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.FieldSinCos;
  35. import org.hipparchus.util.MathArrays;
  36. import org.hipparchus.util.SinCos;
  37. import org.orekit.attitudes.Attitude;
  38. import org.orekit.attitudes.AttitudeProvider;
  39. import org.orekit.attitudes.FieldAttitude;
  40. import org.orekit.forces.ForceModel;
  41. import org.orekit.orbits.EquinoctialOrbit;
  42. import org.orekit.orbits.FieldEquinoctialOrbit;
  43. import org.orekit.orbits.FieldOrbit;
  44. import org.orekit.orbits.Orbit;
  45. import org.orekit.orbits.OrbitType;
  46. import org.orekit.orbits.PositionAngle;
  47. import org.orekit.propagation.FieldSpacecraftState;
  48. import org.orekit.propagation.PropagationType;
  49. import org.orekit.propagation.SpacecraftState;
  50. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  51. import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
  52. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  53. import org.orekit.propagation.semianalytical.dsst.utilities.FieldCjSjCoefficient;
  54. import org.orekit.propagation.semianalytical.dsst.utilities.FieldShortPeriodicsInterpolatedCoefficient;
  55. import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
  56. import org.orekit.time.AbsoluteDate;
  57. import org.orekit.time.FieldAbsoluteDate;
  58. import org.orekit.utils.FieldTimeSpanMap;
  59. import org.orekit.utils.ParameterDriver;
  60. import org.orekit.utils.TimeSpanMap;

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

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

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

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

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

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

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

  127.     /** Contribution to be numerically averaged. */
  128.     private final ForceModel contribution;

  130.     /** Gauss integrator. */
  131.     private final double threshold;

  133.     /** Gauss integrator. */
  134.     private GaussQuadrature integrator;

  136.     /** Flag for Gauss order computation. */
  137.     private boolean isDirty;

  139.     /** Attitude provider. */
  140.     private AttitudeProvider attitudeProvider;

  142.     /** Prefix for coefficients keys. */
  143.     private final String coefficientsKeyPrefix;

  145.     /** Short period terms. */
  146.     private GaussianShortPeriodicCoefficients gaussianSPCoefs;

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

  151.     /** Driver for gravitational parameter. */
  152.     private final ParameterDriver gmParameterDriver;

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

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

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

  175.         gaussianFieldSPCoefs       = new HashMap<>();
  176.     }

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

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

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

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

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

  197.     }

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

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

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

  222.     }

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

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

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

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

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

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

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

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

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

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

  325.         return meanElementRate;
  326.     }

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

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

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

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

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

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

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

  395.         // Field
  396.         final Field<T> field = context.getA().getField();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  482.         }

  484.     }

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

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

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

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

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

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

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

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

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

  516.         }

  518.     }

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

  536.         // (-b)<sup>j</sup>
  537.         double mbtj = 1;

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

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

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

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

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

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

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

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

  595.         /** Current state. */
  596.         private final FieldSpacecraftState<T> state;

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

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

  608.         /** Container for attributes. */
  609.         private final FieldAbstractGaussianContributionContext<T> context;

  611.         /** Auxiliary Elements. */
  612.         private final FieldAuxiliaryElements<T> auxiliaryElements;

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

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

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

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

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

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

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

  672.             // Compute acceleration
  673.             FieldVector3D<T> acc = FieldVector3D.getZero(field);

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

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

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

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

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

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

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

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

  741.             return val;
  742.         }

  744.         /** Converts true longitude to mean longitude.
  745.          * @param x True longitude
  746.          * @return Eccentric longitude
  747.          */
  748.         private T trueToMean (final T x) {
  749.             return eccentricToMean(trueToEccentric(x));
  750.         }

  752.         /** Converts true longitude to eccentric longitude.
  753.          * @param lv True longitude
  754.          * @return Eccentric longitude
  755.          */
  756.         private T trueToEccentric (final T lv) {
  757.             final FieldSinCos<T> sclV = FastMath.sinCos(lv);
  758.             final T cosLv   = sclV.cos();
  759.             final T sinLv   = sclV.sin();
  760.             final T num     = auxiliaryElements.getH().multiply(cosLv).subtract(auxiliaryElements.getK().multiply(sinLv));
  761.             final T den     = auxiliaryElements.getB().add(auxiliaryElements.getK().multiply(cosLv)).add(auxiliaryElements.getH().multiply(sinLv)).add(1.);
  762.             return FastMath.atan(num.divide(den)).multiply(2.).add(lv);
  763.         }

  765.         /** Converts eccentric longitude to mean longitude.
  766.          * @param le Eccentric longitude
  767.          * @return Mean longitude
  768.          */
  769.         private T eccentricToMean (final T le) {
  770.             final FieldSinCos<T> scle = FastMath.sinCos(le);
  771.             return le.subtract(auxiliaryElements.getK().multiply(scle.sin())).add(auxiliaryElements.getH().multiply(scle.cos()));
  772.         }

  774.         /** Compute δa/δv.
  775.          *  @param vel satellite velocity
  776.          *  @return δa/δv
  777.          */
  778.         private FieldVector3D<T> getAoV(final FieldVector3D<T> vel) {
  779.             return new FieldVector3D<>(context.getTon2a(), vel);
  780.         }

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

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

  810.         /** Compute δp/δv.
  811.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  812.          *  @return δp/δv
  813.          */
  814.         private FieldVector3D<T> getPoV(final T Y) {
  815.             return new FieldVector3D<>(context.getCo2AB().multiply(Y), auxiliaryElements.getVectorW());
  816.         }

  818.         /** Compute δq/δv.
  819.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  820.          *  @return δq/δv
  821.          */
  822.         private FieldVector3D<T> getQoV(final T X) {
  823.             return new FieldVector3D<>(context.getCo2AB().multiply(X).multiply(I), auxiliaryElements.getVectorW());
  824.         }

  826.         /** Compute δλ/δv.
  827.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  828.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  829.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  830.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  831.          *  @return δλ/δv
  832.          */
  833.         private FieldVector3D<T> getLoV(final T X, final T Y, final T Xdot, final T Ydot) {
  834.             final FieldVector3D<T> pos = new FieldVector3D<>(X, auxiliaryElements.getVectorF(), Y, auxiliaryElements.getVectorG());
  835.             final FieldVector3D<T> v2  = new FieldVector3D<>(auxiliaryElements.getK(), getHoV(X, Y, Xdot, Ydot), auxiliaryElements.getH().negate(), getKoV(X, Y, Xdot, Ydot));
  836.             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());
  837.         }

  839.     }

  841.     /** Internal class for numerical quadrature. */
  842.     private class IntegrableFunction implements UnivariateVectorFunction {

  844.         /** Current state. */
  845.         private final SpacecraftState state;

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

  851.         /** The j index.
  852.          * <p>
  853.          * Used only for short periodic variation. Ignored for mean elements variation.
  854.          * </p> */
  855.         private final int j;

  857.         /** Container for attributes. */
  858.         private final AbstractGaussianContributionContext context;

  860.         /** Auxiliary Elements. */
  861.         private final AuxiliaryElements auxiliaryElements;

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

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

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

  893.         /** {@inheritDoc} */
  894.         @Override
  895.         public double[] value(final double x) {

  897.             //Compute the time difference from the true longitude difference
  898.             final double shiftedLm = trueToMean(x);
  899.             final double dLm = shiftedLm - auxiliaryElements.getLM();
  900.             final double dt = dLm / context.getMeanMotion();

  902.             final SinCos scL  = FastMath.sinCos(x);
  903.             final double cosL = scL.cos();
  904.             final double sinL = scL.sin();
  905.             final double roa  = auxiliaryElements.getB() * auxiliaryElements.getB() / (1. + auxiliaryElements.getH() * sinL + auxiliaryElements.getK() * cosL);
  906.             final double roa2 = roa * roa;
  907.             final double r    = auxiliaryElements.getSma() * roa;
  908.             final double X    = r * cosL;
  909.             final double Y    = r * sinL;
  910.             final double naob = context.getMeanMotion() * auxiliaryElements.getSma() / auxiliaryElements.getB();
  911.             final double Xdot = -naob * (auxiliaryElements.getH() + sinL);
  912.             final double Ydot =  naob * (auxiliaryElements.getK() + cosL);
  913.             final Vector3D vel = new Vector3D(Xdot, auxiliaryElements.getVectorF(), Ydot, auxiliaryElements.getVectorG());

  915.             // Compute acceleration
  916.             Vector3D acc = Vector3D.ZERO;

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

  921.             // Recompose an orbit with time held fixed to be compliant with DSST theory
  922.             final Orbit recomposedOrbit =
  923.                             new EquinoctialOrbit(shiftedOrbit.getA(),
  924.                                                  shiftedOrbit.getEquinoctialEx(),
  925.                                                  shiftedOrbit.getEquinoctialEy(),
  926.                                                  shiftedOrbit.getHx(),
  927.                                                  shiftedOrbit.getHy(),
  928.                                                  shiftedOrbit.getLv(),
  929.                                                  PositionAngle.TRUE,
  930.                                                  shiftedOrbit.getFrame(),
  931.                                                  state.getDate(),
  932.                                                  context.getMu());

  934.             // Get the corresponding attitude
  935.             final Attitude recomposedAttitude =
  936.                     attitudeProvider.getAttitude(recomposedOrbit,
  937.                                                  recomposedOrbit.getDate(),
  938.                                                  recomposedOrbit.getFrame());

  940.             // create shifted SpacecraftState with attitude at specified time
  941.             final SpacecraftState shiftedState =
  942.                     new SpacecraftState(recomposedOrbit, recomposedAttitude, state.getMass());

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

  946.             //Compute the derivatives of the elements by the speed
  947.             final double[] deriv = new double[6];
  948.             // da/dv
  949.             deriv[0] = getAoV(vel).dotProduct(acc);
  950.             // dex/dv
  951.             deriv[1] = getKoV(X, Y, Xdot, Ydot).dotProduct(acc);
  952.             // dey/dv
  953.             deriv[2] = getHoV(X, Y, Xdot, Ydot).dotProduct(acc);
  954.             // dhx/dv
  955.             deriv[3] = getQoV(X).dotProduct(acc);
  956.             // dhy/dv
  957.             deriv[4] = getPoV(Y).dotProduct(acc);
  958.             // dλ/dv
  959.             deriv[5] = getLoV(X, Y, Xdot, Ydot).dotProduct(acc);

  961.             // Compute mean elements rates
  962.             double[] val = null;
  963.             if (meanMode) {
  964.                 val = new double[6];
  965.                 for (int i = 0; i < 6; i++) {
  966.                     // da<sub>i</sub>/dt
  967.                     val[i] = roa2 * deriv[i];
  968.                 }
  969.             } else {
  970.                 val = new double[12];
  971.                 //Compute cos(j*L) and sin(j*L);
  972.                 final SinCos scjL  = FastMath.sinCos(j * x);
  973.                 final double cosjL = j == 1 ? cosL : scjL.cos();
  974.                 final double sinjL = j == 1 ? sinL : scjL.sin();

  976.                 for (int i = 0; i < 6; i++) {
  977.                     // da<sub>i</sub>/dv * cos(jL)
  978.                     val[i] = cosjL * deriv[i];
  979.                     // da<sub>i</sub>/dv * sin(jL)
  980.                     val[i + 6] = sinjL * deriv[i];
  981.                 }
  982.             }
  983.             return val;
  984.         }

  986.         /** Converts true longitude to eccentric longitude.
  987.          * @param lv True longitude
  988.          * @return Eccentric longitude
  989.          */
  990.         private double trueToEccentric (final double lv) {
  991.             final SinCos scLv    = FastMath.sinCos(lv);
  992.             final double num     = auxiliaryElements.getH() * scLv.cos() - auxiliaryElements.getK() * scLv.sin();
  993.             final double den     = auxiliaryElements.getB() + 1. + auxiliaryElements.getK() * scLv.cos() + auxiliaryElements.getH() * scLv.sin();
  994.             return lv + 2. * FastMath.atan(num / den);
  995.         }

  997.         /** Converts eccentric longitude to mean longitude.
  998.          * @param le Eccentric longitude
  999.          * @return Mean longitude
  1000.          */
  1001.         private double eccentricToMean (final double le) {
  1002.             final SinCos scLe = FastMath.sinCos(le);
  1003.             return le - auxiliaryElements.getK() * scLe.sin() + auxiliaryElements.getH() * scLe.cos();
  1004.         }

  1006.         /** Converts true longitude to mean longitude.
  1007.          * @param lv True longitude
  1008.          * @return Eccentric longitude
  1009.          */
  1010.         private double trueToMean (final double lv) {
  1011.             return eccentricToMean(trueToEccentric(lv));
  1012.         }

  1014.         /** Compute δa/δv.
  1015.          *  @param vel satellite velocity
  1016.          *  @return δa/δv
  1017.          */
  1018.         private Vector3D getAoV(final Vector3D vel) {
  1019.             return new Vector3D(context.getTon2a(), vel);
  1020.         }

  1022.         /** Compute δh/δv.
  1023.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  1024.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  1025.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  1026.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  1027.          *  @return δh/δv
  1028.          */
  1029.         private Vector3D getHoV(final double X, final double Y, final double Xdot, final double Ydot) {
  1030.             final double kf = (2. * Xdot * Y - X * Ydot) * context.getOoMU();
  1031.             final double kg = X * Xdot * context.getOoMU();
  1032.             final double kw = auxiliaryElements.getK() * (I * auxiliaryElements.getQ() * Y - auxiliaryElements.getP() * X) * context.getOOAB();
  1033.             return new Vector3D(kf, auxiliaryElements.getVectorF(), -kg, auxiliaryElements.getVectorG(), kw, auxiliaryElements.getVectorW());
  1034.         }

  1036.         /** Compute δk/δv.
  1037.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  1038.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  1039.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  1040.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  1041.          *  @return δk/δv
  1042.          */
  1043.         private Vector3D getKoV(final double X, final double Y, final double Xdot, final double Ydot) {
  1044.             final double kf = Y * Ydot * context.getOoMU();
  1045.             final double kg = (2. * X * Ydot - Xdot * Y) * context.getOoMU();
  1046.             final double kw = auxiliaryElements.getH() * (I * auxiliaryElements.getQ() * Y - auxiliaryElements.getP() * X) * context.getOOAB();
  1047.             return new Vector3D(-kf, auxiliaryElements.getVectorF(), kg, auxiliaryElements.getVectorG(), -kw, auxiliaryElements.getVectorW());
  1048.         }

  1050.         /** Compute δp/δv.
  1051.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  1052.          *  @return δp/δv
  1053.          */
  1054.         private Vector3D getPoV(final double Y) {
  1055.             return new Vector3D(context.getCo2AB() * Y, auxiliaryElements.getVectorW());
  1056.         }

  1058.         /** Compute δq/δv.
  1059.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  1060.          *  @return δq/δv
  1061.          */
  1062.         private Vector3D getQoV(final double X) {
  1063.             return new Vector3D(I * context.getCo2AB() * X, auxiliaryElements.getVectorW());
  1064.         }

  1066.         /** Compute δλ/δv.
  1067.          *  @param X satellite position component along f, equinoctial reference frame 1st vector
  1068.          *  @param Y satellite position component along g, equinoctial reference frame 2nd vector
  1069.          *  @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
  1070.          *  @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
  1071.          *  @return δλ/δv
  1072.          */
  1073.         private Vector3D getLoV(final double X, final double Y, final double Xdot, final double Ydot) {
  1074.             final Vector3D pos = new Vector3D(X, auxiliaryElements.getVectorF(), Y, auxiliaryElements.getVectorG());
  1075.             final Vector3D v2  = new Vector3D(auxiliaryElements.getK(), getHoV(X, Y, Xdot, Ydot), -auxiliaryElements.getH(), getKoV(X, Y, Xdot, Ydot));
  1076.             return new Vector3D(-2. * context.getOOA(), pos, context.getOoBpo(), v2, (I * auxiliaryElements.getQ() * Y - auxiliaryElements.getP() * X) * context.getOOA(), auxiliaryElements.getVectorW());
  1077.         }

  1079.     }

  1081.     /** Class used to {@link #integrate(UnivariateVectorFunction, double, double) integrate}
  1082.      *  a {@link org.hipparchus.analysis.UnivariateVectorFunction function}
  1083.      *  of the orbital elements using the Gaussian quadrature rule to get the acceleration.
  1084.      */
  1085.     private static class GaussQuadrature {

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

  1089.         /** Points for quadrature of order 12. */
  1090.         private static final double[] P_12 = {
  1091.             -0.98156063424671910000,
  1092.             -0.90411725637047490000,
  1093.             -0.76990267419430470000,
  1094.             -0.58731795428661740000,
  1095.             -0.36783149899818024000,
  1096.             -0.12523340851146890000,
  1097.             0.12523340851146890000,
  1098.             0.36783149899818024000,
  1099.             0.58731795428661740000,
  1100.             0.76990267419430470000,
  1101.             0.90411725637047490000,
  1102.             0.98156063424671910000
  1103.         };

  1105.         /** Weights for quadrature of order 12. */
  1106.         private static final double[] W_12 = {
  1107.             0.04717533638651220000,
  1108.             0.10693932599531830000,
  1109.             0.16007832854334633000,
  1110.             0.20316742672306584000,
  1111.             0.23349253653835478000,
  1112.             0.24914704581340286000,
  1113.             0.24914704581340286000,
  1114.             0.23349253653835478000,
  1115.             0.20316742672306584000,
  1116.             0.16007832854334633000,
  1117.             0.10693932599531830000,
  1118.             0.04717533638651220000
  1119.         };

  1121.         /** Points for quadrature of order 16. */
  1122.         private static final double[] P_16 = {
  1123.             -0.98940093499164990000,
  1124.             -0.94457502307323260000,
  1125.             -0.86563120238783160000,
  1126.             -0.75540440835500310000,
  1127.             -0.61787624440264380000,
  1128.             -0.45801677765722737000,
  1129.             -0.28160355077925890000,
  1130.             -0.09501250983763745000,
  1131.             0.09501250983763745000,
  1132.             0.28160355077925890000,
  1133.             0.45801677765722737000,
  1134.             0.61787624440264380000,
  1135.             0.75540440835500310000,
  1136.             0.86563120238783160000,
  1137.             0.94457502307323260000,
  1138.             0.98940093499164990000
  1139.         };

  1141.         /** Weights for quadrature of order 16. */
  1142.         private static final double[] W_16 = {
  1143.             0.02715245941175405800,
  1144.             0.06225352393864777000,
  1145.             0.09515851168249283000,
  1146.             0.12462897125553388000,
  1147.             0.14959598881657685000,
  1148.             0.16915651939500256000,
  1149.             0.18260341504492360000,
  1150.             0.18945061045506847000,
  1151.             0.18945061045506847000,
  1152.             0.18260341504492360000,
  1153.             0.16915651939500256000,
  1154.             0.14959598881657685000,
  1155.             0.12462897125553388000,
  1156.             0.09515851168249283000,
  1157.             0.06225352393864777000,
  1158.             0.02715245941175405800
  1159.         };

  1161.         /** Points for quadrature of order 20. */
  1162.         private static final double[] P_20 = {
  1163.             -0.99312859918509490000,
  1164.             -0.96397192727791390000,
  1165.             -0.91223442825132600000,
  1166.             -0.83911697182221890000,
  1167.             -0.74633190646015080000,
  1168.             -0.63605368072651510000,
  1169.             -0.51086700195082700000,
  1170.             -0.37370608871541955000,
  1171.             -0.22778585114164507000,
  1172.             -0.07652652113349734000,
  1173.             0.07652652113349734000,
  1174.             0.22778585114164507000,
  1175.             0.37370608871541955000,
  1176.             0.51086700195082700000,
  1177.             0.63605368072651510000,
  1178.             0.74633190646015080000,
  1179.             0.83911697182221890000,
  1180.             0.91223442825132600000,
  1181.             0.96397192727791390000,
  1182.             0.99312859918509490000
  1183.         };

  1185.         /** Weights for quadrature of order 20. */
  1186.         private static final double[] W_20 = {
  1187.             0.01761400713915226400,
  1188.             0.04060142980038684000,
  1189.             0.06267204833410904000,
  1190.             0.08327674157670477000,
  1191.             0.10193011981724048000,
  1192.             0.11819453196151844000,
  1193.             0.13168863844917678000,
  1194.             0.14209610931838212000,
  1195.             0.14917298647260380000,
  1196.             0.15275338713072600000,
  1197.             0.15275338713072600000,
  1198.             0.14917298647260380000,
  1199.             0.14209610931838212000,
  1200.             0.13168863844917678000,
  1201.             0.11819453196151844000,
  1202.             0.10193011981724048000,
  1203.             0.08327674157670477000,
  1204.             0.06267204833410904000,
  1205.             0.04060142980038684000,
  1206.             0.01761400713915226400
  1207.         };

  1209.         /** Points for quadrature of order 24. */
  1210.         private static final double[] P_24 = {
  1211.             -0.99518721999702130000,
  1212.             -0.97472855597130950000,
  1213.             -0.93827455200273270000,
  1214.             -0.88641552700440100000,
  1215.             -0.82000198597390300000,
  1216.             -0.74012419157855440000,
  1217.             -0.64809365193697550000,
  1218.             -0.54542147138883950000,
  1219.             -0.43379350762604520000,
  1220.             -0.31504267969616340000,
  1221.             -0.19111886747361634000,
  1222.             -0.06405689286260563000,
  1223.             0.06405689286260563000,
  1224.             0.19111886747361634000,
  1225.             0.31504267969616340000,
  1226.             0.43379350762604520000,
  1227.             0.54542147138883950000,
  1228.             0.64809365193697550000,
  1229.             0.74012419157855440000,
  1230.             0.82000198597390300000,
  1231.             0.88641552700440100000,
  1232.             0.93827455200273270000,
  1233.             0.97472855597130950000,
  1234.             0.99518721999702130000
  1235.         };

  1237.         /** Weights for quadrature of order 24. */
  1238.         private static final double[] W_24 = {
  1239.             0.01234122979998733500,
  1240.             0.02853138862893380600,
  1241.             0.04427743881741981000,
  1242.             0.05929858491543691500,
  1243.             0.07334648141108027000,
  1244.             0.08619016153195320000,
  1245.             0.09761865210411391000,
  1246.             0.10744427011596558000,
  1247.             0.11550566805372553000,
  1248.             0.12167047292780335000,
  1249.             0.12583745634682825000,
  1250.             0.12793819534675221000,
  1251.             0.12793819534675221000,
  1252.             0.12583745634682825000,
  1253.             0.12167047292780335000,
  1254.             0.11550566805372553000,
  1255.             0.10744427011596558000,
  1256.             0.09761865210411391000,
  1257.             0.08619016153195320000,
  1258.             0.07334648141108027000,
  1259.             0.05929858491543691500,
  1260.             0.04427743881741981000,
  1261.             0.02853138862893380600,
  1262.             0.01234122979998733500
  1263.         };

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

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

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

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

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

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

  1529.         /** Node points. */
  1530.         private final double[] nodePoints;

  1532.         /** Node weights. */
  1533.         private final double[] nodeWeights;

  1535.         /** Number of points. */
  1536.         private final int numberOfPoints;

  1538.         /** Creates a Gauss integrator of the given order.
  1539.          *
  1540.          *  @param numberOfPoints Order of the integration rule.
  1541.          */
  1542.         GaussQuadrature(final int numberOfPoints) {

  1544.             this.numberOfPoints = numberOfPoints;

  1546.             switch(numberOfPoints) {
  1547.                 case 12 :
  1548.                     this.nodePoints  = P_12.clone();
  1549.                     this.nodeWeights = W_12.clone();
  1550.                     break;
  1551.                 case 16 :
  1552.                     this.nodePoints  = P_16.clone();
  1553.                     this.nodeWeights = W_16.clone();
  1554.                     break;
  1555.                 case 20 :
  1556.                     this.nodePoints  = P_20.clone();
  1557.                     this.nodeWeights = W_20.clone();
  1558.                     break;
  1559.                 case 24 :
  1560.                     this.nodePoints  = P_24.clone();
  1561.                     this.nodeWeights = W_24.clone();
  1562.                     break;
  1563.                 case 32 :
  1564.                     this.nodePoints  = P_32.clone();
  1565.                     this.nodeWeights = W_32.clone();
  1566.                     break;
  1567.                 case 40 :
  1568.                     this.nodePoints  = P_40.clone();
  1569.                     this.nodeWeights = W_40.clone();
  1570.                     break;
  1571.                 case 48 :
  1572.                 default :
  1573.                     this.nodePoints  = P_48.clone();
  1574.                     this.nodeWeights = W_48.clone();
  1575.                     break;
  1576.             }

  1578.         }

  1580.         /** Integrates a given function on the given interval.
  1581.         *
  1582.         *  @param f Function to integrate.
  1583.         *  @param lowerBound Lower bound of the integration interval.
  1584.         *  @param upperBound Upper bound of the integration interval.
  1585.         *  @return the integral of the weighted function.
  1586.         */
  1587.         public double[] integrate(final UnivariateVectorFunction f,
  1588.                final double lowerBound, final double upperBound) {

  1590.             final double[] adaptedPoints  = nodePoints.clone();
  1591.             final double[] adaptedWeights = nodeWeights.clone();
  1592.             transform(adaptedPoints, adaptedWeights, lowerBound, upperBound);
  1593.             return basicIntegrate(f, adaptedPoints, adaptedWeights);
  1594.         }

  1596.        /** Integrates a given function on the given interval.
  1597.        *
  1598.        *  @param <T> the type of the field elements
  1599.        *  @param f Function to integrate.
  1600.        *  @param lowerBound Lower bound of the integration interval.
  1601.        *  @param upperBound Upper bound of the integration interval.
  1602.        *  @param field field utilized by default
  1603.        *  @return the integral of the weighted function.
  1604.        */
  1605.         public <T extends RealFieldElement<T>> T[] integrate(final RealFieldUnivariateVectorFunction<T> f,
  1606.                                                              final T lowerBound, final T upperBound,
  1607.                                                              final Field<T> field) {

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

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

  1614.             for (int i = 0; i < numberOfPoints; i++) {
  1615.                 adaptedPoints[i]  = zero.add(nodePoints[i]);
  1616.                 adaptedWeights[i] = zero.add(nodeWeights[i]);
  1617.             }

  1619.             transform(adaptedPoints, adaptedWeights, lowerBound, upperBound);
  1620.             return basicIntegrate(f, adaptedPoints, adaptedWeights, field);
  1621.         }

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

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

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

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

  1719.             T x = points[0];
  1720.             T w = weights[0];
  1721.             T[] v = f.value(x);

  1723.             final T[] y = MathArrays.buildArray(field, v.length);
  1724.             for (int j = 0; j < v.length; j++) {
  1725.                 y[j] = v[j].multiply(w);
  1726.             }
  1727.             final T[] t = y.clone();
  1728.             final T[] c = MathArrays.buildArray(field, v.length);;
  1729.             final T[] s = t.clone();
  1730.             for (int i = 1; i < points.length; i++) {
  1731.                 x = points[i];
  1732.                 w = weights[i];
  1733.                 v = f.value(x);
  1734.                 for (int j = 0; j < v.length; j++) {
  1735.                     y[j] = v[j].multiply(w).subtract(c[j]);
  1736.                     t[j] = y[j].add(s[j]);
  1737.                     c[j] = (t[j].subtract(s[j])).subtract(y[j]);
  1738.                     s[j] = t[j];
  1739.                 }
  1740.             }
  1741.             return s;
  1742.         }

  1744.     }

  1746.     /** Compute the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
  1747.      *  <p>
  1748.      *  Those coefficients are given in Danielson paper by expression 4.4-(6)
  1749.      *  </p>
  1750.      *  @author Petre Bazavan
  1751.      *  @author Lucian Barbulescu
  1752.      */
  1753.     private class FourierCjSjCoefficients {

  1755.         /** Maximum possible value for j. */
  1756.         private final int jMax;

  1758.         /** The C<sub>i</sub><sup>j</sup> coefficients.
  1759.          * <p>
  1760.          * the index i corresponds to the following elements: <br/>
  1761.          * - 0 for a <br>
  1762.          * - 1 for k <br>
  1763.          * - 2 for h <br>
  1764.          * - 3 for q <br>
  1765.          * - 4 for p <br>
  1766.          * - 5 for λ <br>
  1767.          * </p>
  1768.          */
  1769.         private final double[][] cCoef;

  1771.         /** The C<sub>i</sub><sup>j</sup> coefficients.
  1772.          * <p>
  1773.          * the index i corresponds to the following elements: <br/>
  1774.          * - 0 for a <br>
  1775.          * - 1 for k <br>
  1776.          * - 2 for h <br>
  1777.          * - 3 for q <br>
  1778.          * - 4 for p <br>
  1779.          * - 5 for λ <br>
  1780.          * </p>
  1781.          */
  1782.         private final double[][] sCoef;

  1784.         /** Standard constructor.
  1785.          * @param state the current state
  1786.          * @param jMax maximum value for j
  1787.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1788.          * @param parameters values of the force model parameters
  1789.          */
  1790.         FourierCjSjCoefficients(final SpacecraftState state, final int jMax,
  1791.                                 final AuxiliaryElements auxiliaryElements, final double[] parameters) {

  1793.             //Initialise the fields
  1794.             this.jMax = jMax;

  1796.             //Allocate the arrays
  1797.             final int rows = jMax + 1;
  1798.             cCoef = new double[rows][6];
  1799.             sCoef = new double[rows][6];

  1801.             //Compute the coefficients
  1802.             computeCoefficients(state, auxiliaryElements, parameters);
  1803.         }

  1805.         /**
  1806.          * Compute the Fourrier coefficients.
  1807.          * <p>
  1808.          * Only the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients need to be computed
  1809.          * as D<sub>i</sub><sup>m</sup> is always 0.
  1810.          * </p>
  1811.          * @param state the current state
  1812.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1813.          * @param parameters values of the force model parameters
  1814.          */
  1815.         private void computeCoefficients(final SpacecraftState state,
  1816.                                          final AuxiliaryElements auxiliaryElements,
  1817.                                          final double[] parameters) {

  1819.             // Computes the limits for the integral
  1820.             final double[] ll = getLLimits(state, auxiliaryElements);
  1821.             // Computes integrated mean element rates if Llow < Lhigh
  1822.             if (ll[0] < ll[1]) {
  1823.                 //Compute 1 / PI
  1824.                 final double ooPI = 1 / FastMath.PI;

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

  1831.                     //divide by PI and set the values for the coefficients
  1832.                     for (int i = 0; i < 6; i++) {
  1833.                         cCoef[j][i] = ooPI * curentCoefficients[i];
  1834.                         sCoef[j][i] = ooPI * curentCoefficients[i + 6];
  1835.                     }
  1836.                 }
  1837.             }
  1838.         }

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

  1849.         /** Get the coefficient S<sub>i</sub><sup>j</sup>.
  1850.          * @param i i index - corresponds to the required variation
  1851.          * @param j j index
  1852.          * @return the coefficient S<sub>i</sub><sup>j</sup>
  1853.          */
  1854.         public double getSij(final int i, final int j) {
  1855.             return sCoef[j][i];
  1856.         }
  1857.     }

  1859.     /** Compute the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients with field elements.
  1860.      *  <p>
  1861.      *  Those coefficients are given in Danielson paper by expression 4.4-(6)
  1862.      *  </p>
  1863.      *  @author Petre Bazavan
  1864.      *  @author Lucian Barbulescu
  1865.      */
  1866.     private class FieldFourierCjSjCoefficients <T extends RealFieldElement<T>> {

  1868.         /** Maximum possible value for j. */
  1869.         private final int jMax;

  1871.         /** The C<sub>i</sub><sup>j</sup> coefficients.
  1872.          * <p>
  1873.          * the index i corresponds to the following elements: <br/>
  1874.          * - 0 for a <br>
  1875.          * - 1 for k <br>
  1876.          * - 2 for h <br>
  1877.          * - 3 for q <br>
  1878.          * - 4 for p <br>
  1879.          * - 5 for λ <br>
  1880.          * </p>
  1881.          */
  1882.         private final T[][] cCoef;

  1884.         /** The C<sub>i</sub><sup>j</sup> coefficients.
  1885.          * <p>
  1886.          * the index i corresponds to the following elements: <br/>
  1887.          * - 0 for a <br>
  1888.          * - 1 for k <br>
  1889.          * - 2 for h <br>
  1890.          * - 3 for q <br>
  1891.          * - 4 for p <br>
  1892.          * - 5 for λ <br>
  1893.          * </p>
  1894.          */
  1895.         private final T[][] sCoef;

  1897.         /** Standard constructor.
  1898.          * @param state the current state
  1899.          * @param jMax maximum value for j
  1900.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1901.          * @param parameters values of the force model parameters
  1902.          * @param field field used by default
  1903.          */
  1904.         FieldFourierCjSjCoefficients(final FieldSpacecraftState<T> state, final int jMax,
  1905.                                      final FieldAuxiliaryElements<T> auxiliaryElements,
  1906.                                      final T[] parameters,
  1907.                                      final Field<T> field) {
  1908.             //Initialise the fields
  1909.             this.jMax = jMax;

  1911.             //Allocate the arrays
  1912.             final int rows = jMax + 1;
  1913.             cCoef = MathArrays.buildArray(field, rows, 6);
  1914.             sCoef = MathArrays.buildArray(field, rows, 6);

  1916.             //Compute the coefficients
  1917.             computeCoefficients(state, auxiliaryElements, parameters, field);
  1918.         }

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

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

  1949.                     //divide by PI and set the values for the coefficients
  1950.                     for (int i = 0; i < 6; i++) {
  1951.                         cCoef[j][i] = curentCoefficients[i].multiply(ooPI);
  1952.                         sCoef[j][i] = curentCoefficients[i + 6].multiply(ooPI);
  1953.                     }
  1954.                 }
  1955.             }
  1956.         }

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

  1967.         /** Get the coefficient S<sub>i</sub><sup>j</sup>.
  1968.          * @param i i index - corresponds to the required variation
  1969.          * @param j j index
  1970.          * @return the coefficient S<sub>i</sub><sup>j</sup>
  1971.          */
  1972.         public T getSij(final int i, final int j) {
  1973.             return sCoef[j][i];
  1974.         }
  1975.     }

  1977.     /** This class handles the short periodic coefficients described in Danielson 2.5.3-26.
  1978.      *
  1979.      * <p>
  1980.      * The value of M is 0. Also, since the values of the Fourier coefficient D<sub>i</sub><sup>m</sup> is 0
  1981.      * then the values of the coefficients D<sub>i</sub><sup>m</sup> for m &gt; 2 are also 0.
  1982.      * </p>
  1983.      * @author Petre Bazavan
  1984.      * @author Lucian Barbulescu
  1985.      *
  1986.      */
  1987.     private static class GaussianShortPeriodicCoefficients implements ShortPeriodTerms {

  1989.         /** Maximum value for j index. */
  1990.         private final int jMax;

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

  1995.         /** Prefix for coefficients keys. */
  1996.         private final String coefficientsKeyPrefix;

  1998.         /** All coefficients slots. */
  1999.         private final transient TimeSpanMap<Slot> slots;

  2001.         /** Constructor.
  2002.          *  @param coefficientsKeyPrefix prefix for coefficients keys
  2003.          *  @param jMax maximum value for j index
  2004.          *  @param interpolationPoints number of points used in the interpolation process
  2005.          *  @param slots all coefficients slots
  2006.          */
  2007.         GaussianShortPeriodicCoefficients(final String coefficientsKeyPrefix,
  2008.                                           final int jMax, final int interpolationPoints,
  2009.                                           final TimeSpanMap<Slot> slots) {
  2010.             //Initialize fields
  2011.             this.jMax                  = jMax;
  2012.             this.interpolationPoints   = interpolationPoints;
  2013.             this.coefficientsKeyPrefix = coefficientsKeyPrefix;
  2014.             this.slots                 = slots;
  2015.         }

  2017.         /** Get the slot valid for some date.
  2018.          * @param meanStates mean states defining the slot
  2019.          * @return slot valid at the specified date
  2020.          */
  2021.         public Slot createSlot(final SpacecraftState... meanStates) {
  2022.             final Slot         slot  = new Slot(jMax, interpolationPoints);
  2023.             final AbsoluteDate first = meanStates[0].getDate();
  2024.             final AbsoluteDate last  = meanStates[meanStates.length - 1].getDate();
  2025.             if (first.compareTo(last) <= 0) {
  2026.                 slots.addValidAfter(slot, first);
  2027.             } else {
  2028.                 slots.addValidBefore(slot, first);
  2029.             }
  2030.             return slot;
  2031.         }

  2033.         /** Compute the short periodic coefficients.
  2034.          *
  2035.          * @param state current state information: date, kinematics, attitude
  2036.          * @param slot coefficients slot
  2037.          * @param fourierCjSj Fourier coefficients
  2038.          * @param uijvij U and V coefficients
  2039.          * @param n Keplerian mean motion
  2040.          * @param a semi major axis
  2041.          */
  2042.         private void computeCoefficients(final SpacecraftState state, final Slot slot,
  2043.                                          final FourierCjSjCoefficients fourierCjSj,
  2044.                                          final UijVijCoefficients uijvij,
  2045.                                          final double n, final double a) {

  2047.             // get the current date
  2048.             final AbsoluteDate date = state.getDate();

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

  2053.             // 1. / n
  2054.             final double oon = 1. / n;
  2055.             // 3. / (2 * a * n)
  2056.             final double to2an = 1.5 * oon / a;
  2057.             // 3. / (4 * a * n)
  2058.             final double to4an = to2an / 2;

  2060.             // Compute the coefficients for each element
  2061.             final int size = jMax + 1;
  2062.             final double[]   di1        = new double[6];
  2063.             final double[]   di2        = new double[6];
  2064.             final double[][] currentCij = new double[size][6];
  2065.             final double[][] currentSij = new double[size][6];
  2066.             for (int i = 0; i < 6; i++) {

  2068.                 // compute D<sub>i</sub>¹ and D<sub>i</sub>² (all others are 0)
  2069.                 di1[i] = -oon * fourierCjSj.getCij(i, 0);
  2070.                 if (i == 5) {
  2071.                     di1[i] += to2an * uijvij.getU1(0, 0);
  2072.                 }
  2073.                 di2[i] = 0.;
  2074.                 if (i == 5) {
  2075.                     di2[i] += -to4an * fourierCjSj.getCij(0, 0);
  2076.                 }

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

  2081.                 for (int j = 1; j <= jMax; j++) {
  2082.                     // compute the current C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup>
  2083.                     currentCij[j][i] = oon * uijvij.getU1(j, i);
  2084.                     if (i == 5) {
  2085.                         currentCij[j][i] += -to2an * uijvij.getU2(j);
  2086.                     }
  2087.                     currentSij[j][i] = oon * uijvij.getV1(j, i);
  2088.                     if (i == 5) {
  2089.                         currentSij[j][i] += -to2an * uijvij.getV2(j);
  2090.                     }

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

  2096.             }

  2098.             // add the values to the interpolators
  2099.             slot.cij[0].addGridPoint(date, currentCij[0]);
  2100.             slot.dij[1].addGridPoint(date, di1);
  2101.             slot.dij[2].addGridPoint(date, di2);
  2102.             for (int j = 1; j <= jMax; j++) {
  2103.                 slot.cij[j].addGridPoint(date, currentCij[j]);
  2104.                 slot.sij[j].addGridPoint(date, currentSij[j]);
  2105.             }

  2107.         }

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

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

  2128.                 //multiply by 2 / k²
  2129.                 currentTerm *= 2. / (kIndex * kIndex);

  2131.                 // add the term to the result
  2132.                 k20 += currentTerm;
  2133.             }

  2135.             return k20;
  2136.         }

  2138.         /** {@inheritDoc} */
  2139.         @Override
  2140.         public double[] value(final Orbit meanOrbit) {

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

  2145.             // Get the True longitude L
  2146.             final double L = meanOrbit.getLv();

  2148.             // Compute the center (l - λ)
  2149.             final double center =  L - meanOrbit.getLM();
  2150.             // Compute (l - λ)²
  2151.             final double center2 = center * center;

  2153.             // Initialize short periodic variations
  2154.             final double[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  2155.             final double[] d1 = slot.dij[1].value(meanOrbit.getDate());
  2156.             final double[] d2 = slot.dij[2].value(meanOrbit.getDate());
  2157.             for (int i = 0; i < 6; i++) {
  2158.                 shortPeriodicVariation[i] += center * d1[i] + center2 * d2[i];
  2159.             }

  2161.             for (int j = 1; j <= JMAX; j++) {
  2162.                 final double[] c = slot.cij[j].value(meanOrbit.getDate());
  2163.                 final double[] s = slot.sij[j].value(meanOrbit.getDate());
  2164.                 final SinCos sc  = FastMath.sinCos(j * L);
  2165.                 final double cos = sc.cos();
  2166.                 final double sin = sc.sin();
  2167.                 for (int i = 0; i < 6; i++) {
  2168.                     // add corresponding term to the short periodic variation
  2169.                     shortPeriodicVariation[i] += c[i] * cos;
  2170.                     shortPeriodicVariation[i] += s[i] * sin;
  2171.                 }
  2172.             }

  2174.             return shortPeriodicVariation;

  2176.         }

  2178.         /** {@inheritDoc} */
  2179.         public String getCoefficientsKeyPrefix() {
  2180.             return coefficientsKeyPrefix;
  2181.         }

  2183.         /** {@inheritDoc}
  2184.          * <p>
  2185.          * For Gaussian forces, there are JMAX cj coefficients,
  2186.          * JMAX sj coefficients and 3 dj coefficients. As JMAX = 12,
  2187.          * this sums up to 27 coefficients. The j index is the integer
  2188.          * multiplier for the true longitude argument in the cj and sj
  2189.          * coefficients and to the degree in  the polynomial dj coefficients.
  2190.          * </p>
  2191.          */
  2192.         @Override
  2193.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected) {

  2195.             // select the coefficients slot
  2196.             final Slot slot = slots.get(date);

  2198.             final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * JMAX + 3);
  2199.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  2200.             storeIfSelected(coefficients, selected, slot.dij[1].value(date), "d", 1);
  2201.             storeIfSelected(coefficients, selected, slot.dij[2].value(date), "d", 2);
  2202.             for (int j = 1; j <= JMAX; j++) {
  2203.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  2204.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  2205.             }

  2207.             return coefficients;

  2209.         }

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

  2232.     }

  2234.      /** This class handles the short periodic coefficients described in Danielson 2.5.3-26.
  2235.      *
  2236.      * <p>
  2237.      * The value of M is 0. Also, since the values of the Fourier coefficient D<sub>i</sub><sup>m</sup> is 0
  2238.      * then the values of the coefficients D<sub>i</sub><sup>m</sup> for m &gt; 2 are also 0.
  2239.      * </p>
  2240.      * @author Petre Bazavan
  2241.      * @author Lucian Barbulescu
  2242.      *
  2243.      */
  2244.     private static class FieldGaussianShortPeriodicCoefficients <T extends RealFieldElement<T>> implements FieldShortPeriodTerms<T> {

  2246.         /** Maximum value for j index. */
  2247.         private final int jMax;

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

  2252.         /** Prefix for coefficients keys. */
  2253.         private final String coefficientsKeyPrefix;

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

  2258.         /** Constructor.
  2259.          *  @param coefficientsKeyPrefix prefix for coefficients keys
  2260.          *  @param jMax maximum value for j index
  2261.          *  @param interpolationPoints number of points used in the interpolation process
  2262.          *  @param slots all coefficients slots
  2263.          */
  2264.         FieldGaussianShortPeriodicCoefficients(final String coefficientsKeyPrefix,
  2265.                                                final int jMax, final int interpolationPoints,
  2266.                                                final FieldTimeSpanMap<FieldSlot<T>, T> slots) {
  2267.             //Initialize fields
  2268.             this.jMax                  = jMax;
  2269.             this.interpolationPoints   = interpolationPoints;
  2270.             this.coefficientsKeyPrefix = coefficientsKeyPrefix;
  2271.             this.slots                 = slots;
  2272.         }

  2274.         /** Get the slot valid for some date.
  2275.          * @param meanStates mean states defining the slot
  2276.          * @return slot valid at the specified date
  2277.          */
  2278.         @SuppressWarnings("unchecked")
  2279.         public FieldSlot<T> createSlot(final FieldSpacecraftState<T>... meanStates) {
  2280.             final FieldSlot<T>         slot  = new FieldSlot<>(jMax, interpolationPoints);
  2281.             final FieldAbsoluteDate<T> first = meanStates[0].getDate();
  2282.             final FieldAbsoluteDate<T> last  = meanStates[meanStates.length - 1].getDate();
  2283.             if (first.compareTo(last) <= 0) {
  2284.                 slots.addValidAfter(slot, first);
  2285.             } else {
  2286.                 slots.addValidBefore(slot, first);
  2287.             }
  2288.             return slot;
  2289.         }

  2291.         /** Compute the short periodic coefficients.
  2292.          *
  2293.          * @param state current state information: date, kinematics, attitude
  2294.          * @param slot coefficients slot
  2295.          * @param fourierCjSj Fourier coefficients
  2296.          * @param uijvij U and V coefficients
  2297.          * @param n Keplerian mean motion
  2298.          * @param a semi major axis
  2299.          * @param field field used by default
  2300.          */
  2301.         private void computeCoefficients(final FieldSpacecraftState<T> state, final FieldSlot<T> slot,
  2302.                                          final FieldFourierCjSjCoefficients<T> fourierCjSj,
  2303.                                          final FieldUijVijCoefficients<T> uijvij,
  2304.                                          final T n, final T a,
  2305.                                          final Field<T> field) {

  2307.             // Zero
  2308.             final T zero = field.getZero();

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

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

  2316.             // 1. / n
  2317.             final T oon = n.reciprocal();
  2318.             // 3. / (2 * a * n)
  2319.             final T to2an = oon.multiply(1.5).divide(a);
  2320.             // 3. / (4 * a * n)
  2321.             final T to4an = to2an.divide(2.);

  2323.             // Compute the coefficients for each element
  2324.             final int size = jMax + 1;
  2325.             final T[]   di1        = MathArrays.buildArray(field, 6);
  2326.             final T[]   di2        = MathArrays.buildArray(field, 6);
  2327.             final T[][] currentCij = MathArrays.buildArray(field, size, 6);
  2328.             final T[][] currentSij = MathArrays.buildArray(field, size, 6);
  2329.             for (int i = 0; i < 6; i++) {

  2331.                 // compute D<sub>i</sub>¹ and D<sub>i</sub>² (all others are 0)
  2332.                 di1[i] = oon.negate().multiply(fourierCjSj.getCij(i, 0));
  2333.                 if (i == 5) {
  2334.                     di1[i] = di1[i].add(to2an.multiply(uijvij.getU1(0, 0)));
  2335.                 }
  2336.                 di2[i] = zero;
  2337.                 if (i == 5) {
  2338.                     di2[i] = di2[i].add(to4an.negate().multiply(fourierCjSj.getCij(0, 0)));
  2339.                 }

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

  2344.                 for (int j = 1; j <= jMax; j++) {
  2345.                     // compute the current C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup>
  2346.                     currentCij[j][i] = oon.multiply(uijvij.getU1(j, i));
  2347.                     if (i == 5) {
  2348.                         currentCij[j][i] = currentCij[j][i].add(to2an.negate().multiply(uijvij.getU2(j)));
  2349.                     }
  2350.                     currentSij[j][i] = oon.multiply(uijvij.getV1(j, i));
  2351.                     if (i == 5) {
  2352.                         currentSij[j][i] = currentSij[j][i].add(to2an.negate().multiply(uijvij.getV2(j)));
  2353.                     }

  2355.                     // add the computed coefficients to C<sub>i</sub>⁰
  2356.                     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());
  2357.                 }

  2359.             }

  2361.             // add the values to the interpolators
  2362.             slot.cij[0].addGridPoint(date, currentCij[0]);
  2363.             slot.dij[1].addGridPoint(date, di1);
  2364.             slot.dij[2].addGridPoint(date, di2);
  2365.             for (int j = 1; j <= jMax; j++) {
  2366.                 slot.cij[j].addGridPoint(date, currentCij[j]);
  2367.                 slot.sij[j].addGridPoint(date, currentSij[j]);
  2368.             }

  2370.         }

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

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

  2393.                 //multiply by 2 / k²
  2394.                 currentTerm = currentTerm.multiply(2. / (kIndex * kIndex));

  2396.                 // add the term to the result
  2397.                 k20 = k20.add(currentTerm);
  2398.             }

  2400.             return k20;
  2401.         }

  2403.         /** {@inheritDoc} */
  2404.         @Override
  2405.         public T[] value(final FieldOrbit<T> meanOrbit) {

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

  2410.             // Get the True longitude L
  2411.             final T L = meanOrbit.getLv();

  2413.             // Compute the center (l - λ)
  2414.             final T center =  L.subtract(meanOrbit.getLM());
  2415.             // Compute (l - λ)²
  2416.             final T center2 = center.multiply(center);

  2418.             // Initialize short periodic variations
  2419.             final T[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  2420.             final T[] d1 = slot.dij[1].value(meanOrbit.getDate());
  2421.             final T[] d2 = slot.dij[2].value(meanOrbit.getDate());
  2422.             for (int i = 0; i < 6; i++) {
  2423.                 shortPeriodicVariation[i] = shortPeriodicVariation[i].add(center.multiply(d1[i]).add(center2.multiply(d2[i])));
  2424.             }

  2426.             for (int j = 1; j <= JMAX; j++) {
  2427.                 final T[] c = slot.cij[j].value(meanOrbit.getDate());
  2428.                 final T[] s = slot.sij[j].value(meanOrbit.getDate());
  2429.                 final FieldSinCos<T> sc = FastMath.sinCos(L.multiply(j));
  2430.                 final T cos = sc.cos();
  2431.                 final T sin = sc.sin();
  2432.                 for (int i = 0; i < 6; i++) {
  2433.                     // add corresponding term to the short periodic variation
  2434.                     shortPeriodicVariation[i] = shortPeriodicVariation[i].add(c[i].multiply(cos));
  2435.                     shortPeriodicVariation[i] = shortPeriodicVariation[i].add(s[i].multiply(sin));
  2436.                 }
  2437.             }

  2439.             return shortPeriodicVariation;

  2441.         }

  2443.         /** {@inheritDoc} */
  2444.         public String getCoefficientsKeyPrefix() {
  2445.             return coefficientsKeyPrefix;
  2446.         }

  2448.         /** {@inheritDoc}
  2449.          * <p>
  2450.          * For Gaussian forces, there are JMAX cj coefficients,
  2451.          * JMAX sj coefficients and 3 dj coefficients. As JMAX = 12,
  2452.          * this sums up to 27 coefficients. The j index is the integer
  2453.          * multiplier for the true longitude argument in the cj and sj
  2454.          * coefficients and to the degree in  the polynomial dj coefficients.
  2455.          * </p>
  2456.          */
  2457.         @Override
  2458.         public Map<String, T[]> getCoefficients(final FieldAbsoluteDate<T> date, final Set<String> selected) {

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

  2463.             final Map<String, T[]> coefficients = new HashMap<String, T[]>(2 * JMAX + 3);
  2464.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  2465.             storeIfSelected(coefficients, selected, slot.dij[1].value(date), "d", 1);
  2466.             storeIfSelected(coefficients, selected, slot.dij[2].value(date), "d", 2);
  2467.             for (int j = 1; j <= JMAX; j++) {
  2468.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  2469.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  2470.             }

  2472.             return coefficients;

  2474.         }

  2476.         /** Put a coefficient in a map if selected.
  2477.          * @param map map to populate
  2478.          * @param selected set of coefficients that should be put in the map
  2479.          * (empty set means all coefficients are selected)
  2480.          * @param value coefficient value
  2481.          * @param id coefficient identifier
  2482.          * @param indices list of coefficient indices
  2483.          */
  2484.         private void storeIfSelected(final Map<String, T[]> map, final Set<String> selected,
  2485.                                      final T[] value, final String id, final int... indices) {
  2486.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  2487.             keyBuilder.append(id);
  2488.             for (int index : indices) {
  2489.                 keyBuilder.append('[').append(index).append(']');
  2490.             }
  2491.             final String key = keyBuilder.toString();
  2492.             if (selected.isEmpty() || selected.contains(key)) {
  2493.                 map.put(key, value);
  2494.             }
  2495.         }

  2497.     }


  2500.     /** The U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients described by
  2501.      * equations 2.5.3-(21) and 2.5.3-(22) from Danielson.
  2502.      * <p>
  2503.      * The index i takes only the values 1 and 2<br>
  2504.      * For U only the index 0 for j is used.
  2505.      * </p>
  2506.      *
  2507.      * @author Petre Bazavan
  2508.      * @author Lucian Barbulescu
  2509.      */
  2510.     private static class UijVijCoefficients {

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

  2523.         /** The V₁<sup>j</sup> coefficients.
  2524.          * <p>
  2525.          * The first index identifies the Fourier coefficients used<br>
  2526.          * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
  2527.          * for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
  2528.          * to compute the coefficients V₂<sup>j</sup>
  2529.          * </p>
  2530.          */
  2531.         private final double[][] v1ij;

  2533.         /** The U₂<sup>j</sup> coefficients.
  2534.          * <p>
  2535.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
  2536.          * </p>
  2537.          */
  2538.         private final double[] u2ij;

  2540.         /** The V₂<sup>j</sup> coefficients.
  2541.          * <p>
  2542.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
  2543.          * </p>
  2544.          */
  2545.         private final double[] v2ij;

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

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

  2553.         /** The maximum value for j index. */
  2554.         private final int jMax;

  2556.         /** Constructor.
  2557.          * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
  2558.          * @param fourierCjSj the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  2559.          * @param jMax maximum value for j index
  2560.          */
  2561.         UijVijCoefficients(final double[][] currentRhoSigmaj, final FourierCjSjCoefficients fourierCjSj, final int jMax) {
  2562.             this.currentRhoSigmaj = currentRhoSigmaj;
  2563.             this.fourierCjSj = fourierCjSj;
  2564.             this.jMax = jMax;

  2566.             // initialize the internal arrays.
  2567.             this.u1ij = new double[6][2 * jMax + 1];
  2568.             this.v1ij = new double[6][2 * jMax + 1];
  2569.             this.u2ij = new double[jMax + 1];
  2570.             this.v2ij = new double[jMax + 1];

  2572.             //compute the coefficients
  2573.             computeU1V1Coefficients();
  2574.             computeU2V2Coefficients();
  2575.         }

  2577.         /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients. */
  2578.         private void computeU1V1Coefficients() {
  2579.             // generate the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients
  2580.             // for j >= 1
  2581.             // also the U₁⁰ for Fourier index = 1 (i == 0) coefficient will be computed
  2582.             u1ij[0][0] = 0;
  2583.             for (int j = 1; j <= jMax; j++) {
  2584.                 // compute 1 / j
  2585.                 final double ooj = 1. / j;

  2587.                 for (int i = 0; i < 6; i++) {
  2588.                     //j is aready between 1 and J
  2589.                     u1ij[i][j] = fourierCjSj.getSij(i, j);
  2590.                     v1ij[i][j] = fourierCjSj.getCij(i, j);

  2592.                     // 1 - δ<sub>1j</sub> is 1 for all j > 1
  2593.                     if (j > 1) {
  2594.                         // k starts with 1 because j-J is less than or equal to 0
  2595.                         for (int kIndex = 1; kIndex <= j - 1; kIndex++) {
  2596.                             // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
  2597.                             // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
  2598.                             u1ij[i][j] +=   fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[1][kIndex] +
  2599.                                             fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[0][kIndex];

  2601.                             // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  2602.                             // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  2603.                             v1ij[i][j] +=   fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[0][kIndex] -
  2604.                                             fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[1][kIndex];
  2605.                         }
  2606.                     }

  2608.                     // since j must be between 1 and J-1 and is already between 1 and J
  2609.                     // the following sum is skiped only for j = jMax
  2610.                     if (j != jMax) {
  2611.                         for (int kIndex = 1; kIndex <= jMax - j; kIndex++) {
  2612.                             // -C<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub> +
  2613.                             // S<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub>
  2614.                             u1ij[i][j] +=   -fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[1][kIndex] +
  2615.                                             fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[0][kIndex];

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

  2624.                     for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  2625.                         // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  2626.                         // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  2627.                         u1ij[i][j] +=   -fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[1][j + kIndex] -
  2628.                                         fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[0][j + kIndex];

  2630.                         // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  2631.                         // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  2632.                         v1ij[i][j] +=   fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[0][j + kIndex] +
  2633.                                         fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[1][j + kIndex];
  2634.                     }

  2636.                     // divide by 1 / j
  2637.                     u1ij[i][j] *= -ooj;
  2638.                     v1ij[i][j] *= ooj;

  2640.                     // if index = 1 (i == 0) add the computed terms to U₁⁰
  2641.                     if (i == 0) {
  2642.                         //- (U₁<sup>j</sup> * ρ<sub>j</sub> + V₁<sup>j</sup> * σ<sub>j</sub>
  2643.                         u1ij[0][0] += -u1ij[0][j] * currentRhoSigmaj[0][j] - v1ij[0][j] * currentRhoSigmaj[1][j];
  2644.                     }
  2645.                 }
  2646.             }

  2648.             // Terms with j > jMax are required only when computing the coefficients
  2649.             // U₂<sup>j</sup> and V₂<sup>j</sup>
  2650.             // and those coefficients are only required for Fourier index = 1 (i == 0).
  2651.             for (int j = jMax + 1; j <= 2 * jMax; j++) {
  2652.                 // compute 1 / j
  2653.                 final double ooj = 1. / j;
  2654.                 //the value of i is 0
  2655.                 u1ij[0][j] = 0.;
  2656.                 v1ij[0][j] = 0.;

  2658.                 //k starts from j-J as it is always greater than or equal to 1
  2659.                 for (int kIndex = j - jMax; kIndex <= j - 1; kIndex++) {
  2660.                     // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
  2661.                     // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
  2662.                     u1ij[0][j] +=   fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[1][kIndex] +
  2663.                                     fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[0][kIndex];

  2665.                     // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  2666.                     // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  2667.                     v1ij[0][j] +=   fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[0][kIndex] -
  2668.                                     fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[1][kIndex];
  2669.                 }
  2670.                 for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  2671.                     // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  2672.                     // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  2673.                     u1ij[0][j] +=   -fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[1][j + kIndex] -
  2674.                                     fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[0][j + kIndex];

  2676.                     // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  2677.                     // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  2678.                     v1ij[0][j] +=   fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[0][j + kIndex] +
  2679.                                     fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[1][j + kIndex];
  2680.                 }

  2682.                 // divide by 1 / j
  2683.                 u1ij[0][j] *= -ooj;
  2684.                 v1ij[0][j] *= ooj;
  2685.             }
  2686.         }

  2688.         /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
  2689.          * <p>
  2690.          * Only the coefficients for Fourier index = 1 (i == 0) are required.
  2691.          * </p>
  2692.          */
  2693.         private void computeU2V2Coefficients() {
  2694.             for (int j = 1; j <= jMax; j++) {
  2695.                 // compute 1 / j
  2696.                 final double ooj = 1. / j;

  2698.                 // only the values for i == 0 are computed
  2699.                 u2ij[j] = v1ij[0][j];
  2700.                 v2ij[j] = u1ij[0][j];

  2702.                 // 1 - δ<sub>1j</sub> is 1 for all j > 1
  2703.                 if (j > 1) {
  2704.                     for (int l = 1; l <= j - 1; l++) {
  2705.                         // U₁<sup>j-l</sup> * σ<sub>l</sub> +
  2706.                         // V₁<sup>j-l</sup> * ρ<sub>l</sub>
  2707.                         u2ij[j] +=   u1ij[0][j - l] * currentRhoSigmaj[1][l] +
  2708.                                      v1ij[0][j - l] * currentRhoSigmaj[0][l];

  2710.                         // U₁<sup>j-l</sup> * ρ<sub>l</sub> -
  2711.                         // V₁<sup>j-l</sup> * σ<sub>l</sub>
  2712.                         v2ij[j] +=   u1ij[0][j - l] * currentRhoSigmaj[0][l] -
  2713.                                      v1ij[0][j - l] * currentRhoSigmaj[1][l];
  2714.                     }
  2715.                 }

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

  2727.                     // U₁<sup>j+l</sup> * ρ<sub>l</sub> +
  2728.                     // U₁<sup>l</sup> * ρ<sub>j+l</sub> +
  2729.                     // V₁<sup>j+l</sup> * σ<sub>l</sub> +
  2730.                     // V₁<sup>l</sup> * σ<sub>j+l</sub>
  2731.                     u2ij[j] +=   u1ij[0][j + l] * currentRhoSigmaj[0][l] +
  2732.                                  u1ij[0][l] * currentRhoSigmaj[0][j + l] +
  2733.                                  v1ij[0][j + l] * currentRhoSigmaj[1][l] +
  2734.                                  v1ij[0][l] * currentRhoSigmaj[1][j + l];
  2735.                 }

  2737.                 // divide by 1 / j
  2738.                 u2ij[j] *= -ooj;
  2739.                 v2ij[j] *= ooj;
  2740.             }
  2741.         }

  2743.         /** Get the coefficient U₁<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 U₁<sup>j</sup> for the given Fourier index i
  2748.          */
  2749.         public double getU1(final int j, final int i) {
  2750.             return u1ij[i][j];
  2751.         }

  2753.         /** Get the coefficient V₁<sup>j</sup> for Fourier index i.
  2754.          *
  2755.          * @param j j index
  2756.          * @param i Fourier index (starts at 0)
  2757.          * @return the coefficient V₁<sup>j</sup> for the given Fourier index i
  2758.          */
  2759.         public double getV1(final int j, final int i) {
  2760.             return v1ij[i][j];
  2761.         }

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

  2772.         /** Get the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0).
  2773.          *
  2774.          * @param j j index
  2775.          * @return the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0)
  2776.          */
  2777.         public double getV2(final int j) {
  2778.             return v2ij[j];
  2779.         }
  2780.     }

  2782.     /** The U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients described by
  2783.      * equations 2.5.3-(21) and 2.5.3-(22) from Danielson.
  2784.      * <p>
  2785.      * The index i takes only the values 1 and 2<br>
  2786.      * For U only the index 0 for j is used.
  2787.      * </p>
  2788.      *
  2789.      * @author Petre Bazavan
  2790.      * @author Lucian Barbulescu
  2791.      */
  2792.     private static class FieldUijVijCoefficients <T extends RealFieldElement<T>> {

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

  2805.         /** The V₁<sup>j</sup> coefficients.
  2806.          * <p>
  2807.          * The first index identifies the Fourier coefficients used<br>
  2808.          * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
  2809.          * for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
  2810.          * to compute the coefficients V₂<sup>j</sup>
  2811.          * </p>
  2812.          */
  2813.         private final T[][] v1ij;

  2815.         /** The U₂<sup>j</sup> coefficients.
  2816.          * <p>
  2817.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
  2818.          * </p>
  2819.          */
  2820.         private final T[] u2ij;

  2822.         /** The V₂<sup>j</sup> coefficients.
  2823.          * <p>
  2824.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
  2825.          * </p>
  2826.          */
  2827.         private final T[] v2ij;

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

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

  2835.         /** The maximum value for j index. */
  2836.         private final int jMax;

  2838.         /** Constructor.
  2839.          * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
  2840.          * @param fourierCjSj the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  2841.          * @param jMax maximum value for j index
  2842.          * @param field field used by default
  2843.          */
  2844.         FieldUijVijCoefficients(final T[][] currentRhoSigmaj,
  2845.                                 final FieldFourierCjSjCoefficients<T> fourierCjSj,
  2846.                                 final int jMax,
  2847.                                 final Field<T> field) {
  2848.             this.currentRhoSigmaj = currentRhoSigmaj;
  2849.             this.fourierCjSj = fourierCjSj;
  2850.             this.jMax = jMax;

  2852.             // initialize the internal arrays.
  2853.             this.u1ij = MathArrays.buildArray(field, 6, 2 * jMax + 1);
  2854.             this.v1ij = MathArrays.buildArray(field, 6, 2 * jMax + 1);
  2855.             this.u2ij = MathArrays.buildArray(field, jMax + 1);
  2856.             this.v2ij = MathArrays.buildArray(field, jMax + 1);

  2858.             //compute the coefficients
  2859.             computeU1V1Coefficients(field);
  2860.             computeU2V2Coefficients(field);
  2861.         }

  2863.         /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
  2864.          * @param field field used by default
  2865.          */
  2866.         private void computeU1V1Coefficients(final Field<T> field) {
  2867.             // Zero
  2868.             final T zero = field.getZero();

  2870.             // generate the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients
  2871.             // for j >= 1
  2872.             // also the U₁⁰ for Fourier index = 1 (i == 0) coefficient will be computed
  2873.             u1ij[0][0] = zero;
  2874.             for (int j = 1; j <= jMax; j++) {
  2875.                 // compute 1 / j
  2876.                 final double ooj = 1. / j;

  2878.                 for (int i = 0; i < 6; i++) {
  2879.                     //j is aready between 1 and J
  2880.                     u1ij[i][j] = fourierCjSj.getSij(i, j);
  2881.                     v1ij[i][j] = fourierCjSj.getCij(i, j);

  2883.                     // 1 - δ<sub>1j</sub> is 1 for all j > 1
  2884.                     if (j > 1) {
  2885.                         // k starts with 1 because j-J is less than or equal to 0
  2886.                         for (int kIndex = 1; kIndex <= j - 1; kIndex++) {
  2887.                             // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
  2888.                             // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
  2889.                             u1ij[i][j] = u1ij[i][j].add(fourierCjSj.getCij(i, j - kIndex).multiply(currentRhoSigmaj[1][kIndex]).
  2890.                                          add(fourierCjSj.getSij(i, j - kIndex).multiply(currentRhoSigmaj[0][kIndex])));

  2892.                             // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  2893.                             // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  2894.                             v1ij[i][j] = v1ij[i][j].add(fourierCjSj.getCij(i, j - kIndex).multiply(currentRhoSigmaj[0][kIndex]).
  2895.                                          subtract(fourierCjSj.getSij(i, j - kIndex).multiply(currentRhoSigmaj[1][kIndex])));
  2896.                         }
  2897.                     }

  2899.                     // since j must be between 1 and J-1 and is already between 1 and J
  2900.                     // the following sum is skiped only for j = jMax
  2901.                     if (j != jMax) {
  2902.                         for (int kIndex = 1; kIndex <= jMax - j; kIndex++) {
  2903.                             // -C<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub> +
  2904.                             // S<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub>
  2905.                             u1ij[i][j] = u1ij[i][j].add(fourierCjSj.getCij(i, j + kIndex).negate().multiply(currentRhoSigmaj[1][kIndex]).
  2906.                                          add(fourierCjSj.getSij(i, j + kIndex).multiply(currentRhoSigmaj[0][kIndex])));

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

  2915.                     for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  2916.                         // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  2917.                         // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  2918.                         u1ij[i][j] = u1ij[i][j].add(fourierCjSj.getCij(i, kIndex).negate().multiply(currentRhoSigmaj[1][j + kIndex]).
  2919.                                      subtract(fourierCjSj.getSij(i, kIndex).multiply(currentRhoSigmaj[0][j + kIndex])));

  2921.                         // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  2922.                         // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  2923.                         v1ij[i][j] = v1ij[i][j].add(fourierCjSj.getCij(i, kIndex).multiply(currentRhoSigmaj[0][j + kIndex]).
  2924.                                      add(fourierCjSj.getSij(i, kIndex).multiply(currentRhoSigmaj[1][j + kIndex])));
  2925.                     }

  2927.                     // divide by 1 / j
  2928.                     u1ij[i][j] = u1ij[i][j].multiply(-ooj);
  2929.                     v1ij[i][j] = v1ij[i][j].multiply(ooj);

  2931.                     // if index = 1 (i == 0) add the computed terms to U₁⁰
  2932.                     if (i == 0) {
  2933.                         //- (U₁<sup>j</sup> * ρ<sub>j</sub> + V₁<sup>j</sup> * σ<sub>j</sub>
  2934.                         u1ij[0][0] = u1ij[0][0].add(u1ij[0][j].negate().multiply(currentRhoSigmaj[0][j]).subtract(v1ij[0][j].multiply(currentRhoSigmaj[1][j])));
  2935.                     }
  2936.                 }
  2937.             }

  2939.             // Terms with j > jMax are required only when computing the coefficients
  2940.             // U₂<sup>j</sup> and V₂<sup>j</sup>
  2941.             // and those coefficients are only required for Fourier index = 1 (i == 0).
  2942.             for (int j = jMax + 1; j <= 2 * jMax; j++) {
  2943.                 // compute 1 / j
  2944.                 final double ooj = 1. / j;
  2945.                 //the value of i is 0
  2946.                 u1ij[0][j] = zero;
  2947.                 v1ij[0][j] = zero;

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

  2956.                     // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  2957.                     // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  2958.                     v1ij[0][j] = v1ij[0][j].add(fourierCjSj.getCij(0, j - kIndex).multiply(currentRhoSigmaj[0][kIndex]).
  2959.                                  subtract(fourierCjSj.getSij(0, j - kIndex).multiply(currentRhoSigmaj[1][kIndex])));
  2960.                 }
  2961.                 for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  2962.                     // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  2963.                     // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  2964.                     u1ij[0][j] = u1ij[0][j].add(fourierCjSj.getCij(0, kIndex).negate().multiply(currentRhoSigmaj[1][j + kIndex]).
  2965.                                  subtract(fourierCjSj.getSij(0, kIndex).multiply(currentRhoSigmaj[0][j + kIndex])));

  2967.                     // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  2968.                     // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  2969.                     v1ij[0][j] = v1ij[0][j].add(fourierCjSj.getCij(0, kIndex).multiply(currentRhoSigmaj[0][j + kIndex]).add(
  2970.                                  fourierCjSj.getSij(0, kIndex).multiply(currentRhoSigmaj[1][j + kIndex])));
  2971.                 }

  2973.                 // divide by 1 / j
  2974.                 u1ij[0][j] = u1ij[0][j].multiply(-ooj);
  2975.                 v1ij[0][j] = v1ij[0][j].multiply(ooj);
  2976.             }
  2977.         }

  2979.         /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
  2980.          * <p>
  2981.          * Only the coefficients for Fourier index = 1 (i == 0) are required.
  2982.          * </p>
  2983.          * @param field field used by default
  2984.          */
  2985.         private void computeU2V2Coefficients(final Field<T> field) {
  2986.             for (int j = 1; j <= jMax; j++) {
  2987.                 // compute 1 / j
  2988.                 final double ooj = 1. / j;

  2990.                 // only the values for i == 0 are computed
  2991.                 u2ij[j] = v1ij[0][j];
  2992.                 v2ij[j] = u1ij[0][j];

  2994.                 // 1 - δ<sub>1j</sub> is 1 for all j > 1
  2995.                 if (j > 1) {
  2996.                     for (int l = 1; l <= j - 1; l++) {
  2997.                         // U₁<sup>j-l</sup> * σ<sub>l</sub> +
  2998.                         // V₁<sup>j-l</sup> * ρ<sub>l</sub>
  2999.                         u2ij[j] = u2ij[j].add(u1ij[0][j - l].multiply(currentRhoSigmaj[1][l]).
  3000.                                   add(v1ij[0][j - l].multiply(currentRhoSigmaj[0][l])));

  3002.                         // U₁<sup>j-l</sup> * ρ<sub>l</sub> -
  3003.                         // V₁<sup>j-l</sup> * σ<sub>l</sub>
  3004.                         v2ij[j] = v2ij[j].add(u1ij[0][j - l].multiply(currentRhoSigmaj[0][l]).
  3005.                                   subtract(v1ij[0][j - l].multiply(currentRhoSigmaj[1][l])));
  3006.                     }
  3007.                 }

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

  3019.                     // U₁<sup>j+l</sup> * ρ<sub>l</sub> +
  3020.                     // U₁<sup>l</sup> * ρ<sub>j+l</sub> +
  3021.                     // V₁<sup>j+l</sup> * σ<sub>l</sub> +
  3022.                     // V₁<sup>l</sup> * σ<sub>j+l</sub>
  3023.                     u2ij[j] = u2ij[j].add(u1ij[0][j + l].multiply(currentRhoSigmaj[0][l]).
  3024.                               add(u1ij[0][l].multiply(currentRhoSigmaj[0][j + l])).
  3025.                               add(v1ij[0][j + l].multiply(currentRhoSigmaj[1][l])).
  3026.                               add(v1ij[0][l].multiply(currentRhoSigmaj[1][j + l])));
  3027.                 }

  3029.                 // divide by 1 / j
  3030.                 u2ij[j] = u2ij[j].multiply(-ooj);
  3031.                 v2ij[j] = v2ij[j].multiply(ooj);
  3032.             }
  3033.         }

  3035.         /** Get the coefficient U₁<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 U₁<sup>j</sup> for the given Fourier index i
  3040.          */
  3041.         public T getU1(final int j, final int i) {
  3042.             return u1ij[i][j];
  3043.         }

  3045.         /** Get the coefficient V₁<sup>j</sup> for Fourier index i.
  3046.          *
  3047.          * @param j j index
  3048.          * @param i Fourier index (starts at 0)
  3049.          * @return the coefficient V₁<sup>j</sup> for the given Fourier index i
  3050.          */
  3051.         public T getV1(final int j, final int i) {
  3052.             return v1ij[i][j];
  3053.         }

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

  3064.         /** Get the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0).
  3065.          *
  3066.          * @param j j index
  3067.          * @return the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0)
  3068.          */
  3069.         public T getV2(final int j) {
  3070.             return v2ij[j];
  3071.         }
  3072.     }

  3074.     /** Coefficients valid for one time slot. */
  3075.     private static class Slot {

  3077.         /**The coefficients D<sub>i</sub><sup>j</sup>.
  3078.          * <p>
  3079.          * Only for j = 1 and j = 2 the coefficients are not 0. <br>
  3080.          * i corresponds to the equinoctial element, as follows:
  3081.          * - i=0 for a <br/>
  3082.          * - i=1 for k <br/>
  3083.          * - i=2 for h <br/>
  3084.          * - i=3 for q <br/>
  3085.          * - i=4 for p <br/>
  3086.          * - i=5 for λ <br/>
  3087.          * </p>
  3088.          */
  3089.         private final ShortPeriodicsInterpolatedCoefficient[] dij;

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

  3105.         /** The coefficients S<sub>i</sub><sup>j</sup>.
  3106.          * <p>
  3107.          * The index order is sij[j][i] <br/>
  3108.          * i corresponds to the equinoctial element, as follows: <br/>
  3109.          * - i=0 for a <br/>
  3110.          * - i=1 for k <br/>
  3111.          * - i=2 for h <br/>
  3112.          * - i=3 for q <br/>
  3113.          * - i=4 for p <br/>
  3114.          * - i=5 for λ <br/>
  3115.          * </p>
  3116.          */
  3117.         private final ShortPeriodicsInterpolatedCoefficient[] sij;

  3119.         /** Simple constructor.
  3120.          *  @param jMax maximum value for j index
  3121.          *  @param interpolationPoints number of points used in the interpolation process
  3122.          */
  3123.         Slot(final int jMax, final int interpolationPoints) {

  3125.             dij = new ShortPeriodicsInterpolatedCoefficient[3];
  3126.             cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
  3127.             sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];

  3129.             // 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
  3130.             for (int j = 0; j <= jMax; j++) {
  3131.                 cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3132.                 if (j > 0) {
  3133.                     sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3134.                 }
  3135.                 // Initialize only the non-zero D<sub>i</sub><sup>j</sup> coefficients
  3136.                 if (j == 1 || j == 2) {
  3137.                     dij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3138.                 }
  3139.             }

  3141.         }

  3143.     }

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

  3148.         /**The coefficients D<sub>i</sub><sup>j</sup>.
  3149.          * <p>
  3150.          * Only for j = 1 and j = 2 the coefficients are not 0. <br>
  3151.          * i corresponds to the equinoctial element, as follows:
  3152.          * - i=0 for a <br/>
  3153.          * - i=1 for k <br/>
  3154.          * - i=2 for h <br/>
  3155.          * - i=3 for q <br/>
  3156.          * - i=4 for p <br/>
  3157.          * - i=5 for λ <br/>
  3158.          * </p>
  3159.          */
  3160.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] dij;

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

  3176.         /** The coefficients S<sub>i</sub><sup>j</sup>.
  3177.          * <p>
  3178.          * The index order is sij[j][i] <br/>
  3179.          * i corresponds to the equinoctial element, as follows: <br/>
  3180.          * - i=0 for a <br/>
  3181.          * - i=1 for k <br/>
  3182.          * - i=2 for h <br/>
  3183.          * - i=3 for q <br/>
  3184.          * - i=4 for p <br/>
  3185.          * - i=5 for λ <br/>
  3186.          * </p>
  3187.          */
  3188.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] sij;

  3190.         /** Simple constructor.
  3191.          *  @param jMax maximum value for j index
  3192.          *  @param interpolationPoints number of points used in the interpolation process
  3193.          */
  3194.         @SuppressWarnings("unchecked")
  3195.         FieldSlot(final int jMax, final int interpolationPoints) {

  3197.             dij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, 3);
  3198.             cij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);
  3199.             sij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);

  3201.             // 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
  3202.             for (int j = 0; j <= jMax; j++) {
  3203.                 cij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3204.                 if (j > 0) {
  3205.                     sij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3206.                 }
  3207.                 // Initialize only the non-zero D<sub>i</sub><sup>j</sup> coefficients
  3208.                 if (j == 1 || j == 2) {
  3209.                     dij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3210.                 }
  3211.             }

  3213.         }

  3215.     }

  3217. }
Created with JaCoCo 0.8.6.202009150832