1 /* Copyright 2002-2024 Joseph Reed
2 * Licensed to CS GROUP (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * Joseph Reed licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.utils;
18
19 import java.util.Map.Entry;
20 import java.util.TreeMap;
21
22 import org.hipparchus.geometry.euclidean.threed.Vector3D;
23 import org.hipparchus.geometry.spherical.twod.Circle;
24 import org.hipparchus.geometry.spherical.twod.S2Point;
25 import org.hipparchus.util.FastMath;
26 import org.orekit.bodies.GeodeticPoint;
27 import org.orekit.bodies.LoxodromeArc;
28 import org.orekit.bodies.OneAxisEllipsoid;
29 import org.orekit.frames.Frame;
30 import org.orekit.frames.TopocentricFrame;
31 import org.orekit.time.AbsoluteDate;
32
33 /** Builder class, enabling incremental building of an {@link PVCoordinatesProvider}
34 * instance using waypoints defined on an ellipsoid.
35 *
36 * Given a series of waypoints ({@code (date, point)} tuples),
37 * build a {@link PVCoordinatesProvider} representing the path.
38 * The static methods provide implementations for the most common path definitions
39 * (cartesian, great-circle, loxodrome). If these methods are insufficient,
40 * the public constructor provides a way to customize the path definition.
41 *
42 * This class connects the path segments using the {@link AggregatedPVCoordinatesProvider}.
43 * As such, no effort is made to smooth the velocity between segments.
44 * While position is unaffected, the velocity may be discontinuous between adjacent time points.
45 * Thus, care should be taken when modeling paths with abrupt direction changes
46 * (e.g. fast-moving aircraft); understand how the {@link PVCoordinatesProvider}
47 * will be used in the particular application.
48 *
49 * @author Joe Reed
50 * @since 11.3
51 */
52 public class WaypointPVBuilder {
53
54 /** Factory used to create intermediate pv providers between waypoints. */
55 private final InterpolationFactory factory;
56
57 /** Central body, on which the waypoints are defined. */
58 private final OneAxisEllipsoid body;
59
60 /** Set of waypoints, indexed by time. */
61 private final TreeMap<AbsoluteDate, GeodeticPoint> waypoints;
62
63 /** Whether the resulting provider should be invalid or constant prior to the first waypoint. */
64 private boolean invalidBefore;
65
66 /** Whether the resulting provider should be invalid or constant after to the last waypoint. */
67 private boolean invalidAfter;
68
69 /** Create a new instance.
70 * @param factory The factory used to create the intermediate coordinate providers between waypoints.
71 * @param body The central body, on which the way points are defined.
72 */
73 public WaypointPVBuilder(final InterpolationFactory factory, final OneAxisEllipsoid body) {
74 this.factory = factory;
75 this.body = body;
76 this.waypoints = new TreeMap<>();
77 this.invalidBefore = true;
78 this.invalidAfter = true;
79 }
80
81 /** Construct a waypoint builder interpolating points using a linear cartesian interpolation.
82 *
83 * @param body the reference ellipsoid on which the waypoints are defined.
84 * @return the waypoint builder
85 */
86 public static WaypointPVBuilder cartesianBuilder(final OneAxisEllipsoid body) {
87 return new WaypointPVBuilder(CartesianWaypointPVProv::new, body);
88 }
89
90 /** Construct a waypoint builder interpolating points using a loxodrome (or Rhumbline).
91 *
92 * @param body the reference ellipsoid on which the waypoints are defined.
93 * @return the waypoint builder
94 */
95 public static WaypointPVBuilder loxodromeBuilder(final OneAxisEllipsoid body) {
96 return new WaypointPVBuilder(LoxodromeWaypointPVProv::new, body);
97 }
98
99 /** Construct a waypoint builder interpolating points using a great-circle.
100 *
101 * The altitude of the intermediate points is linearly interpolated from the bounding waypoints.
102 * Extrapolating before the first waypoint or after the last waypoint may result in undefined altitudes.
103 *
104 * @param body the reference ellipsoid on which the waypoints are defined.
105 * @return the waypoint builder
106 */
107 public static WaypointPVBuilder greatCircleBuilder(final OneAxisEllipsoid body) {
108 return new WaypointPVBuilder(GreatCircleWaypointPVProv::new, body);
109 }
110
111 /** Add a waypoint.
112 *
113 * @param point the waypoint location
114 * @param date the waypoint time
115 * @return this instance
116 */
117 public WaypointPVBuilder addWaypoint(final GeodeticPoint point, final AbsoluteDate date) {
118 waypoints.put(date, point);
119 return this;
120 }
121
122 /** Indicate the resulting {@link PVCoordinatesProvider} should be invalid before the first waypoint.
123 *
124 * @return this instance
125 */
126 public WaypointPVBuilder invalidBefore() {
127 invalidBefore = true;
128 return this;
129 }
130
131 /** Indicate the resulting {@link PVCoordinatesProvider} provide
132 * a constant location of the first waypoint prior to the first time.
133 *
134 * @return this instance
135 */
136 public WaypointPVBuilder constantBefore() {
137 invalidBefore = false;
138 return this;
139 }
140
141 /** Indicate the resulting {@link PVCoordinatesProvider} should be invalid after the last waypoint.
142 *
143 * @return this instance
144 */
145 public WaypointPVBuilder invalidAfter() {
146 invalidAfter = true;
147 return this;
148 }
149
150 /** Indicate the resulting {@link PVCoordinatesProvider} provide
151 * a constant location of the last waypoint after to the last time.
152 *
153 * @return this instance
154 */
155 public WaypointPVBuilder constantAfter() {
156 invalidAfter = false;
157 return this;
158 }
159
160 /** Build a {@link PVCoordinatesProvider} from the waypoints added to this builder.
161 *
162 * @return the coordinates provider instance.
163 */
164 public PVCoordinatesProvider build() {
165 final PVCoordinatesProvider initialProvider = createInitial(waypoints.firstKey(),
166 waypoints.firstEntry().getValue());
167 final AggregatedPVCoordinatesProvider.Builder builder = new AggregatedPVCoordinatesProvider.Builder(initialProvider);
168
169 Entry<AbsoluteDate, GeodeticPoint> previousEntry = null;
170 for (final Entry<AbsoluteDate, GeodeticPoint> entry: waypoints.entrySet()) {
171 if (previousEntry != null) {
172 builder.addPVProviderAfter(previousEntry.getKey(),
173 factory.create(previousEntry.getKey(),
174 previousEntry.getValue(),
175 entry.getKey(),
176 entry.getValue(),
177 body),
178 true);
179 }
180 previousEntry = entry;
181 }
182 // add the point so we're valid at the final waypoint
183 builder.addPVProviderAfter(previousEntry.getKey(),
184 new ConstantPVCoordinatesProvider(previousEntry.getValue(), body),
185 true);
186 // add the final propagator after the final waypoint
187 builder.addPVProviderAfter(previousEntry.getKey().shiftedBy(Double.MIN_VALUE),
188 createFinal(previousEntry.getKey(), previousEntry.getValue()),
189 true);
190
191 return builder.build();
192 }
193
194 /** Create the initial provider.
195 *
196 * This method uses the internal {@code validBefore} flag to
197 * either return an invalid PVCoordinatesProvider or a constant one.
198 *
199 * @param firstDate the date at which the first waypoint is reached
200 * and this provider will no longer be called
201 * @param firstPoint the first waypoint
202 * @return the coordinate provider
203 */
204 protected PVCoordinatesProvider createInitial(final AbsoluteDate firstDate, final GeodeticPoint firstPoint) {
205 if (invalidBefore) {
206 return new AggregatedPVCoordinatesProvider.InvalidPVProvider();
207 }
208 else {
209 return new ConstantPVCoordinatesProvider(firstPoint, body);
210 }
211 }
212
213 /** Create the final provider.
214 *
215 * This method uses the internal {@code validAfter} flag to
216 * either return an invalid PVCoordinatesProvider or a constant one.
217 *
218 * @param lastDate the date at which the last waypoint is reached
219 * and this provider will be called
220 * @param lastPoint the last waypoint
221 * @return the coordinate provider
222 */
223 protected PVCoordinatesProvider createFinal(final AbsoluteDate lastDate, final GeodeticPoint lastPoint) {
224 if (invalidAfter) {
225 return new AggregatedPVCoordinatesProvider.InvalidPVProvider();
226 }
227 else {
228 return new ConstantPVCoordinatesProvider(lastPoint, body);
229 }
230 }
231
232 /**
233 * Factory interface, creating the {@link PVCoordinatesProvider} instances between the provided waypoints.
234 */
235 @FunctionalInterface
236 public interface InterpolationFactory {
237
238 /** Create a {@link PVCoordinatesProvider} which interpolates between the provided waypoints.
239 *
240 * @param date1 the first waypoint's date
241 * @param point1 the first waypoint's location
242 * @param date2 the second waypoint's date
243 * @param point2 the second waypoint's location
244 * @param body the body on which the waypoints are defined
245 * @return a {@link PVCoordinatesProvider} providing the locations at times between the waypoints.
246 */
247 PVCoordinatesProvider create(AbsoluteDate date1, GeodeticPoint point1,
248 AbsoluteDate date2, GeodeticPoint point2,
249 OneAxisEllipsoid body);
250 }
251
252 /**
253 * Coordinate provider interpolating along the great-circle between two points.
254 */
255 static class GreatCircleWaypointPVProv implements PVCoordinatesProvider {
256
257 /** Great circle estimation. */
258 private final Circle circle;
259 /** Duration between the two points (seconds). */
260 private final double duration;
261 /** Phase along the circle of the first point. */
262 private final double phase0;
263 /** Phase length from the first point to the second. */
264 private final double phaseLength;
265 /** Time at which interpolation results in the initial point. */
266 private final AbsoluteDate t0;
267 /** Body on which the great circle is defined. */
268 private final OneAxisEllipsoid body;
269 /** Phase of one second. */
270 private double oneSecondPhase;
271 /** Altitude of the initial point. */
272 private double initialAltitude;
273 /** Time-derivative of the altitude. */
274 private double altitudeSlope;
275
276 /** Class constructor. Aligns to the {@link InterpolationFactory} functional interface.
277 *
278 * @param date1 the first waypoint's date
279 * @param point1 the first waypoint's location
280 * @param date2 the second waypoint's date
281 * @param point2 the second waypoint's location
282 * @param body the body on which the waypoints are defined
283 * @see InterpolationFactory
284 */
285 GreatCircleWaypointPVProv(final AbsoluteDate date1, final GeodeticPoint point1,
286 final AbsoluteDate date2, final GeodeticPoint point2,
287 final OneAxisEllipsoid body) {
288 this.t0 = date1;
289 this.duration = date2.durationFrom(date1);
290 this.body = body;
291 final S2Point s0 = toSpherical(point1);
292 final S2Point s1 = toSpherical(point2);
293 circle = new Circle(s0, s1, 1e-9);
294
295 phase0 = circle.getPhase(s0.getVector());
296 phaseLength = circle.getPhase(s1.getVector()) - phase0;
297
298 oneSecondPhase = phaseLength / duration;
299 altitudeSlope = (point2.getAltitude() - point1.getAltitude()) / duration;
300 initialAltitude = point1.getAltitude();
301 }
302
303 @Override
304 public Vector3D getPosition(final AbsoluteDate date, final Frame frame) {
305 final double d = date.durationFrom(t0);
306 final double fraction = d / duration;
307 final double phase = fraction * phaseLength;
308
309 final S2Point sp = new S2Point(circle.getPointAt(phase0 + phase));
310 final GeodeticPoint point = toGeodetic(sp, initialAltitude + d * altitudeSlope);
311 final Vector3D p = body.transform(point);
312
313 return body.getBodyFrame().getStaticTransformTo(frame, date).transformPosition(p);
314
315 }
316
317 @Override
318 public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame) {
319 final double d = date.durationFrom(t0);
320 final double fraction = d / duration;
321 final double phase = fraction * phaseLength;
322
323 final S2Point sp = new S2Point(circle.getPointAt(phase0 + phase));
324 final GeodeticPoint point = toGeodetic(sp, initialAltitude + d * altitudeSlope);
325 final Vector3D p = body.transform(point);
326
327 // add 1 second to get another point along the circle, to use for velocity
328 final S2Point sp2 = new S2Point(circle.getPointAt(phase0 + phase + oneSecondPhase));
329 final GeodeticPoint point2 = toGeodetic(sp2, initialAltitude + (d + 1) * altitudeSlope);
330 final Vector3D p2 = body.transform(point2);
331 final Vector3D v = p2.subtract(p);
332
333 final TimeStampedPVCoordinates tpv = new TimeStampedPVCoordinates(date, p, v);
334 return body.getBodyFrame().getTransformTo(frame, date).transformPVCoordinates(tpv);
335 }
336
337 /** Converts the given geodetic point to a point on the 2-sphere.
338 * @param point input geodetic point
339 * @return a point on the 2-sphere
340 */
341 static S2Point toSpherical(final GeodeticPoint point) {
342 return new S2Point(point.getLongitude(), 0.5 * FastMath.PI - point.getLatitude());
343 }
344
345 /** Converts a 2-sphere point to a geodetic point.
346 * @param point point on the 2-sphere
347 * @param alt point altitude
348 * @return a geodetic point
349 */
350 static GeodeticPoint toGeodetic(final S2Point point, final double alt) {
351 return new GeodeticPoint(0.5 * FastMath.PI - point.getPhi(), point.getTheta(), alt);
352 }
353 }
354
355 /**
356 * Coordinate provider interpolating along the loxodrome between two points.
357 */
358 static class LoxodromeWaypointPVProv implements PVCoordinatesProvider {
359
360 /** Arc along which the interpolation occurs. */
361 private final LoxodromeArc arc;
362 /** Time at which the interpolation begins (at arc start). */
363 private final AbsoluteDate t0;
364 /** Total duration to get the length of the arc (seconds). */
365 private final double duration;
366 /** Velocity along the arc (m/s). */
367 private final double velocity;
368
369 /** Class constructor. Aligns to the {@link InterpolationFactory} functional interface.
370 *
371 * @param date1 the first waypoint's date
372 * @param point1 the first waypoint's location
373 * @param date2 the second waypoint's date
374 * @param point2 the second waypoint's location
375 * @param body the body on which the waypoints are defined
376 * @see InterpolationFactory
377 */
378 LoxodromeWaypointPVProv(final AbsoluteDate date1, final GeodeticPoint point1, final AbsoluteDate date2,
379 final GeodeticPoint point2, final OneAxisEllipsoid body) {
380 this.arc = new LoxodromeArc(point1, point2, body);
381 this.t0 = date1;
382 this.duration = date2.durationFrom(date1);
383 this.velocity = arc.getDistance() / duration;
384 }
385
386 @Override
387 public Vector3D getPosition(final AbsoluteDate date, final Frame frame) {
388 final double fraction = date.durationFrom(t0) / duration;
389 final GeodeticPoint point = arc.calculatePointAlongArc(fraction);
390 final Vector3D p = arc.getBody().transform(point);
391
392 return arc.getBody().getBodyFrame().getStaticTransformTo(frame, date).transformPosition(p);
393 }
394
395 @Override
396 public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame) {
397 final double fraction = date.durationFrom(t0) / duration;
398 final GeodeticPoint point = arc.calculatePointAlongArc(fraction);
399 final Vector3D p = arc.getBody().transform(point);
400 final Vector3D vp = arc.getBody().transform(
401 new TopocentricFrame(arc.getBody(), point, "frame")
402 .pointAtDistance(arc.getAzimuth(), 0, velocity));
403
404 final TimeStampedPVCoordinates tpv = new TimeStampedPVCoordinates(date, p, vp.subtract(p));
405 return arc.getBody().getBodyFrame().getTransformTo(frame, date).transformPVCoordinates(tpv);
406 }
407 }
408
409 /**
410 * Coordinate provider interpolating along the cartesian (3-space) line between two points.
411 */
412 static class CartesianWaypointPVProv implements PVCoordinatesProvider {
413
414 /** Date at which the position is valid. */
415 private final AbsoluteDate t0;
416 /** Initial point. */
417 private final Vector3D p0;
418 /** Velocity. */
419 private final Vector3D vel;
420 /** Frame in which the point and velocity are defined. */
421 private final Frame sourceFrame;
422
423 /** Class constructor. Aligns to the {@link InterpolationFactory} functional interface.
424 *
425 * @param date1 the first waypoint's date
426 * @param point1 the first waypoint's location
427 * @param date2 the second waypoint's date
428 * @param point2 the second waypoint's location
429 * @param body the body on which the waypoints are defined
430 * @see InterpolationFactory
431 */
432 CartesianWaypointPVProv(final AbsoluteDate date1, final GeodeticPoint point1,
433 final AbsoluteDate date2, final GeodeticPoint point2,
434 final OneAxisEllipsoid body) {
435 this.t0 = date1;
436 this.p0 = body.transform(point1);
437 this.vel = body.transform(point2).subtract(p0).scalarMultiply(1. / date2.durationFrom(t0));
438 this.sourceFrame = body.getBodyFrame();
439 }
440
441 @Override
442 public Vector3D getPosition(final AbsoluteDate date, final Frame frame) {
443 final double d = date.durationFrom(t0);
444 final Vector3D p = p0.add(vel.scalarMultiply(d));
445 return sourceFrame.getStaticTransformTo(frame, date).transformPosition(p);
446 }
447
448 @Override
449 public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame) {
450 final double d = date.durationFrom(t0);
451 final Vector3D p = p0.add(vel.scalarMultiply(d));
452 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, p, vel);
453 return sourceFrame.getTransformTo(frame, date).transformPVCoordinates(pv);
454 }
455
456 }
457 }