Example usage for org.apache.commons.math3.ode.nonstiff ClassicalRungeKuttaIntegrator ClassicalRungeKuttaIntegrator

List of usage examples for org.apache.commons.math3.ode.nonstiff ClassicalRungeKuttaIntegrator ClassicalRungeKuttaIntegrator

Introduction

In this page you can find the example usage for org.apache.commons.math3.ode.nonstiff ClassicalRungeKuttaIntegrator ClassicalRungeKuttaIntegrator.

Prototype

public ClassicalRungeKuttaIntegrator(final double step) 

Source Link

Document

Simple constructor.

Usage

From source file:fr.cs.examples.propagation.EphemerisMode.java

/** Program entry point.
 * @param args program arguments (unused here)
 *///from w  w w  .j av a 2 s.  co m
public static void main(String[] args) {
    try {

        // configure Orekit
        Autoconfiguration.configureOrekit();

        // Initial orbit parameters
        double a = 24396159; // semi major axis in meters
        double e = 0.72831215; // eccentricity
        double i = FastMath.toRadians(7); // inclination
        double omega = FastMath.toRadians(180); // perigee argument
        double raan = FastMath.toRadians(261); // right ascension of ascending node
        double lM = 0; // mean anomaly

        // Inertial frame
        Frame inertialFrame = FramesFactory.getEME2000();

        // Initial date in UTC time scale
        TimeScale utc = TimeScalesFactory.getUTC();
        AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, utc);

        // gravitation coefficient
        double mu = 3.986004415e+14;

        // Orbit construction as Keplerian
        Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN, inertialFrame,
                initialDate, mu);

        // Initialize state
        SpacecraftState initialState = new SpacecraftState(initialOrbit);

        // Numerical propagation with no perturbation (only keplerian movement)
        // Using a very simple integrator with a fixed step: classical Runge-Kutta
        double stepSize = 10; // the step is ten seconds
        AbstractIntegrator integrator = new ClassicalRungeKuttaIntegrator(stepSize);
        NumericalPropagator propagator = new NumericalPropagator(integrator);

        // Set the propagator to ephemeris mode
        propagator.setEphemerisMode();

        // Initialize propagation
        propagator.setInitialState(initialState);

        // Propagation with storage of the results in an integrated ephemeris
        SpacecraftState finalState = propagator.propagate(initialDate.shiftedBy(6000));

        System.out.println(" Numerical propagation :");
        System.out.println("  Final date : " + finalState.getDate());
        System.out.println("  " + finalState.getOrbit());

        // Getting the integrated ephemeris
        BoundedPropagator ephemeris = propagator.getGeneratedEphemeris();

        System.out
                .println(" Ephemeris defined from " + ephemeris.getMinDate() + " to " + ephemeris.getMaxDate());

        System.out.println(" Ephemeris propagation :");
        AbsoluteDate intermediateDate = initialDate.shiftedBy(3000);
        SpacecraftState intermediateState = ephemeris.propagate(intermediateDate);
        System.out.println("  date :  " + intermediateState.getDate());
        System.out.println("  " + intermediateState.getOrbit());

        intermediateDate = finalState.getDate();
        intermediateState = ephemeris.propagate(intermediateDate);
        System.out.println("  date :  " + intermediateState.getDate());
        System.out.println("  " + intermediateState.getOrbit());

        intermediateDate = initialDate.shiftedBy(-1000);
        System.out.println();
        System.out.println("Attempting to propagate to date " + intermediateDate + " which is OUT OF RANGE");
        System.out.println("This propagation attempt should fail, " + "so an error message shoud appear below, "
                + "this is expected and shows that errors are handled correctly");
        intermediateState = ephemeris.propagate(intermediateDate);

        // these two print should never happen as en exception should have been triggered
        System.out.println("  date :  " + intermediateState.getDate());
        System.out.println("  " + intermediateState.getOrbit());

    } catch (OrekitException oe) {
        System.out.println(oe.getMessage());
    }
}

From source file:de.bund.bfr.math.IntegratorFactory.java

public FirstOrderIntegrator createIntegrator() {
    switch (type) {
    case EULER://from   www.  j av  a  2s. c o  m
        return new EulerIntegrator(step);
    case GILL:
        return new GillIntegrator(step);
    case MIDPOINT:
        return new MidpointIntegrator(step);
    case RUNGE_KUTTA:
        return new ClassicalRungeKuttaIntegrator(step);
    case THREE_EIGHTHES:
        return new ThreeEighthesIntegrator(step);
    default:
        throw new RuntimeException("Unknown type of IntegratorFactory: " + type);
    }
}

From source file:beast.structuredCoalescent.distribution.ExactStructuredCoalescent.java

public double calculateLogP() {
    // Calculate the tree intervals (time between events, which nodes participate at a event etc.)       
    treeIntervalsInput.get().calculateIntervals();
    treeIntervalsInput.get().swap();//w  w w  .  j a v  a  2  s. c om

    // Set up for lineage state probabilities
    activeLineages = new ArrayList<Integer>();
    lineStateProbs = new ArrayList<Double>();

    // Compute likelihood at each integration time and tree event starting at final sampling time and moving backwards
    logP = 0;

    // Initialize the line state probabilities
    // total number of intervals
    final int intervalCount = treeIntervalsInput.get().getIntervalCount();
    // counts in which interval we are in
    int t = 0;
    nr_lineages = 0;
    // Captures the probabilities of lineages being in a state
    double[] p;

    // Initialize the migration rates matrix
    int c = 0;

    for (int k = 0; k < states; k++) {
        for (int l = 0; l < states; l++) {
            if (k != l) {
                migration_rates[c] = migrationRatesInput.get().getArrayValue(c);
                migration_map[k][l] = c;
                c++;
            } else {
                coalescent_rates[k] = coalescentRatesInput.get().getArrayValue(k) / 2;
            }

        }
    }

    boolean first = true;
    // integrate until there are no more tree intervals
    do {
        double nextIntervalTime = treeIntervalsInput.get().getInterval(t);
        // Length of the current interval
        final double duration = nextIntervalTime;// - currTime;
        // if the current interval has a length greater than 0, integrate
        if (duration > 0) {
            p = new double[jointStateProbabilities.size()]; // Captures the probabilities of lineages being in a state

            // convert the array list to double[]
            for (int i = 0; i < jointStateProbabilities.size(); i++)
                p[i] = jointStateProbabilities.get(i);

            double[] p_for_ode = new double[p.length];

            double ts = timeStep;
            if (duration < timeStep)
                ts = duration / 2;

            // initialize integrator
            FirstOrderIntegrator integrator = new ClassicalRungeKuttaIntegrator(ts);
            // set the odes
            FirstOrderDifferentialEquations ode = new ode_integrator(migration_rates, coalescent_rates,
                    nr_lineages, states, connectivity, sums);
            // integrate                   
            integrator.integrate(ode, 0, p, duration, p_for_ode);

            // if the dimension is equal to the max integer, this means that a calculation
            // of a probability of a configuration resulted in a value below 0 and the
            // run will be stopped
            if (ode.getDimension() == Integer.MAX_VALUE) {
                System.out.println("lalalallal");
                return Double.NEGATIVE_INFINITY;
            }

            // set the probabilities of the system being in a configuration again
            for (int i = 0; i < p_for_ode.length; i++)
                jointStateProbabilities.set(i, p_for_ode[i]);

        }
        /*
         *  compute contribution of event to likelihood
         */
        if (treeIntervalsInput.get().getIntervalType(t) == IntervalType.COALESCENT) {
            nr_lineages--;
            logP += coalesce(t);
        }

        /*
         * add new lineage
         */
        if (treeIntervalsInput.get().getIntervalType(t) == IntervalType.SAMPLE) {
            nr_lineages++;
            addLineages(t, first);
            first = false;
        }

        t++;
    } while (t < intervalCount);

    //Compute likelihood of remaining tree intervals (coal events occuring before origin)
    if (Double.isInfinite(logP))
        logP = Double.NEGATIVE_INFINITY;
    if (max_posterior < logP && logP < 0) {
        max_posterior = logP;
        max_mig = new double[states * (states - 1)];
        max_coal = new double[states];
        for (int i = 0; i < 1; i++)
            max_mig[i] = migrationRatesInput.get().getArrayValue(i);
        for (int i = 0; i < 1; i++)
            max_coal[i] = coalescentRatesInput.get().getArrayValue(i);
    }

    return logP;
}

From source file:eu.mihosoft.fx.tutorials.gravity.SolarSystem.java

/**
 * Initializes the ui frame listener.//w w w. j  ava2  s. c  o  m
 *
 * @param ode ode to solve
 * @param y state vector (location and velocity)
 * @param m particle masses
 * @param ignoreFlags ignore flags (used for collisions)
 * @param dt time step size
 */
private void initFrameListener(FirstOrderDifferentialEquations ode, double[] y, final double[] m,
        final boolean[] ignoreFlags, double dt) {

    // local time scale (sec per h * h per day * days per year)
    final double localTimeScale = 3600 * 24 * 365;

    // integrator
    FirstOrderIntegrator integrator = new ClassicalRungeKuttaIntegrator(dt * localTimeScale);

    double[] yPrev = new double[y.length];
    double[] interpolatedY = new double[y.length];

    // create frame listener 
    frameListener = new AnimationTimer() {

        @Override
        public void handle(long now) {

            // thanks to http://gafferongames.com/game-physics/fix-your-timestep/
            // measure elapsed time between last and current pulse (frame)
            double frameDuration = (now - lastTimeStamp) / 1e9;
            lastTimeStamp = now;

            // we don't allow frame durations above 2*dt
            if (frameDuration > 2 * dt) {
                frameDuration = 2 * dt;
            }

            // add elapsed time to remaining simulation interval
            remainingSimulationTime += frameDuration;

            // copy current state to prev state
            System.arraycopy(y, 0, yPrev, 0, yPrev.length);

            // simulate remaining interval
            while (remainingSimulationTime >= dt) {

                double scaledT = time * localTimeScale;
                double scaledDT = dt * localTimeScale * timeScale;

                // integrate one step
                try {
                    integrator.integrate(ode, scaledT, y, scaledT + scaledDT, y);

                } catch (Exception ex) {
                    ex.printStackTrace(System.err);
                }
                // remove integrated interval from remaining simulation time
                remainingSimulationTime -= dt;

                // update t
                time = time + dt;
            }

            // interpolate state
            double alpha = remainingSimulationTime / dt;

            // set interpolated state
            for (int i = 0; i < y.length; i++) {
                interpolatedY[i] = y[i] * alpha + yPrev[i] * (1.0 - alpha);
            }

            // update properties for visualization
            updateView(nodes, interpolatedY, m, ignoreFlags);
        }
    };

    // finally, start the framle listener
    frameListener.start();
}

From source file:beast.structuredCoalescent.distribution.Masco.java

public double calculateLogP() {
    // newly calculate tree intervals
    treeIntervalsInput.get().calculateIntervals();
    // correctly calculate the daughter nodes at coalescent intervals in the case of
    // bifurcation or in case two nodes are at the same height
    treeIntervalsInput.get().swap();//www  .j ava2s  .  c  om

    // Set up ArrayLists for the indices of active lineages and the lineage state probabilities
    activeLineages = new ArrayList<Integer>();
    lineStateProbs = new ArrayList<Double>();

    // Compute likelihood at each integration time and tree event starting at final sampling time and moving backwards
    logP = 0;

    // set the current time
    double currTime = 0.0;
    // total number of intervals
    final int intervalCount = treeIntervalsInput.get().getIntervalCount();
    // interval time counter
    int t = 0;
    // initialize the number of lineages
    nr_lineages = 0;
    // Captures the probabilities of lineages being in a state
    double[] p;

    // Initialize the migration rates matrix
    double[][] migration_rates = new double[states][states];
    int c = 0;

    for (int k = 0; k < states; k++) {
        for (int l = 0; l < states; l++) {
            if (k != l) {
                migration_rates[k][l] = migrationRatesInput.get().getArrayValue(c);
                c++;
            } else { // diagonal
                migration_rates[k][l] = 0.0;
            }

        }
    }

    // Initialize the coalescent rates
    double[] coalescent_rates = new double[states];
    for (int k = 0; k < states; k++) {
        coalescent_rates[k] = coalescentRatesInput.get().getArrayValue(k) / 2;//(epiModelInput.get().getF(t,k,k) / (Y.get(k)*Y.get(k)));
    }

    // integrate until there are no more tree intervals
    do {
        double nextIntervalTime = treeIntervalsInput.get().getInterval(t);

        // Length of the current interval
        final double duration = nextIntervalTime;// - currTime;
        // if the current interval has a length greater than 0, integrate
        if (duration > 0) {
            if (dependentHistory)
                p = new double[lineStateProbs.size()]; // Captures the probabilities of lineages being in a state
            else
                p = new double[lineStateProbs.size() + 1]; // Captures the probabilities of lineages being in a state, last one keeps track of the probability

            // convert the array list to double[]
            for (int i = 0; i < lineStateProbs.size(); i++)
                p[i] = lineStateProbs.get(i);

            // not needed
            if (!dependentHistory)
                p[lineStateProbs.size()] = 1;

            double[] p_for_ode = new double[p.length];
            double ts = 0.0;

            // If proportial time step is true, set the integration time for the given interval 
            // inverse proportional to the number of lineages
            if (propTimeStep)
                ts = timeStep / lineStateProbs.size();
            else
                ts = timeStep;

            // Never choose a longer time step than the integration window
            if (duration < (ts / 2))
                ts = duration / 2;

            FirstOrderIntegrator integrator = new ClassicalRungeKuttaIntegrator(ts);
            // set the odes
            FirstOrderDifferentialEquations ode = new ode_masco(migration_rates, coalescent_rates, nr_lineages,
                    states);
            // integrate                   
            integrator.integrate(ode, 0, p, duration, p_for_ode);

            // If the Dimension is larger than the maximum integer, at least one state prob is below 0 and the step is rejected
            if (ode.getDimension() == Integer.MAX_VALUE) {
                System.out.println(lineStateProbs.size());
                System.out.println("lalalallal");
                return Double.NEGATIVE_INFINITY;
            }

            for (int i = 0; i < lineStateProbs.size(); i++)
                lineStateProbs.set(i, p_for_ode[i]);
        }

        // update the time
        currTime = nextIntervalTime;
        // event is coalescent event
        if (treeIntervalsInput.get().getIntervalType(t) == IntervalType.COALESCENT) {
            logP += coalesce(t);
            nr_lineages--;
        }

        // event is sampling event
        if (treeIntervalsInput.get().getIntervalType(t) == IntervalType.SAMPLE) {
            logP += normalizeLineages();
            addLineages(t);
            nr_lineages++;
        }

        // update the interval number
        t++;
    } while (t < intervalCount);

    //Compute likelihood of remaining tree intervals (coal events occuring before origin)
    if (Double.isInfinite(logP))
        logP = Double.NEGATIVE_INFINITY;
    if (max_posterior < logP && logP < 0) {
        max_posterior = logP;
        max_mig = new double[states * (states - 1)];
        max_coal = new double[states];
        for (int i = 0; i < 1; i++)
            max_mig[i] = migrationRatesInput.get().getArrayValue(i);
        for (int i = 0; i < 1; i++)
            max_coal[i] = coalescentRatesInput.get().getArrayValue(i);
    }
    //        System.exit(0);

    return logP;

}

From source file:beast.structuredCoalescent.distribution.IndependentStructuredCoalescent.java

public double calculateLogP() {
    // newly calculate tree intervals
    treeIntervalsInput.get().calculateIntervals();
    // correctly calculate the daughter nodes at coalescent intervals in the case of
    // bifurcation or in case two nodes are at the same height
    treeIntervalsInput.get().swap();/* w  w w  . j a v  a2 s. c o  m*/

    // Set up ArrayLists for the indices of active lineages and the lineage state probabilities
    activeLineages = new ArrayList<Integer>();
    lineStateProbs = new ArrayList<Double>();

    // Compute likelihood at each integration time and tree event starting at final sampling time and moving backwards
    logP = 0;

    // set the current time
    double currTime = 0.0;
    // total number of intervals
    final int intervalCount = treeIntervalsInput.get().getIntervalCount();
    // interval time counter
    int t = 0;
    // initialize the number of lineages
    nr_lineages = 0;
    // Captures the probabilities of lineages being in a state
    double[] p;

    // Initialize the migration rates matrix
    double[][] migration_rates = new double[states][states];
    int c = 0;

    for (int k = 0; k < states; k++) {
        for (int l = 0; l < states; l++) {
            if (k != l) {
                migration_rates[k][l] = migrationRatesInput.get().getArrayValue(c);
                c++;
            } else { // diagonal
                migration_rates[k][l] = 0.0;
            }

        }
    }

    // Initialize the coalescent rates
    double[] coalescent_rates = new double[states];
    for (int k = 0; k < states; k++) {
        coalescent_rates[k] = coalescentRatesInput.get().getArrayValue(k) / 2;//(epiModelInput.get().getF(t,k,k) / (Y.get(k)*Y.get(k)));
    }

    // integrate until there are no more tree intervals
    do {
        double nextIntervalTime = treeIntervalsInput.get().getInterval(t);

        // Length of the current interval
        final double duration = nextIntervalTime;// - currTime;
        // if the current interval has a length greater than 0, integrate
        if (duration > 0) {
            if (dependentHistory)
                p = new double[lineStateProbs.size()]; // Captures the probabilities of lineages being in a state
            else
                p = new double[lineStateProbs.size() + 1]; // Captures the probabilities of lineages being in a state, last one keeps track of the probability

            // convert the array list to double[]
            for (int i = 0; i < lineStateProbs.size(); i++)
                p[i] = lineStateProbs.get(i);

            // not needed
            if (!dependentHistory)
                p[lineStateProbs.size()] = 1;

            double[] p_for_ode = new double[p.length];
            double ts = 0.0;

            // If proportial time step is true, set the integration time for the given interval 
            // inverse proportional to the number of lineages
            if (propTimeStep)
                ts = timeStep / lineStateProbs.size();
            else
                ts = timeStep;

            // Never choose a longer time step than the integration window
            if (duration < (ts / 2))
                ts = duration / 2;

            FirstOrderIntegrator integrator = new ClassicalRungeKuttaIntegrator(ts);
            // set the odes
            FirstOrderDifferentialEquations ode = new independent_ode_integrator(migration_rates,
                    coalescent_rates, nr_lineages, states);
            // integrate                   
            integrator.integrate(ode, 0, p, duration, p_for_ode);

            // If the Dimension is larger than the maximum integer, at least one state prob is below 0 and the step is rejected
            if (ode.getDimension() == Integer.MAX_VALUE) {
                System.out.println(lineStateProbs.size());
                System.out.println("lalalallal");
                return Double.NEGATIVE_INFINITY;
            }

            for (int i = 0; i < lineStateProbs.size(); i++)
                lineStateProbs.set(i, p_for_ode[i]);
        }

        // update the time
        currTime = nextIntervalTime;
        // event is coalescent event
        if (treeIntervalsInput.get().getIntervalType(t) == IntervalType.COALESCENT) {
            logP += coalesce(t);
            nr_lineages--;
        }

        // event is sampling event
        if (treeIntervalsInput.get().getIntervalType(t) == IntervalType.SAMPLE) {
            addLineages(t);
            nr_lineages++;
        }

        // update the interval number
        t++;
    } while (t < intervalCount);

    //Compute likelihood of remaining tree intervals (coal events occuring before origin)
    if (Double.isInfinite(logP))
        logP = Double.NEGATIVE_INFINITY;
    if (max_posterior < logP && logP < 0) {
        max_posterior = logP;
        max_mig = new double[states * (states - 1)];
        max_coal = new double[states];
        for (int i = 0; i < 1; i++)
            max_mig[i] = migrationRatesInput.get().getArrayValue(i);
        for (int i = 0; i < 1; i++)
            max_coal[i] = coalescentRatesInput.get().getArrayValue(i);
    }

    return logP;

}

From source file:fr.cs.examples.propagation.DSSTPropagation.java

/** Set up the DSST Propagator
 *
 *  @param orbit initial orbit//from   w w  w  .j a va 2 s .co m
 *  @param mass S/C mass (kg)
 *  @param isOsculating if orbital elements are osculating
 *  @param fixedStepSize step size for fixed step integrator (s)
 *  @throws OrekitException
 */
private DSSTPropagator createDSSTProp(final Orbit orbit, final double mass, final boolean isOsculating,
        final double fixedStepSize) throws OrekitException {
    AbstractIntegrator integrator;
    if (fixedStepSize > 0.) {
        integrator = new ClassicalRungeKuttaIntegrator(fixedStepSize);
    } else {
        final double minStep = orbit.getKeplerianPeriod();
        final double maxStep = minStep * 100.;
        final double[][] tol = DSSTPropagator.tolerances(1.0, orbit);
        integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
        ((AdaptiveStepsizeIntegrator) integrator).setInitialStepSize(10. * minStep);
    }

    DSSTPropagator dsstProp = new DSSTPropagator(integrator);
    dsstProp.setInitialState(new SpacecraftState(orbit, mass), isOsculating);

    return dsstProp;
}

From source file:org.orekit.forces.gravity.CunninghamAttractionModelTest.java

@Test
public void testZonalWithDrozinerReference() throws OrekitException, ParseException {
    // initialization
    AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 07, 01), new TimeComponents(13, 59, 27.816),
            TimeScalesFactory.getUTC());
    double i = FastMath.toRadians(98.7);
    double omega = FastMath.toRadians(93.0);
    double OMEGA = FastMath.toRadians(15.0 * 22.5);
    Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i, omega, OMEGA, 0, PositionAngle.MEAN,
            FramesFactory.getEME2000(), date, mu);

    propagator = new NumericalPropagator(new ClassicalRungeKuttaIntegrator(1000));
    propagator.addForceModel(new CunninghamAttractionModel(itrf2008,
            GravityFieldFactory.getUnnormalizedProvider(ae, mu, TideSystem.UNKNOWN,
                    new double[][] { { 0.0 }, { 0.0 }, { c20 }, { c30 }, { c40 }, { c50 }, { c60 }, },
                    new double[][] { { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, })));

    propagator.setInitialState(new SpacecraftState(orbit));
    SpacecraftState cunnOrb = propagator.propagate(date.shiftedBy(Constants.JULIAN_DAY));

    propagator.removeForceModels();//from  w w  w . j a  va2s . co  m

    propagator.addForceModel(new DrozinerAttractionModel(itrf2008,
            GravityFieldFactory.getUnnormalizedProvider(ae, mu, TideSystem.UNKNOWN,
                    new double[][] { { 0.0 }, { 0.0 }, { c20 }, { c30 }, { c40 }, { c50 }, { c60 }, },
                    new double[][] { { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, })));

    propagator.setInitialState(new SpacecraftState(orbit));
    SpacecraftState drozOrb = propagator.propagate(date.shiftedBy(Constants.JULIAN_DAY));

    Vector3D dif = cunnOrb.getPVCoordinates().getPosition().subtract(drozOrb.getPVCoordinates().getPosition());
    Assert.assertEquals(0, dif.getNorm(), 3.1e-7);
    Assert.assertTrue(propagator.getCalls() < 400);
}

From source file:org.orekit.forces.gravity.DrozinerAttractionModelTest.java

@Test
public void testTesserealWithHolmesFeaterstoneReference() throws OrekitException, IOException, ParseException {

    Utils.setDataRoot("regular-data:potential/grgs-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
    UnnormalizedSphericalHarmonicsProvider unnormalized = GravityFieldFactory.getUnnormalizedProvider(10, 10);
    NormalizedSphericalHarmonicsProvider normalized = GravityFieldFactory.getNormalizedProvider(10, 10);

    //  initialization
    AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 07, 01), new TimeComponents(13, 59, 27.816),
            TimeScalesFactory.getUTC());
    double i = FastMath.toRadians(98.7);
    double omega = FastMath.toRadians(93.0);
    double OMEGA = FastMath.toRadians(15.0 * 22.5);
    Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i, omega, OMEGA, 0, PositionAngle.MEAN,
            FramesFactory.getEME2000(), date, mu);
    propagator = new NumericalPropagator(new ClassicalRungeKuttaIntegrator(100));
    propagator.addForceModel(new HolmesFeatherstoneAttractionModel(itrf2008, normalized));
    propagator.setInitialState(new SpacecraftState(orbit));
    SpacecraftState hfOrb = propagator.propagate(date.shiftedBy(Constants.JULIAN_DAY));

    propagator.removeForceModels();//from  w  w  w  .j  a  v  a 2  s  . co m
    propagator.addForceModel(new DrozinerAttractionModel(itrf2008, unnormalized));

    propagator.setInitialState(new SpacecraftState(orbit));
    SpacecraftState drozOrb = propagator.propagate(date.shiftedBy(Constants.JULIAN_DAY));

    Vector3D dif = hfOrb.getPVCoordinates().getPosition().subtract(drozOrb.getPVCoordinates().getPosition());
    Assert.assertEquals(0, dif.getNorm(), 3.1e-7);
    Assert.assertTrue(propagator.getCalls() < 3500);

}

From source file:org.orekit.forces.gravity.HolmesFeatherstoneAttractionModelTest.java

@Test
public void testZonalWithCunninghamReference() throws OrekitException {
    // initialization
    AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 07, 01), new TimeComponents(13, 59, 27.816),
            TimeScalesFactory.getUTC());
    double i = FastMath.toRadians(98.7);
    double omega = FastMath.toRadians(93.0);
    double OMEGA = FastMath.toRadians(15.0 * 22.5);
    Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i, omega, OMEGA, 0, PositionAngle.MEAN,
            FramesFactory.getEME2000(), date, mu);

    propagator = new NumericalPropagator(new ClassicalRungeKuttaIntegrator(1000));
    propagator.addForceModel(new HolmesFeatherstoneAttractionModel(itrf,
            GravityFieldFactory.getNormalizedProvider(ae, mu, TideSystem.UNKNOWN,
                    new double[][] { { 0.0 }, { 0.0 }, { normalizedC20 }, { normalizedC30 }, { normalizedC40 },
                            { normalizedC50 }, { normalizedC60 }, },
                    new double[][] { { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, })));

    propagator.setInitialState(new SpacecraftState(orbit));
    SpacecraftState hfOrb = propagator.propagate(date.shiftedBy(Constants.JULIAN_DAY));

    propagator.removeForceModels();//from  w w w  . j a  va 2  s . co  m

    propagator.addForceModel(new CunninghamAttractionModel(itrf,
            GravityFieldFactory.getUnnormalizedProvider(ae, mu, TideSystem.UNKNOWN,
                    new double[][] { { 0.0 }, { 0.0 }, { unnormalizedC20 }, { unnormalizedC30 },
                            { unnormalizedC40 }, { unnormalizedC50 }, { unnormalizedC60 }, },
                    new double[][] { { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, })));

    propagator.setInitialState(new SpacecraftState(orbit));
    SpacecraftState cOrb = propagator.propagate(date.shiftedBy(Constants.JULIAN_DAY));

    Vector3D dif = hfOrb.getPVCoordinates().getPosition().subtract(cOrb.getPVCoordinates().getPosition());
    Assert.assertEquals(0, dif.getNorm(), 2e-9);
    Assert.assertTrue(propagator.getCalls() < 400);
}