Interesting Plots

Air Density as function of the Solar Cycle

import satkit as sk
import plotly.graph_objects as go
import numpy as np
import math as m

start = sk.time(1995, 1, 1)
stop = sk.time(2022, 12, 31)
duration = stop - start
timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)])

altitude = 400e3
rho_400, _temperature = zip(*np.array([sk.density.nrlmsise(altitude, 0, 0, x) for x in timearray]))

altitude = 500e3
rho_500, emperature = zip(*np.array([sk.density.nrlmsise(altitude, 0, 0, x) for x in timearray]))


fig = go.Figure()
fig.add_trace(go.Scatter(x=[t.datetime() for t in timearray], y=rho_400,
                          mode='lines', name='Altitude = 400km', line=dict(color='black', width=1)))
fig.add_trace(go.Scatter(x=[t.datetime() for t in timearray], y=rho_500,
                         mode='lines', name='Altitude = 500km', line=dict(color='black', width=1, dash='dash')))
fig.update_layout(legend=dict(
    yanchor="top",
    y=0.99,
    xanchor="left",
    x=0.5
))
fig.update_yaxes(type="log", title_font=dict(size=8), tickfont=dict(size=8), title='Density [kg/m<sup>3</sup>]')
fig.update_xaxes(title='Year')
fig.update_layout(yaxis_tickformat=".2e", width=650, height=512,
                  title='Air Density Changes with Solar Cycle',)
fig.update_xaxes(showline=True, linewidth=2, linecolor='black', mirror=True)
fig.update_yaxes(showline=True, linewidth=2, linecolor='black', mirror=True)
fig.update_layout(xaxis=dict(
        gridcolor='#dddddd',
        gridwidth=1,
    ),
    yaxis=dict(
        gridcolor='#dddddd',
        gridwidth=1,
    ),)
fig.show()

Forces acting on satellite as a function of altitude

Note: this is insprired by a very similar plot in Montenbruck & Gill

## Force components acting on satellite vs altitude
import numpy as np
import satkit as sk
import math as m
import plotly.graph_objects as go


N = 1000
range_arr = np.logspace(m.log10(6378.2e3 + 100e3), m.log10(50e6), N)
grav1v = np.array([sk.gravity(np.array([a, 0, 0]), order=1) for a in range_arr])
grav1 = np.linalg.norm(grav1v, axis=1)
grav2v = np.array([sk.gravity(np.array([a, 0, 0]), order=2) for a in range_arr])
grav2 = np.linalg.norm(grav2v-grav1v, axis=1)

grav6v = np.array([sk.gravity(np.array([a, 0, 0]), order=6) for a in range_arr])
grav5v = np.array([sk.gravity(np.array([a, 0, 0]), order=5) for a in range_arr])
grav6 = np.linalg.norm(grav6v-grav5v, axis=1)

aoverm = 0.01
Cd = 2.2
Cr = 1.0

# Drag force at maximum & minimum density
didx = np.argwhere(range_arr - sk.consts.earth_radius < 800e3).flatten()
tm_max = sk.time(2001, 12, 1)
tm_min = sk.time(1996, 12, 1)
rho_max = np.array([sk.density.nrlmsise(a-sk.consts.earth_radius, 0, 0, tm_max)[0] for a in range_arr[didx]])
rho_min = np.array([sk.density.nrlmsise(a-sk.consts.earth_radius, 0, 0, tm_min)[0] for a in range_arr[didx]])
varr = np.sqrt(sk.consts.mu_earth/(range_arr + sk.consts.earth_radius))
drag_max = 0.5 * rho_max * varr[didx]**2 * Cd * aoverm
drag_min = 0.5 * rho_min * varr[didx]**2 * Cd * aoverm

moon_range = np.linalg.norm(sk.jplephem.geocentric_pos(sk.solarsystem.Moon, sk.time(2023, 1, 1)))
moon = sk.consts.mu_moon*((moon_range-range_arr)**(-2) - moon_range**(-2))
sun = sk.consts.mu_sun*((sk.consts.au-range_arr)**(-2) - sk.consts.au**(-2))

a_radiation = 4.56e-6 * 0.5 * Cr * aoverm * np.ones(range_arr.shape)

def add_line(fig, x, y, text, frac=0.5, ax=-20, ay=-30):
    fig.add_scatter(x=x, y=y, mode='lines', name=text, line=dict(width=2, color='black'))
    idx = int(len(x)*frac)
    fig.add_annotation(
        x=np.log10(x[idx]),
        y=np.log10(y[idx]),
        text=text,
        showarrow=True,
        font=dict(
            size=12,
            color='black',
        ),
        align="center",
        ax=ax,
        ay=ay,
        arrowcolor="#636363",
        arrowsize=1,
        arrowwidth=2,
        arrowhead=6,
        borderwidth=2,
        borderpad=4,
        opacity=0.8
    )

fig = go.Figure()
add_line(fig, range_arr/1e3, grav1/1e3, 'Gravity')
add_line(fig, range_arr/1e3, grav2/1e3, 'J2', 0.2, 0, -10)
add_line(fig, range_arr/1e3, grav6/1e3, 'J6', 0.8, 0, -10)
add_line(fig, range_arr[didx]/1e3, drag_max/1e3, 'Drag Max', 0.7, 30, 0)
add_line(fig, range_arr[didx]/1e3, drag_min/1e3, 'Drag Min', 0.8, 10, 30)
add_line(fig, range_arr/1e3, moon/1e3, 'Moon', 0.8, -10, -10)
add_line(fig, range_arr/1e3, sun/1e3, 'Sun', 0.7, -10, 10)
add_line(fig, range_arr/1e3, a_radiation/1e3, 'Radiation\nPressure', 0.3, -10, 10)

fig.update_xaxes(type="log", title='Distance from Earth Origin [km]', range=np.log10([6378.1, 50e3]))
fig.update_yaxes(type="log", title='Acceleration   [km/s<sup>2</sup>]', tickformat=".1e")
fig.update_layout(title='Satellite Forces vs Altitude', width=650, height=650)
fig.update_xaxes(showline=True, linewidth=2, linecolor='black', mirror=True)
fig.update_yaxes(showline=True, linewidth=2, linecolor='black', mirror=True)
fig.update_layout(showlegend=False, xaxis=dict(
        gridcolor='#dddddd',
        gridwidth=1,
    ),
    yaxis=dict(
        gridcolor='#dddddd',
        gridwidth=1,
    ),)
fig.show()

Difference in angle between IAU-2010 reduction and approximate calculation

import plotly.graph_objects as go
import numpy as np
import satkit as sk
import math as m

start = sk.time(2023, 1, 1)
stop = sk.time(1970, 1, 1)
duration = stop-start
timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)])
qexact = sk.frametransform.qgcrf2itrf(timearray)
qapprox = sk.frametransform.qgcrf2itrf_approx(timearray)
qdiff = np.array([q1 * q2.conj for q1, q2 in zip(qexact, qapprox)])
theta = np.array([q.angle for q in qdiff])
fig = go.Figure()
fig.add_trace(go.Scatter(x=[t.datetime() for t in timearray], y=theta*180.0/m.pi*3600,
                         mode='lines', name='Angle Difference',
                         line=dict(width=1, color='black')))
fig.update_layout(title='Angle Error of Approximate ITRF to GCRF Rotation', width=650, height=512,
                  font=dict(size=14))
fig.update_yaxes(title='Error [arcsec]')
fig.update_xaxes(title='Year')

Polar Motion

import plotly.graph_objects as go
import numpy as np
import satkit as sk
import math as m

start = sk.time(2023, 1, 1)
stop = sk.time(1970, 1, 1)
duration = stop-start
timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)])

dut1, xp, yp, lod, dX, dY = zip(*[sk.frametransform.earth_orientation_params(t) for t in timearray])
fig = go.Figure()
fig.add_trace(go.Scatter3d(x=xp, y=yp, z=[t.datetime() for t in timearray],  mode='lines', name='Polar Motion',))
fig.update_scenes(xaxis_title='X [arcsec]', yaxis_title='Y [arcsec]', zaxis_title='Year')
fig.update_layout(title='Polar Motion', width=650, height=600, font=dict(size=12))
fig.show()

Precession / Nutation

import satkit as sk
import numpy as np
import math as m
import plotly.graph_objects as go

# Precession / Nutation is the difference between the IAU-2006 GCRS frame and the CIRS frame
start = sk.time(2000, 1, 1)
stop = sk.time(2020, 1, 1)
duration = stop - start
N = 1000
timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, N)])
qarray = np.array([sk.frametransform.qcirs2gcrf(t) for t in timearray])
theta = np.fromiter((q.angle for q in qarray), float)
rots = np.array([q * np.array([0.0, 0.0, 1.0]) for q in qarray])
pline = np.linspace(0,1,N)
pv = np.polyfit(pline, rots[:,0], 1)
rx0 =rots[:,0] - np.polyval(pv, pline)

rad2asec = 180.0/m.pi * 3600
fig = go.Figure()
fig.add_trace(go.Scatter3d(x=rx0*rad2asec, y=rots[:,1]*rad2asec, z=[t.datetime() for t in timearray], mode='lines',
                          name='Precession/Nutation', line=dict(width=2, color='black')))
fig.update_layout(title='Precession / Nutation', width=650, height=600, font=dict(size=12))

fig.show()