AbstractGaussianContribution.java

  1. /* Copyright 2002-2025 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 org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.Field;
  21. import org.hipparchus.analysis.CalculusFieldUnivariateVectorFunction;
  22. import org.hipparchus.analysis.UnivariateVectorFunction;
  23. import org.hipparchus.geometry.euclidean.threed.FieldRotation;
  24. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  25. import org.hipparchus.geometry.euclidean.threed.Rotation;
  26. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  27. import org.hipparchus.util.FastMath;
  28. import org.hipparchus.util.FieldSinCos;
  29. import org.hipparchus.util.MathArrays;
  30. import org.hipparchus.util.SinCos;
  31. import org.orekit.attitudes.Attitude;
  32. import org.orekit.attitudes.AttitudeProvider;
  33. import org.orekit.attitudes.FieldAttitude;
  34. import org.orekit.forces.ForceModel;
  35. import org.orekit.orbits.EquinoctialOrbit;
  36. import org.orekit.orbits.FieldEquinoctialOrbit;
  37. import org.orekit.orbits.FieldOrbit;
  38. import org.orekit.orbits.Orbit;
  39. import org.orekit.orbits.OrbitType;
  40. import org.orekit.orbits.PositionAngleType;
  41. import org.orekit.propagation.FieldSpacecraftState;
  42. import org.orekit.propagation.PropagationType;
  43. import org.orekit.propagation.SpacecraftState;
  44. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  45. import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
  46. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  47. import org.orekit.propagation.semianalytical.dsst.utilities.FieldCjSjCoefficient;
  48. import org.orekit.propagation.semianalytical.dsst.utilities.FieldShortPeriodicsInterpolatedCoefficient;
  49. import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
  50. import org.orekit.time.AbsoluteDate;
  51. import org.orekit.time.FieldAbsoluteDate;
  52. import org.orekit.utils.FieldTimeSpanMap;
  53. import org.orekit.utils.ParameterDriver;
  54. import org.orekit.utils.TimeSpanMap;

  56. import java.lang.reflect.Array;
  57. import java.util.ArrayList;
  58. import java.util.Collections;
  59. import java.util.HashMap;
  60. import java.util.List;
  61. import java.util.Map;
  62. import java.util.Set;

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

  102.     /**
  103.      * Retrograde factor I.
  104.      * <p>
  105.      * DSST model needs equinoctial orbit as internal representation. Classical
  106.      * equinoctial elements have discontinuities when inclination is close to zero.
  107.      * In this representation, I = +1. <br>
  108.      * To avoid this discontinuity, another representation exists and equinoctial
  109.      * elements can be expressed in a different way, called "retrograde" orbit. This
  110.      * implies I = -1. <br>
  111.      * As Orekit doesn't implement the retrograde orbit, I is always set to +1. But
  112.      * for the sake of consistency with the theory, the retrograde factor has been
  113.      * kept in the formulas.
  114.      * </p>
  115.      */
  116.     private static final int I = 1;

  118.     /**
  119.      * Central attraction scaling factor.
  120.      * <p>
  121.      * We use a power of 2 to avoid numeric noise introduction in the
  122.      * multiplications/divisions sequences.
  123.      * </p>
  124.      */
  125.     private static final double MU_SCALE = FastMath.scalb(1.0, 32);

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

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

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

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

  139.     /** Contribution to be numerically averaged. */
  140.     private final ForceModel contribution;

  142.     /** Gauss integrator. */
  143.     private final double threshold;

  145.     /** Gauss integrator. */
  146.     private GaussQuadrature integrator;

  148.     /** Flag for Gauss order computation. */
  149.     private boolean isDirty;

  151.     /** Attitude provider. */
  152.     private AttitudeProvider attitudeProvider;

  154.     /** Prefix for coefficients keys. */
  155.     private final String coefficientsKeyPrefix;

  157.     /** Short period terms. */
  158.     private GaussianShortPeriodicCoefficients gaussianSPCoefs;

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

  163.     /** Driver for gravitational parameter. */
  164.     private final ParameterDriver gmParameterDriver;

  166.     /**
  167.      * Build a new instance.
  168.      * @param coefficientsKeyPrefix prefix for coefficients keys
  169.      * @param threshold             tolerance for the choice of the Gauss quadrature
  170.      *                              order
  171.      * @param contribution          the {@link ForceModel} to be numerically
  172.      *                              averaged
  173.      * @param mu                    central attraction coefficient
  174.      */
  175.     protected AbstractGaussianContribution(final String coefficientsKeyPrefix, final double threshold,
  176.             final ForceModel contribution, final double mu) {

  178.         gmParameterDriver = new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT, mu, MU_SCALE,
  179.                 0.0, Double.POSITIVE_INFINITY);

  181.         this.coefficientsKeyPrefix = coefficientsKeyPrefix;
  182.         this.contribution = contribution;
  183.         this.threshold = threshold;
  184.         this.integrator = new GaussQuadrature(GAUSS_ORDER[MAX_ORDER_RANK]);
  185.         this.isDirty = true;

  187.         gaussianFieldSPCoefs = new HashMap<>();
  188.     }

  190.     /** {@inheritDoc} */
  191.     @Override
  192.     public void init(final SpacecraftState initialState, final AbsoluteDate target) {
  193.         // Initialize the numerical force model
  194.         contribution.init(initialState, target);
  195.     }

  197.     /** {@inheritDoc} */
  198.     @Override
  199.     public <T extends CalculusFieldElement<T>> void init(final FieldSpacecraftState<T> initialState, final FieldAbsoluteDate<T> target) {
  200.         // Initialize the numerical force model
  201.         contribution.init(initialState, target);
  202.     }

  204.     /** {@inheritDoc} */
  205.     @Override
  206.     public List<ParameterDriver> getParametersDrivers() {
  207.         // Initialize drivers (without central attraction coefficient driver)
  208.         final List<ParameterDriver> drivers = new ArrayList<>(getParametersDriversWithoutMu());
  209.         // We put central attraction coefficient driver at the end of the array
  210.         drivers.add(gmParameterDriver);
  211.         return drivers;
  212.     }

  214.     /**
  215.      * Get the drivers for force model parameters except the one for the central
  216.      * attraction coefficient.
  217.      * <p>
  218.      * The driver for central attraction coefficient is automatically added at the
  219.      * last element of the {@link ParameterDriver} array into
  220.      * {@link #getParametersDrivers()} method.
  221.      * </p>
  222.      * @return drivers for force model parameters
  223.      */
  224.     protected abstract List<ParameterDriver> getParametersDriversWithoutMu();

  226.     /** {@inheritDoc} */
  227.     @Override
  228.     public List<ShortPeriodTerms> initializeShortPeriodTerms(final AuxiliaryElements auxiliaryElements, final PropagationType type,
  229.             final double[] parameters) {

  231.         final List<ShortPeriodTerms> list = new ArrayList<>();
  232.         gaussianSPCoefs = new GaussianShortPeriodicCoefficients(coefficientsKeyPrefix, JMAX, INTERPOLATION_POINTS,
  233.                 new TimeSpanMap<>(new Slot(JMAX, INTERPOLATION_POINTS)));
  234.         list.add(gaussianSPCoefs);
  235.         return list;

  237.     }

  239.     /** {@inheritDoc} */
  240.     @Override
  241.     public <T extends CalculusFieldElement<T>> List<FieldShortPeriodTerms<T>> initializeShortPeriodTerms(
  242.             final FieldAuxiliaryElements<T> auxiliaryElements, final PropagationType type, final T[] parameters) {

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

  246.         final FieldGaussianShortPeriodicCoefficients<T> fgspc = new FieldGaussianShortPeriodicCoefficients<>(
  247.                 coefficientsKeyPrefix, JMAX, INTERPOLATION_POINTS,
  248.                 new FieldTimeSpanMap<>(new FieldSlot<>(JMAX, INTERPOLATION_POINTS), field));
  249.         gaussianFieldSPCoefs.put(field, fgspc);
  250.         return Collections.singletonList(fgspc);
  251.     }

  253.     /**
  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 auxiliaryElements auxiliary elements related to the current orbit
  259.      * @param parameters        parameters values of the force model parameters
  260.      *                          only 1 value for each parameterDriver
  261.      * @return new force model context
  262.      */
  263.     private AbstractGaussianContributionContext initializeStep(final AuxiliaryElements auxiliaryElements,
  264.             final double[] parameters) {
  265.         return new AbstractGaussianContributionContext(auxiliaryElements, parameters);
  266.     }

  268.     /**
  269.      * Performs initialization at each integration step for the current force model.
  270.      * <p>
  271.      * This method aims at being called before mean elements rates computation.
  272.      * </p>
  273.      * @param <T>               type of the elements
  274.      * @param auxiliaryElements auxiliary elements related to the current orbit
  275.      * @param parameters        parameters values of the force model parameters
  276.      *                          (only 1 values for each parameters corresponding
  277.      *                          to state date) by getting the parameters for a specific date.
  278.      * @return new force model context
  279.      */
  280.     private <T extends CalculusFieldElement<T>> FieldAbstractGaussianContributionContext<T> initializeStep(
  281.             final FieldAuxiliaryElements<T> auxiliaryElements, final T[] parameters) {
  282.         return new FieldAbstractGaussianContributionContext<>(auxiliaryElements, parameters);
  283.     }

  285.     /** {@inheritDoc} */
  286.     @Override
  287.     public double[] getMeanElementRate(final SpacecraftState state, final AuxiliaryElements auxiliaryElements,
  288.             final double[] parameters) {

  290.         // Container for attributes

  292.         final AbstractGaussianContributionContext context = initializeStep(auxiliaryElements, parameters);
  293.         double[] meanElementRate = new double[6];
  294.         // Computes the limits for the integral
  295.         final double[] ll = getLLimits(state, auxiliaryElements);
  296.         // Computes integrated mean element rates if Llow < Lhigh
  297.         if (ll[0] < ll[1]) {
  298.             meanElementRate = getMeanElementRate(state, integrator, ll[0], ll[1], context, parameters);
  299.             if (isDirty) {
  300.                 boolean next = true;
  301.                 for (int i = 0; i < MAX_ORDER_RANK && next; i++) {
  302.                     final double[] meanRates = getMeanElementRate(state, new GaussQuadrature(GAUSS_ORDER[i]), ll[0],
  303.                             ll[1], context, parameters);
  304.                     if (getRatesDiff(meanElementRate, meanRates, context) < threshold) {
  305.                         integrator = new GaussQuadrature(GAUSS_ORDER[i]);
  306.                         next = false;
  307.                     }
  308.                 }
  309.                 isDirty = false;
  310.             }
  311.         }
  312.         return meanElementRate;
  313.     }

  315.     /** {@inheritDoc} */
  316.     @Override
  317.     public <T extends CalculusFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> state,
  318.             final FieldAuxiliaryElements<T> auxiliaryElements, final T[] parameters) {

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

  324.         T[] meanElementRate = MathArrays.buildArray(field, 6);
  325.         // Computes the limits for the integral
  326.         final T[] ll = getLLimits(state, auxiliaryElements);
  327.         // Computes integrated mean element rates if Llow < Lhigh
  328.         if (ll[0].getReal() < ll[1].getReal()) {
  329.             meanElementRate = getMeanElementRate(state, integrator, ll[0], ll[1], context, parameters);
  330.             if (isDirty) {
  331.                 boolean next = true;
  332.                 for (int i = 0; i < MAX_ORDER_RANK && next; i++) {
  333.                     final T[] meanRates = getMeanElementRate(state, new GaussQuadrature(GAUSS_ORDER[i]), ll[0], ll[1],
  334.                             context, parameters);
  335.                     if (getRatesDiff(meanElementRate, meanRates, context).getReal() < threshold) {
  336.                         integrator = new GaussQuadrature(GAUSS_ORDER[i]);
  337.                         next = false;
  338.                     }
  339.                 }
  340.                 isDirty = false;
  341.             }
  342.         }

  344.         return meanElementRate;
  345.     }

  347.     /**
  348.      * Compute the limits in L, the true longitude, for integration.
  349.      *
  350.      * @param state             current state information: date, kinematics,
  351.      *                          attitude
  352.      * @param auxiliaryElements auxiliary elements related to the current orbit
  353.      * @return the integration limits in L
  354.      */
  355.     protected abstract double[] getLLimits(SpacecraftState state, AuxiliaryElements auxiliaryElements);

  357.     /**
  358.      * Compute the limits in L, the true longitude, for integration.
  359.      *
  360.      * @param <T>               type of the elements
  361.      * @param state             current state information: date, kinematics,
  362.      *                          attitude
  363.      * @param auxiliaryElements auxiliary elements related to the current orbit
  364.      * @return the integration limits in L
  365.      */
  366.     protected abstract <T extends CalculusFieldElement<T>> T[] getLLimits(FieldSpacecraftState<T> state,
  367.             FieldAuxiliaryElements<T> auxiliaryElements);

  369.     /**
  370.      * Computes the mean equinoctial elements rates da<sub>i</sub> / dt.
  371.      *
  372.      * @param state      current state
  373.      * @param gauss      Gauss quadrature
  374.      * @param low        lower bound of the integral interval
  375.      * @param high       upper bound of the integral interval
  376.      * @param context    container for attributes
  377.      * @param parameters values of the force model parameters
  378.      * at state date (1 values for each parameters)
  379.      * @return the mean element rates
  380.      */
  381.     protected double[] getMeanElementRate(final SpacecraftState state, final GaussQuadrature gauss, final double low,
  382.             final double high, final AbstractGaussianContributionContext context, final double[] parameters) {

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

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

  389.         // Constant multiplier for integral
  390.         final double coef = 1. / (2. * FastMath.PI * auxiliaryElements.getB());
  391.         // Corrects mean element rates
  392.         for (int i = 0; i < 6; i++) {
  393.             meanElementRate[i] *= coef;
  394.         }
  395.         return meanElementRate;
  396.     }

  398.     /**
  399.      * Computes the mean equinoctial elements rates da<sub>i</sub> / dt.
  400.      *
  401.      * @param <T>        type of the elements
  402.      * @param state      current state
  403.      * @param gauss      Gauss quadrature
  404.      * @param low        lower bound of the integral interval
  405.      * @param high       upper bound of the integral interval
  406.      * @param context    container for attributes
  407.      * @param parameters values of the force model parameters(1 values for each parameters)
  408.      * @return the mean element rates
  409.      */
  410.     protected <T extends CalculusFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> state,
  411.             final GaussQuadrature gauss, final T low, final T high,
  412.             final FieldAbstractGaussianContributionContext<T> context, final T[] parameters) {

  414.         // Field
  415.         final Field<T> field = context.getA().getField();

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

  420.         final T[] meanElementRate = gauss.integrate(new FieldIntegrableFunction<>(state, true, 0, parameters, field),
  421.                 low, high, field);
  422.         // Constant multiplier for integral
  423.         final T coef = auxiliaryElements.getB().multiply(low.getPi()).multiply(2.).reciprocal();
  424.         // Corrects mean element rates
  425.         for (int i = 0; i < 6; i++) {
  426.             meanElementRate[i] = meanElementRate[i].multiply(coef);
  427.         }
  428.         return meanElementRate;
  429.     }

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

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

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

  454.     /**
  455.      * Estimates the weighted magnitude of the difference between 2 sets of
  456.      * equinoctial elements rates.
  457.      *
  458.      * @param <T>     type of the elements
  459.      * @param meanRef reference rates
  460.      * @param meanCur current rates
  461.      * @param context container for attributes
  462.      * @return estimated magnitude of weighted differences
  463.      */
  464.     private <T extends CalculusFieldElement<T>> T getRatesDiff(final T[] meanRef, final T[] meanCur,
  465.             final FieldAbstractGaussianContributionContext<T> context) {

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

  470.         T maxDiff = FastMath.abs(meanRef[0].subtract(meanCur[0])).divide(auxiliaryElements.getSma());

  472.         // Corrects mean element rates
  473.         for (int i = 1; i < meanRef.length; i++) {
  474.             maxDiff = FastMath.max(maxDiff, FastMath.abs(meanRef[i].subtract(meanCur[i])));
  475.         }
  476.         return maxDiff;
  477.     }

  479.     /** {@inheritDoc} */
  480.     @Override
  481.     public void registerAttitudeProvider(final AttitudeProvider provider) {
  482.         this.attitudeProvider = provider;
  483.     }

  485.     /** {@inheritDoc} */
  486.     @Override
  487.     public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {

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

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

  495.             // Container of attributes
  496.             // Extract the proper parameters valid for the corresponding meanState date from the input array
  497.             final double[] extractedParameters = this.extractParameters(parameters, auxiliaryElements.getDate());
  498.             final AbstractGaussianContributionContext context = initializeStep(auxiliaryElements, extractedParameters);

  500.             // Compute rhoj and sigmaj
  501.             final double[][] currentRhoSigmaj = computeRhoSigmaCoefficients(auxiliaryElements);

  503.             // Generate the Cij and Sij coefficients
  504.             final FourierCjSjCoefficients fourierCjSj = new FourierCjSjCoefficients(meanState, JMAX, auxiliaryElements,
  505.                                                                                     extractedParameters);

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

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

  513.         }

  515.     }

  517.     /** {@inheritDoc} */
  518.     @Override
  519.     @SuppressWarnings("unchecked")
  520.     public <T extends CalculusFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
  521.             final FieldSpacecraftState<T>... meanStates) {

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

  526.         final FieldGaussianShortPeriodicCoefficients<T> fgspc = (FieldGaussianShortPeriodicCoefficients<T>) gaussianFieldSPCoefs
  527.                 .get(field);
  528.         final FieldSlot<T> slot = fgspc.createSlot(meanStates);
  529.         for (final FieldSpacecraftState<T> meanState : meanStates) {

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

  534.             // Container of attributes
  535.             // Extract the proper parameters valid for the corresponding meanState date from the input array
  536.             final T[] extractedParameters = this.extractParameters(parameters, auxiliaryElements.getDate());
  537.             final FieldAbstractGaussianContributionContext<T> context = initializeStep(auxiliaryElements, extractedParameters);

  539.             // Compute rhoj and sigmaj
  540.             final T[][] currentRhoSigmaj = computeRhoSigmaCoefficients(context, field);

  542.             // Generate the Cij and Sij coefficients
  543.             final FieldFourierCjSjCoefficients<T> fourierCjSj = new FieldFourierCjSjCoefficients<>(meanState, JMAX,
  544.                     auxiliaryElements, extractedParameters, field);

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

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

  553.         }

  555.     }

  557.     /**
  558.      * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
  559.      * <p>
  560.      * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
  561.      * ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
  562.      * σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
  563.      * </p>
  564.      * @param auxiliaryElements auxiliary elements related to the current orbit
  565.      * @return computed coefficients
  566.      */
  567.     private double[][] computeRhoSigmaCoefficients(final AuxiliaryElements auxiliaryElements) {
  568.         final double[][] currentRhoSigmaj = new double[2][3 * JMAX + 1];
  569.         final CjSjCoefficient cjsjKH = new CjSjCoefficient(auxiliaryElements.getK(), auxiliaryElements.getH());
  570.         final double b = 1. / (1 + auxiliaryElements.getB());

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

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

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

  586.     /**
  587.      * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
  588.      * <p>
  589.      * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
  590.      * ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
  591.      * σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
  592.      * </p>
  593.      * @param <T>     type of the elements
  594.      * @param context container for attributes
  595.      * @param field   field used by default
  596.      * @return computed coefficients
  597.      */
  598.     private <T extends CalculusFieldElement<T>> T[][] computeRhoSigmaCoefficients(final FieldAbstractGaussianContributionContext<T> context, final Field<T> field) {
  599.         // zero
  600.         final T zero = field.getZero();

  602.         final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
  603.         final T[][] currentRhoSigmaj = MathArrays.buildArray(field, 2, 3 * JMAX + 1);
  604.         final FieldCjSjCoefficient<T> cjsjKH = new FieldCjSjCoefficient<>(auxiliaryElements.getK(),
  605.                 auxiliaryElements.getH(), field);
  606.         final T b = auxiliaryElements.getB().add(1.).reciprocal();

  608.         // (-b)<sup>j</sup>
  609.         T mbtj = zero.newInstance(1.);

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

  613.             // Compute current rho and sigma;
  614.             mbtj = mbtj.multiply(b.negate());
  615.             final T coef = mbtj.multiply(auxiliaryElements.getB().multiply(j).add(1.));
  616.             currentRhoSigmaj[0][j] = coef.multiply(cjsjKH.getCj(j));
  617.             currentRhoSigmaj[1][j] = coef.multiply(cjsjKH.getSj(j));
  618.         }
  619.         return currentRhoSigmaj;
  620.     }

  622.     /**
  623.      * Internal class for numerical quadrature.
  624.      * <p>
  625.      * This class is a rewrite of {@link IntegrableFunction} for field elements
  626.      * </p>
  627.      * @param <T> type of the field elements
  628.      */
  629.     protected class FieldIntegrableFunction<T extends CalculusFieldElement<T>>
  630.             implements CalculusFieldUnivariateVectorFunction<T> {

  632.         /** Current state. */
  633.         private final FieldSpacecraftState<T> state;

  635.         /**
  636.          * Signal that this class is used to compute the values required by the mean
  637.          * element variations or by the short periodic element variations.
  638.          */
  639.         private final boolean meanMode;

  641.         /**
  642.          * The j index.
  643.          * <p>
  644.          * Used only for short periodic variation. Ignored for mean elements variation.
  645.          * </p>
  646.          */
  647.         private final int j;

  649.         /** Container for attributes. */
  650.         private final FieldAbstractGaussianContributionContext<T> context;

  652.         /** Auxiliary Elements. */
  653.         private final FieldAuxiliaryElements<T> auxiliaryElements;

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

  658.         /**
  659.          * Build a new instance with a new field.
  660.          * @param state      current state information: date, kinematics, attitude
  661.          * @param meanMode   if true return the value associated to the mean elements
  662.          *                   variation, if false return the values associated to the
  663.          *                   short periodic elements variation
  664.          * @param j          the j index. used only for short periodic variation.
  665.          *                   Ignored for mean elements variation.
  666.          * @param parameters values of the force model parameters (only 1 values
  667.          *                   for each parameters corresponding to state date) obtained by
  668.          *                   calling the extract parameter method {@link #extractParameters(double[], AbsoluteDate)}
  669.          *                   to selected the right value for state date or by getting the parameters for a specific date
  670.          * @param field      field utilized by default
  671.          */
  672.         public FieldIntegrableFunction(final FieldSpacecraftState<T> state, final boolean meanMode, final int j,
  673.                 final T[] parameters, final Field<T> field) {

  675.             this.meanMode = meanMode;
  676.             this.j = j;
  677.             this.parameters = parameters.clone();
  678.             this.auxiliaryElements = new FieldAuxiliaryElements<>(state.getOrbit(), I);
  679.             this.context = new FieldAbstractGaussianContributionContext<>(auxiliaryElements, this.parameters);
  680.             // remove derivatives from state
  681.             final T[] stateVector = MathArrays.buildArray(field, 6);
  682.             final PositionAngleType positionAngleType = PositionAngleType.MEAN;
  683.             OrbitType.EQUINOCTIAL.mapOrbitToArray(state.getOrbit(), positionAngleType, stateVector, null);
  684.             final FieldOrbit<T> fixedOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(stateVector, null,
  685.                     positionAngleType, state.getDate(), context.getMu(), state.getFrame());
  686.             this.state = new FieldSpacecraftState<>(fixedOrbit, state.getAttitude()).withMass(state.getMass());
  687.         }

  689.         /** {@inheritDoc} */
  690.         @Override
  691.         public T[] value(final T x) {

  693.             // Parameters for array building
  694.             final Field<T> field = auxiliaryElements.getDate().getField();
  695.             final int dimension = 6;

  697.             // Compute the time difference from the true longitude difference
  698.             final T shiftedLm = trueToMean(x);
  699.             final T dLm = shiftedLm.subtract(auxiliaryElements.getLM());
  700.             final T dt = dLm.divide(context.getMeanMotion());

  702.             final FieldSinCos<T> scL = FastMath.sinCos(x);
  703.             final T cosL = scL.cos();
  704.             final T sinL = scL.sin();
  705.             final T roa  = auxiliaryElements.getB().multiply(auxiliaryElements.getB()).divide(auxiliaryElements.getH().multiply(sinL).add(auxiliaryElements.getK().multiply(cosL)).add(1.));
  706.             final T roa2 = roa.multiply(roa);
  707.             final T r = auxiliaryElements.getSma().multiply(roa);
  708.             final T X = r.multiply(cosL);
  709.             final T Y = r.multiply(sinL);
  710.             final T naob = context.getMeanMotion().multiply(auxiliaryElements.getSma())
  711.                     .divide(auxiliaryElements.getB());
  712.             final T Xdot = naob.multiply(auxiliaryElements.getH().add(sinL)).negate();
  713.             final T Ydot = naob.multiply(auxiliaryElements.getK().add(cosL));
  714.             final FieldVector3D<T> vel = new FieldVector3D<>(Xdot, auxiliaryElements.getVectorF(), Ydot,
  715.                     auxiliaryElements.getVectorG());

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

  720.             // Recompose an orbit with time held fixed to be compliant with DSST theory
  721.             final FieldOrbit<T> recomposedOrbit = new FieldEquinoctialOrbit<>(shiftedOrbit.getA(),
  722.                     shiftedOrbit.getEquinoctialEx(), shiftedOrbit.getEquinoctialEy(), shiftedOrbit.getHx(),
  723.                     shiftedOrbit.getHy(), shiftedOrbit.getLM(), PositionAngleType.MEAN, shiftedOrbit.getFrame(),
  724.                     state.getDate(), context.getMu());

  726.             // Get the corresponding attitude
  727.             final FieldAttitude<T> recomposedAttitude;
  728.             if (contribution.dependsOnAttitudeRate()) {
  729.                 recomposedAttitude = attitudeProvider.getAttitude(recomposedOrbit,
  730.                         recomposedOrbit.getDate(), recomposedOrbit.getFrame());
  731.             } else {
  732.                 final FieldRotation<T> rotation = attitudeProvider.getAttitudeRotation(recomposedOrbit,
  733.                         recomposedOrbit.getDate(), recomposedOrbit.getFrame());
  734.                 final FieldVector3D<T> zeroVector = FieldVector3D.getZero(recomposedOrbit.getA().getField());
  735.                 recomposedAttitude = new FieldAttitude<>(recomposedOrbit.getDate(), recomposedOrbit.getFrame(),
  736.                         rotation, zeroVector, zeroVector);
  737.             }

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

  742.             final FieldVector3D<T> acc = contribution.acceleration(shiftedState, parameters);

  744.             // Compute the derivatives of the elements by the speed
  745.             final T[] deriv = MathArrays.buildArray(field, dimension);
  746.             // da/dv
  747.             deriv[0] = getAoV(vel).dotProduct(acc);
  748.             // dex/dv
  749.             deriv[1] = getKoV(X, Y, Xdot, Ydot).dotProduct(acc);
  750.             // dey/dv
  751.             deriv[2] = getHoV(X, Y, Xdot, Ydot).dotProduct(acc);
  752.             // dhx/dv
  753.             deriv[3] = getQoV(X).dotProduct(acc);
  754.             // dhy/dv
  755.             deriv[4] = getPoV(Y).dotProduct(acc);
  756.             // dλ/dv
  757.             deriv[5] = getLoV(X, Y, Xdot, Ydot).dotProduct(acc);

  759.             // Compute mean elements rates
  760.             final T[] val;
  761.             if (meanMode) {
  762.                 val = MathArrays.buildArray(field, dimension);
  763.                 for (int i = 0; i < 6; i++) {
  764.                     // da<sub>i</sub>/dt
  765.                     val[i] = deriv[i].multiply(roa2);
  766.                 }
  767.             } else {
  768.                 val = MathArrays.buildArray(field, dimension * 2);
  769.                 //Compute cos(j*L) and sin(j*L);
  770.                 final FieldSinCos<T> scjL = FastMath.sinCos(x.multiply(j));
  771.                 final T cosjL = j == 1 ? cosL : scjL.cos();
  772.                 final T sinjL = j == 1 ? sinL : scjL.sin();

  774.                 for (int i = 0; i < 6; i++) {
  775.                     // da<sub>i</sub>/dv * cos(jL)
  776.                     val[i] = deriv[i].multiply(cosjL);
  777.                     // da<sub>i</sub>/dv * sin(jL)
  778.                     val[i + 6] = deriv[i].multiply(sinjL);
  779.                 }
  780.             }

  782.             return val;
  783.         }

  785.         /**
  786.          * Converts true longitude to mean longitude.
  787.          * @param x True longitude
  788.          * @return Eccentric longitude
  789.          */
  790.         private T trueToMean(final T x) {
  791.             return eccentricToMean(trueToEccentric(x));
  792.         }

  794.         /**
  795.          * Converts true longitude to eccentric longitude.
  796.          * @param lv True longitude
  797.          * @return Eccentric longitude
  798.          */
  799.         private T trueToEccentric (final T lv) {
  800.             final FieldSinCos<T> sclV = FastMath.sinCos(lv);
  801.             final T cosLv   = sclV.cos();
  802.             final T sinLv   = sclV.sin();
  803.             final T num     = auxiliaryElements.getH().multiply(cosLv).subtract(auxiliaryElements.getK().multiply(sinLv));
  804.             final T den     = auxiliaryElements.getB().add(auxiliaryElements.getK().multiply(cosLv)).add(auxiliaryElements.getH().multiply(sinLv)).add(1.);
  805.             return FastMath.atan(num.divide(den)).multiply(2.).add(lv);
  806.         }

  808.         /**
  809.          * Converts eccentric longitude to mean longitude.
  810.          * @param le Eccentric longitude
  811.          * @return Mean longitude
  812.          */
  813.         private T eccentricToMean (final T le) {
  814.             final FieldSinCos<T> scle = FastMath.sinCos(le);
  815.             return le.subtract(auxiliaryElements.getK().multiply(scle.sin())).add(auxiliaryElements.getH().multiply(scle.cos()));
  816.         }

  818.         /**
  819.          * Compute δa/δv.
  820.          * @param vel satellite velocity
  821.          * @return δa/δv
  822.          */
  823.         private FieldVector3D<T> getAoV(final FieldVector3D<T> vel) {
  824.             return new FieldVector3D<>(context.getTon2a(), vel);
  825.         }

  827.         /**
  828.          * Compute δh/δv.
  829.          * @param X    satellite position component along f, equinoctial reference frame
  830.          *             1st vector
  831.          * @param Y    satellite position component along g, equinoctial reference frame
  832.          *             2nd vector
  833.          * @param Xdot satellite velocity component along f, equinoctial reference frame
  834.          *             1st vector
  835.          * @param Ydot satellite velocity component along g, equinoctial reference frame
  836.          *             2nd vector
  837.          * @return δh/δv
  838.          */
  839.         private FieldVector3D<T> getHoV(final T X, final T Y, final T Xdot, final T Ydot) {
  840.             final T kf = (Xdot.multiply(Y).multiply(2.).subtract(X.multiply(Ydot))).multiply(context.getOoMU());
  841.             final T kg = X.multiply(Xdot).multiply(context.getOoMU());
  842.             final T kw = auxiliaryElements.getK().multiply(
  843.                     auxiliaryElements.getQ().multiply(Y).multiply(I).subtract(auxiliaryElements.getP().multiply(X)))
  844.                     .multiply(context.getOOAB());
  845.             return new FieldVector3D<>(kf, auxiliaryElements.getVectorF(), kg.negate(), auxiliaryElements.getVectorG(),
  846.                     kw, auxiliaryElements.getVectorW());
  847.         }

  849.         /**
  850.          * Compute δk/δv.
  851.          * @param X    satellite position component along f, equinoctial reference frame
  852.          *             1st vector
  853.          * @param Y    satellite position component along g, equinoctial reference frame
  854.          *             2nd vector
  855.          * @param Xdot satellite velocity component along f, equinoctial reference frame
  856.          *             1st vector
  857.          * @param Ydot satellite velocity component along g, equinoctial reference frame
  858.          *             2nd vector
  859.          * @return δk/δv
  860.          */
  861.         private FieldVector3D<T> getKoV(final T X, final T Y, final T Xdot, final T Ydot) {
  862.             final T kf = Y.multiply(Ydot).multiply(context.getOoMU());
  863.             final T kg = (X.multiply(Ydot).multiply(2.).subtract(Xdot.multiply(Y))).multiply(context.getOoMU());
  864.             final T kw = auxiliaryElements.getH().multiply(
  865.                     auxiliaryElements.getQ().multiply(Y).multiply(I).subtract(auxiliaryElements.getP().multiply(X)))
  866.                     .multiply(context.getOOAB());
  867.             return new FieldVector3D<>(kf.negate(), auxiliaryElements.getVectorF(), kg, auxiliaryElements.getVectorG(),
  868.                     kw.negate(), auxiliaryElements.getVectorW());
  869.         }

  871.         /**
  872.          * Compute δp/δv.
  873.          * @param Y satellite position component along g, equinoctial reference frame
  874.          *          2nd vector
  875.          * @return δp/δv
  876.          */
  877.         private FieldVector3D<T> getPoV(final T Y) {
  878.             return new FieldVector3D<>(context.getCo2AB().multiply(Y), auxiliaryElements.getVectorW());
  879.         }

  881.         /**
  882.          * Compute δq/δv.
  883.          * @param X satellite position component along f, equinoctial reference frame
  884.          *          1st vector
  885.          * @return δq/δv
  886.          */
  887.         private FieldVector3D<T> getQoV(final T X) {
  888.             return new FieldVector3D<>(context.getCo2AB().multiply(X).multiply(I), auxiliaryElements.getVectorW());
  889.         }

  891.         /**
  892.          * Compute δλ/δv.
  893.          * @param X    satellite position component along f, equinoctial reference frame
  894.          *             1st vector
  895.          * @param Y    satellite position component along g, equinoctial reference frame
  896.          *             2nd vector
  897.          * @param Xdot satellite velocity component along f, equinoctial reference frame
  898.          *             1st vector
  899.          * @param Ydot satellite velocity component along g, equinoctial reference frame
  900.          *             2nd vector
  901.          * @return δλ/δv
  902.          */
  903.         private FieldVector3D<T> getLoV(final T X, final T Y, final T Xdot, final T Ydot) {
  904.             final FieldVector3D<T> pos = new FieldVector3D<>(X, auxiliaryElements.getVectorF(), Y,
  905.                     auxiliaryElements.getVectorG());
  906.             final FieldVector3D<T> v2 = new FieldVector3D<>(auxiliaryElements.getK(), getHoV(X, Y, Xdot, Ydot),
  907.                     auxiliaryElements.getH().negate(), getKoV(X, Y, Xdot, Ydot));
  908.             return new FieldVector3D<>(context.getOOA().multiply(-2.), pos, context.getOoBpo(), v2,
  909.                     context.getOOA().multiply(auxiliaryElements.getQ().multiply(Y).multiply(I)
  910.                             .subtract(auxiliaryElements.getP().multiply(X))),
  911.                     auxiliaryElements.getVectorW());
  912.         }

  914.     }

  916.     /** Internal class for numerical quadrature. */
  917.     protected class IntegrableFunction implements UnivariateVectorFunction {

  919.         /** Current state. */
  920.         private final SpacecraftState state;

  922.         /**
  923.          * Signal that this class is used to compute the values required by the mean
  924.          * element variations or by the short periodic element variations.
  925.          */
  926.         private final boolean meanMode;

  928.         /**
  929.          * The j index.
  930.          * <p>
  931.          * Used only for short periodic variation. Ignored for mean elements variation.
  932.          * </p>
  933.          */
  934.         private final int j;

  936.         /** Container for attributes. */
  937.         private final AbstractGaussianContributionContext context;

  939.         /** Auxiliary Elements. */
  940.         private final AuxiliaryElements auxiliaryElements;

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

  945.         /**
  946.          * Build a new instance.
  947.          * @param state      current state information: date, kinematics, attitude
  948.          * @param meanMode   if true return the value associated to the mean elements
  949.          *                   variation, if false return the values associated to the
  950.          *                   short periodic elements variation
  951.          * @param j          the j index. used only for short periodic variation.
  952.          *                   Ignored for mean elements variation.
  953.          * @param parameters list of the estimated values for each driver at state date of the force model parameters
  954.          *                   only 1 value for each parameter
  955.          */
  956.         IntegrableFunction(final SpacecraftState state, final boolean meanMode, final int j,
  957.                 final double[] parameters) {

  959.             this.meanMode = meanMode;
  960.             this.j = j;
  961.             this.parameters = parameters.clone();
  962.             this.auxiliaryElements = new AuxiliaryElements(state.getOrbit(), I);
  963.             this.context = new AbstractGaussianContributionContext(auxiliaryElements, this.parameters);
  964.             // remove derivatives from state
  965.             final double[] stateVector = new double[6];
  966.             final PositionAngleType positionAngleType = PositionAngleType.MEAN;
  967.             OrbitType.EQUINOCTIAL.mapOrbitToArray(state.getOrbit(), positionAngleType, stateVector, null);
  968.             final Orbit fixedOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(stateVector, null, positionAngleType,
  969.                     state.getDate(), context.getMu(), state.getFrame());
  970.             this.state = new SpacecraftState(fixedOrbit, state.getAttitude()).withMass(state.getMass());
  971.         }

  973.         /** {@inheritDoc} */
  974.         @SuppressWarnings("checkstyle:FinalLocalVariable")
  975.         @Override
  976.         public double[] value(final double x) {

  978.             // Compute the time difference from the true longitude difference
  979.             final double shiftedLm = trueToMean(x);
  980.             final double dLm = shiftedLm - auxiliaryElements.getLM();
  981.             final double dt = dLm / context.getMeanMotion();

  983.             final SinCos scL  = FastMath.sinCos(x);
  984.             final double cosL = scL.cos();
  985.             final double sinL = scL.sin();
  986.             final double roa  = auxiliaryElements.getB() * auxiliaryElements.getB() / (1. + auxiliaryElements.getH() * sinL + auxiliaryElements.getK() * cosL);
  987.             final double roa2 = roa * roa;
  988.             final double r = auxiliaryElements.getSma() * roa;
  989.             final double X = r * cosL;
  990.             final double Y = r * sinL;
  991.             final double naob = context.getMeanMotion() * auxiliaryElements.getSma() / auxiliaryElements.getB();
  992.             final double Xdot = -naob * (auxiliaryElements.getH() + sinL);
  993.             final double Ydot = naob * (auxiliaryElements.getK() + cosL);
  994.             final Vector3D vel = new Vector3D(Xdot, auxiliaryElements.getVectorF(), Ydot,
  995.                     auxiliaryElements.getVectorG());

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

  1000.             // Recompose an orbit with time held fixed to be compliant with DSST theory
  1001.             final Orbit recomposedOrbit = new EquinoctialOrbit(shiftedOrbit.getA(), shiftedOrbit.getEquinoctialEx(),
  1002.                     shiftedOrbit.getEquinoctialEy(), shiftedOrbit.getHx(), shiftedOrbit.getHy(), shiftedOrbit.getLM(),
  1003.                     PositionAngleType.MEAN, shiftedOrbit.getFrame(), state.getDate(), context.getMu());

  1005.             // Get the corresponding attitude
  1006.             final Attitude recomposedAttitude;
  1007.             if (contribution.dependsOnAttitudeRate()) {
  1008.                 recomposedAttitude = attitudeProvider.getAttitude(recomposedOrbit,
  1009.                         recomposedOrbit.getDate(), recomposedOrbit.getFrame());
  1010.             } else {
  1011.                 final Rotation rotation = attitudeProvider.getAttitudeRotation(recomposedOrbit,
  1012.                         recomposedOrbit.getDate(), recomposedOrbit.getFrame());
  1013.                 final Vector3D zeroVector = Vector3D.ZERO;
  1014.                 recomposedAttitude = new Attitude(recomposedOrbit.getDate(), recomposedOrbit.getFrame(),
  1015.                         rotation, zeroVector, zeroVector);
  1016.             }

  1018.             // create shifted SpacecraftState with attitude at specified time
  1019.             final SpacecraftState shiftedState = new SpacecraftState(recomposedOrbit, recomposedAttitude).withMass(state.getMass());

  1021.             // here parameters is a list of all span values of each parameter driver
  1022.             final Vector3D acc = contribution.acceleration(shiftedState, parameters);

  1024.             // Compute the derivatives of the elements by the speed
  1025.             final double[] deriv = new double[6];
  1026.             // da/dv
  1027.             deriv[0] = getAoV(vel).dotProduct(acc);
  1028.             // dex/dv
  1029.             deriv[1] = getKoV(X, Y, Xdot, Ydot).dotProduct(acc);
  1030.             // dey/dv
  1031.             deriv[2] = getHoV(X, Y, Xdot, Ydot).dotProduct(acc);
  1032.             // dhx/dv
  1033.             deriv[3] = getQoV(X).dotProduct(acc);
  1034.             // dhy/dv
  1035.             deriv[4] = getPoV(Y).dotProduct(acc);
  1036.             // dλ/dv
  1037.             deriv[5] = getLoV(X, Y, Xdot, Ydot).dotProduct(acc);

  1039.             // Compute mean elements rates
  1040.             final double[] val;
  1041.             if (meanMode) {
  1042.                 val = new double[6];
  1043.                 for (int i = 0; i < 6; i++) {
  1044.                     // da<sub>i</sub>/dt
  1045.                     val[i] = roa2 * deriv[i];
  1046.                 }
  1047.             } else {
  1048.                 val = new double[12];
  1049.                 //Compute cos(j*L) and sin(j*L);
  1050.                 final SinCos scjL  = FastMath.sinCos(j * x);
  1051.                 final double cosjL = j == 1 ? cosL : scjL.cos();
  1052.                 final double sinjL = j == 1 ? sinL : scjL.sin();

  1054.                 for (int i = 0; i < 6; i++) {
  1055.                     // da<sub>i</sub>/dv * cos(jL)
  1056.                     val[i] = cosjL * deriv[i];
  1057.                     // da<sub>i</sub>/dv * sin(jL)
  1058.                     val[i + 6] = sinjL * deriv[i];
  1059.                 }
  1060.             }
  1061.             return val;
  1062.         }

  1064.         /**
  1065.          * Converts true longitude to eccentric longitude.
  1066.          * @param lv True longitude
  1067.          * @return Eccentric longitude
  1068.          */
  1069.         private double trueToEccentric (final double lv) {
  1070.             final SinCos scLv    = FastMath.sinCos(lv);
  1071.             final double num     = auxiliaryElements.getH() * scLv.cos() - auxiliaryElements.getK() * scLv.sin();
  1072.             final double den     = auxiliaryElements.getB() + 1. + auxiliaryElements.getK() * scLv.cos() + auxiliaryElements.getH() * scLv.sin();
  1073.             return lv + 2. * FastMath.atan(num / den);
  1074.         }

  1076.         /**
  1077.          * Converts eccentric longitude to mean longitude.
  1078.          * @param le Eccentric longitude
  1079.          * @return Mean longitude
  1080.          */
  1081.         private double eccentricToMean (final double le) {
  1082.             final SinCos scLe = FastMath.sinCos(le);
  1083.             return le - auxiliaryElements.getK() * scLe.sin() + auxiliaryElements.getH() * scLe.cos();
  1084.         }

  1086.         /**
  1087.          * Converts true longitude to mean longitude.
  1088.          * @param lv True longitude
  1089.          * @return Eccentric longitude
  1090.          */
  1091.         private double trueToMean(final double lv) {
  1092.             return eccentricToMean(trueToEccentric(lv));
  1093.         }

  1095.         /**
  1096.          * Compute δa/δv.
  1097.          * @param vel satellite velocity
  1098.          * @return δa/δv
  1099.          */
  1100.         private Vector3D getAoV(final Vector3D vel) {
  1101.             return new Vector3D(context.getTon2a(), vel);
  1102.         }

  1104.         /**
  1105.          * Compute δh/δv.
  1106.          * @param X    satellite position component along f, equinoctial reference frame
  1107.          *             1st vector
  1108.          * @param Y    satellite position component along g, equinoctial reference frame
  1109.          *             2nd vector
  1110.          * @param Xdot satellite velocity component along f, equinoctial reference frame
  1111.          *             1st vector
  1112.          * @param Ydot satellite velocity component along g, equinoctial reference frame
  1113.          *             2nd vector
  1114.          * @return δh/δv
  1115.          */
  1116.         private Vector3D getHoV(final double X, final double Y, final double Xdot, final double Ydot) {
  1117.             final double kf = (2. * Xdot * Y - X * Ydot) * context.getOoMU();
  1118.             final double kg = X * Xdot * context.getOoMU();
  1119.             final double kw = auxiliaryElements.getK() *
  1120.                     (I * auxiliaryElements.getQ() * Y - auxiliaryElements.getP() * X) * context.getOOAB();
  1121.             return new Vector3D(kf, auxiliaryElements.getVectorF(), -kg, auxiliaryElements.getVectorG(), kw,
  1122.                     auxiliaryElements.getVectorW());
  1123.         }

  1125.         /**
  1126.          * Compute δk/δv.
  1127.          * @param X    satellite position component along f, equinoctial reference frame
  1128.          *             1st vector
  1129.          * @param Y    satellite position component along g, equinoctial reference frame
  1130.          *             2nd vector
  1131.          * @param Xdot satellite velocity component along f, equinoctial reference frame
  1132.          *             1st vector
  1133.          * @param Ydot satellite velocity component along g, equinoctial reference frame
  1134.          *             2nd vector
  1135.          * @return δk/δv
  1136.          */
  1137.         private Vector3D getKoV(final double X, final double Y, final double Xdot, final double Ydot) {
  1138.             final double kf = Y * Ydot * context.getOoMU();
  1139.             final double kg = (2. * X * Ydot - Xdot * Y) * context.getOoMU();
  1140.             final double kw = auxiliaryElements.getH() *
  1141.                     (I * auxiliaryElements.getQ() * Y - auxiliaryElements.getP() * X) * context.getOOAB();
  1142.             return new Vector3D(-kf, auxiliaryElements.getVectorF(), kg, auxiliaryElements.getVectorG(), -kw,
  1143.                     auxiliaryElements.getVectorW());
  1144.         }

  1146.         /**
  1147.          * Compute δp/δv.
  1148.          * @param Y satellite position component along g, equinoctial reference frame
  1149.          *          2nd vector
  1150.          * @return δp/δv
  1151.          */
  1152.         private Vector3D getPoV(final double Y) {
  1153.             return new Vector3D(context.getCo2AB() * Y, auxiliaryElements.getVectorW());
  1154.         }

  1156.         /**
  1157.          * Compute δq/δv.
  1158.          * @param X satellite position component along f, equinoctial reference frame
  1159.          *          1st vector
  1160.          * @return δq/δv
  1161.          */
  1162.         private Vector3D getQoV(final double X) {
  1163.             return new Vector3D(I * context.getCo2AB() * X, auxiliaryElements.getVectorW());
  1164.         }

  1166.         /**
  1167.          * Compute δλ/δv.
  1168.          * @param X    satellite position component along f, equinoctial reference frame
  1169.          *             1st vector
  1170.          * @param Y    satellite position component along g, equinoctial reference frame
  1171.          *             2nd vector
  1172.          * @param Xdot satellite velocity component along f, equinoctial reference frame
  1173.          *             1st vector
  1174.          * @param Ydot satellite velocity component along g, equinoctial reference frame
  1175.          *             2nd vector
  1176.          * @return δλ/δv
  1177.          */
  1178.         private Vector3D getLoV(final double X, final double Y, final double Xdot, final double Ydot) {
  1179.             final Vector3D pos = new Vector3D(X, auxiliaryElements.getVectorF(), Y, auxiliaryElements.getVectorG());
  1180.             final Vector3D v2 = new Vector3D(auxiliaryElements.getK(), getHoV(X, Y, Xdot, Ydot),
  1181.                     -auxiliaryElements.getH(), getKoV(X, Y, Xdot, Ydot));
  1182.             return new Vector3D(-2. * context.getOOA(), pos, context.getOoBpo(), v2,
  1183.                     (I * auxiliaryElements.getQ() * Y - auxiliaryElements.getP() * X) * context.getOOA(),
  1184.                     auxiliaryElements.getVectorW());
  1185.         }

  1187.     }

  1189.     /**
  1190.      * Class used to {@link #integrate(UnivariateVectorFunction, double, double)
  1191.      * integrate} a {@link org.hipparchus.analysis.UnivariateVectorFunction
  1192.      * function} of the orbital elements using the Gaussian quadrature rule to get
  1193.      * the acceleration.
  1194.      */
  1195.     protected static class GaussQuadrature {

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

  1199.         /** Points for quadrature of order 12. */
  1200.         private static final double[] P_12 = { -0.98156063424671910000, -0.90411725637047490000,
  1201.             -0.76990267419430470000, -0.58731795428661740000, -0.36783149899818024000, -0.12523340851146890000,
  1202.             0.12523340851146890000, 0.36783149899818024000, 0.58731795428661740000, 0.76990267419430470000,
  1203.             0.90411725637047490000, 0.98156063424671910000 };

  1205.         /** Weights for quadrature of order 12. */
  1206.         private static final double[] W_12 = { 0.04717533638651220000, 0.10693932599531830000, 0.16007832854334633000,
  1207.             0.20316742672306584000, 0.23349253653835478000, 0.24914704581340286000, 0.24914704581340286000,
  1208.             0.23349253653835478000, 0.20316742672306584000, 0.16007832854334633000, 0.10693932599531830000,
  1209.             0.04717533638651220000 };

  1211.         /** Points for quadrature of order 16. */
  1212.         private static final double[] P_16 = { -0.98940093499164990000, -0.94457502307323260000,
  1213.             -0.86563120238783160000, -0.75540440835500310000, -0.61787624440264380000, -0.45801677765722737000,
  1214.             -0.28160355077925890000, -0.09501250983763745000, 0.09501250983763745000, 0.28160355077925890000,
  1215.             0.45801677765722737000, 0.61787624440264380000, 0.75540440835500310000, 0.86563120238783160000,
  1216.             0.94457502307323260000, 0.98940093499164990000 };

  1218.         /** Weights for quadrature of order 16. */
  1219.         private static final double[] W_16 = { 0.02715245941175405800, 0.06225352393864777000, 0.09515851168249283000,
  1220.             0.12462897125553388000, 0.14959598881657685000, 0.16915651939500256000, 0.18260341504492360000,
  1221.             0.18945061045506847000, 0.18945061045506847000, 0.18260341504492360000, 0.16915651939500256000,
  1222.             0.14959598881657685000, 0.12462897125553388000, 0.09515851168249283000, 0.06225352393864777000,
  1223.             0.02715245941175405800 };

  1225.         /** Points for quadrature of order 20. */
  1226.         private static final double[] P_20 = { -0.99312859918509490000, -0.96397192727791390000,
  1227.             -0.91223442825132600000, -0.83911697182221890000, -0.74633190646015080000, -0.63605368072651510000,
  1228.             -0.51086700195082700000, -0.37370608871541955000, -0.22778585114164507000, -0.07652652113349734000,
  1229.             0.07652652113349734000, 0.22778585114164507000, 0.37370608871541955000, 0.51086700195082700000,
  1230.             0.63605368072651510000, 0.74633190646015080000, 0.83911697182221890000, 0.91223442825132600000,
  1231.             0.96397192727791390000, 0.99312859918509490000 };

  1233.         /** Weights for quadrature of order 20. */
  1234.         private static final double[] W_20 = { 0.01761400713915226400, 0.04060142980038684000, 0.06267204833410904000,
  1235.             0.08327674157670477000, 0.10193011981724048000, 0.11819453196151844000, 0.13168863844917678000,
  1236.             0.14209610931838212000, 0.14917298647260380000, 0.15275338713072600000, 0.15275338713072600000,
  1237.             0.14917298647260380000, 0.14209610931838212000, 0.13168863844917678000, 0.11819453196151844000,
  1238.             0.10193011981724048000, 0.08327674157670477000, 0.06267204833410904000, 0.04060142980038684000,
  1239.             0.01761400713915226400 };

  1241.         /** Points for quadrature of order 24. */
  1242.         private static final double[] P_24 = { -0.99518721999702130000, -0.97472855597130950000,
  1243.             -0.93827455200273270000, -0.88641552700440100000, -0.82000198597390300000, -0.74012419157855440000,
  1244.             -0.64809365193697550000, -0.54542147138883950000, -0.43379350762604520000, -0.31504267969616340000,
  1245.             -0.19111886747361634000, -0.06405689286260563000, 0.06405689286260563000, 0.19111886747361634000,
  1246.             0.31504267969616340000, 0.43379350762604520000, 0.54542147138883950000, 0.64809365193697550000,
  1247.             0.74012419157855440000, 0.82000198597390300000, 0.88641552700440100000, 0.93827455200273270000,
  1248.             0.97472855597130950000, 0.99518721999702130000 };

  1250.         /** Weights for quadrature of order 24. */
  1251.         private static final double[] W_24 = { 0.01234122979998733500, 0.02853138862893380600, 0.04427743881741981000,
  1252.             0.05929858491543691500, 0.07334648141108027000, 0.08619016153195320000, 0.09761865210411391000,
  1253.             0.10744427011596558000, 0.11550566805372553000, 0.12167047292780335000, 0.12583745634682825000,
  1254.             0.12793819534675221000, 0.12793819534675221000, 0.12583745634682825000, 0.12167047292780335000,
  1255.             0.11550566805372553000, 0.10744427011596558000, 0.09761865210411391000, 0.08619016153195320000,
  1256.             0.07334648141108027000, 0.05929858491543691500, 0.04427743881741981000, 0.02853138862893380600,
  1257.             0.01234122979998733500 };

  1259.         /** Points for quadrature of order 32. */
  1260.         private static final double[] P_32 = { -0.99726386184948160000, -0.98561151154526840000,
  1261.             -0.96476225558750640000, -0.93490607593773970000, -0.89632115576605220000, -0.84936761373256990000,
  1262.             -0.79448379596794250000, -0.73218211874028970000, -0.66304426693021520000, -0.58771575724076230000,
  1263.             -0.50689990893222950000, -0.42135127613063540000, -0.33186860228212767000, -0.23928736225213710000,
  1264.             -0.14447196158279646000, -0.04830766568773831000, 0.04830766568773831000, 0.14447196158279646000,
  1265.             0.23928736225213710000, 0.33186860228212767000, 0.42135127613063540000, 0.50689990893222950000,
  1266.             0.58771575724076230000, 0.66304426693021520000, 0.73218211874028970000, 0.79448379596794250000,
  1267.             0.84936761373256990000, 0.89632115576605220000, 0.93490607593773970000, 0.96476225558750640000,
  1268.             0.98561151154526840000, 0.99726386184948160000 };

  1270.         /** Weights for quadrature of order 32. */
  1271.         private static final double[] W_32 = { 0.00701861000947013600, 0.01627439473090571200, 0.02539206530926214200,
  1272.             0.03427386291302141000, 0.04283589802222658600, 0.05099805926237621600, 0.05868409347853559000,
  1273.             0.06582222277636193000, 0.07234579410884862000, 0.07819389578707042000, 0.08331192422694673000,
  1274.             0.08765209300440380000, 0.09117387869576390000, 0.09384439908080441000, 0.09563872007927487000,
  1275.             0.09654008851472784000, 0.09654008851472784000, 0.09563872007927487000, 0.09384439908080441000,
  1276.             0.09117387869576390000, 0.08765209300440380000, 0.08331192422694673000, 0.07819389578707042000,
  1277.             0.07234579410884862000, 0.06582222277636193000, 0.05868409347853559000, 0.05099805926237621600,
  1278.             0.04283589802222658600, 0.03427386291302141000, 0.02539206530926214200, 0.01627439473090571200,
  1279.             0.00701861000947013600 };

  1281.         /** Points for quadrature of order 40. */
  1282.         private static final double[] P_40 = { -0.99823770971055930000, -0.99072623869945710000,
  1283.             -0.97725994998377420000, -0.95791681921379170000, -0.93281280827867660000, -0.90209880696887420000,
  1284.             -0.86595950321225960000, -0.82461223083331170000, -0.77830565142651940000, -0.72731825518992710000,
  1285.             -0.67195668461417960000, -0.61255388966798030000, -0.54946712509512820000, -0.48307580168617870000,
  1286.             -0.41377920437160500000, -0.34199409082575850000, -0.26815218500725370000, -0.19269758070137110000,
  1287.             -0.11608407067525522000, -0.03877241750605081600, 0.03877241750605081600, 0.11608407067525522000,
  1288.             0.19269758070137110000, 0.26815218500725370000, 0.34199409082575850000, 0.41377920437160500000,
  1289.             0.48307580168617870000, 0.54946712509512820000, 0.61255388966798030000, 0.67195668461417960000,
  1290.             0.72731825518992710000, 0.77830565142651940000, 0.82461223083331170000, 0.86595950321225960000,
  1291.             0.90209880696887420000, 0.93281280827867660000, 0.95791681921379170000, 0.97725994998377420000,
  1292.             0.99072623869945710000, 0.99823770971055930000 };

  1294.         /** Weights for quadrature of order 40. */
  1295.         private static final double[] W_40 = { 0.00452127709853309800, 0.01049828453115270400, 0.01642105838190797300,
  1296.             0.02224584919416689000, 0.02793700698002338000, 0.03346019528254786500, 0.03878216797447199000,
  1297.             0.04387090818567333000, 0.04869580763507221000, 0.05322784698393679000, 0.05743976909939157000,
  1298.             0.06130624249292891000, 0.06480401345660108000, 0.06791204581523394000, 0.07061164739128681000,
  1299.             0.07288658239580408000, 0.07472316905796833000, 0.07611036190062619000, 0.07703981816424793000,
  1300.             0.07750594797842482000, 0.07750594797842482000, 0.07703981816424793000, 0.07611036190062619000,
  1301.             0.07472316905796833000, 0.07288658239580408000, 0.07061164739128681000, 0.06791204581523394000,
  1302.             0.06480401345660108000, 0.06130624249292891000, 0.05743976909939157000, 0.05322784698393679000,
  1303.             0.04869580763507221000, 0.04387090818567333000, 0.03878216797447199000, 0.03346019528254786500,
  1304.             0.02793700698002338000, 0.02224584919416689000, 0.01642105838190797300, 0.01049828453115270400,
  1305.             0.00452127709853309800 };

  1307.         /** Points for quadrature of order 48. */
  1308.         private static final double[] P_48 = { -0.99877100725242610000, -0.99353017226635080000,
  1309.             -0.98412458372282700000, -0.97059159254624720000, -0.95298770316043080000, -0.93138669070655440000,
  1310.             -0.90587913671556960000, -0.87657202027424800000, -0.84358826162439350000, -0.80706620402944250000,
  1311.             -0.76715903251574020000, -0.72403413092381470000, -0.67787237963266400000, -0.62886739677651370000,
  1312.             -0.57722472608397270000, -0.52316097472223300000, -0.46690290475095840000, -0.40868648199071680000,
  1313.             -0.34875588629216070000, -0.28736248735545555000, -0.22476379039468908000, -0.16122235606889174000,
  1314.             -0.09700469920946270000, -0.03238017096286937000, 0.03238017096286937000, 0.09700469920946270000,
  1315.             0.16122235606889174000, 0.22476379039468908000, 0.28736248735545555000, 0.34875588629216070000,
  1316.             0.40868648199071680000, 0.46690290475095840000, 0.52316097472223300000, 0.57722472608397270000,
  1317.             0.62886739677651370000, 0.67787237963266400000, 0.72403413092381470000, 0.76715903251574020000,
  1318.             0.80706620402944250000, 0.84358826162439350000, 0.87657202027424800000, 0.90587913671556960000,
  1319.             0.93138669070655440000, 0.95298770316043080000, 0.97059159254624720000, 0.98412458372282700000,
  1320.             0.99353017226635080000, 0.99877100725242610000 };

  1322.         /** Weights for quadrature of order 48. */
  1323.         private static final double[] W_48 = { 0.00315334605230596250, 0.00732755390127620800, 0.01147723457923446900,
  1324.             0.01557931572294386600, 0.01961616045735556700, 0.02357076083932435600, 0.02742650970835688000,
  1325.             0.03116722783279807000, 0.03477722256477045000, 0.03824135106583080600, 0.04154508294346483000,
  1326.             0.04467456085669424000, 0.04761665849249054000, 0.05035903555385448000, 0.05289018948519365000,
  1327.             0.05519950369998416500, 0.05727729210040315000, 0.05911483969839566000, 0.06070443916589384000,
  1328.             0.06203942315989268000, 0.06311419228625403000, 0.06392423858464817000, 0.06446616443595010000,
  1329.             0.06473769681268386000, 0.06473769681268386000, 0.06446616443595010000, 0.06392423858464817000,
  1330.             0.06311419228625403000, 0.06203942315989268000, 0.06070443916589384000, 0.05911483969839566000,
  1331.             0.05727729210040315000, 0.05519950369998416500, 0.05289018948519365000, 0.05035903555385448000,
  1332.             0.04761665849249054000, 0.04467456085669424000, 0.04154508294346483000, 0.03824135106583080600,
  1333.             0.03477722256477045000, 0.03116722783279807000, 0.02742650970835688000, 0.02357076083932435600,
  1334.             0.01961616045735556700, 0.01557931572294386600, 0.01147723457923446900, 0.00732755390127620800,
  1335.             0.00315334605230596250 };

  1337.         /** Node points. */
  1338.         private final double[] nodePoints;

  1340.         /** Node weights. */
  1341.         private final double[] nodeWeights;

  1343.         /** Number of points. */
  1344.         private final int numberOfPoints;

  1346.         /**
  1347.          * Creates a Gauss integrator of the given order.
  1348.          *
  1349.          * @param numberOfPoints Order of the integration rule.
  1350.          */
  1351.         GaussQuadrature(final int numberOfPoints) {

  1353.             this.numberOfPoints = numberOfPoints;

  1355.             switch (numberOfPoints) {
  1356.                 case 12:
  1357.                     this.nodePoints = P_12.clone();
  1358.                     this.nodeWeights = W_12.clone();
  1359.                     break;
  1360.                 case 16:
  1361.                     this.nodePoints = P_16.clone();
  1362.                     this.nodeWeights = W_16.clone();
  1363.                     break;
  1364.                 case 20:
  1365.                     this.nodePoints = P_20.clone();
  1366.                     this.nodeWeights = W_20.clone();
  1367.                     break;
  1368.                 case 24:
  1369.                     this.nodePoints = P_24.clone();
  1370.                     this.nodeWeights = W_24.clone();
  1371.                     break;
  1372.                 case 32:
  1373.                     this.nodePoints = P_32.clone();
  1374.                     this.nodeWeights = W_32.clone();
  1375.                     break;
  1376.                 case 40:
  1377.                     this.nodePoints = P_40.clone();
  1378.                     this.nodeWeights = W_40.clone();
  1379.                     break;
  1380.                 case 48:
  1381.                 default:
  1382.                     this.nodePoints = P_48.clone();
  1383.                     this.nodeWeights = W_48.clone();
  1384.                     break;
  1385.             }

  1387.         }

  1389.         /**
  1390.          * Integrates a given function on the given interval.
  1391.          *
  1392.          * @param f          Function to integrate.
  1393.          * @param lowerBound Lower bound of the integration interval.
  1394.          * @param upperBound Upper bound of the integration interval.
  1395.          * @return the integral of the weighted function.
  1396.          */
  1397.         public double[] integrate(final UnivariateVectorFunction f, final double lowerBound, final double upperBound) {

  1399.             final double[] adaptedPoints = nodePoints.clone();
  1400.             final double[] adaptedWeights = nodeWeights.clone();
  1401.             transform(adaptedPoints, adaptedWeights, lowerBound, upperBound);
  1402.             return basicIntegrate(f, adaptedPoints, adaptedWeights);
  1403.         }

  1405.         /**
  1406.          * Integrates a given function on the given interval.
  1407.          *
  1408.          * @param <T>        the type of the field elements
  1409.          * @param f          Function to integrate.
  1410.          * @param lowerBound Lower bound of the integration interval.
  1411.          * @param upperBound Upper bound of the integration interval.
  1412.          * @param field      field utilized by default
  1413.          * @return the integral of the weighted function.
  1414.          */
  1415.         public <T extends CalculusFieldElement<T>> T[] integrate(final CalculusFieldUnivariateVectorFunction<T> f,
  1416.                 final T lowerBound, final T upperBound, final Field<T> field) {

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

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

  1423.             for (int i = 0; i < numberOfPoints; i++) {
  1424.                 adaptedPoints[i] = zero.newInstance(nodePoints[i]);
  1425.                 adaptedWeights[i] = zero.newInstance(nodeWeights[i]);
  1426.             }

  1428.             transform(adaptedPoints, adaptedWeights, lowerBound, upperBound);
  1429.             return basicIntegrate(f, adaptedPoints, adaptedWeights, field);
  1430.         }

  1432.         /**
  1433.          * Performs a change of variable so that the integration can be performed on an
  1434.          * arbitrary interval {@code [a, b]}.
  1435.          * <p>
  1436.          * It is assumed that the natural interval is {@code [-1, 1]}.
  1437.          * </p>
  1438.          *
  1439.          * @param points  Points to adapt to the new interval.
  1440.          * @param weights Weights to adapt to the new interval.
  1441.          * @param a       Lower bound of the integration interval.
  1442.          * @param b       Lower bound of the integration interval.
  1443.          */
  1444.         private void transform(final double[] points, final double[] weights, final double a, final double b) {
  1445.             // Scaling
  1446.             final double scale = (b - a) / 2;
  1447.             final double shift = a + scale;
  1448.             for (int i = 0; i < points.length; i++) {
  1449.                 points[i] = points[i] * scale + shift;
  1450.                 weights[i] *= scale;
  1451.             }
  1452.         }

  1454.         /**
  1455.          * Performs a change of variable so that the integration can be performed on an
  1456.          * arbitrary interval {@code [a, b]}.
  1457.          * <p>
  1458.          * It is assumed that the natural interval is {@code [-1, 1]}.
  1459.          * </p>
  1460.          * @param <T>     the type of the field elements
  1461.          * @param points  Points to adapt to the new interval.
  1462.          * @param weights Weights to adapt to the new interval.
  1463.          * @param a       Lower bound of the integration interval.
  1464.          * @param b       Lower bound of the integration interval
  1465.          */
  1466.         private <T extends CalculusFieldElement<T>> void transform(final T[] points, final T[] weights, final T a,
  1467.                 final T b) {
  1468.             // Scaling
  1469.             final T scale = (b.subtract(a)).divide(2.);
  1470.             final T shift = a.add(scale);
  1471.             for (int i = 0; i < points.length; i++) {
  1472.                 points[i] = scale.multiply(points[i]).add(shift);
  1473.                 weights[i] = scale.multiply(weights[i]);
  1474.             }
  1475.         }

  1477.         /**
  1478.          * Returns an estimate of the integral of {@code f(x) * w(x)}, where {@code w}
  1479.          * is a weight function that depends on the actual flavor of the Gauss
  1480.          * integration scheme.
  1481.          *
  1482.          * @param f       Function to integrate.
  1483.          * @param points  Nodes.
  1484.          * @param weights Nodes weights.
  1485.          * @return the integral of the weighted function.
  1486.          */
  1487.         private double[] basicIntegrate(final UnivariateVectorFunction f, final double[] points,
  1488.                 final double[] weights) {
  1489.             double x = points[0];
  1490.             double w = weights[0];
  1491.             double[] v = f.value(x);
  1492.             final double[] y = new double[v.length];
  1493.             for (int j = 0; j < v.length; j++) {
  1494.                 y[j] = w * v[j];
  1495.             }
  1496.             final double[] t = y.clone();
  1497.             final double[] c = new double[v.length];
  1498.             final double[] s = t.clone();
  1499.             for (int i = 1; i < points.length; i++) {
  1500.                 x = points[i];
  1501.                 w = weights[i];
  1502.                 v = f.value(x);
  1503.                 for (int j = 0; j < v.length; j++) {
  1504.                     y[j] = w * v[j] - c[j];
  1505.                     t[j] = s[j] + y[j];
  1506.                     c[j] = (t[j] - s[j]) - y[j];
  1507.                     s[j] = t[j];
  1508.                 }
  1509.             }
  1510.             return s;
  1511.         }

  1513.         /**
  1514.          * Returns an estimate of the integral of {@code f(x) * w(x)}, where {@code w}
  1515.          * is a weight function that depends on the actual flavor of the Gauss
  1516.          * integration scheme.
  1517.          *
  1518.          * @param <T>     the type of the field elements.
  1519.          * @param f       Function to integrate.
  1520.          * @param points  Nodes.
  1521.          * @param weights Nodes weight
  1522.          * @param field   field utilized by default
  1523.          * @return the integral of the weighted function.
  1524.          */
  1525.         private <T extends CalculusFieldElement<T>> T[] basicIntegrate(final CalculusFieldUnivariateVectorFunction<T> f,
  1526.                 final T[] points, final T[] weights, final Field<T> field) {

  1528.             T x = points[0];
  1529.             T w = weights[0];
  1530.             T[] v = f.value(x);

  1532.             final T[] y = MathArrays.buildArray(field, v.length);
  1533.             for (int j = 0; j < v.length; j++) {
  1534.                 y[j] = v[j].multiply(w);
  1535.             }
  1536.             final T[] t = y.clone();
  1537.             final T[] c = MathArrays.buildArray(field, v.length);
  1538.             final T[] s = t.clone();
  1539.             for (int i = 1; i < points.length; i++) {
  1540.                 x = points[i];
  1541.                 w = weights[i];
  1542.                 v = f.value(x);
  1543.                 for (int j = 0; j < v.length; j++) {
  1544.                     y[j] = v[j].multiply(w).subtract(c[j]);
  1545.                     t[j] = y[j].add(s[j]);
  1546.                     c[j] = (t[j].subtract(s[j])).subtract(y[j]);
  1547.                     s[j] = t[j];
  1548.                 }
  1549.             }
  1550.             return s;
  1551.         }

  1553.     }

  1555.     /**
  1556.      * Compute the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  1557.      * coefficients.
  1558.      * <p>
  1559.      * Those coefficients are given in Danielson paper by expression 4.4-(6)
  1560.      * </p>
  1561.      * @author Petre Bazavan
  1562.      * @author Lucian Barbulescu
  1563.      */
  1564.     protected class FourierCjSjCoefficients {

  1566.         /** Maximum possible value for j. */
  1567.         private final int jMax;

  1569.         /**
  1570.          * The C<sub>i</sub><sup>j</sup> coefficients.
  1571.          * <p>
  1572.          * the index i corresponds to the following elements: <br/>
  1573.          * - 0 for a <br>
  1574.          * - 1 for k <br>
  1575.          * - 2 for h <br>
  1576.          * - 3 for q <br>
  1577.          * - 4 for p <br>
  1578.          * - 5 for λ <br>
  1579.          * </p>
  1580.          */
  1581.         private final double[][] cCoef;

  1583.         /**
  1584.          * The C<sub>i</sub><sup>j</sup> coefficients.
  1585.          * <p>
  1586.          * the index i corresponds to the following elements: <br/>
  1587.          * - 0 for a <br>
  1588.          * - 1 for k <br>
  1589.          * - 2 for h <br>
  1590.          * - 3 for q <br>
  1591.          * - 4 for p <br>
  1592.          * - 5 for λ <br>
  1593.          * </p>
  1594.          */
  1595.         private final double[][] sCoef;

  1597.         /**
  1598.          * Standard constructor.
  1599.          * @param state             the current state
  1600.          * @param jMax              maximum value for j
  1601.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1602.          * @param parameters        list of parameter values at state date for each driver
  1603.          * of the force model parameters (1 value per parameter)
  1604.          */
  1605.         FourierCjSjCoefficients(final SpacecraftState state, final int jMax, final AuxiliaryElements auxiliaryElements,
  1606.                 final double[] parameters) {

  1608.             // Initialise the fields
  1609.             this.jMax = jMax;

  1611.             // Allocate the arrays
  1612.             final int rows = jMax + 1;
  1613.             cCoef = new double[rows][6];
  1614.             sCoef = new double[rows][6];

  1616.             // Compute the coefficients
  1617.             computeCoefficients(state, auxiliaryElements, parameters);
  1618.         }

  1620.         /**
  1621.          * Compute the Fourrier coefficients.
  1622.          * <p>
  1623.          * Only the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients
  1624.          * need to be computed as D<sub>i</sub><sup>m</sup> is always 0.
  1625.          * </p>
  1626.          * @param state             the current state
  1627.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1628.          * @param parameters        list of parameter values at state date for each driver
  1629.          * of the force model parameters (1 value per parameter)
  1630.          */
  1631.         private void computeCoefficients(final SpacecraftState state, final AuxiliaryElements auxiliaryElements,
  1632.                 final double[] parameters) {

  1634.             // Computes the limits for the integral
  1635.             final double[] ll = getLLimits(state, auxiliaryElements);
  1636.             // Computes integrated mean element rates if Llow < Lhigh
  1637.             if (ll[0] < ll[1]) {
  1638.                 // Compute 1 / PI
  1639.                 final double ooPI = 1 / FastMath.PI;

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

  1646.                     // divide by PI and set the values for the coefficients
  1647.                     for (int i = 0; i < 6; i++) {
  1648.                         cCoef[j][i] = ooPI * curentCoefficients[i];
  1649.                         sCoef[j][i] = ooPI * curentCoefficients[i + 6];
  1650.                     }
  1651.                 }
  1652.             }
  1653.         }

  1655.         /**
  1656.          * Get the coefficient C<sub>i</sub><sup>j</sup>.
  1657.          * @param i i index - corresponds to the required variation
  1658.          * @param j j index
  1659.          * @return the coefficient C<sub>i</sub><sup>j</sup>
  1660.          */
  1661.         public double getCij(final int i, final int j) {
  1662.             return cCoef[j][i];
  1663.         }

  1665.         /**
  1666.          * Get the coefficient S<sub>i</sub><sup>j</sup>.
  1667.          * @param i i index - corresponds to the required variation
  1668.          * @param j j index
  1669.          * @return the coefficient S<sub>i</sub><sup>j</sup>
  1670.          */
  1671.         public double getSij(final int i, final int j) {
  1672.             return sCoef[j][i];
  1673.         }
  1674.     }

  1676.     /**
  1677.      * Compute the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
  1678.      * coefficients with field elements.
  1679.      * <p>
  1680.      * Those coefficients are given in Danielson paper by expression 4.4-(6)
  1681.      * </p>
  1682.      * @author Petre Bazavan
  1683.      * @author Lucian Barbulescu
  1684.      * @param <T> type of the field elements
  1685.      */
  1686.     protected class FieldFourierCjSjCoefficients<T extends CalculusFieldElement<T>> {

  1688.         /** Maximum possible value for j. */
  1689.         private final int jMax;

  1691.         /**
  1692.          * The C<sub>i</sub><sup>j</sup> coefficients.
  1693.          * <p>
  1694.          * the index i corresponds to the following elements: <br/>
  1695.          * - 0 for a <br>
  1696.          * - 1 for k <br>
  1697.          * - 2 for h <br>
  1698.          * - 3 for q <br>
  1699.          * - 4 for p <br>
  1700.          * - 5 for λ <br>
  1701.          * </p>
  1702.          */
  1703.         private final T[][] cCoef;

  1705.         /**
  1706.          * The C<sub>i</sub><sup>j</sup> coefficients.
  1707.          * <p>
  1708.          * the index i corresponds to the following elements: <br/>
  1709.          * - 0 for a <br>
  1710.          * - 1 for k <br>
  1711.          * - 2 for h <br>
  1712.          * - 3 for q <br>
  1713.          * - 4 for p <br>
  1714.          * - 5 for λ <br>
  1715.          * </p>
  1716.          */
  1717.         private final T[][] sCoef;

  1719.         /**
  1720.          * Standard constructor.
  1721.          * @param state             the current state
  1722.          * @param jMax              maximum value for j
  1723.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1724.          * @param parameters        values of the force model parameters
  1725.          * @param field             field used by default
  1726.          */
  1727.         FieldFourierCjSjCoefficients(final FieldSpacecraftState<T> state, final int jMax,
  1728.                 final FieldAuxiliaryElements<T> auxiliaryElements, final T[] parameters, final Field<T> field) {
  1729.             // Initialise the fields
  1730.             this.jMax = jMax;

  1732.             // Allocate the arrays
  1733.             final int rows = jMax + 1;
  1734.             cCoef = MathArrays.buildArray(field, rows, 6);
  1735.             sCoef = MathArrays.buildArray(field, rows, 6);

  1737.             // Compute the coefficients
  1738.             computeCoefficients(state, auxiliaryElements, parameters, field);
  1739.         }

  1741.         /**
  1742.          * Compute the Fourrier coefficients.
  1743.          * <p>
  1744.          * Only the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients
  1745.          * need to be computed as D<sub>i</sub><sup>m</sup> is always 0.
  1746.          * </p>
  1747.          * @param state             the current state
  1748.          * @param auxiliaryElements auxiliary elements related to the current orbit
  1749.          * @param parameters        values of the force model parameters
  1750.          * @param field             field used by default
  1751.          */
  1752.         private void computeCoefficients(final FieldSpacecraftState<T> state,
  1753.                 final FieldAuxiliaryElements<T> auxiliaryElements, final T[] parameters, final Field<T> field) {
  1754.             // Zero
  1755.             final T zero = field.getZero();
  1756.             // Computes the limits for the integral
  1757.             final T[] ll = getLLimits(state, auxiliaryElements);
  1758.             // Computes integrated mean element rates if Llow < Lhigh
  1759.             if (ll[0].getReal() < ll[1].getReal()) {
  1760.                 // Compute 1 / PI
  1761.                 final T ooPI = zero.getPi().reciprocal();

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

  1768.                     // divide by PI and set the values for the coefficients
  1769.                     for (int i = 0; i < 6; i++) {
  1770.                         cCoef[j][i] = curentCoefficients[i].multiply(ooPI);
  1771.                         sCoef[j][i] = curentCoefficients[i + 6].multiply(ooPI);
  1772.                     }
  1773.                 }
  1774.             }
  1775.         }

  1777.         /**
  1778.          * Get the coefficient C<sub>i</sub><sup>j</sup>.
  1779.          * @param i i index - corresponds to the required variation
  1780.          * @param j j index
  1781.          * @return the coefficient C<sub>i</sub><sup>j</sup>
  1782.          */
  1783.         public T getCij(final int i, final int j) {
  1784.             return cCoef[j][i];
  1785.         }

  1787.         /**
  1788.          * Get the coefficient S<sub>i</sub><sup>j</sup>.
  1789.          * @param i i index - corresponds to the required variation
  1790.          * @param j j index
  1791.          * @return the coefficient S<sub>i</sub><sup>j</sup>
  1792.          */
  1793.         public T getSij(final int i, final int j) {
  1794.             return sCoef[j][i];
  1795.         }
  1796.     }

  1798.     /**
  1799.      * This class handles the short periodic coefficients described in Danielson
  1800.      * 2.5.3-26.
  1801.      *
  1802.      * <p>
  1803.      * The value of M is 0. Also, since the values of the Fourier coefficient
  1804.      * D<sub>i</sub><sup>m</sup> is 0 then the values of the coefficients
  1805.      * D<sub>i</sub><sup>m</sup> for m &gt; 2 are also 0.
  1806.      * </p>
  1807.      * @author Petre Bazavan
  1808.      * @author Lucian Barbulescu
  1809.      *
  1810.      */
  1811.     protected static class GaussianShortPeriodicCoefficients implements ShortPeriodTerms {

  1813.         /** Maximum value for j index. */
  1814.         private final int jMax;

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

  1819.         /** Prefix for coefficients keys. */
  1820.         private final String coefficientsKeyPrefix;

  1822.         /** All coefficients slots. */
  1823.         private final transient TimeSpanMap<Slot> slots;

  1825.         /**
  1826.          * Constructor.
  1827.          * @param coefficientsKeyPrefix prefix for coefficients keys
  1828.          * @param jMax                  maximum value for j index
  1829.          * @param interpolationPoints   number of points used in the interpolation
  1830.          *                              process
  1831.          * @param slots                 all coefficients slots
  1832.          */
  1833.         GaussianShortPeriodicCoefficients(final String coefficientsKeyPrefix, final int jMax,
  1834.                 final int interpolationPoints, final TimeSpanMap<Slot> slots) {
  1835.             // Initialize fields
  1836.             this.jMax = jMax;
  1837.             this.interpolationPoints = interpolationPoints;
  1838.             this.coefficientsKeyPrefix = coefficientsKeyPrefix;
  1839.             this.slots = slots;
  1840.         }

  1842.         /**
  1843.          * Get the slot valid for some date.
  1844.          * @param meanStates mean states defining the slot
  1845.          * @return slot valid at the specified date
  1846.          */
  1847.         public Slot createSlot(final SpacecraftState... meanStates) {
  1848.             final Slot slot = new Slot(jMax, interpolationPoints);
  1849.             final AbsoluteDate first = meanStates[0].getDate();
  1850.             final AbsoluteDate last = meanStates[meanStates.length - 1].getDate();
  1851.             final int compare = first.compareTo(last);
  1852.             if (compare < 0) {
  1853.                 slots.addValidAfter(slot, first, false);
  1854.             } else if (compare > 0) {
  1855.                 slots.addValidBefore(slot, first, false);
  1856.             } else {
  1857.                 // single date, valid for all time
  1858.                 slots.addValidAfter(slot, AbsoluteDate.PAST_INFINITY, false);
  1859.             }
  1860.             return slot;
  1861.         }

  1863.         /**
  1864.          * Compute the short periodic coefficients.
  1865.          *
  1866.          * @param state       current state information: date, kinematics, attitude
  1867.          * @param slot        coefficients slot
  1868.          * @param fourierCjSj Fourier coefficients
  1869.          * @param uijvij      U and V coefficients
  1870.          * @param n           Keplerian mean motion
  1871.          * @param a           semi major axis
  1872.          */
  1873.         private void computeCoefficients(final SpacecraftState state, final Slot slot,
  1874.                 final FourierCjSjCoefficients fourierCjSj, final UijVijCoefficients uijvij, final double n,
  1875.                 final double a) {

  1877.             // get the current date
  1878.             final AbsoluteDate date = state.getDate();

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

  1883.             // 1. / n
  1884.             final double oon = 1. / n;
  1885.             // 3. / (2 * a * n)
  1886.             final double to2an = 1.5 * oon / a;
  1887.             // 3. / (4 * a * n)
  1888.             final double to4an = to2an / 2;

  1890.             // Compute the coefficients for each element
  1891.             final int size = jMax + 1;
  1892.             final double[] di1 = new double[6];
  1893.             final double[] di2 = new double[6];
  1894.             final double[][] currentCij = new double[size][6];
  1895.             final double[][] currentSij = new double[size][6];
  1896.             for (int i = 0; i < 6; i++) {

  1898.                 // compute D<sub>i</sub>¹ and D<sub>i</sub>² (all others are 0)
  1899.                 di1[i] = -oon * fourierCjSj.getCij(i, 0);
  1900.                 if (i == 5) {
  1901.                     di1[i] += to2an * uijvij.getU1(0, 0);
  1902.                 }
  1903.                 di2[i] = 0.;
  1904.                 if (i == 5) {
  1905.                     di2[i] += -to4an * fourierCjSj.getCij(0, 0);
  1906.                 }

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

  1911.                 for (int j = 1; j <= jMax; j++) {
  1912.                     // compute the current C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup>
  1913.                     currentCij[j][i] = oon * uijvij.getU1(j, i);
  1914.                     if (i == 5) {
  1915.                         currentCij[j][i] += -to2an * uijvij.getU2(j);
  1916.                     }
  1917.                     currentSij[j][i] = oon * uijvij.getV1(j, i);
  1918.                     if (i == 5) {
  1919.                         currentSij[j][i] += -to2an * uijvij.getV2(j);
  1920.                     }

  1922.                     // add the computed coefficients to C<sub>i</sub>⁰
  1923.                     currentCij[0][i] -= currentCij[j][i] * uijvij.currentRhoSigmaj[0][j] +
  1924.                         currentSij[j][i] * uijvij.currentRhoSigmaj[1][j];
  1925.                 }

  1927.             }

  1929.             // add the values to the interpolators
  1930.             slot.cij[0].addGridPoint(date, currentCij[0]);
  1931.             slot.dij[1].addGridPoint(date, di1);
  1932.             slot.dij[2].addGridPoint(date, di2);
  1933.             for (int j = 1; j <= jMax; j++) {
  1934.                 slot.cij[j].addGridPoint(date, currentCij[j]);
  1935.                 slot.sij[j].addGridPoint(date, currentSij[j]);
  1936.             }

  1938.         }

  1940.         /**
  1941.          * Compute the coefficient k₂⁰ by using the equation 2.5.3-(9a) from Danielson.
  1942.          * <p>
  1943.          * After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:<br>
  1944.          * k₂⁰ = &Sigma;<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² +
  1945.          * ρ<sub>k</sub>²)]
  1946.          * </p>
  1947.          * @param kMax             max value fot k index
  1948.          * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and
  1949.          *                         σ<sub>j</sub> coefficients
  1950.          * @return the coefficient k₂⁰
  1951.          */
  1952.         private double computeK20(final int kMax, final double[][] currentRhoSigmaj) {
  1953.             double k20 = 0.;

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

  1962.                 // multiply by 2 / k²
  1963.                 currentTerm *= 2. / (kIndex * kIndex);

  1965.                 // add the term to the result
  1966.                 k20 += currentTerm;
  1967.             }

  1969.             return k20;
  1970.         }

  1972.         /** {@inheritDoc} */
  1973.         @Override
  1974.         public double[] value(final Orbit meanOrbit) {

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

  1979.             // Get the True longitude L
  1980.             final double L = meanOrbit.getLv();

  1982.             // Compute the center (l - λ)
  1983.             final double center = L - meanOrbit.getLM();
  1984.             // Compute (l - λ)²
  1985.             final double center2 = center * center;

  1987.             // Initialize short periodic variations
  1988.             final double[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  1989.             final double[] d1 = slot.dij[1].value(meanOrbit.getDate());
  1990.             final double[] d2 = slot.dij[2].value(meanOrbit.getDate());
  1991.             for (int i = 0; i < 6; i++) {
  1992.                 shortPeriodicVariation[i] += center * d1[i] + center2 * d2[i];
  1993.             }

  1995.             for (int j = 1; j <= JMAX; j++) {
  1996.                 final double[] c = slot.cij[j].value(meanOrbit.getDate());
  1997.                 final double[] s = slot.sij[j].value(meanOrbit.getDate());
  1998.                 final SinCos sc  = FastMath.sinCos(j * L);
  1999.                 final double cos = sc.cos();
  2000.                 final double sin = sc.sin();
  2001.                 for (int i = 0; i < 6; i++) {
  2002.                     // add corresponding term to the short periodic variation
  2003.                     shortPeriodicVariation[i] += c[i] * cos;
  2004.                     shortPeriodicVariation[i] += s[i] * sin;
  2005.                 }
  2006.             }

  2008.             return shortPeriodicVariation;

  2010.         }

  2012.         /** {@inheritDoc} */
  2013.         public String getCoefficientsKeyPrefix() {
  2014.             return coefficientsKeyPrefix;
  2015.         }

  2017.         /**
  2018.          * {@inheritDoc}
  2019.          * <p>
  2020.          * For Gaussian forces, there are JMAX cj coefficients, JMAX sj coefficients and
  2021.          * 3 dj coefficients. As JMAX = 12, this sums up to 27 coefficients. The j index
  2022.          * is the integer multiplier for the true longitude argument in the cj and sj
  2023.          * coefficients and to the degree in the polynomial dj coefficients.
  2024.          * </p>
  2025.          */
  2026.         @Override
  2027.         public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected) {

  2029.             // select the coefficients slot
  2030.             final Slot slot = slots.get(date);

  2032.             final Map<String, double[]> coefficients = new HashMap<>(2 * JMAX + 3);
  2033.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  2034.             storeIfSelected(coefficients, selected, slot.dij[1].value(date), "d", 1);
  2035.             storeIfSelected(coefficients, selected, slot.dij[2].value(date), "d", 2);
  2036.             for (int j = 1; j <= JMAX; j++) {
  2037.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  2038.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  2039.             }

  2041.             return coefficients;

  2043.         }

  2045.         /**
  2046.          * Put a coefficient in a map if selected.
  2047.          * @param map      map to populate
  2048.          * @param selected set of coefficients that should be put in the map (empty set
  2049.          *                 means all coefficients are selected)
  2050.          * @param value    coefficient value
  2051.          * @param id       coefficient identifier
  2052.          * @param indices  list of coefficient indices
  2053.          */
  2054.         private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected, final double[] value,
  2055.                 final String id, final int... indices) {
  2056.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  2057.             keyBuilder.append(id);
  2058.             for (int index : indices) {
  2059.                 keyBuilder.append('[').append(index).append(']');
  2060.             }
  2061.             final String key = keyBuilder.toString();
  2062.             if (selected.isEmpty() || selected.contains(key)) {
  2063.                 map.put(key, value);
  2064.             }
  2065.         }

  2067.     }

  2069.     /**
  2070.      * This class handles the short periodic coefficients described in Danielson
  2071.      * 2.5.3-26.
  2072.      *
  2073.      * <p>
  2074.      * The value of M is 0. Also, since the values of the Fourier coefficient
  2075.      * D<sub>i</sub><sup>m</sup> is 0 then the values of the coefficients
  2076.      * D<sub>i</sub><sup>m</sup> for m &gt; 2 are also 0.
  2077.      * </p>
  2078.      * @author Petre Bazavan
  2079.      * @author Lucian Barbulescu
  2080.      * @param <T> type of the field elements
  2081.      */
  2082.     protected static class FieldGaussianShortPeriodicCoefficients<T extends CalculusFieldElement<T>>
  2083.             implements FieldShortPeriodTerms<T> {

  2085.         /** Maximum value for j index. */
  2086.         private final int jMax;

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

  2091.         /** Prefix for coefficients keys. */
  2092.         private final String coefficientsKeyPrefix;

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

  2097.         /**
  2098.          * Constructor.
  2099.          * @param coefficientsKeyPrefix prefix for coefficients keys
  2100.          * @param jMax                  maximum value for j index
  2101.          * @param interpolationPoints   number of points used in the interpolation
  2102.          *                              process
  2103.          * @param slots                 all coefficients slots
  2104.          */
  2105.         FieldGaussianShortPeriodicCoefficients(final String coefficientsKeyPrefix, final int jMax,
  2106.                 final int interpolationPoints, final FieldTimeSpanMap<FieldSlot<T>, T> slots) {
  2107.             // Initialize fields
  2108.             this.jMax = jMax;
  2109.             this.interpolationPoints = interpolationPoints;
  2110.             this.coefficientsKeyPrefix = coefficientsKeyPrefix;
  2111.             this.slots = slots;
  2112.         }

  2114.         /**
  2115.          * Get the slot valid for some date.
  2116.          * @param meanStates mean states defining the slot
  2117.          * @return slot valid at the specified date
  2118.          */
  2119.         @SuppressWarnings("unchecked")
  2120.         public FieldSlot<T> createSlot(final FieldSpacecraftState<T>... meanStates) {
  2121.             final FieldSlot<T> slot = new FieldSlot<>(jMax, interpolationPoints);
  2122.             final FieldAbsoluteDate<T> first = meanStates[0].getDate();
  2123.             final FieldAbsoluteDate<T> last = meanStates[meanStates.length - 1].getDate();
  2124.             if (first.compareTo(last) <= 0) {
  2125.                 slots.addValidAfter(slot, first);
  2126.             } else {
  2127.                 slots.addValidBefore(slot, first);
  2128.             }
  2129.             return slot;
  2130.         }

  2132.         /**
  2133.          * Compute the short periodic coefficients.
  2134.          *
  2135.          * @param state       current state information: date, kinematics, attitude
  2136.          * @param slot        coefficients slot
  2137.          * @param fourierCjSj Fourier coefficients
  2138.          * @param uijvij      U and V coefficients
  2139.          * @param n           Keplerian mean motion
  2140.          * @param a           semi major axis
  2141.          * @param field       field used by default
  2142.          */
  2143.         private void computeCoefficients(final FieldSpacecraftState<T> state, final FieldSlot<T> slot,
  2144.                 final FieldFourierCjSjCoefficients<T> fourierCjSj, final FieldUijVijCoefficients<T> uijvij, final T n,
  2145.                 final T a, final Field<T> field) {

  2147.             // Zero
  2148.             final T zero = field.getZero();

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

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

  2156.             // 1. / n
  2157.             final T oon = n.reciprocal();
  2158.             // 3. / (2 * a * n)
  2159.             final T to2an = oon.multiply(1.5).divide(a);
  2160.             // 3. / (4 * a * n)
  2161.             final T to4an = to2an.divide(2.);

  2163.             // Compute the coefficients for each element
  2164.             final int size = jMax + 1;
  2165.             final T[] di1 = MathArrays.buildArray(field, 6);
  2166.             final T[] di2 = MathArrays.buildArray(field, 6);
  2167.             final T[][] currentCij = MathArrays.buildArray(field, size, 6);
  2168.             final T[][] currentSij = MathArrays.buildArray(field, size, 6);
  2169.             for (int i = 0; i < 6; i++) {

  2171.                 // compute D<sub>i</sub>¹ and D<sub>i</sub>² (all others are 0)
  2172.                 di1[i] = oon.negate().multiply(fourierCjSj.getCij(i, 0));
  2173.                 if (i == 5) {
  2174.                     di1[i] = di1[i].add(to2an.multiply(uijvij.getU1(0, 0)));
  2175.                 }
  2176.                 di2[i] = zero;
  2177.                 if (i == 5) {
  2178.                     di2[i] = di2[i].add(to4an.negate().multiply(fourierCjSj.getCij(0, 0)));
  2179.                 }

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

  2184.                 for (int j = 1; j <= jMax; j++) {
  2185.                     // compute the current C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup>
  2186.                     currentCij[j][i] = oon.multiply(uijvij.getU1(j, i));
  2187.                     if (i == 5) {
  2188.                         currentCij[j][i] = currentCij[j][i].add(to2an.negate().multiply(uijvij.getU2(j)));
  2189.                     }
  2190.                     currentSij[j][i] = oon.multiply(uijvij.getV1(j, i));
  2191.                     if (i == 5) {
  2192.                         currentSij[j][i] = currentSij[j][i].add(to2an.negate().multiply(uijvij.getV2(j)));
  2193.                     }

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

  2200.             }

  2202.             // add the values to the interpolators
  2203.             slot.cij[0].addGridPoint(date, currentCij[0]);
  2204.             slot.dij[1].addGridPoint(date, di1);
  2205.             slot.dij[2].addGridPoint(date, di2);
  2206.             for (int j = 1; j <= jMax; j++) {
  2207.                 slot.cij[j].addGridPoint(date, currentCij[j]);
  2208.                 slot.sij[j].addGridPoint(date, currentSij[j]);
  2209.             }

  2211.         }

  2213.         /**
  2214.          * Compute the coefficient k₂⁰ by using the equation 2.5.3-(9a) from Danielson.
  2215.          * <p>
  2216.          * After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:<br>
  2217.          * k₂⁰ = &Sigma;<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² +
  2218.          * ρ<sub>k</sub>²)]
  2219.          * </p>
  2220.          * @param kMax             max value fot k index
  2221.          * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and
  2222.          *                         σ<sub>j</sub> coefficients
  2223.          * @param field            field used by default
  2224.          * @return the coefficient k₂⁰
  2225.          */
  2226.         private T computeK20(final int kMax, final T[][] currentRhoSigmaj, final Field<T> field) {
  2227.             T k20 = field.getZero();

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

  2236.                 // multiply by 2 / k²
  2237.                 currentTerm = currentTerm.multiply(2. / (kIndex * kIndex));

  2239.                 // add the term to the result
  2240.                 k20 = k20.add(currentTerm);
  2241.             }

  2243.             return k20;
  2244.         }

  2246.         /** {@inheritDoc} */
  2247.         @Override
  2248.         public T[] value(final FieldOrbit<T> meanOrbit) {

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

  2253.             // Get the True longitude L
  2254.             final T L = meanOrbit.getLv();

  2256.             // Compute the center (l - λ)
  2257.             final T center = L.subtract(meanOrbit.getLM());
  2258.             // Compute (l - λ)²
  2259.             final T center2 = center.square();

  2261.             // Initialize short periodic variations
  2262.             final T[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
  2263.             final T[] d1 = slot.dij[1].value(meanOrbit.getDate());
  2264.             final T[] d2 = slot.dij[2].value(meanOrbit.getDate());
  2265.             for (int i = 0; i < 6; i++) {
  2266.                 shortPeriodicVariation[i] = shortPeriodicVariation[i]
  2267.                         .add(center.multiply(d1[i]).add(center2.multiply(d2[i])));
  2268.             }

  2270.             for (int j = 1; j <= JMAX; j++) {
  2271.                 final T[] c = slot.cij[j].value(meanOrbit.getDate());
  2272.                 final T[] s = slot.sij[j].value(meanOrbit.getDate());
  2273.                 final FieldSinCos<T> sc = FastMath.sinCos(L.multiply(j));
  2274.                 final T cos = sc.cos();
  2275.                 final T sin = sc.sin();
  2276.                 for (int i = 0; i < 6; i++) {
  2277.                     // add corresponding term to the short periodic variation
  2278.                     shortPeriodicVariation[i] = shortPeriodicVariation[i].add(c[i].multiply(cos));
  2279.                     shortPeriodicVariation[i] = shortPeriodicVariation[i].add(s[i].multiply(sin));
  2280.                 }
  2281.             }

  2283.             return shortPeriodicVariation;

  2285.         }

  2287.         /** {@inheritDoc} */
  2288.         public String getCoefficientsKeyPrefix() {
  2289.             return coefficientsKeyPrefix;
  2290.         }

  2292.         /**
  2293.          * {@inheritDoc}
  2294.          * <p>
  2295.          * For Gaussian forces, there are JMAX cj coefficients, JMAX sj coefficients and
  2296.          * 3 dj coefficients. As JMAX = 12, this sums up to 27 coefficients. The j index
  2297.          * is the integer multiplier for the true longitude argument in the cj and sj
  2298.          * coefficients and to the degree in the polynomial dj coefficients.
  2299.          * </p>
  2300.          */
  2301.         @Override
  2302.         public Map<String, T[]> getCoefficients(final FieldAbsoluteDate<T> date, final Set<String> selected) {

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

  2307.             final Map<String, T[]> coefficients = new HashMap<>(2 * JMAX + 3);
  2308.             storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
  2309.             storeIfSelected(coefficients, selected, slot.dij[1].value(date), "d", 1);
  2310.             storeIfSelected(coefficients, selected, slot.dij[2].value(date), "d", 2);
  2311.             for (int j = 1; j <= JMAX; j++) {
  2312.                 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
  2313.                 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
  2314.             }

  2316.             return coefficients;

  2318.         }

  2320.         /**
  2321.          * Put a coefficient in a map if selected.
  2322.          * @param map      map to populate
  2323.          * @param selected set of coefficients that should be put in the map (empty set
  2324.          *                 means all coefficients are selected)
  2325.          * @param value    coefficient value
  2326.          * @param id       coefficient identifier
  2327.          * @param indices  list of coefficient indices
  2328.          */
  2329.         private void storeIfSelected(final Map<String, T[]> map, final Set<String> selected, final T[] value,
  2330.                 final String id, final int... indices) {
  2331.             final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
  2332.             keyBuilder.append(id);
  2333.             for (int index : indices) {
  2334.                 keyBuilder.append('[').append(index).append(']');
  2335.             }
  2336.             final String key = keyBuilder.toString();
  2337.             if (selected.isEmpty() || selected.contains(key)) {
  2338.                 map.put(key, value);
  2339.             }
  2340.         }

  2342.     }

  2344.     /**
  2345.      * The U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients
  2346.      * described by equations 2.5.3-(21) and 2.5.3-(22) from Danielson.
  2347.      * <p>
  2348.      * The index i takes only the values 1 and 2<br>
  2349.      * For U only the index 0 for j is used.
  2350.      * </p>
  2351.      *
  2352.      * @author Petre Bazavan
  2353.      * @author Lucian Barbulescu
  2354.      */
  2355.     protected static class UijVijCoefficients {

  2357.         /**
  2358.          * The U₁<sup>j</sup> coefficients.
  2359.          * <p>
  2360.          * The first index identifies the Fourier coefficients used<br>
  2361.          * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and
  2362.          * S<sub>i</sub><sup>j</sup><br>
  2363.          * The only exception is when j = 0 when only the coefficient for fourier index
  2364.          * = 1 (i == 0) is needed.<br>
  2365.          * Also, for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are
  2366.          * computed, because are required to compute the coefficients U₂<sup>j</sup>
  2367.          * </p>
  2368.          */
  2369.         private final double[][] u1ij;

  2371.         /**
  2372.          * The V₁<sup>j</sup> coefficients.
  2373.          * <p>
  2374.          * The first index identifies the Fourier coefficients used<br>
  2375.          * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and
  2376.          * S<sub>i</sub><sup>j</sup><br>
  2377.          * for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed,
  2378.          * because are required to compute the coefficients V₂<sup>j</sup>
  2379.          * </p>
  2380.          */
  2381.         private final double[][] v1ij;

  2383.         /**
  2384.          * The U₂<sup>j</sup> coefficients.
  2385.          * <p>
  2386.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as
  2387.          * they are the only ones required.
  2388.          * </p>
  2389.          */
  2390.         private final double[] u2ij;

  2392.         /**
  2393.          * The V₂<sup>j</sup> coefficients.
  2394.          * <p>
  2395.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as
  2396.          * they are the only ones required.
  2397.          * </p>
  2398.          */
  2399.         private final double[] v2ij;

  2401.         /**
  2402.          * The current computed values for the ρ<sub>j</sub> and σ<sub>j</sub>
  2403.          * coefficients.
  2404.          */
  2405.         private final double[][] currentRhoSigmaj;

  2407.         /**
  2408.          * The C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> Fourier
  2409.          * coefficients.
  2410.          */
  2411.         private final FourierCjSjCoefficients fourierCjSj;

  2413.         /** The maximum value for j index. */
  2414.         private final int jMax;

  2416.         /**
  2417.          * Constructor.
  2418.          * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and
  2419.          *                         σ<sub>j</sub> coefficients
  2420.          * @param fourierCjSj      the fourier coefficients C<sub>i</sub><sup>j</sup>
  2421.          *                         and the S<sub>i</sub><sup>j</sup>
  2422.          * @param jMax             maximum value for j index
  2423.          */
  2424.         UijVijCoefficients(final double[][] currentRhoSigmaj, final FourierCjSjCoefficients fourierCjSj,
  2425.                 final int jMax) {
  2426.             this.currentRhoSigmaj = currentRhoSigmaj;
  2427.             this.fourierCjSj = fourierCjSj;
  2428.             this.jMax = jMax;

  2430.             // initialize the internal arrays.
  2431.             this.u1ij = new double[6][2 * jMax + 1];
  2432.             this.v1ij = new double[6][2 * jMax + 1];
  2433.             this.u2ij = new double[jMax + 1];
  2434.             this.v2ij = new double[jMax + 1];

  2436.             // compute the coefficients
  2437.             computeU1V1Coefficients();
  2438.             computeU2V2Coefficients();
  2439.         }

  2441.         /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients. */
  2442.         private void computeU1V1Coefficients() {
  2443.             // generate the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients
  2444.             // for j >= 1
  2445.             // also the U₁⁰ for Fourier index = 1 (i == 0) coefficient will be computed
  2446.             u1ij[0][0] = 0;
  2447.             for (int j = 1; j <= jMax; j++) {
  2448.                 // compute 1 / j
  2449.                 final double ooj = 1. / j;

  2451.                 for (int i = 0; i < 6; i++) {
  2452.                     // j is aready between 1 and J
  2453.                     u1ij[i][j] = fourierCjSj.getSij(i, j);
  2454.                     v1ij[i][j] = fourierCjSj.getCij(i, j);

  2456.                     // 1 - δ<sub>1j</sub> is 1 for all j > 1
  2457.                     if (j > 1) {
  2458.                         // k starts with 1 because j-J is less than or equal to 0
  2459.                         for (int kIndex = 1; kIndex <= j - 1; kIndex++) {
  2460.                             // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
  2461.                             // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
  2462.                             u1ij[i][j] += fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[1][kIndex] +
  2463.                                     fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[0][kIndex];

  2465.                             // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  2466.                             // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  2467.                             v1ij[i][j] += fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[0][kIndex] -
  2468.                                     fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[1][kIndex];
  2469.                         }
  2470.                     }

  2472.                     // since j must be between 1 and J-1 and is already between 1 and J
  2473.                     // the following sum is skiped only for j = jMax
  2474.                     if (j != jMax) {
  2475.                         for (int kIndex = 1; kIndex <= jMax - j; kIndex++) {
  2476.                             // -C<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub> +
  2477.                             // S<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub>
  2478.                             u1ij[i][j] += -fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[1][kIndex] +
  2479.                                     fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[0][kIndex];

  2481.                             // C<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub> +
  2482.                             // S<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub>
  2483.                             v1ij[i][j] += fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[0][kIndex] +
  2484.                                     fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[1][kIndex];
  2485.                         }
  2486.                     }

  2488.                     for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  2489.                         // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  2490.                         // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  2491.                         u1ij[i][j] += -fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[1][j + kIndex] -
  2492.                                 fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[0][j + kIndex];

  2494.                         // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  2495.                         // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  2496.                         v1ij[i][j] += fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[0][j + kIndex] +
  2497.                                 fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[1][j + kIndex];
  2498.                     }

  2500.                     // divide by 1 / j
  2501.                     u1ij[i][j] *= -ooj;
  2502.                     v1ij[i][j] *= ooj;

  2504.                     // if index = 1 (i == 0) add the computed terms to U₁⁰
  2505.                     if (i == 0) {
  2506.                         // - (U₁<sup>j</sup> * ρ<sub>j</sub> + V₁<sup>j</sup> * σ<sub>j</sub>
  2507.                         u1ij[0][0] += -u1ij[0][j] * currentRhoSigmaj[0][j] - v1ij[0][j] * currentRhoSigmaj[1][j];
  2508.                     }
  2509.                 }
  2510.             }

  2512.             // Terms with j > jMax are required only when computing the coefficients
  2513.             // U₂<sup>j</sup> and V₂<sup>j</sup>
  2514.             // and those coefficients are only required for Fourier index = 1 (i == 0).
  2515.             for (int j = jMax + 1; j <= 2 * jMax; j++) {
  2516.                 // compute 1 / j
  2517.                 final double ooj = 1. / j;
  2518.                 // the value of i is 0
  2519.                 u1ij[0][j] = 0.;
  2520.                 v1ij[0][j] = 0.;

  2522.                 // k starts from j-J as it is always greater than or equal to 1
  2523.                 for (int kIndex = j - jMax; kIndex <= j - 1; kIndex++) {
  2524.                     // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
  2525.                     // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
  2526.                     u1ij[0][j] += fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[1][kIndex] +
  2527.                             fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[0][kIndex];

  2529.                     // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  2530.                     // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  2531.                     v1ij[0][j] += fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[0][kIndex] -
  2532.                             fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[1][kIndex];
  2533.                 }
  2534.                 for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  2535.                     // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  2536.                     // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  2537.                     u1ij[0][j] += -fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[1][j + kIndex] -
  2538.                             fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[0][j + kIndex];

  2540.                     // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  2541.                     // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  2542.                     v1ij[0][j] += fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[0][j + kIndex] +
  2543.                             fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[1][j + kIndex];
  2544.                 }

  2546.                 // divide by 1 / j
  2547.                 u1ij[0][j] *= -ooj;
  2548.                 v1ij[0][j] *= ooj;
  2549.             }
  2550.         }

  2552.         /**
  2553.          * Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
  2554.          * <p>
  2555.          * Only the coefficients for Fourier index = 1 (i == 0) are required.
  2556.          * </p>
  2557.          */
  2558.         private void computeU2V2Coefficients() {
  2559.             for (int j = 1; j <= jMax; j++) {
  2560.                 // compute 1 / j
  2561.                 final double ooj = 1. / j;

  2563.                 // only the values for i == 0 are computed
  2564.                 u2ij[j] = v1ij[0][j];
  2565.                 v2ij[j] = u1ij[0][j];

  2567.                 // 1 - δ<sub>1j</sub> is 1 for all j > 1
  2568.                 if (j > 1) {
  2569.                     for (int l = 1; l <= j - 1; l++) {
  2570.                         // U₁<sup>j-l</sup> * σ<sub>l</sub> +
  2571.                         // V₁<sup>j-l</sup> * ρ<sub>l</sub>
  2572.                         u2ij[j] += u1ij[0][j - l] * currentRhoSigmaj[1][l] + v1ij[0][j - l] * currentRhoSigmaj[0][l];

  2574.                         // U₁<sup>j-l</sup> * ρ<sub>l</sub> -
  2575.                         // V₁<sup>j-l</sup> * σ<sub>l</sub>
  2576.                         v2ij[j] += u1ij[0][j - l] * currentRhoSigmaj[0][l] - v1ij[0][j - l] * currentRhoSigmaj[1][l];
  2577.                     }
  2578.                 }

  2580.                 for (int l = 1; l <= jMax; l++) {
  2581.                     // -U₁<sup>j+l</sup> * σ<sub>l</sub> +
  2582.                     // U₁<sup>l</sup> * σ<sub>j+l</sub> +
  2583.                     // V₁<sup>j+l</sup> * ρ<sub>l</sub> -
  2584.                     // V₁<sup>l</sup> * ρ<sub>j+l</sub>
  2585.                     u2ij[j] += -u1ij[0][j + l] * currentRhoSigmaj[1][l] + u1ij[0][l] * currentRhoSigmaj[1][j + l] +
  2586.                             v1ij[0][j + l] * currentRhoSigmaj[0][l] - v1ij[0][l] * currentRhoSigmaj[0][j + l];

  2588.                     // U₁<sup>j+l</sup> * ρ<sub>l</sub> +
  2589.                     // U₁<sup>l</sup> * ρ<sub>j+l</sub> +
  2590.                     // V₁<sup>j+l</sup> * σ<sub>l</sub> +
  2591.                     // V₁<sup>l</sup> * σ<sub>j+l</sub>
  2592.                     u2ij[j] += u1ij[0][j + l] * currentRhoSigmaj[0][l] + u1ij[0][l] * currentRhoSigmaj[0][j + l] +
  2593.                             v1ij[0][j + l] * currentRhoSigmaj[1][l] + v1ij[0][l] * currentRhoSigmaj[1][j + l];
  2594.                 }

  2596.                 // divide by 1 / j
  2597.                 u2ij[j] *= -ooj;
  2598.                 v2ij[j] *= ooj;
  2599.             }
  2600.         }

  2602.         /**
  2603.          * Get the coefficient U₁<sup>j</sup> for Fourier index i.
  2604.          *
  2605.          * @param j j index
  2606.          * @param i Fourier index (starts at 0)
  2607.          * @return the coefficient U₁<sup>j</sup> for the given Fourier index i
  2608.          */
  2609.         public double getU1(final int j, final int i) {
  2610.             return u1ij[i][j];
  2611.         }

  2613.         /**
  2614.          * Get the coefficient V₁<sup>j</sup> for Fourier index i.
  2615.          *
  2616.          * @param j j index
  2617.          * @param i Fourier index (starts at 0)
  2618.          * @return the coefficient V₁<sup>j</sup> for the given Fourier index i
  2619.          */
  2620.         public double getV1(final int j, final int i) {
  2621.             return v1ij[i][j];
  2622.         }

  2624.         /**
  2625.          * Get the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0).
  2626.          *
  2627.          * @param j j index
  2628.          * @return the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0)
  2629.          */
  2630.         public double getU2(final int j) {
  2631.             return u2ij[j];
  2632.         }

  2634.         /**
  2635.          * Get the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0).
  2636.          *
  2637.          * @param j j index
  2638.          * @return the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0)
  2639.          */
  2640.         public double getV2(final int j) {
  2641.             return v2ij[j];
  2642.         }
  2643.     }

  2645.     /**
  2646.      * The U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients
  2647.      * described by equations 2.5.3-(21) and 2.5.3-(22) from Danielson.
  2648.      * <p>
  2649.      * The index i takes only the values 1 and 2<br>
  2650.      * For U only the index 0 for j is used.
  2651.      * </p>
  2652.      *
  2653.      * @author Petre Bazavan
  2654.      * @author Lucian Barbulescu
  2655.      * @param <T> type of the field elements
  2656.      */
  2657.     protected static class FieldUijVijCoefficients<T extends CalculusFieldElement<T>> {

  2659.         /**
  2660.          * The U₁<sup>j</sup> coefficients.
  2661.          * <p>
  2662.          * The first index identifies the Fourier coefficients used<br>
  2663.          * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and
  2664.          * S<sub>i</sub><sup>j</sup><br>
  2665.          * The only exception is when j = 0 when only the coefficient for fourier index
  2666.          * = 1 (i == 0) is needed.<br>
  2667.          * Also, for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are
  2668.          * computed, because are required to compute the coefficients U₂<sup>j</sup>
  2669.          * </p>
  2670.          */
  2671.         private final T[][] u1ij;

  2673.         /**
  2674.          * The V₁<sup>j</sup> coefficients.
  2675.          * <p>
  2676.          * The first index identifies the Fourier coefficients used<br>
  2677.          * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and
  2678.          * S<sub>i</sub><sup>j</sup><br>
  2679.          * for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed,
  2680.          * because are required to compute the coefficients V₂<sup>j</sup>
  2681.          * </p>
  2682.          */
  2683.         private final T[][] v1ij;

  2685.         /**
  2686.          * The U₂<sup>j</sup> coefficients.
  2687.          * <p>
  2688.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as
  2689.          * they are the only ones required.
  2690.          * </p>
  2691.          */
  2692.         private final T[] u2ij;

  2694.         /**
  2695.          * The V₂<sup>j</sup> coefficients.
  2696.          * <p>
  2697.          * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as
  2698.          * they are the only ones required.
  2699.          * </p>
  2700.          */
  2701.         private final T[] v2ij;

  2703.         /**
  2704.          * The current computed values for the ρ<sub>j</sub> and σ<sub>j</sub>
  2705.          * coefficients.
  2706.          */
  2707.         private final T[][] currentRhoSigmaj;

  2709.         /**
  2710.          * The C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> Fourier
  2711.          * coefficients.
  2712.          */
  2713.         private final FieldFourierCjSjCoefficients<T> fourierCjSj;

  2715.         /** The maximum value for j index. */
  2716.         private final int jMax;

  2718.         /**
  2719.          * Constructor.
  2720.          * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and
  2721.          *                         σ<sub>j</sub> coefficients
  2722.          * @param fourierCjSj      the fourier coefficients C<sub>i</sub><sup>j</sup>
  2723.          *                         and the S<sub>i</sub><sup>j</sup>
  2724.          * @param jMax             maximum value for j index
  2725.          * @param field            field used by default
  2726.          */
  2727.         FieldUijVijCoefficients(final T[][] currentRhoSigmaj, final FieldFourierCjSjCoefficients<T> fourierCjSj,
  2728.                 final int jMax, final Field<T> field) {
  2729.             this.currentRhoSigmaj = currentRhoSigmaj;
  2730.             this.fourierCjSj = fourierCjSj;
  2731.             this.jMax = jMax;

  2733.             // initialize the internal arrays.
  2734.             this.u1ij = MathArrays.buildArray(field, 6, 2 * jMax + 1);
  2735.             this.v1ij = MathArrays.buildArray(field, 6, 2 * jMax + 1);
  2736.             this.u2ij = MathArrays.buildArray(field, jMax + 1);
  2737.             this.v2ij = MathArrays.buildArray(field, jMax + 1);

  2739.             // compute the coefficients
  2740.             computeU1V1Coefficients(field);
  2741.             computeU2V2Coefficients();
  2742.         }

  2744.         /**
  2745.          * Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
  2746.          * @param field field used by default
  2747.          */
  2748.         private void computeU1V1Coefficients(final Field<T> field) {
  2749.             // Zero
  2750.             final T zero = field.getZero();

  2752.             // generate the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients
  2753.             // for j >= 1
  2754.             // also the U₁⁰ for Fourier index = 1 (i == 0) coefficient will be computed
  2755.             u1ij[0][0] = zero;
  2756.             for (int j = 1; j <= jMax; j++) {
  2757.                 // compute 1 / j
  2758.                 final double ooj = 1. / j;

  2760.                 for (int i = 0; i < 6; i++) {
  2761.                     // j is aready between 1 and J
  2762.                     u1ij[i][j] = fourierCjSj.getSij(i, j);
  2763.                     v1ij[i][j] = fourierCjSj.getCij(i, j);

  2765.                     // 1 - δ<sub>1j</sub> is 1 for all j > 1
  2766.                     if (j > 1) {
  2767.                         // k starts with 1 because j-J is less than or equal to 0
  2768.                         for (int kIndex = 1; kIndex <= j - 1; kIndex++) {
  2769.                             // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
  2770.                             // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
  2771.                             u1ij[i][j] = u1ij[i][j]
  2772.                                     .add(fourierCjSj.getCij(i, j - kIndex).multiply(currentRhoSigmaj[1][kIndex]).add(
  2773.                                             fourierCjSj.getSij(i, j - kIndex).multiply(currentRhoSigmaj[0][kIndex])));

  2775.                             // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  2776.                             // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  2777.                             v1ij[i][j] = v1ij[i][j].add(
  2778.                                     fourierCjSj.getCij(i, j - kIndex).multiply(currentRhoSigmaj[0][kIndex]).subtract(
  2779.                                             fourierCjSj.getSij(i, j - kIndex).multiply(currentRhoSigmaj[1][kIndex])));
  2780.                         }
  2781.                     }

  2783.                     // since j must be between 1 and J-1 and is already between 1 and J
  2784.                     // the following sum is skiped only for j = jMax
  2785.                     if (j != jMax) {
  2786.                         for (int kIndex = 1; kIndex <= jMax - j; kIndex++) {
  2787.                             // -C<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub> +
  2788.                             // S<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub>
  2789.                             u1ij[i][j] = u1ij[i][j].add(fourierCjSj.getCij(i, j + kIndex).negate()
  2790.                                     .multiply(currentRhoSigmaj[1][kIndex])
  2791.                                     .add(fourierCjSj.getSij(i, j + kIndex).multiply(currentRhoSigmaj[0][kIndex])));

  2793.                             // C<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub> +
  2794.                             // S<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub>
  2795.                             v1ij[i][j] = v1ij[i][j]
  2796.                                     .add(fourierCjSj.getCij(i, j + kIndex).multiply(currentRhoSigmaj[0][kIndex]).add(
  2797.                                             fourierCjSj.getSij(i, j + kIndex).multiply(currentRhoSigmaj[1][kIndex])));
  2798.                         }
  2799.                     }

  2801.                     for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  2802.                         // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  2803.                         // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  2804.                         u1ij[i][j] = u1ij[i][j].add(fourierCjSj.getCij(i, kIndex).negate()
  2805.                                 .multiply(currentRhoSigmaj[1][j + kIndex])
  2806.                                 .subtract(fourierCjSj.getSij(i, kIndex).multiply(currentRhoSigmaj[0][j + kIndex])));

  2808.                         // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  2809.                         // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  2810.                         v1ij[i][j] = v1ij[i][j]
  2811.                                 .add(fourierCjSj.getCij(i, kIndex).multiply(currentRhoSigmaj[0][j + kIndex])
  2812.                                         .add(fourierCjSj.getSij(i, kIndex).multiply(currentRhoSigmaj[1][j + kIndex])));
  2813.                     }

  2815.                     // divide by 1 / j
  2816.                     u1ij[i][j] = u1ij[i][j].multiply(-ooj);
  2817.                     v1ij[i][j] = v1ij[i][j].multiply(ooj);

  2819.                     // if index = 1 (i == 0) add the computed terms to U₁⁰
  2820.                     if (i == 0) {
  2821.                         // - (U₁<sup>j</sup> * ρ<sub>j</sub> + V₁<sup>j</sup> * σ<sub>j</sub>
  2822.                         u1ij[0][0] = u1ij[0][0].add(u1ij[0][j].negate().multiply(currentRhoSigmaj[0][j])
  2823.                                 .subtract(v1ij[0][j].multiply(currentRhoSigmaj[1][j])));
  2824.                     }
  2825.                 }
  2826.             }

  2828.             // Terms with j > jMax are required only when computing the coefficients
  2829.             // U₂<sup>j</sup> and V₂<sup>j</sup>
  2830.             // and those coefficients are only required for Fourier index = 1 (i == 0).
  2831.             for (int j = jMax + 1; j <= 2 * jMax; j++) {
  2832.                 // compute 1 / j
  2833.                 final double ooj = 1. / j;
  2834.                 // the value of i is 0
  2835.                 u1ij[0][j] = zero;
  2836.                 v1ij[0][j] = zero;

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

  2845.                     // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
  2846.                     // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
  2847.                     v1ij[0][j] = v1ij[0][j].add(fourierCjSj.getCij(0, j - kIndex).multiply(currentRhoSigmaj[0][kIndex])
  2848.                             .subtract(fourierCjSj.getSij(0, j - kIndex).multiply(currentRhoSigmaj[1][kIndex])));
  2849.                 }
  2850.                 for (int kIndex = 1; kIndex <= jMax; kIndex++) {
  2851.                     // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
  2852.                     // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
  2853.                     u1ij[0][j] = u1ij[0][j]
  2854.                             .add(fourierCjSj.getCij(0, kIndex).negate().multiply(currentRhoSigmaj[1][j + kIndex])
  2855.                                     .subtract(fourierCjSj.getSij(0, kIndex).multiply(currentRhoSigmaj[0][j + kIndex])));

  2857.                     // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
  2858.                     // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
  2859.                     v1ij[0][j] = v1ij[0][j].add(fourierCjSj.getCij(0, kIndex).multiply(currentRhoSigmaj[0][j + kIndex])
  2860.                             .add(fourierCjSj.getSij(0, kIndex).multiply(currentRhoSigmaj[1][j + kIndex])));
  2861.                 }

  2863.                 // divide by 1 / j
  2864.                 u1ij[0][j] = u1ij[0][j].multiply(-ooj);
  2865.                 v1ij[0][j] = v1ij[0][j].multiply(ooj);
  2866.             }
  2867.         }

  2869.         /**
  2870.          * Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
  2871.          * <p>
  2872.          * Only the coefficients for Fourier index = 1 (i == 0) are required.
  2873.          * </p>
  2874.          */
  2875.         private void computeU2V2Coefficients() {
  2876.             for (int j = 1; j <= jMax; j++) {
  2877.                 // compute 1 / j
  2878.                 final double ooj = 1. / j;

  2880.                 // only the values for i == 0 are computed
  2881.                 u2ij[j] = v1ij[0][j];
  2882.                 v2ij[j] = u1ij[0][j];

  2884.                 // 1 - δ<sub>1j</sub> is 1 for all j > 1
  2885.                 if (j > 1) {
  2886.                     for (int l = 1; l <= j - 1; l++) {
  2887.                         // U₁<sup>j-l</sup> * σ<sub>l</sub> +
  2888.                         // V₁<sup>j-l</sup> * ρ<sub>l</sub>
  2889.                         u2ij[j] = u2ij[j].add(u1ij[0][j - l].multiply(currentRhoSigmaj[1][l])
  2890.                                 .add(v1ij[0][j - l].multiply(currentRhoSigmaj[0][l])));

  2892.                         // U₁<sup>j-l</sup> * ρ<sub>l</sub> -
  2893.                         // V₁<sup>j-l</sup> * σ<sub>l</sub>
  2894.                         v2ij[j] = v2ij[j].add(u1ij[0][j - l].multiply(currentRhoSigmaj[0][l])
  2895.                                 .subtract(v1ij[0][j - l].multiply(currentRhoSigmaj[1][l])));
  2896.                     }
  2897.                 }

  2899.                 for (int l = 1; l <= jMax; l++) {
  2900.                     // -U₁<sup>j+l</sup> * σ<sub>l</sub> +
  2901.                     // U₁<sup>l</sup> * σ<sub>j+l</sub> +
  2902.                     // V₁<sup>j+l</sup> * ρ<sub>l</sub> -
  2903.                     // V₁<sup>l</sup> * ρ<sub>j+l</sub>
  2904.                     u2ij[j] = u2ij[j].add(u1ij[0][j + l].negate().multiply(currentRhoSigmaj[1][l])
  2905.                             .add(u1ij[0][l].multiply(currentRhoSigmaj[1][j + l]))
  2906.                             .add(v1ij[0][j + l].multiply(currentRhoSigmaj[0][l]))
  2907.                             .subtract(v1ij[0][l].multiply(currentRhoSigmaj[0][j + l])));

  2909.                     // U₁<sup>j+l</sup> * ρ<sub>l</sub> +
  2910.                     // U₁<sup>l</sup> * ρ<sub>j+l</sub> +
  2911.                     // V₁<sup>j+l</sup> * σ<sub>l</sub> +
  2912.                     // V₁<sup>l</sup> * σ<sub>j+l</sub>
  2913.                     u2ij[j] = u2ij[j].add(u1ij[0][j + l].multiply(currentRhoSigmaj[0][l])
  2914.                             .add(u1ij[0][l].multiply(currentRhoSigmaj[0][j + l]))
  2915.                             .add(v1ij[0][j + l].multiply(currentRhoSigmaj[1][l]))
  2916.                             .add(v1ij[0][l].multiply(currentRhoSigmaj[1][j + l])));
  2917.                 }

  2919.                 // divide by 1 / j
  2920.                 u2ij[j] = u2ij[j].multiply(-ooj);
  2921.                 v2ij[j] = v2ij[j].multiply(ooj);
  2922.             }
  2923.         }

  2925.         /**
  2926.          * Get the coefficient U₁<sup>j</sup> for Fourier index i.
  2927.          *
  2928.          * @param j j index
  2929.          * @param i Fourier index (starts at 0)
  2930.          * @return the coefficient U₁<sup>j</sup> for the given Fourier index i
  2931.          */
  2932.         public T getU1(final int j, final int i) {
  2933.             return u1ij[i][j];
  2934.         }

  2936.         /**
  2937.          * Get the coefficient V₁<sup>j</sup> for Fourier index i.
  2938.          *
  2939.          * @param j j index
  2940.          * @param i Fourier index (starts at 0)
  2941.          * @return the coefficient V₁<sup>j</sup> for the given Fourier index i
  2942.          */
  2943.         public T getV1(final int j, final int i) {
  2944.             return v1ij[i][j];
  2945.         }

  2947.         /**
  2948.          * Get the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0).
  2949.          *
  2950.          * @param j j index
  2951.          * @return the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0)
  2952.          */
  2953.         public T getU2(final int j) {
  2954.             return u2ij[j];
  2955.         }

  2957.         /**
  2958.          * Get the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0).
  2959.          *
  2960.          * @param j j index
  2961.          * @return the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0)
  2962.          */
  2963.         public T getV2(final int j) {
  2964.             return v2ij[j];
  2965.         }
  2966.     }

  2968.     /** Coefficients valid for one time slot. */
  2969.     protected static class Slot {

  2971.         /**
  2972.          * The coefficients D<sub>i</sub><sup>j</sup>.
  2973.          * <p>
  2974.          * Only for j = 1 and j = 2 the coefficients are not 0. <br>
  2975.          * i corresponds to the equinoctial element, as follows: - i=0 for a <br/>
  2976.          * - i=1 for k <br/>
  2977.          * - i=2 for h <br/>
  2978.          * - i=3 for q <br/>
  2979.          * - i=4 for p <br/>
  2980.          * - i=5 for λ <br/>
  2981.          * </p>
  2982.          */
  2983.         private final ShortPeriodicsInterpolatedCoefficient[] dij;

  2985.         /**
  2986.          * The coefficients C<sub>i</sub><sup>j</sup>.
  2987.          * <p>
  2988.          * The index order is cij[j][i] <br/>
  2989.          * i corresponds to the equinoctial element, as follows: <br/>
  2990.          * - i=0 for a <br/>
  2991.          * - i=1 for k <br/>
  2992.          * - i=2 for h <br/>
  2993.          * - i=3 for q <br/>
  2994.          * - i=4 for p <br/>
  2995.          * - i=5 for λ <br/>
  2996.          * </p>
  2997.          */
  2998.         private final ShortPeriodicsInterpolatedCoefficient[] cij;

  3000.         /**
  3001.          * The coefficients S<sub>i</sub><sup>j</sup>.
  3002.          * <p>
  3003.          * The index order is sij[j][i] <br/>
  3004.          * i corresponds to the equinoctial element, as follows: <br/>
  3005.          * - i=0 for a <br/>
  3006.          * - i=1 for k <br/>
  3007.          * - i=2 for h <br/>
  3008.          * - i=3 for q <br/>
  3009.          * - i=4 for p <br/>
  3010.          * - i=5 for λ <br/>
  3011.          * </p>
  3012.          */
  3013.         private final ShortPeriodicsInterpolatedCoefficient[] sij;

  3015.         /**
  3016.          * Simple constructor.
  3017.          * @param jMax                maximum value for j index
  3018.          * @param interpolationPoints number of points used in the interpolation process
  3019.          */
  3020.         Slot(final int jMax, final int interpolationPoints) {

  3022.             dij = new ShortPeriodicsInterpolatedCoefficient[3];
  3023.             cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
  3024.             sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];

  3026.             // Initialize the C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup> and
  3027.             // D<sub>i</sub><sup>j</sup> coefficients
  3028.             for (int j = 0; j <= jMax; j++) {
  3029.                 cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3030.                 if (j > 0) {
  3031.                     sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3032.                 }
  3033.                 // Initialize only the non-zero D<sub>i</sub><sup>j</sup> coefficients
  3034.                 if (j == 1 || j == 2) {
  3035.                     dij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
  3036.                 }
  3037.             }

  3039.         }

  3041.     }

  3043.     /** Coefficients valid for one time slot.
  3044.      * @param <T> type of the field elements
  3045.      */
  3046.     protected static class FieldSlot<T extends CalculusFieldElement<T>> {

  3048.         /**
  3049.          * The coefficients D<sub>i</sub><sup>j</sup>.
  3050.          * <p>
  3051.          * Only for j = 1 and j = 2 the coefficients are not 0. <br>
  3052.          * i corresponds to the equinoctial element, as follows: - i=0 for a <br/>
  3053.          * - i=1 for k <br/>
  3054.          * - i=2 for h <br/>
  3055.          * - i=3 for q <br/>
  3056.          * - i=4 for p <br/>
  3057.          * - i=5 for λ <br/>
  3058.          * </p>
  3059.          */
  3060.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] dij;

  3062.         /**
  3063.          * The coefficients C<sub>i</sub><sup>j</sup>.
  3064.          * <p>
  3065.          * The index order is cij[j][i] <br/>
  3066.          * i corresponds to the equinoctial element, as follows: <br/>
  3067.          * - i=0 for a <br/>
  3068.          * - i=1 for k <br/>
  3069.          * - i=2 for h <br/>
  3070.          * - i=3 for q <br/>
  3071.          * - i=4 for p <br/>
  3072.          * - i=5 for λ <br/>
  3073.          * </p>
  3074.          */
  3075.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] cij;

  3077.         /**
  3078.          * The coefficients S<sub>i</sub><sup>j</sup>.
  3079.          * <p>
  3080.          * The index order is sij[j][i] <br/>
  3081.          * i corresponds to the equinoctial element, as follows: <br/>
  3082.          * - i=0 for a <br/>
  3083.          * - i=1 for k <br/>
  3084.          * - i=2 for h <br/>
  3085.          * - i=3 for q <br/>
  3086.          * - i=4 for p <br/>
  3087.          * - i=5 for λ <br/>
  3088.          * </p>
  3089.          */
  3090.         private final FieldShortPeriodicsInterpolatedCoefficient<T>[] sij;

  3092.         /**
  3093.          * Simple constructor.
  3094.          * @param jMax                maximum value for j index
  3095.          * @param interpolationPoints number of points used in the interpolation process
  3096.          */
  3097.         @SuppressWarnings("unchecked")
  3098.         FieldSlot(final int jMax, final int interpolationPoints) {

  3100.             dij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array
  3101.                     .newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, 3);
  3102.             cij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array
  3103.                     .newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);
  3104.             sij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array
  3105.                     .newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);

  3107.             // Initialize the C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup> and
  3108.             // D<sub>i</sub><sup>j</sup> coefficients
  3109.             for (int j = 0; j <= jMax; j++) {
  3110.                 cij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3111.                 if (j > 0) {
  3112.                     sij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3113.                 }
  3114.                 // Initialize only the non-zero D<sub>i</sub><sup>j</sup> coefficients
  3115.                 if (j == 1 || j == 2) {
  3116.                     dij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
  3117.                 }
  3118.             }

  3120.         }

  3122.     }

  3124. }
Created with JaCoCo 0.8.12.202403310830