Examples with EcksteinHechlerPropagator org.orekit.propagation.analytical.EcksteinHechlerPropagator used on opensource projects

Example 31 with EcksteinHechlerPropagator

use of org.orekit.propagation.analytical.EcksteinHechlerPropagator in project Orekit by CS-SI.

the class EquinoctialOrbitTest method doTestInterpolation.

private void doTestInterpolation(boolean useDerivatives, double shiftErrorWithin, double interpolationErrorWithin, double shiftErrorSlightlyPast, double interpolationErrorSlightlyPast, double shiftErrorFarPast, double interpolationErrorFarPast) throws OrekitException {
    final double ehMu = 3.9860047e14;
    final double ae = 6.378137e6;
    final double c20 = -1.08263e-3;
    final double c30 = 2.54e-6;
    final double c40 = 1.62e-6;
    final double c50 = 2.3e-7;
    final double c60 = -5.5e-7;
    final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
    final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
    final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
    final EquinoctialOrbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), date, ehMu);
    EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
    // set up a 5 points sample
    List<Orbit> sample = new ArrayList<Orbit>();
    for (double dt = 0; dt < 300.0; dt += 60.0) {
        Orbit orbit = propagator.propagate(date.shiftedBy(dt)).getOrbit();
        if (!useDerivatives) {
            // remove derivatives
            double[] stateVector = new double[6];
            orbit.getType().mapOrbitToArray(orbit, PositionAngle.TRUE, stateVector, null);
            orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngle.TRUE, orbit.getDate(), orbit.getMu(), orbit.getFrame());
        }
        sample.add(orbit);
    }
    // well inside the sample, interpolation should be much better than Keplerian shift
    double maxShiftError = 0;
    double maxInterpolationError = 0;
    for (double dt = 0; dt < 241.0; dt += 1.0) {
        AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
        Vector3D shifted = initialOrbit.shiftedBy(dt).getPVCoordinates().getPosition();
        Vector3D interpolated = initialOrbit.interpolate(t, sample).getPVCoordinates().getPosition();
        Vector3D propagated = propagator.propagate(t).getPVCoordinates().getPosition();
        maxShiftError = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
        maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
    }
    Assert.assertEquals(shiftErrorWithin, maxShiftError, 0.01 * shiftErrorWithin);
    Assert.assertEquals(interpolationErrorWithin, maxInterpolationError, 0.01 * interpolationErrorWithin);
    // slightly past sample end, interpolation should quickly increase, but remain reasonable
    maxShiftError = 0;
    maxInterpolationError = 0;
    for (double dt = 240; dt < 300.0; dt += 1.0) {
        AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
        Vector3D shifted = initialOrbit.shiftedBy(dt).getPVCoordinates().getPosition();
        Vector3D interpolated = initialOrbit.interpolate(t, sample).getPVCoordinates().getPosition();
        Vector3D propagated = propagator.propagate(t).getPVCoordinates().getPosition();
        maxShiftError = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
        maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
    }
    Assert.assertEquals(shiftErrorSlightlyPast, maxShiftError, 0.01 * shiftErrorSlightlyPast);
    Assert.assertEquals(interpolationErrorSlightlyPast, maxInterpolationError, 0.01 * interpolationErrorSlightlyPast);
    // far past sample end, interpolation should become really wrong
    // (in this test case, break even occurs at around 863 seconds, with a 3.9 km error)
    maxShiftError = 0;
    maxInterpolationError = 0;
    for (double dt = 300; dt < 1000; dt += 1.0) {
        AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
        Vector3D shifted = initialOrbit.shiftedBy(dt).getPVCoordinates().getPosition();
        Vector3D interpolated = initialOrbit.interpolate(t, sample).getPVCoordinates().getPosition();
        Vector3D propagated = propagator.propagate(t).getPVCoordinates().getPosition();
        maxShiftError = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
        maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
    }
    Assert.assertEquals(shiftErrorFarPast, maxShiftError, 0.01 * shiftErrorFarPast);
    Assert.assertEquals(interpolationErrorFarPast, maxInterpolationError, 0.01 * interpolationErrorFarPast);
}

Example 32 with EcksteinHechlerPropagator

use of org.orekit.propagation.analytical.EcksteinHechlerPropagator in project Orekit by CS-SI.

the class KeplerianOrbitTest method doTestInterpolation.

private void doTestInterpolation(boolean useDerivatives, double shiftPositionErrorWithin, double interpolationPositionErrorWithin, double shiftEccentricityErrorWithin, double interpolationEccentricityErrorWithin, double shiftPositionErrorSlightlyPast, double interpolationPositionErrorSlightlyPast, double shiftEccentricityErrorSlightlyPast, double interpolationEccentricityErrorSlightlyPast) throws OrekitException {
    final double ehMu = 3.9860047e14;
    final double ae = 6.378137e6;
    final double c20 = -1.08263e-3;
    final double c30 = 2.54e-6;
    final double c40 = 1.62e-6;
    final double c50 = 2.3e-7;
    final double c60 = -5.5e-7;
    final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
    final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
    final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
    final KeplerianOrbit initialOrbit = new KeplerianOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), date, ehMu);
    EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
    // set up a 5 points sample
    List<Orbit> sample = new ArrayList<Orbit>();
    for (double dt = 0; dt < 300.0; dt += 60.0) {
        Orbit orbit = propagator.propagate(date.shiftedBy(dt)).getOrbit();
        if (!useDerivatives) {
            // remove derivatives
            double[] stateVector = new double[6];
            orbit.getType().mapOrbitToArray(orbit, PositionAngle.TRUE, stateVector, null);
            orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngle.TRUE, orbit.getDate(), orbit.getMu(), orbit.getFrame());
        }
        sample.add(orbit);
    }
    // well inside the sample, interpolation should be slightly better than Keplerian shift
    // the relative bad behaviour here is due to eccentricity, which cannot be
    // accurately interpolated with a polynomial in this case
    double maxShiftPositionError = 0;
    double maxInterpolationPositionError = 0;
    double maxShiftEccentricityError = 0;
    double maxInterpolationEccentricityError = 0;
    for (double dt = 0; dt < 241.0; dt += 1.0) {
        AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
        Vector3D shiftedP = initialOrbit.shiftedBy(dt).getPVCoordinates().getPosition();
        Vector3D interpolatedP = initialOrbit.interpolate(t, sample).getPVCoordinates().getPosition();
        Vector3D propagatedP = propagator.propagate(t).getPVCoordinates().getPosition();
        double shiftedE = initialOrbit.shiftedBy(dt).getE();
        double interpolatedE = initialOrbit.interpolate(t, sample).getE();
        double propagatedE = propagator.propagate(t).getE();
        maxShiftPositionError = FastMath.max(maxShiftPositionError, shiftedP.subtract(propagatedP).getNorm());
        maxInterpolationPositionError = FastMath.max(maxInterpolationPositionError, interpolatedP.subtract(propagatedP).getNorm());
        maxShiftEccentricityError = FastMath.max(maxShiftEccentricityError, FastMath.abs(shiftedE - propagatedE));
        maxInterpolationEccentricityError = FastMath.max(maxInterpolationEccentricityError, FastMath.abs(interpolatedE - propagatedE));
    }
    Assert.assertEquals(shiftPositionErrorWithin, maxShiftPositionError, 0.01 * shiftPositionErrorWithin);
    Assert.assertEquals(interpolationPositionErrorWithin, maxInterpolationPositionError, 0.01 * interpolationPositionErrorWithin);
    Assert.assertEquals(shiftEccentricityErrorWithin, maxShiftEccentricityError, 0.01 * shiftEccentricityErrorWithin);
    Assert.assertEquals(interpolationEccentricityErrorWithin, maxInterpolationEccentricityError, 0.01 * interpolationEccentricityErrorWithin);
    // slightly past sample end, bad eccentricity interpolation shows up
    // (in this case, interpolated eccentricity exceeds 1.0 btween 1900
    // and 1910s, while semi-majaxis remains positive, so this is not
    // even a proper hyperbolic orbit...)
    maxShiftPositionError = 0;
    maxInterpolationPositionError = 0;
    maxShiftEccentricityError = 0;
    maxInterpolationEccentricityError = 0;
    for (double dt = 240; dt < 600; dt += 1.0) {
        AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
        Vector3D shiftedP = initialOrbit.shiftedBy(dt).getPVCoordinates().getPosition();
        Vector3D interpolatedP = initialOrbit.interpolate(t, sample).getPVCoordinates().getPosition();
        Vector3D propagatedP = propagator.propagate(t).getPVCoordinates().getPosition();
        double shiftedE = initialOrbit.shiftedBy(dt).getE();
        double interpolatedE = initialOrbit.interpolate(t, sample).getE();
        double propagatedE = propagator.propagate(t).getE();
        maxShiftPositionError = FastMath.max(maxShiftPositionError, shiftedP.subtract(propagatedP).getNorm());
        maxInterpolationPositionError = FastMath.max(maxInterpolationPositionError, interpolatedP.subtract(propagatedP).getNorm());
        maxShiftEccentricityError = FastMath.max(maxShiftEccentricityError, FastMath.abs(shiftedE - propagatedE));
        maxInterpolationEccentricityError = FastMath.max(maxInterpolationEccentricityError, FastMath.abs(interpolatedE - propagatedE));
    }
    Assert.assertEquals(shiftPositionErrorSlightlyPast, maxShiftPositionError, 0.01 * shiftPositionErrorSlightlyPast);
    Assert.assertEquals(interpolationPositionErrorSlightlyPast, maxInterpolationPositionError, 0.01 * interpolationPositionErrorSlightlyPast);
    Assert.assertEquals(shiftEccentricityErrorSlightlyPast, maxShiftEccentricityError, 0.01 * shiftEccentricityErrorSlightlyPast);
    Assert.assertEquals(interpolationEccentricityErrorSlightlyPast, maxInterpolationEccentricityError, 0.01 * interpolationEccentricityErrorSlightlyPast);
}

Example 33 with EcksteinHechlerPropagator

use of org.orekit.propagation.analytical.EcksteinHechlerPropagator in project Orekit by CS-SI.

the class EarthObservation 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));
        final SortedSet<String> output = new TreeSet<String>();
        // Initial state definition : date, orbit
        final AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, TimeScalesFactory.getUTC());
        final Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
        final Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
        final Orbit initialOrbit = new KeplerianOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), initialDate, Constants.EIGEN5C_EARTH_MU);
        // Attitudes sequence definition
        final AttitudeProvider dayObservationLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VVLH, RotationOrder.XYZ, FastMath.toRadians(20), FastMath.toRadians(40), 0);
        final AttitudeProvider nightRestingLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VVLH);
        final PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
        final PVCoordinatesProvider earth = CelestialBodyFactory.getEarth();
        final EventDetector dayNightEvent = new EclipseDetector(sun, 696000000., earth, Constants.WGS84_EARTH_EQUATORIAL_RADIUS).withHandler(new ContinueOnEvent<EclipseDetector>());
        final EventDetector nightDayEvent = new EclipseDetector(sun, 696000000., earth, Constants.WGS84_EARTH_EQUATORIAL_RADIUS).withHandler(new ContinueOnEvent<EclipseDetector>());
        final AttitudesSequence attitudesSequence = new AttitudesSequence();
        final AttitudesSequence.SwitchHandler switchHandler = new AttitudesSequence.SwitchHandler() {

            public void switchOccurred(AttitudeProvider preceding, AttitudeProvider following, SpacecraftState s) {
                if (preceding == dayObservationLaw) {
                    output.add(s.getDate() + ": switching to night law");
                } else {
                    output.add(s.getDate() + ": switching to day law");
                }
            }
        };
        attitudesSequence.addSwitchingCondition(dayObservationLaw, nightRestingLaw, dayNightEvent, false, true, 10.0, AngularDerivativesFilter.USE_R, switchHandler);
        attitudesSequence.addSwitchingCondition(nightRestingLaw, dayObservationLaw, nightDayEvent, true, false, 10.0, AngularDerivativesFilter.USE_R, switchHandler);
        if (dayNightEvent.g(new SpacecraftState(initialOrbit)) >= 0) {
            // initial position is in daytime
            attitudesSequence.resetActiveProvider(dayObservationLaw);
        } else {
            // initial position is in nighttime
            attitudesSequence.resetActiveProvider(nightRestingLaw);
        }
        // Propagator : consider the analytical Eckstein-Hechler model
        final Propagator propagator = new EcksteinHechlerPropagator(initialOrbit, attitudesSequence, Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS, Constants.EIGEN5C_EARTH_MU, Constants.EIGEN5C_EARTH_C20, Constants.EIGEN5C_EARTH_C30, Constants.EIGEN5C_EARTH_C40, Constants.EIGEN5C_EARTH_C50, Constants.EIGEN5C_EARTH_C60);
        // Register the switching events to the propagator
        attitudesSequence.registerSwitchEvents(propagator);
        propagator.setMasterMode(180.0, new OrekitFixedStepHandler() {

            public void init(final SpacecraftState s0, final AbsoluteDate t) {
            }

            public void handleStep(SpacecraftState currentState, boolean isLast) throws OrekitException {
                DecimalFormatSymbols angleDegree = new DecimalFormatSymbols(Locale.US);
                angleDegree.setDecimalSeparator('\u00b0');
                DecimalFormat ad = new DecimalFormat(" 00.000;-00.000", angleDegree);
                // the Earth position in spacecraft frame should be along spacecraft Z axis
                // during nigthtime and away from it during daytime due to roll and pitch offsets
                final Vector3D earth = currentState.toTransform().transformPosition(Vector3D.ZERO);
                final double pointingOffset = Vector3D.angle(earth, Vector3D.PLUS_K);
                // the g function is the eclipse indicator, its an angle between Sun and Earth limb,
                // positive when Sun is outside of Earth limb, negative when Sun is hidden by Earth limb
                final double eclipseAngle = dayNightEvent.g(currentState);
                output.add(currentState.getDate() + " " + ad.format(FastMath.toDegrees(eclipseAngle)) + " " + ad.format(FastMath.toDegrees(pointingOffset)));
            }
        });
        // Propagate from the initial date for the fixed duration
        SpacecraftState finalState = propagator.propagate(initialDate.shiftedBy(12600.));
        // to make sure out of orders calls between step handler and event handlers don't mess things up
        for (final String line : output) {
            System.out.println(line);
        }
        System.out.println("Propagation ended at " + finalState.getDate());
    } catch (OrekitException oe) {
        System.err.println(oe.getMessage());
    }
}

Example 34 with EcksteinHechlerPropagator

use of org.orekit.propagation.analytical.EcksteinHechlerPropagator in project Orekit by CS-SI.

the class Frames3 method main.

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 state definition :
        // ==========================
        // Date
        // ****
        AbsoluteDate initialDate = new AbsoluteDate(2003, 4, 7, 10, 55, 21.575, TimeScalesFactory.getUTC());
        // Orbit
        // *****
        // The Sun is in the orbital plane for raan ~ 202
        // gravitation coefficient
        double mu = 3.986004415e+14;
        // inertial frame
        Frame eme2000 = FramesFactory.getEME2000();
        Orbit orbit = new CircularOrbit(7178000.0, 0.5e-4, -0.5e-4, FastMath.toRadians(50.), FastMath.toRadians(220.), FastMath.toRadians(5.300), PositionAngle.MEAN, eme2000, initialDate, mu);
        // Attitude laws
        // *************
        // Earth
        Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
        BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, earthFrame);
        // Target pointing attitude provider over satellite nadir at date, without yaw compensation
        NadirPointing nadirLaw = new NadirPointing(eme2000, earth);
        // Target pointing attitude provider with yaw compensation
        final PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
        YawSteering yawSteeringLaw = new YawSteering(eme2000, nadirLaw, sun, Vector3D.MINUS_I);
        // Propagator : Eckstein-Hechler analytic propagator
        Propagator propagator = new EcksteinHechlerPropagator(orbit, yawSteeringLaw, Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS, Constants.EIGEN5C_EARTH_MU, Constants.EIGEN5C_EARTH_C20, Constants.EIGEN5C_EARTH_C30, Constants.EIGEN5C_EARTH_C40, Constants.EIGEN5C_EARTH_C50, Constants.EIGEN5C_EARTH_C60);
        // Let's write the results in a file in order to draw some plots.
        propagator.setMasterMode(10, new OrekitFixedStepHandler() {

            PrintStream out = null;

            public void init(SpacecraftState s0, AbsoluteDate t, double step) throws OrekitException {
                try {
                    File file = new File(System.getProperty("user.home"), "XYZ.dat");
                    System.out.println("Results written to file: " + file.getAbsolutePath());
                    out = new PrintStream(file, "UTF-8");
                    out.println("#time X Y Z Wx Wy Wz");
                } catch (IOException ioe) {
                    throw new OrekitException(ioe, LocalizedCoreFormats.SIMPLE_MESSAGE, ioe.getLocalizedMessage());
                }
            }

            public void handleStep(SpacecraftState currentState, boolean isLast) throws OrekitException {
                // get the transform from orbit/attitude reference frame to spacecraft frame
                Transform inertToSpacecraft = currentState.toTransform();
                // get the position of the Sun in orbit/attitude reference frame
                Vector3D sunInert = sun.getPVCoordinates(currentState.getDate(), currentState.getFrame()).getPosition();
                // convert Sun position to spacecraft frame
                Vector3D sunSat = inertToSpacecraft.transformPosition(sunInert);
                // and the spacecraft rotational rate also
                Vector3D spin = inertToSpacecraft.getRotationRate();
                // Lets calculate the reduced coordinates
                double sunX = sunSat.getX() / sunSat.getNorm();
                double sunY = sunSat.getY() / sunSat.getNorm();
                double sunZ = sunSat.getZ() / sunSat.getNorm();
                out.format(Locale.US, "%s %12.3f %12.3f %12.3f %12.7f %12.7f %12.7f%n", currentState.getDate(), sunX, sunY, sunZ, spin.getX(), spin.getY(), spin.getZ());
                if (isLast) {
                    out.close();
                }
            }
        });
        System.out.println("Running...");
        propagator.propagate(initialDate.shiftedBy(6000));
    } catch (OrekitException oe) {
        System.err.println(oe.getMessage());
    }
}