Propagation

The next 4 tutorials detail some elementary usages of the propagation package described in the propagation section of the library architecture documentation.

Steps management

Depending on the needs of the calling application, the propagator can provide intermediate states between the initial time and propagation target time, or it can provide only the final state.

Final state only

This method is used when the user wants to completely control the evolution of time. The application gives a target time and no step handlers at all

Let's define the EME2000 inertial frame as reference frame:

Frame inertialFrame = FramesFactory.getEME2000();

Let's set up an initial state as:

TimeScale utc = TimeScalesFactory.getUTC();
AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, utc);

double mu =  3.986004415e+14;

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

The initial orbit is defined as a KeplerianOrbit. More details on the orbit representation can be found in the orbits section of the library architecture documentation.

Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngleType.MEAN,
                                        inertialFrame, initialDate, mu);

We choose to use a very simple KeplerianPropagator to compute basic Keplerian motion. It could be any of the available propagators.

KeplerianPropagator kepler = new KeplerianPropagator(initialOrbit);

Finally, the propagation features of duration and step size are defined and a propagation loop is performed in order to print the results at each step:

double duration = 600.;
AbsoluteDate finalDate = initialDate.shiftedBy(duration);
double stepT = 60.;
int cpt = 1;
for (AbsoluteDate extrapDate = initialDate;
        extrapDate.compareTo(finalDate) <= 0;
        extrapDate = extrapDate.shiftedBy(stepT))  {
    SpacecraftState currentState = kepler.propagate(extrapDate);
    System.out.format(Locale.US, "step %2d %s %s%n",
                        cpt++, currentState.getDate(), currentState.getOrbit());
}

The printed results are shown below

step  1 2004-01-01T23:30:00.000 Keplerian parameters: {a: 2.4396159E7; e: 0.72831215; i: 7.0; pa: 180.0; raan: 261.0; v: 0.0;}
step  2 2004-01-01T23:31:00.000 Keplerian parameters: {a: 2.4396159E7; e: 0.72831215; i: 7.0; pa: 180.0; raan: 261.0; v: 5.281383633573694;}
step  3 2004-01-01T23:32:00.000 Keplerian parameters: {a: 2.4396159E7; e: 0.72831215; i: 7.0; pa: 180.0; raan: 261.0; v: 10.525261309411585;}
step  4 2004-01-01T23:33:00.000 Keplerian parameters: {a: 2.4396159E7; e: 0.72831215; i: 7.0; pa: 180.0; raan: 261.0; v: 15.695876839823306;}
step  5 2004-01-01T23:34:00.000 Keplerian parameters: {a: 2.4396159E7; e: 0.72831215; i: 7.0; pa: 180.0; raan: 261.0; v: 20.76077035517381;}
step  6 2004-01-01T23:35:00.000 Keplerian parameters: {a: 2.4396159E7; e: 0.72831215; i: 7.0; pa: 180.0; raan: 261.0; v: 25.691961427063767;}
step  7 2004-01-01T23:36:00.000 Keplerian parameters: {a: 2.4396159E7; e: 0.72831215; i: 7.0; pa: 180.0; raan: 261.0; v: 30.466680460539763;}
step  8 2004-01-01T23:37:00.000 Keplerian parameters: {a: 2.4396159E7; e: 0.72831215; i: 7.0; pa: 180.0; raan: 261.0; v: 35.06763907945756;}
step  9 2004-01-01T23:38:00.000 Keplerian parameters: {a: 2.4396159E7; e: 0.72831215; i: 7.0; pa: 180.0; raan: 261.0; v: 39.48289615024968;}
step 10 2004-01-01T23:39:00.000 Keplerian parameters: {a: 2.4396159E7; e: 0.72831215; i: 7.0; pa: 180.0; raan: 261.0; v: 43.70541689946282;}
step 11 2004-01-01T23:40:00.000 Keplerian parameters: {a: 2.4396159E7; e: 0.72831215; i: 7.0; pa: 180.0; raan: 261.0; v: 47.732436294590705;}

The complete code for this example can be found in the source tree of the tutorials, in file src/main/java/org/orekit/tutorials/propagation/KeplerianPropagation.java.

Ephemeris generation

This second mode may be used when the user needs random access to the orbit state at any time between some initial and final dates.

As in the first example above, let's first define some initial state as:

// Inertial frame
Frame inertialFrame = FramesFactory.getEME2000();
// Initial date
TimeScale utc = TimeScalesFactory.getUTC();
AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, utc);
// Central attraction coefficient
double mu =  3.986004415e+14;
// Initial orbit
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 ascention of ascending node
double lM = 0;                          // mean anomaly
Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngleType.MEAN,
                                        inertialFrame, initialDate, mu);
// Initial state definition
SpacecraftState initialState = new SpacecraftState(initialOrbit);

Here we use a simple NumericalPropagator, without perturbation, based on a classical fixed step Runge-Kutta integrator provided by the underlying Hipparchus library.

double stepSize = 10;
AbstractIntegrator integrator = new ClassicalRungeKuttaIntegrator(stepSize);
NumericalPropagator propagator = new NumericalPropagator(integrator);

The initial state is set for this propagator:

propagator.setInitialState(initialState);

Then, the propagator's ephemeris generator is used:

final EphemerisGenerator generator = propagator.getEphemerisGenerator();

And the propagation is performed for a given duration.

SpacecraftState finalState = propagator.propagate(initialDate.shiftedBy(6000));

This finalState can be used for anything, to be printed, for example as shown below:

 Numerical propagation :
  Final date : 2004-01-02T01:10:00.000Z
  equinoctial parameters: {a: 2.4396159E7;
                           ex: 0.11393312156755062; ey: 0.719345418868777;
                           hx: -0.009567941763699867; hy: -0.06040960680288257;
                           lv: 583.1250344407331;}

Throughout the propagation, intermediate states are stored within an ephemeris which can now be recovered with a single call:

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

The ephemeris is defined as a BoundedPropagator, which means that it is valid only between the propagation initial and final dates. The code above give the following result:

Ephemeris defined from 2004-01-01T23:30:00.000Z to 2004-01-02T01:10:00.000

Between these dates, the ephemeris can be used as any propagator to propagate the orbital state towards any intermediate date with a single call:

SpacecraftState intermediateState = ephemeris.propagate(intermediateDate);

Here are results obtained with intermediate dates set to 3000 seconds after start date and to exactly the final date:

Ephemeris propagation :
 date :  2004-01-02T00:20:00.000Z
 equinoctial parameters: {a: 2.4396159E7;
                          ex: 0.11393312156755062; ey: 0.719345418868777;
                          hx: -0.009567941763699867; hy: -0.06040960680288257;
                          lv: 559.0092657655282;}

 date :  2004-01-02T01:10:00.000Z
 equinoctial parameters: {a: 2.4396159E7;
                          ex: 0.11393312156755062; ey: 0.719345418868777;
                          hx: -0.009567941763699867; hy: -0.06040960680288257;
                          lv: 583.1250344407331;}

The following shows the error message we get when we try to use a date outside of the ephemeris range (in this case, the intermediate date was set to 1000 seconds before ephemeris start):

out of range date for ephemerides: 2004-01-01T23:13:20.000Z is 1.0E3 s before [2004-01-01T23:30:00.000Z, 2004-01-02T01:10:00.000Z]

The complete code for this example can be found in the source tree of the library, in file src/main/java/org/orekit/tutorials/propagation/EphemerisMode.java.

Intermediate states

The third mode is the most frequently used and is the recommended way to use propagation. It is intended for the use case of performing some simulation over a single time range with some sampling rate, from start to finish, with some computation performed at regular steps. A naive approach would be to implement a time loop in user code with some small time step and to call propagate within this time loop with a target time set just to the next time step, i.e. to use the first mode described above (dubbed Final State Only), but within a user-controlled loop. This naive approach is not recommended, as it increases the burden on user code (users have to manage a time loop and often have to manage preservation of some meta-data from one iteration to the next) and as it interacts badly with adaptive step size control in numerical propagators.

The recommended way to perform a computation for intermediate states throughout a single time range is to have the intermediate computation encapsulated within a dedicated custom object that implements either OrekitStepHandler or OrekitFixedStephHandler, to register it to the propagator and to run the propagator just once (without a loop), for the complete time range. The propagator will take care of the time loop, it will preserve the intermediate states of the event handlers (see next section), and it will be able to decorrelate the propagation step size from the sampling step size needed for the intermediate computation, hence allowing adaptive step size control to be used to its full potential.

Step handlers can be used with all propagators, analytical, semi-analytical or numerical. We could use the KeplerianPropagator as in the previous examples, or we can use something more accurate. Here, we will use a NumericalPropagator. Setting up a numerical propagator is much more involved than a Keplerian one. Users have to set up an ODE integrator, a propagator and add force models. If some of the force models correspond to surface forces, like atmospheric drag or radiation pressure, then the shape of the spacecraft and the attitude should also be defined.

This preparation phase can be daunting, and it is not the goal of this tutorial to focus on this, so we will rely here on builder methods that are shared among several tutorials to create everything that is needed for setting up a numerical propagator. These builder methods are configured using yaml files. The numerical-propagation.yaml file is an example of such a file, with a fairly complete set of force model. Users can alter this file to change the force models configuration if they want to use this tutorial for different propagation studies. When the tutorial is run without command line arguments, it uses the numerical-propagation.yaml example file from the src/main/resources folder as its input file. When it is run with an argument, it uses this argument as the input file name and generates its output file in the same folder as input, replacing the file extension by .out. This allows this tutorial to also be used as a standalone reference propagator to generate accurate positions.

// Earth model
final OneAxisEllipsoid earth = OrekitBuilder.createBody(inputData.getBody());

// gravity field
final TutorialGravity gravity = inputData.getPropagator().getForceModels().getGravity();
    final NormalizedSphericalHarmonicsProvider normalized =
        GravityFieldFactory.getConstantNormalizedProvider(gravity.getDegree(),
                                                            gravity.getOrder(),
                                                            new AbsoluteDate(inputData.getOrbit().getDate(),
                                                                            TimeScalesFactory.getUTC()));

// Initial orbit
final Orbit orbit = OrekitBuilder.createOrbit(inputData.getOrbit(), normalized.getMu());

// Propagator
final OrbitType type = OrbitType.CARTESIAN;
final ODEIntegrator integrator = OrekitBuilder.createIntegrator(inputData.getPropagator().getIntegrator(),
                                                                orbit, type, -1.0, 0.001, 1000, 0.001);
final AttitudeProvider attitudeProvider = inputData.
                                            getSatellite().
                                            getAttitude().
                                            getProvider(orbit.getFrame(), earth);
final NumericalPropagator propagator =
    OrekitBuilder.createPropagator(orbit, inputData.getSatellite().getMass(),
                                    integrator, attitudeProvider, type,
                                    normalized, earth.getBodyFrame(),
                                    inputData.getPropagator().getForceModels());

Once the propagator has been set up, we register custom TutorialStepHandler that encapsulate our computation needed (here it will just be printing the intermediate states to a file). This object implements both the OrekitFixedStepHandler interface and the Autocloseable interface. The first interface is what allows us to pass this object to the propagator, which will perform the time loop and will call the handleStep method at each fixed step (despite the underlying integrator uses variable step size). The second interface allows us to put the handler in a try-with-resources statement, thus ensuring the output file is properly closed when the program ends.

The really important thing to understand in this example is that there are no explicit time loops. The propagate method is called only once with the target end time. It is the propagator itself that will manage the time loop and will call the step handler.

// Set up a step handler that will autoclose at propagation end
try (TutorialStepHandler stepHandler = new TutorialStepHandler(earth, output)) {

    propagator.getMultiplexer().add(inputData.getOutputStep(), stepHandler);

    // Propagate for the specified duration in one sweep
    final AbsoluteDate targetDate = orbit.getDate().shiftedBy(inputData.getDuration());
    propagator.propagate(targetDate);

}

Our TutorialStepHandler implements the OrekitFixedStepHandler interface, and mainly its handleStep method, as follows (note that the init and finish methods can be omitted if not needed, as they have default implementations that do nothing):

private static class TutorialStepHandler implements OrekitFixedStepHandler, AutoCloseable {

    /** Earth model. */
    private final OneAxisEllipsoid earth;

    /** UTC time scale. */
    private final TimeScale utc;

    /** Output file. */
    private File outputFile;

    /** Output stream. */
    private PrintStream out;

    /** Entries counter. */
    private int count;

    /** Simple constructor.
     * @param earth Earth model
     * @param outputFile output file
     */
    public TutorialStepHandler(final OneAxisEllipsoid earth, final File outputFile) throws IOException{
        this.earth      = earth;
        this.utc        = TimeScalesFactory.getUTC();
        this.outputFile = outputFile;
        this.out        = new PrintStream(outputFile);
    }

    /** {@inheritDoc} */
    @Override
    public void init(final SpacecraftState s0, final AbsoluteDate t, final double step) {
        out.println("         date" +
                    "                PX inert          PY inert          PZ inert" +
                    "       VX inert          VY inert          VZ inert" +
                    "       PX Earth          PY Earth         PZ Earth" +
                    "       VX Earth          VY Earth         VZ Earth" +
                    "     latitude   longitude   altitude");
        count = 0;
    }

    /** {@inheritDoc} */
    @Override
    public void handleStep(final SpacecraftState currentState) {
        final PVCoordinates pvInert     = currentState.getPVCoordinates(currentState.getFrame());
        final Transform     inert2Earth = currentState.getFrame().
                                            getTransformTo(earth.getBodyFrame(), currentState.getDate());
        final PVCoordinates pvEarth     = inert2Earth.transformPVCoordinates(pvInert);
        final GeodeticPoint gp          = earth.transform(pvEarth.getPosition(), earth.getBodyFrame(), currentState.getDate());
        out.format(Locale.US,
                    "%s %16.6f %16.6f %16.6f %16.9f %16.9f %16.9f %16.6f %16.6f %16.6f %16.9f %16.9f %16.9f %8.3f   %8.3f %13.3f%n",
                    currentState.getDate().toStringWithoutUtcOffset(utc, 3),
                    pvInert.getPosition().getX(), pvInert.getPosition().getY(), pvInert.getPosition().getZ(),
                    pvInert.getVelocity().getX(), pvInert.getVelocity().getY(), pvInert.getVelocity().getZ(),
                    pvEarth.getPosition().getX(), pvEarth.getPosition().getY(), pvEarth.getPosition().getZ(),
                    pvEarth.getVelocity().getX(), pvEarth.getVelocity().getY(), pvEarth.getVelocity().getZ(),
                    FastMath.toDegrees(gp.getLatitude()),
                    FastMath.toDegrees(gp.getLongitude()),
                    gp.getAltitude());
            ++count;
    }

    /** {@inheritDoc} */
    @Override
    public void finish(final SpacecraftState finalState) {
        out.flush();
        System.out.format(Locale.US, "%d steps have been written to %s%n", count, outputFile);
    }

    /** {@inheritDoc} */
    @Override
    public void close() {
        if (out != null) {
            out.close();
            out = null;
        }
    }

}

This example is representative of a classical design for Orekit-based applications: they start with a (sometimes big) initialization phase to set up a propagator, they call once a propagate method, and they stop when it completes. The core of the computation is done in the step handler.

When run, this tutorial will create and fill-up a text file and will print a single line on standard output at the end:

11 steps have been written to <some path>/numerical-propagation.out

The numerical-propagation.out file will have the following content:

         date                PX inert          PY inert          PZ inert       VX inert          VY inert          VZ inert       PX Earth          PY Earth         PZ Earth       VX Earth          VY Earth         VZ Earth     latitude   longitude   altitude
2004-01-01T23:30:00.000   1036869.533073   6546536.584961         0.000000  -9994.355195123   1582.950353984  -1242.449124591   6473251.708515  -1424518.060793       557.420126   2070.888720149   9409.982052653  -1246.049874321    0.005    -12.411    250002.988
2004-01-01T23:31:00.000    435145.161985   6625301.741998    -74485.459407 -10054.836792382   1041.961176108  -1239.380166290   6583989.003788   -857387.807832    -74145.281442   1619.437607949   9487.205015972  -1243.017818595   -0.644     -7.419    261859.849
2004-01-01T23:32:00.000   -168717.835613   6671585.513444   -148605.100527 -10065.726801835    501.642499579  -1230.292474265   6667550.869618   -286915.330004   -148483.836008   1166.313921963   9521.389737155  -1233.948925583   -1.283     -2.464    297246.463
2004-01-01T23:33:00.000   -771776.831684   6685676.803337   -222006.598213 -10028.409855529    -29.803599829  -1215.513022549   6724063.110923    284333.350240   -222104.837380    718.991523013   9513.344152793  -1219.170449290   -1.902      2.421    355622.390
2004-01-01T23:34:00.000  -1371221.915221   6668336.314515   -294362.929237  -9945.881145640   -544.926289125  -1195.551423695   6754083.921000    853886.550364   -294680.216735    284.339810432   9465.368753689  -1199.192890149   -2.494      7.205    436124.209
2004-01-01T23:35:00.000  -1964463.458103   6620739.495315   -365381.974160  -9822.426829536  -1037.402744231  -1171.052794413   6758556.199584   1419448.860800   -365916.877998   -131.736933547   9380.985628789  -1174.662754438   -3.052     11.861    537616.316
2004-01-01T23:36:00.000  -2549197.933761   6544401.070930   -434813.038930  -9663.228116929  -1502.273535591  -1142.743812774   6738741.942381   1978964.263401   -435563.245503   -524.534607475   9264.590760551  -1146.308473580   -3.570     16.366    658752.656
2004-01-01T23:37:00.000  -3123449.055188   6441089.984589   -502450.153584  -9473.943493195  -1936.025516757  -1111.380376300   6696146.633011   2530657.303882   -503412.580032   -890.645568268   9121.078227693  -1114.887912158   -4.047     20.703    798042.152
2004-01-01T23:38:00.000  -3685584.661818   6312743.662014   -568132.440681  -9260.326965621  -2336.536385405  -1077.703086651   6632441.853256   3073052.238722   -569303.355755  -1227.925763561   8955.490237875  -1081.143711877   -4.479     24.860    953911.812
2004-01-01T23:39:00.000  -4234312.805188   6161388.817822   -631742.123314  -9027.921390813  -2702.921583535  -1042.403313474   6549392.930236   3604972.950251   -633117.268992  -1535.368391262   8772.731212215  -1045.769239383   -4.868     28.830   1124763.063
2004-01-01T23:40:00.000  -4768662.074253   5989073.201508   -693200.830216  -8781.838911171  -3035.327067886  -1006.100279371   6448795.893118   4125528.045899   -694775.538010  -1812.923846256   8577.359004766  -1009.385585490   -5.214     32.609   1309018.478

The complete code for this example can be found in the source tree of the library, in file src/main/java/org/orekit/tutorials/propagation/NumericalPropagation.java.

Events management

This tutorial aims to demonstrate the power and simplicity of the event-handling mechanism.

One needs to check the visibility between a satellite and a ground station during some time range.

We will use, and extend, the predefined ElevationDetector to perform the task.

First, let's set up an initial state for the satellite defined by a position and a velocity at one date in some inertial frame.

Vector3D position  = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
Vector3D velocity  = new Vector3D(505.8479685, 942.7809215, 7435.922231);
PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
AbsoluteDate initialDate =
    new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, TimeScalesFactory.getUTC());
Frame inertialFrame = FramesFactory.getEME2000();

We also need to set the central attraction coefficient to define the initial orbit as a KeplerianOrbit.

double mu =  3.986004415e+14;
Orbit initialOrbit =
    new KeplerianOrbit(pvCoordinates, inertialFrame, initialDate, mu);

More details on the orbit representation can be found in the orbits section of the library architecture documentation.

As a propagator, we consider a KeplerianPropagator to compute the simple Keplerian motion. It could be more elaborate without modifying the general purpose of this tutorial.

Propagator kepler = new KeplerianPropagator(initialOrbit);

Then, let's define the ground station by its coordinates as a GeodeticPoint:

double longitude = FastMath.toRadians(45.);
double latitude  = FastMath.toRadians(25.);
double altitude  = 0.;
GeodeticPoint station1 = new GeodeticPoint(latitude, longitude, altitude);

And let's associate to it a TopocentricFrame related to a BodyShape in some terrestrial frame.

Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                                        Constants.WGS84_EARTH_FLATTENING,
                                        earthFrame);
TopocentricFrame sta1Frame = new TopocentricFrame(earth, station1, "station1");

More details on BodyShape and GeodeticPoint can be found in the bodies section of the library architecture documentation.
More details on TopocentricFrame can be found in the frames section of the library architecture documentation.

An EventDetector is defined as a an ElevationDetector with constant limit elevation and a dedicated handler, which in this simple case can be written inline as a lambda function that implements the eventOccurred method from the EventHandler interface:

double maxcheck  = 60.0;
double threshold =  0.001;
double elevation = FastMath.toRadians(5.);
EventDetector sta1Visi =
    new ElevationDetector(maxcheck, threshold, sta1Frame).
    withConstantElevation(elevation).
    withHandler((s, detector, increasing) -> {
                    System.out.println(" Visibility on " +
                                        ((ElevationDetector) detector).getTopocentricFrame().getName() +
                                        (increasing ? " begins at " : " ends at ") +
                                        s.getDate().toStringWithoutUtcOffset(utc, 3));
                    return increasing ? Action.CONTINUE : Action.STOP;
                });

This EventDetector is added to the propagator:

kepler.addEventDetector(sta1Visi);

The propagator will propagate from the initial date to the first raising or for the fixed duration according to the behavior implemented in the VisibilityHandler class.

SpacecraftState finalState =
    kepler.propagate(new AbsoluteDate(initialDate, 1500.));
System.out.println(" Final state : " + finalState.getDate().durationFrom(initialDate));

The printed result is shown below. We can see that the propagation stopped just after detecting the setting:

Visibility on station1 begins at 2004-01-01T23:31:52.098
Visibility on station1 ends at 2004-01-01T23:42:48.850
Final state duration from initial date (s) : 768.850267132939

The complete code for this example can be found in the source tree of the tutorials, in file src/main/java/org/orekit/tutorials/propagation/VisibilityCheck.java.