AbstractGaussianContribution.java

/* Copyright 2002-2015 CS Systèmes d'Information
 * Licensed to CS Systèmes d'Information (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * CS licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.propagation.semianalytical.dsst.forces;

import org.apache.commons.math3.analysis.UnivariateVectorFunction;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.apache.commons.math3.util.FastMath;
import org.orekit.attitudes.Attitude;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitExceptionWrapper;
import org.orekit.forces.ForceModel;
import org.orekit.frames.Frame;
import org.orekit.orbits.EquinoctialOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.numerical.TimeDerivativesEquations;
import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
import org.orekit.time.AbsoluteDate;

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

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

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

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

    /** Maximum value for j index. */
    private static final int JMAX = 12;
    // CHECKSTYLE: stop VisibilityModifierCheck

    /** Retrograde factor. */
    protected int I;

    /** a. */
    protected double a;
    /** e<sub>x</sub>. */
    protected double k;
    /** e<sub>y</sub>. */
    protected double h;
    /** h<sub>x</sub>. */
    protected double q;
    /** h<sub>y</sub>. */
    protected double p;

    /** Eccentricity. */
    protected double ecc;

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

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

    /** Equinoctial frame f vector. */
    protected Vector3D f;
    /** Equinoctial frame g vector. */
    protected Vector3D g;
    /** Equinoctial frame w vector. */
    protected Vector3D w;

    /** A = sqrt(μ * a). */
    protected double A;
    /** B = sqrt(1 - h² - k²). */
    protected double B;
    /** C = 1 + p² + q². */
    protected double C;

    /** 2 / (n² * a) . */
    protected double ton2a;
    /** 1 / A .*/
    protected double ooA;
    /** 1 / (A * B) .*/
    protected double ooAB;
    /** C / (2 * A * B) .*/
    protected double co2AB;
    /** 1 / (1 + B) .*/
    protected double ooBpo;
    /** 1 / μ .*/
    protected double ooMu;
    /** μ .*/
    protected double mu;

    // CHECKSTYLE: resume VisibilityModifierCheck

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

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

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

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

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

    /** The C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients used to compute
     * the short-periodic gaussian contribution. */
    private GaussianShortPeriodicCoefficients gaussianSPCoefs;

    /** The frame used to describe the orbits. */
    private Frame frame;

    /** Build a new instance.
     *
     *  @param threshold tolerance for the choice of the Gauss quadrature order
     *  @param contribution the {@link ForceModel} to be numerically averaged
     */
    protected AbstractGaussianContribution(final double threshold,
            final ForceModel contribution) {
        this.contribution = contribution;
        this.threshold  = threshold;
        this.integrator = new GaussQuadrature(GAUSS_ORDER[MAX_ORDER_RANK]);
        this.isDirty    = true;
    }

    /** {@inheritDoc} */
    @Override
    public void initialize(final AuxiliaryElements aux, final boolean meanOnly)
        throws OrekitException {

        // save the frame
        this.frame = aux.getFrame();

        if (!meanOnly) {
            // initialize the short periodic coefficient generator, if needed.
            this.gaussianSPCoefs = new GaussianShortPeriodicCoefficients(JMAX, INTERPOLATION_POINTS);
        }
    }

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

        // Equinoctial elements
        a  = aux.getSma();
        k  = aux.getK();
        h  = aux.getH();
        q  = aux.getQ();
        p  = aux.getP();

        // Retrograde factor
        I = aux.getRetrogradeFactor();

        // Eccentricity
        ecc = aux.getEcc();

        // Equinoctial coefficients
        A = aux.getA();
        B = aux.getB();
        C = aux.getC();

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

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

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

        // 1 / A
        ooA = 1. / A;
        // 1 / AB
        ooAB = ooA / B;
        // C / 2AB
        co2AB = C * ooAB / 2.;
        // 1 / (1 + B)
        ooBpo = 1. / (1. + B);
        // 2 / (n² * a)
        ton2a = 2. / (n * n * a);
        // mu
        mu = aux.getMu();
        // 1 / mu
        ooMu  = 1. / mu;
    }

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

        double[] meanElementRate = new double[6];
        // Computes the limits for the integral
        final double[] ll = getLLimits(state);
        // Computes integrated mean element rates if Llow < Lhigh
        if (ll[0] < ll[1]) {
            meanElementRate = getMeanElementRate(state, integrator, ll[0], ll[1]);
            if (isDirty) {
                boolean next = true;
                for (int i = 0; i < MAX_ORDER_RANK && next; i++) {
                    final double[] meanRates = getMeanElementRate(state, new GaussQuadrature(GAUSS_ORDER[i]), ll[0], ll[1]);
                    if (getRatesDiff(meanElementRate, meanRates) < threshold) {
                        integrator = new GaussQuadrature(GAUSS_ORDER[i]);
                        next = false;
                    }
                }
                isDirty = false;
            }
        }
        return meanElementRate;
    }

    /** Compute the acceleration due to the non conservative perturbing force.
     *
     *  @param state current state information: date, kinematics, attitude
     *  @return the perturbing acceleration
     *  @exception OrekitException if some specific error occurs
     */
    protected Vector3D getAcceleration(final SpacecraftState state)
        throws OrekitException {
        final AccelerationRetriever retriever = new AccelerationRetriever(state);
        contribution.addContribution(state, retriever);

        return retriever.getAcceleration();
    }

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

    /** Computes the mean equinoctial elements rates da<sub>i</sub> / dt.
     *
     *  @param state current state
     *  @param gauss Gauss quadrature
     *  @param low lower bound of the integral interval
     *  @param high upper bound of the integral interval
     *  @return the mean element rates
     *  @throws OrekitException if some specific error occurs
     */
    private double[] getMeanElementRate(final SpacecraftState state,
            final GaussQuadrature gauss,
            final double low,
            final double high) throws OrekitException {
        final double[] meanElementRate = gauss.integrate(new IntegrableFunction(state, true, 0), low, high);
        // Constant multiplier for integral
        final double coef = 1. / (2. * FastMath.PI * B);
        // Corrects mean element rates
        for (int i = 0; i < 6; i++) {
            meanElementRate[i] *= coef;
        }
        return meanElementRate;
    }

    /** Estimates the weighted magnitude of the difference between 2 sets of equinoctial elements rates.
     *
     *  @param meanRef reference rates
     *  @param meanCur current rates
     *  @return estimated magnitude of weighted differences
     */
    private double getRatesDiff(final double[] meanRef, final double[] meanCur) {
        double maxDiff = FastMath.abs(meanRef[0] - meanCur[0]) / a;
        // Corrects mean element rates
        for (int i = 1; i < meanRef.length; i++) {
            final double diff = FastMath.abs(meanRef[i] - meanCur[i]);
            if (maxDiff < diff) maxDiff = diff;
        }
        return maxDiff;
    }

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

    /** {@inheritDoc} */
    @Override
    public double[] getShortPeriodicVariations(final AbsoluteDate date,
            final double[] meanElements) throws OrekitException {

        // Build an Orbit object from the mean elements
        final Orbit meanOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(
                meanElements, PositionAngle.MEAN, date, this.mu, this.frame);

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

        // Compute the center (l - λ)
        final double center =  L - meanElements[5];
        // Compute (l - λ)²
        final double center2 = center * center;

        // Initialize short periodic variations
        final double[] shortPeriodicVariation = new double[6];
        for (int i = 0; i < 6; i++) {
            shortPeriodicVariation[i] = gaussianSPCoefs.getCij(i, 0, date) +
                    center * gaussianSPCoefs.getDij(i, 1, date);
            if (i == 5) {
                shortPeriodicVariation[i] += center2 * gaussianSPCoefs.getDij(i, 2, date);
            }
        }

        for (int j = 1; j <= JMAX; j++) {
            for (int i = 0; i < 6; i++) {
                // Get Cij and Sij
                final double cij = gaussianSPCoefs.getCij(i, j, date);
                final double sij = gaussianSPCoefs.getSij(i, j, date);

                // add corresponding term to the short periodic variation
                shortPeriodicVariation[i] += cij * FastMath.cos(j * L);
                shortPeriodicVariation[i] += sij * FastMath.sin(j * L);
            }
        }

        return shortPeriodicVariation;

    }

    /** {@inheritDoc} */
    @Override
    public void computeShortPeriodicsCoefficients(final SpacecraftState state)
        throws OrekitException {

        //Compute the coefficients
        gaussianSPCoefs.computeCoefficients(state);
    }

    /** {@inheritDoc} */
    @Override
    public void resetShortPeriodicsCoefficients() {
        if (gaussianSPCoefs != null) {
            // reset the coefficients
            gaussianSPCoefs.resetCoefficients();
        }
    }

    /** Internal class for retrieving acceleration from a {@link ForceModel}. */
    private static class AccelerationRetriever implements TimeDerivativesEquations {

        /** acceleration vector. */
        private Vector3D acceleration;

        /** state. */
        private final SpacecraftState state;

        /** Simple constructor.
         *  @param state input state
         */
        public AccelerationRetriever(final SpacecraftState state) {
            this.acceleration = Vector3D.ZERO;
            this.state = state;
        }

        /** {@inheritDoc} */
        @Override
        public void addKeplerContribution(final double mu) {
        }

        /** {@inheritDoc} */
        @Override
        public void addXYZAcceleration(final double x, final double y,
                final double z) {
            //TODO How to be sure we are in the good frame ???
            acceleration = new Vector3D(x, y, z);
        }

        /** {@inheritDoc} */
        @Override
        public void addAcceleration(final Vector3D gamma, final Frame frame)
            throws OrekitException {
            acceleration = frame.getTransformTo(state.getFrame(),
                    state.getDate()).transformVector(gamma);
        }

        /** {@inheritDoc} */
        @Override
        public void addMassDerivative(final double q) {
        }

        /** Get the acceleration vector.
         * @return acceleration vector
         */
        public Vector3D getAcceleration() {
            return acceleration;
        }

    }

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

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

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

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

        /** Build a new instance.
         *  @param  state current state information: date, kinematics, attitude
         *  @param meanMode if true return the value associated to the mean elements variation,
         *                  if false return the values associated to the short periodic elements variation
         * @param j the j index. used only for short periodic variation. Ignored for mean elements variation.
         */
        public IntegrableFunction(final SpacecraftState state, final boolean meanMode, final int j) {
            this.state = state;
            this.meanMode = meanMode;
            this.j = j;
        }

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

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

            final double cosL = FastMath.cos(x);
            final double sinL = FastMath.sin(x);
            final double roa  = B * B / (1. + h * sinL + k * cosL);
            final double roa2 = roa * roa;
            final double r    = a * roa;
            final double X    = r * cosL;
            final double Y    = r * sinL;
            final double naob = n * a / B;
            final double Xdot = -naob * (h + sinL);
            final double Ydot =  naob * (k + cosL);
            final Vector3D vel = new Vector3D(Xdot, f, Ydot, g);

            // Compute acceleration
            Vector3D acc = Vector3D.ZERO;
            try {

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

                // Recompose an orbit with time held fixed to be compliant with DSST theory
                final Orbit recomposedOrbit =
                        new EquinoctialOrbit(shiftedOrbit.getA(),
                                             shiftedOrbit.getEquinoctialEx(),
                                             shiftedOrbit.getEquinoctialEy(),
                                             shiftedOrbit.getHx(),
                                             shiftedOrbit.getHy(),
                                             shiftedOrbit.getLv(),
                                             PositionAngle.TRUE,
                                             shiftedOrbit.getFrame(),
                                             state.getDate(),
                                             shiftedOrbit.getMu());

                // Get the corresponding attitude
                final Attitude recomposedAttitude =
                        attitudeProvider.getAttitude(recomposedOrbit,
                                                     recomposedOrbit.getDate(),
                                                     recomposedOrbit.getFrame());

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

                acc = getAcceleration(shiftedState);

            } catch (OrekitException oe) {
                throw new OrekitExceptionWrapper(oe);
            }
            //Compute the derivatives of the elements by the speed
            final double[] deriv = new double[6];
            // da/dv
            deriv[0] = getAoV(vel).dotProduct(acc);
            // dex/dv
            deriv[1] = getKoV(X, Y, Xdot, Ydot).dotProduct(acc);
            // dey/dv
            deriv[2] = getHoV(X, Y, Xdot, Ydot).dotProduct(acc);
            // dhx/dv
            deriv[3] = getQoV(X).dotProduct(acc);
            // dhy/dv
            deriv[4] = getPoV(Y).dotProduct(acc);
            // dλ/dv
            deriv[5] = getLoV(X, Y, Xdot, Ydot).dotProduct(acc);

            // Compute mean elements rates
            double[] val = null;
            if (meanMode) {
                val = new double[6];
                for (int i = 0; i < 6; i++) {
                    // da<sub>i</sub>/dt
                    val[i] = roa2 * deriv[i];
                }
            } else {
                val = new double [12];
                //Compute cos(j*L) and sin(j*L);
                final double cosjL = j == 1 ? cosL : FastMath.cos(j * x);
                final double sinjL = j == 1 ? sinL : FastMath.sin(j * x);

                for (int i = 0; i < 6; i++) {
                    // da<sub>i</sub>/dv * cos(jL)
                    val[i] = cosjL * deriv[i];
                    // da<sub>i</sub>/dv * sin(jL)
                    val[i + 6] = sinjL * deriv[i];
                }
            }
            return val;
        }

        /** Converts true longitude to eccentric longitude.
         * @param lv True longitude
         * @return Eccentric longitude
         */
        private double trueToEccentric (final double lv) {
            final double cosLv   = FastMath.cos(lv);
            final double sinLv   = FastMath.sin(lv);
            final double num     = h * cosLv - k * sinLv;
            final double den     = B + 1 + k * cosLv + h * sinLv;
            return lv + 2 * FastMath.atan(num / den);
        }

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

        /** Converts true longitude to mean longitude.
         * @param lv True longitude
         * @return Eccentric longitude
         */
        private double trueToMean (final double lv) {
            return eccentricToMean(trueToEccentric(lv));
        }

        /** Compute δa/δv.
         *  @param vel satellite velocity
         *  @return δa/δv
         */
        private Vector3D getAoV(final Vector3D vel) {
            return new Vector3D(ton2a, vel);
        }

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

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

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

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

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

    }

    /** Class used to {@link #integrate(UnivariateVectorFunction, double, double) integrate}
     *  a {@link org.apache.commons.math3.analysis.UnivariateVectorFunction function}
     *  of the orbital elements using the Gaussian quadrature rule to get the acceleration.
     */
    private static class GaussQuadrature {

        // CHECKSTYLE: stop NoWhitespaceAfter

        // Points and weights for the available quadrature orders

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

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

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

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

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

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

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

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

        /** Points for quadrature of order 32. */
        private static final double[] P_32 = {
            -0.99726386184948160000,
            -0.98561151154526840000,
            -0.96476225558750640000,
            -0.93490607593773970000,
            -0.89632115576605220000,
            -0.84936761373256990000,
            -0.79448379596794250000,
            -0.73218211874028970000,
            -0.66304426693021520000,
            -0.58771575724076230000,
            -0.50689990893222950000,
            -0.42135127613063540000,
            -0.33186860228212767000,
            -0.23928736225213710000,
            -0.14447196158279646000,
            -0.04830766568773831000,
            0.04830766568773831000,
            0.14447196158279646000,
            0.23928736225213710000,
            0.33186860228212767000,
            0.42135127613063540000,
            0.50689990893222950000,
            0.58771575724076230000,
            0.66304426693021520000,
            0.73218211874028970000,
            0.79448379596794250000,
            0.84936761373256990000,
            0.89632115576605220000,
            0.93490607593773970000,
            0.96476225558750640000,
            0.98561151154526840000,
            0.99726386184948160000
        };

        /** Weights for quadrature of order 32. */
        private static final double[] W_32 = {
            0.00701861000947013600,
            0.01627439473090571200,
            0.02539206530926214200,
            0.03427386291302141000,
            0.04283589802222658600,
            0.05099805926237621600,
            0.05868409347853559000,
            0.06582222277636193000,
            0.07234579410884862000,
            0.07819389578707042000,
            0.08331192422694673000,
            0.08765209300440380000,
            0.09117387869576390000,
            0.09384439908080441000,
            0.09563872007927487000,
            0.09654008851472784000,
            0.09654008851472784000,
            0.09563872007927487000,
            0.09384439908080441000,
            0.09117387869576390000,
            0.08765209300440380000,
            0.08331192422694673000,
            0.07819389578707042000,
            0.07234579410884862000,
            0.06582222277636193000,
            0.05868409347853559000,
            0.05099805926237621600,
            0.04283589802222658600,
            0.03427386291302141000,
            0.02539206530926214200,
            0.01627439473090571200,
            0.00701861000947013600
        };

        /** Points for quadrature of order 40. */
        private static final double[] P_40 = {
            -0.99823770971055930000,
            -0.99072623869945710000,
            -0.97725994998377420000,
            -0.95791681921379170000,
            -0.93281280827867660000,
            -0.90209880696887420000,
            -0.86595950321225960000,
            -0.82461223083331170000,
            -0.77830565142651940000,
            -0.72731825518992710000,
            -0.67195668461417960000,
            -0.61255388966798030000,
            -0.54946712509512820000,
            -0.48307580168617870000,
            -0.41377920437160500000,
            -0.34199409082575850000,
            -0.26815218500725370000,
            -0.19269758070137110000,
            -0.11608407067525522000,
            -0.03877241750605081600,
            0.03877241750605081600,
            0.11608407067525522000,
            0.19269758070137110000,
            0.26815218500725370000,
            0.34199409082575850000,
            0.41377920437160500000,
            0.48307580168617870000,
            0.54946712509512820000,
            0.61255388966798030000,
            0.67195668461417960000,
            0.72731825518992710000,
            0.77830565142651940000,
            0.82461223083331170000,
            0.86595950321225960000,
            0.90209880696887420000,
            0.93281280827867660000,
            0.95791681921379170000,
            0.97725994998377420000,
            0.99072623869945710000,
            0.99823770971055930000
        };

        /** Weights for quadrature of order 40. */
        private static final double[] W_40 = {
            0.00452127709853309800,
            0.01049828453115270400,
            0.01642105838190797300,
            0.02224584919416689000,
            0.02793700698002338000,
            0.03346019528254786500,
            0.03878216797447199000,
            0.04387090818567333000,
            0.04869580763507221000,
            0.05322784698393679000,
            0.05743976909939157000,
            0.06130624249292891000,
            0.06480401345660108000,
            0.06791204581523394000,
            0.07061164739128681000,
            0.07288658239580408000,
            0.07472316905796833000,
            0.07611036190062619000,
            0.07703981816424793000,
            0.07750594797842482000,
            0.07750594797842482000,
            0.07703981816424793000,
            0.07611036190062619000,
            0.07472316905796833000,
            0.07288658239580408000,
            0.07061164739128681000,
            0.06791204581523394000,
            0.06480401345660108000,
            0.06130624249292891000,
            0.05743976909939157000,
            0.05322784698393679000,
            0.04869580763507221000,
            0.04387090818567333000,
            0.03878216797447199000,
            0.03346019528254786500,
            0.02793700698002338000,
            0.02224584919416689000,
            0.01642105838190797300,
            0.01049828453115270400,
            0.00452127709853309800
        };

        /** Points for quadrature of order 48. */
        private static final double[] P_48 = {
            -0.99877100725242610000,
            -0.99353017226635080000,
            -0.98412458372282700000,
            -0.97059159254624720000,
            -0.95298770316043080000,
            -0.93138669070655440000,
            -0.90587913671556960000,
            -0.87657202027424800000,
            -0.84358826162439350000,
            -0.80706620402944250000,
            -0.76715903251574020000,
            -0.72403413092381470000,
            -0.67787237963266400000,
            -0.62886739677651370000,
            -0.57722472608397270000,
            -0.52316097472223300000,
            -0.46690290475095840000,
            -0.40868648199071680000,
            -0.34875588629216070000,
            -0.28736248735545555000,
            -0.22476379039468908000,
            -0.16122235606889174000,
            -0.09700469920946270000,
            -0.03238017096286937000,
            0.03238017096286937000,
            0.09700469920946270000,
            0.16122235606889174000,
            0.22476379039468908000,
            0.28736248735545555000,
            0.34875588629216070000,
            0.40868648199071680000,
            0.46690290475095840000,
            0.52316097472223300000,
            0.57722472608397270000,
            0.62886739677651370000,
            0.67787237963266400000,
            0.72403413092381470000,
            0.76715903251574020000,
            0.80706620402944250000,
            0.84358826162439350000,
            0.87657202027424800000,
            0.90587913671556960000,
            0.93138669070655440000,
            0.95298770316043080000,
            0.97059159254624720000,
            0.98412458372282700000,
            0.99353017226635080000,
            0.99877100725242610000
        };

        /** Weights for quadrature of order 48. */
        private static final double[] W_48 = {
            0.00315334605230596250,
            0.00732755390127620800,
            0.01147723457923446900,
            0.01557931572294386600,
            0.01961616045735556700,
            0.02357076083932435600,
            0.02742650970835688000,
            0.03116722783279807000,
            0.03477722256477045000,
            0.03824135106583080600,
            0.04154508294346483000,
            0.04467456085669424000,
            0.04761665849249054000,
            0.05035903555385448000,
            0.05289018948519365000,
            0.05519950369998416500,
            0.05727729210040315000,
            0.05911483969839566000,
            0.06070443916589384000,
            0.06203942315989268000,
            0.06311419228625403000,
            0.06392423858464817000,
            0.06446616443595010000,
            0.06473769681268386000,
            0.06473769681268386000,
            0.06446616443595010000,
            0.06392423858464817000,
            0.06311419228625403000,
            0.06203942315989268000,
            0.06070443916589384000,
            0.05911483969839566000,
            0.05727729210040315000,
            0.05519950369998416500,
            0.05289018948519365000,
            0.05035903555385448000,
            0.04761665849249054000,
            0.04467456085669424000,
            0.04154508294346483000,
            0.03824135106583080600,
            0.03477722256477045000,
            0.03116722783279807000,
            0.02742650970835688000,
            0.02357076083932435600,
            0.01961616045735556700,
            0.01557931572294386600,
            0.01147723457923446900,
            0.00732755390127620800,
            0.00315334605230596250
        };
        // CHECKSTYLE: resume NoWhitespaceAfter

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

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

        /** Creates a Gauss integrator of the given order.
         *
         *  @param numberOfPoints Order of the integration rule.
         */
        public GaussQuadrature(final int numberOfPoints) {

            switch(numberOfPoints) {
            case 12 :
                this.nodePoints  = P_12.clone();
                this.nodeWeights = W_12.clone();
                break;
            case 16 :
                this.nodePoints  = P_16.clone();
                this.nodeWeights = W_16.clone();
                break;
            case 20 :
                this.nodePoints  = P_20.clone();
                this.nodeWeights = W_20.clone();
                break;
            case 24 :
                this.nodePoints  = P_24.clone();
                this.nodeWeights = W_24.clone();
                break;
            case 32 :
                this.nodePoints  = P_32.clone();
                this.nodeWeights = W_32.clone();
                break;
            case 40 :
                this.nodePoints  = P_40.clone();
                this.nodeWeights = W_40.clone();
                break;
            case 48 :
            default :
                this.nodePoints  = P_48.clone();
                this.nodeWeights = W_48.clone();
                break;
            }

        }

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

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

        /** Performs a change of variable so that the integration
         *  can be performed on an arbitrary interval {@code [a, b]}.
         *  <p>
         *  It is assumed that the natural interval is {@code [-1, 1]}.
         *  </p>
         *
         * @param points  Points to adapt to the new interval.
         * @param weights Weights to adapt to the new interval.
         * @param a Lower bound of the integration interval.
         * @param b Lower bound of the integration interval.
         */
        private void transform(final double[] points, final double[] weights,
                final double a, final double b) {
            // Scaling
            final double scale = (b - a) / 2;
            final double shift = a + scale;
            for (int i = 0; i < points.length; i++) {
                points[i]   = points[i] * scale + shift;
                weights[i] *= scale;
            }
        }

        /** Returns an estimate of the integral of {@code f(x) * w(x)},
         *  where {@code w} is a weight function that depends on the actual
         *  flavor of the Gauss integration scheme.
         *
         * @param f Function to integrate.
         * @param points  Nodes.
         * @param weights Nodes weights.
         * @return the integral of the weighted function.
         */
        private double[] basicIntegrate(final UnivariateVectorFunction f,
                final double[] points,
                final double[] weights) {
            double x = points[0];
            double w = weights[0];
            double[] v = f.value(x);
            final double[] y = new double[v.length];
            for (int j = 0; j < v.length; j++) {
                y[j] = w * v[j];
            }
            final double[] t = y.clone();
            final double[] c = new double[v.length];
            final double[] s = t.clone();
            for (int i = 1; i < points.length; i++) {
                x = points[i];
                w = weights[i];
                v = f.value(x);
                for (int j = 0; j < v.length; j++) {
                    y[j] = w * v[j] - c[j];
                    t[j] =  s[j] + y[j];
                    c[j] = (t[j] - s[j]) - y[j];
                    s[j] = t[j];
                }
            }
            return s;
        }

    }

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

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

        /** The C<sub>i</sub><sup>j</sup> coefficients.
         * <p>
         * the index i corresponds to the following elements: <br/>
         * - 0 for a <br>
         * - 1 for k <br>
         * - 2 for h <br>
         * - 3 for q <br>
         * - 4 for p <br>
         * - 5 for λ <br>
         * </p>
         */
        private final double[][] cCoef;

        /** The C<sub>i</sub><sup>j</sup> coefficients.
         * <p>
         * the index i corresponds to the following elements: <br/>
         * - 0 for a <br>
         * - 1 for k <br>
         * - 2 for h <br>
         * - 3 for q <br>
         * - 4 for p <br>
         * - 5 for λ <br>
         * </p>
         */
        private final double[][] sCoef;

        /** Standard constructor.
         * @param state the current state
         * @param jMax maximum value for j
         * @throws OrekitException in case of an error
         */
        public FourierCjSjCoefficients(final SpacecraftState state, final int jMax) throws OrekitException {
            //Initialise the fields
            this.jMax = jMax;

            //Allocate the arrays
            cCoef = new double[jMax + 1][6];
            sCoef = new double[jMax + 1][6];

            //Compute the coefficients
            computeCoefficients(state);
        }

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

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

                    //divide by PI and set the values for the coefficients
                    for (int i = 0; i < 6; i++) {
                        cCoef[j][i] = ooPI * curentCoefficients[i];
                        sCoef[j][i] = ooPI * curentCoefficients[i + 6];
                    }
                }
            }
        }

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

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

    /** This class handles the short periodic coefficients described in Danielson 2.5.3-26.
     *
     * <p>
     * The value of M is 0. Also, since the values of the Fourier coefficient D<sub>i</sub><sup>m</sup> is 0
     * then the values of the coefficients D<sub>i</sub><sup>m</sup> for m > 2 are also 0.
     * </p>
     * @author Petre Bazavan
     * @author Lucian Barbulescu
     *
     */
    private class GaussianShortPeriodicCoefficients {

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

        /**The coefficients D<sub>i</sub><sup>j</sup>.
         * <p>
         * Only for j = 1 and j = 2 the coefficients are not 0. <br>
         * i corresponds to the equinoctial element, as follows:
         * - i=0 for a <br/>
         * - i=1 for k <br/>
         * - i=2 for h <br/>
         * - i=3 for q <br/>
         * - i=4 for p <br/>
         * - i=5 for λ <br/>
         * </p>
         */
        private final ShortPeriodicsInterpolatedCoefficient[][] dij;

        /** The coefficients C<sub>i</sub><sup>j</sup>.
         * <p>
         * The index order is cij[j][i] <br/>
         * i corresponds to the equinoctial element, as follows: <br/>
         * - i=0 for a <br/>
         * - i=1 for k <br/>
         * - i=2 for h <br/>
         * - i=3 for q <br/>
         * - i=4 for p <br/>
         * - i=5 for λ <br/>
         * </p>
         */
        private final ShortPeriodicsInterpolatedCoefficient[][] cij;

        /** The coefficients S<sub>i</sub><sup>j</sup>.
         * <p>
         * The index order is sij[j][i] <br/>
         * i corresponds to the equinoctial element, as follows: <br/>
         * - i=0 for a <br/>
         * - i=1 for k <br/>
         * - i=2 for h <br/>
         * - i=3 for q <br/>
         * - i=4 for p <br/>
         * - i=5 for λ <br/>
         * </p>
         */
        private final ShortPeriodicsInterpolatedCoefficient[][] sij;

        /** The current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients.
         * <p>
         * Index 0 corresponds to ρ, index 1 corresponds to σ
         * Used to compute the U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients.
         * </p>
         */
        private final double[][] currentRhoSigmaj;

        /** Constructor.
         *  @param jMax maximum value for j index
         *  @param interpolationPoints number of points used in the interpolation process
         */
        public GaussianShortPeriodicCoefficients(final int jMax, final int interpolationPoints) {
            //Initialise fields
            this.jMax = jMax;

            this.dij = new ShortPeriodicsInterpolatedCoefficient[3][6];

            this.cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1][6];
            this.sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1][6];

            this.currentRhoSigmaj = new double[2][3 * jMax + 1];

            // Initialise the C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup> and D<sub>i</sub><sup>j</sup> coefficients
            for (int j = 0; j <= jMax; j++) {
                for (int i = 0; i < 6; i++) {
                    this.cij[j][i] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
                    if (j > 0) {
                        this.sij[j][i] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
                    }
                    // Initialise only the non-zero D<sub>i</sub><sup>j</sup> coefficients
                    if (j == 1 || (j == 2 && i == 5)) {
                        this.dij[j][i] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
                    }
                }
            }
        }

        /** Compute the short periodic coefficients.
         *
         * @param state current state information: date, kinematics, attitude
         * @throws OrekitException if an error occurs
         */
        public void computeCoefficients(final SpacecraftState state)
            throws OrekitException {

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

            // Compute ρ<sub>j</sub> and σ<sub>j</sub>
            computeRhoSigmaCoefficients(date);

            // Compute the Fourier coefficients
            final FourierCjSjCoefficients fourierCjSj = new FourierCjSjCoefficients(state, jMax);

            // Compute the required U and V coefficients
            final UijVijCoefficients uijvij = new UijVijCoefficients(currentRhoSigmaj, fourierCjSj, jMax);

            // compute the k₂⁰ coefficient
            final double k20 = computeK20(jMax);

            // 1. / n
            final double oon = 1. / n;
            // 3. / (2 * a * n)
            final double to2an = 1.5 * oon / a;
            // 3. / (4 * a * n)
            final double to4an = to2an / 2;

            // Compute the coefficients for each element
            for (int i = 0; i < 6; i++) {

                // compute D<sub>i</sub>¹ and D<sub>i</sub>² (all others are 0)
                double di1 = -oon * fourierCjSj.getCij(i, 0);
                if (i == 5) {
                    di1 += to2an * uijvij.getU1(0, 0);
                }
                double di2 = 0.;
                if (i == 5) {
                    di2 += -to4an * fourierCjSj.getCij(0, 0);
                }

                //the C<sub>i</sub>⁰ is computed based on all others
                double currentCi0 = -di2 * k20;

                for (int j = 1; j <= jMax; j++) {
                    // compute the current C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup>
                    double currentCij = oon * uijvij.getU1(j, i);
                    if (i == 5) {
                        currentCij += -to2an * uijvij.getU2(j);
                    }
                    double currentSij = oon * uijvij.getV1(j, i);
                    if (i == 5) {
                        currentSij += -to2an * uijvij.getV2(j);
                    }

                    // add the computed coefficients to C<sub>i</sub>⁰
                    currentCi0 += -(currentCij * currentRhoSigmaj[0][j] + currentSij * currentRhoSigmaj[1][j]);

                    // add the values to the interpolators
                    cij[j][i].addGridPoint(date, currentCij);
                    sij[j][i].addGridPoint(date, currentSij);
                }

                // add the computed values to the interpolators
                cij[0][i].addGridPoint(date, currentCi0);
                dij[1][i].addGridPoint(date, di1);
                if (i == 5) {
                    dij[2][i].addGridPoint(date, di2);
                }
            }
        }
        /** Reset the coefficients.
         * <p>
         * For each coefficient, clear history of computed points
         * </p>
         */
        public void resetCoefficients() {

            for (int j = 0; j <= jMax; j++) {
                for (int i = 0; i < 6; i++) {
                    this.cij[j][i].clearHistory();
                    if (j > 0) {
                        this.sij[j][i].clearHistory();
                    }
                    if (j == 1 || (j == 2 && i == 5)) {
                        this.dij[j][i].clearHistory();
                    }
                }
            }
        }

        /**
         * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
         * <p>
         * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
         *  ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
         *  σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
         * </p>
         * @param date current date
         */
        private void computeRhoSigmaCoefficients(final AbsoluteDate date) {
            final CjSjCoefficient cjsjKH = new CjSjCoefficient(k, h);
            final double b = 1. / (1 + B);

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

            for (int j = 1; j <= 3 * jMax; j++) {

                //Compute current rho and sigma;
                mbtj *= -b;
                final double coef = (1 + j * B) * mbtj;
                currentRhoSigmaj[0][j] = coef * cjsjKH.getCj(j);
                currentRhoSigmaj[1][j] = coef * cjsjKH.getSj(j);
            }
        }

        /** Compute the coefficient k₂⁰ by using the equation
         * 2.5.3-(9a) from Danielson.
         * <p>
         * After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:<br>
         * k₂⁰ = &Sigma;<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
         * </p>
         * @param kMax max value fot k index
         * @return the coefficient k₂⁰
         */
        private double computeK20(final int kMax) {
            double k20 = 0.;

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

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

                // add the term to the result
                k20 += currentTerm;
            }

            return k20;
        }

        /** Get C<sub>i</sub><sup>j</sup>.
         *
         * @param i i index
         * @param j j index
         * @param date the date
         * @return C<sub>i</sub><sup>j</sup>
         */
        public double getCij(final int i, final int j, final AbsoluteDate date) {
            return cij[j][i].value(date);
        }

        /** Get S<sub>i</sub><sup>j</sup>.
         *
         * @param i i index
         * @param j j index
         * @param date the date
         * @return S<sub>i</sub><sup>j</sup>
         */
        public double getSij(final int i, final int j, final AbsoluteDate date) {
            return sij[j][i].value(date);
        }

        /** Get D<sub>i</sub><sup>j</sup>.
         * @param i i index
         * @param j j index
         * @param date target date
         * @return D<sub>i</sub><sup>j</sup>
         */
        public double getDij(final int i, final int j, final AbsoluteDate date) {
            return dij[j][i].value(date);
        }
    }

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

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

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

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

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

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

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

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

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

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

            //compute the coefficients
            computeU1V1Coefficients();
            computeU2V2Coefficients();
        }

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

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

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

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

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

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

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

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

                    // divide by 1 / j
                    u1ij[i][j] *= -ooj;
                    v1ij[i][j] *= ooj;

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

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

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

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

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

                // divide by 1 / j
                u1ij[0][j] *= -ooj;
                v1ij[0][j] *= ooj;
            }
        }

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

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

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

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

                for (int l = 1; l <= jMax; l++) {
                    // -U₁<sup>j+l</sup> * σ<sub>l</sub> +
                    // U₁<sup>l</sup> * σ<sub>j+l</sub> +
                    // V₁<sup>j+l</sup> * ρ<sub>l</sub> -
                    // V₁<sup>l</sup> * ρ<sub>j+l</sub>
                    u2ij[j] +=   -u1ij[0][j + l] * currentRhoSigmaj[1][l] +
                                  u1ij[0][l] * currentRhoSigmaj[1][j + l] +
                                  v1ij[0][j + l] * currentRhoSigmaj[0][l] -
                                  v1ij[0][l] * currentRhoSigmaj[0][j + l];

                    // U₁<sup>j+l</sup> * ρ<sub>l</sub> +
                    // U₁<sup>l</sup> * ρ<sub>j+l</sub> +
                    // V₁<sup>j+l</sup> * σ<sub>l</sub> +
                    // V₁<sup>l</sup> * σ<sub>j+l</sub>
                    u2ij[j] +=   u1ij[0][j + l] * currentRhoSigmaj[0][l] +
                                 u1ij[0][l] * currentRhoSigmaj[0][j + l] +
                                 v1ij[0][j + l] * currentRhoSigmaj[1][l] +
                                 v1ij[0][l] * currentRhoSigmaj[1][j + l];
                }

                // divide by 1 / j
                u2ij[j] *= -ooj;
                v2ij[j] *= ooj;
            }
        }

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

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

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

        /** Get the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0).
         *
         * @param j j index
         * @return the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0)
         */
        public double getV2(final int j) {
            return v2ij[j];
        }
    }
}