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