DSSTTesseral.java

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

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

  30. import org.hipparchus.CalculusFieldElement;
  31. import org.hipparchus.Field;
  32. import org.hipparchus.analysis.differentiation.FieldGradient;
  33. import org.hipparchus.exception.LocalizedCoreFormats;
  34. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  35. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  36. import org.hipparchus.util.FastMath;
  37. import org.hipparchus.util.FieldSinCos;
  38. import org.hipparchus.util.MathArrays;
  39. import org.hipparchus.util.MathUtils;
  40. import org.hipparchus.util.SinCos;
  41. import org.orekit.attitudes.AttitudeProvider;
  42. import org.orekit.errors.OrekitException;
  43. import org.orekit.errors.OrekitInternalError;
  44. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  45. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
  46. import org.orekit.frames.FieldStaticTransform;
  47. import org.orekit.frames.Frame;
  48. import org.orekit.frames.StaticTransform;
  49. import org.orekit.orbits.FieldOrbit;
  50. import org.orekit.orbits.Orbit;
  51. import org.orekit.propagation.FieldSpacecraftState;
  52. import org.orekit.propagation.PropagationType;
  53. import org.orekit.propagation.SpacecraftState;
  54. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  55. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
  56. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  57. import org.orekit.propagation.semianalytical.dsst.utilities.FieldGHmsjPolynomials;
  58. import org.orekit.propagation.semianalytical.dsst.utilities.FieldGammaMnsFunction;
  59. import org.orekit.propagation.semianalytical.dsst.utilities.FieldShortPeriodicsInterpolatedCoefficient;
  60. import org.orekit.propagation.semianalytical.dsst.utilities.GHmsjPolynomials;
  61. import org.orekit.propagation.semianalytical.dsst.utilities.GammaMnsFunction;
  62. import org.orekit.propagation.semianalytical.dsst.utilities.JacobiPolynomials;
  63. import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
  64. import org.orekit.propagation.semianalytical.dsst.utilities.hansen.FieldHansenTesseralLinear;
  65. import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenTesseralLinear;
  66. import org.orekit.time.AbsoluteDate;
  67. import org.orekit.time.FieldAbsoluteDate;
  68. import org.orekit.utils.FieldTimeSpanMap;
  69. import org.orekit.utils.ParameterDriver;
  70. import org.orekit.utils.TimeSpanMap;

  72. /** Tesseral contribution to the central body gravitational perturbation.
  73.  *  <p>
  74.  *  Only resonant tesserals are considered.
  75.  *  </p>
  76.  *
  77.  *  @author Romain Di Costanzo
  78.  *  @author Pascal Parraud
  79.  *  @author Bryan Cazabonne (field translation)
  80.  */
  81. public class DSSTTesseral implements DSSTForceModel {

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

  86.     /** Identifier for cMm coefficients. */
  87.     public static final String CM_COEFFICIENTS = "cM";

  89.     /** Identifier for sMm coefficients. */
  90.     public static final String SM_COEFFICIENTS = "sM";

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

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

  115.     /** Minimum period for analytically averaged high-order resonant
  116.      *  central body spherical harmonics in seconds.
  117.      */
  118.     private static final double MIN_PERIOD_IN_SECONDS = 864000.;

  120.     /** Minimum period for analytically averaged high-order resonant
  121.      *  central body spherical harmonics in satellite revolutions.
  122.      */
  123.     private static final double MIN_PERIOD_IN_SAT_REV = 10.;

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

  128.     /** Provider for spherical harmonics. */
  129.     private final UnnormalizedSphericalHarmonicsProvider provider;

  131.     /** Central body rotating frame. */
  132.     private final Frame bodyFrame;

  134.     /** Central body rotation rate (rad/s). */
  135.     private final double centralBodyRotationRate;

  137.     /** Central body rotation period (seconds). */
  138.     private final double bodyPeriod;

  140.     /** Maximal degree to consider for harmonics potential. */
  141.     private final int maxDegree;

  143.     /** Maximal degree to consider for short periodics tesseral harmonics potential (without m-daily). */
  144.     private final int maxDegreeTesseralSP;

  146.     /** Maximal degree to consider for short periodics m-daily tesseral harmonics potential . */
  147.     private final int maxDegreeMdailyTesseralSP;

  149.     /** Maximal order to consider for harmonics potential. */
  150.     private final int maxOrder;

  152.     /** Maximal order to consider for short periodics tesseral harmonics potential (without m-daily). */
  153.     private final int maxOrderTesseralSP;

  155.     /** Maximal order to consider for short periodics m-daily tesseral harmonics potential . */
  156.     private final int maxOrderMdailyTesseralSP;

  158.     /** Maximum power of the eccentricity to use in summation over s for
  159.      * short periodic tesseral harmonics (without m-daily). */
  160.     private final int maxEccPowTesseralSP;

  162.     /** Maximum power of the eccentricity to use in summation over s for
  163.      * m-daily tesseral harmonics. */
  164.     private final int maxEccPowMdailyTesseralSP;

  166.     /** Maximum value for j. */
  167.     private final int maxFrequencyShortPeriodics;

  169.     /** Maximum power of the eccentricity to use in summation over s. */
  170.     private int maxEccPow;

  172.     /** Maximum power of the eccentricity to use in Hansen coefficient Kernel expansion. */
  173.     private int maxHansen;

  175.     /** Maximum value between maxOrderMdailyTesseralSP and maxOrderTesseralSP. */
  176.     private int mMax;

  178.     /** List of non resonant orders with j != 0. */
  179.     private final SortedMap<Integer, List<Integer> > nonResOrders;

  181.     /** List of resonant orders. */
  182.     private final List<Integer> resOrders;

  184.     /** Short period terms. */
  185.     private TesseralShortPeriodicCoefficients shortPeriodTerms;

  187.     /** Short period terms. */
  188.     private Map<Field<?>, FieldTesseralShortPeriodicCoefficients<?>> fieldShortPeriodTerms;

  190.     /** Driver for gravitational parameter. */
  191.     private final ParameterDriver gmParameterDriver;

  193.     /** Hansen objects. */
  194.     private HansenObjects hansen;

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

  199.     /** Simple constructor with default reference values.
  200.      * <p>
  201.      * When this constructor is used, maximum allowed values are used
  202.      * for the short periodic coefficients:
  203.      * </p>
  204.      * <ul>
  205.      *    <li> {@link #maxDegreeTesseralSP} is set to {@code provider.getMaxDegree()} </li>
  206.      *    <li> {@link #maxOrderTesseralSP} is set to {@code provider.getMaxOrder()}. </li>
  207.      *    <li> {@link #maxEccPowTesseralSP} is set to {@code min(4, provider.getMaxOrder())} </li>
  208.      *    <li> {@link #maxFrequencyShortPeriodics} is set to {@code min(provider.getMaxDegree() + 4, 12)}.
  209.      *         This parameter should not exceed 12 as higher values will exceed computer capacity </li>
  210.      *    <li> {@link #maxDegreeMdailyTesseralSP} is set to {@code provider.getMaxDegree()} </li>
  211.      *    <li> {@link #maxOrderMdailyTesseralSP} is set to {@code provider.getMaxOrder()} </li>
  212.      *    <li> {@link #maxEccPowMdailyTesseralSP} is set to min(provider.getMaxDegree() - 2, 4).
  213.      *         This parameter should not exceed 4 as higher values will exceed computer capacity </li>
  214.      * </ul>
  215.      * @param centralBodyFrame rotating body frame
  216.      * @param centralBodyRotationRate central body rotation rate (rad/s)
  217.      * @param provider provider for spherical harmonics
  218.      * @since 10.1
  219.      */
  220.     public DSSTTesseral(final Frame centralBodyFrame,
  221.                         final double centralBodyRotationRate,
  222.                         final UnnormalizedSphericalHarmonicsProvider provider) {
  223.         this(centralBodyFrame, centralBodyRotationRate, provider, provider.getMaxDegree(),
  224.              provider.getMaxOrder(), FastMath.min(4, provider.getMaxOrder()),  FastMath.min(12, provider.getMaxDegree() + 4),
  225.              provider.getMaxDegree(), provider.getMaxOrder(), FastMath.min(4, provider.getMaxDegree() - 2));
  226.     }

  228.     /** Simple constructor.
  229.      * @param centralBodyFrame rotating body frame
  230.      * @param centralBodyRotationRate central body rotation rate (rad/s)
  231.      * @param provider provider for spherical harmonics
  232.      * @param maxDegreeTesseralSP maximal degree to consider for short periodics tesseral harmonics potential
  233.      *  (must be between 2 and {@code provider.getMaxDegree()})
  234.      * @param maxOrderTesseralSP maximal order to consider for short periodics tesseral harmonics potential
  235.      *  (must be between 0 and {@code provider.getMaxOrder()})
  236.      * @param maxEccPowTesseralSP maximum power of the eccentricity to use in summation over s for
  237.      * short periodic tesseral harmonics (without m-daily), should typically not exceed 4 as higher
  238.      * values will exceed computer capacity
  239.      * (must be between 0 and {@code provider.getMaxOrder()} though, however if order = 0 the value can be anything
  240.      *  since it won't be used in the code)
  241.      * @param maxFrequencyShortPeriodics maximum frequency in mean longitude for short periodic computations
  242.      * (typically {@code maxDegreeTesseralSP} + {@code maxEccPowTesseralSP and no more than 12})
  243.      * @param maxDegreeMdailyTesseralSP maximal degree to consider for short periodics m-daily tesseral harmonics potential
  244.      *  (must be between 2 and {@code provider.getMaxDegree()})
  245.      * @param maxOrderMdailyTesseralSP maximal order to consider for short periodics m-daily tesseral harmonics potential
  246.      *  (must be between 0 and {@code provider.getMaxOrder()})
  247.      * @param maxEccPowMdailyTesseralSP maximum power of the eccentricity to use in summation over s for
  248.      * m-daily tesseral harmonics, (must be between 0 and {@code maxDegreeMdailyTesseralSP - 2},
  249.      * but should typically not exceed 4 as higher values will exceed computer capacity)
  250.      * @since 7.2
  251.      */
  252.     public DSSTTesseral(final Frame centralBodyFrame,
  253.                         final double centralBodyRotationRate,
  254.                         final UnnormalizedSphericalHarmonicsProvider provider,
  255.                         final int maxDegreeTesseralSP, final int maxOrderTesseralSP,
  256.                         final int maxEccPowTesseralSP, final int maxFrequencyShortPeriodics,
  257.                         final int maxDegreeMdailyTesseralSP, final int maxOrderMdailyTesseralSP,
  258.                         final int maxEccPowMdailyTesseralSP) {

  260.         gmParameterDriver = new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
  261.                                                 provider.getMu(), MU_SCALE,
  262.                                                 0.0, Double.POSITIVE_INFINITY);

  264.         // Central body rotating frame
  265.         this.bodyFrame = centralBodyFrame;

  267.         //Save the rotation rate
  268.         this.centralBodyRotationRate = centralBodyRotationRate;

  270.         // Central body rotation period in seconds
  271.         this.bodyPeriod = MathUtils.TWO_PI / centralBodyRotationRate;

  273.         // Provider for spherical harmonics
  274.         this.provider      = provider;
  275.         this.maxDegree     = provider.getMaxDegree();
  276.         this.maxOrder      = provider.getMaxOrder();

  278.         //set the maximum degree order for short periodics
  279.         checkIndexRange(maxDegreeTesseralSP, 2, maxDegree);
  280.         this.maxDegreeTesseralSP       = maxDegreeTesseralSP;

  282.         checkIndexRange(maxDegreeMdailyTesseralSP, 2, maxDegree);
  283.         this.maxDegreeMdailyTesseralSP = maxDegreeMdailyTesseralSP;

  285.         checkIndexRange(maxOrderTesseralSP, 0, maxOrder);
  286.         this.maxOrderTesseralSP        = maxOrderTesseralSP;

  288.         checkIndexRange(maxOrderMdailyTesseralSP, 0, maxOrder);
  289.         this.maxOrderMdailyTesseralSP  = maxOrderMdailyTesseralSP;

  291.         // set the maximum value for eccentricity power
  292.         if (maxOrder > 0) {
  293.             // Range check can be silently ignored if order = 0
  294.             checkIndexRange(maxEccPowTesseralSP, 0, maxOrder);
  295.         }
  296.         this.maxEccPowTesseralSP       = maxEccPowTesseralSP;

  298.         checkIndexRange(maxEccPowMdailyTesseralSP, 0, maxDegreeMdailyTesseralSP - 2);
  299.         this.maxEccPowMdailyTesseralSP = maxEccPowMdailyTesseralSP;

  301.         // set the maximum value for frequency
  302.         this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;

  304.         // Initialize default values
  305.         this.resOrders    = new ArrayList<Integer>();
  306.         this.nonResOrders = new TreeMap<Integer, List <Integer>>();

  308.         // Initialize default values
  309.         this.fieldShortPeriodTerms = new HashMap<>();
  310.         this.fieldHansen           = new HashMap<>();
  311.         this.maxEccPow             = 0;
  312.         this.maxHansen             = 0;

  314.     }

  316.     /** Check an index range.
  317.      * @param index index value
  318.      * @param min minimum value for index
  319.      * @param max maximum value for index
  320.      */
  321.     private void checkIndexRange(final int index, final int min, final int max) {
  322.         if (index < min || index > max) {
  323.             throw new OrekitException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, index, min, max);
  324.         }
  325.     }

  327.     /** {@inheritDoc} */
  328.     @Override
  329.     public List<ShortPeriodTerms> initializeShortPeriodTerms(final AuxiliaryElements auxiliaryElements,
  330.                                              final PropagationType type,
  331.                                              final double[] parameters) {

  333.         // Initializes specific parameters.

  335.         final DSSTTesseralContext context = initializeStep(auxiliaryElements, parameters);

  337.         // Set the highest power of the eccentricity in the analytical power
  338.         // series expansion for the averaged high order resonant central body
  339.         // spherical harmonic perturbation
  340.         maxEccPow = getMaxEccPow(auxiliaryElements.getEcc());

  342.         // Set the maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
  343.         maxHansen = maxEccPow / 2;

  345.         // The following terms are only used for hansen objects initialization
  346.         final double ratio = context.getRatio();

  348.         // Compute the non resonant tesseral harmonic terms if not set by the user
  349.         getResonantAndNonResonantTerms(type, context.getOrbitPeriod(), ratio);

  351.         hansen = new HansenObjects(ratio, type);

  353.         mMax = FastMath.max(maxOrderTesseralSP, maxOrderMdailyTesseralSP);

  355.         shortPeriodTerms = new TesseralShortPeriodicCoefficients(bodyFrame, maxOrderMdailyTesseralSP,
  356.                                                                  maxDegreeTesseralSP < 0, nonResOrders,
  357.                                                                  mMax, maxFrequencyShortPeriodics, INTERPOLATION_POINTS,
  358.                                                                  new TimeSpanMap<Slot>(new Slot(mMax, maxFrequencyShortPeriodics,
  359.                                                                                                 INTERPOLATION_POINTS)));

  361.         final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
  362.         list.add(shortPeriodTerms);
  363.         return list;

  365.     }

  367.     /** {@inheritDoc} */
  368.     @Override
  369.     public <T extends CalculusFieldElement<T>> List<FieldShortPeriodTerms<T>> initializeShortPeriodTerms(final FieldAuxiliaryElements<T> auxiliaryElements,
  370.                                                                                      final PropagationType type,
  371.                                                                                      final T[] parameters) {

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

  376.         // Initializes specific parameters.
  377.         final FieldDSSTTesseralContext<T> context = initializeStep(auxiliaryElements, parameters);

  379.         // Set the highest power of the eccentricity in the analytical power
  380.         // series expansion for the averaged high order resonant central body
  381.         // spherical harmonic perturbation
  382.         maxEccPow = getMaxEccPow(auxiliaryElements.getEcc().getReal());

  384.         // Set the maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
  385.         maxHansen = maxEccPow / 2;

  387.         // The following terms are only used for hansen objects initialization
  388.         final T ratio = context.getRatio();

  390.         // Compute the non resonant tesseral harmonic terms if not set by the user
  391.         // Field information is not important here
  392.         getResonantAndNonResonantTerms(type, context.getOrbitPeriod().getReal(), ratio.getReal());

  394.         mMax = FastMath.max(maxOrderTesseralSP, maxOrderMdailyTesseralSP);

  396.         fieldHansen.put(field, new FieldHansenObjects<>(ratio, type));

  398.         final FieldTesseralShortPeriodicCoefficients<T> ftspc =
  399.                         new FieldTesseralShortPeriodicCoefficients<>(bodyFrame, maxOrderMdailyTesseralSP,
  400.                                                                      maxDegreeTesseralSP < 0, nonResOrders,
  401.                                                                      mMax, maxFrequencyShortPeriodics, INTERPOLATION_POINTS,
  402.                                                                      new FieldTimeSpanMap<>(new FieldSlot<>(mMax,
  403.                                                                                                             maxFrequencyShortPeriodics,
  404.                                                                                                             INTERPOLATION_POINTS),
  405.                                                                                             field));

  407.         fieldShortPeriodTerms.put(field, ftspc);
  408.         return Collections.singletonList(ftspc);

  410.     }

  412.     /**
  413.      * Get the maximum power of the eccentricity to use in summation over s.
  414.      * @param e eccentricity
  415.      * @return the maximum power of the eccentricity
  416.      */
  417.     private int getMaxEccPow(final double e) {
  418.         // maxEccPow depends on satellite eccentricity
  419.         if (e <= 0.005) {
  420.             return 3;
  421.         } else if (e <= 0.02) {
  422.             return 4;
  423.         } else if (e <= 0.1) {
  424.             return 7;
  425.         } else if (e <= 0.2) {
  426.             return 10;
  427.         } else if (e <= 0.3) {
  428.             return 12;
  429.         } else if (e <= 0.4) {
  430.             return 15;
  431.         } else {
  432.             return 20;
  433.         }
  434.     }

  436.     /** Performs initialization at each integration step for the current force model.
  437.      *  <p>
  438.      *  This method aims at being called before mean elements rates computation.
  439.      *  </p>
  440.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  441.      *  @param parameters values of the force model parameters (only 1 value for each parameter)
  442.      *  that is to say that the extract parameter method {@link #extractParameters(double[], AbsoluteDate)}
  443.      *  should have be called before or the parameters list given in argument must correspond
  444.      *  to the extraction of parameter for a precise date {@link #getParameters(AbsoluteDate)}.
  445.      *  @return new force model context
  446.      */
  447.     private DSSTTesseralContext initializeStep(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
  448.         return new DSSTTesseralContext(auxiliaryElements, bodyFrame, provider, maxFrequencyShortPeriodics, bodyPeriod, parameters);
  449.     }

  451.     /** Performs initialization at each integration step for the current force model.
  452.      *  <p>
  453.      *  This method aims at being called before mean elements rates computation.
  454.      *  </p>
  455.      *  @param <T> type of the elements
  456.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  457.      *  @param parameters list of each estimated values for each driver of the force model parameters
  458.          *                (each span of each driver)
  459.      *  @return new force model context
  460.      */
  461.     private <T extends CalculusFieldElement<T>> FieldDSSTTesseralContext<T> initializeStep(final FieldAuxiliaryElements<T> auxiliaryElements,
  462.                                                                                        final T[] parameters) {
  463.         return new FieldDSSTTesseralContext<>(auxiliaryElements, bodyFrame, provider, maxFrequencyShortPeriodics, bodyPeriod, parameters);
  464.     }

  466.     /** {@inheritDoc} */
  467.     @Override
  468.     public double[] getMeanElementRate(final SpacecraftState spacecraftState,
  469.                                        final AuxiliaryElements auxiliaryElements, final double[] parameters) {

  471.         // Container for attributes

  473.         final DSSTTesseralContext context = initializeStep(auxiliaryElements, parameters);

  475.         // Access to potential U derivatives
  476.         final UAnddU udu = new UAnddU(spacecraftState.getDate(), context, hansen);

  478.         // Compute the cross derivative operator :
  479.         final double UAlphaGamma   = auxiliaryElements.getAlpha() * udu.getdUdGa() - auxiliaryElements.getGamma() * udu.getdUdAl();
  480.         final double UAlphaBeta    = auxiliaryElements.getAlpha() * udu.getdUdBe() - auxiliaryElements.getBeta()  * udu.getdUdAl();
  481.         final double UBetaGamma    = auxiliaryElements.getBeta() * udu.getdUdGa() - auxiliaryElements.getGamma() * udu.getdUdBe();
  482.         final double Uhk           = auxiliaryElements.getH() * udu.getdUdk()  - auxiliaryElements.getK() * udu.getdUdh();
  483.         final double pUagmIqUbgoAB = (auxiliaryElements.getP() * UAlphaGamma - I * auxiliaryElements.getQ() * UBetaGamma) * context.getOoAB();
  484.         final double UhkmUabmdUdl  = Uhk - UAlphaBeta - udu.getdUdl();

  486.         final double da =  context.getAx2oA() * udu.getdUdl();
  487.         final double dh =  context.getBoA() * udu.getdUdk() + auxiliaryElements.getK() * pUagmIqUbgoAB - auxiliaryElements.getH() * context.getBoABpo() * udu.getdUdl();
  488.         final double dk =  -(context.getBoA() * udu.getdUdh() + auxiliaryElements.getH() * pUagmIqUbgoAB + auxiliaryElements.getK() * context.getBoABpo() * udu.getdUdl());
  489.         final double dp =  context.getCo2AB() * (auxiliaryElements.getP() * UhkmUabmdUdl - UBetaGamma);
  490.         final double dq =  context.getCo2AB() * (auxiliaryElements.getQ() * UhkmUabmdUdl - I * UAlphaGamma);
  491.         final double dM = -context.getAx2oA() * udu.getdUda() + context.getBoABpo() * (auxiliaryElements.getH() * udu.getdUdh() + auxiliaryElements.getK() * udu.getdUdk()) + pUagmIqUbgoAB;

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

  496.     /** {@inheritDoc} */
  497.     @Override
  498.     public <T extends CalculusFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> spacecraftState,
  499.                                                                   final FieldAuxiliaryElements<T> auxiliaryElements,
  500.                                                                   final T[] parameters) {

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

  505.         // Container for attributes

  507.         final FieldDSSTTesseralContext<T> context = initializeStep(auxiliaryElements, parameters);

  509.         @SuppressWarnings("unchecked")
  510.         final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);
  511.         // Access to potential U derivatives
  512.         final FieldUAnddU<T> udu = new FieldUAnddU<>(spacecraftState.getDate(), context, fho);

  514.         // Compute the cross derivative operator :
  515.         final T UAlphaGamma   = udu.getdUdGa().multiply(auxiliaryElements.getAlpha()).subtract(udu.getdUdAl().multiply(auxiliaryElements.getGamma()));
  516.         final T UAlphaBeta    = udu.getdUdBe().multiply(auxiliaryElements.getAlpha()).subtract(udu.getdUdAl().multiply(auxiliaryElements.getBeta()));
  517.         final T UBetaGamma    = udu.getdUdGa().multiply(auxiliaryElements.getBeta()).subtract(udu.getdUdBe().multiply(auxiliaryElements.getGamma()));
  518.         final T Uhk           = udu.getdUdk().multiply(auxiliaryElements.getH()).subtract(udu.getdUdh().multiply(auxiliaryElements.getK()));
  519.         final T pUagmIqUbgoAB = (UAlphaGamma.multiply(auxiliaryElements.getP()).subtract(UBetaGamma.multiply(auxiliaryElements.getQ()).multiply(I))).multiply(context.getOoAB());
  520.         final T UhkmUabmdUdl  = Uhk.subtract(UAlphaBeta).subtract(udu.getdUdl());

  522.         final T da = udu.getdUdl().multiply(context.getAx2oA());
  523.         final T dh = udu.getdUdk().multiply(context.getBoA()).add(pUagmIqUbgoAB.multiply(auxiliaryElements.getK())).subtract(udu.getdUdl().multiply(auxiliaryElements.getH()).multiply(context.getBoABpo()));
  524.         final T dk = (udu.getdUdh().multiply(context.getBoA()).add(pUagmIqUbgoAB.multiply(auxiliaryElements.getH())).add(udu.getdUdl().multiply(context.getBoABpo()).multiply(auxiliaryElements.getK()))).negate();
  525.         final T dp = context.getCo2AB().multiply(auxiliaryElements.getP().multiply(UhkmUabmdUdl).subtract(UBetaGamma));
  526.         final T dq = context.getCo2AB().multiply(auxiliaryElements.getQ().multiply(UhkmUabmdUdl).subtract(UAlphaGamma.multiply(I)));
  527.         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()));

  529.         final T[] elements = MathArrays.buildArray(field, 6);
  530.         elements[0] = da;
  531.         elements[1] = dk;
  532.         elements[2] = dh;
  533.         elements[3] = dq;
  534.         elements[4] = dp;
  535.         elements[5] = dM;

  537.         return elements;

  539.     }

  541.     /** {@inheritDoc} */
  542.     @Override
  543.     public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {

  545.         final Slot slot = shortPeriodTerms.createSlot(meanStates);

  547.         for (final SpacecraftState meanState : meanStates) {

  549.             final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(meanState.getOrbit(), I);

  551.             // Extract the proper parameters valid at date from the input array
  552.             final double[] extractedParameters = this.extractParameters(parameters, auxiliaryElements.getDate());
  553.             final DSSTTesseralContext context = initializeStep(auxiliaryElements, extractedParameters);

  555.             // Initialise the Hansen coefficients
  556.             for (int s = -maxDegree; s <= maxDegree; s++) {
  557.                 // coefficients with j == 0 are always needed
  558.                 hansen.computeHansenObjectsInitValues(context, s + maxDegree, 0);
  559.                 if (maxDegreeTesseralSP >= 0) {
  560.                     // initialize other objects only if required
  561.                     for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  562.                         hansen.computeHansenObjectsInitValues(context, s + maxDegree, j);
  563.                     }
  564.                 }
  565.             }

  567.             final FourierCjSjCoefficients cjsjFourier = new FourierCjSjCoefficients(maxFrequencyShortPeriodics, mMax);

  569.             // Compute coefficients
  570.             // Compute only if there is at least one non-resonant tesseral
  571.             if (!nonResOrders.isEmpty() || maxDegreeTesseralSP < 0) {
  572.                 // Generate the fourrier coefficients
  573.                 cjsjFourier.generateCoefficients(meanState.getDate(), context, hansen);

  575.                 // the coefficient 3n / 2a
  576.                 final double tnota = 1.5 * context.getMeanMotion() / auxiliaryElements.getSma();

  578.                 // build the mDaily coefficients
  579.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  580.                     // build the coefficients
  581.                     buildCoefficients(cjsjFourier, meanState.getDate(), slot, m, 0, tnota, context);
  582.                 }

  584.                 if (maxDegreeTesseralSP >= 0) {
  585.                     // generate the other coefficients, if required
  586.                     for (final Map.Entry<Integer, List<Integer>> entry : nonResOrders.entrySet()) {

  588.                         for (int j : entry.getValue()) {
  589.                             // build the coefficients
  590.                             buildCoefficients(cjsjFourier, meanState.getDate(), slot, entry.getKey(), j, tnota, context);
  591.                         }
  592.                     }
  593.                 }
  594.             }

  596.         }

  598.     }

  600.     /** {@inheritDoc} */
  601.     @Override
  602.     @SuppressWarnings("unchecked")
  603.     public <T extends CalculusFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
  604.                                                                        final FieldSpacecraftState<T>... meanStates) {

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

  609.         final FieldTesseralShortPeriodicCoefficients<T> ftspc = (FieldTesseralShortPeriodicCoefficients<T>) fieldShortPeriodTerms.get(field);
  610.         final FieldSlot<T> slot = ftspc.createSlot(meanStates);

  612.         for (final FieldSpacecraftState<T> meanState : meanStates) {

  614.             final FieldAuxiliaryElements<T> auxiliaryElements = new FieldAuxiliaryElements<>(meanState.getOrbit(), I);

  616.             // Extract the proper parameters valid at date from the input array
  617.             final T[] extractedParameters = this.extractParameters(parameters, auxiliaryElements.getDate());
  618.             final FieldDSSTTesseralContext<T> context = initializeStep(auxiliaryElements, extractedParameters);

  620.             final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);
  621.             // Initialise the Hansen coefficients
  622.             for (int s = -maxDegree; s <= maxDegree; s++) {
  623.                 // coefficients with j == 0 are always needed
  624.                 fho.computeHansenObjectsInitValues(context, s + maxDegree, 0);
  625.                 if (maxDegreeTesseralSP >= 0) {
  626.                     // initialize other objects only if required
  627.                     for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  628.                         fho.computeHansenObjectsInitValues(context, s + maxDegree, j);
  629.                     }
  630.                 }
  631.             }

  633.             final FieldFourierCjSjCoefficients<T> cjsjFourier = new FieldFourierCjSjCoefficients<>(maxFrequencyShortPeriodics, mMax, field);

  635.             // Compute coefficients
  636.             // Compute only if there is at least one non-resonant tesseral
  637.             if (!nonResOrders.isEmpty() || maxDegreeTesseralSP < 0) {
  638.                 // Generate the fourrier coefficients
  639.                 cjsjFourier.generateCoefficients(meanState.getDate(), context, fho, field);

  641.                 // the coefficient 3n / 2a
  642.                 final T tnota = context.getMeanMotion().multiply(1.5).divide(auxiliaryElements.getSma());

  644.                 // build the mDaily coefficients
  645.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  646.                     // build the coefficients
  647.                     buildCoefficients(cjsjFourier, meanState.getDate(), slot, m, 0, tnota, context, field);
  648.                 }

  650.                 if (maxDegreeTesseralSP >= 0) {
  651.                     // generate the other coefficients, if required
  652.                     for (final Map.Entry<Integer, List<Integer>> entry : nonResOrders.entrySet()) {

  654.                         for (int j : entry.getValue()) {
  655.                             // build the coefficients
  656.                             buildCoefficients(cjsjFourier, meanState.getDate(), slot, entry.getKey(), j, tnota, context, field);
  657.                         }
  658.                     }
  659.                 }
  660.             }

  662.         }

  664.     }

  666.     /** {@inheritDoc} */
  667.     public List<ParameterDriver> getParametersDrivers() {
  668.         return Collections.singletonList(gmParameterDriver);
  669.     }

  671.     /** Build a set of coefficients.
  672.      * @param cjsjFourier the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  673.      * @param date the current date
  674.      * @param slot slot to which the coefficients belong
  675.      * @param m m index
  676.      * @param j j index
  677.      * @param tnota 3n/2a
  678.      * @param context container for attributes
  679.      */
  680.     private void buildCoefficients(final FourierCjSjCoefficients cjsjFourier,
  681.                                    final AbsoluteDate date, final Slot slot,
  682.                                    final int m, final int j, final double tnota, final DSSTTesseralContext context) {

  684.         // Create local arrays
  685.         final double[] currentCijm = new double[] {0., 0., 0., 0., 0., 0.};
  686.         final double[] currentSijm = new double[] {0., 0., 0., 0., 0., 0.};

  688.         // compute the term 1 / (jn - mθ<sup>.</sup>)
  689.         final double oojnmt = 1. / (j * context.getMeanMotion() - m * centralBodyRotationRate);

  691.         // initialise the coeficients
  692.         for (int i = 0; i < 6; i++) {
  693.             currentCijm[i] = -cjsjFourier.getSijm(i, j, m);
  694.             currentSijm[i] = cjsjFourier.getCijm(i, j, m);
  695.         }
  696.         // Add the separate part for δ<sub>6i</sub>
  697.         currentCijm[5] += tnota * oojnmt * cjsjFourier.getCijm(0, j, m);
  698.         currentSijm[5] += tnota * oojnmt * cjsjFourier.getSijm(0, j, m);

  700.         //Multiply by 1 / (jn - mθ<sup>.</sup>)
  701.         for (int i = 0; i < 6; i++) {
  702.             currentCijm[i] *= oojnmt;
  703.             currentSijm[i] *= oojnmt;
  704.         }

  706.         // Add the coefficients to the interpolation grid
  707.         slot.cijm[m][j + maxFrequencyShortPeriodics].addGridPoint(date, currentCijm);
  708.         slot.sijm[m][j + maxFrequencyShortPeriodics].addGridPoint(date, currentSijm);

  710.     }

  712.      /** Build a set of coefficients.
  713.      * @param <T> the type of the field elements
  714.      * @param cjsjFourier the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  715.      * @param date the current date
  716.      * @param slot slot to which the coefficients belong
  717.      * @param m m index
  718.      * @param j j index
  719.      * @param tnota 3n/2a
  720.      * @param context container for attributes
  721.      * @param field field used by default
  722.      */
  723.     private <T extends CalculusFieldElement<T>> void buildCoefficients(final FieldFourierCjSjCoefficients<T> cjsjFourier,
  724.                                                                    final FieldAbsoluteDate<T> date,
  725.                                                                    final FieldSlot<T> slot,
  726.                                                                    final int m, final int j, final T tnota,
  727.                                                                    final FieldDSSTTesseralContext<T> context,
  728.                                                                    final Field<T> field) {

  730.         // Zero
  731.         final T zero = field.getZero();

  733.         // Create local arrays
  734.         final T[] currentCijm = MathArrays.buildArray(field, 6);
  735.         final T[] currentSijm = MathArrays.buildArray(field, 6);

  737.         Arrays.fill(currentCijm, zero);
  738.         Arrays.fill(currentSijm, zero);

  740.         // compute the term 1 / (jn - mθ<sup>.</sup>)
  741.         final T oojnmt = (context.getMeanMotion().multiply(j).subtract(m * centralBodyRotationRate)).reciprocal();

  743.         // initialise the coeficients
  744.         for (int i = 0; i < 6; i++) {
  745.             currentCijm[i] = cjsjFourier.getSijm(i, j, m).negate();
  746.             currentSijm[i] = cjsjFourier.getCijm(i, j, m);
  747.         }
  748.         // Add the separate part for δ<sub>6i</sub>
  749.         currentCijm[5] = currentCijm[5].add(tnota.multiply(oojnmt).multiply(cjsjFourier.getCijm(0, j, m)));
  750.         currentSijm[5] = currentSijm[5].add(tnota.multiply(oojnmt).multiply(cjsjFourier.getSijm(0, j, m)));

  752.         //Multiply by 1 / (jn - mθ<sup>.</sup>)
  753.         for (int i = 0; i < 6; i++) {
  754.             currentCijm[i] = currentCijm[i].multiply(oojnmt);
  755.             currentSijm[i] = currentSijm[i].multiply(oojnmt);
  756.         }

  758.         // Add the coefficients to the interpolation grid
  759.         slot.cijm[m][j + maxFrequencyShortPeriodics].addGridPoint(date, currentCijm);
  760.         slot.sijm[m][j + maxFrequencyShortPeriodics].addGridPoint(date, currentSijm);

  762.     }

  764.      /**
  765.       * Get the resonant and non-resonant tesseral terms in the central body spherical harmonic field.
  766.       *
  767.       * @param type type of the elements used during the propagation
  768.       * @param orbitPeriod Keplerian period
  769.       * @param ratio ratio of satellite period to central body rotation period
  770.       */
  771.     private void getResonantAndNonResonantTerms(final PropagationType type, final double orbitPeriod,
  772.                                                 final double ratio) {

  774.         // Compute natural resonant terms
  775.         final double tolerance = 1. / FastMath.max(MIN_PERIOD_IN_SAT_REV,
  776.                                                    MIN_PERIOD_IN_SECONDS / orbitPeriod);

  778.         // Search the resonant orders in the tesseral harmonic field
  779.         resOrders.clear();
  780.         nonResOrders.clear();
  781.         for (int m = 1; m <= maxOrder; m++) {
  782.             final double resonance = ratio * m;
  783.             int jRes = 0;
  784.             final int jComputedRes = (int) FastMath.round(resonance);
  785.             if (jComputedRes > 0 && jComputedRes <= maxFrequencyShortPeriodics && FastMath.abs(resonance - jComputedRes) <= tolerance) {
  786.                 // Store each resonant index and order
  787.                 resOrders.add(m);
  788.                 jRes = jComputedRes;
  789.             }

  791.             if (type == PropagationType.OSCULATING && maxDegreeTesseralSP >= 0 && m <= maxOrderTesseralSP) {
  792.                 //compute non resonant orders in the tesseral harmonic field
  793.                 final List<Integer> listJofM = new ArrayList<Integer>();
  794.                 //for the moment we take only the pairs (j,m) with |j| <= maxDegree + maxEccPow (from |s-j| <= maxEccPow and |s| <= maxDegree)
  795.                 for (int j = -maxFrequencyShortPeriodics; j <= maxFrequencyShortPeriodics; j++) {
  796.                     if (j != 0 && j != jRes) {
  797.                         listJofM.add(j);
  798.                     }
  799.                 }

  801.                 nonResOrders.put(m, listJofM);
  802.             }
  803.         }
  804.     }

  806.     /** Compute the n-SUM for potential derivatives components.
  807.      *  @param date current date
  808.      *  @param j resonant index <i>j</i>
  809.      *  @param m resonant order <i>m</i>
  810.      *  @param s d'Alembert characteristic <i>s</i>
  811.      *  @param maxN maximum possible value for <i>n</i> index
  812.      *  @param roaPow powers of R/a up to degree <i>n</i>
  813.      *  @param ghMSJ G<sup>j</sup><sub>m,s</sub> and H<sup>j</sup><sub>m,s</sub> polynomials
  814.      *  @param gammaMNS &Gamma;<sup>m</sup><sub>n,s</sub>(γ) function
  815.      *  @param context container for attributes
  816.      *  @param hansenObjects initialization of hansen objects
  817.      *  @return Components of U<sub>n</sub> derivatives for fixed j, m, s
  818.      */
  819.     private double[][] computeNSum(final AbsoluteDate date,
  820.                                    final int j, final int m, final int s, final int maxN, final double[] roaPow,
  821.                                    final GHmsjPolynomials ghMSJ, final GammaMnsFunction gammaMNS, final DSSTTesseralContext context,
  822.                                    final HansenObjects hansenObjects) {

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

  827.         //spherical harmonics
  828.         final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);

  830.         // Potential derivatives components
  831.         double dUdaCos  = 0.;
  832.         double dUdaSin  = 0.;
  833.         double dUdhCos  = 0.;
  834.         double dUdhSin  = 0.;
  835.         double dUdkCos  = 0.;
  836.         double dUdkSin  = 0.;
  837.         double dUdlCos  = 0.;
  838.         double dUdlSin  = 0.;
  839.         double dUdAlCos = 0.;
  840.         double dUdAlSin = 0.;
  841.         double dUdBeCos = 0.;
  842.         double dUdBeSin = 0.;
  843.         double dUdGaCos = 0.;
  844.         double dUdGaSin = 0.;

  846.         // I^m
  847.         @SuppressWarnings("unused")
  848.         final int Im = I > 0 ? 1 : (m % 2 == 0 ? 1 : -1);

  850.         // jacobi v, w, indices from 2.7.1-(15)
  851.         final int v = FastMath.abs(m - s);
  852.         final int w = FastMath.abs(m + s);

  854.         // Initialise lower degree nmin = (Max(2, m, |s|)) for summation over n
  855.         final int nmin = FastMath.max(FastMath.max(2, m), FastMath.abs(s));

  857.         //Get the corresponding Hansen object
  858.         final int sIndex = maxDegree + (j < 0 ? -s : s);
  859.         final int jIndex = FastMath.abs(j);
  860.         final HansenTesseralLinear hans = hansenObjects.getHansenObjects()[sIndex][jIndex];

  862.         // n-SUM from nmin to N
  863.         for (int n = nmin; n <= maxN; n++) {
  864.             // If (n - s) is odd, the contribution is null because of Vmns
  865.             if ((n - s) % 2 == 0) {

  867.                 // Vmns coefficient
  868.                 final double vMNS   = CoefficientsFactory.getVmns(m, n, s);

  870.                 // Inclination function Gamma and derivative
  871.                 final double gaMNS  = gammaMNS.getValue(m, n, s);
  872.                 final double dGaMNS = gammaMNS.getDerivative(m, n, s);

  874.                 // Hansen kernel value and derivative
  875.                 final double kJNS   = hans.getValue(-n - 1, context.getChi());
  876.                 final double dkJNS  = hans.getDerivative(-n - 1, context.getChi());

  878.                 // Gjms, Hjms polynomials and derivatives
  879.                 final double gMSJ   = ghMSJ.getGmsj(m, s, j);
  880.                 final double hMSJ   = ghMSJ.getHmsj(m, s, j);
  881.                 final double dGdh   = ghMSJ.getdGmsdh(m, s, j);
  882.                 final double dGdk   = ghMSJ.getdGmsdk(m, s, j);
  883.                 final double dGdA   = ghMSJ.getdGmsdAlpha(m, s, j);
  884.                 final double dGdB   = ghMSJ.getdGmsdBeta(m, s, j);
  885.                 final double dHdh   = ghMSJ.getdHmsdh(m, s, j);
  886.                 final double dHdk   = ghMSJ.getdHmsdk(m, s, j);
  887.                 final double dHdA   = ghMSJ.getdHmsdAlpha(m, s, j);
  888.                 final double dHdB   = ghMSJ.getdHmsdBeta(m, s, j);

  890.                 // Jacobi l-index from 2.7.1-(15)
  891.                 final int l = FastMath.min(n - m, n - FastMath.abs(s));
  892.                 // Jacobi polynomial and derivative
  893.                 final double[] jacobi = JacobiPolynomials.getValueAndDerivative(l, v, w, auxiliaryElements.getGamma());

  895.                 // Geopotential coefficients
  896.                 final double cnm = harmonics.getUnnormalizedCnm(n, m);
  897.                 final double snm = harmonics.getUnnormalizedSnm(n, m);

  899.                 // Common factors from expansion of equations 3.3-4
  900.                 final double cf_0      = roaPow[n] * Im * vMNS;
  901.                 final double cf_1      = cf_0 * gaMNS * jacobi[0]; //jacobi.getValue();
  902.                 final double cf_2      = cf_1 * kJNS;
  903.                 final double gcPhs     = gMSJ * cnm + hMSJ * snm;
  904.                 final double gsMhc     = gMSJ * snm - hMSJ * cnm;
  905.                 final double dKgcPhsx2 = 2. * dkJNS * gcPhs;
  906.                 final double dKgsMhcx2 = 2. * dkJNS * gsMhc;
  907.                 final double dUdaCoef  = (n + 1) * cf_2;
  908.                 final double dUdlCoef  = j * cf_2;
  909.                 //final double dUdGaCoef = cf_0 * kJNS * (jacobi.getValue() * dGaMNS + gaMNS * jacobi.getGradient()[0]);
  910.                 final double dUdGaCoef = cf_0 * kJNS * (jacobi[0] * dGaMNS + gaMNS * jacobi[1]);

  912.                 // dU / da components
  913.                 dUdaCos  += dUdaCoef * gcPhs;
  914.                 dUdaSin  += dUdaCoef * gsMhc;

  916.                 // dU / dh components
  917.                 dUdhCos  += cf_1 * (kJNS * (cnm * dGdh + snm * dHdh) + auxiliaryElements.getH() * dKgcPhsx2);
  918.                 dUdhSin  += cf_1 * (kJNS * (snm * dGdh - cnm * dHdh) + auxiliaryElements.getH() * dKgsMhcx2);

  920.                 // dU / dk components
  921.                 dUdkCos  += cf_1 * (kJNS * (cnm * dGdk + snm * dHdk) + auxiliaryElements.getK() * dKgcPhsx2);
  922.                 dUdkSin  += cf_1 * (kJNS * (snm * dGdk - cnm * dHdk) + auxiliaryElements.getK() * dKgsMhcx2);

  924.                 // dU / dLambda components
  925.                 dUdlCos  +=  dUdlCoef * gsMhc;
  926.                 dUdlSin  += -dUdlCoef * gcPhs;

  928.                 // dU / alpha components
  929.                 dUdAlCos += cf_2 * (dGdA * cnm + dHdA * snm);
  930.                 dUdAlSin += cf_2 * (dGdA * snm - dHdA * cnm);

  932.                 // dU / dBeta components
  933.                 dUdBeCos += cf_2 * (dGdB * cnm + dHdB * snm);
  934.                 dUdBeSin += cf_2 * (dGdB * snm - dHdB * cnm);

  936.                 // dU / dGamma components
  937.                 dUdGaCos += dUdGaCoef * gcPhs;
  938.                 dUdGaSin += dUdGaCoef * gsMhc;
  939.             }
  940.         }

  942.         return new double[][] { { dUdaCos,  dUdaSin  },
  943.                                 { dUdhCos,  dUdhSin  },
  944.                                 { dUdkCos,  dUdkSin  },
  945.                                 { dUdlCos,  dUdlSin  },
  946.                                 { dUdAlCos, dUdAlSin },
  947.                                 { dUdBeCos, dUdBeSin },
  948.                                 { dUdGaCos, dUdGaSin } };
  949.     }

  951.     /** Compute the n-SUM for potential derivatives components.
  952.      *  @param <T> the type of the field elements
  953.      *  @param date current date
  954.      *  @param j resonant index <i>j</i>
  955.      *  @param m resonant order <i>m</i>
  956.      *  @param s d'Alembert characteristic <i>s</i>
  957.      *  @param maxN maximum possible value for <i>n</i> index
  958.      *  @param roaPow powers of R/a up to degree <i>n</i>
  959.      *  @param ghMSJ G<sup>j</sup><sub>m,s</sub> and H<sup>j</sup><sub>m,s</sub> polynomials
  960.      *  @param gammaMNS &Gamma;<sup>m</sup><sub>n,s</sub>(γ) function
  961.      *  @param context container for attributes
  962.      *  @param hansenObjects initialization of hansen objects
  963.      *  @return Components of U<sub>n</sub> derivatives for fixed j, m, s
  964.      */
  965.     private <T extends CalculusFieldElement<T>> T[][] computeNSum(final FieldAbsoluteDate<T> date,
  966.                                                               final int j, final int m, final int s, final int maxN,
  967.                                                               final T[] roaPow,
  968.                                                               final FieldGHmsjPolynomials<T> ghMSJ,
  969.                                                               final FieldGammaMnsFunction<T> gammaMNS,
  970.                                                               final FieldDSSTTesseralContext<T> context,
  971.                                                               final FieldHansenObjects<T> hansenObjects) {

  973.         // Auxiliary elements related to the current orbit
  974.         final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
  975.         // Zero for initialization
  976.         final Field<T> field = date.getField();
  977.         final T zero = field.getZero();

  979.         //spherical harmonics
  980.         final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date.toAbsoluteDate());

  982.         // Potential derivatives components
  983.         T dUdaCos  = zero;
  984.         T dUdaSin  = zero;
  985.         T dUdhCos  = zero;
  986.         T dUdhSin  = zero;
  987.         T dUdkCos  = zero;
  988.         T dUdkSin  = zero;
  989.         T dUdlCos  = zero;
  990.         T dUdlSin  = zero;
  991.         T dUdAlCos = zero;
  992.         T dUdAlSin = zero;
  993.         T dUdBeCos = zero;
  994.         T dUdBeSin = zero;
  995.         T dUdGaCos = zero;
  996.         T dUdGaSin = zero;

  998.         // I^m
  999.         @SuppressWarnings("unused")
  1000.         final int Im = I > 0 ? 1 : (m % 2 == 0 ? 1 : -1);

  1002.         // jacobi v, w, indices from 2.7.1-(15)
  1003.         final int v = FastMath.abs(m - s);
  1004.         final int w = FastMath.abs(m + s);

  1006.         // Initialise lower degree nmin = (Max(2, m, |s|)) for summation over n
  1007.         final int nmin = FastMath.max(FastMath.max(2, m), FastMath.abs(s));

  1009.         //Get the corresponding Hansen object
  1010.         final int sIndex = maxDegree + (j < 0 ? -s : s);
  1011.         final int jIndex = FastMath.abs(j);
  1012.         final FieldHansenTesseralLinear<T> hans = hansenObjects.getHansenObjects()[sIndex][jIndex];

  1014.         // n-SUM from nmin to N
  1015.         for (int n = nmin; n <= maxN; n++) {
  1016.             // If (n - s) is odd, the contribution is null because of Vmns
  1017.             if ((n - s) % 2 == 0) {

  1019.                 // Vmns coefficient
  1020.                 final T vMNS   = zero.add(CoefficientsFactory.getVmns(m, n, s));

  1022.                 // Inclination function Gamma and derivative
  1023.                 final T gaMNS  = gammaMNS.getValue(m, n, s);
  1024.                 final T dGaMNS = gammaMNS.getDerivative(m, n, s);

  1026.                 // Hansen kernel value and derivative
  1027.                 final T kJNS   = hans.getValue(-n - 1, context.getChi());
  1028.                 final T dkJNS  = hans.getDerivative(-n - 1, context.getChi());

  1030.                 // Gjms, Hjms polynomials and derivatives
  1031.                 final T gMSJ   = ghMSJ.getGmsj(m, s, j);
  1032.                 final T hMSJ   = ghMSJ.getHmsj(m, s, j);
  1033.                 final T dGdh   = ghMSJ.getdGmsdh(m, s, j);
  1034.                 final T dGdk   = ghMSJ.getdGmsdk(m, s, j);
  1035.                 final T dGdA   = ghMSJ.getdGmsdAlpha(m, s, j);
  1036.                 final T dGdB   = ghMSJ.getdGmsdBeta(m, s, j);
  1037.                 final T dHdh   = ghMSJ.getdHmsdh(m, s, j);
  1038.                 final T dHdk   = ghMSJ.getdHmsdk(m, s, j);
  1039.                 final T dHdA   = ghMSJ.getdHmsdAlpha(m, s, j);
  1040.                 final T dHdB   = ghMSJ.getdHmsdBeta(m, s, j);

  1042.                 // Jacobi l-index from 2.7.1-(15)
  1043.                 final int l = FastMath.min(n - m, n - FastMath.abs(s));
  1044.                 // Jacobi polynomial and derivative
  1045.                 final FieldGradient<T> jacobi =
  1046.                         JacobiPolynomials.getValue(l, v, w, FieldGradient.variable(1, 0, auxiliaryElements.getGamma()));

  1048.                 // Geopotential coefficients
  1049.                 final T cnm = zero.add(harmonics.getUnnormalizedCnm(n, m));
  1050.                 final T snm = zero.add(harmonics.getUnnormalizedSnm(n, m));

  1052.                 // Common factors from expansion of equations 3.3-4
  1053.                 final T cf_0      = roaPow[n].multiply(Im).multiply(vMNS);
  1054.                 final T cf_1      = cf_0.multiply(gaMNS).multiply(jacobi.getValue());
  1055.                 final T cf_2      = cf_1.multiply(kJNS);
  1056.                 final T gcPhs     = gMSJ.multiply(cnm).add(hMSJ.multiply(snm));
  1057.                 final T gsMhc     = gMSJ.multiply(snm).subtract(hMSJ.multiply(cnm));
  1058.                 final T dKgcPhsx2 = dkJNS.multiply(gcPhs).multiply(2.);
  1059.                 final T dKgsMhcx2 = dkJNS.multiply(gsMhc).multiply(2.);
  1060.                 final T dUdaCoef  = cf_2.multiply(n + 1);
  1061.                 final T dUdlCoef  = cf_2.multiply(j);
  1062.                 final T dUdGaCoef = cf_0.multiply(kJNS).multiply(dGaMNS.multiply(jacobi.getValue()).add(gaMNS.multiply(jacobi.getGradient()[0])));

  1064.                 // dU / da components
  1065.                 dUdaCos  = dUdaCos.add(dUdaCoef.multiply(gcPhs));
  1066.                 dUdaSin  = dUdaSin.add(dUdaCoef.multiply(gsMhc));

  1068.                 // dU / dh components
  1069.                 dUdhCos  = dUdhCos.add(cf_1.multiply(kJNS.multiply(cnm.multiply(dGdh).add(snm.multiply(dHdh))).add(dKgcPhsx2.multiply(auxiliaryElements.getH()))));
  1070.                 dUdhSin  = dUdhSin.add(cf_1.multiply(kJNS.multiply(snm.multiply(dGdh).subtract(cnm.multiply(dHdh))).add(dKgsMhcx2.multiply(auxiliaryElements.getH()))));

  1072.                 // dU / dk components
  1073.                 dUdkCos  = dUdkCos.add(cf_1.multiply(kJNS.multiply(cnm.multiply(dGdk).add(snm.multiply(dHdk))).add(dKgcPhsx2.multiply(auxiliaryElements.getK()))));
  1074.                 dUdkSin  = dUdkSin.add(cf_1.multiply(kJNS.multiply(snm.multiply(dGdk).subtract(cnm.multiply(dHdk))).add(dKgsMhcx2.multiply(auxiliaryElements.getK()))));

  1076.                 // dU / dLambda components
  1077.                 dUdlCos  = dUdlCos.add(dUdlCoef.multiply(gsMhc));
  1078.                 dUdlSin  = dUdlSin.add(dUdlCoef.multiply(gcPhs).negate());

  1080.                 // dU / alpha components
  1081.                 dUdAlCos = dUdAlCos.add(cf_2.multiply(dGdA.multiply(cnm).add(dHdA.multiply(snm))));
  1082.                 dUdAlSin = dUdAlSin.add(cf_2.multiply(dGdA.multiply(snm).subtract(dHdA.multiply(cnm))));

  1084.                 // dU / dBeta components
  1085.                 dUdBeCos = dUdBeCos.add(cf_2.multiply(dGdB.multiply(cnm).add(dHdB.multiply(snm))));
  1086.                 dUdBeSin = dUdBeSin.add(cf_2.multiply(dGdB.multiply(snm).subtract(dHdB.multiply(cnm))));

  1088.                 // dU / dGamma components
  1089.                 dUdGaCos = dUdGaCos.add(dUdGaCoef.multiply(gcPhs));
  1090.                 dUdGaSin = dUdGaSin.add(dUdGaCoef.multiply(gsMhc));
  1091.             }
  1092.         }

  1094.         final T[][] derivatives = MathArrays.buildArray(field, 7, 2);
  1095.         derivatives[0][0] = dUdaCos;
  1096.         derivatives[0][1] = dUdaSin;
  1097.         derivatives[1][0] = dUdhCos;
  1098.         derivatives[1][1] = dUdhSin;
  1099.         derivatives[2][0] = dUdkCos;
  1100.         derivatives[2][1] = dUdkSin;
  1101.         derivatives[3][0] = dUdlCos;
  1102.         derivatives[3][1] = dUdlSin;
  1103.         derivatives[4][0] = dUdAlCos;
  1104.         derivatives[4][1] = dUdAlSin;
  1105.         derivatives[5][0] = dUdBeCos;
  1106.         derivatives[5][1] = dUdBeSin;
  1107.         derivatives[6][0] = dUdGaCos;
  1108.         derivatives[6][1] = dUdGaSin;

  1110.         return derivatives;

  1112.     }

  1114.     /** {@inheritDoc} */
  1115.     @Override
  1116.     public void registerAttitudeProvider(final AttitudeProvider attitudeProvider) {
  1117.         //nothing is done since this contribution is not sensitive to attitude
  1118.     }

  1120.     /** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
  1121.      *  <p>
  1122.      *  Those coefficients are given in Danielson paper by substituting the
  1123.      *  disturbing function (2.7.1-16) with m != 0 into (2.2-10)
  1124.      *  </p>
  1125.      */
  1126.     private class FourierCjSjCoefficients {

  1128.         /** Absolute limit for j ( -jMax <= j <= jMax ).  */
  1129.         private final int jMax;

  1131.         /** The C<sub>i</sub><sup>jm</sup> coefficients.
  1132.          * <p>
  1133.          * The index order is [m][j][i] <br/>
  1134.          * The i index corresponds to the C<sub>i</sub><sup>jm</sup> coefficients used to
  1135.          * compute the following: <br/>
  1136.          * - da/dt <br/>
  1137.          * - dk/dt <br/>
  1138.          * - dh/dt / dk <br/>
  1139.          * - dq/dt <br/>
  1140.          * - dp/dt / dα <br/>
  1141.          * - dλ/dt / dβ <br/>
  1142.          * </p>
  1143.          */
  1144.         private final double[][][] cCoef;

  1146.         /** The S<sub>i</sub><sup>jm</sup> coefficients.
  1147.          * <p>
  1148.          * The index order is [m][j][i] <br/>
  1149.          * The i index corresponds to the C<sub>i</sub><sup>jm</sup> coefficients used to
  1150.          * compute the following: <br/>
  1151.          * - da/dt <br/>
  1152.          * - dk/dt <br/>
  1153.          * - dh/dt / dk <br/>
  1154.          * - dq/dt <br/>
  1155.          * - dp/dt / dα <br/>
  1156.          * - dλ/dt / dβ <br/>
  1157.          * </p>
  1158.          */
  1159.         private final double[][][] sCoef;

  1161.         /** G<sub>ms</sub><sup>j</sup> and H<sub>ms</sub><sup>j</sup> polynomials. */
  1162.         private GHmsjPolynomials ghMSJ;

  1164.         /** &Gamma;<sub>ns</sub><sup>m</sup> function. */
  1165.         private GammaMnsFunction gammaMNS;

  1167.         /** R / a up to power degree. */
  1168.         private final double[] roaPow;

  1170.         /** Create a set of C<sub>i</sub><sup>jm</sup> and S<sub>i</sub><sup>jm</sup> coefficients.
  1171.          *  @param jMax absolute limit for j ( -jMax <= j <= jMax )
  1172.          *  @param mMax maximum value for m
  1173.          */
  1174.         FourierCjSjCoefficients(final int jMax, final int mMax) {
  1175.             // initialise fields
  1176.             final int rows    = mMax + 1;
  1177.             final int columns = 2 * jMax + 1;
  1178.             this.jMax         = jMax;
  1179.             this.cCoef        = new double[rows][columns][6];
  1180.             this.sCoef        = new double[rows][columns][6];
  1181.             this.roaPow       = new double[maxDegree + 1];
  1182.             roaPow[0] = 1.;
  1183.         }

  1185.         /**
  1186.          * Generate the coefficients.
  1187.          * @param date the current date
  1188.          * @param context container for attributes
  1189.          * @param hansenObjects initialization of hansen objects
  1190.          */
  1191.         public void generateCoefficients(final AbsoluteDate date, final DSSTTesseralContext context,
  1192.                                          final HansenObjects hansenObjects) {

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

  1196.             // Compute only if there is at least one non-resonant tesseral
  1197.             if (!nonResOrders.isEmpty() || maxDegreeTesseralSP < 0) {
  1198.                 // Gmsj and Hmsj polynomials
  1199.                 ghMSJ = new GHmsjPolynomials(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta(), I);

  1201.                 // GAMMAmns function
  1202.                 gammaMNS = new GammaMnsFunction(maxDegree, auxiliaryElements.getGamma(), I);

  1204.                 final int maxRoaPower = FastMath.max(maxDegreeTesseralSP, maxDegreeMdailyTesseralSP);

  1206.                 // R / a up to power degree
  1207.                 for (int i = 1; i <= maxRoaPower; i++) {
  1208.                     roaPow[i] = context.getRoa() * roaPow[i - 1];
  1209.                 }

  1211.                 //generate the m-daily coefficients
  1212.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1213.                     buildFourierCoefficients(date, m, 0, maxDegreeMdailyTesseralSP, context, hansenObjects);
  1214.                 }

  1216.                 // generate the other coefficients only if required
  1217.                 if (maxDegreeTesseralSP >= 0) {
  1218.                     for (int m: nonResOrders.keySet()) {
  1219.                         final List<Integer> listJ = nonResOrders.get(m);

  1221.                         for (int j: listJ) {
  1222.                             buildFourierCoefficients(date, m, j, maxDegreeTesseralSP, context, hansenObjects);
  1223.                         }
  1224.                     }
  1225.                 }
  1226.             }
  1227.         }

  1229.         /** Build a set of fourier coefficients for a given m and j.
  1230.          *
  1231.          * @param date the date of the coefficients
  1232.          * @param m m index
  1233.          * @param j j index
  1234.          * @param maxN  maximum value for n index
  1235.          * @param context container for attributes
  1236.          * @param hansenObjects initialization of hansen objects
  1237.          */
  1238.         private void buildFourierCoefficients(final AbsoluteDate date,
  1239.                final int m, final int j, final int maxN, final DSSTTesseralContext context,
  1240.                final HansenObjects hansenObjects) {

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

  1244.             // Potential derivatives components for a given non-resonant pair {j,m}
  1245.             double dRdaCos  = 0.;
  1246.             double dRdaSin  = 0.;
  1247.             double dRdhCos  = 0.;
  1248.             double dRdhSin  = 0.;
  1249.             double dRdkCos  = 0.;
  1250.             double dRdkSin  = 0.;
  1251.             double dRdlCos  = 0.;
  1252.             double dRdlSin  = 0.;
  1253.             double dRdAlCos = 0.;
  1254.             double dRdAlSin = 0.;
  1255.             double dRdBeCos = 0.;
  1256.             double dRdBeSin = 0.;
  1257.             double dRdGaCos = 0.;
  1258.             double dRdGaSin = 0.;

  1260.             // s-SUM from -sMin to sMax
  1261.             final int sMin = j == 0 ? maxEccPowMdailyTesseralSP : maxEccPowTesseralSP;
  1262.             final int sMax = j == 0 ? maxEccPowMdailyTesseralSP : maxEccPowTesseralSP;
  1263.             for (int s = 0; s <= sMax; s++) {

  1265.                 // n-SUM for s positive
  1266.                 final double[][] nSumSpos = computeNSum(date, j, m, s, maxN,
  1267.                                                         roaPow, ghMSJ, gammaMNS, context, hansenObjects);
  1268.                 dRdaCos  += nSumSpos[0][0];
  1269.                 dRdaSin  += nSumSpos[0][1];
  1270.                 dRdhCos  += nSumSpos[1][0];
  1271.                 dRdhSin  += nSumSpos[1][1];
  1272.                 dRdkCos  += nSumSpos[2][0];
  1273.                 dRdkSin  += nSumSpos[2][1];
  1274.                 dRdlCos  += nSumSpos[3][0];
  1275.                 dRdlSin  += nSumSpos[3][1];
  1276.                 dRdAlCos += nSumSpos[4][0];
  1277.                 dRdAlSin += nSumSpos[4][1];
  1278.                 dRdBeCos += nSumSpos[5][0];
  1279.                 dRdBeSin += nSumSpos[5][1];
  1280.                 dRdGaCos += nSumSpos[6][0];
  1281.                 dRdGaSin += nSumSpos[6][1];

  1283.                 // n-SUM for s negative
  1284.                 if (s > 0 && s <= sMin) {
  1285.                     final double[][] nSumSneg = computeNSum(date, j, m, -s, maxN,
  1286.                                                             roaPow, ghMSJ, gammaMNS, context, hansenObjects);
  1287.                     dRdaCos  += nSumSneg[0][0];
  1288.                     dRdaSin  += nSumSneg[0][1];
  1289.                     dRdhCos  += nSumSneg[1][0];
  1290.                     dRdhSin  += nSumSneg[1][1];
  1291.                     dRdkCos  += nSumSneg[2][0];
  1292.                     dRdkSin  += nSumSneg[2][1];
  1293.                     dRdlCos  += nSumSneg[3][0];
  1294.                     dRdlSin  += nSumSneg[3][1];
  1295.                     dRdAlCos += nSumSneg[4][0];
  1296.                     dRdAlSin += nSumSneg[4][1];
  1297.                     dRdBeCos += nSumSneg[5][0];
  1298.                     dRdBeSin += nSumSneg[5][1];
  1299.                     dRdGaCos += nSumSneg[6][0];
  1300.                     dRdGaSin += nSumSneg[6][1];
  1301.                 }
  1302.             }
  1303.             dRdaCos  *= -context.getMoa() / auxiliaryElements.getSma();
  1304.             dRdaSin  *= -context.getMoa() / auxiliaryElements.getSma();
  1305.             dRdhCos  *=  context.getMoa();
  1306.             dRdhSin  *=  context.getMoa();
  1307.             dRdkCos  *=  context.getMoa();
  1308.             dRdkSin  *=  context.getMoa();
  1309.             dRdlCos  *=  context.getMoa();
  1310.             dRdlSin  *=  context.getMoa();
  1311.             dRdAlCos *=  context.getMoa();
  1312.             dRdAlSin *=  context.getMoa();
  1313.             dRdBeCos *=  context.getMoa();
  1314.             dRdBeSin *=  context.getMoa();
  1315.             dRdGaCos *=  context.getMoa();
  1316.             dRdGaSin *=  context.getMoa();

  1318.             // Compute the cross derivative operator :
  1319.             final double RAlphaGammaCos   = auxiliaryElements.getAlpha() * dRdGaCos - auxiliaryElements.getGamma() * dRdAlCos;
  1320.             final double RAlphaGammaSin   = auxiliaryElements.getAlpha() * dRdGaSin - auxiliaryElements.getGamma() * dRdAlSin;
  1321.             final double RAlphaBetaCos    = auxiliaryElements.getAlpha() * dRdBeCos - auxiliaryElements.getBeta()  * dRdAlCos;
  1322.             final double RAlphaBetaSin    = auxiliaryElements.getAlpha() * dRdBeSin - auxiliaryElements.getBeta()  * dRdAlSin;
  1323.             final double RBetaGammaCos    =  auxiliaryElements.getBeta() * dRdGaCos - auxiliaryElements.getGamma() * dRdBeCos;
  1324.             final double RBetaGammaSin    =  auxiliaryElements.getBeta() * dRdGaSin - auxiliaryElements.getGamma() * dRdBeSin;
  1325.             final double RhkCos           =     auxiliaryElements.getH() * dRdkCos  -     auxiliaryElements.getK() * dRdhCos;
  1326.             final double RhkSin           =     auxiliaryElements.getH() * dRdkSin  -     auxiliaryElements.getK() * dRdhSin;
  1327.             final double pRagmIqRbgoABCos = (auxiliaryElements.getP() * RAlphaGammaCos - I * auxiliaryElements.getQ() * RBetaGammaCos) * context.getOoAB();
  1328.             final double pRagmIqRbgoABSin = (auxiliaryElements.getP() * RAlphaGammaSin - I * auxiliaryElements.getQ() * RBetaGammaSin) * context.getOoAB();
  1329.             final double RhkmRabmdRdlCos  =  RhkCos - RAlphaBetaCos - dRdlCos;
  1330.             final double RhkmRabmdRdlSin  =  RhkSin - RAlphaBetaSin - dRdlSin;

  1332.             // da/dt
  1333.             cCoef[m][j + jMax][0] = context.getAx2oA() * dRdlCos;
  1334.             sCoef[m][j + jMax][0] = context.getAx2oA() * dRdlSin;

  1336.             // dk/dt
  1337.             cCoef[m][j + jMax][1] = -(context.getBoA() * dRdhCos + auxiliaryElements.getH() * pRagmIqRbgoABCos + auxiliaryElements.getK() * context.getBoABpo() * dRdlCos);
  1338.             sCoef[m][j + jMax][1] = -(context.getBoA() * dRdhSin + auxiliaryElements.getH() * pRagmIqRbgoABSin + auxiliaryElements.getK() * context.getBoABpo() * dRdlSin);

  1340.             // dh/dt
  1341.             cCoef[m][j + jMax][2] = context.getBoA() * dRdkCos + auxiliaryElements.getK() * pRagmIqRbgoABCos - auxiliaryElements.getH() * context.getBoABpo() * dRdlCos;
  1342.             sCoef[m][j + jMax][2] = context.getBoA() * dRdkSin + auxiliaryElements.getK() * pRagmIqRbgoABSin - auxiliaryElements.getH() * context.getBoABpo() * dRdlSin;

  1344.             // dq/dt
  1345.             cCoef[m][j + jMax][3] = context.getCo2AB() * (auxiliaryElements.getQ() * RhkmRabmdRdlCos - I * RAlphaGammaCos);
  1346.             sCoef[m][j + jMax][3] = context.getCo2AB() * (auxiliaryElements.getQ() * RhkmRabmdRdlSin - I * RAlphaGammaSin);

  1348.             // dp/dt
  1349.             cCoef[m][j + jMax][4] = context.getCo2AB() * (auxiliaryElements.getP() * RhkmRabmdRdlCos - RBetaGammaCos);
  1350.             sCoef[m][j + jMax][4] = context.getCo2AB() * (auxiliaryElements.getP() * RhkmRabmdRdlSin - RBetaGammaSin);

  1352.             // dλ/dt
  1353.             cCoef[m][j + jMax][5] = -context.getAx2oA() * dRdaCos + context.getBoABpo() * (auxiliaryElements.getH() * dRdhCos + auxiliaryElements.getK() * dRdkCos) + pRagmIqRbgoABCos;
  1354.             sCoef[m][j + jMax][5] = -context.getAx2oA() * dRdaSin + context.getBoABpo() * (auxiliaryElements.getH() * dRdhSin + auxiliaryElements.getK() * dRdkSin) + pRagmIqRbgoABSin;
  1355.         }

  1357.         /** Get the coefficient C<sub>i</sub><sup>jm</sup>.
  1358.          * @param i i index - corresponds to the required variation
  1359.          * @param j j index
  1360.          * @param m m index
  1361.          * @return the coefficient C<sub>i</sub><sup>jm</sup>
  1362.          */
  1363.         public double getCijm(final int i, final int j, final int m) {
  1364.             return cCoef[m][j + jMax][i];
  1365.         }

  1367.         /** Get the coefficient S<sub>i</sub><sup>jm</sup>.
  1368.          * @param i i index - corresponds to the required variation
  1369.          * @param j j index
  1370.          * @param m m index
  1371.          * @return the coefficient S<sub>i</sub><sup>jm</sup>
  1372.          */
  1373.         public double getSijm(final int i, final int j, final int m) {
  1374.             return sCoef[m][j + jMax][i];
  1375.         }
  1376.     }

  1378.     /** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
  1379.      *  <p>
  1380.      *  Those coefficients are given in Danielson paper by substituting the
  1381.      *  disturbing function (2.7.1-16) with m != 0 into (2.2-10)
  1382.      *  </p>
  1383.      */
  1384.     private class FieldFourierCjSjCoefficients <T extends CalculusFieldElement<T>> {

  1386.         /** Absolute limit for j ( -jMax <= j <= jMax ).  */
  1387.         private final int jMax;

  1389.         /** The C<sub>i</sub><sup>jm</sup> coefficients.
  1390.          * <p>
  1391.          * The index order is [m][j][i] <br/>
  1392.          * The i index corresponds to the C<sub>i</sub><sup>jm</sup> coefficients used to
  1393.          * compute the following: <br/>
  1394.          * - da/dt <br/>
  1395.          * - dk/dt <br/>
  1396.          * - dh/dt / dk <br/>
  1397.          * - dq/dt <br/>
  1398.          * - dp/dt / dα <br/>
  1399.          * - dλ/dt / dβ <br/>
  1400.          * </p>
  1401.          */
  1402.         private final T[][][] cCoef;

  1404.         /** The S<sub>i</sub><sup>jm</sup> coefficients.
  1405.          * <p>
  1406.          * The index order is [m][j][i] <br/>
  1407.          * The i index corresponds to the C<sub>i</sub><sup>jm</sup> coefficients used to
  1408.          * compute the following: <br/>
  1409.          * - da/dt <br/>
  1410.          * - dk/dt <br/>
  1411.          * - dh/dt / dk <br/>
  1412.          * - dq/dt <br/>
  1413.          * - dp/dt / dα <br/>
  1414.          * - dλ/dt / dβ <br/>
  1415.          * </p>
  1416.          */
  1417.         private final T[][][] sCoef;

  1419.         /** G<sub>ms</sub><sup>j</sup> and H<sub>ms</sub><sup>j</sup> polynomials. */
  1420.         private FieldGHmsjPolynomials<T> ghMSJ;

  1422.         /** &Gamma;<sub>ns</sub><sup>m</sup> function. */
  1423.         private FieldGammaMnsFunction<T> gammaMNS;

  1425.         /** R / a up to power degree. */
  1426.         private final T[] roaPow;

  1428.         /** Create a set of C<sub>i</sub><sup>jm</sup> and S<sub>i</sub><sup>jm</sup> coefficients.
  1429.          *  @param jMax absolute limit for j ( -jMax <= j <= jMax )
  1430.          *  @param mMax maximum value for m
  1431.          *  @param field field used by default
  1432.          */
  1433.         FieldFourierCjSjCoefficients(final int jMax, final int mMax, final Field<T> field) {
  1434.             // initialise fields
  1435.             final T zero = field.getZero();
  1436.             final int rows    = mMax + 1;
  1437.             final int columns = 2 * jMax + 1;
  1438.             this.jMax         = jMax;
  1439.             this.cCoef        = MathArrays.buildArray(field, rows, columns, 6);
  1440.             this.sCoef        = MathArrays.buildArray(field, rows, columns, 6);
  1441.             this.roaPow       = MathArrays.buildArray(field, maxDegree + 1);
  1442.             roaPow[0] = zero.add(1.);
  1443.         }

  1445.         /**
  1446.          * Generate the coefficients.
  1447.          * @param date the current date
  1448.          * @param context container for attributes
  1449.          * @param hansenObjects initialization of hansen objects
  1450.          * @param field field used by default
  1451.          */
  1452.         public void generateCoefficients(final FieldAbsoluteDate<T> date,
  1453.                                          final FieldDSSTTesseralContext<T> context,
  1454.                                          final FieldHansenObjects<T> hansenObjects,
  1455.                                          final Field<T> field) {

  1457.             final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
  1458.             // Compute only if there is at least one non-resonant tesseral
  1459.             if (!nonResOrders.isEmpty() || maxDegreeTesseralSP < 0) {
  1460.                 // Gmsj and Hmsj polynomials
  1461.                 ghMSJ = new FieldGHmsjPolynomials<>(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta(), I, field);

  1463.                 // GAMMAmns function
  1464.                 gammaMNS = new FieldGammaMnsFunction<>(maxDegree, auxiliaryElements.getGamma(), I, field);

  1466.                 final int maxRoaPower = FastMath.max(maxDegreeTesseralSP, maxDegreeMdailyTesseralSP);

  1468.                 // R / a up to power degree
  1469.                 for (int i = 1; i <= maxRoaPower; i++) {
  1470.                     roaPow[i] = context.getRoa().multiply(roaPow[i - 1]);
  1471.                 }

  1473.                 //generate the m-daily coefficients
  1474.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1475.                     buildFourierCoefficients(date, m, 0, maxDegreeMdailyTesseralSP, context, hansenObjects, field);
  1476.                 }

  1478.                 // generate the other coefficients only if required
  1479.                 if (maxDegreeTesseralSP >= 0) {
  1480.                     for (int m: nonResOrders.keySet()) {
  1481.                         final List<Integer> listJ = nonResOrders.get(m);

  1483.                         for (int j: listJ) {
  1484.                             buildFourierCoefficients(date, m, j, maxDegreeTesseralSP, context, hansenObjects, field);
  1485.                         }
  1486.                     }
  1487.                 }
  1488.             }
  1489.         }

  1491.         /** Build a set of fourier coefficients for a given m and j.
  1492.          *
  1493.          * @param date the date of the coefficients
  1494.          * @param m m index
  1495.          * @param j j index
  1496.          * @param maxN  maximum value for n index
  1497.          * @param context container for attributes
  1498.          * @param hansenObjects initialization of hansen objects
  1499.          * @param field field used by default
  1500.          */
  1501.         private void buildFourierCoefficients(final FieldAbsoluteDate<T> date,
  1502.                                               final int m, final int j, final int maxN,
  1503.                                               final FieldDSSTTesseralContext<T> context,
  1504.                                               final FieldHansenObjects<T> hansenObjects,
  1505.                                               final Field<T> field) {

  1507.             // Zero
  1508.             final T zero = field.getZero();
  1509.             // Common parameters
  1510.             final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();

  1512.             // Potential derivatives components for a given non-resonant pair {j,m}
  1513.             T dRdaCos  = zero;
  1514.             T dRdaSin  = zero;
  1515.             T dRdhCos  = zero;
  1516.             T dRdhSin  = zero;
  1517.             T dRdkCos  = zero;
  1518.             T dRdkSin  = zero;
  1519.             T dRdlCos  = zero;
  1520.             T dRdlSin  = zero;
  1521.             T dRdAlCos = zero;
  1522.             T dRdAlSin = zero;
  1523.             T dRdBeCos = zero;
  1524.             T dRdBeSin = zero;
  1525.             T dRdGaCos = zero;
  1526.             T dRdGaSin = zero;

  1528.             // s-SUM from -sMin to sMax
  1529.             final int sMin = j == 0 ? maxEccPowMdailyTesseralSP : maxEccPowTesseralSP;
  1530.             final int sMax = j == 0 ? maxEccPowMdailyTesseralSP : maxEccPowTesseralSP;
  1531.             for (int s = 0; s <= sMax; s++) {

  1533.                 // n-SUM for s positive
  1534.                 final T[][] nSumSpos = computeNSum(date, j, m, s, maxN,
  1535.                                                         roaPow, ghMSJ, gammaMNS, context, hansenObjects);
  1536.                 dRdaCos  =  dRdaCos.add(nSumSpos[0][0]);
  1537.                 dRdaSin  =  dRdaSin.add(nSumSpos[0][1]);
  1538.                 dRdhCos  =  dRdhCos.add(nSumSpos[1][0]);
  1539.                 dRdhSin  =  dRdhSin.add(nSumSpos[1][1]);
  1540.                 dRdkCos  =  dRdkCos.add(nSumSpos[2][0]);
  1541.                 dRdkSin  =  dRdkSin.add(nSumSpos[2][1]);
  1542.                 dRdlCos  =  dRdlCos.add(nSumSpos[3][0]);
  1543.                 dRdlSin  =  dRdlSin.add(nSumSpos[3][1]);
  1544.                 dRdAlCos = dRdAlCos.add(nSumSpos[4][0]);
  1545.                 dRdAlSin = dRdAlSin.add(nSumSpos[4][1]);
  1546.                 dRdBeCos = dRdBeCos.add(nSumSpos[5][0]);
  1547.                 dRdBeSin = dRdBeSin.add(nSumSpos[5][1]);
  1548.                 dRdGaCos = dRdGaCos.add(nSumSpos[6][0]);
  1549.                 dRdGaSin = dRdGaSin.add(nSumSpos[6][1]);

  1551.                 // n-SUM for s negative
  1552.                 if (s > 0 && s <= sMin) {
  1553.                     final T[][] nSumSneg = computeNSum(date, j, m, -s, maxN,
  1554.                                                             roaPow, ghMSJ, gammaMNS, context, hansenObjects);
  1555.                     dRdaCos  =  dRdaCos.add(nSumSneg[0][0]);
  1556.                     dRdaSin  =  dRdaSin.add(nSumSneg[0][1]);
  1557.                     dRdhCos  =  dRdhCos.add(nSumSneg[1][0]);
  1558.                     dRdhSin  =  dRdhSin.add(nSumSneg[1][1]);
  1559.                     dRdkCos  =  dRdkCos.add(nSumSneg[2][0]);
  1560.                     dRdkSin  =  dRdkSin.add(nSumSneg[2][1]);
  1561.                     dRdlCos  =  dRdlCos.add(nSumSneg[3][0]);
  1562.                     dRdlSin  =  dRdlSin.add(nSumSneg[3][1]);
  1563.                     dRdAlCos = dRdAlCos.add(nSumSneg[4][0]);
  1564.                     dRdAlSin = dRdAlSin.add(nSumSneg[4][1]);
  1565.                     dRdBeCos = dRdBeCos.add(nSumSneg[5][0]);
  1566.                     dRdBeSin = dRdBeSin.add(nSumSneg[5][1]);
  1567.                     dRdGaCos = dRdGaCos.add(nSumSneg[6][0]);
  1568.                     dRdGaSin = dRdGaSin.add(nSumSneg[6][1]);
  1569.                 }
  1570.             }
  1571.             dRdaCos  =  dRdaCos.multiply(context.getMoa().negate().divide(auxiliaryElements.getSma()));
  1572.             dRdaSin  =  dRdaSin.multiply(context.getMoa().negate().divide(auxiliaryElements.getSma()));
  1573.             dRdhCos  =  dRdhCos.multiply(context.getMoa());
  1574.             dRdhSin  =  dRdhSin.multiply(context.getMoa());
  1575.             dRdkCos  =  dRdkCos.multiply(context.getMoa());
  1576.             dRdkSin  =  dRdkSin.multiply(context.getMoa());
  1577.             dRdlCos  =  dRdlCos.multiply(context.getMoa());
  1578.             dRdlSin  =  dRdlSin.multiply(context.getMoa());
  1579.             dRdAlCos = dRdAlCos.multiply(context.getMoa());
  1580.             dRdAlSin = dRdAlSin.multiply(context.getMoa());
  1581.             dRdBeCos = dRdBeCos.multiply(context.getMoa());
  1582.             dRdBeSin = dRdBeSin.multiply(context.getMoa());
  1583.             dRdGaCos = dRdGaCos.multiply(context.getMoa());
  1584.             dRdGaSin = dRdGaSin.multiply(context.getMoa());

  1586.             // Compute the cross derivative operator :
  1587.             final T RAlphaGammaCos   = auxiliaryElements.getAlpha().multiply(dRdGaCos).subtract(auxiliaryElements.getGamma().multiply(dRdAlCos));
  1588.             final T RAlphaGammaSin   = auxiliaryElements.getAlpha().multiply(dRdGaSin).subtract(auxiliaryElements.getGamma().multiply(dRdAlSin));
  1589.             final T RAlphaBetaCos    = auxiliaryElements.getAlpha().multiply(dRdBeCos).subtract(auxiliaryElements.getBeta().multiply(dRdAlCos));
  1590.             final T RAlphaBetaSin    = auxiliaryElements.getAlpha().multiply(dRdBeSin).subtract(auxiliaryElements.getBeta().multiply(dRdAlSin));
  1591.             final T RBetaGammaCos    =  auxiliaryElements.getBeta().multiply(dRdGaCos).subtract(auxiliaryElements.getGamma().multiply(dRdBeCos));
  1592.             final T RBetaGammaSin    =  auxiliaryElements.getBeta().multiply(dRdGaSin).subtract(auxiliaryElements.getGamma().multiply(dRdBeSin));
  1593.             final T RhkCos           =     auxiliaryElements.getH().multiply(dRdkCos).subtract(auxiliaryElements.getK().multiply(dRdhCos));
  1594.             final T RhkSin           =     auxiliaryElements.getH().multiply(dRdkSin).subtract(auxiliaryElements.getK().multiply(dRdhSin));
  1595.             final T pRagmIqRbgoABCos = (auxiliaryElements.getP().multiply(RAlphaGammaCos).subtract(auxiliaryElements.getQ().multiply(RBetaGammaCos).multiply(I))).multiply(context.getOoAB());
  1596.             final T pRagmIqRbgoABSin = (auxiliaryElements.getP().multiply(RAlphaGammaSin).subtract(auxiliaryElements.getQ().multiply(RBetaGammaSin).multiply(I))).multiply(context.getOoAB());
  1597.             final T RhkmRabmdRdlCos  =  RhkCos.subtract(RAlphaBetaCos).subtract(dRdlCos);
  1598.             final T RhkmRabmdRdlSin  =  RhkSin.subtract(RAlphaBetaSin).subtract(dRdlSin);

  1600.             // da/dt
  1601.             cCoef[m][j + jMax][0] = context.getAx2oA().multiply(dRdlCos);
  1602.             sCoef[m][j + jMax][0] = context.getAx2oA().multiply(dRdlSin);

  1604.             // dk/dt
  1605.             cCoef[m][j + jMax][1] = (context.getBoA().multiply(dRdhCos).add(auxiliaryElements.getH().multiply(pRagmIqRbgoABCos)).add(auxiliaryElements.getK().multiply(context.getBoABpo()).multiply(dRdlCos))).negate();
  1606.             sCoef[m][j + jMax][1] = (context.getBoA().multiply(dRdhSin).add(auxiliaryElements.getH().multiply(pRagmIqRbgoABSin)).add(auxiliaryElements.getK().multiply(context.getBoABpo()).multiply(dRdlSin))).negate();

  1608.             // dh/dt
  1609.             cCoef[m][j + jMax][2] = context.getBoA().multiply(dRdkCos).add(auxiliaryElements.getK().multiply(pRagmIqRbgoABCos)).subtract(auxiliaryElements.getH().multiply(context.getBoABpo()).multiply(dRdlCos));
  1610.             sCoef[m][j + jMax][2] = context.getBoA().multiply(dRdkSin).add(auxiliaryElements.getK().multiply(pRagmIqRbgoABSin)).subtract(auxiliaryElements.getH().multiply(context.getBoABpo()).multiply(dRdlSin));

  1612.             // dq/dt
  1613.             cCoef[m][j + jMax][3] = context.getCo2AB().multiply(auxiliaryElements.getQ().multiply(RhkmRabmdRdlCos).subtract(RAlphaGammaCos.multiply(I)));
  1614.             sCoef[m][j + jMax][3] = context.getCo2AB().multiply(auxiliaryElements.getQ().multiply(RhkmRabmdRdlSin).subtract(RAlphaGammaSin.multiply(I)));

  1616.             // dp/dt
  1617.             cCoef[m][j + jMax][4] = context.getCo2AB().multiply(auxiliaryElements.getP().multiply(RhkmRabmdRdlCos).subtract(RBetaGammaCos));
  1618.             sCoef[m][j + jMax][4] = context.getCo2AB().multiply(auxiliaryElements.getP().multiply(RhkmRabmdRdlSin).subtract(RBetaGammaSin));

  1620.             // dλ/dt
  1621.             cCoef[m][j + jMax][5] = context.getAx2oA().negate().multiply(dRdaCos).add(context.getBoABpo().multiply(auxiliaryElements.getH().multiply(dRdhCos).add(auxiliaryElements.getK().multiply(dRdkCos)))).add(pRagmIqRbgoABCos);
  1622.             sCoef[m][j + jMax][5] = context.getAx2oA().negate().multiply(dRdaSin).add(context.getBoABpo().multiply(auxiliaryElements.getH().multiply(dRdhSin).add(auxiliaryElements.getK().multiply(dRdkSin)))).add(pRagmIqRbgoABSin);
  1623.         }

  1625.         /** Get the coefficient C<sub>i</sub><sup>jm</sup>.
  1626.          * @param i i index - corresponds to the required variation
  1627.          * @param j j index
  1628.          * @param m m index
  1629.          * @return the coefficient C<sub>i</sub><sup>jm</sup>
  1630.          */
  1631.         public T getCijm(final int i, final int j, final int m) {
  1632.             return cCoef[m][j + jMax][i];
  1633.         }

  1635.         /** Get the coefficient S<sub>i</sub><sup>jm</sup>.
  1636.          * @param i i index - corresponds to the required variation
  1637.          * @param j j index
  1638.          * @param m m index
  1639.          * @return the coefficient S<sub>i</sub><sup>jm</sup>
  1640.          */
  1641.         public T getSijm(final int i, final int j, final int m) {
  1642.             return sCoef[m][j + jMax][i];
  1643.         }
  1644.     }

  1646.     /** The C<sup>i</sup><sub>m</sub><sub>t</sub> and S<sup>i</sup><sub>m</sub><sub>t</sub> coefficients used to compute
  1647.      * the short-periodic zonal contribution.
  1648.      *   <p>
  1649.      *  Those coefficients are given by expression 2.5.4-5 from the Danielson paper.
  1650.      *   </p>
  1651.      *
  1652.      * @author Sorin Scortan
  1653.      *
  1654.      */
  1655.     private static class TesseralShortPeriodicCoefficients implements ShortPeriodTerms {

  1657.         /** Retrograde factor I.
  1658.          *  <p>
  1659.          *  DSST model needs equinoctial orbit as internal representation.
  1660.          *  Classical equinoctial elements have discontinuities when inclination
  1661.          *  is close to zero. In this representation, I = +1. <br>
  1662.          *  To avoid this discontinuity, another representation exists and equinoctial
  1663.          *  elements can be expressed in a different way, called "retrograde" orbit.
  1664.          *  This implies I = -1. <br>
  1665.          *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  1666.          *  But for the sake of consistency with the theory, the retrograde factor
  1667.          *  has been kept in the formulas.
  1668.          *  </p>
  1669.          */
  1670.         private static final int I = 1;

  1672.         /** Central body rotating frame. */
  1673.         private final Frame bodyFrame;

  1675.         /** Maximal order to consider for short periodics m-daily tesseral harmonics potential. */
  1676.         private final int maxOrderMdailyTesseralSP;

  1678.         /** Flag to take into account only M-dailies harmonic tesserals for short periodic perturbations.  */
  1679.         private final boolean mDailiesOnly;

  1681.         /** List of non resonant orders with j != 0. */
  1682.         private final SortedMap<Integer, List<Integer> > nonResOrders;

  1684.         /** Maximum value for m index. */
  1685.         private final int mMax;

  1687.         /** Maximum value for j index. */
  1688.         private final int jMax;

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

  1693.         /** All coefficients slots. */
  1694.         private final transient TimeSpanMap<Slot> slots;

  1696.         /** Constructor.
  1697.          * @param bodyFrame central body rotating frame
  1698.          * @param maxOrderMdailyTesseralSP maximal order to consider for short periodics m-daily tesseral harmonics potential
  1699.          * @param mDailiesOnly flag to take into account only M-dailies harmonic tesserals for short periodic perturbations
  1700.          * @param nonResOrders lst of non resonant orders with j != 0
  1701.          * @param mMax maximum value for m index
  1702.          * @param jMax maximum value for j index
  1703.          * @param interpolationPoints number of points used in the interpolation process
  1704.          * @param slots all coefficients slots
  1705.          */
  1706.         TesseralShortPeriodicCoefficients(final Frame bodyFrame, final int maxOrderMdailyTesseralSP,
  1707.                                           final boolean mDailiesOnly, final SortedMap<Integer, List<Integer> > nonResOrders,
  1708.                                           final int mMax, final int jMax, final int interpolationPoints,
  1709.                                           final TimeSpanMap<Slot> slots) {
  1710.             this.bodyFrame                = bodyFrame;
  1711.             this.maxOrderMdailyTesseralSP = maxOrderMdailyTesseralSP;
  1712.             this.mDailiesOnly             = mDailiesOnly;
  1713.             this.nonResOrders             = nonResOrders;
  1714.             this.mMax                     = mMax;
  1715.             this.jMax                     = jMax;
  1716.             this.interpolationPoints      = interpolationPoints;
  1717.             this.slots                    = slots;
  1718.         }

  1720.         /** Get the slot valid for some date.
  1721.          * @param meanStates mean states defining the slot
  1722.          * @return slot valid at the specified date
  1723.          */
  1724.         public Slot createSlot(final SpacecraftState... meanStates) {
  1725.             final Slot         slot  = new Slot(mMax, jMax, interpolationPoints);
  1726.             final AbsoluteDate first = meanStates[0].getDate();
  1727.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  1728.             final int compare = first.compareTo(last);
  1729.             if (compare < 0) {
  1730.                 slots.addValidAfter(slot, first, false);
  1731.             } else if (compare > 0) {
  1732.                 slots.addValidBefore(slot, first, false);
  1733.             } else {
  1734.                 // single date, valid for all time
  1735.                 slots.addValidAfter(slot, AbsoluteDate.PAST_INFINITY, false);
  1736.             }
  1737.             return slot;
  1738.         }

  1740.         /** {@inheritDoc} */
  1741.         @Override
  1742.         public double[] value(final Orbit meanOrbit) {

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

  1747.             // Initialise the short periodic variations
  1748.             final double[] shortPeriodicVariation = new double[6];

  1750.             // Compute only if there is at least one non-resonant tesseral or
  1751.             // only the m-daily tesseral should be taken into account
  1752.             if (!nonResOrders.isEmpty() || mDailiesOnly) {

  1754.                 //Build an auxiliary object
  1755.                 final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(meanOrbit, I);

  1757.                 // Central body rotation angle from equation 2.7.1-(3)(4).
  1758.                 final StaticTransform t = bodyFrame.getStaticTransformTo(
  1759.                         auxiliaryElements.getFrame(),
  1760.                         auxiliaryElements.getDate());
  1761.                 final Vector3D xB = t.transformVector(Vector3D.PLUS_I);
  1762.                 final Vector3D yB = t.transformVector(Vector3D.PLUS_J);
  1763.                 final Vector3D  f = auxiliaryElements.getVectorF();
  1764.                 final Vector3D  g = auxiliaryElements.getVectorG();
  1765.                 final double currentTheta = FastMath.atan2(-f.dotProduct(yB) + I * g.dotProduct(xB),
  1766.                                                             f.dotProduct(xB) + I * g.dotProduct(yB));

  1768.                 //Add the m-daily contribution
  1769.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1770.                     // Phase angle
  1771.                     final double jlMmt  = -m * currentTheta;
  1772.                     final SinCos scPhi  = FastMath.sinCos(jlMmt);
  1773.                     final double sinPhi = scPhi.sin();
  1774.                     final double cosPhi = scPhi.cos();

  1776.                     // compute contribution for each element
  1777.                     final double[] c = slot.getCijm(0, m, meanOrbit.getDate());
  1778.                     final double[] s = slot.getSijm(0, m, meanOrbit.getDate());
  1779.                     for (int i = 0; i < 6; i++) {
  1780.                         shortPeriodicVariation[i] += c[i] * cosPhi + s[i] * sinPhi;
  1781.                     }
  1782.                 }

  1784.                 // loop through all non-resonant (j,m) pairs
  1785.                 for (final Map.Entry<Integer, List<Integer>> entry : nonResOrders.entrySet()) {
  1786.                     final int           m     = entry.getKey();
  1787.                     final List<Integer> listJ = entry.getValue();

  1789.                     for (int j : listJ) {
  1790.                         // Phase angle
  1791.                         final double jlMmt  = j * meanOrbit.getLM() - m * currentTheta;
  1792.                         final SinCos scPhi  = FastMath.sinCos(jlMmt);
  1793.                         final double sinPhi = scPhi.sin();
  1794.                         final double cosPhi = scPhi.cos();

  1796.                         // compute contribution for each element
  1797.                         final double[] c = slot.getCijm(j, m, meanOrbit.getDate());
  1798.                         final double[] s = slot.getSijm(j, m, meanOrbit.getDate());
  1799.                         for (int i = 0; i < 6; i++) {
  1800.                             shortPeriodicVariation[i] += c[i] * cosPhi + s[i] * sinPhi;
  1801.                         }

  1803.                     }
  1804.                 }
  1805.             }

  1807.             return shortPeriodicVariation;

  1809.         }

  1811.         /** {@inheritDoc} */
  1812.         @Override
  1813.         public String getCoefficientsKeyPrefix() {
  1814.             return DSSTTesseral.SHORT_PERIOD_PREFIX;
  1815.         }

  1817.         /** {@inheritDoc}
  1818.          * <p>
  1819.          * For tesseral terms contributions,there are maxOrderMdailyTesseralSP
  1820.          * m-daily cMm coefficients, maxOrderMdailyTesseralSP m-daily sMm
  1821.          * coefficients, nbNonResonant cjm coefficients and nbNonResonant
  1822.          * sjm coefficients, where maxOrderMdailyTesseralSP and nbNonResonant both
  1823.          * depend on the orbit. The j index is the integer multiplier for the true
  1824.          * longitude argument and the m index is the integer multiplier for m-dailies.
  1825.          * </p>
  1826.          */
  1827.         @Override
  1828.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected) {

  1830.             // select the coefficients slot
  1831.             final Slot slot = slots.get(date);

  1833.             if (!nonResOrders.isEmpty() || mDailiesOnly) {
  1834.                 final Map<String, double[]> coefficients =
  1835.                                 new HashMap<String, double[]>(12 * maxOrderMdailyTesseralSP +
  1836.                                                               12 * nonResOrders.size());

  1838.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1839.                     storeIfSelected(coefficients, selected, slot.getCijm(0, m, date), DSSTTesseral.CM_COEFFICIENTS, m);
  1840.                     storeIfSelected(coefficients, selected, slot.getSijm(0, m, date), DSSTTesseral.SM_COEFFICIENTS, m);
  1841.                 }

  1843.                 for (final Map.Entry<Integer, List<Integer>> entry : nonResOrders.entrySet()) {
  1844.                     final int           m     = entry.getKey();
  1845.                     final List<Integer> listJ = entry.getValue();

  1847.                     for (int j : listJ) {
  1848.                         for (int i = 0; i < 6; ++i) {
  1849.                             storeIfSelected(coefficients, selected, slot.getCijm(j, m, date), "c", j, m);
  1850.                             storeIfSelected(coefficients, selected, slot.getSijm(j, m, date), "s", j, m);
  1851.                         }
  1852.                     }
  1853.                 }

  1855.                 return coefficients;

  1857.             } else {
  1858.                 return Collections.emptyMap();
  1859.             }

  1861.         }

  1863.         /** Put a coefficient in a map if selected.
  1864.          * @param map map to populate
  1865.          * @param selected set of coefficients that should be put in the map
  1866.          * (empty set means all coefficients are selected)
  1867.          * @param value coefficient value
  1868.          * @param id coefficient identifier
  1869.          * @param indices list of coefficient indices
  1870.          */
  1871.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
  1872.                                      final double[] value, final String id, final int... indices) {
  1873.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  1874.             keyBuilder.append(id);
  1875.             for (int index : indices) {
  1876.                 keyBuilder.append('[').append(index).append(']');
  1877.             }
  1878.             final String key = keyBuilder.toString();
  1879.             if (selected.isEmpty() || selected.contains(key)) {
  1880.                 map.put(key, value);
  1881.             }
  1882.         }

  1884.     }

  1886.     /** The C<sup>i</sup><sub>m</sub><sub>t</sub> and S<sup>i</sup><sub>m</sub><sub>t</sub> coefficients used to compute
  1887.      * the short-periodic zonal contribution.
  1888.      *   <p>
  1889.      *  Those coefficients are given by expression 2.5.4-5 from the Danielson paper.
  1890.      *   </p>
  1891.      *
  1892.      * @author Sorin Scortan
  1893.      *
  1894.      */
  1895.     private static class FieldTesseralShortPeriodicCoefficients <T extends CalculusFieldElement<T>> implements FieldShortPeriodTerms<T> {

  1897.         /** Retrograde factor I.
  1898.          *  <p>
  1899.          *  DSST model needs equinoctial orbit as internal representation.
  1900.          *  Classical equinoctial elements have discontinuities when inclination
  1901.          *  is close to zero. In this representation, I = +1. <br>
  1902.          *  To avoid this discontinuity, another representation exists and equinoctial
  1903.          *  elements can be expressed in a different way, called "retrograde" orbit.
  1904.          *  This implies I = -1. <br>
  1905.          *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  1906.          *  But for the sake of consistency with the theory, the retrograde factor
  1907.          *  has been kept in the formulas.
  1908.          *  </p>
  1909.          */
  1910.         private static final int I = 1;

  1912.         /** Central body rotating frame. */
  1913.         private final Frame bodyFrame;

  1915.         /** Maximal order to consider for short periodics m-daily tesseral harmonics potential. */
  1916.         private final int maxOrderMdailyTesseralSP;

  1918.         /** Flag to take into account only M-dailies harmonic tesserals for short periodic perturbations.  */
  1919.         private final boolean mDailiesOnly;

  1921.         /** List of non resonant orders with j != 0. */
  1922.         private final SortedMap<Integer, List<Integer> > nonResOrders;

  1924.         /** Maximum value for m index. */
  1925.         private final int mMax;

  1927.         /** Maximum value for j index. */
  1928.         private final int jMax;

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

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

  1936.         /** Constructor.
  1937.          * @param bodyFrame central body rotating frame
  1938.          * @param maxOrderMdailyTesseralSP maximal order to consider for short periodics m-daily tesseral harmonics potential
  1939.          * @param mDailiesOnly flag to take into account only M-dailies harmonic tesserals for short periodic perturbations
  1940.          * @param nonResOrders lst of non resonant orders with j != 0
  1941.          * @param mMax maximum value for m index
  1942.          * @param jMax maximum value for j index
  1943.          * @param interpolationPoints number of points used in the interpolation process
  1944.          * @param slots all coefficients slots
  1945.          */
  1946.         FieldTesseralShortPeriodicCoefficients(final Frame bodyFrame, final int maxOrderMdailyTesseralSP,
  1947.                                                final boolean mDailiesOnly, final SortedMap<Integer, List<Integer> > nonResOrders,
  1948.                                                final int mMax, final int jMax, final int interpolationPoints,
  1949.                                                final FieldTimeSpanMap<FieldSlot<T>, T> slots) {
  1950.             this.bodyFrame                = bodyFrame;
  1951.             this.maxOrderMdailyTesseralSP = maxOrderMdailyTesseralSP;
  1952.             this.mDailiesOnly             = mDailiesOnly;
  1953.             this.nonResOrders             = nonResOrders;
  1954.             this.mMax                     = mMax;
  1955.             this.jMax                     = jMax;
  1956.             this.interpolationPoints      = interpolationPoints;
  1957.             this.slots                    = slots;
  1958.         }

  1960.         /** Get the slot valid for some date.
  1961.          * @param meanStates mean states defining the slot
  1962.          * @return slot valid at the specified date
  1963.          */
  1964.         @SuppressWarnings("unchecked")
  1965.         public FieldSlot<T> createSlot(final FieldSpacecraftState<T>... meanStates) {
  1966.             final FieldSlot<T>         slot  = new FieldSlot<>(mMax, jMax, interpolationPoints);
  1967.             final FieldAbsoluteDate<T> first = meanStates[0].getDate();
  1968.             final FieldAbsoluteDate<T> last  = meanStates[meanStates.length - 1].getDate();
  1969.             if (first.compareTo(last) <= 0) {
  1970.                 slots.addValidAfter(slot, first);
  1971.             } else {
  1972.                 slots.addValidBefore(slot, first);
  1973.             }
  1974.             return slot;
  1975.         }

  1977.         /** {@inheritDoc} */
  1978.         @Override
  1979.         public T[] value(final FieldOrbit<T> meanOrbit) {

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

  1984.             // Initialise the short periodic variations
  1985.             final T[] shortPeriodicVariation = MathArrays.buildArray(meanOrbit.getDate().getField(), 6);

  1987.             // Compute only if there is at least one non-resonant tesseral or
  1988.             // only the m-daily tesseral should be taken into account
  1989.             if (!nonResOrders.isEmpty() || mDailiesOnly) {

  1991.                 //Build an auxiliary object
  1992.                 final FieldAuxiliaryElements<T> auxiliaryElements = new FieldAuxiliaryElements<>(meanOrbit, I);

  1994.                 // Central body rotation angle from equation 2.7.1-(3)(4).
  1995.                 final FieldStaticTransform<T> t = bodyFrame.getStaticTransformTo(auxiliaryElements.getFrame(), auxiliaryElements.getDate());
  1996.                 final FieldVector3D<T> xB = t.transformVector(Vector3D.PLUS_I);
  1997.                 final FieldVector3D<T> yB = t.transformVector(Vector3D.PLUS_J);
  1998.                 final FieldVector3D<T>  f = auxiliaryElements.getVectorF();
  1999.                 final FieldVector3D<T>  g = auxiliaryElements.getVectorG();
  2000.                 final T currentTheta = FastMath.atan2(f.dotProduct(yB).negate().add(g.dotProduct(xB).multiply(I)),
  2001.                                                       f.dotProduct(xB).add(g.dotProduct(yB).multiply(I)));

  2003.                 //Add the m-daily contribution
  2004.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  2005.                     // Phase angle
  2006.                     final T jlMmt              = currentTheta.multiply(-m);
  2007.                     final FieldSinCos<T> scPhi = FastMath.sinCos(jlMmt);
  2008.                     final T sinPhi             = scPhi.sin();
  2009.                     final T cosPhi             = scPhi.cos();

  2011.                     // compute contribution for each element
  2012.                     final T[] c = slot.getCijm(0, m, meanOrbit.getDate());
  2013.                     final T[] s = slot.getSijm(0, m, meanOrbit.getDate());
  2014.                     for (int i = 0; i < 6; i++) {
  2015.                         shortPeriodicVariation[i] = shortPeriodicVariation[i].add(c[i].multiply(cosPhi).add(s[i].multiply(sinPhi)));
  2016.                     }
  2017.                 }

  2019.                 // loop through all non-resonant (j,m) pairs
  2020.                 for (final Map.Entry<Integer, List<Integer>> entry : nonResOrders.entrySet()) {
  2021.                     final int           m     = entry.getKey();
  2022.                     final List<Integer> listJ = entry.getValue();

  2024.                     for (int j : listJ) {
  2025.                         // Phase angle
  2026.                         final T jlMmt              = meanOrbit.getLM().multiply(j).subtract(currentTheta.multiply(m));
  2027.                         final FieldSinCos<T> scPhi = FastMath.sinCos(jlMmt);
  2028.                         final T sinPhi             = scPhi.sin();
  2029.                         final T cosPhi             = scPhi.cos();

  2031.                         // compute contribution for each element
  2032.                         final T[] c = slot.getCijm(j, m, meanOrbit.getDate());
  2033.                         final T[] s = slot.getSijm(j, m, meanOrbit.getDate());
  2034.                         for (int i = 0; i < 6; i++) {
  2035.                             shortPeriodicVariation[i] = shortPeriodicVariation[i].add(c[i].multiply(cosPhi).add(s[i].multiply(sinPhi)));
  2036.                         }

  2038.                     }
  2039.                 }
  2040.             }

  2042.             return shortPeriodicVariation;

  2044.         }

  2046.         /** {@inheritDoc} */
  2047.         @Override
  2048.         public String getCoefficientsKeyPrefix() {
  2049.             return DSSTTesseral.SHORT_PERIOD_PREFIX;
  2050.         }

  2052.         /** {@inheritDoc}
  2053.          * <p>
  2054.          * For tesseral terms contributions,there are maxOrderMdailyTesseralSP
  2055.          * m-daily cMm coefficients, maxOrderMdailyTesseralSP m-daily sMm
  2056.          * coefficients, nbNonResonant cjm coefficients and nbNonResonant
  2057.          * sjm coefficients, where maxOrderMdailyTesseralSP and nbNonResonant both
  2058.          * depend on the orbit. The j index is the integer multiplier for the true
  2059.          * longitude argument and the m index is the integer multiplier for m-dailies.
  2060.          * </p>
  2061.          */
  2062.         @Override
  2063.         public Map<String, T[]> getCoefficients(final FieldAbsoluteDate<T> date, final Set<String> selected) {

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

  2068.             if (!nonResOrders.isEmpty() || mDailiesOnly) {
  2069.                 final Map<String, T[]> coefficients =
  2070.                                 new HashMap<String, T[]>(12 * maxOrderMdailyTesseralSP +
  2071.                                                          12 * nonResOrders.size());

  2073.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  2074.                     storeIfSelected(coefficients, selected, slot.getCijm(0, m, date), DSSTTesseral.CM_COEFFICIENTS, m);
  2075.                     storeIfSelected(coefficients, selected, slot.getSijm(0, m, date), DSSTTesseral.SM_COEFFICIENTS, m);
  2076.                 }

  2078.                 for (final Map.Entry<Integer, List<Integer>> entry : nonResOrders.entrySet()) {
  2079.                     final int           m     = entry.getKey();
  2080.                     final List<Integer> listJ = entry.getValue();

  2082.                     for (int j : listJ) {
  2083.                         for (int i = 0; i < 6; ++i) {
  2084.                             storeIfSelected(coefficients, selected, slot.getCijm(j, m, date), "c", j, m);
  2085.                             storeIfSelected(coefficients, selected, slot.getSijm(j, m, date), "s", j, m);
  2086.                         }
  2087.                     }
  2088.                 }

  2090.                 return coefficients;

  2092.             } else {
  2093.                 return Collections.emptyMap();
  2094.             }

  2096.         }

  2098.         /** Put a coefficient in a map if selected.
  2099.          * @param map map to populate
  2100.          * @param selected set of coefficients that should be put in the map
  2101.          * (empty set means all coefficients are selected)
  2102.          * @param value coefficient value
  2103.          * @param id coefficient identifier
  2104.          * @param indices list of coefficient indices
  2105.          */
  2106.         private void storeIfSelected(final Map<String, T[]> map, final Set<String> selected,
  2107.                                      final T[] value, final String id, final int... indices) {
  2108.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  2109.             keyBuilder.append(id);
  2110.             for (int index : indices) {
  2111.                 keyBuilder.append('[').append(index).append(']');
  2112.             }
  2113.             final String key = keyBuilder.toString();
  2114.             if (selected.isEmpty() || selected.contains(key)) {
  2115.                 map.put(key, value);
  2116.             }
  2117.         }
  2118.     }

  2120.     /** Coefficients valid for one time slot. */
  2121.     private static class Slot {

  2123.         /** The coefficients C<sub>i</sub><sup>j</sup><sup>m</sup>.
  2124.          * <p>
  2125.          * The index order is cijm[m][j][i] <br/>
  2126.          * i corresponds to the equinoctial element, as follows: <br/>
  2127.          * - i=0 for a <br/>
  2128.          * - i=1 for k <br/>
  2129.          * - i=2 for h <br/>
  2130.          * - i=3 for q <br/>
  2131.          * - i=4 for p <br/>
  2132.          * - i=5 for λ <br/>
  2133.          * </p>
  2134.          */
  2135.         private final ShortPeriodicsInterpolatedCoefficient[][] cijm;

  2137.         /** The coefficients S<sub>i</sub><sup>j</sup><sup>m</sup>.
  2138.          * <p>
  2139.          * The index order is sijm[m][j][i] <br/>
  2140.          * i corresponds to the equinoctial element, as follows: <br/>
  2141.          * - i=0 for a <br/>
  2142.          * - i=1 for k <br/>
  2143.          * - i=2 for h <br/>
  2144.          * - i=3 for q <br/>
  2145.          * - i=4 for p <br/>
  2146.          * - i=5 for λ <br/>
  2147.          * </p>
  2148.          */
  2149.         private final ShortPeriodicsInterpolatedCoefficient[][] sijm;

  2151.         /** Simple constructor.
  2152.          *  @param mMax maximum value for m index
  2153.          *  @param jMax maximum value for j index
  2154.          *  @param interpolationPoints number of points used in the interpolation process
  2155.          */
  2156.         Slot(final int mMax, final int jMax, final int interpolationPoints) {

  2158.             final int rows    = mMax + 1;
  2159.             final int columns = 2 * jMax + 1;
  2160.             cijm = new ShortPeriodicsInterpolatedCoefficient[rows][columns];
  2161.             sijm = new ShortPeriodicsInterpolatedCoefficient[rows][columns];
  2162.             for (int m = 1; m <= mMax; m++) {
  2163.                 for (int j = -jMax; j <= jMax; j++) {
  2164.                     cijm[m][j + jMax] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  2165.                     sijm[m][j + jMax] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  2166.                 }
  2167.             }

  2169.         }

  2171.         /** Get C<sub>i</sub><sup>j</sup><sup>m</sup>.
  2172.          *
  2173.          * @param j j index
  2174.          * @param m m index
  2175.          * @param date the date
  2176.          * @return C<sub>i</sub><sup>j</sup><sup>m</sup>
  2177.          */
  2178.         double[] getCijm(final int j, final int m, final AbsoluteDate date) {
  2179.             final int jMax = (cijm[m].length - 1) / 2;
  2180.             return cijm[m][j + jMax].value(date);
  2181.         }

  2183.         /** Get S<sub>i</sub><sup>j</sup><sup>m</sup>.
  2184.          *
  2185.          * @param j j index
  2186.          * @param m m index
  2187.          * @param date the date
  2188.          * @return S<sub>i</sub><sup>j</sup><sup>m</sup>
  2189.          */
  2190.         double[] getSijm(final int j, final int m, final AbsoluteDate date) {
  2191.             final int jMax = (cijm[m].length - 1) / 2;
  2192.             return sijm[m][j + jMax].value(date);
  2193.         }

  2195.     }

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

  2200.         /** The coefficients C<sub>i</sub><sup>j</sup><sup>m</sup>.
  2201.          * <p>
  2202.          * The index order is cijm[m][j][i] <br/>
  2203.          * i corresponds to the equinoctial element, as follows: <br/>
  2204.          * - i=0 for a <br/>
  2205.          * - i=1 for k <br/>
  2206.          * - i=2 for h <br/>
  2207.          * - i=3 for q <br/>
  2208.          * - i=4 for p <br/>
  2209.          * - i=5 for λ <br/>
  2210.          * </p>
  2211.          */
  2212.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[][] cijm;

  2214.         /** The coefficients S<sub>i</sub><sup>j</sup><sup>m</sup>.
  2215.          * <p>
  2216.          * The index order is sijm[m][j][i] <br/>
  2217.          * i corresponds to the equinoctial element, as follows: <br/>
  2218.          * - i=0 for a <br/>
  2219.          * - i=1 for k <br/>
  2220.          * - i=2 for h <br/>
  2221.          * - i=3 for q <br/>
  2222.          * - i=4 for p <br/>
  2223.          * - i=5 for λ <br/>
  2224.          * </p>
  2225.          */
  2226.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[][] sijm;

  2228.         /** Simple constructor.
  2229.          *  @param mMax maximum value for m index
  2230.          *  @param jMax maximum value for j index
  2231.          *  @param interpolationPoints number of points used in the interpolation process
  2232.          */
  2233.         @SuppressWarnings("unchecked")
  2234.         FieldSlot(final int mMax, final int jMax, final int interpolationPoints) {

  2236.             final int rows    = mMax + 1;
  2237.             final int columns = 2 * jMax + 1;
  2238.             cijm = (FieldShortPeriodicsInterpolatedCoefficient<T>[][]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows, columns);
  2239.             sijm = (FieldShortPeriodicsInterpolatedCoefficient<T>[][]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows, columns);
  2240.             for (int m = 1; m <= mMax; m++) {
  2241.                 for (int j = -jMax; j <= jMax; j++) {
  2242.                     cijm[m][j + jMax] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  2243.                     sijm[m][j + jMax] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  2244.                 }
  2245.             }

  2247.         }

  2249.         /** Get C<sub>i</sub><sup>j</sup><sup>m</sup>.
  2250.          *
  2251.          * @param j j index
  2252.          * @param m m index
  2253.          * @param date the date
  2254.          * @return C<sub>i</sub><sup>j</sup><sup>m</sup>
  2255.          */
  2256.         T[] getCijm(final int j, final int m, final FieldAbsoluteDate<T> date) {
  2257.             final int jMax = (cijm[m].length - 1) / 2;
  2258.             return cijm[m][j + jMax].value(date);
  2259.         }

  2261.         /** Get S<sub>i</sub><sup>j</sup><sup>m</sup>.
  2262.          *
  2263.          * @param j j index
  2264.          * @param m m index
  2265.          * @param date the date
  2266.          * @return S<sub>i</sub><sup>j</sup><sup>m</sup>
  2267.          */
  2268.         T[] getSijm(final int j, final int m, final FieldAbsoluteDate<T> date) {
  2269.             final int jMax = (cijm[m].length - 1) / 2;
  2270.             return sijm[m][j + jMax].value(date);
  2271.         }

  2273.     }

  2275.     /** Compute potential and potential derivatives with respect to orbital parameters.
  2276.      *  <p>The following elements are computed from expression 3.3 - (4).
  2277.      *  <pre>
  2278.      *  dU / da
  2279.      *  dU / dh
  2280.      *  dU / dk
  2281.      *  dU / dλ
  2282.      *  dU / dα
  2283.      *  dU / dβ
  2284.      *  dU / dγ
  2285.      *  </pre>
  2286.      *  </p>
  2287.      */
  2288.     private class UAnddU {

  2290.         /** dU / da. */
  2291.         private  double dUda;

  2293.         /** dU / dk. */
  2294.         private double dUdk;

  2296.         /** dU / dh. */
  2297.         private double dUdh;

  2299.         /** dU / dl. */
  2300.         private double dUdl;

  2302.         /** dU / dAlpha. */
  2303.         private double dUdAl;

  2305.         /** dU / dBeta. */
  2306.         private double dUdBe;

  2308.         /** dU / dGamma. */
  2309.         private double dUdGa;

  2311.         /** Simple constuctor.
  2312.          * @param date current date
  2313.          * @param context container for attributes
  2314.          * @param hansen hansen objects
  2315.          */
  2316.         UAnddU(final AbsoluteDate date, final DSSTTesseralContext context, final HansenObjects hansen) {

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

  2321.             // Potential derivatives
  2322.             dUda  = 0.;
  2323.             dUdh  = 0.;
  2324.             dUdk  = 0.;
  2325.             dUdl  = 0.;
  2326.             dUdAl = 0.;
  2327.             dUdBe = 0.;
  2328.             dUdGa = 0.;

  2330.             // Compute only if there is at least one resonant tesseral
  2331.             if (!resOrders.isEmpty()) {
  2332.                 // Gmsj and Hmsj polynomials
  2333.                 final GHmsjPolynomials ghMSJ = new GHmsjPolynomials(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta(), I);

  2335.                 // GAMMAmns function
  2336.                 final GammaMnsFunction gammaMNS = new GammaMnsFunction(maxDegree, auxiliaryElements.getGamma(), I);

  2338.                 // R / a up to power degree
  2339.                 final double[] roaPow = new double[maxDegree + 1];
  2340.                 roaPow[0] = 1.;
  2341.                 for (int i = 1; i <= maxDegree; i++) {
  2342.                     roaPow[i] = context.getRoa() * roaPow[i - 1];
  2343.                 }

  2345.                 // SUM over resonant terms {j,m}
  2346.                 for (int m : resOrders) {

  2348.                     // Resonant index for the current resonant order
  2349.                     final int j = FastMath.max(1, (int) FastMath.round(context.getRatio() * m));

  2351.                     // Phase angle
  2352.                     final double jlMmt  = j * auxiliaryElements.getLM() - m * context.getTheta();
  2353.                     final SinCos scPhi  = FastMath.sinCos(jlMmt);
  2354.                     final double sinPhi = scPhi.sin();
  2355.                     final double cosPhi = scPhi.cos();

  2357.                     // Potential derivatives components for a given resonant pair {j,m}
  2358.                     double dUdaCos  = 0.;
  2359.                     double dUdaSin  = 0.;
  2360.                     double dUdhCos  = 0.;
  2361.                     double dUdhSin  = 0.;
  2362.                     double dUdkCos  = 0.;
  2363.                     double dUdkSin  = 0.;
  2364.                     double dUdlCos  = 0.;
  2365.                     double dUdlSin  = 0.;
  2366.                     double dUdAlCos = 0.;
  2367.                     double dUdAlSin = 0.;
  2368.                     double dUdBeCos = 0.;
  2369.                     double dUdBeSin = 0.;
  2370.                     double dUdGaCos = 0.;
  2371.                     double dUdGaSin = 0.;

  2373.                     // s-SUM from -sMin to sMax
  2374.                     final int sMin = FastMath.min(maxEccPow - j, maxDegree);
  2375.                     final int sMax = FastMath.min(maxEccPow + j, maxDegree);
  2376.                     for (int s = 0; s <= sMax; s++) {

  2378.                         //Compute the initial values for Hansen coefficients using newComb operators
  2379.                         hansen.computeHansenObjectsInitValues(context, s + maxDegree, j);

  2381.                         // n-SUM for s positive
  2382.                         final double[][] nSumSpos = computeNSum(date, j, m, s, maxDegree,
  2383.                                                                 roaPow, ghMSJ, gammaMNS, context, hansen);
  2384.                         dUdaCos  += nSumSpos[0][0];
  2385.                         dUdaSin  += nSumSpos[0][1];
  2386.                         dUdhCos  += nSumSpos[1][0];
  2387.                         dUdhSin  += nSumSpos[1][1];
  2388.                         dUdkCos  += nSumSpos[2][0];
  2389.                         dUdkSin  += nSumSpos[2][1];
  2390.                         dUdlCos  += nSumSpos[3][0];
  2391.                         dUdlSin  += nSumSpos[3][1];
  2392.                         dUdAlCos += nSumSpos[4][0];
  2393.                         dUdAlSin += nSumSpos[4][1];
  2394.                         dUdBeCos += nSumSpos[5][0];
  2395.                         dUdBeSin += nSumSpos[5][1];
  2396.                         dUdGaCos += nSumSpos[6][0];
  2397.                         dUdGaSin += nSumSpos[6][1];

  2399.                         // n-SUM for s negative
  2400.                         if (s > 0 && s <= sMin) {
  2401.                             //Compute the initial values for Hansen coefficients using newComb operators
  2402.                             hansen.computeHansenObjectsInitValues(context, maxDegree - s, j);

  2404.                             final double[][] nSumSneg = computeNSum(date, j, m, -s, maxDegree,
  2405.                                                                     roaPow, ghMSJ, gammaMNS, context, hansen);
  2406.                             dUdaCos  += nSumSneg[0][0];
  2407.                             dUdaSin  += nSumSneg[0][1];
  2408.                             dUdhCos  += nSumSneg[1][0];
  2409.                             dUdhSin  += nSumSneg[1][1];
  2410.                             dUdkCos  += nSumSneg[2][0];
  2411.                             dUdkSin  += nSumSneg[2][1];
  2412.                             dUdlCos  += nSumSneg[3][0];
  2413.                             dUdlSin  += nSumSneg[3][1];
  2414.                             dUdAlCos += nSumSneg[4][0];
  2415.                             dUdAlSin += nSumSneg[4][1];
  2416.                             dUdBeCos += nSumSneg[5][0];
  2417.                             dUdBeSin += nSumSneg[5][1];
  2418.                             dUdGaCos += nSumSneg[6][0];
  2419.                             dUdGaSin += nSumSneg[6][1];
  2420.                         }
  2421.                     }

  2423.                     // Assembly of potential derivatives componants
  2424.                     dUda  += cosPhi * dUdaCos  + sinPhi * dUdaSin;
  2425.                     dUdh  += cosPhi * dUdhCos  + sinPhi * dUdhSin;
  2426.                     dUdk  += cosPhi * dUdkCos  + sinPhi * dUdkSin;
  2427.                     dUdl  += cosPhi * dUdlCos  + sinPhi * dUdlSin;
  2428.                     dUdAl += cosPhi * dUdAlCos + sinPhi * dUdAlSin;
  2429.                     dUdBe += cosPhi * dUdBeCos + sinPhi * dUdBeSin;
  2430.                     dUdGa += cosPhi * dUdGaCos + sinPhi * dUdGaSin;
  2431.                 }

  2433.                 this.dUda  = dUda * (-context.getMoa() / auxiliaryElements.getSma());
  2434.                 this.dUdh  = dUdh * context.getMoa();
  2435.                 this.dUdk  = dUdk * context.getMoa();
  2436.                 this.dUdl  = dUdl * context.getMoa();
  2437.                 this.dUdAl = dUdAl * context.getMoa();
  2438.                 this.dUdBe = dUdBe * context.getMoa();
  2439.                 this.dUdGa = dUdGa * context.getMoa();
  2440.             }

  2442.         }

  2444.         /** Return value of dU / da.
  2445.          * @return dUda
  2446.          */
  2447.         public double getdUda() {
  2448.             return dUda;
  2449.         }

  2451.         /** Return value of dU / dk.
  2452.          * @return dUdk
  2453.          */
  2454.         public double getdUdk() {
  2455.             return dUdk;
  2456.         }

  2458.         /** Return value of dU / dh.
  2459.          * @return dUdh
  2460.          */
  2461.         public double getdUdh() {
  2462.             return dUdh;
  2463.         }

  2465.         /** Return value of dU / dl.
  2466.          * @return dUdl
  2467.          */
  2468.         public double getdUdl() {
  2469.             return dUdl;
  2470.         }

  2472.         /** Return value of dU / dAlpha.
  2473.          * @return dUdAl
  2474.          */
  2475.         public double getdUdAl() {
  2476.             return dUdAl;
  2477.         }

  2479.         /** Return value of dU / dBeta.
  2480.          * @return dUdBe
  2481.          */
  2482.         public double getdUdBe() {
  2483.             return dUdBe;
  2484.         }

  2486.         /** Return value of dU / dGamma.
  2487.          * @return dUdGa
  2488.          */
  2489.         public double getdUdGa() {
  2490.             return dUdGa;
  2491.         }

  2493.     }

  2495.     /**  Computes the potential U derivatives.
  2496.      *  <p>The following elements are computed from expression 3.3 - (4).
  2497.      *  <pre>
  2498.      *  dU / da
  2499.      *  dU / dh
  2500.      *  dU / dk
  2501.      *  dU / dλ
  2502.      *  dU / dα
  2503.      *  dU / dβ
  2504.      *  dU / dγ
  2505.      *  </pre>
  2506.      *  </p>
  2507.      */
  2508.     private class FieldUAnddU <T extends CalculusFieldElement<T>> {

  2510.         /** dU / da. */
  2511.         private T dUda;

  2513.         /** dU / dk. */
  2514.         private T dUdk;

  2516.         /** dU / dh. */
  2517.         private T dUdh;

  2519.         /** dU / dl. */
  2520.         private T dUdl;

  2522.         /** dU / dAlpha. */
  2523.         private T dUdAl;

  2525.         /** dU / dBeta. */
  2526.         private T dUdBe;

  2528.         /** dU / dGamma. */
  2529.         private T dUdGa;

  2531.         /** Simple constuctor.
  2532.          * @param date current date
  2533.          * @param context container for attributes
  2534.          * @param hansen hansen objects
  2535.          */
  2536.         FieldUAnddU(final FieldAbsoluteDate<T> date, final FieldDSSTTesseralContext<T> context,
  2537.                     final FieldHansenObjects<T> hansen) {

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

  2542.             // Zero for initialization
  2543.             final Field<T> field = date.getField();
  2544.             final T zero = field.getZero();

  2546.             // Potential derivatives
  2547.             dUda  = zero;
  2548.             dUdh  = zero;
  2549.             dUdk  = zero;
  2550.             dUdl  = zero;
  2551.             dUdAl = zero;
  2552.             dUdBe = zero;
  2553.             dUdGa = zero;

  2555.             // Compute only if there is at least one resonant tesseral
  2556.             if (!resOrders.isEmpty()) {
  2557.                 // Gmsj and Hmsj polynomials
  2558.                 final FieldGHmsjPolynomials<T> ghMSJ = new FieldGHmsjPolynomials<>(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta(), I, field);

  2560.                 // GAMMAmns function
  2561.                 final FieldGammaMnsFunction<T> gammaMNS = new FieldGammaMnsFunction<>(maxDegree, auxiliaryElements.getGamma(), I, field);

  2563.                 // R / a up to power degree
  2564.                 final T[] roaPow = MathArrays.buildArray(field, maxDegree + 1);
  2565.                 roaPow[0] = zero.add(1.);
  2566.                 for (int i = 1; i <= maxDegree; i++) {
  2567.                     roaPow[i] = roaPow[i - 1].multiply(context.getRoa());
  2568.                 }

  2570.                 // SUM over resonant terms {j,m}
  2571.                 for (int m : resOrders) {

  2573.                     // Resonant index for the current resonant order
  2574.                     final int j = FastMath.max(1, (int) FastMath.round(context.getRatio().multiply(m)));

  2576.                     // Phase angle
  2577.                     final T jlMmt              = auxiliaryElements.getLM().multiply(j).subtract(context.getTheta().multiply(m));
  2578.                     final FieldSinCos<T> scPhi = FastMath.sinCos(jlMmt);
  2579.                     final T sinPhi             = scPhi.sin();
  2580.                     final T cosPhi             = scPhi.cos();

  2582.                     // Potential derivatives components for a given resonant pair {j,m}
  2583.                     T dUdaCos  = zero;
  2584.                     T dUdaSin  = zero;
  2585.                     T dUdhCos  = zero;
  2586.                     T dUdhSin  = zero;
  2587.                     T dUdkCos  = zero;
  2588.                     T dUdkSin  = zero;
  2589.                     T dUdlCos  = zero;
  2590.                     T dUdlSin  = zero;
  2591.                     T dUdAlCos = zero;
  2592.                     T dUdAlSin = zero;
  2593.                     T dUdBeCos = zero;
  2594.                     T dUdBeSin = zero;
  2595.                     T dUdGaCos = zero;
  2596.                     T dUdGaSin = zero;

  2598.                     // s-SUM from -sMin to sMax
  2599.                     final int sMin = FastMath.min(maxEccPow - j, maxDegree);
  2600.                     final int sMax = FastMath.min(maxEccPow + j, maxDegree);
  2601.                     for (int s = 0; s <= sMax; s++) {

  2603.                         //Compute the initial values for Hansen coefficients using newComb operators
  2604.                         hansen.computeHansenObjectsInitValues(context, s + maxDegree, j);

  2606.                         // n-SUM for s positive
  2607.                         final T[][] nSumSpos = computeNSum(date, j, m, s, maxDegree,
  2608.                                                                 roaPow, ghMSJ, gammaMNS, context, hansen);
  2609.                         dUdaCos  = dUdaCos.add(nSumSpos[0][0]);
  2610.                         dUdaSin  = dUdaSin.add(nSumSpos[0][1]);
  2611.                         dUdhCos  = dUdhCos.add(nSumSpos[1][0]);
  2612.                         dUdhSin  = dUdhSin.add(nSumSpos[1][1]);
  2613.                         dUdkCos  = dUdkCos.add(nSumSpos[2][0]);
  2614.                         dUdkSin  = dUdkSin.add(nSumSpos[2][1]);
  2615.                         dUdlCos  = dUdlCos.add(nSumSpos[3][0]);
  2616.                         dUdlSin  = dUdlSin.add(nSumSpos[3][1]);
  2617.                         dUdAlCos = dUdAlCos.add(nSumSpos[4][0]);
  2618.                         dUdAlSin = dUdAlSin.add(nSumSpos[4][1]);
  2619.                         dUdBeCos = dUdBeCos.add(nSumSpos[5][0]);
  2620.                         dUdBeSin = dUdBeSin.add(nSumSpos[5][1]);
  2621.                         dUdGaCos = dUdGaCos.add(nSumSpos[6][0]);
  2622.                         dUdGaSin = dUdGaSin.add(nSumSpos[6][1]);

  2624.                         // n-SUM for s negative
  2625.                         if (s > 0 && s <= sMin) {
  2626.                             //Compute the initial values for Hansen coefficients using newComb operators
  2627.                             hansen.computeHansenObjectsInitValues(context, maxDegree - s, j);

  2629.                             final T[][] nSumSneg = computeNSum(date, j, m, -s, maxDegree,
  2630.                                                                     roaPow, ghMSJ, gammaMNS, context, hansen);
  2631.                             dUdaCos  = dUdaCos.add(nSumSneg[0][0]);
  2632.                             dUdaSin  = dUdaSin.add(nSumSneg[0][1]);
  2633.                             dUdhCos  = dUdhCos.add(nSumSneg[1][0]);
  2634.                             dUdhSin  = dUdhSin.add(nSumSneg[1][1]);
  2635.                             dUdkCos  = dUdkCos.add(nSumSneg[2][0]);
  2636.                             dUdkSin  = dUdkSin.add(nSumSneg[2][1]);
  2637.                             dUdlCos  = dUdlCos.add(nSumSneg[3][0]);
  2638.                             dUdlSin  = dUdlSin.add(nSumSneg[3][1]);
  2639.                             dUdAlCos = dUdAlCos.add(nSumSneg[4][0]);
  2640.                             dUdAlSin = dUdAlSin.add(nSumSneg[4][1]);
  2641.                             dUdBeCos = dUdBeCos.add(nSumSneg[5][0]);
  2642.                             dUdBeSin = dUdBeSin.add(nSumSneg[5][1]);
  2643.                             dUdGaCos = dUdGaCos.add(nSumSneg[6][0]);
  2644.                             dUdGaSin = dUdGaSin.add(nSumSneg[6][1]);
  2645.                         }
  2646.                     }

  2648.                     // Assembly of potential derivatives componants
  2649.                     dUda  = dUda.add(dUdaCos.multiply(cosPhi).add(dUdaSin.multiply(sinPhi)));
  2650.                     dUdh  = dUdh.add(dUdhCos.multiply(cosPhi).add(dUdhSin.multiply(sinPhi)));
  2651.                     dUdk  = dUdk.add(dUdkCos.multiply(cosPhi).add(dUdkSin.multiply(sinPhi)));
  2652.                     dUdl  = dUdl.add(dUdlCos.multiply(cosPhi).add(dUdlSin.multiply(sinPhi)));
  2653.                     dUdAl = dUdAl.add(dUdAlCos.multiply(cosPhi).add(dUdAlSin.multiply(sinPhi)));
  2654.                     dUdBe = dUdBe.add(dUdBeCos.multiply(cosPhi).add(dUdBeSin.multiply(sinPhi)));
  2655.                     dUdGa = dUdGa.add(dUdGaCos.multiply(cosPhi).add(dUdGaSin.multiply(sinPhi)));
  2656.                 }

  2658.                 dUda  =  dUda.multiply(context.getMoa().divide(auxiliaryElements.getSma())).negate();
  2659.                 dUdh  =  dUdh.multiply(context.getMoa());
  2660.                 dUdk  =  dUdk.multiply(context.getMoa());
  2661.                 dUdl  =  dUdl.multiply(context.getMoa());
  2662.                 dUdAl =  dUdAl.multiply(context.getMoa());
  2663.                 dUdBe =  dUdBe.multiply(context.getMoa());
  2664.                 dUdGa =  dUdGa.multiply(context.getMoa());

  2666.             }

  2668.         }

  2670.         /** Return value of dU / da.
  2671.          * @return dUda
  2672.          */
  2673.         public T getdUda() {
  2674.             return dUda;
  2675.         }

  2677.         /** Return value of dU / dk.
  2678.          * @return dUdk
  2679.          */
  2680.         public T getdUdk() {
  2681.             return dUdk;
  2682.         }

  2684.         /** Return value of dU / dh.
  2685.          * @return dUdh
  2686.          */
  2687.         public T getdUdh() {
  2688.             return dUdh;
  2689.         }

  2691.         /** Return value of dU / dl.
  2692.          * @return dUdl
  2693.          */
  2694.         public T getdUdl() {
  2695.             return dUdl;
  2696.         }

  2698.         /** Return value of dU / dAlpha.
  2699.          * @return dUdAl
  2700.          */
  2701.         public T getdUdAl() {
  2702.             return dUdAl;
  2703.         }

  2705.         /** Return value of dU / dBeta.
  2706.          * @return dUdBe
  2707.          */
  2708.         public T getdUdBe() {
  2709.             return dUdBe;
  2710.         }

  2712.         /** Return value of dU / dGamma.
  2713.          * @return dUdGa
  2714.          */
  2715.         public T getdUdGa() {
  2716.             return dUdGa;
  2717.         }

  2719.     }

  2721.     /** Computes init values of the Hansen Objects. */
  2722.     private class HansenObjects {

  2724.         /** A two dimensional array that contains the objects needed to build the Hansen coefficients. <br/>
  2725.          * The indexes are s + maxDegree and j */
  2726.         private HansenTesseralLinear[][] hansenObjects;

  2728.         /** Simple constructor.
  2729.          * @param ratio Ratio of satellite period to central body rotation period
  2730.          * @param type type of the elements used during the propagation
  2731.          */
  2732.         HansenObjects(final double ratio,
  2733.                       final PropagationType type) {

  2735.             //Allocate the two dimensional array
  2736.             final int rows     = 2 * maxDegree + 1;
  2737.             final int columns  = maxFrequencyShortPeriodics + 1;
  2738.             this.hansenObjects = new HansenTesseralLinear[rows][columns];

  2740.             switch (type) {
  2741.                 case MEAN:
  2742.                     // loop through the resonant orders
  2743.                     for (int m : resOrders) {
  2744.                         //Compute the corresponding j term
  2745.                         final int j = FastMath.max(1, (int) FastMath.round(ratio * m));

  2747.                         //Compute the sMin and sMax values
  2748.                         final int sMin = FastMath.min(maxEccPow - j, maxDegree);
  2749.                         final int sMax = FastMath.min(maxEccPow + j, maxDegree);

  2751.                         //loop through the s values
  2752.                         for (int s = 0; s <= sMax; s++) {
  2753.                             //Compute the n0 value
  2754.                             final int n0 = FastMath.max(FastMath.max(2, m), s);

  2756.                             //Create the object for the pair j, s
  2757.                             this.hansenObjects[s + maxDegree][j] = new HansenTesseralLinear(maxDegree, s, j, n0, maxHansen);

  2759.                             if (s > 0 && s <= sMin) {
  2760.                                 //Also create the object for the pair j, -s
  2761.                                 this.hansenObjects[maxDegree - s][j] =  new HansenTesseralLinear(maxDegree, -s, j, n0, maxHansen);
  2762.                             }
  2763.                         }
  2764.                     }
  2765.                     break;

  2767.                 case OSCULATING:
  2768.                     // create all objects
  2769.                     for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
  2770.                         for (int s = -maxDegree; s <= maxDegree; s++) {
  2771.                             //Compute the n0 value
  2772.                             final int n0 = FastMath.max(2, FastMath.abs(s));
  2773.                             this.hansenObjects[s + maxDegree][j] = new HansenTesseralLinear(maxDegree, s, j, n0, maxHansen);
  2774.                         }
  2775.                     }
  2776.                     break;

  2778.                 default:
  2779.                     throw new OrekitInternalError(null);
  2780.             }

  2782.         }

  2784.         /** Compute init values for hansen objects.
  2785.          * @param context container for attributes
  2786.          * @param rows number of rows of the hansen matrix
  2787.          * @param columns columns number of columns of the hansen matrix
  2788.          */
  2789.         public void computeHansenObjectsInitValues(final DSSTTesseralContext context, final int rows, final int columns) {
  2790.             hansenObjects[rows][columns].computeInitValues(context.getE2(), context.getChi(), context.getChi2());
  2791.         }

  2793.         /** Get hansen object.
  2794.          * @return hansenObjects
  2795.          */
  2796.         public HansenTesseralLinear[][] getHansenObjects() {
  2797.             return hansenObjects;
  2798.         }

  2800.     }

  2802.     /** Computes init values of the Hansen Objects. */
  2803.     private class FieldHansenObjects<T extends CalculusFieldElement<T>> {

  2805.         /** A two dimensional array that contains the objects needed to build the Hansen coefficients. <br/>
  2806.          * The indexes are s + maxDegree and j */
  2807.         private FieldHansenTesseralLinear<T>[][] hansenObjects;

  2809.         /** Simple constructor.
  2810.          * @param ratio Ratio of satellite period to central body rotation period
  2811.          * @param type type of the elements used during the propagation
  2812.          */
  2813.         @SuppressWarnings("unchecked")
  2814.         FieldHansenObjects(final T ratio,
  2815.                            final PropagationType type) {

  2817.             // Set the maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
  2818.             maxHansen = maxEccPow / 2;

  2820.             //Allocate the two dimensional array
  2821.             final int rows     = 2 * maxDegree + 1;
  2822.             final int columns  = maxFrequencyShortPeriodics + 1;
  2823.             this.hansenObjects = (FieldHansenTesseralLinear<T>[][]) Array.newInstance(FieldHansenTesseralLinear.class, rows, columns);

  2825.             switch (type) {
  2826.                 case MEAN:
  2827.                  // loop through the resonant orders
  2828.                     for (int m : resOrders) {
  2829.                         //Compute the corresponding j term
  2830.                         final int j = FastMath.max(1, (int) FastMath.round(ratio.multiply(m)));

  2832.                         //Compute the sMin and sMax values
  2833.                         final int sMin = FastMath.min(maxEccPow - j, maxDegree);
  2834.                         final int sMax = FastMath.min(maxEccPow + j, maxDegree);

  2836.                         //loop through the s values
  2837.                         for (int s = 0; s <= sMax; s++) {
  2838.                             //Compute the n0 value
  2839.                             final int n0 = FastMath.max(FastMath.max(2, m), s);

  2841.                             //Create the object for the pair j, s
  2842.                             this.hansenObjects[s + maxDegree][j] = new FieldHansenTesseralLinear<>(maxDegree, s, j, n0, maxHansen, ratio.getField());

  2844.                             if (s > 0 && s <= sMin) {
  2845.                                 //Also create the object for the pair j, -s
  2846.                                 this.hansenObjects[maxDegree - s][j] =  new FieldHansenTesseralLinear<>(maxDegree, -s, j, n0, maxHansen, ratio.getField());
  2847.                             }
  2848.                         }
  2849.                     }
  2850.                     break;

  2852.                 case OSCULATING:
  2853.                     // create all objects
  2854.                     for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
  2855.                         for (int s = -maxDegree; s <= maxDegree; s++) {
  2856.                             //Compute the n0 value
  2857.                             final int n0 = FastMath.max(2, FastMath.abs(s));
  2858.                             this.hansenObjects[s + maxDegree][j] = new FieldHansenTesseralLinear<>(maxDegree, s, j, n0, maxHansen, ratio.getField());
  2859.                         }
  2860.                     }
  2861.                     break;

  2863.                 default:
  2864.                     throw new OrekitInternalError(null);
  2865.             }

  2867.         }

  2869.         /** Compute init values for hansen objects.
  2870.          * @param context container for attributes
  2871.          * @param rows number of rows of the hansen matrix
  2872.          * @param columns columns number of columns of the hansen matrix
  2873.          */
  2874.         public void computeHansenObjectsInitValues(final FieldDSSTTesseralContext<T> context,
  2875.                                                    final int rows, final int columns) {
  2876.             hansenObjects[rows][columns].computeInitValues(context.getE2(), context.getChi(), context.getChi2());
  2877.         }

  2879.         /** Get hansen object.
  2880.          * @return hansenObjects
  2881.          */
  2882.         public FieldHansenTesseralLinear<T>[][] getHansenObjects() {
  2883.             return hansenObjects;
  2884.         }

  2886.     }

  2888. }
Created with JaCoCo 0.8.11.202310140853