STK Components for .NET 2020 r2

## 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 |

Configuring Context

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

C#

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 JplDE421 ephem = new JplDE421(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.

C#

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.

C#

SphericalHarmonicGravity gravity = new SphericalHarmonicGravity(); gravity.TargetPoint = position.IntegrationPoint; int degree = 41; int order = 41; SphericalHarmonicGravityModel gravityModel = SphericalHarmonicGravityModel.ReadFrom(gravityEGM96FilePath); gravity.GravityField = new SphericalHarmonicGravityField(gravityModel, degree, order, true, SphericalHarmonicsTideType.None); 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.

C#

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 STK Components. 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.

C#

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:

C#

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

C#

// 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