DSSTTesseral.java

  1. /* Copyright 2002-2019 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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.DSFactory;
  33. import org.hipparchus.analysis.differentiation.DerivativeStructure;
  34. import org.hipparchus.analysis.differentiation.FDSFactory;
  35. import org.hipparchus.analysis.differentiation.FieldDerivativeStructure;
  36. import org.hipparchus.exception.LocalizedCoreFormats;
  37. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  38. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  39. import org.hipparchus.util.FastMath;
  40. import org.hipparchus.util.MathArrays;
  41. import org.hipparchus.util.MathUtils;
  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.     /** Factory for the DerivativeStructure instances. */
  194.     private DSFactory factory;

  196.     /** Factory for the DerivativeStructure instances. */
  197.     private Map<Field<?>, FDSFactory<?>> fieldFactory;

  199.     /** Flag for force model initialization with field elements. */
  200.     private boolean pendingInitialization;

  202.     /** Simple constructor.
  203.      * @param centralBodyFrame rotating body frame
  204.      * @param centralBodyRotationRate central body rotation rate (rad/s)
  205.      * @param provider provider for spherical harmonics
  206.      * @param maxDegreeTesseralSP maximal degree to consider for short periodics tesseral harmonics potential
  207.      *  (must be between 2 and {@code provider.getMaxDegree()})
  208.      * @param maxOrderTesseralSP maximal order to consider for short periodics tesseral harmonics potential
  209.      *  (must be between 0 and {@code provider.getMaxOrder()})
  210.      * @param maxEccPowTesseralSP maximum power of the eccentricity to use in summation over s for
  211.      * short periodic tesseral harmonics (without m-daily), should typically not exceed 4 as higher
  212.      * values will exceed computer capacity
  213.      * @param maxFrequencyShortPeriodics maximum frequency in mean longitude for short periodic computations
  214.      * (typically {@code maxDegreeTesseralSP} + {@code maxEccPowTesseralSP and no more than 12})
  215.      * @param maxDegreeMdailyTesseralSP maximal degree to consider for short periodics m-daily tesseral harmonics potential
  216.      *  (must be between 2 and {@code provider.getMaxDegree()})
  217.      * @param maxOrderMdailyTesseralSP maximal order to consider for short periodics m-daily tesseral harmonics potential
  218.      *  (must be between 0 and {@code provider.getMaxOrder()})
  219.      * @param maxEccPowMdailyTesseralSP maximum power of the eccentricity to use in summation over s for
  220.      * m-daily tesseral harmonics, (must be between 0 and {@code maxDegreeMdailyTesseralSP - 2},
  221.      * but should typically not exceed 4 as higher values will exceed computer capacity)
  222.      * @since 7.2
  223.      */
  224.     public DSSTTesseral(final Frame centralBodyFrame,
  225.                         final double centralBodyRotationRate,
  226.                         final UnnormalizedSphericalHarmonicsProvider provider,
  227.                         final int maxDegreeTesseralSP, final int maxOrderTesseralSP,
  228.                         final int maxEccPowTesseralSP, final int maxFrequencyShortPeriodics,
  229.                         final int maxDegreeMdailyTesseralSP, final int maxOrderMdailyTesseralSP,
  230.                         final int maxEccPowMdailyTesseralSP) {

  232.         gmParameterDriver = new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
  233.                                                 provider.getMu(), MU_SCALE,
  234.                                                 0.0, Double.POSITIVE_INFINITY);

  236.         // Central body rotating frame
  237.         this.bodyFrame = centralBodyFrame;

  239.         //Save the rotation rate
  240.         this.centralBodyRotationRate = centralBodyRotationRate;

  242.         // Central body rotation period in seconds
  243.         this.bodyPeriod = MathUtils.TWO_PI / centralBodyRotationRate;

  245.         // Provider for spherical harmonics
  246.         this.provider      = provider;
  247.         this.maxDegree     = provider.getMaxDegree();
  248.         this.maxOrder      = provider.getMaxOrder();

  250.         //set the maximum degree order for short periodics
  251.         checkIndexRange(maxDegreeTesseralSP, 2, maxDegree);
  252.         this.maxDegreeTesseralSP       = maxDegreeTesseralSP;

  254.         checkIndexRange(maxDegreeMdailyTesseralSP, 2, maxDegree);
  255.         this.maxDegreeMdailyTesseralSP = maxDegreeMdailyTesseralSP;

  257.         checkIndexRange(maxOrderTesseralSP, 0, maxOrder);
  258.         this.maxOrderTesseralSP        = maxOrderTesseralSP;

  260.         checkIndexRange(maxOrderMdailyTesseralSP, 0, maxOrder);
  261.         this.maxOrderMdailyTesseralSP  = maxOrderMdailyTesseralSP;

  263.         // set the maximum value for eccentricity power
  264.         this.maxEccPowTesseralSP       = maxEccPowTesseralSP;

  266.         checkIndexRange(maxEccPowMdailyTesseralSP, 0, maxDegreeMdailyTesseralSP - 2);
  267.         this.maxEccPowMdailyTesseralSP = maxEccPowMdailyTesseralSP;

  269.         // set the maximum value for frequency
  270.         this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;

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

  275.         pendingInitialization = true;

  277.         fieldShortPeriodTerms = new HashMap<>();
  278.         fieldHansen           = new HashMap<>();
  279.         fieldFactory          = new HashMap<>();

  281.     }

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

  294.     /** {@inheritDoc} */
  295.     @Override
  296.     public List<ShortPeriodTerms> initialize(final AuxiliaryElements auxiliaryElements,
  297.                                              final PropagationType type,
  298.                                              final double[] parameters) {

  300.         // Initializes specific parameters.
  301.         final DSSTTesseralContext context = initializeStep(auxiliaryElements, parameters);

  303.         // The following terms are only used for hansen objects initialization
  304.         final double ratio            = context.getRatio();
  305.         final int maxEccPow           = context.getMaxEccPow();
  306.         final List<Integer> resOrders = context.getResOrders();

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

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

  313.         factory = new DSFactory(1, 1);

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

  317.         shortPeriodTerms = new TesseralShortPeriodicCoefficients(bodyFrame, maxOrderMdailyTesseralSP,
  318.                                                                  maxDegreeTesseralSP < 0, nonResOrders,
  319.                                                                  mMax, maxFrequencyShortPeriodics, INTERPOLATION_POINTS,
  320.                                                                  new TimeSpanMap<Slot>(new Slot(mMax, maxFrequencyShortPeriodics,
  321.                                                                                                 INTERPOLATION_POINTS)));

  323.         final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
  324.         list.add(shortPeriodTerms);
  325.         return list;

  327.     }

  329.     /** {@inheritDoc} */
  330.     @Override
  331.     public <T extends RealFieldElement<T>> List<FieldShortPeriodTerms<T>> initialize(final FieldAuxiliaryElements<T> auxiliaryElements,
  332.                                                                                      final PropagationType type,
  333.                                                                                      final T[] parameters) {

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

  338.         if (pendingInitialization == true) {

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

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

  346.             // The following terms are only used for hansen objects initialization
  347.             final T      ratio            = context.getRatio();
  348.             final int maxEccPow           = context.getMaxEccPow();
  349.             final List<Integer> resOrders = context.getResOrders();


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

  354.             fieldFactory.put(field, new FDSFactory<>(field, 1, 1));
  355.             fieldHansen.put(field, new FieldHansenObjects<>(ratio, maxEccPow, resOrders, type, field));

  357.             pendingInitialization = false;
  358.         }

  360.         final FieldTesseralShortPeriodicCoefficients<T> ftspc =
  361.                         new FieldTesseralShortPeriodicCoefficients<>(bodyFrame, maxOrderMdailyTesseralSP,
  362.                                                                      maxDegreeTesseralSP < 0, nonResOrders,
  363.                                                                      mMax, maxFrequencyShortPeriodics, INTERPOLATION_POINTS,
  364.                                                                      new FieldTimeSpanMap<>(new FieldSlot<>(mMax,
  365.                                                                                                             maxFrequencyShortPeriodics,
  366.                                                                                                             INTERPOLATION_POINTS),
  367.                                                                                             field));

  369.         fieldShortPeriodTerms.put(field, ftspc);
  370.         return Collections.singletonList(ftspc);

  372.     }

  374.     /** Performs initialization at each integration step for the current force model.
  375.      *  <p>
  376.      *  This method aims at being called before mean elements rates computation.
  377.      *  </p>
  378.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  379.      *  @param parameters values of the force model parameters
  380.      *  @return new force model context
  381.      */
  382.     private DSSTTesseralContext initializeStep(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
  383.         return new DSSTTesseralContext(auxiliaryElements, bodyFrame, provider, maxFrequencyShortPeriodics, bodyPeriod, parameters);
  384.     }

  386.     /** Performs initialization at each integration step for the current force model.
  387.      *  <p>
  388.      *  This method aims at being called before mean elements rates computation.
  389.      *  </p>
  390.      *  @param <T> type of the elements
  391.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  392.      *  @param parameters values of the force model parameters
  393.      *  @return new force model context
  394.      */
  395.     private <T extends RealFieldElement<T>> FieldDSSTTesseralContext<T> initializeStep(final FieldAuxiliaryElements<T> auxiliaryElements,
  396.                                                                                        final T[] parameters) {
  397.         return new FieldDSSTTesseralContext<>(auxiliaryElements, bodyFrame, provider, maxFrequencyShortPeriodics, bodyPeriod, parameters);
  398.     }

  400.     /** {@inheritDoc} */
  401.     @Override
  402.     public double[] getMeanElementRate(final SpacecraftState spacecraftState,
  403.                                        final AuxiliaryElements auxiliaryElements, final double[] parameters) {

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

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

  411.         // Compute the cross derivative operator :
  412.         final double UAlphaGamma   = auxiliaryElements.getAlpha() * udu.getdUdGa() - auxiliaryElements.getGamma() * udu.getdUdAl();
  413.         final double UAlphaBeta    = auxiliaryElements.getAlpha() * udu.getdUdBe() - auxiliaryElements.getBeta()  * udu.getdUdAl();
  414.         final double UBetaGamma    = auxiliaryElements.getBeta() * udu.getdUdGa() - auxiliaryElements.getGamma() * udu.getdUdBe();
  415.         final double Uhk           = auxiliaryElements.getH() * udu.getdUdk()  - auxiliaryElements.getK() * udu.getdUdh();
  416.         final double pUagmIqUbgoAB = (auxiliaryElements.getP() * UAlphaGamma - I * auxiliaryElements.getQ() * UBetaGamma) * context.getOoAB();
  417.         final double UhkmUabmdUdl  = Uhk - UAlphaBeta - udu.getdUdl();

  419.         final double da =  context.getAx2oA() * udu.getdUdl();
  420.         final double dh =  context.getBoA() * udu.getdUdk() + auxiliaryElements.getK() * pUagmIqUbgoAB - auxiliaryElements.getH() * context.getBoABpo() * udu.getdUdl();
  421.         final double dk =  -(context.getBoA() * udu.getdUdh() + auxiliaryElements.getH() * pUagmIqUbgoAB + auxiliaryElements.getK() * context.getBoABpo() * udu.getdUdl());
  422.         final double dp =  context.getCo2AB() * (auxiliaryElements.getP() * UhkmUabmdUdl - UBetaGamma);
  423.         final double dq =  context.getCo2AB() * (auxiliaryElements.getQ() * UhkmUabmdUdl - I * UAlphaGamma);
  424.         final double dM = -context.getAx2oA() * udu.getdUda() + context.getBoABpo() * (auxiliaryElements.getH() * udu.getdUdh() + auxiliaryElements.getK() * udu.getdUdk()) + pUagmIqUbgoAB;

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

  429.     /** {@inheritDoc} */
  430.     @Override
  431.     public <T extends RealFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> spacecraftState,
  432.                                                                   final FieldAuxiliaryElements<T> auxiliaryElements,
  433.                                                                   final T[] parameters) {

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

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

  441.         @SuppressWarnings("unchecked")
  442.         final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);
  443.         // Access to potential U derivatives
  444.         final FieldUAnddU<T> udu = new FieldUAnddU<>(spacecraftState.getDate(), context, fho);

  446.         // Compute the cross derivative operator :
  447.         final T UAlphaGamma   = udu.getdUdGa().multiply(auxiliaryElements.getAlpha()).subtract(udu.getdUdAl().multiply(auxiliaryElements.getGamma()));
  448.         final T UAlphaBeta    = udu.getdUdBe().multiply(auxiliaryElements.getAlpha()).subtract(udu.getdUdAl().multiply(auxiliaryElements.getBeta()));
  449.         final T UBetaGamma    = udu.getdUdGa().multiply(auxiliaryElements.getBeta()).subtract(udu.getdUdBe().multiply(auxiliaryElements.getGamma()));
  450.         final T Uhk           = udu.getdUdk().multiply(auxiliaryElements.getH()).subtract(udu.getdUdh().multiply(auxiliaryElements.getK()));
  451.         final T pUagmIqUbgoAB = (UAlphaGamma.multiply(auxiliaryElements.getP()).subtract(UBetaGamma.multiply(auxiliaryElements.getQ()).multiply(I))).multiply(context.getOoAB());
  452.         final T UhkmUabmdUdl  = Uhk.subtract(UAlphaBeta).subtract(udu.getdUdl());

  454.         final T da = udu.getdUdl().multiply(context.getAx2oA());
  455.         final T dh = udu.getdUdk().multiply(context.getBoA()).add(pUagmIqUbgoAB.multiply(auxiliaryElements.getK())).subtract(udu.getdUdl().multiply(auxiliaryElements.getH()).multiply(context.getBoABpo()));
  456.         final T dk = (udu.getdUdh().multiply(context.getBoA()).add(pUagmIqUbgoAB.multiply(auxiliaryElements.getH())).add(udu.getdUdl().multiply(context.getBoABpo()).multiply(auxiliaryElements.getK()))).negate();
  457.         final T dp = context.getCo2AB().multiply(auxiliaryElements.getP().multiply(UhkmUabmdUdl).subtract(UBetaGamma));
  458.         final T dq = context.getCo2AB().multiply(auxiliaryElements.getQ().multiply(UhkmUabmdUdl).subtract(UAlphaGamma.multiply(I)));
  459.         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()));

  461.         final T[] elements = MathArrays.buildArray(field, 6);
  462.         elements[0] = da;
  463.         elements[1] = dk;
  464.         elements[2] = dh;
  465.         elements[3] = dq;
  466.         elements[4] = dp;
  467.         elements[5] = dM;

  469.         return elements;

  471.     }

  473.     /** {@inheritDoc} */
  474.     @Override
  475.     public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {

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

  479.         for (final SpacecraftState meanState : meanStates) {

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

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

  485.             // Initialise the Hansen coefficients
  486.             for (int s = -maxDegree; s <= maxDegree; s++) {
  487.                 // coefficients with j == 0 are always needed
  488.                 hansen.computeHansenObjectsInitValues(context, s + maxDegree, 0);
  489.                 if (maxDegreeTesseralSP >= 0) {
  490.                     // initialize other objects only if required
  491.                     for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  492.                         hansen.computeHansenObjectsInitValues(context, s + maxDegree, j);
  493.                     }
  494.                 }
  495.             }

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

  499.             // Compute coefficients
  500.             // Compute only if there is at least one non-resonant tesseral
  501.             if (!nonResOrders.isEmpty() || maxDegreeTesseralSP < 0) {
  502.                 // Generate the fourrier coefficients
  503.                 cjsjFourier.generateCoefficients(meanState.getDate(), context, hansen);

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

  508.                 // build the mDaily coefficients
  509.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  510.                     // build the coefficients
  511.                     buildCoefficients(cjsjFourier, meanState.getDate(), slot, m, 0, tnota, context);
  512.                 }

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

  518.                         for (int j : entry.getValue()) {
  519.                             // build the coefficients
  520.                             buildCoefficients(cjsjFourier, meanState.getDate(), slot, entry.getKey(), j, tnota, context);
  521.                         }
  522.                     }
  523.                 }
  524.             }

  526.         }

  528.     }

  530.     /** {@inheritDoc} */
  531.     @Override
  532.     @SuppressWarnings("unchecked")
  533.     public <T extends RealFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
  534.                                                                        final FieldSpacecraftState<T>... meanStates) {

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

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

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

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

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

  548.             final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);
  549.             // Initialise the Hansen coefficients
  550.             for (int s = -maxDegree; s <= maxDegree; s++) {
  551.                 // coefficients with j == 0 are always needed
  552.                 fho.computeHansenObjectsInitValues(context, s + maxDegree, 0);
  553.                 if (maxDegreeTesseralSP >= 0) {
  554.                     // initialize other objects only if required
  555.                     for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  556.                         fho.computeHansenObjectsInitValues(context, s + maxDegree, j);
  557.                     }
  558.                 }
  559.             }

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

  563.             // Compute coefficients
  564.             // Compute only if there is at least one non-resonant tesseral
  565.             if (!nonResOrders.isEmpty() || maxDegreeTesseralSP < 0) {
  566.                 // Generate the fourrier coefficients
  567.                 cjsjFourier.generateCoefficients(meanState.getDate(), context, fho, field);

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

  572.                 // build the mDaily coefficients
  573.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  574.                     // build the coefficients
  575.                     buildCoefficients(cjsjFourier, meanState.getDate(), slot, m, 0, tnota, context, field);
  576.                 }

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

  582.                         for (int j : entry.getValue()) {
  583.                             // build the coefficients
  584.                             buildCoefficients(cjsjFourier, meanState.getDate(), slot, entry.getKey(), j, tnota, context, field);
  585.                         }
  586.                     }
  587.                 }
  588.             }

  590.         }

  592.     }

  594.     /** {@inheritDoc} */
  595.     public ParameterDriver[] getParametersDrivers() {
  596.         return new ParameterDriver[] {
  597.             gmParameterDriver
  598.         };
  599.     }

  601.     /** Build a set of coefficients.
  602.      * @param cjsjFourier the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  603.      * @param date the current date
  604.      * @param slot slot to which the coefficients belong
  605.      * @param m m index
  606.      * @param j j index
  607.      * @param tnota 3n/2a
  608.      * @param context container for attributes
  609.      */
  610.     private void buildCoefficients(final FourierCjSjCoefficients cjsjFourier,
  611.                                    final AbsoluteDate date, final Slot slot,
  612.                                    final int m, final int j, final double tnota, final DSSTTesseralContext context) {

  614.         // Create local arrays
  615.         final double[] currentCijm = new double[] {0., 0., 0., 0., 0., 0.};
  616.         final double[] currentSijm = new double[] {0., 0., 0., 0., 0., 0.};

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

  621.         // initialise the coeficients
  622.         for (int i = 0; i < 6; i++) {
  623.             currentCijm[i] = -cjsjFourier.getSijm(i, j, m);
  624.             currentSijm[i] = cjsjFourier.getCijm(i, j, m);
  625.         }
  626.         // Add the separate part for δ<sub>6i</sub>
  627.         currentCijm[5] += tnota * oojnmt * cjsjFourier.getCijm(0, j, m);
  628.         currentSijm[5] += tnota * oojnmt * cjsjFourier.getSijm(0, j, m);

  630.         //Multiply by 1 / (jn - mθ<sup>.</sup>)
  631.         for (int i = 0; i < 6; i++) {
  632.             currentCijm[i] *= oojnmt;
  633.             currentSijm[i] *= oojnmt;
  634.         }

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

  640.     }

  642.      /** Build a set of coefficients.
  643.      * @param <T> the type of the field elements
  644.      * @param cjsjFourier the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  645.      * @param date the current date
  646.      * @param slot slot to which the coefficients belong
  647.      * @param m m index
  648.      * @param j j index
  649.      * @param tnota 3n/2a
  650.      * @param context container for attributes
  651.      * @param field field used by default
  652.      */
  653.     private <T extends RealFieldElement<T>> void buildCoefficients(final FieldFourierCjSjCoefficients<T> cjsjFourier,
  654.                                                                    final FieldAbsoluteDate<T> date,
  655.                                                                    final FieldSlot<T> slot,
  656.                                                                    final int m, final int j, final T tnota,
  657.                                                                    final FieldDSSTTesseralContext<T> context,
  658.                                                                    final Field<T> field) {

  660.         // Zero
  661.         final T zero = field.getZero();

  663.         // Create local arrays
  664.         final T[] currentCijm = MathArrays.buildArray(field, 6);
  665.         final T[] currentSijm = MathArrays.buildArray(field, 6);

  667.         Arrays.fill(currentCijm, zero);
  668.         Arrays.fill(currentSijm, zero);

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

  673.         // initialise the coeficients
  674.         for (int i = 0; i < 6; i++) {
  675.             currentCijm[i] = cjsjFourier.getSijm(i, j, m).negate();
  676.             currentSijm[i] = cjsjFourier.getCijm(i, j, m);
  677.         }
  678.         // Add the separate part for δ<sub>6i</sub>
  679.         currentCijm[5] = currentCijm[5].add(tnota.multiply(oojnmt).multiply(cjsjFourier.getCijm(0, j, m)));
  680.         currentSijm[5] = currentSijm[5].add(tnota.multiply(oojnmt).multiply(cjsjFourier.getSijm(0, j, m)));

  682.         //Multiply by 1 / (jn - mθ<sup>.</sup>)
  683.         for (int i = 0; i < 6; i++) {
  684.             currentCijm[i] = currentCijm[i].multiply(oojnmt);
  685.             currentSijm[i] = currentSijm[i].multiply(oojnmt);
  686.         }

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

  692.     }

  694.     /** {@inheritDoc} */
  695.     @Override
  696.     public EventDetector[] getEventsDetectors() {
  697.         return null;
  698.     }

  700.     /** {@inheritDoc} */
  701.     @Override
  702.     public <T extends RealFieldElement<T>> FieldEventDetector<T>[] getFieldEventsDetectors(final Field<T> field) {
  703.         return null;
  704.     }

  706.      /**
  707.       * Get the non-resonant tesseral terms in the central body spherical harmonic field.
  708.       *
  709.       * @param type type of the elements used during the propagation
  710.       * @param context container for attributes
  711.       */
  712.     private void getNonResonantTerms(final PropagationType type, final DSSTTesseralContext context) {

  714.         // Compute natural resonant terms
  715.         final double tolerance = 1. / FastMath.max(MIN_PERIOD_IN_SAT_REV,
  716.                                                    MIN_PERIOD_IN_SECONDS / context.getOrbitPeriod());

  718.         nonResOrders.clear();
  719.         for (int m = 1; m <= maxOrder; m++) {
  720.             final double resonance = context.getRatio() * m;
  721.             int jRes = 0;
  722.             final int jComputedRes = (int) FastMath.round(resonance);
  723.             if (jComputedRes > 0 && jComputedRes <= maxFrequencyShortPeriodics && FastMath.abs(resonance - jComputedRes) <= tolerance) {
  724.                 // Store each resonant index and order
  725.                 jRes = jComputedRes;
  726.             }

  728.             if (type == PropagationType.OSCULATING && maxDegreeTesseralSP >= 0 && m <= maxOrderTesseralSP) {
  729.                 //compute non resonant orders in the tesseral harmonic field
  730.                 final List<Integer> listJofM = new ArrayList<Integer>();
  731.                 //for the moment we take only the pairs (j,m) with |j| <= maxDegree + maxEccPow (from |s-j| <= maxEccPow and |s| <= maxDegree)
  732.                 for (int j = -maxFrequencyShortPeriodics; j <= maxFrequencyShortPeriodics; j++) {
  733.                     if (j != 0 && j != jRes) {
  734.                         listJofM.add(j);
  735.                     }
  736.                 }

  738.                 nonResOrders.put(m, listJofM);
  739.             }
  740.         }
  741.     }

  743.     /**
  744.      * Get the non-resonant tesseral terms in the central body spherical harmonic field.
  745.      *
  746.      * @param <T> type of the elements
  747.      * @param type type of the elements used during the propagation
  748.      * @param context container for attributes
  749.      * @param field field used by default
  750.      */
  751.     private <T extends RealFieldElement<T>> void getNonResonantTerms(final PropagationType type,
  752.                                                                      final FieldDSSTTesseralContext<T> context,
  753.                                                                      final Field<T> field) {

  755.         final T zero = field.getZero();
  756.         // Compute natural resonant terms
  757.         final T tolerance = FastMath.max(zero.add(MIN_PERIOD_IN_SAT_REV),
  758.                                          context.getOrbitPeriod().divide(MIN_PERIOD_IN_SECONDS).reciprocal()).reciprocal();

  760.         nonResOrders.clear();
  761.         for (int m = 1; m <= maxOrder; m++) {
  762.             final T resonance = context.getRatio().multiply(m);
  763.             int jRes = 0;
  764.             final int jComputedRes = (int) FastMath.round(resonance);
  765.             if (jComputedRes > 0 && jComputedRes <= maxFrequencyShortPeriodics && FastMath.abs(resonance.subtract(jComputedRes)).getReal() <= tolerance.getReal()) {
  766.                 // Store each resonant index and order
  767.                 jRes = jComputedRes;
  768.             }

  770.             if (type == PropagationType.OSCULATING && maxDegreeTesseralSP >= 0 && m <= maxOrderTesseralSP) {
  771.                 //compute non resonant orders in the tesseral harmonic field
  772.                 final List<Integer> listJofM = new ArrayList<Integer>();
  773.                 //for the moment we take only the pairs (j,m) with |j| <= maxDegree + maxEccPow (from |s-j| <= maxEccPow and |s| <= maxDegree)
  774.                 for (int j = -maxFrequencyShortPeriodics; j <= maxFrequencyShortPeriodics; j++) {
  775.                     if (j != 0 && j != jRes) {
  776.                         listJofM.add(j);
  777.                     }
  778.                 }

  780.                 nonResOrders.put(m, listJofM);
  781.             }
  782.         }
  783.     }

  785.     /** Compute the n-SUM for potential derivatives components.
  786.      *  @param date current date
  787.      *  @param j resonant index <i>j</i>
  788.      *  @param m resonant order <i>m</i>
  789.      *  @param s d'Alembert characteristic <i>s</i>
  790.      *  @param maxN maximum possible value for <i>n</i> index
  791.      *  @param roaPow powers of R/a up to degree <i>n</i>
  792.      *  @param ghMSJ G<sup>j</sup><sub>m,s</sub> and H<sup>j</sup><sub>m,s</sub> polynomials
  793.      *  @param gammaMNS &Gamma;<sup>m</sup><sub>n,s</sub>(γ) function
  794.      *  @param context container for attributes
  795.      *  @param hansenObjects initialization of hansen objects
  796.      *  @return Components of U<sub>n</sub> derivatives for fixed j, m, s
  797.      */
  798.     private double[][] computeNSum(final AbsoluteDate date,
  799.                                    final int j, final int m, final int s, final int maxN, final double[] roaPow,
  800.                                    final GHmsjPolynomials ghMSJ, final GammaMnsFunction gammaMNS, final DSSTTesseralContext context,
  801.                                    final HansenObjects hansenObjects) {

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

  806.         //spherical harmonics
  807.         final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);

  809.         // Potential derivatives components
  810.         double dUdaCos  = 0.;
  811.         double dUdaSin  = 0.;
  812.         double dUdhCos  = 0.;
  813.         double dUdhSin  = 0.;
  814.         double dUdkCos  = 0.;
  815.         double dUdkSin  = 0.;
  816.         double dUdlCos  = 0.;
  817.         double dUdlSin  = 0.;
  818.         double dUdAlCos = 0.;
  819.         double dUdAlSin = 0.;
  820.         double dUdBeCos = 0.;
  821.         double dUdBeSin = 0.;
  822.         double dUdGaCos = 0.;
  823.         double dUdGaSin = 0.;

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

  829.         // jacobi v, w, indices from 2.7.1-(15)
  830.         final int v = FastMath.abs(m - s);
  831.         final int w = FastMath.abs(m + s);

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

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

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

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

  849.                 // Inclination function Gamma and derivative
  850.                 final double gaMNS  = gammaMNS.getValue(m, n, s);
  851.                 final double dGaMNS = gammaMNS.getDerivative(m, n, s);

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

  857.                 // Gjms, Hjms polynomials and derivatives
  858.                 final double gMSJ   = ghMSJ.getGmsj(m, s, j);
  859.                 final double hMSJ   = ghMSJ.getHmsj(m, s, j);
  860.                 final double dGdh   = ghMSJ.getdGmsdh(m, s, j);
  861.                 final double dGdk   = ghMSJ.getdGmsdk(m, s, j);
  862.                 final double dGdA   = ghMSJ.getdGmsdAlpha(m, s, j);
  863.                 final double dGdB   = ghMSJ.getdGmsdBeta(m, s, j);
  864.                 final double dHdh   = ghMSJ.getdHmsdh(m, s, j);
  865.                 final double dHdk   = ghMSJ.getdHmsdk(m, s, j);
  866.                 final double dHdA   = ghMSJ.getdHmsdAlpha(m, s, j);
  867.                 final double dHdB   = ghMSJ.getdHmsdBeta(m, s, j);

  869.                 // Jacobi l-index from 2.7.1-(15)
  870.                 final int l = FastMath.min(n - m, n - FastMath.abs(s));
  871.                 // Jacobi polynomial and derivative
  872.                 final DerivativeStructure jacobi =
  873.                         JacobiPolynomials.getValue(l, v, w, factory.variable(0, auxiliaryElements.getGamma()));

  875.                 // Geopotential coefficients
  876.                 final double cnm = harmonics.getUnnormalizedCnm(n, m);
  877.                 final double snm = harmonics.getUnnormalizedSnm(n, m);

  879.                 // Common factors from expansion of equations 3.3-4
  880.                 final double cf_0      = roaPow[n] * Im * vMNS;
  881.                 final double cf_1      = cf_0 * gaMNS * jacobi.getValue();
  882.                 final double cf_2      = cf_1 * kJNS;
  883.                 final double gcPhs     = gMSJ * cnm + hMSJ * snm;
  884.                 final double gsMhc     = gMSJ * snm - hMSJ * cnm;
  885.                 final double dKgcPhsx2 = 2. * dkJNS * gcPhs;
  886.                 final double dKgsMhcx2 = 2. * dkJNS * gsMhc;
  887.                 final double dUdaCoef  = (n + 1) * cf_2;
  888.                 final double dUdlCoef  = j * cf_2;
  889.                 final double dUdGaCoef = cf_0 * kJNS * (jacobi.getValue() * dGaMNS + gaMNS * jacobi.getPartialDerivative(1));

  891.                 // dU / da components
  892.                 dUdaCos  += dUdaCoef * gcPhs;
  893.                 dUdaSin  += dUdaCoef * gsMhc;

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

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

  903.                 // dU / dLambda components
  904.                 dUdlCos  +=  dUdlCoef * gsMhc;
  905.                 dUdlSin  += -dUdlCoef * gcPhs;

  907.                 // dU / alpha components
  908.                 dUdAlCos += cf_2 * (dGdA * cnm + dHdA * snm);
  909.                 dUdAlSin += cf_2 * (dGdA * snm - dHdA * cnm);

  911.                 // dU / dBeta components
  912.                 dUdBeCos += cf_2 * (dGdB * cnm + dHdB * snm);
  913.                 dUdBeSin += cf_2 * (dGdB * snm - dHdB * cnm);

  915.                 // dU / dGamma components
  916.                 dUdGaCos += dUdGaCoef * gcPhs;
  917.                 dUdGaSin += dUdGaCoef * gsMhc;
  918.             }
  919.         }

  921.         return new double[][] {{dUdaCos,  dUdaSin},
  922.                                {dUdhCos,  dUdhSin},
  923.                                {dUdkCos,  dUdkSin},
  924.                                {dUdlCos,  dUdlSin},
  925.                                {dUdAlCos, dUdAlSin},
  926.                                {dUdBeCos, dUdBeSin},
  927.                                {dUdGaCos, dUdGaSin}};
  928.     }

  930.     /** Compute the n-SUM for potential derivatives components.
  931.      *  @param <T> the type of the field elements
  932.      *  @param date current date
  933.      *  @param j resonant index <i>j</i>
  934.      *  @param m resonant order <i>m</i>
  935.      *  @param s d'Alembert characteristic <i>s</i>
  936.      *  @param maxN maximum possible value for <i>n</i> index
  937.      *  @param roaPow powers of R/a up to degree <i>n</i>
  938.      *  @param ghMSJ G<sup>j</sup><sub>m,s</sub> and H<sup>j</sup><sub>m,s</sub> polynomials
  939.      *  @param gammaMNS &Gamma;<sup>m</sup><sub>n,s</sub>(γ) function
  940.      *  @param context container for attributes
  941.      *  @param hansenObjects initialization of hansen objects
  942.      *  @return Components of U<sub>n</sub> derivatives for fixed j, m, s
  943.      */
  944.     private <T extends RealFieldElement<T>> T[][] computeNSum(final FieldAbsoluteDate<T> date,
  945.                                                               final int j, final int m, final int s, final int maxN,
  946.                                                               final T[] roaPow,
  947.                                                               final FieldGHmsjPolynomials<T> ghMSJ,
  948.                                                               final FieldGammaMnsFunction<T> gammaMNS,
  949.                                                               final FieldDSSTTesseralContext<T> context,
  950.                                                               final FieldHansenObjects<T> hansenObjects) {

  952.         // Auxiliary elements related to the current orbit
  953.         final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
  954.         // Zero for initialization
  955.         final Field<T> field = date.getField();
  956.         final T zero = field.getZero();

  958.         //spherical harmonics
  959.         final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date.toAbsoluteDate());

  961.         // Potential derivatives components
  962.         T dUdaCos  = zero;
  963.         T dUdaSin  = zero;
  964.         T dUdhCos  = zero;
  965.         T dUdhSin  = zero;
  966.         T dUdkCos  = zero;
  967.         T dUdkSin  = zero;
  968.         T dUdlCos  = zero;
  969.         T dUdlSin  = zero;
  970.         T dUdAlCos = zero;
  971.         T dUdAlSin = zero;
  972.         T dUdBeCos = zero;
  973.         T dUdBeSin = zero;
  974.         T dUdGaCos = zero;
  975.         T dUdGaSin = zero;

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

  981.         // jacobi v, w, indices from 2.7.1-(15)
  982.         final int v = FastMath.abs(m - s);
  983.         final int w = FastMath.abs(m + s);

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

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

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

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

  1001.                 // Inclination function Gamma and derivative
  1002.                 final T gaMNS  = gammaMNS.getValue(m, n, s);
  1003.                 final T dGaMNS = gammaMNS.getDerivative(m, n, s);

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

  1009.                 // Gjms, Hjms polynomials and derivatives
  1010.                 final T gMSJ   = ghMSJ.getGmsj(m, s, j);
  1011.                 final T hMSJ   = ghMSJ.getHmsj(m, s, j);
  1012.                 final T dGdh   = ghMSJ.getdGmsdh(m, s, j);
  1013.                 final T dGdk   = ghMSJ.getdGmsdk(m, s, j);
  1014.                 final T dGdA   = ghMSJ.getdGmsdAlpha(m, s, j);
  1015.                 final T dGdB   = ghMSJ.getdGmsdBeta(m, s, j);
  1016.                 final T dHdh   = ghMSJ.getdHmsdh(m, s, j);
  1017.                 final T dHdk   = ghMSJ.getdHmsdk(m, s, j);
  1018.                 final T dHdA   = ghMSJ.getdHmsdAlpha(m, s, j);
  1019.                 final T dHdB   = ghMSJ.getdHmsdBeta(m, s, j);

  1021.                 // Jacobi l-index from 2.7.1-(15)
  1022.                 final int l = FastMath.min(n - m, n - FastMath.abs(s));
  1023.                 // Jacobi polynomial and derivative
  1024.                 @SuppressWarnings("unchecked")
  1025.                 final FDSFactory<T> fdsf = (FDSFactory<T>) fieldFactory.get(field);
  1026.                 final FieldDerivativeStructure<T> jacobi =
  1027.                         JacobiPolynomials.getValue(l, v, w, fdsf.variable(0, auxiliaryElements.getGamma()));

  1029.                 // Geopotential coefficients
  1030.                 final T cnm = zero.add(harmonics.getUnnormalizedCnm(n, m));
  1031.                 final T snm = zero.add(harmonics.getUnnormalizedSnm(n, m));

  1033.                 // Common factors from expansion of equations 3.3-4
  1034.                 final T cf_0      = roaPow[n].multiply(Im).multiply(vMNS);
  1035.                 final T cf_1      = cf_0.multiply(gaMNS).multiply(jacobi.getValue());
  1036.                 final T cf_2      = cf_1.multiply(kJNS);
  1037.                 final T gcPhs     = gMSJ.multiply(cnm).add(hMSJ.multiply(snm));
  1038.                 final T gsMhc     = gMSJ.multiply(snm).subtract(hMSJ.multiply(cnm));
  1039.                 final T dKgcPhsx2 = dkJNS.multiply(gcPhs).multiply(2.);
  1040.                 final T dKgsMhcx2 = dkJNS.multiply(gsMhc).multiply(2.);
  1041.                 final T dUdaCoef  = cf_2.multiply(n + 1);
  1042.                 final T dUdlCoef  = cf_2.multiply(j);
  1043.                 final T dUdGaCoef = cf_0.multiply(kJNS).multiply(dGaMNS.multiply(jacobi.getValue()).add(gaMNS.multiply(jacobi.getPartialDerivative(1))));

  1045.                 // dU / da components
  1046.                 dUdaCos  = dUdaCos.add(dUdaCoef.multiply(gcPhs));
  1047.                 dUdaSin  = dUdaSin.add(dUdaCoef.multiply(gsMhc));

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

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

  1057.                 // dU / dLambda components
  1058.                 dUdlCos  = dUdlCos.add(dUdlCoef.multiply(gsMhc));
  1059.                 dUdlSin  = dUdlSin.add(dUdlCoef.multiply(gcPhs).negate());

  1061.                 // dU / alpha components
  1062.                 dUdAlCos = dUdAlCos.add(cf_2.multiply(dGdA.multiply(cnm).add(dHdA.multiply(snm))));
  1063.                 dUdAlSin = dUdAlSin.add(cf_2.multiply(dGdA.multiply(snm).subtract(dHdA.multiply(cnm))));

  1065.                 // dU / dBeta components
  1066.                 dUdBeCos = dUdBeCos.add(cf_2.multiply(dGdB.multiply(cnm).add(dHdB.multiply(snm))));
  1067.                 dUdBeSin = dUdBeSin.add(cf_2.multiply(dGdB.multiply(snm).subtract(dHdB.multiply(cnm))));

  1069.                 // dU / dGamma components
  1070.                 dUdGaCos = dUdGaCos.add(dUdGaCoef.multiply(gcPhs));
  1071.                 dUdGaSin = dUdGaSin.add(dUdGaCoef.multiply(gsMhc));
  1072.             }
  1073.         }

  1075.         final T[][] derivatives = MathArrays.buildArray(field, 7, 2);
  1076.         derivatives[0][0] = dUdaCos;
  1077.         derivatives[0][1] = dUdaSin;
  1078.         derivatives[1][0] = dUdhCos;
  1079.         derivatives[1][1] = dUdhSin;
  1080.         derivatives[2][0] = dUdkCos;
  1081.         derivatives[2][1] = dUdkSin;
  1082.         derivatives[3][0] = dUdlCos;
  1083.         derivatives[3][1] = dUdlSin;
  1084.         derivatives[4][0] = dUdAlCos;
  1085.         derivatives[4][1] = dUdAlSin;
  1086.         derivatives[5][0] = dUdBeCos;
  1087.         derivatives[5][1] = dUdBeSin;
  1088.         derivatives[6][0] = dUdGaCos;
  1089.         derivatives[6][1] = dUdGaSin;

  1091.         return derivatives;

  1093.     }

  1095.     /** {@inheritDoc} */
  1096.     @Override
  1097.     public void registerAttitudeProvider(final AttitudeProvider attitudeProvider) {
  1098.         //nothing is done since this contribution is not sensitive to attitude
  1099.     }

  1101.     /** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
  1102.      *  <p>
  1103.      *  Those coefficients are given in Danielson paper by substituting the
  1104.      *  disturbing function (2.7.1-16) with m != 0 into (2.2-10)
  1105.      *  </p>
  1106.      */
  1107.     private class FourierCjSjCoefficients {

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

  1112.         /** The C<sub>i</sub><sup>jm</sup> coefficients.
  1113.          * <p>
  1114.          * The index order is [m][j][i] <br/>
  1115.          * The i index corresponds to the C<sub>i</sub><sup>jm</sup> coefficients used to
  1116.          * compute the following: <br/>
  1117.          * - da/dt <br/>
  1118.          * - dk/dt <br/>
  1119.          * - dh/dt / dk <br/>
  1120.          * - dq/dt <br/>
  1121.          * - dp/dt / dα <br/>
  1122.          * - dλ/dt / dβ <br/>
  1123.          * </p>
  1124.          */
  1125.         private final double[][][] cCoef;

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

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

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

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

  1151.         /** Create a set of C<sub>i</sub><sup>jm</sup> and S<sub>i</sub><sup>jm</sup> coefficients.
  1152.          *  @param jMax absolute limit for j ( -jMax <= j <= jMax )
  1153.          *  @param mMax maximum value for m
  1154.          */
  1155.         FourierCjSjCoefficients(final int jMax, final int mMax) {
  1156.             // initialise fields
  1157.             final int rows    = mMax + 1;
  1158.             final int columns = 2 * jMax + 1;
  1159.             this.jMax         = jMax;
  1160.             this.cCoef        = new double[rows][columns][6];
  1161.             this.sCoef        = new double[rows][columns][6];
  1162.             this.roaPow       = new double[maxDegree + 1];
  1163.             roaPow[0] = 1.;
  1164.         }

  1166.         /**
  1167.          * Generate the coefficients.
  1168.          * @param date the current date
  1169.          * @param context container for attributes
  1170.          * @param hansenObjects initialization of hansen objects
  1171.          */
  1172.         public void generateCoefficients(final AbsoluteDate date, final DSSTTesseralContext context,
  1173.                                          final HansenObjects hansenObjects) {

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

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

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

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

  1187.                 // R / a up to power degree
  1188.                 for (int i = 1; i <= maxRoaPower; i++) {
  1189.                     roaPow[i] = context.getRoa() * roaPow[i - 1];
  1190.                 }

  1192.                 //generate the m-daily coefficients
  1193.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1194.                     buildFourierCoefficients(date, m, 0, maxDegreeMdailyTesseralSP, context, hansenObjects);
  1195.                 }

  1197.                 // generate the other coefficients only if required
  1198.                 if (maxDegreeTesseralSP >= 0) {
  1199.                     for (int m: nonResOrders.keySet()) {
  1200.                         final List<Integer> listJ = nonResOrders.get(m);

  1202.                         for (int j: listJ) {
  1203.                             buildFourierCoefficients(date, m, j, maxDegreeTesseralSP, context, hansenObjects);
  1204.                         }
  1205.                     }
  1206.                 }
  1207.             }
  1208.         }

  1210.         /** Build a set of fourier coefficients for a given m and j.
  1211.          *
  1212.          * @param date the date of the coefficients
  1213.          * @param m m index
  1214.          * @param j j index
  1215.          * @param maxN  maximum value for n index
  1216.          * @param context container for attributes
  1217.          * @param hansenObjects initialization of hansen objects
  1218.          */
  1219.         private void buildFourierCoefficients(final AbsoluteDate date,
  1220.                final int m, final int j, final int maxN, final DSSTTesseralContext context,
  1221.                final HansenObjects hansenObjects) {

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

  1225.             // Potential derivatives components for a given non-resonant pair {j,m}
  1226.             double dRdaCos  = 0.;
  1227.             double dRdaSin  = 0.;
  1228.             double dRdhCos  = 0.;
  1229.             double dRdhSin  = 0.;
  1230.             double dRdkCos  = 0.;
  1231.             double dRdkSin  = 0.;
  1232.             double dRdlCos  = 0.;
  1233.             double dRdlSin  = 0.;
  1234.             double dRdAlCos = 0.;
  1235.             double dRdAlSin = 0.;
  1236.             double dRdBeCos = 0.;
  1237.             double dRdBeSin = 0.;
  1238.             double dRdGaCos = 0.;
  1239.             double dRdGaSin = 0.;

  1241.             // s-SUM from -sMin to sMax
  1242.             final int sMin = j == 0 ? maxEccPowMdailyTesseralSP : maxEccPowTesseralSP;
  1243.             final int sMax = j == 0 ? maxEccPowMdailyTesseralSP : maxEccPowTesseralSP;
  1244.             for (int s = 0; s <= sMax; s++) {

  1246.                 // n-SUM for s positive
  1247.                 final double[][] nSumSpos = computeNSum(date, j, m, s, maxN,
  1248.                                                         roaPow, ghMSJ, gammaMNS, context, hansenObjects);
  1249.                 dRdaCos  += nSumSpos[0][0];
  1250.                 dRdaSin  += nSumSpos[0][1];
  1251.                 dRdhCos  += nSumSpos[1][0];
  1252.                 dRdhSin  += nSumSpos[1][1];
  1253.                 dRdkCos  += nSumSpos[2][0];
  1254.                 dRdkSin  += nSumSpos[2][1];
  1255.                 dRdlCos  += nSumSpos[3][0];
  1256.                 dRdlSin  += nSumSpos[3][1];
  1257.                 dRdAlCos += nSumSpos[4][0];
  1258.                 dRdAlSin += nSumSpos[4][1];
  1259.                 dRdBeCos += nSumSpos[5][0];
  1260.                 dRdBeSin += nSumSpos[5][1];
  1261.                 dRdGaCos += nSumSpos[6][0];
  1262.                 dRdGaSin += nSumSpos[6][1];

  1264.                 // n-SUM for s negative
  1265.                 if (s > 0 && s <= sMin) {
  1266.                     final double[][] nSumSneg = computeNSum(date, j, m, -s, maxN,
  1267.                                                             roaPow, ghMSJ, gammaMNS, context, hansenObjects);
  1268.                     dRdaCos  += nSumSneg[0][0];
  1269.                     dRdaSin  += nSumSneg[0][1];
  1270.                     dRdhCos  += nSumSneg[1][0];
  1271.                     dRdhSin  += nSumSneg[1][1];
  1272.                     dRdkCos  += nSumSneg[2][0];
  1273.                     dRdkSin  += nSumSneg[2][1];
  1274.                     dRdlCos  += nSumSneg[3][0];
  1275.                     dRdlSin  += nSumSneg[3][1];
  1276.                     dRdAlCos += nSumSneg[4][0];
  1277.                     dRdAlSin += nSumSneg[4][1];
  1278.                     dRdBeCos += nSumSneg[5][0];
  1279.                     dRdBeSin += nSumSneg[5][1];
  1280.                     dRdGaCos += nSumSneg[6][0];
  1281.                     dRdGaSin += nSumSneg[6][1];
  1282.                 }
  1283.             }
  1284.             dRdaCos  *= -context.getMoa() / auxiliaryElements.getSma();
  1285.             dRdaSin  *= -context.getMoa() / auxiliaryElements.getSma();
  1286.             dRdhCos  *=  context.getMoa();
  1287.             dRdhSin  *=  context.getMoa();
  1288.             dRdkCos  *=  context.getMoa();
  1289.             dRdkSin  *=  context.getMoa();
  1290.             dRdlCos  *=  context.getMoa();
  1291.             dRdlSin  *=  context.getMoa();
  1292.             dRdAlCos *=  context.getMoa();
  1293.             dRdAlSin *=  context.getMoa();
  1294.             dRdBeCos *=  context.getMoa();
  1295.             dRdBeSin *=  context.getMoa();
  1296.             dRdGaCos *=  context.getMoa();
  1297.             dRdGaSin *=  context.getMoa();

  1299.             // Compute the cross derivative operator :
  1300.             final double RAlphaGammaCos   = auxiliaryElements.getAlpha() * dRdGaCos - auxiliaryElements.getGamma() * dRdAlCos;
  1301.             final double RAlphaGammaSin   = auxiliaryElements.getAlpha() * dRdGaSin - auxiliaryElements.getGamma() * dRdAlSin;
  1302.             final double RAlphaBetaCos    = auxiliaryElements.getAlpha() * dRdBeCos - auxiliaryElements.getBeta()  * dRdAlCos;
  1303.             final double RAlphaBetaSin    = auxiliaryElements.getAlpha() * dRdBeSin - auxiliaryElements.getBeta()  * dRdAlSin;
  1304.             final double RBetaGammaCos    =  auxiliaryElements.getBeta() * dRdGaCos - auxiliaryElements.getGamma() * dRdBeCos;
  1305.             final double RBetaGammaSin    =  auxiliaryElements.getBeta() * dRdGaSin - auxiliaryElements.getGamma() * dRdBeSin;
  1306.             final double RhkCos           =     auxiliaryElements.getH() * dRdkCos  -     auxiliaryElements.getK() * dRdhCos;
  1307.             final double RhkSin           =     auxiliaryElements.getH() * dRdkSin  -     auxiliaryElements.getK() * dRdhSin;
  1308.             final double pRagmIqRbgoABCos = (auxiliaryElements.getP() * RAlphaGammaCos - I * auxiliaryElements.getQ() * RBetaGammaCos) * context.getOoAB();
  1309.             final double pRagmIqRbgoABSin = (auxiliaryElements.getP() * RAlphaGammaSin - I * auxiliaryElements.getQ() * RBetaGammaSin) * context.getOoAB();
  1310.             final double RhkmRabmdRdlCos  =  RhkCos - RAlphaBetaCos - dRdlCos;
  1311.             final double RhkmRabmdRdlSin  =  RhkSin - RAlphaBetaSin - dRdlSin;

  1313.             // da/dt
  1314.             cCoef[m][j + jMax][0] = context.getAx2oA() * dRdlCos;
  1315.             sCoef[m][j + jMax][0] = context.getAx2oA() * dRdlSin;

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

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

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

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

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

  1338.         /** Get the coefficient C<sub>i</sub><sup>jm</sup>.
  1339.          * @param i i index - corresponds to the required variation
  1340.          * @param j j index
  1341.          * @param m m index
  1342.          * @return the coefficient C<sub>i</sub><sup>jm</sup>
  1343.          */
  1344.         public double getCijm(final int i, final int j, final int m) {
  1345.             return cCoef[m][j + jMax][i];
  1346.         }

  1348.         /** Get the coefficient S<sub>i</sub><sup>jm</sup>.
  1349.          * @param i i index - corresponds to the required variation
  1350.          * @param j j index
  1351.          * @param m m index
  1352.          * @return the coefficient S<sub>i</sub><sup>jm</sup>
  1353.          */
  1354.         public double getSijm(final int i, final int j, final int m) {
  1355.             return sCoef[m][j + jMax][i];
  1356.         }
  1357.     }

  1359.     /** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
  1360.      *  <p>
  1361.      *  Those coefficients are given in Danielson paper by substituting the
  1362.      *  disturbing function (2.7.1-16) with m != 0 into (2.2-10)
  1363.      *  </p>
  1364.      */
  1365.     private class FieldFourierCjSjCoefficients <T extends RealFieldElement<T>> {

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

  1370.         /** The C<sub>i</sub><sup>jm</sup> coefficients.
  1371.          * <p>
  1372.          * The index order is [m][j][i] <br/>
  1373.          * The i index corresponds to the C<sub>i</sub><sup>jm</sup> coefficients used to
  1374.          * compute the following: <br/>
  1375.          * - da/dt <br/>
  1376.          * - dk/dt <br/>
  1377.          * - dh/dt / dk <br/>
  1378.          * - dq/dt <br/>
  1379.          * - dp/dt / dα <br/>
  1380.          * - dλ/dt / dβ <br/>
  1381.          * </p>
  1382.          */
  1383.         private final T[][][] cCoef;

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

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

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

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

  1409.         /** Create a set of C<sub>i</sub><sup>jm</sup> and S<sub>i</sub><sup>jm</sup> coefficients.
  1410.          *  @param jMax absolute limit for j ( -jMax <= j <= jMax )
  1411.          *  @param mMax maximum value for m
  1412.          *  @param field field used by default
  1413.          */
  1414.         FieldFourierCjSjCoefficients(final int jMax, final int mMax, final Field<T> field) {
  1415.             // initialise fields
  1416.             final T zero = field.getZero();
  1417.             final int rows    = mMax + 1;
  1418.             final int columns = 2 * jMax + 1;
  1419.             this.jMax         = jMax;
  1420.             this.cCoef        = MathArrays.buildArray(field, rows, columns, 6);
  1421.             this.sCoef        = MathArrays.buildArray(field, rows, columns, 6);
  1422.             this.roaPow       = MathArrays.buildArray(field, maxDegree + 1);
  1423.             roaPow[0] = zero.add(1.);
  1424.         }

  1426.         /**
  1427.          * Generate the coefficients.
  1428.          * @param date the current date
  1429.          * @param context container for attributes
  1430.          * @param hansenObjects initialization of hansen objects
  1431.          * @param field field used by default
  1432.          */
  1433.         public void generateCoefficients(final FieldAbsoluteDate<T> date,
  1434.                                          final FieldDSSTTesseralContext<T> context,
  1435.                                          final FieldHansenObjects<T> hansenObjects,
  1436.                                          final Field<T> field) {

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

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

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

  1449.                 // R / a up to power degree
  1450.                 for (int i = 1; i <= maxRoaPower; i++) {
  1451.                     roaPow[i] = context.getRoa().multiply(roaPow[i - 1]);
  1452.                 }

  1454.                 //generate the m-daily coefficients
  1455.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1456.                     buildFourierCoefficients(date, m, 0, maxDegreeMdailyTesseralSP, context, hansenObjects, field);
  1457.                 }

  1459.                 // generate the other coefficients only if required
  1460.                 if (maxDegreeTesseralSP >= 0) {
  1461.                     for (int m: nonResOrders.keySet()) {
  1462.                         final List<Integer> listJ = nonResOrders.get(m);

  1464.                         for (int j: listJ) {
  1465.                             buildFourierCoefficients(date, m, j, maxDegreeTesseralSP, context, hansenObjects, field);
  1466.                         }
  1467.                     }
  1468.                 }
  1469.             }
  1470.         }

  1472.         /** Build a set of fourier coefficients for a given m and j.
  1473.          *
  1474.          * @param date the date of the coefficients
  1475.          * @param m m index
  1476.          * @param j j index
  1477.          * @param maxN  maximum value for n index
  1478.          * @param context container for attributes
  1479.          * @param hansenObjects initialization of hansen objects
  1480.          * @param field field used by default
  1481.          */
  1482.         private void buildFourierCoefficients(final FieldAbsoluteDate<T> date,
  1483.                                               final int m, final int j, final int maxN,
  1484.                                               final FieldDSSTTesseralContext<T> context,
  1485.                                               final FieldHansenObjects<T> hansenObjects,
  1486.                                               final Field<T> field) {

  1488.             // Zero
  1489.             final T zero = field.getZero();
  1490.             // Common parameters
  1491.             final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();

  1493.             // Potential derivatives components for a given non-resonant pair {j,m}
  1494.             T dRdaCos  = zero;
  1495.             T dRdaSin  = zero;
  1496.             T dRdhCos  = zero;
  1497.             T dRdhSin  = zero;
  1498.             T dRdkCos  = zero;
  1499.             T dRdkSin  = zero;
  1500.             T dRdlCos  = zero;
  1501.             T dRdlSin  = zero;
  1502.             T dRdAlCos = zero;
  1503.             T dRdAlSin = zero;
  1504.             T dRdBeCos = zero;
  1505.             T dRdBeSin = zero;
  1506.             T dRdGaCos = zero;
  1507.             T dRdGaSin = zero;

  1509.             // s-SUM from -sMin to sMax
  1510.             final int sMin = j == 0 ? maxEccPowMdailyTesseralSP : maxEccPowTesseralSP;
  1511.             final int sMax = j == 0 ? maxEccPowMdailyTesseralSP : maxEccPowTesseralSP;
  1512.             for (int s = 0; s <= sMax; s++) {

  1514.                 // n-SUM for s positive
  1515.                 final T[][] nSumSpos = computeNSum(date, j, m, s, maxN,
  1516.                                                         roaPow, ghMSJ, gammaMNS, context, hansenObjects);
  1517.                 dRdaCos  =  dRdaCos.add(nSumSpos[0][0]);
  1518.                 dRdaSin  =  dRdaSin.add(nSumSpos[0][1]);
  1519.                 dRdhCos  =  dRdhCos.add(nSumSpos[1][0]);
  1520.                 dRdhSin  =  dRdhSin.add(nSumSpos[1][1]);
  1521.                 dRdkCos  =  dRdkCos.add(nSumSpos[2][0]);
  1522.                 dRdkSin  =  dRdkSin.add(nSumSpos[2][1]);
  1523.                 dRdlCos  =  dRdlCos.add(nSumSpos[3][0]);
  1524.                 dRdlSin  =  dRdlSin.add(nSumSpos[3][1]);
  1525.                 dRdAlCos = dRdAlCos.add(nSumSpos[4][0]);
  1526.                 dRdAlSin = dRdAlSin.add(nSumSpos[4][1]);
  1527.                 dRdBeCos = dRdBeCos.add(nSumSpos[5][0]);
  1528.                 dRdBeSin = dRdBeSin.add(nSumSpos[5][1]);
  1529.                 dRdGaCos = dRdGaCos.add(nSumSpos[6][0]);
  1530.                 dRdGaSin = dRdGaSin.add(nSumSpos[6][1]);

  1532.                 // n-SUM for s negative
  1533.                 if (s > 0 && s <= sMin) {
  1534.                     final T[][] nSumSneg = computeNSum(date, j, m, -s, maxN,
  1535.                                                             roaPow, ghMSJ, gammaMNS, context, hansenObjects);
  1536.                     dRdaCos  =  dRdaCos.add(nSumSneg[0][0]);
  1537.                     dRdaSin  =  dRdaSin.add(nSumSneg[0][1]);
  1538.                     dRdhCos  =  dRdhCos.add(nSumSneg[1][0]);
  1539.                     dRdhSin  =  dRdhSin.add(nSumSneg[1][1]);
  1540.                     dRdkCos  =  dRdkCos.add(nSumSneg[2][0]);
  1541.                     dRdkSin  =  dRdkSin.add(nSumSneg[2][1]);
  1542.                     dRdlCos  =  dRdlCos.add(nSumSneg[3][0]);
  1543.                     dRdlSin  =  dRdlSin.add(nSumSneg[3][1]);
  1544.                     dRdAlCos = dRdAlCos.add(nSumSneg[4][0]);
  1545.                     dRdAlSin = dRdAlSin.add(nSumSneg[4][1]);
  1546.                     dRdBeCos = dRdBeCos.add(nSumSneg[5][0]);
  1547.                     dRdBeSin = dRdBeSin.add(nSumSneg[5][1]);
  1548.                     dRdGaCos = dRdGaCos.add(nSumSneg[6][0]);
  1549.                     dRdGaSin = dRdGaSin.add(nSumSneg[6][1]);
  1550.                 }
  1551.             }
  1552.             dRdaCos  =  dRdaCos.multiply(context.getMoa().negate().divide(auxiliaryElements.getSma()));
  1553.             dRdaSin  =  dRdaSin.multiply(context.getMoa().negate().divide(auxiliaryElements.getSma()));
  1554.             dRdhCos  =  dRdhCos.multiply(context.getMoa());
  1555.             dRdhSin  =  dRdhSin.multiply(context.getMoa());
  1556.             dRdkCos  =  dRdkCos.multiply(context.getMoa());
  1557.             dRdkSin  =  dRdkSin.multiply(context.getMoa());
  1558.             dRdlCos  =  dRdlCos.multiply(context.getMoa());
  1559.             dRdlSin  =  dRdlSin.multiply(context.getMoa());
  1560.             dRdAlCos = dRdAlCos.multiply(context.getMoa());
  1561.             dRdAlSin = dRdAlSin.multiply(context.getMoa());
  1562.             dRdBeCos = dRdBeCos.multiply(context.getMoa());
  1563.             dRdBeSin = dRdBeSin.multiply(context.getMoa());
  1564.             dRdGaCos = dRdGaCos.multiply(context.getMoa());
  1565.             dRdGaSin = dRdGaSin.multiply(context.getMoa());

  1567.             // Compute the cross derivative operator :
  1568.             final T RAlphaGammaCos   = auxiliaryElements.getAlpha().multiply(dRdGaCos).subtract(auxiliaryElements.getGamma().multiply(dRdAlCos));
  1569.             final T RAlphaGammaSin   = auxiliaryElements.getAlpha().multiply(dRdGaSin).subtract(auxiliaryElements.getGamma().multiply(dRdAlSin));
  1570.             final T RAlphaBetaCos    = auxiliaryElements.getAlpha().multiply(dRdBeCos).subtract(auxiliaryElements.getBeta().multiply(dRdAlCos));
  1571.             final T RAlphaBetaSin    = auxiliaryElements.getAlpha().multiply(dRdBeSin).subtract(auxiliaryElements.getBeta().multiply(dRdAlSin));
  1572.             final T RBetaGammaCos    =  auxiliaryElements.getBeta().multiply(dRdGaCos).subtract(auxiliaryElements.getGamma().multiply(dRdBeCos));
  1573.             final T RBetaGammaSin    =  auxiliaryElements.getBeta().multiply(dRdGaSin).subtract(auxiliaryElements.getGamma().multiply(dRdBeSin));
  1574.             final T RhkCos           =     auxiliaryElements.getH().multiply(dRdkCos).subtract(auxiliaryElements.getK().multiply(dRdhCos));
  1575.             final T RhkSin           =     auxiliaryElements.getH().multiply(dRdkSin).subtract(auxiliaryElements.getK().multiply(dRdhSin));
  1576.             final T pRagmIqRbgoABCos = (auxiliaryElements.getP().multiply(RAlphaGammaCos).subtract(auxiliaryElements.getQ().multiply(RBetaGammaCos).multiply(I))).multiply(context.getOoAB());
  1577.             final T pRagmIqRbgoABSin = (auxiliaryElements.getP().multiply(RAlphaGammaSin).subtract(auxiliaryElements.getQ().multiply(RBetaGammaSin).multiply(I))).multiply(context.getOoAB());
  1578.             final T RhkmRabmdRdlCos  =  RhkCos.subtract(RAlphaBetaCos).subtract(dRdlCos);
  1579.             final T RhkmRabmdRdlSin  =  RhkSin.subtract(RAlphaBetaSin).subtract(dRdlSin);

  1581.             // da/dt
  1582.             cCoef[m][j + jMax][0] = context.getAx2oA().multiply(dRdlCos);
  1583.             sCoef[m][j + jMax][0] = context.getAx2oA().multiply(dRdlSin);

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

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

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

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

  1601.             // dλ/dt
  1602.             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);
  1603.             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);
  1604.         }

  1606.         /** Get the coefficient C<sub>i</sub><sup>jm</sup>.
  1607.          * @param i i index - corresponds to the required variation
  1608.          * @param j j index
  1609.          * @param m m index
  1610.          * @return the coefficient C<sub>i</sub><sup>jm</sup>
  1611.          */
  1612.         public T getCijm(final int i, final int j, final int m) {
  1613.             return cCoef[m][j + jMax][i];
  1614.         }

  1616.         /** Get the coefficient S<sub>i</sub><sup>jm</sup>.
  1617.          * @param i i index - corresponds to the required variation
  1618.          * @param j j index
  1619.          * @param m m index
  1620.          * @return the coefficient S<sub>i</sub><sup>jm</sup>
  1621.          */
  1622.         public T getSijm(final int i, final int j, final int m) {
  1623.             return sCoef[m][j + jMax][i];
  1624.         }
  1625.     }

  1627.     /** 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
  1628.      * the short-periodic zonal contribution.
  1629.      *   <p>
  1630.      *  Those coefficients are given by expression 2.5.4-5 from the Danielson paper.
  1631.      *   </p>
  1632.      *
  1633.      * @author Sorin Scortan
  1634.      *
  1635.      */
  1636.     private static class TesseralShortPeriodicCoefficients implements ShortPeriodTerms {

  1638.         /** Retrograde factor I.
  1639.          *  <p>
  1640.          *  DSST model needs equinoctial orbit as internal representation.
  1641.          *  Classical equinoctial elements have discontinuities when inclination
  1642.          *  is close to zero. In this representation, I = +1. <br>
  1643.          *  To avoid this discontinuity, another representation exists and equinoctial
  1644.          *  elements can be expressed in a different way, called "retrograde" orbit.
  1645.          *  This implies I = -1. <br>
  1646.          *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  1647.          *  But for the sake of consistency with the theory, the retrograde factor
  1648.          *  has been kept in the formulas.
  1649.          *  </p>
  1650.          */
  1651.         private static final int I = 1;

  1653.         /** Central body rotating frame. */
  1654.         private final Frame bodyFrame;

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

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

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

  1665.         /** Maximum value for m index. */
  1666.         private final int mMax;

  1668.         /** Maximum value for j index. */
  1669.         private final int jMax;

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

  1674.         /** All coefficients slots. */
  1675.         private final transient TimeSpanMap<Slot> slots;

  1677.         /** Constructor.
  1678.          * @param bodyFrame central body rotating frame
  1679.          * @param maxOrderMdailyTesseralSP maximal order to consider for short periodics m-daily tesseral harmonics potential
  1680.          * @param mDailiesOnly flag to take into account only M-dailies harmonic tesserals for short periodic perturbations
  1681.          * @param nonResOrders lst of non resonant orders with j != 0
  1682.          * @param mMax maximum value for m index
  1683.          * @param jMax maximum value for j index
  1684.          * @param interpolationPoints number of points used in the interpolation process
  1685.          * @param slots all coefficients slots
  1686.          */
  1687.         TesseralShortPeriodicCoefficients(final Frame bodyFrame, final int maxOrderMdailyTesseralSP,
  1688.                                           final boolean mDailiesOnly, final SortedMap<Integer, List<Integer> > nonResOrders,
  1689.                                           final int mMax, final int jMax, final int interpolationPoints,
  1690.                                           final TimeSpanMap<Slot> slots) {
  1691.             this.bodyFrame                = bodyFrame;
  1692.             this.maxOrderMdailyTesseralSP = maxOrderMdailyTesseralSP;
  1693.             this.mDailiesOnly             = mDailiesOnly;
  1694.             this.nonResOrders             = nonResOrders;
  1695.             this.mMax                     = mMax;
  1696.             this.jMax                     = jMax;
  1697.             this.interpolationPoints      = interpolationPoints;
  1698.             this.slots                    = slots;
  1699.         }

  1701.         /** Get the slot valid for some date.
  1702.          * @param meanStates mean states defining the slot
  1703.          * @return slot valid at the specified date
  1704.          */
  1705.         public Slot createSlot(final SpacecraftState... meanStates) {
  1706.             final Slot         slot  = new Slot(mMax, jMax, interpolationPoints);
  1707.             final AbsoluteDate first = meanStates[0].getDate();
  1708.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  1709.             if (first.compareTo(last) <= 0) {
  1710.                 slots.addValidAfter(slot, first);
  1711.             } else {
  1712.                 slots.addValidBefore(slot, first);
  1713.             }
  1714.             return slot;
  1715.         }

  1717.         /** {@inheritDoc} */
  1718.         @Override
  1719.         public double[] value(final Orbit meanOrbit) {

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

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

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

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

  1734.                 // Central body rotation angle from equation 2.7.1-(3)(4).
  1735.                 final Transform t = bodyFrame.getTransformTo(auxiliaryElements.getFrame(), auxiliaryElements.getDate());
  1736.                 final Vector3D xB = t.transformVector(Vector3D.PLUS_I);
  1737.                 final Vector3D yB = t.transformVector(Vector3D.PLUS_J);
  1738.                 final Vector3D  f = auxiliaryElements.getVectorF();
  1739.                 final Vector3D  g = auxiliaryElements.getVectorG();
  1740.                 final double currentTheta = FastMath.atan2(-f.dotProduct(yB) + I * g.dotProduct(xB),
  1741.                                                             f.dotProduct(xB) + I * g.dotProduct(yB));

  1743.                 //Add the m-daily contribution
  1744.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1745.                     // Phase angle
  1746.                     final double jlMmt  = -m * currentTheta;
  1747.                     final double sinPhi = FastMath.sin(jlMmt);
  1748.                     final double cosPhi = FastMath.cos(jlMmt);

  1750.                     // compute contribution for each element
  1751.                     final double[] c = slot.getCijm(0, m, meanOrbit.getDate());
  1752.                     final double[] s = slot.getSijm(0, m, meanOrbit.getDate());
  1753.                     for (int i = 0; i < 6; i++) {
  1754.                         shortPeriodicVariation[i] += c[i] * cosPhi + s[i] * sinPhi;
  1755.                     }
  1756.                 }

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

  1763.                     for (int j : listJ) {
  1764.                         // Phase angle
  1765.                         final double jlMmt  = j * meanOrbit.getLM() - m * currentTheta;
  1766.                         final double sinPhi = FastMath.sin(jlMmt);
  1767.                         final double cosPhi = FastMath.cos(jlMmt);

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

  1776.                     }
  1777.                 }
  1778.             }

  1780.             return shortPeriodicVariation;

  1782.         }

  1784.         /** {@inheritDoc} */
  1785.         @Override
  1786.         public String getCoefficientsKeyPrefix() {
  1787.             return DSSTTesseral.SHORT_PERIOD_PREFIX;
  1788.         }

  1790.         /** {@inheritDoc}
  1791.          * <p>
  1792.          * For tesseral terms contributions,there are maxOrderMdailyTesseralSP
  1793.          * m-daily cMm coefficients, maxOrderMdailyTesseralSP m-daily sMm
  1794.          * coefficients, nbNonResonant cjm coefficients and nbNonResonant
  1795.          * sjm coefficients, where maxOrderMdailyTesseralSP and nbNonResonant both
  1796.          * depend on the orbit. The j index is the integer multiplier for the true
  1797.          * longitude argument and the m index is the integer multiplier for m-dailies.
  1798.          * </p>
  1799.          */
  1800.         @Override
  1801.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected) {

  1803.             // select the coefficients slot
  1804.             final Slot slot = slots.get(date);

  1806.             if (!nonResOrders.isEmpty() || mDailiesOnly) {
  1807.                 final Map<String, double[]> coefficients =
  1808.                                 new HashMap<String, double[]>(12 * maxOrderMdailyTesseralSP +
  1809.                                                               12 * nonResOrders.size());

  1811.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1812.                     storeIfSelected(coefficients, selected, slot.getCijm(0, m, date), DSSTTesseral.CM_COEFFICIENTS, m);
  1813.                     storeIfSelected(coefficients, selected, slot.getSijm(0, m, date), DSSTTesseral.SM_COEFFICIENTS, m);
  1814.                 }

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

  1820.                     for (int j : listJ) {
  1821.                         for (int i = 0; i < 6; ++i) {
  1822.                             storeIfSelected(coefficients, selected, slot.getCijm(j, m, date), "c", j, m);
  1823.                             storeIfSelected(coefficients, selected, slot.getSijm(j, m, date), "s", j, m);
  1824.                         }
  1825.                     }
  1826.                 }

  1828.                 return coefficients;

  1830.             } else {
  1831.                 return Collections.emptyMap();
  1832.             }

  1834.         }

  1836.         /** Put a coefficient in a map if selected.
  1837.          * @param map map to populate
  1838.          * @param selected set of coefficients that should be put in the map
  1839.          * (empty set means all coefficients are selected)
  1840.          * @param value coefficient value
  1841.          * @param id coefficient identifier
  1842.          * @param indices list of coefficient indices
  1843.          */
  1844.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
  1845.                                      final double[] value, final String id, final int... indices) {
  1846.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  1847.             keyBuilder.append(id);
  1848.             for (int index : indices) {
  1849.                 keyBuilder.append('[').append(index).append(']');
  1850.             }
  1851.             final String key = keyBuilder.toString();
  1852.             if (selected.isEmpty() || selected.contains(key)) {
  1853.                 map.put(key, value);
  1854.             }
  1855.         }

  1857.     }

  1859.     /** 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
  1860.      * the short-periodic zonal contribution.
  1861.      *   <p>
  1862.      *  Those coefficients are given by expression 2.5.4-5 from the Danielson paper.
  1863.      *   </p>
  1864.      *
  1865.      * @author Sorin Scortan
  1866.      *
  1867.      */
  1868.     private static class FieldTesseralShortPeriodicCoefficients <T extends RealFieldElement<T>> implements FieldShortPeriodTerms<T> {

  1870.         /** Retrograde factor I.
  1871.          *  <p>
  1872.          *  DSST model needs equinoctial orbit as internal representation.
  1873.          *  Classical equinoctial elements have discontinuities when inclination
  1874.          *  is close to zero. In this representation, I = +1. <br>
  1875.          *  To avoid this discontinuity, another representation exists and equinoctial
  1876.          *  elements can be expressed in a different way, called "retrograde" orbit.
  1877.          *  This implies I = -1. <br>
  1878.          *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  1879.          *  But for the sake of consistency with the theory, the retrograde factor
  1880.          *  has been kept in the formulas.
  1881.          *  </p>
  1882.          */
  1883.         private static final int I = 1;

  1885.         /** Central body rotating frame. */
  1886.         private final Frame bodyFrame;

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

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

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

  1897.         /** Maximum value for m index. */
  1898.         private final int mMax;

  1900.         /** Maximum value for j index. */
  1901.         private final int jMax;

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

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

  1909.         /** Constructor.
  1910.          * @param bodyFrame central body rotating frame
  1911.          * @param maxOrderMdailyTesseralSP maximal order to consider for short periodics m-daily tesseral harmonics potential
  1912.          * @param mDailiesOnly flag to take into account only M-dailies harmonic tesserals for short periodic perturbations
  1913.          * @param nonResOrders lst of non resonant orders with j != 0
  1914.          * @param mMax maximum value for m index
  1915.          * @param jMax maximum value for j index
  1916.          * @param interpolationPoints number of points used in the interpolation process
  1917.          * @param slots all coefficients slots
  1918.          */
  1919.         FieldTesseralShortPeriodicCoefficients(final Frame bodyFrame, final int maxOrderMdailyTesseralSP,
  1920.                                                final boolean mDailiesOnly, final SortedMap<Integer, List<Integer> > nonResOrders,
  1921.                                                final int mMax, final int jMax, final int interpolationPoints,
  1922.                                                final FieldTimeSpanMap<FieldSlot<T>, T> slots) {
  1923.             this.bodyFrame                = bodyFrame;
  1924.             this.maxOrderMdailyTesseralSP = maxOrderMdailyTesseralSP;
  1925.             this.mDailiesOnly             = mDailiesOnly;
  1926.             this.nonResOrders             = nonResOrders;
  1927.             this.mMax                     = mMax;
  1928.             this.jMax                     = jMax;
  1929.             this.interpolationPoints      = interpolationPoints;
  1930.             this.slots                    = slots;
  1931.         }

  1933.         /** Get the slot valid for some date.
  1934.          * @param meanStates mean states defining the slot
  1935.          * @return slot valid at the specified date
  1936.          */
  1937.         @SuppressWarnings("unchecked")
  1938.         public FieldSlot<T> createSlot(final FieldSpacecraftState<T>... meanStates) {
  1939.             final FieldSlot<T>         slot  = new FieldSlot<>(mMax, jMax, interpolationPoints);
  1940.             final FieldAbsoluteDate<T> first = meanStates[0].getDate();
  1941.             final FieldAbsoluteDate<T> last  = meanStates[meanStates.length - 1].getDate();
  1942.             if (first.compareTo(last) <= 0) {
  1943.                 slots.addValidAfter(slot, first);
  1944.             } else {
  1945.                 slots.addValidBefore(slot, first);
  1946.             }
  1947.             return slot;
  1948.         }

  1950.         /** {@inheritDoc} */
  1951.         @Override
  1952.         public T[] value(final FieldOrbit<T> meanOrbit) {

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

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

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

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

  1967.                 // Central body rotation angle from equation 2.7.1-(3)(4).
  1968.                 final FieldTransform<T> t = bodyFrame.getTransformTo(auxiliaryElements.getFrame(), auxiliaryElements.getDate());
  1969.                 final FieldVector3D<T> xB = t.transformVector(Vector3D.PLUS_I);
  1970.                 final FieldVector3D<T> yB = t.transformVector(Vector3D.PLUS_J);
  1971.                 final FieldVector3D<T>  f = auxiliaryElements.getVectorF();
  1972.                 final FieldVector3D<T>  g = auxiliaryElements.getVectorG();
  1973.                 final T currentTheta = FastMath.atan2(f.dotProduct(yB).negate().add(g.dotProduct(xB).multiply(I)),
  1974.                                                       f.dotProduct(xB).add(g.dotProduct(yB).multiply(I)));

  1976.                 //Add the m-daily contribution
  1977.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1978.                     // Phase angle
  1979.                     final T jlMmt  = currentTheta.multiply(-m);
  1980.                     final T sinPhi = FastMath.sin(jlMmt);
  1981.                     final T cosPhi = FastMath.cos(jlMmt);

  1983.                     // compute contribution for each element
  1984.                     final T[] c = slot.getCijm(0, m, meanOrbit.getDate());
  1985.                     final T[] s = slot.getSijm(0, m, meanOrbit.getDate());
  1986.                     for (int i = 0; i < 6; i++) {
  1987.                         shortPeriodicVariation[i] = shortPeriodicVariation[i].add(c[i].multiply(cosPhi).add(s[i].multiply(sinPhi)));
  1988.                     }
  1989.                 }

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

  1996.                     for (int j : listJ) {
  1997.                         // Phase angle
  1998.                         final T jlMmt  = meanOrbit.getLM().multiply(j).subtract(currentTheta.multiply(m));
  1999.                         final T sinPhi = FastMath.sin(jlMmt);
  2000.                         final T cosPhi = FastMath.cos(jlMmt);

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

  2009.                     }
  2010.                 }
  2011.             }

  2013.             return shortPeriodicVariation;

  2015.         }

  2017.         /** {@inheritDoc} */
  2018.         @Override
  2019.         public String getCoefficientsKeyPrefix() {
  2020.             return DSSTTesseral.SHORT_PERIOD_PREFIX;
  2021.         }

  2023.         /** {@inheritDoc}
  2024.          * <p>
  2025.          * For tesseral terms contributions,there are maxOrderMdailyTesseralSP
  2026.          * m-daily cMm coefficients, maxOrderMdailyTesseralSP m-daily sMm
  2027.          * coefficients, nbNonResonant cjm coefficients and nbNonResonant
  2028.          * sjm coefficients, where maxOrderMdailyTesseralSP and nbNonResonant both
  2029.          * depend on the orbit. The j index is the integer multiplier for the true
  2030.          * longitude argument and the m index is the integer multiplier for m-dailies.
  2031.          * </p>
  2032.          */
  2033.         @Override
  2034.         public Map<String, T[]> getCoefficients(final FieldAbsoluteDate<T> date, final Set<String> selected) {

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

  2039.             if (!nonResOrders.isEmpty() || mDailiesOnly) {
  2040.                 final Map<String, T[]> coefficients =
  2041.                                 new HashMap<String, T[]>(12 * maxOrderMdailyTesseralSP +
  2042.                                                          12 * nonResOrders.size());

  2044.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  2045.                     storeIfSelected(coefficients, selected, slot.getCijm(0, m, date), DSSTTesseral.CM_COEFFICIENTS, m);
  2046.                     storeIfSelected(coefficients, selected, slot.getSijm(0, m, date), DSSTTesseral.SM_COEFFICIENTS, m);
  2047.                 }

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

  2053.                     for (int j : listJ) {
  2054.                         for (int i = 0; i < 6; ++i) {
  2055.                             storeIfSelected(coefficients, selected, slot.getCijm(j, m, date), "c", j, m);
  2056.                             storeIfSelected(coefficients, selected, slot.getSijm(j, m, date), "s", j, m);
  2057.                         }
  2058.                     }
  2059.                 }

  2061.                 return coefficients;

  2063.             } else {
  2064.                 return Collections.emptyMap();
  2065.             }

  2067.         }

  2069.         /** Put a coefficient in a map if selected.
  2070.          * @param map map to populate
  2071.          * @param selected set of coefficients that should be put in the map
  2072.          * (empty set means all coefficients are selected)
  2073.          * @param value coefficient value
  2074.          * @param id coefficient identifier
  2075.          * @param indices list of coefficient indices
  2076.          */
  2077.         private void storeIfSelected(final Map<String, T[]> map, final Set<String> selected,
  2078.                                      final T[] value, final String id, final int... indices) {
  2079.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  2080.             keyBuilder.append(id);
  2081.             for (int index : indices) {
  2082.                 keyBuilder.append('[').append(index).append(']');
  2083.             }
  2084.             final String key = keyBuilder.toString();
  2085.             if (selected.isEmpty() || selected.contains(key)) {
  2086.                 map.put(key, value);
  2087.             }
  2088.         }
  2089.     }

  2091.     /** Coefficients valid for one time slot. */
  2092.     private static class Slot {

  2094.         /** The coefficients C<sub>i</sub><sup>j</sup><sup>m</sup>.
  2095.          * <p>
  2096.          * The index order is cijm[m][j][i] <br/>
  2097.          * i corresponds to the equinoctial element, as follows: <br/>
  2098.          * - i=0 for a <br/>
  2099.          * - i=1 for k <br/>
  2100.          * - i=2 for h <br/>
  2101.          * - i=3 for q <br/>
  2102.          * - i=4 for p <br/>
  2103.          * - i=5 for λ <br/>
  2104.          * </p>
  2105.          */
  2106.         private final ShortPeriodicsInterpolatedCoefficient[][] cijm;

  2108.         /** The coefficients S<sub>i</sub><sup>j</sup><sup>m</sup>.
  2109.          * <p>
  2110.          * The index order is sijm[m][j][i] <br/>
  2111.          * i corresponds to the equinoctial element, as follows: <br/>
  2112.          * - i=0 for a <br/>
  2113.          * - i=1 for k <br/>
  2114.          * - i=2 for h <br/>
  2115.          * - i=3 for q <br/>
  2116.          * - i=4 for p <br/>
  2117.          * - i=5 for λ <br/>
  2118.          * </p>
  2119.          */
  2120.         private final ShortPeriodicsInterpolatedCoefficient[][] sijm;

  2122.         /** Simple constructor.
  2123.          *  @param mMax maximum value for m index
  2124.          *  @param jMax maximum value for j index
  2125.          *  @param interpolationPoints number of points used in the interpolation process
  2126.          */
  2127.         Slot(final int mMax, final int jMax, final int interpolationPoints) {

  2129.             final int rows    = mMax + 1;
  2130.             final int columns = 2 * jMax + 1;
  2131.             cijm = new ShortPeriodicsInterpolatedCoefficient[rows][columns];
  2132.             sijm = new ShortPeriodicsInterpolatedCoefficient[rows][columns];
  2133.             for (int m = 1; m <= mMax; m++) {
  2134.                 for (int j = -jMax; j <= jMax; j++) {
  2135.                     cijm[m][j + jMax] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  2136.                     sijm[m][j + jMax] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  2137.                 }
  2138.             }

  2140.         }

  2142.         /** Get C<sub>i</sub><sup>j</sup><sup>m</sup>.
  2143.          *
  2144.          * @param j j index
  2145.          * @param m m index
  2146.          * @param date the date
  2147.          * @return C<sub>i</sub><sup>j</sup><sup>m</sup>
  2148.          */
  2149.         double[] getCijm(final int j, final int m, final AbsoluteDate date) {
  2150.             final int jMax = (cijm[m].length - 1) / 2;
  2151.             return cijm[m][j + jMax].value(date);
  2152.         }

  2154.         /** Get S<sub>i</sub><sup>j</sup><sup>m</sup>.
  2155.          *
  2156.          * @param j j index
  2157.          * @param m m index
  2158.          * @param date the date
  2159.          * @return S<sub>i</sub><sup>j</sup><sup>m</sup>
  2160.          */
  2161.         double[] getSijm(final int j, final int m, final AbsoluteDate date) {
  2162.             final int jMax = (cijm[m].length - 1) / 2;
  2163.             return sijm[m][j + jMax].value(date);
  2164.         }

  2166.     }

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

  2171.         /** The coefficients C<sub>i</sub><sup>j</sup><sup>m</sup>.
  2172.          * <p>
  2173.          * The index order is cijm[m][j][i] <br/>
  2174.          * i corresponds to the equinoctial element, as follows: <br/>
  2175.          * - i=0 for a <br/>
  2176.          * - i=1 for k <br/>
  2177.          * - i=2 for h <br/>
  2178.          * - i=3 for q <br/>
  2179.          * - i=4 for p <br/>
  2180.          * - i=5 for λ <br/>
  2181.          * </p>
  2182.          */
  2183.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[][] cijm;

  2185.         /** The coefficients S<sub>i</sub><sup>j</sup><sup>m</sup>.
  2186.          * <p>
  2187.          * The index order is sijm[m][j][i] <br/>
  2188.          * i corresponds to the equinoctial element, as follows: <br/>
  2189.          * - i=0 for a <br/>
  2190.          * - i=1 for k <br/>
  2191.          * - i=2 for h <br/>
  2192.          * - i=3 for q <br/>
  2193.          * - i=4 for p <br/>
  2194.          * - i=5 for λ <br/>
  2195.          * </p>
  2196.          */
  2197.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[][] sijm;

  2199.         /** Simple constructor.
  2200.          *  @param mMax maximum value for m index
  2201.          *  @param jMax maximum value for j index
  2202.          *  @param interpolationPoints number of points used in the interpolation process
  2203.          */
  2204.         @SuppressWarnings("unchecked")
  2205.         FieldSlot(final int mMax, final int jMax, final int interpolationPoints) {

  2207.             final int rows    = mMax + 1;
  2208.             final int columns = 2 * jMax + 1;
  2209.             cijm = (FieldShortPeriodicsInterpolatedCoefficient<T>[][]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows, columns);
  2210.             sijm = (FieldShortPeriodicsInterpolatedCoefficient<T>[][]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows, columns);
  2211.             for (int m = 1; m <= mMax; m++) {
  2212.                 for (int j = -jMax; j <= jMax; j++) {
  2213.                     cijm[m][j + jMax] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  2214.                     sijm[m][j + jMax] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  2215.                 }
  2216.             }

  2218.         }

  2220.         /** Get C<sub>i</sub><sup>j</sup><sup>m</sup>.
  2221.          *
  2222.          * @param j j index
  2223.          * @param m m index
  2224.          * @param date the date
  2225.          * @return C<sub>i</sub><sup>j</sup><sup>m</sup>
  2226.          */
  2227.         T[] getCijm(final int j, final int m, final FieldAbsoluteDate<T> date) {
  2228.             final int jMax = (cijm[m].length - 1) / 2;
  2229.             return cijm[m][j + jMax].value(date);
  2230.         }

  2232.         /** Get S<sub>i</sub><sup>j</sup><sup>m</sup>.
  2233.          *
  2234.          * @param j j index
  2235.          * @param m m index
  2236.          * @param date the date
  2237.          * @return S<sub>i</sub><sup>j</sup><sup>m</sup>
  2238.          */
  2239.         T[] getSijm(final int j, final int m, final FieldAbsoluteDate<T> date) {
  2240.             final int jMax = (cijm[m].length - 1) / 2;
  2241.             return sijm[m][j + jMax].value(date);
  2242.         }

  2244.     }

  2246.     /** Compute potential and potential derivatives with respect to orbital parameters.
  2247.      *  <p>The following elements are computed from expression 3.3 - (4).
  2248.      *  <pre>
  2249.      *  dU / da
  2250.      *  dU / dh
  2251.      *  dU / dk
  2252.      *  dU / dλ
  2253.      *  dU / dα
  2254.      *  dU / dβ
  2255.      *  dU / dγ
  2256.      *  </pre>
  2257.      *  </p>
  2258.      */
  2259.     private class UAnddU {

  2261.         /** dU / da. */
  2262.         private  double dUda;

  2264.         /** dU / dk. */
  2265.         private double dUdk;

  2267.         /** dU / dh. */
  2268.         private double dUdh;

  2270.         /** dU / dl. */
  2271.         private double dUdl;

  2273.         /** dU / dAlpha. */
  2274.         private double dUdAl;

  2276.         /** dU / dBeta. */
  2277.         private double dUdBe;

  2279.         /** dU / dGamma. */
  2280.         private double dUdGa;

  2282.         /** Simple constuctor.
  2283.          * @param date current date
  2284.          * @param context container for attributes
  2285.          * @param hansen hansen objects
  2286.          */
  2287.         UAnddU(final AbsoluteDate date, final DSSTTesseralContext context, final HansenObjects hansen) {

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

  2292.             // Potential derivatives
  2293.             dUda  = 0.;
  2294.             dUdh  = 0.;
  2295.             dUdk  = 0.;
  2296.             dUdl  = 0.;
  2297.             dUdAl = 0.;
  2298.             dUdBe = 0.;
  2299.             dUdGa = 0.;

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

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

  2309.                 // R / a up to power degree
  2310.                 final double[] roaPow = new double[maxDegree + 1];
  2311.                 roaPow[0] = 1.;
  2312.                 for (int i = 1; i <= maxDegree; i++) {
  2313.                     roaPow[i] = context.getRoa() * roaPow[i - 1];
  2314.                 }

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

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

  2322.                     // Phase angle
  2323.                     final double jlMmt  = j * auxiliaryElements.getLM() - m * context.getTheta();
  2324.                     final double sinPhi = FastMath.sin(jlMmt);
  2325.                     final double cosPhi = FastMath.cos(jlMmt);

  2327.                     // Potential derivatives components for a given resonant pair {j,m}
  2328.                     double dUdaCos  = 0.;
  2329.                     double dUdaSin  = 0.;
  2330.                     double dUdhCos  = 0.;
  2331.                     double dUdhSin  = 0.;
  2332.                     double dUdkCos  = 0.;
  2333.                     double dUdkSin  = 0.;
  2334.                     double dUdlCos  = 0.;
  2335.                     double dUdlSin  = 0.;
  2336.                     double dUdAlCos = 0.;
  2337.                     double dUdAlSin = 0.;
  2338.                     double dUdBeCos = 0.;
  2339.                     double dUdBeSin = 0.;
  2340.                     double dUdGaCos = 0.;
  2341.                     double dUdGaSin = 0.;

  2343.                     // s-SUM from -sMin to sMax
  2344.                     final int sMin = FastMath.min(context.getMaxEccPow() - j, maxDegree);
  2345.                     final int sMax = FastMath.min(context.getMaxEccPow() + j, maxDegree);
  2346.                     for (int s = 0; s <= sMax; s++) {

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

  2351.                         // n-SUM for s positive
  2352.                         final double[][] nSumSpos = computeNSum(date, j, m, s, maxDegree,
  2353.                                                                 roaPow, ghMSJ, gammaMNS, context, hansen);
  2354.                         dUdaCos  += nSumSpos[0][0];
  2355.                         dUdaSin  += nSumSpos[0][1];
  2356.                         dUdhCos  += nSumSpos[1][0];
  2357.                         dUdhSin  += nSumSpos[1][1];
  2358.                         dUdkCos  += nSumSpos[2][0];
  2359.                         dUdkSin  += nSumSpos[2][1];
  2360.                         dUdlCos  += nSumSpos[3][0];
  2361.                         dUdlSin  += nSumSpos[3][1];
  2362.                         dUdAlCos += nSumSpos[4][0];
  2363.                         dUdAlSin += nSumSpos[4][1];
  2364.                         dUdBeCos += nSumSpos[5][0];
  2365.                         dUdBeSin += nSumSpos[5][1];
  2366.                         dUdGaCos += nSumSpos[6][0];
  2367.                         dUdGaSin += nSumSpos[6][1];

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

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

  2393.                     // Assembly of potential derivatives componants
  2394.                     dUda  += cosPhi * dUdaCos  + sinPhi * dUdaSin;
  2395.                     dUdh  += cosPhi * dUdhCos  + sinPhi * dUdhSin;
  2396.                     dUdk  += cosPhi * dUdkCos  + sinPhi * dUdkSin;
  2397.                     dUdl  += cosPhi * dUdlCos  + sinPhi * dUdlSin;
  2398.                     dUdAl += cosPhi * dUdAlCos + sinPhi * dUdAlSin;
  2399.                     dUdBe += cosPhi * dUdBeCos + sinPhi * dUdBeSin;
  2400.                     dUdGa += cosPhi * dUdGaCos + sinPhi * dUdGaSin;
  2401.                 }

  2403.                 this.dUda  = dUda * (-context.getMoa() / auxiliaryElements.getSma());
  2404.                 this.dUdh  = dUdh * context.getMoa();
  2405.                 this.dUdk  = dUdk * context.getMoa();
  2406.                 this.dUdl  = dUdl * context.getMoa();
  2407.                 this.dUdAl = dUdAl * context.getMoa();
  2408.                 this.dUdBe = dUdBe * context.getMoa();
  2409.                 this.dUdGa = dUdGa * context.getMoa();
  2410.             }

  2412.         }

  2414.         /** Return value of dU / da.
  2415.          * @return dUda
  2416.          */
  2417.         public double getdUda() {
  2418.             return dUda;
  2419.         }

  2421.         /** Return value of dU / dk.
  2422.          * @return dUdk
  2423.          */
  2424.         public double getdUdk() {
  2425.             return dUdk;
  2426.         }

  2428.         /** Return value of dU / dh.
  2429.          * @return dUdh
  2430.          */
  2431.         public double getdUdh() {
  2432.             return dUdh;
  2433.         }

  2435.         /** Return value of dU / dl.
  2436.          * @return dUdl
  2437.          */
  2438.         public double getdUdl() {
  2439.             return dUdl;
  2440.         }

  2442.         /** Return value of dU / dAlpha.
  2443.          * @return dUdAl
  2444.          */
  2445.         public double getdUdAl() {
  2446.             return dUdAl;
  2447.         }

  2449.         /** Return value of dU / dBeta.
  2450.          * @return dUdBe
  2451.          */
  2452.         public double getdUdBe() {
  2453.             return dUdBe;
  2454.         }

  2456.         /** Return value of dU / dGamma.
  2457.          * @return dUdGa
  2458.          */
  2459.         public double getdUdGa() {
  2460.             return dUdGa;
  2461.         }

  2463.     }

  2465.     /**  Computes the potential U derivatives.
  2466.      *  <p>The following elements are computed from expression 3.3 - (4).
  2467.      *  <pre>
  2468.      *  dU / da
  2469.      *  dU / dh
  2470.      *  dU / dk
  2471.      *  dU / dλ
  2472.      *  dU / dα
  2473.      *  dU / dβ
  2474.      *  dU / dγ
  2475.      *  </pre>
  2476.      *  </p>
  2477.      */
  2478.     private class FieldUAnddU <T extends RealFieldElement<T>> {

  2480.         /** dU / da. */
  2481.         private T dUda;

  2483.         /** dU / dk. */
  2484.         private T dUdk;

  2486.         /** dU / dh. */
  2487.         private T dUdh;

  2489.         /** dU / dl. */
  2490.         private T dUdl;

  2492.         /** dU / dAlpha. */
  2493.         private T dUdAl;

  2495.         /** dU / dBeta. */
  2496.         private T dUdBe;

  2498.         /** dU / dGamma. */
  2499.         private T dUdGa;

  2501.         /** Simple constuctor.
  2502.          * @param date current date
  2503.          * @param context container for attributes
  2504.          * @param hansen hansen objects
  2505.          */
  2506.         FieldUAnddU(final FieldAbsoluteDate<T> date, final FieldDSSTTesseralContext<T> context,
  2507.                     final FieldHansenObjects<T> hansen) {

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

  2512.             // Zero for initialization
  2513.             final Field<T> field = date.getField();
  2514.             final T zero = field.getZero();

  2516.             // Potential derivatives
  2517.             dUda  = zero;
  2518.             dUdh  = zero;
  2519.             dUdk  = zero;
  2520.             dUdl  = zero;
  2521.             dUdAl = zero;
  2522.             dUdBe = zero;
  2523.             dUdGa = zero;

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

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

  2533.                 // R / a up to power degree
  2534.                 final T[] roaPow = MathArrays.buildArray(field, maxDegree + 1);
  2535.                 roaPow[0] = zero.add(1.);
  2536.                 for (int i = 1; i <= maxDegree; i++) {
  2537.                     roaPow[i] = roaPow[i - 1].multiply(context.getRoa());
  2538.                 }

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

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

  2546.                     // Phase angle
  2547.                     final T jlMmt  = auxiliaryElements.getLM().multiply(j).subtract(context.getTheta().multiply(m));
  2548.                     final T sinPhi = FastMath.sin(jlMmt);
  2549.                     final T cosPhi = FastMath.cos(jlMmt);

  2551.                     // Potential derivatives components for a given resonant pair {j,m}
  2552.                     T dUdaCos  = zero;
  2553.                     T dUdaSin  = zero;
  2554.                     T dUdhCos  = zero;
  2555.                     T dUdhSin  = zero;
  2556.                     T dUdkCos  = zero;
  2557.                     T dUdkSin  = zero;
  2558.                     T dUdlCos  = zero;
  2559.                     T dUdlSin  = zero;
  2560.                     T dUdAlCos = zero;
  2561.                     T dUdAlSin = zero;
  2562.                     T dUdBeCos = zero;
  2563.                     T dUdBeSin = zero;
  2564.                     T dUdGaCos = zero;
  2565.                     T dUdGaSin = zero;

  2567.                     // s-SUM from -sMin to sMax
  2568.                     final int sMin = FastMath.min(context.getMaxEccPow() - j, maxDegree);
  2569.                     final int sMax = FastMath.min(context.getMaxEccPow() + j, maxDegree);
  2570.                     for (int s = 0; s <= sMax; s++) {

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

  2575.                         // n-SUM for s positive
  2576.                         final T[][] nSumSpos = computeNSum(date, j, m, s, maxDegree,
  2577.                                                                 roaPow, ghMSJ, gammaMNS, context, hansen);
  2578.                         dUdaCos  = dUdaCos.add(nSumSpos[0][0]);
  2579.                         dUdaSin  = dUdaSin.add(nSumSpos[0][1]);
  2580.                         dUdhCos  = dUdhCos.add(nSumSpos[1][0]);
  2581.                         dUdhSin  = dUdhSin.add(nSumSpos[1][1]);
  2582.                         dUdkCos  = dUdkCos.add(nSumSpos[2][0]);
  2583.                         dUdkSin  = dUdkSin.add(nSumSpos[2][1]);
  2584.                         dUdlCos  = dUdlCos.add(nSumSpos[3][0]);
  2585.                         dUdlSin  = dUdlSin.add(nSumSpos[3][1]);
  2586.                         dUdAlCos = dUdAlCos.add(nSumSpos[4][0]);
  2587.                         dUdAlSin = dUdAlSin.add(nSumSpos[4][1]);
  2588.                         dUdBeCos = dUdBeCos.add(nSumSpos[5][0]);
  2589.                         dUdBeSin = dUdBeSin.add(nSumSpos[5][1]);
  2590.                         dUdGaCos = dUdGaCos.add(nSumSpos[6][0]);
  2591.                         dUdGaSin = dUdGaSin.add(nSumSpos[6][1]);

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

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

  2617.                     // Assembly of potential derivatives componants
  2618.                     dUda  = dUda.add(dUdaCos.multiply(cosPhi).add(dUdaSin.multiply(sinPhi)));
  2619.                     dUdh  = dUdh.add(dUdhCos.multiply(cosPhi).add(dUdhSin.multiply(sinPhi)));
  2620.                     dUdk  = dUdk.add(dUdkCos.multiply(cosPhi).add(dUdkSin.multiply(sinPhi)));
  2621.                     dUdl  = dUdl.add(dUdlCos.multiply(cosPhi).add(dUdlSin.multiply(sinPhi)));
  2622.                     dUdAl = dUdAl.add(dUdAlCos.multiply(cosPhi).add(dUdAlSin.multiply(sinPhi)));
  2623.                     dUdBe = dUdBe.add(dUdBeCos.multiply(cosPhi).add(dUdBeSin.multiply(sinPhi)));
  2624.                     dUdGa = dUdGa.add(dUdGaCos.multiply(cosPhi).add(dUdGaSin.multiply(sinPhi)));
  2625.                 }

  2627.                 dUda  =  dUda.multiply(context.getMoa().divide(auxiliaryElements.getSma())).negate();
  2628.                 dUdh  =  dUdh.multiply(context.getMoa());
  2629.                 dUdk  =  dUdk.multiply(context.getMoa());
  2630.                 dUdl  =  dUdl.multiply(context.getMoa());
  2631.                 dUdAl =  dUdAl.multiply(context.getMoa());
  2632.                 dUdBe =  dUdBe.multiply(context.getMoa());
  2633.                 dUdGa =  dUdGa.multiply(context.getMoa());

  2635.             }

  2637.         }

  2639.         /** Return value of dU / da.
  2640.          * @return dUda
  2641.          */
  2642.         public T getdUda() {
  2643.             return dUda;
  2644.         }

  2646.         /** Return value of dU / dk.
  2647.          * @return dUdk
  2648.          */
  2649.         public T getdUdk() {
  2650.             return dUdk;
  2651.         }

  2653.         /** Return value of dU / dh.
  2654.          * @return dUdh
  2655.          */
  2656.         public T getdUdh() {
  2657.             return dUdh;
  2658.         }

  2660.         /** Return value of dU / dl.
  2661.          * @return dUdl
  2662.          */
  2663.         public T getdUdl() {
  2664.             return dUdl;
  2665.         }

  2667.         /** Return value of dU / dAlpha.
  2668.          * @return dUdAl
  2669.          */
  2670.         public T getdUdAl() {
  2671.             return dUdAl;
  2672.         }

  2674.         /** Return value of dU / dBeta.
  2675.          * @return dUdBe
  2676.          */
  2677.         public T getdUdBe() {
  2678.             return dUdBe;
  2679.         }

  2681.         /** Return value of dU / dGamma.
  2682.          * @return dUdGa
  2683.          */
  2684.         public T getdUdGa() {
  2685.             return dUdGa;
  2686.         }

  2688.     }

  2690.     /** Computes init values of the Hansen Objects. */
  2691.     private class HansenObjects {

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

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

  2700.         /** Simple constructor.
  2701.          * @param ratio Ratio of satellite period to central body rotation period
  2702.          * @param maxEccPow Maximum power of the eccentricity to use in summation over s.
  2703.          * @param resOrders List of resonant orders
  2704.          * @param type type of the elements used during the propagation
  2705.          */
  2706.         HansenObjects(final double ratio,
  2707.                       final int maxEccPow,
  2708.                       final List<Integer> resOrders,
  2709.                       final PropagationType type) {

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

  2714.             //Allocate the two dimensional array
  2715.             final int rows     = 2 * maxDegree + 1;
  2716.             final int columns  = maxFrequencyShortPeriodics + 1;
  2717.             this.hansenObjects = new HansenTesseralLinear[rows][columns];

  2719.             switch (type) {
  2720.                 case MEAN:
  2721.                     // loop through the resonant orders
  2722.                     for (int m : resOrders) {
  2723.                         //Compute the corresponding j term
  2724.                         final int j = FastMath.max(1, (int) FastMath.round(ratio * m));

  2726.                         //Compute the sMin and sMax values
  2727.                         final int sMin = FastMath.min(maxEccPow - j, maxDegree);
  2728.                         final int sMax = FastMath.min(maxEccPow + j, maxDegree);

  2730.                         //loop through the s values
  2731.                         for (int s = 0; s <= sMax; s++) {
  2732.                             //Compute the n0 value
  2733.                             final int n0 = FastMath.max(FastMath.max(2, m), s);

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

  2738.                             if (s > 0 && s <= sMin) {
  2739.                                 //Also create the object for the pair j, -s
  2740.                                 this.hansenObjects[maxDegree - s][j] =  new HansenTesseralLinear(maxDegree, -s, j, n0, maxHansen);
  2741.                             }
  2742.                         }
  2743.                     }
  2744.                     break;

  2746.                 case OSCULATING:
  2747.                     // create all objects
  2748.                     for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
  2749.                         for (int s = -maxDegree; s <= maxDegree; s++) {
  2750.                             //Compute the n0 value
  2751.                             final int n0 = FastMath.max(2, FastMath.abs(s));
  2752.                             this.hansenObjects[s + maxDegree][j] = new HansenTesseralLinear(maxDegree, s, j, n0, maxHansen);
  2753.                         }
  2754.                     }
  2755.                     break;

  2757.                 default:
  2758.                     throw new OrekitInternalError(null);
  2759.             }

  2761.         }

  2763.         /** Compute init values for hansen objects.
  2764.          * @param context container for attributes
  2765.          * @param rows number of rows of the hansen matrix
  2766.          * @param columns columns number of columns of the hansen matrix
  2767.          */
  2768.         public void computeHansenObjectsInitValues(final DSSTTesseralContext context, final int rows, final int columns) {
  2769.             hansenObjects[rows][columns].computeInitValues(context.getE2(), context.getChi(), context.getChi2());
  2770.         }

  2772.         /** Get hansen object.
  2773.          * @return hansenObjects
  2774.          */
  2775.         public HansenTesseralLinear[][] getHansenObjects() {
  2776.             return hansenObjects;
  2777.         }

  2779.     }

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

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

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

  2791.         /** Simple constructor.
  2792.          * @param ratio Ratio of satellite period to central body rotation period
  2793.          * @param maxEccPow Maximum power of the eccentricity to use in summation over s.
  2794.          * @param resOrders List of resonant orders
  2795.          * @param type type of the elements used during the propagation
  2796.          * @param field field used by default.
  2797.          */
  2798.         @SuppressWarnings("unchecked")
  2799.         FieldHansenObjects(final T ratio,
  2800.                            final int maxEccPow,
  2801.                            final List<Integer> resOrders,
  2802.                            final PropagationType type,
  2803.                            final Field<T> field) {

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

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

  2813.             switch (type) {
  2814.                 case MEAN:
  2815.                  // loop through the resonant orders
  2816.                     for (int m : resOrders) {
  2817.                         //Compute the corresponding j term
  2818.                         final int j = FastMath.max(1, (int) FastMath.round(ratio.multiply(m)));

  2820.                         //Compute the sMin and sMax values
  2821.                         final int sMin = FastMath.min(maxEccPow - j, maxDegree);
  2822.                         final int sMax = FastMath.min(maxEccPow + j, maxDegree);

  2824.                         //loop through the s values
  2825.                         for (int s = 0; s <= sMax; s++) {
  2826.                             //Compute the n0 value
  2827.                             final int n0 = FastMath.max(FastMath.max(2, m), s);

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

  2832.                             if (s > 0 && s <= sMin) {
  2833.                                 //Also create the object for the pair j, -s
  2834.                                 this.hansenObjects[maxDegree - s][j] =  new FieldHansenTesseralLinear<>(maxDegree, -s, j, n0, maxHansen, field);
  2835.                             }
  2836.                         }
  2837.                     }
  2838.                     break;

  2840.                 case OSCULATING:
  2841.                     // create all objects
  2842.                     for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
  2843.                         for (int s = -maxDegree; s <= maxDegree; s++) {
  2844.                             //Compute the n0 value
  2845.                             final int n0 = FastMath.max(2, FastMath.abs(s));
  2846.                             this.hansenObjects[s + maxDegree][j] = new FieldHansenTesseralLinear<>(maxDegree, s, j, n0, maxHansen, field);
  2847.                         }
  2848.                     }
  2849.                     break;

  2851.                 default:
  2852.                     throw new OrekitInternalError(null);
  2853.             }

  2855.         }

  2857.         /** Compute init values for hansen objects.
  2858.          * @param context container for attributes
  2859.          * @param rows number of rows of the hansen matrix
  2860.          * @param columns columns number of columns of the hansen matrix
  2861.          */
  2862.         public void computeHansenObjectsInitValues(final FieldDSSTTesseralContext<T> context,
  2863.                                                    final int rows, final int columns) {
  2864.             hansenObjects[rows][columns].computeInitValues(context.getE2(), context.getChi(), context.getChi2());
  2865.         }

  2867.         /** Get hansen object.
  2868.          * @return hansenObjects
  2869.          */
  2870.         public FieldHansenTesseralLinear<T>[][] getHansenObjects() {
  2871.             return hansenObjects;
  2872.         }

  2874.     }

  2876. }
Created with JaCoCo 0.8.4.201905082037