In [1]:
!which python; python -V;

# This makes the diagrams to more reliably appear in Jupyter environment
import plotly.io as pio
pio.renderers.default = "notebook_connected"

# This will cause the ephemerides to be imported from JPL horizons system
from astropy.coordinates import solar_system_ephemeris
solar_system_ephemeris.set("jpl")

/home/thomson/devel/perylune/venv/bin/python
Python 3.8.5
Out[1]:
<ScienceState solar_system_ephemeris: 'jpl'>
In [2]:
from poliastro.bodies import Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Sun, Pluto
from poliastro.twobody import Orbit

from poliastro.constants import GM_earth

from astropy import units as u
import numpy as np

from perylune.orbit_tools import *
from perylune.interplanetary import *
from perylune.constants import *
In [3]:
# Let's create two departure orbits. Typically the departure orbit is circular (orb1), but we can make some
# interesting observations on eliptical orbit (orb2) and study how the escape velocity changes between apoapsis and periapsis.
orb1 = Orbit.circular(Earth, 250*u.km)
orb2 = Orbit.from_classical(Earth, 350*u.km + Earth.R, 0.14285714285714285*u.one, 0*u.deg, 0*u.deg, 0*u.deg, 0*u.deg)
In [4]:
print_orb(orb1)
print_orb(orb2)

6628 x 6628 km x 0.0 deg (GCRS) orbit around Earth (♁) at epoch J2000.000 (TT)
a(𝑎)=6628.1366km, b=6628.1366km, e=0.00, i=0.00deg raan(Ω)=0.00deg argp(𝜔)=0.00deg nu(𝜈)=0.00deg
period=5370.30s perapis=6628.1366km(250.00km) apoapsis=6628.1366km(250.00km)
5767 x 7689 km x 0.0 deg (GCRS) orbit around Earth (♁) at epoch J2000.000 (TT)
a(𝑎)=6728.1366km, b=6659.1282km, e=0.14, i=0.00deg raan(Ω)=0.00deg argp(𝜔)=0.00deg nu(𝜈)=0.00deg
period=5492.29s perapis=5766.9742km(-611.16km) apoapsis=7689.2990km(1311.16km)
In [5]:
# Calculate delta-v necessary to reach escape velocity for circular orbit. The values returned are current position, periapsis, apoapis.
# No surprises here - it's perfectly circular, so it's always the same.
escape_vel(orb1, False)
Out[5]:
(<Quantity 10967.00800742 m / s>,
 <Quantity 10967.00800742 m / s>,
 <Quantity 10967.00800742 m / s>)
In [6]:
# The situation for eliptical orbit is more interesting. The delta-v for periapsis burn is 2869 m/s, which is the smallest value
# for this orbit and up to 3516 for apoapsis (when the velocity is smallest, so extra delta-v is needed to escape)
escape_vel(orb2, False)[1]
Out[6]:
$11757.362 \; \mathrm{\frac{m}{s}}$
In [7]:
# To obtain data, head to Horizons (https://ssd.jpl.nasa.gov) and use the following parameters:
# Ephemepris type: ELEMENTS
# Target body: <whatever planet needed, e.g. Mars>
# Center: Sun
# Time span: pick any date
# Table settings: output inits=KM-S
# Display/Output: plain text

# For example, Earth returns the following data for 2020-Nov-03
#2459156.500000000 = A.D. 2020-Nov-03 00:00:00.0000 TDB 
# EC= 1.623243125148644E-02 QR= 1.470544578776156E+08 IN= 4.337517056126222E-03
# OM= 1.956969904857147E+02 W = 2.693660952667993E+02 Tp=  2459219.993764857296
# N = 1.142090613016157E-05 MA= 2.973464932363198E+02 TA= 2.956789660740816E+02
# A = 1.494808962493945E+08 AD= 1.519073346211734E+08 PR= 3.152114165874045E+07

a = 1.494808962493945E+08 * u.km
e = 1.623243125148644E-02 * u.one
inc = 4.337517056126222E-03 * u.deg
raan = 1.956969904857147E+02 * u.deg
argp = 2.693660952667993E+02 * u.deg
nu = 2.956789660740816E+02 * u.deg

# Orbital data from Horizons
orbit_earth = Orbit.from_classical(Sun, a, e, inc, raan, argp, nu)

# MARS
# EC= 9.333873836923647E-02 QR= 2.066557040253017E+08 IN= 1.847950769090515E+00
# OM= 4.949400395133841E+01 W = 2.866239949629700E+02 Tp=  2459064.901374706533
# N = 6.065617422097314E-06 MA= 4.800403158507560E+01 TA= 5.659109019293314E+01
# A = 2.279304441149290E+08 AD= 2.492051842045562E+08 PR= 5.935092422553787E+07
a = 2.279304441149290E+08 * u.km
e = 1.623243125148644E-02 * u.one
inc = 1.847950769090515E+00 * u.deg
raan = 4.949400395133841E+01 * u.deg
argp = 2.866239949629700E+02  * u.deg
nu = 5.659109019293314E+01 * u.deg

# TODO: insert actual date
orbit_mars = Orbit.from_classical(Sun, a, e, inc, raan, argp, nu)

hoh = Maneuver.hohmann(orbit_earth, orbit_mars.a)
orb3,orb4 = orbit_earth.apply_maneuver(hoh, intermediate=True)
In [8]:
from poliastro.plotting import OrbitPlotter3D

op = OrbitPlotter3D()
op.plot(orbit_earth, label="Earth")
op.plot(orbit_mars, label="Mars")
op.plot(orb3, label="transfer")
#op.plot(orb4, label="arrival")
#op.plot(orb3, label="Transfer orbit - Hohmann transfer")
#op.plot(orb4, label="Final orbit - GEO")
In [9]:
# Show cost of the first burn
get_cost(hoh, 0)
Out[9]:
$2838.7196 \; \mathrm{\frac{m}{s}}$
In [10]:
# Show cost of the second burn
get_cost(hoh,1)
Out[10]:
$2759.9828 \; \mathrm{\frac{m}{s}}$
In [11]:
# and now show a summary of the transfer. This is the total cost, assuming the craft enters orbit around Mars.
hoh
Out[11]:
Number of impulses: 2, Total cost: 5.598702 km / s
In [12]:
from astropy import time
from poliastro.ephem import Ephem
from poliastro.util import time_range

# Mars 2020
date_launch = time.Time("2020-07-20 11:50", scale="utc").tdb
date_arrival = time.Time("2021-02-18 12:00", scale="utc").tdb

# Ideal, read from the porkchop
date_launch = time.Time("2020-07-15 11:50", scale="utc").tdb
date_arrival = time.Time("2021-01-20 12:00", scale="utc").tdb

earth = Ephem.from_body(Earth, time_range(date_launch, end=date_arrival))
mars = Ephem.from_body(Mars, time_range(date_launch, end=date_arrival))

# Solve for departure and target orbits
ss_earth = Orbit.from_ephem(Sun, earth, date_launch)
ss_mars = Orbit.from_ephem(Sun, mars, date_arrival)

# Solve for the transfer maneuver
man_lambert = Maneuver.lambert(ss_earth, ss_mars)

# Get the transfer and final orbits
ss_trans, ss_target = ss_earth.apply_maneuver(man_lambert, intermediate=True)

from poliastro.plotting import OrbitPlotter3D
plotter = OrbitPlotter3D()
plotter.set_attractor(Sun)

plotter.plot_ephem(earth, date_launch, label="Earth at launch (2020-07-20)")
plotter.plot_ephem(mars, date_arrival, label="Mars at arrival (2021-02-18)")
plotter.plot_trajectory(ss_trans.sample(max_anomaly=140 * u.deg), color="black", label="Transfer orbit")
plotter.set_view(30 * u.deg, 260 * u.deg, distance=3 * u.km)

print(man_lambert)

Number of impulses: 2, Total cost: 6.884949 km / s
In [13]:
plotter._layout.autosize = False
plotter._layout.width = 1200
plotter._layout.height = 800
plotter._layout.margin=dict(l=10, r=10, b=10, t=10, pad=4 )
plotter.show()
In [14]:
print_orb(ss_trans)

1 x 2 AU x 24.0 deg (HCRS) orbit around Sun (☉) at epoch 2020-07-15 11:51:09.184 (TDB)
a(𝑎)=1.3303AU, b=1.2906AU, e=0.24, i=24.01deg raan(Ω)=356.87deg argp(𝜔)=298.85deg nu(𝜈)=-2.94deg
period=560.45d perapis=1.0078AU apoapsis=1.6529AU
In [ ]: