The goal of this tutorial is to introduce orbital integration using SingleBodyAttraction class.
This class can replace all the different kinds of point mass interactions (ThirdBodyAttraction, NewtonianAttraction).
Using SingleBodyAttraction and InertialForces enables a richer modelling, allowing the user to compute the motion of a satellite in a reference frame that is not necessarily centered on the main attractor and does not necessarily possess inertial axis.
We will propagate the trajectory of a satellite departing from the Lagrange point L2 of the Earth-Moon system.
The equations of motion will be computed in a reference Frame centered on L2, its X-axis continuously oriented toward Earth and Moon.
The rotation of this frame means that we cannot simply use the fundamental principle of dynamics (Newton's second law).
Since we will consider only the gravitational attractions of the Earth and the Moon, the inertial reference frame most closely
related to our problem is a frame that is centered on the Earth-Moon barycenter, with inertial axis.
Let’s initialize our program. First, the time settings : the initial moment of integration, the length of integration (here in seconds), and the time interval between each output.
final AbsoluteDate initialDate = new AbsoluteDate(2000, 01, 01, 0, 0, 00.000, TimeScalesFactory.getUTC());
double integrationTime = 600000.;
double outputStep = 600.0;
Then the initial position and velocity of the satellite, located exactly at L2 point. We still haven't introduced a reference Frame, so the user must keep in mind in which reference Frame he wants to define the initial conditions of its satellite.
final PVCoordinates initialConditions = new PVCoordinates(new Vector3D(0.0, 0.0, 0.0), new Vector3D(0.0, 0.0, 0.0));
We need to load a few Celestial bodies, those we consider in gravitational interaction with our satellite, as well as the Earth-Moon barycenter, since we will compute our inertia forces with respect to its attached frame. We then create the reference frame attached to the L2 point, and the inertial reference frame attached to the Earth-Moon barycenter.
final CelestialBody earth = CelestialBodyFactory.getEarth();
final CelestialBody moon = CelestialBodyFactory.getMoon();
final CelestialBody earthMoonBary = CelestialBodyFactory.getEarthMoonBarycenter();
final Frame l2Frame = new L2Frame(earth, moon);
final Frame earthMoonBaryFrame = earthMoonBary.getInertiallyOrientedFrame();
final Frame inertiaFrame = earthMoonBaryFrame;
final Frame integrationFrame = l2Frame;
final Frame outputFrame = l2Frame;
We now transform our PVCoordinates into AbsolutePVCoordinates, define
the satellite attitude, and use all of it to build the SpacecraftState
final AbsolutePVCoordinates initialAbsPV = new AbsolutePVCoordinates(integrationFrame, initialDate, initialConditions);
Attitude arbitraryAttitude = new Attitude(integrationFrame,
new TimeStampedAngularCoordinates(initialDate,
new PVCoordinates(Vector3D.PLUS_I, Vector3D.PLUS_J),
new PVCoordinates(Vector3D.PLUS_I, Vector3D.PLUS_J)));
final SpacecraftState initialState = new SpacecraftState(initialAbsPV, arbitraryAttitude);
We will use a variable-step 8(5,3) Dormand-Prince integrator. This integrator needs a few initialization parameters. First, let's set boundaries on the integration steps, to prevent a too long computing or a too long integration step.
final double minStep = 0.001;
final double maxstep = 3600.0;
We also need to set the acceptable error of our integration. This error is used to adjust the integration step and does not stand for the overall error of the propagation.
final double positionTolerance = 0.001;
final double velocityTolerance = 0.00001;
final double massTolerance = 1.0e-6;
final double[] vecAbsoluteTolerances = {
positionTolerance, positionTolerance, positionTolerance,
velocityTolerance, velocityTolerance, velocityTolerance,
massTolerance
};
final double[] vecRelativeTolerances = new double[vecAbsoluteTolerances.length];
We can now define the numerical integrator of our propagator.
AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxstep,
vecAbsoluteTolerances,
vecRelativeTolerances);
And finally, we can build our propagator. First we use the NumericalPropagator
constructor with the previously defined numerical integrator, and then we specify the
properties of the evolution model. The use of SingleBodyAttraction means that
we do not constrain the type of orbit, and that we need to use setOrbitType(null)
as well as setIgnoreCentralAttraction(true).
The non-inertial aspect of our integration reference frame is handled by InertialForces, which will compute the
inertial accelerations that appear in this frame. Each attracting body will be provided
to the propagator by calling SingleBodyAttraction for each body in gravitational
interaction with our spacecraft.
We then only need to set the inital spacecraft state,
and the propagator mode, more details about propagator modes can be found in the Propagation tutorial .
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(null);
propagator.setIgnoreCentralAttraction(true);
propagator.addForceModel(new InertialForces(earthMoonBaryFrame));
propagator.addForceModel(new SingleBodyAbsoluteAttraction(earth));
propagator.addForceModel(new SingleBodyAbsoluteAttraction(moon));
propagator.setInitialState(initialState);
propagator.getMultiplexer().add(outputStep, new TutorialStepHandler("test.dat", outputFrame));
In the end we can do the propagation itself, from an initial date, for a given duration.
SpacecraftState finalState = propagator.propagate(initialDate.shiftedBy(integrationTime));
final PVCoordinates pv = finalState.getPVCoordinates(outputFrame);
System.out.println("initial conditions: " + initialConditions);
System.out.println("final conditions: " + pv);
Now, we can deal with our stephandler, it will print the position, velocities and acceleration at each output step of the integration.
private static class TutorialStepHandler implements OrekitFixedStepHandler {
private Frame outputFrame;
private TutorialStepHandler( final Frame frame) {
outputFrame = frame;
}
public void init(final SpacecraftState s0, final AbsoluteDate t) {
System.out.format(Locale.US,
"%s %s %s %s %s %s %s %s %s %s %n",
"date", " X", " Y",
" Z", " Vx", " Vy",
" Vz", " ax", " ay",
" az");
}
public void handleStep(SpacecraftState currentState) {
final AbsoluteDate d = currentState.getDate();
final PVCoordinates pv = currentState.getPVCoordinates(outputFrame);
System.out.format(Locale.US,
"%s %18.12f %18.12f %18.12f %18.12f %18.12f %18.12f %18.12f %18.12f %18.12f%n",
d, pv.getPosition().getX(),
pv.getPosition().getY(), pv.getPosition().getZ(),
pv.getVelocity().getX(), pv.getVelocity().getY(),
pv.getVelocity().getZ(), pv.getAcceleration().getX(),
pv.getAcceleration().getY(),
pv.getAcceleration().getZ());
}
public void finish(final SpacecraftState finalState) {
final PVCoordinates finalPv =
finalState.getPVCoordinates(outputFrame);
System.out.println();
System.out.format(Locale.US,
"%s %12.0f %12.0f %12.0f %12.0f %12.0f %12.0f%n",
finalState.getDate(), finalPv.getPosition().getX(),
finalPv.getPosition().getY(),
finalPv.getPosition().getZ(),
finalPv.getVelocity().getX(),
finalPv.getVelocity().getY(),
finalPv.getVelocity().getZ());
System.out.println();
}
}
The printed results are shown below:
date X Y Z Vx Vy Vz ax ay az
2000-01-01T00:00:00.000Z 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000007203264 -0.000023617334 0.000000988884
2000-01-01T00:10:00.000Z 1.310181730693 -4.509937975829 0.177995372228 0.004389441334 -0.015464490879 0.000593311180 0.000007424953 -0.000027930976 0.000000988817
2000-01-01T00:20:00.000Z 5.292760978714 -19.075027472324 0.711964496463 0.008906039665 -0.033517172720 0.001186577761 0.000007627122 -0.000032244627 0.000000988736
2000-01-01T00:30:00.000Z 12.020527119451 -45.248282259272 1.601878011878 0.013538105089 -0.054158285678 0.001779791113 0.000007690269 -0.000035156300 0.000000989424
2000-01-01T00:40:00.000Z 21.178259937228 -80.021547310428 2.850254956138 0.017011615221 -0.062141479011 0.002381481685 0.000005918995 -0.000015387558 0.000001002821
2000-01-01T00:50:00.000Z 32.466255913066 -120.334453824629 4.459646951922 0.020640563371 -0.072665310900 0.002983149974 0.000006174249 -0.000019691905 0.000001002738
2000-01-01T01:00:00.000Z 45.976396611751 -167.736449656047 6.430024082490 0.024416814768 -0.085771787993 0.003584764051 0.000006410007 -0.000023996366 0.000001002640
2000-01-01T01:10:00.000Z 61.793555128366 -223.777138815637 8.761351268544 0.028328673735 -0.101460960959 0.004186315403 0.000006626275 -0.000028300882 0.000001002529
2000-01-01T01:20:00.000Z 79.995588058317 -290.006145592032 11.453588310832 0.032364447630 -0.119732846726 0.004787795477 0.000006823058 -0.000032605399 0.000001002403
2000-01-01T01:30:00.000Z 100.653337431216 -367.973095343246 14.506689867546 0.036512309222 -0.140585791722 0.005389196595 0.000006837664 -0.000034976714 0.000001003336
2000-01-01T01:40:00.000Z 123.453665764425 -454.667651202580 17.923140739980 0.039513023305 -0.148793728322 0.005998979319 0.000005129687 -0.000015774366 0.000001016276
2000-01-01T01:50:00.000Z 148.100098098847 -547.040990119446 21.705450574838 0.042666722324 -0.159546931615 0.006608707351 0.000005379399 -0.000020069667 0.000001016148
2000-01-01T02:00:00.000Z 174.682530302316 -646.639412217862 25.853573367492 0.045964410184 -0.172877349174 0.007218354413 0.000005609651 -0.000024365071 0.000001016006
2000-01-01T02:10:00.000Z 203.283853139330 -755.009261136276 30.367457953699 0.049394412241 -0.188785025652 0.007827911891 0.000005820447 -0.000028660521 0.000001015850
2000-01-01T02:20:00.000Z 233.979953476934 -873.696897209172 35.247047991524 0.052945056828 -0.207269972054 0.008437371140 0.000006011793 -0.000032955961 0.000001015679
2000-01-01T02:30:00.000Z 266.839718231943 -1004.248702962090 40.492281925576 0.056604717930 -0.228332679264 0.009046723193 0.000006192961 -0.000037365335 0.000001015430
2000-01-01T02:40:00.000Z 301.553619978550 -1143.666361515510 46.105600609180 0.059130769542 -0.236774827689 0.009664348966 0.000004333782 -0.000016135290 0.000001029339
2000-01-01T02:50:00.000Z 337.827105187818 -1288.892787478942 52.089480854137 0.061805266666 -0.247741890570 0.010281901244 0.000004577970 -0.000020421606 0.000001029166
2000-01-01T03:00:00.000Z 375.748076996488 -1441.470994477480 58.443860505298 0.064620448977 -0.261280773723 0.010899345538 0.000004802733 -0.000024708015 0.000001028979
2000-01-01T03:10:00.000Z 415.397450190460 -1602.944087970172 65.168672162654 0.067564662535 -0.277391515874 0.011516673144 0.000005008075 -0.000028994460 0.000001028777
2000-01-01T03:20:00.000Z 456.849148063114 -1774.855186584263 72.263843192313 0.070626256317 -0.296074122196 0.012133875320 0.000005194002 -0.000033280886 0.000001028561
2000-01-01T03:20:00.000Z 457 -1775 72 0 -0 0
initial conditions: {P(0.0, 0.0, 0.0), V(0.0, 0.0, 0.0), A(0.0, 0.0, 0.0)}
final conditions : {2000-01-01T03:20:00.000, P(456.8491480631136, -1774.8551865842628, 72.26384319231306), V(0.07062625631702495, -0.2960741221959937, 0.012133875319548914), A(5.194001612691349E-6, -3.3280886249059965E-5, 1.0285608966641585E-6)}
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/PropagationInNonInertialFrame.java.