FieldHansenThirdBodyLinear.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-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.  * @param <T> type of the field elements
  37.  */
  38. public class FieldHansenThirdBodyLinear <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 number of slices needed to compute the coefficients. */
  59.     private int numSlices;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  121.     }

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

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

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

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

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

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

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

  206.         int sliceCounter = 0;

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

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

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

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

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

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

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

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

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

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

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

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

  312.     /**
  313.      * Compute the value of the Hansen coefficient K₀<sup>n, s</sup>.
  314.      *
  315.      * @param n n value
  316.      * @param chitm1 χ<sup>-1</sup>
  317.      * @return the coefficient K₀<sup>n, s</sup>
  318.      */
  319.     public T getValue(final int n, final T chitm1) {
  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 hansenRoot[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 2.7.3-(6c)/Collins 4-242 and Petre's paper
  337.         final PolynomialFunction[] v = mpvec[n];
  338.         T ret = v[1].value(chitm1).multiply(hansenRoot[sliceNo][1]);
  339.         if (hansenRoot[sliceNo][0].getReal() != 0) {
  340.             ret = ret.add(v[0].value(chitm1).multiply(hansenRoot[sliceNo][0]));
  341.         }

  343.         return ret;

  345.     }

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

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

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

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

  385.     }

  387. }
Created with JaCoCo 0.8.11.202310140853