1 /* Copyright 2002-2023 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.data;
18
19 import java.util.Arrays;
20
21 import org.hipparchus.CalculusFieldElement;
22 import org.hipparchus.util.FastMath;
23 import org.hipparchus.util.FieldSinCos;
24 import org.hipparchus.util.MathArrays;
25 import org.hipparchus.util.SinCos;
26 import org.orekit.errors.OrekitInternalError;
27 import org.orekit.utils.Constants;
28
29 /** Base class for nutation series terms.
30 * @author Luc Maisonobe
31 * @see PoissonSeries
32 */
33 abstract class SeriesTerm {
34
35 /** Coefficients for the sine part. */
36 private double[][] sinCoeff;
37
38 /** Coefficients for the cosine part. */
39 private double[][] cosCoeff;
40
41 /** Simple constructor for the base class.
42 */
43 protected SeriesTerm() {
44 this.sinCoeff = new double[0][0];
45 this.cosCoeff = new double[0][0];
46 }
47
48 /** Get the degree of the function component.
49 * @param index index of the function component (must be less than dimension)
50 * @return degree of the function component
51 */
52 public int getDegree(final int index) {
53 return sinCoeff[index].length - 1;
54 }
55
56 /** Add a pair of values to existing term coefficients.
57 * <p>
58 * Despite it would seem logical to simply set coefficients
59 * rather than add to them, this does not work for some IERS
60 * files. As an example in table 5.3a in IERS conventions 2003,
61 * the coefficients for luni-solar term for 2F+Ω with period
62 * 13.633 days appears twice with different coefficients, as
63 * well as term for 2(F+D+Ω)+l with period 5.643 days, term for
64 * 2(F+D+Ω)-l with period 9.557 days, term for 2(Ω-l') with
65 * period -173.318 days, term for 2D-l with period 31.812 days ...
66 * 35 different duplicated terms have been identified in the
67 * tables 5.3a and 5.3b in IERS conventions 2003.
68 * The coefficients read in lines duplicating a term must be
69 * added together.
70 * </p>
71 * @param index index of the components (will automatically
72 * increase dimension if needed)
73 * @param degree degree of the coefficients, may be negative if
74 * the term does not contribute to component (will automatically
75 * increase {@link #getDegree() degree} of the component if needed)
76 * @param sinID coefficient for the sine part, at index and degree
77 * @param cosID coefficient for the cosine part, at index and degree
78 */
79 public void add(final int index, final int degree,
80 final double sinID, final double cosID) {
81 sinCoeff = extendArray(index, degree, sinCoeff);
82 cosCoeff = extendArray(index, degree, cosCoeff);
83 if (degree >= 0) {
84 sinCoeff[index][degree] += sinID;
85 cosCoeff[index][degree] += cosID;
86 }
87 }
88
89 /** Get a coefficient for the sine part.
90 * @param index index of the function component (must be less than dimension)
91 * @param degree degree of the coefficients
92 * (must be less than {@link #getDegree() degree} for the component)
93 * @return coefficient for the sine part, at index and degree
94 */
95 public double getSinCoeff(final int index, final int degree) {
96 return sinCoeff[index][degree];
97 }
98
99 /** Get a coefficient for the cosine part.
100 * @param index index of the function component (must be less than dimension)
101 * @param degree degree of the coefficients
102 * (must be less than {@link #getDegree() degree} for the component)
103 * @return coefficient for the cosine part, at index and degree
104 */
105 public double getCosCoeff(final int index, final int degree) {
106 return cosCoeff[index][degree];
107 }
108
109 /** Evaluate the value of the series term.
110 * @param elements bodies elements for nutation
111 * @return value of the series term
112 */
113 public double[] value(final BodiesElements elements) {
114
115 // preliminary computation
116 final double tc = elements.getTC();
117 final double a = argument(elements);
118 final SinCos sc = FastMath.sinCos(a);
119
120 // compute each function
121 final double[] values = new double[sinCoeff.length];
122 for (int i = 0; i < values.length; ++i) {
123 double s = 0;
124 double c = 0;
125 for (int j = sinCoeff[i].length - 1; j >= 0; --j) {
126 s = s * tc + sinCoeff[i][j];
127 c = c * tc + cosCoeff[i][j];
128 }
129 values[i] = s * sc.sin() + c * sc.cos();
130 }
131
132 return values;
133
134 }
135
136 /** Evaluate the time derivative of the series term.
137 * @param elements bodies elements for nutation
138 * @return time derivative of the series term
139 */
140 public double[] derivative(final BodiesElements elements) {
141
142 // preliminary computation
143 final double tc = elements.getTC();
144 final double a = argument(elements);
145 final double aDot = argumentDerivative(elements);
146 final SinCos sc = FastMath.sinCos(a);
147
148 // compute each function
149 final double[] derivatives = new double[sinCoeff.length];
150 for (int i = 0; i < derivatives.length; ++i) {
151 double s = 0;
152 double c = 0;
153 double sDot = 0;
154 double cDot = 0;
155 if (sinCoeff[i].length > 0) {
156 for (int j = sinCoeff[i].length - 1; j > 0; --j) {
157 s = s * tc + sinCoeff[i][j];
158 c = c * tc + cosCoeff[i][j];
159 sDot = sDot * tc + j * sinCoeff[i][j];
160 cDot = cDot * tc + j * cosCoeff[i][j];
161 }
162 s = s * tc + sinCoeff[i][0];
163 c = c * tc + cosCoeff[i][0];
164 sDot /= Constants.JULIAN_CENTURY;
165 cDot /= Constants.JULIAN_CENTURY;
166 }
167 derivatives[i] = (sDot - c * aDot) * sc.sin() + (cDot + s * aDot) * sc.cos();
168 }
169
170 return derivatives;
171
172 }
173
174 /** Compute the argument for the current date.
175 * @param elements luni-solar and planetary elements for the current date
176 * @return current value of the argument
177 */
178 protected abstract double argument(BodiesElements elements);
179
180 /** Compute the time derivative of the argument for the current date.
181 * @param elements luni-solar and planetary elements for the current date
182 * @return current time derivative of the argument
183 */
184 protected abstract double argumentDerivative(BodiesElements elements);
185
186 /** Evaluate the value of the series term.
187 * @param elements bodies elements for nutation
188 * @param <T> the type of the field elements
189 * @return value of the series term
190 */
191 public <T extends CalculusFieldElement<T>> T[] value(final FieldBodiesElements<T> elements) {
192
193 // preliminary computation
194 final T tc = elements.getTC();
195 final T a = argument(elements);
196 final FieldSinCos<T> sc = FastMath.sinCos(a);
197
198 // compute each function
199 final T[] values = MathArrays.buildArray(tc.getField(), sinCoeff.length);
200 for (int i = 0; i < values.length; ++i) {
201 T s = tc.getField().getZero();
202 T c = tc.getField().getZero();
203 for (int j = sinCoeff[i].length - 1; j >= 0; --j) {
204 s = s.multiply(tc).add(sinCoeff[i][j]);
205 c = c.multiply(tc).add(cosCoeff[i][j]);
206 }
207 values[i] = s.multiply(sc.sin()).add(c.multiply(sc.cos()));
208 }
209
210 return values;
211
212 }
213
214 /** Evaluate the time derivative of the series term.
215 * @param elements bodies elements for nutation
216 * @param <T> the type of the field elements
217 * @return time derivative of the series term
218 */
219 public <T extends CalculusFieldElement<T>> T[] derivative(final FieldBodiesElements<T> elements) {
220
221 // preliminary computation
222 final T tc = elements.getTC();
223 final T a = argument(elements);
224 final T aDot = argumentDerivative(elements);
225 final FieldSinCos<T> sc = FastMath.sinCos(a);
226
227 // compute each function
228 final T[] derivatives = MathArrays.buildArray(tc.getField(), sinCoeff.length);
229 for (int i = 0; i < derivatives.length; ++i) {
230 T s = tc.getField().getZero();
231 T c = tc.getField().getZero();
232 T sDot = tc.getField().getZero();
233 T cDot = tc.getField().getZero();
234 if (sinCoeff[i].length > 0) {
235 for (int j = sinCoeff[i].length - 1; j > 0; --j) {
236 s = s.multiply(tc).add(sinCoeff[i][j]);
237 c = c.multiply(tc).add(cosCoeff[i][j]);
238 sDot = sDot.multiply(tc).add(j * sinCoeff[i][j]);
239 cDot = cDot.multiply(tc).add(j * cosCoeff[i][j]);
240 }
241 s = s.multiply(tc).add(sinCoeff[i][0]);
242 c = c.multiply(tc).add(cosCoeff[i][0]);
243 sDot = sDot.divide(Constants.JULIAN_CENTURY);
244 cDot = cDot.divide(Constants.JULIAN_CENTURY);
245 }
246 derivatives[i] = sDot.subtract(c.multiply(aDot)).multiply(sc.sin()).
247 add(cDot.add(s.multiply(aDot)).multiply(sc.cos()));
248 }
249
250 return derivatives;
251
252 }
253
254 /** Compute the argument for the current date.
255 * @param elements luni-solar and planetary elements for the current date
256 * @param <T> the type of the field elements
257 * @return current value of the argument
258 */
259 protected abstract <T extends CalculusFieldElement<T>> T argument(FieldBodiesElements<T> elements);
260
261 /** Compute the time derivative of the argument for the current date.
262 * @param elements luni-solar and planetary elements for the current date
263 * @param <T> the type of the field elements
264 * @return current time derivative of the argument
265 */
266 protected abstract <T extends CalculusFieldElement<T>> T argumentDerivative(FieldBodiesElements<T> elements);
267
268 /** Factory method for building the appropriate object.
269 * <p>The method checks the null coefficients and build an instance
270 * of an appropriate type to avoid too many unnecessary multiplications
271 * by zero coefficients.</p>
272 * @param cGamma coefficient for γ = GMST + π tide parameter
273 * @param cL coefficient for mean anomaly of the Moon
274 * @param cLPrime coefficient for mean anomaly of the Sun
275 * @param cF coefficient for L - Ω where L is the mean longitude of the Moon
276 * @param cD coefficient for mean elongation of the Moon from the Sun
277 * @param cOmega coefficient for mean longitude of the ascending node of the Moon
278 * @param cMe coefficient for mean Mercury longitude
279 * @param cVe coefficient for mean Venus longitude
280 * @param cE coefficient for mean Earth longitude
281 * @param cMa coefficient for mean Mars longitude
282 * @param cJu coefficient for mean Jupiter longitude
283 * @param cSa coefficient for mean Saturn longitude
284 * @param cUr coefficient for mean Uranus longitude
285 * @param cNe coefficient for mean Neptune longitude
286 * @param cPa coefficient for general accumulated precession in longitude
287 * @return a nutation serie term instance well suited for the set of coefficients
288 */
289 public static SeriesTerm buildTerm(final int cGamma,
290 final int cL, final int cLPrime, final int cF,
291 final int cD, final int cOmega,
292 final int cMe, final int cVe, final int cE,
293 final int cMa, final int cJu, final int cSa,
294 final int cUr, final int cNe, final int cPa) {
295 if (cGamma == 0 && cL == 0 && cLPrime == 0 && cF == 0 && cD == 0 && cOmega == 0) {
296 return new PlanetaryTerm(cMe, cVe, cE, cMa, cJu, cSa, cUr, cNe, cPa);
297 } else if (cGamma == 0 &&
298 cMe == 0 && cVe == 0 && cE == 0 && cMa == 0 && cJu == 0 &&
299 cSa == 0 && cUr == 0 && cNe == 0 && cPa == 0) {
300 return new LuniSolarTerm(cL, cLPrime, cF, cD, cOmega);
301 } else if (cGamma != 0 &&
302 cMe == 0 && cVe == 0 && cE == 0 && cMa == 0 && cJu == 0 &&
303 cSa == 0 && cUr == 0 && cNe == 0 && cPa == 0) {
304 return new TideTerm(cGamma, cL, cLPrime, cF, cD, cOmega);
305 } else if (cGamma == 0 && cLPrime == 0 && cUr == 0 && cNe == 0 && cPa == 0) {
306 return new NoFarPlanetsTerm(cL, cF, cD, cOmega,
307 cMe, cVe, cE, cMa, cJu, cSa);
308 } else if (cGamma == 0) {
309 return new GeneralTerm(cL, cLPrime, cF, cD, cOmega,
310 cMe, cVe, cE, cMa, cJu, cSa, cUr, cNe, cPa);
311 } else {
312 throw new OrekitInternalError(null);
313 }
314
315 }
316
317 /** Extend an array to hold at least index and degree.
318 * @param index index of the function
319 * @param degree degree of the coefficients
320 * @param array to extend
321 * @return extended array
322 */
323 private static double[][] extendArray(final int index, final int degree,
324 final double[][] array) {
325
326 // extend the number of rows if needed
327 final double[][] extended;
328 if (array.length > index) {
329 extended = array;
330 } else {
331 final int rows = index + 1;
332 extended = new double[rows][];
333 System.arraycopy(array, 0, extended, 0, array.length);
334 Arrays.fill(extended, array.length, index + 1, new double[0]);
335 }
336
337 // extend the number of elements in the row if needed
338 extended[index] = extendArray(degree, extended[index]);
339
340 return extended;
341
342 }
343
344 /** Extend an array to hold at least degree.
345 * @param degree degree of the coefficients
346 * @param array to extend
347 * @return extended array
348 */
349 private static double[] extendArray(final int degree, final double[] array) {
350 // extend the number of elements if needed
351 if (array.length > degree) {
352 return array;
353 } else {
354 final double[] extended = new double[degree + 1];
355 System.arraycopy(array, 0, extended, 0, array.length);
356 return extended;
357 }
358 }
359
360 }