SolarSystem.EphemerisSource = 'DE421';
Luna.EquatorialRadius = 1738.2e-30;
Earth.EquatorialRadius = 1738.2e-30;
   

Create CoordinateSystem SunOrigin_EME;
SunOrigin_EME.Origin = Sun;
SunOrigin_EME.Axes = MJ2000Eq;

Create Spacecraft Fermi_cb;  %cb = canon ball
Fermi_cb.DateFormat         = UTCGregorian;
Fermi_cb.Epoch              = '01 Mar 2013 00:00:00.000';
Fermi_cb.CoordinateSystem = SunOrigin_EME;
Fermi_cb.DisplayStateType = Keplerian;
Fermi_cb.SMA = 145864620
Fermi_cb.ECC = 0.001
Fermi_cb.INC = 0;
Fermi_cb.RAAN = 0;
Fermi_cb.AOP = 0;
Fermi_cb.TA = 0;
Fermi_cb.DryMass            = 1000
Fermi_cb.Cr                 = 1
Fermi_cb.SRPArea            = 500


Create Spacecraft Fermi_NPlate;
Fermi_NPlate.DateFormat         = UTCGregorian;
Fermi_NPlate.Epoch              = '01 Mar 2013 00:00:00.000';
Fermi_NPlate.CoordinateSystem = SunOrigin_EME;
Fermi_NPlate.DisplayStateType = Keplerian;
Fermi_NPlate.SMA = 145864620
Fermi_NPlate.ECC = 0.001
Fermi_NPlate.INC = 0;
Fermi_NPlate.RAAN = 0;
Fermi_NPlate.AOP = 0;
Fermi_NPlate.TA = 0;
Fermi_NPlate.DryMass            = 1000
Fermi_NPlate.Cr                 = 1e70  %this value not used
Fermi_NPlate.SRPArea            = 1e70  %this value not used
Fermi_NPlate.Attitude 			 = NadirPointing;
Fermi_NPlate.AttitudeReferenceBody = Sun;
Fermi_NPlate.AttitudeConstraintType = 'OrbitNormal';
Fermi_NPlate.BodyAlignmentVectorX = 1;
Fermi_NPlate.BodyAlignmentVectorY = 0;
Fermi_NPlate.BodyAlignmentVectorZ = 0;
Fermi_NPlate.BodyConstraintVectorX = 0;
Fermi_NPlate.BodyConstraintVectorY = 0;
Fermi_NPlate.BodyConstraintVectorZ = 1;
%above attitude should point X body axis towards Sun.  

Create Spacecraft Fermi_SPAD;
Fermi_SPAD.DateFormat         = UTCGregorian;
Fermi_SPAD.Epoch              = '01 Mar 2013 00:00:00.000';
Fermi_SPAD.CoordinateSystem = SunOrigin_EME;
Fermi_SPAD.DisplayStateType = Keplerian;
Fermi_SPAD.SMA = 145864620
Fermi_SPAD.ECC = 0.001
Fermi_SPAD.INC = 0;
Fermi_SPAD.RAAN = 0;
Fermi_SPAD.AOP = 0;
Fermi_SPAD.TA = 0;
Fermi_SPAD.DryMass            = 1000
Fermi_SPAD.Cr                 = 1e70  %this value not used
Fermi_SPAD.SRPArea            = 1e70  %this value not used
Fermi_SPAD.SPADSRPFile = 'C:\Users\dcooley\Documents\GMAT\Nav\Tickets\GMT-7004_NPlate_Compare\OneHalfDeg_500m2.spo';
Fermi_SPAD.SPADSRPScaleFactor = 1;
Fermi_SPAD.SPADSRPInterpolationMethod = Bicubic


Create Plate BP1 BP2 
% 
BP1.Type = FixedInBody
BP1.PlateNormal = [1.0, 0.0, 0.0] 					
BP1.LitFraction = 1.0                          				
BP1.AreaCoefficient = 1.0 					
BP1.Area = 500
BP1.SpecularFraction = 0
BP1.DiffuseFraction = 0


BP2.Type = FixedInBody
BP2.PlateNormal = [-1.0, 0.0, 0.0] 					
BP2.LitFraction = 1.0                          				
BP2.AreaCoefficient = 1.0			
BP2.Area = 1e70    %shouldnt be used since Sun not shining on this panel
BP2.SpecularFraction = 0.5;
BP2.DiffuseFraction = 0.1;


Fermi_NPlate.AddPlates    = {BP1, BP2}; 

%The Fermi_NPlate S/c attitude points the X axis towards the Sun.  Fermi s/c uses NPlate model.
%The orbit obtained should be identical to the SimFermi S/C canon-ball
%SRP model, Cr=2.0, with a cross sectional area of 500 m^2
%and a mass of 1000 kg.



%
%   Force model and propagator
%

Create ForceModel FM_NPlate FM_cb FM_SPAD
Create Propagator Prop_NPlate Prop_cb Prop_SPAD;

FM_cb.CentralBody                       = Sun;
FM_cb.PointMasses                       = {Sun};
FM_cb.SRP                               = On;
FM_cb.SRP.SRPModel 			 			     = Spherical
FM_cb.Drag.AtmosphereModel              = 'None'
FM_cb.ErrorControl                      = 'None';
FM_cb.SRP.Flux 								  = 1358;
FM_cb.SRP.Nominal_Sun 						  = 149597870.691;

Prop_cb.FM                              = FM_cb;
Prop_cb.Type                            = RungeKutta89;
Prop_cb.InitialStepSize                 = 5;  %set small to make sure attitude is more accurate
Prop_cb.Accuracy                        = 1e-13;
Prop_cb.MinStep                         = 0;
Prop_cb.MaxStep                         = 60;
Prop_cb.MaxStepAttempts                 = 50;

FM_NPlate.CentralBody                       = Sun;
FM_NPlate.PointMasses                       = {Sun};
FM_NPlate.SRP                               = On;
FM_NPlate.SRP.SRPModel 			 			    = NPlate
FM_NPlate.Drag.AtmosphereModel              = 'None'
FM_NPlate.ErrorControl                      = 'None';
FM_NPlate.SRP.Flux 								 = 1358;
FM_NPlate.SRP.Nominal_Sun 						 = 149597870.691;


Prop_NPlate.FM                              = FM_NPlate;
Prop_NPlate.Type                            = RungeKutta89;
Prop_NPlate.InitialStepSize                 = 5;  %set small to make sure attitude is more accurate
Prop_NPlate.Accuracy                        = 1e-13;
Prop_NPlate.MinStep                         = 0;
Prop_NPlate.MaxStep                         = 60;
Prop_NPlate.MaxStepAttempts                 = 50;

FM_SPAD.CentralBody                       = Sun;
FM_SPAD.PointMasses                       = {Sun};
FM_SPAD.SRP                               = On;
FM_SPAD.SRP.SRPModel 			 			   = 'SPADFile' 
FM_SPAD.Drag.AtmosphereModel              = 'None'
FM_SPAD.ErrorControl                      = 'None';
FM_SPAD.SRP.Flux 									= 1358;
FM_SPAD.SRP.Nominal_Sun 						= 149597870.691;

Prop_SPAD.FM                              = FM_SPAD;
Prop_SPAD.Type                            = RungeKutta89;
Prop_SPAD.InitialStepSize                 = 5;  %set small to make sure attitude is more accurate
Prop_SPAD.Accuracy                        = 1e-13;
Prop_SPAD.MinStep                         = 0;
Prop_SPAD.MaxStep                         = 60;
Prop_SPAD.MaxStepAttempts                 = 50;




Create ReportFile rf
Create Variable Pos_Error_Nplate_meters Pos_Error_SPAD_meters

BeginMissionSequence

% 
Propagate Prop_cb(Fermi_cb) {Fermi_cb.ElapsedSecs = 86400};
Propagate Prop_NPlate(Fermi_NPlate) {Fermi_NPlate.ElapsedSecs = 86400};
Propagate Prop_SPAD(Fermi_SPAD) {Fermi_SPAD.ElapsedSecs = 86400};
% 
% 
Pos_Error_Nplate_meters = 1000*sqrt ( (Fermi_cb.X-Fermi_NPlate.X)^2  + (Fermi_cb.Y-Fermi_NPlate.Y)^2  + (Fermi_cb.Z-Fermi_NPlate.Z)^2  );

Pos_Error_SPAD_meters = 1000*sqrt ( (Fermi_cb.X-Fermi_SPAD.X)^2  + (Fermi_cb.Y-Fermi_SPAD.Y)^2  + (Fermi_cb.Z-Fermi_SPAD.Z)^2  );


Report rf Pos_Error_Nplate_meters Pos_Error_SPAD_meters
%