DSSTTesseral.java

  1. /* Copyright 2002-2018 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.io.NotSerializableException;
  20. import java.io.Serializable;
  21. import java.util.ArrayList;
  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.SortedSet;
  29. import java.util.TreeMap;

  31. import org.hipparchus.analysis.differentiation.DSFactory;
  32. import org.hipparchus.analysis.differentiation.DerivativeStructure;
  33. import org.hipparchus.exception.LocalizedCoreFormats;
  34. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  35. import org.hipparchus.util.FastMath;
  36. import org.hipparchus.util.MathUtils;
  37. import org.orekit.attitudes.AttitudeProvider;
  38. import org.orekit.errors.OrekitException;
  39. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  40. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
  41. import org.orekit.frames.Frame;
  42. import org.orekit.frames.Transform;
  43. import org.orekit.orbits.Orbit;
  44. import org.orekit.propagation.SpacecraftState;
  45. import org.orekit.propagation.events.EventDetector;
  46. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  47. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
  48. import org.orekit.propagation.semianalytical.dsst.utilities.GHmsjPolynomials;
  49. import org.orekit.propagation.semianalytical.dsst.utilities.GammaMnsFunction;
  50. import org.orekit.propagation.semianalytical.dsst.utilities.JacobiPolynomials;
  51. import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
  52. import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenTesseralLinear;
  53. import org.orekit.time.AbsoluteDate;
  54. import org.orekit.utils.TimeSpanMap;

  56. /** Tesseral contribution to the central body gravitational perturbation.
  57.  *  <p>
  58.  *  Only resonant tesserals are considered.
  59.  *  </p>
  60.  *
  61.  *  @author Romain Di Costanzo
  62.  *  @author Pascal Parraud
  63.  */
  64. public class DSSTTesseral implements DSSTForceModel {

  66.     /** Minimum period for analytically averaged high-order resonant
  67.      *  central body spherical harmonics in seconds.
  68.      */
  69.     private static final double MIN_PERIOD_IN_SECONDS = 864000.;

  71.     /** Minimum period for analytically averaged high-order resonant
  72.      *  central body spherical harmonics in satellite revolutions.
  73.      */
  74.     private static final double MIN_PERIOD_IN_SAT_REV = 10.;

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

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

  94.     /** Provider for spherical harmonics. */
  95.     private final UnnormalizedSphericalHarmonicsProvider provider;

  97.     /** Central body rotating frame. */
  98.     private final Frame bodyFrame;

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

  103.     /** Central body rotation period (seconds). */
  104.     private final double bodyPeriod;

  106.     /** Maximal degree to consider for harmonics potential. */
  107.     private final int maxDegree;

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

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

  115.     /** Maximal order to consider for harmonics potential. */
  116.     private final int maxOrder;

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

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

  124.     /** List of resonant orders. */
  125.     private final List<Integer> resOrders;

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

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

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

  138.     /** Maximum value for j. */
  139.     private final int maxFrequencyShortPeriodics;

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

  144.     /** Keplerian period. */
  145.     private double orbitPeriod;

  147.     /** Ratio of satellite period to central body rotation period. */
  148.     private double ratio;

  150.     // Equinoctial elements (according to DSST notation)
  151.     /** a. */
  152.     private double a;

  154.     /** ex. */
  155.     private double k;

  157.     /** ey. */
  158.     private double h;

  160.     /** hx. */
  161.     private double q;

  163.     /** hy. */
  164.     private double p;

  166.     /** lm. */
  167.     private double lm;

  169.     /** Eccentricity. */
  170.     private double ecc;

  172.     // Common factors for potential computation
  173.     /** &Chi; = 1 / sqrt(1 - e²) = 1 / B. */
  174.     private double chi;

  176.     /** &Chi;². */
  177.     private double chi2;

  179.     // Equinoctial reference frame vectors (according to DSST notation)
  180.     /** Equinoctial frame f vector. */
  181.     private Vector3D f;

  183.     /** Equinoctial frame g vector. */
  184.     private Vector3D g;

  186.     /** Central body rotation angle θ. */
  187.     private double theta;

  189.     /** Direction cosine α. */
  190.     private double alpha;

  192.     /** Direction cosine β. */
  193.     private double beta;

  195.     /** Direction cosine γ. */
  196.     private double gamma;

  198.     // Common factors from equinoctial coefficients
  199.     /** 2 * a / A .*/
  200.     private double ax2oA;

  202.     /** 1 / (A * B) .*/
  203.     private double ooAB;

  205.     /** B / A .*/
  206.     private double BoA;

  208.     /** B / (A * (1 + B)) .*/
  209.     private double BoABpo;

  211.     /** C / (2 * A * B) .*/
  212.     private double Co2AB;

  214.     /** μ / a .*/
  215.     private double moa;

  217.     /** R / a .*/
  218.     private double roa;

  220.     /** ecc². */
  221.     private double e2;

  223.     /** The satellite mean motion. */
  224.     private double meanMotion;

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

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

  233.     /** Fourier coefficients. */
  234.     private FourierCjSjCoefficients cjsjFourier;

  236.     /** Short period terms. */
  237.     private TesseralShortPeriodicCoefficients shortPeriodTerms;

  239.     /** Factory for the DerivativeStructure instances. */
  240.     private final DSFactory factory;

  242.     /** Simple constructor.
  243.      * @param centralBodyFrame rotating body frame
  244.      * @param centralBodyRotationRate central body rotation rate (rad/s)
  245.      * @param provider provider for spherical harmonics
  246.      * @param maxDegreeTesseralSP maximal degree to consider for short periodics tesseral harmonics potential
  247.      *  (must be between 2 and {@code provider.getMaxDegree()})
  248.      * @param maxOrderTesseralSP maximal order to consider for short periodics tesseral harmonics potential
  249.      *  (must be between 0 and {@code provider.getMaxOrder()})
  250.      * @param maxEccPowTesseralSP maximum power of the eccentricity to use in summation over s for
  251.      * short periodic tesseral harmonics (without m-daily), should typically not exceed 4 as higher
  252.      * values will exceed computer capacity
  253.      * @param maxFrequencyShortPeriodics maximum frequency in mean longitude for short periodic computations
  254.      * (typically {@code maxDegreeTesseralSP} + {@code maxEccPowTesseralSP and no more than 12})
  255.      * @param maxDegreeMdailyTesseralSP maximal degree to consider for short periodics m-daily tesseral harmonics potential
  256.      *  (must be between 2 and {@code provider.getMaxDegree()})
  257.      * @param maxOrderMdailyTesseralSP maximal order to consider for short periodics m-daily tesseral harmonics potential
  258.      *  (must be between 0 and {@code provider.getMaxOrder()})
  259.      * @param maxEccPowMdailyTesseralSP maximum power of the eccentricity to use in summation over s for
  260.      * m-daily tesseral harmonics, (must be between 0 and {@code maxDegreeMdailyTesseralSP - 2},
  261.      * but should typically not exceed 4 as higher values will exceed computer capacity)
  262.      * @exception OrekitException if degrees or powers are out of range
  263.      * @since 7.2
  264.      */
  265.     public DSSTTesseral(final Frame centralBodyFrame,
  266.                         final double centralBodyRotationRate,
  267.                         final UnnormalizedSphericalHarmonicsProvider provider,
  268.                         final int maxDegreeTesseralSP, final int maxOrderTesseralSP,
  269.                         final int maxEccPowTesseralSP, final int maxFrequencyShortPeriodics,
  270.                         final int maxDegreeMdailyTesseralSP, final int maxOrderMdailyTesseralSP,
  271.                         final int maxEccPowMdailyTesseralSP)
  272.         throws OrekitException {

  274.         // Central body rotating frame
  275.         this.bodyFrame = centralBodyFrame;

  277.         //Save the rotation rate
  278.         this.centralBodyRotationRate = centralBodyRotationRate;

  280.         // Central body rotation period in seconds
  281.         this.bodyPeriod = MathUtils.TWO_PI / centralBodyRotationRate;

  283.         // Provider for spherical harmonics
  284.         this.provider  = provider;
  285.         this.maxDegree = provider.getMaxDegree();
  286.         this.maxOrder  = provider.getMaxOrder();

  288.         //set the maximum degree order for short periodics
  289.         checkIndexRange(maxDegreeTesseralSP, 2, maxDegree);
  290.         this.maxDegreeTesseralSP       = maxDegreeTesseralSP;

  292.         checkIndexRange(maxDegreeMdailyTesseralSP, 2, maxDegree);
  293.         this.maxDegreeMdailyTesseralSP = maxDegreeMdailyTesseralSP;

  295.         checkIndexRange(maxOrderTesseralSP, 0, maxOrder);
  296.         this.maxOrderTesseralSP        = maxOrderTesseralSP;

  298.         checkIndexRange(maxOrderMdailyTesseralSP, 0, maxOrder);
  299.         this.maxOrderMdailyTesseralSP  = maxOrderMdailyTesseralSP;

  301.         // set the maximum value for eccentricity power
  302.         this.maxEccPowTesseralSP       = maxEccPowTesseralSP;

  304.         checkIndexRange(maxEccPowMdailyTesseralSP, 0, maxDegreeMdailyTesseralSP - 2);
  305.         this.maxEccPowMdailyTesseralSP = maxEccPowMdailyTesseralSP;

  307.         // set the maximum value for frequency
  308.         this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;

  310.         // Initialize default values
  311.         this.resOrders = new ArrayList<Integer>();
  312.         this.nonResOrders = new TreeMap<Integer, List <Integer> >();
  313.         this.maxEccPow = 0;
  314.         this.maxHansen = 0;

  316.         this.factory = new DSFactory(1, 1);

  318.     }

  320.     /** Check an index range.
  321.      * @param index index value
  322.      * @param min minimum value for index
  323.      * @param max maximum value for index
  324.      * @exception OrekitException if index is out of range
  325.      */
  326.     private void checkIndexRange(final int index, final int min, final int max)
  327.         throws OrekitException {
  328.         if (index < min || index > max) {
  329.             throw new OrekitException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, index, min, max);
  330.         }
  331.     }

  333.     /** {@inheritDoc} */
  334.     @Override
  335.     public List<ShortPeriodTerms> initialize(final AuxiliaryElements aux, final boolean meanOnly)
  336.         throws OrekitException {

  338.         // Keplerian period
  339.         orbitPeriod = aux.getKeplerianPeriod();

  341.         // Set the highest power of the eccentricity in the analytical power
  342.         // series expansion for the averaged high order resonant central body
  343.         // spherical harmonic perturbation
  344.         final double e = aux.getEcc();
  345.         if (e <= 0.005) {
  346.             maxEccPow = 3;
  347.         } else if (e <= 0.02) {
  348.             maxEccPow = 4;
  349.         } else if (e <= 0.1) {
  350.             maxEccPow = 7;
  351.         } else if (e <= 0.2) {
  352.             maxEccPow = 10;
  353.         } else if (e <= 0.3) {
  354.             maxEccPow = 12;
  355.         } else if (e <= 0.4) {
  356.             maxEccPow = 15;
  357.         } else {
  358.             maxEccPow = 20;
  359.         }

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

  364.         // Ratio of satellite to central body periods to define resonant terms
  365.         ratio = orbitPeriod / bodyPeriod;

  367.         // Compute the resonant tesseral harmonic terms if not set by the user
  368.         getResonantAndNonResonantTerms(meanOnly);

  370.         //initialize the HansenTesseralLinear objects needed
  371.         createHansenObjects(meanOnly);

  373.         final int mMax = FastMath.max(maxOrderTesseralSP, maxOrderMdailyTesseralSP);
  374.         cjsjFourier = new FourierCjSjCoefficients(maxFrequencyShortPeriodics, mMax);

  376.         shortPeriodTerms = new TesseralShortPeriodicCoefficients(bodyFrame, maxOrderMdailyTesseralSP,
  377.                                                                  maxDegreeTesseralSP < 0, nonResOrders,
  378.                                                                  mMax, maxFrequencyShortPeriodics, INTERPOLATION_POINTS,
  379.                                                                  new TimeSpanMap<Slot>(new Slot(mMax, maxFrequencyShortPeriodics,
  380.                                                                                                 INTERPOLATION_POINTS)));

  382.         final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
  383.         list.add(shortPeriodTerms);
  384.         return list;

  386.     }

  388.     /** Create the objects needed for linear transformation.
  389.      *
  390.      * <p>
  391.      * Each {@link org.orekit.propagation.semianalytical.dsst.utilities.hansenHansenTesseralLinear HansenTesseralLinear} uses
  392.      * a fixed value for s and j. Since j varies from -maxJ to +maxJ and s varies from -maxDegree to +maxDegree,
  393.      * a 2 * maxDegree + 1 x 2 * maxJ + 1 matrix of objects should be created. The size of this matrix can be reduced
  394.      * by taking into account the expression (2.7.3-2). This means that is enough to create the objects for  positive
  395.      * values of j and all values of s.
  396.      * </p>
  397.      *
  398.      * @param meanOnly create only the objects required for the mean contribution
  399.      */
  400.     private void createHansenObjects(final boolean meanOnly) {
  401.         //Allocate the two dimensional array
  402.         final int rows     = 2 * maxDegree + 1;
  403.         final int columns  = maxFrequencyShortPeriodics + 1;
  404.         this.hansenObjects = new HansenTesseralLinear[rows][columns];

  406.         if (meanOnly) {
  407.             // loop through the resonant orders
  408.             for (int m : resOrders) {
  409.                 //Compute the corresponding j term
  410.                 final int j = FastMath.max(1, (int) FastMath.round(ratio * m));

  412.                 //Compute the sMin and sMax values
  413.                 final int sMin = FastMath.min(maxEccPow - j, maxDegree);
  414.                 final int sMax = FastMath.min(maxEccPow + j, maxDegree);

  416.                 //loop through the s values
  417.                 for (int s = 0; s <= sMax; s++) {
  418.                     //Compute the n0 value
  419.                     final int n0 = FastMath.max(FastMath.max(2, m), s);

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

  424.                     if (s > 0 && s <= sMin) {
  425.                         //Also create the object for the pair j, -s
  426.                         this.hansenObjects[maxDegree - s][j] =  new HansenTesseralLinear(maxDegree, -s, j, n0, maxHansen);
  427.                     }
  428.                 }
  429.             }
  430.         } else {
  431.             // create all objects
  432.             for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
  433.                 for (int s = -maxDegree; s <= maxDegree; s++) {
  434.                     //Compute the n0 value
  435.                     final int n0 = FastMath.max(2, FastMath.abs(s));

  437.                     this.hansenObjects[s + maxDegree][j] = new HansenTesseralLinear(maxDegree, s, j, n0, maxHansen);
  438.                 }
  439.             }
  440.         }
  441.     }

  443.     /** {@inheritDoc} */
  444.     @Override
  445.     public void initializeStep(final AuxiliaryElements aux) throws OrekitException {

  447.         // Equinoctial elements
  448.         a  = aux.getSma();
  449.         k  = aux.getK();
  450.         h  = aux.getH();
  451.         q  = aux.getQ();
  452.         p  = aux.getP();
  453.         lm = aux.getLM();

  455.         // Eccentricity
  456.         ecc = aux.getEcc();
  457.         e2 = ecc * ecc;

  459.         // Equinoctial frame vectors
  460.         f = aux.getVectorF();
  461.         g = aux.getVectorG();

  463.         // Central body rotation angle from equation 2.7.1-(3)(4).
  464.         final Transform t = bodyFrame.getTransformTo(aux.getFrame(), aux.getDate());
  465.         final Vector3D xB = t.transformVector(Vector3D.PLUS_I);
  466.         final Vector3D yB = t.transformVector(Vector3D.PLUS_J);
  467.         theta = FastMath.atan2(-f.dotProduct(yB) + I * g.dotProduct(xB),
  468.                                 f.dotProduct(xB) + I * g.dotProduct(yB));

  470.         // Direction cosines
  471.         alpha = aux.getAlpha();
  472.         beta  = aux.getBeta();
  473.         gamma = aux.getGamma();

  475.         // Equinoctial coefficients
  476.         // A = sqrt(μ * a)
  477.         final double A = aux.getA();
  478.         // B = sqrt(1 - h² - k²)
  479.         final double B = aux.getB();
  480.         // C = 1 + p² + q²
  481.         final double C = aux.getC();
  482.         // Common factors from equinoctial coefficients
  483.         // 2 * a / A
  484.         ax2oA  = 2. * a / A;
  485.         // B / A
  486.         BoA  = B / A;
  487.         // 1 / AB
  488.         ooAB = 1. / (A * B);
  489.         // C / 2AB
  490.         Co2AB = C * ooAB / 2.;
  491.         // B / (A * (1 + B))
  492.         BoABpo = BoA / (1. + B);
  493.         // &mu / a
  494.         moa = provider.getMu() / a;
  495.         // R / a
  496.         roa = provider.getAe() / a;

  498.         // &Chi; = 1 / B
  499.         chi = 1. / B;
  500.         chi2 = chi * chi;

  502.         //mean motion n
  503.         meanMotion = aux.getMeanMotion();
  504.     }

  506.     /** {@inheritDoc} */
  507.     @Override
  508.     public double[] getMeanElementRate(final SpacecraftState spacecraftState) throws OrekitException {

  510.         // Compute potential derivatives
  511.         final double[] dU  = computeUDerivatives(spacecraftState.getDate());
  512.         final double dUda  = dU[0];
  513.         final double dUdh  = dU[1];
  514.         final double dUdk  = dU[2];
  515.         final double dUdl  = dU[3];
  516.         final double dUdAl = dU[4];
  517.         final double dUdBe = dU[5];
  518.         final double dUdGa = dU[6];

  520.         // Compute the cross derivative operator :
  521.         final double UAlphaGamma   = alpha * dUdGa - gamma * dUdAl;
  522.         final double UAlphaBeta    = alpha * dUdBe - beta  * dUdAl;
  523.         final double UBetaGamma    =  beta * dUdGa - gamma * dUdBe;
  524.         final double Uhk           =     h * dUdk  -     k * dUdh;
  525.         final double pUagmIqUbgoAB = (p * UAlphaGamma - I * q * UBetaGamma) * ooAB;
  526.         final double UhkmUabmdUdl  =  Uhk - UAlphaBeta - dUdl;

  528.         final double da =  ax2oA * dUdl;
  529.         final double dh =    BoA * dUdk + k * pUagmIqUbgoAB - h * BoABpo * dUdl;
  530.         final double dk =  -(BoA * dUdh + h * pUagmIqUbgoAB + k * BoABpo * dUdl);
  531.         final double dp =  Co2AB * (p * UhkmUabmdUdl - UBetaGamma);
  532.         final double dq =  Co2AB * (q * UhkmUabmdUdl - I * UAlphaGamma);
  533.         final double dM = -ax2oA * dUda + BoABpo * (h * dUdh + k * dUdk) + pUagmIqUbgoAB;

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

  538.     /** {@inheritDoc} */
  539.     @Override
  540.     public void updateShortPeriodTerms(final SpacecraftState... meanStates)
  541.         throws OrekitException {

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

  545.         for (final SpacecraftState meanState : meanStates) {

  547.             initializeStep(new AuxiliaryElements(meanState.getOrbit(), I));

  549.             // Initialise the Hansen coefficients
  550.             for (int s = -maxDegree; s <= maxDegree; s++) {
  551.                 // coefficients with j == 0 are always needed
  552.                 this.hansenObjects[s + maxDegree][0].computeInitValues(e2, chi, chi2);
  553.                 if (maxDegreeTesseralSP >= 0) {
  554.                     // initialize other objects only if required
  555.                     for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  556.                         this.hansenObjects[s + maxDegree][j].computeInitValues(e2, chi, chi2);
  557.                     }
  558.                 }
  559.             }

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

  567.                 // the coefficient 3n / 2a
  568.                 final double tnota = 1.5 * meanMotion / a;

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

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

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

  588.         }

  590.     }

  592.     /** Build a set of coefficients.
  593.      *
  594.      * @param date the current date
  595.      * @param slot slot to which the coefficients belong
  596.      * @param m m index
  597.      * @param j j index
  598.      * @param tnota 3n/2a
  599.      */
  600.     private void buildCoefficients(final AbsoluteDate date, final Slot slot,
  601.                                    final int m, final int j, final double tnota) {
  602.         // Create local arrays
  603.         final double[] currentCijm = new double[] {0., 0., 0., 0., 0., 0.};
  604.         final double[] currentSijm = new double[] {0., 0., 0., 0., 0., 0.};

  606.         // compute the term 1 / (jn - mθ<sup>.</sup>)
  607.         final double oojnmt = 1. / (j * meanMotion - m * centralBodyRotationRate);

  609.         // initialise the coeficients
  610.         for (int i = 0; i < 6; i++) {
  611.             currentCijm[i] = -cjsjFourier.getSijm(i, j, m);
  612.             currentSijm[i] = cjsjFourier.getCijm(i, j, m);
  613.         }
  614.         // Add the separate part for δ<sub>6i</sub>
  615.         currentCijm[5] += tnota * oojnmt * cjsjFourier.getCijm(0, j, m);
  616.         currentSijm[5] += tnota * oojnmt * cjsjFourier.getSijm(0, j, m);

  618.         //Multiply by 1 / (jn - mθ<sup>.</sup>)
  619.         for (int i = 0; i < 6; i++) {
  620.             currentCijm[i] *= oojnmt;
  621.             currentSijm[i] *= oojnmt;
  622.         }

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

  628.     }

  630.     /** {@inheritDoc} */
  631.     @Override
  632.     public EventDetector[] getEventsDetectors() {
  633.         return null;
  634.     }

  636.      /**
  637.       * Get the resonant and non-resonant tesseral terms in the central body spherical harmonic field.
  638.       *
  639.       * @param resonantOnly extract only resonant terms
  640.       */
  641.     private void getResonantAndNonResonantTerms(final boolean resonantOnly) {

  643.         // Compute natural resonant terms
  644.         final double tolerance = 1. / FastMath.max(MIN_PERIOD_IN_SAT_REV,
  645.                                                    MIN_PERIOD_IN_SECONDS / orbitPeriod);

  647.         // Search the resonant orders in the tesseral harmonic field
  648.         resOrders.clear();
  649.         nonResOrders.clear();
  650.         for (int m = 1; m <= maxOrder; m++) {
  651.             final double resonance = ratio * m;
  652.             int jRes = 0;
  653.             final int jComputedRes = (int) FastMath.round(resonance);
  654.             if (jComputedRes > 0 && jComputedRes <= maxFrequencyShortPeriodics && FastMath.abs(resonance - jComputedRes) <= tolerance) {
  655.                 // Store each resonant index and order
  656.                 resOrders.add(m);
  657.                 jRes = jComputedRes;
  658.             }

  660.             if (!resonantOnly && maxDegreeTesseralSP >= 0 && m <= maxOrderTesseralSP) {
  661.                 //compute non resonant orders in the tesseral harmonic field
  662.                 final List<Integer> listJofM = new ArrayList<Integer>();
  663.                 //for the moment we take only the pairs (j,m) with |j| <= maxDegree + maxEccPow (from |s-j| <= maxEccPow and |s| <= maxDegree)
  664.                 for (int j = -maxFrequencyShortPeriodics; j <= maxFrequencyShortPeriodics; j++) {
  665.                     if (j != 0 && j != jRes) {
  666.                         listJofM.add(j);
  667.                     }
  668.                 }

  670.                 nonResOrders.put(m, listJofM);
  671.             }
  672.         }
  673.     }

  675.     /** Computes the potential U derivatives.
  676.      *  <p>The following elements are computed from expression 3.3 - (4).
  677.      *  <pre>
  678.      *  dU / da
  679.      *  dU / dh
  680.      *  dU / dk
  681.      *  dU / dλ
  682.      *  dU / dα
  683.      *  dU / dβ
  684.      *  dU / dγ
  685.      *  </pre>
  686.      *  </p>
  687.      *
  688.      *  @param date current date
  689.      *  @return potential derivatives
  690.      *  @throws OrekitException if an error occurs
  691.      */
  692.     private double[] computeUDerivatives(final AbsoluteDate date) throws OrekitException {

  694.         // Potential derivatives
  695.         double dUda  = 0.;
  696.         double dUdh  = 0.;
  697.         double dUdk  = 0.;
  698.         double dUdl  = 0.;
  699.         double dUdAl = 0.;
  700.         double dUdBe = 0.;
  701.         double dUdGa = 0.;

  703.         // Compute only if there is at least one resonant tesseral
  704.         if (!resOrders.isEmpty()) {
  705.             // Gmsj and Hmsj polynomials
  706.             final GHmsjPolynomials ghMSJ = new GHmsjPolynomials(k, h, alpha, beta, I);

  708.             // GAMMAmns function
  709.             final GammaMnsFunction gammaMNS = new GammaMnsFunction(maxDegree, gamma, I);

  711.             // R / a up to power degree
  712.             final double[] roaPow = new double[maxDegree + 1];
  713.             roaPow[0] = 1.;
  714.             for (int i = 1; i <= maxDegree; i++) {
  715.                 roaPow[i] = roa * roaPow[i - 1];
  716.             }

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

  721.                 // Resonant index for the current resonant order
  722.                 final int j = FastMath.max(1, (int) FastMath.round(ratio * m));

  724.                 // Phase angle
  725.                 final double jlMmt  = j * lm - m * theta;
  726.                 final double sinPhi = FastMath.sin(jlMmt);
  727.                 final double cosPhi = FastMath.cos(jlMmt);

  729.                 // Potential derivatives components for a given resonant pair {j,m}
  730.                 double dUdaCos  = 0.;
  731.                 double dUdaSin  = 0.;
  732.                 double dUdhCos  = 0.;
  733.                 double dUdhSin  = 0.;
  734.                 double dUdkCos  = 0.;
  735.                 double dUdkSin  = 0.;
  736.                 double dUdlCos  = 0.;
  737.                 double dUdlSin  = 0.;
  738.                 double dUdAlCos = 0.;
  739.                 double dUdAlSin = 0.;
  740.                 double dUdBeCos = 0.;
  741.                 double dUdBeSin = 0.;
  742.                 double dUdGaCos = 0.;
  743.                 double dUdGaSin = 0.;

  745.                 // s-SUM from -sMin to sMax
  746.                 final int sMin = FastMath.min(maxEccPow - j, maxDegree);
  747.                 final int sMax = FastMath.min(maxEccPow + j, maxDegree);
  748.                 for (int s = 0; s <= sMax; s++) {

  750.                     //Compute the initial values for Hansen coefficients using newComb operators
  751.                     this.hansenObjects[s + maxDegree][j].computeInitValues(e2, chi, chi2);

  753.                     // n-SUM for s positive
  754.                     final double[][] nSumSpos = computeNSum(date, j, m, s, maxDegree,
  755.                                                             roaPow, ghMSJ, gammaMNS);
  756.                     dUdaCos  += nSumSpos[0][0];
  757.                     dUdaSin  += nSumSpos[0][1];
  758.                     dUdhCos  += nSumSpos[1][0];
  759.                     dUdhSin  += nSumSpos[1][1];
  760.                     dUdkCos  += nSumSpos[2][0];
  761.                     dUdkSin  += nSumSpos[2][1];
  762.                     dUdlCos  += nSumSpos[3][0];
  763.                     dUdlSin  += nSumSpos[3][1];
  764.                     dUdAlCos += nSumSpos[4][0];
  765.                     dUdAlSin += nSumSpos[4][1];
  766.                     dUdBeCos += nSumSpos[5][0];
  767.                     dUdBeSin += nSumSpos[5][1];
  768.                     dUdGaCos += nSumSpos[6][0];
  769.                     dUdGaSin += nSumSpos[6][1];

  771.                     // n-SUM for s negative
  772.                     if (s > 0 && s <= sMin) {
  773.                         //Compute the initial values for Hansen coefficients using newComb operators
  774.                         this.hansenObjects[maxDegree - s][j].computeInitValues(e2, chi, chi2);

  776.                         final double[][] nSumSneg = computeNSum(date, j, m, -s, maxDegree,
  777.                                                                 roaPow, ghMSJ, gammaMNS);
  778.                         dUdaCos  += nSumSneg[0][0];
  779.                         dUdaSin  += nSumSneg[0][1];
  780.                         dUdhCos  += nSumSneg[1][0];
  781.                         dUdhSin  += nSumSneg[1][1];
  782.                         dUdkCos  += nSumSneg[2][0];
  783.                         dUdkSin  += nSumSneg[2][1];
  784.                         dUdlCos  += nSumSneg[3][0];
  785.                         dUdlSin  += nSumSneg[3][1];
  786.                         dUdAlCos += nSumSneg[4][0];
  787.                         dUdAlSin += nSumSneg[4][1];
  788.                         dUdBeCos += nSumSneg[5][0];
  789.                         dUdBeSin += nSumSneg[5][1];
  790.                         dUdGaCos += nSumSneg[6][0];
  791.                         dUdGaSin += nSumSneg[6][1];
  792.                     }
  793.                 }

  795.                 // Assembly of potential derivatives componants
  796.                 dUda  += cosPhi * dUdaCos  + sinPhi * dUdaSin;
  797.                 dUdh  += cosPhi * dUdhCos  + sinPhi * dUdhSin;
  798.                 dUdk  += cosPhi * dUdkCos  + sinPhi * dUdkSin;
  799.                 dUdl  += cosPhi * dUdlCos  + sinPhi * dUdlSin;
  800.                 dUdAl += cosPhi * dUdAlCos + sinPhi * dUdAlSin;
  801.                 dUdBe += cosPhi * dUdBeCos + sinPhi * dUdBeSin;
  802.                 dUdGa += cosPhi * dUdGaCos + sinPhi * dUdGaSin;
  803.             }

  805.             dUda  *= -moa / a;
  806.             dUdh  *=  moa;
  807.             dUdk  *=  moa;
  808.             dUdl  *=  moa;
  809.             dUdAl *=  moa;
  810.             dUdBe *=  moa;
  811.             dUdGa *=  moa;
  812.         }

  814.         return new double[] {dUda, dUdh, dUdk, dUdl, dUdAl, dUdBe, dUdGa};
  815.     }

  817.     /** Compute the n-SUM for potential derivatives components.
  818.      *  @param date current date
  819.      *  @param j resonant index <i>j</i>
  820.      *  @param m resonant order <i>m</i>
  821.      *  @param s d'Alembert characteristic <i>s</i>
  822.      *  @param maxN maximum possible value for <i>n</i> index
  823.      *  @param roaPow powers of R/a up to degree <i>n</i>
  824.      *  @param ghMSJ G<sup>j</sup><sub>m,s</sub> and H<sup>j</sup><sub>m,s</sub> polynomials
  825.      *  @param gammaMNS &Gamma;<sup>m</sup><sub>n,s</sub>(γ) function
  826.      *  @return Components of U<sub>n</sub> derivatives for fixed j, m, s
  827.      * @throws OrekitException if some error occurred
  828.      */
  829.     private double[][] computeNSum(final AbsoluteDate date,
  830.                                    final int j, final int m, final int s, final int maxN, final double[] roaPow,
  831.                                    final GHmsjPolynomials ghMSJ, final GammaMnsFunction gammaMNS)
  832.         throws OrekitException {

  834.         //spherical harmonics
  835.         final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);

  837.         // Potential derivatives components
  838.         double dUdaCos  = 0.;
  839.         double dUdaSin  = 0.;
  840.         double dUdhCos  = 0.;
  841.         double dUdhSin  = 0.;
  842.         double dUdkCos  = 0.;
  843.         double dUdkSin  = 0.;
  844.         double dUdlCos  = 0.;
  845.         double dUdlSin  = 0.;
  846.         double dUdAlCos = 0.;
  847.         double dUdAlSin = 0.;
  848.         double dUdBeCos = 0.;
  849.         double dUdBeSin = 0.;
  850.         double dUdGaCos = 0.;
  851.         double dUdGaSin = 0.;

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

  857.         // jacobi v, w, indices from 2.7.1-(15)
  858.         final int v = FastMath.abs(m - s);
  859.         final int w = FastMath.abs(m + s);

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

  864.         //Get the corresponding Hansen object
  865.         final int sIndex = maxDegree + (j < 0 ? -s : s);
  866.         final int jIndex = FastMath.abs(j);
  867.         final HansenTesseralLinear hans = this.hansenObjects[sIndex][jIndex];

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

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

  877.                 // Inclination function Gamma and derivative
  878.                 final double gaMNS  = gammaMNS.getValue(m, n, s);
  879.                 final double dGaMNS = gammaMNS.getDerivative(m, n, s);

  881.                 // Hansen kernel value and derivative
  882.                 final double kJNS   = hans.getValue(-n - 1, chi);
  883.                 final double dkJNS  = hans.getDerivative(-n - 1, chi);

  885.                 // Gjms, Hjms polynomials and derivatives
  886.                 final double gMSJ   = ghMSJ.getGmsj(m, s, j);
  887.                 final double hMSJ   = ghMSJ.getHmsj(m, s, j);
  888.                 final double dGdh   = ghMSJ.getdGmsdh(m, s, j);
  889.                 final double dGdk   = ghMSJ.getdGmsdk(m, s, j);
  890.                 final double dGdA   = ghMSJ.getdGmsdAlpha(m, s, j);
  891.                 final double dGdB   = ghMSJ.getdGmsdBeta(m, s, j);
  892.                 final double dHdh   = ghMSJ.getdHmsdh(m, s, j);
  893.                 final double dHdk   = ghMSJ.getdHmsdk(m, s, j);
  894.                 final double dHdA   = ghMSJ.getdHmsdAlpha(m, s, j);
  895.                 final double dHdB   = ghMSJ.getdHmsdBeta(m, s, j);

  897.                 // Jacobi l-index from 2.7.1-(15)
  898.                 final int l = FastMath.min(n - m, n - FastMath.abs(s));
  899.                 // Jacobi polynomial and derivative
  900.                 final DerivativeStructure jacobi =
  901.                         JacobiPolynomials.getValue(l, v, w, factory.variable(0, gamma));

  903.                 // Geopotential coefficients
  904.                 final double cnm = harmonics.getUnnormalizedCnm(n, m);
  905.                 final double snm = harmonics.getUnnormalizedSnm(n, m);

  907.                 // Common factors from expansion of equations 3.3-4
  908.                 final double cf_0      = roaPow[n] * Im * vMNS;
  909.                 final double cf_1      = cf_0 * gaMNS * jacobi.getValue();
  910.                 final double cf_2      = cf_1 * kJNS;
  911.                 final double gcPhs     = gMSJ * cnm + hMSJ * snm;
  912.                 final double gsMhc     = gMSJ * snm - hMSJ * cnm;
  913.                 final double dKgcPhsx2 = 2. * dkJNS * gcPhs;
  914.                 final double dKgsMhcx2 = 2. * dkJNS * gsMhc;
  915.                 final double dUdaCoef  = (n + 1) * cf_2;
  916.                 final double dUdlCoef  = j * cf_2;
  917.                 final double dUdGaCoef = cf_0 * kJNS * (jacobi.getValue() * dGaMNS + gaMNS * jacobi.getPartialDerivative(1));

  919.                 // dU / da components
  920.                 dUdaCos  += dUdaCoef * gcPhs;
  921.                 dUdaSin  += dUdaCoef * gsMhc;

  923.                 // dU / dh components
  924.                 dUdhCos  += cf_1 * (kJNS * (cnm * dGdh + snm * dHdh) + h * dKgcPhsx2);
  925.                 dUdhSin  += cf_1 * (kJNS * (snm * dGdh - cnm * dHdh) + h * dKgsMhcx2);

  927.                 // dU / dk components
  928.                 dUdkCos  += cf_1 * (kJNS * (cnm * dGdk + snm * dHdk) + k * dKgcPhsx2);
  929.                 dUdkSin  += cf_1 * (kJNS * (snm * dGdk - cnm * dHdk) + k * dKgsMhcx2);

  931.                 // dU / dLambda components
  932.                 dUdlCos  +=  dUdlCoef * gsMhc;
  933.                 dUdlSin  += -dUdlCoef * gcPhs;

  935.                 // dU / alpha components
  936.                 dUdAlCos += cf_2 * (dGdA * cnm + dHdA * snm);
  937.                 dUdAlSin += cf_2 * (dGdA * snm - dHdA * cnm);

  939.                 // dU / dBeta components
  940.                 dUdBeCos += cf_2 * (dGdB * cnm + dHdB * snm);
  941.                 dUdBeSin += cf_2 * (dGdB * snm - dHdB * cnm);

  943.                 // dU / dGamma components
  944.                 dUdGaCos += dUdGaCoef * gcPhs;
  945.                 dUdGaSin += dUdGaCoef * gsMhc;
  946.             }
  947.         }

  949.         return new double[][] {{dUdaCos,  dUdaSin},
  950.                                {dUdhCos,  dUdhSin},
  951.                                {dUdkCos,  dUdkSin},
  952.                                {dUdlCos,  dUdlSin},
  953.                                {dUdAlCos, dUdAlSin},
  954.                                {dUdBeCos, dUdBeSin},
  955.                                {dUdGaCos, dUdGaSin}};
  956.     }

  958.     /** {@inheritDoc} */
  959.     @Override
  960.     public void registerAttitudeProvider(final AttitudeProvider attitudeProvider) {
  961.         //nothing is done since this contribution is not sensitive to attitude
  962.     }

  964.     /** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
  965.      *  <p>
  966.      *  Those coefficients are given in Danielson paper by substituting the
  967.      *  disturbing function (2.7.1-16) with m != 0 into (2.2-10)
  968.      *  </p>
  969.      */
  970.     private class FourierCjSjCoefficients {

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

  975.         /** The C<sub>i</sub><sup>jm</sup> coefficients.
  976.          * <p>
  977.          * The index order is [m][j][i] <br/>
  978.          * The i index corresponds to the C<sub>i</sub><sup>jm</sup> coefficients used to
  979.          * compute the following: <br/>
  980.          * - da/dt <br/>
  981.          * - dk/dt <br/>
  982.          * - dh/dt / dk <br/>
  983.          * - dq/dt <br/>
  984.          * - dp/dt / dα <br/>
  985.          * - dλ/dt / dβ <br/>
  986.          * </p>
  987.          */
  988.         private final double[][][] cCoef;

  990.         /** The S<sub>i</sub><sup>jm</sup> coefficients.
  991.          * <p>
  992.          * The index order is [m][j][i] <br/>
  993.          * The i index corresponds to the C<sub>i</sub><sup>jm</sup> coefficients used to
  994.          * compute the following: <br/>
  995.          * - da/dt <br/>
  996.          * - dk/dt <br/>
  997.          * - dh/dt / dk <br/>
  998.          * - dq/dt <br/>
  999.          * - dp/dt / dα <br/>
  1000.          * - dλ/dt / dβ <br/>
  1001.          * </p>
  1002.          */
  1003.         private final double[][][] sCoef;

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

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

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

  1014.         /** Create a set of C<sub>i</sub><sup>jm</sup> and S<sub>i</sub><sup>jm</sup> coefficients.
  1015.          *  @param jMax absolute limit for j ( -jMax <= j <= jMax )
  1016.          *  @param mMax maximum value for m
  1017.          */
  1018.         FourierCjSjCoefficients(final int jMax, final int mMax) {
  1019.             // initialise fields
  1020.             final int rows    = mMax + 1;
  1021.             final int columns = 2 * jMax + 1;
  1022.             this.jMax         = jMax;
  1023.             this.cCoef        = new double[rows][columns][6];
  1024.             this.sCoef        = new double[rows][columns][6];
  1025.             this.roaPow       = new double[maxDegree + 1];
  1026.             roaPow[0] = 1.;
  1027.         }

  1029.         /**
  1030.          * Generate the coefficients.
  1031.          * @param date the current date
  1032.          * @throws OrekitException if an error occurs while generating the coefficients
  1033.          */
  1034.         public void generateCoefficients(final AbsoluteDate date) throws OrekitException {
  1035.             // Compute only if there is at least one non-resonant tesseral
  1036.             if (!nonResOrders.isEmpty() || maxDegreeTesseralSP < 0) {
  1037.                 // Gmsj and Hmsj polynomials
  1038.                 ghMSJ = new GHmsjPolynomials(k, h, alpha, beta, I);

  1040.                 // GAMMAmns function
  1041.                 gammaMNS = new GammaMnsFunction(maxDegree, gamma, I);

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

  1045.                 // R / a up to power degree
  1046.                 for (int i = 1; i <= maxRoaPower; i++) {
  1047.                     roaPow[i] = roa * roaPow[i - 1];
  1048.                 }

  1050.                 //generate the m-daily coefficients
  1051.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1052.                     buildFourierCoefficients(date, m, 0, maxDegreeMdailyTesseralSP);
  1053.                 }

  1055.                 // generate the other coefficients only if required
  1056.                 if (maxDegreeTesseralSP >= 0) {
  1057.                     for (int m: nonResOrders.keySet()) {
  1058.                         final List<Integer> listJ = nonResOrders.get(m);

  1060.                         for (int j: listJ) {
  1061.                             buildFourierCoefficients(date, m, j, maxDegreeTesseralSP);
  1062.                         }
  1063.                     }
  1064.                 }
  1065.             }
  1066.         }

  1068.         /** Build a set of fourier coefficients for a given m and j.
  1069.          *
  1070.          * @param date the date of the coefficients
  1071.          * @param m m index
  1072.          * @param j j index
  1073.          * @param maxN  maximum value for n index
  1074.          * @throws OrekitException in case of Hansen kernel generation error
  1075.          */
  1076.         private void buildFourierCoefficients(final AbsoluteDate date,
  1077.                final int m, final int j, final int maxN) throws OrekitException {
  1078.             // Potential derivatives components for a given non-resonant pair {j,m}
  1079.             double dRdaCos  = 0.;
  1080.             double dRdaSin  = 0.;
  1081.             double dRdhCos  = 0.;
  1082.             double dRdhSin  = 0.;
  1083.             double dRdkCos  = 0.;
  1084.             double dRdkSin  = 0.;
  1085.             double dRdlCos  = 0.;
  1086.             double dRdlSin  = 0.;
  1087.             double dRdAlCos = 0.;
  1088.             double dRdAlSin = 0.;
  1089.             double dRdBeCos = 0.;
  1090.             double dRdBeSin = 0.;
  1091.             double dRdGaCos = 0.;
  1092.             double dRdGaSin = 0.;

  1094.             // s-SUM from -sMin to sMax
  1095.             final int sMin = j == 0 ? maxEccPowMdailyTesseralSP : maxEccPowTesseralSP;
  1096.             final int sMax = j == 0 ? maxEccPowMdailyTesseralSP : maxEccPowTesseralSP;
  1097.             for (int s = 0; s <= sMax; s++) {

  1099.                 // n-SUM for s positive
  1100.                 final double[][] nSumSpos = computeNSum(date, j, m, s, maxN,
  1101.                                                         roaPow, ghMSJ, gammaMNS);
  1102.                 dRdaCos  += nSumSpos[0][0];
  1103.                 dRdaSin  += nSumSpos[0][1];
  1104.                 dRdhCos  += nSumSpos[1][0];
  1105.                 dRdhSin  += nSumSpos[1][1];
  1106.                 dRdkCos  += nSumSpos[2][0];
  1107.                 dRdkSin  += nSumSpos[2][1];
  1108.                 dRdlCos  += nSumSpos[3][0];
  1109.                 dRdlSin  += nSumSpos[3][1];
  1110.                 dRdAlCos += nSumSpos[4][0];
  1111.                 dRdAlSin += nSumSpos[4][1];
  1112.                 dRdBeCos += nSumSpos[5][0];
  1113.                 dRdBeSin += nSumSpos[5][1];
  1114.                 dRdGaCos += nSumSpos[6][0];
  1115.                 dRdGaSin += nSumSpos[6][1];

  1117.                 // n-SUM for s negative
  1118.                 if (s > 0 && s <= sMin) {
  1119.                     final double[][] nSumSneg = computeNSum(date, j, m, -s, maxN,
  1120.                                                             roaPow, ghMSJ, gammaMNS);
  1121.                     dRdaCos  += nSumSneg[0][0];
  1122.                     dRdaSin  += nSumSneg[0][1];
  1123.                     dRdhCos  += nSumSneg[1][0];
  1124.                     dRdhSin  += nSumSneg[1][1];
  1125.                     dRdkCos  += nSumSneg[2][0];
  1126.                     dRdkSin  += nSumSneg[2][1];
  1127.                     dRdlCos  += nSumSneg[3][0];
  1128.                     dRdlSin  += nSumSneg[3][1];
  1129.                     dRdAlCos += nSumSneg[4][0];
  1130.                     dRdAlSin += nSumSneg[4][1];
  1131.                     dRdBeCos += nSumSneg[5][0];
  1132.                     dRdBeSin += nSumSneg[5][1];
  1133.                     dRdGaCos += nSumSneg[6][0];
  1134.                     dRdGaSin += nSumSneg[6][1];
  1135.                 }
  1136.             }
  1137.             dRdaCos  *= -moa / a;
  1138.             dRdaSin  *= -moa / a;
  1139.             dRdhCos  *=  moa;
  1140.             dRdhSin  *=  moa;
  1141.             dRdkCos  *=  moa;
  1142.             dRdkSin *=  moa;
  1143.             dRdlCos *=  moa;
  1144.             dRdlSin *=  moa;
  1145.             dRdAlCos *=  moa;
  1146.             dRdAlSin *=  moa;
  1147.             dRdBeCos *=  moa;
  1148.             dRdBeSin *=  moa;
  1149.             dRdGaCos *=  moa;
  1150.             dRdGaSin *=  moa;

  1152.             // Compute the cross derivative operator :
  1153.             final double RAlphaGammaCos   = alpha * dRdGaCos - gamma * dRdAlCos;
  1154.             final double RAlphaGammaSin   = alpha * dRdGaSin - gamma * dRdAlSin;
  1155.             final double RAlphaBetaCos    = alpha * dRdBeCos - beta  * dRdAlCos;
  1156.             final double RAlphaBetaSin    = alpha * dRdBeSin - beta  * dRdAlSin;
  1157.             final double RBetaGammaCos    =  beta * dRdGaCos - gamma * dRdBeCos;
  1158.             final double RBetaGammaSin    =  beta * dRdGaSin - gamma * dRdBeSin;
  1159.             final double RhkCos           =     h * dRdkCos  -     k * dRdhCos;
  1160.             final double RhkSin           =     h * dRdkSin  -     k * dRdhSin;
  1161.             final double pRagmIqRbgoABCos = (p * RAlphaGammaCos - I * q * RBetaGammaCos) * ooAB;
  1162.             final double pRagmIqRbgoABSin = (p * RAlphaGammaSin - I * q * RBetaGammaSin) * ooAB;
  1163.             final double RhkmRabmdRdlCos  =  RhkCos - RAlphaBetaCos - dRdlCos;
  1164.             final double RhkmRabmdRdlSin  =  RhkSin - RAlphaBetaSin - dRdlSin;

  1166.             // da/dt
  1167.             cCoef[m][j + jMax][0] = ax2oA * dRdlCos;
  1168.             sCoef[m][j + jMax][0] = ax2oA * dRdlSin;

  1170.             // dk/dt
  1171.             cCoef[m][j + jMax][1] = -(BoA * dRdhCos + h * pRagmIqRbgoABCos + k * BoABpo * dRdlCos);
  1172.             sCoef[m][j + jMax][1] = -(BoA * dRdhSin + h * pRagmIqRbgoABSin + k * BoABpo * dRdlSin);

  1174.             // dh/dt
  1175.             cCoef[m][j + jMax][2] = BoA * dRdkCos + k * pRagmIqRbgoABCos - h * BoABpo * dRdlCos;
  1176.             sCoef[m][j + jMax][2] = BoA * dRdkSin + k * pRagmIqRbgoABSin - h * BoABpo * dRdlSin;

  1178.             // dq/dt
  1179.             cCoef[m][j + jMax][3] = Co2AB * (q * RhkmRabmdRdlCos - I * RAlphaGammaCos);
  1180.             sCoef[m][j + jMax][3] = Co2AB * (q * RhkmRabmdRdlSin - I * RAlphaGammaSin);

  1182.             // dp/dt
  1183.             cCoef[m][j + jMax][4] = Co2AB * (p * RhkmRabmdRdlCos - RBetaGammaCos);
  1184.             sCoef[m][j + jMax][4] = Co2AB * (p * RhkmRabmdRdlSin - RBetaGammaSin);

  1186.             // dλ/dt
  1187.             cCoef[m][j + jMax][5] = -ax2oA * dRdaCos + BoABpo * (h * dRdhCos + k * dRdkCos) + pRagmIqRbgoABCos;
  1188.             sCoef[m][j + jMax][5] = -ax2oA * dRdaSin + BoABpo * (h * dRdhSin + k * dRdkSin) + pRagmIqRbgoABSin;
  1189.         }

  1191.         /** Get the coefficient C<sub>i</sub><sup>jm</sup>.
  1192.          * @param i i index - corresponds to the required variation
  1193.          * @param j j index
  1194.          * @param m m index
  1195.          * @return the coefficient C<sub>i</sub><sup>jm</sup>
  1196.          */
  1197.         public double getCijm(final int i, final int j, final int m) {
  1198.             return cCoef[m][j + jMax][i];
  1199.         }

  1201.         /** Get the coefficient S<sub>i</sub><sup>jm</sup>.
  1202.          * @param i i index - corresponds to the required variation
  1203.          * @param j j index
  1204.          * @param m m index
  1205.          * @return the coefficient S<sub>i</sub><sup>jm</sup>
  1206.          */
  1207.         public double getSijm(final int i, final int j, final int m) {
  1208.             return sCoef[m][j + jMax][i];
  1209.         }
  1210.     }

  1212.     /** 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
  1213.      * the short-periodic zonal contribution.
  1214.      *   <p>
  1215.      *  Those coefficients are given by expression 2.5.4-5 from the Danielson paper.
  1216.      *   </p>
  1217.      *
  1218.      * @author Sorin Scortan
  1219.      *
  1220.      */
  1221.     private static class TesseralShortPeriodicCoefficients implements ShortPeriodTerms {

  1223.         /** Serializable UID. */
  1224.         private static final long serialVersionUID = 20151119L;

  1226.         /** Retrograde factor I.
  1227.          *  <p>
  1228.          *  DSST model needs equinoctial orbit as internal representation.
  1229.          *  Classical equinoctial elements have discontinuities when inclination
  1230.          *  is close to zero. In this representation, I = +1. <br>
  1231.          *  To avoid this discontinuity, another representation exists and equinoctial
  1232.          *  elements can be expressed in a different way, called "retrograde" orbit.
  1233.          *  This implies I = -1. <br>
  1234.          *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  1235.          *  But for the sake of consistency with the theory, the retrograde factor
  1236.          *  has been kept in the formulas.
  1237.          *  </p>
  1238.          */
  1239.         private static final int I = 1;

  1241.         /** Central body rotating frame. */
  1242.         private final Frame bodyFrame;

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

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

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

  1253.         /** Maximum value for m index. */
  1254.         private final int mMax;

  1256.         /** Maximum value for j index. */
  1257.         private final int jMax;

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

  1262.         /** All coefficients slots. */
  1263.         private final transient TimeSpanMap<Slot> slots;

  1265.         /** Constructor.
  1266.          * @param bodyFrame central body rotating frame
  1267.          * @param maxOrderMdailyTesseralSP maximal order to consider for short periodics m-daily tesseral harmonics potential
  1268.          * @param mDailiesOnly flag to take into account only M-dailies harmonic tesserals for short periodic perturbations
  1269.          * @param nonResOrders lst of non resonant orders with j != 0
  1270.          * @param mMax maximum value for m index
  1271.          * @param jMax maximum value for j index
  1272.          * @param interpolationPoints number of points used in the interpolation process
  1273.          * @param slots all coefficients slots
  1274.          */
  1275.         TesseralShortPeriodicCoefficients(final Frame bodyFrame, final int maxOrderMdailyTesseralSP,
  1276.                                           final boolean mDailiesOnly, final SortedMap<Integer, List<Integer> > nonResOrders,
  1277.                                           final int mMax, final int jMax, final int interpolationPoints,
  1278.                                           final TimeSpanMap<Slot> slots) {
  1279.             this.bodyFrame                = bodyFrame;
  1280.             this.maxOrderMdailyTesseralSP = maxOrderMdailyTesseralSP;
  1281.             this.mDailiesOnly             = mDailiesOnly;
  1282.             this.nonResOrders             = nonResOrders;
  1283.             this.mMax                     = mMax;
  1284.             this.jMax                     = jMax;
  1285.             this.interpolationPoints      = interpolationPoints;
  1286.             this.slots                    = slots;
  1287.         }

  1289.         /** Get the slot valid for some date.
  1290.          * @param meanStates mean states defining the slot
  1291.          * @return slot valid at the specified date
  1292.          */
  1293.         public Slot createSlot(final SpacecraftState... meanStates) {
  1294.             final Slot         slot  = new Slot(mMax, jMax, interpolationPoints);
  1295.             final AbsoluteDate first = meanStates[0].getDate();
  1296.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  1297.             if (first.compareTo(last) <= 0) {
  1298.                 slots.addValidAfter(slot, first);
  1299.             } else {
  1300.                 slots.addValidBefore(slot, first);
  1301.             }
  1302.             return slot;
  1303.         }

  1305.         /** {@inheritDoc} */
  1306.         @Override
  1307.         public double[] value(final Orbit meanOrbit) throws OrekitException {

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

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

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

  1319.                 //Build an auxiliary object
  1320.                 final AuxiliaryElements aux = new AuxiliaryElements(meanOrbit, I);

  1322.                 // Central body rotation angle from equation 2.7.1-(3)(4).
  1323.                 final Transform t = bodyFrame.getTransformTo(aux.getFrame(), aux.getDate());
  1324.                 final Vector3D xB = t.transformVector(Vector3D.PLUS_I);
  1325.                 final Vector3D yB = t.transformVector(Vector3D.PLUS_J);
  1326.                 final Vector3D  f = aux.getVectorF();
  1327.                 final Vector3D  g = aux.getVectorG();
  1328.                 final double currentTheta = FastMath.atan2(-f.dotProduct(yB) + I * g.dotProduct(xB),
  1329.                                                             f.dotProduct(xB) + I * g.dotProduct(yB));

  1331.                 //Add the m-daily contribution
  1332.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1333.                     // Phase angle
  1334.                     final double jlMmt  = -m * currentTheta;
  1335.                     final double sinPhi = FastMath.sin(jlMmt);
  1336.                     final double cosPhi = FastMath.cos(jlMmt);

  1338.                     // compute contribution for each element
  1339.                     final double[] c = slot.getCijm(0, m, meanOrbit.getDate());
  1340.                     final double[] s = slot.getSijm(0, m, meanOrbit.getDate());
  1341.                     for (int i = 0; i < 6; i++) {
  1342.                         shortPeriodicVariation[i] += c[i] * cosPhi + s[i] * sinPhi;
  1343.                     }
  1344.                 }

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

  1351.                     for (int j : listJ) {
  1352.                         // Phase angle
  1353.                         final double jlMmt  = j * meanOrbit.getLM() - m * currentTheta;
  1354.                         final double sinPhi = FastMath.sin(jlMmt);
  1355.                         final double cosPhi = FastMath.cos(jlMmt);

  1357.                         // compute contribution for each element
  1358.                         final double[] c = slot.getCijm(j, m, meanOrbit.getDate());
  1359.                         final double[] s = slot.getSijm(j, m, meanOrbit.getDate());
  1360.                         for (int i = 0; i < 6; i++) {
  1361.                             shortPeriodicVariation[i] += c[i] * cosPhi + s[i] * sinPhi;
  1362.                         }

  1364.                     }
  1365.                 }
  1366.             }

  1368.             return shortPeriodicVariation;

  1370.         }

  1372.         /** {@inheritDoc} */
  1373.         @Override
  1374.         public String getCoefficientsKeyPrefix() {
  1375.             return "DSST-central-body-tesseral-";
  1376.         }

  1378.         /** {@inheritDoc}
  1379.          * <p>
  1380.          * For tesseral terms contributions,there are maxOrderMdailyTesseralSP
  1381.          * m-daily cMm coefficients, maxOrderMdailyTesseralSP m-daily sMm
  1382.          * coefficients, nbNonResonant cjm coefficients and nbNonResonant
  1383.          * sjm coefficients, where maxOrderMdailyTesseralSP and nbNonResonant both
  1384.          * depend on the orbit. The j index is the integer multiplier for the true
  1385.          * longitude argument and the m index is the integer multiplier for m-dailies.
  1386.          * </p>
  1387.          */
  1388.         @Override
  1389.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected)
  1390.             throws OrekitException {

  1392.             // select the coefficients slot
  1393.             final Slot slot = slots.get(date);

  1395.             if (!nonResOrders.isEmpty() || mDailiesOnly) {
  1396.                 final Map<String, double[]> coefficients =
  1397.                                 new HashMap<String, double[]>(12 * maxOrderMdailyTesseralSP +
  1398.                                                               12 * nonResOrders.size());

  1400.                 for (int m = 1; m <= maxOrderMdailyTesseralSP; m++) {
  1401.                     storeIfSelected(coefficients, selected, slot.getCijm(0, m, date), "cM", m);
  1402.                     storeIfSelected(coefficients, selected, slot.getSijm(0, m, date), "sM", m);
  1403.                 }

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

  1409.                     for (int j : listJ) {
  1410.                         for (int i = 0; i < 6; ++i) {
  1411.                             storeIfSelected(coefficients, selected, slot.getCijm(j, m, date), "c", j, m);
  1412.                             storeIfSelected(coefficients, selected, slot.getSijm(j, m, date), "s", j, m);
  1413.                         }
  1414.                     }
  1415.                 }

  1417.                 return coefficients;

  1419.             } else {
  1420.                 return Collections.emptyMap();
  1421.             }

  1423.         }

  1425.         /** Put a coefficient in a map if selected.
  1426.          * @param map map to populate
  1427.          * @param selected set of coefficients that should be put in the map
  1428.          * (empty set means all coefficients are selected)
  1429.          * @param value coefficient value
  1430.          * @param id coefficient identifier
  1431.          * @param indices list of coefficient indices
  1432.          */
  1433.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
  1434.                                      final double[] value, final String id, final int... indices) {
  1435.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  1436.             keyBuilder.append(id);
  1437.             for (int index : indices) {
  1438.                 keyBuilder.append('[').append(index).append(']');
  1439.             }
  1440.             final String key = keyBuilder.toString();
  1441.             if (selected.isEmpty() || selected.contains(key)) {
  1442.                 map.put(key, value);
  1443.             }
  1444.         }

  1446.         /** Replace the instance with a data transfer object for serialization.
  1447.          * @return data transfer object that will be serialized
  1448.          * @exception NotSerializableException if an additional state provider is not serializable
  1449.          */
  1450.         private Object writeReplace() throws NotSerializableException {

  1452.             // slots transitions
  1453.             final SortedSet<TimeSpanMap.Transition<Slot>> transitions     = slots.getTransitions();
  1454.             final AbsoluteDate[]                          transitionDates = new AbsoluteDate[transitions.size()];
  1455.             final Slot[]                                  allSlots        = new Slot[transitions.size() + 1];
  1456.             int i = 0;
  1457.             for (final TimeSpanMap.Transition<Slot> transition : transitions) {
  1458.                 if (i == 0) {
  1459.                     // slot before the first transition
  1460.                     allSlots[i] = transition.getBefore();
  1461.                 }
  1462.                 if (i < transitionDates.length) {
  1463.                     transitionDates[i] = transition.getDate();
  1464.                     allSlots[++i]      = transition.getAfter();
  1465.                 }
  1466.             }

  1468.             return new DataTransferObject(bodyFrame, maxOrderMdailyTesseralSP,
  1469.                                           mDailiesOnly, nonResOrders,
  1470.                                           mMax, jMax, interpolationPoints,
  1471.                                           transitionDates, allSlots);

  1473.         }


  1476.         /** Internal class used only for serialization. */
  1477.         private static class DataTransferObject implements Serializable {

  1479.             /** Serializable UID. */
  1480.             private static final long serialVersionUID = 20160319L;

  1482.             /** Central body rotating frame. */
  1483.             private final Frame bodyFrame;

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

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

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

  1494.             /** Maximum value for m index. */
  1495.             private final int mMax;

  1497.             /** Maximum value for j index. */
  1498.             private final int jMax;

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

  1503.             /** Transitions dates. */
  1504.             private final AbsoluteDate[] transitionDates;

  1506.             /** All slots. */
  1507.             private final Slot[] allSlots;

  1509.             /** Simple constructor.
  1510.              * @param bodyFrame central body rotating frame
  1511.              * @param maxOrderMdailyTesseralSP maximal order to consider for short periodics m-daily tesseral harmonics potential
  1512.              * @param mDailiesOnly flag to take into account only M-dailies harmonic tesserals for short periodic perturbations
  1513.              * @param nonResOrders lst of non resonant orders with j != 0
  1514.              * @param mMax maximum value for m index
  1515.              * @param jMax maximum value for j index
  1516.              * @param interpolationPoints number of points used in the interpolation process
  1517.              * @param transitionDates transitions dates
  1518.              * @param allSlots all slots
  1519.              */
  1520.             DataTransferObject(final Frame bodyFrame, final int maxOrderMdailyTesseralSP,
  1521.                                final boolean mDailiesOnly, final SortedMap<Integer, List<Integer> > nonResOrders,
  1522.                                final int mMax, final int jMax, final int interpolationPoints,
  1523.                                final AbsoluteDate[] transitionDates, final Slot[] allSlots) {
  1524.                 this.bodyFrame                = bodyFrame;
  1525.                 this.maxOrderMdailyTesseralSP = maxOrderMdailyTesseralSP;
  1526.                 this.mDailiesOnly             = mDailiesOnly;
  1527.                 this.nonResOrders             = nonResOrders;
  1528.                 this.mMax                     = mMax;
  1529.                 this.jMax                     = jMax;
  1530.                 this.interpolationPoints      = interpolationPoints;
  1531.                 this.transitionDates          = transitionDates;
  1532.                 this.allSlots                 = allSlots;
  1533.             }

  1535.             /** Replace the deserialized data transfer object with a {@link TesseralShortPeriodicCoefficients}.
  1536.              * @return replacement {@link TesseralShortPeriodicCoefficients}
  1537.              */
  1538.             private Object readResolve() {

  1540.                 final TimeSpanMap<Slot> slots = new TimeSpanMap<Slot>(allSlots[0]);
  1541.                 for (int i = 0; i < transitionDates.length; ++i) {
  1542.                     slots.addValidAfter(allSlots[i + 1], transitionDates[i]);
  1543.                 }

  1545.                 return new TesseralShortPeriodicCoefficients(bodyFrame, maxOrderMdailyTesseralSP, mDailiesOnly,
  1546.                                                              nonResOrders, mMax, jMax, interpolationPoints,
  1547.                                                              slots);

  1549.             }

  1551.         }

  1553.     }

  1555.     /** Coefficients valid for one time slot. */
  1556.     private static class Slot implements Serializable {

  1558.         /** Serializable UID. */
  1559.         private static final long serialVersionUID = 20160319L;

  1561.         /** The coefficients C<sub>i</sub><sup>j</sup><sup>m</sup>.
  1562.          * <p>
  1563.          * The index order is cijm[m][j][i] <br/>
  1564.          * i corresponds to the equinoctial element, as follows: <br/>
  1565.          * - i=0 for a <br/>
  1566.          * - i=1 for k <br/>
  1567.          * - i=2 for h <br/>
  1568.          * - i=3 for q <br/>
  1569.          * - i=4 for p <br/>
  1570.          * - i=5 for λ <br/>
  1571.          * </p>
  1572.          */
  1573.         private final ShortPeriodicsInterpolatedCoefficient[][] cijm;

  1575.         /** The coefficients S<sub>i</sub><sup>j</sup><sup>m</sup>.
  1576.          * <p>
  1577.          * The index order is sijm[m][j][i] <br/>
  1578.          * i corresponds to the equinoctial element, as follows: <br/>
  1579.          * - i=0 for a <br/>
  1580.          * - i=1 for k <br/>
  1581.          * - i=2 for h <br/>
  1582.          * - i=3 for q <br/>
  1583.          * - i=4 for p <br/>
  1584.          * - i=5 for λ <br/>
  1585.          * </p>
  1586.          */
  1587.         private final ShortPeriodicsInterpolatedCoefficient[][] sijm;

  1589.         /** Simple constructor.
  1590.          *  @param mMax maximum value for m index
  1591.          *  @param jMax maximum value for j index
  1592.          *  @param interpolationPoints number of points used in the interpolation process
  1593.          */
  1594.         Slot(final int mMax, final int jMax, final int interpolationPoints) {

  1596.             final int rows    = mMax + 1;
  1597.             final int columns = 2 * jMax + 1;
  1598.             cijm = new ShortPeriodicsInterpolatedCoefficient[rows][columns];
  1599.             sijm = new ShortPeriodicsInterpolatedCoefficient[rows][columns];
  1600.             for (int m = 1; m <= mMax; m++) {
  1601.                 for (int j = -jMax; j <= jMax; j++) {
  1602.                     cijm[m][j + jMax] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  1603.                     sijm[m][j + jMax] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  1604.                 }
  1605.             }

  1607.         }

  1609.         /** Get C<sub>i</sub><sup>j</sup><sup>m</sup>.
  1610.          *
  1611.          * @param j j index
  1612.          * @param m m index
  1613.          * @param date the date
  1614.          * @return C<sub>i</sub><sup>j</sup><sup>m</sup>
  1615.          */
  1616.         double[] getCijm(final int j, final int m, final AbsoluteDate date) {
  1617.             final int jMax = (cijm[m].length - 1) / 2;
  1618.             return cijm[m][j + jMax].value(date);
  1619.         }

  1621.         /** Get S<sub>i</sub><sup>j</sup><sup>m</sup>.
  1622.          *
  1623.          * @param j j index
  1624.          * @param m m index
  1625.          * @param date the date
  1626.          * @return S<sub>i</sub><sup>j</sup><sup>m</sup>
  1627.          */
  1628.         double[] getSijm(final int j, final int m, final AbsoluteDate date) {
  1629.             final int jMax = (cijm[m].length - 1) / 2;
  1630.             return sijm[m][j + jMax].value(date);
  1631.         }

  1633.     }

  1635. }
Created with JaCoCo 0.7.9.201702052155