HansenTesseralLinear.java

  1. /* Copyright 2002-2024 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 org.hipparchus.analysis.differentiation.Gradient;
  19. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  20. import org.hipparchus.util.FastMath;
  21. import org.orekit.propagation.semianalytical.dsst.utilities.NewcombOperators;

  22. /**
  23.  * Hansen coefficients K(t,n,s) for t!=0 and n < 0.
  24.  * <p>
  25.  * Implements Collins 4-236 or Danielson 2.7.3-(9) for Hansen Coefficients and
  26.  * Collins 4-240 for derivatives. The recursions are transformed into
  27.  * composition of linear transformations to obtain the associated polynomials
  28.  * for coefficients and their derivatives - see Petre's paper
  29.  *
  30.  * @author Petre Bazavan
  31.  * @author Lucian Barbulescu
  32.  */
  33. public class HansenTesseralLinear {

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

  36.     /**
  37.      * The first vector of polynomials associated to Hansen coefficients and
  38.      * derivatives.
  39.      */
  40.     private PolynomialFunction[][] mpvec;

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

  43.     /** The Hansen coefficients used as roots. */
  44.     private double[][] hansenRoot;

  45.     /** The derivatives of the Hansen coefficients used as roots. */
  46.     private double[][] hansenDerivRoot;

  47.     /** The minimum value for the order. */
  48.     private int Nmin;

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

  51.     /** The s coefficient. */
  52.     private int s;

  53.     /** The j coefficient. */
  54.     private int j;

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

  57.     /**
  58.      * The offset used to identify the polynomial that corresponds to a negative.
  59.      * value of n in the internal array that starts at 0
  60.      */
  61.     private int offset;

  62.     /** The objects used to calculate initial data by means of Newcomb operators. */
  63.     private HansenCoefficientsBySeries[] hansenInit;

  64.     /**
  65.      * Constructor.
  66.      *
  67.      * @param nMax the maximum (absolute) value of n parameter
  68.      * @param s s parameter
  69.      * @param j j parameter
  70.      * @param n0 the minimum (absolute) value of n
  71.      * @param maxHansen maximum power of e2 in Hansen expansion
  72.      */
  73.     public HansenTesseralLinear(final int nMax, final int s, final int j, final int n0, final int maxHansen) {
  74.         //Initialize the fields
  75.         this.offset = nMax + 1;
  76.         this.Nmin = -nMax - 1;
  77.         this.N0 = -n0 - 4;
  78.         this.s = s;
  79.         this.j = j;

  80.         //Ensure that only the needed terms are computed
  81.         final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
  82.         this.hansenInit = new HansenCoefficientsBySeries[maxRoots];
  83.         for (int i = 0; i < maxRoots; i++) {
  84.             this.hansenInit[i] = new HansenCoefficientsBySeries(N0 - i + 3, s, j, maxHansen);
  85.         }

  86.         // The first 4 values are computed with series. No linear combination is needed.
  87.         final int size = N0 - Nmin;
  88.         this.numSlices = (int) FastMath.max(FastMath.ceil(((double) size) / SLICE), 1);
  89.         hansenRoot = new double[numSlices][4];
  90.         hansenDerivRoot = new double[numSlices][4];
  91.         if (size > 0) {
  92.             mpvec = new PolynomialFunction[size][];
  93.             mpvecDeriv = new PolynomialFunction[size][];

  94.             // Prepare the database of the associated polynomials
  95.             generatePolynomials();
  96.         }

  97.     }

  98.     /**
  99.      * Compute polynomial coefficient a.
  100.      *
  101.      *  <p>
  102.      *  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>
  103.      *  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²
  104.      *  </p>
  105.      *
  106.      *  <p>
  107.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  108.      *  </p>
  109.      *
  110.      * @param mnm1 -n-1
  111.      * @return the polynomial
  112.      */
  113.     private PolynomialFunction a(final int mnm1) {
  114.         // Collins 4-236, Danielson 2.7.3-(9)
  115.         final double r1 = (mnm1 + 2.) * (2. * mnm1 + 5.);
  116.         final double r2 = (2. + mnm1 + s) * (2. + mnm1 - s);
  117.         return new PolynomialFunction(new double[] {
  118.             0.0, 0.0, r1 / r2
  119.         });
  120.     }

  121.     /**
  122.      * Compute polynomial coefficient b.
  123.      *
  124.      *  <p>
  125.      *  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>
  126.      *  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²
  127.      *  </p>
  128.      *
  129.      *  <p>
  130.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  131.      *  </p>
  132.      *
  133.      * @param mnm1 -n-1
  134.      * @return the polynomial
  135.      */
  136.     private PolynomialFunction b(final int mnm1) {
  137.         // Collins 4-236, Danielson 2.7.3-(9)
  138.         final double r2 = (2. + mnm1 + s) * (2. + mnm1 - s);
  139.         final double d1 = (mnm1 + 3.) * 2. * j * s / (r2 * (mnm1 + 4.));
  140.         final double d2 = (mnm1 + 3.) * (mnm1 + 2.) / r2;
  141.         return new PolynomialFunction(new double[] {
  142.             0.0, -d1, -d2
  143.         });
  144.     }

  145.     /**
  146.      * Compute polynomial coefficient c.
  147.      *
  148.      *  <p>
  149.      *  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>
  150.      *  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²
  151.      *  </p>
  152.      *
  153.      *  <p>
  154.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  155.      *  </p>
  156.      *
  157.      * @param mnm1 -n-1
  158.      * @return the polynomial
  159.      */
  160.     private PolynomialFunction c(final int mnm1) {
  161.         // Collins 4-236, Danielson 2.7.3-(9)
  162.         final double r1 = j * j * (mnm1 + 2.);
  163.         final double r2 = (mnm1 + 4.) * (2. + mnm1 + s) * (2. + mnm1 - s);

  164.         return new PolynomialFunction(new double[] {
  165.             0.0, 0.0, r1 / r2
  166.         });
  167.     }

  168.     /**
  169.      * Compute polynomial coefficient d.
  170.      *
  171.      *  <p>
  172.      *  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²
  173.      *  </p>
  174.      *
  175.      *  <p>
  176.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  177.      *  </p>
  178.      *
  179.      * @param mnm1 -n-1
  180.      * @return the polynomial
  181.      */
  182.     private PolynomialFunction d(final int mnm1) {
  183.         // Collins 4-236, Danielson 2.7.3-(9)
  184.         return new PolynomialFunction(new double[] {
  185.             0.0, 0.0, 1.0
  186.         });
  187.     }

  188.     /**
  189.      * Compute polynomial coefficient f.
  190.      *
  191.      *  <p>
  192.      *  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²
  193.      *  </p>
  194.      *
  195.      *  <p>
  196.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  197.      *  </p>
  198.      *
  199.      * @param n index
  200.      * @return the polynomial
  201.      */
  202.     private PolynomialFunction f(final int n) {
  203.         // Collins 4-236, Danielson 2.7.3-(9)
  204.         final double r1 = (n + 3.0) * j * s;
  205.         final double r2 = (n + 4.0) * (2.0 + n + s) * (2.0 + n - s);
  206.         return new PolynomialFunction(new double[] {
  207.             0.0, 0.0, 0.0, r1 / r2
  208.         });
  209.     }

  210.     /**
  211.      * Generate the polynomials needed in the linear transformation.
  212.      *
  213.      * <p>
  214.      * See Petre's paper
  215.      * </p>
  216.      */
  217.     private void generatePolynomials() {


  218.         // Initialization of the matrices for linear transformations
  219.         // The final configuration of these matrices are obtained by composition
  220.         // of linear transformations

  221.         // The matrix of polynomials associated to Hansen coefficients, Petre's
  222.         // paper
  223.         PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix4();

  224.         // The matrix of polynomials associated to derivatives, Petre's paper
  225.         final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix4();
  226.         PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix4();
  227.         final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix4();

  228.         // The matrix of the current linear transformation
  229.         a.setMatrixLine(0, new PolynomialFunction[] {
  230.             HansenUtilities.ZERO, HansenUtilities.ONE, HansenUtilities.ZERO, HansenUtilities.ZERO
  231.         });
  232.         a.setMatrixLine(1, new PolynomialFunction[] {
  233.             HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ONE, HansenUtilities.ZERO
  234.         });
  235.         a.setMatrixLine(2, new PolynomialFunction[] {
  236.             HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ONE
  237.         });
  238.         // The generation process
  239.         int index;
  240.         int sliceCounter = 0;
  241.         for (int i = N0 - 1; i > Nmin - 1; i--) {
  242.             index = i + this.offset;
  243.             // The matrix of the current linear transformation is updated
  244.             // Petre's paper
  245.             a.setMatrixLine(3, new PolynomialFunction[] {
  246.                     c(i), HansenUtilities.ZERO, b(i), a(i)
  247.             });

  248.             // composition of the linear transformations to calculate
  249.             // the polynomials associated to Hansen coefficients
  250.             // Petre's paper
  251.             A = A.multiply(a);
  252.             // store the polynomials for Hansen coefficients
  253.             mpvec[index] = A.getMatrixLine(3);
  254.             // composition of the linear transformations to calculate
  255.             // the polynomials associated to derivatives
  256.             // Petre's paper
  257.             D = D.multiply(a);

  258.             //Update the B matrix
  259.             B.setMatrixLine(3, new PolynomialFunction[] {
  260.                 HansenUtilities.ZERO, f(i),
  261.                 HansenUtilities.ZERO, d(i)
  262.             });
  263.             D = D.add(A.multiply(B));

  264.             // store the polynomials for Hansen coefficients from the
  265.             // expressions of derivatives
  266.             mpvecDeriv[index] = D.getMatrixLine(3);

  267.             if (++sliceCounter % SLICE == 0) {
  268.                 // Re-Initialisation of matrix for linear transformmations
  269.                 // The final configuration of these matrix are obtained by composition
  270.                 // of linear transformations
  271.                 A = HansenUtilities.buildIdentityMatrix4();
  272.                 D = HansenUtilities.buildZeroMatrix4();
  273.             }
  274.         }
  275.     }

  276.     /**
  277.      * Compute the values for the first four coefficients and their derivatives by means of series.
  278.      *
  279.      * @param e2 e²
  280.      * @param chi &Chi;
  281.      * @param chi2 &Chi;²
  282.      */
  283.     public void computeInitValues(final double e2, final double chi, final double chi2) {
  284.         // compute the values for n, n+1, n+2 and n+3 by series
  285.         // See Danielson 2.7.3-(10)
  286.         //Ensure that only the needed terms are computed
  287.         final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
  288.         for (int i = 0; i < maxRoots; i++) {
  289.             final Gradient hansenKernel = hansenInit[i].getValueGradient(e2, chi, chi2);
  290.             this.hansenRoot[0][i] = hansenKernel.getValue();
  291.             this.hansenDerivRoot[0][i] = hansenKernel.getPartialDerivative(0);
  292.         }

  293.         for (int i = 1; i < numSlices; i++) {
  294.             for (int k = 0; k < 4; k++) {
  295.                 final PolynomialFunction[] mv = mpvec[N0 - (i * SLICE) - k + 3 + offset];
  296.                 final PolynomialFunction[] sv = mpvecDeriv[N0 - (i * SLICE) - k + 3 + offset];

  297.                 hansenDerivRoot[i][k] = mv[3].value(chi) * hansenDerivRoot[i - 1][3] +
  298.                                         mv[2].value(chi) * hansenDerivRoot[i - 1][2] +
  299.                                         mv[1].value(chi) * hansenDerivRoot[i - 1][1] +
  300.                                         mv[0].value(chi) * hansenDerivRoot[i - 1][0] +
  301.                                         sv[3].value(chi) * hansenRoot[i - 1][3] +
  302.                                         sv[2].value(chi) * hansenRoot[i - 1][2] +
  303.                                         sv[1].value(chi) * hansenRoot[i - 1][1] +
  304.                                         sv[0].value(chi) * hansenRoot[i - 1][0];

  305.                 hansenRoot[i][k] =  mv[3].value(chi) * hansenRoot[i - 1][3] +
  306.                                     mv[2].value(chi) * hansenRoot[i - 1][2] +
  307.                                     mv[1].value(chi) * hansenRoot[i - 1][1] +
  308.                                     mv[0].value(chi) * hansenRoot[i - 1][0];
  309.             }
  310.         }
  311.     }

  312.     /**
  313.      * Compute the value of the Hansen coefficient K<sub>j</sub><sup>-n-1, s</sup>.
  314.      *
  315.      * @param mnm1 -n-1
  316.      * @param chi χ
  317.      * @return the coefficient K<sub>j</sub><sup>-n-1, s</sup>
  318.      */
  319.     public double getValue(final int mnm1, final double chi) {
  320.         //Compute n
  321.         final int n = -mnm1 - 1;

  322.         //Compute the potential slice
  323.         int sliceNo = (n + N0 + 4) / SLICE;
  324.         if (sliceNo < numSlices) {
  325.             //Compute the index within the slice
  326.             final int indexInSlice = (n + N0 + 4) % SLICE;

  327.             //Check if a root must be returned
  328.             if (indexInSlice <= 3) {
  329.                 return hansenRoot[sliceNo][indexInSlice];
  330.             }
  331.         } else {
  332.             //the value was a potential root for a slice, but that slice was not required
  333.             //Decrease the slice number
  334.             sliceNo--;
  335.         }

  336.         // Computes the coefficient by linear transformation
  337.         // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
  338.         final PolynomialFunction[] v = mpvec[mnm1 + offset];
  339.         return v[3].value(chi) * hansenRoot[sliceNo][3] +
  340.                v[2].value(chi) * hansenRoot[sliceNo][2] +
  341.                v[1].value(chi) * hansenRoot[sliceNo][1] +
  342.                v[0].value(chi) * hansenRoot[sliceNo][0];

  343.     }

  344.     /**
  345.      * Compute the value of the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de².
  346.      *
  347.      * @param mnm1 -n-1
  348.      * @param chi χ
  349.      * @return the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de²
  350.      */
  351.     public double getDerivative(final int mnm1, final double chi) {

  352.         //Compute n
  353.         final int n = -mnm1 - 1;

  354.         //Compute the potential slice
  355.         int sliceNo = (n + N0 + 4) / SLICE;
  356.         if (sliceNo < numSlices) {
  357.             //Compute the index within the slice
  358.             final int indexInSlice = (n + N0 + 4) % SLICE;

  359.             //Check if a root must be returned
  360.             if (indexInSlice <= 3) {
  361.                 return hansenDerivRoot[sliceNo][indexInSlice];
  362.             }
  363.         } else {
  364.             //the value was a potential root for a slice, but that slice was not required
  365.             //Decrease the slice number
  366.             sliceNo--;
  367.         }

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

  372.         return v[3].value(chi)  * hansenDerivRoot[sliceNo][3] +
  373.                v[2].value(chi)  * hansenDerivRoot[sliceNo][2] +
  374.                v[1].value(chi)  * hansenDerivRoot[sliceNo][1] +
  375.                v[0].value(chi)  * hansenDerivRoot[sliceNo][0] +
  376.                vv[3].value(chi) * hansenRoot[sliceNo][3] +
  377.                vv[2].value(chi) * hansenRoot[sliceNo][2] +
  378.                vv[1].value(chi) * hansenRoot[sliceNo][1] +
  379.                vv[0].value(chi) * hansenRoot[sliceNo][0];

  380.     }

  381.     /**
  382.      * Compute a Hansen coefficient with series.
  383.      * <p>
  384.      * This class implements the computation of the Hansen kernels
  385.      * through a power series in e² and that is using
  386.      * modified Newcomb operators. The reference formulae can be found
  387.      * in Danielson 2.7.3-10 and 3.3-5
  388.      * </p>
  389.      */
  390.     private static class HansenCoefficientsBySeries {

  391.         /** -n-1 coefficient. */
  392.         private final int mnm1;

  393.         /** s coefficient. */
  394.         private final int s;

  395.         /** j coefficient. */
  396.         private final int j;

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

  399.         /** Polynomial representing the serie. */
  400.         private PolynomialFunction polynomial;

  401.         /**
  402.          * Class constructor.
  403.          *
  404.          * @param mnm1 -n-1 value
  405.          * @param s s value
  406.          * @param j j value
  407.          * @param maxHansen max power of e² in series expansion
  408.          */
  409.         HansenCoefficientsBySeries(final int mnm1, final int s,
  410.                                           final int j, final int maxHansen) {
  411.             this.mnm1 = mnm1;
  412.             this.s = s;
  413.             this.j = j;
  414.             this.maxNewcomb = maxHansen;
  415.             this.polynomial = generatePolynomial();
  416.         }

  417.         /** Computes the value of Hansen kernel and its derivative at e².
  418.          * <p>
  419.          * The formulae applied are described in Danielson 2.7.3-10 and
  420.          * 3.3-5
  421.          * </p>
  422.          * @param e2 e²
  423.          * @param chi &Chi;
  424.          * @param chi2 &Chi;²
  425.          * @return the value of the Hansen coefficient and its derivative for e²
  426.          */
  427.         public Gradient getValueGradient(final double e2, final double chi, final double chi2) {

  428.             //Estimation of the serie expansion at e2
  429.             final Gradient serie = polynomial.value(Gradient.variable(1, 0, e2));

  430.             final double value      =  FastMath.pow(chi2, -mnm1 - 1) * serie.getValue() / chi;
  431.             final double coef       = -(mnm1 + 1.5);
  432.             final double derivative = coef * chi2 * value +
  433.                                       FastMath.pow(chi2, -mnm1 - 1) * serie.getPartialDerivative(0) / chi;
  434.             return new Gradient(value, derivative);
  435.         }

  436.         /** Generate the serie expansion in e².
  437.          * <p>
  438.          * Generate the series expansion in e² used in the formulation
  439.          * of the Hansen kernel (see Danielson 2.7.3-10):
  440.          * &Sigma; Y<sup>ns</sup><sub>α+a,α+b</sub>
  441.          * *e<sup>2α</sup>
  442.          * </p>
  443.          * @return polynomial representing the power serie expansion
  444.          */
  445.         private PolynomialFunction generatePolynomial() {
  446.             // Initialization
  447.             final int aHT = FastMath.max(j - s, 0);
  448.             final int bHT = FastMath.max(s - j, 0);

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

  450.             //Loop for getting the Newcomb operators
  451.             for (int alphaHT = 0; alphaHT <= maxNewcomb; alphaHT++) {
  452.                 coefficients[alphaHT] =
  453.                         NewcombOperators.getValue(alphaHT + aHT, alphaHT + bHT, mnm1, s);
  454.             }

  455.             //Creation of the polynomial
  456.             return new PolynomialFunction(coefficients);
  457.         }
  458.     }

  459. }