In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import xarray as xr
import datetime
In [2]:
mjd_unixtimestamp_offset = 10587.5
seconds_in_day = 3600 * 24
c = 299792458
def mjd2unixtimestamp(m):
return (m - mjd_unixtimestamp_offset) * seconds_in_day
def load_gmat_report(path):
GMAT_COLS = 7
data = np.fromfile(path, sep=' ').reshape((-1, GMAT_COLS))
t = [np.datetime64(datetime.datetime.utcfromtimestamp(s)) for s in mjd2unixtimestamp(data[:,0])]
xr_data = xr.DataArray(data[:,1:].reshape((-1, 2 ,3)), dims = ('time', 'rv', 'xyz'),\
coords = {'time' : t, 'rv' : ['r', 'v'], 'xyz' : ['x', 'y', 'z']})
return xr_data
def doppler_ppb(data):
return -1e12 / c * ((data.sel(rv = 'r') * data.sel(rv = 'v')).sum('xyz') / np.sqrt((data.sel(rv = 'r')**2).sum('xyz'))).drop('rv')
In [3]:
ecef = load_gmat_report('ecef_report.txt')
ea4gpz = load_gmat_report('ea4gpz_report.txt')
In [4]:
plt.figure(figsize = (12,6), facecolor = 'w')
np.sqrt((ea4gpz.sel(rv = 'r')**2).sum('xyz')).plot()
plt.ylabel('Distance (km)')
plt.xlabel('UTC time')
plt.title('Distance from Es\'hail 2 to EA4GPZ');
In [5]:
plt.figure(figsize = (12,6), facecolor = 'w')
doppler_ppb(ea4gpz).plot()
plt.ylabel('Doppler (ppb)')
plt.xlabel('UTC time')
plt.title('Doppler from Es\'hail 2 to EA4GPZ');
In [6]:
plt.figure(figsize = (12,6), facecolor = 'w')
ecef.sel(rv = 'r', xyz = 'z').plot()
plt.ylabel('ECEF z (km)')
plt.xlabel('UTC time')
plt.title('Es\'hail 2 position: ECEF z coordinate');
In [7]:
plt.figure(figsize = (12,6), facecolor = 'w')
np.rad2deg(np.arctan(ecef.sel(rv = 'r', xyz = 'y')/ecef.sel(rv = 'r', xyz = 'x'))).plot()
plt.ylabel('Longitude (deg)')
plt.xlabel('UTC time')
plt.title('Es\'hail 2 position: longitude');
In [8]:
plt.figure(figsize = (12,6), facecolor = 'w')
ax1 = plt.gca()
ax2 = ax1.twinx()
ax1.plot(ecef.coords['time'], ecef.sel(rv = 'r', xyz = 'x'), color = 'C0', label = 'ECEF x')
ax2.plot(ecef.coords['time'], ecef.sel(rv = 'r', xyz = 'y'), color = 'C1', label = 'ECEF y')
ax1.set_ylabel('ECEF x (km)', color = 'C0')
ax2.set_ylabel('ECEF y (km)', color = 'C1')
ax1.set_xlabel('UTC time')
plt.title('Es\'hail 2 position: ECEF x and y coordinates');
In [12]:
plt.figure(figsize = (12,6), facecolor = 'w')
np.sqrt((ecef.sel(rv = 'r')**2).sum('xyz')).plot()
plt.ylabel('ECEF radius (km)')
plt.xlabel('UTC time')
plt.title('Es\'hail 2 position: ECEF radius coordinate');