HansenTesseralLinear.java

  1. /* Copyright 2002-2018 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;

  19. import org.hipparchus.analysis.differentiation.DSFactory;
  20. import org.hipparchus.analysis.differentiation.DerivativeStructure;
  21. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  22. import org.hipparchus.util.FastMath;
  23. import org.orekit.propagation.semianalytical.dsst.utilities.NewcombOperators;

  25. /**
  26.  * Hansen coefficients K(t,n,s) for t!=0 and n < 0.
  27.  * <p>
  28.  * Implements Collins 4-236 or Danielson 2.7.3-(9) for Hansen Coefficients and
  29.  * Collins 4-240 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.  */
  36. public class HansenTesseralLinear {

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

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

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

  50.     /** The Hansen coefficients used as roots. */
  51.     private double[][] hansenRoot;

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

  56.     /** The minimum value for the order. */
  57.     private int Nmin;

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

  62.     /** The s coefficient. */
  63.     private int s;

  65.     /** The j coefficient. */
  66.     private int j;

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

  71.     /**
  72.      * The offset used to identify the polynomial that corresponds to a negative.
  73.      * value of n in the internal array that starts at 0
  74.      */
  75.     private int offset;

  77.     /** The objects used to calculate initial data by means of Newcomb operators. */
  78.     private HansenCoefficientsBySeries[] hansenInit;

  80.     /**
  81.      * Constructor.
  82.      *
  83.      * @param nMax the maximum (absolute) value of n parameter
  84.      * @param s s parameter
  85.      * @param j j parameter
  86.      * @param n0 the minimum (absolute) value of n
  87.      * @param maxHansen maximum power of e2 in Hansen expansion
  88.      */
  89.     public HansenTesseralLinear(final int nMax, final int s, final int j, final int n0, final int maxHansen) {
  90.         //Initialize the fields
  91.         this.offset = nMax + 1;
  92.         this.Nmin = -nMax - 1;
  93.         this.N0 = -n0 - 4;
  94.         this.s = s;
  95.         this.j = j;

  97.         //Ensure that only the needed terms are computed
  98.         final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
  99.         this.hansenInit = new HansenCoefficientsBySeries[maxRoots];
  100.         for (int i = 0; i < maxRoots; i++) {
  101.             this.hansenInit[i] = new HansenCoefficientsBySeries(N0 - i + 3, s, j, maxHansen);
  102.         }

  104.         // The first 4 values are computed with series. No linear combination is needed.
  105.         final int size = N0 - Nmin;
  106.         this.numSlices = (int) FastMath.max(FastMath.ceil(((double) size) / SLICE), 1);
  107.         hansenRoot = new double[numSlices][4];
  108.         hansenDerivRoot = new double[numSlices][4];
  109.         if (size > 0) {
  110.             mpvec = new PolynomialFunction[size][];
  111.             mpvecDeriv = new PolynomialFunction[size][];

  113.             // Prepare the database of the associated polynomials
  114.             generatePolynomials();
  115.         }

  117.     }

  119.     /**
  120.      * Compute polynomial coefficient a.
  121.      *
  122.      *  <p>
  123.      *  It is used to generate the coefficient for K<sub>j</sub><sup>-n, s</sup> when computing K<sub>j</sub><sup>-n-1, s</sup>
  124.      *  and the coefficient for dK<sub>j</sub><sup>-n, s</sup> / de² when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
  125.      *  </p>
  126.      *
  127.      *  <p>
  128.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  129.      *  </p>
  130.      *
  131.      * @param mnm1 -n-1
  132.      * @return the polynomial
  133.      */
  134.     private PolynomialFunction a(final int mnm1) {
  135.         // Collins 4-236, Danielson 2.7.3-(9)
  136.         final double r1 = (mnm1 + 2.) * (2. * mnm1 + 5.);
  137.         final double r2 = (2. + mnm1 + s) * (2. + mnm1 - s);
  138.         return new PolynomialFunction(new double[] {
  139.             0.0, 0.0, r1 / r2
  140.         });
  141.     }

  143.     /**
  144.      * Compute polynomial coefficient b.
  145.      *
  146.      *  <p>
  147.      *  It is used to generate the coefficient for K<sub>j</sub><sup>-n+1, s</sup> when computing K<sub>j</sub><sup>-n-1, s</sup>
  148.      *  and the coefficient for dK<sub>j</sub><sup>-n+1, s</sup> / de² when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
  149.      *  </p>
  150.      *
  151.      *  <p>
  152.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  153.      *  </p>
  154.      *
  155.      * @param mnm1 -n-1
  156.      * @return the polynomial
  157.      */
  158.     private PolynomialFunction b(final int mnm1) {
  159.         // Collins 4-236, Danielson 2.7.3-(9)
  160.         final double r2 = (2. + mnm1 + s) * (2. + mnm1 - s);
  161.         final double d1 = (mnm1 + 3.) * 2. * j * s / (r2 * (mnm1 + 4.));
  162.         final double d2 = (mnm1 + 3.) * (mnm1 + 2.) / r2;
  163.         return new PolynomialFunction(new double[] {
  164.             0.0, -d1, -d2
  165.         });
  166.     }

  168.     /**
  169.      * Compute polynomial coefficient c.
  170.      *
  171.      *  <p>
  172.      *  It is used to generate the coefficient for K<sub>j</sub><sup>-n+3, s</sup> when computing K<sub>j</sub><sup>-n-1, s</sup>
  173.      *  and the coefficient for dK<sub>j</sub><sup>-n+3, s</sup> / de² when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
  174.      *  </p>
  175.      *
  176.      *  <p>
  177.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  178.      *  </p>
  179.      *
  180.      * @param mnm1 -n-1
  181.      * @return the polynomial
  182.      */
  183.     private PolynomialFunction c(final int mnm1) {
  184.         // Collins 4-236, Danielson 2.7.3-(9)
  185.         final double r1 = j * j * (mnm1 + 2.);
  186.         final double r2 = (mnm1 + 4.) * (2. + mnm1 + s) * (2. + mnm1 - s);

  188.         return new PolynomialFunction(new double[] {
  189.             0.0, 0.0, r1 / r2
  190.         });
  191.     }

  193.     /**
  194.      * Compute polynomial coefficient d.
  195.      *
  196.      *  <p>
  197.      *  It is used to generate the coefficient for K<sub>j</sub><sup>-n-1, s</sup> / dχ when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
  198.      *  </p>
  199.      *
  200.      *  <p>
  201.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  202.      *  </p>
  203.      *
  204.      * @param mnm1 -n-1
  205.      * @return the polynomial
  206.      */
  207.     private PolynomialFunction d(final int mnm1) {
  208.         // Collins 4-236, Danielson 2.7.3-(9)
  209.         return new PolynomialFunction(new double[] {
  210.             0.0, 0.0, 1.0
  211.         });
  212.     }

  214.     /**
  215.      * Compute polynomial coefficient f.
  216.      *
  217.      *  <p>
  218.      *  It is used to generate the coefficient for K<sub>j</sub><sup>-n+1, s</sup> / dχ when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
  219.      *  </p>
  220.      *
  221.      *  <p>
  222.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  223.      *  </p>
  224.      *
  225.      * @param n index
  226.      * @return the polynomial
  227.      */
  228.     private PolynomialFunction f(final int n) {
  229.         // Collins 4-236, Danielson 2.7.3-(9)
  230.         final double r1 = (n + 3.0) * j * s;
  231.         final double r2 = (n + 4.0) * (2.0 + n + s) * (2.0 + n - s);
  232.         return new PolynomialFunction(new double[] {
  233.             0.0, 0.0, 0.0, r1 / r2
  234.         });
  235.     }

  237.     /**
  238.      * Generate the polynomials needed in the linear transformation.
  239.      *
  240.      * <p>
  241.      * See Petre's paper
  242.      * </p>
  243.      */
  244.     private void generatePolynomials() {


  247.         // Initialization of the matrices for linear transformations
  248.         // The final configuration of these matrices are obtained by composition
  249.         // of linear transformations

  251.         // The matrix of polynomials associated to Hansen coefficients, Petre's
  252.         // paper
  253.         PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix4();

  255.         // The matrix of polynomials associated to derivatives, Petre's paper
  256.         final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix4();
  257.         PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix4();
  258.         final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix4();

  260.         // The matrix of the current linear transformation
  261.         a.setMatrixLine(0, new PolynomialFunction[] {
  262.             HansenUtilities.ZERO, HansenUtilities.ONE, HansenUtilities.ZERO, HansenUtilities.ZERO
  263.         });
  264.         a.setMatrixLine(1, new PolynomialFunction[] {
  265.             HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ONE, HansenUtilities.ZERO
  266.         });
  267.         a.setMatrixLine(2, new PolynomialFunction[] {
  268.             HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ONE
  269.         });
  270.         // The generation process
  271.         int index;
  272.         int sliceCounter = 0;
  273.         for (int i = N0 - 1; i > Nmin - 1; i--) {
  274.             index = i + this.offset;
  275.             // The matrix of the current linear transformation is updated
  276.             // Petre's paper
  277.             a.setMatrixLine(3, new PolynomialFunction[] {
  278.                     c(i), HansenUtilities.ZERO, b(i), a(i)
  279.             });

  281.             // composition of the linear transformations to calculate
  282.             // the polynomials associated to Hansen coefficients
  283.             // Petre's paper
  284.             A = A.multiply(a);
  285.             // store the polynomials for Hansen coefficients
  286.             mpvec[index] = A.getMatrixLine(3);
  287.             // composition of the linear transformations to calculate
  288.             // the polynomials associated to derivatives
  289.             // Petre's paper
  290.             D = D.multiply(a);

  292.             //Update the B matrix
  293.             B.setMatrixLine(3, new PolynomialFunction[] {
  294.                 HansenUtilities.ZERO, f(i),
  295.                 HansenUtilities.ZERO, d(i)
  296.             });
  297.             D = D.add(A.multiply(B));

  299.             // store the polynomials for Hansen coefficients from the
  300.             // expressions of derivatives
  301.             mpvecDeriv[index] = D.getMatrixLine(3);

  303.             if (++sliceCounter % SLICE == 0) {
  304.                 // Re-Initialisation of matrix for linear transformmations
  305.                 // The final configuration of these matrix are obtained by composition
  306.                 // of linear transformations
  307.                 A = HansenUtilities.buildIdentityMatrix4();
  308.                 D = HansenUtilities.buildZeroMatrix4();
  309.             }
  310.         }
  311.     }

  313.     /**
  314.      * Compute the values for the first four coefficients and their derivatives by means of series.
  315.      *
  316.      * @param e2 e²
  317.      * @param chi &Chi;
  318.      * @param chi2 &Chi;²
  319.      */
  320.     public void computeInitValues(final double e2, final double chi, final double chi2) {
  321.         // compute the values for n, n+1, n+2 and n+3 by series
  322.         // See Danielson 2.7.3-(10)
  323.         //Ensure that only the needed terms are computed
  324.         final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
  325.         for (int i = 0; i < maxRoots; i++) {
  326.             final DerivativeStructure hansenKernel = hansenInit[i].getValue(e2, chi, chi2);
  327.             this.hansenRoot[0][i] = hansenKernel.getValue();
  328.             this.hansenDerivRoot[0][i] = hansenKernel.getPartialDerivative(1);
  329.         }

  331.         for (int i = 1; i < numSlices; i++) {
  332.             for (int k = 0; k < 4; k++) {
  333.                 final PolynomialFunction[] mv = mpvec[N0 - (i * SLICE) - k + 3 + offset];
  334.                 final PolynomialFunction[] sv = mpvecDeriv[N0 - (i * SLICE) - k + 3 + offset];

  336.                 hansenDerivRoot[i][k] = mv[3].value(chi) * hansenDerivRoot[i - 1][3] +
  337.                                         mv[2].value(chi) * hansenDerivRoot[i - 1][2] +
  338.                                         mv[1].value(chi) * hansenDerivRoot[i - 1][1] +
  339.                                         mv[0].value(chi) * hansenDerivRoot[i - 1][0] +
  340.                                         sv[3].value(chi) * hansenRoot[i - 1][3] +
  341.                                         sv[2].value(chi) * hansenRoot[i - 1][2] +
  342.                                         sv[1].value(chi) * hansenRoot[i - 1][1] +
  343.                                         sv[0].value(chi) * hansenRoot[i - 1][0];

  345.                 hansenRoot[i][k] =  mv[3].value(chi) * hansenRoot[i - 1][3] +
  346.                                     mv[2].value(chi) * hansenRoot[i - 1][2] +
  347.                                     mv[1].value(chi) * hansenRoot[i - 1][1] +
  348.                                     mv[0].value(chi) * hansenRoot[i - 1][0];
  349.             }
  350.         }
  351.     }

  353.     /**
  354.      * Compute the value of the Hansen coefficient K<sub>j</sub><sup>-n-1, s</sup>.
  355.      *
  356.      * @param mnm1 -n-1
  357.      * @param chi χ
  358.      * @return the coefficient K<sub>j</sub><sup>-n-1, s</sup>
  359.      */
  360.     public double getValue(final int mnm1, final double chi) {
  361.         //Compute n
  362.         final int n = -mnm1 - 1;

  364.         //Compute the potential slice
  365.         int sliceNo = (n + N0 + 4) / SLICE;
  366.         if (sliceNo < numSlices) {
  367.             //Compute the index within the slice
  368.             final int indexInSlice = (n + N0 + 4) % SLICE;

  370.             //Check if a root must be returned
  371.             if (indexInSlice <= 3) {
  372.                 return hansenRoot[sliceNo][indexInSlice];
  373.             }
  374.         } else {
  375.             //the value was a potential root for a slice, but that slice was not required
  376.             //Decrease the slice number
  377.             sliceNo--;
  378.         }

  380.         // Computes the coefficient by linear transformation
  381.         // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
  382.         final PolynomialFunction[] v = mpvec[mnm1 + offset];
  383.         return v[3].value(chi) * hansenRoot[sliceNo][3] +
  384.                v[2].value(chi) * hansenRoot[sliceNo][2] +
  385.                v[1].value(chi) * hansenRoot[sliceNo][1] +
  386.                v[0].value(chi) * hansenRoot[sliceNo][0];

  388.     }

  390.     /**
  391.      * Compute the value of the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de².
  392.      *
  393.      * @param mnm1 -n-1
  394.      * @param chi χ
  395.      * @return the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de²
  396.      */
  397.     public double getDerivative(final int mnm1, final double chi) {

  399.         //Compute n
  400.         final int n = -mnm1 - 1;

  402.         //Compute the potential slice
  403.         int sliceNo = (n + N0 + 4) / SLICE;
  404.         if (sliceNo < numSlices) {
  405.             //Compute the index within the slice
  406.             final int indexInSlice = (n + N0 + 4) % SLICE;

  408.             //Check if a root must be returned
  409.             if (indexInSlice <= 3) {
  410.                 return hansenDerivRoot[sliceNo][indexInSlice];
  411.             }
  412.         } else {
  413.             //the value was a potential root for a slice, but that slice was not required
  414.             //Decrease the slice number
  415.             sliceNo--;
  416.         }

  418.         // Computes the coefficient by linear transformation
  419.         // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
  420.         final PolynomialFunction[] v = mpvec[mnm1 + this.offset];
  421.         final PolynomialFunction[] vv = mpvecDeriv[mnm1 + this.offset];

  423.         return v[3].value(chi)  * hansenDerivRoot[sliceNo][3] +
  424.                v[2].value(chi)  * hansenDerivRoot[sliceNo][2] +
  425.                v[1].value(chi)  * hansenDerivRoot[sliceNo][1] +
  426.                v[0].value(chi)  * hansenDerivRoot[sliceNo][0] +
  427.                vv[3].value(chi) * hansenRoot[sliceNo][3] +
  428.                vv[2].value(chi) * hansenRoot[sliceNo][2] +
  429.                vv[1].value(chi) * hansenRoot[sliceNo][1] +
  430.                vv[0].value(chi) * hansenRoot[sliceNo][0];

  432.     }

  434.     /**
  435.      * Compute a Hansen coefficient with series.
  436.      * <p>
  437.      * This class implements the computation of the Hansen kernels
  438.      * through a power series in e² and that is using
  439.      * modified Newcomb operators. The reference formulae can be found
  440.      * in Danielson 2.7.3-10 and 3.3-5
  441.      * </p>
  442.      */
  443.     private static class HansenCoefficientsBySeries {

  445.         /** -n-1 coefficient. */
  446.         private final int mnm1;

  448.         /** s coefficient. */
  449.         private final int s;

  451.         /** j coefficient. */
  452.         private final int j;

  454.         /** Max power in e² for the Newcomb's series expansion. */
  455.         private final int maxNewcomb;

  457.         /** Polynomial representing the serie. */
  458.         private PolynomialFunction polynomial;

  460.         /** Factory for the DerivativeStructure instances. */
  461.         private final DSFactory factory;

  463.         /**
  464.          * Class constructor.
  465.          *
  466.          * @param mnm1 -n-1 value
  467.          * @param s s value
  468.          * @param j j value
  469.          * @param maxHansen max power of e² in series expansion
  470.          */
  471.         HansenCoefficientsBySeries(final int mnm1, final int s,
  472.                                           final int j, final int maxHansen) {
  473.             this.mnm1 = mnm1;
  474.             this.s = s;
  475.             this.j = j;
  476.             this.maxNewcomb = maxHansen;
  477.             this.polynomial = generatePolynomial();
  478.             this.factory = new DSFactory(1, 1);
  479.         }

  481.         /** Computes the value of Hansen kernel and its derivative at e².
  482.          * <p>
  483.          * The formulae applied are described in Danielson 2.7.3-10 and
  484.          * 3.3-5
  485.          * </p>
  486.          * @param e2 e²
  487.          * @param chi &Chi;
  488.          * @param chi2 &Chi;²
  489.          * @return the value of the Hansen coefficient and its derivative for e²
  490.          */
  491.         public DerivativeStructure getValue(final double e2, final double chi, final double chi2) {

  493.             //Estimation of the serie expansion at e2
  494.             final DerivativeStructure serie = polynomial.value(factory.variable(0, e2));

  496.             final double value      =  FastMath.pow(chi2, -mnm1 - 1) * serie.getValue() / chi;
  497.             final double coef       = -(mnm1 + 1.5);
  498.             final double derivative = coef * chi2 * value +
  499.                                       FastMath.pow(chi2, -mnm1 - 1) * serie.getPartialDerivative(1) / chi;
  500.             return factory.build(value, derivative);
  501.         }

  503.         /** Generate the serie expansion in e².
  504.          * <p>
  505.          * Generate the series expansion in e² used in the formulation
  506.          * of the Hansen kernel (see Danielson 2.7.3-10):
  507.          * &Sigma; Y<sup>ns</sup><sub>α+a,α+b</sub>
  508.          * *e<sup>2α</sup>
  509.          * </p>
  510.          * @return polynomial representing the power serie expansion
  511.          */
  512.         private PolynomialFunction generatePolynomial() {
  513.             // Initialization
  514.             final int aHT = FastMath.max(j - s, 0);
  515.             final int bHT = FastMath.max(s - j, 0);

  517.             final double[] coefficients = new double[maxNewcomb + 1];

  519.             //Loop for getting the Newcomb operators
  520.             for (int alphaHT = 0; alphaHT <= maxNewcomb; alphaHT++) {
  521.                 coefficients[alphaHT] =
  522.                         NewcombOperators.getValue(alphaHT + aHT, alphaHT + bHT, mnm1, s);
  523.             }

  525.             //Creation of the polynomial
  526.             return new PolynomialFunction(coefficients);
  527.         }
  528.     }

  530. }
Created with JaCoCo 0.7.9.201702052155