FieldHansenZonalLinear.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.Field;
  20. import org.hipparchus.CalculusFieldElement;
  21. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  22. import org.hipparchus.util.FastMath;
  23. import org.hipparchus.util.MathArrays;

  25. /**
  26.  * Hansen coefficients K(t,n,s) for t=0 and n < 0.
  27.  * <p>
  28.  *Implements Collins 4-242 or echivalently, Danielson 2.7.3-(6) for Hansen Coefficients and
  29.  * Collins 4-245 or Danielson 3.1-(7) for derivatives. The recursions are transformed into
  30.  * composition of linear transformations to obtain the associated polynomials
  31.  * for coefficients and their derivatives - see Petre's paper
  32.  *
  33.  * @author Petre Bazavan
  34.  * @author Lucian Barbulescu
  35.  * @author Bryan Cazabonne
  36.  * @param <T> type of the field elements
  37.  */
  38. public class FieldHansenZonalLinear <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;

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

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

  52.     /** The Hansen coefficients used as roots. */
  53.     private T[][] hansenRoot;

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

  58.     /** The minimum value for the order. */
  59.     private int Nmin;


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

  65.     /** The s coefficient. */
  66.     private int s;

  68.     /**
  69.      * The offset used to identify the polynomial that corresponds to a negative
  70.      * value of n in the internal array that starts at 0.
  71.      */
  72.     private int offset;

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

  77.     /** 2<sup>s</sup>. */
  78.     private double twots;

  80.     /** 2*s+1. */
  81.     private int twosp1;

  83.     /** 2*s. */
  84.     private int twos;

  86.     /** (2*s+1) / 2<sup>s</sup>. */
  87.     private double twosp1otwots;

  89.     /**
  90.      * Constructor.
  91.      *
  92.      * @param nMax the maximum (absolute) value of n coefficient
  93.      * @param s s coefficient
  94.      * @param field field used by default
  95.      */
  96.     public FieldHansenZonalLinear(final int nMax, final int s, final Field<T> field) {
  97.         //Initialize fields
  98.         this.offset = nMax + 1;
  99.         this.Nmin = -nMax - 1;
  100.         N0 = -(s + 2);
  101.         this.s = s;
  102.         this.twots = FastMath.pow(2., s);
  103.         this.twos = 2 * s;
  104.         this.twosp1 = this.twos + 1;
  105.         this.twosp1otwots = (double) this.twosp1 / this.twots;

  107.         // prepare structures for stored data
  108.         final int size = nMax - s - 1;
  109.         mpvec      = new PolynomialFunction[size][];
  110.         mpvecDeriv = new PolynomialFunction[size][];

  112.         this.numSlices  = FastMath.max((int) FastMath.ceil(((double) size) / SLICE), 1);
  113.         hansenRoot      = MathArrays.buildArray(field, numSlices, 2);
  114.         hansenDerivRoot = MathArrays.buildArray(field, numSlices, 2);

  116.         // Prepare the data base of associated polynomials
  117.         generatePolynomials();

  119.     }

  121.     /**
  122.      * Compute polynomial coefficient a.
  123.      *
  124.      *  <p>
  125.      *  It is used to generate the coefficient for K₀<sup>-n, s</sup> when computing K₀<sup>-n-1, s</sup>
  126.      *  and the coefficient for dK₀<sup>-n, s</sup> / de² when computing dK₀<sup>-n-1, s</sup> / de²
  127.      *  </p>
  128.      *
  129.      *  <p>
  130.      *  See Danielson 2.7.3-(6) and Collins 4-242 and 4-245
  131.      *  </p>
  132.      *
  133.      * @param mnm1 -n-1 value
  134.      * @return the polynomial
  135.      */
  136.     private PolynomialFunction a(final int mnm1) {
  137.         // from recurence Collins 4-242
  138.         final double d1 = (mnm1 + 2) * (2 * mnm1 + 5);
  139.         final double d2 = (mnm1 + 2 - s) * (mnm1 + 2 + s);
  140.         return new PolynomialFunction(new double[] {
  141.             0.0, 0.0, d1 / d2
  142.         });
  143.     }

  145.     /**
  146.      * Compute polynomial coefficient b.
  147.      *
  148.      *  <p>
  149.      *  It is used to generate the coefficient for K₀<sup>-n+1, s</sup> when computing K₀<sup>-n-1, s</sup>
  150.      *  and the coefficient for dK₀<sup>-n+1, s</sup> / de² when computing dK₀<sup>-n-1, s</sup> / de²
  151.      *  </p>
  152.      *
  153.      *  <p>
  154.      *  See Danielson 2.7.3-(6) and Collins 4-242 and 4-245
  155.      *  </p>
  156.      *
  157.      * @param mnm1 -n-1 value
  158.      * @return the polynomial
  159.      */
  160.     private PolynomialFunction b(final int mnm1) {
  161.         // from recurence Collins 4-242
  162.         final double d1 = (mnm1 + 2) * (mnm1 + 3);
  163.         final double d2 = (mnm1 + 2 - s) * (mnm1 + 2 + s);
  164.         return new PolynomialFunction(new double[] {
  165.             0.0, 0.0, -d1 / d2
  166.         });
  167.     }

  169.     /**
  170.      * Generate the polynomials needed in the linear transformation.
  171.      *
  172.      * <p>
  173.      * See Petre's paper
  174.      * </p>
  175.      */
  176.     private void generatePolynomials() {

  178.         int sliceCounter = 0;
  179.         int index;

  181.         // Initialisation of matrix for linear transformmations
  182.         // The final configuration of these matrix are obtained by composition
  183.         // of linear transformations
  184.         PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix2();
  185.         PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix2();
  186.         PolynomialFunctionMatrix E = HansenUtilities.buildIdentityMatrix2();

  188.         // generation of polynomials associated to Hansen coefficients and to
  189.         // their derivatives
  190.         final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix2();
  191.         a.setElem(0, 1, HansenUtilities.ONE);

  193.         //The B matrix is constant.
  194.         final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix2();
  195.         // from Collins 4-245 and Petre's paper
  196.         B.setElem(1, 1, new PolynomialFunction(new double[] {
  197.             2.0
  198.         }));

  200.         for (int i = N0 - 1; i > Nmin - 1; i--) {
  201.             index = i + offset;
  202.             // Matrix of the current linear transformation
  203.             // Petre's paper
  204.             a.setMatrixLine(1, new PolynomialFunction[] {
  205.                 b(i), a(i)
  206.             });
  207.             // composition of linear transformations
  208.             // see Petre's paper
  209.             A = A.multiply(a);
  210.             // store the polynomials for Hansen coefficients
  211.             mpvec[index] = A.getMatrixLine(1);

  213.             D = D.multiply(a);
  214.             E = E.multiply(a);
  215.             D = D.add(E.multiply(B));

  217.             // store the polynomials for Hansen coefficients from the expressions
  218.             // of derivatives
  219.             mpvecDeriv[index] = D.getMatrixLine(1);

  221.             if (++sliceCounter % SLICE == 0) {
  222.                 // Re-Initialisation of matrix for linear transformmations
  223.                 // The final configuration of these matrix are obtained by composition
  224.                 // of linear transformations
  225.                 A = HansenUtilities.buildIdentityMatrix2();
  226.                 D = HansenUtilities.buildZeroMatrix2();
  227.                 E = HansenUtilities.buildIdentityMatrix2();
  228.             }

  230.         }
  231.     }

  233.     /**
  234.      * Compute the roots for the Hansen coefficients and their derivatives.
  235.      *
  236.      * @param chi 1 / sqrt(1 - e²)
  237.      */
  238.     public void computeInitValues(final T chi) {
  239.         final Field<T> field = chi.getField();
  240.         final T zero = field.getZero();
  241.         // compute the values for n=s and n=s+1
  242.         // See Danielson 2.7.3-(6a,b)
  243.         hansenRoot[0][0] = zero;
  244.         hansenRoot[0][1] = FastMath.pow(chi, this.twosp1).divide(this.twots);
  245.         hansenDerivRoot[0][0] = zero;
  246.         hansenDerivRoot[0][1] = FastMath.pow(chi, this.twos).multiply(this.twosp1otwots);

  248.         final int st = -s - 1;
  249.         for (int i = 1; i < numSlices; i++) {
  250.             for (int j = 0; j < 2; j++) {
  251.                 // Get the required polynomials
  252.                 final PolynomialFunction[] mv = mpvec[st - (i * SLICE) - j + offset];
  253.                 final PolynomialFunction[] sv = mpvecDeriv[st - (i * SLICE) - j + offset];

  255.                 //Compute the root derivatives
  256.                 hansenDerivRoot[i][j] = mv[1].value(chi).multiply(hansenDerivRoot[i - 1][1]).
  257.                                         add(mv[0].value(chi).multiply(hansenDerivRoot[i - 1][0])).
  258.                                         add((sv[1].value(chi).multiply(hansenRoot[i - 1][1]).
  259.                                         add(sv[0].value(chi).multiply(hansenRoot[i - 1][0]))
  260.                                         ).divide(chi));
  261.                 hansenRoot[i][j] =     mv[1].value(chi).multiply(hansenRoot[i - 1][1]).
  262.                                        add(mv[0].value(chi).multiply(hansenRoot[i - 1][0]));

  264.             }

  266.         }
  267.     }

  269.     /**
  270.      * Get the K₀<sup>-n-1,s</sup> coefficient value.
  271.      *
  272.      * <p> The s value is given in the class constructor
  273.      *
  274.      * @param mnm1 (-n-1) coefficient
  275.      * @param chi The value of χ
  276.      * @return K₀<sup>-n-1,s</sup>
  277.      */
  278.     public T getValue(final int mnm1, final T chi) {
  279.         //Compute n
  280.         final int n = -mnm1 - 1;

  282.         //Compute the potential slice
  283.         int sliceNo = (n - s) / SLICE;
  284.         if (sliceNo < numSlices) {
  285.             //Compute the index within the slice
  286.             final int indexInSlice = (n - s) % SLICE;

  288.             //Check if a root must be returned
  289.             if (indexInSlice <= 1) {
  290.                 return hansenRoot[sliceNo][indexInSlice];
  291.             }
  292.         } else {
  293.             //the value was a potential root for a slice, but that slice was not required
  294.             //Decrease the slice number
  295.             sliceNo--;
  296.         }

  298.         // Danielson 2.7.3-(6c)/Collins 4-242 and Petre's paper
  299.         final PolynomialFunction[] v = mpvec[mnm1 + offset];
  300.         T ret = v[1].value(chi).multiply(hansenRoot[sliceNo][1]);
  301.         if (hansenRoot[sliceNo][0].getReal() != 0) {
  302.             ret = ret.add(v[0].value(chi).multiply(hansenRoot[sliceNo][0]));
  303.         }
  304.         return  ret;
  305.     }

  307.     /**
  308.      * Get the dK₀<sup>-n-1,s</sup> / d&Chi; coefficient derivative.
  309.      *
  310.      * <p> The s value is given in the class constructor.
  311.      *
  312.      * @param mnm1 (-n-1) coefficient
  313.      * @param chi The value of χ
  314.      * @return dK₀<sup>-n-1,s</sup> / d&Chi;
  315.      */
  316.     public T getDerivative(final int mnm1, final T chi) {
  317.         //Compute n
  318.         final int n = -mnm1 - 1;

  320.         //Compute the potential slice
  321.         int sliceNo = (n - s) / SLICE;
  322.         if (sliceNo < numSlices) {
  323.             //Compute the index within the slice
  324.             final int indexInSlice = (n - s) % SLICE;

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

  336.         // Danielson 3.1-(7c) and Petre's paper
  337.         final PolynomialFunction[] v = mpvec[mnm1 + offset];
  338.         T ret = v[1].value(chi).multiply(hansenDerivRoot[sliceNo][1]);
  339.         if (hansenDerivRoot[sliceNo][0].getReal() != 0) {
  340.             ret = ret.add(v[0].value(chi).multiply(hansenDerivRoot[sliceNo][0]));
  341.         }

  343.         // Danielson 2.7.3-(6b)
  344.         final PolynomialFunction[] v1 = mpvecDeriv[mnm1 + offset];
  345.         T hret = v1[1].value(chi).multiply(hansenRoot[sliceNo][1]);
  346.         if (hansenRoot[sliceNo][0].getReal() != 0) {
  347.             hret = hret.add(v1[0].value(chi).multiply(hansenRoot[sliceNo][0]));
  348.         }
  349.         ret = ret.add(hret.divide(chi));

  351.         return ret;

  353.     }

  355. }
Created with JaCoCo 0.8.11.202310140853