FieldHansenZonalLinear.java

  1. /* Copyright 2002-2020 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 org.hipparchus.Field;
  19. import org.hipparchus.RealFieldElement;
  20. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  21. import org.hipparchus.util.FastMath;
  22. import org.hipparchus.util.MathArrays;

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

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

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

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

  45.     /** The Hansen coefficients used as roots. */
  46.     private T[][] hansenRoot;

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

  49.     /** The minimum value for the order. */
  50.     private int Nmin;


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

  53.     /** The s coefficient. */
  54.     private int s;

  55.     /**
  56.      * The offset used to identify the polynomial that corresponds to a negative
  57.      * value of n in the internal array that starts at 0.
  58.      */
  59.     private int offset;

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

  62.     /** 2<sup>s</sup>. */
  63.     private double twots;

  64.     /** 2*s+1. */
  65.     private int twosp1;

  66.     /** 2*s. */
  67.     private int twos;

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

  70.     /**
  71.      * Constructor.
  72.      *
  73.      * @param nMax the maximum (absolute) value of n coefficient
  74.      * @param s s coefficient
  75.      * @param field field used by default
  76.      */
  77.     public FieldHansenZonalLinear(final int nMax, final int s, final Field<T> field) {
  78.         //Initialize fields
  79.         this.offset = nMax + 1;
  80.         this.Nmin = -nMax - 1;
  81.         N0 = -(s + 2);
  82.         this.s = s;
  83.         this.twots = FastMath.pow(2., s);
  84.         this.twos = 2 * s;
  85.         this.twosp1 = this.twos + 1;
  86.         this.twosp1otwots = (double) this.twosp1 / this.twots;

  87.         // prepare structures for stored data
  88.         final int size = nMax - s - 1;
  89.         mpvec      = new PolynomialFunction[size][];
  90.         mpvecDeriv = new PolynomialFunction[size][];

  91.         this.numSlices  = FastMath.max((int) FastMath.ceil(((double) size) / SLICE), 1);
  92.         hansenRoot      = MathArrays.buildArray(field, numSlices, 2);
  93.         hansenDerivRoot = MathArrays.buildArray(field, numSlices, 2);

  94.         // Prepare the data base of associated polynomials
  95.         generatePolynomials();

  96.     }

  97.     /**
  98.      * Compute polynomial coefficient a.
  99.      *
  100.      *  <p>
  101.      *  It is used to generate the coefficient for K₀<sup>-n, s</sup> when computing K₀<sup>-n-1, s</sup>
  102.      *  and the coefficient for dK₀<sup>-n, s</sup> / de² when computing dK₀<sup>-n-1, s</sup> / de²
  103.      *  </p>
  104.      *
  105.      *  <p>
  106.      *  See Danielson 2.7.3-(6) and Collins 4-242 and 4-245
  107.      *  </p>
  108.      *
  109.      * @param mnm1 -n-1 value
  110.      * @return the polynomial
  111.      */
  112.     private PolynomialFunction a(final int mnm1) {
  113.         // from recurence Collins 4-242
  114.         final double d1 = (mnm1 + 2) * (2 * mnm1 + 5);
  115.         final double d2 = (mnm1 + 2 - s) * (mnm1 + 2 + s);
  116.         return new PolynomialFunction(new double[] {
  117.             0.0, 0.0, d1 / d2
  118.         });
  119.     }

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

  143.     /**
  144.      * Generate the polynomials needed in the linear transformation.
  145.      *
  146.      * <p>
  147.      * See Petre's paper
  148.      * </p>
  149.      */
  150.     private void generatePolynomials() {

  151.         int sliceCounter = 0;
  152.         int index;

  153.         // Initialisation of matrix for linear transformmations
  154.         // The final configuration of these matrix are obtained by composition
  155.         // of linear transformations
  156.         PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix2();
  157.         PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix2();
  158.         PolynomialFunctionMatrix E = HansenUtilities.buildIdentityMatrix2();

  159.         // generation of polynomials associated to Hansen coefficients and to
  160.         // their derivatives
  161.         final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix2();
  162.         a.setElem(0, 1, HansenUtilities.ONE);

  163.         //The B matrix is constant.
  164.         final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix2();
  165.         // from Collins 4-245 and Petre's paper
  166.         B.setElem(1, 1, new PolynomialFunction(new double[] {
  167.             2.0
  168.         }));

  169.         for (int i = N0 - 1; i > Nmin - 1; i--) {
  170.             index = i + offset;
  171.             // Matrix of the current linear transformation
  172.             // Petre's paper
  173.             a.setMatrixLine(1, new PolynomialFunction[] {
  174.                 b(i), a(i)
  175.             });
  176.             // composition of linear transformations
  177.             // see Petre's paper
  178.             A = A.multiply(a);
  179.             // store the polynomials for Hansen coefficients
  180.             mpvec[index] = A.getMatrixLine(1);

  181.             D = D.multiply(a);
  182.             E = E.multiply(a);
  183.             D = D.add(E.multiply(B));

  184.             // store the polynomials for Hansen coefficients from the expressions
  185.             // of derivatives
  186.             mpvecDeriv[index] = D.getMatrixLine(1);

  187.             if (++sliceCounter % SLICE == 0) {
  188.                 // Re-Initialisation of matrix for linear transformmations
  189.                 // The final configuration of these matrix are obtained by composition
  190.                 // of linear transformations
  191.                 A = HansenUtilities.buildIdentityMatrix2();
  192.                 D = HansenUtilities.buildZeroMatrix2();
  193.                 E = HansenUtilities.buildIdentityMatrix2();
  194.             }

  195.         }
  196.     }

  197.     /**
  198.      * Compute the roots for the Hansen coefficients and their derivatives.
  199.      *
  200.      * @param chi 1 / sqrt(1 - e²)
  201.      */
  202.     public void computeInitValues(final T chi) {
  203.         final Field<T> field = chi.getField();
  204.         final T zero = field.getZero();
  205.         // compute the values for n=s and n=s+1
  206.         // See Danielson 2.7.3-(6a,b)
  207.         hansenRoot[0][0] = zero;
  208.         hansenRoot[0][1] = FastMath.pow(chi, this.twosp1).divide(this.twots);
  209.         hansenDerivRoot[0][0] = zero;
  210.         hansenDerivRoot[0][1] = FastMath.pow(chi, this.twos).multiply(this.twosp1otwots);

  211.         final int st = -s - 1;
  212.         for (int i = 1; i < numSlices; i++) {
  213.             for (int j = 0; j < 2; j++) {
  214.                 // Get the required polynomials
  215.                 final PolynomialFunction[] mv = mpvec[st - (i * SLICE) - j + offset];
  216.                 final PolynomialFunction[] sv = mpvecDeriv[st - (i * SLICE) - j + offset];

  217.                 //Compute the root derivatives
  218.                 hansenDerivRoot[i][j] = mv[1].value(chi).multiply(hansenDerivRoot[i - 1][1]).
  219.                                         add(mv[0].value(chi).multiply(hansenDerivRoot[i - 1][0])).
  220.                                         add((sv[1].value(chi).multiply(hansenRoot[i - 1][1]).
  221.                                         add(sv[0].value(chi).multiply(hansenRoot[i - 1][0]))
  222.                                         ).divide(chi));
  223.                 hansenRoot[i][j] =     mv[1].value(chi).multiply(hansenRoot[i - 1][1]).
  224.                                        add(mv[0].value(chi).multiply(hansenRoot[i - 1][0]));

  225.             }

  226.         }
  227.     }

  228.     /**
  229.      * Get the K₀<sup>-n-1,s</sup> coefficient value.
  230.      *
  231.      * <p> The s value is given in the class constructor
  232.      *
  233.      * @param mnm1 (-n-1) coefficient
  234.      * @param chi The value of χ
  235.      * @return K₀<sup>-n-1,s</sup>
  236.      */
  237.     public T getValue(final int mnm1, final T chi) {
  238.         //Compute n
  239.         final int n = -mnm1 - 1;

  240.         //Compute the potential slice
  241.         int sliceNo = (n - s) / SLICE;
  242.         if (sliceNo < numSlices) {
  243.             //Compute the index within the slice
  244.             final int indexInSlice = (n - s) % SLICE;

  245.             //Check if a root must be returned
  246.             if (indexInSlice <= 1) {
  247.                 return hansenRoot[sliceNo][indexInSlice];
  248.             }
  249.         } else {
  250.             //the value was a potential root for a slice, but that slice was not required
  251.             //Decrease the slice number
  252.             sliceNo--;
  253.         }

  254.         // Danielson 2.7.3-(6c)/Collins 4-242 and Petre's paper
  255.         final PolynomialFunction[] v = mpvec[mnm1 + offset];
  256.         T ret = v[1].value(chi).multiply(hansenRoot[sliceNo][1]);
  257.         if (hansenRoot[sliceNo][0].getReal() != 0) {
  258.             ret = ret.add(v[0].value(chi).multiply(hansenRoot[sliceNo][0]));
  259.         }
  260.         return  ret;
  261.     }

  262.     /**
  263.      * Get the dK₀<sup>-n-1,s</sup> / d&Chi; coefficient derivative.
  264.      *
  265.      * <p> The s value is given in the class constructor.
  266.      *
  267.      * @param mnm1 (-n-1) coefficient
  268.      * @param chi The value of χ
  269.      * @return dK₀<sup>-n-1,s</sup> / d&Chi;
  270.      */
  271.     public T getDerivative(final int mnm1, final T chi) {
  272.         //Compute n
  273.         final int n = -mnm1 - 1;

  274.         //Compute the potential slice
  275.         int sliceNo = (n - s) / SLICE;
  276.         if (sliceNo < numSlices) {
  277.             //Compute the index within the slice
  278.             final int indexInSlice = (n - s) % SLICE;

  279.             //Check if a root must be returned
  280.             if (indexInSlice <= 1) {
  281.                 return hansenDerivRoot[sliceNo][indexInSlice];
  282.             }
  283.         } else {
  284.             //the value was a potential root for a slice, but that slice was not required
  285.             //Decrease the slice number
  286.             sliceNo--;
  287.         }

  288.         // Danielson 3.1-(7c) and Petre's paper
  289.         final PolynomialFunction[] v = mpvec[mnm1 + offset];
  290.         T ret = v[1].value(chi).multiply(hansenDerivRoot[sliceNo][1]);
  291.         if (hansenDerivRoot[sliceNo][0].getReal() != 0) {
  292.             ret = ret.add(v[0].value(chi).multiply(hansenDerivRoot[sliceNo][0]));
  293.         }

  294.         // Danielson 2.7.3-(6b)
  295.         final PolynomialFunction[] v1 = mpvecDeriv[mnm1 + offset];
  296.         T hret = v1[1].value(chi).multiply(hansenRoot[sliceNo][1]);
  297.         if (hansenRoot[sliceNo][0].getReal() != 0) {
  298.             hret = hret.add(v1[0].value(chi).multiply(hansenRoot[sliceNo][0]));
  299.         }
  300.         ret = ret.add(hret.divide(chi));

  301.         return ret;

  302.     }

  303. }