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.HashMap;
20 import java.util.Map;
21
22 import org.hipparchus.CalculusFieldElement;
23 import org.hipparchus.util.MathArrays;
24 import org.hipparchus.util.MathUtils;
25 import org.hipparchus.util.MathUtils.SumAndResidual;
26
27 /**
28 * Class representing a Poisson series for nutation or ephemeris computations.
29 * <p>
30 * A Poisson series is composed of a time polynomial part and a non-polynomial
31 * part which consist in summation series. The {@link SeriesTerm series terms}
32 * are harmonic functions (combination of sines and cosines) of polynomial
33 * <em>arguments</em>. The polynomial arguments are combinations of luni-solar or
34 * planetary {@link BodiesElements elements}.
35 * </p>
36 * @author Luc Maisonobe
37 * @see PoissonSeriesParser
38 * @see SeriesTerm
39 * @see PolynomialNutation
40 */
41 public class PoissonSeries {
42
43 /** Polynomial part. */
44 private final PolynomialNutation polynomial;
45
46 /** Non-polynomial series. */
47 private final Map<Long, SeriesTerm> series;
48
49 /** Build a Poisson series from an IERS table file.
50 * @param polynomial polynomial part (may be null)
51 * @param series non-polynomial part
52 */
53 public PoissonSeries(final PolynomialNutation polynomial, final Map<Long, SeriesTerm> series) {
54 this.polynomial = polynomial;
55 this.series = series;
56 }
57
58 /** Get the polynomial part of the series.
59 * @return polynomial part of the series.
60 */
61 public PolynomialNutation getPolynomial() {
62 return polynomial;
63 }
64
65 /** Get the number of different terms in the non-polynomial part.
66 * @return number of different terms in the non-polynomial part
67 */
68 public int getNonPolynomialSize() {
69 return series.size();
70 }
71
72 /** Evaluate the value of the series.
73 * @param elements bodies elements for nutation
74 * @return value of the series
75 */
76 public double value(final BodiesElements elements) {
77
78 // polynomial part
79 final double p = polynomial.value(elements.getTC());
80
81 // non-polynomial part
82 double npHigh = 0;
83 double npLow = 0;
84 for (final Map.Entry<Long, SeriesTerm> entry : series.entrySet()) {
85 final double v = entry.getValue().value(elements)[0];
86 // Use 2Sum for high precision.
87 final SumAndResidual sumAndResidual = MathUtils.twoSum(npHigh, v);
88 npHigh = sumAndResidual.getSum();
89 npLow += sumAndResidual.getResidual();
90 }
91
92 // add the polynomial and the non-polynomial parts
93 return p + (npHigh + npLow);
94
95 }
96
97 /** Evaluate the value of the series.
98 * @param elements bodies elements for nutation
99 * @param <T> type of the field elements
100 * @return value of the series
101 */
102 public <T extends CalculusFieldElement<T>> T value(final FieldBodiesElements<T> elements) {
103
104 // polynomial part
105 final T tc = elements.getTC();
106 final T p = polynomial.value(tc);
107
108 // non-polynomial part
109 T sum = tc.getField().getZero();
110 for (final Map.Entry<Long, SeriesTerm> entry : series.entrySet()) {
111 sum = sum.add(entry.getValue().value(elements)[0]);
112 }
113
114 // add the polynomial and the non-polynomial parts
115 return p.add(sum);
116
117 }
118
119 /** This interface represents a fast evaluator for Poisson series.
120 * @see PoissonSeries#compile(PoissonSeries...)
121 * @since 6.1
122 */
123 public interface CompiledSeries {
124
125 /** Evaluate a set of Poisson series.
126 * @param elements bodies elements for nutation
127 * @return value of the series
128 */
129 double[] value(BodiesElements elements);
130
131 /** Evaluate time derivative of a set of Poisson series.
132 * @param elements bodies elements for nutation
133 * @return time derivative of the series
134 */
135 double[] derivative(BodiesElements elements);
136
137 /** Evaluate a set of Poisson series.
138 * @param elements bodies elements for nutation
139 * @param <S> the type of the field elements
140 * @return value of the series
141 */
142 <S extends CalculusFieldElement<S>> S[] value(FieldBodiesElements<S> elements);
143
144 /** Evaluate time derivative of a set of Poisson series.
145 * @param elements bodies elements for nutation
146 * @param <S> the type of the field elements
147 * @return time derivative of the series
148 */
149 <S extends CalculusFieldElement<S>> S[] derivative(FieldBodiesElements<S> elements);
150
151 }
152
153 /** Join several nutation series, for fast simultaneous evaluation.
154 * @param poissonSeries Poisson series to join
155 * @return a single function that evaluates all series together
156 * @since 6.1
157 */
158 @SafeVarargs
159 public static CompiledSeries compile(final PoissonSeries... poissonSeries) {
160
161 // store all polynomials
162 final PolynomialNutation[] polynomials = new PolynomialNutation[poissonSeries.length];
163 for (int i = 0; i < polynomials.length; ++i) {
164 polynomials[i] = poissonSeries[i].polynomial;
165 }
166
167 // gather all series terms
168 final Map<Long, SeriesTerm> joinedMap = new HashMap<Long, SeriesTerm>();
169 for (final PoissonSeries ps : poissonSeries) {
170 for (Map.Entry<Long, SeriesTerm> entry : ps.series.entrySet()) {
171 final long key = entry.getKey();
172 if (!joinedMap.containsKey(key)) {
173
174 // retrieve all Delaunay and planetary multipliers from the key
175 final int[] m = NutationCodec.decode(key);
176
177 // prepare a new term, ready to handle the required dimension
178 final SeriesTerm term =
179 SeriesTerm.buildTerm(m[0],
180 m[1], m[2], m[3], m[4], m[5],
181 m[6], m[7], m[8], m[9], m[10], m[11], m[12], m[13], m[14]);
182 term.add(poissonSeries.length - 1, -1, Double.NaN, Double.NaN);
183
184 // store it
185 joinedMap.put(key, term);
186
187 }
188 }
189 }
190
191 // join series by sharing terms, in order to speed up evaluation
192 // which is dominated by the computation of sine/cosine in each term
193 for (int i = 0; i < poissonSeries.length; ++i) {
194 for (final Map.Entry<Long, SeriesTerm> entry : poissonSeries[i].series.entrySet()) {
195 final SeriesTerm singleTerm = entry.getValue();
196 final SeriesTerm joinedTerm = joinedMap.get(entry.getKey());
197 for (int degree = 0; degree <= singleTerm.getDegree(0); ++degree) {
198 joinedTerm.add(i, degree,
199 singleTerm.getSinCoeff(0, degree),
200 singleTerm.getCosCoeff(0, degree));
201 }
202 }
203 }
204
205 // use a single array for faster access
206 final SeriesTerm[] joinedTerms = new SeriesTerm[joinedMap.size()];
207 int index = 0;
208 for (final Map.Entry<Long, SeriesTerm> entry : joinedMap.entrySet()) {
209 joinedTerms[index++] = entry.getValue();
210 }
211
212 return new CompiledSeries() {
213
214 /** {@inheritDoc} */
215 @Override
216 public double[] value(final BodiesElements elements) {
217
218 // non-polynomial part
219 final double[] npHigh = new double[polynomials.length];
220 final double[] npLow = new double[polynomials.length];
221 for (final SeriesTerm term : joinedTerms) {
222 final double[] termValue = term.value(elements);
223 for (int i = 0; i < termValue.length; ++i) {
224 // Use 2Sum for high precision.
225 final SumAndResidual sumAndResidual = MathUtils.twoSum(npHigh[i], termValue[i]);
226 npHigh[i] = sumAndResidual.getSum();
227 npLow[i] += sumAndResidual.getResidual();
228 }
229 }
230
231 // add residual and polynomial part
232 for (int i = 0; i < npHigh.length; ++i) {
233 npHigh[i] += npLow[i] + polynomials[i].value(elements.getTC());
234 }
235 return npHigh;
236
237 }
238
239 /** {@inheritDoc} */
240 @Override
241 public double[] derivative(final BodiesElements elements) {
242
243 // non-polynomial part
244 final double[] v = new double[polynomials.length];
245 for (final SeriesTerm term : joinedTerms) {
246 final double[] termDerivative = term.derivative(elements);
247 for (int i = 0; i < termDerivative.length; ++i) {
248 v[i] += termDerivative[i];
249 }
250 }
251
252 // add polynomial part
253 for (int i = 0; i < v.length; ++i) {
254 v[i] += polynomials[i].derivative(elements.getTC());
255 }
256 return v;
257
258 }
259
260 /** {@inheritDoc} */
261 @Override
262 public <S extends CalculusFieldElement<S>> S[] value(final FieldBodiesElements<S> elements) {
263
264 // non-polynomial part
265 final S[] v = MathArrays.buildArray(elements.getTC().getField(), polynomials.length);
266 for (final SeriesTerm term : joinedTerms) {
267 final S[] termValue = term.value(elements);
268 for (int i = 0; i < termValue.length; ++i) {
269 v[i] = v[i].add(termValue[i]);
270 }
271 }
272
273 // add polynomial part
274 final S tc = elements.getTC();
275 for (int i = 0; i < v.length; ++i) {
276 v[i] = v[i].add(polynomials[i].value(tc));
277 }
278 return v;
279
280 }
281
282 /** {@inheritDoc} */
283 @Override
284 public <S extends CalculusFieldElement<S>> S[] derivative(final FieldBodiesElements<S> elements) {
285
286 // non-polynomial part
287 final S[] v = MathArrays.buildArray(elements.getTC().getField(), polynomials.length);
288 for (final SeriesTerm term : joinedTerms) {
289 final S[] termDerivative = term.derivative(elements);
290 for (int i = 0; i < termDerivative.length; ++i) {
291 v[i] = v[i].add(termDerivative[i]);
292 }
293 }
294
295 // add polynomial part
296 final S tc = elements.getTC();
297 for (int i = 0; i < v.length; ++i) {
298 v[i] = v[i].add(polynomials[i].derivative(tc));
299 }
300 return v;
301
302 }
303
304 };
305
306 }
307
308 }