GLONASSNumericalPropagator.java

/* Copyright 2002-2024 CS GROUP
 * Licensed to CS GROUP (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.numerical;

import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.ode.ODEIntegrator;
import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;
import org.hipparchus.util.SinCos;
import org.orekit.attitudes.Attitude;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.data.DataContext;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.frames.Frame;
import org.orekit.orbits.CartesianOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.PropagationType;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.analytical.gnss.data.GLONASSAlmanac;
import org.orekit.propagation.analytical.gnss.data.GLONASSNavigationMessage;
import org.orekit.propagation.analytical.gnss.data.GLONASSOrbitalElements;
import org.orekit.propagation.analytical.gnss.data.GNSSConstants;
import org.orekit.propagation.integration.AbstractIntegratedPropagator;
import org.orekit.propagation.integration.StateMapper;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.GLONASSDate;
import org.orekit.utils.AbsolutePVCoordinates;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;

import java.util.Arrays;

/**
 * This class propagates GLONASS orbits using numerical integration.
 * <p>
 * As recommended by the GLONASS Interface Control Document (ICD),
 * a {@link ClassicalRungeKuttaIntegrator  4th order Runge-Kutta technique}
 * shall be used to integrate the equations.
 * </p>
 * <p>
 * Classical used of this orbit propagator is to compute GLONASS satellite
 * coordinates from the navigation message.
 * </p>
 * <p>
 * If the projections of luni-solar accelerations to axes of
 * Greenwich geocentric coordinates {@link GLONASSOrbitalElements#getXDotDot() X''(tb)},
 * {@link GLONASSOrbitalElements#getYDotDot() Y''(tb)} and {@link GLONASSOrbitalElements#getZDotDot() Z''(tb)}
 * are available in the navigation message; a transformation is performed to convert these
 * accelerations into the correct coordinate system. In the case where they are not
 * available into the navigation message, these accelerations are computed.
 * </p>
 * <p>
 * <b>Caution:</b> The Glonass numerical propagator can only be used with {@link GLONASSNavigationMessage}.
 * Using this propagator with a {@link GLONASSAlmanac} is prone to error.
 * </p>
 *
 * @see <a href="http://russianspacesystems.ru/wp-content/uploads/2016/08/ICD-GLONASS-CDMA-General.-Edition-1.0-2016.pdf">
 *       GLONASS Interface Control Document</a>
 *
 * @author Bryan Cazabonne
 *
 */
public class GLONASSNumericalPropagator extends AbstractIntegratedPropagator {

    /** Second degree coefficient of normal potential. */
    private static final double GLONASS_J20 = 1.08262575e-3;

    /** Equatorial radius of Earth (m). */
    private static final double GLONASS_EARTH_EQUATORIAL_RADIUS = 6378136;

    /** Value of the Earth's rotation rate in rad/s (See Ref). */
    private static final double GLONASS_AV = 7.2921151467e-5;

    // Data used to solve Kepler's equation
    /** First coefficient to compute Kepler equation solver starter. */
    private static final double A;

    /** Second coefficient to compute Kepler equation solver starter. */
    private static final double B;

    static {
        final double k1 = 3 * FastMath.PI + 2;
        final double k2 = FastMath.PI - 1;
        final double k3 = 6 * FastMath.PI - 1;
        A  = 3 * k2 * k2 / k1;
        B  = k3 * k3 / (6 * k1);
    }

    /** The GLONASS orbital elements used. */
    private final GLONASSOrbitalElements glonassOrbit;

    /** Initial date in GLONASS form. */
    private final GLONASSDate initDate;

    /** The spacecraft mass (kg). */
    private final double mass;

    /** The ECI frame used for GLONASS propagation. */
    private final Frame eci;

    /** Direction cosines and distance of perturbing body: Moon.
     * <p>
     * <ul>
     * <li>double[0] = ξ<sub>m</sub></li>
     * <li>double[1] = η<sub>m</sub></li>
     * <li>double[2] = ψ<sub>m</sub></li>
     * <li>double[3] = r<sub>m</sub></li>
     * </ul>
     * </p>
     */
    private double[] moonElements;

    /** Direction cosines and distance of perturbing body: Sun.
     * <p>
     * <ul>
     * <li>double[0] = ξ<sub>s</sub></li>
     * <li>double[1] = η<sub>s</sub></li>
     * <li>double[2] = ψ<sub>s</sub></li>
     * <li>double[3] = r<sub>s</sub></li>
     * </ul>
     * </p>
     */
    private double[] sunElements;

    /** Flag for availability of projections of acceleration transmitted within the navigation message. */
    private final boolean isAccAvailable;

    /** Data context used for propagation. */
    private final DataContext dataContext;

    /**
     * Private constructor.
     * @param integrator Runge-Kutta integrator
     * @param glonassOrbit Glonass orbital elements
     * @param eci Earth Centered Inertial frame
     * @param provider Attitude provider
     * @param mass Satellite mass (kg)
     * @param context Data context
     * @param isAccAvailable true if the acceleration  is transmitted within the navigation message
     */
    public GLONASSNumericalPropagator(final ClassicalRungeKuttaIntegrator integrator,
                                      final GLONASSOrbitalElements glonassOrbit,
                                      final Frame eci, final AttitudeProvider provider,
                                      final double mass, final DataContext context,
                                      final boolean isAccAvailable) {
        super(integrator, PropagationType.OSCULATING);
        this.dataContext = context;
        this.isAccAvailable = isAccAvailable;
        // Stores the GLONASS orbital elements
        this.glonassOrbit = glonassOrbit;
        // Sets the Earth Centered Inertial frame
        this.eci = eci;
        // Sets the mass
        this.mass = mass;
        this.initDate = new GLONASSDate(
                glonassOrbit.getDate(),
                dataContext.getTimeScales().getGLONASS());

        // Initialize state mapper
        initMapper();
        setInitialState();
        setAttitudeProvider(provider);
        setOrbitType(OrbitType.CARTESIAN);
        // It is not meaningful for propagation in Cartesian parameters
        setPositionAngleType(PositionAngleType.TRUE);
        setMu(GNSSConstants.GLONASS_MU);

        // As recommended by GLONASS ICD (2016), the direction cosines and distance
        // of perturbing body are calculated one time (at tb).
        if (!isAccAvailable) {
            computeMoonElements(initDate);
            computeSunElements(initDate);
        }
    }

    /**
     * Gets the underlying GLONASS orbital elements.
     *
     * @return the underlying GLONASS orbital elements
     */
    public GLONASSOrbitalElements getGLONASSOrbitalElements() {
        return glonassOrbit;
    }

    /** {@inheritDoc} */
    @Override
    public SpacecraftState propagate(final AbsoluteDate date) {
        // Spacecraft state in inertial frame
        final SpacecraftState stateInInertial = super.propagate(date);

        // Build the spacecraft state in inertial frame
        final PVCoordinates pvInPZ90 = getPVInPZ90(stateInInertial);
        final AbsolutePVCoordinates absPV = new AbsolutePVCoordinates(
                dataContext.getFrames().getPZ9011(IERSConventions.IERS_2010, true),
                stateInInertial.getDate(), pvInPZ90);
        final TimeStampedPVCoordinates pvInInertial = absPV.getPVCoordinates(eci);
        return new SpacecraftState(new CartesianOrbit(pvInInertial, eci, pvInInertial.getDate(), GNSSConstants.GLONASS_MU),
                                                      stateInInertial.getAttitude(),
                                                      stateInInertial.getMass(),
                                                      stateInInertial.getAdditionalStatesValues(),
                                                      stateInInertial.getAdditionalStatesDerivatives());
    }

    /**
     * Set the initial state.
     * <p>
     * The initial conditions on position and velocity are in the ECEF coordinate system PZ-90.
     * Previous to orbit integration, they must be transformed to an absolute inertial coordinate system.
     * </p>
     */
    private void setInitialState() {

        // Transform initial PV coordinates to an absolute inertial coordinate system.
        final PVCoordinates pvInInertial = getPVInInertial(initDate);

        // Create a new orbit
        final Orbit orbit = new CartesianOrbit(pvInInertial,
                                               eci, initDate.getDate(),
                                               GNSSConstants.GLONASS_MU);

        // Reset the initial state to apply the transformation
        resetInitialState(new SpacecraftState(orbit, mass));
    }

    /**
     * This method computes the direction cosines and the distance used to
     * compute the gravitational perturbations of the Moon.
     *
     * @param date the computation date in GLONASS scale
     */
    private void computeMoonElements(final GLONASSDate date) {

        moonElements = new double[4];

        // Constants
        // Semi-major axis of the Moon's orbit (m)
        final double am = 3.84385243e8;
        // The Moon's orbit eccentricity
        final double em = 0.054900489;
        // Mean inclination of the Moon's orbit to the ecliptic (rad)
        final double im = 0.0898041080;

        // Computed parameters
        // Time from epoch 2000 to the instant tb of GLONASS ephemeris broadcast
        final double dtoJD = (glonassOrbit.getTime() - 10800.) / Constants.JULIAN_DAY;
        final double t  = (date.getJD0() + dtoJD - 2451545.0) / 36525;
        final double t2 = t * t;
        // Mean inclination of Earth equator to ecliptic (rad)
        final double eps = 0.4090926006 - 0.0002270711 * t;
        // Mean longitude of the Moon's orbit perigee (rad)
        final double gammaM = 1.4547885346 + 71.0176852437 * t - 0.0001801481 * t2;
        // Mean longitude of the ascending node of the Moon (rad)
        final double omegaM = 2.1824391966 - 33.7570459536 * t + 0.0000362262 * t2;
        // Mean anomaly of the Moon (rad)
        final double qm = 2.3555557435 + 8328.6914257190 * t + 0.0001545547 * t2;

        // Commons parameters
        final SinCos scOm  = FastMath.sinCos(omegaM);
        final SinCos scIm  = FastMath.sinCos(im);
        final SinCos scEs  = FastMath.sinCos(eps);
        final SinCos scGm  = FastMath.sinCos(gammaM);
        final double cosOm = scOm.cos();
        final double sinOm = scOm.sin();
        final double cosIm = scIm.cos();
        final double sinIm = scIm.sin();
        final double cosEs = scEs.cos();
        final double sinEs = scEs.sin();
        final double cosGm = scGm.cos();
        final double sinGm = scGm.sin();

        // Intermediate parameters
        final double psiStar = cosOm * sinIm;
        final double etaStar = sinOm * sinIm;
        final double epsStar = 1. - cosOm * cosOm * (1. - cosIm);
        final double eps11 = sinOm * cosOm * (1. - cosIm);
        final double eps12 = 1. - sinOm * sinOm * (1. - cosIm);
        final double eta11 = epsStar * cosEs - psiStar * sinEs;
        final double eta12 = eps11 * cosEs + etaStar * sinEs;
        final double psi11 = epsStar * sinEs + psiStar * cosEs;
        final double psi12 = eps11 * sinEs - etaStar * cosEs;

        // Eccentric Anomaly
        final double ek = getEccentricAnomaly(qm, em);

        // True Anomaly
        final double vk    = getTrueAnomaly(ek, em);
        final SinCos scVk  = FastMath.sinCos(vk);
        final double sinVk = scVk.sin();
        final double cosVk = scVk.cos();

        // Direction cosine
        final double epsM = eps11 * (sinVk * cosGm + cosVk * sinGm) + eps12 * (cosVk * cosGm - sinVk * sinGm);
        final double etaM = eta11 * (sinVk * cosGm + cosVk * sinGm) + eta12 * (cosVk * cosGm - sinVk * sinGm);
        final double psiM = psi11 * (sinVk * cosGm + cosVk * sinGm) + psi12 * (cosVk * cosGm - sinVk * sinGm);

        // Distance
        final double rm = am * (1. - em * FastMath.cos(ek));

        moonElements[0] = epsM;
        moonElements[1] = etaM;
        moonElements[2] = psiM;
        moonElements[3] = rm;

    }

    /**
     * This method computes the direction cosines and the distance used to
     * compute the gravitational perturbations of the Sun.
     *
     * @param date the computation date in GLONASS scale
     */
    private void computeSunElements(final GLONASSDate date) {

        sunElements = new double[4];

        // Constants
        //  Major semi-axis of the Earth’s orbit around the Sun (m)
        final double as = 1.49598e11;
        // The eccentricity of the Earth’s orbit around the Sun
        final double es = 0.016719;

        // Computed parameters
        // Time from epoch 2000 to the instant tb of GLONASS ephemeris broadcast
        final double dtoJD = (glonassOrbit.getTime() - 10800.) / Constants.JULIAN_DAY;
        final double t  = (date.getJD0() + dtoJD - 2451545.0) / 36525;
        final double t2 = t * t;
        // Mean inclination of Earth equator to ecliptic (rad)
        final double eps = 0.4090926006 - 0.0002270711 * t;
        // Mean tropic longitude of the Sun orbit perigee (rad)
        final double ws = -7.6281824375 + 0.0300101976 * t + 0.0000079741 * t2;
        // Mean anomaly of the Sun (rad)
        final double qs = 6.2400601269 + 628.3019551714 * t - 0.0000026820 * t2;

        // Eccentric Anomaly
        final double ek = getEccentricAnomaly(qs, es);

        // True Anomaly
        final double vk    =  getTrueAnomaly(ek, es);
        final SinCos scVk  = FastMath.sinCos(vk);
        final double sinVk = scVk.sin();
        final double cosVk = scVk.cos();

        // Commons parameters
        final SinCos scWs  = FastMath.sinCos(ws);
        final SinCos scEs  = FastMath.sinCos(eps);
        final double cosWs = scWs.cos();
        final double sinWs = scWs.sin();
        final double cosEs = scEs.cos();
        final double sinEs = scEs.sin();

        // Direction cosine
        final double epsS = cosVk * cosWs - sinVk * sinWs;
        final double etaS = cosEs * (sinVk * cosWs + cosVk * sinWs);
        final double psiS = sinEs * (sinVk * cosWs + cosVk * sinWs);

        // Distance
        final double rs = as * (1. - es * FastMath.cos(ek));

        sunElements[0] = epsS;
        sunElements[1] = etaS;
        sunElements[2] = psiS;
        sunElements[3] = rs;

    }

    /**
     * Computes the elliptic eccentric anomaly from the mean anomaly.
     * <p>
     * The algorithm used here for solving Kepler equation has been published
     * in: "Procedures for  solving Kepler's Equation", A. W. Odell and
     * R. H. Gooding, Celestial Mechanics 38 (1986) 307-334
     * </p>
     * <p>It has been copied from the OREKIT library (KeplerianOrbit class).</p>
     *
     * @param M mean anomaly (rad)
     * @param e eccentricity
     * @return E the eccentric anomaly
     */
    private double getEccentricAnomaly(final double M, final double e) {

        // reduce M to [-PI PI) interval
        final double reducedM = MathUtils.normalizeAngle(M, 0.0);

        // compute start value according to A. W. Odell and R. H. Gooding S12 starter
        double E;
        if (FastMath.abs(reducedM) < 1.0 / 6.0) {
            E = reducedM + e * (FastMath.cbrt(6 * reducedM) - reducedM);
        } else {
            if (reducedM < 0) {
                final double w = FastMath.PI + reducedM;
                E = reducedM + e * (A * w / (B - w) - FastMath.PI - reducedM);
            } else {
                final double w = FastMath.PI - reducedM;
                E = reducedM + e * (FastMath.PI - A * w / (B - w) - reducedM);
            }
        }

        final double e1 = 1 - e;
        final boolean noCancellationRisk = (e1 + E * E / 6) >= 0.1;

        // perform two iterations, each consisting of one Halley step and one Newton-Raphson step
        for (int j = 0; j < 2; ++j) {
            final double f;
            double fd;
            final SinCos scE  = FastMath.sinCos(E);
            final double fdd  = e * scE.sin();
            final double fddd = e * scE.cos();
            if (noCancellationRisk) {
                f  = (E - fdd) - reducedM;
                fd = 1 - fddd;
            } else {
                f  = eMeSinE(E, e) - reducedM;
                final double s = FastMath.sin(0.5 * E);
                fd = e1 + 2 * e * s * s;
            }
            final double dee = f * fd / (0.5 * f * fdd - fd * fd);

            // update eccentric anomaly, using expressions that limit underflow problems
            final double w = fd + 0.5 * dee * (fdd + dee * fddd / 3);
            fd += dee * (fdd + 0.5 * dee * fddd);
            E  -= (f - dee * (fd - w)) / fd;

        }

        // expand the result back to original range
        E += M - reducedM;

        return E;

    }

    /**
     * Accurate computation of E - e sin(E).
     *
     * @param E eccentric anomaly
     * @param e eccentricity
     * @return E - e sin(E)
     */
    private static double eMeSinE(final double E, final double e) {
        double x = (1 - e) * FastMath.sin(E);
        final double mE2 = -E * E;
        double term = E;
        double d    = 0;
        // the inequality test below IS intentional and should NOT be replaced by a check with a small tolerance
        for (double x0 = Double.NaN; !Double.valueOf(x).equals(x0);) {
            d += 2;
            term *= mE2 / (d * (d + 1));
            x0 = x;
            x = x - term;
        }
        return x;
    }

    /**
     * Get true anomaly from eccentric anomaly and eccentricity.
     *
     * @param ek the eccentric anomaly (rad)
     * @param ecc the eccentricity
     * @return the true anomaly (rad)
     */
    private double getTrueAnomaly(final double ek, final double ecc) {
        final SinCos scek = FastMath.sinCos(ek);
        final double svk  = FastMath.sqrt(1. - ecc * ecc) * scek.sin();
        final double cvk  = scek.cos() - ecc;
        return FastMath.atan2(svk, cvk);
    }

    /**
     * This method transforms the PV coordinates obtained after the numerical
     * integration in the ECEF PZ-90.
     *
     * @param state spacecraft state after integration
     * @return the PV coordinates in the ECEF PZ-90.
     */
    private PVCoordinates getPVInPZ90(final SpacecraftState state) {

        // Compute time difference between start date and end date
        final double dt = state.getDate().durationFrom(initDate.getDate());

        // Position and velocity vectors
        final PVCoordinates pv = state.getPVCoordinates();
        final Vector3D pos = pv.getPosition();
        final Vector3D vel = pv.getVelocity();

        // Components of position and velocity vectors
        final double x0 = pos.getX();
        final double y0 = pos.getY();
        final double z0 = pos.getZ();
        final double vx0 = vel.getX();
        final double vy0 = vel.getY();
        final double vz0 = vel.getZ();

        // Greenwich Mean Sidereal Time (GMST)
        final GLONASSDate gloDate = new GLONASSDate(
                state.getDate(),
                dataContext.getTimeScales().getGLONASS());
        final double gmst = gloDate.getGMST();

        final double ti = glonassOrbit.getTime() + dt;
        // We use the GMST instead of the GMT as it is recommended into GLONASS ICD (2016)
        final double s = gmst + GLONASS_AV * (ti - 10800.);

        // Commons Parameters
        final SinCos scS  = FastMath.sinCos(s);
        final double cosS = scS.cos();
        final double sinS = scS.sin();

        // Transformed coordinates
        final double x = x0 * cosS + y0 * sinS;
        final double y = -x0 * sinS + y0 * cosS;
        final double z = z0;
        final double vx = vx0 * cosS + vy0 * sinS + GLONASS_AV * y;
        final double vy = -vx0 * sinS + vy0 * cosS - GLONASS_AV * x;
        final double vz = vz0;

        // Transformed orbit
        return new PVCoordinates(new Vector3D(x, y, z),
                                 new Vector3D(vx, vy, vz));
    }

    /**
     * This method computes the PV coordinates of the spacecraft center of mass.
     * The returned PV are expressed in inertial coordinates system at the instant tb.
     *
     * @param date the computation date in GLONASS scale
     * @return the PV Coordinates in inertial coordinates system
     */
    private PVCoordinates getPVInInertial(final GLONASSDate date) {

        // Greenwich Mean Sidereal Time (GMST)
        final double gmst = date.getGMST();

        final double time = glonassOrbit.getTime();
        final double dt   = time - 10800.;
        // We use the GMST instead of the GMT as it is recommended into GLONASS ICD (2016)
        final double s = gmst + GLONASS_AV * dt;

        // Commons Parameters
        final SinCos scS  = FastMath.sinCos(s);
        final double cosS = scS.cos();
        final double sinS = scS.sin();

        // PV coordinates in inertial frame
        final double x0  = glonassOrbit.getX() * cosS - glonassOrbit.getY() * sinS;
        final double y0  = glonassOrbit.getX() * sinS + glonassOrbit.getY() * cosS;
        final double z0  = glonassOrbit.getZ();
        final double vx0 = glonassOrbit.getXDot() * cosS - glonassOrbit.getYDot() * sinS - GLONASS_AV * y0;
        final double vy0 = glonassOrbit.getXDot() * sinS + glonassOrbit.getYDot() * cosS + GLONASS_AV * x0;
        final double vz0 = glonassOrbit.getZDot();
        return new PVCoordinates(new Vector3D(x0, y0, z0),
                                 new Vector3D(vx0, vy0, vz0));
    }

    @Override
    protected StateMapper createMapper(final AbsoluteDate referenceDate, final double mu,
                                       final OrbitType orbitType, final PositionAngleType positionAngleType,
                                       final AttitudeProvider attitudeProvider, final Frame frame) {
        return new Mapper(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
    }

    /** Internal mapper. */
    private static class Mapper extends StateMapper {

        /**
         * Simple constructor.
         *
         * @param referenceDate reference date
         * @param mu central attraction coefficient (m³/s²)
         * @param orbitType orbit type to use for mapping
         * @param positionAngleType angle type to use for propagation
         * @param attitudeProvider attitude provider
         * @param frame inertial frame
         */
        Mapper(final AbsoluteDate referenceDate, final double mu,
               final OrbitType orbitType, final PositionAngleType positionAngleType,
               final AttitudeProvider attitudeProvider, final Frame frame) {
            super(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
        }

        @Override
        public SpacecraftState mapArrayToState(final AbsoluteDate date, final double[] y,
                                               final double[] yDot, final PropagationType type) {
            // The parameter meanOnly is ignored for the GLONASS Propagator
            final double mass = y[6];
            if (mass <= 0.0) {
                throw new OrekitException(OrekitMessages.NOT_POSITIVE_SPACECRAFT_MASS, mass);
            }

            final Orbit orbit       = getOrbitType().mapArrayToOrbit(y, yDot, getPositionAngleType(), date, getMu(), getFrame());
            final Attitude attitude = getAttitudeProvider().getAttitude(orbit, date, getFrame());

            return new SpacecraftState(orbit, attitude, mass);
        }

        @Override
        public void mapStateToArray(final SpacecraftState state, final double[] y,
                                    final double[] yDot) {
            getOrbitType().mapOrbitToArray(state.getOrbit(), getPositionAngleType(), y, yDot);
            y[6] = state.getMass();
        }

    }

    @Override
    protected MainStateEquations getMainStateEquations(final ODEIntegrator integ) {
        return new Main();
    }

    /** Internal class for orbital parameters integration. */
    private class Main implements MainStateEquations {

        /** Derivatives array. */
        private final double[] yDot;

        /**
         * Simple constructor.
         */
        Main() {
            yDot = new double[7];
        }

        @Override
        public double[] computeDerivatives(final SpacecraftState state) {

            // Date in Glonass form
            final GLONASSDate gloDate = new GLONASSDate(
                    state.getDate(),
                    dataContext.getTimeScales().getGLONASS());

            // Position and Velocity vectors
            final Vector3D vel = state.getPVCoordinates().getVelocity();
            final Vector3D pos = state.getPosition();

            Arrays.fill(yDot, 0.0);

            // dPos/dt = Vel
            yDot[0] += vel.getX();
            yDot[1] += vel.getY();
            yDot[2] += vel.getZ();

            // Components of position and velocity vectors
            final double x0 = pos.getX();
            final double y0 = pos.getY();
            final double z0 = pos.getZ();

            // Normalized values
            final double r  = pos.getNorm();
            final double r2 = r * r;
            final double oor = 1. / r;
            final double oor2 = 1. / r2;
            final double x = x0 * oor;
            final double y = y0 * oor;
            final double z = z0 * oor;
            final double g = GNSSConstants.GLONASS_MU * oor2;
            final double ro = GLONASS_EARTH_EQUATORIAL_RADIUS * oor;

            yDot[3] += x * (-g + (-1.5 * GLONASS_J20 * g * ro * ro * (1. - 5. * z * z)));
            yDot[4] += y * (-g + (-1.5 * GLONASS_J20 * g * ro * ro * (1. - 5. * z * z)));
            yDot[5] += z * (-g + (-1.5 * GLONASS_J20 * g * ro * ro * (3. - 5. * z * z)));

            // Luni-Solar contribution
            final Vector3D acc;
            if (isAccAvailable) {
                acc = getLuniSolarPerturbations(gloDate);
            } else {
                final Vector3D accMoon = computeLuniSolarPerturbations(
                        state, moonElements[0], moonElements[1], moonElements[2],
                        moonElements[3],
                        dataContext.getCelestialBodies().getMoon().getGM());
                final Vector3D accSun = computeLuniSolarPerturbations(
                        state,
                        sunElements[0], sunElements[1], sunElements[2],
                        sunElements[3],
                        dataContext.getCelestialBodies().getSun().getGM());
                acc = accMoon.add(accSun);
            }

            yDot[3] += acc.getX();
            yDot[4] += acc.getY();
            yDot[5] += acc.getZ();

            return yDot.clone();
        }

        /**
         * This method computes the accelerations induced by gravitational
         * perturbations of the Sun and the Moon if they are not available into
         * the navigation message data.
         *
         * @param state current state
         * @param eps first direction cosine
         * @param eta second direction cosine
         * @param psi third direction cosine
         * @param r distance of perturbing body
         * @param g body gravitational field constant
         * @return a vector containing the accelerations
         */
        private Vector3D computeLuniSolarPerturbations(final SpacecraftState state, final double eps,
                                                       final double eta, final double psi,
                                                       final double r, final double g) {

            // Current pv coordinates
            final PVCoordinates pv = state.getPVCoordinates();

            final double oor = 1. / r;
            final double oor2 = oor * oor;

            // Normalized variable
            final double x = pv.getPosition().getX() * oor;
            final double y = pv.getPosition().getY() * oor;
            final double z = pv.getPosition().getZ() * oor;
            final double gm = g * oor2;

            final double epsmX  = eps - x;
            final double etamY  = eta - y;
            final double psimZ  = psi - z;
            final Vector3D vector = new Vector3D(epsmX, etamY, psimZ);
            final double d2 = vector.getNormSq();
            final double deltaM = FastMath.sqrt(d2) * d2;

            // Accelerations
            final double accX = gm * ((epsmX / deltaM) - eps);
            final double accY = gm * ((etamY / deltaM) - eta);
            final double accZ = gm * ((psimZ / deltaM) - psi);

            return new Vector3D(accX, accY, accZ);
        }

        /**
         * Get the accelerations induced by gravitational
         * perturbations of the Sun and the Moon in a geocentric
         * coordinate system.
         * <p>
         * The accelerations are obtained using projections of accelerations
         * transmitted within navigation message data.
         * </p>
         *
         * @param date the computation date in GLONASS scale
         * @return a vector containing the sum of both accelerations
         */
        private Vector3D getLuniSolarPerturbations(final GLONASSDate date) {

            // Greenwich Mean Sidereal Time (GMST)
            final double gmst = date.getGMST();

            final double time = glonassOrbit.getTime();
            final double dt   = time - 10800.;
            // We use the GMST instead of the GMT as it is recommended into GLONASS ICD (see Ref)
            final double s = gmst + GLONASS_AV * dt;

            // Commons Parameters
            final SinCos scS  = FastMath.sinCos(s);
            final double cosS = scS.cos();
            final double sinS = scS.sin();

            // Accelerations
            final double accX = glonassOrbit.getXDotDot() * cosS - glonassOrbit.getYDotDot() * sinS;
            final double accY = glonassOrbit.getXDotDot() * sinS + glonassOrbit.getYDotDot() * cosS;
            final double accZ = glonassOrbit.getZDotDot();

            return new Vector3D(accX, accY, accZ);
        }

    }

}