In [1]:
!which python; python -V;
from astropy import units as u
import numpy as np

from poliastro.bodies import Earth, Mars, Sun
from poliastro.twobody import Orbit

# Needed to handle TLE into Poliastro's Orbit
from tletools import TLE

import plotly.io as pio
pio.renderers.default = "notebook_connected"

# Needed for defining manouvers
from poliastro.maneuver import Maneuver
from poliastro.twobody import thrust
from poliastro.twobody.propagation import cowell
from poliastro.plotting import OrbitPlotter2D, OrbitPlotter3D

from perylune.orbit_tools import *

/home/thomson/devel/perylune/venv/bin/python
Python 3.8.5
In [2]:
# STEP 0: DEFINE TARGET ORBIT
# Let's get the TLE orbital data for the target satellite. As an example, let's pick a dead NOAA-17 sat.
# Its TLE data can be obtained from many places, such as celestrak or n2yo (https://www.n2yo.com/satellite/?s=27453)
tle_text = """NOAA-17
1 27453U 02032A   20263.80942421 -.00000011 +00000-0 +13583-4 0  9998
2 27453 098.5909 208.3215 0011096 327.5463 032.5033 14.25072668948324"""
tle_lines = tle_text.strip().splitlines()
tle_target = TLE.from_lines(*tle_lines)
orb_target = tle_target.to_orbit()
In [3]:
# STEP 1: LAUNCH
# Let's use ANDESITE sat as an example. It was launched as shareride by RocketLab's Electron on
# the latest rideshare launch on 2020 June 13. NORAD ID is 45726. 
# TLE elements taken from n2yo website, details https://www.n2yo.com/satellite/?s=45726

tle_text = """GDASAT-1
1 45726U 20037D   20278.45278018  .00000608  00000-0  65390-4 0  9991
2 45726  97.7132  96.2906 0012962 283.6573  76.3213 14.92011802 14275"""
tle_lines = tle_text.strip().splitlines()
tle1 = TLE.from_lines(*tle_lines)
orb1 = tle1.to_orbit()
In [4]:
# Let's get the details about initial orbit
print_orb(orb1)

6961 x 6979 km x 97.7 deg (GCRS) orbit around Earth (♁) at epoch 2020-10-04T10:52:00.207552000 (UTC)
a(𝑎)=6969.80km, b=6969.79km, e=0.00, i=97.71deg raan(Ω)=96.29deg argp(𝜔)=283.66deg nu(𝜈)=76.47deg
period=5790.84s perapis=6961km(582.63km) apoapsis=6979km(600.70km)
In [5]:
# Let's take a look at the target orbit
print_orb(orb_target)

7178 x 7194 km x 98.6 deg (GCRS) orbit around Earth (♁) at epoch 2020-09-19T19:25:34.251744000 (UTC)
a(𝑎)=7186.39km, b=7186.38km, e=0.00, i=98.59deg raan(Ω)=208.32deg argp(𝜔)=327.55deg nu(𝜈)=32.57deg
period=6062.85s perapis=7178km(800.27km) apoapsis=7194km(816.22km)
In [6]:
# Let's visualize the orbits

# Obtain the orbit after applying Hohmann transfer.
plot = OrbitPlotter3D(dark=False)
#plot._layout['autosize'] = False
#print(plot._layout['template'])
plot.plot(orb1, color="orange", label="GDASAT-1, after orbital insertion")
plot.plot(orb_target, color="red", label="NOAA-17, target sat")
#plot._layout.width=1200
#plot._layout.height=800
plot.show();
In [7]:
# STEP 2: Synchronize Right Ascension of the Ascending Node. Changing RAAN as any other out of plane maneuver
# is in general a costly operation. However, there's a great saving trick that allows performing
# such change almost for free - the J2 perturbation (oblate earth)

# First, lets estimate how fast the orbit precesses in the course of 24 hours.

from poliastro.core.perturbations import J2_perturbation
tof = (24.0 * u.h).to(u.s)
orb24h = orb1.propagate(tof, method=cowell, ad=J2_perturbation, J2=Earth.J2.value, R=Earth.R.to(u.km).value)
print("Orbit precessed from RAAN=%s to %s in 24h, giving the rate of %s/24h." % (orb1.raan, orb24h.raan.to(u.deg), (orb24h.raan - orb1.raan).to(u.deg)))

# The target orbit has RAAN of 208.32deg. Our initial orbit has RANN of 96.29deg. The orbital plane needs to rotate by 112.03deg.
# Since the perfromance is 0.988 deg/24h, the rotation will take roughly 113 days. Let's do the 113 days in one simulation...
tof = (113 * 24 * u.h).to(u.s)
orb2 = orb1.propagate(tof, method=cowell, ad=J2_perturbation, J2=Earth.J2.value, R=Earth.R.to(u.km).value)

print("RAAN after 113 days is %s" % orb2.raan.to(u.deg))

# ... and then continue simulating one hour at a time until the RAAN is at least the same as that of target's orbit:
cnt = 0
while orb2.raan < orb_target.raan:
    tof = (10*u.min).to(u.s)
    cnt +=1
    orb2 = orb2.propagate(tof, method=cowell, ad=J2_perturbation, J2=Earth.J2.value, R=Earth.R.to(u.km).value)

print("Orbit after J2 perturbation (%d days, %d mins):" % (113, cnt*10))
print_orb(orb2)

print("Target orbit:")
print_orb(orb_target)

Orbit precessed from RAAN=96.2906 deg to 97.27942040624983 deg in 24h, giving the rate of 0.9888204062498249 deg/24h.
RAAN after 113 days is 207.57302661216775 deg
Orbit after J2 perturbation (113 days, 1100 mins):
6953 x 6967 km x 97.7 deg (GCRS) orbit around Earth (♁) at epoch 2021-01-26T05:12:00.207552000 (UTC)
a(𝑎)=6960.24km, b=6960.24km, e=0.00, i=97.72deg raan(Ω)=208.32deg argp(𝜔)=282.47deg nu(𝜈)=-56.73deg
period=5778.93s perapis=6953km(575.25km) apoapsis=6967km(588.96km)
Target orbit:
7178 x 7194 km x 98.6 deg (GCRS) orbit around Earth (♁) at epoch 2020-09-19T19:25:34.251744000 (UTC)
a(𝑎)=7186.39km, b=7186.38km, e=0.00, i=98.59deg raan(Ω)=208.32deg argp(𝜔)=327.55deg nu(𝜈)=32.57deg
period=6062.85s perapis=7178km(800.27km) apoapsis=7194km(816.22km)
In [8]:
# Let's visualize the orbits. The virtually identical RAAN value makes the orbits almost co-planar. However,
# the only difference is their slightly different inclination: 98.59deg vs 97.71 deg.
plot = OrbitPlotter3D()
plot.plot(orb2, color="orange", label="GDASAT-1, after RAAN drift")
plot.plot(orb_target, color="red", label="NOAA-17, target sat")
plot.show()
In [9]:
# STEP 3: Change orbital inclination.
# Our orbit has inclination of 97.7 deg, while the target is 98.6 deg. Let's calculate the delta inclination
delta_theta = orb_target.inc - orb1.inc
print("Need to change inclination by %s" % delta_theta)

# Now let's wait till the spacecraft flies to the descending node.
orb3 = propagate_to_desc_node(orb2)

man1 = plane_change_maneuver(orb2, delta_theta)
print("Plane change maneuver: %s" % man1)
print("maneuver 1 is %s" % repr(man1[0]))
orb4 = orb3.apply_maneuver(man1)

print("== Current orbit (after inclination change) ==")
print_orb(orb4)

print("== Target orbit ==")
print_orb(orb_target)

# The orbits are now fully co-planar!
# We can use Hohmann transfer now!

Need to change inclination by 0.8777000000000044 deg
Plane change maneuver: Number of impulses: 1, Total cost: 0.115987 km / s
maneuver 1 is (<Quantity 0. s>, <Quantity [ 0.05504824, -0.1007999 ,  0.01619127] km / s>)
== Current orbit (after inclination change) ==
6955 x 6967 km x 98.6 deg (GCRS) orbit around Earth (♁) at epoch 2021-01-26T06:36:05.106362871 (UTC)
a(𝑎)=6960.74km, b=6960.73km, e=0.00, i=98.60deg raan(Ω)=208.32deg argp(𝜔)=279.16deg nu(𝜈)=-99.16deg
period=5779.54s perapis=6955km(576.41km) apoapsis=6967km(588.79km)
== Target orbit ==
7178 x 7194 km x 98.6 deg (GCRS) orbit around Earth (♁) at epoch 2020-09-19T19:25:34.251744000 (UTC)
a(𝑎)=7186.39km, b=7186.38km, e=0.00, i=98.59deg raan(Ω)=208.32deg argp(𝜔)=327.55deg nu(𝜈)=32.57deg
period=6062.85s perapis=7178km(800.27km) apoapsis=7194km(816.22km)
In [10]:
plot = OrbitPlotter2D()
plot.plot(orb4, color="orange", label="GDASAT-1, after inclination change")
plot.plot(orb_target, color="red", label="NOAA-17, target sat")
plot._layout.width=1200
plot._layout.height=800
plot.show()
In [11]:
# STEP 4: Hohmann transfer
# This is generally simple maneuver, however the classic Hohmann maneuver is between two circular orbits.
# Nevertheless we can reasonably use it as the orbits are nearly circular. Hohmann manevuer consists
# of two burns.
man23 = Maneuver.hohmann(orb4, 7178 * u.km)
print(man23)
print("Maneuver 2: %s" % repr(man23[0]))
print("Maneuver 3: %s" % repr(man23[1]))
print("== Current orbit (after Hohmann transfer) ==")
orb5 = orb4.apply_maneuver(man23)
print_orb(orb5)

Number of impulses: 2, Total cost: 0.115398 km / s
Maneuver 2: (<Quantity 1590.30720989 s>, <Quantity [-49.52185833, -25.16882297,   8.85292832] m / s>)
Maneuver 3: (<Quantity 2955.74292612 s>, <Quantity [52.07033882, 26.46405413, -9.30851532] m / s>)
== Current orbit (after Hohmann transfer) ==
7178 x 7178 km x 98.6 deg (GCRS) orbit around Earth (♁) at epoch 2021-01-26T07:51:51.156498878 (UTC)
a(𝑎)=7178.00km, b=7178.00km, e=0.00, i=98.60deg raan(Ω)=208.32deg argp(𝜔)=0.00deg nu(𝜈)=99.16deg
period=6052.24s perapis=7178km(799.86km) apoapsis=7178km(799.86km)
In [12]:
# STEP 5: Final corrections.
# The final touches needed to synchronize orbits. Our orbit is 800x800km, but the target is 800x816km.
# The elliptic burn can be done at the periapsis (which is equal to our altitude), so the raised 
# point (apoapsis) will be higher.

# First we need to calculate the difference of orbital velocity of 800x800km orbit and
# the 800x816km orbit at apoapsis.
from poliastro.constants import GM_earth

r_i = 800000 *u.m + Earth.R
v800 = np.sqrt(GM_earth*(1/r_i))

print("Orbital velocity of 800x800km orbit is %s" % v800)

r_1 = 800000 *u.m + Earth.R
r_2 = 816000 *u.m + Earth.R

v816 = np.sqrt(GM_earth*(2/r_1 - 2/(r_1 + r_2)))

print("Orbital velocity of 800x816km orbit at periapsis is %s" % v816)

# This was first experimentally chosen value is 4.3 m/s. Later on, I've calculated
# the value (see astro-dynamics.ods, Hohmann sheet, M11), also the v800, v816 calculations above)
dv = v816 - v800

print("Required delta-v is %s" % dv)

## TODO: fix this to conduct the burn at a time that will sync argp
orb6 = propagate_to_periapsis(orb5)

man4 = prograde_maneuver(orb6, dv.value, None)
print("Maneuver 4 is %s" % repr(man4[0]))

orb7 = orb6.apply_maneuver(man4)

Orbital velocity of 800x800km orbit is 7451.831541112051 m / s
Orbital velocity of 800x816km orbit at periapsis is 7455.978280403477 m / s
Required delta-v is 4.146739291425547 m / s
Maneuver 4 is (<Quantity 0. s>, <Quantity [-0.00029408,  0.00054566,  0.00410015] km / s>)
In [13]:
# SUMMARY: Orbits are now synced. There are minor differences, but the delta-v required to
# correct it is nagligible (e.g. the cost of raising apoapsis was in the order of 4m/s)
# 
print("== Here's our orbit ==")
print_orb(orb7)

print("== Here's the target orbit ==")
print_orb(orb_target)

== Here's our orbit ==
7178 x 7194 km x 98.6 deg (GCRS) orbit around Earth (♁) at epoch 2021-01-26T09:04:56.363403879 (UTC)
a(𝑎)=7186.00km, b=7186.00km, e=0.00, i=98.60deg raan(Ω)=208.32deg argp(𝜔)=360.00deg nu(𝜈)=0.00deg
period=6062.36s perapis=7178km(799.86km) apoapsis=7194km(815.86km)
== Here's the target orbit ==
7178 x 7194 km x 98.6 deg (GCRS) orbit around Earth (♁) at epoch 2020-09-19T19:25:34.251744000 (UTC)
a(𝑎)=7186.39km, b=7186.38km, e=0.00, i=98.59deg raan(Ω)=208.32deg argp(𝜔)=327.55deg nu(𝜈)=32.57deg
period=6062.85s perapis=7178km(800.27km) apoapsis=7194km(816.22km)
In [14]:
# The total delta-v of all maneuvers combined is as follows
delta_v = man1.get_total_cost() + man23.get_total_cost() + man4.get_total_cost()
delta_v.to(u.m/u.s)
Out[14]:
$235.53247 \; \mathrm{\frac{m}{s}}$
In [15]:
from poliastro.util import norm
print(man1.get_total_cost().to(u.m/u.s))
print(norm(man23._dvs[0]))
print(norm(man23._dvs[1]))
print(man4.get_total_cost().to(u.m/u.s))

115.98743740537127 m / s
56.25174167464684 m / s
59.14655360541718 m / s
4.146739291425546 m / s
In [ ]: