FieldHansenTesseralLinear.java

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

  18. import java.lang.reflect.Array;

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

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

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

  42.     /**
  43.      * The first vector of polynomials associated to Hansen coefficients and
  44.      * derivatives.
  45.      */
  46.     private PolynomialFunction[][] mpvec;

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

  49.     /** The Hansen coefficients used as roots. */
  50.     private T[][] hansenRoot;

  51.     /** The derivatives of the Hansen coefficients used as roots. */
  52.     private T[][] hansenDerivRoot;

  53.     /** The minimum value for the order. */
  54.     private int Nmin;

  55.     /** The index of the initial condition, Petre's paper. */
  56.     private int N0;

  57.     /** The s coefficient. */
  58.     private int s;

  59.     /** The j coefficient. */
  60.     private int j;

  61.     /** The number of slices needed to compute the coefficients. */
  62.     private int numSlices;

  63.     /**
  64.      * The offset used to identify the polynomial that corresponds to a negative.
  65.      * value of n in the internal array that starts at 0
  66.      */
  67.     private int offset;

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

  70.     /**
  71.      * Constructor.
  72.      *
  73.      * @param nMax the maximum (absolute) value of n parameter
  74.      * @param s s parameter
  75.      * @param j j parameter
  76.      * @param n0 the minimum (absolute) value of n
  77.      * @param maxHansen maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
  78.      * @param field field used by default
  79.      */
  80.     @SuppressWarnings("unchecked")
  81.     public FieldHansenTesseralLinear(final int nMax, final int s, final int j, final int n0,
  82.                                      final int maxHansen, final Field<T> field) {
  83.         //Initialize the fields
  84.         this.offset = nMax + 1;
  85.         this.Nmin = -nMax - 1;
  86.         this.N0 = -n0 - 4;
  87.         this.s = s;
  88.         this.j = j;

  89.         final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
  90.         //Ensure that only the needed terms are computed
  91.         this.hansenInit = (FieldHansenCoefficientsBySeries<T>[]) Array.newInstance(FieldHansenCoefficientsBySeries.class, maxRoots);
  92.         for (int i = 0; i < maxRoots; i++) {
  93.             this.hansenInit[i] = new FieldHansenCoefficientsBySeries<>(N0 - i + 3, s, j, maxHansen, field);
  94.         }

  95.         // The first 4 values are computed with series. No linear combination is needed.
  96.         final int size = N0 - Nmin;
  97.         this.numSlices = (int) FastMath.max(FastMath.ceil(((double) size) / SLICE), 1);
  98.         hansenRoot = MathArrays.buildArray(field, numSlices, 4);
  99.         hansenDerivRoot = MathArrays.buildArray(field, numSlices, 4);
  100.         if (size > 0) {
  101.             mpvec = new PolynomialFunction[size][];
  102.             mpvecDeriv = new PolynomialFunction[size][];

  103.             // Prepare the database of the associated polynomials
  104.             generatePolynomials();
  105.         }

  106.     }

  107.     /**
  108.      * Compute polynomial coefficient a.
  109.      *
  110.      *  <p>
  111.      *  It is used to generate the coefficient for K<sub>j</sub><sup>-n, s</sup> when computing K<sub>j</sub><sup>-n-1, s</sup>
  112.      *  and the coefficient for dK<sub>j</sub><sup>-n, s</sup> / de² when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
  113.      *  </p>
  114.      *
  115.      *  <p>
  116.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  117.      *  </p>
  118.      *
  119.      * @param mnm1 -n-1
  120.      * @return the polynomial
  121.      */
  122.     private PolynomialFunction a(final int mnm1) {
  123.         // Collins 4-236, Danielson 2.7.3-(9)
  124.         final double r1 = (mnm1 + 2.) * (2. * mnm1 + 5.);
  125.         final double r2 = (2. + mnm1 + s) * (2. + mnm1 - s);
  126.         return new PolynomialFunction(new double[] {
  127.             0.0, 0.0, r1 / r2
  128.         });
  129.     }

  130.     /**
  131.      * Compute polynomial coefficient b.
  132.      *
  133.      *  <p>
  134.      *  It is used to generate the coefficient for K<sub>j</sub><sup>-n+1, s</sup> when computing K<sub>j</sub><sup>-n-1, s</sup>
  135.      *  and the coefficient for dK<sub>j</sub><sup>-n+1, s</sup> / de² when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
  136.      *  </p>
  137.      *
  138.      *  <p>
  139.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  140.      *  </p>
  141.      *
  142.      * @param mnm1 -n-1
  143.      * @return the polynomial
  144.      */
  145.     private PolynomialFunction b(final int mnm1) {
  146.         // Collins 4-236, Danielson 2.7.3-(9)
  147.         final double r2 = (2. + mnm1 + s) * (2. + mnm1 - s);
  148.         final double d1 = (mnm1 + 3.) * 2. * j * s / (r2 * (mnm1 + 4.));
  149.         final double d2 = (mnm1 + 3.) * (mnm1 + 2.) / r2;
  150.         return new PolynomialFunction(new double[] {
  151.             0.0, -d1, -d2
  152.         });
  153.     }

  154.     /**
  155.      * Compute polynomial coefficient c.
  156.      *
  157.      *  <p>
  158.      *  It is used to generate the coefficient for K<sub>j</sub><sup>-n+3, s</sup> when computing K<sub>j</sub><sup>-n-1, s</sup>
  159.      *  and the coefficient for dK<sub>j</sub><sup>-n+3, s</sup> / de² when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
  160.      *  </p>
  161.      *
  162.      *  <p>
  163.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  164.      *  </p>
  165.      *
  166.      * @param mnm1 -n-1
  167.      * @return the polynomial
  168.      */
  169.     private PolynomialFunction c(final int mnm1) {
  170.         // Collins 4-236, Danielson 2.7.3-(9)
  171.         final double r1 = j * j * (mnm1 + 2.);
  172.         final double r2 = (mnm1 + 4.) * (2. + mnm1 + s) * (2. + mnm1 - s);

  173.         return new PolynomialFunction(new double[] {
  174.             0.0, 0.0, r1 / r2
  175.         });
  176.     }

  177.     /**
  178.      * Compute polynomial coefficient d.
  179.      *
  180.      *  <p>
  181.      *  It is used to generate the coefficient for K<sub>j</sub><sup>-n-1, s</sup> / dχ when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
  182.      *  </p>
  183.      *
  184.      *  <p>
  185.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  186.      *  </p>
  187.      *
  188.      * @param mnm1 -n-1
  189.      * @return the polynomial
  190.      */
  191.     private PolynomialFunction d(final int mnm1) {
  192.         // Collins 4-236, Danielson 2.7.3-(9)
  193.         return new PolynomialFunction(new double[] {
  194.             0.0, 0.0, 1.0
  195.         });
  196.     }

  197.     /**
  198.      * Compute polynomial coefficient f.
  199.      *
  200.      *  <p>
  201.      *  It is used to generate the coefficient for K<sub>j</sub><sup>-n+1, s</sup> / dχ when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
  202.      *  </p>
  203.      *
  204.      *  <p>
  205.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  206.      *  </p>
  207.      *
  208.      * @param n index
  209.      * @return the polynomial
  210.      */
  211.     private PolynomialFunction f(final int n) {
  212.         // Collins 4-236, Danielson 2.7.3-(9)
  213.         final double r1 = (n + 3.0) * j * s;
  214.         final double r2 = (n + 4.0) * (2.0 + n + s) * (2.0 + n - s);
  215.         return new PolynomialFunction(new double[] {
  216.             0.0, 0.0, 0.0, r1 / r2
  217.         });
  218.     }

  219.     /**
  220.      * Generate the polynomials needed in the linear transformation.
  221.      *
  222.      * <p>
  223.      * See Petre's paper
  224.      * </p>
  225.      */
  226.     private void generatePolynomials() {


  227.         // Initialization of the matrices for linear transformations
  228.         // The final configuration of these matrices are obtained by composition
  229.         // of linear transformations

  230.         // The matrix of polynomials associated to Hansen coefficients, Petre's
  231.         // paper
  232.         PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix4();

  233.         // The matrix of polynomials associated to derivatives, Petre's paper
  234.         final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix4();
  235.         PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix4();
  236.         final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix4();

  237.         // The matrix of the current linear transformation
  238.         a.setMatrixLine(0, new PolynomialFunction[] {
  239.             HansenUtilities.ZERO, HansenUtilities.ONE, HansenUtilities.ZERO, HansenUtilities.ZERO
  240.         });
  241.         a.setMatrixLine(1, new PolynomialFunction[] {
  242.             HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ONE, HansenUtilities.ZERO
  243.         });
  244.         a.setMatrixLine(2, new PolynomialFunction[] {
  245.             HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ONE
  246.         });
  247.         // The generation process
  248.         int index;
  249.         int sliceCounter = 0;
  250.         for (int i = N0 - 1; i > Nmin - 1; i--) {
  251.             index = i + this.offset;
  252.             // The matrix of the current linear transformation is updated
  253.             // Petre's paper
  254.             a.setMatrixLine(3, new PolynomialFunction[] {
  255.                     c(i), HansenUtilities.ZERO, b(i), a(i)
  256.             });

  257.             // composition of the linear transformations to calculate
  258.             // the polynomials associated to Hansen coefficients
  259.             // Petre's paper
  260.             A = A.multiply(a);
  261.             // store the polynomials for Hansen coefficients
  262.             mpvec[index] = A.getMatrixLine(3);
  263.             // composition of the linear transformations to calculate
  264.             // the polynomials associated to derivatives
  265.             // Petre's paper
  266.             D = D.multiply(a);

  267.             //Update the B matrix
  268.             B.setMatrixLine(3, new PolynomialFunction[] {
  269.                 HansenUtilities.ZERO, f(i),
  270.                 HansenUtilities.ZERO, d(i)
  271.             });
  272.             D = D.add(A.multiply(B));

  273.             // store the polynomials for Hansen coefficients from the
  274.             // expressions of derivatives
  275.             mpvecDeriv[index] = D.getMatrixLine(3);

  276.             if (++sliceCounter % SLICE == 0) {
  277.                 // Re-Initialisation of matrix for linear transformmations
  278.                 // The final configuration of these matrix are obtained by composition
  279.                 // of linear transformations
  280.                 A = HansenUtilities.buildIdentityMatrix4();
  281.                 D = HansenUtilities.buildZeroMatrix4();
  282.             }
  283.         }
  284.     }

  285.     /**
  286.      * Compute the values for the first four coefficients and their derivatives by means of series.
  287.      *
  288.      * @param e2 e²
  289.      * @param chi &Chi;
  290.      * @param chi2 &Chi;²
  291.      */
  292.     public void computeInitValues(final T e2, final T chi, final T chi2) {
  293.         // compute the values for n, n+1, n+2 and n+3 by series
  294.         // See Danielson 2.7.3-(10)
  295.         //Ensure that only the needed terms are computed
  296.         final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
  297.         for (int i = 0; i < maxRoots; i++) {
  298.             final FieldGradient<T> hansenKernel = hansenInit[i].getValueGradient(e2, chi, chi2);
  299.             this.hansenRoot[0][i] = hansenKernel.getValue();
  300.             this.hansenDerivRoot[0][i] = hansenKernel.getPartialDerivative(0);
  301.         }

  302.         for (int i = 1; i < numSlices; i++) {
  303.             for (int k = 0; k < 4; k++) {
  304.                 final PolynomialFunction[] mv = mpvec[N0 - (i * SLICE) - k + 3 + offset];
  305.                 final PolynomialFunction[] sv = mpvecDeriv[N0 - (i * SLICE) - k + 3 + offset];

  306.                 hansenDerivRoot[i][k] = mv[3].value(chi).multiply(hansenDerivRoot[i - 1][3]).
  307.                                         add(mv[2].value(chi).multiply(hansenDerivRoot[i - 1][2])).
  308.                                         add(mv[1].value(chi).multiply(hansenDerivRoot[i - 1][1])).
  309.                                         add(mv[0].value(chi).multiply(hansenDerivRoot[i - 1][0])).
  310.                                         add(sv[3].value(chi).multiply(hansenRoot[i - 1][3])).
  311.                                         add(sv[2].value(chi).multiply(hansenRoot[i - 1][2])).
  312.                                         add(sv[1].value(chi).multiply(hansenRoot[i - 1][1])).
  313.                                         add(sv[0].value(chi).multiply(hansenRoot[i - 1][0]));

  314.                 hansenRoot[i][k] =  mv[3].value(chi).multiply(hansenRoot[i - 1][3]).
  315.                                     add(mv[2].value(chi).multiply(hansenRoot[i - 1][2])).
  316.                                     add(mv[1].value(chi).multiply(hansenRoot[i - 1][1])).
  317.                                     add(mv[0].value(chi).multiply(hansenRoot[i - 1][0]));
  318.             }
  319.         }
  320.     }

  321.     /**
  322.      * Compute the value of the Hansen coefficient K<sub>j</sub><sup>-n-1, s</sup>.
  323.      *
  324.      * @param mnm1 -n-1
  325.      * @param chi χ
  326.      * @return the coefficient K<sub>j</sub><sup>-n-1, s</sup>
  327.      */
  328.     public T getValue(final int mnm1, final T chi) {
  329.         //Compute n
  330.         final int n = -mnm1 - 1;

  331.         //Compute the potential slice
  332.         int sliceNo = (n + N0 + 4) / SLICE;
  333.         if (sliceNo < numSlices) {
  334.             //Compute the index within the slice
  335.             final int indexInSlice = (n + N0 + 4) % SLICE;

  336.             //Check if a root must be returned
  337.             if (indexInSlice <= 3) {
  338.                 return hansenRoot[sliceNo][indexInSlice];
  339.             }
  340.         } else {
  341.             //the value was a potential root for a slice, but that slice was not required
  342.             //Decrease the slice number
  343.             sliceNo--;
  344.         }

  345.         // Computes the coefficient by linear transformation
  346.         // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
  347.         final PolynomialFunction[] v = mpvec[mnm1 + offset];
  348.         return v[3].value(chi).multiply(hansenRoot[sliceNo][3]).
  349.                add(v[2].value(chi).multiply(hansenRoot[sliceNo][2])).
  350.                add(v[1].value(chi).multiply(hansenRoot[sliceNo][1])).
  351.                add(v[0].value(chi).multiply(hansenRoot[sliceNo][0]));

  352.     }

  353.     /**
  354.      * Compute the value of the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de².
  355.      *
  356.      * @param mnm1 -n-1
  357.      * @param chi χ
  358.      * @return the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de²
  359.      */
  360.     public T getDerivative(final int mnm1, final T chi) {

  361.         //Compute n
  362.         final int n = -mnm1 - 1;

  363.         //Compute the potential slice
  364.         int sliceNo = (n + N0 + 4) / SLICE;
  365.         if (sliceNo < numSlices) {
  366.             //Compute the index within the slice
  367.             final int indexInSlice = (n + N0 + 4) % SLICE;

  368.             //Check if a root must be returned
  369.             if (indexInSlice <= 3) {
  370.                 return hansenDerivRoot[sliceNo][indexInSlice];
  371.             }
  372.         } else {
  373.             //the value was a potential root for a slice, but that slice was not required
  374.             //Decrease the slice number
  375.             sliceNo--;
  376.         }

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

  381.         return v[3].value(chi).multiply(hansenDerivRoot[sliceNo][3]).
  382.                add(v[2].value(chi).multiply(hansenDerivRoot[sliceNo][2])).
  383.                add(v[1].value(chi).multiply(hansenDerivRoot[sliceNo][1])).
  384.                add(v[0].value(chi).multiply(hansenDerivRoot[sliceNo][0])).
  385.                add(vv[3].value(chi).multiply(hansenRoot[sliceNo][3])).
  386.                add(vv[2].value(chi).multiply(hansenRoot[sliceNo][2])).
  387.                add( vv[1].value(chi).multiply(hansenRoot[sliceNo][1])).
  388.                add(vv[0].value(chi).multiply(hansenRoot[sliceNo][0]));

  389.     }

  390.     /**
  391.      * Compute a Hansen coefficient with series.
  392.      * <p>
  393.      * This class implements the computation of the Hansen kernels
  394.      * through a power series in e² and that is using
  395.      * modified Newcomb operators. The reference formulae can be found
  396.      * in Danielson 2.7.3-10 and 3.3-5
  397.      * </p>
  398.      */
  399.     private static class FieldHansenCoefficientsBySeries <T extends CalculusFieldElement<T>> {

  400.         /** -n-1 coefficient. */
  401.         private final int mnm1;

  402.         /** s coefficient. */
  403.         private final int s;

  404.         /** j coefficient. */
  405.         private final int j;

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

  408.         /** Polynomial representing the serie. */
  409.         private PolynomialFunction polynomial;

  410.         /**
  411.          * Class constructor.
  412.          *
  413.          * @param mnm1 -n-1 value
  414.          * @param s s value
  415.          * @param j j value
  416.          * @param maxHansen max power of e² in series expansion
  417.          * @param field field for the function parameters and value
  418.          */
  419.         FieldHansenCoefficientsBySeries(final int mnm1, final int s,
  420.                                           final int j, final int maxHansen, final Field<T> field) {
  421.             this.mnm1 = mnm1;
  422.             this.s = s;
  423.             this.j = j;
  424.             this.maxNewcomb = maxHansen;
  425.             this.polynomial = generatePolynomial();
  426.         }

  427.         /** Computes the value of Hansen kernel and its derivative at e².
  428.          * <p>
  429.          * The formulae applied are described in Danielson 2.7.3-10 and
  430.          * 3.3-5
  431.          * </p>
  432.          * @param e2 e²
  433.          * @param chi &Chi;
  434.          * @param chi2 &Chi;²
  435.          * @return the value of the Hansen coefficient and its derivative for e²
  436.          */
  437.         private FieldGradient<T> getValueGradient(final T e2, final T chi, final T chi2) {

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

  441.             final T value      =  FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getValue()).divide(chi);
  442.             final T coef       = zero.subtract(mnm1 + 1.5);
  443.             final T derivative = coef.multiply(chi2).multiply(value).
  444.                             add(FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getPartialDerivative(0)).divide(chi));
  445.             return new FieldGradient<T>(value, derivative);
  446.         }

  447.         /** Generate the serie expansion in e².
  448.          * <p>
  449.          * Generate the series expansion in e² used in the formulation
  450.          * of the Hansen kernel (see Danielson 2.7.3-10):
  451.          * &Sigma; Y<sup>ns</sup><sub>α+a,α+b</sub>
  452.          * *e<sup>2α</sup>
  453.          * </p>
  454.          * @return polynomial representing the power serie expansion
  455.          */
  456.         private PolynomialFunction generatePolynomial() {
  457.             // Initialization
  458.             final int aHT = FastMath.max(j - s, 0);
  459.             final int bHT = FastMath.max(s - j, 0);

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

  461.             //Loop for getting the Newcomb operators
  462.             for (int alphaHT = 0; alphaHT <= maxNewcomb; alphaHT++) {
  463.                 coefficients[alphaHT] =
  464.                         NewcombOperators.getValue(alphaHT + aHT, alphaHT + bHT, mnm1, s);
  465.             }

  466.             //Creation of the polynomial
  467.             return new PolynomialFunction(coefficients);
  468.         }
  469.     }

  470. }