In [1]:

    

import tsyganenko as tsyg
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from spacepy import coordinates as coord
import spacepy.time as spt
from spacepy.time import Ticktock
import datetime as dt
from mpl_toolkits.mplot3d import Axes3D
import sys
Re = 6371
earth_radius_ax = 1.5*Re #km
#adding the year data here so I don't have to crush my github repo
#pathname = '../../data-se3-path-planner/yearData/batch2019/'

# ORIGINALLY RAN ON YEAR2019 DATA
# pathname = '../../data-se3-path-planner/yearData/year2019/'
pathname = '../../batch2019/'
sys.path.append(pathname)


    

In [2]:

    

df = pd.read_csv(pathname + 'Jan65_results.csv')
GMAT_MJD_OFFSET = 29999.5
t = df['DefaultSC.A1ModJulian'] + GMAT_MJD_OFFSET
x = df['DefaultSC.gse.X']
y = df['DefaultSC.gse.Y']
z = df['DefaultSC.gse.Z']

#
xa = np.array(x)
ya = np.array(y)
za = np.array(z)
print(za)
ta = np.array(t)

c_system = 'SM'
spacecraft = coord.Coords([[i,j,k] for i,j,k in zip(x,y,z)], 'GSE', 'car')
spacecraft.ticks = Ticktock(t,'MJD')
spacecraft = spacecraft.convert('SM','car')
points = 10000
# this figure validates what I already expected
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.plot(spacecraft.x[:points],spacecraft.y[:points],spacecraft.z[:points])
plt.title('SM Orbit')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()


    

    




[ 5626.54749147  5700.06199539  5772.07267905 ...,  1244.12738047
  1836.66693301  1984.04644116]






    

In [3]:

    

# adding this important term to generate the -30, -15, 0, +15, +30 plots
# psi = tsyg.getTilt(t)
r = np.linspace(0,10)
psi = np.deg2rad(-30)
psi = np.array(psi)
psiAlt = np.deg2rad(tsyg.getTilt(t))
np.array(psiAlt)
print(psi)
phi_0 = 0.24
alpha1 = 0.1287
alpha2 = 0.0314
phi_0 = 0.24
phi_1 = phi_0 - (alpha1*psi + alpha2*psi**2)

phi_1Alt = phi_0 - (alpha1*psiAlt + alpha2*psiAlt**2)

plt.plot(np.rad2deg(psiAlt),np.rad2deg(phi_1Alt))
plt.xlabel('psi')
plt.ylabel('phi_1')
plt.show()

plt.plot(t,np.rad2deg(phi_1Alt))
plt.xlabel('time')
plt.ylabel('phi_1Alt')
plt.show()

plt.plot(t,np.rad2deg(psiAlt))
plt.ylabel('psiAlt')
plt.xlabel('time')
plt.show()


    

    




-0.5235987755982988






    













    













    

In [4]:

    

# getting phi_c

phi_c  = np.arcsin( (np.sqrt(r))/(np.sqrt(r + (1/np.sin(phi_1))**2 -1)  )) + psi
# convert to degrees here if desired

zgsm = r*np.cos(phi_c)
xgsm = r*np.sin(phi_c)
plt.figure()
plt.plot(xgsm, zgsm)
# plt.axis('equal')
# plt.ylim([3,10])
# plt.xlim([-3,4])
plt.title('gsm coordinates for 10 Re and varying $\psi$')
plt.show()


    

In [5]:

    

#let's get all the plots below

psiVector_d = np.array([-30, -15, 0, 15, 30])
psiVector_r = np.deg2rad(psiVector_d)
plt.figure()
for psi, psi_d in zip(psiVector_r, psiVector_d):
    phi_1 = phi_0 - (alpha1*psi + alpha2*psi**2)
    phi_c  = np.arcsin( (np.sqrt(r))/(np.sqrt(r + (1/np.sin(phi_1))**2 -1)  )) + psi
    # convert to degrees here if desired

    zgsm = r*np.cos(phi_c)
    xgsm = r*np.sin(phi_c)
    plt.plot(xgsm, zgsm, label='%s'%(str(psi_d)))

# add in a circle to represent the earth
param = np.linspace(0,2*np.pi)
xEarth = np.cos(param)
yEarth = np.sin(param)
plt.plot(xEarth,yEarth,label='Earth')

# fancy stuff
plt.legend()
plt.title('$\psi$ effect on cusp location')
plt.xlabel('$x_{gsm}$')
plt.ylabel('$z_{gsm}$')

# axis limits to make it "isotropic"
plt.axis('equal')
# plt.ylim([-1,10])
# plt.xlim([-1,10])
plt.show()


    

In [ ]: