DSSTZonal.java

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

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

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

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

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

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

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

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

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

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

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

  115.     /** Central body rotating frame (fixed with respect to body). */
  116.     private final Frame bodyFixedFrame;

  118.     /** Maximal degree to consider for harmonics potential. */
  119.     private final int maxDegree;

  121.     /** Maximal degree to consider for harmonics potential. */
  122.     private final int maxDegreeShortPeriodics;

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

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

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

  133.     /** Highest power of the eccentricity. */
  134.     private int maxEccPow;

  136.     /** Short period terms. */
  137.     private ZonalShortPeriodicCoefficients zonalSPCoefs;

  139.     /** Short period terms. */
  140.     private final Map<Field<?>, FieldZonalShortPeriodicCoefficients<?>> zonalFieldSPCoefs;

  142.     /** Driver for gravitational parameter. */
  143.     private final ParameterDriver gmParameterDriver;

  145.     /** Hansen objects. */
  146.     private HansenObjects hansen;

  148.     /** Hansen objects for field elements. */
  149.     private final Map<Field<?>, FieldHansenObjects<?>> fieldHansen;

  151.     /** Constructor with default reference values.
  152.      * <p>
  153.      * When this constructor is used, maximum allowed values are used
  154.      * for the short periodic coefficients:
  155.      * </p>
  156.      * <ul>
  157.      *    <li> {@link #maxDegreeShortPeriodics} is set to {@code provider.getMaxDegree()} </li>
  158.      *    <li> {@link #maxEccPowShortPeriodics} is set to {@code min(provider.getMaxDegree() - 1, 4)}.
  159.      *         This parameter should not exceed 4 as higher values will exceed computer capacity </li>
  160.      *    <li> {@link #maxFrequencyShortPeriodics} is set to {@code 2 * provider.getMaxDegree() + 1} </li>
  161.      * </ul>
  162.      * @param bodyFixedFrame rotating body frame
  163.      * <p>If set, this frame will be used to compute the Earth pole direction and the zonal contribution.
  164.      * If not, the inertial frame from the propagator will be used.<br>
  165.      * Ideally this frame should be always set to ITRF.
  166.      * However, using TOD is advised since it's generally a good trade-off between performance and precision.
  167.      *
  168.      * @param provider provider for spherical harmonics
  169.      * @since 12.1
  170.      */
  171.     public DSSTZonal(final Frame bodyFixedFrame, final UnnormalizedSphericalHarmonicsProvider provider) {
  172.         this(bodyFixedFrame, provider, provider.getMaxDegree(), FastMath.min(4, provider.getMaxDegree() - 1), 2 * provider.getMaxDegree() + 1);
  173.     }

  175.     /**
  176.      * Constructor with bodyFixedFrame = orbit frame, and default reference values.
  177.      * <p>
  178.      * Kept for compatibility with anterior versions.
  179.      * </p>
  180.      * <p>
  181.      * Setting bodyFixedFrame to null will lead to large errors if the orbit frame is far from Earth rotating frame (GCRF, EME2000...).
  182.      * The error gets smaller as the orbit frame gets closer to Earth rotating frame (MOD, then TOD).
  183.      * @see <a href="https://gitlab.orekit.org/orekit/orekit/-/issues/1104">issue-1104 on the forge</a>
  184.      * <p>
  185.      * When this constructor is used, maximum allowed values are used
  186.      * for the short periodic coefficients:
  187.      * </p>
  188.      * <ul>
  189.      *    <li> {@link #maxDegreeShortPeriodics} is set to {@code provider.getMaxDegree()} </li>
  190.      *    <li> {@link #maxEccPowShortPeriodics} is set to {@code min(provider.getMaxDegree() - 1, 4)}.
  191.      *         This parameter should not exceed 4 as higher values will exceed computer capacity </li>
  192.      *    <li> {@link #maxFrequencyShortPeriodics} is set to {@code 2 * provider.getMaxDegree() + 1} </li>
  193.      * </ul>
  194.      * @param provider provider for spherical harmonics
  195.      * @see <a href="https://gitlab.orekit.org/orekit/orekit/-/issues/1104">issue-1104 on the forge</a>
  196.      * @since 10.1
  197.      */
  198.     public DSSTZonal(final UnnormalizedSphericalHarmonicsProvider provider) {
  199.         this(null, provider);
  200.     }

  202.     /** Constructor.
  203.      * @param bodyFixedFrame rotating body frame
  204.      * <p>If set, this frame will be used to compute the Earth pole direction and the zonal contribution.
  205.      * If not, the inertial frame from the propagator will be used.<br>
  206.      * Ideally this frame should be always set to ITRF.
  207.      * However, using TOD is advised since it's a good trade-off between performance and precision.
  208.      *
  209.      * @param provider provider for spherical harmonics
  210.      * @param maxDegreeShortPeriodics maximum degree to consider for short periodics zonal harmonics potential
  211.      * (must be between 2 and {@code provider.getMaxDegree()})
  212.      * @param maxEccPowShortPeriodics maximum power of the eccentricity to be used in short periodic computations
  213.      * (must be between 0 and {@code maxDegreeShortPeriodics - 1}, but should typically not exceed 4 as higher
  214.      * values will exceed computer capacity)
  215.      * @param maxFrequencyShortPeriodics maximum frequency in true longitude for short periodic computations
  216.      * (must be between 1 and {@code 2 * maxDegreeShortPeriodics + 1})
  217.      * @since 12.1
  218.      */
  219.     public DSSTZonal(final Frame bodyFixedFrame,
  220.                      final UnnormalizedSphericalHarmonicsProvider provider,
  221.                      final int maxDegreeShortPeriodics,
  222.                      final int maxEccPowShortPeriodics,
  223.                      final int maxFrequencyShortPeriodics) {

  225.         gmParameterDriver = new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
  226.                                                 provider.getMu(), MU_SCALE,
  227.                                                 0.0, Double.POSITIVE_INFINITY);

  229.         // Central body rotating frame
  230.         this.bodyFixedFrame = bodyFixedFrame;

  232.         // Vns coefficients
  233.         this.Vns = CoefficientsFactory.computeVns(provider.getMaxDegree() + 1);

  235.         this.provider  = provider;
  236.         this.maxDegree = provider.getMaxDegree();

  238.         checkIndexRange(maxDegreeShortPeriodics, 2, provider.getMaxDegree());
  239.         this.maxDegreeShortPeriodics = maxDegreeShortPeriodics;

  241.         checkIndexRange(maxEccPowShortPeriodics, 0, maxDegreeShortPeriodics - 1);
  242.         this.maxEccPowShortPeriodics = maxEccPowShortPeriodics;

  244.         checkIndexRange(maxFrequencyShortPeriodics, 1, 2 * maxDegreeShortPeriodics + 1);
  245.         this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;

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

  250.         zonalFieldSPCoefs = new HashMap<>();
  251.         fieldHansen       = new HashMap<>();
  252.     }

  254.     /**
  255.      * Constructor with bodyFixedFrame = orbit frame.
  256.      * <p>
  257.      * Added for retro-compatibility with anterior versions for initialization.
  258.      * </p>
  259.      * <p>
  260.      * Setting bodyFixedFrame to null will lead to large errors if the orbit frame is far from Earth rotating frame (GCRF, EME2000...).
  261.      * The error gets smaller as the orbit frame gets closer to Earth rotating frame (MOD, then TOD).
  262.      * @see <a href="https://gitlab.orekit.org/orekit/orekit/-/issues/1104">issue-1104 on the forge</a>
  263.      * @param provider provider for spherical harmonics
  264.      * @param maxDegreeShortPeriodics maximum degree to consider for short periodics zonal harmonics potential
  265.      * (must be between 2 and {@code provider.getMaxDegree()})
  266.      * @param maxEccPowShortPeriodics maximum power of the eccentricity to be used in short periodic computations
  267.      * (must be between 0 and {@code maxDegreeShortPeriodics - 1}, but should typically not exceed 4 as higher
  268.      * values will exceed computer capacity)
  269.      * @param maxFrequencyShortPeriodics maximum frequency in true longitude for short periodic computations
  270.      * (must be between 1 and {@code 2 * maxDegreeShortPeriodics + 1})
  271.      * @see <a href="https://gitlab.orekit.org/orekit/orekit/-/issues/1104">issue-1104 on the forge</a>
  272.      * @since 7.2
  273.      */
  274.     public DSSTZonal(final UnnormalizedSphericalHarmonicsProvider provider,
  275.                      final int maxDegreeShortPeriodics,
  276.                      final int maxEccPowShortPeriodics,
  277.                      final int maxFrequencyShortPeriodics) {
  278.         this(null, provider, maxDegreeShortPeriodics, maxEccPowShortPeriodics, maxFrequencyShortPeriodics);
  279.     }

  281.     /** Check an index range.
  282.      * @param index index value
  283.      * @param min minimum value for index
  284.      * @param max maximum value for index
  285.      */
  286.     private void checkIndexRange(final int index, final int min, final int max) {
  287.         if (index < min || index > max) {
  288.             throw new OrekitException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, index, min, max);
  289.         }
  290.     }

  292.     /** Get the spherical harmonics provider.
  293.      *  @return the spherical harmonics provider
  294.      */
  295.     public UnnormalizedSphericalHarmonicsProvider getProvider() {
  296.         return provider;
  297.     }

  299.     /** {@inheritDoc}
  300.      *  <p>
  301.      *  Computes the highest power of the eccentricity to appear in the truncated
  302.      *  analytical power series expansion.
  303.      *  </p>
  304.      *  <p>
  305.      *  This method computes the upper value for the central body potential and
  306.      *  determines the maximal power for the eccentricity producing potential
  307.      *  terms bigger than a defined tolerance.
  308.      *  </p>
  309.      */
  310.     @Override
  311.     public List<ShortPeriodTerms> initializeShortPeriodTerms(final AuxiliaryElements auxiliaryElements,
  312.                                              final PropagationType type,
  313.                                              final double[] parameters) {

  315.         computeMeanElementsTruncations(auxiliaryElements, parameters);

  317.         switch (type) {
  318.             case MEAN:
  319.                 maxEccPow = maxEccPowMeanElements;
  320.                 break;
  321.             case OSCULATING:
  322.                 maxEccPow = FastMath.max(maxEccPowMeanElements, maxEccPowShortPeriodics);
  323.                 break;
  324.             default:
  325.                 throw new OrekitInternalError(null);
  326.         }

  328.         hansen = new HansenObjects();

  330.         final List<ShortPeriodTerms> list = new ArrayList<>();
  331.         zonalSPCoefs = new ZonalShortPeriodicCoefficients(maxFrequencyShortPeriodics,
  332.                                                           INTERPOLATION_POINTS,
  333.                                                           new TimeSpanMap<>(new Slot(maxFrequencyShortPeriodics, INTERPOLATION_POINTS)));
  334.         list.add(zonalSPCoefs);
  335.         return list;

  337.     }

  339.     /** {@inheritDoc}
  340.      *  <p>
  341.      *  Computes the highest power of the eccentricity to appear in the truncated
  342.      *  analytical power series expansion.
  343.      *  </p>
  344.      *  <p>
  345.      *  This method computes the upper value for the central body potential and
  346.      *  determines the maximal power for the eccentricity producing potential
  347.      *  terms bigger than a defined tolerance.
  348.      *  </p>
  349.      */
  350.     @Override
  351.     public <T extends CalculusFieldElement<T>> List<FieldShortPeriodTerms<T>> initializeShortPeriodTerms(final FieldAuxiliaryElements<T> auxiliaryElements,
  352.                                                                                      final PropagationType type,
  353.                                                                                      final T[] parameters) {

  355.         // Field used by default
  356.         final Field<T> field = auxiliaryElements.getDate().getField();
  357.         computeMeanElementsTruncations(auxiliaryElements, parameters, field);

  359.         switch (type) {
  360.             case MEAN:
  361.                 maxEccPow = maxEccPowMeanElements;
  362.                 break;
  363.             case OSCULATING:
  364.                 maxEccPow = FastMath.max(maxEccPowMeanElements, maxEccPowShortPeriodics);
  365.                 break;
  366.             default:
  367.                 throw new OrekitInternalError(null);
  368.         }

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

  372.         final FieldZonalShortPeriodicCoefficients<T> fzspc =
  373.                         new FieldZonalShortPeriodicCoefficients<>(maxFrequencyShortPeriodics,
  374.                                                                   INTERPOLATION_POINTS,
  375.                                                                   new FieldTimeSpanMap<>(new FieldSlot<>(maxFrequencyShortPeriodics,
  376.                                                                                                          INTERPOLATION_POINTS),
  377.                                                                                          field));
  378.         zonalFieldSPCoefs.put(field, fzspc);
  379.         return Collections.singletonList(fzspc);

  381.     }

  383.     /** Compute indices truncations for mean elements computations.
  384.      * @param auxiliaryElements auxiliary elements
  385.      * @param parameters values of the force model parameters for state date (only 1 value for each parameter)
  386.      */
  387.     private void computeMeanElementsTruncations(final AuxiliaryElements auxiliaryElements, final double[] parameters) {

  389.         final DSSTZonalContext context = new DSSTZonalContext(auxiliaryElements, bodyFixedFrame, provider, parameters);
  390.         //Compute the max eccentricity power for the mean element rate expansion
  391.         if (maxDegree == 2) {
  392.             maxEccPowMeanElements = 0;
  393.         } else {
  394.             // Initializes specific parameters.
  395.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(auxiliaryElements.getDate());

  397.             // Utilities for truncation
  398.             final double ax2or = 2. * auxiliaryElements.getSma() / provider.getAe();
  399.             double xmuran = parameters[0] / auxiliaryElements.getSma();
  400.             // Set a lower bound for eccentricity
  401.             final double eo2  = FastMath.max(0.0025, auxiliaryElements.getEcc() / 2.);
  402.             final double x2o2 = context.getChi2() / 2.;
  403.             final double[] eccPwr = new double[maxDegree + 1];
  404.             final double[] chiPwr = new double[maxDegree + 1];
  405.             final double[] hafPwr = new double[maxDegree + 1];
  406.             eccPwr[0] = 1.;
  407.             chiPwr[0] = context.getChi();
  408.             hafPwr[0] = 1.;
  409.             for (int i = 1; i <= maxDegree; i++) {
  410.                 eccPwr[i] = eccPwr[i - 1] * eo2;
  411.                 chiPwr[i] = chiPwr[i - 1] * x2o2;
  412.                 hafPwr[i] = hafPwr[i - 1] * 0.5;
  413.                 xmuran  /= ax2or;
  414.             }

  416.             // Set highest power of e and degree of current spherical harmonic.
  417.             maxEccPowMeanElements = 0;
  418.             int maxDeg = maxDegree;
  419.             // Loop over n
  420.             do {
  421.                 // Set order of current spherical harmonic.
  422.                 int m = 0;
  423.                 // Loop over m
  424.                 do {
  425.                     // Compute magnitude of current spherical harmonic coefficient.
  426.                     final double cnm = harmonics.getUnnormalizedCnm(maxDeg, m);
  427.                     final double snm = harmonics.getUnnormalizedSnm(maxDeg, m);
  428.                     final double csnm = FastMath.hypot(cnm, snm);
  429.                     if (csnm == 0.) {
  430.                         break;
  431.                     }
  432.                     // Set magnitude of last spherical harmonic term.
  433.                     double lastTerm = 0.;
  434.                     // Set current power of e and related indices.
  435.                     int nsld2 = (maxDeg - maxEccPowMeanElements - 1) / 2;
  436.                     int l = maxDeg - 2 * nsld2;
  437.                     // Loop over l
  438.                     double term = 0.;
  439.                     do {
  440.                         // Compute magnitude of current spherical harmonic term.
  441.                         if (m < l) {
  442.                             term =  csnm * xmuran *
  443.                                     (CombinatoricsUtils.factorialDouble(maxDeg - l) / (CombinatoricsUtils.factorialDouble(maxDeg - m))) *
  444.                                     (CombinatoricsUtils.factorialDouble(maxDeg + l) / (CombinatoricsUtils.factorialDouble(nsld2) * CombinatoricsUtils.factorialDouble(nsld2 + l))) *
  445.                                     eccPwr[l] * UpperBounds.getDnl(context.getChi2(), chiPwr[l], maxDeg, l) *
  446.                                     (UpperBounds.getRnml(context.getGamma(), maxDeg, l, m, 1, I) + UpperBounds.getRnml(context.getGamma(), maxDeg, l, m, -1, I));
  447.                         } else {
  448.                             term =  csnm * xmuran *
  449.                                     (CombinatoricsUtils.factorialDouble(maxDeg + m) / (CombinatoricsUtils.factorialDouble(nsld2) * CombinatoricsUtils.factorialDouble(nsld2 + l))) *
  450.                                     eccPwr[l] * hafPwr[m - l] * UpperBounds.getDnl(context.getChi2(), chiPwr[l], maxDeg, l) *
  451.                                     (UpperBounds.getRnml(context.getGamma(), maxDeg, m, l, 1, I) + UpperBounds.getRnml(context.getGamma(), maxDeg, m, l, -1, I));
  452.                         }
  453.                         // Is the current spherical harmonic term bigger than the truncation tolerance ?
  454.                         if (term >= TRUNCATION_TOLERANCE) {
  455.                             maxEccPowMeanElements = l;
  456.                         } else {
  457.                             // Is the current term smaller than the last term ?
  458.                             if (term < lastTerm) {
  459.                                 break;
  460.                             }
  461.                         }
  462.                         // Proceed to next power of e.
  463.                         lastTerm = term;
  464.                         l += 2;
  465.                         nsld2--;
  466.                     } while (l < maxDeg);
  467.                     // Is the current spherical harmonic term bigger than the truncation tolerance ?
  468.                     if (term >= TRUNCATION_TOLERANCE) {
  469.                         maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
  470.                         return;
  471.                     }
  472.                     // Proceed to next order.
  473.                     m++;
  474.                 } while (m <= FastMath.min(maxDeg, provider.getMaxOrder()));
  475.                 // Proceed to next degree.
  476.                 xmuran *= ax2or;
  477.                 maxDeg--;
  478.             } while (maxDeg > maxEccPowMeanElements + 2);

  480.             maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
  481.         }
  482.     }

  484.     /** Compute indices truncations for mean elements computations.
  485.      * @param <T> type of the elements
  486.      * @param auxiliaryElements auxiliary elements
  487.      * @param parameters values of the force model parameters
  488.      * @param field field used by default
  489.      */
  490.     private <T extends CalculusFieldElement<T>> void computeMeanElementsTruncations(final FieldAuxiliaryElements<T> auxiliaryElements,
  491.                                                                                     final T[] parameters,
  492.                                                                                     final Field<T> field) {

  494.         final T zero = field.getZero();
  495.         final FieldDSSTZonalContext<T> context = new FieldDSSTZonalContext<>(auxiliaryElements, bodyFixedFrame, provider, parameters);
  496.         //Compute the max eccentricity power for the mean element rate expansion
  497.         if (maxDegree == 2) {
  498.             maxEccPowMeanElements = 0;
  499.         } else {
  500.             // Initializes specific parameters.
  501.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(auxiliaryElements.getDate().toAbsoluteDate());

  503.             // Utilities for truncation
  504.             final T ax2or = auxiliaryElements.getSma().multiply(2.).divide(provider.getAe());
  505.             T xmuran = parameters[0].divide(auxiliaryElements.getSma());
  506.             // Set a lower bound for eccentricity
  507.             final T eo2  = FastMath.max(zero.newInstance(0.0025), auxiliaryElements.getEcc().divide(2.));
  508.             final T x2o2 = context.getChi2().divide(2.);
  509.             final T[] eccPwr = MathArrays.buildArray(field, maxDegree + 1);
  510.             final T[] chiPwr = MathArrays.buildArray(field, maxDegree + 1);
  511.             final T[] hafPwr = MathArrays.buildArray(field, maxDegree + 1);
  512.             eccPwr[0] = zero.newInstance(1.);
  513.             chiPwr[0] = context.getChi();
  514.             hafPwr[0] = zero.newInstance(1.);
  515.             for (int i = 1; i <= maxDegree; i++) {
  516.                 eccPwr[i] = eccPwr[i - 1].multiply(eo2);
  517.                 chiPwr[i] = chiPwr[i - 1].multiply(x2o2);
  518.                 hafPwr[i] = hafPwr[i - 1].multiply(0.5);
  519.                 xmuran  = xmuran.divide(ax2or);
  520.             }

  522.             // Set highest power of e and degree of current spherical harmonic.
  523.             maxEccPowMeanElements = 0;
  524.             int maxDeg = maxDegree;
  525.             // Loop over n
  526.             do {
  527.                 // Set order of current spherical harmonic.
  528.                 int m = 0;
  529.                 // Loop over m
  530.                 do {
  531.                     // Compute magnitude of current spherical harmonic coefficient.
  532.                     final T cnm = zero.newInstance(harmonics.getUnnormalizedCnm(maxDeg, m));
  533.                     final T snm = zero.newInstance(harmonics.getUnnormalizedSnm(maxDeg, m));
  534.                     final T csnm = FastMath.hypot(cnm, snm);
  535.                     if (csnm.getReal() == 0.) {
  536.                         break;
  537.                     }
  538.                     // Set magnitude of last spherical harmonic term.
  539.                     T lastTerm = zero;
  540.                     // Set current power of e and related indices.
  541.                     int nsld2 = (maxDeg - maxEccPowMeanElements - 1) / 2;
  542.                     int l = maxDeg - 2 * nsld2;
  543.                     // Loop over l
  544.                     T term = zero;
  545.                     do {
  546.                         // Compute magnitude of current spherical harmonic term.
  547.                         if (m < l) {
  548.                             term = csnm.multiply(xmuran).
  549.                                    multiply((CombinatoricsUtils.factorialDouble(maxDeg - l) / (CombinatoricsUtils.factorialDouble(maxDeg - m))) *
  550.                                    (CombinatoricsUtils.factorialDouble(maxDeg + l) / (CombinatoricsUtils.factorialDouble(nsld2) * CombinatoricsUtils.factorialDouble(nsld2 + l)))).
  551.                                    multiply(eccPwr[l]).multiply(UpperBounds.getDnl(context.getChi2(), chiPwr[l], maxDeg, l)).
  552.                                    multiply(UpperBounds.getRnml(context.getGamma(), maxDeg, l, m, 1, I).add(UpperBounds.getRnml(context.getGamma(), maxDeg, l, m, -1, I)));
  553.                         } else {
  554.                             term = csnm.multiply(xmuran).
  555.                                    multiply(CombinatoricsUtils.factorialDouble(maxDeg + m) / (CombinatoricsUtils.factorialDouble(nsld2) * CombinatoricsUtils.factorialDouble(nsld2 + l))).
  556.                                    multiply(eccPwr[l]).multiply(hafPwr[m - l]).multiply(UpperBounds.getDnl(context.getChi2(), chiPwr[l], maxDeg, l)).
  557.                                    multiply(UpperBounds.getRnml(context.getGamma(), maxDeg, m, l, 1, I).add(UpperBounds.getRnml(context.getGamma(), maxDeg, m, l, -1, I)));
  558.                         }
  559.                         // Is the current spherical harmonic term bigger than the truncation tolerance ?
  560.                         if (term.getReal() >= TRUNCATION_TOLERANCE) {
  561.                             maxEccPowMeanElements = l;
  562.                         } else {
  563.                             // Is the current term smaller than the last term ?
  564.                             if (term.getReal() < lastTerm.getReal()) {
  565.                                 break;
  566.                             }
  567.                         }
  568.                         // Proceed to next power of e.
  569.                         lastTerm = term;
  570.                         l += 2;
  571.                         nsld2--;
  572.                     } while (l < maxDeg);
  573.                     // Is the current spherical harmonic term bigger than the truncation tolerance ?
  574.                     if (term.getReal() >= TRUNCATION_TOLERANCE) {
  575.                         maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
  576.                         return;
  577.                     }
  578.                     // Proceed to next order.
  579.                     m++;
  580.                 } while (m <= FastMath.min(maxDeg, provider.getMaxOrder()));
  581.                 // Proceed to next degree.
  582.                 xmuran = xmuran.multiply(ax2or);
  583.                 maxDeg--;
  584.             } while (maxDeg > maxEccPowMeanElements + 2);

  586.             maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
  587.         }
  588.     }

  590.     /** Performs initialization at each integration step for the current force model.
  591.      *  <p>
  592.      *  This method aims at being called before mean elements rates computation.
  593.      *  </p>
  594.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  595.      *  @param parameters values of the force model parameters
  596.      *  @return new force model context
  597.      */
  598.     private DSSTZonalContext initializeStep(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
  599.         return new DSSTZonalContext(auxiliaryElements, bodyFixedFrame, provider, parameters);
  600.     }

  602.     /** Performs initialization at each integration step for the current force model.
  603.      *  <p>
  604.      *  This method aims at being called before mean elements rates computation.
  605.      *  </p>
  606.      *  @param <T> type of the elements
  607.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  608.      *  @param parameters values of the force model parameters
  609.      *  @return new force model context
  610.      */
  611.     private <T extends CalculusFieldElement<T>> FieldDSSTZonalContext<T> initializeStep(final FieldAuxiliaryElements<T> auxiliaryElements,
  612.                                                                                     final T[] parameters) {
  613.         return new FieldDSSTZonalContext<>(auxiliaryElements, bodyFixedFrame, provider, parameters);
  614.     }

  616.     /** {@inheritDoc} */
  617.     @Override
  618.     public double[] getMeanElementRate(final SpacecraftState spacecraftState,
  619.                                        final AuxiliaryElements auxiliaryElements, final double[] parameters) {

  621.         // Container of attributes

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

  627.         return computeMeanElementRates(context, udu);

  629.     }

  631.     /** {@inheritDoc} */
  632.     @Override
  633.     public <T extends CalculusFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> spacecraftState,
  634.                                                                   final FieldAuxiliaryElements<T> auxiliaryElements,
  635.                                                                   final T[] parameters) {

  637.         // Field used by default
  638.         final Field<T> field = auxiliaryElements.getDate().getField();
  639.         // Container of attributes
  640.         final FieldDSSTZonalContext<T> context = initializeStep(auxiliaryElements, parameters);

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

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

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

  650.     }

  652.     /** Compute the mean element rates.
  653.      * @param context container for attributes
  654.      * @param udu derivatives of the gravitational potential U
  655.      * @return the mean element rates
  656.      */
  657.     private double[] computeMeanElementRates(final DSSTZonalContext context,
  658.                                              final UAnddU udu) {

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

  663.         // Compute cross derivatives [Eq. 2.2-(8)]
  664.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  665.         final double UAlphaGamma   = context.getAlpha() * udu.getdUdGa() - context.getGamma() * udu.getdUdAl();
  666.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  667.         final double UBetaGamma    =  context.getBeta() * udu.getdUdGa() - context.getGamma() * udu.getdUdBe();
  668.         // Common factor
  669.         final double pUAGmIqUBGoAB = (auxiliaryElements.getP() * UAlphaGamma - I * auxiliaryElements.getQ() * UBetaGamma) * context.getOoAB();

  671.         // Compute mean elements rates [Eq. 3.1-(1)]
  672.         final double da =  0.;
  673.         final double dh =  context.getBoA() * udu.getdUdk() + auxiliaryElements.getK() * pUAGmIqUBGoAB;
  674.         final double dk = -context.getBoA() * udu.getdUdh() - auxiliaryElements.getH() * pUAGmIqUBGoAB;
  675.         final double dp = -context.getCo2AB() * UBetaGamma;
  676.         final double dq = -context.getCo2AB() * UAlphaGamma * I;
  677.         final double dM = -context.getAx2oA() * udu.getdUda() + context.getBoABpo() * (auxiliaryElements.getH() * udu.getdUdh() + auxiliaryElements.getK() * udu.getdUdk()) + pUAGmIqUBGoAB;

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

  682.     /** Compute the mean element rates.
  683.      * @param <T> type of the elements
  684.      * @param date current date
  685.      * @param context container for attributes
  686.      * @param udu derivatives of the gravitational potential U
  687.      * @return the mean element rates
  688.      */
  689.     private <T extends CalculusFieldElement<T>> T[] computeMeanElementRates(final FieldAbsoluteDate<T> date,
  690.                                                                         final FieldDSSTZonalContext<T> context,
  691.                                                                         final FieldUAnddU<T> udu) {

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

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

  699.         // Compute cross derivatives [Eq. 2.2-(8)]
  700.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  701.         final T UAlphaGamma   = udu.getdUdGa().multiply(context.getAlpha()).subtract(udu.getdUdAl().multiply(context.getGamma()));
  702.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  703.         final T UBetaGamma    =  udu.getdUdGa().multiply(context.getBeta()).subtract(udu.getdUdBe().multiply(context.getGamma()));
  704.         // Common factor
  705.         final T pUAGmIqUBGoAB = (UAlphaGamma.multiply(auxiliaryElements.getP()).subtract(UBetaGamma.multiply(I).multiply(auxiliaryElements.getQ()))).multiply(context.getOoAB());

  707.         // Compute mean elements rates [Eq. 3.1-(1)]
  708.         final T da =  field.getZero();
  709.         final T dh =  udu.getdUdk().multiply(context.getBoA()).add(pUAGmIqUBGoAB.multiply(auxiliaryElements.getK()));
  710.         final T dk =  (udu.getdUdh().multiply(context.getBoA()).negate()).subtract(pUAGmIqUBGoAB.multiply(auxiliaryElements.getH()));
  711.         final T dp =  UBetaGamma.multiply(context.getCo2AB().negate());
  712.         final T dq =  UAlphaGamma.multiply(context.getCo2AB().negate()).multiply(I);
  713.         final T dM =  pUAGmIqUBGoAB.add(udu.getdUda().multiply(context.getAx2oA().negate())).add((udu.getdUdh().multiply(auxiliaryElements.getH()).add(udu.getdUdk().multiply(auxiliaryElements.getK()))).multiply(context.getBoABpo()));

  715.         final T[] elements =  MathArrays.buildArray(field, 6);
  716.         elements[0] = da;
  717.         elements[1] = dk;
  718.         elements[2] = dh;
  719.         elements[3] = dq;
  720.         elements[4] = dp;
  721.         elements[5] = dM;

  723.         return elements;
  724.     }

  726.     /** {@inheritDoc} */
  727.     @Override
  728.     public void registerAttitudeProvider(final AttitudeProvider attitudeProvider) {
  729.         //nothing is done since this contribution is not sensitive to attitude
  730.     }

  732.     /** Check if an index is within the accepted interval.
  733.      *
  734.      * @param index the index to check
  735.      * @param lowerBound the lower bound of the interval
  736.      * @param upperBound the upper bound of the interval
  737.      * @return true if the index is between the lower and upper bounds, false otherwise
  738.      */
  739.     private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
  740.         return index >= lowerBound && index <= upperBound;
  741.     }

  743.     /** {@inheritDoc} */
  744.     @Override
  745.     public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {

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

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

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

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

  761.             // Compute rhoj and sigmaj
  762.             final double[][] rhoSigma = computeRhoSigmaCoefficients(slot, auxiliaryElements);

  764.             // Compute Di
  765.             computeDiCoefficients(meanState.getDate(), slot, context, udu);

  767.             // generate the Cij and Sij coefficients
  768.             final FourierCjSjCoefficients cjsj = new FourierCjSjCoefficients(meanState.getDate(),
  769.                                                                              maxDegreeShortPeriodics,
  770.                                                                              maxEccPowShortPeriodics,
  771.                                                                              maxFrequencyShortPeriodics,
  772.                                                                              context,
  773.                                                                              hansen);

  775.             computeCijSijCoefficients(meanState.getDate(), slot, cjsj, rhoSigma, context, auxiliaryElements, udu);
  776.         }

  778.     }

  780.     /** {@inheritDoc} */
  781.     @Override
  782.     @SuppressWarnings("unchecked")
  783.     public <T extends CalculusFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
  784.                                                                        final FieldSpacecraftState<T>... meanStates) {

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

  789.         final FieldZonalShortPeriodicCoefficients<T> fzspc = (FieldZonalShortPeriodicCoefficients<T>) zonalFieldSPCoefs.get(field);
  790.         final FieldSlot<T> slot = fzspc.createSlot(meanStates);
  791.         for (final FieldSpacecraftState<T> meanState : meanStates) {

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

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

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

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

  806.             // Compute rhoj and sigmaj
  807.             final T[][] rhoSigma = computeRhoSigmaCoefficients(slot, auxiliaryElements, field);

  809.             // Compute Di
  810.             computeDiCoefficients(meanState.getDate(), slot, context, field, udu);

  812.             // generate the Cij and Sij coefficients
  813.             final FieldFourierCjSjCoefficients<T> cjsj = new FieldFourierCjSjCoefficients<>(meanState.getDate(),
  814.                                                                                             maxDegreeShortPeriodics,
  815.                                                                                             maxEccPowShortPeriodics,
  816.                                                                                             maxFrequencyShortPeriodics,
  817.                                                                                             context,
  818.                                                                                             fho);


  821.             computeCijSijCoefficients(meanState.getDate(), slot, cjsj, rhoSigma, context, auxiliaryElements, field, udu);
  822.         }
  823.     }

  825.     /** {@inheritDoc} */
  826.     public List<ParameterDriver> getParametersDrivers() {
  827.         return Collections.singletonList(gmParameterDriver);
  828.     }

  830.     /** Generate the values for the D<sub>i</sub> coefficients.
  831.      * @param date target date
  832.      * @param slot slot to which the coefficients belong
  833.      * @param context container for attributes
  834.      * @param udu derivatives of the gravitational potential U
  835.      */
  836.     private void computeDiCoefficients(final AbsoluteDate date,
  837.                                        final Slot slot,
  838.                                        final DSSTZonalContext context,
  839.                                        final UAnddU udu) {

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

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

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

  849.             if (i == 5) {
  850.                 currentDi[i] += -1.5 * 2 * udu.getU() * context.getOON2A2();
  851.             }

  853.         }

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

  857.     }

  859.     /** Generate the values for the D<sub>i</sub> coefficients.
  860.      * @param <T> type of the elements
  861.      * @param date target date
  862.      * @param slot slot to which the coefficients belong
  863.      * @param context container for attributes
  864.      * @param field field used by default
  865.      * @param udu derivatives of the gravitational potential U
  866.      */
  867.     private <T extends CalculusFieldElement<T>> void computeDiCoefficients(final FieldAbsoluteDate<T> date,
  868.                                                                        final FieldSlot<T> slot,
  869.                                                                        final FieldDSSTZonalContext<T> context,
  870.                                                                        final Field<T> field,
  871.                                                                        final FieldUAnddU<T> udu) {

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

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

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

  881.             if (i == 5) {
  882.                 currentDi[i] = currentDi[i].add(context.getOON2A2().multiply(udu.getU()).multiply(2.).multiply(-1.5));
  883.             }

  885.         }

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

  889.     }

  891.     /**
  892.      * Generate the values for the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
  893.      * @param date date of computation
  894.      * @param slot slot to which the coefficients belong
  895.      * @param cjsj Fourier coefficients
  896.      * @param rhoSigma ρ<sub>j</sub> and σ<sub>j</sub>
  897.      * @param context container for attributes
  898.      * @param auxiliaryElements auxiliary elements related to the current orbit
  899.      * @param udu derivatives of the gravitational potential U
  900.      */
  901.     private void computeCijSijCoefficients(final AbsoluteDate date, final Slot slot,
  902.                                            final FourierCjSjCoefficients cjsj,
  903.                                            final double[][] rhoSigma, final DSSTZonalContext context,
  904.                                            final AuxiliaryElements auxiliaryElements,
  905.                                            final UAnddU udu) {

  907.         final int nMax = maxDegreeShortPeriodics;

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

  913.             // Create local arrays
  914.             final double[] currentCij = new double[] {0., 0., 0., 0., 0., 0.};
  915.             final double[] currentSij = new double[] {0., 0., 0., 0., 0., 0.};

  917.             // j == 1
  918.             if (j == 1) {
  919.                 final double coef1 = 4 * auxiliaryElements.getK() * udu.getU() - context.getHK() * cjsj.getCj(1) + context.getK2MH2O2() * cjsj.getSj(1);
  920.                 final double coef2 = 4 * auxiliaryElements.getH() * udu.getU() + context.getK2MH2O2() * cjsj.getCj(1) + context.getHK() * cjsj.getSj(1);
  921.                 final double coef3 = (auxiliaryElements.getK() * cjsj.getCj(1) + auxiliaryElements.getH() * cjsj.getSj(1)) / 4.;
  922.                 final double coef4 = (8 * udu.getU() - auxiliaryElements.getH() * cjsj.getCj(1) + auxiliaryElements.getK() * cjsj.getSj(1)) / 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 == 2
  942.             if (j == 2) {
  943.                 final double coef1 = context.getK2MH2() * udu.getU();
  944.                 final double coef2 = 2 * context.getHK() * udu.getU();
  945.                 final double coef3 = auxiliaryElements.getH() * udu.getU() / 2;
  946.                 final double coef4 = auxiliaryElements.getK() * udu.getU() / 2;

  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 1 and 2N-3
  966.             if (isBetween(j, 1, 2 * nMax - 3) && j + 2 < cjsj.jMax) {
  967.                 final double coef1 = ( j + 2 ) * (-context.getHK() * cjsj.getCj(j + 2) + context.getK2MH2O2() * cjsj.getSj(j + 2));
  968.                 final double coef2 = ( j + 2 ) * (context.getK2MH2O2() * cjsj.getCj(j + 2) + context.getHK() * cjsj.getSj(j + 2));
  969.                 final double coef3 = ( j + 2 ) * (auxiliaryElements.getK() * cjsj.getCj(j + 2) + auxiliaryElements.getH() * cjsj.getSj(j + 2)) / 4;
  970.                 final double coef4 = ( j + 2 ) * (auxiliaryElements.getH() * cjsj.getCj(j + 2) - auxiliaryElements.getK() * cjsj.getSj(j + 2)) / 4;

  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 1 and 2N-2
  990.             if (isBetween(j, 1, 2 * nMax - 2) && j + 1 < cjsj.jMax) {
  991.                 final double coef1 = 2 * ( j + 1 ) * (-auxiliaryElements.getH() * cjsj.getCj(j + 1) + auxiliaryElements.getK() * cjsj.getSj(j + 1));
  992.                 final double coef2 = 2 * ( j + 1 ) * (auxiliaryElements.getK() * cjsj.getCj(j + 1) + auxiliaryElements.getH() * cjsj.getSj(j + 1));
  993.                 final double coef3 = ( j + 1 ) * cjsj.getCj(j + 1);
  994.                 final double coef4 = ( j + 1 ) * cjsj.getSj(j + 1);

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

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

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

  1008.                 //Coefficients for λ
  1009.                 currentCij[5] -= coef2 / 2;
  1010.                 currentSij[5] += -coef1 / 2;
  1011.             }

  1013.             // j between 2 and 2N
  1014.             if (isBetween(j, 2, 2 * nMax) && j - 1 < cjsj.jMax) {
  1015.                 final double coef1 = 2 * ( j - 1 ) * (auxiliaryElements.getH() * cjsj.getCj(j - 1) + auxiliaryElements.getK() * cjsj.getSj(j - 1));
  1016.                 final double coef2 = 2 * ( j - 1 ) * (auxiliaryElements.getK() * cjsj.getCj(j - 1) - auxiliaryElements.getH() * cjsj.getSj(j - 1));
  1017.                 final double coef3 = ( j - 1 ) * cjsj.getCj(j - 1);
  1018.                 final double coef4 = ( j - 1 ) * cjsj.getSj(j - 1);

  1020.                 //Coefficients for a
  1021.                 currentCij[0] += coef1;
  1022.                 currentSij[0] -= coef2;

  1024.                 //Coefficients for k
  1025.                 currentCij[1] += coef4;
  1026.                 currentSij[1] -= coef3;

  1028.                 //Coefficients for h
  1029.                 currentCij[2] += coef3;
  1030.                 currentSij[2] += coef4;

  1032.                 //Coefficients for λ
  1033.                 currentCij[5] += coef2 / 2;
  1034.                 currentSij[5] += coef1 / 2;
  1035.             }

  1037.             // j between 3 and 2N + 1
  1038.             if (isBetween(j, 3, 2 * nMax + 1) && j - 2 < cjsj.jMax) {
  1039.                 final double coef1 = ( j - 2 ) * (context.getHK() * cjsj.getCj(j - 2) + context.getK2MH2O2() * cjsj.getSj(j - 2));
  1040.                 final double coef2 = ( j - 2 ) * (-context.getK2MH2O2() * cjsj.getCj(j - 2) + context.getHK() * cjsj.getSj(j - 2));
  1041.                 final double coef3 = ( j - 2 ) * (auxiliaryElements.getK() * cjsj.getCj(j - 2) - auxiliaryElements.getH() * cjsj.getSj(j - 2)) / 4;
  1042.                 final double coef4 = ( j - 2 ) * (auxiliaryElements.getH() * cjsj.getCj(j - 2) + auxiliaryElements.getK() * cjsj.getSj(j - 2)) / 4;
  1043.                 final double coef5 = ( j - 2 ) * (context.getK2MH2O2() * cjsj.getCj(j - 2) - context.getHK() * cjsj.getSj(j - 2));

  1045.                 //Coefficients for a
  1046.                 currentCij[0] += coef1;
  1047.                 currentSij[0] += coef2;

  1049.                 //Coefficients for k
  1050.                 currentCij[1] += coef4;
  1051.                 currentSij[1] -= coef3;

  1053.                 //Coefficients for h
  1054.                 currentCij[2] += coef3;
  1055.                 currentSij[2] += coef4;

  1057.                 //Coefficients for λ
  1058.                 currentCij[5] += coef5 / 2;
  1059.                 currentSij[5] += coef1 / 2;
  1060.             }

  1062.             //multiply by the common factor
  1063.             //for a (i == 0) -> χ³ / (n² * a)
  1064.             currentCij[0] *= context.getX3ON2A();
  1065.             currentSij[0] *= context.getX3ON2A();
  1066.             //for k (i == 1) -> χ / (n² * a²)
  1067.             currentCij[1] *= context.getXON2A2();
  1068.             currentSij[1] *= context.getXON2A2();
  1069.             //for h (i == 2) -> χ / (n² * a²)
  1070.             currentCij[2] *= context.getXON2A2();
  1071.             currentSij[2] *= context.getXON2A2();
  1072.             //for λ (i == 5) -> (χ²) / (n² * a² * (χ + 1 ) )
  1073.             currentCij[5] *= context.getX2ON2A2XP1();
  1074.             currentSij[5] *= context.getX2ON2A2XP1();

  1076.             // j is between 1 and 2 * N - 1
  1077.             if (isBetween(j, 1, 2 * nMax - 1) && j < cjsj.jMax) {
  1078.                 // Compute cross derivatives
  1079.                 // Cj(alpha,gamma) = alpha * dC/dgamma - gamma * dC/dalpha
  1080.                 final double CjAlphaGamma   = context.getAlpha() * cjsj.getdCjdGamma(j) - context.getGamma() * cjsj.getdCjdAlpha(j);
  1081.                 // Cj(alpha,beta) = alpha * dC/dbeta - beta * dC/dalpha
  1082.                 final double CjAlphaBeta   = context.getAlpha() * cjsj.getdCjdBeta(j) - context.getBeta() * cjsj.getdCjdAlpha(j);
  1083.                 // Cj(beta,gamma) = beta * dC/dgamma - gamma * dC/dbeta
  1084.                 final double CjBetaGamma    =  context.getBeta() * cjsj.getdCjdGamma(j) - context.getGamma() * cjsj.getdCjdBeta(j);
  1085.                 // Cj(h,k) = h * dC/dk - k * dC/dh
  1086.                 final double CjHK   = auxiliaryElements.getH() * cjsj.getdCjdK(j) - auxiliaryElements.getK() * cjsj.getdCjdH(j);
  1087.                 // Sj(alpha,gamma) = alpha * dS/dgamma - gamma * dS/dalpha
  1088.                 final double SjAlphaGamma   = context.getAlpha() * cjsj.getdSjdGamma(j) - context.getGamma() * cjsj.getdSjdAlpha(j);
  1089.                 // Sj(alpha,beta) = alpha * dS/dbeta - beta * dS/dalpha
  1090.                 final double SjAlphaBeta   = context.getAlpha() * cjsj.getdSjdBeta(j) - context.getBeta() * cjsj.getdSjdAlpha(j);
  1091.                 // Sj(beta,gamma) = beta * dS/dgamma - gamma * dS/dbeta
  1092.                 final double SjBetaGamma    =  context.getBeta() * cjsj.getdSjdGamma(j) - context.getGamma() * cjsj.getdSjdBeta(j);
  1093.                 // Sj(h,k) = h * dS/dk - k * dS/dh
  1094.                 final double SjHK   = auxiliaryElements.getH() * cjsj.getdSjdK(j) - auxiliaryElements.getK() * cjsj.getdSjdH(j);

  1096.                 //Coefficients for a
  1097.                 final double coef1 = context.getX3ON2A() * (3 - context.getBB()) * j;
  1098.                 currentCij[0] += coef1 * cjsj.getSj(j);
  1099.                 currentSij[0] -= coef1 * cjsj.getCj(j);

  1101.                 //Coefficients for k and h
  1102.                 final double coef2 = auxiliaryElements.getP() * CjAlphaGamma - I * auxiliaryElements.getQ() * CjBetaGamma;
  1103.                 final double coef3 = auxiliaryElements.getP() * SjAlphaGamma - I * auxiliaryElements.getQ() * SjBetaGamma;
  1104.                 currentCij[1] -= context.getXON2A2() * (auxiliaryElements.getH() * coef2 + context.getBB() * cjsj.getdCjdH(j) - 1.5 * auxiliaryElements.getK() * j * cjsj.getSj(j));
  1105.                 currentSij[1] -= context.getXON2A2() * (auxiliaryElements.getH() * coef3 + context.getBB() * cjsj.getdSjdH(j) + 1.5 * auxiliaryElements.getK() * j * cjsj.getCj(j));
  1106.                 currentCij[2] += context.getXON2A2() * (auxiliaryElements.getK() * coef2 + context.getBB() * cjsj.getdCjdK(j) + 1.5 * auxiliaryElements.getH() * j * cjsj.getSj(j));
  1107.                 currentSij[2] += context.getXON2A2() * (auxiliaryElements.getK() * coef3 + context.getBB() * cjsj.getdSjdK(j) - 1.5 * auxiliaryElements.getH() * j * cjsj.getCj(j));

  1109.                 //Coefficients for q and p
  1110.                 final double coef4 = CjHK - CjAlphaBeta - j * cjsj.getSj(j);
  1111.                 final double coef5 = SjHK - SjAlphaBeta + j * cjsj.getCj(j);
  1112.                 currentCij[3] = context.getCXO2N2A2() * (-I * CjAlphaGamma + auxiliaryElements.getQ() * coef4);
  1113.                 currentSij[3] = context.getCXO2N2A2() * (-I * SjAlphaGamma + auxiliaryElements.getQ() * coef5);
  1114.                 currentCij[4] = context.getCXO2N2A2() * (-CjBetaGamma + auxiliaryElements.getP() * coef4);
  1115.                 currentSij[4] = context.getCXO2N2A2() * (-SjBetaGamma + auxiliaryElements.getP() * coef5);

  1117.                 //Coefficients for λ
  1118.                 final double coef6 = auxiliaryElements.getH() * cjsj.getdCjdH(j) + auxiliaryElements.getK() * cjsj.getdCjdK(j);
  1119.                 final double coef7 = auxiliaryElements.getH() * cjsj.getdSjdH(j) + auxiliaryElements.getK() * cjsj.getdSjdK(j);
  1120.                 currentCij[5] += context.getOON2A2() * (-2 * auxiliaryElements.getSma() * cjsj.getdCjdA(j) + coef6 / (context.getChi() + 1) + context.getChi() * coef2 - 3 * cjsj.getCj(j));
  1121.                 currentSij[5] += context.getOON2A2() * (-2 * auxiliaryElements.getSma() * cjsj.getdSjdA(j) + coef7 / (context.getChi() + 1) + context.getChi() * coef3 - 3 * cjsj.getSj(j));
  1122.             }

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

  1129.             // Add the coefficients to the interpolation grid
  1130.             slot.cij[j].addGridPoint(date, currentCij);
  1131.             slot.sij[j].addGridPoint(date, currentSij);

  1133.         }

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

  1138.     }

  1140.     /**
  1141.      * Generate the values for the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
  1142.      * @param <T> type of the elements
  1143.      * @param date date of computation
  1144.      * @param slot slot to which the coefficients belong
  1145.      * @param cjsj Fourier coefficients
  1146.      * @param rhoSigma ρ<sub>j</sub> and σ<sub>j</sub>
  1147.      * @param context container for attributes
  1148.      * @param auxiliaryElements auxiliary elements related to the current orbit
  1149.      * @param field field used by default
  1150.      * @param udu derivatives of the gravitational potential U
  1151.      */
  1152.     private <T extends CalculusFieldElement<T>> void computeCijSijCoefficients(final FieldAbsoluteDate<T> date,
  1153.                                                                            final FieldSlot<T> slot,
  1154.                                                                            final FieldFourierCjSjCoefficients<T> cjsj,
  1155.                                                                            final T[][] rhoSigma,
  1156.                                                                            final FieldDSSTZonalContext<T> context,
  1157.                                                                            final FieldAuxiliaryElements<T> auxiliaryElements,
  1158.                                                                            final Field<T> field,
  1159.                                                                            final FieldUAnddU<T> udu) {

  1161.         // Zero
  1162.         final T zero = field.getZero();

  1164.         final int nMax = maxDegreeShortPeriodics;

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

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

  1172.             // Create local arrays
  1173.             final T[] currentCij = MathArrays.buildArray(field, 6);
  1174.             final T[] currentSij = MathArrays.buildArray(field, 6);

  1176.             Arrays.fill(currentCij, zero);
  1177.             Arrays.fill(currentSij, zero);

  1179.             // j == 1
  1180.             if (j == 1) {
  1181.                 final T coef1 = auxiliaryElements.getK().multiply(udu.getU()).multiply(4.).subtract(context.getHK().multiply(cjsj.getCj(1))).add(context.getK2MH2O2().multiply(cjsj.getSj(1)));
  1182.                 final T coef2 = auxiliaryElements.getH().multiply(udu.getU()).multiply(4.).add(context.getK2MH2O2().multiply(cjsj.getCj(1))).add(context.getHK().multiply(cjsj.getSj(1)));
  1183.                 final T coef3 = auxiliaryElements.getK().multiply(cjsj.getCj(1)).add(auxiliaryElements.getH().multiply(cjsj.getSj(1))).divide(4.);
  1184.                 final T coef4 = udu.getU().multiply(8.).subtract(auxiliaryElements.getH().multiply(cjsj.getCj(1))).add(auxiliaryElements.getK().multiply(cjsj.getSj(1))).divide(4.);

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

  1190.                 //Coefficients for k
  1191.                 currentCij[1] = currentCij[1].add(coef4);
  1192.                 currentSij[1] = currentSij[1].add(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.divide(2.));
  1201.             }

  1203.             // j == 2
  1204.             if (j == 2) {
  1205.                 final T coef1 = context.getK2MH2().multiply(udu.getU());
  1206.                 final T coef2 = context.getHK().multiply(udu.getU()).multiply(2.);
  1207.                 final T coef3 = auxiliaryElements.getH().multiply(udu.getU()).divide(2.);
  1208.                 final T coef4 = auxiliaryElements.getK().multiply(udu.getU()).divide(2.);

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

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

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

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

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

  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.negate());
  1240.                 currentSij[1] = currentSij[1].subtract(coef3);

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

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

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

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

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

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

  1270.                 //Coefficients for λ
  1271.                 currentCij[5] = currentCij[5].subtract(coef2.divide(2.));
  1272.                 currentSij[5] = currentSij[5].add(coef1.negate().divide(2.));
  1273.             }

  1275.             // j between 2 and 2N
  1276.             if (isBetween(j, 2, 2 * nMax) && j - 1 < cjsj.jMax) {
  1277.                 final T coef1 = auxiliaryElements.getH().multiply(cjsj.getCj(j - 1)).add(auxiliaryElements.getK().multiply(cjsj.getSj(j - 1))).multiply(2 * ( j - 1 ));
  1278.                 final T coef2 = auxiliaryElements.getK().multiply(cjsj.getCj(j - 1)).subtract(auxiliaryElements.getH().multiply(cjsj.getSj(j - 1))).multiply(2 * ( j - 1 ));
  1279.                 final T coef3 = cjsj.getCj(j - 1).multiply(j - 1);
  1280.                 final T coef4 = cjsj.getSj(j - 1).multiply(j - 1);

  1282.                 //Coefficients for a
  1283.                 currentCij[0] = currentCij[0].add(coef1);
  1284.                 currentSij[0] = currentSij[0].subtract(coef2);

  1286.                 //Coefficients for k
  1287.                 currentCij[1] = currentCij[1].add(coef4);
  1288.                 currentSij[1] = currentSij[1].subtract(coef3);

  1290.                 //Coefficients for h
  1291.                 currentCij[2] = currentCij[2].add(coef3);
  1292.                 currentSij[2] = currentSij[2].add(coef4);

  1294.                 //Coefficients for λ
  1295.                 currentCij[5] = currentCij[5].add(coef2.divide(2.));
  1296.                 currentSij[5] = currentSij[5].add(coef1.divide(2.));
  1297.             }

  1299.             // j between 3 and 2N + 1
  1300.             if (isBetween(j, 3, 2 * nMax + 1) && j - 2 < cjsj.jMax) {
  1301.                 final T coef1 = context.getHK().multiply(cjsj.getCj(j - 2)).add(context.getK2MH2O2().multiply(cjsj.getSj(j - 2))).multiply(j - 2);
  1302.                 final T coef2 = context.getK2MH2O2().negate().multiply(cjsj.getCj(j - 2)).add(context.getHK().multiply(cjsj.getSj(j - 2))).multiply(j - 2);
  1303.                 final T coef3 = auxiliaryElements.getK().multiply(cjsj.getCj(j - 2)).subtract(auxiliaryElements.getH().multiply(cjsj.getSj(j - 2))).multiply(j - 2).divide(4.);
  1304.                 final T coef4 = auxiliaryElements.getH().multiply(cjsj.getCj(j - 2)).add(auxiliaryElements.getK().multiply(cjsj.getSj(j - 2))).multiply(j - 2).divide(4.);
  1305.                 final T coef5 = context.getK2MH2O2().multiply(cjsj.getCj(j - 2)).subtract(context.getHK().multiply(cjsj.getSj(j - 2))).multiply(j - 2);

  1307.                 //Coefficients for a
  1308.                 currentCij[0] = currentCij[0].add(coef1);
  1309.                 currentSij[0] = currentSij[0].add(coef2);

  1311.                 //Coefficients for k
  1312.                 currentCij[1] = currentCij[1].add(coef4);
  1313.                 currentSij[1] = currentSij[1].add(coef3.negate());

  1315.                 //Coefficients for h
  1316.                 currentCij[2] = currentCij[2].add(coef3);
  1317.                 currentSij[2] = currentSij[2].add(coef4);

  1319.                 //Coefficients for λ
  1320.                 currentCij[5] = currentCij[5].add(coef5.divide(2.));
  1321.                 currentSij[5] = currentSij[5].add(coef1.divide(2.));
  1322.             }

  1324.             //multiply by the common factor
  1325.             //for a (i == 0) -> χ³ / (n² * a)
  1326.             currentCij[0] = currentCij[0].multiply(context.getX3ON2A());
  1327.             currentSij[0] = currentSij[0].multiply(context.getX3ON2A());
  1328.             //for k (i == 1) -> χ / (n² * a²)
  1329.             currentCij[1] = currentCij[1].multiply(context.getXON2A2());
  1330.             currentSij[1] = currentSij[1].multiply(context.getXON2A2());
  1331.             //for h (i == 2) -> χ / (n² * a²)
  1332.             currentCij[2] = currentCij[2].multiply(context.getXON2A2());
  1333.             currentSij[2] = currentSij[2].multiply(context.getXON2A2());
  1334.             //for λ (i == 5) -> (χ²) / (n² * a² * (χ + 1 ) )
  1335.             currentCij[5] = currentCij[5].multiply(context.getX2ON2A2XP1());
  1336.             currentSij[5] = currentSij[5].multiply(context.getX2ON2A2XP1());

  1338.             // j is between 1 and 2 * N - 1
  1339.             if (isBetween(j, 1, 2 * nMax - 1) && j < cjsj.jMax) {
  1340.                 // Compute cross derivatives
  1341.                 // Cj(alpha,gamma) = alpha * dC/dgamma - gamma * dC/dalpha
  1342.                 final T CjAlphaGamma = context.getAlpha().multiply(cjsj.getdCjdGamma(j)).subtract(context.getGamma().multiply(cjsj.getdCjdAlpha(j)));
  1343.                 // Cj(alpha,beta) = alpha * dC/dbeta - beta * dC/dalpha
  1344.                 final T CjAlphaBeta  = context.getAlpha().multiply(cjsj.getdCjdBeta(j)).subtract(context.getBeta().multiply(cjsj.getdCjdAlpha(j)));
  1345.                 // Cj(beta,gamma) = beta * dC/dgamma - gamma * dC/dbeta
  1346.                 final T CjBetaGamma  = context.getBeta().multiply(cjsj.getdCjdGamma(j)).subtract(context.getGamma().multiply(cjsj.getdCjdBeta(j)));
  1347.                 // Cj(h,k) = h * dC/dk - k * dC/dh
  1348.                 final T CjHK         = auxiliaryElements.getH().multiply(cjsj.getdCjdK(j)).subtract(auxiliaryElements.getK().multiply(cjsj.getdCjdH(j)));
  1349.                 // Sj(alpha,gamma) = alpha * dS/dgamma - gamma * dS/dalpha
  1350.                 final T SjAlphaGamma = context.getAlpha().multiply(cjsj.getdSjdGamma(j)).subtract(context.getGamma().multiply(cjsj.getdSjdAlpha(j)));
  1351.                 // Sj(alpha,beta) = alpha * dS/dbeta - beta * dS/dalpha
  1352.                 final T SjAlphaBeta  = context.getAlpha().multiply(cjsj.getdSjdBeta(j)).subtract(context.getBeta().multiply(cjsj.getdSjdAlpha(j)));
  1353.                 // Sj(beta,gamma) = beta * dS/dgamma - gamma * dS/dbeta
  1354.                 final T SjBetaGamma  = context.getBeta().multiply(cjsj.getdSjdGamma(j)).subtract(context.getGamma().multiply(cjsj.getdSjdBeta(j)));
  1355.                 // Sj(h,k) = h * dS/dk - k * dS/dh
  1356.                 final T SjHK         = auxiliaryElements.getH().multiply(cjsj.getdSjdK(j)).subtract(auxiliaryElements.getK().multiply(cjsj.getdSjdH(j)));

  1358.                 //Coefficients for a
  1359.                 final T coef1 = context.getX3ON2A().multiply(context.getBB().negate().add(3.)).multiply(j);
  1360.                 currentCij[0] = currentCij[0].add(coef1.multiply(cjsj.getSj(j)));
  1361.                 currentSij[0] = currentSij[0].subtract(coef1.multiply(cjsj.getCj(j)));

  1363.                 //Coefficients for k and h
  1364.                 final T coef2 = auxiliaryElements.getP().multiply(CjAlphaGamma).subtract(auxiliaryElements.getQ().multiply(CjBetaGamma).multiply(I));
  1365.                 final T coef3 = auxiliaryElements.getP().multiply(SjAlphaGamma).subtract(auxiliaryElements.getQ().multiply(SjBetaGamma).multiply(I));
  1366.                 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)))));
  1367.                 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)))));
  1368.                 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)))));
  1369.                 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)))));

  1371.                 //Coefficients for q and p
  1372.                 final T coef4 = CjHK.subtract(CjAlphaBeta).subtract(cjsj.getSj(j).multiply(j));
  1373.                 final T coef5 = SjHK.subtract(SjAlphaBeta).add(cjsj.getCj(j).multiply(j));
  1374.                 currentCij[3] = context.getCXO2N2A2().multiply(CjAlphaGamma.multiply(-I).add(auxiliaryElements.getQ().multiply(coef4)));
  1375.                 currentSij[3] = context.getCXO2N2A2().multiply(SjAlphaGamma.multiply(-I).add(auxiliaryElements.getQ().multiply(coef5)));
  1376.                 currentCij[4] = context.getCXO2N2A2().multiply(CjBetaGamma.negate().add(auxiliaryElements.getP().multiply(coef4)));
  1377.                 currentSij[4] = context.getCXO2N2A2().multiply(SjBetaGamma.negate().add(auxiliaryElements.getP().multiply(coef5)));

  1379.                 //Coefficients for λ
  1380.                 final T coef6 = auxiliaryElements.getH().multiply(cjsj.getdCjdH(j)).add(auxiliaryElements.getK().multiply(cjsj.getdCjdK(j)));
  1381.                 final T coef7 = auxiliaryElements.getH().multiply(cjsj.getdSjdH(j)).add(auxiliaryElements.getK().multiply(cjsj.getdSjdK(j)));
  1382.                 currentCij[5] = currentCij[5].add(context.getOON2A2().multiply(auxiliaryElements.getSma().multiply(-2.).multiply(cjsj.getdCjdA(j)).add(coef6.divide(context.getChi().add(1.))).add(context.getChi().multiply(coef2)).subtract(cjsj.getCj(j).multiply(3.))));
  1383.                 currentSij[5] = currentSij[5].add(context.getOON2A2().multiply(auxiliaryElements.getSma().multiply(-2.).multiply(cjsj.getdSjdA(j)).add(coef7.divide(context.getChi().add(1.))).add(context.getChi().multiply(coef3)).subtract(cjsj.getSj(j).multiply(3.))));
  1384.             }

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

  1391.             // Add the coefficients to the interpolation grid
  1392.             slot.cij[j].addGridPoint(date, currentCij);
  1393.             slot.sij[j].addGridPoint(date, currentSij);

  1395.         }

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

  1400.     }

  1402.     /**
  1403.      * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
  1404.      * <p>
  1405.      * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
  1406.      *  ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
  1407.      *  σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
  1408.      * </p>
  1409.      * @param slot slot to which the coefficients belong
  1410.      * @param auxiliaryElements auxiliary elements related to the current orbit
  1411.      * @return array containing ρ<sub>j</sub> and σ<sub>j</sub>
  1412.      */
  1413.     private double[][] computeRhoSigmaCoefficients(final Slot slot,
  1414.                                                    final AuxiliaryElements auxiliaryElements) {

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

  1419.         // (-b)<sup>j</sup>
  1420.         double mbtj = 1;

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

  1425.             //Compute current rho and sigma;
  1426.             mbtj *= -b;
  1427.             final double coef  = (1 + j * auxiliaryElements.getB()) * mbtj;
  1428.             final double rho   = coef * cjsjKH.getCj(j);
  1429.             final double sigma = coef * cjsjKH.getSj(j);

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

  1436.         return rhoSigma;

  1438.     }

  1440.     /**
  1441.      * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
  1442.      * <p>
  1443.      * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
  1444.      *  ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
  1445.      *  σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
  1446.      * </p>
  1447.      * @param <T> type of the elements
  1448.      * @param slot slot to which the coefficients belong
  1449.      * @param auxiliaryElements auxiliary elements related to the current orbit
  1450.      * @param field field used by default
  1451.      * @return array containing ρ<sub>j</sub> and σ<sub>j</sub>
  1452.      */
  1453.     private <T extends CalculusFieldElement<T>> T[][] computeRhoSigmaCoefficients(final FieldSlot<T> slot,
  1454.                                                                                   final FieldAuxiliaryElements<T> auxiliaryElements,
  1455.                                                                                   final Field<T> field) {
  1456.         final T zero = field.getZero();

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

  1461.         // (-b)<sup>j</sup>
  1462.         T mbtj = zero.newInstance(1.);

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

  1467.             //Compute current rho and sigma;
  1468.             mbtj = mbtj.multiply(b.negate());
  1469.             final T coef  = mbtj.multiply(auxiliaryElements.getB().multiply(j).add(1.));
  1470.             final T rho   = coef.multiply(cjsjKH.getCj(j));
  1471.             final T sigma = coef.multiply(cjsjKH.getSj(j));

  1473.             // Add the coefficients to the interpolation grid
  1474.             rhoSigma[j][0] = rho;
  1475.             rhoSigma[j][1] = sigma;
  1476.         }

  1478.         return rhoSigma;

  1480.     }

  1482.     /** The coefficients used to compute the short-periodic zonal contribution.
  1483.      *
  1484.      * <p>
  1485.      * Those coefficients are given in Danielson paper by expressions 4.1-(20) to 4.1.-(25)
  1486.      * </p>
  1487.      * <p>
  1488.      * The coefficients are: <br>
  1489.      * - C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> <br>
  1490.      * - ρ<sub>j</sub> and σ<sub>j</sub> <br>
  1491.      * - C<sub>i</sub>⁰
  1492.      * </p>
  1493.      *
  1494.      * @author Lucian Barbulescu
  1495.      */
  1496.     private static class ZonalShortPeriodicCoefficients implements ShortPeriodTerms {

  1498.         /** Maximum value for j index. */
  1499.         private final int maxFrequencyShortPeriodics;

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

  1504.         /** All coefficients slots. */
  1505.         private final transient TimeSpanMap<Slot> slots;

  1507.         /** Constructor.
  1508.          * @param maxFrequencyShortPeriodics maximum value for j index
  1509.          * @param interpolationPoints number of points used in the interpolation process
  1510.          * @param slots all coefficients slots
  1511.          */
  1512.         ZonalShortPeriodicCoefficients(final int maxFrequencyShortPeriodics, final int interpolationPoints,
  1513.                                        final TimeSpanMap<Slot> slots) {

  1515.             // Save parameters
  1516.             this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
  1517.             this.interpolationPoints        = interpolationPoints;
  1518.             this.slots                      = slots;

  1520.         }

  1522.         /** Get the slot valid for some date.
  1523.          * @param meanStates mean states defining the slot
  1524.          * @return slot valid at the specified date
  1525.          */
  1526.         public Slot createSlot(final SpacecraftState... meanStates) {
  1527.             final Slot         slot  = new Slot(maxFrequencyShortPeriodics, interpolationPoints);
  1528.             final AbsoluteDate first = meanStates[0].getDate();
  1529.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  1530.             final int compare = first.compareTo(last);
  1531.             if (compare < 0) {
  1532.                 slots.addValidAfter(slot, first, false);
  1533.             } else if (compare > 0) {
  1534.                 slots.addValidBefore(slot, first, false);
  1535.             } else {
  1536.                 // single date, valid for all time
  1537.                 slots.addValidAfter(slot, AbsoluteDate.PAST_INFINITY, false);
  1538.             }
  1539.             return slot;
  1540.         }

  1542.         /** {@inheritDoc} */
  1543.         @Override
  1544.         public double[] value(final Orbit meanOrbit) {

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

  1549.             // Get the True longitude L
  1550.             final double L = meanOrbit.getLv();

  1552.             // Compute the center
  1553.             final double center = L - meanOrbit.getLM();

  1555.             // Initialize short periodic variations
  1556.             final double[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  1557.             final double[] d = slot.di.value(meanOrbit.getDate());
  1558.             for (int i = 0; i < 6; i++) {
  1559.                 shortPeriodicVariation[i] +=  center * d[i];
  1560.             }

  1562.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1563.                 final double[] c = slot.cij[j].value(meanOrbit.getDate());
  1564.                 final double[] s = slot.sij[j].value(meanOrbit.getDate());
  1565.                 final SinCos sc  = FastMath.sinCos(j * L);
  1566.                 for (int i = 0; i < 6; i++) {
  1567.                     // add corresponding term to the short periodic variation
  1568.                     shortPeriodicVariation[i] += c[i] * sc.cos();
  1569.                     shortPeriodicVariation[i] += s[i] * sc.sin();
  1570.                 }
  1571.             }

  1573.             return shortPeriodicVariation;
  1574.         }

  1576.         /** {@inheritDoc} */
  1577.         @Override
  1578.         public String getCoefficientsKeyPrefix() {
  1579.             return DSSTZonal.SHORT_PERIOD_PREFIX;
  1580.         }

  1582.         /** {@inheritDoc}
  1583.          * <p>
  1584.          * For zonal terms contributions,there are maxJ cj coefficients,
  1585.          * maxJ sj coefficients and 2 dj coefficients, where maxJ depends
  1586.          * on the orbit. The j index is the integer multiplier for the true
  1587.          * longitude argument in the cj and sj coefficients and the degree
  1588.          * in the polynomial dj coefficients.
  1589.          * </p>
  1590.          */
  1591.         @Override
  1592.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected) {

  1594.             // select the coefficients slot
  1595.             final Slot slot = slots.get(date);

  1597.             final Map<String, double[]> coefficients = new HashMap<>(2 * maxFrequencyShortPeriodics + 2);
  1598.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  1599.             storeIfSelected(coefficients, selected, slot.di.value(date), "d", 1);
  1600.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1601.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  1602.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  1603.             }
  1604.             return coefficients;

  1606.         }

  1608.         /** Put a coefficient in a map if selected.
  1609.          * @param map map to populate
  1610.          * @param selected set of coefficients that should be put in the map
  1611.          * (empty set means all coefficients are selected)
  1612.          * @param value coefficient value
  1613.          * @param id coefficient identifier
  1614.          * @param indices list of coefficient indices
  1615.          */
  1616.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
  1617.                                      final double[] value, final String id, final int... indices) {
  1618.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  1619.             keyBuilder.append(id);
  1620.             for (int index : indices) {
  1621.                 keyBuilder.append('[').append(index).append(']');
  1622.             }
  1623.             final String key = keyBuilder.toString();
  1624.             if (selected.isEmpty() || selected.contains(key)) {
  1625.                 map.put(key, value);
  1626.             }
  1627.         }

  1629.     }

  1631.     /** The coefficients used to compute the short-periodic zonal contribution.
  1632.     *
  1633.     * <p>
  1634.     * Those coefficients are given in Danielson paper by expressions 4.1-(20) to 4.1.-(25)
  1635.     * </p>
  1636.     * <p>
  1637.     * The coefficients are: <br>
  1638.     * - C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> <br>
  1639.     * - ρ<sub>j</sub> and σ<sub>j</sub> <br>
  1640.     * - C<sub>i</sub>⁰
  1641.     * </p>
  1642.     *
  1643.     * @author Lucian Barbulescu
  1644.     */
  1645.     private static class FieldZonalShortPeriodicCoefficients <T extends CalculusFieldElement<T>> implements FieldShortPeriodTerms<T> {

  1647.         /** Maximum value for j index. */
  1648.         private final int maxFrequencyShortPeriodics;

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

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

  1656.        /** Constructor.
  1657.         * @param maxFrequencyShortPeriodics maximum value for j index
  1658.         * @param interpolationPoints number of points used in the interpolation process
  1659.         * @param slots all coefficients slots
  1660.         */
  1661.         FieldZonalShortPeriodicCoefficients(final int maxFrequencyShortPeriodics, final int interpolationPoints,
  1662.                                             final FieldTimeSpanMap<FieldSlot<T>, T> slots) {

  1664.             // Save parameters
  1665.             this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
  1666.             this.interpolationPoints        = interpolationPoints;
  1667.             this.slots                      = slots;

  1669.         }

  1671.        /** Get the slot valid for some date.
  1672.         * @param meanStates mean states defining the slot
  1673.         * @return slot valid at the specified date
  1674.         */
  1675.         @SuppressWarnings("unchecked")
  1676.         public FieldSlot<T> createSlot(final FieldSpacecraftState<T>... meanStates) {
  1677.             final FieldSlot<T>         slot  = new FieldSlot<>(maxFrequencyShortPeriodics, interpolationPoints);
  1678.             final FieldAbsoluteDate<T> first = meanStates[0].getDate();
  1679.             final FieldAbsoluteDate<T> last  = meanStates[meanStates.length - 1].getDate();
  1680.             if (first.compareTo(last) <= 0) {
  1681.                 slots.addValidAfter(slot, first);
  1682.             } else {
  1683.                 slots.addValidBefore(slot, first);
  1684.             }
  1685.             return slot;
  1686.         }

  1688.         /** {@inheritDoc} */
  1689.         @Override
  1690.         public T[] value(final FieldOrbit<T> meanOrbit) {

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

  1695.             // Get the True longitude L
  1696.             final T L = meanOrbit.getLv();

  1698.             // Compute the center
  1699.             final T center = L.subtract(meanOrbit.getLM());

  1701.             // Initialize short periodic variations
  1702.             final T[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  1703.             final T[] d = slot.di.value(meanOrbit.getDate());
  1704.             for (int i = 0; i < 6; i++) {
  1705.                 shortPeriodicVariation[i] = shortPeriodicVariation[i].add(center.multiply(d[i]));
  1706.             }

  1708.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1709.                 final T[]            c   = slot.cij[j].value(meanOrbit.getDate());
  1710.                 final T[]            s   = slot.sij[j].value(meanOrbit.getDate());
  1711.                 final FieldSinCos<T> sc  = FastMath.sinCos(L.multiply(j));
  1712.                 for (int i = 0; i < 6; i++) {
  1713.                     // add corresponding term to the short periodic variation
  1714.                     shortPeriodicVariation[i] = shortPeriodicVariation[i].add(c[i].multiply(sc.cos()));
  1715.                     shortPeriodicVariation[i] = shortPeriodicVariation[i].add(s[i].multiply(sc.sin()));
  1716.                 }
  1717.             }

  1719.             return shortPeriodicVariation;
  1720.         }

  1722.         /** {@inheritDoc} */
  1723.         @Override
  1724.         public String getCoefficientsKeyPrefix() {
  1725.             return DSSTZonal.SHORT_PERIOD_PREFIX;
  1726.         }

  1728.        /** {@inheritDoc}
  1729.         * <p>
  1730.         * For zonal terms contributions,there are maxJ cj coefficients,
  1731.         * maxJ sj coefficients and 2 dj coefficients, where maxJ depends
  1732.         * on the orbit. The j index is the integer multiplier for the true
  1733.         * longitude argument in the cj and sj coefficients and the degree
  1734.         * in the polynomial dj coefficients.
  1735.         * </p>
  1736.         */
  1737.         @Override
  1738.         public Map<String, T[]> getCoefficients(final FieldAbsoluteDate<T> date, final Set<String> selected) {

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

  1743.             final Map<String, T[]> coefficients = new HashMap<>(2 * maxFrequencyShortPeriodics + 2);
  1744.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  1745.             storeIfSelected(coefficients, selected, slot.di.value(date), "d", 1);
  1746.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1747.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  1748.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  1749.             }
  1750.             return coefficients;

  1752.         }

  1754.        /** Put a coefficient in a map if selected.
  1755.         * @param map map to populate
  1756.         * @param selected set of coefficients that should be put in the map
  1757.         * (empty set means all coefficients are selected)
  1758.         * @param value coefficient value
  1759.         * @param id coefficient identifier
  1760.         * @param indices list of coefficient indices
  1761.         */
  1762.         private void storeIfSelected(final Map<String, T[]> map, final Set<String> selected,
  1763.                                      final T[] value, final String id, final int... indices) {
  1764.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  1765.             keyBuilder.append(id);
  1766.             for (int index : indices) {
  1767.                 keyBuilder.append('[').append(index).append(']');
  1768.             }
  1769.             final String key = keyBuilder.toString();
  1770.             if (selected.isEmpty() || selected.contains(key)) {
  1771.                 map.put(key, value);
  1772.             }
  1773.         }

  1775.     }

  1777.     /** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
  1778.      *  <p>
  1779.      *  Those coefficients are given in Danielson paper by expressions 4.1-(13) to 4.1.-(16b)
  1780.      *  </p>
  1781.      */
  1782.     private class FourierCjSjCoefficients {

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

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

  1790.         /** Maximum possible value for n. */
  1791.         private final int nMax;

  1793.         /** Maximum possible value for s. */
  1794.         private final int sMax;

  1796.         /** Maximum possible value for j. */
  1797.         private final int jMax;

  1799.         /** The C<sup>j</sup> coefficients and their derivatives.
  1800.          * <p>
  1801.          * Each column of the matrix contains the following values: <br/>
  1802.          * - C<sup>j</sup> <br/>
  1803.          * - dC<sup>j</sup> / da <br/>
  1804.          * - dC<sup>j</sup> / dh <br/>
  1805.          * - dC<sup>j</sup> / dk <br/>
  1806.          * - dC<sup>j</sup> / dα <br/>
  1807.          * - dC<sup>j</sup> / dβ <br/>
  1808.          * - dC<sup>j</sup> / dγ <br/>
  1809.          * </p>
  1810.          */
  1811.         private final double[][] cCoef;

  1813.         /** The S<sup>j</sup> coefficients and their derivatives.
  1814.          * <p>
  1815.          * Each column of the matrix contains the following values: <br/>
  1816.          * - S<sup>j</sup> <br/>
  1817.          * - dS<sup>j</sup> / da <br/>
  1818.          * - dS<sup>j</sup> / dh <br/>
  1819.          * - dS<sup>j</sup> / dk <br/>
  1820.          * - dS<sup>j</sup> / dα <br/>
  1821.          * - dS<sup>j</sup> / dβ <br/>
  1822.          * - dS<sup>j</sup> / dγ <br/>
  1823.          * </p>
  1824.          */
  1825.         private final double[][] sCoef;

  1827.         /** h * &Chi;³. */
  1828.         private final double hXXX;
  1829.         /** k * &Chi;³. */
  1830.         private final double kXXX;

  1832.         /** Create a set of C<sup>j</sup> and the S<sup>j</sup> coefficients.
  1833.          *  @param date the current date
  1834.          *  @param nMax maximum possible value for n
  1835.          *  @param sMax maximum possible value for s
  1836.          *  @param jMax maximum possible value for j
  1837.          *  @param context container for attributes
  1838.          *  @param hansenObjects initialization of hansen objects
  1839.          */
  1840.         FourierCjSjCoefficients(final AbsoluteDate date,
  1841.                                 final int nMax, final int sMax, final int jMax, final DSSTZonalContext context,
  1842.                                 final HansenObjects hansenObjects) {

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

  1846.             this.ghijCoef = new GHIJjsPolynomials(auxiliaryElements.getK(), auxiliaryElements.getH(), context.getAlpha(), context.getBeta());
  1847.             // Qns coefficients
  1848.             final double[][] Qns  = CoefficientsFactory.computeQns(context.getGamma(), nMax, nMax);

  1850.             this.lnsCoef = new LnsCoefficients(nMax, nMax, Qns, Vns, context.getRoa());
  1851.             this.nMax = nMax;
  1852.             this.sMax = sMax;
  1853.             this.jMax = jMax;

  1855.             // compute the common factors that depends on the mean elements
  1856.             this.hXXX = auxiliaryElements.getH() * context.getChi3();
  1857.             this.kXXX = auxiliaryElements.getK() * context.getChi3();

  1859.             this.cCoef = new double[7][jMax + 1];
  1860.             this.sCoef = new double[7][jMax + 1];

  1862.             for (int s = 0; s <= sMax; s++) {
  1863.                 //Initialise the Hansen roots
  1864.                 hansenObjects.computeHansenObjectsInitValues(context, s);
  1865.             }
  1866.             generateCoefficients(date, context, auxiliaryElements, hansenObjects);
  1867.         }

  1869.         /** Generate all coefficients.
  1870.          * @param date the current date
  1871.          * @param context container for attributes
  1872.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1873.          * @param hansenObjects initialization of hansen objects
  1874.          */
  1875.         private void generateCoefficients(final AbsoluteDate date,
  1876.                                           final DSSTZonalContext context,
  1877.                                           final AuxiliaryElements auxiliaryElements,
  1878.                                           final HansenObjects hansenObjects) {

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

  1883.                 //init arrays
  1884.                 for (int i = 0; i <= 6; i++) {
  1885.                     cCoef[i][j] = 0.;
  1886.                     sCoef[i][j] = 0.;
  1887.                 }

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

  1891.                     //compute first double sum where s: j -> N-1 and n: s+1 -> N
  1892.                     for (int s = j; s <= FastMath.min(nMax - 1, sMax); s++) {
  1893.                         // j - s
  1894.                         final int jms = j - s;
  1895.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  1896.                         final int d0smj = (s == j) ? 1 : 2;

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

  1905.                                 final double hjs = ghijCoef.getHjs(s, -jms);
  1906.                                 final double dHjsdh = ghijCoef.getdHjsdh(s, -jms);
  1907.                                 final double dHjsdk = ghijCoef.getdHjsdk(s, -jms);
  1908.                                 final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, -jms);
  1909.                                 final double dHjsdBeta = ghijCoef.getdHjsdBeta(s, -jms);

  1911.                                 final double gjs = ghijCoef.getGjs(s, -jms);
  1912.                                 final double dGjsdh = ghijCoef.getdGjsdh(s, -jms);
  1913.                                 final double dGjsdk = ghijCoef.getdGjsdk(s, -jms);
  1914.                                 final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, -jms);
  1915.                                 final double dGjsdBeta = ghijCoef.getdGjsdBeta(s, -jms);

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

  1920.                                 // K₀<sup>-n-1,s</sup>
  1921.                                 final double kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  1922.                                 final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  1924.                                 final double coef0 = d0smj * jn;
  1925.                                 final double coef1 = coef0 * lns;
  1926.                                 final double coef2 = coef1 * kns;
  1927.                                 final double coef3 = coef2 * hjs;
  1928.                                 final double coef4 = coef2 * gjs;

  1930.                                 // Add the term to the coefficients
  1931.                                 cCoef[0][j] += coef3;
  1932.                                 cCoef[1][j] += coef3 * (n + 1);
  1933.                                 cCoef[2][j] += coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
  1934.                                 cCoef[3][j] += coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
  1935.                                 cCoef[4][j] += coef2 * dHjsdAlpha;
  1936.                                 cCoef[5][j] += coef2 * dHjsdBeta;
  1937.                                 cCoef[6][j] += coef0 * dlns * kns * hjs;

  1939.                                 sCoef[0][j] += coef4;
  1940.                                 sCoef[1][j] += coef4 * (n + 1);
  1941.                                 sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
  1942.                                 sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
  1943.                                 sCoef[4][j] += coef2 * dGjsdAlpha;
  1944.                                 sCoef[5][j] += coef2 * dGjsdBeta;
  1945.                                 sCoef[6][j] += coef0 * dlns * kns * gjs;
  1946.                             }
  1947.                         }
  1948.                     }

  1950.                     //compute second double sum where s: 0 -> N-j and n: max(j+s, j+1) -> N
  1951.                     for (int s = 0; s <= FastMath.min(nMax - j, sMax); s++) {
  1952.                         // j + s
  1953.                         final int jps = j + s;
  1954.                         // Kronecker symbols (2 - delta(0,j+s))
  1955.                         final double d0spj = (s == -j) ? 1 : 2;

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

  1964.                                 final double hjs = ghijCoef.getHjs(s, jps);
  1965.                                 final double dHjsdh = ghijCoef.getdHjsdh(s, jps);
  1966.                                 final double dHjsdk = ghijCoef.getdHjsdk(s, jps);
  1967.                                 final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, jps);
  1968.                                 final double dHjsdBeta = ghijCoef.getdHjsdBeta(s, jps);

  1970.                                 final double gjs = ghijCoef.getGjs(s, jps);
  1971.                                 final double dGjsdh = ghijCoef.getdGjsdh(s, jps);
  1972.                                 final double dGjsdk = ghijCoef.getdGjsdk(s, jps);
  1973.                                 final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, jps);
  1974.                                 final double dGjsdBeta = ghijCoef.getdGjsdBeta(s, jps);

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

  1979.                                 // K₀<sup>-n-1,s</sup>
  1980.                                 final double kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  1981.                                 final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  1983.                                 final double coef0 = d0spj * jn;
  1984.                                 final double coef1 = coef0 * lns;
  1985.                                 final double coef2 = coef1 * kns;

  1987.                                 final double coef3 = coef2 * hjs;
  1988.                                 final double coef4 = coef2 * gjs;

  1990.                                 // Add the term to the coefficients
  1991.                                 cCoef[0][j] -= coef3;
  1992.                                 cCoef[1][j] -= coef3 * (n + 1);
  1993.                                 cCoef[2][j] -= coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
  1994.                                 cCoef[3][j] -= coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
  1995.                                 cCoef[4][j] -= coef2 * dHjsdAlpha;
  1996.                                 cCoef[5][j] -= coef2 * dHjsdBeta;
  1997.                                 cCoef[6][j] -= coef0 * dlns * kns * hjs;

  1999.                                 sCoef[0][j] += coef4;
  2000.                                 sCoef[1][j] += coef4 * (n + 1);
  2001.                                 sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
  2002.                                 sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
  2003.                                 sCoef[4][j] += coef2 * dGjsdAlpha;
  2004.                                 sCoef[5][j] += coef2 * dGjsdBeta;
  2005.                                 sCoef[6][j] += coef0 * dlns * kns * gjs;
  2006.                             }
  2007.                         }
  2008.                     }

  2010.                     //compute third double sum where s: 1 -> j and  n: j+1 -> N
  2011.                     for (int s = 1; s <= FastMath.min(j, sMax); s++) {
  2012.                         // j - s
  2013.                         final int jms = j - s;
  2014.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2015.                         final int d0smj = (s == j) ? 1 : 2;

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

  2024.                                 final double ijs = ghijCoef.getIjs(s, jms);
  2025.                                 final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2026.                                 final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2027.                                 final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2028.                                 final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2030.                                 final double jjs = ghijCoef.getJjs(s, jms);
  2031.                                 final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2032.                                 final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2033.                                 final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2034.                                 final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

  2039.                                 // K₀<sup>-n-1,s</sup>
  2040.                                 final double kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  2041.                                 final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  2043.                                 final double coef0 = d0smj * jn;
  2044.                                 final double coef1 = coef0 * lns;
  2045.                                 final double coef2 = coef1 * kns;

  2047.                                 final double coef3 = coef2 * ijs;
  2048.                                 final double coef4 = coef2 * jjs;

  2050.                                 // Add the term to the coefficients
  2051.                                 cCoef[0][j] -= coef3;
  2052.                                 cCoef[1][j] -= coef3 * (n + 1);
  2053.                                 cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2054.                                 cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2055.                                 cCoef[4][j] -= coef2 * dIjsdAlpha;
  2056.                                 cCoef[5][j] -= coef2 * dIjsdBeta;
  2057.                                 cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2059.                                 sCoef[0][j] += coef4;
  2060.                                 sCoef[1][j] += coef4 * (n + 1);
  2061.                                 sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2062.                                 sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2063.                                 sCoef[4][j] += coef2 * dJjsdAlpha;
  2064.                                 sCoef[5][j] += coef2 * dJjsdBeta;
  2065.                                 sCoef[6][j] += coef0 * dlns * kns * jjs;
  2066.                             }
  2067.                         }
  2068.                     }
  2069.                 }

  2071.                 if (isBetween(j, 2, nMax)) {
  2072.                     //add first term
  2073.                     // J<sub>j</sub>
  2074.                     final double jj = -harmonics.getUnnormalizedCnm(j, 0);
  2075.                     double kns = hansenObjects.getHansenObjects()[0].getValue(-j - 1, context.getChi());
  2076.                     double dkns = hansenObjects.getHansenObjects()[0].getDerivative(-j - 1, context.getChi());

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

  2081.                     final double hjs = ghijCoef.getHjs(0, j);
  2082.                     final double dHjsdh = ghijCoef.getdHjsdh(0, j);
  2083.                     final double dHjsdk = ghijCoef.getdHjsdk(0, j);
  2084.                     final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(0, j);
  2085.                     final double dHjsdBeta = ghijCoef.getdHjsdBeta(0, j);

  2087.                     final double gjs = ghijCoef.getGjs(0, j);
  2088.                     final double dGjsdh = ghijCoef.getdGjsdh(0, j);
  2089.                     final double dGjsdk = ghijCoef.getdGjsdk(0, j);
  2090.                     final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(0, j);
  2091.                     final double dGjsdBeta = ghijCoef.getdGjsdBeta(0, j);

  2093.                     // 2 * J<sub>j</sub> * K₀<sup>-j-1,0</sup> * L<sub>j</sub><sup>j</sup>
  2094.                     double coef0 = 2 * jj;
  2095.                     double coef1 = coef0 * lns;
  2096.                     double coef2 = coef1 * kns;

  2098.                     double coef3 = coef2 * hjs;
  2099.                     double coef4 = coef2 * gjs;

  2101.                     // Add the term to the coefficients
  2102.                     cCoef[0][j] -= coef3;
  2103.                     cCoef[1][j] -= coef3 * (j + 1);
  2104.                     cCoef[2][j] -= coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
  2105.                     cCoef[3][j] -= coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
  2106.                     cCoef[4][j] -= coef2 * dHjsdAlpha;
  2107.                     cCoef[5][j] -= coef2 * dHjsdBeta;
  2108.                     //no contribution to cCoef[6][j] because dlns is 0

  2110.                     sCoef[0][j] += coef4;
  2111.                     sCoef[1][j] += coef4 * (j + 1);
  2112.                     sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
  2113.                     sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
  2114.                     sCoef[4][j] += coef2 * dGjsdAlpha;
  2115.                     sCoef[5][j] += coef2 * dGjsdBeta;
  2116.                     //no contribution to sCoef[6][j] because dlns is 0

  2118.                     //compute simple sum where s: 1 -> j-1
  2119.                     for (int s = 1; s <= FastMath.min(j - 1, sMax); s++) {
  2120.                         // j - s
  2121.                         final int jms = j - s;
  2122.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2123.                         final int d0smj = (s == j) ? 1 : 2;

  2125.                         // if s is odd, then the term is equal to zero due to the factor Vj,s-j
  2126.                         if (s % 2 == 0) {
  2127.                             // (2 - delta(0,j-s)) * J<sub>j</sub> * K₀<sup>-j-1,s</sup> * L<sub>j</sub><sup>j-s</sup>
  2128.                             kns   = hansenObjects.getHansenObjects()[s].getValue(-j - 1, context.getChi());
  2129.                             dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-j - 1, context.getChi());

  2131.                             lns = lnsCoef.getLns(j, jms);
  2132.                             final double dlns = lnsCoef.getdLnsdGamma(j, jms);

  2134.                             final double ijs = ghijCoef.getIjs(s, jms);
  2135.                             final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2136.                             final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2137.                             final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2138.                             final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2140.                             final double jjs = ghijCoef.getJjs(s, jms);
  2141.                             final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2142.                             final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2143.                             final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2144.                             final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

  2146.                             coef0 = d0smj * jj;
  2147.                             coef1 = coef0 * lns;
  2148.                             coef2 = coef1 * kns;

  2150.                             coef3 = coef2 * ijs;
  2151.                             coef4 = coef2 * jjs;

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

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

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

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

  2179.                     //if j is even
  2180.                     if (j % 2 == 0) {
  2181.                         //compute first double sum where s: j-min(j-1,N) -> j/2-1 and n: j-s -> min(j-1,N)
  2182.                         for (int s = j - minjm1on; s <= FastMath.min(j / 2 - 1, sMax); s++) {
  2183.                             // j - s
  2184.                             final int jms = j - s;
  2185.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2186.                             final int d0smj = (s == j) ? 1 : 2;

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

  2195.                                     final double ijs = ghijCoef.getIjs(s, jms);
  2196.                                     final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2197.                                     final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2198.                                     final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2199.                                     final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2201.                                     final double jjs = ghijCoef.getJjs(s, jms);
  2202.                                     final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2203.                                     final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2204.                                     final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2205.                                     final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

  2210.                                     // K₀<sup>-n-1,s</sup>
  2211.                                     final double kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  2212.                                     final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  2214.                                     final double coef0 = d0smj * jn;
  2215.                                     final double coef1 = coef0 * lns;
  2216.                                     final double coef2 = coef1 * kns;

  2218.                                     final double coef3 = coef2 * ijs;
  2219.                                     final double coef4 = coef2 * jjs;

  2221.                                     // Add the term to the coefficients
  2222.                                     cCoef[0][j] -= coef3;
  2223.                                     cCoef[1][j] -= coef3 * (n + 1);
  2224.                                     cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2225.                                     cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2226.                                     cCoef[4][j] -= coef2 * dIjsdAlpha;
  2227.                                     cCoef[5][j] -= coef2 * dIjsdBeta;
  2228.                                     cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2230.                                     sCoef[0][j] += coef4;
  2231.                                     sCoef[1][j] += coef4 * (n + 1);
  2232.                                     sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2233.                                     sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2234.                                     sCoef[4][j] += coef2 * dJjsdAlpha;
  2235.                                     sCoef[5][j] += coef2 * dJjsdBeta;
  2236.                                     sCoef[6][j] += coef0 * dlns * kns * jjs;
  2237.                                 }
  2238.                             }
  2239.                         }

  2241.                         //compute second double sum where s: j/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  2242.                         for (int s = j / 2; s <=  FastMath.min(minjm1on - 1, sMax); s++) {
  2243.                             // j - s
  2244.                             final int jms = j - s;
  2245.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2246.                             final int d0smj = (s == j) ? 1 : 2;

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

  2255.                                     final double ijs = ghijCoef.getIjs(s, jms);
  2256.                                     final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2257.                                     final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2258.                                     final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2259.                                     final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2261.                                     final double jjs = ghijCoef.getJjs(s, jms);
  2262.                                     final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2263.                                     final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2264.                                     final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2265.                                     final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

  2270.                                     // K₀<sup>-n-1,s</sup>
  2271.                                     final double kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  2272.                                     final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  2274.                                     final double coef0 = d0smj * jn;
  2275.                                     final double coef1 = coef0 * lns;
  2276.                                     final double coef2 = coef1 * kns;

  2278.                                     final double coef3 = coef2 * ijs;
  2279.                                     final double coef4 = coef2 * jjs;

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

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

  2302.                     //if j is odd
  2303.                     else {
  2304.                         //compute first double sum where s: (j-1)/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  2305.                         for (int s = (j - 1) / 2; s <= FastMath.min(minjm1on - 1, sMax); s++) {
  2306.                             // j - s
  2307.                             final int jms = j - s;
  2308.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2309.                             final int d0smj = (s == j) ? 1 : 2;

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

  2318.                                     final double ijs = ghijCoef.getIjs(s, jms);
  2319.                                     final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2320.                                     final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2321.                                     final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2322.                                     final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2324.                                     final double jjs = ghijCoef.getJjs(s, jms);
  2325.                                     final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2326.                                     final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2327.                                     final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2328.                                     final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

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

  2335.                                     final double kns = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  2336.                                     final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  2338.                                     final double coef0 = d0smj * jn;
  2339.                                     final double coef1 = coef0 * lns;
  2340.                                     final double coef2 = coef1 * kns;

  2342.                                     final double coef3 = coef2 * ijs;
  2343.                                     final double coef4 = coef2 * jjs;

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

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

  2365.                         //the second double sum is added only if N >= 4 and j between 5 and 2*N-3
  2366.                         if (nMax >= 4 && isBetween(j, 5, 2 * nMax - 3)) {
  2367.                             //compute second double sum where s: j-min(j-1,N) -> (j-3)/2 and n: j-s -> min(j-1,N)
  2368.                             for (int s = j - minjm1on; s <= FastMath.min((j - 3) / 2, sMax); s++) {
  2369.                                 // j - s
  2370.                                 final int jms = j - s;
  2371.                                 // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2372.                                 final int d0smj = (s == j) ? 1 : 2;

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

  2381.                                         final double ijs = ghijCoef.getIjs(s, jms);
  2382.                                         final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  2383.                                         final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  2384.                                         final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2385.                                         final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  2387.                                         final double jjs = ghijCoef.getJjs(s, jms);
  2388.                                         final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  2389.                                         final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  2390.                                         final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2391.                                         final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

  2396.                                         // K₀<sup>-n-1,s</sup>
  2397.                                         final double kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  2398.                                         final double dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  2400.                                         final double coef0 = d0smj * jn;
  2401.                                         final double coef1 = coef0 * lns;
  2402.                                         final double coef2 = coef1 * kns;

  2404.                                         final double coef3 = coef2 * ijs;
  2405.                                         final double coef4 = coef2 * jjs;

  2407.                                         // Add the term to the coefficients
  2408.                                         cCoef[0][j] -= coef3;
  2409.                                         cCoef[1][j] -= coef3 * (n + 1);
  2410.                                         cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  2411.                                         cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  2412.                                         cCoef[4][j] -= coef2 * dIjsdAlpha;
  2413.                                         cCoef[5][j] -= coef2 * dIjsdBeta;
  2414.                                         cCoef[6][j] -= coef0 * dlns * kns * ijs;

  2416.                                         sCoef[0][j] += coef4;
  2417.                                         sCoef[1][j] += coef4 * (n + 1);
  2418.                                         sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  2419.                                         sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  2420.                                         sCoef[4][j] += coef2 * dJjsdAlpha;
  2421.                                         sCoef[5][j] += coef2 * dJjsdBeta;
  2422.                                         sCoef[6][j] += coef0 * dlns * kns * jjs;
  2423.                                     }
  2424.                                 }
  2425.                             }
  2426.                         }
  2427.                     }
  2428.                 }

  2430.                 cCoef[0][j] *= -context.getMuoa() / j;
  2431.                 cCoef[1][j] *=  context.getMuoa() / ( j * auxiliaryElements.getSma() );
  2432.                 cCoef[2][j] *= -context.getMuoa() / j;
  2433.                 cCoef[3][j] *= -context.getMuoa() / j;
  2434.                 cCoef[4][j] *= -context.getMuoa() / j;
  2435.                 cCoef[5][j] *= -context.getMuoa() / j;
  2436.                 cCoef[6][j] *= -context.getMuoa() / j;

  2438.                 sCoef[0][j] *= -context.getMuoa() / j;
  2439.                 sCoef[1][j] *=  context.getMuoa() / ( j * auxiliaryElements.getSma() );
  2440.                 sCoef[2][j] *= -context.getMuoa() / j;
  2441.                 sCoef[3][j] *= -context.getMuoa() / j;
  2442.                 sCoef[4][j] *= -context.getMuoa() / j;
  2443.                 sCoef[5][j] *= -context.getMuoa() / j;
  2444.                 sCoef[6][j] *= -context.getMuoa() / j;

  2446.             }
  2447.         }

  2449.         /** Check if an index is within the accepted interval.
  2450.          *
  2451.          * @param index the index to check
  2452.          * @param lowerBound the lower bound of the interval
  2453.          * @param upperBound the upper bound of the interval
  2454.          * @return true if the index is between the lower and upper bounds, false otherwise
  2455.          */
  2456.         private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
  2457.             return index >= lowerBound && index <= upperBound;
  2458.         }

  2460.         /**Get the value of C<sup>j</sup>.
  2461.          *
  2462.          * @param j j index
  2463.          * @return C<sup>j</sup>
  2464.          */
  2465.         public double getCj(final int j) {
  2466.             return cCoef[0][j];
  2467.         }

  2469.         /**Get the value of dC<sup>j</sup> / da.
  2470.          *
  2471.          * @param j j index
  2472.          * @return dC<sup>j</sup> / da
  2473.          */
  2474.         public double getdCjdA(final int j) {
  2475.             return cCoef[1][j];
  2476.         }

  2478.         /**Get the value of dC<sup>j</sup> / dh.
  2479.          *
  2480.          * @param j j index
  2481.          * @return dC<sup>j</sup> / dh
  2482.          */
  2483.         public double getdCjdH(final int j) {
  2484.             return cCoef[2][j];
  2485.         }

  2487.         /**Get the value of dC<sup>j</sup> / dk.
  2488.          *
  2489.          * @param j j index
  2490.          * @return dC<sup>j</sup> / dk
  2491.          */
  2492.         public double getdCjdK(final int j) {
  2493.             return cCoef[3][j];
  2494.         }

  2496.         /**Get the value of dC<sup>j</sup> / dα.
  2497.          *
  2498.          * @param j j index
  2499.          * @return dC<sup>j</sup> / dα
  2500.          */
  2501.         public double getdCjdAlpha(final int j) {
  2502.             return cCoef[4][j];
  2503.         }

  2505.         /**Get the value of dC<sup>j</sup> / dβ.
  2506.          *
  2507.          * @param j j index
  2508.          * @return dC<sup>j</sup> / dβ
  2509.          */
  2510.         public double getdCjdBeta(final int j) {
  2511.             return cCoef[5][j];
  2512.         }

  2514.         /**Get the value of dC<sup>j</sup> / dγ.
  2515.          *
  2516.          * @param j j index
  2517.          * @return dC<sup>j</sup> / dγ
  2518.          */
  2519.         public double getdCjdGamma(final int j) {
  2520.             return cCoef[6][j];
  2521.         }

  2523.         /**Get the value of S<sup>j</sup>.
  2524.          *
  2525.          * @param j j index
  2526.          * @return S<sup>j</sup>
  2527.          */
  2528.         public double getSj(final int j) {
  2529.             return sCoef[0][j];
  2530.         }

  2532.         /**Get the value of dS<sup>j</sup> / da.
  2533.          *
  2534.          * @param j j index
  2535.          * @return dS<sup>j</sup> / da
  2536.          */
  2537.         public double getdSjdA(final int j) {
  2538.             return sCoef[1][j];
  2539.         }

  2541.         /**Get the value of dS<sup>j</sup> / dh.
  2542.          *
  2543.          * @param j j index
  2544.          * @return dS<sup>j</sup> / dh
  2545.          */
  2546.         public double getdSjdH(final int j) {
  2547.             return sCoef[2][j];
  2548.         }

  2550.         /**Get the value of dS<sup>j</sup> / dk.
  2551.          *
  2552.          * @param j j index
  2553.          * @return dS<sup>j</sup> / dk
  2554.          */
  2555.         public double getdSjdK(final int j) {
  2556.             return sCoef[3][j];
  2557.         }

  2559.         /**Get the value of dS<sup>j</sup> / dα.
  2560.          *
  2561.          * @param j j index
  2562.          * @return dS<sup>j</sup> / dα
  2563.          */
  2564.         public double getdSjdAlpha(final int j) {
  2565.             return sCoef[4][j];
  2566.         }

  2568.         /**Get the value of dS<sup>j</sup> / dβ.
  2569.          *
  2570.          * @param j j index
  2571.          * @return dS<sup>j</sup> / dβ
  2572.          */
  2573.         public double getdSjdBeta(final int j) {
  2574.             return sCoef[5][j];
  2575.         }

  2577.         /**Get the value of dS<sup>j</sup> /  dγ.
  2578.          *
  2579.          * @param j j index
  2580.          * @return dS<sup>j</sup> /  dγ
  2581.          */
  2582.         public double getdSjdGamma(final int j) {
  2583.             return sCoef[6][j];
  2584.         }
  2585.     }

  2587.     /** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
  2588.      *  <p>
  2589.      *  Those coefficients are given in Danielson paper by expressions 4.1-(13) to 4.1.-(16b)
  2590.      *  </p>
  2591.      */
  2592.     private class FieldFourierCjSjCoefficients <T extends CalculusFieldElement<T>> {

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

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

  2600.         /** Maximum possible value for n. */
  2601.         private final int nMax;

  2603.         /** Maximum possible value for s. */
  2604.         private final int sMax;

  2606.         /** Maximum possible value for j. */
  2607.         private final int jMax;

  2609.         /** The C<sup>j</sup> coefficients and their derivatives.
  2610.          * <p>
  2611.          * Each column of the matrix contains the following values: <br/>
  2612.          * - C<sup>j</sup> <br/>
  2613.          * - dC<sup>j</sup> / da <br/>
  2614.          * - dC<sup>j</sup> / dh <br/>
  2615.          * - dC<sup>j</sup> / dk <br/>
  2616.          * - dC<sup>j</sup> / dα <br/>
  2617.          * - dC<sup>j</sup> / dβ <br/>
  2618.          * - dC<sup>j</sup> / dγ <br/>
  2619.          * </p>
  2620.          */
  2621.         private final T[][] cCoef;

  2623.         /** The S<sup>j</sup> coefficients and their derivatives.
  2624.          * <p>
  2625.          * Each column of the matrix contains the following values: <br/>
  2626.          * - S<sup>j</sup> <br/>
  2627.          * - dS<sup>j</sup> / da <br/>
  2628.          * - dS<sup>j</sup> / dh <br/>
  2629.          * - dS<sup>j</sup> / dk <br/>
  2630.          * - dS<sup>j</sup> / dα <br/>
  2631.          * - dS<sup>j</sup> / dβ <br/>
  2632.          * - dS<sup>j</sup> / dγ <br/>
  2633.          * </p>
  2634.          */
  2635.         private final T[][] sCoef;

  2637.         /** h * &Chi;³. */
  2638.         private final T hXXX;
  2639.         /** k * &Chi;³. */
  2640.         private final T kXXX;

  2642.         /** Create a set of C<sup>j</sup> and the S<sup>j</sup> coefficients.
  2643.          *  @param date the current date
  2644.          *  @param nMax maximum possible value for n
  2645.          *  @param sMax maximum possible value for s
  2646.          *  @param jMax maximum possible value for j
  2647.          *  @param context container for attributes
  2648.          *  @param hansenObjects initialization of hansen objects
  2649.          */
  2650.         FieldFourierCjSjCoefficients(final FieldAbsoluteDate<T> date,
  2651.                                      final int nMax, final int sMax, final int jMax,
  2652.                                      final FieldDSSTZonalContext<T> context,
  2653.                                      final FieldHansenObjects<T> hansenObjects) {

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

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

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

  2664.             this.lnsCoef = new FieldLnsCoefficients<>(nMax, nMax, Qns, Vns, context.getRoa(), field);
  2665.             this.nMax = nMax;
  2666.             this.sMax = sMax;
  2667.             this.jMax = jMax;

  2669.             // compute the common factors that depends on the mean elements
  2670.             this.hXXX = auxiliaryElements.getH().multiply(context.getChi3());
  2671.             this.kXXX = auxiliaryElements.getK().multiply(context.getChi3());

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

  2676.             for (int s = 0; s <= sMax; s++) {
  2677.                 //Initialise the Hansen roots
  2678.                 hansenObjects.computeHansenObjectsInitValues(context, s);
  2679.             }
  2680.             generateCoefficients(date, context, auxiliaryElements, hansenObjects, field);
  2681.         }

  2683.         /** Generate all coefficients.
  2684.          * @param date the current date
  2685.          * @param context container for attributes
  2686.          * @param hansenObjects initialization of hansen objects
  2687.          * @param auxiliaryElements auxiliary elements related to the current orbit
  2688.          * @param field field used by default
  2689.          */
  2690.         private void generateCoefficients(final FieldAbsoluteDate<T> date,
  2691.                                           final FieldDSSTZonalContext<T> context,
  2692.                                           final FieldAuxiliaryElements<T> auxiliaryElements,
  2693.                                           final FieldHansenObjects<T> hansenObjects,
  2694.                                           final Field<T> field) {

  2696.             //Zero
  2697.             final T zero = field.getZero();

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

  2702.                 //init arrays
  2703.                 for (int i = 0; i <= 6; i++) {
  2704.                     cCoef[i][j] = zero;
  2705.                     sCoef[i][j] = zero;
  2706.                 }

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

  2710.                     //compute first double sum where s: j -> N-1 and n: s+1 -> N
  2711.                     for (int s = j; s <= FastMath.min(nMax - 1, sMax); s++) {
  2712.                         // j - s
  2713.                         final int jms = j - s;
  2714.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2715.                         final int d0smj = (s == j) ? 1 : 2;

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

  2724.                                 final T hjs        = ghijCoef.getHjs(s, -jms);
  2725.                                 final T dHjsdh     = ghijCoef.getdHjsdh(s, -jms);
  2726.                                 final T dHjsdk     = ghijCoef.getdHjsdk(s, -jms);
  2727.                                 final T dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, -jms);
  2728.                                 final T dHjsdBeta  = ghijCoef.getdHjsdBeta(s, -jms);

  2730.                                 final T gjs        = ghijCoef.getGjs(s, -jms);
  2731.                                 final T dGjsdh     = ghijCoef.getdGjsdh(s, -jms);
  2732.                                 final T dGjsdk     = ghijCoef.getdGjsdk(s, -jms);
  2733.                                 final T dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, -jms);
  2734.                                 final T dGjsdBeta  = ghijCoef.getdGjsdBeta(s, -jms);

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

  2739.                                 // K₀<sup>-n-1,s</sup>
  2740.                                 final T kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  2741.                                 final T dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  2743.                                 final T coef0 = jn.multiply(d0smj);
  2744.                                 final T coef1 = coef0.multiply(lns);
  2745.                                 final T coef2 = coef1.multiply(kns);
  2746.                                 final T coef3 = coef2.multiply(hjs);
  2747.                                 final T coef4 = coef2.multiply(gjs);

  2749.                                 // Add the term to the coefficients
  2750.                                 cCoef[0][j] = cCoef[0][j].add(coef3);
  2751.                                 cCoef[1][j] = cCoef[1][j].add(coef3.multiply(n + 1));
  2752.                                 cCoef[2][j] = cCoef[2][j].add(coef1.multiply(kns.multiply(dHjsdh).add(hjs.multiply(hXXX).multiply(dkns))));
  2753.                                 cCoef[3][j] = cCoef[3][j].add(coef1.multiply(kns.multiply(dHjsdk).add(hjs.multiply(kXXX).multiply(dkns))));
  2754.                                 cCoef[4][j] = cCoef[4][j].add(coef2.multiply(dHjsdAlpha));
  2755.                                 cCoef[5][j] = cCoef[5][j].add(coef2.multiply(dHjsdBeta));
  2756.                                 cCoef[6][j] = cCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(hjs));

  2758.                                 sCoef[0][j] = sCoef[0][j].add(coef4);
  2759.                                 sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  2760.                                 sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dGjsdh).add(gjs.multiply(hXXX).multiply(dkns))));
  2761.                                 sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dGjsdk).add(gjs.multiply(kXXX).multiply(dkns))));
  2762.                                 sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dGjsdAlpha));
  2763.                                 sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dGjsdBeta));
  2764.                                 sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(gjs));
  2765.                             }
  2766.                         }
  2767.                     }

  2769.                     //compute second double sum where s: 0 -> N-j and n: max(j+s, j+1) -> N
  2770.                     for (int s = 0; s <= FastMath.min(nMax - j, sMax); s++) {
  2771.                         // j + s
  2772.                         final int jps = j + s;
  2773.                         // Kronecker symbols (2 - delta(0,j+s))
  2774.                         final double d0spj = (s == -j) ? 1 : 2;

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

  2783.                                 final T hjs        = ghijCoef.getHjs(s, jps);
  2784.                                 final T dHjsdh     = ghijCoef.getdHjsdh(s, jps);
  2785.                                 final T dHjsdk     = ghijCoef.getdHjsdk(s, jps);
  2786.                                 final T dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, jps);
  2787.                                 final T dHjsdBeta  = ghijCoef.getdHjsdBeta(s, jps);

  2789.                                 final T gjs        = ghijCoef.getGjs(s, jps);
  2790.                                 final T dGjsdh     = ghijCoef.getdGjsdh(s, jps);
  2791.                                 final T dGjsdk     = ghijCoef.getdGjsdk(s, jps);
  2792.                                 final T dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, jps);
  2793.                                 final T dGjsdBeta  = ghijCoef.getdGjsdBeta(s, jps);

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

  2798.                                 // K₀<sup>-n-1,s</sup>
  2799.                                 final T kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  2800.                                 final T dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  2802.                                 final T coef0 = jn.multiply(d0spj);
  2803.                                 final T coef1 = coef0.multiply(lns);
  2804.                                 final T coef2 = coef1.multiply(kns);

  2806.                                 final T coef3 = coef2.multiply(hjs);
  2807.                                 final T coef4 = coef2.multiply(gjs);

  2809.                                 // Add the term to the coefficients
  2810.                                 cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2811.                                 cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  2812.                                 cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dHjsdh).add(hjs.multiply(hXXX).multiply(dkns))));
  2813.                                 cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dHjsdk).add(hjs.multiply(kXXX).multiply(dkns))));
  2814.                                 cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dHjsdAlpha));
  2815.                                 cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dHjsdBeta));
  2816.                                 cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(hjs));

  2818.                                 sCoef[0][j] = sCoef[0][j].add(coef4);
  2819.                                 sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  2820.                                 sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dGjsdh).add(gjs.multiply(hXXX).multiply(dkns))));
  2821.                                 sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dGjsdk).add(gjs.multiply(kXXX).multiply(dkns))));
  2822.                                 sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dGjsdAlpha));
  2823.                                 sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dGjsdBeta));
  2824.                                 sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(gjs));
  2825.                             }
  2826.                         }
  2827.                     }

  2829.                     //compute third double sum where s: 1 -> j and  n: j+1 -> N
  2830.                     for (int s = 1; s <= FastMath.min(j, sMax); s++) {
  2831.                         // j - s
  2832.                         final int jms = j - s;
  2833.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2834.                         final int d0smj = (s == j) ? 1 : 2;

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

  2843.                                 final T ijs        = ghijCoef.getIjs(s, jms);
  2844.                                 final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  2845.                                 final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  2846.                                 final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2847.                                 final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  2849.                                 final T jjs        = ghijCoef.getJjs(s, jms);
  2850.                                 final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  2851.                                 final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  2852.                                 final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2853.                                 final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

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

  2858.                                 // K₀<sup>-n-1,s</sup>
  2859.                                 final T kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  2860.                                 final T dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  2862.                                 final T coef0 = jn.multiply(d0smj);
  2863.                                 final T coef1 = coef0.multiply(lns);
  2864.                                 final T coef2 = coef1.multiply(kns);

  2866.                                 final T coef3 = coef2.multiply(ijs);
  2867.                                 final T coef4 = coef2.multiply(jjs);

  2869.                                 // Add the term to the coefficients
  2870.                                 cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2871.                                 cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  2872.                                 cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  2873.                                 cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  2874.                                 cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  2875.                                 cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  2876.                                 cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  2878.                                 sCoef[0][j] = sCoef[0][j].add(coef4);
  2879.                                 sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  2880.                                 sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  2881.                                 sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  2882.                                 sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  2883.                                 sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  2884.                                 sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  2885.                             }
  2886.                         }
  2887.                     }
  2888.                 }

  2890.                 if (isBetween(j, 2, nMax)) {
  2891.                     //add first term
  2892.                     // J<sub>j</sub>
  2893.                     final T jj = zero.subtract(harmonics.getUnnormalizedCnm(j, 0));
  2894.                     T kns  = hansenObjects.getHansenObjects()[0].getValue(-j - 1, context.getChi());
  2895.                     T dkns = hansenObjects.getHansenObjects()[0].getDerivative(-j - 1, context.getChi());

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

  2900.                     final T hjs        = ghijCoef.getHjs(0, j);
  2901.                     final T dHjsdh     = ghijCoef.getdHjsdh(0, j);
  2902.                     final T dHjsdk     = ghijCoef.getdHjsdk(0, j);
  2903.                     final T dHjsdAlpha = ghijCoef.getdHjsdAlpha(0, j);
  2904.                     final T dHjsdBeta  = ghijCoef.getdHjsdBeta(0, j);

  2906.                     final T gjs        = ghijCoef.getGjs(0, j);
  2907.                     final T dGjsdh     = ghijCoef.getdGjsdh(0, j);
  2908.                     final T dGjsdk     = ghijCoef.getdGjsdk(0, j);
  2909.                     final T dGjsdAlpha = ghijCoef.getdGjsdAlpha(0, j);
  2910.                     final T dGjsdBeta  = ghijCoef.getdGjsdBeta(0, j);

  2912.                     // 2 * J<sub>j</sub> * K₀<sup>-j-1,0</sup> * L<sub>j</sub><sup>j</sup>
  2913.                     T coef0 = jj.multiply(2.);
  2914.                     T coef1 = coef0.multiply(lns);
  2915.                     T coef2 = coef1.multiply(kns);

  2917.                     T coef3 = coef2.multiply(hjs);
  2918.                     T coef4 = coef2.multiply(gjs);

  2920.                     // Add the term to the coefficients
  2921.                     cCoef[0][j] = cCoef[0][j].subtract(coef3);
  2922.                     cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(j + 1));
  2923.                     cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dHjsdh).add(hjs.multiply(hXXX).multiply(dkns))));
  2924.                     cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dHjsdk).add(hjs.multiply(kXXX).multiply(dkns))));
  2925.                     cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dHjsdAlpha));
  2926.                     cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dHjsdBeta));
  2927.                     //no contribution to cCoef[6][j] because dlns is 0

  2929.                     sCoef[0][j] = sCoef[0][j].add(coef4);
  2930.                     sCoef[1][j] = sCoef[1][j].add(coef4.multiply(j + 1));
  2931.                     sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dGjsdh).add(gjs.multiply(hXXX).multiply(dkns))));
  2932.                     sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dGjsdk).add(gjs.multiply(kXXX).multiply(dkns))));
  2933.                     sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dGjsdAlpha));
  2934.                     sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dGjsdBeta));
  2935.                     //no contribution to sCoef[6][j] because dlns is 0

  2937.                     //compute simple sum where s: 1 -> j-1
  2938.                     for (int s = 1; s <= FastMath.min(j - 1, sMax); s++) {
  2939.                         // j - s
  2940.                         final int jms = j - s;
  2941.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  2942.                         final int d0smj = (s == j) ? 1 : 2;

  2944.                         // if s is odd, then the term is equal to zero due to the factor Vj,s-j
  2945.                         if (s % 2 == 0) {
  2946.                             // (2 - delta(0,j-s)) * J<sub>j</sub> * K₀<sup>-j-1,s</sup> * L<sub>j</sub><sup>j-s</sup>
  2947.                             kns   = hansenObjects.getHansenObjects()[s].getValue(-j - 1, context.getChi());
  2948.                             dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-j - 1, context.getChi());

  2950.                             lns = lnsCoef.getLns(j, jms);
  2951.                             final T dlns = lnsCoef.getdLnsdGamma(j, jms);

  2953.                             final T ijs        = ghijCoef.getIjs(s, jms);
  2954.                             final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  2955.                             final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  2956.                             final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  2957.                             final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  2959.                             final T jjs        = ghijCoef.getJjs(s, jms);
  2960.                             final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  2961.                             final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  2962.                             final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  2963.                             final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

  2965.                             coef0 = jj.multiply(d0smj);
  2966.                             coef1 = coef0.multiply(lns);
  2967.                             coef2 = coef1.multiply(kns);

  2969.                             coef3 = coef2.multiply(ijs);
  2970.                             coef4 = coef2.multiply(jjs);

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

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

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

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

  2998.                     //if j is even
  2999.                     if (j % 2 == 0) {
  3000.                         //compute first double sum where s: j-min(j-1,N) -> j/2-1 and n: j-s -> min(j-1,N)
  3001.                         for (int s = j - minjm1on; s <= FastMath.min(j / 2 - 1, sMax); s++) {
  3002.                             // j - s
  3003.                             final int jms = j - s;
  3004.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  3005.                             final int d0smj = (s == j) ? 1 : 2;

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

  3014.                                     final T ijs        = ghijCoef.getIjs(s, jms);
  3015.                                     final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  3016.                                     final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  3017.                                     final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  3018.                                     final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  3020.                                     final T jjs        = ghijCoef.getJjs(s, jms);
  3021.                                     final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  3022.                                     final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  3023.                                     final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  3024.                                     final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

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

  3029.                                     // K₀<sup>-n-1,s</sup>
  3030.                                     final T kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  3031.                                     final T dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  3033.                                     final T coef0 = jn.multiply(d0smj);
  3034.                                     final T coef1 = coef0.multiply(lns);
  3035.                                     final T coef2 = coef1.multiply(kns);

  3037.                                     final T coef3 = coef2.multiply(ijs);
  3038.                                     final T coef4 = coef2.multiply(jjs);

  3040.                                     // Add the term to the coefficients
  3041.                                     cCoef[0][j] = cCoef[0][j].subtract(coef3);
  3042.                                     cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  3043.                                     cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  3044.                                     cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  3045.                                     cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  3046.                                     cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  3047.                                     cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  3049.                                     sCoef[0][j] = sCoef[0][j].add(coef4);
  3050.                                     sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  3051.                                     sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  3052.                                     sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  3053.                                     sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  3054.                                     sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  3055.                                     sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  3056.                                 }
  3057.                             }
  3058.                         }

  3060.                         //compute second double sum where s: j/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  3061.                         for (int s = j / 2; s <=  FastMath.min(minjm1on - 1, sMax); s++) {
  3062.                             // j - s
  3063.                             final int jms = j - s;
  3064.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  3065.                             final int d0smj = (s == j) ? 1 : 2;

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

  3074.                                     final T ijs        = ghijCoef.getIjs(s, jms);
  3075.                                     final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  3076.                                     final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  3077.                                     final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  3078.                                     final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  3080.                                     final T jjs        = ghijCoef.getJjs(s, jms);
  3081.                                     final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  3082.                                     final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  3083.                                     final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  3084.                                     final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

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

  3089.                                     // K₀<sup>-n-1,s</sup>
  3090.                                     final T kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  3091.                                     final T dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

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

  3097.                                     final T coef3 = coef2.multiply(ijs);
  3098.                                     final T coef4 = coef2.multiply(jjs);

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

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

  3121.                     //if j is odd
  3122.                     else {
  3123.                         //compute first double sum where s: (j-1)/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  3124.                         for (int s = (j - 1) / 2; s <= FastMath.min(minjm1on - 1, sMax); s++) {
  3125.                             // j - s
  3126.                             final int jms = j - s;
  3127.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  3128.                             final int d0smj = (s == j) ? 1 : 2;

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

  3137.                                     final T ijs        = ghijCoef.getIjs(s, jms);
  3138.                                     final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  3139.                                     final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  3140.                                     final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  3141.                                     final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  3143.                                     final T jjs        = ghijCoef.getJjs(s, jms);
  3144.                                     final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  3145.                                     final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  3146.                                     final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  3147.                                     final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

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

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

  3154.                                     final T kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  3155.                                     final T dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  3157.                                     final T coef0 = jn.multiply(d0smj);
  3158.                                     final T coef1 = coef0.multiply(lns);
  3159.                                     final T coef2 = coef1.multiply(kns);

  3161.                                     final T coef3 = coef2.multiply(ijs);
  3162.                                     final T coef4 = coef2.multiply(jjs);

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

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

  3184.                         //the second double sum is added only if N >= 4 and j between 5 and 2*N-3
  3185.                         if (nMax >= 4 && isBetween(j, 5, 2 * nMax - 3)) {
  3186.                             //compute second double sum where s: j-min(j-1,N) -> (j-3)/2 and n: j-s -> min(j-1,N)
  3187.                             for (int s = j - minjm1on; s <= FastMath.min((j - 3) / 2, sMax); s++) {
  3188.                                 // j - s
  3189.                                 final int jms = j - s;
  3190.                                 // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  3191.                                 final int d0smj = (s == j) ? 1 : 2;

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

  3200.                                         final T ijs        = ghijCoef.getIjs(s, jms);
  3201.                                         final T dIjsdh     = ghijCoef.getdIjsdh(s, jms);
  3202.                                         final T dIjsdk     = ghijCoef.getdIjsdk(s, jms);
  3203.                                         final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  3204.                                         final T dIjsdBeta  = ghijCoef.getdIjsdBeta(s, jms);

  3206.                                         final T jjs        = ghijCoef.getJjs(s, jms);
  3207.                                         final T dJjsdh     = ghijCoef.getdJjsdh(s, jms);
  3208.                                         final T dJjsdk     = ghijCoef.getdJjsdk(s, jms);
  3209.                                         final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  3210.                                         final T dJjsdBeta  = ghijCoef.getdJjsdBeta(s, jms);

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

  3215.                                         // K₀<sup>-n-1,s</sup>
  3216.                                         final T kns   = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  3217.                                         final T dkns  = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  3219.                                         final T coef0 = jn.multiply(d0smj);
  3220.                                         final T coef1 = coef0.multiply(lns);
  3221.                                         final T coef2 = coef1.multiply(kns);

  3223.                                         final T coef3 = coef2.multiply(ijs);
  3224.                                         final T coef4 = coef2.multiply(jjs);

  3226.                                         // Add the term to the coefficients
  3227.                                         cCoef[0][j] = cCoef[0][j].subtract(coef3);
  3228.                                         cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
  3229.                                         cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
  3230.                                         cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
  3231.                                         cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
  3232.                                         cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
  3233.                                         cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));

  3235.                                         sCoef[0][j] = sCoef[0][j].add(coef4);
  3236.                                         sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
  3237.                                         sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
  3238.                                         sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
  3239.                                         sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
  3240.                                         sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
  3241.                                         sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
  3242.                                     }
  3243.                                 }
  3244.                             }
  3245.                         }
  3246.                     }
  3247.                 }

  3249.                 cCoef[0][j] = cCoef[0][j].multiply(context.getMuoa().divide(j).negate());
  3250.                 cCoef[1][j] = cCoef[1][j].multiply(context.getMuoa().divide(auxiliaryElements.getSma().multiply(j)));
  3251.                 cCoef[2][j] = cCoef[2][j].multiply(context.getMuoa().divide(j).negate());
  3252.                 cCoef[3][j] = cCoef[3][j].multiply(context.getMuoa().divide(j).negate());
  3253.                 cCoef[4][j] = cCoef[4][j].multiply(context.getMuoa().divide(j).negate());
  3254.                 cCoef[5][j] = cCoef[5][j].multiply(context.getMuoa().divide(j).negate());
  3255.                 cCoef[6][j] = cCoef[6][j].multiply(context.getMuoa().divide(j).negate());

  3257.                 sCoef[0][j] = sCoef[0][j].multiply(context.getMuoa().divide(j).negate());
  3258.                 sCoef[1][j] = sCoef[1][j].multiply(context.getMuoa().divide(auxiliaryElements.getSma().multiply(j)));
  3259.                 sCoef[2][j] = sCoef[2][j].multiply(context.getMuoa().divide(j).negate());
  3260.                 sCoef[3][j] = sCoef[3][j].multiply(context.getMuoa().divide(j).negate());
  3261.                 sCoef[4][j] = sCoef[4][j].multiply(context.getMuoa().divide(j).negate());
  3262.                 sCoef[5][j] = sCoef[5][j].multiply(context.getMuoa().divide(j).negate());
  3263.                 sCoef[6][j] = sCoef[6][j].multiply(context.getMuoa().divide(j).negate());

  3265.             }
  3266.         }

  3268.         /** Check if an index is within the accepted interval.
  3269.          *
  3270.          * @param index the index to check
  3271.          * @param lowerBound the lower bound of the interval
  3272.          * @param upperBound the upper bound of the interval
  3273.          * @return true if the index is between the lower and upper bounds, false otherwise
  3274.          */
  3275.         private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
  3276.             return index >= lowerBound && index <= upperBound;
  3277.         }

  3279.         /**Get the value of C<sup>j</sup>.
  3280.          *
  3281.          * @param j j index
  3282.          * @return C<sup>j</sup>
  3283.          */
  3284.         public T getCj(final int j) {
  3285.             return cCoef[0][j];
  3286.         }

  3288.         /**Get the value of dC<sup>j</sup> / da.
  3289.          *
  3290.          * @param j j index
  3291.          * @return dC<sup>j</sup> / da
  3292.          */
  3293.         public T getdCjdA(final int j) {
  3294.             return cCoef[1][j];
  3295.         }

  3297.         /**Get the value of dC<sup>j</sup> / dh.
  3298.          *
  3299.          * @param j j index
  3300.          * @return dC<sup>j</sup> / dh
  3301.          */
  3302.         public T getdCjdH(final int j) {
  3303.             return cCoef[2][j];
  3304.         }

  3306.         /**Get the value of dC<sup>j</sup> / dk.
  3307.          *
  3308.          * @param j j index
  3309.          * @return dC<sup>j</sup> / dk
  3310.          */
  3311.         public T getdCjdK(final int j) {
  3312.             return cCoef[3][j];
  3313.         }

  3315.         /**Get the value of dC<sup>j</sup> / dα.
  3316.          *
  3317.          * @param j j index
  3318.          * @return dC<sup>j</sup> / dα
  3319.          */
  3320.         public T getdCjdAlpha(final int j) {
  3321.             return cCoef[4][j];
  3322.         }

  3324.         /**Get the value of dC<sup>j</sup> / dβ.
  3325.          *
  3326.          * @param j j index
  3327.          * @return dC<sup>j</sup> / dβ
  3328.          */
  3329.         public T getdCjdBeta(final int j) {
  3330.             return cCoef[5][j];
  3331.         }

  3333.         /**Get the value of dC<sup>j</sup> / dγ.
  3334.          *
  3335.          * @param j j index
  3336.          * @return dC<sup>j</sup> / dγ
  3337.          */
  3338.         public T getdCjdGamma(final int j) {
  3339.             return cCoef[6][j];
  3340.         }

  3342.         /**Get the value of S<sup>j</sup>.
  3343.          *
  3344.          * @param j j index
  3345.          * @return S<sup>j</sup>
  3346.          */
  3347.         public T getSj(final int j) {
  3348.             return sCoef[0][j];
  3349.         }

  3351.         /**Get the value of dS<sup>j</sup> / da.
  3352.          *
  3353.          * @param j j index
  3354.          * @return dS<sup>j</sup> / da
  3355.          */
  3356.         public T getdSjdA(final int j) {
  3357.             return sCoef[1][j];
  3358.         }

  3360.         /**Get the value of dS<sup>j</sup> / dh.
  3361.          *
  3362.          * @param j j index
  3363.          * @return dS<sup>j</sup> / dh
  3364.          */
  3365.         public T getdSjdH(final int j) {
  3366.             return sCoef[2][j];
  3367.         }

  3369.         /**Get the value of dS<sup>j</sup> / dk.
  3370.          *
  3371.          * @param j j index
  3372.          * @return dS<sup>j</sup> / dk
  3373.          */
  3374.         public T getdSjdK(final int j) {
  3375.             return sCoef[3][j];
  3376.         }

  3378.         /**Get the value of dS<sup>j</sup> / dα.
  3379.          *
  3380.          * @param j j index
  3381.          * @return dS<sup>j</sup> / dα
  3382.          */
  3383.         public T getdSjdAlpha(final int j) {
  3384.             return sCoef[4][j];
  3385.         }

  3387.         /**Get the value of dS<sup>j</sup> / dβ.
  3388.          *
  3389.          * @param j j index
  3390.          * @return dS<sup>j</sup> / dβ
  3391.          */
  3392.         public T getdSjdBeta(final int j) {
  3393.             return sCoef[5][j];
  3394.         }

  3396.         /**Get the value of dS<sup>j</sup> /  dγ.
  3397.          *
  3398.          * @param j j index
  3399.          * @return dS<sup>j</sup> /  dγ
  3400.          */
  3401.         public T getdSjdGamma(final int j) {
  3402.             return sCoef[6][j];
  3403.         }
  3404.     }

  3406.     /** Coefficients valid for one time slot. */
  3407.     private static class Slot {

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

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

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

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

  3457.             final int rows = maxFrequencyShortPeriodics + 1;
  3458.             di  = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3459.             cij = new ShortPeriodicsInterpolatedCoefficient[rows];
  3460.             sij = new ShortPeriodicsInterpolatedCoefficient[rows];

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

  3468.         }

  3470.     }

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

  3475.         /**The coefficients D<sub>i</sub>.
  3476.          * <p>
  3477.          * i corresponds to the equinoctial element, as follows:
  3478.          * - i=0 for a <br/>
  3479.          * - i=1 for k <br/>
  3480.          * - i=2 for h <br/>
  3481.          * - i=3 for q <br/>
  3482.          * - i=4 for p <br/>
  3483.          * - i=5 for λ <br/>
  3484.          * </p>
  3485.          */
  3486.         private final FieldShortPeriodicsInterpolatedCoefficient<T> di;

  3488.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  3489.          * <p>
  3490.          * The constant term C<sub>i</sub>⁰ is also stored in this variable at index j = 0 <br>
  3491.          * The index order is cij[j][i] <br/>
  3492.          * i corresponds to the equinoctial element, as follows: <br/>
  3493.          * - i=0 for a <br/>
  3494.          * - i=1 for k <br/>
  3495.          * - i=2 for h <br/>
  3496.          * - i=3 for q <br/>
  3497.          * - i=4 for p <br/>
  3498.          * - i=5 for λ <br/>
  3499.          * </p>
  3500.          */
  3501.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] cij;

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

  3517.         /** Simple constructor.
  3518.          *  @param maxFrequencyShortPeriodics maximum value for j index
  3519.          *  @param interpolationPoints number of points used in the interpolation process
  3520.          */
  3521.         @SuppressWarnings("unchecked")
  3522.         FieldSlot(final int maxFrequencyShortPeriodics, final int interpolationPoints) {

  3524.             final int rows = maxFrequencyShortPeriodics + 1;
  3525.             di  = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3526.             cij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows);
  3527.             sij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows);

  3529.             //Initialize the arrays
  3530.             for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
  3531.                 cij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3532.                 sij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3533.             }

  3535.         }

  3537.     }

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

  3542.         /** The current value of the U function. <br/>
  3543.          * Needed for the short periodic contribution */
  3544.         private double U;

  3546.         /** dU / da. */
  3547.         private double dUda;

  3549.         /** dU / dk. */
  3550.         private double dUdk;

  3552.         /** dU / dh. */
  3553.         private double dUdh;

  3555.         /** dU / dAlpha. */
  3556.         private double dUdAl;

  3558.         /** dU / dBeta. */
  3559.         private double dUdBe;

  3561.         /** dU / dGamma. */
  3562.         private double dUdGa;

  3564.         /** Simple constuctor.
  3565.          *  @param date current date
  3566.          *  @param context container for attributes
  3567.          *  @param auxiliaryElements auxiliary elements related to the current orbit
  3568.          *  @param hansen initialization of hansen objects
  3569.          */
  3570.         UAnddU(final AbsoluteDate date,
  3571.                final DSSTZonalContext context,
  3572.                final AuxiliaryElements auxiliaryElements,
  3573.                final HansenObjects hansen) {

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

  3577.             //Reset U
  3578.             U = 0.;

  3580.             // Gs and Hs coefficients
  3581.             final double[][] GsHs = CoefficientsFactory.computeGsHs(auxiliaryElements.getK(), auxiliaryElements.getH(), context.getAlpha(), context.getBeta(), maxEccPowMeanElements);
  3582.             // Qns coefficients
  3583.             final double[][] Qns  = CoefficientsFactory.computeQns(context.getGamma(), maxDegree, maxEccPowMeanElements);

  3585.             final double[] roaPow = new double[maxDegree + 1];
  3586.             roaPow[0] = 1.;
  3587.             for (int i = 1; i <= maxDegree; i++) {
  3588.                 roaPow[i] = context.getRoa() * roaPow[i - 1];
  3589.             }

  3591.             // Potential derivatives
  3592.             dUda  = 0.;
  3593.             dUdk  = 0.;
  3594.             dUdh  = 0.;
  3595.             dUdAl = 0.;
  3596.             dUdBe = 0.;
  3597.             dUdGa = 0.;

  3599.             for (int s = 0; s <= maxEccPowMeanElements; s++) {
  3600.                 //Initialize the Hansen roots
  3601.                 hansen.computeHansenObjectsInitValues(context, s);

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

  3606.                 // Compute Gs partial derivatives from 3.1-(9)
  3607.                 double dGsdh  = 0.;
  3608.                 double dGsdk  = 0.;
  3609.                 double dGsdAl = 0.;
  3610.                 double dGsdBe = 0.;
  3611.                 if (s > 0) {
  3612.                     // First get the G(s-1) and the H(s-1) coefficients
  3613.                     final double sxgsm1 = s * GsHs[0][s - 1];
  3614.                     final double sxhsm1 = s * GsHs[1][s - 1];
  3615.                     // Then compute derivatives
  3616.                     dGsdh  = context.getBeta()  * sxgsm1 - context.getAlpha() * sxhsm1;
  3617.                     dGsdk  = context.getAlpha() * sxgsm1 + context.getBeta()  * sxhsm1;
  3618.                     dGsdAl = auxiliaryElements.getK() * sxgsm1 - auxiliaryElements.getH() * sxhsm1;
  3619.                     dGsdBe = auxiliaryElements.getH() * sxgsm1 + auxiliaryElements.getK() * sxhsm1;
  3620.                 }

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

  3625.                 for (int n = s + 2; n <= maxDegree; n++) {
  3626.                     // (n - s) must be even
  3627.                     if ((n - s) % 2 == 0) {

  3629.                         //Extract data from previous computation :
  3630.                         final double kns   = hansen.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  3631.                         final double dkns  = hansen.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  3633.                         final double vns   = Vns.get(new NSKey(n, s));
  3634.                         final double coef0 = d0s * roaPow[n] * vns * -harmonics.getUnnormalizedCnm(n, 0);
  3635.                         final double coef1 = coef0 * Qns[n][s];
  3636.                         final double coef2 = coef1 * kns;
  3637.                         final double coef3 = coef2 * gs;
  3638.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  3639.                         final double dqns  = Qns[n][s + 1];

  3641.                         // Compute U
  3642.                         U += coef3;
  3643.                         // Compute dU / da :
  3644.                         dUda += coef3 * (n + 1);
  3645.                         // Compute dU / dEx
  3646.                         dUdk += coef1 * (kns * dGsdk + auxiliaryElements.getK() * context.getChi3() * gs * dkns);
  3647.                         // Compute dU / dEy
  3648.                         dUdh += coef1 * (kns * dGsdh + auxiliaryElements.getH() * context.getChi3() * gs * dkns);
  3649.                         // Compute dU / dAlpha
  3650.                         dUdAl += coef2 * dGsdAl;
  3651.                         // Compute dU / dBeta
  3652.                         dUdBe += coef2 * dGsdBe;
  3653.                         // Compute dU / dGamma
  3654.                         dUdGa += coef0 * kns * dqns * gs;

  3656.                     }
  3657.                 }
  3658.             }

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

  3663.             this.dUda = dUda *  context.getMuoa() / auxiliaryElements.getSma();
  3664.             this.dUdk = dUdk * -context.getMuoa();
  3665.             this.dUdh = dUdh * -context.getMuoa();
  3666.             this.dUdAl = dUdAl * -context.getMuoa();
  3667.             this.dUdBe = dUdBe * -context.getMuoa();
  3668.             this.dUdGa = dUdGa * -context.getMuoa();

  3670.         }

  3672.         /** Return value of U.
  3673.          * @return U
  3674.          */
  3675.         public double getU() {
  3676.             return U;
  3677.         }

  3679.         /** Return value of dU / da.
  3680.          * @return dUda
  3681.          */
  3682.         public double getdUda() {
  3683.             return dUda;
  3684.         }

  3686.         /** Return value of dU / dk.
  3687.          * @return dUdk
  3688.          */
  3689.         public double getdUdk() {
  3690.             return dUdk;
  3691.         }

  3693.         /** Return value of dU / dh.
  3694.          * @return dUdh
  3695.          */
  3696.         public double getdUdh() {
  3697.             return dUdh;
  3698.         }

  3700.         /** Return value of dU / dAlpha.
  3701.          * @return dUdAl
  3702.          */
  3703.         public double getdUdAl() {
  3704.             return dUdAl;
  3705.         }

  3707.         /** Return value of dU / dBeta.
  3708.          * @return dUdBe
  3709.          */
  3710.         public double getdUdBe() {
  3711.             return dUdBe;
  3712.         }

  3714.         /** Return value of dU / dGamma.
  3715.          * @return dUdGa
  3716.          */
  3717.         public double getdUdGa() {
  3718.             return dUdGa;
  3719.         }

  3721.     }

  3723.     /** Compute the derivatives of the gravitational potential U [Eq. 3.1-(6)].
  3724.      *  <p>
  3725.      *  The result is the array
  3726.      *  [dU/da, dU/dk, dU/dh, dU/dα, dU/dβ, dU/dγ]
  3727.      *  </p>
  3728.      */
  3729.     private class FieldUAnddU <T extends CalculusFieldElement<T>> {

  3731.          /** The current value of the U function. <br/>
  3732.           * Needed for the short periodic contribution */
  3733.         private T U;

  3735.          /** dU / da. */
  3736.         private T dUda;

  3738.          /** dU / dk. */
  3739.         private T dUdk;

  3741.          /** dU / dh. */
  3742.         private T dUdh;

  3744.          /** dU / dAlpha. */
  3745.         private T dUdAl;

  3747.          /** dU / dBeta. */
  3748.         private T dUdBe;

  3750.          /** dU / dGamma. */
  3751.         private T dUdGa;

  3753.          /** Simple constuctor.
  3754.           *  @param date current date
  3755.           *  @param context container for attributes
  3756.           *  @param auxiliaryElements auxiliary elements related to the current orbit
  3757.           *  @param hansen initialization of hansen objects
  3758.           */
  3759.         FieldUAnddU(final FieldAbsoluteDate<T> date,
  3760.                     final FieldDSSTZonalContext<T> context,
  3761.                     final FieldAuxiliaryElements<T> auxiliaryElements,
  3762.                     final FieldHansenObjects<T> hansen) {

  3764.             // Zero for initialization
  3765.             final Field<T> field = date.getField();
  3766.             final T zero = field.getZero();

  3768.             //spherical harmonics
  3769.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date.toAbsoluteDate());

  3771.             //Reset U
  3772.             U = zero;

  3774.             // Gs and Hs coefficients
  3775.             final T[][] GsHs = CoefficientsFactory.computeGsHs(auxiliaryElements.getK(), auxiliaryElements.getH(),
  3776.                                                                context.getAlpha(), context.getBeta(),
  3777.                                                                maxEccPowMeanElements, field);
  3778.             // Qns coefficients
  3779.             final T[][] Qns  = CoefficientsFactory.computeQns(context.getGamma(), maxDegree, maxEccPowMeanElements);

  3781.             final T[] roaPow = MathArrays.buildArray(field, maxDegree + 1);
  3782.             roaPow[0] = zero.newInstance(1.);
  3783.             for (int i = 1; i <= maxDegree; i++) {
  3784.                 roaPow[i] = roaPow[i - 1].multiply(context.getRoa());
  3785.             }

  3787.             // Potential derivatives
  3788.             dUda  = zero;
  3789.             dUdk  = zero;
  3790.             dUdh  = zero;
  3791.             dUdAl = zero;
  3792.             dUdBe = zero;
  3793.             dUdGa = zero;

  3795.             for (int s = 0; s <= maxEccPowMeanElements; s++) {
  3796.                 //Initialize the Hansen roots
  3797.                 hansen.computeHansenObjectsInitValues(context, s);

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

  3802.                 // Compute Gs partial derivatives from 3.1-(9)
  3803.                 T dGsdh  = zero;
  3804.                 T dGsdk  = zero;
  3805.                 T dGsdAl = zero;
  3806.                 T dGsdBe = zero;
  3807.                 if (s > 0) {
  3808.                     // First get the G(s-1) and the H(s-1) coefficients
  3809.                     final T sxgsm1 = GsHs[0][s - 1].multiply(s);
  3810.                     final T sxhsm1 = GsHs[1][s - 1].multiply(s);
  3811.                     // Then compute derivatives
  3812.                     dGsdh  = sxgsm1.multiply(context.getBeta()).subtract(sxhsm1.multiply(context.getAlpha()));
  3813.                     dGsdk  = sxgsm1.multiply(context.getAlpha()).add(sxhsm1.multiply(context.getBeta()));
  3814.                     dGsdAl = sxgsm1.multiply(auxiliaryElements.getK()).subtract(sxhsm1.multiply(auxiliaryElements.getH()));
  3815.                     dGsdBe = sxgsm1.multiply(auxiliaryElements.getH()).add(sxhsm1.multiply(auxiliaryElements.getK()));
  3816.                 }

  3818.                 // Kronecker symbol (2 - delta(0,s))
  3819.                 final T d0s = zero.newInstance((s == 0) ? 1 : 2);

  3821.                 for (int n = s + 2; n <= maxDegree; n++) {
  3822.                     // (n - s) must be even
  3823.                     if ((n - s) % 2 == 0) {

  3825.                         //Extract data from previous computation :
  3826.                         final T kns   = hansen.getHansenObjects()[s].getValue(-n - 1, context.getChi());
  3827.                         final T dkns  = hansen.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());

  3829.                         final double vns   = Vns.get(new NSKey(n, s));
  3830.                         final T coef0 = d0s.multiply(roaPow[n]).multiply(vns).multiply(-harmonics.getUnnormalizedCnm(n, 0));
  3831.                         final T coef1 = coef0.multiply(Qns[n][s]);
  3832.                         final T coef2 = coef1.multiply(kns);
  3833.                         final T coef3 = coef2.multiply(gs);
  3834.                         // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  3835.                         final T dqns  = Qns[n][s + 1];

  3837.                         // Compute U
  3838.                         U = U.add(coef3);
  3839.                         // Compute dU / da :
  3840.                         dUda  = dUda.add(coef3.multiply(n + 1));
  3841.                         // Compute dU / dEx
  3842.                         dUdk  = dUdk.add(coef1.multiply(dGsdk.multiply(kns).add(auxiliaryElements.getK().multiply(context.getChi3()).multiply(dkns).multiply(gs))));
  3843.                         // Compute dU / dEy
  3844.                         dUdh  = dUdh.add(coef1.multiply(dGsdh.multiply(kns).add(auxiliaryElements.getH().multiply(context.getChi3()).multiply(dkns).multiply(gs))));
  3845.                         // Compute dU / dAlpha
  3846.                         dUdAl = dUdAl.add(coef2.multiply(dGsdAl));
  3847.                         // Compute dU / dBeta
  3848.                         dUdBe = dUdBe.add(coef2.multiply(dGsdBe));
  3849.                         // Compute dU / dGamma
  3850.                         dUdGa = dUdGa.add(coef0.multiply(kns).multiply(dqns).multiply(gs));
  3851.                     }
  3852.                 }
  3853.             }

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

  3858.             dUda  = dUda.multiply(context.getMuoa().divide(auxiliaryElements.getSma()));
  3859.             dUdk  = dUdk.multiply(context.getMuoa()).negate();
  3860.             dUdh  = dUdh.multiply(context.getMuoa()).negate();
  3861.             dUdAl = dUdAl.multiply(context.getMuoa()).negate();
  3862.             dUdBe = dUdBe.multiply(context.getMuoa()).negate();
  3863.             dUdGa = dUdGa.multiply(context.getMuoa()).negate();

  3865.         }

  3867.         /** Return value of U.
  3868.          * @return U
  3869.          */
  3870.         public T getU() {
  3871.             return U;
  3872.         }

  3874.          /** Return value of dU / da.
  3875.           * @return dUda
  3876.           */
  3877.         public T getdUda() {
  3878.             return dUda;
  3879.         }

  3881.          /** Return value of dU / dk.
  3882.           * @return dUdk
  3883.           */
  3884.         public T getdUdk() {
  3885.             return dUdk;
  3886.         }

  3888.          /** Return value of dU / dh.
  3889.           * @return dUdh
  3890.           */
  3891.         public T getdUdh() {
  3892.             return dUdh;
  3893.         }

  3895.          /** Return value of dU / dAlpha.
  3896.           * @return dUdAl
  3897.           */
  3898.         public T getdUdAl() {
  3899.             return dUdAl;
  3900.         }

  3902.          /** Return value of dU / dBeta.
  3903.           * @return dUdBe
  3904.           */
  3905.         public T getdUdBe() {
  3906.             return dUdBe;
  3907.         }

  3909.          /** Return value of dU / dGamma.
  3910.           * @return dUdGa
  3911.           */
  3912.         public T getdUdGa() {
  3913.             return dUdGa;
  3914.         }

  3916.     }

  3918.     /** Computes init values of the Hansen Objects. */
  3919.     private class HansenObjects {

  3921.         /** An array that contains the objects needed to build the Hansen coefficients. <br/>
  3922.          * The index is s*/
  3923.         private final HansenZonalLinear[] hansenObjects;

  3925.         /** Simple constructor. */
  3926.         HansenObjects() {
  3927.             this.hansenObjects = new HansenZonalLinear[maxEccPow + 1];
  3928.             for (int s = 0; s <= maxEccPow; s++) {
  3929.                 this.hansenObjects[s] = new HansenZonalLinear(maxDegree, s);
  3930.             }
  3931.         }

  3933.         /** Compute init values for hansen objects.
  3934.          * @param context container for attributes
  3935.          * @param element element of the array to compute the init values
  3936.          */
  3937.         public void computeHansenObjectsInitValues(final DSSTZonalContext context, final int element) {
  3938.             hansenObjects[element].computeInitValues(context.getChi());
  3939.         }

  3941.         /** Get the Hansen Objects.
  3942.          * @return hansenObjects
  3943.          */
  3944.         public HansenZonalLinear[] getHansenObjects() {
  3945.             return hansenObjects;
  3946.         }

  3948.     }

  3950.     /** Computes init values of the Hansen Objects.
  3951.      * @param <T> type of the elements
  3952.      */
  3953.     private class FieldHansenObjects<T extends CalculusFieldElement<T>> {

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

  3959.         /** Simple constructor.
  3960.          * @param field field used by default
  3961.          */
  3962.         @SuppressWarnings("unchecked")
  3963.         FieldHansenObjects(final Field<T> field) {
  3964.             this.hansenObjects = (FieldHansenZonalLinear<T>[]) Array.newInstance(FieldHansenZonalLinear.class, maxEccPow + 1);
  3965.             for (int s = 0; s <= maxEccPow; s++) {
  3966.                 this.hansenObjects[s] = new FieldHansenZonalLinear<>(maxDegree, s, field);
  3967.             }
  3968.         }

  3970.         /** Compute init values for hansen objects.
  3971.          * @param context container for attributes
  3972.          * @param element element of the array to compute the init values
  3973.          */
  3974.         public void computeHansenObjectsInitValues(final FieldDSSTZonalContext<T> context, final int element) {
  3975.             hansenObjects[element].computeInitValues(context.getChi());
  3976.         }

  3978.         /** Get the Hansen Objects.
  3979.          * @return hansenObjects
  3980.          */
  3981.         public FieldHansenZonalLinear<T>[] getHansenObjects() {
  3982.             return hansenObjects;
  3983.         }

  3985.     }

  3987. }
Created with JaCoCo 0.8.12.202403310830