DSSTThirdBody.java

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

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

  65. /** Third body attraction perturbation to the
  66.  *  {@link org.orekit.propagation.semianalytical.dsst.DSSTPropagator DSSTPropagator}.
  67.  *
  68.  *  @author Romain Di Costanzo
  69.  *  @author Pascal Parraud
  70.  *  @author Bryan Cazabonne (field translation)
  71.  */
  72. public class DSSTThirdBody implements DSSTForceModel {

  74.     /**  Name of the prefix for short period coefficients keys. */
  75.     public static final String SHORT_PERIOD_PREFIX = "DSST-3rd-body-";

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

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

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

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

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

  109.     /** Max power for summation. */
  110.     private static final int    MAX_POWER = 22;

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

  115.     /** Max frequency of F. */
  116.     private int    maxFreqF;

  118.     /** Max frequency of F for field initialization. */
  119.     private int    maxFieldFreqF;

  121.     /** The 3rd body to consider. */
  122.     private final CelestialBody    body;

  124.     /** Short period terms. */
  125.     private ThirdBodyShortPeriodicCoefficients shortPeriods;

  127.     /** Short period terms. */
  128.     private Map<Field<?>, FieldThirdBodyShortPeriodicCoefficients<?>> fieldShortPeriods;

  130.     /** Drivers for third body attraction coefficient. */
  131.     private final ParameterDriver bodyParameterDriver;

  133.     /** Driver for gravitational parameter. */
  134.     private final ParameterDriver gmParameterDriver;

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

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

  142.     /** Flag for force model initialization with field elements. */
  143.     private boolean pendingInitialization;

  145.     /** Complete constructor.
  146.      *  @param body the 3rd body to consider
  147.      *  @param mu central attraction coefficient
  148.      *  @see CelestialBodies
  149.      */
  150.     public DSSTThirdBody(final CelestialBody body, final double mu) {
  151.         bodyParameterDriver = new ParameterDriver(body.getName() + DSSTThirdBody.ATTRACTION_COEFFICIENT,
  152.                                                   body.getGM(), MU_SCALE,
  153.                                                   0.0, Double.POSITIVE_INFINITY);
  154.         gmParameterDriver = new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
  155.                                                 mu, MU_SCALE,
  156.                                                 0.0, Double.POSITIVE_INFINITY);

  158.         this.body = body;
  159.         this.Vns  = CoefficientsFactory.computeVns(MAX_POWER);

  161.         pendingInitialization = true;

  163.         fieldShortPeriods = new HashMap<>();
  164.         fieldHansen       = new HashMap<>();
  165.     }

  167.     /** Get third body.
  168.      *  @return third body
  169.      */
  170.     public CelestialBody getBody() {
  171.         return body;
  172.     }

  174.     /** Computes the highest power of the eccentricity and the highest power
  175.      *  of a/R3 to appear in the truncated analytical power series expansion.
  176.      *  <p>
  177.      *  This method computes the upper value for the 3rd body potential and
  178.      *  determines the maximal powers for the eccentricity and a/R3 producing
  179.      *  potential terms bigger than a defined tolerance.
  180.      *  </p>
  181.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  182.      *  @param type type of the elements used during the propagation
  183.      *  @param parameters values of the force model parameters
  184.      */
  185.     @Override
  186.     public List<ShortPeriodTerms> initialize(final AuxiliaryElements auxiliaryElements,
  187.                                              final PropagationType type,
  188.                                              final double[] parameters) {

  190.         // Initializes specific parameters.
  191.         final DSSTThirdBodyContext context = initializeStep(auxiliaryElements, parameters);

  193.         maxFreqF = context.getMaxFreqF();

  195.         hansen = new HansenObjects();

  197.         final int jMax = maxFreqF;
  198.         shortPeriods = new ThirdBodyShortPeriodicCoefficients(jMax, INTERPOLATION_POINTS,
  199.                                                               maxFreqF, body.getName(),
  200.                                                               new TimeSpanMap<Slot>(new Slot(jMax, INTERPOLATION_POINTS)));

  202.         final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
  203.         list.add(shortPeriods);
  204.         return list;

  206.     }

  208.     /** {@inheritDoc} */
  209.     @Override
  210.     public <T extends RealFieldElement<T>> List<FieldShortPeriodTerms<T>> initialize(final FieldAuxiliaryElements<T> auxiliaryElements,
  211.                                                                                      final PropagationType type,
  212.                                                                                      final T[] parameters) {

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

  217.         if (pendingInitialization == true) {
  218.             // Initializes specific parameters.
  219.             final FieldDSSTThirdBodyContext<T> context = initializeStep(auxiliaryElements, parameters);

  221.             maxFieldFreqF = context.getMaxFreqF();

  223.             fieldHansen.put(field, new FieldHansenObjects<>(field));

  225.             pendingInitialization = false;
  226.         }

  228.         final int jMax = maxFieldFreqF;
  229.         final FieldThirdBodyShortPeriodicCoefficients<T> ftbspc =
  230.                         new FieldThirdBodyShortPeriodicCoefficients<>(jMax, INTERPOLATION_POINTS,
  231.                                                                       maxFieldFreqF, body.getName(),
  232.                                                                       new FieldTimeSpanMap<>(new FieldSlot<>(jMax,
  233.                                                                                                              INTERPOLATION_POINTS),
  234.                                                                                              field));
  235.         fieldShortPeriods.put(field, ftbspc);
  236.         return Collections.singletonList(ftbspc);
  237.     }

  239.     /** Performs initialization at each integration step for the current force model.
  240.      *  <p>
  241.      *  This method aims at being called before mean elements rates computation.
  242.      *  </p>
  243.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  244.      *  @param parameters values of the force model parameters
  245.      *  @return new force model context
  246.      */
  247.     private DSSTThirdBodyContext initializeStep(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
  248.         return new DSSTThirdBodyContext(auxiliaryElements, body, parameters);
  249.     }

  251.     /** Performs initialization at each integration step for the current force model.
  252.      *  <p>
  253.      *  This method aims at being called before mean elements rates computation.
  254.      *  </p>
  255.      *  @param <T> type of the elements
  256.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  257.      *  @param parameters values of the force model parameters
  258.      *  @return new force model context
  259.      */
  260.     private <T extends RealFieldElement<T>> FieldDSSTThirdBodyContext<T> initializeStep(final FieldAuxiliaryElements<T> auxiliaryElements,
  261.                                                                                         final T[] parameters) {
  262.         return new FieldDSSTThirdBodyContext<>(auxiliaryElements, body, parameters);
  263.     }

  265.     /** {@inheritDoc} */
  266.     @Override
  267.     public double[] getMeanElementRate(final SpacecraftState currentState,
  268.                                        final AuxiliaryElements auxiliaryElements, final double[] parameters) {

  270.         // Container for attributes
  271.         final DSSTThirdBodyContext context = initializeStep(auxiliaryElements, parameters);
  272.         // Access to potential U derivatives
  273.         final UAnddU udu = new UAnddU(context, hansen);

  275.         // Compute cross derivatives [Eq. 2.2-(8)]
  276.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  277.         final double UAlphaGamma   = context.getAlpha() * udu.getdUdGa() - context.getGamma() * udu.getdUdAl();
  278.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  279.         final double UBetaGamma    =  context.getBeta() * udu.getdUdGa() - context.getGamma() * udu.getdUdBe();
  280.         // Common factor
  281.         final double pUAGmIqUBGoAB = (auxiliaryElements.getP() * UAlphaGamma - I * auxiliaryElements.getQ() * UBetaGamma) * context.getOoAB();

  283.         // Compute mean elements rates [Eq. 3.1-(1)]
  284.         final double da =  0.;
  285.         final double dh =  context.getBoA() * udu.getdUdk() + auxiliaryElements.getK() * pUAGmIqUBGoAB;
  286.         final double dk = -context.getBoA() * udu.getdUdh() - auxiliaryElements.getH() * pUAGmIqUBGoAB;
  287.         final double dp =  context.getMCo2AB() * UBetaGamma;
  288.         final double dq =  context.getMCo2AB() * UAlphaGamma * I;
  289.         final double dM =  context.getM2aoA() * udu.getdUda() + context.getBoABpo() * (auxiliaryElements.getH() * udu.getdUdh() + auxiliaryElements.getK() * udu.getdUdk()) + pUAGmIqUBGoAB;

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

  293.     }

  295.     /** {@inheritDoc} */
  296.     @Override
  297.     public <T extends RealFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> currentState,
  298.                                                                   final FieldAuxiliaryElements<T> auxiliaryElements,
  299.                                                                   final T[] parameters) {

  301.         // Parameters for array building
  302.         final Field<T> field = currentState.getDate().getField();
  303.         final T        zero  = field.getZero();

  305.         // Container for attributes
  306.         final FieldDSSTThirdBodyContext<T> context = initializeStep(auxiliaryElements, parameters);

  308.         @SuppressWarnings("unchecked")
  309.         final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);

  311.         // Access to potential U derivatives
  312.         final FieldUAnddU<T> udu = new FieldUAnddU<>(context, fho);

  314.         // Compute cross derivatives [Eq. 2.2-(8)]
  315.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  316.         final T UAlphaGamma   = udu.getdUdGa().multiply(context.getAlpha()).subtract(udu.getdUdAl().multiply(context.getGamma()));
  317.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  318.         final T UBetaGamma    = udu.getdUdGa().multiply(context.getBeta()).subtract(udu.getdUdBe().multiply(context.getGamma()));
  319.         // Common factor
  320.         final T pUAGmIqUBGoAB = (UAlphaGamma.multiply(auxiliaryElements.getP()).subtract(UBetaGamma.multiply(auxiliaryElements.getQ()).multiply(I))).multiply(context.getOoAB());

  322.         // Compute mean elements rates [Eq. 3.1-(1)]
  323.         final T da =  zero;
  324.         final T dh =  udu.getdUdk().multiply(context.getBoA()).add(pUAGmIqUBGoAB.multiply(auxiliaryElements.getK()));
  325.         final T dk =  ((udu.getdUdh().multiply(context.getBoA())).negate()).subtract(pUAGmIqUBGoAB.multiply(auxiliaryElements.getH()));
  326.         final T dp =  UBetaGamma.multiply(context.getMCo2AB());
  327.         final T dq =  UAlphaGamma.multiply(I).multiply(context.getMCo2AB());
  328.         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()));

  330.         final T[] elements = MathArrays.buildArray(field, 6);
  331.         elements[0] = da;
  332.         elements[1] = dk;
  333.         elements[2] = dh;
  334.         elements[3] = dq;
  335.         elements[4] = dp;
  336.         elements[5] = dM;

  338.         return elements;

  340.     }

  342.     /** {@inheritDoc} */
  343.     @Override
  344.     public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {

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

  348.         for (final SpacecraftState meanState : meanStates) {

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

  353.             // Container of attributes
  354.             final DSSTThirdBodyContext context = initializeStep(auxiliaryElements, parameters);

  356.             final GeneratingFunctionCoefficients gfCoefs =
  357.                             new GeneratingFunctionCoefficients(context.getMaxAR3Pow(), MAX_ECCPOWER_SP, context.getMaxAR3Pow() + 1, context, hansen);

  359.             //Compute additional quantities
  360.             // 2 * a / An
  361.             final double ax2oAn  = -context.getM2aoA() / context.getMeanMotion();
  362.             // B / An
  363.             final double BoAn    = context.getBoA() / context.getMeanMotion();
  364.             // 1 / ABn
  365.             final double ooABn   = context.getOoAB() / context.getMeanMotion();
  366.             // C / 2ABn
  367.             final double Co2ABn  = -context.getMCo2AB() / context.getMeanMotion();
  368.             // B / (A * (1 + B) * n)
  369.             final double BoABpon = context.getBoABpo() / context.getMeanMotion();
  370.             // -3 / n²a² = -3 / nA
  371.             final double m3onA   = -3 / (context.getA() * context.getMeanMotion());

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

  378.                 // Compute the cross derivatives operator :
  379.                 final double SAlphaGammaCj    = context.getAlpha() * gfCoefs.getdSdgammaCj(j) - context.getGamma() * gfCoefs.getdSdalphaCj(j);
  380.                 final double SAlphaBetaCj     = context.getAlpha() * gfCoefs.getdSdbetaCj(j)  - context.getBeta()  * gfCoefs.getdSdalphaCj(j);
  381.                 final double SBetaGammaCj     = context.getBeta() * gfCoefs.getdSdgammaCj(j) - context.getGamma() * gfCoefs.getdSdbetaCj(j);
  382.                 final double ShkCj            = auxiliaryElements.getH() * gfCoefs.getdSdkCj(j)     -  auxiliaryElements.getK()    * gfCoefs.getdSdhCj(j);
  383.                 final double pSagmIqSbgoABnCj = (auxiliaryElements.getP() * SAlphaGammaCj - I * auxiliaryElements.getQ() * SBetaGammaCj) * ooABn;
  384.                 final double ShkmSabmdSdlCj   = ShkCj - SAlphaBetaCj - gfCoefs.getdSdlambdaCj(j);

  386.                 currentCij[0] =  ax2oAn * gfCoefs.getdSdlambdaCj(j);
  387.                 currentCij[1] =  -(BoAn * gfCoefs.getdSdhCj(j) + auxiliaryElements.getH() * pSagmIqSbgoABnCj + auxiliaryElements.getK() * BoABpon * gfCoefs.getdSdlambdaCj(j));
  388.                 currentCij[2] =    BoAn * gfCoefs.getdSdkCj(j) + auxiliaryElements.getK() * pSagmIqSbgoABnCj - auxiliaryElements.getH() * BoABpon * gfCoefs.getdSdlambdaCj(j);
  389.                 currentCij[3] =  Co2ABn * (auxiliaryElements.getQ() * ShkmSabmdSdlCj - I * SAlphaGammaCj);
  390.                 currentCij[4] =  Co2ABn * (auxiliaryElements.getP() * ShkmSabmdSdlCj - SBetaGammaCj);
  391.                 currentCij[5] = -ax2oAn * gfCoefs.getdSdaCj(j) + BoABpon * (auxiliaryElements.getH() * gfCoefs.getdSdhCj(j) + auxiliaryElements.getK() * gfCoefs.getdSdkCj(j)) + pSagmIqSbgoABnCj + m3onA * gfCoefs.getSCj(j);

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

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

  399.                 // Compute the cross derivatives operator :
  400.                 final double SAlphaGammaSj    = context.getAlpha() * gfCoefs.getdSdgammaSj(j) - context.getGamma() * gfCoefs.getdSdalphaSj(j);
  401.                 final double SAlphaBetaSj     = context.getAlpha() * gfCoefs.getdSdbetaSj(j)  - context.getBeta()  * gfCoefs.getdSdalphaSj(j);
  402.                 final double SBetaGammaSj     =  context.getBeta() * gfCoefs.getdSdgammaSj(j) - context.getGamma() * gfCoefs.getdSdbetaSj(j);
  403.                 final double ShkSj            =     auxiliaryElements.getH() * gfCoefs.getdSdkSj(j)     -  auxiliaryElements.getK()    * gfCoefs.getdSdhSj(j);
  404.                 final double pSagmIqSbgoABnSj = (auxiliaryElements.getP() * SAlphaGammaSj - I * auxiliaryElements.getQ() * SBetaGammaSj) * ooABn;
  405.                 final double ShkmSabmdSdlSj   =  ShkSj - SAlphaBetaSj - gfCoefs.getdSdlambdaSj(j);

  407.                 currentSij[0] =  ax2oAn * gfCoefs.getdSdlambdaSj(j);
  408.                 currentSij[1] =  -(BoAn * gfCoefs.getdSdhSj(j) + auxiliaryElements.getH() * pSagmIqSbgoABnSj + auxiliaryElements.getK() * BoABpon * gfCoefs.getdSdlambdaSj(j));
  409.                 currentSij[2] =    BoAn * gfCoefs.getdSdkSj(j) + auxiliaryElements.getK() * pSagmIqSbgoABnSj - auxiliaryElements.getH() * BoABpon * gfCoefs.getdSdlambdaSj(j);
  410.                 currentSij[3] =  Co2ABn * (auxiliaryElements.getQ() * ShkmSabmdSdlSj - I * SAlphaGammaSj);
  411.                 currentSij[4] =  Co2ABn * (auxiliaryElements.getP() * ShkmSabmdSdlSj - SBetaGammaSj);
  412.                 currentSij[5] = -ax2oAn * gfCoefs.getdSdaSj(j) + BoABpon * (auxiliaryElements.getH() * gfCoefs.getdSdhSj(j) + auxiliaryElements.getK() * gfCoefs.getdSdkSj(j)) + pSagmIqSbgoABnSj + m3onA * gfCoefs.getSSj(j);

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

  417.                 if (j == 1) {
  418.                     //Compute the C⁰ coefficients using Danielson 2.5.2-15a.
  419.                     final double[] value = new double[6];
  420.                     for (int i = 0; i < 6; ++i) {
  421.                         value[i] = currentCij[i] * auxiliaryElements.getK() / 2. + currentSij[i] * auxiliaryElements.getH() / 2.;
  422.                     }
  423.                     slot.cij[0].addGridPoint(meanState.getDate(), value);
  424.                 }
  425.             }
  426.         }
  427.     }

  429.     /** {@inheritDoc} */
  430.     @Override
  431.     @SuppressWarnings("unchecked")
  432.     public <T extends RealFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
  433.                                                                        final FieldSpacecraftState<T>... meanStates) {

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

  438.         final FieldThirdBodyShortPeriodicCoefficients<T> ftbspc = (FieldThirdBodyShortPeriodicCoefficients<T>) fieldShortPeriods.get(field);
  439.         final FieldSlot<T> slot = ftbspc.createSlot(meanStates);
  440.         for (final FieldSpacecraftState<T> meanState : meanStates) {

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

  445.             // Container of attributes
  446.             final FieldDSSTThirdBodyContext<T> context = initializeStep(auxiliaryElements, parameters);

  448.             final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);

  450.             final FieldGeneratingFunctionCoefficients<T> gfCoefs =
  451.                             new FieldGeneratingFunctionCoefficients<>(context.getMaxAR3Pow(), MAX_ECCPOWER_SP, context.getMaxAR3Pow() + 1, context, fho, field);

  453.             //Compute additional quantities
  454.             // 2 * a / An
  455.             final T ax2oAn  = context.getM2aoA().negate().divide(context.getMeanMotion());
  456.             // B / An
  457.             final T BoAn    = context.getBoA().divide(context.getMeanMotion());
  458.             // 1 / ABn
  459.             final T ooABn   = context.getOoAB().divide(context.getMeanMotion());
  460.             // C / 2ABn
  461.             final T Co2ABn  = context.getMCo2AB().negate().divide(context.getMeanMotion());
  462.             // B / (A * (1 + B) * n)
  463.             final T BoABpon = context.getBoABpo().divide(context.getMeanMotion());
  464.             // -3 / n²a² = -3 / nA
  465.             final T m3onA   = context.getA().multiply(context.getMeanMotion()).divide(-3.).reciprocal();

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

  472.                 // Compute the cross derivatives operator :
  473.                 final T SAlphaGammaCj    = context.getAlpha().multiply(gfCoefs.getdSdgammaCj(j)).subtract(context.getGamma().multiply(gfCoefs.getdSdalphaCj(j)));
  474.                 final T SAlphaBetaCj     = context.getAlpha().multiply(gfCoefs.getdSdbetaCj(j)).subtract(context.getBeta().multiply(gfCoefs.getdSdalphaCj(j)));
  475.                 final T SBetaGammaCj     = context.getBeta().multiply(gfCoefs.getdSdgammaCj(j)).subtract(context.getGamma().multiply(gfCoefs.getdSdbetaCj(j)));
  476.                 final T ShkCj            = auxiliaryElements.getH().multiply(gfCoefs.getdSdkCj(j)).subtract(auxiliaryElements.getK().multiply(gfCoefs.getdSdhCj(j)));
  477.                 final T pSagmIqSbgoABnCj = ooABn.multiply(auxiliaryElements.getP().multiply(SAlphaGammaCj).subtract(auxiliaryElements.getQ().multiply(SBetaGammaCj).multiply(I)));
  478.                 final T ShkmSabmdSdlCj   = ShkCj.subtract(SAlphaBetaCj).subtract(gfCoefs.getdSdlambdaCj(j));

  480.                 currentCij[0] = ax2oAn.multiply(gfCoefs.getdSdlambdaCj(j));
  481.                 currentCij[1] = BoAn.multiply(gfCoefs.getdSdhCj(j)).add(auxiliaryElements.getH().multiply(pSagmIqSbgoABnCj)).add(auxiliaryElements.getK().multiply(BoABpon).multiply(gfCoefs.getdSdlambdaCj(j))).negate();
  482.                 currentCij[2] = BoAn.multiply(gfCoefs.getdSdkCj(j)).add(auxiliaryElements.getK().multiply(pSagmIqSbgoABnCj)).subtract(auxiliaryElements.getH().multiply(BoABpon).multiply(gfCoefs.getdSdlambdaCj(j)));
  483.                 currentCij[3] = Co2ABn.multiply(auxiliaryElements.getQ().multiply(ShkmSabmdSdlCj).subtract(SAlphaGammaCj.multiply(I)));
  484.                 currentCij[4] = Co2ABn.multiply(auxiliaryElements.getP().multiply(ShkmSabmdSdlCj).subtract(SBetaGammaCj));
  485.                 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)));

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

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

  493.                 // Compute the cross derivatives operator :
  494.                 final T SAlphaGammaSj    = context.getAlpha().multiply(gfCoefs.getdSdgammaSj(j)).subtract(context.getGamma().multiply(gfCoefs.getdSdalphaSj(j)));
  495.                 final T SAlphaBetaSj     = context.getAlpha().multiply(gfCoefs.getdSdbetaSj(j)).subtract(context.getBeta().multiply(gfCoefs.getdSdalphaSj(j)));
  496.                 final T SBetaGammaSj     = context.getBeta().multiply(gfCoefs.getdSdgammaSj(j)).subtract(context.getGamma().multiply(gfCoefs.getdSdbetaSj(j)));
  497.                 final T ShkSj            = auxiliaryElements.getH().multiply(gfCoefs.getdSdkSj(j)).subtract(auxiliaryElements.getK().multiply(gfCoefs.getdSdhSj(j)));
  498.                 final T pSagmIqSbgoABnSj = ooABn.multiply(auxiliaryElements.getP().multiply(SAlphaGammaSj).subtract(auxiliaryElements.getQ().multiply(SBetaGammaSj).multiply(I)));
  499.                 final T ShkmSabmdSdlSj   = ShkSj.subtract(SAlphaBetaSj).subtract(gfCoefs.getdSdlambdaSj(j));

  501.                 currentSij[0] = ax2oAn.multiply(gfCoefs.getdSdlambdaSj(j));
  502.                 currentSij[1] = BoAn.multiply(gfCoefs.getdSdhSj(j)).add(auxiliaryElements.getH().multiply(pSagmIqSbgoABnSj)).add(auxiliaryElements.getK().multiply(BoABpon).multiply(gfCoefs.getdSdlambdaSj(j))).negate();
  503.                 currentSij[2] = BoAn.multiply(gfCoefs.getdSdkSj(j)).add(auxiliaryElements.getK().multiply(pSagmIqSbgoABnSj)).subtract(auxiliaryElements.getH().multiply(BoABpon).multiply(gfCoefs.getdSdlambdaSj(j)));
  504.                 currentSij[3] = Co2ABn.multiply(auxiliaryElements.getQ().multiply(ShkmSabmdSdlSj).subtract(SAlphaGammaSj.multiply(I)));
  505.                 currentSij[4] = Co2ABn.multiply(auxiliaryElements.getP().multiply(ShkmSabmdSdlSj).subtract(SBetaGammaSj));
  506.                 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)));

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

  511.                 if (j == 1) {
  512.                     //Compute the C⁰ coefficients using Danielson 2.5.2-15a.
  513.                     final T[] value = MathArrays.buildArray(field, 6);
  514.                     for (int i = 0; i < 6; ++i) {
  515.                         value[i] = currentCij[i].multiply(auxiliaryElements.getK()).divide(2.).add(currentSij[i].multiply(auxiliaryElements.getH()).divide(2.));
  516.                     }
  517.                     slot.cij[0].addGridPoint(meanState.getDate(), value);
  518.                 }
  519.             }
  520.         }
  521.     }

  523.     /** {@inheritDoc} */
  524.     @Override
  525.     public EventDetector[] getEventsDetectors() {
  526.         return null;
  527.     }

  529.     /** {@inheritDoc} */
  530.     @Override
  531.     public <T extends RealFieldElement<T>> FieldEventDetector<T>[] getFieldEventsDetectors(final Field<T> field) {
  532.         return null;
  533.     }

  535.     /** {@inheritDoc} */
  536.     @Override
  537.     public void registerAttitudeProvider(final AttitudeProvider provider) {
  538.         //nothing is done since this contribution is not sensitive to attitude
  539.     }

  541.     /** {@inheritDoc} */
  542.     @Override
  543.     public ParameterDriver[] getParametersDrivers() {
  544.         return new ParameterDriver[] {
  545.             bodyParameterDriver,
  546.             gmParameterDriver
  547.         };
  548.     }

  550.     /** Computes the C<sup>j</sup> and S<sup>j</sup> coefficients Danielson 4.2-(15,16)
  551.      * and their derivatives.
  552.      *  <p>
  553.      *  CS Mathematical Report $3.5.3.2
  554.      *  </p>
  555.      */
  556.     private class FourierCjSjCoefficients {

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

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

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

  567.         /** The Fourier coefficients C<sup>j</sup> and their derivatives.
  568.          * <p>
  569.          * Each column of the matrix contains the following values: <br/>
  570.          * - C<sup>j</sup> <br/>
  571.          * - dC<sup>j</sup> / da <br/>
  572.          * - dC<sup>j</sup> / dk <br/>
  573.          * - dC<sup>j</sup> / dh <br/>
  574.          * - dC<sup>j</sup> / dα <br/>
  575.          * - dC<sup>j</sup> / dβ <br/>
  576.          * - dC<sup>j</sup> / dγ <br/>
  577.          * </p>
  578.          */
  579.         private final double[][] cj;

  581.         /** The S<sup>j</sup> coefficients and their derivatives.
  582.          * <p>
  583.          * Each column of the matrix contains the following values: <br/>
  584.          * - S<sup>j</sup> <br/>
  585.          * - dS<sup>j</sup> / da <br/>
  586.          * - dS<sup>j</sup> / dk <br/>
  587.          * - dS<sup>j</sup> / dh <br/>
  588.          * - dS<sup>j</sup> / dα <br/>
  589.          * - dS<sup>j</sup> / dβ <br/>
  590.          * - dS<sup>j</sup> / dγ <br/>
  591.          * </p>
  592.          */
  593.         private final double[][] sj;

  595.         /** The Coefficients C<sup>j</sup><sub>,λ</sub>.
  596.          * <p>
  597.          * See Danielson 4.2-21
  598.          * </p>
  599.          */
  600.         private final double[] cjlambda;

  602.         /** The Coefficients S<sup>j</sup><sub>,λ</sub>.
  603.         * <p>
  604.         * See Danielson 4.2-21
  605.         * </p>
  606.         */
  607.         private final double[] sjlambda;

  609.         /** Maximum value for n. */
  610.         private final int nMax;

  612.         /** Maximum value for s. */
  613.         private final int sMax;

  615.         /** Maximum value for j. */
  616.         private final int jMax;

  618.         /**
  619.          * Private constructor.
  620.          *
  621.          * @param nMax maximum value for n index
  622.          * @param sMax maximum value for s index
  623.          * @param jMax maximum value for j index
  624.          * @param context container for attributes
  625.          */
  626.         FourierCjSjCoefficients(final int nMax, final int sMax, final int jMax, final DSSTThirdBodyContext context) {
  627.             //Save parameters
  628.             this.nMax = nMax;
  629.             this.sMax = sMax;
  630.             this.jMax = jMax;

  632.             //Create objects
  633.             wnsjEtomjmsCoefficient = new WnsjEtomjmsCoefficient(context);
  634.             ABDECoefficients = new CjSjAlphaBetaKH(context);
  635.             gns = new GnsCoefficients(nMax, sMax, context);

  637.             //create arays
  638.             this.cj = new double[7][jMax + 1];
  639.             this.sj = new double[7][jMax + 1];
  640.             this.cjlambda = new double[jMax];
  641.             this.sjlambda = new double[jMax];

  643.             computeCoefficients(context);
  644.         }

  646.         /**
  647.          * Compute all coefficients.
  648.          * @param context container for attributes
  649.          */
  650.         private void computeCoefficients(final DSSTThirdBodyContext context) {

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

  654.             for (int j = 1; j <= jMax; j++) {
  655.                 // initialise the coefficients
  656.                 for (int i = 0; i <= 6; i++) {
  657.                     cj[i][j] = 0.;
  658.                     sj[i][j] = 0.;
  659.                 }
  660.                 if (j < jMax) {
  661.                     // initialise the C<sup>j</sup><sub>,λ</sub> and S<sup>j</sup><sub>,λ</sub> coefficients
  662.                     cjlambda[j] = 0.;
  663.                     sjlambda[j] = 0.;
  664.                 }
  665.                 for (int s = 0; s <= sMax; s++) {

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

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

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

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

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

  686.                             //Compute C<sup>j</sup>
  687.                             cj[0][j] += gns.getGns(n, s) * coef1;
  688.                             //Compute dC<sup>j</sup> / da
  689.                             cj[1][j] += gns.getdGnsda(n, s) * coef1;
  690.                             //Compute dC<sup>j</sup> / dk
  691.                             cj[2][j] += -gns.getGns(n, s) *
  692.                                         (
  693.                                             wjnp1semjms[1] * ABDECoefficients.getCoefA() +
  694.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefAdk() +
  695.                                             wmjnp1semjms[1] * ABDECoefficients.getCoefB() +
  696.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefBdk()
  697.                                          );
  698.                             //Compute dC<sup>j</sup> / dh
  699.                             cj[3][j] += -gns.getGns(n, s) *
  700.                                         (
  701.                                             wjnp1semjms[2] * ABDECoefficients.getCoefA() +
  702.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefAdh() +
  703.                                             wmjnp1semjms[2] * ABDECoefficients.getCoefB() +
  704.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefBdh()
  705.                                          );
  706.                             //Compute dC<sup>j</sup> / dα
  707.                             cj[4][j] += -gns.getGns(n, s) *
  708.                                         (
  709.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefAdalpha() +
  710.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefBdalpha()
  711.                                         );
  712.                             //Compute dC<sup>j</sup> / dβ
  713.                             cj[5][j] += -gns.getGns(n, s) *
  714.                                         (
  715.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefAdbeta() +
  716.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefBdbeta()
  717.                                         );
  718.                             //Compute dC<sup>j</sup> / dγ
  719.                             cj[6][j] += gns.getdGnsdgamma(n, s) * coef1;

  721.                             //Compute S<sup>j</sup>
  722.                             sj[0][j] += gns.getGns(n, s) * coef2;
  723.                             //Compute dS<sup>j</sup> / da
  724.                             sj[1][j] += gns.getdGnsda(n, s) * coef2;
  725.                             //Compute dS<sup>j</sup> / dk
  726.                             sj[2][j] += gns.getGns(n, s) *
  727.                                         (
  728.                                             wjnp1semjms[1] * ABDECoefficients.getCoefD() +
  729.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefDdk() +
  730.                                             wmjnp1semjms[1] * ABDECoefficients.getCoefE() +
  731.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefEdk()
  732.                                          );
  733.                             //Compute dS<sup>j</sup> / dh
  734.                             sj[3][j] += gns.getGns(n, s) *
  735.                                         (
  736.                                             wjnp1semjms[2] * ABDECoefficients.getCoefD() +
  737.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefDdh() +
  738.                                             wmjnp1semjms[2] * ABDECoefficients.getCoefE() +
  739.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefEdh()
  740.                                          );
  741.                             //Compute dS<sup>j</sup> / dα
  742.                             sj[4][j] += gns.getGns(n, s) *
  743.                                         (
  744.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefDdalpha() +
  745.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefEdalpha()
  746.                                         );
  747.                             //Compute dS<sup>j</sup> / dβ
  748.                             sj[5][j] += gns.getGns(n, s) *
  749.                                         (
  750.                                             wjnp1semjms[0] * ABDECoefficients.getdCoefDdbeta() +
  751.                                             wmjnp1semjms[0] * ABDECoefficients.getdCoefEdbeta()
  752.                                         );
  753.                             //Compute dS<sup>j</sup> / dγ
  754.                             sj[6][j] += gns.getdGnsdgamma(n, s) * coef2;

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

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

  764.                                 //Compute C<sup>j</sup><sub>,λ</sub>
  765.                                 cjlambda[j] += gns.getGns(n, s) * (wjnsemjms[0] * ABDECoefficients.getCoefD() + wmjnsemjms[0] * ABDECoefficients.getCoefE());
  766.                                 //Compute S<sup>j</sup><sub>,λ</sub>
  767.                                 sjlambda[j] += gns.getGns(n, s) * (wjnsemjms[0] * ABDECoefficients.getCoefA() + wmjnsemjms[0] * ABDECoefficients.getCoefB());
  768.                             }
  769.                         }
  770.                     }
  771.                 }
  772.                 // Divide by j
  773.                 for (int i = 0; i <= 6; i++) {
  774.                     cj[i][j] /= j;
  775.                     sj[i][j] /= j;
  776.                 }
  777.             }
  778.             //The C⁰ coefficients are not computed here.
  779.             //They are evaluated at the final point.

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

  785.         /** Get the Fourier coefficient C<sup>j</sup>.
  786.          * @param j j index
  787.          * @return C<sup>j</sup>
  788.          */
  789.         public double getCj(final int j) {
  790.             return cj[0][j];
  791.         }

  793.         /** Get the derivative dC<sup>j</sup>/da.
  794.          * @param j j index
  795.          * @return dC<sup>j</sup>/da
  796.          */
  797.         public double getdCjda(final int j) {
  798.             return cj[1][j];
  799.         }

  801.         /** Get the derivative dC<sup>j</sup>/dk.
  802.          * @param j j index
  803.          * @return dC<sup>j</sup>/dk
  804.          */
  805.         public double getdCjdk(final int j) {
  806.             return cj[2][j];
  807.         }

  809.         /** Get the derivative dC<sup>j</sup>/dh.
  810.          * @param j j index
  811.          * @return dC<sup>j</sup>/dh
  812.          */
  813.         public double getdCjdh(final int j) {
  814.             return cj[3][j];
  815.         }

  817.         /** Get the derivative dC<sup>j</sup>/dα.
  818.          * @param j j index
  819.          * @return dC<sup>j</sup>/dα
  820.          */
  821.         public double getdCjdalpha(final int j) {
  822.             return cj[4][j];
  823.         }

  825.         /** Get the derivative dC<sup>j</sup>/dβ.
  826.          * @param j j index
  827.          * @return dC<sup>j</sup>/dβ
  828.          */
  829.         public double getdCjdbeta(final int j) {
  830.             return cj[5][j];
  831.         }

  833.         /** Get the derivative dC<sup>j</sup>/dγ.
  834.          * @param j j index
  835.          * @return dC<sup>j</sup>/dγ
  836.          */
  837.         public double getdCjdgamma(final int j) {
  838.             return cj[6][j];
  839.         }

  841.         /** Get the Fourier coefficient S<sup>j</sup>.
  842.          * @param j j index
  843.          * @return S<sup>j</sup>
  844.          */
  845.         public double getSj(final int j) {
  846.             return sj[0][j];
  847.         }

  849.         /** Get the derivative dS<sup>j</sup>/da.
  850.          * @param j j index
  851.          * @return dS<sup>j</sup>/da
  852.          */
  853.         public double getdSjda(final int j) {
  854.             return sj[1][j];
  855.         }

  857.         /** Get the derivative dS<sup>j</sup>/dk.
  858.          * @param j j index
  859.          * @return dS<sup>j</sup>/dk
  860.          */
  861.         public double getdSjdk(final int j) {
  862.             return sj[2][j];
  863.         }

  865.         /** Get the derivative dS<sup>j</sup>/dh.
  866.          * @param j j index
  867.          * @return dS<sup>j</sup>/dh
  868.          */
  869.         public double getdSjdh(final int j) {
  870.             return sj[3][j];
  871.         }

  873.         /** Get the derivative dS<sup>j</sup>/dα.
  874.          * @param j j index
  875.          * @return dS<sup>j</sup>/dα
  876.          */
  877.         public double getdSjdalpha(final int j) {
  878.             return sj[4][j];
  879.         }

  881.         /** Get the derivative dS<sup>j</sup>/dβ.
  882.          * @param j j index
  883.          * @return dS<sup>j</sup>/dβ
  884.          */
  885.         public double getdSjdbeta(final int j) {
  886.             return sj[5][j];
  887.         }

  889.         /** Get the derivative dS<sup>j</sup>/dγ.
  890.          * @param j j index
  891.          * @return dS<sup>j</sup>/dγ
  892.          */
  893.         public double getdSjdgamma(final int j) {
  894.             return sj[6][j];
  895.         }

  897.         /** Get the coefficient C⁰<sub>,λ</sub>.
  898.          * @return C⁰<sub>,λ</sub>
  899.          */
  900.         public double getC0Lambda() {
  901.             return cjlambda[0];
  902.         }

  904.         /** Get the coefficient C<sup>j</sup><sub>,λ</sub>.
  905.          * @param j j index
  906.          * @return C<sup>j</sup><sub>,λ</sub>
  907.          */
  908.         public double getCjLambda(final int j) {
  909.             if (j < 1 || j >= jMax) {
  910.                 return 0.;
  911.             }
  912.             return cjlambda[j];
  913.         }

  915.         /** Get the coefficient S<sup>j</sup><sub>,λ</sub>.
  916.          * @param j j index
  917.          * @return S<sup>j</sup><sub>,λ</sub>
  918.          */
  919.         public double getSjLambda(final int j) {
  920.             if (j < 1 || j >= jMax) {
  921.                 return 0.;
  922.             }
  923.             return sjlambda[j];
  924.         }
  925.     }

  927.     /** Computes the C<sup>j</sup> and S<sup>j</sup> coefficients Danielson 4.2-(15,16)
  928.      * and their derivatives.
  929.      *  <p>
  930.      *  CS Mathematical Report $3.5.3.2
  931.      *  </p>
  932.      */
  933.     private class FieldFourierCjSjCoefficients <T extends RealFieldElement<T>> {

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

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

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

  944.         /** The Fourier coefficients C<sup>j</sup> and their derivatives.
  945.          * <p>
  946.          * Each column of the matrix contains the following values: <br/>
  947.          * - C<sup>j</sup> <br/>
  948.          * - dC<sup>j</sup> / da <br/>
  949.          * - dC<sup>j</sup> / dk <br/>
  950.          * - dC<sup>j</sup> / dh <br/>
  951.          * - dC<sup>j</sup> / dα <br/>
  952.          * - dC<sup>j</sup> / dβ <br/>
  953.          * - dC<sup>j</sup> / dγ <br/>
  954.          * </p>
  955.          */
  956.         private final T[][] cj;

  958.         /** The S<sup>j</sup> coefficients and their derivatives.
  959.          * <p>
  960.          * Each column of the matrix contains the following values: <br/>
  961.          * - S<sup>j</sup> <br/>
  962.          * - dS<sup>j</sup> / da <br/>
  963.          * - dS<sup>j</sup> / dk <br/>
  964.          * - dS<sup>j</sup> / dh <br/>
  965.          * - dS<sup>j</sup> / dα <br/>
  966.          * - dS<sup>j</sup> / dβ <br/>
  967.          * - dS<sup>j</sup> / dγ <br/>
  968.          * </p>
  969.          */
  970.         private final T[][] sj;

  972.         /** The Coefficients C<sup>j</sup><sub>,λ</sub>.
  973.          * <p>
  974.          * See Danielson 4.2-21
  975.          * </p>
  976.          */
  977.         private final T[] cjlambda;

  979.         /** The Coefficients S<sup>j</sup><sub>,λ</sub>.
  980.         * <p>
  981.         * See Danielson 4.2-21
  982.         * </p>
  983.         */
  984.         private final T[] sjlambda;

  986.         /** Zero. */
  987.         private final T zero;

  989.         /** Maximum value for n. */
  990.         private final int nMax;

  992.         /** Maximum value for s. */
  993.         private final int sMax;

  995.         /** Maximum value for j. */
  996.         private final int jMax;

  998.         /**
  999.          * Private constructor.
  1000.          *
  1001.          * @param nMax maximum value for n index
  1002.          * @param sMax maximum value for s index
  1003.          * @param jMax maximum value for j index
  1004.          * @param context container for attributes
  1005.          * @param field field used by default
  1006.          */
  1007.         FieldFourierCjSjCoefficients(final int nMax, final int sMax, final int jMax,
  1008.                                      final FieldDSSTThirdBodyContext<T> context,
  1009.                                      final Field<T> field) {
  1010.             //Zero
  1011.             this.zero = field.getZero();

  1013.             //Save parameters
  1014.             this.nMax = nMax;
  1015.             this.sMax = sMax;
  1016.             this.jMax = jMax;

  1018.             //Create objects
  1019.             wnsjEtomjmsCoefficient = new FieldWnsjEtomjmsCoefficient<>(context, field);
  1020.             ABDECoefficients       = new FieldCjSjAlphaBetaKH<>(context, field);
  1021.             gns                    = new FieldGnsCoefficients<>(nMax, sMax, context, field);

  1023.             //create arays
  1024.             this.cj = MathArrays.buildArray(field, 7, jMax + 1);
  1025.             this.sj = MathArrays.buildArray(field, 7, jMax + 1);
  1026.             this.cjlambda = MathArrays.buildArray(field, jMax);
  1027.             this.sjlambda = MathArrays.buildArray(field, jMax);

  1029.             computeCoefficients(context, field);
  1030.         }

  1032.         /**
  1033.          * Compute all coefficients.
  1034.          * @param context container for attributes
  1035.          * @param field field used by default
  1036.          */
  1037.         private void computeCoefficients(final FieldDSSTThirdBodyContext<T> context,
  1038.                                          final Field<T> field) {

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

  1042.             for (int j = 1; j <= jMax; j++) {
  1043.                 // initialise the coefficients
  1044.                 for (int i = 0; i <= 6; i++) {
  1045.                     cj[i][j] = zero;
  1046.                     sj[i][j] = zero;
  1047.                 }
  1048.                 if (j < jMax) {
  1049.                     // initialise the C<sup>j</sup><sub>,λ</sub> and S<sup>j</sup><sub>,λ</sub> coefficients
  1050.                     cjlambda[j] = zero;
  1051.                     sjlambda[j] = zero;
  1052.                 }
  1053.                 for (int s = 0; s <= sMax; s++) {

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

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

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

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

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

  1074.                             //Compute C<sup>j</sup>
  1075.                             cj[0][j] = cj[0][j].add(gns.getGns(n, s).multiply(coef1));
  1076.                             //Compute dC<sup>j</sup> / da
  1077.                             cj[1][j] = cj[1][j].add(gns.getdGnsda(n, s).multiply(coef1));
  1078.                             //Compute dC<sup>j</sup> / dk
  1079.                             cj[2][j] = cj[2][j].add(gns.getGns(n, s).negate().
  1080.                                        multiply(
  1081.                                             wjnp1semjms[1].multiply(ABDECoefficients.getCoefA()).
  1082.                                             add(wjnp1semjms[0].multiply(ABDECoefficients.getdCoefAdk())).
  1083.                                             add(wmjnp1semjms[1].multiply(ABDECoefficients.getCoefB())).
  1084.                                             add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefBdk()))
  1085.                                          ));
  1086.                             //Compute dC<sup>j</sup> / dh
  1087.                             cj[3][j] = cj[3][j].add(gns.getGns(n, s).negate().
  1088.                                        multiply(
  1089.                                              wjnp1semjms[2].multiply(ABDECoefficients.getCoefA()).
  1090.                                              add(wjnp1semjms[0].multiply(ABDECoefficients.getdCoefAdh())).
  1091.                                              add(wmjnp1semjms[2].multiply(ABDECoefficients.getCoefB())).
  1092.                                              add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefBdh()))
  1093.                                              ));
  1094.                             //Compute dC<sup>j</sup> / dα
  1095.                             cj[4][j] = cj[4][j].add(gns.getGns(n, s).negate().
  1096.                                        multiply(
  1097.                                            wjnp1semjms[0].multiply(ABDECoefficients.getdCoefAdalpha()).
  1098.                                            add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefBdalpha()))
  1099.                                        ));
  1100.                             //Compute dC<sup>j</sup> / dβ
  1101.                             cj[5][j] = cj[5][j].add(gns.getGns(n, s).negate().
  1102.                                        multiply(
  1103.                                            wjnp1semjms[0].multiply(ABDECoefficients.getdCoefAdbeta()).
  1104.                                            add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefBdbeta()))
  1105.                                        ));
  1106.                             //Compute dC<sup>j</sup> / dγ
  1107.                             cj[6][j] = cj[6][j].add(gns.getdGnsdgamma(n, s).multiply(coef1));

  1109.                             //Compute S<sup>j</sup>
  1110.                             sj[0][j] = sj[0][j].add(gns.getGns(n, s).multiply(coef2));
  1111.                             //Compute dS<sup>j</sup> / da
  1112.                             sj[1][j] = sj[1][j].add(gns.getdGnsda(n, s).multiply(coef2));
  1113.                             //Compute dS<sup>j</sup> / dk
  1114.                             sj[2][j] = sj[2][j].add(gns.getGns(n, s).
  1115.                                        multiply(
  1116.                                            wjnp1semjms[1].multiply(ABDECoefficients.getCoefD()).
  1117.                                            add(wjnp1semjms[0].multiply(ABDECoefficients.getdCoefDdk())).
  1118.                                            add(wmjnp1semjms[1].multiply(ABDECoefficients.getCoefE())).
  1119.                                            add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefEdk()))
  1120.                                        ));
  1121.                             //Compute dS<sup>j</sup> / dh
  1122.                             sj[3][j] = sj[3][j].add(gns.getGns(n, s).
  1123.                                        multiply(
  1124.                                            wjnp1semjms[2].multiply(ABDECoefficients.getCoefD()).
  1125.                                            add(wjnp1semjms[0].multiply(ABDECoefficients.getdCoefDdh())).
  1126.                                            add(wmjnp1semjms[2].multiply(ABDECoefficients.getCoefE())).
  1127.                                            add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefEdh()))
  1128.                                        ));
  1129.                             //Compute dS<sup>j</sup> / dα
  1130.                             sj[4][j] = sj[4][j].add(gns.getGns(n, s).
  1131.                                        multiply(
  1132.                                             wjnp1semjms[0].multiply(ABDECoefficients.getdCoefDdalpha()).
  1133.                                             add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefEdalpha()))
  1134.                                        ));
  1135.                             //Compute dS<sup>j</sup> / dβ
  1136.                             sj[5][j] = sj[5][j].add(gns.getGns(n, s).
  1137.                                         multiply(
  1138.                                             wjnp1semjms[0].multiply(ABDECoefficients.getdCoefDdbeta()).
  1139.                                             add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefEdbeta()))
  1140.                                        ));
  1141.                             //Compute dS<sup>j</sup> / dγ
  1142.                             sj[6][j] = sj[6][j].add(gns.getdGnsdgamma(n, s).multiply(coef2));

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

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

  1152.                                 //Compute C<sup>j</sup><sub>,λ</sub>
  1153.                                 cjlambda[j] = cjlambda[j].add(gns.getGns(n, s).multiply(wjnsemjms[0].multiply(ABDECoefficients.getCoefD()).add(wmjnsemjms[0].multiply(ABDECoefficients.getCoefE()))));
  1154.                                 //Compute S<sup>j</sup><sub>,λ</sub>
  1155.                                 sjlambda[j] = sjlambda[j].add(gns.getGns(n, s).multiply(wjnsemjms[0].multiply(ABDECoefficients.getCoefA()).add(wmjnsemjms[0].multiply(ABDECoefficients.getCoefB()))));
  1156.                             }
  1157.                         }
  1158.                     }
  1159.                 }
  1160.                 // Divide by j
  1161.                 for (int i = 0; i <= 6; i++) {
  1162.                     cj[i][j] = cj[i][j].divide(j);
  1163.                     sj[i][j] = sj[i][j].divide(j);
  1164.                 }
  1165.             }
  1166.             //The C⁰ coefficients are not computed here.
  1167.             //They are evaluated at the final point.

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

  1173.         /** Get the Fourier coefficient C<sup>j</sup>.
  1174.          * @param j j index
  1175.          * @return C<sup>j</sup>
  1176.          */
  1177.         public T getCj(final int j) {
  1178.             return cj[0][j];
  1179.         }

  1181.         /** Get the derivative dC<sup>j</sup>/da.
  1182.          * @param j j index
  1183.          * @return dC<sup>j</sup>/da
  1184.          */
  1185.         public T getdCjda(final int j) {
  1186.             return cj[1][j];
  1187.         }

  1189.         /** Get the derivative dC<sup>j</sup>/dk.
  1190.          * @param j j index
  1191.          * @return dC<sup>j</sup>/dk
  1192.          */
  1193.         public T getdCjdk(final int j) {
  1194.             return cj[2][j];
  1195.         }

  1197.         /** Get the derivative dC<sup>j</sup>/dh.
  1198.          * @param j j index
  1199.          * @return dC<sup>j</sup>/dh
  1200.          */
  1201.         public T getdCjdh(final int j) {
  1202.             return cj[3][j];
  1203.         }

  1205.         /** Get the derivative dC<sup>j</sup>/dα.
  1206.          * @param j j index
  1207.          * @return dC<sup>j</sup>/dα
  1208.          */
  1209.         public T getdCjdalpha(final int j) {
  1210.             return cj[4][j];
  1211.         }

  1213.         /** Get the derivative dC<sup>j</sup>/dβ.
  1214.          * @param j j index
  1215.          * @return dC<sup>j</sup>/dβ
  1216.          */
  1217.         public T getdCjdbeta(final int j) {
  1218.             return cj[5][j];
  1219.         }

  1221.         /** Get the derivative dC<sup>j</sup>/dγ.
  1222.          * @param j j index
  1223.          * @return dC<sup>j</sup>/dγ
  1224.          */
  1225.         public T getdCjdgamma(final int j) {
  1226.             return cj[6][j];
  1227.         }

  1229.         /** Get the Fourier coefficient S<sup>j</sup>.
  1230.          * @param j j index
  1231.          * @return S<sup>j</sup>
  1232.          */
  1233.         public T getSj(final int j) {
  1234.             return sj[0][j];
  1235.         }

  1237.         /** Get the derivative dS<sup>j</sup>/da.
  1238.          * @param j j index
  1239.          * @return dS<sup>j</sup>/da
  1240.          */
  1241.         public T getdSjda(final int j) {
  1242.             return sj[1][j];
  1243.         }

  1245.         /** Get the derivative dS<sup>j</sup>/dk.
  1246.          * @param j j index
  1247.          * @return dS<sup>j</sup>/dk
  1248.          */
  1249.         public T getdSjdk(final int j) {
  1250.             return sj[2][j];
  1251.         }

  1253.         /** Get the derivative dS<sup>j</sup>/dh.
  1254.          * @param j j index
  1255.          * @return dS<sup>j</sup>/dh
  1256.          */
  1257.         public T getdSjdh(final int j) {
  1258.             return sj[3][j];
  1259.         }

  1261.         /** Get the derivative dS<sup>j</sup>/dα.
  1262.          * @param j j index
  1263.          * @return dS<sup>j</sup>/dα
  1264.          */
  1265.         public T getdSjdalpha(final int j) {
  1266.             return sj[4][j];
  1267.         }

  1269.         /** Get the derivative dS<sup>j</sup>/dβ.
  1270.          * @param j j index
  1271.          * @return dS<sup>j</sup>/dβ
  1272.          */
  1273.         public T getdSjdbeta(final int j) {
  1274.             return sj[5][j];
  1275.         }

  1277.         /** Get the derivative dS<sup>j</sup>/dγ.
  1278.          * @param j j index
  1279.          * @return dS<sup>j</sup>/dγ
  1280.          */
  1281.         public T getdSjdgamma(final int j) {
  1282.             return sj[6][j];
  1283.         }

  1285.         /** Get the coefficient C⁰<sub>,λ</sub>.
  1286.          * @return C⁰<sub>,λ</sub>
  1287.          */
  1288.         public T getC0Lambda() {
  1289.             return cjlambda[0];
  1290.         }

  1292.         /** Get the coefficient C<sup>j</sup><sub>,λ</sub>.
  1293.          * @param j j index
  1294.          * @return C<sup>j</sup><sub>,λ</sub>
  1295.          */
  1296.         public T getCjLambda(final int j) {
  1297.             if (j < 1 || j >= jMax) {
  1298.                 return zero;
  1299.             }
  1300.             return cjlambda[j];
  1301.         }

  1303.         /** Get the coefficient S<sup>j</sup><sub>,λ</sub>.
  1304.          * @param j j index
  1305.          * @return S<sup>j</sup><sub>,λ</sub>
  1306.          */
  1307.         public T getSjLambda(final int j) {
  1308.             if (j < 1 || j >= jMax) {
  1309.                 return zero;
  1310.             }
  1311.             return sjlambda[j];
  1312.         }
  1313.     }

  1315.     /** 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.
  1316.      *
  1317.      * <p>
  1318.      * Starting from Danielson 4.2-9,10,11 and taking into account that fact that: <br />
  1319.      * c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
  1320.      * the expression e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup>
  1321.      * can be written as: <br/ >
  1322.      * - for |s| > |j| <br />
  1323.      * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
  1324.      *          (((n + s)!(n - s)!)/((n + j)!(n - j)!)) *
  1325.      *          (-b)<sup>|j-s|</sup> *
  1326.      *          ((1 - c²)<sup>n-|s|</sup>/(1 + c²)<sup>n</sup>) *
  1327.      *          P<sub>n-|s|</sub><sup>|j-s|, |j+s|</sup>(χ) <br />
  1328.      * <br />
  1329.      * - for |s| <= |j| <br />
  1330.      * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
  1331.      *          (-b)<sup>|j-s|</sup> *
  1332.      *          ((1 - c²)<sup>n-|j|</sup>/(1 + c²)<sup>n</sup>) *
  1333.      *          P<sub>n-|j|</sub><sup>|j-s|, |j+s|</sup>(χ)
  1334.      * </p>
  1335.      *
  1336.      * @author Lucian Barbulescu
  1337.      */
  1338.     private static class WnsjEtomjmsCoefficient {

  1340.         /** The value c.
  1341.          * <p>
  1342.          *  c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
  1343.          * </p>
  1344.          *  */
  1345.         private final double c;

  1347.         /** db / dh. */
  1348.         private final double dbdh;

  1350.         /** db / dk. */
  1351.         private final double dbdk;

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

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

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

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

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

  1370.         /**
  1371.          * Standard constructor.
  1372.          * @param context container for attributes
  1373.          */
  1374.         WnsjEtomjmsCoefficient(final DSSTThirdBodyContext context) {

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

  1378.             //initialise fields
  1379.             c = auxiliaryElements.getEcc() * context.getb();
  1380.             final double c2 = c * c;

  1382.             //b² * χ
  1383.             final double b2Chi = context.getb() * context.getb() * context.getX();
  1384.             //Compute derivatives of b
  1385.             dbdh = auxiliaryElements.getH() * b2Chi;
  1386.             dbdk = auxiliaryElements.getK() * b2Chi;

  1388.             //Compute derivatives of c
  1389.             if (auxiliaryElements.getEcc() == 0.0) {
  1390.                 // we are at a perfectly circular orbit singularity here
  1391.                 // we arbitrarily consider the perigee is along the X axis,
  1392.                 // i.e cos(ω + Ω) = h/ecc 1 and sin(ω + Ω) = k/ecc = 0
  1393.                 dcdh = auxiliaryElements.getEcc() * dbdh + context.getb();
  1394.                 dcdk = auxiliaryElements.getEcc() * dbdk;
  1395.             } else {
  1396.                 dcdh = auxiliaryElements.getEcc() * dbdh + context.getb() * auxiliaryElements.getH() / auxiliaryElements.getEcc();
  1397.                 dcdk = auxiliaryElements.getEcc() * dbdk + context.getb() * auxiliaryElements.getK() / auxiliaryElements.getEcc();
  1398.             }

  1400.             //Compute the powers (1 - c²)<sup>n</sup> and (1 + c²)<sup>n</sup>
  1401.             omc2tn = new double[context.getMaxAR3Pow() + context.getMaxFreqF() + 2];
  1402.             opc2tn = new double[context.getMaxAR3Pow() + context.getMaxFreqF() + 2];
  1403.             final double omc2 = 1. - c2;
  1404.             final double opc2 = 1. + c2;
  1405.             omc2tn[0] = 1.;
  1406.             opc2tn[0] = 1.;
  1407.             for (int i = 1; i <= context.getMaxAR3Pow() + context.getMaxFreqF() + 1; i++) {
  1408.                 omc2tn[i] = omc2tn[i - 1] * omc2;
  1409.                 opc2tn[i] = opc2tn[i - 1] * opc2;
  1410.             }

  1412.             //Compute the powers of b
  1413.             btjms = new double[context.getMaxAR3Pow() + context.getMaxFreqF() + 1];
  1414.             btjms[0] = 1.;
  1415.             for (int i = 1; i <= context.getMaxAR3Pow() + context.getMaxFreqF(); i++) {
  1416.                 btjms[i] = btjms[i - 1] * context.getb();
  1417.             }
  1418.         }

  1420.         /** 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 />
  1421.          *
  1422.          * @param j j index
  1423.          * @param s s index
  1424.          * @param n n index
  1425.          * @param context container for attributes
  1426.          * @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
  1427.          */
  1428.         public double[] computeWjnsEmjmsAndDeriv(final int j, final int s, final int n, final DSSTThirdBodyContext context) {
  1429.             final double[] wjnsemjms = new double[] {0., 0., 0.};

  1431.             // |j|
  1432.             final int absJ = FastMath.abs(j);
  1433.             // |s|
  1434.             final int absS = FastMath.abs(s);
  1435.             // |j - s|
  1436.             final int absJmS = FastMath.abs(j - s);
  1437.             // |j + s|
  1438.             final int absJpS = FastMath.abs(j + s);

  1440.             //The lower index of P. Also the power of (1 - c²)
  1441.             final int l;
  1442.             // the factorial ratio coefficient or 1. if |s| <= |j|
  1443.             final double factCoef;
  1444.             if (absS > absJ) {
  1445.                 //factCoef = (fact[n + s] / fact[n + j]) * (fact[n - s] / fact[n - j]);
  1446.                 factCoef = (CombinatoricsUtils.factorialDouble(n + s) / CombinatoricsUtils.factorialDouble(n + j)) * (CombinatoricsUtils.factorialDouble(n - s) / CombinatoricsUtils.factorialDouble(n - j));
  1447.                 l = n - absS;
  1448.             } else {
  1449.                 factCoef = 1.;
  1450.                 l = n - absJ;
  1451.             }

  1453.             // (-1)<sup>|j-s|</sup>
  1454.             final double sign = absJmS % 2 != 0 ? -1. : 1.;
  1455.             //(1 - c²)<sup>n-|s|</sup> / (1 + c²)<sup>n</sup>
  1456.             final double coef1 = omc2tn[l] / opc2tn[n];
  1457.             //-b<sup>|j-s|</sup>
  1458.             final double coef2 = sign * btjms[absJmS];
  1459.             // P<sub>l</sub><sup>|j-s|, |j+s|</sup>(χ)
  1460.             final Gradient jac =
  1461.                     JacobiPolynomials.getValue(l, absJmS, absJpS, Gradient.variable(1, 0, context.getX()));

  1463.             // the derivative of coef1 by c
  1464.             final double dcoef1dc = -coef1 * 2. * c * (((double) n) / opc2tn[1] + ((double) l) / omc2tn[1]);
  1465.             // the derivative of coef1 by h
  1466.             final double dcoef1dh = dcoef1dc * dcdh;
  1467.             // the derivative of coef1 by k
  1468.             final double dcoef1dk = dcoef1dc * dcdk;

  1470.             // the derivative of coef2 by b
  1471.             final double dcoef2db = absJmS == 0 ? 0 : sign * (double) absJmS * btjms[absJmS - 1];
  1472.             // the derivative of coef2 by h
  1473.             final double dcoef2dh = dcoef2db * dbdh;
  1474.             // the derivative of coef2 by k
  1475.             final double dcoef2dk = dcoef2db * dbdk;

  1477.             // the jacobi polynomial value
  1478.             final double jacobi = jac.getValue();
  1479.             // the derivative of the Jacobi polynomial by h
  1480.             final double djacobidh = jac.getGradient()[0] * context.getHXXX();
  1481.             // the derivative of the Jacobi polynomial by k
  1482.             final double djacobidk = jac.getGradient()[0] * context.getKXXX();

  1484.             //group the above coefficients to limit the mathematical operations
  1485.             final double term1 = factCoef * coef1 * coef2;
  1486.             final double term2 = factCoef * coef1 * jacobi;
  1487.             final double term3 = factCoef * coef2 * jacobi;

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

  1494.             return wjnsemjms;
  1495.         }
  1496.     }

  1498.     /** 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.
  1499.     *
  1500.     * <p>
  1501.     * Starting from Danielson 4.2-9,10,11 and taking into account that fact that: <br />
  1502.     * c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
  1503.     * the expression e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup>
  1504.     * can be written as: <br/ >
  1505.     * - for |s| > |j| <br />
  1506.     * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
  1507.     *          (((n + s)!(n - s)!)/((n + j)!(n - j)!)) *
  1508.     *          (-b)<sup>|j-s|</sup> *
  1509.     *          ((1 - c²)<sup>n-|s|</sup>/(1 + c²)<sup>n</sup>) *
  1510.     *          P<sub>n-|s|</sub><sup>|j-s|, |j+s|</sup>(χ) <br />
  1511.     * <br />
  1512.     * - for |s| <= |j| <br />
  1513.     * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
  1514.     *          (-b)<sup>|j-s|</sup> *
  1515.     *          ((1 - c²)<sup>n-|j|</sup>/(1 + c²)<sup>n</sup>) *
  1516.     *          P<sub>n-|j|</sub><sup>|j-s|, |j+s|</sup>(χ)
  1517.     * </p>
  1518.     *
  1519.     * @author Lucian Barbulescu
  1520.     */
  1521.     private static class FieldWnsjEtomjmsCoefficient <T extends RealFieldElement<T>> {

  1523.         /** The value c.
  1524.          * <p>
  1525.          *  c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
  1526.          * </p>
  1527.          *  */
  1528.         private final T c;

  1530.         /** db / dh. */
  1531.         private final T dbdh;

  1533.         /** db / dk. */
  1534.         private final T dbdk;

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

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

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

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

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

  1553.         /**
  1554.          * Standard constructor.
  1555.          * @param context container for attributes
  1556.          * @param field field used by default
  1557.          */
  1558.         FieldWnsjEtomjmsCoefficient(final FieldDSSTThirdBodyContext<T> context, final Field<T> field) {

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

  1562.             //Zero
  1563.             final T zero = field.getZero();

  1565.             //initialise fields
  1566.             c = auxiliaryElements.getEcc().multiply(context.getb());
  1567.             final T c2 = c.multiply(c);

  1569.             //b² * χ
  1570.             final T b2Chi = context.getb().multiply(context.getb()).multiply(context.getX());
  1571.             //Compute derivatives of b
  1572.             dbdh = auxiliaryElements.getH().multiply(b2Chi);
  1573.             dbdk = auxiliaryElements.getK().multiply(b2Chi);

  1575.             //Compute derivatives of c
  1576.             if (auxiliaryElements.getEcc().getReal() == 0.0) {
  1577.                 // we are at a perfectly circular orbit singularity here
  1578.                 // we arbitrarily consider the perigee is along the X axis,
  1579.                 // i.e cos(ω + Ω) = h/ecc 1 and sin(ω + Ω) = k/ecc = 0
  1580.                 dcdh = auxiliaryElements.getEcc().multiply(dbdh).add(context.getb());
  1581.                 dcdk = auxiliaryElements.getEcc().multiply(dbdk);
  1582.             } else {
  1583.                 dcdh = auxiliaryElements.getEcc().multiply(dbdh).add(context.getb().multiply(auxiliaryElements.getH()).divide(auxiliaryElements.getEcc()));
  1584.                 dcdk = auxiliaryElements.getEcc().multiply(dbdk).add(context.getb().multiply(auxiliaryElements.getK()).divide(auxiliaryElements.getEcc()));
  1585.             }

  1587.             //Compute the powers (1 - c²)<sup>n</sup> and (1 + c²)<sup>n</sup>
  1588.             omc2tn = MathArrays.buildArray(field, context.getMaxAR3Pow() + context.getMaxFreqF() + 2);
  1589.             opc2tn = MathArrays.buildArray(field, context.getMaxAR3Pow() + context.getMaxFreqF() + 2);
  1590.             final T omc2 = c2.negate().add(1.);
  1591.             final T opc2 = c2.add(1.);
  1592.             omc2tn[0] = zero.add(1.);
  1593.             opc2tn[0] = zero.add(1.);
  1594.             for (int i = 1; i <= context.getMaxAR3Pow() + context.getMaxFreqF() + 1; i++) {
  1595.                 omc2tn[i] = omc2tn[i - 1].multiply(omc2);
  1596.                 opc2tn[i] = opc2tn[i - 1].multiply(opc2);
  1597.             }

  1599.             //Compute the powers of b
  1600.             btjms = MathArrays.buildArray(field, context.getMaxAR3Pow() + context.getMaxFreqF() + 1);
  1601.             btjms[0] = zero.add(1.);
  1602.             for (int i = 1; i <= context.getMaxAR3Pow() + context.getMaxFreqF(); i++) {
  1603.                 btjms[i] = btjms[i - 1].multiply(context.getb());
  1604.             }
  1605.         }

  1607.         /** 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 />
  1608.          *
  1609.          * @param j j index
  1610.          * @param s s index
  1611.          * @param n n index
  1612.          * @param context container for attributes
  1613.          * @param field field used by default
  1614.          * @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
  1615.          */
  1616.         public T[] computeWjnsEmjmsAndDeriv(final int j, final int s, final int n,
  1617.                                             final FieldDSSTThirdBodyContext<T> context,
  1618.                                             final Field<T> field) {
  1619.             //Zero
  1620.             final T zero = field.getZero();

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

  1625.             // |j|
  1626.             final int absJ = FastMath.abs(j);
  1627.             // |s|
  1628.             final int absS = FastMath.abs(s);
  1629.             // |j - s|
  1630.             final int absJmS = FastMath.abs(j - s);
  1631.             // |j + s|
  1632.             final int absJpS = FastMath.abs(j + s);

  1634.             //The lower index of P. Also the power of (1 - c²)
  1635.             final int l;
  1636.             // the factorial ratio coefficient or 1. if |s| <= |j|
  1637.             final T factCoef;
  1638.             if (absS > absJ) {
  1639.                 //factCoef = (fact[n + s] / fact[n + j]) * (fact[n - s] / fact[n - j]);
  1640.                 factCoef = zero.add((CombinatoricsUtils.factorialDouble(n + s) / CombinatoricsUtils.factorialDouble(n + j)) * (CombinatoricsUtils.factorialDouble(n - s) / CombinatoricsUtils.factorialDouble(n - j)));
  1641.                 l = n - absS;
  1642.             } else {
  1643.                 factCoef = zero.add(1.);
  1644.                 l = n - absJ;
  1645.             }

  1647.             // (-1)<sup>|j-s|</sup>
  1648.             final T sign = absJmS % 2 != 0 ? zero.add(-1.) : zero.add(1.);
  1649.             //(1 - c²)<sup>n-|s|</sup> / (1 + c²)<sup>n</sup>
  1650.             final T coef1 = omc2tn[l].divide(opc2tn[n]);
  1651.             //-b<sup>|j-s|</sup>
  1652.             final T coef2 = btjms[absJmS].multiply(sign);
  1653.             // P<sub>l</sub><sup>|j-s|, |j+s|</sup>(χ)
  1654.             final FieldGradient<T> jac =
  1655.                     JacobiPolynomials.getValue(l, absJmS, absJpS, FieldGradient.variable(1, 0, context.getX()));

  1657.             // the derivative of coef1 by c
  1658.             final T dcoef1dc = coef1.negate().multiply(2.).multiply(c).multiply(opc2tn[1].reciprocal().multiply(n).add(omc2tn[1].reciprocal().multiply(l)));
  1659.             // the derivative of coef1 by h
  1660.             final T dcoef1dh = dcoef1dc.multiply(dcdh);
  1661.             // the derivative of coef1 by k
  1662.             final T dcoef1dk = dcoef1dc.multiply(dcdk);

  1664.             // the derivative of coef2 by b
  1665.             final T dcoef2db = absJmS == 0 ? zero : sign.multiply(absJmS).multiply(btjms[absJmS - 1]);
  1666.             // the derivative of coef2 by h
  1667.             final T dcoef2dh = dcoef2db.multiply(dbdh);
  1668.             // the derivative of coef2 by k
  1669.             final T dcoef2dk = dcoef2db.multiply(dbdk);

  1671.             // the jacobi polynomial value
  1672.             final T jacobi = jac.getValue();
  1673.             // the derivative of the Jacobi polynomial by h
  1674.             final T djacobidh = jac.getGradient()[0].multiply(context.getHXXX());
  1675.             // the derivative of the Jacobi polynomial by k
  1676.             final T djacobidk = jac.getGradient()[0].multiply(context.getKXXX());

  1678.             //group the above coefficients to limit the mathematical operations
  1679.             final T term1 = factCoef.multiply(coef1).multiply(coef2);
  1680.             final T term2 = factCoef.multiply(coef1).multiply(jacobi);
  1681.             final T term3 = factCoef.multiply(coef2).multiply(jacobi);

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

  1688.             return wjnsemjms;
  1689.         }
  1690.     }

  1692.     /** The G<sub>n,s</sub> coefficients and their derivatives.
  1693.      * <p>
  1694.      * See Danielson, 4.2-17
  1695.      *
  1696.      * @author Lucian Barbulescu
  1697.      */
  1698.     private class GnsCoefficients {

  1700.         /** Maximum value for n index. */
  1701.         private final int nMax;

  1703.         /** Maximum value for s index. */
  1704.         private final int sMax;

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

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

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

  1715.         /** Standard constructor.
  1716.          *
  1717.          * @param nMax maximum value for n indes
  1718.          * @param sMax maximum value for s index
  1719.          * @param context container for attributes
  1720.          */
  1721.         GnsCoefficients(final int nMax, final int sMax, final DSSTThirdBodyContext context) {
  1722.             this.nMax = nMax;
  1723.             this.sMax = sMax;

  1725.             final int rows    = nMax + 1;
  1726.             final int columns = sMax + 1;
  1727.             this.gns          = new double[rows][columns];
  1728.             this.dgnsda       = new double[rows][columns];
  1729.             this.dgnsdgamma   = new double[rows][columns];

  1731.             // Generate the coefficients
  1732.             generateCoefficients(context);
  1733.         }
  1734.         /**
  1735.          * Compute the coefficient G<sub>n,s</sub> and its derivatives.
  1736.          * <p>
  1737.          * Only the derivatives by a and γ are computed as all others are 0
  1738.          * </p>
  1739.          * @param context container for attributes
  1740.          */
  1741.         private void generateCoefficients(final DSSTThirdBodyContext context) {

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

  1745.             for (int s = 0; s <= sMax; s++) {
  1746.                 // The n index is always at least the maximum between 2 and s
  1747.                 final int minN = FastMath.max(2, s);
  1748.                 for (int n = minN; n <= nMax; n++) {
  1749.                     // compute the coefficients only if (n - s) % 2 == 0
  1750.                     if ( (n - s) % 2 == 0 ) {
  1751.                         // Kronecker symbol (2 - delta(0,s))
  1752.                         final double delta0s = (s == 0) ? 1. : 2.;
  1753.                         final double vns   = Vns.get(new NSKey(n, s));
  1754.                         final double coef0 = delta0s * context.getAoR3Pow()[n] * vns * context.getMuoR3();
  1755.                         final double coef1 = coef0 * context.getQns()[n][s];
  1756.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  1757.                         // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
  1758.                         final double dqns = (n == s) ? 0. : context.getQns()[n][s + 1];

  1760.                         //Compute the coefficient and its derivatives.
  1761.                         this.gns[n][s] = coef1;
  1762.                         this.dgnsda[n][s] = coef1 * n / auxiliaryElements.getSma();
  1763.                         this.dgnsdgamma[n][s] = coef0 * dqns;
  1764.                     } else {
  1765.                         // the coefficient and its derivatives is 0
  1766.                         this.gns[n][s] = 0.;
  1767.                         this.dgnsda[n][s] = 0.;
  1768.                         this.dgnsdgamma[n][s] = 0.;
  1769.                     }
  1770.                 }
  1771.             }
  1772.         }

  1774.         /** Get the coefficient G<sub>n,s</sub>.
  1775.          *
  1776.          * @param n n index
  1777.          * @param s s index
  1778.          * @return the coefficient G<sub>n,s</sub>
  1779.          */
  1780.         public double getGns(final int n, final int s) {
  1781.             return this.gns[n][s];
  1782.         }

  1784.         /** Get the derivative dG<sub>n,s</sub> / da.
  1785.          *
  1786.          * @param n n index
  1787.          * @param s s index
  1788.          * @return the derivative dG<sub>n,s</sub> / da
  1789.          */
  1790.         public double getdGnsda(final int n, final int s) {
  1791.             return this.dgnsda[n][s];
  1792.         }

  1794.         /** Get the derivative dG<sub>n,s</sub> / dγ.
  1795.          *
  1796.          * @param n n index
  1797.          * @param s s index
  1798.          * @return the derivative dG<sub>n,s</sub> / dγ
  1799.          */
  1800.         public double getdGnsdgamma(final int n, final int s) {
  1801.             return this.dgnsdgamma[n][s];
  1802.         }
  1803.     }

  1805.     /** The G<sub>n,s</sub> coefficients and their derivatives.
  1806.      * <p>
  1807.      * See Danielson, 4.2-17
  1808.      *
  1809.      * @author Lucian Barbulescu
  1810.      */
  1811.     private class FieldGnsCoefficients  <T extends RealFieldElement<T>> {

  1813.         /** Maximum value for n index. */
  1814.         private final int nMax;

  1816.         /** Maximum value for s index. */
  1817.         private final int sMax;

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

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

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

  1828.         /** Standard constructor.
  1829.          *
  1830.          * @param nMax maximum value for n indes
  1831.          * @param sMax maximum value for s index
  1832.          * @param context container for attributes
  1833.          * @param field field used by default
  1834.          */
  1835.         FieldGnsCoefficients(final int nMax, final int sMax,
  1836.                              final FieldDSSTThirdBodyContext<T> context,
  1837.                              final Field<T> field) {
  1838.             this.nMax = nMax;
  1839.             this.sMax = sMax;

  1841.             final int rows    = nMax + 1;
  1842.             final int columns = sMax + 1;
  1843.             this.gns          = MathArrays.buildArray(field, rows, columns);
  1844.             this.dgnsda       = MathArrays.buildArray(field, rows, columns);
  1845.             this.dgnsdgamma   = MathArrays.buildArray(field, rows, columns);

  1847.             // Generate the coefficients
  1848.             generateCoefficients(context, field);
  1849.         }
  1850.         /**
  1851.          * Compute the coefficient G<sub>n,s</sub> and its derivatives.
  1852.          * <p>
  1853.          * Only the derivatives by a and γ are computed as all others are 0
  1854.          * </p>
  1855.          * @param context container for attributes
  1856.          * @param field field used by default
  1857.          */
  1858.         private void generateCoefficients(final FieldDSSTThirdBodyContext<T> context,
  1859.                                           final Field<T> field) {

  1861.             //Zero
  1862.             final T zero = field.getZero();

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

  1866.             for (int s = 0; s <= sMax; s++) {
  1867.                 // The n index is always at least the maximum between 2 and s
  1868.                 final int minN = FastMath.max(2, s);
  1869.                 for (int n = minN; n <= nMax; n++) {
  1870.                     // compute the coefficients only if (n - s) % 2 == 0
  1871.                     if ( (n - s) % 2 == 0 ) {
  1872.                         // Kronecker symbol (2 - delta(0,s))
  1873.                         final T delta0s = (s == 0) ? zero.add(1.) : zero.add(2.);
  1874.                         final double vns = Vns.get(new NSKey(n, s));
  1875.                         final T coef0 = context.getAoR3Pow()[n].multiply(vns).multiply(context.getMuoR3()).multiply(delta0s);
  1876.                         final T coef1 = coef0.multiply(context.getQns()[n][s]);
  1877.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  1878.                         // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
  1879.                         final T dqns = (n == s) ? zero : context.getQns()[n][s + 1];

  1881.                         //Compute the coefficient and its derivatives.
  1882.                         this.gns[n][s] = coef1;
  1883.                         this.dgnsda[n][s] = coef1.multiply(n).divide(auxiliaryElements.getSma());
  1884.                         this.dgnsdgamma[n][s] = coef0.multiply(dqns);
  1885.                     } else {
  1886.                         // the coefficient and its derivatives is 0
  1887.                         this.gns[n][s] = zero;
  1888.                         this.dgnsda[n][s] = zero;
  1889.                         this.dgnsdgamma[n][s] = zero;
  1890.                     }
  1891.                 }
  1892.             }
  1893.         }

  1895.         /** Get the coefficient G<sub>n,s</sub>.
  1896.          *
  1897.          * @param n n index
  1898.          * @param s s index
  1899.          * @return the coefficient G<sub>n,s</sub>
  1900.          */
  1901.         public T getGns(final int n, final int s) {
  1902.             return this.gns[n][s];
  1903.         }

  1905.         /** Get the derivative dG<sub>n,s</sub> / da.
  1906.          *
  1907.          * @param n n index
  1908.          * @param s s index
  1909.          * @return the derivative dG<sub>n,s</sub> / da
  1910.          */
  1911.         public T getdGnsda(final int n, final int s) {
  1912.             return this.dgnsda[n][s];
  1913.         }

  1915.         /** Get the derivative dG<sub>n,s</sub> / dγ.
  1916.          *
  1917.          * @param n n index
  1918.          * @param s s index
  1919.          * @return the derivative dG<sub>n,s</sub> / dγ
  1920.          */
  1921.         public T getdGnsdgamma(final int n, final int s) {
  1922.             return this.dgnsdgamma[n][s];
  1923.         }
  1924.     }

  1926.     /** This class computes the terms containing the coefficients C<sub>j</sub> and S<sub>j</sub> of (α, β) or (k, h).
  1927.      *
  1928.      * <p>
  1929.      * The following terms and their derivatives by k, h, alpha and beta are considered: <br/ >
  1930.      * - 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 />
  1931.      * - C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) <br />
  1932.      * - 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 />
  1933.      * - C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) <br />
  1934.      * 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 />
  1935.      * See the CS Mathematical report $3.5.3.2 for more details
  1936.      * </p>
  1937.      * @author Lucian Barbulescu
  1938.      */
  1939.     private static class CjSjAlphaBetaKH {

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

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

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

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

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

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

  1963.         /**
  1964.          * Standard constructor.
  1965.          * @param context container for attributes
  1966.          */
  1967.         CjSjAlphaBetaKH(final DSSTThirdBodyContext context) {

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

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

  1974.             coefAandDeriv = new double[5];
  1975.             coefBandDeriv = new double[5];
  1976.             coefDandDeriv = new double[5];
  1977.             coefEandDeriv = new double[5];
  1978.         }

  1980.         /** Compute the coefficients and their derivatives for a given (j,s) pair.
  1981.          * @param j j index
  1982.          * @param s s index
  1983.          */
  1984.         public void computeCoefficients(final int j, final int s) {
  1985.             // sign of j-s
  1986.             final int sign = j < s ? -1 : 1;

  1988.             //|j-s|
  1989.             final int absJmS = FastMath.abs(j - s);

  1991.             //j+s
  1992.             final int jps = j + s;

  1994.             //Compute the coefficient A and its derivatives
  1995.             coefAandDeriv[0] = sign * cjsjalbe.getCj(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getCj(absJmS);
  1996.             coefAandDeriv[1] = sign * cjsjalbe.getCj(s) * cjsjkh.getDsjDk(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getDcjDk(absJmS);
  1997.             coefAandDeriv[2] = sign * cjsjalbe.getCj(s) * cjsjkh.getDsjDh(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getDcjDh(absJmS);
  1998.             coefAandDeriv[3] = sign * cjsjalbe.getDcjDk(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getDsjDk(s) * cjsjkh.getCj(absJmS);
  1999.             coefAandDeriv[4] = sign * cjsjalbe.getDcjDh(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getDsjDh(s) * cjsjkh.getCj(absJmS);

  2001.             //Compute the coefficient B and its derivatives
  2002.             coefBandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getSj(jps) - cjsjalbe.getSj(s) * cjsjkh.getCj(jps);
  2003.             coefBandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDsjDk(jps) - cjsjalbe.getSj(s) * cjsjkh.getDcjDk(jps);
  2004.             coefBandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDsjDh(jps) - cjsjalbe.getSj(s) * cjsjkh.getDcjDh(jps);
  2005.             coefBandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getSj(jps) - cjsjalbe.getDsjDk(s) * cjsjkh.getCj(jps);
  2006.             coefBandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getSj(jps) - cjsjalbe.getDsjDh(s) * cjsjkh.getCj(jps);

  2008.             //Compute the coefficient D and its derivatives
  2009.             coefDandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getSj(absJmS);
  2010.             coefDandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDcjDk(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getDsjDk(absJmS);
  2011.             coefDandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDcjDh(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getDsjDh(absJmS);
  2012.             coefDandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getDsjDk(s) * cjsjkh.getSj(absJmS);
  2013.             coefDandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getDsjDh(s) * cjsjkh.getSj(absJmS);

  2015.             //Compute the coefficient E and its derivatives
  2016.             coefEandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getCj(jps) + cjsjalbe.getSj(s) * cjsjkh.getSj(jps);
  2017.             coefEandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDcjDk(jps) + cjsjalbe.getSj(s) * cjsjkh.getDsjDk(jps);
  2018.             coefEandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDcjDh(jps) + cjsjalbe.getSj(s) * cjsjkh.getDsjDh(jps);
  2019.             coefEandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getCj(jps) + cjsjalbe.getDsjDk(s) * cjsjkh.getSj(jps);
  2020.             coefEandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getCj(jps) + cjsjalbe.getDsjDh(s) * cjsjkh.getSj(jps);
  2021.         }

  2023.         /** Get the value of coefficient A<sub>j,s</sub>.
  2024.          *
  2025.          * @return the coefficient A<sub>j,s</sub>
  2026.          */
  2027.         public double getCoefA() {
  2028.             return coefAandDeriv[0];
  2029.         }

  2031.         /** Get the value of coefficient dA<sub>j,s</sub>/dk.
  2032.          *
  2033.          * @return the coefficient dA<sub>j,s</sub>/dk
  2034.          */
  2035.         public double getdCoefAdk() {
  2036.             return coefAandDeriv[1];
  2037.         }

  2039.         /** Get the value of coefficient dA<sub>j,s</sub>/dh.
  2040.          *
  2041.          * @return the coefficient dA<sub>j,s</sub>/dh
  2042.          */
  2043.         public double getdCoefAdh() {
  2044.             return coefAandDeriv[2];
  2045.         }

  2047.         /** Get the value of coefficient dA<sub>j,s</sub>/dα.
  2048.          *
  2049.          * @return the coefficient dA<sub>j,s</sub>/dα
  2050.          */
  2051.         public double getdCoefAdalpha() {
  2052.             return coefAandDeriv[3];
  2053.         }

  2055.         /** Get the value of coefficient dA<sub>j,s</sub>/dβ.
  2056.          *
  2057.          * @return the coefficient dA<sub>j,s</sub>/dβ
  2058.          */
  2059.         public double getdCoefAdbeta() {
  2060.             return coefAandDeriv[4];
  2061.         }

  2063.         /** Get the value of coefficient B<sub>j,s</sub>.
  2064.          *
  2065.          * @return the coefficient B<sub>j,s</sub>
  2066.          */
  2067.         public double getCoefB() {
  2068.             return coefBandDeriv[0];
  2069.         }

  2071.         /** Get the value of coefficient dB<sub>j,s</sub>/dk.
  2072.          *
  2073.          * @return the coefficient dB<sub>j,s</sub>/dk
  2074.          */
  2075.         public double getdCoefBdk() {
  2076.             return coefBandDeriv[1];
  2077.         }

  2079.         /** Get the value of coefficient dB<sub>j,s</sub>/dh.
  2080.          *
  2081.          * @return the coefficient dB<sub>j,s</sub>/dh
  2082.          */
  2083.         public double getdCoefBdh() {
  2084.             return coefBandDeriv[2];
  2085.         }

  2087.         /** Get the value of coefficient dB<sub>j,s</sub>/dα.
  2088.          *
  2089.          * @return the coefficient dB<sub>j,s</sub>/dα
  2090.          */
  2091.         public double getdCoefBdalpha() {
  2092.             return coefBandDeriv[3];
  2093.         }

  2095.         /** Get the value of coefficient dB<sub>j,s</sub>/dβ.
  2096.          *
  2097.          * @return the coefficient dB<sub>j,s</sub>/dβ
  2098.          */
  2099.         public double getdCoefBdbeta() {
  2100.             return coefBandDeriv[4];
  2101.         }

  2103.         /** Get the value of coefficient D<sub>j,s</sub>.
  2104.          *
  2105.          * @return the coefficient D<sub>j,s</sub>
  2106.          */
  2107.         public double getCoefD() {
  2108.             return coefDandDeriv[0];
  2109.         }

  2111.         /** Get the value of coefficient dD<sub>j,s</sub>/dk.
  2112.          *
  2113.          * @return the coefficient dD<sub>j,s</sub>/dk
  2114.          */
  2115.         public double getdCoefDdk() {
  2116.             return coefDandDeriv[1];
  2117.         }

  2119.         /** Get the value of coefficient dD<sub>j,s</sub>/dh.
  2120.          *
  2121.          * @return the coefficient dD<sub>j,s</sub>/dh
  2122.          */
  2123.         public double getdCoefDdh() {
  2124.             return coefDandDeriv[2];
  2125.         }

  2127.         /** Get the value of coefficient dD<sub>j,s</sub>/dα.
  2128.          *
  2129.          * @return the coefficient dD<sub>j,s</sub>/dα
  2130.          */
  2131.         public double getdCoefDdalpha() {
  2132.             return coefDandDeriv[3];
  2133.         }

  2135.         /** Get the value of coefficient dD<sub>j,s</sub>/dβ.
  2136.          *
  2137.          * @return the coefficient dD<sub>j,s</sub>/dβ
  2138.          */
  2139.         public double getdCoefDdbeta() {
  2140.             return coefDandDeriv[4];
  2141.         }

  2143.         /** Get the value of coefficient E<sub>j,s</sub>.
  2144.          *
  2145.          * @return the coefficient E<sub>j,s</sub>
  2146.          */
  2147.         public double getCoefE() {
  2148.             return coefEandDeriv[0];
  2149.         }

  2151.         /** Get the value of coefficient dE<sub>j,s</sub>/dk.
  2152.          *
  2153.          * @return the coefficient dE<sub>j,s</sub>/dk
  2154.          */
  2155.         public double getdCoefEdk() {
  2156.             return coefEandDeriv[1];
  2157.         }

  2159.         /** Get the value of coefficient dE<sub>j,s</sub>/dh.
  2160.          *
  2161.          * @return the coefficient dE<sub>j,s</sub>/dh
  2162.          */
  2163.         public double getdCoefEdh() {
  2164.             return coefEandDeriv[2];
  2165.         }

  2167.         /** Get the value of coefficient dE<sub>j,s</sub>/dα.
  2168.          *
  2169.          * @return the coefficient dE<sub>j,s</sub>/dα
  2170.          */
  2171.         public double getdCoefEdalpha() {
  2172.             return coefEandDeriv[3];
  2173.         }

  2175.         /** Get the value of coefficient dE<sub>j,s</sub>/dβ.
  2176.          *
  2177.          * @return the coefficient dE<sub>j,s</sub>/dβ
  2178.          */
  2179.         public double getdCoefEdbeta() {
  2180.             return coefEandDeriv[4];
  2181.         }
  2182.     }

  2184.      /** This class computes the terms containing the coefficients C<sub>j</sub> and S<sub>j</sub> of (α, β) or (k, h).
  2185.      *
  2186.      * <p>
  2187.      * The following terms and their derivatives by k, h, alpha and beta are considered: <br/ >
  2188.      * - 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 />
  2189.      * - C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) <br />
  2190.      * - 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 />
  2191.      * - C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) <br />
  2192.      * 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 />
  2193.      * See the CS Mathematical report $3.5.3.2 for more details
  2194.      * </p>
  2195.      * @author Lucian Barbulescu
  2196.      */
  2197.     private static class FieldCjSjAlphaBetaKH <T extends RealFieldElement<T>> {

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

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

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

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

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

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

  2221.         /**
  2222.          * Standard constructor.
  2223.          * @param context container for attributes
  2224.          * @param field field used by default
  2225.          */
  2226.         FieldCjSjAlphaBetaKH(final FieldDSSTThirdBodyContext<T> context, final Field<T> field) {

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

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

  2233.             coefAandDeriv = MathArrays.buildArray(field, 5);
  2234.             coefBandDeriv = MathArrays.buildArray(field, 5);
  2235.             coefDandDeriv = MathArrays.buildArray(field, 5);
  2236.             coefEandDeriv = MathArrays.buildArray(field, 5);
  2237.         }

  2239.         /** Compute the coefficients and their derivatives for a given (j,s) pair.
  2240.          * @param j j index
  2241.          * @param s s index
  2242.          */
  2243.         public void computeCoefficients(final int j, final int s) {
  2244.             // sign of j-s
  2245.             final int sign = j < s ? -1 : 1;

  2247.             //|j-s|
  2248.             final int absJmS = FastMath.abs(j - s);

  2250.             //j+s
  2251.             final int jps = j + s;

  2253.             //Compute the coefficient A and its derivatives
  2254.             coefAandDeriv[0] = cjsjalbe.getCj(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign).add(cjsjalbe.getSj(s).multiply(cjsjkh.getCj(absJmS)));
  2255.             coefAandDeriv[1] = cjsjalbe.getCj(s).multiply(cjsjkh.getDsjDk(absJmS)).multiply(sign).add(cjsjalbe.getSj(s).multiply(cjsjkh.getDcjDk(absJmS)));
  2256.             coefAandDeriv[2] = cjsjalbe.getCj(s).multiply(cjsjkh.getDsjDh(absJmS)).multiply(sign).add(cjsjalbe.getSj(s).multiply(cjsjkh.getDcjDh(absJmS)));
  2257.             coefAandDeriv[3] = cjsjalbe.getDcjDk(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign).add(cjsjalbe.getDsjDk(s).multiply(cjsjkh.getCj(absJmS)));
  2258.             coefAandDeriv[4] = cjsjalbe.getDcjDh(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign).add(cjsjalbe.getDsjDh(s).multiply(cjsjkh.getCj(absJmS)));

  2260.             //Compute the coefficient B and its derivatives
  2261.             coefBandDeriv[0] = cjsjalbe.getCj(s).multiply(cjsjkh.getSj(jps)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getCj(jps)));
  2262.             coefBandDeriv[1] = cjsjalbe.getCj(s).multiply(cjsjkh.getDsjDk(jps)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getDcjDk(jps)));
  2263.             coefBandDeriv[2] = cjsjalbe.getCj(s).multiply(cjsjkh.getDsjDh(jps)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getDcjDh(jps)));
  2264.             coefBandDeriv[3] = cjsjalbe.getDcjDk(s).multiply(cjsjkh.getSj(jps)).subtract(cjsjalbe.getDsjDk(s).multiply(cjsjkh.getCj(jps)));
  2265.             coefBandDeriv[4] = cjsjalbe.getDcjDh(s).multiply(cjsjkh.getSj(jps)).subtract(cjsjalbe.getDsjDh(s).multiply(cjsjkh.getCj(jps)));

  2267.             //Compute the coefficient D and its derivatives
  2268.             coefDandDeriv[0] = cjsjalbe.getCj(s).multiply(cjsjkh.getCj(absJmS)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign));
  2269.             coefDandDeriv[1] = cjsjalbe.getCj(s).multiply(cjsjkh.getDcjDk(absJmS)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getDsjDk(absJmS)).multiply(sign));
  2270.             coefDandDeriv[2] = cjsjalbe.getCj(s).multiply(cjsjkh.getDcjDh(absJmS)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getDsjDh(absJmS)).multiply(sign));
  2271.             coefDandDeriv[3] = cjsjalbe.getDcjDk(s).multiply(cjsjkh.getCj(absJmS)).subtract(cjsjalbe.getDsjDk(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign));
  2272.             coefDandDeriv[4] = cjsjalbe.getDcjDh(s).multiply(cjsjkh.getCj(absJmS)).subtract(cjsjalbe.getDsjDh(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign));

  2274.             //Compute the coefficient E and its derivatives
  2275.             coefEandDeriv[0] = cjsjalbe.getCj(s).multiply(cjsjkh.getCj(jps)).add(cjsjalbe.getSj(s).multiply(cjsjkh.getSj(jps)));
  2276.             coefEandDeriv[1] = cjsjalbe.getCj(s).multiply(cjsjkh.getDcjDk(jps)).add(cjsjalbe.getSj(s).multiply(cjsjkh.getDsjDk(jps)));
  2277.             coefEandDeriv[2] = cjsjalbe.getCj(s).multiply(cjsjkh.getDcjDh(jps)).add(cjsjalbe.getSj(s).multiply(cjsjkh.getDsjDh(jps)));
  2278.             coefEandDeriv[3] = cjsjalbe.getDcjDk(s).multiply(cjsjkh.getCj(jps)).add(cjsjalbe.getDsjDk(s).multiply(cjsjkh.getSj(jps)));
  2279.             coefEandDeriv[4] = cjsjalbe.getDcjDh(s).multiply(cjsjkh.getCj(jps)).add(cjsjalbe.getDsjDh(s).multiply(cjsjkh.getSj(jps)));
  2280.         }

  2282.         /** Get the value of coefficient A<sub>j,s</sub>.
  2283.          *
  2284.          * @return the coefficient A<sub>j,s</sub>
  2285.          */
  2286.         public T getCoefA() {
  2287.             return coefAandDeriv[0];
  2288.         }

  2290.         /** Get the value of coefficient dA<sub>j,s</sub>/dk.
  2291.          *
  2292.          * @return the coefficient dA<sub>j,s</sub>/dk
  2293.          */
  2294.         public T getdCoefAdk() {
  2295.             return coefAandDeriv[1];
  2296.         }

  2298.         /** Get the value of coefficient dA<sub>j,s</sub>/dh.
  2299.          *
  2300.          * @return the coefficient dA<sub>j,s</sub>/dh
  2301.          */
  2302.         public T getdCoefAdh() {
  2303.             return coefAandDeriv[2];
  2304.         }

  2306.         /** Get the value of coefficient dA<sub>j,s</sub>/dα.
  2307.          *
  2308.          * @return the coefficient dA<sub>j,s</sub>/dα
  2309.          */
  2310.         public T getdCoefAdalpha() {
  2311.             return coefAandDeriv[3];
  2312.         }

  2314.         /** Get the value of coefficient dA<sub>j,s</sub>/dβ.
  2315.          *
  2316.          * @return the coefficient dA<sub>j,s</sub>/dβ
  2317.          */
  2318.         public T getdCoefAdbeta() {
  2319.             return coefAandDeriv[4];
  2320.         }

  2322.        /** Get the value of coefficient B<sub>j,s</sub>.
  2323.         *
  2324.         * @return the coefficient B<sub>j,s</sub>
  2325.         */
  2326.         public T getCoefB() {
  2327.             return coefBandDeriv[0];
  2328.         }

  2330.         /** Get the value of coefficient dB<sub>j,s</sub>/dk.
  2331.          *
  2332.          * @return the coefficient dB<sub>j,s</sub>/dk
  2333.          */
  2334.         public T getdCoefBdk() {
  2335.             return coefBandDeriv[1];
  2336.         }

  2338.         /** Get the value of coefficient dB<sub>j,s</sub>/dh.
  2339.          *
  2340.          * @return the coefficient dB<sub>j,s</sub>/dh
  2341.          */
  2342.         public T getdCoefBdh() {
  2343.             return coefBandDeriv[2];
  2344.         }

  2346.         /** Get the value of coefficient dB<sub>j,s</sub>/dα.
  2347.          *
  2348.          * @return the coefficient dB<sub>j,s</sub>/dα
  2349.          */
  2350.         public T getdCoefBdalpha() {
  2351.             return coefBandDeriv[3];
  2352.         }

  2354.         /** Get the value of coefficient dB<sub>j,s</sub>/dβ.
  2355.          *
  2356.          * @return the coefficient dB<sub>j,s</sub>/dβ
  2357.          */
  2358.         public T getdCoefBdbeta() {
  2359.             return coefBandDeriv[4];
  2360.         }

  2362.         /** Get the value of coefficient D<sub>j,s</sub>.
  2363.          *
  2364.          * @return the coefficient D<sub>j,s</sub>
  2365.          */
  2366.         public T getCoefD() {
  2367.             return coefDandDeriv[0];
  2368.         }

  2370.         /** Get the value of coefficient dD<sub>j,s</sub>/dk.
  2371.          *
  2372.          * @return the coefficient dD<sub>j,s</sub>/dk
  2373.          */
  2374.         public T getdCoefDdk() {
  2375.             return coefDandDeriv[1];
  2376.         }

  2378.         /** Get the value of coefficient dD<sub>j,s</sub>/dh.
  2379.          *
  2380.          * @return the coefficient dD<sub>j,s</sub>/dh
  2381.          */
  2382.         public T getdCoefDdh() {
  2383.             return coefDandDeriv[2];
  2384.         }

  2386.         /** Get the value of coefficient dD<sub>j,s</sub>/dα.
  2387.          *
  2388.          * @return the coefficient dD<sub>j,s</sub>/dα
  2389.          */
  2390.         public T getdCoefDdalpha() {
  2391.             return coefDandDeriv[3];
  2392.         }

  2394.         /** Get the value of coefficient dD<sub>j,s</sub>/dβ.
  2395.          *
  2396.          * @return the coefficient dD<sub>j,s</sub>/dβ
  2397.          */
  2398.         public T getdCoefDdbeta() {
  2399.             return coefDandDeriv[4];
  2400.         }

  2402.         /** Get the value of coefficient E<sub>j,s</sub>.
  2403.          *
  2404.          * @return the coefficient E<sub>j,s</sub>
  2405.          */
  2406.         public T getCoefE() {
  2407.             return coefEandDeriv[0];
  2408.         }

  2410.         /** Get the value of coefficient dE<sub>j,s</sub>/dk.
  2411.          *
  2412.          * @return the coefficient dE<sub>j,s</sub>/dk
  2413.          */
  2414.         public T getdCoefEdk() {
  2415.             return coefEandDeriv[1];
  2416.         }

  2418.         /** Get the value of coefficient dE<sub>j,s</sub>/dh.
  2419.          *
  2420.          * @return the coefficient dE<sub>j,s</sub>/dh
  2421.          */
  2422.         public T getdCoefEdh() {
  2423.             return coefEandDeriv[2];
  2424.         }

  2426.         /** Get the value of coefficient dE<sub>j,s</sub>/dα.
  2427.          *
  2428.          * @return the coefficient dE<sub>j,s</sub>/dα
  2429.          */
  2430.         public T getdCoefEdalpha() {
  2431.             return coefEandDeriv[3];
  2432.         }

  2434.         /** Get the value of coefficient dE<sub>j,s</sub>/dβ.
  2435.          *
  2436.          * @return the coefficient dE<sub>j,s</sub>/dβ
  2437.          */
  2438.         public T getdCoefEdbeta() {
  2439.             return coefEandDeriv[4];
  2440.         }
  2441.     }

  2443.     /** This class computes the coefficients for the generating function S and its derivatives.
  2444.      * <p>
  2445.      * The form of the generating functions is: <br>
  2446.      *  S = C⁰ + &Sigma;<sub>j=1</sub><sup>N+1</sup>(C<sup>j</sup> * cos(jF) + S<sup>j</sup> * sin(jF)) <br>
  2447.      *  The coefficients C⁰, C<sup>j</sup>, S<sup>j</sup> are the Fourrier coefficients
  2448.      *  presented in Danielson 4.2-14,15 except for the case j=1 where
  2449.      *  C¹ = C¹<sub>Fourier</sub> - hU and
  2450.      *  S¹ = S¹<sub>Fourier</sub> + kU <br>
  2451.      *  Also the coefficients of the derivatives of S by a, k, h, α, β, γ and λ
  2452.      *  are computed end expressed in a similar manner. The formulas used are 4.2-19, 20, 23, 24
  2453.      * </p>
  2454.      * @author Lucian Barbulescu
  2455.      */
  2456.     private class GeneratingFunctionCoefficients {

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

  2461.         /** Maximum value of j index. */
  2462.         private final int jMax;

  2464.         /** The coefficients C<sup>j</sup> of the function S and its derivatives.
  2465.          * <p>
  2466.          * The index j belongs to the interval [0,jMax]. The coefficient C⁰ is the free coefficient.<br>
  2467.          * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
  2468.          * - S <br/>
  2469.          * - dS / da <br/>
  2470.          * - dS / dk <br/>
  2471.          * - dS / dh <br/>
  2472.          * - dS / dα <br/>
  2473.          * - dS / dβ <br/>
  2474.          * - dS / dγ <br/>
  2475.          * - dS / dλ
  2476.          * </p>
  2477.          */
  2478.         private final double[][] cjCoefs;

  2480.         /** The coefficients S<sup>j</sup> of the function S and its derivatives.
  2481.          * <p>
  2482.          * The index j belongs to the interval [0,jMax].<br>
  2483.          * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
  2484.          * - S <br/>
  2485.          * - dS / da <br/>
  2486.          * - dS / dk <br/>
  2487.          * - dS / dh <br/>
  2488.          * - dS / dα <br/>
  2489.          * - dS / dβ <br/>
  2490.          * - dS / dγ <br/>
  2491.          * - dS / dλ
  2492.          * </p>
  2493.          */
  2494.         private final double[][] sjCoefs;

  2496.         /**
  2497.          * Standard constructor.
  2498.          *
  2499.          * @param nMax maximum value of n index
  2500.          * @param sMax maximum value of s index
  2501.          * @param jMax maximum value of j index
  2502.          * @param context container for attributes
  2503.          * @param hansen hansen objects
  2504.          */
  2505.         GeneratingFunctionCoefficients(final int nMax, final int sMax, final int jMax,
  2506.                                        final DSSTThirdBodyContext context, final HansenObjects hansen) {
  2507.             this.jMax = jMax;
  2508.             this.cjsjFourier = new FourierCjSjCoefficients(nMax, sMax, jMax, context);
  2509.             this.cjCoefs = new double[8][jMax + 1];
  2510.             this.sjCoefs = new double[8][jMax + 1];

  2512.             computeGeneratingFunctionCoefficients(context, hansen);
  2513.         }

  2515.         /**
  2516.          * Compute the coefficients for the generating function S and its derivatives.
  2517.          * @param context container for attributes
  2518.          * @param hansenObjects hansen objects
  2519.          */
  2520.         private void computeGeneratingFunctionCoefficients(final DSSTThirdBodyContext context, final HansenObjects hansenObjects) {

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

  2524.             // Access to potential U derivatives
  2525.             final UAnddU udu = new UAnddU(context, hansenObjects);

  2527.             //Compute the C<sup>j</sup> coefficients
  2528.             for (int j = 1; j <= jMax; j++) {
  2529.                 //Compute the C<sup>j</sup> coefficients
  2530.                 cjCoefs[0][j] = cjsjFourier.getCj(j);
  2531.                 cjCoefs[1][j] = cjsjFourier.getdCjda(j);
  2532.                 cjCoefs[2][j] = cjsjFourier.getdCjdk(j) - (cjsjFourier.getSjLambda(j - 1) - cjsjFourier.getSjLambda(j + 1)) / 2;
  2533.                 cjCoefs[3][j] = cjsjFourier.getdCjdh(j) - (cjsjFourier.getCjLambda(j - 1) + cjsjFourier.getCjLambda(j + 1)) / 2;
  2534.                 cjCoefs[4][j] = cjsjFourier.getdCjdalpha(j);
  2535.                 cjCoefs[5][j] = cjsjFourier.getdCjdbeta(j);
  2536.                 cjCoefs[6][j] = cjsjFourier.getdCjdgamma(j);
  2537.                 cjCoefs[7][j] = cjsjFourier.getCjLambda(j);

  2539.                 //Compute the S<sup>j</sup> coefficients
  2540.                 sjCoefs[0][j] = cjsjFourier.getSj(j);
  2541.                 sjCoefs[1][j] = cjsjFourier.getdSjda(j);
  2542.                 sjCoefs[2][j] = cjsjFourier.getdSjdk(j) + (cjsjFourier.getCjLambda(j - 1) - cjsjFourier.getCjLambda(j + 1)) / 2;
  2543.                 sjCoefs[3][j] = cjsjFourier.getdSjdh(j) - (cjsjFourier.getSjLambda(j - 1) + cjsjFourier.getSjLambda(j + 1)) / 2;
  2544.                 sjCoefs[4][j] = cjsjFourier.getdSjdalpha(j);
  2545.                 sjCoefs[5][j] = cjsjFourier.getdSjdbeta(j);
  2546.                 sjCoefs[6][j] = cjsjFourier.getdSjdgamma(j);
  2547.                 sjCoefs[7][j] = cjsjFourier.getSjLambda(j);

  2549.                 //In the special case j == 1 there are some additional terms to be added
  2550.                 if (j == 1) {
  2551.                     //Additional terms for C<sup>j</sup> coefficients
  2552.                     cjCoefs[0][j] += -auxiliaryElements.getH() * udu.getU();
  2553.                     cjCoefs[1][j] += -auxiliaryElements.getH() * udu.getdUda();
  2554.                     cjCoefs[2][j] += -auxiliaryElements.getH() * udu.getdUdk();
  2555.                     cjCoefs[3][j] += -(auxiliaryElements.getH() * udu.getdUdh() + udu.getU() + cjsjFourier.getC0Lambda());
  2556.                     cjCoefs[4][j] += -auxiliaryElements.getH() * udu.getdUdAl();
  2557.                     cjCoefs[5][j] += -auxiliaryElements.getH() * udu.getdUdBe();
  2558.                     cjCoefs[6][j] += -auxiliaryElements.getH() * udu.getdUdGa();

  2560.                     //Additional terms for S<sup>j</sup> coefficients
  2561.                     sjCoefs[0][j] += auxiliaryElements.getK() * udu.getU();
  2562.                     sjCoefs[1][j] += auxiliaryElements.getK() * udu.getdUda();
  2563.                     sjCoefs[2][j] += auxiliaryElements.getK() * udu.getdUdk() + udu.getU() + cjsjFourier.getC0Lambda();
  2564.                     sjCoefs[3][j] += auxiliaryElements.getK() * udu.getdUdh();
  2565.                     sjCoefs[4][j] += auxiliaryElements.getK() * udu.getdUdAl();
  2566.                     sjCoefs[5][j] += auxiliaryElements.getK() * udu.getdUdBe();
  2567.                     sjCoefs[6][j] += auxiliaryElements.getK() * udu.getdUdGa();
  2568.                 }
  2569.             }
  2570.         }

  2572.         /** Get the coefficient C<sup>j</sup> for the function S.
  2573.          * <br>
  2574.          * Possible values for j are within the interval [0,jMax].
  2575.          * The value 0 is used to obtain the free coefficient C⁰
  2576.          * @param j j index
  2577.          * @return C<sup>j</sup> for the function S
  2578.          */
  2579.         public double getSCj(final int j) {
  2580.             return cjCoefs[0][j];
  2581.         }

  2583.         /** Get the coefficient S<sup>j</sup> for the function S.
  2584.          * <br>
  2585.          * Possible values for j are within the interval [1,jMax].
  2586.          * @param j j index
  2587.          * @return S<sup>j</sup> for the function S
  2588.          */
  2589.         public double getSSj(final int j) {
  2590.             return sjCoefs[0][j];
  2591.         }

  2593.         /** Get the coefficient C<sup>j</sup> for the derivative dS/da.
  2594.          * <br>
  2595.          * Possible values for j are within the interval [0,jMax].
  2596.          * The value 0 is used to obtain the free coefficient C⁰
  2597.          * @param j j index
  2598.          * @return C<sup>j</sup> for the function dS/da
  2599.          */
  2600.         public double getdSdaCj(final int j) {
  2601.             return cjCoefs[1][j];
  2602.         }

  2604.         /** Get the coefficient S<sup>j</sup> for the derivative dS/da.
  2605.          * <br>
  2606.          * Possible values for j are within the interval [1,jMax].
  2607.          * @param j j index
  2608.          * @return S<sup>j</sup> for the derivative dS/da
  2609.          */
  2610.         public double getdSdaSj(final int j) {
  2611.             return sjCoefs[1][j];
  2612.         }

  2614.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dk
  2615.          * <br>
  2616.          * Possible values for j are within the interval [0,jMax].
  2617.          * The value 0 is used to obtain the free coefficient C⁰
  2618.          * @param j j index
  2619.          * @return C<sup>j</sup> for the function dS/dk
  2620.          */
  2621.         public double getdSdkCj(final int j) {
  2622.             return cjCoefs[2][j];
  2623.         }

  2625.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dk.
  2626.          * <br>
  2627.          * Possible values for j are within the interval [1,jMax].
  2628.          * @param j j index
  2629.          * @return S<sup>j</sup> for the derivative dS/dk
  2630.          */
  2631.         public double getdSdkSj(final int j) {
  2632.             return sjCoefs[2][j];
  2633.         }

  2635.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dh
  2636.          * <br>
  2637.          * Possible values for j are within the interval [0,jMax].
  2638.          * The value 0 is used to obtain the free coefficient C⁰
  2639.          * @param j j index
  2640.          * @return C<sup>j</sup> for the function dS/dh
  2641.          */
  2642.         public double getdSdhCj(final int j) {
  2643.             return cjCoefs[3][j];
  2644.         }

  2646.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dh.
  2647.          * <br>
  2648.          * Possible values for j are within the interval [1,jMax].
  2649.          * @param j j index
  2650.          * @return S<sup>j</sup> for the derivative dS/dh
  2651.          */
  2652.         public double getdSdhSj(final int j) {
  2653.             return sjCoefs[3][j];
  2654.         }

  2656.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dα
  2657.          * <br>
  2658.          * Possible values for j are within the interval [0,jMax].
  2659.          * The value 0 is used to obtain the free coefficient C⁰
  2660.          * @param j j index
  2661.          * @return C<sup>j</sup> for the function dS/dα
  2662.          */
  2663.         public double getdSdalphaCj(final int j) {
  2664.             return cjCoefs[4][j];
  2665.         }

  2667.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dα.
  2668.          * <br>
  2669.          * Possible values for j are within the interval [1,jMax].
  2670.          * @param j j index
  2671.          * @return S<sup>j</sup> for the derivative dS/dα
  2672.          */
  2673.         public double getdSdalphaSj(final int j) {
  2674.             return sjCoefs[4][j];
  2675.         }

  2677.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dβ
  2678.          * <br>
  2679.          * Possible values for j are within the interval [0,jMax].
  2680.          * The value 0 is used to obtain the free coefficient C⁰
  2681.          * @param j j index
  2682.          * @return C<sup>j</sup> for the function dS/dβ
  2683.          */
  2684.         public double getdSdbetaCj(final int j) {
  2685.             return cjCoefs[5][j];
  2686.         }

  2688.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dβ.
  2689.          * <br>
  2690.          * Possible values for j are within the interval [1,jMax].
  2691.          * @param j j index
  2692.          * @return S<sup>j</sup> for the derivative dS/dβ
  2693.          */
  2694.         public double getdSdbetaSj(final int j) {
  2695.             return sjCoefs[5][j];
  2696.         }

  2698.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dγ
  2699.          * <br>
  2700.          * Possible values for j are within the interval [0,jMax].
  2701.          * The value 0 is used to obtain the free coefficient C⁰
  2702.          * @param j j index
  2703.          * @return C<sup>j</sup> for the function dS/dγ
  2704.          */
  2705.         public double getdSdgammaCj(final int j) {
  2706.             return cjCoefs[6][j];
  2707.         }

  2709.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dγ.
  2710.          * <br>
  2711.          * Possible values for j are within the interval [1,jMax].
  2712.          * @param j j index
  2713.          * @return S<sup>j</sup> for the derivative dS/dγ
  2714.          */
  2715.         public double getdSdgammaSj(final int j) {
  2716.             return sjCoefs[6][j];
  2717.         }

  2719.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dλ
  2720.          * <br>
  2721.          * Possible values for j are within the interval [0,jMax].
  2722.          * The value 0 is used to obtain the free coefficient C⁰
  2723.          * @param j j index
  2724.          * @return C<sup>j</sup> for the function dS/dλ
  2725.          */
  2726.         public double getdSdlambdaCj(final int j) {
  2727.             return cjCoefs[7][j];
  2728.         }

  2730.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dλ.
  2731.          * <br>
  2732.          * Possible values for j are within the interval [1,jMax].
  2733.          * @param j j index
  2734.          * @return S<sup>j</sup> for the derivative dS/dλ
  2735.          */
  2736.         public double getdSdlambdaSj(final int j) {
  2737.             return sjCoefs[7][j];
  2738.         }
  2739.     }

  2741.     /** This class computes the coefficients for the generating function S and its derivatives.
  2742.      * <p>
  2743.      * The form of the generating functions is: <br>
  2744.      *  S = C⁰ + &Sigma;<sub>j=1</sub><sup>N+1</sup>(C<sup>j</sup> * cos(jF) + S<sup>j</sup> * sin(jF)) <br>
  2745.      *  The coefficients C⁰, C<sup>j</sup>, S<sup>j</sup> are the Fourrier coefficients
  2746.      *  presented in Danielson 4.2-14,15 except for the case j=1 where
  2747.      *  C¹ = C¹<sub>Fourier</sub> - hU and
  2748.      *  S¹ = S¹<sub>Fourier</sub> + kU <br>
  2749.      *  Also the coefficients of the derivatives of S by a, k, h, α, β, γ and λ
  2750.      *  are computed end expressed in a similar manner. The formulas used are 4.2-19, 20, 23, 24
  2751.      * </p>
  2752.      * @author Lucian Barbulescu
  2753.      */
  2754.     private class FieldGeneratingFunctionCoefficients <T extends RealFieldElement<T>> {

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

  2759.         /** Maximum value of j index. */
  2760.         private final int jMax;

  2762.         /** The coefficients C<sup>j</sup> of the function S and its derivatives.
  2763.          * <p>
  2764.          * The index j belongs to the interval [0,jMax]. The coefficient C⁰ is the free coefficient.<br>
  2765.          * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
  2766.          * - S <br/>
  2767.          * - dS / da <br/>
  2768.          * - dS / dk <br/>
  2769.          * - dS / dh <br/>
  2770.          * - dS / dα <br/>
  2771.          * - dS / dβ <br/>
  2772.          * - dS / dγ <br/>
  2773.          * - dS / dλ
  2774.          * </p>
  2775.          */
  2776.         private final T[][] cjCoefs;

  2778.         /** The coefficients S<sup>j</sup> of the function S and its derivatives.
  2779.          * <p>
  2780.          * The index j belongs to the interval [0,jMax].<br>
  2781.          * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
  2782.          * - S <br/>
  2783.          * - dS / da <br/>
  2784.          * - dS / dk <br/>
  2785.          * - dS / dh <br/>
  2786.          * - dS / dα <br/>
  2787.          * - dS / dβ <br/>
  2788.          * - dS / dγ <br/>
  2789.          * - dS / dλ
  2790.          * </p>
  2791.          */
  2792.         private final T[][] sjCoefs;

  2794.         /**
  2795.          * Standard constructor.
  2796.          *
  2797.          * @param nMax maximum value of n index
  2798.          * @param sMax maximum value of s index
  2799.          * @param jMax maximum value of j index
  2800.          * @param context container for attributes
  2801.          * @param hansen hansen objects
  2802.          * @param field field used by default
  2803.          */
  2804.         FieldGeneratingFunctionCoefficients(final int nMax, final int sMax, final int jMax,
  2805.                                             final FieldDSSTThirdBodyContext<T> context,
  2806.                                             final FieldHansenObjects<T> hansen,
  2807.                                             final Field<T> field) {
  2808.             this.jMax = jMax;
  2809.             this.cjsjFourier = new FieldFourierCjSjCoefficients<>(nMax, sMax, jMax, context, field);
  2810.             this.cjCoefs     = MathArrays.buildArray(field, 8, jMax + 1);
  2811.             this.sjCoefs     = MathArrays.buildArray(field, 8, jMax + 1);

  2813.             computeGeneratingFunctionCoefficients(context, hansen);
  2814.         }

  2816.         /**
  2817.          * Compute the coefficients for the generating function S and its derivatives.
  2818.          * @param context container for attributes
  2819.          * @param hansenObjects hansen objects
  2820.          */
  2821.         private void computeGeneratingFunctionCoefficients(final FieldDSSTThirdBodyContext<T> context,
  2822.                                                            final FieldHansenObjects<T> hansenObjects) {

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

  2826.             // Access to potential U derivatives
  2827.             final FieldUAnddU<T> udu = new FieldUAnddU<>(context, hansenObjects);

  2829.             //Compute the C<sup>j</sup> coefficients
  2830.             for (int j = 1; j <= jMax; j++) {
  2831.                 //Compute the C<sup>j</sup> coefficients
  2832.                 cjCoefs[0][j] = cjsjFourier.getCj(j);
  2833.                 cjCoefs[1][j] = cjsjFourier.getdCjda(j);
  2834.                 cjCoefs[2][j] = cjsjFourier.getdCjdk(j).subtract((cjsjFourier.getSjLambda(j - 1).subtract(cjsjFourier.getSjLambda(j + 1))).divide(2.));
  2835.                 cjCoefs[3][j] = cjsjFourier.getdCjdh(j).subtract((cjsjFourier.getCjLambda(j - 1).add(cjsjFourier.getCjLambda(j + 1))).divide(2.));
  2836.                 cjCoefs[4][j] = cjsjFourier.getdCjdalpha(j);
  2837.                 cjCoefs[5][j] = cjsjFourier.getdCjdbeta(j);
  2838.                 cjCoefs[6][j] = cjsjFourier.getdCjdgamma(j);
  2839.                 cjCoefs[7][j] = cjsjFourier.getCjLambda(j);

  2841.                 //Compute the S<sup>j</sup> coefficients
  2842.                 sjCoefs[0][j] = cjsjFourier.getSj(j);
  2843.                 sjCoefs[1][j] = cjsjFourier.getdSjda(j);
  2844.                 sjCoefs[2][j] = cjsjFourier.getdSjdk(j).add((cjsjFourier.getCjLambda(j - 1).subtract(cjsjFourier.getCjLambda(j + 1))).divide(2.));
  2845.                 sjCoefs[3][j] = cjsjFourier.getdSjdh(j).subtract((cjsjFourier.getSjLambda(j - 1).add(cjsjFourier.getSjLambda(j + 1))).divide(2.));
  2846.                 sjCoefs[4][j] = cjsjFourier.getdSjdalpha(j);
  2847.                 sjCoefs[5][j] = cjsjFourier.getdSjdbeta(j);
  2848.                 sjCoefs[6][j] = cjsjFourier.getdSjdgamma(j);
  2849.                 sjCoefs[7][j] = cjsjFourier.getSjLambda(j);

  2851.                 //In the special case j == 1 there are some additional terms to be added
  2852.                 if (j == 1) {
  2853.                     //Additional terms for C<sup>j</sup> coefficients
  2854.                     cjCoefs[0][j] = cjCoefs[0][j].add(auxiliaryElements.getH().negate().multiply(udu.getU()));
  2855.                     cjCoefs[1][j] = cjCoefs[1][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUda()));
  2856.                     cjCoefs[2][j] = cjCoefs[2][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUdk()));
  2857.                     cjCoefs[3][j] = cjCoefs[3][j].add(auxiliaryElements.getH().multiply(udu.getdUdh()).add(udu.getU()).add(cjsjFourier.getC0Lambda()).negate());
  2858.                     cjCoefs[4][j] = cjCoefs[4][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUdAl()));
  2859.                     cjCoefs[5][j] = cjCoefs[5][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUdBe()));
  2860.                     cjCoefs[6][j] = cjCoefs[6][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUdGa()));

  2862.                     //Additional terms for S<sup>j</sup> coefficients
  2863.                     sjCoefs[0][j] = sjCoefs[0][j].add(auxiliaryElements.getK().multiply(udu.getU()));
  2864.                     sjCoefs[1][j] = sjCoefs[1][j].add(auxiliaryElements.getK().multiply(udu.getdUda()));
  2865.                     sjCoefs[2][j] = sjCoefs[2][j].add(auxiliaryElements.getK().multiply(udu.getdUdk()).add(udu.getU()).add(cjsjFourier.getC0Lambda()));
  2866.                     sjCoefs[3][j] = sjCoefs[3][j].add(auxiliaryElements.getK().multiply(udu.getdUdh()));
  2867.                     sjCoefs[4][j] = sjCoefs[4][j].add(auxiliaryElements.getK().multiply(udu.getdUdAl()));
  2868.                     sjCoefs[5][j] = sjCoefs[5][j].add(auxiliaryElements.getK().multiply(udu.getdUdBe()));
  2869.                     sjCoefs[6][j] = sjCoefs[6][j].add(auxiliaryElements.getK().multiply(udu.getdUdGa()));
  2870.                 }
  2871.             }
  2872.         }

  2874.         /** Get the coefficient C<sup>j</sup> for the function S.
  2875.          * <br>
  2876.          * Possible values for j are within the interval [0,jMax].
  2877.          * The value 0 is used to obtain the free coefficient C⁰
  2878.          * @param j j index
  2879.          * @return C<sup>j</sup> for the function S
  2880.          */
  2881.         public T getSCj(final int j) {
  2882.             return cjCoefs[0][j];
  2883.         }

  2885.         /** Get the coefficient S<sup>j</sup> for the function S.
  2886.          * <br>
  2887.          * Possible values for j are within the interval [1,jMax].
  2888.          * @param j j index
  2889.          * @return S<sup>j</sup> for the function S
  2890.          */
  2891.         public T getSSj(final int j) {
  2892.             return sjCoefs[0][j];
  2893.         }

  2895.         /** Get the coefficient C<sup>j</sup> for the derivative dS/da.
  2896.          * <br>
  2897.          * Possible values for j are within the interval [0,jMax].
  2898.          * The value 0 is used to obtain the free coefficient C⁰
  2899.          * @param j j index
  2900.          * @return C<sup>j</sup> for the function dS/da
  2901.          */
  2902.         public T getdSdaCj(final int j) {
  2903.             return cjCoefs[1][j];
  2904.         }

  2906.         /** Get the coefficient S<sup>j</sup> for the derivative dS/da.
  2907.          * <br>
  2908.          * Possible values for j are within the interval [1,jMax].
  2909.          * @param j j index
  2910.          * @return S<sup>j</sup> for the derivative dS/da
  2911.          */
  2912.         public T getdSdaSj(final int j) {
  2913.             return sjCoefs[1][j];
  2914.         }

  2916.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dk
  2917.          * <br>
  2918.          * Possible values for j are within the interval [0,jMax].
  2919.          * The value 0 is used to obtain the free coefficient C⁰
  2920.          * @param j j index
  2921.          * @return C<sup>j</sup> for the function dS/dk
  2922.          */
  2923.         public T getdSdkCj(final int j) {
  2924.             return cjCoefs[2][j];
  2925.         }

  2927.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dk.
  2928.          * <br>
  2929.          * Possible values for j are within the interval [1,jMax].
  2930.          * @param j j index
  2931.          * @return S<sup>j</sup> for the derivative dS/dk
  2932.          */
  2933.         public T getdSdkSj(final int j) {
  2934.             return sjCoefs[2][j];
  2935.         }

  2937.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dh
  2938.          * <br>
  2939.          * Possible values for j are within the interval [0,jMax].
  2940.          * The value 0 is used to obtain the free coefficient C⁰
  2941.          * @param j j index
  2942.          * @return C<sup>j</sup> for the function dS/dh
  2943.          */
  2944.         public T getdSdhCj(final int j) {
  2945.             return cjCoefs[3][j];
  2946.         }

  2948.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dh.
  2949.          * <br>
  2950.          * Possible values for j are within the interval [1,jMax].
  2951.          * @param j j index
  2952.          * @return S<sup>j</sup> for the derivative dS/dh
  2953.          */
  2954.         public T getdSdhSj(final int j) {
  2955.             return sjCoefs[3][j];
  2956.         }

  2958.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dα
  2959.          * <br>
  2960.          * Possible values for j are within the interval [0,jMax].
  2961.          * The value 0 is used to obtain the free coefficient C⁰
  2962.          * @param j j index
  2963.          * @return C<sup>j</sup> for the function dS/dα
  2964.          */
  2965.         public T getdSdalphaCj(final int j) {
  2966.             return cjCoefs[4][j];
  2967.         }

  2969.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dα.
  2970.          * <br>
  2971.          * Possible values for j are within the interval [1,jMax].
  2972.          * @param j j index
  2973.          * @return S<sup>j</sup> for the derivative dS/dα
  2974.          */
  2975.         public T getdSdalphaSj(final int j) {
  2976.             return sjCoefs[4][j];
  2977.         }

  2979.         /** Get the coefficient C<sup>j</sup> for the derivative dS/dβ
  2980.          * <br>
  2981.          * Possible values for j are within the interval [0,jMax].
  2982.          * The value 0 is used to obtain the free coefficient C⁰
  2983.          * @param j j index
  2984.          * @return C<sup>j</sup> for the function dS/dβ
  2985.          */
  2986.         public T getdSdbetaCj(final int j) {
  2987.             return cjCoefs[5][j];
  2988.         }

  2990.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dβ.
  2991.          * <br>
  2992.          * Possible values for j are within the interval [1,jMax].
  2993.          * @param j j index
  2994.          * @return S<sup>j</sup> for the derivative dS/dβ
  2995.          */
  2996.         public T getdSdbetaSj(final int j) {
  2997.             return sjCoefs[5][j];
  2998.         }

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

  3011.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dγ.
  3012.          * <br>
  3013.          * Possible values for j are within the interval [1,jMax].
  3014.          * @param j j index
  3015.          * @return S<sup>j</sup> for the derivative dS/dγ
  3016.          */
  3017.         public T getdSdgammaSj(final int j) {
  3018.             return sjCoefs[6][j];
  3019.         }

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

  3032.         /** Get the coefficient S<sup>j</sup> for the derivative dS/dλ.
  3033.          * <br>
  3034.          * Possible values for j are within the interval [1,jMax].
  3035.          * @param j j index
  3036.          * @return S<sup>j</sup> for the derivative dS/dλ
  3037.          */
  3038.         public T getdSdlambdaSj(final int j) {
  3039.             return sjCoefs[7][j];
  3040.         }
  3041.     }

  3043.     /**
  3044.      * The coefficients used to compute the short periodic contribution for the Third body perturbation.
  3045.      * <p>
  3046.      * The short periodic contribution for the Third Body is expressed in Danielson 4.2-25.<br>
  3047.      * The coefficients C<sub>i</sub>⁰, C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup>
  3048.      * are computed by replacing the corresponding values in formula 2.5.5-10.
  3049.      * </p>
  3050.      * @author Lucian Barbulescu
  3051.      */
  3052.     private static class ThirdBodyShortPeriodicCoefficients implements ShortPeriodTerms {

  3054.         /** Maximal value for j. */
  3055.         private final int jMax;

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

  3060.         /** Max frequency of F. */
  3061.         private final int    maxFreqF;

  3063.         /** Coefficients prefix. */
  3064.         private final String prefix;

  3066.         /** All coefficients slots. */
  3067.         private final transient TimeSpanMap<Slot> slots;

  3069.         /**
  3070.          * Standard constructor.
  3071.          *  @param interpolationPoints number of points used in the interpolation process
  3072.          * @param jMax maximal value for j
  3073.          * @param maxFreqF Max frequency of F
  3074.          * @param bodyName third body name
  3075.          * @param slots all coefficients slots
  3076.          */
  3077.         ThirdBodyShortPeriodicCoefficients(final int jMax, final int interpolationPoints,
  3078.                                            final int maxFreqF, final String bodyName,
  3079.                                            final TimeSpanMap<Slot> slots) {
  3080.             this.jMax                = jMax;
  3081.             this.interpolationPoints = interpolationPoints;
  3082.             this.maxFreqF            = maxFreqF;
  3083.             this.prefix              = DSSTThirdBody.SHORT_PERIOD_PREFIX + bodyName + "-";
  3084.             this.slots               = slots;
  3085.         }

  3087.         /** Get the slot valid for some date.
  3088.          * @param meanStates mean states defining the slot
  3089.          * @return slot valid at the specified date
  3090.          */
  3091.         public Slot createSlot(final SpacecraftState... meanStates) {
  3092.             final Slot         slot  = new Slot(jMax, interpolationPoints);
  3093.             final AbsoluteDate first = meanStates[0].getDate();
  3094.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  3095.             if (first.compareTo(last) <= 0) {
  3096.                 slots.addValidAfter(slot, first);
  3097.             } else {
  3098.                 slots.addValidBefore(slot, first);
  3099.             }
  3100.             return slot;
  3101.         }

  3103.         /** {@inheritDoc} */
  3104.         @Override
  3105.         public double[] value(final Orbit meanOrbit) {

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

  3110.             // the current eccentric longitude
  3111.             final double F = meanOrbit.getLE();

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

  3116.             // Add the cos and sin dependent terms
  3117.             for (int j = 1; j <= maxFreqF; j++) {
  3118.                 //compute cos and sin
  3119.                 final SinCos scjF  = FastMath.sinCos(j * F);

  3121.                 final double[] c = slot.cij[j].value(meanOrbit.getDate());
  3122.                 final double[] s = slot.sij[j].value(meanOrbit.getDate());
  3123.                 for (int i = 0; i < 6; i++) {
  3124.                     shortPeriodic[i] += c[i] * scjF.cos() + s[i] * scjF.sin();
  3125.                 }
  3126.             }

  3128.             return shortPeriodic;

  3130.         }

  3132.         /** {@inheritDoc} */
  3133.         @Override
  3134.         public String getCoefficientsKeyPrefix() {
  3135.             return prefix;
  3136.         }

  3138.         /** {@inheritDoc}
  3139.          * <p>
  3140.          * For third body attraction forces,there are maxFreqF + 1 cj coefficients,
  3141.          * maxFreqF sj coefficients where maxFreqF depends on the orbit.
  3142.          * The j index is the integer multiplier for the eccentric longitude argument
  3143.          * in the cj and sj coefficients.
  3144.          * </p>
  3145.          */
  3146.         @Override
  3147.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected) {

  3149.             // select the coefficients slot
  3150.             final Slot slot = slots.get(date);

  3152.             final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * maxFreqF + 1);
  3153.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "c", 0);
  3154.             for (int j = 1; j <= maxFreqF; j++) {
  3155.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  3156.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  3157.             }
  3158.             return coefficients;

  3160.         }

  3162.         /** Put a coefficient in a map if selected.
  3163.          * @param map map to populate
  3164.          * @param selected set of coefficients that should be put in the map
  3165.          * (empty set means all coefficients are selected)
  3166.          * @param value coefficient value
  3167.          * @param id coefficient identifier
  3168.          * @param indices list of coefficient indices
  3169.          */
  3170.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
  3171.                                      final double[] value, final String id, final int... indices) {
  3172.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  3173.             keyBuilder.append(id);
  3174.             for (int index : indices) {
  3175.                 keyBuilder.append('[').append(index).append(']');
  3176.             }
  3177.             final String key = keyBuilder.toString();
  3178.             if (selected.isEmpty() || selected.contains(key)) {
  3179.                 map.put(key, value);
  3180.             }
  3181.         }

  3183.     }

  3185.     /**
  3186.      * The coefficients used to compute the short periodic contribution for the Third body perturbation.
  3187.      * <p>
  3188.      * The short periodic contribution for the Third Body is expressed in Danielson 4.2-25.<br>
  3189.      * The coefficients C<sub>i</sub>⁰, C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup>
  3190.      * are computed by replacing the corresponding values in formula 2.5.5-10.
  3191.      * </p>
  3192.      * @author Lucian Barbulescu
  3193.      */
  3194.     private static class FieldThirdBodyShortPeriodicCoefficients <T extends RealFieldElement<T>> implements FieldShortPeriodTerms<T> {

  3196.         /** Maximal value for j. */
  3197.         private final int jMax;

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

  3202.         /** Max frequency of F. */
  3203.         private final int    maxFreqF;

  3205.         /** Coefficients prefix. */
  3206.         private final String prefix;

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

  3211.         /**
  3212.          * Standard constructor.
  3213.          * @param interpolationPoints number of points used in the interpolation process
  3214.          * @param jMax maximal value for j
  3215.          * @param maxFreqF Max frequency of F
  3216.          * @param bodyName third body name
  3217.          * @param slots all coefficients slots
  3218.          */
  3219.         FieldThirdBodyShortPeriodicCoefficients(final int jMax, final int interpolationPoints,
  3220.                                                 final int maxFreqF, final String bodyName,
  3221.                                                 final FieldTimeSpanMap<FieldSlot<T>, T> slots) {
  3222.             this.jMax                = jMax;
  3223.             this.interpolationPoints = interpolationPoints;
  3224.             this.maxFreqF            = maxFreqF;
  3225.             this.prefix              = DSSTThirdBody.SHORT_PERIOD_PREFIX + bodyName + "-";
  3226.             this.slots               = slots;
  3227.         }

  3229.         /** Get the slot valid for some date.
  3230.          * @param meanStates mean states defining the slot
  3231.          * @return slot valid at the specified date
  3232.          */
  3233.         @SuppressWarnings("unchecked")
  3234.         public FieldSlot<T> createSlot(final FieldSpacecraftState<T>... meanStates) {
  3235.             final FieldSlot<T>         slot  = new FieldSlot<>(jMax, interpolationPoints);
  3236.             final FieldAbsoluteDate<T> first = meanStates[0].getDate();
  3237.             final FieldAbsoluteDate<T> last  = meanStates[meanStates.length - 1].getDate();
  3238.             if (first.compareTo(last) <= 0) {
  3239.                 slots.addValidAfter(slot, first);
  3240.             } else {
  3241.                 slots.addValidBefore(slot, first);
  3242.             }
  3243.             return slot;
  3244.         }

  3246.         /** {@inheritDoc} */
  3247.         @Override
  3248.         public T[] value(final FieldOrbit<T> meanOrbit) {

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

  3253.             // the current eccentric longitude
  3254.             final T F = meanOrbit.getLE();

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

  3259.             // Add the cos and sin dependent terms
  3260.             for (int j = 1; j <= maxFreqF; j++) {
  3261.                 //compute cos and sin
  3262.                 final FieldSinCos<T> scjF = FastMath.sinCos(F.multiply(j));

  3264.                 final T[] c = (T[]) slot.cij[j].value(meanOrbit.getDate());
  3265.                 final T[] s = (T[]) slot.sij[j].value(meanOrbit.getDate());
  3266.                 for (int i = 0; i < 6; i++) {
  3267.                     shortPeriodic[i] = shortPeriodic[i].add(c[i].multiply(scjF.cos()).add(s[i].multiply(scjF.sin())));
  3268.                 }
  3269.             }

  3271.             return shortPeriodic;

  3273.         }

  3275.         /** {@inheritDoc} */
  3276.         @Override
  3277.         public String getCoefficientsKeyPrefix() {
  3278.             return prefix;
  3279.         }

  3281.         /** {@inheritDoc}
  3282.          * <p>
  3283.          * For third body attraction forces,there are maxFreqF + 1 cj coefficients,
  3284.          * maxFreqF sj coefficients where maxFreqF depends on the orbit.
  3285.          * The j index is the integer multiplier for the eccentric longitude argument
  3286.          * in the cj and sj coefficients.
  3287.          * </p>
  3288.          */
  3289.         @Override
  3290.         public Map<String, T[]> getCoefficients(final FieldAbsoluteDate<T> date, final Set<String> selected) {

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

  3295.             final Map<String, T[]> coefficients = new HashMap<String, T[]>(2 * maxFreqF + 1);
  3296.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "c", 0);
  3297.             for (int j = 1; j <= maxFreqF; j++) {
  3298.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  3299.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  3300.             }
  3301.             return coefficients;

  3303.         }

  3305.         /** Put a coefficient in a map if selected.
  3306.          * @param map map to populate
  3307.          * @param selected set of coefficients that should be put in the map
  3308.          * (empty set means all coefficients are selected)
  3309.          * @param value coefficient value
  3310.          * @param id coefficient identifier
  3311.          * @param indices list of coefficient indices
  3312.          */
  3313.         private void storeIfSelected(final Map<String, T[]> map, final Set<String> selected,
  3314.                                      final T[] value, final String id, final int... indices) {
  3315.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  3316.             keyBuilder.append(id);
  3317.             for (int index : indices) {
  3318.                 keyBuilder.append('[').append(index).append(']');
  3319.             }
  3320.             final String key = keyBuilder.toString();
  3321.             if (selected.isEmpty() || selected.contains(key)) {
  3322.                 map.put(key, value);
  3323.             }
  3324.         }

  3326.     }

  3328.     /** Coefficients valid for one time slot. */
  3329.     private static class Slot {

  3331.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  3332.          * <p>
  3333.          * The index order is cij[j][i] <br/>
  3334.          * i corresponds to the equinoctial element, as follows: <br/>
  3335.          * - i=0 for a <br/>
  3336.          * - i=1 for k <br/>
  3337.          * - i=2 for h <br/>
  3338.          * - i=3 for q <br/>
  3339.          * - i=4 for p <br/>
  3340.          * - i=5 for λ <br/>
  3341.          * </p>
  3342.          */
  3343.         private final ShortPeriodicsInterpolatedCoefficient[] cij;

  3345.         /** The coefficients S<sub>i</sub><sup>j</sup>.
  3346.          * <p>
  3347.          * The index order is sij[j][i] <br/>
  3348.          * i corresponds to the equinoctial element, as follows: <br/>
  3349.          * - i=0 for a <br/>
  3350.          * - i=1 for k <br/>
  3351.          * - i=2 for h <br/>
  3352.          * - i=3 for q <br/>
  3353.          * - i=4 for p <br/>
  3354.          * - i=5 for λ <br/>
  3355.          * </p>
  3356.          */
  3357.         private final ShortPeriodicsInterpolatedCoefficient[] sij;

  3359.         /** Simple constructor.
  3360.          *  @param jMax maximum value for j index
  3361.          *  @param interpolationPoints number of points used in the interpolation process
  3362.          */
  3363.         Slot(final int jMax, final int interpolationPoints) {
  3364.             // allocate the coefficients arrays
  3365.             cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
  3366.             sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
  3367.             for (int j = 0; j <= jMax; j++) {
  3368.                 cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3369.                 sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3370.             }


  3373.         }
  3374.     }

  3376.     /** Coefficients valid for one time slot. */
  3377.     private static class FieldSlot <T extends RealFieldElement<T>> {

  3379.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  3380.          * <p>
  3381.          * The index order is cij[j][i] <br/>
  3382.          * i corresponds to the equinoctial element, as follows: <br/>
  3383.          * - i=0 for a <br/>
  3384.          * - i=1 for k <br/>
  3385.          * - i=2 for h <br/>
  3386.          * - i=3 for q <br/>
  3387.          * - i=4 for p <br/>
  3388.          * - i=5 for λ <br/>
  3389.          * </p>
  3390.          */
  3391.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] cij;

  3393.         /** The coefficients S<sub>i</sub><sup>j</sup>.
  3394.          * <p>
  3395.          * The index order is sij[j][i] <br/>
  3396.          * i corresponds to the equinoctial element, as follows: <br/>
  3397.          * - i=0 for a <br/>
  3398.          * - i=1 for k <br/>
  3399.          * - i=2 for h <br/>
  3400.          * - i=3 for q <br/>
  3401.          * - i=4 for p <br/>
  3402.          * - i=5 for λ <br/>
  3403.          * </p>
  3404.          */
  3405.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] sij;

  3407.         /** Simple constructor.
  3408.          *  @param jMax maximum value for j index
  3409.          *  @param interpolationPoints number of points used in the interpolation process
  3410.          */
  3411.         @SuppressWarnings("unchecked")
  3412.         FieldSlot(final int jMax, final int interpolationPoints) {
  3413.             // allocate the coefficients arrays
  3414.             cij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);
  3415.             sij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);
  3416.             for (int j = 0; j <= jMax; j++) {
  3417.                 cij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3418.                 sij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3419.             }


  3422.         }
  3423.     }

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

  3428.         /** The current value of the U function. <br/>
  3429.          * Needed for the short periodic contribution */
  3430.         private double U;

  3432.         /** dU / da. */
  3433.         private  double dUda;

  3435.         /** dU / dk. */
  3436.         private double dUdk;

  3438.         /** dU / dh. */
  3439.         private double dUdh;

  3441.         /** dU / dAlpha. */
  3442.         private double dUdAl;

  3444.         /** dU / dBeta. */
  3445.         private double dUdBe;

  3447.         /** dU / dGamma. */
  3448.         private double dUdGa;

  3450.         /** Simple constuctor.
  3451.          * @param context container for attributes
  3452.          * @param hansen hansen objects
  3453.          */
  3454.         UAnddU(final DSSTThirdBodyContext context,
  3455.                final HansenObjects hansen) {
  3456.             // Auxiliary elements related to the current orbit
  3457.             final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();

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

  3462.             // Initialise U.
  3463.             U = 0.;

  3465.             // Potential derivatives
  3466.             dUda  = 0.;
  3467.             dUdk  = 0.;
  3468.             dUdh  = 0.;
  3469.             dUdAl = 0.;
  3470.             dUdBe = 0.;
  3471.             dUdGa = 0.;

  3473.             for (int s = 0; s <= context.getMaxEccPow(); s++) {

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

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

  3481.                 // Compute Gs partial derivatives from 3.1-(9)
  3482.                 double dGsdh  = 0.;
  3483.                 double dGsdk  = 0.;
  3484.                 double dGsdAl = 0.;
  3485.                 double dGsdBe = 0.;
  3486.                 if (s > 0) {
  3487.                     // First get the G(s-1) and the H(s-1) coefficients
  3488.                     final double sxGsm1 = s * GsHs[0][s - 1];
  3489.                     final double sxHsm1 = s * GsHs[1][s - 1];
  3490.                     // Then compute derivatives
  3491.                     dGsdh  = context.getBeta()  * sxGsm1 - context.getAlpha() * sxHsm1;
  3492.                     dGsdk  = context.getAlpha() * sxGsm1 + context.getBeta()  * sxHsm1;
  3493.                     dGsdAl = auxiliaryElements.getK() * sxGsm1 - auxiliaryElements.getH() * sxHsm1;
  3494.                     dGsdBe = auxiliaryElements.getH() * sxGsm1 + auxiliaryElements.getK() * sxHsm1;
  3495.                 }

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

  3500.                 for (int n = FastMath.max(2, s); n <= context.getMaxAR3Pow(); n++) {
  3501.                     // (n - s) must be even
  3502.                     if ((n - s) % 2 == 0) {
  3503.                         // Extract data from previous computation :
  3504.                         final double kns   = hansen.getHansenObjects()[s].getValue(n, auxiliaryElements.getB());
  3505.                         final double dkns  = hansen.getHansenObjects()[s].getDerivative(n, auxiliaryElements.getB());

  3507.                         final double vns   = Vns.get(new NSKey(n, s));
  3508.                         final double coef0 = delta0s * context.getAoR3Pow()[n] * vns;
  3509.                         final double coef1 = coef0 * context.getQns()[n][s];
  3510.                         final double coef2 = coef1 * kns;
  3511.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  3512.                         // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
  3513.                         final double dqns = (n == s) ? 0. : context.getQns()[n][s + 1];

  3515.                         //Compute U:
  3516.                         U += coef2 * gs;

  3518.                         // Compute dU / da :
  3519.                         dUda  += coef2 * n * gs;
  3520.                         // Compute dU / dh
  3521.                         dUdh  += coef1 * (kns * dGsdh + context.getHXXX() * gs * dkns);
  3522.                         // Compute dU / dk
  3523.                         dUdk  += coef1 * (kns * dGsdk + context.getKXXX() * gs * dkns);
  3524.                         // Compute dU / dAlpha
  3525.                         dUdAl += coef2 * dGsdAl;
  3526.                         // Compute dU / dBeta
  3527.                         dUdBe += coef2 * dGsdBe;
  3528.                         // Compute dU / dGamma
  3529.                         dUdGa += coef0 * kns * dqns * gs;
  3530.                     }
  3531.                 }
  3532.             }

  3534.             // multiply by mu3 / R3
  3535.             this.U = U * context.getMuoR3();

  3537.             this.dUda  = dUda  * context.getMuoR3() / auxiliaryElements.getSma();
  3538.             this.dUdk  = dUdk  * context.getMuoR3();
  3539.             this.dUdh  = dUdh  * context.getMuoR3();
  3540.             this.dUdAl = dUdAl * context.getMuoR3();
  3541.             this.dUdBe = dUdBe * context.getMuoR3();
  3542.             this.dUdGa = dUdGa * context.getMuoR3();

  3544.         }

  3546.         /** Return value of U.
  3547.          * @return U
  3548.          */
  3549.         public double getU() {
  3550.             return U;
  3551.         }

  3553.         /** Return value of dU / da.
  3554.          * @return dUda
  3555.          */
  3556.         public double getdUda() {
  3557.             return dUda;
  3558.         }

  3560.         /** Return value of dU / dk.
  3561.          * @return dUdk
  3562.          */
  3563.         public double getdUdk() {
  3564.             return dUdk;
  3565.         }

  3567.         /** Return value of dU / dh.
  3568.          * @return dUdh
  3569.          */
  3570.         public double getdUdh() {
  3571.             return dUdh;
  3572.         }

  3574.         /** Return value of dU / dAlpha.
  3575.          * @return dUdAl
  3576.          */
  3577.         public double getdUdAl() {
  3578.             return dUdAl;
  3579.         }

  3581.         /** Return value of dU / dBeta.
  3582.          * @return dUdBe
  3583.          */
  3584.         public double getdUdBe() {
  3585.             return dUdBe;
  3586.         }

  3588.         /** Return value of dU / dGamma.
  3589.          * @return dUdGa
  3590.          */
  3591.         public double getdUdGa() {
  3592.             return dUdGa;
  3593.         }

  3595.     }

  3597.     /** Compute potential and potential derivatives with respect to orbital parameters. */
  3598.     private class FieldUAnddU <T extends RealFieldElement<T>> {

  3600.         /** The current value of the U function. <br/>
  3601.          * Needed for the short periodic contribution */
  3602.         private T U;

  3604.         /** dU / da. */
  3605.         private T dUda;

  3607.         /** dU / dk. */
  3608.         private T dUdk;

  3610.         /** dU / dh. */
  3611.         private T dUdh;

  3613.         /** dU / dAlpha. */
  3614.         private T dUdAl;

  3616.         /** dU / dBeta. */
  3617.         private T dUdBe;

  3619.         /** dU / dGamma. */
  3620.         private T dUdGa;

  3622.         /** Simple constuctor.
  3623.          * @param context container for attributes
  3624.          * @param hansen hansen objects
  3625.          */
  3626.         FieldUAnddU(final FieldDSSTThirdBodyContext<T> context,
  3627.                     final FieldHansenObjects<T> hansen) {

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

  3632.             // Field for array building
  3633.             final Field<T> field = auxiliaryElements.getDate().getField();
  3634.             // Zero for initialization
  3635.             final T zero = field.getZero();

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

  3640.             // Initialise U.
  3641.             U = zero;

  3643.             // Potential derivatives
  3644.             dUda  = zero;
  3645.             dUdk  = zero;
  3646.             dUdh  = zero;
  3647.             dUdAl = zero;
  3648.             dUdBe = zero;
  3649.             dUdGa = zero;

  3651.             for (int s = 0; s <= context.getMaxEccPow(); s++) {
  3652.                 // initialise the Hansen roots
  3653.                 hansen.computeHansenObjectsInitValues(context, auxiliaryElements.getB(), s);

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

  3658.                 // Compute Gs partial derivatives from 3.1-(9)
  3659.                 T dGsdh  = zero;
  3660.                 T dGsdk  = zero;
  3661.                 T dGsdAl = zero;
  3662.                 T dGsdBe = zero;
  3663.                 if (s > 0) {
  3664.                     // First get the G(s-1) and the H(s-1) coefficients
  3665.                     final T sxGsm1 = GsHs[0][s - 1].multiply(s);
  3666.                     final T sxHsm1 = GsHs[1][s - 1].multiply(s);
  3667.                     // Then compute derivatives
  3668.                     dGsdh  = sxGsm1.multiply(context.getBeta()).subtract(sxHsm1.multiply(context.getAlpha()));
  3669.                     dGsdk  = sxGsm1.multiply(context.getAlpha()).add(sxHsm1.multiply(context.getBeta()));
  3670.                     dGsdAl = sxGsm1.multiply(auxiliaryElements.getK()).subtract(sxHsm1.multiply(auxiliaryElements.getH()));
  3671.                     dGsdBe = sxGsm1.multiply(auxiliaryElements.getH()).add(sxHsm1.multiply(auxiliaryElements.getK()));
  3672.                 }

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

  3677.                 for (int n = FastMath.max(2, s); n <= context.getMaxAR3Pow(); n++) {
  3678.                     // (n - s) must be even
  3679.                     if ((n - s) % 2 == 0) {
  3680.                         // Extract data from previous computation :
  3681.                         final T kns   = (T) hansen.getHansenObjects()[s].getValue(n, auxiliaryElements.getB());
  3682.                         final T dkns  = (T) hansen.getHansenObjects()[s].getDerivative(n, auxiliaryElements.getB());

  3684.                         final double vns = Vns.get(new NSKey(n, s));
  3685.                         final T coef0 = delta0s.multiply(vns).multiply(context.getAoR3Pow()[n]);
  3686.                         final T coef1 = coef0.multiply(context.getQns()[n][s]);
  3687.                         final T coef2 = coef1.multiply(kns);
  3688.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  3689.                         // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
  3690.                         final T dqns = (n == s) ? zero : context.getQns()[n][s + 1];

  3692.                         //Compute U:
  3693.                         U = U.add(coef2.multiply(gs));

  3695.                         // Compute dU / da :
  3696.                         dUda  = dUda.add(coef2.multiply(n).multiply(gs));
  3697.                         // Compute dU / dh
  3698.                         dUdh  = dUdh.add(coef1.multiply(dGsdh.multiply(kns).add(context.getHXXX().multiply(gs).multiply(dkns))));
  3699.                         // Compute dU / dk
  3700.                         dUdk  = dUdk.add(coef1.multiply(dGsdk.multiply(kns).add(context.getKXXX().multiply(gs).multiply(dkns))));
  3701.                         // Compute dU / dAlpha
  3702.                         dUdAl = dUdAl.add(coef2.multiply(dGsdAl));
  3703.                         // Compute dU / dBeta
  3704.                         dUdBe = dUdBe.add(coef2.multiply(dGsdBe));
  3705.                         // Compute dU / dGamma
  3706.                         dUdGa = dUdGa.add(coef0.multiply(kns).multiply(dqns).multiply(gs));
  3707.                     }
  3708.                 }
  3709.             }

  3711.             // multiply by mu3 / R3
  3712.             this.U = U.multiply(context.getMuoR3());

  3714.             this.dUda  = dUda.multiply(context.getMuoR3().divide(auxiliaryElements.getSma()));
  3715.             this.dUdk  = dUdk.multiply(context.getMuoR3());
  3716.             this.dUdh  = dUdh.multiply(context.getMuoR3());
  3717.             this.dUdAl = dUdAl.multiply(context.getMuoR3());
  3718.             this.dUdBe = dUdBe.multiply(context.getMuoR3());
  3719.             this.dUdGa = dUdGa.multiply(context.getMuoR3());

  3721.         }

  3723.         /** Return value of U.
  3724.          * @return U
  3725.          */
  3726.         public T getU() {
  3727.             return U;
  3728.         }

  3730.         /** Return value of dU / da.
  3731.          * @return dUda
  3732.          */
  3733.         public T getdUda() {
  3734.             return dUda;
  3735.         }

  3737.         /** Return value of dU / dk.
  3738.          * @return dUdk
  3739.          */
  3740.         public T getdUdk() {
  3741.             return dUdk;
  3742.         }

  3744.         /** Return value of dU / dh.
  3745.          * @return dUdh
  3746.          */
  3747.         public T getdUdh() {
  3748.             return dUdh;
  3749.         }

  3751.         /** Return value of dU / dAlpha.
  3752.          * @return dUdAl
  3753.          */
  3754.         public T getdUdAl() {
  3755.             return dUdAl;
  3756.         }

  3758.         /** Return value of dU / dBeta.
  3759.          * @return dUdBe
  3760.          */
  3761.         public T getdUdBe() {
  3762.             return dUdBe;
  3763.         }

  3765.         /** Return value of dU / dGamma.
  3766.          * @return dUdGa
  3767.          */
  3768.         public T getdUdGa() {
  3769.             return dUdGa;
  3770.         }

  3772.     }

  3774.     /** Computes init values of the Hansen Objects. */
  3775.     private static class HansenObjects {

  3777.         /** Max power for summation. */
  3778.         private static final int    MAX_POWER = 22;

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

  3784.         /** Simple constructor. */
  3785.         HansenObjects() {
  3786.             this.hansenObjects = new HansenThirdBodyLinear[MAX_POWER + 1];
  3787.             for (int s = 0; s <= MAX_POWER; s++) {
  3788.                 this.hansenObjects[s] = new HansenThirdBodyLinear(MAX_POWER, s);
  3789.             }
  3790.         }

  3792.         /** Compute init values for hansen objects.
  3793.          * @param context container for attributes
  3794.          * @param B = sqrt(1 - e²).
  3795.          * @param element element of the array to compute the init values
  3796.          */
  3797.         public void computeHansenObjectsInitValues(final DSSTThirdBodyContext context, final double B, final int element) {
  3798.             hansenObjects[element].computeInitValues(B, context.getBB(), context.getBBB());
  3799.         }

  3801.         /** Get the Hansen Objects.
  3802.          * @return hansenObjects
  3803.          */
  3804.         public HansenThirdBodyLinear[] getHansenObjects() {
  3805.             return hansenObjects;
  3806.         }

  3808.     }

  3810.     /** Computes init values of the Hansen Objects. */
  3811.     private static class FieldHansenObjects<T extends RealFieldElement<T>> {

  3813.         /** Max power for summation. */
  3814.         private static final int    MAX_POWER = 22;

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

  3820.         /** Simple constructor.
  3821.          * @param field field used by default
  3822.          */
  3823.         @SuppressWarnings("unchecked")
  3824.         FieldHansenObjects(final Field<T> field) {
  3825.             this.hansenObjects = (FieldHansenThirdBodyLinear<T>[]) Array.newInstance(FieldHansenThirdBodyLinear.class, MAX_POWER + 1);
  3826.             for (int s = 0; s <= MAX_POWER; s++) {
  3827.                 this.hansenObjects[s] = new FieldHansenThirdBodyLinear<>(MAX_POWER, s, field);
  3828.             }
  3829.         }

  3831.         /** Initialise the Hansen roots for third body problem.
  3832.          * @param context container for attributes
  3833.          * @param B = sqrt(1 - e²).
  3834.          * @param element element of the array to compute the init values
  3835.          */
  3836.         public void computeHansenObjectsInitValues(final FieldDSSTThirdBodyContext<T> context,
  3837.                                                    final T B, final int element) {
  3838.             hansenObjects[element].computeInitValues(B, context.getBB(), context.getBBB());
  3839.         }

  3841.         /** Get the Hansen Objects.
  3842.          * @return hansenObjects
  3843.          */
  3844.         public FieldHansenThirdBodyLinear<T>[] getHansenObjects() {
  3845.             return hansenObjects;
  3846.         }

  3848.     }

  3850. }
Created with JaCoCo 0.8.6.202009150832