Examples with AbstractIntegrator org.hipparchus.ode.AbstractIntegrator used on opensource projects

Example 6 with AbstractIntegrator

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

the class DragForceTest method testIssue229.

@Test
public void testIssue229() throws OrekitException {
    AbsoluteDate initialDate = new AbsoluteDate(2004, 1, 1, 0, 0, 0., TimeScalesFactory.getUTC());
    Frame frame = FramesFactory.getEME2000();
    double rpe = 160.e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
    double rap = 2000.e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
    double inc = FastMath.toRadians(0.);
    double aop = FastMath.toRadians(0.);
    double raan = FastMath.toRadians(0.);
    double mean = FastMath.toRadians(180.);
    double mass = 100.;
    KeplerianOrbit orbit = new KeplerianOrbit(0.5 * (rpe + rap), (rap - rpe) / (rpe + rap), inc, aop, raan, mean, PositionAngle.MEAN, frame, initialDate, Constants.EIGEN5C_EARTH_MU);
    IsotropicDrag shape = new IsotropicDrag(10., 2.2);
    Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
    BodyShape earthShape = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf);
    Atmosphere atmosphere = new SimpleExponentialAtmosphere(earthShape, 2.6e-10, 200000, 26000);
    double[][] tolerance = NumericalPropagator.tolerances(0.1, orbit, OrbitType.CARTESIAN);
    AbstractIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 300, tolerance[0], tolerance[1]);
    NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.setOrbitType(OrbitType.CARTESIAN);
    propagator.setMu(orbit.getMu());
    propagator.addForceModel(new DragForce(atmosphere, shape));
    PartialDerivativesEquations partials = new PartialDerivativesEquations("partials", propagator);
    propagator.setInitialState(partials.setInitialJacobians(new SpacecraftState(orbit, mass)));
    SpacecraftState state = propagator.propagate(new AbsoluteDate(2004, 1, 1, 1, 30, 0., TimeScalesFactory.getUTC()));
    double delta = 0.1;
    Orbit shifted = new CartesianOrbit(new TimeStampedPVCoordinates(orbit.getDate(), orbit.getPVCoordinates().getPosition().add(new Vector3D(delta, 0, 0)), orbit.getPVCoordinates().getVelocity()), orbit.getFrame(), orbit.getMu());
    propagator.setInitialState(partials.setInitialJacobians(new SpacecraftState(shifted, mass)));
    SpacecraftState newState = propagator.propagate(new AbsoluteDate(2004, 1, 1, 1, 30, 0., TimeScalesFactory.getUTC()));
    double[] dPVdX = new double[] { (newState.getPVCoordinates().getPosition().getX() - state.getPVCoordinates().getPosition().getX()) / delta, (newState.getPVCoordinates().getPosition().getY() - state.getPVCoordinates().getPosition().getY()) / delta, (newState.getPVCoordinates().getPosition().getZ() - state.getPVCoordinates().getPosition().getZ()) / delta, (newState.getPVCoordinates().getVelocity().getX() - state.getPVCoordinates().getVelocity().getX()) / delta, (newState.getPVCoordinates().getVelocity().getY() - state.getPVCoordinates().getVelocity().getY()) / delta, (newState.getPVCoordinates().getVelocity().getZ() - state.getPVCoordinates().getVelocity().getZ()) / delta };
    double[][] dYdY0 = new double[6][6];
    partials.getMapper().getStateJacobian(state, dYdY0);
    for (int i = 0; i < 6; ++i) {
        Assert.assertEquals(dPVdX[i], dYdY0[i][0], 6.2e-6 * FastMath.abs(dPVdX[i]));
    }
}

Example 7 with AbstractIntegrator

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

the class PropagatorConversion 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));
        // gravity field
        NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(2, 0);
        double mu = provider.getMu();
        // inertial frame
        Frame inertialFrame = FramesFactory.getEME2000();
        // Initial date
        AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, TimeScalesFactory.getUTC());
        // Initial orbit (GTO)
        // semi major axis in meters
        final double a = 24396159;
        // eccentricity
        final double e = 0.72831215;
        // inclination
        final double i = FastMath.toRadians(7);
        // perigee argument
        final double omega = FastMath.toRadians(180);
        // right ascention of ascending node
        final double raan = FastMath.toRadians(261);
        // mean anomaly
        final double lM = 0;
        Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN, inertialFrame, initialDate, mu);
        final double period = initialOrbit.getKeplerianPeriod();
        // Initial state definition
        final 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.;
        final double dP = 1.e-2;
        final OrbitType orbType = OrbitType.CARTESIAN;
        final double[][] tol = NumericalPropagator.tolerances(dP, initialOrbit, orbType);
        final AbstractIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
        // Propagator
        NumericalPropagator numProp = new NumericalPropagator(integrator);
        numProp.setInitialState(initialState);
        numProp.setOrbitType(orbType);
        // Force Models:
        // 1 - Perturbing gravity field (only J2 is considered here)
        ForceModel gravity = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
        // Add force models to the propagator
        numProp.addForceModel(gravity);
        // Propagator factory
        PropagatorBuilder builder = new KeplerianPropagatorBuilder(initialOrbit, PositionAngle.TRUE, dP);
        // Propagator converter
        PropagatorConverter fitter = new FiniteDifferencePropagatorConverter(builder, 1.e-6, 5000);
        // Resulting propagator
        KeplerianPropagator kepProp = (KeplerianPropagator) fitter.convert(numProp, 2 * period, 251);
        // Step handlers
        StatesHandler numStepHandler = new StatesHandler();
        StatesHandler kepStepHandler = new StatesHandler();
        // Set up operating mode for the propagator as master mode
        // with fixed step and specialized step handler
        numProp.setMasterMode(60., numStepHandler);
        kepProp.setMasterMode(60., kepStepHandler);
        // Extrapolate from the initial to the final date
        numProp.propagate(initialDate.shiftedBy(10. * period));
        kepProp.propagate(initialDate.shiftedBy(10. * period));
        // retrieve the states
        List<SpacecraftState> numStates = numStepHandler.getStates();
        List<SpacecraftState> kepStates = kepStepHandler.getStates();
        // Print the results on the output file
        File output = new File(new File(System.getProperty("user.home")), "elements.dat");
        try (final PrintStream stream = new PrintStream(output, "UTF-8")) {
            stream.println("# date Anum Akep Enum Ekep Inum Ikep LMnum LMkep");
            for (SpacecraftState numState : numStates) {
                for (SpacecraftState kepState : kepStates) {
                    if (numState.getDate().compareTo(kepState.getDate()) == 0) {
                        stream.println(numState.getDate() + " " + numState.getA() + " " + kepState.getA() + " " + numState.getE() + " " + kepState.getE() + " " + FastMath.toDegrees(numState.getI()) + " " + FastMath.toDegrees(kepState.getI()) + " " + FastMath.toDegrees(MathUtils.normalizeAngle(numState.getLM(), FastMath.PI)) + " " + FastMath.toDegrees(MathUtils.normalizeAngle(kepState.getLM(), FastMath.PI)));
                        break;
                    }
                }
            }
        }
        System.out.println("Results saved as file " + output);
        File output1 = new File(new File(System.getProperty("user.home")), "elts_pv.dat");
        try (final PrintStream stream = new PrintStream(output1, "UTF-8")) {
            stream.println("# date pxn pyn pzn vxn vyn vzn pxk pyk pzk vxk vyk vzk");
            for (SpacecraftState numState : numStates) {
                for (SpacecraftState kepState : kepStates) {
                    if (numState.getDate().compareTo(kepState.getDate()) == 0) {
                        final double pxn = numState.getPVCoordinates().getPosition().getX();
                        final double pyn = numState.getPVCoordinates().getPosition().getY();
                        final double pzn = numState.getPVCoordinates().getPosition().getZ();
                        final double vxn = numState.getPVCoordinates().getVelocity().getX();
                        final double vyn = numState.getPVCoordinates().getVelocity().getY();
                        final double vzn = numState.getPVCoordinates().getVelocity().getZ();
                        final double pxk = kepState.getPVCoordinates().getPosition().getX();
                        final double pyk = kepState.getPVCoordinates().getPosition().getY();
                        final double pzk = kepState.getPVCoordinates().getPosition().getZ();
                        final double vxk = kepState.getPVCoordinates().getVelocity().getX();
                        final double vyk = kepState.getPVCoordinates().getVelocity().getY();
                        final double vzk = kepState.getPVCoordinates().getVelocity().getZ();
                        stream.println(numState.getDate() + " " + pxn + " " + pyn + " " + pzn + " " + vxn + " " + vyn + " " + vzn + " " + pxk + " " + pyk + " " + pzk + " " + vxk + " " + vyk + " " + vzk);
                        break;
                    }
                }
            }
        }
        System.out.println("Results saved as file " + output1);
    } catch (OrekitException oe) {
        System.err.println(oe.getLocalizedMessage());
        System.exit(1);
    } catch (IOException ioe) {
        System.err.println(ioe.getLocalizedMessage());
        System.exit(1);
    }
}

Example 8 with AbstractIntegrator

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

the class DSSTPropagation method createDSSTProp.

/**
 * Set up the DSST Propagator
 *
 *  @param orbit initial orbit
 *  @param mass S/C mass (kg)
 *  @param initialIsOsculating if initial orbital elements are osculating
 *  @param outputIsOsculating if we want to output osculating parameters
 *  @param fixedStepSize step size for fixed step integrator (s)
 *  @param minStep minimum step size, if step is not fixed (s)
 *  @param maxStep maximum step size, if step is not fixed (s)
 *  @param dP position tolerance for step size control, if step is not fixed (m)
 *  @param shortPeriodCoefficients list of short periodic coefficients
 *  to output (null means no coefficients at all, empty list means all
 *  possible coefficients)
 *  @throws OrekitException
 */
private DSSTPropagator createDSSTProp(final Orbit orbit, final double mass, final boolean initialIsOsculating, final boolean outputIsOsculating, final double fixedStepSize, final double minStep, final double maxStep, final double dP, final List<String> shortPeriodCoefficients) throws OrekitException {
    AbstractIntegrator integrator;
    if (fixedStepSize > 0.) {
        integrator = new ClassicalRungeKuttaIntegrator(fixedStepSize);
    } else {
        final double[][] tol = DSSTPropagator.tolerances(dP, orbit);
        integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
        ((AdaptiveStepsizeIntegrator) integrator).setInitialStepSize(10. * minStep);
    }
    DSSTPropagator dsstProp = new DSSTPropagator(integrator, !outputIsOsculating);
    dsstProp.setInitialState(new SpacecraftState(orbit, mass), initialIsOsculating);
    dsstProp.setSelectedCoefficients(shortPeriodCoefficients == null ? null : new HashSet<String>(shortPeriodCoefficients));
    return dsstProp;
}

Example 9 with AbstractIntegrator

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

the class EphemerisMode 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));
        // Initial orbit parameters
        // 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 ascension of ascending node
        double raan = FastMath.toRadians(261);
        // mean anomaly
        double lM = 0;
        // 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
        // the step is ten seconds
        double stepSize = 10;
        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());
    }
}

Example 10 with AbstractIntegrator

use of org.hipparchus.ode.AbstractIntegrator 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;
}