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.     @SuppressWarnings("unchecked")
  91.     public FieldHansenTesseralLinear(final int nMax, final int s, final int j, final int n0,
  92.                                      final int maxHansen, final Field<T> field) {
  93.         //Initialize the fields
  94.         this.offset = nMax + 1;
  95.         this.Nmin = -nMax - 1;
  96.         this.N0 = -n0 - 4;

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

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

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

  119.     }

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

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

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

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

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

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

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

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

  196.     }

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

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

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

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

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

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

  240.     }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  336. }
Created with JaCoCo 0.8.12.202403310830