{ "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", "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", "\n", "Note: parts of this software were adapted from Cees Bassa (ASTRON);\n", " see https://github.com/cbassa/satellite_analysis\n", "```" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import datetime\n", "import numpy as np\n", "\n", "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d import Axes3D\n", "import matplotlib\n", "from matplotlib import animation, rc\n", "\n", "from astropy import time\n", "from astropy import units as u\n", "from astropy import constants as const\n", "from astropy.coordinates import EarthLocation\n", "from pycraf import satellite\n", "\n", "rc('animation', html='html5')\n", "matplotlib.rcParams['animation.embed_limit'] = 2 ** 128" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Name = Nominal Earth mass parameter\n", " Value = 398600400000000.0\n", " Uncertainty = 0.0\n", " Unit = m3 / s2\n", " Reference = IAU 2015 Resolution B 3\n", " Name = Nominal Earth equatorial radius\n", " Value = 6378100.0\n", " Uncertainty = 0.0\n", " Unit = m\n", " Reference = IAU 2015 Resolution B 3\n" ] } ], "source": [ "print(const.GM_earth, const.R_earth, sep='\\n')" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def satellite_mean_motion(altitude, mu=const.GM_earth, r_earth=const.R_earth):\n", " '''\n", " Compute mean motion of satellite at altitude in Earth's gravitational field.\n", " \n", " See https://en.wikipedia.org/wiki/Mean_motion#Formulae\n", " '''\n", " no = np.sqrt(4.0 * np.pi ** 2 * (altitude + r_earth) ** 3 / mu).to(u.day)\n", " return 1 / no" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$15.219487 \\; \\mathrm{\\frac{1}{d}}$" ], "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "satellite_mean_motion(500 * u.km)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def tle_from_orbital_parameters(\n", " sat_name, sat_nr, epoch, inclination, raan, mean_anomaly, mean_motion\n", " ):\n", " '''\n", " Generate TLE strings from orbital parameters.\n", " \n", " Note: epoch has a very strange format: first two digits are the year, next three\n", " digits are the day from beginning of year, then fraction of a day is given, e.g.\n", " 20180.25 would be 2020, day 180, 6 hours (UT?)\n", " '''\n", " \n", " # Note: RAAN = right ascention (or longitude) of ascending node\n", " \n", " def checksum(line):\n", " s = 0\n", " for c in line[:-1]:\n", " if c.isdigit():\n", " s += int(c)\n", " if c == \"-\":\n", " s += 1\n", " return '{:s}{:1d}'.format(line[:-1], s % 10)\n", " \n", " tle0 = sat_name\n", " tle1 = checksum(\n", " '1 {:05d}U 20001A {:14.8f} .00000000 00000-0 50000-4 '\n", " '0 0X'.format(sat_nr, epoch)\n", " )\n", " tle2 = checksum(\n", " '2 {:05d} {:8.4f} {:8.4f} 0001000 0.0000 {:8.4f} '\n", " '{:11.8f} 0X'.format(\n", " sat_nr, inclination.to_value(u.deg), raan.to_value(u.deg),\n", " mean_anomaly.to_value(u.deg), mean_motion.to_value(1 / u.day)\n", " ))\n", " \n", " return '\\n'.join([tle0, tle1, tle2])\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "TEST\n", "1 00099U 20001A 19050.10000000 .00000000 00000-0 50000-4 0 09\n", "2 00099 10.0000 20.0000 0001000 0.0000 45.0000 15.21948688 05\n" ] } ], "source": [ "my_sat_tle = tle_from_orbital_parameters(\n", " 'TEST', 99, 19050.1, 10 * u.deg, 20 * u.deg, 45 * u.deg,\n", " satellite_mean_motion(500 * u.km)\n", " )\n", "print(my_sat_tle)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "('TEST', )" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "my_sat = satellite.get_sat(my_sat_tle)\n", "my_sat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Starlink constellation" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "altitudes = np.array([550, 1110, 1130, 1275, 1325, 345.6, 340.8, 335.9]) * u.km\n", "inclinations = np.array([53.0, 53.8, 74.0, 81.0, 70.0, 53.0, 48.0, 42.0]) * u.deg\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": 9, "metadata": {}, "outputs": [], "source": [ "def create_constellation(altitudes, inclinations, nplanes, sats_per_plane):\n", " \n", " my_sat_tles = []\n", " sat_nr = 80000\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) * u.deg\n", " raans = np.random.uniform(0, 360, n) * u.deg\n", " else:\n", " mas = np.linspace(0.0, 360.0, s, endpoint=False) * u.deg\n", " mas += np.random.uniform(0, 360, 1) * u.deg\n", " raans = np.linspace(0.0, 360.0, n, endpoint=False) * u.deg\n", " mas, raans = np.meshgrid(mas, raans)\n", " mas, raans = mas.flatten(), raans.flatten()\n", " \n", " mm = satellite_mean_motion(alt)\n", " for ma, raan in zip(mas, raans):\n", " my_sat_tles.append(\n", " tle_from_orbital_parameters(\n", " 'TEST {:d}'.format(sat_nr), sat_nr, 19050.1,\n", " inc, raan, ma, mm\n", " ))\n", " sat_nr += 1\n", " \n", " return my_sat_tles" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "my_sat_tles = create_constellation(altitudes, inclinations, nplanes, sats_per_plane)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "11927" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(my_sat_tles)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "my_sats = [satellite.get_sat(sat_tle) for sat_tle in my_sat_tles]" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "def _propagate(sat, dt):\n", " '''\n", " True equator mean equinox (TEME) position from `sgp4` at given time.\n", " Parameters\n", " ----------\n", " sat : `sgp4.io.Satellite` instance\n", " Satellite object filled from TLE\n", " dt : `~datetime.datetime`\n", " Time\n", " Returns\n", " -------\n", " xs, ys, zs : float\n", " TEME (=True equator mean equinox) position of satellite [km]\n", " '''\n", "\n", " # pos [km], vel [km/s]\n", " position, velocity = sat.propagate(\n", " dt.year, dt.month, dt.day,\n", " dt.hour, dt.minute, dt.second + dt.microsecond / 1e6\n", " )\n", "\n", " if position is None:\n", " raise ValueError('Satellite propagation error')\n", "\n", " return position\n", "\n", "\n", "vec_propagate = np.vectorize(\n", " _propagate, excluded=['sat'], otypes=[np.float64] * 3\n", " )" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Start time: 2019-10-23 20:35:23.514158\n" ] } ], "source": [ "start_time = time.Time(datetime.datetime.utcnow(), 'utc')\n", "print('Start time:', start_time)\n", "td = time.TimeDelta(np.arange(0, 600, 1), format='sec')\n", "times = start_time + td # 10 min in steps of 1 s" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "pos_eci = np.array([vec_propagate(sat[1], times.datetime) for sat in my_sats])\n", "# pos_eci" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plim = 8000\n", "\n", "fig = plt.figure(figsize=(12, 9))\n", "ax = Axes3D(fig)\n", "ax.view_init(azim=60, elev=30)\n", "ax.set_xlim((-plim, plim))\n", "ax.set_ylim((-plim, plim))\n", "ax.set_zlim((-plim, plim))\n", "ax.set_aspect('equal')\n", "ax.set_xlabel('x [km]')\n", "ax.set_ylabel('y [km]')\n", "ax.set_zlabel('z [km]')\n", "\n", "rads = np.sqrt(pos_eci[:, 0, 0] ** 2 + pos_eci[:, 1, 0] ** 2 + pos_eci[:, 2, 0] ** 2)\n", "points = ax.scatter(\n", " pos_eci[:, 0, 0], pos_eci[:, 1, 0], pos_eci[:, 2, 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(times[0].datetime), loc='center', fontsize=20)\n", "\n", "def init():\n", " points._offsets3d = (pos_eci[:, 0, 0], pos_eci[:, 1, 0], pos_eci[:, 2, 0])\n", " title.set_text('{:%y/%m/%d %H:%M:%S}'.format(times[0].datetime))\n", " return points, title\n", "\n", "def animate(i):\n", " points._offsets3d = (pos_eci[:, 0, i], pos_eci[:, 1, i], pos_eci[:, 2, i])\n", " title.set_text('{:%y/%m/%d %H:%M:%S}'.format(times[i].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=pos_eci.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": 17, "metadata": {}, "outputs": [], "source": [ "obs_loc_parkes = EarthLocation(148.25738, -32.9933, 414.8)\n", "obs_loc_effbg = EarthLocation(6.88375, 50.525, 366.)\n", "\n", "# create a SatelliteObserver instance\n", "sat_obs_parkes = satellite.SatelliteObserver(obs_loc_parkes)\n", "sat_obs_effbg = satellite.SatelliteObserver(obs_loc_effbg)\n", "pos_horiz_parkes = np.array([\n", " sat_obs_parkes.azel_from_sat(sat[1], times)\n", " for sat in my_sats\n", " ])\n", "pos_horiz_effbg = np.array([\n", " sat_obs_effbg.azel_from_sat(sat[1], times)\n", " for sat in my_sats\n", " ])" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "vmin, vmax = np.log10(100), np.log10(10000)\n", "\n", "fig = plt.figure(figsize=(12, 7))\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", " pos_horiz_parkes[:, 0, 0], pos_horiz_parkes[:, 1, 0],\n", " c=np.log10(pos_horiz_parkes[:, 2, 0]),\n", " cmap='viridis', vmin=vmin, vmax=vmax,\n", " )\n", "points2 = ax2.scatter(\n", " pos_horiz_effbg[:, 0, 0], pos_horiz_effbg[:, 1, 0],\n", " c=np.log10(pos_horiz_effbg[:, 2, 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(times[0].datetime),\n", " fontsize=15, ha='right'\n", " )\n", "\n", "def init():\n", " points1.set_offsets(pos_horiz_parkes[:, 0:2, 0])\n", " points1.set_array(np.log10(pos_horiz_parkes[:, 2, 0]))\n", " points2.set_offsets(pos_horiz_effbg[:, 0:2, 0])\n", " points2.set_array(np.log10(pos_horiz_effbg[:, 2, 0]))\n", " title.set_text('{:%y/%m/%d %H:%M:%S}'.format(times[0].datetime))\n", " return points, title\n", "\n", "def animate(i):\n", " points1.set_offsets(pos_horiz_parkes[:, 0:2, i])\n", " points1.set_array(np.log10(pos_horiz_parkes[:, 2, i]))\n", " points2.set_offsets(pos_horiz_effbg[:, 0:2, i])\n", " points2.set_array(np.log10(pos_horiz_effbg[:, 2, i]))\n", " title.set_text('{:%y/%m/%d %H:%M:%S}'.format(times[i].datetime))\n", " return points, title\n", "\n", "anim = animation.FuncAnimation(\n", " fig, animate, init_func=init, frames=pos_horiz_parkes.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": 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 }