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, PositionAngle.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.

Intermediate states

In order to get intermediate states, users need to register some custom object implementation of the OrekitStepHandler or OrekitFixedStephHandler

As in the “Final state only” example above, let's define some initial state with:

// 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, PositionAngle.MEAN,
                                        inertialFrame, initialDate, mu);
// Initial state definition
SpacecraftState initialState = new SpacecraftState(initialOrbit);

Here we use a more sophisticated NumericalPropagator based on an adaptive step integrator provided by the underlying Hipparchus library, but it doesn't matter which integrator is used.

// Adaptive step integrator
// with a minimum step of 0.001 and a maximum step of 1000
double minStep = 0.001;
double maxstep = 1000.0;
double positionTolerance = 10.0;
OrbitType propagationType = OrbitType.KEPLERIAN;
double[][] tolerances =
    NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);
AdaptiveStepsizeIntegrator integrator =
    new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);

We set up the integrator, and force it to use Keplerian parameters for propagation.

NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(propagationType);

A force model, reduced here to a single perturbing gravity field, is taken into account. More details on force models can be found in the forces section of the library architecture documentation.

NormalizedSphericalHarmonicsProvider provider =
    GravityFieldFactory.getNormalizedProvider(10, 10);
ForceModel holmesFeatherstone =
    new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010,
                                                                true),
                                          provider)

This force model is simply added to the propagator:

propagator.addForceModel(holmesFeatherstone);

Then, the initial state is set for the propagator:

propagator.setInitialState(initialState);

The propagator configuration is completed with a fixed step TutorialStepHandler which implements the interface OrekitFixedStepHandler in order to fulfill the handleStep method to be called within the loop. For the purpose of this tutorial, the handleStep method will just print the current state at the moment.

propagator.getMultiplexer().add(60., new TutorialStepHandler());

Finally, the propagator is just asked to propagate, from the initial state, for a given duration.

SpacecraftState finalState =
    propagator.propagate(new AbsoluteDate(initialDate, 630.));

Clearly, with a few lines of code, the main application delegates to the propagator the care of handling regular outputs through a variable step integration loop.

All that is needed is to derive some class from the interface OrekitFixedStepHandler to realize a handleStep method, as follows:

private static class TutorialStepHandler implements OrekitFixedStepHandler {

    public void init(final SpacecraftState s0, final AbsoluteDate t, final double step) {
        System.out.println("          date                a           e" +
                           "           i         \u03c9          \u03a9" +
                           "          \u03bd");
    }

    public void handleStep(SpacecraftState currentState) {
        KeplerianOrbit o = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(currentState.getOrbit());
        System.out.format(Locale.US, "%s %12.3f %10.8f %10.6f %10.6f %10.6f %10.6f%n",
                          currentState.getDate(),
                          o.getA(), o.getE(),
                          FastMath.toDegrees(o.getI()),
                          FastMath.toDegrees(o.getPerigeeArgument()),
                          FastMath.toDegrees(o.getRightAscensionOfAscendingNode()),
                          FastMath.toDegrees(o.getTrueAnomaly()));
    }

    public void finish(final SpacecraftState finalState) {
        System.out.println("this was the last step ");
        System.out.println();
    }

}

The printed results are shown below:

          date                a           e           i         ω          Ω          ν
2004-01-01T23:30:00.000 24396159.000 0.72831215   7.000000 180.000000 261.000000   0.000000
2004-01-01T23:31:00.000 24395672.948 0.72830573   6.999927 180.010992 260.999966   5.270439
2004-01-01T23:32:00.000 24394149.110 0.72828580   6.999758 180.021885 260.999767  10.503710
2004-01-01T23:33:00.000 24391676.555 0.72825352   6.999506 180.032557 260.999281  15.664370
2004-01-01T23:34:00.000 24388396.501 0.72821073   6.999190 180.042889 260.998422  20.720220
2004-01-01T23:35:00.000 24384487.471 0.72815976   6.998827 180.052785 260.997139  25.643463
2004-01-01T23:36:00.000 24380142.510 0.72810308   6.998434 180.062182 260.995416  30.411433
2004-01-01T23:37:00.000 24375546.779 0.72804310   6.998029 180.071043 260.993268  35.006866
2004-01-01T23:38:00.000 24370861.668 0.72798191   6.997624 180.079354 260.990738  39.417786
2004-01-01T23:39:00.000 24366216.959 0.72792119   6.997234 180.087112 260.987893  43.637076
2004-01-01T23:40:00.000 24361709.735 0.72786221   6.996868 180.094320 260.984809  47.661865
this was the last step

Final state:
2004-01-01T23:40:30.000 24359529.816 0.72783366   6.996697 180.097720 260.983203  49.601407

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.

Ephemeris generation

This third 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 two examples 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, PositionAngle.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.000
  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.000 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.000
 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.000
 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.000, [2004-01-01T23:30:00.000, 2004-01-02T01:10:00.000]

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.

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 " +
                                       detector.getTopocentricFrame().getName() +
                                       (increasing ? " begins at " : " ends at ") +
                                       s.getDate());
                    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 : 768.850266477389

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.

Orekit Tutorials

Markdown Tutorials