PoissonSeries.java

  1. /* Copyright 2002-2018 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;

  19. import java.util.HashMap;
  20. import java.util.Map;

  22. import org.hipparchus.RealFieldElement;
  23. import org.hipparchus.util.MathArrays;

  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 {

  41.     /** Polynomial part. */
  42.     private final PolynomialNutation polynomial;

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

  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.     }

  56.     /** Get the polynomial part of the series.
  57.      * @return polynomial part of the series.
  58.      */
  59.     public PolynomialNutation getPolynomial() {
  60.         return polynomial;
  61.     }

  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.     }

  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) {

  76.         // polynomial part
  77.         final double p = polynomial.value(elements.getTC());

  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.         }

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

  99.     }

  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) {

  108.         // polynomial part
  109.         final T tc = elements.getTC();
  110.         final T p  = polynomial.value(tc);

  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.         }

  118.         // add the polynomial and the non-polynomial parts
  119.         return p.add(sum);

  121.     }

  123.     /** This interface represents a fast evaluator for Poisson series.
  124.      * @see PoissonSeries#compile(PoissonSeries...)
  125.      * @since 6.1
  126.      */
  127.     public interface  CompiledSeries {

  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);

  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);

  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);

  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);

  155.     }

  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) {

  165.         // store all polynomials
  166.         final PolynomialNutation[] polynomials = new PolynomialNutation[poissonSeries.length];
  167.         for (int i = 0; i < polynomials.length; ++i) {
  168.             polynomials[i] = poissonSeries[i].polynomial;
  169.         }

  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)) {

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

  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);

  188.                     // store it
  189.                     joinedMap.put(key, term);

  191.                 }
  192.             }
  193.         }

  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.         }

  209.         // use a single array for faster access
  210.         final 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.         }

  216.         return new CompiledSeries() {

  218.             /** {@inheritDoc} */
  219.             @Override
  220.             public double[] value(final BodiesElements elements) {

  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.                 }

  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;

  248.             }

  250.             /** {@inheritDoc} */
  251.             @Override
  252.             public double[] derivative(final BodiesElements elements) {

  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.                 }

  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;

  269.             }

  271.             /** {@inheritDoc} */
  272.             @Override
  273.             public <S extends RealFieldElement<S>> S[] value(final FieldBodiesElements<S> elements) {

  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.                 }

  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;

  291.             }

  293.             /** {@inheritDoc} */
  294.             @Override
  295.             public <S extends RealFieldElement<S>> S[] derivative(final FieldBodiesElements<S> elements) {

  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.                 }

  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;

  313.             }

  315.         };

  317.     }

  319. }
Created with JaCoCo 0.7.9.201702052155