FieldHansenThirdBodyLinear.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;

  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-254 or Danielson 2.7.3-(7) for Hansen Coefficients and
  29.  * Danielson 3.2-(3) 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.  */
  37. public class FieldHansenThirdBodyLinear <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;

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

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

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

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

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

  60.     /** The maximum order of n indexes. */
  61.     private int nMax;

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

  66.     /** The s index. */
  67.     private int s;

  69.     /** (-1)<sup>s</sup> * (2*s + 1)!! / (s+1)!  */
  70.     private double twosp1dfosp1f;

  72.     /** (-1)<sup>s</sup> * (2*s + 1)!! / (s+2)!  */
  73.     private double twosp1dfosp2f;

  75.     /** (-1)<sup>s</sup> * 2 * (2*s + 1)!! / (s+2)!  */
  76.     private double two2sp1dfosp2f;

  78.     /** (2*s + 3). */
  79.     private double twosp3;

  81.     /**
  82.      * Constructor.
  83.      *
  84.      * @param nMax the maximum value of n
  85.      * @param s the value of s
  86.      * @param field field used by default
  87.      */
  88.     public FieldHansenThirdBodyLinear(final int nMax, final int s, final Field<T> field) {
  89.         // initialise fields
  90.         this.nMax = nMax;
  91.         N0 = s;
  92.         this.s = s;

  94.         // initialization of structures for stored data
  95.         mpvec = new PolynomialFunction[this.nMax + 1][];
  96.         mpvecDeriv = new PolynomialFunction[this.nMax + 1][];

  98.         //Compute the fields that will be used to determine the initial values for the coefficients
  99.         this.twosp1dfosp1f = (s % 2 == 0) ? 1.0 : -1.0;
  100.         for (int i = s; i >= 1; i--) {
  101.             this.twosp1dfosp1f *= (2.0 * i + 1.0) / (i + 1.0);
  102.         }

  104.         this.twosp1dfosp2f = this.twosp1dfosp1f / (s + 2.);
  105.         this.twosp3 = 2 * s + 3;
  106.         this.two2sp1dfosp2f = 2 * this.twosp1dfosp2f;

  108.         // initialization of structures for stored data
  109.         mpvec = new PolynomialFunction[this.nMax + 1][];
  110.         mpvecDeriv = new PolynomialFunction[this.nMax + 1][];

  112.         this.numSlices  = FastMath.max(1, (nMax - s + SLICE - 2) / SLICE);

  114.         hansenRoot      = MathArrays.buildArray(field, numSlices, 2);
  115.         hansenDerivRoot = MathArrays.buildArray(field, numSlices, 2);

  117.         // Prepare the database of the associated polynomials
  118.         generatePolynomials();

  120.     }

  122.     /**
  123.      * Compute polynomial coefficient a.
  124.      *
  125.      *  <p>
  126.      *  It is used to generate the coefficient for K₀<sup>n-1, s</sup> when computing K₀<sup>n, s</sup>
  127.      *  and the coefficient for dK₀<sup>n-1, s</sup> / d&Chi; when computing dK₀<sup>n, s</sup> / d&Chi;
  128.      *  </p>
  129.      *
  130.      *  <p>
  131.      *  See Danielson 2.7.3-(7c) and Collins 4-254 and 4-257
  132.      *  </p>
  133.      *
  134.      * @param n n value
  135.      * @return the polynomial
  136.      */
  137.     private PolynomialFunction a(final int n) {
  138.         // from recurrence Danielson 2.7.3-(7c), Collins 4-254
  139.         final double r1 = 2 * n + 1;
  140.         final double r2 = n + 1;

  142.         return new PolynomialFunction(new double[] {
  143.             r1 / r2
  144.         });
  145.     }

  147.     /**
  148.      * Compute polynomial coefficient b.
  149.      *
  150.      *  <p>
  151.      *  It is used to generate the coefficient for K₀<sup>n-2, s</sup> when computing K₀<sup>n, s</sup>
  152.      *  and the coefficient for dK₀<sup>n-2, s</sup> / d&Chi; when computing dK₀<sup>n, s</sup> / d&Chi;
  153.      *  </p>
  154.      *
  155.      *  <p>
  156.      *  See Danielson 2.7.3-(7c) and Collins 4-254 and 4-257
  157.      *  </p>
  158.      *
  159.      * @param n n value
  160.      * @return the polynomial
  161.      */
  162.     private PolynomialFunction b(final int n) {
  163.         // from recurrence Danielson 2.7.3-(7c), Collins 4-254
  164.         final double r1 = (n + s) * (n - s);
  165.         final double r2 = n * (n + 1);

  167.         return new PolynomialFunction(new double[] {
  168.             0.0, 0.0, -r1 / r2
  169.         });
  170.     }

  172.     /**
  173.      * Compute polynomial coefficient d.
  174.      *
  175.      *  <p>
  176.      *  It is used to generate the coefficient for K₀<sup>n-2, s</sup> when computing dK₀<sup>n, s</sup> / d&Chi;
  177.      *  </p>
  178.      *
  179.      *  <p>
  180.      *  See Danielson 2.7.3-(7c) and Collins 4-254 and 4-257
  181.      *  </p>
  182.      *
  183.      * @param n n value
  184.      * @return the polynomial
  185.      */
  186.     private PolynomialFunction d(final int n) {
  187.         // from Danielson 3.2-(3b)
  188.         final double r1 = 2 * (n + s) * (n - s);
  189.         final double r2 = n * (n + 1);

  191.         return new PolynomialFunction(new double[] {
  192.             0.0, 0.0, 0.0, r1 / r2
  193.         });
  194.     }

  196.     /**
  197.      * Generate the polynomials needed in the linear transformation.
  198.      *
  199.      * <p>
  200.      * See Petre's paper
  201.      * </p>
  202.      */
  203.     private void generatePolynomials() {

  205.         int sliceCounter = 0;

  207.         // Initialization of the matrices for linear transformations
  208.         // The final configuration of these matrices are obtained by composition
  209.         // of linear transformations

  211.         // the matrix A for the polynomials associated
  212.         // to Hansen coefficients, Petre's pater
  213.         PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix2();

  215.         // the matrix D for the polynomials associated
  216.         // to derivatives, Petre's paper
  217.         final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix2();
  218.         PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix2();
  219.         PolynomialFunctionMatrix E = HansenUtilities.buildIdentityMatrix2();

  221.         // The matrix that contains the coefficients at each step
  222.         final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix2();
  223.         a.setElem(0, 1, HansenUtilities.ONE);

  225.         // The generation process
  226.         for (int i = N0 + 2; i <= nMax; i++) {
  227.             // Collins 4-254 or Danielson 2.7.3-(7)
  228.             // Petre's paper
  229.             // The matrix of the current linear transformation is actualized
  230.             a.setMatrixLine(1, new PolynomialFunction[] {
  231.                 b(i), a(i)
  232.             });

  234.             // composition of the linear transformations to calculate
  235.             // the polynomials associated to Hansen coefficients
  236.             A = A.multiply(a);
  237.             // store the polynomials associated to Hansen coefficients
  238.             this.mpvec[i] = A.getMatrixLine(1);
  239.             // composition of the linear transformations to calculate
  240.             // the polynomials associated to derivatives
  241.             // Danielson 3.2-(3b) and Petre's paper
  242.             D = D.multiply(a);
  243.             if (sliceCounter % SLICE != 0) {
  244.                 a.setMatrixLine(1, new PolynomialFunction[] {
  245.                     b(i - 1), a(i - 1)
  246.                 });
  247.                 E = E.multiply(a);
  248.             }

  250.             B.setElem(1, 0, d(i));
  251.             // F = E.prod(B);
  252.             D = D.add(E.multiply(B));
  253.             // store the polynomials associated to the derivatives
  254.             this.mpvecDeriv[i] = D.getMatrixLine(1);

  256.             if (++sliceCounter % SLICE == 0) {
  257.                 // Re-Initialization of the matrices for linear transformations
  258.                 // The final configuration of these matrices are obtained by composition
  259.                 // of linear transformations
  260.                 A = HansenUtilities.buildIdentityMatrix2();
  261.                 D = HansenUtilities.buildZeroMatrix2();
  262.                 E = HansenUtilities.buildIdentityMatrix2();
  263.             }
  264.         }
  265.     }

  267.     /**
  268.      * Compute the initial values (see Collins, 4-255, 4-256 and 4.259)
  269.      * <p>
  270.      * K₀<sup>s, s</sup> = (-1)<sup>s</sup> * ( (2*s+1)!! / (s+1)! )
  271.      * </p>
  272.      * <p>
  273.      * K₀<sup>s+1, s</sup> = (-1)<sup>s</sup> * ( (2*s+1)!! / (s+2)!
  274.      * ) * (2*s+3 - χ<sup>-2</sup>)
  275.      * </p>
  276.      * <p>
  277.      * dK₀<sup>s+1, s</sup> / dχ = = (-1)<sup>s</sup> * 2 * (
  278.      * (2*s+1)!! / (s+2)! ) * χ<sup>-3</sup>
  279.      * </p>
  280.      * @param chitm1 sqrt(1 - e²)
  281.      * @param chitm2 sqrt(1 - e²)²
  282.      * @param chitm3 sqrt(1 - e²)³
  283.      */
  284.     public void computeInitValues(final T chitm1, final T chitm2, final T chitm3) {
  285.         final Field<T> field = chitm2.getField();
  286.         final T zero = field.getZero();
  287.         this.hansenRoot[0][0] = zero.add(twosp1dfosp1f);
  288.         this.hansenRoot[0][1] = (chitm2.negate().add(this.twosp3)).multiply(this.twosp1dfosp2f);
  289.         this.hansenDerivRoot[0][0] = zero;
  290.         this.hansenDerivRoot[0][1] = chitm3.multiply(two2sp1dfosp2f);

  292.         for (int i = 1; i < numSlices; i++) {
  293.             for (int j = 0; j < 2; j++) {
  294.                 // Get the required polynomials
  295.                 final PolynomialFunction[] mv = mpvec[s + (i * SLICE) + j];
  296.                 final PolynomialFunction[] sv = mpvecDeriv[s + (i * SLICE) + j];

  298.                 //Compute the root derivatives
  299.                 hansenDerivRoot[i][j] = mv[1].value(chitm1).multiply(hansenDerivRoot[i - 1][1]).
  300.                                     add(mv[0].value(chitm1).multiply(hansenDerivRoot[i - 1][0])).
  301.                                     add(sv[1].value(chitm1).multiply(hansenRoot[i - 1][1])).
  302.                                     add(sv[0].value(chitm1).multiply(hansenRoot[i - 1][0]));

  304.                 //Compute the root Hansen coefficients
  305.                 hansenRoot[i][j] =  mv[1].value(chitm1).multiply(hansenRoot[i - 1][1]).
  306.                                 add(mv[0].value(chitm1).multiply(hansenRoot[i - 1][0]));
  307.             }
  308.         }
  309.     }

  311.     /**
  312.      * Compute the value of the Hansen coefficient K₀<sup>n, s</sup>.
  313.      *
  314.      * @param n n value
  315.      * @param chitm1 χ<sup>-1</sup>
  316.      * @return the coefficient K₀<sup>n, s</sup>
  317.      */
  318.     public T getValue(final int n, final T chitm1) {
  319.         //Compute the potential slice
  320.         int sliceNo = (n - s) / SLICE;
  321.         if (sliceNo < numSlices) {
  322.             //Compute the index within the slice
  323.             final int indexInSlice = (n - s) % SLICE;

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

  335.         // Danielson 2.7.3-(6c)/Collins 4-242 and Petre's paper
  336.         final PolynomialFunction[] v = mpvec[n];
  337.         T ret = v[1].value(chitm1).multiply(hansenRoot[sliceNo][1]);
  338.         if (hansenRoot[sliceNo][0].getReal() != 0) {
  339.             ret = ret.add(v[0].value(chitm1).multiply(hansenRoot[sliceNo][0]));
  340.         }

  342.         return ret;

  344.     }

  346.     /**
  347.      * Compute the value of the Hansen coefficient dK₀<sup>n, s</sup> / d&Chi;.
  348.      *
  349.      * @param n n value
  350.      * @param chitm1 χ<sup>-1</sup>
  351.      * @return the coefficient dK₀<sup>n, s</sup> / d&Chi;
  352.      */
  353.     public T getDerivative(final int n, final T chitm1) {
  354.         //Compute the potential slice
  355.         int sliceNo = (n - s) / SLICE;
  356.         if (sliceNo < numSlices) {
  357.             //Compute the index within the slice
  358.             final int indexInSlice = (n - s) % SLICE;

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

  370.         final PolynomialFunction[] v = mpvec[n];
  371.         T ret = v[1].value(chitm1).multiply(hansenDerivRoot[sliceNo][1]);
  372.         if (hansenDerivRoot[sliceNo][0].getReal() != 0) {
  373.             ret = ret.add(v[0].value(chitm1).multiply(hansenDerivRoot[sliceNo][0]));
  374.         }

  376.         // Danielson 2.7.3-(7c)/Collins 4-254 and Petre's paper
  377.         final PolynomialFunction[] v1 = mpvecDeriv[n];
  378.         ret = ret.add(v1[1].value(chitm1).multiply(hansenRoot[sliceNo][1]));
  379.         if (hansenRoot[sliceNo][0].getReal() != 0) {
  380.             ret = ret.add(v1[0].value(chitm1).multiply(hansenRoot[sliceNo][0]));
  381.         }
  382.         return ret;

  384.     }

  386. }
Created with JaCoCo 0.8.8.202204050719