{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## License\n",
"```\n",
"SpaceX/Starlink satellite constellation: orbit visualization.\n",
"Copyright (C) 2019+ Benjamin Winkel (bwinkel@mpifr.de)\n",
"\n",
"Note: parts of this software were adapted from Cees Bassa (ASTRON);\n",
" see https://github.com/cbassa/satellite_analysis\n",
"\n",
"This program is free software: you can redistribute it and/or modify\n",
"it under the terms of the GNU General Public License as published by\n",
"the Free Software Foundation, either version 3 of the License, or\n",
"(at your option) any later version.\n",
"\n",
"This program is distributed in the hope that it will be useful,\n",
"but WITHOUT ANY WARRANTY; without even the implied warranty of\n",
"MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n",
"GNU General Public License for more details.\n",
"\n",
"You should have received a copy of the GNU General Public License\n",
"along with this program. If not, see .\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"# This import registers the 3D projection, but is otherwise unused.\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"# for animations, you may need to install \"ffmpeg\" and/or \"imagemagick\"\n",
"from matplotlib import animation, rc\n",
"\n",
"import cysgp4\n",
"\n",
"rc('animation', html='html5')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you plan to make big movies (i.e., large file sizes), you may need to start Jupyter with the command line option ``--NotebookApp.iopub_data_rate_limit=1.0e10`` and do the following:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# matplotlib.rcParams['animation.embed_limit'] = 2 ** 128"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From the FCC filing of the Starlink constellation, we know at which altitudes, inclinations etc. the satellites will be operated. From this information, we need to produce TLEs, in order to process them with the SGP4 software. For example:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# want epoch for the following time\n",
"mjd_epoch = 58813.5\n",
"pydt = cysgp4.PyDateTime.from_mjd(mjd_epoch)\n",
"pydt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note, the TLE uses a special format for the date/time:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"19330.5"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# the TLE uses a special format for the date/time:\n",
"pydt.tle_epoch"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But the `tle_linestrings_from_orbital_parameters` works with standard MJD:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"('MYSAT',\n",
" '1 00001U 20001A 19330.50000000 .00000000 00000-0 50000-4 0 07',\n",
" '2 00001 10.0000 35.0000 0001000 0.0000 112.0000 9.55934723 04')"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sat_name, sat_nr = 'MYSAT', 1\n",
"alt_km = 3000. # satellite altitude\n",
"mean_motion = cysgp4.satellite_mean_motion(alt_km)\n",
"inclination = 10. # deg\n",
"raan = 35. # deg\n",
"eccentricity = 0.0001\n",
"argument_of_perigee = 0. # deg\n",
"mean_anomaly = 112. # deg\n",
"\n",
"tle_tuple = cysgp4.tle_linestrings_from_orbital_parameters(\n",
" sat_name,\n",
" sat_nr,\n",
" mjd_epoch,\n",
" inclination,\n",
" raan,\n",
" eccentricity,\n",
" argument_of_perigee,\n",
" mean_anomaly,\n",
" mean_motion\n",
" )\n",
"\n",
"tle_tuple"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To check, that this is working as intended:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"my_tle = cysgp4.PyTle(*tle_tuple)\n",
"# want position at a given time, which can of course differ from the epoch\n",
"obs_dt = cysgp4.PyDateTime.from_mjd(58814.23)\n",
"my_sat = cysgp4.Satellite(my_tle)\n",
"my_sat.pydt = obs_dt\n",
"my_sat.eci_pos() # in Geodetic frame"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now lets define the Starlink constellation."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"altitudes = np.array([550, 1110, 1130, 1275, 1325, 345.6, 340.8, 335.9])\n",
"inclinations = np.array([53.0, 53.8, 74.0, 81.0, 70.0, 53.0, 48.0, 42.0])\n",
"nplanes = np.array([72, 32, 8, 5, 6, 2547, 2478, 2493])\n",
"sats_per_plane = np.array([22, 50, 50, 75, 75, 1, 1, 1])"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total number of satellites 11927\n"
]
}
],
"source": [
"print('total number of satellites', np.sum(nplanes * sats_per_plane))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note, that Starlink now even seems to ask for permission to launch almost 40000 satellites."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def create_constellation(mjd_epoch, altitudes, inclinations, nplanes, sats_per_plane):\n",
" \n",
" my_sat_tles = []\n",
" sat_nr = 1\n",
" for alt, inc, n, s in zip(\n",
" altitudes, inclinations, nplanes, sats_per_plane\n",
" ):\n",
" \n",
" if s == 1:\n",
" # random placement for lower orbits\n",
" mas = np.random.uniform(0, 360, n)\n",
" raans = np.random.uniform(0, 360, n)\n",
" else:\n",
" mas = np.linspace(0.0, 360.0, s, endpoint=False)\n",
" mas += np.random.uniform(0, 360, 1)\n",
" raans = np.linspace(0.0, 360.0, n, endpoint=False)\n",
" mas, raans = np.meshgrid(mas, raans)\n",
" mas, raans = mas.flatten(), raans.flatten()\n",
" \n",
" mm = cysgp4.satellite_mean_motion(alt)\n",
" for ma, raan in zip(mas, raans):\n",
" my_sat_tles.append(\n",
" cysgp4.tle_linestrings_from_orbital_parameters(\n",
" 'TEST {:d}'.format(sat_nr), sat_nr, mjd_epoch,\n",
" inc, raan, 0.001, 0., ma, mm\n",
" ))\n",
" \n",
" sat_nr += 1\n",
" \n",
" return my_sat_tles"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"11927"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"starlink_tle_tuples = create_constellation(\n",
" mjd_epoch, altitudes, inclinations, nplanes, sats_per_plane\n",
" )\n",
"len(starlink_tle_tuples)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"starlink_tles = np.array([\n",
" cysgp4.PyTle(*tle)\n",
" for tle in starlink_tle_tuples\n",
" ])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now plot similarly as in the first notebook."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"start_mjd = mjd_epoch\n",
"td = np.arange(0, 600, 5) / 86400. # 1 d in steps of 10 s\n",
"mjds = start_mjd + td"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# Effelsberg 100-m radio telescope\n",
"effbg_observer = cysgp4.PyObserver(6.88375, 50.525, 0.366)\n",
"# Parkes telescope (\"The Dish\")\n",
"parkes_observer = cysgp4.PyObserver(148.25738, -32.9933, 414.8)\n",
"observers = np.array([effbg_observer, parkes_observer])"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"result = cysgp4.propagate_many(\n",
" mjds[np.newaxis, np.newaxis, :],\n",
" starlink_tles[:, np.newaxis, np.newaxis],\n",
" observers[np.newaxis, :, np.newaxis]\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(120, 11927, 2, (11927, 2, 120, 3), (11927, 2, 120, 4))"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eci_pos = result['eci_pos']\n",
"topo_pos = result['topo']\n",
"len(mjds), len(starlink_tles), len(observers), eci_pos.shape, topo_pos.shape"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"eci_pos_x, eci_pos_y, eci_pos_z = (eci_pos[..., i] for i in range(3))\n",
"topo_pos_az, topo_pos_el, topo_pos_dist, _ = (topo_pos[..., i] for i in range(4))\n",
"topo_pos_az = (topo_pos_az + 180.) % 360. - 180."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"my_time = cysgp4.PyDateTime()\n",
"my_time.mjd = mjds[0]\n",
"plim = 8000\n",
"\n",
"# The figure size should make such that one gets a nice pixel canvas\n",
"# that fits the standard movie sizes (at given dpi):\n",
"# 854 x 480 (480p) --> figsize=(8.54, 4.8), dpi=100\n",
"# 1280 x 720 (720p) --> figsize=(12.8, 7.2), dpi=100\n",
"# 1920 x 1080 (1080p) --> figsize=(12.8, 7.2), dpi=150\n",
"# 3840 x 2160 (4K) --> figsize=(12.8, 7.2), dpi=300\n",
"# so basically, divide desired width and height with dpi\n",
"# (beware, 4K videos get large and need a lot of RAM!)\n",
"fig = plt.figure(figsize=(12.8, 7.2), dpi=100)\n",
"ax = fig.add_subplot(111, projection='3d') # 3D axes\n",
"ax.view_init(azim=60, elev=30)\n",
"# Aspect ratio is not implemented; \n",
"# see https://github.com/matplotlib/matplotlib/issues/1077/\n",
"# ax.set_aspect('equal')\n",
"# need to manually stretch to make it approx. right\n",
"ax.set_xlim((-2 * plim, 2 * plim))\n",
"ax.set_ylim((-2 * plim, 2 * plim))\n",
"ax.set_zlim((-plim, plim))\n",
"# ax.auto_scale_xyz([-2 * plim, 2 * plim], [-2 * plim, 2 * plim], [-plim, plim])\n",
"# axisEqual3D(ax)\n",
"ax.set_xlabel('x [km]')\n",
"ax.set_ylabel('y [km]')\n",
"ax.set_zlabel('z [km]')\n",
"\n",
"rads = np.sqrt(eci_pos_x[:, 0, 0] ** 2 + eci_pos_y[:, 0, 0] ** 2 + eci_pos_z[:, 0, 0] ** 2)\n",
"points = ax.scatter(\n",
" eci_pos_x[:, 0, 0], eci_pos_y[:, 0, 0], eci_pos_z[:, 0, 0],\n",
" c=rads, s=1, vmin=rads.min(), vmax=rads.max(), marker='o'\n",
" )\n",
"title = ax.set_title('{:%y/%m/%d %H:%M:%S}'.format(my_time.datetime), loc='center', fontsize=20)\n",
"\n",
"def init():\n",
" points._offsets3d = (eci_pos_x[:, 0, 0], eci_pos_y[:, 0, 0], eci_pos_z[:, 0, 0])\n",
" my_time.mjd = mjds[0]\n",
" title.set_text('{:%y/%m/%d %H:%M:%S}'.format(my_time.datetime))\n",
" return points, title\n",
"\n",
"def animate(i):\n",
" points._offsets3d = (eci_pos_x[:, 0, i], eci_pos_y[:, 0, i], eci_pos_z[:, 0, i])\n",
" my_time.mjd = mjds[i]\n",
" title.set_text('{:%y/%m/%d %H:%M:%S}'.format(my_time.datetime))\n",
" return points, title\n",
"\n",
"# call the animator. blit=True means only re-draw the parts that have changed.\n",
"anim = animation.FuncAnimation(\n",
" fig, animate, init_func=init, frames=eci_pos_x.shape[2], interval=20, blit=True\n",
" )\n",
"\n",
"# this takes a while!\n",
"plt.close(anim._fig)\n",
"# anim"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"FFwriter = animation.FFMpegWriter(\n",
" fps=30, bitrate=8000,\n",
" extra_args=['-vcodec', 'libx264'],\n",
" )\n",
"anim.save('starlink_constellation_eci_720p.mp4', writer=FFwriter)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"my_time = cysgp4.PyDateTime()\n",
"my_time.mjd = mjds[0]\n",
"vmin, vmax = np.log10(100), np.log10(10000)\n",
"\n",
"fig = plt.figure(figsize=(12.8, 7.2), dpi=100)\n",
"ax1 = fig.add_axes((0.1, 0.5, 0.8, 0.35))\n",
"ax2 = fig.add_axes((0.1, 0.1, 0.8, 0.35))\n",
"cax = fig.add_axes((0.91, 0.2, 0.02, 0.5))\n",
"ax2.set_xlabel('Azimuth [deg]')\n",
"ax1.set_ylabel('Elevation [deg]')\n",
"for ax in [ax1, ax2]:\n",
" ax.set_xlim((-180, 180))\n",
" ax.set_ylim((0, 90))\n",
" ax.set_xticks(range(-150, 180, 30))\n",
" ax.set_yticks(range(0, 91, 30))\n",
" ax.set_aspect('equal')\n",
" ax.grid()\n",
"\n",
"points1 = ax1.scatter(\n",
" topo_pos_az[:, 1, 0], topo_pos_el[:, 1, 0],\n",
" c=np.log10(topo_pos_dist[:, 1, 0]),\n",
" cmap='viridis', vmin=vmin, vmax=vmax,\n",
" )\n",
"points2 = ax2.scatter(\n",
" topo_pos_az[:, 0, 0], topo_pos_el[:, 0, 0],\n",
" c=np.log10(topo_pos_dist[:, 0, 0]),\n",
" cmap='viridis', vmin=vmin, vmax=vmax,\n",
" )\n",
"cbar = fig.colorbar(points1, cax=cax)\n",
"cbar.set_label('Distance (km)')\n",
"cbar.set_ticks([2, 3, 4])\n",
"cbar.set_ticklabels([100, 1000, 10000])\n",
"\n",
"ax1.text(-170, 75, 'Parkes 64-m', fontsize=16)\n",
"ax2.text(-170, 75, 'Effelsberg 100-m', fontsize=16)\n",
"title = ax1.text(\n",
" 174, 75, '{:%y/%m/%d %H:%M:%S}'.format(my_time.datetime),\n",
" fontsize=15, ha='right'\n",
" )\n",
"\n",
"def init():\n",
" points1.set_offsets(np.column_stack([topo_pos_az[:, 1, 0], topo_pos_el[:, 1, 0]]))\n",
" points1.set_array(np.log10(topo_pos_dist[:, 1, 0]))\n",
" points2.set_offsets(np.column_stack([topo_pos_az[:, 0, 0], topo_pos_el[:, 0, 0]]))\n",
" points2.set_array(np.log10(topo_pos_dist[:, 0, 0]))\n",
" my_time.mjd = mjds[0]\n",
" title.set_text('{:%y/%m/%d %H:%M:%S}'.format(my_time.datetime))\n",
" return points, title\n",
"\n",
"def animate(i):\n",
" points1.set_offsets(np.column_stack([topo_pos_az[:, 1, i], topo_pos_el[:, 1, i]]))\n",
" points1.set_array(np.log10(topo_pos_dist[:, 1, i]))\n",
" points2.set_offsets(np.column_stack([topo_pos_az[:, 0, i], topo_pos_el[:, 0, i]]))\n",
" points2.set_array(np.log10(topo_pos_dist[:, 0, i]))\n",
" my_time.mjd = mjds[i]\n",
" title.set_text('{:%y/%m/%d %H:%M:%S}'.format(my_time.datetime))\n",
" return points, title\n",
"\n",
"anim = animation.FuncAnimation(\n",
" fig, animate, init_func=init, frames=topo_pos_az.shape[2], interval=20, blit=True\n",
" )\n",
"\n",
"# this takes a while!\n",
"plt.close(anim._fig)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"FFwriter = animation.FFMpegWriter(\n",
" fps=30, bitrate=8000,\n",
" extra_args=['-vcodec', 'libx264'],\n",
" )\n",
"anim.save('starlink_constellation_horizon_720p.mp4', writer=FFwriter)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}