Examples of utilisation

Ground track

#!/usr/bin/env python

"""Script showing the position of the ISS at the time of the TLE
and the ground track for the previous and the next orbit
"""

import sys
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

from beyond.io.tle import Tle
from beyond.dates import Date, timedelta


# Parsing of TLE
tle = Tle("""ISS (ZARYA)
1 25544U 98067A   19004.59354167  .00000715  00000-0  18267-4 0  9995
2 25544  51.6416  95.0104 0002419 236.2184 323.8248 15.53730729149833""")

# Conversion into `Orbit` object
orb = tle.orbit()

# Tables containing the positions of the ground track
latitudes, longitudes = [], []
prev_lon, prev_lat = None, None

period = orb.infos.period
start = orb.date - period
stop = 2 * period
step = period / 100

for point in orb.ephemeris(start=start, stop=stop, step=step):

    # Conversion to earth rotating frame
    point.frame = 'ITRF'

    # Conversion from cartesian to spherical coordinates (range, latitude, longitude)
    point.form = 'spherical'

    # Conversion from radians to degrees
    lon, lat = np.degrees(point[1:3])

    # Creation of multiple segments in order to not have a ground track
    # doing impossible paths
    if prev_lon is None:
        lons = []
        lats = []
        longitudes.append(lons)
        latitudes.append(lats)
    elif orb.i < np.pi /2 and (np.sign(prev_lon) == 1 and np.sign(lon) == -1):
        lons.append(lon + 360)
        lats.append(lat)
        lons = [prev_lon - 360]
        lats = [prev_lat]
        longitudes.append(lons)
        latitudes.append(lats)
    elif orb.i > np.pi/2 and (np.sign(prev_lon) == -1 and np.sign(lon) == 1):
        lons.append(lon - 360)
        lats.append(lat)
        lons = [prev_lon + 360]
        lats = [prev_lat]
        longitudes.append(lons)
        latitudes.append(lats)

    lons.append(lon)
    lats.append(lat)
    prev_lon = lon
    prev_lat = lat

img = Path(__file__).parent / "earth.png"

im = plt.imread(str(img))
plt.figure(figsize=(15.2, 8.2))
plt.imshow(im, extent=[-180, 180, -90, 90])

for lons, lats in zip(longitudes, latitudes):
    plt.plot(lons, lats, 'r')

lon, lat = np.degrees(orb.copy(frame='ITRF', form='spherical')[1:3])
plt.plot([lon], [lat], 'ro')

plt.xlim([-180, 180])
plt.ylim([-90, 90])
plt.grid(True, color='w', linestyle=":", alpha=0.4)
plt.xticks(range(-180, 181, 30))
plt.yticks(range(-90, 91, 30))
plt.tight_layout()

if "no-display" not in sys.argv:
    plt.show()
_images/ground-track-ISS.png

Station pointings

#!/usr/bin/env python

import sys
import numpy as np
import matplotlib.pyplot as plt

from beyond.dates import Date, timedelta
from beyond.io.tle import Tle
from beyond.frames import create_station
from beyond.config import config


tle = Tle("""ISS (ZARYA)
1 25544U 98067A   16086.49419020  .00003976  00000-0  66962-4 0  9998
2 25544  51.6423 110.4590 0001967   0.7896 153.8407 15.54256299992114""").orbit()

# Station definition
station = create_station('TLS', (43.428889, 1.497778, 178.0))
azims, elevs = [], []

print("    Time      Azim    Elev    Distance   Radial Velocity")
print("=========================================================")

for orb in station.visibility(tle, start=Date.now(), stop=timedelta(hours=24), step=timedelta(seconds=30), events=True):
    elev = np.degrees(orb.phi)
    # Radians are counterclockwise and azimuth is clockwise
    azim = np.degrees(-orb.theta) % 360

    # Archive for plotting
    azims.append(azim)
    # Matplotlib actually force 0 to be at the center of the polar plot,
    # so we trick it by inverting the values
    elevs.append(90 - elev)

    r = orb.r / 1000.
    print("{event:7} {orb.date:%H:%M:%S} {azim:7.2f} {elev:7.2f} {r:10.2f} {orb.r_dot:10.2f}".format(
        orb=orb, r=r, azim=azim, elev=elev, event=orb.event.info if orb.event is not None else ""
    ))

    if orb.event and orb.event.info.startswith("LOS"):
        # We stop at the end of the first pass
        print()
        break

plt.figure()
ax = plt.subplot(111, projection='polar')
ax.set_theta_direction(-1)
ax.set_theta_zero_location('N')
plt.plot(np.radians(azims), elevs, '.')
ax.set_yticks(range(0, 90, 20))
ax.set_yticklabels(map(str, range(90, 0, -20)))
ax.set_rmax(90)

if "no-display" not in sys.argv:
    plt.show()
    Time      Azim    Elev    Distance   Radial Velocity
=========================================================
AOS 02:01:23  302.10   -0.00    2312.91   -6906.48
    02:01:39  302.30    1.00    2204.46   -6904.71
    02:02:09  302.71    3.07    1997.46   -6894.10
    02:02:39  303.20    5.39    1790.94   -6870.86
    02:03:09  303.79    8.08    1585.39   -6829.11
    02:03:39  304.54   11.30    1381.49   -6758.37
    02:04:09  305.53   15.29    1180.37   -6638.85
    02:04:39  306.94   20.51     984.03   -6430.17
    02:05:09  309.17   27.78     796.31   -6042.60
    02:05:39  313.28   38.63     625.31   -5262.66
    02:06:09  323.59   55.36     489.17   -3610.39
    02:06:39   12.81   74.28     423.23    -565.54
MAX 02:06:43   30.39   75.00     421.90       0.00
    02:07:09   90.23   61.99     458.89    2806.46
    02:07:39  105.74   43.24     577.62    4871.54
    02:08:09  111.00   30.79     740.40    5857.82
    02:08:39  113.63   22.61     924.12    6336.73
    02:09:09  115.23   16.86    1118.38    6588.72
    02:09:39  116.32   12.54    1318.37    6730.93
    02:10:09  117.12    9.12    1521.66    6814.75
    02:10:39  117.74    6.29    1726.92    6864.71
    02:11:09  118.24    3.86    1933.33    6893.49
    02:11:39  118.65    1.72    2140.39    6908.11
LOS 02:12:05  118.97    0.01    2323.42    6912.59
_images/station.png

Jupiter and its moons

"""Follow the informations provided in the `beyond.env.jpl` module about configuration
in order to supply the data needed.
"""

import numpy as np
import matplotlib.pyplot as plt

from beyond.dates import Date
from beyond.env.jpl import get_orbit
from beyond.frames import create_station
from beyond.config import config

# Load the ".bsp" file
config.update({
    "env": {
        "jpl": [
            "/home/jules/.space/jpl/jup310.bsp"
        ]
    }
})


date = Date.now()

# Definition of the location of observation
station = create_station('TLS', (43.428889, 1.497778, 178.0))

# Retrieve Jupiter and its moons state-vectors
jupiter = get_orbit('Jupiter', date)
io = get_orbit('Io', date)
europa = get_orbit('Europa', date)
ganymede = get_orbit('Ganymede', date)
callisto = get_orbit('Callisto', date)

# Convert them to the observer frame
jupiter.frame = station
io.frame = station
europa.frame = station
ganymede.frame = station
callisto.frame = station

# Change the orbit form in order to get azimuth and elevation
europa.form = 'spherical'
io.form = 'spherical'
ganymede.form = 'spherical'
callisto.form = 'spherical'
jupiter.form = 'spherical'

plt.figure()

# The size of the dots is not on scale
plt.plot([-np.degrees(io.theta)], [np.degrees(io.phi)], 'o', zorder=-io.r, label="Io")
plt.plot([-np.degrees(europa.theta)], [np.degrees(europa.phi)], 'o', zorder=-europa.r, label="Europa")
plt.plot([-np.degrees(ganymede.theta)], [np.degrees(ganymede.phi)], 'o', zorder=-ganymede.r, label="Ganymede")
plt.plot([-np.degrees(callisto.theta)], [np.degrees(callisto.phi)], 'o', zorder=-callisto.r, label="Callisto")
plt.plot([-np.degrees(jupiter.theta)], [np.degrees(jupiter.phi)], 'o', zorder=-jupiter.r, ms=12, label="Jupiter")

plt.suptitle("Jupiter as seen from Toulouse")
plt.title(f"{date:%Y-%m-%d %H:%M:%S UTC}")
plt.legend()
plt.axis('equal')
plt.show()
_images/jupiter.png

Listeners

#!/usr/bin/env python

from beyond.dates import Date, timedelta
from beyond.io.tle import Tle
from beyond.frames import create_station
from beyond.propagators.listeners import stations_listeners, NodeListener, ApsideListener, LightListener


tle = Tle("""ISS (ZARYA)
1 25544U 98067A   17153.89608442  .00001425  00000-0  28940-4 0  9997
2 25544  51.6419 109.5559 0004850 223.1783 180.8272 15.53969766 59532""").orbit()

# Station definition
station = create_station('TLS', (43.428889, 1.497778, 178.0))

# Listeners declaration
listeners = stations_listeners(station)  # AOS, LOS and MAX elevation
listeners.append(NodeListener())         # Ascending and Descending Node
listeners.append(ApsideListener())       # Apogee and Perigee
listeners.append(LightListener())        # Illumination events

start = Date.now()
stop = timedelta(minutes=100)
step = timedelta(seconds=180)

for orb in tle.iter(start=start, stop=stop, step=step, listeners=listeners):
    event = orb.event if orb.event is not None else ""
    print(f"{orb.date:%Y-%m-%d %H:%M:%S} {event}")
2017-06-03 21:39:11
2017-06-03 21:42:11
2017-06-03 21:45:11
2017-06-03 21:48:11
2017-06-03 21:51:11
2017-06-03 21:52:37 Apogee
2017-06-03 21:54:11
2017-06-03 21:57:11
2017-06-03 22:00:11
2017-06-03 22:00:43 AscNode
2017-06-03 22:03:11
2017-06-03 22:06:11
2017-06-03 22:09:11
2017-06-03 22:12:11
2017-06-03 22:15:11
2017-06-03 22:18:11
2017-06-03 22:21:11
2017-06-03 22:24:11
2017-06-03 22:25:14 AOS TLS
2017-06-03 22:27:11
2017-06-03 22:28:16 Perigee
2017-06-03 22:30:11
2017-06-03 22:30:25 Umbra in
2017-06-03 22:30:36 MAX TLS
2017-06-03 22:33:11
2017-06-03 22:35:57 LOS TLS
2017-06-03 22:36:11
2017-06-03 22:39:11
2017-06-03 22:42:11
2017-06-03 22:45:11
2017-06-03 22:46:59 DscNode
2017-06-03 22:48:11
2017-06-03 22:51:11
2017-06-03 22:54:11
2017-06-03 22:57:11
2017-06-03 23:00:11
2017-06-03 23:03:11
2017-06-03 23:03:18 Umbra out
2017-06-03 23:06:11
2017-06-03 23:09:11
2017-06-03 23:12:11
2017-06-03 23:15:11
2017-06-03 23:18:11

Maneuvers


"""Example of Hohmann transfer

The orbit we are starting with is a Tle of the ISS. The amplitude of the maneuver is greatly
exagerated regarding the ISS's capability, but has the convenience to be particularly visual.
"""

import sys
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from beyond.io.tle import Tle
from beyond.dates import timedelta
from beyond.propagators.keplernum import KeplerNum
from beyond.env.solarsystem import get_body
from beyond.orbits.man import ImpulsiveMan
from beyond.propagators.listeners import ApsideListener, find_event


orb = Tle("""ISS (ZARYA)
1 25544U 98067A   18124.55610684  .00001524  00000-0  30197-4 0  9997
2 25544  51.6421 236.2139 0003381  47.8509  47.6767 15.54198229111731""").orbit()

start = orb.date
stop = timedelta(minutes=300)
step = timedelta(seconds=60)

# Changing the propagator to Keplerian, as SGP4 is not able to perform maneuvers
orb.propagator = KeplerNum(step, bodies=get_body("Earth"))

# Research for the next perigee
perigee = find_event(orb.iter(stop=stop, listeners=ApsideListener()), 'Periapsis')

man1 = ImpulsiveMan(perigee.date, (280, 0, 0), frame="TNW")
orb.maneuvers = [man1]

dates1, alt1 = [], []

# Research for the next apogee after the first maneuver
apogee = find_event(orb.iter(start=perigee.date - step * 10, stop=stop, listeners=ApsideListener()), 'Apoapsis')
# apogee = find_event(orb.iter(stop=stop, listeners=ApsideListener()), 'Apoapsis', offset=1)

# Adding the second maneuver to the orbit
man2 = ImpulsiveMan(apogee.date, (270, 0, 0), frame="TNW")
orb.maneuvers.append(man2)

print(man1.date)
print(man2.date)

# Propagation throught the two maneuvers
ephem = orb.ephem(start=start, stop=stop, step=step)

# graphs
plt.figure()

data = np.array(ephem)
dates = [x.date for x in ephem]
# Altitude in km
alt = (np.linalg.norm(data[:, :3], axis=1) - orb.frame.center.body.r) / 1000
events_dates = [perigee.date, apogee.date]
events_alt = (np.linalg.norm([perigee[:3], apogee[:3]], axis=1) - orb.frame.center.body.r) / 1000

plt.plot(dates, alt)
plt.plot([events_dates[0]], [events_alt[0]], 'ro', label="perigee")
plt.plot([events_dates[1]], [events_alt[1]], 'ko', label="apogee")

plt.ylabel("altitude (km)")
plt.legend()
plt.grid(linestyle=':', alpha=0.4)
plt.tight_layout()

fig = plt.figure()
ax = plt.gca(projection='3d')
ax.view_init(elev=52, azim=140)

x, y, z = zip(perigee[:3], apogee[:3])

plt.plot(data[:, 0], data[:, 1], data[:, 2])
plt.plot([perigee[0]], [perigee[1]], [perigee[2]], 'ro')
plt.plot([apogee[0]], [apogee[1]], [apogee[2]], 'ko')

if "no-display" not in sys.argv:
    plt.show()
2018-05-04T14:41:09.984493 UTC
2018-05-04T15:33:11.374985 UTC
_images/hohmann_1.png _images/hohmann_2.png

Docking

"""Example script for docking to the ISS, starting from an orbit 3000 m below

This script uses the V-bar approach strategy, with an approach ellipse at 2000
meters and a keep-out sphere at 200 meters.
"""

import sys
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import matplotlib.patches as mpatches

from beyond.dates import Date, timedelta
from beyond.orbits import Orbit
from beyond.propagators.cw import ClohessyWiltshire
from beyond.utils.cwhelper import CWHelper
from beyond.propagators.listeners import Listener, Event, find_event
from beyond.orbits.man import ImpulsiveMan, ContinuousMan


class TangentialDistanceListener(Listener):
    """This listener is useful to detect the exact time at which to
    start the Hohmann transfer in order to place the chaser spacecraft
    at the first hold point
    """

    def __init__(self, delta, orientation="QSW"):
        self.delta = delta
        self.orientation = orientation

    def info(self, orb):
        return Event(self, "go")

    def __call__(self, orb):
        idx = 1 if self.orientation == "QSW" else 0
        return orb[idx] - self.delta


propagator = ClohessyWiltshire(6378000 + 410000)
helper = CWHelper(propagator)

radial = -3000  # radial distance at the beginning of the simulation
date = Date(Date.now().d)
step = timedelta(seconds=10)
hold_points = [
    -2500,  # First hold point well outside the approach ellipse
    -200,  # Second hold point just outside the keep-out sphere
    -10,  # Last hold point close to the docking port for the last checks
]
hold = timedelta(minutes=10)  # hold duration
linear_dv = 0.5  # Proximity linear velocity
final_dv = 0.1  # Final linear velocity for docking
duration = helper.period * 3

# Tangential distance necessary for a Hohmann transfer
tangential = helper.hohmann_distance(radial)

# Initial orbit
orb = helper.coelliptic(date, radial, tangential * 1.5)

# Search for the point of desired maneuver in order to arrive at the
# first hold point
ephem = orb.ephem(
    step=step,
    stop=helper.period * 2,
    listeners=TangentialDistanceListener(
        tangential + hold_points[0], helper.frame.orientation
    ),
)
start_man = find_event(ephem, "go")

# Creation of the two impulsive maneuvers that place the spacecraft at the
# first hold point (at 2500 m)
orb.maneuvers.extend(helper.hohmann(-radial, start_man.date))

# Eccentric boost in order to go to the next hold point (at 200 m)
orb.maneuvers.extend(
    helper.eccentric_boost(
        hold_points[1] - hold_points[0], orb.maneuvers[-1].date + hold
    )
)

# Linear maneuver at 0.5 m/s to reach the next hold point (at 10 m)
orb.maneuvers.extend(
    helper.vbar_linear(
        hold_points[2] - hold_points[1], orb.maneuvers[-1].date + hold, linear_dv
    )
)

# Final maneuver at 0.1 m/s to reach the docking port
orb.maneuvers.extend(
    helper.vbar_linear(-hold_points[2], orb.maneuvers[-1].date + hold, final_dv)
)

chaser_ephem = orb.ephem(step=step, stop=duration)
dates = list(chaser_ephem.dates)
result = np.asarray(chaser_ephem)

# Display the trajectory
print("Final position")
print(f"Q = {chaser_ephem[-1, 0]: 7.3f}")
print(f"S = {chaser_ephem[-1, 1]: 7.3f}")
print(f"W = {chaser_ephem[-1, 2]: 7.3f}")

man_dates = [man.date for man in orb.maneuvers if isinstance(man, ImpulsiveMan)]
cont_man = [man for man in orb.maneuvers if isinstance(man, ContinuousMan)]
mans = np.array([chaser_ephem.propagate(d) for d in man_dates])

plt.figure(figsize=(12.8, 4.8))
ax = plt.subplot(111)

# Approach sphere (2 km)
circle = mpatches.Ellipse((0, 0), 4000, 2000, ec="none", fc="#aaaaaa60")
ax.add_patch(circle)

# Keep out sphere (200 m) with corridors
w1 = mpatches.Wedge((0, 0), 200, theta1=-170, theta2=-100, ec="none", fc="#aaaaaa")
w2 = mpatches.Wedge((0, 0), 200, theta1=-80, theta2=-10, ec="none", fc="#aaaaaa")
w3 = mpatches.Wedge((0, 0), 200, theta1=10, theta2=170, ec="none", fc="#aaaaaa")
ax.add_patch(w1)
ax.add_patch(w2)
ax.add_patch(w3)

plt.text(0, -750, "Approach ellipse", ha="center")
plt.text(0, 200, "Keep-out sphere", ha="center")

plt.plot(chaser_ephem[:, 1], chaser_ephem[:, 0])
if len(mans):
    plt.plot(mans[:, 1], mans[:, 0], "o")

for man in cont_man:
    subeph = chaser_ephem.ephem(start=man.start, stop=man.stop, step=step)
    plt.plot(subeph[:, 1], subeph[:, 0], "r")

plt.xlabel("S (m)")
plt.ylabel("Q (m)")
ax.set_aspect("equal")
plt.grid(ls=":")
plt.xlim(-11000, 1000)
plt.ylim(-3200, 500)
plt.tight_layout()

if "no-display" not in sys.argv:
    plt.show()
Final position
Q =   0.000
S =   0.000
W =   0.000
_images/docking.png