FieldHansenTesseralLinear.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.propagation.semianalytical.dsst.utilities.hansen;

  19. import java.lang.reflect.Array;

  21. import org.hipparchus.CalculusFieldElement;
  22. import org.hipparchus.Field;
  23. import org.hipparchus.analysis.differentiation.FieldGradient;
  24. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  25. import org.hipparchus.util.FastMath;
  26. import org.hipparchus.util.MathArrays;
  27. import org.orekit.propagation.semianalytical.dsst.utilities.NewcombOperators;

  29. /**
  30.  * Hansen coefficients K(t,n,s) for t!=0 and n < 0.
  31.  * <p>
  32.  * Implements Collins 4-236 or Danielson 2.7.3-(9) for Hansen Coefficients and
  33.  * Collins 4-240 for derivatives. The recursions are transformed into
  34.  * composition of linear transformations to obtain the associated polynomials
  35.  * for coefficients and their derivatives - see Petre's paper
  36.  *
  37.  * @author Petre Bazavan
  38.  * @author Lucian Barbulescu
  39.  * @author Bryan Cazabonne
  40.  * @param <T> type of the field elements
  41.  */
  42. public class FieldHansenTesseralLinear <T extends CalculusFieldElement<T>> {

  44.     /** The number of coefficients that will be computed with a set of roots. */
  45.     private static final int SLICE = 10;

  47.     /**
  48.      * The first vector of polynomials associated to Hansen coefficients and
  49.      * derivatives.
  50.      */
  51.     private PolynomialFunction[][] mpvec;

  53.     /** The second vector of polynomials associated only to derivatives. */
  54.     private PolynomialFunction[][] mpvecDeriv;

  56.     /** The Hansen coefficients used as roots. */
  57.     private final T[][] hansenRoot;

  59.     /** The derivatives of the Hansen coefficients used as roots. */
  60.     private final T[][] hansenDerivRoot;

  62.     /** The minimum value for the order. */
  63.     private final int Nmin;

  65.     /** The index of the initial condition, Petre's paper. */
  66.     private final int N0;

  68.     /** The number of slices needed to compute the coefficients. */
  69.     private final int numSlices;

  71.     /**
  72.      * The offset used to identify the polynomial that corresponds to a negative.
  73.      * value of n in the internal array that starts at 0
  74.      */
  75.     private final int offset;

  77.     /** The objects used to calculate initial data by means of Newcomb operators. */
  78.     private final FieldHansenCoefficientsBySeries<T>[] hansenInit;

  80.     /**
  81.      * Constructor.
  82.      *
  83.      * @param nMax the maximum (absolute) value of n parameter
  84.      * @param s s parameter
  85.      * @param j j parameter
  86.      * @param n0 the minimum (absolute) value of n
  87.      * @param maxHansen maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
  88.      * @param field field used by default
  89.      */
  90.     public FieldHansenTesseralLinear(final int nMax, final int s, final int j, final int n0,
  91.                                      final int maxHansen, final Field<T> field) {
  92.         //Initialize the fields
  93.         this.offset = nMax + 1;
  94.         this.Nmin = -nMax - 1;
  95.         this.N0 = -n0 - 4;

  97.         final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
  98.         //Ensure that only the needed terms are computed
  99.         this.hansenInit = (FieldHansenCoefficientsBySeries<T>[]) Array.newInstance(FieldHansenCoefficientsBySeries.class, maxRoots);
  100.         for (int i = 0; i < maxRoots; i++) {
  101.             this.hansenInit[i] = new FieldHansenCoefficientsBySeries<>(N0 - i + 3, s, j, maxHansen);
  102.         }

  104.         // The first 4 values are computed with series. No linear combination is needed.
  105.         final int size = N0 - Nmin;
  106.         this.numSlices = (int) FastMath.max(FastMath.ceil(((double) size) / SLICE), 1);
  107.         hansenRoot = MathArrays.buildArray(field, numSlices, 4);
  108.         hansenDerivRoot = MathArrays.buildArray(field, numSlices, 4);
  109.         if (size > 0) {
  110.             mpvec = new PolynomialFunction[size][];
  111.             mpvecDeriv = new PolynomialFunction[size][];

  113.             // Prepare the database of the associated polynomials
  114.             HansenUtilities.generateTesseralPolynomials(N0, Nmin, offset, SLICE, j, s,
  115.                                                         mpvec, mpvecDeriv);
  116.         }

  118.     }

  120.     /**
  121.      * Compute the values for the first four coefficients and their derivatives by means of series.
  122.      *
  123.      * @param e2 e²
  124.      * @param chi &Chi;
  125.      * @param chi2 &Chi;²
  126.      */
  127.     public void computeInitValues(final T e2, final T chi, final T chi2) {
  128.         // compute the values for n, n+1, n+2 and n+3 by series
  129.         // See Danielson 2.7.3-(10)
  130.         //Ensure that only the needed terms are computed
  131.         final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
  132.         for (int i = 0; i < maxRoots; i++) {
  133.             final FieldGradient<T> hansenKernel = hansenInit[i].getValueGradient(e2, chi, chi2);
  134.             this.hansenRoot[0][i] = hansenKernel.getValue();
  135.             this.hansenDerivRoot[0][i] = hansenKernel.getPartialDerivative(0);
  136.         }

  138.         for (int i = 1; i < numSlices; i++) {
  139.             for (int k = 0; k < 4; k++) {
  140.                 final PolynomialFunction[] mv = mpvec[N0 - (i * SLICE) - k + 3 + offset];
  141.                 final PolynomialFunction[] sv = mpvecDeriv[N0 - (i * SLICE) - k + 3 + offset];

  143.                 hansenDerivRoot[i][k] = mv[3].value(chi).multiply(hansenDerivRoot[i - 1][3]).
  144.                                         add(mv[2].value(chi).multiply(hansenDerivRoot[i - 1][2])).
  145.                                         add(mv[1].value(chi).multiply(hansenDerivRoot[i - 1][1])).
  146.                                         add(mv[0].value(chi).multiply(hansenDerivRoot[i - 1][0])).
  147.                                         add(sv[3].value(chi).multiply(hansenRoot[i - 1][3])).
  148.                                         add(sv[2].value(chi).multiply(hansenRoot[i - 1][2])).
  149.                                         add(sv[1].value(chi).multiply(hansenRoot[i - 1][1])).
  150.                                         add(sv[0].value(chi).multiply(hansenRoot[i - 1][0]));

  152.                 hansenRoot[i][k] =  mv[3].value(chi).multiply(hansenRoot[i - 1][3]).
  153.                                     add(mv[2].value(chi).multiply(hansenRoot[i - 1][2])).
  154.                                     add(mv[1].value(chi).multiply(hansenRoot[i - 1][1])).
  155.                                     add(mv[0].value(chi).multiply(hansenRoot[i - 1][0]));
  156.             }
  157.         }
  158.     }

  160.     /**
  161.      * Compute the value of the Hansen coefficient K<sub>j</sub><sup>-n-1, s</sup>.
  162.      *
  163.      * @param mnm1 -n-1
  164.      * @param chi χ
  165.      * @return the coefficient K<sub>j</sub><sup>-n-1, s</sup>
  166.      */
  167.     public T getValue(final int mnm1, final T chi) {
  168.         //Compute n
  169.         final int n = -mnm1 - 1;

  171.         //Compute the potential slice
  172.         int sliceNo = (n + N0 + 4) / SLICE;
  173.         if (sliceNo < numSlices) {
  174.             //Compute the index within the slice
  175.             final int indexInSlice = (n + N0 + 4) % SLICE;

  177.             //Check if a root must be returned
  178.             if (indexInSlice <= 3) {
  179.                 return hansenRoot[sliceNo][indexInSlice];
  180.             }
  181.         } else {
  182.             //the value was a potential root for a slice, but that slice was not required
  183.             //Decrease the slice number
  184.             sliceNo--;
  185.         }

  187.         // Computes the coefficient by linear transformation
  188.         // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
  189.         final PolynomialFunction[] v = mpvec[mnm1 + offset];
  190.         return v[3].value(chi).multiply(hansenRoot[sliceNo][3]).
  191.                add(v[2].value(chi).multiply(hansenRoot[sliceNo][2])).
  192.                add(v[1].value(chi).multiply(hansenRoot[sliceNo][1])).
  193.                add(v[0].value(chi).multiply(hansenRoot[sliceNo][0]));

  195.     }

  197.     /**
  198.      * Compute the value of the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de².
  199.      *
  200.      * @param mnm1 -n-1
  201.      * @param chi χ
  202.      * @return the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de²
  203.      */
  204.     public T getDerivative(final int mnm1, final T chi) {

  206.         //Compute n
  207.         final int n = -mnm1 - 1;

  209.         //Compute the potential slice
  210.         int sliceNo = (n + N0 + 4) / SLICE;
  211.         if (sliceNo < numSlices) {
  212.             //Compute the index within the slice
  213.             final int indexInSlice = (n + N0 + 4) % SLICE;

  215.             //Check if a root must be returned
  216.             if (indexInSlice <= 3) {
  217.                 return hansenDerivRoot[sliceNo][indexInSlice];
  218.             }
  219.         } else {
  220.             //the value was a potential root for a slice, but that slice was not required
  221.             //Decrease the slice number
  222.             sliceNo--;
  223.         }

  225.         // Computes the coefficient by linear transformation
  226.         // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
  227.         final PolynomialFunction[] v = mpvec[mnm1 + this.offset];
  228.         final PolynomialFunction[] vv = mpvecDeriv[mnm1 + this.offset];

  230.         return v[3].value(chi).multiply(hansenDerivRoot[sliceNo][3]).
  231.                add(v[2].value(chi).multiply(hansenDerivRoot[sliceNo][2])).
  232.                add(v[1].value(chi).multiply(hansenDerivRoot[sliceNo][1])).
  233.                add(v[0].value(chi).multiply(hansenDerivRoot[sliceNo][0])).
  234.                add(vv[3].value(chi).multiply(hansenRoot[sliceNo][3])).
  235.                add(vv[2].value(chi).multiply(hansenRoot[sliceNo][2])).
  236.                add( vv[1].value(chi).multiply(hansenRoot[sliceNo][1])).
  237.                add(vv[0].value(chi).multiply(hansenRoot[sliceNo][0]));

  239.     }

  241.     /**
  242.      * Compute a Hansen coefficient with series.
  243.      * <p>
  244.      * This class implements the computation of the Hansen kernels
  245.      * through a power series in e² and that is using
  246.      * modified Newcomb operators. The reference formulae can be found
  247.      * in Danielson 2.7.3-10 and 3.3-5
  248.      * </p>
  249.      */
  250.     private static class FieldHansenCoefficientsBySeries <T extends CalculusFieldElement<T>> {

  252.         /** -n-1 coefficient. */
  253.         private final int mnm1;

  255.         /** s coefficient. */
  256.         private final int s;

  258.         /** j coefficient. */
  259.         private final int j;

  261.         /** Max power in e² for the Newcomb's series expansion. */
  262.         private final int maxNewcomb;

  264.         /** Polynomial representing the serie. */
  265.         private final PolynomialFunction polynomial;

  267.         /**
  268.          * Class constructor.
  269.          *
  270.          * @param mnm1      -n-1 value
  271.          * @param s         s value
  272.          * @param j         j value
  273.          * @param maxHansen max power of e² in series expansion
  274.          */
  275.         FieldHansenCoefficientsBySeries(final int mnm1, final int s,
  276.                                         final int j, final int maxHansen) {
  277.             this.mnm1 = mnm1;
  278.             this.s = s;
  279.             this.j = j;
  280.             this.maxNewcomb = maxHansen;
  281.             this.polynomial = generatePolynomial();
  282.         }

  284.         /** Computes the value of Hansen kernel and its derivative at e².
  285.          * <p>
  286.          * The formulae applied are described in Danielson 2.7.3-10 and
  287.          * 3.3-5
  288.          * </p>
  289.          * @param e2 e²
  290.          * @param chi &Chi;
  291.          * @param chi2 &Chi;²
  292.          * @return the value of the Hansen coefficient and its derivative for e²
  293.          */
  294.         private FieldGradient<T> getValueGradient(final T e2, final T chi, final T chi2) {

  296.             final T zero = e2.getField().getZero();
  297.             //Estimation of the serie expansion at e2
  298.             final FieldGradient<T> serie = polynomial.value(FieldGradient.variable(1, 0, e2));

  300.             final T value      =  FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getValue()).divide(chi);
  301.             final T coef       = zero.subtract(mnm1 + 1.5);
  302.             final T derivative = coef.multiply(chi2).multiply(value).
  303.                             add(FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getPartialDerivative(0)).divide(chi));
  304.             return new FieldGradient<>(value, derivative);
  305.         }

  307.         /** Generate the serie expansion in e².
  308.          * <p>
  309.          * Generate the series expansion in e² used in the formulation
  310.          * of the Hansen kernel (see Danielson 2.7.3-10):
  311.          * &Sigma; Y<sup>ns</sup><sub>α+a,α+b</sub>
  312.          * *e<sup>2α</sup>
  313.          * </p>
  314.          * @return polynomial representing the power serie expansion
  315.          */
  316.         private PolynomialFunction generatePolynomial() {
  317.             // Initialization
  318.             final int aHT = FastMath.max(j - s, 0);
  319.             final int bHT = FastMath.max(s - j, 0);

  321.             final double[] coefficients = new double[maxNewcomb + 1];

  323.             //Loop for getting the Newcomb operators
  324.             for (int alphaHT = 0; alphaHT <= maxNewcomb; alphaHT++) {
  325.                 coefficients[alphaHT] =
  326.                         NewcombOperators.getValue(alphaHT + aHT, alphaHT + bHT, mnm1, s);
  327.             }

  329.             //Creation of the polynomial
  330.             return new PolynomialFunction(coefficients);
  331.         }
  332.     }

  334. }
Created with JaCoCo 0.8.12.202403310830