HansenThirdBodyLinear.java

  1. /* Copyright 2002-2016 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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.apache.commons.math3.analysis.polynomials.PolynomialFunction;
  19. import org.apache.commons.math3.util.FastMath;

  20. /**
  21.  * Hansen coefficients K(t,n,s) for t=0 and n > 0.
  22.  * <p>
  23.  * Implements Collins 4-254 or Danielson 2.7.3-(7) for Hansen Coefficients and
  24.  * Danielson 3.2-(3) 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 HansenThirdBodyLinear {

  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 number of slices needed to compute the coefficients. */
  46.     private int numSlices;

  47.     /** The maximum order of n indexes. */
  48.     private int nMax;

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

  51.     /** The s index. */
  52.     private int s;

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

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

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

  59.     /** (2*s + 3). */
  60.     private double twosp3;

  61.     /**
  62.      * Constructor.
  63.      *
  64.      * @param nMax the maximum value of n
  65.      * @param s the value of s
  66.      */
  67.     public HansenThirdBodyLinear(final int nMax, final int s) {

  68.         // initialise fields
  69.         this.nMax = nMax;
  70.         N0 = s;
  71.         this.s = s;

  72.         // initialization of structures for stored data
  73.         mpvec = new PolynomialFunction[this.nMax + 1][];
  74.         mpvecDeriv = new PolynomialFunction[this.nMax + 1][];

  75.         //Compute the fields that will be used to determine the initial values for the coefficients
  76.         this.twosp1dfosp1f = (s % 2 == 0) ? 1.0 : -1.0;
  77.         for (int i = s; i >= 1; i--) {
  78.             this.twosp1dfosp1f *= (2.0 * i + 1.0) / (i + 1.0);
  79.         }

  80.         this.twosp1dfosp2f = this.twosp1dfosp1f / (s + 2.);
  81.         this.twosp3 = 2 * s + 3;
  82.         this.two2sp1dfosp2f = 2 * this.twosp1dfosp2f;

  83.         // initialization of structures for stored data
  84.         mpvec = new PolynomialFunction[this.nMax + 1][];
  85.         mpvecDeriv = new PolynomialFunction[this.nMax + 1][];

  86.         this.numSlices = (int) FastMath.ceil(((double) nMax - s - 1) / SLICE);
  87.         hansenRoot = new double[numSlices][2];
  88.         hansenDerivRoot = new double[numSlices][2];

  89.         // Prepare the database of the 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-1, s</sup> when computing K₀<sup>n, s</sup>
  97.      *  and the coefficient for dK₀<sup>n-1, s</sup> / d&Chi; when computing dK₀<sup>n, s</sup> / d&Chi;
  98.      *  </p>
  99.      *
  100.      *  <p>
  101.      *  See Danielson 2.7.3-(7c) and Collins 4-254 and 4-257
  102.      *  </p>
  103.      *
  104.      * @param n n value
  105.      * @return the polynomial
  106.      */
  107.     private PolynomialFunction a(final int n) {
  108.         // from recurrence Danielson 2.7.3-(7c), Collins 4-254
  109.         final double r1 = 2 * n + 1;
  110.         final double r2 = n + 1;

  111.         return new PolynomialFunction(new double[] {
  112.             r1 / r2
  113.         });
  114.     }

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

  134.         return new PolynomialFunction(new double[] {
  135.             0.0, 0.0, -r1 / r2
  136.         });
  137.     }

  138.     /**
  139.      * Compute polynomial coefficient d.
  140.      *
  141.      *  <p>
  142.      *  It is used to generate the coefficient for K₀<sup>n-2, s</sup> when computing dK₀<sup>n, s</sup> / d&Chi;
  143.      *  </p>
  144.      *
  145.      *  <p>
  146.      *  See Danielson 2.7.3-(7c) and Collins 4-254 and 4-257
  147.      *  </p>
  148.      *
  149.      * @param n n value
  150.      * @return the polynomial
  151.      */
  152.     private PolynomialFunction d(final int n) {
  153.         // from Danielson 3.2-(3b)
  154.         final double r1 = 2 * (n + s) * (n - s);
  155.         final double r2 = n * (n + 1);

  156.         return new PolynomialFunction(new double[] {
  157.             0.0, 0.0, 0.0, r1 / r2
  158.         });
  159.     }

  160.     /**
  161.      * Generate the polynomials needed in the linear transformation.
  162.      *
  163.      * <p>
  164.      * See Petre's paper
  165.      * </p>
  166.      */
  167.     private void generatePolynomials() {

  168.         int sliceCounter = 0;

  169.         // Initialization of the matrices for linear transformations
  170.         // The final configuration of these matrices are obtained by composition
  171.         // of linear transformations

  172.         // the matrix A for the polynomials associated
  173.         // to Hansen coefficients, Petre's pater
  174.         PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix2();

  175.         // the matrix D for the polynomials associated
  176.         // to derivatives, Petre's paper
  177.         final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix2();
  178.         PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix2();
  179.         PolynomialFunctionMatrix E = HansenUtilities.buildIdentityMatrix2();

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

  183.         // The generation process
  184.         for (int i = N0 + 2; i <= nMax; i++) {
  185.             // Collins 4-254 or Danielson 2.7.3-(7)
  186.             // Petre's paper
  187.             // The matrix of the current linear transformation is actualized
  188.             a.setMatrixLine(1, new PolynomialFunction[] {
  189.                 b(i), a(i)
  190.             });

  191.             // composition of the linear transformations to calculate
  192.             // the polynomials associated to Hansen coefficients
  193.             A = A.multiply(a);
  194.             // store the polynomials associated to Hansen coefficients
  195.             this.mpvec[i] = A.getMatrixLine(1);
  196.             // composition of the linear transformations to calculate
  197.             // the polynomials associated to derivatives
  198.             // Danielson 3.2-(3b) and Petre's paper
  199.             D = D.multiply(a);
  200.             if (sliceCounter % SLICE != 0) {
  201.                 a.setMatrixLine(1, new PolynomialFunction[] {
  202.                     b(i - 1), a(i - 1)
  203.                 });
  204.                 E = E.multiply(a);
  205.             }

  206.             B.setElem(1, 0, d(i));
  207.             // F = E.prod(B);
  208.             D = D.add(E.multiply(B));
  209.             // store the polynomials associated to the derivatives
  210.             this.mpvecDeriv[i] = D.getMatrixLine(1);

  211.             if (++sliceCounter % SLICE == 0) {
  212.                 // Re-Initialization of the matrices for linear transformations
  213.                 // The final configuration of these matrices are obtained by composition
  214.                 // of linear transformations
  215.                 A = HansenUtilities.buildIdentityMatrix2();
  216.                 D = HansenUtilities.buildZeroMatrix2();
  217.                 E = HansenUtilities.buildIdentityMatrix2();
  218.             }
  219.         }
  220.     }

  221.     /**
  222.      * Compute the initial values (see Collins, 4-255, 4-256 and 4.259)
  223.      * <p>
  224.      * K₀<sup>s, s</sup> = (-1)<sup>s</sup> * ( (2*s+1)!! / (s+1)! )
  225.      * </p>
  226.      * <p>
  227.      * K₀<sup>s+1, s</sup> = (-1)<sup>s</sup> * ( (2*s+1)!! / (s+2)!
  228.      * ) * (2*s+3 - χ<sup>-2</sup>)
  229.      * </p>
  230.      * <p>
  231.      * dK₀<sup>s+1, s</sup> / dχ = = (-1)<sup>s</sup> * 2 * (
  232.      * (2*s+1)!! / (s+2)! ) * χ<sup>-3</sup>
  233.      * </p>
  234.      * @param chitm1 sqrt(1 - e²)
  235.      * @param chitm2 sqrt(1 - e²)²
  236.      * @param chitm3 sqrt(1 - e²)³
  237.      */
  238.     public void computeInitValues(final double chitm1, final double chitm2, final double chitm3) {
  239.         this.hansenRoot[0][0] = this.twosp1dfosp1f;
  240.         this.hansenRoot[0][1] = this.twosp1dfosp2f * (this.twosp3 - chitm2);
  241.         this.hansenDerivRoot[0][0] = 0;
  242.         this.hansenDerivRoot[0][1] = this.two2sp1dfosp2f * chitm3;

  243.         for (int i = 1; i < numSlices; i++) {
  244.             for (int j = 0; j < 2; j++) {
  245.                 // Get the required polynomials
  246.                 final PolynomialFunction[] mv = mpvec[s + (i * SLICE) + j];
  247.                 final PolynomialFunction[] sv = mpvecDeriv[s + (i * SLICE) + j];

  248.                 //Compute the root derivatives
  249.                 hansenDerivRoot[i][j] = mv[1].value(chitm1) * hansenDerivRoot[i - 1][1] +
  250.                                         mv[0].value(chitm1) * hansenDerivRoot[i - 1][0] +
  251.                                         sv[1].value(chitm1) * hansenRoot[i - 1][1] +
  252.                                         sv[0].value(chitm1) * hansenRoot[i - 1][0];

  253.                 //Compute the root Hansen coefficients
  254.                 hansenRoot[i][j] =  mv[1].value(chitm1) * hansenRoot[i - 1][1] +
  255.                                     mv[0].value(chitm1) * hansenRoot[i - 1][0];
  256.             }
  257.         }
  258.     }

  259.     /**
  260.      * Compute the value of the Hansen coefficient K₀<sup>n, s</sup>.
  261.      *
  262.      * @param n n value
  263.      * @param chitm1 χ<sup>-1</sup>
  264.      * @return the coefficient K₀<sup>n, s</sup>
  265.      */
  266.     public double getValue(final int n, final double chitm1) {

  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 hansenRoot[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 2.7.3-(6c)/Collins 4-242 and Petre's paper
  282.         final PolynomialFunction[] v = mpvec[n];
  283.         double ret = v[1].value(chitm1) * hansenRoot[sliceNo][1];
  284.         if (hansenRoot[sliceNo][0] != 0) {
  285.             ret += v[0].value(chitm1) * hansenRoot[sliceNo][0];
  286.         }

  287.         return ret;

  288.     }

  289.     /**
  290.      * Compute the value of the Hansen coefficient dK₀<sup>n, s</sup> / d&Chi;.
  291.      *
  292.      * @param n n value
  293.      * @param chitm1 χ<sup>-1</sup>
  294.      * @return the coefficient dK₀<sup>n, s</sup> / d&Chi;
  295.      */
  296.     public double getDerivative(final int n, final double chitm1) {
  297.         //Compute the potential slice
  298.         int sliceNo = (n - s) / SLICE;
  299.         if (sliceNo < numSlices) {
  300.             //Compute the index within the slice
  301.             final int indexInSlice = (n - s) % SLICE;

  302.             //Check if a root must be returned
  303.             if (indexInSlice <= 1) {
  304.                 return hansenDerivRoot[sliceNo][indexInSlice];
  305.             }
  306.         } else {
  307.             //the value was a potential root for a slice, but that slice was not required
  308.             //Decrease the slice number
  309.             sliceNo--;
  310.         }

  311.         final PolynomialFunction[] v = mpvec[n];
  312.         double ret = v[1].value(chitm1) * hansenDerivRoot[sliceNo][1];
  313.         if (hansenDerivRoot[sliceNo][0] != 0) {
  314.             ret += v[0].value(chitm1) * hansenDerivRoot[sliceNo][0];
  315.         }

  316.         // Danielson 2.7.3-(7c)/Collins 4-254 and Petre's paper
  317.         final PolynomialFunction[] v1 = mpvecDeriv[n];
  318.         ret += v1[1].value(chitm1) * hansenRoot[sliceNo][1];
  319.         if (hansenRoot[sliceNo][0] != 0) {
  320.             ret += v1[0].value(chitm1) * hansenRoot[sliceNo][0];
  321.         }
  322.         return ret;

  323.     }

  324. }