Lunar Free Return Code Sample

The two primary classes of Lunar free return trajectories are categorized in terms of whether their Lunar flyby occurs between the Earth and the Moon (cislunar) or past the Moon (circumlunar). Their design is so similar that most of the code in this example is shared between the two classes. The main practical difference between the two classes is that the cislunar free return has a much longer time of flight: about 13 days from the Earth to the Moon and about 13 days to return to the Earth. The circumlunar free return takes about 3 days to get to the Moon and 3 days to return to the Earth.

Free return trajectories are symmetric in the circular, restricted three-body problem, but the the Solar gravitational perturbations and spherical harmonic gravity models used in this example break that symmetry. However, this symmetry can still be used to inform a mission design strategy.

A rotating frame can be defined using the Earth-Moon line and the normal vector of the Moon's orbit about the Earth. Initial guesses for the flybys can be defined easily in terms of this rotating frame. These initial guesses can be used to backward target initial guesses for Earth parking orbit and translunar injection maneuver conditions. The Earth conditions can be varied in order to forward target the desired Earth re-entry conditions with an implicit Lunar flyby. In this example, the Lunar flyby altitudes are higher than the initial guess which is not necessarily going to be true in the general case. The flyby initial guesses may need to have their altitudes increased to ensure safe lunar flyby conditions occur.

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

Free Return Initial Setup

The initial setup code includes common setup that is needed for both circumlunar and cislunar free return targeting.

Copy

// Pick names.
const string motionID = "Satellite";
const string dryMassID = "Dry_Mass";

// Use JPL DE 440 ephemerides for positions of Sun, Moon, and Earth.
JplDE440 ephem = new JplDE440(pathToJplDe440);
ephem.UseForCentralBodyPositions(CentralBodiesFacet.GetFromContext());

// Store commonly used central bodies and constants.
CentralBody moon = CentralBodiesFacet.GetFromContext().Moon;
SphericalHarmonicGravityModel moonGravityModel = SphericalHarmonicGravityModel.ReadFrom(pathToMoonGravFile);
double moonGravitationalParameter = moonGravityModel.GravitationalParameter; // m^3/s^2.
double moonPhysicalRadius = moonGravityModel.ReferenceDistance; // Meters.

CentralBody sun = CentralBodiesFacet.GetFromContext().Sun;
const double sunGravitationalParameter = 1.3271244004193938e20; // m^3/s^2.

CentralBody earth = CentralBodiesFacet.GetFromContext().Earth;
SphericalHarmonicGravityModel earthGravityModel = SphericalHarmonicGravityModel.ReadFrom(pathToEarthGravFile);
double earthGravitationalParameter = earthGravityModel.GravitationalParameter;
double earthPhysicalRadius = earthGravityModel.ReferenceDistance; // Meters.

// Set up rotating reference frames.
VectorTrueDisplacement earthMoonVector = new VectorTrueDisplacement(earth.CenterOfMassPoint, moon.CenterOfMassPoint);
VectorVelocity moonOrbitalVelocity = new VectorVelocity(moon.CenterOfMassPoint, Earth.InertialFrame);
VectorCrossProduct moonOrbitalAngularMomentum = new VectorCrossProduct(earthMoonVector, moonOrbitalVelocity);
AxesAlignedConstrained rotatingAxes = new AxesAlignedConstrained(earthMoonVector, AxisIndicator.First, moonOrbitalAngularMomentum, AxisIndicator.Third);

ReferenceFrame moonCenteredRotatingFrame = new ReferenceFrame(moon.CenterOfMassPoint, rotatingAxes);

double earthMoonMassRatio = earthGravitationalParameter / moonGravitationalParameter;
double scaleFactor = 1.0 / (1.0 + earthMoonMassRatio);
VectorScaled earthToBarycenterVector = new VectorScaled(earthMoonVector, scaleFactor);
Point barycenter = new PointVectorToPoint(earthToBarycenterVector, earth.CenterOfMassPoint);
ReferenceFrame barycentricRotatingFrame = new ReferenceFrame(barycenter, rotatingAxes);

// Get reference frame transformation from Moon rotating frame to Moon inertial frame.
EvaluatorGroup group = new EvaluatorGroup();
ReferenceFrameEvaluator fromMoonRotatingToMoonInertial = GeometryTransformer.GetReferenceFrameTransformation(moonCenteredRotatingFrame, moon.InertialFrame, group);

// Initial guess date for Lunar flyby.
JulianDate epoch = new JulianDate(new GregorianDate(2023, 1, 1, 0, 0, 0.0), TimeStandard.CoordinatedUniversalTime);

// Get transformations at epoch. Order 1 allows velocity information to be transformed.
KinematicTransformation fromMoonRotatingAtEpoch = fromMoonRotatingToMoonInertial.Evaluate(epoch, 1);
KinematicTransformation toMoonRotatingAtEpoch = fromMoonRotatingAtEpoch.Invert(1);

Cartesian initialPositionMoonRotating;
Cartesian initialVelocityMoonRotating;

Circumlunar Initial Guess

The initial guess for the circumlunar free return trajectory is a Lunar flyby that occurs 100 km past the Moon (along the line from the Earth to the Moon).

Copy

// Initial state is in moon-centered rotating coordinates to ensure free return conditions.
// Units are meters, so this is consistent with a 100 km circumlunar Lunar flyby.
initialPositionMoonRotating = new Cartesian(moonPhysicalRadius + 100000.0, 0.0, 0.0);
initialVelocityMoonRotating = new Cartesian(0.0, -2400.0, 0.0);

Cislunar Initial Guess

The initial guess for the cislunar free return trajectory is a Lunar flyby that occurs 100 km before the Moon (along the line from the Earth to the Moon).

Copy

// Initial state is in moon-centered rotating coordinates to ensure free return conditions.
// Units are meters, so this is consistent with a 100 km cislunar Lunar flyby.
initialPositionMoonRotating = new Cartesian(-moonPhysicalRadius - 100000.0, 0.0, 0.0);
initialVelocityMoonRotating = new Cartesian(0.0, -2400.0, 0.0);

Free Return Backward Propagation Setup

This code is common to both the cislunar and the circumlunar examples. It sets up the force models, propagators, segments, and differential correctors that are used to perform the backward targeting of an initial guess for an Earth parking orbit from the given initial guess for the Lunar flyby conditions. It also runs those differential correctors to actually compute the desired initial guess.

Copy

Motion<Cartesian> initialStateMoonRotating = new Motion<Cartesian>(initialPositionMoonRotating, initialVelocityMoonRotating);

Motion<Cartesian> initialStateMoonInertial = fromMoonRotatingAtEpoch.Transform(initialStateMoonRotating, 1);

// Dry mass only because translunar injection maneuver is modeled without
// using any propellant mass expenditure model.
PropagationScalar dryMass = new PropagationScalar(1000.0);
dryMass.Identification = dryMassID;
dryMass.ScalarDerivative = new ScalarFixed(0.0);

// Create the integrator.
var integrator = new RungeKuttaFehlberg78Integrator();
integrator.InitialStepSize = 60.0;
integrator.MinimumStepSize = 1.0;
integrator.MaximumStepSize = 86400.0;

// Need to define a propagator for the Moon's sphere of influence.
PropagationNewtonianPoint propagationPointAroundMoon = new PropagationNewtonianPoint();
propagationPointAroundMoon.Identification = motionID;

propagationPointAroundMoon.InitialPosition = initialStateMoonInertial.Value;
propagationPointAroundMoon.InitialVelocity = initialStateMoonInertial.FirstDerivative;
propagationPointAroundMoon.Mass = dryMass.IntegrationValue;

propagationPointAroundMoon.IntegrationFrame = moon.InertialFrame;

// Moon is the primary gravity field.
SphericalHarmonicGravityField moonGravityField = new SphericalHarmonicGravityField(moonGravityModel, 4, 4, true, null);

SphericalHarmonicGravity moonPrimaryGravity = new SphericalHarmonicGravity(propagationPointAroundMoon.IntegrationPoint, moonGravityField);

propagationPointAroundMoon.AppliedForces.Add(moonPrimaryGravity);

// The Earth and Sun produce third body gravity for Moon propagator.
ThirdBodyGravity thirdBodyForMoonPropagator = new ThirdBodyGravity(propagationPointAroundMoon.IntegrationPoint);
thirdBodyForMoonPropagator.AddThirdBody(earth.Name, earth.CenterOfMassPoint, earthGravitationalParameter);
thirdBodyForMoonPropagator.AddThirdBody(sun.Name, sun.CenterOfMassPoint, sunGravitationalParameter);
thirdBodyForMoonPropagator.CentralBody = moon;

propagationPointAroundMoon.AppliedForces.Add(thirdBodyForMoonPropagator);

// Make the actual propagator.
NumericalPropagatorDefinition nearMoonNumericalPropagator = new NumericalPropagatorDefinition();

// Add the elements.
nearMoonNumericalPropagator.IntegrationElements.Add(propagationPointAroundMoon);
nearMoonNumericalPropagator.IntegrationElements.Add(dryMass);

// Set the integrator.
nearMoonNumericalPropagator.Integrator = integrator;

// Set an epoch.
nearMoonNumericalPropagator.Epoch = epoch;

// Make an initial state segment that uses this propagator.
NumericalInitialStateSegment initialStateSegment = new NumericalInitialStateSegment();
initialStateSegment.Name = "Lunar_Closest_Approach";
initialStateSegment.PropagatorDefinition = nearMoonNumericalPropagator;

// Make a segment that propagates backwards to edge of Lunar sphere of influence (SOI).
NumericalPropagatorSegment propagateBackwardToSoi = new NumericalPropagatorSegment();
propagateBackwardToSoi.Name = "Propagate_Backward_ToSOI";
propagateBackwardToSoi.PropagatorDefinition = nearMoonNumericalPropagator;
propagateBackwardToSoi.PropagationDirection = IntegrationSense.Decreasing;
propagateBackwardToSoi.StopOnMaximumDurationBehavior = MaximumDurationBehavior.Throw;

// Stop propagating this segment at the Lunar SOI boundary.
const double moonSemimajorAxis = 384400000.0; // Meters from ssd.jpl.nasa.gov/sats/elem/

// Calculated value is about 66182924 meters.
double moonSphereOfInfluence = moonSemimajorAxis * Math.Pow(moonGravitationalParameter / earthGravitationalParameter, 0.4);

ScalarStoppingCondition moonSoiStoppingCondition = new ScalarStoppingCondition(new ScalarPointElement(propagationPointAroundMoon.IntegrationPoint,
                                                                                                      CartesianElement.Magnitude,
                                                                                                      moon.InertialFrame),
                                                                               moonSphereOfInfluence,
                                                                               0.0001, // Value tolerance, meters.
                                                                               StopType.AnyThreshold);
moonSoiStoppingCondition.Name = "Moon_SOI";

propagateBackwardToSoi.StoppingConditions.Add(moonSoiStoppingCondition);

// Make an Earth-centered propagator for any propagation outside the Lunar SOI.
// The initial positions and velocity will be provided by the preceding segment,
// so it is okay to use Cartesian.Undefined here.
PropagationNewtonianPoint propagationPointAroundEarth = new PropagationNewtonianPoint(motionID, earth.InertialFrame, Cartesian.Undefined, Cartesian.Undefined);

// Total mass.
propagationPointAroundEarth.Mass = dryMass.IntegrationValue;

// Gravitational field for Earth.
SphericalHarmonicGravityField earthGravityField = new SphericalHarmonicGravityField(earthGravityModel, 4, 4, true, new PermanentSolidTideModel());
SphericalHarmonicGravity earthGravityAsPrimary = new SphericalHarmonicGravity(propagationPointAroundEarth.IntegrationPoint, earthGravityField);
propagationPointAroundEarth.AppliedForces.Add(earthGravityAsPrimary);

// The gravity of the Sun and Moon are modeled as third bodies for the Earth-centered propagator.
ThirdBodyGravity thirdBodyForEarthPropagator = new ThirdBodyGravity();
thirdBodyForEarthPropagator.TargetPoint = propagationPointAroundEarth.IntegrationPoint;
thirdBodyForEarthPropagator.AddThirdBody(moon.Name, moon.CenterOfMassPoint, moonGravitationalParameter);
thirdBodyForEarthPropagator.AddThirdBody(sun.Name, sun.CenterOfMassPoint, sunGravitationalParameter);
propagationPointAroundEarth.AppliedForces.Add(thirdBodyForEarthPropagator);

// Make the actual propagator.
NumericalPropagatorDefinition earthCenteredNumericalPropagator = new NumericalPropagatorDefinition();

// Add the elements.
earthCenteredNumericalPropagator.IntegrationElements.Add(propagationPointAroundEarth);
earthCenteredNumericalPropagator.IntegrationElements.Add(dryMass);

// Set the integrator.
earthCenteredNumericalPropagator.Integrator = integrator;

// Set epoch.
earthCenteredNumericalPropagator.Epoch = epoch;

// Make a segment that propagates backward two days using Earth-centered propagator.
NumericalPropagatorSegment propagateBackwardTwoDays = new NumericalPropagatorSegment();
propagateBackwardTwoDays.Name = "Propagate_Backward_TwoDays";
propagateBackwardTwoDays.PropagatorDefinition = earthCenteredNumericalPropagator;
propagateBackwardTwoDays.PropagationDirection = IntegrationSense.Decreasing;
propagateBackwardTwoDays.StopOnMaximumDurationBehavior = MaximumDurationBehavior.Throw;
propagateBackwardTwoDays.StoppingConditions.Add(new DurationStoppingCondition(Duration.FromDays(-2.0)));

// Make a segment that propagates backward to perigee.
NumericalPropagatorSegment propagateBackwardToPerigee = new NumericalPropagatorSegment();
propagateBackwardToPerigee.Name = "Propagate_Backward_ToPerigee";
propagateBackwardToPerigee.PropagatorDefinition = earthCenteredNumericalPropagator;
propagateBackwardToPerigee.PropagationDirection = IntegrationSense.Decreasing;
propagateBackwardToPerigee.StopOnMaximumDurationBehavior = MaximumDurationBehavior.Throw;

ScalarStoppingCondition perigeeStoppingCondition1 = new ScalarStoppingCondition(new ScalarPointElement(propagationPointAroundEarth.IntegrationPoint,
                                                                                                       CartesianElement.Magnitude,
                                                                                                       earth.InertialFrame),
                                                                                0.0001, // Value tolerance, meters.
                                                                                StopType.LocalMinimum);
perigeeStoppingCondition1.Name = "Earth_Perigee_Backwards_Targeting";

propagateBackwardToPerigee.StoppingConditions.Add(perigeeStoppingCondition1);

// Set up segments to use backward propagation.
TargetedSegmentList backwardSegmentList = new TargetedSegmentList();
backwardSegmentList.Name = "Backward_Segment_List";

backwardSegmentList.Segments.Add(initialStateSegment);
backwardSegmentList.Segments.Add(propagateBackwardToSoi);
backwardSegmentList.Segments.Add(propagateBackwardTwoDays);
backwardSegmentList.Segments.Add(propagateBackwardToPerigee);

// Two variables are rotating frame y and z velocity directions.
// Need to transform from initial to rotating and back inside setter for variable.
SetVariableCallback<NumericalPropagatorState> setterVy = (variable, configuration) =>
{
    var currentMotion = configuration.GetMotion<Cartesian>(motionID);
    var rotatingMotion = toMoonRotatingAtEpoch.Transform(currentMotion, 1);
    double updatedVy = rotatingMotion.FirstDerivative.Y + variable;
    var updatedRotating = new Motion<Cartesian>(rotatingMotion.Value,
                                                new Cartesian(rotatingMotion.FirstDerivative.X, updatedVy, rotatingMotion.FirstDerivative.Z));
    var updatedInertial = fromMoonRotatingAtEpoch.Transform(updatedRotating, 1);
    configuration.ModifyMotion(motionID, updatedInertial);
};

SegmentPropagatorVariable velocityYVariable = initialStateSegment.CreateVariable(10.0, // Maximum step, m/s.
                                                                                 0.1, // Perturbation, m/s.
                                                                                 setterVy);

// Need to transform from initial to rotating and back inside setter for variable.
SetVariableCallback<NumericalPropagatorState> setterVz = (variable, configuration) =>
{
    var currentMotion = configuration.GetMotion<Cartesian>(motionID);
    var rotatingMotion = toMoonRotatingAtEpoch.Transform(currentMotion, 1);
    double updatedVz = rotatingMotion.FirstDerivative.Z + variable;
    var updatedRotating = new Motion<Cartesian>(rotatingMotion.Value,
                                                new Cartesian(rotatingMotion.FirstDerivative.X, rotatingMotion.FirstDerivative.Y, updatedVz));
    var updatedInertial = fromMoonRotatingAtEpoch.Transform(updatedRotating, 1);
    configuration.ModifyMotion(motionID, updatedInertial);
};

SegmentPropagatorVariable velocityZVariable = initialStateSegment.CreateVariable(10.0, // Maximum step, m/s.
                                                                                 0.1, // Perturbation, m/s.
                                                                                 setterVz);

// Constraints are Earth departure altitude of 241 km and inclination of 28.5 degrees.
ScalarPointElement backwardPerigeeRadius = new ScalarPointElement(propagationPointAroundEarth.IntegrationPoint,
                                                                  CartesianElement.Magnitude,
                                                                  earth.InertialFrame);
ScalarAtEndOfNumericalSegmentConstraint backwardPerigeeRadiusConstraint = new ScalarAtEndOfNumericalSegmentConstraint(backwardPerigeeRadius,
                                                                                                                      propagateBackwardToPerigee,
                                                                                                                      earthPhysicalRadius + 241000.0, // Desired value, meters.
                                                                                                                      0.001, // Tolerance, meters.
                                                                                                                      EndsOfSegment.Final);

ScalarModifiedKeplerianElement backwardInclination = new ScalarModifiedKeplerianElement(earthGravitationalParameter,
                                                                                        propagationPointAroundEarth.IntegrationPoint,
                                                                                        KeplerianElement.Inclination,
                                                                                        earth.InertialFrame);
ScalarAtEndOfNumericalSegmentConstraint backwardInclinationConstraint = new ScalarAtEndOfNumericalSegmentConstraint(backwardInclination,
                                                                                                                    propagateBackwardToPerigee,
                                                                                                                    Trig.DegreesToRadians(28.5), // Desired value, radians.
                                                                                                                    Trig.DegreesToRadians(0.1), // Tolerance, radians.
                                                                                                                    EndsOfSegment.Final);

var correctorBackwardTrajectory = new TargetedSegmentListDifferentialCorrector();
correctorBackwardTrajectory.Name = "Correct_Backward_Trajectory";
correctorBackwardTrajectory.MaximumIterations = 75;
correctorBackwardTrajectory.OnlyStoreFinalResults = false;

correctorBackwardTrajectory.Constraints.Add(backwardPerigeeRadiusConstraint);
correctorBackwardTrajectory.Constraints.Add(backwardInclinationConstraint);

correctorBackwardTrajectory.Variables.Add(velocityYVariable);
correctorBackwardTrajectory.Variables.Add(velocityZVariable);

// Add the differential corrector to the targeted segment list.
backwardSegmentList.Operators.Add(correctorBackwardTrajectory);

backwardSegmentList.OperatorAction = TargetedSegmentListOperatorBehavior.RunActiveOperators;

// Run backwardSegmentList.
SegmentPropagator propagator = backwardSegmentList.GetSegmentPropagator(new EvaluatorGroup());
TargetedSegmentListResults listResults = (TargetedSegmentListResults)propagator.Propagate();

DateMotionCollection<Cartesian> barycentricRotatingEphemeris = listResults.GetDateMotionCollectionOfOverallTrajectory(motionID, barycentricRotatingFrame);

// Barycentric or moon centered inertial ephemeris could also be created here if desired.
DateMotionCollection<Cartesian> earthCenteredInertialEphemeris = listResults.GetDateMotionCollectionOfOverallTrajectory(motionID, earth.InertialFrame);

// Compute translunar injection and parking orbit parameters.
Motion<Cartesian> translunarInjectionState = earthCenteredInertialEphemeris.Motions[0];
ModifiedKeplerianElements translunarKeplerianElements = new ModifiedKeplerianElements(translunarInjectionState, earthGravitationalParameter);

double translunarInjectionVelocity = translunarInjectionState.FirstDerivative.Magnitude;
double parkingOrbitVelocity = Math.Sqrt(earthGravitationalParameter / translunarInjectionState.Value.Magnitude);
double translunarInjectionDeltaV = translunarInjectionVelocity - parkingOrbitVelocity;

Cartesian parkingOrbitVelocityCartesian = translunarInjectionState.FirstDerivative * (parkingOrbitVelocity / translunarInjectionVelocity);

Motion<Cartesian> parkingOrbitState = new Motion<Cartesian>(translunarInjectionState.Value, parkingOrbitVelocityCartesian);
ModifiedKeplerianElements parkingOrbitElements = new ModifiedKeplerianElements(parkingOrbitState, earthGravitationalParameter);

Free Return Forward Propagation Setup

This code is common to both the cislunar and the circumlunar examples. It sets up the force models, propagators, segments, and differential correctors that are used to perform the forward targeting of the Earth re-entry conditions that occur after the Lunar flyby. The parking orbit guess determined in the previous section is varied to target those Earth re-entry conditions.

Copy

// Use Earth parking orbit right ascension of the ascending node, true anomaly, and translunar injection delta-v as variables
// to target Earth re-entry after the lunar flyby.
NumericalInitialStateSegment forwardInitialStateSegment = new NumericalInitialStateSegment();
forwardInitialStateSegment.Name = "Forward_Initial_State_Segment";

// Use cloned version of the Earth centered backward propagator that has modified initial state and epoch.
var forwardEarthCenteredPropagator = (NumericalPropagatorDefinition)earthCenteredNumericalPropagator.Clone(new CopyContext());
forwardEarthCenteredPropagator.Epoch = earthCenteredInertialEphemeris.Dates[0];
var forwardEarthCenteredPoint = (PropagationNewtonianPoint)((DefinitionalObjectCollection<PropagationStateElement>)forwardEarthCenteredPropagator.IntegrationElements)[0];

forwardEarthCenteredPoint.InitialPosition = parkingOrbitState.Value;
forwardEarthCenteredPoint.InitialVelocity = parkingOrbitState.FirstDerivative;

forwardInitialStateSegment.PropagatorDefinition = forwardEarthCenteredPropagator;

// The translunar injection from the parking orbit.
ImpulsiveManeuverSegment transLunarInjectionSegment = new ImpulsiveManeuverSegment();
transLunarInjectionSegment.Name = "Trans-Lunar_Injection";

// The maneuver.
ImpulsiveManeuverInformation transLunarBurnDetails = new ImpulsiveManeuverInformation(motionID, new Cartesian(translunarInjectionDeltaV, 0.0, 0.0));
transLunarBurnDetails.Orientation = new AxesVelocityOrbitNormal(transLunarBurnDetails.PropagationPoint, earth);
transLunarInjectionSegment.Maneuver = transLunarBurnDetails;

// Right ascension variable. Perturbation values for each of these variables
// may need to be modified using trial and error if convergence fails.
SetVariableCallback<NumericalPropagatorState> setterRightAscension = (variable, configuration) =>
{
    Motion<Cartesian> currentMotion = configuration.GetMotion<Cartesian>(motionID);

    ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter);
    ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(oldElements.RadiusOfPeriapsis,
                                                                          oldElements.InverseSemimajorAxis,
                                                                          oldElements.Inclination,
                                                                          oldElements.ArgumentOfPeriapsis,
                                                                          oldElements.RightAscensionOfAscendingNode + variable,
                                                                          oldElements.TrueAnomaly,
                                                                          earthGravitationalParameter);
    configuration.ModifyMotion(motionID, newElements.ToCartesian());
};

SegmentPropagatorVariable rightAscensionVariable = forwardInitialStateSegment.CreateVariable(Trig.DegreesToRadians(0.1), // Maximum step, radians.
                                                                                             Trig.DegreesToRadians(0.001), // Perturbation, radians.
                                                                                             setterRightAscension);
rightAscensionVariable.Name = "InitialState_RightAscension";

// True anomaly variable. (This is really argument of latitude since the orbit is circular and inclined.) 
SetVariableCallback<NumericalPropagatorState> setterTrueAnomaly = (variable, configuration) =>
{
    Motion<Cartesian> currentMotion = configuration.GetMotion<Cartesian>(motionID);
    ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter);
    ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(oldElements.RadiusOfPeriapsis,
                                                                          oldElements.InverseSemimajorAxis,
                                                                          oldElements.Inclination,
                                                                          oldElements.ArgumentOfPeriapsis,
                                                                          oldElements.RightAscensionOfAscendingNode,
                                                                          oldElements.TrueAnomaly + variable,
                                                                          earthGravitationalParameter);
    configuration.ModifyMotion(motionID, newElements.ToCartesian());
};

SegmentPropagatorVariable trueAnomalyVariable = forwardInitialStateSegment.CreateVariable(Trig.DegreesToRadians(0.1), // Maximum step, radians.
                                                                                          Trig.DegreesToRadians(0.001), // Perturbation, radians.
                                                                                          setterTrueAnomaly);
trueAnomalyVariable.Name = "InitialState_TrueAnomaly";

// Maneuver delta-v variable.
SegmentPropagatorVariable translunarInjectionDeltaVxVariable = transLunarInjectionSegment.CreateVariable(100.0, // Maximum step, m/s.
                                                                                                         1.0, // Perturbation, m/s.
                                                                                                         (variable, configuration) => { configuration.Maneuver.X += variable; });
translunarInjectionDeltaVxVariable.Name = "TranslunarInjection_ThrustVector_X";

// Make a segment that propagates forward to edge of Lunar SOI.
NumericalPropagatorSegment propagateForwardToSoi1 = new NumericalPropagatorSegment();
propagateForwardToSoi1.Name = "Propagate_Forward_ToSOI_BeforeFlyby";
propagateForwardToSoi1.PropagatorDefinition = forwardEarthCenteredPropagator;
propagateForwardToSoi1.PropagationDirection = IntegrationSense.Increasing;
propagateForwardToSoi1.StopOnMaximumDurationBehavior = MaximumDurationBehavior.Throw;

// Stopping condition can be reused here.
propagateForwardToSoi1.StoppingConditions.Add(moonSoiStoppingCondition);

// Make a segment that propagates forward until perilune (closest approach distance to Moon)
// using Moon-centered propagator.
var propagateForwardToPerilune = new NumericalPropagatorSegment();
propagateForwardToPerilune.Name = "Propagate_Forward_ToPerilune";

// No need to change initial conditions here because the state and time will be overriden.
propagateForwardToPerilune.PropagatorDefinition = nearMoonNumericalPropagator;
propagateForwardToPerilune.PropagationDirection = IntegrationSense.Increasing;
propagateForwardToPerilune.StopOnMaximumDurationBehavior = MaximumDurationBehavior.Throw;

ScalarStoppingCondition periluneStoppingCondition = new ScalarStoppingCondition(new ScalarPointElement(propagationPointAroundMoon.IntegrationPoint,
                                                                                                       CartesianElement.Magnitude,
                                                                                                       moon.InertialFrame),
                                                                                0.0001, // Value tolerance, meters.
                                                                                StopType.LocalMinimum);
periluneStoppingCondition.Name = "Moon_Perilune";

propagateForwardToPerilune.StoppingConditions.Add(periluneStoppingCondition);

// Make a segment that propagates forward to edge of Lunar SOI.
NumericalPropagatorSegment propagateForwardToSoi2 = new NumericalPropagatorSegment();
propagateForwardToSoi2.Name = "Propagate_Forward_ToSOI_AfterFlyby";
propagateForwardToSoi2.PropagatorDefinition = nearMoonNumericalPropagator;
propagateForwardToSoi2.PropagationDirection = IntegrationSense.Increasing;
propagateForwardToSoi2.StopOnMaximumDurationBehavior = MaximumDurationBehavior.Throw;

// Stopping condition can be reused here.
propagateForwardToSoi2.StoppingConditions.Add(moonSoiStoppingCondition);

// Make a segment that propagates forward two days.
NumericalPropagatorSegment propagateForwardTwoDays = new NumericalPropagatorSegment();
propagateForwardTwoDays.Name = "Propagate_Forward_Two_Days";
propagateForwardTwoDays.PropagatorDefinition = forwardEarthCenteredPropagator;
propagateForwardTwoDays.StoppingConditions.Add(new DurationStoppingCondition(Duration.FromDays(2.0)));

// Make a segment that propagates forward to perigee or Earth altitude of 122.5 km.
NumericalPropagatorSegment propagateForwardToPerigee = new NumericalPropagatorSegment();
propagateForwardToPerigee.Name = "Propagate_Backward_ToPerigee";
propagateForwardToPerigee.PropagatorDefinition = forwardEarthCenteredPropagator;
propagateForwardToPerigee.PropagationDirection = IntegrationSense.Increasing;
propagateForwardToPerigee.StopOnMaximumDurationBehavior = MaximumDurationBehavior.Throw;

ScalarStoppingCondition perigeeStoppingCondition2 = new ScalarStoppingCondition(new ScalarPointElement(forwardEarthCenteredPoint.IntegrationPoint,
                                                                                                       CartesianElement.Magnitude,
                                                                                                       earth.InertialFrame),
                                                                                0.0001, // Value tolerance, meters.
                                                                                StopType.LocalMinimum);
perigeeStoppingCondition2.Name = "Earth_Perigee_After_Return";

ScalarStoppingCondition altitudeStoppingCondition =  new ScalarStoppingCondition(new ScalarPointElement(forwardEarthCenteredPoint.IntegrationPoint,
                                                                                                        CartesianElement.Magnitude,
                                                                                                        earth.InertialFrame),
                                                                                 earthPhysicalRadius + 122500.0, // Threshold, meters.
                                                                                 0.0001, // Value tolerance, meters.
                                                                                 StopType.AnyThreshold);
altitudeStoppingCondition.Name = "Earth_Entry_Altitude";

propagateForwardToPerigee.StoppingConditions.Add(perigeeStoppingCondition2);
propagateForwardToPerigee.StoppingConditions.Add(altitudeStoppingCondition);

// A re-entry altitude of 120 - 125 km with a flight path angle of between -6 and -7 degrees
// are the target conditions for manned re-entry used by Mark Jesick, "AAS 11-452, Abort Options for Human
// Missions to Earth-Moon Halo orbits" in Proceedings of the 2011 AAS/AIAA Astrodynamics Conference,
// Girdwood, AK, August 2011.

// Choosing the midpoint of those ranges gives altitude constraint 125.5 km and flight path angle
// constraint of -6.5 degree for re-entry conditions.
ScalarPointElement radiusScalar = new ScalarPointElement(forwardEarthCenteredPoint.IntegrationPoint, CartesianElement.Magnitude, earth.InertialFrame);
ScalarAtEndOfNumericalSegmentConstraint altitudeConstraint = new ScalarAtEndOfNumericalSegmentConstraint(radiusScalar,
                                                                                                         propagateForwardToPerigee,
                                                                                                         earthPhysicalRadius + 122500.0, // Desired value, meters.
                                                                                                         100.0, // Tolerance, meters.
                                                                                                         EndsOfSegment.Final);

// Calculating Earth's J2 is needed to get a radial velocity scalar from
// the Kozai-Izsak mean elements.
double earthJ2 = earthGravityModel.ZonalCoefficients[2];
earthJ2 = earthGravityModel.Normalized ? earthJ2 * Math.Sqrt(5.0) : earthJ2;

ScalarKozaiIzsakMeanElement radialVelocityScalar = new ScalarKozaiIzsakMeanElement(earthGravitationalParameter,
                                                                                   forwardEarthCenteredPoint.IntegrationPoint,
                                                                                   KozaiIzsakMeanElement.RadialVelocity,
                                                                                   earth.InertialFrame,
                                                                                   earthJ2,
                                                                                   earthPhysicalRadius);

ScalarVectorElement velocityMagnitudeScalar = new ScalarVectorElement(forwardEarthCenteredPoint.IntegrationPoint.CreateVectorVelocity(earth.InertialFrame),
                                                                      CartesianElement.Magnitude,
                                                                      earth.InertialFrame.Axes);

Scalar sineFlightPathAngleScalar = radialVelocityScalar / velocityMagnitudeScalar;

double desiredSineFlightPathAngle = Math.Sin(Trig.DegreesToRadians(-6.5));

ScalarAtEndOfNumericalSegmentConstraint flightPathAngleConstraint = new ScalarAtEndOfNumericalSegmentConstraint(sineFlightPathAngleScalar,
                                                                                                                propagateForwardToPerigee,
                                                                                                                desiredSineFlightPathAngle, // Desired value, radians.
                                                                                                                0.001, // Tolerance, radians.
                                                                                                                EndsOfSegment.Final);

// Set up segments to use forward propagation.
TargetedSegmentList forwardSegmentList = new TargetedSegmentList();
forwardSegmentList.Name = "Forward_Segment_List";

forwardSegmentList.Segments.Add(forwardInitialStateSegment);
forwardSegmentList.Segments.Add(transLunarInjectionSegment);
forwardSegmentList.Segments.Add(propagateForwardToSoi1);
forwardSegmentList.Segments.Add(propagateForwardToPerilune);
forwardSegmentList.Segments.Add(propagateForwardToSoi2);
forwardSegmentList.Segments.Add(propagateForwardTwoDays);
forwardSegmentList.Segments.Add(propagateForwardToPerigee);

TargetedSegmentListDifferentialCorrector altitudeCorrectorForwardTrajectory = new TargetedSegmentListDifferentialCorrector();

altitudeCorrectorForwardTrajectory.Name = "Correct_Forward_Trajectory_Altitude";
altitudeCorrectorForwardTrajectory.Variables.Add(rightAscensionVariable);
altitudeCorrectorForwardTrajectory.Variables.Add(trueAnomalyVariable);
altitudeCorrectorForwardTrajectory.Variables.Add(translunarInjectionDeltaVxVariable);

altitudeCorrectorForwardTrajectory.Constraints.Add(altitudeConstraint);

altitudeCorrectorForwardTrajectory.MaximumIterations = 25;
altitudeCorrectorForwardTrajectory.OnlyStoreFinalResults = false;

forwardSegmentList.Operators.Add(altitudeCorrectorForwardTrajectory);

TargetedSegmentListDifferentialCorrector flightPathCorrectorForwardTrajectory = new TargetedSegmentListDifferentialCorrector();

flightPathCorrectorForwardTrajectory.Name = "Correct_Forward_Trajectory_FlightPath";
flightPathCorrectorForwardTrajectory.Variables.Add(rightAscensionVariable);
flightPathCorrectorForwardTrajectory.Variables.Add(trueAnomalyVariable);
flightPathCorrectorForwardTrajectory.Variables.Add(translunarInjectionDeltaVxVariable);

flightPathCorrectorForwardTrajectory.Constraints.Add(altitudeConstraint);
flightPathCorrectorForwardTrajectory.Constraints.Add(flightPathAngleConstraint);

flightPathCorrectorForwardTrajectory.MaximumIterations = 25;
flightPathCorrectorForwardTrajectory.OnlyStoreFinalResults = false;

forwardSegmentList.Operators.Add(flightPathCorrectorForwardTrajectory);

Results

This section demonstrates how to run the circumlunar and cislunar free-return forward targeters that were created above. After the runs are complete, the results are extracted and tabulated below.

Copy

// Run forwardSegmentList.
SegmentPropagator forwardPropagator = forwardSegmentList.GetSegmentPropagator(new EvaluatorGroup());
TargetedSegmentListResults forwardListResults = (TargetedSegmentListResults)forwardPropagator.Propagate();

DateMotionCollection<Cartesian> forwardEarthCenteredInertialEphemeris = forwardListResults.GetDateMotionCollectionOfOverallTrajectory(motionID, earth.InertialFrame);

DateMotionCollection<Cartesian> forwardEarthFixedEphemeris = forwardListResults.GetDateMotionCollectionOfOverallTrajectory(motionID, earth.FixedFrame);

// Note that it is also possible to create ephemeris in the barycentric rotating frame
// and the Moon-centered inertial frame.

// Get updated Earth parking orbit state and elements.
Motion<Cartesian> updatedParkingOrbitState = forwardEarthCenteredInertialEphemeris.Motions[0];
ModifiedKeplerianElements updatedParkingOrbitElements = new ModifiedKeplerianElements(updatedParkingOrbitState, earthGravitationalParameter);

JulianDate parkingOrbitEpoch = forwardEarthCenteredInertialEphemeris.Dates[0];
double parkingOrbitSemimajorAxis = 1.0 / updatedParkingOrbitElements.InverseSemimajorAxis;
double parkingOrbitEccentricity = updatedParkingOrbitElements.ComputeEccentricity();
double parkingOrbitInclination = Trig.RadiansToDegrees(updatedParkingOrbitElements.Inclination);
double parkingOrbitArgumentPerigee = Trig.RadiansToDegrees(updatedParkingOrbitElements.ArgumentOfPeriapsis);
double updatedRaan = Trig.RadiansToDegrees(updatedParkingOrbitElements.RightAscensionOfAscendingNode);
double updatedTrueAnomaly = Trig.RadiansToDegrees(updatedParkingOrbitElements.TrueAnomaly);

Motion<Cartesian> updatedTranslunarInjectionState = forwardEarthCenteredInertialEphemeris.Motions[1];
Cartesian updatedTranslunarInjectionDeltaV = updatedTranslunarInjectionState.FirstDerivative - updatedParkingOrbitState.FirstDerivative;
double updatedDeltaVMagnitude = updatedTranslunarInjectionDeltaV.Magnitude;

// Get lunar flyby state and elements.
SegmentResults propagateToPeriluneResults = forwardListResults.GetResultsOfSegment(propagateForwardToPerilune);

if (propagateToPeriluneResults != null)
{
    JulianDate periluneDate = propagateToPeriluneResults.FinalPropagatedState.CurrentDate;

    Motion<Cartesian> periluneState = propagateToPeriluneResults.FinalPropagatedState.GetMotion<Cartesian>(motionID);
    ModifiedKeplerianElements periluneElements = new ModifiedKeplerianElements(periluneState, moonGravitationalParameter);

    double lunarFlybyAltitude = periluneElements.RadiusOfPeriapsis - moonPhysicalRadius;
}

// Get Earth re-entry state and elements.
int forwardCount = forwardEarthCenteredInertialEphemeris.Count;
Motion<Cartesian> earthReentryState = forwardEarthCenteredInertialEphemeris.Motions[forwardCount - 1];
Cartesian earthFixedPosition = forwardEarthFixedEphemeris.Values[forwardCount - 1];
ModifiedKeplerianElements earthReentryElements = new ModifiedKeplerianElements(earthReentryState, earthGravitationalParameter);

JulianDate earthReentryDate = forwardEarthCenteredInertialEphemeris.Dates[forwardCount - 1];

Cartographic earthReentryCartographic = earth.Shape.CartesianToCartographic(earthFixedPosition);

double altitude = earthReentryCartographic.Height;
double latitude = Trig.RadiansToDegrees(earthReentryCartographic.Latitude);
double longitude = Trig.RadiansToDegrees(earthReentryCartographic.Longitude);

double eccentricity = earthReentryElements.ComputeEccentricity();
double trueAnomaly = earthReentryElements.TrueAnomaly;
double sineTrueAnomaly = Math.Sin(trueAnomaly);
double cosineTrueAnomaly = Math.Cos(trueAnomaly);

double sineFlightPathAngle = eccentricity * sineTrueAnomaly / Math.Sqrt(1.0 + 2.0 * eccentricity * cosineTrueAnomaly + eccentricity * eccentricity);

double flightPathAngle = Trig.RadiansToDegrees(Math.Asin(sineFlightPathAngle));

Results Table

Parameters	Circumlunar Results	Cislunar Results
Parking Orbit Epoch	12/29/2022 4:05:14 UTCG	12/19/2022 10:44:28 UTCG
Parking Orbit Inclination	28.5 degrees	28.5 degrees
Parking Orbit RAAN	8.59 degrees	8.37 degrees
Parking Orbit True Anomaly	90.00 degrees	-89.94 degrees
Translunar Injection Delta-v	3154.8 m/s	3148.0 m/s
Lunar Flyby Epoch	1/1/2023 00:01:20 UTCG	1/1/2023 00:13:01 UTCG
Lunar Flyby Altitude	145 km	373 km
Earth Re-entry Epoch	1/4/2023 00:29:47 UTCG	1/14/2023 4:34:32 UTCG
Earth Re-entry Altitude	122.5 km	122.9 km
Earth Re-entry Latitude	1.6 degrees South	8.1 degrees South
Earth Re-entry Longitude	81.32 degrees East	23.08 degrees East
Earth Re-entry Flight Path Angle	-6.46 degrees	-6.54 degrees

Discussion

Any high-fidelity mission design would require spacecraft-specific atmospheric drag and solar radiation pressure parameters. These parameters would be used as inputs for the solar radiation force and atmospheric drag models that are available in DME Component Libraries. Any propagator with atmospheric drag should only be active within 1000 km of Earth's surface. Solar radiation pressure should be added to all propagators because it is likely the largest un-modeled force in the above example code.

Spacecraft-specific dry mass, propellant mass, and specific impulse would also need to be provided to model the propellant mass expenditure of translunar injection maneuver.

Targeting and optimization are tricky for free returns because both the launch conditions and the re-entry conditions are often heavily constrained. This example assumed a 90-degree azimuth Cape Canaveral launch with an inclination of 28.5 degrees. The geographical areas where re-entry occurs may also be important if a manned vehicle can only land in the ocean or on specific land areas.