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 }