SeriesTerm.java

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

  19. import java.util.Arrays;

  21. import org.hipparchus.CalculusFieldElement;
  22. import org.hipparchus.util.FastMath;
  23. import org.hipparchus.util.FieldSinCos;
  24. import org.hipparchus.util.MathArrays;
  25. import org.hipparchus.util.SinCos;
  26. import org.orekit.errors.OrekitInternalError;
  27. import org.orekit.utils.Constants;

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

  35.     /** Coefficients for the sine part. */
  36.     private double[][] sinCoeff;

  38.     /** Coefficients for the cosine part. */
  39.     private double[][] cosCoeff;

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

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

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

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

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

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

  115.         // preliminary computation
  116.         final double tc  = elements.getTC();
  117.         final double a   = argument(elements);
  118.         final SinCos sc  = FastMath.sinCos(a);

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

  132.         return values;

  134.     }

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

  142.         // preliminary computation
  143.         final double tc   = elements.getTC();
  144.         final double a    = argument(elements);
  145.         final double aDot = argumentDerivative(elements);
  146.         final SinCos sc   = FastMath.sinCos(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) * sc.sin() + (cDot + s * aDot) * sc.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 CalculusFieldElement<T>> T[] value(final FieldBodiesElements<T> elements) {

  193.         // preliminary computation
  194.         final T tc  = elements.getTC();
  195.         final T a   = argument(elements);
  196.         final FieldSinCos<T> sc = FastMath.sinCos(a);

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

  210.         return values;

  212.     }

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

  221.         // preliminary computation
  222.         final T tc   = elements.getTC();
  223.         final T a    = argument(elements);
  224.         final T aDot = argumentDerivative(elements);
  225.         final FieldSinCos<T> sc = FastMath.sinCos(a);

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

  250.         return derivatives;

  252.     }

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

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

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

  315.     }

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

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

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

  340.         return extended;

  342.     }

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

  360. }
Created with JaCoCo 0.8.12.202403310830