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