FieldHansenTesseralLinear.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 java.lang.reflect.Array;

  19. import org.hipparchus.Field;
  20. import org.hipparchus.RealFieldElement;
  21. import org.hipparchus.analysis.differentiation.FDSFactory;
  22. import org.hipparchus.analysis.differentiation.FieldDerivativeStructure;
  23. import org.hipparchus.analysis.differentiation.FieldGradient;
  24. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  25. import org.hipparchus.util.FastMath;
  26. import org.hipparchus.util.MathArrays;
  27. import org.orekit.propagation.semianalytical.dsst.utilities.NewcombOperators;

  28. /**
  29.  * Hansen coefficients K(t,n,s) for t!=0 and n < 0.
  30.  * <p>
  31.  * Implements Collins 4-236 or Danielson 2.7.3-(9) for Hansen Coefficients and
  32.  * Collins 4-240 for derivatives. The recursions are transformed into
  33.  * composition of linear transformations to obtain the associated polynomials
  34.  * for coefficients and their derivatives - see Petre's paper
  35.  *
  36.  * @author Petre Bazavan
  37.  * @author Lucian Barbulescu
  38.  * @author Bryan Cazabonne
  39.  */
  40. public class FieldHansenTesseralLinear <T extends RealFieldElement<T>> {

  41.     /** The number of coefficients that will be computed with a set of roots. */
  42.     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;

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

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

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

  54.     /** The minimum value for the order. */
  55.     private int Nmin;

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

  58.     /** The s coefficient. */
  59.     private int s;

  60.     /** The j coefficient. */
  61.     private int j;

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

  64.     /**
  65.      * The offset used to identify the polynomial that corresponds to a negative.
  66.      * value of n in the internal array that starts at 0
  67.      */
  68.     private int offset;

  69.     /** The objects used to calculate initial data by means of Newcomb operators. */
  70.     private FieldHansenCoefficientsBySeries<T>[] hansenInit;

  71.     /**
  72.      * Constructor.
  73.      *
  74.      * @param nMax the maximum (absolute) value of n parameter
  75.      * @param s s parameter
  76.      * @param j j parameter
  77.      * @param n0 the minimum (absolute) value of n
  78.      * @param maxHansen maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
  79.      * @param field field used by default
  80.      */
  81.     @SuppressWarnings("unchecked")
  82.     public FieldHansenTesseralLinear(final int nMax, final int s, final int j, final int n0,
  83.                                      final int maxHansen, final Field<T> field) {
  84.         //Initialize the fields
  85.         this.offset = nMax + 1;
  86.         this.Nmin = -nMax - 1;
  87.         this.N0 = -n0 - 4;
  88.         this.s = s;
  89.         this.j = j;

  90.         final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
  91.         //Ensure that only the needed terms are computed
  92.         this.hansenInit = (FieldHansenCoefficientsBySeries<T>[]) Array.newInstance(FieldHansenCoefficientsBySeries.class, maxRoots);
  93.         for (int i = 0; i < maxRoots; i++) {
  94.             this.hansenInit[i] = new FieldHansenCoefficientsBySeries<>(N0 - i + 3, s, j, maxHansen, field);
  95.         }

  96.         // The first 4 values are computed with series. No linear combination is needed.
  97.         final int size = N0 - Nmin;
  98.         this.numSlices = (int) FastMath.max(FastMath.ceil(((double) size) / SLICE), 1);
  99.         hansenRoot = MathArrays.buildArray(field, numSlices, 4);
  100.         hansenDerivRoot = MathArrays.buildArray(field, numSlices, 4);
  101.         if (size > 0) {
  102.             mpvec = new PolynomialFunction[size][];
  103.             mpvecDeriv = new PolynomialFunction[size][];

  104.             // Prepare the database of the associated polynomials
  105.             generatePolynomials();
  106.         }

  107.     }

  108.     /**
  109.      * Compute polynomial coefficient a.
  110.      *
  111.      *  <p>
  112.      *  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>
  113.      *  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²
  114.      *  </p>
  115.      *
  116.      *  <p>
  117.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  118.      *  </p>
  119.      *
  120.      * @param mnm1 -n-1
  121.      * @return the polynomial
  122.      */
  123.     private PolynomialFunction a(final int mnm1) {
  124.         // Collins 4-236, Danielson 2.7.3-(9)
  125.         final double r1 = (mnm1 + 2.) * (2. * mnm1 + 5.);
  126.         final double r2 = (2. + mnm1 + s) * (2. + mnm1 - s);
  127.         return new PolynomialFunction(new double[] {
  128.             0.0, 0.0, r1 / r2
  129.         });
  130.     }

  131.     /**
  132.      * Compute polynomial coefficient b.
  133.      *
  134.      *  <p>
  135.      *  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>
  136.      *  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²
  137.      *  </p>
  138.      *
  139.      *  <p>
  140.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  141.      *  </p>
  142.      *
  143.      * @param mnm1 -n-1
  144.      * @return the polynomial
  145.      */
  146.     private PolynomialFunction b(final int mnm1) {
  147.         // Collins 4-236, Danielson 2.7.3-(9)
  148.         final double r2 = (2. + mnm1 + s) * (2. + mnm1 - s);
  149.         final double d1 = (mnm1 + 3.) * 2. * j * s / (r2 * (mnm1 + 4.));
  150.         final double d2 = (mnm1 + 3.) * (mnm1 + 2.) / r2;
  151.         return new PolynomialFunction(new double[] {
  152.             0.0, -d1, -d2
  153.         });
  154.     }

  155.     /**
  156.      * Compute polynomial coefficient c.
  157.      *
  158.      *  <p>
  159.      *  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>
  160.      *  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²
  161.      *  </p>
  162.      *
  163.      *  <p>
  164.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  165.      *  </p>
  166.      *
  167.      * @param mnm1 -n-1
  168.      * @return the polynomial
  169.      */
  170.     private PolynomialFunction c(final int mnm1) {
  171.         // Collins 4-236, Danielson 2.7.3-(9)
  172.         final double r1 = j * j * (mnm1 + 2.);
  173.         final double r2 = (mnm1 + 4.) * (2. + mnm1 + s) * (2. + mnm1 - s);

  174.         return new PolynomialFunction(new double[] {
  175.             0.0, 0.0, r1 / r2
  176.         });
  177.     }

  178.     /**
  179.      * Compute polynomial coefficient d.
  180.      *
  181.      *  <p>
  182.      *  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²
  183.      *  </p>
  184.      *
  185.      *  <p>
  186.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  187.      *  </p>
  188.      *
  189.      * @param mnm1 -n-1
  190.      * @return the polynomial
  191.      */
  192.     private PolynomialFunction d(final int mnm1) {
  193.         // Collins 4-236, Danielson 2.7.3-(9)
  194.         return new PolynomialFunction(new double[] {
  195.             0.0, 0.0, 1.0
  196.         });
  197.     }

  198.     /**
  199.      * Compute polynomial coefficient f.
  200.      *
  201.      *  <p>
  202.      *  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²
  203.      *  </p>
  204.      *
  205.      *  <p>
  206.      *  See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
  207.      *  </p>
  208.      *
  209.      * @param n index
  210.      * @return the polynomial
  211.      */
  212.     private PolynomialFunction f(final int n) {
  213.         // Collins 4-236, Danielson 2.7.3-(9)
  214.         final double r1 = (n + 3.0) * j * s;
  215.         final double r2 = (n + 4.0) * (2.0 + n + s) * (2.0 + n - s);
  216.         return new PolynomialFunction(new double[] {
  217.             0.0, 0.0, 0.0, r1 / r2
  218.         });
  219.     }

  220.     /**
  221.      * Generate the polynomials needed in the linear transformation.
  222.      *
  223.      * <p>
  224.      * See Petre's paper
  225.      * </p>
  226.      */
  227.     private void generatePolynomials() {


  228.         // Initialization of the matrices for linear transformations
  229.         // The final configuration of these matrices are obtained by composition
  230.         // of linear transformations

  231.         // The matrix of polynomials associated to Hansen coefficients, Petre's
  232.         // paper
  233.         PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix4();

  234.         // The matrix of polynomials associated to derivatives, Petre's paper
  235.         final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix4();
  236.         PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix4();
  237.         final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix4();

  238.         // The matrix of the current linear transformation
  239.         a.setMatrixLine(0, new PolynomialFunction[] {
  240.             HansenUtilities.ZERO, HansenUtilities.ONE, HansenUtilities.ZERO, HansenUtilities.ZERO
  241.         });
  242.         a.setMatrixLine(1, new PolynomialFunction[] {
  243.             HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ONE, HansenUtilities.ZERO
  244.         });
  245.         a.setMatrixLine(2, new PolynomialFunction[] {
  246.             HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ONE
  247.         });
  248.         // The generation process
  249.         int index;
  250.         int sliceCounter = 0;
  251.         for (int i = N0 - 1; i > Nmin - 1; i--) {
  252.             index = i + this.offset;
  253.             // The matrix of the current linear transformation is updated
  254.             // Petre's paper
  255.             a.setMatrixLine(3, new PolynomialFunction[] {
  256.                     c(i), HansenUtilities.ZERO, b(i), a(i)
  257.             });

  258.             // composition of the linear transformations to calculate
  259.             // the polynomials associated to Hansen coefficients
  260.             // Petre's paper
  261.             A = A.multiply(a);
  262.             // store the polynomials for Hansen coefficients
  263.             mpvec[index] = A.getMatrixLine(3);
  264.             // composition of the linear transformations to calculate
  265.             // the polynomials associated to derivatives
  266.             // Petre's paper
  267.             D = D.multiply(a);

  268.             //Update the B matrix
  269.             B.setMatrixLine(3, new PolynomialFunction[] {
  270.                 HansenUtilities.ZERO, f(i),
  271.                 HansenUtilities.ZERO, d(i)
  272.             });
  273.             D = D.add(A.multiply(B));

  274.             // store the polynomials for Hansen coefficients from the
  275.             // expressions of derivatives
  276.             mpvecDeriv[index] = D.getMatrixLine(3);

  277.             if (++sliceCounter % SLICE == 0) {
  278.                 // Re-Initialisation of matrix for linear transformmations
  279.                 // The final configuration of these matrix are obtained by composition
  280.                 // of linear transformations
  281.                 A = HansenUtilities.buildIdentityMatrix4();
  282.                 D = HansenUtilities.buildZeroMatrix4();
  283.             }
  284.         }
  285.     }

  286.     /**
  287.      * Compute the values for the first four coefficients and their derivatives by means of series.
  288.      *
  289.      * @param e2 e²
  290.      * @param chi &Chi;
  291.      * @param chi2 &Chi;²
  292.      */
  293.     public void computeInitValues(final T e2, final T chi, final T chi2) {
  294.         // compute the values for n, n+1, n+2 and n+3 by series
  295.         // See Danielson 2.7.3-(10)
  296.         //Ensure that only the needed terms are computed
  297.         final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
  298.         for (int i = 0; i < maxRoots; i++) {
  299.             final FieldGradient<T> hansenKernel = hansenInit[i].getValueGradient(e2, chi, chi2);
  300.             this.hansenRoot[0][i] = hansenKernel.getValue();
  301.             this.hansenDerivRoot[0][i] = hansenKernel.getPartialDerivative(0);
  302.         }

  303.         for (int i = 1; i < numSlices; i++) {
  304.             for (int k = 0; k < 4; k++) {
  305.                 final PolynomialFunction[] mv = mpvec[N0 - (i * SLICE) - k + 3 + offset];
  306.                 final PolynomialFunction[] sv = mpvecDeriv[N0 - (i * SLICE) - k + 3 + offset];

  307.                 hansenDerivRoot[i][k] = mv[3].value(chi).multiply(hansenDerivRoot[i - 1][3]).
  308.                                         add(mv[2].value(chi).multiply(hansenDerivRoot[i - 1][2])).
  309.                                         add(mv[1].value(chi).multiply(hansenDerivRoot[i - 1][1])).
  310.                                         add(mv[0].value(chi).multiply(hansenDerivRoot[i - 1][0])).
  311.                                         add(sv[3].value(chi).multiply(hansenRoot[i - 1][3])).
  312.                                         add(sv[2].value(chi).multiply(hansenRoot[i - 1][2])).
  313.                                         add(sv[1].value(chi).multiply(hansenRoot[i - 1][1])).
  314.                                         add(sv[0].value(chi).multiply(hansenRoot[i - 1][0]));

  315.                 hansenRoot[i][k] =  mv[3].value(chi).multiply(hansenRoot[i - 1][3]).
  316.                                     add(mv[2].value(chi).multiply(hansenRoot[i - 1][2])).
  317.                                     add(mv[1].value(chi).multiply(hansenRoot[i - 1][1])).
  318.                                     add(mv[0].value(chi).multiply(hansenRoot[i - 1][0]));
  319.             }
  320.         }
  321.     }

  322.     /**
  323.      * Compute the value of the Hansen coefficient K<sub>j</sub><sup>-n-1, s</sup>.
  324.      *
  325.      * @param mnm1 -n-1
  326.      * @param chi χ
  327.      * @return the coefficient K<sub>j</sub><sup>-n-1, s</sup>
  328.      */
  329.     public T getValue(final int mnm1, final T chi) {
  330.         //Compute n
  331.         final int n = -mnm1 - 1;

  332.         //Compute the potential slice
  333.         int sliceNo = (n + N0 + 4) / SLICE;
  334.         if (sliceNo < numSlices) {
  335.             //Compute the index within the slice
  336.             final int indexInSlice = (n + N0 + 4) % SLICE;

  337.             //Check if a root must be returned
  338.             if (indexInSlice <= 3) {
  339.                 return hansenRoot[sliceNo][indexInSlice];
  340.             }
  341.         } else {
  342.             //the value was a potential root for a slice, but that slice was not required
  343.             //Decrease the slice number
  344.             sliceNo--;
  345.         }

  346.         // Computes the coefficient by linear transformation
  347.         // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
  348.         final PolynomialFunction[] v = mpvec[mnm1 + offset];
  349.         return v[3].value(chi).multiply(hansenRoot[sliceNo][3]).
  350.                add(v[2].value(chi).multiply(hansenRoot[sliceNo][2])).
  351.                add(v[1].value(chi).multiply(hansenRoot[sliceNo][1])).
  352.                add(v[0].value(chi).multiply(hansenRoot[sliceNo][0]));

  353.     }

  354.     /**
  355.      * Compute the value of the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de².
  356.      *
  357.      * @param mnm1 -n-1
  358.      * @param chi χ
  359.      * @return the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de²
  360.      */
  361.     public T getDerivative(final int mnm1, final T chi) {

  362.         //Compute n
  363.         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;

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

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

  382.         return v[3].value(chi).multiply(hansenDerivRoot[sliceNo][3]).
  383.                add(v[2].value(chi).multiply(hansenDerivRoot[sliceNo][2])).
  384.                add(v[1].value(chi).multiply(hansenDerivRoot[sliceNo][1])).
  385.                add(v[0].value(chi).multiply(hansenDerivRoot[sliceNo][0])).
  386.                add(vv[3].value(chi).multiply(hansenRoot[sliceNo][3])).
  387.                add(vv[2].value(chi).multiply(hansenRoot[sliceNo][2])).
  388.                add( vv[1].value(chi).multiply(hansenRoot[sliceNo][1])).
  389.                add(vv[0].value(chi).multiply(hansenRoot[sliceNo][0]));

  390.     }

  391.     /**
  392.      * Compute a Hansen coefficient with series.
  393.      * <p>
  394.      * This class implements the computation of the Hansen kernels
  395.      * through a power series in e² and that is using
  396.      * modified Newcomb operators. The reference formulae can be found
  397.      * in Danielson 2.7.3-10 and 3.3-5
  398.      * </p>
  399.      */
  400.     private static class FieldHansenCoefficientsBySeries <T extends RealFieldElement<T>> {

  401.         /** -n-1 coefficient. */
  402.         private final int mnm1;

  403.         /** s coefficient. */
  404.         private final int s;

  405.         /** j coefficient. */
  406.         private final int j;

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

  409.         /** Polynomial representing the serie. */
  410.         private PolynomialFunction polynomial;

  411.         /** Factory for the DerivativeStructure instances. */
  412.         private final FDSFactory<T> factory;

  413.         /**
  414.          * Class constructor.
  415.          *
  416.          * @param mnm1 -n-1 value
  417.          * @param s s value
  418.          * @param j j value
  419.          * @param maxHansen max power of e² in series expansion
  420.          * @param field field for the function parameters and value
  421.          */
  422.         FieldHansenCoefficientsBySeries(final int mnm1, final int s,
  423.                                           final int j, final int maxHansen, final Field<T> field) {
  424.             this.mnm1 = mnm1;
  425.             this.s = s;
  426.             this.j = j;
  427.             this.maxNewcomb = maxHansen;
  428.             this.polynomial = generatePolynomial();
  429.             this.factory = new FDSFactory<>(field, 1, 1);
  430.         }

  431.         /** Computes the value of Hansen kernel and its derivative at e².
  432.          * <p>
  433.          * The formulae applied are described in Danielson 2.7.3-10 and
  434.          * 3.3-5
  435.          * </p>
  436.          * @param e2 e²
  437.          * @param chi &Chi;
  438.          * @param chi2 &Chi;²
  439.          * @return the value of the Hansen coefficient and its derivative for e²
  440.          * @deprecated as of 10.2, replaced by {@link #getValueGradient(RealFieldElement, RealFieldElement, RealFieldElement)}
  441.          */
  442.         @SuppressWarnings("unused")
  443.         @Deprecated
  444.         public FieldDerivativeStructure<T> getValue(final T e2, final T chi, final T chi2) {

  445.             final T zero = e2.getField().getZero();
  446.             //Estimation of the serie expansion at e2
  447.             final FieldDerivativeStructure<T> serie = polynomial.value(factory.variable(0, e2));

  448.             final T value      =  FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getValue()).divide(chi);
  449.             final T coef       = zero.subtract(mnm1 + 1.5);
  450.             final T derivative = coef.multiply(chi2).multiply(value).
  451.                             add(FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getPartialDerivative(1)).divide(chi));
  452.             return factory.build(value, derivative);
  453.         }

  454.         /** Computes the value of Hansen kernel and its derivative at e².
  455.          * <p>
  456.          * The formulae applied are described in Danielson 2.7.3-10 and
  457.          * 3.3-5
  458.          * </p>
  459.          * @param e2 e²
  460.          * @param chi &Chi;
  461.          * @param chi2 &Chi;²
  462.          * @return the value of the Hansen coefficient and its derivative for e²
  463.          */
  464.         private FieldGradient<T> getValueGradient(final T e2, final T chi, final T chi2) {

  465.             final T zero = e2.getField().getZero();
  466.             //Estimation of the serie expansion at e2
  467.             final FieldGradient<T> serie = polynomial.value(FieldGradient.variable(1, 0, e2));

  468.             final T value      =  FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getValue()).divide(chi);
  469.             final T coef       = zero.subtract(mnm1 + 1.5);
  470.             final T derivative = coef.multiply(chi2).multiply(value).
  471.                             add(FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getPartialDerivative(0)).divide(chi));
  472.             return new FieldGradient<T>(value, derivative);
  473.         }

  474.         /** Generate the serie expansion in e².
  475.          * <p>
  476.          * Generate the series expansion in e² used in the formulation
  477.          * of the Hansen kernel (see Danielson 2.7.3-10):
  478.          * &Sigma; Y<sup>ns</sup><sub>α+a,α+b</sub>
  479.          * *e<sup>2α</sup>
  480.          * </p>
  481.          * @return polynomial representing the power serie expansion
  482.          */
  483.         private PolynomialFunction generatePolynomial() {
  484.             // Initialization
  485.             final int aHT = FastMath.max(j - s, 0);
  486.             final int bHT = FastMath.max(s - j, 0);

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

  488.             //Loop for getting the Newcomb operators
  489.             for (int alphaHT = 0; alphaHT <= maxNewcomb; alphaHT++) {
  490.                 coefficients[alphaHT] =
  491.                         NewcombOperators.getValue(alphaHT + aHT, alphaHT + bHT, mnm1, s);
  492.             }

  493.             //Creation of the polynomial
  494.             return new PolynomialFunction(coefficients);
  495.         }
  496.     }

  497. }