<commit_before>package org.orekit.propagation.semianalytical.dsst.dsstforcemodel;

import org.apache.commons.math.analysis.UnivariateFunction;
import org.apache.commons.math.analysis.integration.SimpsonIntegrator;
import org.apache.commons.math.analysis.integration.UnivariateRealIntegrator;
import org.apache.commons.math.geometry.euclidean.threed.Vector3D;
import org.apache.commons.math.util.FastMath;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitExceptionWrapper;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.SpacecraftState;

/** 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 (i.e. atmospheric drag and solar radiation pressure).
 * </p><p>
 * Gaussian contributions can be expressed as: da<sub>i</sub>/dt = &delta;a<sub>i</sub>/&delta;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>
 * Only two methods must be implemented by derived classes:
 * {@link #getAcceleration(SpacecraftState)} and {@link #getLLimits(SpacecraftState)}.
 * </p>
 * @author <NAME>
 */
public abstract class AbstractDSSTGaussianContribution implements DSSTForceModel {

    /** Propagation orbit type. */
    protected static final OrbitType ORBIT_TYPE = OrbitType.EQUINOCTIAL;

    /** Position angle type. */
    protected static final PositionAngle ANGLE_TYPE = PositionAngle.MEAN;

    /** DSST model needs equinoctial orbit as internal representation.
     *  Classical equinoctial elements have discontinuities when inclination is close to zero.
     *  In this representation, I = +1. <br>
     *  To avoid this discontinuity, another representation exists and equinoctial elements can
     *  be expressed in a different way, called "retrograde" orbit. This implies I = -1.
     *  As Orekit doesn't implement the retrograde orbit, I = +1 here.
     */
    protected double I = 1;

    // Equinoctial elements (according to DSST notation)
    /** a */
    protected double a;
    /** ex */
    protected double k;
    /** ey */
    protected double h;
    /** hx */
    protected double q;
    /** hy */
    protected double p;

    // Kepler mean motion
    /** n = sqrt(&mu; / a<sup>3</sup>) */
    protected double n;

    // Equinoctial reference frame vectors (according to DSST notation)
    /** f */
    protected Vector3D f;
    /** g */
    protected Vector3D g;
    /** w */
    protected Vector3D w;

    // Useful equinoctial coefficients
    /** A = sqrt(&mu; * a) */
    protected double A;
    /** B = sqrt(1 - h<sup>2</sup> - k<sup>2</sup>) */
    protected double B;
    /** C = 1 + p<sup>2</sup> + q<sup>2</sup> */
    protected double C;

    // Common factors
    /** 2. / (n<sup>2</sup> * 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. / &mu; */
    protected double ooMu;

    /** Build a new instance. */
    protected AbstractDSSTGaussianContribution() {
    }

    /** Get the current retrograde factor I.
     *  @return the retrograde factor I
     */
    public double getRetrogradeFactor() {
        return I;
    }

    /** {@inheritDoc} */
    public double[] getMeanElementRate(final SpacecraftState state) throws OrekitException {
        final double[] meanElementRate = new double[6];
        // Constant multiplier for integral
        final double h = state.getOrbit().getEquinoctialEy();
        final double k = state.getOrbit().getEquinoctialEx();
        final double coef = 1. / (2. * FastMath.PI * FastMath.sqrt(1. - h * h - k * k));
        final double[] ll = getLLimits(state);
        // Define integrable functions
        final IntegrableFunction iFct = new IntegrableFunction(state);
        // Compute mean element rates
        for (int i = 0; i < 6; i++) {
            // Select element
            iFct.setElement(i);
            // Define numerical quadrature operator parameters
//            final int    nbPts = getNbPoints(i);
//            final double relAc = getRelativeAccuracy(i);
//            final double absAc = getAbsoluteAccuracy(i);
            final int    maxEv = getMaxEval(i);
            // Define numerical quadrature operator
//            final UnivariateRealIntegrator numQuad = new LegendreGaussIntegrator(nbPts, relAc, absAc);
            final UnivariateRealIntegrator numQuad = new SimpsonIntegrator();
            try {
                // Numerical quadrature
                meanElementRate[i] = coef * numQuad.integrate(maxEv, iFct, ll[0], ll[1]);
            } catch (OrekitExceptionWrapper oew) {
                throw oew.getException();
            }
        }
        return meanElementRate;
    }

    /** Compute useful equinoctial parameters: A, B, C, f, g, w.
     *  @param  state current state information: date, kinematics, attitude
     */
    protected void computeParameters(final SpacecraftState state) {
        // Initialisation of A, B, C coefficients, f, g, w basis
    
        // Get current state vector (equinoctial elements):
        double[] stateVector = new double[6];
        ORBIT_TYPE.mapOrbitToArray(state.getOrbit(), ANGLE_TYPE, stateVector);
    
        // Equinoctial elements
        a = stateVector[0];
        k = stateVector[1];
        h = stateVector[2];
        q = stateVector[3];
        p = stateVector[4];
    
        // Factors
        final double k2 = k * k;<commit_msg>Updated orekit after an API change in Apache Commons Math.<commit_after>package org.orekit.propagation.semianalytical.dsst.dsstforcemodel;

import org.apache.commons.math.analysis.UnivariateFunction;
import org.apache.commons.math.analysis.integration.SimpsonIntegrator;
import org.apache.commons.math.analysis.integration.UnivariateIntegrator;
import org.apache.commons.math.geometry.euclidean.threed.Vector3D;
import org.apache.commons.math.util.FastMath;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitExceptionWrapper;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.SpacecraftState;

/** 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 (i.e. atmospheric drag and solar radiation pressure).
 * </p><p>
 * Gaussian contributions can be expressed as: da<sub>i</sub>/dt = &delta;a<sub>i</sub>/&delta;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>
 * Only two methods must be implemented by derived classes:
 * {@link #getAcceleration(SpacecraftState)} and {@link #getLLimits(SpacecraftState)}.
 * </p>
 * @author <NAME>
 */
public abstract class AbstractDSSTGaussianContribution implements DSSTForceModel {

    /** Propagation orbit type. */
    protected static final OrbitType ORBIT_TYPE = OrbitType.EQUINOCTIAL;

    /** Position angle type. */
    protected static final PositionAngle ANGLE_TYPE = PositionAngle.MEAN;

    /** DSST model needs equinoctial orbit as internal representation.
     *  Classical equinoctial elements have discontinuities when inclination is close to zero.
     *  In this representation, I = +1. <br>
     *  To avoid this discontinuity, another representation exists and equinoctial elements can
     *  be expressed in a different way, called "retrograde" orbit. This implies I = -1.
     *  As Orekit doesn't implement the retrograde orbit, I = +1 here.
     */
    protected double I = 1;

    // Equinoctial elements (according to DSST notation)
    /** a */
    protected double a;
    /** ex */
    protected double k;
    /** ey */
    protected double h;
    /** hx */
    protected double q;
    /** hy */
    protected double p;

    // Kepler mean motion
    /** n = sqrt(&mu; / a<sup>3</sup>) */
    protected double n;

    // Equinoctial reference frame vectors (according to DSST notation)
    /** f */
    protected Vector3D f;
    /** g */
    protected Vector3D g;
    /** w */
    protected Vector3D w;

    // Useful equinoctial coefficients
    /** A = sqrt(&mu; * a) */
    protected double A;
    /** B = sqrt(1 - h<sup>2</sup> - k<sup>2</sup>) */
    protected double B;
    /** C = 1 + p<sup>2</sup> + q<sup>2</sup> */
    protected double C;

    // Common factors
    /** 2. / (n<sup>2</sup> * 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. / &mu; */
    protected double ooMu;

    /** Build a new instance. */
    protected AbstractDSSTGaussianContribution() {
    }

    /** Get the current retrograde factor I.
     *  @return the retrograde factor I
     */
    public double getRetrogradeFactor() {
        return I;
    }

    /** {@inheritDoc} */
    public double[] getMeanElementRate(final SpacecraftState state) throws OrekitException {
        final double[] meanElementRate = new double[6];
        // Constant multiplier for integral
        final double h = state.getOrbit().getEquinoctialEy();
        final double k = state.getOrbit().getEquinoctialEx();
        final double coef = 1. / (2. * FastMath.PI * FastMath.sqrt(1. - h * h - k * k));
        final double[] ll = getLLimits(state);
        // Define integrable functions
        final IntegrableFunction iFct = new IntegrableFunction(state);
        // Compute mean element rates
        for (int i = 0; i < 6; i++) {
            // Select element
            iFct.setElement(i);
            // Define numerical quadrature operator parameters
//            final int    nbPts = getNbPoints(i);
//            final double relAc = getRelativeAccuracy(i);
//            final double absAc = getAbsoluteAccuracy(i);
            final int    maxEv = getMaxEval(i);
            // Define numerical quadrature operator
//            final UnivariateRealIntegrator numQuad = new LegendreGaussIntegrator(nbPts, relAc, absAc);
            final UnivariateIntegrator numQuad = new SimpsonIntegrator();
            try {
                // Numerical quadrature
                meanElementRate[i] = coef * numQuad.integrate(maxEv, iFct, ll[0], ll[1]);
            } catch (OrekitExceptionWrapper oew) {
                throw oew.getException();
            }
        }
        return meanElementRate;
    }

    /** Compute useful equinoctial parameters: A, B, C, f, g, w.
     *  @param  state current state information: date, kinematics, attitude
     */
    protected void computeParameters(final SpacecraftState state) {
        // Initialisation of A, B, C coefficients, f, g, w basis
    
        // Get current state vector (equinoctial elements):
        double[] stateVector = new double[6];
        ORBIT_TYPE.mapOrbitToArray(state.getOrbit(), ANGLE_TYPE, stateVector);
    
        // Equinoctial elements
        a = stateVector[0];
        k = stateVector[1];
        h = stateVector[2];
        q = stateVector[3];
        p = stateVector[4];
    
        // Factors
        final double k2 = k * k;
{'fraction_non_alphanumeric': 0.11286438003488082, 'fraction_numerical': 0.006643966447969438, 'mean_word_length': 2.6086305064429127, 'pattern_counts': {'":': 0, '<': 75, '<?xml version=': 0, '>': 73, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3116055', 'n_tokens_mistral': 3623, 'n_tokens_neox': 3457, 'n_words': 1080}