AbstractGaussianContribution.java

  1. /* Copyright 2002-2017 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;

  18. import java.io.NotSerializableException;
  19. import java.io.Serializable;
  20. import java.util.ArrayList;
  21. import java.util.HashMap;
  22. import java.util.List;
  23. import java.util.Map;
  24. import java.util.Set;
  25. import java.util.SortedSet;

  26. import org.hipparchus.analysis.UnivariateVectorFunction;
  27. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  28. import org.hipparchus.util.FastMath;
  29. import org.orekit.attitudes.Attitude;
  30. import org.orekit.attitudes.AttitudeProvider;
  31. import org.orekit.errors.OrekitException;
  32. import org.orekit.errors.OrekitExceptionWrapper;
  33. import org.orekit.forces.ForceModel;
  34. import org.orekit.orbits.EquinoctialOrbit;
  35. import org.orekit.orbits.Orbit;
  36. import org.orekit.orbits.OrbitType;
  37. import org.orekit.orbits.PositionAngle;
  38. import org.orekit.propagation.SpacecraftState;
  39. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  40. import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
  41. import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
  42. import org.orekit.time.AbsoluteDate;
  43. import org.orekit.utils.TimeSpanMap;

  44. /** Common handling of {@link DSSTForceModel} methods for Gaussian contributions to DSST propagation.
  45.  * <p>
  46.  * This abstract class allows to provide easily a subset of {@link DSSTForceModel} methods
  47.  * for specific Gaussian contributions.
  48.  * </p><p>
  49.  * This class implements the notion of numerical averaging of the DSST theory.
  50.  * Numerical averaging is mainly used for non-conservative disturbing forces such as
  51.  * atmospheric drag and solar radiation pressure.
  52.  * </p><p>
  53.  * Gaussian contributions can be expressed as: da<sub>i</sub>/dt = δa<sub>i</sub>/δv . q<br>
  54.  * where:
  55.  * <ul>
  56.  * <li>a<sub>i</sub> are the six equinoctial elements</li>
  57.  * <li>v is the velocity vector</li>
  58.  * <li>q is the perturbing acceleration due to the considered force</li>
  59.  * </ul>
  60.  *
  61.  * <p> The averaging process and other considerations lead to integrate this contribution
  62.  * over the true longitude L possibly taking into account some limits.
  63.  *
  64.  * <p> To create a numerically averaged contribution, one needs only to provide a
  65.  * {@link ForceModel} and to implement in the derived class the method:
  66.  * {@link #getLLimits(SpacecraftState)}.
  67.  * </p>
  68.  * @author Pascal Parraud
  69.  */
  70. public abstract class AbstractGaussianContribution implements DSSTForceModel {

  71.     /** Available orders for Gauss quadrature. */
  72.     private static final int[] GAUSS_ORDER = {12, 16, 20, 24, 32, 40, 48};

  73.     /** Max rank in Gauss quadrature orders array. */
  74.     private static final int MAX_ORDER_RANK = GAUSS_ORDER.length - 1;

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

  77.     /** Maximum value for j index. */
  78.     private static final int JMAX = 12;

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

  93.     // CHECKSTYLE: stop VisibilityModifierCheck

  94.     /** a. */
  95.     protected double a;
  96.     /** e<sub>x</sub>. */
  97.     protected double k;
  98.     /** e<sub>y</sub>. */
  99.     protected double h;
  100.     /** h<sub>x</sub>. */
  101.     protected double q;
  102.     /** h<sub>y</sub>. */
  103.     protected double p;

  104.     /** Eccentricity. */
  105.     protected double ecc;

  106.     /** Kepler mean motion: n = sqrt(μ / a³). */
  107.     protected double n;

  108.     /** Mean longitude. */
  109.     protected double lm;

  110.     /** Equinoctial frame f vector. */
  111.     protected Vector3D f;
  112.     /** Equinoctial frame g vector. */
  113.     protected Vector3D g;
  114.     /** Equinoctial frame w vector. */
  115.     protected Vector3D w;

  116.     /** A = sqrt(μ * a). */
  117.     protected double A;
  118.     /** B = sqrt(1 - h² - k²). */
  119.     protected double B;
  120.     /** C = 1 + p² + q². */
  121.     protected double C;

  122.     /** 2 / (n² * a) . */
  123.     protected double ton2a;
  124.     /** 1 / A .*/
  125.     protected double ooA;
  126.     /** 1 / (A * B) .*/
  127.     protected double ooAB;
  128.     /** C / (2 * A * B) .*/
  129.     protected double co2AB;
  130.     /** 1 / (1 + B) .*/
  131.     protected double ooBpo;
  132.     /** 1 / μ .*/
  133.     protected double ooMu;
  134.     /** μ .*/
  135.     protected double mu;

  136.     // CHECKSTYLE: resume VisibilityModifierCheck

  137.     /** Contribution to be numerically averaged. */
  138.     private final ForceModel contribution;

  139.     /** Gauss integrator. */
  140.     private final double threshold;

  141.     /** Gauss integrator. */
  142.     private GaussQuadrature integrator;

  143.     /** Flag for Gauss order computation. */
  144.     private boolean isDirty;

  145.     /** Attitude provider. */
  146.     private AttitudeProvider attitudeProvider;

  147.     /** Prefix for coefficients keys. */
  148.     private final String coefficientsKeyPrefix;

  149.     /** Short period terms. */
  150.     private GaussianShortPeriodicCoefficients gaussianSPCoefs;

  151.     /** Build a new instance.
  152.      *  @param coefficientsKeyPrefix prefix for coefficients keys
  153.      *  @param threshold tolerance for the choice of the Gauss quadrature order
  154.      *  @param contribution the {@link ForceModel} to be numerically averaged
  155.      */
  156.     protected AbstractGaussianContribution(final String coefficientsKeyPrefix,
  157.                                            final double threshold,
  158.                                            final ForceModel contribution) {
  159.         this.coefficientsKeyPrefix = coefficientsKeyPrefix;
  160.         this.contribution          = contribution;
  161.         this.threshold             = threshold;
  162.         this.integrator            = new GaussQuadrature(GAUSS_ORDER[MAX_ORDER_RANK]);
  163.         this.isDirty               = true;
  164.     }

  165.     /** {@inheritDoc} */
  166.     @Override
  167.     public List<ShortPeriodTerms> initialize(final AuxiliaryElements aux, final boolean meanOnly) {

  168.         final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
  169.         gaussianSPCoefs = new GaussianShortPeriodicCoefficients(coefficientsKeyPrefix,
  170.                                                                 JMAX, INTERPOLATION_POINTS,
  171.                                                                 new TimeSpanMap<Slot>(new Slot(JMAX, INTERPOLATION_POINTS)));
  172.         list.add(gaussianSPCoefs);
  173.         return list;

  174.     }

  175.     /** {@inheritDoc} */
  176.     @Override
  177.     public void initializeStep(final AuxiliaryElements aux)
  178.         throws OrekitException {

  179.         // Equinoctial elements
  180.         a  = aux.getSma();
  181.         k  = aux.getK();
  182.         h  = aux.getH();
  183.         q  = aux.getQ();
  184.         p  = aux.getP();

  185.         // Eccentricity
  186.         ecc = aux.getEcc();

  187.         // Equinoctial coefficients
  188.         A = aux.getA();
  189.         B = aux.getB();
  190.         C = aux.getC();

  191.         // Equinoctial frame vectors
  192.         f = aux.getVectorF();
  193.         g = aux.getVectorG();
  194.         w = aux.getVectorW();

  195.         // Kepler mean motion
  196.         n = aux.getMeanMotion();

  197.         // Mean longitude
  198.         lm = aux.getLM();

  199.         // 1 / A
  200.         ooA = 1. / A;
  201.         // 1 / AB
  202.         ooAB = ooA / B;
  203.         // C / 2AB
  204.         co2AB = C * ooAB / 2.;
  205.         // 1 / (1 + B)
  206.         ooBpo = 1. / (1. + B);
  207.         // 2 / (n² * a)
  208.         ton2a = 2. / (n * n * a);
  209.         // mu
  210.         mu = aux.getMu();
  211.         // 1 / mu
  212.         ooMu  = 1. / mu;
  213.     }

  214.     /** {@inheritDoc} */
  215.     @Override
  216.     public double[] getMeanElementRate(final SpacecraftState state) throws OrekitException {

  217.         double[] meanElementRate = new double[6];
  218.         // Computes the limits for the integral
  219.         final double[] ll = getLLimits(state);
  220.         // Computes integrated mean element rates if Llow < Lhigh
  221.         if (ll[0] < ll[1]) {
  222.             meanElementRate = getMeanElementRate(state, integrator, ll[0], ll[1]);
  223.             if (isDirty) {
  224.                 boolean next = true;
  225.                 for (int i = 0; i < MAX_ORDER_RANK && next; i++) {
  226.                     final double[] meanRates = getMeanElementRate(state, new GaussQuadrature(GAUSS_ORDER[i]), ll[0], ll[1]);
  227.                     if (getRatesDiff(meanElementRate, meanRates) < threshold) {
  228.                         integrator = new GaussQuadrature(GAUSS_ORDER[i]);
  229.                         next = false;
  230.                     }
  231.                 }
  232.                 isDirty = false;
  233.             }
  234.         }
  235.         return meanElementRate;
  236.     }

  237.     /** Compute the limits in L, the true longitude, for integration.
  238.      *
  239.      *  @param  state current state information: date, kinematics, attitude
  240.      *  @return the integration limits in L
  241.      *  @exception OrekitException if some specific error occurs
  242.      */
  243.     protected abstract double[] getLLimits(SpacecraftState state) throws OrekitException;

  244.     /** Computes the mean equinoctial elements rates da<sub>i</sub> / dt.
  245.      *
  246.      *  @param state current state
  247.      *  @param gauss Gauss quadrature
  248.      *  @param low lower bound of the integral interval
  249.      *  @param high upper bound of the integral interval
  250.      *  @return the mean element rates
  251.      *  @throws OrekitException if some specific error occurs
  252.      */
  253.     private double[] getMeanElementRate(final SpacecraftState state,
  254.             final GaussQuadrature gauss,
  255.             final double low,
  256.             final double high) throws OrekitException {
  257.         final double[] meanElementRate = gauss.integrate(new IntegrableFunction(state, true, 0), low, high);
  258.         // Constant multiplier for integral
  259.         final double coef = 1. / (2. * FastMath.PI * B);
  260.         // Corrects mean element rates
  261.         for (int i = 0; i < 6; i++) {
  262.             meanElementRate[i] *= coef;
  263.         }
  264.         return meanElementRate;
  265.     }

  266.     /** Estimates the weighted magnitude of the difference between 2 sets of equinoctial elements rates.
  267.      *
  268.      *  @param meanRef reference rates
  269.      *  @param meanCur current rates
  270.      *  @return estimated magnitude of weighted differences
  271.      */
  272.     private double getRatesDiff(final double[] meanRef, final double[] meanCur) {
  273.         double maxDiff = FastMath.abs(meanRef[0] - meanCur[0]) / a;
  274.         // Corrects mean element rates
  275.         for (int i = 1; i < meanRef.length; i++) {
  276.             final double diff = FastMath.abs(meanRef[i] - meanCur[i]);
  277.             if (maxDiff < diff) maxDiff = diff;
  278.         }
  279.         return maxDiff;
  280.     }

  281.     /** {@inheritDoc} */
  282.     @Override
  283.     public void registerAttitudeProvider(final AttitudeProvider provider) {
  284.         this.attitudeProvider = provider;
  285.     }

  286.     /** {@inheritDoc} */
  287.     @Override
  288.     public void updateShortPeriodTerms(final SpacecraftState... meanStates)
  289.         throws OrekitException {

  290.         final Slot slot = gaussianSPCoefs.createSlot(meanStates);
  291.         for (final SpacecraftState meanState : meanStates) {
  292.             initializeStep(new AuxiliaryElements(meanState.getOrbit(), I));
  293.             final double[][] currentRhoSigmaj = computeRhoSigmaCoefficients(meanState.getDate());
  294.             final FourierCjSjCoefficients fourierCjSj = new FourierCjSjCoefficients(meanState, JMAX);
  295.             final UijVijCoefficients uijvij = new UijVijCoefficients(currentRhoSigmaj, fourierCjSj, JMAX);
  296.             gaussianSPCoefs.computeCoefficients(meanState, slot, fourierCjSj, uijvij, n, a);
  297.         }

  298.     }

  299.     /**
  300.      * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
  301.      * <p>
  302.      * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
  303.      *  ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
  304.      *  σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
  305.      * </p>
  306.      * @param date current date
  307.      * @return computed coefficients
  308.      */
  309.     private double[][] computeRhoSigmaCoefficients(final AbsoluteDate date) {
  310.         final double[][] currentRhoSigmaj = new double[2][3 * JMAX + 1];
  311.         final CjSjCoefficient cjsjKH = new CjSjCoefficient(k, h);
  312.         final double b = 1. / (1 + B);

  313.         // (-b)<sup>j</sup>
  314.         double mbtj = 1;

  315.         for (int j = 1; j <= 3 * JMAX; j++) {

  316.             //Compute current rho and sigma;
  317.             mbtj *= -b;
  318.             final double coef = (1 + j * B) * mbtj;
  319.             currentRhoSigmaj[0][j] = coef * cjsjKH.getCj(j);
  320.             currentRhoSigmaj[1][j] = coef * cjsjKH.getSj(j);
  321.         }
  322.         return currentRhoSigmaj;
  323.     }

  324.     /** Internal class for numerical quadrature. */
  325.     private class IntegrableFunction implements UnivariateVectorFunction {

  326.         /** Current state. */
  327.         private final SpacecraftState state;

  328.         /** Signal that this class is used to compute the values required by the mean element variations
  329.          * or by the short periodic element variations. */
  330.         private final boolean meanMode;

  331.         /** The j index.
  332.          * <p>
  333.          * Used only for short periodic variation. Ignored for mean elements variation.
  334.          * </p> */
  335.         private final int j;

  336.         /** Build a new instance.
  337.          *  @param  state current state information: date, kinematics, attitude
  338.          *  @param meanMode if true return the value associated to the mean elements variation,
  339.          *                  if false return the values associated to the short periodic elements variation
  340.          * @param j the j index. used only for short periodic variation. Ignored for mean elements variation.
  341.          */
  342.         IntegrableFunction(final SpacecraftState state, final boolean meanMode, final int j) {

  343.             // remove derivatives from state
  344.             final double[] stateVector = new double[6];
  345.             OrbitType.EQUINOCTIAL.mapOrbitToArray(state.getOrbit(), PositionAngle.TRUE, stateVector, null);
  346.             final Orbit fixedOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(stateVector, null, PositionAngle.TRUE,
  347.                                                                            state.getDate(),
  348.                                                                            state.getMu(),
  349.                                                                            state.getFrame());
  350.             this.state = new SpacecraftState(fixedOrbit, state.getAttitude(), state.getMass());
  351.             this.meanMode = meanMode;
  352.             this.j = j;
  353.         }

  354.         /** {@inheritDoc} */
  355.         @Override
  356.         public double[] value(final double x) {

  357.             //Compute the time difference from the true longitude difference
  358.             final double shiftedLm = trueToMean(x);
  359.             final double dLm = shiftedLm - lm;
  360.             final double dt = dLm / n;

  361.             final double cosL = FastMath.cos(x);
  362.             final double sinL = FastMath.sin(x);
  363.             final double roa  = B * B / (1. + h * sinL + k * cosL);
  364.             final double roa2 = roa * roa;
  365.             final double r    = a * roa;
  366.             final double X    = r * cosL;
  367.             final double Y    = r * sinL;
  368.             final double naob = n * a / B;
  369.             final double Xdot = -naob * (h + sinL);
  370.             final double Ydot =  naob * (k + cosL);
  371.             final Vector3D vel = new Vector3D(Xdot, f, Ydot, g);

  372.             // Compute acceleration
  373.             Vector3D acc = Vector3D.ZERO;
  374.             try {

  375.                 // shift the orbit to dt
  376.                 final Orbit shiftedOrbit = state.getOrbit().shiftedBy(dt);

  377.                 // Recompose an orbit with time held fixed to be compliant with DSST theory
  378.                 final Orbit recomposedOrbit =
  379.                         new EquinoctialOrbit(shiftedOrbit.getA(),
  380.                                              shiftedOrbit.getEquinoctialEx(),
  381.                                              shiftedOrbit.getEquinoctialEy(),
  382.                                              shiftedOrbit.getHx(),
  383.                                              shiftedOrbit.getHy(),
  384.                                              shiftedOrbit.getLv(),
  385.                                              PositionAngle.TRUE,
  386.                                              shiftedOrbit.getFrame(),
  387.                                              state.getDate(),
  388.                                              shiftedOrbit.getMu());

  389.                 // Get the corresponding attitude
  390.                 final Attitude recomposedAttitude =
  391.                         attitudeProvider.getAttitude(recomposedOrbit,
  392.                                                      recomposedOrbit.getDate(),
  393.                                                      recomposedOrbit.getFrame());

  394.                 // create shifted SpacecraftState with attitude at specified time
  395.                 final SpacecraftState shiftedState =
  396.                         new SpacecraftState(recomposedOrbit, recomposedAttitude, state.getMass());

  397.                 acc = contribution.acceleration(shiftedState, contribution.getParameters());

  398.             } catch (OrekitException oe) {
  399.                 throw new OrekitExceptionWrapper(oe);
  400.             }
  401.             //Compute the derivatives of the elements by the speed
  402.             final double[] deriv = new double[6];
  403.             // da/dv
  404.             deriv[0] = getAoV(vel).dotProduct(acc);
  405.             // dex/dv
  406.             deriv[1] = getKoV(X, Y, Xdot, Ydot).dotProduct(acc);
  407.             // dey/dv
  408.             deriv[2] = getHoV(X, Y, Xdot, Ydot).dotProduct(acc);
  409.             // dhx/dv
  410.             deriv[3] = getQoV(X).dotProduct(acc);
  411.             // dhy/dv
  412.             deriv[4] = getPoV(Y).dotProduct(acc);
  413.             // dλ/dv
  414.             deriv[5] = getLoV(X, Y, Xdot, Ydot).dotProduct(acc);

  415.             // Compute mean elements rates
  416.             double[] val = null;
  417.             if (meanMode) {
  418.                 val = new double[6];
  419.                 for (int i = 0; i < 6; i++) {
  420.                     // da<sub>i</sub>/dt
  421.                     val[i] = roa2 * deriv[i];
  422.                 }
  423.             } else {
  424.                 val = new double[12];
  425.                 //Compute cos(j*L) and sin(j*L);
  426.                 final double cosjL = j == 1 ? cosL : FastMath.cos(j * x);
  427.                 final double sinjL = j == 1 ? sinL : FastMath.sin(j * x);

  428.                 for (int i = 0; i < 6; i++) {
  429.                     // da<sub>i</sub>/dv * cos(jL)
  430.                     val[i] = cosjL * deriv[i];
  431.                     // da<sub>i</sub>/dv * sin(jL)
  432.                     val[i + 6] = sinjL * deriv[i];
  433.                 }
  434.             }
  435.             return val;
  436.         }

  437.         /** Converts true longitude to eccentric longitude.
  438.          * @param lv True longitude
  439.          * @return Eccentric longitude
  440.          */
  441.         private double trueToEccentric (final double lv) {
  442.             final double cosLv   = FastMath.cos(lv);
  443.             final double sinLv   = FastMath.sin(lv);
  444.             final double num     = h * cosLv - k * sinLv;
  445.             final double den     = B + 1 + k * cosLv + h * sinLv;
  446.             return lv + 2 * FastMath.atan(num / den);
  447.         }

  448.         /** Converts eccentric longitude to mean longitude.
  449.          * @param le Eccentric longitude
  450.          * @return Mean longitude
  451.          */
  452.         private double eccentricToMean (final double le) {
  453.             return le - k * FastMath.sin(le) + h * FastMath.cos(le);
  454.         }

  455.         /** Converts true longitude to mean longitude.
  456.          * @param lv True longitude
  457.          * @return Eccentric longitude
  458.          */
  459.         private double trueToMean (final double lv) {
  460.             return eccentricToMean(trueToEccentric(lv));
  461.         }

  462.         /** Compute δa/δv.
  463.          *  @param vel satellite velocity
  464.          *  @return δa/δv
  465.          */
  466.         private Vector3D getAoV(final Vector3D vel) {
  467.             return new Vector3D(ton2a, vel);
  468.         }

  469.         /** Compute δh/δv.
  470.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  471.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  472.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  473.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  474.          *  @return δh/δv
  475.          */
  476.         private Vector3D getHoV(final double X, final double Y, final double Xdot, final double Ydot) {
  477.             final double kf = (2. * Xdot * Y - X * Ydot) * ooMu;
  478.             final double kg = X * Xdot * ooMu;
  479.             final double kw = k * (I * q * Y - p * X) * ooAB;
  480.             return new Vector3D(kf, f, -kg, g, kw, w);
  481.         }

  482.         /** Compute δk/δv.
  483.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  484.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  485.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  486.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  487.          *  @return δk/δv
  488.          */
  489.         private Vector3D getKoV(final double X, final double Y, final double Xdot, final double Ydot) {
  490.             final double kf = Y * Ydot * ooMu;
  491.             final double kg = (2. * X * Ydot - Xdot * Y) * ooMu;
  492.             final double kw = h * (I * q * Y - p * X) * ooAB;
  493.             return new Vector3D(-kf, f, kg, g, -kw, w);
  494.         }

  495.         /** Compute δp/δv.
  496.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  497.          *  @return δp/δv
  498.          */
  499.         private Vector3D getPoV(final double Y) {
  500.             return new Vector3D(co2AB * Y, w);
  501.         }

  502.         /** Compute δq/δv.
  503.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  504.          *  @return δq/δv
  505.          */
  506.         private Vector3D getQoV(final double X) {
  507.             return new Vector3D(I * co2AB * X, w);
  508.         }

  509.         /** Compute δλ/δv.
  510.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  511.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  512.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  513.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  514.          *  @return δλ/δv
  515.          */
  516.         private Vector3D getLoV(final double X, final double Y, final double Xdot, final double Ydot) {
  517.             final Vector3D pos = new Vector3D(X, f, Y, g);
  518.             final Vector3D v2  = new Vector3D(k, getHoV(X, Y, Xdot, Ydot), -h, getKoV(X, Y, Xdot, Ydot));
  519.             return new Vector3D(-2. * ooA, pos, ooBpo, v2, (I * q * Y - p * X) * ooA, w);
  520.         }

  521.     }

  522.     /** Class used to {@link #integrate(UnivariateVectorFunction, double, double) integrate}
  523.      *  a {@link org.hipparchus.analysis.UnivariateVectorFunction function}
  524.      *  of the orbital elements using the Gaussian quadrature rule to get the acceleration.
  525.      */
  526.     private static class GaussQuadrature {

  527.         // CHECKSTYLE: stop NoWhitespaceAfter

  528.         // Points and weights for the available quadrature orders

  529.         /** Points for quadrature of order 12. */
  530.         private static final double[] P_12 = {
  531.             -0.98156063424671910000,
  532.             -0.90411725637047490000,
  533.             -0.76990267419430470000,
  534.             -0.58731795428661740000,
  535.             -0.36783149899818024000,
  536.             -0.12523340851146890000,
  537.             0.12523340851146890000,
  538.             0.36783149899818024000,
  539.             0.58731795428661740000,
  540.             0.76990267419430470000,
  541.             0.90411725637047490000,
  542.             0.98156063424671910000
  543.         };

  544.         /** Weights for quadrature of order 12. */
  545.         private static final double[] W_12 = {
  546.             0.04717533638651220000,
  547.             0.10693932599531830000,
  548.             0.16007832854334633000,
  549.             0.20316742672306584000,
  550.             0.23349253653835478000,
  551.             0.24914704581340286000,
  552.             0.24914704581340286000,
  553.             0.23349253653835478000,
  554.             0.20316742672306584000,
  555.             0.16007832854334633000,
  556.             0.10693932599531830000,
  557.             0.04717533638651220000
  558.         };

  559.         /** Points for quadrature of order 16. */
  560.         private static final double[] P_16 = {
  561.             -0.98940093499164990000,
  562.             -0.94457502307323260000,
  563.             -0.86563120238783160000,
  564.             -0.75540440835500310000,
  565.             -0.61787624440264380000,
  566.             -0.45801677765722737000,
  567.             -0.28160355077925890000,
  568.             -0.09501250983763745000,
  569.             0.09501250983763745000,
  570.             0.28160355077925890000,
  571.             0.45801677765722737000,
  572.             0.61787624440264380000,
  573.             0.75540440835500310000,
  574.             0.86563120238783160000,
  575.             0.94457502307323260000,
  576.             0.98940093499164990000
  577.         };

  578.         /** Weights for quadrature of order 16. */
  579.         private static final double[] W_16 = {
  580.             0.02715245941175405800,
  581.             0.06225352393864777000,
  582.             0.09515851168249283000,
  583.             0.12462897125553388000,
  584.             0.14959598881657685000,
  585.             0.16915651939500256000,
  586.             0.18260341504492360000,
  587.             0.18945061045506847000,
  588.             0.18945061045506847000,
  589.             0.18260341504492360000,
  590.             0.16915651939500256000,
  591.             0.14959598881657685000,
  592.             0.12462897125553388000,
  593.             0.09515851168249283000,
  594.             0.06225352393864777000,
  595.             0.02715245941175405800
  596.         };

  597.         /** Points for quadrature of order 20. */
  598.         private static final double[] P_20 = {
  599.             -0.99312859918509490000,
  600.             -0.96397192727791390000,
  601.             -0.91223442825132600000,
  602.             -0.83911697182221890000,
  603.             -0.74633190646015080000,
  604.             -0.63605368072651510000,
  605.             -0.51086700195082700000,
  606.             -0.37370608871541955000,
  607.             -0.22778585114164507000,
  608.             -0.07652652113349734000,
  609.             0.07652652113349734000,
  610.             0.22778585114164507000,
  611.             0.37370608871541955000,
  612.             0.51086700195082700000,
  613.             0.63605368072651510000,
  614.             0.74633190646015080000,
  615.             0.83911697182221890000,
  616.             0.91223442825132600000,
  617.             0.96397192727791390000,
  618.             0.99312859918509490000
  619.         };

  620.         /** Weights for quadrature of order 20. */
  621.         private static final double[] W_20 = {
  622.             0.01761400713915226400,
  623.             0.04060142980038684000,
  624.             0.06267204833410904000,
  625.             0.08327674157670477000,
  626.             0.10193011981724048000,
  627.             0.11819453196151844000,
  628.             0.13168863844917678000,
  629.             0.14209610931838212000,
  630.             0.14917298647260380000,
  631.             0.15275338713072600000,
  632.             0.15275338713072600000,
  633.             0.14917298647260380000,
  634.             0.14209610931838212000,
  635.             0.13168863844917678000,
  636.             0.11819453196151844000,
  637.             0.10193011981724048000,
  638.             0.08327674157670477000,
  639.             0.06267204833410904000,
  640.             0.04060142980038684000,
  641.             0.01761400713915226400
  642.         };

  643.         /** Points for quadrature of order 24. */
  644.         private static final double[] P_24 = {
  645.             -0.99518721999702130000,
  646.             -0.97472855597130950000,
  647.             -0.93827455200273270000,
  648.             -0.88641552700440100000,
  649.             -0.82000198597390300000,
  650.             -0.74012419157855440000,
  651.             -0.64809365193697550000,
  652.             -0.54542147138883950000,
  653.             -0.43379350762604520000,
  654.             -0.31504267969616340000,
  655.             -0.19111886747361634000,
  656.             -0.06405689286260563000,
  657.             0.06405689286260563000,
  658.             0.19111886747361634000,
  659.             0.31504267969616340000,
  660.             0.43379350762604520000,
  661.             0.54542147138883950000,
  662.             0.64809365193697550000,
  663.             0.74012419157855440000,
  664.             0.82000198597390300000,
  665.             0.88641552700440100000,
  666.             0.93827455200273270000,
  667.             0.97472855597130950000,
  668.             0.99518721999702130000
  669.         };

  670.         /** Weights for quadrature of order 24. */
  671.         private static final double[] W_24 = {
  672.             0.01234122979998733500,
  673.             0.02853138862893380600,
  674.             0.04427743881741981000,
  675.             0.05929858491543691500,
  676.             0.07334648141108027000,
  677.             0.08619016153195320000,
  678.             0.09761865210411391000,
  679.             0.10744427011596558000,
  680.             0.11550566805372553000,
  681.             0.12167047292780335000,
  682.             0.12583745634682825000,
  683.             0.12793819534675221000,
  684.             0.12793819534675221000,
  685.             0.12583745634682825000,
  686.             0.12167047292780335000,
  687.             0.11550566805372553000,
  688.             0.10744427011596558000,
  689.             0.09761865210411391000,
  690.             0.08619016153195320000,
  691.             0.07334648141108027000,
  692.             0.05929858491543691500,
  693.             0.04427743881741981000,
  694.             0.02853138862893380600,
  695.             0.01234122979998733500
  696.         };

  697.         /** Points for quadrature of order 32. */
  698.         private static final double[] P_32 = {
  699.             -0.99726386184948160000,
  700.             -0.98561151154526840000,
  701.             -0.96476225558750640000,
  702.             -0.93490607593773970000,
  703.             -0.89632115576605220000,
  704.             -0.84936761373256990000,
  705.             -0.79448379596794250000,
  706.             -0.73218211874028970000,
  707.             -0.66304426693021520000,
  708.             -0.58771575724076230000,
  709.             -0.50689990893222950000,
  710.             -0.42135127613063540000,
  711.             -0.33186860228212767000,
  712.             -0.23928736225213710000,
  713.             -0.14447196158279646000,
  714.             -0.04830766568773831000,
  715.             0.04830766568773831000,
  716.             0.14447196158279646000,
  717.             0.23928736225213710000,
  718.             0.33186860228212767000,
  719.             0.42135127613063540000,
  720.             0.50689990893222950000,
  721.             0.58771575724076230000,
  722.             0.66304426693021520000,
  723.             0.73218211874028970000,
  724.             0.79448379596794250000,
  725.             0.84936761373256990000,
  726.             0.89632115576605220000,
  727.             0.93490607593773970000,
  728.             0.96476225558750640000,
  729.             0.98561151154526840000,
  730.             0.99726386184948160000
  731.         };

  732.         /** Weights for quadrature of order 32. */
  733.         private static final double[] W_32 = {
  734.             0.00701861000947013600,
  735.             0.01627439473090571200,
  736.             0.02539206530926214200,
  737.             0.03427386291302141000,
  738.             0.04283589802222658600,
  739.             0.05099805926237621600,
  740.             0.05868409347853559000,
  741.             0.06582222277636193000,
  742.             0.07234579410884862000,
  743.             0.07819389578707042000,
  744.             0.08331192422694673000,
  745.             0.08765209300440380000,
  746.             0.09117387869576390000,
  747.             0.09384439908080441000,
  748.             0.09563872007927487000,
  749.             0.09654008851472784000,
  750.             0.09654008851472784000,
  751.             0.09563872007927487000,
  752.             0.09384439908080441000,
  753.             0.09117387869576390000,
  754.             0.08765209300440380000,
  755.             0.08331192422694673000,
  756.             0.07819389578707042000,
  757.             0.07234579410884862000,
  758.             0.06582222277636193000,
  759.             0.05868409347853559000,
  760.             0.05099805926237621600,
  761.             0.04283589802222658600,
  762.             0.03427386291302141000,
  763.             0.02539206530926214200,
  764.             0.01627439473090571200,
  765.             0.00701861000947013600
  766.         };

  767.         /** Points for quadrature of order 40. */
  768.         private static final double[] P_40 = {
  769.             -0.99823770971055930000,
  770.             -0.99072623869945710000,
  771.             -0.97725994998377420000,
  772.             -0.95791681921379170000,
  773.             -0.93281280827867660000,
  774.             -0.90209880696887420000,
  775.             -0.86595950321225960000,
  776.             -0.82461223083331170000,
  777.             -0.77830565142651940000,
  778.             -0.72731825518992710000,
  779.             -0.67195668461417960000,
  780.             -0.61255388966798030000,
  781.             -0.54946712509512820000,
  782.             -0.48307580168617870000,
  783.             -0.41377920437160500000,
  784.             -0.34199409082575850000,
  785.             -0.26815218500725370000,
  786.             -0.19269758070137110000,
  787.             -0.11608407067525522000,
  788.             -0.03877241750605081600,
  789.             0.03877241750605081600,
  790.             0.11608407067525522000,
  791.             0.19269758070137110000,
  792.             0.26815218500725370000,
  793.             0.34199409082575850000,
  794.             0.41377920437160500000,
  795.             0.48307580168617870000,
  796.             0.54946712509512820000,
  797.             0.61255388966798030000,
  798.             0.67195668461417960000,
  799.             0.72731825518992710000,
  800.             0.77830565142651940000,
  801.             0.82461223083331170000,
  802.             0.86595950321225960000,
  803.             0.90209880696887420000,
  804.             0.93281280827867660000,
  805.             0.95791681921379170000,
  806.             0.97725994998377420000,
  807.             0.99072623869945710000,
  808.             0.99823770971055930000
  809.         };

  810.         /** Weights for quadrature of order 40. */
  811.         private static final double[] W_40 = {
  812.             0.00452127709853309800,
  813.             0.01049828453115270400,
  814.             0.01642105838190797300,
  815.             0.02224584919416689000,
  816.             0.02793700698002338000,
  817.             0.03346019528254786500,
  818.             0.03878216797447199000,
  819.             0.04387090818567333000,
  820.             0.04869580763507221000,
  821.             0.05322784698393679000,
  822.             0.05743976909939157000,
  823.             0.06130624249292891000,
  824.             0.06480401345660108000,
  825.             0.06791204581523394000,
  826.             0.07061164739128681000,
  827.             0.07288658239580408000,
  828.             0.07472316905796833000,
  829.             0.07611036190062619000,
  830.             0.07703981816424793000,
  831.             0.07750594797842482000,
  832.             0.07750594797842482000,
  833.             0.07703981816424793000,
  834.             0.07611036190062619000,
  835.             0.07472316905796833000,
  836.             0.07288658239580408000,
  837.             0.07061164739128681000,
  838.             0.06791204581523394000,
  839.             0.06480401345660108000,
  840.             0.06130624249292891000,
  841.             0.05743976909939157000,
  842.             0.05322784698393679000,
  843.             0.04869580763507221000,
  844.             0.04387090818567333000,
  845.             0.03878216797447199000,
  846.             0.03346019528254786500,
  847.             0.02793700698002338000,
  848.             0.02224584919416689000,
  849.             0.01642105838190797300,
  850.             0.01049828453115270400,
  851.             0.00452127709853309800
  852.         };

  853.         /** Points for quadrature of order 48. */
  854.         private static final double[] P_48 = {
  855.             -0.99877100725242610000,
  856.             -0.99353017226635080000,
  857.             -0.98412458372282700000,
  858.             -0.97059159254624720000,
  859.             -0.95298770316043080000,
  860.             -0.93138669070655440000,
  861.             -0.90587913671556960000,
  862.             -0.87657202027424800000,
  863.             -0.84358826162439350000,
  864.             -0.80706620402944250000,
  865.             -0.76715903251574020000,
  866.             -0.72403413092381470000,
  867.             -0.67787237963266400000,
  868.             -0.62886739677651370000,
  869.             -0.57722472608397270000,
  870.             -0.52316097472223300000,
  871.             -0.46690290475095840000,
  872.             -0.40868648199071680000,
  873.             -0.34875588629216070000,
  874.             -0.28736248735545555000,
  875.             -0.22476379039468908000,
  876.             -0.16122235606889174000,
  877.             -0.09700469920946270000,
  878.             -0.03238017096286937000,
  879.             0.03238017096286937000,
  880.             0.09700469920946270000,
  881.             0.16122235606889174000,
  882.             0.22476379039468908000,
  883.             0.28736248735545555000,
  884.             0.34875588629216070000,
  885.             0.40868648199071680000,
  886.             0.46690290475095840000,
  887.             0.52316097472223300000,
  888.             0.57722472608397270000,
  889.             0.62886739677651370000,
  890.             0.67787237963266400000,
  891.             0.72403413092381470000,
  892.             0.76715903251574020000,
  893.             0.80706620402944250000,
  894.             0.84358826162439350000,
  895.             0.87657202027424800000,
  896.             0.90587913671556960000,
  897.             0.93138669070655440000,
  898.             0.95298770316043080000,
  899.             0.97059159254624720000,
  900.             0.98412458372282700000,
  901.             0.99353017226635080000,
  902.             0.99877100725242610000
  903.         };

  904.         /** Weights for quadrature of order 48. */
  905.         private static final double[] W_48 = {
  906.             0.00315334605230596250,
  907.             0.00732755390127620800,
  908.             0.01147723457923446900,
  909.             0.01557931572294386600,
  910.             0.01961616045735556700,
  911.             0.02357076083932435600,
  912.             0.02742650970835688000,
  913.             0.03116722783279807000,
  914.             0.03477722256477045000,
  915.             0.03824135106583080600,
  916.             0.04154508294346483000,
  917.             0.04467456085669424000,
  918.             0.04761665849249054000,
  919.             0.05035903555385448000,
  920.             0.05289018948519365000,
  921.             0.05519950369998416500,
  922.             0.05727729210040315000,
  923.             0.05911483969839566000,
  924.             0.06070443916589384000,
  925.             0.06203942315989268000,
  926.             0.06311419228625403000,
  927.             0.06392423858464817000,
  928.             0.06446616443595010000,
  929.             0.06473769681268386000,
  930.             0.06473769681268386000,
  931.             0.06446616443595010000,
  932.             0.06392423858464817000,
  933.             0.06311419228625403000,
  934.             0.06203942315989268000,
  935.             0.06070443916589384000,
  936.             0.05911483969839566000,
  937.             0.05727729210040315000,
  938.             0.05519950369998416500,
  939.             0.05289018948519365000,
  940.             0.05035903555385448000,
  941.             0.04761665849249054000,
  942.             0.04467456085669424000,
  943.             0.04154508294346483000,
  944.             0.03824135106583080600,
  945.             0.03477722256477045000,
  946.             0.03116722783279807000,
  947.             0.02742650970835688000,
  948.             0.02357076083932435600,
  949.             0.01961616045735556700,
  950.             0.01557931572294386600,
  951.             0.01147723457923446900,
  952.             0.00732755390127620800,
  953.             0.00315334605230596250
  954.         };
  955.         // CHECKSTYLE: resume NoWhitespaceAfter

  956.         /** Node points. */
  957.         private final double[] nodePoints;

  958.         /** Node weights. */
  959.         private final double[] nodeWeights;

  960.         /** Creates a Gauss integrator of the given order.
  961.          *
  962.          *  @param numberOfPoints Order of the integration rule.
  963.          */
  964.         GaussQuadrature(final int numberOfPoints) {

  965.             switch(numberOfPoints) {
  966.                 case 12 :
  967.                     this.nodePoints  = P_12.clone();
  968.                     this.nodeWeights = W_12.clone();
  969.                     break;
  970.                 case 16 :
  971.                     this.nodePoints  = P_16.clone();
  972.                     this.nodeWeights = W_16.clone();
  973.                     break;
  974.                 case 20 :
  975.                     this.nodePoints  = P_20.clone();
  976.                     this.nodeWeights = W_20.clone();
  977.                     break;
  978.                 case 24 :
  979.                     this.nodePoints  = P_24.clone();
  980.                     this.nodeWeights = W_24.clone();
  981.                     break;
  982.                 case 32 :
  983.                     this.nodePoints  = P_32.clone();
  984.                     this.nodeWeights = W_32.clone();
  985.                     break;
  986.                 case 40 :
  987.                     this.nodePoints  = P_40.clone();
  988.                     this.nodeWeights = W_40.clone();
  989.                     break;
  990.                 case 48 :
  991.                 default :
  992.                     this.nodePoints  = P_48.clone();
  993.                     this.nodeWeights = W_48.clone();
  994.                     break;
  995.             }

  996.         }

  997.         /** Integrates a given function on the given interval.
  998.          *
  999.          *  @param f Function to integrate.
  1000.          *  @param lowerBound Lower bound of the integration interval.
  1001.          *  @param upperBound Upper bound of the integration interval.
  1002.          *  @return the integral of the weighted function.
  1003.          */
  1004.         public double[] integrate(final UnivariateVectorFunction f,
  1005.                 final double lowerBound, final double upperBound) {

  1006.             final double[] adaptedPoints  = nodePoints.clone();
  1007.             final double[] adaptedWeights = nodeWeights.clone();
  1008.             transform(adaptedPoints, adaptedWeights, lowerBound, upperBound);
  1009.             return basicIntegrate(f, adaptedPoints, adaptedWeights);
  1010.         }

  1011.         /** Performs a change of variable so that the integration
  1012.          *  can be performed on an arbitrary interval {@code [a, b]}.
  1013.          *  <p>
  1014.          *  It is assumed that the natural interval is {@code [-1, 1]}.
  1015.          *  </p>
  1016.          *
  1017.          * @param points  Points to adapt to the new interval.
  1018.          * @param weights Weights to adapt to the new interval.
  1019.          * @param a Lower bound of the integration interval.
  1020.          * @param b Lower bound of the integration interval.
  1021.          */
  1022.         private void transform(final double[] points, final double[] weights,
  1023.                 final double a, final double b) {
  1024.             // Scaling
  1025.             final double scale = (b - a) / 2;
  1026.             final double shift = a + scale;
  1027.             for (int i = 0; i < points.length; i++) {
  1028.                 points[i]   = points[i] * scale + shift;
  1029.                 weights[i] *= scale;
  1030.             }
  1031.         }

  1032.         /** Returns an estimate of the integral of {@code f(x) * w(x)},
  1033.          *  where {@code w} is a weight function that depends on the actual
  1034.          *  flavor of the Gauss integration scheme.
  1035.          *
  1036.          * @param f Function to integrate.
  1037.          * @param points  Nodes.
  1038.          * @param weights Nodes weights.
  1039.          * @return the integral of the weighted function.
  1040.          */
  1041.         private double[] basicIntegrate(final UnivariateVectorFunction f,
  1042.                 final double[] points,
  1043.                 final double[] weights) {
  1044.             double x = points[0];
  1045.             double w = weights[0];
  1046.             double[] v = f.value(x);
  1047.             final double[] y = new double[v.length];
  1048.             for (int j = 0; j < v.length; j++) {
  1049.                 y[j] = w * v[j];
  1050.             }
  1051.             final double[] t = y.clone();
  1052.             final double[] c = new double[v.length];
  1053.             final double[] s = t.clone();
  1054.             for (int i = 1; i < points.length; i++) {
  1055.                 x = points[i];
  1056.                 w = weights[i];
  1057.                 v = f.value(x);
  1058.                 for (int j = 0; j < v.length; j++) {
  1059.                     y[j] = w * v[j] - c[j];
  1060.                     t[j] =  s[j] + y[j];
  1061.                     c[j] = (t[j] - s[j]) - y[j];
  1062.                     s[j] = t[j];
  1063.                 }
  1064.             }
  1065.             return s;
  1066.         }

  1067.     }

  1068.     /** Compute the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
  1069.      *  <p>
  1070.      *  Those coefficients are given in Danielson paper by expression 4.4-(6)
  1071.      *  </p>
  1072.      *  @author Petre Bazavan
  1073.      *  @author Lucian Barbulescu
  1074.      */
  1075.     private class FourierCjSjCoefficients {

  1076.         /** Maximum possible value for j. */
  1077.         private final int jMax;

  1078.         /** The C<sub>i</sub><sup>j</sup> coefficients.
  1079.          * <p>
  1080.          * the index i corresponds to the following elements: <br/>
  1081.          * - 0 for a <br>
  1082.          * - 1 for k <br>
  1083.          * - 2 for h <br>
  1084.          * - 3 for q <br>
  1085.          * - 4 for p <br>
  1086.          * - 5 for λ <br>
  1087.          * </p>
  1088.          */
  1089.         private final double[][] cCoef;

  1090.         /** The C<sub>i</sub><sup>j</sup> coefficients.
  1091.          * <p>
  1092.          * the index i corresponds to the following elements: <br/>
  1093.          * - 0 for a <br>
  1094.          * - 1 for k <br>
  1095.          * - 2 for h <br>
  1096.          * - 3 for q <br>
  1097.          * - 4 for p <br>
  1098.          * - 5 for λ <br>
  1099.          * </p>
  1100.          */
  1101.         private final double[][] sCoef;

  1102.         /** Standard constructor.
  1103.          * @param state the current state
  1104.          * @param jMax maximum value for j
  1105.          * @throws OrekitException in case of an error
  1106.          */
  1107.         FourierCjSjCoefficients(final SpacecraftState state, final int jMax)
  1108.             throws OrekitException {
  1109.             //Initialise the fields
  1110.             this.jMax = jMax;

  1111.             //Allocate the arrays
  1112.             final int rows = jMax + 1;
  1113.             cCoef = new double[rows][6];
  1114.             sCoef = new double[rows][6];

  1115.             //Compute the coefficients
  1116.             computeCoefficients(state);
  1117.         }

  1118.         /**
  1119.          * Compute the Fourrier coefficients.
  1120.          * <p>
  1121.          * Only the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients need to be computed
  1122.          * as D<sub>i</sub><sup>m</sup> is always 0.
  1123.          * </p>
  1124.          * @param state the current state
  1125.          * @throws OrekitException in case of an error
  1126.          */
  1127.         private void computeCoefficients(final SpacecraftState state)
  1128.             throws OrekitException {
  1129.             // Computes the limits for the integral
  1130.             final double[] ll = getLLimits(state);
  1131.             // Computes integrated mean element rates if Llow < Lhigh
  1132.             if (ll[0] < ll[1]) {
  1133.                 //Compute 1 / PI
  1134.                 final double ooPI = 1 / FastMath.PI;

  1135.                 // loop through all values of j
  1136.                 for (int j = 0; j <= jMax; j++) {
  1137.                     final double[] curentCoefficients =
  1138.                             integrator.integrate(new IntegrableFunction(state, false, j), ll[0], ll[1]);

  1139.                     //divide by PI and set the values for the coefficients
  1140.                     for (int i = 0; i < 6; i++) {
  1141.                         cCoef[j][i] = ooPI * curentCoefficients[i];
  1142.                         sCoef[j][i] = ooPI * curentCoefficients[i + 6];
  1143.                     }
  1144.                 }
  1145.             }
  1146.         }

  1147.         /** Get the coefficient C<sub>i</sub><sup>j</sup>.
  1148.          * @param i i index - corresponds to the required variation
  1149.          * @param j j index
  1150.          * @return the coefficient C<sub>i</sub><sup>j</sup>
  1151.          */
  1152.         public double getCij(final int i, final int j) {
  1153.             return cCoef[j][i];
  1154.         }

  1155.         /** Get the coefficient S<sub>i</sub><sup>j</sup>.
  1156.          * @param i i index - corresponds to the required variation
  1157.          * @param j j index
  1158.          * @return the coefficient S<sub>i</sub><sup>j</sup>
  1159.          */
  1160.         public double getSij(final int i, final int j) {
  1161.             return sCoef[j][i];
  1162.         }
  1163.     }

  1164.     /** This class handles the short periodic coefficients described in Danielson 2.5.3-26.
  1165.      *
  1166.      * <p>
  1167.      * The value of M is 0. Also, since the values of the Fourier coefficient D<sub>i</sub><sup>m</sup> is 0
  1168.      * then the values of the coefficients D<sub>i</sub><sup>m</sup> for m &gt; 2 are also 0.
  1169.      * </p>
  1170.      * @author Petre Bazavan
  1171.      * @author Lucian Barbulescu
  1172.      *
  1173.      */
  1174.     private static class GaussianShortPeriodicCoefficients implements ShortPeriodTerms {

  1175.         /** Serializable UID. */
  1176.         private static final long serialVersionUID = 20151118L;

  1177.         /** Maximum value for j index. */
  1178.         private final int jMax;

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

  1181.         /** Prefix for coefficients keys. */
  1182.         private final String coefficientsKeyPrefix;

  1183.         /** All coefficients slots. */
  1184.         private final transient TimeSpanMap<Slot> slots;

  1185.         /** Constructor.
  1186.          *  @param coefficientsKeyPrefix prefix for coefficients keys
  1187.          *  @param jMax maximum value for j index
  1188.          *  @param interpolationPoints number of points used in the interpolation process
  1189.          *  @param slots all coefficients slots
  1190.          */
  1191.         GaussianShortPeriodicCoefficients(final String coefficientsKeyPrefix,
  1192.                                           final int jMax, final int interpolationPoints,
  1193.                                           final TimeSpanMap<Slot> slots) {
  1194.             //Initialize fields
  1195.             this.jMax                  = jMax;
  1196.             this.interpolationPoints   = interpolationPoints;
  1197.             this.coefficientsKeyPrefix = coefficientsKeyPrefix;
  1198.             this.slots                 = slots;
  1199.         }

  1200.         /** Get the slot valid for some date.
  1201.          * @param meanStates mean states defining the slot
  1202.          * @return slot valid at the specified date
  1203.          */
  1204.         public Slot createSlot(final SpacecraftState... meanStates) {
  1205.             final Slot         slot  = new Slot(jMax, interpolationPoints);
  1206.             final AbsoluteDate first = meanStates[0].getDate();
  1207.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  1208.             if (first.compareTo(last) <= 0) {
  1209.                 slots.addValidAfter(slot, first);
  1210.             } else {
  1211.                 slots.addValidBefore(slot, first);
  1212.             }
  1213.             return slot;
  1214.         }

  1215.         /** Compute the short periodic coefficients.
  1216.          *
  1217.          * @param state current state information: date, kinematics, attitude
  1218.          * @param slot coefficients slot
  1219.          * @param fourierCjSj Fourier coefficients
  1220.          * @param uijvij U and V coefficients
  1221.          * @param n Keplerian mean motion
  1222.          * @param a semi major axis
  1223.          * @throws OrekitException if an error occurs
  1224.          */
  1225.         private void computeCoefficients(final SpacecraftState state, final Slot slot,
  1226.                                          final FourierCjSjCoefficients fourierCjSj,
  1227.                                          final UijVijCoefficients uijvij,
  1228.                                          final double n, final double a)
  1229.             throws OrekitException {

  1230.             // get the current date
  1231.             final AbsoluteDate date = state.getDate();

  1232.             // compute the k₂⁰ coefficient
  1233.             final double k20 = computeK20(jMax, uijvij.currentRhoSigmaj);

  1234.             // 1. / n
  1235.             final double oon = 1. / n;
  1236.             // 3. / (2 * a * n)
  1237.             final double to2an = 1.5 * oon / a;
  1238.             // 3. / (4 * a * n)
  1239.             final double to4an = to2an / 2;

  1240.             // Compute the coefficients for each element
  1241.             final int size = jMax + 1;
  1242.             final double[]   di1        = new double[6];
  1243.             final double[]   di2        = new double[6];
  1244.             final double[][] currentCij = new double[size][6];
  1245.             final double[][] currentSij = new double[size][6];
  1246.             for (int i = 0; i < 6; i++) {

  1247.                 // compute D<sub>i</sub>¹ and D<sub>i</sub>² (all others are 0)
  1248.                 di1[i] = -oon * fourierCjSj.getCij(i, 0);
  1249.                 if (i == 5) {
  1250.                     di1[i] += to2an * uijvij.getU1(0, 0);
  1251.                 }
  1252.                 di2[i] = 0.;
  1253.                 if (i == 5) {
  1254.                     di2[i] += -to4an * fourierCjSj.getCij(0, 0);
  1255.                 }

  1256.                 //the C<sub>i</sub>⁰ is computed based on all others
  1257.                 currentCij[0][i] = -di2[i] * k20;

  1258.                 for (int j = 1; j <= jMax; j++) {
  1259.                     // compute the current C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup>
  1260.                     currentCij[j][i] = oon * uijvij.getU1(j, i);
  1261.                     if (i == 5) {
  1262.                         currentCij[j][i] += -to2an * uijvij.getU2(j);
  1263.                     }
  1264.                     currentSij[j][i] = oon * uijvij.getV1(j, i);
  1265.                     if (i == 5) {
  1266.                         currentSij[j][i] += -to2an * uijvij.getV2(j);
  1267.                     }

  1268.                     // add the computed coefficients to C<sub>i</sub>⁰
  1269.                     currentCij[0][i] += -(currentCij[j][i] * uijvij.currentRhoSigmaj[0][j] + currentSij[j][i] * uijvij.currentRhoSigmaj[1][j]);
  1270.                 }

  1271.             }

  1272.             // add the values to the interpolators
  1273.             slot.cij[0].addGridPoint(date, currentCij[0]);
  1274.             slot.dij[1].addGridPoint(date, di1);
  1275.             slot.dij[2].addGridPoint(date, di2);
  1276.             for (int j = 1; j <= jMax; j++) {
  1277.                 slot.cij[j].addGridPoint(date, currentCij[j]);
  1278.                 slot.sij[j].addGridPoint(date, currentSij[j]);
  1279.             }

  1280.         }

  1281.         /** Compute the coefficient k₂⁰ by using the equation
  1282.          * 2.5.3-(9a) from Danielson.
  1283.          * <p>
  1284.          * After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:<br>
  1285.          * k₂⁰ = &Sigma;<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
  1286.          * </p>
  1287.          * @param kMax max value fot k index
  1288.          * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
  1289.          * @return the coefficient k₂⁰
  1290.          */
  1291.         private double computeK20(final int kMax, final double[][] currentRhoSigmaj) {
  1292.             double k20 = 0.;

  1293.             for (int kIndex = 1; kIndex <= kMax; kIndex++) {
  1294.                 // After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:
  1295.                 //k₂⁰ = &Sigma;<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
  1296.                 double currentTerm = currentRhoSigmaj[1][kIndex] * currentRhoSigmaj[1][kIndex] +
  1297.                                      currentRhoSigmaj[0][kIndex] * currentRhoSigmaj[0][kIndex];

  1298.                 //multiply by 2 / k²
  1299.                 currentTerm *= 2. / (kIndex * kIndex);

  1300.                 // add the term to the result
  1301.                 k20 += currentTerm;
  1302.             }

  1303.             return k20;
  1304.         }

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

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

  1310.             // Get the True longitude L
  1311.             final double L = meanOrbit.getLv();

  1312.             // Compute the center (l - λ)
  1313.             final double center =  L - meanOrbit.getLM();
  1314.             // Compute (l - λ)²
  1315.             final double center2 = center * center;

  1316.             // Initialize short periodic variations
  1317.             final double[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  1318.             final double[] d1 = slot.dij[1].value(meanOrbit.getDate());
  1319.             final double[] d2 = slot.dij[2].value(meanOrbit.getDate());
  1320.             for (int i = 0; i < 6; i++) {
  1321.                 shortPeriodicVariation[i] += center * d1[i] + center2 * d2[i];
  1322.             }

  1323.             for (int j = 1; j <= JMAX; j++) {
  1324.                 final double[] c = slot.cij[j].value(meanOrbit.getDate());
  1325.                 final double[] s = slot.sij[j].value(meanOrbit.getDate());
  1326.                 final double cos = FastMath.cos(j * L);
  1327.                 final double sin = FastMath.sin(j * L);
  1328.                 for (int i = 0; i < 6; i++) {
  1329.                     // add corresponding term to the short periodic variation
  1330.                     shortPeriodicVariation[i] += c[i] * cos;
  1331.                     shortPeriodicVariation[i] += s[i] * sin;
  1332.                 }
  1333.             }

  1334.             return shortPeriodicVariation;

  1335.         }

  1336.         /** {@inheritDoc} */
  1337.         public String getCoefficientsKeyPrefix() {
  1338.             return coefficientsKeyPrefix;
  1339.         }

  1340.         /** {@inheritDoc}
  1341.          * <p>
  1342.          * For Gaussian forces, there are JMAX cj coefficients,
  1343.          * JMAX sj coefficients and 3 dj coefficients. As JMAX = 12,
  1344.          * this sums up to 27 coefficients. The j index is the integer
  1345.          * multiplier for the true longitude argument in the cj and sj
  1346.          * coefficients and to the degree in  the polynomial dj coefficients.
  1347.          * </p>
  1348.          */
  1349.         @Override
  1350.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected)
  1351.             throws OrekitException {

  1352.             // select the coefficients slot
  1353.             final Slot slot = slots.get(date);

  1354.             final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * JMAX + 3);
  1355.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  1356.             storeIfSelected(coefficients, selected, slot.dij[1].value(date), "d", 1);
  1357.             storeIfSelected(coefficients, selected, slot.dij[2].value(date), "d", 2);
  1358.             for (int j = 1; j <= JMAX; j++) {
  1359.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  1360.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  1361.             }

  1362.             return coefficients;

  1363.         }

  1364.         /** Put a coefficient in a map if selected.
  1365.          * @param map map to populate
  1366.          * @param selected set of coefficients that should be put in the map
  1367.          * (empty set means all coefficients are selected)
  1368.          * @param value coefficient value
  1369.          * @param id coefficient identifier
  1370.          * @param indices list of coefficient indices
  1371.          */
  1372.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
  1373.                                      final double[] value, final String id, final int... indices) {
  1374.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  1375.             keyBuilder.append(id);
  1376.             for (int index : indices) {
  1377.                 keyBuilder.append('[').append(index).append(']');
  1378.             }
  1379.             final String key = keyBuilder.toString();
  1380.             if (selected.isEmpty() || selected.contains(key)) {
  1381.                 map.put(key, value);
  1382.             }
  1383.         }

  1384.         /** Replace the instance with a data transfer object for serialization.
  1385.          * @return data transfer object that will be serialized
  1386.          * @exception NotSerializableException if an additional state provider is not serializable
  1387.          */
  1388.         private Object writeReplace() throws NotSerializableException {

  1389.             // slots transitions
  1390.             final SortedSet<TimeSpanMap.Transition<Slot>> transitions     = slots.getTransitions();
  1391.             final AbsoluteDate[]                          transitionDates = new AbsoluteDate[transitions.size()];
  1392.             final Slot[]                                  allSlots        = new Slot[transitions.size() + 1];
  1393.             int i = 0;
  1394.             for (final TimeSpanMap.Transition<Slot> transition : transitions) {
  1395.                 if (i == 0) {
  1396.                     // slot before the first transition
  1397.                     allSlots[i] = transition.getBefore();
  1398.                 }
  1399.                 if (i < transitionDates.length) {
  1400.                     transitionDates[i] = transition.getDate();
  1401.                     allSlots[++i]      = transition.getAfter();
  1402.                 }
  1403.             }

  1404.             return new DataTransferObject(jMax, interpolationPoints, coefficientsKeyPrefix,
  1405.                                           transitionDates, allSlots);

  1406.         }


  1407.         /** Internal class used only for serialization. */
  1408.         private static class DataTransferObject implements Serializable {

  1409.             /** Serializable UID. */
  1410.             private static final long serialVersionUID = 20160319L;

  1411.             /** Maximum value for j index. */
  1412.             private final int jMax;

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

  1415.             /** Prefix for coefficients keys. */
  1416.             private final String coefficientsKeyPrefix;

  1417.             /** Transitions dates. */
  1418.             private final AbsoluteDate[] transitionDates;

  1419.             /** All slots. */
  1420.             private final Slot[] allSlots;

  1421.             /** Simple constructor.
  1422.              * @param jMax maximum value for j index
  1423.              * @param interpolationPoints number of points used in the interpolation process
  1424.              * @param coefficientsKeyPrefix prefix for coefficients keys
  1425.              * @param transitionDates transitions dates
  1426.              * @param allSlots all slots
  1427.              */
  1428.             DataTransferObject(final int jMax, final int interpolationPoints,
  1429.                                final String coefficientsKeyPrefix,
  1430.                                final AbsoluteDate[] transitionDates, final Slot[] allSlots) {
  1431.                 this.jMax                  = jMax;
  1432.                 this.interpolationPoints   = interpolationPoints;
  1433.                 this.coefficientsKeyPrefix = coefficientsKeyPrefix;
  1434.                 this.transitionDates       = transitionDates;
  1435.                 this.allSlots              = allSlots;
  1436.             }

  1437.             /** Replace the deserialized data transfer object with a {@link GaussianShortPeriodicCoefficients}.
  1438.              * @return replacement {@link GaussianShortPeriodicCoefficients}
  1439.              */
  1440.             private Object readResolve() {

  1441.                 final TimeSpanMap<Slot> slots = new TimeSpanMap<Slot>(allSlots[0]);
  1442.                 for (int i = 0; i < transitionDates.length; ++i) {
  1443.                     slots.addValidAfter(allSlots[i + 1], transitionDates[i]);
  1444.                 }

  1445.                 return new GaussianShortPeriodicCoefficients(coefficientsKeyPrefix, jMax, interpolationPoints, slots);

  1446.             }

  1447.         }

  1448.     }

  1449.     /** The U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients described by
  1450.      * equations 2.5.3-(21) and 2.5.3-(22) from Danielson.
  1451.      * <p>
  1452.      * The index i takes only the values 1 and 2<br>
  1453.      * For U only the index 0 for j is used.
  1454.      * </p>
  1455.      *
  1456.      * @author Petre Bazavan
  1457.      * @author Lucian Barbulescu
  1458.      */
  1459.     private static class UijVijCoefficients {

  1460.         /** The U₁<sup>j</sup> coefficients.
  1461.          * <p>
  1462.          * The first index identifies the Fourier coefficients used<br>
  1463.          * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
  1464.          * The only exception is when j = 0 when only the coefficient for fourier index = 1 (i == 0) is needed.<br>
  1465.          * Also, for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
  1466.          * to compute the coefficients U₂<sup>j</sup>
  1467.          * </p>
  1468.          */
  1469.         private final double[][] u1ij;

  1470.         /** The V₁<sup>j</sup> coefficients.
  1471.          * <p>
  1472.          * The first index identifies the Fourier coefficients used<br>
  1473.          * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
  1474.          * for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
  1475.          * to compute the coefficients V₂<sup>j</sup>
  1476.          * </p>
  1477.          */
  1478.         private final double[][] v1ij;

  1479.         /** The U₂<sup>j</sup> coefficients.
  1480.          * <p>
  1481.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
  1482.          * </p>
  1483.          */
  1484.         private final double[] u2ij;

  1485.         /** The V₂<sup>j</sup> coefficients.
  1486.          * <p>
  1487.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
  1488.          * </p>
  1489.          */
  1490.         private final double[] v2ij;

  1491.         /** The current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients. */
  1492.         private final double[][] currentRhoSigmaj;

  1493.         /** The C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> Fourier coefficients. */
  1494.         private final FourierCjSjCoefficients fourierCjSj;

  1495.         /** The maximum value for j index. */
  1496.         private final int jMax;

  1497.         /** Constructor.
  1498.          * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
  1499.          * @param fourierCjSj the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  1500.          * @param jMax maximum value for j index
  1501.          */
  1502.         UijVijCoefficients(final double[][] currentRhoSigmaj, final FourierCjSjCoefficients fourierCjSj, final int jMax) {
  1503.             this.currentRhoSigmaj = currentRhoSigmaj;
  1504.             this.fourierCjSj = fourierCjSj;
  1505.             this.jMax = jMax;

  1506.             // initialize the internal arrays.
  1507.             this.u1ij = new double[6][2 * jMax + 1];
  1508.             this.v1ij = new double[6][2 * jMax + 1];
  1509.             this.u2ij = new double[jMax + 1];
  1510.             this.v2ij = new double[jMax + 1];

  1511.             //compute the coefficients
  1512.             computeU1V1Coefficients();
  1513.             computeU2V2Coefficients();
  1514.         }

  1515.         /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients. */
  1516.         private void computeU1V1Coefficients() {
  1517.             // generate the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients
  1518.             // for j >= 1
  1519.             // also the U₁⁰ for Fourier index = 1 (i == 0) coefficient will be computed
  1520.             u1ij[0][0] = 0;
  1521.             for (int j = 1; j <= jMax; j++) {
  1522.                 // compute 1 / j
  1523.                 final double ooj = 1. / j;

  1524.                 for (int i = 0; i < 6; i++) {
  1525.                     //j is aready between 1 and J
  1526.                     u1ij[i][j] = fourierCjSj.getSij(i, j);
  1527.                     v1ij[i][j] = fourierCjSj.getCij(i, j);

  1528.                     // 1 - δ<sub>1j</sub> is 1 for all j > 1
  1529.                     if (j > 1) {
  1530.                         // k starts with 1 because j-J is less than or equal to 0
  1531.                         for (int kIndex = 1; kIndex <= j - 1; kIndex++) {
  1532.                             // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
  1533.                             // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
  1534.                             u1ij[i][j] +=   fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[1][kIndex] +
  1535.                                             fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[0][kIndex];

  1536.                             // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  1537.                             // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  1538.                             v1ij[i][j] +=   fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[0][kIndex] -
  1539.                                             fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[1][kIndex];
  1540.                         }
  1541.                     }

  1542.                     // since j must be between 1 and J-1 and is already between 1 and J
  1543.                     // the following sum is skiped only for j = jMax
  1544.                     if (j != jMax) {
  1545.                         for (int kIndex = 1; kIndex <= jMax - j; kIndex++) {
  1546.                             // -C<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub> +
  1547.                             // S<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub>
  1548.                             u1ij[i][j] +=   -fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[1][kIndex] +
  1549.                                             fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[0][kIndex];

  1550.                             // C<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub> +
  1551.                             // S<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub>
  1552.                             v1ij[i][j] +=   fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[0][kIndex] +
  1553.                                             fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[1][kIndex];
  1554.                         }
  1555.                     }

  1556.                     for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  1557.                         // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  1558.                         // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  1559.                         u1ij[i][j] +=   -fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[1][j + kIndex] -
  1560.                                         fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[0][j + kIndex];

  1561.                         // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  1562.                         // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  1563.                         v1ij[i][j] +=   fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[0][j + kIndex] +
  1564.                                         fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[1][j + kIndex];
  1565.                     }

  1566.                     // divide by 1 / j
  1567.                     u1ij[i][j] *= -ooj;
  1568.                     v1ij[i][j] *= ooj;

  1569.                     // if index = 1 (i == 0) add the computed terms to U₁⁰
  1570.                     if (i == 0) {
  1571.                         //- (U₁<sup>j</sup> * ρ<sub>j</sub> + V₁<sup>j</sup> * σ<sub>j</sub>
  1572.                         u1ij[0][0] += -u1ij[0][j] * currentRhoSigmaj[0][j] - v1ij[0][j] * currentRhoSigmaj[1][j];
  1573.                     }
  1574.                 }
  1575.             }

  1576.             // Terms with j > jMax are required only when computing the coefficients
  1577.             // U₂<sup>j</sup> and V₂<sup>j</sup>
  1578.             // and those coefficients are only required for Fourier index = 1 (i == 0).
  1579.             for (int j = jMax + 1; j <= 2 * jMax; j++) {
  1580.                 // compute 1 / j
  1581.                 final double ooj = 1. / j;
  1582.                 //the value of i is 0
  1583.                 u1ij[0][j] = 0.;
  1584.                 v1ij[0][j] = 0.;

  1585.                 //k starts from j-J as it is always greater than or equal to 1
  1586.                 for (int kIndex = j - jMax; kIndex <= j - 1; kIndex++) {
  1587.                     // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
  1588.                     // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
  1589.                     u1ij[0][j] +=   fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[1][kIndex] +
  1590.                                     fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[0][kIndex];

  1591.                     // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  1592.                     // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  1593.                     v1ij[0][j] +=   fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[0][kIndex] -
  1594.                                     fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[1][kIndex];
  1595.                 }
  1596.                 for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  1597.                     // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  1598.                     // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  1599.                     u1ij[0][j] +=   -fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[1][j + kIndex] -
  1600.                                     fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[0][j + kIndex];

  1601.                     // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  1602.                     // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  1603.                     v1ij[0][j] +=   fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[0][j + kIndex] +
  1604.                                     fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[1][j + kIndex];
  1605.                 }

  1606.                 // divide by 1 / j
  1607.                 u1ij[0][j] *= -ooj;
  1608.                 v1ij[0][j] *= ooj;
  1609.             }
  1610.         }

  1611.         /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
  1612.          * <p>
  1613.          * Only the coefficients for Fourier index = 1 (i == 0) are required.
  1614.          * </p>
  1615.          */
  1616.         private void computeU2V2Coefficients() {
  1617.             for (int j = 1; j <= jMax; j++) {
  1618.                 // compute 1 / j
  1619.                 final double ooj = 1. / j;

  1620.                 // only the values for i == 0 are computed
  1621.                 u2ij[j] = v1ij[0][j];
  1622.                 v2ij[j] = u1ij[0][j];

  1623.                 // 1 - δ<sub>1j</sub> is 1 for all j > 1
  1624.                 if (j > 1) {
  1625.                     for (int l = 1; l <= j - 1; l++) {
  1626.                         // U₁<sup>j-l</sup> * σ<sub>l</sub> +
  1627.                         // V₁<sup>j-l</sup> * ρ<sub>l</sub>
  1628.                         u2ij[j] +=   u1ij[0][j - l] * currentRhoSigmaj[1][l] +
  1629.                                      v1ij[0][j - l] * currentRhoSigmaj[0][l];

  1630.                         // U₁<sup>j-l</sup> * ρ<sub>l</sub> -
  1631.                         // V₁<sup>j-l</sup> * σ<sub>l</sub>
  1632.                         v2ij[j] +=   u1ij[0][j - l] * currentRhoSigmaj[0][l] -
  1633.                                      v1ij[0][j - l] * currentRhoSigmaj[1][l];
  1634.                     }
  1635.                 }

  1636.                 for (int l = 1; l <= jMax; l++) {
  1637.                     // -U₁<sup>j+l</sup> * σ<sub>l</sub> +
  1638.                     // U₁<sup>l</sup> * σ<sub>j+l</sub> +
  1639.                     // V₁<sup>j+l</sup> * ρ<sub>l</sub> -
  1640.                     // V₁<sup>l</sup> * ρ<sub>j+l</sub>
  1641.                     u2ij[j] +=   -u1ij[0][j + l] * currentRhoSigmaj[1][l] +
  1642.                                   u1ij[0][l] * currentRhoSigmaj[1][j + l] +
  1643.                                   v1ij[0][j + l] * currentRhoSigmaj[0][l] -
  1644.                                   v1ij[0][l] * currentRhoSigmaj[0][j + l];

  1645.                     // U₁<sup>j+l</sup> * ρ<sub>l</sub> +
  1646.                     // U₁<sup>l</sup> * ρ<sub>j+l</sub> +
  1647.                     // V₁<sup>j+l</sup> * σ<sub>l</sub> +
  1648.                     // V₁<sup>l</sup> * σ<sub>j+l</sub>
  1649.                     u2ij[j] +=   u1ij[0][j + l] * currentRhoSigmaj[0][l] +
  1650.                                  u1ij[0][l] * currentRhoSigmaj[0][j + l] +
  1651.                                  v1ij[0][j + l] * currentRhoSigmaj[1][l] +
  1652.                                  v1ij[0][l] * currentRhoSigmaj[1][j + l];
  1653.                 }

  1654.                 // divide by 1 / j
  1655.                 u2ij[j] *= -ooj;
  1656.                 v2ij[j] *= ooj;
  1657.             }
  1658.         }

  1659.         /** Get the coefficient U₁<sup>j</sup> for Fourier index i.
  1660.          *
  1661.          * @param j j index
  1662.          * @param i Fourier index (starts at 0)
  1663.          * @return the coefficient U₁<sup>j</sup> for the given Fourier index i
  1664.          */
  1665.         public double getU1(final int j, final int i) {
  1666.             return u1ij[i][j];
  1667.         }

  1668.         /** Get the coefficient V₁<sup>j</sup> for Fourier index i.
  1669.          *
  1670.          * @param j j index
  1671.          * @param i Fourier index (starts at 0)
  1672.          * @return the coefficient V₁<sup>j</sup> for the given Fourier index i
  1673.          */
  1674.         public double getV1(final int j, final int i) {
  1675.             return v1ij[i][j];
  1676.         }

  1677.         /** Get the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0).
  1678.          *
  1679.          * @param j j index
  1680.          * @return the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0)
  1681.          */
  1682.         public double getU2(final int j) {
  1683.             return u2ij[j];
  1684.         }

  1685.         /** Get the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0).
  1686.          *
  1687.          * @param j j index
  1688.          * @return the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0)
  1689.          */
  1690.         public double getV2(final int j) {
  1691.             return v2ij[j];
  1692.         }
  1693.     }

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

  1696.         /** Serializable UID. */
  1697.         private static final long serialVersionUID = 20160319L;

  1698.         /**The coefficients D<sub>i</sub><sup>j</sup>.
  1699.          * <p>
  1700.          * Only for j = 1 and j = 2 the coefficients are not 0. <br>
  1701.          * i corresponds to the equinoctial element, as follows:
  1702.          * - i=0 for a <br/>
  1703.          * - i=1 for k <br/>
  1704.          * - i=2 for h <br/>
  1705.          * - i=3 for q <br/>
  1706.          * - i=4 for p <br/>
  1707.          * - i=5 for λ <br/>
  1708.          * </p>
  1709.          */
  1710.         private final ShortPeriodicsInterpolatedCoefficient[] dij;

  1711.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  1712.          * <p>
  1713.          * The index order is cij[j][i] <br/>
  1714.          * i corresponds to the equinoctial element, as follows: <br/>
  1715.          * - i=0 for a <br/>
  1716.          * - i=1 for k <br/>
  1717.          * - i=2 for h <br/>
  1718.          * - i=3 for q <br/>
  1719.          * - i=4 for p <br/>
  1720.          * - i=5 for λ <br/>
  1721.          * </p>
  1722.          */
  1723.         private final ShortPeriodicsInterpolatedCoefficient[] cij;

  1724.         /** The coefficients S<sub>i</sub><sup>j</sup>.
  1725.          * <p>
  1726.          * The index order is sij[j][i] <br/>
  1727.          * i corresponds to the equinoctial element, as follows: <br/>
  1728.          * - i=0 for a <br/>
  1729.          * - i=1 for k <br/>
  1730.          * - i=2 for h <br/>
  1731.          * - i=3 for q <br/>
  1732.          * - i=4 for p <br/>
  1733.          * - i=5 for λ <br/>
  1734.          * </p>
  1735.          */
  1736.         private final ShortPeriodicsInterpolatedCoefficient[] sij;

  1737.         /** Simple constructor.
  1738.          *  @param jMax maximum value for j index
  1739.          *  @param interpolationPoints number of points used in the interpolation process
  1740.          */
  1741.         Slot(final int jMax, final int interpolationPoints) {

  1742.             dij = new ShortPeriodicsInterpolatedCoefficient[3];
  1743.             cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
  1744.             sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];

  1745.             // Initialize the C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup> and D<sub>i</sub><sup>j</sup> coefficients
  1746.             for (int j = 0; j <= jMax; j++) {
  1747.                 cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  1748.                 if (j > 0) {
  1749.                     sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  1750.                 }
  1751.                 // Initialize only the non-zero D<sub>i</sub><sup>j</sup> coefficients
  1752.                 if (j == 1 || j == 2) {
  1753.                     dij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  1754.                 }
  1755.             }

  1756.         }

  1757.     }

  1758. }