Code Sample

The following code samples illustrate the various pieces required to set up, configure, and numerically propagate a state.

Note
The functionality described in this topic requires a license for the Orbit Propagation Library.

Configuring Context

First, we configure the context in DME Component Libraries to include EOP data, using EarthOrientationParameters, and loading JPL ephemerides using JplDE440. The EOP and JPL data are not required to run, but they will affect the accuracy of the force models during propagation.

Copy

EarthCentralBody earth = CentralBodiesFacet.GetFromContext().Earth;
SunCentralBody sun = CentralBodiesFacet.GetFromContext().Sun;
MoonCentralBody moon = CentralBodiesFacet.GetFromContext().Moon;

// Load EOP table into the library
earth.OrientationParameters = EarthOrientationParametersFile.ReadData(eopFilePath);

// Load JPL Ephemerides into the library
JplDE440 ephem = new JplDE440(deFilePath);
// The following updates the "CenterOfMassPoint" on each of the CentralBodies
ephem.UseForCentralBodyPositions(CentralBodiesFacet.GetFromContext());
// The following is needed to accurately model the precise orientation of the Moon
moon.FixedFrame = ephem.GetMoonTopographicFixedFrame();

Defining State Elements

Next, we decide which elements need to be integrated in the state and which should be evaluated as auxiliary data during propagation. A position which depends on complex force models needs to be integrated to predict its values. In contrast, a density is simply a function of position and has no differential equation associated with it. The differential equations will be filled in next.

Copy

NumericalPropagatorDefinition state = new NumericalPropagatorDefinition();

string positionID = "MyVehiclePosition";
PropagationNewtonianPoint position = new PropagationNewtonianPoint(positionID, earth.InertialFrame,
                                                                   initialPosition, initialVelocity);

state.Epoch = epoch;
state.IntegrationElements.Add(position);
// state.IntegrationElements.Add(angleOfAttack);
// etc.

// Add auxiliary output
AuxiliaryStateScalar densityOutput = new AuxiliaryStateScalar();
string densityID = "MyAtmosphericDensity";
densityOutput.Identification = densityID;
state.AuxiliaryElements.Add(densityOutput);

Configuring Force Models

In general, each PropagationStateElement will have some way of defining its derivative for the integrator. In the case of a PropagationVector, the derivative is simply a Vector whose value represents the derivative of the highest order of the represented vector state. However, the most common element will likely be the PropagationNewtonianPoint which models the position and velocity of a body subject to a list of ForceModels.

Copy

SphericalHarmonicGravity gravity = new SphericalHarmonicGravity();
gravity.TargetPoint = position.IntegrationPoint;
int degree = 41;
int order = 41;
SphericalHarmonicGravityModel gravityModel = SphericalHarmonicGravityModel.ReadFrom(gravityEGM96FilePath);
SolidTideModel tideModel = null;
gravity.GravityField = new SphericalHarmonicGravityField(gravityModel, degree, order,
                                                         true, tideModel);

ThirdBodyGravity thirdBody = new ThirdBodyGravity();
thirdBody.CentralBody = earth;
const double sunGm = 1.327122e20; // m^3 / sec^2
const double moonGm = 4.902801076e12; // m^3 / sec^2
thirdBody.AddThirdBody(sun.Name, sun.CenterOfMassPoint, sunGm);
thirdBody.AddThirdBody(moon.Name, moon.CenterOfMassPoint, moonGm);
thirdBody.TargetPoint = position.IntegrationPoint;

SimpleSolarRadiationForce srp = new SimpleSolarRadiationForce();
ScalarOccultation occultationFactor = new ScalarOccultationDualCone();
occultationFactor.IlluminatingBody = sun;
occultationFactor.OccludingBodies.Add(earth);
occultationFactor.OccludingBodies.Add(moon);
occultationFactor.TargetPoint = position.IntegrationPoint;
occultationFactor.VectorType = RadiationVectorType.TruePosition;
srp.OccultationFactor = occultationFactor;
srp.CoefficientOfReflectivity = 1.0;
srp.ReferenceArea = 20.0;

AtmosphericDragForce drag = new AtmosphericDragForce();
ConstantSolarGeophysicalData geoData = new ConstantSolarGeophysicalData(150.0, 3.00);
ScalarDensityJacchiaRoberts density = new ScalarDensityJacchiaRoberts(position.IntegrationPoint, geoData);
drag.Density = density;
drag.CoefficientOfDrag = 2.2;
drag.ReferenceArea = 20.0;

position.AppliedForces.Add(gravity);
position.AppliedForces.Add(thirdBody);
position.AppliedForces.Add(srp);
position.AppliedForces.Add(drag);
position.Mass = mass;

densityOutput.AuxiliaryScalar = density;

Configuring Integrator

The last step before creating the propagator is to choose which integration scheme to use in order to propagate the state and to configure its settings.

Copy

RungeKuttaFehlberg78Integrator integrator = new RungeKuttaFehlberg78Integrator();
integrator.InitialStepSize = stepSize;
integrator.StepSizeBehavior = KindOfStepSize.Relative;
integrator.MaximumStepSize = double.MaxValue;
integrator.MinimumStepSize = 0.001;
// This is STK's default and truncates each step to 3 decimal places
integrator.StepTruncationOrder = -3;

state.Integrator = integrator;

Propagating

After creating a NumericalPropagator, we use it to generate ephemeris over a given span of time, then create a PointInterpolator to represent the propagated ephemeris for further analysis in DME Component Libraries. Note that this point should be used to define the LocationPoint of a platform, rather than the IntegrationPoint which only represents the position during the integration process.

Copy

int outputSparsity = 1; // 1 indicates to output after every integration step, 2 every other step, etc, etc
Duration propagationTime = Duration.FromSeconds(totalTime);
NumericalPropagationStateHistory rawData = propagator.Propagate(propagationTime, outputSparsity);

// Recover the position, velocity, (and acceleration if specified on the position before creating the propagator)
DateMotionCollection<Cartesian> ephemeris = rawData.GetDateMotionCollection<Cartesian>(positionID);

DateMotionCollection<double> densityData = rawData.GetDateMotionCollection<double>(densityID);

Point platformPosition = new PointInterpolator(position.IntegrationFrame, InterpolationAlgorithmType.Hermite, 6, ephemeris);

// Perform analysis using other libraries...

The above code snippet illustrates a simple propagation over time, however, it is also possible to have the duration of propagation determined not by an overall time but by some other stopping condition which is a function of time and the state. This can be accomplished by attaching an event handler to the StepTaken event, which will fire every time the propagator takes a step, as shown below:

Copy

// The following allows us to propagate all at once or step the state manually
NumericalPropagator propagator = state.CreatePropagator(group);

// Note: the "IntegrationPoint" from the PropagationNewtonianPoint cannot be evaluated
// outside of the context of the NumericalPropagator.  However, the StepTaken event
// occurs inside the scope of propagation and the point is created in the same EvaluatorGroup.
// So, the evaluation of the point will represent the current position and velocity in the state.
// Similarly, any evaluator which was built off of the point will evaluate using the current
// position and velocity.

PointEvaluator activePointEvaluator =
    GeometryTransformer.ObservePoint(position.IntegrationPoint, position.IntegrationFrame, group);

// Add an event handler to check each step during propagation.
propagator.StepTaken += (sender, e) =>
{
    JulianDate t = e.CurrentTime;
    double[] x = e.CurrentState;

    StepSizeInformation info = e.StepSizeInformation;

    // The following two methods of retrieving the position are equally valid...
    Motion<Cartesian> posVel = propagator.Converter.ConvertState<Cartesian>(positionID, x);
    // ...but this one is faster and eliminates the need to perform a transformation manually
    // if the position is needed in a different frame
    posVel = activePointEvaluator.Evaluate(e.CurrentTime, 1);

    // Check for a "stopping condition" and if reached stop propagation
    if (posVel.Value.Magnitude > 1e8)
    {
        e.Indication = PropagationEventIndication.StopPropagationAfterStep;
    }
    else
    {
        e.Indication = PropagationEventIndication.ContinuePropagation;
    }
};

It is also possible to use the propagator in situations where surrounding code is iterating (or "shooting") to find some desired final state based on an initial guess. It's always possible to create a new propagator, but if the only things that change are the initial epoch and initial state, it's also possible to reset the propagator to a new set of initial conditions and try again. To do this, use the PropagationStateConverter in reverse to specify new initial conditions and update the initial state accordingly, then call the Reset method on the propagator.

Without calling Reset, each call to Propagate will start where the previous call left off. This is designed to allow systems which propagate for extended periods of time to stop intermittently to verify the ephemeris, write or commit the data, or provide a means of specifying segments of propagation. However, note that if any of the differential equations defined on the state elements need to change, a new propagator will need to be created based on the new equations. However, to start the new propagator where the last one left off, simply set the new propagator's InitialState to a copy of the old propagator's CurrentState.

Copy

// If conducting an iterative analysis with the same configuration (same state, force models, correctors, etc)...
// Reinitialize the state and repropagate

double[] newInitialState = propagator.InitialState; // Gets a copy

Motion<Cartesian> newInitialConditions = new Motion<Cartesian>(nextInitialPosition, nextInitialVelocity);
propagator.Converter.AdjustState(positionID, newInitialState, newInitialConditions);

propagator.Reset(epoch, newInitialState);
rawData = propagator.Propagate(Duration.FromSeconds(totalTime), 1);

// To use a different state or different force models, create a new NumericalPropagatorDefinition with new StateElements and ForceModels