PoissonSeries.java

  1. /* Copyright 2002-2022 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. import java.util.HashMap;
  19. import java.util.Map;

  20. import org.hipparchus.CalculusFieldElement;
  21. import org.hipparchus.util.MathArrays;
  22. import org.hipparchus.util.MathUtils;
  23. import org.hipparchus.util.MathUtils.SumAndResidual;

  24. /**
  25.  * Class representing a Poisson series for nutation or ephemeris computations.
  26.  * <p>
  27.  * A Poisson series is composed of a time polynomial part and a non-polynomial
  28.  * part which consist in summation series. The {@link SeriesTerm series terms}
  29.  * are harmonic functions (combination of sines and cosines) of polynomial
  30.  * <em>arguments</em>. The polynomial arguments are combinations of luni-solar or
  31.  * planetary {@link BodiesElements elements}.
  32.  * </p>
  33.  * @author Luc Maisonobe
  34.  * @see PoissonSeriesParser
  35.  * @see SeriesTerm
  36.  * @see PolynomialNutation
  37.  */
  38. public class PoissonSeries {

  39.     /** Polynomial part. */
  40.     private final PolynomialNutation polynomial;

  41.     /** Non-polynomial series. */
  42.     private final Map<Long, SeriesTerm> series;

  43.     /** Build a Poisson series from an IERS table file.
  44.      * @param polynomial polynomial part (may be null)
  45.      * @param series non-polynomial part
  46.      */
  47.     public PoissonSeries(final PolynomialNutation polynomial, final Map<Long, SeriesTerm> series) {
  48.         this.polynomial = polynomial;
  49.         this.series     = series;
  50.     }

  51.     /** Get the polynomial part of the series.
  52.      * @return polynomial part of the series.
  53.      */
  54.     public PolynomialNutation getPolynomial() {
  55.         return polynomial;
  56.     }

  57.     /** Get the number of different terms in the non-polynomial part.
  58.      * @return number of different terms in the non-polynomial part
  59.      */
  60.     public int getNonPolynomialSize() {
  61.         return series.size();
  62.     }

  63.     /** Evaluate the value of the series.
  64.      * @param elements bodies elements for nutation
  65.      * @return value of the series
  66.      */
  67.     public double value(final BodiesElements elements) {

  68.         // polynomial part
  69.         final double p = polynomial.value(elements.getTC());

  70.         // non-polynomial part
  71.         double npHigh = 0;
  72.         double npLow  = 0;
  73.         for (final Map.Entry<Long, SeriesTerm> entry : series.entrySet()) {
  74.             final double v       = entry.getValue().value(elements)[0];
  75.             // Use 2Sum for high precision.
  76.             final SumAndResidual sumAndResidual = MathUtils.twoSum(npHigh, v);
  77.             npHigh = sumAndResidual.getSum();
  78.             npLow += sumAndResidual.getResidual();
  79.         }

  80.         // add the polynomial and the non-polynomial parts
  81.         return p + (npHigh + npLow);

  82.     }

  83.     /** Evaluate the value of the series.
  84.      * @param elements bodies elements for nutation
  85.      * @param <T> type fo the field elements
  86.      * @return value of the series
  87.      */
  88.     public <T extends CalculusFieldElement<T>> T value(final FieldBodiesElements<T> elements) {

  89.         // polynomial part
  90.         final T tc = elements.getTC();
  91.         final T p  = polynomial.value(tc);

  92.         // non-polynomial part
  93.         T sum = tc.getField().getZero();
  94.         for (final Map.Entry<Long, SeriesTerm> entry : series.entrySet()) {
  95.             sum = sum.add(entry.getValue().value(elements)[0]);
  96.         }

  97.         // add the polynomial and the non-polynomial parts
  98.         return p.add(sum);

  99.     }

  100.     /** This interface represents a fast evaluator for Poisson series.
  101.      * @see PoissonSeries#compile(PoissonSeries...)
  102.      * @since 6.1
  103.      */
  104.     public interface  CompiledSeries {

  105.         /** Evaluate a set of Poisson series.
  106.          * @param elements bodies elements for nutation
  107.          * @return value of the series
  108.          */
  109.         double[] value(BodiesElements elements);

  110.         /** Evaluate time derivative of a set of Poisson series.
  111.          * @param elements bodies elements for nutation
  112.          * @return time derivative of the series
  113.          */
  114.         double[] derivative(BodiesElements elements);

  115.         /** Evaluate a set of Poisson series.
  116.          * @param elements bodies elements for nutation
  117.          * @param <S> the type of the field elements
  118.          * @return value of the series
  119.          */
  120.         <S extends CalculusFieldElement<S>> S[] value(FieldBodiesElements<S> elements);

  121.         /** Evaluate time derivative of a set of Poisson series.
  122.          * @param elements bodies elements for nutation
  123.          * @param <S> the type of the field elements
  124.          * @return time derivative of the series
  125.          */
  126.         <S extends CalculusFieldElement<S>> S[] derivative(FieldBodiesElements<S> elements);

  127.     }

  128.     /** Join several nutation series, for fast simultaneous evaluation.
  129.      * @param poissonSeries Poisson series to join
  130.      * @return a single function that evaluates all series together
  131.      * @since 6.1
  132.      */
  133.     @SafeVarargs
  134.     public static CompiledSeries compile(final PoissonSeries... poissonSeries) {

  135.         // store all polynomials
  136.         final PolynomialNutation[] polynomials = new PolynomialNutation[poissonSeries.length];
  137.         for (int i = 0; i < polynomials.length; ++i) {
  138.             polynomials[i] = poissonSeries[i].polynomial;
  139.         }

  140.         // gather all series terms
  141.         final Map<Long, SeriesTerm> joinedMap = new HashMap<Long, SeriesTerm>();
  142.         for (final PoissonSeries ps : poissonSeries) {
  143.             for (Map.Entry<Long, SeriesTerm> entry : ps.series.entrySet()) {
  144.                 final long key = entry.getKey();
  145.                 if (!joinedMap.containsKey(key)) {

  146.                     // retrieve all Delaunay and planetary multipliers from the key
  147.                     final int[] m = NutationCodec.decode(key);

  148.                     // prepare a new term, ready to handle the required dimension
  149.                     final SeriesTerm term =
  150.                             SeriesTerm.buildTerm(m[0],
  151.                                                  m[1], m[2], m[3], m[4], m[5],
  152.                                                  m[6], m[7], m[8], m[9], m[10], m[11], m[12], m[13], m[14]);
  153.                     term.add(poissonSeries.length - 1, -1, Double.NaN, Double.NaN);

  154.                     // store it
  155.                     joinedMap.put(key, term);

  156.                 }
  157.             }
  158.         }

  159.         // join series by sharing terms, in order to speed up evaluation
  160.         // which is dominated by the computation of sine/cosine in each term
  161.         for (int i = 0; i < poissonSeries.length; ++i) {
  162.             for (final Map.Entry<Long, SeriesTerm> entry : poissonSeries[i].series.entrySet()) {
  163.                 final SeriesTerm singleTerm = entry.getValue();
  164.                 final SeriesTerm joinedTerm = joinedMap.get(entry.getKey());
  165.                 for (int degree = 0; degree <= singleTerm.getDegree(0); ++degree) {
  166.                     joinedTerm.add(i, degree,
  167.                                    singleTerm.getSinCoeff(0, degree),
  168.                                    singleTerm.getCosCoeff(0, degree));
  169.                 }
  170.             }
  171.         }

  172.         // use a single array for faster access
  173.         final SeriesTerm[] joinedTerms = new SeriesTerm[joinedMap.size()];
  174.         int index = 0;
  175.         for (final Map.Entry<Long, SeriesTerm> entry : joinedMap.entrySet()) {
  176.             joinedTerms[index++] = entry.getValue();
  177.         }

  178.         return new CompiledSeries() {

  179.             /** {@inheritDoc} */
  180.             @Override
  181.             public double[] value(final BodiesElements elements) {

  182.                 // non-polynomial part
  183.                 final double[] npHigh = new double[polynomials.length];
  184.                 final double[] npLow  = new double[polynomials.length];
  185.                 for (final SeriesTerm term : joinedTerms) {
  186.                     final double[] termValue = term.value(elements);
  187.                     for (int i = 0; i < termValue.length; ++i) {
  188.                         // Use 2Sum for high precision.
  189.                         final SumAndResidual sumAndResidual = MathUtils.twoSum(npHigh[i], termValue[i]);
  190.                         npHigh[i] = sumAndResidual.getSum();
  191.                         npLow[i] += sumAndResidual.getResidual();
  192.                     }
  193.                 }

  194.                 // add residual and polynomial part
  195.                 for (int i = 0; i < npHigh.length; ++i) {
  196.                     npHigh[i] += npLow[i] + polynomials[i].value(elements.getTC());
  197.                 }
  198.                 return npHigh;

  199.             }

  200.             /** {@inheritDoc} */
  201.             @Override
  202.             public double[] derivative(final BodiesElements elements) {

  203.                 // non-polynomial part
  204.                 final double[] v = new double[polynomials.length];
  205.                 for (final SeriesTerm term : joinedTerms) {
  206.                     final double[] termDerivative = term.derivative(elements);
  207.                     for (int i = 0; i < termDerivative.length; ++i) {
  208.                         v[i] += termDerivative[i];
  209.                     }
  210.                 }

  211.                 // add polynomial part
  212.                 for (int i = 0; i < v.length; ++i) {
  213.                     v[i] += polynomials[i].derivative(elements.getTC());
  214.                 }
  215.                 return v;

  216.             }

  217.             /** {@inheritDoc} */
  218.             @Override
  219.             public <S extends CalculusFieldElement<S>> S[] value(final FieldBodiesElements<S> elements) {

  220.                // non-polynomial part
  221.                 final S[] v = MathArrays.buildArray(elements.getTC().getField(), polynomials.length);
  222.                 for (final SeriesTerm term : joinedTerms) {
  223.                     final S[] termValue = term.value(elements);
  224.                     for (int i = 0; i < termValue.length; ++i) {
  225.                         v[i] = v[i].add(termValue[i]);
  226.                     }
  227.                 }

  228.                 // add polynomial part
  229.                 final S tc = elements.getTC();
  230.                 for (int i = 0; i < v.length; ++i) {
  231.                     v[i] = v[i].add(polynomials[i].value(tc));
  232.                 }
  233.                 return v;

  234.             }

  235.             /** {@inheritDoc} */
  236.             @Override
  237.             public <S extends CalculusFieldElement<S>> S[] derivative(final FieldBodiesElements<S> elements) {

  238.                // non-polynomial part
  239.                 final S[] v = MathArrays.buildArray(elements.getTC().getField(), polynomials.length);
  240.                 for (final SeriesTerm term : joinedTerms) {
  241.                     final S[] termDerivative = term.derivative(elements);
  242.                     for (int i = 0; i < termDerivative.length; ++i) {
  243.                         v[i] = v[i].add(termDerivative[i]);
  244.                     }
  245.                 }

  246.                 // add polynomial part
  247.                 final S tc = elements.getTC();
  248.                 for (int i = 0; i < v.length; ++i) {
  249.                     v[i] = v[i].add(polynomials[i].derivative(tc));
  250.                 }
  251.                 return v;

  252.             }

  253.         };

  254.     }

  255. }