Notebook

License¶

SpaceX/Starlink satellite constellation: orbit visualization.
Copyright (C) 2019+  Benjamin Winkel (bwinkel@mpifr.de)

Note: parts of this software were adapted from Cees Bassa (ASTRON);
      see https://github.com/cbassa/satellite_analysis

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
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/>.

In [1]:

import numpy as np

import matplotlib
import matplotlib.pyplot as plt
# This import registers the 3D projection, but is otherwise unused.
from mpl_toolkits.mplot3d import Axes3D
# for animations, you may need to install "ffmpeg" and/or "imagemagick"
from matplotlib import animation, rc

import cysgp4

rc('animation', html='html5')

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:

In [2]:

# matplotlib.rcParams['animation.embed_limit'] = 2 ** 128

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:

In [3]:

# want epoch for the following time
mjd_epoch = 58813.5
pydt = cysgp4.PyDateTime.from_mjd(mjd_epoch)
pydt

Out[3]:

<PyDateTime: 2019-11-26 12:00:00.000000 UTC>

Note, the TLE uses a special format for the date/time:

In [4]:

# the TLE uses a special format for the date/time:
pydt.tle_epoch

Out[4]:

19330.5

But the tle_linestrings_from_orbital_parameters works with standard MJD:

In [5]:

sat_name, sat_nr = 'MYSAT', 1
alt_km = 3000.  # satellite altitude
mean_motion = cysgp4.satellite_mean_motion(alt_km)
inclination = 10.  # deg
raan = 35.  # deg
eccentricity = 0.0001
argument_of_perigee = 0.  # deg
mean_anomaly = 112.  # deg

tle_tuple = cysgp4.tle_linestrings_from_orbital_parameters(
    sat_name,
    sat_nr,
    mjd_epoch,
    inclination,
    raan,
    eccentricity,
    argument_of_perigee,
    mean_anomaly,
    mean_motion
    )

tle_tuple

Out[5]:

('MYSAT',
 '1 00001U 20001A   19330.50000000  .00000000  00000-0  50000-4 0    07',
 '2 00001  10.0000  35.0000 0001000   0.0000 112.0000  9.55934723    04')

To check, that this is working as intended:

In [6]:

my_tle = cysgp4.PyTle(*tle_tuple)
# want position at a given time, which can of course differ from the epoch
obs_dt = cysgp4.PyDateTime.from_mjd(58814.23)
my_sat = cysgp4.Satellite(my_tle)
my_sat.pydt = obs_dt
my_sat.eci_pos()  # in Geodetic frame

Out[6]:

<PyEci: -7.4040d, 9.5532d, 2997.3144km 2019-11-27 05:31:12.000000 UTC>

Now lets define the Starlink constellation.

In [7]:

altitudes = np.array([550, 1110, 1130, 1275, 1325, 345.6, 340.8, 335.9])
inclinations = np.array([53.0, 53.8, 74.0, 81.0, 70.0, 53.0, 48.0, 42.0])
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 [8]:

print('total number of satellites', np.sum(nplanes * sats_per_plane))

total number of satellites 11927

Note, that Starlink now even seems to ask for permission to launch almost 40000 satellites.

In [9]:

def create_constellation(mjd_epoch, altitudes, inclinations, nplanes, sats_per_plane):
    
    my_sat_tles = []
    sat_nr = 1
    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)
            raans = np.random.uniform(0, 360, n)
        else:
            mas = np.linspace(0.0, 360.0, s, endpoint=False)
            mas += np.random.uniform(0, 360, 1)
            raans = np.linspace(0.0, 360.0, n, endpoint=False)
            mas, raans = np.meshgrid(mas, raans)
            mas, raans = mas.flatten(), raans.flatten()
        
        mm = cysgp4.satellite_mean_motion(alt)
        for ma, raan in zip(mas, raans):
            my_sat_tles.append(
                cysgp4.tle_linestrings_from_orbital_parameters(
                    'TEST {:d}'.format(sat_nr), sat_nr, mjd_epoch,
                    inc, raan, 0.001, 0., ma, mm
                    ))
                
            sat_nr += 1
    
    return my_sat_tles

In [10]:

starlink_tle_tuples = create_constellation(
    mjd_epoch, altitudes, inclinations, nplanes, sats_per_plane
    )
len(starlink_tle_tuples)

Out[10]:

In [11]:

starlink_tles = np.array([
    cysgp4.PyTle(*tle)
    for tle in starlink_tle_tuples
    ])

Now plot similarly as in the first notebook.

In [12]:

start_mjd = mjd_epoch
td = np.arange(0, 600, 5) / 86400.  # 1 d in steps of 10 s
mjds = start_mjd + td

In [13]:

# Effelsberg 100-m radio telescope
effbg_observer = cysgp4.PyObserver(6.88375, 50.525, 0.366)
# Parkes telescope ("The Dish")
parkes_observer = cysgp4.PyObserver(148.25738, -32.9933, 414.8)
observers = np.array([effbg_observer, parkes_observer])

In [14]:

result = cysgp4.propagate_many(
    mjds[np.newaxis, np.newaxis, :],
    starlink_tles[:, np.newaxis, np.newaxis],
    observers[np.newaxis, :, np.newaxis]
    )

In [15]:

eci_pos = result['eci_pos']
topo_pos = result['topo']
len(mjds), len(starlink_tles), len(observers), eci_pos.shape, topo_pos.shape

Out[15]:

(120, 11927, 2, (11927, 2, 120, 3), (11927, 2, 120, 4))

In [16]:

eci_pos_x, eci_pos_y, eci_pos_z = (eci_pos[..., i] for i in range(3))
topo_pos_az, topo_pos_el, topo_pos_dist, _ = (topo_pos[..., i] for i in range(4))
topo_pos_az = (topo_pos_az + 180.) % 360. - 180.

In [17]:

my_time = cysgp4.PyDateTime()
my_time.mjd = mjds[0]
plim = 8000

# The figure size should make such that one gets a nice pixel canvas
# that fits the standard movie sizes (at given dpi):
#    854 x  480  (480p) --> figsize=(8.54, 4.8), dpi=100
#   1280 x  720  (720p) --> figsize=(12.8, 7.2), dpi=100
#   1920 x 1080 (1080p) --> figsize=(12.8, 7.2), dpi=150
#   3840 x 2160    (4K) --> figsize=(12.8, 7.2), dpi=300
# so basically, divide desired width and height with dpi
# (beware, 4K videos get large and need a lot of RAM!)
fig = plt.figure(figsize=(12.8, 7.2), dpi=100)
ax = fig.add_subplot(111, projection='3d')  # 3D axes
ax.view_init(azim=60, elev=30)
# Aspect ratio is  not implemented; 
# see https://github.com/matplotlib/matplotlib/issues/1077/
# ax.set_aspect('equal')
# need to manually stretch to make it approx. right
ax.set_xlim((-2 * plim, 2 * plim))
ax.set_ylim((-2 * plim, 2 * plim))
ax.set_zlim((-plim, plim))
# ax.auto_scale_xyz([-2 * plim, 2 * plim], [-2 * plim, 2 * plim], [-plim, plim])
# axisEqual3D(ax)
ax.set_xlabel('x [km]')
ax.set_ylabel('y [km]')
ax.set_zlabel('z [km]')

rads = np.sqrt(eci_pos_x[:, 0, 0] ** 2 + eci_pos_y[:, 0, 0] ** 2 + eci_pos_z[:, 0, 0] ** 2)
points = ax.scatter(
    eci_pos_x[:, 0, 0], eci_pos_y[:, 0, 0], eci_pos_z[:, 0, 0],
    c=rads, s=1, vmin=rads.min(), vmax=rads.max(), marker='o'
    )
title = ax.set_title('{:%y/%m/%d %H:%M:%S}'.format(my_time.datetime), loc='center', fontsize=20)

def init():
    points._offsets3d = (eci_pos_x[:, 0, 0], eci_pos_y[:, 0, 0], eci_pos_z[:, 0, 0])
    my_time.mjd = mjds[0]
    title.set_text('{:%y/%m/%d %H:%M:%S}'.format(my_time.datetime))
    return points, title

def animate(i):
    points._offsets3d = (eci_pos_x[:, 0, i], eci_pos_y[:, 0, i], eci_pos_z[:, 0, i])
    my_time.mjd = mjds[i]
    title.set_text('{:%y/%m/%d %H:%M:%S}'.format(my_time.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=eci_pos_x.shape[2], interval=20, blit=True
    )

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

In [18]:

FFwriter = animation.FFMpegWriter(
    fps=30, bitrate=8000,
    extra_args=['-vcodec', 'libx264'],
    )
anim.save('starlink_constellation_eci_720p.mp4', writer=FFwriter)

In [19]:

my_time = cysgp4.PyDateTime()
my_time.mjd = mjds[0]
vmin, vmax = np.log10(100), np.log10(10000)

fig = plt.figure(figsize=(12.8, 7.2), dpi=100)
ax1 = fig.add_axes((0.1, 0.5, 0.8, 0.35))
ax2 = fig.add_axes((0.1, 0.1, 0.8, 0.35))
cax = fig.add_axes((0.91, 0.2, 0.02, 0.5))
ax2.set_xlabel('Azimuth [deg]')
ax1.set_ylabel('Elevation [deg]')
for ax in [ax1, ax2]:
    ax.set_xlim((-180, 180))
    ax.set_ylim((0, 90))
    ax.set_xticks(range(-150, 180, 30))
    ax.set_yticks(range(0, 91, 30))
    ax.set_aspect('equal')
    ax.grid()

points1 = ax1.scatter(
    topo_pos_az[:, 1, 0], topo_pos_el[:, 1, 0],
    c=np.log10(topo_pos_dist[:, 1, 0]),
    cmap='viridis', vmin=vmin, vmax=vmax,
    )
points2 = ax2.scatter(
    topo_pos_az[:, 0, 0], topo_pos_el[:, 0, 0],
    c=np.log10(topo_pos_dist[:, 0, 0]),
    cmap='viridis', vmin=vmin, vmax=vmax,
    )
cbar = fig.colorbar(points1, cax=cax)
cbar.set_label('Distance (km)')
cbar.set_ticks([2, 3, 4])
cbar.set_ticklabels([100, 1000, 10000])

ax1.text(-170, 75, 'Parkes 64-m', fontsize=16)
ax2.text(-170, 75, 'Effelsberg 100-m', fontsize=16)
title = ax1.text(
    174, 75, '{:%y/%m/%d %H:%M:%S}'.format(my_time.datetime),
    fontsize=15, ha='right'
    )

def init():
    points1.set_offsets(np.column_stack([topo_pos_az[:, 1, 0], topo_pos_el[:, 1, 0]]))
    points1.set_array(np.log10(topo_pos_dist[:, 1, 0]))
    points2.set_offsets(np.column_stack([topo_pos_az[:, 0, 0], topo_pos_el[:, 0, 0]]))
    points2.set_array(np.log10(topo_pos_dist[:, 0, 0]))
    my_time.mjd = mjds[0]
    title.set_text('{:%y/%m/%d %H:%M:%S}'.format(my_time.datetime))
    return points, title

def animate(i):
    points1.set_offsets(np.column_stack([topo_pos_az[:, 1, i], topo_pos_el[:, 1, i]]))
    points1.set_array(np.log10(topo_pos_dist[:, 1, i]))
    points2.set_offsets(np.column_stack([topo_pos_az[:, 0, i], topo_pos_el[:, 0, i]]))
    points2.set_array(np.log10(topo_pos_dist[:, 0, i]))
    my_time.mjd = mjds[i]
    title.set_text('{:%y/%m/%d %H:%M:%S}'.format(my_time.datetime))
    return points, title

anim = animation.FuncAnimation(
    fig, animate, init_func=init, frames=topo_pos_az.shape[2], interval=20, blit=True
    )

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

In [20]:

FFwriter = animation.FFMpegWriter(
    fps=30, bitrate=8000,
    extra_args=['-vcodec', 'libx264'],
    )
anim.save('starlink_constellation_horizon_720p.mp4', writer=FFwriter)

In [ ]: