Click or drag to resize

Force Models

A ForceModel defines a particular force acting on a body based on a model of the physics in the body's environment. When combined together with a PropagationNewtonianPoint, the force models define a body's equations of motion according to Newton's second law. The Newtonian position uses the force models along with a specified mass to determine the acceleration used to integrate the position and velocity.

Note Note

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

Gravity Models

There are several different gravity models available in DME Component Libraries:

  • Spherical Harmonic Gravity

    Spherical harmonic gravity models the gravitational effects produced by the mass distribution of a large body. The gravity is handled by three types: SphericalHarmonicGravityModel, SphericalHarmonicGravityField, and SphericalHarmonicGravity.

    The gravity model represents the full set of coefficients and can be read in from a ".grv" file. This is then used to construct the immutable gravity field by downselecting to the desired fidelity in the degree and order of the represented field, as well as configuring other options, such as the inclusion of tidal data. This gravity field is then used to define the gravity force at a given position.

    C#
    Copy
    // The partial configuration of a PropagationNewtonianPoint was demonstrated in an earlier recipe.
PropagationNewtonianPoint satellite = SeeRepresentASatellitePosition();

// Configuring a SphericalHarmonicGravity ForceModel requires the Point of the object whose force is being calculated,
// as well as the SphericalHarmonicGravityField with the subset of desired coefficients describing the gravity field,
// which are pulled from the full SphericalHarmonicGravityModel, which is generally read in from a file.
SphericalHarmonicGravityModel model = SphericalHarmonicGravityModel.ReadFrom(pathToGravFile);

// For performance reasons the SphericalHarmonicGravity ForceModel does not operate directly
// on the gravity model. Instead an immutable object with the subset of the coefficients
// at the desired degree and order is created.
SphericalHarmonicGravityField field = new SphericalHarmonicGravityField(model, 21, 21, true, new PermanentSolidTideModel());
// The degree and order of the gravity field determine the accuracy of the force. Higher
// degrees and orders will lead to a more accurate result, but are more computationally intensive.
// At the bare minimum the degree should be 2, and can go as high as the maximum degree of the model
// (which is in the .grv file). The order can be equal or less than the degree.
// The SphericalHarmonicsGravity force can include or exclude the two body gravitational force. If it is excluded
// then a separate TwoBodyGravity force should be added to the PropagationNewtonianPoint.
// The last parameter adds the Earth's permanent tides to the gravity field.

// Either of the above fields are fully configured and can be used to create a SphericalHarmonicsGravity object
// (or set to the GravityField property of an existing SphericalHarmonicsGravity object).
SphericalHarmonicGravity harmonicGrav = new SphericalHarmonicGravity(satellite.IntegrationPoint, field);

// The gravity force is now fully configured and can be added as an AppliedForce.

    In addition, there is the simple TwoBodyGravity force, which treats the central gravitational body as a point mass, and the ThirdBodyGravity perturbation force which models the addition of the point-mass gravity from bodies other than the primary gravitational body.

  • Two Body Gravity

    C#
    Copy
    // The partial configuration of a PropagationNewtonianPoint was demonstrated in an earlier recipe.
PropagationNewtonianPoint satellite = SeeRepresentASatellitePosition();

CentralBody earth = CentralBodiesFacet.GetFromContext().Earth;
CentralBody moon = CentralBodiesFacet.GetFromContext().Moon;

// A TwoBodyGravityForce requires three things, the point of the object the force is applied to, the
// CentralBody of the gravitational body, and a double of the gravitational parameter
// of the body, or mu. As always, the value of the gravitational parameter must be in SI Notation.
// In this case that means m^3 / s^2.
TwoBodyGravity grav = new TwoBodyGravity(satellite.IntegrationPoint, earth, WorldGeodeticSystem1984.GravitationalParameter);

// The majority of the time a TwoBodyGravity force is created, it will represent a rotation around earth, so the
// default constructors initialize the force with the above values. The following processes create equivalent forces:
TwoBodyGravity gravIdentical = new TwoBodyGravity(satellite.IntegrationPoint);
TwoBodyGravity gravIdentical2 = new TwoBodyGravity();
gravIdentical2.TargetPoint = satellite.IntegrationPoint;

// However, any body can serve as the focal point for a TwoBodyGravity ForceModel
TwoBodyGravity tbGravMoon = new TwoBodyGravity(satellite.IntegrationPoint, moon, 4.902795e12);
// Note that the above force model is for a satellite orbiting the moon. If you desire to model the perturbing effects
// of the moon's (or other body's) gravity on an earth-orbiting satellite use "ThirdBodyGravity" which accounts for the
// effect of the moon's gravity on the Earth.

// All of the above example forces are now completely configured and ready to be added to a PropagationNewtonianPoint.

  • Third Body Gravity

    C#
    Copy
    // Getting the necessary bodies from the CentralBodiesFacet.
CentralBody moon = CentralBodiesFacet.GetFromContext().Moon;
CentralBody sun = CentralBodiesFacet.GetFromContext().Sun;

// The partial configuration of a PropagationNewtonianPoint was demonstrated in an earlier recipe.
PropagationNewtonianPoint satellite = SeeRepresentASatellitePosition();

// Construct the ThirdBodyGravity force model with the Point of the target object.
ThirdBodyGravity thirdBodyGrav = new ThirdBodyGravity(satellite.IntegrationPoint);

// The CentralBody of the force is automatically set to the Earth, however the desired
// "third bodies" must still be added.
ThirdBodyGravity.ThirdBodyInfo moonInfo = new ThirdBodyGravity.ThirdBodyInfo("moon", moon.CenterOfMassPoint, 4.902795e12);
ThirdBodyGravity.ThirdBodyInfo sunInfo = new ThirdBodyGravity.ThirdBodyInfo("sun", sun.CenterOfMassPoint, 1.327122e20);
thirdBodyGrav.ThirdBodies.Add(moonInfo);
thirdBodyGrav.ThirdBodies.Add(sunInfo);

// The ThirdBodyGravity force is now completely configured and ready to be added to the PropagationNewtonianPoint.
Simple Atmospheric Drag

The AtmosphericDragForce is constructed using one of several ScalarAtmosphericDensity models, as well as Scalars representing the object's CoefficientOfDrag and ReferenceArea. Since the position and velocity are required to compute density and in order to avoid possible inconsistency, the drag force acquires the position and velocity from the scalar density model in order to compute the force.

It's important to note that there are many different models for the Earth's atmosphere and some are more accurate than others. Ultimately the variability and uncertainty in modeling the atmosphere makes atmospheric drag the least physically accurate force model in the provided set. However, for low-altitude satellites having an approximation for the atmospheric drag is crucial for determining an estimation for its trajectory.

C#
Copy
// The partial configuration of a PropagationNewtonianPoint was demonstrated in an earlier recipe.
PropagationNewtonianPoint satellite = SeeRepresentASatellitePosition();

// First we will create an atmospheric drag force from the default constructor.
AtmosphericDragForce drag = new AtmosphericDragForce();

// An AtmophericDrag force needs three properties set before it can be used in an evaluation:
// A ScalarAtmosphericDensity (which also contains the object's position data), a
// Scalar Coefficient of Drag, and a Scalar Reference Area. Note that doubles
// are implicitly converted into ScalarFixed objects.
drag.CoefficientOfDrag = 2.2;
drag.ReferenceArea = 20.0;

// Lastly we need to set the density property, there are currently five Density models that come with DME Component Libraries,
// and more can be added by the user. These five models are Jacchia 1970, Jacchia-Roberts, NRLMSISE 2000, MSISE 1990,
// and MSIS 1986. Each of these models takes into account the solar and magnetic flux represented by a
// SolarGeophysicalData object. We will use the a default instance of ConstantSolarGeophysicalData for this purpose.
SolarGeophysicalData spaceWeather = new ConstantSolarGeophysicalData();

ScalarAtmosphericDensity densityJ1970 = new ScalarDensityJacchia1970(satellite.IntegrationPoint, spaceWeather);
ScalarAtmosphericDensity densityJR = new ScalarDensityJacchiaRoberts(satellite.IntegrationPoint, spaceWeather);
ScalarAtmosphericDensity densityMSIS2000 = new ScalarDensityMsis2000(satellite.IntegrationPoint, spaceWeather);
ScalarAtmosphericDensity densityMSIS90 = new ScalarDensityMsis90(satellite.IntegrationPoint, spaceWeather);
ScalarAtmosphericDensity densityMSIS86 = new ScalarDensityMsis86(satellite.IntegrationPoint, spaceWeather);

// All of the five previous density models are configured and ready to add to an AtmosphericDragForce.
drag.Density = densityJ1970;

// The drag force is completely configured and can be added as an AppliedForce.
Simple Solar Radiation Force

  • Solar Radiation Force

    SimpleSolarRadiationForce models the effect of pressure exerted by radiation from the Sun. The force is modeled as being projected along the vector from the Sun to the object. The entire ReferenceArea is taken to be orthogonal to the incident radiation.

    SimpleSolarRadiationForce is constructed using a ScalarOccultation model, as well as Scalars representing the object's coefficient of reflectivity, normalized area, and the Sun's luminosity. Similar to atmospheric drag force, the position is used by the occultation factor and also provides the position for the force.

    By default, a constant value (ConstantsSolarLuminosity2015IauResolutionB3) is used for luminosity, but a time-varying Scalar can be used instead.

    Three basic occultation factor models are available: ScalarOccultationDualCone, ScalarOccultationCylindrical, and ScalarOccultationNoShadow. These models take an observing Point, and CentralBody instances representing the illuminating body and occluding bodies. The VectorType property can be used to specify whether the calculation uses the true position for all bodies involved, or the apparent position that adjusts for light time delay and aberration. Also note that while an occultation factor model is needed to compute the SimpleSolarRadiationForce, it can also be used independently as a metric to determine the percentage of illumination at a particular (potentially time-varying) location in space.

    C#
    Copy
    // Getting the necessary bodies from the CentralBodiesFacet.
CentralBody earth = CentralBodiesFacet.GetFromContext().Earth;
CentralBody moon = CentralBodiesFacet.GetFromContext().Moon;
CentralBody sun = CentralBodiesFacet.GetFromContext().Sun;

// The partial configuration of a PropagationNewtonianPoint was demonstrated in an earlier recipe.
PropagationNewtonianPoint satellite = SeeRepresentASatellitePosition();

// There are two main types of Occultation model, Dual-Cone and Cylindrical. Of the two Dual-Cone
// is far more accurate, and should nearly always be used. Cylindrical has mostly been provided as
// a way to verify results against older models. For convenience there is also a No-Shadow model, and
// a No-Shadow condition can also be attained with the other models by not adding any occulting bodies.
// Lastly there is a Regulated Dual-Cone variety, which should only be used in conjunction with
// SolarRadiationBoundaryMitigation, which is covered in a later recipe.
ScalarOccultation occultationDualCone = new ScalarOccultationDualCone(sun, satellite.IntegrationPoint, earth, moon);
ScalarOccultation occultationCylindrical = new ScalarOccultationCylindrical(sun, satellite.IntegrationPoint, earth, moon);
// A No-Shadow occultation model does not require eclipsing bodies, although they can be added,
// and are ignored if added).
ScalarOccultation occultationNone = new ScalarOccultationNoShadow(sun, satellite.IntegrationPoint);

// All of the above ScalarOccultation models are fully configured and ready to be used within
// a SolarRadiationForce, however there is also an advanced property that can be set, the VectorType.
// By default all ScalarOccultation models' will correct for light time delay and velocity aberration
// when finding the vectors from the TargetPoint to the illuminating and occluding bodies for maximum accuracy.
// If a higher performance is desirable at the expense of this correction the VectorType can be sent to
// simply use the true positions of the bodies without corrections.
occultationCylindrical.VectorType = RadiationVectorType.TruePosition;

// The other two required pieces of information are Scalars defining the Coefficient of Reflectivity
// and the normalized area of the satellite. Note that doubles
// are implicitly converted into ScalarFixed objects.
SimpleSolarRadiationForce srf = new SimpleSolarRadiationForce(occultationDualCone, 1.0, 20.0);
// The force object now models the solar radiation force on the satellite assuming a Dual-Cone shadow model
// the Earth and Moon as eclipsing bodies. The modeled satellite properties are of a black-body
// made up of a 20 square meter surface.

// The above force is now ready to be added to the satellite's AppliedForces, but there are some additional
// optional properties.

// There are two possible formats for the Coefficient of Reflectivity, BlackBody and RadiationPressure.
// These are related by the following simple formula:
// RadiationPressure = BlackBody + 1
// The BlackBody format is more common among astronomers, physicists, or other professions interested mainly in the
// optical properties of a surface. 0.0 represents a perfectly absorbent surface, and 1.0 a perfectly reflective one.
// When calculating RadiationPressure a perfectly absorbent surface receives an impulse equal to that of the impacting
// light, so is in this format a perfectly absorbent body is given a value of 1.0. A perfectly reflective body
// receives twice this impulse as the photon is rebounded, so the coefficient is defined as 2.0.
// The type of coefficient reflectivity is set to RadiationPressure by default.
srf.CoefficientOfReflectivityType = CoefficientOfReflectivityType.BlackBody;
srf.CoefficientOfReflectivity = 0.0;
// This force will return an identical result as when it was first created with a reflectivity type of RadiationPressure
// and a coefficient value of 1.0.

// Lastly the Scalar of the Luminosity value of the illuminating body can be changed. By default it is a fixed scalar with
// the value of Constants.SolarLuminosity;

// The above radiation force is completely configured and ready to add as an AppliedForce on the satellite.

// The most common application of a solar force will be the situation we have manually configured above:
// using Dual-Cone occultation with the sun as the illuminating body and
// the Earth and Moon as eclipsing bodies. That scenario can easily be attained with the following constructor, which
// only requires a point, coefficient of reflectivity, and a reference area.
SimpleSolarRadiationForce srfIdentical = new SimpleSolarRadiationForce(satellite.IntegrationPoint, 1.0, 20.0);

  • Solar Radiation Boundary Mitigation

    A fourth occultation factor is also available: ScalarOccultationRegulatedDualCone. This model must be used when correcting the solar radiation force with SolarRadiationBoundaryMitigation during propagation. If the corrector is not being used, then the normal dual-cone shadow model should be used instead.

    While the dual-cone model will most frequently be used in conjunction with the solar radiation force, it can also be used in more general cases, such as the occlusion of a target body by its moon(s), or as a simple means of determining whether a vehicle is in direct, partial, or no sunlight.

    C#
    Copy
    EarthCentralBody earth = CentralBodiesFacet.GetFromContext().Earth;

// We begin by configuring our NumericalPropagatorDefinition as shown in "Create Numerical Propagator"
// With the exception of not adding the solar radiation force.
NumericalPropagatorDefinition definition = new NumericalPropagatorDefinition();

PropagationNewtonianPoint satellite = SeeRepresentASatellitePosition();

satellite.AppliedForces.Add(SeeCreateSphericalHarmonicsGravityForce());
satellite.AppliedForces.Add(SeeCreateThirdBodyGravityForce());
satellite.AppliedForces.Add(SeeCreateAtmosphericDragForce());

definition.Add(satellite);

definition.Epoch = new GregorianDate(2011, 10, 2).ToJulianDate();

definition.Integrator = SeeCreateNumericalIntegrator();

// Boundary mitigation requires a two step process. A SolarRadiationBoundaryMitigation object
// must be created and added as a PropagationStateCorrector to the definition.
// However boundary mitigation can only be performed on a solar radiation force made with a
// regulated dual cone occultation model.

// Start with the SimpleSolarRadiationForce we created in "Create Solar Radiation Force"
SimpleSolarRadiationForce srf = SeeCreateSolarRadiationForce();

// You can either construct a ScalarOccultationRegulatedDualCone from scratch, or from an
// existing ScalarOccultationDualCone.
srf.OccultationFactor = new ScalarOccultationRegulatedDualCone((ScalarOccultationDualCone)srf.OccultationFactor);

// The solar radiation force can now be added to the satellite and used to create the corrector.
satellite.AppliedForces.Add(srf);

// The corrector must be constructed with the satellite, the force to correct, the satellite's Mass, and
// since this is an Encke corrector it also needs enough information to create a TwoBodyGravity force.
SolarRadiationBoundaryMitigation boundaryMitigation = new SolarRadiationBoundaryMitigation(satellite, earth,
            WorldGeodeticSystem1984.GravitationalParameter, srf, satellite.Mass);

// Now we add the corrector to the NumericalPropagatorDefinition.
definition.StateCorrectors.Add(boundaryMitigation);

// The following propagator will now propagate a satellite with (in addition to the other forces)
// a solar radiation force that has been corrected for boundary mitigation.
definition.CreatePropagator();
N-Plate Models for Atmospheric Drag and Solar Radiation Pressure

  • N-Plate Atmospheric Drag

    In addition to the AtmosphericDragForce model that is low-fidelity because it assumes a spherical spacecraft, a medium-fidelity model for atmospheric drag called NPlateAtmosphericDragForce is also available.

    This model splits the spacecraft into individual plates that have distinct areas, pointing directions, and drag coefficients. These individual plates are modeled by DragPlates that are stored in a NPlateModel. There are several different concrete DragPlate types available:

    • A DragBodyPlate points in a specific direction defined by its UnitNormalCartesian, which is defined in the spacecraft's BodyAxes. Drag is only computed for a body plate if it is facing toward the incoming atmosphere. Drag is not reduced if one plate is supposed to be in front of another, but both are pointing toward the incoming atmosphere. This is a limitation on the fidelity of the model.

    • A DragTwoDofSolarPanelPlate points in the apparent direction of the Sun and its pointing is independent of the spacecraft's BodyAxes.

      If the solar panels are facing away from the incoming atmosphere, the shadowed sides of the solar panels will be modeled with the same area and drag coefficient as the Sun-facing sides, but the plate will be modeled as pointing directly away from the Sun.

    • A DragOneDofSolarPanelPlate points as close to the apparent direction of the Sun as it can while being constrained to only rotate around its UnitAxisCartesian, which is defined in the spacecraft's BodyAxes.

      If the solar panels are facing away from the incoming atmosphere, the shadowed sides of the solar panels will be modeled with the same area and drag coefficient as the Sun-facing sides, but the plate will be modeled as pointing in the opposing direction from its typical pointing.

    A code example of how to build a NPlateAtmosphericDragForce instance is given as follows:

    C#
    Copy
    // Getting the necessary bodies from the CentralBodiesFacet.
EarthCentralBody earth = CentralBodiesFacet.GetFromContext().Earth;
SunCentralBody sun = CentralBodiesFacet.GetFromContext().Sun;

// The partial configuration of a PropagationNewtonianPoint was demonstrated in an earlier recipe.
PropagationNewtonianPoint satellite = SeeRepresentASatellitePosition();

// An N-Plate model can be created from individual DragPlates such as DragBodyPlate, which
// model plates that point in specific directions in the spacecraft's body frame.
NPlateModel nPlateModel = new NPlateModel();
nPlateModel.Add(new DragBodyPlate("PlusZ", "Body", 5.0, UnitCartesian.UnitZ, 1.5));
nPlateModel.Add(new DragBodyPlate("PlusY", "Body", 5.0, UnitCartesian.UnitY, 1.5));
nPlateModel.Add(new DragBodyPlate("PlusX", "Body", 5.0, UnitCartesian.UnitX, 1.5));
nPlateModel.Add(new DragBodyPlate("MinusZ", "Body", 5.0, -UnitCartesian.UnitZ, 1.5));
nPlateModel.Add(new DragBodyPlate("MinusY", "Body", 5.0, -UnitCartesian.UnitY, 1.5));
nPlateModel.Add(new DragBodyPlate("MinusX", "Body", 5.0, -UnitCartesian.UnitX, 1.5));

// Another DragPlate is the DragTwoDofSolarPanelPlate, that models a solar panel that is able to always point
// directly toward the apparent Sun direction.
nPlateModel.Add(new DragTwoDofSolarPanelPlate("SolarPanel", "SolarPanels", 60.0, 2.5));

// The final DragPlate is the DragOneDofSolarPanelPlate, that models a solar panel that is only able to point toward the
// apparent Sun direction to the extent that its single rotation axis allows it to do so.
nPlateModel.Add(new DragOneDofSolarPanelPlate("AxisPanel", "SolarPanels", 5.0,
                                              UnitCartesian.UnitZ, 2.3));

// However, N-Plate models are usually read from .nplate files.
nPlateModel = NPlateModel.ReadFrom(pathToNPlateFile);

// N-Plate models need a model of a spacecraft's attitude, which is parameterized using
// body axes. The following defines the attitude of a nadir pointing satellite.
Vector nadirPointingVector = new VectorTrueDisplacement(satellite.IntegrationPoint,
                                                        Earth.CenterOfMassPoint);
Vector velocityVector = satellite.IntegrationPoint.CreateVectorVelocity(earth.InertialFrame);
Axes bodyAxes = new AxesAlignedConstrained(nadirPointingVector, AxisIndicator.Third,
                                           velocityVector, AxisIndicator.First);

// The density is defined the same way as the simpler AtmosphericDragForce model.
var density = new ScalarDensityJacchia1970(satellite.IntegrationPoint,
                                           new ConstantSolarGeophysicalData(150, 150.0, 3));

// The N-Plate drag model can now be constructed. It needs the Sun to be defined
// as its illuminating body in order to define the direction in which any
// solar panel plates are facing.
var nPlateDrag = new NPlateAtmosphericDragForce(nPlateModel, bodyAxes, density, sun);

  • N-Plate Solar Radiation Pressure

    In addition to the SimpleSolarRadiationForce model that is low-fidelity because it assumes a spherical spacecraft, a medium-fidelity model for solar radiation pressure called NPlateSolarRadiationForce is also available.

    This model splits the spacecraft into individual plates that have distinct areas, pointing directions, specular reflectivity coefficients, and diffuse reflectivity coefficients. These individual plates are modeled by SolarRadiationPlates that are stored in a NPlateModel. There are several different concrete SolarRadiationPlate types available:

    • A SolarRadiationBodyPlate points in a specific direction defined by its UnitNormalCartesian, which is defined in the spacecraft's BodyAxes. Solar radiation force is only computed for a body plate if it is facing toward the Sun. Self-shadowing effects from other plates are not taken into account if one plate is in front of another, but both are pointing toward the Sun. This is a limitation on the fidelity of the model.

    • A SolarRadiationTwoDofSolarPanelPlate points in the apparent direction of the Sun and its pointing is independent of the spacecraft's BodyAxes.

    • A SolarRadiationOneDofSolarPanelPlate points as close to the apparent direction of the Sun as it can while being constrained to only rotate around its UnitAxisCartesian, which is defined in the spacecraft's BodyAxes.

    A code example of how to build a NPlateSolarRadiationForce instance is given as follows:

    C#
    Copy
    // Getting the necessary bodies from the CentralBodiesFacet.
EarthCentralBody earth = CentralBodiesFacet.GetFromContext().Earth;
SunCentralBody sun = CentralBodiesFacet.GetFromContext().Sun;

// The partial configuration of a PropagationNewtonianPoint was demonstrated in an earlier recipe.
PropagationNewtonianPoint satellite = SeeRepresentASatellitePosition();

// An N-Plate model can be created from individual SolarRadiationPlates such as SolarRadiationBodyPlate, which
// model plates that point in specific directions in the spacecraft's body frame.
NPlateModel nPlateModel = new NPlateModel();
nPlateModel.Add(new SolarRadiationBodyPlate("PlusZ", "Body", 5.0, UnitCartesian.UnitZ, 0.1, 0.7));
nPlateModel.Add(new SolarRadiationBodyPlate("PlusY", "Body", 5.0, UnitCartesian.UnitY, 0.1, 0.7));
nPlateModel.Add(new SolarRadiationBodyPlate("PlusX", "Body", 5.0, UnitCartesian.UnitX, 0.1, 0.7));
nPlateModel.Add(new SolarRadiationBodyPlate("MinusZ", "Body", 5.0, -UnitCartesian.UnitZ, 0.1, 0.7));
nPlateModel.Add(new SolarRadiationBodyPlate("MinusY", "Body", 5.0, -UnitCartesian.UnitY, 0.1, 0.7));
nPlateModel.Add(new SolarRadiationBodyPlate("MinusX", "Body", 5.0, -UnitCartesian.UnitX, 0.1, 0.7));

// Another SolarRadiationPlate is the SolarRadiationTwoDofSolarPanelPlate, that models a solar panel that is
// able to always point directly toward the apparent Sun direction.
nPlateModel.Add(new SolarRadiationTwoDofSolarPanelPlate("SolarPanel", "SolarPanels", 60.0, 0.3, 0.5));

// The final SolarRadiationPlate is the SolarRadiationOneDofSolarPanelPlate, that models a solar panel that is
// only able to point toward the apparent Sun direction to the extent that its single rotation axis allows
// it to do so.
nPlateModel.Add(new SolarRadiationOneDofSolarPanelPlate("AxisPanel", "SolarPanels", 5.0,
                                                        UnitCartesian.UnitZ, 0.4, 0.4));

// However, N-Plate models are usually read from .nplate files.
nPlateModel = NPlateModel.ReadFrom(pathToNPlateFile);

// N-Plate models need a model of a spacecraft's attitude, which is parameterized using
// body axes. The following defines the attitude of a nadir pointing satellite.
Vector nadirPointingVector = new VectorTrueDisplacement(satellite.IntegrationPoint,
                                                        Earth.CenterOfMassPoint);
Vector velocityVector = satellite.IntegrationPoint.CreateVectorVelocity(earth.InertialFrame);
Axes bodyAxes = new AxesAlignedConstrained(nadirPointingVector, AxisIndicator.Third,
                                           velocityVector, AxisIndicator.First);

// The occultation is defined the same way as in the SimpleSolarRadiationForce,
// except that it cannot be corrected using SolarRadiationBoundaryMitigation.
// Thus, ScalarOccultationRegulatedDualCone should not be used.
var occultationFactor = new ScalarOccultationDualCone(Sun, satellite.IntegrationPoint,
                                                      Earth, Moon);

// The N-Plate solar radiation force model can now be constructed.
var nPlateSolarRadiationForce = new NPlateSolarRadiationForce(nPlateModel, bodyAxes, occultationFactor);
Variable Area Files

The content of variable area files, which are typically .dat files, can be read by the ReadFrom method of ScalarVariableArea. Depending on the type of variable area file that is read, this operation will result in either a ScalarVariableAreaArgumentLatitude instance or a ScalarVariableAreaJulianDate instance.

The ScalarVariableAreaArgumentLatitude type loosely models the differences in the atmosphere-facing or Sun-facing area of a spacecraft over the course of its orbit by interpolating its area based on the argument of latitude vs. area data given by the variable area file.

The ScalarVariableAreaJulianDate type loosely models the differences in the atmosphere-facing or Sun-facing area of a spacecraft over time by interpolating its area based on the time vs. area data given by the variable area file. There are two possible behaviors that are available for handling times before or after the range of the data in the file:

  • If CycleRepeats is true, the interval of the time data will be added or subtracted from the time until the updated time is inside the interval of the time data. Then, the area will be interpolated from the updated time.

  • If CycleRepeats is false, the area at the beginning of the time data will be used for all times before that time, and the area at the end of the time data will be used for all times after that time.

Either concrete ScalarVariableArea type can be used in the low-fidelity AtmosphericDragForce or SimpleSolarRadiationForce as the reference area property.

Custom Forces

Custom forces can be created by extending the ForceModel class and defining a ForceEvaluator which computes the force. See the Evaluators And Evaluator Groups topic for more information. Also, see the Torque Models topic for an example of how to create torque models, which are very similar to force models.

Some forces are specific, indicating that they already account for mass in their calculation. For instance, most gravity models usually express the specific force (acceleration) on a body, since the mass simply cancels. For the sake of coordinate transformations, specific forces are still treated as forces even though the PropagationNewtonianPoint does not divide them by mass in order to compute the acceleration.

There are some custom forces that will require some care when using them to propagate a trajectory, for example, thrust. It's simple for a user to create a thrust model by specifying a force and determining a mass flow rate (or vice versa) and then specifying a position element and a mass element to integrate in the state. However, care needs to be taken that the mass doesn't become negative after the fuel is exhausted. Generally, forces that only apply for short periods of time need to be handled slightly differently from forces that occur as a result of the environment. Such forces should only be applied in a NumericalPropagator over the span when the force is active and it will need to target the time at which the engine stops, after which a new propagator should be used. Otherwise, it's possible to adjust results with a PropagationStateCorrector, but in some cases this method can cause inaccuracies and other problems during propagation.

Force Models and Validation

DME Component Libraries is validated using an automated test suite that uses the most recent version of STK. There are a few things that users should know when attempting to compare values from the numerical propagator in DME Component Libraries against HPOP in STK:

  • Generally, it is best to start by integrating a trajectory with TwoBodyGravity as the only force model, to verify that the initial state, frame, and time standards are all identical before proceeding to more complex force models. In STK, you can do this by creating an STK HPOP satellite and turning off all of the force models except for gravity and then decreasing the degree and order to zero with tides turned off.

  • Both STK 11 and DME Component Libraries use InternationalCelestialReferenceFrame (ICRF) by default to define the Earth's InertialFrame. If you are comparing against earlier versions of STK, you may need to reassign the Earth's InertialFrame to be the J2000Frame. It is also important to note that propagation should generally occur in an inertial frame (either ICRF or J2000), unless the user intends to model coriolis and centrifugal effects separately.

  • By default, STK uses Earth Orientation Parameters (EOP) to define polar corrections to the orientation of the Earth's FixedFrame. To match exactly, DME Component Libraries must use the same EOP data as STK. See the External Data topic for more information on obtaining and configuring EOP data.

  • By default, STK 11 uses JPL DE430 ephemerides for the positions of the Sun, Moon, and other celestial bodies. To match with STK, you should use JplDE430 to update the central bodies in the calculation context, as well as set the MoonCentralBodyFixedFrame to be the topographic fixed frame, available from the GetMoonTopographicFixedFrame method. Without setting the lunar fixed frame, third body occlusion calculations may be slightly inaccurate when radiation is obscured by the Moon. See the External Data topic for more information on obtaining JPL Ephemerides and configuring them in DME Component Libraries.

  • In STK 11, the default constant value for solar luminosity is 3.839e+26 watts. Old versions of STK used a different default value. SimpleSolarRadiationForce in DME Component Libraries uses the same value as STK 11 as its default.

  • In STK 11, the default radius of the solar body is 6.95508e+8 meters. Old versions of STK used a different default value. SunCentralBody in DME Component Libraries uses the same value as STK 11 as its default.

  • Small differences in how the Earth's geoid Shape is represented and used to transform between geodetic and cartesian positions may cause slight discrepancies in ephemeris between DME Component Libraries and STK. The iteration involved in computing geodetic positions causes numerical noise which will create small discrepancies between two different algorithms and will grow over several integration steps.

  • Care should be taken when noting settings on types requiring vectors to the Sun or other celestial bodies. Force models will often provide options indicating whether to use apparent or true vectors. STK 11 and DME Component Libraries share the same defaults and in most cases use apparent vectors accounting for light time delay and aberration where radiation is involved.

  • STK uses a ratio of area/mass to denote constants for drag and solar radiation. Due to the modular nature of DME Component Libraries, the constants are specified directly, leaving the choice of whether to use a scalar ratio up to the user. So, for the STK defaults (a ratio of 0.02 m2/kg for a 1000 kg spacecraft), the corresponding reference area would be 20 m2. You could also create a Scalar representing the STK ratio value and multiply by the Scalar representing mass.

  • Note that when comparing values between any two sources, care needs to be taken that interpolation schemes are not inserting additional numerical error above and beyond that associated with propagation. To match precisely, it's important that the values be compared at the same time, in the same time standard, and that time needs to be EXACTLY at an ephemeris point.

  • Using a RungeKuttaFehlberg78Integrator using a fixed 60 second step size, with SphericalHarmonicGravity, ThirdBodyGravity to model solar and lunar perturbations, ScalarDensityJacchiaRoberts, SimpleSolarRadiationForce using ScalarOccultationDualCone, and a lot of careful management of numerical noise in initial state and time, the correspondence to STK 11 over one day is well below a millimeter of error and in some cases is on the order of a micrometer, depending on the test orbit used.

  • Lastly, if you do decide to attempt your own validation of DME Component Libraries and run into discrepancies, please contact AGI support. Chances are there are settings that are causing the discrepancy. However, in many cases it is important to remember that small "discrepancies" between DME Component Libraries and STK may simply be the accumulation of numerical noise and still be well below the accuracy of the force models themselves, most especially ones involving models of the Earth's atmosphere.