SpaceX/Starlink satellite constellation: orbit visualization.
Copyright (C) 2019+  Benjamin Winkel ([email protected])

This program is free software: you can redistribute it and/or modify
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <https://www.gnu.org/licenses/>.

Note: parts of this software were adapted from Cees Bassa (ASTRON);
see https://github.com/cbassa/satellite_analysis
In [1]:
import datetime
import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib
from matplotlib import animation, rc

from astropy import time
from astropy import units as u
from astropy import constants as const
from astropy.coordinates import EarthLocation
from pycraf import satellite

rc('animation', html='html5')
matplotlib.rcParams['animation.embed_limit'] = 2 ** 128

In [2]:
print(const.GM_earth, const.R_earth, sep='\n')

  Name   = Nominal Earth mass parameter
Value  = 398600400000000.0
Uncertainty  = 0.0
Unit  = m3 / s2
Reference = IAU 2015 Resolution B 3
Name   = Nominal Earth equatorial radius
Value  = 6378100.0
Uncertainty  = 0.0
Unit  = m
Reference = IAU 2015 Resolution B 3

In [3]:
def satellite_mean_motion(altitude, mu=const.GM_earth, r_earth=const.R_earth):
'''
Compute mean motion of satellite at altitude in Earth's gravitational field.

See https://en.wikipedia.org/wiki/Mean_motion#Formulae
'''
no = np.sqrt(4.0 * np.pi ** 2 * (altitude + r_earth) ** 3 / mu).to(u.day)
return 1 / no

In [4]:
satellite_mean_motion(500 * u.km)

Out[4]:
$15.219487 \; \mathrm{\frac{1}{d}}$
In [5]:
def tle_from_orbital_parameters(
sat_name, sat_nr, epoch, inclination, raan, mean_anomaly, mean_motion
):
'''
Generate TLE strings from orbital parameters.

Note: epoch has a very strange format: first two digits are the year, next three
digits are the day from beginning of year, then fraction of a day is given, e.g.
20180.25 would be 2020, day 180, 6 hours (UT?)
'''

# Note: RAAN = right ascention (or longitude) of ascending node

def checksum(line):
s = 0
for c in line[:-1]:
if c.isdigit():
s += int(c)
if c == "-":
s += 1
return '{:s}{:1d}'.format(line[:-1], s % 10)

tle0 = sat_name
tle1 = checksum(
'1 {:05d}U 20001A   {:14.8f}  .00000000  00000-0  50000-4 '
'0    0X'.format(sat_nr, epoch)
)
tle2 = checksum(
'2 {:05d} {:8.4f} {:8.4f} 0001000   0.0000 {:8.4f} '
'{:11.8f}    0X'.format(
sat_nr, inclination.to_value(u.deg), raan.to_value(u.deg),
mean_anomaly.to_value(u.deg), mean_motion.to_value(1 / u.day)
))

return '\n'.join([tle0, tle1, tle2])

In [6]:
my_sat_tle = tle_from_orbital_parameters(
'TEST', 99, 19050.1, 10 * u.deg, 20 * u.deg, 45 * u.deg,
satellite_mean_motion(500 * u.km)
)
print(my_sat_tle)

TEST
1 00099U 20001A   19050.10000000  .00000000  00000-0  50000-4 0    09
2 00099  10.0000  20.0000 0001000   0.0000  45.0000 15.21948688    05

In [7]:
my_sat = satellite.get_sat(my_sat_tle)
my_sat

Out[7]:
('TEST', <sgp4.model.Satellite at 0x7f14529f7b00>)

In [8]:
altitudes = np.array([550, 1110, 1130, 1275, 1325, 345.6, 340.8, 335.9]) * u.km
inclinations = np.array([53.0, 53.8, 74.0, 81.0, 70.0, 53.0, 48.0, 42.0]) * u.deg
nplanes = np.array([72, 32, 8, 5, 6, 2547, 2478, 2493])
sats_per_plane = np.array([22, 50, 50, 75, 75, 1, 1, 1])

In [9]:
def create_constellation(altitudes, inclinations, nplanes, sats_per_plane):

my_sat_tles = []
sat_nr = 80000
for alt, inc, n, s in zip(
altitudes, inclinations, nplanes, sats_per_plane
):

if s == 1:
# random placement for lower orbits
mas = np.random.uniform(0, 360, n) * u.deg
raans = np.random.uniform(0, 360, n) * u.deg
else:
mas = np.linspace(0.0, 360.0, s, endpoint=False) * u.deg
mas += np.random.uniform(0, 360, 1) * u.deg
raans = np.linspace(0.0, 360.0, n, endpoint=False) * u.deg
mas, raans = np.meshgrid(mas, raans)
mas, raans = mas.flatten(), raans.flatten()

mm = satellite_mean_motion(alt)
for ma, raan in zip(mas, raans):
my_sat_tles.append(
tle_from_orbital_parameters(
'TEST {:d}'.format(sat_nr), sat_nr, 19050.1,
inc, raan, ma, mm
))
sat_nr += 1

return my_sat_tles

In [10]:
my_sat_tles = create_constellation(altitudes, inclinations, nplanes, sats_per_plane)

In [11]:
len(my_sat_tles)

Out[11]:
11927
In [12]:
my_sats = [satellite.get_sat(sat_tle) for sat_tle in my_sat_tles]

In [13]:
def _propagate(sat, dt):
'''
True equator mean equinox (TEME) position from sgp4 at given time.
Parameters
----------
sat : sgp4.io.Satellite instance
Satellite object filled from TLE
dt : ~datetime.datetime
Time
Returns
-------
xs, ys, zs : float
TEME (=True equator mean equinox) position of satellite [km]
'''

# pos [km], vel [km/s]
position, velocity = sat.propagate(
dt.year, dt.month, dt.day,
dt.hour, dt.minute, dt.second + dt.microsecond / 1e6
)

if position is None:
raise ValueError('Satellite propagation error')

return position

vec_propagate = np.vectorize(
_propagate, excluded=['sat'], otypes=[np.float64] * 3
)

In [14]:
start_time = time.Time(datetime.datetime.utcnow(), 'utc')
print('Start time:', start_time)
td = time.TimeDelta(np.arange(0, 600, 1), format='sec')
times = start_time + td  # 10 min in steps of 1 s

Start time: 2019-10-23 20:35:23.514158

In [15]:
pos_eci = np.array([vec_propagate(sat[1], times.datetime) for sat in my_sats])
# pos_eci

In [16]:
plim = 8000

fig = plt.figure(figsize=(12, 9))
ax = Axes3D(fig)
ax.view_init(azim=60, elev=30)
ax.set_xlim((-plim, plim))
ax.set_ylim((-plim, plim))
ax.set_zlim((-plim, plim))
ax.set_aspect('equal')
ax.set_xlabel('x [km]')
ax.set_ylabel('y [km]')
ax.set_zlabel('z [km]')

rads = np.sqrt(pos_eci[:, 0, 0] ** 2 + pos_eci[:, 1, 0] ** 2 + pos_eci[:, 2, 0] ** 2)
points = ax.scatter(
pos_eci[:, 0, 0], pos_eci[:, 1, 0], pos_eci[:, 2, 0],
)
title = ax.set_title('{:%y/%m/%d %H:%M:%S}'.format(times[0].datetime), loc='center', fontsize=20)

def init():
points._offsets3d = (pos_eci[:, 0, 0], pos_eci[:, 1, 0], pos_eci[:, 2, 0])
title.set_text('{:%y/%m/%d %H:%M:%S}'.format(times[0].datetime))
return points, title

def animate(i):
points._offsets3d = (pos_eci[:, 0, i], pos_eci[:, 1, i], pos_eci[:, 2, i])
title.set_text('{:%y/%m/%d %H:%M:%S}'.format(times[i].datetime))
return points, title

# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(
fig, animate, init_func=init, frames=pos_eci.shape[2], interval=20, blit=True
)

# this takes a while!
plt.close(anim._fig)
anim

Out[16]: