DSSTThirdBody.java

  1. /* Copyright 2002-2023 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.propagation.semianalytical.dsst.forces;

  19. import java.lang.reflect.Array;
  20. import java.util.ArrayList;
  21. import java.util.Arrays;
  22. import java.util.Collections;
  23. import java.util.HashMap;
  24. import java.util.List;
  25. import java.util.Map;
  26. import java.util.Set;
  27. import java.util.SortedMap;

  29. import org.hipparchus.CalculusFieldElement;
  30. import org.hipparchus.Field;
  31. import org.hipparchus.analysis.differentiation.FieldGradient;
  32. import org.hipparchus.util.CombinatoricsUtils;
  33. import org.hipparchus.util.FastMath;
  34. import org.hipparchus.util.FieldSinCos;
  35. import org.hipparchus.util.MathArrays;
  36. import org.hipparchus.util.SinCos;
  37. import org.orekit.attitudes.AttitudeProvider;
  38. import org.orekit.bodies.CelestialBodies;
  39. import org.orekit.bodies.CelestialBody;
  40. import org.orekit.orbits.FieldOrbit;
  41. import org.orekit.orbits.Orbit;
  42. import org.orekit.propagation.FieldSpacecraftState;
  43. import org.orekit.propagation.PropagationType;
  44. import org.orekit.propagation.SpacecraftState;
  45. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  46. import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
  47. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
  48. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory.NSKey;
  49. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  50. import org.orekit.propagation.semianalytical.dsst.utilities.FieldCjSjCoefficient;
  51. import org.orekit.propagation.semianalytical.dsst.utilities.FieldShortPeriodicsInterpolatedCoefficient;
  52. import org.orekit.propagation.semianalytical.dsst.utilities.JacobiPolynomials;
  53. import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
  54. import org.orekit.propagation.semianalytical.dsst.utilities.hansen.FieldHansenThirdBodyLinear;
  55. import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenThirdBodyLinear;
  56. import org.orekit.time.AbsoluteDate;
  57. import org.orekit.time.FieldAbsoluteDate;
  58. import org.orekit.utils.FieldTimeSpanMap;
  59. import org.orekit.utils.ParameterDriver;
  60. import org.orekit.utils.TimeSpanMap;

  62. /** Third body attraction perturbation to the
  63.  *  {@link org.orekit.propagation.semianalytical.dsst.DSSTPropagator DSSTPropagator}.
  64.  *
  65.  *  @author Romain Di Costanzo
  66.  *  @author Pascal Parraud
  67.  *  @author Bryan Cazabonne (field translation)
  68.  */
  69. public class DSSTThirdBody implements DSSTForceModel {
  70.     /**  Name of the prefix for short period coefficients keys. */
  71.     public static final String SHORT_PERIOD_PREFIX = "DSST-3rd-body-";

  73.     /** Name of the single parameter of this model: the attraction coefficient. */
  74.     public static final String ATTRACTION_COEFFICIENT = " attraction coefficient";

  76.     /** Max power for summation. */
  77.     public static final int MAX_POWER = 22;

  79.     /** Truncation tolerance for big, eccentric orbits. */
  80.     public static final double BIG_TRUNCATION_TOLERANCE = 1.e-1;

  82.     /** Truncation tolerance for small orbits. */
  83.     public static final double SMALL_TRUNCATION_TOLERANCE = 1.9e-6;

  85.     /** Central attraction scaling factor.
  86.      * <p>
  87.      * We use a power of 2 to avoid numeric noise introduction
  88.      * in the multiplications/divisions sequences.
  89.      * </p>
  90.      */
  91.     private static final double MU_SCALE = FastMath.scalb(1.0, 32);

  93.     /** Retrograde factor I.
  94.      *  <p>
  95.      *  DSST model needs equinoctial orbit as internal representation.
  96.      *  Classical equinoctial elements have discontinuities when inclination
  97.      *  is close to zero. In this representation, I = +1. <br>
  98.      *  To avoid this discontinuity, another representation exists and equinoctial
  99.      *  elements can be expressed in a different way, called "retrograde" orbit.
  100.      *  This implies I = -1. <br>
  101.      *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  102.      *  But for the sake of consistency with the theory, the retrograde factor
  103.      *  has been kept in the formulas.
  104.      *  </p>
  105.      */
  106.     private static final int    I = 1;

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

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

  114.     /** V<sub>ns</sub> coefficients. */
  115.     private final SortedMap<NSKey, Double> Vns;

  117.     /** Force model static context. Initialized with short period terms. */
  118.     private DSSTThirdBodyStaticContext staticContext;

  120.     /** If false, the static context needs to be initialized. */
  121.     private boolean doesStaticContextNeedsInitialization;

  123.     /** The 3rd body to consider. */
  124.     private final CelestialBody    body;

  126.     /** Short period terms. */
  127.     private ThirdBodyShortPeriodicCoefficients shortPeriods;

  129.     /** "Field" Short period terms. */
  130.     private Map<Field<?>, FieldThirdBodyShortPeriodicCoefficients<?>> fieldShortPeriods;

  132.     /** Drivers for third body attraction coefficient and gravitational parameter. */
  133.     private final List<ParameterDriver> parameterDrivers;

  135.     /** Hansen objects. */
  136.     private HansenObjects hansen;

  138.     /** Hansen objects for field elements. */
  139.     private Map<Field<?>, FieldHansenObjects<?>> fieldHansen;

  141.     /** Complete constructor.
  142.      *  @param body the 3rd body to consider
  143.      *  @param mu central attraction coefficient
  144.      *            (<b>i.e., attraction coefficient of the central body, not the one of the 3rd body</b>)
  145.      *  @see CelestialBodies
  146.      */
  147.     public DSSTThirdBody(final CelestialBody body, final double mu) {
  148.         parameterDrivers = new ArrayList<>(2);
  149.         parameterDrivers.add(new ParameterDriver(body.getName() + DSSTThirdBody.ATTRACTION_COEFFICIENT,
  150.                                                  body.getGM(), MU_SCALE,
  151.                                                  0.0, Double.POSITIVE_INFINITY));
  152.         parameterDrivers.add(new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
  153.                                                  mu, MU_SCALE,
  154.                                                  0.0, Double.POSITIVE_INFINITY));

  156.         this.body = body;
  157.         this.Vns  = CoefficientsFactory.computeVns(MAX_POWER);
  158.         this.doesStaticContextNeedsInitialization = true;

  160.         fieldShortPeriods  = new HashMap<>();
  161.         fieldHansen        = new HashMap<>();
  162.     }

  164.     /** Get third body.
  165.      *  @return third body
  166.      */
  167.     public CelestialBody getBody() {
  168.         return body;
  169.     }

  171.     /** Initializes the static 3rd body context if needed.
  172.      * @param auxiliaryElements auxiliary elements
  173.      * @param x DSST Chi element
  174.      * @param r3 distance from center of mass of the central body to the 3rd body
  175.      * @param parameters force model parameters
  176.      */
  177.     private void initializeStaticContextIfNeeded(final AuxiliaryElements auxiliaryElements,
  178.                                                  final double x, final double r3,
  179.                                                  final double[] parameters) {
  180.         if (doesStaticContextNeedsInitialization) {
  181.             staticContext = new DSSTThirdBodyStaticContext(auxiliaryElements, x, r3, parameters);
  182.             doesStaticContextNeedsInitialization = false;
  183.         }
  184.     }

  186.     /** Computes the highest power of the eccentricity and the highest power
  187.      *  of a/R3 to appear in the truncated analytical power series expansion.
  188.      *  <p>
  189.      *  This method computes the upper value for the 3rd body potential and
  190.      *  determines the maximal powers for the eccentricity and a/R3 producing
  191.      *  potential terms bigger than a defined tolerance.
  192.      *  </p>
  193.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  194.      *  @param type type of the elements used during the propagation
  195.      *  @param parameters values of the force model parameters for state date (1 value for each parameters)
  196.      */
  197.     @Override
  198.     public List<ShortPeriodTerms> initializeShortPeriodTerms(final AuxiliaryElements auxiliaryElements,
  199.                                                              final PropagationType type,
  200.                                                              final double[] parameters) {

  202.         // Initializes specific parameters.
  203.         final DSSTThirdBodyDynamicContext context = initializeStep(auxiliaryElements, parameters);

  205.         // Static context
  206.         initializeStaticContextIfNeeded(auxiliaryElements, context.getX(), context.getR3(), parameters);

  208.         // Hansen objects
  209.         hansen = new HansenObjects();

  211.         // Initialize short period terms
  212.         final int jMax = staticContext.getMaxAR3Pow() + 1;
  213.         shortPeriods = new ThirdBodyShortPeriodicCoefficients(jMax, INTERPOLATION_POINTS,
  214.                                                               staticContext.getMaxFreqF(), body.getName(),
  215.                                                               new TimeSpanMap<Slot>(new Slot(jMax, INTERPOLATION_POINTS)));

  217.         final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
  218.         list.add(shortPeriods);
  219.         return list;

  221.     }

  223.     /** {@inheritDoc} */
  224.     @Override
  225.     public <T extends CalculusFieldElement<T>> List<FieldShortPeriodTerms<T>> initializeShortPeriodTerms(final FieldAuxiliaryElements<T> auxiliaryElements,
  226.                                                                                                          final PropagationType type,
  227.                                                                                                          final T[] parameters) {

  229.         // Field used by default
  230.         final Field<T> field = auxiliaryElements.getDate().getField();

  232.         // Initializes specific parameters.
  233.         final FieldDSSTThirdBodyDynamicContext<T> context = initializeStep(auxiliaryElements, parameters);

  235.         // Static context (only provide integers. So, derivatives are not taken into account and getReal() method is accepted)
  236.         initializeStaticContextIfNeeded(auxiliaryElements.toAuxiliaryElements(), context.getX().getReal(), context.getR3().getReal(), getParameters());

  238.         // Hansen object
  239.         fieldHansen.put(field, new FieldHansenObjects<>(field));

  241.         // Initialize short period terms
  242.         final int jMax = staticContext.getMaxAR3Pow() + 1;
  243.         final FieldThirdBodyShortPeriodicCoefficients<T> ftbspc =
  244.                         new FieldThirdBodyShortPeriodicCoefficients<>(jMax, INTERPOLATION_POINTS,
  245.                                         staticContext.getMaxFreqF(), body.getName(),
  246.                                                                       new FieldTimeSpanMap<>(new FieldSlot<>(jMax,
  247.                                                                                                              INTERPOLATION_POINTS),
  248.                                                                                              field));
  249.         fieldShortPeriods.put(field, ftbspc);
  250.         return Collections.singletonList(ftbspc);
  251.     }

  253.     /** Performs initialization at each integration step for the current force model.
  254.      *  <p>
  255.      *  This method aims at being called before mean elements rates computation.
  256.      *  </p>
  257.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  258.      *  @param parameters values of the force model parameters  (only 1 values for each parameters corresponding
  259.      * to state date) obtained by calling the extract parameter method {@link #extractParameters(double[], AbsoluteDate)}
  260.      * to selected the right value for state date or by getting the parameters for a specific date
  261.      * {@link #getParameters(AbsoluteDate)}.
  262.      *  @return new force model context
  263.      */
  264.     private DSSTThirdBodyDynamicContext initializeStep(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
  265.         return new DSSTThirdBodyDynamicContext(auxiliaryElements, body, parameters);
  266.     }

  268.     /** Performs initialization at each integration step for the current force model.
  269.      *  <p>
  270.      *  This method aims at being called before mean elements rates computation.
  271.      *  </p>
  272.      *  @param <T> type of the elements
  273.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  274.      *  @param parameters values of the force model parameters at state date (1 value per parameter driver)
  275.      *  @return new force model context
  276.      */
  277.     private <T extends CalculusFieldElement<T>> FieldDSSTThirdBodyDynamicContext<T> initializeStep(final FieldAuxiliaryElements<T> auxiliaryElements,
  278.                                                                                         final T[] parameters) {
  279.         return new FieldDSSTThirdBodyDynamicContext<>(auxiliaryElements, body, parameters);
  280.     }

  282.     /** {@inheritDoc} */
  283.     @Override
  284.     public double[] getMeanElementRate(final SpacecraftState currentState,
  285.                                        final AuxiliaryElements auxiliaryElements, final double[] parameters) {


  288.         // Container for attributes
  289.         final DSSTThirdBodyDynamicContext context = initializeStep(auxiliaryElements, parameters);

  291.         // a / R3 up to power maxAR3Pow
  292.         final double[] aoR3Pow = computeAoR3Pow(context);

  294.         // Qns coefficients
  295.         final double[][] Qns = CoefficientsFactory.computeQns(context.getGamma(), staticContext.getMaxAR3Pow(), FastMath.max(staticContext.getMaxEccPow(), MAX_ECCPOWER_SP));

  297.         // Access to potential U derivatives
  298.         final UAnddU udu = new UAnddU(context, hansen, aoR3Pow, Qns);

  300.         // Compute cross derivatives [Eq. 2.2-(8)]
  301.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  302.         final double UAlphaGamma   = context.getAlpha() * udu.getdUdGa() - context.getGamma() * udu.getdUdAl();
  303.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  304.         final double UBetaGamma    =  context.getBeta() * udu.getdUdGa() - context.getGamma() * udu.getdUdBe();
  305.         // Common factor
  306.         final double pUAGmIqUBGoAB = (auxiliaryElements.getP() * UAlphaGamma - I * auxiliaryElements.getQ() * UBetaGamma) * context.getOoAB();

  308.         // Compute mean elements rates [Eq. 3.1-(1)]
  309.         final double da =  0.;
  310.         final double dh =  context.getBoA() * udu.getdUdk() + auxiliaryElements.getK() * pUAGmIqUBGoAB;
  311.         final double dk = -context.getBoA() * udu.getdUdh() - auxiliaryElements.getH() * pUAGmIqUBGoAB;
  312.         final double dp =  context.getMCo2AB() * UBetaGamma;
  313.         final double dq =  context.getMCo2AB() * UAlphaGamma * I;
  314.         final double dM =  context.getM2aoA() * udu.getdUda() + context.getBoABpo() * (auxiliaryElements.getH() * udu.getdUdh() + auxiliaryElements.getK() * udu.getdUdk()) + pUAGmIqUBGoAB;

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

  318.     }

  320.     /** {@inheritDoc} */
  321.     @Override
  322.     public <T extends CalculusFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> currentState,
  323.                                                                   final FieldAuxiliaryElements<T> auxiliaryElements,
  324.                                                                   final T[] parameters) {

  326.         // Parameters for array building
  327.         final Field<T> field = currentState.getDate().getField();
  328.         final T        zero  = field.getZero();

  330.         // Container for attributes
  331.         final FieldDSSTThirdBodyDynamicContext<T> context = initializeStep(auxiliaryElements, parameters);

  333.         // a / R3 up to power maxAR3Pow
  334.         final T[] aoR3Pow = computeAoR3Pow(context);

  336.         // Qns coefficients
  337.         final T[][] Qns = CoefficientsFactory.computeQns(context.getGamma(), staticContext.getMaxAR3Pow(), FastMath.max(staticContext.getMaxEccPow(), MAX_ECCPOWER_SP));

  339.         // Hansen objects
  340.         @SuppressWarnings("unchecked")
  341.         final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);

  343.         // Access to potential U derivatives
  344.         final FieldUAnddU<T> udu = new FieldUAnddU<>(context, fho, aoR3Pow, Qns);

  346.         // Compute cross derivatives [Eq. 2.2-(8)]
  347.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  348.         final T UAlphaGamma   = udu.getdUdGa().multiply(context.getAlpha()).subtract(udu.getdUdAl().multiply(context.getGamma()));
  349.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  350.         final T UBetaGamma    = udu.getdUdGa().multiply(context.getBeta()).subtract(udu.getdUdBe().multiply(context.getGamma()));
  351.         // Common factor
  352.         final T pUAGmIqUBGoAB = (UAlphaGamma.multiply(auxiliaryElements.getP()).subtract(UBetaGamma.multiply(auxiliaryElements.getQ()).multiply(I))).multiply(context.getOoAB());

  354.         // Compute mean elements rates [Eq. 3.1-(1)]
  355.         final T da =  zero;
  356.         final T dh =  udu.getdUdk().multiply(context.getBoA()).add(pUAGmIqUBGoAB.multiply(auxiliaryElements.getK()));
  357.         final T dk =  ((udu.getdUdh().multiply(context.getBoA())).negate()).subtract(pUAGmIqUBGoAB.multiply(auxiliaryElements.getH()));
  358.         final T dp =  UBetaGamma.multiply(context.getMCo2AB());
  359.         final T dq =  UAlphaGamma.multiply(I).multiply(context.getMCo2AB());
  360.         final T dM =  pUAGmIqUBGoAB.add(udu.getdUda().multiply(context.getM2aoA())).add((udu.getdUdh().multiply(auxiliaryElements.getH()).add(udu.getdUdk().multiply(auxiliaryElements.getK()))).multiply(context.getBoABpo()));

  362.         final T[] elements = MathArrays.buildArray(field, 6);
  363.         elements[0] = da;
  364.         elements[1] = dk;
  365.         elements[2] = dh;
  366.         elements[3] = dq;
  367.         elements[4] = dp;
  368.         elements[5] = dM;

  370.         return elements;

  372.     }

  374.     /** {@inheritDoc} */
  375.     @Override
  376.     public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {

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

  380.         for (final SpacecraftState meanState : meanStates) {

  382.             // Auxiliary elements related to the current orbit
  383.             final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(meanState.getOrbit(), I);

  385.             // Extract the proper parameters valid for the corresponding meanState date from the input array
  386.             final double[] extractedParameters = this.extractParameters(parameters, auxiliaryElements.getDate());

  388.             // Container of attributes
  389.             final DSSTThirdBodyDynamicContext context = initializeStep(auxiliaryElements, extractedParameters);

  391.             // a / R3 up to power maxAR3Pow
  392.             final double[] aoR3Pow = computeAoR3Pow(context);

  394.             // Qns coefficients
  395.             final double[][] Qns = CoefficientsFactory.computeQns(context.getGamma(), staticContext.getMaxAR3Pow(), FastMath.max(staticContext.getMaxEccPow(), MAX_ECCPOWER_SP));

  397.             final GeneratingFunctionCoefficients gfCoefs =
  398.                             new GeneratingFunctionCoefficients(staticContext.getMaxAR3Pow(), MAX_ECCPOWER_SP, staticContext.getMaxAR3Pow() + 1, context, hansen, aoR3Pow, Qns);

  400.             //Compute additional quantities
  401.             // 2 * a / An
  402.             final double ax2oAn  = -context.getM2aoA() / context.getMeanMotion();
  403.             // B / An
  404.             final double BoAn    = context.getBoA() / context.getMeanMotion();
  405.             // 1 / ABn
  406.             final double ooABn   = context.getOoAB() / context.getMeanMotion();
  407.             // C / 2ABn
  408.             final double Co2ABn  = -context.getMCo2AB() / context.getMeanMotion();
  409.             // B / (A * (1 + B) * n)
  410.             final double BoABpon = context.getBoABpo() / context.getMeanMotion();
  411.             // -3 / n²a² = -3 / nA
  412.             final double m3onA   = -3 / (context.getA() * context.getMeanMotion());

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

  419.                 // Compute the cross derivatives operator :
  420.                 final double SAlphaGammaCj    = context.getAlpha() * gfCoefs.getdSdgammaCj(j) - context.getGamma() * gfCoefs.getdSdalphaCj(j);
  421.                 final double SAlphaBetaCj     = context.getAlpha() * gfCoefs.getdSdbetaCj(j)  - context.getBeta()  * gfCoefs.getdSdalphaCj(j);
  422.                 final double SBetaGammaCj     = context.getBeta() * gfCoefs.getdSdgammaCj(j) - context.getGamma() * gfCoefs.getdSdbetaCj(j);
  423.                 final double ShkCj            = auxiliaryElements.getH() * gfCoefs.getdSdkCj(j)     -  auxiliaryElements.getK()    * gfCoefs.getdSdhCj(j);
  424.                 final double pSagmIqSbgoABnCj = (auxiliaryElements.getP() * SAlphaGammaCj - I * auxiliaryElements.getQ() * SBetaGammaCj) * ooABn;
  425.                 final double ShkmSabmdSdlCj   = ShkCj - SAlphaBetaCj - gfCoefs.getdSdlambdaCj(j);

  427.                 currentCij[0] =  ax2oAn * gfCoefs.getdSdlambdaCj(j);
  428.                 currentCij[1] =  -(BoAn * gfCoefs.getdSdhCj(j) + auxiliaryElements.getH() * pSagmIqSbgoABnCj + auxiliaryElements.getK() * BoABpon * gfCoefs.getdSdlambdaCj(j));
  429.                 currentCij[2] =    BoAn * gfCoefs.getdSdkCj(j) + auxiliaryElements.getK() * pSagmIqSbgoABnCj - auxiliaryElements.getH() * BoABpon * gfCoefs.getdSdlambdaCj(j);
  430.                 currentCij[3] =  Co2ABn * (auxiliaryElements.getQ() * ShkmSabmdSdlCj - I * SAlphaGammaCj);
  431.                 currentCij[4] =  Co2ABn * (auxiliaryElements.getP() * ShkmSabmdSdlCj - SBetaGammaCj);
  432.                 currentCij[5] = -ax2oAn * gfCoefs.getdSdaCj(j) + BoABpon * (auxiliaryElements.getH() * gfCoefs.getdSdhCj(j) + auxiliaryElements.getK() * gfCoefs.getdSdkCj(j)) + pSagmIqSbgoABnCj + m3onA * gfCoefs.getSCj(j);

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

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

  440.                 // Compute the cross derivatives operator :
  441.                 final double SAlphaGammaSj    = context.getAlpha() * gfCoefs.getdSdgammaSj(j) - context.getGamma() * gfCoefs.getdSdalphaSj(j);
  442.                 final double SAlphaBetaSj     = context.getAlpha() * gfCoefs.getdSdbetaSj(j)  - context.getBeta()  * gfCoefs.getdSdalphaSj(j);
  443.                 final double SBetaGammaSj     =  context.getBeta() * gfCoefs.getdSdgammaSj(j) - context.getGamma() * gfCoefs.getdSdbetaSj(j);
  444.                 final double ShkSj            =     auxiliaryElements.getH() * gfCoefs.getdSdkSj(j)     -  auxiliaryElements.getK()    * gfCoefs.getdSdhSj(j);
  445.                 final double pSagmIqSbgoABnSj = (auxiliaryElements.getP() * SAlphaGammaSj - I * auxiliaryElements.getQ() * SBetaGammaSj) * ooABn;
  446.                 final double ShkmSabmdSdlSj   =  ShkSj - SAlphaBetaSj - gfCoefs.getdSdlambdaSj(j);

  448.                 currentSij[0] =  ax2oAn * gfCoefs.getdSdlambdaSj(j);
  449.                 currentSij[1] =  -(BoAn * gfCoefs.getdSdhSj(j) + auxiliaryElements.getH() * pSagmIqSbgoABnSj + auxiliaryElements.getK() * BoABpon * gfCoefs.getdSdlambdaSj(j));
  450.                 currentSij[2] =    BoAn * gfCoefs.getdSdkSj(j) + auxiliaryElements.getK() * pSagmIqSbgoABnSj - auxiliaryElements.getH() * BoABpon * gfCoefs.getdSdlambdaSj(j);
  451.                 currentSij[3] =  Co2ABn * (auxiliaryElements.getQ() * ShkmSabmdSdlSj - I * SAlphaGammaSj);
  452.                 currentSij[4] =  Co2ABn * (auxiliaryElements.getP() * ShkmSabmdSdlSj - SBetaGammaSj);
  453.                 currentSij[5] = -ax2oAn * gfCoefs.getdSdaSj(j) + BoABpon * (auxiliaryElements.getH() * gfCoefs.getdSdhSj(j) + auxiliaryElements.getK() * gfCoefs.getdSdkSj(j)) + pSagmIqSbgoABnSj + m3onA * gfCoefs.getSSj(j);

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

  458.                 if (j == 1) {
  459.                     //Compute the C⁰ coefficients using Danielson 2.5.2-15a.
  460.                     final double[] value = new double[6];
  461.                     for (int i = 0; i < 6; ++i) {
  462.                         value[i] = currentCij[i] * auxiliaryElements.getK() / 2. + currentSij[i] * auxiliaryElements.getH() / 2.;
  463.                     }
  464.                     slot.cij[0].addGridPoint(meanState.getDate(), value);
  465.                 }
  466.             }
  467.         }
  468.     }

  470.     /** {@inheritDoc} */
  471.     @Override
  472.     @SuppressWarnings("unchecked")
  473.     public <T extends CalculusFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
  474.                                                                        final FieldSpacecraftState<T>... meanStates) {

  476.         // Field used by default
  477.         final Field<T> field = meanStates[0].getDate().getField();

  479.         final FieldThirdBodyShortPeriodicCoefficients<T> ftbspc = (FieldThirdBodyShortPeriodicCoefficients<T>) fieldShortPeriods.get(field);
  480.         final FieldSlot<T> slot = ftbspc.createSlot(meanStates);
  481.         for (final FieldSpacecraftState<T> meanState : meanStates) {

  483.             // Auxiliary elements related to the current orbit
  484.             final FieldAuxiliaryElements<T> auxiliaryElements = new FieldAuxiliaryElements<>(meanState.getOrbit(), I);

  486.             // Extract the proper parameters valid for the corresponding meanState date from the input array
  487.             final T[] extractedParameters = this.extractParameters(parameters, auxiliaryElements.getDate());

  489.             // Container of attributes
  490.             final FieldDSSTThirdBodyDynamicContext<T> context = initializeStep(auxiliaryElements, extractedParameters);

  492.             // a / R3 up to power maxAR3Pow
  493.             final T[] aoR3Pow = computeAoR3Pow(context);

  495.             // Qns coefficients
  496.             final T[][] Qns = CoefficientsFactory.computeQns(context.getGamma(), staticContext.getMaxAR3Pow(), FastMath.max(staticContext.getMaxEccPow(), MAX_ECCPOWER_SP));

  498.             final FieldGeneratingFunctionCoefficients<T> gfCoefs =
  499.                             new FieldGeneratingFunctionCoefficients<>(staticContext.getMaxAR3Pow(), MAX_ECCPOWER_SP, staticContext.getMaxAR3Pow() + 1,
  500.                                                                       context, (FieldHansenObjects<T>) fieldHansen.get(field), field, aoR3Pow, Qns);

  502.             //Compute additional quantities
  503.             // 2 * a / An
  504.             final T ax2oAn  = context.getM2aoA().negate().divide(context.getMeanMotion());
  505.             // B / An
  506.             final T BoAn    = context.getBoA().divide(context.getMeanMotion());
  507.             // 1 / ABn
  508.             final T ooABn   = context.getOoAB().divide(context.getMeanMotion());
  509.             // C / 2ABn
  510.             final T Co2ABn  = context.getMCo2AB().negate().divide(context.getMeanMotion());
  511.             // B / (A * (1 + B) * n)
  512.             final T BoABpon = context.getBoABpo().divide(context.getMeanMotion());
  513.             // -3 / n²a² = -3 / nA
  514.             final T m3onA   = context.getA().multiply(context.getMeanMotion()).divide(-3.).reciprocal();

  516.             //Compute the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients.
  517.             for (int j = 1; j < slot.cij.length; j++) {
  518.                 // First compute the C<sub>i</sub><sup>j</sup> coefficients
  519.                 final T[] currentCij = MathArrays.buildArray(field, 6);

  521.                 // Compute the cross derivatives operator :
  522.                 final T SAlphaGammaCj    = context.getAlpha().multiply(gfCoefs.getdSdgammaCj(j)).subtract(context.getGamma().multiply(gfCoefs.getdSdalphaCj(j)));
  523.                 final T SAlphaBetaCj     = context.getAlpha().multiply(gfCoefs.getdSdbetaCj(j)).subtract(context.getBeta().multiply(gfCoefs.getdSdalphaCj(j)));
  524.                 final T SBetaGammaCj     = context.getBeta().multiply(gfCoefs.getdSdgammaCj(j)).subtract(context.getGamma().multiply(gfCoefs.getdSdbetaCj(j)));
  525.                 final T ShkCj            = auxiliaryElements.getH().multiply(gfCoefs.getdSdkCj(j)).subtract(auxiliaryElements.getK().multiply(gfCoefs.getdSdhCj(j)));
  526.                 final T pSagmIqSbgoABnCj = ooABn.multiply(auxiliaryElements.getP().multiply(SAlphaGammaCj).subtract(auxiliaryElements.getQ().multiply(SBetaGammaCj).multiply(I)));
  527.                 final T ShkmSabmdSdlCj   = ShkCj.subtract(SAlphaBetaCj).subtract(gfCoefs.getdSdlambdaCj(j));

  529.                 currentCij[0] = ax2oAn.multiply(gfCoefs.getdSdlambdaCj(j));
  530.                 currentCij[1] = BoAn.multiply(gfCoefs.getdSdhCj(j)).add(auxiliaryElements.getH().multiply(pSagmIqSbgoABnCj)).add(auxiliaryElements.getK().multiply(BoABpon).multiply(gfCoefs.getdSdlambdaCj(j))).negate();
  531.                 currentCij[2] = BoAn.multiply(gfCoefs.getdSdkCj(j)).add(auxiliaryElements.getK().multiply(pSagmIqSbgoABnCj)).subtract(auxiliaryElements.getH().multiply(BoABpon).multiply(gfCoefs.getdSdlambdaCj(j)));
  532.                 currentCij[3] = Co2ABn.multiply(auxiliaryElements.getQ().multiply(ShkmSabmdSdlCj).subtract(SAlphaGammaCj.multiply(I)));
  533.                 currentCij[4] = Co2ABn.multiply(auxiliaryElements.getP().multiply(ShkmSabmdSdlCj).subtract(SBetaGammaCj));
  534.                 currentCij[5] = ax2oAn.negate().multiply(gfCoefs.getdSdaCj(j)).add(BoABpon.multiply(auxiliaryElements.getH().multiply(gfCoefs.getdSdhCj(j)).add(auxiliaryElements.getK().multiply(gfCoefs.getdSdkCj(j))))).add(pSagmIqSbgoABnCj).add(m3onA.multiply(gfCoefs.getSCj(j)));

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

  539.                 // Compute the S<sub>i</sub><sup>j</sup> coefficients
  540.                 final T[] currentSij = MathArrays.buildArray(field, 6);

  542.                 // Compute the cross derivatives operator :
  543.                 final T SAlphaGammaSj    = context.getAlpha().multiply(gfCoefs.getdSdgammaSj(j)).subtract(context.getGamma().multiply(gfCoefs.getdSdalphaSj(j)));
  544.                 final T SAlphaBetaSj     = context.getAlpha().multiply(gfCoefs.getdSdbetaSj(j)).subtract(context.getBeta().multiply(gfCoefs.getdSdalphaSj(j)));
  545.                 final T SBetaGammaSj     = context.getBeta().multiply(gfCoefs.getdSdgammaSj(j)).subtract(context.getGamma().multiply(gfCoefs.getdSdbetaSj(j)));
  546.                 final T ShkSj            = auxiliaryElements.getH().multiply(gfCoefs.getdSdkSj(j)).subtract(auxiliaryElements.getK().multiply(gfCoefs.getdSdhSj(j)));
  547.                 final T pSagmIqSbgoABnSj = ooABn.multiply(auxiliaryElements.getP().multiply(SAlphaGammaSj).subtract(auxiliaryElements.getQ().multiply(SBetaGammaSj).multiply(I)));
  548.                 final T ShkmSabmdSdlSj   = ShkSj.subtract(SAlphaBetaSj).subtract(gfCoefs.getdSdlambdaSj(j));

  550.                 currentSij[0] = ax2oAn.multiply(gfCoefs.getdSdlambdaSj(j));
  551.                 currentSij[1] = BoAn.multiply(gfCoefs.getdSdhSj(j)).add(auxiliaryElements.getH().multiply(pSagmIqSbgoABnSj)).add(auxiliaryElements.getK().multiply(BoABpon).multiply(gfCoefs.getdSdlambdaSj(j))).negate();
  552.                 currentSij[2] = BoAn.multiply(gfCoefs.getdSdkSj(j)).add(auxiliaryElements.getK().multiply(pSagmIqSbgoABnSj)).subtract(auxiliaryElements.getH().multiply(BoABpon).multiply(gfCoefs.getdSdlambdaSj(j)));
  553.                 currentSij[3] = Co2ABn.multiply(auxiliaryElements.getQ().multiply(ShkmSabmdSdlSj).subtract(SAlphaGammaSj.multiply(I)));
  554.                 currentSij[4] = Co2ABn.multiply(auxiliaryElements.getP().multiply(ShkmSabmdSdlSj).subtract(SBetaGammaSj));
  555.                 currentSij[5] = ax2oAn.negate().multiply(gfCoefs.getdSdaSj(j)).add(BoABpon.multiply(auxiliaryElements.getH().multiply(gfCoefs.getdSdhSj(j)).add(auxiliaryElements.getK().multiply(gfCoefs.getdSdkSj(j))))).add(pSagmIqSbgoABnSj).add(m3onA.multiply(gfCoefs.getSSj(j)));

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

  560.                 if (j == 1) {
  561.                     //Compute the C⁰ coefficients using Danielson 2.5.2-15a.
  562.                     final T[] value = MathArrays.buildArray(field, 6);
  563.                     for (int i = 0; i < 6; ++i) {
  564.                         value[i] = currentCij[i].multiply(auxiliaryElements.getK()).divide(2.).add(currentSij[i].multiply(auxiliaryElements.getH()).divide(2.));
  565.                     }
  566.                     slot.cij[0].addGridPoint(meanState.getDate(), value);
  567.                 }
  568.             }
  569.         }
  570.     }

  572.     /**
  573.      * Computes a / R3 to the power maxAR3Pow.
  574.      * @param context force model dynamic context
  575.      * @return aoR3Pow array
  576.      */
  577.     private double[] computeAoR3Pow(final DSSTThirdBodyDynamicContext context) {
  578.         // a / R3 up to power maxAR3Pow
  579.         final double aoR3 = context.getAuxiliaryElements().getSma() / context.getR3();
  580.         final double[] aoR3Pow = new double[staticContext.getMaxAR3Pow() + 1];
  581.         aoR3Pow[0] = 1.;
  582.         for (int i = 1; i <= staticContext.getMaxAR3Pow(); i++) {
  583.             aoR3Pow[i] = aoR3 * aoR3Pow[i - 1];
  584.         }
  585.         return aoR3Pow;
  586.     }

  588.     /**
  589.      * Computes a / R3 to the power maxAR3Pow.
  590.      * @param context force model dynamic context
  591.      * @param <T> type of the elements
  592.      * @return aoR3Pow array
  593.      */
  594.     private <T extends CalculusFieldElement<T>> T[] computeAoR3Pow(final FieldDSSTThirdBodyDynamicContext<T> context) {
  595.         final Field<T> field = context.getA().getField();
  596.         // a / R3 up to power maxAR3Pow
  597.         final T aoR3 = context.getFieldAuxiliaryElements().getSma().divide(context.getR3());
  598.         final T[] aoR3Pow = MathArrays.buildArray(field, staticContext.getMaxAR3Pow() + 1);
  599.         aoR3Pow[0] = field.getOne();
  600.         for (int i = 1; i <= staticContext.getMaxAR3Pow(); i++) {
  601.             aoR3Pow[i] = aoR3.multiply(aoR3Pow[i - 1]);
  602.         }
  603.         return aoR3Pow;
  604.     }

  606.     /** {@inheritDoc} */
  607.     @Override
  608.     public void registerAttitudeProvider(final AttitudeProvider provider) {
  609.         //nothing is done since this contribution is not sensitive to attitude
  610.     }

  612.     /** {@inheritDoc} */
  613.     @Override
  614.     public List<ParameterDriver> getParametersDrivers() {
  615.         return Collections.unmodifiableList(parameterDrivers);
  616.     }

  618.     /** Computes the C<sup>j</sup> and S<sup>j</sup> coefficients Danielson 4.2-(15,16)
  619.      * and their derivatives.
  620.      *  <p>
  621.      *  CS Mathematical Report $3.5.3.2
  622.      *  </p>
  623.      */
  624.     private class FourierCjSjCoefficients {

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

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

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

  635.         /** The Fourier coefficients C<sup>j</sup> and their derivatives.
  636.          * <p>
  637.          * Each column of the matrix contains the following values: <br/>
  638.          * - C<sup>j</sup> <br/>
  639.          * - dC<sup>j</sup> / da <br/>
  640.          * - dC<sup>j</sup> / dk <br/>
  641.          * - dC<sup>j</sup> / dh <br/>
  642.          * - dC<sup>j</sup> / dα <br/>
  643.          * - dC<sup>j</sup> / dβ <br/>
  644.          * - dC<sup>j</sup> / dγ <br/>
  645.          * </p>
  646.          */
  647.         private final double[][] cj;

  649.         /** The S<sup>j</sup> coefficients and their derivatives.
  650.          * <p>
  651.          * Each column of the matrix contains the following values: <br/>
  652.          * - S<sup>j</sup> <br/>
  653.          * - dS<sup>j</sup> / da <br/>
  654.          * - dS<sup>j</sup> / dk <br/>
  655.          * - dS<sup>j</sup> / dh <br/>
  656.          * - dS<sup>j</sup> / dα <br/>
  657.          * - dS<sup>j</sup> / dβ <br/>
  658.          * - dS<sup>j</sup> / dγ <br/>
  659.          * </p>
  660.          */
  661.         private final double[][] sj;

  663.         /** The Coefficients C<sup>j</sup><sub>,λ</sub>.
  664.          * <p>
  665.          * See Danielson 4.2-21
  666.          * </p>
  667.          */
  668.         private final double[] cjlambda;

  670.         /** The Coefficients S<sup>j</sup><sub>,λ</sub>.
  671.         * <p>
  672.         * See Danielson 4.2-21
  673.         * </p>
  674.         */
  675.         private final double[] sjlambda;

  677.         /** Maximum value for n. */
  678.         private final int nMax;

  680.         /** Maximum value for s. */
  681.         private final int sMax;

  683.         /** Maximum value for j. */
  684.         private final int jMax;

  686.         /**
  687.          * Private constructor.
  688.          *
  689.          * @param nMax maximum value for n index
  690.          * @param sMax maximum value for s index
  691.          * @param jMax maximum value for j index
  692.          * @param context container for dynamic force model attributes
  693.          * @param aoR3Pow a / R3 up to power maxAR3Pow
  694.          * @param qns Qns coefficients
  695.          */
  696.         FourierCjSjCoefficients(final int nMax, final int sMax, final int jMax,
  697.                                 final DSSTThirdBodyDynamicContext context,
  698.                                 final double[] aoR3Pow, final double[][] qns) {
  699.             //Save parameters
  700.             this.nMax = nMax;
  701.             this.sMax = sMax;
  702.             this.jMax = jMax;

  704.             //Create objects
  705.             wnsjEtomjmsCoefficient = new WnsjEtomjmsCoefficient(context);
  706.             ABDECoefficients = new CjSjAlphaBetaKH(context);
  707.             gns = new GnsCoefficients(nMax, sMax, context, aoR3Pow, qns);

  709.             //create arays
  710.             this.cj = new double[7][jMax + 1];
  711.             this.sj = new double[7][jMax + 1];
  712.             this.cjlambda = new double[jMax];
  713.             this.sjlambda = new double[jMax];

  715.             computeCoefficients(context);
  716.         }

  718.         /**
  719.          * Compute all coefficients.
  720.          * @param context container for attributes
  721.          */
  722.         private void computeCoefficients(final DSSTThirdBodyDynamicContext context) {

  724.             final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();

  726.             for (int j = 1; j <= jMax; j++) {
  727.                 // initialise the coefficients
  728.                 for (int i = 0; i <= 6; i++) {
  729.                     cj[i][j] = 0.;
  730.                     sj[i][j] = 0.;
  731.                 }
  732.                 if (j < jMax) {
  733.                     // initialise the C<sup>j</sup><sub>,λ</sub> and S<sup>j</sup><sub>,λ</sub> coefficients
  734.                     cjlambda[j] = 0.;
  735.                     sjlambda[j] = 0.;
  736.                 }
  737.                 for (int s = 0; s <= sMax; s++) {

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

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

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

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

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

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

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

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

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

  836.                                 //Compute C<sup>j</sup><sub>,λ</sub>
  837.                                 cjlambda[j] += gns.getGns(n, s) * (wjnsemjms[0] * ABDECoefficients.getCoefD() + wmjnsemjms[0] * ABDECoefficients.getCoefE());
  838.                                 //Compute S<sup>j</sup><sub>,λ</sub>
  839.                                 sjlambda[j] += gns.getGns(n, s) * (wjnsemjms[0] * ABDECoefficients.getCoefA() + wmjnsemjms[0] * ABDECoefficients.getCoefB());
  840.                             }
  841.                         }
  842.                     }
  843.                 }
  844.                 // Divide by j
  845.                 for (int i = 0; i <= 6; i++) {
  846.                     cj[i][j] /= j;
  847.                     sj[i][j] /= j;
  848.                 }
  849.             }
  850.             //The C⁰ coefficients are not computed here.
  851.             //They are evaluated at the final point.

  853.             //C⁰<sub>,λ</sub>
  854.             cjlambda[0] = auxiliaryElements.getK() * cjlambda[1] / 2. + auxiliaryElements.getH() * sjlambda[1] / 2.;
  855.         }

  857.         /** Get the Fourier coefficient C<sup>j</sup>.
  858.          * @param j j index
  859.          * @return C<sup>j</sup>
  860.          */
  861.         public double getCj(final int j) {
  862.             return cj[0][j];
  863.         }

  865.         /** Get the derivative dC<sup>j</sup>/da.
  866.          * @param j j index
  867.          * @return dC<sup>j</sup>/da
  868.          */
  869.         public double getdCjda(final int j) {
  870.             return cj[1][j];
  871.         }

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

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

  889.         /** Get the derivative dC<sup>j</sup>/dα.
  890.          * @param j j index
  891.          * @return dC<sup>j</sup>/dα
  892.          */
  893.         public double getdCjdalpha(final int j) {
  894.             return cj[4][j];
  895.         }

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

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

  913.         /** Get the Fourier coefficient S<sup>j</sup>.
  914.          * @param j j index
  915.          * @return S<sup>j</sup>
  916.          */
  917.         public double getSj(final int j) {
  918.             return sj[0][j];
  919.         }

  921.         /** Get the derivative dS<sup>j</sup>/da.
  922.          * @param j j index
  923.          * @return dS<sup>j</sup>/da
  924.          */
  925.         public double getdSjda(final int j) {
  926.             return sj[1][j];
  927.         }

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

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

  945.         /** Get the derivative dS<sup>j</sup>/dα.
  946.          * @param j j index
  947.          * @return dS<sup>j</sup>/dα
  948.          */
  949.         public double getdSjdalpha(final int j) {
  950.             return sj[4][j];
  951.         }

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

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

  969.         /** Get the coefficient C⁰<sub>,λ</sub>.
  970.          * @return C⁰<sub>,λ</sub>
  971.          */
  972.         public double getC0Lambda() {
  973.             return cjlambda[0];
  974.         }

  976.         /** Get the coefficient C<sup>j</sup><sub>,λ</sub>.
  977.          * @param j j index
  978.          * @return C<sup>j</sup><sub>,λ</sub>
  979.          */
  980.         public double getCjLambda(final int j) {
  981.             if (j < 1 || j >= jMax) {
  982.                 return 0.;
  983.             }
  984.             return cjlambda[j];
  985.         }

  987.         /** Get the coefficient S<sup>j</sup><sub>,λ</sub>.
  988.          * @param j j index
  989.          * @return S<sup>j</sup><sub>,λ</sub>
  990.          */
  991.         public double getSjLambda(final int j) {
  992.             if (j < 1 || j >= jMax) {
  993.                 return 0.;
  994.             }
  995.             return sjlambda[j];
  996.         }
  997.     }

  999.     /** Computes the C<sup>j</sup> and S<sup>j</sup> coefficients Danielson 4.2-(15,16)
  1000.      * and their derivatives.
  1001.      *  <p>
  1002.      *  CS Mathematical Report $3.5.3.2
  1003.      *  </p>
  1004.      */
  1005.     private class FieldFourierCjSjCoefficients <T extends CalculusFieldElement<T>> {

  1007.         /** The coefficients G<sub>n, s</sub> and their derivatives. */
  1008.         private final FieldGnsCoefficients<T> gns;

  1010.         /** the coefficients e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and their derivatives by h and k. */
  1011.         private final FieldWnsjEtomjmsCoefficient<T> wnsjEtomjmsCoefficient;

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

  1016.         /** The Fourier coefficients C<sup>j</sup> and their derivatives.
  1017.          * <p>
  1018.          * Each column of the matrix contains the following values: <br/>
  1019.          * - C<sup>j</sup> <br/>
  1020.          * - dC<sup>j</sup> / da <br/>
  1021.          * - dC<sup>j</sup> / dk <br/>
  1022.          * - dC<sup>j</sup> / dh <br/>
  1023.          * - dC<sup>j</sup> / dα <br/>
  1024.          * - dC<sup>j</sup> / dβ <br/>
  1025.          * - dC<sup>j</sup> / dγ <br/>
  1026.          * </p>
  1027.          */
  1028.         private final T[][] cj;

  1030.         /** The S<sup>j</sup> coefficients and their derivatives.
  1031.          * <p>
  1032.          * Each column of the matrix contains the following values: <br/>
  1033.          * - S<sup>j</sup> <br/>
  1034.          * - dS<sup>j</sup> / da <br/>
  1035.          * - dS<sup>j</sup> / dk <br/>
  1036.          * - dS<sup>j</sup> / dh <br/>
  1037.          * - dS<sup>j</sup> / dα <br/>
  1038.          * - dS<sup>j</sup> / dβ <br/>
  1039.          * - dS<sup>j</sup> / dγ <br/>
  1040.          * </p>
  1041.          */
  1042.         private final T[][] sj;

  1044.         /** The Coefficients C<sup>j</sup><sub>,λ</sub>.
  1045.          * <p>
  1046.          * See Danielson 4.2-21
  1047.          * </p>
  1048.          */
  1049.         private final T[] cjlambda;

  1051.         /** The Coefficients S<sup>j</sup><sub>,λ</sub>.
  1052.         * <p>
  1053.         * See Danielson 4.2-21
  1054.         * </p>
  1055.         */
  1056.         private final T[] sjlambda;

  1058.         /** Zero. */
  1059.         private final T zero;

  1061.         /** Maximum value for n. */
  1062.         private final int nMax;

  1064.         /** Maximum value for s. */
  1065.         private final int sMax;

  1067.         /** Maximum value for j. */
  1068.         private final int jMax;

  1070.         /**
  1071.          * Private constructor.
  1072.          *
  1073.          * @param nMax maximum value for n index
  1074.          * @param sMax maximum value for s index
  1075.          * @param jMax maximum value for j index
  1076.          * @param context container for dynamic force model attributes
  1077.          * @param aoR3Pow a / R3 up to power maxAR3Pow
  1078.          * @param qns Qns coefficients
  1079.          * @param field field used by default
  1080.          */
  1081.         FieldFourierCjSjCoefficients(final int nMax, final int sMax, final int jMax,
  1082.                                      final FieldDSSTThirdBodyDynamicContext<T> context,
  1083.                                      final T[] aoR3Pow, final T[][] qns,
  1084.                                      final Field<T> field) {
  1085.             //Zero
  1086.             this.zero = field.getZero();

  1088.             //Save parameters
  1089.             this.nMax = nMax;
  1090.             this.sMax = sMax;
  1091.             this.jMax = jMax;

  1093.             //Create objects
  1094.             wnsjEtomjmsCoefficient = new FieldWnsjEtomjmsCoefficient<>(context, field);
  1095.             ABDECoefficients       = new FieldCjSjAlphaBetaKH<>(context, field);
  1096.             gns                    = new FieldGnsCoefficients<>(nMax, sMax, context, aoR3Pow, qns, field);

  1098.             //create arays
  1099.             this.cj = MathArrays.buildArray(field, 7, jMax + 1);
  1100.             this.sj = MathArrays.buildArray(field, 7, jMax + 1);
  1101.             this.cjlambda = MathArrays.buildArray(field, jMax);
  1102.             this.sjlambda = MathArrays.buildArray(field, jMax);

  1104.             computeCoefficients(context, field);
  1105.         }

  1107.         /**
  1108.          * Compute all coefficients.
  1109.          * @param context container for attributes
  1110.          * @param field field used by default
  1111.          */
  1112.         private void computeCoefficients(final FieldDSSTThirdBodyDynamicContext<T> context,
  1113.                                          final Field<T> field) {

  1115.             final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();

  1117.             for (int j = 1; j <= jMax; j++) {
  1118.                 // initialise the coefficients
  1119.                 for (int i = 0; i <= 6; i++) {
  1120.                     cj[i][j] = zero;
  1121.                     sj[i][j] = zero;
  1122.                 }
  1123.                 if (j < jMax) {
  1124.                     // initialise the C<sup>j</sup><sub>,λ</sub> and S<sup>j</sup><sub>,λ</sub> coefficients
  1125.                     cjlambda[j] = zero;
  1126.                     sjlambda[j] = zero;
  1127.                 }
  1128.                 for (int s = 0; s <= sMax; s++) {

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

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

  1136.                     for (int n = minN; n <= nMax; n++) {
  1137.                         // check if n-s is even
  1138.                         if ((n - s) % 2 == 0) {
  1139.                             // compute the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n+1, s</sup> and its derivatives
  1140.                             final T[] wjnp1semjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(j, s, n + 1, context, field);

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

  1145.                             // compute common factors
  1146.                             final T coef1 = (wjnp1semjms[0].multiply(ABDECoefficients.getCoefA()).add(wmjnp1semjms[0].multiply(ABDECoefficients.getCoefB()))).negate();
  1147.                             final T coef2 =  wjnp1semjms[0].multiply(ABDECoefficients.getCoefD()).add(wmjnp1semjms[0].multiply(ABDECoefficients.getCoefE()));

  1149.                             //Compute C<sup>j</sup>
  1150.                             cj[0][j] = cj[0][j].add(gns.getGns(n, s).multiply(coef1));
  1151.                             //Compute dC<sup>j</sup> / da
  1152.                             cj[1][j] = cj[1][j].add(gns.getdGnsda(n, s).multiply(coef1));
  1153.                             //Compute dC<sup>j</sup> / dk
  1154.                             cj[2][j] = cj[2][j].add(gns.getGns(n, s).negate().
  1155.                                        multiply(
  1156.                                             wjnp1semjms[1].multiply(ABDECoefficients.getCoefA()).
  1157.                                             add(wjnp1semjms[0].multiply(ABDECoefficients.getdCoefAdk())).
  1158.                                             add(wmjnp1semjms[1].multiply(ABDECoefficients.getCoefB())).
  1159.                                             add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefBdk()))
  1160.                                          ));
  1161.                             //Compute dC<sup>j</sup> / dh
  1162.                             cj[3][j] = cj[3][j].add(gns.getGns(n, s).negate().
  1163.                                        multiply(
  1164.                                              wjnp1semjms[2].multiply(ABDECoefficients.getCoefA()).
  1165.                                              add(wjnp1semjms[0].multiply(ABDECoefficients.getdCoefAdh())).
  1166.                                              add(wmjnp1semjms[2].multiply(ABDECoefficients.getCoefB())).
  1167.                                              add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefBdh()))
  1168.                                              ));
  1169.                             //Compute dC<sup>j</sup> / dα
  1170.                             cj[4][j] = cj[4][j].add(gns.getGns(n, s).negate().
  1171.                                        multiply(
  1172.                                            wjnp1semjms[0].multiply(ABDECoefficients.getdCoefAdalpha()).
  1173.                                            add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefBdalpha()))
  1174.                                        ));
  1175.                             //Compute dC<sup>j</sup> / dβ
  1176.                             cj[5][j] = cj[5][j].add(gns.getGns(n, s).negate().
  1177.                                        multiply(
  1178.                                            wjnp1semjms[0].multiply(ABDECoefficients.getdCoefAdbeta()).
  1179.                                            add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefBdbeta()))
  1180.                                        ));
  1181.                             //Compute dC<sup>j</sup> / dγ
  1182.                             cj[6][j] = cj[6][j].add(gns.getdGnsdgamma(n, s).multiply(coef1));

  1184.                             //Compute S<sup>j</sup>
  1185.                             sj[0][j] = sj[0][j].add(gns.getGns(n, s).multiply(coef2));
  1186.                             //Compute dS<sup>j</sup> / da
  1187.                             sj[1][j] = sj[1][j].add(gns.getdGnsda(n, s).multiply(coef2));
  1188.                             //Compute dS<sup>j</sup> / dk
  1189.                             sj[2][j] = sj[2][j].add(gns.getGns(n, s).
  1190.                                        multiply(
  1191.                                            wjnp1semjms[1].multiply(ABDECoefficients.getCoefD()).
  1192.                                            add(wjnp1semjms[0].multiply(ABDECoefficients.getdCoefDdk())).
  1193.                                            add(wmjnp1semjms[1].multiply(ABDECoefficients.getCoefE())).
  1194.                                            add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefEdk()))
  1195.                                        ));
  1196.                             //Compute dS<sup>j</sup> / dh
  1197.                             sj[3][j] = sj[3][j].add(gns.getGns(n, s).
  1198.                                        multiply(
  1199.                                            wjnp1semjms[2].multiply(ABDECoefficients.getCoefD()).
  1200.                                            add(wjnp1semjms[0].multiply(ABDECoefficients.getdCoefDdh())).
  1201.                                            add(wmjnp1semjms[2].multiply(ABDECoefficients.getCoefE())).
  1202.                                            add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefEdh()))
  1203.                                        ));
  1204.                             //Compute dS<sup>j</sup> / dα
  1205.                             sj[4][j] = sj[4][j].add(gns.getGns(n, s).
  1206.                                        multiply(
  1207.                                             wjnp1semjms[0].multiply(ABDECoefficients.getdCoefDdalpha()).
  1208.                                             add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefEdalpha()))
  1209.                                        ));
  1210.                             //Compute dS<sup>j</sup> / dβ
  1211.                             sj[5][j] = sj[5][j].add(gns.getGns(n, s).
  1212.                                         multiply(
  1213.                                             wjnp1semjms[0].multiply(ABDECoefficients.getdCoefDdbeta()).
  1214.                                             add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefEdbeta()))
  1215.                                        ));
  1216.                             //Compute dS<sup>j</sup> / dγ
  1217.                             sj[6][j] = sj[6][j].add(gns.getdGnsdgamma(n, s).multiply(coef2));

  1219.                             //Check if n is greater or equal to j and j is at most jMax-1
  1220.                             if (n >= j && j < jMax) {
  1221.                                 // compute the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives
  1222.                                 final T[] wjnsemjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(j, s, n, context, field);

  1224.                                 // compute the coefficient e<sup>-|j-s|</sup>*w<sub>-j</sub><sup>n, s</sup> and its derivatives
  1225.                                 final T[] wmjnsemjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(-j, s, n, context, field);

  1227.                                 //Compute C<sup>j</sup><sub>,λ</sub>
  1228.                                 cjlambda[j] = cjlambda[j].add(gns.getGns(n, s).multiply(wjnsemjms[0].multiply(ABDECoefficients.getCoefD()).add(wmjnsemjms[0].multiply(ABDECoefficients.getCoefE()))));
  1229.                                 //Compute S<sup>j</sup><sub>,λ</sub>
  1230.                                 sjlambda[j] = sjlambda[j].add(gns.getGns(n, s).multiply(wjnsemjms[0].multiply(ABDECoefficients.getCoefA()).add(wmjnsemjms[0].multiply(ABDECoefficients.getCoefB()))));
  1231.                             }
  1232.                         }
  1233.                     }
  1234.                 }
  1235.                 // Divide by j
  1236.                 for (int i = 0; i <= 6; i++) {
  1237.                     cj[i][j] = cj[i][j].divide(j);
  1238.                     sj[i][j] = sj[i][j].divide(j);
  1239.                 }
  1240.             }
  1241.             //The C⁰ coefficients are not computed here.
  1242.             //They are evaluated at the final point.

  1244.             //C⁰<sub>,λ</sub>
  1245.             cjlambda[0] = auxiliaryElements.getK().multiply(cjlambda[1]).divide(2.).add(auxiliaryElements.getH().multiply(sjlambda[1]).divide(2.));
  1246.         }

  1248.         /** Get the Fourier coefficient C<sup>j</sup>.
  1249.          * @param j j index
  1250.          * @return C<sup>j</sup>
  1251.          */
  1252.         public T getCj(final int j) {
  1253.             return cj[0][j];
  1254.         }

  1256.         /** Get the derivative dC<sup>j</sup>/da.
  1257.          * @param j j index
  1258.          * @return dC<sup>j</sup>/da
  1259.          */
  1260.         public T getdCjda(final int j) {
  1261.             return cj[1][j];
  1262.         }

  1264.         /** Get the derivative dC<sup>j</sup>/dk.
  1265.          * @param j j index
  1266.          * @return dC<sup>j</sup>/dk
  1267.          */
  1268.         public T getdCjdk(final int j) {
  1269.             return cj[2][j];
  1270.         }

  1272.         /** Get the derivative dC<sup>j</sup>/dh.
  1273.          * @param j j index
  1274.          * @return dC<sup>j</sup>/dh
  1275.          */
  1276.         public T getdCjdh(final int j) {
  1277.             return cj[3][j];
  1278.         }

  1280.         /** Get the derivative dC<sup>j</sup>/dα.
  1281.          * @param j j index
  1282.          * @return dC<sup>j</sup>/dα
  1283.          */
  1284.         public T getdCjdalpha(final int j) {
  1285.             return cj[4][j];
  1286.         }

  1288.         /** Get the derivative dC<sup>j</sup>/dβ.
  1289.          * @param j j index
  1290.          * @return dC<sup>j</sup>/dβ
  1291.          */
  1292.         public T getdCjdbeta(final int j) {
  1293.             return cj[5][j];
  1294.         }

  1296.         /** Get the derivative dC<sup>j</sup>/dγ.
  1297.          * @param j j index
  1298.          * @return dC<sup>j</sup>/dγ
  1299.          */
  1300.         public T getdCjdgamma(final int j) {
  1301.             return cj[6][j];
  1302.         }

  1304.         /** Get the Fourier coefficient S<sup>j</sup>.
  1305.          * @param j j index
  1306.          * @return S<sup>j</sup>
  1307.          */
  1308.         public T getSj(final int j) {
  1309.             return sj[0][j];
  1310.         }

  1312.         /** Get the derivative dS<sup>j</sup>/da.
  1313.          * @param j j index
  1314.          * @return dS<sup>j</sup>/da
  1315.          */
  1316.         public T getdSjda(final int j) {
  1317.             return sj[1][j];
  1318.         }

  1320.         /** Get the derivative dS<sup>j</sup>/dk.
  1321.          * @param j j index
  1322.          * @return dS<sup>j</sup>/dk
  1323.          */
  1324.         public T getdSjdk(final int j) {
  1325.             return sj[2][j];
  1326.         }

  1328.         /** Get the derivative dS<sup>j</sup>/dh.
  1329.          * @param j j index
  1330.          * @return dS<sup>j</sup>/dh
  1331.          */
  1332.         public T getdSjdh(final int j) {
  1333.             return sj[3][j];
  1334.         }

  1336.         /** Get the derivative dS<sup>j</sup>/dα.
  1337.          * @param j j index
  1338.          * @return dS<sup>j</sup>/dα
  1339.          */
  1340.         public T getdSjdalpha(final int j) {
  1341.             return sj[4][j];
  1342.         }

  1344.         /** Get the derivative dS<sup>j</sup>/dβ.
  1345.          * @param j j index
  1346.          * @return dS<sup>j</sup>/dβ
  1347.          */
  1348.         public T getdSjdbeta(final int j) {
  1349.             return sj[5][j];
  1350.         }

  1352.         /** Get the derivative dS<sup>j</sup>/dγ.
  1353.          * @param j j index
  1354.          * @return dS<sup>j</sup>/dγ
  1355.          */
  1356.         public T getdSjdgamma(final int j) {
  1357.             return sj[6][j];
  1358.         }

  1360.         /** Get the coefficient C⁰<sub>,λ</sub>.
  1361.          * @return C⁰<sub>,λ</sub>
  1362.          */
  1363.         public T getC0Lambda() {
  1364.             return cjlambda[0];
  1365.         }

  1367.         /** Get the coefficient C<sup>j</sup><sub>,λ</sub>.
  1368.          * @param j j index
  1369.          * @return C<sup>j</sup><sub>,λ</sub>
  1370.          */
  1371.         public T getCjLambda(final int j) {
  1372.             if (j < 1 || j >= jMax) {
  1373.                 return zero;
  1374.             }
  1375.             return cjlambda[j];
  1376.         }

  1378.         /** Get the coefficient S<sup>j</sup><sub>,λ</sub>.
  1379.          * @param j j index
  1380.          * @return S<sup>j</sup><sub>,λ</sub>
  1381.          */
  1382.         public T getSjLambda(final int j) {
  1383.             if (j < 1 || j >= jMax) {
  1384.                 return zero;
  1385.             }
  1386.             return sjlambda[j];
  1387.         }
  1388.     }

  1390.     /** 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.
  1391.      *
  1392.      * <p>
  1393.      * Starting from Danielson 4.2-9,10,11 and taking into account that fact that: <br />
  1394.      * c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
  1395.      * the expression e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup>
  1396.      * can be written as: <br/ >
  1397.      * - for |s| > |j| <br />
  1398.      * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
  1399.      *          (((n + s)!(n - s)!)/((n + j)!(n - j)!)) *
  1400.      *          (-b)<sup>|j-s|</sup> *
  1401.      *          ((1 - c²)<sup>n-|s|</sup>/(1 + c²)<sup>n</sup>) *
  1402.      *          P<sub>n-|s|</sub><sup>|j-s|, |j+s|</sup>(χ) <br />
  1403.      * <br />
  1404.      * - for |s| <= |j| <br />
  1405.      * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
  1406.      *          (-b)<sup>|j-s|</sup> *
  1407.      *          ((1 - c²)<sup>n-|j|</sup>/(1 + c²)<sup>n</sup>) *
  1408.      *          P<sub>n-|j|</sub><sup>|j-s|, |j+s|</sup>(χ)
  1409.      * </p>
  1410.      *
  1411.      * @author Lucian Barbulescu
  1412.      */
  1413.     private class WnsjEtomjmsCoefficient {

  1415.         /** The value c.
  1416.          * <p>
  1417.          *  c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
  1418.          * </p>
  1419.          *  */
  1420.         private final double c;

  1422.         /** db / dh. */
  1423.         private final double dbdh;

  1425.         /** db / dk. */
  1426.         private final double dbdk;

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

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

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

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

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

  1445.         /**
  1446.          * Standard constructor.
  1447.          * @param context container for attributes
  1448.          */
  1449.         WnsjEtomjmsCoefficient(final DSSTThirdBodyDynamicContext context) {

  1451.             final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();

  1453.             //initialise fields
  1454.             c = auxiliaryElements.getEcc() * context.getb();
  1455.             final double c2 = c * c;

  1457.             //b² * χ
  1458.             final double b2Chi = context.getb() * context.getb() * context.getX();
  1459.             //Compute derivatives of b
  1460.             dbdh = auxiliaryElements.getH() * b2Chi;
  1461.             dbdk = auxiliaryElements.getK() * b2Chi;

  1463.             //Compute derivatives of c
  1464.             if (auxiliaryElements.getEcc() == 0.0) {
  1465.                 // we are at a perfectly circular orbit singularity here
  1466.                 // we arbitrarily consider the perigee is along the X axis,
  1467.                 // i.e cos(ω + Ω) = h/ecc 1 and sin(ω + Ω) = k/ecc = 0
  1468.                 dcdh = auxiliaryElements.getEcc() * dbdh + context.getb();
  1469.                 dcdk = auxiliaryElements.getEcc() * dbdk;
  1470.             } else {
  1471.                 dcdh = auxiliaryElements.getEcc() * dbdh + context.getb() * auxiliaryElements.getH() / auxiliaryElements.getEcc();
  1472.                 dcdk = auxiliaryElements.getEcc() * dbdk + context.getb() * auxiliaryElements.getK() / auxiliaryElements.getEcc();
  1473.             }

  1475.             //Compute the powers (1 - c²)<sup>n</sup> and (1 + c²)<sup>n</sup>
  1476.             omc2tn = new double[staticContext.getMaxAR3Pow() + staticContext.getMaxFreqF() + 2];
  1477.             opc2tn = new double[staticContext.getMaxAR3Pow() + staticContext.getMaxFreqF() + 2];
  1478.             final double omc2 = 1. - c2;
  1479.             final double opc2 = 1. + c2;
  1480.             omc2tn[0] = 1.;
  1481.             opc2tn[0] = 1.;
  1482.             for (int i = 1; i <= staticContext.getMaxAR3Pow() + staticContext.getMaxFreqF() + 1; i++) {
  1483.                 omc2tn[i] = omc2tn[i - 1] * omc2;
  1484.                 opc2tn[i] = opc2tn[i - 1] * opc2;
  1485.             }

  1487.             //Compute the powers of b
  1488.             btjms = new double[staticContext.getMaxAR3Pow() + staticContext.getMaxFreqF() + 1];
  1489.             btjms[0] = 1.;
  1490.             for (int i = 1; i <= staticContext.getMaxAR3Pow() + staticContext.getMaxFreqF(); i++) {
  1491.                 btjms[i] = btjms[i - 1] * context.getb();
  1492.             }
  1493.         }

  1495.         /** 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 />
  1496.          *
  1497.          * @param j j index
  1498.          * @param s s index
  1499.          * @param n n index
  1500.          * @param context container for attributes
  1501.          * @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
  1502.          */
  1503.         public double[] computeWjnsEmjmsAndDeriv(final int j, final int s, final int n, final DSSTThirdBodyDynamicContext context) {
  1504.             final double[] wjnsemjms = new double[] {0., 0., 0.};

  1506.             // |j|
  1507.             final int absJ = FastMath.abs(j);
  1508.             // |s|
  1509.             final int absS = FastMath.abs(s);
  1510.             // |j - s|
  1511.             final int absJmS = FastMath.abs(j - s);
  1512.             // |j + s|
  1513.             final int absJpS = FastMath.abs(j + s);

  1515.             //The lower index of P. Also the power of (1 - c²)
  1516.             final int l;
  1517.             // the factorial ratio coefficient or 1. if |s| <= |j|
  1518.             final double factCoef;
  1519.             if (absS > absJ) {
  1520.                 //factCoef = (fact[n + s] / fact[n + j]) * (fact[n - s] / fact[n - j]);
  1521.                 factCoef = (CombinatoricsUtils.factorialDouble(n + s) / CombinatoricsUtils.factorialDouble(n + j)) * (CombinatoricsUtils.factorialDouble(n - s) / CombinatoricsUtils.factorialDouble(n - j));
  1522.                 l = n - absS;
  1523.             } else {
  1524.                 factCoef = 1.;
  1525.                 l = n - absJ;
  1526.             }

  1528.             // (-1)<sup>|j-s|</sup>
  1529.             final double sign = absJmS % 2 != 0 ? -1. : 1.;
  1530.             //(1 - c²)<sup>n-|s|</sup> / (1 + c²)<sup>n</sup>
  1531.             final double coef1 = omc2tn[l] / opc2tn[n];
  1532.             //-b<sup>|j-s|</sup>
  1533.             final double coef2 = sign * btjms[absJmS];
  1534.             // P<sub>l</sub><sup>|j-s|, |j+s|</sup>(χ)
  1535.             // Jacobi polynomial value (0) and first-order derivative (1)
  1536.             final double[] jac = JacobiPolynomials.getValueAndDerivative(l, absJmS, absJpS, context.getX());

  1538.             // the derivative of coef1 by c
  1539.             final double dcoef1dc = -coef1 * 2. * c * (((double) n) / opc2tn[1] + ((double) l) / omc2tn[1]);
  1540.             // the derivative of coef1 by h
  1541.             final double dcoef1dh = dcoef1dc * dcdh;
  1542.             // the derivative of coef1 by k
  1543.             final double dcoef1dk = dcoef1dc * dcdk;

  1545.             // the derivative of coef2 by b
  1546.             final double dcoef2db = absJmS == 0 ? 0 : sign * (double) absJmS * btjms[absJmS - 1];
  1547.             // the derivative of coef2 by h
  1548.             final double dcoef2dh = dcoef2db * dbdh;
  1549.             // the derivative of coef2 by k
  1550.             final double dcoef2dk = dcoef2db * dbdk;

  1552.             // the jacobi polynomial value
  1553.             // final double jacobi = jac.getValue();
  1554.             final double jacobi = jac[0];
  1555.             // the derivative of the Jacobi polynomial by h
  1556.             //final double djacobidh = jac.getGradient()[0] * context.getHXXX();
  1557.             final double djacobidh = jac[1] * context.getHXXX();
  1558.             // the derivative of the Jacobi polynomial by k
  1559.             //final double djacobidk = jac.getGradient()[0] * context.getKXXX();
  1560.             final double djacobidk = jac[1] * context.getKXXX();

  1562.             //group the above coefficients to limit the mathematical operations
  1563.             final double term1 = factCoef * coef1 * coef2;
  1564.             final double term2 = factCoef * coef1 * jacobi;
  1565.             final double term3 = factCoef * coef2 * jacobi;

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

  1572.             return wjnsemjms;
  1573.         }
  1574.     }

  1576.     /** 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.
  1577.     *
  1578.     * <p>
  1579.     * Starting from Danielson 4.2-9,10,11 and taking into account that fact that: <br />
  1580.     * c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
  1581.     * the expression e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup>
  1582.     * can be written as: <br/ >
  1583.     * - for |s| > |j| <br />
  1584.     * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
  1585.     *          (((n + s)!(n - s)!)/((n + j)!(n - j)!)) *
  1586.     *          (-b)<sup>|j-s|</sup> *
  1587.     *          ((1 - c²)<sup>n-|s|</sup>/(1 + c²)<sup>n</sup>) *
  1588.     *          P<sub>n-|s|</sub><sup>|j-s|, |j+s|</sup>(χ) <br />
  1589.     * <br />
  1590.     * - for |s| <= |j| <br />
  1591.     * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
  1592.     *          (-b)<sup>|j-s|</sup> *
  1593.     *          ((1 - c²)<sup>n-|j|</sup>/(1 + c²)<sup>n</sup>) *
  1594.     *          P<sub>n-|j|</sub><sup>|j-s|, |j+s|</sup>(χ)
  1595.     * </p>
  1596.     *
  1597.     * @author Lucian Barbulescu
  1598.     */
  1599.     private class FieldWnsjEtomjmsCoefficient <T extends CalculusFieldElement<T>> {

  1601.         /** The value c.
  1602.          * <p>
  1603.          *  c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
  1604.          * </p>
  1605.          *  */
  1606.         private final T c;

  1608.         /** db / dh. */
  1609.         private final T dbdh;

  1611.         /** db / dk. */
  1612.         private final T dbdk;

  1614.         /** dc / dh = e * db/dh + b * de/dh. */
  1615.         private final T dcdh;

  1617.         /** dc / dk = e * db/dk + b * de/dk. */
  1618.         private final T dcdk;

  1620.         /** The values (1 - c²)<sup>n</sup>. <br />
  1621.          * The maximum possible value for the power is N + 1 */
  1622.         private final T[] omc2tn;

  1624.         /** The values (1 + c²)<sup>n</sup>. <br />
  1625.          * The maximum possible value for the power is N + 1 */
  1626.         private final T[] opc2tn;

  1628.         /** The values b<sup>|j-s|</sup>. */
  1629.         private final T[] btjms;

  1631.         /**
  1632.          * Standard constructor.
  1633.          * @param context container for attributes
  1634.          * @param field field used by default
  1635.          */
  1636.         FieldWnsjEtomjmsCoefficient(final FieldDSSTThirdBodyDynamicContext<T> context, final Field<T> field) {

  1638.             final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();

  1640.             //Zero
  1641.             final T zero = field.getZero();

  1643.             //initialise fields
  1644.             c = auxiliaryElements.getEcc().multiply(context.getb());
  1645.             final T c2 = c.multiply(c);

  1647.             //b² * χ
  1648.             final T b2Chi = context.getb().multiply(context.getb()).multiply(context.getX());
  1649.             //Compute derivatives of b
  1650.             dbdh = auxiliaryElements.getH().multiply(b2Chi);
  1651.             dbdk = auxiliaryElements.getK().multiply(b2Chi);

  1653.             //Compute derivatives of c
  1654.             if (auxiliaryElements.getEcc().getReal() == 0.0) {
  1655.                 // we are at a perfectly circular orbit singularity here
  1656.                 // we arbitrarily consider the perigee is along the X axis,
  1657.                 // i.e cos(ω + Ω) = h/ecc 1 and sin(ω + Ω) = k/ecc = 0
  1658.                 dcdh = auxiliaryElements.getEcc().multiply(dbdh).add(context.getb());
  1659.                 dcdk = auxiliaryElements.getEcc().multiply(dbdk);
  1660.             } else {
  1661.                 dcdh = auxiliaryElements.getEcc().multiply(dbdh).add(context.getb().multiply(auxiliaryElements.getH()).divide(auxiliaryElements.getEcc()));
  1662.                 dcdk = auxiliaryElements.getEcc().multiply(dbdk).add(context.getb().multiply(auxiliaryElements.getK()).divide(auxiliaryElements.getEcc()));
  1663.             }

  1665.             //Compute the powers (1 - c²)<sup>n</sup> and (1 + c²)<sup>n</sup>
  1666.             omc2tn = MathArrays.buildArray(field, staticContext.getMaxAR3Pow() + staticContext.getMaxFreqF() + 2);
  1667.             opc2tn = MathArrays.buildArray(field, staticContext.getMaxAR3Pow() + staticContext.getMaxFreqF() + 2);
  1668.             final T omc2 = c2.negate().add(1.);
  1669.             final T opc2 = c2.add(1.);
  1670.             omc2tn[0] = zero.add(1.);
  1671.             opc2tn[0] = zero.add(1.);
  1672.             for (int i = 1; i <= staticContext.getMaxAR3Pow() + staticContext.getMaxFreqF() + 1; i++) {
  1673.                 omc2tn[i] = omc2tn[i - 1].multiply(omc2);
  1674.                 opc2tn[i] = opc2tn[i - 1].multiply(opc2);
  1675.             }

  1677.             //Compute the powers of b
  1678.             btjms = MathArrays.buildArray(field, staticContext.getMaxAR3Pow() + staticContext.getMaxFreqF() + 1);
  1679.             btjms[0] = zero.add(1.);
  1680.             for (int i = 1; i <= staticContext.getMaxAR3Pow() + staticContext.getMaxFreqF(); i++) {
  1681.                 btjms[i] = btjms[i - 1].multiply(context.getb());
  1682.             }
  1683.         }

  1685.         /** 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 />
  1686.          *
  1687.          * @param j j index
  1688.          * @param s s index
  1689.          * @param n n index
  1690.          * @param context container for attributes
  1691.          * @param field field used by default
  1692.          * @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
  1693.          */
  1694.         public T[] computeWjnsEmjmsAndDeriv(final int j, final int s, final int n,
  1695.                                             final FieldDSSTThirdBodyDynamicContext<T> context,
  1696.                                             final Field<T> field) {
  1697.             //Zero
  1698.             final T zero = field.getZero();

  1700.             final T[] wjnsemjms = MathArrays.buildArray(field, 3);
  1701.             Arrays.fill(wjnsemjms, zero);

  1703.             // |j|
  1704.             final int absJ = FastMath.abs(j);
  1705.             // |s|
  1706.             final int absS = FastMath.abs(s);
  1707.             // |j - s|
  1708.             final int absJmS = FastMath.abs(j - s);
  1709.             // |j + s|
  1710.             final int absJpS = FastMath.abs(j + s);

  1712.             //The lower index of P. Also the power of (1 - c²)
  1713.             final int l;
  1714.             // the factorial ratio coefficient or 1. if |s| <= |j|
  1715.             final T factCoef;
  1716.             if (absS > absJ) {
  1717.                 //factCoef = (fact[n + s] / fact[n + j]) * (fact[n - s] / fact[n - j]);
  1718.                 factCoef = zero.add((CombinatoricsUtils.factorialDouble(n + s) / CombinatoricsUtils.factorialDouble(n + j)) * (CombinatoricsUtils.factorialDouble(n - s) / CombinatoricsUtils.factorialDouble(n - j)));
  1719.                 l = n - absS;
  1720.             } else {
  1721.                 factCoef = zero.add(1.);
  1722.                 l = n - absJ;
  1723.             }

  1725.             // (-1)<sup>|j-s|</sup>
  1726.             final T sign = absJmS % 2 != 0 ? zero.add(-1.) : zero.add(1.);
  1727.             //(1 - c²)<sup>n-|s|</sup> / (1 + c²)<sup>n</sup>
  1728.             final T coef1 = omc2tn[l].divide(opc2tn[n]);
  1729.             //-b<sup>|j-s|</sup>
  1730.             final T coef2 = btjms[absJmS].multiply(sign);
  1731.             // P<sub>l</sub><sup>|j-s|, |j+s|</sup>(χ)
  1732.             final FieldGradient<T> jac =
  1733.                     JacobiPolynomials.getValue(l, absJmS, absJpS, FieldGradient.variable(1, 0, context.getX()));

  1735.             // the derivative of coef1 by c
  1736.             final T dcoef1dc = coef1.negate().multiply(2.).multiply(c).multiply(opc2tn[1].reciprocal().multiply(n).add(omc2tn[1].reciprocal().multiply(l)));
  1737.             // the derivative of coef1 by h
  1738.             final T dcoef1dh = dcoef1dc.multiply(dcdh);
  1739.             // the derivative of coef1 by k
  1740.             final T dcoef1dk = dcoef1dc.multiply(dcdk);

  1742.             // the derivative of coef2 by b
  1743.             final T dcoef2db = absJmS == 0 ? zero : sign.multiply(absJmS).multiply(btjms[absJmS - 1]);
  1744.             // the derivative of coef2 by h
  1745.             final T dcoef2dh = dcoef2db.multiply(dbdh);
  1746.             // the derivative of coef2 by k
  1747.             final T dcoef2dk = dcoef2db.multiply(dbdk);

  1749.             // the jacobi polynomial value
  1750.             final T jacobi = jac.getValue();
  1751.             // the derivative of the Jacobi polynomial by h
  1752.             final T djacobidh = jac.getGradient()[0].multiply(context.getHXXX());
  1753.             // the derivative of the Jacobi polynomial by k
  1754.             final T djacobidk = jac.getGradient()[0].multiply(context.getKXXX());

  1756.             //group the above coefficients to limit the mathematical operations
  1757.             final T term1 = factCoef.multiply(coef1).multiply(coef2);
  1758.             final T term2 = factCoef.multiply(coef1).multiply(jacobi);
  1759.             final T term3 = factCoef.multiply(coef2).multiply(jacobi);

  1761.             //compute e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives by k and h
  1762.             wjnsemjms[0] = term1.multiply(jacobi);
  1763.             wjnsemjms[1] = dcoef1dk.multiply(term3).add(dcoef2dk.multiply(term2)).add(djacobidk.multiply(term1));
  1764.             wjnsemjms[2] = dcoef1dh.multiply(term3).add(dcoef2dh.multiply(term2)).add(djacobidh.multiply(term1));

  1766.             return wjnsemjms;
  1767.         }
  1768.     }

  1770.     /** The G<sub>n,s</sub> coefficients and their derivatives.
  1771.      * <p>
  1772.      * See Danielson, 4.2-17
  1773.      *
  1774.      * @author Lucian Barbulescu
  1775.      */
  1776.     private class GnsCoefficients {

  1778.         /** Maximum value for n index. */
  1779.         private final int nMax;

  1781.         /** Maximum value for s index. */
  1782.         private final int sMax;

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

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

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

  1793.         /** Standard constructor.
  1794.          *
  1795.          * @param nMax maximum value for n indes
  1796.          * @param sMax maximum value for s index
  1797.          * @param context container for attributes
  1798.          * @param aoR3Pow a / R3 up to power maxAR3Pow
  1799.          * @param qns Qns coefficients
  1800.          */
  1801.         GnsCoefficients(final int nMax, final int sMax,
  1802.                         final DSSTThirdBodyDynamicContext context,
  1803.                         final double[] aoR3Pow, final double[][] qns) {
  1804.             this.nMax = nMax;
  1805.             this.sMax = sMax;

  1807.             final int rows    = nMax + 1;
  1808.             final int columns = sMax + 1;
  1809.             this.gns          = new double[rows][columns];
  1810.             this.dgnsda       = new double[rows][columns];
  1811.             this.dgnsdgamma   = new double[rows][columns];

  1813.             // Generate the coefficients
  1814.             generateCoefficients(context, aoR3Pow, qns);
  1815.         }
  1816.         /**
  1817.          * Compute the coefficient G<sub>n,s</sub> and its derivatives.
  1818.          * <p>
  1819.          * Only the derivatives by a and γ are computed as all others are 0
  1820.          * </p>
  1821.          * @param context container for attributes
  1822.          * @param aoR3Pow a / R3 up to power maxAR3Pow
  1823.          * @param qns Qns coefficients
  1824.          */
  1825.         private void generateCoefficients(final DSSTThirdBodyDynamicContext context, final double[] aoR3Pow, final double[][] qns) {

  1827.             final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();

  1829.             for (int s = 0; s <= sMax; s++) {
  1830.                 // The n index is always at least the maximum between 2 and s
  1831.                 final int minN = FastMath.max(2, s);
  1832.                 for (int n = minN; n <= nMax; n++) {
  1833.                     // compute the coefficients only if (n - s) % 2 == 0
  1834.                     if ( (n - s) % 2 == 0 ) {
  1835.                         // Kronecker symbol (2 - delta(0,s))
  1836.                         final double delta0s = (s == 0) ? 1. : 2.;
  1837.                         final double vns   = Vns.get(new NSKey(n, s));
  1838.                         final double coef0 = delta0s * aoR3Pow[n] * vns * context.getMuoR3();
  1839.                         final double coef1 = coef0 * qns[n][s];
  1840.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  1841.                         // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
  1842.                         final double dqns = (n == s) ? 0. : qns[n][s + 1];

  1844.                         //Compute the coefficient and its derivatives.
  1845.                         this.gns[n][s] = coef1;
  1846.                         this.dgnsda[n][s] = coef1 * n / auxiliaryElements.getSma();
  1847.                         this.dgnsdgamma[n][s] = coef0 * dqns;
  1848.                     } else {
  1849.                         // the coefficient and its derivatives is 0
  1850.                         this.gns[n][s] = 0.;
  1851.                         this.dgnsda[n][s] = 0.;
  1852.                         this.dgnsdgamma[n][s] = 0.;
  1853.                     }
  1854.                 }
  1855.             }
  1856.         }

  1858.         /** Get the coefficient G<sub>n,s</sub>.
  1859.          *
  1860.          * @param n n index
  1861.          * @param s s index
  1862.          * @return the coefficient G<sub>n,s</sub>
  1863.          */
  1864.         public double getGns(final int n, final int s) {
  1865.             return this.gns[n][s];
  1866.         }

  1868.         /** Get the derivative dG<sub>n,s</sub> / da.
  1869.          *
  1870.          * @param n n index
  1871.          * @param s s index
  1872.          * @return the derivative dG<sub>n,s</sub> / da
  1873.          */
  1874.         public double getdGnsda(final int n, final int s) {
  1875.             return this.dgnsda[n][s];
  1876.         }

  1878.         /** Get the derivative dG<sub>n,s</sub> / dγ.
  1879.          *
  1880.          * @param n n index
  1881.          * @param s s index
  1882.          * @return the derivative dG<sub>n,s</sub> / dγ
  1883.          */
  1884.         public double getdGnsdgamma(final int n, final int s) {
  1885.             return this.dgnsdgamma[n][s];
  1886.         }
  1887.     }

  1889.     /** The G<sub>n,s</sub> coefficients and their derivatives.
  1890.      * <p>
  1891.      * See Danielson, 4.2-17
  1892.      *
  1893.      * @author Lucian Barbulescu
  1894.      */
  1895.     private class FieldGnsCoefficients  <T extends CalculusFieldElement<T>> {

  1897.         /** Maximum value for n index. */
  1898.         private final int nMax;

  1900.         /** Maximum value for s index. */
  1901.         private final int sMax;

  1903.         /** The coefficients G<sub>n,s</sub>. */
  1904.         private final T gns[][];

  1906.         /** The derivatives of the coefficients G<sub>n,s</sub> by a. */
  1907.         private final T dgnsda[][];

  1909.         /** The derivatives of the coefficients G<sub>n,s</sub> by γ. */
  1910.         private final T dgnsdgamma[][];

  1912.         /** Standard constructor.
  1913.          *
  1914.          * @param nMax maximum value for n indes
  1915.          * @param sMax maximum value for s index
  1916.          * @param context container for attributes
  1917.          * @param aoR3Pow a / R3 up to power maxAR3Pow
  1918.          * @param qns Qns coefficients
  1919.          * @param field field used by default
  1920.          */
  1921.         FieldGnsCoefficients(final int nMax, final int sMax,
  1922.                              final FieldDSSTThirdBodyDynamicContext<T> context,
  1923.                              final T[] aoR3Pow, final T[][] qns,
  1924.                              final Field<T> field) {
  1925.             this.nMax = nMax;
  1926.             this.sMax = sMax;

  1928.             final int rows    = nMax + 1;
  1929.             final int columns = sMax + 1;
  1930.             this.gns          = MathArrays.buildArray(field, rows, columns);
  1931.             this.dgnsda       = MathArrays.buildArray(field, rows, columns);
  1932.             this.dgnsdgamma   = MathArrays.buildArray(field, rows, columns);

  1934.             // Generate the coefficients
  1935.             generateCoefficients(context, aoR3Pow, qns, field);
  1936.         }
  1937.         /**
  1938.          * Compute the coefficient G<sub>n,s</sub> and its derivatives.
  1939.          * <p>
  1940.          * Only the derivatives by a and γ are computed as all others are 0
  1941.          * </p>
  1942.          * @param context container for attributes
  1943.          * @param aoR3Pow a / R3 up to power maxAR3Pow
  1944.          * @param qns Qns coefficients
  1945.          * @param field field used by default
  1946.          */
  1947.         private void generateCoefficients(final FieldDSSTThirdBodyDynamicContext<T> context,
  1948.                                           final T[] aoR3Pow, final T[][] qns,
  1949.                                           final Field<T> field) {

  1951.             //Zero
  1952.             final T zero = field.getZero();

  1954.             final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();

  1956.             for (int s = 0; s <= sMax; s++) {
  1957.                 // The n index is always at least the maximum between 2 and s
  1958.                 final int minN = FastMath.max(2, s);
  1959.                 for (int n = minN; n <= nMax; n++) {
  1960.                     // compute the coefficients only if (n - s) % 2 == 0
  1961.                     if ( (n - s) % 2 == 0 ) {
  1962.                         // Kronecker symbol (2 - delta(0,s))
  1963.                         final T delta0s = (s == 0) ? zero.add(1.) : zero.add(2.);
  1964.                         final double vns = Vns.get(new NSKey(n, s));
  1965.                         final T coef0 = aoR3Pow[n].multiply(vns).multiply(context.getMuoR3()).multiply(delta0s);
  1966.                         final T coef1 = coef0.multiply(qns[n][s]);
  1967.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  1968.                         // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
  1969.                         final T dqns = (n == s) ? zero : qns[n][s + 1];

  1971.                         //Compute the coefficient and its derivatives.
  1972.                         this.gns[n][s] = coef1;
  1973.                         this.dgnsda[n][s] = coef1.multiply(n).divide(auxiliaryElements.getSma());
  1974.                         this.dgnsdgamma[n][s] = coef0.multiply(dqns);
  1975.                     } else {
  1976.                         // the coefficient and its derivatives is 0
  1977.                         this.gns[n][s] = zero;
  1978.                         this.dgnsda[n][s] = zero;
  1979.                         this.dgnsdgamma[n][s] = zero;
  1980.                     }
  1981.                 }
  1982.             }
  1983.         }

  1985.         /** Get the coefficient G<sub>n,s</sub>.
  1986.          *
  1987.          * @param n n index
  1988.          * @param s s index
  1989.          * @return the coefficient G<sub>n,s</sub>
  1990.          */
  1991.         public T getGns(final int n, final int s) {
  1992.             return this.gns[n][s];
  1993.         }

  1995.         /** Get the derivative dG<sub>n,s</sub> / da.
  1996.          *
  1997.          * @param n n index
  1998.          * @param s s index
  1999.          * @return the derivative dG<sub>n,s</sub> / da
  2000.          */
  2001.         public T getdGnsda(final int n, final int s) {
  2002.             return this.dgnsda[n][s];
  2003.         }

  2005.         /** Get the derivative dG<sub>n,s</sub> / dγ.
  2006.          *
  2007.          * @param n n index
  2008.          * @param s s index
  2009.          * @return the derivative dG<sub>n,s</sub> / dγ
  2010.          */
  2011.         public T getdGnsdgamma(final int n, final int s) {
  2012.             return this.dgnsdgamma[n][s];
  2013.         }
  2014.     }

  2016.     /** This class computes the terms containing the coefficients C<sub>j</sub> and S<sub>j</sub> of (α, β) or (k, h).
  2017.      *
  2018.      * <p>
  2019.      * The following terms and their derivatives by k, h, alpha and beta are considered: <br/ >
  2020.      * - 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 />
  2021.      * - C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) <br />
  2022.      * - 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 />
  2023.      * - C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) <br />
  2024.      * 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 />
  2025.      * See the CS Mathematical report $3.5.3.2 for more details
  2026.      * </p>
  2027.      * @author Lucian Barbulescu
  2028.      */
  2029.     private static class CjSjAlphaBetaKH {

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

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

  2037.         /** 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)
  2038.          * and its derivative by k, h, α and β. */
  2039.         private final double coefAandDeriv[];

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

  2045.         /** 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)
  2046.          * and its derivative by k, h, α and β. */
  2047.         private final double coefDandDeriv[];

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

  2053.         /**
  2054.          * Standard constructor.
  2055.          * @param context container for attributes
  2056.          */
  2057.         CjSjAlphaBetaKH(final DSSTThirdBodyDynamicContext context) {

  2059.             final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();

  2061.             cjsjkh = new CjSjCoefficient(auxiliaryElements.getK(), auxiliaryElements.getH());
  2062.             cjsjalbe = new CjSjCoefficient(context.getAlpha(), context.getBeta());

  2064.             coefAandDeriv = new double[5];
  2065.             coefBandDeriv = new double[5];
  2066.             coefDandDeriv = new double[5];
  2067.             coefEandDeriv = new double[5];
  2068.         }

  2070.         /** Compute the coefficients and their derivatives for a given (j,s) pair.
  2071.          * @param j j index
  2072.          * @param s s index
  2073.          */
  2074.         public void computeCoefficients(final int j, final int s) {
  2075.             // sign of j-s
  2076.             final int sign = j < s ? -1 : 1;

  2078.             //|j-s|
  2079.             final int absJmS = FastMath.abs(j - s);

  2081.             //j+s
  2082.             final int jps = j + s;

  2084.             //Compute the coefficient A and its derivatives
  2085.             coefAandDeriv[0] = sign * cjsjalbe.getCj(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getCj(absJmS);
  2086.             coefAandDeriv[1] = sign * cjsjalbe.getCj(s) * cjsjkh.getDsjDk(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getDcjDk(absJmS);
  2087.             coefAandDeriv[2] = sign * cjsjalbe.getCj(s) * cjsjkh.getDsjDh(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getDcjDh(absJmS);
  2088.             coefAandDeriv[3] = sign * cjsjalbe.getDcjDk(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getDsjDk(s) * cjsjkh.getCj(absJmS);
  2089.             coefAandDeriv[4] = sign * cjsjalbe.getDcjDh(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getDsjDh(s) * cjsjkh.getCj(absJmS);

  2091.             //Compute the coefficient B and its derivatives
  2092.             coefBandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getSj(jps) - cjsjalbe.getSj(s) * cjsjkh.getCj(jps);
  2093.             coefBandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDsjDk(jps) - cjsjalbe.getSj(s) * cjsjkh.getDcjDk(jps);
  2094.             coefBandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDsjDh(jps) - cjsjalbe.getSj(s) * cjsjkh.getDcjDh(jps);
  2095.             coefBandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getSj(jps) - cjsjalbe.getDsjDk(s) * cjsjkh.getCj(jps);
  2096.             coefBandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getSj(jps) - cjsjalbe.getDsjDh(s) * cjsjkh.getCj(jps);

  2098.             //Compute the coefficient D and its derivatives
  2099.             coefDandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getSj(absJmS);
  2100.             coefDandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDcjDk(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getDsjDk(absJmS);
  2101.             coefDandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDcjDh(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getDsjDh(absJmS);
  2102.             coefDandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getDsjDk(s) * cjsjkh.getSj(absJmS);
  2103.             coefDandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getDsjDh(s) * cjsjkh.getSj(absJmS);

  2105.             //Compute the coefficient E and its derivatives
  2106.             coefEandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getCj(jps) + cjsjalbe.getSj(s) * cjsjkh.getSj(jps);
  2107.             coefEandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDcjDk(jps) + cjsjalbe.getSj(s) * cjsjkh.getDsjDk(jps);
  2108.             coefEandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDcjDh(jps) + cjsjalbe.getSj(s) * cjsjkh.getDsjDh(jps);
  2109.             coefEandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getCj(jps) + cjsjalbe.getDsjDk(s) * cjsjkh.getSj(jps);
  2110.             coefEandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getCj(jps) + cjsjalbe.getDsjDh(s) * cjsjkh.getSj(jps);
  2111.         }

  2113.         /** Get the value of coefficient A<sub>j,s</sub>.
  2114.          *
  2115.          * @return the coefficient A<sub>j,s</sub>
  2116.          */
  2117.         public double getCoefA() {
  2118.             return coefAandDeriv[0];
  2119.         }

  2121.         /** Get the value of coefficient dA<sub>j,s</sub>/dk.
  2122.          *
  2123.          * @return the coefficient dA<sub>j,s</sub>/dk
  2124.          */
  2125.         public double getdCoefAdk() {
  2126.             return coefAandDeriv[1];
  2127.         }

  2129.         /** Get the value of coefficient dA<sub>j,s</sub>/dh.
  2130.          *
  2131.          * @return the coefficient dA<sub>j,s</sub>/dh
  2132.          */
  2133.         public double getdCoefAdh() {
  2134.             return coefAandDeriv[2];
  2135.         }

  2137.         /** Get the value of coefficient dA<sub>j,s</sub>/dα.
  2138.          *
  2139.          * @return the coefficient dA<sub>j,s</sub>/dα
  2140.          */
  2141.         public double getdCoefAdalpha() {
  2142.             return coefAandDeriv[3];
  2143.         }

  2145.         /** Get the value of coefficient dA<sub>j,s</sub>/dβ.
  2146.          *
  2147.          * @return the coefficient dA<sub>j,s</sub>/dβ
  2148.          */
  2149.         public double getdCoefAdbeta() {
  2150.             return coefAandDeriv[4];
  2151.         }

  2153.         /** Get the value of coefficient B<sub>j,s</sub>.
  2154.          *
  2155.          * @return the coefficient B<sub>j,s</sub>
  2156.          */
  2157.         public double getCoefB() {
  2158.             return coefBandDeriv[0];
  2159.         }

  2161.         /** Get the value of coefficient dB<sub>j,s</sub>/dk.
  2162.          *
  2163.          * @return the coefficient dB<sub>j,s</sub>/dk
  2164.          */
  2165.         public double getdCoefBdk() {
  2166.             return coefBandDeriv[1];
  2167.         }

  2169.         /** Get the value of coefficient dB<sub>j,s</sub>/dh.
  2170.          *
  2171.          * @return the coefficient dB<sub>j,s</sub>/dh
  2172.          */
  2173.         public double getdCoefBdh() {
  2174.             return coefBandDeriv[2];
  2175.         }

  2177.         /** Get the value of coefficient dB<sub>j,s</sub>/dα.
  2178.          *
  2179.          * @return the coefficient dB<sub>j,s</sub>/dα
  2180.          */
  2181.         public double getdCoefBdalpha() {
  2182.             return coefBandDeriv[3];
  2183.         }

  2185.         /** Get the value of coefficient dB<sub>j,s</sub>/dβ.
  2186.          *
  2187.          * @return the coefficient dB<sub>j,s</sub>/dβ
  2188.          */
  2189.         public double getdCoefBdbeta() {
  2190.             return coefBandDeriv[4];
  2191.         }

  2193.         /** Get the value of coefficient D<sub>j,s</sub>.
  2194.          *
  2195.          * @return the coefficient D<sub>j,s</sub>
  2196.          */
  2197.         public double getCoefD() {
  2198.             return coefDandDeriv[0];
  2199.         }

  2201.         /** Get the value of coefficient dD<sub>j,s</sub>/dk.
  2202.          *
  2203.          * @return the coefficient dD<sub>j,s</sub>/dk
  2204.          */
  2205.         public double getdCoefDdk() {
  2206.             return coefDandDeriv[1];
  2207.         }

  2209.         /** Get the value of coefficient dD<sub>j,s</sub>/dh.
  2210.          *
  2211.          * @return the coefficient dD<sub>j,s</sub>/dh
  2212.          */
  2213.         public double getdCoefDdh() {
  2214.             return coefDandDeriv[2];
  2215.         }

  2217.         /** Get the value of coefficient dD<sub>j,s</sub>/dα.
  2218.          *
  2219.          * @return the coefficient dD<sub>j,s</sub>/dα
  2220.          */
  2221.         public double getdCoefDdalpha() {
  2222.             return coefDandDeriv[3];
  2223.         }

  2225.         /** Get the value of coefficient dD<sub>j,s</sub>/dβ.
  2226.          *
  2227.          * @return the coefficient dD<sub>j,s</sub>/dβ
  2228.          */
  2229.         public double getdCoefDdbeta() {
  2230.             return coefDandDeriv[4];
  2231.         }

  2233.         /** Get the value of coefficient E<sub>j,s</sub>.
  2234.          *
  2235.          * @return the coefficient E<sub>j,s</sub>
  2236.          */
  2237.         public double getCoefE() {
  2238.             return coefEandDeriv[0];
  2239.         }

  2241.         /** Get the value of coefficient dE<sub>j,s</sub>/dk.
  2242.          *
  2243.          * @return the coefficient dE<sub>j,s</sub>/dk
  2244.          */
  2245.         public double getdCoefEdk() {
  2246.             return coefEandDeriv[1];
  2247.         }

  2249.         /** Get the value of coefficient dE<sub>j,s</sub>/dh.
  2250.          *
  2251.          * @return the coefficient dE<sub>j,s</sub>/dh
  2252.          */
  2253.         public double getdCoefEdh() {
  2254.             return coefEandDeriv[2];
  2255.         }

  2257.         /** Get the value of coefficient dE<sub>j,s</sub>/dα.
  2258.          *
  2259.          * @return the coefficient dE<sub>j,s</sub>/dα
  2260.          */
  2261.         public double getdCoefEdalpha() {
  2262.             return coefEandDeriv[3];
  2263.         }

  2265.         /** Get the value of coefficient dE<sub>j,s</sub>/dβ.
  2266.          *
  2267.          * @return the coefficient dE<sub>j,s</sub>/dβ
  2268.          */
  2269.         public double getdCoefEdbeta() {
  2270.             return coefEandDeriv[4];
  2271.         }
  2272.     }

  2274.      /** This class computes the terms containing the coefficients C<sub>j</sub> and S<sub>j</sub> of (α, β) or (k, h).
  2275.      *
  2276.      * <p>
  2277.      * The following terms and their derivatives by k, h, alpha and beta are considered: <br/ >
  2278.      * - 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 />
  2279.      * - C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) <br />
  2280.      * - 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 />
  2281.      * - C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) <br />
  2282.      * 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 />
  2283.      * See the CS Mathematical report $3.5.3.2 for more details
  2284.      * </p>
  2285.      * @author Lucian Barbulescu
  2286.      */
  2287.     private static class FieldCjSjAlphaBetaKH <T extends CalculusFieldElement<T>> {

  2289.         /** The C<sub>j</sub>(k, h) and the S<sub>j</sub>(k, h) series. */
  2290.         private final FieldCjSjCoefficient<T> cjsjkh;

  2292.         /** The C<sub>j</sub>(α, β) and the S<sub>j</sub>(α, β) series. */
  2293.         private final FieldCjSjCoefficient<T> cjsjalbe;

  2295.         /** 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)
  2296.          * and its derivative by k, h, α and β. */
  2297.         private final T coefAandDeriv[];

  2299.         /** The coeficient C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h)
  2300.          * and its derivative by k, h, α and β. */
  2301.         private final T coefBandDeriv[];

  2303.         /** 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)
  2304.          * and its derivative by k, h, α and β. */
  2305.         private final T coefDandDeriv[];

  2307.         /** The coeficient C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h)
  2308.          * and its derivative by k, h, α and β. */
  2309.         private final T coefEandDeriv[];

  2311.         /**
  2312.          * Standard constructor.
  2313.          * @param context container for attributes
  2314.          * @param field field used by default
  2315.          */
  2316.         FieldCjSjAlphaBetaKH(final FieldDSSTThirdBodyDynamicContext<T> context, final Field<T> field) {

  2318.             final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();

  2320.             cjsjkh   = new FieldCjSjCoefficient<>(auxiliaryElements.getK(), auxiliaryElements.getH(), field);
  2321.             cjsjalbe = new FieldCjSjCoefficient<>(context.getAlpha(), context.getBeta(), field);

  2323.             coefAandDeriv = MathArrays.buildArray(field, 5);
  2324.             coefBandDeriv = MathArrays.buildArray(field, 5);
  2325.             coefDandDeriv = MathArrays.buildArray(field, 5);
  2326.             coefEandDeriv = MathArrays.buildArray(field, 5);
  2327.         }

  2329.         /** Compute the coefficients and their derivatives for a given (j,s) pair.
  2330.          * @param j j index
  2331.          * @param s s index
  2332.          */
  2333.         public void computeCoefficients(final int j, final int s) {
  2334.             // sign of j-s
  2335.             final int sign = j < s ? -1 : 1;

  2337.             //|j-s|
  2338.             final int absJmS = FastMath.abs(j - s);

  2340.             //j+s
  2341.             final int jps = j + s;

  2343.             //Compute the coefficient A and its derivatives
  2344.             coefAandDeriv[0] = cjsjalbe.getCj(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign).add(cjsjalbe.getSj(s).multiply(cjsjkh.getCj(absJmS)));
  2345.             coefAandDeriv[1] = cjsjalbe.getCj(s).multiply(cjsjkh.getDsjDk(absJmS)).multiply(sign).add(cjsjalbe.getSj(s).multiply(cjsjkh.getDcjDk(absJmS)));
  2346.             coefAandDeriv[2] = cjsjalbe.getCj(s).multiply(cjsjkh.getDsjDh(absJmS)).multiply(sign).add(cjsjalbe.getSj(s).multiply(cjsjkh.getDcjDh(absJmS)));
  2347.             coefAandDeriv[3] = cjsjalbe.getDcjDk(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign).add(cjsjalbe.getDsjDk(s).multiply(cjsjkh.getCj(absJmS)));
  2348.             coefAandDeriv[4] = cjsjalbe.getDcjDh(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign).add(cjsjalbe.getDsjDh(s).multiply(cjsjkh.getCj(absJmS)));

  2350.             //Compute the coefficient B and its derivatives
  2351.             coefBandDeriv[0] = cjsjalbe.getCj(s).multiply(cjsjkh.getSj(jps)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getCj(jps)));
  2352.             coefBandDeriv[1] = cjsjalbe.getCj(s).multiply(cjsjkh.getDsjDk(jps)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getDcjDk(jps)));
  2353.             coefBandDeriv[2] = cjsjalbe.getCj(s).multiply(cjsjkh.getDsjDh(jps)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getDcjDh(jps)));
  2354.             coefBandDeriv[3] = cjsjalbe.getDcjDk(s).multiply(cjsjkh.getSj(jps)).subtract(cjsjalbe.getDsjDk(s).multiply(cjsjkh.getCj(jps)));
  2355.             coefBandDeriv[4] = cjsjalbe.getDcjDh(s).multiply(cjsjkh.getSj(jps)).subtract(cjsjalbe.getDsjDh(s).multiply(cjsjkh.getCj(jps)));

  2357.             //Compute the coefficient D and its derivatives
  2358.             coefDandDeriv[0] = cjsjalbe.getCj(s).multiply(cjsjkh.getCj(absJmS)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign));
  2359.             coefDandDeriv[1] = cjsjalbe.getCj(s).multiply(cjsjkh.getDcjDk(absJmS)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getDsjDk(absJmS)).multiply(sign));
  2360.             coefDandDeriv[2] = cjsjalbe.getCj(s).multiply(cjsjkh.getDcjDh(absJmS)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getDsjDh(absJmS)).multiply(sign));
  2361.             coefDandDeriv[3] = cjsjalbe.getDcjDk(s).multiply(cjsjkh.getCj(absJmS)).subtract(cjsjalbe.getDsjDk(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign));
  2362.             coefDandDeriv[4] = cjsjalbe.getDcjDh(s).multiply(cjsjkh.getCj(absJmS)).subtract(cjsjalbe.getDsjDh(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign));

  2364.             //Compute the coefficient E and its derivatives
  2365.             coefEandDeriv[0] = cjsjalbe.getCj(s).multiply(cjsjkh.getCj(jps)).add(cjsjalbe.getSj(s).multiply(cjsjkh.getSj(jps)));
  2366.             coefEandDeriv[1] = cjsjalbe.getCj(s).multiply(cjsjkh.getDcjDk(jps)).add(cjsjalbe.getSj(s).multiply(cjsjkh.getDsjDk(jps)));
  2367.             coefEandDeriv[2] = cjsjalbe.getCj(s).multiply(cjsjkh.getDcjDh(jps)).add(cjsjalbe.getSj(s).multiply(cjsjkh.getDsjDh(jps)));
  2368.             coefEandDeriv[3] = cjsjalbe.getDcjDk(s).multiply(cjsjkh.getCj(jps)).add(cjsjalbe.getDsjDk(s).multiply(cjsjkh.getSj(jps)));
  2369.             coefEandDeriv[4] = cjsjalbe.getDcjDh(s).multiply(cjsjkh.getCj(jps)).add(cjsjalbe.getDsjDh(s).multiply(cjsjkh.getSj(jps)));
  2370.         }

  2372.         /** Get the value of coefficient A<sub>j,s</sub>.
  2373.          *
  2374.          * @return the coefficient A<sub>j,s</sub>
  2375.          */
  2376.         public T getCoefA() {
  2377.             return coefAandDeriv[0];
  2378.         }

  2380.         /** Get the value of coefficient dA<sub>j,s</sub>/dk.
  2381.          *
  2382.          * @return the coefficient dA<sub>j,s</sub>/dk
  2383.          */
  2384.         public T getdCoefAdk() {
  2385.             return coefAandDeriv[1];
  2386.         }

  2388.         /** Get the value of coefficient dA<sub>j,s</sub>/dh.
  2389.          *
  2390.          * @return the coefficient dA<sub>j,s</sub>/dh
  2391.          */
  2392.         public T getdCoefAdh() {
  2393.             return coefAandDeriv[2];
  2394.         }

  2396.         /** Get the value of coefficient dA<sub>j,s</sub>/dα.
  2397.          *
  2398.          * @return the coefficient dA<sub>j,s</sub>/dα
  2399.          */
  2400.         public T getdCoefAdalpha() {
  2401.             return coefAandDeriv[3];
  2402.         }

  2404.         /** Get the value of coefficient dA<sub>j,s</sub>/dβ.
  2405.          *
  2406.          * @return the coefficient dA<sub>j,s</sub>/dβ
  2407.          */
  2408.         public T getdCoefAdbeta() {
  2409.             return coefAandDeriv[4];
  2410.         }

  2412.        /** Get the value of coefficient B<sub>j,s</sub>.
  2413.         *
  2414.         * @return the coefficient B<sub>j,s</sub>
  2415.         */
  2416.         public T getCoefB() {
  2417.             return coefBandDeriv[0];
  2418.         }

  2420.         /** Get the value of coefficient dB<sub>j,s</sub>/dk.
  2421.          *
  2422.          * @return the coefficient dB<sub>j,s</sub>/dk
  2423.          */
  2424.         public T getdCoefBdk() {
  2425.             return coefBandDeriv[1];
  2426.         }

  2428.         /** Get the value of coefficient dB<sub>j,s</sub>/dh.
  2429.          *
  2430.          * @return the coefficient dB<sub>j,s</sub>/dh
  2431.          */
  2432.         public T getdCoefBdh() {
  2433.             return coefBandDeriv[2];
  2434.         }

  2436.         /** Get the value of coefficient dB<sub>j,s</sub>/dα.
  2437.          *
  2438.          * @return the coefficient dB<sub>j,s</sub>/dα
  2439.          */
  2440.         public T getdCoefBdalpha() {
  2441.             return coefBandDeriv[3];
  2442.         }

  2444.         /** Get the value of coefficient dB<sub>j,s</sub>/dβ.
  2445.          *
  2446.          * @return the coefficient dB<sub>j,s</sub>/dβ
  2447.          */
  2448.         public T getdCoefBdbeta() {
  2449.             return coefBandDeriv[4];
  2450.         }

  2452.         /** Get the value of coefficient D<sub>j,s</sub>.
  2453.          *
  2454.          * @return the coefficient D<sub>j,s</sub>
  2455.          */
  2456.         public T getCoefD() {
  2457.             return coefDandDeriv[0];
  2458.         }

  2460.         /** Get the value of coefficient dD<sub>j,s</sub>/dk.
  2461.          *
  2462.          * @return the coefficient dD<sub>j,s</sub>/dk
  2463.          */
  2464.         public T getdCoefDdk() {
  2465.             return coefDandDeriv[1];
  2466.         }

  2468.         /** Get the value of coefficient dD<sub>j,s</sub>/dh.
  2469.          *
  2470.          * @return the coefficient dD<sub>j,s</sub>/dh
  2471.          */
  2472.         public T getdCoefDdh() {
  2473.             return coefDandDeriv[2];
  2474.         }

  2476.         /** Get the value of coefficient dD<sub>j,s</sub>/dα.
  2477.          *
  2478.          * @return the coefficient dD<sub>j,s</sub>/dα
  2479.          */
  2480.         public T getdCoefDdalpha() {
  2481.             return coefDandDeriv[3];
  2482.         }

  2484.         /** Get the value of coefficient dD<sub>j,s</sub>/dβ.
  2485.          *
  2486.          * @return the coefficient dD<sub>j,s</sub>/dβ
  2487.          */
  2488.         public T getdCoefDdbeta() {
  2489.             return coefDandDeriv[4];
  2490.         }

  2492.         /** Get the value of coefficient E<sub>j,s</sub>.
  2493.          *
  2494.          * @return the coefficient E<sub>j,s</sub>
  2495.          */
  2496.         public T getCoefE() {
  2497.             return coefEandDeriv[0];
  2498.         }

  2500.         /** Get the value of coefficient dE<sub>j,s</sub>/dk.
  2501.          *
  2502.          * @return the coefficient dE<sub>j,s</sub>/dk
  2503.          */
  2504.         public T getdCoefEdk() {
  2505.             return coefEandDeriv[1];
  2506.         }

  2508.         /** Get the value of coefficient dE<sub>j,s</sub>/dh.
  2509.          *
  2510.          * @return the coefficient dE<sub>j,s</sub>/dh
  2511.          */
  2512.         public T getdCoefEdh() {
  2513.             return coefEandDeriv[2];
  2514.         }

  2516.         /** Get the value of coefficient dE<sub>j,s</sub>/dα.
  2517.          *
  2518.          * @return the coefficient dE<sub>j,s</sub>/dα
  2519.          */
  2520.         public T getdCoefEdalpha() {
  2521.             return coefEandDeriv[3];
  2522.         }

  2524.         /** Get the value of coefficient dE<sub>j,s</sub>/dβ.
  2525.          *
  2526.          * @return the coefficient dE<sub>j,s</sub>/dβ
  2527.          */
  2528.         public T getdCoefEdbeta() {
  2529.             return coefEandDeriv[4];
  2530.         }
  2531.     }

  2533.     /** This class computes the coefficients for the generating function S and its derivatives.
  2534.      * <p>
  2535.      * The form of the generating functions is: <br>
  2536.      *  S = C⁰ + &Sigma;<sub>j=1</sub><sup>N+1</sup>(C<sup>j</sup> * cos(jF) + S<sup>j</sup> * sin(jF)) <br>
  2537.      *  The coefficients C⁰, C<sup>j</sup>, S<sup>j</sup> are the Fourrier coefficients
  2538.      *  presented in Danielson 4.2-14,15 except for the case j=1 where
  2539.      *  C¹ = C¹<sub>Fourier</sub> - hU and
  2540.      *  S¹ = S¹<sub>Fourier</sub> + kU <br>
  2541.      *  Also the coefficients of the derivatives of S by a, k, h, α, β, γ and λ
  2542.      *  are computed end expressed in a similar manner. The formulas used are 4.2-19, 20, 23, 24
  2543.      * </p>
  2544.      * @author Lucian Barbulescu
  2545.      */
  2546.     private class GeneratingFunctionCoefficients {

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

  2551.         /** Maximum value of j index. */
  2552.         private final int jMax;

  2554.         /** The coefficients C<sup>j</sup> of the function S and its derivatives.
  2555.          * <p>
  2556.          * The index j belongs to the interval [0,jMax]. The coefficient C⁰ is the free coefficient.<br>
  2557.          * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
  2558.          * - S <br/>
  2559.          * - dS / da <br/>
  2560.          * - dS / dk <br/>
  2561.          * - dS / dh <br/>
  2562.          * - dS / dα <br/>
  2563.          * - dS / dβ <br/>
  2564.          * - dS / dγ <br/>
  2565.          * - dS / dλ
  2566.          * </p>
  2567.          */
  2568.         private final double[][] cjCoefs;

  2570.         /** The coefficients S<sup>j</sup> of the function S and its derivatives.
  2571.          * <p>
  2572.          * The index j belongs to the interval [0,jMax].<br>
  2573.          * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
  2574.          * - S <br/>
  2575.          * - dS / da <br/>
  2576.          * - dS / dk <br/>
  2577.          * - dS / dh <br/>
  2578.          * - dS / dα <br/>
  2579.          * - dS / dβ <br/>
  2580.          * - dS / dγ <br/>
  2581.          * - dS / dλ
  2582.          * </p>
  2583.          */
  2584.         private final double[][] sjCoefs;

  2586.         /**
  2587.          * Standard constructor.
  2588.          *
  2589.          * @param nMax maximum value of n index
  2590.          * @param sMax maximum value of s index
  2591.          * @param jMax maximum value of j index
  2592.          * @param context container for attributes
  2593.          * @param hansen hansen objects
  2594.          * @param aoR3Pow a / R3 up to power maxAR3Pow
  2595.          * @param qns Qns coefficients
  2596.          */
  2597.         GeneratingFunctionCoefficients(final int nMax, final int sMax, final int jMax,
  2598.                                        final DSSTThirdBodyDynamicContext context,
  2599.                                        final HansenObjects hansen,
  2600.                                        final double[] aoR3Pow, final double[][] qns) {
  2601.             this.jMax = jMax;
  2602.             this.cjsjFourier = new FourierCjSjCoefficients(nMax, sMax, jMax, context, aoR3Pow, qns);
  2603.             this.cjCoefs = new double[8][jMax + 1];
  2604.             this.sjCoefs = new double[8][jMax + 1];

  2606.             computeGeneratingFunctionCoefficients(context, hansen, aoR3Pow, qns);
  2607.         }

  2609.         /**
  2610.          * Compute the coefficients for the generating function S and its derivatives.
  2611.          * @param context container for attributes
  2612.          * @param hansenObjects hansen objects
  2613.          * @param aoR3Pow a / R3 up to power maxAR3Pow
  2614.          * @param qns Qns coefficients
  2615.          */
  2616.         private void computeGeneratingFunctionCoefficients(final DSSTThirdBodyDynamicContext context, final HansenObjects hansenObjects,
  2617.                                                            final double[] aoR3Pow, final double[][] qns) {

  2619.             final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();

  2621.             // Access to potential U derivatives
  2622.             final UAnddU udu = new UAnddU(context, hansenObjects, aoR3Pow, qns);

  2624.             //Compute the C<sup>j</sup> coefficients
  2625.             for (int j = 1; j <= jMax; j++) {
  2626.                 //Compute the C<sup>j</sup> coefficients
  2627.                 cjCoefs[0][j] = cjsjFourier.getCj(j);
  2628.                 cjCoefs[1][j] = cjsjFourier.getdCjda(j);
  2629.                 cjCoefs[2][j] = cjsjFourier.getdCjdk(j) - (cjsjFourier.getSjLambda(j - 1) - cjsjFourier.getSjLambda(j + 1)) / 2;
  2630.                 cjCoefs[3][j] = cjsjFourier.getdCjdh(j) - (cjsjFourier.getCjLambda(j - 1) + cjsjFourier.getCjLambda(j + 1)) / 2;
  2631.                 cjCoefs[4][j] = cjsjFourier.getdCjdalpha(j);
  2632.                 cjCoefs[5][j] = cjsjFourier.getdCjdbeta(j);
  2633.                 cjCoefs[6][j] = cjsjFourier.getdCjdgamma(j);
  2634.                 cjCoefs[7][j] = cjsjFourier.getCjLambda(j);

  2636.                 //Compute the S<sup>j</sup> coefficients
  2637.                 sjCoefs[0][j] = cjsjFourier.getSj(j);
  2638.                 sjCoefs[1][j] = cjsjFourier.getdSjda(j);
  2639.                 sjCoefs[2][j] = cjsjFourier.getdSjdk(j) + (cjsjFourier.getCjLambda(j - 1) - cjsjFourier.getCjLambda(j + 1)) / 2;
  2640.                 sjCoefs[3][j] = cjsjFourier.getdSjdh(j) - (cjsjFourier.getSjLambda(j - 1) + cjsjFourier.getSjLambda(j + 1)) / 2;
  2641.                 sjCoefs[4][j] = cjsjFourier.getdSjdalpha(j);
  2642.                 sjCoefs[5][j] = cjsjFourier.getdSjdbeta(j);
  2643.                 sjCoefs[6][j] = cjsjFourier.getdSjdgamma(j);
  2644.                 sjCoefs[7][j] = cjsjFourier.getSjLambda(j);

  2646.                 //In the special case j == 1 there are some additional terms to be added
  2647.                 if (j == 1) {
  2648.                     //Additional terms for C<sup>j</sup> coefficients
  2649.                     cjCoefs[0][j] += -auxiliaryElements.getH() * udu.getU();
  2650.                     cjCoefs[1][j] += -auxiliaryElements.getH() * udu.getdUda();
  2651.                     cjCoefs[2][j] += -auxiliaryElements.getH() * udu.getdUdk();
  2652.                     cjCoefs[3][j] += -(auxiliaryElements.getH() * udu.getdUdh() + udu.getU() + cjsjFourier.getC0Lambda());
  2653.                     cjCoefs[4][j] += -auxiliaryElements.getH() * udu.getdUdAl();
  2654.                     cjCoefs[5][j] += -auxiliaryElements.getH() * udu.getdUdBe();
  2655.                     cjCoefs[6][j] += -auxiliaryElements.getH() * udu.getdUdGa();

  2657.                     //Additional terms for S<sup>j</sup> coefficients
  2658.                     sjCoefs[0][j] += auxiliaryElements.getK() * udu.getU();
  2659.                     sjCoefs[1][j] += auxiliaryElements.getK() * udu.getdUda();
  2660.                     sjCoefs[2][j] += auxiliaryElements.getK() * udu.getdUdk() + udu.getU() + cjsjFourier.getC0Lambda();
  2661.                     sjCoefs[3][j] += auxiliaryElements.getK() * udu.getdUdh();
  2662.                     sjCoefs[4][j] += auxiliaryElements.getK() * udu.getdUdAl();
  2663.                     sjCoefs[5][j] += auxiliaryElements.getK() * udu.getdUdBe();
  2664.                     sjCoefs[6][j] += auxiliaryElements.getK() * udu.getdUdGa();
  2665.                 }
  2666.             }
  2667.         }

  2669.         /** Get the coefficient C<sup>j</sup> for the function S.
  2670.          * <br>
  2671.          * Possible values for j are within the interval [0,jMax].
  2672.          * The value 0 is used to obtain the free coefficient C⁰
  2673.          * @param j j index
  2674.          * @return C<sup>j</sup> for the function S
  2675.          */
  2676.         public double getSCj(final int j) {
  2677.             return cjCoefs[0][j];
  2678.         }

  2680.         /** Get the coefficient S<sup>j</sup> for the function S.
  2681.          * <br>
  2682.          * Possible values for j are within the interval [1,jMax].
  2683.          * @param j j index
  2684.          * @return S<sup>j</sup> for the function S
  2685.          */
  2686.         public double getSSj(final int j) {
  2687.             return sjCoefs[0][j];
  2688.         }

  2690.         /** Get the coefficient C<sup>j</sup> for the derivative dS/da.
  2691.          * <br>
  2692.          * Possible values for j are within the interval [0,jMax].
  2693.          * The value 0 is used to obtain the free coefficient C⁰
  2694.          * @param j j index
  2695.          * @return C<sup>j</sup> for the function dS/da
  2696.          */
  2697.         public double getdSdaCj(final int j) {
  2698.             return cjCoefs[1][j];
  2699.         }

  2701.         /** Get the coefficient S<sup>j</sup> for the derivative dS/da.
  2702.          * <br>
  2703.          * Possible values for j are within the interval [1,jMax].
  2704.          * @param j j index
  2705.          * @return S<sup>j</sup> for the derivative dS/da
  2706.          */
  2707.         public double getdSdaSj(final int j) {
  2708.             return sjCoefs[1][j];
  2709.         }

  2711.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dk
  2712.          * <br>
  2713.          * Possible values for j are within the interval [0,jMax].
  2714.          * The value 0 is used to obtain the free coefficient C⁰
  2715.          * @param j j index
  2716.          * @return C<sup>j</sup> for the function dS/dk
  2717.          */
  2718.         public double getdSdkCj(final int j) {
  2719.             return cjCoefs[2][j];
  2720.         }

  2722.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dk.
  2723.          * <br>
  2724.          * Possible values for j are within the interval [1,jMax].
  2725.          * @param j j index
  2726.          * @return S<sup>j</sup> for the derivative dS/dk
  2727.          */
  2728.         public double getdSdkSj(final int j) {
  2729.             return sjCoefs[2][j];
  2730.         }

  2732.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dh
  2733.          * <br>
  2734.          * Possible values for j are within the interval [0,jMax].
  2735.          * The value 0 is used to obtain the free coefficient C⁰
  2736.          * @param j j index
  2737.          * @return C<sup>j</sup> for the function dS/dh
  2738.          */
  2739.         public double getdSdhCj(final int j) {
  2740.             return cjCoefs[3][j];
  2741.         }

  2743.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dh.
  2744.          * <br>
  2745.          * Possible values for j are within the interval [1,jMax].
  2746.          * @param j j index
  2747.          * @return S<sup>j</sup> for the derivative dS/dh
  2748.          */
  2749.         public double getdSdhSj(final int j) {
  2750.             return sjCoefs[3][j];
  2751.         }

  2753.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dα
  2754.          * <br>
  2755.          * Possible values for j are within the interval [0,jMax].
  2756.          * The value 0 is used to obtain the free coefficient C⁰
  2757.          * @param j j index
  2758.          * @return C<sup>j</sup> for the function dS/dα
  2759.          */
  2760.         public double getdSdalphaCj(final int j) {
  2761.             return cjCoefs[4][j];
  2762.         }

  2764.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dα.
  2765.          * <br>
  2766.          * Possible values for j are within the interval [1,jMax].
  2767.          * @param j j index
  2768.          * @return S<sup>j</sup> for the derivative dS/dα
  2769.          */
  2770.         public double getdSdalphaSj(final int j) {
  2771.             return sjCoefs[4][j];
  2772.         }

  2774.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dβ
  2775.          * <br>
  2776.          * Possible values for j are within the interval [0,jMax].
  2777.          * The value 0 is used to obtain the free coefficient C⁰
  2778.          * @param j j index
  2779.          * @return C<sup>j</sup> for the function dS/dβ
  2780.          */
  2781.         public double getdSdbetaCj(final int j) {
  2782.             return cjCoefs[5][j];
  2783.         }

  2785.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dβ.
  2786.          * <br>
  2787.          * Possible values for j are within the interval [1,jMax].
  2788.          * @param j j index
  2789.          * @return S<sup>j</sup> for the derivative dS/dβ
  2790.          */
  2791.         public double getdSdbetaSj(final int j) {
  2792.             return sjCoefs[5][j];
  2793.         }

  2795.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dγ
  2796.          * <br>
  2797.          * Possible values for j are within the interval [0,jMax].
  2798.          * The value 0 is used to obtain the free coefficient C⁰
  2799.          * @param j j index
  2800.          * @return C<sup>j</sup> for the function dS/dγ
  2801.          */
  2802.         public double getdSdgammaCj(final int j) {
  2803.             return cjCoefs[6][j];
  2804.         }

  2806.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dγ.
  2807.          * <br>
  2808.          * Possible values for j are within the interval [1,jMax].
  2809.          * @param j j index
  2810.          * @return S<sup>j</sup> for the derivative dS/dγ
  2811.          */
  2812.         public double getdSdgammaSj(final int j) {
  2813.             return sjCoefs[6][j];
  2814.         }

  2816.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dλ
  2817.          * <br>
  2818.          * Possible values for j are within the interval [0,jMax].
  2819.          * The value 0 is used to obtain the free coefficient C⁰
  2820.          * @param j j index
  2821.          * @return C<sup>j</sup> for the function dS/dλ
  2822.          */
  2823.         public double getdSdlambdaCj(final int j) {
  2824.             return cjCoefs[7][j];
  2825.         }

  2827.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dλ.
  2828.          * <br>
  2829.          * Possible values for j are within the interval [1,jMax].
  2830.          * @param j j index
  2831.          * @return S<sup>j</sup> for the derivative dS/dλ
  2832.          */
  2833.         public double getdSdlambdaSj(final int j) {
  2834.             return sjCoefs[7][j];
  2835.         }
  2836.     }

  2838.     /** This class computes the coefficients for the generating function S and its derivatives.
  2839.      * <p>
  2840.      * The form of the generating functions is: <br>
  2841.      *  S = C⁰ + &Sigma;<sub>j=1</sub><sup>N+1</sup>(C<sup>j</sup> * cos(jF) + S<sup>j</sup> * sin(jF)) <br>
  2842.      *  The coefficients C⁰, C<sup>j</sup>, S<sup>j</sup> are the Fourrier coefficients
  2843.      *  presented in Danielson 4.2-14,15 except for the case j=1 where
  2844.      *  C¹ = C¹<sub>Fourier</sub> - hU and
  2845.      *  S¹ = S¹<sub>Fourier</sub> + kU <br>
  2846.      *  Also the coefficients of the derivatives of S by a, k, h, α, β, γ and λ
  2847.      *  are computed end expressed in a similar manner. The formulas used are 4.2-19, 20, 23, 24
  2848.      * </p>
  2849.      * @author Lucian Barbulescu
  2850.      */
  2851.     private class FieldGeneratingFunctionCoefficients <T extends CalculusFieldElement<T>> {

  2853.         /** The Fourier coefficients as presented in Danielson 4.2-14,15. */
  2854.         private final FieldFourierCjSjCoefficients<T> cjsjFourier;

  2856.         /** Maximum value of j index. */
  2857.         private final int jMax;

  2859.         /** The coefficients C<sup>j</sup> of the function S and its derivatives.
  2860.          * <p>
  2861.          * The index j belongs to the interval [0,jMax]. The coefficient C⁰ is the free coefficient.<br>
  2862.          * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
  2863.          * - S <br/>
  2864.          * - dS / da <br/>
  2865.          * - dS / dk <br/>
  2866.          * - dS / dh <br/>
  2867.          * - dS / dα <br/>
  2868.          * - dS / dβ <br/>
  2869.          * - dS / dγ <br/>
  2870.          * - dS / dλ
  2871.          * </p>
  2872.          */
  2873.         private final T[][] cjCoefs;

  2875.         /** The coefficients S<sup>j</sup> of the function S and its derivatives.
  2876.          * <p>
  2877.          * The index j belongs to the interval [0,jMax].<br>
  2878.          * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
  2879.          * - S <br/>
  2880.          * - dS / da <br/>
  2881.          * - dS / dk <br/>
  2882.          * - dS / dh <br/>
  2883.          * - dS / dα <br/>
  2884.          * - dS / dβ <br/>
  2885.          * - dS / dγ <br/>
  2886.          * - dS / dλ
  2887.          * </p>
  2888.          */
  2889.         private final T[][] sjCoefs;

  2891.         /**
  2892.          * Standard constructor.
  2893.          *
  2894.          * @param nMax maximum value of n index
  2895.          * @param sMax maximum value of s index
  2896.          * @param jMax maximum value of j index
  2897.          * @param context container for attributes
  2898.          * @param hansen hansen objects
  2899.          * @param field field used by default
  2900.          * @param aoR3Pow a / R3 up to power maxAR3Pow
  2901.          * @param qns Qns coefficients
  2902.          */
  2903.         FieldGeneratingFunctionCoefficients(final int nMax, final int sMax, final int jMax,
  2904.                                             final FieldDSSTThirdBodyDynamicContext<T> context,
  2905.                                             final FieldHansenObjects<T> hansen, final Field<T> field,
  2906.                                             final T[] aoR3Pow, final T[][] qns) {
  2907.             this.jMax = jMax;
  2908.             this.cjsjFourier = new FieldFourierCjSjCoefficients<>(nMax, sMax, jMax, context, aoR3Pow, qns, field);
  2909.             this.cjCoefs     = MathArrays.buildArray(field, 8, jMax + 1);
  2910.             this.sjCoefs     = MathArrays.buildArray(field, 8, jMax + 1);

  2912.             computeGeneratingFunctionCoefficients(context, hansen, aoR3Pow, qns);
  2913.         }

  2915.         /**
  2916.          * Compute the coefficients for the generating function S and its derivatives.
  2917.          * @param context container for attributes
  2918.          * @param hansenObjects hansen objects
  2919.          * @param aoR3Pow a / R3 up to power maxAR3Pow
  2920.          * @param qns Qns coefficients
  2921.          */
  2922.         private void computeGeneratingFunctionCoefficients(final FieldDSSTThirdBodyDynamicContext<T> context,
  2923.                                                            final FieldHansenObjects<T> hansenObjects,
  2924.                                                            final T[] aoR3Pow, final T[][] qns) {

  2926.             final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();

  2928.             // Access to potential U derivatives
  2929.             final FieldUAnddU<T> udu = new FieldUAnddU<>(context, hansenObjects, aoR3Pow, qns);

  2931.             //Compute the C<sup>j</sup> coefficients
  2932.             for (int j = 1; j <= jMax; j++) {
  2933.                 //Compute the C<sup>j</sup> coefficients
  2934.                 cjCoefs[0][j] = cjsjFourier.getCj(j);
  2935.                 cjCoefs[1][j] = cjsjFourier.getdCjda(j);
  2936.                 cjCoefs[2][j] = cjsjFourier.getdCjdk(j).subtract((cjsjFourier.getSjLambda(j - 1).subtract(cjsjFourier.getSjLambda(j + 1))).divide(2.));
  2937.                 cjCoefs[3][j] = cjsjFourier.getdCjdh(j).subtract((cjsjFourier.getCjLambda(j - 1).add(cjsjFourier.getCjLambda(j + 1))).divide(2.));
  2938.                 cjCoefs[4][j] = cjsjFourier.getdCjdalpha(j);
  2939.                 cjCoefs[5][j] = cjsjFourier.getdCjdbeta(j);
  2940.                 cjCoefs[6][j] = cjsjFourier.getdCjdgamma(j);
  2941.                 cjCoefs[7][j] = cjsjFourier.getCjLambda(j);

  2943.                 //Compute the S<sup>j</sup> coefficients
  2944.                 sjCoefs[0][j] = cjsjFourier.getSj(j);
  2945.                 sjCoefs[1][j] = cjsjFourier.getdSjda(j);
  2946.                 sjCoefs[2][j] = cjsjFourier.getdSjdk(j).add((cjsjFourier.getCjLambda(j - 1).subtract(cjsjFourier.getCjLambda(j + 1))).divide(2.));
  2947.                 sjCoefs[3][j] = cjsjFourier.getdSjdh(j).subtract((cjsjFourier.getSjLambda(j - 1).add(cjsjFourier.getSjLambda(j + 1))).divide(2.));
  2948.                 sjCoefs[4][j] = cjsjFourier.getdSjdalpha(j);
  2949.                 sjCoefs[5][j] = cjsjFourier.getdSjdbeta(j);
  2950.                 sjCoefs[6][j] = cjsjFourier.getdSjdgamma(j);
  2951.                 sjCoefs[7][j] = cjsjFourier.getSjLambda(j);

  2953.                 //In the special case j == 1 there are some additional terms to be added
  2954.                 if (j == 1) {
  2955.                     //Additional terms for C<sup>j</sup> coefficients
  2956.                     cjCoefs[0][j] = cjCoefs[0][j].add(auxiliaryElements.getH().negate().multiply(udu.getU()));
  2957.                     cjCoefs[1][j] = cjCoefs[1][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUda()));
  2958.                     cjCoefs[2][j] = cjCoefs[2][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUdk()));
  2959.                     cjCoefs[3][j] = cjCoefs[3][j].add(auxiliaryElements.getH().multiply(udu.getdUdh()).add(udu.getU()).add(cjsjFourier.getC0Lambda()).negate());
  2960.                     cjCoefs[4][j] = cjCoefs[4][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUdAl()));
  2961.                     cjCoefs[5][j] = cjCoefs[5][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUdBe()));
  2962.                     cjCoefs[6][j] = cjCoefs[6][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUdGa()));

  2964.                     //Additional terms for S<sup>j</sup> coefficients
  2965.                     sjCoefs[0][j] = sjCoefs[0][j].add(auxiliaryElements.getK().multiply(udu.getU()));
  2966.                     sjCoefs[1][j] = sjCoefs[1][j].add(auxiliaryElements.getK().multiply(udu.getdUda()));
  2967.                     sjCoefs[2][j] = sjCoefs[2][j].add(auxiliaryElements.getK().multiply(udu.getdUdk()).add(udu.getU()).add(cjsjFourier.getC0Lambda()));
  2968.                     sjCoefs[3][j] = sjCoefs[3][j].add(auxiliaryElements.getK().multiply(udu.getdUdh()));
  2969.                     sjCoefs[4][j] = sjCoefs[4][j].add(auxiliaryElements.getK().multiply(udu.getdUdAl()));
  2970.                     sjCoefs[5][j] = sjCoefs[5][j].add(auxiliaryElements.getK().multiply(udu.getdUdBe()));
  2971.                     sjCoefs[6][j] = sjCoefs[6][j].add(auxiliaryElements.getK().multiply(udu.getdUdGa()));
  2972.                 }
  2973.             }
  2974.         }

  2976.         /** Get the coefficient C<sup>j</sup> for the function S.
  2977.          * <br>
  2978.          * Possible values for j are within the interval [0,jMax].
  2979.          * The value 0 is used to obtain the free coefficient C⁰
  2980.          * @param j j index
  2981.          * @return C<sup>j</sup> for the function S
  2982.          */
  2983.         public T getSCj(final int j) {
  2984.             return cjCoefs[0][j];
  2985.         }

  2987.         /** Get the coefficient S<sup>j</sup> for the function S.
  2988.          * <br>
  2989.          * Possible values for j are within the interval [1,jMax].
  2990.          * @param j j index
  2991.          * @return S<sup>j</sup> for the function S
  2992.          */
  2993.         public T getSSj(final int j) {
  2994.             return sjCoefs[0][j];
  2995.         }

  2997.         /** Get the coefficient C<sup>j</sup> for the derivative dS/da.
  2998.          * <br>
  2999.          * Possible values for j are within the interval [0,jMax].
  3000.          * The value 0 is used to obtain the free coefficient C⁰
  3001.          * @param j j index
  3002.          * @return C<sup>j</sup> for the function dS/da
  3003.          */
  3004.         public T getdSdaCj(final int j) {
  3005.             return cjCoefs[1][j];
  3006.         }

  3008.         /** Get the coefficient S<sup>j</sup> for the derivative dS/da.
  3009.          * <br>
  3010.          * Possible values for j are within the interval [1,jMax].
  3011.          * @param j j index
  3012.          * @return S<sup>j</sup> for the derivative dS/da
  3013.          */
  3014.         public T getdSdaSj(final int j) {
  3015.             return sjCoefs[1][j];
  3016.         }

  3018.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dk
  3019.          * <br>
  3020.          * Possible values for j are within the interval [0,jMax].
  3021.          * The value 0 is used to obtain the free coefficient C⁰
  3022.          * @param j j index
  3023.          * @return C<sup>j</sup> for the function dS/dk
  3024.          */
  3025.         public T getdSdkCj(final int j) {
  3026.             return cjCoefs[2][j];
  3027.         }

  3029.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dk.
  3030.          * <br>
  3031.          * Possible values for j are within the interval [1,jMax].
  3032.          * @param j j index
  3033.          * @return S<sup>j</sup> for the derivative dS/dk
  3034.          */
  3035.         public T getdSdkSj(final int j) {
  3036.             return sjCoefs[2][j];
  3037.         }

  3039.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dh
  3040.          * <br>
  3041.          * Possible values for j are within the interval [0,jMax].
  3042.          * The value 0 is used to obtain the free coefficient C⁰
  3043.          * @param j j index
  3044.          * @return C<sup>j</sup> for the function dS/dh
  3045.          */
  3046.         public T getdSdhCj(final int j) {
  3047.             return cjCoefs[3][j];
  3048.         }

  3050.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dh.
  3051.          * <br>
  3052.          * Possible values for j are within the interval [1,jMax].
  3053.          * @param j j index
  3054.          * @return S<sup>j</sup> for the derivative dS/dh
  3055.          */
  3056.         public T getdSdhSj(final int j) {
  3057.             return sjCoefs[3][j];
  3058.         }

  3060.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dα
  3061.          * <br>
  3062.          * Possible values for j are within the interval [0,jMax].
  3063.          * The value 0 is used to obtain the free coefficient C⁰
  3064.          * @param j j index
  3065.          * @return C<sup>j</sup> for the function dS/dα
  3066.          */
  3067.         public T getdSdalphaCj(final int j) {
  3068.             return cjCoefs[4][j];
  3069.         }

  3071.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dα.
  3072.          * <br>
  3073.          * Possible values for j are within the interval [1,jMax].
  3074.          * @param j j index
  3075.          * @return S<sup>j</sup> for the derivative dS/dα
  3076.          */
  3077.         public T getdSdalphaSj(final int j) {
  3078.             return sjCoefs[4][j];
  3079.         }

  3081.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dβ
  3082.          * <br>
  3083.          * Possible values for j are within the interval [0,jMax].
  3084.          * The value 0 is used to obtain the free coefficient C⁰
  3085.          * @param j j index
  3086.          * @return C<sup>j</sup> for the function dS/dβ
  3087.          */
  3088.         public T getdSdbetaCj(final int j) {
  3089.             return cjCoefs[5][j];
  3090.         }

  3092.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dβ.
  3093.          * <br>
  3094.          * Possible values for j are within the interval [1,jMax].
  3095.          * @param j j index
  3096.          * @return S<sup>j</sup> for the derivative dS/dβ
  3097.          */
  3098.         public T getdSdbetaSj(final int j) {
  3099.             return sjCoefs[5][j];
  3100.         }

  3102.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dγ
  3103.          * <br>
  3104.          * Possible values for j are within the interval [0,jMax].
  3105.          * The value 0 is used to obtain the free coefficient C⁰
  3106.          * @param j j index
  3107.          * @return C<sup>j</sup> for the function dS/dγ
  3108.          */
  3109.         public T getdSdgammaCj(final int j) {
  3110.             return cjCoefs[6][j];
  3111.         }

  3113.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dγ.
  3114.          * <br>
  3115.          * Possible values for j are within the interval [1,jMax].
  3116.          * @param j j index
  3117.          * @return S<sup>j</sup> for the derivative dS/dγ
  3118.          */
  3119.         public T getdSdgammaSj(final int j) {
  3120.             return sjCoefs[6][j];
  3121.         }

  3123.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dλ
  3124.          * <br>
  3125.          * Possible values for j are within the interval [0,jMax].
  3126.          * The value 0 is used to obtain the free coefficient C⁰
  3127.          * @param j j index
  3128.          * @return C<sup>j</sup> for the function dS/dλ
  3129.          */
  3130.         public T getdSdlambdaCj(final int j) {
  3131.             return cjCoefs[7][j];
  3132.         }

  3134.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dλ.
  3135.          * <br>
  3136.          * Possible values for j are within the interval [1,jMax].
  3137.          * @param j j index
  3138.          * @return S<sup>j</sup> for the derivative dS/dλ
  3139.          */
  3140.         public T getdSdlambdaSj(final int j) {
  3141.             return sjCoefs[7][j];
  3142.         }
  3143.     }

  3145.     /**
  3146.      * The coefficients used to compute the short periodic contribution for the Third body perturbation.
  3147.      * <p>
  3148.      * The short periodic contribution for the Third Body is expressed in Danielson 4.2-25.<br>
  3149.      * The coefficients C<sub>i</sub>⁰, C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup>
  3150.      * are computed by replacing the corresponding values in formula 2.5.5-10.
  3151.      * </p>
  3152.      * @author Lucian Barbulescu
  3153.      */
  3154.     private static class ThirdBodyShortPeriodicCoefficients implements ShortPeriodTerms {

  3156.         /** Maximal value for j. */
  3157.         private final int jMax;

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

  3162.         /** Max frequency of F. */
  3163.         private final int    maxFreqF;

  3165.         /** Coefficients prefix. */
  3166.         private final String prefix;

  3168.         /** All coefficients slots. */
  3169.         private final transient TimeSpanMap<Slot> slots;

  3171.         /**
  3172.          * Standard constructor.
  3173.          *  @param interpolationPoints number of points used in the interpolation process
  3174.          * @param jMax maximal value for j
  3175.          * @param maxFreqF Max frequency of F
  3176.          * @param bodyName third body name
  3177.          * @param slots all coefficients slots
  3178.          */
  3179.         ThirdBodyShortPeriodicCoefficients(final int jMax, final int interpolationPoints,
  3180.                                            final int maxFreqF, final String bodyName,
  3181.                                            final TimeSpanMap<Slot> slots) {
  3182.             this.jMax                = jMax;
  3183.             this.interpolationPoints = interpolationPoints;
  3184.             this.maxFreqF            = maxFreqF;
  3185.             this.prefix              = DSSTThirdBody.SHORT_PERIOD_PREFIX + bodyName + "-";
  3186.             this.slots               = slots;
  3187.         }

  3189.         /** Get the slot valid for some date.
  3190.          * @param meanStates mean states defining the slot
  3191.          * @return slot valid at the specified date
  3192.          */
  3193.         public Slot createSlot(final SpacecraftState... meanStates) {
  3194.             final Slot         slot  = new Slot(jMax, interpolationPoints);
  3195.             final AbsoluteDate first = meanStates[0].getDate();
  3196.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  3197.             final int compare = first.compareTo(last);
  3198.             if (compare < 0) {
  3199.                 slots.addValidAfter(slot, first, false);
  3200.             } else if (compare > 0) {
  3201.                 slots.addValidBefore(slot, first, false);
  3202.             } else {
  3203.                 // single date, valid for all time
  3204.                 slots.addValidAfter(slot, AbsoluteDate.PAST_INFINITY, false);
  3205.             }
  3206.             return slot;
  3207.         }

  3209.         /** {@inheritDoc} */
  3210.         @Override
  3211.         public double[] value(final Orbit meanOrbit) {

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

  3216.             // the current eccentric longitude
  3217.             final double F = meanOrbit.getLE();

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

  3222.             // Add the cos and sin dependent terms
  3223.             for (int j = 1; j <= maxFreqF; j++) {
  3224.                 //compute cos and sin
  3225.                 final SinCos scjF  = FastMath.sinCos(j * F);

  3227.                 final double[] c = slot.cij[j].value(meanOrbit.getDate());
  3228.                 final double[] s = slot.sij[j].value(meanOrbit.getDate());
  3229.                 for (int i = 0; i < 6; i++) {
  3230.                     shortPeriodic[i] += c[i] * scjF.cos() + s[i] * scjF.sin();
  3231.                 }
  3232.             }

  3234.             return shortPeriodic;

  3236.         }

  3238.         /** {@inheritDoc} */
  3239.         @Override
  3240.         public String getCoefficientsKeyPrefix() {
  3241.             return prefix;
  3242.         }

  3244.         /** {@inheritDoc}
  3245.          * <p>
  3246.          * For third body attraction forces,there are maxFreqF + 1 cj coefficients,
  3247.          * maxFreqF sj coefficients where maxFreqF depends on the orbit.
  3248.          * The j index is the integer multiplier for the eccentric longitude argument
  3249.          * in the cj and sj coefficients.
  3250.          * </p>
  3251.          */
  3252.         @Override
  3253.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected) {

  3255.             // select the coefficients slot
  3256.             final Slot slot = slots.get(date);

  3258.             final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * maxFreqF + 1);
  3259.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "c", 0);
  3260.             for (int j = 1; j <= maxFreqF; j++) {
  3261.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  3262.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  3263.             }
  3264.             return coefficients;

  3266.         }

  3268.         /** Put a coefficient in a map if selected.
  3269.          * @param map map to populate
  3270.          * @param selected set of coefficients that should be put in the map
  3271.          * (empty set means all coefficients are selected)
  3272.          * @param value coefficient value
  3273.          * @param id coefficient identifier
  3274.          * @param indices list of coefficient indices
  3275.          */
  3276.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
  3277.                                      final double[] value, final String id, final int... indices) {
  3278.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  3279.             keyBuilder.append(id);
  3280.             for (int index : indices) {
  3281.                 keyBuilder.append('[').append(index).append(']');
  3282.             }
  3283.             final String key = keyBuilder.toString();
  3284.             if (selected.isEmpty() || selected.contains(key)) {
  3285.                 map.put(key, value);
  3286.             }
  3287.         }

  3289.     }

  3291.     /**
  3292.      * The coefficients used to compute the short periodic contribution for the Third body perturbation.
  3293.      * <p>
  3294.      * The short periodic contribution for the Third Body is expressed in Danielson 4.2-25.<br>
  3295.      * The coefficients C<sub>i</sub>⁰, C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup>
  3296.      * are computed by replacing the corresponding values in formula 2.5.5-10.
  3297.      * </p>
  3298.      * @author Lucian Barbulescu
  3299.      */
  3300.     private static class FieldThirdBodyShortPeriodicCoefficients <T extends CalculusFieldElement<T>> implements FieldShortPeriodTerms<T> {

  3302.         /** Maximal value for j. */
  3303.         private final int jMax;

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

  3308.         /** Max frequency of F. */
  3309.         private final int    maxFreqF;

  3311.         /** Coefficients prefix. */
  3312.         private final String prefix;

  3314.         /** All coefficients slots. */
  3315.         private final transient FieldTimeSpanMap<FieldSlot<T>, T> slots;

  3317.         /**
  3318.          * Standard constructor.
  3319.          * @param interpolationPoints number of points used in the interpolation process
  3320.          * @param jMax maximal value for j
  3321.          * @param maxFreqF Max frequency of F
  3322.          * @param bodyName third body name
  3323.          * @param slots all coefficients slots
  3324.          */
  3325.         FieldThirdBodyShortPeriodicCoefficients(final int jMax, final int interpolationPoints,
  3326.                                                 final int maxFreqF, final String bodyName,
  3327.                                                 final FieldTimeSpanMap<FieldSlot<T>, T> slots) {
  3328.             this.jMax                = jMax;
  3329.             this.interpolationPoints = interpolationPoints;
  3330.             this.maxFreqF            = maxFreqF;
  3331.             this.prefix              = DSSTThirdBody.SHORT_PERIOD_PREFIX + bodyName + "-";
  3332.             this.slots               = slots;
  3333.         }

  3335.         /** Get the slot valid for some date.
  3336.          * @param meanStates mean states defining the slot
  3337.          * @return slot valid at the specified date
  3338.          */
  3339.         @SuppressWarnings("unchecked")
  3340.         public FieldSlot<T> createSlot(final FieldSpacecraftState<T>... meanStates) {
  3341.             final FieldSlot<T>         slot  = new FieldSlot<>(jMax, interpolationPoints);
  3342.             final FieldAbsoluteDate<T> first = meanStates[0].getDate();
  3343.             final FieldAbsoluteDate<T> last  = meanStates[meanStates.length - 1].getDate();
  3344.             if (first.compareTo(last) <= 0) {
  3345.                 slots.addValidAfter(slot, first);
  3346.             } else {
  3347.                 slots.addValidBefore(slot, first);
  3348.             }
  3349.             return slot;
  3350.         }

  3352.         /** {@inheritDoc} */
  3353.         @Override
  3354.         public T[] value(final FieldOrbit<T> meanOrbit) {

  3356.             // select the coefficients slot
  3357.             final FieldSlot<T> slot = slots.get(meanOrbit.getDate());

  3359.             // the current eccentric longitude
  3360.             final T F = meanOrbit.getLE();

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

  3365.             // Add the cos and sin dependent terms
  3366.             for (int j = 1; j <= maxFreqF; j++) {
  3367.                 //compute cos and sin
  3368.                 final FieldSinCos<T> scjF = FastMath.sinCos(F.multiply(j));

  3370.                 final T[] c = (T[]) slot.cij[j].value(meanOrbit.getDate());
  3371.                 final T[] s = (T[]) slot.sij[j].value(meanOrbit.getDate());
  3372.                 for (int i = 0; i < 6; i++) {
  3373.                     shortPeriodic[i] = shortPeriodic[i].add(c[i].multiply(scjF.cos()).add(s[i].multiply(scjF.sin())));
  3374.                 }
  3375.             }

  3377.             return shortPeriodic;

  3379.         }

  3381.         /** {@inheritDoc} */
  3382.         @Override
  3383.         public String getCoefficientsKeyPrefix() {
  3384.             return prefix;
  3385.         }

  3387.         /** {@inheritDoc}
  3388.          * <p>
  3389.          * For third body attraction forces,there are maxFreqF + 1 cj coefficients,
  3390.          * maxFreqF sj coefficients where maxFreqF depends on the orbit.
  3391.          * The j index is the integer multiplier for the eccentric longitude argument
  3392.          * in the cj and sj coefficients.
  3393.          * </p>
  3394.          */
  3395.         @Override
  3396.         public Map<String, T[]> getCoefficients(final FieldAbsoluteDate<T> date, final Set<String> selected) {

  3398.             // select the coefficients slot
  3399.             final FieldSlot<T> slot = slots.get(date);

  3401.             final Map<String, T[]> coefficients = new HashMap<String, T[]>(2 * maxFreqF + 1);
  3402.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "c", 0);
  3403.             for (int j = 1; j <= maxFreqF; j++) {
  3404.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  3405.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  3406.             }
  3407.             return coefficients;

  3409.         }

  3411.         /** Put a coefficient in a map if selected.
  3412.          * @param map map to populate
  3413.          * @param selected set of coefficients that should be put in the map
  3414.          * (empty set means all coefficients are selected)
  3415.          * @param value coefficient value
  3416.          * @param id coefficient identifier
  3417.          * @param indices list of coefficient indices
  3418.          */
  3419.         private void storeIfSelected(final Map<String, T[]> map, final Set<String> selected,
  3420.                                      final T[] value, final String id, final int... indices) {
  3421.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  3422.             keyBuilder.append(id);
  3423.             for (int index : indices) {
  3424.                 keyBuilder.append('[').append(index).append(']');
  3425.             }
  3426.             final String key = keyBuilder.toString();
  3427.             if (selected.isEmpty() || selected.contains(key)) {
  3428.                 map.put(key, value);
  3429.             }
  3430.         }

  3432.     }

  3434.     /** Coefficients valid for one time slot. */
  3435.     private static class Slot {

  3437.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  3438.          * <p>
  3439.          * The index order is cij[j][i] <br/>
  3440.          * i corresponds to the equinoctial element, as follows: <br/>
  3441.          * - i=0 for a <br/>
  3442.          * - i=1 for k <br/>
  3443.          * - i=2 for h <br/>
  3444.          * - i=3 for q <br/>
  3445.          * - i=4 for p <br/>
  3446.          * - i=5 for λ <br/>
  3447.          * </p>
  3448.          */
  3449.         private final ShortPeriodicsInterpolatedCoefficient[] cij;

  3451.         /** The coefficients S<sub>i</sub><sup>j</sup>.
  3452.          * <p>
  3453.          * The index order is sij[j][i] <br/>
  3454.          * i corresponds to the equinoctial element, as follows: <br/>
  3455.          * - i=0 for a <br/>
  3456.          * - i=1 for k <br/>
  3457.          * - i=2 for h <br/>
  3458.          * - i=3 for q <br/>
  3459.          * - i=4 for p <br/>
  3460.          * - i=5 for λ <br/>
  3461.          * </p>
  3462.          */
  3463.         private final ShortPeriodicsInterpolatedCoefficient[] sij;

  3465.         /** Simple constructor.
  3466.          *  @param jMax maximum value for j index
  3467.          *  @param interpolationPoints number of points used in the interpolation process
  3468.          */
  3469.         Slot(final int jMax, final int interpolationPoints) {
  3470.             // allocate the coefficients arrays
  3471.             cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
  3472.             sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
  3473.             for (int j = 0; j <= jMax; j++) {
  3474.                 cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3475.                 sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3476.             }


  3479.         }
  3480.     }

  3482.     /** Coefficients valid for one time slot. */
  3483.     private static class FieldSlot <T extends CalculusFieldElement<T>> {

  3485.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  3486.          * <p>
  3487.          * The index order is cij[j][i] <br/>
  3488.          * i corresponds to the equinoctial element, as follows: <br/>
  3489.          * - i=0 for a <br/>
  3490.          * - i=1 for k <br/>
  3491.          * - i=2 for h <br/>
  3492.          * - i=3 for q <br/>
  3493.          * - i=4 for p <br/>
  3494.          * - i=5 for λ <br/>
  3495.          * </p>
  3496.          */
  3497.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] cij;

  3499.         /** The coefficients S<sub>i</sub><sup>j</sup>.
  3500.          * <p>
  3501.          * The index order is sij[j][i] <br/>
  3502.          * i corresponds to the equinoctial element, as follows: <br/>
  3503.          * - i=0 for a <br/>
  3504.          * - i=1 for k <br/>
  3505.          * - i=2 for h <br/>
  3506.          * - i=3 for q <br/>
  3507.          * - i=4 for p <br/>
  3508.          * - i=5 for λ <br/>
  3509.          * </p>
  3510.          */
  3511.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] sij;

  3513.         /** Simple constructor.
  3514.          *  @param jMax maximum value for j index
  3515.          *  @param interpolationPoints number of points used in the interpolation process
  3516.          */
  3517.         @SuppressWarnings("unchecked")
  3518.         FieldSlot(final int jMax, final int interpolationPoints) {
  3519.             // allocate the coefficients arrays
  3520.             cij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);
  3521.             sij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);
  3522.             for (int j = 0; j <= jMax; j++) {
  3523.                 cij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3524.                 sij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3525.             }


  3528.         }
  3529.     }

  3531.     /** Compute potential and potential derivatives with respect to orbital parameters. */
  3532.     private class UAnddU {

  3534.         /** The current value of the U function. <br/>
  3535.          * Needed for the short periodic contribution */
  3536.         private double U;

  3538.         /** dU / da. */
  3539.         private  double dUda;

  3541.         /** dU / dk. */
  3542.         private double dUdk;

  3544.         /** dU / dh. */
  3545.         private double dUdh;

  3547.         /** dU / dAlpha. */
  3548.         private double dUdAl;

  3550.         /** dU / dBeta. */
  3551.         private double dUdBe;

  3553.         /** dU / dGamma. */
  3554.         private double dUdGa;

  3556.         /** Simple constuctor.
  3557.          * @param context container for attributes
  3558.          * @param hansen hansen objects
  3559.          * @param aoR3Pow a / R3 up to power maxAR3Pow
  3560.          * @param qns Qns coefficients
  3561.          */
  3562.         UAnddU(final DSSTThirdBodyDynamicContext context, final HansenObjects hansen,
  3563.                final double[] aoR3Pow, final double[][] qns) {
  3564.             // Auxiliary elements related to the current orbit
  3565.             final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();

  3567.             // Gs and Hs coefficients
  3568.             final double[][] GsHs = CoefficientsFactory.computeGsHs(auxiliaryElements.getK(), auxiliaryElements.getH(), context.getAlpha(), context.getBeta(), staticContext.getMaxEccPow());

  3570.             // Initialise U.
  3571.             U = 0.;

  3573.             // Potential derivatives
  3574.             dUda  = 0.;
  3575.             dUdk  = 0.;
  3576.             dUdh  = 0.;
  3577.             dUdAl = 0.;
  3578.             dUdBe = 0.;
  3579.             dUdGa = 0.;

  3581.             for (int s = 0; s <= staticContext.getMaxEccPow(); s++) {

  3583.                 // initialise the Hansen roots
  3584.                 hansen.computeHansenObjectsInitValues(context, auxiliaryElements.getB(), s);

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

  3589.                 // Compute Gs partial derivatives from 3.1-(9)
  3590.                 double dGsdh  = 0.;
  3591.                 double dGsdk  = 0.;
  3592.                 double dGsdAl = 0.;
  3593.                 double dGsdBe = 0.;
  3594.                 if (s > 0) {
  3595.                     // First get the G(s-1) and the H(s-1) coefficients
  3596.                     final double sxGsm1 = s * GsHs[0][s - 1];
  3597.                     final double sxHsm1 = s * GsHs[1][s - 1];
  3598.                     // Then compute derivatives
  3599.                     dGsdh  = context.getBeta()  * sxGsm1 - context.getAlpha() * sxHsm1;
  3600.                     dGsdk  = context.getAlpha() * sxGsm1 + context.getBeta()  * sxHsm1;
  3601.                     dGsdAl = auxiliaryElements.getK() * sxGsm1 - auxiliaryElements.getH() * sxHsm1;
  3602.                     dGsdBe = auxiliaryElements.getH() * sxGsm1 + auxiliaryElements.getK() * sxHsm1;
  3603.                 }

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

  3608.                 for (int n = FastMath.max(2, s); n <= staticContext.getMaxAR3Pow(); n++) {
  3609.                     // (n - s) must be even
  3610.                     if ((n - s) % 2 == 0) {
  3611.                         // Extract data from previous computation :
  3612.                         final double kns   = hansen.getHansenObjects()[s].getValue(n, auxiliaryElements.getB());
  3613.                         final double dkns  = hansen.getHansenObjects()[s].getDerivative(n, auxiliaryElements.getB());

  3615.                         final double vns   = Vns.get(new NSKey(n, s));
  3616.                         final double coef0 = delta0s * aoR3Pow[n] * vns;
  3617.                         final double coef1 = coef0 * qns[n][s];
  3618.                         final double coef2 = coef1 * kns;
  3619.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  3620.                         // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
  3621.                         final double dqns = (n == s) ? 0. : qns[n][s + 1];

  3623.                         //Compute U:
  3624.                         U += coef2 * gs;

  3626.                         // Compute dU / da :
  3627.                         dUda  += coef2 * n * gs;
  3628.                         // Compute dU / dh
  3629.                         dUdh  += coef1 * (kns * dGsdh + context.getHXXX() * gs * dkns);
  3630.                         // Compute dU / dk
  3631.                         dUdk  += coef1 * (kns * dGsdk + context.getKXXX() * gs * dkns);
  3632.                         // Compute dU / dAlpha
  3633.                         dUdAl += coef2 * dGsdAl;
  3634.                         // Compute dU / dBeta
  3635.                         dUdBe += coef2 * dGsdBe;
  3636.                         // Compute dU / dGamma
  3637.                         dUdGa += coef0 * kns * dqns * gs;
  3638.                     }
  3639.                 }
  3640.             }

  3642.             // multiply by mu3 / R3
  3643.             this.U = U * context.getMuoR3();

  3645.             this.dUda  = dUda  * context.getMuoR3() / auxiliaryElements.getSma();
  3646.             this.dUdk  = dUdk  * context.getMuoR3();
  3647.             this.dUdh  = dUdh  * context.getMuoR3();
  3648.             this.dUdAl = dUdAl * context.getMuoR3();
  3649.             this.dUdBe = dUdBe * context.getMuoR3();
  3650.             this.dUdGa = dUdGa * context.getMuoR3();

  3652.         }

  3654.         /** Return value of U.
  3655.          * @return U
  3656.          */
  3657.         public double getU() {
  3658.             return U;
  3659.         }

  3661.         /** Return value of dU / da.
  3662.          * @return dUda
  3663.          */
  3664.         public double getdUda() {
  3665.             return dUda;
  3666.         }

  3668.         /** Return value of dU / dk.
  3669.          * @return dUdk
  3670.          */
  3671.         public double getdUdk() {
  3672.             return dUdk;
  3673.         }

  3675.         /** Return value of dU / dh.
  3676.          * @return dUdh
  3677.          */
  3678.         public double getdUdh() {
  3679.             return dUdh;
  3680.         }

  3682.         /** Return value of dU / dAlpha.
  3683.          * @return dUdAl
  3684.          */
  3685.         public double getdUdAl() {
  3686.             return dUdAl;
  3687.         }

  3689.         /** Return value of dU / dBeta.
  3690.          * @return dUdBe
  3691.          */
  3692.         public double getdUdBe() {
  3693.             return dUdBe;
  3694.         }

  3696.         /** Return value of dU / dGamma.
  3697.          * @return dUdGa
  3698.          */
  3699.         public double getdUdGa() {
  3700.             return dUdGa;
  3701.         }

  3703.     }

  3705.     /** Compute potential and potential derivatives with respect to orbital parameters. */
  3706.     private class FieldUAnddU <T extends CalculusFieldElement<T>> {

  3708.         /** The current value of the U function. <br/>
  3709.          * Needed for the short periodic contribution */
  3710.         private T U;

  3712.         /** dU / da. */
  3713.         private T dUda;

  3715.         /** dU / dk. */
  3716.         private T dUdk;

  3718.         /** dU / dh. */
  3719.         private T dUdh;

  3721.         /** dU / dAlpha. */
  3722.         private T dUdAl;

  3724.         /** dU / dBeta. */
  3725.         private T dUdBe;

  3727.         /** dU / dGamma. */
  3728.         private T dUdGa;

  3730.         /** Simple constuctor.
  3731.          * @param context container for attributes
  3732.          * @param hansen hansen objects
  3733.          * @param aoR3Pow a / R3 up to power maxAR3Pow
  3734.          * @param qns Qns coefficients
  3735.          */
  3736.         FieldUAnddU(final FieldDSSTThirdBodyDynamicContext<T> context, final FieldHansenObjects<T> hansen,
  3737.                     final T[] aoR3Pow, final T[][] qns) {

  3739.             // Auxiliary elements related to the current orbit
  3740.             final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();

  3742.             // Field for array building
  3743.             final Field<T> field = auxiliaryElements.getDate().getField();
  3744.             // Zero for initialization
  3745.             final T zero = field.getZero();

  3747.             // Gs and Hs coefficients
  3748.             final T[][] GsHs = CoefficientsFactory.computeGsHs(auxiliaryElements.getK(), auxiliaryElements.getH(), context.getAlpha(), context.getBeta(), staticContext.getMaxEccPow(), field);

  3750.             // Initialise U.
  3751.             U = zero;

  3753.             // Potential derivatives
  3754.             dUda  = zero;
  3755.             dUdk  = zero;
  3756.             dUdh  = zero;
  3757.             dUdAl = zero;
  3758.             dUdBe = zero;
  3759.             dUdGa = zero;

  3761.             for (int s = 0; s <= staticContext.getMaxEccPow(); s++) {
  3762.                 // initialise the Hansen roots
  3763.                 hansen.computeHansenObjectsInitValues(context, auxiliaryElements.getB(), s);

  3765.                 // Get the current Gs coefficient
  3766.                 final T gs = GsHs[0][s];

  3768.                 // Compute Gs partial derivatives from 3.1-(9)
  3769.                 T dGsdh  = zero;
  3770.                 T dGsdk  = zero;
  3771.                 T dGsdAl = zero;
  3772.                 T dGsdBe = zero;
  3773.                 if (s > 0) {
  3774.                     // First get the G(s-1) and the H(s-1) coefficients
  3775.                     final T sxGsm1 = GsHs[0][s - 1].multiply(s);
  3776.                     final T sxHsm1 = GsHs[1][s - 1].multiply(s);
  3777.                     // Then compute derivatives
  3778.                     dGsdh  = sxGsm1.multiply(context.getBeta()).subtract(sxHsm1.multiply(context.getAlpha()));
  3779.                     dGsdk  = sxGsm1.multiply(context.getAlpha()).add(sxHsm1.multiply(context.getBeta()));
  3780.                     dGsdAl = sxGsm1.multiply(auxiliaryElements.getK()).subtract(sxHsm1.multiply(auxiliaryElements.getH()));
  3781.                     dGsdBe = sxGsm1.multiply(auxiliaryElements.getH()).add(sxHsm1.multiply(auxiliaryElements.getK()));
  3782.                 }

  3784.                 // Kronecker symbol (2 - delta(0,s))
  3785.                 final T delta0s = zero.add((s == 0) ? 1. : 2.);

  3787.                 for (int n = FastMath.max(2, s); n <= staticContext.getMaxAR3Pow(); n++) {
  3788.                     // (n - s) must be even
  3789.                     if ((n - s) % 2 == 0) {
  3790.                         // Extract data from previous computation :
  3791.                         final T kns   = (T) hansen.getHansenObjects()[s].getValue(n, auxiliaryElements.getB());
  3792.                         final T dkns  = (T) hansen.getHansenObjects()[s].getDerivative(n, auxiliaryElements.getB());

  3794.                         final double vns = Vns.get(new NSKey(n, s));
  3795.                         final T coef0 = delta0s.multiply(vns).multiply(aoR3Pow[n]);
  3796.                         final T coef1 = coef0.multiply(qns[n][s]);
  3797.                         final T coef2 = coef1.multiply(kns);
  3798.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  3799.                         // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
  3800.                         final T dqns = (n == s) ? zero : qns[n][s + 1];

  3802.                         //Compute U:
  3803.                         U = U.add(coef2.multiply(gs));

  3805.                         // Compute dU / da :
  3806.                         dUda  = dUda.add(coef2.multiply(n).multiply(gs));
  3807.                         // Compute dU / dh
  3808.                         dUdh  = dUdh.add(coef1.multiply(dGsdh.multiply(kns).add(context.getHXXX().multiply(gs).multiply(dkns))));
  3809.                         // Compute dU / dk
  3810.                         dUdk  = dUdk.add(coef1.multiply(dGsdk.multiply(kns).add(context.getKXXX().multiply(gs).multiply(dkns))));
  3811.                         // Compute dU / dAlpha
  3812.                         dUdAl = dUdAl.add(coef2.multiply(dGsdAl));
  3813.                         // Compute dU / dBeta
  3814.                         dUdBe = dUdBe.add(coef2.multiply(dGsdBe));
  3815.                         // Compute dU / dGamma
  3816.                         dUdGa = dUdGa.add(coef0.multiply(kns).multiply(dqns).multiply(gs));
  3817.                     }
  3818.                 }
  3819.             }

  3821.             // multiply by mu3 / R3
  3822.             this.U = U.multiply(context.getMuoR3());

  3824.             this.dUda  = dUda.multiply(context.getMuoR3().divide(auxiliaryElements.getSma()));
  3825.             this.dUdk  = dUdk.multiply(context.getMuoR3());
  3826.             this.dUdh  = dUdh.multiply(context.getMuoR3());
  3827.             this.dUdAl = dUdAl.multiply(context.getMuoR3());
  3828.             this.dUdBe = dUdBe.multiply(context.getMuoR3());
  3829.             this.dUdGa = dUdGa.multiply(context.getMuoR3());

  3831.         }

  3833.         /** Return value of U.
  3834.          * @return U
  3835.          */
  3836.         public T getU() {
  3837.             return U;
  3838.         }

  3840.         /** Return value of dU / da.
  3841.          * @return dUda
  3842.          */
  3843.         public T getdUda() {
  3844.             return dUda;
  3845.         }

  3847.         /** Return value of dU / dk.
  3848.          * @return dUdk
  3849.          */
  3850.         public T getdUdk() {
  3851.             return dUdk;
  3852.         }

  3854.         /** Return value of dU / dh.
  3855.          * @return dUdh
  3856.          */
  3857.         public T getdUdh() {
  3858.             return dUdh;
  3859.         }

  3861.         /** Return value of dU / dAlpha.
  3862.          * @return dUdAl
  3863.          */
  3864.         public T getdUdAl() {
  3865.             return dUdAl;
  3866.         }

  3868.         /** Return value of dU / dBeta.
  3869.          * @return dUdBe
  3870.          */
  3871.         public T getdUdBe() {
  3872.             return dUdBe;
  3873.         }

  3875.         /** Return value of dU / dGamma.
  3876.          * @return dUdGa
  3877.          */
  3878.         public T getdUdGa() {
  3879.             return dUdGa;
  3880.         }

  3882.     }

  3884.     /** Computes init values of the Hansen Objects. */
  3885.     private static class HansenObjects {

  3887.         /** Max power for summation. */
  3888.         private static final int    MAX_POWER = 22;

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

  3894.         /** Simple constructor. */
  3895.         HansenObjects() {
  3896.             this.hansenObjects = new HansenThirdBodyLinear[MAX_POWER + 1];
  3897.             for (int s = 0; s <= MAX_POWER; s++) {
  3898.                 this.hansenObjects[s] = new HansenThirdBodyLinear(MAX_POWER, s);
  3899.             }
  3900.         }

  3902.         /** Compute init values for hansen objects.
  3903.          * @param context container for attributes
  3904.          * @param B = sqrt(1 - e²).
  3905.          * @param element element of the array to compute the init values
  3906.          */
  3907.         public void computeHansenObjectsInitValues(final DSSTThirdBodyDynamicContext context, final double B, final int element) {
  3908.             hansenObjects[element].computeInitValues(B, context.getBB(), context.getBBB());
  3909.         }

  3911.         /** Get the Hansen Objects.
  3912.          * @return hansenObjects
  3913.          */
  3914.         public HansenThirdBodyLinear[] getHansenObjects() {
  3915.             return hansenObjects;
  3916.         }

  3918.     }

  3920.     /** Computes init values of the Hansen Objects. */
  3921.     private static class FieldHansenObjects<T extends CalculusFieldElement<T>> {

  3923.         /** Max power for summation. */
  3924.         private static final int    MAX_POWER = 22;

  3926.         /** An array that contains the objects needed to build the Hansen coefficients. <br/>
  3927.          * The index is s */
  3928.         private final FieldHansenThirdBodyLinear<T>[] hansenObjects;

  3930.         /** Simple constructor.
  3931.          * @param field field used by default
  3932.          */
  3933.         @SuppressWarnings("unchecked")
  3934.         FieldHansenObjects(final Field<T> field) {
  3935.             this.hansenObjects = (FieldHansenThirdBodyLinear<T>[]) Array.newInstance(FieldHansenThirdBodyLinear.class, MAX_POWER + 1);
  3936.             for (int s = 0; s <= MAX_POWER; s++) {
  3937.                 this.hansenObjects[s] = new FieldHansenThirdBodyLinear<>(MAX_POWER, s, field);
  3938.             }
  3939.         }

  3941.         /** Initialise the Hansen roots for third body problem.
  3942.          * @param context container for attributes
  3943.          * @param B = sqrt(1 - e²).
  3944.          * @param element element of the array to compute the init values
  3945.          */
  3946.         public void computeHansenObjectsInitValues(final FieldDSSTThirdBodyDynamicContext<T> context,
  3947.                                                    final T B, final int element) {
  3948.             hansenObjects[element].computeInitValues(B, context.getBB(), context.getBBB());
  3949.         }

  3951.         /** Get the Hansen Objects.
  3952.          * @return hansenObjects
  3953.          */
  3954.         public FieldHansenThirdBodyLinear<T>[] getHansenObjects() {
  3955.             return hansenObjects;
  3956.         }

  3958.     }

  3960. }
Created with JaCoCo 0.8.11.202310140853