# GMAT: Finite burn targeting, solving for burn duration and quaternion

I'm trying to using GMAT to design a finite burn with a fixed attitude throughout the burn. I do not know the attitude needed at the start, nor do I know the burn duration needed.

I set up six variables: BurnDuration, q1, q2, q3, and q4, respectively representing the burn duration in seconds and the components of the quaternion (I don't know which component corresponds to the scalar value). In the differential correcter, I ask GMAT to vary these variables, propagate to periapse, turn on the thruster for the burn duration, and then propagate until the next apoapse, where the objectives should be met. In the script below, I'm simply trying to increase the SMA by 100km.

In my setup, GMAT does not seem to modify the quaternions representing the attitude of the spacecraft (I believe that I fixed the thruster to the spacecraft), and only tries to modify the burn duration. After a number of iterations, the differential corrector tried to modify the burn duration by the max step in each direction (i.e. +30 seconds at iteration i and then immediately by -30 seconds, but never tries a smaller step). The differential correction does not converge at all.

Question: Why does GMAT not try to vary the quaternions? If they do not change the final value, then GMAT typically says that a given variable is ineffective.

I suspect that I have incorrectly set up something, but I cannot seem to figure out what it is. Any hint is appreciated, thanks.

The script:

%General Mission Analysis Tool(GMAT) Script
%Created: 2021-12-08 17:44:08

%----------------------------------------
%---------- Spacecraft
%----------------------------------------

Create Spacecraft DefaultSC;
GMAT DefaultSC.DateFormat = UTCGregorian;
GMAT DefaultSC.Epoch = '25 Nov 2023 14:11:26.789';
GMAT DefaultSC.CoordinateSystem = MoonJ2K;
GMAT DefaultSC.DisplayStateType = Cartesian;
GMAT DefaultSC.X = -406.8997620000155;
GMAT DefaultSC.Y = -111.8177819999983;
GMAT DefaultSC.Z = -1697.662681000002;
GMAT DefaultSC.VX = 0.620734;
GMAT DefaultSC.VY = 1.591183;
GMAT DefaultSC.VZ = -0.252232;
GMAT DefaultSC.DryMass = 850;
GMAT DefaultSC.Cd = 2.2;
GMAT DefaultSC.Cr = 1.8;
GMAT DefaultSC.DragArea = 15;
GMAT DefaultSC.SRPArea = 1;
GMAT DefaultSC.Tanks = {ChemicalTank1};
GMAT DefaultSC.Thrusters = {ChemicalThruster1};
GMAT DefaultSC.NAIFId = -10000001;
GMAT DefaultSC.NAIFIdReferenceFrame = -9000001;
GMAT DefaultSC.OrbitColor = Red;
GMAT DefaultSC.TargetColor = Teal;
GMAT DefaultSC.OrbitErrorCovariance = [ 1e+70 0 0 0 0 0 ; 0 1e+70 0 0 0 0 ; 0 0 1e+70 0 0 0 ; 0 0 0 1e+70 0 0 ; 0 0 0 0 1e+70 0 ; 0 0 0 0 0 1e+70 ];
GMAT DefaultSC.CdSigma = 1e+70;
GMAT DefaultSC.CrSigma = 1e+70;
GMAT DefaultSC.Id = 'SatId';
GMAT DefaultSC.Attitude = CoordinateSystemFixed;
GMAT DefaultSC.ModelFile = 'aura.3ds';
GMAT DefaultSC.ModelOffsetX = 0;
GMAT DefaultSC.ModelOffsetY = 0;
GMAT DefaultSC.ModelOffsetZ = 0;
GMAT DefaultSC.ModelRotationX = 0;
GMAT DefaultSC.ModelRotationY = 0;
GMAT DefaultSC.ModelRotationZ = 0;
GMAT DefaultSC.ModelScale = 1;
GMAT DefaultSC.AttitudeDisplayStateType = 'Quaternion';
GMAT DefaultSC.AttitudeRateDisplayStateType = 'AngularVelocity';
GMAT DefaultSC.AttitudeCoordinateSystem = EarthMJ2000Eq;
GMAT DefaultSC.EulerAngleSequence = '321';

%----------------------------------------
%---------- Hardware Components
%----------------------------------------

Create ChemicalThruster ChemicalThruster1;
GMAT ChemicalThruster1.CoordinateSystem = Local;
GMAT ChemicalThruster1.Origin = Luna;
GMAT ChemicalThruster1.Axes = SpacecraftBody;
GMAT ChemicalThruster1.ThrustDirection1 = 1;
GMAT ChemicalThruster1.ThrustDirection2 = 0;
GMAT ChemicalThruster1.ThrustDirection3 = 0;
GMAT ChemicalThruster1.DutyCycle = 1;
GMAT ChemicalThruster1.ThrustScaleFactor = 1;
GMAT ChemicalThruster1.DecrementMass = false;
GMAT ChemicalThruster1.Tank = {ChemicalTank1};
GMAT ChemicalThruster1.MixRatio = [ 1 ];
GMAT ChemicalThruster1.GravitationalAccel = 9.81;
GMAT ChemicalThruster1.C1 = 10;
GMAT ChemicalThruster1.C2 = 0;
GMAT ChemicalThruster1.C3 = 0;
GMAT ChemicalThruster1.C4 = 0;
GMAT ChemicalThruster1.C5 = 0;
GMAT ChemicalThruster1.C6 = 0;
GMAT ChemicalThruster1.C7 = 0;
GMAT ChemicalThruster1.C8 = 0;
GMAT ChemicalThruster1.C9 = 0;
GMAT ChemicalThruster1.C10 = 0;
GMAT ChemicalThruster1.C11 = 0;
GMAT ChemicalThruster1.C12 = 0;
GMAT ChemicalThruster1.C13 = 0;
GMAT ChemicalThruster1.C14 = 0;
GMAT ChemicalThruster1.C15 = 0;
GMAT ChemicalThruster1.C16 = 0;
GMAT ChemicalThruster1.K1 = 300;
GMAT ChemicalThruster1.K2 = 0;
GMAT ChemicalThruster1.K3 = 0;
GMAT ChemicalThruster1.K4 = 0;
GMAT ChemicalThruster1.K5 = 0;
GMAT ChemicalThruster1.K6 = 0;
GMAT ChemicalThruster1.K7 = 0;
GMAT ChemicalThruster1.K8 = 0;
GMAT ChemicalThruster1.K9 = 0;
GMAT ChemicalThruster1.K10 = 0;
GMAT ChemicalThruster1.K11 = 0;
GMAT ChemicalThruster1.K12 = 0;
GMAT ChemicalThruster1.K13 = 0;
GMAT ChemicalThruster1.K14 = 0;
GMAT ChemicalThruster1.K15 = 0;
GMAT ChemicalThruster1.K16 = 0;

Create ChemicalTank ChemicalTank1;
GMAT ChemicalTank1.AllowNegativeFuelMass = false;
GMAT ChemicalTank1.FuelMass = 756;
GMAT ChemicalTank1.Pressure = 1500;
GMAT ChemicalTank1.Temperature = 20;
GMAT ChemicalTank1.RefTemperature = 20;
GMAT ChemicalTank1.Volume = 0.75;
GMAT ChemicalTank1.FuelDensity = 1260;
GMAT ChemicalTank1.PressureModel = PressureRegulated;

%----------------------------------------
%---------- ForceModels
%----------------------------------------

Create ForceModel DefaultProp_ForceModel;
GMAT DefaultProp_ForceModel.CentralBody = Luna;
GMAT DefaultProp_ForceModel.PrimaryBodies = {Luna};
GMAT DefaultProp_ForceModel.Drag = None;
GMAT DefaultProp_ForceModel.SRP = Off;
GMAT DefaultProp_ForceModel.RelativisticCorrection = Off;
GMAT DefaultProp_ForceModel.GravityField.Luna.Degree = 0;
GMAT DefaultProp_ForceModel.GravityField.Luna.Order = 0;
GMAT DefaultProp_ForceModel.GravityField.Luna.StmLimit = 100;
GMAT DefaultProp_ForceModel.GravityField.Luna.PotentialFile = 'LP165P.cof';
GMAT DefaultProp_ForceModel.GravityField.Luna.TideModel = 'None';

%----------------------------------------
%---------- Propagators
%----------------------------------------

Create Propagator DefaultProp;
GMAT DefaultProp.FM = DefaultProp_ForceModel;
GMAT DefaultProp.Type = RungeKutta89;
GMAT DefaultProp.InitialStepSize = 60;
GMAT DefaultProp.Accuracy = 9.999999999999999e-12;
GMAT DefaultProp.MinStep = 0.001;
GMAT DefaultProp.MaxStep = 2700;
GMAT DefaultProp.MaxStepAttempts = 50;
GMAT DefaultProp.StopIfAccuracyIsViolated = true;

%----------------------------------------
%---------- Burns
%----------------------------------------

Create FiniteBurn DefaultFB;
GMAT DefaultFB.Thrusters = {ChemicalThruster1};
GMAT DefaultFB.ThrottleLogicAlgorithm = 'MaxNumberOfThrusters';

Create ImpulsiveBurn DefaultIB;
GMAT DefaultIB.CoordinateSystem = Local;
GMAT DefaultIB.Origin = Earth;
GMAT DefaultIB.Axes = VNB;
GMAT DefaultIB.Element1 = 0;
GMAT DefaultIB.Element2 = 0;
GMAT DefaultIB.Element3 = 0;
GMAT DefaultIB.DecrementMass = false;
GMAT DefaultIB.Isp = 300;
GMAT DefaultIB.GravitationalAccel = 9.81;

%----------------------------------------
%---------- Coordinate Systems
%----------------------------------------

Create CoordinateSystem MoonJ2K;
GMAT MoonJ2K.Origin = Luna;
GMAT MoonJ2K.Axes = MJ2000Eq;

Create CoordinateSystem MoonIAU;
GMAT MoonIAU.Origin = Luna;
GMAT MoonIAU.Axes = BodyFixed;

%----------------------------------------
%---------- Solvers
%----------------------------------------

Create DifferentialCorrector DefaultDC;
GMAT DefaultDC.ShowProgress = true;
GMAT DefaultDC.ReportStyle = Normal;
GMAT DefaultDC.ReportFile = 'DifferentialCorrectorDefaultDC.data';
GMAT DefaultDC.MaximumIterations = 75;
GMAT DefaultDC.DerivativeMethod = ForwardDifference;
GMAT DefaultDC.Algorithm = NewtonRaphson;

%----------------------------------------
%---------- Subscribers
%----------------------------------------

Create ReportFile ReportFile1;
GMAT ReportFile1.SolverIterations = Current;
GMAT ReportFile1.UpperLeft = [ 0 0 ];
GMAT ReportFile1.Size = [ 0 0 ];
GMAT ReportFile1.RelativeZOrder = 0;
GMAT ReportFile1.Maximized = false;
GMAT ReportFile1.Filename = 'ReportFile1.txt';
GMAT ReportFile1.Precision = 16;
GMAT ReportFile1.Add = {DefaultSC.UTCGregorian, DefaultSC.MoonIAU.X, DefaultSC.MoonIAU.Y, DefaultSC.MoonIAU.Z, DefaultSC.MoonIAU.VX, DefaultSC.MoonIAU.VY, DefaultSC.MoonIAU.VZ, DefaultSC.MoonIAU.VMAG, DefaultSC.MoonJ2K.X, DefaultSC.MoonJ2K.Y, DefaultSC.MoonJ2K.Z, DefaultSC.MoonJ2K.VX, DefaultSC.MoonJ2K.VY, DefaultSC.MoonJ2K.VZ, DefaultSC.MoonJ2K.VMAG};
GMAT ReportFile1.LeftJustify = On;
GMAT ReportFile1.ZeroFill = Off;
GMAT ReportFile1.FixedWidth = true;
GMAT ReportFile1.Delimiter = ' ';
GMAT ReportFile1.ColumnWidth = 23;
GMAT ReportFile1.WriteReport = true;

%----------------------------------------
%---------- Arrays, Variables, Strings
%----------------------------------------
Create Variable BurnDuration q1 q2 q3 q4 AchieveTime;
GMAT BurnDuration = 600;
GMAT q1 = 0;
GMAT q2 = 0;
GMAT q3 = 0;
GMAT q4 = 0;

%----------------------------------------
%---------- Mission Sequence
%----------------------------------------

BeginMissionSequence;
% Set the thurst direction
GMAT DefaultSC.Q1 = q1;
GMAT DefaultSC.Q2 = q2;
GMAT DefaultSC.Q3 = q3;
GMAT DefaultSC.Q4 = q4;
Target DefaultDC {SolveMode = Solve, ExitMode = DiscardAndContinue, ShowProgressWindow = true};
Vary 'Vary burn duration' DefaultDC(BurnDuration = 600, {Perturbation = 0.5, Lower = 0.0, Upper = 9000, MaxStep = 60, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
Vary 'Vary q1' DefaultDC(q1 = 1, {Perturbation = 0.1, Lower = -1.0, Upper = 1.0, MaxStep = 0.2, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
Vary 'Vary q2' DefaultDC(q2 = 0, {Perturbation = 0.1, Lower = -1.0, Upper = 1.0, MaxStep = 0.2, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
Vary 'Vary q3' DefaultDC(q3 = 0, {Perturbation = 0.1, Lower = -1.0, Upper = 1.0, MaxStep = 0.2, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
Vary 'Vary q4' DefaultDC(q4 = 0, {Perturbation = 0.1, Lower = -1.0, Upper = 1.0, MaxStep = 0.2, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
GMAT DefaultSC.Quaternion = [q1 q2 q3 q4];
Propagate DefaultProp(DefaultSC) {DefaultSC.Luna.Periapsis};
BeginFiniteBurn DefaultFB(DefaultSC);
Propagate DefaultProp(DefaultSC) {DefaultSC.ElapsedSecs = BurnDuration};
EndFiniteBurn DefaultFB(DefaultSC);
Propagate DefaultProp(DefaultSC) {DefaultSC.Luna.Apoapsis};
Achieve 'Achieve SMA' DefaultDC(DefaultSC.Luna.SMA = 1968.0, {Tolerance = 0.1});
%Achieve 'Achieve VMAG' DefaultDC(DefaultSC.MoonJ2K.VMAG = 0.1, {Tolerance = 0.1});
%Achieve 'Achieve Rx' DefaultDC(DefaultSC.MoonJ2K.X = -77.018285, {Tolerance = 0.1});
%Achieve 'Achieve Ry' DefaultDC(DefaultSC.MoonJ2K.Y = 626.836085, {Tolerance = 0.1});
%Achieve 'Achieve Rz' DefaultDC(DefaultSC.MoonJ2K.Z = -1617.498663, {Tolerance = 0.1});
EndTarget;  % For targeter DefaultDC

$$$$


First off: Chris you may already know some/all of my response due to your background. I will attempt to both answer your GMAT question and provide insight to others who may not be as familiar with astrodynamics.

TL;DR Finite burns have a direction, not an attitude. Also, targeting (e.g. via a differential corrector) does best with a square system where the number of variables and constraints are equal.

## Using GMAT

You must set the spacecraft attitude model to Spinner if you want to set its attitude directly. You currently have it set to CoordinateSystemFixed, which as the name implies fixes the spacecraft attitude to align with the selected coordinate system. (and ignores any attitude you specify)

But your actual issue here is much higher up than that implementation detail. I'll explain more in the sections below, but for now I suggest that you not link the burn direction to the spacecraft attitude. That adds complexity to the problem that you don't need.

Instead, vary the Chemical Thruster thrust direction directly. See the Finite Burns section below for reasoning.

1. Set the ChemicalThruster1 coordinate system to MoonJ2K. Make sure you have its C1 set to the desired thrust magnitude.

2. Create two GMAT variables, $$\alpha$$ and $$\beta$$, that you will use to compute the burn direction. Vary those two variables along with the burn direction.

Vary 'Vary burn duration' DefaultDC(BurnDuration = 600, {Perturbation = 0.5, Lower = 0.0, Upper = 9000, MaxStep = 60, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
Vary 'Vary Thrust Alpha' DefaultDC(alpha = 0.0, {Perturbation = 0.1, Lower = 0, Upper = 360.0, MaxStep = 1.0, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
Vary 'Vary Thrust Beta' DefaultDC(beta = 0.0, {Perturbation = 0.1, Lower = 0, Upper = 360.0, MaxStep = 1.0, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});

GMAT 'Compute Thrust Z' ChemicalThruster1.ThrustDirection3 = sin(Deg2Rad(beta));


Optional You should have an equal number of constraints and variables to do targeting. So either add two more constraints, or switch to optimization, or be ok with a non-unique solution. (See Targeting and Optimization section below)

If I take your script and make this change, and also increase thrust magnitude C1 to 200N, the targeter converges in 7 iterations. Increasing the thrust magnitude is important since the default 10N thrust on a 850kg spacecraft won't do much.

## Finite Burns

Finite-burn maneuvers have a duration, a magnitude, and a unit direction. FB maneuvers do not have an attitude. This may seem like nitpicking, but it actually speaks to the crux of the issue here. Attitude describes how the spacecraft is oriented (e.g. with a quaternion), and is specific to the spacecraft (e.g. its thruster configuration, how it was modelled, etc.). In other words, attitude is more than just direction and is therefore an inappropriate/overkill way to parameterize a finite burn during targeting or optimization. Using it forces the targeter/optimizer to contend with variables that do not affect the finite burn, which increases the solver's risk of failure.

So to properly formulate the problem, you should parameterize the maneuver using duration and two angles to represent the burn direction: $$\alpha$$ (in-plane angle) and $$\beta$$ (out-of-plane angle). I assume the thrust magnitude (which GMAT calls C1) is fixed since you did not say anything about it.

## Targeting and Optimization

Targeting is the process by which you vary certain input parameters in order to drive a system to a desired state. Because targeting is basically solving a system of nonlinear equations via linearization, it needs an equal numbers of variables and constraints to work properly (e.g. a square system). However, you are trying to vary 3 variables (duration, $$\alpha$$, $$\beta$$) in order to achieve 1 constraint (semimajor axis). This is an underdetermined system with an infinite number of solutions, so the targeter will either not find a solution, or find a nonunique solution. As long as you go in with eyes wide open, that's fine, but know that targeting is usually done with an equal number of variables and constraints.

Alternatively, you could use an optimizer instead of a targeter. That would allow you to keep your underdetermined parameterization, but you would have to add a cost function. e.g. you could minimize thrust duration. Then you could find at least a locally optimial solution. GMAT comes with several samples scripts of optimization that you can look into.

## The Modified GMAT Script

%General Mission Analysis Tool(GMAT) Script
%Created: 2021-12-08 17:44:08
%Modified: 2022-03-27

%----------------------------------------
%---------- Spacecraft
%----------------------------------------

Create Spacecraft DefaultSC;
GMAT DefaultSC.DateFormat = UTCGregorian;
GMAT DefaultSC.Epoch = '25 Nov 2023 14:11:26.789';
GMAT DefaultSC.CoordinateSystem = MoonJ2K;
GMAT DefaultSC.DisplayStateType = Cartesian;
GMAT DefaultSC.X = -406.8997620000155;
GMAT DefaultSC.Y = -111.8177819999983;
GMAT DefaultSC.Z = -1697.662681000002;
GMAT DefaultSC.VX = 0.620734;
GMAT DefaultSC.VY = 1.591183;
GMAT DefaultSC.VZ = -0.252232;
GMAT DefaultSC.DryMass = 850;
GMAT DefaultSC.Cd = 2.2;
GMAT DefaultSC.Cr = 1.8;
GMAT DefaultSC.DragArea = 15;
GMAT DefaultSC.SRPArea = 1;
GMAT DefaultSC.Tanks = {ChemicalTank1};
GMAT DefaultSC.Thrusters = {ChemicalThruster1};
GMAT DefaultSC.NAIFId = -10000001;
GMAT DefaultSC.NAIFIdReferenceFrame = -9000001;
GMAT DefaultSC.OrbitColor = Red;
GMAT DefaultSC.TargetColor = Teal;
GMAT DefaultSC.OrbitErrorCovariance = [ 1e+70 0 0 0 0 0 ; 0 1e+70 0 0 0 0 ; 0 0 1e+70 0 0 0 ; 0 0 0 1e+70 0 0 ; 0 0 0 0 1e+70 0 ; 0 0 0 0 0 1e+70 ];
GMAT DefaultSC.CdSigma = 1e+70;
GMAT DefaultSC.CrSigma = 1e+70;
GMAT DefaultSC.Id = 'SatId';
GMAT DefaultSC.Attitude = CoordinateSystemFixed;
GMAT DefaultSC.ModelFile = 'aura.3ds';
GMAT DefaultSC.ModelOffsetX = 0;
GMAT DefaultSC.ModelOffsetY = 0;
GMAT DefaultSC.ModelOffsetZ = 0;
GMAT DefaultSC.ModelRotationX = 0;
GMAT DefaultSC.ModelRotationY = 0;
GMAT DefaultSC.ModelRotationZ = 0;
GMAT DefaultSC.ModelScale = 1;
GMAT DefaultSC.AttitudeDisplayStateType = 'Quaternion';
GMAT DefaultSC.AttitudeRateDisplayStateType = 'AngularVelocity';
GMAT DefaultSC.AttitudeCoordinateSystem = EarthMJ2000Eq;
GMAT DefaultSC.EulerAngleSequence = '321';

%----------------------------------------
%---------- Hardware Components
%----------------------------------------

Create ChemicalThruster ChemicalThruster1;
GMAT ChemicalThruster1.CoordinateSystem = MoonJ2K;
GMAT ChemicalThruster1.ThrustDirection1 = 1;
GMAT ChemicalThruster1.ThrustDirection2 = 0;
GMAT ChemicalThruster1.ThrustDirection3 = 0;
GMAT ChemicalThruster1.DutyCycle = 1;
GMAT ChemicalThruster1.ThrustScaleFactor = 1;
GMAT ChemicalThruster1.DecrementMass = false;
GMAT ChemicalThruster1.Tank = {ChemicalTank1};
GMAT ChemicalThruster1.MixRatio = [ 1 ];
GMAT ChemicalThruster1.GravitationalAccel = 9.81;
GMAT ChemicalThruster1.C1 = 200;
GMAT ChemicalThruster1.C2 = 0;
GMAT ChemicalThruster1.C3 = 0;
GMAT ChemicalThruster1.C4 = 0;
GMAT ChemicalThruster1.C5 = 0;
GMAT ChemicalThruster1.C6 = 0;
GMAT ChemicalThruster1.C7 = 0;
GMAT ChemicalThruster1.C8 = 0;
GMAT ChemicalThruster1.C9 = 0;
GMAT ChemicalThruster1.C10 = 0;
GMAT ChemicalThruster1.C11 = 0;
GMAT ChemicalThruster1.C12 = 0;
GMAT ChemicalThruster1.C13 = 0;
GMAT ChemicalThruster1.C14 = 0;
GMAT ChemicalThruster1.C15 = 0;
GMAT ChemicalThruster1.C16 = 0;
GMAT ChemicalThruster1.K1 = 300;
GMAT ChemicalThruster1.K2 = 0;
GMAT ChemicalThruster1.K3 = 0;
GMAT ChemicalThruster1.K4 = 0;
GMAT ChemicalThruster1.K5 = 0;
GMAT ChemicalThruster1.K6 = 0;
GMAT ChemicalThruster1.K7 = 0;
GMAT ChemicalThruster1.K8 = 0;
GMAT ChemicalThruster1.K9 = 0;
GMAT ChemicalThruster1.K10 = 0;
GMAT ChemicalThruster1.K11 = 0;
GMAT ChemicalThruster1.K12 = 0;
GMAT ChemicalThruster1.K13 = 0;
GMAT ChemicalThruster1.K14 = 0;
GMAT ChemicalThruster1.K15 = 0;
GMAT ChemicalThruster1.K16 = 0;

Create ChemicalTank ChemicalTank1;
GMAT ChemicalTank1.AllowNegativeFuelMass = false;
GMAT ChemicalTank1.FuelMass = 756;
GMAT ChemicalTank1.Pressure = 1500;
GMAT ChemicalTank1.Temperature = 20;
GMAT ChemicalTank1.RefTemperature = 20;
GMAT ChemicalTank1.Volume = 0.75;
GMAT ChemicalTank1.FuelDensity = 1260;
GMAT ChemicalTank1.PressureModel = PressureRegulated;

%----------------------------------------
%---------- ForceModels
%----------------------------------------

Create ForceModel DefaultProp_ForceModel;
GMAT DefaultProp_ForceModel.CentralBody = Luna;
GMAT DefaultProp_ForceModel.PointMasses = {Luna};
GMAT DefaultProp_ForceModel.Drag = None;
GMAT DefaultProp_ForceModel.SRP = Off;
GMAT DefaultProp_ForceModel.RelativisticCorrection = Off;

%----------------------------------------
%---------- Propagators
%----------------------------------------

Create Propagator DefaultProp;
GMAT DefaultProp.FM = DefaultProp_ForceModel;
GMAT DefaultProp.Type = RungeKutta89;
GMAT DefaultProp.InitialStepSize = 60;
GMAT DefaultProp.Accuracy = 9.999999999999999e-12;
GMAT DefaultProp.MinStep = 0.001;
GMAT DefaultProp.MaxStep = 2700;
GMAT DefaultProp.MaxStepAttempts = 50;
GMAT DefaultProp.StopIfAccuracyIsViolated = true;

%----------------------------------------
%---------- Burns
%----------------------------------------

Create FiniteBurn DefaultFB;
GMAT DefaultFB.Thrusters = {ChemicalThruster1};
GMAT DefaultFB.ThrottleLogicAlgorithm = 'MaxNumberOfThrusters';

%----------------------------------------
%---------- Coordinate Systems
%----------------------------------------

Create CoordinateSystem MoonJ2K;
GMAT MoonJ2K.Origin = Luna;
GMAT MoonJ2K.Axes = MJ2000Eq;

Create CoordinateSystem MoonIAU;
GMAT MoonIAU.Origin = Luna;
GMAT MoonIAU.Axes = BodyFixed;

%----------------------------------------
%---------- Solvers
%----------------------------------------

Create DifferentialCorrector DefaultDC;
GMAT DefaultDC.ShowProgress = true;
GMAT DefaultDC.ReportStyle = Normal;
GMAT DefaultDC.ReportFile = 'DifferentialCorrectorDefaultDC.data';
GMAT DefaultDC.MaximumIterations = 1000;
GMAT DefaultDC.DerivativeMethod = ForwardDifference;
GMAT DefaultDC.Algorithm = NewtonRaphson;

%----------------------------------------
%---------- Subscribers
%----------------------------------------

Create ReportFile ReportFile1;
GMAT ReportFile1.SolverIterations = Current;
GMAT ReportFile1.UpperLeft = [ 0 0 ];
GMAT ReportFile1.Size = [ 0 0 ];
GMAT ReportFile1.RelativeZOrder = 0;
GMAT ReportFile1.Maximized = false;
GMAT ReportFile1.Filename = 'ReportFile1.txt';
GMAT ReportFile1.Precision = 16;
GMAT ReportFile1.Add = {DefaultSC.UTCGregorian, DefaultSC.MoonIAU.X, DefaultSC.MoonIAU.Y, DefaultSC.MoonIAU.Z, DefaultSC.MoonIAU.VX, DefaultSC.MoonIAU.VY, DefaultSC.MoonIAU.VZ, DefaultSC.MoonIAU.VMAG, DefaultSC.MoonJ2K.X, DefaultSC.MoonJ2K.Y, DefaultSC.MoonJ2K.Z, DefaultSC.MoonJ2K.VX, DefaultSC.MoonJ2K.VY, DefaultSC.MoonJ2K.VZ, DefaultSC.MoonJ2K.VMAG};
GMAT ReportFile1.LeftJustify = On;
GMAT ReportFile1.ZeroFill = Off;
GMAT ReportFile1.FixedWidth = true;
GMAT ReportFile1.Delimiter = ' ';
GMAT ReportFile1.ColumnWidth = 23;
GMAT ReportFile1.WriteReport = true;

Create OpenFramesInterface OpenFrames1;
GMAT OpenFrames1.SolverIterations = Current;
GMAT OpenFrames1.UpperLeft = [ 0.3869047619047619 0.1314285714285714 ];
GMAT OpenFrames1.Size = [ 0.4452380952380952 0.4838095238095238 ];
GMAT OpenFrames1.RelativeZOrder = 71;
GMAT OpenFrames1.Maximized = false;
GMAT OpenFrames1.View = {LunaView1, DefaultSCView1};
GMAT OpenFrames1.CoordinateSystem = MoonJ2K;
GMAT OpenFrames1.DrawObject = [ true true ];
GMAT OpenFrames1.DrawTrajectory = [ true true ];
GMAT OpenFrames1.DrawAxes = [ false false ];
GMAT OpenFrames1.DrawXYPlane = [ false false ];
GMAT OpenFrames1.DrawLabel = [ true true ];
GMAT OpenFrames1.DrawUsePropLabel = [ false false ];
GMAT OpenFrames1.DrawCenterPoint = [ true false ];
GMAT OpenFrames1.DrawEndPoints = [ true false ];
GMAT OpenFrames1.DrawVelocity = [ false false ];
GMAT OpenFrames1.DrawGrid = [ false false ];
GMAT OpenFrames1.DrawLineWidth = [ 2 2 ];
GMAT OpenFrames1.DrawMarkerSize = [ 10 10 ];
GMAT OpenFrames1.DrawFontSize = [ 14 14 ];
GMAT OpenFrames1.Axes = On;
GMAT OpenFrames1.AxesLength = 1;
GMAT OpenFrames1.AxesLabels = On;
GMAT OpenFrames1.FrameLabel = Off;
GMAT OpenFrames1.XYPlane = On;
GMAT OpenFrames1.EclipticPlane = Off;
GMAT OpenFrames1.EnableStars = On;
GMAT OpenFrames1.StarCatalog = 'inp_StarsHYGv3.txt';
GMAT OpenFrames1.StarCount = 40000;
GMAT OpenFrames1.MinStarMag = -2;
GMAT OpenFrames1.MaxStarMag = 6;
GMAT OpenFrames1.MinStarPixels = 1;
GMAT OpenFrames1.MaxStarPixels = 10;
GMAT OpenFrames1.MinStarDimRatio = 0.5;
GMAT OpenFrames1.ShowPlot = true;
GMAT OpenFrames1.ShowToolbar = true;
GMAT OpenFrames1.SolverIterLastN = 1;
GMAT OpenFrames1.ShowVR = false;
GMAT OpenFrames1.PlaybackTimeScale = 3600;
GMAT OpenFrames1.MultisampleAntiAliasing = On;
GMAT OpenFrames1.MSAASamples = 2;
GMAT OpenFrames1.DrawFontPosition = {'Top-Right', 'Top-Right'};

%----------------------------------------
%---------- Arrays, Variables, Strings
%----------------------------------------
Create Variable BurnDuration alpha beta AchieveTime;
GMAT BurnDuration = 600;
GMAT alpha = 0;
GMAT beta = 0;

%----------------------------------------
%---------- User Objects
%----------------------------------------

Create OpenFramesView LunaView1;
GMAT LunaView1.ViewFrame = Luna;
GMAT LunaView1.ViewTrajectory = Off;
GMAT LunaView1.InertialFrame = Off;
GMAT LunaView1.SetDefaultLocation = Off;
GMAT LunaView1.SetCurrentLocation = On;
GMAT LunaView1.CurrentEye = [ 4393.341861034918 -9387.820084998433 -1225.40538297525 ];
GMAT LunaView1.CurrentCenter = [ 0 0 -1.591615728102624e-12 ];
GMAT LunaView1.CurrentUp = [ -0.3526234018341157 -0.2815103104031884 0.8924173248069744 ];
GMAT LunaView1.FOVy = 45;

Create OpenFramesView DefaultSCView1;
GMAT DefaultSCView1.ViewFrame = DefaultSC;
GMAT DefaultSCView1.ViewTrajectory = Off;
GMAT DefaultSCView1.InertialFrame = Off;
GMAT DefaultSCView1.SetDefaultLocation = Off;
GMAT DefaultSCView1.SetCurrentLocation = Off;
GMAT DefaultSCView1.FOVy = 45;

%----------------------------------------
%---------- Mission Sequence
%----------------------------------------

BeginMissionSequence;
Target DefaultDC {SolveMode = Solve, ExitMode = DiscardAndContinue, ShowProgressWindow = true};
Vary 'Vary burn duration' DefaultDC(BurnDuration = 600, {Perturbation = 0.5, Lower = 0.0, Upper = 9000, MaxStep = 60, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
Vary 'Vary Thrust Alpha' DefaultDC(alpha = 0, {Perturbation = 0.1, Lower = 0.0, Upper = 360.0, MaxStep = 1.0, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
Vary 'Vary Thrust Beta' DefaultDC(beta = 0, {Perturbation = 0.1, Lower = 0.0, Upper = 360.0, MaxStep = 1.0, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});

GMAT 'Compute Thrust Z' ChemicalThruster1.ThrustDirection3 = sin(Deg2Rad(beta));

Propagate 'Prop to Perilune' DefaultProp(DefaultSC) {DefaultSC.Luna.Periapsis};
BeginFiniteBurn DefaultFB(DefaultSC);
Propagate 'Prop Burn' DefaultProp(DefaultSC) {DefaultSC.ElapsedSecs = BurnDuration};
EndFiniteBurn DefaultFB(DefaultSC);
Propagate 'Prop to Apolune' DefaultProp(DefaultSC) {DefaultSC.Luna.Apoapsis};
Achieve 'Achieve SMA' DefaultDC(DefaultSC.Luna.SMA = 1968.0, {Tolerance = 0.1});
%Achieve 'Achieve VMAG' DefaultDC(DefaultSC.MoonJ2K.VMAG = 0.1, {Tolerance = 0.1});
%Achieve 'Achieve Rx' DefaultDC(DefaultSC.MoonJ2K.X = -77.018285, {Tolerance = 0.1});
%Achieve 'Achieve Ry' DefaultDC(DefaultSC.MoonJ2K.Y = 626.836085, {Tolerance = 0.1});
%Achieve 'Achieve Rz' DefaultDC(DefaultSC.MoonJ2K.Z = -1617.498663, {Tolerance = 0.1});
EndTarget;  % For targeter DefaultDC
$$$$