DSSTZonal.java

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

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

  29. import org.hipparchus.Field;
  30. import org.hipparchus.RealFieldElement;
  31. import org.hipparchus.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.events.EventDetector;
  48. import org.orekit.propagation.events.FieldEventDetector;
  49. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  50. import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
  51. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
  52. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory.NSKey;
  53. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  54. import org.orekit.propagation.semianalytical.dsst.utilities.FieldCjSjCoefficient;
  55. import org.orekit.propagation.semianalytical.dsst.utilities.FieldGHIJjsPolynomials;
  56. import org.orekit.propagation.semianalytical.dsst.utilities.FieldLnsCoefficients;
  57. import org.orekit.propagation.semianalytical.dsst.utilities.FieldShortPeriodicsInterpolatedCoefficient;
  58. import org.orekit.propagation.semianalytical.dsst.utilities.GHIJjsPolynomials;
  59. import org.orekit.propagation.semianalytical.dsst.utilities.LnsCoefficients;
  60. import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
  61. import org.orekit.propagation.semianalytical.dsst.utilities.UpperBounds;
  62. import org.orekit.propagation.semianalytical.dsst.utilities.hansen.FieldHansenZonalLinear;
  63. import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenZonalLinear;
  64. import org.orekit.time.AbsoluteDate;
  65. import org.orekit.time.FieldAbsoluteDate;
  66. import org.orekit.utils.FieldTimeSpanMap;
  67. import org.orekit.utils.ParameterDriver;
  68. import org.orekit.utils.TimeSpanMap;

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

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

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

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

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

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

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

  113.     /** Provider for spherical harmonics. */
  114.     private final UnnormalizedSphericalHarmonicsProvider provider;

  116.     /** Maximal degree to consider for harmonics potential. */
  117.     private final int maxDegree;

  119.     /** Maximal degree to consider for harmonics potential. */
  120.     private final int maxDegreeShortPeriodics;

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

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

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

  131.     /** Highest power of the eccentricity. */
  132.     private int maxEccPow;

  134.     /** Short period terms. */
  135.     private ZonalShortPeriodicCoefficients zonalSPCoefs;

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

  140.     /** Driver for gravitational parameter. */
  141.     private final ParameterDriver gmParameterDriver;

  143.     /** Hansen objects. */
  144.     private HansenObjects hansen;

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

  149.     /** Flag for force model initialization with field elements. */
  150.     private boolean pendingInitialization;

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

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

  186.         gmParameterDriver = new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
  187.                                                 provider.getMu(), MU_SCALE,
  188.                                                 0.0, Double.POSITIVE_INFINITY);

  190.         // Vns coefficients
  191.         this.Vns = CoefficientsFactory.computeVns(provider.getMaxDegree() + 1);

  193.         this.provider  = provider;
  194.         this.maxDegree = provider.getMaxDegree();

  196.         checkIndexRange(maxDegreeShortPeriodics, 2, provider.getMaxDegree());
  197.         this.maxDegreeShortPeriodics = maxDegreeShortPeriodics;

  199.         checkIndexRange(maxEccPowShortPeriodics, 0, maxDegreeShortPeriodics - 1);
  200.         this.maxEccPowShortPeriodics = maxEccPowShortPeriodics;

  202.         checkIndexRange(maxFrequencyShortPeriodics, 1, 2 * maxDegreeShortPeriodics + 1);
  203.         this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;

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


  209.         pendingInitialization = true;

  211.         zonalFieldSPCoefs = new HashMap<>();
  212.         fieldHansen       = new HashMap<>();

  214.     }

  216.     /** Check an index range.
  217.      * @param index index value
  218.      * @param min minimum value for index
  219.      * @param max maximum value for index
  220.      */
  221.     private void checkIndexRange(final int index, final int min, final int max) {
  222.         if (index < min || index > max) {
  223.             throw new OrekitException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, index, min, max);
  224.         }
  225.     }

  227.     /** Get the spherical harmonics provider.
  228.      *  @return the spherical harmonics provider
  229.      */
  230.     public UnnormalizedSphericalHarmonicsProvider getProvider() {
  231.         return provider;
  232.     }

  234.     /** {@inheritDoc}
  235.      *  <p>
  236.      *  Computes the highest power of the eccentricity to appear in the truncated
  237.      *  analytical power series expansion.
  238.      *  </p>
  239.      *  <p>
  240.      *  This method computes the upper value for the central body potential and
  241.      *  determines the maximal power for the eccentricity producing potential
  242.      *  terms bigger than a defined tolerance.
  243.      *  </p>
  244.      */
  245.     @Override
  246.     public List<ShortPeriodTerms> initialize(final AuxiliaryElements auxiliaryElements,
  247.                                              final PropagationType type,
  248.                                              final double[] parameters) {

  250.         computeMeanElementsTruncations(auxiliaryElements, parameters);

  252.         switch (type) {
  253.             case MEAN:
  254.                 maxEccPow = maxEccPowMeanElements;
  255.                 break;
  256.             case OSCULATING:
  257.                 maxEccPow = FastMath.max(maxEccPowMeanElements, maxEccPowShortPeriodics);
  258.                 break;
  259.             default:
  260.                 throw new OrekitInternalError(null);
  261.         }

  263.         hansen = new HansenObjects();

  265.         final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
  266.         zonalSPCoefs = new ZonalShortPeriodicCoefficients(maxFrequencyShortPeriodics,
  267.                                                           INTERPOLATION_POINTS,
  268.                                                           new TimeSpanMap<Slot>(new Slot(maxFrequencyShortPeriodics,
  269.                                                                                          INTERPOLATION_POINTS)));
  270.         list.add(zonalSPCoefs);
  271.         return list;

  273.     }

  275.     /** {@inheritDoc}
  276.      *  <p>
  277.      *  Computes the highest power of the eccentricity to appear in the truncated
  278.      *  analytical power series expansion.
  279.      *  </p>
  280.      *  <p>
  281.      *  This method computes the upper value for the central body potential and
  282.      *  determines the maximal power for the eccentricity producing potential
  283.      *  terms bigger than a defined tolerance.
  284.      *  </p>
  285.      */
  286.     @Override
  287.     public <T extends RealFieldElement<T>> List<FieldShortPeriodTerms<T>> initialize(final FieldAuxiliaryElements<T> auxiliaryElements,
  288.                                                                                      final PropagationType type,
  289.                                                                                      final T[] parameters) {

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

  294.         if (pendingInitialization == true) {
  295.             computeMeanElementsTruncations(auxiliaryElements, parameters, field);

  297.             switch (type) {
  298.                 case MEAN:
  299.                     maxEccPow = maxEccPowMeanElements;
  300.                     break;
  301.                 case OSCULATING:
  302.                     maxEccPow = FastMath.max(maxEccPowMeanElements, maxEccPowShortPeriodics);
  303.                     break;
  304.                 default:
  305.                     throw new OrekitInternalError(null);
  306.             }

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

  310.             pendingInitialization = false;
  311.         }

  313.         final FieldZonalShortPeriodicCoefficients<T> fzspc =
  314.                         new FieldZonalShortPeriodicCoefficients<>(maxFrequencyShortPeriodics,
  315.                                                                   INTERPOLATION_POINTS,
  316.                                                                   new FieldTimeSpanMap<>(new FieldSlot<>(maxFrequencyShortPeriodics,
  317.                                                                                                          INTERPOLATION_POINTS),
  318.                                                                                          field));
  319.         zonalFieldSPCoefs.put(field, fzspc);
  320.         return Collections.singletonList(fzspc);

  322.     }

  324.     /** Compute indices truncations for mean elements computations.
  325.      * @param auxiliaryElements auxiliary elements
  326.      * @param parameters values of the force model parameters
  327.      */
  328.     private void computeMeanElementsTruncations(final AuxiliaryElements auxiliaryElements, final double[] parameters) {

  330.         final DSSTZonalContext context = new DSSTZonalContext(auxiliaryElements, provider, parameters);
  331.         //Compute the max eccentricity power for the mean element rate expansion
  332.         if (maxDegree == 2) {
  333.             maxEccPowMeanElements = 0;
  334.         } else {
  335.             // Initializes specific parameters.
  336.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(auxiliaryElements.getDate());

  338.             // Utilities for truncation
  339.             final double ax2or = 2. * auxiliaryElements.getSma() / provider.getAe();
  340.             double xmuran = parameters[0] / auxiliaryElements.getSma();
  341.             // Set a lower bound for eccentricity
  342.             final double eo2  = FastMath.max(0.0025, auxiliaryElements.getEcc() / 2.);
  343.             final double x2o2 = context.getXX() / 2.;
  344.             final double[] eccPwr = new double[maxDegree + 1];
  345.             final double[] chiPwr = new double[maxDegree + 1];
  346.             final double[] hafPwr = new double[maxDegree + 1];
  347.             eccPwr[0] = 1.;
  348.             chiPwr[0] = context.getX();
  349.             hafPwr[0] = 1.;
  350.             for (int i = 1; i <= maxDegree; i++) {
  351.                 eccPwr[i] = eccPwr[i - 1] * eo2;
  352.                 chiPwr[i] = chiPwr[i - 1] * x2o2;
  353.                 hafPwr[i] = hafPwr[i - 1] * 0.5;
  354.                 xmuran  /= ax2or;
  355.             }

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

  421.             maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
  422.         }
  423.     }

  425.     /** Compute indices truncations for mean elements computations.
  426.      * @param <T> type of the elements
  427.      * @param auxiliaryElements auxiliary elements
  428.      * @param parameters values of the force model parameters
  429.      * @param field field used by default
  430.      */
  431.     private <T extends RealFieldElement<T>> void computeMeanElementsTruncations(final FieldAuxiliaryElements<T> auxiliaryElements,
  432.                                                                                 final T[] parameters,
  433.                                                                                 final Field<T> field) {

  435.         final T zero = field.getZero();
  436.         final FieldDSSTZonalContext<T> context = new FieldDSSTZonalContext<>(auxiliaryElements, provider, parameters);
  437.         //Compute the max eccentricity power for the mean element rate expansion
  438.         if (maxDegree == 2) {
  439.             maxEccPowMeanElements = 0;
  440.         } else {
  441.             // Initializes specific parameters.
  442.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(auxiliaryElements.getDate().toAbsoluteDate());

  444.             // Utilities for truncation
  445.             final T ax2or = auxiliaryElements.getSma().multiply(2.).divide(provider.getAe());
  446.             T xmuran = parameters[0].divide(auxiliaryElements.getSma());
  447.             // Set a lower bound for eccentricity
  448.             final T eo2  = FastMath.max(zero.add(0.0025), auxiliaryElements.getEcc().divide(2.));
  449.             final T x2o2 = context.getXX().divide(2.);
  450.             final T[] eccPwr = MathArrays.buildArray(field, maxDegree + 1);
  451.             final T[] chiPwr = MathArrays.buildArray(field, maxDegree + 1);
  452.             final T[] hafPwr = MathArrays.buildArray(field, maxDegree + 1);
  453.             eccPwr[0] = zero.add(1.);
  454.             chiPwr[0] = context.getX();
  455.             hafPwr[0] = zero.add(1.);
  456.             for (int i = 1; i <= maxDegree; i++) {
  457.                 eccPwr[i] = eccPwr[i - 1].multiply(eo2);
  458.                 chiPwr[i] = chiPwr[i - 1].multiply(x2o2);
  459.                 hafPwr[i] = hafPwr[i - 1].multiply(0.5);
  460.                 xmuran  = xmuran.divide(ax2or);
  461.             }

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

  527.             maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
  528.         }
  529.     }

  531.     /** Performs initialization at each integration step for the current force model.
  532.      *  <p>
  533.      *  This method aims at being called before mean elements rates computation.
  534.      *  </p>
  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 DSSTZonalContext initializeStep(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
  540.         return new DSSTZonalContext(auxiliaryElements, provider, parameters);
  541.     }

  543.     /** Performs initialization at each integration step for the current force model.
  544.      *  <p>
  545.      *  This method aims at being called before mean elements rates computation.
  546.      *  </p>
  547.      *  @param <T> type of the elements
  548.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  549.      *  @param parameters values of the force model parameters
  550.      *  @return new force model context
  551.      */
  552.     private <T extends RealFieldElement<T>> FieldDSSTZonalContext<T> initializeStep(final FieldAuxiliaryElements<T> auxiliaryElements,
  553.                                                                                     final T[] parameters) {
  554.         return new FieldDSSTZonalContext<>(auxiliaryElements, provider, parameters);
  555.     }

  557.     /** {@inheritDoc} */
  558.     @Override
  559.     public double[] getMeanElementRate(final SpacecraftState spacecraftState,
  560.                                        final AuxiliaryElements auxiliaryElements, final double[] parameters) {

  562.         // Container of attributes
  563.         final DSSTZonalContext context = initializeStep(auxiliaryElements, parameters);
  564.         // Access to potential U derivatives
  565.         final UAnddU udu = new UAnddU(spacecraftState.getDate(), context, auxiliaryElements, hansen);

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

  569.     }

  571.     /** {@inheritDoc} */
  572.     @Override
  573.     public <T extends RealFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> spacecraftState,
  574.                                                                   final FieldAuxiliaryElements<T> auxiliaryElements,
  575.                                                                   final T[] parameters) {

  577.         // Field used by default
  578.         final Field<T> field = auxiliaryElements.getDate().getField();
  579.         // Container of attributes
  580.         final FieldDSSTZonalContext<T> context = initializeStep(auxiliaryElements, parameters);

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

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

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

  590.     }

  592.     /** {@inheritDoc} */
  593.     @Override
  594.     public EventDetector[] getEventsDetectors() {
  595.         return null;
  596.     }

  598.     /** {@inheritDoc} */
  599.     @Override
  600.     public <T extends RealFieldElement<T>> FieldEventDetector<T>[] getFieldEventsDetectors(final Field<T> field) {
  601.         return null;
  602.     }

  604.     /** Compute the mean element rates.
  605.      * @param date current date
  606.      * @param context container for attributes
  607.      * @param udu derivatives of the gravitational potential U
  608.      * @return the mean element rates
  609.      */
  610.     private double[] computeMeanElementRates(final AbsoluteDate date,
  611.                                              final DSSTZonalContext context,
  612.                                              final UAnddU udu) {

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

  617.         // Compute cross derivatives [Eq. 2.2-(8)]
  618.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  619.         final double UAlphaGamma   = auxiliaryElements.getAlpha() * udu.getdUdGa() - auxiliaryElements.getGamma() * udu.getdUdAl();
  620.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  621.         final double UBetaGamma    =  auxiliaryElements.getBeta() * udu.getdUdGa() - auxiliaryElements.getGamma() * udu.getdUdBe();
  622.         // Common factor
  623.         final double pUAGmIqUBGoAB = (auxiliaryElements.getP() * UAlphaGamma - I * auxiliaryElements.getQ() * UBetaGamma) * context.getOoAB();

  625.         // Compute mean elements rates [Eq. 3.1-(1)]
  626.         final double da =  0.;
  627.         final double dh =  context.getBoA() * udu.getdUdk() + auxiliaryElements.getK() * pUAGmIqUBGoAB;
  628.         final double dk = -context.getBoA() * udu.getdUdh() - auxiliaryElements.getH() * pUAGmIqUBGoAB;
  629.         final double dp =  context.getMCo2AB() * UBetaGamma;
  630.         final double dq =  context.getMCo2AB() * UAlphaGamma * I;
  631.         final double dM =  context.getM2aoA() * udu.getdUda() + context.getBoABpo() * (auxiliaryElements.getH() * udu.getdUdh() + auxiliaryElements.getK() * udu.getdUdk()) + pUAGmIqUBGoAB;

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

  636.     /** Compute the mean element rates.
  637.      * @param <T> type of the elements
  638.      * @param date current date
  639.      * @param context container for attributes
  640.      * @param udu derivatives of the gravitational potential U
  641.      * @return the mean element rates
  642.      */
  643.     private <T extends RealFieldElement<T>> T[] computeMeanElementRates(final FieldAbsoluteDate<T> date,
  644.                                                                         final FieldDSSTZonalContext<T> context,
  645.                                                                         final FieldUAnddU<T> udu) {

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

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

  653.         // Compute cross derivatives [Eq. 2.2-(8)]
  654.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  655.         final T UAlphaGamma   = udu.getdUdGa().multiply(auxiliaryElements.getAlpha()).subtract(udu.getdUdAl().multiply(auxiliaryElements.getGamma()));
  656.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  657.         final T UBetaGamma    =  udu.getdUdGa().multiply(auxiliaryElements.getBeta()).subtract(udu.getdUdBe().multiply(auxiliaryElements.getGamma()));
  658.         // Common factor
  659.         final T pUAGmIqUBGoAB = (UAlphaGamma.multiply(auxiliaryElements.getP()).subtract(UBetaGamma.multiply(I).multiply(auxiliaryElements.getQ()))).multiply(context.getOoAB());

  661.         // Compute mean elements rates [Eq. 3.1-(1)]
  662.         final T da =  field.getZero();
  663.         final T dh =  udu.getdUdk().multiply(context.getBoA()).add(pUAGmIqUBGoAB.multiply(auxiliaryElements.getK()));
  664.         final T dk =  (udu.getdUdh().multiply(context.getBoA()).negate()).subtract(pUAGmIqUBGoAB.multiply(auxiliaryElements.getH()));
  665.         final T dp =  UBetaGamma.multiply(context.getMCo2AB());
  666.         final T dq =  UAlphaGamma.multiply(context.getMCo2AB()).multiply(I);
  667.         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()));

  669.         final T[] elements =  MathArrays.buildArray(field, 6);
  670.         elements[0] = da;
  671.         elements[1] = dk;
  672.         elements[2] = dh;
  673.         elements[3] = dq;
  674.         elements[4] = dp;
  675.         elements[5] = dM;

  677.         return elements;
  678.     }

  680.     /** {@inheritDoc} */
  681.     @Override
  682.     public void registerAttitudeProvider(final AttitudeProvider attitudeProvider) {
  683.         //nothing is done since this contribution is not sensitive to attitude
  684.     }

  686.     /** Check if an index is within the accepted interval.
  687.      *
  688.      * @param index the index to check
  689.      * @param lowerBound the lower bound of the interval
  690.      * @param upperBound the upper bound of the interval
  691.      * @return true if the index is between the lower and upper bounds, false otherwise
  692.      */
  693.     private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
  694.         return index >= lowerBound && index <= upperBound;
  695.     }

  697.     /** {@inheritDoc} */
  698.     @Override
  699.     public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {

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

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

  707.             // Container of attributes
  708.             final DSSTZonalContext context = initializeStep(auxiliaryElements, parameters);

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

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

  716.             // Compute Di
  717.             computeDiCoefficients(meanState.getDate(), slot, context, udu);

  719.             // generate the Cij and Sij coefficients
  720.             final FourierCjSjCoefficients cjsj = new FourierCjSjCoefficients(meanState.getDate(),
  721.                                                                              maxDegreeShortPeriodics,
  722.                                                                              maxEccPowShortPeriodics,
  723.                                                                              maxFrequencyShortPeriodics,
  724.                                                                              context,
  725.                                                                              hansen);

  727.             computeCijSijCoefficients(meanState.getDate(), slot, cjsj, rhoSigma, context, auxiliaryElements, udu);
  728.         }

  730.     }

  732.     /** {@inheritDoc} */
  733.     @Override
  734.     @SuppressWarnings("unchecked")
  735.     public <T extends RealFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
  736.                                                                        final FieldSpacecraftState<T>... meanStates) {

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

  741.         final FieldZonalShortPeriodicCoefficients<T> fzspc = (FieldZonalShortPeriodicCoefficients<T>) zonalFieldSPCoefs.get(field);
  742.         final FieldSlot<T> slot = fzspc.createSlot(meanStates);
  743.         for (final FieldSpacecraftState<T> meanState : meanStates) {

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

  748.             // Container of attributes
  749.             final FieldDSSTZonalContext<T> context = initializeStep(auxiliaryElements, parameters);

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

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

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

  759.             // Compute Di
  760.             computeDiCoefficients(meanState.getDate(), slot, context, field, udu);

  762.             // generate the Cij and Sij coefficients
  763.             final FieldFourierCjSjCoefficients<T> cjsj = new FieldFourierCjSjCoefficients<>(meanState.getDate(),
  764.                                                                                             maxDegreeShortPeriodics,
  765.                                                                                             maxEccPowShortPeriodics,
  766.                                                                                             maxFrequencyShortPeriodics,
  767.                                                                                             context,
  768.                                                                                             fho);


  771.             computeCijSijCoefficients(meanState.getDate(), slot, cjsj, rhoSigma, context, auxiliaryElements, field, udu);
  772.         }
  773.     }

  775.     /** {@inheritDoc} */
  776.     public ParameterDriver[] getParametersDrivers() {
  777.         return new ParameterDriver[] {
  778.             gmParameterDriver
  779.         };
  780.     }

  782.     /** Generate the values for the D<sub>i</sub> coefficients.
  783.      * @param date target date
  784.      * @param slot slot to which the coefficients belong
  785.      * @param context container for attributes
  786.      * @param udu derivatives of the gravitational potential U
  787.      */
  788.     private void computeDiCoefficients(final AbsoluteDate date,
  789.                                        final Slot slot,
  790.                                        final DSSTZonalContext context,
  791.                                        final UAnddU udu) {

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

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

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

  801.             if (i == 5) {
  802.                 currentDi[i] += -1.5 * 2 * udu.getU() * context.getOON2A2();
  803.             }

  805.         }

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

  809.     }

  811.     /** Generate the values for the D<sub>i</sub> coefficients.
  812.      * @param <T> type of the elements
  813.      * @param date target date
  814.      * @param slot slot to which the coefficients belong
  815.      * @param context container for attributes
  816.      * @param field field used by default
  817.      * @param udu derivatives of the gravitational potential U
  818.      */
  819.     private <T extends RealFieldElement<T>> void computeDiCoefficients(final FieldAbsoluteDate<T> date,
  820.                                                                        final FieldSlot<T> slot,
  821.                                                                        final FieldDSSTZonalContext<T> context,
  822.                                                                        final Field<T> field,
  823.                                                                        final FieldUAnddU<T> udu) {

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

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

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

  833.             if (i == 5) {
  834.                 currentDi[i] = currentDi[i].add(context.getOON2A2().multiply(udu.getU()).multiply(2.).multiply(-1.5));
  835.             }

  837.         }

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

  841.     }

  843.     /**
  844.      * Generate the values for the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
  845.      * @param date date of computation
  846.      * @param slot slot to which the coefficients belong
  847.      * @param cjsj Fourier coefficients
  848.      * @param rhoSigma ρ<sub>j</sub> and σ<sub>j</sub>
  849.      * @param context container for attributes
  850.      * @param auxiliaryElements auxiliary elements related to the current orbit
  851.      * @param udu derivatives of the gravitational potential U
  852.      */
  853.     private void computeCijSijCoefficients(final AbsoluteDate date, final Slot slot,
  854.                                            final FourierCjSjCoefficients cjsj,
  855.                                            final double[][] rhoSigma, final DSSTZonalContext context,
  856.                                            final AuxiliaryElements auxiliaryElements,
  857.                                            final UAnddU udu) {

  859.         final int nMax = maxDegreeShortPeriodics;

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

  865.             // Create local arrays
  866.             final double[] currentCij = new double[] {0., 0., 0., 0., 0., 0.};
  867.             final double[] currentSij = new double[] {0., 0., 0., 0., 0., 0.};

  869.             // j == 1
  870.             if (j == 1) {
  871.                 final double coef1 = 4 * auxiliaryElements.getK() * udu.getU() - context.getHK() * cjsj.getCj(1) + context.getK2MH2O2() * cjsj.getSj(1);
  872.                 final double coef2 = 4 * auxiliaryElements.getH() * udu.getU() + context.getK2MH2O2() * cjsj.getCj(1) + context.getHK() * cjsj.getSj(1);
  873.                 final double coef3 = (auxiliaryElements.getK() * cjsj.getCj(1) + auxiliaryElements.getH() * cjsj.getSj(1)) / 4.;
  874.                 final double coef4 = (8 * udu.getU() - auxiliaryElements.getH() * cjsj.getCj(1) + auxiliaryElements.getK() * cjsj.getSj(1)) / 4.;

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

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

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

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

  893.             // j == 2
  894.             if (j == 2) {
  895.                 final double coef1 = context.getK2MH2() * udu.getU();
  896.                 final double coef2 = 2 * context.getHK() * udu.getU();
  897.                 final double coef3 = auxiliaryElements.getH() * udu.getU() / 2;
  898.                 final double coef4 = auxiliaryElements.getK() * udu.getU() / 2;

  900.                 //Coefficients for a
  901.                 currentCij[0] += coef1;
  902.                 currentSij[0] += coef2;

  904.                 //Coefficients for k
  905.                 currentCij[1] += coef4;
  906.                 currentSij[1] += coef3;

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

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

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

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

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

  932.                 //Coefficients for h
  933.                 currentCij[2] -= coef3;
  934.                 currentSij[2] += coef4;

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

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

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

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

  956.                 //Coefficients for h
  957.                 currentCij[2] -= coef3;
  958.                 currentSij[2] -= coef4;

  960.                 //Coefficients for λ
  961.                 currentCij[5] -= coef2 / 2;
  962.                 currentSij[5] += -coef1 / 2;
  963.             }

  965.             // j between 2 and 2N
  966.             if (isBetween(j, 2, 2 * nMax) && j - 1 < cjsj.jMax) {
  967.                 final double coef1 = 2 * ( j - 1 ) * (auxiliaryElements.getH() * cjsj.getCj(j - 1) + auxiliaryElements.getK() * cjsj.getSj(j - 1));
  968.                 final double coef2 = 2 * ( j - 1 ) * (auxiliaryElements.getK() * cjsj.getCj(j - 1) - auxiliaryElements.getH() * cjsj.getSj(j - 1));
  969.                 final double coef3 = ( j - 1 ) * cjsj.getCj(j - 1);
  970.                 final double coef4 = ( j - 1 ) * cjsj.getSj(j - 1);

  972.                 //Coefficients for a
  973.                 currentCij[0] += coef1;
  974.                 currentSij[0] -= coef2;

  976.                 //Coefficients for k
  977.                 currentCij[1] += coef4;
  978.                 currentSij[1] -= coef3;

  980.                 //Coefficients for h
  981.                 currentCij[2] += coef3;
  982.                 currentSij[2] += coef4;

  984.                 //Coefficients for λ
  985.                 currentCij[5] += coef2 / 2;
  986.                 currentSij[5] += coef1 / 2;
  987.             }

  989.             // j between 3 and 2N + 1
  990.             if (isBetween(j, 3, 2 * nMax + 1) && j - 2 < cjsj.jMax) {
  991.                 final double coef1 = ( j - 2 ) * (context.getHK() * cjsj.getCj(j - 2) + context.getK2MH2O2() * cjsj.getSj(j - 2));
  992.                 final double coef2 = ( j - 2 ) * (-context.getK2MH2O2() * cjsj.getCj(j - 2) + context.getHK() * cjsj.getSj(j - 2));
  993.                 final double coef3 = ( j - 2 ) * (auxiliaryElements.getK() * cjsj.getCj(j - 2) - auxiliaryElements.getH() * cjsj.getSj(j - 2)) / 4;
  994.                 final double coef4 = ( j - 2 ) * (auxiliaryElements.getH() * cjsj.getCj(j - 2) + auxiliaryElements.getK() * cjsj.getSj(j - 2)) / 4;
  995.                 final double coef5 = ( j - 2 ) * (context.getK2MH2O2() * cjsj.getCj(j - 2) - context.getHK() * cjsj.getSj(j - 2));

  997.                 //Coefficients for a
  998.                 currentCij[0] += coef1;
  999.                 currentSij[0] += coef2;

  1001.                 //Coefficients for k
  1002.                 currentCij[1] += coef4;
  1003.                 currentSij[1] += -coef3;

  1005.                 //Coefficients for h
  1006.                 currentCij[2] += coef3;
  1007.                 currentSij[2] += coef4;

  1009.                 //Coefficients for λ
  1010.                 currentCij[5] += coef5 / 2;
  1011.                 currentSij[5] += coef1 / 2;
  1012.             }

  1014.             //multiply by the common factor
  1015.             //for a (i == 0) -> χ³ / (n² * a)
  1016.             currentCij[0] *= context.getX3ON2A();
  1017.             currentSij[0] *= context.getX3ON2A();
  1018.             //for k (i == 1) -> χ / (n² * a²)
  1019.             currentCij[1] *= context.getXON2A2();
  1020.             currentSij[1] *= context.getXON2A2();
  1021.             //for h (i == 2) -> χ / (n² * a²)
  1022.             currentCij[2] *= context.getXON2A2();
  1023.             currentSij[2] *= context.getXON2A2();
  1024.             //for λ (i == 5) -> (χ²) / (n² * a² * (χ + 1 ) )
  1025.             currentCij[5] *= context.getX2ON2A2XP1();
  1026.             currentSij[5] *= context.getX2ON2A2XP1();

  1028.             // j is between 1 and 2 * N - 1
  1029.             if (isBetween(j, 1, 2 * nMax - 1) && j < cjsj.jMax) {
  1030.                 // Compute cross derivatives
  1031.                 // Cj(alpha,gamma) = alpha * dC/dgamma - gamma * dC/dalpha
  1032.                 final double CjAlphaGamma   = auxiliaryElements.getAlpha() * cjsj.getdCjdGamma(j) - auxiliaryElements.getGamma() * cjsj.getdCjdAlpha(j);
  1033.                 // Cj(alpha,beta) = alpha * dC/dbeta - beta * dC/dalpha
  1034.                 final double CjAlphaBeta   = auxiliaryElements.getAlpha() * cjsj.getdCjdBeta(j) - auxiliaryElements.getBeta() * cjsj.getdCjdAlpha(j);
  1035.                 // Cj(beta,gamma) = beta * dC/dgamma - gamma * dC/dbeta
  1036.                 final double CjBetaGamma    =  auxiliaryElements.getBeta() * cjsj.getdCjdGamma(j) - auxiliaryElements.getGamma() * cjsj.getdCjdBeta(j);
  1037.                 // Cj(h,k) = h * dC/dk - k * dC/dh
  1038.                 final double CjHK   = auxiliaryElements.getH() * cjsj.getdCjdK(j) - auxiliaryElements.getK() * cjsj.getdCjdH(j);
  1039.                 // Sj(alpha,gamma) = alpha * dS/dgamma - gamma * dS/dalpha
  1040.                 final double SjAlphaGamma   = auxiliaryElements.getAlpha() * cjsj.getdSjdGamma(j) - auxiliaryElements.getGamma() * cjsj.getdSjdAlpha(j);
  1041.                 // Sj(alpha,beta) = alpha * dS/dbeta - beta * dS/dalpha
  1042.                 final double SjAlphaBeta   = auxiliaryElements.getAlpha() * cjsj.getdSjdBeta(j) - auxiliaryElements.getBeta() * cjsj.getdSjdAlpha(j);
  1043.                 // Sj(beta,gamma) = beta * dS/dgamma - gamma * dS/dbeta
  1044.                 final double SjBetaGamma    =  auxiliaryElements.getBeta() * cjsj.getdSjdGamma(j) - auxiliaryElements.getGamma() * cjsj.getdSjdBeta(j);
  1045.                 // Sj(h,k) = h * dS/dk - k * dS/dh
  1046.                 final double SjHK   = auxiliaryElements.getH() * cjsj.getdSjdK(j) - auxiliaryElements.getK() * cjsj.getdSjdH(j);

  1048.                 //Coefficients for a
  1049.                 final double coef1 = context.getX3ON2A() * (3 - context.getBB()) * j;
  1050.                 currentCij[0] += coef1 * cjsj.getSj(j);
  1051.                 currentSij[0] -= coef1 * cjsj.getCj(j);

  1053.                 //Coefficients for k and h
  1054.                 final double coef2 = auxiliaryElements.getP() * CjAlphaGamma - I * auxiliaryElements.getQ() * CjBetaGamma;
  1055.                 final double coef3 = auxiliaryElements.getP() * SjAlphaGamma - I * auxiliaryElements.getQ() * SjBetaGamma;
  1056.                 currentCij[1] -= context.getXON2A2() * (auxiliaryElements.getH() * coef2 + context.getBB() * cjsj.getdCjdH(j) - 1.5 * auxiliaryElements.getK() * j * cjsj.getSj(j));
  1057.                 currentSij[1] -= context.getXON2A2() * (auxiliaryElements.getH() * coef3 + context.getBB() * cjsj.getdSjdH(j) + 1.5 * auxiliaryElements.getK() * j * cjsj.getCj(j));
  1058.                 currentCij[2] += context.getXON2A2() * (auxiliaryElements.getK() * coef2 + context.getBB() * cjsj.getdCjdK(j) + 1.5 * auxiliaryElements.getH() * j * cjsj.getSj(j));
  1059.                 currentSij[2] += context.getXON2A2() * (auxiliaryElements.getK() * coef3 + context.getBB() * cjsj.getdSjdK(j) - 1.5 * auxiliaryElements.getH() * j * cjsj.getCj(j));

  1061.                 //Coefficients for q and p
  1062.                 final double coef4 = CjHK - CjAlphaBeta - j * cjsj.getSj(j);
  1063.                 final double coef5 = SjHK - SjAlphaBeta + j * cjsj.getCj(j);
  1064.                 currentCij[3] = context.getCXO2N2A2() * (-I * CjAlphaGamma + auxiliaryElements.getQ() * coef4);
  1065.                 currentSij[3] = context.getCXO2N2A2() * (-I * SjAlphaGamma + auxiliaryElements.getQ() * coef5);
  1066.                 currentCij[4] = context.getCXO2N2A2() * (-CjBetaGamma + auxiliaryElements.getP() * coef4);
  1067.                 currentSij[4] = context.getCXO2N2A2() * (-SjBetaGamma + auxiliaryElements.getP() * coef5);

  1069.                 //Coefficients for λ
  1070.                 final double coef6 = auxiliaryElements.getH() * cjsj.getdCjdH(j) + auxiliaryElements.getK() * cjsj.getdCjdK(j);
  1071.                 final double coef7 = auxiliaryElements.getH() * cjsj.getdSjdH(j) + auxiliaryElements.getK() * cjsj.getdSjdK(j);
  1072.                 currentCij[5] += context.getOON2A2() * (-2 * auxiliaryElements.getSma() * cjsj.getdCjdA(j) + coef6 / (context.getX() + 1) + context.getX() * coef2 - 3 * cjsj.getCj(j));
  1073.                 currentSij[5] += context.getOON2A2() * (-2 * auxiliaryElements.getSma() * cjsj.getdSjdA(j) + coef7 / (context.getX() + 1) + context.getX() * coef3 - 3 * cjsj.getSj(j));
  1074.             }

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

  1081.             // Add the coefficients to the interpolation grid
  1082.             slot.cij[j].addGridPoint(date, currentCij);
  1083.             slot.sij[j].addGridPoint(date, currentSij);

  1085.         }

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

  1090.     }

  1092.     /**
  1093.      * Generate the values for the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
  1094.      * @param <T> type of the elements
  1095.      * @param date date of computation
  1096.      * @param slot slot to which the coefficients belong
  1097.      * @param cjsj Fourier coefficients
  1098.      * @param rhoSigma ρ<sub>j</sub> and σ<sub>j</sub>
  1099.      * @param context container for attributes
  1100.      * @param auxiliaryElements auxiliary elements related to the current orbit
  1101.      * @param field field used by default
  1102.      * @param udu derivatives of the gravitational potential U
  1103.      */
  1104.     private <T extends RealFieldElement<T>> void computeCijSijCoefficients(final FieldAbsoluteDate<T> date,
  1105.                                                                            final FieldSlot<T> slot,
  1106.                                                                            final FieldFourierCjSjCoefficients<T> cjsj,
  1107.                                                                            final T[][] rhoSigma,
  1108.                                                                            final FieldDSSTZonalContext<T> context,
  1109.                                                                            final FieldAuxiliaryElements<T> auxiliaryElements,
  1110.                                                                            final Field<T> field,
  1111.                                                                            final FieldUAnddU<T> udu) {

  1113.         // Zero
  1114.         final T zero = field.getZero();

  1116.         final int nMax = maxDegreeShortPeriodics;

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

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

  1124.             // Create local arrays
  1125.             final T[] currentCij = MathArrays.buildArray(field, 6);
  1126.             final T[] currentSij = MathArrays.buildArray(field, 6);

  1128.             Arrays.fill(currentCij, zero);
  1129.             Arrays.fill(currentSij, zero);

  1131.             // j == 1
  1132.             if (j == 1) {
  1133.                 final T coef1 = auxiliaryElements.getK().multiply(udu.getU()).multiply(4.).subtract(context.getHK().multiply(cjsj.getCj(1))).add(context.getK2MH2O2().multiply(cjsj.getSj(1)));
  1134.                 final T coef2 = auxiliaryElements.getH().multiply(udu.getU()).multiply(4.).add(context.getK2MH2O2().multiply(cjsj.getCj(1))).add(context.getHK().multiply(cjsj.getSj(1)));
  1135.                 final T coef3 = auxiliaryElements.getK().multiply(cjsj.getCj(1)).add(auxiliaryElements.getH().multiply(cjsj.getSj(1))).divide(4.);
  1136.                 final T coef4 = udu.getU().multiply(8.).subtract(auxiliaryElements.getH().multiply(cjsj.getCj(1))).add(auxiliaryElements.getK().multiply(cjsj.getSj(1))).divide(4.);

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

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

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

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

  1155.             // j == 2
  1156.             if (j == 2) {
  1157.                 final T coef1 = context.getK2MH2().multiply(udu.getU());
  1158.                 final T coef2 = context.getHK().multiply(udu.getU()).multiply(2.);
  1159.                 final T coef3 = auxiliaryElements.getH().multiply(udu.getU()).divide(2.);
  1160.                 final T coef4 = auxiliaryElements.getK().multiply(udu.getU()).divide(2.);

  1162.                 //Coefficients for a
  1163.                 currentCij[0] = currentCij[0].add(coef1);
  1164.                 currentSij[0] = currentSij[0].add(coef2);

  1166.                 //Coefficients for k
  1167.                 currentCij[1] = currentCij[1].add(coef4);
  1168.                 currentSij[1] = currentSij[1].add(coef3);

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

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

  1179.             // j between 1 and 2N-3
  1180.             if (isBetween(j, 1, 2 * nMax - 3) && j + 2 < cjsj.jMax) {
  1181.                 final T coef1 = context.getHK().negate().multiply(cjsj.getCj(j + 2)).add(context.getK2MH2O2().multiply(cjsj.getSj(j + 2))).multiply(j + 2);
  1182.                 final T coef2 = context.getK2MH2O2().multiply(cjsj.getCj(j + 2)).add(context.getHK().multiply(cjsj.getSj(j + 2))).multiply(j + 2);
  1183.                 final T coef3 = auxiliaryElements.getK().multiply(cjsj.getCj(j + 2)).add(auxiliaryElements.getH().multiply(cjsj.getSj(j + 2))).multiply(j + 2).divide(4.);
  1184.                 final T coef4 = auxiliaryElements.getH().multiply(cjsj.getCj(j + 2)).subtract(auxiliaryElements.getK().multiply(cjsj.getSj(j + 2))).multiply(j + 2).divide(4.);

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

  1190.                 //Coefficients for k
  1191.                 currentCij[1] = currentCij[1].add(coef4.negate());
  1192.                 currentSij[1] = currentSij[1].subtract(coef3);

  1194.                 //Coefficients for h
  1195.                 currentCij[2] = currentCij[2].subtract(coef3);
  1196.                 currentSij[2] = currentSij[2].add(coef4);

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

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

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

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

  1218.                 //Coefficients for h
  1219.                 currentCij[2] = currentCij[2].subtract(coef3);
  1220.                 currentSij[2] = currentSij[2].subtract(coef4);

  1222.                 //Coefficients for λ
  1223.                 currentCij[5] = currentCij[5].subtract(coef2.divide(2.));
  1224.                 currentSij[5] = currentSij[5].add(coef1.negate().divide(2.));
  1225.             }

  1227.             // j between 2 and 2N
  1228.             if (isBetween(j, 2, 2 * nMax) && j - 1 < cjsj.jMax) {
  1229.                 final T coef1 = auxiliaryElements.getH().multiply(cjsj.getCj(j - 1)).add(auxiliaryElements.getK().multiply(cjsj.getSj(j - 1))).multiply(2 * ( j - 1 ));
  1230.                 final T coef2 = auxiliaryElements.getK().multiply(cjsj.getCj(j - 1)).subtract(auxiliaryElements.getH().multiply(cjsj.getSj(j - 1))).multiply(2 * ( j - 1 ));
  1231.                 final T coef3 = cjsj.getCj(j - 1).multiply(j - 1);
  1232.                 final T coef4 = cjsj.getSj(j - 1).multiply(j - 1);

  1234.                 //Coefficients for a
  1235.                 currentCij[0] = currentCij[0].add(coef1);
  1236.                 currentSij[0] = currentSij[0].subtract(coef2);

  1238.                 //Coefficients for k
  1239.                 currentCij[1] = currentCij[1].add(coef4);
  1240.                 currentSij[1] = currentSij[1].subtract(coef3);

  1242.                 //Coefficients for h
  1243.                 currentCij[2] = currentCij[2].add(coef3);
  1244.                 currentSij[2] = currentSij[2].add(coef4);

  1246.                 //Coefficients for λ
  1247.                 currentCij[5] = currentCij[5].add(coef2.divide(2.));
  1248.                 currentSij[5] = currentSij[5].add(coef1.divide(2.));
  1249.             }

  1251.             // j between 3 and 2N + 1
  1252.             if (isBetween(j, 3, 2 * nMax + 1) && j - 2 < cjsj.jMax) {
  1253.                 final T coef1 = context.getHK().multiply(cjsj.getCj(j - 2)).add(context.getK2MH2O2().multiply(cjsj.getSj(j - 2))).multiply(j - 2);
  1254.                 final T coef2 = context.getK2MH2O2().negate().multiply(cjsj.getCj(j - 2)).add(context.getHK().multiply(cjsj.getSj(j - 2))).multiply(j - 2);
  1255.                 final T coef3 = auxiliaryElements.getK().multiply(cjsj.getCj(j - 2)).subtract(auxiliaryElements.getH().multiply(cjsj.getSj(j - 2))).multiply(j - 2).divide(4.);
  1256.                 final T coef4 = auxiliaryElements.getH().multiply(cjsj.getCj(j - 2)).add(auxiliaryElements.getK().multiply(cjsj.getSj(j - 2))).multiply(j - 2).divide(4.);
  1257.                 final T coef5 = context.getK2MH2O2().multiply(cjsj.getCj(j - 2)).subtract(context.getHK().multiply(cjsj.getSj(j - 2))).multiply(j - 2);

  1259.                 //Coefficients for a
  1260.                 currentCij[0] = currentCij[0].add(coef1);
  1261.                 currentSij[0] = currentSij[0].add(coef2);

  1263.                 //Coefficients for k
  1264.                 currentCij[1] = currentCij[1].add(coef4);
  1265.                 currentSij[1] = currentSij[1].add(coef3.negate());

  1267.                 //Coefficients for h
  1268.                 currentCij[2] = currentCij[2].add(coef3);
  1269.                 currentSij[2] = currentSij[2].add(coef4);

  1271.                 //Coefficients for λ
  1272.                 currentCij[5] = currentCij[5].add(coef5.divide(2.));
  1273.                 currentSij[5] = currentSij[5].add(coef1.divide(2.));
  1274.             }

  1276.             //multiply by the common factor
  1277.             //for a (i == 0) -> χ³ / (n² * a)
  1278.             currentCij[0] = currentCij[0].multiply(context.getX3ON2A());
  1279.             currentSij[0] = currentSij[0].multiply(context.getX3ON2A());
  1280.             //for k (i == 1) -> χ / (n² * a²)
  1281.             currentCij[1] = currentCij[1].multiply(context.getXON2A2());
  1282.             currentSij[1] = currentSij[1].multiply(context.getXON2A2());
  1283.             //for h (i == 2) -> χ / (n² * a²)
  1284.             currentCij[2] = currentCij[2].multiply(context.getXON2A2());
  1285.             currentSij[2] = currentSij[2].multiply(context.getXON2A2());
  1286.             //for λ (i == 5) -> (χ²) / (n² * a² * (χ + 1 ) )
  1287.             currentCij[5] = currentCij[5].multiply(context.getX2ON2A2XP1());
  1288.             currentSij[5] = currentSij[5].multiply(context.getX2ON2A2XP1());

  1290.             // j is between 1 and 2 * N - 1
  1291.             if (isBetween(j, 1, 2 * nMax - 1) && j < cjsj.jMax) {
  1292.                 // Compute cross derivatives
  1293.                 // Cj(alpha,gamma) = alpha * dC/dgamma - gamma * dC/dalpha
  1294.                 final T CjAlphaGamma = auxiliaryElements.getAlpha().multiply(cjsj.getdCjdGamma(j)).subtract(auxiliaryElements.getGamma().multiply(cjsj.getdCjdAlpha(j)));
  1295.                 // Cj(alpha,beta) = alpha * dC/dbeta - beta * dC/dalpha
  1296.                 final T CjAlphaBeta  = auxiliaryElements.getAlpha().multiply(cjsj.getdCjdBeta(j)).subtract(auxiliaryElements.getBeta().multiply(cjsj.getdCjdAlpha(j)));
  1297.                 // Cj(beta,gamma) = beta * dC/dgamma - gamma * dC/dbeta
  1298.                 final T CjBetaGamma  = auxiliaryElements.getBeta().multiply(cjsj.getdCjdGamma(j)).subtract(auxiliaryElements.getGamma().multiply(cjsj.getdCjdBeta(j)));
  1299.                 // Cj(h,k) = h * dC/dk - k * dC/dh
  1300.                 final T CjHK         = auxiliaryElements.getH().multiply(cjsj.getdCjdK(j)).subtract(auxiliaryElements.getK().multiply(cjsj.getdCjdH(j)));
  1301.                 // Sj(alpha,gamma) = alpha * dS/dgamma - gamma * dS/dalpha
  1302.                 final T SjAlphaGamma = auxiliaryElements.getAlpha().multiply(cjsj.getdSjdGamma(j)).subtract(auxiliaryElements.getGamma().multiply(cjsj.getdSjdAlpha(j)));
  1303.                 // Sj(alpha,beta) = alpha * dS/dbeta - beta * dS/dalpha
  1304.                 final T SjAlphaBeta  = auxiliaryElements.getAlpha().multiply(cjsj.getdSjdBeta(j)).subtract(auxiliaryElements.getBeta().multiply(cjsj.getdSjdAlpha(j)));
  1305.                 // Sj(beta,gamma) = beta * dS/dgamma - gamma * dS/dbeta
  1306.                 final T SjBetaGamma  = auxiliaryElements.getBeta().multiply(cjsj.getdSjdGamma(j)).subtract(auxiliaryElements.getGamma().multiply(cjsj.getdSjdBeta(j)));
  1307.                 // Sj(h,k) = h * dS/dk - k * dS/dh
  1308.                 final T SjHK         = auxiliaryElements.getH().multiply(cjsj.getdSjdK(j)).subtract(auxiliaryElements.getK().multiply(cjsj.getdSjdH(j)));

  1310.                 //Coefficients for a
  1311.                 final T coef1 = context.getX3ON2A().multiply(context.getBB().negate().add(3.)).multiply(j);
  1312.                 currentCij[0] = currentCij[0].add(coef1.multiply(cjsj.getSj(j)));
  1313.                 currentSij[0] = currentSij[0].subtract(coef1.multiply(cjsj.getCj(j)));

  1315.                 //Coefficients for k and h
  1316.                 final T coef2 = auxiliaryElements.getP().multiply(CjAlphaGamma).subtract(auxiliaryElements.getQ().multiply(CjBetaGamma).multiply(I));
  1317.                 final T coef3 = auxiliaryElements.getP().multiply(SjAlphaGamma).subtract(auxiliaryElements.getQ().multiply(SjBetaGamma).multiply(I));
  1318.                 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)))));
  1319.                 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)))));
  1320.                 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)))));
  1321.                 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)))));

  1323.                 //Coefficients for q and p
  1324.                 final T coef4 = CjHK.subtract(CjAlphaBeta).subtract(cjsj.getSj(j).multiply(j));
  1325.                 final T coef5 = SjHK.subtract(SjAlphaBeta).add(cjsj.getCj(j).multiply(j));
  1326.                 currentCij[3] = context.getCXO2N2A2().multiply(CjAlphaGamma.multiply(-I).add(auxiliaryElements.getQ().multiply(coef4)));
  1327.                 currentSij[3] = context.getCXO2N2A2().multiply(SjAlphaGamma.multiply(-I).add(auxiliaryElements.getQ().multiply(coef5)));
  1328.                 currentCij[4] = context.getCXO2N2A2().multiply(CjBetaGamma.negate().add(auxiliaryElements.getP().multiply(coef4)));
  1329.                 currentSij[4] = context.getCXO2N2A2().multiply(SjBetaGamma.negate().add(auxiliaryElements.getP().multiply(coef5)));

  1331.                 //Coefficients for λ
  1332.                 final T coef6 = auxiliaryElements.getH().multiply(cjsj.getdCjdH(j)).add(auxiliaryElements.getK().multiply(cjsj.getdCjdK(j)));
  1333.                 final T coef7 = auxiliaryElements.getH().multiply(cjsj.getdSjdH(j)).add(auxiliaryElements.getK().multiply(cjsj.getdSjdK(j)));
  1334.                 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.))));
  1335.                 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.))));
  1336.             }

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

  1343.             // Add the coefficients to the interpolation grid
  1344.             slot.cij[j].addGridPoint(date, currentCij);
  1345.             slot.sij[j].addGridPoint(date, currentSij);

  1347.         }

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

  1352.     }

  1354.     /**
  1355.      * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
  1356.      * <p>
  1357.      * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
  1358.      *  ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
  1359.      *  σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
  1360.      * </p>
  1361.      * @param date target date
  1362.      * @param slot slot to which the coefficients belong
  1363.      * @param auxiliaryElements auxiliary elements related to the current orbit
  1364.      * @return array containing ρ<sub>j</sub> and σ<sub>j</sub>
  1365.      */
  1366.     private double[][] computeRhoSigmaCoefficients(final AbsoluteDate date,
  1367.                                                    final Slot slot,
  1368.                                                    final AuxiliaryElements auxiliaryElements) {

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

  1373.         // (-b)<sup>j</sup>
  1374.         double mbtj = 1;

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

  1379.             //Compute current rho and sigma;
  1380.             mbtj *= -b;
  1381.             final double coef  = (1 + j * auxiliaryElements.getB()) * mbtj;
  1382.             final double rho   = coef * cjsjKH.getCj(j);
  1383.             final double sigma = coef * cjsjKH.getSj(j);

  1385.             // Add the coefficients to the interpolation grid
  1386.             rhoSigma[j][0] = rho;
  1387.             rhoSigma[j][1] = sigma;
  1388.         }

  1390.         return rhoSigma;

  1392.     }

  1394.     /**
  1395.      * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
  1396.      * <p>
  1397.      * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
  1398.      *  ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
  1399.      *  σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
  1400.      * </p>
  1401.      * @param <T> type of the elements
  1402.      * @param date target date
  1403.      * @param slot slot to which the coefficients belong
  1404.      * @param auxiliaryElements auxiliary elements related to the current orbit
  1405.      * @param field field used by default
  1406.      * @return array containing ρ<sub>j</sub> and σ<sub>j</sub>
  1407.      */
  1408.     private <T extends RealFieldElement<T>> T[][] computeRhoSigmaCoefficients(final FieldAbsoluteDate<T> date,
  1409.                                                                               final FieldSlot<T> slot,
  1410.                                                                               final FieldAuxiliaryElements<T> auxiliaryElements,
  1411.                                                                               final Field<T> field) {
  1412.         final T zero = field.getZero();

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

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

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

  1423.             //Compute current rho and sigma;
  1424.             mbtj = mbtj.multiply(b.negate());
  1425.             final T coef  = mbtj.multiply(auxiliaryElements.getB().multiply(j).add(1.));
  1426.             final T rho   = coef.multiply(cjsjKH.getCj(j));
  1427.             final T sigma = coef.multiply(cjsjKH.getSj(j));

  1429.             // Add the coefficients to the interpolation grid
  1430.             rhoSigma[j][0] = rho;
  1431.             rhoSigma[j][1] = sigma;
  1432.         }

  1434.         return rhoSigma;

  1436.     }

  1438.     /** The coefficients used to compute the short-periodic zonal contribution.
  1439.      *
  1440.      * <p>
  1441.      * Those coefficients are given in Danielson paper by expressions 4.1-(20) to 4.1.-(25)
  1442.      * </p>
  1443.      * <p>
  1444.      * The coefficients are: <br>
  1445.      * - C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> <br>
  1446.      * - ρ<sub>j</sub> and σ<sub>j</sub> <br>
  1447.      * - C<sub>i</sub>⁰
  1448.      * </p>
  1449.      *
  1450.      * @author Lucian Barbulescu
  1451.      */
  1452.     private static class ZonalShortPeriodicCoefficients implements ShortPeriodTerms {

  1454.         /** Maximum value for j index. */
  1455.         private final int maxFrequencyShortPeriodics;

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

  1460.         /** All coefficients slots. */
  1461.         private final transient TimeSpanMap<Slot> slots;

  1463.         /** Constructor.
  1464.          * @param maxFrequencyShortPeriodics maximum value for j index
  1465.          * @param interpolationPoints number of points used in the interpolation process
  1466.          * @param slots all coefficients slots
  1467.          */
  1468.         ZonalShortPeriodicCoefficients(final int maxFrequencyShortPeriodics, final int interpolationPoints,
  1469.                                        final TimeSpanMap<Slot> slots) {

  1471.             // Save parameters
  1472.             this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
  1473.             this.interpolationPoints        = interpolationPoints;
  1474.             this.slots                      = slots;

  1476.         }

  1478.         /** Get the slot valid for some date.
  1479.          * @param meanStates mean states defining the slot
  1480.          * @return slot valid at the specified date
  1481.          */
  1482.         public Slot createSlot(final SpacecraftState... meanStates) {
  1483.             final Slot         slot  = new Slot(maxFrequencyShortPeriodics, interpolationPoints);
  1484.             final AbsoluteDate first = meanStates[0].getDate();
  1485.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  1486.             if (first.compareTo(last) <= 0) {
  1487.                 slots.addValidAfter(slot, first);
  1488.             } else {
  1489.                 slots.addValidBefore(slot, first);
  1490.             }
  1491.             return slot;
  1492.         }

  1494.         /** {@inheritDoc} */
  1495.         @Override
  1496.         public double[] value(final Orbit meanOrbit) {

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

  1501.             // Get the True longitude L
  1502.             final double L = meanOrbit.getLv();

  1504.             // Compute the center
  1505.             final double center = L - meanOrbit.getLM();

  1507.             // Initialize short periodic variations
  1508.             final double[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  1509.             final double[] d = slot.di.value(meanOrbit.getDate());
  1510.             for (int i = 0; i < 6; i++) {
  1511.                 shortPeriodicVariation[i] +=  center * d[i];
  1512.             }

  1514.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1515.                 final double[] c = slot.cij[j].value(meanOrbit.getDate());
  1516.                 final double[] s = slot.sij[j].value(meanOrbit.getDate());
  1517.                 final SinCos sc  = FastMath.sinCos(j * L);
  1518.                 for (int i = 0; i < 6; i++) {
  1519.                     // add corresponding term to the short periodic variation
  1520.                     shortPeriodicVariation[i] += c[i] * sc.cos();
  1521.                     shortPeriodicVariation[i] += s[i] * sc.sin();
  1522.                 }
  1523.             }

  1525.             return shortPeriodicVariation;
  1526.         }

  1528.         /** {@inheritDoc} */
  1529.         @Override
  1530.         public String getCoefficientsKeyPrefix() {
  1531.             return DSSTZonal.SHORT_PERIOD_PREFIX;
  1532.         }

  1534.         /** {@inheritDoc}
  1535.          * <p>
  1536.          * For zonal terms contributions,there are maxJ cj coefficients,
  1537.          * maxJ sj coefficients and 2 dj coefficients, where maxJ depends
  1538.          * on the orbit. The j index is the integer multiplier for the true
  1539.          * longitude argument in the cj and sj coefficients and the degree
  1540.          * in the polynomial dj coefficients.
  1541.          * </p>
  1542.          */
  1543.         @Override
  1544.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected) {

  1546.             // select the coefficients slot
  1547.             final Slot slot = slots.get(date);

  1549.             final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * maxFrequencyShortPeriodics + 2);
  1550.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  1551.             storeIfSelected(coefficients, selected, slot.di.value(date), "d", 1);
  1552.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1553.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  1554.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  1555.             }
  1556.             return coefficients;

  1558.         }

  1560.         /** Put a coefficient in a map if selected.
  1561.          * @param map map to populate
  1562.          * @param selected set of coefficients that should be put in the map
  1563.          * (empty set means all coefficients are selected)
  1564.          * @param value coefficient value
  1565.          * @param id coefficient identifier
  1566.          * @param indices list of coefficient indices
  1567.          */
  1568.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
  1569.                                      final double[] value, final String id, final int... indices) {
  1570.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  1571.             keyBuilder.append(id);
  1572.             for (int index : indices) {
  1573.                 keyBuilder.append('[').append(index).append(']');
  1574.             }
  1575.             final String key = keyBuilder.toString();
  1576.             if (selected.isEmpty() || selected.contains(key)) {
  1577.                 map.put(key, value);
  1578.             }
  1579.         }

  1581.     }

  1583.     /** The coefficients used to compute the short-periodic zonal contribution.
  1584.     *
  1585.     * <p>
  1586.     * Those coefficients are given in Danielson paper by expressions 4.1-(20) to 4.1.-(25)
  1587.     * </p>
  1588.     * <p>
  1589.     * The coefficients are: <br>
  1590.     * - C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> <br>
  1591.     * - ρ<sub>j</sub> and σ<sub>j</sub> <br>
  1592.     * - C<sub>i</sub>⁰
  1593.     * </p>
  1594.     *
  1595.     * @author Lucian Barbulescu
  1596.     */
  1597.     private static class FieldZonalShortPeriodicCoefficients <T extends RealFieldElement<T>> implements FieldShortPeriodTerms<T> {

  1599.         /** Maximum value for j index. */
  1600.         private final int maxFrequencyShortPeriodics;

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

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

  1608.        /** Constructor.
  1609.         * @param maxFrequencyShortPeriodics maximum value for j index
  1610.         * @param interpolationPoints number of points used in the interpolation process
  1611.         * @param slots all coefficients slots
  1612.         */
  1613.         FieldZonalShortPeriodicCoefficients(final int maxFrequencyShortPeriodics, final int interpolationPoints,
  1614.                                             final FieldTimeSpanMap<FieldSlot<T>, T> slots) {

  1616.             // Save parameters
  1617.             this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
  1618.             this.interpolationPoints        = interpolationPoints;
  1619.             this.slots                      = slots;

  1621.         }

  1623.        /** Get the slot valid for some date.
  1624.         * @param meanStates mean states defining the slot
  1625.         * @return slot valid at the specified date
  1626.         */
  1627.         @SuppressWarnings("unchecked")
  1628.         public FieldSlot<T> createSlot(final FieldSpacecraftState<T>... meanStates) {
  1629.             final FieldSlot<T>         slot  = new FieldSlot<>(maxFrequencyShortPeriodics, interpolationPoints);
  1630.             final FieldAbsoluteDate<T> first = meanStates[0].getDate();
  1631.             final FieldAbsoluteDate<T> last  = meanStates[meanStates.length - 1].getDate();
  1632.             if (first.compareTo(last) <= 0) {
  1633.                 slots.addValidAfter(slot, first);
  1634.             } else {
  1635.                 slots.addValidBefore(slot, first);
  1636.             }
  1637.             return slot;
  1638.         }

  1640.         /** {@inheritDoc} */
  1641.         @Override
  1642.         public T[] value(final FieldOrbit<T> meanOrbit) {

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

  1647.             // Get the True longitude L
  1648.             final T L = meanOrbit.getLv();

  1650.             // Compute the center
  1651.             final T center = L.subtract(meanOrbit.getLM());

  1653.             // Initialize short periodic variations
  1654.             final T[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  1655.             final T[] d = slot.di.value(meanOrbit.getDate());
  1656.             for (int i = 0; i < 6; i++) {
  1657.                 shortPeriodicVariation[i] = shortPeriodicVariation[i].add(center.multiply(d[i]));
  1658.             }

  1660.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1661.                 final T[]            c   = slot.cij[j].value(meanOrbit.getDate());
  1662.                 final T[]            s   = slot.sij[j].value(meanOrbit.getDate());
  1663.                 final FieldSinCos<T> sc  = FastMath.sinCos(L.multiply(j));
  1664.                 for (int i = 0; i < 6; i++) {
  1665.                     // add corresponding term to the short periodic variation
  1666.                     shortPeriodicVariation[i] = shortPeriodicVariation[i].add(c[i].multiply(sc.cos()));
  1667.                     shortPeriodicVariation[i] = shortPeriodicVariation[i].add(s[i].multiply(sc.sin()));
  1668.                 }
  1669.             }

  1671.             return shortPeriodicVariation;
  1672.         }

  1674.         /** {@inheritDoc} */
  1675.         @Override
  1676.         public String getCoefficientsKeyPrefix() {
  1677.             return DSSTZonal.SHORT_PERIOD_PREFIX;
  1678.         }

  1680.        /** {@inheritDoc}
  1681.         * <p>
  1682.         * For zonal terms contributions,there are maxJ cj coefficients,
  1683.         * maxJ sj coefficients and 2 dj coefficients, where maxJ depends
  1684.         * on the orbit. The j index is the integer multiplier for the true
  1685.         * longitude argument in the cj and sj coefficients and the degree
  1686.         * in the polynomial dj coefficients.
  1687.         * </p>
  1688.         */
  1689.         @Override
  1690.         public Map<String, T[]> getCoefficients(final FieldAbsoluteDate<T> date, final Set<String> selected) {

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

  1695.             final Map<String, T[]> coefficients = new HashMap<String, T[]>(2 * maxFrequencyShortPeriodics + 2);
  1696.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  1697.             storeIfSelected(coefficients, selected, slot.di.value(date), "d", 1);
  1698.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1699.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  1700.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  1701.             }
  1702.             return coefficients;

  1704.         }

  1706.        /** Put a coefficient in a map if selected.
  1707.         * @param map map to populate
  1708.         * @param selected set of coefficients that should be put in the map
  1709.         * (empty set means all coefficients are selected)
  1710.         * @param value coefficient value
  1711.         * @param id coefficient identifier
  1712.         * @param indices list of coefficient indices
  1713.         */
  1714.         private void storeIfSelected(final Map<String, T[]> map, final Set<String> selected,
  1715.                                      final T[] value, final String id, final int... indices) {
  1716.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  1717.             keyBuilder.append(id);
  1718.             for (int index : indices) {
  1719.                 keyBuilder.append('[').append(index).append(']');
  1720.             }
  1721.             final String key = keyBuilder.toString();
  1722.             if (selected.isEmpty() || selected.contains(key)) {
  1723.                 map.put(key, value);
  1724.             }
  1725.         }

  1727.     }

  1729.     /** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
  1730.      *  <p>
  1731.      *  Those coefficients are given in Danielson paper by expressions 4.1-(13) to 4.1.-(16b)
  1732.      *  </p>
  1733.      */
  1734.     private class FourierCjSjCoefficients {

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

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

  1742.         /** Maximum possible value for n. */
  1743.         private final int nMax;

  1745.         /** Maximum possible value for s. */
  1746.         private final int sMax;

  1748.         /** Maximum possible value for j. */
  1749.         private final int jMax;

  1751.         /** The C<sup>j</sup> coefficients and their derivatives.
  1752.          * <p>
  1753.          * Each column of the matrix contains the following values: <br/>
  1754.          * - C<sup>j</sup> <br/>
  1755.          * - dC<sup>j</sup> / da <br/>
  1756.          * - dC<sup>j</sup> / dh <br/>
  1757.          * - dC<sup>j</sup> / dk <br/>
  1758.          * - dC<sup>j</sup> / dα <br/>
  1759.          * - dC<sup>j</sup> / dβ <br/>
  1760.          * - dC<sup>j</sup> / dγ <br/>
  1761.          * </p>
  1762.          */
  1763.         private final double[][] cCoef;

  1765.         /** The S<sup>j</sup> coefficients and their derivatives.
  1766.          * <p>
  1767.          * Each column of the matrix contains the following values: <br/>
  1768.          * - S<sup>j</sup> <br/>
  1769.          * - dS<sup>j</sup> / da <br/>
  1770.          * - dS<sup>j</sup> / dh <br/>
  1771.          * - dS<sup>j</sup> / dk <br/>
  1772.          * - dS<sup>j</sup> / dα <br/>
  1773.          * - dS<sup>j</sup> / dβ <br/>
  1774.          * - dS<sup>j</sup> / dγ <br/>
  1775.          * </p>
  1776.          */
  1777.         private final double[][] sCoef;

  1779.         /** h * &Chi;³. */
  1780.         private final double hXXX;
  1781.         /** k * &Chi;³. */
  1782.         private final double kXXX;

  1784.         /** Create a set of C<sup>j</sup> and the S<sup>j</sup> coefficients.
  1785.          *  @param date the current date
  1786.          *  @param nMax maximum possible value for n
  1787.          *  @param sMax maximum possible value for s
  1788.          *  @param jMax maximum possible value for j
  1789.          *  @param context container for attributes
  1790.          *  @param hansenObjects initialization of hansen objects
  1791.          */
  1792.         FourierCjSjCoefficients(final AbsoluteDate date,
  1793.                                 final int nMax, final int sMax, final int jMax, final DSSTZonalContext context,
  1794.                                 final HansenObjects hansenObjects) {

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

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

  1802.             this.lnsCoef = new LnsCoefficients(nMax, nMax, Qns, Vns, context.getRoa());
  1803.             this.nMax = nMax;
  1804.             this.sMax = sMax;
  1805.             this.jMax = jMax;

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

  1811.             this.cCoef = new double[7][jMax + 1];
  1812.             this.sCoef = new double[7][jMax + 1];

  1814.             for (int s = 0; s <= sMax; s++) {
  1815.                 //Initialise the Hansen roots
  1816.                 hansenObjects.computeHansenObjectsInitValues(context, s);
  1817.             }
  1818.             generateCoefficients(date, context, auxiliaryElements, hansenObjects);
  1819.         }

  1821.         /** Generate all coefficients.
  1822.          * @param date the current date
  1823.          * @param context container for attributes
  1824.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1825.          * @param hansenObjects initialization of hansen objects
  1826.          */
  1827.         private void generateCoefficients(final AbsoluteDate date,
  1828.                                           final DSSTZonalContext context,
  1829.                                           final AuxiliaryElements auxiliaryElements,
  1830.                                           final HansenObjects hansenObjects) {

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

  1835.                 //init arrays
  1836.                 for (int i = 0; i <= 6; i++) {
  1837.                     cCoef[i][j] = 0.;
  1838.                     sCoef[i][j] = 0.;
  1839.                 }

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

  1843.                     //compute first double sum where s: j -> N-1 and n: s+1 -> N
  1844.                     for (int s = j; s <= FastMath.min(nMax - 1, sMax); s++) {
  1845.                         // j - s
  1846.                         final int jms = j - s;
  1847.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  1848.                         final int d0smj = (s == j) ? 1 : 2;

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

  1857.                                 final double hjs = ghijCoef.getHjs(s, -jms);
  1858.                                 final double dHjsdh = ghijCoef.getdHjsdh(s, -jms);
  1859.                                 final double dHjsdk = ghijCoef.getdHjsdk(s, -jms);
  1860.                                 final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, -jms);
  1861.                                 final double dHjsdBeta = ghijCoef.getdHjsdBeta(s, -jms);

  1863.                                 final double gjs = ghijCoef.getGjs(s, -jms);
  1864.                                 final double dGjsdh = ghijCoef.getdGjsdh(s, -jms);
  1865.                                 final double dGjsdk = ghijCoef.getdGjsdk(s, -jms);
  1866.                                 final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, -jms);
  1867.                                 final double dGjsdBeta = ghijCoef.getdGjsdBeta(s, -jms);

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

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

  1876.                                 final double coef0 = d0smj * jn;
  1877.                                 final double coef1 = coef0 * lns;
  1878.                                 final double coef2 = coef1 * kns;
  1879.                                 final double coef3 = coef2 * hjs;
  1880.                                 final double coef4 = coef2 * gjs;

  1882.                                 // Add the term to the coefficients
  1883.                                 cCoef[0][j] += coef3;
  1884.                                 cCoef[1][j] += coef3 * (n + 1);
  1885.                                 cCoef[2][j] += coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
  1886.                                 cCoef[3][j] += coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
  1887.                                 cCoef[4][j] += coef2 * dHjsdAlpha;
  1888.                                 cCoef[5][j] += coef2 * dHjsdBeta;
  1889.                                 cCoef[6][j] += coef0 * dlns * kns * hjs;

  1891.                                 sCoef[0][j] += coef4;
  1892.                                 sCoef[1][j] += coef4 * (n + 1);
  1893.                                 sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
  1894.                                 sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
  1895.                                 sCoef[4][j] += coef2 * dGjsdAlpha;
  1896.                                 sCoef[5][j] += coef2 * dGjsdBeta;
  1897.                                 sCoef[6][j] += coef0 * dlns * kns * gjs;
  1898.                             }
  1899.                         }
  1900.                     }

  1902.                     //compute second double sum where s: 0 -> N-j and n: max(j+s, j+1) -> N
  1903.                     for (int s = 0; s <= FastMath.min(nMax - j, sMax); s++) {
  1904.                         // j + s
  1905.                         final int jps = j + s;
  1906.                         // Kronecker symbols (2 - delta(0,j+s))
  1907.                         final double d0spj = (s == -j) ? 1 : 2;

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

  1916.                                 final double hjs = ghijCoef.getHjs(s, jps);
  1917.                                 final double dHjsdh = ghijCoef.getdHjsdh(s, jps);
  1918.                                 final double dHjsdk = ghijCoef.getdHjsdk(s, jps);
  1919.                                 final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, jps);
  1920.                                 final double dHjsdBeta = ghijCoef.getdHjsdBeta(s, jps);

  1922.                                 final double gjs = ghijCoef.getGjs(s, jps);
  1923.                                 final double dGjsdh = ghijCoef.getdGjsdh(s, jps);
  1924.                                 final double dGjsdk = ghijCoef.getdGjsdk(s, jps);
  1925.                                 final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, jps);
  1926.                                 final double dGjsdBeta = ghijCoef.getdGjsdBeta(s, jps);

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

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

  1935.                                 final double coef0 = d0spj * jn;
  1936.                                 final double coef1 = coef0 * lns;
  1937.                                 final double coef2 = coef1 * kns;

  1939.                                 final double coef3 = coef2 * hjs;
  1940.                                 final double coef4 = coef2 * gjs;

  1942.                                 // Add the term to the coefficients
  1943.                                 cCoef[0][j] -= coef3;
  1944.                                 cCoef[1][j] -= coef3 * (n + 1);
  1945.                                 cCoef[2][j] -= coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
  1946.                                 cCoef[3][j] -= coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
  1947.                                 cCoef[4][j] -= coef2 * dHjsdAlpha;
  1948.                                 cCoef[5][j] -= coef2 * dHjsdBeta;
  1949.                                 cCoef[6][j] -= coef0 * dlns * kns * hjs;

  1951.                                 sCoef[0][j] += coef4;
  1952.                                 sCoef[1][j] += coef4 * (n + 1);
  1953.                                 sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
  1954.                                 sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
  1955.                                 sCoef[4][j] += coef2 * dGjsdAlpha;
  1956.                                 sCoef[5][j] += coef2 * dGjsdBeta;
  1957.                                 sCoef[6][j] += coef0 * dlns * kns * gjs;
  1958.                             }
  1959.                         }
  1960.                     }

  1962.                     //compute third double sum where s: 1 -> j and  n: j+1 -> N
  1963.                     for (int s = 1; s <= FastMath.min(j, sMax); s++) {
  1964.                         // j - s
  1965.                         final int jms = j - s;
  1966.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  1967.                         final int d0smj = (s == j) ? 1 : 2;

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

  1976.                                 final double ijs = ghijCoef.getIjs(s, jms);
  1977.                                 final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  1978.                                 final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  1979.                                 final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  1980.                                 final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  1982.                                 final double jjs = ghijCoef.getJjs(s, jms);
  1983.                                 final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  1984.                                 final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  1985.                                 final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  1986.                                 final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

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

  1995.                                 final double coef0 = d0smj * jn;
  1996.                                 final double coef1 = coef0 * lns;
  1997.                                 final double coef2 = coef1 * kns;

  1999.                                 final double coef3 = coef2 * ijs;
  2000.                                 final double coef4 = coef2 * jjs;

  2002.                                 // Add the term to the coefficients
  2003.                                 cCoef[0][j] -= coef3;
  2004.                                 cCoef[1][j] -= coef3 * (n + 1);
  2005.                                 cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2006.                                 cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2007.                                 cCoef[4][j] -= coef2 * dIjsdAlpha;
  2008.                                 cCoef[5][j] -= coef2 * dIjsdBeta;
  2009.                                 cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2011.                                 sCoef[0][j] += coef4;
  2012.                                 sCoef[1][j] += coef4 * (n + 1);
  2013.                                 sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2014.                                 sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2015.                                 sCoef[4][j] += coef2 * dJjsdAlpha;
  2016.                                 sCoef[5][j] += coef2 * dJjsdBeta;
  2017.                                 sCoef[6][j] += coef0 * dlns * kns * jjs;
  2018.                             }
  2019.                         }
  2020.                     }
  2021.                 }

  2023.                 if (isBetween(j, 2, nMax)) {
  2024.                     //add first term
  2025.                     // J<sub>j</sub>
  2026.                     final double jj = -harmonics.getUnnormalizedCnm(j, 0);
  2027.                     double kns = hansenObjects.getHansenObjects()[0].getValue(-j - 1, context.getX());
  2028.                     double dkns = hansenObjects.getHansenObjects()[0].getDerivative(-j - 1, context.getX());

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

  2033.                     final double hjs = ghijCoef.getHjs(0, j);
  2034.                     final double dHjsdh = ghijCoef.getdHjsdh(0, j);
  2035.                     final double dHjsdk = ghijCoef.getdHjsdk(0, j);
  2036.                     final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(0, j);
  2037.                     final double dHjsdBeta = ghijCoef.getdHjsdBeta(0, j);

  2039.                     final double gjs = ghijCoef.getGjs(0, j);
  2040.                     final double dGjsdh = ghijCoef.getdGjsdh(0, j);
  2041.                     final double dGjsdk = ghijCoef.getdGjsdk(0, j);
  2042.                     final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(0, j);
  2043.                     final double dGjsdBeta = ghijCoef.getdGjsdBeta(0, j);

  2045.                     // 2 * J<sub>j</sub> * K₀<sup>-j-1,0</sup> * L<sub>j</sub><sup>j</sup>
  2046.                     double coef0 = 2 * jj;
  2047.                     double coef1 = coef0 * lns;
  2048.                     double coef2 = coef1 * kns;

  2050.                     double coef3 = coef2 * hjs;
  2051.                     double coef4 = coef2 * gjs;

  2053.                     // Add the term to the coefficients
  2054.                     cCoef[0][j] -= coef3;
  2055.                     cCoef[1][j] -= coef3 * (j + 1);
  2056.                     cCoef[2][j] -= coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
  2057.                     cCoef[3][j] -= coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
  2058.                     cCoef[4][j] -= coef2 * dHjsdAlpha;
  2059.                     cCoef[5][j] -= coef2 * dHjsdBeta;
  2060.                     //no contribution to cCoef[6][j] because dlns is 0

  2062.                     sCoef[0][j] += coef4;
  2063.                     sCoef[1][j] += coef4 * (j + 1);
  2064.                     sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
  2065.                     sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
  2066.                     sCoef[4][j] += coef2 * dGjsdAlpha;
  2067.                     sCoef[5][j] += coef2 * dGjsdBeta;
  2068.                     //no contribution to sCoef[6][j] because dlns is 0

  2070.                     //compute simple sum where s: 1 -> j-1
  2071.                     for (int s = 1; s <= FastMath.min(j - 1, sMax); s++) {
  2072.                         // j - s
  2073.                         final int jms = j - s;
  2074.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2075.                         final int d0smj = (s == j) ? 1 : 2;

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

  2083.                             lns = lnsCoef.getLns(j, jms);
  2084.                             final double dlns = lnsCoef.getdLnsdGamma(j, jms);

  2086.                             final double ijs = ghijCoef.getIjs(s, jms);
  2087.                             final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2088.                             final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2089.                             final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2090.                             final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2092.                             final double jjs = ghijCoef.getJjs(s, jms);
  2093.                             final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2094.                             final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2095.                             final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2096.                             final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

  2098.                             coef0 = d0smj * jj;
  2099.                             coef1 = coef0 * lns;
  2100.                             coef2 = coef1 * kns;

  2102.                             coef3 = coef2 * ijs;
  2103.                             coef4 = coef2 * jjs;

  2105.                             // Add the term to the coefficients
  2106.                             cCoef[0][j] -= coef3;
  2107.                             cCoef[1][j] -= coef3 * (j + 1);
  2108.                             cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2109.                             cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2110.                             cCoef[4][j] -= coef2 * dIjsdAlpha;
  2111.                             cCoef[5][j] -= coef2 * dIjsdBeta;
  2112.                             cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2114.                             sCoef[0][j] += coef4;
  2115.                             sCoef[1][j] += coef4 * (j + 1);
  2116.                             sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2117.                             sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2118.                             sCoef[4][j] += coef2 * dJjsdAlpha;
  2119.                             sCoef[5][j] += coef2 * dJjsdBeta;
  2120.                             sCoef[6][j] += coef0 * dlns * kns * jjs;
  2121.                         }
  2122.                     }
  2123.                 }

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

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

  2131.                     //if j is even
  2132.                     if (j % 2 == 0) {
  2133.                         //compute first double sum where s: j-min(j-1,N) -> j/2-1 and n: j-s -> min(j-1,N)
  2134.                         for (int s = j - minjm1on; s <= FastMath.min(j / 2 - 1, sMax); s++) {
  2135.                             // j - s
  2136.                             final int jms = j - s;
  2137.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2138.                             final int d0smj = (s == j) ? 1 : 2;

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

  2147.                                     final double ijs = ghijCoef.getIjs(s, jms);
  2148.                                     final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2149.                                     final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2150.                                     final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2151.                                     final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2153.                                     final double jjs = ghijCoef.getJjs(s, jms);
  2154.                                     final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2155.                                     final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2156.                                     final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2157.                                     final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

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

  2166.                                     final double coef0 = d0smj * jn;
  2167.                                     final double coef1 = coef0 * lns;
  2168.                                     final double coef2 = coef1 * kns;

  2170.                                     final double coef3 = coef2 * ijs;
  2171.                                     final double coef4 = coef2 * jjs;

  2173.                                     // Add the term to the coefficients
  2174.                                     cCoef[0][j] -= coef3;
  2175.                                     cCoef[1][j] -= coef3 * (n + 1);
  2176.                                     cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2177.                                     cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2178.                                     cCoef[4][j] -= coef2 * dIjsdAlpha;
  2179.                                     cCoef[5][j] -= coef2 * dIjsdBeta;
  2180.                                     cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2182.                                     sCoef[0][j] += coef4;
  2183.                                     sCoef[1][j] += coef4 * (n + 1);
  2184.                                     sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2185.                                     sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2186.                                     sCoef[4][j] += coef2 * dJjsdAlpha;
  2187.                                     sCoef[5][j] += coef2 * dJjsdBeta;
  2188.                                     sCoef[6][j] += coef0 * dlns * kns * jjs;
  2189.                                 }
  2190.                             }
  2191.                         }

  2193.                         //compute second double sum where s: j/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  2194.                         for (int s = j / 2; s <=  FastMath.min(minjm1on - 1, sMax); s++) {
  2195.                             // j - s
  2196.                             final int jms = j - s;
  2197.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2198.                             final int d0smj = (s == j) ? 1 : 2;

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

  2207.                                     final double ijs = ghijCoef.getIjs(s, jms);
  2208.                                     final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2209.                                     final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2210.                                     final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2211.                                     final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2213.                                     final double jjs = ghijCoef.getJjs(s, jms);
  2214.                                     final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2215.                                     final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2216.                                     final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2217.                                     final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

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

  2226.                                     final double coef0 = d0smj * jn;
  2227.                                     final double coef1 = coef0 * lns;
  2228.                                     final double coef2 = coef1 * kns;

  2230.                                     final double coef3 = coef2 * ijs;
  2231.                                     final double coef4 = coef2 * jjs;

  2233.                                     // Add the term to the coefficients
  2234.                                     cCoef[0][j] -= coef3;
  2235.                                     cCoef[1][j] -= coef3 * (n + 1);
  2236.                                     cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2237.                                     cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2238.                                     cCoef[4][j] -= coef2 * dIjsdAlpha;
  2239.                                     cCoef[5][j] -= coef2 * dIjsdBeta;
  2240.                                     cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2242.                                     sCoef[0][j] += coef4;
  2243.                                     sCoef[1][j] += coef4 * (n + 1);
  2244.                                     sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2245.                                     sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2246.                                     sCoef[4][j] += coef2 * dJjsdAlpha;
  2247.                                     sCoef[5][j] += coef2 * dJjsdBeta;
  2248.                                     sCoef[6][j] += coef0 * dlns * kns * jjs;
  2249.                                 }
  2250.                             }
  2251.                         }
  2252.                     }

  2254.                     //if j is odd
  2255.                     else {
  2256.                         //compute first double sum where s: (j-1)/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  2257.                         for (int s = (j - 1) / 2; s <= FastMath.min(minjm1on - 1, sMax); s++) {
  2258.                             // j - s
  2259.                             final int jms = j - s;
  2260.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2261.                             final int d0smj = (s == j) ? 1 : 2;

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

  2270.                                     final double ijs = ghijCoef.getIjs(s, jms);
  2271.                                     final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2272.                                     final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2273.                                     final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2274.                                     final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2276.                                     final double jjs = ghijCoef.getJjs(s, jms);
  2277.                                     final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2278.                                     final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2279.                                     final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2280.                                     final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

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

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

  2290.                                     final double coef0 = d0smj * jn;
  2291.                                     final double coef1 = coef0 * lns;
  2292.                                     final double coef2 = coef1 * kns;

  2294.                                     final double coef3 = coef2 * ijs;
  2295.                                     final double coef4 = coef2 * jjs;

  2297.                                     // Add the term to the coefficients
  2298.                                     cCoef[0][j] -= coef3;
  2299.                                     cCoef[1][j] -= coef3 * (n + 1);
  2300.                                     cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2301.                                     cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2302.                                     cCoef[4][j] -= coef2 * dIjsdAlpha;
  2303.                                     cCoef[5][j] -= coef2 * dIjsdBeta;
  2304.                                     cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2306.                                     sCoef[0][j] += coef4;
  2307.                                     sCoef[1][j] += coef4 * (n + 1);
  2308.                                     sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2309.                                     sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2310.                                     sCoef[4][j] += coef2 * dJjsdAlpha;
  2311.                                     sCoef[5][j] += coef2 * dJjsdBeta;
  2312.                                     sCoef[6][j] += coef0 * dlns * kns * jjs;
  2313.                                 }
  2314.                             }
  2315.                         }

  2317.                         //the second double sum is added only if N >= 4 and j between 5 and 2*N-3
  2318.                         if (nMax >= 4 && isBetween(j, 5, 2 * nMax - 3)) {
  2319.                             //compute second double sum where s: j-min(j-1,N) -> (j-3)/2 and n: j-s -> min(j-1,N)
  2320.                             for (int s = j - minjm1on; s <= FastMath.min((j - 3) / 2, sMax); s++) {
  2321.                                 // j - s
  2322.                                 final int jms = j - s;
  2323.                                 // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2324.                                 final int d0smj = (s == j) ? 1 : 2;

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

  2333.                                         final double ijs = ghijCoef.getIjs(s, jms);
  2334.                                         final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2335.                                         final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2336.                                         final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2337.                                         final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2339.                                         final double jjs = ghijCoef.getJjs(s, jms);
  2340.                                         final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2341.                                         final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2342.                                         final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2343.                                         final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

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

  2352.                                         final double coef0 = d0smj * jn;
  2353.                                         final double coef1 = coef0 * lns;
  2354.                                         final double coef2 = coef1 * kns;

  2356.                                         final double coef3 = coef2 * ijs;
  2357.                                         final double coef4 = coef2 * jjs;

  2359.                                         // Add the term to the coefficients
  2360.                                         cCoef[0][j] -= coef3;
  2361.                                         cCoef[1][j] -= coef3 * (n + 1);
  2362.                                         cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2363.                                         cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2364.                                         cCoef[4][j] -= coef2 * dIjsdAlpha;
  2365.                                         cCoef[5][j] -= coef2 * dIjsdBeta;
  2366.                                         cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2368.                                         sCoef[0][j] += coef4;
  2369.                                         sCoef[1][j] += coef4 * (n + 1);
  2370.                                         sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2371.                                         sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2372.                                         sCoef[4][j] += coef2 * dJjsdAlpha;
  2373.                                         sCoef[5][j] += coef2 * dJjsdBeta;
  2374.                                         sCoef[6][j] += coef0 * dlns * kns * jjs;
  2375.                                     }
  2376.                                 }
  2377.                             }
  2378.                         }
  2379.                     }
  2380.                 }

  2382.                 cCoef[0][j] *= -context.getMuoa() / j;
  2383.                 cCoef[1][j] *=  context.getMuoa() / ( j * auxiliaryElements.getSma() );
  2384.                 cCoef[2][j] *= -context.getMuoa() / j;
  2385.                 cCoef[3][j] *= -context.getMuoa() / j;
  2386.                 cCoef[4][j] *= -context.getMuoa() / j;
  2387.                 cCoef[5][j] *= -context.getMuoa() / j;
  2388.                 cCoef[6][j] *= -context.getMuoa() / j;

  2390.                 sCoef[0][j] *= -context.getMuoa() / j;
  2391.                 sCoef[1][j] *=  context.getMuoa() / ( j * auxiliaryElements.getSma() );
  2392.                 sCoef[2][j] *= -context.getMuoa() / j;
  2393.                 sCoef[3][j] *= -context.getMuoa() / j;
  2394.                 sCoef[4][j] *= -context.getMuoa() / j;
  2395.                 sCoef[5][j] *= -context.getMuoa() / j;
  2396.                 sCoef[6][j] *= -context.getMuoa() / j;

  2398.             }
  2399.         }

  2401.         /** Check if an index is within the accepted interval.
  2402.          *
  2403.          * @param index the index to check
  2404.          * @param lowerBound the lower bound of the interval
  2405.          * @param upperBound the upper bound of the interval
  2406.          * @return true if the index is between the lower and upper bounds, false otherwise
  2407.          */
  2408.         private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
  2409.             return index >= lowerBound && index <= upperBound;
  2410.         }

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

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

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

  2439.         /**Get the value of dC<sup>j</sup> / dk.
  2440.          *
  2441.          * @param j j index
  2442.          * @return dC<sup>j</sup> / dk
  2443.          */
  2444.         public double getdCjdK(final int j) {
  2445.             return cCoef[3][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 getdCjdAlpha(final int j) {
  2454.             return cCoef[4][j];
  2455.         }

  2457.         /**Get the value of dC<sup>j</sup> / dβ.
  2458.          *
  2459.          * @param j j index
  2460.          * @return dC<sup>j</sup> / dβ
  2461.          */
  2462.         public double getdCjdBeta(final int j) {
  2463.             return cCoef[5][j];
  2464.         }

  2466.         /**Get the value of dC<sup>j</sup> / dγ.
  2467.          *
  2468.          * @param j j index
  2469.          * @return dC<sup>j</sup> / dγ
  2470.          */
  2471.         public double getdCjdGamma(final int j) {
  2472.             return cCoef[6][j];
  2473.         }

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

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

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

  2502.         /**Get the value of dS<sup>j</sup> / dk.
  2503.          *
  2504.          * @param j j index
  2505.          * @return dS<sup>j</sup> / dk
  2506.          */
  2507.         public double getdSjdK(final int j) {
  2508.             return sCoef[3][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 getdSjdAlpha(final int j) {
  2517.             return sCoef[4][j];
  2518.         }

  2520.         /**Get the value of dS<sup>j</sup> / dβ.
  2521.          *
  2522.          * @param j j index
  2523.          * @return dS<sup>j</sup> / dβ
  2524.          */
  2525.         public double getdSjdBeta(final int j) {
  2526.             return sCoef[5][j];
  2527.         }

  2529.         /**Get the value of dS<sup>j</sup> /  dγ.
  2530.          *
  2531.          * @param j j index
  2532.          * @return dS<sup>j</sup> /  dγ
  2533.          */
  2534.         public double getdSjdGamma(final int j) {
  2535.             return sCoef[6][j];
  2536.         }
  2537.     }

  2539.     /** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
  2540.      *  <p>
  2541.      *  Those coefficients are given in Danielson paper by expressions 4.1-(13) to 4.1.-(16b)
  2542.      *  </p>
  2543.      */
  2544.     private class FieldFourierCjSjCoefficients <T extends RealFieldElement<T>> {

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

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

  2552.         /** Maximum possible value for n. */
  2553.         private final int nMax;

  2555.         /** Maximum possible value for s. */
  2556.         private final int sMax;

  2558.         /** Maximum possible value for j. */
  2559.         private final int jMax;

  2561.         /** The C<sup>j</sup> coefficients and their derivatives.
  2562.          * <p>
  2563.          * Each column of the matrix contains the following values: <br/>
  2564.          * - C<sup>j</sup> <br/>
  2565.          * - dC<sup>j</sup> / da <br/>
  2566.          * - dC<sup>j</sup> / dh <br/>
  2567.          * - dC<sup>j</sup> / dk <br/>
  2568.          * - dC<sup>j</sup> / dα <br/>
  2569.          * - dC<sup>j</sup> / dβ <br/>
  2570.          * - dC<sup>j</sup> / dγ <br/>
  2571.          * </p>
  2572.          */
  2573.         private final T[][] cCoef;

  2575.         /** The S<sup>j</sup> coefficients and their derivatives.
  2576.          * <p>
  2577.          * Each column of the matrix contains the following values: <br/>
  2578.          * - S<sup>j</sup> <br/>
  2579.          * - dS<sup>j</sup> / da <br/>
  2580.          * - dS<sup>j</sup> / dh <br/>
  2581.          * - dS<sup>j</sup> / dk <br/>
  2582.          * - dS<sup>j</sup> / dα <br/>
  2583.          * - dS<sup>j</sup> / dβ <br/>
  2584.          * - dS<sup>j</sup> / dγ <br/>
  2585.          * </p>
  2586.          */
  2587.         private final T[][] sCoef;

  2589.         /** h * &Chi;³. */
  2590.         private final T hXXX;
  2591.         /** k * &Chi;³. */
  2592.         private final T kXXX;

  2594.         /** Create a set of C<sup>j</sup> and the S<sup>j</sup> coefficients.
  2595.          *  @param date the current date
  2596.          *  @param nMax maximum possible value for n
  2597.          *  @param sMax maximum possible value for s
  2598.          *  @param jMax maximum possible value for j
  2599.          *  @param context container for attributes
  2600.          *  @param hansenObjects initialization of hansen objects
  2601.          */
  2602.         FieldFourierCjSjCoefficients(final FieldAbsoluteDate<T> date,
  2603.                                      final int nMax, final int sMax, final int jMax,
  2604.                                      final FieldDSSTZonalContext<T> context,
  2605.                                      final FieldHansenObjects<T> hansenObjects) {

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

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

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

  2616.             this.lnsCoef = new FieldLnsCoefficients<>(nMax, nMax, Qns, Vns, context.getRoa(), field);
  2617.             this.nMax = nMax;
  2618.             this.sMax = sMax;
  2619.             this.jMax = jMax;

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

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

  2628.             for (int s = 0; s <= sMax; s++) {
  2629.                 //Initialise the Hansen roots
  2630.                 hansenObjects.computeHansenObjectsInitValues(context, s);
  2631.             }
  2632.             generateCoefficients(date, context, auxiliaryElements, hansenObjects, field);
  2633.         }

  2635.         /** Generate all coefficients.
  2636.          * @param date the current date
  2637.          * @param context container for attributes
  2638.          * @param hansenObjects initialization of hansen objects
  2639.          * @param auxiliaryElements auxiliary elements related to the current orbit
  2640.          * @param field field used by default
  2641.          */
  2642.         private void generateCoefficients(final FieldAbsoluteDate<T> date,
  2643.                                           final FieldDSSTZonalContext<T> context,
  2644.                                           final FieldAuxiliaryElements<T> auxiliaryElements,
  2645.                                           final FieldHansenObjects<T> hansenObjects,
  2646.                                           final Field<T> field) {

  2648.             //Zero
  2649.             final T zero = field.getZero();

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

  2654.                 //init arrays
  2655.                 for (int i = 0; i <= 6; i++) {
  2656.                     cCoef[i][j] = zero;
  2657.                     sCoef[i][j] = zero;
  2658.                 }

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

  2662.                     //compute first double sum where s: j -> N-1 and n: s+1 -> N
  2663.                     for (int s = j; s <= FastMath.min(nMax - 1, sMax); s++) {
  2664.                         // j - s
  2665.                         final int jms = j - s;
  2666.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2667.                         final int d0smj = (s == j) ? 1 : 2;

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

  2676.                                 final T hjs        = ghijCoef.getHjs(s, -jms);
  2677.                                 final T dHjsdh     = ghijCoef.getdHjsdh(s, -jms);
  2678.                                 final T dHjsdk     = ghijCoef.getdHjsdk(s, -jms);
  2679.                                 final T dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, -jms);
  2680.                                 final T dHjsdBeta  = ghijCoef.getdHjsdBeta(s, -jms);

  2682.                                 final T gjs        = ghijCoef.getGjs(s, -jms);
  2683.                                 final T dGjsdh     = ghijCoef.getdGjsdh(s, -jms);
  2684.                                 final T dGjsdk     = ghijCoef.getdGjsdk(s, -jms);
  2685.                                 final T dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, -jms);
  2686.                                 final T dGjsdBeta  = ghijCoef.getdGjsdBeta(s, -jms);

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

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

  2695.                                 final T coef0 = jn.multiply(d0smj);
  2696.                                 final T coef1 = coef0.multiply(lns);
  2697.                                 final T coef2 = coef1.multiply(kns);
  2698.                                 final T coef3 = coef2.multiply(hjs);
  2699.                                 final T coef4 = coef2.multiply(gjs);

  2701.                                 // Add the term to the coefficients
  2702.                                 cCoef[0][j] = cCoef[0][j].add(coef3);
  2703.                                 cCoef[1][j] = cCoef[1][j].add(coef3.multiply(n + 1));
  2704.                                 cCoef[2][j] = cCoef[2][j].add(coef1.multiply(kns.multiply(dHjsdh).add(hjs.multiply(hXXX).multiply(dkns))));
  2705.                                 cCoef[3][j] = cCoef[3][j].add(coef1.multiply(kns.multiply(dHjsdk).add(hjs.multiply(kXXX).multiply(dkns))));
  2706.                                 cCoef[4][j] = cCoef[4][j].add(coef2.multiply(dHjsdAlpha));
  2707.                                 cCoef[5][j] = cCoef[5][j].add(coef2.multiply(dHjsdBeta));
  2708.                                 cCoef[6][j] = cCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(hjs));

  2710.                                 sCoef[0][j] = sCoef[0][j].add(coef4);
  2711.                                 sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  2712.                                 sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dGjsdh).add(gjs.multiply(hXXX).multiply(dkns))));
  2713.                                 sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dGjsdk).add(gjs.multiply(kXXX).multiply(dkns))));
  2714.                                 sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dGjsdAlpha));
  2715.                                 sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dGjsdBeta));
  2716.                                 sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(gjs));
  2717.                             }
  2718.                         }
  2719.                     }

  2721.                     //compute second double sum where s: 0 -> N-j and n: max(j+s, j+1) -> N
  2722.                     for (int s = 0; s <= FastMath.min(nMax - j, sMax); s++) {
  2723.                         // j + s
  2724.                         final int jps = j + s;
  2725.                         // Kronecker symbols (2 - delta(0,j+s))
  2726.                         final double d0spj = (s == -j) ? 1 : 2;

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

  2735.                                 final T hjs        = ghijCoef.getHjs(s, jps);
  2736.                                 final T dHjsdh     = ghijCoef.getdHjsdh(s, jps);
  2737.                                 final T dHjsdk     = ghijCoef.getdHjsdk(s, jps);
  2738.                                 final T dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, jps);
  2739.                                 final T dHjsdBeta  = ghijCoef.getdHjsdBeta(s, jps);

  2741.                                 final T gjs        = ghijCoef.getGjs(s, jps);
  2742.                                 final T dGjsdh     = ghijCoef.getdGjsdh(s, jps);
  2743.                                 final T dGjsdk     = ghijCoef.getdGjsdk(s, jps);
  2744.                                 final T dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, jps);
  2745.                                 final T dGjsdBeta  = ghijCoef.getdGjsdBeta(s, jps);

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

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

  2754.                                 final T coef0 = jn.multiply(d0spj);
  2755.                                 final T coef1 = coef0.multiply(lns);
  2756.                                 final T coef2 = coef1.multiply(kns);

  2758.                                 final T coef3 = coef2.multiply(hjs);
  2759.                                 final T coef4 = coef2.multiply(gjs);

  2761.                                 // Add the term to the coefficients
  2762.                                 cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2763.                                 cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  2764.                                 cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dHjsdh).add(hjs.multiply(hXXX).multiply(dkns))));
  2765.                                 cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dHjsdk).add(hjs.multiply(kXXX).multiply(dkns))));
  2766.                                 cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dHjsdAlpha));
  2767.                                 cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dHjsdBeta));
  2768.                                 cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(hjs));

  2770.                                 sCoef[0][j] = sCoef[0][j].add(coef4);
  2771.                                 sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  2772.                                 sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dGjsdh).add(gjs.multiply(hXXX).multiply(dkns))));
  2773.                                 sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dGjsdk).add(gjs.multiply(kXXX).multiply(dkns))));
  2774.                                 sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dGjsdAlpha));
  2775.                                 sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dGjsdBeta));
  2776.                                 sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(gjs));
  2777.                             }
  2778.                         }
  2779.                     }

  2781.                     //compute third double sum where s: 1 -> j and  n: j+1 -> N
  2782.                     for (int s = 1; s <= FastMath.min(j, sMax); s++) {
  2783.                         // j - s
  2784.                         final int jms = j - s;
  2785.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2786.                         final int d0smj = (s == j) ? 1 : 2;

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

  2795.                                 final T ijs        = ghijCoef.getIjs(s, jms);
  2796.                                 final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  2797.                                 final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  2798.                                 final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2799.                                 final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  2801.                                 final T jjs        = ghijCoef.getJjs(s, jms);
  2802.                                 final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  2803.                                 final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  2804.                                 final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2805.                                 final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

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

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

  2814.                                 final T coef0 = jn.multiply(d0smj);
  2815.                                 final T coef1 = coef0.multiply(lns);
  2816.                                 final T coef2 = coef1.multiply(kns);

  2818.                                 final T coef3 = coef2.multiply(ijs);
  2819.                                 final T coef4 = coef2.multiply(jjs);

  2821.                                 // Add the term to the coefficients
  2822.                                 cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2823.                                 cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  2824.                                 cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  2825.                                 cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  2826.                                 cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  2827.                                 cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  2828.                                 cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  2830.                                 sCoef[0][j] = sCoef[0][j].add(coef4);
  2831.                                 sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  2832.                                 sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  2833.                                 sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  2834.                                 sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  2835.                                 sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  2836.                                 sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  2837.                             }
  2838.                         }
  2839.                     }
  2840.                 }

  2842.                 if (isBetween(j, 2, nMax)) {
  2843.                     //add first term
  2844.                     // J<sub>j</sub>
  2845.                     final T jj = zero.subtract(harmonics.getUnnormalizedCnm(j, 0));
  2846.                     T kns  = (T) hansenObjects.getHansenObjects()[0].getValue(-j - 1, context.getX());
  2847.                     T dkns = (T) hansenObjects.getHansenObjects()[0].getDerivative(-j - 1, context.getX());

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

  2852.                     final T hjs        = ghijCoef.getHjs(0, j);
  2853.                     final T dHjsdh     = ghijCoef.getdHjsdh(0, j);
  2854.                     final T dHjsdk     = ghijCoef.getdHjsdk(0, j);
  2855.                     final T dHjsdAlpha = ghijCoef.getdHjsdAlpha(0, j);
  2856.                     final T dHjsdBeta  = ghijCoef.getdHjsdBeta(0, j);

  2858.                     final T gjs        = ghijCoef.getGjs(0, j);
  2859.                     final T dGjsdh     = ghijCoef.getdGjsdh(0, j);
  2860.                     final T dGjsdk     = ghijCoef.getdGjsdk(0, j);
  2861.                     final T dGjsdAlpha = ghijCoef.getdGjsdAlpha(0, j);
  2862.                     final T dGjsdBeta  = ghijCoef.getdGjsdBeta(0, j);

  2864.                     // 2 * J<sub>j</sub> * K₀<sup>-j-1,0</sup> * L<sub>j</sub><sup>j</sup>
  2865.                     T coef0 = jj.multiply(2.);
  2866.                     T coef1 = coef0.multiply(lns);
  2867.                     T coef2 = coef1.multiply(kns);

  2869.                     T coef3 = coef2.multiply(hjs);
  2870.                     T coef4 = coef2.multiply(gjs);

  2872.                     // Add the term to the coefficients
  2873.                     cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2874.                     cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(j + 1));
  2875.                     cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dHjsdh).add(hjs.multiply(hXXX).multiply(dkns))));
  2876.                     cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dHjsdk).add(hjs.multiply(kXXX).multiply(dkns))));
  2877.                     cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dHjsdAlpha));
  2878.                     cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dHjsdBeta));
  2879.                     //no contribution to cCoef[6][j] because dlns is 0

  2881.                     sCoef[0][j] = sCoef[0][j].add(coef4);
  2882.                     sCoef[1][j] = sCoef[1][j].add(coef4.multiply(j + 1));
  2883.                     sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dGjsdh).add(gjs.multiply(hXXX).multiply(dkns))));
  2884.                     sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dGjsdk).add(gjs.multiply(kXXX).multiply(dkns))));
  2885.                     sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dGjsdAlpha));
  2886.                     sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dGjsdBeta));
  2887.                     //no contribution to sCoef[6][j] because dlns is 0

  2889.                     //compute simple sum where s: 1 -> j-1
  2890.                     for (int s = 1; s <= FastMath.min(j - 1, sMax); s++) {
  2891.                         // j - s
  2892.                         final int jms = j - s;
  2893.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2894.                         final int d0smj = (s == j) ? 1 : 2;

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

  2902.                             lns = lnsCoef.getLns(j, jms);
  2903.                             final T dlns = lnsCoef.getdLnsdGamma(j, jms);

  2905.                             final T ijs        = ghijCoef.getIjs(s, jms);
  2906.                             final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  2907.                             final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  2908.                             final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2909.                             final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  2911.                             final T jjs        = ghijCoef.getJjs(s, jms);
  2912.                             final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  2913.                             final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  2914.                             final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2915.                             final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

  2917.                             coef0 = jj.multiply(d0smj);
  2918.                             coef1 = coef0.multiply(lns);
  2919.                             coef2 = coef1.multiply(kns);

  2921.                             coef3 = coef2.multiply(ijs);
  2922.                             coef4 = coef2.multiply(jjs);

  2924.                             // Add the term to the coefficients
  2925.                             cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2926.                             cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(j + 1));
  2927.                             cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  2928.                             cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  2929.                             cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  2930.                             cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  2931.                             cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  2933.                             sCoef[0][j] = sCoef[0][j].add(coef4);
  2934.                             sCoef[1][j] = sCoef[1][j].add(coef4.multiply(j + 1));
  2935.                             sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  2936.                             sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  2937.                             sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  2938.                             sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  2939.                             sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  2940.                         }
  2941.                     }
  2942.                 }

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

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

  2950.                     //if j is even
  2951.                     if (j % 2 == 0) {
  2952.                         //compute first double sum where s: j-min(j-1,N) -> j/2-1 and n: j-s -> min(j-1,N)
  2953.                         for (int s = j - minjm1on; s <= FastMath.min(j / 2 - 1, sMax); s++) {
  2954.                             // j - s
  2955.                             final int jms = j - s;
  2956.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2957.                             final int d0smj = (s == j) ? 1 : 2;

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

  2966.                                     final T ijs        = ghijCoef.getIjs(s, jms);
  2967.                                     final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  2968.                                     final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  2969.                                     final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2970.                                     final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  2972.                                     final T jjs        = ghijCoef.getJjs(s, jms);
  2973.                                     final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  2974.                                     final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  2975.                                     final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2976.                                     final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

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

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

  2985.                                     final T coef0 = jn.multiply(d0smj);
  2986.                                     final T coef1 = coef0.multiply(lns);
  2987.                                     final T coef2 = coef1.multiply(kns);

  2989.                                     final T coef3 = coef2.multiply(ijs);
  2990.                                     final T coef4 = coef2.multiply(jjs);

  2992.                                     // Add the term to the coefficients
  2993.                                     cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2994.                                     cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  2995.                                     cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  2996.                                     cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  2997.                                     cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  2998.                                     cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  2999.                                     cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  3001.                                     sCoef[0][j] = sCoef[0][j].add(coef4);
  3002.                                     sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  3003.                                     sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  3004.                                     sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  3005.                                     sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  3006.                                     sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  3007.                                     sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  3008.                                 }
  3009.                             }
  3010.                         }

  3012.                         //compute second double sum where s: j/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  3013.                         for (int s = j / 2; s <=  FastMath.min(minjm1on - 1, sMax); s++) {
  3014.                             // j - s
  3015.                             final int jms = j - s;
  3016.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  3017.                             final int d0smj = (s == j) ? 1 : 2;

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

  3026.                                     final T ijs        = ghijCoef.getIjs(s, jms);
  3027.                                     final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  3028.                                     final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  3029.                                     final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  3030.                                     final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  3032.                                     final T jjs        = ghijCoef.getJjs(s, jms);
  3033.                                     final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  3034.                                     final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  3035.                                     final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  3036.                                     final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

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

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

  3045.                                     final T coef0 = jn.multiply(d0smj);
  3046.                                     final T coef1 = coef0.multiply(lns);
  3047.                                     final T coef2 = coef1.multiply(kns);

  3049.                                     final T coef3 = coef2.multiply(ijs);
  3050.                                     final T coef4 = coef2.multiply(jjs);

  3052.                                     // Add the term to the coefficients
  3053.                                     cCoef[0][j] = cCoef[0][j].subtract(coef3);
  3054.                                     cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  3055.                                     cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  3056.                                     cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  3057.                                     cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  3058.                                     cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  3059.                                     cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  3061.                                     sCoef[0][j] = sCoef[0][j].add(coef4);
  3062.                                     sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  3063.                                     sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  3064.                                     sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  3065.                                     sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  3066.                                     sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  3067.                                     sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  3068.                                 }
  3069.                             }
  3070.                         }
  3071.                     }

  3073.                     //if j is odd
  3074.                     else {
  3075.                         //compute first double sum where s: (j-1)/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  3076.                         for (int s = (j - 1) / 2; s <= FastMath.min(minjm1on - 1, sMax); s++) {
  3077.                             // j - s
  3078.                             final int jms = j - s;
  3079.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  3080.                             final int d0smj = (s == j) ? 1 : 2;

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

  3089.                                     final T ijs        = ghijCoef.getIjs(s, jms);
  3090.                                     final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  3091.                                     final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  3092.                                     final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  3093.                                     final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  3095.                                     final T jjs        = ghijCoef.getJjs(s, jms);
  3096.                                     final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  3097.                                     final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  3098.                                     final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  3099.                                     final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

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

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

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

  3109.                                     final T coef0 = jn.multiply(d0smj);
  3110.                                     final T coef1 = coef0.multiply(lns);
  3111.                                     final T coef2 = coef1.multiply(kns);

  3113.                                     final T coef3 = coef2.multiply(ijs);
  3114.                                     final T coef4 = coef2.multiply(jjs);

  3116.                                     // Add the term to the coefficients
  3117.                                     cCoef[0][j] = cCoef[0][j].subtract(coef3);
  3118.                                     cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  3119.                                     cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  3120.                                     cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  3121.                                     cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  3122.                                     cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  3123.                                     cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  3125.                                     sCoef[0][j] = sCoef[0][j].add(coef4);
  3126.                                     sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  3127.                                     sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  3128.                                     sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  3129.                                     sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  3130.                                     sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  3131.                                     sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  3132.                                 }
  3133.                             }
  3134.                         }

  3136.                         //the second double sum is added only if N >= 4 and j between 5 and 2*N-3
  3137.                         if (nMax >= 4 && isBetween(j, 5, 2 * nMax - 3)) {
  3138.                             //compute second double sum where s: j-min(j-1,N) -> (j-3)/2 and n: j-s -> min(j-1,N)
  3139.                             for (int s = j - minjm1on; s <= FastMath.min((j - 3) / 2, sMax); s++) {
  3140.                                 // j - s
  3141.                                 final int jms = j - s;
  3142.                                 // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  3143.                                 final int d0smj = (s == j) ? 1 : 2;

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

  3152.                                         final T ijs        = ghijCoef.getIjs(s, jms);
  3153.                                         final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  3154.                                         final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  3155.                                         final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  3156.                                         final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  3158.                                         final T jjs        = ghijCoef.getJjs(s, jms);
  3159.                                         final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  3160.                                         final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  3161.                                         final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  3162.                                         final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

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

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

  3171.                                         final T coef0 = jn.multiply(d0smj);
  3172.                                         final T coef1 = coef0.multiply(lns);
  3173.                                         final T coef2 = coef1.multiply(kns);

  3175.                                         final T coef3 = coef2.multiply(ijs);
  3176.                                         final T coef4 = coef2.multiply(jjs);

  3178.                                         // Add the term to the coefficients
  3179.                                         cCoef[0][j] = cCoef[0][j].subtract(coef3);
  3180.                                         cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  3181.                                         cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  3182.                                         cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  3183.                                         cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  3184.                                         cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  3185.                                         cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  3187.                                         sCoef[0][j] = sCoef[0][j].add(coef4);
  3188.                                         sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  3189.                                         sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  3190.                                         sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  3191.                                         sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  3192.                                         sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  3193.                                         sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  3194.                                     }
  3195.                                 }
  3196.                             }
  3197.                         }
  3198.                     }
  3199.                 }

  3201.                 cCoef[0][j] = cCoef[0][j].multiply(context.getMuoa().divide(j).negate());
  3202.                 cCoef[1][j] = cCoef[1][j].multiply(context.getMuoa().divide(auxiliaryElements.getSma().multiply(j)));
  3203.                 cCoef[2][j] = cCoef[2][j].multiply(context.getMuoa().divide(j).negate());
  3204.                 cCoef[3][j] = cCoef[3][j].multiply(context.getMuoa().divide(j).negate());
  3205.                 cCoef[4][j] = cCoef[4][j].multiply(context.getMuoa().divide(j).negate());
  3206.                 cCoef[5][j] = cCoef[5][j].multiply(context.getMuoa().divide(j).negate());
  3207.                 cCoef[6][j] = cCoef[6][j].multiply(context.getMuoa().divide(j).negate());

  3209.                 sCoef[0][j] = sCoef[0][j].multiply(context.getMuoa().divide(j).negate());
  3210.                 sCoef[1][j] = sCoef[1][j].multiply(context.getMuoa().divide(auxiliaryElements.getSma().multiply(j)));
  3211.                 sCoef[2][j] = sCoef[2][j].multiply(context.getMuoa().divide(j).negate());
  3212.                 sCoef[3][j] = sCoef[3][j].multiply(context.getMuoa().divide(j).negate());
  3213.                 sCoef[4][j] = sCoef[4][j].multiply(context.getMuoa().divide(j).negate());
  3214.                 sCoef[5][j] = sCoef[5][j].multiply(context.getMuoa().divide(j).negate());
  3215.                 sCoef[6][j] = sCoef[6][j].multiply(context.getMuoa().divide(j).negate());

  3217.             }
  3218.         }

  3220.         /** Check if an index is within the accepted interval.
  3221.          *
  3222.          * @param index the index to check
  3223.          * @param lowerBound the lower bound of the interval
  3224.          * @param upperBound the upper bound of the interval
  3225.          * @return true if the index is between the lower and upper bounds, false otherwise
  3226.          */
  3227.         private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
  3228.             return index >= lowerBound && index <= upperBound;
  3229.         }

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

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

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

  3258.         /**Get the value of dC<sup>j</sup> / dk.
  3259.          *
  3260.          * @param j j index
  3261.          * @return dC<sup>j</sup> / dk
  3262.          */
  3263.         public T getdCjdK(final int j) {
  3264.             return cCoef[3][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 getdCjdAlpha(final int j) {
  3273.             return cCoef[4][j];
  3274.         }

  3276.         /**Get the value of dC<sup>j</sup> / dβ.
  3277.          *
  3278.          * @param j j index
  3279.          * @return dC<sup>j</sup> / dβ
  3280.          */
  3281.         public T getdCjdBeta(final int j) {
  3282.             return cCoef[5][j];
  3283.         }

  3285.         /**Get the value of dC<sup>j</sup> / dγ.
  3286.          *
  3287.          * @param j j index
  3288.          * @return dC<sup>j</sup> / dγ
  3289.          */
  3290.         public T getdCjdGamma(final int j) {
  3291.             return cCoef[6][j];
  3292.         }

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

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

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

  3321.         /**Get the value of dS<sup>j</sup> / dk.
  3322.          *
  3323.          * @param j j index
  3324.          * @return dS<sup>j</sup> / dk
  3325.          */
  3326.         public T getdSjdK(final int j) {
  3327.             return sCoef[3][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 getdSjdAlpha(final int j) {
  3336.             return sCoef[4][j];
  3337.         }

  3339.         /**Get the value of dS<sup>j</sup> / dβ.
  3340.          *
  3341.          * @param j j index
  3342.          * @return dS<sup>j</sup> / dβ
  3343.          */
  3344.         public T getdSjdBeta(final int j) {
  3345.             return sCoef[5][j];
  3346.         }

  3348.         /**Get the value of dS<sup>j</sup> /  dγ.
  3349.          *
  3350.          * @param j j index
  3351.          * @return dS<sup>j</sup> /  dγ
  3352.          */
  3353.         public T getdSjdGamma(final int j) {
  3354.             return sCoef[6][j];
  3355.         }
  3356.     }

  3358.     /** Coefficients valid for one time slot. */
  3359.     private static class Slot {

  3361.         /**The coefficients D<sub>i</sub>.
  3362.          * <p>
  3363.          * i corresponds to the equinoctial element, as follows:
  3364.          * - i=0 for a <br/>
  3365.          * - i=1 for k <br/>
  3366.          * - i=2 for h <br/>
  3367.          * - i=3 for q <br/>
  3368.          * - i=4 for p <br/>
  3369.          * - i=5 for λ <br/>
  3370.          * </p>
  3371.          */
  3372.         private final ShortPeriodicsInterpolatedCoefficient di;

  3374.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  3375.          * <p>
  3376.          * The constant term C<sub>i</sub>⁰ is also stored in this variable at index j = 0 <br>
  3377.          * The index order is cij[j][i] <br/>
  3378.          * i corresponds to the equinoctial element, as follows: <br/>
  3379.          * - i=0 for a <br/>
  3380.          * - i=1 for k <br/>
  3381.          * - i=2 for h <br/>
  3382.          * - i=3 for q <br/>
  3383.          * - i=4 for p <br/>
  3384.          * - i=5 for λ <br/>
  3385.          * </p>
  3386.          */
  3387.         private final ShortPeriodicsInterpolatedCoefficient[] cij;

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

  3403.         /** Simple constructor.
  3404.          *  @param maxFrequencyShortPeriodics maximum value for j index
  3405.          *  @param interpolationPoints number of points used in the interpolation process
  3406.          */
  3407.         Slot(final int maxFrequencyShortPeriodics, final int interpolationPoints) {

  3409.             final int rows = maxFrequencyShortPeriodics + 1;
  3410.             di  = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3411.             cij = new ShortPeriodicsInterpolatedCoefficient[rows];
  3412.             sij = new ShortPeriodicsInterpolatedCoefficient[rows];

  3414.             //Initialize the arrays
  3415.             for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
  3416.                 cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3417.                 sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3418.             }

  3420.         }

  3422.     }

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

  3427.         /**The coefficients D<sub>i</sub>.
  3428.          * <p>
  3429.          * i corresponds to the equinoctial element, as follows:
  3430.          * - i=0 for a <br/>
  3431.          * - i=1 for k <br/>
  3432.          * - i=2 for h <br/>
  3433.          * - i=3 for q <br/>
  3434.          * - i=4 for p <br/>
  3435.          * - i=5 for λ <br/>
  3436.          * </p>
  3437.          */
  3438.         private final FieldShortPeriodicsInterpolatedCoefficient<T> di;

  3440.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  3441.          * <p>
  3442.          * The constant term C<sub>i</sub>⁰ is also stored in this variable at index j = 0 <br>
  3443.          * The index order is cij[j][i] <br/>
  3444.          * i corresponds to the equinoctial element, as follows: <br/>
  3445.          * - i=0 for a <br/>
  3446.          * - i=1 for k <br/>
  3447.          * - i=2 for h <br/>
  3448.          * - i=3 for q <br/>
  3449.          * - i=4 for p <br/>
  3450.          * - i=5 for λ <br/>
  3451.          * </p>
  3452.          */
  3453.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] cij;

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

  3469.         /** Simple constructor.
  3470.          *  @param maxFrequencyShortPeriodics maximum value for j index
  3471.          *  @param interpolationPoints number of points used in the interpolation process
  3472.          */
  3473.         @SuppressWarnings("unchecked")
  3474.         FieldSlot(final int maxFrequencyShortPeriodics, final int interpolationPoints) {

  3476.             final int rows = maxFrequencyShortPeriodics + 1;
  3477.             di  = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3478.             cij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows);
  3479.             sij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows);

  3481.             //Initialize the arrays
  3482.             for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
  3483.                 cij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3484.                 sij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3485.             }

  3487.         }

  3489.     }

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

  3494.         /** The current value of the U function. <br/>
  3495.          * Needed for the short periodic contribution */
  3496.         private double U;

  3498.         /** dU / da. */
  3499.         private double dUda;

  3501.         /** dU / dk. */
  3502.         private double dUdk;

  3504.         /** dU / dh. */
  3505.         private double dUdh;

  3507.         /** dU / dAlpha. */
  3508.         private double dUdAl;

  3510.         /** dU / dBeta. */
  3511.         private double dUdBe;

  3513.         /** dU / dGamma. */
  3514.         private double dUdGa;

  3516.         /** Simple constuctor.
  3517.          *  @param date current date
  3518.          *  @param context container for attributes
  3519.          *  @param auxiliaryElements auxiliary elements related to the current orbit
  3520.          *  @param hansen initialization of hansen objects
  3521.          */
  3522.         UAnddU(final AbsoluteDate date,
  3523.                final DSSTZonalContext context,
  3524.                final AuxiliaryElements auxiliaryElements,
  3525.                final HansenObjects hansen) {

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

  3529.             //Reset U
  3530.             U = 0.;

  3532.             // Gs and Hs coefficients
  3533.             final double[][] GsHs = CoefficientsFactory.computeGsHs(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta(), maxEccPowMeanElements);
  3534.             // Qns coefficients
  3535.             final double[][] Qns  = CoefficientsFactory.computeQns(auxiliaryElements.getGamma(), maxDegree, maxEccPowMeanElements);

  3537.             final double[] roaPow = new double[maxDegree + 1];
  3538.             roaPow[0] = 1.;
  3539.             for (int i = 1; i <= maxDegree; i++) {
  3540.                 roaPow[i] = context.getRoa() * roaPow[i - 1];
  3541.             }

  3543.             // Potential derivatives
  3544.             dUda  = 0.;
  3545.             dUdk  = 0.;
  3546.             dUdh  = 0.;
  3547.             dUdAl = 0.;
  3548.             dUdBe = 0.;
  3549.             dUdGa = 0.;

  3551.             for (int s = 0; s <= maxEccPowMeanElements; s++) {
  3552.                 //Initialize the Hansen roots
  3553.                 hansen.computeHansenObjectsInitValues(context, s);

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

  3558.                 // Compute Gs partial derivatives from 3.1-(9)
  3559.                 double dGsdh  = 0.;
  3560.                 double dGsdk  = 0.;
  3561.                 double dGsdAl = 0.;
  3562.                 double dGsdBe = 0.;
  3563.                 if (s > 0) {
  3564.                     // First get the G(s-1) and the H(s-1) coefficients
  3565.                     final double sxgsm1 = s * GsHs[0][s - 1];
  3566.                     final double sxhsm1 = s * GsHs[1][s - 1];
  3567.                     // Then compute derivatives
  3568.                     dGsdh  = auxiliaryElements.getBeta()  * sxgsm1 - auxiliaryElements.getAlpha() * sxhsm1;
  3569.                     dGsdk  = auxiliaryElements.getAlpha() * sxgsm1 + auxiliaryElements.getBeta()  * sxhsm1;
  3570.                     dGsdAl = auxiliaryElements.getK() * sxgsm1 - auxiliaryElements.getH() * sxhsm1;
  3571.                     dGsdBe = auxiliaryElements.getH() * sxgsm1 + auxiliaryElements.getK() * sxhsm1;
  3572.                 }

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

  3577.                 for (int n = s + 2; n <= maxDegree; n++) {
  3578.                     // (n - s) must be even
  3579.                     if ((n - s) % 2 == 0) {

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

  3585.                         final double vns   = Vns.get(new NSKey(n, s));
  3586.                         final double coef0 = d0s * roaPow[n] * vns * -harmonics.getUnnormalizedCnm(n, 0);
  3587.                         final double coef1 = coef0 * Qns[n][s];
  3588.                         final double coef2 = coef1 * kns;
  3589.                         final double coef3 = coef2 * gs;
  3590.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  3591.                         final double dqns  = Qns[n][s + 1];

  3593.                         // Compute U
  3594.                         U += coef3;
  3595.                         // Compute dU / da :
  3596.                         dUda += coef3 * (n + 1);
  3597.                         // Compute dU / dEx
  3598.                         dUdk += coef1 * (kns * dGsdk + auxiliaryElements.getK() * context.getXXX() * gs * dkns);
  3599.                         // Compute dU / dEy
  3600.                         dUdh += coef1 * (kns * dGsdh + auxiliaryElements.getH() * context.getXXX() * gs * dkns);
  3601.                         // Compute dU / dAlpha
  3602.                         dUdAl += coef2 * dGsdAl;
  3603.                         // Compute dU / dBeta
  3604.                         dUdBe += coef2 * dGsdBe;
  3605.                         // Compute dU / dGamma
  3606.                         dUdGa += coef0 * kns * dqns * gs;

  3608.                     }
  3609.                 }
  3610.             }

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

  3615.             this.dUda = dUda *  context.getMuoa() / auxiliaryElements.getSma();
  3616.             this.dUdk = dUdk * -context.getMuoa();
  3617.             this.dUdh = dUdh * -context.getMuoa();
  3618.             this.dUdAl = dUdAl * -context.getMuoa();
  3619.             this.dUdBe = dUdBe * -context.getMuoa();
  3620.             this.dUdGa = dUdGa * -context.getMuoa();

  3622.         }

  3624.         /** Return value of U.
  3625.          * @return U
  3626.          */
  3627.         public double getU() {
  3628.             return U;
  3629.         }

  3631.         /** Return value of dU / da.
  3632.          * @return dUda
  3633.          */
  3634.         public double getdUda() {
  3635.             return dUda;
  3636.         }

  3638.         /** Return value of dU / dk.
  3639.          * @return dUdk
  3640.          */
  3641.         public double getdUdk() {
  3642.             return dUdk;
  3643.         }

  3645.         /** Return value of dU / dh.
  3646.          * @return dUdh
  3647.          */
  3648.         public double getdUdh() {
  3649.             return dUdh;
  3650.         }

  3652.         /** Return value of dU / dAlpha.
  3653.          * @return dUdAl
  3654.          */
  3655.         public double getdUdAl() {
  3656.             return dUdAl;
  3657.         }

  3659.         /** Return value of dU / dBeta.
  3660.          * @return dUdBe
  3661.          */
  3662.         public double getdUdBe() {
  3663.             return dUdBe;
  3664.         }

  3666.         /** Return value of dU / dGamma.
  3667.          * @return dUdGa
  3668.          */
  3669.         public double getdUdGa() {
  3670.             return dUdGa;
  3671.         }

  3673.     }

  3675.     /** Compute the derivatives of the gravitational potential U [Eq. 3.1-(6)].
  3676.      *  <p>
  3677.      *  The result is the array
  3678.      *  [dU/da, dU/dk, dU/dh, dU/dα, dU/dβ, dU/dγ]
  3679.      *  </p>
  3680.      */
  3681.     private class FieldUAnddU <T extends RealFieldElement<T>> {

  3683.          /** The current value of the U function. <br/>
  3684.           * Needed for the short periodic contribution */
  3685.         private T U;

  3687.          /** dU / da. */
  3688.         private T dUda;

  3690.          /** dU / dk. */
  3691.         private T dUdk;

  3693.          /** dU / dh. */
  3694.         private T dUdh;

  3696.          /** dU / dAlpha. */
  3697.         private T dUdAl;

  3699.          /** dU / dBeta. */
  3700.         private T dUdBe;

  3702.          /** dU / dGamma. */
  3703.         private T dUdGa;

  3705.          /** Simple constuctor.
  3706.           *  @param date current date
  3707.           *  @param context container for attributes
  3708.           *  @param auxiliaryElements auxiliary elements related to the current orbit
  3709.           *  @param hansen initialization of hansen objects
  3710.           */
  3711.         FieldUAnddU(final FieldAbsoluteDate<T> date,
  3712.                     final FieldDSSTZonalContext<T> context,
  3713.                     final FieldAuxiliaryElements<T> auxiliaryElements,
  3714.                     final FieldHansenObjects<T> hansen) {

  3716.             // Zero for initialization
  3717.             final Field<T> field = date.getField();
  3718.             final T zero = field.getZero();

  3720.             //spherical harmonics
  3721.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date.toAbsoluteDate());

  3723.             //Reset U
  3724.             U = zero;

  3726.             // Gs and Hs coefficients
  3727.             final T[][] GsHs = CoefficientsFactory.computeGsHs(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta(), maxEccPowMeanElements, field);
  3728.             // Qns coefficients
  3729.             final T[][] Qns  = CoefficientsFactory.computeQns(auxiliaryElements.getGamma(), maxDegree, maxEccPowMeanElements);

  3731.             final T[] roaPow = MathArrays.buildArray(field, maxDegree + 1);
  3732.             roaPow[0] = zero.add(1.);
  3733.             for (int i = 1; i <= maxDegree; i++) {
  3734.                 roaPow[i] = roaPow[i - 1].multiply(context.getRoa());
  3735.             }

  3737.             // Potential derivatives
  3738.             dUda  = zero;
  3739.             dUdk  = zero;
  3740.             dUdh  = zero;
  3741.             dUdAl = zero;
  3742.             dUdBe = zero;
  3743.             dUdGa = zero;

  3745.             for (int s = 0; s <= maxEccPowMeanElements; s++) {
  3746.                 //Initialize the Hansen roots
  3747.                 hansen.computeHansenObjectsInitValues(context, s);

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

  3752.                 // Compute Gs partial derivatives from 3.1-(9)
  3753.                 T dGsdh  = zero;
  3754.                 T dGsdk  = zero;
  3755.                 T dGsdAl = zero;
  3756.                 T dGsdBe = zero;
  3757.                 if (s > 0) {
  3758.                     // First get the G(s-1) and the H(s-1) coefficients
  3759.                     final T sxgsm1 = GsHs[0][s - 1].multiply(s);
  3760.                     final T sxhsm1 = GsHs[1][s - 1].multiply(s);
  3761.                     // Then compute derivatives
  3762.                     dGsdh  = sxgsm1.multiply(auxiliaryElements.getBeta()).subtract(sxhsm1.multiply(auxiliaryElements.getAlpha()));
  3763.                     dGsdk  = sxgsm1.multiply(auxiliaryElements.getAlpha()).add(sxhsm1.multiply(auxiliaryElements.getBeta()));
  3764.                     dGsdAl = sxgsm1.multiply(auxiliaryElements.getK()).subtract(sxhsm1.multiply(auxiliaryElements.getH()));
  3765.                     dGsdBe = sxgsm1.multiply(auxiliaryElements.getH()).add(sxhsm1.multiply(auxiliaryElements.getK()));
  3766.                 }

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

  3771.                 for (int n = s + 2; n <= maxDegree; n++) {
  3772.                     // (n - s) must be even
  3773.                     if ((n - s) % 2 == 0) {

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

  3779.                         final double vns   = Vns.get(new NSKey(n, s));
  3780.                         final T coef0 = d0s.multiply(roaPow[n]).multiply(vns).multiply(-harmonics.getUnnormalizedCnm(n, 0));
  3781.                         final T coef1 = coef0.multiply(Qns[n][s]);
  3782.                         final T coef2 = coef1.multiply(kns);
  3783.                         final T coef3 = coef2.multiply(gs);
  3784.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  3785.                         final T dqns  = Qns[n][s + 1];

  3787.                         // Compute U
  3788.                         U = U.add(coef3);
  3789.                         // Compute dU / da :
  3790.                         dUda  = dUda.add(coef3.multiply(n + 1));
  3791.                         // Compute dU / dEx
  3792.                         dUdk  = dUdk.add(coef1.multiply(dGsdk.multiply(kns).add(auxiliaryElements.getK().multiply(context.getXXX()).multiply(dkns).multiply(gs))));
  3793.                         // Compute dU / dEy
  3794.                         dUdh  = dUdh.add(coef1.multiply(dGsdh.multiply(kns).add(auxiliaryElements.getH().multiply(context.getXXX()).multiply(dkns).multiply(gs))));
  3795.                         // Compute dU / dAlpha
  3796.                         dUdAl = dUdAl.add(coef2.multiply(dGsdAl));
  3797.                         // Compute dU / dBeta
  3798.                         dUdBe = dUdBe.add(coef2.multiply(dGsdBe));
  3799.                         // Compute dU / dGamma
  3800.                         dUdGa = dUdGa.add(coef0.multiply(kns).multiply(dqns).multiply(gs));
  3801.                     }
  3802.                 }
  3803.             }

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

  3808.             dUda  = dUda.multiply(context.getMuoa().divide(auxiliaryElements.getSma()));
  3809.             dUdk  = dUdk.multiply(context.getMuoa()).negate();
  3810.             dUdh  = dUdh.multiply(context.getMuoa()).negate();
  3811.             dUdAl = dUdAl.multiply(context.getMuoa()).negate();
  3812.             dUdBe = dUdBe.multiply(context.getMuoa()).negate();
  3813.             dUdGa = dUdGa.multiply(context.getMuoa()).negate();

  3815.         }

  3817.         /** Return value of U.
  3818.          * @return U
  3819.          */
  3820.         public T getU() {
  3821.             return U;
  3822.         }

  3824.          /** Return value of dU / da.
  3825.           * @return dUda
  3826.           */
  3827.         public T getdUda() {
  3828.             return dUda;
  3829.         }

  3831.          /** Return value of dU / dk.
  3832.           * @return dUdk
  3833.           */
  3834.         public T getdUdk() {
  3835.             return dUdk;
  3836.         }

  3838.          /** Return value of dU / dh.
  3839.           * @return dUdh
  3840.           */
  3841.         public T getdUdh() {
  3842.             return dUdh;
  3843.         }

  3845.          /** Return value of dU / dAlpha.
  3846.           * @return dUdAl
  3847.           */
  3848.         public T getdUdAl() {
  3849.             return dUdAl;
  3850.         }

  3852.          /** Return value of dU / dBeta.
  3853.           * @return dUdBe
  3854.           */
  3855.         public T getdUdBe() {
  3856.             return dUdBe;
  3857.         }

  3859.          /** Return value of dU / dGamma.
  3860.           * @return dUdGa
  3861.           */
  3862.         public T getdUdGa() {
  3863.             return dUdGa;
  3864.         }

  3866.     }

  3868.     /** Computes init values of the Hansen Objects. */
  3869.     private class HansenObjects {

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

  3875.         /** Simple constructor. */
  3876.         HansenObjects() {
  3877.             this.hansenObjects = new HansenZonalLinear[maxEccPow + 1];
  3878.             for (int s = 0; s <= maxEccPow; s++) {
  3879.                 this.hansenObjects[s] = new HansenZonalLinear(maxDegree, s);
  3880.             }
  3881.         }

  3883.         /** Compute init values for hansen objects.
  3884.          * @param context container for attributes
  3885.          * @param element element of the array to compute the init values
  3886.          */
  3887.         public void computeHansenObjectsInitValues(final DSSTZonalContext context, final int element) {
  3888.             hansenObjects[element].computeInitValues(context.getX());
  3889.         }

  3891.         /** Get the Hansen Objects.
  3892.          * @return hansenObjects
  3893.          */
  3894.         public HansenZonalLinear[] getHansenObjects() {
  3895.             return hansenObjects;
  3896.         }

  3898.     }

  3900.     /** Computes init values of the Hansen Objects.
  3901.      * @param <T> type of the elements
  3902.      */
  3903.     private class FieldHansenObjects<T extends RealFieldElement<T>> {

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

  3909.         /** Simple constructor.
  3910.          * @param field field used by default
  3911.          */
  3912.         @SuppressWarnings("unchecked")
  3913.         FieldHansenObjects(final Field<T> field) {
  3914.             this.hansenObjects = (FieldHansenZonalLinear<T>[]) Array.newInstance(FieldHansenZonalLinear.class, maxEccPow + 1);
  3915.             for (int s = 0; s <= maxEccPow; s++) {
  3916.                 this.hansenObjects[s] = new FieldHansenZonalLinear<>(maxDegree, s, field);
  3917.             }
  3918.         }

  3920.         /** Compute init values for hansen objects.
  3921.          * @param context container for attributes
  3922.          * @param element element of the array to compute the init values
  3923.          */
  3924.         public void computeHansenObjectsInitValues(final FieldDSSTZonalContext<T> context, final int element) {
  3925.             hansenObjects[element].computeInitValues(context.getX());
  3926.         }

  3928.         /** Get the Hansen Objects.
  3929.          * @return hansenObjects
  3930.          */
  3931.         public FieldHansenZonalLinear<T>[] getHansenObjects() {
  3932.             return hansenObjects;
  3933.         }

  3935.     }

  3937. }
Created with JaCoCo 0.8.6.202009150832