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