Example usage for org.apache.commons.math.util MathUtils EPSILON

Introduction

In this page you can find the example usage for org.apache.commons.math.util MathUtils EPSILON.

Prototype

double EPSILON

To view the source code for org.apache.commons.math.util MathUtils EPSILON.

Click Source Link

Document

Smallest positive number such that 1 - EPSILON is not numerically equal to 1.

Usage

From source file:org.orekit.propagation.numerical.NumericalPropagatorWithJacobians.java

/** Propagate towards a target date and compute partial derivatives.
 * <p>Propagation is the same as the
 * {@link NumericalPropagator#propagate(AbsoluteDate) basic one}.</p>
 * <p>Jacobian for orbit parameters is given as a 7x7 array.</p>
 * <p>Partial derivatives will be computed as a 7xn array
 * when n parameters have been {@link #selectParameters(String[]) selected}
 * (n may be 0).</p>/*from  ww  w.  j av a2 s.c  o  m*/
 * <p>Those parameters are related to some {@link ForceModelWithJacobians force models}
 * which must have been added elsewhere.</p>
 * @param finalDate target date towards which orbit state should be propagated
 * @param dFdY equinoctial orbit parameters + mass jacobian (7x7 array)
 * @param dFdP partial derivatives with respect to selected parameters (7xn)
 * @return propagated state
 * @exception PropagationException if state cannot be propagated
 */
public SpacecraftState propagate(final AbsoluteDate finalDate, final double[][] dFdY, final double[][] dFdP)
        throws PropagationException {
    try {

        if (initialState == null) {
            throw new PropagationException("initial state not specified for orbit propagation");
        }
        if (initialState.getMass() <= 0.0) {
            throw new IllegalArgumentException("Mass is null or negative");
        }
        if (initialState.getDate().equals(finalDate)) {
            // don't extrapolate
            return initialState;
        }
        if (integrator == null) {
            throw new PropagationException("ODE integrator not set for orbit propagation");
        }

        // space dynamics view
        startDate = initialState.getDate();
        if (modeHandler != null) {
            modeHandler.initialize(startDate, initialState.getFrame(), mu, attitudeLaw);
        }

        final EquinoctialOrbit initialOrbit = new EquinoctialOrbit(initialState.getOrbit());

        currentState = new SpacecraftState(initialOrbit, initialState.getAttitude(), initialState.getMass());

        adder = new TimeDerivativesEquationsWithJacobians(initialOrbit);

        integrator.clearEventHandlers();

        // set up events related to force models
        for (final ForceModelWithJacobians forceModel : forceModelsWJ) {
            final EventDetector[] modelDetectors = forceModel.getEventsDetectors();
            if (modelDetectors != null) {
                for (final EventDetector detector : modelDetectors) {
                    setUpEventDetector(detector);
                }
            }
        }

        // set up events added by user
        for (final EventDetector detector : detectors) {
            setUpEventDetector(detector);
        }

        // mathematical view
        final double t0 = 0;
        final double t1 = finalDate.durationFrom(startDate);

        // Map state to array
        stateVector[0] = initialOrbit.getA();
        stateVector[1] = initialOrbit.getEquinoctialEx();
        stateVector[2] = initialOrbit.getEquinoctialEy();
        stateVector[3] = initialOrbit.getHx();
        stateVector[4] = initialOrbit.getHy();
        stateVector[5] = initialOrbit.getLv();
        stateVector[6] = initialState.getMass();

        // set up parameters for jacobian computation
        int noParam = 0;
        final int nbParam = selectedParameters.length;
        final double[] paramWJ = new double[nbParam];
        final double[] hP = new double[nbParam];

        for (final String parameter : selectedParameters) {
            boolean found = false;
            for (final ForceModelWithJacobians fmwj : forceModelsWJ) {
                for (String parFMWJ : fmwj.getParametersNames()) {
                    if (parFMWJ.matches(parameter)) {
                        found = true;
                        paramWJ[noParam] = fmwj.getParameter(parFMWJ);
                        hP[noParam] = paramWJ[noParam] * Math.sqrt(MathUtils.EPSILON);
                        paramPairs.add(new ParameterPair(parameter, fmwj));
                        noParam++;
                    }
                }
            }
            if (!found) {
                throw new PropagationException("unknown parameter {0}", parameter);
            }
        }

        // if selectParameters was not invoked and then no parameter selected
        if (DY0DP == null) {
            DY0DP = new double[7][nbParam];
            for (final double[] row : DY0DP) {
                Arrays.fill(row, 0.0);
            }
        }

        // get hY from integrator tolerance array
        final double[] hY = getHy(integrator);

        // resize integrator tolerance array
        expandToleranceArray(integrator);

        final FirstOrderIntegratorWithJacobians integratorWJ = new FirstOrderIntegratorWithJacobians(integrator,
                new DifferentialEquations(), paramWJ, hY, hP);

        try {
            // mathematical integration
            final double stopTime = integratorWJ.integrate(t0, stateVector, DY0DP, t1, stateVector, dFdY,
                    DY0DP);
            // fill in jacobian
            for (int i = 0; i < DY0DP.length; i++) {
                System.arraycopy(DY0DP[i], 0, dFdP[i], 0, DY0DP[i].length);
            }

            // back to space dynamics view
            final AbsoluteDate date = startDate.shiftedBy(stopTime);

            final EquinoctialOrbit orbit = new EquinoctialOrbit(stateVector[0], stateVector[1], stateVector[2],
                    stateVector[3], stateVector[4], stateVector[5], EquinoctialOrbit.TRUE_LATITUDE_ARGUMENT,
                    initialOrbit.getFrame(), date, mu);

            resetInitialState(new SpacecraftState(orbit, attitudeLaw.getAttitude(orbit), stateVector[6]));
        } finally {
            if (integrator != null) {
                resetToleranceArray(integrator);
            }
        }

        return initialState;

    } catch (OrekitException oe) {
        throw new PropagationException(oe.getMessage(), oe);
    } catch (DerivativeException de) {

        // recover a possible embedded PropagationException
        for (Throwable t = de; t != null; t = t.getCause()) {
            if (t instanceof PropagationException) {
                throw (PropagationException) t;
            }
        }

        throw new PropagationException(de.getMessage(), de);

    } catch (IntegratorException ie) {

        // recover a possible embedded PropagationException
        for (Throwable t = ie; t != null; t = t.getCause()) {
            if (t instanceof PropagationException) {
                throw (PropagationException) t;
            }
        }

        throw new PropagationException(ie.getMessage(), ie);

    }
}

---