1 /* Copyright 2002-2025 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.bodies;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
21 import org.hipparchus.geometry.euclidean.threed.Vector3D;
22 import org.orekit.time.AbsoluteDate;
23 import org.orekit.time.FieldAbsoluteDate;
24 import org.orekit.time.TimeScale;
25 import org.orekit.time.TimeStamped;
26 import org.orekit.utils.FieldPVCoordinates;
27 import org.orekit.utils.PVCoordinates;
28
29
30 /** Position-Velocity model based on Chebyshev polynomials.
31 * <p>This class represent the most basic element of the piecewise ephemerides
32 * for solar system bodies like JPL DE 405 ephemerides.</p>
33 * @see JPLEphemeridesLoader
34 * @author Luc Maisonobe
35 */
36 class PosVelChebyshev implements TimeStamped {
37
38 /** Time scale in which the ephemeris is defined. */
39 private final TimeScale timeScale;
40
41 /** Start of the validity range of the instance. */
42 private final AbsoluteDate start;
43
44 /** Duration of validity range of the instance. */
45 private final double duration;
46
47 /** Chebyshev polynomials coefficients for the X component. */
48 private final double[] xCoeffs;
49
50 /** Chebyshev polynomials coefficients for the Y component. */
51 private final double[] yCoeffs;
52
53 /** Chebyshev polynomials coefficients for the Z component. */
54 private final double[] zCoeffs;
55
56 /** Velocity scale for internal use. */
57 private final double vScale;
58 /** Acceleration scale for internal use. */
59 private final double aScale;
60
61 /** Simple constructor.
62 * @param start start of the validity range of the instance
63 * @param timeScale time scale in which the ephemeris is defined
64 * @param duration duration of the validity range of the instance
65 * @param xCoeffs Chebyshev polynomials coefficients for the X component
66 * (a reference to the array will be stored in the instance)
67 * @param yCoeffs Chebyshev polynomials coefficients for the Y component
68 * (a reference to the array will be stored in the instance)
69 * @param zCoeffs Chebyshev polynomials coefficients for the Z component
70 * (a reference to the array will be stored in the instance)
71 */
72 PosVelChebyshev(final AbsoluteDate start, final TimeScale timeScale, final double duration,
73 final double[] xCoeffs, final double[] yCoeffs, final double[] zCoeffs) {
74 this.start = start;
75 this.timeScale = timeScale;
76 this.duration = duration;
77 this.xCoeffs = xCoeffs;
78 this.yCoeffs = yCoeffs;
79 this.zCoeffs = zCoeffs;
80 this.vScale = 2 / duration;
81 this.aScale = this.vScale * this.vScale;
82 }
83
84 /** {@inheritDoc} */
85 public AbsoluteDate getDate() {
86 return start;
87 }
88
89 /** Compute value of Chebyshev's polynomial independent variable.
90 * @param date date
91 * @return double independent variable value
92 */
93 private double computeValueIndependentVariable(final AbsoluteDate date) {
94 return (2 * date.offsetFrom(start, timeScale) - duration) / duration;
95 }
96
97 /** Compute value of Chebyshev's polynomial independent variable.
98 * @param date date
99 * @param <T> type of the field elements
100 * @return <T> independent variable value
101 */
102 private <T extends CalculusFieldElement<T>> T computeValueIndependentVariable(final FieldAbsoluteDate<T> date) {
103 return date.offsetFrom(new FieldAbsoluteDate<>(date.getField(), start), timeScale).multiply(2).subtract(duration).divide(duration);
104 }
105
106 /** Check if a date is in validity range.
107 * @param date date to check
108 * @return true if date is in validity range
109 */
110 public boolean inRange(final AbsoluteDate date) {
111 final double dt = date.offsetFrom(start, timeScale);
112 return dt >= -0.001 && dt <= duration + 0.001;
113 }
114
115 /** Get the position at a specified date.
116 * @param date date at which position is requested
117 * @return position at specified date
118 */
119 Vector3D getPosition(final AbsoluteDate date) {
120
121 // normalize date
122 final double t = computeValueIndependentVariable(date);
123 final double twoT = 2 * t;
124
125 // initialize Chebyshev polynomials recursion
126 double pKm1 = 1;
127 double pK = t;
128 double xP = xCoeffs[0];
129 double yP = yCoeffs[0];
130 double zP = zCoeffs[0];
131
132 // combine polynomials by applying coefficients
133 for (int k = 1; k < xCoeffs.length; ++k) {
134
135 // consider last computed polynomials on position
136 xP += xCoeffs[k] * pK;
137 yP += yCoeffs[k] * pK;
138 zP += zCoeffs[k] * pK;
139
140 // compute next Chebyshev polynomial value
141 final double pKm2 = pKm1;
142 pKm1 = pK;
143 pK = twoT * pKm1 - pKm2;
144
145 }
146
147 return new Vector3D(xP, yP, zP);
148 }
149
150 /** Get the position at a specified date.
151 * @param date date at which position is requested
152 * @param <T> type of the field elements
153 * @return position at specified date
154 */
155 <T extends CalculusFieldElement<T>> FieldVector3D<T> getPosition(final FieldAbsoluteDate<T> date) {
156
157 final T zero = date.getField().getZero();
158 final T one = date.getField().getOne();
159
160 // normalize date
161 final T t = computeValueIndependentVariable(date);
162 final T twoT = t.add(t);
163
164 // initialize Chebyshev polynomials recursion
165 T pKm1 = one;
166 T pK = t;
167 T xP = zero.newInstance(xCoeffs[0]);
168 T yP = zero.newInstance(yCoeffs[0]);
169 T zP = zero.newInstance(zCoeffs[0]);
170
171 // combine polynomials by applying coefficients
172 for (int k = 1; k < xCoeffs.length; ++k) {
173
174 // consider last computed polynomials on position
175 xP = xP.add(pK.multiply(xCoeffs[k]));
176 yP = yP.add(pK.multiply(yCoeffs[k]));
177 zP = zP.add(pK.multiply(zCoeffs[k]));
178
179 // compute next Chebyshev polynomial value
180 final T pKm2 = pKm1;
181 pKm1 = pK;
182 pK = twoT.multiply(pKm1).subtract(pKm2);
183
184 }
185
186 return new FieldVector3D<>(xP, yP, zP);
187
188 }
189
190 /** Get the position-velocity-acceleration at a specified date.
191 * @param date date at which position-velocity-acceleration is requested
192 * @return position-velocity-acceleration at specified date
193 */
194 PVCoordinates getPositionVelocityAcceleration(final AbsoluteDate date) {
195
196 // normalize date
197 final double t = computeValueIndependentVariable(date);
198 final double twoT = 2 * t;
199
200 // initialize Chebyshev polynomials recursion
201 double pKm1 = 1;
202 double pK = t;
203 double xP = xCoeffs[0];
204 double yP = yCoeffs[0];
205 double zP = zCoeffs[0];
206
207 // initialize Chebyshev polynomials derivatives recursion
208 double qKm1 = 0;
209 double qK = 1;
210 double xV = 0;
211 double yV = 0;
212 double zV = 0;
213
214 // initialize Chebyshev polynomials second derivatives recursion
215 double rKm1 = 0;
216 double rK = 0;
217 double xA = 0;
218 double yA = 0;
219 double zA = 0;
220
221 // combine polynomials by applying coefficients
222 for (int k = 1; k < xCoeffs.length; ++k) {
223
224 // consider last computed polynomials on position
225 xP += xCoeffs[k] * pK;
226 yP += yCoeffs[k] * pK;
227 zP += zCoeffs[k] * pK;
228
229 // consider last computed polynomials on velocity
230 xV += xCoeffs[k] * qK;
231 yV += yCoeffs[k] * qK;
232 zV += zCoeffs[k] * qK;
233
234 // consider last computed polynomials on acceleration
235 xA += xCoeffs[k] * rK;
236 yA += yCoeffs[k] * rK;
237 zA += zCoeffs[k] * rK;
238
239 // compute next Chebyshev polynomial value
240 final double pKm2 = pKm1;
241 pKm1 = pK;
242 pK = twoT * pKm1 - pKm2;
243
244 // compute next Chebyshev polynomial derivative
245 final double qKm2 = qKm1;
246 qKm1 = qK;
247 qK = twoT * qKm1 + 2 * pKm1 - qKm2;
248
249 // compute next Chebyshev polynomial second derivative
250 final double rKm2 = rKm1;
251 rKm1 = rK;
252 rK = twoT * rKm1 + 4 * qKm1 - rKm2;
253
254 }
255
256 return new PVCoordinates(new Vector3D(xP, yP, zP),
257 new Vector3D(xV * vScale, yV * vScale, zV * vScale),
258 new Vector3D(xA * aScale, yA * aScale, zA * aScale));
259
260 }
261
262 /** Get the position-velocity-acceleration at a specified date.
263 * @param date date at which position-velocity-acceleration is requested
264 * @param <T> type of the field elements
265 * @return position-velocity-acceleration at specified date
266 */
267 <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getPositionVelocityAcceleration(final FieldAbsoluteDate<T> date) {
268
269 final T zero = date.getField().getZero();
270 final T one = date.getField().getOne();
271
272 // normalize date
273 final T t = computeValueIndependentVariable(date);
274 final T twoT = t.add(t);
275
276 // initialize Chebyshev polynomials recursion
277 T pKm1 = one;
278 T pK = t;
279 T xP = zero.newInstance(xCoeffs[0]);
280 T yP = zero.newInstance(yCoeffs[0]);
281 T zP = zero.newInstance(zCoeffs[0]);
282
283 // initialize Chebyshev polynomials derivatives recursion
284 T qKm1 = zero;
285 T qK = one;
286 T xV = zero;
287 T yV = zero;
288 T zV = zero;
289
290 // initialize Chebyshev polynomials second derivatives recursion
291 T rKm1 = zero;
292 T rK = zero;
293 T xA = zero;
294 T yA = zero;
295 T zA = zero;
296
297 // combine polynomials by applying coefficients
298 for (int k = 1; k < xCoeffs.length; ++k) {
299
300 // consider last computed polynomials on position
301 xP = xP.add(pK.multiply(xCoeffs[k]));
302 yP = yP.add(pK.multiply(yCoeffs[k]));
303 zP = zP.add(pK.multiply(zCoeffs[k]));
304
305 // consider last computed polynomials on velocity
306 xV = xV.add(qK.multiply(xCoeffs[k]));
307 yV = yV.add(qK.multiply(yCoeffs[k]));
308 zV = zV.add(qK.multiply(zCoeffs[k]));
309
310 // consider last computed polynomials on acceleration
311 xA = xA.add(rK.multiply(xCoeffs[k]));
312 yA = yA.add(rK.multiply(yCoeffs[k]));
313 zA = zA.add(rK.multiply(zCoeffs[k]));
314
315 // compute next Chebyshev polynomial value
316 final T pKm2 = pKm1;
317 pKm1 = pK;
318 pK = twoT.multiply(pKm1).subtract(pKm2);
319
320 // compute next Chebyshev polynomial derivative
321 final T qKm2 = qKm1;
322 qKm1 = qK;
323 qK = twoT.multiply(qKm1).add(pKm1.multiply(2)).subtract(qKm2);
324
325 // compute next Chebyshev polynomial second derivative
326 final T rKm2 = rKm1;
327 rKm1 = rK;
328 rK = twoT.multiply(rKm1).add(qKm1.multiply(4)).subtract(rKm2);
329
330 }
331
332 return new FieldPVCoordinates<>(new FieldVector3D<>(xP, yP, zP),
333 new FieldVector3D<>(xV.multiply(vScale), yV.multiply(vScale), zV.multiply(vScale)),
334 new FieldVector3D<>(xA.multiply(aScale), yA.multiply(aScale), zA.multiply(aScale)));
335
336 }
337
338 }