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