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

  18. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  19. import org.hipparchus.util.FastMath;

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

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

  34.     /**
  35.      * The first vector of polynomials associated to Hansen coefficients and
  36.      * derivatives.
  37.      */
  38.     private PolynomialFunction[][] mpvec;

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

  41.     /** The Hansen coefficients used as roots. */
  42.     private double[][] hansenRoot;

  43.     /** The derivatives of the Hansen coefficients used as roots. */
  44.     private double[][] hansenDerivRoot;

  45.     /** The minimum value for the order. */
  46.     private int Nmin;


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

  49.     /** The s coefficient. */
  50.     private int s;

  51.     /**
  52.      * The offset used to identify the polynomial that corresponds to a negative
  53.      * value of n in the internal array that starts at 0.
  54.      */
  55.     private int offset;

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

  58.     /** 2<sup>s</sup>. */
  59.     private double twots;

  60.     /** 2*s+1. */
  61.     private int twosp1;

  62.     /** 2*s. */
  63.     private int twos;

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

  66.     /**
  67.      * Constructor.
  68.      *
  69.      * @param nMax the maximum (absolute) value of n coefficient
  70.      * @param s s coefficient
  71.      */
  72.     public HansenZonalLinear(final int nMax, final int s) {

  73.         //Initialize fields
  74.         this.offset = nMax + 1;
  75.         this.Nmin = -nMax - 1;
  76.         N0 = -(s + 2);
  77.         this.s = s;
  78.         this.twots = FastMath.pow(2., s);
  79.         this.twos = 2 * s;
  80.         this.twosp1 = this.twos + 1;
  81.         this.twosp1otwots = (double) this.twosp1 / this.twots;

  82.         // prepare structures for stored data
  83.         final int size = nMax - s - 1;
  84.         mpvec      = new PolynomialFunction[size][];
  85.         mpvecDeriv = new PolynomialFunction[size][];

  86.         this.numSlices  = FastMath.max((int) FastMath.ceil(((double) size) / SLICE), 1);
  87.         hansenRoot      = new double[numSlices][2];
  88.         hansenDerivRoot = new double[numSlices][2];

  89.         // Prepare the data base of associated polynomials
  90.         generatePolynomials();

  91.     }

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

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

  138.     /**
  139.      * Generate the polynomials needed in the linear transformation.
  140.      *
  141.      * <p>
  142.      * See Petre's paper
  143.      * </p>
  144.      */
  145.     private void generatePolynomials() {

  146.         int sliceCounter = 0;
  147.         int index;

  148.         // Initialisation of matrix for linear transformmations
  149.         // The final configuration of these matrix are obtained by composition
  150.         // of linear transformations
  151.         PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix2();
  152.         PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix2();
  153.         PolynomialFunctionMatrix E = HansenUtilities.buildIdentityMatrix2();

  154.         // generation of polynomials associated to Hansen coefficients and to
  155.         // their derivatives
  156.         final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix2();
  157.         a.setElem(0, 1, HansenUtilities.ONE);

  158.         //The B matrix is constant.
  159.         final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix2();
  160.         // from Collins 4-245 and Petre's paper
  161.         B.setElem(1, 1, new PolynomialFunction(new double[] {
  162.             2.0
  163.         }));

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

  176.             D = D.multiply(a);
  177.             E = E.multiply(a);
  178.             D = D.add(E.multiply(B));

  179.             // store the polynomials for Hansen coefficients from the expressions
  180.             // of derivatives
  181.             mpvecDeriv[index] = D.getMatrixLine(1);

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

  190.         }
  191.     }

  192.     /**
  193.      * Compute the roots for the Hansen coefficients and their derivatives.
  194.      *
  195.      * @param chi 1 / sqrt(1 - e²)
  196.      */
  197.     public void computeInitValues(final double chi) {
  198.         // compute the values for n=s and n=s+1
  199.         // See Danielson 2.7.3-(6a,b)
  200.         hansenRoot[0][0] = 0;
  201.         hansenRoot[0][1] = FastMath.pow(chi, this.twosp1) / this.twots;
  202.         hansenDerivRoot[0][0] = 0;
  203.         hansenDerivRoot[0][1] = this.twosp1otwots * FastMath.pow(chi, this.twos);

  204.         final int st = -s - 1;
  205.         for (int i = 1; i < numSlices; i++) {
  206.             for (int j = 0; j < 2; j++) {
  207.                 // Get the required polynomials
  208.                 final PolynomialFunction[] mv = mpvec[st - (i * SLICE) - j + offset];
  209.                 final PolynomialFunction[] sv = mpvecDeriv[st - (i * SLICE) - j + offset];

  210.                 //Compute the root derivatives
  211.                 hansenDerivRoot[i][j] = mv[1].value(chi) * hansenDerivRoot[i - 1][1] +
  212.                                         mv[0].value(chi) * hansenDerivRoot[i - 1][0] +
  213.                                         (sv[1].value(chi) * hansenRoot[i - 1][1] +
  214.                                          sv[0].value(chi) * hansenRoot[i - 1][0]
  215.                                         ) / chi;
  216.                 hansenRoot[i][j] =     mv[1].value(chi) * hansenRoot[i - 1][1] +
  217.                                        mv[0].value(chi) * hansenRoot[i - 1][0];

  218.             }

  219.         }
  220.     }

  221.     /**
  222.      * Get the K₀<sup>-n-1,s</sup> coefficient value.
  223.      *
  224.      * <p> The s value is given in the class constructor
  225.      *
  226.      * @param mnm1 (-n-1) coefficient
  227.      * @param chi The value of χ
  228.      * @return K₀<sup>-n-1,s</sup>
  229.      */
  230.     public double getValue(final int mnm1, final double chi) {

  231.         //Compute n
  232.         final int n = -mnm1 - 1;

  233.         //Compute the potential slice
  234.         int sliceNo = (n - s) / SLICE;
  235.         if (sliceNo < numSlices) {
  236.             //Compute the index within the slice
  237.             final int indexInSlice = (n - s) % SLICE;

  238.             //Check if a root must be returned
  239.             if (indexInSlice <= 1) {
  240.                 return hansenRoot[sliceNo][indexInSlice];
  241.             }
  242.         } else {
  243.             //the value was a potential root for a slice, but that slice was not required
  244.             //Decrease the slice number
  245.             sliceNo--;
  246.         }

  247.         // Danielson 2.7.3-(6c)/Collins 4-242 and Petre's paper
  248.         final PolynomialFunction[] v = mpvec[mnm1 + offset];
  249.         double ret = v[1].value(chi) * hansenRoot[sliceNo][1];
  250.         if (hansenRoot[sliceNo][0] != 0) {
  251.             ret += v[0].value(chi) * hansenRoot[sliceNo][0];
  252.         }
  253.         return  ret;
  254.     }

  255.     /**
  256.      * Get the dK₀<sup>-n-1,s</sup> / d&Chi; coefficient derivative.
  257.      *
  258.      * <p> The s value is given in the class constructor.
  259.      *
  260.      * @param mnm1 (-n-1) coefficient
  261.      * @param chi The value of χ
  262.      * @return dK₀<sup>-n-1,s</sup> / d&Chi;
  263.      */
  264.     public double getDerivative(final int mnm1, final double chi) {

  265.         //Compute n
  266.         final int n = -mnm1 - 1;

  267.         //Compute the potential slice
  268.         int sliceNo = (n - s) / SLICE;
  269.         if (sliceNo < numSlices) {
  270.             //Compute the index within the slice
  271.             final int indexInSlice = (n - s) % SLICE;

  272.             //Check if a root must be returned
  273.             if (indexInSlice <= 1) {
  274.                 return hansenDerivRoot[sliceNo][indexInSlice];
  275.             }
  276.         } else {
  277.             //the value was a potential root for a slice, but that slice was not required
  278.             //Decrease the slice number
  279.             sliceNo--;
  280.         }

  281.         // Danielson 3.1-(7c) and Petre's paper
  282.         final PolynomialFunction[] v = mpvec[mnm1 + offset];
  283.         double ret = v[1].value(chi) * hansenDerivRoot[sliceNo][1];
  284.         if (hansenDerivRoot[sliceNo][0] != 0) {
  285.             ret += v[0].value(chi) * hansenDerivRoot[sliceNo][0];
  286.         }

  287.         // Danielson 2.7.3-(6b)
  288.         final PolynomialFunction[] v1 = mpvecDeriv[mnm1 + offset];
  289.         double hret = v1[1].value(chi) * hansenRoot[sliceNo][1];
  290.         if (hansenRoot[sliceNo][0] != 0) {
  291.             hret += v1[0].value(chi) * hansenRoot[sliceNo][0];
  292.         }
  293.         ret += hret / chi;

  294.         return ret;

  295.     }

  296. }