Click or drag to resize

Code Sample

This topic presents a complete code example using Segment Propagation Library to compute a satellite's trajectory to the Moon. We state the problem as follows:

Calculate a trajectory to the Moon starting around 1 Aug 2009 15:00:00.000. Begin in a low Earth orbit. The first maneuver should have a delta-v of 4.0 km/s, but the delta-v for the capture maneuver around the moon is to be determined. Propagate around the Moon in a polar orbit for 3 days.

Note Note

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

Designing Segments

To start with, we break up our trajectory into conceptually separate segments and think about what we know about our problem and what we need to solve for:

  • Initial State - Must determine best starting values

  • Parking Orbit - Stop after an hour

  • Impulsive Maneuver into Transfer Orbit

  • Transfer orbit away from the Earth - Stop at 320,000 km

  • Continue the transfer orbit in the Moon reference frame - Ideally stop at perilune

  • Impulsive Maneuver into a circular lunar orbit

  • Propagate for 3 days around the Moon

We will start with defining several constants and types that will be used throughout the example:

Java
Copy
// pick names
final String motionID = "Satellite";
final String fuelMassName = "Fuel_Mass";
final String dryMassName = "Dry_Mass";

// store commonly used central bodies and constants
CentralBody moon = CentralBodiesFacet.getFromContext().getMoon();
double moonGravitationalParameter = 4902801076000.0;
CentralBody earth = CentralBodiesFacet.getFromContext().getEarth();
final double earthGravitationalParameter = WorldGeodeticSystem1984.GravitationalParameter;

// initial date
JulianDate epoch = new JulianDate(new GregorianDate(2009, 7, 1, 15, 0, 0.0), TimeStandard.getCoordinatedUniversalTime());

// initial state
Motion1<Cartesian> initialState =
        new ModifiedKeplerianElements(6674314.0,
                                      1.0 / 6810520.0,
                                      Trig.degreesToRadians(15.0),
                                      Trig.degreesToRadians(300.0),
                                      Trig.degreesToRadians(290.0),
                                      Trig.degreesToRadians(1.808),
                                      earthGravitationalParameter).toCartesian();

double exhaustVelocity = 4500.0; // meters/second

We then create several of the items for the NumericalPropagatorDefinition that we will be using:

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

// dry mass
PropagationScalar dryMass = new PropagationScalar(100.0);
dryMass.setIdentification(dryMassName);
dryMass.setScalarDerivative(new ScalarFixed(0.0));

// Create the integrator
RungeKutta4Integrator integrator = new RungeKutta4Integrator();
integrator.setInitialStepSize(120.0);

Next, we create the actual NumericalPropagatorDefinition that we will use for our Earth-centered segments:

Java
Copy
// the propagation point
PropagationNewtonianPoint propagationPointAroundEarth =
        new PropagationNewtonianPoint(motionID,
                                      earth.getInertialFrame(),
                                      initialState.getValue(),
                                      initialState.getFirstDerivative());

// total mass
propagationPointAroundEarth.setMass(Scalar.add(fuelMass.getIntegrationValue(), dryMass.getIntegrationValue()));

// simple gravity for Earth
TwoBodyGravity earthGravityAsPrimary =
        new TwoBodyGravity(propagationPointAroundEarth.getIntegrationPoint(),
                           earth,
                           earthGravitationalParameter);
propagationPointAroundEarth.getAppliedForces().add(earthGravityAsPrimary);

// and the Moons gravity
ThirdBodyGravity moonGravityAsThirdBody = new ThirdBodyGravity();
moonGravityAsThirdBody.setTargetPoint(propagationPointAroundEarth.getIntegrationPoint());
moonGravityAsThirdBody.addThirdBody(moon.getName(), moon.getCenterOfMassPoint(), moonGravitationalParameter);
propagationPointAroundEarth.getAppliedForces().add(moonGravityAsThirdBody);

// make the actual propagator
NumericalPropagatorDefinition earthCenteredNumericalPropagator = new NumericalPropagatorDefinition();

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

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

// set epoch
earthCenteredNumericalPropagator.setEpoch(epoch);
Initial Segments

We need to specify a whole set of orbital elements for the initial state. However, because we will need to have a TargetedSegmentList to edit some of those initial values, we cannot just start propagating forward with a NumericalPropagatorSegment. Only the InitialStateSegment<T> allows initial state values to be modified in a TargetedSegmentList:

Java
Copy
// Simply the initial state.  We will need it since the differential corrector will operate
// on the initial values.
NumericalInitialStateSegment initialStateSegment = new NumericalInitialStateSegment();
initialStateSegment.setName("Initial_State_Segment");
initialStateSegment.setPropagatorDefinition(earthCenteredNumericalPropagator);

Next, we propagate for a span of time. As long as we are propagating the same elements and the force models won't change, we can keep on using the same NumericalPropagatorDefinition.

Java
Copy
// configure the segment
NumericalPropagatorSegment parkingOrbitSegment = new NumericalPropagatorSegment();
parkingOrbitSegment.setName("Parking_Orbit");
parkingOrbitSegment.setPropagatorDefinition(earthCenteredNumericalPropagator);

// propagate for almost an hour
DurationStoppingCondition timeOnParkingOrbit = new DurationStoppingCondition(Duration.fromSeconds(3200.0));
timeOnParkingOrbit.setName("Duration_On_Parking_Orbit");
parkingOrbitSegment.getStoppingConditions().add(timeOnParkingOrbit);
Lunar Transfer

Next, we create our first maneuver getting us onto a lunar transfer orbit. We want the delta-v in the velocity-orbit-normal axes.

Java
Copy
ImpulsiveManeuverSegment transLunarInjectionSegment = new ImpulsiveManeuverSegment();
transLunarInjectionSegment.setName("Trans-Lunar_Injection");

// the maneuver
ImpulsiveManeuverInformation transLunarBurnDetails =
        new ImpulsiveManeuverInformation(motionID,
                                         new Cartesian(4000.0, 0.0, 0.0), // 4000 meters/second is actually more than we need
                                         fuelMassName,
                                         dryMassName,
                                         Scalar.toScalar(exhaustVelocity),
                                         InvalidFuelStateBehavior.THROW_EXCEPTION);
transLunarBurnDetails.setOrientation(new AxesVelocityOrbitNormal(transLunarBurnDetails.getPropagationPoint(), earth));
transLunarInjectionSegment.setManeuver(transLunarBurnDetails);

Now we propagate away from the Earth. Again, we use the same NumericalPropagatorDefinition as we used in the previous segments. Remember that at this point, if we just propagate, we have no idea where our satellite will end up. But while defining the problem, we can assume we are on the correct trajectory for now.

Deciding when to stop is an interesting question. We could just propagate all the way to the Moon, but having an Earth-centered integration frame for a lunar satellite doesn't make much sense. Our propagator has both the Earth and Moon as gravity, so we don't need to stop exactly on the sphere of influence between the two. Let's pick 320,000 kilometers from the Earth as our cut-off point (close to the Moon's sphere of influence radius). Before we reach that distance, we are leaving the Earth. After that distance, we are approaching the Moon.

Java
Copy
// propagate until we are inside the moons sphere of influence
NumericalPropagatorSegment propagateOutOfEarthsSoiSegment = new NumericalPropagatorSegment();
propagateOutOfEarthsSoiSegment.setName("Propagate_To_Moons_SOI");
propagateOutOfEarthsSoiSegment.setPropagatorDefinition(earthCenteredNumericalPropagator);

// stop propagating this segment at the SOI boundary
ScalarStoppingCondition awayFromEarthStoppingCondition =
        new ScalarStoppingCondition(new ScalarPointElement(propagationPointAroundEarth.getIntegrationPoint(),
                                                           CartesianElement.MAGNITUDE,
                                                           earth.getInertialFrame()),
                320000000.0, // approximate SOI for moon, meters
                0.0001, // value tolerance, meters
                StopType.ANY_THRESHOLD);
awayFromEarthStoppingCondition.setName("Moon_SOI");

propagateOutOfEarthsSoiSegment.getStoppingConditions().add(awayFromEarthStoppingCondition);

Since we are changing the force models, we need to create a new NumericalPropagatorDefinition. Notice that we are using the same names for all of the propagation elements.

Java
Copy
PropagationNewtonianPoint propagationPointAroundMoon = new PropagationNewtonianPoint();
propagationPointAroundMoon.setIdentification(motionID);

// NOTE: That because this segment is not the first segment, its initial values
// will be initialized with the final state of the previous segment.  So it
// does not matter what values we put in here; they will get updated at
// propagation time.
propagationPointAroundMoon.setInitialPosition(initialState.getValue());
propagationPointAroundMoon.setInitialVelocity(initialState.getFirstDerivative());

// This is the true reason why we are changing propagators, so that we are
// producing ephemeris in the Moon's frame instead of Earth's.
propagationPointAroundMoon.setIntegrationFrame(moon.getInertialFrame());

// Configure fuel.  Again, the values will get initialized from the previous
// segment's final state.
propagationPointAroundMoon.setMass(Scalar.add(dryMass.getIntegrationValue(), fuelMass.getIntegrationValue()));

// Moon's primary gravity
TwoBodyGravity moonPrimaryGravity =
        new TwoBodyGravity(propagationPointAroundMoon.getIntegrationPoint(),
                           moon,
                           moonGravitationalParameter);

propagationPointAroundMoon.getAppliedForces().add(moonPrimaryGravity);

// Earth gravity as third body gravity
ThirdBodyGravity earthGravityAsThirdBody = new ThirdBodyGravity(propagationPointAroundMoon.getIntegrationPoint());
earthGravityAsThirdBody.addThirdBody(earth.getName(), earth.getCenterOfMassPoint(), earthGravitationalParameter);
earthGravityAsThirdBody.setCentralBody(moon);

propagationPointAroundMoon.getAppliedForces().add(earthGravityAsThirdBody);

// make the actual propagator
NumericalPropagatorDefinition nearMoonNumericalPropagator = new NumericalPropagatorDefinition();

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

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

// set an epoch
nearMoonNumericalPropagator.setEpoch(epoch);
Lunar Orbit

Finally, we are approaching the Moon. The transformation of the motion from the Earth's inertial frame to the Moon's inertial frame is handled automatically. Ideally we want to stop at perilune, and we want perilune to be along a polar orbit. However, as we solve for our perfect orbit the intermediate guesses may be very far off. As such, we include several other stopping conditions to ensure the extreme early guesses don't go forever. We will add an altitude stopping condition over the moon to avoid disturbing the lunar environment, and a two day duration so we don't propagate forever if we don't get anywhere close to the Moon.

Java
Copy
NumericalPropagatorSegment transferOrbitWithinMoonsSoiSegment = new NumericalPropagatorSegment();
transferOrbitWithinMoonsSoiSegment.setName("Transfer_Orbit_Approaching_Moon");

// set the NumericalPropagator on the propagate segment
transferOrbitWithinMoonsSoiSegment.setPropagatorDefinition(nearMoonNumericalPropagator);

// lunar altitude stopping condition to avoid a rapid unplanned descent
// The tolerances are more relaxed because this condition is not the intended
// one to stop on and we don't need to be exact in such a case.
ScalarStoppingCondition altitudeOfMoonStoppingCondition =
        new ScalarStoppingCondition(new ScalarCartographicElement(moon,
                                                                  propagationPointAroundMoon.getIntegrationPoint(),
                                                                  CartographicElement.HEIGHT),
                1000.0, // threshold, meters
                0.1, // value tolerance, meters
                StopType.ANY_THRESHOLD);
altitudeOfMoonStoppingCondition.setName("Lunar_Altitude");

transferOrbitWithinMoonsSoiSegment.getStoppingConditions().add(altitudeOfMoonStoppingCondition);

// periapsis stopping condition
ScalarStoppingCondition periluneStoppingCondition =
        new ScalarStoppingCondition(new ScalarModifiedKeplerianElement(moonGravitationalParameter,
                                                                       propagationPointAroundMoon.getIntegrationPoint(),
                                                                       KeplerianElement.TRUE_ANOMALY,
                                                                       moon.getInertialFrame()),
                0.0, // threshold, radians
                0.000001, // tolerance, radians
                StopType.ANY_THRESHOLD);
// leave the threshold 0.0 radians
periluneStoppingCondition.setAngularSetting(CircularRange.NEGATIVE_PI_TO_PI);
periluneStoppingCondition.setName("Lunar_Periapsis");

transferOrbitWithinMoonsSoiSegment.getStoppingConditions().add(periluneStoppingCondition);

// and stop on Duration just in case we completely miss
DurationStoppingCondition durationStoppingConditionWhenApproachingMoon =
        new DurationStoppingCondition(Duration.fromDays(2.0));
transferOrbitWithinMoonsSoiSegment.getStoppingConditions().add(durationStoppingConditionWhenApproachingMoon);

Notice that because we are defining our point's propagation element in a different ReferenceFrame, the ephemeris for this point element in the EphemerisForOverallTrajectory (get) will be in the Moon's inertial frame. For the overall SegmentListResults, its entire ephemeris will be in the frame of the final segment propagated, which in this case is the Moon's inertial frame.

Ignoring that we need to perform some form of targeting to even get this far, we continue designing our segments as if we were on our ideal orbit. We perform a maneuver to capture our satellite around the moon:

Java
Copy
// to get captured around the moon
ImpulsiveManeuverSegment lunarOrbitInsertionManeuver = new ImpulsiveManeuverSegment();

lunarOrbitInsertionManeuver.setName("Lunar_Orbit_Insertion");

// the maneuver details
ImpulsiveManeuverInformation lunarOrbitInsertionManeuverDetails =
        new ImpulsiveManeuverInformation(motionID,
                                         new Cartesian(-1000.0, 0.0, 0.0),
                                         fuelMassName,
                                         dryMassName,
                                         Scalar.toScalar(exhaustVelocity),
                                         InvalidFuelStateBehavior.THROW_EXCEPTION);

lunarOrbitInsertionManeuverDetails.setOrientation(
        new AxesVelocityOrbitNormal(lunarOrbitInsertionManeuverDetails.getPropagationPoint(), moon));

lunarOrbitInsertionManeuver.setManeuver(lunarOrbitInsertionManeuverDetails);

Finally, we propagate around the moon for three days:

Java
Copy
// final segment for this example, we made it to the Moon!
NumericalPropagatorSegment aroundTheMoonSegment = new NumericalPropagatorSegment();
aroundTheMoonSegment.setName("Lunar_Orbit");

// same propagator
aroundTheMoonSegment.setPropagatorDefinition(nearMoonNumericalPropagator);

// propagate and do science
aroundTheMoonSegment.getStoppingConditions().add(new DurationStoppingCondition(Duration.fromDays(3.0)));
Segment Lists

Before diving into the targeting problem, we create our overall SegmentList and the TargetedSegmentLists that will organize and run all of the individual segments we have made:

Java
Copy
// make the main list
SegmentList mainSegmentList = new SegmentList();
mainSegmentList.setName("Main_Segment_List");

// we want to target items in the initial state, so we need a TargetedSegmentList
TargetedSegmentList targetedSegmentToGetToTheMoon = new TargetedSegmentList();
targetedSegmentToGetToTheMoon.setName("Targeted_Segment_To_Get_To_Moon");

targetedSegmentToGetToTheMoon.getSegments().add(initialStateSegment);
targetedSegmentToGetToTheMoon.getSegments().add(parkingOrbitSegment);
targetedSegmentToGetToTheMoon.getSegments().add(transLunarInjectionSegment);
targetedSegmentToGetToTheMoon.getSegments().add(propagateOutOfEarthsSoiSegment);
targetedSegmentToGetToTheMoon.getSegments().add(transferOrbitWithinMoonsSoiSegment);

// add that TargetedSegmentList to the overall segment list
mainSegmentList.getSegments().add(targetedSegmentToGetToTheMoon);

// since circularizing the final orbit is done after the previous TargetedSegmentList
// is done, we can put it into its own TargetedSegmentList.
TargetedSegmentList circularizeLunarOrbitTargetedSegment = new TargetedSegmentList();
circularizeLunarOrbitTargetedSegment.setName("Lunar_Orbit_Insertion_Targeted_Segment");

// adding the single maneuver
circularizeLunarOrbitTargetedSegment.getSegments().add(lunarOrbitInsertionManeuver);

// add the circularizing TargetedSegmentList to the overall segment list
mainSegmentList.getSegments().add(circularizeLunarOrbitTargetedSegment);

// and the final propagate segment
mainSegmentList.getSegments().add(aroundTheMoonSegment);
Targeting

Now that we have our segments defined, we turn to figuring out how to solve for an orbit that actually gets us to the Moon. Our satellite will get nowhere close to the Moon if we just propagate using our rough initial guesses. We need to hone in on the exact solution to get to the Moon, after which we have the subsequent problem of getting into a stable orbit around the Moon. This implies that two target segments are needed.

When picking variables and constraints for a given problem, you must consider how much those constraints will change for a given change with the variable. For this example, ultimately we want to be in a lunar polar orbit with an altitude of 500 km. To get into that exact orbit, we need to perform very minor changes to our initial conditions. Otherwise, we are just going to jump into the middle of nowhere. But, to get anywhere near the Moon at all requires fairly large changes to the conditions.

This indicates that we should use multiple differential correctors with different levels of fidelity. For this problem, it makes sense to constrain our satellite's declination and right ascension to within a degree or two of the Moon relative to the Earth, then use a second solver that constrains the b-plane targeting to get the satellite into the polar orbit, and finally target the lunar inclination and periapsis altitude.

For this problem, the variables that make the most sense to vary are the initial right ascension of the ascending node, initial inclination. This will effectively swivel and tilt the outgoing orbit until it gets close to the moon. Note that this and all of the other inclination variables in this topic make use of a helper method that properly handles the inclination, as described in the Multivariable Function Solvers topic.

Java
Copy
// we want to modify the initial right ascension and Inclination so that our orbit swivels around to
// approach the moon

// right ascension variable
SegmentPropagatorVariable firstRightAscensionVariable = initialStateSegment.createVariable(
        Trig.degreesToRadians(10.0), // maximum step, radians
        Trig.degreesToRadians(0.2), // perturbation, radians
        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());
        }));
firstRightAscensionVariable.setName("Initial_State_InitialState_Keplerian_RightAscension");

// inclination variable
SegmentPropagatorVariable firstInclinationVariable = initialStateSegment.createVariable(
        Trig.degreesToRadians(10.0), // maximum step, radians
        Trig.degreesToRadians(0.2), // perturbation, radians
        SetVariableCallback.of((variable, configuration) -> {
            Motion1<Cartesian> currentMotion = configuration.getMotion(motionID);

            ModifiedKeplerianElements oldElements =
                    new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter);
            ModifiedKeplerianElements newElements =
                    computeElementsWithProperNewInclination(oldElements, variable);

            configuration.modifyMotion(motionID, newElements.toCartesian());
        }));
firstInclinationVariable.setName("Initial_State_InitialState_Keplerian_Inclination");
Java
Copy
public static ModifiedKeplerianElements computeElementsWithProperNewInclination(ModifiedKeplerianElements oldElements, double newInclinationDelta) {
    double rightAscensionCorrection = 0.0;
    double argumentOfPeriapsisCorrection = 0.0;
    double inclination = oldElements.getInclination() + newInclinationDelta;
    if (inclination > Math.PI) {
        // these 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 right ascension.
        inclination -= Math.PI;
        rightAscensionCorrection = -Math.PI;
        argumentOfPeriapsisCorrection = -Math.PI;
    } else if (inclination < 0) {
        // This is correct. If the inclination becomes less and less,
        // to avoid discontinuities in position and velocity,
        // the right ascension and argument of periapsis need to flip too.
        inclination += Math.PI;
        rightAscensionCorrection = Math.PI;
        argumentOfPeriapsisCorrection = Math.PI;
    }

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

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

    return new ModifiedKeplerianElements(
            oldElements.getRadiusOfPeriapsis(),
            oldElements.getInverseSemimajorAxis(),
            inclination,
            newArgumentOfPeriapsis,
            newRightAscension,
            oldElements.getTrueAnomaly(),
            oldElements.getGravitationalParameter());
}

Next, we make delta-declination and delta-right ascension constraints to line up our satellite with the Moon at the end of the transfer orbit:

Java
Copy
// for the first differential corrector just getting us in the vicinity of the moon,
// the delta-declination and delta-right ascension behave very well on the large scale
// we are initially solving for

// delta declination constraint
ScalarAtEndOfSegmentConstraint deltaDeclinationConstraint =
        new ScalarAtEndOfSegmentConstraint(transferOrbitWithinMoonsSoiSegment,
                                           0.0, // desired value, radians
                                           Trig.degreesToRadians(0.0001)); // tolerance, radians
ScalarDeltaSphericalElement deltaDeclinationScalar =
        new ScalarDeltaSphericalElement(new ParameterizedOnStatePoint(deltaDeclinationConstraint.getParameter(),
                                                                      moon.getInertialFrame(),
                                                                      motionID),
                                        earth.getCenterOfMassPoint(),
                                        moon.getCenterOfMassPoint(),
                                        earth.getInertialFrame(),
                                        SphericalElement.CONE);
deltaDeclinationConstraint.setScalar(deltaDeclinationScalar);
deltaDeclinationConstraint.setName("Delta-Declination");

// delta right ascension constraint
ScalarAtEndOfSegmentConstraint deltaRightAscensionConstraint =
        new ScalarAtEndOfSegmentConstraint(transferOrbitWithinMoonsSoiSegment,
                                           0.0, // desired value, radians
                                           Trig.degreesToRadians(0.0001)); // tolerance, radians
ScalarDeltaSphericalElement deltaRightAscensionScalar =
        new ScalarDeltaSphericalElement(new ParameterizedOnStatePoint(deltaRightAscensionConstraint.getParameter(),
                                                                      moon.getInertialFrame(),
                                                                      motionID),
                                        earth.getCenterOfMassPoint(),
                                        moon.getCenterOfMassPoint(),
                                        earth.getInertialFrame(),
                                        SphericalElement.CLOCK);
deltaRightAscensionConstraint.setScalar(deltaRightAscensionScalar);
deltaRightAscensionConstraint.setName("Delta-Right_Ascension");

Then, we create our first TargetedSegmentListDifferentialCorrector and add it to the first target segment:

Java
Copy
// make the differential corrector
TargetedSegmentListDifferentialCorrector firstApproachDifferentialCorrector =
        new TargetedSegmentListDifferentialCorrector();
firstApproachDifferentialCorrector.setName("Coarse_Solving_For_Lunar_Approach");

// add the constraints
firstApproachDifferentialCorrector.getConstraints().add(deltaDeclinationConstraint);
firstApproachDifferentialCorrector.getConstraints().add(deltaRightAscensionConstraint);

// add the variables
firstApproachDifferentialCorrector.getVariables().add(firstInclinationVariable);
firstApproachDifferentialCorrector.getVariables().add(firstRightAscensionVariable);

// add the differential corrector to the correct targeted segment list
targetedSegmentToGetToTheMoon.getOperators().add(firstApproachDifferentialCorrector);

We create our second TargetedSegmentListDifferentialCorrector, with conceptually the same variables, but because we want a finer precision than we had before, we must make new ones with a smaller perturbation and max step:

Java
Copy
// for the second differential corrector, we will use similar variables,
// but with a smaller maximum steps and perturbations to target with finer
// constraints, BDotT and BDotR

// second right ascension variable
SegmentPropagatorVariable secondRightAscensionVariable = initialStateSegment.createVariable(
        Trig.degreesToRadians(2.0), // smaller maximum step, radians
        Trig.degreesToRadians(0.1), // smaller perturbation, radians

        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());
        }));
secondRightAscensionVariable.setName("Initial_State_InitialState_Keplerian_RightAscension");

// second inclination variable
SegmentPropagatorVariable secondInclinationVariable = initialStateSegment.createVariable(
        Trig.degreesToRadians(2.0), // smaller maximum step, radians
        Trig.degreesToRadians(0.1), // smaller perturbation, radians

        SetVariableCallback.of((variable, configuration) -> {
            Motion1<Cartesian> currentMotion = configuration.getMotion(motionID);

            ModifiedKeplerianElements oldElements =
                    new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter);
            ModifiedKeplerianElements newElements =
                    computeElementsWithProperNewInclination(oldElements, variable);

            configuration.modifyMotion(motionID, newElements.toCartesian());
        }));
secondInclinationVariable.setName("Initial_State_InitialState_Keplerian_Inclination");

Here we will use values derived from a VectorBPlane as our constraints. The B-Plane applies to a spacecraft on a hyperbolic trajectory approaching a central body. The B-Plane is defined as the plane that is perpendicular to the incoming asymptote and passes through the center of the central body. The B-Plane vector is the vector on that plane from the central body to where the asymptote meets that plane. By targeting the X and Y elements of that vector (getBDotT and getBDotR, respectively), we can control our incoming orbit to a much better degree than we could with just the difference between the declination and right ascension of our satellite and the moon.

With the B-plane constraints, we assemble these into another TargetedSegmentListDifferentialCorrector and add that to the TargetedSegmentList:

Java
Copy
// the B-Plane constraints will allow us to do more precise targeting of
// our trajectory than the delta-declination/right ascension would.  But
// they still won't get us into our desired orbit

// a B dot T of 0 and a B dot R of 5000.0 will put us in a polar orbit.

// BDotT constraint
ScalarAtEndOfSegmentConstraint bDotTConstraint =
        new ScalarAtEndOfSegmentConstraint(transferOrbitWithinMoonsSoiSegment,
                                           0.0, // desired value, meters
                                           0.00001); // tolerance, meters
Scalar bDotTScalar =
        new VectorBPlane(new ParameterizedOnStatePoint(bDotTConstraint.getParameter(),
                                                       moon.getInertialFrame(),
                                                       motionID),
                         moon,
                         moonGravitationalParameter).getBDotT();
bDotTConstraint.setScalar(bDotTScalar);
bDotTConstraint.setName("BDotT");

ScalarAtEndOfSegmentConstraint bDotRConstraint =
        new ScalarAtEndOfSegmentConstraint(transferOrbitWithinMoonsSoiSegment,
                                           5000000.0, // desired value, meters
                                           0.00001);  // tolerance, meters
Scalar bDotRScalar = new VectorBPlane(new ParameterizedOnStatePoint(bDotRConstraint.getParameter(),
                                                                    moon.getInertialFrame(),
                                                                    motionID),
                                      moon,
                                      moonGravitationalParameter).getBDotR();
bDotRConstraint.setScalar(bDotRScalar);
bDotRConstraint.setName("BDotR");

We then create the second differential corrector, similar to the first:

Java
Copy
// make the second corrector
TargetedSegmentListDifferentialCorrector secondApproachDifferentialCorrector =
        new TargetedSegmentListDifferentialCorrector();
secondApproachDifferentialCorrector.setName("Fine_Solving_For_Lunar_Approach");

// add the constraints
secondApproachDifferentialCorrector.getConstraints().add(bDotRConstraint);
secondApproachDifferentialCorrector.getConstraints().add(bDotTConstraint);

// and the variables
secondApproachDifferentialCorrector.getVariables().add(secondRightAscensionVariable);
secondApproachDifferentialCorrector.getVariables().add(secondInclinationVariable);

// add to the targeted segment list
targetedSegmentToGetToTheMoon.getOperators().add(secondApproachDifferentialCorrector);

We need a third differential corrector to determine a unique orbit around the moon. The same variables will be used, but with again a smaller maximum step and perturbation. The constraints this time will be the lunar altitude and inclination. We start again with the variables:

Java
Copy
// The two previous differential correctors have put us in such an orbit that
// we can reliably target a specific inclination and lunar altitude at periapsis.

// third RightAscension variable
SegmentPropagatorVariable thirdRightAscensionVariable = initialStateSegment.createVariable(
        Trig.degreesToRadians(0.01),   // even smaller maximum step
        Trig.degreesToRadians(0.0005), // and a smaller perturbation
        SetVariableCallback.of((variable, configuration) -> {
            // get motion
            Motion1<Cartesian> currentMotion = configuration.getMotion(motionID);

            // get old elements
            ModifiedKeplerianElements oldElements =
                    new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter);
            // create new elements
            ModifiedKeplerianElements newElements =
                    new ModifiedKeplerianElements(
                            oldElements.getRadiusOfPeriapsis(),
                            oldElements.getInverseSemimajorAxis(),
                            oldElements.getInclination(),
                            oldElements.getArgumentOfPeriapsis(),
                            oldElements.getRightAscensionOfAscendingNode() + variable, // adding the variable
                            oldElements.getTrueAnomaly(),
                            earthGravitationalParameter);

            // apply the new motion
            configuration.modifyMotion(motionID, newElements.toCartesian());
        }));
thirdRightAscensionVariable.setName("Initial_State_InitialState_Keplerian_RightAscension");

SegmentPropagatorVariable thirdInclinationVariable = initialStateSegment.createVariable(
        Trig.degreesToRadians(0.01),   // even smaller maximum step
        Trig.degreesToRadians(0.0005), // and a smaller perturbation
        SetVariableCallback.of((variable, configuration) -> {
            Motion1<Cartesian> currentMotion = configuration.getMotion(motionID);

            ModifiedKeplerianElements oldElements =
                    new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter);
            ModifiedKeplerianElements newElements =
                    computeElementsWithProperNewInclination(oldElements, variable);

            configuration.modifyMotion(motionID, newElements.toCartesian());
        }));
thirdInclinationVariable.setName("Initial_State_InitialState_Keplerian_Inclination");

Now we add constraints for the lunar altitude and inclination at the end of the transfer orbit:

Java
Copy
// The propagate segment that ends at the moon has three stopping conditions:
// A duration of 2 days just so we don't propagate forever.
// A lunar altitude of 1 km
// And periapsis around the moon
// 
// We are targeting a specific lunar altitude with this first constraint, but
// since it is higher than the altitude stopping condition, we are in effect
// targeting the periapsis altitude.  In this case it shouldn't matter which of
// those you constraint the trajectory too, but be aware that you have many choice
// like that.

// lunar altitude constraint
ScalarAtEndOfSegmentConstraint lunarAltitudeConstraint =
        new ScalarAtEndOfSegmentConstraint(transferOrbitWithinMoonsSoiSegment,
                                           500000.0,  // desired value, meters
                                           0.000001); // tolerance, meters
ScalarCartographicElement lunarAltitudeScalar =
        new ScalarCartographicElement(moon,
                                      new ParameterizedOnStatePoint(lunarAltitudeConstraint.getParameter(),
                                                                    moon.getInertialFrame(),
                                                                    motionID),
                                      CartographicElement.HEIGHT);
lunarAltitudeConstraint.setScalar(lunarAltitudeScalar);
lunarAltitudeConstraint.setName("Lunar_Altitude_Constraint");

// lunar inclination constraint, we want a polar orbit
ScalarAtEndOfSegmentConstraint lunarInclinationConstraint =
        new ScalarAtEndOfSegmentConstraint(transferOrbitWithinMoonsSoiSegment,
                                           Math.PI / 2.0, // desired value, radians
                                           0.00000001); // tolerance, radians
ScalarModifiedKeplerianElement lunarInclinationScalar =
        new ScalarModifiedKeplerianElement(moonGravitationalParameter,
                                           new ParameterizedOnStatePoint(lunarInclinationConstraint.getParameter(),
                                                                         moon.getInertialFrame(),
                                                                         motionID),
                                           KeplerianElement.INCLINATION,
                                           moon.getInertialFrame());
lunarInclinationConstraint.setScalar(lunarInclinationScalar);
lunarInclinationConstraint.setName("Lunar_Inclination_Constraint");

And similar to the previous differential correctors, we make the third:

Java
Copy
// make the differential corrector
TargetedSegmentListDifferentialCorrector thirdApproachDifferentialCorrector =
        new TargetedSegmentListDifferentialCorrector();
thirdApproachDifferentialCorrector.setName("Precise_Differential_Corrector");

// add the variables
thirdApproachDifferentialCorrector.getVariables().add(thirdInclinationVariable);
thirdApproachDifferentialCorrector.getVariables().add(thirdRightAscensionVariable);

// add the constraints
thirdApproachDifferentialCorrector.getConstraints().add(lunarAltitudeConstraint);
thirdApproachDifferentialCorrector.getConstraints().add(lunarInclinationConstraint);

// add the differential corrector to the targeted segment list
targetedSegmentToGetToTheMoon.getOperators().add(thirdApproachDifferentialCorrector);

After that, the final differential corrector is rather straight forward. We vary the x component of the delta-v in the lunar VNC axes such that our orbital eccentricity is 0.

Java
Copy
SegmentPropagatorVariable lunarCaptureBurnVariable = lunarOrbitInsertionManeuver.createVariable(
        1000.0, // maximum step, meters/second
        0.1,    // perturbation, meters/second
        SetVariableCallback.of((variable, configuration) -> {
            configuration.getManeuver().setX(configuration.getManeuver().getX() + variable);
        }));
lunarCaptureBurnVariable.setName("ImpulsiveManeuver_ThrustVector_X");

ScalarAtEndOfSegmentConstraint eccentricityAroundMoonConstraint =
        new ScalarAtEndOfSegmentConstraint(lunarOrbitInsertionManeuver,
                                           0.0, // desired value, unitless
                                           0.00003); // tolerance, unitless
ScalarModifiedKeplerianElement eccentricityAroundMoonScalar =
        new ScalarModifiedKeplerianElement(moonGravitationalParameter,
                                           new ParameterizedOnStatePoint(eccentricityAroundMoonConstraint.getParameter(),
                                                                         moon.getInertialFrame(),
                                                                         motionID),
                                           KeplerianElement.ECCENTRICITY,
                                           moon.getInertialFrame());
eccentricityAroundMoonConstraint.setScalar(eccentricityAroundMoonScalar);
eccentricityAroundMoonConstraint.setName("Eccentricity");

TargetedSegmentListDifferentialCorrector lunarCaptureDifferentialCorrector =
        new TargetedSegmentListDifferentialCorrector();
lunarCaptureDifferentialCorrector.setName("Circularize_Lunar_Orbit");

lunarCaptureDifferentialCorrector.getVariables().add(lunarCaptureBurnVariable);
lunarCaptureDifferentialCorrector.getConstraints().add(eccentricityAroundMoonConstraint);

circularizeLunarOrbitTargetedSegment.getOperators().add(lunarCaptureDifferentialCorrector);
Propagating

Before we propagate, we set some general settings on all of our differential correctors.

Java
Copy
firstApproachDifferentialCorrector.setMaximumIterations(20);
secondApproachDifferentialCorrector.setMaximumIterations(20);
thirdApproachDifferentialCorrector.setMaximumIterations(50); // since we are taking so small steps
lunarCaptureDifferentialCorrector.setMaximumIterations(20);

firstApproachDifferentialCorrector.getSolver().setMultithreaded(isMultithreaded);
secondApproachDifferentialCorrector.getSolver().setMultithreaded(isMultithreaded);
thirdApproachDifferentialCorrector.getSolver().setMultithreaded(isMultithreaded);
lunarCaptureDifferentialCorrector.getSolver().setMultithreaded(isMultithreaded);

Finally, we get the SegmentPropagator from the overall SegmentList, and propagate. When it finishes, you can use the SegmentResults to perform further analysis.

Java
Copy
SegmentPropagator propagator = mainSegmentList.getSegmentPropagator(new EvaluatorGroup());
SegmentListResults listResults = (SegmentListResults) propagator.propagate();
DateMotionCollection1<Cartesian> overallEphemeris =
        listResults.getDateMotionCollectionOfOverallTrajectory(motionID,
                                                               earth.getInertialFrame());