DSSTThirdBody.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.forces;

  19. import java.io.NotSerializableException;
  20. import java.io.Serializable;
  21. import java.util.ArrayList;
  22. import java.util.HashMap;
  23. import java.util.List;
  24. import java.util.Map;
  25. import java.util.Set;
  26. import java.util.SortedSet;
  27. import java.util.TreeMap;

  29. import org.hipparchus.analysis.differentiation.DSFactory;
  30. import org.hipparchus.analysis.differentiation.DerivativeStructure;
  31. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  32. import org.hipparchus.util.FastMath;
  33. import org.orekit.attitudes.AttitudeProvider;
  34. import org.orekit.bodies.CelestialBody;
  35. import org.orekit.errors.OrekitException;
  36. import org.orekit.orbits.Orbit;
  37. import org.orekit.propagation.SpacecraftState;
  38. import org.orekit.propagation.events.EventDetector;
  39. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  40. import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
  41. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
  42. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory.NSKey;
  43. import org.orekit.propagation.semianalytical.dsst.utilities.JacobiPolynomials;
  44. import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
  45. import org.orekit.propagation.semianalytical.dsst.utilities.UpperBounds;
  46. import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenThirdBodyLinear;
  47. import org.orekit.time.AbsoluteDate;
  48. import org.orekit.utils.TimeSpanMap;

  50. /** Third body attraction perturbation to the
  51.  *  {@link org.orekit.propagation.semianalytical.dsst.DSSTPropagator DSSTPropagator}.
  52.  *
  53.  *  @author Romain Di Costanzo
  54.  *  @author Pascal Parraud
  55.  */
  56. public class DSSTThirdBody implements DSSTForceModel {

  58.     /** Max power for summation. */
  59.     private static final int    MAX_POWER = 22;

  61.     /** Truncation tolerance for big, eccentric  orbits. */
  62.     private static final double BIG_TRUNCATION_TOLERANCE = 1.e-1;

  64.     /** Truncation tolerance for small orbits. */
  65.     private static final double SMALL_TRUNCATION_TOLERANCE = 1.9e-6;

  67.     /** Number of points for interpolation. */
  68.     private static final int    INTERPOLATION_POINTS = 3;

  70.     /** Maximum power for eccentricity used in short periodic computation. */
  71.     private static final int    MAX_ECCPOWER_SP = 4;

  73.     /** Retrograde factor I.
  74.      *  <p>
  75.      *  DSST model needs equinoctial orbit as internal representation.
  76.      *  Classical equinoctial elements have discontinuities when inclination
  77.      *  is close to zero. In this representation, I = +1. <br>
  78.      *  To avoid this discontinuity, another representation exists and equinoctial
  79.      *  elements can be expressed in a different way, called "retrograde" orbit.
  80.      *  This implies I = -1. <br>
  81.      *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  82.      *  But for the sake of consistency with the theory, the retrograde factor
  83.      *  has been kept in the formulas.
  84.      *  </p>
  85.      */
  86.     private static final int    I = 1;

  88.     /** The 3rd body to consider. */
  89.     private final CelestialBody    body;

  91.     /** Standard gravitational parameter μ for the body in m³/s². */
  92.     private final double           gm;

  94.     /** Factorial. */
  95.     private final double[]         fact;

  97.     /** V<sub>ns</sub> coefficients. */
  98.     private final TreeMap<NSKey, Double> Vns;

  100.     /** Distance from center of mass of the central body to the 3rd body. */
  101.     private double R3;

  103.     /** Short period terms. */
  104.     private ThirdBodyShortPeriodicCoefficients shortPeriods;

  106.     // Equinoctial elements (according to DSST notation)
  107.     /** a. */
  108.     private double a;
  109.     /** e<sub>x</sub>. */
  110.     private double k;
  111.     /** e<sub>y</sub>. */
  112.     private double h;
  113.     /** h<sub>x</sub>. */
  114.     private double q;
  115.     /** h<sub>y</sub>. */
  116.     private double p;

  118.     /** Eccentricity. */
  119.     private double ecc;

  121.     // Direction cosines of the symmetry axis
  122.     /** α. */
  123.     private double alpha;
  124.     /** β. */
  125.     private double beta;
  126.     /** γ. */
  127.     private double gamma;

  129.     // Common factors for potential computation
  130.     /** A = n * a². */
  131.     private double A;
  132.     /** B = sqrt(1 - e²). */
  133.     private double B;
  134.     /** C = 1 + p² + q². */
  135.     private double C;
  136.     /** B². */
  137.     private double BB;
  138.     /** B³. */
  139.     private double BBB;

  141.     /** The mean motion (n). */
  142.     private double meanMotion;

  144.     /** &Chi; = 1 / sqrt(1 - e²) = 1 / B. */
  145.     private double X;
  146.     /** &Chi;². */
  147.     private double XX;
  148.     /** &Chi;³. */
  149.     private double XXX;
  150.     /** -2 * a / A. */
  151.     private double m2aoA;
  152.     /** B / A. */
  153.     private double BoA;
  154.     /** 1 / (A * B). */
  155.     private double ooAB;
  156.     /** -C / (2 * A * B). */
  157.     private double mCo2AB;
  158.     /** B / A(1 + B). */
  159.     private double BoABpo;

  161.     /** Max power for a/R3 in the serie expansion. */
  162.     private int    maxAR3Pow;

  164.     /** Max power for e in the serie expansion. */
  165.     private int    maxEccPow;

  167.     /** Max power for e in the serie expansion (for short periodics). */
  168.     private int    maxEccPowShort;

  170.     /** Max frequency of F. */
  171.     private int    maxFreqF;

  173.     /** An array that contains the objects needed to build the Hansen coefficients. <br/>
  174.      * The index is s */
  175.     private final HansenThirdBodyLinear[] hansenObjects;

  177.     /** The current value of the U function. <br/>
  178.      * Needed for the short periodic contribution */
  179.     private double U;

  181.     /** Qns coefficients. */
  182.     private double[][] Qns;
  183.     /** a / R3 up to power maxAR3Pow. */
  184.     private double[] aoR3Pow;
  185.     /** mu3 / R3. */
  186.     private double muoR3;

  188.     /** b = 1 / (1 + sqrt(1 - e²)) = 1 / (1 + B).*/
  189.     private double b;

  191.     /** h * &Chi;³. */
  192.     private double hXXX;
  193.     /** k * &Chi;³. */
  194.     private double kXXX;

  196.     /** Factory for the DerivativeStructure instances. */
  197.     private final DSFactory factory;

  199.     /** Complete constructor.
  200.      *  @param body the 3rd body to consider
  201.      *  @see org.orekit.bodies.CelestialBodyFactory
  202.      */
  203.     public DSSTThirdBody(final CelestialBody body) {
  204.         this.body = body;
  205.         this.gm   = body.getGM();

  207.         this.maxAR3Pow = Integer.MIN_VALUE;
  208.         this.maxEccPow = Integer.MIN_VALUE;

  210.         this.Vns = CoefficientsFactory.computeVns(MAX_POWER);

  212.         // Factorials computation
  213.         final int dim = 2 * MAX_POWER;
  214.         this.fact = new double[dim];
  215.         fact[0] = 1.;
  216.         for (int i = 1; i < dim; i++) {
  217.             fact[i] = i * fact[i - 1];
  218.         }

  220.         //Initialise the HansenCoefficient generator
  221.         this.hansenObjects = new HansenThirdBodyLinear[MAX_POWER + 1];
  222.         for (int s = 0; s <= MAX_POWER; s++) {
  223.             this.hansenObjects[s] = new HansenThirdBodyLinear(MAX_POWER, s);
  224.         }

  226.         this.factory = new DSFactory(1, 1);

  228.     }

  230.     /** Get third body.
  231.      *  @return third body
  232.      */
  233.     public CelestialBody getBody() {
  234.         return body;
  235.     }

  237.     /** Computes the highest power of the eccentricity and the highest power
  238.      *  of a/R3 to appear in the truncated analytical power series expansion.
  239.      *  <p>
  240.      *  This method computes the upper value for the 3rd body potential and
  241.      *  determines the maximal powers for the eccentricity and a/R3 producing
  242.      *  potential terms bigger than a defined tolerance.
  243.      *  </p>
  244.      *  @param aux auxiliary elements related to the current orbit
  245.      *  @param meanOnly only mean elements will be used for the propagation
  246.      *  @throws OrekitException if some specific error occurs
  247.      */
  248.     @Override
  249.     public List<ShortPeriodTerms> initialize(final AuxiliaryElements aux, final boolean meanOnly)
  250.         throws OrekitException {

  252.         // Initializes specific parameters.
  253.         initializeStep(aux);

  255.         // Truncation tolerance.
  256.         final double aor = a / R3;
  257.         final double tol = ( aor > .3 || (aor > .15  && ecc > .25) ) ? BIG_TRUNCATION_TOLERANCE : SMALL_TRUNCATION_TOLERANCE;

  259.         // Utilities for truncation
  260.         // Set a lower bound for eccentricity
  261.         final double eo2  = FastMath.max(0.0025, ecc / 2.);
  262.         final double x2o2 = XX / 2.;
  263.         final double[] eccPwr = new double[MAX_POWER];
  264.         final double[] chiPwr = new double[MAX_POWER];
  265.         eccPwr[0] = 1.;
  266.         chiPwr[0] = X;
  267.         for (int i = 1; i < MAX_POWER; i++) {
  268.             eccPwr[i] = eccPwr[i - 1] * eo2;
  269.             chiPwr[i] = chiPwr[i - 1] * x2o2;
  270.         }

  272.         // Auxiliary quantities.
  273.         final double ao2rxx = aor / (2. * XX);
  274.         double xmuarn       = ao2rxx * ao2rxx * gm / (X * R3);
  275.         double term         = 0.;

  277.         // Compute max power for a/R3 and e.
  278.         maxAR3Pow = 2;
  279.         maxEccPow = 0;
  280.         int n     = 2;
  281.         int m     = 2;
  282.         int nsmd2 = 0;

  284.         do {
  285.             // Upper bound for Tnm.
  286.             term =  xmuarn *
  287.                    (fact[n + m] / (fact[nsmd2] * fact[nsmd2 + m])) *
  288.                    (fact[n + m + 1] / (fact[m] * fact[n + 1])) *
  289.                    (fact[n - m + 1] / fact[n + 1]) *
  290.                    eccPwr[m] * UpperBounds.getDnl(XX, chiPwr[m], n + 2, m);

  292.             if (term < tol) {
  293.                 if (m == 0) {
  294.                     break;
  295.                 } else if (m < 2) {
  296.                     xmuarn *= ao2rxx;
  297.                     m = 0;
  298.                     n++;
  299.                     nsmd2++;
  300.                 } else {
  301.                     m -= 2;
  302.                     nsmd2++;
  303.                 }
  304.             } else {
  305.                 maxAR3Pow = n;
  306.                 maxEccPow = FastMath.max(m, maxEccPow);
  307.                 xmuarn *= ao2rxx;
  308.                 m++;
  309.                 n++;
  310.             }
  311.         } while (n < MAX_POWER);

  313.         maxEccPow = FastMath.min(maxAR3Pow, maxEccPow);

  315.         // allocate the array aoR3Pow
  316.         aoR3Pow = new double[maxAR3Pow + 1];

  318.         maxFreqF = maxAR3Pow + 1;
  319.         maxEccPowShort = MAX_ECCPOWER_SP;

  321.         Qns = CoefficientsFactory.computeQns(gamma, maxAR3Pow, FastMath.max(maxEccPow, maxEccPowShort));
  322.         final int jMax = maxAR3Pow + 1;
  323.         shortPeriods = new ThirdBodyShortPeriodicCoefficients(jMax, INTERPOLATION_POINTS,
  324.                                                               maxFreqF, body.getName(),
  325.                                                               new TimeSpanMap<Slot>(new Slot(jMax, INTERPOLATION_POINTS)));

  327.         final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
  328.         list.add(shortPeriods);
  329.         return list;

  331.     }

  333.     /** {@inheritDoc} */
  334.     @Override
  335.     public void initializeStep(final AuxiliaryElements aux) throws OrekitException {

  337.         // Equinoctial elements
  338.         a = aux.getSma();
  339.         k = aux.getK();
  340.         h = aux.getH();
  341.         q = aux.getQ();
  342.         p = aux.getP();

  344.         // Eccentricity
  345.         ecc = aux.getEcc();

  347.         // Distance from center of mass of the central body to the 3rd body
  348.         final Vector3D bodyPos = body.getPVCoordinates(aux.getDate(), aux.getFrame()).getPosition();
  349.         R3 = bodyPos.getNorm();

  351.         // Direction cosines
  352.         final Vector3D bodyDir = bodyPos.normalize();
  353.         alpha = bodyDir.dotProduct(aux.getVectorF());
  354.         beta  = bodyDir.dotProduct(aux.getVectorG());
  355.         gamma = bodyDir.dotProduct(aux.getVectorW());

  357.         // Equinoctial coefficients
  358.         A = aux.getA();
  359.         B = aux.getB();
  360.         C = aux.getC();
  361.         meanMotion = aux.getMeanMotion();

  363.         //&Chi;<sup>-2</sup>.
  364.         BB = B * B;
  365.         //&Chi;<sup>-3</sup>.
  366.         BBB = BB * B;

  368.         //b = 1 / (1 + B)
  369.         b = 1. / (1. + B);

  371.         // &Chi;
  372.         X = 1. / B;
  373.         XX = X * X;
  374.         XXX = X * XX;
  375.         // -2 * a / A
  376.         m2aoA = -2. * a / A;
  377.         // B / A
  378.         BoA = B / A;
  379.         // 1 / AB
  380.         ooAB = 1. / (A * B);
  381.         // -C / 2AB
  382.         mCo2AB = -C * ooAB / 2.;
  383.         // B / A(1 + B)
  384.         BoABpo = BoA / (1. + B);

  386.         // mu3 / R3
  387.         muoR3 = gm / R3;

  389.         //h * &Chi;³
  390.         hXXX = h * XXX;
  391.         //k * &Chi;³
  392.         kXXX = k * XXX;
  393.     }

  395.     /** {@inheritDoc} */
  396.     @Override
  397.     public double[] getMeanElementRate(final SpacecraftState currentState) {

  399.         // Qns coefficients
  400.         Qns = CoefficientsFactory.computeQns(gamma, maxAR3Pow, maxEccPow);
  401.         // a / R3 up to power maxAR3Pow
  402.         final double aoR3 = a / R3;
  403.         aoR3Pow[0] = 1.;
  404.         for (int i = 1; i <= maxAR3Pow; i++) {
  405.             aoR3Pow[i] = aoR3 * aoR3Pow[i - 1];
  406.         }

  408.         // Compute potential U derivatives
  409.         final double[] dU  = computeUDerivatives();
  410.         final double dUda  = dU[0];
  411.         final double dUdk  = dU[1];
  412.         final double dUdh  = dU[2];
  413.         final double dUdAl = dU[3];
  414.         final double dUdBe = dU[4];
  415.         final double dUdGa = dU[5];

  417.         // Compute cross derivatives [Eq. 2.2-(8)]
  418.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  419.         final double UAlphaGamma   = alpha * dUdGa - gamma * dUdAl;
  420.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  421.         final double UBetaGamma    =  beta * dUdGa - gamma * dUdBe;
  422.         // Common factor
  423.         final double pUAGmIqUBGoAB = (p * UAlphaGamma - I * q * UBetaGamma) * ooAB;

  425.         // Compute mean elements rates [Eq. 3.1-(1)]
  426.         final double da =  0.;
  427.         final double dh =  BoA * dUdk + k * pUAGmIqUBGoAB;
  428.         final double dk = -BoA * dUdh - h * pUAGmIqUBGoAB;
  429.         final double dp =  mCo2AB * UBetaGamma;
  430.         final double dq =  mCo2AB * UAlphaGamma * I;
  431.         final double dM =  m2aoA * dUda + BoABpo * (h * dUdh + k * dUdk) + pUAGmIqUBGoAB;

  433.         return new double[] {da, dk, dh, dq, dp, dM};

  435.     }

  437.     /** {@inheritDoc} */
  438.     @Override
  439.     public void updateShortPeriodTerms(final SpacecraftState... meanStates)
  440.         throws OrekitException {

  442.         final Slot slot = shortPeriods.createSlot(meanStates);

  444.         for (final SpacecraftState meanState : meanStates) {

  446.             initializeStep(new AuxiliaryElements(meanState.getOrbit(), I));

  448.             // a / R3 up to power maxAR3Pow
  449.             final double aoR3 = a / R3;
  450.             aoR3Pow[0] = 1.;
  451.             for (int i = 1; i <= maxAR3Pow; i++) {
  452.                 aoR3Pow[i] = aoR3 * aoR3Pow[i - 1];
  453.             }

  455.             // Qns coefficients
  456.             Qns = CoefficientsFactory.computeQns(gamma, maxAR3Pow, FastMath.max(maxEccPow, maxEccPowShort));
  457.             final GeneratingFunctionCoefficients gfCoefs =
  458.                             new GeneratingFunctionCoefficients(maxAR3Pow, MAX_ECCPOWER_SP, maxAR3Pow + 1);

  460.             //Compute additional quantities
  461.             // 2 * a / An
  462.             final double ax2oAn  = -m2aoA / meanMotion;
  463.             // B / An
  464.             final double BoAn  = BoA / meanMotion;
  465.             // 1 / ABn
  466.             final double ooABn = ooAB / meanMotion;
  467.             // C / 2ABn
  468.             final double Co2ABn = -mCo2AB / meanMotion;
  469.             // B / (A * (1 + B) * n)
  470.             final double BoABpon = BoABpo / meanMotion;
  471.             // -3 / n²a² = -3 / nA
  472.             final double m3onA = -3 / (A * meanMotion);

  474.             //Compute the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients.
  475.             for (int j = 1; j < slot.cij.length; j++) {
  476.                 // First compute the C<sub>i</sub><sup>j</sup> coefficients
  477.                 final double[] currentCij = new double[6];

  479.                 // Compute the cross derivatives operator :
  480.                 final double SAlphaGammaCj   = alpha * gfCoefs.getdSdgammaCj(j) - gamma * gfCoefs.getdSdalphaCj(j);
  481.                 final double SAlphaBetaCj    = alpha * gfCoefs.getdSdbetaCj(j)  - beta  * gfCoefs.getdSdalphaCj(j);
  482.                 final double SBetaGammaCj    =  beta * gfCoefs.getdSdgammaCj(j) - gamma * gfCoefs.getdSdbetaCj(j);
  483.                 final double ShkCj           =     h * gfCoefs.getdSdkCj(j)     -  k    * gfCoefs.getdSdhCj(j);
  484.                 final double pSagmIqSbgoABnCj = (p * SAlphaGammaCj - I * q * SBetaGammaCj) * ooABn;
  485.                 final double ShkmSabmdSdlCj  =  ShkCj - SAlphaBetaCj - gfCoefs.getdSdlambdaCj(j);

  487.                 currentCij[0] =  ax2oAn * gfCoefs.getdSdlambdaCj(j);
  488.                 currentCij[1] =  -(BoAn * gfCoefs.getdSdhCj(j) + h * pSagmIqSbgoABnCj + k * BoABpon * gfCoefs.getdSdlambdaCj(j));
  489.                 currentCij[2] =    BoAn * gfCoefs.getdSdkCj(j) + k * pSagmIqSbgoABnCj - h * BoABpon * gfCoefs.getdSdlambdaCj(j);
  490.                 currentCij[3] =  Co2ABn * (q * ShkmSabmdSdlCj - I * SAlphaGammaCj);
  491.                 currentCij[4] =  Co2ABn * (p * ShkmSabmdSdlCj - SBetaGammaCj);
  492.                 currentCij[5] = -ax2oAn * gfCoefs.getdSdaCj(j) + BoABpon * (h * gfCoefs.getdSdhCj(j) + k * gfCoefs.getdSdkCj(j)) + pSagmIqSbgoABnCj + m3onA * gfCoefs.getSCj(j);

  494.                 // add the computed coefficients to the interpolators
  495.                 slot.cij[j].addGridPoint(meanState.getDate(), currentCij);

  497.                 // Compute the S<sub>i</sub><sup>j</sup> coefficients
  498.                 final double[] currentSij = new double[6];

  500.                 // Compute the cross derivatives operator :
  501.                 final double SAlphaGammaSj   = alpha * gfCoefs.getdSdgammaSj(j) - gamma * gfCoefs.getdSdalphaSj(j);
  502.                 final double SAlphaBetaSj    = alpha * gfCoefs.getdSdbetaSj(j)  - beta  * gfCoefs.getdSdalphaSj(j);
  503.                 final double SBetaGammaSj    =  beta * gfCoefs.getdSdgammaSj(j) - gamma * gfCoefs.getdSdbetaSj(j);
  504.                 final double ShkSj           =     h * gfCoefs.getdSdkSj(j)     -  k    * gfCoefs.getdSdhSj(j);
  505.                 final double pSagmIqSbgoABnSj = (p * SAlphaGammaSj - I * q * SBetaGammaSj) * ooABn;
  506.                 final double ShkmSabmdSdlSj  =  ShkSj - SAlphaBetaSj - gfCoefs.getdSdlambdaSj(j);

  508.                 currentSij[0] =  ax2oAn * gfCoefs.getdSdlambdaSj(j);
  509.                 currentSij[1] =  -(BoAn * gfCoefs.getdSdhSj(j) + h * pSagmIqSbgoABnSj + k * BoABpon * gfCoefs.getdSdlambdaSj(j));
  510.                 currentSij[2] =    BoAn * gfCoefs.getdSdkSj(j) + k * pSagmIqSbgoABnSj - h * BoABpon * gfCoefs.getdSdlambdaSj(j);
  511.                 currentSij[3] =  Co2ABn * (q * ShkmSabmdSdlSj - I * SAlphaGammaSj);
  512.                 currentSij[4] =  Co2ABn * (p * ShkmSabmdSdlSj - SBetaGammaSj);
  513.                 currentSij[5] = -ax2oAn * gfCoefs.getdSdaSj(j) + BoABpon * (h * gfCoefs.getdSdhSj(j) + k * gfCoefs.getdSdkSj(j)) + pSagmIqSbgoABnSj + m3onA * gfCoefs.getSSj(j);

  515.                 // add the computed coefficients to the interpolators
  516.                 slot.sij[j].addGridPoint(meanState.getDate(), currentSij);

  518.                 if (j == 1) {
  519.                     //Compute the C⁰ coefficients using Danielson 2.5.2-15a.
  520.                     final double[] value = new double[6];
  521.                     for (int i = 0; i < 6; ++i) {
  522.                         value[i] = currentCij[i] * k / 2. + currentSij[i] * h / 2.;
  523.                     }
  524.                     slot.cij[0].addGridPoint(meanState.getDate(), value);
  525.                 }
  526.             }
  527.         }
  528.     }

  530.     /** {@inheritDoc} */
  531.     @Override
  532.     public EventDetector[] getEventsDetectors() {
  533.         return null;
  534.     }

  536.     /** Compute potential derivatives.
  537.      *  @return derivatives of the potential with respect to orbital parameters
  538.      */
  539.     private double[] computeUDerivatives() {

  541.         // Gs and Hs coefficients
  542.         final double[][] GsHs = CoefficientsFactory.computeGsHs(k, h, alpha, beta, maxEccPow);

  544.         // Initialise U.
  545.         U = 0.;

  547.         // Potential derivatives
  548.         double dUda  = 0.;
  549.         double dUdk  = 0.;
  550.         double dUdh  = 0.;
  551.         double dUdAl = 0.;
  552.         double dUdBe = 0.;
  553.         double dUdGa = 0.;

  555.         for (int s = 0; s <= maxEccPow; s++) {
  556.             // initialise the Hansen roots
  557.             this.hansenObjects[s].computeInitValues(B, BB, BBB);

  559.             // Get the current Gs coefficient
  560.             final double gs = GsHs[0][s];

  562.             // Compute Gs partial derivatives from 3.1-(9)
  563.             double dGsdh  = 0.;
  564.             double dGsdk  = 0.;
  565.             double dGsdAl = 0.;
  566.             double dGsdBe = 0.;
  567.             if (s > 0) {
  568.                 // First get the G(s-1) and the H(s-1) coefficients
  569.                 final double sxGsm1 = s * GsHs[0][s - 1];
  570.                 final double sxHsm1 = s * GsHs[1][s - 1];
  571.                 // Then compute derivatives
  572.                 dGsdh  = beta  * sxGsm1 - alpha * sxHsm1;
  573.                 dGsdk  = alpha * sxGsm1 + beta  * sxHsm1;
  574.                 dGsdAl = k * sxGsm1 - h * sxHsm1;
  575.                 dGsdBe = h * sxGsm1 + k * sxHsm1;
  576.             }

  578.             // Kronecker symbol (2 - delta(0,s))
  579.             final double delta0s = (s == 0) ? 1. : 2.;

  581.             for (int n = FastMath.max(2, s); n <= maxAR3Pow; n++) {
  582.                 // (n - s) must be even
  583.                 if ((n - s) % 2 == 0) {
  584.                     // Extract data from previous computation :
  585.                     final double kns   = this.hansenObjects[s].getValue(n, B);
  586.                     final double dkns  = this.hansenObjects[s].getDerivative(n, B);

  588.                     final double vns   = Vns.get(new NSKey(n, s));
  589.                     final double coef0 = delta0s * aoR3Pow[n] * vns;
  590.                     final double coef1 = coef0 * Qns[n][s];
  591.                     final double coef2 = coef1 * kns;
  592.                     // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  593.                     // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
  594.                     final double dqns = (n == s) ? 0. : Qns[n][s + 1];

  596.                     //Compute U:
  597.                     U += coef2 * gs;

  599.                     // Compute dU / da :
  600.                     dUda += coef2 * n * gs;
  601.                     // Compute dU / dh
  602.                     dUdh += coef1 * (kns * dGsdh + hXXX * gs * dkns);
  603.                     // Compute dU / dk
  604.                     dUdk += coef1 * (kns * dGsdk + kXXX * gs * dkns);
  605.                     // Compute dU / dAlpha
  606.                     dUdAl += coef2 * dGsdAl;
  607.                     // Compute dU / dBeta
  608.                     dUdBe += coef2 * dGsdBe;
  609.                     // Compute dU / dGamma
  610.                     dUdGa += coef0 * kns * dqns * gs;
  611.                 }
  612.             }
  613.         }

  615.         // multiply by mu3 / R3
  616.         U *= muoR3;

  618.         return new double[] {
  619.             dUda  * muoR3 / a,
  620.             dUdk  * muoR3,
  621.             dUdh  * muoR3,
  622.             dUdAl * muoR3,
  623.             dUdBe * muoR3,
  624.             dUdGa * muoR3
  625.         };

  627.     }

  629.     /** {@inheritDoc} */
  630.     @Override
  631.     public void registerAttitudeProvider(final AttitudeProvider provider) {
  632.         //nothing is done since this contribution is not sensitive to attitude
  633.     }

  635.     /** Computes the C<sup>j</sup> and S<sup>j</sup> coefficients Danielson 4.2-(15,16)
  636.      * and their derivatives.
  637.      *  <p>
  638.      *  CS Mathematical Report $3.5.3.2
  639.      *  </p>
  640.      */
  641.     private class FourierCjSjCoefficients {

  643.         /** The coefficients G<sub>n, s</sub> and their derivatives. */
  644.         private final GnsCoefficients gns;

  646.         /** the coefficients e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and their derivatives by h and k. */
  647.         private final WnsjEtomjmsCoefficient wnsjEtomjmsCoefficient;

  649.         /** The terms containing the coefficients C<sub>j</sub> and S<sub>j</sub> of (α, β) or (k, h). */
  650.         private final CjSjAlphaBetaKH ABDECoefficients;

  652.         /** The Fourier coefficients C<sup>j</sup> and their derivatives.
  653.          * <p>
  654.          * Each column of the matrix contains the following values: <br/>
  655.          * - C<sup>j</sup> <br/>
  656.          * - dC<sup>j</sup> / da <br/>
  657.          * - dC<sup>j</sup> / dk <br/>
  658.          * - dC<sup>j</sup> / dh <br/>
  659.          * - dC<sup>j</sup> / dα <br/>
  660.          * - dC<sup>j</sup> / dβ <br/>
  661.          * - dC<sup>j</sup> / dγ <br/>
  662.          * </p>
  663.          */
  664.         private final double[][] cj;

  666.         /** The S<sup>j</sup> coefficients and their derivatives.
  667.          * <p>
  668.          * Each column of the matrix contains the following values: <br/>
  669.          * - S<sup>j</sup> <br/>
  670.          * - dS<sup>j</sup> / da <br/>
  671.          * - dS<sup>j</sup> / dk <br/>
  672.          * - dS<sup>j</sup> / dh <br/>
  673.          * - dS<sup>j</sup> / dα <br/>
  674.          * - dS<sup>j</sup> / dβ <br/>
  675.          * - dS<sup>j</sup> / dγ <br/>
  676.          * </p>
  677.          */
  678.         private final double[][] sj;

  680.         /** The Coefficients C<sup>j</sup><sub>,λ</sub>.
  681.          * <p>
  682.          * See Danielson 4.2-21
  683.          * </p>
  684.          */
  685.         private final double[] cjlambda;

  687.         /** The Coefficients S<sup>j</sup><sub>,λ</sub>.
  688.         * <p>
  689.         * See Danielson 4.2-21
  690.         * </p>
  691.         */
  692.         private final double[] sjlambda;

  694.         /** Maximum value for n. */
  695.         private final int nMax;

  697.         /** Maximum value for s. */
  698.         private final int sMax;

  700.         /** Maximum value for j. */
  701.         private final int jMax;

  703.         /**
  704.          * Private constructor.
  705.          *
  706.          * @param nMax maximum value for n index
  707.          * @param sMax maximum value for s index
  708.          * @param jMax maximum value for j index
  709.          */
  710.         FourierCjSjCoefficients(final int nMax, final int sMax, final int jMax) {
  711.             //Save parameters
  712.             this.nMax = nMax;
  713.             this.sMax = sMax;
  714.             this.jMax = jMax;

  716.             //Create objects
  717.             wnsjEtomjmsCoefficient = new WnsjEtomjmsCoefficient();
  718.             ABDECoefficients = new CjSjAlphaBetaKH();
  719.             gns = new GnsCoefficients(nMax, sMax);

  721.             //create arays
  722.             this.cj = new double[7][jMax + 1];
  723.             this.sj = new double[7][jMax + 1];
  724.             this.cjlambda = new double[jMax];
  725.             this.sjlambda = new double[jMax];

  727.             computeCoefficients();
  728.         }

  730.         /**
  731.          * Compute all coefficients.
  732.          */
  733.         private void computeCoefficients() {

  735.             for (int j = 1; j <= jMax; j++) {
  736.                 // initialise the coefficients
  737.                 for (int i = 0; i <= 6; i++) {
  738.                     cj[i][j] = 0.;
  739.                     sj[i][j] = 0.;
  740.                 }
  741.                 if (j < jMax) {
  742.                     // initialise the C<sup>j</sup><sub>,λ</sub> and S<sup>j</sup><sub>,λ</sub> coefficients
  743.                     cjlambda[j] = 0.;
  744.                     sjlambda[j] = 0.;
  745.                 }
  746.                 for (int s = 0; s <= sMax; s++) {

  748.                     // Compute the coefficients A<sub>j, s</sub>, B<sub>j, s</sub>, D<sub>j, s</sub> and E<sub>j, s</sub>
  749.                     ABDECoefficients.computeCoefficients(j, s);

  751.                     // compute starting value for n
  752.                     final int minN = FastMath.max(2, FastMath.max(j - 1, s));

  754.                     for (int n = minN; n <= nMax; n++) {
  755.                         // check if n-s is even
  756.                         if ((n - s) % 2 == 0) {
  757.                             // compute the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n+1, s</sup> and its derivatives
  758.                             final double[] wjnp1semjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(j, s, n + 1);

  760.                             // compute the coefficient e<sup>-|j-s|</sup>*w<sub>-j</sub><sup>n+1, s</sup> and its derivatives
  761.                             final double[] wmjnp1semjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(-j, s, n + 1);

  763.                             // compute common factors
  764.                             final double coef1 = -(wjnp1semjms[0] * ABDECoefficients.getCoefA() + wmjnp1semjms[0] * ABDECoefficients.getCoefB());
  765.                             final double coef2 =   wjnp1semjms[0] * ABDECoefficients.getCoefD() + wmjnp1semjms[0] * ABDECoefficients.getCoefE();

  767.                             //Compute C<sup>j</sup>
  768.                             cj[0][j] += gns.getGns(n, s) * coef1;
  769.                             //Compute dC<sup>j</sup> / da
  770.                             cj[1][j] += gns.getdGnsda(n, s) * coef1;
  771.                             //Compute dC<sup>j</sup> / dk
  772.                             cj[2][j] += -gns.getGns(n, s) *
  773.                                         (
  774.                                             wjnp1semjms[1] * ABDECoefficients.getCoefA() +
  775.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefAdk() +
  776.                                             wmjnp1semjms[1] * ABDECoefficients.getCoefB() +
  777.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefBdk()
  778.                                          );
  779.                             //Compute dC<sup>j</sup> / dh
  780.                             cj[3][j] += -gns.getGns(n, s) *
  781.                                         (
  782.                                             wjnp1semjms[2] * ABDECoefficients.getCoefA() +
  783.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefAdh() +
  784.                                             wmjnp1semjms[2] * ABDECoefficients.getCoefB() +
  785.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefBdh()
  786.                                          );
  787.                             //Compute dC<sup>j</sup> / dα
  788.                             cj[4][j] += -gns.getGns(n, s) *
  789.                                         (
  790.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefAdalpha() +
  791.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefBdalpha()
  792.                                         );
  793.                             //Compute dC<sup>j</sup> / dβ
  794.                             cj[5][j] += -gns.getGns(n, s) *
  795.                                         (
  796.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefAdbeta() +
  797.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefBdbeta()
  798.                                         );
  799.                             //Compute dC<sup>j</sup> / dγ
  800.                             cj[6][j] += gns.getdGnsdgamma(n, s) * coef1;

  802.                             //Compute S<sup>j</sup>
  803.                             sj[0][j] += gns.getGns(n, s) * coef2;
  804.                             //Compute dS<sup>j</sup> / da
  805.                             sj[1][j] += gns.getdGnsda(n, s) * coef2;
  806.                             //Compute dS<sup>j</sup> / dk
  807.                             sj[2][j] += gns.getGns(n, s) *
  808.                                         (
  809.                                             wjnp1semjms[1] * ABDECoefficients.getCoefD() +
  810.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefDdk() +
  811.                                             wmjnp1semjms[1] * ABDECoefficients.getCoefE() +
  812.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefEdk()
  813.                                          );
  814.                             //Compute dS<sup>j</sup> / dh
  815.                             sj[3][j] += gns.getGns(n, s) *
  816.                                         (
  817.                                             wjnp1semjms[2] * ABDECoefficients.getCoefD() +
  818.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefDdh() +
  819.                                             wmjnp1semjms[2] * ABDECoefficients.getCoefE() +
  820.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefEdh()
  821.                                          );
  822.                             //Compute dS<sup>j</sup> / dα
  823.                             sj[4][j] += gns.getGns(n, s) *
  824.                                         (
  825.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefDdalpha() +
  826.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefEdalpha()
  827.                                         );
  828.                             //Compute dS<sup>j</sup> / dβ
  829.                             sj[5][j] += gns.getGns(n, s) *
  830.                                         (
  831.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefDdbeta() +
  832.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefEdbeta()
  833.                                         );
  834.                             //Compute dS<sup>j</sup> / dγ
  835.                             sj[6][j] += gns.getdGnsdgamma(n, s) * coef2;

  837.                             //Check if n is greater or equal to j and j is at most jMax-1
  838.                             if (n >= j && j < jMax) {
  839.                                 // compute the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives
  840.                                 final double[] wjnsemjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(j, s, n);

  842.                                 // compute the coefficient e<sup>-|j-s|</sup>*w<sub>-j</sub><sup>n, s</sup> and its derivatives
  843.                                 final double[] wmjnsemjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(-j, s, n);

  845.                                 //Compute C<sup>j</sup><sub>,λ</sub>
  846.                                 cjlambda[j] += gns.getGns(n, s) * (wjnsemjms[0] * ABDECoefficients.getCoefD() + wmjnsemjms[0] * ABDECoefficients.getCoefE());
  847.                                 //Compute S<sup>j</sup><sub>,λ</sub>
  848.                                 sjlambda[j] += gns.getGns(n, s) * (wjnsemjms[0] * ABDECoefficients.getCoefA() + wmjnsemjms[0] * ABDECoefficients.getCoefB());
  849.                             }
  850.                         }
  851.                     }
  852.                 }
  853.                 // Divide by j
  854.                 for (int i = 0; i <= 6; i++) {
  855.                     cj[i][j] /= j;
  856.                     sj[i][j] /= j;
  857.                 }
  858.             }
  859.             //The C⁰ coefficients are not computed here.
  860.             //They are evaluated at the final point.

  862.             //C⁰<sub>,λ</sub>
  863.             cjlambda[0] = k * cjlambda[1] / 2. + h * sjlambda[1] / 2.;
  864.         }

  866.         /** Get the Fourier coefficient C<sup>j</sup>.
  867.          * @param j j index
  868.          * @return C<sup>j</sup>
  869.          */
  870.         public double getCj(final int j) {
  871.             return cj[0][j];
  872.         }

  874.         /** Get the derivative dC<sup>j</sup>/da.
  875.          * @param j j index
  876.          * @return dC<sup>j</sup>/da
  877.          */
  878.         public double getdCjda(final int j) {
  879.             return cj[1][j];
  880.         }

  882.         /** Get the derivative dC<sup>j</sup>/dk.
  883.          * @param j j index
  884.          * @return dC<sup>j</sup>/dk
  885.          */
  886.         public double getdCjdk(final int j) {
  887.             return cj[2][j];
  888.         }

  890.         /** Get the derivative dC<sup>j</sup>/dh.
  891.          * @param j j index
  892.          * @return dC<sup>j</sup>/dh
  893.          */
  894.         public double getdCjdh(final int j) {
  895.             return cj[3][j];
  896.         }

  898.         /** Get the derivative dC<sup>j</sup>/dα.
  899.          * @param j j index
  900.          * @return dC<sup>j</sup>/dα
  901.          */
  902.         public double getdCjdalpha(final int j) {
  903.             return cj[4][j];
  904.         }

  906.         /** Get the derivative dC<sup>j</sup>/dβ.
  907.          * @param j j index
  908.          * @return dC<sup>j</sup>/dβ
  909.          */
  910.         public double getdCjdbeta(final int j) {
  911.             return cj[5][j];
  912.         }

  914.         /** Get the derivative dC<sup>j</sup>/dγ.
  915.          * @param j j index
  916.          * @return dC<sup>j</sup>/dγ
  917.          */
  918.         public double getdCjdgamma(final int j) {
  919.             return cj[6][j];
  920.         }

  922.         /** Get the Fourier coefficient S<sup>j</sup>.
  923.          * @param j j index
  924.          * @return S<sup>j</sup>
  925.          */
  926.         public double getSj(final int j) {
  927.             return sj[0][j];
  928.         }

  930.         /** Get the derivative dS<sup>j</sup>/da.
  931.          * @param j j index
  932.          * @return dS<sup>j</sup>/da
  933.          */
  934.         public double getdSjda(final int j) {
  935.             return sj[1][j];
  936.         }

  938.         /** Get the derivative dS<sup>j</sup>/dk.
  939.          * @param j j index
  940.          * @return dS<sup>j</sup>/dk
  941.          */
  942.         public double getdSjdk(final int j) {
  943.             return sj[2][j];
  944.         }

  946.         /** Get the derivative dS<sup>j</sup>/dh.
  947.          * @param j j index
  948.          * @return dS<sup>j</sup>/dh
  949.          */
  950.         public double getdSjdh(final int j) {
  951.             return sj[3][j];
  952.         }

  954.         /** Get the derivative dS<sup>j</sup>/dα.
  955.          * @param j j index
  956.          * @return dS<sup>j</sup>/dα
  957.          */
  958.         public double getdSjdalpha(final int j) {
  959.             return sj[4][j];
  960.         }

  962.         /** Get the derivative dS<sup>j</sup>/dβ.
  963.          * @param j j index
  964.          * @return dS<sup>j</sup>/dβ
  965.          */
  966.         public double getdSjdbeta(final int j) {
  967.             return sj[5][j];
  968.         }

  970.         /** Get the derivative dS<sup>j</sup>/dγ.
  971.          * @param j j index
  972.          * @return dS<sup>j</sup>/dγ
  973.          */
  974.         public double getdSjdgamma(final int j) {
  975.             return sj[6][j];
  976.         }

  978.         /** Get the coefficient C⁰<sub>,λ</sub>.
  979.          * @return C⁰<sub>,λ</sub>
  980.          */
  981.         public double getC0Lambda() {
  982.             return cjlambda[0];
  983.         }

  985.         /** Get the coefficient C<sup>j</sup><sub>,λ</sub>.
  986.          * @param j j index
  987.          * @return C<sup>j</sup><sub>,λ</sub>
  988.          */
  989.         public double getCjLambda(final int j) {
  990.             if (j < 1 || j >= jMax) {
  991.                 return 0.;
  992.             }
  993.             return cjlambda[j];
  994.         }

  996.         /** Get the coefficient S<sup>j</sup><sub>,λ</sub>.
  997.          * @param j j index
  998.          * @return S<sup>j</sup><sub>,λ</sub>
  999.          */
  1000.         public double getSjLambda(final int j) {
  1001.             if (j < 1 || j >= jMax) {
  1002.                 return 0.;
  1003.             }
  1004.             return sjlambda[j];
  1005.         }
  1006.     }

  1008.     /** This class covers the coefficients e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and their derivatives by h and k.
  1009.      *
  1010.      * <p>
  1011.      * Starting from Danielson 4.2-9,10,11 and taking into account that fact that: <br />
  1012.      * c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
  1013.      * the expression e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup>
  1014.      * can be written as: <br/ >
  1015.      * - for |s| > |j| <br />
  1016.      * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
  1017.      *          (((n + s)!(n - s)!)/((n + j)!(n - j)!)) *
  1018.      *          (-b)<sup>|j-s|</sup> *
  1019.      *          ((1 - c²)<sup>n-|s|</sup>/(1 + c²)<sup>n</sup>) *
  1020.      *          P<sub>n-|s|</sub><sup>|j-s|, |j+s|</sup>(χ) <br />
  1021.      * <br />
  1022.      * - for |s| <= |j| <br />
  1023.      * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
  1024.      *          (-b)<sup>|j-s|</sup> *
  1025.      *          ((1 - c²)<sup>n-|j|</sup>/(1 + c²)<sup>n</sup>) *
  1026.      *          P<sub>n-|j|</sub><sup>|j-s|, |j+s|</sup>(χ)
  1027.      * </p>
  1028.      *
  1029.      * @author Lucian Barbulescu
  1030.      */
  1031.     private class WnsjEtomjmsCoefficient {

  1033.         /** The value c.
  1034.          * <p>
  1035.          *  c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
  1036.          * </p>
  1037.          *  */
  1038.         private final double c;

  1040.         /** db / dh. */
  1041.         private final double dbdh;

  1043.         /** db / dk. */
  1044.         private final double dbdk;

  1046.         /** dc / dh = e * db/dh + b * de/dh. */
  1047.         private final double dcdh;

  1049.         /** dc / dk = e * db/dk + b * de/dk. */
  1050.         private final double dcdk;

  1052.         /** The values (1 - c²)<sup>n</sup>. <br />
  1053.          * The maximum possible value for the power is N + 1 */
  1054.         private final double[] omc2tn;

  1056.         /** The values (1 + c²)<sup>n</sup>. <br />
  1057.          * The maximum possible value for the power is N + 1 */
  1058.         private final double[] opc2tn;

  1060.         /** The values b<sup>|j-s|</sup>. */
  1061.         private final double[] btjms;

  1063.         /**
  1064.          * Standard constructor.
  1065.          */
  1066.         WnsjEtomjmsCoefficient() {
  1067.             //initialise fields
  1068.             c = ecc * b;
  1069.             final double c2 = c * c;

  1071.             //b² * χ
  1072.             final double b2Chi = b * b * X;
  1073.             //Compute derivatives of b
  1074.             dbdh = h * b2Chi;
  1075.             dbdk = k * b2Chi;

  1077.             //Compute derivatives of c
  1078.             if (ecc == 0.0) {
  1079.                 // we are at a perfectly circular orbit singularity here
  1080.                 // we arbitrarily consider the perigee is along the X axis,
  1081.                 // i.e cos(ω + Ω) = h/ecc 1 and sin(ω + Ω) = k/ecc = 0
  1082.                 dcdh = ecc * dbdh + b;
  1083.                 dcdk = ecc * dbdk;
  1084.             } else {
  1085.                 dcdh = ecc * dbdh + b * h / ecc;
  1086.                 dcdk = ecc * dbdk + b * k / ecc;
  1087.             }

  1089.             //Compute the powers (1 - c²)<sup>n</sup> and (1 + c²)<sup>n</sup>
  1090.             omc2tn = new double[maxAR3Pow + maxFreqF + 2];
  1091.             opc2tn = new double[maxAR3Pow + maxFreqF + 2];
  1092.             final double omc2 = 1. - c2;
  1093.             final double opc2 = 1. + c2;
  1094.             omc2tn[0] = 1.;
  1095.             opc2tn[0] = 1.;
  1096.             for (int i = 1; i <= maxAR3Pow + maxFreqF + 1; i++) {
  1097.                 omc2tn[i] = omc2tn[i - 1] * omc2;
  1098.                 opc2tn[i] = opc2tn[i - 1] * opc2;
  1099.             }

  1101.             //Compute the powers of b
  1102.             btjms = new double[maxAR3Pow + maxFreqF + 1];
  1103.             btjms[0] = 1.;
  1104.             for (int i = 1; i <= maxAR3Pow + maxFreqF; i++) {
  1105.                 btjms[i] = btjms[i - 1] * b;
  1106.             }
  1107.         }

  1109.         /** Compute the value of the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives by h and k. <br />
  1110.          *
  1111.          * @param j j index
  1112.          * @param s s index
  1113.          * @param n n index
  1114.          * @return an array containing the value of the coefficient at index 0, the derivative by k at index 1 and the derivative by h at index 2
  1115.          */
  1116.         public double[] computeWjnsEmjmsAndDeriv(final int j, final int s, final int n) {
  1117.             final double[] wjnsemjms = new double[] {0., 0., 0.};

  1119.             // |j|
  1120.             final int absJ = FastMath.abs(j);
  1121.             // |s|
  1122.             final int absS = FastMath.abs(s);
  1123.             // |j - s|
  1124.             final int absJmS = FastMath.abs(j - s);
  1125.             // |j + s|
  1126.             final int absJpS = FastMath.abs(j + s);

  1128.             //The lower index of P. Also the power of (1 - c²)
  1129.             final int l;
  1130.             // the factorial ratio coefficient or 1. if |s| <= |j|
  1131.             final double factCoef;
  1132.             if (absS > absJ) {
  1133.                 factCoef = (fact[n + s] / fact[n + j]) * (fact[n - s] / fact[n - j]);
  1134.                 l = n - absS;
  1135.             } else {
  1136.                 factCoef = 1.;
  1137.                 l = n - absJ;
  1138.             }

  1140.             // (-1)<sup>|j-s|</sup>
  1141.             final double sign = absJmS % 2 != 0 ? -1. : 1.;
  1142.             //(1 - c²)<sup>n-|s|</sup> / (1 + c²)<sup>n</sup>
  1143.             final double coef1 = omc2tn[l] / opc2tn[n];
  1144.             //-b<sup>|j-s|</sup>
  1145.             final double coef2 = sign * btjms[absJmS];
  1146.             // P<sub>l</sub><sup>|j-s|, |j+s|</sup>(χ)
  1147.             final DerivativeStructure jac =
  1148.                     JacobiPolynomials.getValue(l, absJmS, absJpS, factory.variable(0, X));

  1150.             // the derivative of coef1 by c
  1151.             final double dcoef1dc = -coef1 * 2. * c * (((double) n) / opc2tn[1] + ((double) l) / omc2tn[1]);
  1152.             // the derivative of coef1 by h
  1153.             final double dcoef1dh = dcoef1dc * dcdh;
  1154.             // the derivative of coef1 by k
  1155.             final double dcoef1dk = dcoef1dc * dcdk;

  1157.             // the derivative of coef2 by b
  1158.             final double dcoef2db = absJmS == 0 ? 0 : sign * (double) absJmS * btjms[absJmS - 1];
  1159.             // the derivative of coef2 by h
  1160.             final double dcoef2dh = dcoef2db * dbdh;
  1161.             // the derivative of coef2 by k
  1162.             final double dcoef2dk = dcoef2db * dbdk;

  1164.             // the jacobi polynomial value
  1165.             final double jacobi = jac.getValue();
  1166.             // the derivative of the Jacobi polynomial by h
  1167.             final double djacobidh = jac.getPartialDerivative(1) * hXXX;
  1168.             // the derivative of the Jacobi polynomial by k
  1169.             final double djacobidk = jac.getPartialDerivative(1) * kXXX;

  1171.             //group the above coefficients to limit the mathematical operations
  1172.             final double term1 = factCoef * coef1 * coef2;
  1173.             final double term2 = factCoef * coef1 * jacobi;
  1174.             final double term3 = factCoef * coef2 * jacobi;

  1176.             //compute e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives by k and h
  1177.             wjnsemjms[0] = term1 * jacobi;
  1178.             wjnsemjms[1] = dcoef1dk * term3 + dcoef2dk * term2 + djacobidk * term1;
  1179.             wjnsemjms[2] = dcoef1dh * term3 + dcoef2dh * term2 + djacobidh * term1;

  1181.             return wjnsemjms;
  1182.         }
  1183.     }

  1185.     /** The G<sub>n,s</sub> coefficients and their derivatives.
  1186.      * <p>
  1187.      * See Danielson, 4.2-17
  1188.      *
  1189.      * @author Lucian Barbulescu
  1190.      */
  1191.     private class GnsCoefficients {

  1193.         /** Maximum value for n index. */
  1194.         private final int nMax;

  1196.         /** Maximum value for s index. */
  1197.         private final int sMax;

  1199.         /** The coefficients G<sub>n,s</sub>. */
  1200.         private final double gns[][];

  1202.         /** The derivatives of the coefficients G<sub>n,s</sub> by a. */
  1203.         private final double dgnsda[][];

  1205.         /** The derivatives of the coefficients G<sub>n,s</sub> by γ. */
  1206.         private final double dgnsdgamma[][];

  1208.         /** Standard constructor.
  1209.          *
  1210.          * @param nMax maximum value for n indes
  1211.          * @param sMax maximum value for s index
  1212.          */
  1213.         GnsCoefficients(final int nMax, final int sMax) {
  1214.             this.nMax = nMax;
  1215.             this.sMax = sMax;

  1217.             final int rows    = nMax + 1;
  1218.             final int columns = sMax + 1;
  1219.             this.gns          = new double[rows][columns];
  1220.             this.dgnsda       = new double[rows][columns];
  1221.             this.dgnsdgamma   = new double[rows][columns];

  1223.             // Generate the coefficients
  1224.             generateCoefficients();
  1225.         }
  1226.         /**
  1227.          * Compute the coefficient G<sub>n,s</sub> and its derivatives.
  1228.          * <p>
  1229.          * Only the derivatives by a and γ are computed as all others are 0
  1230.          * </p>
  1231.          */
  1232.         private void generateCoefficients() {
  1233.             for (int s = 0; s <= sMax; s++) {
  1234.                 // The n index is always at least the maximum between 2 and s
  1235.                 final int minN = FastMath.max(2, s);
  1236.                 for (int n = minN; n <= nMax; n++) {
  1237.                     // compute the coefficients only if (n - s) % 2 == 0
  1238.                     if ( (n - s) % 2 == 0 ) {
  1239.                         // Kronecker symbol (2 - delta(0,s))
  1240.                         final double delta0s = (s == 0) ? 1. : 2.;
  1241.                         final double vns   = Vns.get(new NSKey(n, s));
  1242.                         final double coef0 = delta0s * aoR3Pow[n] * vns * muoR3;
  1243.                         final double coef1 = coef0 * Qns[n][s];
  1244.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  1245.                         // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
  1246.                         final double dqns = (n == s) ? 0. : Qns[n][s + 1];

  1248.                         //Compute the coefficient and its derivatives.
  1249.                         this.gns[n][s] = coef1;
  1250.                         this.dgnsda[n][s] = coef1 * n / a;
  1251.                         this.dgnsdgamma[n][s] = coef0 * dqns;
  1252.                     } else {
  1253.                         // the coefficient and its derivatives is 0
  1254.                         this.gns[n][s] = 0.;
  1255.                         this.dgnsda[n][s] = 0.;
  1256.                         this.dgnsdgamma[n][s] = 0.;
  1257.                     }
  1258.                 }
  1259.             }
  1260.         }

  1262.         /** Get the coefficient G<sub>n,s</sub>.
  1263.          *
  1264.          * @param n n index
  1265.          * @param s s index
  1266.          * @return the coefficient G<sub>n,s</sub>
  1267.          */
  1268.         public double getGns(final int n, final int s) {
  1269.             return this.gns[n][s];
  1270.         }

  1272.         /** Get the derivative dG<sub>n,s</sub> / da.
  1273.          *
  1274.          * @param n n index
  1275.          * @param s s index
  1276.          * @return the derivative dG<sub>n,s</sub> / da
  1277.          */
  1278.         public double getdGnsda(final int n, final int s) {
  1279.             return this.dgnsda[n][s];
  1280.         }

  1282.         /** Get the derivative dG<sub>n,s</sub> / dγ.
  1283.          *
  1284.          * @param n n index
  1285.          * @param s s index
  1286.          * @return the derivative dG<sub>n,s</sub> / dγ
  1287.          */
  1288.         public double getdGnsdgamma(final int n, final int s) {
  1289.             return this.dgnsdgamma[n][s];
  1290.         }
  1291.     }

  1293.     /** This class computes the terms containing the coefficients C<sub>j</sub> and S<sub>j</sub> of (α, β) or (k, h).
  1294.      *
  1295.      * <p>
  1296.      * The following terms and their derivatives by k, h, alpha and beta are considered: <br/ >
  1297.      * - sign(j-s) * C<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h) + S<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h) <br />
  1298.      * - C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) <br />
  1299.      * - C<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h) - sign(j-s) * S<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h) <br />
  1300.      * - C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) <br />
  1301.      * For the ease of usage the above terms are renamed A<sub>js</sub>, B<sub>js</sub>, D<sub>js</sub> and E<sub>js</sub> respectively <br />
  1302.      * See the CS Mathematical report $3.5.3.2 for more details
  1303.      * </p>
  1304.      * @author Lucian Barbulescu
  1305.      */
  1306.     private class CjSjAlphaBetaKH {

  1308.         /** The C<sub>j</sub>(k, h) and the S<sub>j</sub>(k, h) series. */
  1309.         private final CjSjCoefficient cjsjkh;

  1311.         /** The C<sub>j</sub>(α, β) and the S<sub>j</sub>(α, β) series. */
  1312.         private final CjSjCoefficient cjsjalbe;

  1314.         /** The coeficient sign(j-s) * C<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h) + S<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h)
  1315.          * and its derivative by k, h, α and β. */
  1316.         private final double coefAandDeriv[];

  1318.         /** The coeficient C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h)
  1319.          * and its derivative by k, h, α and β. */
  1320.         private final double coefBandDeriv[];

  1322.         /** The coeficient C<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h) - sign(j-s) * S<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h)
  1323.          * and its derivative by k, h, α and β. */
  1324.         private final double coefDandDeriv[];

  1326.         /** The coeficient C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h)
  1327.          * and its derivative by k, h, α and β. */
  1328.         private final double coefEandDeriv[];

  1330.         /**
  1331.          * Standard constructor.
  1332.          */
  1333.         CjSjAlphaBetaKH() {
  1334.             cjsjkh = new CjSjCoefficient(k, h);
  1335.             cjsjalbe = new CjSjCoefficient(alpha, beta);

  1337.             coefAandDeriv = new double[5];
  1338.             coefBandDeriv = new double[5];
  1339.             coefDandDeriv = new double[5];
  1340.             coefEandDeriv = new double[5];
  1341.         }

  1343.         /** Compute the coefficients and their derivatives for a given (j,s) pair.
  1344.          * @param j j index
  1345.          * @param s s index
  1346.          */
  1347.         public void computeCoefficients(final int j, final int s) {
  1348.             // sign of j-s
  1349.             final int sign = j < s ? -1 : 1;

  1351.             //|j-s|
  1352.             final int absJmS = FastMath.abs(j - s);

  1354.             //j+s
  1355.             final int jps = j + s;

  1357.             //Compute the coefficient A and its derivatives
  1358.             coefAandDeriv[0] = sign * cjsjalbe.getCj(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getCj(absJmS);
  1359.             coefAandDeriv[1] = sign * cjsjalbe.getCj(s) * cjsjkh.getDsjDk(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getDcjDk(absJmS);
  1360.             coefAandDeriv[2] = sign * cjsjalbe.getCj(s) * cjsjkh.getDsjDh(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getDcjDh(absJmS);
  1361.             coefAandDeriv[3] = sign * cjsjalbe.getDcjDk(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getDsjDk(s) * cjsjkh.getCj(absJmS);
  1362.             coefAandDeriv[4] = sign * cjsjalbe.getDcjDh(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getDsjDh(s) * cjsjkh.getCj(absJmS);

  1364.             //Compute the coefficient B and its derivatives
  1365.             coefBandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getSj(jps) - cjsjalbe.getSj(s) * cjsjkh.getCj(jps);
  1366.             coefBandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDsjDk(jps) - cjsjalbe.getSj(s) * cjsjkh.getDcjDk(jps);
  1367.             coefBandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDsjDh(jps) - cjsjalbe.getSj(s) * cjsjkh.getDcjDh(jps);
  1368.             coefBandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getSj(jps) - cjsjalbe.getDsjDk(s) * cjsjkh.getCj(jps);
  1369.             coefBandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getSj(jps) - cjsjalbe.getDsjDh(s) * cjsjkh.getCj(jps);

  1371.             //Compute the coefficient D and its derivatives
  1372.             coefDandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getSj(absJmS);
  1373.             coefDandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDcjDk(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getDsjDk(absJmS);
  1374.             coefDandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDcjDh(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getDsjDh(absJmS);
  1375.             coefDandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getDsjDk(s) * cjsjkh.getSj(absJmS);
  1376.             coefDandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getDsjDh(s) * cjsjkh.getSj(absJmS);

  1378.             //Compute the coefficient E and its derivatives
  1379.             coefEandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getCj(jps) + cjsjalbe.getSj(s) * cjsjkh.getSj(jps);
  1380.             coefEandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDcjDk(jps) + cjsjalbe.getSj(s) * cjsjkh.getDsjDk(jps);
  1381.             coefEandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDcjDh(jps) + cjsjalbe.getSj(s) * cjsjkh.getDsjDh(jps);
  1382.             coefEandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getCj(jps) + cjsjalbe.getDsjDk(s) * cjsjkh.getSj(jps);
  1383.             coefEandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getCj(jps) + cjsjalbe.getDsjDh(s) * cjsjkh.getSj(jps);
  1384.         }

  1386.         /** Get the value of coefficient A<sub>j,s</sub>.
  1387.          *
  1388.          * @return the coefficient A<sub>j,s</sub>
  1389.          */
  1390.         public double getCoefA() {
  1391.             return coefAandDeriv[0];
  1392.         }

  1394.         /** Get the value of coefficient dA<sub>j,s</sub>/dk.
  1395.          *
  1396.          * @return the coefficient dA<sub>j,s</sub>/dk
  1397.          */
  1398.         public double getdCoefAdk() {
  1399.             return coefAandDeriv[1];
  1400.         }

  1402.         /** Get the value of coefficient dA<sub>j,s</sub>/dh.
  1403.          *
  1404.          * @return the coefficient dA<sub>j,s</sub>/dh
  1405.          */
  1406.         public double getdCoefAdh() {
  1407.             return coefAandDeriv[2];
  1408.         }

  1410.         /** Get the value of coefficient dA<sub>j,s</sub>/dα.
  1411.          *
  1412.          * @return the coefficient dA<sub>j,s</sub>/dα
  1413.          */
  1414.         public double getdCoefAdalpha() {
  1415.             return coefAandDeriv[3];
  1416.         }

  1418.         /** Get the value of coefficient dA<sub>j,s</sub>/dβ.
  1419.          *
  1420.          * @return the coefficient dA<sub>j,s</sub>/dβ
  1421.          */
  1422.         public double getdCoefAdbeta() {
  1423.             return coefAandDeriv[4];
  1424.         }

  1426.         /** Get the value of coefficient B<sub>j,s</sub>.
  1427.          *
  1428.          * @return the coefficient B<sub>j,s</sub>
  1429.          */
  1430.         public double getCoefB() {
  1431.             return coefBandDeriv[0];
  1432.         }

  1434.         /** Get the value of coefficient dB<sub>j,s</sub>/dk.
  1435.          *
  1436.          * @return the coefficient dB<sub>j,s</sub>/dk
  1437.          */
  1438.         public double getdCoefBdk() {
  1439.             return coefBandDeriv[1];
  1440.         }

  1442.         /** Get the value of coefficient dB<sub>j,s</sub>/dh.
  1443.          *
  1444.          * @return the coefficient dB<sub>j,s</sub>/dh
  1445.          */
  1446.         public double getdCoefBdh() {
  1447.             return coefBandDeriv[2];
  1448.         }

  1450.         /** Get the value of coefficient dB<sub>j,s</sub>/dα.
  1451.          *
  1452.          * @return the coefficient dB<sub>j,s</sub>/dα
  1453.          */
  1454.         public double getdCoefBdalpha() {
  1455.             return coefBandDeriv[3];
  1456.         }

  1458.         /** Get the value of coefficient dB<sub>j,s</sub>/dβ.
  1459.          *
  1460.          * @return the coefficient dB<sub>j,s</sub>/dβ
  1461.          */
  1462.         public double getdCoefBdbeta() {
  1463.             return coefBandDeriv[4];
  1464.         }

  1466.         /** Get the value of coefficient D<sub>j,s</sub>.
  1467.          *
  1468.          * @return the coefficient D<sub>j,s</sub>
  1469.          */
  1470.         public double getCoefD() {
  1471.             return coefDandDeriv[0];
  1472.         }

  1474.         /** Get the value of coefficient dD<sub>j,s</sub>/dk.
  1475.          *
  1476.          * @return the coefficient dD<sub>j,s</sub>/dk
  1477.          */
  1478.         public double getdCoefDdk() {
  1479.             return coefDandDeriv[1];
  1480.         }

  1482.         /** Get the value of coefficient dD<sub>j,s</sub>/dh.
  1483.          *
  1484.          * @return the coefficient dD<sub>j,s</sub>/dh
  1485.          */
  1486.         public double getdCoefDdh() {
  1487.             return coefDandDeriv[2];
  1488.         }

  1490.         /** Get the value of coefficient dD<sub>j,s</sub>/dα.
  1491.          *
  1492.          * @return the coefficient dD<sub>j,s</sub>/dα
  1493.          */
  1494.         public double getdCoefDdalpha() {
  1495.             return coefDandDeriv[3];
  1496.         }

  1498.         /** Get the value of coefficient dD<sub>j,s</sub>/dβ.
  1499.          *
  1500.          * @return the coefficient dD<sub>j,s</sub>/dβ
  1501.          */
  1502.         public double getdCoefDdbeta() {
  1503.             return coefDandDeriv[4];
  1504.         }

  1506.         /** Get the value of coefficient E<sub>j,s</sub>.
  1507.          *
  1508.          * @return the coefficient E<sub>j,s</sub>
  1509.          */
  1510.         public double getCoefE() {
  1511.             return coefEandDeriv[0];
  1512.         }

  1514.         /** Get the value of coefficient dE<sub>j,s</sub>/dk.
  1515.          *
  1516.          * @return the coefficient dE<sub>j,s</sub>/dk
  1517.          */
  1518.         public double getdCoefEdk() {
  1519.             return coefEandDeriv[1];
  1520.         }

  1522.         /** Get the value of coefficient dE<sub>j,s</sub>/dh.
  1523.          *
  1524.          * @return the coefficient dE<sub>j,s</sub>/dh
  1525.          */
  1526.         public double getdCoefEdh() {
  1527.             return coefEandDeriv[2];
  1528.         }

  1530.         /** Get the value of coefficient dE<sub>j,s</sub>/dα.
  1531.          *
  1532.          * @return the coefficient dE<sub>j,s</sub>/dα
  1533.          */
  1534.         public double getdCoefEdalpha() {
  1535.             return coefEandDeriv[3];
  1536.         }

  1538.         /** Get the value of coefficient dE<sub>j,s</sub>/dβ.
  1539.          *
  1540.          * @return the coefficient dE<sub>j,s</sub>/dβ
  1541.          */
  1542.         public double getdCoefEdbeta() {
  1543.             return coefEandDeriv[4];
  1544.         }
  1545.     }

  1547.     /** This class computes the coefficients for the generating function S and its derivatives.
  1548.      * <p>
  1549.      * The form of the generating functions is: <br>
  1550.      *  S = C⁰ + &Sigma;<sub>j=1</sub><sup>N+1</sup>(C<sup>j</sup> * cos(jF) + S<sup>j</sup> * sin(jF)) <br>
  1551.      *  The coefficients C⁰, C<sup>j</sup>, S<sup>j</sup> are the Fourrier coefficients
  1552.      *  presented in Danielson 4.2-14,15 except for the case j=1 where
  1553.      *  C¹ = C¹<sub>Fourier</sub> - hU and
  1554.      *  S¹ = S¹<sub>Fourier</sub> + kU <br>
  1555.      *  Also the coefficients of the derivatives of S by a, k, h, α, β, γ and λ
  1556.      *  are computed end expressed in a similar manner. The formulas used are 4.2-19, 20, 23, 24
  1557.      * </p>
  1558.      * @author Lucian Barbulescu
  1559.      */
  1560.     private class GeneratingFunctionCoefficients {

  1562.         /** The Fourier coefficients as presented in Danielson 4.2-14,15. */
  1563.         private final FourierCjSjCoefficients cjsjFourier;

  1565.         /** Maximum value of j index. */
  1566.         private final int jMax;

  1568.         /** The coefficients C<sup>j</sup> of the function S and its derivatives.
  1569.          * <p>
  1570.          * The index j belongs to the interval [0,jMax]. The coefficient C⁰ is the free coefficient.<br>
  1571.          * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
  1572.          * - S <br/>
  1573.          * - dS / da <br/>
  1574.          * - dS / dk <br/>
  1575.          * - dS / dh <br/>
  1576.          * - dS / dα <br/>
  1577.          * - dS / dβ <br/>
  1578.          * - dS / dγ <br/>
  1579.          * - dS / dλ
  1580.          * </p>
  1581.          */
  1582.         private final double[][] cjCoefs;

  1584.         /** The coefficients S<sup>j</sup> of the function S and its derivatives.
  1585.          * <p>
  1586.          * The index j belongs to the interval [0,jMax].<br>
  1587.          * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
  1588.          * - S <br/>
  1589.          * - dS / da <br/>
  1590.          * - dS / dk <br/>
  1591.          * - dS / dh <br/>
  1592.          * - dS / dα <br/>
  1593.          * - dS / dβ <br/>
  1594.          * - dS / dγ <br/>
  1595.          * - dS / dλ
  1596.          * </p>
  1597.          */
  1598.         private final double[][] sjCoefs;

  1600.         /**
  1601.          * Standard constructor.
  1602.          *
  1603.          * @param nMax maximum value of n index
  1604.          * @param sMax maximum value of s index
  1605.          * @param jMax maximum value of j index
  1606.          */
  1607.         GeneratingFunctionCoefficients(final int nMax, final int sMax, final int jMax) {
  1608.             this.jMax = jMax;
  1609.             this.cjsjFourier = new FourierCjSjCoefficients(nMax, sMax, jMax);
  1610.             this.cjCoefs = new double[8][jMax + 1];
  1611.             this.sjCoefs = new double[8][jMax + 1];

  1613.             computeGeneratingFunctionCoefficients();
  1614.         }

  1616.         /**
  1617.          * Compute the coefficients for the generating function S and its derivatives.
  1618.          */
  1619.         private void computeGeneratingFunctionCoefficients() {

  1621.             // Compute potential U and derivatives
  1622.             final double[] dU  = computeUDerivatives();

  1624.             //Compute the C<sup>j</sup> coefficients
  1625.             for (int j = 1; j <= jMax; j++) {
  1626.                 //Compute the C<sup>j</sup> coefficients
  1627.                 cjCoefs[0][j] = cjsjFourier.getCj(j);
  1628.                 cjCoefs[1][j] = cjsjFourier.getdCjda(j);
  1629.                 cjCoefs[2][j] = cjsjFourier.getdCjdk(j) - (cjsjFourier.getSjLambda(j - 1) - cjsjFourier.getSjLambda(j + 1)) / 2;
  1630.                 cjCoefs[3][j] = cjsjFourier.getdCjdh(j) - (cjsjFourier.getCjLambda(j - 1) + cjsjFourier.getCjLambda(j + 1)) / 2;
  1631.                 cjCoefs[4][j] = cjsjFourier.getdCjdalpha(j);
  1632.                 cjCoefs[5][j] = cjsjFourier.getdCjdbeta(j);
  1633.                 cjCoefs[6][j] = cjsjFourier.getdCjdgamma(j);
  1634.                 cjCoefs[7][j] = cjsjFourier.getCjLambda(j);

  1636.                 //Compute the S<sup>j</sup> coefficients
  1637.                 sjCoefs[0][j] = cjsjFourier.getSj(j);
  1638.                 sjCoefs[1][j] = cjsjFourier.getdSjda(j);
  1639.                 sjCoefs[2][j] = cjsjFourier.getdSjdk(j) + (cjsjFourier.getCjLambda(j - 1) - cjsjFourier.getCjLambda(j + 1)) / 2;
  1640.                 sjCoefs[3][j] = cjsjFourier.getdSjdh(j) - (cjsjFourier.getSjLambda(j - 1) + cjsjFourier.getSjLambda(j + 1)) / 2;
  1641.                 sjCoefs[4][j] = cjsjFourier.getdSjdalpha(j);
  1642.                 sjCoefs[5][j] = cjsjFourier.getdSjdbeta(j);
  1643.                 sjCoefs[6][j] = cjsjFourier.getdSjdgamma(j);
  1644.                 sjCoefs[7][j] = cjsjFourier.getSjLambda(j);

  1646.                 //In the special case j == 1 there are some additional terms to be added
  1647.                 if (j == 1) {
  1648.                     //Additional terms for C<sup>j</sup> coefficients
  1649.                     cjCoefs[0][j] += -h * U;
  1650.                     cjCoefs[1][j] += -h * dU[0];
  1651.                     cjCoefs[2][j] += -h * dU[1];
  1652.                     cjCoefs[3][j] += -(h * dU[2] + U + cjsjFourier.getC0Lambda());
  1653.                     cjCoefs[4][j] += -h * dU[3];
  1654.                     cjCoefs[5][j] += -h * dU[4];
  1655.                     cjCoefs[6][j] += -h * dU[5];

  1657.                     //Additional terms for S<sup>j</sup> coefficients
  1658.                     sjCoefs[0][j] += k * U;
  1659.                     sjCoefs[1][j] += k * dU[0];
  1660.                     sjCoefs[2][j] += k * dU[1] + U + cjsjFourier.getC0Lambda();
  1661.                     sjCoefs[3][j] += k * dU[2];
  1662.                     sjCoefs[4][j] += k * dU[3];
  1663.                     sjCoefs[5][j] += k * dU[4];
  1664.                     sjCoefs[6][j] += k * dU[5];
  1665.                 }
  1666.             }
  1667.         }

  1669.         /** Get the coefficient C<sup>j</sup> for the function S.
  1670.          * <br>
  1671.          * Possible values for j are within the interval [0,jMax].
  1672.          * The value 0 is used to obtain the free coefficient C⁰
  1673.          * @param j j index
  1674.          * @return C<sup>j</sup> for the function S
  1675.          */
  1676.         public double getSCj(final int j) {
  1677.             return cjCoefs[0][j];
  1678.         }

  1680.         /** Get the coefficient S<sup>j</sup> for the function S.
  1681.          * <br>
  1682.          * Possible values for j are within the interval [1,jMax].
  1683.          * @param j j index
  1684.          * @return S<sup>j</sup> for the function S
  1685.          */
  1686.         public double getSSj(final int j) {
  1687.             return sjCoefs[0][j];
  1688.         }

  1690.         /** Get the coefficient C<sup>j</sup> for the derivative dS/da.
  1691.          * <br>
  1692.          * Possible values for j are within the interval [0,jMax].
  1693.          * The value 0 is used to obtain the free coefficient C⁰
  1694.          * @param j j index
  1695.          * @return C<sup>j</sup> for the function dS/da
  1696.          */
  1697.         public double getdSdaCj(final int j) {
  1698.             return cjCoefs[1][j];
  1699.         }

  1701.         /** Get the coefficient S<sup>j</sup> for the derivative dS/da.
  1702.          * <br>
  1703.          * Possible values for j are within the interval [1,jMax].
  1704.          * @param j j index
  1705.          * @return S<sup>j</sup> for the derivative dS/da
  1706.          */
  1707.         public double getdSdaSj(final int j) {
  1708.             return sjCoefs[1][j];
  1709.         }

  1711.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dk
  1712.          * <br>
  1713.          * Possible values for j are within the interval [0,jMax].
  1714.          * The value 0 is used to obtain the free coefficient C⁰
  1715.          * @param j j index
  1716.          * @return C<sup>j</sup> for the function dS/dk
  1717.          */
  1718.         public double getdSdkCj(final int j) {
  1719.             return cjCoefs[2][j];
  1720.         }

  1722.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dk.
  1723.          * <br>
  1724.          * Possible values for j are within the interval [1,jMax].
  1725.          * @param j j index
  1726.          * @return S<sup>j</sup> for the derivative dS/dk
  1727.          */
  1728.         public double getdSdkSj(final int j) {
  1729.             return sjCoefs[2][j];
  1730.         }

  1732.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dh
  1733.          * <br>
  1734.          * Possible values for j are within the interval [0,jMax].
  1735.          * The value 0 is used to obtain the free coefficient C⁰
  1736.          * @param j j index
  1737.          * @return C<sup>j</sup> for the function dS/dh
  1738.          */
  1739.         public double getdSdhCj(final int j) {
  1740.             return cjCoefs[3][j];
  1741.         }

  1743.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dh.
  1744.          * <br>
  1745.          * Possible values for j are within the interval [1,jMax].
  1746.          * @param j j index
  1747.          * @return S<sup>j</sup> for the derivative dS/dh
  1748.          */
  1749.         public double getdSdhSj(final int j) {
  1750.             return sjCoefs[3][j];
  1751.         }

  1753.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dα
  1754.          * <br>
  1755.          * Possible values for j are within the interval [0,jMax].
  1756.          * The value 0 is used to obtain the free coefficient C⁰
  1757.          * @param j j index
  1758.          * @return C<sup>j</sup> for the function dS/dα
  1759.          */
  1760.         public double getdSdalphaCj(final int j) {
  1761.             return cjCoefs[4][j];
  1762.         }

  1764.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dα.
  1765.          * <br>
  1766.          * Possible values for j are within the interval [1,jMax].
  1767.          * @param j j index
  1768.          * @return S<sup>j</sup> for the derivative dS/dα
  1769.          */
  1770.         public double getdSdalphaSj(final int j) {
  1771.             return sjCoefs[4][j];
  1772.         }

  1774.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dβ
  1775.          * <br>
  1776.          * Possible values for j are within the interval [0,jMax].
  1777.          * The value 0 is used to obtain the free coefficient C⁰
  1778.          * @param j j index
  1779.          * @return C<sup>j</sup> for the function dS/dβ
  1780.          */
  1781.         public double getdSdbetaCj(final int j) {
  1782.             return cjCoefs[5][j];
  1783.         }

  1785.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dβ.
  1786.          * <br>
  1787.          * Possible values for j are within the interval [1,jMax].
  1788.          * @param j j index
  1789.          * @return S<sup>j</sup> for the derivative dS/dβ
  1790.          */
  1791.         public double getdSdbetaSj(final int j) {
  1792.             return sjCoefs[5][j];
  1793.         }

  1795.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dγ
  1796.          * <br>
  1797.          * Possible values for j are within the interval [0,jMax].
  1798.          * The value 0 is used to obtain the free coefficient C⁰
  1799.          * @param j j index
  1800.          * @return C<sup>j</sup> for the function dS/dγ
  1801.          */
  1802.         public double getdSdgammaCj(final int j) {
  1803.             return cjCoefs[6][j];
  1804.         }

  1806.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dγ.
  1807.          * <br>
  1808.          * Possible values for j are within the interval [1,jMax].
  1809.          * @param j j index
  1810.          * @return S<sup>j</sup> for the derivative dS/dγ
  1811.          */
  1812.         public double getdSdgammaSj(final int j) {
  1813.             return sjCoefs[6][j];
  1814.         }

  1816.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dλ
  1817.          * <br>
  1818.          * Possible values for j are within the interval [0,jMax].
  1819.          * The value 0 is used to obtain the free coefficient C⁰
  1820.          * @param j j index
  1821.          * @return C<sup>j</sup> for the function dS/dλ
  1822.          */
  1823.         public double getdSdlambdaCj(final int j) {
  1824.             return cjCoefs[7][j];
  1825.         }

  1827.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dλ.
  1828.          * <br>
  1829.          * Possible values for j are within the interval [1,jMax].
  1830.          * @param j j index
  1831.          * @return S<sup>j</sup> for the derivative dS/dλ
  1832.          */
  1833.         public double getdSdlambdaSj(final int j) {
  1834.             return sjCoefs[7][j];
  1835.         }
  1836.     }

  1838.     /**
  1839.      * The coefficients used to compute the short periodic contribution for the Third body perturbation.
  1840.      * <p>
  1841.      * The short periodic contribution for the Third Body is expressed in Danielson 4.2-25.<br>
  1842.      * The coefficients C<sub>i</sub>⁰, C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup>
  1843.      * are computed by replacing the corresponding values in formula 2.5.5-10.
  1844.      * </p>
  1845.      * @author Lucian Barbulescu
  1846.      */
  1847.     private static class ThirdBodyShortPeriodicCoefficients implements ShortPeriodTerms {

  1849.         /** Serializable UID. */
  1850.         private static final long serialVersionUID = 20151119L;

  1852.         /** Maximal value for j. */
  1853.         private final int jMax;

  1855.         /** Number of points used in the interpolation process. */
  1856.         private final int interpolationPoints;

  1858.         /** Max frequency of F. */
  1859.         private final int    maxFreqF;

  1861.         /** Coefficients prefix. */
  1862.         private final String prefix;

  1864.         /** All coefficients slots. */
  1865.         private final transient TimeSpanMap<Slot> slots;

  1867.         /**
  1868.          * Standard constructor.
  1869.          *  @param interpolationPoints number of points used in the interpolation process
  1870.          * @param jMax maximal value for j
  1871.          * @param maxFreqF Max frequency of F
  1872.          * @param bodyName third body name
  1873.          * @param slots all coefficients slots
  1874.          */
  1875.         ThirdBodyShortPeriodicCoefficients(final int jMax, final int interpolationPoints,
  1876.                                            final int maxFreqF, final String bodyName,
  1877.                                            final TimeSpanMap<Slot> slots) {
  1878.             this.jMax                = jMax;
  1879.             this.interpolationPoints = interpolationPoints;
  1880.             this.maxFreqF            = maxFreqF;
  1881.             this.prefix              = "DSST-3rd-body-" + bodyName + "-";
  1882.             this.slots               = slots;
  1883.         }

  1885.         /** Get the slot valid for some date.
  1886.          * @param meanStates mean states defining the slot
  1887.          * @return slot valid at the specified date
  1888.          */
  1889.         public Slot createSlot(final SpacecraftState... meanStates) {
  1890.             final Slot         slot  = new Slot(jMax, interpolationPoints);
  1891.             final AbsoluteDate first = meanStates[0].getDate();
  1892.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  1893.             if (first.compareTo(last) <= 0) {
  1894.                 slots.addValidAfter(slot, first);
  1895.             } else {
  1896.                 slots.addValidBefore(slot, first);
  1897.             }
  1898.             return slot;
  1899.         }

  1901.         /** {@inheritDoc} */
  1902.         @Override
  1903.         public double[] value(final Orbit meanOrbit) {

  1905.             // select the coefficients slot
  1906.             final Slot slot = slots.get(meanOrbit.getDate());

  1908.             // the current eccentric longitude
  1909.             final double F = meanOrbit.getLE();

  1911.             //initialize the short periodic contribution with the corresponding C⁰ coeficient
  1912.             final double[] shortPeriodic = slot.cij[0].value(meanOrbit.getDate());

  1914.             // Add the cos and sin dependent terms
  1915.             for (int j = 1; j <= maxFreqF; j++) {
  1916.                 //compute cos and sin
  1917.                 final double cosjF = FastMath.cos(j * F);
  1918.                 final double sinjF = FastMath.sin(j * F);

  1920.                 final double[] c = slot.cij[j].value(meanOrbit.getDate());
  1921.                 final double[] s = slot.sij[j].value(meanOrbit.getDate());
  1922.                 for (int i = 0; i < 6; i++) {
  1923.                     shortPeriodic[i] += c[i] * cosjF + s[i] * sinjF;
  1924.                 }
  1925.             }

  1927.             return shortPeriodic;

  1929.         }

  1931.         /** {@inheritDoc} */
  1932.         @Override
  1933.         public String getCoefficientsKeyPrefix() {
  1934.             return prefix;
  1935.         }

  1937.         /** {@inheritDoc}
  1938.          * <p>
  1939.          * For third body attraction forces,there are maxFreqF + 1 cj coefficients,
  1940.          * maxFreqF sj coefficients where maxFreqF depends on the orbit.
  1941.          * The j index is the integer multiplier for the eccentric longitude argument
  1942.          * in the cj and sj coefficients.
  1943.          * </p>
  1944.          */
  1945.         @Override
  1946.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected)
  1947.             throws OrekitException {

  1949.             // select the coefficients slot
  1950.             final Slot slot = slots.get(date);

  1952.             final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * maxFreqF + 1);
  1953.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "c", 0);
  1954.             for (int j = 1; j <= maxFreqF; j++) {
  1955.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  1956.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  1957.             }
  1958.             return coefficients;

  1960.         }

  1962.         /** Put a coefficient in a map if selected.
  1963.          * @param map map to populate
  1964.          * @param selected set of coefficients that should be put in the map
  1965.          * (empty set means all coefficients are selected)
  1966.          * @param value coefficient value
  1967.          * @param id coefficient identifier
  1968.          * @param indices list of coefficient indices
  1969.          */
  1970.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
  1971.                                      final double[] value, final String id, final int... indices) {
  1972.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  1973.             keyBuilder.append(id);
  1974.             for (int index : indices) {
  1975.                 keyBuilder.append('[').append(index).append(']');
  1976.             }
  1977.             final String key = keyBuilder.toString();
  1978.             if (selected.isEmpty() || selected.contains(key)) {
  1979.                 map.put(key, value);
  1980.             }
  1981.         }

  1983.         /** Replace the instance with a data transfer object for serialization.
  1984.          * @return data transfer object that will be serialized
  1985.          * @exception NotSerializableException if an additional state provider is not serializable
  1986.          */
  1987.         private Object writeReplace() throws NotSerializableException {

  1989.             // slots transitions
  1990.             final SortedSet<TimeSpanMap.Transition<Slot>> transitions     = slots.getTransitions();
  1991.             final AbsoluteDate[]                          transitionDates = new AbsoluteDate[transitions.size()];
  1992.             final Slot[]                                  allSlots        = new Slot[transitions.size() + 1];
  1993.             int i = 0;
  1994.             for (final TimeSpanMap.Transition<Slot> transition : transitions) {
  1995.                 if (i == 0) {
  1996.                     // slot before the first transition
  1997.                     allSlots[i] = transition.getBefore();
  1998.                 }
  1999.                 if (i < transitionDates.length) {
  2000.                     transitionDates[i] = transition.getDate();
  2001.                     allSlots[++i]      = transition.getAfter();
  2002.                 }
  2003.             }

  2005.             return new DataTransferObject(jMax, interpolationPoints, maxFreqF, prefix,
  2006.                                           transitionDates, allSlots);

  2008.         }


  2011.         /** Internal class used only for serialization. */
  2012.         private static class DataTransferObject implements Serializable {

  2014.             /** Serializable UID. */
  2015.             private static final long serialVersionUID = 20160319L;

  2017.             /** Maximum value for j index. */
  2018.             private final int jMax;

  2020.             /** Number of points used in the interpolation process. */
  2021.             private final int interpolationPoints;

  2023.             /** Max frequency of F. */
  2024.             private final int    maxFreqF;

  2026.             /** Coefficients prefix. */
  2027.             private final String prefix;

  2029.             /** Transitions dates. */
  2030.             private final AbsoluteDate[] transitionDates;

  2032.             /** All slots. */
  2033.             private final Slot[] allSlots;

  2035.             /** Simple constructor.
  2036.              * @param jMax maximum value for j index
  2037.              * @param interpolationPoints number of points used in the interpolation process
  2038.              * @param maxFreqF max frequency of F
  2039.              * @param prefix prefix for coefficients keys
  2040.              * @param transitionDates transitions dates
  2041.              * @param allSlots all slots
  2042.              */
  2043.             DataTransferObject(final int jMax, final int interpolationPoints,
  2044.                                final int maxFreqF, final String prefix,
  2045.                                final AbsoluteDate[] transitionDates, final Slot[] allSlots) {
  2046.                 this.jMax                  = jMax;
  2047.                 this.interpolationPoints   = interpolationPoints;
  2048.                 this.maxFreqF              = maxFreqF;
  2049.                 this.prefix                = prefix;
  2050.                 this.transitionDates       = transitionDates;
  2051.                 this.allSlots              = allSlots;
  2052.             }

  2054.             /** Replace the deserialized data transfer object with a {@link ThirdBodyShortPeriodicCoefficients}.
  2055.              * @return replacement {@link ThirdBodyShortPeriodicCoefficients}
  2056.              */
  2057.             private Object readResolve() {

  2059.                 final TimeSpanMap<Slot> slots = new TimeSpanMap<Slot>(allSlots[0]);
  2060.                 for (int i = 0; i < transitionDates.length; ++i) {
  2061.                     slots.addValidAfter(allSlots[i + 1], transitionDates[i]);
  2062.                 }

  2064.                 return new ThirdBodyShortPeriodicCoefficients(jMax, interpolationPoints, maxFreqF, prefix, slots);

  2066.             }

  2068.         }

  2070.     }

  2072.     /** Coefficients valid for one time slot. */
  2073.     private static class Slot implements Serializable {

  2075.         /** Serializable UID. */
  2076.         private static final long serialVersionUID = 20160319L;

  2078.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  2079.          * <p>
  2080.          * The index order is cij[j][i] <br/>
  2081.          * i corresponds to the equinoctial element, as follows: <br/>
  2082.          * - i=0 for a <br/>
  2083.          * - i=1 for k <br/>
  2084.          * - i=2 for h <br/>
  2085.          * - i=3 for q <br/>
  2086.          * - i=4 for p <br/>
  2087.          * - i=5 for λ <br/>
  2088.          * </p>
  2089.          */
  2090.         private final ShortPeriodicsInterpolatedCoefficient[] cij;

  2092.         /** The coefficients S<sub>i</sub><sup>j</sup>.
  2093.          * <p>
  2094.          * The index order is sij[j][i] <br/>
  2095.          * i corresponds to the equinoctial element, as follows: <br/>
  2096.          * - i=0 for a <br/>
  2097.          * - i=1 for k <br/>
  2098.          * - i=2 for h <br/>
  2099.          * - i=3 for q <br/>
  2100.          * - i=4 for p <br/>
  2101.          * - i=5 for λ <br/>
  2102.          * </p>
  2103.          */
  2104.         private final ShortPeriodicsInterpolatedCoefficient[] sij;

  2106.         /** Simple constructor.
  2107.          *  @param jMax maximum value for j index
  2108.          *  @param interpolationPoints number of points used in the interpolation process
  2109.          */
  2110.         Slot(final int jMax, final int interpolationPoints) {
  2111.             // allocate the coefficients arrays
  2112.             cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
  2113.             sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
  2114.             for (int j = 0; j <= jMax; j++) {
  2115.                 cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  2116.                 sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  2117.             }


  2120.         }
  2121.     }

  2123. }
Created with JaCoCo 0.7.9.201702052155