Examples with DormandPrince853Integrator org.hipparchus.ode.nonstiff.DormandPrince853Integrator used on opensource projects

Example 81 with DormandPrince853Integrator

use of org.hipparchus.ode.nonstiff.DormandPrince853Integrator in project Orekit by CS-SI.

the class MasterMode method main.

/**
 * Program entry point.
 * @param args program arguments (unused here)
 */
public static void main(String[] args) {
    try {
        // configure Orekit
        File home = new File(System.getProperty("user.home"));
        File orekitData = new File(home, "orekit-data");
        if (!orekitData.exists()) {
            System.err.format(Locale.US, "Failed to find %s folder%n", orekitData.getAbsolutePath());
            System.err.format(Locale.US, "You need to download %s from the %s page and unzip it in %s for this tutorial to work%n", "orekit-data.zip", "https://www.orekit.org/forge/projects/orekit/files", home.getAbsolutePath());
            System.exit(1);
        }
        DataProvidersManager manager = DataProvidersManager.getInstance();
        manager.addProvider(new DirectoryCrawler(orekitData));
        // gravitation coefficient
        double mu = 3.986004415e+14;
        // inertial frame
        Frame inertialFrame = FramesFactory.getEME2000();
        // Initial date
        AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, TimeScalesFactory.getUTC());
        // Initial orbit
        // semi major axis in meters
        double a = 24396159;
        // eccentricity
        double e = 0.72831215;
        // inclination
        double i = FastMath.toRadians(7);
        // perigee argument
        double omega = FastMath.toRadians(180);
        // right ascention of ascending node
        double raan = FastMath.toRadians(261);
        // mean anomaly
        double lM = 0;
        Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN, inertialFrame, initialDate, mu);
        // Initial state definition
        SpacecraftState initialState = new SpacecraftState(initialOrbit);
        // Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000
        final double minStep = 0.001;
        final double maxstep = 1000.0;
        final double positionTolerance = 10.0;
        final OrbitType propagationType = OrbitType.KEPLERIAN;
        final double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);
        AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
        // Propagator
        NumericalPropagator propagator = new NumericalPropagator(integrator);
        propagator.setOrbitType(propagationType);
        // Force Model (reduced to perturbing gravity field)
        final NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(10, 10);
        ForceModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
        // Add force model to the propagator
        propagator.addForceModel(holmesFeatherstone);
        // Set up initial state in the propagator
        propagator.setInitialState(initialState);
        // Set up operating mode for the propagator as master mode
        // with fixed step and specialized step handler
        propagator.setMasterMode(60., new TutorialStepHandler());
        // Extrapolate from the initial to the final date
        SpacecraftState finalState = propagator.propagate(initialDate.shiftedBy(630.));
        KeplerianOrbit o = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalState.getOrbit());
        System.out.format(Locale.US, "Final state:%n%s %12.3f %10.8f %10.6f %10.6f %10.6f %10.6f%n", finalState.getDate(), o.getA(), o.getE(), FastMath.toDegrees(o.getI()), FastMath.toDegrees(o.getPerigeeArgument()), FastMath.toDegrees(o.getRightAscensionOfAscendingNode()), FastMath.toDegrees(o.getTrueAnomaly()));
    } catch (OrekitException oe) {
        System.err.println(oe.getMessage());
    }
}

Example 82 with DormandPrince853Integrator

use of org.hipparchus.ode.nonstiff.DormandPrince853Integrator in project Orekit by CS-SI.

the class AngularCoordinatesTest method testShiftWithAcceleration.

@Test
public void testShiftWithAcceleration() throws OrekitException {
    double rate = 2 * FastMath.PI / (12 * 60);
    double acc = 0.001;
    double dt = 1.0;
    int n = 2000;
    final AngularCoordinates quadratic = new AngularCoordinates(Rotation.IDENTITY, new Vector3D(rate, Vector3D.PLUS_K), new Vector3D(acc, Vector3D.PLUS_J));
    final AngularCoordinates linear = new AngularCoordinates(quadratic.getRotation(), quadratic.getRotationRate(), Vector3D.ZERO);
    final OrdinaryDifferentialEquation ode = new OrdinaryDifferentialEquation() {

        public int getDimension() {
            return 4;
        }

        public double[] computeDerivatives(final double t, final double[] q) {
            final double omegaX = quadratic.getRotationRate().getX() + t * quadratic.getRotationAcceleration().getX();
            final double omegaY = quadratic.getRotationRate().getY() + t * quadratic.getRotationAcceleration().getY();
            final double omegaZ = quadratic.getRotationRate().getZ() + t * quadratic.getRotationAcceleration().getZ();
            return new double[] { 0.5 * MathArrays.linearCombination(-q[1], omegaX, -q[2], omegaY, -q[3], omegaZ), 0.5 * MathArrays.linearCombination(q[0], omegaX, -q[3], omegaY, q[2], omegaZ), 0.5 * MathArrays.linearCombination(q[3], omegaX, q[0], omegaY, -q[1], omegaZ), 0.5 * MathArrays.linearCombination(-q[2], omegaX, q[1], omegaY, q[0], omegaZ) };
        }
    };
    ODEIntegrator integrator = new DormandPrince853Integrator(1.0e-6, 1.0, 1.0e-12, 1.0e-12);
    integrator.addStepHandler(new StepNormalizer(dt / n, new ODEFixedStepHandler() {

        public void handleStep(ODEStateAndDerivative s, boolean isLast) {
            final double t = s.getTime();
            final double[] y = s.getPrimaryState();
            Rotation reference = new Rotation(y[0], y[1], y[2], y[3], true);
            // the error in shiftedBy taking acceleration into account is cubic
            double expectedCubicError = 1.4544e-6 * t * t * t;
            Assert.assertEquals(expectedCubicError, Rotation.distance(reference, quadratic.shiftedBy(t).getRotation()), 0.0001 * expectedCubicError);
            // the error in shiftedBy not taking acceleration into account is quadratic
            double expectedQuadraticError = 5.0e-4 * t * t;
            Assert.assertEquals(expectedQuadraticError, Rotation.distance(reference, linear.shiftedBy(t).getRotation()), 0.00001 * expectedQuadraticError);
        }
    }));
    double[] y = new double[] { quadratic.getRotation().getQ0(), quadratic.getRotation().getQ1(), quadratic.getRotation().getQ2(), quadratic.getRotation().getQ3() };
    integrator.integrate(ode, new ODEState(0, y), dt);
}

Example 83 with DormandPrince853Integrator

use of org.hipparchus.ode.nonstiff.DormandPrince853Integrator in project symja_android_library by axkr.

the class NDSolve method evaluate.

@Override
public IExpr evaluate(final IAST ast, EvalEngine engine) {
    if (!ToggleFeature.DSOLVE) {
        return F.NIL;
    }
    if (ast.arg3().isList()) {
        final IAST tRangeList = (IAST) ast.arg3();
        if (!(tRangeList.isAST2() || tRangeList.isAST3())) {
            return F.NIL;
        }
        try {
            final IAST listOfVariables = Validate.checkIsVariableOrVariableList(ast, 2, ast.topHead(), engine);
            if (!listOfVariables.isPresent()) {
                return F.NIL;
            }
            final int numberOfVariables = listOfVariables.argSize();
            final ISymbol timeVar = (ISymbol) tRangeList.arg1();
            IExpr tMinExpr = F.C0;
            IExpr tMaxExpr = tRangeList.arg2();
            if (tRangeList.isAST3()) {
                tMinExpr = tRangeList.arg2();
                tMaxExpr = tRangeList.arg3();
            }
            final double tMin = tMinExpr.evalDouble();
            final double tMax = tMaxExpr.evalDouble();
            final double tStep = 0.1;
            IASTAppendable listOfEquations = Validate.checkEquations(ast, 1).copyAppendable();
            IExpr[][] boundaryCondition = new IExpr[2][numberOfVariables];
            int i = 1;
            while (i < listOfEquations.size()) {
                IExpr equation = listOfEquations.get(i);
                if (equation.isFree(timeVar)) {
                    if (determineSingleBoundary(equation, listOfVariables, timeVar, boundaryCondition, engine)) {
                        listOfEquations.remove(i);
                        continue;
                    }
                }
                i++;
            }
            IExpr[] dotEquations = new IExpr[numberOfVariables];
            i = 1;
            while (i < listOfEquations.size()) {
                IExpr equation = listOfEquations.get(i);
                if (!equation.isFree(timeVar)) {
                    if (determineSingleDotEquation(equation, listOfVariables, timeVar, dotEquations, engine)) {
                        listOfEquations.remove(i);
                        continue;
                    }
                }
                i++;
            }
            if (listOfVariables.isList()) {
                AbstractIntegrator abstractIntegrator = new DormandPrince853Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
                // AbstractIntegrator abstractIntegrator = new ClassicalRungeKuttaIntegrator(1.0);
                double[] primaryState = new double[numberOfVariables];
                for (int j = 0; j < numberOfVariables; j++) {
                    primaryState[j] = ((INum) engine.evalN(boundaryCondition[1][j])).doubleValue();
                }
                OrdinaryDifferentialEquation ode = new FirstODE(engine, dotEquations, listOfVariables, timeVar);
                if (listOfVariables.size() > 1) {
                    IASTAppendable[] resultLists = new IASTAppendable[numberOfVariables];
                    for (int j = 0; j < primaryState.length; j++) {
                        resultLists[j] = F.ListAlloc();
                    }
                    for (double time = tMin; time < tMax; time += tStep) {
                        final ODEStateAndDerivative finalstate = abstractIntegrator.integrate(ode, new ODEState(time, primaryState), time + tStep);
                        primaryState = finalstate.getPrimaryState();
                        for (int j = 0; j < primaryState.length; j++) {
                            resultLists[j].append(F.list(F.num(time), F.num(primaryState[j])));
                        }
                    }
                    IASTAppendable primaryList = F.ListAlloc();
                    IASTAppendable secondaryList = F.ListAlloc();
                    for (int j = 1; j < listOfVariables.size(); j++) {
                        secondaryList.append(F.Rule(listOfVariables.get(j), engine.evaluate(F.Interpolation(resultLists[j - 1]))));
                    }
                    primaryList.append(secondaryList);
                    return primaryList;
                }
            }
        } catch (LimitException le) {
            throw le;
        } catch (RuntimeException rex) {
            LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
        }
    }
    return F.NIL;
}