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 |
|---|
The functionality described in this topic requires a license for the Segment Propagation Library. |
Designing SegmentsTo 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:
// pick names const string motionID = "Satellite"; const string fuelMassName = "Fuel_Mass"; const string dryMassName = "Dry_Mass"; // store commonly used central bodies and constants CentralBody moon = CentralBodiesFacet.GetFromContext().Moon; double moonGravitationalParameter = 4902801076000.0; CentralBody earth = CentralBodiesFacet.GetFromContext().Earth; double earthGravitationalParameter = WorldGeodeticSystem1984.GravitationalParameter; // initial date JulianDate epoch = new JulianDate(new GregorianDate(2009, 7, 1, 15, 0, 0.0), TimeStandard.CoordinatedUniversalTime); // initial state Motion<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:
// fuel mass PropagationScalar fuelMass = new PropagationScalar(2000.0); fuelMass.Identification = fuelMassName; fuelMass.ScalarDerivative = new ScalarFixed(0.0); // dry mass PropagationScalar dryMass = new PropagationScalar(100.0); dryMass.Identification = dryMassName; dryMass.ScalarDerivative = new ScalarFixed(0.0); // Create the integrator RungeKutta4Integrator integrator = new RungeKutta4Integrator(); integrator.InitialStepSize = 120.0;
Next, we create the actual NumericalPropagatorDefinition that we will use for our Earth-centered segments:
// the propagation point PropagationNewtonianPoint propagationPointAroundEarth = new PropagationNewtonianPoint(motionID, earth.InertialFrame, initialState.Value, initialState.FirstDerivative); // total mass propagationPointAroundEarth.Mass = fuelMass.IntegrationValue + dryMass.IntegrationValue; // simple gravity for Earth TwoBodyGravity earthGravityAsPrimary = new TwoBodyGravity(propagationPointAroundEarth.IntegrationPoint, earth, earthGravitationalParameter); propagationPointAroundEarth.AppliedForces.Add(earthGravityAsPrimary); // and the Moons gravity ThirdBodyGravity moonGravityAsThirdBody = new ThirdBodyGravity(); moonGravityAsThirdBody.TargetPoint = propagationPointAroundEarth.IntegrationPoint; moonGravityAsThirdBody.AddThirdBody(moon.Name, moon.CenterOfMassPoint, moonGravitationalParameter); propagationPointAroundEarth.AppliedForces.Add(moonGravityAsThirdBody); // make the actual propagator NumericalPropagatorDefinition earthCenteredNumericalPropagator = new NumericalPropagatorDefinition(); // add the elements earthCenteredNumericalPropagator.IntegrationElements.Add(propagationPointAroundEarth); earthCenteredNumericalPropagator.IntegrationElements.Add(fuelMass); earthCenteredNumericalPropagator.IntegrationElements.Add(dryMass); // set the integrator earthCenteredNumericalPropagator.Integrator = integrator; // set epoch earthCenteredNumericalPropagator.Epoch = epoch;
Initial SegmentsWe 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:
// Simply the initial state. We will need it since the differential corrector will operate // on the initial values. NumericalInitialStateSegment initialStateSegment = new NumericalInitialStateSegment(); initialStateSegment.Name = "Initial_State_Segment"; initialStateSegment.PropagatorDefinition = 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.
// configure the segment NumericalPropagatorSegment parkingOrbitSegment = new NumericalPropagatorSegment(); parkingOrbitSegment.Name = "Parking_Orbit"; parkingOrbitSegment.PropagatorDefinition = earthCenteredNumericalPropagator; // propagate for almost an hour DurationStoppingCondition timeOnParkingOrbit = new DurationStoppingCondition(Duration.FromSeconds(3200.0)); timeOnParkingOrbit.Name = "Duration_On_Parking_Orbit"; parkingOrbitSegment.StoppingConditions.Add(timeOnParkingOrbit);
Lunar TransferNext, we create our first maneuver getting us onto a lunar transfer orbit. We want the delta-v in the velocity-orbit-normal axes.
ImpulsiveManeuverSegment transLunarInjectionSegment = new ImpulsiveManeuverSegment(); transLunarInjectionSegment.Name = "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, exhaustVelocity, InvalidFuelStateBehavior.ThrowException); transLunarBurnDetails.Orientation = new AxesVelocityOrbitNormal(transLunarBurnDetails.PropagationPoint, earth); transLunarInjectionSegment.Maneuver = 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.
// propagate until we are inside the moons sphere of influence NumericalPropagatorSegment propagateOutOfEarthsSoiSegment = new NumericalPropagatorSegment(); propagateOutOfEarthsSoiSegment.Name = "Propagate_To_Moons_SOI"; propagateOutOfEarthsSoiSegment.PropagatorDefinition = earthCenteredNumericalPropagator; // stop propagating this segment at the SOI boundary ScalarStoppingCondition awayFromEarthStoppingCondition = new ScalarStoppingCondition(new ScalarPointElement(propagationPointAroundEarth.IntegrationPoint, CartesianElement.Magnitude, earth.InertialFrame), 320000000, // approximate SOI for moon, meters 0.0001, // value tolerance, meters StopType.AnyThreshold); awayFromEarthStoppingCondition.Name = "Moon_SOI"; propagateOutOfEarthsSoiSegment.StoppingConditions.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.
PropagationNewtonianPoint propagationPointAroundMoon = new PropagationNewtonianPoint(); propagationPointAroundMoon.Identification = 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.InitialPosition = initialState.Value; propagationPointAroundMoon.InitialVelocity = initialState.FirstDerivative; // 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.IntegrationFrame = moon.InertialFrame; // Configure fuel. Again, the values will get initialized from the previous // segment's final state. propagationPointAroundMoon.Mass = dryMass.IntegrationValue + fuelMass.IntegrationValue; // Moon's primary gravity TwoBodyGravity moonPrimaryGravity = new TwoBodyGravity(propagationPointAroundMoon.IntegrationPoint, moon, moonGravitationalParameter); propagationPointAroundMoon.AppliedForces.Add(moonPrimaryGravity); // Earth gravity as third body gravity ThirdBodyGravity earthGravityAsThirdBody = new ThirdBodyGravity(propagationPointAroundMoon.IntegrationPoint); earthGravityAsThirdBody.AddThirdBody(earth.Name, earth.CenterOfMassPoint, earthGravitationalParameter); earthGravityAsThirdBody.CentralBody = moon; propagationPointAroundMoon.AppliedForces.Add(earthGravityAsThirdBody); // make the actual propagator NumericalPropagatorDefinition nearMoonNumericalPropagator = new NumericalPropagatorDefinition(); // add the elements nearMoonNumericalPropagator.IntegrationElements.Add(propagationPointAroundMoon); nearMoonNumericalPropagator.IntegrationElements.Add(dryMass); nearMoonNumericalPropagator.IntegrationElements.Add(fuelMass); // set the integrator nearMoonNumericalPropagator.Integrator = integrator; // set an epoch nearMoonNumericalPropagator.Epoch = epoch;
Lunar OrbitFinally, 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.
NumericalPropagatorSegment transferOrbitWithinMoonsSoiSegment = new NumericalPropagatorSegment(); transferOrbitWithinMoonsSoiSegment.Name = "Transfer_Orbit_Approaching_Moon"; // set the NumericalPropagator on the propagate segment transferOrbitWithinMoonsSoiSegment.PropagatorDefinition = 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.IntegrationPoint, CartographicElement.Height), 1000.0, // threshold, meters 0.1, // value tolerance, meters StopType.AnyThreshold); altitudeOfMoonStoppingCondition.Name = "Lunar_Altitude"; transferOrbitWithinMoonsSoiSegment.StoppingConditions.Add(altitudeOfMoonStoppingCondition); // periapsis stopping condition ScalarStoppingCondition periluneStoppingCondition = new ScalarStoppingCondition(new ScalarModifiedKeplerianElement(moonGravitationalParameter, propagationPointAroundMoon.IntegrationPoint, KeplerianElement.TrueAnomaly, moon.InertialFrame), 0.0, // threshold, radians 0.000001, // tolerance, radians StopType.AnyThreshold); // leave the threshold 0.0 radians periluneStoppingCondition.AngularSetting = CircularRange.NegativePiToPi; periluneStoppingCondition.Name = "Lunar_Periapsis"; transferOrbitWithinMoonsSoiSegment.StoppingConditions.Add(periluneStoppingCondition); // and stop on Duration just in case we completely miss DurationStoppingCondition durationStoppingConditionWhenApproachingMoon = new DurationStoppingCondition(Duration.FromDays(2.0)); transferOrbitWithinMoonsSoiSegment.StoppingConditions.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 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:
ImpulsiveManeuverSegment lunarOrbitInsertionManeuver = new ImpulsiveManeuverSegment(); lunarOrbitInsertionManeuver.Name = "Lunar_Orbit_Insertion"; // the maneuver details ImpulsiveManeuverInformation lunarOrbitInsertionManeuverDetails = new ImpulsiveManeuverInformation(motionID, new Cartesian(-1000.0, 0.0, 0.0), fuelMassName, dryMassName, exhaustVelocity, InvalidFuelStateBehavior.ThrowException); lunarOrbitInsertionManeuverDetails.Orientation = new AxesVelocityOrbitNormal(lunarOrbitInsertionManeuverDetails.PropagationPoint, moon); lunarOrbitInsertionManeuver.Maneuver = lunarOrbitInsertionManeuverDetails;
Finally, we propagate around the moon for three days:
// final segment for this example, we made it to the Moon! NumericalPropagatorSegment aroundTheMoonSegment = new NumericalPropagatorSegment(); aroundTheMoonSegment.Name = "Lunar_Orbit"; // same propagator aroundTheMoonSegment.PropagatorDefinition = nearMoonNumericalPropagator; // propagate and do science aroundTheMoonSegment.StoppingConditions.Add(new DurationStoppingCondition(Duration.FromDays(3.0)));
Segment ListsBefore 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:
// make the main list SegmentList mainSegmentList = new SegmentList(); mainSegmentList.Name = "Main_Segment_List"; // we want to target items in the initial state, so we need a TargetedSegmentList TargetedSegmentList targetedSegmentToGetToTheMoon = new TargetedSegmentList(); targetedSegmentToGetToTheMoon.Name = "Targeted_Segment_To_Get_To_Moon"; targetedSegmentToGetToTheMoon.Segments.Add(initialStateSegment); targetedSegmentToGetToTheMoon.Segments.Add(parkingOrbitSegment); targetedSegmentToGetToTheMoon.Segments.Add(transLunarInjectionSegment); targetedSegmentToGetToTheMoon.Segments.Add(propagateOutOfEarthsSoiSegment); targetedSegmentToGetToTheMoon.Segments.Add(transferOrbitWithinMoonsSoiSegment); // add that TargetedSegmentList to the overall segment list mainSegmentList.Segments.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.Name = "Lunar_Orbit_Insertion_Targeted_Segment"; // adding the single maneuver circularizeLunarOrbitTargetedSegment.Segments.Add(lunarOrbitInsertionManeuver); // add the circularizing TargetedSegmentList to the overall segment list mainSegmentList.Segments.Add(circularizeLunarOrbitTargetedSegment); // and the final propagate segment mainSegmentList.Segments.Add(aroundTheMoonSegment);
TargetingNow 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.
// 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 (variable, configuration) => { Motion<Cartesian> currentMotion = configuration.GetMotion<Cartesian>(motionID); ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter); ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(oldElements.RadiusOfPeriapsis, oldElements.InverseSemimajorAxis, oldElements.Inclination, oldElements.ArgumentOfPeriapsis, oldElements.RightAscensionOfAscendingNode + variable, oldElements.TrueAnomaly, earthGravitationalParameter); configuration.ModifyMotion(motionID, newElements.ToCartesian()); }); firstRightAscensionVariable.Name = "Initial_State_InitialState_Keplerian_RightAscension"; // inclination variable SegmentPropagatorVariable firstInclinationVariable = initialStateSegment.CreateVariable( Trig.DegreesToRadians(10.0), // maximum step, radians Trig.DegreesToRadians(0.2), // perturbation, radians (variable, configuration) => { Motion<Cartesian> currentMotion = configuration.GetMotion<Cartesian>(motionID); ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter); ModifiedKeplerianElements newElements = ComputeElementsWithProperNewInclination(oldElements, variable); configuration.ModifyMotion(motionID, newElements.ToCartesian()); }); firstInclinationVariable.Name = "Initial_State_InitialState_Keplerian_Inclination";
public static ModifiedKeplerianElements ComputeElementsWithProperNewInclination( ModifiedKeplerianElements oldElements, double newInclinationDelta) { double rightAscensionCorrection = 0.0; double argumentOfPeriapsisCorrection = 0.0; double inclination = oldElements.Inclination + 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.ArgumentOfPeriapsis + argumentOfPeriapsisCorrection); double newRightAscension = Trig.ZeroToTwoPi(oldElements.RightAscensionOfAscendingNode + rightAscensionCorrection); return new ModifiedKeplerianElements(oldElements.RadiusOfPeriapsis, oldElements.InverseSemimajorAxis, inclination, newArgumentOfPeriapsis, newRightAscension, oldElements.TrueAnomaly, oldElements.GravitationalParameter); }
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:
// 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.Parameter, moon.InertialFrame, motionID), earth.CenterOfMassPoint, moon.CenterOfMassPoint, earth.InertialFrame, SphericalElement.Cone); deltaDeclinationConstraint.Scalar = deltaDeclinationScalar; deltaDeclinationConstraint.Name = "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.Parameter, moon.InertialFrame, motionID), earth.CenterOfMassPoint, moon.CenterOfMassPoint, earth.InertialFrame, SphericalElement.Clock); deltaRightAscensionConstraint.Scalar = deltaRightAscensionScalar; deltaRightAscensionConstraint.Name = "Delta-Right_Ascension";
Then, we create our first TargetedSegmentListDifferentialCorrector and add it to the first target segment:
// make the differential corrector TargetedSegmentListDifferentialCorrector firstApproachDifferentialCorrector = new TargetedSegmentListDifferentialCorrector(); firstApproachDifferentialCorrector.Name = "Coarse_Solving_For_Lunar_Approach"; // add the constraints firstApproachDifferentialCorrector.Constraints.Add(deltaDeclinationConstraint); firstApproachDifferentialCorrector.Constraints.Add(deltaRightAscensionConstraint); // add the variables firstApproachDifferentialCorrector.Variables.Add(firstInclinationVariable); firstApproachDifferentialCorrector.Variables.Add(firstRightAscensionVariable); // add the differential corrector to the correct targeted segment list targetedSegmentToGetToTheMoon.Operators.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:
// 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 (variable, configuration) => { Motion<Cartesian> currentMotion = configuration.GetMotion<Cartesian>(motionID); ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter); ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(oldElements.RadiusOfPeriapsis, oldElements.InverseSemimajorAxis, oldElements.Inclination, oldElements.ArgumentOfPeriapsis, oldElements.RightAscensionOfAscendingNode + variable, oldElements.TrueAnomaly, earthGravitationalParameter); configuration.ModifyMotion(motionID, newElements.ToCartesian()); }); secondRightAscensionVariable.Name = "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 (variable, configuration) => { Motion<Cartesian> currentMotion = configuration.GetMotion<Cartesian>(motionID); ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter); ModifiedKeplerianElements newElements = ComputeElementsWithProperNewInclination(oldElements, variable); configuration.ModifyMotion(motionID, newElements.ToCartesian()); }); secondInclinationVariable.Name = "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:
// 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.Parameter, moon.InertialFrame, motionID), moon, moonGravitationalParameter).GetBDotT(); bDotTConstraint.Scalar = bDotTScalar; bDotTConstraint.Name = "BDotT"; ScalarAtEndOfSegmentConstraint bDotRConstraint = new ScalarAtEndOfSegmentConstraint(transferOrbitWithinMoonsSoiSegment, 5000000, // desired value, meters 0.00001); // tolerance, meters Scalar bDotRScalar = new VectorBPlane(new ParameterizedOnStatePoint(bDotRConstraint.Parameter, moon.InertialFrame, motionID), moon, moonGravitationalParameter).GetBDotR(); bDotRConstraint.Scalar = bDotRScalar; bDotRConstraint.Name = "BDotR";
We then create the second differential corrector, similar to the first:
// make the second corrector TargetedSegmentListDifferentialCorrector secondApproachDifferentialCorrector = new TargetedSegmentListDifferentialCorrector(); secondApproachDifferentialCorrector.Name = "Fine_Solving_For_Lunar_Approach"; // add the constraints secondApproachDifferentialCorrector.Constraints.Add(bDotRConstraint); secondApproachDifferentialCorrector.Constraints.Add(bDotTConstraint); // and the variables secondApproachDifferentialCorrector.Variables.Add(secondRightAscensionVariable); secondApproachDifferentialCorrector.Variables.Add(secondInclinationVariable); // add to the targeted segment list targetedSegmentToGetToTheMoon.Operators.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:
// 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 (variable, configuration) => { // get motion Motion<Cartesian> currentMotion = configuration.GetMotion<Cartesian>(motionID); // get old elements ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter); // create new elements ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(oldElements.RadiusOfPeriapsis, oldElements.InverseSemimajorAxis, oldElements.Inclination, oldElements.ArgumentOfPeriapsis, oldElements.RightAscensionOfAscendingNode + variable, // adding the variable oldElements.TrueAnomaly, earthGravitationalParameter); // apply the new motion configuration.ModifyMotion(motionID, newElements.ToCartesian()); }); thirdRightAscensionVariable.Name = "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 (variable, configuration) => { Motion<Cartesian> currentMotion = configuration.GetMotion<Cartesian>(motionID); ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, earthGravitationalParameter); ModifiedKeplerianElements newElements = ComputeElementsWithProperNewInclination(oldElements, variable); configuration.ModifyMotion(motionID, newElements.ToCartesian()); }); thirdInclinationVariable.Name = "Initial_State_InitialState_Keplerian_Inclination";
Now we add constraints for the lunar altitude and inclination at the end of the transfer orbit:
// 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.Parameter, moon.InertialFrame, motionID), CartographicElement.Height); lunarAltitudeConstraint.Scalar = lunarAltitudeScalar; lunarAltitudeConstraint.Name = "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.Parameter, moon.InertialFrame, motionID), KeplerianElement.Inclination, moon.InertialFrame); lunarInclinationConstraint.Scalar = lunarInclinationScalar; lunarInclinationConstraint.Name = "Lunar_Inclination_Constraint";
And similar to the previous differential correctors, we make the third:
// make the differential corrector TargetedSegmentListDifferentialCorrector thirdApproachDifferentialCorrector = new TargetedSegmentListDifferentialCorrector(); thirdApproachDifferentialCorrector.Name = "Precise_Differential_Corrector"; // add the variables thirdApproachDifferentialCorrector.Variables.Add(thirdInclinationVariable); thirdApproachDifferentialCorrector.Variables.Add(thirdRightAscensionVariable); // add the constraints thirdApproachDifferentialCorrector.Constraints.Add(lunarAltitudeConstraint); thirdApproachDifferentialCorrector.Constraints.Add(lunarInclinationConstraint); // add the differential corrector to the targeted segment list targetedSegmentToGetToTheMoon.Operators.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.
SegmentPropagatorVariable lunarCaptureBurnVariable = lunarOrbitInsertionManeuver.CreateVariable(
1000.0, // maximum step, meters/second
0.1, // perturbation, meters/second
(variable, configuration) =>
{
configuration.Maneuver.X += variable;
});
lunarCaptureBurnVariable.Name = "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.Parameter,
moon.InertialFrame,
motionID),
KeplerianElement.Eccentricity,
moon.InertialFrame);
eccentricityAroundMoonConstraint.Scalar = eccentricityAroundMoonScalar;
eccentricityAroundMoonConstraint.Name = "Eccentricity";
TargetedSegmentListDifferentialCorrector lunarCaptureDifferentialCorrector =
new TargetedSegmentListDifferentialCorrector();
lunarCaptureDifferentialCorrector.Name = "Circularize_Lunar_Orbit";
lunarCaptureDifferentialCorrector.Variables.Add(lunarCaptureBurnVariable);
lunarCaptureDifferentialCorrector.Constraints.Add(eccentricityAroundMoonConstraint);
circularizeLunarOrbitTargetedSegment.Operators.Add(lunarCaptureDifferentialCorrector);PropagatingBefore we propagate, we set some general settings on all of our differential correctors.
firstApproachDifferentialCorrector.MaximumIterations = 20; secondApproachDifferentialCorrector.MaximumIterations = 20; thirdApproachDifferentialCorrector.MaximumIterations = 50; // since we are taking so small steps lunarCaptureDifferentialCorrector.MaximumIterations = 20; firstApproachDifferentialCorrector.Solver.Multithreaded = isMultithreaded; secondApproachDifferentialCorrector.Solver.Multithreaded = isMultithreaded; thirdApproachDifferentialCorrector.Solver.Multithreaded = isMultithreaded; lunarCaptureDifferentialCorrector.Solver.Multithreaded = isMultithreaded;
Finally, we get the SegmentPropagator from the overall SegmentList, and propagate. When it finishes, you can use the SegmentResults to perform further analysis.
SegmentPropagator propagator = mainSegmentList.GetSegmentPropagator(new EvaluatorGroup());
SegmentListResults listResults = (SegmentListResults)propagator.Propagate();
DateMotionCollection<Cartesian> overallEphemeris =
listResults.GetDateMotionCollectionOfOverallTrajectory(motionID,
earth.InertialFrame);