{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Interesting Plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Air Density as function of the Solar Cycle"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import plotly.graph_objects as go\n",
"import numpy as np\n",
"import math as m\n",
"\n",
"start = sk.time(1995, 1, 1)\n",
"end = sk.time(2022, 12, 31)\n",
"duration = end - start\n",
"timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)])\n",
"\n",
"altitude = 400e3\n",
"rho_400, _temperature = zip(*np.array([sk.density.nrlmsise(altitude, 0, 0, x) for x in timearray]))\n",
"\n",
"altitude = 500e3\n",
"rho_500, emperature = zip(*np.array([sk.density.nrlmsise(altitude, 0, 0, x) for x in timearray]))\n",
"\n",
"\n",
"fig = go.Figure()\n",
"fig.add_trace(go.Scatter(x=[t.datetime() for t in timearray], y=rho_400,\n",
" mode='lines', name='Altitude = 400km', line=dict(color='black', width=1)))\n",
"fig.add_trace(go.Scatter(x=[t.datetime() for t in timearray], y=rho_500,\n",
" mode='lines', name='Altitude = 500km', line=dict(color='black', width=1, dash='dash')))\n",
"fig.update_layout(legend=dict(\n",
" yanchor=\"top\",\n",
" y=0.99,\n",
" xanchor=\"left\",\n",
" x=0.5\n",
"))\n",
"fig.update_yaxes(type=\"log\", title_font=dict(size=8), tickfont=dict(size=8), title='Density [kg/m3]')\n",
"fig.update_xaxes(title='Year')\n",
"fig.update_layout(yaxis_tickformat=\".2e\", width=650, height=512,\n",
" title='Air Density Changes with Solar Cycle',)\n",
"fig.update_xaxes(showline=True, linewidth=2, linecolor='black', mirror=True)\n",
"fig.update_yaxes(showline=True, linewidth=2, linecolor='black', mirror=True)\n",
"fig.update_layout(xaxis=dict(\n",
" gridcolor='#dddddd',\n",
" gridwidth=1,\n",
" ),\n",
" yaxis=dict(\n",
" gridcolor='#dddddd',\n",
" gridwidth=1,\n",
" ),)\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Forces acting on satellite as a function of altitude\n",
"\n",
"Note: this is insprired by a very similar plot in Montenbruck & Gill"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"## Force components acting on satellite vs altitude\n",
"import numpy as np\n",
"import satkit as sk\n",
"import math as m\n",
"import plotly.graph_objects as go\n",
"\n",
"\n",
"N = 1000\n",
"range_arr = np.logspace(m.log10(6378.2e3 + 100e3), m.log10(50e6), N)\n",
"grav1v = np.array([sk.gravity(np.array([a, 0, 0]), order=1) for a in range_arr])\n",
"grav1 = np.linalg.norm(grav1v, axis=1)\n",
"grav2v = np.array([sk.gravity(np.array([a, 0, 0]), order=2) for a in range_arr])\n",
"grav2 = np.linalg.norm(grav2v-grav1v, axis=1)\n",
"\n",
"grav6v = np.array([sk.gravity(np.array([a, 0, 0]), order=6) for a in range_arr])\n",
"grav5v = np.array([sk.gravity(np.array([a, 0, 0]), order=5) for a in range_arr])\n",
"grav6 = np.linalg.norm(grav6v-grav5v, axis=1)\n",
"\n",
"aoverm = 0.01\n",
"Cd = 2.2\n",
"Cr = 1.0\n",
"\n",
"# Drag force at maximum & minimum density\n",
"didx = np.argwhere(range_arr - sk.consts.earth_radius < 800e3).flatten()\n",
"tm_max = sk.time(2001, 12, 1)\n",
"tm_min = sk.time(1996, 12, 1)\n",
"rho_max = np.array([sk.density.nrlmsise(a-sk.consts.earth_radius, 0, 0, tm_max)[0] for a in range_arr[didx]])\n",
"rho_min = np.array([sk.density.nrlmsise(a-sk.consts.earth_radius, 0, 0, tm_min)[0] for a in range_arr[didx]])\n",
"varr = np.sqrt(sk.consts.mu_earth/(range_arr + sk.consts.earth_radius))\n",
"drag_max = 0.5 * rho_max * varr[didx]**2 * Cd * aoverm\n",
"drag_min = 0.5 * rho_min * varr[didx]**2 * Cd * aoverm\n",
"\n",
"moon_range = np.linalg.norm(sk.jplephem.geocentric_pos(sk.solarsystem.Moon, sk.time(2023, 1, 1)))\n",
"moon = sk.consts.mu_moon*((moon_range-range_arr)**(-2) - moon_range**(-2))\n",
"sun = sk.consts.mu_sun*((sk.consts.au-range_arr)**(-2) - sk.consts.au**(-2))\n",
"\n",
"a_radiation = 4.56e-6 * 0.5 * Cr * aoverm * np.ones(range_arr.shape)\n",
"\n",
"def add_line(fig, x, y, text, frac=0.5, ax=-20, ay=-30):\n",
" fig.add_scatter(x=x, y=y, mode='lines', name=text, line=dict(width=2, color='black'))\n",
" idx = int(len(x)*frac)\n",
" fig.add_annotation(\n",
" x=np.log10(x[idx]),\n",
" y=np.log10(y[idx]),\n",
" text=text,\n",
" showarrow=True,\n",
" font=dict(\n",
" size=12,\n",
" color='black',\n",
" ),\n",
" align=\"center\",\n",
" ax=ax,\n",
" ay=ay,\n",
" arrowcolor=\"#636363\",\n",
" arrowsize=1,\n",
" arrowwidth=2,\n",
" arrowhead=6,\n",
" borderwidth=2,\n",
" borderpad=4,\n",
" opacity=0.8\n",
" )\n",
"\n",
"fig = go.Figure()\n",
"add_line(fig, range_arr/1e3, grav1/1e3, 'Gravity')\n",
"add_line(fig, range_arr/1e3, grav2/1e3, 'J2', 0.2, 0, -10)\n",
"add_line(fig, range_arr/1e3, grav6/1e3, 'J6', 0.8, 0, -10)\n",
"add_line(fig, range_arr[didx]/1e3, drag_max/1e3, 'Drag Max', 0.7, 30, 0)\n",
"add_line(fig, range_arr[didx]/1e3, drag_min/1e3, 'Drag Min', 0.8, 10, 30)\n",
"add_line(fig, range_arr/1e3, moon/1e3, 'Moon', 0.8, -10, -10)\n",
"add_line(fig, range_arr/1e3, sun/1e3, 'Sun', 0.7, -10, 10)\n",
"add_line(fig, range_arr/1e3, a_radiation/1e3, 'Radiation\\nPressure', 0.3, -10, 10)\n",
"\n",
"fig.update_xaxes(type=\"log\", title='Distance from Earth Origin [km]', range=np.log10([6378.1, 50e3]))\n",
"fig.update_yaxes(type=\"log\", title='Acceleration [km/s2]', tickformat=\".1e\")\n",
"fig.update_layout(title='Satellite Forces vs Altitude', width=650, height=650)\n",
"fig.update_xaxes(showline=True, linewidth=2, linecolor='black', mirror=True)\n",
"fig.update_yaxes(showline=True, linewidth=2, linecolor='black', mirror=True)\n",
"fig.update_layout(showlegend=False, xaxis=dict(\n",
" gridcolor='#dddddd',\n",
" gridwidth=1,\n",
" ),\n",
" yaxis=dict(\n",
" gridcolor='#dddddd',\n",
" gridwidth=1,\n",
" ),)\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Difference in angle between IAU-2010 reduction and approximate calculation\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import plotly.graph_objects as go\n",
"import numpy as np\n",
"import satkit as sk\n",
"import math as m\n",
"\n",
"start = sk.time(2023, 1, 1)\n",
"end = sk.time(1970, 1, 1)\n",
"duration = end - start\n",
"timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)])\n",
"qexact = sk.frametransform.qgcrf2itrf(timearray)\n",
"qapprox = sk.frametransform.qgcrf2itrf_approx(timearray)\n",
"qdiff = np.array([q1 * q2.conj for q1, q2 in zip(qexact, qapprox)]) # type: ignore\n",
"theta = np.array([q.angle for q in qdiff])\n",
"fig = go.Figure()\n",
"fig.add_trace(go.Scatter(x=[t.datetime() for t in timearray], y=theta*180.0/m.pi*3600,\n",
" mode='lines', name='Angle Difference',\n",
" line=dict(width=1, color='black')))\n",
"fig.update_layout(title='Angle Error of Approximate ITRF to GCRF Rotation', width=650, height=512,\n",
" font=dict(size=14))\n",
"fig.update_yaxes(title='Error [arcsec]')\n",
"fig.update_xaxes(title='Year')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Polar Motion"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import plotly.graph_objects as go\n",
"import numpy as np\n",
"import satkit as sk\n",
"import math as m\n",
"\n",
"start = sk.time(2023, 1, 1)\n",
"end = sk.time(1970, 1, 1)\n",
"duration = end - start\n",
"timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, 4000)])\n",
"\n",
"dut1, xp, yp, lod, dX, dY = zip(*[sk.frametransform.earth_orientation_params(t) for t in timearray])\n",
"fig = go.Figure()\n",
"fig.add_trace(go.Scatter3d(x=xp, y=yp, z=[t.datetime() for t in timearray], mode='lines', name='Polar Motion',))\n",
"fig.update_scenes(xaxis_title='X [arcsec]', yaxis_title='Y [arcsec]', zaxis_title='Year')\n",
"fig.update_layout(title='Polar Motion', width=650, height=600, font=dict(size=12))\n",
"fig.show()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Precession / Nutation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"import math as m\n",
"import plotly.graph_objects as go\n",
"\n",
"# Precession / Nutation is the difference between the IAU-2006 GCRS frame and the CIRS frame\n",
"start = sk.time(2000, 1, 1)\n",
"end = sk.time(2020, 1, 1)\n",
"duration = end - start\n",
"N = 1000\n",
"timearray = np.array([start + sk.duration(days=x) for x in np.linspace(0, duration.days, N)])\n",
"qarray = np.array([sk.frametransform.qcirs2gcrf(t) for t in timearray])\n",
"theta = np.fromiter((q.angle for q in qarray), float)\n",
"rots = np.array([q * np.array([0.0, 0.0, 1.0]) for q in qarray])\n",
"pline = np.linspace(0,1,N)\n",
"pv = np.polyfit(pline, rots[:,0], 1)\n",
"rx0 =rots[:,0] - np.polyval(pv, pline)\n",
"\n",
"rad2asec = 180.0/m.pi * 3600\n",
"fig = go.Figure()\n",
"fig.add_trace(go.Scatter3d(x=rx0*rad2asec, y=rots[:,1]*rad2asec, z=[t.datetime() for t in timearray], mode='lines',\n",
" name='Precession/Nutation', line=dict(width=2, color='black')))\n",
"fig.update_layout(title='Precession / Nutation', width=650, height=600, font=dict(size=12))\n",
"\n",
"fig.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}