## PropagationTutorial » History » Version 4

« Previous - Version 4/10 (diff) - Next » - Current version
Michael Turner, 2012-04-26 16:28
minor fixes

h1. Propagation

{{toc}}

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

h2. Propagation modes

Three different operational modes are available for all propagators. They are mutually exclusive.

h3. Slave mode

This is the default mode. It doesn't need to be explicitly set, but it can be, to recover from any other mode.

This tutorial shows how to propagate from an initial state to a target date.

In this case, the calling application coordinates all the tasks, the propagator just propagates.

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

Frame inertialFrame = FramesFactory.getEME2000();

Let's set up an initial state as:

• a date in some time scale
• a central attraction coefficient
• an orbit defined by its keplerian parameters.

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 = Math.toRadians(7); // inclination
double omega = Math.toRadians(180); // perigee argument
double raan = Math.toRadians(261); // right ascension of ascending node
double lM = 0; // mean anomaly

The initial orbit is defined as a KeplerianOrbit.
More detail 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);

Let's set the propagator to slave mode for the purpose of this tutorial, but keep in mind that it can be omitted here as it is the default mode.

kepler.setSlaveMode();

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 = new AbsoluteDate(initialDate, duration, utc);
double stepT = 60.;
int cpt = 1;
AbsoluteDate extrapDate = initialDate;
while (extrapDate.compareTo(finalDate) <= 0) {
SpacecraftState currentState = kepler.propagate(extrapDate);
System.out.println("step " + cpt++);
System.out.println(" time : " + currentState.getDate());
System.out.println(" " + currentState.getOrbit());
extrapDate = new AbsoluteDate(extrapDate, stepT, utc);
}

The printed results are shown below (lines have been wrapped for better presentation)

step 1
time : 2004-01-01T23:30:00.000
equinoctial parameters: {a: 2.4396159E7;
ex: 0.11393312156755062; ey: 0.719345418868777;
hx: -0.009567941763699867; hy: -0.06040960680288257;
lv: 441.0;}
step 2
time : 2004-01-01T23:31:00.000
equinoctial parameters: {a: 2.4396159E7;
ex: 0.11393312156755062; ey: 0.719345418868777;
hx: -0.009567941763699867; hy: -0.06040960680288257;
lv: 446.2813836335737;}
step 3
time : 2004-01-01T23:32:00.000
equinoctial parameters: {a: 2.4396159E7;
ex: 0.11393312156755062; ey: 0.719345418868777;
hx: -0.009567941763699867; hy: -0.06040960680288257;
lv: 451.5252613094112;}
step 4
time : 2004-01-01T23:33:00.000
equinoctial parameters: {a: 2.4396159E7;
ex: 0.11393312156755062; ey: 0.719345418868777;
hx: -0.009567941763699867; hy: -0.06040960680288257;
lv: 456.6958768398235;}
step 5
time : 2004-01-01T23:34:00.000
equinoctial parameters: {a: 2.4396159E7;
ex: 0.11393312156755062; ey: 0.719345418868777;
hx: -0.009567941763699867; hy: -0.06040960680288257;
lv: 461.7607703551739;}
step 6
time : 2004-01-01T23:35:00.000
equinoctial parameters: {a: 2.4396159E7;
ex: 0.11393312156755062; ey: 0.719345418868777;
hx: -0.009567941763699867; hy: -0.06040960680288257;
lv: 466.6919614270639;}
step 7
time : 2004-01-01T23:36:00.000
equinoctial parameters: {a: 2.4396159E7;
ex: 0.11393312156755062; ey: 0.719345418868777;
hx: -0.009567941763699867; hy: -0.06040960680288257;
lv: 471.4666804605399;}
step 8
time : 2004-01-01T23:37:00.000
equinoctial parameters: {a: 2.4396159E7;
ex: 0.11393312156755062; ey: 0.719345418868777;
hx: -0.009567941763699867; hy: -0.06040960680288257;
lv: 476.0676390794578;}
step 9
time : 2004-01-01T23:38:00.000
equinoctial parameters: {a: 2.4396159E7;
ex: 0.11393312156755062; ey: 0.719345418868777;
hx: -0.009567941763699867; hy: -0.06040960680288257;
lv: 480.4828961502498;}
step 10
time : 2004-01-01T23:39:00.000
equinoctial parameters: {a: 2.4396159E7;
ex: 0.11393312156755062; ey: 0.719345418868777;
hx: -0.009567941763699867; hy: -0.06040960680288257;
lv: 484.70541689946305;}
step 11
time : 2004-01-01T23:40:00.000
equinoctial parameters: {a: 2.4396159E7;
ex: 0.11393312156755062; ey: 0.719345418868777;
hx: -0.009567941763699867; hy: -0.06040960680288257;
lv: 488.7324362945905;}

Source code

The complete code for this example can be found in the source tree of the library, in file @src/tutorials/fr/cs/examples/propagation/SlaveMode.java@.

h3. Master mode

In this mode, the propagator, as a master, calls some callback methods provided by the application throughout its internal integration loop.

As in the slave mode tutorial above, let's define some initial state with:

• an inertial frame
• a date in some time scale
• a central attraction coefficient
• an orbit defined by its keplerian parameters.

// 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 = Math.toRadians(7); // inclination
double omega = Math.toRadians(180); // perigee argument
double raan = Math.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 @commons-math@ library, but it doesn't matter which integrator is used.

// with a minimum step of 0.001 and a maximum step of 1000
double minStep = 0.001;
double maxstep = 1000.0;
double[] absoluteTolerance = {
0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001
};
double[] relativeTolerance = {
1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7
};
new DormandPrince853Integrator(minStep, maxstep,
absoluteTolerance, relativeTolerance);
NumericalPropagator propagator = new NumericalPropagator(integrator);

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

Frame ITRF2005 = FramesFactory.getITRF2005(); // terrestrial frame
double ae = 6378137.; // equatorial radius in meter
double c20 = -1.08262631303e-3; // J2 potential coefficient

// potential coefficients arrays (only J2 is considered here)
double[][] c = new double[3][1];
c[0][0] = 0.0;
c[2][0] = c20;
double[][] s = new double[3][1];
ForceModel cunningham = new CunninghamAttractionModel(ITRF2005, ae, mu, c, s);

This force model is simply added to the propagator:

The propagator operating mode is set to master mode with fixed step and a 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.setMasterMode(60., new TutorialStepHandler());

Then, the initial state is set for the propagator:

propagator.setInitialState(initialState);

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 derived some class from the interface OrekitFixedStepHandler to realize a handleStep method, as follows:

private static class TutorialStepHandler implements OrekitFixedStepHandler {

```public void handleStep(SpacecraftState currentState, boolean isLast) {
System.out.println(" time : " + currentState.getDate());
System.out.println(" " + currentState.getOrbit());
if (isLast) {
System.out.println(" this was the last step ");
}
}
```

}

The same result, with the slave mode, would have required much more programming.

The printed results are shown below:

time : 2004-01-01T23:30:00.000
equinoctial parameters: {a: 2.4396159E7;
ex: 0.11393312156755062; ey: 0.719345418868777;
hx: -0.009567941763699867; hy: -0.06040960680288257;
lv: 441.0;}
time : 2004-01-01T23:31:00.000
equinoctial parameters: {a: 2.43956442447147E7;
ex: 0.11379555463005957; ey: 0.7193603780367466;
hx: -0.009567899603538911; hy: -0.06040916627598484;
lv: 446.281396252324;}
time : 2004-01-01T23:32:00.000
equinoctial parameters: {a: 2.4394130981128443E7;
ex: 0.1136592209017635; ey: 0.7193618943111063;
hx: -0.009567880963307815; hy: -0.06040784657447202;
lv: 451.5253608293239;}
time : 2004-01-01T23:33:00.000
equinoctial parameters: {a: 2.43917106643025E7;
ex: 0.11352728523102393; ey: 0.7193506733385568;
hx: -0.009568040191295672; hy: -0.0604057162729131;
lv: 456.6962051211052;}
time : 2004-01-01T23:34:00.000
equinoctial parameters: {a: 2.438852139218653E7;
ex: 0.11340245896040695; ey: 0.719328108923199;
hx: -0.009568504032782268; hy: -0.060402895325632394;
lv: 461.76152443974274;}
time : 2004-01-01T23:35:00.000
equinoctial parameters: {a: 2.4384730095000014E7;
ex: 0.11328681468263016; ey: 0.7192960781495514;
hx: -0.009569360025726937; hy: -0.060399538106672035;
lv: 466.69337761753326;}
time : 2004-01-01T23:36:00.000
equinoctial parameters: {a: 2.4380513630602702E7;
ex: 0.1131817075460741; ey: 0.719256704937947;
hx: -0.009570652226141946; hy: -0.06039581413595108;
lv: 471.4690168499321;}
time : 2004-01-01T23:37:00.000
equinoctial parameters: {a: 2.4376042191173874E7;
ex: 0.11308779626346281; ey: 0.7192121378362317;
hx: -0.009572383735759913; hy: -0.06039189021718959;
lv: 476.0711590450501;}
time : 2004-01-01T23:38:00.000
equinoctial parameters: {a: 2.437146715705459E7;
ex: 0.11300513750428136; ey: 0.7191643744054912;
hx: -0.0095745241076911; hy: -0.06038791660587941;
lv: 480.48785948784604;}
time : 2004-01-01T23:39:00.000
equinoctial parameters: {a: 2.436691418568236E7;
ex: 0.11293331934941137; ey: 0.7191151481891875;
hx: -0.009577019067983016; hy: -0.06038401850862352;
lv: 484.71206392270483;}
time : 2004-01-01T23:40:00.000
equinoctial parameters: {a: 2.4362480759924155E7;
ex: 0.11287160608332707; ey: 0.7190658726134559;
hx: -0.009579800635231736; hy: -0.060380292373369754;
lv: 488.7409705510215;}
this was the last step
Final date : 2004-01-01T23:40:30.000
Final state : equinoctial parameters: {a: 2.4360331842829864E7;
ex: 0.11284425291236964;
ey: 0.7190415678613115;
hx: -0.009581276101643383;
hy: -0.06037851608173257;
lv: 490.68230357067466;}

Source code

The complete code for this example can be found in the source tree of the library, in file @src/tutorials/fr/cs/examples/propagation/MasterMode.java@.

h3. Ephemeris mode

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 tutorials above, let's first define some initial state with:

• an inertial frame
• a date in some time scale
• a central attraction coefficient
• an orbit defined by its keplerian parameters.

// 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 = Math.toRadians(7); // inclination
double omega = Math.toRadians(180); // perigee argument
double raan = Math.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 @commons-math@ library.

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

The initial state is set for this propagator:

propagator.setInitialState(initialState);

Then, the propagator operating mode is simply set to ephemeris mode:

propagator.setEphemerisMode();

And the propagation is performed for a given duration.

SpacecraftState finalState =
propagator.propagate(new AbsoluteDate(initialDate, 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 = propagator.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.0092657655284;}
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

Source code

The complete code for this example can be found in the source tree of the library, in file @src/tutorials/fr/cs/examples/propagation/EphemerisMode.java@.

h2. 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 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.

double ae = 6378137.0; // equatorial radius in meter
double f = 1.0 / 298.257223563; // flattening
Frame ITRF2005 = FramesFactory.getITRF2005(); // terrestrial frame at an arbitrary date
BodyShape earth = new OneAxisEllipsoid(ae, f, ITRF2005);
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 VisibilityDetector which is derived from an ElevationDetector with the same specific parameters.

double maxcheck = 1.;
EventDetector sta1Visi = new VisibilityDetector(maxcheck, elevation, sta1Frame);

This EventDetector is added to the propagator:

Finally, the propagator is simply asked to perform until some final date, in slave mode by default.

It will propagate from the initial date to the first raising or for the fixed duration according to the behavior implemented in the VisibilityDetector class.

SpacecraftState finalState =
kepler.propagate(new AbsoluteDate(initialDate, 1500.));

The main application code is very simple, all the work is done inside the propagator thanks to the VisibilityDetector class especially created for the purpose.

This class extends the ElevationDetector class, furnished by OREKIT, by overriding the eventOccurred method with the special ability to print the results of the visibility check, both the raising and the setting time, and to stop the propagation just after the setting detection.

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

Visibility on station1 begins at 2004-01-01T23:31:52.097
Visibility on station1 ends at 2004-01-01T23:42:48.850
Final state : 2004-01-01T23:42:48.850

Source code

The complete code for this example can be found in the source tree of the library, in file @src/tutorials/fr/cs/examples/propagation/VisibilityCheck.java@.