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.utils;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.analysis.interpolation.FieldHermiteInterpolator;
21 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
22 import org.orekit.errors.OrekitInternalError;
23 import org.orekit.frames.Frame;
24 import org.orekit.time.AbstractFieldTimeInterpolator;
25 import org.orekit.time.FieldAbsoluteDate;
26
27 import java.util.List;
28
29 /**
30 * Class using a Hermite interpolator to interpolate absolute position-velocity-acceleration coordinates.
31 * <p>
32 * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation points
33 * (about 10-20 points) in order to avoid <a href="https://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
34 * and numerical problems (including NaN appearing).
35 *
36 * @author Luc Maisonobe
37 * @author Vincent Cucchietti
38 * @see FieldHermiteInterpolator
39 * @see FieldAbsolutePVCoordinates
40 * @param <KK> type of the field elements
41 */
42 public class FieldAbsolutePVCoordinatesHermiteInterpolator<KK extends CalculusFieldElement<KK>>
43 extends AbstractFieldTimeInterpolator<FieldAbsolutePVCoordinates<KK>, KK> {
44
45 /** Filter for derivatives from the sample to use in interpolation. */
46 private final CartesianDerivativesFilter filter;
47
48 /** Output frame for the interpolated instance. */
49 private final Frame outputFrame;
50
51 /**
52 * Constructor with :
53 * <ul>
54 * <li>Default number of interpolation points of {@code DEFAULT_INTERPOLATION_POINTS}</li>
55 * <li>Default extrapolation threshold value ({@code DEFAULT_EXTRAPOLATION_THRESHOLD_SEC} s)</li>
56 * <li>Use of position and two time derivatives during interpolation</li>
57 * </ul>
58 * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
59 * points (about 10-20 points) in order to avoid <a href="https://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
60 * phenomenon</a> and numerical problems (including NaN appearing).
61 *
62 * @param outputFrame frame for the interpolated instance
63 */
64 public FieldAbsolutePVCoordinatesHermiteInterpolator(final Frame outputFrame) {
65 this(DEFAULT_INTERPOLATION_POINTS, outputFrame);
66 }
67
68 /**
69 * Constructor with :
70 * <ul>
71 * <li>Default extrapolation threshold value ({@code DEFAULT_EXTRAPOLATION_THRESHOLD_SEC} s)</li>
72 * <li>Use of position and two time derivatives during interpolation</li>
73 * </ul>
74 * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
75 * points (about 10-20 points) in order to avoid <a href="https://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
76 * phenomenon</a> and numerical problems (including NaN appearing).
77 *
78 * @param interpolationPoints number of interpolation points
79 * @param outputFrame frame for the interpolated instance
80 */
81 public FieldAbsolutePVCoordinatesHermiteInterpolator(final int interpolationPoints, final Frame outputFrame) {
82 this(interpolationPoints, outputFrame, CartesianDerivativesFilter.USE_PVA);
83 }
84
85 /**
86 * Constructor with default extrapolation threshold value ({@code DEFAULT_EXTRAPOLATION_THRESHOLD_SEC} s).
87 * <p>
88 * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
89 * points (about 10-20 points) in order to avoid <a href="https://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
90 * phenomenon</a> and numerical problems (including NaN appearing).
91 *
92 * @param interpolationPoints number of interpolation points
93 * @param outputFrame frame for the interpolated instance
94 * @param filter filter for derivatives from the sample to use in interpolation
95 */
96 public FieldAbsolutePVCoordinatesHermiteInterpolator(final int interpolationPoints, final Frame outputFrame,
97 final CartesianDerivativesFilter filter) {
98 super(interpolationPoints, DEFAULT_EXTRAPOLATION_THRESHOLD_SEC);
99 this.outputFrame = outputFrame;
100 this.filter = filter;
101 }
102
103 /**
104 * Constructor.
105 * <p>
106 * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
107 * points (about 10-20 points) in order to avoid <a href="https://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
108 * phenomenon</a> and numerical problems (including NaN appearing).
109 *
110 * @param interpolationPoints number of interpolation points
111 * @param extrapolationThreshold extrapolation threshold beyond which the propagation will fail
112 * @param outputFrame frame for the interpolated instance
113 * @param filter filter for derivatives from the sample to use in interpolation
114 */
115 public FieldAbsolutePVCoordinatesHermiteInterpolator(final int interpolationPoints, final double extrapolationThreshold,
116 final Frame outputFrame, final CartesianDerivativesFilter filter) {
117 super(interpolationPoints, extrapolationThreshold);
118 this.outputFrame = outputFrame;
119 this.filter = filter;
120 }
121
122 /** Get filter for derivatives from the sample to use in interpolation.
123 * @return filter for derivatives from the sample to use in interpolation
124 */
125 public CartesianDerivativesFilter getFilter() {
126 return filter;
127 }
128
129 /** Get output frame for the interpolated instance.
130 * @return output frame for the interpolated instance
131 */
132 public Frame getOutputFrame() {
133 return outputFrame;
134 }
135
136 /**
137 * {@inheritDoc}
138 * <p>
139 * The interpolated instance is created by polynomial Hermite interpolation ensuring velocity remains the exact
140 * derivative of position.
141 * <p>
142 * Note that even if first time derivatives (velocities) from sample can be ignored, the interpolated instance always
143 * includes interpolated derivatives. This feature can be used explicitly to compute these derivatives when it would be
144 * too complex to compute them from an analytical formula: just compute a few sample points from the explicit formula and
145 * set the derivatives to zero in these sample points, then use interpolation to add derivatives consistent with the
146 * positions.
147 */
148 @Override
149 protected FieldAbsolutePVCoordinates<KK> interpolate(final InterpolationData interpolationData) {
150
151 // Get interpolation date
152 final FieldAbsoluteDate<KK> date = interpolationData.getInterpolationDate();
153
154 // Get sample
155 final List<FieldAbsolutePVCoordinates<KK>> sample = interpolationData.getNeighborList();
156
157 // Extract input frame from sample
158 final Frame inputFrame = sample.getFirst().getFrame();
159
160 // Set up an interpolator taking derivatives into account
161 final FieldHermiteInterpolator<KK> interpolator = new FieldHermiteInterpolator<>();
162
163 // Add sample points
164 switch (filter) {
165 case USE_P:
166 // Populate sample with position data, ignoring velocity
167 sample.forEach(pv -> {
168 final FieldVector3D<KK> position = pv.getPosition();
169 interpolator.addSamplePoint(pv.getDate().durationFrom(date),
170 position.toArray());
171 });
172 break;
173 case USE_PV:
174 // Populate sample with position and velocity data
175 sample.forEach(pv -> {
176 final FieldVector3D<KK> position = pv.getPosition();
177 final FieldVector3D<KK> velocity = pv.getVelocity();
178 interpolator.addSamplePoint(pv.getDate().durationFrom(date),
179 position.toArray(), velocity.toArray());
180 });
181 break;
182 case USE_PVA:
183 // Populate sample with position, velocity and acceleration data
184 sample.forEach(pv -> {
185 final FieldVector3D<KK> position = pv.getPosition();
186 final FieldVector3D<KK> velocity = pv.getVelocity();
187 final FieldVector3D<KK> acceleration = pv.getAcceleration();
188 interpolator.addSamplePoint(pv.getDate().durationFrom(date),
189 position.toArray(), velocity.toArray(), acceleration.toArray());
190 });
191 break;
192 default:
193 // this should never happen
194 throw new OrekitInternalError(null);
195 }
196
197 // interpolate
198 final KK zero = interpolationData.getZero();
199 final KK[][] p = interpolator.derivatives(zero, 2);
200
201 // build a new interpolated instance
202 final FieldAbsolutePVCoordinates<KK> interpolated = new FieldAbsolutePVCoordinates<>(inputFrame, date,
203 new FieldVector3D<>(p[0]),
204 new FieldVector3D<>(p[1]),
205 new FieldVector3D<>(p[2]));
206
207 // return the interpolated instance as is
208 if (inputFrame == outputFrame) {
209 return interpolated;
210 }
211
212 // convert to output frame
213 return new FieldAbsolutePVCoordinates<>(outputFrame, interpolated.getPVCoordinates(outputFrame));
214 }
215 }