Python Lab Reference

Every library, method, and pattern you need for the Orbital Mechanics labs.

65 entries across 8 libraries

🔢NumPy

The numerical workhorse. Every lab uses it for vectors, trigonometry and roots. Angles in NumPy are always in radians — convert deliberately.

import numpy as np

Creating & shaping arrays

np.array→ ndarray

np.array(object)

Build an array (a vector or matrix) from a Python list.

r = np.array([7000.0, 0.0, 0.0])   # position vector, km
v = np.array([0.0, 7.5, 1.0])       # velocity vector, km/s

Use floats (7000.0, not 7000) so later division stays continuous.

np.concatenate→ ndarray

np.concatenate([a, b])

Join arrays end-to-end — the standard way to pack a state vector.

y0 = np.concatenate([r0, v0])   # 6-element state [x,y,z, vx,vy,vz]

Pass a *list* of arrays. solve_ivp wants a flat 1-D state, and this is how you build it.

np.zeros→ ndarray

np.zeros(shape)

An array of zeros — handy for pre-allocating results.

accel = np.zeros(3)   # [0., 0., 0.]

np.linspace→ ndarray

np.linspace(start, stop, num)

num evenly spaced values from start to stop (both ends included).

times = np.linspace(0, T, 500)   # 500 sample times over one period

Pair with solve_ivp’s t_eval to get a smooth trajectory to plot.

Powers, roots & exponentials

np.sqrt→ float or ndarray

np.sqrt(x)

Square root. Works on scalars and whole arrays.

v = np.sqrt(mu * (2/r - 1/a))   # vis-viva speed

np.pi

The constant π (≈ 3.14159). Not a function — no parentheses.

T = 2 * np.pi * np.sqrt(a**3 / mu)   # orbital period

** (power)

Python’s exponent operator — use it, not np.power, for readability.

a_cubed = a ** 3

a**0.5 is another way to write a square root.

np.exp / np.log

np.exp(x) · np.log(x)

e^x and the natural logarithm (base e).

rho = rho0 * np.exp(-h / H)   # exponential atmosphere density

np.log is ln. For base-10 use np.log10.

Trigonometry & angles

np.radians / np.degrees

np.radians(deg) · np.degrees(rad)

Convert between degrees and radians. Your bridge to every trig call.

theta = np.radians(45)      # 0.785...
deg   = np.degrees(theta)   # back to 45.0

GOLDEN RULE: np.sin/np.cos expect radians. Convert on the way in, convert back on the way out.

np.sin / np.cos / np.tan→ float or ndarray

np.sin(x) — x in radians

Standard trig functions.

x = r * np.cos(nu)
y = r * np.sin(nu)   # perifocal position from true anomaly

np.arctan2→ radians in (−π, π]

np.arctan2(y, x)

Angle of the vector (x, y), correct in all four quadrants. Prefer over arctan.

Plain np.arctan(y/x) loses the quadrant and blows up when x=0. arctan2 takes y and x separately and returns the true angle in (−π, π].

lon = np.degrees(np.arctan2(r[1], r[0]))

Argument order is (y, x) — y first. A classic source of sign bugs.

np.arccos / np.arcsin / np.arctan

Inverse trig — return radians.

inc = np.degrees(np.arccos(h[2] / np.linalg.norm(h)))   # inclination

arccos/arcsin need input in [−1, 1]. Floating-point drift can give 1.0000001 → NaN. Guard with np.clip(x, -1, 1).

Vectors & linear algebra

np.linalg.norm→ float

np.linalg.norm(v)

Magnitude (length) of a vector — the single most-used call in these labs.

r_mag = np.linalg.norm(r)     # distance from Earth centre
speed = np.linalg.norm(v)     # scalar speed from velocity vector

For a matrix, pass axis=1 to get the norm of each row.

np.dot→ float

np.dot(a, b)

Dot product — projects one vector onto another; zero means perpendicular.

radial_speed = np.dot(r, v) / np.linalg.norm(r)

a @ b is the same thing and reads nicely: r @ v.

np.cross→ ndarray (3-vector)

np.cross(a, b)

Cross product — gives a vector perpendicular to both. The gateway to orbital elements.

The specific angular momentum h = r × v defines the orbital plane. Its direction gives inclination and node; its magnitude gives the orbit’s shape.

h = np.cross(r, v)                 # angular momentum vector
h_mag = np.linalg.norm(h)

📐Astropy Units

poliastro speaks in *quantities* — a number glued to a physical unit. This kills “km vs m” bugs, but you must attach units on the way in and strip them on the way out.

from astropy import units as u

Attaching units

value * u.<unit>→ Quantity

Multiply a number by a unit to make a Quantity.

a   = 7000 * u.km
inc = 51.6 * u.deg
dt  = 90 * u.min

This is the idiom poliastro constructors expect for every argument.

Common units

The units the labs use most.

u.km   u.m         # length
u.s    u.min  u.hour  # time
u.deg  u.rad        # angle
u.kg                # mass
u.one               # dimensionless (eccentricity!)

Eccentricity is unitless: pass ecc * u.one.

Compound units→ Quantity

Build derived units by dividing/multiplying.

speed = 7.5 * u.km / u.s
mu    = 398600 * u.km**3 / u.s**2

Converting & extracting

.to(unit)→ Quantity (in the new unit)

quantity.to(u.km / u.s)

Convert a Quantity to different units. Errors loudly if units are incompatible.

v_kms = np.linalg.norm(orb.v).to(u.km / u.s)
T_hr  = orb.period.to(u.hour)

.value→ float or ndarray

quantity.value

Strip the unit and get the raw float/array back.

e = orb.ecc.value          # 0.0006 as a plain float
print(f"e = {e:.4f}")

You almost always need .value before an f-string format spec like :.4f.

.to_value(unit)→ float or ndarray

quantity.to_value(u.km)

Convert *and* strip in one step — .to(unit).value combined.

r_km = orb.r_p.to_value(u.km)   # plain number in km

🌍Astropy Constants

Physical constants that already carry units and cited values — no more copy-pasting 398600.

from astropy.constants import GM_earth, R_earth

The constants you’ll use

GM_earth

Earth’s gravitational parameter μ = G·M (≈ 3.986×10¹⁴ m³/s²).

mu = GM_earth.to(u.km**3 / u.s**2)   # convert to km-based units
mu = GM_earth.value                  # raw number (SI, m³/s²)

Ships in SI (m³/s²). Convert to km³/s² if the rest of your math is in km.

R_earth

Earth’s equatorial radius (≈ 6378 km).

altitude = orb.r_p.to(u.km) - R_earth.to(u.km)

Distinguish radius-from-centre (r) vs altitude-above-surface (r − R_earth).

.value / .to(...)

Constants are Quantities — same rules as astropy units apply.

print(GM_earth)            # 3.986e+14 m3 / s2  (value + unit + reference)
print(GM_earth.value)      # 3.986004418e+14

🛰️poliastro

The high-level orbit toolkit. An Orbit object bundles the six classical elements and lets you read any derived property, propagate in time, and apply maneuvers — all unit-safe.

from poliastro.bodies import Earth
from poliastro.twobody.orbit import Orbit
from poliastro.maneuver import Maneuver
from poliastro.iod import izzo

Creating an orbit

Orbit.from_classical→ Orbit

Orbit.from_classical(Earth, a, ecc, inc, raan, argp, nu)

Build an orbit from the six classical orbital elements (COEs).

Arguments in order: semi-major axis a, eccentricity ecc, inclination inc, right ascension of ascending node (RAAN), argument of periapsis argp, true anomaly nu. Every argument is a Quantity.

orb = Orbit.from_classical(
    Earth,
    7000 * u.km,     # a
    0.01 * u.one,    # ecc
    51.6 * u.deg,    # inc
    0 * u.deg,       # raan
    0 * u.deg,       # argp
    0 * u.deg,       # nu (true anomaly)
)

nu = 0° is periapsis, nu = 180° is apoapsis.

Orbit.circular→ Orbit

Orbit.circular(Earth, alt, inc=0*u.deg)

Shortcut for a circular orbit at a given altitude above the surface.

iss = Orbit.circular(Earth, alt=420 * u.km, inc=51.6 * u.deg)

"alt" is altitude above Earth’s surface, not radius from centre.

Orbit.from_vectors→ Orbit

Orbit.from_vectors(Earth, r, v)

Build an orbit directly from position and velocity vectors.

r = [7000, 0, 0] * u.km
v = [0, 7.5, 1] * u.km / u.s
orb = Orbit.from_vectors(Earth, r, v)

The inverse of orb.rv(): vectors → orbit instead of orbit → vectors.

Reading orbit properties

orb.r · orb.v→ Quantity vectors

Current position and velocity vectors (as Quantities).

speed = np.linalg.norm(orb.v).to(u.km / u.s)

orb.a · orb.ecc · orb.inc

The classical elements: semi-major axis, eccentricity, inclination.

print(f"a = {orb.a:.1f}, e = {orb.ecc.value:.4f}")

Also available: orb.raan, orb.argp, orb.nu.

orb.r_p · orb.r_a→ Quantity

Periapsis and apoapsis radius (closest / farthest distance from centre).

print(f"perigee {orb.r_p.to(u.km):.0f}, apogee {orb.r_a.to(u.km):.0f}")

orb.period (orb.T)→ Quantity (time)

Time to complete one full orbit.

print(f"period = {orb.period.to(u.min):.1f}")

orb.energy→ Quantity

Specific orbital energy — negative for bound orbits, zero at escape.

print(f"energy = {orb.energy:.3f}")

Propagating & maneuvering

orb.propagate→ Orbit (a new one — the original is unchanged)

orb.propagate(time)

Advance the orbit forward (or back) in time; returns a new Orbit.

later = orb.propagate(45 * u.min)
print(later.nu.to(u.deg))   # where the satellite is now

Orbits are immutable: propagate/apply_maneuver return copies, they don’t mutate.

Maneuver.hohmann→ Maneuver

Maneuver.hohmann(orb, r_final)

Build the classic two-burn Hohmann transfer to a target circular radius.

man = Maneuver.hohmann(orb, 42164 * u.km)   # LEO → GEO
print(man.get_total_cost().to(u.km / u.s))  # total Δv
print(man.get_total_time().to(u.hour))      # transfer time

Maneuver.impulse→ Maneuver

Maneuver.impulse(dv)

A single instantaneous velocity change (Δv burn).

burn = Maneuver.impulse([0, 0.1, 0] * u.km / u.s)

orb.apply_maneuver→ Orbit

orb.apply_maneuver(maneuver)

Apply a maneuver and get the resulting orbit.

geo = orb.apply_maneuver(man)

izzo.lambert→ (v1, v2) velocity Quantities

izzo.lambert(k, r1, r2, tof)

Solve Lambert’s problem: the two velocities that fly from r1 to r2 in a given time.

k is the gravitational parameter (Earth.k), r1/r2 are position vectors, tof is time of flight. Returns the departure and arrival velocity vectors.

v1, v2 = izzo.lambert(Earth.k, r1, r2, tof)

This is how you plan an intercept or rendezvous between two known points.

∫SciPy · solve_ivp

The numerical integrator behind “propagate it yourself” labs. You write the equations of motion as a function; solve_ivp marches them through time.

from scipy.integrate import solve_ivp

Integrating equations of motion

solve_ivp→ a solution object (see below)

solve_ivp(fun, t_span, y0, t_eval=None, rtol=1e-3, atol=1e-6)

Integrate dy/dt = fun(t, y) from t_span[0] to t_span[1], starting at y0.

fun(t, y) must return the derivative of the state. For two-body motion the state is [x,y,z, vx,vy,vz] and its derivative is [vx,vy,vz, ax,ay,az].

def two_body(t, y):
    r = y[:3]
    v = y[3:]
    a = -mu * r / np.linalg.norm(r)**3
    return np.concatenate([v, a])

sol = solve_ivp(two_body, [0, T], y0,
                t_eval=times, rtol=1e-10, atol=1e-10)

sol.t · sol.y

The results: sol.t is the array of times, sol.y is the state at each time.

x = sol.y[0]      # x-position over time
vx = sol.y[3]     # x-velocity over time
final = sol.y[:, -1]   # full state at the last instant

sol.y is shaped (n_states, n_times) — each ROW is one variable through time.

t_eval

The specific times you want the answer sampled at (for smooth plots).

t_eval = np.linspace(0, T, 1000)

Without it you get whatever adaptive steps the solver chose — usually too few to plot.

rtol / atol

Relative and absolute error tolerances — smaller = more accurate but slower.

rtol=1e-10, atol=1e-10   # tight, good for orbit energy conservation

Loose tolerances make orbits visibly drift/spiral. For orbital mechanics, go tight.

📡SGP4 & Skyfield

Real-satellite tracking from TLEs (Two-Line Element sets). sgp4 is the low-level propagator; skyfield wraps it with friendly time, coordinate and ground-track handling.

from sgp4.api import Satrec
from skyfield.api import load, EarthSatellite, wgs84

sgp4 — raw TLE propagation

Satrec.twoline2rv→ Satrec

Satrec.twoline2rv(line1, line2)

Parse a two-line element set into a propagatable satellite record.

sat = Satrec.twoline2rv(line1, line2)

sat.sgp4→ (error_code, r, v) — error 0 means success

sat.sgp4(jd, fr)

Propagate to a Julian date (jd + fractional day fr); returns error, r, v.

e, r, v = sat.sgp4(jd, fr)
if e == 0:
    print("position (km):", r)

r and v are plain tuples in the TEME frame, km and km/s.

skyfield — friendly tracking

load.timescale→ Timescale

ts = load.timescale()

Create a timescale, then build times from it.

ts = load.timescale()
t = ts.utc(2026, 7, 9, 12, 0, 0)

EarthSatellite→ EarthSatellite

EarthSatellite(line1, line2, name, ts)

A satellite built from a TLE, tied to a timescale.

sat = EarthSatellite(line1, line2, "ISS", ts)

sat.at(t)→ Geocentric position

The satellite’s geocentric position at time t.

geo = sat.at(t)
print(geo.position.km)   # [x, y, z] in km

wgs84.subpoint→ GeographicPosition (latitude, longitude, elevation)

wgs84.subpoint(geo)

The point on Earth directly below the satellite — its ground track.

sp = wgs84.subpoint(sat.at(t))
print(sp.latitude.degrees, sp.longitude.degrees)

🐍Python Essentials

The plain-Python patterns that show up in every lab — defining functions, printing clean output, and formatting numbers so the auto-grader sees what it expects.

f-strings & number formatting

f"..."

Format string — drop a variable straight into text with {var}.

a = 7000.0
print(f"semi-major axis = {a} km")

The f before the quote is required. Without it {a} prints literally.

:.2f (fixed decimals)

Show a number with a set number of decimal places.

print(f"speed = {v:.2f} km/s")   # 7.55
print(f"e = {ecc:.4f}")          # 0.0006

Rounds for display only — the underlying value is unchanged.

:.3e (scientific)

Scientific notation with N significant digits.

print(f"mu = {mu:.3e}")   # 3.986e+05

:>10 / :<10 (alignment)

Pad to a width for tidy columns — right (>) or left (<) aligned.

print(f"{name:<12}{value:>8.1f}")

Great for building the aligned tables the labs often ask for.

Functions & structure

def→ whatever you `return`

def name(arg, arg2=default):

Define a reusable function. Indentation (4 spaces) marks the body.

def period(a, mu=398600):
    return 2 * np.pi * np.sqrt(a**3 / mu)

print(period(7000))

A function with no return gives None. Remember to return your result.

return

Hand a value back to the caller and exit the function.

return v_periapsis, v_apoapsis   # returns a tuple of two values

print

Write to stdout — this is what the auto-grader reads to check your answer.

print(f"Result: {answer:.2f}")

A lab passes when the expected text appears in your printed output, so match wording/units.

for loop

for item in sequence:

Repeat a block over a list, range, or array.

for alt in [200, 400, 800]:
    orb = Orbit.circular(Earth, alt * u.km)
    print(alt, orb.period.to(u.min))

📖Orbital Cookbook

Recurring recipes and the traps that cost students the most time. Skim this before you start debugging — the bug is usually here.

Go-to recipes

Vis-viva (speed at any point)

The single most useful orbital equation: speed from radius and semi-major axis.

v = np.sqrt(mu * (2 / r - 1 / a))

r = current distance from centre, a = semi-major axis, mu = GM. Circular orbit ⇒ r = a.

Angular momentum vector

h = r × v defines the orbital plane; its z-component sets inclination.

h = np.cross(r, v)
inc = np.degrees(np.arccos(h[2] / np.linalg.norm(h)))

Circular orbit speed

Speed to hold a circular orbit at radius r.

v_circ = np.sqrt(mu / r)

Escape velocity

Speed needed to just barely escape from radius r.

v_esc = np.sqrt(2 * mu / r)   # = sqrt(2) * v_circ

Period (Kepler’s third law)

Orbital period from semi-major axis.

T = 2 * np.pi * np.sqrt(a**3 / mu)

Common traps

Degrees vs radians

NumPy trig is radians; poliastro wants u.deg. Mixing them silently wrecks answers.

# WRONG: np.sin(45)      → sine of 45 radians
# RIGHT: np.sin(np.radians(45))

Symptom: numbers off by a weird factor with no error message.

arccos domain (NaN)

Rounding can push arccos input just past ±1, producing NaN.

inc = np.degrees(np.arccos(np.clip(h[2] / h_mag, -1, 1)))

Symptom: result is nan. Fix: wrap the argument in np.clip(x, -1, 1).

Forgetting .value

A Quantity won’t accept an f-string format spec like :.4f.

# ERROR: f"{orb.ecc:.4f}"       (unless it happens to work)
print(f"{orb.ecc.value:.4f}")   # strip the unit first

radius vs altitude

r is measured from Earth’s centre; altitude is above the surface.

altitude = r - R_earth.to(u.km)   # subtract the radius

A "400 km orbit" means 400 km altitude ⇒ r ≈ 6778 km.

km vs m consistency

Pick one length unit and keep mu in matching units.

mu = 398600        # km³/s²  → use km everywhere
mu = 3.986e14      # m³/s²   → use m everywhere

Mixing them is the #1 cause of answers off by 10³ or 10⁶.