DSSTZonal.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.HashMap;
  23. import java.util.List;
  24. import java.util.Map;
  25. import java.util.Set;
  26. import java.util.SortedSet;
  27. import java.util.TreeMap;

  29. import org.hipparchus.exception.LocalizedCoreFormats;
  30. import org.hipparchus.util.FastMath;
  31. import org.orekit.attitudes.AttitudeProvider;
  32. import org.orekit.errors.OrekitException;
  33. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  34. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
  35. import org.orekit.orbits.Orbit;
  36. import org.orekit.propagation.SpacecraftState;
  37. import org.orekit.propagation.events.EventDetector;
  38. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  39. import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
  40. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
  41. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory.NSKey;
  42. import org.orekit.propagation.semianalytical.dsst.utilities.GHIJjsPolynomials;
  43. import org.orekit.propagation.semianalytical.dsst.utilities.LnsCoefficients;
  44. import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
  45. import org.orekit.propagation.semianalytical.dsst.utilities.UpperBounds;
  46. import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenZonalLinear;
  47. import org.orekit.time.AbsoluteDate;
  48. import org.orekit.utils.TimeSpanMap;

  50. /** Zonal contribution to the central body gravitational perturbation.
  51.  *
  52.  *   @author Romain Di Costanzo
  53.  *   @author Pascal Parraud
  54.  */
  55. public class DSSTZonal implements DSSTForceModel {

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

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

  63.     /** Retrograde factor I.
  64.      *  <p>
  65.      *  DSST model needs equinoctial orbit as internal representation.
  66.      *  Classical equinoctial elements have discontinuities when inclination
  67.      *  is close to zero. In this representation, I = +1. <br>
  68.      *  To avoid this discontinuity, another representation exists and equinoctial
  69.      *  elements can be expressed in a different way, called "retrograde" orbit.
  70.      *  This implies I = -1. <br>
  71.      *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  72.      *  But for the sake of consistency with the theory, the retrograde factor
  73.      *  has been kept in the formulas.
  74.      *  </p>
  75.      */
  76.     private static final int I = 1;

  78.     /** Provider for spherical harmonics. */
  79.     private final UnnormalizedSphericalHarmonicsProvider provider;

  81.     /** Maximal degree to consider for harmonics potential. */
  82.     private final int maxDegree;

  84.     /** Maximal degree to consider for harmonics potential. */
  85.     private final int maxDegreeShortPeriodics;

  87.     /** Maximal degree to consider for harmonics potential in short periodic computations. */
  88.     private final int maxOrder;

  90.     /** Factorial. */
  91.     private final double[] fact;

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

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

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

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

  105.     /** Short period terms. */
  106.     private ZonalShortPeriodicCoefficients zonalSPCoefs;

  108.     // Equinoctial elements (according to DSST notation)
  109.     /** a. */
  110.     private double a;
  111.     /** ex. */
  112.     private double k;
  113.     /** ey. */
  114.     private double h;
  115.     /** hx. */
  116.     private double q;
  117.     /** hy. */
  118.     private double p;

  120.     /** Eccentricity. */
  121.     private double ecc;

  123.     /** Direction cosine &alpha. */
  124.     private double alpha;
  125.     /** Direction cosine &beta. */
  126.     private double beta;
  127.     /** Direction cosine &gamma. */
  128.     private double gamma;

  130.     // Common factors for potential computation
  131.     /** &Chi; = 1 / sqrt(1 - e²) = 1 / B. */
  132.     private double X;
  133.     /** &Chi;². */
  134.     private double XX;
  135.     /** &Chi;³. */
  136.     private double XXX;
  137.     /** 1 / (A * B) .*/
  138.     private double ooAB;
  139.     /** B / A .*/
  140.     private double BoA;
  141.     /** B / A(1 + B) .*/
  142.     private double BoABpo;
  143.     /** -C / (2 * A * B) .*/
  144.     private double mCo2AB;
  145.     /** -2 * a / A .*/
  146.     private double m2aoA;
  147.     /** μ / a .*/
  148.     private double muoa;
  149.     /** R / a .*/
  150.     private double roa;

  152.     /** An array that contains the objects needed to build the Hansen coefficients. <br/>
  153.      * The index is s*/
  154.     private HansenZonalLinear[] hansenObjects;

  156.     /** The current value of the U function. <br/>
  157.      * Needed when computed the short periodic contribution */
  158.     private double U;

  160.     /** A = sqrt( μ * a ) = n * a². */
  161.     private double A;
  162.     /** B = sqrt( 1 - k² - h² ). */
  163.     private double B;
  164.     /** C = 1 + p² + Q². */
  165.     private double C;

  167.     /** The mean motion (n). */
  168.     private double meanMotion;

  170.     /** h * k. */
  171.     private double hk;
  172.     /** k² - h². */
  173.     private double k2mh2;
  174.     /** (k² - h²) / 2. */
  175.     private double k2mh2o2;
  176.     /** 1 / (n² * a²). */
  177.     private double oon2a2;
  178.     /** 1 / (n² * a) . */
  179.     private double oon2a;
  180.     /** χ³ / (n² * a). */
  181.     private double x3on2a;
  182.     /** χ / (n² * a²). */
  183.     private double xon2a2;
  184.     /** (C * χ) / ( 2 * n² * a² ). */
  185.     private double cxo2n2a2;
  186.     /** (χ²) / (n² * a² * (χ + 1 ) ). */
  187.     private double x2on2a2xp1;
  188.     /** B * B.*/
  189.     private double BB;

  191.     /** Simple constructor.
  192.      * @param provider provider for spherical harmonics
  193.      * @param maxDegreeShortPeriodics maximum degree to consider for short periodics zonal harmonics potential
  194.      * (must be between 2 and {@code provider.getMaxDegree()})
  195.      * @param maxEccPowShortPeriodics maximum power of the eccentricity to be used in short periodic computations
  196.      * (must be between 0 and {@code maxDegreeShortPeriodics - 1}, but should typically not exceed 4 as higher
  197.      * values will exceed computer capacity)
  198.      * @param maxFrequencyShortPeriodics maximum frequency in true longitude for short periodic computations
  199.      * (must be between 1 and {@code 2 * maxDegreeShortPeriodics + 1})
  200.      * @exception OrekitException if degrees or powers are out of range
  201.      * @since 7.2
  202.      */
  203.     public DSSTZonal(final UnnormalizedSphericalHarmonicsProvider provider,
  204.                      final int maxDegreeShortPeriodics,
  205.                      final int maxEccPowShortPeriodics,
  206.                      final int maxFrequencyShortPeriodics)
  207.         throws OrekitException {

  209.         this.provider  = provider;
  210.         this.maxDegree = provider.getMaxDegree();
  211.         this.maxOrder  = provider.getMaxOrder();

  213.         checkIndexRange(maxDegreeShortPeriodics, 2, maxDegree);
  214.         this.maxDegreeShortPeriodics = maxDegreeShortPeriodics;

  216.         checkIndexRange(maxEccPowShortPeriodics, 0, maxDegreeShortPeriodics - 1);
  217.         this.maxEccPowShortPeriodics = maxEccPowShortPeriodics;

  219.         checkIndexRange(maxFrequencyShortPeriodics, 1, 2 * maxDegreeShortPeriodics + 1);
  220.         this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;

  222.         // Vns coefficients
  223.         this.Vns = CoefficientsFactory.computeVns(maxDegree + 1);

  225.         // Factorials computation
  226.         final int maxFact = 2 * maxDegree + 1;
  227.         this.fact = new double[maxFact];
  228.         fact[0] = 1.;
  229.         for (int i = 1; i < maxFact; i++) {
  230.             fact[i] = i * fact[i - 1];
  231.         }

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

  236.     }

  238.     /** Check an index range.
  239.      * @param index index value
  240.      * @param min minimum value for index
  241.      * @param max maximum value for index
  242.      * @exception OrekitException if index is out of range
  243.      */
  244.     private void checkIndexRange(final int index, final int min, final int max)
  245.         throws OrekitException {
  246.         if (index < min || index > max) {
  247.             throw new OrekitException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, index, min, max);
  248.         }
  249.     }

  251.     /** Get the spherical harmonics provider.
  252.      *  @return the spherical harmonics provider
  253.      */
  254.     public UnnormalizedSphericalHarmonicsProvider getProvider() {
  255.         return provider;
  256.     }

  258.     /** {@inheritDoc}
  259.      *  <p>
  260.      *  Computes the highest power of the eccentricity to appear in the truncated
  261.      *  analytical power series expansion.
  262.      *  </p>
  263.      *  <p>
  264.      *  This method computes the upper value for the central body potential and
  265.      *  determines the maximal power for the eccentricity producing potential
  266.      *  terms bigger than a defined tolerance.
  267.      *  </p>
  268.      */
  269.     @Override
  270.     public List<ShortPeriodTerms> initialize(final AuxiliaryElements aux, final boolean meanOnly)
  271.         throws OrekitException {

  273.         computeMeanElementsTruncations(aux);

  275.         final int maxEccPow;
  276.         if (!meanOnly) {
  277.             maxEccPow = FastMath.max(maxEccPowMeanElements, maxEccPowShortPeriodics);
  278.         } else {
  279.             maxEccPow = maxEccPowMeanElements;
  280.         }

  282.         //Initialize the HansenCoefficient generator
  283.         this.hansenObjects = new HansenZonalLinear[maxEccPow + 1];

  285.         for (int s = 0; s <= maxEccPow; s++) {
  286.             this.hansenObjects[s] = new HansenZonalLinear(maxDegree, s);
  287.         }

  289.         final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
  290.         zonalSPCoefs = new ZonalShortPeriodicCoefficients(maxFrequencyShortPeriodics,
  291.                                                           INTERPOLATION_POINTS,
  292.                                                           new TimeSpanMap<Slot>(new Slot(maxFrequencyShortPeriodics,
  293.                                                                                          INTERPOLATION_POINTS)));
  294.         list.add(zonalSPCoefs);
  295.         return list;

  297.     }

  299.     /** Compute indices truncations for mean elements computations.
  300.      * @param aux auxiliary elements
  301.      * @throws OrekitException if an error occurs
  302.      */
  303.     private void computeMeanElementsTruncations(final AuxiliaryElements aux) throws OrekitException {

  305.         //Compute the max eccentricity power for the mean element rate expansion
  306.         if (maxDegree == 2) {
  307.             maxEccPowMeanElements = 0;
  308.         } else {
  309.             // Initializes specific parameters.
  310.             initializeStep(aux);
  311.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(aux.getDate());

  313.             // Utilities for truncation
  314.             final double ax2or = 2. * a / provider.getAe();
  315.             double xmuran = provider.getMu() / a;
  316.             // Set a lower bound for eccentricity
  317.             final double eo2  = FastMath.max(0.0025, ecc / 2.);
  318.             final double x2o2 = XX / 2.;
  319.             final double[] eccPwr = new double[maxDegree + 1];
  320.             final double[] chiPwr = new double[maxDegree + 1];
  321.             final double[] hafPwr = new double[maxDegree + 1];
  322.             eccPwr[0] = 1.;
  323.             chiPwr[0] = X;
  324.             hafPwr[0] = 1.;
  325.             for (int i = 1; i <= maxDegree; i++) {
  326.                 eccPwr[i] = eccPwr[i - 1] * eo2;
  327.                 chiPwr[i] = chiPwr[i - 1] * x2o2;
  328.                 hafPwr[i] = hafPwr[i - 1] * 0.5;
  329.                 xmuran  /= ax2or;
  330.             }

  332.             // Set highest power of e and degree of current spherical harmonic.
  333.             maxEccPowMeanElements = 0;
  334.             int n = maxDegree;
  335.             // Loop over n
  336.             do {
  337.                 // Set order of current spherical harmonic.
  338.                 int m = 0;
  339.                 // Loop over m
  340.                 do {
  341.                     // Compute magnitude of current spherical harmonic coefficient.
  342.                     final double cnm = harmonics.getUnnormalizedCnm(n, m);
  343.                     final double snm = harmonics.getUnnormalizedSnm(n, m);
  344.                     final double csnm = FastMath.hypot(cnm, snm);
  345.                     if (csnm == 0.) break;
  346.                     // Set magnitude of last spherical harmonic term.
  347.                     double lastTerm = 0.;
  348.                     // Set current power of e and related indices.
  349.                     int nsld2 = (n - maxEccPowMeanElements - 1) / 2;
  350.                     int l = n - 2 * nsld2;
  351.                     // Loop over l
  352.                     double term = 0.;
  353.                     do {
  354.                         // Compute magnitude of current spherical harmonic term.
  355.                         if (m < l) {
  356.                             term = csnm * xmuran *
  357.                                     (fact[n - l] / (fact[n - m])) *
  358.                                     (fact[n + l] / (fact[nsld2] * fact[nsld2 + l])) *
  359.                                     eccPwr[l] * UpperBounds.getDnl(XX, chiPwr[l], n, l) *
  360.                                     (UpperBounds.getRnml(gamma, n, l, m, 1, I) + UpperBounds.getRnml(gamma, n, l, m, -1, I));
  361.                         } else {
  362.                             term = csnm * xmuran *
  363.                                     (fact[n + m] / (fact[nsld2] * fact[nsld2 + l])) *
  364.                                     eccPwr[l] * hafPwr[m - l] * UpperBounds.getDnl(XX, chiPwr[l], n, l) *
  365.                                     (UpperBounds.getRnml(gamma, n, m, l, 1, I) + UpperBounds.getRnml(gamma, n, m, l, -1, I));
  366.                         }
  367.                         // Is the current spherical harmonic term bigger than the truncation tolerance ?
  368.                         if (term >= TRUNCATION_TOLERANCE) {
  369.                             maxEccPowMeanElements = l;
  370.                         } else {
  371.                             // Is the current term smaller than the last term ?
  372.                             if (term < lastTerm) {
  373.                                 break;
  374.                             }
  375.                         }
  376.                         // Proceed to next power of e.
  377.                         lastTerm = term;
  378.                         l += 2;
  379.                         nsld2--;
  380.                     } while (l < n);
  381.                     // Is the current spherical harmonic term bigger than the truncation tolerance ?
  382.                     if (term >= TRUNCATION_TOLERANCE) {
  383.                         maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
  384.                         return;
  385.                     }
  386.                     // Proceed to next order.
  387.                     m++;
  388.                 } while (m <= FastMath.min(n, maxOrder));
  389.                 // Proceed to next degree.
  390.                 xmuran *= ax2or;
  391.                 n--;
  392.             } while (n > maxEccPowMeanElements + 2);

  394.             maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
  395.         }
  396.     }

  398.     /** {@inheritDoc} */
  399.     @Override
  400.     public void initializeStep(final AuxiliaryElements aux) throws OrekitException {

  402.         // Equinoctial elements
  403.         a = aux.getSma();
  404.         k = aux.getK();
  405.         h = aux.getH();
  406.         q = aux.getQ();
  407.         p = aux.getP();

  409.         // Eccentricity
  410.         ecc = aux.getEcc();

  412.         // Direction cosines
  413.         alpha = aux.getAlpha();
  414.         beta  = aux.getBeta();
  415.         gamma = aux.getGamma();

  417.         // Equinoctial coefficients
  418.         A = aux.getA();
  419.         B = aux.getB();
  420.         C = aux.getC();

  422.         // &Chi; = 1 / B
  423.         X = 1. / B;
  424.         XX = X * X;
  425.         XXX = X * XX;
  426.         // 1 / AB
  427.         ooAB = 1. / (A * B);
  428.         // B / A
  429.         BoA = B / A;
  430.         // -C / 2AB
  431.         mCo2AB = -C * ooAB / 2.;
  432.         // B / A(1 + B)
  433.         BoABpo = BoA / (1. + B);
  434.         // -2 * a / A
  435.         m2aoA = -2 * a / A;
  436.         // μ / a
  437.         muoa = provider.getMu() / a;
  438.         // R / a
  439.         roa = provider.getAe() / a;

  441.         // Mean motion
  442.         meanMotion = aux.getMeanMotion();
  443.     }

  445.     /** {@inheritDoc} */
  446.     @Override
  447.     public double[] getMeanElementRate(final SpacecraftState spacecraftState) throws OrekitException {
  448.         return computeMeanElementRates(spacecraftState.getDate());
  449.     }

  451.     /** {@inheritDoc} */
  452.     @Override
  453.     public EventDetector[] getEventsDetectors() {
  454.         return null;
  455.     }

  457.     /** Compute the mean element rates.
  458.      * @param date current date
  459.      * @return the mean element rates
  460.      * @throws OrekitException if an error occurs in hansen computation
  461.      */
  462.     private double[] computeMeanElementRates(final AbsoluteDate date) throws OrekitException {
  463.         // Compute potential derivative
  464.         final double[] dU  = computeUDerivatives(date);
  465.         final double dUda  = dU[0];
  466.         final double dUdk  = dU[1];
  467.         final double dUdh  = dU[2];
  468.         final double dUdAl = dU[3];
  469.         final double dUdBe = dU[4];
  470.         final double dUdGa = dU[5];

  472.         // Compute cross derivatives [Eq. 2.2-(8)]
  473.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  474.         final double UAlphaGamma   = alpha * dUdGa - gamma * dUdAl;
  475.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  476.         final double UBetaGamma    =  beta * dUdGa - gamma * dUdBe;
  477.         // Common factor
  478.         final double pUAGmIqUBGoAB = (p * UAlphaGamma - I * q * UBetaGamma) * ooAB;

  480.         // Compute mean elements rates [Eq. 3.1-(1)]
  481.         final double da =  0.;
  482.         final double dh =  BoA * dUdk + k * pUAGmIqUBGoAB;
  483.         final double dk = -BoA * dUdh - h * pUAGmIqUBGoAB;
  484.         final double dp =  mCo2AB * UBetaGamma;
  485.         final double dq =  mCo2AB * UAlphaGamma * I;
  486.         final double dM =  m2aoA * dUda + BoABpo * (h * dUdh + k * dUdk) + pUAGmIqUBGoAB;

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

  491.     /** Compute the derivatives of the gravitational potential U [Eq. 3.1-(6)].
  492.      *  <p>
  493.      *  The result is the array
  494.      *  [dU/da, dU/dk, dU/dh, dU/dα, dU/dβ, dU/dγ]
  495.      *  </p>
  496.      *  @param date current date
  497.      *  @return potential derivatives
  498.      *  @throws OrekitException if an error occurs in hansen computation
  499.      */
  500.     private double[] computeUDerivatives(final AbsoluteDate date) throws OrekitException {

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

  504.         //Reset U
  505.         U = 0.;

  507.         // Gs and Hs coefficients
  508.         final double[][] GsHs = CoefficientsFactory.computeGsHs(k, h, alpha, beta, maxEccPowMeanElements);
  509.         // Qns coefficients
  510.         final double[][] Qns  = CoefficientsFactory.computeQns(gamma, maxDegree, maxEccPowMeanElements);

  512.         final double[] roaPow = new double[maxDegree + 1];
  513.         roaPow[0] = 1.;
  514.         for (int i = 1; i <= maxDegree; i++) {
  515.             roaPow[i] = roa * roaPow[i - 1];
  516.         }

  518.         // Potential derivatives
  519.         double dUda  = 0.;
  520.         double dUdk  = 0.;
  521.         double dUdh  = 0.;
  522.         double dUdAl = 0.;
  523.         double dUdBe = 0.;
  524.         double dUdGa = 0.;

  526.         for (int s = 0; s <= maxEccPowMeanElements; s++) {
  527.             //Initialize the Hansen roots
  528.             this.hansenObjects[s].computeInitValues(X);

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

  533.             // Compute Gs partial derivatives from 3.1-(9)
  534.             double dGsdh  = 0.;
  535.             double dGsdk  = 0.;
  536.             double dGsdAl = 0.;
  537.             double dGsdBe = 0.;
  538.             if (s > 0) {
  539.                 // First get the G(s-1) and the H(s-1) coefficients
  540.                 final double sxgsm1 = s * GsHs[0][s - 1];
  541.                 final double sxhsm1 = s * GsHs[1][s - 1];
  542.                 // Then compute derivatives
  543.                 dGsdh  = beta  * sxgsm1 - alpha * sxhsm1;
  544.                 dGsdk  = alpha * sxgsm1 + beta  * sxhsm1;
  545.                 dGsdAl = k * sxgsm1 - h * sxhsm1;
  546.                 dGsdBe = h * sxgsm1 + k * sxhsm1;
  547.             }

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

  552.             for (int n = s + 2; n <= maxDegree; n++) {
  553.                 // (n - s) must be even
  554.                 if ((n - s) % 2 == 0) {

  556.                     //Extract data from previous computation :
  557.                     final double kns   = this.hansenObjects[s].getValue(-n - 1, X);
  558.                     final double dkns  = this.hansenObjects[s].getDerivative(-n - 1, X);

  560.                     final double vns   = Vns.get(new NSKey(n, s));
  561.                     final double coef0 = d0s * roaPow[n] * vns * -harmonics.getUnnormalizedCnm(n, 0);
  562.                     final double coef1 = coef0 * Qns[n][s];
  563.                     final double coef2 = coef1 * kns;
  564.                     final double coef3 = coef2 * gs;
  565.                     // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  566.                     final double dqns  = Qns[n][s + 1];

  568.                     // Compute U
  569.                     U += coef3;
  570.                     // Compute dU / da :
  571.                     dUda += coef3 * (n + 1);
  572.                     // Compute dU / dEx
  573.                     dUdk += coef1 * (kns * dGsdk + k * XXX * gs * dkns);
  574.                     // Compute dU / dEy
  575.                     dUdh += coef1 * (kns * dGsdh + h * XXX * gs * dkns);
  576.                     // Compute dU / dAlpha
  577.                     dUdAl += coef2 * dGsdAl;
  578.                     // Compute dU / dBeta
  579.                     dUdBe += coef2 * dGsdBe;
  580.                     // Compute dU / dGamma
  581.                     dUdGa += coef0 * kns * dqns * gs;

  583.                 }
  584.             }
  585.         }

  587.         // Multiply by -(μ / a)
  588.         U *= -muoa;

  590.         return new double[] {
  591.             dUda  *  muoa / a,
  592.             dUdk  * -muoa,
  593.             dUdh  * -muoa,
  594.             dUdAl * -muoa,
  595.             dUdBe * -muoa,
  596.             dUdGa * -muoa
  597.         };
  598.     }

  600.     /** {@inheritDoc} */
  601.     @Override
  602.     public void registerAttitudeProvider(final AttitudeProvider attitudeProvider) {
  603.         //nothing is done since this contribution is not sensitive to attitude
  604.     }

  606.     /** Check if an index is within the accepted interval.
  607.      *
  608.      * @param index the index to check
  609.      * @param lowerBound the lower bound of the interval
  610.      * @param upperBound the upper bound of the interval
  611.      * @return true if the index is between the lower and upper bounds, false otherwise
  612.      */
  613.     private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
  614.         return index >= lowerBound && index <= upperBound;
  615.     }

  617.     /** {@inheritDoc} */
  618.     @Override
  619.     public void updateShortPeriodTerms(final SpacecraftState... meanStates)
  620.         throws OrekitException {

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

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

  627.             // h * k.
  628.             this.hk = h * k;
  629.             // k² - h².
  630.             this.k2mh2 = k * k - h * h;
  631.             // (k² - h²) / 2.
  632.             this.k2mh2o2 = k2mh2 / 2.;
  633.             // 1 / (n² * a²) = 1 / (n * A)
  634.             this.oon2a2 = 1 / (A * meanMotion);
  635.             // 1 / (n² * a) = a / (n * A)
  636.             this.oon2a = a * oon2a2;
  637.             // χ³ / (n² * a)
  638.             this.x3on2a = XXX * oon2a;
  639.             // χ / (n² * a²)
  640.             this.xon2a2 = X * oon2a2;
  641.             // (C * χ) / ( 2 * n² * a² )
  642.             this.cxo2n2a2 = xon2a2 * C / 2;
  643.             // (χ²) / (n² * a² * (χ + 1 ) )
  644.             this.x2on2a2xp1 = xon2a2 * X / (X + 1);
  645.             // B * B
  646.             this.BB = B * B;

  648.             // Compute rhoj and sigmaj
  649.             final double[][] rhoSigma = computeRhoSigmaCoefficients(meanState.getDate(), slot);

  651.             // Compute Di
  652.             computeDiCoefficients(meanState.getDate(), slot);

  654.             // generate the Cij and Sij coefficients
  655.             final FourierCjSjCoefficients cjsj = new FourierCjSjCoefficients(meanState.getDate(),
  656.                                                                              maxDegreeShortPeriodics,
  657.                                                                              maxEccPowShortPeriodics,
  658.                                                                              maxFrequencyShortPeriodics);
  659.             computeCijSijCoefficients(meanState.getDate(), slot, cjsj, rhoSigma);
  660.         }

  662.     }

  664.     /** Generate the values for the D<sub>i</sub> coefficients.
  665.      * @param date target date
  666.      * @param slot slot to which the coefficients belong
  667.      * @throws OrekitException if an error occurs during the coefficient computation
  668.      */
  669.     private void computeDiCoefficients(final AbsoluteDate date, final Slot slot)
  670.         throws OrekitException {
  671.         final double[] meanElementRates = computeMeanElementRates(date);
  672.         final double[] currentDi = new double[6];

  674.         // Add the coefficients to the interpolation grid
  675.         for (int i = 0; i < 6; i++) {
  676.             currentDi[i] = meanElementRates[i] / meanMotion;

  678.             if (i == 5) {
  679.                 currentDi[i] += -1.5 * 2 * U * oon2a2;
  680.             }

  682.         }

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

  686.     }

  688.     /**
  689.      * Generate the values for the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
  690.      * @param date date of computation
  691.      * @param slot slot to which the coefficients belong
  692.      * @param cjsj Fourier coefficients
  693.      * @param rhoSigma ρ<sub>j</sub> and σ<sub>j</sub>
  694.      */
  695.     private void computeCijSijCoefficients(final AbsoluteDate date, final Slot slot,
  696.                                            final FourierCjSjCoefficients cjsj,
  697.                                            final double[][] rhoSigma) {

  699.         final int nMax = maxDegreeShortPeriodics;

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

  705.             // Create local arrays
  706.             final double[] currentCij = new double[] {0., 0., 0., 0., 0., 0.};
  707.             final double[] currentSij = new double[] {0., 0., 0., 0., 0., 0.};

  709.             // j == 1
  710.             if (j == 1) {
  711.                 final double coef1 = 4 * k * U - hk * cjsj.getCj(1) + k2mh2o2 * cjsj.getSj(1);
  712.                 final double coef2 = 4 * h * U + k2mh2o2 * cjsj.getCj(1) + hk * cjsj.getSj(1);
  713.                 final double coef3 = (k * cjsj.getCj(1) + h * cjsj.getSj(1)) / 4.;
  714.                 final double coef4 = (8 * U - h * cjsj.getCj(1) + k * cjsj.getSj(1)) / 4.;

  716.                 //Coefficients for a
  717.                 currentCij[0] += coef1;
  718.                 currentSij[0] += coef2;

  720.                 //Coefficients for k
  721.                 currentCij[1] += coef4;
  722.                 currentSij[1] += coef3;

  724.                 //Coefficients for h
  725.                 currentCij[2] -= coef3;
  726.                 currentSij[2] += coef4;

  728.                 //Coefficients for λ
  729.                 currentCij[5] -= coef2 / 2;
  730.                 currentSij[5] += coef1 / 2;
  731.             }

  733.             // j == 2
  734.             if (j == 2) {
  735.                 final double coef1 = k2mh2 * U;
  736.                 final double coef2 = 2 * hk * U;
  737.                 final double coef3 = h * U / 2;
  738.                 final double coef4 = k * U / 2;

  740.                 //Coefficients for a
  741.                 currentCij[0] += coef1;
  742.                 currentSij[0] += coef2;

  744.                 //Coefficients for k
  745.                 currentCij[1] += coef4;
  746.                 currentSij[1] += coef3;

  748.                 //Coefficients for h
  749.                 currentCij[2] -= coef3;
  750.                 currentSij[2] += coef4;

  752.                 //Coefficients for λ
  753.                 currentCij[5] -= coef2 / 2;
  754.                 currentSij[5] += coef1 / 2;
  755.             }

  757.             // j between 1 and 2N-3
  758.             if (isBetween(j, 1, 2 * nMax - 3) && j + 2 < cjsj.jMax) {
  759.                 final double coef1 = ( j + 2 ) * (-hk * cjsj.getCj(j + 2) + k2mh2o2 * cjsj.getSj(j + 2));
  760.                 final double coef2 = ( j + 2 ) * (k2mh2o2 * cjsj.getCj(j + 2) + hk * cjsj.getSj(j + 2));
  761.                 final double coef3 = ( j + 2 ) * (k * cjsj.getCj(j + 2) + h * cjsj.getSj(j + 2)) / 4;
  762.                 final double coef4 = ( j + 2 ) * (h * cjsj.getCj(j + 2) - k * cjsj.getSj(j + 2)) / 4;

  764.                 //Coefficients for a
  765.                 currentCij[0] += coef1;
  766.                 currentSij[0] -= coef2;

  768.                 //Coefficients for k
  769.                 currentCij[1] += -coef4;
  770.                 currentSij[1] -= coef3;

  772.                 //Coefficients for h
  773.                 currentCij[2] -= coef3;
  774.                 currentSij[2] += coef4;

  776.                 //Coefficients for λ
  777.                 currentCij[5] -= coef2 / 2;
  778.                 currentSij[5] += -coef1 / 2;
  779.             }

  781.             // j between 1 and 2N-2
  782.             if (isBetween(j, 1, 2 * nMax - 2) && j + 1 < cjsj.jMax) {
  783.                 final double coef1 = 2 * ( j + 1 ) * (-h * cjsj.getCj(j + 1) + k * cjsj.getSj(j + 1));
  784.                 final double coef2 = 2 * ( j + 1 ) * (k * cjsj.getCj(j + 1) + h * cjsj.getSj(j + 1));
  785.                 final double coef3 = ( j + 1 ) * cjsj.getCj(j + 1);
  786.                 final double coef4 = ( j + 1 ) * cjsj.getSj(j + 1);

  788.                 //Coefficients for a
  789.                 currentCij[0] += coef1;
  790.                 currentSij[0] -= coef2;

  792.                 //Coefficients for k
  793.                 currentCij[1] += coef4;
  794.                 currentSij[1] -= coef3;

  796.                 //Coefficients for h
  797.                 currentCij[2] -= coef3;
  798.                 currentSij[2] -= coef4;

  800.                 //Coefficients for λ
  801.                 currentCij[5] -= coef2 / 2;
  802.                 currentSij[5] += -coef1 / 2;
  803.             }

  805.             // j between 2 and 2N
  806.             if (isBetween(j, 2, 2 * nMax) && j - 1 < cjsj.jMax) {
  807.                 final double coef1 = 2 * ( j - 1 ) * (h * cjsj.getCj(j - 1) + k * cjsj.getSj(j - 1));
  808.                 final double coef2 = 2 * ( j - 1 ) * (k * cjsj.getCj(j - 1) - h * cjsj.getSj(j - 1));
  809.                 final double coef3 = ( j - 1 ) * cjsj.getCj(j - 1);
  810.                 final double coef4 = ( j - 1 ) * cjsj.getSj(j - 1);

  812.                 //Coefficients for a
  813.                 currentCij[0] += coef1;
  814.                 currentSij[0] -= coef2;

  816.                 //Coefficients for k
  817.                 currentCij[1] += coef4;
  818.                 currentSij[1] -= coef3;

  820.                 //Coefficients for h
  821.                 currentCij[2] += coef3;
  822.                 currentSij[2] += coef4;

  824.                 //Coefficients for λ
  825.                 currentCij[5] += coef2 / 2;
  826.                 currentSij[5] += coef1 / 2;
  827.             }

  829.             // j between 3 and 2N + 1
  830.             if (isBetween(j, 3, 2 * nMax + 1) && j - 2 < cjsj.jMax) {
  831.                 final double coef1 = ( j - 2 ) * (hk * cjsj.getCj(j - 2) + k2mh2o2 * cjsj.getSj(j - 2));
  832.                 final double coef2 = ( j - 2 ) * (-k2mh2o2 * cjsj.getCj(j - 2) + hk * cjsj.getSj(j - 2));
  833.                 final double coef3 = ( j - 2 ) * (k * cjsj.getCj(j - 2) - h * cjsj.getSj(j - 2)) / 4;
  834.                 final double coef4 = ( j - 2 ) * (h * cjsj.getCj(j - 2) + k * cjsj.getSj(j - 2)) / 4;
  835.                 final double coef5 = ( j - 2 ) * (k2mh2o2 * cjsj.getCj(j - 2) - hk * cjsj.getSj(j - 2));

  837.                 //Coefficients for a
  838.                 currentCij[0] += coef1;
  839.                 currentSij[0] += coef2;

  841.                 //Coefficients for k
  842.                 currentCij[1] += coef4;
  843.                 currentSij[1] += -coef3;

  845.                 //Coefficients for h
  846.                 currentCij[2] += coef3;
  847.                 currentSij[2] += coef4;

  849.                 //Coefficients for λ
  850.                 currentCij[5] += coef5 / 2;
  851.                 currentSij[5] += coef1 / 2;
  852.             }

  854.             //multiply by the common factor
  855.             //for a (i == 0) -> χ³ / (n² * a)
  856.             currentCij[0] *= this.x3on2a;
  857.             currentSij[0] *= this.x3on2a;
  858.             //for k (i == 1) -> χ / (n² * a²)
  859.             currentCij[1] *= this.xon2a2;
  860.             currentSij[1] *= this.xon2a2;
  861.             //for h (i == 2) -> χ / (n² * a²)
  862.             currentCij[2] *= this.xon2a2;
  863.             currentSij[2] *= this.xon2a2;
  864.             //for λ (i == 5) -> (χ²) / (n² * a² * (χ + 1 ) )
  865.             currentCij[5] *= this.x2on2a2xp1;
  866.             currentSij[5] *= this.x2on2a2xp1;

  868.             // j is between 1 and 2 * N - 1
  869.             if (isBetween(j, 1, 2 * nMax - 1) && j < cjsj.jMax) {
  870.                 // Compute cross derivatives
  871.                 // Cj(alpha,gamma) = alpha * dC/dgamma - gamma * dC/dalpha
  872.                 final double CjAlphaGamma   = alpha * cjsj.getdCjdGamma(j) - gamma * cjsj.getdCjdAlpha(j);
  873.                 // Cj(alpha,beta) = alpha * dC/dbeta - beta * dC/dalpha
  874.                 final double CjAlphaBeta   = alpha * cjsj.getdCjdBeta(j) - beta * cjsj.getdCjdAlpha(j);
  875.                 // Cj(beta,gamma) = beta * dC/dgamma - gamma * dC/dbeta
  876.                 final double CjBetaGamma    =  beta * cjsj.getdCjdGamma(j) - gamma * cjsj.getdCjdBeta(j);
  877.                 // Cj(h,k) = h * dC/dk - k * dC/dh
  878.                 final double CjHK   = h * cjsj.getdCjdK(j) - k * cjsj.getdCjdH(j);
  879.                 // Sj(alpha,gamma) = alpha * dS/dgamma - gamma * dS/dalpha
  880.                 final double SjAlphaGamma   = alpha * cjsj.getdSjdGamma(j) - gamma * cjsj.getdSjdAlpha(j);
  881.                 // Sj(alpha,beta) = alpha * dS/dbeta - beta * dS/dalpha
  882.                 final double SjAlphaBeta   = alpha * cjsj.getdSjdBeta(j) - beta * cjsj.getdSjdAlpha(j);
  883.                 // Sj(beta,gamma) = beta * dS/dgamma - gamma * dS/dbeta
  884.                 final double SjBetaGamma    =  beta * cjsj.getdSjdGamma(j) - gamma * cjsj.getdSjdBeta(j);
  885.                 // Sj(h,k) = h * dS/dk - k * dS/dh
  886.                 final double SjHK   = h * cjsj.getdSjdK(j) - k * cjsj.getdSjdH(j);

  888.                 //Coefficients for a
  889.                 final double coef1 = this.x3on2a * (3 - BB) * j;
  890.                 currentCij[0] += coef1 * cjsj.getSj(j);
  891.                 currentSij[0] -= coef1 * cjsj.getCj(j);

  893.                 //Coefficients for k and h
  894.                 final double coef2 = p * CjAlphaGamma - I * q * CjBetaGamma;
  895.                 final double coef3 = p * SjAlphaGamma - I * q * SjBetaGamma;
  896.                 currentCij[1] -= this.xon2a2 * (h * coef2 + BB * cjsj.getdCjdH(j) - 1.5 * k * j * cjsj.getSj(j));
  897.                 currentSij[1] -= this.xon2a2 * (h * coef3 + BB * cjsj.getdSjdH(j) + 1.5 * k * j * cjsj.getCj(j));
  898.                 currentCij[2] += this.xon2a2 * (k * coef2 + BB * cjsj.getdCjdK(j) + 1.5 * h * j * cjsj.getSj(j));
  899.                 currentSij[2] += this.xon2a2 * (k * coef3 + BB * cjsj.getdSjdK(j) - 1.5 * h * j * cjsj.getCj(j));

  901.                 //Coefficients for q and p
  902.                 final double coef4 = CjHK - CjAlphaBeta - j * cjsj.getSj(j);
  903.                 final double coef5 = SjHK - SjAlphaBeta + j * cjsj.getCj(j);
  904.                 currentCij[3] = this.cxo2n2a2 * (-I * CjAlphaGamma + q * coef4);
  905.                 currentSij[3] = this.cxo2n2a2 * (-I * SjAlphaGamma + q * coef5);
  906.                 currentCij[4] = this.cxo2n2a2 * (-CjBetaGamma + p * coef4);
  907.                 currentSij[4] = this.cxo2n2a2 * (-SjBetaGamma + p * coef5);

  909.                 //Coefficients for λ
  910.                 final double coef6 = h * cjsj.getdCjdH(j) + k * cjsj.getdCjdK(j);
  911.                 final double coef7 = h * cjsj.getdSjdH(j) + k * cjsj.getdSjdK(j);
  912.                 currentCij[5] += this.oon2a2 * (-2 * a * cjsj.getdCjdA(j) + coef6 / (X + 1) + X * coef2 - 3 * cjsj.getCj(j));
  913.                 currentSij[5] += this.oon2a2 * (-2 * a * cjsj.getdSjdA(j) + coef7 / (X + 1) + X * coef3 - 3 * cjsj.getSj(j));
  914.             }

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

  921.             // Add the coefficients to the interpolation grid
  922.             slot.cij[j].addGridPoint(date, currentCij);
  923.             slot.sij[j].addGridPoint(date, currentSij);

  925.         }

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

  930.     }

  932.     /**
  933.      * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
  934.      * <p>
  935.      * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
  936.      *  ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
  937.      *  σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
  938.      * </p>
  939.      * @param date target date
  940.      * @param slot slot to which the coefficients belong
  941.      * @return array containing ρ<sub>j</sub> and σ<sub>j</sub>
  942.      */
  943.     private double[][] computeRhoSigmaCoefficients(final AbsoluteDate date, final Slot slot) {
  944.         final CjSjCoefficient cjsjKH = new CjSjCoefficient(k, h);
  945.         final double b = 1. / (1 + B);

  947.         // (-b)<sup>j</sup>
  948.         double mbtj = 1;

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

  953.             //Compute current rho and sigma;
  954.             mbtj *= -b;
  955.             final double coef  = (1 + j * B) * mbtj;
  956.             final double rho   = coef * cjsjKH.getCj(j);
  957.             final double sigma = coef * cjsjKH.getSj(j);

  959.             // Add the coefficients to the interpolation grid
  960.             rhoSigma[j][0] = rho;
  961.             rhoSigma[j][1] = sigma;
  962.         }

  964.         return rhoSigma;

  966.     }

  968.     /** The coefficients used to compute the short-periodic zonal contribution.
  969.      *
  970.      * <p>
  971.      * Those coefficients are given in Danielson paper by expressions 4.1-(20) to 4.1.-(25)
  972.      * </p>
  973.      * <p>
  974.      * The coefficients are: <br>
  975.      * - C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> <br>
  976.      * - ρ<sub>j</sub> and σ<sub>j</sub> <br>
  977.      * - C<sub>i</sub>⁰
  978.      * </p>
  979.      *
  980.      * @author Lucian Barbulescu
  981.      */
  982.     private static class ZonalShortPeriodicCoefficients implements ShortPeriodTerms {

  984.         /** Serializable UID. */
  985.         private static final long serialVersionUID = 20151118L;

  987.         /** Maximum value for j index. */
  988.         private final int maxFrequencyShortPeriodics;

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

  993.         /** All coefficients slots. */
  994.         private final transient TimeSpanMap<Slot> slots;

  996.         /** Constructor.
  997.          * @param maxFrequencyShortPeriodics maximum value for j index
  998.          * @param interpolationPoints number of points used in the interpolation process
  999.          * @param slots all coefficients slots
  1000.          */
  1001.         ZonalShortPeriodicCoefficients(final int maxFrequencyShortPeriodics, final int interpolationPoints,
  1002.                                        final TimeSpanMap<Slot> slots) {

  1004.             // Save parameters
  1005.             this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
  1006.             this.interpolationPoints        = interpolationPoints;
  1007.             this.slots                      = slots;

  1009.         }

  1011.         /** Get the slot valid for some date.
  1012.          * @param meanStates mean states defining the slot
  1013.          * @return slot valid at the specified date
  1014.          */
  1015.         public Slot createSlot(final SpacecraftState... meanStates) {
  1016.             final Slot         slot  = new Slot(maxFrequencyShortPeriodics, interpolationPoints);
  1017.             final AbsoluteDate first = meanStates[0].getDate();
  1018.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  1019.             if (first.compareTo(last) <= 0) {
  1020.                 slots.addValidAfter(slot, first);
  1021.             } else {
  1022.                 slots.addValidBefore(slot, first);
  1023.             }
  1024.             return slot;
  1025.         }

  1027.         /** {@inheritDoc} */
  1028.         @Override
  1029.         public double[] value(final Orbit meanOrbit) {

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

  1034.             // Get the True longitude L
  1035.             final double L = meanOrbit.getLv();

  1037.             // Compute the center
  1038.             final double center = L - meanOrbit.getLM();

  1040.             // Initialize short periodic variations
  1041.             final double[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  1042.             final double[] d = slot.di.value(meanOrbit.getDate());
  1043.             for (int i = 0; i < 6; i++) {
  1044.                 shortPeriodicVariation[i] +=  center * d[i];
  1045.             }

  1047.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1048.                 final double[] c = slot.cij[j].value(meanOrbit.getDate());
  1049.                 final double[] s = slot.sij[j].value(meanOrbit.getDate());
  1050.                 final double cos = FastMath.cos(j * L);
  1051.                 final double sin = FastMath.sin(j * L);
  1052.                 for (int i = 0; i < 6; i++) {
  1053.                     // add corresponding term to the short periodic variation
  1054.                     shortPeriodicVariation[i] += c[i] * cos;
  1055.                     shortPeriodicVariation[i] += s[i] * sin;
  1056.                 }
  1057.             }

  1059.             return shortPeriodicVariation;
  1060.         }

  1062.         /** {@inheritDoc} */
  1063.         @Override
  1064.         public String getCoefficientsKeyPrefix() {
  1065.             return "DSST-central-body-zonal-";
  1066.         }

  1068.         /** {@inheritDoc}
  1069.          * <p>
  1070.          * For zonal terms contributions,there are maxJ cj coefficients,
  1071.          * maxJ sj coefficients and 2 dj coefficients, where maxJ depends
  1072.          * on the orbit. The j index is the integer multiplier for the true
  1073.          * longitude argument in the cj and sj coefficients and the degree
  1074.          * in the polynomial dj coefficients.
  1075.          * </p>
  1076.          */
  1077.         @Override
  1078.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected)
  1079.                 throws OrekitException {

  1081.             // select the coefficients slot
  1082.             final Slot slot = slots.get(date);

  1084.             final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * maxFrequencyShortPeriodics + 2);
  1085.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  1086.             storeIfSelected(coefficients, selected, slot.di.value(date), "d", 1);
  1087.             for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
  1088.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  1089.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  1090.             }
  1091.             return coefficients;

  1093.         }

  1095.         /** Put a coefficient in a map if selected.
  1096.          * @param map map to populate
  1097.          * @param selected set of coefficients that should be put in the map
  1098.          * (empty set means all coefficients are selected)
  1099.          * @param value coefficient value
  1100.          * @param id coefficient identifier
  1101.          * @param indices list of coefficient indices
  1102.          */
  1103.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
  1104.                                      final double[] value, final String id, final int... indices) {
  1105.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  1106.             keyBuilder.append(id);
  1107.             for (int index : indices) {
  1108.                 keyBuilder.append('[').append(index).append(']');
  1109.             }
  1110.             final String key = keyBuilder.toString();
  1111.             if (selected.isEmpty() || selected.contains(key)) {
  1112.                 map.put(key, value);
  1113.             }
  1114.         }

  1116.         /** Replace the instance with a data transfer object for serialization.
  1117.          * @return data transfer object that will be serialized
  1118.          * @exception NotSerializableException if an additional state provider is not serializable
  1119.          */
  1120.         private Object writeReplace() throws NotSerializableException {

  1122.             // slots transitions
  1123.             final SortedSet<TimeSpanMap.Transition<Slot>> transitions     = slots.getTransitions();
  1124.             final AbsoluteDate[]                          transitionDates = new AbsoluteDate[transitions.size()];
  1125.             final Slot[]                                  allSlots        = new Slot[transitions.size() + 1];
  1126.             int i = 0;
  1127.             for (final TimeSpanMap.Transition<Slot> transition : transitions) {
  1128.                 if (i == 0) {
  1129.                     // slot before the first transition
  1130.                     allSlots[i] = transition.getBefore();
  1131.                 }
  1132.                 if (i < transitionDates.length) {
  1133.                     transitionDates[i] = transition.getDate();
  1134.                     allSlots[++i]      = transition.getAfter();
  1135.                 }
  1136.             }

  1138.             return new DataTransferObject(maxFrequencyShortPeriodics, interpolationPoints,
  1139.                                           transitionDates, allSlots);

  1141.         }


  1144.         /** Internal class used only for serialization. */
  1145.         private static class DataTransferObject implements Serializable {

  1147.             /** Serializable UID. */
  1148.             private static final long serialVersionUID = 20170420L;

  1150.             /** Maximum value for j index. */
  1151.             private final int maxFrequencyShortPeriodics;

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

  1156.             /** Transitions dates. */
  1157.             private final AbsoluteDate[] transitionDates;

  1159.             /** All slots. */
  1160.             private final Slot[] allSlots;

  1162.             /** Simple constructor.
  1163.              * @param maxFrequencyShortPeriodics maximum value for j index
  1164.              * @param interpolationPoints number of points used in the interpolation process
  1165.              * @param transitionDates transitions dates
  1166.              * @param allSlots all slots
  1167.              */
  1168.             DataTransferObject(final int maxFrequencyShortPeriodics, final int interpolationPoints,
  1169.                                final AbsoluteDate[] transitionDates, final Slot[] allSlots) {
  1170.                 this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
  1171.                 this.interpolationPoints        = interpolationPoints;
  1172.                 this.transitionDates            = transitionDates;
  1173.                 this.allSlots                   = allSlots;
  1174.             }

  1176.             /** Replace the deserialized data transfer object with a {@link ZonalShortPeriodicCoefficients}.
  1177.              * @return replacement {@link ZonalShortPeriodicCoefficients}
  1178.              */
  1179.             private Object readResolve() {

  1181.                 final TimeSpanMap<Slot> slots = new TimeSpanMap<Slot>(allSlots[0]);
  1182.                 for (int i = 0; i < transitionDates.length; ++i) {
  1183.                     slots.addValidAfter(allSlots[i + 1], transitionDates[i]);
  1184.                 }

  1186.                 return new ZonalShortPeriodicCoefficients(maxFrequencyShortPeriodics,
  1187.                                                           interpolationPoints,
  1188.                                                           slots);

  1190.             }

  1192.         }

  1194.     }

  1196.     /** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
  1197.      *  <p>
  1198.      *  Those coefficients are given in Danielson paper by expressions 4.1-(13) to 4.1.-(16b)
  1199.      *  </p>
  1200.      */
  1201.     private class FourierCjSjCoefficients {

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

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

  1209.         /** Maximum possible value for n. */
  1210.         private final int nMax;

  1212.         /** Maximum possible value for s. */
  1213.         private final int sMax;

  1215.         /** Maximum possible value for j. */
  1216.         private final int jMax;

  1218.         /** The C<sup>j</sup> coefficients and their derivatives.
  1219.          * <p>
  1220.          * Each column of the matrix contains the following values: <br/>
  1221.          * - C<sup>j</sup> <br/>
  1222.          * - dC<sup>j</sup> / da <br/>
  1223.          * - dC<sup>j</sup> / dh <br/>
  1224.          * - dC<sup>j</sup> / dk <br/>
  1225.          * - dC<sup>j</sup> / dα <br/>
  1226.          * - dC<sup>j</sup> / dβ <br/>
  1227.          * - dC<sup>j</sup> / dγ <br/>
  1228.          * </p>
  1229.          */
  1230.         private final double[][] cCoef;

  1232.         /** The S<sup>j</sup> coefficients and their derivatives.
  1233.          * <p>
  1234.          * Each column of the matrix contains the following values: <br/>
  1235.          * - S<sup>j</sup> <br/>
  1236.          * - dS<sup>j</sup> / da <br/>
  1237.          * - dS<sup>j</sup> / dh <br/>
  1238.          * - dS<sup>j</sup> / dk <br/>
  1239.          * - dS<sup>j</sup> / dα <br/>
  1240.          * - dS<sup>j</sup> / dβ <br/>
  1241.          * - dS<sup>j</sup> / dγ <br/>
  1242.          * </p>
  1243.          */
  1244.         private final double[][] sCoef;

  1246.         /** h * &Chi;³. */
  1247.         private final double hXXX;
  1248.         /** k * &Chi;³. */
  1249.         private final double kXXX;

  1251.         /** Create a set of C<sup>j</sup> and the S<sup>j</sup> coefficients.
  1252.          *  @param date the current date
  1253.          *  @param nMax maximum possible value for n
  1254.          *  @param sMax maximum possible value for s
  1255.          *  @param jMax maximum possible value for j
  1256.          * @throws OrekitException if an error occurs while generating the coefficients
  1257.          */
  1258.         FourierCjSjCoefficients(final AbsoluteDate date,
  1259.                                 final int nMax, final int sMax, final int jMax)
  1260.                 throws OrekitException {
  1261.             this.ghijCoef = new GHIJjsPolynomials(k, h, alpha, beta);
  1262.             // Qns coefficients
  1263.             final double[][] Qns  = CoefficientsFactory.computeQns(gamma, nMax, nMax);

  1265.             this.lnsCoef = new LnsCoefficients(nMax, nMax, Qns, Vns, roa);
  1266.             this.nMax = nMax;
  1267.             this.sMax = sMax;
  1268.             this.jMax = jMax;

  1270.             // compute the common factors that depends on the mean elements
  1271.             this.hXXX = h * XXX;
  1272.             this.kXXX = k * XXX;

  1274.             this.cCoef = new double[7][jMax + 1];
  1275.             this.sCoef = new double[7][jMax + 1];

  1277.             for (int s = 0; s <= sMax; s++) {
  1278.                 //Initialise the Hansen roots
  1279.                 hansenObjects[s].computeInitValues(X);
  1280.             }
  1281.             generateCoefficients(date);
  1282.         }

  1284.         /** Generate all coefficients.
  1285.          * @param date the current date
  1286.          * @throws OrekitException if an error occurs while generating the coefficients
  1287.          */
  1288.         private void generateCoefficients(final AbsoluteDate date) throws OrekitException {
  1289.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);
  1290.             for (int j = 1; j <= jMax; j++) {

  1292.                 //init arrays
  1293.                 for (int i = 0; i <= 6; i++) {
  1294.                     cCoef[i][j] = 0.;
  1295.                     sCoef[i][j] = 0.;
  1296.                 }

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

  1300.                     //compute first double sum where s: j -> N-1 and n: s+1 -> N
  1301.                     for (int s = j; s <= FastMath.min(nMax - 1, sMax); s++) {
  1302.                         // j - s
  1303.                         final int jms = j - s;
  1304.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  1305.                         final int d0smj = (s == j) ? 1 : 2;

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

  1314.                                 final double hjs = ghijCoef.getHjs(s, -jms);
  1315.                                 final double dHjsdh = ghijCoef.getdHjsdh(s, -jms);
  1316.                                 final double dHjsdk = ghijCoef.getdHjsdk(s, -jms);
  1317.                                 final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, -jms);
  1318.                                 final double dHjsdBeta = ghijCoef.getdHjsdBeta(s, -jms);

  1320.                                 final double gjs = ghijCoef.getGjs(s, -jms);
  1321.                                 final double dGjsdh = ghijCoef.getdGjsdh(s, -jms);
  1322.                                 final double dGjsdk = ghijCoef.getdGjsdk(s, -jms);
  1323.                                 final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, -jms);
  1324.                                 final double dGjsdBeta = ghijCoef.getdGjsdBeta(s, -jms);

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

  1329.                                 // K₀<sup>-n-1,s</sup>
  1330.                                 final double kns   = hansenObjects[s].getValue(-n - 1, X);
  1331.                                 final double dkns  = hansenObjects[s].getDerivative(-n - 1, X);

  1333.                                 final double coef0 = d0smj * jn;
  1334.                                 final double coef1 = coef0 * lns;
  1335.                                 final double coef2 = coef1 * kns;
  1336.                                 final double coef3 = coef2 * hjs;
  1337.                                 final double coef4 = coef2 * gjs;

  1339.                                 // Add the term to the coefficients
  1340.                                 cCoef[0][j] += coef3;
  1341.                                 cCoef[1][j] += coef3 * (n + 1);
  1342.                                 cCoef[2][j] += coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
  1343.                                 cCoef[3][j] += coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
  1344.                                 cCoef[4][j] += coef2 * dHjsdAlpha;
  1345.                                 cCoef[5][j] += coef2 * dHjsdBeta;
  1346.                                 cCoef[6][j] += coef0 * dlns * kns * hjs;

  1348.                                 sCoef[0][j] += coef4;
  1349.                                 sCoef[1][j] += coef4 * (n + 1);
  1350.                                 sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
  1351.                                 sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
  1352.                                 sCoef[4][j] += coef2 * dGjsdAlpha;
  1353.                                 sCoef[5][j] += coef2 * dGjsdBeta;
  1354.                                 sCoef[6][j] += coef0 * dlns * kns * gjs;
  1355.                             }
  1356.                         }
  1357.                     }

  1359.                     //compute second double sum where s: 0 -> N-j and n: max(j+s, j+1) -> N
  1360.                     for (int s = 0; s <= FastMath.min(nMax - j, sMax); s++) {
  1361.                         // j + s
  1362.                         final int jps = j + s;
  1363.                         // Kronecker symbols (2 - delta(0,j+s))
  1364.                         final double d0spj = (s == -j) ? 1 : 2;

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

  1373.                                 final double hjs = ghijCoef.getHjs(s, jps);
  1374.                                 final double dHjsdh = ghijCoef.getdHjsdh(s, jps);
  1375.                                 final double dHjsdk = ghijCoef.getdHjsdk(s, jps);
  1376.                                 final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, jps);
  1377.                                 final double dHjsdBeta = ghijCoef.getdHjsdBeta(s, jps);

  1379.                                 final double gjs = ghijCoef.getGjs(s, jps);
  1380.                                 final double dGjsdh = ghijCoef.getdGjsdh(s, jps);
  1381.                                 final double dGjsdk = ghijCoef.getdGjsdk(s, jps);
  1382.                                 final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, jps);
  1383.                                 final double dGjsdBeta = ghijCoef.getdGjsdBeta(s, jps);

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

  1388.                                 // K₀<sup>-n-1,s</sup>
  1389.                                 final double kns   = hansenObjects[s].getValue(-n - 1, X);
  1390.                                 final double dkns  = hansenObjects[s].getDerivative(-n - 1, X);

  1392.                                 final double coef0 = d0spj * jn;
  1393.                                 final double coef1 = coef0 * lns;
  1394.                                 final double coef2 = coef1 * kns;

  1396.                                 final double coef3 = coef2 * hjs;
  1397.                                 final double coef4 = coef2 * gjs;

  1399.                                 // Add the term to the coefficients
  1400.                                 cCoef[0][j] -= coef3;
  1401.                                 cCoef[1][j] -= coef3 * (n + 1);
  1402.                                 cCoef[2][j] -= coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
  1403.                                 cCoef[3][j] -= coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
  1404.                                 cCoef[4][j] -= coef2 * dHjsdAlpha;
  1405.                                 cCoef[5][j] -= coef2 * dHjsdBeta;
  1406.                                 cCoef[6][j] -= coef0 * dlns * kns * hjs;

  1408.                                 sCoef[0][j] += coef4;
  1409.                                 sCoef[1][j] += coef4 * (n + 1);
  1410.                                 sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
  1411.                                 sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
  1412.                                 sCoef[4][j] += coef2 * dGjsdAlpha;
  1413.                                 sCoef[5][j] += coef2 * dGjsdBeta;
  1414.                                 sCoef[6][j] += coef0 * dlns * kns * gjs;
  1415.                             }
  1416.                         }
  1417.                     }

  1419.                     //compute third double sum where s: 1 -> j and  n: j+1 -> N
  1420.                     for (int s = 1; s <= FastMath.min(j, sMax); s++) {
  1421.                         // j - s
  1422.                         final int jms = j - s;
  1423.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  1424.                         final int d0smj = (s == j) ? 1 : 2;

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

  1433.                                 final double ijs = ghijCoef.getIjs(s, jms);
  1434.                                 final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  1435.                                 final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  1436.                                 final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  1437.                                 final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  1439.                                 final double jjs = ghijCoef.getJjs(s, jms);
  1440.                                 final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  1441.                                 final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  1442.                                 final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  1443.                                 final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

  1448.                                 // K₀<sup>-n-1,s</sup>
  1449.                                 final double kns   = hansenObjects[s].getValue(-n - 1, X);
  1450.                                 final double dkns  = hansenObjects[s].getDerivative(-n - 1, X);

  1452.                                 final double coef0 = d0smj * jn;
  1453.                                 final double coef1 = coef0 * lns;
  1454.                                 final double coef2 = coef1 * kns;

  1456.                                 final double coef3 = coef2 * ijs;
  1457.                                 final double coef4 = coef2 * jjs;

  1459.                                 // Add the term to the coefficients
  1460.                                 cCoef[0][j] -= coef3;
  1461.                                 cCoef[1][j] -= coef3 * (n + 1);
  1462.                                 cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  1463.                                 cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  1464.                                 cCoef[4][j] -= coef2 * dIjsdAlpha;
  1465.                                 cCoef[5][j] -= coef2 * dIjsdBeta;
  1466.                                 cCoef[6][j] -= coef0 * dlns * kns * ijs;

  1468.                                 sCoef[0][j] += coef4;
  1469.                                 sCoef[1][j] += coef4 * (n + 1);
  1470.                                 sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  1471.                                 sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  1472.                                 sCoef[4][j] += coef2 * dJjsdAlpha;
  1473.                                 sCoef[5][j] += coef2 * dJjsdBeta;
  1474.                                 sCoef[6][j] += coef0 * dlns * kns * jjs;
  1475.                             }
  1476.                         }
  1477.                     }
  1478.                 }

  1480.                 if (isBetween(j, 2, nMax)) {
  1481.                     //add first term
  1482.                     // J<sub>j</sub>
  1483.                     final double jj = -harmonics.getUnnormalizedCnm(j, 0);
  1484.                     double kns = hansenObjects[0].getValue(-j - 1, X);
  1485.                     double dkns = hansenObjects[0].getDerivative(-j - 1, X);

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

  1490.                     final double hjs = ghijCoef.getHjs(0, j);
  1491.                     final double dHjsdh = ghijCoef.getdHjsdh(0, j);
  1492.                     final double dHjsdk = ghijCoef.getdHjsdk(0, j);
  1493.                     final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(0, j);
  1494.                     final double dHjsdBeta = ghijCoef.getdHjsdBeta(0, j);

  1496.                     final double gjs = ghijCoef.getGjs(0, j);
  1497.                     final double dGjsdh = ghijCoef.getdGjsdh(0, j);
  1498.                     final double dGjsdk = ghijCoef.getdGjsdk(0, j);
  1499.                     final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(0, j);
  1500.                     final double dGjsdBeta = ghijCoef.getdGjsdBeta(0, j);

  1502.                     // 2 * J<sub>j</sub> * K₀<sup>-j-1,0</sup> * L<sub>j</sub><sup>j</sup>
  1503.                     double coef0 = 2 * jj;
  1504.                     double coef1 = coef0 * lns;
  1505.                     double coef2 = coef1 * kns;

  1507.                     double coef3 = coef2 * hjs;
  1508.                     double coef4 = coef2 * gjs;

  1510.                     // Add the term to the coefficients
  1511.                     cCoef[0][j] -= coef3;
  1512.                     cCoef[1][j] -= coef3 * (j + 1);
  1513.                     cCoef[2][j] -= coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
  1514.                     cCoef[3][j] -= coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
  1515.                     cCoef[4][j] -= coef2 * dHjsdAlpha;
  1516.                     cCoef[5][j] -= coef2 * dHjsdBeta;
  1517.                     //no contribution to cCoef[6][j] because dlns is 0

  1519.                     sCoef[0][j] += coef4;
  1520.                     sCoef[1][j] += coef4 * (j + 1);
  1521.                     sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
  1522.                     sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
  1523.                     sCoef[4][j] += coef2 * dGjsdAlpha;
  1524.                     sCoef[5][j] += coef2 * dGjsdBeta;
  1525.                     //no contribution to sCoef[6][j] because dlns is 0

  1527.                     //compute simple sum where s: 1 -> j-1
  1528.                     for (int s = 1; s <= FastMath.min(j - 1, sMax); s++) {
  1529.                         // j - s
  1530.                         final int jms = j - s;
  1531.                         // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  1532.                         final int d0smj = (s == j) ? 1 : 2;

  1534.                         // if s is odd, then the term is equal to zero due to the factor Vj,s-j
  1535.                         if (s % 2 == 0) {
  1536.                             // (2 - delta(0,j-s)) * J<sub>j</sub> * K₀<sup>-j-1,s</sup> * L<sub>j</sub><sup>j-s</sup>
  1537.                             kns   = hansenObjects[s].getValue(-j - 1, X);
  1538.                             dkns  = hansenObjects[s].getDerivative(-j - 1, X);

  1540.                             lns = lnsCoef.getLns(j, jms);
  1541.                             final double dlns = lnsCoef.getdLnsdGamma(j, jms);

  1543.                             final double ijs = ghijCoef.getIjs(s, jms);
  1544.                             final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  1545.                             final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  1546.                             final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  1547.                             final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  1549.                             final double jjs = ghijCoef.getJjs(s, jms);
  1550.                             final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  1551.                             final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  1552.                             final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  1553.                             final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

  1555.                             coef0 = d0smj * jj;
  1556.                             coef1 = coef0 * lns;
  1557.                             coef2 = coef1 * kns;

  1559.                             coef3 = coef2 * ijs;
  1560.                             coef4 = coef2 * jjs;

  1562.                             // Add the term to the coefficients
  1563.                             cCoef[0][j] -= coef3;
  1564.                             cCoef[1][j] -= coef3 * (j + 1);
  1565.                             cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  1566.                             cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  1567.                             cCoef[4][j] -= coef2 * dIjsdAlpha;
  1568.                             cCoef[5][j] -= coef2 * dIjsdBeta;
  1569.                             cCoef[6][j] -= coef0 * dlns * kns * ijs;

  1571.                             sCoef[0][j] += coef4;
  1572.                             sCoef[1][j] += coef4 * (j + 1);
  1573.                             sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  1574.                             sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  1575.                             sCoef[4][j] += coef2 * dJjsdAlpha;
  1576.                             sCoef[5][j] += coef2 * dJjsdBeta;
  1577.                             sCoef[6][j] += coef0 * dlns * kns * jjs;
  1578.                         }
  1579.                     }
  1580.                 }

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

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

  1588.                     //if j is even
  1589.                     if (j % 2 == 0) {
  1590.                         //compute first double sum where s: j-min(j-1,N) -> j/2-1 and n: j-s -> min(j-1,N)
  1591.                         for (int s = j - minjm1on; s <= FastMath.min(j / 2 - 1, sMax); s++) {
  1592.                             // j - s
  1593.                             final int jms = j - s;
  1594.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  1595.                             final int d0smj = (s == j) ? 1 : 2;

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

  1604.                                     final double ijs = ghijCoef.getIjs(s, jms);
  1605.                                     final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  1606.                                     final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  1607.                                     final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  1608.                                     final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  1610.                                     final double jjs = ghijCoef.getJjs(s, jms);
  1611.                                     final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  1612.                                     final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  1613.                                     final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  1614.                                     final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

  1619.                                     // K₀<sup>-n-1,s</sup>
  1620.                                     final double kns   = hansenObjects[s].getValue(-n - 1, X);
  1621.                                     final double dkns  = hansenObjects[s].getDerivative(-n - 1, X);

  1623.                                     final double coef0 = d0smj * jn;
  1624.                                     final double coef1 = coef0 * lns;
  1625.                                     final double coef2 = coef1 * kns;

  1627.                                     final double coef3 = coef2 * ijs;
  1628.                                     final double coef4 = coef2 * jjs;

  1630.                                     // Add the term to the coefficients
  1631.                                     cCoef[0][j] -= coef3;
  1632.                                     cCoef[1][j] -= coef3 * (n + 1);
  1633.                                     cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  1634.                                     cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  1635.                                     cCoef[4][j] -= coef2 * dIjsdAlpha;
  1636.                                     cCoef[5][j] -= coef2 * dIjsdBeta;
  1637.                                     cCoef[6][j] -= coef0 * dlns * kns * ijs;

  1639.                                     sCoef[0][j] += coef4;
  1640.                                     sCoef[1][j] += coef4 * (n + 1);
  1641.                                     sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  1642.                                     sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  1643.                                     sCoef[4][j] += coef2 * dJjsdAlpha;
  1644.                                     sCoef[5][j] += coef2 * dJjsdBeta;
  1645.                                     sCoef[6][j] += coef0 * dlns * kns * jjs;
  1646.                                 }
  1647.                             }
  1648.                         }

  1650.                         //compute second double sum where s: j/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  1651.                         for (int s = j / 2; s <=  FastMath.min(minjm1on - 1, sMax); s++) {
  1652.                             // j - s
  1653.                             final int jms = j - s;
  1654.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  1655.                             final int d0smj = (s == j) ? 1 : 2;

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

  1664.                                     final double ijs = ghijCoef.getIjs(s, jms);
  1665.                                     final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  1666.                                     final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  1667.                                     final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  1668.                                     final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  1670.                                     final double jjs = ghijCoef.getJjs(s, jms);
  1671.                                     final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  1672.                                     final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  1673.                                     final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  1674.                                     final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

  1679.                                     // K₀<sup>-n-1,s</sup>
  1680.                                     final double kns   = hansenObjects[s].getValue(-n - 1, X);
  1681.                                     final double dkns  = hansenObjects[s].getDerivative(-n - 1, X);

  1683.                                     final double coef0 = d0smj * jn;
  1684.                                     final double coef1 = coef0 * lns;
  1685.                                     final double coef2 = coef1 * kns;

  1687.                                     final double coef3 = coef2 * ijs;
  1688.                                     final double coef4 = coef2 * jjs;

  1690.                                     // Add the term to the coefficients
  1691.                                     cCoef[0][j] -= coef3;
  1692.                                     cCoef[1][j] -= coef3 * (n + 1);
  1693.                                     cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  1694.                                     cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  1695.                                     cCoef[4][j] -= coef2 * dIjsdAlpha;
  1696.                                     cCoef[5][j] -= coef2 * dIjsdBeta;
  1697.                                     cCoef[6][j] -= coef0 * dlns * kns * ijs;

  1699.                                     sCoef[0][j] += coef4;
  1700.                                     sCoef[1][j] += coef4 * (n + 1);
  1701.                                     sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  1702.                                     sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  1703.                                     sCoef[4][j] += coef2 * dJjsdAlpha;
  1704.                                     sCoef[5][j] += coef2 * dJjsdBeta;
  1705.                                     sCoef[6][j] += coef0 * dlns * kns * jjs;
  1706.                                 }
  1707.                             }
  1708.                         }
  1709.                     }

  1711.                     //if j is odd
  1712.                     else {
  1713.                         //compute first double sum where s: (j-1)/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
  1714.                         for (int s = (j - 1) / 2; s <= FastMath.min(minjm1on - 1, sMax); s++) {
  1715.                             // j - s
  1716.                             final int jms = j - s;
  1717.                             // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  1718.                             final int d0smj = (s == j) ? 1 : 2;

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

  1727.                                     final double ijs = ghijCoef.getIjs(s, jms);
  1728.                                     final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  1729.                                     final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  1730.                                     final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  1731.                                     final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  1733.                                     final double jjs = ghijCoef.getJjs(s, jms);
  1734.                                     final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  1735.                                     final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  1736.                                     final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  1737.                                     final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

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

  1744.                                     final double kns = hansenObjects[s].getValue(-n - 1, X);
  1745.                                     final double dkns  = hansenObjects[s].getDerivative(-n - 1, X);

  1747.                                     final double coef0 = d0smj * jn;
  1748.                                     final double coef1 = coef0 * lns;
  1749.                                     final double coef2 = coef1 * kns;

  1751.                                     final double coef3 = coef2 * ijs;
  1752.                                     final double coef4 = coef2 * jjs;

  1754.                                     // Add the term to the coefficients
  1755.                                     cCoef[0][j] -= coef3;
  1756.                                     cCoef[1][j] -= coef3 * (n + 1);
  1757.                                     cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  1758.                                     cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  1759.                                     cCoef[4][j] -= coef2 * dIjsdAlpha;
  1760.                                     cCoef[5][j] -= coef2 * dIjsdBeta;
  1761.                                     cCoef[6][j] -= coef0 * dlns * kns * ijs;

  1763.                                     sCoef[0][j] += coef4;
  1764.                                     sCoef[1][j] += coef4 * (n + 1);
  1765.                                     sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  1766.                                     sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  1767.                                     sCoef[4][j] += coef2 * dJjsdAlpha;
  1768.                                     sCoef[5][j] += coef2 * dJjsdBeta;
  1769.                                     sCoef[6][j] += coef0 * dlns * kns * jjs;
  1770.                                 }
  1771.                             }
  1772.                         }

  1774.                         //the second double sum is added only if N >= 4 and j between 5 and 2*N-3
  1775.                         if (nMax >= 4 && isBetween(j, 5, 2 * nMax - 3)) {
  1776.                             //compute second double sum where s: j-min(j-1,N) -> (j-3)/2 and n: j-s -> min(j-1,N)
  1777.                             for (int s = j - minjm1on; s <= FastMath.min((j - 3) / 2, sMax); s++) {
  1778.                                 // j - s
  1779.                                 final int jms = j - s;
  1780.                                 // Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
  1781.                                 final int d0smj = (s == j) ? 1 : 2;

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

  1790.                                         final double ijs = ghijCoef.getIjs(s, jms);
  1791.                                         final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
  1792.                                         final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
  1793.                                         final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
  1794.                                         final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);

  1796.                                         final double jjs = ghijCoef.getJjs(s, jms);
  1797.                                         final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
  1798.                                         final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
  1799.                                         final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
  1800.                                         final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);

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

  1805.                                         // K₀<sup>-n-1,s</sup>
  1806.                                         final double kns   = hansenObjects[s].getValue(-n - 1, X);
  1807.                                         final double dkns  = hansenObjects[s].getDerivative(-n - 1, X);

  1809.                                         final double coef0 = d0smj * jn;
  1810.                                         final double coef1 = coef0 * lns;
  1811.                                         final double coef2 = coef1 * kns;

  1813.                                         final double coef3 = coef2 * ijs;
  1814.                                         final double coef4 = coef2 * jjs;

  1816.                                         // Add the term to the coefficients
  1817.                                         cCoef[0][j] -= coef3;
  1818.                                         cCoef[1][j] -= coef3 * (n + 1);
  1819.                                         cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
  1820.                                         cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
  1821.                                         cCoef[4][j] -= coef2 * dIjsdAlpha;
  1822.                                         cCoef[5][j] -= coef2 * dIjsdBeta;
  1823.                                         cCoef[6][j] -= coef0 * dlns * kns * ijs;

  1825.                                         sCoef[0][j] += coef4;
  1826.                                         sCoef[1][j] += coef4 * (n + 1);
  1827.                                         sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
  1828.                                         sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
  1829.                                         sCoef[4][j] += coef2 * dJjsdAlpha;
  1830.                                         sCoef[5][j] += coef2 * dJjsdBeta;
  1831.                                         sCoef[6][j] += coef0 * dlns * kns * jjs;
  1832.                                     }
  1833.                                 }
  1834.                             }
  1835.                         }
  1836.                     }
  1837.                 }

  1839.                 cCoef[0][j] *= -muoa / j;
  1840.                 cCoef[1][j] *=  muoa / ( j * a );
  1841.                 cCoef[2][j] *= -muoa / j;
  1842.                 cCoef[3][j] *= -muoa / j;
  1843.                 cCoef[4][j] *= -muoa / j;
  1844.                 cCoef[5][j] *= -muoa / j;
  1845.                 cCoef[6][j] *= -muoa / j;

  1847.                 sCoef[0][j] *= -muoa / j;
  1848.                 sCoef[1][j] *=  muoa / ( j * a );
  1849.                 sCoef[2][j] *= -muoa / j;
  1850.                 sCoef[3][j] *= -muoa / j;
  1851.                 sCoef[4][j] *= -muoa / j;
  1852.                 sCoef[5][j] *= -muoa / j;
  1853.                 sCoef[6][j] *= -muoa / j;

  1855.             }
  1856.         }

  1858.         /** Check if an index is within the accepted interval.
  1859.          *
  1860.          * @param index the index to check
  1861.          * @param lowerBound the lower bound of the interval
  1862.          * @param upperBound the upper bound of the interval
  1863.          * @return true if the index is between the lower and upper bounds, false otherwise
  1864.          */
  1865.         private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
  1866.             return index >= lowerBound && index <= upperBound;
  1867.         }

  1869.         /**Get the value of C<sup>j</sup>.
  1870.          *
  1871.          * @param j j index
  1872.          * @return C<sup>j</sup>
  1873.          */
  1874.         public double getCj(final int j) {
  1875.             return cCoef[0][j];
  1876.         }

  1878.         /**Get the value of dC<sup>j</sup> / da.
  1879.          *
  1880.          * @param j j index
  1881.          * @return dC<sup>j</sup> / da
  1882.          */
  1883.         public double getdCjdA(final int j) {
  1884.             return cCoef[1][j];
  1885.         }

  1887.         /**Get the value of dC<sup>j</sup> / dh.
  1888.          *
  1889.          * @param j j index
  1890.          * @return dC<sup>j</sup> / dh
  1891.          */
  1892.         public double getdCjdH(final int j) {
  1893.             return cCoef[2][j];
  1894.         }

  1896.         /**Get the value of dC<sup>j</sup> / dk.
  1897.          *
  1898.          * @param j j index
  1899.          * @return dC<sup>j</sup> / dk
  1900.          */
  1901.         public double getdCjdK(final int j) {
  1902.             return cCoef[3][j];
  1903.         }

  1905.         /**Get the value of dC<sup>j</sup> / dα.
  1906.          *
  1907.          * @param j j index
  1908.          * @return dC<sup>j</sup> / dα
  1909.          */
  1910.         public double getdCjdAlpha(final int j) {
  1911.             return cCoef[4][j];
  1912.         }

  1914.         /**Get the value of dC<sup>j</sup> / dβ.
  1915.          *
  1916.          * @param j j index
  1917.          * @return dC<sup>j</sup> / dβ
  1918.          */
  1919.         public double getdCjdBeta(final int j) {
  1920.             return cCoef[5][j];
  1921.         }

  1923.         /**Get the value of dC<sup>j</sup> / dγ.
  1924.          *
  1925.          * @param j j index
  1926.          * @return dC<sup>j</sup> / dγ
  1927.          */
  1928.         public double getdCjdGamma(final int j) {
  1929.             return cCoef[6][j];
  1930.         }

  1932.         /**Get the value of S<sup>j</sup>.
  1933.          *
  1934.          * @param j j index
  1935.          * @return S<sup>j</sup>
  1936.          */
  1937.         public double getSj(final int j) {
  1938.             return sCoef[0][j];
  1939.         }

  1941.         /**Get the value of dS<sup>j</sup> / da.
  1942.          *
  1943.          * @param j j index
  1944.          * @return dS<sup>j</sup> / da
  1945.          */
  1946.         public double getdSjdA(final int j) {
  1947.             return sCoef[1][j];
  1948.         }

  1950.         /**Get the value of dS<sup>j</sup> / dh.
  1951.          *
  1952.          * @param j j index
  1953.          * @return dS<sup>j</sup> / dh
  1954.          */
  1955.         public double getdSjdH(final int j) {
  1956.             return sCoef[2][j];
  1957.         }

  1959.         /**Get the value of dS<sup>j</sup> / dk.
  1960.          *
  1961.          * @param j j index
  1962.          * @return dS<sup>j</sup> / dk
  1963.          */
  1964.         public double getdSjdK(final int j) {
  1965.             return sCoef[3][j];
  1966.         }

  1968.         /**Get the value of dS<sup>j</sup> / dα.
  1969.          *
  1970.          * @param j j index
  1971.          * @return dS<sup>j</sup> / dα
  1972.          */
  1973.         public double getdSjdAlpha(final int j) {
  1974.             return sCoef[4][j];
  1975.         }

  1977.         /**Get the value of dS<sup>j</sup> / dβ.
  1978.          *
  1979.          * @param j j index
  1980.          * @return dS<sup>j</sup> / dβ
  1981.          */
  1982.         public double getdSjdBeta(final int j) {
  1983.             return sCoef[5][j];
  1984.         }

  1986.         /**Get the value of dS<sup>j</sup> /  dγ.
  1987.          *
  1988.          * @param j j index
  1989.          * @return dS<sup>j</sup> /  dγ
  1990.          */
  1991.         public double getdSjdGamma(final int j) {
  1992.             return sCoef[6][j];
  1993.         }
  1994.     }

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

  1999.         /** Serializable UID. */
  2000.         private static final long serialVersionUID = 20160319L;

  2002.         /**The coefficients D<sub>i</sub>.
  2003.          * <p>
  2004.          * i corresponds to the equinoctial element, as follows:
  2005.          * - i=0 for a <br/>
  2006.          * - i=1 for k <br/>
  2007.          * - i=2 for h <br/>
  2008.          * - i=3 for q <br/>
  2009.          * - i=4 for p <br/>
  2010.          * - i=5 for λ <br/>
  2011.          * </p>
  2012.          */
  2013.         private final ShortPeriodicsInterpolatedCoefficient di;

  2015.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  2016.          * <p>
  2017.          * The constant term C<sub>i</sub>⁰ is also stored in this variable at index j = 0 <br>
  2018.          * The index order is cij[j][i] <br/>
  2019.          * i corresponds to the equinoctial element, as follows: <br/>
  2020.          * - i=0 for a <br/>
  2021.          * - i=1 for k <br/>
  2022.          * - i=2 for h <br/>
  2023.          * - i=3 for q <br/>
  2024.          * - i=4 for p <br/>
  2025.          * - i=5 for λ <br/>
  2026.          * </p>
  2027.          */
  2028.         private final ShortPeriodicsInterpolatedCoefficient[] cij;

  2030.         /** The coefficients S<sub>i</sub><sup>j</sup>.
  2031.          * <p>
  2032.          * The index order is sij[j][i] <br/>
  2033.          * i corresponds to the equinoctial element, as follows: <br/>
  2034.          * - i=0 for a <br/>
  2035.          * - i=1 for k <br/>
  2036.          * - i=2 for h <br/>
  2037.          * - i=3 for q <br/>
  2038.          * - i=4 for p <br/>
  2039.          * - i=5 for λ <br/>
  2040.          * </p>
  2041.          */
  2042.         private final ShortPeriodicsInterpolatedCoefficient[] sij;

  2044.         /** Simple constructor.
  2045.          *  @param maxFrequencyShortPeriodics maximum value for j index
  2046.          *  @param interpolationPoints number of points used in the interpolation process
  2047.          */
  2048.         Slot(final int maxFrequencyShortPeriodics, final int interpolationPoints) {

  2050.             final int rows = maxFrequencyShortPeriodics + 1;
  2051.             di  = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  2052.             cij = new ShortPeriodicsInterpolatedCoefficient[rows];
  2053.             sij = new ShortPeriodicsInterpolatedCoefficient[rows];

  2055.             //Initialize the arrays
  2056.             for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
  2057.                 cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  2058.                 sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  2059.             }

  2061.         }

  2063.     }

  2065. }
Created with JaCoCo 0.7.9.201702052155