Slide 1

Slide 1 text

poliastro: an Astrodynamics library written in Python poliastro: an Astrodynamics library written in Python poli Astrodynamics in Python astro Juan Luis Cano - 2017-11-23 OSCW @ ESOC Juan Luis Cano - 2017-11-23 OSCW @ ESOC

Slide 2

Slide 2 text

1. Introduction (2 minutes) 1. Introduction (2 minutes) What is poliastro? What is poliastro? Pure Python library for Astrodynamics Orbital element conversion, analytical and numerical propagation, interface to NEOs databases, orbital maneuvers Physical unit handling, astronomical time scales Cool documentation Version 0.8.0 released few days ago! MIT License (permissive, i.e. commercial-friendly) http://docs.poliastro.space/ (http://docs.poliastro.space/) http://docs.poliastro.space/en/v0.8.0 /changelog.html#poliastro-0-8-0-2017-11-18 (http://docs.poliastro.space /en/v0.8.0/changelog.html#poliastro-0-8-0-2017-11-18)

Slide 3

Slide 3 text

Brief history Brief history 2013: First version: Octave + FORTRAN + Python 2014: Refactor of the API, much friendlier 2015: Replace FORTRAN algorithms by Python + numba 2016: Izzo algorithm (Lambert's problem), presentation at 6th ICATT 2017: SOCIS (ESA) grant, OpenAstronomy & Astropy membership

Slide 4

Slide 4 text

Why another... How does it compare to...? How does it compare to...? Orekit: not restricted to Earth orbits, more general NASA SPICE: more lightweight, focused on interactive use, not backed by NASA

Slide 5

Slide 5 text

2. Main features (7½ minutes) 2. Main features (7½ minutes) 2.1. First contact 2.1. First contact Interactive usage on Jupyter notebook Physical units, astronomical time scales and reference frame handling (thanks to Astropy) Orbital elements conversion (cartesian, keplerian, equinoctial) 2D and 3D visualization

Slide 6

Slide 6 text

In [2]: from poliastro.bodies import Earth from poliastro.twobody import Orbit In [3]: r = [-6045, -3490, 2500] * u.km v = [-3.457, 6.618, 2.533] * u.km / u.s ss = Orbit.from_vectors(Earth, r, v, Time.now()) ss Out[3]: 7283 x 10293 km x 153.2 deg orbit around Earth (♁) In [4]: ss.epoch Out[4]:

Slide 7

Slide 7 text

In [6]: from poliastro.plotting import plot, plot3d plot(ss, label="Sample orbit");

Slide 8

Slide 8 text

In [9]: plot3d(ss, label="Sample orbit"); Out[9]:

Slide 9

Slide 9 text

2.2. External data 2.2. External data Planetary ephemerides SPICE kernels Near Earth Objects (NEOs) NeoWs DASTCOM5 database ftp://ssd.jpl.nasa.gov/pub/ssd/README.txt https://api.nasa.gov/neo/?api_key=DEMO_KEY (https://api.nasa.gov/neo/?api_key=DEMO_KEY) In [10]: from astropy.coordinates import solar_system_ephemeris solar_system_ephemeris.set("jpl") Orbit.from_body_ephem(Earth) Out[10]: 1 x 1 AU x 23.4 deg orbit around Sun (☉) In [11]: from poliastro.neos import neows florence = neows.orbit_from_name("Florence") florence Out[11]: 1 x 3 AU x 22.2 deg orbit around Sun (☉) In [12]: from poliastro.neos import dastcom5 halley_1835 = dastcom5.orbit_from_name('1P')[-3] halley_1835 Out[12]: 1 x 35 AU x 162.3 deg orbit around Sun (☉)

Slide 10

Slide 10 text

In [13]: from poliastro.plotting import OrbitPlotter from poliastro.bodies import Mercury, Venus, Earth, Mars frame = OrbitPlotter(num_points=300) # Florence close approach frame.plot(Orbit.from_body_ephem(Earth, epoch=Time("2017-09-01 12:05:50", scale="td b")), label=Earth) for planet in Mercury, Venus, Mars: frame.plot(Orbit.from_body_ephem(planet), label=planet) frame.plot(halley_1835, label='Halley', color='#666666') frame.plot(florence, label='Florence', color='#000000') plt.title("Inner solar system + Florence + Halley") plt.xlim(-.3e9, .3e9) plt.ylim(-.3e9, .3e9) Out[13]: (-300000000.0, 300000000.0)

Slide 11

Slide 11 text

No content

Slide 12

Slide 12 text

2.3. Core algorithms 2.3. Core algorithms Analytical propagation, aka Kepler's problem (Battin universal variables approach) Numerical propagation (Cowell's method) Boundary-value problem, aka Lambert's problem (Izzo algorithm with multiple revolution) Orbital maneuvers (Hohmann's transfer, bielliptic transfer)

Slide 13

Slide 13 text

Analytical propagation: In [14]: florence.nu Out[14]: −0.40059292 rad In [15]: florence.propagate(1 * u.day).nu Out[15]: 5.9016222 rad

Slide 14

Slide 14 text

Numerical propagation: In [17]: from poliastro.twobody.propagation import propagate, cowell from poliastro.twobody.decorators import state_from_vector from poliastro.util import norm In [18]: @state_from_vector def accel(t0, ss): v = ss.v.value norm_v = norm(v) return 1e-5 * v / norm_v

Slide 15

Slide 15 text

In [19]: ss Out[19]: 7283 x 10293 km x 153.2 deg orbit around Earth (♁) In [20]: propagate(ss, 3 * u.day, method=cowell, rtol=1e-9, ad=accel, callback=results) Out[20]: 21188 x 24857 km x 153.2 deg orbit around Earth (♁)

Slide 16

Slide 16 text

In [24]: frame = OrbitPlotter3D() frame.set_attractor(Earth) frame.plot_trajectory(results.r) frame.show() Out[24]:

Slide 17

Slide 17 text

Lambert's problem: In [26]: date_launch = Time('2011-11-26 15:02', scale='utc').tdb date_arrival = Time('2012-08-06 05:17', scale='utc').tdb ss_e0 = Orbit.from_body_ephem(Earth, date_launch) ss_mf = Orbit.from_body_ephem(Mars, date_arrival) tof = date_arrival - date_launch (v0, v), = lambert(Sun.k, ss_e0.r, ss_mf.r, tof) In [27]: v0 Out[27]: [−29.291373, 14.532221, 5.4164177] km s In [29]: v Out[29]: [17.615607, − 10.998499, − 4.2080014] km s

Slide 18

Slide 18 text

Orbital maneuvers: In [30]: from poliastro.maneuver import Maneuver In [31]: ss_i = Orbit.circular(Earth, alt=700 * u.km) ss_i Out[31]: 7078 x 7078 km x 0.0 deg orbit around Earth (♁) In [32]: hoh = Maneuver.hohmann(ss_i, 36000 * u.km) hoh.get_total_cost() Out[32]: 3.6173999 km s

Slide 19

Slide 19 text

In [33]: frame = OrbitPlotter() ss_a, ss_f = ss_i.apply_maneuver(hoh, intermediate=True) frame.plot(ss_i, label="Initial orbit") frame.plot(ss_a, label="Transfer orbit") frame.plot(ss_f, label="Final orbit");

Slide 20

Slide 20 text

Good performance thanks to numba (Runs per second of the Izzo algorithm compiled or interpreted in several ways. Figure from ) https://indico.esa.int/indico/event/111/session/32/contribution/5 (https://indico.esa.int/indico/event/111/session/32/contribution/5)

Slide 21

Slide 21 text

3. Future work (2 minutes) 3. Future work (2 minutes) Continuous-thrust maneuvers and perturbations: Several optimal or quasi-optimal transfer maneuvers (studied as part of my �nal MSc project ) Solar pressure, atmospheric drag, third body effects... https://github.com/Juanlu001/pfc-uc3m (https://github.com/Juanlu001/pfc-uc3m) Improvements on propagation: Support for events (pericenter and apocenter passage, eclipses, custom events) Better convergence Storage of intermediate results (plotting, analysis) Better validation: So far, results from textbook assignments and scholarly papers: scarce and often not available in full precision Next step: validate against SPICE

Slide 22

Slide 22 text

...And, perhaps more importantly, attract contributors!

Slide 23

Slide 23 text

4. Conclusions 4. Conclusions poliastro offers an attractive, easy to use API while being "fast enough" Lots of functionality provided by a strong ecosystem of packages Interactive visualization is useful in preliminary analysis There is some work to do to add more features and improve stability We would love to have you as a contributor!

Slide 24

Slide 24 text

Per Python ad astra ➭ ✉ https://www.linkedin.com/in/juanluiscanor (https://www.linkedin.com /in/juanluiscanor) [email protected] (mailto:[email protected]) https://github.com/poliastro/oscw2017-talk (https://github.com/poliastro /oscw2017-talk)