DSSTZonal.java

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

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

  29. import org.hipparchus.CalculusFieldElement;
  30. import org.hipparchus.Field;
  31. import org.hipparchus.exception.LocalizedCoreFormats;
  32. import org.hipparchus.util.CombinatoricsUtils;
  33. import org.hipparchus.util.FastMath;
  34. import org.hipparchus.util.FieldSinCos;
  35. import org.hipparchus.util.MathArrays;
  36. import org.hipparchus.util.SinCos;
  37. import org.orekit.attitudes.AttitudeProvider;
  38. import org.orekit.errors.OrekitException;
  39. import org.orekit.errors.OrekitInternalError;
  40. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  41. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
  42. import org.orekit.orbits.FieldOrbit;
  43. import org.orekit.orbits.Orbit;
  44. import org.orekit.propagation.FieldSpacecraftState;
  45. import org.orekit.propagation.PropagationType;
  46. import org.orekit.propagation.SpacecraftState;
  47. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  48. import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
  49. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
  50. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory.NSKey;
  51. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  52. import org.orekit.propagation.semianalytical.dsst.utilities.FieldCjSjCoefficient;
  53. import org.orekit.propagation.semianalytical.dsst.utilities.FieldGHIJjsPolynomials;
  54. import org.orekit.propagation.semianalytical.dsst.utilities.FieldLnsCoefficients;
  55. import org.orekit.propagation.semianalytical.dsst.utilities.FieldShortPeriodicsInterpolatedCoefficient;
  56. import org.orekit.propagation.semianalytical.dsst.utilities.GHIJjsPolynomials;
  57. import org.orekit.propagation.semianalytical.dsst.utilities.LnsCoefficients;
  58. import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
  59. import org.orekit.propagation.semianalytical.dsst.utilities.UpperBounds;
  60. import org.orekit.propagation.semianalytical.dsst.utilities.hansen.FieldHansenZonalLinear;
  61. import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenZonalLinear;
  62. import org.orekit.time.AbsoluteDate;
  63. import org.orekit.time.FieldAbsoluteDate;
  64. import org.orekit.utils.FieldTimeSpanMap;
  65. import org.orekit.utils.ParameterDriver;
  66. import org.orekit.utils.TimeSpanMap;

  68. /** Zonal contribution to the central body gravitational perturbation.
  69.  *
  70.  *   @author Romain Di Costanzo
  71.  *   @author Pascal Parraud
  72.  *   @author Bryan Cazabonne (field translation)
  73.  */
  74. public class DSSTZonal implements DSSTForceModel {

  76.     /**  Name of the prefix for short period coefficients keys. */
  77.     public static final String SHORT_PERIOD_PREFIX = "DSST-central-body-zonal-";

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

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

  97.     /** Truncation tolerance. */
  98.     private static final double TRUNCATION_TOLERANCE = 1e-4;

  100.     /** Central attraction scaling factor.
  101.      * <p>
  102.      * We use a power of 2 to avoid numeric noise introduction
  103.      * in the multiplications/divisions sequences.
  104.      * </p>
  105.      */
  106.     private static final double MU_SCALE = FastMath.scalb(1.0, 32);

  108.     /** Coefficient used to define the mean disturbing function V<sub>ns</sub> coefficient. */
  109.     private final SortedMap<NSKey, Double> Vns;

  111.     /** Provider for spherical harmonics. */
  112.     private final UnnormalizedSphericalHarmonicsProvider provider;

  114.     /** Maximal degree to consider for harmonics potential. */
  115.     private final int maxDegree;

  117.     /** Maximal degree to consider for harmonics potential. */
  118.     private final int maxDegreeShortPeriodics;

  120.     /** Highest power of the eccentricity to be used in short periodic computations. */
  121.     private final int maxEccPowShortPeriodics;

  123.     /** Maximum frequency in true longitude for short periodic computations. */
  124.     private final int maxFrequencyShortPeriodics;

  126.     /** Highest power of the eccentricity to be used in mean elements computations. */
  127.     private int maxEccPowMeanElements;

  129.     /** Highest power of the eccentricity. */
  130.     private int maxEccPow;

  132.     /** Short period terms. */
  133.     private ZonalShortPeriodicCoefficients zonalSPCoefs;

  135.     /** Short period terms. */
  136.     private Map<Field<?>, FieldZonalShortPeriodicCoefficients<?>> zonalFieldSPCoefs;

  138.     /** Driver for gravitational parameter. */
  139.     private final ParameterDriver gmParameterDriver;

  141.     /** Hansen objects. */
  142.     private HansenObjects hansen;

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

  147.     /** Constructor with default reference values.
  148.      * <p>
  149.      * When this constructor is used, maximum allowed values are used
  150.      * for the short periodic coefficients:
  151.      * </p>
  152.      * <ul>
  153.      *    <li> {@link #maxDegreeShortPeriodics} is set to {@code provider.getMaxDegree()} </li>
  154.      *    <li> {@link #maxEccPowShortPeriodics} is set to {@code min(provider.getMaxDegree() - 1, 4)}.
  155.      *         This parameter should not exceed 4 as higher values will exceed computer capacity </li>
  156.      *    <li> {@link #maxFrequencyShortPeriodics} is set to {@code 2 * provider.getMaxDegree() + 1} </li>
  157.      * </ul>
  158.      * @param provider provider for spherical harmonics
  159.      * @since 10.1
  160.      */
  161.     public DSSTZonal(final UnnormalizedSphericalHarmonicsProvider provider) {
  162.         this(provider, provider.getMaxDegree(), FastMath.min(4, provider.getMaxDegree() - 1), 2 * provider.getMaxDegree() + 1);
  163.     }

  165.     /** Simple constructor.
  166.      * @param provider provider for spherical harmonics
  167.      * @param maxDegreeShortPeriodics maximum degree to consider for short periodics zonal harmonics potential
  168.      * (must be between 2 and {@code provider.getMaxDegree()})
  169.      * @param maxEccPowShortPeriodics maximum power of the eccentricity to be used in short periodic computations
  170.      * (must be between 0 and {@code maxDegreeShortPeriodics - 1}, but should typically not exceed 4 as higher
  171.      * values will exceed computer capacity)
  172.      * @param maxFrequencyShortPeriodics maximum frequency in true longitude for short periodic computations
  173.      * (must be between 1 and {@code 2 * maxDegreeShortPeriodics + 1})
  174.      * @since 7.2
  175.      */
  176.     public DSSTZonal(final UnnormalizedSphericalHarmonicsProvider provider,
  177.                      final int maxDegreeShortPeriodics,
  178.                      final int maxEccPowShortPeriodics,
  179.                      final int maxFrequencyShortPeriodics) {

  181.         gmParameterDriver = new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
  182.                                                 provider.getMu(), MU_SCALE,
  183.                                                 0.0, Double.POSITIVE_INFINITY);

  185.         // Vns coefficients
  186.         this.Vns = CoefficientsFactory.computeVns(provider.getMaxDegree() + 1);

  188.         this.provider  = provider;
  189.         this.maxDegree = provider.getMaxDegree();

  191.         checkIndexRange(maxDegreeShortPeriodics, 2, provider.getMaxDegree());
  192.         this.maxDegreeShortPeriodics = maxDegreeShortPeriodics;

  194.         checkIndexRange(maxEccPowShortPeriodics, 0, maxDegreeShortPeriodics - 1);
  195.         this.maxEccPowShortPeriodics = maxEccPowShortPeriodics;

  197.         checkIndexRange(maxFrequencyShortPeriodics, 1, 2 * maxDegreeShortPeriodics + 1);
  198.         this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;

  200.         // Initialize default values
  201.         this.maxEccPowMeanElements = (maxDegree == 2) ? 0 : Integer.MIN_VALUE;

  203.         zonalFieldSPCoefs = new HashMap<>();
  204.         fieldHansen       = new HashMap<>();

  206.     }

  208.     /** Check an index range.
  209.      * @param index index value
  210.      * @param min minimum value for index
  211.      * @param max maximum value for index
  212.      */
  213.     private void checkIndexRange(final int index, final int min, final int max) {
  214.         if (index < min || index > max) {
  215.             throw new OrekitException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, index, min, max);
  216.         }
  217.     }

  219.     /** Get the spherical harmonics provider.
  220.      *  @return the spherical harmonics provider
  221.      */
  222.     public UnnormalizedSphericalHarmonicsProvider getProvider() {
  223.         return provider;
  224.     }

  226.     /** {@inheritDoc}
  227.      *  <p>
  228.      *  Computes the highest power of the eccentricity to appear in the truncated
  229.      *  analytical power series expansion.
  230.      *  </p>
  231.      *  <p>
  232.      *  This method computes the upper value for the central body potential and
  233.      *  determines the maximal power for the eccentricity producing potential
  234.      *  terms bigger than a defined tolerance.
  235.      *  </p>
  236.      */
  237.     @Override
  238.     public List<ShortPeriodTerms> initializeShortPeriodTerms(final AuxiliaryElements auxiliaryElements,
  239.                                              final PropagationType type,
  240.                                              final double[] parameters) {

  242.         computeMeanElementsTruncations(auxiliaryElements, parameters);

  244.         switch (type) {
  245.             case MEAN:
  246.                 maxEccPow = maxEccPowMeanElements;
  247.                 break;
  248.             case OSCULATING:
  249.                 maxEccPow = FastMath.max(maxEccPowMeanElements, maxEccPowShortPeriodics);
  250.                 break;
  251.             default:
  252.                 throw new OrekitInternalError(null);
  253.         }

  255.         hansen = new HansenObjects();

  257.         final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
  258.         zonalSPCoefs = new ZonalShortPeriodicCoefficients(maxFrequencyShortPeriodics,
  259.                                                           INTERPOLATION_POINTS,
  260.                                                           new TimeSpanMap<Slot>(new Slot(maxFrequencyShortPeriodics,
  261.                                                                                          INTERPOLATION_POINTS)));
  262.         list.add(zonalSPCoefs);
  263.         return list;

  265.     }

  267.     /** {@inheritDoc}
  268.      *  <p>
  269.      *  Computes the highest power of the eccentricity to appear in the truncated
  270.      *  analytical power series expansion.
  271.      *  </p>
  272.      *  <p>
  273.      *  This method computes the upper value for the central body potential and
  274.      *  determines the maximal power for the eccentricity producing potential
  275.      *  terms bigger than a defined tolerance.
  276.      *  </p>
  277.      */
  278.     @Override
  279.     public <T extends CalculusFieldElement<T>> List<FieldShortPeriodTerms<T>> initializeShortPeriodTerms(final FieldAuxiliaryElements<T> auxiliaryElements,
  280.                                                                                      final PropagationType type,
  281.                                                                                      final T[] parameters) {

  283.         // Field used by default
  284.         final Field<T> field = auxiliaryElements.getDate().getField();
  285.         computeMeanElementsTruncations(auxiliaryElements, parameters, field);

  287.         switch (type) {
  288.             case MEAN:
  289.                 maxEccPow = maxEccPowMeanElements;
  290.                 break;
  291.             case OSCULATING:
  292.                 maxEccPow = FastMath.max(maxEccPowMeanElements, maxEccPowShortPeriodics);
  293.                 break;
  294.             default:
  295.                 throw new OrekitInternalError(null);
  296.         }

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

  300.         final FieldZonalShortPeriodicCoefficients<T> fzspc =
  301.                         new FieldZonalShortPeriodicCoefficients<>(maxFrequencyShortPeriodics,
  302.                                                                   INTERPOLATION_POINTS,
  303.                                                                   new FieldTimeSpanMap<>(new FieldSlot<>(maxFrequencyShortPeriodics,
  304.                                                                                                          INTERPOLATION_POINTS),
  305.                                                                                          field));
  306.         zonalFieldSPCoefs.put(field, fzspc);
  307.         return Collections.singletonList(fzspc);

  309.     }

  311.     /** Compute indices truncations for mean elements computations.
  312.      * @param auxiliaryElements auxiliary elements
  313.      * @param parameters values of the force model parameters for state date (only 1 value for each parameter)
  314.      */
  315.     private void computeMeanElementsTruncations(final AuxiliaryElements auxiliaryElements, final double[] parameters) {

  317.         final DSSTZonalContext context = new DSSTZonalContext(auxiliaryElements, provider, parameters);
  318.         //Compute the max eccentricity power for the mean element rate expansion
  319.         if (maxDegree == 2) {
  320.             maxEccPowMeanElements = 0;
  321.         } else {
  322.             // Initializes specific parameters.
  323.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(auxiliaryElements.getDate());

  325.             // Utilities for truncation
  326.             final double ax2or = 2. * auxiliaryElements.getSma() / provider.getAe();
  327.             double xmuran = parameters[0] / auxiliaryElements.getSma();
  328.             // Set a lower bound for eccentricity
  329.             final double eo2  = FastMath.max(0.0025, auxiliaryElements.getEcc() / 2.);
  330.             final double x2o2 = context.getXX() / 2.;
  331.             final double[] eccPwr = new double[maxDegree + 1];
  332.             final double[] chiPwr = new double[maxDegree + 1];
  333.             final double[] hafPwr = new double[maxDegree + 1];
  334.             eccPwr[0] = 1.;
  335.             chiPwr[0] = context.getX();
  336.             hafPwr[0] = 1.;
  337.             for (int i = 1; i <= maxDegree; i++) {
  338.                 eccPwr[i] = eccPwr[i - 1] * eo2;
  339.                 chiPwr[i] = chiPwr[i - 1] * x2o2;
  340.                 hafPwr[i] = hafPwr[i - 1] * 0.5;
  341.                 xmuran  /= ax2or;
  342.             }

  344.             // Set highest power of e and degree of current spherical harmonic.
  345.             maxEccPowMeanElements = 0;
  346.             int maxDeg = maxDegree;
  347.             // Loop over n
  348.             do {
  349.                 // Set order of current spherical harmonic.
  350.                 int m = 0;
  351.                 // Loop over m
  352.                 do {
  353.                     // Compute magnitude of current spherical harmonic coefficient.
  354.                     final double cnm = harmonics.getUnnormalizedCnm(maxDeg, m);
  355.                     final double snm = harmonics.getUnnormalizedSnm(maxDeg, m);
  356.                     final double csnm = FastMath.hypot(cnm, snm);
  357.                     if (csnm == 0.) {
  358.                         break;
  359.                     }
  360.                     // Set magnitude of last spherical harmonic term.
  361.                     double lastTerm = 0.;
  362.                     // Set current power of e and related indices.
  363.                     int nsld2 = (maxDeg - maxEccPowMeanElements - 1) / 2;
  364.                     int l = maxDeg - 2 * nsld2;
  365.                     // Loop over l
  366.                     double term = 0.;
  367.                     do {
  368.                         // Compute magnitude of current spherical harmonic term.
  369.                         if (m < l) {
  370.                             term =  csnm * xmuran *
  371.                                     (CombinatoricsUtils.factorialDouble(maxDeg - l) / (CombinatoricsUtils.factorialDouble(maxDeg - m))) *
  372.                                     (CombinatoricsUtils.factorialDouble(maxDeg + l) / (CombinatoricsUtils.factorialDouble(nsld2) * CombinatoricsUtils.factorialDouble(nsld2 + l))) *
  373.                                     eccPwr[l] * UpperBounds.getDnl(context.getXX(), chiPwr[l], maxDeg, l) *
  374.                                     (UpperBounds.getRnml(auxiliaryElements.getGamma(), maxDeg, l, m, 1, I) + UpperBounds.getRnml(auxiliaryElements.getGamma(), maxDeg, l, m, -1, I));
  375.                         } else {
  376.                             term =  csnm * xmuran *
  377.                                     (CombinatoricsUtils.factorialDouble(maxDeg + m) / (CombinatoricsUtils.factorialDouble(nsld2) * CombinatoricsUtils.factorialDouble(nsld2 + l))) *
  378.                                     eccPwr[l] * hafPwr[m - l] * UpperBounds.getDnl(context.getXX(), chiPwr[l], maxDeg, l) *
  379.                                     (UpperBounds.getRnml(auxiliaryElements.getGamma(), maxDeg, m, l, 1, I) + UpperBounds.getRnml(auxiliaryElements.getGamma(), maxDeg, m, l, -1, I));
  380.                         }
  381.                         // Is the current spherical harmonic term bigger than the truncation tolerance ?
  382.                         if (term >= TRUNCATION_TOLERANCE) {
  383.                             maxEccPowMeanElements = l;
  384.                         } else {
  385.                             // Is the current term smaller than the last term ?
  386.                             if (term < lastTerm) {
  387.                                 break;
  388.                             }
  389.                         }
  390.                         // Proceed to next power of e.
  391.                         lastTerm = term;
  392.                         l += 2;
  393.                         nsld2--;
  394.                     } while (l < maxDeg);
  395.                     // Is the current spherical harmonic term bigger than the truncation tolerance ?
  396.                     if (term >= TRUNCATION_TOLERANCE) {
  397.                         maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
  398.                         return;
  399.                     }
  400.                     // Proceed to next order.
  401.                     m++;
  402.                 } while (m <= FastMath.min(maxDeg, provider.getMaxOrder()));
  403.                 // Proceed to next degree.
  404.                 xmuran *= ax2or;
  405.                 maxDeg--;
  406.             } while (maxDeg > maxEccPowMeanElements + 2);

  408.             maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
  409.         }
  410.     }

  412.     /** Compute indices truncations for mean elements computations.
  413.      * @param <T> type of the elements
  414.      * @param auxiliaryElements auxiliary elements
  415.      * @param parameters values of the force model parameters
  416.      * @param field field used by default
  417.      */
  418.     private <T extends CalculusFieldElement<T>> void computeMeanElementsTruncations(final FieldAuxiliaryElements<T> auxiliaryElements,
  419.                                                                                 final T[] parameters,
  420.                                                                                 final Field<T> field) {

  422.         final T zero = field.getZero();
  423.         final FieldDSSTZonalContext<T> context = new FieldDSSTZonalContext<>(auxiliaryElements, provider, parameters);
  424.         //Compute the max eccentricity power for the mean element rate expansion
  425.         if (maxDegree == 2) {
  426.             maxEccPowMeanElements = 0;
  427.         } else {
  428.             // Initializes specific parameters.
  429.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(auxiliaryElements.getDate().toAbsoluteDate());

  431.             // Utilities for truncation
  432.             final T ax2or = auxiliaryElements.getSma().multiply(2.).divide(provider.getAe());
  433.             T xmuran = parameters[0].divide(auxiliaryElements.getSma());
  434.             // Set a lower bound for eccentricity
  435.             final T eo2  = FastMath.max(zero.add(0.0025), auxiliaryElements.getEcc().divide(2.));
  436.             final T x2o2 = context.getXX().divide(2.);
  437.             final T[] eccPwr = MathArrays.buildArray(field, maxDegree + 1);
  438.             final T[] chiPwr = MathArrays.buildArray(field, maxDegree + 1);
  439.             final T[] hafPwr = MathArrays.buildArray(field, maxDegree + 1);
  440.             eccPwr[0] = zero.add(1.);
  441.             chiPwr[0] = context.getX();
  442.             hafPwr[0] = zero.add(1.);
  443.             for (int i = 1; i <= maxDegree; i++) {
  444.                 eccPwr[i] = eccPwr[i - 1].multiply(eo2);
  445.                 chiPwr[i] = chiPwr[i - 1].multiply(x2o2);
  446.                 hafPwr[i] = hafPwr[i - 1].multiply(0.5);
  447.                 xmuran  = xmuran.divide(ax2or);
  448.             }

  450.             // Set highest power of e and degree of current spherical harmonic.
  451.             maxEccPowMeanElements = 0;
  452.             int maxDeg = maxDegree;
  453.             // Loop over n
  454.             do {
  455.                 // Set order of current spherical harmonic.
  456.                 int m = 0;
  457.                 // Loop over m
  458.                 do {
  459.                     // Compute magnitude of current spherical harmonic coefficient.
  460.                     final T cnm = zero.add(harmonics.getUnnormalizedCnm(maxDeg, m));
  461.                     final T snm = zero.add(harmonics.getUnnormalizedSnm(maxDeg, m));
  462.                     final T csnm = FastMath.hypot(cnm, snm);
  463.                     if (csnm.getReal() == 0.) {
  464.                         break;
  465.                     }
  466.                     // Set magnitude of last spherical harmonic term.
  467.                     T lastTerm = zero;
  468.                     // Set current power of e and related indices.
  469.                     int nsld2 = (maxDeg - maxEccPowMeanElements - 1) / 2;
  470.                     int l = maxDeg - 2 * nsld2;
  471.                     // Loop over l
  472.                     T term = zero;
  473.                     do {
  474.                         // Compute magnitude of current spherical harmonic term.
  475.                         if (m < l) {
  476.                             term = csnm.multiply(xmuran).
  477.                                    multiply((CombinatoricsUtils.factorialDouble(maxDeg - l) / (CombinatoricsUtils.factorialDouble(maxDeg - m))) *
  478.                                    (CombinatoricsUtils.factorialDouble(maxDeg + l) / (CombinatoricsUtils.factorialDouble(nsld2) * CombinatoricsUtils.factorialDouble(nsld2 + l)))).
  479.                                    multiply(eccPwr[l]).multiply(UpperBounds.getDnl(context.getXX(), chiPwr[l], maxDeg, l)).
  480.                                    multiply(UpperBounds.getRnml(auxiliaryElements.getGamma(), maxDeg, l, m, 1, I).add(UpperBounds.getRnml(auxiliaryElements.getGamma(), maxDeg, l, m, -1, I)));
  481.                         } else {
  482.                             term = csnm.multiply(xmuran).
  483.                                    multiply(CombinatoricsUtils.factorialDouble(maxDeg + m) / (CombinatoricsUtils.factorialDouble(nsld2) * CombinatoricsUtils.factorialDouble(nsld2 + l))).
  484.                                    multiply(eccPwr[l]).multiply(hafPwr[m - l]).multiply(UpperBounds.getDnl(context.getXX(), chiPwr[l], maxDeg, l)).
  485.                                    multiply(UpperBounds.getRnml(auxiliaryElements.getGamma(), maxDeg, m, l, 1, I).add(UpperBounds.getRnml(auxiliaryElements.getGamma(), maxDeg, m, l, -1, I)));
  486.                         }
  487.                         // Is the current spherical harmonic term bigger than the truncation tolerance ?
  488.                         if (term.getReal() >= TRUNCATION_TOLERANCE) {
  489.                             maxEccPowMeanElements = l;
  490.                         } else {
  491.                             // Is the current term smaller than the last term ?
  492.                             if (term.getReal() < lastTerm.getReal()) {
  493.                                 break;
  494.                             }
  495.                         }
  496.                         // Proceed to next power of e.
  497.                         lastTerm = term;
  498.                         l += 2;
  499.                         nsld2--;
  500.                     } while (l < maxDeg);
  501.                     // Is the current spherical harmonic term bigger than the truncation tolerance ?
  502.                     if (term.getReal() >= TRUNCATION_TOLERANCE) {
  503.                         maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
  504.                         return;
  505.                     }
  506.                     // Proceed to next order.
  507.                     m++;
  508.                 } while (m <= FastMath.min(maxDeg, provider.getMaxOrder()));
  509.                 // Proceed to next degree.
  510.                 xmuran = xmuran.multiply(ax2or);
  511.                 maxDeg--;
  512.             } while (maxDeg > maxEccPowMeanElements + 2);

  514.             maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
  515.         }
  516.     }

  518.     /** Performs initialization at each integration step for the current force model.
  519.      *  <p>
  520.      *  This method aims at being called before mean elements rates computation.
  521.      *  </p>
  522.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  523.      *  @param parameters values of the force model parameters
  524.      *  @return new force model context
  525.      */
  526.     private DSSTZonalContext initializeStep(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
  527.         return new DSSTZonalContext(auxiliaryElements, provider, parameters);
  528.     }

  530.     /** Performs initialization at each integration step for the current force model.
  531.      *  <p>
  532.      *  This method aims at being called before mean elements rates computation.
  533.      *  </p>
  534.      *  @param <T> type of the elements
  535.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  536.      *  @param parameters values of the force model parameters
  537.      *  @return new force model context
  538.      */
  539.     private <T extends CalculusFieldElement<T>> FieldDSSTZonalContext<T> initializeStep(final FieldAuxiliaryElements<T> auxiliaryElements,
  540.                                                                                     final T[] parameters) {
  541.         return new FieldDSSTZonalContext<>(auxiliaryElements, provider, parameters);
  542.     }

  544.     /** {@inheritDoc} */
  545.     @Override
  546.     public double[] getMeanElementRate(final SpacecraftState spacecraftState,
  547.                                        final AuxiliaryElements auxiliaryElements, final double[] parameters) {

  549.         // Container of attributes

  551.         final DSSTZonalContext context = initializeStep(auxiliaryElements, parameters);
  552.         // Access to potential U derivatives
  553.         final UAnddU udu = new UAnddU(spacecraftState.getDate(), context, auxiliaryElements, hansen);

  555.         return computeMeanElementRates(spacecraftState.getDate(), context, udu);

  557.     }

  559.     /** {@inheritDoc} */
  560.     @Override
  561.     public <T extends CalculusFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> spacecraftState,
  562.                                                                   final FieldAuxiliaryElements<T> auxiliaryElements,
  563.                                                                   final T[] parameters) {

  565.         // Field used by default
  566.         final Field<T> field = auxiliaryElements.getDate().getField();
  567.         // Container of attributes
  568.         final FieldDSSTZonalContext<T> context = initializeStep(auxiliaryElements, parameters);

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

  573.         // Access to potential U derivatives
  574.         final FieldUAnddU<T> udu = new FieldUAnddU<>(spacecraftState.getDate(), context, auxiliaryElements, fho);

  576.         return computeMeanElementRates(spacecraftState.getDate(), context, udu);

  578.     }

  580.     /** Compute the mean element rates.
  581.      * @param date current date
  582.      * @param context container for attributes
  583.      * @param udu derivatives of the gravitational potential U
  584.      * @return the mean element rates
  585.      */
  586.     private double[] computeMeanElementRates(final AbsoluteDate date,
  587.                                              final DSSTZonalContext context,
  588.                                              final UAnddU udu) {

  590.         // Auxiliary elements related to the current orbit
  591.         final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();

  593.         // Compute cross derivatives [Eq. 2.2-(8)]
  594.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  595.         final double UAlphaGamma   = auxiliaryElements.getAlpha() * udu.getdUdGa() - auxiliaryElements.getGamma() * udu.getdUdAl();
  596.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  597.         final double UBetaGamma    =  auxiliaryElements.getBeta() * udu.getdUdGa() - auxiliaryElements.getGamma() * udu.getdUdBe();
  598.         // Common factor
  599.         final double pUAGmIqUBGoAB = (auxiliaryElements.getP() * UAlphaGamma - I * auxiliaryElements.getQ() * UBetaGamma) * context.getOoAB();

  601.         // Compute mean elements rates [Eq. 3.1-(1)]
  602.         final double da =  0.;
  603.         final double dh =  context.getBoA() * udu.getdUdk() + auxiliaryElements.getK() * pUAGmIqUBGoAB;
  604.         final double dk = -context.getBoA() * udu.getdUdh() - auxiliaryElements.getH() * pUAGmIqUBGoAB;
  605.         final double dp =  context.getMCo2AB() * UBetaGamma;
  606.         final double dq =  context.getMCo2AB() * UAlphaGamma * I;
  607.         final double dM =  context.getM2aoA() * udu.getdUda() + context.getBoABpo() * (auxiliaryElements.getH() * udu.getdUdh() + auxiliaryElements.getK() * udu.getdUdk()) + pUAGmIqUBGoAB;

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

  612.     /** Compute the mean element rates.
  613.      * @param <T> type of the elements
  614.      * @param date current date
  615.      * @param context container for attributes
  616.      * @param udu derivatives of the gravitational potential U
  617.      * @return the mean element rates
  618.      */
  619.     private <T extends CalculusFieldElement<T>> T[] computeMeanElementRates(final FieldAbsoluteDate<T> date,
  620.                                                                         final FieldDSSTZonalContext<T> context,
  621.                                                                         final FieldUAnddU<T> udu) {

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

  626.         // Parameter for array building
  627.         final Field<T> field = date.getField();

  629.         // Compute cross derivatives [Eq. 2.2-(8)]
  630.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  631.         final T UAlphaGamma   = udu.getdUdGa().multiply(auxiliaryElements.getAlpha()).subtract(udu.getdUdAl().multiply(auxiliaryElements.getGamma()));
  632.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  633.         final T UBetaGamma    =  udu.getdUdGa().multiply(auxiliaryElements.getBeta()).subtract(udu.getdUdBe().multiply(auxiliaryElements.getGamma()));
  634.         // Common factor
  635.         final T pUAGmIqUBGoAB = (UAlphaGamma.multiply(auxiliaryElements.getP()).subtract(UBetaGamma.multiply(I).multiply(auxiliaryElements.getQ()))).multiply(context.getOoAB());

  637.         // Compute mean elements rates [Eq. 3.1-(1)]
  638.         final T da =  field.getZero();
  639.         final T dh =  udu.getdUdk().multiply(context.getBoA()).add(pUAGmIqUBGoAB.multiply(auxiliaryElements.getK()));
  640.         final T dk =  (udu.getdUdh().multiply(context.getBoA()).negate()).subtract(pUAGmIqUBGoAB.multiply(auxiliaryElements.getH()));
  641.         final T dp =  UBetaGamma.multiply(context.getMCo2AB());
  642.         final T dq =  UAlphaGamma.multiply(context.getMCo2AB()).multiply(I);
  643.         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()));

  645.         final T[] elements =  MathArrays.buildArray(field, 6);
  646.         elements[0] = da;
  647.         elements[1] = dk;
  648.         elements[2] = dh;
  649.         elements[3] = dq;
  650.         elements[4] = dp;
  651.         elements[5] = dM;

  653.         return elements;
  654.     }

  656.     /** {@inheritDoc} */
  657.     @Override
  658.     public void registerAttitudeProvider(final AttitudeProvider attitudeProvider) {
  659.         //nothing is done since this contribution is not sensitive to attitude
  660.     }

  662.     /** Check if an index is within the accepted interval.
  663.      *
  664.      * @param index the index to check
  665.      * @param lowerBound the lower bound of the interval
  666.      * @param upperBound the upper bound of the interval
  667.      * @return true if the index is between the lower and upper bounds, false otherwise
  668.      */
  669.     private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
  670.         return index >= lowerBound && index <= upperBound;
  671.     }

  673.     /** {@inheritDoc} */
  674.     @Override
  675.     public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {

  677.         final Slot slot = zonalSPCoefs.createSlot(meanStates);
  678.         for (final SpacecraftState meanState : meanStates) {

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

  683.             // Container of attributes
  684.             // Extract the proper parameters valid for the corresponding meanState date from the input array
  685.             final double[] extractedParameters = this.extractParameters(parameters, auxiliaryElements.getDate());
  686.             final DSSTZonalContext context = initializeStep(auxiliaryElements, extractedParameters);

  688.             // Access to potential U derivatives
  689.             final UAnddU udu = new UAnddU(meanState.getDate(), context, auxiliaryElements, hansen);

  691.             // Compute rhoj and sigmaj
  692.             final double[][] rhoSigma = computeRhoSigmaCoefficients(meanState.getDate(), slot, auxiliaryElements);

  694.             // Compute Di
  695.             computeDiCoefficients(meanState.getDate(), slot, context, udu);

  697.             // generate the Cij and Sij coefficients
  698.             final FourierCjSjCoefficients cjsj = new FourierCjSjCoefficients(meanState.getDate(),
  699.                                                                              maxDegreeShortPeriodics,
  700.                                                                              maxEccPowShortPeriodics,
  701.                                                                              maxFrequencyShortPeriodics,
  702.                                                                              context,
  703.                                                                              hansen);

  705.             computeCijSijCoefficients(meanState.getDate(), slot, cjsj, rhoSigma, context, auxiliaryElements, udu);
  706.         }

  708.     }

  710.     /** {@inheritDoc} */
  711.     @Override
  712.     @SuppressWarnings("unchecked")
  713.     public <T extends CalculusFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
  714.                                                                        final FieldSpacecraftState<T>... meanStates) {

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

  719.         final FieldZonalShortPeriodicCoefficients<T> fzspc = (FieldZonalShortPeriodicCoefficients<T>) zonalFieldSPCoefs.get(field);
  720.         final FieldSlot<T> slot = fzspc.createSlot(meanStates);
  721.         for (final FieldSpacecraftState<T> meanState : meanStates) {

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

  726.             // Container of attributes
  727.             // Extract the proper parameters valid for the corresponding meanState date from the input array
  728.             final T[] extractedParameters = this.extractParameters(parameters, auxiliaryElements.getDate());
  729.             final FieldDSSTZonalContext<T> context = initializeStep(auxiliaryElements, extractedParameters);

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

  733.             // Access to potential U derivatives
  734.             final FieldUAnddU<T> udu = new FieldUAnddU<>(meanState.getDate(), context, auxiliaryElements, fho);

  736.             // Compute rhoj and sigmaj
  737.             final T[][] rhoSigma = computeRhoSigmaCoefficients(meanState.getDate(), slot, auxiliaryElements, field);

  739.             // Compute Di
  740.             computeDiCoefficients(meanState.getDate(), slot, context, field, udu);

  742.             // generate the Cij and Sij coefficients
  743.             final FieldFourierCjSjCoefficients<T> cjsj = new FieldFourierCjSjCoefficients<>(meanState.getDate(),
  744.                                                                                             maxDegreeShortPeriodics,
  745.                                                                                             maxEccPowShortPeriodics,
  746.                                                                                             maxFrequencyShortPeriodics,
  747.                                                                                             context,
  748.                                                                                             fho);


  751.             computeCijSijCoefficients(meanState.getDate(), slot, cjsj, rhoSigma, context, auxiliaryElements, field, udu);
  752.         }
  753.     }

  755.     /** {@inheritDoc} */
  756.     public List<ParameterDriver> getParametersDrivers() {
  757.         return Collections.singletonList(gmParameterDriver);
  758.     }

  760.     /** Generate the values for the D<sub>i</sub> coefficients.
  761.      * @param date target date
  762.      * @param slot slot to which the coefficients belong
  763.      * @param context container for attributes
  764.      * @param udu derivatives of the gravitational potential U
  765.      */
  766.     private void computeDiCoefficients(final AbsoluteDate date,
  767.                                        final Slot slot,
  768.                                        final DSSTZonalContext context,
  769.                                        final UAnddU udu) {

  771.         final double[] meanElementRates = computeMeanElementRates(date, context, udu);

  773.         final double[] currentDi = new double[6];

  775.         // Add the coefficients to the interpolation grid
  776.         for (int i = 0; i < 6; i++) {
  777.             currentDi[i] = meanElementRates[i] / context.getMeanMotion();

  779.             if (i == 5) {
  780.                 currentDi[i] += -1.5 * 2 * udu.getU() * context.getOON2A2();
  781.             }

  783.         }

  785.         slot.di.addGridPoint(date, currentDi);

  787.     }

  789.     /** Generate the values for the D<sub>i</sub> coefficients.
  790.      * @param <T> type of the elements
  791.      * @param date target date
  792.      * @param slot slot to which the coefficients belong
  793.      * @param context container for attributes
  794.      * @param field field used by default
  795.      * @param udu derivatives of the gravitational potential U
  796.      */
  797.     private <T extends CalculusFieldElement<T>> void computeDiCoefficients(final FieldAbsoluteDate<T> date,
  798.                                                                        final FieldSlot<T> slot,
  799.                                                                        final FieldDSSTZonalContext<T> context,
  800.                                                                        final Field<T> field,
  801.                                                                        final FieldUAnddU<T> udu) {

  803.         final T[] meanElementRates = computeMeanElementRates(date, context, udu);

  805.         final T[] currentDi = MathArrays.buildArray(field, 6);

  807.         // Add the coefficients to the interpolation grid
  808.         for (int i = 0; i < 6; i++) {
  809.             currentDi[i] = meanElementRates[i].divide(context.getMeanMotion());

  811.             if (i == 5) {
  812.                 currentDi[i] = currentDi[i].add(context.getOON2A2().multiply(udu.getU()).multiply(2.).multiply(-1.5));
  813.             }

  815.         }

  817.         slot.di.addGridPoint(date, currentDi);

  819.     }

  821.     /**
  822.      * Generate the values for the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
  823.      * @param date date of computation
  824.      * @param slot slot to which the coefficients belong
  825.      * @param cjsj Fourier coefficients
  826.      * @param rhoSigma ρ<sub>j</sub> and σ<sub>j</sub>
  827.      * @param context container for attributes
  828.      * @param auxiliaryElements auxiliary elements related to the current orbit
  829.      * @param udu derivatives of the gravitational potential U
  830.      */
  831.     private void computeCijSijCoefficients(final AbsoluteDate date, final Slot slot,
  832.                                            final FourierCjSjCoefficients cjsj,
  833.                                            final double[][] rhoSigma, final DSSTZonalContext context,
  834.                                            final AuxiliaryElements auxiliaryElements,
  835.                                            final UAnddU udu) {

  837.         final int nMax = maxDegreeShortPeriodics;

  839.         // The C<sub>i</sub>⁰ coefficients
  840.         final double[] currentCi0 = new double[] {0., 0., 0., 0., 0., 0.};
  841.         for (int j = 1; j < slot.cij.length; j++) {

  843.             // Create local arrays
  844.             final double[] currentCij = new double[] {0., 0., 0., 0., 0., 0.};
  845.             final double[] currentSij = new double[] {0., 0., 0., 0., 0., 0.};

  847.             // j == 1
  848.             if (j == 1) {
  849.                 final double coef1 = 4 * auxiliaryElements.getK() * udu.getU() - context.getHK() * cjsj.getCj(1) + context.getK2MH2O2() * cjsj.getSj(1);
  850.                 final double coef2 = 4 * auxiliaryElements.getH() * udu.getU() + context.getK2MH2O2() * cjsj.getCj(1) + context.getHK() * cjsj.getSj(1);
  851.                 final double coef3 = (auxiliaryElements.getK() * cjsj.getCj(1) + auxiliaryElements.getH() * cjsj.getSj(1)) / 4.;
  852.                 final double coef4 = (8 * udu.getU() - auxiliaryElements.getH() * cjsj.getCj(1) + auxiliaryElements.getK() * cjsj.getSj(1)) / 4.;

  854.                 //Coefficients for a
  855.                 currentCij[0] += coef1;
  856.                 currentSij[0] += coef2;

  858.                 //Coefficients for k
  859.                 currentCij[1] += coef4;
  860.                 currentSij[1] += coef3;

  862.                 //Coefficients for h
  863.                 currentCij[2] -= coef3;
  864.                 currentSij[2] += coef4;

  866.                 //Coefficients for λ
  867.                 currentCij[5] -= coef2 / 2;
  868.                 currentSij[5] += coef1 / 2;
  869.             }

  871.             // j == 2
  872.             if (j == 2) {
  873.                 final double coef1 = context.getK2MH2() * udu.getU();
  874.                 final double coef2 = 2 * context.getHK() * udu.getU();
  875.                 final double coef3 = auxiliaryElements.getH() * udu.getU() / 2;
  876.                 final double coef4 = auxiliaryElements.getK() * udu.getU() / 2;

  878.                 //Coefficients for a
  879.                 currentCij[0] += coef1;
  880.                 currentSij[0] += coef2;

  882.                 //Coefficients for k
  883.                 currentCij[1] += coef4;
  884.                 currentSij[1] += coef3;

  886.                 //Coefficients for h
  887.                 currentCij[2] -= coef3;
  888.                 currentSij[2] += coef4;

  890.                 //Coefficients for λ
  891.                 currentCij[5] -= coef2 / 2;
  892.                 currentSij[5] += coef1 / 2;
  893.             }

  895.             // j between 1 and 2N-3
  896.             if (isBetween(j, 1, 2 * nMax - 3) && j + 2 < cjsj.jMax) {
  897.                 final double coef1 = ( j + 2 ) * (-context.getHK() * cjsj.getCj(j + 2) + context.getK2MH2O2() * cjsj.getSj(j + 2));
  898.                 final double coef2 = ( j + 2 ) * (context.getK2MH2O2() * cjsj.getCj(j + 2) + context.getHK() * cjsj.getSj(j + 2));
  899.                 final double coef3 = ( j + 2 ) * (auxiliaryElements.getK() * cjsj.getCj(j + 2) + auxiliaryElements.getH() * cjsj.getSj(j + 2)) / 4;
  900.                 final double coef4 = ( j + 2 ) * (auxiliaryElements.getH() * cjsj.getCj(j + 2) - auxiliaryElements.getK() * cjsj.getSj(j + 2)) / 4;

  902.                 //Coefficients for a
  903.                 currentCij[0] += coef1;
  904.                 currentSij[0] -= coef2;

  906.                 //Coefficients for k
  907.                 currentCij[1] += -coef4;
  908.                 currentSij[1] -= coef3;

  910.                 //Coefficients for h
  911.                 currentCij[2] -= coef3;
  912.                 currentSij[2] += coef4;

  914.                 //Coefficients for λ
  915.                 currentCij[5] -= coef2 / 2;
  916.                 currentSij[5] += -coef1 / 2;
  917.             }

  919.             // j between 1 and 2N-2
  920.             if (isBetween(j, 1, 2 * nMax - 2) && j + 1 < cjsj.jMax) {
  921.                 final double coef1 = 2 * ( j + 1 ) * (-auxiliaryElements.getH() * cjsj.getCj(j + 1) + auxiliaryElements.getK() * cjsj.getSj(j + 1));
  922.                 final double coef2 = 2 * ( j + 1 ) * (auxiliaryElements.getK() * cjsj.getCj(j + 1) + auxiliaryElements.getH() * cjsj.getSj(j + 1));
  923.                 final double coef3 = ( j + 1 ) * cjsj.getCj(j + 1);
  924.                 final double coef4 = ( j + 1 ) * cjsj.getSj(j + 1);

  926.                 //Coefficients for a
  927.                 currentCij[0] += coef1;
  928.                 currentSij[0] -= coef2;

  930.                 //Coefficients for k
  931.                 currentCij[1] += coef4;
  932.                 currentSij[1] -= coef3;

  934.                 //Coefficients for h
  935.                 currentCij[2] -= coef3;
  936.                 currentSij[2] -= coef4;

  938.                 //Coefficients for λ
  939.                 currentCij[5] -= coef2 / 2;
  940.                 currentSij[5] += -coef1 / 2;
  941.             }

  943.             // j between 2 and 2N
  944.             if (isBetween(j, 2, 2 * nMax) && j - 1 < cjsj.jMax) {
  945.                 final double coef1 = 2 * ( j - 1 ) * (auxiliaryElements.getH() * cjsj.getCj(j - 1) + auxiliaryElements.getK() * cjsj.getSj(j - 1));
  946.                 final double coef2 = 2 * ( j - 1 ) * (auxiliaryElements.getK() * cjsj.getCj(j - 1) - auxiliaryElements.getH() * cjsj.getSj(j - 1));
  947.                 final double coef3 = ( j - 1 ) * cjsj.getCj(j - 1);
  948.                 final double coef4 = ( j - 1 ) * cjsj.getSj(j - 1);

  950.                 //Coefficients for a
  951.                 currentCij[0] += coef1;
  952.                 currentSij[0] -= coef2;

  954.                 //Coefficients for k
  955.                 currentCij[1] += coef4;
  956.                 currentSij[1] -= coef3;

  958.                 //Coefficients for h
  959.                 currentCij[2] += coef3;
  960.                 currentSij[2] += coef4;

  962.                 //Coefficients for λ
  963.                 currentCij[5] += coef2 / 2;
  964.                 currentSij[5] += coef1 / 2;
  965.             }

  967.             // j between 3 and 2N + 1
  968.             if (isBetween(j, 3, 2 * nMax + 1) && j - 2 < cjsj.jMax) {
  969.                 final double coef1 = ( j - 2 ) * (context.getHK() * cjsj.getCj(j - 2) + context.getK2MH2O2() * cjsj.getSj(j - 2));
  970.                 final double coef2 = ( j - 2 ) * (-context.getK2MH2O2() * cjsj.getCj(j - 2) + context.getHK() * cjsj.getSj(j - 2));
  971.                 final double coef3 = ( j - 2 ) * (auxiliaryElements.getK() * cjsj.getCj(j - 2) - auxiliaryElements.getH() * cjsj.getSj(j - 2)) / 4;
  972.                 final double coef4 = ( j - 2 ) * (auxiliaryElements.getH() * cjsj.getCj(j - 2) + auxiliaryElements.getK() * cjsj.getSj(j - 2)) / 4;
  973.                 final double coef5 = ( j - 2 ) * (context.getK2MH2O2() * cjsj.getCj(j - 2) - context.getHK() * cjsj.getSj(j - 2));

  975.                 //Coefficients for a
  976.                 currentCij[0] += coef1;
  977.                 currentSij[0] += coef2;

  979.                 //Coefficients for k
  980.                 currentCij[1] += coef4;
  981.                 currentSij[1] += -coef3;

  983.                 //Coefficients for h
  984.                 currentCij[2] += coef3;
  985.                 currentSij[2] += coef4;

  987.                 //Coefficients for λ
  988.                 currentCij[5] += coef5 / 2;
  989.                 currentSij[5] += coef1 / 2;
  990.             }

  992.             //multiply by the common factor
  993.             //for a (i == 0) -> χ³ / (n² * a)
  994.             currentCij[0] *= context.getX3ON2A();
  995.             currentSij[0] *= context.getX3ON2A();
  996.             //for k (i == 1) -> χ / (n² * a²)
  997.             currentCij[1] *= context.getXON2A2();
  998.             currentSij[1] *= context.getXON2A2();
  999.             //for h (i == 2) -> χ / (n² * a²)
  1000.             currentCij[2] *= context.getXON2A2();
  1001.             currentSij[2] *= context.getXON2A2();
  1002.             //for λ (i == 5) -> (χ²) / (n² * a² * (χ + 1 ) )
  1003.             currentCij[5] *= context.getX2ON2A2XP1();
  1004.             currentSij[5] *= context.getX2ON2A2XP1();

  1006.             // j is between 1 and 2 * N - 1
  1007.             if (isBetween(j, 1, 2 * nMax - 1) && j < cjsj.jMax) {
  1008.                 // Compute cross derivatives
  1009.                 // Cj(alpha,gamma) = alpha * dC/dgamma - gamma * dC/dalpha
  1010.                 final double CjAlphaGamma   = auxiliaryElements.getAlpha() * cjsj.getdCjdGamma(j) - auxiliaryElements.getGamma() * cjsj.getdCjdAlpha(j);
  1011.                 // Cj(alpha,beta) = alpha * dC/dbeta - beta * dC/dalpha
  1012.                 final double CjAlphaBeta   = auxiliaryElements.getAlpha() * cjsj.getdCjdBeta(j) - auxiliaryElements.getBeta() * cjsj.getdCjdAlpha(j);
  1013.                 // Cj(beta,gamma) = beta * dC/dgamma - gamma * dC/dbeta
  1014.                 final double CjBetaGamma    =  auxiliaryElements.getBeta() * cjsj.getdCjdGamma(j) - auxiliaryElements.getGamma() * cjsj.getdCjdBeta(j);
  1015.                 // Cj(h,k) = h * dC/dk - k * dC/dh
  1016.                 final double CjHK   = auxiliaryElements.getH() * cjsj.getdCjdK(j) - auxiliaryElements.getK() * cjsj.getdCjdH(j);
  1017.                 // Sj(alpha,gamma) = alpha * dS/dgamma - gamma * dS/dalpha
  1018.                 final double SjAlphaGamma   = auxiliaryElements.getAlpha() * cjsj.getdSjdGamma(j) - auxiliaryElements.getGamma() * cjsj.getdSjdAlpha(j);
  1019.                 // Sj(alpha,beta) = alpha * dS/dbeta - beta * dS/dalpha
  1020.                 final double SjAlphaBeta   = auxiliaryElements.getAlpha() * cjsj.getdSjdBeta(j) - auxiliaryElements.getBeta() * cjsj.getdSjdAlpha(j);
  1021.                 // Sj(beta,gamma) = beta * dS/dgamma - gamma * dS/dbeta
  1022.                 final double SjBetaGamma    =  auxiliaryElements.getBeta() * cjsj.getdSjdGamma(j) - auxiliaryElements.getGamma() * cjsj.getdSjdBeta(j);
  1023.                 // Sj(h,k) = h * dS/dk - k * dS/dh
  1024.                 final double SjHK   = auxiliaryElements.getH() * cjsj.getdSjdK(j) - auxiliaryElements.getK() * cjsj.getdSjdH(j);

  1026.                 //Coefficients for a
  1027.                 final double coef1 = context.getX3ON2A() * (3 - context.getBB()) * j;
  1028.                 currentCij[0] += coef1 * cjsj.getSj(j);
  1029.                 currentSij[0] -= coef1 * cjsj.getCj(j);

  1031.                 //Coefficients for k and h
  1032.                 final double coef2 = auxiliaryElements.getP() * CjAlphaGamma - I * auxiliaryElements.getQ() * CjBetaGamma;
  1033.                 final double coef3 = auxiliaryElements.getP() * SjAlphaGamma - I * auxiliaryElements.getQ() * SjBetaGamma;
  1034.                 currentCij[1] -= context.getXON2A2() * (auxiliaryElements.getH() * coef2 + context.getBB() * cjsj.getdCjdH(j) - 1.5 * auxiliaryElements.getK() * j * cjsj.getSj(j));
  1035.                 currentSij[1] -= context.getXON2A2() * (auxiliaryElements.getH() * coef3 + context.getBB() * cjsj.getdSjdH(j) + 1.5 * auxiliaryElements.getK() * j * cjsj.getCj(j));
  1036.                 currentCij[2] += context.getXON2A2() * (auxiliaryElements.getK() * coef2 + context.getBB() * cjsj.getdCjdK(j) + 1.5 * auxiliaryElements.getH() * j * cjsj.getSj(j));
  1037.                 currentSij[2] += context.getXON2A2() * (auxiliaryElements.getK() * coef3 + context.getBB() * cjsj.getdSjdK(j) - 1.5 * auxiliaryElements.getH() * j * cjsj.getCj(j));

  1039.                 //Coefficients for q and p
  1040.                 final double coef4 = CjHK - CjAlphaBeta - j * cjsj.getSj(j);
  1041.                 final double coef5 = SjHK - SjAlphaBeta + j * cjsj.getCj(j);
  1042.                 currentCij[3] = context.getCXO2N2A2() * (-I * CjAlphaGamma + auxiliaryElements.getQ() * coef4);
  1043.                 currentSij[3] = context.getCXO2N2A2() * (-I * SjAlphaGamma + auxiliaryElements.getQ() * coef5);
  1044.                 currentCij[4] = context.getCXO2N2A2() * (-CjBetaGamma + auxiliaryElements.getP() * coef4);
  1045.                 currentSij[4] = context.getCXO2N2A2() * (-SjBetaGamma + auxiliaryElements.getP() * coef5);

  1047.                 //Coefficients for λ
  1048.                 final double coef6 = auxiliaryElements.getH() * cjsj.getdCjdH(j) + auxiliaryElements.getK() * cjsj.getdCjdK(j);
  1049.                 final double coef7 = auxiliaryElements.getH() * cjsj.getdSjdH(j) + auxiliaryElements.getK() * cjsj.getdSjdK(j);
  1050.                 currentCij[5] += context.getOON2A2() * (-2 * auxiliaryElements.getSma() * cjsj.getdCjdA(j) + coef6 / (context.getX() + 1) + context.getX() * coef2 - 3 * cjsj.getCj(j));
  1051.                 currentSij[5] += context.getOON2A2() * (-2 * auxiliaryElements.getSma() * cjsj.getdSjdA(j) + coef7 / (context.getX() + 1) + context.getX() * coef3 - 3 * cjsj.getSj(j));
  1052.             }

  1054.             for (int i = 0; i < 6; i++) {
  1055.                 //Add the current coefficients contribution to C<sub>i</sub>⁰
  1056.                 currentCi0[i] -= currentCij[i] * rhoSigma[j][0] + currentSij[i] * rhoSigma[j][1];
  1057.             }

  1059.             // Add the coefficients to the interpolation grid
  1060.             slot.cij[j].addGridPoint(date, currentCij);
  1061.             slot.sij[j].addGridPoint(date, currentSij);

  1063.         }

  1065.         //Add C<sub>i</sub>⁰ to the interpolation grid
  1066.         slot.cij[0].addGridPoint(date, currentCi0);

  1068.     }

  1070.     /**
  1071.      * Generate the values for the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
  1072.      * @param <T> type of the elements
  1073.      * @param date date of computation
  1074.      * @param slot slot to which the coefficients belong
  1075.      * @param cjsj Fourier coefficients
  1076.      * @param rhoSigma ρ<sub>j</sub> and σ<sub>j</sub>
  1077.      * @param context container for attributes
  1078.      * @param auxiliaryElements auxiliary elements related to the current orbit
  1079.      * @param field field used by default
  1080.      * @param udu derivatives of the gravitational potential U
  1081.      */
  1082.     private <T extends CalculusFieldElement<T>> void computeCijSijCoefficients(final FieldAbsoluteDate<T> date,
  1083.                                                                            final FieldSlot<T> slot,
  1084.                                                                            final FieldFourierCjSjCoefficients<T> cjsj,
  1085.                                                                            final T[][] rhoSigma,
  1086.                                                                            final FieldDSSTZonalContext<T> context,
  1087.                                                                            final FieldAuxiliaryElements<T> auxiliaryElements,
  1088.                                                                            final Field<T> field,
  1089.                                                                            final FieldUAnddU<T> udu) {

  1091.         // Zero
  1092.         final T zero = field.getZero();

  1094.         final int nMax = maxDegreeShortPeriodics;

  1096.         // The C<sub>i</sub>⁰ coefficients
  1097.         final T[] currentCi0 = MathArrays.buildArray(field, 6);
  1098.         Arrays.fill(currentCi0, zero);

  1100.         for (int j = 1; j < slot.cij.length; j++) {

  1102.             // Create local arrays
  1103.             final T[] currentCij = MathArrays.buildArray(field, 6);
  1104.             final T[] currentSij = MathArrays.buildArray(field, 6);

  1106.             Arrays.fill(currentCij, zero);
  1107.             Arrays.fill(currentSij, zero);

  1109.             // j == 1
  1110.             if (j == 1) {
  1111.                 final T coef1 = auxiliaryElements.getK().multiply(udu.getU()).multiply(4.).subtract(context.getHK().multiply(cjsj.getCj(1))).add(context.getK2MH2O2().multiply(cjsj.getSj(1)));
  1112.                 final T coef2 = auxiliaryElements.getH().multiply(udu.getU()).multiply(4.).add(context.getK2MH2O2().multiply(cjsj.getCj(1))).add(context.getHK().multiply(cjsj.getSj(1)));
  1113.                 final T coef3 = auxiliaryElements.getK().multiply(cjsj.getCj(1)).add(auxiliaryElements.getH().multiply(cjsj.getSj(1))).divide(4.);
  1114.                 final T coef4 = udu.getU().multiply(8.).subtract(auxiliaryElements.getH().multiply(cjsj.getCj(1))).add(auxiliaryElements.getK().multiply(cjsj.getSj(1))).divide(4.);

  1116.                 //Coefficients for a
  1117.                 currentCij[0] = currentCij[0].add(coef1);
  1118.                 currentSij[0] = currentSij[0].add(coef2);

  1120.                 //Coefficients for k
  1121.                 currentCij[1] = currentCij[1].add(coef4);
  1122.                 currentSij[1] = currentSij[1].add(coef3);

  1124.                 //Coefficients for h
  1125.                 currentCij[2] = currentCij[2].subtract(coef3);
  1126.                 currentSij[2] = currentSij[2].add(coef4);

  1128.                 //Coefficients for λ
  1129.                 currentCij[5] = currentCij[5].subtract(coef2.divide(2.));
  1130.                 currentSij[5] = currentSij[5].add(coef1.divide(2.));
  1131.             }

  1133.             // j == 2
  1134.             if (j == 2) {
  1135.                 final T coef1 = context.getK2MH2().multiply(udu.getU());
  1136.                 final T coef2 = context.getHK().multiply(udu.getU()).multiply(2.);
  1137.                 final T coef3 = auxiliaryElements.getH().multiply(udu.getU()).divide(2.);
  1138.                 final T coef4 = auxiliaryElements.getK().multiply(udu.getU()).divide(2.);

  1140.                 //Coefficients for a
  1141.                 currentCij[0] = currentCij[0].add(coef1);
  1142.                 currentSij[0] = currentSij[0].add(coef2);

  1144.                 //Coefficients for k
  1145.                 currentCij[1] = currentCij[1].add(coef4);
  1146.                 currentSij[1] = currentSij[1].add(coef3);

  1148.                 //Coefficients for h
  1149.                 currentCij[2] = currentCij[2].subtract(coef3);
  1150.                 currentSij[2] = currentSij[2].add(coef4);

  1152.                 //Coefficients for λ
  1153.                 currentCij[5] = currentCij[5].subtract(coef2.divide(2.));
  1154.                 currentSij[5] = currentSij[5].add(coef1.divide(2.));
  1155.             }

  1157.             // j between 1 and 2N-3
  1158.             if (isBetween(j, 1, 2 * nMax - 3) && j + 2 < cjsj.jMax) {
  1159.                 final T coef1 = context.getHK().negate().multiply(cjsj.getCj(j + 2)).add(context.getK2MH2O2().multiply(cjsj.getSj(j + 2))).multiply(j + 2);
  1160.                 final T coef2 = context.getK2MH2O2().multiply(cjsj.getCj(j + 2)).add(context.getHK().multiply(cjsj.getSj(j + 2))).multiply(j + 2);
  1161.                 final T coef3 = auxiliaryElements.getK().multiply(cjsj.getCj(j + 2)).add(auxiliaryElements.getH().multiply(cjsj.getSj(j + 2))).multiply(j + 2).divide(4.);
  1162.                 final T coef4 = auxiliaryElements.getH().multiply(cjsj.getCj(j + 2)).subtract(auxiliaryElements.getK().multiply(cjsj.getSj(j + 2))).multiply(j + 2).divide(4.);

  1164.                 //Coefficients for a
  1165.                 currentCij[0] = currentCij[0].add(coef1);
  1166.                 currentSij[0] = currentSij[0].subtract(coef2);

  1168.                 //Coefficients for k
  1169.                 currentCij[1] = currentCij[1].add(coef4.negate());
  1170.                 currentSij[1] = currentSij[1].subtract(coef3);

  1172.                 //Coefficients for h
  1173.                 currentCij[2] = currentCij[2].subtract(coef3);
  1174.                 currentSij[2] = currentSij[2].add(coef4);

  1176.                 //Coefficients for λ
  1177.                 currentCij[5] = currentCij[5].subtract(coef2.divide(2.));
  1178.                 currentSij[5] = currentSij[5].add(coef1.negate().divide(2.));
  1179.             }

  1181.             // j between 1 and 2N-2
  1182.             if (isBetween(j, 1, 2 * nMax - 2) && j + 1 < cjsj.jMax) {
  1183.                 final T coef1 = auxiliaryElements.getH().negate().multiply(cjsj.getCj(j + 1)).add(auxiliaryElements.getK().multiply(cjsj.getSj(j + 1))).multiply(2. * (j + 1));
  1184.                 final T coef2 = auxiliaryElements.getK().multiply(cjsj.getCj(j + 1)).add(auxiliaryElements.getH().multiply(cjsj.getSj(j + 1))).multiply(2. * (j + 1));
  1185.                 final T coef3 = cjsj.getCj(j + 1).multiply(j + 1);
  1186.                 final T coef4 = cjsj.getSj(j + 1).multiply(j + 1);

  1188.                 //Coefficients for a
  1189.                 currentCij[0] = currentCij[0].add(coef1);
  1190.                 currentSij[0] = currentSij[0].subtract(coef2);

  1192.                 //Coefficients for k
  1193.                 currentCij[1] = currentCij[1].add(coef4);
  1194.                 currentSij[1] = currentSij[1].subtract(coef3);

  1196.                 //Coefficients for h
  1197.                 currentCij[2] = currentCij[2].subtract(coef3);
  1198.                 currentSij[2] = currentSij[2].subtract(coef4);

  1200.                 //Coefficients for λ
  1201.                 currentCij[5] = currentCij[5].subtract(coef2.divide(2.));
  1202.                 currentSij[5] = currentSij[5].add(coef1.negate().divide(2.));
  1203.             }

  1205.             // j between 2 and 2N
  1206.             if (isBetween(j, 2, 2 * nMax) && j - 1 < cjsj.jMax) {
  1207.                 final T coef1 = auxiliaryElements.getH().multiply(cjsj.getCj(j - 1)).add(auxiliaryElements.getK().multiply(cjsj.getSj(j - 1))).multiply(2 * ( j - 1 ));
  1208.                 final T coef2 = auxiliaryElements.getK().multiply(cjsj.getCj(j - 1)).subtract(auxiliaryElements.getH().multiply(cjsj.getSj(j - 1))).multiply(2 * ( j - 1 ));
  1209.                 final T coef3 = cjsj.getCj(j - 1).multiply(j - 1);
  1210.                 final T coef4 = cjsj.getSj(j - 1).multiply(j - 1);

  1212.                 //Coefficients for a
  1213.                 currentCij[0] = currentCij[0].add(coef1);
  1214.                 currentSij[0] = currentSij[0].subtract(coef2);

  1216.                 //Coefficients for k
  1217.                 currentCij[1] = currentCij[1].add(coef4);
  1218.                 currentSij[1] = currentSij[1].subtract(coef3);

  1220.                 //Coefficients for h
  1221.                 currentCij[2] = currentCij[2].add(coef3);
  1222.                 currentSij[2] = currentSij[2].add(coef4);

  1224.                 //Coefficients for λ
  1225.                 currentCij[5] = currentCij[5].add(coef2.divide(2.));
  1226.                 currentSij[5] = currentSij[5].add(coef1.divide(2.));
  1227.             }

  1229.             // j between 3 and 2N + 1
  1230.             if (isBetween(j, 3, 2 * nMax + 1) && j - 2 < cjsj.jMax) {
  1231.                 final T coef1 = context.getHK().multiply(cjsj.getCj(j - 2)).add(context.getK2MH2O2().multiply(cjsj.getSj(j - 2))).multiply(j - 2);
  1232.                 final T coef2 = context.getK2MH2O2().negate().multiply(cjsj.getCj(j - 2)).add(context.getHK().multiply(cjsj.getSj(j - 2))).multiply(j - 2);
  1233.                 final T coef3 = auxiliaryElements.getK().multiply(cjsj.getCj(j - 2)).subtract(auxiliaryElements.getH().multiply(cjsj.getSj(j - 2))).multiply(j - 2).divide(4.);
  1234.                 final T coef4 = auxiliaryElements.getH().multiply(cjsj.getCj(j - 2)).add(auxiliaryElements.getK().multiply(cjsj.getSj(j - 2))).multiply(j - 2).divide(4.);
  1235.                 final T coef5 = context.getK2MH2O2().multiply(cjsj.getCj(j - 2)).subtract(context.getHK().multiply(cjsj.getSj(j - 2))).multiply(j - 2);

  1237.                 //Coefficients for a
  1238.                 currentCij[0] = currentCij[0].add(coef1);
  1239.                 currentSij[0] = currentSij[0].add(coef2);

  1241.                 //Coefficients for k
  1242.                 currentCij[1] = currentCij[1].add(coef4);
  1243.                 currentSij[1] = currentSij[1].add(coef3.negate());

  1245.                 //Coefficients for h
  1246.                 currentCij[2] = currentCij[2].add(coef3);
  1247.                 currentSij[2] = currentSij[2].add(coef4);

  1249.                 //Coefficients for λ
  1250.                 currentCij[5] = currentCij[5].add(coef5.divide(2.));
  1251.                 currentSij[5] = currentSij[5].add(coef1.divide(2.));
  1252.             }

  1254.             //multiply by the common factor
  1255.             //for a (i == 0) -> χ³ / (n² * a)
  1256.             currentCij[0] = currentCij[0].multiply(context.getX3ON2A());
  1257.             currentSij[0] = currentSij[0].multiply(context.getX3ON2A());
  1258.             //for k (i == 1) -> χ / (n² * a²)
  1259.             currentCij[1] = currentCij[1].multiply(context.getXON2A2());
  1260.             currentSij[1] = currentSij[1].multiply(context.getXON2A2());
  1261.             //for h (i == 2) -> χ / (n² * a²)
  1262.             currentCij[2] = currentCij[2].multiply(context.getXON2A2());
  1263.             currentSij[2] = currentSij[2].multiply(context.getXON2A2());
  1264.             //for λ (i == 5) -> (χ²) / (n² * a² * (χ + 1 ) )
  1265.             currentCij[5] = currentCij[5].multiply(context.getX2ON2A2XP1());
  1266.             currentSij[5] = currentSij[5].multiply(context.getX2ON2A2XP1());

  1268.             // j is between 1 and 2 * N - 1
  1269.             if (isBetween(j, 1, 2 * nMax - 1) && j < cjsj.jMax) {
  1270.                 // Compute cross derivatives
  1271.                 // Cj(alpha,gamma) = alpha * dC/dgamma - gamma * dC/dalpha
  1272.                 final T CjAlphaGamma = auxiliaryElements.getAlpha().multiply(cjsj.getdCjdGamma(j)).subtract(auxiliaryElements.getGamma().multiply(cjsj.getdCjdAlpha(j)));
  1273.                 // Cj(alpha,beta) = alpha * dC/dbeta - beta * dC/dalpha
  1274.                 final T CjAlphaBeta  = auxiliaryElements.getAlpha().multiply(cjsj.getdCjdBeta(j)).subtract(auxiliaryElements.getBeta().multiply(cjsj.getdCjdAlpha(j)));
  1275.                 // Cj(beta,gamma) = beta * dC/dgamma - gamma * dC/dbeta
  1276.                 final T CjBetaGamma  = auxiliaryElements.getBeta().multiply(cjsj.getdCjdGamma(j)).subtract(auxiliaryElements.getGamma().multiply(cjsj.getdCjdBeta(j)));
  1277.                 // Cj(h,k) = h * dC/dk - k * dC/dh
  1278.                 final T CjHK         = auxiliaryElements.getH().multiply(cjsj.getdCjdK(j)).subtract(auxiliaryElements.getK().multiply(cjsj.getdCjdH(j)));
  1279.                 // Sj(alpha,gamma) = alpha * dS/dgamma - gamma * dS/dalpha
  1280.                 final T SjAlphaGamma = auxiliaryElements.getAlpha().multiply(cjsj.getdSjdGamma(j)).subtract(auxiliaryElements.getGamma().multiply(cjsj.getdSjdAlpha(j)));
  1281.                 // Sj(alpha,beta) = alpha * dS/dbeta - beta * dS/dalpha
  1282.                 final T SjAlphaBeta  = auxiliaryElements.getAlpha().multiply(cjsj.getdSjdBeta(j)).subtract(auxiliaryElements.getBeta().multiply(cjsj.getdSjdAlpha(j)));
  1283.                 // Sj(beta,gamma) = beta * dS/dgamma - gamma * dS/dbeta
  1284.                 final T SjBetaGamma  = auxiliaryElements.getBeta().multiply(cjsj.getdSjdGamma(j)).subtract(auxiliaryElements.getGamma().multiply(cjsj.getdSjdBeta(j)));
  1285.                 // Sj(h,k) = h * dS/dk - k * dS/dh
  1286.                 final T SjHK         = auxiliaryElements.getH().multiply(cjsj.getdSjdK(j)).subtract(auxiliaryElements.getK().multiply(cjsj.getdSjdH(j)));

  1288.                 //Coefficients for a
  1289.                 final T coef1 = context.getX3ON2A().multiply(context.getBB().negate().add(3.)).multiply(j);
  1290.                 currentCij[0] = currentCij[0].add(coef1.multiply(cjsj.getSj(j)));
  1291.                 currentSij[0] = currentSij[0].subtract(coef1.multiply(cjsj.getCj(j)));

  1293.                 //Coefficients for k and h
  1294.                 final T coef2 = auxiliaryElements.getP().multiply(CjAlphaGamma).subtract(auxiliaryElements.getQ().multiply(CjBetaGamma).multiply(I));
  1295.                 final T coef3 = auxiliaryElements.getP().multiply(SjAlphaGamma).subtract(auxiliaryElements.getQ().multiply(SjBetaGamma).multiply(I));
  1296.                 currentCij[1] = currentCij[1].subtract(context.getXON2A2().multiply(auxiliaryElements.getH().multiply(coef2).add(context.getBB().multiply(cjsj.getdCjdH(j))).subtract(auxiliaryElements.getK().multiply(1.5).multiply(j).multiply(cjsj.getSj(j)))));
  1297.                 currentSij[1] = currentSij[1].subtract(context.getXON2A2().multiply(auxiliaryElements.getH().multiply(coef3).add(context.getBB().multiply(cjsj.getdSjdH(j))).add(auxiliaryElements.getK().multiply(1.5).multiply(j).multiply(cjsj.getCj(j)))));
  1298.                 currentCij[2] = currentCij[2].add(context.getXON2A2().multiply(auxiliaryElements.getK().multiply(coef2).add(context.getBB().multiply(cjsj.getdCjdK(j))).add(auxiliaryElements.getH().multiply(1.5).multiply(j).multiply(cjsj.getSj(j)))));
  1299.                 currentSij[2] = currentSij[2].add(context.getXON2A2().multiply(auxiliaryElements.getK().multiply(coef3).add(context.getBB().multiply(cjsj.getdSjdK(j))).subtract(auxiliaryElements.getH().multiply(1.5).multiply(j).multiply(cjsj.getCj(j)))));

  1301.                 //Coefficients for q and p
  1302.                 final T coef4 = CjHK.subtract(CjAlphaBeta).subtract(cjsj.getSj(j).multiply(j));
  1303.                 final T coef5 = SjHK.subtract(SjAlphaBeta).add(cjsj.getCj(j).multiply(j));
  1304.                 currentCij[3] = context.getCXO2N2A2().multiply(CjAlphaGamma.multiply(-I).add(auxiliaryElements.getQ().multiply(coef4)));
  1305.                 currentSij[3] = context.getCXO2N2A2().multiply(SjAlphaGamma.multiply(-I).add(auxiliaryElements.getQ().multiply(coef5)));
  1306.                 currentCij[4] = context.getCXO2N2A2().multiply(CjBetaGamma.negate().add(auxiliaryElements.getP().multiply(coef4)));
  1307.                 currentSij[4] = context.getCXO2N2A2().multiply(SjBetaGamma.negate().add(auxiliaryElements.getP().multiply(coef5)));

  1309.                 //Coefficients for λ
  1310.                 final T coef6 = auxiliaryElements.getH().multiply(cjsj.getdCjdH(j)).add(auxiliaryElements.getK().multiply(cjsj.getdCjdK(j)));
  1311.                 final T coef7 = auxiliaryElements.getH().multiply(cjsj.getdSjdH(j)).add(auxiliaryElements.getK().multiply(cjsj.getdSjdK(j)));
  1312.                 currentCij[5] = currentCij[5].add(context.getOON2A2().multiply(auxiliaryElements.getSma().multiply(-2.).multiply(cjsj.getdCjdA(j)).add(coef6.divide(context.getX().add(1.))).add(context.getX().multiply(coef2)).subtract(cjsj.getCj(j).multiply(3.))));
  1313.                 currentSij[5] = currentSij[5].add(context.getOON2A2().multiply(auxiliaryElements.getSma().multiply(-2.).multiply(cjsj.getdSjdA(j)).add(coef7.divide(context.getX().add(1.))).add(context.getX().multiply(coef3)).subtract(cjsj.getSj(j).multiply(3.))));
  1314.             }

  1316.             for (int i = 0; i < 6; i++) {
  1317.                 //Add the current coefficients contribution to C<sub>i</sub>⁰
  1318.                 currentCi0[i] = currentCi0[i].subtract(currentCij[i].multiply(rhoSigma[j][0]).add(currentSij[i].multiply(rhoSigma[j][1])));
  1319.             }

  1321.             // Add the coefficients to the interpolation grid
  1322.             slot.cij[j].addGridPoint(date, currentCij);
  1323.             slot.sij[j].addGridPoint(date, currentSij);

  1325.         }

  1327.         //Add C<sub>i</sub>⁰ to the interpolation grid
  1328.         slot.cij[0].addGridPoint(date, currentCi0);

  1330.     }

  1332.     /**
  1333.      * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
  1334.      * <p>
  1335.      * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
  1336.      *  ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
  1337.      *  σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
  1338.      * </p>
  1339.      * @param date target date
  1340.      * @param slot slot to which the coefficients belong
  1341.      * @param auxiliaryElements auxiliary elements related to the current orbit
  1342.      * @return array containing ρ<sub>j</sub> and σ<sub>j</sub>
  1343.      */
  1344.     private double[][] computeRhoSigmaCoefficients(final AbsoluteDate date,
  1345.                                                    final Slot slot,
  1346.                                                    final AuxiliaryElements auxiliaryElements) {

  1348.         final CjSjCoefficient cjsjKH = new CjSjCoefficient(auxiliaryElements.getK(), auxiliaryElements.getH());
  1349.         final double b = 1. / (1 + auxiliaryElements.getB());

  1351.         // (-b)<sup>j</sup>
  1352.         double mbtj = 1;

  1354.         final double[][] rhoSigma = new double[slot.cij.length][2];
  1355.         for (int j = 1; j < rhoSigma.length; j++) {

  1357.             //Compute current rho and sigma;
  1358.             mbtj *= -b;
  1359.             final double coef  = (1 + j * auxiliaryElements.getB()) * mbtj;
  1360.             final double rho   = coef * cjsjKH.getCj(j);
  1361.             final double sigma = coef * cjsjKH.getSj(j);

  1363.             // Add the coefficients to the interpolation grid
  1364.             rhoSigma[j][0] = rho;
  1365.             rhoSigma[j][1] = sigma;
  1366.         }

  1368.         return rhoSigma;

  1370.     }

  1372.     /**
  1373.      * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
  1374.      * <p>
  1375.      * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
  1376.      *  ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
  1377.      *  σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
  1378.      * </p>
  1379.      * @param <T> type of the elements
  1380.      * @param date target date
  1381.      * @param slot slot to which the coefficients belong
  1382.      * @param auxiliaryElements auxiliary elements related to the current orbit
  1383.      * @param field field used by default
  1384.      * @return array containing ρ<sub>j</sub> and σ<sub>j</sub>
  1385.      */
  1386.     private <T extends CalculusFieldElement<T>> T[][] computeRhoSigmaCoefficients(final FieldAbsoluteDate<T> date,
  1387.                                                                               final FieldSlot<T> slot,
  1388.                                                                               final FieldAuxiliaryElements<T> auxiliaryElements,
  1389.                                                                               final Field<T> field) {
  1390.         final T zero = field.getZero();

  1392.         final FieldCjSjCoefficient<T> cjsjKH = new FieldCjSjCoefficient<>(auxiliaryElements.getK(), auxiliaryElements.getH(), field);
  1393.         final T b = auxiliaryElements.getB().add(1.).reciprocal();

  1395.         // (-b)<sup>j</sup>
  1396.         T mbtj = zero.add(1.);

  1398.         final T[][] rhoSigma = MathArrays.buildArray(field, slot.cij.length, 2);
  1399.         for (int j = 1; j < rhoSigma.length; j++) {

  1401.             //Compute current rho and sigma;
  1402.             mbtj = mbtj.multiply(b.negate());
  1403.             final T coef  = mbtj.multiply(auxiliaryElements.getB().multiply(j).add(1.));
  1404.             final T rho   = coef.multiply(cjsjKH.getCj(j));
  1405.             final T sigma = coef.multiply(cjsjKH.getSj(j));

  1407.             // Add the coefficients to the interpolation grid
  1408.             rhoSigma[j][0] = rho;
  1409.             rhoSigma[j][1] = sigma;
  1410.         }

  1412.         return rhoSigma;

  1414.     }

  1416.     /** The coefficients used to compute the short-periodic zonal contribution.
  1417.      *
  1418.      * <p>
  1419.      * Those coefficients are given in Danielson paper by expressions 4.1-(20) to 4.1.-(25)
  1420.      * </p>
  1421.      * <p>
  1422.      * The coefficients are: <br>
  1423.      * - C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> <br>
  1424.      * - ρ<sub>j</sub> and σ<sub>j</sub> <br>
  1425.      * - C<sub>i</sub>⁰
  1426.      * </p>
  1427.      *
  1428.      * @author Lucian Barbulescu
  1429.      */
  1430.     private static class ZonalShortPeriodicCoefficients implements ShortPeriodTerms {

  1432.         /** Maximum value for j index. */
  1433.         private final int maxFrequencyShortPeriodics;

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

  1438.         /** All coefficients slots. */
  1439.         private final transient TimeSpanMap<Slot> slots;

  1441.         /** Constructor.
  1442.          * @param maxFrequencyShortPeriodics maximum value for j index
  1443.          * @param interpolationPoints number of points used in the interpolation process
  1444.          * @param slots all coefficients slots
  1445.          */
  1446.         ZonalShortPeriodicCoefficients(final int maxFrequencyShortPeriodics, final int interpolationPoints,
  1447.                                        final TimeSpanMap<Slot> slots) {

  1449.             // Save parameters
  1450.             this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
  1451.             this.interpolationPoints        = interpolationPoints;
  1452.             this.slots                      = slots;

  1454.         }

  1456.         /** Get the slot valid for some date.
  1457.          * @param meanStates mean states defining the slot
  1458.          * @return slot valid at the specified date
  1459.          */
  1460.         public Slot createSlot(final SpacecraftState... meanStates) {
  1461.             final Slot         slot  = new Slot(maxFrequencyShortPeriodics, interpolationPoints);
  1462.             final AbsoluteDate first = meanStates[0].getDate();
  1463.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  1464.             final int compare = first.compareTo(last);
  1465.             if (compare < 0) {
  1466.                 slots.addValidAfter(slot, first, false);
  1467.             } else if (compare > 0) {
  1468.                 slots.addValidBefore(slot, first, false);
  1469.             } else {
  1470.                 // single date, valid for all time
  1471.                 slots.addValidAfter(slot, AbsoluteDate.PAST_INFINITY, false);
  1472.             }
  1473.             return slot;
  1474.         }

  1476.         /** {@inheritDoc} */
  1477.         @Override
  1478.         public double[] value(final Orbit meanOrbit) {

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

  1483.             // Get the True longitude L
  1484.             final double L = meanOrbit.getLv();

  1486.             // Compute the center
  1487.             final double center = L - meanOrbit.getLM();

  1489.             // Initialize short periodic variations
  1490.             final double[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  1491.             final double[] d = slot.di.value(meanOrbit.getDate());
  1492.             for (int i = 0; i < 6; i++) {
  1493.                 shortPeriodicVariation[i] +=  center * d[i];
  1494.             }

  1496.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1497.                 final double[] c = slot.cij[j].value(meanOrbit.getDate());
  1498.                 final double[] s = slot.sij[j].value(meanOrbit.getDate());
  1499.                 final SinCos sc  = FastMath.sinCos(j * L);
  1500.                 for (int i = 0; i < 6; i++) {
  1501.                     // add corresponding term to the short periodic variation
  1502.                     shortPeriodicVariation[i] += c[i] * sc.cos();
  1503.                     shortPeriodicVariation[i] += s[i] * sc.sin();
  1504.                 }
  1505.             }

  1507.             return shortPeriodicVariation;
  1508.         }

  1510.         /** {@inheritDoc} */
  1511.         @Override
  1512.         public String getCoefficientsKeyPrefix() {
  1513.             return DSSTZonal.SHORT_PERIOD_PREFIX;
  1514.         }

  1516.         /** {@inheritDoc}
  1517.          * <p>
  1518.          * For zonal terms contributions,there are maxJ cj coefficients,
  1519.          * maxJ sj coefficients and 2 dj coefficients, where maxJ depends
  1520.          * on the orbit. The j index is the integer multiplier for the true
  1521.          * longitude argument in the cj and sj coefficients and the degree
  1522.          * in the polynomial dj coefficients.
  1523.          * </p>
  1524.          */
  1525.         @Override
  1526.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected) {

  1528.             // select the coefficients slot
  1529.             final Slot slot = slots.get(date);

  1531.             final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * maxFrequencyShortPeriodics + 2);
  1532.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  1533.             storeIfSelected(coefficients, selected, slot.di.value(date), "d", 1);
  1534.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1535.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  1536.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  1537.             }
  1538.             return coefficients;

  1540.         }

  1542.         /** Put a coefficient in a map if selected.
  1543.          * @param map map to populate
  1544.          * @param selected set of coefficients that should be put in the map
  1545.          * (empty set means all coefficients are selected)
  1546.          * @param value coefficient value
  1547.          * @param id coefficient identifier
  1548.          * @param indices list of coefficient indices
  1549.          */
  1550.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
  1551.                                      final double[] value, final String id, final int... indices) {
  1552.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  1553.             keyBuilder.append(id);
  1554.             for (int index : indices) {
  1555.                 keyBuilder.append('[').append(index).append(']');
  1556.             }
  1557.             final String key = keyBuilder.toString();
  1558.             if (selected.isEmpty() || selected.contains(key)) {
  1559.                 map.put(key, value);
  1560.             }
  1561.         }

  1563.     }

  1565.     /** The coefficients used to compute the short-periodic zonal contribution.
  1566.     *
  1567.     * <p>
  1568.     * Those coefficients are given in Danielson paper by expressions 4.1-(20) to 4.1.-(25)
  1569.     * </p>
  1570.     * <p>
  1571.     * The coefficients are: <br>
  1572.     * - C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> <br>
  1573.     * - ρ<sub>j</sub> and σ<sub>j</sub> <br>
  1574.     * - C<sub>i</sub>⁰
  1575.     * </p>
  1576.     *
  1577.     * @author Lucian Barbulescu
  1578.     */
  1579.     private static class FieldZonalShortPeriodicCoefficients <T extends CalculusFieldElement<T>> implements FieldShortPeriodTerms<T> {

  1581.         /** Maximum value for j index. */
  1582.         private final int maxFrequencyShortPeriodics;

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

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

  1590.        /** Constructor.
  1591.         * @param maxFrequencyShortPeriodics maximum value for j index
  1592.         * @param interpolationPoints number of points used in the interpolation process
  1593.         * @param slots all coefficients slots
  1594.         */
  1595.         FieldZonalShortPeriodicCoefficients(final int maxFrequencyShortPeriodics, final int interpolationPoints,
  1596.                                             final FieldTimeSpanMap<FieldSlot<T>, T> slots) {

  1598.             // Save parameters
  1599.             this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
  1600.             this.interpolationPoints        = interpolationPoints;
  1601.             this.slots                      = slots;

  1603.         }

  1605.        /** Get the slot valid for some date.
  1606.         * @param meanStates mean states defining the slot
  1607.         * @return slot valid at the specified date
  1608.         */
  1609.         @SuppressWarnings("unchecked")
  1610.         public FieldSlot<T> createSlot(final FieldSpacecraftState<T>... meanStates) {
  1611.             final FieldSlot<T>         slot  = new FieldSlot<>(maxFrequencyShortPeriodics, interpolationPoints);
  1612.             final FieldAbsoluteDate<T> first = meanStates[0].getDate();
  1613.             final FieldAbsoluteDate<T> last  = meanStates[meanStates.length - 1].getDate();
  1614.             if (first.compareTo(last) <= 0) {
  1615.                 slots.addValidAfter(slot, first);
  1616.             } else {
  1617.                 slots.addValidBefore(slot, first);
  1618.             }
  1619.             return slot;
  1620.         }

  1622.         /** {@inheritDoc} */
  1623.         @Override
  1624.         public T[] value(final FieldOrbit<T> meanOrbit) {

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

  1629.             // Get the True longitude L
  1630.             final T L = meanOrbit.getLv();

  1632.             // Compute the center
  1633.             final T center = L.subtract(meanOrbit.getLM());

  1635.             // Initialize short periodic variations
  1636.             final T[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  1637.             final T[] d = slot.di.value(meanOrbit.getDate());
  1638.             for (int i = 0; i < 6; i++) {
  1639.                 shortPeriodicVariation[i] = shortPeriodicVariation[i].add(center.multiply(d[i]));
  1640.             }

  1642.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1643.                 final T[]            c   = slot.cij[j].value(meanOrbit.getDate());
  1644.                 final T[]            s   = slot.sij[j].value(meanOrbit.getDate());
  1645.                 final FieldSinCos<T> sc  = FastMath.sinCos(L.multiply(j));
  1646.                 for (int i = 0; i < 6; i++) {
  1647.                     // add corresponding term to the short periodic variation
  1648.                     shortPeriodicVariation[i] = shortPeriodicVariation[i].add(c[i].multiply(sc.cos()));
  1649.                     shortPeriodicVariation[i] = shortPeriodicVariation[i].add(s[i].multiply(sc.sin()));
  1650.                 }
  1651.             }

  1653.             return shortPeriodicVariation;
  1654.         }

  1656.         /** {@inheritDoc} */
  1657.         @Override
  1658.         public String getCoefficientsKeyPrefix() {
  1659.             return DSSTZonal.SHORT_PERIOD_PREFIX;
  1660.         }

  1662.        /** {@inheritDoc}
  1663.         * <p>
  1664.         * For zonal terms contributions,there are maxJ cj coefficients,
  1665.         * maxJ sj coefficients and 2 dj coefficients, where maxJ depends
  1666.         * on the orbit. The j index is the integer multiplier for the true
  1667.         * longitude argument in the cj and sj coefficients and the degree
  1668.         * in the polynomial dj coefficients.
  1669.         * </p>
  1670.         */
  1671.         @Override
  1672.         public Map<String, T[]> getCoefficients(final FieldAbsoluteDate<T> date, final Set<String> selected) {

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

  1677.             final Map<String, T[]> coefficients = new HashMap<String, T[]>(2 * maxFrequencyShortPeriodics + 2);
  1678.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  1679.             storeIfSelected(coefficients, selected, slot.di.value(date), "d", 1);
  1680.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1681.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  1682.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  1683.             }
  1684.             return coefficients;

  1686.         }

  1688.        /** Put a coefficient in a map if selected.
  1689.         * @param map map to populate
  1690.         * @param selected set of coefficients that should be put in the map
  1691.         * (empty set means all coefficients are selected)
  1692.         * @param value coefficient value
  1693.         * @param id coefficient identifier
  1694.         * @param indices list of coefficient indices
  1695.         */
  1696.         private void storeIfSelected(final Map<String, T[]> map, final Set<String> selected,
  1697.                                      final T[] value, final String id, final int... indices) {
  1698.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  1699.             keyBuilder.append(id);
  1700.             for (int index : indices) {
  1701.                 keyBuilder.append('[').append(index).append(']');
  1702.             }
  1703.             final String key = keyBuilder.toString();
  1704.             if (selected.isEmpty() || selected.contains(key)) {
  1705.                 map.put(key, value);
  1706.             }
  1707.         }

  1709.     }

  1711.     /** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
  1712.      *  <p>
  1713.      *  Those coefficients are given in Danielson paper by expressions 4.1-(13) to 4.1.-(16b)
  1714.      *  </p>
  1715.      */
  1716.     private class FourierCjSjCoefficients {

  1718.         /** The G<sub>js</sub>, H<sub>js</sub>, I<sub>js</sub> and J<sub>js</sub> polynomials. */
  1719.         private final GHIJjsPolynomials ghijCoef;

  1721.         /** L<sub>n</sub><sup>s</sup>(γ). */
  1722.         private final LnsCoefficients lnsCoef;

  1724.         /** Maximum possible value for n. */
  1725.         private final int nMax;

  1727.         /** Maximum possible value for s. */
  1728.         private final int sMax;

  1730.         /** Maximum possible value for j. */
  1731.         private final int jMax;

  1733.         /** The C<sup>j</sup> coefficients and their derivatives.
  1734.          * <p>
  1735.          * Each column of the matrix contains the following values: <br/>
  1736.          * - C<sup>j</sup> <br/>
  1737.          * - dC<sup>j</sup> / da <br/>
  1738.          * - dC<sup>j</sup> / dh <br/>
  1739.          * - dC<sup>j</sup> / dk <br/>
  1740.          * - dC<sup>j</sup> / dα <br/>
  1741.          * - dC<sup>j</sup> / dβ <br/>
  1742.          * - dC<sup>j</sup> / dγ <br/>
  1743.          * </p>
  1744.          */
  1745.         private final double[][] cCoef;

  1747.         /** The S<sup>j</sup> coefficients and their derivatives.
  1748.          * <p>
  1749.          * Each column of the matrix contains the following values: <br/>
  1750.          * - S<sup>j</sup> <br/>
  1751.          * - dS<sup>j</sup> / da <br/>
  1752.          * - dS<sup>j</sup> / dh <br/>
  1753.          * - dS<sup>j</sup> / dk <br/>
  1754.          * - dS<sup>j</sup> / dα <br/>
  1755.          * - dS<sup>j</sup> / dβ <br/>
  1756.          * - dS<sup>j</sup> / dγ <br/>
  1757.          * </p>
  1758.          */
  1759.         private final double[][] sCoef;

  1761.         /** h * &Chi;³. */
  1762.         private final double hXXX;
  1763.         /** k * &Chi;³. */
  1764.         private final double kXXX;

  1766.         /** Create a set of C<sup>j</sup> and the S<sup>j</sup> coefficients.
  1767.          *  @param date the current date
  1768.          *  @param nMax maximum possible value for n
  1769.          *  @param sMax maximum possible value for s
  1770.          *  @param jMax maximum possible value for j
  1771.          *  @param context container for attributes
  1772.          *  @param hansenObjects initialization of hansen objects
  1773.          */
  1774.         FourierCjSjCoefficients(final AbsoluteDate date,
  1775.                                 final int nMax, final int sMax, final int jMax, final DSSTZonalContext context,
  1776.                                 final HansenObjects hansenObjects) {

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

  1780.             this.ghijCoef = new GHIJjsPolynomials(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta());
  1781.             // Qns coefficients
  1782.             final double[][] Qns  = CoefficientsFactory.computeQns(auxiliaryElements.getGamma(), nMax, nMax);

  1784.             this.lnsCoef = new LnsCoefficients(nMax, nMax, Qns, Vns, context.getRoa());
  1785.             this.nMax = nMax;
  1786.             this.sMax = sMax;
  1787.             this.jMax = jMax;

  1789.             // compute the common factors that depends on the mean elements
  1790.             this.hXXX = auxiliaryElements.getH() * context.getXXX();
  1791.             this.kXXX = auxiliaryElements.getK() * context.getXXX();

  1793.             this.cCoef = new double[7][jMax + 1];
  1794.             this.sCoef = new double[7][jMax + 1];

  1796.             for (int s = 0; s <= sMax; s++) {
  1797.                 //Initialise the Hansen roots
  1798.                 hansenObjects.computeHansenObjectsInitValues(context, s);
  1799.             }
  1800.             generateCoefficients(date, context, auxiliaryElements, hansenObjects);
  1801.         }

  1803.         /** Generate all coefficients.
  1804.          * @param date the current date
  1805.          * @param context container for attributes
  1806.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1807.          * @param hansenObjects initialization of hansen objects
  1808.          */
  1809.         private void generateCoefficients(final AbsoluteDate date,
  1810.                                           final DSSTZonalContext context,
  1811.                                           final AuxiliaryElements auxiliaryElements,
  1812.                                           final HansenObjects hansenObjects) {

  1814.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);
  1815.             for (int j = 1; j <= jMax; j++) {

  1817.                 //init arrays
  1818.                 for (int i = 0; i <= 6; i++) {
  1819.                     cCoef[i][j] = 0.;
  1820.                     sCoef[i][j] = 0.;
  1821.                 }

  1823.                 if (isBetween(j, 1, nMax - 1)) {

  1825.                     //compute first double sum where s: j -> N-1 and n: s+1 -> N
  1826.                     for (int s = j; s <= FastMath.min(nMax - 1, sMax); s++) {
  1827.                         // j - s
  1828.                         final int jms = j - s;
  1829.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  1830.                         final int d0smj = (s == j) ? 1 : 2;

  1832.                         for (int n = s + 1; n <= nMax; n++) {
  1833.                             // if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
  1834.                             if ((n + jms) % 2 == 0) {
  1835.                                 // (2 - delta(0,s-j)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
  1836.                                 final double lns = lnsCoef.getLns(n, -jms);
  1837.                                 final double dlns = lnsCoef.getdLnsdGamma(n, -jms);

  1839.                                 final double hjs = ghijCoef.getHjs(s, -jms);
  1840.                                 final double dHjsdh = ghijCoef.getdHjsdh(s, -jms);
  1841.                                 final double dHjsdk = ghijCoef.getdHjsdk(s, -jms);
  1842.                                 final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, -jms);
  1843.                                 final double dHjsdBeta = ghijCoef.getdHjsdBeta(s, -jms);

  1845.                                 final double gjs = ghijCoef.getGjs(s, -jms);
  1846.                                 final double dGjsdh = ghijCoef.getdGjsdh(s, -jms);
  1847.                                 final double dGjsdk = ghijCoef.getdGjsdk(s, -jms);
  1848.                                 final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, -jms);
  1849.                                 final double dGjsdBeta = ghijCoef.getdGjsdBeta(s, -jms);

  1851.                                 // J<sub>n</sub>
  1852.                                 final double jn = -harmonics.getUnnormalizedCnm(n, 0);

  1854.                                 // K₀<sup>-n-1,s</sup>
  1855.                                 final double kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  1856.                                 final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  1858.                                 final double coef0 = d0smj * jn;
  1859.                                 final double coef1 = coef0 * lns;
  1860.                                 final double coef2 = coef1 * kns;
  1861.                                 final double coef3 = coef2 * hjs;
  1862.                                 final double coef4 = coef2 * gjs;

  1864.                                 // Add the term to the coefficients
  1865.                                 cCoef[0][j] += coef3;
  1866.                                 cCoef[1][j] += coef3 * (n + 1);
  1867.                                 cCoef[2][j] += coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
  1868.                                 cCoef[3][j] += coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
  1869.                                 cCoef[4][j] += coef2 * dHjsdAlpha;
  1870.                                 cCoef[5][j] += coef2 * dHjsdBeta;
  1871.                                 cCoef[6][j] += coef0 * dlns * kns * hjs;

  1873.                                 sCoef[0][j] += coef4;
  1874.                                 sCoef[1][j] += coef4 * (n + 1);
  1875.                                 sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
  1876.                                 sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
  1877.                                 sCoef[4][j] += coef2 * dGjsdAlpha;
  1878.                                 sCoef[5][j] += coef2 * dGjsdBeta;
  1879.                                 sCoef[6][j] += coef0 * dlns * kns * gjs;
  1880.                             }
  1881.                         }
  1882.                     }

  1884.                     //compute second double sum where s: 0 -> N-j and n: max(j+s, j+1) -> N
  1885.                     for (int s = 0; s <= FastMath.min(nMax - j, sMax); s++) {
  1886.                         // j + s
  1887.                         final int jps = j + s;
  1888.                         // Kronecker symbols (2 - delta(0,j+s))
  1889.                         final double d0spj = (s == -j) ? 1 : 2;

  1891.                         for (int n = FastMath.max(j + s, j + 1); n <= nMax; n++) {
  1892.                             // if n + (j+s) is odd, then the term is equal to zero due to the factor Vn,s+j
  1893.                             if ((n + jps) % 2 == 0) {
  1894.                                 // (2 - delta(0,s+j)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j+s</sup>
  1895.                                 final double lns = lnsCoef.getLns(n, jps);
  1896.                                 final double dlns = lnsCoef.getdLnsdGamma(n, jps);

  1898.                                 final double hjs = ghijCoef.getHjs(s, jps);
  1899.                                 final double dHjsdh = ghijCoef.getdHjsdh(s, jps);
  1900.                                 final double dHjsdk = ghijCoef.getdHjsdk(s, jps);
  1901.                                 final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, jps);
  1902.                                 final double dHjsdBeta = ghijCoef.getdHjsdBeta(s, jps);

  1904.                                 final double gjs = ghijCoef.getGjs(s, jps);
  1905.                                 final double dGjsdh = ghijCoef.getdGjsdh(s, jps);
  1906.                                 final double dGjsdk = ghijCoef.getdGjsdk(s, jps);
  1907.                                 final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, jps);
  1908.                                 final double dGjsdBeta = ghijCoef.getdGjsdBeta(s, jps);

  1910.                                 // J<sub>n</sub>
  1911.                                 final double jn = -harmonics.getUnnormalizedCnm(n, 0);

  1913.                                 // K₀<sup>-n-1,s</sup>
  1914.                                 final double kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  1915.                                 final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  1917.                                 final double coef0 = d0spj * jn;
  1918.                                 final double coef1 = coef0 * lns;
  1919.                                 final double coef2 = coef1 * kns;

  1921.                                 final double coef3 = coef2 * hjs;
  1922.                                 final double coef4 = coef2 * gjs;

  1924.                                 // Add the term to the coefficients
  1925.                                 cCoef[0][j] -= coef3;
  1926.                                 cCoef[1][j] -= coef3 * (n + 1);
  1927.                                 cCoef[2][j] -= coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
  1928.                                 cCoef[3][j] -= coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
  1929.                                 cCoef[4][j] -= coef2 * dHjsdAlpha;
  1930.                                 cCoef[5][j] -= coef2 * dHjsdBeta;
  1931.                                 cCoef[6][j] -= coef0 * dlns * kns * hjs;

  1933.                                 sCoef[0][j] += coef4;
  1934.                                 sCoef[1][j] += coef4 * (n + 1);
  1935.                                 sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
  1936.                                 sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
  1937.                                 sCoef[4][j] += coef2 * dGjsdAlpha;
  1938.                                 sCoef[5][j] += coef2 * dGjsdBeta;
  1939.                                 sCoef[6][j] += coef0 * dlns * kns * gjs;
  1940.                             }
  1941.                         }
  1942.                     }

  1944.                     //compute third double sum where s: 1 -> j and  n: j+1 -> N
  1945.                     for (int s = 1; s <= FastMath.min(j, sMax); s++) {
  1946.                         // j - s
  1947.                         final int jms = j - s;
  1948.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  1949.                         final int d0smj = (s == j) ? 1 : 2;

  1951.                         for (int n = j + 1; n <= nMax; n++) {
  1952.                             // if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
  1953.                             if ((n + jms) % 2 == 0) {
  1954.                                 // (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
  1955.                                 final double lns = lnsCoef.getLns(n, jms);
  1956.                                 final double dlns = lnsCoef.getdLnsdGamma(n, jms);

  1958.                                 final double ijs = ghijCoef.getIjs(s, jms);
  1959.                                 final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  1960.                                 final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  1961.                                 final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  1962.                                 final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  1964.                                 final double jjs = ghijCoef.getJjs(s, jms);
  1965.                                 final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  1966.                                 final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  1967.                                 final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  1968.                                 final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

  1970.                                 // J<sub>n</sub>
  1971.                                 final double jn = -harmonics.getUnnormalizedCnm(n, 0);

  1973.                                 // K₀<sup>-n-1,s</sup>
  1974.                                 final double kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  1975.                                 final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  1977.                                 final double coef0 = d0smj * jn;
  1978.                                 final double coef1 = coef0 * lns;
  1979.                                 final double coef2 = coef1 * kns;

  1981.                                 final double coef3 = coef2 * ijs;
  1982.                                 final double coef4 = coef2 * jjs;

  1984.                                 // Add the term to the coefficients
  1985.                                 cCoef[0][j] -= coef3;
  1986.                                 cCoef[1][j] -= coef3 * (n + 1);
  1987.                                 cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  1988.                                 cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  1989.                                 cCoef[4][j] -= coef2 * dIjsdAlpha;
  1990.                                 cCoef[5][j] -= coef2 * dIjsdBeta;
  1991.                                 cCoef[6][j] -= coef0 * dlns * kns * ijs;

  1993.                                 sCoef[0][j] += coef4;
  1994.                                 sCoef[1][j] += coef4 * (n + 1);
  1995.                                 sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  1996.                                 sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  1997.                                 sCoef[4][j] += coef2 * dJjsdAlpha;
  1998.                                 sCoef[5][j] += coef2 * dJjsdBeta;
  1999.                                 sCoef[6][j] += coef0 * dlns * kns * jjs;
  2000.                             }
  2001.                         }
  2002.                     }
  2003.                 }

  2005.                 if (isBetween(j, 2, nMax)) {
  2006.                     //add first term
  2007.                     // J<sub>j</sub>
  2008.                     final double jj = -harmonics.getUnnormalizedCnm(j, 0);
  2009.                     double kns = hansenObjects.getHansenObjects()[0].getValue(-j - 1, context.getX());
  2010.                     double dkns = hansenObjects.getHansenObjects()[0].getDerivative(-j - 1, context.getX());

  2012.                     double lns = lnsCoef.getLns(j, j);
  2013.                     //dlns is 0 because n == s == j

  2015.                     final double hjs = ghijCoef.getHjs(0, j);
  2016.                     final double dHjsdh = ghijCoef.getdHjsdh(0, j);
  2017.                     final double dHjsdk = ghijCoef.getdHjsdk(0, j);
  2018.                     final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(0, j);
  2019.                     final double dHjsdBeta = ghijCoef.getdHjsdBeta(0, j);

  2021.                     final double gjs = ghijCoef.getGjs(0, j);
  2022.                     final double dGjsdh = ghijCoef.getdGjsdh(0, j);
  2023.                     final double dGjsdk = ghijCoef.getdGjsdk(0, j);
  2024.                     final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(0, j);
  2025.                     final double dGjsdBeta = ghijCoef.getdGjsdBeta(0, j);

  2027.                     // 2 * J<sub>j</sub> * K₀<sup>-j-1,0</sup> * L<sub>j</sub><sup>j</sup>
  2028.                     double coef0 = 2 * jj;
  2029.                     double coef1 = coef0 * lns;
  2030.                     double coef2 = coef1 * kns;

  2032.                     double coef3 = coef2 * hjs;
  2033.                     double coef4 = coef2 * gjs;

  2035.                     // Add the term to the coefficients
  2036.                     cCoef[0][j] -= coef3;
  2037.                     cCoef[1][j] -= coef3 * (j + 1);
  2038.                     cCoef[2][j] -= coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
  2039.                     cCoef[3][j] -= coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
  2040.                     cCoef[4][j] -= coef2 * dHjsdAlpha;
  2041.                     cCoef[5][j] -= coef2 * dHjsdBeta;
  2042.                     //no contribution to cCoef[6][j] because dlns is 0

  2044.                     sCoef[0][j] += coef4;
  2045.                     sCoef[1][j] += coef4 * (j + 1);
  2046.                     sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
  2047.                     sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
  2048.                     sCoef[4][j] += coef2 * dGjsdAlpha;
  2049.                     sCoef[5][j] += coef2 * dGjsdBeta;
  2050.                     //no contribution to sCoef[6][j] because dlns is 0

  2052.                     //compute simple sum where s: 1 -> j-1
  2053.                     for (int s = 1; s <= FastMath.min(j - 1, sMax); s++) {
  2054.                         // j - s
  2055.                         final int jms = j - s;
  2056.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2057.                         final int d0smj = (s == j) ? 1 : 2;

  2059.                         // if s is odd, then the term is equal to zero due to the factor Vj,s-j
  2060.                         if (s % 2 == 0) {
  2061.                             // (2 - delta(0,j-s)) * J<sub>j</sub> * K₀<sup>-j-1,s</sup> * L<sub>j</sub><sup>j-s</sup>
  2062.                             kns   = hansenObjects.getHansenObjects()[s].getValue(-j - 1, context.getX());
  2063.                             dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-j - 1, context.getX());

  2065.                             lns = lnsCoef.getLns(j, jms);
  2066.                             final double dlns = lnsCoef.getdLnsdGamma(j, jms);

  2068.                             final double ijs = ghijCoef.getIjs(s, jms);
  2069.                             final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2070.                             final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2071.                             final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2072.                             final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2074.                             final double jjs = ghijCoef.getJjs(s, jms);
  2075.                             final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2076.                             final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2077.                             final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2078.                             final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

  2080.                             coef0 = d0smj * jj;
  2081.                             coef1 = coef0 * lns;
  2082.                             coef2 = coef1 * kns;

  2084.                             coef3 = coef2 * ijs;
  2085.                             coef4 = coef2 * jjs;

  2087.                             // Add the term to the coefficients
  2088.                             cCoef[0][j] -= coef3;
  2089.                             cCoef[1][j] -= coef3 * (j + 1);
  2090.                             cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2091.                             cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2092.                             cCoef[4][j] -= coef2 * dIjsdAlpha;
  2093.                             cCoef[5][j] -= coef2 * dIjsdBeta;
  2094.                             cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2096.                             sCoef[0][j] += coef4;
  2097.                             sCoef[1][j] += coef4 * (j + 1);
  2098.                             sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2099.                             sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2100.                             sCoef[4][j] += coef2 * dJjsdAlpha;
  2101.                             sCoef[5][j] += coef2 * dJjsdBeta;
  2102.                             sCoef[6][j] += coef0 * dlns * kns * jjs;
  2103.                         }
  2104.                     }
  2105.                 }

  2107.                 if (isBetween(j, 3, 2 * nMax - 1)) {
  2108.                     //compute uppercase sigma expressions

  2110.                     //min(j-1,N)
  2111.                     final int minjm1on = FastMath.min(j - 1, nMax);

  2113.                     //if j is even
  2114.                     if (j % 2 == 0) {
  2115.                         //compute first double sum where s: j-min(j-1,N) -> j/2-1 and n: j-s -> min(j-1,N)
  2116.                         for (int s = j - minjm1on; s <= FastMath.min(j / 2 - 1, sMax); s++) {
  2117.                             // j - s
  2118.                             final int jms = j - s;
  2119.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2120.                             final int d0smj = (s == j) ? 1 : 2;

  2122.                             for (int n = j - s; n <= minjm1on; n++) {
  2123.                                 // if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
  2124.                                 if ((n + jms) % 2 == 0) {
  2125.                                     // (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
  2126.                                     final double lns = lnsCoef.getLns(n, jms);
  2127.                                     final double dlns = lnsCoef.getdLnsdGamma(n, jms);

  2129.                                     final double ijs = ghijCoef.getIjs(s, jms);
  2130.                                     final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2131.                                     final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2132.                                     final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2133.                                     final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2135.                                     final double jjs = ghijCoef.getJjs(s, jms);
  2136.                                     final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2137.                                     final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2138.                                     final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2139.                                     final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

  2141.                                     // J<sub>n</sub>
  2142.                                     final double jn = -harmonics.getUnnormalizedCnm(n, 0);

  2144.                                     // K₀<sup>-n-1,s</sup>
  2145.                                     final double kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  2146.                                     final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  2148.                                     final double coef0 = d0smj * jn;
  2149.                                     final double coef1 = coef0 * lns;
  2150.                                     final double coef2 = coef1 * kns;

  2152.                                     final double coef3 = coef2 * ijs;
  2153.                                     final double coef4 = coef2 * jjs;

  2155.                                     // Add the term to the coefficients
  2156.                                     cCoef[0][j] -= coef3;
  2157.                                     cCoef[1][j] -= coef3 * (n + 1);
  2158.                                     cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2159.                                     cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2160.                                     cCoef[4][j] -= coef2 * dIjsdAlpha;
  2161.                                     cCoef[5][j] -= coef2 * dIjsdBeta;
  2162.                                     cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2164.                                     sCoef[0][j] += coef4;
  2165.                                     sCoef[1][j] += coef4 * (n + 1);
  2166.                                     sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2167.                                     sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2168.                                     sCoef[4][j] += coef2 * dJjsdAlpha;
  2169.                                     sCoef[5][j] += coef2 * dJjsdBeta;
  2170.                                     sCoef[6][j] += coef0 * dlns * kns * jjs;
  2171.                                 }
  2172.                             }
  2173.                         }

  2175.                         //compute second double sum where s: j/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  2176.                         for (int s = j / 2; s <=  FastMath.min(minjm1on - 1, sMax); s++) {
  2177.                             // j - s
  2178.                             final int jms = j - s;
  2179.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2180.                             final int d0smj = (s == j) ? 1 : 2;

  2182.                             for (int n = s + 1; n <= minjm1on; n++) {
  2183.                                 // if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
  2184.                                 if ((n + jms) % 2 == 0) {
  2185.                                     // (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
  2186.                                     final double lns = lnsCoef.getLns(n, jms);
  2187.                                     final double dlns = lnsCoef.getdLnsdGamma(n, jms);

  2189.                                     final double ijs = ghijCoef.getIjs(s, jms);
  2190.                                     final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2191.                                     final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2192.                                     final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2193.                                     final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2195.                                     final double jjs = ghijCoef.getJjs(s, jms);
  2196.                                     final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2197.                                     final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2198.                                     final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2199.                                     final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

  2201.                                     // J<sub>n</sub>
  2202.                                     final double jn = -harmonics.getUnnormalizedCnm(n, 0);

  2204.                                     // K₀<sup>-n-1,s</sup>
  2205.                                     final double kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  2206.                                     final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  2208.                                     final double coef0 = d0smj * jn;
  2209.                                     final double coef1 = coef0 * lns;
  2210.                                     final double coef2 = coef1 * kns;

  2212.                                     final double coef3 = coef2 * ijs;
  2213.                                     final double coef4 = coef2 * jjs;

  2215.                                     // Add the term to the coefficients
  2216.                                     cCoef[0][j] -= coef3;
  2217.                                     cCoef[1][j] -= coef3 * (n + 1);
  2218.                                     cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2219.                                     cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2220.                                     cCoef[4][j] -= coef2 * dIjsdAlpha;
  2221.                                     cCoef[5][j] -= coef2 * dIjsdBeta;
  2222.                                     cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2224.                                     sCoef[0][j] += coef4;
  2225.                                     sCoef[1][j] += coef4 * (n + 1);
  2226.                                     sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2227.                                     sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2228.                                     sCoef[4][j] += coef2 * dJjsdAlpha;
  2229.                                     sCoef[5][j] += coef2 * dJjsdBeta;
  2230.                                     sCoef[6][j] += coef0 * dlns * kns * jjs;
  2231.                                 }
  2232.                             }
  2233.                         }
  2234.                     }

  2236.                     //if j is odd
  2237.                     else {
  2238.                         //compute first double sum where s: (j-1)/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  2239.                         for (int s = (j - 1) / 2; s <= FastMath.min(minjm1on - 1, sMax); s++) {
  2240.                             // j - s
  2241.                             final int jms = j - s;
  2242.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2243.                             final int d0smj = (s == j) ? 1 : 2;

  2245.                             for (int n = s + 1; n <= minjm1on; n++) {
  2246.                                 // if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
  2247.                                 if ((n + jms) % 2 == 0) {
  2248.                                     // (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
  2249.                                     final double lns = lnsCoef.getLns(n, jms);
  2250.                                     final double dlns = lnsCoef.getdLnsdGamma(n, jms);

  2252.                                     final double ijs = ghijCoef.getIjs(s, jms);
  2253.                                     final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2254.                                     final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2255.                                     final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2256.                                     final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2258.                                     final double jjs = ghijCoef.getJjs(s, jms);
  2259.                                     final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2260.                                     final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2261.                                     final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2262.                                     final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

  2264.                                     // J<sub>n</sub>
  2265.                                     final double jn = -harmonics.getUnnormalizedCnm(n, 0);

  2267.                                     // K₀<sup>-n-1,s</sup>

  2269.                                     final double kns = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  2270.                                     final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  2272.                                     final double coef0 = d0smj * jn;
  2273.                                     final double coef1 = coef0 * lns;
  2274.                                     final double coef2 = coef1 * kns;

  2276.                                     final double coef3 = coef2 * ijs;
  2277.                                     final double coef4 = coef2 * jjs;

  2279.                                     // Add the term to the coefficients
  2280.                                     cCoef[0][j] -= coef3;
  2281.                                     cCoef[1][j] -= coef3 * (n + 1);
  2282.                                     cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2283.                                     cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2284.                                     cCoef[4][j] -= coef2 * dIjsdAlpha;
  2285.                                     cCoef[5][j] -= coef2 * dIjsdBeta;
  2286.                                     cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2288.                                     sCoef[0][j] += coef4;
  2289.                                     sCoef[1][j] += coef4 * (n + 1);
  2290.                                     sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2291.                                     sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2292.                                     sCoef[4][j] += coef2 * dJjsdAlpha;
  2293.                                     sCoef[5][j] += coef2 * dJjsdBeta;
  2294.                                     sCoef[6][j] += coef0 * dlns * kns * jjs;
  2295.                                 }
  2296.                             }
  2297.                         }

  2299.                         //the second double sum is added only if N >= 4 and j between 5 and 2*N-3
  2300.                         if (nMax >= 4 && isBetween(j, 5, 2 * nMax - 3)) {
  2301.                             //compute second double sum where s: j-min(j-1,N) -> (j-3)/2 and n: j-s -> min(j-1,N)
  2302.                             for (int s = j - minjm1on; s <= FastMath.min((j - 3) / 2, sMax); s++) {
  2303.                                 // j - s
  2304.                                 final int jms = j - s;
  2305.                                 // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2306.                                 final int d0smj = (s == j) ? 1 : 2;

  2308.                                 for (int n = j - s; n <= minjm1on; n++) {
  2309.                                     // if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
  2310.                                     if ((n + jms) % 2 == 0) {
  2311.                                         // (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
  2312.                                         final double lns = lnsCoef.getLns(n, jms);
  2313.                                         final double dlns = lnsCoef.getdLnsdGamma(n, jms);

  2315.                                         final double ijs = ghijCoef.getIjs(s, jms);
  2316.                                         final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2317.                                         final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2318.                                         final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2319.                                         final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2321.                                         final double jjs = ghijCoef.getJjs(s, jms);
  2322.                                         final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2323.                                         final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2324.                                         final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2325.                                         final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

  2327.                                         // J<sub>n</sub>
  2328.                                         final double jn = -harmonics.getUnnormalizedCnm(n, 0);

  2330.                                         // K₀<sup>-n-1,s</sup>
  2331.                                         final double kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  2332.                                         final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  2334.                                         final double coef0 = d0smj * jn;
  2335.                                         final double coef1 = coef0 * lns;
  2336.                                         final double coef2 = coef1 * kns;

  2338.                                         final double coef3 = coef2 * ijs;
  2339.                                         final double coef4 = coef2 * jjs;

  2341.                                         // Add the term to the coefficients
  2342.                                         cCoef[0][j] -= coef3;
  2343.                                         cCoef[1][j] -= coef3 * (n + 1);
  2344.                                         cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2345.                                         cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2346.                                         cCoef[4][j] -= coef2 * dIjsdAlpha;
  2347.                                         cCoef[5][j] -= coef2 * dIjsdBeta;
  2348.                                         cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2350.                                         sCoef[0][j] += coef4;
  2351.                                         sCoef[1][j] += coef4 * (n + 1);
  2352.                                         sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2353.                                         sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2354.                                         sCoef[4][j] += coef2 * dJjsdAlpha;
  2355.                                         sCoef[5][j] += coef2 * dJjsdBeta;
  2356.                                         sCoef[6][j] += coef0 * dlns * kns * jjs;
  2357.                                     }
  2358.                                 }
  2359.                             }
  2360.                         }
  2361.                     }
  2362.                 }

  2364.                 cCoef[0][j] *= -context.getMuoa() / j;
  2365.                 cCoef[1][j] *=  context.getMuoa() / ( j * auxiliaryElements.getSma() );
  2366.                 cCoef[2][j] *= -context.getMuoa() / j;
  2367.                 cCoef[3][j] *= -context.getMuoa() / j;
  2368.                 cCoef[4][j] *= -context.getMuoa() / j;
  2369.                 cCoef[5][j] *= -context.getMuoa() / j;
  2370.                 cCoef[6][j] *= -context.getMuoa() / j;

  2372.                 sCoef[0][j] *= -context.getMuoa() / j;
  2373.                 sCoef[1][j] *=  context.getMuoa() / ( j * auxiliaryElements.getSma() );
  2374.                 sCoef[2][j] *= -context.getMuoa() / j;
  2375.                 sCoef[3][j] *= -context.getMuoa() / j;
  2376.                 sCoef[4][j] *= -context.getMuoa() / j;
  2377.                 sCoef[5][j] *= -context.getMuoa() / j;
  2378.                 sCoef[6][j] *= -context.getMuoa() / j;

  2380.             }
  2381.         }

  2383.         /** Check if an index is within the accepted interval.
  2384.          *
  2385.          * @param index the index to check
  2386.          * @param lowerBound the lower bound of the interval
  2387.          * @param upperBound the upper bound of the interval
  2388.          * @return true if the index is between the lower and upper bounds, false otherwise
  2389.          */
  2390.         private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
  2391.             return index >= lowerBound && index <= upperBound;
  2392.         }

  2394.         /**Get the value of C<sup>j</sup>.
  2395.          *
  2396.          * @param j j index
  2397.          * @return C<sup>j</sup>
  2398.          */
  2399.         public double getCj(final int j) {
  2400.             return cCoef[0][j];
  2401.         }

  2403.         /**Get the value of dC<sup>j</sup> / da.
  2404.          *
  2405.          * @param j j index
  2406.          * @return dC<sup>j</sup> / da
  2407.          */
  2408.         public double getdCjdA(final int j) {
  2409.             return cCoef[1][j];
  2410.         }

  2412.         /**Get the value of dC<sup>j</sup> / dh.
  2413.          *
  2414.          * @param j j index
  2415.          * @return dC<sup>j</sup> / dh
  2416.          */
  2417.         public double getdCjdH(final int j) {
  2418.             return cCoef[2][j];
  2419.         }

  2421.         /**Get the value of dC<sup>j</sup> / dk.
  2422.          *
  2423.          * @param j j index
  2424.          * @return dC<sup>j</sup> / dk
  2425.          */
  2426.         public double getdCjdK(final int j) {
  2427.             return cCoef[3][j];
  2428.         }

  2430.         /**Get the value of dC<sup>j</sup> / dα.
  2431.          *
  2432.          * @param j j index
  2433.          * @return dC<sup>j</sup> / dα
  2434.          */
  2435.         public double getdCjdAlpha(final int j) {
  2436.             return cCoef[4][j];
  2437.         }

  2439.         /**Get the value of dC<sup>j</sup> / dβ.
  2440.          *
  2441.          * @param j j index
  2442.          * @return dC<sup>j</sup> / dβ
  2443.          */
  2444.         public double getdCjdBeta(final int j) {
  2445.             return cCoef[5][j];
  2446.         }

  2448.         /**Get the value of dC<sup>j</sup> / dγ.
  2449.          *
  2450.          * @param j j index
  2451.          * @return dC<sup>j</sup> / dγ
  2452.          */
  2453.         public double getdCjdGamma(final int j) {
  2454.             return cCoef[6][j];
  2455.         }

  2457.         /**Get the value of S<sup>j</sup>.
  2458.          *
  2459.          * @param j j index
  2460.          * @return S<sup>j</sup>
  2461.          */
  2462.         public double getSj(final int j) {
  2463.             return sCoef[0][j];
  2464.         }

  2466.         /**Get the value of dS<sup>j</sup> / da.
  2467.          *
  2468.          * @param j j index
  2469.          * @return dS<sup>j</sup> / da
  2470.          */
  2471.         public double getdSjdA(final int j) {
  2472.             return sCoef[1][j];
  2473.         }

  2475.         /**Get the value of dS<sup>j</sup> / dh.
  2476.          *
  2477.          * @param j j index
  2478.          * @return dS<sup>j</sup> / dh
  2479.          */
  2480.         public double getdSjdH(final int j) {
  2481.             return sCoef[2][j];
  2482.         }

  2484.         /**Get the value of dS<sup>j</sup> / dk.
  2485.          *
  2486.          * @param j j index
  2487.          * @return dS<sup>j</sup> / dk
  2488.          */
  2489.         public double getdSjdK(final int j) {
  2490.             return sCoef[3][j];
  2491.         }

  2493.         /**Get the value of dS<sup>j</sup> / dα.
  2494.          *
  2495.          * @param j j index
  2496.          * @return dS<sup>j</sup> / dα
  2497.          */
  2498.         public double getdSjdAlpha(final int j) {
  2499.             return sCoef[4][j];
  2500.         }

  2502.         /**Get the value of dS<sup>j</sup> / dβ.
  2503.          *
  2504.          * @param j j index
  2505.          * @return dS<sup>j</sup> / dβ
  2506.          */
  2507.         public double getdSjdBeta(final int j) {
  2508.             return sCoef[5][j];
  2509.         }

  2511.         /**Get the value of dS<sup>j</sup> /  dγ.
  2512.          *
  2513.          * @param j j index
  2514.          * @return dS<sup>j</sup> /  dγ
  2515.          */
  2516.         public double getdSjdGamma(final int j) {
  2517.             return sCoef[6][j];
  2518.         }
  2519.     }

  2521.     /** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
  2522.      *  <p>
  2523.      *  Those coefficients are given in Danielson paper by expressions 4.1-(13) to 4.1.-(16b)
  2524.      *  </p>
  2525.      */
  2526.     private class FieldFourierCjSjCoefficients <T extends CalculusFieldElement<T>> {

  2528.         /** The G<sub>js</sub>, H<sub>js</sub>, I<sub>js</sub> and J<sub>js</sub> polynomials. */
  2529.         private final FieldGHIJjsPolynomials<T> ghijCoef;

  2531.         /** L<sub>n</sub><sup>s</sup>(γ). */
  2532.         private final FieldLnsCoefficients<T> lnsCoef;

  2534.         /** Maximum possible value for n. */
  2535.         private final int nMax;

  2537.         /** Maximum possible value for s. */
  2538.         private final int sMax;

  2540.         /** Maximum possible value for j. */
  2541.         private final int jMax;

  2543.         /** The C<sup>j</sup> coefficients and their derivatives.
  2544.          * <p>
  2545.          * Each column of the matrix contains the following values: <br/>
  2546.          * - C<sup>j</sup> <br/>
  2547.          * - dC<sup>j</sup> / da <br/>
  2548.          * - dC<sup>j</sup> / dh <br/>
  2549.          * - dC<sup>j</sup> / dk <br/>
  2550.          * - dC<sup>j</sup> / dα <br/>
  2551.          * - dC<sup>j</sup> / dβ <br/>
  2552.          * - dC<sup>j</sup> / dγ <br/>
  2553.          * </p>
  2554.          */
  2555.         private final T[][] cCoef;

  2557.         /** The S<sup>j</sup> coefficients and their derivatives.
  2558.          * <p>
  2559.          * Each column of the matrix contains the following values: <br/>
  2560.          * - S<sup>j</sup> <br/>
  2561.          * - dS<sup>j</sup> / da <br/>
  2562.          * - dS<sup>j</sup> / dh <br/>
  2563.          * - dS<sup>j</sup> / dk <br/>
  2564.          * - dS<sup>j</sup> / dα <br/>
  2565.          * - dS<sup>j</sup> / dβ <br/>
  2566.          * - dS<sup>j</sup> / dγ <br/>
  2567.          * </p>
  2568.          */
  2569.         private final T[][] sCoef;

  2571.         /** h * &Chi;³. */
  2572.         private final T hXXX;
  2573.         /** k * &Chi;³. */
  2574.         private final T kXXX;

  2576.         /** Create a set of C<sup>j</sup> and the S<sup>j</sup> coefficients.
  2577.          *  @param date the current date
  2578.          *  @param nMax maximum possible value for n
  2579.          *  @param sMax maximum possible value for s
  2580.          *  @param jMax maximum possible value for j
  2581.          *  @param context container for attributes
  2582.          *  @param hansenObjects initialization of hansen objects
  2583.          */
  2584.         FieldFourierCjSjCoefficients(final FieldAbsoluteDate<T> date,
  2585.                                      final int nMax, final int sMax, final int jMax,
  2586.                                      final FieldDSSTZonalContext<T> context,
  2587.                                      final FieldHansenObjects<T> hansenObjects) {

  2589.             //Field used by default
  2590.             final Field<T> field = date.getField();

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

  2594.             this.ghijCoef = new FieldGHIJjsPolynomials<>(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta());
  2595.             // Qns coefficients
  2596.             final T[][] Qns = CoefficientsFactory.computeQns(auxiliaryElements.getGamma(), nMax, nMax);

  2598.             this.lnsCoef = new FieldLnsCoefficients<>(nMax, nMax, Qns, Vns, context.getRoa(), field);
  2599.             this.nMax = nMax;
  2600.             this.sMax = sMax;
  2601.             this.jMax = jMax;

  2603.             // compute the common factors that depends on the mean elements
  2604.             this.hXXX = auxiliaryElements.getH().multiply(context.getXXX());
  2605.             this.kXXX = auxiliaryElements.getK().multiply(context.getXXX());

  2607.             this.cCoef = MathArrays.buildArray(field, 7, jMax + 1);
  2608.             this.sCoef = MathArrays.buildArray(field, 7, jMax + 1);

  2610.             for (int s = 0; s <= sMax; s++) {
  2611.                 //Initialise the Hansen roots
  2612.                 hansenObjects.computeHansenObjectsInitValues(context, s);
  2613.             }
  2614.             generateCoefficients(date, context, auxiliaryElements, hansenObjects, field);
  2615.         }

  2617.         /** Generate all coefficients.
  2618.          * @param date the current date
  2619.          * @param context container for attributes
  2620.          * @param hansenObjects initialization of hansen objects
  2621.          * @param auxiliaryElements auxiliary elements related to the current orbit
  2622.          * @param field field used by default
  2623.          */
  2624.         private void generateCoefficients(final FieldAbsoluteDate<T> date,
  2625.                                           final FieldDSSTZonalContext<T> context,
  2626.                                           final FieldAuxiliaryElements<T> auxiliaryElements,
  2627.                                           final FieldHansenObjects<T> hansenObjects,
  2628.                                           final Field<T> field) {

  2630.             //Zero
  2631.             final T zero = field.getZero();

  2633.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date.toAbsoluteDate());
  2634.             for (int j = 1; j <= jMax; j++) {

  2636.                 //init arrays
  2637.                 for (int i = 0; i <= 6; i++) {
  2638.                     cCoef[i][j] = zero;
  2639.                     sCoef[i][j] = zero;
  2640.                 }

  2642.                 if (isBetween(j, 1, nMax - 1)) {

  2644.                     //compute first double sum where s: j -> N-1 and n: s+1 -> N
  2645.                     for (int s = j; s <= FastMath.min(nMax - 1, sMax); s++) {
  2646.                         // j - s
  2647.                         final int jms = j - s;
  2648.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2649.                         final int d0smj = (s == j) ? 1 : 2;

  2651.                         for (int n = s + 1; n <= nMax; n++) {
  2652.                             // if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
  2653.                             if ((n + jms) % 2 == 0) {
  2654.                                 // (2 - delta(0,s-j)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
  2655.                                 final T lns  = lnsCoef.getLns(n, -jms);
  2656.                                 final T dlns = lnsCoef.getdLnsdGamma(n, -jms);

  2658.                                 final T hjs        = ghijCoef.getHjs(s, -jms);
  2659.                                 final T dHjsdh     = ghijCoef.getdHjsdh(s, -jms);
  2660.                                 final T dHjsdk     = ghijCoef.getdHjsdk(s, -jms);
  2661.                                 final T dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, -jms);
  2662.                                 final T dHjsdBeta  = ghijCoef.getdHjsdBeta(s, -jms);

  2664.                                 final T gjs        = ghijCoef.getGjs(s, -jms);
  2665.                                 final T dGjsdh     = ghijCoef.getdGjsdh(s, -jms);
  2666.                                 final T dGjsdk     = ghijCoef.getdGjsdk(s, -jms);
  2667.                                 final T dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, -jms);
  2668.                                 final T dGjsdBeta  = ghijCoef.getdGjsdBeta(s, -jms);

  2670.                                 // J<sub>n</sub>
  2671.                                 final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));

  2673.                                 // K₀<sup>-n-1,s</sup>
  2674.                                 final T kns   = (T) hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  2675.                                 final T dkns  = (T) hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  2677.                                 final T coef0 = jn.multiply(d0smj);
  2678.                                 final T coef1 = coef0.multiply(lns);
  2679.                                 final T coef2 = coef1.multiply(kns);
  2680.                                 final T coef3 = coef2.multiply(hjs);
  2681.                                 final T coef4 = coef2.multiply(gjs);

  2683.                                 // Add the term to the coefficients
  2684.                                 cCoef[0][j] = cCoef[0][j].add(coef3);
  2685.                                 cCoef[1][j] = cCoef[1][j].add(coef3.multiply(n + 1));
  2686.                                 cCoef[2][j] = cCoef[2][j].add(coef1.multiply(kns.multiply(dHjsdh).add(hjs.multiply(hXXX).multiply(dkns))));
  2687.                                 cCoef[3][j] = cCoef[3][j].add(coef1.multiply(kns.multiply(dHjsdk).add(hjs.multiply(kXXX).multiply(dkns))));
  2688.                                 cCoef[4][j] = cCoef[4][j].add(coef2.multiply(dHjsdAlpha));
  2689.                                 cCoef[5][j] = cCoef[5][j].add(coef2.multiply(dHjsdBeta));
  2690.                                 cCoef[6][j] = cCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(hjs));

  2692.                                 sCoef[0][j] = sCoef[0][j].add(coef4);
  2693.                                 sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  2694.                                 sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dGjsdh).add(gjs.multiply(hXXX).multiply(dkns))));
  2695.                                 sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dGjsdk).add(gjs.multiply(kXXX).multiply(dkns))));
  2696.                                 sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dGjsdAlpha));
  2697.                                 sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dGjsdBeta));
  2698.                                 sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(gjs));
  2699.                             }
  2700.                         }
  2701.                     }

  2703.                     //compute second double sum where s: 0 -> N-j and n: max(j+s, j+1) -> N
  2704.                     for (int s = 0; s <= FastMath.min(nMax - j, sMax); s++) {
  2705.                         // j + s
  2706.                         final int jps = j + s;
  2707.                         // Kronecker symbols (2 - delta(0,j+s))
  2708.                         final double d0spj = (s == -j) ? 1 : 2;

  2710.                         for (int n = FastMath.max(j + s, j + 1); n <= nMax; n++) {
  2711.                             // if n + (j+s) is odd, then the term is equal to zero due to the factor Vn,s+j
  2712.                             if ((n + jps) % 2 == 0) {
  2713.                                 // (2 - delta(0,s+j)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j+s</sup>
  2714.                                 final T lns  = lnsCoef.getLns(n, jps);
  2715.                                 final T dlns = lnsCoef.getdLnsdGamma(n, jps);

  2717.                                 final T hjs        = ghijCoef.getHjs(s, jps);
  2718.                                 final T dHjsdh     = ghijCoef.getdHjsdh(s, jps);
  2719.                                 final T dHjsdk     = ghijCoef.getdHjsdk(s, jps);
  2720.                                 final T dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, jps);
  2721.                                 final T dHjsdBeta  = ghijCoef.getdHjsdBeta(s, jps);

  2723.                                 final T gjs        = ghijCoef.getGjs(s, jps);
  2724.                                 final T dGjsdh     = ghijCoef.getdGjsdh(s, jps);
  2725.                                 final T dGjsdk     = ghijCoef.getdGjsdk(s, jps);
  2726.                                 final T dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, jps);
  2727.                                 final T dGjsdBeta  = ghijCoef.getdGjsdBeta(s, jps);

  2729.                                 // J<sub>n</sub>
  2730.                                 final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));

  2732.                                 // K₀<sup>-n-1,s</sup>
  2733.                                 final T kns   = (T) hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  2734.                                 final T dkns  = (T) hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  2736.                                 final T coef0 = jn.multiply(d0spj);
  2737.                                 final T coef1 = coef0.multiply(lns);
  2738.                                 final T coef2 = coef1.multiply(kns);

  2740.                                 final T coef3 = coef2.multiply(hjs);
  2741.                                 final T coef4 = coef2.multiply(gjs);

  2743.                                 // Add the term to the coefficients
  2744.                                 cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2745.                                 cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  2746.                                 cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dHjsdh).add(hjs.multiply(hXXX).multiply(dkns))));
  2747.                                 cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dHjsdk).add(hjs.multiply(kXXX).multiply(dkns))));
  2748.                                 cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dHjsdAlpha));
  2749.                                 cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dHjsdBeta));
  2750.                                 cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(hjs));

  2752.                                 sCoef[0][j] = sCoef[0][j].add(coef4);
  2753.                                 sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  2754.                                 sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dGjsdh).add(gjs.multiply(hXXX).multiply(dkns))));
  2755.                                 sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dGjsdk).add(gjs.multiply(kXXX).multiply(dkns))));
  2756.                                 sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dGjsdAlpha));
  2757.                                 sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dGjsdBeta));
  2758.                                 sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(gjs));
  2759.                             }
  2760.                         }
  2761.                     }

  2763.                     //compute third double sum where s: 1 -> j and  n: j+1 -> N
  2764.                     for (int s = 1; s <= FastMath.min(j, sMax); s++) {
  2765.                         // j - s
  2766.                         final int jms = j - s;
  2767.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2768.                         final int d0smj = (s == j) ? 1 : 2;

  2770.                         for (int n = j + 1; n <= nMax; n++) {
  2771.                             // if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
  2772.                             if ((n + jms) % 2 == 0) {
  2773.                                 // (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
  2774.                                 final T lns  = lnsCoef.getLns(n, jms);
  2775.                                 final T dlns = lnsCoef.getdLnsdGamma(n, jms);

  2777.                                 final T ijs        = ghijCoef.getIjs(s, jms);
  2778.                                 final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  2779.                                 final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  2780.                                 final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2781.                                 final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  2783.                                 final T jjs        = ghijCoef.getJjs(s, jms);
  2784.                                 final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  2785.                                 final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  2786.                                 final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2787.                                 final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

  2789.                                 // J<sub>n</sub>
  2790.                                 final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));

  2792.                                 // K₀<sup>-n-1,s</sup>
  2793.                                 final T kns   = (T) hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  2794.                                 final T dkns  = (T) hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  2796.                                 final T coef0 = jn.multiply(d0smj);
  2797.                                 final T coef1 = coef0.multiply(lns);
  2798.                                 final T coef2 = coef1.multiply(kns);

  2800.                                 final T coef3 = coef2.multiply(ijs);
  2801.                                 final T coef4 = coef2.multiply(jjs);

  2803.                                 // Add the term to the coefficients
  2804.                                 cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2805.                                 cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  2806.                                 cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  2807.                                 cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  2808.                                 cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  2809.                                 cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  2810.                                 cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  2812.                                 sCoef[0][j] = sCoef[0][j].add(coef4);
  2813.                                 sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  2814.                                 sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  2815.                                 sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  2816.                                 sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  2817.                                 sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  2818.                                 sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  2819.                             }
  2820.                         }
  2821.                     }
  2822.                 }

  2824.                 if (isBetween(j, 2, nMax)) {
  2825.                     //add first term
  2826.                     // J<sub>j</sub>
  2827.                     final T jj = zero.subtract(harmonics.getUnnormalizedCnm(j, 0));
  2828.                     T kns  = (T) hansenObjects.getHansenObjects()[0].getValue(-j - 1, context.getX());
  2829.                     T dkns = (T) hansenObjects.getHansenObjects()[0].getDerivative(-j - 1, context.getX());

  2831.                     T lns = lnsCoef.getLns(j, j);
  2832.                     //dlns is 0 because n == s == j

  2834.                     final T hjs        = ghijCoef.getHjs(0, j);
  2835.                     final T dHjsdh     = ghijCoef.getdHjsdh(0, j);
  2836.                     final T dHjsdk     = ghijCoef.getdHjsdk(0, j);
  2837.                     final T dHjsdAlpha = ghijCoef.getdHjsdAlpha(0, j);
  2838.                     final T dHjsdBeta  = ghijCoef.getdHjsdBeta(0, j);

  2840.                     final T gjs        = ghijCoef.getGjs(0, j);
  2841.                     final T dGjsdh     = ghijCoef.getdGjsdh(0, j);
  2842.                     final T dGjsdk     = ghijCoef.getdGjsdk(0, j);
  2843.                     final T dGjsdAlpha = ghijCoef.getdGjsdAlpha(0, j);
  2844.                     final T dGjsdBeta  = ghijCoef.getdGjsdBeta(0, j);

  2846.                     // 2 * J<sub>j</sub> * K₀<sup>-j-1,0</sup> * L<sub>j</sub><sup>j</sup>
  2847.                     T coef0 = jj.multiply(2.);
  2848.                     T coef1 = coef0.multiply(lns);
  2849.                     T coef2 = coef1.multiply(kns);

  2851.                     T coef3 = coef2.multiply(hjs);
  2852.                     T coef4 = coef2.multiply(gjs);

  2854.                     // Add the term to the coefficients
  2855.                     cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2856.                     cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(j + 1));
  2857.                     cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dHjsdh).add(hjs.multiply(hXXX).multiply(dkns))));
  2858.                     cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dHjsdk).add(hjs.multiply(kXXX).multiply(dkns))));
  2859.                     cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dHjsdAlpha));
  2860.                     cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dHjsdBeta));
  2861.                     //no contribution to cCoef[6][j] because dlns is 0

  2863.                     sCoef[0][j] = sCoef[0][j].add(coef4);
  2864.                     sCoef[1][j] = sCoef[1][j].add(coef4.multiply(j + 1));
  2865.                     sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dGjsdh).add(gjs.multiply(hXXX).multiply(dkns))));
  2866.                     sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dGjsdk).add(gjs.multiply(kXXX).multiply(dkns))));
  2867.                     sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dGjsdAlpha));
  2868.                     sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dGjsdBeta));
  2869.                     //no contribution to sCoef[6][j] because dlns is 0

  2871.                     //compute simple sum where s: 1 -> j-1
  2872.                     for (int s = 1; s <= FastMath.min(j - 1, sMax); s++) {
  2873.                         // j - s
  2874.                         final int jms = j - s;
  2875.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2876.                         final int d0smj = (s == j) ? 1 : 2;

  2878.                         // if s is odd, then the term is equal to zero due to the factor Vj,s-j
  2879.                         if (s % 2 == 0) {
  2880.                             // (2 - delta(0,j-s)) * J<sub>j</sub> * K₀<sup>-j-1,s</sup> * L<sub>j</sub><sup>j-s</sup>
  2881.                             kns   = (T) hansenObjects.getHansenObjects()[s].getValue(-j - 1, context.getX());
  2882.                             dkns  = (T) hansenObjects.getHansenObjects()[s].getDerivative(-j - 1, context.getX());

  2884.                             lns = lnsCoef.getLns(j, jms);
  2885.                             final T dlns = lnsCoef.getdLnsdGamma(j, jms);

  2887.                             final T ijs        = ghijCoef.getIjs(s, jms);
  2888.                             final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  2889.                             final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  2890.                             final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2891.                             final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  2893.                             final T jjs        = ghijCoef.getJjs(s, jms);
  2894.                             final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  2895.                             final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  2896.                             final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2897.                             final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

  2899.                             coef0 = jj.multiply(d0smj);
  2900.                             coef1 = coef0.multiply(lns);
  2901.                             coef2 = coef1.multiply(kns);

  2903.                             coef3 = coef2.multiply(ijs);
  2904.                             coef4 = coef2.multiply(jjs);

  2906.                             // Add the term to the coefficients
  2907.                             cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2908.                             cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(j + 1));
  2909.                             cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  2910.                             cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  2911.                             cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  2912.                             cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  2913.                             cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  2915.                             sCoef[0][j] = sCoef[0][j].add(coef4);
  2916.                             sCoef[1][j] = sCoef[1][j].add(coef4.multiply(j + 1));
  2917.                             sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  2918.                             sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  2919.                             sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  2920.                             sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  2921.                             sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  2922.                         }
  2923.                     }
  2924.                 }

  2926.                 if (isBetween(j, 3, 2 * nMax - 1)) {
  2927.                     //compute uppercase sigma expressions

  2929.                     //min(j-1,N)
  2930.                     final int minjm1on = FastMath.min(j - 1, nMax);

  2932.                     //if j is even
  2933.                     if (j % 2 == 0) {
  2934.                         //compute first double sum where s: j-min(j-1,N) -> j/2-1 and n: j-s -> min(j-1,N)
  2935.                         for (int s = j - minjm1on; s <= FastMath.min(j / 2 - 1, sMax); s++) {
  2936.                             // j - s
  2937.                             final int jms = j - s;
  2938.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2939.                             final int d0smj = (s == j) ? 1 : 2;

  2941.                             for (int n = j - s; n <= minjm1on; n++) {
  2942.                                 // if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
  2943.                                 if ((n + jms) % 2 == 0) {
  2944.                                     // (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
  2945.                                     final T lns  = lnsCoef.getLns(n, jms);
  2946.                                     final T dlns = lnsCoef.getdLnsdGamma(n, jms);

  2948.                                     final T ijs        = ghijCoef.getIjs(s, jms);
  2949.                                     final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  2950.                                     final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  2951.                                     final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2952.                                     final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  2954.                                     final T jjs        = ghijCoef.getJjs(s, jms);
  2955.                                     final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  2956.                                     final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  2957.                                     final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2958.                                     final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

  2960.                                     // J<sub>n</sub>
  2961.                                     final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));

  2963.                                     // K₀<sup>-n-1,s</sup>
  2964.                                     final T kns   = (T) hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  2965.                                     final T dkns  = (T) hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  2967.                                     final T coef0 = jn.multiply(d0smj);
  2968.                                     final T coef1 = coef0.multiply(lns);
  2969.                                     final T coef2 = coef1.multiply(kns);

  2971.                                     final T coef3 = coef2.multiply(ijs);
  2972.                                     final T coef4 = coef2.multiply(jjs);

  2974.                                     // Add the term to the coefficients
  2975.                                     cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2976.                                     cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  2977.                                     cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  2978.                                     cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  2979.                                     cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  2980.                                     cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  2981.                                     cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  2983.                                     sCoef[0][j] = sCoef[0][j].add(coef4);
  2984.                                     sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  2985.                                     sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  2986.                                     sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  2987.                                     sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  2988.                                     sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  2989.                                     sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  2990.                                 }
  2991.                             }
  2992.                         }

  2994.                         //compute second double sum where s: j/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  2995.                         for (int s = j / 2; s <=  FastMath.min(minjm1on - 1, sMax); s++) {
  2996.                             // j - s
  2997.                             final int jms = j - s;
  2998.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2999.                             final int d0smj = (s == j) ? 1 : 2;

  3001.                             for (int n = s + 1; n <= minjm1on; n++) {
  3002.                                 // if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
  3003.                                 if ((n + jms) % 2 == 0) {
  3004.                                     // (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
  3005.                                     final T lns  = lnsCoef.getLns(n, jms);
  3006.                                     final T dlns = lnsCoef.getdLnsdGamma(n, jms);

  3008.                                     final T ijs        = ghijCoef.getIjs(s, jms);
  3009.                                     final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  3010.                                     final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  3011.                                     final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  3012.                                     final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  3014.                                     final T jjs        = ghijCoef.getJjs(s, jms);
  3015.                                     final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  3016.                                     final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  3017.                                     final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  3018.                                     final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

  3020.                                     // J<sub>n</sub>
  3021.                                     final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));

  3023.                                     // K₀<sup>-n-1,s</sup>
  3024.                                     final T kns   = (T) hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  3025.                                     final T dkns  = (T) hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  3027.                                     final T coef0 = jn.multiply(d0smj);
  3028.                                     final T coef1 = coef0.multiply(lns);
  3029.                                     final T coef2 = coef1.multiply(kns);

  3031.                                     final T coef3 = coef2.multiply(ijs);
  3032.                                     final T coef4 = coef2.multiply(jjs);

  3034.                                     // Add the term to the coefficients
  3035.                                     cCoef[0][j] = cCoef[0][j].subtract(coef3);
  3036.                                     cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  3037.                                     cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  3038.                                     cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  3039.                                     cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  3040.                                     cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  3041.                                     cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  3043.                                     sCoef[0][j] = sCoef[0][j].add(coef4);
  3044.                                     sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  3045.                                     sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  3046.                                     sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  3047.                                     sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  3048.                                     sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  3049.                                     sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  3050.                                 }
  3051.                             }
  3052.                         }
  3053.                     }

  3055.                     //if j is odd
  3056.                     else {
  3057.                         //compute first double sum where s: (j-1)/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  3058.                         for (int s = (j - 1) / 2; s <= FastMath.min(minjm1on - 1, sMax); s++) {
  3059.                             // j - s
  3060.                             final int jms = j - s;
  3061.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  3062.                             final int d0smj = (s == j) ? 1 : 2;

  3064.                             for (int n = s + 1; n <= minjm1on; n++) {
  3065.                                 // if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
  3066.                                 if ((n + jms) % 2 == 0) {
  3067.                                     // (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
  3068.                                     final T lns  = lnsCoef.getLns(n, jms);
  3069.                                     final T dlns = lnsCoef.getdLnsdGamma(n, jms);

  3071.                                     final T ijs        = ghijCoef.getIjs(s, jms);
  3072.                                     final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  3073.                                     final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  3074.                                     final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  3075.                                     final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  3077.                                     final T jjs        = ghijCoef.getJjs(s, jms);
  3078.                                     final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  3079.                                     final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  3080.                                     final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  3081.                                     final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

  3083.                                     // J<sub>n</sub>
  3084.                                     final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));

  3086.                                     // K₀<sup>-n-1,s</sup>

  3088.                                     final T kns   = (T) hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  3089.                                     final T dkns  = (T) hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  3091.                                     final T coef0 = jn.multiply(d0smj);
  3092.                                     final T coef1 = coef0.multiply(lns);
  3093.                                     final T coef2 = coef1.multiply(kns);

  3095.                                     final T coef3 = coef2.multiply(ijs);
  3096.                                     final T coef4 = coef2.multiply(jjs);

  3098.                                     // Add the term to the coefficients
  3099.                                     cCoef[0][j] = cCoef[0][j].subtract(coef3);
  3100.                                     cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  3101.                                     cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  3102.                                     cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  3103.                                     cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  3104.                                     cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  3105.                                     cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  3107.                                     sCoef[0][j] = sCoef[0][j].add(coef4);
  3108.                                     sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  3109.                                     sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  3110.                                     sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  3111.                                     sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  3112.                                     sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  3113.                                     sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  3114.                                 }
  3115.                             }
  3116.                         }

  3118.                         //the second double sum is added only if N >= 4 and j between 5 and 2*N-3
  3119.                         if (nMax >= 4 && isBetween(j, 5, 2 * nMax - 3)) {
  3120.                             //compute second double sum where s: j-min(j-1,N) -> (j-3)/2 and n: j-s -> min(j-1,N)
  3121.                             for (int s = j - minjm1on; s <= FastMath.min((j - 3) / 2, sMax); s++) {
  3122.                                 // j - s
  3123.                                 final int jms = j - s;
  3124.                                 // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  3125.                                 final int d0smj = (s == j) ? 1 : 2;

  3127.                                 for (int n = j - s; n <= minjm1on; n++) {
  3128.                                     // if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
  3129.                                     if ((n + jms) % 2 == 0) {
  3130.                                         // (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
  3131.                                         final T lns  = lnsCoef.getLns(n, jms);
  3132.                                         final T dlns = lnsCoef.getdLnsdGamma(n, jms);

  3134.                                         final T ijs        = ghijCoef.getIjs(s, jms);
  3135.                                         final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  3136.                                         final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  3137.                                         final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  3138.                                         final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  3140.                                         final T jjs        = ghijCoef.getJjs(s, jms);
  3141.                                         final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  3142.                                         final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  3143.                                         final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  3144.                                         final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

  3146.                                         // J<sub>n</sub>
  3147.                                         final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));

  3149.                                         // K₀<sup>-n-1,s</sup>
  3150.                                         final T kns   = (T) hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getX());
  3151.                                         final T dkns  = (T) hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  3153.                                         final T coef0 = jn.multiply(d0smj);
  3154.                                         final T coef1 = coef0.multiply(lns);
  3155.                                         final T coef2 = coef1.multiply(kns);

  3157.                                         final T coef3 = coef2.multiply(ijs);
  3158.                                         final T coef4 = coef2.multiply(jjs);

  3160.                                         // Add the term to the coefficients
  3161.                                         cCoef[0][j] = cCoef[0][j].subtract(coef3);
  3162.                                         cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  3163.                                         cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  3164.                                         cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  3165.                                         cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  3166.                                         cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  3167.                                         cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  3169.                                         sCoef[0][j] = sCoef[0][j].add(coef4);
  3170.                                         sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  3171.                                         sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  3172.                                         sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  3173.                                         sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  3174.                                         sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  3175.                                         sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  3176.                                     }
  3177.                                 }
  3178.                             }
  3179.                         }
  3180.                     }
  3181.                 }

  3183.                 cCoef[0][j] = cCoef[0][j].multiply(context.getMuoa().divide(j).negate());
  3184.                 cCoef[1][j] = cCoef[1][j].multiply(context.getMuoa().divide(auxiliaryElements.getSma().multiply(j)));
  3185.                 cCoef[2][j] = cCoef[2][j].multiply(context.getMuoa().divide(j).negate());
  3186.                 cCoef[3][j] = cCoef[3][j].multiply(context.getMuoa().divide(j).negate());
  3187.                 cCoef[4][j] = cCoef[4][j].multiply(context.getMuoa().divide(j).negate());
  3188.                 cCoef[5][j] = cCoef[5][j].multiply(context.getMuoa().divide(j).negate());
  3189.                 cCoef[6][j] = cCoef[6][j].multiply(context.getMuoa().divide(j).negate());

  3191.                 sCoef[0][j] = sCoef[0][j].multiply(context.getMuoa().divide(j).negate());
  3192.                 sCoef[1][j] = sCoef[1][j].multiply(context.getMuoa().divide(auxiliaryElements.getSma().multiply(j)));
  3193.                 sCoef[2][j] = sCoef[2][j].multiply(context.getMuoa().divide(j).negate());
  3194.                 sCoef[3][j] = sCoef[3][j].multiply(context.getMuoa().divide(j).negate());
  3195.                 sCoef[4][j] = sCoef[4][j].multiply(context.getMuoa().divide(j).negate());
  3196.                 sCoef[5][j] = sCoef[5][j].multiply(context.getMuoa().divide(j).negate());
  3197.                 sCoef[6][j] = sCoef[6][j].multiply(context.getMuoa().divide(j).negate());

  3199.             }
  3200.         }

  3202.         /** Check if an index is within the accepted interval.
  3203.          *
  3204.          * @param index the index to check
  3205.          * @param lowerBound the lower bound of the interval
  3206.          * @param upperBound the upper bound of the interval
  3207.          * @return true if the index is between the lower and upper bounds, false otherwise
  3208.          */
  3209.         private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
  3210.             return index >= lowerBound && index <= upperBound;
  3211.         }

  3213.         /**Get the value of C<sup>j</sup>.
  3214.          *
  3215.          * @param j j index
  3216.          * @return C<sup>j</sup>
  3217.          */
  3218.         public T getCj(final int j) {
  3219.             return cCoef[0][j];
  3220.         }

  3222.         /**Get the value of dC<sup>j</sup> / da.
  3223.          *
  3224.          * @param j j index
  3225.          * @return dC<sup>j</sup> / da
  3226.          */
  3227.         public T getdCjdA(final int j) {
  3228.             return cCoef[1][j];
  3229.         }

  3231.         /**Get the value of dC<sup>j</sup> / dh.
  3232.          *
  3233.          * @param j j index
  3234.          * @return dC<sup>j</sup> / dh
  3235.          */
  3236.         public T getdCjdH(final int j) {
  3237.             return cCoef[2][j];
  3238.         }

  3240.         /**Get the value of dC<sup>j</sup> / dk.
  3241.          *
  3242.          * @param j j index
  3243.          * @return dC<sup>j</sup> / dk
  3244.          */
  3245.         public T getdCjdK(final int j) {
  3246.             return cCoef[3][j];
  3247.         }

  3249.         /**Get the value of dC<sup>j</sup> / dα.
  3250.          *
  3251.          * @param j j index
  3252.          * @return dC<sup>j</sup> / dα
  3253.          */
  3254.         public T getdCjdAlpha(final int j) {
  3255.             return cCoef[4][j];
  3256.         }

  3258.         /**Get the value of dC<sup>j</sup> / dβ.
  3259.          *
  3260.          * @param j j index
  3261.          * @return dC<sup>j</sup> / dβ
  3262.          */
  3263.         public T getdCjdBeta(final int j) {
  3264.             return cCoef[5][j];
  3265.         }

  3267.         /**Get the value of dC<sup>j</sup> / dγ.
  3268.          *
  3269.          * @param j j index
  3270.          * @return dC<sup>j</sup> / dγ
  3271.          */
  3272.         public T getdCjdGamma(final int j) {
  3273.             return cCoef[6][j];
  3274.         }

  3276.         /**Get the value of S<sup>j</sup>.
  3277.          *
  3278.          * @param j j index
  3279.          * @return S<sup>j</sup>
  3280.          */
  3281.         public T getSj(final int j) {
  3282.             return sCoef[0][j];
  3283.         }

  3285.         /**Get the value of dS<sup>j</sup> / da.
  3286.          *
  3287.          * @param j j index
  3288.          * @return dS<sup>j</sup> / da
  3289.          */
  3290.         public T getdSjdA(final int j) {
  3291.             return sCoef[1][j];
  3292.         }

  3294.         /**Get the value of dS<sup>j</sup> / dh.
  3295.          *
  3296.          * @param j j index
  3297.          * @return dS<sup>j</sup> / dh
  3298.          */
  3299.         public T getdSjdH(final int j) {
  3300.             return sCoef[2][j];
  3301.         }

  3303.         /**Get the value of dS<sup>j</sup> / dk.
  3304.          *
  3305.          * @param j j index
  3306.          * @return dS<sup>j</sup> / dk
  3307.          */
  3308.         public T getdSjdK(final int j) {
  3309.             return sCoef[3][j];
  3310.         }

  3312.         /**Get the value of dS<sup>j</sup> / dα.
  3313.          *
  3314.          * @param j j index
  3315.          * @return dS<sup>j</sup> / dα
  3316.          */
  3317.         public T getdSjdAlpha(final int j) {
  3318.             return sCoef[4][j];
  3319.         }

  3321.         /**Get the value of dS<sup>j</sup> / dβ.
  3322.          *
  3323.          * @param j j index
  3324.          * @return dS<sup>j</sup> / dβ
  3325.          */
  3326.         public T getdSjdBeta(final int j) {
  3327.             return sCoef[5][j];
  3328.         }

  3330.         /**Get the value of dS<sup>j</sup> /  dγ.
  3331.          *
  3332.          * @param j j index
  3333.          * @return dS<sup>j</sup> /  dγ
  3334.          */
  3335.         public T getdSjdGamma(final int j) {
  3336.             return sCoef[6][j];
  3337.         }
  3338.     }

  3340.     /** Coefficients valid for one time slot. */
  3341.     private static class Slot {

  3343.         /**The coefficients D<sub>i</sub>.
  3344.          * <p>
  3345.          * i corresponds to the equinoctial element, as follows:
  3346.          * - i=0 for a <br/>
  3347.          * - i=1 for k <br/>
  3348.          * - i=2 for h <br/>
  3349.          * - i=3 for q <br/>
  3350.          * - i=4 for p <br/>
  3351.          * - i=5 for λ <br/>
  3352.          * </p>
  3353.          */
  3354.         private final ShortPeriodicsInterpolatedCoefficient di;

  3356.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  3357.          * <p>
  3358.          * The constant term C<sub>i</sub>⁰ is also stored in this variable at index j = 0 <br>
  3359.          * The index order is cij[j][i] <br/>
  3360.          * i corresponds to the equinoctial element, as follows: <br/>
  3361.          * - i=0 for a <br/>
  3362.          * - i=1 for k <br/>
  3363.          * - i=2 for h <br/>
  3364.          * - i=3 for q <br/>
  3365.          * - i=4 for p <br/>
  3366.          * - i=5 for λ <br/>
  3367.          * </p>
  3368.          */
  3369.         private final ShortPeriodicsInterpolatedCoefficient[] cij;

  3371.         /** The coefficients S<sub>i</sub><sup>j</sup>.
  3372.          * <p>
  3373.          * The index order is sij[j][i] <br/>
  3374.          * i corresponds to the equinoctial element, as follows: <br/>
  3375.          * - i=0 for a <br/>
  3376.          * - i=1 for k <br/>
  3377.          * - i=2 for h <br/>
  3378.          * - i=3 for q <br/>
  3379.          * - i=4 for p <br/>
  3380.          * - i=5 for λ <br/>
  3381.          * </p>
  3382.          */
  3383.         private final ShortPeriodicsInterpolatedCoefficient[] sij;

  3385.         /** Simple constructor.
  3386.          *  @param maxFrequencyShortPeriodics maximum value for j index
  3387.          *  @param interpolationPoints number of points used in the interpolation process
  3388.          */
  3389.         Slot(final int maxFrequencyShortPeriodics, final int interpolationPoints) {

  3391.             final int rows = maxFrequencyShortPeriodics + 1;
  3392.             di  = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3393.             cij = new ShortPeriodicsInterpolatedCoefficient[rows];
  3394.             sij = new ShortPeriodicsInterpolatedCoefficient[rows];

  3396.             //Initialize the arrays
  3397.             for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
  3398.                 cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3399.                 sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3400.             }

  3402.         }

  3404.     }

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

  3409.         /**The coefficients D<sub>i</sub>.
  3410.          * <p>
  3411.          * i corresponds to the equinoctial element, as follows:
  3412.          * - i=0 for a <br/>
  3413.          * - i=1 for k <br/>
  3414.          * - i=2 for h <br/>
  3415.          * - i=3 for q <br/>
  3416.          * - i=4 for p <br/>
  3417.          * - i=5 for λ <br/>
  3418.          * </p>
  3419.          */
  3420.         private final FieldShortPeriodicsInterpolatedCoefficient<T> di;

  3422.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  3423.          * <p>
  3424.          * The constant term C<sub>i</sub>⁰ is also stored in this variable at index j = 0 <br>
  3425.          * The index order is cij[j][i] <br/>
  3426.          * i corresponds to the equinoctial element, as follows: <br/>
  3427.          * - i=0 for a <br/>
  3428.          * - i=1 for k <br/>
  3429.          * - i=2 for h <br/>
  3430.          * - i=3 for q <br/>
  3431.          * - i=4 for p <br/>
  3432.          * - i=5 for λ <br/>
  3433.          * </p>
  3434.          */
  3435.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] cij;

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

  3451.         /** Simple constructor.
  3452.          *  @param maxFrequencyShortPeriodics maximum value for j index
  3453.          *  @param interpolationPoints number of points used in the interpolation process
  3454.          */
  3455.         @SuppressWarnings("unchecked")
  3456.         FieldSlot(final int maxFrequencyShortPeriodics, final int interpolationPoints) {

  3458.             final int rows = maxFrequencyShortPeriodics + 1;
  3459.             di  = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3460.             cij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows);
  3461.             sij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows);

  3463.             //Initialize the arrays
  3464.             for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
  3465.                 cij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3466.                 sij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3467.             }

  3469.         }

  3471.     }

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

  3476.         /** The current value of the U function. <br/>
  3477.          * Needed for the short periodic contribution */
  3478.         private double U;

  3480.         /** dU / da. */
  3481.         private double dUda;

  3483.         /** dU / dk. */
  3484.         private double dUdk;

  3486.         /** dU / dh. */
  3487.         private double dUdh;

  3489.         /** dU / dAlpha. */
  3490.         private double dUdAl;

  3492.         /** dU / dBeta. */
  3493.         private double dUdBe;

  3495.         /** dU / dGamma. */
  3496.         private double dUdGa;

  3498.         /** Simple constuctor.
  3499.          *  @param date current date
  3500.          *  @param context container for attributes
  3501.          *  @param auxiliaryElements auxiliary elements related to the current orbit
  3502.          *  @param hansen initialization of hansen objects
  3503.          */
  3504.         UAnddU(final AbsoluteDate date,
  3505.                final DSSTZonalContext context,
  3506.                final AuxiliaryElements auxiliaryElements,
  3507.                final HansenObjects hansen) {

  3509.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);

  3511.             //Reset U
  3512.             U = 0.;

  3514.             // Gs and Hs coefficients
  3515.             final double[][] GsHs = CoefficientsFactory.computeGsHs(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta(), maxEccPowMeanElements);
  3516.             // Qns coefficients
  3517.             final double[][] Qns  = CoefficientsFactory.computeQns(auxiliaryElements.getGamma(), maxDegree, maxEccPowMeanElements);

  3519.             final double[] roaPow = new double[maxDegree + 1];
  3520.             roaPow[0] = 1.;
  3521.             for (int i = 1; i <= maxDegree; i++) {
  3522.                 roaPow[i] = context.getRoa() * roaPow[i - 1];
  3523.             }

  3525.             // Potential derivatives
  3526.             dUda  = 0.;
  3527.             dUdk  = 0.;
  3528.             dUdh  = 0.;
  3529.             dUdAl = 0.;
  3530.             dUdBe = 0.;
  3531.             dUdGa = 0.;

  3533.             for (int s = 0; s <= maxEccPowMeanElements; s++) {
  3534.                 //Initialize the Hansen roots
  3535.                 hansen.computeHansenObjectsInitValues(context, s);

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

  3540.                 // Compute Gs partial derivatives from 3.1-(9)
  3541.                 double dGsdh  = 0.;
  3542.                 double dGsdk  = 0.;
  3543.                 double dGsdAl = 0.;
  3544.                 double dGsdBe = 0.;
  3545.                 if (s > 0) {
  3546.                     // First get the G(s-1) and the H(s-1) coefficients
  3547.                     final double sxgsm1 = s * GsHs[0][s - 1];
  3548.                     final double sxhsm1 = s * GsHs[1][s - 1];
  3549.                     // Then compute derivatives
  3550.                     dGsdh  = auxiliaryElements.getBeta()  * sxgsm1 - auxiliaryElements.getAlpha() * sxhsm1;
  3551.                     dGsdk  = auxiliaryElements.getAlpha() * sxgsm1 + auxiliaryElements.getBeta()  * sxhsm1;
  3552.                     dGsdAl = auxiliaryElements.getK() * sxgsm1 - auxiliaryElements.getH() * sxhsm1;
  3553.                     dGsdBe = auxiliaryElements.getH() * sxgsm1 + auxiliaryElements.getK() * sxhsm1;
  3554.                 }

  3556.                 // Kronecker symbol (2 - delta(0,s))
  3557.                 final double d0s = (s == 0) ? 1 : 2;

  3559.                 for (int n = s + 2; n <= maxDegree; n++) {
  3560.                     // (n - s) must be even
  3561.                     if ((n - s) % 2 == 0) {

  3563.                         //Extract data from previous computation :
  3564.                         final double kns   = hansen.getHansenObjects()[s].getValue(-n - 1, context.getX());
  3565.                         final double dkns  = hansen.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  3567.                         final double vns   = Vns.get(new NSKey(n, s));
  3568.                         final double coef0 = d0s * roaPow[n] * vns * -harmonics.getUnnormalizedCnm(n, 0);
  3569.                         final double coef1 = coef0 * Qns[n][s];
  3570.                         final double coef2 = coef1 * kns;
  3571.                         final double coef3 = coef2 * gs;
  3572.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  3573.                         final double dqns  = Qns[n][s + 1];

  3575.                         // Compute U
  3576.                         U += coef3;
  3577.                         // Compute dU / da :
  3578.                         dUda += coef3 * (n + 1);
  3579.                         // Compute dU / dEx
  3580.                         dUdk += coef1 * (kns * dGsdk + auxiliaryElements.getK() * context.getXXX() * gs * dkns);
  3581.                         // Compute dU / dEy
  3582.                         dUdh += coef1 * (kns * dGsdh + auxiliaryElements.getH() * context.getXXX() * gs * dkns);
  3583.                         // Compute dU / dAlpha
  3584.                         dUdAl += coef2 * dGsdAl;
  3585.                         // Compute dU / dBeta
  3586.                         dUdBe += coef2 * dGsdBe;
  3587.                         // Compute dU / dGamma
  3588.                         dUdGa += coef0 * kns * dqns * gs;

  3590.                     }
  3591.                 }
  3592.             }

  3594.             // Multiply by -(μ / a)
  3595.             this.U = -context.getMuoa() * U;

  3597.             this.dUda = dUda *  context.getMuoa() / auxiliaryElements.getSma();
  3598.             this.dUdk = dUdk * -context.getMuoa();
  3599.             this.dUdh = dUdh * -context.getMuoa();
  3600.             this.dUdAl = dUdAl * -context.getMuoa();
  3601.             this.dUdBe = dUdBe * -context.getMuoa();
  3602.             this.dUdGa = dUdGa * -context.getMuoa();

  3604.         }

  3606.         /** Return value of U.
  3607.          * @return U
  3608.          */
  3609.         public double getU() {
  3610.             return U;
  3611.         }

  3613.         /** Return value of dU / da.
  3614.          * @return dUda
  3615.          */
  3616.         public double getdUda() {
  3617.             return dUda;
  3618.         }

  3620.         /** Return value of dU / dk.
  3621.          * @return dUdk
  3622.          */
  3623.         public double getdUdk() {
  3624.             return dUdk;
  3625.         }

  3627.         /** Return value of dU / dh.
  3628.          * @return dUdh
  3629.          */
  3630.         public double getdUdh() {
  3631.             return dUdh;
  3632.         }

  3634.         /** Return value of dU / dAlpha.
  3635.          * @return dUdAl
  3636.          */
  3637.         public double getdUdAl() {
  3638.             return dUdAl;
  3639.         }

  3641.         /** Return value of dU / dBeta.
  3642.          * @return dUdBe
  3643.          */
  3644.         public double getdUdBe() {
  3645.             return dUdBe;
  3646.         }

  3648.         /** Return value of dU / dGamma.
  3649.          * @return dUdGa
  3650.          */
  3651.         public double getdUdGa() {
  3652.             return dUdGa;
  3653.         }

  3655.     }

  3657.     /** Compute the derivatives of the gravitational potential U [Eq. 3.1-(6)].
  3658.      *  <p>
  3659.      *  The result is the array
  3660.      *  [dU/da, dU/dk, dU/dh, dU/dα, dU/dβ, dU/dγ]
  3661.      *  </p>
  3662.      */
  3663.     private class FieldUAnddU <T extends CalculusFieldElement<T>> {

  3665.          /** The current value of the U function. <br/>
  3666.           * Needed for the short periodic contribution */
  3667.         private T U;

  3669.          /** dU / da. */
  3670.         private T dUda;

  3672.          /** dU / dk. */
  3673.         private T dUdk;

  3675.          /** dU / dh. */
  3676.         private T dUdh;

  3678.          /** dU / dAlpha. */
  3679.         private T dUdAl;

  3681.          /** dU / dBeta. */
  3682.         private T dUdBe;

  3684.          /** dU / dGamma. */
  3685.         private T dUdGa;

  3687.          /** Simple constuctor.
  3688.           *  @param date current date
  3689.           *  @param context container for attributes
  3690.           *  @param auxiliaryElements auxiliary elements related to the current orbit
  3691.           *  @param hansen initialization of hansen objects
  3692.           */
  3693.         FieldUAnddU(final FieldAbsoluteDate<T> date,
  3694.                     final FieldDSSTZonalContext<T> context,
  3695.                     final FieldAuxiliaryElements<T> auxiliaryElements,
  3696.                     final FieldHansenObjects<T> hansen) {

  3698.             // Zero for initialization
  3699.             final Field<T> field = date.getField();
  3700.             final T zero = field.getZero();

  3702.             //spherical harmonics
  3703.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date.toAbsoluteDate());

  3705.             //Reset U
  3706.             U = zero;

  3708.             // Gs and Hs coefficients
  3709.             final T[][] GsHs = CoefficientsFactory.computeGsHs(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta(), maxEccPowMeanElements, field);
  3710.             // Qns coefficients
  3711.             final T[][] Qns  = CoefficientsFactory.computeQns(auxiliaryElements.getGamma(), maxDegree, maxEccPowMeanElements);

  3713.             final T[] roaPow = MathArrays.buildArray(field, maxDegree + 1);
  3714.             roaPow[0] = zero.add(1.);
  3715.             for (int i = 1; i <= maxDegree; i++) {
  3716.                 roaPow[i] = roaPow[i - 1].multiply(context.getRoa());
  3717.             }

  3719.             // Potential derivatives
  3720.             dUda  = zero;
  3721.             dUdk  = zero;
  3722.             dUdh  = zero;
  3723.             dUdAl = zero;
  3724.             dUdBe = zero;
  3725.             dUdGa = zero;

  3727.             for (int s = 0; s <= maxEccPowMeanElements; s++) {
  3728.                 //Initialize the Hansen roots
  3729.                 hansen.computeHansenObjectsInitValues(context, s);

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

  3734.                 // Compute Gs partial derivatives from 3.1-(9)
  3735.                 T dGsdh  = zero;
  3736.                 T dGsdk  = zero;
  3737.                 T dGsdAl = zero;
  3738.                 T dGsdBe = zero;
  3739.                 if (s > 0) {
  3740.                     // First get the G(s-1) and the H(s-1) coefficients
  3741.                     final T sxgsm1 = GsHs[0][s - 1].multiply(s);
  3742.                     final T sxhsm1 = GsHs[1][s - 1].multiply(s);
  3743.                     // Then compute derivatives
  3744.                     dGsdh  = sxgsm1.multiply(auxiliaryElements.getBeta()).subtract(sxhsm1.multiply(auxiliaryElements.getAlpha()));
  3745.                     dGsdk  = sxgsm1.multiply(auxiliaryElements.getAlpha()).add(sxhsm1.multiply(auxiliaryElements.getBeta()));
  3746.                     dGsdAl = sxgsm1.multiply(auxiliaryElements.getK()).subtract(sxhsm1.multiply(auxiliaryElements.getH()));
  3747.                     dGsdBe = sxgsm1.multiply(auxiliaryElements.getH()).add(sxhsm1.multiply(auxiliaryElements.getK()));
  3748.                 }

  3750.                 // Kronecker symbol (2 - delta(0,s))
  3751.                 final T d0s = zero.add((s == 0) ? 1 : 2);

  3753.                 for (int n = s + 2; n <= maxDegree; n++) {
  3754.                     // (n - s) must be even
  3755.                     if ((n - s) % 2 == 0) {

  3757.                         //Extract data from previous computation :
  3758.                         final T kns   = (T) hansen.getHansenObjects()[s].getValue(-n - 1, context.getX());
  3759.                         final T dkns  = (T) hansen.getHansenObjects()[s].getDerivative(-n - 1, context.getX());

  3761.                         final double vns   = Vns.get(new NSKey(n, s));
  3762.                         final T coef0 = d0s.multiply(roaPow[n]).multiply(vns).multiply(-harmonics.getUnnormalizedCnm(n, 0));
  3763.                         final T coef1 = coef0.multiply(Qns[n][s]);
  3764.                         final T coef2 = coef1.multiply(kns);
  3765.                         final T coef3 = coef2.multiply(gs);
  3766.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  3767.                         final T dqns  = Qns[n][s + 1];

  3769.                         // Compute U
  3770.                         U = U.add(coef3);
  3771.                         // Compute dU / da :
  3772.                         dUda  = dUda.add(coef3.multiply(n + 1));
  3773.                         // Compute dU / dEx
  3774.                         dUdk  = dUdk.add(coef1.multiply(dGsdk.multiply(kns).add(auxiliaryElements.getK().multiply(context.getXXX()).multiply(dkns).multiply(gs))));
  3775.                         // Compute dU / dEy
  3776.                         dUdh  = dUdh.add(coef1.multiply(dGsdh.multiply(kns).add(auxiliaryElements.getH().multiply(context.getXXX()).multiply(dkns).multiply(gs))));
  3777.                         // Compute dU / dAlpha
  3778.                         dUdAl = dUdAl.add(coef2.multiply(dGsdAl));
  3779.                         // Compute dU / dBeta
  3780.                         dUdBe = dUdBe.add(coef2.multiply(dGsdBe));
  3781.                         // Compute dU / dGamma
  3782.                         dUdGa = dUdGa.add(coef0.multiply(kns).multiply(dqns).multiply(gs));
  3783.                     }
  3784.                 }
  3785.             }

  3787.             // Multiply by -(μ / a)
  3788.             U = U.multiply(context.getMuoa().negate());

  3790.             dUda  = dUda.multiply(context.getMuoa().divide(auxiliaryElements.getSma()));
  3791.             dUdk  = dUdk.multiply(context.getMuoa()).negate();
  3792.             dUdh  = dUdh.multiply(context.getMuoa()).negate();
  3793.             dUdAl = dUdAl.multiply(context.getMuoa()).negate();
  3794.             dUdBe = dUdBe.multiply(context.getMuoa()).negate();
  3795.             dUdGa = dUdGa.multiply(context.getMuoa()).negate();

  3797.         }

  3799.         /** Return value of U.
  3800.          * @return U
  3801.          */
  3802.         public T getU() {
  3803.             return U;
  3804.         }

  3806.          /** Return value of dU / da.
  3807.           * @return dUda
  3808.           */
  3809.         public T getdUda() {
  3810.             return dUda;
  3811.         }

  3813.          /** Return value of dU / dk.
  3814.           * @return dUdk
  3815.           */
  3816.         public T getdUdk() {
  3817.             return dUdk;
  3818.         }

  3820.          /** Return value of dU / dh.
  3821.           * @return dUdh
  3822.           */
  3823.         public T getdUdh() {
  3824.             return dUdh;
  3825.         }

  3827.          /** Return value of dU / dAlpha.
  3828.           * @return dUdAl
  3829.           */
  3830.         public T getdUdAl() {
  3831.             return dUdAl;
  3832.         }

  3834.          /** Return value of dU / dBeta.
  3835.           * @return dUdBe
  3836.           */
  3837.         public T getdUdBe() {
  3838.             return dUdBe;
  3839.         }

  3841.          /** Return value of dU / dGamma.
  3842.           * @return dUdGa
  3843.           */
  3844.         public T getdUdGa() {
  3845.             return dUdGa;
  3846.         }

  3848.     }

  3850.     /** Computes init values of the Hansen Objects. */
  3851.     private class HansenObjects {

  3853.         /** An array that contains the objects needed to build the Hansen coefficients. <br/>
  3854.          * The index is s*/
  3855.         private HansenZonalLinear[] hansenObjects;

  3857.         /** Simple constructor. */
  3858.         HansenObjects() {
  3859.             this.hansenObjects = new HansenZonalLinear[maxEccPow + 1];
  3860.             for (int s = 0; s <= maxEccPow; s++) {
  3861.                 this.hansenObjects[s] = new HansenZonalLinear(maxDegree, s);
  3862.             }
  3863.         }

  3865.         /** Compute init values for hansen objects.
  3866.          * @param context container for attributes
  3867.          * @param element element of the array to compute the init values
  3868.          */
  3869.         public void computeHansenObjectsInitValues(final DSSTZonalContext context, final int element) {
  3870.             hansenObjects[element].computeInitValues(context.getX());
  3871.         }

  3873.         /** Get the Hansen Objects.
  3874.          * @return hansenObjects
  3875.          */
  3876.         public HansenZonalLinear[] getHansenObjects() {
  3877.             return hansenObjects;
  3878.         }

  3880.     }

  3882.     /** Computes init values of the Hansen Objects.
  3883.      * @param <T> type of the elements
  3884.      */
  3885.     private class FieldHansenObjects<T extends CalculusFieldElement<T>> {

  3887.         /** An array that contains the objects needed to build the Hansen coefficients. <br/>
  3888.          * The index is s*/
  3889.         private FieldHansenZonalLinear<T>[] hansenObjects;

  3891.         /** Simple constructor.
  3892.          * @param field field used by default
  3893.          */
  3894.         @SuppressWarnings("unchecked")
  3895.         FieldHansenObjects(final Field<T> field) {
  3896.             this.hansenObjects = (FieldHansenZonalLinear<T>[]) Array.newInstance(FieldHansenZonalLinear.class, maxEccPow + 1);
  3897.             for (int s = 0; s <= maxEccPow; s++) {
  3898.                 this.hansenObjects[s] = new FieldHansenZonalLinear<>(maxDegree, s, field);
  3899.             }
  3900.         }

  3902.         /** Compute init values for hansen objects.
  3903.          * @param context container for attributes
  3904.          * @param element element of the array to compute the init values
  3905.          */
  3906.         public void computeHansenObjectsInitValues(final FieldDSSTZonalContext<T> context, final int element) {
  3907.             hansenObjects[element].computeInitValues(context.getX());
  3908.         }

  3910.         /** Get the Hansen Objects.
  3911.          * @return hansenObjects
  3912.          */
  3913.         public FieldHansenZonalLinear<T>[] getHansenObjects() {
  3914.             return hansenObjects;
  3915.         }

  3917.     }

  3919. }
Created with JaCoCo 0.8.11.202310140853