<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 DSLWP-B orbital parameter analysis
In this notebook we analyse the Keplerian orbital parameters derived from the DSLWP-B tracking files published by Wei Mingchuan BG2BHC using GMAT. The ECEF cartesian state is loaded from the first line of the tracking file and then the orbit is propagated and some Keplerian state parameters are plotted to compare the elliptical orbits in different tracking files.
Step1:  Set this to the path of your GMAT executable
Step3:  The GMAT script template contains fields ready to be filled in using Python's format() function.
Step4:  Conversion between UNIX timestamp (used by the tracking files) and GMAT Modified Julian Day.
Step5:  Utility function to load the first row from a tracking file and to load the Keplerian state report from GMAT.
Step6:  Keys for each of the columns in the orbit (Keplerian state) report.
Step7:  The function below takes the data from a tracking file, generates a GMAT script and executes it. GMAT is closed automatically after the script has run unless do_not_close is set to True. This can be useful to examine the simulation output in more detail.
Step8:  Load the cartesian state from each tracking file, propagate the orbit and write a Keplerian state report using GMAT.
Step9:  Plot the orbital parameters which are vary significantly between different tracking files.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
Explanation: DSLWP-B orbital parameter analysis
In this notebook we analyse the Keplerian orbital parameters derived from the DSLWP-B tracking files published by Wei Mingchuan BG2BHC using GMAT. The ECEF cartesian state is loaded from the first line of the tracking file and then the orbit is propagated and some Keplerian state parameters are plotted to compare the elliptical orbits in different tracking files.
End of explanation
GMAT_PATH = '/home/daniel/GMAT/R2018a/bin/GMAT-R2018a'
import numpy as np
import matplotlib.pyplot as plt
from astropy.time import Time
import subprocess
# Larger figure size
fig_size = [10, 6]
plt.rcParams['figure.figsize'] = fig_size
Explanation: Set this to the path of your GMAT executable:
End of explanation
gmat_script_template = 
%----------------------------------------
%---------- Spacecraft
%----------------------------------------
Create Spacecraft DSLWP_B;
DSLWP_B.DateFormat = UTCModJulian;
DSLWP_B.Epoch = '{epoch}';
DSLWP_B.CoordinateSystem = EarthFixed;
DSLWP_B.DisplayStateType = Cartesian;
DSLWP_B.X = {x};
DSLWP_B.Y = {y};
DSLWP_B.Z = {z};
DSLWP_B.VX = {vx};
DSLWP_B.VY = {vy};
DSLWP_B.VZ = {vz};
DSLWP_B.DryMass = 45;
DSLWP_B.DragArea = 0.25;
DSLWP_B.SRPArea = 0.25;
%----------------------------------------
%---------- ForceModels
%----------------------------------------
Create ForceModel LunaProp_ForceModel;
LunaProp_ForceModel.CentralBody = Luna;
LunaProp_ForceModel.PrimaryBodies = {{Luna}};
LunaProp_ForceModel.PointMasses = {{Earth, Jupiter, Mars, Neptune, Saturn, Sun, Uranus, Venus}};
LunaProp_ForceModel.Drag = None;
LunaProp_ForceModel.SRP = On;
LunaProp_ForceModel.RelativisticCorrection = On;
LunaProp_ForceModel.ErrorControl = RSSStep;
LunaProp_ForceModel.GravityField.Luna.Degree = 10;
LunaProp_ForceModel.GravityField.Luna.Order = 10;
LunaProp_ForceModel.GravityField.Luna.StmLimit = 100;
LunaProp_ForceModel.GravityField.Luna.PotentialFile = 'LP165P.cof';
LunaProp_ForceModel.GravityField.Luna.TideModel = 'None';
%----------------------------------------
%---------- Propagators
%----------------------------------------
Create Propagator LunaProp;
LunaProp.FM = LunaProp_ForceModel;
LunaProp.Type = PrinceDormand78;
LunaProp.InitialStepSize = 1;
LunaProp.Accuracy = 1e-13;
LunaProp.MinStep = 0;
LunaProp.MaxStep = 600;
%----------------------------------------
%---------- Coordinate Systems
%----------------------------------------
Create CoordinateSystem LunaInertial;
LunaInertial.Origin = Luna;
LunaInertial.Axes = BodyInertial;
%----------------------------------------
%---------- Subscribers
%----------------------------------------
Create OrbitView LunaOrbitView;
GMAT LunaOrbitView.SolverIterations = None;
GMAT LunaOrbitView.UpperLeft = [ 0.1801470588235294 0.04190751445086705 ];
GMAT LunaOrbitView.Size = [ 0.9926470588235294 0.9552023121387283 ];
GMAT LunaOrbitView.RelativeZOrder = 126;
GMAT LunaOrbitView.Maximized = true;
GMAT LunaOrbitView.Add = {{DSLWP_B, Earth, Luna, Sun}};
GMAT LunaOrbitView.CoordinateSystem = LunaInertial;
GMAT LunaOrbitView.DrawObject = [ true true true true ];
GMAT LunaOrbitView.DataCollectFrequency = 1;
GMAT LunaOrbitView.UpdatePlotFrequency = 50;
GMAT LunaOrbitView.NumPointsToRedraw = 0;
GMAT LunaOrbitView.ShowPlot = true;
GMAT LunaOrbitView.MaxPlotPoints = 20000;
GMAT LunaOrbitView.ShowLabels = true;
GMAT LunaOrbitView.ViewPointReference = Luna;
GMAT LunaOrbitView.ViewPointVector = [ 30000 0 0 ];
GMAT LunaOrbitView.ViewDirection = Luna;
GMAT LunaOrbitView.ViewScaleFactor = 1;
GMAT LunaOrbitView.ViewUpCoordinateSystem = LunaInertial;
GMAT LunaOrbitView.ViewUpAxis = Z;
GMAT LunaOrbitView.EclipticPlane = Off;
GMAT LunaOrbitView.XYPlane = On;
GMAT LunaOrbitView.WireFrame = Off;
GMAT LunaOrbitView.Axes = On;
GMAT LunaOrbitView.Grid = Off;
GMAT LunaOrbitView.SunLine = Off;
GMAT LunaOrbitView.UseInitialView = On;
GMAT LunaOrbitView.StarCount = 7000;
GMAT LunaOrbitView.EnableStars = On;
GMAT LunaOrbitView.EnableConstellations = Off;
Create ReportFile OrbitReport;
OrbitReport.Filename = '/home/daniel/jupyter_notebooks/dslwp/OrbitReport_{label}.txt';
OrbitReport.Add = {{DSLWP_B.UTCModJulian, DSLWP_B.Luna.SMA, DSLWP_B.Luna.ECC, DSLWP_B.LunaInertial.INC, DSLWP_B.LunaInertial.RAAN, DSLWP_B.LunaInertial.AOP, DSLWP_B.Luna.MA, DSLWP_B.Luna.TA}};
OrbitReport.WriteHeaders = false;
OrbitReport.WriteReport = true;
%----------------------------------------
%---------- Mission Sequence
%----------------------------------------
BeginMissionSequence;
Toggle OrbitReport Off
If DSLWP_B.UTCModJulian <= {start}
    Propagate LunaProp(DSLWP_B) {{DSLWP_B.UTCModJulian = {start}}}
Else
    Propagate BackProp LunaProp(DSLWP_B) {{DSLWP_B.UTCModJulian = {start}}}
EndIf
Toggle OrbitReport On
Propagate LunaProp(DSLWP_B) {{DSLWP_B.UTCModJulian = {end}}}
Explanation: The GMAT script template contains fields ready to be filled in using Python's format() function.
End of explanation
mjd_unixtimestamp_offset = 10587.5
seconds_in_day = 3600 * 24
def mjd2unixtimestamp(m):
    return (m - mjd_unixtimestamp_offset) * seconds_in_day
def unixtimestamp2mjd(u):
    return u / seconds_in_day + mjd_unixtimestamp_offset
unixtimestamp2mjd(1528607994)
Explanation: Conversion between UNIX timestamp (used by the tracking files) and GMAT Modified Julian Day.
End of explanation
def load_tracking_file(path):
    ncols = 7
    data = np.fromfile(path, sep=' ', count=ncols)
    return data
def load_orbit_file(path):
    ncols = 8
    data = np.fromfile(path, sep=' ')
    return data.reshape((data.size // ncols, ncols))
Explanation: Utility function to load the first row from a tracking file and to load the Keplerian state report from GMAT.
End of explanation
utc = 0
sma = 1
ecc = 2
inc = 3
raan = 4
aop = 5
ma = 6
ta = 7
Explanation: Keys for each of the columns in the orbit (Keplerian state) report.
End of explanation
def gmat_propagate_tracking(track, start, end, label = '', do_not_close = False):
    data = {'label' : label, 'start' : start, 'end' : end}
    data['epoch'] = unixtimestamp2mjd(track[0])
    data['x'], data['y'], data['z'] = track[1:4]
    data['vx'], data['vy'], data['vz'] = track[4:7]
    SCRIPT_PATH = '/tmp/gmat.script'
    with open(SCRIPT_PATH, 'w') as f:
        f.write(gmat_script_template.format(**data))
    subprocess.call([GMAT_PATH, '-r', SCRIPT_PATH] + (['-x'] if not do_not_close else []))
Explanation: The function below takes the data from a tracking file, generates a GMAT script and executes it. GMAT is closed automatically after the script has run unless do_not_close is set to True. This can be useful to examine the simulation output in more detail.
End of explanation
#parts = ['20180526', '20180528', '20180529', '20180531', '20180601', '20180602', '20180603', '20180607', '20180609']
#parts = ['20180602', '20180603', '20180607', '20180609', '20180610', '20180615', '20180619', '20180622']
#parts = ['20180610', '20180615', '20180619', '20180622']
parts = ['20180629', '20180714', '20180727a']
parts = ['20180727a']
parts = ['20180622']
parts = ['20180727a', '20180803', '20180812', '20180814', '20180816', '20180818']
parts = ['20180812', '20180912', '20180914', '20180916', '20180930', '20181004', '20181006']
parts = ['20180812', '20180912', '20180914', '20180916', '20180930', '20181004', '20181006']
parts = ['20181006', '20181008', '20181010', '20181013', '20181015', '20181017', '20181019', '20181021']
parts = ['20181019']
parts = ['20190317']
parts = ['20190426']
parts = ['20190520']
parts = ['20190603']
parts = ['20190630']
for part in parts:
    tracking = load_tracking_file('tracking_files/program_tracking_dslwp-b_{}.txt'.format(part))
    gmat_propagate_tracking(tracking, start = '28560', end = '28570', label = part)
Explanation: Load the cartesian state from each tracking file, propagate the orbit and write a Keplerian state report using GMAT.
End of explanation
fig1 = plt.figure(figsize = [15,8], facecolor='w')
fig2 = plt.figure(figsize = [15,8], facecolor='w')
fig3 = plt.figure(figsize = [15,8], facecolor='w')
fig4 = plt.figure(figsize = [15,8], facecolor='w')
fig5 = plt.figure(figsize = [15,8], facecolor='w')
fig6 = plt.figure(figsize = [15,8], facecolor='w')
sub1 = fig1.add_subplot(111)
sub2 = fig2.add_subplot(111)
sub3 = fig3.add_subplot(111)
sub4 = fig4.add_subplot(111)
sub5 = fig5.add_subplot(111)
sub6 = fig6.add_subplot(111)
for part in parts:
    orbit = load_orbit_file('OrbitReport_{}.txt'.format(part))
    t = Time(mjd2unixtimestamp(orbit[:,utc]), format='unix')
    sub1.plot(t.datetime, orbit[:,sma])
    sub2.plot(t.datetime, orbit[:,ma])
    sub3.plot(t.datetime, orbit[:,ecc])
    sub4.plot(t.datetime, orbit[:,aop])
    sub5.plot(t.datetime, orbit[:,inc])
    sub6.plot(orbit[:,utc], orbit[:,raan])
sub1.legend(parts)
sub2.legend(parts)
sub3.legend(parts)
sub4.legend(parts)
sub5.legend(parts)
sub6.legend(parts)
sub1.set_xlabel('UTC time')
sub2.set_xlabel('UTC time')
sub3.set_xlabel('UTC time')
sub4.set_xlabel('UTC time')
sub5.set_xlabel('UTC time')
sub6.set_xlabel('UTC time')
sub1.set_ylabel('SMA (km)')
sub2.set_ylabel('MA (deg)')
sub3.set_ylabel('ECC')
sub4.set_ylabel('AOP (deg)')
sub5.set_ylabel('INC (deg)')
sub6.set_ylabel('RAAN (deg)')
sub1.set_title('Semi-major axis')
sub2.set_title('Mean anomaly')
sub3.set_title('Eccentricity')
sub4.set_title('Argument of periapsis')
sub5.set_title('Inclination')
sub6.set_title('Right ascension of ascending node');
plt.figure(figsize = [15,8], facecolor='w')
for part in parts:
    orbit = load_orbit_file('OrbitReport_{}.txt'.format(part))
    t = Time(mjd2unixtimestamp(orbit[:,utc]), format='unix')
    plt.plot(t[:200].datetime, orbit[:200,ma])
plt.legend(parts)
plt.xlabel('UTC time')
plt.ylabel('MA (deg)')
plt.title('Mean anomaly');
Explanation: Plot the orbital parameters which are vary significantly between different tracking files.
End of explanation
<END_TASK>