Decay of the ISS Orbit.

The ISS operates within the LEO environment and, given its size, at-mospheric drag results in noticeable orbit decay that requires periodic corrections via thrusters. Toexamine this in more detail, suppose that the ISS is at an altitude of 350 km and in a circular orbit.If the ballistic coefficient for the ISS is estimated as B= 165 kg/m^2, estimate the orbit decay over aperiod of six months.

In [1]:
#Rex Calabrese
import numpy as np
import numpy.linalg as la

import scipy as sp
from scipy import stats
from scipy.integrate import odeint
from scipy.integrate import solve_ivp
from scipy.signal import find_peaks

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import Image

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

from matplotlib import rc
rc('text', usetex=False) 

fs_default = (11,4)
In [3]:
Image("images/hw8_diagram_1.png",width=400)
Out[3]:

Video of results!

In [4]:
%%html
<video width="600" autoplay loop muted>
  <source src="Drag_animation_2.mp4" type="video/mp4" />
</video>

Physical Parameters

Universal Gravatational Constant

$ G= 6.67408 \cdot 10^{-20} \frac{km^3}{kg \cdot s^2}$

Mass Earth

$m_{e} = 5.974 \cdot 10^{24} kg$

Radius of earth

$r_{e} = 6378 \cdot 10^{24} km$

Grav Constant Earth

$\mu = 3.986 \cdot 10^{5} \frac{{kg}^3}{s^{2}}$

In [6]:
me = 5.974e24 #kg
G = 6.67259e-20 #km^3 · kg^-1 · s^-2
re = 6378.0   #km
mu = G*me # km^3 · s^-2
Out[6]:
398620.52660000004

Satellite Parameters

$B=\frac{m}{C_{D} A}$

$\mathbf{a}_{d r a g}=-\frac{1}{2} \frac{\rho}{B}\left|\mathbf{v}_{r e l}\right| \mathbf{v}_{r e l}$

$\mathbf{v}_{r e l}=\mathbf{v}-\omega \times \mathbf{r}$

In [7]:
A = 1. * (1. / 1000)**2  #km
Cd = 2.  #unitless
m = 100 #kg

#Ballistic Coefficent
B = m / (Cd * A) #kg / km^2

print('Ballistic Coefficent: B = {0:0.1e} kg / km^2'.format(B))

Ballistic Coefficent: B = 5.0e+07 kg / km^2

Initial conditions

$\mathbf{r}_{0} = [5873.40, -658.52, 3007.49] km $

$\mathbf{v}_{0} = [-2.90, 4.09, 6.14] \ \ \frac{km}{s} $

In [8]:
r0 = [5873.40, -658.52, 3007.49] #km 
v0 = [-2.90, 4.09, 6.14] #km/s

u0 = np.hstack((r0, v0))
In [22]:
import orbital
orbital_elements_r0 = orbital.elements_from_state_vector(np.array(r0), np.array(v0), mu)
print('ASIDE:\nuse library to check my functions for calculating orbital elements (given position & velocity)\n')
print(orbital_elements_r0)

ASIDE:
use library to check my functions for calculating orbital elements (given position & velocity)

OrbitalElements(a=6946.493413751562, e=0.05129398476251154, i=1.1364997307253533, raan=5.933256915059585, arg_pe=1.033814758799821, f=5.772892295021246)

Density Modelling

$\text{for z < 100 km} : \rho = 1 \cdot 10^9 \frac{kg}{{km}^3} $

$\text{for z >= 100 km} : \rho = 2 \cdot 10^9 \cdot z ^ {-8.284} \frac{kg}{{km}^3}$

In [23]:
def rho_fn(z):
    if z < 100:
        return False
    if z >= 100: 
        return 2e9 * z**(-8.284)
    
fig00, ax00 = plt.subplots(figsize=fs_default)
z_plot = np.linspace(5000,7000,500)
rho_plot = np.zeros(np.size(z_plot))

for i in range(np.size(z_plot)):
    rho_plot[i] = rho_fn(z_plot[i])
    
ax00.plot(z_plot,rho_plot,c='grey')

plt.xlabel('Altitude (km)')
plt.ylabel(r'$Density (\frac{kg}{km^{3}})$')

plt.show()

Relative Velocity

In [24]:
omega_earth = [0., 0., 7.27e-5] #rad/s  

def v_rel_fn(r,v,omega_earth):
    return v - np.cross(omega_earth,r)

def a_drag_fn(r,r_mag,v):
    v_rel = v - np.cross(omega_earth,r)
    v_rel_mag = la.norm(v_rel)
    
    rho = rho_fn(r_mag - re)
    return   -0.5 * (rho / B) * v_rel * v_rel_mag
In [25]:
# %% State Vector Derivative Function

def deriv(u, dt):
    r = u[0:3]
    r_mag = la.norm(r)
    v = u[3:6]
    n = -mu / r_mag**3
    a_d = a_drag_fn(r,r_mag,v) 
    return [u[3],     #  dotu[0] = u[3]'
            u[4],     #  dotu[1] = u[4]'
            u[5],     #  dotu[2] = u[5]'
            u[0] * n + a_d[0],       #  dotu[3] = u[0] * n
            u[1] * n + a_d[1],       #  dotu[4] = u[1] * n
            u[2] * n + a_d[2]]       #  dotu[5] = u[2] * n

Two Body Model

In [26]:
def twobody(dt,u0):
    u = odeint(deriv, u0, dt, atol=1.0E-8) #1E-8
    r = u[:,0:3]
    v = u[:,3:6]
    return r,v
In [27]:
Y2S = 365 * 24 * 3600
D2S = 24*3600 
Y2D = 365 

#ORIGINAL: 5 years
nDays = 5 * Y2D
#SMALLER
#nDays = 1 * Y2D


ElapsedTime = nDays*D2S 


TimePeriod = int(ElapsedTime / 5)

CALCULATE_HIGH_RESOLUTION = True
if CALCULATE_HIGH_RESOLUTION == True:
    
    dt_hires = np.arange(0.0, ElapsedTime, 1)  
    r1_hires = np.zeros((dt_hires.size,3))
    v1_hires = np.zeros((dt_hires.size,3))
    
    u0_i = u0
    for i in range(5): #5
        
        i_L = i * TimePeriod
        i_R = (i+1) * TimePeriod
        
        dt_i = dt_hires[i_L:i_R]
        
        r1_hires[i_L:i_R],v1_hires[i_L:i_R] = twobody(dt_i,u0_i)
        u0_i = np.hstack((r1_hires[i_R-1,:],v1_hires[i_R-1,:]))
   
    
LOAD_HIGH_RESOLUTION = False        #127 million time steps
if LOAD_HIGH_RESOLUTION == True:
    r1_hires = np.load('./results_export/full/2/r1.npy')
    v1_hires = np.load('./results_export/full/2/v1.npy')
    dt1_hires = np.load('./results_export/full/2/dt.npy')

#fix this for solving vs loading data
dt1_hires = dt_hires
del dt_hires

Down-Sampling Data

In [28]:
# %% Chop off the inside of the high resolution trajectory data
ends = 10000
r1_ends = np.vstack((r1_hires[0:ends,:],r1_hires[-ends:,:]))
v1_ends = np.vstack((v1_hires[0:ends,:],v1_hires[-ends:,:]))
dt1_ends = np.vstack((dt1_hires[0:ends],dt1_hires[-ends:]))
NU = np.zeros(dt1_ends.size)
In [29]:
# %% Keep the tail end of hires trajectory data.
r1 = r1_hires[-100000:]  #100k
v1 = v1_hires[-100000:]  #100k
dt1 = dt1_hires[-100000:]  #100k
dt1_days = dt1 / D2S

r1_mag = la.norm(r1,axis=1)
v1_mag = la.norm(v1,axis=1)
In [30]:
# %% Downsample hires trajectory data.
sample_interval = int(dt1_hires.size / 100) #int(D2S / 10)
sample_interval = 97
r2 = r1_hires[0::sample_interval,:].copy()
del r1_hires

v2 = v1_hires[0::sample_interval,:].copy()
del v1_hires

dt2_days = dt1_hires[0::sample_interval].copy() / D2S
del dt1_hires

#Will be plotting with time in years.
dt2 = dt2_days / Y2D
In [ ]:
# Subscript 2 data is downsampled.
# Subscript 1 data is clipped.

Calculate Orbital Elements

In [31]:
h2 = np.cross(r2,v2,axis=1)
r2_mag = np.linalg.norm(r2, axis=1)
v2_mag = np.linalg.norm(v2, axis=1)

eccvec2 = np.cross(v2, h2)/mu   -    np.divide(r2,r2_mag[:,None])
eccmag2 = la.norm(eccvec2,axis=1)


energy2 = (v2_mag**2 / 2) - (mu/r2_mag)
a2 = -mu/(2*energy2)

r2_a = a2 * (1 + eccmag2)
r2_p = a2 * (1 - eccmag2)

Plot Semi-Major axis over time

In [32]:
fig0, ax0 = plt.subplots(figsize=fs_default)
ax0.plot(dt2,a2,c='grey')

slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,a2)
line = slope*dt2+intercept
ax0.plot(dt2,line,'r--',lw=1)

ax0.set_title('Semi-Major axis :  $\Delta a$/$\Delta t$={0:0.2f} km/yr'.format(slope,intercept))  
ax0.set_xlabel('Time (yr)')
ax0.set_ylabel('Semi-Major Axis (km)')
ax0.legend(['e','linear regression'],loc='upper right')
ax0.grid(True)
plt.show()
print('semi major axis decay over {0:0.2f} years'.format(dt2[-1] ))
print('Δa = {0:0.3f} km'.format(a2[-1] - a2[0] ))

semi major axis decay over 5.00 years
Δa = -8.606 km

Plot eccentricity over time

In [33]:
fig1, ax1 = plt.subplots(figsize=fs_default)
ax1.plot(dt2,eccmag2,c='grey')

slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,eccmag2)
line = slope*dt2+intercept
ax1.plot(dt2,line,'r--',lw=1)

ax1.set_title('Eccentricity :  $\Delta e$/$\Delta t$={0:0.2e} yr$^-1$'.format(slope,intercept))  
ax1.set_xlabel('Time (yr)')
ax1.set_ylabel('Eccentricity (km/km)')
ax1.legend(['e','linear regression'],loc='upper right')
plt.show()
print('Eccentricity decay over {0:0.2f} years'.format(dt2[-1] ))
print('Δe = {0:0.4f}'.format(eccmag2[-1] - eccmag2[0] ))

Eccentricity decay over 5.00 years
Δe = -0.0013

Plot Perigee and Apoapsis Values over Time

In [37]:
fig3, ax3 = plt.subplots(figsize=fs_default)

c_rp = 'orange'
c_ra = 'purple'

#apoapsis
ax3.plot(dt2,r2_a - re,c=c_ra,lw=3,label='$apogee$')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,r2_a - re)
line = slope*dt2+intercept
ax3.plot(dt2,line,'r--',lw=1,label='$\Delta r_a$/$\Delta t$ = {0:0.1f} km/yr'.format(slope))

#perigee
ax3.plot(dt2,r2_p  - re,c=c_rp,lw=3,label='$perigee$')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,r2_p - re)
line = slope*dt2+intercept
ax3.plot(dt2,line,'b--',lw=1,label='$\Delta r_p$/$\Delta t$ = {0:0.1f} km/yr'.format(slope))


ax3.set_title('Atmospheric Drag: Perigee and Apoapsis Values')
ax3.set_ylabel('$Altitude$   (km)')
ax3.set_xlabel('Time  (years)')
ax3.legend(loc='center right')
ax3.grid(True)
plt.show()

Re-frame the Above Plot.

In [39]:
# %% Plot flat orbit with r_a and r_p

r2_a_shift = r2_a-r2_a[0]
r2_p_shift = r2_p-r2_p[0]

fig4, ax4 = plt.subplots(figsize=fs_default)

#apoapsis
ax4.plot(dt2,r2_a_shift,c=c_ra,lw=3,label='$r_a$')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,r2_a_shift)
line = slope*dt2+intercept
ax4.plot(dt2,line,'r--',lw=1,label='$\Delta r_a$/$\Delta t$ = {0:0.1f} km/yr'.format(slope))

#perigee
ax4.plot(dt2,r2_p_shift,c=c_rp,lw=3,label='$r_p$')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,r2_p_shift)
line = slope*dt2+intercept
ax4.plot(dt2,line,'b--',lw=1,label='$\Delta r_p$/$\Delta t$ = {0:0.1f} km/yr'.format(slope))


ax4.set_title('Atmospheric Drag: Perigee and Apoapsis Values (arbitrary initial values)')
ax4.set_ylabel('$r$   (km)')
ax4.set_xlabel('Time  (years)')
ax4.legend(loc='lower left')
plt.show()

The above plot shows a relatively steady perigee. We know the drag force is largest at perigee due to increased atmospheric density, similar to an instantaneous ΔV maneuver at perigee, this location remains 'pinned'.

'Manually' find r_a, r_p index locations.

In [40]:
# %% 'Micro View': FIND r_a and r_p locations, and plot over r vs t.

r1_apo = find_peaks(r1_mag)
r1_peri = find_peaks(-r1_mag)

plt.figure(figsize=fs_default)

#Plot r
plt.plot(dt1_days,r1_mag,color='k')

#Plot the found r_a,t_p locations
for i in range(r1_peri[0].size - 1): #
    j = r1_peri[0][i]
    j1 = r1_apo[0][i]
    plt.scatter(dt1_days[j],r1_mag[j], c=c_ra)
    plt.scatter(dt1_days[j1],r1_mag[j1], c=c_rp)
plt.legend(['|r|','Apogee','Perigee'],loc=1)
plt.title('Test: Find Perigee and Apoasis Values Through Peak Searching')
plt.show()

Recalculate Deceleration and Velocity Values (were not saved from inside 2body loop)

In [42]:
# RECALCULATE ACCEL (#Using r1, v1 ,dt1 )
# For the tail end of the high resolution data!!!
a1 = np.zeros((dt1.size,3))
for i in range(dt1.size):
    a1[i,:] = a_drag_fn(r1[i,:],r1_mag[i],v1[i,:])
    
a1_mag =  la.norm(a1,axis=1)
In [43]:
# Recalculate V_rel
v1rel = np.zeros((dt1.size,3))
for i in range(dt1.size):
    v1rel[i,:] = v_rel_fn(r1[i,:],v1[i,:],omega_earth)
    
v1rel_mag = la.norm(v1rel,axis=1)

Trying to understand where the negative deceleration is coming from.

In [44]:
dt1_yr = dt1 / Y2S
dt1_day = dt1 / D2S

fig5, ax5 = plt.subplots(3, figsize=(11,8), sharex=True)
ax5[0].plot(dt1_day, r1[:,0],c='r',label='$r_x$')
ax5[0].plot(dt1_day, r1[:,1],c='g', label='$r_y$')
ax5[0].plot(dt1_day, r1[:,2],c='b', label='$r_z$')
ax5[0].plot(dt1_day, r1_mag[:],c='k', label='|r|')

ax5[1].plot(dt1_day, v1rel[:,0],c='r', label='$v_x$')
ax5[1].plot(dt1_day, v1rel[:,1],c='g', label='$v_y$')
ax5[1].plot(dt1_day, v1rel[:,2],c='b', label='$v_z$')
ax5[1].plot(dt1_day, v1rel_mag,c='k', label='|v|')

ax5[2].plot(dt1_day, a1[:,0],c='r', label='$a_x$')
ax5[2].plot(dt1_day, a1[:,1],c='g', label='$a_y$')
ax5[2].plot(dt1_day, a1[:,2],c='b', label='$a_z$')
ax5[2].plot(dt1_day, a1_mag,c='k', label='|a|')

#r_a and r_p
for i in range(r1_apo[0].size): #
    j1 = r1_apo[0][i]
    if i ==0:
        ax5[0].scatter(dt1_day[j1],r1_mag[j1], c=c_rp, label='A')
        ax5[1].scatter(dt1_day[j1],v1rel_mag[j1], c=c_rp,label='A')    
        ax5[2].scatter(dt1_day[j1],a1_mag[j1], c=c_rp,label='A')
    else:
        ax5[0].scatter(dt1_day[j1],r1_mag[j1], c=c_rp)
        ax5[1].scatter(dt1_day[j1],v1rel_mag[j1], c=c_rp)    
        ax5[2].scatter(dt1_day[j1],a1_mag[j1], c=c_rp)
    
    
for i in range(r1_peri[0].size): #
    j = r1_peri[0][i]
    if i ==0:
        ax5[0].scatter(dt1_day[j],r1_mag[j], c=c_ra, label='P')
        ax5[1].scatter(dt1_day[j],v1rel_mag[j], c=c_ra, label='P')   
        ax5[2].scatter(dt1_day[j],a1_mag[j], c=c_ra, label='P')
    else:
        ax5[0].scatter(dt1_day[j],r1_mag[j], c=c_ra)
        ax5[1].scatter(dt1_day[j],v1rel_mag[j], c=c_ra)   
        ax5[2].scatter(dt1_day[j],a1_mag[j], c=c_ra)



for i in range(3): ax5[i].legend(loc=1)
ax5[0].set_title('Position, Relative Velocity, Deceleration (final orbit)')

for i in range(3): ax5[i].set_xlim((5*Y2D - 0.18,5*Y2D))
ax5[2].set_xlabel('Time (days)')

ax5[0].set_ylabel('Position (km)')
ax5[1].set_ylabel('Relative Velocity (km/s)')
ax5[2].set_ylabel('Deceleration (km/s$^2$)')
plt.show()

The above plot shows the deceleration values are inversely proportional to altitude.

Negative deceleration may be due to the orbit and resulting relative velocities w.r.t earth rotation.

In [46]:
%matplotlib notebook
#fig = plt.figure(dpi=100) #
fig = plt.figure(figsize=(6,6),dpi=135)
#fig = plt.figure(dpi=100)
#ax = fig.add_subplot(111, projection = '3d')
ax = Axes3D(fig)
ax.set_xlim([-8000,8000])
ax.set_ylim([-8000,8000])
ax.set_zlim([-8000,8000])

###########################################
#Add Earth (Bottom Half)
ue1, ve1 = np.mgrid[0:2*np.pi:60j, 0:np.pi:60j]
xe=re*np.cos(ue1)*np.sin(ve1)
ye=re*np.sin(ue1)*np.sin(ve1)
ze=re*np.cos(ve1)
clip0 = int(xe.size / 10)
ax.plot_surface(xe, ye, ze, lw=0.1,alpha=0.2,color='deepskyblue',zorder = -10)
ax.plot_wireframe(xe, ye, ze, lw=0.1,alpha=0.99,color='deepskyblue',zorder = -10)

#Show  Earth Axis of rotation.
ax.plot([0,0],[0,0],[-re*1.2,re*1.2],color='k',lw=4,ls='-',alpha=0.2,label='Axis of Rotation: Earth')
ax.scatter(0,0,0,color='k',alpha=0.8)

###########################################
#Plot First Full Orbit
clip1 = 0
clip2 = 7500
ax.plot(r1[clip1:clip2,0], r1[clip1:clip2,1], r1[clip1:clip2,2],'k',lw=1,zorder = 2.5,label='Initial Orbit')

###########################################
#apoapsis
c_v = 'mediumspringgreen'
c_a = 'hotpink'

i=0
j1 = r1_apo[0][i]
ax.scatter(r1[j1,0],r1[j1,1],r1[j1,2], c=c_ra,s=30) #show location of apoapsis
ax.quiver(0,0,0,r1[j1,0],r1[j1,1],r1[j1,2], color=c_ra,label='$r_{Apoapsis}$') # show vector to apoapsis

#Show velocity vector at apoapsis
scale1 = 1000
ax.quiver(r1[j1,0],r1[j1,1],r1[j1,2],v1rel[j1,0]*scale1,v1rel[j1,1]*scale1,v1rel[j1,2]*scale1, color=c_v)

#Show decel vector at apoapsis
scale2 = 1e15
ax.quiver(r1[j1,0],r1[j1,1],r1[j1,2],a1[j1,0]*scale2,a1[j1,1]*scale2,a1[j1,2]*scale2, color=c_a)

###########################################
#perigee
i=0
j1 = r1_peri[0][i]
ax.scatter(r1[j1,0],r1[j1,1],r1[j1,2], c=c_rp,s=30)#show location of perigee
ax.quiver(0,0,0,r1[j1,0],r1[j1,1],r1[j1,2], color=c_rp,label='$r_{Periapsis}$') # show vector to perigee

#Show velocity vector at perigee
ax.quiver(r1[j1,0],r1[j1,1],r1[j1,2],v1rel[j1,0]*scale1,v1rel[j1,1]*scale1,v1rel[j1,2]*scale1, color=c_v,label='$\\vec v$')

#Show deceleration vector at perigee (Scales of deceleration are not equal)
scale3 = 1e15
ax.quiver(r1[j1,0],r1[j1,1],r1[j1,2],a1[j1,0]*scale3,a1[j1,1]*scale3,a1[j1,2]*scale3, color=c_a,label='$ \\vec a$ (not to scale)')

###########################################
#Label
ax.set_xlabel('x (km)')
ax.set_ylabel('y (km)')
ax.set_zlabel('z (km)')

#Set ticks
ax.set_xticks([-5000,0,5000])
ax.set_yticks([-5000,0,5000])
ax.set_zticks([-5000,0,5000])


plt.legend()
ax.set_axis_off()
  
plt.show()    

animate_rotate = False
if animate_rotate==True:
    for angle in range(0, 360):
        ax.view_init(0, angle)
        fig.canvas.draw()
        plt.pause(.001)

ax.view_init(-22,70)

The above figure is another representation of the inverse proportionality between deceleration and altitude.

DO AGAIN BUT ANIMATE THE TRAJECTORY

In [47]:
%matplotlib notebook
#fig = plt.figure(dpi=100) #
#fig = plt.figure(dpi=100)


#Very good
fig = plt.figure(figsize=(6,6),dpi=135)



#ax = fig.add_subplot(111, projection = '3d')
ax = Axes3D(fig)
axlim = 7000
ax.set_xlim([-axlim,axlim])
ax.set_ylim([-axlim,axlim])
ax.set_zlim([-axlim,axlim])
ax.view_init(-22,70)

###########################################
#Add Earth (Bottom Half)
ue1, ve1 = np.mgrid[0:2*np.pi:60j, 0:np.pi:60j]
xe=re*np.cos(ue1)*np.sin(ve1)
ye=re*np.sin(ue1)*np.sin(ve1)
ze=re*np.cos(ve1)
clip0 = int(xe.size / 10)
ax.plot_surface(xe, ye, ze, lw=0.1,alpha=0.2,color='deepskyblue',zorder = -10)
ax.plot_wireframe(xe, ye, ze, lw=0.1,alpha=0.99,color='deepskyblue',zorder = -10)

#Show  Earth Axis of rotation.
ax.plot([0,0],[0,0],[-re*1.2,re*1.2],color='k',lw=4,ls='-',alpha=0.2,label='Axis of Rotation')
ax.scatter(0,0,0,color='k',alpha=0.8)

###########################################
#Plot First Full Orbit
#clip1 = 0
#clip2 = 7500
#

###########################################
#apoapsis
c_v = 'mediumspringgreen'
c_a = 'hotpink'


#Label
ax.set_xlabel('x (km)')
ax.set_ylabel('y (km)')
ax.set_zlabel('z (km)')

#Set ticks
ax.set_xticks([-5000,0,5000])
ax.set_yticks([-5000,0,5000])
ax.set_zticks([-5000,0,5000])


#ax.set_axis_off()
def clear_lines():
    quiv_accel.remove()
    quiv_vel.remove()
    quiv_radius.remove()
    

plt.show()

LW_vec = 3


#For parameter output in plot
line = -20000
ypos_0 = 0.7 #-0.1
xpos_0 =0.02
    
    
animation_range = 1000
for j in range(animation_range):    

    
    i = j * 30   # * 80 for display
    
    #ORBIT
    ax.plot((r1[i,0],r1[i+1,0]), (r1[i,1],r1[i+1,1]), (r1[i,2],r1[i+1,2]),'k',lw=1,color='k',zorder = 2.5,label='Orbit')
    
    #RADIUS VECTOR
    quiv_radius = ax.quiver(0,0,0,r1[i,0],r1[i,1],r1[i,2], color=c_rp,linewidths=LW_vec, label='Vector to sat') 
    
    #VELOCITY VECTOR
    s1 = 1000
    quiv_vel = ax.quiver(r1[i,0],r1[i,1],r1[i,2],v1rel[i,0]*s1,v1rel[i,1]*s1,v1rel[i,2]*s1,color=c_v,linewidths=LW_vec,label='Velocity')
    
    #DECEL VECTOR
    s2 = 1e14
    quiv_accel = ax.quiver(r1[i,0],r1[i,1],r1[i,2],a1[i,0]*s2,a1[i,1]*s2,a1[i,2]*s2, color=c_a, linewidths=LW_vec,label = 'Drag')
    
    
    
    text_below = '$altitude =$'+'{0:.0f} km \ndrag = {1:.1e} $km/s^2$'.format(r1_mag[i] - re,a1_mag[i]  )      
    text_parameter_1 = ax.text(xpos_0, ypos_0 + 0*line ,0,text_below,
        horizontalalignment='left',
        verticalalignment='bottom',
        size='large',
        bbox=dict(facecolor='gray', alpha=0.3),
        transform = ax.transAxes)
    
    

    
    #SHOW
    fig.canvas.draw()
    plt.pause(0.0001)

    if j == 0:
        plt.legend(loc='lower right',  borderaxespad=0.) #bbox_to_anchor=(1.05, 1),
    
    #Save frames for video
    export_animation = False
    if export_animation == True:
        frame_name = 'animation/frame_{}'.format(j)
        plt.savefig(frame_name)
    
    #REMOVE PREVIOUS
    if j < animation_range-1:
        quiv_accel.remove()
        quiv_vel.remove()
        quiv_radius.remove()
        text_parameter_1.remove()

        del quiv_accel
        del quiv_vel
        del quiv_radius
        del text_parameter_1

    #quiv_vel.set_alpha(0)
    #quiv_radius.set_alpha(0)
    #quiv_accel.set_alpha(0)
###########################################
  
plt.show()
In [ ]: