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

Create Spacecraft LISA1;
GMAT LISA1.DateFormat = TAIModJulian;
GMAT LISA1.Epoch = 21545.000000000;
GMAT LISA1.CoordinateSystem = SunMJ2000Ec;
GMAT LISA1.StateType = Keplerian;
LISA1.AnomalyType = MA;
GMAT LISA1.DryMass = 850;
GMAT LISA1.Cd = 2.2;
GMAT LISA1.Cr = 1.8;
GMAT LISA1.DragArea = 15;
GMAT LISA1.SRPArea = 1;

Create Spacecraft LISA2;
GMAT LISA2.DateFormat = TAIModJulian;
GMAT LISA2.Epoch = 21545.000000000;
GMAT LISA2.CoordinateSystem = SunMJ2000Ec;
GMAT LISA2.StateType = Keplerian;
LISA2.AnomalyType = MA;
GMAT LISA2.DryMass = 850;
GMAT LISA2.Cd = 2.2;
GMAT LISA2.Cr = 1.8;
GMAT LISA2.DragArea = 15;
GMAT LISA2.SRPArea = 1;

Create Spacecraft LISA3;
GMAT LISA3.DateFormat = TAIModJulian;
GMAT LISA3.Epoch = 21545.000000000;
GMAT LISA3.CoordinateSystem = SunMJ2000Ec;
GMAT LISA3.StateType = Keplerian;
LISA3.AnomalyType = MA;
GMAT LISA3.DryMass = 850;
GMAT LISA3.Cd = 2.2;
GMAT LISA3.Cr = 1.8;
GMAT LISA3.DragArea = 15;
GMAT LISA3.SRPArea = 1;

Create Spacecraft EarthSat;
GMAT EarthSat.DateFormat = TAIModJulian;
GMAT EarthSat.Epoch = 21545.000000000;
GMAT EarthSat.CoordinateSystem = EarthMJ2000Eq;
GMAT EarthSat.StateType = Cartesian;
GMAT EarthSat.DryMass = 850;
GMAT EarthSat.Cd = 2.2;
GMAT EarthSat.Cr = 1.8;
GMAT EarthSat.DragArea = 15;
GMAT EarthSat.SRPArea = 1;


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

Create ForceModel EarthProp_ForceModel;
GMAT EarthProp_ForceModel.CentralBody = Earth;
GMAT EarthProp_ForceModel.PrimaryBodies = {Earth};
GMAT EarthProp_ForceModel.Drag = None;
GMAT EarthProp_ForceModel.SRP = Off;
GMAT EarthProp_ForceModel.ErrorControl = RSSStep;
GMAT EarthProp_ForceModel.Gravity.Earth.Degree = 4;
GMAT EarthProp_ForceModel.Gravity.Earth.Order = 4;
GMAT EarthProp_ForceModel.Gravity.Earth.PotentialFile = JGM2.grv;

Create ForceModel SunProp_ForceModel;
GMAT SunProp_ForceModel.CentralBody = Sun;
GMAT SunProp_ForceModel.PointMasses = {Sun, Earth, Luna, Jupiter};
GMAT SunProp_ForceModel.Drag = None;
GMAT SunProp_ForceModel.SRP = Off;
GMAT SunProp_ForceModel.ErrorControl = RSSStep;

SolarSystem.Ephemeris = {SLP};

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

Create Propagator EarthProp;
GMAT EarthProp.FM = EarthProp_ForceModel;
GMAT EarthProp.Type = RungeKutta89;
GMAT EarthProp.InitialStepSize = 60;
GMAT EarthProp.Accuracy = 9.999999999999999e-012;
GMAT EarthProp.MinStep = 0.001;
GMAT EarthProp.MaxStep = 2700;
GMAT EarthProp.MaxStepAttempts = 50;

Create Propagator SunProp;
GMAT SunProp.FM = SunProp_ForceModel;
GMAT SunProp.Type = PrinceDormand78;
GMAT SunProp.InitialStepSize = 45000;
GMAT SunProp.Accuracy = 1e-010;
GMAT SunProp.MinStep = 0.001;
GMAT SunProp.MaxStep = 160000;
GMAT SunProp.MaxStepAttempts = 50; 


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

Create ImpulsiveBurn DefaultIB;
GMAT DefaultIB.Origin = Earth;
GMAT DefaultIB.Axes = VNB;
GMAT DefaultIB.VectorFormat = Cartesian;
GMAT DefaultIB.Element1 = 0;
GMAT DefaultIB.Element2 = 0;
GMAT DefaultIB.Element3 = 0;




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

Create CoordinateSystem EarthMJ2000Eq;
GMAT EarthMJ2000Eq.Origin = Earth;
GMAT EarthMJ2000Eq.Axes = MJ2000Eq;
GMAT EarthMJ2000Eq.UpdateInterval = 60;
GMAT EarthMJ2000Eq.OverrideOriginInterval = false;

Create CoordinateSystem EarthMJ2000Ec;
GMAT EarthMJ2000Ec.Origin = Earth;
GMAT EarthMJ2000Ec.Axes = MJ2000Ec;
GMAT EarthMJ2000Ec.UpdateInterval = 60;
GMAT EarthMJ2000Ec.OverrideOriginInterval = false;

Create CoordinateSystem EarthFixed;
GMAT EarthFixed.Origin = Earth;
GMAT EarthFixed.Axes = BodyFixed;
GMAT EarthFixed.UpdateInterval = 60;
GMAT EarthFixed.OverrideOriginInterval = false;

Create CoordinateSystem SunMJ2000Ec;
GMAT SunMJ2000Ec.Origin = Sun;
GMAT SunMJ2000Ec.Axes = MJ2000Ec;
GMAT SunMJ2000Ec.UpdateInterval = 60;
GMAT SunMJ2000Ec.OverrideOriginInterval = false;


%----------------------------------------
%---------- Plots and Reports
%----------------------------------------

Create OpenGLPlot SolarSystemView;
GMAT SolarSystemView.Add = {LISA1 Earth, Sun, Venus, Mars};
GMAT SolarSystemView.OrbitColor = [ 255 32768 4227327 ];
GMAT SolarSystemView.CoordinateSystem = SunMJ2000Ec;
GMAT SolarSystemView.ViewPointRef = Sun;
GMAT SolarSystemView.ViewPointRefVector = [ 0 0 0];
GMAT SolarSystemView.ViewPointVector = Vector;
GMAT SolarSystemView.ViewPointVectorVector = [ 0 0 30000];
GMAT SolarSystemView.ViewDirection = Sun;
GMAT SolarSystemView.ViewDirectionVector = [ 0 0 -1];
GMAT SolarSystemView.ViewScaleFactor = 10000;
GMAT SolarSystemView.FixedFovAngle = 45;
GMAT SolarSystemView.ViewUpCoordinateSystem = SunMJ2000Ec;
GMAT SolarSystemView.ViewUpAxis = X; 
GMAT SolarSystemView.CelestialPlane = Off;
GMAT SolarSystemView.XYPlane = Off;
GMAT SolarSystemView.WireFrame = Off;
GMAT SolarSystemView.TargetStatus = Off;
GMAT SolarSystemView.Axes = On;
GMAT SolarSystemView.Grid = On;
GMAT SolarSystemView.EarthSunLines = Off;
GMAT SolarSystemView.UseInitialView = On;
GMAT SolarSystemView.PerspectiveMode = Off;
GMAT SolarSystemView.UseFixedFov = Off;
GMAT SolarSystemView.DataCollectFrequency = 1;
GMAT SolarSystemView.UpdatePlotFrequency = 50;
GMAT SolarSystemView.NumPointsToRedraw = 0;
GMAT SolarSystemView.ShowPlot = true;

Create ReportFile LISAGeometry;
LISAGeometry.Filename = LISAGeometry;



%----------------------------------------
%---------- Parameters
%----------------------------------------

Create Variable pi SMA ECC INC AOP RAAN MA 
Create Variable d thetaE EarthAngle I
Create Variable S1 S2 S3 FinalEpoch;
Create Array Side1[3,1] Side2[3,1] Side3[3,1] ;

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

%  Define constants and initializations
pi = 3.141592653589793;
d = 5000000;
SMA = EarthSat.Sun.SMA;
ECC = d/(2*sqrt(3)*SMA);
INC = d/2/SMA*180/pi;
AOP = 90;
MA  = 90;
EarthAngle = 180; 

%  Set up spacecraft initial conditions
thetaE = atan2(EarthSat.SunMJ2000Ec.Y,EarthSat.SunMJ2000Ec.X)*180/pi; 

If thetaE < 0
   thetaE = thetaE + 360;
EndIf;
RAAN = thetaE - AOP - MA - EarthAngle;


LISA1.SMA = SMA;
LISA1.ECC = ECC;
LISA1.INC = INC;
LISA1.AOP = AOP;
LISA1.RAAN = RAAN;
LISA1.MA   = MA; 

LISA2.SMA = SMA;
LISA2.ECC = ECC;
LISA2.INC = INC;
LISA2.AOP = AOP;
LISA2.RAAN = RAAN + 120;
LISA2.MA   = MA - 120;

LISA3.SMA = SMA;
LISA3.ECC = ECC;
LISA3.INC = INC;
LISA3.AOP = AOP;
LISA3.RAAN = RAAN - 120;
LISA3.MA   = MA + 120;

LISA2 = LISA1;


%  ---------------------------------------
%  ------ Propagation Type One -----------
%  ---------------------------------------
%  First Propagate for about a year, using 
%  a stopping condition of about 3 days, 
%  using a while loop

While LISA1.ElapsedDays < 365

   Propagate SunProp(LISA1){LISA1.ElapsedSecs = 250000}; 

EndWhile

Report  LISAGeometry  LISA1.A1ModJulian LISA1.X LISA1.Y LISA1.Z
FinalEpoch = LISA1.A1ModJulian;



%  ---------------------------------------
%  ------ Propagation Type Two -----------
%  ---------------------------------------
LISA1 = LISA2;

While LISA1.ElapsedDays < 365

   Propagate SunProp(LISA1); 

EndWhile

Propagate SunProp(LISA1){LISA1.A1ModJulian = FinalEpoch} ; 
Report  LISAGeometry  LISA1.A1ModJulian LISA1.X LISA1.Y LISA1.Z



%  ---------------------------------------
%  ------ Propagation Type Three ---------
%  ---------------------------------------
LISA1 = LISA2
Propagate SunProp(LISA1){LISA1.A1ModJulian = FinalEpoch} ; 
Report  LISAGeometry  LISA1.A1ModJulian LISA1.X LISA1.Y LISA1.Z