DSSTTesseral.java

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

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

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

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

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

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

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

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

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

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

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

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

  131.     /** Provider for spherical harmonics. */
  132.     private final UnnormalizedSphericalHarmonicsProvider provider;

  134.     /** Central body rotating frame. */
  135.     private final Frame bodyFrame;

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

  140.     /** Central body rotation period (seconds). */
  141.     private final double bodyPeriod;

  143.     /** Maximal degree to consider for harmonics potential. */
  144.     private final int maxDegree;

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

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

  152.     /** Maximal order to consider for harmonics potential. */
  153.     private final int maxOrder;

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

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

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

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

  169.     /** Maximum value for j. */
  170.     private final int maxFrequencyShortPeriodics;

  172.     /** Maximum value between maxOrderMdailyTesseralSP and maxOrderTesseralSP. */
  173.     private int mMax;

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

  178.     /** Short period terms. */
  179.     private TesseralShortPeriodicCoefficients shortPeriodTerms;

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

  184.     /** Driver for gravitational parameter. */
  185.     private final ParameterDriver gmParameterDriver;

  187.     /** Hansen objects. */
  188.     private HansenObjects hansen;

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

  193.     /** Flag for force model initialization with field elements. */
  194.     private boolean pendingInitialization;

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

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

  255.         gmParameterDriver = new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
  256.                                                 provider.getMu(), MU_SCALE,
  257.                                                 0.0, Double.POSITIVE_INFINITY);

  259.         // Central body rotating frame
  260.         this.bodyFrame = centralBodyFrame;

  262.         //Save the rotation rate
  263.         this.centralBodyRotationRate = centralBodyRotationRate;

  265.         // Central body rotation period in seconds
  266.         this.bodyPeriod = MathUtils.TWO_PI / centralBodyRotationRate;

  268.         // Provider for spherical harmonics
  269.         this.provider      = provider;
  270.         this.maxDegree     = provider.getMaxDegree();
  271.         this.maxOrder      = provider.getMaxOrder();

  273.         //set the maximum degree order for short periodics
  274.         checkIndexRange(maxDegreeTesseralSP, 2, maxDegree);
  275.         this.maxDegreeTesseralSP       = maxDegreeTesseralSP;

  277.         checkIndexRange(maxDegreeMdailyTesseralSP, 2, maxDegree);
  278.         this.maxDegreeMdailyTesseralSP = maxDegreeMdailyTesseralSP;

  280.         checkIndexRange(maxOrderTesseralSP, 0, maxOrder);
  281.         this.maxOrderTesseralSP        = maxOrderTesseralSP;

  283.         checkIndexRange(maxOrderMdailyTesseralSP, 0, maxOrder);
  284.         this.maxOrderMdailyTesseralSP  = maxOrderMdailyTesseralSP;

  286.         // set the maximum value for eccentricity power
  287.         this.maxEccPowTesseralSP       = maxEccPowTesseralSP;

  289.         checkIndexRange(maxEccPowMdailyTesseralSP, 0, maxDegreeMdailyTesseralSP - 2);
  290.         this.maxEccPowMdailyTesseralSP = maxEccPowMdailyTesseralSP;

  292.         // set the maximum value for frequency
  293.         this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;

  295.         // Initialize default values
  296.         this.nonResOrders = new TreeMap<Integer, List <Integer> >();

  298.         pendingInitialization = true;

  300.         fieldShortPeriodTerms = new HashMap<>();
  301.         fieldHansen           = new HashMap<>();

  303.     }

  305.     /** Check an index range.
  306.      * @param index index value
  307.      * @param min minimum value for index
  308.      * @param max maximum value for index
  309.      */
  310.     private void checkIndexRange(final int index, final int min, final int max) {
  311.         if (index < min || index > max) {
  312.             throw new OrekitException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, index, min, max);
  313.         }
  314.     }

  316.     /** {@inheritDoc} */
  317.     @Override
  318.     public List<ShortPeriodTerms> initialize(final AuxiliaryElements auxiliaryElements,
  319.                                              final PropagationType type,
  320.                                              final double[] parameters) {

  322.         // Initializes specific parameters.
  323.         final DSSTTesseralContext context = initializeStep(auxiliaryElements, parameters);

  325.         // The following terms are only used for hansen objects initialization
  326.         final double ratio            = context.getRatio();
  327.         final int maxEccPow           = context.getMaxEccPow();
  328.         final List<Integer> resOrders = context.getResOrders();

  330.         // Compute the non resonant tesseral harmonic terms if not set by the user
  331.         getNonResonantTerms(type, context);

  333.         hansen = new HansenObjects(ratio, maxEccPow, resOrders, type);

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

  337.         shortPeriodTerms = new TesseralShortPeriodicCoefficients(bodyFrame, maxOrderMdailyTesseralSP,
  338.                                                                  maxDegreeTesseralSP < 0, nonResOrders,
  339.                                                                  mMax, maxFrequencyShortPeriodics, INTERPOLATION_POINTS,
  340.                                                                  new TimeSpanMap<Slot>(new Slot(mMax, maxFrequencyShortPeriodics,
  341.                                                                                                 INTERPOLATION_POINTS)));

  343.         final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
  344.         list.add(shortPeriodTerms);
  345.         return list;

  347.     }

  349.     /** {@inheritDoc} */
  350.     @Override
  351.     public <T extends RealFieldElement<T>> List<FieldShortPeriodTerms<T>> initialize(final FieldAuxiliaryElements<T> auxiliaryElements,
  352.                                                                                      final PropagationType type,
  353.                                                                                      final T[] parameters) {

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

  358.         if (pendingInitialization == true) {

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

  363.             // Compute the non resonant tesseral harmonic terms if not set by the user
  364.             getNonResonantTerms(type, context, field);

  366.             // The following terms are only used for hansen objects initialization
  367.             final T      ratio            = context.getRatio();
  368.             final int maxEccPow           = context.getMaxEccPow();
  369.             final List<Integer> resOrders = context.getResOrders();


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

  374.             fieldHansen.put(field, new FieldHansenObjects<>(ratio, maxEccPow, resOrders, type, field));

  376.             pendingInitialization = false;
  377.         }

  379.         final FieldTesseralShortPeriodicCoefficients<T> ftspc =
  380.                         new FieldTesseralShortPeriodicCoefficients<>(bodyFrame, maxOrderMdailyTesseralSP,
  381.                                                                      maxDegreeTesseralSP < 0, nonResOrders,
  382.                                                                      mMax, maxFrequencyShortPeriodics, INTERPOLATION_POINTS,
  383.                                                                      new FieldTimeSpanMap<>(new FieldSlot<>(mMax,
  384.                                                                                                             maxFrequencyShortPeriodics,
  385.                                                                                                             INTERPOLATION_POINTS),
  386.                                                                                             field));

  388.         fieldShortPeriodTerms.put(field, ftspc);
  389.         return Collections.singletonList(ftspc);

  391.     }

  393.     /** Performs initialization at each integration step for the current force model.
  394.      *  <p>
  395.      *  This method aims at being called before mean elements rates computation.
  396.      *  </p>
  397.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  398.      *  @param parameters values of the force model parameters
  399.      *  @return new force model context
  400.      */
  401.     private DSSTTesseralContext initializeStep(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
  402.         return new DSSTTesseralContext(auxiliaryElements, bodyFrame, provider, maxFrequencyShortPeriodics, bodyPeriod, parameters);
  403.     }

  405.     /** Performs initialization at each integration step for the current force model.
  406.      *  <p>
  407.      *  This method aims at being called before mean elements rates computation.
  408.      *  </p>
  409.      *  @param <T> type of the elements
  410.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  411.      *  @param parameters values of the force model parameters
  412.      *  @return new force model context
  413.      */
  414.     private <T extends RealFieldElement<T>> FieldDSSTTesseralContext<T> initializeStep(final FieldAuxiliaryElements<T> auxiliaryElements,
  415.                                                                                        final T[] parameters) {
  416.         return new FieldDSSTTesseralContext<>(auxiliaryElements, bodyFrame, provider, maxFrequencyShortPeriodics, bodyPeriod, parameters);
  417.     }

  419.     /** {@inheritDoc} */
  420.     @Override
  421.     public double[] getMeanElementRate(final SpacecraftState spacecraftState,
  422.                                        final AuxiliaryElements auxiliaryElements, final double[] parameters) {

  424.         // Container for attributes
  425.         final DSSTTesseralContext context = initializeStep(auxiliaryElements, parameters);

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

  430.         // Compute the cross derivative operator :
  431.         final double UAlphaGamma   = auxiliaryElements.getAlpha() * udu.getdUdGa() - auxiliaryElements.getGamma() * udu.getdUdAl();
  432.         final double UAlphaBeta    = auxiliaryElements.getAlpha() * udu.getdUdBe() - auxiliaryElements.getBeta()  * udu.getdUdAl();
  433.         final double UBetaGamma    = auxiliaryElements.getBeta() * udu.getdUdGa() - auxiliaryElements.getGamma() * udu.getdUdBe();
  434.         final double Uhk           = auxiliaryElements.getH() * udu.getdUdk()  - auxiliaryElements.getK() * udu.getdUdh();
  435.         final double pUagmIqUbgoAB = (auxiliaryElements.getP() * UAlphaGamma - I * auxiliaryElements.getQ() * UBetaGamma) * context.getOoAB();
  436.         final double UhkmUabmdUdl  = Uhk - UAlphaBeta - udu.getdUdl();

  438.         final double da =  context.getAx2oA() * udu.getdUdl();
  439.         final double dh =  context.getBoA() * udu.getdUdk() + auxiliaryElements.getK() * pUagmIqUbgoAB - auxiliaryElements.getH() * context.getBoABpo() * udu.getdUdl();
  440.         final double dk =  -(context.getBoA() * udu.getdUdh() + auxiliaryElements.getH() * pUagmIqUbgoAB + auxiliaryElements.getK() * context.getBoABpo() * udu.getdUdl());
  441.         final double dp =  context.getCo2AB() * (auxiliaryElements.getP() * UhkmUabmdUdl - UBetaGamma);
  442.         final double dq =  context.getCo2AB() * (auxiliaryElements.getQ() * UhkmUabmdUdl - I * UAlphaGamma);
  443.         final double dM = -context.getAx2oA() * udu.getdUda() + context.getBoABpo() * (auxiliaryElements.getH() * udu.getdUdh() + auxiliaryElements.getK() * udu.getdUdk()) + pUagmIqUbgoAB;

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

  448.     /** {@inheritDoc} */
  449.     @Override
  450.     public <T extends RealFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> spacecraftState,
  451.                                                                   final FieldAuxiliaryElements<T> auxiliaryElements,
  452.                                                                   final T[] parameters) {

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

  457.         // Container for attributes
  458.         final FieldDSSTTesseralContext<T> context = initializeStep(auxiliaryElements, parameters);

  460.         @SuppressWarnings("unchecked")
  461.         final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);
  462.         // Access to potential U derivatives
  463.         final FieldUAnddU<T> udu = new FieldUAnddU<>(spacecraftState.getDate(), context, fho);

  465.         // Compute the cross derivative operator :
  466.         final T UAlphaGamma   = udu.getdUdGa().multiply(auxiliaryElements.getAlpha()).subtract(udu.getdUdAl().multiply(auxiliaryElements.getGamma()));
  467.         final T UAlphaBeta    = udu.getdUdBe().multiply(auxiliaryElements.getAlpha()).subtract(udu.getdUdAl().multiply(auxiliaryElements.getBeta()));
  468.         final T UBetaGamma    = udu.getdUdGa().multiply(auxiliaryElements.getBeta()).subtract(udu.getdUdBe().multiply(auxiliaryElements.getGamma()));
  469.         final T Uhk           = udu.getdUdk().multiply(auxiliaryElements.getH()).subtract(udu.getdUdh().multiply(auxiliaryElements.getK()));
  470.         final T pUagmIqUbgoAB = (UAlphaGamma.multiply(auxiliaryElements.getP()).subtract(UBetaGamma.multiply(auxiliaryElements.getQ()).multiply(I))).multiply(context.getOoAB());
  471.         final T UhkmUabmdUdl  = Uhk.subtract(UAlphaBeta).subtract(udu.getdUdl());

  473.         final T da = udu.getdUdl().multiply(context.getAx2oA());
  474.         final T dh = udu.getdUdk().multiply(context.getBoA()).add(pUagmIqUbgoAB.multiply(auxiliaryElements.getK())).subtract(udu.getdUdl().multiply(auxiliaryElements.getH()).multiply(context.getBoABpo()));
  475.         final T dk = (udu.getdUdh().multiply(context.getBoA()).add(pUagmIqUbgoAB.multiply(auxiliaryElements.getH())).add(udu.getdUdl().multiply(context.getBoABpo()).multiply(auxiliaryElements.getK()))).negate();
  476.         final T dp = context.getCo2AB().multiply(auxiliaryElements.getP().multiply(UhkmUabmdUdl).subtract(UBetaGamma));
  477.         final T dq = context.getCo2AB().multiply(auxiliaryElements.getQ().multiply(UhkmUabmdUdl).subtract(UAlphaGamma.multiply(I)));
  478.         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()));

  480.         final T[] elements = MathArrays.buildArray(field, 6);
  481.         elements[0] = da;
  482.         elements[1] = dk;
  483.         elements[2] = dh;
  484.         elements[3] = dq;
  485.         elements[4] = dp;
  486.         elements[5] = dM;

  488.         return elements;

  490.     }

  492.     /** {@inheritDoc} */
  493.     @Override
  494.     public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {

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

  498.         for (final SpacecraftState meanState : meanStates) {

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

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

  504.             // Initialise the Hansen coefficients
  505.             for (int s = -maxDegree; s <= maxDegree; s++) {
  506.                 // coefficients with j == 0 are always needed
  507.                 hansen.computeHansenObjectsInitValues(context, s + maxDegree, 0);
  508.                 if (maxDegreeTesseralSP >= 0) {
  509.                     // initialize other objects only if required
  510.                     for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  511.                         hansen.computeHansenObjectsInitValues(context, s + maxDegree, j);
  512.                     }
  513.                 }
  514.             }

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

  518.             // Compute coefficients
  519.             // Compute only if there is at least one non-resonant tesseral
  520.             if (!nonResOrders.isEmpty() || maxDegreeTesseralSP < 0) {
  521.                 // Generate the fourrier coefficients
  522.                 cjsjFourier.generateCoefficients(meanState.getDate(), context, hansen);

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

  527.                 // build the mDaily coefficients
  528.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  529.                     // build the coefficients
  530.                     buildCoefficients(cjsjFourier, meanState.getDate(), slot, m, 0, tnota, context);
  531.                 }

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

  537.                         for (int j : entry.getValue()) {
  538.                             // build the coefficients
  539.                             buildCoefficients(cjsjFourier, meanState.getDate(), slot, entry.getKey(), j, tnota, context);
  540.                         }
  541.                     }
  542.                 }
  543.             }

  545.         }

  547.     }

  549.     /** {@inheritDoc} */
  550.     @Override
  551.     @SuppressWarnings("unchecked")
  552.     public <T extends RealFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
  553.                                                                        final FieldSpacecraftState<T>... meanStates) {

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

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

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

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

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

  567.             final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);
  568.             // Initialise the Hansen coefficients
  569.             for (int s = -maxDegree; s <= maxDegree; s++) {
  570.                 // coefficients with j == 0 are always needed
  571.                 fho.computeHansenObjectsInitValues(context, s + maxDegree, 0);
  572.                 if (maxDegreeTesseralSP >= 0) {
  573.                     // initialize other objects only if required
  574.                     for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  575.                         fho.computeHansenObjectsInitValues(context, s + maxDegree, j);
  576.                     }
  577.                 }
  578.             }

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

  582.             // Compute coefficients
  583.             // Compute only if there is at least one non-resonant tesseral
  584.             if (!nonResOrders.isEmpty() || maxDegreeTesseralSP < 0) {
  585.                 // Generate the fourrier coefficients
  586.                 cjsjFourier.generateCoefficients(meanState.getDate(), context, fho, field);

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

  591.                 // build the mDaily coefficients
  592.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  593.                     // build the coefficients
  594.                     buildCoefficients(cjsjFourier, meanState.getDate(), slot, m, 0, tnota, context, field);
  595.                 }

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

  601.                         for (int j : entry.getValue()) {
  602.                             // build the coefficients
  603.                             buildCoefficients(cjsjFourier, meanState.getDate(), slot, entry.getKey(), j, tnota, context, field);
  604.                         }
  605.                     }
  606.                 }
  607.             }

  609.         }

  611.     }

  613.     /** {@inheritDoc} */
  614.     public ParameterDriver[] getParametersDrivers() {
  615.         return new ParameterDriver[] {
  616.             gmParameterDriver
  617.         };
  618.     }

  620.     /** Build a set of coefficients.
  621.      * @param cjsjFourier the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  622.      * @param date the current date
  623.      * @param slot slot to which the coefficients belong
  624.      * @param m m index
  625.      * @param j j index
  626.      * @param tnota 3n/2a
  627.      * @param context container for attributes
  628.      */
  629.     private void buildCoefficients(final FourierCjSjCoefficients cjsjFourier,
  630.                                    final AbsoluteDate date, final Slot slot,
  631.                                    final int m, final int j, final double tnota, final DSSTTesseralContext context) {

  633.         // Create local arrays
  634.         final double[] currentCijm = new double[] {0., 0., 0., 0., 0., 0.};
  635.         final double[] currentSijm = new double[] {0., 0., 0., 0., 0., 0.};

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

  640.         // initialise the coeficients
  641.         for (int i = 0; i < 6; i++) {
  642.             currentCijm[i] = -cjsjFourier.getSijm(i, j, m);
  643.             currentSijm[i] = cjsjFourier.getCijm(i, j, m);
  644.         }
  645.         // Add the separate part for δ<sub>6i</sub>
  646.         currentCijm[5] += tnota * oojnmt * cjsjFourier.getCijm(0, j, m);
  647.         currentSijm[5] += tnota * oojnmt * cjsjFourier.getSijm(0, j, m);

  649.         //Multiply by 1 / (jn - mθ<sup>.</sup>)
  650.         for (int i = 0; i < 6; i++) {
  651.             currentCijm[i] *= oojnmt;
  652.             currentSijm[i] *= oojnmt;
  653.         }

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

  659.     }

  661.      /** Build a set of coefficients.
  662.      * @param <T> the type of the field elements
  663.      * @param cjsjFourier the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  664.      * @param date the current date
  665.      * @param slot slot to which the coefficients belong
  666.      * @param m m index
  667.      * @param j j index
  668.      * @param tnota 3n/2a
  669.      * @param context container for attributes
  670.      * @param field field used by default
  671.      */
  672.     private <T extends RealFieldElement<T>> void buildCoefficients(final FieldFourierCjSjCoefficients<T> cjsjFourier,
  673.                                                                    final FieldAbsoluteDate<T> date,
  674.                                                                    final FieldSlot<T> slot,
  675.                                                                    final int m, final int j, final T tnota,
  676.                                                                    final FieldDSSTTesseralContext<T> context,
  677.                                                                    final Field<T> field) {

  679.         // Zero
  680.         final T zero = field.getZero();

  682.         // Create local arrays
  683.         final T[] currentCijm = MathArrays.buildArray(field, 6);
  684.         final T[] currentSijm = MathArrays.buildArray(field, 6);

  686.         Arrays.fill(currentCijm, zero);
  687.         Arrays.fill(currentSijm, zero);

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

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

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

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

  711.     }

  713.     /** {@inheritDoc} */
  714.     @Override
  715.     public EventDetector[] getEventsDetectors() {
  716.         return null;
  717.     }

  719.     /** {@inheritDoc} */
  720.     @Override
  721.     public <T extends RealFieldElement<T>> FieldEventDetector<T>[] getFieldEventsDetectors(final Field<T> field) {
  722.         return null;
  723.     }

  725.      /**
  726.       * Get the non-resonant tesseral terms in the central body spherical harmonic field.
  727.       *
  728.       * @param type type of the elements used during the propagation
  729.       * @param context container for attributes
  730.       */
  731.     private void getNonResonantTerms(final PropagationType type, final DSSTTesseralContext context) {

  733.         // Compute natural resonant terms
  734.         final double tolerance = 1. / FastMath.max(MIN_PERIOD_IN_SAT_REV,
  735.                                                    MIN_PERIOD_IN_SECONDS / context.getOrbitPeriod());

  737.         nonResOrders.clear();
  738.         for (int m = 1; m <= maxOrder; m++) {
  739.             final double resonance = context.getRatio() * m;
  740.             int jRes = 0;
  741.             final int jComputedRes = (int) FastMath.round(resonance);
  742.             if (jComputedRes > 0 && jComputedRes <= maxFrequencyShortPeriodics && FastMath.abs(resonance - jComputedRes) <= tolerance) {
  743.                 // Store each resonant index and order
  744.                 jRes = jComputedRes;
  745.             }

  747.             if (type == PropagationType.OSCULATING && maxDegreeTesseralSP >= 0 && m <= maxOrderTesseralSP) {
  748.                 //compute non resonant orders in the tesseral harmonic field
  749.                 final List<Integer> listJofM = new ArrayList<Integer>();
  750.                 //for the moment we take only the pairs (j,m) with |j| <= maxDegree + maxEccPow (from |s-j| <= maxEccPow and |s| <= maxDegree)
  751.                 for (int j = -maxFrequencyShortPeriodics; j <= maxFrequencyShortPeriodics; j++) {
  752.                     if (j != 0 && j != jRes) {
  753.                         listJofM.add(j);
  754.                     }
  755.                 }

  757.                 nonResOrders.put(m, listJofM);
  758.             }
  759.         }
  760.     }

  762.     /**
  763.      * Get the non-resonant tesseral terms in the central body spherical harmonic field.
  764.      *
  765.      * @param <T> type of the elements
  766.      * @param type type of the elements used during the propagation
  767.      * @param context container for attributes
  768.      * @param field field used by default
  769.      */
  770.     private <T extends RealFieldElement<T>> void getNonResonantTerms(final PropagationType type,
  771.                                                                      final FieldDSSTTesseralContext<T> context,
  772.                                                                      final Field<T> field) {

  774.         final T zero = field.getZero();
  775.         // Compute natural resonant terms
  776.         final T tolerance = FastMath.max(zero.add(MIN_PERIOD_IN_SAT_REV),
  777.                                          context.getOrbitPeriod().divide(MIN_PERIOD_IN_SECONDS).reciprocal()).reciprocal();

  779.         nonResOrders.clear();
  780.         for (int m = 1; m <= maxOrder; m++) {
  781.             final T resonance = context.getRatio().multiply(m);
  782.             int jRes = 0;
  783.             final int jComputedRes = (int) FastMath.round(resonance);
  784.             if (jComputedRes > 0 && jComputedRes <= maxFrequencyShortPeriodics && FastMath.abs(resonance.subtract(jComputedRes)).getReal() <= tolerance.getReal()) {
  785.                 // Store each resonant index and order
  786.                 jRes = jComputedRes;
  787.             }

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

  799.                 nonResOrders.put(m, listJofM);
  800.             }
  801.         }
  802.     }

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

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

  825.         //spherical harmonics
  826.         final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);

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

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

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

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

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

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

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

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

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

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

  888.                 // Jacobi l-index from 2.7.1-(15)
  889.                 final int l = FastMath.min(n - m, n - FastMath.abs(s));
  890.                 // Jacobi polynomial and derivative
  891.                 final Gradient jacobi =
  892.                         JacobiPolynomials.getValue(l, v, w, Gradient.variable(1, 0, auxiliaryElements.getGamma()));

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

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

  910.                 // dU / da components
  911.                 dUdaCos  += dUdaCoef * gcPhs;
  912.                 dUdaSin  += dUdaCoef * gsMhc;

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

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

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

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

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

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

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

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

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

  977.         //spherical harmonics
  978.         final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date.toAbsoluteDate());

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  1108.         return derivatives;

  1110.     }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  1618.             // dλ/dt
  1619.             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);
  1620.             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);
  1621.         }

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

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

  1644.     /** 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
  1645.      * the short-periodic zonal contribution.
  1646.      *   <p>
  1647.      *  Those coefficients are given by expression 2.5.4-5 from the Danielson paper.
  1648.      *   </p>
  1649.      *
  1650.      * @author Sorin Scortan
  1651.      *
  1652.      */
  1653.     private static class TesseralShortPeriodicCoefficients implements ShortPeriodTerms {

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

  1670.         /** Central body rotating frame. */
  1671.         private final Frame bodyFrame;

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

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

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

  1682.         /** Maximum value for m index. */
  1683.         private final int mMax;

  1685.         /** Maximum value for j index. */
  1686.         private final int jMax;

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

  1691.         /** All coefficients slots. */
  1692.         private final transient TimeSpanMap<Slot> slots;

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

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

  1734.         /** {@inheritDoc} */
  1735.         @Override
  1736.         public double[] value(final Orbit meanOrbit) {

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

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

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

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

  1751.                 // Central body rotation angle from equation 2.7.1-(3)(4).
  1752.                 final Transform t = bodyFrame.getTransformTo(auxiliaryElements.getFrame(), auxiliaryElements.getDate());
  1753.                 final Vector3D xB = t.transformVector(Vector3D.PLUS_I);
  1754.                 final Vector3D yB = t.transformVector(Vector3D.PLUS_J);
  1755.                 final Vector3D  f = auxiliaryElements.getVectorF();
  1756.                 final Vector3D  g = auxiliaryElements.getVectorG();
  1757.                 final double currentTheta = FastMath.atan2(-f.dotProduct(yB) + I * g.dotProduct(xB),
  1758.                                                             f.dotProduct(xB) + I * g.dotProduct(yB));

  1760.                 //Add the m-daily contribution
  1761.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1762.                     // Phase angle
  1763.                     final double jlMmt  = -m * currentTheta;
  1764.                     final SinCos scPhi  = FastMath.sinCos(jlMmt);
  1765.                     final double sinPhi = scPhi.sin();
  1766.                     final double cosPhi = scPhi.cos();

  1768.                     // compute contribution for each element
  1769.                     final double[] c = slot.getCijm(0, m, meanOrbit.getDate());
  1770.                     final double[] s = slot.getSijm(0, m, meanOrbit.getDate());
  1771.                     for (int i = 0; i < 6; i++) {
  1772.                         shortPeriodicVariation[i] += c[i] * cosPhi + s[i] * sinPhi;
  1773.                     }
  1774.                 }

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

  1781.                     for (int j : listJ) {
  1782.                         // Phase angle
  1783.                         final double jlMmt  = j * meanOrbit.getLM() - m * currentTheta;
  1784.                         final SinCos scPhi  = FastMath.sinCos(jlMmt);
  1785.                         final double sinPhi = scPhi.sin();
  1786.                         final double cosPhi = scPhi.cos();

  1788.                         // compute contribution for each element
  1789.                         final double[] c = slot.getCijm(j, m, meanOrbit.getDate());
  1790.                         final double[] s = slot.getSijm(j, m, meanOrbit.getDate());
  1791.                         for (int i = 0; i < 6; i++) {
  1792.                             shortPeriodicVariation[i] += c[i] * cosPhi + s[i] * sinPhi;
  1793.                         }

  1795.                     }
  1796.                 }
  1797.             }

  1799.             return shortPeriodicVariation;

  1801.         }

  1803.         /** {@inheritDoc} */
  1804.         @Override
  1805.         public String getCoefficientsKeyPrefix() {
  1806.             return DSSTTesseral.SHORT_PERIOD_PREFIX;
  1807.         }

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

  1822.             // select the coefficients slot
  1823.             final Slot slot = slots.get(date);

  1825.             if (!nonResOrders.isEmpty() || mDailiesOnly) {
  1826.                 final Map<String, double[]> coefficients =
  1827.                                 new HashMap<String, double[]>(12 * maxOrderMdailyTesseralSP +
  1828.                                                               12 * nonResOrders.size());

  1830.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1831.                     storeIfSelected(coefficients, selected, slot.getCijm(0, m, date), DSSTTesseral.CM_COEFFICIENTS, m);
  1832.                     storeIfSelected(coefficients, selected, slot.getSijm(0, m, date), DSSTTesseral.SM_COEFFICIENTS, m);
  1833.                 }

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

  1839.                     for (int j : listJ) {
  1840.                         for (int i = 0; i < 6; ++i) {
  1841.                             storeIfSelected(coefficients, selected, slot.getCijm(j, m, date), "c", j, m);
  1842.                             storeIfSelected(coefficients, selected, slot.getSijm(j, m, date), "s", j, m);
  1843.                         }
  1844.                     }
  1845.                 }

  1847.                 return coefficients;

  1849.             } else {
  1850.                 return Collections.emptyMap();
  1851.             }

  1853.         }

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

  1876.     }

  1878.     /** 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
  1879.      * the short-periodic zonal contribution.
  1880.      *   <p>
  1881.      *  Those coefficients are given by expression 2.5.4-5 from the Danielson paper.
  1882.      *   </p>
  1883.      *
  1884.      * @author Sorin Scortan
  1885.      *
  1886.      */
  1887.     private static class FieldTesseralShortPeriodicCoefficients <T extends RealFieldElement<T>> implements FieldShortPeriodTerms<T> {

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

  1904.         /** Central body rotating frame. */
  1905.         private final Frame bodyFrame;

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

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

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

  1916.         /** Maximum value for m index. */
  1917.         private final int mMax;

  1919.         /** Maximum value for j index. */
  1920.         private final int jMax;

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

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

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

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

  1969.         /** {@inheritDoc} */
  1970.         @Override
  1971.         public T[] value(final FieldOrbit<T> meanOrbit) {

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

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

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

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

  1986.                 // Central body rotation angle from equation 2.7.1-(3)(4).
  1987.                 final FieldTransform<T> t = bodyFrame.getTransformTo(auxiliaryElements.getFrame(), auxiliaryElements.getDate());
  1988.                 final FieldVector3D<T> xB = t.transformVector(Vector3D.PLUS_I);
  1989.                 final FieldVector3D<T> yB = t.transformVector(Vector3D.PLUS_J);
  1990.                 final FieldVector3D<T>  f = auxiliaryElements.getVectorF();
  1991.                 final FieldVector3D<T>  g = auxiliaryElements.getVectorG();
  1992.                 final T currentTheta = FastMath.atan2(f.dotProduct(yB).negate().add(g.dotProduct(xB).multiply(I)),
  1993.                                                       f.dotProduct(xB).add(g.dotProduct(yB).multiply(I)));

  1995.                 //Add the m-daily contribution
  1996.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1997.                     // Phase angle
  1998.                     final T jlMmt              = currentTheta.multiply(-m);
  1999.                     final FieldSinCos<T> scPhi = FastMath.sinCos(jlMmt);
  2000.                     final T sinPhi             = scPhi.sin();
  2001.                     final T cosPhi             = scPhi.cos();

  2003.                     // compute contribution for each element
  2004.                     final T[] c = slot.getCijm(0, m, meanOrbit.getDate());
  2005.                     final T[] s = slot.getSijm(0, m, meanOrbit.getDate());
  2006.                     for (int i = 0; i < 6; i++) {
  2007.                         shortPeriodicVariation[i] = shortPeriodicVariation[i].add(c[i].multiply(cosPhi).add(s[i].multiply(sinPhi)));
  2008.                     }
  2009.                 }

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

  2016.                     for (int j : listJ) {
  2017.                         // Phase angle
  2018.                         final T jlMmt              = meanOrbit.getLM().multiply(j).subtract(currentTheta.multiply(m));
  2019.                         final FieldSinCos<T> scPhi = FastMath.sinCos(jlMmt);
  2020.                         final T sinPhi             = scPhi.sin();
  2021.                         final T cosPhi             = scPhi.cos();

  2023.                         // compute contribution for each element
  2024.                         final T[] c = slot.getCijm(j, m, meanOrbit.getDate());
  2025.                         final T[] s = slot.getSijm(j, m, meanOrbit.getDate());
  2026.                         for (int i = 0; i < 6; i++) {
  2027.                             shortPeriodicVariation[i] = shortPeriodicVariation[i].add(c[i].multiply(cosPhi).add(s[i].multiply(sinPhi)));
  2028.                         }

  2030.                     }
  2031.                 }
  2032.             }

  2034.             return shortPeriodicVariation;

  2036.         }

  2038.         /** {@inheritDoc} */
  2039.         @Override
  2040.         public String getCoefficientsKeyPrefix() {
  2041.             return DSSTTesseral.SHORT_PERIOD_PREFIX;
  2042.         }

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

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

  2060.             if (!nonResOrders.isEmpty() || mDailiesOnly) {
  2061.                 final Map<String, T[]> coefficients =
  2062.                                 new HashMap<String, T[]>(12 * maxOrderMdailyTesseralSP +
  2063.                                                          12 * nonResOrders.size());

  2065.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  2066.                     storeIfSelected(coefficients, selected, slot.getCijm(0, m, date), DSSTTesseral.CM_COEFFICIENTS, m);
  2067.                     storeIfSelected(coefficients, selected, slot.getSijm(0, m, date), DSSTTesseral.SM_COEFFICIENTS, m);
  2068.                 }

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

  2074.                     for (int j : listJ) {
  2075.                         for (int i = 0; i < 6; ++i) {
  2076.                             storeIfSelected(coefficients, selected, slot.getCijm(j, m, date), "c", j, m);
  2077.                             storeIfSelected(coefficients, selected, slot.getSijm(j, m, date), "s", j, m);
  2078.                         }
  2079.                     }
  2080.                 }

  2082.                 return coefficients;

  2084.             } else {
  2085.                 return Collections.emptyMap();
  2086.             }

  2088.         }

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

  2112.     /** Coefficients valid for one time slot. */
  2113.     private static class Slot {

  2115.         /** The coefficients C<sub>i</sub><sup>j</sup><sup>m</sup>.
  2116.          * <p>
  2117.          * The index order is cijm[m][j][i] <br/>
  2118.          * i corresponds to the equinoctial element, as follows: <br/>
  2119.          * - i=0 for a <br/>
  2120.          * - i=1 for k <br/>
  2121.          * - i=2 for h <br/>
  2122.          * - i=3 for q <br/>
  2123.          * - i=4 for p <br/>
  2124.          * - i=5 for λ <br/>
  2125.          * </p>
  2126.          */
  2127.         private final ShortPeriodicsInterpolatedCoefficient[][] cijm;

  2129.         /** The coefficients S<sub>i</sub><sup>j</sup><sup>m</sup>.
  2130.          * <p>
  2131.          * The index order is sijm[m][j][i] <br/>
  2132.          * i corresponds to the equinoctial element, as follows: <br/>
  2133.          * - i=0 for a <br/>
  2134.          * - i=1 for k <br/>
  2135.          * - i=2 for h <br/>
  2136.          * - i=3 for q <br/>
  2137.          * - i=4 for p <br/>
  2138.          * - i=5 for λ <br/>
  2139.          * </p>
  2140.          */
  2141.         private final ShortPeriodicsInterpolatedCoefficient[][] sijm;

  2143.         /** Simple constructor.
  2144.          *  @param mMax maximum value for m index
  2145.          *  @param jMax maximum value for j index
  2146.          *  @param interpolationPoints number of points used in the interpolation process
  2147.          */
  2148.         Slot(final int mMax, final int jMax, final int interpolationPoints) {

  2150.             final int rows    = mMax + 1;
  2151.             final int columns = 2 * jMax + 1;
  2152.             cijm = new ShortPeriodicsInterpolatedCoefficient[rows][columns];
  2153.             sijm = new ShortPeriodicsInterpolatedCoefficient[rows][columns];
  2154.             for (int m = 1; m <= mMax; m++) {
  2155.                 for (int j = -jMax; j <= jMax; j++) {
  2156.                     cijm[m][j + jMax] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  2157.                     sijm[m][j + jMax] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  2158.                 }
  2159.             }

  2161.         }

  2163.         /** Get C<sub>i</sub><sup>j</sup><sup>m</sup>.
  2164.          *
  2165.          * @param j j index
  2166.          * @param m m index
  2167.          * @param date the date
  2168.          * @return C<sub>i</sub><sup>j</sup><sup>m</sup>
  2169.          */
  2170.         double[] getCijm(final int j, final int m, final AbsoluteDate date) {
  2171.             final int jMax = (cijm[m].length - 1) / 2;
  2172.             return cijm[m][j + jMax].value(date);
  2173.         }

  2175.         /** Get S<sub>i</sub><sup>j</sup><sup>m</sup>.
  2176.          *
  2177.          * @param j j index
  2178.          * @param m m index
  2179.          * @param date the date
  2180.          * @return S<sub>i</sub><sup>j</sup><sup>m</sup>
  2181.          */
  2182.         double[] getSijm(final int j, final int m, final AbsoluteDate date) {
  2183.             final int jMax = (cijm[m].length - 1) / 2;
  2184.             return sijm[m][j + jMax].value(date);
  2185.         }

  2187.     }

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

  2192.         /** The coefficients C<sub>i</sub><sup>j</sup><sup>m</sup>.
  2193.          * <p>
  2194.          * The index order is cijm[m][j][i] <br/>
  2195.          * i corresponds to the equinoctial element, as follows: <br/>
  2196.          * - i=0 for a <br/>
  2197.          * - i=1 for k <br/>
  2198.          * - i=2 for h <br/>
  2199.          * - i=3 for q <br/>
  2200.          * - i=4 for p <br/>
  2201.          * - i=5 for λ <br/>
  2202.          * </p>
  2203.          */
  2204.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[][] cijm;

  2206.         /** The coefficients S<sub>i</sub><sup>j</sup><sup>m</sup>.
  2207.          * <p>
  2208.          * The index order is sijm[m][j][i] <br/>
  2209.          * i corresponds to the equinoctial element, as follows: <br/>
  2210.          * - i=0 for a <br/>
  2211.          * - i=1 for k <br/>
  2212.          * - i=2 for h <br/>
  2213.          * - i=3 for q <br/>
  2214.          * - i=4 for p <br/>
  2215.          * - i=5 for λ <br/>
  2216.          * </p>
  2217.          */
  2218.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[][] sijm;

  2220.         /** Simple constructor.
  2221.          *  @param mMax maximum value for m index
  2222.          *  @param jMax maximum value for j index
  2223.          *  @param interpolationPoints number of points used in the interpolation process
  2224.          */
  2225.         @SuppressWarnings("unchecked")
  2226.         FieldSlot(final int mMax, final int jMax, final int interpolationPoints) {

  2228.             final int rows    = mMax + 1;
  2229.             final int columns = 2 * jMax + 1;
  2230.             cijm = (FieldShortPeriodicsInterpolatedCoefficient<T>[][]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows, columns);
  2231.             sijm = (FieldShortPeriodicsInterpolatedCoefficient<T>[][]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows, columns);
  2232.             for (int m = 1; m <= mMax; m++) {
  2233.                 for (int j = -jMax; j <= jMax; j++) {
  2234.                     cijm[m][j + jMax] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  2235.                     sijm[m][j + jMax] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  2236.                 }
  2237.             }

  2239.         }

  2241.         /** Get C<sub>i</sub><sup>j</sup><sup>m</sup>.
  2242.          *
  2243.          * @param j j index
  2244.          * @param m m index
  2245.          * @param date the date
  2246.          * @return C<sub>i</sub><sup>j</sup><sup>m</sup>
  2247.          */
  2248.         T[] getCijm(final int j, final int m, final FieldAbsoluteDate<T> date) {
  2249.             final int jMax = (cijm[m].length - 1) / 2;
  2250.             return cijm[m][j + jMax].value(date);
  2251.         }

  2253.         /** Get S<sub>i</sub><sup>j</sup><sup>m</sup>.
  2254.          *
  2255.          * @param j j index
  2256.          * @param m m index
  2257.          * @param date the date
  2258.          * @return S<sub>i</sub><sup>j</sup><sup>m</sup>
  2259.          */
  2260.         T[] getSijm(final int j, final int m, final FieldAbsoluteDate<T> date) {
  2261.             final int jMax = (cijm[m].length - 1) / 2;
  2262.             return sijm[m][j + jMax].value(date);
  2263.         }

  2265.     }

  2267.     /** Compute potential and potential derivatives with respect to orbital parameters.
  2268.      *  <p>The following elements are computed from expression 3.3 - (4).
  2269.      *  <pre>
  2270.      *  dU / da
  2271.      *  dU / dh
  2272.      *  dU / dk
  2273.      *  dU / dλ
  2274.      *  dU / dα
  2275.      *  dU / dβ
  2276.      *  dU / dγ
  2277.      *  </pre>
  2278.      *  </p>
  2279.      */
  2280.     private class UAnddU {

  2282.         /** dU / da. */
  2283.         private  double dUda;

  2285.         /** dU / dk. */
  2286.         private double dUdk;

  2288.         /** dU / dh. */
  2289.         private double dUdh;

  2291.         /** dU / dl. */
  2292.         private double dUdl;

  2294.         /** dU / dAlpha. */
  2295.         private double dUdAl;

  2297.         /** dU / dBeta. */
  2298.         private double dUdBe;

  2300.         /** dU / dGamma. */
  2301.         private double dUdGa;

  2303.         /** Simple constuctor.
  2304.          * @param date current date
  2305.          * @param context container for attributes
  2306.          * @param hansen hansen objects
  2307.          */
  2308.         UAnddU(final AbsoluteDate date, final DSSTTesseralContext context, final HansenObjects hansen) {

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

  2313.             // Potential derivatives
  2314.             dUda  = 0.;
  2315.             dUdh  = 0.;
  2316.             dUdk  = 0.;
  2317.             dUdl  = 0.;
  2318.             dUdAl = 0.;
  2319.             dUdBe = 0.;
  2320.             dUdGa = 0.;

  2322.             // Compute only if there is at least one resonant tesseral
  2323.             if (!context.getResOrders().isEmpty()) {
  2324.                 // Gmsj and Hmsj polynomials
  2325.                 final GHmsjPolynomials ghMSJ = new GHmsjPolynomials(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta(), I);

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

  2330.                 // R / a up to power degree
  2331.                 final double[] roaPow = new double[maxDegree + 1];
  2332.                 roaPow[0] = 1.;
  2333.                 for (int i = 1; i <= maxDegree; i++) {
  2334.                     roaPow[i] = context.getRoa() * roaPow[i - 1];
  2335.                 }

  2337.                 // SUM over resonant terms {j,m}
  2338.                 for (int m : context.getResOrders()) {

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

  2343.                     // Phase angle
  2344.                     final double jlMmt  = j * auxiliaryElements.getLM() - m * context.getTheta();
  2345.                     final SinCos scPhi  = FastMath.sinCos(jlMmt);
  2346.                     final double sinPhi = scPhi.sin();
  2347.                     final double cosPhi = scPhi.cos();

  2349.                     // Potential derivatives components for a given resonant pair {j,m}
  2350.                     double dUdaCos  = 0.;
  2351.                     double dUdaSin  = 0.;
  2352.                     double dUdhCos  = 0.;
  2353.                     double dUdhSin  = 0.;
  2354.                     double dUdkCos  = 0.;
  2355.                     double dUdkSin  = 0.;
  2356.                     double dUdlCos  = 0.;
  2357.                     double dUdlSin  = 0.;
  2358.                     double dUdAlCos = 0.;
  2359.                     double dUdAlSin = 0.;
  2360.                     double dUdBeCos = 0.;
  2361.                     double dUdBeSin = 0.;
  2362.                     double dUdGaCos = 0.;
  2363.                     double dUdGaSin = 0.;

  2365.                     // s-SUM from -sMin to sMax
  2366.                     final int sMin = FastMath.min(context.getMaxEccPow() - j, maxDegree);
  2367.                     final int sMax = FastMath.min(context.getMaxEccPow() + j, maxDegree);
  2368.                     for (int s = 0; s <= sMax; s++) {

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

  2373.                         // n-SUM for s positive
  2374.                         final double[][] nSumSpos = computeNSum(date, j, m, s, maxDegree,
  2375.                                                                 roaPow, ghMSJ, gammaMNS, context, hansen);
  2376.                         dUdaCos  += nSumSpos[0][0];
  2377.                         dUdaSin  += nSumSpos[0][1];
  2378.                         dUdhCos  += nSumSpos[1][0];
  2379.                         dUdhSin  += nSumSpos[1][1];
  2380.                         dUdkCos  += nSumSpos[2][0];
  2381.                         dUdkSin  += nSumSpos[2][1];
  2382.                         dUdlCos  += nSumSpos[3][0];
  2383.                         dUdlSin  += nSumSpos[3][1];
  2384.                         dUdAlCos += nSumSpos[4][0];
  2385.                         dUdAlSin += nSumSpos[4][1];
  2386.                         dUdBeCos += nSumSpos[5][0];
  2387.                         dUdBeSin += nSumSpos[5][1];
  2388.                         dUdGaCos += nSumSpos[6][0];
  2389.                         dUdGaSin += nSumSpos[6][1];

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

  2396.                             final double[][] nSumSneg = computeNSum(date, j, m, -s, maxDegree,
  2397.                                                                     roaPow, ghMSJ, gammaMNS, context, hansen);
  2398.                             dUdaCos  += nSumSneg[0][0];
  2399.                             dUdaSin  += nSumSneg[0][1];
  2400.                             dUdhCos  += nSumSneg[1][0];
  2401.                             dUdhSin  += nSumSneg[1][1];
  2402.                             dUdkCos  += nSumSneg[2][0];
  2403.                             dUdkSin  += nSumSneg[2][1];
  2404.                             dUdlCos  += nSumSneg[3][0];
  2405.                             dUdlSin  += nSumSneg[3][1];
  2406.                             dUdAlCos += nSumSneg[4][0];
  2407.                             dUdAlSin += nSumSneg[4][1];
  2408.                             dUdBeCos += nSumSneg[5][0];
  2409.                             dUdBeSin += nSumSneg[5][1];
  2410.                             dUdGaCos += nSumSneg[6][0];
  2411.                             dUdGaSin += nSumSneg[6][1];
  2412.                         }
  2413.                     }

  2415.                     // Assembly of potential derivatives componants
  2416.                     dUda  += cosPhi * dUdaCos  + sinPhi * dUdaSin;
  2417.                     dUdh  += cosPhi * dUdhCos  + sinPhi * dUdhSin;
  2418.                     dUdk  += cosPhi * dUdkCos  + sinPhi * dUdkSin;
  2419.                     dUdl  += cosPhi * dUdlCos  + sinPhi * dUdlSin;
  2420.                     dUdAl += cosPhi * dUdAlCos + sinPhi * dUdAlSin;
  2421.                     dUdBe += cosPhi * dUdBeCos + sinPhi * dUdBeSin;
  2422.                     dUdGa += cosPhi * dUdGaCos + sinPhi * dUdGaSin;
  2423.                 }

  2425.                 this.dUda  = dUda * (-context.getMoa() / auxiliaryElements.getSma());
  2426.                 this.dUdh  = dUdh * context.getMoa();
  2427.                 this.dUdk  = dUdk * context.getMoa();
  2428.                 this.dUdl  = dUdl * context.getMoa();
  2429.                 this.dUdAl = dUdAl * context.getMoa();
  2430.                 this.dUdBe = dUdBe * context.getMoa();
  2431.                 this.dUdGa = dUdGa * context.getMoa();
  2432.             }

  2434.         }

  2436.         /** Return value of dU / da.
  2437.          * @return dUda
  2438.          */
  2439.         public double getdUda() {
  2440.             return dUda;
  2441.         }

  2443.         /** Return value of dU / dk.
  2444.          * @return dUdk
  2445.          */
  2446.         public double getdUdk() {
  2447.             return dUdk;
  2448.         }

  2450.         /** Return value of dU / dh.
  2451.          * @return dUdh
  2452.          */
  2453.         public double getdUdh() {
  2454.             return dUdh;
  2455.         }

  2457.         /** Return value of dU / dl.
  2458.          * @return dUdl
  2459.          */
  2460.         public double getdUdl() {
  2461.             return dUdl;
  2462.         }

  2464.         /** Return value of dU / dAlpha.
  2465.          * @return dUdAl
  2466.          */
  2467.         public double getdUdAl() {
  2468.             return dUdAl;
  2469.         }

  2471.         /** Return value of dU / dBeta.
  2472.          * @return dUdBe
  2473.          */
  2474.         public double getdUdBe() {
  2475.             return dUdBe;
  2476.         }

  2478.         /** Return value of dU / dGamma.
  2479.          * @return dUdGa
  2480.          */
  2481.         public double getdUdGa() {
  2482.             return dUdGa;
  2483.         }

  2485.     }

  2487.     /**  Computes the potential U derivatives.
  2488.      *  <p>The following elements are computed from expression 3.3 - (4).
  2489.      *  <pre>
  2490.      *  dU / da
  2491.      *  dU / dh
  2492.      *  dU / dk
  2493.      *  dU / dλ
  2494.      *  dU / dα
  2495.      *  dU / dβ
  2496.      *  dU / dγ
  2497.      *  </pre>
  2498.      *  </p>
  2499.      */
  2500.     private class FieldUAnddU <T extends RealFieldElement<T>> {

  2502.         /** dU / da. */
  2503.         private T dUda;

  2505.         /** dU / dk. */
  2506.         private T dUdk;

  2508.         /** dU / dh. */
  2509.         private T dUdh;

  2511.         /** dU / dl. */
  2512.         private T dUdl;

  2514.         /** dU / dAlpha. */
  2515.         private T dUdAl;

  2517.         /** dU / dBeta. */
  2518.         private T dUdBe;

  2520.         /** dU / dGamma. */
  2521.         private T dUdGa;

  2523.         /** Simple constuctor.
  2524.          * @param date current date
  2525.          * @param context container for attributes
  2526.          * @param hansen hansen objects
  2527.          */
  2528.         FieldUAnddU(final FieldAbsoluteDate<T> date, final FieldDSSTTesseralContext<T> context,
  2529.                     final FieldHansenObjects<T> hansen) {

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

  2534.             // Zero for initialization
  2535.             final Field<T> field = date.getField();
  2536.             final T zero = field.getZero();

  2538.             // Potential derivatives
  2539.             dUda  = zero;
  2540.             dUdh  = zero;
  2541.             dUdk  = zero;
  2542.             dUdl  = zero;
  2543.             dUdAl = zero;
  2544.             dUdBe = zero;
  2545.             dUdGa = zero;

  2547.             // Compute only if there is at least one resonant tesseral
  2548.             if (!context.getResOrders().isEmpty()) {
  2549.                 // Gmsj and Hmsj polynomials
  2550.                 final FieldGHmsjPolynomials<T> ghMSJ = new FieldGHmsjPolynomials<>(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta(), I, field);

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

  2555.                 // R / a up to power degree
  2556.                 final T[] roaPow = MathArrays.buildArray(field, maxDegree + 1);
  2557.                 roaPow[0] = zero.add(1.);
  2558.                 for (int i = 1; i <= maxDegree; i++) {
  2559.                     roaPow[i] = roaPow[i - 1].multiply(context.getRoa());
  2560.                 }

  2562.                 // SUM over resonant terms {j,m}
  2563.                 for (int m : context.getResOrders()) {

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

  2568.                     // Phase angle
  2569.                     final T jlMmt              = auxiliaryElements.getLM().multiply(j).subtract(context.getTheta().multiply(m));
  2570.                     final FieldSinCos<T> scPhi = FastMath.sinCos(jlMmt);
  2571.                     final T sinPhi             = scPhi.sin();
  2572.                     final T cosPhi             = scPhi.cos();

  2574.                     // Potential derivatives components for a given resonant pair {j,m}
  2575.                     T dUdaCos  = zero;
  2576.                     T dUdaSin  = zero;
  2577.                     T dUdhCos  = zero;
  2578.                     T dUdhSin  = zero;
  2579.                     T dUdkCos  = zero;
  2580.                     T dUdkSin  = zero;
  2581.                     T dUdlCos  = zero;
  2582.                     T dUdlSin  = zero;
  2583.                     T dUdAlCos = zero;
  2584.                     T dUdAlSin = zero;
  2585.                     T dUdBeCos = zero;
  2586.                     T dUdBeSin = zero;
  2587.                     T dUdGaCos = zero;
  2588.                     T dUdGaSin = zero;

  2590.                     // s-SUM from -sMin to sMax
  2591.                     final int sMin = FastMath.min(context.getMaxEccPow() - j, maxDegree);
  2592.                     final int sMax = FastMath.min(context.getMaxEccPow() + j, maxDegree);
  2593.                     for (int s = 0; s <= sMax; s++) {

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

  2598.                         // n-SUM for s positive
  2599.                         final T[][] nSumSpos = computeNSum(date, j, m, s, maxDegree,
  2600.                                                                 roaPow, ghMSJ, gammaMNS, context, hansen);
  2601.                         dUdaCos  = dUdaCos.add(nSumSpos[0][0]);
  2602.                         dUdaSin  = dUdaSin.add(nSumSpos[0][1]);
  2603.                         dUdhCos  = dUdhCos.add(nSumSpos[1][0]);
  2604.                         dUdhSin  = dUdhSin.add(nSumSpos[1][1]);
  2605.                         dUdkCos  = dUdkCos.add(nSumSpos[2][0]);
  2606.                         dUdkSin  = dUdkSin.add(nSumSpos[2][1]);
  2607.                         dUdlCos  = dUdlCos.add(nSumSpos[3][0]);
  2608.                         dUdlSin  = dUdlSin.add(nSumSpos[3][1]);
  2609.                         dUdAlCos = dUdAlCos.add(nSumSpos[4][0]);
  2610.                         dUdAlSin = dUdAlSin.add(nSumSpos[4][1]);
  2611.                         dUdBeCos = dUdBeCos.add(nSumSpos[5][0]);
  2612.                         dUdBeSin = dUdBeSin.add(nSumSpos[5][1]);
  2613.                         dUdGaCos = dUdGaCos.add(nSumSpos[6][0]);
  2614.                         dUdGaSin = dUdGaSin.add(nSumSpos[6][1]);

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

  2621.                             final T[][] nSumSneg = computeNSum(date, j, m, -s, maxDegree,
  2622.                                                                     roaPow, ghMSJ, gammaMNS, context, hansen);
  2623.                             dUdaCos  = dUdaCos.add(nSumSneg[0][0]);
  2624.                             dUdaSin  = dUdaSin.add(nSumSneg[0][1]);
  2625.                             dUdhCos  = dUdhCos.add(nSumSneg[1][0]);
  2626.                             dUdhSin  = dUdhSin.add(nSumSneg[1][1]);
  2627.                             dUdkCos  = dUdkCos.add(nSumSneg[2][0]);
  2628.                             dUdkSin  = dUdkSin.add(nSumSneg[2][1]);
  2629.                             dUdlCos  = dUdlCos.add(nSumSneg[3][0]);
  2630.                             dUdlSin  = dUdlSin.add(nSumSneg[3][1]);
  2631.                             dUdAlCos = dUdAlCos.add(nSumSneg[4][0]);
  2632.                             dUdAlSin = dUdAlSin.add(nSumSneg[4][1]);
  2633.                             dUdBeCos = dUdBeCos.add(nSumSneg[5][0]);
  2634.                             dUdBeSin = dUdBeSin.add(nSumSneg[5][1]);
  2635.                             dUdGaCos = dUdGaCos.add(nSumSneg[6][0]);
  2636.                             dUdGaSin = dUdGaSin.add(nSumSneg[6][1]);
  2637.                         }
  2638.                     }

  2640.                     // Assembly of potential derivatives componants
  2641.                     dUda  = dUda.add(dUdaCos.multiply(cosPhi).add(dUdaSin.multiply(sinPhi)));
  2642.                     dUdh  = dUdh.add(dUdhCos.multiply(cosPhi).add(dUdhSin.multiply(sinPhi)));
  2643.                     dUdk  = dUdk.add(dUdkCos.multiply(cosPhi).add(dUdkSin.multiply(sinPhi)));
  2644.                     dUdl  = dUdl.add(dUdlCos.multiply(cosPhi).add(dUdlSin.multiply(sinPhi)));
  2645.                     dUdAl = dUdAl.add(dUdAlCos.multiply(cosPhi).add(dUdAlSin.multiply(sinPhi)));
  2646.                     dUdBe = dUdBe.add(dUdBeCos.multiply(cosPhi).add(dUdBeSin.multiply(sinPhi)));
  2647.                     dUdGa = dUdGa.add(dUdGaCos.multiply(cosPhi).add(dUdGaSin.multiply(sinPhi)));
  2648.                 }

  2650.                 dUda  =  dUda.multiply(context.getMoa().divide(auxiliaryElements.getSma())).negate();
  2651.                 dUdh  =  dUdh.multiply(context.getMoa());
  2652.                 dUdk  =  dUdk.multiply(context.getMoa());
  2653.                 dUdl  =  dUdl.multiply(context.getMoa());
  2654.                 dUdAl =  dUdAl.multiply(context.getMoa());
  2655.                 dUdBe =  dUdBe.multiply(context.getMoa());
  2656.                 dUdGa =  dUdGa.multiply(context.getMoa());

  2658.             }

  2660.         }

  2662.         /** Return value of dU / da.
  2663.          * @return dUda
  2664.          */
  2665.         public T getdUda() {
  2666.             return dUda;
  2667.         }

  2669.         /** Return value of dU / dk.
  2670.          * @return dUdk
  2671.          */
  2672.         public T getdUdk() {
  2673.             return dUdk;
  2674.         }

  2676.         /** Return value of dU / dh.
  2677.          * @return dUdh
  2678.          */
  2679.         public T getdUdh() {
  2680.             return dUdh;
  2681.         }

  2683.         /** Return value of dU / dl.
  2684.          * @return dUdl
  2685.          */
  2686.         public T getdUdl() {
  2687.             return dUdl;
  2688.         }

  2690.         /** Return value of dU / dAlpha.
  2691.          * @return dUdAl
  2692.          */
  2693.         public T getdUdAl() {
  2694.             return dUdAl;
  2695.         }

  2697.         /** Return value of dU / dBeta.
  2698.          * @return dUdBe
  2699.          */
  2700.         public T getdUdBe() {
  2701.             return dUdBe;
  2702.         }

  2704.         /** Return value of dU / dGamma.
  2705.          * @return dUdGa
  2706.          */
  2707.         public T getdUdGa() {
  2708.             return dUdGa;
  2709.         }

  2711.     }

  2713.     /** Computes init values of the Hansen Objects. */
  2714.     private class HansenObjects {

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

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

  2723.         /** Simple constructor.
  2724.          * @param ratio Ratio of satellite period to central body rotation period
  2725.          * @param maxEccPow Maximum power of the eccentricity to use in summation over s.
  2726.          * @param resOrders List of resonant orders
  2727.          * @param type type of the elements used during the propagation
  2728.          */
  2729.         HansenObjects(final double ratio,
  2730.                       final int maxEccPow,
  2731.                       final List<Integer> resOrders,
  2732.                       final PropagationType type) {

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

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

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

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

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

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

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

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

  2780.                 default:
  2781.                     throw new OrekitInternalError(null);
  2782.             }

  2784.         }

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

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

  2802.     }

  2804.     /** Computes init values of the Hansen Objects. */
  2805.     private class FieldHansenObjects<T extends RealFieldElement<T>> {

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

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

  2814.         /** Simple constructor.
  2815.          * @param ratio Ratio of satellite period to central body rotation period
  2816.          * @param maxEccPow Maximum power of the eccentricity to use in summation over s.
  2817.          * @param resOrders List of resonant orders
  2818.          * @param type type of the elements used during the propagation
  2819.          * @param field field used by default.
  2820.          */
  2821.         @SuppressWarnings("unchecked")
  2822.         FieldHansenObjects(final T ratio,
  2823.                            final int maxEccPow,
  2824.                            final List<Integer> resOrders,
  2825.                            final PropagationType type,
  2826.                            final Field<T> field) {

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

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

  2836.             switch (type) {
  2837.                 case MEAN:
  2838.                  // loop through the resonant orders
  2839.                     for (int m : resOrders) {
  2840.                         //Compute the corresponding j term
  2841.                         final int j = FastMath.max(1, (int) FastMath.round(ratio.multiply(m)));

  2843.                         //Compute the sMin and sMax values
  2844.                         final int sMin = FastMath.min(maxEccPow - j, maxDegree);
  2845.                         final int sMax = FastMath.min(maxEccPow + j, maxDegree);

  2847.                         //loop through the s values
  2848.                         for (int s = 0; s <= sMax; s++) {
  2849.                             //Compute the n0 value
  2850.                             final int n0 = FastMath.max(FastMath.max(2, m), s);

  2852.                             //Create the object for the pair j, s
  2853.                             this.hansenObjects[s + maxDegree][j] = new FieldHansenTesseralLinear<>(maxDegree, s, j, n0, maxHansen, field);

  2855.                             if (s > 0 && s <= sMin) {
  2856.                                 //Also create the object for the pair j, -s
  2857.                                 this.hansenObjects[maxDegree - s][j] =  new FieldHansenTesseralLinear<>(maxDegree, -s, j, n0, maxHansen, field);
  2858.                             }
  2859.                         }
  2860.                     }
  2861.                     break;

  2863.                 case OSCULATING:
  2864.                     // create all objects
  2865.                     for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
  2866.                         for (int s = -maxDegree; s <= maxDegree; s++) {
  2867.                             //Compute the n0 value
  2868.                             final int n0 = FastMath.max(2, FastMath.abs(s));
  2869.                             this.hansenObjects[s + maxDegree][j] = new FieldHansenTesseralLinear<>(maxDegree, s, j, n0, maxHansen, field);
  2870.                         }
  2871.                     }
  2872.                     break;

  2874.                 default:
  2875.                     throw new OrekitInternalError(null);
  2876.             }

  2878.         }

  2880.         /** Compute init values for hansen objects.
  2881.          * @param context container for attributes
  2882.          * @param rows number of rows of the hansen matrix
  2883.          * @param columns columns number of columns of the hansen matrix
  2884.          */
  2885.         public void computeHansenObjectsInitValues(final FieldDSSTTesseralContext<T> context,
  2886.                                                    final int rows, final int columns) {
  2887.             hansenObjects[rows][columns].computeInitValues(context.getE2(), context.getChi(), context.getChi2());
  2888.         }

  2890.         /** Get hansen object.
  2891.          * @return hansenObjects
  2892.          */
  2893.         public FieldHansenTesseralLinear<T>[][] getHansenObjects() {
  2894.             return hansenObjects;
  2895.         }

  2897.     }

  2899. }
Created with JaCoCo 0.8.6.202009150832