1 /* Copyright 2022-2025 Luc Maisonobe
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.files.sp3;
18
19 import java.util.stream.Stream;
20
21 import org.hipparchus.analysis.interpolation.HermiteInterpolator;
22 import org.hipparchus.geometry.euclidean.threed.Vector3D;
23 import org.orekit.time.AbsoluteDate;
24 import org.orekit.time.AbstractTimeInterpolator;
25
26 /** Interpolator for {@link SP3Coordinate SP3 coordinates}.
27 * <p>
28 * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation points
29 * (about 10-20 points) in order to avoid <a href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
30 * and numerical problems (including NaN appearing).
31 * </p>
32 * <p>
33 * If some clock or clock rate are present in the SP3 files as default values (999999.999999), then they
34 * are replaced by {@code Double.NaN} during parsing, so the interpolation will exhibit NaNs, but the
35 * positions will be properly interpolated.
36 * </p>
37 *
38 * @author Luc Maisonobe
39 * @see HermiteInterpolator
40 * @see SP3Coordinate
41 * @since 12.0
42 */
43 public class SP3CoordinateHermiteInterpolator extends AbstractTimeInterpolator<SP3Coordinate> {
44
45 /** Flag for using velocity and clock rate. */
46 private final boolean useRates;
47
48 /**
49 * Constructor.
50 * <p>
51 * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
52 * points (about 10-20 points) in order to avoid <a href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
53 * phenomenon</a> and numerical problems (including NaN appearing).
54 * </p>
55 *
56 * @param interpolationPoints number of interpolation points
57 * @param extrapolationThreshold extrapolation threshold beyond which the propagation will fail
58 * @param useRates if true, use velocity and clock rates for interpolation
59 */
60 public SP3CoordinateHermiteInterpolator(final int interpolationPoints,
61 final double extrapolationThreshold,
62 final boolean useRates) {
63 super(interpolationPoints, extrapolationThreshold);
64 this.useRates = useRates;
65 }
66
67 /**
68 * {@inheritDoc}
69 * <p>
70 * The interpolated instance is created by polynomial Hermite interpolation ensuring velocity remains the exact
71 * derivative of position.
72 * <p>
73 * Note that even if first time derivatives (velocities) from sample can be ignored, the interpolated instance always
74 * includes interpolated derivatives. This feature can be used explicitly to compute these derivatives when it would be
75 * too complex to compute them from an analytical formula: just compute a few sample points from the explicit formula and
76 * set the derivatives to zero in these sample points, then use interpolation to add derivatives consistent with the
77 * positions.
78 */
79 @Override
80 protected SP3Coordinate interpolate(final InterpolationData interpolationData) {
81
82 // Get date
83 final AbsoluteDate date = interpolationData.getInterpolationDate();
84
85 // Convert sample to stream
86 final Stream<SP3Coordinate> sample = interpolationData.getNeighborList().stream();
87
88 // Set up an interpolator taking derivatives into account
89 final HermiteInterpolator interpolator = new HermiteInterpolator();
90
91 // Add sample points
92 if (useRates) {
93 // populate sample with position, clock, velocity and clock rate data
94 sample.forEach(c -> interpolator.addSamplePoint(c.getDate().durationFrom(date),
95 new double[] {
96 c.getPosition().getX(),
97 c.getPosition().getY(),
98 c.getPosition().getZ(),
99 c.getClockCorrection(),
100 },
101 new double[] {
102 c.getVelocity().getX(),
103 c.getVelocity().getY(),
104 c.getVelocity().getZ(),
105 c.getClockRateChange(),
106 }));
107 } else {
108 // populate sample with position and clock data, ignoring velocity and clock rate
109 sample.forEach(c -> interpolator.addSamplePoint(c.getDate().durationFrom(date),
110 new double[] {
111 c.getPosition().getX(),
112 c.getPosition().getY(),
113 c.getPosition().getZ(),
114 c.getClockCorrection(),
115 }));
116 }
117
118 // Interpolate
119 final double[][] interpolated = interpolator.derivatives(0.0, 1);
120
121 // Build a new interpolated instance
122 return new SP3Coordinate(date,
123 new Vector3D(interpolated[0][0], interpolated[0][1], interpolated[0][2]), null,
124 new Vector3D(interpolated[1][0], interpolated[1][1], interpolated[1][2]), null,
125 interpolated[0][3], Double.NaN,
126 interpolated[1][3], Double.NaN,
127 false, false, false, false);
128
129 }
130
131 }