Appendix I: Mission Design Toolkit — the Complete astrotools Package

If you followed the optional "ambitious reader" track of the Design Your Mission project, you built a small Python package called astrotools one piece at a time, adding a module (or growing an existing one) in each chapter's Mission Design Checkpoint. This appendix reproduces the finished package in full, organized by module, so that a reader who skipped the build-it-yourself track can still read the completed code — or copy it whole and start experimenting. It is an illustrative teaching toolkit, not production software: like every listing in this book it was never executed at build time, and the hand-traced # Expected output: values that verify each function live in the chapters where the code was introduced (each chapter's code/project-checkpoint.py). The code targets Python 3.10+ and depends only on the four libraries pinned in requirements.txt (numpy, scipy, matplotlib, astropy); most modules need only numpy or the standard-library math.

A standing convention: every listing below is the importable body of a module — its constants and functions, with docstrings. The interactive demonstrations and their hand-derived # Expected output: comments (the if __name__ == "__main__": blocks) stay in each chapter, where they are worked and checked against reality. Modules import one another by flat name, exactly as the chapters do (for example from rocket import delta_v), and every module that needs the standard gravity defines the same constant, G0 = 9.80665 m/s². Functions that take a gravitational parameter mu are unit-agnostic: feed them consistent units (SI throughout, or km and km/s throughout) and the result comes back in the same system.

How the package is organized

The eleven canonical modules below are the spine of the toolkit; their function signatures were frozen at the start of the build so that a piece written in one chapter composes cleanly with a piece from another. A set of additional utilities — smaller modules and standalone Mission-Design helpers contributed by chapters off the core schedule — follows.

rocket.py

The rocket equation and its relatives: specific impulse to exhaust velocity, delta-v and mass ratio, a launch-ascent estimate, and the staging helpers. Built in Chapter 3 (core), extended in Chapter 4 (ascent_delta_v) and Chapter 22 (staging); the $I_{sp}\leftrightarrow v_e$ relationship is also formalized in Chapter 16's propulsion.py.

"""astrotools/rocket.py -- the rocket equation, Isp <-> v_e, and staging.

Contributed across Chapters 3 (core), 4 (ascent_delta_v), and 22 (staging).
Illustrative teaching code; never executed at build (see Appendix I).
"""
import math

G0 = 9.80665  # standard gravity, m/s^2


def isp_to_ve(isp):
    """Specific impulse (s) -> effective exhaust velocity (m/s)."""
    return isp * G0


def ve_to_isp(ve):
    """Effective exhaust velocity (m/s) -> specific impulse (s)."""
    return ve / G0


def delta_v(ve, m0, mf):
    """Tsiolkovsky rocket equation. Returns delta-v (m/s)."""
    return ve * math.log(m0 / mf)


def mass_ratio(dv, ve):
    """Mass ratio m0/mf needed for a target delta-v (m/s) at exhaust velocity ve (m/s)."""
    return math.exp(dv / ve)


def propellant_fraction(dv, ve):
    """Propellant fraction m_p / m_0 needed for a target delta-v (m/s)."""
    return 1 - math.exp(-dv / ve)


def stage_delta_v(stages):
    """Total delta-v (m/s) of a multistage rocket. `stages` is a list of (ve, m0, mf)
    tuples, each including everything above it. (Chapter 22 inlines the same sum.)"""
    return sum(delta_v(ve, m0, mf) for (ve, m0, mf) in stages)


# --- ascent increment (Chapter 4) --------------------------------------------
def ascent_delta_v(v_orbit, gravity_loss=1500.0, drag_loss=100.0, rotation_bonus=0.0):
    """Launch delta-v (m/s) to reach an orbit of tangential speed v_orbit (m/s).

    Adds representative gravity and drag losses; subtracts an optional Earth-rotation
    credit (up to ~465 m/s, due east from the equator). The losses are order-of-magnitude
    figures -- refine with a trajectory model for a real design.
    """
    return v_orbit + gravity_loss + drag_loss - rotation_bonus


# --- staging increment (Chapter 22) ------------------------------------------
def payload_fraction(dv_total, ve, eps, n):
    """Overall payload fraction of n IDENTICAL stages (exhaust velocity ve, structural
    coefficient eps) splitting dv_total equally. Returns nan if a stage is unbuildable."""
    R = math.exp((dv_total / n) / ve)
    if R >= 1.0 / eps:                      # exceeds the structural ceiling 1/eps
        return float("nan")
    pi = (1 - eps * R) / (R * (1 - eps))
    return pi ** n


def stage_mass(m_above, dv_stage, ve, eps):
    """Loaded mass (structure + propellant, kg) of one stage carrying m_above on top:
    m_stage = m_above (R - 1) / (1 - eps*R), R = exp(dv/ve)."""
    R = math.exp(dv_stage / ve)
    return m_above * (R - 1) / (1 - eps * R)

orbits.py

Orbital energy and speeds, Kepler's third law and equation, the conversion from classical elements to inertial state vectors, and the secular $J_2$ perturbation helpers. Built in Chapter 6 (energy) and Chapter 8 (Kepler); Chapter 2's preview contributes escape_velocity and surface_gravity; the perifocal-to-inertial rotation that finishes elements_to_rv was scheduled for Chapter 9; and Chapter 12 adds the $J_2$ helpers.

"""astrotools/orbits.py -- energy, Kepler, elements <-> state vectors, and J2.

Contributed across Chapters 2 (escape_velocity, surface_gravity -- the Chapter 2
'astrotools preview'), 6 (energy), 8 (Kepler III, Kepler's equation, the perifocal
sketch of elements_to_rv), 9 (the perifocal->inertial rotation), and 12 (secular J2).
The energy/Kepler functions are unit-agnostic; the J2 helpers default to Earth in km.
"""
import math
import numpy as np

G0 = 9.80665  # standard gravity, m/s^2


# --- Chapter 2 (astrotools preview) ------------------------------------------
def escape_velocity(mu, r):
    """Escape velocity from radius r for gravitational parameter mu: sqrt(2*mu/r).
    (Introduced in the Chapter 2 preview, which assigns it to orbits.py.)"""
    return math.sqrt(2 * mu / r)


def surface_gravity(mu, r):
    """Gravitational acceleration at radius r: mu/r**2. (Chapter 2 preview.)"""
    return mu / r**2


# --- Chapter 6 (orbital energy) ----------------------------------------------
def circular_velocity(mu, r):
    """Circular orbital speed at radius r: sqrt(mu/r)."""
    return math.sqrt(mu / r)


def specific_energy(mu, a):
    """Specific orbital energy: -mu/(2a). <0 bound, 0 parabolic, >0 hyperbolic."""
    return -mu / (2 * a)


def vis_viva(mu, r, a):
    """Speed at radius r on an orbit of semi-major axis a: sqrt(mu*(2/r - 1/a)).
    The workhorse of orbital mechanics."""
    return math.sqrt(mu * (2 / r - 1 / a))


# --- Chapter 8 (Kepler's third law + time of flight) -------------------------
def period(mu, a):
    """Orbital period from Kepler's third law: T = 2*pi*sqrt(a**3/mu)."""
    return 2 * math.pi * math.sqrt(a**3 / mu)


def mean_motion(mu, a):
    """Mean motion (rad/time): n = 2*pi/T = sqrt(mu/a**3)."""
    return math.sqrt(mu / a**3)


def solve_kepler(M, e, tol=1e-10):
    """Solve Kepler's equation M = E - e*sin(E) for the eccentric anomaly E (rad) by
    Newton's method, starting from E0 = M + e*sin(M)."""
    E = M + e * math.sin(M)
    for _ in range(50):
        dE = (E - e * math.sin(E) - M) / (1 - e * math.cos(E))
        E -= dE
        if abs(dE) < tol:
            break
    return E


# --- Chapters 8 + 9 (classical elements -> inertial state vectors) -----------
def _rot_z(theta):
    """Active rotation matrix about the z-axis by theta (rad)."""
    c, s = math.cos(theta), math.sin(theta)
    return np.array([[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]])


def _rot_x(theta):
    """Active rotation matrix about the x-axis by theta (rad)."""
    c, s = math.cos(theta), math.sin(theta)
    return np.array([[1.0, 0.0, 0.0], [0.0, c, -s], [0.0, s, c]])


def elements_to_rv(mu, a, e, i, raan, argp, nu):
    """Convert classical orbital elements to inertial position and velocity vectors
    (angles in rad; lengths consistent with mu). Chapter 8 builds the perifocal (PQW)
    r and v from the orbit equation and its derivative; Chapter 9 rotates PQW -> inertial
    with the 3-1-3 sequence R = Rz(raan) @ Rx(i) @ Rz(argp). Returns (r, v) as arrays."""
    p = a * (1 - e**2)                                   # semi-latus rectum
    r = p / (1 + e * math.cos(nu))                       # radius (orbit equation, Ch. 8)
    r_pqw = np.array([r * math.cos(nu), r * math.sin(nu), 0.0])
    s = math.sqrt(mu / p)
    v_pqw = np.array([-s * math.sin(nu), s * (e + math.cos(nu)), 0.0])
    R = _rot_z(raan) @ _rot_x(i) @ _rot_z(argp)          # perifocal -> inertial (Ch. 9)
    return R @ r_pqw, R @ v_pqw


# --- Chapter 12 (secular J2 perturbations; Earth defaults in km) --------------
MU_EARTH = 3.986e5      # km^3/s^2  (default for the J2 helpers below)
RE_EARTH = 6378.0       # km, equatorial radius (used with J2)
J2_EARTH = 1.0826e-3


def j2_nodal_rate(a, e, i_deg, mu=MU_EARTH, R=RE_EARTH, J2=J2_EARTH):
    """Secular nodal-regression rate d(RAAN)/dt in deg/day (Section 12.2)."""
    n = math.sqrt(mu / a**3)                                       # mean motion, rad/s
    rate = -1.5 * n * J2 * (R / a)**2 / (1 - e * e)**2 * math.cos(math.radians(i_deg))
    return math.degrees(rate) * 86400.0


def sun_sync_inclination(a, mu=MU_EARTH, R=RE_EARTH, J2=J2_EARTH):
    """Inclination (deg) making a circular orbit sun-synchronous (+0.9856 deg/day)."""
    n = math.sqrt(mu / a**3)
    raan_req = math.radians(360.0 / 365.2422) / 86400.0           # rad/s, eastward
    return math.degrees(math.acos(-raan_req / (1.5 * n * J2 * (R / a)**2)))


def station_keeping_dv(ns_per_year, ew_per_year, years):
    """Total station-keeping delta-v (m/s) over a mission life (Section 12.6)."""
    return (ns_per_year + ew_per_year) * years


def propellant_mass(m0, dv, isp):
    """Propellant mass (kg) for a delta-v (m/s) via the rocket equation, from wet mass m0
    (kg) and Isp (s). Mirrors rocket.propellant_fraction, returning kilograms."""
    ve = isp * G0
    return m0 * (1 - math.exp(-dv / ve))

maneuvers.py

Impulsive maneuvers between orbits: the two-burn Hohmann transfer, a pure plane change, and a helper that rolls a list of legs into a total. Built in Chapter 10; it builds on the orbital speeds of orbits.py and feeds mission.py (Chapter 29).

"""astrotools/maneuvers.py -- impulsive maneuvers: Hohmann, plane change, budgets.

Built in Chapter 10. Unit-agnostic: work in km/km-per-s or in SI throughout.
"""
import math


def hohmann(mu, r1, r2):
    """Two-burn Hohmann transfer between circular orbits r1 -> r2 (same length unit).
    Returns (dv1, dv2, dv_total) in the speed unit implied by mu and r."""
    a_t = (r1 + r2) / 2.0
    v_c1 = math.sqrt(mu / r1)
    v_c2 = math.sqrt(mu / r2)
    v_p = math.sqrt(mu * (2.0 / r1 - 1.0 / a_t))   # transfer-ellipse perigee speed
    v_a = math.sqrt(mu * (2.0 / r2 - 1.0 / a_t))   # transfer-ellipse apogee speed
    dv1 = abs(v_p - v_c1)                           # leave inner circle
    dv2 = abs(v_c2 - v_a)                           # join outer circle
    return dv1, dv2, dv1 + dv2


def plane_change(v, di_deg):
    """Delta-v for a pure plane change of di_deg degrees at orbital speed v:
    dv = 2 v sin(di/2)."""
    return 2.0 * v * math.sin(math.radians(di_deg) / 2.0)


def dv_budget(legs):
    """Roll up a list of (label, dv) maneuver legs into a total delta-v budget."""
    return sum(dv for _, dv in legs)

interplanetary.py

Heliocentric transfers: the interplanetary Hohmann (returning the hyperbolic-excess speeds a spacecraft needs at each planet plus the coast time), the synodic period that sets launch windows, and the departure characteristic energy $C_3$. Built in Chapter 11; it extends maneuvers.py from geocentric to heliocentric transfers.

"""astrotools/interplanetary.py -- heliocentric transfers, windows, and C3.

Built in Chapter 11. Distances in km, speeds in km/s, times in seconds.
"""
import math

MU_SUN = 1.327e11    # km^3/s^2, Sun's gravitational parameter
AU = 1.496e8         # km per astronomical unit


def hohmann_transfer(mu, r1, r2):
    """Heliocentric Hohmann transfer r1 -> r2 (km). Returns
    (v_inf_depart, v_inf_arrive, tof_seconds): the hyperbolic-excess speeds (km/s)
    the spacecraft needs relative to each planet, and the coast time (s)."""
    v1 = math.sqrt(mu / r1)                  # departure-planet circular speed
    v2 = math.sqrt(mu / r2)                  # arrival-planet   circular speed
    a = (r1 + r2) / 2                        # transfer-ellipse semi-major axis
    vp = math.sqrt(mu * (2 / r1 - 1 / a))    # perihelion (departure) speed
    va = math.sqrt(mu * (2 / r2 - 1 / a))    # aphelion   (arrival)   speed
    return abs(vp - v1), abs(v2 - va), math.pi * math.sqrt(a ** 3 / mu)


def synodic_period(T1, T2):
    """Time between launch windows for two orbits of periods T1, T2 (same units)."""
    return 1.0 / abs(1.0 / T1 - 1.0 / T2)


def c3_required(v_inf):
    """Characteristic energy C3 (km^2/s^2) from hyperbolic excess speed v_inf (km/s):
    C3 = v_inf**2. The number a launch vehicle is rated against."""
    return v_inf ** 2

attitude.py

Spacecraft orientation with quaternions and direction-cosine matrices. Built in Chapter 14. Conventions are frozen: Hamilton quaternions, scalar-first $q = [w, x, y, z]$ with unit norm; a passive DCM $R_{B/N}$ that carries reference-frame components into body-frame components; and a 3-2-1 (yaw-pitch-roll) Euler sequence for euler_to_quat.

"""astrotools/attitude.py -- quaternions and direction-cosine matrices.

Built in Chapter 14. CONVENTIONS: Hamilton quaternions, SCALAR-FIRST q = [w, x, y, z],
unit norm; DCM is the PASSIVE transform R_{B/N}; euler_to_quat uses a 3-2-1 sequence.
"""
import numpy as np


def quat_multiply(q1, q2):
    """Hamilton product q1 (x) q2 of two scalar-first quaternions. Composes rotations."""
    w1, x1, y1, z1 = q1
    w2, x2, y2, z2 = q2
    return np.array([
        w1*w2 - x1*x2 - y1*y2 - z1*z2,
        w1*x2 + x1*w2 + y1*z2 - z1*y2,
        w1*y2 - x1*z2 + y1*w2 + z1*x2,
        w1*z2 + x1*y2 - y1*x2 + z1*w2,
    ])


def dcm_from_quat(q):
    """Passive DCM R_{B/N} from a unit, scalar-first quaternion."""
    w, x, y, z = q
    return np.array([
        [1 - 2*(y*y + z*z),     2*(x*y + w*z),     2*(x*z - w*y)],
        [    2*(x*y - w*z), 1 - 2*(x*x + z*z),     2*(y*z + w*x)],
        [    2*(x*z + w*y),     2*(y*z - w*x), 1 - 2*(x*x + y*y)],
    ])


def _axis_quat(axis, angle):
    """Unit quaternion for a rotation of `angle` (rad) about a unit `axis`."""
    half = angle / 2.0
    c, sn = np.cos(half), np.sin(half)
    return np.array([c, sn*axis[0], sn*axis[1], sn*axis[2]])


def euler_to_quat(yaw, pitch, roll):
    """Quaternion for a 3-2-1 (yaw-pitch-roll) Euler sequence, angles in radians."""
    q_yaw = _axis_quat((0.0, 0.0, 1.0), yaw)
    q_pitch = _axis_quat((0.0, 1.0, 0.0), pitch)
    q_roll = _axis_quat((1.0, 0.0, 0.0), roll)
    return quat_multiply(q_yaw, quat_multiply(q_pitch, q_roll))

propulsion.py

Where thrust and specific impulse come from: the thrust equation, the ideal exhaust velocity set by chamber temperature and molecular weight, the nozzle expansion ratio, a chemical-propellant selection table, and the electric-propulsion trio. Built across Chapters 16 (thrust, $I_{sp}$), 17 (propellant selection), 18–19 (exit_velocity, expansion_ratio), and 20 (electric propulsion).

"""astrotools/propulsion.py -- thrust, exit velocity, nozzles, electric thrusters.

Contributed across Chapters 16 (thrust, specific impulse), 17 (propellant selection),
18-19 (exit_velocity, expansion_ratio), and 20 (electric propulsion). Chapter 21 (a
survey) adds no code; its thrust-from-power relation is ep_thrust below.
"""
import math

G0 = 9.80665  # standard gravity, m/s^2
R_U = 8.314   # universal gas constant, J/(mol*K)


# --- Chapter 16 (thrust and specific impulse) --------------------------------
def thrust(mdot, ve, pe, pa, ae):
    """Thrust (N) from the thrust equation: the mdot*ve momentum term plus the
    (pe - pa)*ae pressure term. mdot kg/s, ve m/s, pe/pa Pa, ae m^2."""
    return mdot * ve + (pe - pa) * ae


def effective_exhaust_velocity(isp):
    """Specific impulse (s) -> effective exhaust velocity c (m/s).
    (Same relation as rocket.isp_to_ve, kept here for propulsion self-containment.)"""
    return isp * G0


def specific_impulse(F, mdot):
    """Rigorous specific impulse (s): thrust per unit propellant weight-flow, F/(mdot*g0)."""
    return F / (mdot * G0)


def thrust_to_weight(F, m, g=9.81):
    """Thrust-to-weight ratio (dimensionless); liftoff from a surface needs > 1.
    g is the LOCAL gravitational acceleration, not the Isp standard gravity G0."""
    return F / (m * g)


def total_impulse(ve, mp):
    """Total impulse (N.s) = effective exhaust velocity (m/s) * propellant mass (kg)."""
    return ve * mp


# --- Chapter 17 (chemical-propellant selection) ------------------------------
# name: (vac Isp low, vac Isp high [s], bulk density [kg/m^3], storable?, restartable?)
PROPELLANTS = {
    "hydrolox":   (450, 465,  360, False, True),
    "methalox":   (350, 380,  830, False, True),
    "kerolox":    (300, 353, 1020, False, True),
    "hypergolic": (315, 340, 1180, True,  True),
    "solid":      (250, 285, 1750, True,  False),
}


def select(storable_needed, restart_needed):
    """Propellant classes meeting the mission's constraints, best (high) Isp first."""
    ok = []
    for name, (isp_lo, isp_hi, rho, storable, restart) in PROPELLANTS.items():
        if storable_needed and not storable:
            continue
        if restart_needed and not restart:
            continue
        ok.append((isp_hi, name))
    return [name for _isp, name in sorted(ok, reverse=True)]


# --- Chapters 18-19 (nozzle flow) --------------------------------------------
def exit_velocity(gamma, Tc, M, pe, pc):
    """Ideal (isentropic) exhaust velocity (m/s):
    sqrt( (2*gamma/(gamma-1)) * (R_u*Tc/M) * (1 - (pe/pc)**((gamma-1)/gamma)) ).
    gamma = ratio of specific heats, Tc = chamber temperature (K), M = mean exhaust
    molar mass (kg/mol), pe/pc = exit/chamber pressure (shared unit). Introduced in
    Chapter 18; the nozzle-theory home is Chapter 19 (identical formula)."""
    R = R_U / M
    bracket = 1.0 - (pe / pc) ** ((gamma - 1) / gamma)
    return math.sqrt((2 * gamma / (gamma - 1)) * R * Tc * bracket)


def expansion_ratio(gamma, pe, pc):
    """Nozzle area ratio A_e/A_t to expand from chamber pressure pc to exit pressure pe,
    for isentropic flow choked at the throat."""
    g = gamma
    num = ((2 / (g + 1)) ** (1 / (g - 1))) * (pc / pe) ** (1 / g)
    den = math.sqrt(((g + 1) / (g - 1)) * (1 - (pe / pc) ** ((g - 1) / g)))
    return num / den


# --- Chapter 20 (electric propulsion) ----------------------------------------
def ep_thrust(power_w, isp_s, efficiency):
    """Thrust (N) of an electric thruster from input power (W), Isp (s), and total
    efficiency (0-1): F = 2*eta*P / v_e, with v_e = Isp*g0. At fixed power, higher Isp
    means lower thrust. (Chapter 21's thrust_from_power is the same relation.)"""
    ve = isp_s * G0
    return 2 * efficiency * power_w / ve


def ep_mass_flow(thrust_n, isp_s):
    """Propellant mass flow (kg/s): mdot = F / v_e."""
    return thrust_n / (isp_s * G0)


def ep_burn_time(prop_kg, thrust_n, isp_s):
    """Continuous-thrust time (s) to expend prop_kg of propellant: t = m_p / mdot."""
    return prop_kg / ep_mass_flow(thrust_n, isp_s)

power.py

Sizing the electrical power subsystem: solar flux with distance, the array output needed to run the loads and recharge through eclipse, the resulting array area at end of life, and the battery capacity. Built in Chapter 25, alongside its siblings thermal.py (Chapter 24) and comms.py (Chapter 26).

"""astrotools/power.py -- solar-array sizing, eclipse, and battery.

Built in Chapter 25. Efficiency defaults follow SMAD conventions.
"""

S_1AU = 1361.0  # solar constant at 1 AU, W/m^2


def solar_flux(au):
    """Solar flux (W/m^2) at a distance of `au` astronomical units.
    (Shared verbatim with thermal.solar_flux.)"""
    return S_1AU / au**2


def required_array_power(p_day, t_day, p_ecl, t_ecl, x_day=0.85, x_ecl=0.65):
    """Array output power (W) needed to run the loads and recharge for eclipse.
    Times in consistent units; x_day/x_ecl are power-path efficiencies (SMAD)."""
    return (p_day * t_day / x_day + p_ecl * t_ecl / x_ecl) / t_day


def array_area(p_required, flux=S_1AU, eff=0.30, cos_theta=0.90, life_factor=1.0):
    """Solar-array area (m^2) to produce p_required watts, given cell efficiency,
    sun-angle cosine loss, and an end-of-life degradation factor (<=1)."""
    return p_required / (flux * eff * cos_theta * life_factor)


def battery_capacity_wh(p_ecl, t_ecl_h, dod=0.25, eff=0.90):
    """Battery capacity (W*h) to deliver an eclipse, given depth of discharge and
    discharge-path efficiency."""
    return p_ecl * t_ecl_h / (dod * eff)

thermal.py

Radiative thermal balance: solar flux with distance, the equilibrium temperature of a surface from its absorptivity, emissivity, areas, and internal power, and the radiator area needed to reject a given heat load. Built in Chapter 24; it pairs with power.py, since nearly every watt the power system delivers becomes waste heat this module must reject.

"""astrotools/thermal.py -- radiative equilibrium temperature and radiator sizing.

Built in Chapter 24. Pairs with power.py (Chapter 25).
"""

SIGMA = 5.670e-8    # Stefan-Boltzmann constant, W m^-2 K^-4
S_1AU = 1361.0      # solar flux at 1 AU, W/m^2


def solar_flux(au):
    """Solar flux (W/m^2) at heliocentric distance `au` in astronomical units.
    (Shared verbatim with power.solar_flux.)"""
    return S_1AU / au**2


def equilibrium_temperature(alpha, eps, a_sun, a_rad, flux=S_1AU, q_internal=0.0):
    """Equilibrium temperature (K) from radiative balance.

    alpha: solar absorptivity; eps: IR emissivity; a_sun: sunlit area (m^2);
    a_rad: radiating area (m^2); flux: solar flux (W/m^2); q_internal: internal power (W).
    """
    absorbed = alpha * flux * a_sun + q_internal
    return (absorbed / (eps * SIGMA * a_rad)) ** 0.25


def radiator_area(q_reject, t_rad, eps=0.85, alpha=0.0, flux=0.0):
    """Radiator area (m^2) to reject q_reject (W) at t_rad (K), net of absorbed sunlight."""
    net_per_m2 = eps * SIGMA * t_rad**4 - alpha * flux
    return q_reject / net_per_m2

comms.py

The link budget, done entirely in decibels so that a chain of gains and losses is just a sum: transmit power, free-space path loss, parabolic-dish gain, noise spectral density, and the received power they combine into. Built in Chapter 26.

"""astrotools/comms.py -- the decibel link budget.

Built in Chapter 26; works entirely in decibels so a link budget is just a sum.
"""
import math

C = 2.998e8                  # speed of light, m/s
K_BOLTZMANN = 1.380649e-23   # J/K


def watts_to_dbw(power_w):
    """Power in watts -> dBW."""
    return 10 * math.log10(power_w)


def fspl_db(dist_m, freq_hz):
    """Free-space path loss (dB). The ratio is squared -> a factor of 20."""
    wavelength = C / freq_hz
    return 20 * math.log10(4 * math.pi * dist_m / wavelength)


def dish_gain_db(diameter_m, freq_hz, efficiency=0.6):
    """Parabolic-dish gain (dBi) from diameter and frequency."""
    wavelength = C / freq_hz
    return 10 * math.log10(efficiency * (math.pi * diameter_m / wavelength) ** 2)


def noise_density_dbw_hz(ts_k):
    """Noise spectral density N0 = k*Ts, in dBW/Hz."""
    return 10 * math.log10(K_BOLTZMANN) + 10 * math.log10(ts_k)


def link_budget(pt_w, gt_db, gr_db, freq_hz, dist_m, extra_loss_db=0.0):
    """Received power (dBW) for a one-way link -- add the decibel column."""
    eirp = watts_to_dbw(pt_w) + gt_db
    return eirp - fspl_db(dist_m, freq_hz) + gr_db - extra_loss_db

structures.py

Mass budgeting — the discipline of the enemy — and a structural margin check: subsystem dry mass, a growth margin, the combined budget, and the margin of safety. Built in Chapter 23; its mass budget feeds power sizing (Chapter 25) and the synthesis in Chapter 29.

"""astrotools/structures.py -- mass budgeting and structural margins.

Built in Chapter 23. Signatures are called by later chapters, so they stay stable.
"""


def dry_mass(subsystems):
    """Total dry mass (kg) from a dict {subsystem_name: mass_kg}."""
    return sum(subsystems.values())


def apply_margin(mass_kg, margin_frac):
    """Mass carried after a growth margin (e.g. 0.25 for 25%)."""
    return mass_kg * (1.0 + margin_frac)


def mass_budget(subsystems, margin_frac=0.25):
    """Return (dry_kg, allocated_kg) for a subsystem mass dict and a margin."""
    dry = dry_mass(subsystems)
    return dry, apply_margin(dry, margin_frac)


def margin_of_safety(sigma_allow, sigma_limit, fos):
    """Structural margin of safety; positive => passes with reserve. sigma_allow = yield
    or ultimate strength; sigma_limit = limit (max) stress; fos = factor of safety."""
    return sigma_allow / (fos * sigma_limit) - 1.0

mission.py

The module that turns a mission into a vehicle: roll up a delta-v budget with per-leg margins, invert the rocket equation to size the vehicle, confirm the result closes on a launcher, lay out the phase timeline, and weigh reuse economics. Built in Chapter 29 (roll_up_dv, size_vehicle), extended in Chapter 38 (reuse economics) and Chapter 40 (closes, timeline). At the capstone it composes the whole package, importing circular_velocity from orbits.py and hohmann from maneuvers.py.

"""astrotools/mission.py -- roll up the budget, size the vehicle, check closure, time it.

Built in Chapter 29 (roll_up_dv, size_vehicle); extended in Chapter 38 (reuse economics)
and Chapter 40 (closes, timeline). The Chapter 40 capstone composes this module with
orbits.circular_velocity and maneuvers.hohmann to size a mission end to end.
"""
import math

G0 = 9.80665  # standard gravity, m/s^2


def roll_up_dv(budget):
    """Roll up a delta-v budget. `budget` is a list of (label, dv, margin) tuples:
    dv in m/s, margin a fraction (0.05 = 5%). Returns (ideal_total, margined_total)."""
    ideal = sum(dv for _, dv, _ in budget)
    margined = sum(dv * (1.0 + m) for _, dv, m in budget)
    return ideal, margined


def size_vehicle(dv, isp, payload):
    """Invert the rocket equation to size a vehicle. dv (m/s), isp (s),
    payload = delivered dry mass m_f (kg). Returns (m_prop, m_wet) in kg."""
    ve = isp * G0
    m_wet = math.exp(dv / ve) * payload      # m0 = (m0/mf) * mf
    return m_wet - payload, m_wet


def closes(m_wet, capacity):
    """Closure check (Chapters 30/40): does the sized wet mass fit under a launcher's
    capacity? Returns (fits, how_many_fit_per_launch)."""
    return m_wet < capacity, int(capacity // m_wet)


def timeline(phases):
    """Build a mission timeline from a list of (phase_name, duration) pairs. Returns a
    list of (phase_name, start, end) with cumulative times in the durations' own unit.
    (Canonical mission.py signature; completed at the Chapter 40 capstone.)"""
    schedule, t = [], 0.0
    for name, duration in phases:
        schedule.append((name, t, t + duration))
        t += duration
    return schedule


# --- Chapter 38 (launch economics under reuse) -------------------------------
def reuse_cost_per_flight(M, R, F, N):
    """Cost per flight ($M) of a reusable stage flown N times: build cost M amortized over
    N flights, plus per-flight refurbishment R and the fixed cost F reuse never recovers
    (expended upper stage, propellant, range fees, operations)."""
    return M / N + R + F


def cost_per_kg(cost_per_flight_musd, payload_kg):
    """Launch cost in dollars per kilogram."""
    return cost_per_flight_musd * 1e6 / payload_kg


def reuse_break_even(M, R, F):
    """Flights N at which reuse ties an expendable vehicle costing (M + F) per flight:
    solve M/N + R + F = M + F  ->  N = M / (M - R)."""
    return M / (M - R)

Additional utilities

Several chapters contributed code beyond the eleven canonical modules. The first group below are small, self-standing modules that compose with the package the same way — a launch-vehicle screen, a light-time calculator, reliability block math, a coverage-geometry helper, an end-of-life disposal module, and an orbit-regime lookup. The second group, standalone Mission-Design helpers, are tools the chapters explicitly flagged as not part of the importable flight package: they belong in your Mission Design Review folder, and are collected here so the appendix reproduces every line of code the project produced. (Headings for these use a role and chapter rather than a module filename, to keep the distinction honest.)

launch.py

A first-cut screen listing the Appendix-H vehicles whose approximate capability clears a payload plus margin to a destination. Built in Chapter 30; feeds the launch-vehicle-selection step of mission.py. Capacities are Tier-2 approximate — planning intuition, never a value to fly a mission on.

"""astrotools/launch.py -- first-cut launch-vehicle screen.

Built in Chapter 30. Capabilities are Tier-2 approximate figures from Appendix H.
"""

# Approximate capability in kg (Appendix H, Tier 2).
VEHICLES = {            # (LEO kg, GTO kg)
    "Electron":       (300,    0),
    "LVM3":           (10000,  4000),
    "Ariane 6":       (21600,  11500),
    "Falcon 9":       (22800,  8300),
    "Vulcan Centaur": (27000,  14000),
    "New Glenn":      (45000,  13000),
    "Falcon Heavy":   (63800,  26700),
}


def select_launcher(payload_kg, destination, margin=0.10):
    """Sorted names of vehicles whose capability clears payload*(1+margin).
    `destination` is 'LEO' or 'GTO'. A first screen, not a design tool."""
    col = {"LEO": 0, "GTO": 1}[destination]
    need = payload_kg * (1 + margin)
    return sorted(name for name, cap in VEHICLES.items() if cap[col] >= need)

operations.py

Turns a mission's distance into the round-trip light time that sets its autonomy level. Built in Chapter 31.

"""astrotools/operations.py -- light-time and autonomy level.

Built in Chapter 31.
"""

C_KM_S = 2.998e5    # speed of light, km/s
AU_KM = 1.496e8     # astronomical unit, km


def light_time(distance_km):
    """One-way light-travel time in seconds."""
    return distance_km / C_KM_S


def round_trip_light_time(distance_km):
    """Round-trip (command up + telemetry back) light time in seconds."""
    return 2 * light_time(distance_km)


def autonomy_level(distance_km):
    """A rough operations rule of thumb keyed to round-trip light time."""
    rtlt = round_trip_light_time(distance_km)
    if rtlt < 1:
        return "real-time control feasible"
    if rtlt < 10:
        return "supervised; time-tag critical events"
    return "autonomous; sequence-based commanding"

reliability.py

Reliability block math for the risk-assessment checkpoint: series and parallel combinations, redundant units, and the k-of-n survival probability behind engine-out capability. Built in Chapter 32; used in the Chapter 40 capstone.

"""astrotools/reliability.py -- reliability block math.

Built in Chapter 32; used in the Chapter 40 capstone risk assessment.
"""
from math import comb


def series(rs):
    """Reliability of independent components all required (in series): the product of R_i."""
    p = 1.0
    for r in rs:
        p *= r
    return p


def parallel(rs):
    """Reliability of independent redundant units (works if at least one works)."""
    q = 1.0
    for r in rs:
        q *= (1.0 - r)
    return 1.0 - q


def redundant(r, n):
    """n identical parallel units: probability that at least one of n survives = 1-(1-r)^n."""
    return 1.0 - (1.0 - r) ** n


def k_of_n(r, n, k):
    """Probability that at least k of n identical independent units survive (engine-out)."""
    return sum(comb(n, i) * r ** i * (1 - r) ** (n - i) for i in range(k, n + 1))

constellation.py

Coverage geometry for a satellite constellation: the Earth-central angle of a footprint, the fraction of the surface one satellite sees, and a floor on the number of satellites for instantaneous global coverage. Added in Chapter 33 — a convenience helper; orbits.py remains the workhorse.

"""astrotools/constellation.py -- coverage geometry (convenience helper).

Added in Chapter 33. A design convenience, not one of the canonical modules.
"""
import math

R_EARTH = 6371.0  # km, mean radius


def earth_central_angle(alt_km, min_elev_deg):
    """Earth-central angle (deg) of the ground footprint a satellite at altitude alt_km
    can serve down to a minimum elevation of min_elev_deg above the horizon."""
    r = R_EARTH + alt_km
    eps = math.radians(min_elev_deg)
    eta = math.asin((R_EARTH / r) * math.cos(eps))   # nadir angle
    lam = math.pi / 2 - eps - eta                    # Earth-central angle
    return math.degrees(lam)


def coverage_fraction(alt_km, min_elev_deg):
    """Fraction of Earth's surface inside one satellite's footprint (spherical cap)."""
    lam = math.radians(earth_central_angle(alt_km, min_elev_deg))
    return (1 - math.cos(lam)) / 2


def min_satellites_for_global(alt_km, min_elev_deg):
    """Floor on satellites for instantaneous global coverage (ignores tiling overlap; the
    real number is several times larger, and larger still for capacity)."""
    return math.ceil(1 / coverage_fraction(alt_km, min_elev_deg))

debris.py

End-of-life disposal and sustainability: collision probability from a debris flux, the natural-decay lifetime of a circular orbit, and the delta-v to raise to a graveyard. Added in Chapter 35. The atmospheric density rho0 is solar-cycle-dependent — a modeling caveat the chapter flags.

"""astrotools/debris.py -- end-of-life disposal and sustainability.

Added in Chapter 35. rho0 is solar-cycle-dependent (a flagged modeling caveat).
"""
import math

MU_EARTH = 3.986e14     # m^3/s^2
RE_EARTH = 6.371e6      # m (mean radius)
SEC_PER_YEAR = 3.156e7


def collision_probability(flux, area_m2, years):
    """P(>=1 impact) for a debris flux (impacts/m^2/yr), area (m^2), and exposure (yr)."""
    n = flux * area_m2 * years
    return 1.0 - math.exp(-n)


def orbit_lifetime(alt_km, beta=100.0, rho0=2e-13, H_km=70.0):
    """Natural circular-orbit decay lifetime (years): tau ~ H*beta/(rho0*a*v)."""
    a = RE_EARTH + alt_km * 1e3
    v = math.sqrt(MU_EARTH / a)
    tau_s = (H_km * 1e3) * beta / (rho0 * a * v)
    return tau_s / SEC_PER_YEAR


def reorbit_dv(v, dh_km, a_km):
    """Two-burn delta-v (m/s) to raise a circular orbit by dh (small-change approx)."""
    return 0.5 * v * (dh_km / a_km)

orbit_catalog.py

A lookup table of the standard Earth-orbit regimes with a circular-period calculator, so you can sanity-check any altitude you consider. Added in Chapter 9 — a convenience helper; orbits.py (Chapter 6) stays the real workhorse.

"""astrotools/orbit_catalog.py -- Earth-orbit regime lookup (convenience helper).

Added in Chapter 9. A lookup table with a period calculator; not a canonical module.
"""
import math

MU_EARTH = 3.986e5      # km^3/s^2
R_EARTH = 6371.0        # km (mean radius; use 6378 km equatorial for the exact GEO altitude)

# regime: (typical altitude km, typical inclination deg, one-line use)
CATALOG = {
    "LEO": (400,   "28-98", "imaging, ISS, broadband; cheapest to reach"),
    "SSO": (800,   "98",    "Earth observation at constant local lighting"),
    "MEO": (20200, "55",    "navigation: GPS / Galileo / GLONASS"),
    "GEO": (35786, "0",     "comms & weather; fixed over one longitude"),
}


def circular_period_min(alt_km):
    """Period (minutes) of a circular Earth orbit at altitude alt_km (km)."""
    a = R_EARTH + alt_km
    return 2 * math.pi * math.sqrt(a**3 / MU_EARTH) / 60.0


def classify(alt_km):
    """Rough regime name for a circular orbit at altitude alt_km (km)."""
    if alt_km < 2000:
        return "LEO"
    if abs(alt_km - 35786) < 50:
        return "GEO"
    if alt_km < 35786:
        return "MEO"
    return "above GEO"

Standalone Mission-Design helpers

These come from chapters off the module schedule and were flagged in the text as teaching aids for your Mission Design Review, not part of the importable flight package. They are reproduced here for completeness; each is labeled with the chapter that wrote it.

Ascent aerodynamics (Chapter 5)

Dynamic pressure, max-Q, and the constant-q speed cap for a structural limit — the inputs to the MDR "Ascent Loads and Environment" note.

"""Ascent aerodynamics helper (Chapter 5) -- dynamic pressure and max-Q.

Standalone MDR helper, off the canonical module schedule.
"""
import math

RHO0, H = 1.225, 8000.0   # sea-level density (kg/m^3), scale height (m)


def density(h):
    """Exponential-atmosphere density (kg/m^3) at altitude h (m)."""
    return RHO0 * math.exp(-h / H)


def dynamic_pressure(h, v):
    """Dynamic pressure q (Pa) at altitude h (m), speed v (m/s)."""
    return 0.5 * density(h) * v**2


def find_max_q(profile):
    """Given [(h_m, v_ms), ...], return (q_max_Pa, h_at_max_m)."""
    best = max(profile, key=lambda hv: dynamic_pressure(*hv))
    return dynamic_pressure(*best), best[0]


def v_cap(h, q_lim):
    """Max speed (m/s) at altitude h (m) keeping q at the limit q_lim (Pa)."""
    return math.sqrt(2 * q_lim / density(h))

Ballistic re-entry (Chapter 7)

Entry kinetic energy, the Allen-Eggers peak-g (and its inversion to an entry angle), and the plasma frequency that sets radio blackout — the MDR entry / disposal note.

"""Ballistic re-entry helper (Chapter 7) -- Allen-Eggers peak-g and blackout.

Standalone MDR helper; re-entry has no canonical module in the toolkit.
"""
import math


def specific_ke(v):
    """Kinetic energy per kilogram (J/kg) at entry speed v (m/s)."""
    return 0.5 * v**2


def peak_g(v_entry, gamma_deg, scale_height=7200.0):
    """Peak deceleration (Earth g's) of a ballistic entry (Allen-Eggers), independent of
    ballistic coefficient. gamma_deg is the entry flight-path angle below horizontal."""
    a_max = v_entry**2 * math.sin(math.radians(gamma_deg)) / (2 * math.e * scale_height)
    return a_max / 9.81


def entry_angle_for_g(g_limit, v_entry, scale_height=7200.0):
    """Ballistic entry angle (deg) yielding a target peak-g limit (Allen-Eggers, inverted)."""
    sin_g = g_limit * 9.81 * (2 * math.e * scale_height) / v_entry**2
    return math.degrees(math.asin(sin_g))


def plasma_frequency_ghz(n_e):
    """Plasma frequency (GHz) for electron density n_e (electrons/m^3); radio links below
    this frequency are blacked out."""
    return 8.98 * math.sqrt(n_e) / 1e9

Navigation and control: scalar Kalman update and PID (Chapters 13, 27)

A scalar Kalman update that fuses a prediction with one measurement (Chapter 13) and a minimal discrete PID controller (Chapter 27) — together they sketch a navigate-then-control loop for the GN&C note. The flight versions are a full extended Kalman filter and attitude controller, beyond this book's linear-algebra level.

"""Navigation/control sketches (Chapters 13 and 27) -- scalar Kalman + PID.

Standalone teaching aids, NOT part of the stable flight API.
"""


def kalman_update(x_pred, var_pred, z, var_meas):
    """Optimally fuse a prediction with one measurement of the same quantity (Chapter 13).
    Returns (updated estimate, updated variance)."""
    K = var_pred / (var_pred + var_meas)     # gain: relative trust in the measurement
    x_upd = x_pred + K * (z - x_pred)         # move from prediction toward measurement
    var_upd = (1 - K) * var_pred             # information added -> uncertainty shrinks
    return x_upd, var_upd


class PID:
    """A minimal discrete PID controller (Chapter 27): u = Kp*e + Ki*sum(e*dt) + Kd*de/dt."""
    def __init__(self, kp, ki, kd, dt):
        self.kp, self.ki, self.kd, self.dt = kp, ki, kd, dt
        self.integral, self.prev_error = 0.0, 0.0

    def step(self, error):
        self.integral += error * self.dt
        derivative = (error - self.prev_error) / self.dt
        self.prev_error = error
        return self.kp * error + self.ki * self.integral + self.kd * derivative

Crew consumables / life support (Chapter 28)

Open-loop consumable mass and the mission day at which a closed-loop recycler undercuts it — for a crewed mission or a crewed variant of an otherwise-uncrewed track.

"""Crew-consumables helper (Chapter 28) -- open- vs closed-loop life support.

Standalone helper; life support has no canonical module in the flight toolkit.
"""


def crew_consumables(crew, days, o2=0.84, water=2.5, food=1.8):
    """Open-loop supplied mass (kg) for a crew over a mission, by commodity (kg/person-day
    defaults). Returns a dict of o2, water, food, and total."""
    per_day = o2 + water + food                 # kg per person-day
    return {"o2": crew * days * o2, "water": crew * days * water,
            "food": crew * days * food, "total": crew * days * per_day}


def recycle_breakeven(hw_mass, crew, open_rate=5.14, closed_rate=2.47):
    """Mission day at which closed-loop mass undercuts open-loop, for a fixed crew."""
    return hw_mass / (crew * (open_rate - closed_rate))

Launch economics: flight rate and price per kilogram (Chapters 37, 39)

Two economics one-liners: the cost per flight when a fixed program cost is amortized over a flight rate (Chapter 37), and the launch cost of a given mass at a price per kilogram (Chapter 39). The reuse-specific cost model lives in mission.py (Chapter 38).

"""Launch-economics helpers (Chapters 37 and 39).

Standalone Tier-2/Tier-3 illustrative helpers. Reuse economics live in mission.py.
"""


def cost_per_flight(marginal, fixed_annual, flights_per_year):
    """Cost per flight = marginal cost + fixed annual program cost / flight rate
    (Chapter 37). Consistent units, e.g. $B."""
    return marginal + fixed_annual / flights_per_year


def launch_cost(wet_mass_kg, dollars_per_kg):
    """Cost to launch wet_mass_kg at a given price per kg to LEO (Chapter 39)."""
    return wet_mass_kg * dollars_per_kg

Mission-selection warm-up (Chapter 1)

The very first checkpoint: record the chosen track and a rough, order-of-magnitude delta-v beyond low Earth orbit, read off the delta-v map. Every number here is refined by later modules; the package proper begins with rocket.py in Chapter 3.

"""Mission-selection warm-up (Chapter 1) -- choose and record a mission.

The seed of the Design Your Mission project; the package proper begins in Chapter 3.
"""

DV_TO_LEO = 9.4  # km/s, the admission price of orbit

# track: (name, rough delta-v beyond LEO in km/s -- order of magnitude only)
MISSIONS = {
    "A": ("GEO communications satellite", 4.3),
    "B": ("Lunar lander (cargo)",         6.0),
    "C": ("Mars orbiter (science)",       5.3),
    "D": ("Asteroid rendezvous",          5.5),
}


def summarize(track):
    """Return (name, delta-v beyond LEO, rough total from the ground) for a track."""
    name, dv_beyond_leo = MISSIONS[track]
    return name, dv_beyond_leo, DV_TO_LEO + dv_beyond_leo

A closing note

This is a teaching toolkit, not production software. Its virtue is that you can read every line and see exactly which physics it encodes; its limitation is that it trades the robustness, edge-case handling, and current-ephemeris precision of a real astrodynamics library for clarity. For mission-grade work, reach past these round numbers to a maintained package — poliastro, astropy, GMAT, or an agency's own tools — exactly as Appendix B and Appendix C advise. But for building intuition, for checking a back-of-the-envelope number, and for watching your own Mission Design Review close, astrotools is enough. To run it, drop these modules into a folder alongside a requirements.txt and install the four libraries in one line:

pip install -r requirements.txt

Then point the same functions that sized the Track-A comsat at your mission's numbers, and let the package size your vehicle too. The symbols and units the code uses are cataloged in Appendix A; the constants it hard-codes are tabulated, at full precision, in Appendix B. Now go run something, break it, and fix it. That is how rockets get built.