PoissonSeries.java

/* Copyright 2002-2016 CS Systèmes d'Information
 * Licensed to CS Systèmes d'Information (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * CS licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.data;

import java.lang.reflect.Array;
import java.util.HashMap;
import java.util.Map;

import org.apache.commons.math3.RealFieldElement;
import org.apache.commons.math3.util.MathArrays;

/**
 * Class representing a Poisson series for nutation or ephemeris computations.
 * <p>
 * A Poisson series is composed of a time polynomial part and a non-polynomial
 * part which consist in summation series. The {@link SeriesTerm series terms}
 * are harmonic functions (combination of sines and cosines) of polynomial
 * <em>arguments</em>. The polynomial arguments are combinations of luni-solar or
 * planetary {@link BodiesElements elements}.
 * </p>
 * @param <T> the type of the field elements
 * @author Luc Maisonobe
 * @see PoissonSeriesParser
 * @see SeriesTerm
 * @see PolynomialNutation
 */
public class PoissonSeries<T extends RealFieldElement<T>> {

    /** Polynomial part. */
    private final PolynomialNutation<T> polynomial;

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

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

    /** Get the polynomial part of the series.
     * @return polynomial part of the series.
     */
    public PolynomialNutation<T> getPolynomial() {
        return polynomial;
    }

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

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

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

        // non-polynomial part
        // compute sum accurately, using Møller-Knuth TwoSum algorithm without branching
        // the following statements must NOT be simplified, they rely on floating point
        // arithmetic properties (rounding and representable numbers)
        double npHigh = 0;
        double npLow  = 0;
        for (final SeriesTerm<T> term : series.values()) {
            final double v       = term.value(elements)[0];
            final double sum     = npHigh + v;
            final double sPrime  = sum - v;
            final double tPrime  = sum - sPrime;
            final double deltaS  = npHigh  - sPrime;
            final double deltaT  = v - tPrime;
            npLow  += deltaS   + deltaT;
            npHigh  = sum;
        }

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

    }

    /** Evaluate the value of the series.
     * @param elements bodies elements for nutation
     * @return value of the series
     */
    public T value(final FieldBodiesElements<T> elements) {

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

        // non-polynomial part
        T sum = tc.getField().getZero();
        for (final SeriesTerm<T> term : series.values()) {
            sum = sum.add(term.value(elements)[0]);
        }

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

    }

    /** This interface represents a fast evaluator for Poisson series.
     * @see PoissonSeries#compile(PoissonSeries...)
     * @param <S> the type of the field elements
     * @since 6.1
     */
    public interface  CompiledSeries<S extends RealFieldElement<S>> {

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

        /** Evaluate a set of Poisson series.
         * @param elements bodies elements for nutation
         * @return value of the series
         */
        S[] value(FieldBodiesElements<S> elements);

    }

    /** Join several nutation series, for fast simultaneous evaluation.
     * @param poissonSeries Poisson series to join
     * @return a single function that evaluates all series together
     * @param <S> the type of the field elements
     * @since 6.1
     */
    public static <S extends RealFieldElement<S>> CompiledSeries<S> compile(final PoissonSeries<S> ... poissonSeries) {

        // store all polynomials
        @SuppressWarnings("unchecked")
        final PolynomialNutation<S>[] polynomials =
                (PolynomialNutation<S>[]) Array.newInstance(PolynomialNutation.class, poissonSeries.length);
        for (int i = 0; i < polynomials.length; ++i) {
            polynomials[i] = poissonSeries[i].polynomial;
        }

        // gather all series terms
        final Map<Long, SeriesTerm<S>> joinedMap = new HashMap<Long, SeriesTerm<S>>();
        for (final PoissonSeries<S> ps : poissonSeries) {
            for (long key : ps.series.keySet()) {
                if (!joinedMap.containsKey(key)) {

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

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

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

                }
            }
        }

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

        // use a single array for faster access
        @SuppressWarnings("unchecked")
        final SeriesTerm<S>[] joinedTerms =
                joinedMap.values().toArray((SeriesTerm<S>[]) Array.newInstance(SeriesTerm.class, joinedMap.size()));

        return new CompiledSeries<S>() {

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

                // non-polynomial part
                // compute sum accurately, using Møller-Knuth TwoSum algorithm without branching
                // the following statements must NOT be simplified, they rely on floating point
                // arithmetic properties (rounding and representable numbers)
                final double[] npHigh = new double[polynomials.length];
                final double[] npLow  = new double[polynomials.length];
                for (final SeriesTerm<S> term : joinedTerms) {
                    final double[] termValue = term.value(elements);
                    for (int i = 0; i < termValue.length; ++i) {
                        final double v       = termValue[i];
                        final double sum     = npHigh[i] + v;
                        final double sPrime  = sum - v;
                        final double tPrime  = sum - sPrime;
                        final double deltaS  = npHigh[i]  - sPrime;
                        final double deltaT  = v - tPrime;
                        npLow[i]  += deltaS   + deltaT;
                        npHigh[i]  = sum;
                    }
                }

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

            }

            /** {@inheritDoc} */
            @Override
            public S[] value(final FieldBodiesElements<S> elements) {

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

                // add residual and polynomial part
                final S tc = elements.getTC();
                for (int i = 0; i < v.length; ++i) {
                    v[i] = v[i].add(polynomials[i].value(tc));
                }
                return v;

            }

        };

    }

}