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;

  19. import org.hipparchus.analysis.differentiation.Gradient;
  20. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  21. import org.hipparchus.util.FastMath;
  22. import org.orekit.propagation.semianalytical.dsst.utilities.NewcombOperators;

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

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

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

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

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

  52.     /** The derivatives of the Hansen coefficients used as roots. */
  53.     private final double[][] hansenDerivRoot;

  55.     /** The minimum value for the order. */
  56.     private final int Nmin;

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

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

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

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

  73.     /**
  74.      * Constructor.
  75.      *
  76.      * @param nMax the maximum (absolute) value of n parameter
  77.      * @param s s parameter
  78.      * @param j j parameter
  79.      * @param n0 the minimum (absolute) value of n
  80.      * @param maxHansen maximum power of e2 in Hansen expansion
  81.      */
  82.     public HansenTesseralLinear(final int nMax, final int s, final int j, final int n0, final int maxHansen) {
  83.         //Initialize the fields
  84.         this.offset = nMax + 1;
  85.         this.Nmin = -nMax - 1;
  86.         this.N0 = -n0 - 4;

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

  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 = new double[numSlices][4];
  99.         hansenDerivRoot = new double[numSlices][4];
  100.         if (size > 0) {
  101.             mpvec = new PolynomialFunction[size][];
  102.             mpvecDeriv = new PolynomialFunction[size][];

  104.             // Prepare the database of the associated polynomials
  105.             HansenUtilities.generateTesseralPolynomials(N0, Nmin, offset, SLICE, j, s,
  106.                                                         mpvec, mpvecDeriv);
  107.         }

  109.     }

  111.     /**
  112.      * Compute the values for the first four coefficients and their derivatives by means of series.
  113.      *
  114.      * @param e2 e²
  115.      * @param chi &Chi;
  116.      * @param chi2 &Chi;²
  117.      */
  118.     public void computeInitValues(final double e2, final double chi, final double chi2) {
  119.         // compute the values for n, n+1, n+2 and n+3 by series
  120.         // See Danielson 2.7.3-(10)
  121.         //Ensure that only the needed terms are computed
  122.         final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
  123.         for (int i = 0; i < maxRoots; i++) {
  124.             final Gradient hansenKernel = hansenInit[i].getValueGradient(e2, chi, chi2);
  125.             this.hansenRoot[0][i] = hansenKernel.getValue();
  126.             this.hansenDerivRoot[0][i] = hansenKernel.getPartialDerivative(0);
  127.         }

  129.         for (int i = 1; i < numSlices; i++) {
  130.             for (int k = 0; k < 4; k++) {
  131.                 final PolynomialFunction[] mv = mpvec[N0 - (i * SLICE) - k + 3 + offset];
  132.                 final PolynomialFunction[] sv = mpvecDeriv[N0 - (i * SLICE) - k + 3 + offset];

  134.                 hansenDerivRoot[i][k] = mv[3].value(chi) * hansenDerivRoot[i - 1][3] +
  135.                                         mv[2].value(chi) * hansenDerivRoot[i - 1][2] +
  136.                                         mv[1].value(chi) * hansenDerivRoot[i - 1][1] +
  137.                                         mv[0].value(chi) * hansenDerivRoot[i - 1][0] +
  138.                                         sv[3].value(chi) * hansenRoot[i - 1][3] +
  139.                                         sv[2].value(chi) * hansenRoot[i - 1][2] +
  140.                                         sv[1].value(chi) * hansenRoot[i - 1][1] +
  141.                                         sv[0].value(chi) * hansenRoot[i - 1][0];

  143.                 hansenRoot[i][k] =  mv[3].value(chi) * hansenRoot[i - 1][3] +
  144.                                     mv[2].value(chi) * hansenRoot[i - 1][2] +
  145.                                     mv[1].value(chi) * hansenRoot[i - 1][1] +
  146.                                     mv[0].value(chi) * hansenRoot[i - 1][0];
  147.             }
  148.         }
  149.     }

  151.     /**
  152.      * Compute the value of the Hansen coefficient K<sub>j</sub><sup>-n-1, s</sup>.
  153.      *
  154.      * @param mnm1 -n-1
  155.      * @param chi χ
  156.      * @return the coefficient K<sub>j</sub><sup>-n-1, s</sup>
  157.      */
  158.     public double getValue(final int mnm1, final double chi) {
  159.         //Compute n
  160.         final int n = -mnm1 - 1;

  162.         //Compute the potential slice
  163.         int sliceNo = (n + N0 + 4) / SLICE;
  164.         if (sliceNo < numSlices) {
  165.             //Compute the index within the slice
  166.             final int indexInSlice = (n + N0 + 4) % SLICE;

  168.             //Check if a root must be returned
  169.             if (indexInSlice <= 3) {
  170.                 return hansenRoot[sliceNo][indexInSlice];
  171.             }
  172.         } else {
  173.             //the value was a potential root for a slice, but that slice was not required
  174.             //Decrease the slice number
  175.             sliceNo--;
  176.         }

  178.         // Computes the coefficient by linear transformation
  179.         // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
  180.         final PolynomialFunction[] v = mpvec[mnm1 + offset];
  181.         return v[3].value(chi) * hansenRoot[sliceNo][3] +
  182.                v[2].value(chi) * hansenRoot[sliceNo][2] +
  183.                v[1].value(chi) * hansenRoot[sliceNo][1] +
  184.                v[0].value(chi) * hansenRoot[sliceNo][0];

  186.     }

  188.     /**
  189.      * Compute the value of the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de².
  190.      *
  191.      * @param mnm1 -n-1
  192.      * @param chi χ
  193.      * @return the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de²
  194.      */
  195.     public double getDerivative(final int mnm1, final double chi) {

  197.         //Compute n
  198.         final int n = -mnm1 - 1;

  200.         //Compute the potential slice
  201.         int sliceNo = (n + N0 + 4) / SLICE;
  202.         if (sliceNo < numSlices) {
  203.             //Compute the index within the slice
  204.             final int indexInSlice = (n + N0 + 4) % SLICE;

  206.             //Check if a root must be returned
  207.             if (indexInSlice <= 3) {
  208.                 return hansenDerivRoot[sliceNo][indexInSlice];
  209.             }
  210.         } else {
  211.             //the value was a potential root for a slice, but that slice was not required
  212.             //Decrease the slice number
  213.             sliceNo--;
  214.         }

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

  221.         return v[3].value(chi)  * hansenDerivRoot[sliceNo][3] +
  222.                v[2].value(chi)  * hansenDerivRoot[sliceNo][2] +
  223.                v[1].value(chi)  * hansenDerivRoot[sliceNo][1] +
  224.                v[0].value(chi)  * hansenDerivRoot[sliceNo][0] +
  225.                vv[3].value(chi) * hansenRoot[sliceNo][3] +
  226.                vv[2].value(chi) * hansenRoot[sliceNo][2] +
  227.                vv[1].value(chi) * hansenRoot[sliceNo][1] +
  228.                vv[0].value(chi) * hansenRoot[sliceNo][0];

  230.     }

  232.     /**
  233.      * Compute a Hansen coefficient with series.
  234.      * <p>
  235.      * This class implements the computation of the Hansen kernels
  236.      * through a power series in e² and that is using
  237.      * modified Newcomb operators. The reference formulae can be found
  238.      * in Danielson 2.7.3-10 and 3.3-5
  239.      * </p>
  240.      */
  241.     private static class HansenCoefficientsBySeries {

  243.         /** -n-1 coefficient. */
  244.         private final int mnm1;

  246.         /** s coefficient. */
  247.         private final int s;

  249.         /** j coefficient. */
  250.         private final int j;

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

  255.         /** Polynomial representing the serie. */
  256.         private final PolynomialFunction polynomial;

  258.         /**
  259.          * Class constructor.
  260.          *
  261.          * @param mnm1 -n-1 value
  262.          * @param s s value
  263.          * @param j j value
  264.          * @param maxHansen max power of e² in series expansion
  265.          */
  266.         HansenCoefficientsBySeries(final int mnm1, final int s,
  267.                                           final int j, final int maxHansen) {
  268.             this.mnm1 = mnm1;
  269.             this.s = s;
  270.             this.j = j;
  271.             this.maxNewcomb = maxHansen;
  272.             this.polynomial = generatePolynomial();
  273.         }

  275.         /** Computes the value of Hansen kernel and its derivative at e².
  276.          * <p>
  277.          * The formulae applied are described in Danielson 2.7.3-10 and
  278.          * 3.3-5
  279.          * </p>
  280.          * @param e2 e²
  281.          * @param chi &Chi;
  282.          * @param chi2 &Chi;²
  283.          * @return the value of the Hansen coefficient and its derivative for e²
  284.          */
  285.         public Gradient getValueGradient(final double e2, final double chi, final double chi2) {

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

  290.             final double value      =  FastMath.pow(chi2, -mnm1 - 1) * serie.getValue() / chi;
  291.             final double coef       = -(mnm1 + 1.5);
  292.             final double derivative = coef * chi2 * value +
  293.                                       FastMath.pow(chi2, -mnm1 - 1) * serie.getPartialDerivative(0) / chi;
  294.             return new Gradient(value, derivative);
  295.         }

  297.         /** Generate the serie expansion in e².
  298.          * <p>
  299.          * Generate the series expansion in e² used in the formulation
  300.          * of the Hansen kernel (see Danielson 2.7.3-10):
  301.          * &Sigma; Y<sup>ns</sup><sub>α+a,α+b</sub>
  302.          * *e<sup>2α</sup>
  303.          * </p>
  304.          * @return polynomial representing the power serie expansion
  305.          */
  306.         private PolynomialFunction generatePolynomial() {
  307.             // Initialization
  308.             final int aHT = FastMath.max(j - s, 0);
  309.             final int bHT = FastMath.max(s - j, 0);

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

  313.             //Loop for getting the Newcomb operators
  314.             for (int alphaHT = 0; alphaHT <= maxNewcomb; alphaHT++) {
  315.                 coefficients[alphaHT] =
  316.                         NewcombOperators.getValue(alphaHT + aHT, alphaHT + bHT, mnm1, s);
  317.             }

  319.             //Creation of the polynomial
  320.             return new PolynomialFunction(coefficients);
  321.         }
  322.     }

  324. }
Created with JaCoCo 0.8.12.202403310830