AbstractGaussianContribution.java

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

  19. import java.lang.reflect.Array;
  20. import java.util.ArrayList;
  21. import java.util.Collections;
  22. import java.util.HashMap;
  23. import java.util.List;
  24. import java.util.Map;
  25. import java.util.Set;

  27. import org.hipparchus.Field;
  28. import org.hipparchus.RealFieldElement;
  29. import org.hipparchus.analysis.RealFieldUnivariateVectorFunction;
  30. import org.hipparchus.analysis.UnivariateVectorFunction;
  31. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  32. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  33. import org.hipparchus.util.FastMath;
  34. import org.hipparchus.util.MathArrays;
  35. import org.orekit.attitudes.Attitude;
  36. import org.orekit.attitudes.AttitudeProvider;
  37. import org.orekit.attitudes.FieldAttitude;
  38. import org.orekit.forces.ForceModel;
  39. import org.orekit.orbits.EquinoctialOrbit;
  40. import org.orekit.orbits.FieldEquinoctialOrbit;
  41. import org.orekit.orbits.FieldOrbit;
  42. import org.orekit.orbits.Orbit;
  43. import org.orekit.orbits.OrbitType;
  44. import org.orekit.orbits.PositionAngle;
  45. import org.orekit.propagation.FieldSpacecraftState;
  46. import org.orekit.propagation.PropagationType;
  47. import org.orekit.propagation.SpacecraftState;
  48. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  49. import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
  50. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  51. import org.orekit.propagation.semianalytical.dsst.utilities.FieldCjSjCoefficient;
  52. import org.orekit.propagation.semianalytical.dsst.utilities.FieldShortPeriodicsInterpolatedCoefficient;
  53. import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
  54. import org.orekit.time.AbsoluteDate;
  55. import org.orekit.time.FieldAbsoluteDate;
  56. import org.orekit.utils.FieldTimeSpanMap;
  57. import org.orekit.utils.ParameterDriver;
  58. import org.orekit.utils.TimeSpanMap;

  60. /** Common handling of {@link DSSTForceModel} methods for Gaussian contributions to DSST propagation.
  61.  * <p>
  62.  * This abstract class allows to provide easily a subset of {@link DSSTForceModel} methods
  63.  * for specific Gaussian contributions.
  64.  * </p><p>
  65.  * This class implements the notion of numerical averaging of the DSST theory.
  66.  * Numerical averaging is mainly used for non-conservative disturbing forces such as
  67.  * atmospheric drag and solar radiation pressure.
  68.  * </p><p>
  69.  * Gaussian contributions can be expressed as: da<sub>i</sub>/dt = δa<sub>i</sub>/δv . q<br>
  70.  * where:
  71.  * <ul>
  72.  * <li>a<sub>i</sub> are the six equinoctial elements</li>
  73.  * <li>v is the velocity vector</li>
  74.  * <li>q is the perturbing acceleration due to the considered force</li>
  75.  * </ul>
  76.  *
  77.  * <p> The averaging process and other considerations lead to integrate this contribution
  78.  * over the true longitude L possibly taking into account some limits.
  79.  *
  80.  * <p> To create a numerically averaged contribution, one needs only to provide a
  81.  * {@link ForceModel} and to implement in the derived class the methods:
  82.  * {@link #getLLimits(SpacecraftState, AuxiliaryElements)}
  83.  * and {@link #getParametersDriversWithoutMu()}.
  84.  * </p>
  85.  * @author Pascal Parraud
  86.  * @author Bryan Cazabonne (field translation)
  87.  */
  88. public abstract class AbstractGaussianContribution implements DSSTForceModel {

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

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

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

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

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

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

  125.     /** Contribution to be numerically averaged. */
  126.     private final ForceModel contribution;

  128.     /** Gauss integrator. */
  129.     private final double threshold;

  131.     /** Gauss integrator. */
  132.     private GaussQuadrature integrator;

  134.     /** Flag for Gauss order computation. */
  135.     private boolean isDirty;

  137.     /** Attitude provider. */
  138.     private AttitudeProvider attitudeProvider;

  140.     /** Prefix for coefficients keys. */
  141.     private final String coefficientsKeyPrefix;

  143.     /** Short period terms. */
  144.     private GaussianShortPeriodicCoefficients gaussianSPCoefs;

  146.     /** Short period terms. */
  147.     private Map<Field<?>, FieldGaussianShortPeriodicCoefficients<?>> gaussianFieldSPCoefs;

  149.     /** Driver for gravitational parameter. */
  150.     private final ParameterDriver gmParameterDriver;

  152.     /** Build a new instance.
  153.      *  @param coefficientsKeyPrefix prefix for coefficients keys
  154.      *  @param threshold tolerance for the choice of the Gauss quadrature order
  155.      *  @param contribution the {@link ForceModel} to be numerically averaged
  156.      *  @param mu central attraction coefficient
  157.      */
  158.     protected AbstractGaussianContribution(final String coefficientsKeyPrefix,
  159.                                            final double threshold,
  160.                                            final ForceModel contribution,
  161.                                            final double mu) {

  163.         gmParameterDriver = new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
  164.                                                 mu, MU_SCALE,
  165.                                                 0.0, Double.POSITIVE_INFINITY);

  167.         this.coefficientsKeyPrefix = coefficientsKeyPrefix;
  168.         this.contribution          = contribution;
  169.         this.threshold             = threshold;
  170.         this.integrator            = new GaussQuadrature(GAUSS_ORDER[MAX_ORDER_RANK]);
  171.         this.isDirty               = true;

  173.         gaussianFieldSPCoefs       = new HashMap<>();
  174.     }

  176.     /** {@inheritDoc} */
  177.     @Override
  178.     public ParameterDriver[] getParametersDrivers() {

  180.         final ParameterDriver[] driversWithoutMu = getParametersDriversWithoutMu();

  182.         // + 1 for central attraction coefficient driver
  183.         final ParameterDriver[] drivers = new ParameterDriver[driversWithoutMu.length + 1];

  185.         int index = 0;
  186.         for (final ParameterDriver driver : driversWithoutMu) {
  187.             drivers[index] = driver;
  188.             index++;
  189.         }

  191.         // We put central attraction coefficient driver at the end of the array
  192.         drivers[driversWithoutMu.length] = gmParameterDriver;
  193.         return drivers;

  195.     }

  197.     /** Get the drivers for force model parameters except the one for the central attraction coefficient.
  198.      * <p>
  199.      * The driver for central attraction coefficient is automatically
  200.      * added at the last element of the {@link ParameterDriver} array
  201.      * into {@link #getParametersDrivers()} method.
  202.      * </p>
  203.      * @return drivers for force model parameters
  204.      */
  205.     protected abstract ParameterDriver[] getParametersDriversWithoutMu();

  207.     /** {@inheritDoc} */
  208.     @Override
  209.     public List<ShortPeriodTerms> initialize(final AuxiliaryElements auxiliaryElements,
  210.                                              final PropagationType type,
  211.                                              final double[] parameters) {

  213.         final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
  214.         gaussianSPCoefs = new GaussianShortPeriodicCoefficients(coefficientsKeyPrefix,
  215.                                                                 JMAX, INTERPOLATION_POINTS,
  216.                                                                 new TimeSpanMap<Slot>(new Slot(JMAX, INTERPOLATION_POINTS)));
  217.         list.add(gaussianSPCoefs);
  218.         return list;

  220.     }

  222.     /** {@inheritDoc} */
  223.     @Override
  224.     public <T extends RealFieldElement<T>> List<FieldShortPeriodTerms<T>> initialize(final FieldAuxiliaryElements<T> auxiliaryElements,
  225.                                                                                      final PropagationType type,
  226.                                                                                      final T[] parameters) {

  228.         final Field<T> field = auxiliaryElements.getDate().getField();

  230.         final FieldGaussianShortPeriodicCoefficients<T> fgspc =
  231.                         new FieldGaussianShortPeriodicCoefficients<>(coefficientsKeyPrefix,
  232.                                                                      JMAX, INTERPOLATION_POINTS,
  233.                                                                      new FieldTimeSpanMap<>(new FieldSlot<>(JMAX,
  234.                                                                                                             INTERPOLATION_POINTS),
  235.                                                                                             field));
  236.         gaussianFieldSPCoefs.put(field, fgspc);
  237.         return Collections.singletonList(fgspc);
  238.     }

  240.     /** Performs initialization at each integration step for the current force model.
  241.      *  <p>
  242.      *  This method aims at being called before mean elements rates computation.
  243.      *  </p>
  244.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  245.      *  @param parameters parameters values of the force model parameters
  246.      *  @return new force model context
  247.      */
  248.     private AbstractGaussianContributionContext initializeStep(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
  249.         return new AbstractGaussianContributionContext(auxiliaryElements, parameters);
  250.     }

  252.     /** Performs initialization at each integration step for the current force model.
  253.      *  <p>
  254.      *  This method aims at being called before mean elements rates computation.
  255.      *  </p>
  256.      *  @param <T> type of the elements
  257.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  258.      *  @param parameters parameters values of the force model parameters
  259.      *  @return new force model context
  260.      */
  261.     private <T extends RealFieldElement<T>> FieldAbstractGaussianContributionContext<T> initializeStep(final FieldAuxiliaryElements<T> auxiliaryElements,
  262.                                                                                                        final T[] parameters) {
  263.         return new FieldAbstractGaussianContributionContext<>(auxiliaryElements, parameters);
  264.     }

  266.     /** {@inheritDoc} */
  267.     @Override
  268.     public double[] getMeanElementRate(final SpacecraftState state,
  269.                                        final AuxiliaryElements auxiliaryElements, final double[] parameters) {

  271.         // Container for attributes
  272.         final AbstractGaussianContributionContext context = initializeStep(auxiliaryElements, parameters);
  273.         double[] meanElementRate = new double[6];
  274.         // Computes the limits for the integral
  275.         final double[] ll = getLLimits(state, auxiliaryElements);
  276.         // Computes integrated mean element rates if Llow < Lhigh
  277.         if (ll[0] < ll[1]) {
  278.             meanElementRate = getMeanElementRate(state, integrator, ll[0], ll[1], context, parameters);
  279.             if (isDirty) {
  280.                 boolean next = true;
  281.                 for (int i = 0; i < MAX_ORDER_RANK && next; i++) {
  282.                     final double[] meanRates = getMeanElementRate(state, new GaussQuadrature(GAUSS_ORDER[i]), ll[0], ll[1], context, parameters);
  283.                     if (getRatesDiff(meanElementRate, meanRates, context) < threshold) {
  284.                         integrator = new GaussQuadrature(GAUSS_ORDER[i]);
  285.                         next = false;
  286.                     }
  287.                 }
  288.                 isDirty = false;
  289.             }
  290.         }
  291.         return meanElementRate;
  292.     }

  294.     /** {@inheritDoc} */
  295.     @Override
  296.     public <T extends RealFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> state,
  297.                                                                   final FieldAuxiliaryElements<T> auxiliaryElements,
  298.                                                                   final T[] parameters) {

  300.         // Container for attributes
  301.         final FieldAbstractGaussianContributionContext<T> context = initializeStep(auxiliaryElements, parameters);
  302.         final Field<T> field = state.getDate().getField();

  304.         T[] meanElementRate = MathArrays.buildArray(field, 6);
  305.         // Computes the limits for the integral
  306.         final T[] ll = getLLimits(state, auxiliaryElements);
  307.         // Computes integrated mean element rates if Llow < Lhigh
  308.         if (ll[0].getReal() < ll[1].getReal()) {
  309.             meanElementRate = getMeanElementRate(state, integrator, ll[0], ll[1], context, parameters);
  310.             if (isDirty) {
  311.                 boolean next = true;
  312.                 for (int i = 0; i < MAX_ORDER_RANK && next; i++) {
  313.                     final T[] meanRates = getMeanElementRate(state, new GaussQuadrature(GAUSS_ORDER[i]), ll[0], ll[1], context, parameters);
  314.                     if (getRatesDiff(meanElementRate, meanRates, context).getReal() < threshold) {
  315.                         integrator = new GaussQuadrature(GAUSS_ORDER[i]);
  316.                         next = false;
  317.                     }
  318.                 }
  319.                 isDirty = false;
  320.             }
  321.         }

  323.         return meanElementRate;
  324.     }

  326.     /** Compute the limits in L, the true longitude, for integration.
  327.      *
  328.      *  @param  state current state information: date, kinematics, attitude
  329.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  330.      *  @return the integration limits in L
  331.      */
  332.     protected abstract double[] getLLimits(SpacecraftState state, AuxiliaryElements auxiliaryElements);

  334.     /** Compute the limits in L, the true longitude, for integration.
  335.      *
  336.      *  @param <T> type of the elements
  337.      *  @param state current state information: date, kinematics, attitude
  338.      *  @param auxiliaryElements auxiliary elements related to the current orbit
  339.      *  @return the integration limits in L
  340.      */
  341.     protected abstract <T extends RealFieldElement<T>> T[] getLLimits(FieldSpacecraftState<T> state,
  342.                                                                       FieldAuxiliaryElements<T> auxiliaryElements);

  344.     /** Computes the mean equinoctial elements rates da<sub>i</sub> / dt.
  345.     *
  346.     *  @param state current state
  347.     *  @param gauss Gauss quadrature
  348.     *  @param low lower bound of the integral interval
  349.     *  @param high upper bound of the integral interval
  350.     *  @param context container for attributes
  351.     *  @param parameters values of the force model parameters
  352.     *  @return the mean element rates
  353.     */
  354.     private double[] getMeanElementRate(final SpacecraftState state,
  355.                                         final GaussQuadrature gauss,
  356.                                         final double low,
  357.                                         final double high,
  358.                                         final AbstractGaussianContributionContext context,
  359.                                         final double[] parameters) {

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

  364.         final double[] meanElementRate = gauss.integrate(new IntegrableFunction(state, true, 0, parameters), low, high);

  366.         // Constant multiplier for integral
  367.         final double coef = 1. / (2. * FastMath.PI * auxiliaryElements.getB());
  368.         // Corrects mean element rates
  369.         for (int i = 0; i < 6; i++) {
  370.             meanElementRate[i] *= coef;
  371.         }
  372.         return meanElementRate;
  373.     }

  375.     /** Computes the mean equinoctial elements rates da<sub>i</sub> / dt.
  376.     *
  377.     *  @param <T> type of the elements
  378.     *  @param state current state
  379.     *  @param gauss Gauss quadrature
  380.     *  @param low lower bound of the integral interval
  381.     *  @param high upper bound of the integral interval
  382.     *  @param context container for attributes
  383.     *  @param parameters values of the force model parameters
  384.     *  @return the mean element rates
  385.     */
  386.     private <T extends RealFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> state,
  387.                                                                    final GaussQuadrature gauss,
  388.                                                                    final T low,
  389.                                                                    final T high,
  390.                                                                    final FieldAbstractGaussianContributionContext<T> context,
  391.                                                                    final T[] parameters) {

  393.         // Field
  394.         final Field<T> field = context.getA().getField();

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

  399.         final T[] meanElementRate = gauss.integrate(new FieldIntegrableFunction<>(state, true, 0, parameters, field), low, high, field);
  400.         // Constant multiplier for integral
  401.         final T coef = auxiliaryElements.getB().multiply(FastMath.PI).multiply(2.).reciprocal();
  402.         // Corrects mean element rates
  403.         for (int i = 0; i < 6; i++) {
  404.             meanElementRate[i] = meanElementRate[i].multiply(coef);
  405.         }
  406.         return meanElementRate;
  407.     }

  409.     /** Estimates the weighted magnitude of the difference between 2 sets of equinoctial elements rates.
  410.      *
  411.      *  @param meanRef reference rates
  412.      *  @param meanCur current rates
  413.      *  @param context container for attributes
  414.      *  @return estimated magnitude of weighted differences
  415.      */
  416.     private double getRatesDiff(final double[] meanRef, final double[] meanCur, final AbstractGaussianContributionContext context) {

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

  421.         double maxDiff = FastMath.abs(meanRef[0] - meanCur[0]) / auxiliaryElements.getSma();
  422.         // Corrects mean element rates
  423.         for (int i = 1; i < meanRef.length; i++) {
  424.             maxDiff = FastMath.max(maxDiff, FastMath.abs(meanRef[i] - meanCur[i]));
  425.         }
  426.         return maxDiff;
  427.     }

  429.     /** Estimates the weighted magnitude of the difference between 2 sets of equinoctial elements rates.
  430.     *
  431.     *  @param <T> type of the elements
  432.     *  @param meanRef reference rates
  433.     *  @param meanCur current rates
  434.     *  @param context container for attributes
  435.     *  @return estimated magnitude of weighted differences
  436.     */
  437.     private <T extends RealFieldElement<T>> T getRatesDiff(final T[] meanRef, final T[] meanCur, final FieldAbstractGaussianContributionContext<T> context) {

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

  442.         T maxDiff = FastMath.abs(meanRef[0].subtract(meanCur[0])).divide(auxiliaryElements.getSma());;
  443.         // Corrects mean element rates
  444.         for (int i = 1; i < meanRef.length; i++) {
  445.             maxDiff = FastMath.max(maxDiff, FastMath.abs(meanRef[i].subtract(meanCur[i])));
  446.         }
  447.         return maxDiff;
  448.     }

  450.     /** {@inheritDoc} */
  451.     @Override
  452.     public void registerAttitudeProvider(final AttitudeProvider provider) {
  453.         this.attitudeProvider = provider;
  454.     }

  456.     /** {@inheritDoc} */
  457.     @Override
  458.     public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {

  460.         final Slot slot = gaussianSPCoefs.createSlot(meanStates);
  461.         for (final SpacecraftState meanState : meanStates) {

  463.             // Auxiliary elements related to the current orbit
  464.             final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(meanState.getOrbit(), I);

  466.             // Container of attributes
  467.             final AbstractGaussianContributionContext context = initializeStep(auxiliaryElements, parameters);

  469.             // Compute rhoj and sigmaj
  470.             final double[][] currentRhoSigmaj = computeRhoSigmaCoefficients(meanState.getDate(), auxiliaryElements);

  472.             // Generate the Cij and Sij coefficients
  473.             final FourierCjSjCoefficients fourierCjSj = new FourierCjSjCoefficients(meanState, JMAX, auxiliaryElements, parameters);

  475.             // Generate the Uij and Vij coefficients
  476.             final UijVijCoefficients uijvij = new UijVijCoefficients(currentRhoSigmaj, fourierCjSj, JMAX);

  478.             gaussianSPCoefs.computeCoefficients(meanState, slot, fourierCjSj, uijvij, context.getMeanMotion(), auxiliaryElements.getSma());

  480.         }

  482.     }

  484.     /** {@inheritDoc} */
  485.     @Override
  486.     @SuppressWarnings("unchecked")
  487.     public <T extends RealFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
  488.                                                                        final FieldSpacecraftState<T>... meanStates) {

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

  493.         final FieldGaussianShortPeriodicCoefficients<T> fgspc = (FieldGaussianShortPeriodicCoefficients<T>) gaussianFieldSPCoefs.get(field);
  494.         final FieldSlot<T> slot = fgspc.createSlot(meanStates);
  495.         for (final FieldSpacecraftState<T> meanState : meanStates) {

  497.             // Auxiliary elements related to the current orbit
  498.             final FieldAuxiliaryElements<T> auxiliaryElements = new FieldAuxiliaryElements<>(meanState.getOrbit(), I);

  500.             // Container of attributes
  501.             final FieldAbstractGaussianContributionContext<T> context = initializeStep(auxiliaryElements, parameters);

  503.             // Compute rhoj and sigmaj
  504.             final T[][] currentRhoSigmaj = computeRhoSigmaCoefficients(meanState.getDate(), context, field);

  506.             // Generate the Cij and Sij coefficients
  507.             final FieldFourierCjSjCoefficients<T> fourierCjSj = new FieldFourierCjSjCoefficients<>(meanState, JMAX, auxiliaryElements, parameters, field);

  509.             // Generate the Uij and Vij coefficients
  510.             final FieldUijVijCoefficients<T> uijvij = new FieldUijVijCoefficients<>(currentRhoSigmaj, fourierCjSj, JMAX, field);

  512.             fgspc.computeCoefficients(meanState, slot, fourierCjSj, uijvij, context.getMeanMotion(), auxiliaryElements.getSma(), field);

  514.         }

  516.     }

  518.     /**
  519.      * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
  520.      * <p>
  521.      * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
  522.      *  ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
  523.      *  σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
  524.      * </p>
  525.      * @param date current date
  526.      * @param auxiliaryElements auxiliary elements related to the current orbit
  527.      * @return computed coefficients
  528.      */
  529.     private double[][] computeRhoSigmaCoefficients(final AbsoluteDate date, final AuxiliaryElements auxiliaryElements) {
  530.         final double[][] currentRhoSigmaj = new double[2][3 * JMAX + 1];
  531.         final CjSjCoefficient cjsjKH = new CjSjCoefficient(auxiliaryElements.getK(), auxiliaryElements.getH());
  532.         final double b = 1. / (1 + auxiliaryElements.getB());

  534.         // (-b)<sup>j</sup>
  535.         double mbtj = 1;

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

  539.             //Compute current rho and sigma;
  540.             mbtj *= -b;
  541.             final double coef = (1 + j * auxiliaryElements.getB()) * mbtj;
  542.             currentRhoSigmaj[0][j] = coef * cjsjKH.getCj(j);
  543.             currentRhoSigmaj[1][j] = coef * cjsjKH.getSj(j);
  544.         }
  545.         return currentRhoSigmaj;
  546.     }

  548.     /**
  549.      * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
  550.      * <p>
  551.      * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
  552.      *  ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
  553.      *  σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
  554.      * </p>
  555.      * @param <T> type of the elements
  556.      * @param date current date
  557.      * @param context container for attributes
  558.      * @param field field used by default
  559.      * @return computed coefficients
  560.      */
  561.     private <T extends RealFieldElement<T>> T[][] computeRhoSigmaCoefficients(final FieldAbsoluteDate<T> date,
  562.                                                                               final FieldAbstractGaussianContributionContext<T> context,
  563.                                                                               final Field<T> field) {
  564.         // zero
  565.         final T zero = field.getZero();

  567.         final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
  568.         final T[][] currentRhoSigmaj = MathArrays.buildArray(field, 2, 3 * JMAX + 1);
  569.         final FieldCjSjCoefficient<T> cjsjKH = new FieldCjSjCoefficient<>(auxiliaryElements.getK(), auxiliaryElements.getH(), field);
  570.         final T b = auxiliaryElements.getB().add(1.).reciprocal();

  572.         // (-b)<sup>j</sup>
  573.         T mbtj = zero.add(1.);

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

  577.             //Compute current rho and sigma;
  578.             mbtj = mbtj.multiply(b.negate());
  579.             final T coef = mbtj.multiply(auxiliaryElements.getB().multiply(j).add(1.));
  580.             currentRhoSigmaj[0][j] = coef.multiply(cjsjKH.getCj(j));
  581.             currentRhoSigmaj[1][j] = coef.multiply(cjsjKH.getSj(j));
  582.         }
  583.         return currentRhoSigmaj;
  584.     }

  586.     /** Internal class for numerical quadrature.
  587.      *  <p>
  588.      *  This class is a rewrite of {@link IntegrableFunction} for field elements
  589.      *  </p>
  590.      */
  591.     private class FieldIntegrableFunction <T extends RealFieldElement<T>> implements RealFieldUnivariateVectorFunction<T> {

  593.         /** Current state. */
  594.         private final FieldSpacecraftState<T> state;

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

  600.         /** The j index.
  601.          * <p>
  602.          * Used only for short periodic variation. Ignored for mean elements variation.
  603.          * </p> */
  604.         private final int j;

  606.         /** Container for attributes. */
  607.         private final FieldAbstractGaussianContributionContext<T> context;

  609.         /** Auxiliary Elements. */
  610.         private final FieldAuxiliaryElements<T> auxiliaryElements;

  612.         /** Drivers for solar radiation and atmospheric drag forces. */
  613.         private final T[] parameters;

  615.         /** Build a new instance with a new field.
  616.          *  @param  state current state information: date, kinematics, attitude
  617.          *  @param meanMode if true return the value associated to the mean elements variation,
  618.          *                  if false return the values associated to the short periodic elements variation
  619.          *  @param j the j index. used only for short periodic variation. Ignored for mean elements variation.
  620.          *  @param parameters values of the force model parameters
  621.          *  @param field field utilized by default
  622.          */
  623.         FieldIntegrableFunction(final FieldSpacecraftState<T> state,
  624.                                 final boolean meanMode,
  625.                                 final int j,
  626.                                 final T[] parameters,
  627.                                 final Field<T> field) {

  629.             this.meanMode = meanMode;
  630.             this.j = j;
  631.             this.parameters = parameters;
  632.             this.auxiliaryElements = new FieldAuxiliaryElements<>(state.getOrbit(), I);
  633.             this.context = new FieldAbstractGaussianContributionContext<>(auxiliaryElements, parameters);
  634.             // remove derivatives from state
  635.             final T[] stateVector = MathArrays.buildArray(field, 6);
  636.             OrbitType.EQUINOCTIAL.mapOrbitToArray(state.getOrbit(), PositionAngle.TRUE, stateVector, null);
  637.             final FieldOrbit<T> fixedOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(stateVector, null, PositionAngle.TRUE,
  638.                                                                            state.getDate(),
  639.                                                                            context.getMu(),
  640.                                                                            state.getFrame());
  641.             this.state = new FieldSpacecraftState<>(fixedOrbit, state.getAttitude(), state.getMass());
  642.         }

  644.         /** {@inheritDoc} */
  645.         @Override
  646.         public T[] value(final T x) {

  648.             // Parameters for array building
  649.             final Field<T> field = auxiliaryElements.getDate().getField();
  650.             final int dimension = 6;

  652.             //Compute the time difference from the true longitude difference
  653.             final T shiftedLm = trueToMean(x);
  654.             final T dLm = shiftedLm.subtract(auxiliaryElements.getLM());
  655.             final T dt = dLm.divide(context.getMeanMotion());

  657.             final T cosL = FastMath.cos(x);
  658.             final T sinL = FastMath.sin(x);
  659.             final T roa  = auxiliaryElements.getB().multiply(auxiliaryElements.getB()).divide(auxiliaryElements.getH().multiply(sinL).add(auxiliaryElements.getK().multiply(cosL)).add(1.));
  660.             final T roa2 = roa.multiply(roa);
  661.             final T r    = auxiliaryElements.getSma().multiply(roa);
  662.             final T X    = r.multiply(cosL);
  663.             final T Y    = r.multiply(sinL);
  664.             final T naob = context.getMeanMotion().multiply(auxiliaryElements.getSma()).divide(auxiliaryElements.getB());
  665.             final T Xdot = naob.multiply(auxiliaryElements.getH().add(sinL)).negate();
  666.             final T Ydot = naob.multiply(auxiliaryElements.getK().add(cosL));
  667.             final FieldVector3D<T> vel = new FieldVector3D<>(Xdot, auxiliaryElements.getVectorF(), Ydot, auxiliaryElements.getVectorG());

  669.             // Compute acceleration
  670.             FieldVector3D<T> acc = FieldVector3D.getZero(field);

  672.             // shift the orbit to dt
  673.             final FieldOrbit<T> shiftedOrbit = state.getOrbit().shiftedBy(dt);

  675.             // Recompose an orbit with time held fixed to be compliant with DSST theory
  676.             final FieldOrbit<T> recomposedOrbit =
  677.                             new FieldEquinoctialOrbit<>(shiftedOrbit.getA(),
  678.                                             shiftedOrbit.getEquinoctialEx(),
  679.                                             shiftedOrbit.getEquinoctialEy(),
  680.                                             shiftedOrbit.getHx(),
  681.                                             shiftedOrbit.getHy(),
  682.                                             shiftedOrbit.getLv(),
  683.                                             PositionAngle.TRUE,
  684.                                             shiftedOrbit.getFrame(),
  685.                                             state.getDate(),
  686.                                             context.getMu());

  688.             // Get the corresponding attitude
  689.             final FieldAttitude<T> recomposedAttitude =
  690.                             attitudeProvider.getAttitude(recomposedOrbit,
  691.                                                          recomposedOrbit.getDate(),
  692.                                                          recomposedOrbit.getFrame());

  694.             // create shifted SpacecraftState with attitude at specified time
  695.             final FieldSpacecraftState<T> shiftedState =
  696.                             new FieldSpacecraftState<>(recomposedOrbit, recomposedAttitude, state.getMass());

  698.             acc = contribution.acceleration(shiftedState, parameters);

  700.             //Compute the derivatives of the elements by the speed
  701.             final T[] deriv = MathArrays.buildArray(field, dimension);
  702.             // da/dv
  703.             deriv[0] = getAoV(vel).dotProduct(acc);
  704.             // dex/dv
  705.             deriv[1] = getKoV(X, Y, Xdot, Ydot).dotProduct(acc);
  706.             // dey/dv
  707.             deriv[2] = getHoV(X, Y, Xdot, Ydot).dotProduct(acc);
  708.             // dhx/dv
  709.             deriv[3] = getQoV(X).dotProduct(acc);
  710.             // dhy/dv
  711.             deriv[4] = getPoV(Y).dotProduct(acc);
  712.             // dλ/dv
  713.             deriv[5] = getLoV(X, Y, Xdot, Ydot).dotProduct(acc);

  715.             // Compute mean elements rates
  716.             T[] val = null;
  717.             if (meanMode) {
  718.                 val = MathArrays.buildArray(field, dimension);
  719.                 for (int i = 0; i < 6; i++) {
  720.                     // da<sub>i</sub>/dt
  721.                     val[i] = deriv[i].multiply(roa2);
  722.                 }
  723.             } else {
  724.                 val = MathArrays.buildArray(field, dimension * 2);
  725.                 //Compute cos(j*L) and sin(j*L);
  726.                 final T cosjL = j == 1 ? cosL : FastMath.cos(x.multiply(j));
  727.                 final T sinjL = j == 1 ? sinL : FastMath.sin(x.multiply(j));

  729.                 for (int i = 0; i < 6; i++) {
  730.                     // da<sub>i</sub>/dv * cos(jL)
  731.                     val[i] = deriv[i].multiply(cosjL);
  732.                     // da<sub>i</sub>/dv * sin(jL)
  733.                     val[i + 6] = deriv[i].multiply(sinjL);
  734.                 }
  735.             }

  737.             return val;
  738.         }

  740.         /** Converts true longitude to mean longitude.
  741.          * @param x True longitude
  742.          * @return Eccentric longitude
  743.          */
  744.         private T trueToMean (final T x) {
  745.             return eccentricToMean(trueToEccentric(x));
  746.         }

  748.         /** Converts true longitude to eccentric longitude.
  749.          * @param lv True longitude
  750.          * @return Eccentric longitude
  751.          */
  752.         private T trueToEccentric (final T lv) {
  753.             final T cosLv   = FastMath.cos(lv);
  754.             final T sinLv   = FastMath.sin(lv);
  755.             final T num     = auxiliaryElements.getH().multiply(cosLv).subtract(auxiliaryElements.getK().multiply(sinLv));
  756.             final T den     = auxiliaryElements.getB().add(auxiliaryElements.getK().multiply(cosLv)).add(auxiliaryElements.getH().multiply(sinLv)).add(1.);
  757.             return FastMath.atan(num.divide(den)).multiply(2.).add(lv);
  758.         }

  760.         /** Converts eccentric longitude to mean longitude.
  761.          * @param le Eccentric longitude
  762.          * @return Mean longitude
  763.          */
  764.         private T eccentricToMean (final T le) {
  765.             return le.subtract(auxiliaryElements.getK().multiply(FastMath.sin(le))).add(auxiliaryElements.getH().multiply(FastMath.cos(le)));
  766.         }

  768.         /** Compute δa/δv.
  769.          *  @param vel satellite velocity
  770.          *  @return δa/δv
  771.          */
  772.         private FieldVector3D<T> getAoV(final FieldVector3D<T> vel) {
  773.             return new FieldVector3D<>(context.getTon2a(), vel);
  774.         }

  776.         /** Compute δh/δv.
  777.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  778.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  779.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  780.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  781.          *  @return δh/δv
  782.          */
  783.         private FieldVector3D<T> getHoV(final T X, final T Y, final T Xdot, final T Ydot) {
  784.             final T kf = (Xdot.multiply(Y).multiply(2.).subtract(X.multiply(Ydot))).multiply(context.getOoMU());
  785.             final T kg = X.multiply(Xdot).multiply(context.getOoMU());
  786.             final T kw = auxiliaryElements.getK().multiply(auxiliaryElements.getQ().multiply(Y).multiply(I).subtract(auxiliaryElements.getP().multiply(X))).multiply(context.getOOAB());
  787.             return new FieldVector3D<>(kf, auxiliaryElements.getVectorF(), kg.negate(), auxiliaryElements.getVectorG(), kw, auxiliaryElements.getVectorW());
  788.         }

  790.         /** Compute δk/δv.
  791.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  792.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  793.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  794.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  795.          *  @return δk/δv
  796.          */
  797.         private FieldVector3D<T> getKoV(final T X, final T Y, final T Xdot, final T Ydot) {
  798.             final T kf = Y.multiply(Ydot).multiply(context.getOoMU());
  799.             final T kg = (X.multiply(Ydot).multiply(2.).subtract(Xdot.multiply(Y))).multiply(context.getOoMU());
  800.             final T kw = auxiliaryElements.getH().multiply(auxiliaryElements.getQ().multiply(Y).multiply(I).subtract(auxiliaryElements.getP().multiply(X))).multiply(context.getOOAB());
  801.             return new FieldVector3D<>(kf.negate(), auxiliaryElements.getVectorF(), kg, auxiliaryElements.getVectorG(), kw.negate(), auxiliaryElements.getVectorW());
  802.         }

  804.         /** Compute δp/δv.
  805.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  806.          *  @return δp/δv
  807.          */
  808.         private FieldVector3D<T> getPoV(final T Y) {
  809.             return new FieldVector3D<>(context.getCo2AB().multiply(Y), auxiliaryElements.getVectorW());
  810.         }

  812.         /** Compute δq/δv.
  813.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  814.          *  @return δq/δv
  815.          */
  816.         private FieldVector3D<T> getQoV(final T X) {
  817.             return new FieldVector3D<>(context.getCo2AB().multiply(X).multiply(I), auxiliaryElements.getVectorW());
  818.         }

  820.         /** Compute δλ/δv.
  821.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  822.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  823.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  824.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  825.          *  @return δλ/δv
  826.          */
  827.         private FieldVector3D<T> getLoV(final T X, final T Y, final T Xdot, final T Ydot) {
  828.             final FieldVector3D<T> pos = new FieldVector3D<>(X, auxiliaryElements.getVectorF(), Y, auxiliaryElements.getVectorG());
  829.             final FieldVector3D<T> v2  = new FieldVector3D<>(auxiliaryElements.getK(), getHoV(X, Y, Xdot, Ydot), auxiliaryElements.getH().negate(), getKoV(X, Y, Xdot, Ydot));
  830.             return new FieldVector3D<>(context.getOOA().multiply(-2.), pos, context.getOoBpo(), v2, context.getOOA().multiply(auxiliaryElements.getQ().multiply(Y).multiply(I).subtract(auxiliaryElements.getP().multiply(X))), auxiliaryElements.getVectorW());
  831.         }

  833.     }

  835.     /** Internal class for numerical quadrature. */
  836.     private class IntegrableFunction implements UnivariateVectorFunction {

  838.         /** Current state. */
  839.         private final SpacecraftState state;

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

  845.         /** The j index.
  846.          * <p>
  847.          * Used only for short periodic variation. Ignored for mean elements variation.
  848.          * </p> */
  849.         private final int j;

  851.         /** Container for attributes. */
  852.         private final AbstractGaussianContributionContext context;

  854.         /** Auxiliary Elements. */
  855.         private final AuxiliaryElements auxiliaryElements;

  857.         /** Drivers for solar radiation and atmospheric drag forces. */
  858.         private final double[] parameters;

  860.         /** Build a new instance.
  861.          *  @param  state current state information: date, kinematics, attitude
  862.          *  @param meanMode if true return the value associated to the mean elements variation,
  863.          *                  if false return the values associated to the short periodic elements variation
  864.          *  @param j the j index. used only for short periodic variation. Ignored for mean elements variation.
  865.          *  @param parameters values of the force model parameters
  866.          */
  867.         IntegrableFunction(final SpacecraftState state,
  868.                            final boolean meanMode,
  869.                            final int j,
  870.                            final double[] parameters) {

  872.             this.meanMode = meanMode;
  873.             this.j = j;
  874.             this.parameters = parameters;
  875.             this.auxiliaryElements = new AuxiliaryElements(state.getOrbit(), I);
  876.             this.context = new AbstractGaussianContributionContext(auxiliaryElements, parameters);
  877.             // remove derivatives from state
  878.             final double[] stateVector = new double[6];
  879.             OrbitType.EQUINOCTIAL.mapOrbitToArray(state.getOrbit(), PositionAngle.TRUE, stateVector, null);
  880.             final Orbit fixedOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(stateVector, null, PositionAngle.TRUE,
  881.                                                                            state.getDate(),
  882.                                                                            context.getMu(),
  883.                                                                            state.getFrame());
  884.             this.state = new SpacecraftState(fixedOrbit, state.getAttitude(), state.getMass());
  885.         }

  887.         /** {@inheritDoc} */
  888.         @Override
  889.         public double[] value(final double x) {

  891.             //Compute the time difference from the true longitude difference
  892.             final double shiftedLm = trueToMean(x);
  893.             final double dLm = shiftedLm - auxiliaryElements.getLM();
  894.             final double dt = dLm / context.getMeanMotion();

  896.             final double cosL = FastMath.cos(x);
  897.             final double sinL = FastMath.sin(x);
  898.             final double roa  = auxiliaryElements.getB() * auxiliaryElements.getB() / (1. + auxiliaryElements.getH() * sinL + auxiliaryElements.getK() * cosL);
  899.             final double roa2 = roa * roa;
  900.             final double r    = auxiliaryElements.getSma() * roa;
  901.             final double X    = r * cosL;
  902.             final double Y    = r * sinL;
  903.             final double naob = context.getMeanMotion() * auxiliaryElements.getSma() / auxiliaryElements.getB();
  904.             final double Xdot = -naob * (auxiliaryElements.getH() + sinL);
  905.             final double Ydot =  naob * (auxiliaryElements.getK() + cosL);
  906.             final Vector3D vel = new Vector3D(Xdot, auxiliaryElements.getVectorF(), Ydot, auxiliaryElements.getVectorG());

  908.             // Compute acceleration
  909.             Vector3D acc = Vector3D.ZERO;

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

  914.             // Recompose an orbit with time held fixed to be compliant with DSST theory
  915.             final Orbit recomposedOrbit =
  916.                             new EquinoctialOrbit(shiftedOrbit.getA(),
  917.                                                  shiftedOrbit.getEquinoctialEx(),
  918.                                                  shiftedOrbit.getEquinoctialEy(),
  919.                                                  shiftedOrbit.getHx(),
  920.                                                  shiftedOrbit.getHy(),
  921.                                                  shiftedOrbit.getLv(),
  922.                                                  PositionAngle.TRUE,
  923.                                                  shiftedOrbit.getFrame(),
  924.                                                  state.getDate(),
  925.                                                  context.getMu());

  927.             // Get the corresponding attitude
  928.             final Attitude recomposedAttitude =
  929.                     attitudeProvider.getAttitude(recomposedOrbit,
  930.                                                  recomposedOrbit.getDate(),
  931.                                                  recomposedOrbit.getFrame());

  933.             // create shifted SpacecraftState with attitude at specified time
  934.             final SpacecraftState shiftedState =
  935.                     new SpacecraftState(recomposedOrbit, recomposedAttitude, state.getMass());

  937.             acc = contribution.acceleration(shiftedState, parameters);

  939.             //Compute the derivatives of the elements by the speed
  940.             final double[] deriv = new double[6];
  941.             // da/dv
  942.             deriv[0] = getAoV(vel).dotProduct(acc);
  943.             // dex/dv
  944.             deriv[1] = getKoV(X, Y, Xdot, Ydot).dotProduct(acc);
  945.             // dey/dv
  946.             deriv[2] = getHoV(X, Y, Xdot, Ydot).dotProduct(acc);
  947.             // dhx/dv
  948.             deriv[3] = getQoV(X).dotProduct(acc);
  949.             // dhy/dv
  950.             deriv[4] = getPoV(Y).dotProduct(acc);
  951.             // dλ/dv
  952.             deriv[5] = getLoV(X, Y, Xdot, Ydot).dotProduct(acc);

  954.             // Compute mean elements rates
  955.             double[] val = null;
  956.             if (meanMode) {
  957.                 val = new double[6];
  958.                 for (int i = 0; i < 6; i++) {
  959.                     // da<sub>i</sub>/dt
  960.                     val[i] = roa2 * deriv[i];
  961.                 }
  962.             } else {
  963.                 val = new double[12];
  964.                 //Compute cos(j*L) and sin(j*L);
  965.                 final double cosjL = j == 1 ? cosL : FastMath.cos(j * x);
  966.                 final double sinjL = j == 1 ? sinL : FastMath.sin(j * x);

  968.                 for (int i = 0; i < 6; i++) {
  969.                     // da<sub>i</sub>/dv * cos(jL)
  970.                     val[i] = cosjL * deriv[i];
  971.                     // da<sub>i</sub>/dv * sin(jL)
  972.                     val[i + 6] = sinjL * deriv[i];
  973.                 }
  974.             }
  975.             return val;
  976.         }

  978.         /** Converts true longitude to eccentric longitude.
  979.          * @param lv True longitude
  980.          * @return Eccentric longitude
  981.          */
  982.         private double trueToEccentric (final double lv) {
  983.             final double cosLv   = FastMath.cos(lv);
  984.             final double sinLv   = FastMath.sin(lv);
  985.             final double num     = auxiliaryElements.getH() * cosLv - auxiliaryElements.getK() * sinLv;
  986.             final double den     = auxiliaryElements.getB() + 1. + auxiliaryElements.getK() * cosLv + auxiliaryElements.getH() * sinLv;
  987.             return lv + 2. * FastMath.atan(num / den);
  988.         }

  990.         /** Converts eccentric longitude to mean longitude.
  991.          * @param le Eccentric longitude
  992.          * @return Mean longitude
  993.          */
  994.         private double eccentricToMean (final double le) {
  995.             return le - auxiliaryElements.getK() * FastMath.sin(le) + auxiliaryElements.getH() * FastMath.cos(le);
  996.         }

  998.         /** Converts true longitude to mean longitude.
  999.          * @param lv True longitude
  1000.          * @return Eccentric longitude
  1001.          */
  1002.         private double trueToMean (final double lv) {
  1003.             return eccentricToMean(trueToEccentric(lv));
  1004.         }

  1006.         /** Compute δa/δv.
  1007.          *  @param vel satellite velocity
  1008.          *  @return δa/δv
  1009.          */
  1010.         private Vector3D getAoV(final Vector3D vel) {
  1011.             return new Vector3D(context.getTon2a(), vel);
  1012.         }

  1014.         /** Compute δh/δv.
  1015.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  1016.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  1017.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  1018.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  1019.          *  @return δh/δv
  1020.          */
  1021.         private Vector3D getHoV(final double X, final double Y, final double Xdot, final double Ydot) {
  1022.             final double kf = (2. * Xdot * Y - X * Ydot) * context.getOoMU();
  1023.             final double kg = X * Xdot * context.getOoMU();
  1024.             final double kw = auxiliaryElements.getK() * (I * auxiliaryElements.getQ() * Y - auxiliaryElements.getP() * X) * context.getOOAB();
  1025.             return new Vector3D(kf, auxiliaryElements.getVectorF(), -kg, auxiliaryElements.getVectorG(), kw, auxiliaryElements.getVectorW());
  1026.         }

  1028.         /** Compute δk/δv.
  1029.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  1030.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  1031.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  1032.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  1033.          *  @return δk/δv
  1034.          */
  1035.         private Vector3D getKoV(final double X, final double Y, final double Xdot, final double Ydot) {
  1036.             final double kf = Y * Ydot * context.getOoMU();
  1037.             final double kg = (2. * X * Ydot - Xdot * Y) * context.getOoMU();
  1038.             final double kw = auxiliaryElements.getH() * (I * auxiliaryElements.getQ() * Y - auxiliaryElements.getP() * X) * context.getOOAB();
  1039.             return new Vector3D(-kf, auxiliaryElements.getVectorF(), kg, auxiliaryElements.getVectorG(), -kw, auxiliaryElements.getVectorW());
  1040.         }

  1042.         /** Compute δp/δv.
  1043.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  1044.          *  @return δp/δv
  1045.          */
  1046.         private Vector3D getPoV(final double Y) {
  1047.             return new Vector3D(context.getCo2AB() * Y, auxiliaryElements.getVectorW());
  1048.         }

  1050.         /** Compute δq/δv.
  1051.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  1052.          *  @return δq/δv
  1053.          */
  1054.         private Vector3D getQoV(final double X) {
  1055.             return new Vector3D(I * context.getCo2AB() * X, auxiliaryElements.getVectorW());
  1056.         }

  1058.         /** Compute δλ/δv.
  1059.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  1060.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  1061.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  1062.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  1063.          *  @return δλ/δv
  1064.          */
  1065.         private Vector3D getLoV(final double X, final double Y, final double Xdot, final double Ydot) {
  1066.             final Vector3D pos = new Vector3D(X, auxiliaryElements.getVectorF(), Y, auxiliaryElements.getVectorG());
  1067.             final Vector3D v2  = new Vector3D(auxiliaryElements.getK(), getHoV(X, Y, Xdot, Ydot), -auxiliaryElements.getH(), getKoV(X, Y, Xdot, Ydot));
  1068.             return new Vector3D(-2. * context.getOOA(), pos, context.getOoBpo(), v2, (I * auxiliaryElements.getQ() * Y - auxiliaryElements.getP() * X) * context.getOOA(), auxiliaryElements.getVectorW());
  1069.         }

  1071.     }

  1073.     /** Class used to {@link #integrate(UnivariateVectorFunction, double, double) integrate}
  1074.      *  a {@link org.hipparchus.analysis.UnivariateVectorFunction function}
  1075.      *  of the orbital elements using the Gaussian quadrature rule to get the acceleration.
  1076.      */
  1077.     private static class GaussQuadrature {

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

  1081.         /** Points for quadrature of order 12. */
  1082.         private static final double[] P_12 = {
  1083.             -0.98156063424671910000,
  1084.             -0.90411725637047490000,
  1085.             -0.76990267419430470000,
  1086.             -0.58731795428661740000,
  1087.             -0.36783149899818024000,
  1088.             -0.12523340851146890000,
  1089.             0.12523340851146890000,
  1090.             0.36783149899818024000,
  1091.             0.58731795428661740000,
  1092.             0.76990267419430470000,
  1093.             0.90411725637047490000,
  1094.             0.98156063424671910000
  1095.         };

  1097.         /** Weights for quadrature of order 12. */
  1098.         private static final double[] W_12 = {
  1099.             0.04717533638651220000,
  1100.             0.10693932599531830000,
  1101.             0.16007832854334633000,
  1102.             0.20316742672306584000,
  1103.             0.23349253653835478000,
  1104.             0.24914704581340286000,
  1105.             0.24914704581340286000,
  1106.             0.23349253653835478000,
  1107.             0.20316742672306584000,
  1108.             0.16007832854334633000,
  1109.             0.10693932599531830000,
  1110.             0.04717533638651220000
  1111.         };

  1113.         /** Points for quadrature of order 16. */
  1114.         private static final double[] P_16 = {
  1115.             -0.98940093499164990000,
  1116.             -0.94457502307323260000,
  1117.             -0.86563120238783160000,
  1118.             -0.75540440835500310000,
  1119.             -0.61787624440264380000,
  1120.             -0.45801677765722737000,
  1121.             -0.28160355077925890000,
  1122.             -0.09501250983763745000,
  1123.             0.09501250983763745000,
  1124.             0.28160355077925890000,
  1125.             0.45801677765722737000,
  1126.             0.61787624440264380000,
  1127.             0.75540440835500310000,
  1128.             0.86563120238783160000,
  1129.             0.94457502307323260000,
  1130.             0.98940093499164990000
  1131.         };

  1133.         /** Weights for quadrature of order 16. */
  1134.         private static final double[] W_16 = {
  1135.             0.02715245941175405800,
  1136.             0.06225352393864777000,
  1137.             0.09515851168249283000,
  1138.             0.12462897125553388000,
  1139.             0.14959598881657685000,
  1140.             0.16915651939500256000,
  1141.             0.18260341504492360000,
  1142.             0.18945061045506847000,
  1143.             0.18945061045506847000,
  1144.             0.18260341504492360000,
  1145.             0.16915651939500256000,
  1146.             0.14959598881657685000,
  1147.             0.12462897125553388000,
  1148.             0.09515851168249283000,
  1149.             0.06225352393864777000,
  1150.             0.02715245941175405800
  1151.         };

  1153.         /** Points for quadrature of order 20. */
  1154.         private static final double[] P_20 = {
  1155.             -0.99312859918509490000,
  1156.             -0.96397192727791390000,
  1157.             -0.91223442825132600000,
  1158.             -0.83911697182221890000,
  1159.             -0.74633190646015080000,
  1160.             -0.63605368072651510000,
  1161.             -0.51086700195082700000,
  1162.             -0.37370608871541955000,
  1163.             -0.22778585114164507000,
  1164.             -0.07652652113349734000,
  1165.             0.07652652113349734000,
  1166.             0.22778585114164507000,
  1167.             0.37370608871541955000,
  1168.             0.51086700195082700000,
  1169.             0.63605368072651510000,
  1170.             0.74633190646015080000,
  1171.             0.83911697182221890000,
  1172.             0.91223442825132600000,
  1173.             0.96397192727791390000,
  1174.             0.99312859918509490000
  1175.         };

  1177.         /** Weights for quadrature of order 20. */
  1178.         private static final double[] W_20 = {
  1179.             0.01761400713915226400,
  1180.             0.04060142980038684000,
  1181.             0.06267204833410904000,
  1182.             0.08327674157670477000,
  1183.             0.10193011981724048000,
  1184.             0.11819453196151844000,
  1185.             0.13168863844917678000,
  1186.             0.14209610931838212000,
  1187.             0.14917298647260380000,
  1188.             0.15275338713072600000,
  1189.             0.15275338713072600000,
  1190.             0.14917298647260380000,
  1191.             0.14209610931838212000,
  1192.             0.13168863844917678000,
  1193.             0.11819453196151844000,
  1194.             0.10193011981724048000,
  1195.             0.08327674157670477000,
  1196.             0.06267204833410904000,
  1197.             0.04060142980038684000,
  1198.             0.01761400713915226400
  1199.         };

  1201.         /** Points for quadrature of order 24. */
  1202.         private static final double[] P_24 = {
  1203.             -0.99518721999702130000,
  1204.             -0.97472855597130950000,
  1205.             -0.93827455200273270000,
  1206.             -0.88641552700440100000,
  1207.             -0.82000198597390300000,
  1208.             -0.74012419157855440000,
  1209.             -0.64809365193697550000,
  1210.             -0.54542147138883950000,
  1211.             -0.43379350762604520000,
  1212.             -0.31504267969616340000,
  1213.             -0.19111886747361634000,
  1214.             -0.06405689286260563000,
  1215.             0.06405689286260563000,
  1216.             0.19111886747361634000,
  1217.             0.31504267969616340000,
  1218.             0.43379350762604520000,
  1219.             0.54542147138883950000,
  1220.             0.64809365193697550000,
  1221.             0.74012419157855440000,
  1222.             0.82000198597390300000,
  1223.             0.88641552700440100000,
  1224.             0.93827455200273270000,
  1225.             0.97472855597130950000,
  1226.             0.99518721999702130000
  1227.         };

  1229.         /** Weights for quadrature of order 24. */
  1230.         private static final double[] W_24 = {
  1231.             0.01234122979998733500,
  1232.             0.02853138862893380600,
  1233.             0.04427743881741981000,
  1234.             0.05929858491543691500,
  1235.             0.07334648141108027000,
  1236.             0.08619016153195320000,
  1237.             0.09761865210411391000,
  1238.             0.10744427011596558000,
  1239.             0.11550566805372553000,
  1240.             0.12167047292780335000,
  1241.             0.12583745634682825000,
  1242.             0.12793819534675221000,
  1243.             0.12793819534675221000,
  1244.             0.12583745634682825000,
  1245.             0.12167047292780335000,
  1246.             0.11550566805372553000,
  1247.             0.10744427011596558000,
  1248.             0.09761865210411391000,
  1249.             0.08619016153195320000,
  1250.             0.07334648141108027000,
  1251.             0.05929858491543691500,
  1252.             0.04427743881741981000,
  1253.             0.02853138862893380600,
  1254.             0.01234122979998733500
  1255.         };

  1257.         /** Points for quadrature of order 32. */
  1258.         private static final double[] P_32 = {
  1259.             -0.99726386184948160000,
  1260.             -0.98561151154526840000,
  1261.             -0.96476225558750640000,
  1262.             -0.93490607593773970000,
  1263.             -0.89632115576605220000,
  1264.             -0.84936761373256990000,
  1265.             -0.79448379596794250000,
  1266.             -0.73218211874028970000,
  1267.             -0.66304426693021520000,
  1268.             -0.58771575724076230000,
  1269.             -0.50689990893222950000,
  1270.             -0.42135127613063540000,
  1271.             -0.33186860228212767000,
  1272.             -0.23928736225213710000,
  1273.             -0.14447196158279646000,
  1274.             -0.04830766568773831000,
  1275.             0.04830766568773831000,
  1276.             0.14447196158279646000,
  1277.             0.23928736225213710000,
  1278.             0.33186860228212767000,
  1279.             0.42135127613063540000,
  1280.             0.50689990893222950000,
  1281.             0.58771575724076230000,
  1282.             0.66304426693021520000,
  1283.             0.73218211874028970000,
  1284.             0.79448379596794250000,
  1285.             0.84936761373256990000,
  1286.             0.89632115576605220000,
  1287.             0.93490607593773970000,
  1288.             0.96476225558750640000,
  1289.             0.98561151154526840000,
  1290.             0.99726386184948160000
  1291.         };

  1293.         /** Weights for quadrature of order 32. */
  1294.         private static final double[] W_32 = {
  1295.             0.00701861000947013600,
  1296.             0.01627439473090571200,
  1297.             0.02539206530926214200,
  1298.             0.03427386291302141000,
  1299.             0.04283589802222658600,
  1300.             0.05099805926237621600,
  1301.             0.05868409347853559000,
  1302.             0.06582222277636193000,
  1303.             0.07234579410884862000,
  1304.             0.07819389578707042000,
  1305.             0.08331192422694673000,
  1306.             0.08765209300440380000,
  1307.             0.09117387869576390000,
  1308.             0.09384439908080441000,
  1309.             0.09563872007927487000,
  1310.             0.09654008851472784000,
  1311.             0.09654008851472784000,
  1312.             0.09563872007927487000,
  1313.             0.09384439908080441000,
  1314.             0.09117387869576390000,
  1315.             0.08765209300440380000,
  1316.             0.08331192422694673000,
  1317.             0.07819389578707042000,
  1318.             0.07234579410884862000,
  1319.             0.06582222277636193000,
  1320.             0.05868409347853559000,
  1321.             0.05099805926237621600,
  1322.             0.04283589802222658600,
  1323.             0.03427386291302141000,
  1324.             0.02539206530926214200,
  1325.             0.01627439473090571200,
  1326.             0.00701861000947013600
  1327.         };

  1329.         /** Points for quadrature of order 40. */
  1330.         private static final double[] P_40 = {
  1331.             -0.99823770971055930000,
  1332.             -0.99072623869945710000,
  1333.             -0.97725994998377420000,
  1334.             -0.95791681921379170000,
  1335.             -0.93281280827867660000,
  1336.             -0.90209880696887420000,
  1337.             -0.86595950321225960000,
  1338.             -0.82461223083331170000,
  1339.             -0.77830565142651940000,
  1340.             -0.72731825518992710000,
  1341.             -0.67195668461417960000,
  1342.             -0.61255388966798030000,
  1343.             -0.54946712509512820000,
  1344.             -0.48307580168617870000,
  1345.             -0.41377920437160500000,
  1346.             -0.34199409082575850000,
  1347.             -0.26815218500725370000,
  1348.             -0.19269758070137110000,
  1349.             -0.11608407067525522000,
  1350.             -0.03877241750605081600,
  1351.             0.03877241750605081600,
  1352.             0.11608407067525522000,
  1353.             0.19269758070137110000,
  1354.             0.26815218500725370000,
  1355.             0.34199409082575850000,
  1356.             0.41377920437160500000,
  1357.             0.48307580168617870000,
  1358.             0.54946712509512820000,
  1359.             0.61255388966798030000,
  1360.             0.67195668461417960000,
  1361.             0.72731825518992710000,
  1362.             0.77830565142651940000,
  1363.             0.82461223083331170000,
  1364.             0.86595950321225960000,
  1365.             0.90209880696887420000,
  1366.             0.93281280827867660000,
  1367.             0.95791681921379170000,
  1368.             0.97725994998377420000,
  1369.             0.99072623869945710000,
  1370.             0.99823770971055930000
  1371.         };

  1373.         /** Weights for quadrature of order 40. */
  1374.         private static final double[] W_40 = {
  1375.             0.00452127709853309800,
  1376.             0.01049828453115270400,
  1377.             0.01642105838190797300,
  1378.             0.02224584919416689000,
  1379.             0.02793700698002338000,
  1380.             0.03346019528254786500,
  1381.             0.03878216797447199000,
  1382.             0.04387090818567333000,
  1383.             0.04869580763507221000,
  1384.             0.05322784698393679000,
  1385.             0.05743976909939157000,
  1386.             0.06130624249292891000,
  1387.             0.06480401345660108000,
  1388.             0.06791204581523394000,
  1389.             0.07061164739128681000,
  1390.             0.07288658239580408000,
  1391.             0.07472316905796833000,
  1392.             0.07611036190062619000,
  1393.             0.07703981816424793000,
  1394.             0.07750594797842482000,
  1395.             0.07750594797842482000,
  1396.             0.07703981816424793000,
  1397.             0.07611036190062619000,
  1398.             0.07472316905796833000,
  1399.             0.07288658239580408000,
  1400.             0.07061164739128681000,
  1401.             0.06791204581523394000,
  1402.             0.06480401345660108000,
  1403.             0.06130624249292891000,
  1404.             0.05743976909939157000,
  1405.             0.05322784698393679000,
  1406.             0.04869580763507221000,
  1407.             0.04387090818567333000,
  1408.             0.03878216797447199000,
  1409.             0.03346019528254786500,
  1410.             0.02793700698002338000,
  1411.             0.02224584919416689000,
  1412.             0.01642105838190797300,
  1413.             0.01049828453115270400,
  1414.             0.00452127709853309800
  1415.         };

  1417.         /** Points for quadrature of order 48. */
  1418.         private static final double[] P_48 = {
  1419.             -0.99877100725242610000,
  1420.             -0.99353017226635080000,
  1421.             -0.98412458372282700000,
  1422.             -0.97059159254624720000,
  1423.             -0.95298770316043080000,
  1424.             -0.93138669070655440000,
  1425.             -0.90587913671556960000,
  1426.             -0.87657202027424800000,
  1427.             -0.84358826162439350000,
  1428.             -0.80706620402944250000,
  1429.             -0.76715903251574020000,
  1430.             -0.72403413092381470000,
  1431.             -0.67787237963266400000,
  1432.             -0.62886739677651370000,
  1433.             -0.57722472608397270000,
  1434.             -0.52316097472223300000,
  1435.             -0.46690290475095840000,
  1436.             -0.40868648199071680000,
  1437.             -0.34875588629216070000,
  1438.             -0.28736248735545555000,
  1439.             -0.22476379039468908000,
  1440.             -0.16122235606889174000,
  1441.             -0.09700469920946270000,
  1442.             -0.03238017096286937000,
  1443.             0.03238017096286937000,
  1444.             0.09700469920946270000,
  1445.             0.16122235606889174000,
  1446.             0.22476379039468908000,
  1447.             0.28736248735545555000,
  1448.             0.34875588629216070000,
  1449.             0.40868648199071680000,
  1450.             0.46690290475095840000,
  1451.             0.52316097472223300000,
  1452.             0.57722472608397270000,
  1453.             0.62886739677651370000,
  1454.             0.67787237963266400000,
  1455.             0.72403413092381470000,
  1456.             0.76715903251574020000,
  1457.             0.80706620402944250000,
  1458.             0.84358826162439350000,
  1459.             0.87657202027424800000,
  1460.             0.90587913671556960000,
  1461.             0.93138669070655440000,
  1462.             0.95298770316043080000,
  1463.             0.97059159254624720000,
  1464.             0.98412458372282700000,
  1465.             0.99353017226635080000,
  1466.             0.99877100725242610000
  1467.         };

  1469.         /** Weights for quadrature of order 48. */
  1470.         private static final double[] W_48 = {
  1471.             0.00315334605230596250,
  1472.             0.00732755390127620800,
  1473.             0.01147723457923446900,
  1474.             0.01557931572294386600,
  1475.             0.01961616045735556700,
  1476.             0.02357076083932435600,
  1477.             0.02742650970835688000,
  1478.             0.03116722783279807000,
  1479.             0.03477722256477045000,
  1480.             0.03824135106583080600,
  1481.             0.04154508294346483000,
  1482.             0.04467456085669424000,
  1483.             0.04761665849249054000,
  1484.             0.05035903555385448000,
  1485.             0.05289018948519365000,
  1486.             0.05519950369998416500,
  1487.             0.05727729210040315000,
  1488.             0.05911483969839566000,
  1489.             0.06070443916589384000,
  1490.             0.06203942315989268000,
  1491.             0.06311419228625403000,
  1492.             0.06392423858464817000,
  1493.             0.06446616443595010000,
  1494.             0.06473769681268386000,
  1495.             0.06473769681268386000,
  1496.             0.06446616443595010000,
  1497.             0.06392423858464817000,
  1498.             0.06311419228625403000,
  1499.             0.06203942315989268000,
  1500.             0.06070443916589384000,
  1501.             0.05911483969839566000,
  1502.             0.05727729210040315000,
  1503.             0.05519950369998416500,
  1504.             0.05289018948519365000,
  1505.             0.05035903555385448000,
  1506.             0.04761665849249054000,
  1507.             0.04467456085669424000,
  1508.             0.04154508294346483000,
  1509.             0.03824135106583080600,
  1510.             0.03477722256477045000,
  1511.             0.03116722783279807000,
  1512.             0.02742650970835688000,
  1513.             0.02357076083932435600,
  1514.             0.01961616045735556700,
  1515.             0.01557931572294386600,
  1516.             0.01147723457923446900,
  1517.             0.00732755390127620800,
  1518.             0.00315334605230596250
  1519.         };

  1521.         /** Node points. */
  1522.         private final double[] nodePoints;

  1524.         /** Node weights. */
  1525.         private final double[] nodeWeights;

  1527.         /** Number of points. */
  1528.         private final int numberOfPoints;

  1530.         /** Creates a Gauss integrator of the given order.
  1531.          *
  1532.          *  @param numberOfPoints Order of the integration rule.
  1533.          */
  1534.         GaussQuadrature(final int numberOfPoints) {

  1536.             this.numberOfPoints = numberOfPoints;

  1538.             switch(numberOfPoints) {
  1539.                 case 12 :
  1540.                     this.nodePoints  = P_12.clone();
  1541.                     this.nodeWeights = W_12.clone();
  1542.                     break;
  1543.                 case 16 :
  1544.                     this.nodePoints  = P_16.clone();
  1545.                     this.nodeWeights = W_16.clone();
  1546.                     break;
  1547.                 case 20 :
  1548.                     this.nodePoints  = P_20.clone();
  1549.                     this.nodeWeights = W_20.clone();
  1550.                     break;
  1551.                 case 24 :
  1552.                     this.nodePoints  = P_24.clone();
  1553.                     this.nodeWeights = W_24.clone();
  1554.                     break;
  1555.                 case 32 :
  1556.                     this.nodePoints  = P_32.clone();
  1557.                     this.nodeWeights = W_32.clone();
  1558.                     break;
  1559.                 case 40 :
  1560.                     this.nodePoints  = P_40.clone();
  1561.                     this.nodeWeights = W_40.clone();
  1562.                     break;
  1563.                 case 48 :
  1564.                 default :
  1565.                     this.nodePoints  = P_48.clone();
  1566.                     this.nodeWeights = W_48.clone();
  1567.                     break;
  1568.             }

  1570.         }

  1572.         /** Integrates a given function on the given interval.
  1573.         *
  1574.         *  @param f Function to integrate.
  1575.         *  @param lowerBound Lower bound of the integration interval.
  1576.         *  @param upperBound Upper bound of the integration interval.
  1577.         *  @return the integral of the weighted function.
  1578.         */
  1579.         public double[] integrate(final UnivariateVectorFunction f,
  1580.                final double lowerBound, final double upperBound) {

  1582.             final double[] adaptedPoints  = nodePoints.clone();
  1583.             final double[] adaptedWeights = nodeWeights.clone();
  1584.             transform(adaptedPoints, adaptedWeights, lowerBound, upperBound);
  1585.             return basicIntegrate(f, adaptedPoints, adaptedWeights);
  1586.         }

  1588.        /** Integrates a given function on the given interval.
  1589.        *
  1590.        *  @param <T> the type of the field elements
  1591.        *  @param f Function to integrate.
  1592.        *  @param lowerBound Lower bound of the integration interval.
  1593.        *  @param upperBound Upper bound of the integration interval.
  1594.        *  @param field field utilized by default
  1595.        *  @return the integral of the weighted function.
  1596.        */
  1597.         public <T extends RealFieldElement<T>> T[] integrate(final RealFieldUnivariateVectorFunction<T> f,
  1598.                                                              final T lowerBound, final T upperBound,
  1599.                                                              final Field<T> field) {

  1601.             final T zero = field.getZero();

  1603.             final T[] adaptedPoints  = MathArrays.buildArray(field, numberOfPoints);
  1604.             final T[] adaptedWeights = MathArrays.buildArray(field, numberOfPoints);;

  1606.             for (int i = 0; i < numberOfPoints; i++) {
  1607.                 adaptedPoints[i]  = zero.add(nodePoints[i]);
  1608.                 adaptedWeights[i] = zero.add(nodeWeights[i]);
  1609.             }

  1611.             transform(adaptedPoints, adaptedWeights, lowerBound, upperBound);
  1612.             return basicIntegrate(f, adaptedPoints, adaptedWeights, field);
  1613.         }

  1615.         /** Performs a change of variable so that the integration
  1616.          *  can be performed on an arbitrary interval {@code [a, b]}.
  1617.          *  <p>
  1618.          *  It is assumed that the natural interval is {@code [-1, 1]}.
  1619.          *  </p>
  1620.          *
  1621.          * @param points  Points to adapt to the new interval.
  1622.          * @param weights Weights to adapt to the new interval.
  1623.          * @param a Lower bound of the integration interval.
  1624.          * @param b Lower bound of the integration interval.
  1625.          */
  1626.         private void transform(final double[] points, final double[] weights,
  1627.                 final double a, final double b) {
  1628.             // Scaling
  1629.             final double scale = (b - a) / 2;
  1630.             final double shift = a + scale;
  1631.             for (int i = 0; i < points.length; i++) {
  1632.                 points[i]   = points[i] * scale + shift;
  1633.                 weights[i] *= scale;
  1634.             }
  1635.         }

  1637.         /** Performs a change of variable so that the integration
  1638.          *  can be performed on an arbitrary interval {@code [a, b]}.
  1639.          *  <p>
  1640.          *  It is assumed that the natural interval is {@code [-1, 1]}.
  1641.          *  </p>
  1642.          * @param <T> the type of the field elements
  1643.          * @param points  Points to adapt to the new interval.
  1644.          * @param weights Weights to adapt to the new interval.
  1645.          * @param a Lower bound of the integration interval.
  1646.          * @param b Lower bound of the integration interval
  1647.          */
  1648.         private <T extends RealFieldElement<T>> void transform(final T[] points, final T[] weights,
  1649.                 final T a, final T b) {
  1650.             // Scaling
  1651.             final T scale = (b.subtract(a)).divide(2.);
  1652.             final T shift = a.add(scale);
  1653.             for (int i = 0; i < points.length; i++) {
  1654.                 points[i]   = scale.multiply(points[i]).add(shift);
  1655.                 weights[i]  = scale.multiply(weights[i]);
  1656.             }
  1657.         }

  1659.         /** Returns an estimate of the integral of {@code f(x) * w(x)},
  1660.          *  where {@code w} is a weight function that depends on the actual
  1661.          *  flavor of the Gauss integration scheme.
  1662.          *
  1663.          * @param f Function to integrate.
  1664.          * @param points  Nodes.
  1665.          * @param weights Nodes weights.
  1666.          * @return the integral of the weighted function.
  1667.          */
  1668.         private double[] basicIntegrate(final UnivariateVectorFunction f,
  1669.                 final double[] points,
  1670.                 final double[] weights) {
  1671.             double x = points[0];
  1672.             double w = weights[0];
  1673.             double[] v = f.value(x);
  1674.             final double[] y = new double[v.length];
  1675.             for (int j = 0; j < v.length; j++) {
  1676.                 y[j] = w * v[j];
  1677.             }
  1678.             final double[] t = y.clone();
  1679.             final double[] c = new double[v.length];
  1680.             final double[] s = t.clone();
  1681.             for (int i = 1; i < points.length; i++) {
  1682.                 x = points[i];
  1683.                 w = weights[i];
  1684.                 v = f.value(x);
  1685.                 for (int j = 0; j < v.length; j++) {
  1686.                     y[j] = w * v[j] - c[j];
  1687.                     t[j] =  s[j] + y[j];
  1688.                     c[j] = (t[j] - s[j]) - y[j];
  1689.                     s[j] = t[j];
  1690.                 }
  1691.             }
  1692.             return s;
  1693.         }

  1695.         /** Returns an estimate of the integral of {@code f(x) * w(x)},
  1696.          *  where {@code w} is a weight function that depends on the actual
  1697.          *  flavor of the Gauss integration scheme.
  1698.          *
  1699.          * @param <T> the type of the field elements.
  1700.          * @param f Function to integrate.
  1701.          * @param points  Nodes.
  1702.          * @param weights Nodes weight
  1703.          * @param field field utilized by default
  1704.          * @return the integral of the weighted function.
  1705.          */
  1706.         private <T extends RealFieldElement<T>> T[] basicIntegrate(final RealFieldUnivariateVectorFunction<T> f,
  1707.                 final T[] points,
  1708.                 final T[] weights,
  1709.                 final Field<T> field) {

  1711.             T x = points[0];
  1712.             T w = weights[0];
  1713.             T[] v = f.value(x);

  1715.             final T[] y = MathArrays.buildArray(field, v.length);
  1716.             for (int j = 0; j < v.length; j++) {
  1717.                 y[j] = v[j].multiply(w);
  1718.             }
  1719.             final T[] t = y.clone();
  1720.             final T[] c = MathArrays.buildArray(field, v.length);;
  1721.             final T[] s = t.clone();
  1722.             for (int i = 1; i < points.length; i++) {
  1723.                 x = points[i];
  1724.                 w = weights[i];
  1725.                 v = f.value(x);
  1726.                 for (int j = 0; j < v.length; j++) {
  1727.                     y[j] = v[j].multiply(w).subtract(c[j]);
  1728.                     t[j] = y[j].add(s[j]);
  1729.                     c[j] = (t[j].subtract(s[j])).subtract(y[j]);
  1730.                     s[j] = t[j];
  1731.                 }
  1732.             }
  1733.             return s;
  1734.         }

  1736.     }

  1738.     /** Compute the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
  1739.      *  <p>
  1740.      *  Those coefficients are given in Danielson paper by expression 4.4-(6)
  1741.      *  </p>
  1742.      *  @author Petre Bazavan
  1743.      *  @author Lucian Barbulescu
  1744.      */
  1745.     private class FourierCjSjCoefficients {

  1747.         /** Maximum possible value for j. */
  1748.         private final int jMax;

  1750.         /** The C<sub>i</sub><sup>j</sup> coefficients.
  1751.          * <p>
  1752.          * the index i corresponds to the following elements: <br/>
  1753.          * - 0 for a <br>
  1754.          * - 1 for k <br>
  1755.          * - 2 for h <br>
  1756.          * - 3 for q <br>
  1757.          * - 4 for p <br>
  1758.          * - 5 for λ <br>
  1759.          * </p>
  1760.          */
  1761.         private final double[][] cCoef;

  1763.         /** The C<sub>i</sub><sup>j</sup> coefficients.
  1764.          * <p>
  1765.          * the index i corresponds to the following elements: <br/>
  1766.          * - 0 for a <br>
  1767.          * - 1 for k <br>
  1768.          * - 2 for h <br>
  1769.          * - 3 for q <br>
  1770.          * - 4 for p <br>
  1771.          * - 5 for λ <br>
  1772.          * </p>
  1773.          */
  1774.         private final double[][] sCoef;

  1776.         /** Standard constructor.
  1777.          * @param state the current state
  1778.          * @param jMax maximum value for j
  1779.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1780.          * @param parameters values of the force model parameters
  1781.          */
  1782.         FourierCjSjCoefficients(final SpacecraftState state, final int jMax,
  1783.                                 final AuxiliaryElements auxiliaryElements, final double[] parameters) {

  1785.             //Initialise the fields
  1786.             this.jMax = jMax;

  1788.             //Allocate the arrays
  1789.             final int rows = jMax + 1;
  1790.             cCoef = new double[rows][6];
  1791.             sCoef = new double[rows][6];

  1793.             //Compute the coefficients
  1794.             computeCoefficients(state, auxiliaryElements, parameters);
  1795.         }

  1797.         /**
  1798.          * Compute the Fourrier coefficients.
  1799.          * <p>
  1800.          * Only the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients need to be computed
  1801.          * as D<sub>i</sub><sup>m</sup> is always 0.
  1802.          * </p>
  1803.          * @param state the current state
  1804.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1805.          * @param parameters values of the force model parameters
  1806.          */
  1807.         private void computeCoefficients(final SpacecraftState state,
  1808.                                          final AuxiliaryElements auxiliaryElements,
  1809.                                          final double[] parameters) {

  1811.             // Computes the limits for the integral
  1812.             final double[] ll = getLLimits(state, auxiliaryElements);
  1813.             // Computes integrated mean element rates if Llow < Lhigh
  1814.             if (ll[0] < ll[1]) {
  1815.                 //Compute 1 / PI
  1816.                 final double ooPI = 1 / FastMath.PI;

  1818.                 // loop through all values of j
  1819.                 for (int j = 0; j <= jMax; j++) {
  1820.                     final double[] curentCoefficients =
  1821.                             integrator.integrate(new IntegrableFunction(state, false, j, parameters), ll[0], ll[1]);

  1823.                     //divide by PI and set the values for the coefficients
  1824.                     for (int i = 0; i < 6; i++) {
  1825.                         cCoef[j][i] = ooPI * curentCoefficients[i];
  1826.                         sCoef[j][i] = ooPI * curentCoefficients[i + 6];
  1827.                     }
  1828.                 }
  1829.             }
  1830.         }

  1832.         /** Get the coefficient C<sub>i</sub><sup>j</sup>.
  1833.          * @param i i index - corresponds to the required variation
  1834.          * @param j j index
  1835.          * @return the coefficient C<sub>i</sub><sup>j</sup>
  1836.          */
  1837.         public double getCij(final int i, final int j) {
  1838.             return cCoef[j][i];
  1839.         }

  1841.         /** Get the coefficient S<sub>i</sub><sup>j</sup>.
  1842.          * @param i i index - corresponds to the required variation
  1843.          * @param j j index
  1844.          * @return the coefficient S<sub>i</sub><sup>j</sup>
  1845.          */
  1846.         public double getSij(final int i, final int j) {
  1847.             return sCoef[j][i];
  1848.         }
  1849.     }

  1851.     /** Compute the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients with field elements.
  1852.      *  <p>
  1853.      *  Those coefficients are given in Danielson paper by expression 4.4-(6)
  1854.      *  </p>
  1855.      *  @author Petre Bazavan
  1856.      *  @author Lucian Barbulescu
  1857.      */
  1858.     private class FieldFourierCjSjCoefficients <T extends RealFieldElement<T>> {

  1860.         /** Maximum possible value for j. */
  1861.         private final int jMax;

  1863.         /** The C<sub>i</sub><sup>j</sup> coefficients.
  1864.          * <p>
  1865.          * the index i corresponds to the following elements: <br/>
  1866.          * - 0 for a <br>
  1867.          * - 1 for k <br>
  1868.          * - 2 for h <br>
  1869.          * - 3 for q <br>
  1870.          * - 4 for p <br>
  1871.          * - 5 for λ <br>
  1872.          * </p>
  1873.          */
  1874.         private final T[][] cCoef;

  1876.         /** The C<sub>i</sub><sup>j</sup> coefficients.
  1877.          * <p>
  1878.          * the index i corresponds to the following elements: <br/>
  1879.          * - 0 for a <br>
  1880.          * - 1 for k <br>
  1881.          * - 2 for h <br>
  1882.          * - 3 for q <br>
  1883.          * - 4 for p <br>
  1884.          * - 5 for λ <br>
  1885.          * </p>
  1886.          */
  1887.         private final T[][] sCoef;

  1889.         /** Standard constructor.
  1890.          * @param state the current state
  1891.          * @param jMax maximum value for j
  1892.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1893.          * @param parameters values of the force model parameters
  1894.          * @param field field used by default
  1895.          */
  1896.         FieldFourierCjSjCoefficients(final FieldSpacecraftState<T> state, final int jMax,
  1897.                                      final FieldAuxiliaryElements<T> auxiliaryElements,
  1898.                                      final T[] parameters,
  1899.                                      final Field<T> field) {
  1900.             //Initialise the fields
  1901.             this.jMax = jMax;

  1903.             //Allocate the arrays
  1904.             final int rows = jMax + 1;
  1905.             cCoef = MathArrays.buildArray(field, rows, 6);
  1906.             sCoef = MathArrays.buildArray(field, rows, 6);

  1908.             //Compute the coefficients
  1909.             computeCoefficients(state, auxiliaryElements, parameters, field);
  1910.         }

  1912.         /**
  1913.          * Compute the Fourrier coefficients.
  1914.          * <p>
  1915.          * Only the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients need to be computed
  1916.          * as D<sub>i</sub><sup>m</sup> is always 0.
  1917.          * </p>
  1918.          * @param state the current state
  1919.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1920.          * @param parameters values of the force model parameters
  1921.          * @param field field used by default
  1922.          */
  1923.         private void computeCoefficients(final FieldSpacecraftState<T> state,
  1924.                                          final FieldAuxiliaryElements<T> auxiliaryElements,
  1925.                                          final T[] parameters,
  1926.                                          final Field<T> field) {
  1927.             // Zero
  1928.             final T zero = field.getZero();
  1929.             // Computes the limits for the integral
  1930.             final T[] ll = getLLimits(state, auxiliaryElements);
  1931.             // Computes integrated mean element rates if Llow < Lhigh
  1932.             if (ll[0].getReal() < ll[1].getReal()) {
  1933.                 //Compute 1 / PI
  1934.                 final T ooPI = zero.add(FastMath.PI).reciprocal();

  1936.                 // loop through all values of j
  1937.                 for (int j = 0; j <= jMax; j++) {
  1938.                     final T[] curentCoefficients =
  1939.                             integrator.integrate(new FieldIntegrableFunction<>(state, false, j, parameters, field), ll[0], ll[1], field);

  1941.                     //divide by PI and set the values for the coefficients
  1942.                     for (int i = 0; i < 6; i++) {
  1943.                         cCoef[j][i] = curentCoefficients[i].multiply(ooPI);
  1944.                         sCoef[j][i] = curentCoefficients[i + 6].multiply(ooPI);
  1945.                     }
  1946.                 }
  1947.             }
  1948.         }

  1950.         /** Get the coefficient C<sub>i</sub><sup>j</sup>.
  1951.          * @param i i index - corresponds to the required variation
  1952.          * @param j j index
  1953.          * @return the coefficient C<sub>i</sub><sup>j</sup>
  1954.          */
  1955.         public T getCij(final int i, final int j) {
  1956.             return cCoef[j][i];
  1957.         }

  1959.         /** Get the coefficient S<sub>i</sub><sup>j</sup>.
  1960.          * @param i i index - corresponds to the required variation
  1961.          * @param j j index
  1962.          * @return the coefficient S<sub>i</sub><sup>j</sup>
  1963.          */
  1964.         public T getSij(final int i, final int j) {
  1965.             return sCoef[j][i];
  1966.         }
  1967.     }

  1969.     /** This class handles the short periodic coefficients described in Danielson 2.5.3-26.
  1970.      *
  1971.      * <p>
  1972.      * The value of M is 0. Also, since the values of the Fourier coefficient D<sub>i</sub><sup>m</sup> is 0
  1973.      * then the values of the coefficients D<sub>i</sub><sup>m</sup> for m &gt; 2 are also 0.
  1974.      * </p>
  1975.      * @author Petre Bazavan
  1976.      * @author Lucian Barbulescu
  1977.      *
  1978.      */
  1979.     private static class GaussianShortPeriodicCoefficients implements ShortPeriodTerms {

  1981.         /** Maximum value for j index. */
  1982.         private final int jMax;

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

  1987.         /** Prefix for coefficients keys. */
  1988.         private final String coefficientsKeyPrefix;

  1990.         /** All coefficients slots. */
  1991.         private final transient TimeSpanMap<Slot> slots;

  1993.         /** Constructor.
  1994.          *  @param coefficientsKeyPrefix prefix for coefficients keys
  1995.          *  @param jMax maximum value for j index
  1996.          *  @param interpolationPoints number of points used in the interpolation process
  1997.          *  @param slots all coefficients slots
  1998.          */
  1999.         GaussianShortPeriodicCoefficients(final String coefficientsKeyPrefix,
  2000.                                           final int jMax, final int interpolationPoints,
  2001.                                           final TimeSpanMap<Slot> slots) {
  2002.             //Initialize fields
  2003.             this.jMax                  = jMax;
  2004.             this.interpolationPoints   = interpolationPoints;
  2005.             this.coefficientsKeyPrefix = coefficientsKeyPrefix;
  2006.             this.slots                 = slots;
  2007.         }

  2009.         /** Get the slot valid for some date.
  2010.          * @param meanStates mean states defining the slot
  2011.          * @return slot valid at the specified date
  2012.          */
  2013.         public Slot createSlot(final SpacecraftState... meanStates) {
  2014.             final Slot         slot  = new Slot(jMax, interpolationPoints);
  2015.             final AbsoluteDate first = meanStates[0].getDate();
  2016.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  2017.             if (first.compareTo(last) <= 0) {
  2018.                 slots.addValidAfter(slot, first);
  2019.             } else {
  2020.                 slots.addValidBefore(slot, first);
  2021.             }
  2022.             return slot;
  2023.         }

  2025.         /** Compute the short periodic coefficients.
  2026.          *
  2027.          * @param state current state information: date, kinematics, attitude
  2028.          * @param slot coefficients slot
  2029.          * @param fourierCjSj Fourier coefficients
  2030.          * @param uijvij U and V coefficients
  2031.          * @param n Keplerian mean motion
  2032.          * @param a semi major axis
  2033.          */
  2034.         private void computeCoefficients(final SpacecraftState state, final Slot slot,
  2035.                                          final FourierCjSjCoefficients fourierCjSj,
  2036.                                          final UijVijCoefficients uijvij,
  2037.                                          final double n, final double a) {

  2039.             // get the current date
  2040.             final AbsoluteDate date = state.getDate();

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

  2045.             // 1. / n
  2046.             final double oon = 1. / n;
  2047.             // 3. / (2 * a * n)
  2048.             final double to2an = 1.5 * oon / a;
  2049.             // 3. / (4 * a * n)
  2050.             final double to4an = to2an / 2;

  2052.             // Compute the coefficients for each element
  2053.             final int size = jMax + 1;
  2054.             final double[]   di1        = new double[6];
  2055.             final double[]   di2        = new double[6];
  2056.             final double[][] currentCij = new double[size][6];
  2057.             final double[][] currentSij = new double[size][6];
  2058.             for (int i = 0; i < 6; i++) {

  2060.                 // compute D<sub>i</sub>¹ and D<sub>i</sub>² (all others are 0)
  2061.                 di1[i] = -oon * fourierCjSj.getCij(i, 0);
  2062.                 if (i == 5) {
  2063.                     di1[i] += to2an * uijvij.getU1(0, 0);
  2064.                 }
  2065.                 di2[i] = 0.;
  2066.                 if (i == 5) {
  2067.                     di2[i] += -to4an * fourierCjSj.getCij(0, 0);
  2068.                 }

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

  2073.                 for (int j = 1; j <= jMax; j++) {
  2074.                     // compute the current C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup>
  2075.                     currentCij[j][i] = oon * uijvij.getU1(j, i);
  2076.                     if (i == 5) {
  2077.                         currentCij[j][i] += -to2an * uijvij.getU2(j);
  2078.                     }
  2079.                     currentSij[j][i] = oon * uijvij.getV1(j, i);
  2080.                     if (i == 5) {
  2081.                         currentSij[j][i] += -to2an * uijvij.getV2(j);
  2082.                     }

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

  2088.             }

  2090.             // add the values to the interpolators
  2091.             slot.cij[0].addGridPoint(date, currentCij[0]);
  2092.             slot.dij[1].addGridPoint(date, di1);
  2093.             slot.dij[2].addGridPoint(date, di2);
  2094.             for (int j = 1; j <= jMax; j++) {
  2095.                 slot.cij[j].addGridPoint(date, currentCij[j]);
  2096.                 slot.sij[j].addGridPoint(date, currentSij[j]);
  2097.             }

  2099.         }

  2101.         /** Compute the coefficient k₂⁰ by using the equation
  2102.          * 2.5.3-(9a) from Danielson.
  2103.          * <p>
  2104.          * After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:<br>
  2105.          * k₂⁰ = &Sigma;<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
  2106.          * </p>
  2107.          * @param kMax max value fot k index
  2108.          * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
  2109.          * @return the coefficient k₂⁰
  2110.          */
  2111.         private double computeK20(final int kMax, final double[][] currentRhoSigmaj) {
  2112.             double k20 = 0.;

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

  2120.                 //multiply by 2 / k²
  2121.                 currentTerm *= 2. / (kIndex * kIndex);

  2123.                 // add the term to the result
  2124.                 k20 += currentTerm;
  2125.             }

  2127.             return k20;
  2128.         }

  2130.         /** {@inheritDoc} */
  2131.         @Override
  2132.         public double[] value(final Orbit meanOrbit) {

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

  2137.             // Get the True longitude L
  2138.             final double L = meanOrbit.getLv();

  2140.             // Compute the center (l - λ)
  2141.             final double center =  L - meanOrbit.getLM();
  2142.             // Compute (l - λ)²
  2143.             final double center2 = center * center;

  2145.             // Initialize short periodic variations
  2146.             final double[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  2147.             final double[] d1 = slot.dij[1].value(meanOrbit.getDate());
  2148.             final double[] d2 = slot.dij[2].value(meanOrbit.getDate());
  2149.             for (int i = 0; i < 6; i++) {
  2150.                 shortPeriodicVariation[i] += center * d1[i] + center2 * d2[i];
  2151.             }

  2153.             for (int j = 1; j <= JMAX; j++) {
  2154.                 final double[] c = slot.cij[j].value(meanOrbit.getDate());
  2155.                 final double[] s = slot.sij[j].value(meanOrbit.getDate());
  2156.                 final double cos = FastMath.cos(j * L);
  2157.                 final double sin = FastMath.sin(j * L);
  2158.                 for (int i = 0; i < 6; i++) {
  2159.                     // add corresponding term to the short periodic variation
  2160.                     shortPeriodicVariation[i] += c[i] * cos;
  2161.                     shortPeriodicVariation[i] += s[i] * sin;
  2162.                 }
  2163.             }

  2165.             return shortPeriodicVariation;

  2167.         }

  2169.         /** {@inheritDoc} */
  2170.         public String getCoefficientsKeyPrefix() {
  2171.             return coefficientsKeyPrefix;
  2172.         }

  2174.         /** {@inheritDoc}
  2175.          * <p>
  2176.          * For Gaussian forces, there are JMAX cj coefficients,
  2177.          * JMAX sj coefficients and 3 dj coefficients. As JMAX = 12,
  2178.          * this sums up to 27 coefficients. The j index is the integer
  2179.          * multiplier for the true longitude argument in the cj and sj
  2180.          * coefficients and to the degree in  the polynomial dj coefficients.
  2181.          * </p>
  2182.          */
  2183.         @Override
  2184.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected) {

  2186.             // select the coefficients slot
  2187.             final Slot slot = slots.get(date);

  2189.             final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * JMAX + 3);
  2190.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  2191.             storeIfSelected(coefficients, selected, slot.dij[1].value(date), "d", 1);
  2192.             storeIfSelected(coefficients, selected, slot.dij[2].value(date), "d", 2);
  2193.             for (int j = 1; j <= JMAX; j++) {
  2194.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  2195.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  2196.             }

  2198.             return coefficients;

  2200.         }

  2202.         /** Put a coefficient in a map if selected.
  2203.          * @param map map to populate
  2204.          * @param selected set of coefficients that should be put in the map
  2205.          * (empty set means all coefficients are selected)
  2206.          * @param value coefficient value
  2207.          * @param id coefficient identifier
  2208.          * @param indices list of coefficient indices
  2209.          */
  2210.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
  2211.                                      final double[] value, final String id, final int... indices) {
  2212.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  2213.             keyBuilder.append(id);
  2214.             for (int index : indices) {
  2215.                 keyBuilder.append('[').append(index).append(']');
  2216.             }
  2217.             final String key = keyBuilder.toString();
  2218.             if (selected.isEmpty() || selected.contains(key)) {
  2219.                 map.put(key, value);
  2220.             }
  2221.         }

  2223.     }

  2225.      /** This class handles the short periodic coefficients described in Danielson 2.5.3-26.
  2226.      *
  2227.      * <p>
  2228.      * The value of M is 0. Also, since the values of the Fourier coefficient D<sub>i</sub><sup>m</sup> is 0
  2229.      * then the values of the coefficients D<sub>i</sub><sup>m</sup> for m &gt; 2 are also 0.
  2230.      * </p>
  2231.      * @author Petre Bazavan
  2232.      * @author Lucian Barbulescu
  2233.      *
  2234.      */
  2235.     private static class FieldGaussianShortPeriodicCoefficients <T extends RealFieldElement<T>> implements FieldShortPeriodTerms<T> {

  2237.         /** Maximum value for j index. */
  2238.         private final int jMax;

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

  2243.         /** Prefix for coefficients keys. */
  2244.         private final String coefficientsKeyPrefix;

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

  2249.         /** Constructor.
  2250.          *  @param coefficientsKeyPrefix prefix for coefficients keys
  2251.          *  @param jMax maximum value for j index
  2252.          *  @param interpolationPoints number of points used in the interpolation process
  2253.          *  @param slots all coefficients slots
  2254.          */
  2255.         FieldGaussianShortPeriodicCoefficients(final String coefficientsKeyPrefix,
  2256.                                                final int jMax, final int interpolationPoints,
  2257.                                                final FieldTimeSpanMap<FieldSlot<T>, T> slots) {
  2258.             //Initialize fields
  2259.             this.jMax                  = jMax;
  2260.             this.interpolationPoints   = interpolationPoints;
  2261.             this.coefficientsKeyPrefix = coefficientsKeyPrefix;
  2262.             this.slots                 = slots;
  2263.         }

  2265.         /** Get the slot valid for some date.
  2266.          * @param meanStates mean states defining the slot
  2267.          * @return slot valid at the specified date
  2268.          */
  2269.         @SuppressWarnings("unchecked")
  2270.         public FieldSlot<T> createSlot(final FieldSpacecraftState<T>... meanStates) {
  2271.             final FieldSlot<T>         slot  = new FieldSlot<>(jMax, interpolationPoints);
  2272.             final FieldAbsoluteDate<T> first = meanStates[0].getDate();
  2273.             final FieldAbsoluteDate<T> last  = meanStates[meanStates.length - 1].getDate();
  2274.             if (first.compareTo(last) <= 0) {
  2275.                 slots.addValidAfter(slot, first);
  2276.             } else {
  2277.                 slots.addValidBefore(slot, first);
  2278.             }
  2279.             return slot;
  2280.         }

  2282.         /** Compute the short periodic coefficients.
  2283.          *
  2284.          * @param state current state information: date, kinematics, attitude
  2285.          * @param slot coefficients slot
  2286.          * @param fourierCjSj Fourier coefficients
  2287.          * @param uijvij U and V coefficients
  2288.          * @param n Keplerian mean motion
  2289.          * @param a semi major axis
  2290.          * @param field field used by default
  2291.          */
  2292.         private void computeCoefficients(final FieldSpacecraftState<T> state, final FieldSlot<T> slot,
  2293.                                          final FieldFourierCjSjCoefficients<T> fourierCjSj,
  2294.                                          final FieldUijVijCoefficients<T> uijvij,
  2295.                                          final T n, final T a,
  2296.                                          final Field<T> field) {

  2298.             // Zero
  2299.             final T zero = field.getZero();

  2301.             // get the current date
  2302.             final FieldAbsoluteDate<T> date = state.getDate();

  2304.             // compute the k₂⁰ coefficient
  2305.             final T k20 = computeK20(jMax, uijvij.currentRhoSigmaj, field);

  2307.             // 1. / n
  2308.             final T oon = n.reciprocal();
  2309.             // 3. / (2 * a * n)
  2310.             final T to2an = oon.multiply(1.5).divide(a);
  2311.             // 3. / (4 * a * n)
  2312.             final T to4an = to2an.divide(2.);

  2314.             // Compute the coefficients for each element
  2315.             final int size = jMax + 1;
  2316.             final T[]   di1        = MathArrays.buildArray(field, 6);
  2317.             final T[]   di2        = MathArrays.buildArray(field, 6);
  2318.             final T[][] currentCij = MathArrays.buildArray(field, size, 6);
  2319.             final T[][] currentSij = MathArrays.buildArray(field, size, 6);
  2320.             for (int i = 0; i < 6; i++) {

  2322.                 // compute D<sub>i</sub>¹ and D<sub>i</sub>² (all others are 0)
  2323.                 di1[i] = oon.negate().multiply(fourierCjSj.getCij(i, 0));
  2324.                 if (i == 5) {
  2325.                     di1[i] = di1[i].add(to2an.multiply(uijvij.getU1(0, 0)));
  2326.                 }
  2327.                 di2[i] = zero;
  2328.                 if (i == 5) {
  2329.                     di2[i] = di2[i].add(to4an.negate().multiply(fourierCjSj.getCij(0, 0)));
  2330.                 }

  2332.                 //the C<sub>i</sub>⁰ is computed based on all others
  2333.                 currentCij[0][i] = di2[i].negate().multiply(k20);

  2335.                 for (int j = 1; j <= jMax; j++) {
  2336.                     // compute the current C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup>
  2337.                     currentCij[j][i] = oon.multiply(uijvij.getU1(j, i));
  2338.                     if (i == 5) {
  2339.                         currentCij[j][i] = currentCij[j][i].add(to2an.negate().multiply(uijvij.getU2(j)));
  2340.                     }
  2341.                     currentSij[j][i] = oon.multiply(uijvij.getV1(j, i));
  2342.                     if (i == 5) {
  2343.                         currentSij[j][i] = currentSij[j][i].add(to2an.negate().multiply(uijvij.getV2(j)));
  2344.                     }

  2346.                     // add the computed coefficients to C<sub>i</sub>⁰
  2347.                     currentCij[0][i] = currentCij[0][i].add(currentCij[j][i].multiply(uijvij.currentRhoSigmaj[0][j]).add(currentSij[j][i].multiply(uijvij.currentRhoSigmaj[1][j])).negate());
  2348.                 }

  2350.             }

  2352.             // add the values to the interpolators
  2353.             slot.cij[0].addGridPoint(date, currentCij[0]);
  2354.             slot.dij[1].addGridPoint(date, di1);
  2355.             slot.dij[2].addGridPoint(date, di2);
  2356.             for (int j = 1; j <= jMax; j++) {
  2357.                 slot.cij[j].addGridPoint(date, currentCij[j]);
  2358.                 slot.sij[j].addGridPoint(date, currentSij[j]);
  2359.             }

  2361.         }

  2363.         /** Compute the coefficient k₂⁰ by using the equation
  2364.          * 2.5.3-(9a) from Danielson.
  2365.          * <p>
  2366.          * After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:<br>
  2367.          * k₂⁰ = &Sigma;<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
  2368.          * </p>
  2369.          * @param kMax max value fot k index
  2370.          * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
  2371.          * @param field field used by default
  2372.          * @return the coefficient k₂⁰
  2373.          */
  2374.         private T computeK20(final int kMax, final T[][] currentRhoSigmaj, final Field<T> field) {
  2375.             final T zero = field.getZero();
  2376.             T k20 = zero;

  2378.             for (int kIndex = 1; kIndex <= kMax; kIndex++) {
  2379.                 // After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:
  2380.                 //k₂⁰ = &Sigma;<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
  2381.                 T currentTerm = currentRhoSigmaj[1][kIndex].multiply(currentRhoSigmaj[1][kIndex]).
  2382.                                 add(currentRhoSigmaj[0][kIndex].multiply(currentRhoSigmaj[0][kIndex]));

  2384.                 //multiply by 2 / k²
  2385.                 currentTerm = currentTerm.multiply(2. / (kIndex * kIndex));

  2387.                 // add the term to the result
  2388.                 k20 = k20.add(currentTerm);
  2389.             }

  2391.             return k20;
  2392.         }

  2394.         /** {@inheritDoc} */
  2395.         @Override
  2396.         public T[] value(final FieldOrbit<T> meanOrbit) {

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

  2401.             // Get the True longitude L
  2402.             final T L = meanOrbit.getLv();

  2404.             // Compute the center (l - λ)
  2405.             final T center =  L.subtract(meanOrbit.getLM());
  2406.             // Compute (l - λ)²
  2407.             final T center2 = center.multiply(center);

  2409.             // Initialize short periodic variations
  2410.             final T[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  2411.             final T[] d1 = slot.dij[1].value(meanOrbit.getDate());
  2412.             final T[] d2 = slot.dij[2].value(meanOrbit.getDate());
  2413.             for (int i = 0; i < 6; i++) {
  2414.                 shortPeriodicVariation[i] = shortPeriodicVariation[i].add(center.multiply(d1[i]).add(center2.multiply(d2[i])));
  2415.             }

  2417.             for (int j = 1; j <= JMAX; j++) {
  2418.                 final T[] c = slot.cij[j].value(meanOrbit.getDate());
  2419.                 final T[] s = slot.sij[j].value(meanOrbit.getDate());
  2420.                 final T cos = FastMath.cos(L.multiply(j));
  2421.                 final T sin = FastMath.sin(L.multiply(j));
  2422.                 for (int i = 0; i < 6; i++) {
  2423.                     // add corresponding term to the short periodic variation
  2424.                     shortPeriodicVariation[i] = shortPeriodicVariation[i].add(c[i].multiply(cos));
  2425.                     shortPeriodicVariation[i] = shortPeriodicVariation[i].add(s[i].multiply(sin));
  2426.                 }
  2427.             }

  2429.             return shortPeriodicVariation;

  2431.         }

  2433.         /** {@inheritDoc} */
  2434.         public String getCoefficientsKeyPrefix() {
  2435.             return coefficientsKeyPrefix;
  2436.         }

  2438.         /** {@inheritDoc}
  2439.          * <p>
  2440.          * For Gaussian forces, there are JMAX cj coefficients,
  2441.          * JMAX sj coefficients and 3 dj coefficients. As JMAX = 12,
  2442.          * this sums up to 27 coefficients. The j index is the integer
  2443.          * multiplier for the true longitude argument in the cj and sj
  2444.          * coefficients and to the degree in  the polynomial dj coefficients.
  2445.          * </p>
  2446.          */
  2447.         @Override
  2448.         public Map<String, T[]> getCoefficients(final FieldAbsoluteDate<T> date, final Set<String> selected) {

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

  2453.             final Map<String, T[]> coefficients = new HashMap<String, T[]>(2 * JMAX + 3);
  2454.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  2455.             storeIfSelected(coefficients, selected, slot.dij[1].value(date), "d", 1);
  2456.             storeIfSelected(coefficients, selected, slot.dij[2].value(date), "d", 2);
  2457.             for (int j = 1; j <= JMAX; j++) {
  2458.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  2459.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  2460.             }

  2462.             return coefficients;

  2464.         }

  2466.         /** Put a coefficient in a map if selected.
  2467.          * @param map map to populate
  2468.          * @param selected set of coefficients that should be put in the map
  2469.          * (empty set means all coefficients are selected)
  2470.          * @param value coefficient value
  2471.          * @param id coefficient identifier
  2472.          * @param indices list of coefficient indices
  2473.          */
  2474.         private void storeIfSelected(final Map<String, T[]> map, final Set<String> selected,
  2475.                                      final T[] value, final String id, final int... indices) {
  2476.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  2477.             keyBuilder.append(id);
  2478.             for (int index : indices) {
  2479.                 keyBuilder.append('[').append(index).append(']');
  2480.             }
  2481.             final String key = keyBuilder.toString();
  2482.             if (selected.isEmpty() || selected.contains(key)) {
  2483.                 map.put(key, value);
  2484.             }
  2485.         }

  2487.     }


  2490.     /** The U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients described by
  2491.      * equations 2.5.3-(21) and 2.5.3-(22) from Danielson.
  2492.      * <p>
  2493.      * The index i takes only the values 1 and 2<br>
  2494.      * For U only the index 0 for j is used.
  2495.      * </p>
  2496.      *
  2497.      * @author Petre Bazavan
  2498.      * @author Lucian Barbulescu
  2499.      */
  2500.     private static class UijVijCoefficients {

  2502.         /** The U₁<sup>j</sup> coefficients.
  2503.          * <p>
  2504.          * The first index identifies the Fourier coefficients used<br>
  2505.          * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
  2506.          * The only exception is when j = 0 when only the coefficient for fourier index = 1 (i == 0) is needed.<br>
  2507.          * Also, for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
  2508.          * to compute the coefficients U₂<sup>j</sup>
  2509.          * </p>
  2510.          */
  2511.         private final double[][] u1ij;

  2513.         /** The V₁<sup>j</sup> coefficients.
  2514.          * <p>
  2515.          * The first index identifies the Fourier coefficients used<br>
  2516.          * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
  2517.          * for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
  2518.          * to compute the coefficients V₂<sup>j</sup>
  2519.          * </p>
  2520.          */
  2521.         private final double[][] v1ij;

  2523.         /** The U₂<sup>j</sup> coefficients.
  2524.          * <p>
  2525.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
  2526.          * </p>
  2527.          */
  2528.         private final double[] u2ij;

  2530.         /** The V₂<sup>j</sup> coefficients.
  2531.          * <p>
  2532.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
  2533.          * </p>
  2534.          */
  2535.         private final double[] v2ij;

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

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

  2543.         /** The maximum value for j index. */
  2544.         private final int jMax;

  2546.         /** Constructor.
  2547.          * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
  2548.          * @param fourierCjSj the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  2549.          * @param jMax maximum value for j index
  2550.          */
  2551.         UijVijCoefficients(final double[][] currentRhoSigmaj, final FourierCjSjCoefficients fourierCjSj, final int jMax) {
  2552.             this.currentRhoSigmaj = currentRhoSigmaj;
  2553.             this.fourierCjSj = fourierCjSj;
  2554.             this.jMax = jMax;

  2556.             // initialize the internal arrays.
  2557.             this.u1ij = new double[6][2 * jMax + 1];
  2558.             this.v1ij = new double[6][2 * jMax + 1];
  2559.             this.u2ij = new double[jMax + 1];
  2560.             this.v2ij = new double[jMax + 1];

  2562.             //compute the coefficients
  2563.             computeU1V1Coefficients();
  2564.             computeU2V2Coefficients();
  2565.         }

  2567.         /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients. */
  2568.         private void computeU1V1Coefficients() {
  2569.             // generate the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients
  2570.             // for j >= 1
  2571.             // also the U₁⁰ for Fourier index = 1 (i == 0) coefficient will be computed
  2572.             u1ij[0][0] = 0;
  2573.             for (int j = 1; j <= jMax; j++) {
  2574.                 // compute 1 / j
  2575.                 final double ooj = 1. / j;

  2577.                 for (int i = 0; i < 6; i++) {
  2578.                     //j is aready between 1 and J
  2579.                     u1ij[i][j] = fourierCjSj.getSij(i, j);
  2580.                     v1ij[i][j] = fourierCjSj.getCij(i, j);

  2582.                     // 1 - δ<sub>1j</sub> is 1 for all j > 1
  2583.                     if (j > 1) {
  2584.                         // k starts with 1 because j-J is less than or equal to 0
  2585.                         for (int kIndex = 1; kIndex <= j - 1; kIndex++) {
  2586.                             // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
  2587.                             // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
  2588.                             u1ij[i][j] +=   fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[1][kIndex] +
  2589.                                             fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[0][kIndex];

  2591.                             // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  2592.                             // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  2593.                             v1ij[i][j] +=   fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[0][kIndex] -
  2594.                                             fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[1][kIndex];
  2595.                         }
  2596.                     }

  2598.                     // since j must be between 1 and J-1 and is already between 1 and J
  2599.                     // the following sum is skiped only for j = jMax
  2600.                     if (j != jMax) {
  2601.                         for (int kIndex = 1; kIndex <= jMax - j; kIndex++) {
  2602.                             // -C<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub> +
  2603.                             // S<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub>
  2604.                             u1ij[i][j] +=   -fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[1][kIndex] +
  2605.                                             fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[0][kIndex];

  2607.                             // C<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub> +
  2608.                             // S<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub>
  2609.                             v1ij[i][j] +=   fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[0][kIndex] +
  2610.                                             fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[1][kIndex];
  2611.                         }
  2612.                     }

  2614.                     for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  2615.                         // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  2616.                         // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  2617.                         u1ij[i][j] +=   -fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[1][j + kIndex] -
  2618.                                         fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[0][j + kIndex];

  2620.                         // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  2621.                         // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  2622.                         v1ij[i][j] +=   fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[0][j + kIndex] +
  2623.                                         fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[1][j + kIndex];
  2624.                     }

  2626.                     // divide by 1 / j
  2627.                     u1ij[i][j] *= -ooj;
  2628.                     v1ij[i][j] *= ooj;

  2630.                     // if index = 1 (i == 0) add the computed terms to U₁⁰
  2631.                     if (i == 0) {
  2632.                         //- (U₁<sup>j</sup> * ρ<sub>j</sub> + V₁<sup>j</sup> * σ<sub>j</sub>
  2633.                         u1ij[0][0] += -u1ij[0][j] * currentRhoSigmaj[0][j] - v1ij[0][j] * currentRhoSigmaj[1][j];
  2634.                     }
  2635.                 }
  2636.             }

  2638.             // Terms with j > jMax are required only when computing the coefficients
  2639.             // U₂<sup>j</sup> and V₂<sup>j</sup>
  2640.             // and those coefficients are only required for Fourier index = 1 (i == 0).
  2641.             for (int j = jMax + 1; j <= 2 * jMax; j++) {
  2642.                 // compute 1 / j
  2643.                 final double ooj = 1. / j;
  2644.                 //the value of i is 0
  2645.                 u1ij[0][j] = 0.;
  2646.                 v1ij[0][j] = 0.;

  2648.                 //k starts from j-J as it is always greater than or equal to 1
  2649.                 for (int kIndex = j - jMax; kIndex <= j - 1; kIndex++) {
  2650.                     // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
  2651.                     // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
  2652.                     u1ij[0][j] +=   fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[1][kIndex] +
  2653.                                     fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[0][kIndex];

  2655.                     // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  2656.                     // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  2657.                     v1ij[0][j] +=   fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[0][kIndex] -
  2658.                                     fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[1][kIndex];
  2659.                 }
  2660.                 for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  2661.                     // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  2662.                     // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  2663.                     u1ij[0][j] +=   -fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[1][j + kIndex] -
  2664.                                     fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[0][j + kIndex];

  2666.                     // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  2667.                     // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  2668.                     v1ij[0][j] +=   fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[0][j + kIndex] +
  2669.                                     fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[1][j + kIndex];
  2670.                 }

  2672.                 // divide by 1 / j
  2673.                 u1ij[0][j] *= -ooj;
  2674.                 v1ij[0][j] *= ooj;
  2675.             }
  2676.         }

  2678.         /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
  2679.          * <p>
  2680.          * Only the coefficients for Fourier index = 1 (i == 0) are required.
  2681.          * </p>
  2682.          */
  2683.         private void computeU2V2Coefficients() {
  2684.             for (int j = 1; j <= jMax; j++) {
  2685.                 // compute 1 / j
  2686.                 final double ooj = 1. / j;

  2688.                 // only the values for i == 0 are computed
  2689.                 u2ij[j] = v1ij[0][j];
  2690.                 v2ij[j] = u1ij[0][j];

  2692.                 // 1 - δ<sub>1j</sub> is 1 for all j > 1
  2693.                 if (j > 1) {
  2694.                     for (int l = 1; l <= j - 1; l++) {
  2695.                         // U₁<sup>j-l</sup> * σ<sub>l</sub> +
  2696.                         // V₁<sup>j-l</sup> * ρ<sub>l</sub>
  2697.                         u2ij[j] +=   u1ij[0][j - l] * currentRhoSigmaj[1][l] +
  2698.                                      v1ij[0][j - l] * currentRhoSigmaj[0][l];

  2700.                         // U₁<sup>j-l</sup> * ρ<sub>l</sub> -
  2701.                         // V₁<sup>j-l</sup> * σ<sub>l</sub>
  2702.                         v2ij[j] +=   u1ij[0][j - l] * currentRhoSigmaj[0][l] -
  2703.                                      v1ij[0][j - l] * currentRhoSigmaj[1][l];
  2704.                     }
  2705.                 }

  2707.                 for (int l = 1; l <= jMax; l++) {
  2708.                     // -U₁<sup>j+l</sup> * σ<sub>l</sub> +
  2709.                     // U₁<sup>l</sup> * σ<sub>j+l</sub> +
  2710.                     // V₁<sup>j+l</sup> * ρ<sub>l</sub> -
  2711.                     // V₁<sup>l</sup> * ρ<sub>j+l</sub>
  2712.                     u2ij[j] +=   -u1ij[0][j + l] * currentRhoSigmaj[1][l] +
  2713.                                   u1ij[0][l] * currentRhoSigmaj[1][j + l] +
  2714.                                   v1ij[0][j + l] * currentRhoSigmaj[0][l] -
  2715.                                   v1ij[0][l] * currentRhoSigmaj[0][j + l];

  2717.                     // U₁<sup>j+l</sup> * ρ<sub>l</sub> +
  2718.                     // U₁<sup>l</sup> * ρ<sub>j+l</sub> +
  2719.                     // V₁<sup>j+l</sup> * σ<sub>l</sub> +
  2720.                     // V₁<sup>l</sup> * σ<sub>j+l</sub>
  2721.                     u2ij[j] +=   u1ij[0][j + l] * currentRhoSigmaj[0][l] +
  2722.                                  u1ij[0][l] * currentRhoSigmaj[0][j + l] +
  2723.                                  v1ij[0][j + l] * currentRhoSigmaj[1][l] +
  2724.                                  v1ij[0][l] * currentRhoSigmaj[1][j + l];
  2725.                 }

  2727.                 // divide by 1 / j
  2728.                 u2ij[j] *= -ooj;
  2729.                 v2ij[j] *= ooj;
  2730.             }
  2731.         }

  2733.         /** Get the coefficient U₁<sup>j</sup> for Fourier index i.
  2734.          *
  2735.          * @param j j index
  2736.          * @param i Fourier index (starts at 0)
  2737.          * @return the coefficient U₁<sup>j</sup> for the given Fourier index i
  2738.          */
  2739.         public double getU1(final int j, final int i) {
  2740.             return u1ij[i][j];
  2741.         }

  2743.         /** Get the coefficient V₁<sup>j</sup> for Fourier index i.
  2744.          *
  2745.          * @param j j index
  2746.          * @param i Fourier index (starts at 0)
  2747.          * @return the coefficient V₁<sup>j</sup> for the given Fourier index i
  2748.          */
  2749.         public double getV1(final int j, final int i) {
  2750.             return v1ij[i][j];
  2751.         }

  2753.         /** Get the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0).
  2754.          *
  2755.          * @param j j index
  2756.          * @return the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0)
  2757.          */
  2758.         public double getU2(final int j) {
  2759.             return u2ij[j];
  2760.         }

  2762.         /** Get the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0).
  2763.          *
  2764.          * @param j j index
  2765.          * @return the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0)
  2766.          */
  2767.         public double getV2(final int j) {
  2768.             return v2ij[j];
  2769.         }
  2770.     }

  2772.     /** The U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients described by
  2773.      * equations 2.5.3-(21) and 2.5.3-(22) from Danielson.
  2774.      * <p>
  2775.      * The index i takes only the values 1 and 2<br>
  2776.      * For U only the index 0 for j is used.
  2777.      * </p>
  2778.      *
  2779.      * @author Petre Bazavan
  2780.      * @author Lucian Barbulescu
  2781.      */
  2782.     private static class FieldUijVijCoefficients <T extends RealFieldElement<T>> {

  2784.         /** The U₁<sup>j</sup> coefficients.
  2785.          * <p>
  2786.          * The first index identifies the Fourier coefficients used<br>
  2787.          * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
  2788.          * The only exception is when j = 0 when only the coefficient for fourier index = 1 (i == 0) is needed.<br>
  2789.          * Also, for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
  2790.          * to compute the coefficients U₂<sup>j</sup>
  2791.          * </p>
  2792.          */
  2793.         private final T[][] u1ij;

  2795.         /** The V₁<sup>j</sup> coefficients.
  2796.          * <p>
  2797.          * The first index identifies the Fourier coefficients used<br>
  2798.          * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
  2799.          * for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
  2800.          * to compute the coefficients V₂<sup>j</sup>
  2801.          * </p>
  2802.          */
  2803.         private final T[][] v1ij;

  2805.         /** The U₂<sup>j</sup> coefficients.
  2806.          * <p>
  2807.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
  2808.          * </p>
  2809.          */
  2810.         private final T[] u2ij;

  2812.         /** The V₂<sup>j</sup> coefficients.
  2813.          * <p>
  2814.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
  2815.          * </p>
  2816.          */
  2817.         private final T[] v2ij;

  2819.         /** The current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients. */
  2820.         private final T[][] currentRhoSigmaj;

  2822.         /** The C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> Fourier coefficients. */
  2823.         private final FieldFourierCjSjCoefficients<T> fourierCjSj;

  2825.         /** The maximum value for j index. */
  2826.         private final int jMax;

  2828.         /** Constructor.
  2829.          * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
  2830.          * @param fourierCjSj the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  2831.          * @param jMax maximum value for j index
  2832.          * @param field field used by default
  2833.          */
  2834.         FieldUijVijCoefficients(final T[][] currentRhoSigmaj,
  2835.                                 final FieldFourierCjSjCoefficients<T> fourierCjSj,
  2836.                                 final int jMax,
  2837.                                 final Field<T> field) {
  2838.             this.currentRhoSigmaj = currentRhoSigmaj;
  2839.             this.fourierCjSj = fourierCjSj;
  2840.             this.jMax = jMax;

  2842.             // initialize the internal arrays.
  2843.             this.u1ij = MathArrays.buildArray(field, 6, 2 * jMax + 1);
  2844.             this.v1ij = MathArrays.buildArray(field, 6, 2 * jMax + 1);
  2845.             this.u2ij = MathArrays.buildArray(field, jMax + 1);
  2846.             this.v2ij = MathArrays.buildArray(field, jMax + 1);

  2848.             //compute the coefficients
  2849.             computeU1V1Coefficients(field);
  2850.             computeU2V2Coefficients(field);
  2851.         }

  2853.         /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
  2854.          * @param field field used by default
  2855.          */
  2856.         private void computeU1V1Coefficients(final Field<T> field) {
  2857.             // Zero
  2858.             final T zero = field.getZero();

  2860.             // generate the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients
  2861.             // for j >= 1
  2862.             // also the U₁⁰ for Fourier index = 1 (i == 0) coefficient will be computed
  2863.             u1ij[0][0] = zero;
  2864.             for (int j = 1; j <= jMax; j++) {
  2865.                 // compute 1 / j
  2866.                 final double ooj = 1. / j;

  2868.                 for (int i = 0; i < 6; i++) {
  2869.                     //j is aready between 1 and J
  2870.                     u1ij[i][j] = fourierCjSj.getSij(i, j);
  2871.                     v1ij[i][j] = fourierCjSj.getCij(i, j);

  2873.                     // 1 - δ<sub>1j</sub> is 1 for all j > 1
  2874.                     if (j > 1) {
  2875.                         // k starts with 1 because j-J is less than or equal to 0
  2876.                         for (int kIndex = 1; kIndex <= j - 1; kIndex++) {
  2877.                             // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
  2878.                             // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
  2879.                             u1ij[i][j] = u1ij[i][j].add(fourierCjSj.getCij(i, j - kIndex).multiply(currentRhoSigmaj[1][kIndex]).
  2880.                                          add(fourierCjSj.getSij(i, j - kIndex).multiply(currentRhoSigmaj[0][kIndex])));

  2882.                             // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  2883.                             // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  2884.                             v1ij[i][j] = v1ij[i][j].add(fourierCjSj.getCij(i, j - kIndex).multiply(currentRhoSigmaj[0][kIndex]).
  2885.                                          subtract(fourierCjSj.getSij(i, j - kIndex).multiply(currentRhoSigmaj[1][kIndex])));
  2886.                         }
  2887.                     }

  2889.                     // since j must be between 1 and J-1 and is already between 1 and J
  2890.                     // the following sum is skiped only for j = jMax
  2891.                     if (j != jMax) {
  2892.                         for (int kIndex = 1; kIndex <= jMax - j; kIndex++) {
  2893.                             // -C<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub> +
  2894.                             // S<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub>
  2895.                             u1ij[i][j] = u1ij[i][j].add(fourierCjSj.getCij(i, j + kIndex).negate().multiply(currentRhoSigmaj[1][kIndex]).
  2896.                                          add(fourierCjSj.getSij(i, j + kIndex).multiply(currentRhoSigmaj[0][kIndex])));

  2898.                             // C<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub> +
  2899.                             // S<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub>
  2900.                             v1ij[i][j] = v1ij[i][j].add(fourierCjSj.getCij(i, j + kIndex).multiply(currentRhoSigmaj[0][kIndex]).
  2901.                                          add(fourierCjSj.getSij(i, j + kIndex).multiply(currentRhoSigmaj[1][kIndex])));
  2902.                         }
  2903.                     }

  2905.                     for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  2906.                         // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  2907.                         // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  2908.                         u1ij[i][j] = u1ij[i][j].add(fourierCjSj.getCij(i, kIndex).negate().multiply(currentRhoSigmaj[1][j + kIndex]).
  2909.                                      subtract(fourierCjSj.getSij(i, kIndex).multiply(currentRhoSigmaj[0][j + kIndex])));

  2911.                         // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  2912.                         // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  2913.                         v1ij[i][j] = v1ij[i][j].add(fourierCjSj.getCij(i, kIndex).multiply(currentRhoSigmaj[0][j + kIndex]).
  2914.                                      add(fourierCjSj.getSij(i, kIndex).multiply(currentRhoSigmaj[1][j + kIndex])));
  2915.                     }

  2917.                     // divide by 1 / j
  2918.                     u1ij[i][j] = u1ij[i][j].multiply(-ooj);
  2919.                     v1ij[i][j] = v1ij[i][j].multiply(ooj);

  2921.                     // if index = 1 (i == 0) add the computed terms to U₁⁰
  2922.                     if (i == 0) {
  2923.                         //- (U₁<sup>j</sup> * ρ<sub>j</sub> + V₁<sup>j</sup> * σ<sub>j</sub>
  2924.                         u1ij[0][0] = u1ij[0][0].add(u1ij[0][j].negate().multiply(currentRhoSigmaj[0][j]).subtract(v1ij[0][j].multiply(currentRhoSigmaj[1][j])));
  2925.                     }
  2926.                 }
  2927.             }

  2929.             // Terms with j > jMax are required only when computing the coefficients
  2930.             // U₂<sup>j</sup> and V₂<sup>j</sup>
  2931.             // and those coefficients are only required for Fourier index = 1 (i == 0).
  2932.             for (int j = jMax + 1; j <= 2 * jMax; j++) {
  2933.                 // compute 1 / j
  2934.                 final double ooj = 1. / j;
  2935.                 //the value of i is 0
  2936.                 u1ij[0][j] = zero;
  2937.                 v1ij[0][j] = zero;

  2939.                 //k starts from j-J as it is always greater than or equal to 1
  2940.                 for (int kIndex = j - jMax; kIndex <= j - 1; kIndex++) {
  2941.                     // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
  2942.                     // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
  2943.                     u1ij[0][j] = u1ij[0][j].add(fourierCjSj.getCij(0, j - kIndex).multiply(currentRhoSigmaj[1][kIndex]).
  2944.                                  add(fourierCjSj.getSij(0, j - kIndex).multiply(currentRhoSigmaj[0][kIndex])));

  2946.                     // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  2947.                     // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  2948.                     v1ij[0][j] = v1ij[0][j].add(fourierCjSj.getCij(0, j - kIndex).multiply(currentRhoSigmaj[0][kIndex]).
  2949.                                  subtract(fourierCjSj.getSij(0, j - kIndex).multiply(currentRhoSigmaj[1][kIndex])));
  2950.                 }
  2951.                 for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  2952.                     // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  2953.                     // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  2954.                     u1ij[0][j] = u1ij[0][j].add(fourierCjSj.getCij(0, kIndex).negate().multiply(currentRhoSigmaj[1][j + kIndex]).
  2955.                                  subtract(fourierCjSj.getSij(0, kIndex).multiply(currentRhoSigmaj[0][j + kIndex])));

  2957.                     // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  2958.                     // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  2959.                     v1ij[0][j] = v1ij[0][j].add(fourierCjSj.getCij(0, kIndex).multiply(currentRhoSigmaj[0][j + kIndex]).add(
  2960.                                  fourierCjSj.getSij(0, kIndex).multiply(currentRhoSigmaj[1][j + kIndex])));
  2961.                 }

  2963.                 // divide by 1 / j
  2964.                 u1ij[0][j] = u1ij[0][j].multiply(-ooj);
  2965.                 v1ij[0][j] = v1ij[0][j].multiply(ooj);
  2966.             }
  2967.         }

  2969.         /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
  2970.          * <p>
  2971.          * Only the coefficients for Fourier index = 1 (i == 0) are required.
  2972.          * </p>
  2973.          * @param field field used by default
  2974.          */
  2975.         private void computeU2V2Coefficients(final Field<T> field) {
  2976.             for (int j = 1; j <= jMax; j++) {
  2977.                 // compute 1 / j
  2978.                 final double ooj = 1. / j;

  2980.                 // only the values for i == 0 are computed
  2981.                 u2ij[j] = v1ij[0][j];
  2982.                 v2ij[j] = u1ij[0][j];

  2984.                 // 1 - δ<sub>1j</sub> is 1 for all j > 1
  2985.                 if (j > 1) {
  2986.                     for (int l = 1; l <= j - 1; l++) {
  2987.                         // U₁<sup>j-l</sup> * σ<sub>l</sub> +
  2988.                         // V₁<sup>j-l</sup> * ρ<sub>l</sub>
  2989.                         u2ij[j] = u2ij[j].add(u1ij[0][j - l].multiply(currentRhoSigmaj[1][l]).
  2990.                                   add(v1ij[0][j - l].multiply(currentRhoSigmaj[0][l])));

  2992.                         // U₁<sup>j-l</sup> * ρ<sub>l</sub> -
  2993.                         // V₁<sup>j-l</sup> * σ<sub>l</sub>
  2994.                         v2ij[j] = v2ij[j].add(u1ij[0][j - l].multiply(currentRhoSigmaj[0][l]).
  2995.                                   subtract(v1ij[0][j - l].multiply(currentRhoSigmaj[1][l])));
  2996.                     }
  2997.                 }

  2999.                 for (int l = 1; l <= jMax; l++) {
  3000.                     // -U₁<sup>j+l</sup> * σ<sub>l</sub> +
  3001.                     // U₁<sup>l</sup> * σ<sub>j+l</sub> +
  3002.                     // V₁<sup>j+l</sup> * ρ<sub>l</sub> -
  3003.                     // V₁<sup>l</sup> * ρ<sub>j+l</sub>
  3004.                     u2ij[j] = u2ij[j].add(u1ij[0][j + l].negate().multiply(currentRhoSigmaj[1][l]).
  3005.                               add(u1ij[0][l].multiply(currentRhoSigmaj[1][j + l])).
  3006.                               add(v1ij[0][j + l].multiply(currentRhoSigmaj[0][l])).
  3007.                               subtract(v1ij[0][l].multiply(currentRhoSigmaj[0][j + l])));

  3009.                     // U₁<sup>j+l</sup> * ρ<sub>l</sub> +
  3010.                     // U₁<sup>l</sup> * ρ<sub>j+l</sub> +
  3011.                     // V₁<sup>j+l</sup> * σ<sub>l</sub> +
  3012.                     // V₁<sup>l</sup> * σ<sub>j+l</sub>
  3013.                     u2ij[j] = u2ij[j].add(u1ij[0][j + l].multiply(currentRhoSigmaj[0][l]).
  3014.                               add(u1ij[0][l].multiply(currentRhoSigmaj[0][j + l])).
  3015.                               add(v1ij[0][j + l].multiply(currentRhoSigmaj[1][l])).
  3016.                               add(v1ij[0][l].multiply(currentRhoSigmaj[1][j + l])));
  3017.                 }

  3019.                 // divide by 1 / j
  3020.                 u2ij[j] = u2ij[j].multiply(-ooj);
  3021.                 v2ij[j] = v2ij[j].multiply(ooj);
  3022.             }
  3023.         }

  3025.         /** Get the coefficient U₁<sup>j</sup> for Fourier index i.
  3026.          *
  3027.          * @param j j index
  3028.          * @param i Fourier index (starts at 0)
  3029.          * @return the coefficient U₁<sup>j</sup> for the given Fourier index i
  3030.          */
  3031.         public T getU1(final int j, final int i) {
  3032.             return u1ij[i][j];
  3033.         }

  3035.         /** Get the coefficient V₁<sup>j</sup> for Fourier index i.
  3036.          *
  3037.          * @param j j index
  3038.          * @param i Fourier index (starts at 0)
  3039.          * @return the coefficient V₁<sup>j</sup> for the given Fourier index i
  3040.          */
  3041.         public T getV1(final int j, final int i) {
  3042.             return v1ij[i][j];
  3043.         }

  3045.         /** Get the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0).
  3046.          *
  3047.          * @param j j index
  3048.          * @return the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0)
  3049.          */
  3050.         public T getU2(final int j) {
  3051.             return u2ij[j];
  3052.         }

  3054.         /** Get the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0).
  3055.          *
  3056.          * @param j j index
  3057.          * @return the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0)
  3058.          */
  3059.         public T getV2(final int j) {
  3060.             return v2ij[j];
  3061.         }
  3062.     }

  3064.     /** Coefficients valid for one time slot. */
  3065.     private static class Slot {

  3067.         /**The coefficients D<sub>i</sub><sup>j</sup>.
  3068.          * <p>
  3069.          * Only for j = 1 and j = 2 the coefficients are not 0. <br>
  3070.          * i corresponds to the equinoctial element, as follows:
  3071.          * - i=0 for a <br/>
  3072.          * - i=1 for k <br/>
  3073.          * - i=2 for h <br/>
  3074.          * - i=3 for q <br/>
  3075.          * - i=4 for p <br/>
  3076.          * - i=5 for λ <br/>
  3077.          * </p>
  3078.          */
  3079.         private final ShortPeriodicsInterpolatedCoefficient[] dij;

  3081.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  3082.          * <p>
  3083.          * The index order is cij[j][i] <br/>
  3084.          * i corresponds to the equinoctial element, as follows: <br/>
  3085.          * - i=0 for a <br/>
  3086.          * - i=1 for k <br/>
  3087.          * - i=2 for h <br/>
  3088.          * - i=3 for q <br/>
  3089.          * - i=4 for p <br/>
  3090.          * - i=5 for λ <br/>
  3091.          * </p>
  3092.          */
  3093.         private final ShortPeriodicsInterpolatedCoefficient[] cij;

  3095.         /** The coefficients S<sub>i</sub><sup>j</sup>.
  3096.          * <p>
  3097.          * The index order is sij[j][i] <br/>
  3098.          * i corresponds to the equinoctial element, as follows: <br/>
  3099.          * - i=0 for a <br/>
  3100.          * - i=1 for k <br/>
  3101.          * - i=2 for h <br/>
  3102.          * - i=3 for q <br/>
  3103.          * - i=4 for p <br/>
  3104.          * - i=5 for λ <br/>
  3105.          * </p>
  3106.          */
  3107.         private final ShortPeriodicsInterpolatedCoefficient[] sij;

  3109.         /** Simple constructor.
  3110.          *  @param jMax maximum value for j index
  3111.          *  @param interpolationPoints number of points used in the interpolation process
  3112.          */
  3113.         Slot(final int jMax, final int interpolationPoints) {

  3115.             dij = new ShortPeriodicsInterpolatedCoefficient[3];
  3116.             cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
  3117.             sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];

  3119.             // 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
  3120.             for (int j = 0; j <= jMax; j++) {
  3121.                 cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3122.                 if (j > 0) {
  3123.                     sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3124.                 }
  3125.                 // Initialize only the non-zero D<sub>i</sub><sup>j</sup> coefficients
  3126.                 if (j == 1 || j == 2) {
  3127.                     dij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3128.                 }
  3129.             }

  3131.         }

  3133.     }

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

  3138.         /**The coefficients D<sub>i</sub><sup>j</sup>.
  3139.          * <p>
  3140.          * Only for j = 1 and j = 2 the coefficients are not 0. <br>
  3141.          * i corresponds to the equinoctial element, as follows:
  3142.          * - i=0 for a <br/>
  3143.          * - i=1 for k <br/>
  3144.          * - i=2 for h <br/>
  3145.          * - i=3 for q <br/>
  3146.          * - i=4 for p <br/>
  3147.          * - i=5 for λ <br/>
  3148.          * </p>
  3149.          */
  3150.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] dij;

  3152.         /** The coefficients C<sub>i</sub><sup>j</sup>.
  3153.          * <p>
  3154.          * The index order is cij[j][i] <br/>
  3155.          * i corresponds to the equinoctial element, as follows: <br/>
  3156.          * - i=0 for a <br/>
  3157.          * - i=1 for k <br/>
  3158.          * - i=2 for h <br/>
  3159.          * - i=3 for q <br/>
  3160.          * - i=4 for p <br/>
  3161.          * - i=5 for λ <br/>
  3162.          * </p>
  3163.          */
  3164.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] cij;

  3166.         /** The coefficients S<sub>i</sub><sup>j</sup>.
  3167.          * <p>
  3168.          * The index order is sij[j][i] <br/>
  3169.          * i corresponds to the equinoctial element, as follows: <br/>
  3170.          * - i=0 for a <br/>
  3171.          * - i=1 for k <br/>
  3172.          * - i=2 for h <br/>
  3173.          * - i=3 for q <br/>
  3174.          * - i=4 for p <br/>
  3175.          * - i=5 for λ <br/>
  3176.          * </p>
  3177.          */
  3178.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] sij;

  3180.         /** Simple constructor.
  3181.          *  @param jMax maximum value for j index
  3182.          *  @param interpolationPoints number of points used in the interpolation process
  3183.          */
  3184.         @SuppressWarnings("unchecked")
  3185.         FieldSlot(final int jMax, final int interpolationPoints) {

  3187.             dij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, 3);
  3188.             cij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);
  3189.             sij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);

  3191.             // 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
  3192.             for (int j = 0; j <= jMax; j++) {
  3193.                 cij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3194.                 if (j > 0) {
  3195.                     sij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3196.                 }
  3197.                 // Initialize only the non-zero D<sub>i</sub><sup>j</sup> coefficients
  3198.                 if (j == 1 || j == 2) {
  3199.                     dij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3200.                 }
  3201.             }

  3203.         }

  3205.     }

  3207. }
Created with JaCoCo 0.8.5.201910111838