DME Component Libraries for Java 2024 r1

DME Component Libraries for Java 2024 r1

Programmer's Guide

Segment Propagation

Stopping Conditions

Targeted Segment List Operators

Multivariable Function Solvers

Lambert and Primer Vector Example

Trajectory Optimization Example

Lunar Free Return Code Sample

Click or drag to resize

Multivariable Function Solvers

Function solvers allow a system of multivariable equations to be driven to some set of values. DME Component Libraries includes one such solver type, the NewtonRaphsonMultivariableFunctionSolver. This solver will sequentially evaluate the objective function at the current variable vector, compute the derivative of the function, and iterate towards the solution. This process will be repeated until the difference between the computed function value and the constraint's target values are equal within tolerance. To assist with convergence the input variables have several settings such as, MaximumStep (get / set), VariableTolerance (get / set), and, Scaling (get / set). For more information on the available settings see the settings class.

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

Newton Raphson Function Solver

Suppose you want to solve the following pair of equations for x and y:

x² + y² = 25

x + 3.5 sin(x) = y

First it would help to get both equations to be implicit equations. The first equation already has the constant on one side, but the second should be rewritten as:

x - y + 3.5 sin(x) = 0

Eventually, we will need to define a SolvableMultivariableFunction, but first we start by building up the variables and constraints. The following code sample creates two variables, one for x and one for y:

Java

SolverVariableSettings xVariable = new SolverVariableSettings(0.3); // maximum step
// initial value is left at the default 0
xVariable.setName("x");

SolverVariableSettings yVariable = new SolverVariableSettings(0.3);
yVariable.setName("y");

Next, we create two constraints, each corresponding to a value computed by the function. These define the criteria that the function will be solved for. In this case, the constraints are the answer for the individual equations, with the desired value corresponding to the answer of each of the equations.

Java

SolverConstraintSettings firstEquation = new SolverConstraintSettings(25, // desired value
                                                                      0.001); // tolerance
firstEquation.setName("firstEquation");

SolverConstraintSettings secondEquation = new SolverConstraintSettings(0.0, 0.001);
secondEquation.setName("secondEquation");

These variable and constraint settings are then added to the function solver, as seen in the code sample below. The order in which the variables and constraint settings are added to the function solver must match the order that the variables are used and the constraints are computed in the evaluate method of the function.

Java

NewtonRaphsonMultivariableFunctionSolver solver = new NewtonRaphsonMultivariableFunctionSolver();
solver.getConstraints().add(firstEquation);
solver.getConstraints().add(secondEquation);
solver.getVariables().add(xVariable);
solver.getVariables().add(yVariable);

With the variables and the constraints created, we now define the function. The function must derive from the SolvableMultivariableFunction abstract base class, and implement the evaluate method. This method will be given an array of the current values of the variables, as computed by the NewtonRaphsonMultivariableFunctionSolver. Then, it must compute each equation, and assign those values to the results, in the same order as the constraints on the NewtonRaphsonMultivariableFunctionSolver. The following code sample shows the implementation of the evaluate method:

Java

@Override
public SolvableMultivariableFunctionResults evaluate(double[] variables, ITrackCalculationProgress progressTracker) {
    double[] equationAnswers = new double[2];
    double x = variables[0]; // x is always the first variable in this type
    double y = variables[1]; // y is always the second variable in this type

    // the order of the constraint values must match the order of the constraints in this type
    equationAnswers[0] = x * x + y * y;
    equationAnswers[1] = x - y + 3.5 * Math.sin(x);

    return new SolvableMultivariableFunctionResults(variables, equationAnswers);
}

public final void setPerturbationValues(double xPerturbation, double yPerturbation) {
    setPerturbationValues(new double[] {
        xPerturbation,
        yPerturbation
    });
}

The following code shows how to create and run the NewtonRaphsonMultivariableFunctionSolver:

Java

// SimpleFunction is our example class derived from SolvableMultivariableFunction
SimpleFunction function = new SimpleFunction();
function.setPerturbationValues(0.1, 0.1);
solver.setFunction(function);

solver.findSolution(25);
MultivariableFunctionSolverResults<MultivariableFunctionSolverIterationResults> results = solver.getLastRunsResults();

double[] variablesUsed = results.getFinalIteration().getFunctionResult().getVariablesUsed();
double foundX = variablesUsed[0];
double foundY = variablesUsed[1];
// 1.3744 4.8073

You can also check the variables directly to see that x = 1.37448 and y = 4.8073. However, if you graph the functions, you will see that there are multiple solutions:

Multivariable Function Plot

The initial guess will influence the results that the solver finds. Because of this fact, if you change your initial guess, it will find different answers:

Java

// recreate the variables with new initial values
xVariable = new SolverVariableSettings(0.3, 5.0, 0.001);
yVariable = new SolverVariableSettings(0.3, 5.0, 0.001);

function = new SimpleFunction();
function.setPerturbationValues(0.1, 0.1);

// recreate the solver
solver = new NewtonRaphsonMultivariableFunctionSolver();
solver.getVariables().add(xVariable);
solver.getVariables().add(yVariable);
solver.getConstraints().add(firstEquation);
solver.getConstraints().add(secondEquation);
solver.setFunction(function);
solver.findSolution(25);

results = solver.getLastRunsResults();
variablesUsed = results.getFinalIteration().getFunctionResult().getVariablesUsed();

foundX = variablesUsed[0];
foundY = variablesUsed[1];
// 4.8178 1.3372

One limitation of the Newton Raphson method is that it assumes a linear function. Consider if the slope of the function was near 0 when one of the variable is perturbed. The variable would be asked to take a large step that, if taken, could take it well beyond the parts of the function you are interested in. The MaximumStep (get / set) property on the variables prevents those large jumps.

Segmented Propagation and Function Solvers

The function you are solving can be very abstract. For example, you can set up a list of segments to compute a spacecraft's ephemeris, compute some value based on that ephemeris (such as the eccentricity and inclination at a point in time), and set up a variable that modifies the delta-v of some maneuver. In effect, you have the following equation, which can be run through the TargetedSegmentListDifferentialCorrector:

f(delta-v-x, delta-v-z) = [inclination, eccentricity]

Targeted Segment List Differential Corrector

TargetedSegmentListDifferentialCorrector requires that at least one SegmentPropagatorConstraint and one SegmentPropagatorVariable are configured on it, as shown below:

Java

// create the differential corrector
TargetedSegmentListDifferentialCorrector differentialCorrector = new TargetedSegmentListDifferentialCorrector();
differentialCorrector.setName("Corrector");

// add the differential corrector to the TargetedSegmentList
targetedSegmentList.getOperators().add(differentialCorrector);

// create a variable
DelegateBasedVariable<ImpulsiveManeuverSegmentConfiguration> velocityXVariable = impulsiveManeuverSegment.createVariable(
        200.0, // maximum step, meters/second 
        1.0,   // perturbation, meters/second
        SetVariableCallback.of((variable, configuration) -> {
            configuration.getManeuver().setX(configuration.getManeuver().getX() + variable);
        }));
velocityXVariable.setName("Velocity_X_Variable");

// add the variable to the differential corrector
differentialCorrector.getVariables().add(velocityXVariable);

// create a constraint
DelegateBasedConstraint periodConstraint = new DelegateBasedConstraint(
        DelegateBasedConstraintCallback.of(segmentResults -> {
            Motion1<Cartesian> objectsMotion = segmentResults.getStateForNextSegment().getMotion(SatelliteMotionIdentification);
            return new ModifiedKeplerianElements(
                    objectsMotion, 
                    WorldGeodeticSystem1984.GravitationalParameter).computePeriod();
        }),
        impulsiveManeuverSegment, 
        TimeConstants.SecondsPerDay, // desired value, seconds 
        0.01); // tolerance, seconds

periodConstraint.setName("Period_Constraint");

// add the constraint to the differential corrector
differentialCorrector.getConstraints().add(periodConstraint);

Segment Propagation Constraints

There are several constraints included in the Segment Propagation Library. The first is DelegateBasedConstraint, which allows you to specify a delegate that returns some computed value from the ephemeris of the segments propagated. For example, if you wanted to constrain an orbit such that its period is 24 hours, you could do so as shown below:

Java

DelegateBasedConstraint constraint = new DelegateBasedConstraint(
        DelegateBasedConstraintCallback.of(segmentResults -> {
            Motion1<Cartesian> objectsMotion = segmentResults.getStateForNextSegment().getMotion(satelliteID);
            return new ModifiedKeplerianElements(objectsMotion, WorldGeodeticSystem1984.GravitationalParameter).computePeriod();
        }), 
        segment, 
        TimeConstants.SecondsPerDay, // a desired value of 1 day in seconds
        0.01); // tolerance, seconds
solver.getConstraints().add(constraint);

The other constraints generally consist of evaluating a Scalar of some sort. For example, the ScalarAtEndOfSegmentConstraint will evaluate a Scalar at the end of a segment in the TargetedSegmentList. For example, if you want your segment final state to be at a specific altitude, you could specify that as shown below:

Java

ScalarAtEndOfSegmentConstraint constraint =
        new ScalarAtEndOfSegmentConstraint(segment, 
                                           500000.0, // desired value of 500 km
                                           1.0); // 1 meter tolerance
ScalarCartographicElement altitudeScalar = 
        new ScalarCartographicElement(earth, 
                                      new ParameterizedOnStatePoint(constraint.getParameter(), 
                                                                    earth.getFixedFrame(), 
                                                                    SatelliteMotionIdentification),
                                      CartographicElement.HEIGHT);
constraint.setScalar(altitudeScalar);

solver.getConstraints().add(constraint);

Segment Propagation Variables

Many segments have values in them that can be used as variables. Although each segment has different properties that can be used as variables, the process used to create variables is similar across the different types of segments.

In the first example, where we showed solving a simple system of multivariable equations, the evaluate method of the function was created such that the variables were the familiar X and Y used in mathematics literature. In more complicated functions present within the Segment Propagation Library, the variables can be parameterized primitives or elements within the propagation state. These variables may be trivially altered, require auxiliary calculation, or even require that they are applied sequentially. To handle the more complex behavior involved in modeling changes to the propagation elements and segment primitives, we include a class that handles most of the heavy lifting- TargetedSegmentListFunction.

To illustrate how the complex behavior of variables could be a problem: consider an example where the initial state of propagation has two variables being solved for, the inclination and the right ascension of the ascending node. As we will explore shortly, the inclination may also need to change the right ascension (if the function solver tries to push the inclination to be less than 0 or more than 180 degrees, other orbital elements must flip otherwise the actual position of the satellite will be on the other side of the central body). This complex behavior requires the use of DelegateBasedVariable<T>. See the example under the appropriate section for the details.

Now let's dive a little further into the specifics of the underlying mechanics in SegmentPropagatorVariables. For the delegate type variable one typically modifies a value directly present in the state or performs a complex calculation which then modifies the state in some manner. If you look at the below example variable in the impulsive maneuver segment, you'll notice that when we call our variable setter we add the input value to the configuration's maneuver element. Each time the TargetedSegmentListFunction is evaluated in the solver it passes a clean configuration containing any initial state either passed-in or set when constructing the propagation sequence. Thus the passed in configuration functions as an initial value for the variable that the delegate alters. The solver itself maintains a record of the variable's current value with no knowledge of the initial value supplied by the configuration.

This behavior is similar for the other types of SegmentPropagatorVariables. For each of these types the initial value set on the definitional object itself (the variable) is used in propagation/function evaluation but is not present in iteration output. What this means is that when querying the values of the variables the iteration results will all be relative to the initial value on the definitional object or found within the clean configuration. If inspecting state elements directly the initial value will be incorporated and thus be the true value used in the computations and propagation.

Delegate Based Variable

DelegateBasedVariable<T> works with a setter delegate. The setter will take in the segments configuration and the actual value of the variable, and this delegate must modify the appropriate value in the configuration.

InitialStateSegment<T> can have any value in its state edited. If, for example, you don't know what inclination your satellite should start at, you can create such variables by calling the createVariable method on the segment:

Java

DelegateBasedVariable<InitialStateSegmentConfiguration> inclinationVariable = initialStateSegment.createVariable(
        0.2,  // maximum step, radians
        0.01, // perturbation, radians
        SetVariableCallback.of((variable, configuration) -> {
            Motion1<Cartesian> currentMotion = configuration.getMotion(SatelliteMotionIdentification);
            ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, gravitationalParameter);
            double rightAscensionCorrection = 0.0;
            double argumentOfPeriapsisCorrection = 0.0;
            double inclination = oldElements.getInclination() + variable;
            if (inclination > Math.PI) {
                // this checks are in here because it is possible that the differential corrector could
                // push the inclination to be greater than PI, and if it does we need to flip the
                // inclination and RightAscension.
                inclination = inclination - Math.PI;
                rightAscensionCorrection = -Math.PI;
                argumentOfPeriapsisCorrection = -Math.PI;
            } else if (inclination < 0) {
                // This is correct. If the inclination becomes less than 0, to avoid discontinuities
                // in position and velocity, the right ascension and AOP need to flip too.
                inclination = inclination + Math.PI;
                rightAscensionCorrection = Math.PI;
                argumentOfPeriapsisCorrection = Math.PI;
            }

            double newArgumentOfPeriapsis = Trig.zeroToTwoPi(oldElements.getArgumentOfPeriapsis() + argumentOfPeriapsisCorrection);
            double newRightAscension = Trig.zeroToTwoPi(oldElements.getRightAscensionOfAscendingNode() + rightAscensionCorrection);

            ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(
                    oldElements.getRadiusOfPeriapsis(),
                    oldElements.getInverseSemimajorAxis(), 
                    inclination, 
                    newArgumentOfPeriapsis, 
                    newRightAscension,
                    oldElements.getTrueAnomaly(), 
                    oldElements.getGravitationalParameter());

            configuration.modifyMotion(SatelliteMotionIdentification, newElements.toCartesian());
        }));

When making the delegates for variables, you must remember that the Newton Raphson solver does not know anything about the physics of the problem. It doesn't know that the inclination can only be between 0 and 180 degrees. When it computes how much the variable must change in order to meet its constraints, it may try to push the inclination beyond its bounds. It is up to the variable to handle that properly, hence the extra checks in this setter delegate.

Notice how by modifying the inclination you may also need to modify other values too. This is perfectly acceptable as long as conceptually one variable is being modified.

However, if you want to solve for both inclination and for the right ascension of the ascending node, you would need a second variable. Do not be tempted to try to apply two separate variables in one delegate. Do not do this:

Java

// DO NOT DO THIS

DelegateBasedVariable<InitialStateSegmentConfiguration> incorrectVariable = initialStateSegment.createVariable(
        0.2,  // maximum step, radians
        0.01, // perturbation, radians
        SetVariableCallback.of((variable, configuration) -> {
            Motion1<Cartesian> currentMotion = configuration.getMotion(SatelliteMotionIdentification);
            ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, gravitationalParameter);
            ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(
                    oldElements.getRadiusOfPeriapsis(), 
                    oldElements.getInverseSemimajorAxis(),
                    oldElements.getInclination() + variable,  // this will not work properly
                    oldElements.getRightAscensionOfAscendingNode() + variable,  // this will not work properly
                    oldElements.getArgumentOfPeriapsis(),
                    oldElements.getTrueAnomaly(),
                    gravitationalParameter);
            configuration.modifyMotion(SatelliteMotionIdentification, newElements.toCartesian());
        }));

The function solver will try to solve it and it may even converge, but the results you get will not be what you expect.

Any of the values in an ImpulsiveManeuverInformation can be directly edited, except for the Orientation (get / set). Varying a delta-v is a common use case:

Java

DelegateBasedVariable<ImpulsiveManeuverSegmentConfiguration> velocityXVariable = firstBurn.createVariable(
        100.0, // maximum step, meters/second
        0.1,   // perturbation, meters/second
        SetVariableCallback.of((variable, configuration) -> {
            configuration.getManeuver().setX(configuration.getManeuver().getX() + variable);
        }));

Parameterized Variables

Throughout the library, there are many settings that you might want to vary in a function solver. For those, consider using ParameterizedScalarVariable, ParameterizedDoubleVariable, ParameterizedDurationVariable, or ParameterizedDateVariable. These variables all work the same way, by providing a parameterized placeholder for the value in question.

For example, any Scalar in any segment can be a variable by using the ParameterizedScalarVariable. This variable requires that you use its Value (get) property, in place of the Scalar, as shown below:

Java

// What drag coefficient will give my satellite an altitude of 250 km after 1 day of propagation?

// common items
EarthCentralBody earth = CentralBodiesFacet.getFromContext().getEarth();
double gravityParameter = WorldGeodeticSystem1984.GravitationalParameter;

SolarGeophysicalData solarData = new ConstantSolarGeophysicalData(150.0, 3.0);

// make the segment list
TargetedSegmentList targetedSegment = new TargetedSegmentList();

// the propagate segment
NumericalPropagatorSegment propagationSegment = new NumericalPropagatorSegment();

// configure the numerical propagator
NumericalPropagatorDefinition propagatorDefinition = new NumericalPropagatorDefinition();
propagationSegment.setPropagatorDefinition(propagatorDefinition);
propagatorDefinition.setEpoch(TimeConstants.J2000);
RungeKuttaFehlberg78Integrator integrator = new RungeKuttaFehlberg78Integrator();
integrator.setRelativeTolerance(Constants.Epsilon13);
integrator.setMinimumStepSize(1.0);
integrator.setMaximumStepSize(86400.0);
integrator.setInitialStepSize(300.0);
integrator.setStepSizeBehavior(KindOfStepSize.RELATIVE);
propagatorDefinition.setIntegrator(integrator);
propagationSegment.getStoppingConditions().add(new DurationStoppingCondition(Duration.fromDays(1.0)));

// create the numerical propagator point
PropagationNewtonianPoint satellitePoint = new PropagationNewtonianPoint();
propagatorDefinition.getIntegrationElements().add(satellitePoint);
satellitePoint.setIdentification(SatelliteMotionIdentification);
satellitePoint.setMass(Scalar.toScalar(500.0));

//Values of an orbit starting at 255 km
satellitePoint.setInitialPosition(new Cartesian(6633140.0, 0.0, 0.0));
satellitePoint.setInitialVelocity(new Cartesian(0.0, 6812.52, 3689.9));

// add gravity
satellitePoint.getAppliedForces().add(new TwoBodyGravity(satellitePoint.getIntegrationPoint(), 
                                                         earth, 
                                                         gravityParameter));

// make the differential corrector
TargetedSegmentListDifferentialCorrector differentialCorrector = 
    new TargetedSegmentListDifferentialCorrector();
differentialCorrector.setName("Differential_Corrector");
differentialCorrector.getSolver().setMultithreaded(false);
differentialCorrector.setMaximumIterations(25);

// the variable for the differential corrector
ParameterizedScalarVariable dragCoefficientVariable = 
        new ParameterizedScalarVariable(0.1, // initial value of the drag coefficient
                                        0.3, // maximum step
                                        0.001, // perturbation
                                        propagationSegment);
dragCoefficientVariable.setName("Drag Coefficient");
dragCoefficientVariable.getSettings().setVariableTolerance(0.0001);

// the constraint
ScalarAtEndOfSegmentConstraint altitudeConstraint = 
        new ScalarAtEndOfSegmentConstraint(propagationSegment, 
                                           250000.0, // desired value, meters
                                           0.1); // tolerance, meters

ScalarCartographicElement altitudeScalar = 
        new ScalarCartographicElement(earth,
                                      new ParameterizedOnStatePoint(altitudeConstraint.getParameter(),
                                                                    satellitePoint.getIntegrationFrame(),
                                                                    satellitePoint.getIdentification()), 
                                      CartographicElement.HEIGHT);
altitudeConstraint.setScalar(altitudeScalar);

// add the variable and the constraint
differentialCorrector.getVariables().add(dragCoefficientVariable);
differentialCorrector.getConstraints().add(altitudeConstraint);

// add the differential corrector
targetedSegment.getOperators().add(differentialCorrector);

// create the drag force
ScalarAtmosphericDensity density = 
        new ScalarDensityMsis90(satellitePoint.getIntegrationPoint(), solarData);
AtmosphericDragForce dragForce = 
        new AtmosphericDragForce(density, 
                                 dragCoefficientVariable.getValue(), // here we use the parameterized scalar
                                 Scalar.toScalar(50.0)); // surface area, meters squared

// add the force
satellitePoint.getAppliedForces().add(dragForce);

// add the segment
targetedSegment.getSegments().add(propagationSegment);

// propagate
SegmentListPropagator propagator = targetedSegment.getSegmentListPropagator();

// since we know we were propagating a TargetedSegmentList, we can cast the results safely
TargetedSegmentListResults results = (TargetedSegmentListResults)propagator.propagateSegmentList();

TargetedSegmentListDifferentialCorrectorResults correctorResults = 
    (TargetedSegmentListDifferentialCorrectorResults)results.getOperatorResults().get(0);

SolvableMultivariableFunctionResults finalFunctionResults = 
    correctorResults.getFunctionSolverResults().getFinalIteration().getFunctionResult();

double finalAltitude = finalFunctionResults.getConstraintValues()[0];
double dragCoefficient = finalFunctionResults.getVariablesUsed()[0];
// final altitude will be 250000.0 meters + or - 0.1 meters
// the drag coefficient will be 0.03803

In fact, the DesiredValue (get / set) of a constraint in one differential corrector can be varied by an outer differential corrector using a ParameterizedDoubleVariable, as shown below:

Java

TargetedSegmentList outerTargetedSegmentList = new TargetedSegmentList();
outerTargetedSegmentList.setName("Outer_Target_List");

TargetedSegmentListDifferentialCorrector solverFor24HourPeriod = new TargetedSegmentListDifferentialCorrector();
solverFor24HourPeriod.setName("Solving_For_24_Hour_Period");

// 24 hour period constraint
ScalarAtEndOfSegmentConstraint periodConstraint = new ScalarAtEndOfSegmentConstraint(
        dv2, 
        TimeConstants.SecondsPerDay, // desired value, seconds
        0.0001);                     // tolerance, seconds
ScalarModifiedKeplerianElement periodScalar = new ScalarModifiedKeplerianElement(
        gravitationalParameter, 
        new ParameterizedOnStatePoint(
                periodConstraint.getParameter(), 
                frame, 
                motionID),
        KeplerianElement.PERIOD, 
        frame);
periodConstraint.setScalar(periodScalar);
periodConstraint.setName("Orbital_Period");

// vary the desired value of the nested Radius of Apoapsis
ParameterizedDoubleVariable nestedDesiredRadiusOfApoapsisVariable = new ParameterizedDoubleVariable(
        42000000.0, // initial value, meters
        5000000.0,  // maximum step, meters 
        500000.0,   // perturbation, meters
        transferOrbit);
nestedDesiredRadiusOfApoapsisVariable.setName("Desired_Radius_Of_Apoapsis");
radiusOfApoapsisPostDv1Constraint.setDesiredValue(nestedDesiredRadiusOfApoapsisVariable.getValue());

// add the variable and constraint
solverFor24HourPeriod.getConstraints().add(periodConstraint);
solverFor24HourPeriod.getVariables().add(nestedDesiredRadiusOfApoapsisVariable);

// add the differentialCorrector to the segment
outerTargetedSegmentList.getOperators().add(solverFor24HourPeriod);

Below is a complete example using a nested TargetedSegmentList, demonstrating many of the features described above:

Java

    String motionID = "Satellite";
    String fuelMassName = "Fuel_Mass";
    String dryMassName = "Dry_Mass";
    EarthCentralBody earth = CentralBodiesFacet.getFromContext().getEarth();
    ReferenceFrame frame = earth.getInertialFrame();
    double gravitationalParameter = WorldGeodeticSystem1984.GravitationalParameter;
    JulianDate epoch = TimeConstants.J2000.toTimeStandard(TimeStandard.getCoordinatedUniversalTime());
    Motion1<Cartesian> initialConditions = new Motion1<Cartesian>(
        new Cartesian(8000000.0, 0.0, 0.0), 
        new Cartesian(1500.0, 7500.0, 0.0));



    NumericalPropagatorDefinition numericalPropagatorDefinition = new NumericalPropagatorDefinition();
    numericalPropagatorDefinition.setEpoch(epoch);

    PropagationNewtonianPoint integrationPoint = new PropagationNewtonianPoint(
        motionID, 
        frame, 
        initialConditions.getValue(), 
        initialConditions.getFirstDerivative());

    // simple gravity
    TwoBodyGravity gravity = new TwoBodyGravity(
        integrationPoint.getIntegrationPoint(), 
        earth, 
        gravitationalParameter);
    integrationPoint.getAppliedForces().add(gravity);

    // fuel mass
    PropagationScalar fuelMass = new PropagationScalar(900.0);
    fuelMass.setIdentification(fuelMassName);
    fuelMass.setScalarDerivative(new ScalarFixed(0.0));

    // dry mass
    AuxiliaryStateScalar dryMass = new AuxiliaryStateScalar(Scalar.toScalar(100.0));
    dryMass.setIdentification(dryMassName);
    integrationPoint.setMass(Scalar.add(fuelMass.getIntegrationValue(), dryMass.getAuxiliaryScalar()));

    // configure the propagator
    RungeKuttaFehlberg78Integrator integrator = new RungeKuttaFehlberg78Integrator();
    integrator.setRelativeTolerance(Constants.Epsilon13);
    integrator.setMinimumStepSize(1.0);
    integrator.setMaximumStepSize(86400.0);
    integrator.setInitialStepSize(300.0);
    integrator.setStepSizeBehavior(KindOfStepSize.RELATIVE);

    numericalPropagatorDefinition.setIntegrator(integrator);
    numericalPropagatorDefinition.getIntegrationElements().add(integrationPoint);
    numericalPropagatorDefinition.getIntegrationElements().add(fuelMass);
    numericalPropagatorDefinition.getAuxiliaryElements().add(dryMass);



    NumericalInitialStateSegment initial = new NumericalInitialStateSegment();
    initial.setName("Initial_State_Segment");
    initial.setPropagatorDefinition(numericalPropagatorDefinition);



    NumericalPropagatorSegment propagateToPerigeeSegment = new NumericalPropagatorSegment();
    propagateToPerigeeSegment.setName("Propagate_To_Perigee");
    propagateToPerigeeSegment.setPropagatorDefinition(numericalPropagatorDefinition);

    // perigee stopping condition

    ScalarStoppingCondition perStop = new ScalarStoppingCondition(
            new ScalarModifiedKeplerianElement(
                    gravitationalParameter,
                    integrationPoint.getIntegrationPoint(),
                    KeplerianElement.TRUE_ANOMALY,
                    frame),
            0.0,          // threshold in radians
            0.0000000001, // value tolerance, radians
            StopType.ANY_THRESHOLD);
    perStop.setAngularSetting(CircularRange.NEGATIVE_PI_TO_PI); // to avoid the branch cut
    perStop.setName("Periapsis");

    propagateToPerigeeSegment.getStoppingConditions().add(perStop);



    TargetedSegmentList innerTargetSegment = new TargetedSegmentList();
    innerTargetSegment.setName("Inner_Target_List");


    // maneuver info

    ImpulsiveManeuverSegment dv1 = new ImpulsiveManeuverSegment();

    dv1.setName("First_Maneuver");

    // maneuver info
    ImpulsiveManeuverInformation dv1Info = new ImpulsiveManeuverInformation(
            motionID, 
            new Cartesian(200.0, 0.0, 0.0), // first guess
            fuelMassName, 
            dryMassName, 
            Scalar.toScalar(4500.0), // exhaust velocity, meters/second
            InvalidFuelStateBehavior.THROW_EXCEPTION);

    dv1Info.setOrientation(new AxesVelocityOrbitNormal(
            dv1Info.getPropagationPoint(), 
            earth));

    dv1.setManeuver(dv1Info);



    NumericalPropagatorSegment transferOrbit = new NumericalPropagatorSegment();
    transferOrbit.setName("Transfer_Orbit");
    transferOrbit.setPropagatorDefinition(numericalPropagatorDefinition);

    // apogee stopping condition
    ScalarStoppingCondition apoapsisStoppingCondition = new ScalarStoppingCondition(
            new ScalarModifiedKeplerianElement(
                    gravitationalParameter, 
                    integrationPoint.getIntegrationPoint(),
                    KeplerianElement.TRUE_ANOMALY, 
                    frame), 
            Math.PI,    // threshold, radians
            0.00000001, // value tolerance, radians
            StopType.ANY_THRESHOLD);
    apoapsisStoppingCondition.setAngularSetting(CircularRange.ZERO_TO_TWO_PI); // let the system know to avoid the branch cut
    apoapsisStoppingCondition.setName("Apoapsis_Stopping_Condition");
    transferOrbit.getStoppingConditions().add(apoapsisStoppingCondition);



    ImpulsiveManeuverSegment dv2 = new ImpulsiveManeuverSegment();

    dv2.setName("Second_Maneuver");
    ImpulsiveManeuverInformation dv2Info = new ImpulsiveManeuverInformation(
            motionID, 
            new Cartesian(500.0, 0.0, 0.0), 
            fuelMassName, 
            dryMassName, 
            Scalar.toScalar(4500.0),
            InvalidFuelStateBehavior.DO_FULL_DELTA_V);
    dv2.setManeuver(dv2Info);

    dv2Info.setOrientation(new AxesVelocityOrbitNormal(
            dv2Info.getPropagationPoint(), 
            earth));


    innerTargetSegment.getSegments().add(dv1);
    innerTargetSegment.getSegments().add(transferOrbit);
    innerTargetSegment.getSegments().add(dv2);


    TargetedSegmentListDifferentialCorrector hohmannTransferDifferentialCorrector = 
        new TargetedSegmentListDifferentialCorrector();
    hohmannTransferDifferentialCorrector.setName("Solve_For_Hohmann_Transfer");

    // variable 1, dv1's x velocity
    DelegateBasedVariable<ImpulsiveManeuverSegmentConfiguration> dv1VelocityXVariable = dv1.createVariable(
            200.0, // maximum step, meters/second
            0.1,   // perturbation, meters/second
            SetVariableCallback.of((variable, configuration) -> {
                configuration.getManeuver().setX(configuration.getManeuver().getX() + variable);
            }));
    dv1VelocityXVariable.setName("ImpulsiveManeuver_ThrustVector_X");

    // constraint 1, dv1's

    ScalarAtEndOfSegmentConstraint radiusOfApoapsisPostDv1Constraint = new ScalarAtEndOfSegmentConstraint(
            transferOrbit, 
            0, // desired value, meters, leaving as null since a outer variable will set this
            0.0001); // tolerance, meters
    ScalarModifiedKeplerianElement apoapsisScalar = new ScalarModifiedKeplerianElement(
            gravitationalParameter, 
            new ParameterizedOnStatePoint(
                    radiusOfApoapsisPostDv1Constraint.getParameter(),
                    frame, 
                    motionID),
            KeplerianElement.RADIUS_OF_APOAPSIS, 
            frame);

    radiusOfApoapsisPostDv1Constraint.setScalar(apoapsisScalar);
    radiusOfApoapsisPostDv1Constraint.setName("Radius_Of_Apoapsis");

    // variable 2, dv2's x velocity
    DelegateBasedVariable<ImpulsiveManeuverSegmentConfiguration> dv2VelocityXVariable = dv2.createVariable(
            200.0, // maximum step, meters/second
            0.1,   // perturbation, meters/second
            SetVariableCallback.of((variable, configuration) -> {
                configuration.getManeuver().setX(configuration.getManeuver().getX() + variable);
            }));
    dv2VelocityXVariable.setName("ImpulsiveManeuver_ThrustVector_X");

    // constraint 2, eccentricity constraint after the second maneuver
    ScalarAtEndOfSegmentConstraint eccentricityPostDv2Constraint = new ScalarAtEndOfSegmentConstraint(
            dv2, 
            0.1,       // desired value, unitless
            0.000001); // tolerance, unitless

    ScalarModifiedKeplerianElement eccentricityScalar = new ScalarModifiedKeplerianElement(
            gravitationalParameter, 
            new ParameterizedOnStatePoint(
                    eccentricityPostDv2Constraint.getParameter(),
                    frame, 
                    motionID), 
            KeplerianElement.ECCENTRICITY, 
            frame);
    eccentricityPostDv2Constraint.setScalar(eccentricityScalar);
    eccentricityPostDv2Constraint.setName("Eccentricity");

    // add the variables
    hohmannTransferDifferentialCorrector.getVariables().add(dv1VelocityXVariable);
    hohmannTransferDifferentialCorrector.getVariables().add(dv2VelocityXVariable);

    // add the constraints
    hohmannTransferDifferentialCorrector.getConstraints().add(radiusOfApoapsisPostDv1Constraint);
    hohmannTransferDifferentialCorrector.getConstraints().add(eccentricityPostDv2Constraint);

    // add the corrector to the targeted segment
    innerTargetSegment.getOperators().add(hohmannTransferDifferentialCorrector);



    TargetedSegmentList outerTargetedSegmentList = new TargetedSegmentList();
    outerTargetedSegmentList.setName("Outer_Target_List");

    TargetedSegmentListDifferentialCorrector solverFor24HourPeriod = new TargetedSegmentListDifferentialCorrector();
    solverFor24HourPeriod.setName("Solving_For_24_Hour_Period");

    // 24 hour period constraint
    ScalarAtEndOfSegmentConstraint periodConstraint = new ScalarAtEndOfSegmentConstraint(
            dv2, 
            TimeConstants.SecondsPerDay, // desired value, seconds
            0.0001);                     // tolerance, seconds
    ScalarModifiedKeplerianElement periodScalar = new ScalarModifiedKeplerianElement(
            gravitationalParameter, 
            new ParameterizedOnStatePoint(
                    periodConstraint.getParameter(), 
                    frame, 
                    motionID),
            KeplerianElement.PERIOD, 
            frame);
    periodConstraint.setScalar(periodScalar);
    periodConstraint.setName("Orbital_Period");

    // vary the desired value of the nested Radius of Apoapsis
    ParameterizedDoubleVariable nestedDesiredRadiusOfApoapsisVariable = new ParameterizedDoubleVariable(
            42000000.0, // initial value, meters
            5000000.0,  // maximum step, meters 
            500000.0,   // perturbation, meters
            transferOrbit);
    nestedDesiredRadiusOfApoapsisVariable.setName("Desired_Radius_Of_Apoapsis");
    radiusOfApoapsisPostDv1Constraint.setDesiredValue(nestedDesiredRadiusOfApoapsisVariable.getValue());

    // add the variable and constraint
    solverFor24HourPeriod.getConstraints().add(periodConstraint);
    solverFor24HourPeriod.getVariables().add(nestedDesiredRadiusOfApoapsisVariable);

    // add the differentialCorrector to the segment
    outerTargetedSegmentList.getOperators().add(solverFor24HourPeriod);


    NumericalPropagatorSegment finalPropagator = new NumericalPropagatorSegment();
    finalPropagator.setName("Final_Orbit");
    finalPropagator.setPropagatorDefinition(numericalPropagatorDefinition);
    finalPropagator.getStoppingConditions().add(new DurationStoppingCondition(Duration.fromDays(1.0)));
    finalPropagator.getStoppingConditions().get(0).setName("Duration");


    SegmentList overallList = new SegmentList();
    overallList.setName("Overall_List");
    overallList.getSegments().add(initial);
    overallList.getSegments().add(propagateToPerigeeSegment);

    // add the inner targeted segment list to the outer targeted segment list
    outerTargetedSegmentList.getSegments().add(innerTargetSegment);

    // and add the outer list
    overallList.getSegments().add(outerTargetedSegmentList);
    overallList.getSegments().add(finalPropagator);



    SegmentPropagator propagator = overallList.getSegmentPropagator(new EvaluatorGroup());
    SegmentListResults propagatorResults = (SegmentListResults) propagator.propagate();


    // double check that everything worked
    double finalPointRadius = 1.0
            / new ModifiedKeplerianElements(propagatorResults.getFinalSegmentResult().getStateForNextSegment().<Cartesian> getMotion(motionID), gravitationalParameter).getInverseSemimajorAxis();
    AssertHelper.assertEquals(Math.pow(Math.pow(TimeConstants.SecondsPerDay / Constants.TwoPi, 2.0) * WorldGeodeticSystem1984.GravitationalParameter, 1.0 / 3.0), finalPointRadius, 1.0);
}

@Test
public final void impulsiveManeuverVariableExample() {
    EarthCentralBody earth = getEarth();
    final String satelliteID = NumericalPropagatorHelperMethods.Id;
    final ImpulsiveManeuverSegment tempObj$2 = new ImpulsiveManeuverSegment();
    tempObj$2.setName("Trans-Lunar Injection");
    tempObj$2.setManeuver(new ImpulsiveManeuverInformation(satelliteID, new Cartesian(3300.0, 0.0, 0.0), earth.getInertialFrame().getAxes()));
    ImpulsiveManeuverSegment firstBurn = tempObj$2;

    DelegateBasedVariable<ImpulsiveManeuverSegmentConfiguration> velocityXVariable = firstBurn.createVariable(
            100.0, // maximum step, meters/second
            0.1,   // perturbation, meters/second
            SetVariableCallback.of((variable, configuration) -> {
                configuration.getManeuver().setX(configuration.getManeuver().getX() + variable);
            }));

}

@Test
public final void stopSegmentExample() {

    // define constants
    String satelliteName = "Satellite";
    EarthCentralBody earth = CentralBodiesFacet.getFromContext().getEarth();

    // make the segment list
    SegmentList segmentList = new SegmentList();

    // initial state
    NumericalInitialStateSegment initialStateSegment = configureNumericalInitialStateSegment();

    // make the propagation segment
    // we will need the integration point to make the altitude stopping condition
    PropagationNewtonianPoint propagationPoint = createPropagatedNewtonianPoint(satelliteName);
    NumericalPropagatorSegment propagationSegment = new NumericalPropagatorSegment();
    propagationSegment.setPropagatorDefinition(createNumericalPropagator(propagationPoint));

    // make a stopping condition
    ScalarStoppingCondition altitudeStoppingCondition = new ScalarStoppingCondition(
            new ScalarCartographicElement(
                    earth,
                    propagationPoint.getIntegrationPoint(),
                    CartographicElement.HEIGHT),
            100000.0, // 100 kilometers
            1.0,      // stop within a meter
            StopType.THRESHOLD_DECREASING); // the satellite won't come back up

    // the satellite won't come back up
    propagationSegment.getStoppingConditions().add(altitudeStoppingCondition);

    // if the satellite goes below 100 km, then stop propagation
    StopSegment stopIfAltitudeIsTooLow = new StopSegment(true);
    propagationSegment.setStoppingConditionAutoSegment(altitudeStoppingCondition, stopIfAltitudeIsTooLow, 1);

    // or our orbit might not reenter in which case we want to continue propagating another segment.
    StoppingCondition anotherStoppingCondition = createSomeOtherStoppingCondition();
    propagationSegment.getStoppingConditions().add(anotherStoppingCondition);

    // if we don't stop, this segment will propagate after the previous propagation segment finishes
    NumericalPropagatorSegment someOtherPropagationSegment = configureNumericalPropagatorSegment();

    // configure the segment list
    segmentList.getSegments().add(initialStateSegment);
    segmentList.getSegments().add(propagationSegment);
    segmentList.getSegments().add(someOtherPropagationSegment);

}

@Test
public final void returnSegmentExample() {

    // define constants
    String satelliteName = "Satellite";
    EarthCentralBody earth = CentralBodiesFacet.getFromContext().getEarth();

    // make the segment list (create the differential corrector too)
    TargetedSegmentList segmentList = createTargetedSegmentList();

    // initial state
    NumericalInitialStateSegment initialStateSegment = configureNumericalInitialStateSegment();
    ImpulsiveManeuverSegment firstManeuver = createFirstManeuver();

    // make the propagation segment
    // we will need the integration point to make the altitude stopping condition
    PropagationNewtonianPoint propagationPoint = createPropagatedNewtonianPoint(satelliteName);
    NumericalPropagatorSegment propagationSegment = new NumericalPropagatorSegment();
    propagationSegment.setPropagatorDefinition(createNumericalPropagator(propagationPoint));

    // make a stopping condition
    ScalarStoppingCondition altitudeStoppingCondition = new ScalarStoppingCondition(
            new ScalarCartographicElement(
                    earth,
                    propagationPoint.getIntegrationPoint(),
                    CartographicElement.HEIGHT),
            420000000.0, // 42,000 kilometers
            1.0,         // stop within a meter
            StopType.ANY_THRESHOLD);
    propagationSegment.getStoppingConditions().add(altitudeStoppingCondition);

    // if the satellite goes above 42 000 km, then we don't want to do the second maneuver
    ReturnSegment stopIfAltitudeIsTooLow = new ReturnSegment(segmentList, ReturnSegmentBehavior.ENABLED);
    propagationSegment.setStoppingConditionAutoSegment(altitudeStoppingCondition, stopIfAltitudeIsTooLow, 1);

    // or our orbit might still be too low in which case we want to do another maneuver
    StoppingCondition anotherStoppingCondition = createSomeOtherStoppingCondition();
    propagationSegment.getStoppingConditions().add(anotherStoppingCondition);
    ImpulsiveManeuverSegment impulsiveManeuver = createSecondManeuver();

    // configure the segment list
    segmentList.getSegments().add(initialStateSegment);
    segmentList.getSegments().add(firstManeuver);
    segmentList.getSegments().add(propagationSegment);
    segmentList.getSegments().add(impulsiveManeuver);

}

public final void changePropagateSegmentOperatorExample() {

    EarthCentralBody earth = CentralBodiesFacet.getFromContext().getEarth();

    // configure the first low precision propagator
    TwoBodyStoppablePropagator firstLowPrecisionPropagator = new TwoBodyStoppablePropagator(
            WorldGeodeticSystem1984.GravitationalParameter, 
            SatelliteMotionIdentification, 
            earth.getInertialFrame(), 
            Duration.fromMinutes(5.0));
    firstLowPrecisionPropagator.setInitialDate(TimeConstants.J2000);
    firstLowPrecisionPropagator.setInitialMotion(new Motion1<Cartesian>(new Cartesian(42000000.0, 0.0, 0.0), new Cartesian(0.0, 3000.0, 0.0)));

    // set its stopping condition
    DurationStoppingCondition firstCondition = new DurationStoppingCondition();
    firstCondition.setThreshold(ValueDefinition.toValueDefinition(Duration.fromMinutes(10.0)));
    firstLowPrecisionPropagator.getStoppingConditions().add(firstCondition);

    // create the segment
    PropagateSegment segment = new PropagateSegment(firstLowPrecisionPropagator);
    ReferenceFrameAdapter fullyConfiguredAdapter = new ReferenceFrameAdapter(SatelliteMotionIdentification, earth.getInertialFrame());
    fullyConfiguredAdapter.setPreviousReferenceFrame(earth.getInertialFrame()); // note that the previous defined in must be set
    segment.setElementAndAdapter(fullyConfiguredAdapter); // set the adapter

    // create the high precision propagator that will replace the existing propagator
    StoppableNumericalPropagator secondHighPrecisionPropagator = new StoppableNumericalPropagator();
    NumericalPropagatorDefinition numericalPropagator = createNumericalPropagator();
    secondHighPrecisionPropagator.setPropagatorDefinition(numericalPropagator);

    // this needs its own stopping condition too
    DurationStoppingCondition secondCondition = new DurationStoppingCondition();
    secondCondition.setThreshold(firstCondition.getThreshold());
    secondHighPrecisionPropagator.getStoppingConditions().add(secondCondition);

    // create the targeted segment list
    TargetedSegmentList list = new TargetedSegmentList();
    list.getSegments().add(segment);

    // add the operator
    ChangePropagateSegmentsPropagatorOperator changeOperator = new ChangePropagateSegmentsPropagatorOperator(segment, secondHighPrecisionPropagator);
    list.getOperators().add(createDifferentialCorrector1(segment, segment));
    list.getOperators().add(changeOperator);
    list.getOperators().add(createDifferentialCorrector2(segment, segment));

    // propagate
    SegmentListPropagator propagator = list.getSegmentListPropagator();
    TargetedSegmentListResults results = (TargetedSegmentListResults) propagator.propagateSegmentList();
    PropagateSegmentResults highPrecisionResults = (PropagateSegmentResults) results.getResultsOfSegment(segment);

}

private static TargetedSegmentListDifferentialCorrector createDifferentialCorrector1(SegmentDefinition segmentToMakeVariableFrom, SegmentDefinition segmentToMakeConstraintFrom) {
    TargetedSegmentListDifferentialCorrector difeq = new TargetedSegmentListDifferentialCorrector();
    difeq.getVariables().add(new DelegateBasedVariable<SegmentConfiguration>(1.0, 0.1, segmentToMakeVariableFrom, new SetVariableCallback<SegmentConfiguration>() {
        @Override
        public void invoke(double a, SegmentConfiguration b) {}
    }));
    difeq.getConstraints().add(new DelegateBasedConstraint(new DelegateBasedConstraintCallback() {
        @Override
        public double invoke(SegmentResults results) {
            return (JulianDate.subtract(results.getInitialState().getCurrentDate(), TimeConstants.J2000)).getTotalSeconds();
        }
    }, segmentToMakeConstraintFrom, 100.0, 1.0));
    return difeq;
}

private static TargetedSegmentListDifferentialCorrector createDifferentialCorrector2(SegmentDefinition segmentToMakeVariableFrom, SegmentDefinition segmentToMakeConstraintFrom) {
    return createDifferentialCorrector1(segmentToMakeVariableFrom, segmentToMakeConstraintFrom);
}

private static ImpulsiveManeuverSegment createFirstManeuver() {
    return new ImpulsiveManeuverSegment();
}

private static ImpulsiveManeuverSegment createSecondManeuver() {
    return new ImpulsiveManeuverSegment();
}

private static TargetedSegmentList createTargetedSegmentList() {
    return new TargetedSegmentList();
}

private static ScalarStoppingCondition createRangeOfZeroStoppingCondition() {
    return new ScalarStoppingCondition();
}

private static StkEphemerisFile getEphemerisFile() {
    return StkEphemerisFile.readFrom(getRelativeDataPath("Data\\Covariance6x6.e"));
}

private static SegmentList getVerticalCorrectionSegmentList() {
    return new SegmentList();
}

private static SegmentList getHorizontalCorrectionSegmentList() {
    return new SegmentList();
}

public final FiniteManeuverSegment basicFiniteManeuver(@Nonnull Scalar[] fuelFlowRate, @Nonnull DurationStoppingCondition[] primaryDurationStoppingCondition) {

    // CreatePropagatedNewtonianPoint also creates fuel and dry mass elements
    PropagationNewtonianPoint propagationPoint = createPropagatedNewtonianPoint(SatelliteMotionIdentification);
    NumericalPropagatorDefinition basePropagator = createNumericalPropagator(propagationPoint);
    PropagationScalar fuelMassScalar = new PropagationScalar(500.0);
    fuelMassScalar.setIdentification("Fuel Mass");
    AuxiliaryStateScalar dryMass = new AuxiliaryStateScalar(Scalar.toScalar(500.0));
    dryMass.setIdentification("Dry Mass");
    propagationPoint.setMass(Scalar.add(fuelMassScalar.getIntegrationValue(), dryMass.getAuxiliaryScalar()));
    basePropagator.getIntegrationElements().add(fuelMassScalar);
    basePropagator.getAuxiliaryElements().add(dryMass);

    double nominalFuelFlow = -0.02; // kg/sec, negative because fuel will be consumed
    ScalarFixed fuelFlowRateScalar = new ScalarFixed(nominalFuelFlow);
    fuelMassScalar.setScalarDerivative(fuelFlowRateScalar);

    double isp = 250.0; // seconds
    Cartesian thrustDirection = new Cartesian(0.0, 0.0, 1.0);   // unit vector
    Vector thrustVector = ContinuousThrustForceModel.createIspThrustVector(
            Scalar.toScalar(isp), 
            fuelFlowRateScalar, 
            new AxesVelocityOrbitNormal(propagationPoint.getIntegrationPoint()),
            thrustDirection);
    ContinuousThrustForceModel thrustForce = new ContinuousThrustForceModel(thrustVector, propagationPoint.getIntegrationFrame().getAxes());
    propagationPoint.getAppliedForces().add(thrustForce);



    FiniteManeuverSegment segment = new FiniteManeuverSegment();
    segment.setName("Finite Maneuver");
    segment.setPropagatorDefinition(basePropagator);

    DurationStoppingCondition tenMinuteBurn = new DurationStoppingCondition(Duration.fromMinutes(10.0));
    segment.getStoppingConditions().add(tenMinuteBurn);

    // add a stopping condition to run out of fuel
    StoppingCondition outOfFuelStoppingCondition = 
            FiniteManeuverSegment.createOutOfFuelStoppingCondition(fuelMassScalar.getIntegrationValue(), "Out of Fuel", true);
    segment.getStoppingConditions().add(outOfFuelStoppingCondition);

    // propagate
    SegmentPropagator propagator = segment.getSegmentPropagator();
    SegmentResults results = propagator.propagate();

    primaryDurationStoppingCondition[0] = tenMinuteBurn;
    fuelFlowRate[0] = fuelFlowRateScalar;
    AssertHelper.assertNotNull(results, "sanity check");
    return segment;
}

public final FiniteManeuverSegment basicFiniteManeuver() {
    Scalar fuelFlow = null;
    DurationStoppingCondition durationStop = null;
    final Scalar[] out$fuelFlow$4 = new Scalar[1];
    final DurationStoppingCondition[] out$durationStop$5 = new DurationStoppingCondition[1];
    final FiniteManeuverSegment temp$3 = basicFiniteManeuver(out$fuelFlow$4, out$durationStop$5);
    durationStop = out$durationStop$5[0];
    fuelFlow = out$fuelFlow$4[0];
    return temp$3;
}

@Test
public final void burnCenteringExample() {

    InitialStateSegment<BasicState> initialState = configureInitialStateSegment();
    PropagateSegment initialPropagation = createAndConfigurePropagateSegment();

    FiniteManeuverSegment finiteManeuverSegment = basicFiniteManeuver();

    DurationStoppingCondition durationCondition = new DurationStoppingCondition(Duration.fromMinutes(5.0));
    finiteManeuverSegment.getStoppingConditions().add(durationCondition);
    finiteManeuverSegment.setBurnCenteringDuration(durationCondition.getThreshold());

    PropagateSegment finalPropagation = createAndConfigurePropagateSegment();

    SegmentList mcs = new SegmentList();
    mcs.getSegments().add(initialState);
    mcs.getSegments().add(initialPropagation);
    mcs.getSegments().add(finiteManeuverSegment);
    mcs.getSegments().add(finalPropagation);

    // propagate
    SegmentListResults results = mcs.getSegmentListPropagator().propagateSegmentList();
    PropagateSegmentResults firstPropagateSegmentResults = (PropagateSegmentResults) results.getResultsOfSegment(initialPropagation);
    FiniteManeuverSegmentResults maneuverResults = (FiniteManeuverSegmentResults) results.getResultsOfSegment(finiteManeuverSegment);


    Assert.assertNotNull(firstPropagateSegmentResults);
    Assert.assertNotNull(maneuverResults);
    ITimeBasedState finalEvaluatedStateOfFirstPropagate = firstPropagateSegmentResults.getEntireComputedEphemeris().get(firstPropagateSegmentResults.getEntireComputedEphemeris().size() - 1);
    AssertHelper.assertEquals(JulianDate.subtract(finalEvaluatedStateOfFirstPropagate.getCurrentDate(), maneuverResults.getStateAtEngineIgnition().getCurrentDate()), Duration.fromMinutes(2.5));
}

private static PropagationScalar configureDefaultMassesAndReturnFuelElement(PropagationNewtonianPoint propagationPoint, NumericalPropagatorDefinition basePropagator, double fuelMass,
        double dryMass) {

    PropagationScalar fuelMassScalar = new PropagationScalar(fuelMass);
    fuelMassScalar.setIdentification("Fuel Mass");

    AuxiliaryStateScalar dryMassScalar = new AuxiliaryStateScalar(Scalar.toScalar(dryMass));
    dryMassScalar.setIdentification("Dry Mass");

    propagationPoint.setMass(Scalar.add(fuelMassScalar.getIntegrationValue(), dryMassScalar.getAuxiliaryScalar()));

    basePropagator.getIntegrationElements().add(fuelMassScalar);
    basePropagator.getAuxiliaryElements().add(dryMassScalar);

    return fuelMassScalar;
}

@Test
public final void fixedAtTimeAxesExample() {

    FiniteManeuverSegment segment = new FiniteManeuverSegment();
    segment.setName("Finite Maneuver");

    // createPropagatedNewtonianPoint also creates fuel and dry mass elements
    PropagationNewtonianPoint propagationPoint = createPropagatedNewtonianPoint(SatelliteMotionIdentification);
    NumericalPropagatorDefinition basePropagator = createNumericalPropagator(propagationPoint);
    segment.setPropagatorDefinition(basePropagator);
    PropagationScalar fuelMassElement = configureDefaultMassesAndReturnFuelElement(propagationPoint, basePropagator, 500.0, 500.0);

    double nominalFuelFlow = -0.02; // kg/sec, negative because fuel will be consumed
    ScalarFixed fuelFlowRateScalar = new ScalarFixed(nominalFuelFlow);
    fuelMassElement.setScalarDerivative(fuelFlowRateScalar);
    double isp = 250.0; // seconds
    Cartesian thrustDirection = new Cartesian(0.0, 0.0, 1.0); // unit vector

    AxesVelocityOrbitNormal axesThatWillGetFixed = new AxesVelocityOrbitNormal(propagationPoint.getIntegrationPoint());
    AxesFixedAtJulianDate axesFixedAtEngineIgnition = new AxesFixedAtJulianDate(axesThatWillGetFixed, new TimeFromStateValueDefinition(segment.getIgnitionState()));
    Vector thrustVector = ContinuousThrustForceModel.createIspThrustVector(Scalar.toScalar(isp), fuelFlowRateScalar, axesFixedAtEngineIgnition, thrustDirection);
    ContinuousThrustForceModel thrustForce = new ContinuousThrustForceModel(thrustVector, propagationPoint.getIntegrationFrame().getAxes());
    propagationPoint.getAppliedForces().add(thrustForce);

    StoppingCondition tenMinuteBurn = new DurationStoppingCondition(Duration.fromMinutes(10.0));
    segment.getStoppingConditions().add(tenMinuteBurn);

    // add a stopping condition to run out of fuel
    StoppingCondition outOfFuelStoppingCondition = FiniteManeuverSegment.createOutOfFuelStoppingCondition(fuelMassElement.getIntegrationValue(), "Out of Fuel", true);
    segment.getStoppingConditions().add(outOfFuelStoppingCondition);

    SegmentPropagator propagator = segment.getSegmentPropagator();
    SegmentResults results = propagator.propagate();

}

@Test
public final void burnToZeroInclinationExample() {
    double earthGm = WorldGeodeticSystem1984.GravitationalParameter;
    Motion1<Cartesian> initialConditions = new ModifiedKeplerianElements(17000000.0, 1.0 / 17200000.0, Trig.degreesToRadians(23.0), 
            0.0, 0.0, Trig.degreesToRadians(5.0), earthGm).toCartesian();

    PropagationNewtonianPoint propagationPoint = createPropagatedNewtonianPoint(SatelliteMotionIdentification, initialConditions);
    NumericalPropagatorDefinition basePropagator = createNumericalPropagator(propagationPoint);

    double initialActualFuelMass = 5000.0;
    double initialDryMass = 1000.0;
    PropagationScalar fuelMassElement = configureDefaultMassesAndReturnFuelElement(propagationPoint, basePropagator, initialActualFuelMass, initialDryMass);

    ScalarStoppingCondition outOfFuelStoppingCondition = new ScalarStoppingCondition(fuelMassElement.getIntegrationValue(), 0.00000001, 0.000000001,
            StopType.THRESHOLD_DECREASING);
    outOfFuelStoppingCondition.setEnabled(StoppingConditionEnabled.ENABLED);
    ScalarFixed fuelFlow = new ScalarFixed(-0.005);
    fuelMassElement.setScalarDerivative(fuelFlow);



    // create the continuous thrust force model
    ContinuousThrustForceModel thrust = new ContinuousThrustForceModel();
    thrust.setIntegrationAxes(propagationPoint.getIntegrationFrame().getAxes());

    // use the constant thrust/isp thrust vector
    double isp = 1000.0;
    Vector thrustVector = ContinuousThrustForceModel.createIspThrustVector(Scalar.toScalar(isp), 
            fuelFlow,
            propagationPoint.getIntegrationFrame().getAxes(),
            new Cartesian(0.0, 0.0, 1.0)); // thrust along the negative Z
    thrust.setThrustVector(thrustVector);

    Vector velocityVector = new VectorVelocity(propagationPoint.getIntegrationPoint(), propagationPoint.getIntegrationFrame());

    ScalarSwitchingOnVelocityVector switchingScalar = new ScalarSwitchingOnVelocityVector(velocityVector);
    thrust.setThrustVector(Vector.multiply(thrustVector, switchingScalar)); // switch the thrust back and forth

    propagationPoint.getAppliedForces().add(thrust);



    FiniteManeuverSegment segment = new FiniteManeuverSegment();
    segment.setPropagatorDefinition(basePropagator);

    // a time condition just in case
    segment.getStoppingConditions().add(new DurationStoppingCondition(Duration.fromDays(10.0)));

    // don't forget the out of fuel stopping condition we made
    segment.getStoppingConditions().add(outOfFuelStoppingCondition);

    // don't keep going once we've hit our flat orbit
    ScalarStoppingCondition inclinationStoppingCondition = new ScalarStoppingCondition(
            new ScalarModifiedKeplerianElement(
                earthGm, 
                propagationPoint.getIntegrationPoint(),
                KeplerianElement.INCLINATION, 
                propagationPoint.getIntegrationFrame()), 
            0.0, Trig.degreesToRadians(0.05), StopType.ANY_THRESHOLD);
    segment.getStoppingConditions().add(inclinationStoppingCondition);

    SegmentList list = new SegmentList();
    list.getSegments().add(segment);

    DateMotionCollection1<Cartesian> ephemeris = list.getSegmentListPropagator().propagate().getDateMotionCollectionOfOverallTrajectory(SatelliteMotionIdentification,
            propagationPoint.getIntegrationFrame());
    AssertHelper.assertNotEqual(0, ephemeris.getCount(), "sanity check");
}

@Test
public final void seedFiniteManeuverExample() {
    EarthCentralBody earth = getEarth();
    Scalar fuelFlowRate = null;
    DurationStoppingCondition durationStoppingCondition = null;
    final Scalar[] out$fuelFlowRate$6 = new Scalar[1];
    final DurationStoppingCondition[] out$durationStoppingCondition$7 = new DurationStoppingCondition[1];
    FiniteManeuverSegment finiteManeuver = basicFiniteManeuver(out$fuelFlowRate$6, out$durationStoppingCondition$7);
    durationStoppingCondition = out$durationStoppingCondition$7[0];
    fuelFlowRate = out$fuelFlowRate$6[0];


    TargetedSegmentList mcs = new TargetedSegmentList();
    PropagateSegment propToSomePoint = createAndConfigurePropagateSegment();

    ImpulsiveManeuverSegment impulsiveManeuver = createFirstManeuver();

    ImpulsiveManeuverInformation impulsiveManeuverInformation = 
            new ImpulsiveManeuverInformation(SatelliteMotionIdentification, new Cartesian(0.0, 0.0, 100.0));
    impulsiveManeuverInformation.setOrientation(earth.getInertialFrame().getAxes());
    impulsiveManeuver.setManeuver(impulsiveManeuverInformation);

    durationStoppingCondition.setThreshold(new ParameterizedValueDefinition<Duration>());

    SwitchableSegment overallManeuverSegment = new SwitchableSegment(finiteManeuver, impulsiveManeuver);
    PropagateSegment propagateSomeMore = createAndConfigurePropagateSegment();

    mcs.getSegments().add(propToSomePoint);
    mcs.getSegments().add(overallManeuverSegment);
    mcs.getSegments().add(propagateSomeMore);

    TargetedSegmentListDifferentialCorrector firstDifferentialCorrector = createDifferentialCorrector1(impulsiveManeuver, impulsiveManeuver);

    SeedFiniteManeuverOperator seedManeuverOperator = new SeedFiniteManeuverOperator(overallManeuverSegment, SatelliteMotionIdentification);
    seedManeuverOperator.setFuelFlowRateScalar(fuelFlowRate);
    seedManeuverOperator.setStoppingConditionToUpdate(durationStoppingCondition);

    TargetedSegmentListDifferentialCorrector secondDifferentialCorrector = createDifferentialCorrector2(finiteManeuver, finiteManeuver);

    mcs.getOperators().add(firstDifferentialCorrector);
    mcs.getOperators().add(seedManeuverOperator);
    mcs.getOperators().add(secondDifferentialCorrector);

    Assert.assertNotNull(mcs.getSegmentListPropagator().propagateSegmentList());
}

@Test
public final void swapSegmentExample() {
    TargetedSegmentList mcs = new TargetedSegmentList();

    PropagateSegment propagateWithoutManeuver = createAndConfigurePropagateSegment();
    FiniteManeuverSegment propagateWithManeuver = basicFiniteManeuver();

    SwitchableSegment overallManeuverSegment = new SwitchableSegment(propagateWithManeuver, propagateWithoutManeuver);
    PropagateSegment propagateSomeMore = createAndConfigurePropagateSegment();

    mcs.getSegments().add(overallManeuverSegment);
    mcs.getSegments().add(propagateSomeMore);

    TargetedSegmentListDifferentialCorrector firstDifferentialCorrector = createDifferentialCorrector1(propagateWithoutManeuver, propagateWithoutManeuver);

    SwapSegmentOperator swapSegmentOperator = new SwapSegmentOperator(overallManeuverSegment);

    TargetedSegmentListDifferentialCorrector secondDifferentialCorrector = createDifferentialCorrector2(propagateWithManeuver, propagateWithManeuver);

    mcs.getOperators().add(firstDifferentialCorrector);
    mcs.getOperators().add(swapSegmentOperator);
    mcs.getOperators().add(secondDifferentialCorrector);

    Assert.assertNotNull(mcs.getSegmentListPropagator().propagateSegmentList());
}

@Test
public final void primerVectorExample() {
    // We want orbits with less than 1 complete revolution.
    final int numberOfRevolutions = 0;

    // Define central body parameters, reference frame, and initial date.
    EarthCentralBody earth = CentralBodiesFacet.getFromContext().getEarth();
    ReferenceFrame icrfFrame = earth.getInternationalCelestialReferenceFrame();
    double gravitationalParameter = EarthGravitationalModel2008.GravitationalParameter; // 3.986004415e+14 m^3/s^2
    double earthRadius = EarthGravitationalModel2008.SemimajorAxis; // 6378136.3 meters
    JulianDate initialDate = new GregorianDate(2022, 3, 16, 0, 0, 0.0).toJulianDate();

    // Parameters specific to initial LEO.
    double initialSemimajorAxis = 300000.0 + earthRadius; // Circular orbit with 300 km altitude. 
    double initialEccentricity = 0.0;
    double initialInclination = 0.0;
    double initialRightAscensionAscendingNode = 0.0;
    double initialArgumentOfPerigee = 0.0;
    double initialTrueAnomaly = 0.0;
    Motion1<Cartesian> initialLeo = new KeplerianElements(initialSemimajorAxis, initialEccentricity,
            initialInclination, initialArgumentOfPerigee, initialRightAscensionAscendingNode,
            initialTrueAnomaly, gravitationalParameter).toCartesian();

    // Parameters specific to final GEO.
    Motion1<Cartesian> finalGeo = new Motion1<Cartesian>(new Cartesian(-42165000.0, 0.0, 0.0), // position vector, meters.
            new Cartesian(0.0, -3074.6, 0.0)); // velocity vector, m/s.
    // Time of flight of GTO transfer orbit.
    double gtoSemimajorAxis = (initialSemimajorAxis + finalGeo.getValue().getMagnitude()) / 2.0; // meters.
    double gtoTimeOfFlight = 0.5 * OrbitalElements.computePeriod(1.0 / gtoSemimajorAxis, gravitationalParameter); // seconds.

    // Find solution using Lambert solver.
    LambertOrbitSolver lambertSolver = new LambertOrbitSolver(gravitationalParameter);
    LambertResult lambertSolution = lambertSolver.solveFixedDurationTransfer(initialLeo.getValue(), finalGeo.getValue(),
            Duration.fromSeconds(gtoTimeOfFlight), numberOfRevolutions, LambertPathType.SHORT_PATH,
            OrbitDirectionType.PROGRADE, initialLeo.getFirstDerivative());

    // Calculate delta-v's required to complete Lambert solution.
    Cartesian deltaV1 = Cartesian.subtract(lambertSolution.getInitialPositionMotion().getFirstDerivative(), initialLeo.getFirstDerivative());
    Cartesian deltaV2 = Cartesian.subtract(finalGeo.getFirstDerivative(), lambertSolution.getFinalPositionMotion().getFirstDerivative());

    double deltaV1Magnitude = deltaV1.getMagnitude();
    double deltaV2Magnitude = deltaV2.getMagnitude();
    UnitCartesian initialPrimerVector = deltaV1.normalize();
    UnitCartesian finalPrimerVector = deltaV2.normalize();

    // Setup preliminary propagator to calculate state transition matrix.
    NumericalPropagatorDefinition stmPropagatorDefinition = new NumericalPropagatorDefinition();
    stmPropagatorDefinition.setEpoch(initialDate);

    double initialMass = 1000.0; // kg.
    String massElementId = "Mass";
    PropagationScalar mass = new PropagationScalar(initialMass);
    mass.setIdentification(massElementId);
    mass.setScalarDerivative(Scalar.toScalar(0.0));

    // Define states using PropagationNewtonianPoint.
    final String id = "PrimaryMotion";
    PropagationNewtonianPoint point = new PropagationNewtonianPoint();
    point.setIdentification(id);
    point.setMass(mass.getIntegrationValue());
    point.setInitialPosition(initialLeo.getValue());
    point.setInitialVelocity(Cartesian.add(initialLeo.getFirstDerivative(), deltaV1));
    point.setIntegrationFrame(icrfFrame);

    // Configure two-body gravity.
    TwoBodyGravity twoBodyGravity = new TwoBodyGravity(point.getIntegrationPoint(), earth, gravitationalParameter);
    point.getAppliedForces().add(twoBodyGravity);
    stmPropagatorDefinition.getIntegrationElements().add(mass);
    stmPropagatorDefinition.getIntegrationElements().add(point);

    // Create state transition matrix and add it as an integration element.
    final String stmId = "StateTransitionMatrix";
    IPartialDifferentiable position = (IPartialDifferentiable)point.getIntegrationPoint();
    IPartialDifferentiable velocity = (IPartialDifferentiable)point.getIntegrationPoint().createVectorVelocity(point.getIntegrationFrame());
    // point is the IPartialDifferentiable for the acceleration.

    StateTransitionMatrix stateMatrix = new StateTransitionMatrix();
    stateMatrix.setIdentification(stmId);
    stateMatrix.addStateParameter(position, velocity);
    stateMatrix.addStateParameter(velocity, point);
    stmPropagatorDefinition.getIntegrationElements().add(stateMatrix);

    // Configure numerical integrator.
    RungeKuttaFehlberg78Integrator integrator = new RungeKuttaFehlberg78Integrator();
    integrator.setRelativeTolerance(Constants.Epsilon13);
    integrator.setMinimumStepSize(1.0);
    integrator.setMaximumStepSize(6.0 * 3600.0);
    integrator.setInitialStepSize(60.0);
    integrator.setAbsoluteTolerance(Constants.Epsilon13);
    integrator.setStepSizeBehavior(KindOfStepSize.RELATIVE);
    stmPropagatorDefinition.setIntegrator(integrator);

    // Propagate and extract results.
    NumericalPropagator stmPropagator = stmPropagatorDefinition.createPropagator();
    NumericalPropagationStateHistory stmPropagatorResults = stmPropagator.propagate(Duration.fromSeconds(gtoTimeOfFlight), 1);
    List<Matrix> stmHistory = stmPropagatorResults.<Matrix> getDateMotionCollection(stmId).getValues();
    int stmHistoryLength = stmHistory.size();
    Matrix finalStm = stmHistory.get(stmHistoryLength - 1);

    // Extract dPosition/dPosition0 submatrix.
    Matrix finalStmSubmatrix1 = finalStm.getMatrix(0, 3, 0, 3);
    Matrix3By3 convertedSubmatrix1 
        = new Matrix3By3(finalStmSubmatrix1.get(0, 0), finalStmSubmatrix1.get(0, 1), finalStmSubmatrix1.get(0, 2),
                         finalStmSubmatrix1.get(1, 0), finalStmSubmatrix1.get(1, 1), finalStmSubmatrix1.get(1, 2),
                         finalStmSubmatrix1.get(2, 0), finalStmSubmatrix1.get(2, 1), finalStmSubmatrix1.get(2, 2));

    // Extract dPosition/dVelocity0 submatrix.
    Matrix finalStmSubmatrix2 = finalStm.getMatrix(0, 3, 3, 3);
    Matrix3By3 convertedSubmatrix2 
        = new Matrix3By3(finalStmSubmatrix2.get(0, 0), finalStmSubmatrix2.get(0, 1), finalStmSubmatrix2.get(0, 2), 
                         finalStmSubmatrix2.get(1, 0), finalStmSubmatrix2.get(1, 1), finalStmSubmatrix2.get(1, 2),
                         finalStmSubmatrix2.get(2, 0), finalStmSubmatrix2.get(2, 1), finalStmSubmatrix2.get(2, 2));

    // Solve for initial primer vector rate.
    Cartesian rightHandSide = finalPrimerVector.subtract(
            Matrix3By3.multiply(convertedSubmatrix1, Cartesian.toCartesian(initialPrimerVector)));
    Cartesian initialPrimerVectorRate = Matrix3By3.multiply(convertedSubmatrix2.invert(), rightHandSide);

    // Now set-up main propagator that includes primer-vector propagation and maneuvers.
    point.setInitialVelocity(initialLeo.getFirstDerivative()); // Needs to be reset to no longer include deltaV1.

    // Define PrimerVector as a costate PropagationVector
    String primerVectorId = "PrimerVector";
    PropagationVector primerVector = new PropagationVector(Cartesian.toCartesian(initialPrimerVector), initialPrimerVectorRate);
    primerVector.setIdentification(primerVectorId);
    primerVector.setPropagationAxes(icrfFrame.getAxes());
    // The VectorDerivative is the second derivative and is set in the following code.

    // Get partial derivatives of the gravity.
    Scalar radius = point.getIntegrationPoint().getScalarElement(CartesianElement.MAGNITUDE, point.getIntegrationFrame());
    Scalar radiusX = point.getIntegrationPoint().getScalarElement(CartesianElement.X, point.getIntegrationFrame());
    Scalar radiusY = point.getIntegrationPoint().getScalarElement(CartesianElement.Y, point.getIntegrationFrame());
    Scalar radiusZ = point.getIntegrationPoint().getScalarElement(CartesianElement.Z, point.getIntegrationFrame());
    Scalar radiusSquared = Scalar.multiply(radius, radius);
    Scalar factor = Scalar.divide(gravitationalParameter, (Scalar.multiply(radiusSquared, radius)));

    // Two-body gravity has simple partial derivatives. Even then, these operations only work because
    // all of the states and co-states are in the same reference frame.
    Scalar dGravityX_dX = Scalar.multiply(factor, (Scalar.subtract(Scalar.divide(Scalar.multiply(
            Scalar.multiply(3.0, radiusX), radiusX), radiusSquared), Scalar.toScalar(1.0))));
    Scalar dGravityY_dY = Scalar.multiply(factor, (Scalar.subtract(Scalar.divide(Scalar.multiply(
            Scalar.multiply(3.0, radiusY), radiusY), radiusSquared), Scalar.toScalar(1.0))));
    Scalar dGravityZ_dZ = Scalar.multiply(factor, (Scalar.subtract(Scalar.divide(Scalar.multiply(
            Scalar.multiply(3.0, radiusZ), radiusZ), radiusSquared), Scalar.toScalar(1.0))));
    // Equal to dGravityY_dX
    Scalar dGravityX_dY = Scalar.divide(Scalar.multiply(Scalar.multiply(Scalar.multiply(factor, 3.0), radiusX), radiusY), radiusSquared);
    // Equal to dGravityZ_dX
    Scalar dGravityX_dZ = Scalar.divide(Scalar.multiply(Scalar.multiply(Scalar.multiply(factor, 3.0), radiusX), radiusZ), radiusSquared);
    // Equal to dGravityZ_dY
    Scalar dGravityY_dZ = Scalar.divide(Scalar.multiply(Scalar.multiply(Scalar.multiply(factor, 3.0), radiusY), radiusZ), radiusSquared);

    // Primer vector is multiplied by partials matrix to get derivatives.
    Scalar primerVectorX = primerVector.getIntegrationValue().getScalarElement(CartesianElement.X, icrfFrame.getAxes());
    Scalar primerVectorY = primerVector.getIntegrationValue().getScalarElement(CartesianElement.Y, icrfFrame.getAxes());
    Scalar primerVectorZ = primerVector.getIntegrationValue().getScalarElement(CartesianElement.Z, icrfFrame.getAxes());
    Scalar primerSecondDerivativeX = Scalar.add(Scalar.add(Scalar.multiply(primerVectorX, dGravityX_dX),
                                                           Scalar.multiply(primerVectorY, dGravityX_dY)),
                                                           Scalar.multiply(primerVectorZ, dGravityX_dZ));
    Scalar primerSecondDerivativeY = Scalar.add(Scalar.add(Scalar.multiply(primerVectorX, dGravityX_dY),
                                                           Scalar.multiply(primerVectorY, dGravityY_dY)),
                                                           Scalar.multiply(primerVectorZ, dGravityY_dZ));
    Scalar primerSecondDerivativeZ = Scalar.add(Scalar.add(Scalar.multiply(primerVectorX, dGravityX_dZ),
                                                           Scalar.multiply(primerVectorY, dGravityY_dZ)),
                                                           Scalar.multiply(primerVectorZ, dGravityZ_dZ));
    Vector primerSecondDerivative = Vector.add(
            Vector.add(Vector.multiply(primerSecondDerivativeX, VectorFixed.getUnitXVector(icrfFrame.getAxes())),
                    Vector.multiply(primerSecondDerivativeY, VectorFixed.getUnitYVector(icrfFrame.getAxes()))),
            Vector.multiply(primerSecondDerivativeZ, VectorFixed.getUnitZVector(icrfFrame.getAxes())));
    primerVector.setVectorDerivative(primerSecondDerivative);

    // Set-up main propagator.
    NumericalPropagatorDefinition mainPropagatorDefinition = new NumericalPropagatorDefinition();
    mainPropagatorDefinition.setEpoch(initialDate);
    mainPropagatorDefinition.setIntegrator(integrator);
    mainPropagatorDefinition.getIntegrationElements().add(mass);
    mainPropagatorDefinition.getIntegrationElements().add(point);
    mainPropagatorDefinition.getIntegrationElements().add(primerVector);

    // Set-up segment propagation.
    NumericalInitialStateSegment initial = new NumericalInitialStateSegment();
    initial.setName("Initial State");
    initial.setPropagatorDefinition(mainPropagatorDefinition);

    // To illustrate the use of primer vectors, the delta-v will be purely in the x-direction of the maneuver's orientation axes.
    ImpulsiveManeuverSegment maneuver1 = new ImpulsiveManeuverSegment();
    maneuver1.setManeuver(new ImpulsiveManeuverInformation(id, new Cartesian(deltaV1Magnitude, 0.0, 0.0)));
    maneuver1.getManeuver().setCostateElementIdentification(primerVectorId);

    // The maneuver's orientation axes are defined using the costate vector of the maneuver, which is the primer vector.
    AxesAlignedConstrained axes1 = new AxesAlignedConstrained(maneuver1.getManeuver().getCostateVector(),
            new VectorTrueDisplacement(earth.getCenterOfMassPoint(), maneuver1.getManeuver().getPropagationPoint()));
    maneuver1.getManeuver().setOrientation(axes1);

    // The GTO transfer orbit is defined next.
    NumericalPropagatorSegment gtoSegment = new NumericalPropagatorSegment();
    gtoSegment.setName("GTO coast segment");
    gtoSegment.setPropagatorDefinition(mainPropagatorDefinition);
    gtoSegment.setMaximumDuration(Duration.fromDays(1.0));

    // Stop after the GTO is finished.
    DurationStoppingCondition durationStop = new DurationStoppingCondition(Duration.fromSeconds(gtoTimeOfFlight));
    gtoSegment.getStoppingConditions().add(durationStop);

    // To illustrate the use of primer vectors, the delta-v will be purely in the x-direction of the maneuver's orientation axes.
    ImpulsiveManeuverSegment maneuver2 = new ImpulsiveManeuverSegment();
    maneuver2.setManeuver(new ImpulsiveManeuverInformation(id, new Cartesian(deltaV2Magnitude, 0.0, 0.0)));
    maneuver2.getManeuver().setCostateElementIdentification(primerVectorId);

    // The maneuver's orientation axes are defined using the costate vector of the maneuver, which is the primer vector.
    AxesAlignedConstrained axes2 = new AxesAlignedConstrained(maneuver2.getManeuver().getCostateVector(),
            new VectorTrueDisplacement(earth.getCenterOfMassPoint(), maneuver2.getManeuver().getPropagationPoint()));
    maneuver2.getManeuver().setOrientation(axes2);

    // Make a segment list, populate the segment list, and propagate the segment list.
    SegmentList mainSegmentList = new SegmentList();
    mainSegmentList.getSegments().add(initial);
    mainSegmentList.getSegments().add(maneuver1);
    mainSegmentList.getSegments().add(gtoSegment);
    mainSegmentList.getSegments().add(maneuver2);

    // Define propagator and propagate it.
    SegmentListPropagator mainPropagator = mainSegmentList.getSegmentListPropagator();
    SegmentListResults mainResults = mainPropagator.propagateSegmentList();

    // Extract results for each segment inside segment list.
    SegmentResults initialSegmentResults = mainResults.getResultsOfSegment(initial);
    ImpulsiveManeuverSegmentResults maneuver1Results = (ImpulsiveManeuverSegmentResults) mainResults.getResultsOfSegment(maneuver1);
    SegmentResults gtoSegmentResults = mainResults.getResultsOfSegment(gtoSegment);
    ImpulsiveManeuverSegmentResults maneuver2Results = (ImpulsiveManeuverSegmentResults) mainResults.getResultsOfSegment(maneuver2);

    Motion1<Cartesian> propagationInitialState = initialSegmentResults.getFirstPropagatedState().<Cartesian> getMotion(id);

    Motion1<Cartesian> maneuver1InitialState = maneuver1Results.getFirstPropagatedState().<Cartesian> getMotion(id);
    Motion1<Cartesian> maneuver1FinalState = maneuver1Results.getFinalPropagatedState().<Cartesian> getMotion(id);

    Motion1<Cartesian> gtoInitialState = gtoSegmentResults.getFirstPropagatedState().<Cartesian> getMotion(id);
    Motion1<Cartesian> gtoFinalState = gtoSegmentResults.getFinalPropagatedState().<Cartesian> getMotion(id);

    Motion1<Cartesian> gtoInitialCostate = gtoSegmentResults.getFirstPropagatedState().<Cartesian> getMotion(primerVectorId);
    Motion1<Cartesian> gtoFinalCostate = gtoSegmentResults.getFinalPropagatedState().<Cartesian> getMotion(primerVectorId);

    Motion1<Cartesian> maneuver2InitialState = maneuver2Results.getFirstPropagatedState().<Cartesian> getMotion(id);
    Motion1<Cartesian> maneuver2FinalState = maneuver2Results.getFinalPropagatedState().<Cartesian> getMotion(id);
    AssertHelper.assertEquals(initialLeo.getValue().getX(), propagationInitialState.getValue().getX());
    AssertHelper.assertEquals(initialLeo.getValue().getY(), propagationInitialState.getValue().getY());
    AssertHelper.assertEquals(initialLeo.getValue().getZ(), propagationInitialState.getValue().getZ());
    AssertHelper.assertEquals(initialLeo.getFirstDerivative().getX(), propagationInitialState.getFirstDerivative().getX());
    AssertHelper.assertEquals(initialLeo.getFirstDerivative().getY(), propagationInitialState.getFirstDerivative().getY());
    AssertHelper.assertEquals(initialLeo.getFirstDerivative().getZ(), propagationInitialState.getFirstDerivative().getZ());
    AssertHelper.assertEquals(initialLeo.getValue().getX(), maneuver1InitialState.getValue().getX());
    AssertHelper.assertEquals(initialLeo.getValue().getY(), maneuver1InitialState.getValue().getY());
    AssertHelper.assertEquals(initialLeo.getValue().getZ(), maneuver1InitialState.getValue().getZ());
    AssertHelper.assertEquals(initialLeo.getFirstDerivative().getX(), maneuver1InitialState.getFirstDerivative().getX());
    AssertHelper.assertEquals(initialLeo.getFirstDerivative().getY(), maneuver1InitialState.getFirstDerivative().getY());
    AssertHelper.assertEquals(initialLeo.getFirstDerivative().getZ(), maneuver1InitialState.getFirstDerivative().getZ());
    AssertHelper.assertEquals(initialLeo.getValue().getX(), maneuver1FinalState.getValue().getX());
    AssertHelper.assertEquals(initialLeo.getValue().getY(), maneuver1FinalState.getValue().getY());
    AssertHelper.assertEquals(initialLeo.getValue().getZ(), maneuver1FinalState.getValue().getZ());
    AssertHelper.assertEquals(initialLeo.getFirstDerivative().getX() + deltaV1.getX(), maneuver1FinalState.getFirstDerivative().getX(), Constants.Epsilon12);
    AssertHelper.assertEquals(initialLeo.getFirstDerivative().getY() + deltaV1.getY(), maneuver1FinalState.getFirstDerivative().getY(), Constants.Epsilon12);
    AssertHelper.assertEquals(initialLeo.getFirstDerivative().getZ() + deltaV1.getZ(), maneuver1FinalState.getFirstDerivative().getZ(), Constants.Epsilon12);
    AssertHelper.assertEquals(initialLeo.getValue().getX(), gtoInitialState.getValue().getX());
    AssertHelper.assertEquals(initialLeo.getValue().getY(), gtoInitialState.getValue().getY());
    AssertHelper.assertEquals(initialLeo.getValue().getZ(), gtoInitialState.getValue().getZ());
    AssertHelper.assertEquals(initialLeo.getFirstDerivative().getX() + deltaV1.getX(), gtoInitialState.getFirstDerivative().getX(), Constants.Epsilon12);
    AssertHelper.assertEquals(initialLeo.getFirstDerivative().getY() + deltaV1.getY(), gtoInitialState.getFirstDerivative().getY(), Constants.Epsilon12);
    AssertHelper.assertEquals(initialLeo.getFirstDerivative().getZ() + deltaV1.getZ(), gtoInitialState.getFirstDerivative().getZ(), Constants.Epsilon12);
    AssertHelper.assertEquals(finalGeo.getValue().getX(), gtoFinalState.getValue().getX(), Constants.Epsilon4);
    AssertHelper.assertEquals(finalGeo.getValue().getY(), gtoFinalState.getValue().getY(), 60.0);
    AssertHelper.assertEquals(finalGeo.getValue().getZ(), gtoFinalState.getValue().getZ(), Constants.Epsilon14);
    AssertHelper.assertEquals(finalGeo.getFirstDerivative().getX() - deltaV2.getX(), gtoFinalState.getFirstDerivative().getX(), Constants.Epsilon2);
    AssertHelper.assertEquals(finalGeo.getFirstDerivative().getY() - deltaV2.getY(), gtoFinalState.getFirstDerivative().getY(), Constants.Epsilon8);
    AssertHelper.assertEquals(finalGeo.getFirstDerivative().getZ() - deltaV2.getZ(), gtoFinalState.getFirstDerivative().getZ(), Constants.Epsilon12);
    AssertHelper.assertEquals(initialPrimerVector.getX(), gtoInitialCostate.getValue().getX());
    AssertHelper.assertEquals(initialPrimerVector.getY(), gtoInitialCostate.getValue().getY());
    AssertHelper.assertEquals(initialPrimerVector.getZ(), gtoInitialCostate.getValue().getZ());
    AssertHelper.assertEquals(initialPrimerVectorRate.getX(), gtoInitialCostate.getFirstDerivative().getX());
    AssertHelper.assertEquals(initialPrimerVectorRate.getY(), gtoInitialCostate.getFirstDerivative().getY());
    AssertHelper.assertEquals(initialPrimerVectorRate.getZ(), gtoInitialCostate.getFirstDerivative().getZ());
    AssertHelper.assertEquals(finalPrimerVector.getX(), gtoFinalCostate.getValue().getX(), Constants.Epsilon10);
    AssertHelper.assertEquals(finalPrimerVector.getY(), gtoFinalCostate.getValue().getY(), Constants.Epsilon11);
    AssertHelper.assertEquals(finalPrimerVector.getZ(), gtoFinalCostate.getValue().getZ(), Constants.Epsilon15);
    AssertHelper.assertEquals(finalGeo.getValue().getX(), maneuver2InitialState.getValue().getX(), Constants.Epsilon4);
    AssertHelper.assertEquals(finalGeo.getValue().getY(), maneuver2InitialState.getValue().getY(), 60.0);
    AssertHelper.assertEquals(finalGeo.getValue().getZ(), maneuver2InitialState.getValue().getZ(), Constants.Epsilon14);
    AssertHelper.assertEquals(finalGeo.getFirstDerivative().getX() - deltaV2.getX(), maneuver2InitialState.getFirstDerivative().getX(), Constants.Epsilon2);
    AssertHelper.assertEquals(finalGeo.getFirstDerivative().getY() - deltaV2.getY(), maneuver2InitialState.getFirstDerivative().getY(), Constants.Epsilon8);
    AssertHelper.assertEquals(finalGeo.getFirstDerivative().getZ() - deltaV2.getZ(), maneuver2InitialState.getFirstDerivative().getZ(), Constants.Epsilon13);
    AssertHelper.assertEquals(finalGeo.getValue().getX(), maneuver2FinalState.getValue().getX(), Constants.Epsilon4);
    AssertHelper.assertEquals(finalGeo.getValue().getY(), maneuver2FinalState.getValue().getY(), 60.0);
    AssertHelper.assertEquals(finalGeo.getValue().getZ(), maneuver2FinalState.getValue().getZ(), Constants.Epsilon14);
    AssertHelper.assertEquals(finalGeo.getFirstDerivative().getX(), maneuver2FinalState.getFirstDerivative().getX(), Constants.Epsilon2);
    AssertHelper.assertEquals(finalGeo.getFirstDerivative().getY(), maneuver2FinalState.getFirstDerivative().getY(), Constants.Epsilon8);
    AssertHelper.assertEquals(finalGeo.getFirstDerivative().getZ(), maneuver2FinalState.getFirstDerivative().getZ(), Constants.Epsilon13);
}


public static class ScalarSwitchingOnVelocityVector extends Scalar {
    public ScalarSwitchingOnVelocityVector() {}

    public ScalarSwitchingOnVelocityVector(Vector vector) {
        m_vector = vector;
    }

    protected ScalarSwitchingOnVelocityVector(ScalarSwitchingOnVelocityVector existingInstance, CopyContext context) {
        super(existingInstance, context);
        m_vector = context.updateReference(existingInstance.m_vector);
    }

    @Override
    public Object clone(CopyContext context) {
        return new ScalarSwitchingOnVelocityVector(this, context);
    }

    @Override
    protected boolean checkForSameDefinition(Scalar other) {
        ScalarSwitchingOnVelocityVector o = other instanceof ScalarSwitchingOnVelocityVector ? (ScalarSwitchingOnVelocityVector) other : null;
        return o != null && areSameDefinition(m_vector, o.m_vector) && checkForSameDefinition(o);
    }

    protected boolean checkForSameDefinition(ScalarSwitchingOnVelocityVector other) {
        return other != null && other.getClass() == ScalarSwitchingOnVelocityVector.class;
    }

    @Override
    protected int computeCurrentDefinitionHashCode() {
        return HashCode.combine(ScalarSwitchingOnVelocityVector.class.hashCode(),
                                super.computeCurrentDefinitionHashCode(),
                                getDefinitionHashCode(m_vector));
    }

    @Override
    public void enumerateDependencies(DependencyEnumerator enumerator) {
        super.enumerateDependencies(enumerator);
        enumerator.enumerate(m_vector);
    }

    @Override
    public ScalarEvaluator getEvaluator(EvaluatorGroup group) {
        if (group == null) 
            throw new ArgumentNullException("group");

        PropertyInvalidException.validateNonNull(m_vector, "Vector");
        return group.createEvaluator(createEvaluatorCallback);
    }

    private final EvaluatorGroup.Callback0<ScalarEvaluator> createEvaluatorCallback = EvaluatorGroup.Callback0.of(this::createEvaluator);

    private final ScalarEvaluator createEvaluator(EvaluatorGroup group) {
        return new Evaluator(group, m_vector.getEvaluator(group));
    }

    public final Vector getVector() {
        return m_vector;
    }

    public final void setVector(Vector value) {
        throwIfFrozen();
        m_vector = value;
    }

    private Vector m_vector;

    private static final class Evaluator extends ScalarEvaluator {
        public Evaluator(EvaluatorGroup group, VectorEvaluator evaluator) {
            super(group);
            m_evaluator = evaluator;
        }

        private Evaluator(Evaluator existingInstance, CopyContext context) {
            super(existingInstance, context);
            m_evaluator = existingInstance.m_evaluator;
            updateEvaluatorReferences(context);
        }

        @Override
        public Object clone(CopyContext context) {
            return new Evaluator(this, context);
        }

        @Override
        protected void dispose(boolean disposing) {
            if (m_isDisposed || !disposing) {
                return;
            }
            m_isDisposed = true;
            m_evaluator.dispose();
        }

        @Override
        public boolean getIsThreadSafe() {
            return false;
        }

        @Override
        public void updateEvaluatorReferences(CopyContext context) {
            m_evaluator = context.updateReference(m_evaluator);
        }

        @Override
        public TimeIntervalCollection getAvailabilityIntervals(TimeIntervalCollection consideredIntervals) {
            return m_evaluator.getAvailabilityIntervals(consideredIntervals);
        }

        @Override
        public boolean isAvailable(JulianDate date) {
            return m_evaluator.isAvailable(date);
        }

        @Override
        public boolean getIsTimeVarying() {
            return m_evaluator.getIsTimeVarying();
        }

        @Override
        public Motion1<Double> evaluate(JulianDate date, int order) {
            Motion1<Cartesian> motion = m_evaluator.evaluate(date, order + 1);

            double value = motion.getValue().getZ();
            if (value == 0) {
                return new Motion1<Double>(0.0);
            }

            // throttle down as our Z speed gets low
            double multiple = -1.0;
            if (Math.abs(value) < 10.0) {
                multiple = multiple * value * 0.1;
            }

            return new Motion1<Double>(multiple * Math.signum(value));
        }

        private boolean m_isDisposed;
        private VectorEvaluator m_evaluator;
    }
}


public final void freeReturnExample(@Nonnull LunarTrajectoryType type) {
    String pathToEarthGravFile = EGM2008GravityPath.getFullPath();
    String pathToJplDe440 = JplDE440Path.getFullPath();
    String pathToMoonGravFile = new DataPath(DataPathRoot.STK_INSTALL_HOME, "STKData\\CentralBodies\\Moon\\GL0660B.grv").getFullPath();

    // Pick names.
    final String motionID = "Satellite";
    final 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().getMoon();
    SphericalHarmonicGravityModel moonGravityModel = SphericalHarmonicGravityModel.readFrom(pathToMoonGravFile);
    double moonGravitationalParameter = moonGravityModel.getGravitationalParameter(); // m^3/s^2.
    double moonPhysicalRadius = moonGravityModel.getReferenceDistance(); // Meters.

    CentralBody sun = CentralBodiesFacet.getFromContext().getSun();
    final double sunGravitationalParameter = 1.3271244004193938e20; // m^3/s^2.

    CentralBody earth = CentralBodiesFacet.getFromContext().getEarth();
    SphericalHarmonicGravityModel earthGravityModel = SphericalHarmonicGravityModel.readFrom(pathToEarthGravFile);
    final double earthGravitationalParameter = earthGravityModel.getGravitationalParameter();
    double earthPhysicalRadius = earthGravityModel.getReferenceDistance(); // Meters.

    // Set up rotating reference frames.
    VectorTrueDisplacement earthMoonVector = new VectorTrueDisplacement(earth.getCenterOfMassPoint(), moon.getCenterOfMassPoint());
    VectorVelocity moonOrbitalVelocity = new VectorVelocity(moon.getCenterOfMassPoint(), getEarth().getInertialFrame());
    VectorCrossProduct moonOrbitalAngularMomentum = new VectorCrossProduct(earthMoonVector, moonOrbitalVelocity);
    AxesAlignedConstrained rotatingAxes = new AxesAlignedConstrained(earthMoonVector, AxisIndicator.FIRST, moonOrbitalAngularMomentum, AxisIndicator.THIRD);

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

    double earthMoonMassRatio = earthGravitationalParameter / moonGravitationalParameter;
    double scaleFactor = 1.0 / (1.0 + earthMoonMassRatio);
    VectorScaled earthToBarycenterVector = new VectorScaled(earthMoonVector, Scalar.toScalar(scaleFactor));
    Point barycenter = new PointVectorToPoint(earthToBarycenterVector, earth.getCenterOfMassPoint());
    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.getInertialFrame(), group);

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

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

    Cartesian initialPositionMoonRotating = new Cartesian();
    Cartesian initialVelocityMoonRotating = new Cartesian();

    switch (type) {
        case CIRCUMLUNAR: {
            // 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);
            break;
        }
        case CISLUNAR: {
            // 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);
            break;
        }
        default: {
            throw new ArgumentOutOfRangeException();
        }
    }

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

    Motion1<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.setIdentification(dryMassID);
    dryMass.setScalarDerivative(new ScalarFixed(0.0));

    // Create the integrator.
    RungeKuttaFehlberg78Integrator integrator = new RungeKuttaFehlberg78Integrator();
    integrator.setInitialStepSize(60.0);
    integrator.setMinimumStepSize(1.0);
    integrator.setMaximumStepSize(86400.0);

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

    propagationPointAroundMoon.setInitialPosition(initialStateMoonInertial.getValue());
    propagationPointAroundMoon.setInitialVelocity(initialStateMoonInertial.getFirstDerivative());
    propagationPointAroundMoon.setMass(dryMass.getIntegrationValue());

    propagationPointAroundMoon.setIntegrationFrame(moon.getInertialFrame());

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

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

    propagationPointAroundMoon.getAppliedForces().add(moonPrimaryGravity);

    // The Earth and Sun produce third body gravity for Moon propagator.
    ThirdBodyGravity thirdBodyForMoonPropagator = new ThirdBodyGravity(propagationPointAroundMoon.getIntegrationPoint());
    thirdBodyForMoonPropagator.addThirdBody(earth.getName(), earth.getCenterOfMassPoint(), earthGravitationalParameter);
    thirdBodyForMoonPropagator.addThirdBody(sun.getName(), sun.getCenterOfMassPoint(), sunGravitationalParameter);
    thirdBodyForMoonPropagator.setCentralBody(moon);

    propagationPointAroundMoon.getAppliedForces().add(thirdBodyForMoonPropagator);

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

    // Add the elements.
    nearMoonNumericalPropagator.getIntegrationElements().add(propagationPointAroundMoon);
    nearMoonNumericalPropagator.getIntegrationElements().add(dryMass);

    // Set the integrator.
    nearMoonNumericalPropagator.setIntegrator(integrator);

    // Set an epoch.
    nearMoonNumericalPropagator.setEpoch(epoch);

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

    // Make a segment that propagates backwards to edge of Lunar sphere of influence (SOI).
    NumericalPropagatorSegment propagateBackwardToSoi = new NumericalPropagatorSegment();
    propagateBackwardToSoi.setName("Propagate_Backward_ToSOI");
    propagateBackwardToSoi.setPropagatorDefinition(nearMoonNumericalPropagator);
    propagateBackwardToSoi.setPropagationDirection(IntegrationSense.DECREASING);
    propagateBackwardToSoi.setStopOnMaximumDurationBehavior(MaximumDurationBehavior.THROW);

    // Stop propagating this segment at the Lunar SOI boundary.
    final 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.getIntegrationPoint(),
                                                                                                          CartesianElement.MAGNITUDE,
                                                                                                          moon.getInertialFrame()),
                                                                                   moonSphereOfInfluence,
                                                                                   0.0001, // Value tolerance, meters.
                                                                                   StopType.ANY_THRESHOLD);
    moonSoiStoppingCondition.setName("Moon_SOI");

    propagateBackwardToSoi.getStoppingConditions().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.getInertialFrame(), Cartesian.getUndefined(), Cartesian.getUndefined());

    // Total mass.
    propagationPointAroundEarth.setMass(dryMass.getIntegrationValue());

    // Gravitational field for Earth.
    SphericalHarmonicGravityField earthGravityField = new SphericalHarmonicGravityField(earthGravityModel, 4, 4, true, new PermanentSolidTideModel());
    SphericalHarmonicGravity earthGravityAsPrimary = new SphericalHarmonicGravity(propagationPointAroundEarth.getIntegrationPoint(), earthGravityField);
    propagationPointAroundEarth.getAppliedForces().add(earthGravityAsPrimary);

    // The gravity of the Sun and Moon are modeled as third bodies for the Earth-centered propagator.
    ThirdBodyGravity thirdBodyForEarthPropagator = new ThirdBodyGravity();
    thirdBodyForEarthPropagator.setTargetPoint(propagationPointAroundEarth.getIntegrationPoint());
    thirdBodyForEarthPropagator.addThirdBody(moon.getName(), moon.getCenterOfMassPoint(), moonGravitationalParameter);
    thirdBodyForEarthPropagator.addThirdBody(sun.getName(), sun.getCenterOfMassPoint(), sunGravitationalParameter);
    propagationPointAroundEarth.getAppliedForces().add(thirdBodyForEarthPropagator);

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

    // Add the elements.
    earthCenteredNumericalPropagator.getIntegrationElements().add(propagationPointAroundEarth);
    earthCenteredNumericalPropagator.getIntegrationElements().add(dryMass);

    // Set the integrator.
    earthCenteredNumericalPropagator.setIntegrator(integrator);

    // Set epoch.
    earthCenteredNumericalPropagator.setEpoch(epoch);

    // Make a segment that propagates backward two days using Earth-centered propagator.
    NumericalPropagatorSegment propagateBackwardTwoDays = new NumericalPropagatorSegment();
    propagateBackwardTwoDays.setName("Propagate_Backward_TwoDays");
    propagateBackwardTwoDays.setPropagatorDefinition(earthCenteredNumericalPropagator);
    propagateBackwardTwoDays.setPropagationDirection(IntegrationSense.DECREASING);
    propagateBackwardTwoDays.setStopOnMaximumDurationBehavior(MaximumDurationBehavior.THROW);
    propagateBackwardTwoDays.getStoppingConditions().add(new DurationStoppingCondition(Duration.fromDays(-2.0)));

    // Make a segment that propagates backward to perigee.
    NumericalPropagatorSegment propagateBackwardToPerigee = new NumericalPropagatorSegment();
    propagateBackwardToPerigee.setName("Propagate_Backward_ToPerigee");
    propagateBackwardToPerigee.setPropagatorDefinition(earthCenteredNumericalPropagator);
    propagateBackwardToPerigee.setPropagationDirection(IntegrationSense.DECREASING);
    propagateBackwardToPerigee.setStopOnMaximumDurationBehavior(MaximumDurationBehavior.THROW);

    ScalarStoppingCondition perigeeStoppingCondition1 = new ScalarStoppingCondition(new ScalarPointElement(propagationPointAroundEarth.getIntegrationPoint(),
                                                                                                           CartesianElement.MAGNITUDE,
                                                                                                           earth.getInertialFrame()),
                                                                                    0.0001, // Value tolerance, meters.
                                                                                    StopType.LOCAL_MINIMUM);
    perigeeStoppingCondition1.setName("Earth_Perigee_Backwards_Targeting");

    propagateBackwardToPerigee.getStoppingConditions().add(perigeeStoppingCondition1);

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

    backwardSegmentList.getSegments().add(initialStateSegment);
    backwardSegmentList.getSegments().add(propagateBackwardToSoi);
    backwardSegmentList.getSegments().add(propagateBackwardTwoDays);
    backwardSegmentList.getSegments().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 = SetVariableCallback.of((variable, configuration) ->  {
            Motion1<Cartesian> currentMotion = configuration.getMotion(motionID);
            Motion1<Cartesian> rotatingMotion = toMoonRotatingAtEpoch.transform(currentMotion, 1);
            double updatedVy = rotatingMotion.getFirstDerivative().getY() + variable;
            Motion1<Cartesian> updatedRotating = new Motion1<Cartesian>(rotatingMotion.getValue(),
                                                                        new Cartesian(rotatingMotion.getFirstDerivative().getX(), updatedVy, rotatingMotion.getFirstDerivative().getZ()));
            Motion1<Cartesian> 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 = SetVariableCallback.of((variable, configuration) -> {
            Motion1<Cartesian> currentMotion = configuration.getMotion(motionID);
            Motion1<Cartesian> rotatingMotion = toMoonRotatingAtEpoch.transform(currentMotion, 1);
            double updatedVz = rotatingMotion.getFirstDerivative().getZ() + variable;
            Motion1<Cartesian> updatedRotating = new Motion1<Cartesian>(rotatingMotion.getValue(),
                                                                        new Cartesian(rotatingMotion.getFirstDerivative().getX(), rotatingMotion.getFirstDerivative().getY(), updatedVz));
            Motion1<Cartesian> 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.getIntegrationPoint(), 
                                                                      CartesianElement.MAGNITUDE,
                                                                      earth.getInertialFrame());
    ScalarAtEndOfNumericalSegmentConstraint backwardPerigeeRadiusConstraint = new ScalarAtEndOfNumericalSegmentConstraint(backwardPerigeeRadius,
                                                                                                                          propagateBackwardToPerigee,
                                                                                                                          earthPhysicalRadius + 241000.0, // Desired value, meters.
                                                                                                                          0.001, // Tolerance, meters.
                                                                                                                          EndsOfSegment.FINAL);

    ScalarModifiedKeplerianElement backwardInclination = new ScalarModifiedKeplerianElement(earthGravitationalParameter, 
                                                                                            propagationPointAroundEarth.getIntegrationPoint(),
                                                                                            KeplerianElement.INCLINATION,
                                                                                            earth.getInertialFrame());
    ScalarAtEndOfNumericalSegmentConstraint backwardInclinationConstraint = new ScalarAtEndOfNumericalSegmentConstraint(backwardInclination,
                                                                                                                        propagateBackwardToPerigee,
                                                                                                                        Trig.degreesToRadians(28.5),
                                                                                                                        Trig.degreesToRadians(0.1),
                                                                                                                        EndsOfSegment.FINAL);

    TargetedSegmentListDifferentialCorrector correctorBackwardTrajectory = new TargetedSegmentListDifferentialCorrector();
    correctorBackwardTrajectory.setName("Correct_Backward_Trajectory");
    correctorBackwardTrajectory.setMaximumIterations(75);
    correctorBackwardTrajectory.setOnlyStoreFinalResults(false);
    correctorBackwardTrajectory.getSolver().setMultithreaded(false);

    correctorBackwardTrajectory.getConstraints().add(backwardPerigeeRadiusConstraint);
    correctorBackwardTrajectory.getConstraints().add(backwardInclinationConstraint);

    correctorBackwardTrajectory.getVariables().add(velocityYVariable);
    correctorBackwardTrajectory.getVariables().add(velocityZVariable);

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

    backwardSegmentList.setOperatorAction(TargetedSegmentListOperatorBehavior.RUN_ACTIVE_OPERATORS);

    // Run backwardSegmentList.
    SegmentPropagator propagator = backwardSegmentList.getSegmentPropagator(new EvaluatorGroup());
    TargetedSegmentListResults listResults = (TargetedSegmentListResults)propagator.propagate();

    DateMotionCollection1<Cartesian> barycentricRotatingEphemeris = listResults.getDateMotionCollectionOfOverallTrajectory(motionID, barycentricRotatingFrame);

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

    // Compute translunar injection and parking orbit parameters.
    Motion1<Cartesian> translunarInjectionState = earthCenteredInertialEphemeris.getMotions().get(0);
    ModifiedKeplerianElements translunarKeplerianElements = new ModifiedKeplerianElements(translunarInjectionState, earthGravitationalParameter);

    double translunarInjectionVelocity = translunarInjectionState.getFirstDerivative().getMagnitude();
    double parkingOrbitVelocity = Math.sqrt(earthGravitationalParameter / translunarInjectionState.getValue().getMagnitude());
    double translunarInjectionDeltaV = translunarInjectionVelocity - parkingOrbitVelocity;

    Cartesian parkingOrbitVelocityCartesian = Cartesian.multiply(translunarInjectionState.getFirstDerivative(), (parkingOrbitVelocity / translunarInjectionVelocity));

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

    // 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.setName("Forward_Initial_State_Segment");

    // Use cloned version of the Earth centered backward propagator that has modified initial state and epoch.
    NumericalPropagatorDefinition forwardEarthCenteredPropagator = (NumericalPropagatorDefinition) earthCenteredNumericalPropagator.clone(new CopyContext());
    forwardEarthCenteredPropagator.setEpoch(earthCenteredInertialEphemeris.getDates().get(0));
    PropagationNewtonianPoint forwardEarthCenteredPoint = (PropagationNewtonianPoint)((DefinitionalObjectCollection<PropagationStateElement>)forwardEarthCenteredPropagator.getIntegrationElements()).get(0);

    forwardEarthCenteredPoint.setInitialPosition(parkingOrbitState.getValue());
    forwardEarthCenteredPoint.setInitialVelocity(parkingOrbitState.getFirstDerivative());

    forwardInitialStateSegment.setPropagatorDefinition(forwardEarthCenteredPropagator);

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

    // The maneuver.
    ImpulsiveManeuverInformation transLunarBurnDetails = new ImpulsiveManeuverInformation(motionID, new Cartesian(translunarInjectionDeltaV, 0.0, 0.0));
    transLunarBurnDetails.setOrientation(new AxesVelocityOrbitNormal(transLunarBurnDetails.getPropagationPoint(), earth));
    transLunarInjectionSegment.setManeuver(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 = SetVariableCallback.of((variable, configuration) -> {
            Motion1<Cartesian> currentMotion = configuration.getMotion(motionID);
            ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter);
            ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(oldElements.getRadiusOfPeriapsis(), 
                                                                                  oldElements.getInverseSemimajorAxis(),
                                                                                  oldElements.getInclination(),
                                                                                  oldElements.getArgumentOfPeriapsis(),
                                                                                  oldElements.getRightAscensionOfAscendingNode() + variable,
                                                                                  oldElements.getTrueAnomaly(),
                                                                                  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.setName("InitialState_RightAscension");

    // True anomaly variable. (This is really argument of latitude since the orbit is circular and inclined.) 
    SetVariableCallback<NumericalPropagatorState> setterTrueAnomaly = SetVariableCallback.of((variable, configuration) -> {
            Motion1<Cartesian> currentMotion = configuration.getMotion(motionID);
            ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter);
            ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(oldElements.getRadiusOfPeriapsis(), 
                                                                                  oldElements.getInverseSemimajorAxis(),
                                                                                  oldElements.getInclination(),
                                                                                  oldElements.getArgumentOfPeriapsis(),
                                                                                  oldElements.getRightAscensionOfAscendingNode(),
                                                                                  oldElements.getTrueAnomaly() + 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.setName("InitialState_TrueAnomaly");

    // Maneuver delta-v variable.
    SegmentPropagatorVariable translunarInjectionDeltaVxVariable = transLunarInjectionSegment.createVariable(100.0, // Maximum step, m/s.
                                                                                                             1.0, // Perturbation, m/s.
                                                                                                             SetVariableCallback.of((variable, configuration) -> {
            configuration.getManeuver().setX(configuration.getManeuver().getX() + variable);
    }));
    translunarInjectionDeltaVxVariable.setName("TranslunarInjection_ThrustVector_X");

    // Make a segment that propagates forward to edge of Lunar SOI.
    NumericalPropagatorSegment propagateForwardToSoi1 = new NumericalPropagatorSegment();
    propagateForwardToSoi1.setName("Propagate_Forward_ToSOI_BeforeFlyby");
    propagateForwardToSoi1.setPropagatorDefinition(forwardEarthCenteredPropagator);
    propagateForwardToSoi1.setPropagationDirection(IntegrationSense.INCREASING);
    propagateForwardToSoi1.setStopOnMaximumDurationBehavior(MaximumDurationBehavior.THROW);

    // Stopping condition can be reused here.
    propagateForwardToSoi1.getStoppingConditions().add(moonSoiStoppingCondition);

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

    // No need to change initial conditions here because the state and time will be overriden.
    propagateForwardToPerilune.setPropagatorDefinition(nearMoonNumericalPropagator);
    propagateForwardToPerilune.setPropagationDirection(IntegrationSense.INCREASING);
    propagateForwardToPerilune.setStopOnMaximumDurationBehavior(MaximumDurationBehavior.THROW);

    ScalarStoppingCondition periluneStoppingCondition = new ScalarStoppingCondition(new ScalarPointElement(propagationPointAroundMoon.getIntegrationPoint(),
                                                                                                           CartesianElement.MAGNITUDE,
                                                                                                           moon.getInertialFrame()),
                                                                                    0.0001, // Value tolerance, meters.
                                                                                    StopType.LOCAL_MINIMUM);
    periluneStoppingCondition.setName("Moon_Perilune");

    propagateForwardToPerilune.getStoppingConditions().add(periluneStoppingCondition);

    // Make a segment that propagates forward to edge of Lunar SOI.
    NumericalPropagatorSegment propagateForwardToSoi2 = new NumericalPropagatorSegment();
    propagateForwardToSoi2.setName("Propagate_Forward_ToSOI_AfterFlyby");
    propagateForwardToSoi2.setPropagatorDefinition(nearMoonNumericalPropagator);
    propagateForwardToSoi2.setPropagationDirection(IntegrationSense.INCREASING);
    propagateForwardToSoi2.setStopOnMaximumDurationBehavior(MaximumDurationBehavior.THROW);

    // Stopping condition can be reused here.
    propagateForwardToSoi2.getStoppingConditions().add(moonSoiStoppingCondition);

    // Make a segment that propagates forward two days.
    NumericalPropagatorSegment propagateForwardTwoDays = new NumericalPropagatorSegment();
    propagateForwardTwoDays.setName("Propagate_Forward_Two_Days");
    propagateForwardTwoDays.setPropagatorDefinition(forwardEarthCenteredPropagator);
    propagateForwardTwoDays.getStoppingConditions().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.setName("Propagate_Backward_ToPerigee");
    propagateForwardToPerigee.setPropagatorDefinition(forwardEarthCenteredPropagator);
    propagateForwardToPerigee.setPropagationDirection(IntegrationSense.INCREASING);
    propagateForwardToPerigee.setStopOnMaximumDurationBehavior(MaximumDurationBehavior.THROW);

    ScalarStoppingCondition perigeeStoppingCondition2 = new ScalarStoppingCondition(new ScalarPointElement(forwardEarthCenteredPoint.getIntegrationPoint(),
                                                                                                           CartesianElement.MAGNITUDE,
                                                                                                           earth.getInertialFrame()),
                                                                                    0.0001, // Value tolerance, meters.
                                                                                    StopType.LOCAL_MINIMUM);
    perigeeStoppingCondition2.setName("Earth_Perigee_After_Return");

    ScalarStoppingCondition altitudeStoppingCondition = new ScalarStoppingCondition(new ScalarPointElement(forwardEarthCenteredPoint.getIntegrationPoint(),
                                                                                                           CartesianElement.MAGNITUDE,
                                                                                                           earth.getInertialFrame()), 
                                                                                    earthPhysicalRadius + 122500.0, // Threshold, meters.
                                                                                    0.0001, // Value tolerance, meters.
                                                                                    StopType.ANY_THRESHOLD);
    altitudeStoppingCondition.setName("Earth_Entry_Altitude");

    propagateForwardToPerigee.getStoppingConditions().add(perigeeStoppingCondition2);
    propagateForwardToPerigee.getStoppingConditions().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.getIntegrationPoint(), CartesianElement.MAGNITUDE, earth.getInertialFrame());
    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.getZonalCoefficients()[2];
    earthJ2 = earthGravityModel.getNormalized() ? earthJ2 * Math.sqrt(5.0) : earthJ2;

    ScalarKozaiIzsakMeanElement radialVelocityScalar = new ScalarKozaiIzsakMeanElement(earthGravitationalParameter,
                                                                                       forwardEarthCenteredPoint.getIntegrationPoint(),
                                                                                       KozaiIzsakMeanElement.RADIAL_VELOCITY,
                                                                                       earth.getInertialFrame(),
                                                                                       earthJ2,
                                                                                       earthPhysicalRadius);

    ScalarVectorElement velocityMagnitudeScalar = new ScalarVectorElement(forwardEarthCenteredPoint.getIntegrationPoint().createVectorVelocity(earth.getInertialFrame()),
                                                                        CartesianElement.MAGNITUDE,
                                                                        earth.getInertialFrame().getAxes());

    Scalar sineFlightPathAngleScalar = Scalar.divide(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.setName("Forward_Segment_List");

    forwardSegmentList.getSegments().add(forwardInitialStateSegment);
    forwardSegmentList.getSegments().add(transLunarInjectionSegment);
    forwardSegmentList.getSegments().add(propagateForwardToSoi1);
    forwardSegmentList.getSegments().add(propagateForwardToPerilune);
    forwardSegmentList.getSegments().add(propagateForwardToSoi2);
    forwardSegmentList.getSegments().add(propagateForwardTwoDays);
    forwardSegmentList.getSegments().add(propagateForwardToPerigee);

    TargetedSegmentListDifferentialCorrector altitudeCorrectorForwardTrajectory = new TargetedSegmentListDifferentialCorrector();

    altitudeCorrectorForwardTrajectory.setName("Correct_Forward_Trajectory_Altitude");
    altitudeCorrectorForwardTrajectory.getVariables().add(rightAscensionVariable);
    altitudeCorrectorForwardTrajectory.getVariables().add(trueAnomalyVariable);
    altitudeCorrectorForwardTrajectory.getVariables().add(translunarInjectionDeltaVxVariable);

    altitudeCorrectorForwardTrajectory.getConstraints().add(altitudeConstraint);

    altitudeCorrectorForwardTrajectory.setMaximumIterations(25);
    altitudeCorrectorForwardTrajectory.setOnlyStoreFinalResults(false);

    forwardSegmentList.getOperators().add(altitudeCorrectorForwardTrajectory);

    TargetedSegmentListDifferentialCorrector flightPathCorrectorForwardTrajectory = new TargetedSegmentListDifferentialCorrector();

    flightPathCorrectorForwardTrajectory.setName("Correct_Forward_Trajectory_FlightPath");
    flightPathCorrectorForwardTrajectory.getVariables().add(rightAscensionVariable);
    flightPathCorrectorForwardTrajectory.getVariables().add(trueAnomalyVariable);
    flightPathCorrectorForwardTrajectory.getVariables().add(translunarInjectionDeltaVxVariable);

    flightPathCorrectorForwardTrajectory.getConstraints().add(altitudeConstraint);
    flightPathCorrectorForwardTrajectory.getConstraints().add(flightPathAngleConstraint);

    flightPathCorrectorForwardTrajectory.setMaximumIterations(25);
    flightPathCorrectorForwardTrajectory.setOnlyStoreFinalResults(false);

    forwardSegmentList.getOperators().add(flightPathCorrectorForwardTrajectory);

    // Run forwardSegmentList.
    SegmentPropagator forwardPropagator = forwardSegmentList.getSegmentPropagator(new EvaluatorGroup());
    TargetedSegmentListResults forwardListResults = (TargetedSegmentListResults) forwardPropagator.propagate();

    DateMotionCollection1<Cartesian> forwardEarthCenteredInertialEphemeris = forwardListResults.getDateMotionCollectionOfOverallTrajectory(motionID, earth.getInertialFrame());

    DateMotionCollection1<Cartesian> forwardEarthFixedEphemeris = forwardListResults.getDateMotionCollectionOfOverallTrajectory(motionID, earth.getFixedFrame());

    // 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.
    Motion1<Cartesian> updatedParkingOrbitState = forwardEarthCenteredInertialEphemeris.getMotions().get(0);
    ModifiedKeplerianElements updatedParkingOrbitElements = new ModifiedKeplerianElements(updatedParkingOrbitState, earthGravitationalParameter);

    JulianDate parkingOrbitEpoch = forwardEarthCenteredInertialEphemeris.getDates().get(0);
    double parkingOrbitSemimajorAxis = 1.0 / updatedParkingOrbitElements.getInverseSemimajorAxis();
    double parkingOrbitEccentricity = updatedParkingOrbitElements.computeEccentricity();
    double parkingOrbitInclination = Trig.radiansToDegrees(updatedParkingOrbitElements.getInclination());
    double parkingOrbitArgumentPerigee = Trig.radiansToDegrees(updatedParkingOrbitElements.getArgumentOfPeriapsis());
    double updatedRaan = Trig.radiansToDegrees(updatedParkingOrbitElements.getRightAscensionOfAscendingNode());
    double updatedTrueAnomaly = Trig.radiansToDegrees(updatedParkingOrbitElements.getTrueAnomaly());

    Motion1<Cartesian> updatedTranslunarInjectionState = forwardEarthCenteredInertialEphemeris.getMotions().get(1);
    Cartesian updatedTranslunarInjectionDeltaV = Cartesian.subtract(updatedTranslunarInjectionState.getFirstDerivative(), updatedParkingOrbitState.getFirstDerivative());
    double updatedDeltaVMagnitude = updatedTranslunarInjectionDeltaV.getMagnitude();

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

    if (propagateToPeriluneResults != null) {
        JulianDate periluneDate = propagateToPeriluneResults.getFinalPropagatedState().getCurrentDate();

        Motion1<Cartesian> periluneState = propagateToPeriluneResults.getFinalPropagatedState().<Cartesian> getMotion(motionID);
        ModifiedKeplerianElements periluneElements = new ModifiedKeplerianElements(periluneState, moonGravitationalParameter);

        double lunarFlybyAltitude = periluneElements.getRadiusOfPeriapsis() - moonPhysicalRadius;
    }

    // Get Earth re-entry state and elements.
    int forwardCount = forwardEarthCenteredInertialEphemeris.getCount();
    Motion1<Cartesian> earthReentryState = forwardEarthCenteredInertialEphemeris.getMotions().get(forwardCount - 1);
    Cartesian earthFixedPosition = forwardEarthFixedEphemeris.getValues().get(forwardCount - 1);
    ModifiedKeplerianElements earthReentryElements = new ModifiedKeplerianElements(earthReentryState, earthGravitationalParameter);

    JulianDate earthReentryDate = forwardEarthCenteredInertialEphemeris.getDates().get(forwardCount - 1);

    Cartographic earthReentryCartographic = earth.getShape().cartesianToCartographic(earthFixedPosition);

    double altitude = earthReentryCartographic.getHeight();
    double latitude = Trig.radiansToDegrees(earthReentryCartographic.getLatitude());
    double longitude = Trig.radiansToDegrees(earthReentryCartographic.getLongitude());

    double eccentricity = earthReentryElements.computeEccentricity();
    double trueAnomaly = earthReentryElements.getTrueAnomaly();
    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));

Copyright © 2007-2023 ANSYS, Inc. and/or its affiliated companies. All rights reserved. Unauthorized use, distribution, or reproduction is prohibited.

Send comments on this topic to support@agi.com