1 /* Copyright 2002-2026 CS GROUP
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 * CS 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.orbits;
18
19 import org.hipparchus.analysis.interpolation.HermiteInterpolator;
20 import org.hipparchus.util.MathUtils;
21 import org.orekit.errors.OrekitInternalError;
22 import org.orekit.frames.Frame;
23 import org.orekit.time.AbsoluteDate;
24 import org.orekit.time.TimeInterpolator;
25 import org.orekit.utils.CartesianDerivativesFilter;
26 import org.orekit.utils.TimeStampedPVCoordinates;
27 import org.orekit.utils.TimeStampedPVCoordinatesHermiteInterpolator;
28
29 import java.util.List;
30 import java.util.stream.Stream;
31
32 /**
33 * Class using a Hermite interpolator to interpolate orbits. It is assumed that provided samples have the same
34 * gravitational parameter and frame. Otherwise, an error will be thrown
35 * <p>
36 * Depending on given sample orbit type, the interpolation may differ :
37 * <ul>
38 * <li>For Keplerian, Circular and Equinoctial orbits, the interpolated instance is created by polynomial Hermite
39 * interpolation, using derivatives when available. </li>
40 * <li>For Cartesian orbits, the interpolated instance is created using the cartesian derivatives filter given at
41 * instance construction. Hence, it will fall back to Lagrange interpolation if this instance has been designed to not
42 * use derivatives.
43 * </ul>
44 * <p>
45 * In any case, it should be used only with small number of interpolation points (about 10-20 points) in order to avoid
46 * <a href="https://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a> and numerical problems
47 * (including NaN appearing).
48 *
49 * @author Luc Maisonobe
50 * @author Vincent Cucchietti
51 * @see Orbit
52 * @see HermiteInterpolator
53 */
54 public class OrbitHermiteInterpolator extends AbstractOrbitInterpolator {
55
56 /** Filter for derivatives from the sample to use in position-velocity-acceleration interpolation. */
57 private final CartesianDerivativesFilter pvaFilter;
58
59 /**
60 * Constructor with :
61 * <ul>
62 * <li>Default number of interpolation points of {@code DEFAULT_INTERPOLATION_POINTS}</li>
63 * <li>Default extrapolation threshold value ({@code DEFAULT_EXTRAPOLATION_THRESHOLD_SEC} s)</li>
64 * <li>Use of position and two time derivatives during interpolation</li>
65 * </ul>
66 * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
67 * points (about 10-20 points) in order to avoid <a href="https://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
68 * phenomenon</a> and numerical problems (including NaN appearing).
69 *
70 * @param outputInertialFrame output inertial frame
71 */
72 public OrbitHermiteInterpolator(final Frame outputInertialFrame) {
73 this(DEFAULT_INTERPOLATION_POINTS, outputInertialFrame);
74 }
75
76 /**
77 * Constructor with :
78 * <ul>
79 * <li>Default extrapolation threshold value ({@code DEFAULT_EXTRAPOLATION_THRESHOLD_SEC} s)</li>
80 * <li>Use of position and two time derivatives during interpolation</li>
81 * </ul>
82 * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
83 * points (about 10-20 points) in order to avoid <a href="https://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
84 * phenomenon</a> and numerical problems (including NaN appearing).
85 *
86 * @param interpolationPoints number of interpolation points
87 * @param outputInertialFrame output inertial frame
88 */
89 public OrbitHermiteInterpolator(final int interpolationPoints, final Frame outputInertialFrame) {
90 this(interpolationPoints, outputInertialFrame, CartesianDerivativesFilter.USE_PVA);
91 }
92
93 /**
94 * Constructor with default extrapolation threshold value ({@code DEFAULT_EXTRAPOLATION_THRESHOLD_SEC} s).
95 * <p>
96 * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
97 * points (about 10-20 points) in order to avoid <a href="https://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
98 * phenomenon</a> and numerical problems (including NaN appearing).
99 *
100 * @param interpolationPoints number of interpolation points
101 * @param outputInertialFrame output inertial frame
102 * @param pvaFilter filter for derivatives from the sample to use in position-velocity-acceleration interpolation. Used
103 * only when interpolating Cartesian orbits.
104 */
105 public OrbitHermiteInterpolator(final int interpolationPoints, final Frame outputInertialFrame,
106 final CartesianDerivativesFilter pvaFilter) {
107 this(interpolationPoints, DEFAULT_EXTRAPOLATION_THRESHOLD_SEC, outputInertialFrame, pvaFilter);
108 }
109
110 /**
111 * Constructor.
112 * <p>
113 * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
114 * points (about 10-20 points) in order to avoid <a href="https://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
115 * phenomenon</a> and numerical problems (including NaN appearing).
116 *
117 * @param interpolationPoints number of interpolation points
118 * @param extrapolationThreshold extrapolation threshold beyond which the propagation will fail
119 * @param outputInertialFrame output inertial frame
120 * @param pvaFilter filter for derivatives from the sample to use in position-velocity-acceleration interpolation. Used
121 * only when interpolating Cartesian orbits.
122 */
123 public OrbitHermiteInterpolator(final int interpolationPoints, final double extrapolationThreshold,
124 final Frame outputInertialFrame, final CartesianDerivativesFilter pvaFilter) {
125 super(interpolationPoints, extrapolationThreshold, outputInertialFrame);
126 this.pvaFilter = pvaFilter;
127 }
128
129 /** Get filter for derivatives from the sample to use in position-velocity-acceleration interpolation.
130 * @return filter for derivatives from the sample to use in position-velocity-acceleration interpolation
131 */
132 public CartesianDerivativesFilter getPVAFilter() {
133 return pvaFilter;
134 }
135
136 /**
137 * {@inheritDoc}
138 * <p>
139 * Depending on given sample orbit type, the interpolation may differ :
140 * <ul>
141 * <li>For Keplerian, Circular and Equinoctial orbits, the interpolated instance is created by polynomial Hermite
142 * interpolation, using derivatives when available. </li>
143 * <li>For Cartesian orbits, the interpolated instance is created using the cartesian derivatives filter given at
144 * instance construction. Hence, it will fall back to Lagrange interpolation if this instance has been designed to not
145 * use derivatives.
146 * </ul>
147 * If orbit interpolation on large samples is needed, using the {@link
148 * org.orekit.propagation.analytical.Ephemeris} class is a better way than using this
149 * low-level interpolation. The Ephemeris class automatically handles selection of
150 * a neighboring sub-sample with a predefined number of point from a large global sample
151 * in a thread-safe way.
152 *
153 * @param interpolationData interpolation data
154 *
155 * @return interpolated instance at given date
156 */
157 @Override
158 protected Orbit interpolate(final InterpolationData interpolationData) {
159 // Get orbit sample
160 final List<Orbit> sample = interpolationData.getNeighborList();
161
162 // Get information for interpolation
163 final AbsoluteDate interpolationDate = interpolationData.getInterpolationDate();
164 final Orbit firstEntry = sample.getFirst();
165 final OrbitType orbitType = firstEntry.getType();
166 final Frame inputFrame = firstEntry.getFrame();
167 final Frame outputFrame = getOutputInertialFrame();
168
169 final Orbit interpolated;
170 if (orbitType == OrbitType.CARTESIAN) {
171 interpolated = interpolateCartesian(interpolationDate, inputFrame, sample);
172 }
173 else {
174 interpolated = interpolateCommon(interpolationDate, inputFrame, sample, orbitType);
175 }
176
177 // Return interpolated if input and output frame are the same
178 if (inputFrame.equals(outputFrame)) {
179 return interpolated;
180 }
181
182 // Otherwise, express in the output frame
183 return interpolated.inFrame(outputFrame);
184 }
185
186 /**
187 * Interpolate Cartesian orbit using a specific method for Cartesian orbit. Assumes that the input frame is the same
188 * throughout the sample.
189 *
190 * @param interpolationDate interpolation date
191 * @param inputFrame input frame
192 * @param sample orbits sample
193 * @return interpolated Cartesian orbit
194 */
195 private CartesianOrbit interpolateCartesian(final AbsoluteDate interpolationDate,
196 final Frame inputFrame,
197 final List<Orbit> sample) {
198
199 // Create time stamped position-velocity-acceleration Hermite interpolator
200 final TimeInterpolator<TimeStampedPVCoordinates> interpolator =
201 new TimeStampedPVCoordinatesHermiteInterpolator(getNbInterpolationPoints(),
202 getExtrapolationThreshold(),
203 pvaFilter);
204
205 // Convert sample to stream
206 final Stream<Orbit> sampleStream = sample.stream();
207
208 // Map time stamped position-velocity-acceleration coordinates
209 final Stream<TimeStampedPVCoordinates> sampleTimeStampedPV = sampleStream.map(Orbit::getPVCoordinates);
210
211 // Interpolate PVA
212 final TimeStampedPVCoordinates interpolated = interpolator.interpolate(interpolationDate, sampleTimeStampedPV);
213
214 // Use first entry gravitational parameter (initially checked that it is consistent throughout the sample)
215 final double mu = sample.getFirst().getMu();
216
217 return new CartesianOrbit(interpolated,
218 inputFrame,
219 interpolationDate,
220 mu);
221 }
222
223 /**
224 * Method gathering common parts of interpolation between circular, equinoctial and keplerian orbit.
225 *
226 * @param interpolationDate interpolation date
227 * @param inputFrame input frame
228 * @param orbits orbits sample
229 * @param orbitType interpolation method to use
230 *
231 * @return interpolated orbit
232 */
233 private Orbit interpolateCommon(final AbsoluteDate interpolationDate,
234 final Frame inputFrame,
235 final List<Orbit> orbits,
236 final OrbitType orbitType) {
237
238 // First pass to check if derivatives are available throughout the sample
239 boolean useDerivatives = true;
240 for (final Orbit orbit : orbits) {
241 useDerivatives = useDerivatives && orbit.hasNonKeplerianAcceleration();
242 }
243
244 // Use first entry gravitational parameter
245 final double mu = orbits.getFirst().getMu();
246
247 // Interpolate and build a new instance
248 final double[][] interpolated;
249 switch (orbitType) {
250 case CIRCULAR:
251 interpolated = interpolateCircular(interpolationDate, orbits, useDerivatives);
252 return new CircularOrbit(interpolated[0][0], interpolated[0][1], interpolated[0][2],
253 interpolated[0][3], interpolated[0][4], interpolated[0][5],
254 interpolated[1][0], interpolated[1][1], interpolated[1][2],
255 interpolated[1][3], interpolated[1][4], interpolated[1][5],
256 PositionAngleType.MEAN, inputFrame, interpolationDate, mu);
257 case KEPLERIAN:
258 interpolated = interpolateKeplerian(interpolationDate, orbits, useDerivatives);
259 return new KeplerianOrbit(interpolated[0][0], interpolated[0][1], interpolated[0][2],
260 interpolated[0][3], interpolated[0][4], interpolated[0][5],
261 interpolated[1][0], interpolated[1][1], interpolated[1][2],
262 interpolated[1][3], interpolated[1][4], interpolated[1][5],
263 PositionAngleType.MEAN, inputFrame, interpolationDate, mu);
264 case EQUINOCTIAL:
265 interpolated = interpolateEquinoctial(interpolationDate, orbits, useDerivatives);
266 return new EquinoctialOrbit(interpolated[0][0], interpolated[0][1], interpolated[0][2],
267 interpolated[0][3], interpolated[0][4], interpolated[0][5],
268 interpolated[1][0], interpolated[1][1], interpolated[1][2],
269 interpolated[1][3], interpolated[1][4], interpolated[1][5],
270 PositionAngleType.MEAN, inputFrame, interpolationDate, mu);
271 default:
272 // Should never happen
273 throw new OrekitInternalError(null);
274 }
275
276 }
277
278 /**
279 * Build interpolating functions for circular orbit parameters.
280 *
281 * @param interpolationDate interpolation date
282 * @param orbits orbits sample
283 * @param useDerivatives flag defining if derivatives are available throughout the sample
284 *
285 * @return interpolating functions for circular orbit parameters
286 */
287 private double[][] interpolateCircular(final AbsoluteDate interpolationDate, final List<Orbit> orbits,
288 final boolean useDerivatives) {
289
290 // Set up an interpolator
291 final HermiteInterpolator interpolator = new HermiteInterpolator();
292
293 // Second pass to feed interpolator
294 AbsoluteDate previousDate = null;
295 double previousRAAN = Double.NaN;
296 double previousAlphaM = Double.NaN;
297 for (final Orbit orbit : orbits) {
298 final CircularOrbit circ = (CircularOrbit) OrbitType.CIRCULAR.convertType(orbit);
299 final double continuousRAAN;
300 final double continuousAlphaM;
301 if (previousDate == null) {
302 continuousRAAN = circ.getRightAscensionOfAscendingNode();
303 continuousAlphaM = circ.getAlphaM();
304 }
305 else {
306 final double dt = circ.getDate().durationFrom(previousDate);
307 final double keplerAM = previousAlphaM + circ.getKeplerianMeanMotion() * dt;
308 continuousRAAN = MathUtils.normalizeAngle(circ.getRightAscensionOfAscendingNode(), previousRAAN);
309 continuousAlphaM = MathUtils.normalizeAngle(circ.getAlphaM(), keplerAM);
310 }
311 previousDate = circ.getDate();
312 previousRAAN = continuousRAAN;
313 previousAlphaM = continuousAlphaM;
314 if (useDerivatives) {
315 interpolator.addSamplePoint(circ.getDate().durationFrom(interpolationDate),
316 new double[] { circ.getA(),
317 circ.getCircularEx(),
318 circ.getCircularEy(),
319 circ.getI(),
320 continuousRAAN,
321 continuousAlphaM },
322 new double[] { circ.getADot(),
323 circ.getCircularExDot(),
324 circ.getCircularEyDot(),
325 circ.getIDot(),
326 circ.getRightAscensionOfAscendingNodeDot(),
327 circ.getAlphaMDot() });
328 }
329 else {
330 interpolator.addSamplePoint(circ.getDate().durationFrom(interpolationDate),
331 new double[] { circ.getA(),
332 circ.getCircularEx(),
333 circ.getCircularEy(),
334 circ.getI(),
335 continuousRAAN,
336 continuousAlphaM });
337 }
338 }
339
340 return interpolator.derivatives(0.0, 1);
341 }
342
343 /**
344 * Build interpolating functions for keplerian orbit parameters.
345 *
346 * @param interpolationDate interpolation date
347 * @param orbits orbits sample
348 * @param useDerivatives flag defining if derivatives are available throughout the sample
349 *
350 * @return interpolating functions for keplerian orbit parameters
351 */
352 private double[][] interpolateKeplerian(final AbsoluteDate interpolationDate, final List<Orbit> orbits,
353 final boolean useDerivatives) {
354
355 // Set up an interpolator
356 final HermiteInterpolator interpolator = new HermiteInterpolator();
357
358 // Second pass to feed interpolator
359 AbsoluteDate previousDate = null;
360 double previousPA = Double.NaN;
361 double previousRAAN = Double.NaN;
362 double previousM = Double.NaN;
363 for (final Orbit orbit : orbits) {
364 final KeplerianOrbit kep = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(orbit);
365 final double continuousPA;
366 final double continuousRAAN;
367 final double continuousM;
368 if (previousDate == null) {
369 continuousPA = kep.getPerigeeArgument();
370 continuousRAAN = kep.getRightAscensionOfAscendingNode();
371 continuousM = kep.getMeanAnomaly();
372 }
373 else {
374 final double dt = kep.getDate().durationFrom(previousDate);
375 final double keplerM = previousM + kep.getKeplerianMeanMotion() * dt;
376 continuousPA = MathUtils.normalizeAngle(kep.getPerigeeArgument(), previousPA);
377 continuousRAAN = MathUtils.normalizeAngle(kep.getRightAscensionOfAscendingNode(), previousRAAN);
378 continuousM = MathUtils.normalizeAngle(kep.getMeanAnomaly(), keplerM);
379 }
380 previousDate = kep.getDate();
381 previousPA = continuousPA;
382 previousRAAN = continuousRAAN;
383 previousM = continuousM;
384 if (useDerivatives) {
385 interpolator.addSamplePoint(kep.getDate().durationFrom(interpolationDate),
386 new double[] { kep.getA(),
387 kep.getE(),
388 kep.getI(),
389 continuousPA,
390 continuousRAAN,
391 continuousM },
392 new double[] { kep.getADot(),
393 kep.getEDot(),
394 kep.getIDot(),
395 kep.getPerigeeArgumentDot(),
396 kep.getRightAscensionOfAscendingNodeDot(),
397 kep.getMeanAnomalyDot() });
398 }
399 else {
400 interpolator.addSamplePoint(kep.getDate().durationFrom(interpolationDate),
401 new double[] { kep.getA(),
402 kep.getE(),
403 kep.getI(),
404 continuousPA,
405 continuousRAAN,
406 continuousM });
407 }
408 }
409
410 return interpolator.derivatives(0.0, 1);
411 }
412
413 /**
414 * Build interpolating functions for equinoctial orbit parameters.
415 *
416 * @param interpolationDate interpolation date
417 * @param orbits orbits sample
418 * @param useDerivatives flag defining if derivatives are available throughout the sample
419 *
420 * @return interpolating functions for equinoctial orbit parameters
421 */
422 private double[][] interpolateEquinoctial(final AbsoluteDate interpolationDate, final List<Orbit> orbits,
423 final boolean useDerivatives) {
424
425 // Set up an interpolator
426 final HermiteInterpolator interpolator = new HermiteInterpolator();
427
428 // Second pass to feed interpolator
429 AbsoluteDate previousDate = null;
430 double previousLm = Double.NaN;
431 for (final Orbit orbit : orbits) {
432 final EquinoctialOrbit equi = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(orbit);
433 final double continuousLm;
434 if (previousDate == null) {
435 continuousLm = equi.getLM();
436 }
437 else {
438 final double dt = equi.getDate().durationFrom(previousDate);
439 final double keplerLm = previousLm + equi.getKeplerianMeanMotion() * dt;
440 continuousLm = MathUtils.normalizeAngle(equi.getLM(), keplerLm);
441 }
442 previousDate = equi.getDate();
443 previousLm = continuousLm;
444 if (useDerivatives) {
445 interpolator.addSamplePoint(equi.getDate().durationFrom(interpolationDate),
446 new double[] { equi.getA(),
447 equi.getEquinoctialEx(),
448 equi.getEquinoctialEy(),
449 equi.getHx(),
450 equi.getHy(),
451 continuousLm },
452 new double[] {
453 equi.getADot(),
454 equi.getEquinoctialExDot(),
455 equi.getEquinoctialEyDot(),
456 equi.getHxDot(),
457 equi.getHyDot(),
458 equi.getLMDot() });
459 }
460 else {
461 interpolator.addSamplePoint(equi.getDate().durationFrom(interpolationDate),
462 new double[] { equi.getA(),
463 equi.getEquinoctialEx(),
464 equi.getEquinoctialEy(),
465 equi.getHx(),
466 equi.getHy(),
467 continuousLm });
468 }
469 }
470
471 return interpolator.derivatives(0.0, 1);
472 }
473 }