Code Sample

The following code samples illustrate the various pieces required to setup, configure, and numerically propagate an aircraft through multiple maneuver segments. The propagation sequence has the aircraft hold a constant heading, pull up to a desired rate of climb, reach an altitude at which to begin leveling off, push over to zero rate of climb, and turn left at a constant rate of change in heading.

Note
The functionality described in this topic requires a license for the Aircraft Propagation Library and the Segment Propagation Library.

Height Reference

Some of the objects used in the Aircraft Propagation Library require a height reference surface, which defaults to EarthCentralBodyMeanSeaLevel, which is the Mean Sea Level (MSL) surface of the Earth. This surface is initially undefined and must be configured before using the library. See the Terrain topic for more information on obtaining and loading MSL data. Alternatively, an EllipsoidTerrainProvider can be assigned to the MeanSeaLevel property to use the Earth ellipsoidal shape as an inaccurate approximation of the mean sea level surface, without additional data files.

Atmosphere, Wind, and Performance Models

Next, we define the atmosphere and wind models to be used. USStandardAtmosphere1976 is used as the atmospheric model in this example. For wind, we use a ConstantWindModel and define the local horizontal wind speed and direction (from North) and up/down draft. It is important to note that the wind direction indicates where the wind is coming from, as would the arrowhead of a weather vane.

Copy

USStandardAtmosphere1976 atmosphere = new USStandardAtmosphere1976();

ConstantWindModel winds = new ConstantWindModel
{
    CentralBody = earth,
    HorizontalWindSpeed = 10.0, // meters per second
    WindDirection = Trig.DegreesToRadians(90.0), // convert to radians, winds are coming from the East
};

We also configure the set of performance models corresponding to the various flight phases. These performance models govern the aircraft's motion during the maneuvers. The load limits are specified as G's with the EarthSurfaceGravity used to represent standard gravity for conversion to accelerations.

Copy

// Construct a cruise performance model which varies with height relative to mean sea level.
BandedCruisePerformanceModel cruise = new BandedCruisePerformanceModel(earth.MeanSeaLevel);
cruise.AddBand(10000.0, new CruiseCommandedValues(150.0, 0.0755987));
cruise.AddBand(11000.0, new CruiseCommandedValues(175.0, 0.0755987));

PerformanceModels performance = new PerformanceModels
{
    Acceleration = new SimpleAccelerationPerformanceModel
    {
        LongitudinalLoadLimits = new Bounds(-0.5, 0.5),  // Deceleration and acceleration in G's
        LateralLoadLimits = new Bounds(-1.1547, 1.1547), // Left and right turn acceleration in G's
        VerticalLoadLimits = new Bounds(0.75, 1.1547),   // Push over and pull up acceleration in G's
        GravitationalAcceleration = Constants.EarthSurfaceGravity,
    },
    Cruise = cruise,
    Climb = new SimpleClimbPerformanceModel(92.6, 10.16, 0.0629989),
    Descent = new SimpleDescentPerformanceModel(92.6, -10.16, 0.0629989),
};

Defining the First Segment

We define the initial position, velocity, and mass of the aircraft. Also, we configure a numerical integrator for use in propagating the state.

Copy

double longitude = Trig.DegreesToRadians(0.0); // convert to radians
double latitude = Trig.DegreesToRadians(0.0); // convert to radians
double height = 10000.0; // meters

double initialMass = 1500.0; // kilograms

AircraftReferenceState initialState = new AircraftReferenceState(epoch)
{
    HeightReferenceSurface = earth.MeanSeaLevel,              // This is the default, but we'll set it here for clarity.
    Position = new Cartographic(longitude, latitude, height), // Then, this initial height is altitude with respect to mean sea level.
    VelocityReferenceAxes = ManeuverReferenceAxes.MovingAtmosphere,                    // The velocity is relative to the moving atmosphere,
    Velocity = new AzimuthHorizontalVertical(Trig.DegreesToRadians(45.0), 150.0, 0.0), // so these values are true heading and true airspeed.
};

NumericalIntegrator integrator = new RungeKutta4Integrator
{
    InitialStepSize = 1.0,
};

Our first maneuver will hold a constant angle and speed using a MaintainCourseOrHeadingBehavior as its horizontal behavior and a MaintainHorizontalAndVerticalSpeedsBehavior as its vertical behavior. The ManeuverReferenceAxes is set to the moving atmosphere, so the true heading of the aircraft is the angle that is maintained. If the ManeuverReferenceAxes were set to the static atmosphere, then the true course of the aircraft would be the angle that is maintained.

Copy

CompositeManeuver flyStraight = new CompositeManeuver
{
    PerformanceModels = performance,
    WindModel = winds,
    HorizontalBehavior = new MaintainCourseOrHeadingBehavior(ManeuverReferenceAxes.MovingAtmosphere),           // Maintain true heading.
    VerticalBehavior = new MaintainHorizontalAndVerticalSpeedsBehavior(ManeuverReferenceAxes.MovingAtmosphere), // Maintain true air velocities.
};

Next, some configuration is needed to set up the propagator. The state elements to be propagated are the aircraft's Cartesian motion and mass. Each are configured with an initial value and an identifier for use later in recovering the results of the propagation. The aircraft motion element's mass property is set to the mass state element's integrated value so that the correct value is used during numerical integration.

The propagator definition is configured with the numerical integration method and the aircraft motion element and mass state element are added to the list of elements to be integrated. The ApplyManeuver method creates the state element which will be the source of the derivatives for the motion element during numerical integration.

The GetFuelFlowRateScalar helper function provides the scalar derivative for the rate of change of mass. The details of this function will be discussed next, after we cover the final details of setting up the first propagation segment.

Copy

// Define the mass state element for the propagator definition.
string massID = "MyAircraftMass";
PropagationScalar massStateElement = new PropagationScalar
{
    Identification = massID,
    InitialState = new Motion<double>(initialMass),
    IncludeHighestDerivativeInOutput = true, // Include the mass flow rate information.
};

// Define the position state element for the propagator definition.
string positionID = "MyAircraftPosition";
AircraftMotionElement aircraftMotion = new AircraftMotionElement(positionID, initialState, winds)
{
    Mass = massStateElement.IntegrationValue,
};

// Configure the propagator definition with the state elements.
NumericalPropagatorDefinition firstPropagatorDefinition = new NumericalPropagatorDefinition
{
    Epoch = epoch,
    Integrator = integrator,
};
firstPropagatorDefinition.IntegrationElements.Add(aircraftMotion.ApplyManeuver(flyStraight));
firstPropagatorDefinition.IntegrationElements.Add(massStateElement);

// Define the derivative function for the mass state element.
massStateElement.ScalarDerivative = GetFuelFlowRateScalar(atmosphere, winds, aircraftMotion);

The first segment is configured with the propagator definition and a stopping condition which is triggered once the aircraft has flown for 5 minutes.

Copy

NumericalPropagatorSegment flyStraightFor5Minutes = new NumericalPropagatorSegment
{
    Name = "FlyStraightFor5Minutes",
    PropagatorDefinition = firstPropagatorDefinition
};
flyStraightFor5Minutes.StoppingConditions.Add(new DurationStoppingCondition(Duration.FromMinutes(5))
{
    Name = "StopIn5Minutes",
});

Fuel Flow Rate

We now return to the mass flow rate mentioned above. In this example, the mass flow rate is derived from a simple propeller propulsion model. In order to determine the fuel flow rate, the applied thrust must be determined, which in turn means that the propulsion and aerodynamic forces acting on the aircraft must be computed, which requires that the orientation of the aircraft body axes must be assessed. Fortunately, the SimpleFixedWingCoordinatedFlight type can be used to facilitate these computations.

The GetFuelFlowRateScalar helper function configures the aerodynamic model to be used with the required parameters and a ScalarDynamicPressure. The SimpleFixedWingCoordinatedFlight instance uses the aerodynamic model to determine the orientation of the aircraft (bank angle and angle of attack) and aerodynamic and propulsion forces required to follow the motion prescribed by the maneuver. The required thrust determined by the coordinated flight condition is used to configure the propulsion model so that the corresponding fuel flow rate can be determined. The scalar which can produce this fuel flow rate is returned from the GetFuelFlowRateScalar helper function for use as the derivative of the mass state element.

Copy

private Scalar GetFuelFlowRateScalar(USStandardAtmosphere1976 atmosphere, WindModel winds, AircraftMotionElement element)
{
    // The relative velocity vector accounts for the local winds at the desired location.
    // The dynamic pressure is needed to compute the aerodynamic forces.
    ScalarDynamicPressure dynamicPressure = new ScalarDynamicPressure(atmosphere.Density, winds.RelativeVelocityVector);

    SimpleFixedWingForwardFlightAerodynamics fixedWingAerodynamics = new SimpleFixedWingForwardFlightAerodynamics
    {
        ReferenceArea = 5.0,
        AspectRatio = 4.0,
        LiftCoefficientAtZeroAngleOfAttack = 0.0,
        LiftCoefficientSlope = 0.065 * Constants.DegreesPerRadian,
        OswaldEfficiencyFactor = 0.8,
        DragCoefficientAtMinimum = 0.02,
        DragIndex = 0.0,
        DynamicPressure = dynamicPressure.ApplyServiceProvider(element.IntegrationPoint),
    };

    SimpleFixedWingCoordinatedFlight flightCondition = new SimpleFixedWingCoordinatedFlight(element.IntegrationPoint, winds, fixedWingAerodynamics, element.Mass, Constants.EarthSurfaceGravity);

    SimpleForwardFlightPropellerPropulsion propulsionModel = new SimpleForwardFlightPropellerPropulsion
    {
        Density = atmosphere.Density,
        DensityRatio = atmosphere.DensityRatio,
        FuelFlowRateLimits = new Bounds(0.022167, 0.22167),
        PowerAvailableLimits = new Bounds(0.0, 1.70064e+006),
        PowerDropOffParameter = 0.117,
        PropellerCount = 1,
        PropellerRadius = 3.048,
        ReferencePoint = element.IntegrationPoint,
        RelativeVelocity = winds.RelativeVelocityVector,
        ThrustRequired = flightCondition.ThrustRequired,
    };

    return propulsionModel.FuelFlowRate;
}

Defining the Remaining Segments

Returning to the propagation set-up, the second segment is configured similarly to the first. This segment performs a pull-up maneuver while holding a constant heading. The segment stops once the rate of climb is 50 meters per second.

Copy

CompositeManeuver pullUp = new CompositeManeuver
{
    PerformanceModels = performance,
    WindModel = winds,
    HorizontalBehavior = new MaintainCourseOrHeadingBehavior(ManeuverReferenceAxes.MovingAtmosphere), // Maintain true heading.
    VerticalBehavior = new PushOverOrPullUpBehavior
    {
        TransverseLoadFactor = 1.15, // This is a pull-up since the load factor is greater than 1 G.
        LongitudinalLoadFactor = 0.0, // The velocity does not increase during the pull-up.
        GravitationalAcceleration = Constants.EarthSurfaceGravity,
    },
};

NumericalPropagatorDefinition secondPropagatorDefinition = new NumericalPropagatorDefinition
{
    Epoch = epoch,
    Integrator = integrator,
};
secondPropagatorDefinition.IntegrationElements.Add(aircraftMotion.ApplyManeuver(pullUp));
secondPropagatorDefinition.IntegrationElements.Add(massStateElement);

NumericalPropagatorSegment pullUpTo50MpsClimb = new NumericalPropagatorSegment
{
    Name = "PullUpTo50MpsClimb",
    PropagatorDefinition = secondPropagatorDefinition
};
pullUpTo50MpsClimb.StoppingConditions.Add(new ScalarStoppingCondition
{
    Name = "StopAt50MpsClimb",
    Scalar = new ScalarDerivative(new ScalarCartographicElement(earth, aircraftMotion.IntegrationPoint, CartographicElement.Height), 1),
    Threshold = 50.0,
    TypeOfStoppingCondition = StopType.AnyThreshold,
    FunctionTolerance = 0.001,
});

A constant rate of climb is held for the third segment until the aircraft reaches an altitude of 11000 meters.

Copy

CompositeManeuver climb = new CompositeManeuver
{
    PerformanceModels = performance,
    WindModel = winds,
    HorizontalBehavior = new MaintainCourseOrHeadingBehavior(ManeuverReferenceAxes.MovingAtmosphere),           // Maintain true heading.
    VerticalBehavior = new MaintainHorizontalAndVerticalSpeedsBehavior(ManeuverReferenceAxes.MovingAtmosphere), // Maintain true air velocities.
};

NumericalPropagatorDefinition thirdPropagatorDefinition = new NumericalPropagatorDefinition
{
    Epoch = epoch,
    Integrator = integrator,
};
thirdPropagatorDefinition.IntegrationElements.Add(aircraftMotion.ApplyManeuver(climb));
thirdPropagatorDefinition.IntegrationElements.Add(massStateElement);

NumericalPropagatorSegment climbTo11KmAltitude = new NumericalPropagatorSegment
{
    Name = "ClimbTo11KmAltitude",
    PropagatorDefinition = thirdPropagatorDefinition
};
climbTo11KmAltitude.StoppingConditions.Add(new ScalarStoppingCondition
{
    Name = "StopAt11KmAltitude",
    Scalar = new ScalarCartographicElement(earth, aircraftMotion.IntegrationPoint, CartographicElement.Height),
    Threshold = 11000.0,
    TypeOfStoppingCondition = StopType.AnyThreshold,
    FunctionTolerance = 0.001,
});

The aircraft then performs a push-over maneuver until a zero rate of climb is attained.

Copy

CompositeManeuver pushOver = new CompositeManeuver
{
    PerformanceModels = performance,
    WindModel = winds,
    HorizontalBehavior = new MaintainCourseOrHeadingBehavior(ManeuverReferenceAxes.MovingAtmosphere), // Maintain true heading.
    VerticalBehavior = new PushOverOrPullUpBehavior
    {
        TransverseLoadFactor = 0.75, // This is a push-over since the load factor is less than 1 G.
        LongitudinalLoadFactor = 0.0, // The velocity does not increase during the push-over.
        GravitationalAcceleration = Constants.EarthSurfaceGravity,
    },
};

NumericalPropagatorDefinition fourthPropagatorDefinition = new NumericalPropagatorDefinition
{
    Epoch = epoch,
    Integrator = integrator,
};
fourthPropagatorDefinition.IntegrationElements.Add(aircraftMotion.ApplyManeuver(pushOver));
fourthPropagatorDefinition.IntegrationElements.Add(massStateElement);

NumericalPropagatorSegment pushOverTo0MpsClimb = new NumericalPropagatorSegment
{
    Name = "PushOverTo0MpsClimb",
    PropagatorDefinition = fourthPropagatorDefinition
};
pushOverTo0MpsClimb.StoppingConditions.Add(new ScalarStoppingCondition
{
    Name = "StopAt0MpsClimb",
    Scalar = new ScalarDerivative(new ScalarCartographicElement(earth, aircraftMotion.IntegrationPoint, CartographicElement.Height), 1),
    Threshold = 0.0,
    TypeOfStoppingCondition = StopType.AnyThreshold,
    FunctionTolerance = 0.001,
});

Finally, the aircraft enters a left turn where the heading changes by 1 degree per second. This turn is maintained for 5 minutes.

Copy

CompositeManeuver turnLeft = new CompositeManeuver
{
    PerformanceModels = performance,
    WindModel = winds,
    HorizontalBehavior = new TurnAtConstantRateBehavior
    {
        ManeuverReferenceAxes = ManeuverReferenceAxes.MovingAtmosphere, // The true heading will change at the constant rate.
        TurnRate = -Trig.DegreesToRadians(1.0),                         // Left turn since the rate is negative.
    },
    VerticalBehavior = new CruiseBehavior
    {
        HeightReferenceSurface = earth.MeanSeaLevel, // This is the default, but we'll set it here for clarity.
        TimeConstant = 1.0,
        Altitude = 11000.0,                          // Then, this cruise height is altitude with respect to mean sea level.
        LevelOffAltitudeBandAboveCruise = 100.0,     // Above this band the descent performance model is used.
        LevelOffAltitudeBandBelowCruise = 100.0,     // Below this band the climb performance model is used.
    },
};

NumericalPropagatorDefinition fifthPropagatorDefinition = new NumericalPropagatorDefinition
{
    Epoch = epoch,
    Integrator = integrator,
};
fifthPropagatorDefinition.IntegrationElements.Add(aircraftMotion.ApplyManeuver(turnLeft));
fifthPropagatorDefinition.IntegrationElements.Add(massStateElement);

NumericalPropagatorSegment turnLeftFor5Minutes = new NumericalPropagatorSegment
{
    Name = "TurnLeftFor5Minutes",
    PropagatorDefinition = fifthPropagatorDefinition
};
turnLeftFor5Minutes.StoppingConditions.Add(new DurationStoppingCondition(Duration.FromMinutes(5))
{
    Name = "StopIn5Minutes",
});

Propagating

All five segments are added to a SegmentList and the sequence is propagated.

The position element identifier is used to retrieve the aircraft's motion from the propagation results. Propagation of the maneuver segments produces ephemeris for the aircraft, which defines its position, velocity, and acceleration at discrete times. In order to provide information at an arbitrary time over the propagation, the motion must be interpolated. The code below uses 5th-order Hermite interpolation.

Copy

SegmentList segmentList = new SegmentList
{
    Name = "StraightThenClimbThenTurn"
};

segmentList.Segments.Add(flyStraightFor5Minutes);
segmentList.Segments.Add(pullUpTo50MpsClimb);
segmentList.Segments.Add(climbTo11KmAltitude);
segmentList.Segments.Add(pushOverTo0MpsClimb);
segmentList.Segments.Add(turnLeftFor5Minutes);

SegmentListResults results = segmentList.GetSegmentListPropagator().PropagateSegmentList();

// Recover the position, velocity, and acceleration.
DateMotionCollection<Cartesian> ephemeris = results.GetDateMotionCollectionOfOverallTrajectory<Cartesian>(positionID);

// Recover the mass and flow rate.
DateMotionCollection<double> mass = results.GetDateMotionCollectionOfOverallTrajectory<double>(massID);

Some additional scalars are defined from the interpolated point in order to compute the aircraft's altitude, true airspeed, heading, and course.

Copy

EvaluatorGroup group = new EvaluatorGroup();
PointInterpolator location = new PointInterpolator(earth.FixedFrame, new HermitePolynomialApproximation(), 5, ephemeris);

// The fixed velocity vector is needed to compute the true course of the aircraft.
Vector fixedVelocity = new VectorVelocity(location, earth.FixedFrame);

// The relative velocity vector accounts for the local winds at the desired location.
// It will be used to compute the true heading and true airspeed.
Vector relativeVelocity = winds.RelativeVelocityVector.ApplyServiceProvider(location);

MotionEvaluator<Cartographic> cartographicEvaluator = earth.ObserveCartographicPoint(location, group);

Axes eastNorthUpAxes = new AxesEastNorthUp(earth, location);

ScalarEvaluator trueAirspeedEvaluator = new ScalarVectorElement(relativeVelocity, CartesianElement.Magnitude, eastNorthUpAxes).GetEvaluator(group);

ScalarEvaluator trueHeadingEvaluator = new ScalarInverseTangent(new ScalarVectorElement(relativeVelocity, CartesianElement.Y, eastNorthUpAxes),
                                                                new ScalarVectorElement(relativeVelocity, CartesianElement.X, eastNorthUpAxes)).GetEvaluator(group);

ScalarEvaluator trueCourseEvaluator = new ScalarInverseTangent(new ScalarVectorElement(fixedVelocity, CartesianElement.Y, eastNorthUpAxes),
                                                               new ScalarVectorElement(fixedVelocity, CartesianElement.X, eastNorthUpAxes)).GetEvaluator(group);

ScalarInterpolator interpolatedMass = new ScalarInterpolator(new HermitePolynomialApproximation(), 1, mass);

ScalarEvaluator massEvaluator = interpolatedMass.GetEvaluator(group);

With some assumptions about how the aircraft is configured and flown, its orientation can be determined. Specifically, assuming that the thrust acts along the longitudinal axis (x-axis) of the vehicle and that there is no slip angle (and correspondingly, no aerodynamic side force), the bank angle and angle of attack can be determined from an aerodynamic force model.

In the following code sample, a SimpleFixedWingForwardFlightAerodynamics model is applied to define the orientation of the body axes of the aircraft. The force model is configured with a ScalarDynamicPressure derived from a density and velocity relative to the moving atmosphere.

Copy

// The dynamic pressure is needed to compute the aerodynamic forces.
ScalarDynamicPressure q = new ScalarDynamicPressure(atmosphere.Density, relativeVelocity);

SimpleFixedWingForwardFlightAerodynamics aerodynamics = new SimpleFixedWingForwardFlightAerodynamics
{
    ReferenceArea = 5.0,
    AspectRatio = 4.0,
    LiftCoefficientAtZeroAngleOfAttack = 0.0,
    LiftCoefficientSlope = 0.065 * Constants.DegreesPerRadian,
    OswaldEfficiencyFactor = 0.8,
    DragCoefficientAtMinimum = 0.02,
    DragIndex = 0.0,
    DynamicPressure = q.ApplyServiceProvider(location),
};

SimpleFixedWingCoordinatedFlight coordinatedFlight = new SimpleFixedWingCoordinatedFlight(location, winds, aerodynamics, interpolatedMass, Constants.EarthSurfaceGravity);

Axes bodyAxes = coordinatedFlight.BodyAxes;