SeriesTerm.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.Arrays;

  21. import org.hipparchus.RealFieldElement;
  22. import org.hipparchus.util.FastMath;
  23. import org.hipparchus.util.MathArrays;
  24. import org.orekit.errors.OrekitInternalError;
  25. import org.orekit.utils.Constants;

  27. /** Base class for nutation series terms.
  28.  * @author Luc Maisonobe
  29.  * @see PoissonSeries
  30.  */
  31. abstract class SeriesTerm {

  33.     /** Coefficients for the sine part. */
  34.     private double[][] sinCoeff;

  36.     /** Coefficients for the cosine part. */
  37.     private double[][] cosCoeff;

  39.     /** Simple constructor for the base class.
  40.      */
  41.     protected SeriesTerm() {
  42.         this.sinCoeff = new double[0][0];
  43.         this.cosCoeff = new double[0][0];
  44.     }

  46.     /** Get the degree of the function component.
  47.      * @param index index of the function component (must be less than dimension)
  48.      * @return degree of the function component
  49.      */
  50.     public int getDegree(final int index) {
  51.         return  sinCoeff[index].length - 1;
  52.     }

  54.     /** Add a pair of values to existing term coefficients.
  55.      * <p>
  56.      * Despite it would seem logical to simply set coefficients
  57.      * rather than add to them, this does not work for some IERS
  58.      * files. As an example in table 5.3a in IERS conventions 2003,
  59.      * the coefficients for luni-solar term for 2F+Ω with period
  60.      * 13.633 days appears twice with different coefficients, as
  61.      * well as term for 2(F+D+Ω)+l with period 5.643 days, term for
  62.      * 2(F+D+Ω)-l with period 9.557 days, term for 2(Ω-l') with
  63.      * period -173.318 days, term for 2D-l with period 31.812 days ...
  64.      * 35 different duplicated terms have been identified in the
  65.      * tables 5.3a and 5.3b in IERS conventions 2003.
  66.      * The coefficients read in lines duplicating a term must be
  67.      * added together.
  68.      * </p>
  69.      * @param index index of the components (will automatically
  70.      * increase dimension if needed)
  71.      * @param degree degree of the coefficients, may be negative if
  72.      * the term does not contribute to component (will automatically
  73.      * increase {@link #getDegree() degree} of the component if needed)
  74.      * @param sinID coefficient for the sine part, at index and degree
  75.      * @param cosID coefficient for the cosine part, at index and degree
  76.      */
  77.     public void add(final int index, final int degree,
  78.                     final double sinID, final double cosID) {
  79.         sinCoeff = extendArray(index, degree, sinCoeff);
  80.         cosCoeff = extendArray(index, degree, cosCoeff);
  81.         if (degree >= 0) {
  82.             sinCoeff[index][degree] += sinID;
  83.             cosCoeff[index][degree] += cosID;
  84.         }
  85.     }

  87.     /** Get a coefficient for the sine part.
  88.      * @param index index of the function component (must be less than dimension)
  89.      * @param degree degree of the coefficients
  90.      * (must be less than {@link #getDegree() degree} for the component)
  91.      * @return coefficient for the sine part, at index and degree
  92.      */
  93.     public double getSinCoeff(final int index, final int degree) {
  94.         return sinCoeff[index][degree];
  95.     }

  97.     /** Get a coefficient for the cosine part.
  98.      * @param index index of the function component (must be less than dimension)
  99.      * @param degree degree of the coefficients
  100.      * (must be less than {@link #getDegree() degree} for the component)
  101.      * @return coefficient for the cosine part, at index and degree
  102.      */
  103.     public double getCosCoeff(final int index, final int degree) {
  104.         return cosCoeff[index][degree];
  105.     }

  107.     /** Evaluate the value of the series term.
  108.      * @param elements bodies elements for nutation
  109.      * @return value of the series term
  110.      */
  111.     public double[] value(final BodiesElements elements) {

  113.         // preliminary computation
  114.         final double tc  = elements.getTC();
  115.         final double a   = argument(elements);
  116.         final double sin = FastMath.sin(a);
  117.         final double cos = FastMath.cos(a);

  119.         // compute each function
  120.         final double[] values = new double[sinCoeff.length];
  121.         for (int i = 0; i < values.length; ++i) {
  122.             double s = 0;
  123.             double c = 0;
  124.             for (int j = sinCoeff[i].length - 1; j >= 0; --j) {
  125.                 s = s * tc + sinCoeff[i][j];
  126.                 c = c * tc + cosCoeff[i][j];
  127.             }
  128.             values[i] = s * sin + c * cos;
  129.         }

  131.         return values;

  133.     }

  135.     /** Evaluate the time derivative of the series term.
  136.      * @param elements bodies elements for nutation
  137.      * @return time derivative of the series term
  138.      */
  139.     public double[] derivative(final BodiesElements elements) {

  141.         // preliminary computation
  142.         final double tc   = elements.getTC();
  143.         final double a    = argument(elements);
  144.         final double aDot = argumentDerivative(elements);
  145.         final double sin  = FastMath.sin(a);
  146.         final double cos  = FastMath.cos(a);

  148.         // compute each function
  149.         final double[] derivatives = new double[sinCoeff.length];
  150.         for (int i = 0; i < derivatives.length; ++i) {
  151.             double s    = 0;
  152.             double c    = 0;
  153.             double sDot = 0;
  154.             double cDot = 0;
  155.             if (sinCoeff[i].length > 0) {
  156.                 for (int j = sinCoeff[i].length - 1; j > 0; --j) {
  157.                     s    = s    * tc +     sinCoeff[i][j];
  158.                     c    = c    * tc +     cosCoeff[i][j];
  159.                     sDot = sDot * tc + j * sinCoeff[i][j];
  160.                     cDot = cDot * tc + j * cosCoeff[i][j];
  161.                 }
  162.                 s     = s * tc + sinCoeff[i][0];
  163.                 c     = c * tc + cosCoeff[i][0];
  164.                 sDot /= Constants.JULIAN_CENTURY;
  165.                 cDot /= Constants.JULIAN_CENTURY;
  166.             }
  167.             derivatives[i] = (sDot - c * aDot) * sin + (cDot + s * aDot) * cos;
  168.         }

  170.         return derivatives;

  172.     }

  174.     /** Compute the argument for the current date.
  175.      * @param elements luni-solar and planetary elements for the current date
  176.      * @return current value of the argument
  177.      */
  178.     protected abstract double argument(BodiesElements elements);

  180.     /** Compute the time derivative of the argument for the current date.
  181.      * @param elements luni-solar and planetary elements for the current date
  182.      * @return current time derivative of the argument
  183.      */
  184.     protected abstract double argumentDerivative(BodiesElements elements);

  186.     /** Evaluate the value of the series term.
  187.      * @param elements bodies elements for nutation
  188.      * @param <T> the type of the field elements
  189.      * @return value of the series term
  190.      */
  191.     public <T extends RealFieldElement<T>> T[] value(final FieldBodiesElements<T> elements) {

  193.         // preliminary computation
  194.         final T tc  = elements.getTC();
  195.         final T a   = argument(elements);
  196.         final T sin = a.sin();
  197.         final T cos = a.cos();

  199.         // compute each function
  200.         final T[] values = MathArrays.buildArray(tc.getField(), sinCoeff.length);
  201.         for (int i = 0; i < values.length; ++i) {
  202.             T s = tc.getField().getZero();
  203.             T c = tc.getField().getZero();
  204.             for (int j = sinCoeff[i].length - 1; j >= 0; --j) {
  205.                 s = s.multiply(tc).add(sinCoeff[i][j]);
  206.                 c = c.multiply(tc).add(cosCoeff[i][j]);
  207.             }
  208.             values[i] = s.multiply(sin).add(c.multiply(cos));
  209.         }

  211.         return values;

  213.     }

  215.     /** Evaluate the time derivative of the series term.
  216.      * @param elements bodies elements for nutation
  217.      * @param <T> the type of the field elements
  218.      * @return time derivative of the series term
  219.      */
  220.     public <T extends RealFieldElement<T>> T[] derivative(final FieldBodiesElements<T> elements) {

  222.         // preliminary computation
  223.         final T tc   = elements.getTC();
  224.         final T a    = argument(elements);
  225.         final T aDot = argumentDerivative(elements);
  226.         final T sin  = a.sin();
  227.         final T cos  = a.cos();

  229.         // compute each function
  230.         final T[] derivatives = MathArrays.buildArray(tc.getField(), sinCoeff.length);
  231.         for (int i = 0; i < derivatives.length; ++i) {
  232.             T s    = tc.getField().getZero();
  233.             T c    = tc.getField().getZero();
  234.             T sDot = tc.getField().getZero();
  235.             T cDot = tc.getField().getZero();
  236.             if (sinCoeff[i].length > 0) {
  237.                 for (int j = sinCoeff[i].length - 1; j > 0; --j) {
  238.                     s    = s.multiply(tc).add(sinCoeff[i][j]);
  239.                     c    = c.multiply(tc).add(cosCoeff[i][j]);
  240.                     sDot = sDot.multiply(tc).add(j * sinCoeff[i][j]);
  241.                     cDot = cDot.multiply(tc).add(j * cosCoeff[i][j]);
  242.                 }
  243.                 s    = s.multiply(tc).add(sinCoeff[i][0]);
  244.                 c    = c.multiply(tc).add(cosCoeff[i][0]);
  245.                 sDot = sDot.divide(Constants.JULIAN_CENTURY);
  246.                 cDot = cDot.divide(Constants.JULIAN_CENTURY);
  247.             }
  248.             derivatives[i] = sDot.subtract(c.multiply(aDot)).multiply(sin).
  249.                              add(cDot.add(s.multiply(aDot)).multiply(cos));
  250.         }

  252.         return derivatives;

  254.     }

  256.     /** Compute the argument for the current date.
  257.      * @param elements luni-solar and planetary elements for the current date
  258.      * @param <T> the type of the field elements
  259.      * @return current value of the argument
  260.      */
  261.     protected abstract <T extends RealFieldElement<T>> T argument(FieldBodiesElements<T> elements);

  263.     /** Compute the time derivative of the argument for the current date.
  264.      * @param elements luni-solar and planetary elements for the current date
  265.      * @param <T> the type of the field elements
  266.      * @return current time derivative of the argument
  267.      */
  268.     protected abstract <T extends RealFieldElement<T>> T argumentDerivative(FieldBodiesElements<T> elements);

  270.     /** Factory method for building the appropriate object.
  271.      * <p>The method checks the null coefficients and build an instance
  272.      * of an appropriate type to avoid too many unnecessary multiplications
  273.      * by zero coefficients.</p>
  274.      * @param cGamma coefficient for γ = GMST + π tide parameter
  275.      * @param cL coefficient for mean anomaly of the Moon
  276.      * @param cLPrime coefficient for mean anomaly of the Sun
  277.      * @param cF coefficient for L - Ω where L is the mean longitude of the Moon
  278.      * @param cD coefficient for mean elongation of the Moon from the Sun
  279.      * @param cOmega coefficient for mean longitude of the ascending node of the Moon
  280.      * @param cMe coefficient for mean Mercury longitude
  281.      * @param cVe coefficient for mean Venus longitude
  282.      * @param cE coefficient for mean Earth longitude
  283.      * @param cMa coefficient for mean Mars longitude
  284.      * @param cJu coefficient for mean Jupiter longitude
  285.      * @param cSa coefficient for mean Saturn longitude
  286.      * @param cUr coefficient for mean Uranus longitude
  287.      * @param cNe coefficient for mean Neptune longitude
  288.      * @param cPa coefficient for general accumulated precession in longitude
  289.      * @return a nutation serie term instance well suited for the set of coefficients
  290.      */
  291.     public static  SeriesTerm buildTerm(final int cGamma,
  292.                                         final int cL, final int cLPrime, final int cF,
  293.                                         final int cD, final int cOmega,
  294.                                         final int cMe, final int cVe, final int cE,
  295.                                         final int cMa, final int cJu, final int cSa,
  296.                                         final int cUr, final int cNe, final int cPa) {
  297.         if (cGamma == 0 && cL == 0 && cLPrime == 0 && cF == 0 && cD == 0 && cOmega == 0) {
  298.             return new PlanetaryTerm(cMe, cVe, cE, cMa, cJu, cSa, cUr, cNe, cPa);
  299.         } else if (cGamma == 0 &&
  300.                    cMe == 0 && cVe == 0 && cE == 0 && cMa == 0 && cJu == 0 &&
  301.                    cSa == 0 && cUr == 0 && cNe == 0 && cPa == 0) {
  302.             return new LuniSolarTerm(cL, cLPrime, cF, cD, cOmega);
  303.         } else if (cGamma != 0 &&
  304.                    cMe == 0 && cVe == 0 && cE == 0 && cMa == 0 && cJu == 0 &&
  305.                    cSa == 0 && cUr == 0 && cNe == 0 && cPa == 0) {
  306.             return new TideTerm(cGamma, cL, cLPrime, cF, cD, cOmega);
  307.         } else if (cGamma == 0 && cLPrime == 0 && cUr == 0 && cNe == 0 && cPa == 0) {
  308.             return new NoFarPlanetsTerm(cL, cF, cD, cOmega,
  309.                                         cMe, cVe, cE, cMa, cJu, cSa);
  310.         } else if (cGamma == 0) {
  311.             return new GeneralTerm(cL, cLPrime, cF, cD, cOmega,
  312.                                    cMe, cVe, cE, cMa, cJu, cSa, cUr, cNe, cPa);
  313.         } else {
  314.             throw new OrekitInternalError(null);
  315.         }

  317.     }

  319.     /** Extend an array to hold at least index and degree.
  320.      * @param index index of the function
  321.      * @param degree degree of the coefficients
  322.      * @param array to extend
  323.      * @return extended array
  324.      */
  325.     private static double[][] extendArray(final int index, final int degree,
  326.                                           final double[][] array) {

  328.         // extend the number of rows if needed
  329.         final double[][] extended;
  330.         if (array.length > index) {
  331.             extended = array;
  332.         } else {
  333.             final int rows =  index + 1;
  334.             extended = new double[rows][];
  335.             System.arraycopy(array, 0, extended, 0, array.length);
  336.             Arrays.fill(extended, array.length, index + 1, new double[0]);
  337.         }

  339.         // extend the number of elements in the row if needed
  340.         extended[index] = extendArray(degree, extended[index]);

  342.         return extended;

  344.     }

  346.     /** Extend an array to hold at least degree.
  347.      * @param degree degree of the coefficients
  348.      * @param array to extend
  349.      * @return extended array
  350.      */
  351.     private static double[] extendArray(final int degree, final double[] array) {
  352.         // extend the number of elements if needed
  353.         if (array.length > degree) {
  354.             return array;
  355.         } else {
  356.             final double[] extended = new double[degree + 1];
  357.             System.arraycopy(array, 0, extended, 0, array.length);
  358.             return extended;
  359.         }
  360.     }

  362. }
Created with JaCoCo 0.7.9.201702052155