1   /* Copyright 2002-2024 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 }