1 /* Copyright 2002-2022 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 /** 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 }
81
82 /** {@inheritDoc} */
83 public AbsoluteDate getDate() {
84 return start;
85 }
86
87 /** Check if a date is in validity range.
88 * @param date date to check
89 * @return true if date is in validity range
90 */
91 public boolean inRange(final AbsoluteDate date) {
92 final double dt = date.offsetFrom(start, timeScale);
93 return dt >= -0.001 && dt <= duration + 0.001;
94 }
95
96 /** Get the position-velocity-acceleration at a specified date.
97 * @param date date at which position-velocity-acceleration is requested
98 * @return position-velocity-acceleration at specified date
99 */
100 public PVCoordinates getPositionVelocityAcceleration(final AbsoluteDate date) {
101
102 // normalize date
103 final double t = (2 * date.offsetFrom(start, timeScale) - duration) / duration;
104 final double twoT = 2 * t;
105
106 // initialize Chebyshev polynomials recursion
107 double pKm1 = 1;
108 double pK = t;
109 double xP = xCoeffs[0];
110 double yP = yCoeffs[0];
111 double zP = zCoeffs[0];
112
113 // initialize Chebyshev polynomials derivatives recursion
114 double qKm1 = 0;
115 double qK = 1;
116 double xV = 0;
117 double yV = 0;
118 double zV = 0;
119
120 // initialize Chebyshev polynomials second derivatives recursion
121 double rKm1 = 0;
122 double rK = 0;
123 double xA = 0;
124 double yA = 0;
125 double zA = 0;
126
127 // combine polynomials by applying coefficients
128 for (int k = 1; k < xCoeffs.length; ++k) {
129
130 // consider last computed polynomials on position
131 xP += xCoeffs[k] * pK;
132 yP += yCoeffs[k] * pK;
133 zP += zCoeffs[k] * pK;
134
135 // consider last computed polynomials on velocity
136 xV += xCoeffs[k] * qK;
137 yV += yCoeffs[k] * qK;
138 zV += zCoeffs[k] * qK;
139
140 // consider last computed polynomials on acceleration
141 xA += xCoeffs[k] * rK;
142 yA += yCoeffs[k] * rK;
143 zA += zCoeffs[k] * rK;
144
145 // compute next Chebyshev polynomial value
146 final double pKm2 = pKm1;
147 pKm1 = pK;
148 pK = twoT * pKm1 - pKm2;
149
150 // compute next Chebyshev polynomial derivative
151 final double qKm2 = qKm1;
152 qKm1 = qK;
153 qK = twoT * qKm1 + 2 * pKm1 - qKm2;
154
155 // compute next Chebyshev polynomial second derivative
156 final double rKm2 = rKm1;
157 rKm1 = rK;
158 rK = twoT * rKm1 + 4 * qKm1 - rKm2;
159
160 }
161
162 final double vScale = 2 / duration;
163 final double aScale = vScale * vScale;
164 return new PVCoordinates(new Vector3D(xP, yP, zP),
165 new Vector3D(xV * vScale, yV * vScale, zV * vScale),
166 new Vector3D(xA * aScale, yA * aScale, zA * aScale));
167
168 }
169
170 /** Get the position-velocity-acceleration at a specified date.
171 * @param date date at which position-velocity-acceleration is requested
172 * @param <T> type fo the field elements
173 * @return position-velocity-acceleration at specified date
174 */
175 public <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getPositionVelocityAcceleration(final FieldAbsoluteDate<T> date) {
176
177 final T zero = date.getField().getZero();
178 final T one = date.getField().getOne();
179
180 // normalize date
181 final T t = date.offsetFrom(new FieldAbsoluteDate<>(date.getField(), start), timeScale).multiply(2).subtract(duration).divide(duration);
182 final T twoT = t.add(t);
183
184 // initialize Chebyshev polynomials recursion
185 T pKm1 = one;
186 T pK = t;
187 T xP = zero.add(xCoeffs[0]);
188 T yP = zero.add(yCoeffs[0]);
189 T zP = zero.add(zCoeffs[0]);
190
191 // initialize Chebyshev polynomials derivatives recursion
192 T qKm1 = zero;
193 T qK = one;
194 T xV = zero;
195 T yV = zero;
196 T zV = zero;
197
198 // initialize Chebyshev polynomials second derivatives recursion
199 T rKm1 = zero;
200 T rK = zero;
201 T xA = zero;
202 T yA = zero;
203 T zA = zero;
204
205 // combine polynomials by applying coefficients
206 for (int k = 1; k < xCoeffs.length; ++k) {
207
208 // consider last computed polynomials on position
209 xP = xP.add(pK.multiply(xCoeffs[k]));
210 yP = yP.add(pK.multiply(yCoeffs[k]));
211 zP = zP.add(pK.multiply(zCoeffs[k]));
212
213 // consider last computed polynomials on velocity
214 xV = xV.add(qK.multiply(xCoeffs[k]));
215 yV = yV.add(qK.multiply(yCoeffs[k]));
216 zV = zV.add(qK.multiply(zCoeffs[k]));
217
218 // consider last computed polynomials on acceleration
219 xA = xA.add(rK.multiply(xCoeffs[k]));
220 yA = yA.add(rK.multiply(yCoeffs[k]));
221 zA = zA.add(rK.multiply(zCoeffs[k]));
222
223 // compute next Chebyshev polynomial value
224 final T pKm2 = pKm1;
225 pKm1 = pK;
226 pK = twoT.multiply(pKm1).subtract(pKm2);
227
228 // compute next Chebyshev polynomial derivative
229 final T qKm2 = qKm1;
230 qKm1 = qK;
231 qK = twoT.multiply(qKm1).add(pKm1.multiply(2)).subtract(qKm2);
232
233 // compute next Chebyshev polynomial second derivative
234 final T rKm2 = rKm1;
235 rKm1 = rK;
236 rK = twoT.multiply(rKm1).add(qKm1.multiply(4)).subtract(rKm2);
237
238 }
239
240 final double vScale = 2 / duration;
241 final double aScale = vScale * vScale;
242 return new FieldPVCoordinates<>(new FieldVector3D<>(xP, yP, zP),
243 new FieldVector3D<>(xV.multiply(vScale), yV.multiply(vScale), zV.multiply(vScale)),
244 new FieldVector3D<>(xA.multiply(aScale), yA.multiply(aScale), zA.multiply(aScale)));
245
246 }
247
248 }