FieldHansenTesseralLinear.java

  1. /* Copyright 2002-2022 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.  */
  38. public class FieldHansenTesseralLinear <T extends CalculusFieldElement<T>> {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  105.     }

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

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

  153.     /**
  154.      * Compute polynomial coefficient c.
  155.      *
  156.      *  <p>
  157.      *  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>
  158.      *  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²
  159.      *  </p>
  160.      *
  161.      *  <p>
  162.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  163.      *  </p>
  164.      *
  165.      * @param mnm1 -n-1
  166.      * @return the polynomial
  167.      */
  168.     private PolynomialFunction c(final int mnm1) {
  169.         // Collins 4-236, Danielson 2.7.3-(9)
  170.         final double r1 = j * j * (mnm1 + 2.);
  171.         final double r2 = (mnm1 + 4.) * (2. + mnm1 + s) * (2. + mnm1 - s);

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

  176.     /**
  177.      * Compute polynomial coefficient d.
  178.      *
  179.      *  <p>
  180.      *  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²
  181.      *  </p>
  182.      *
  183.      *  <p>
  184.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  185.      *  </p>
  186.      *
  187.      * @param mnm1 -n-1
  188.      * @return the polynomial
  189.      */
  190.     private PolynomialFunction d(final int mnm1) {
  191.         // Collins 4-236, Danielson 2.7.3-(9)
  192.         return new PolynomialFunction(new double[] {
  193.             0.0, 0.0, 1.0
  194.         });
  195.     }

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  351.     }

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

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

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

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

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

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

  388.     }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  469. }