Running the Simulation

Running the Simulation

All functionality is in the orbital_thermal package. The modules you will use most:

Module	Purpose
equilibrium	Steady radiator temperature and fixed-temperature capacity
radiation	Gray-body net flux, required area
bounds	The analytic thermodynamic bounds (Theorems 1–5)
environment	Orbital period/radius and the exact tilted-plate-to-sphere view factor
sink	Orbit-varying effective sink and the subpoint-albedo model
transient	One-node RK4 radiator solver, convergence certificate, averaging-load bias
mccalip_replication / mccalip_exact_vf	Replication of McCalip's model and the exact-view-factor correction
fluids	CoolProp-backed ammonia coolant screen (optional)

Required flag: the orbit-coupled functions (sink.*, transient.simulate, transient.averaging_bias) take a mandatory keyword assume_sun_shielded=True. It is the single point where the model asserts that direct solar load on the panel face is omitted; there is no default, so you must pass it explicitly.

---

1. Static radiator sizing

from orbital_thermal import equilibrium, radiation

# Equilibrium temperature to reject 120 kW through 220 m^2 at emissivity 0.91, sink 220 K
T = equilibrium.equilibrium_temperature(Q=120e3, area=220.0, emissivity=0.91, T_sink=220.0)
print(T)            # 337.1 K  (the AI1 primary operating point)

# Inverse: heat rejection capacity (W) of a panel held at a temperature
Q = equilibrium.radiative_capacity(T=337.1, area=220.0, emissivity=0.91, T_sink=220.0)

---

2. The McCalip exact-view-factor correction (headline result)

from orbital_thermal import mccalip_exact_vf as mx

# Correction table across orbit beta angle (exact per-face VF vs McCalip's heuristic)
for row in mx.correction_table_vs_beta():
    print(row["beta_deg"], round(row["delta_K"], 2))
# ... 90 6.35   <- +6.35 K at the edge-on default

# Recompute McCalip's default equilibrium with the exact view factor
print(round(mx.eqtemp_exact_vf({}), 4))     # 342.0992 K  (vs his coded 335.7495 K)

The exact tilted-plate-to-sphere Earth view factor itself:

from orbital_thermal import environment as env
print(round(env.sphere_view_factor(altitude_km=550, tilt_deg=90.0), 6))   # 0.257773 (edge-on, per face)

---

3. Orbit-coupled effective sink

from orbital_thermal import sink

# Radiatively-weighted orbit-average sink (K) for a nadir-facing panel
print(round(sink.orbit_averaged_sink(550, 0.0, tilt_deg=0, assume_sun_shielded=True), 2))   # 250.99 K

# The full effective-sink profile around one orbit (u in degrees, T_s_eff in K)
u, Ts = sink.sink_profile(550, 30.0, tilt_deg=0, assume_sun_shielded=True)

---

4. Transient simulation and the convergence certificate

The one-node model C dT/dt = q_load − εσ(T⁴ − T_sink_eff(t)⁴) is marched with RK4 to a periodic steady state. Use simulate(...); pass return_diagnostics=True and check_time_resolution=True to get the full certificate.

from orbital_thermal import transient
from orbital_thermal.constants import SIGMA_SB

q_load = 0.91 * SIGMA_SB * (337.1**4 - 220.0**4)     # ~545.5 W/m^2

t, T, T_sink, diag = transient.simulate(
    altitude_km=550, beta_deg=30, q_load=q_load,
    areal_heat_capacity=8000.0,       # J/m^2/K
    tilt_deg=0, assume_sun_shielded=True,
    return_diagnostics=True, check_time_resolution=True,
)
print(diag["converged"], diag["time_residual_K"])    # True  (certified)

diag reports three independent convergence axes — periodic_converged (closure + scale-aware energy balance), time_discretization_converged (N/2N/4N grid refinement + forcing certificate), and the combined converged — plus component residuals (forcing_residual_K, pointwise_n_to_4n_K, peak_phase_residual_deg, …) and refined_orbits_used.

Averaging-load bias

averaging_bias(...) runs simulate and compares the transient mean to the steady, averaged-sink solution. It raises RuntimeError unless the result is fully certified (so you never get a bias from an unconverged orbit). Pass require_convergence=False to inspect an uncertified run.

b = transient.averaging_bias(550, 30, q_load, 8000.0, tilt_deg=0, assume_sun_shielded=True)
print(round(b["transient_mean_K"], 2),           # 346.89 K
      round(b["peak_excess_over_steady_K"], 2),   # 3.01 K  (the operational penalty)
      round(b["bias_K"], 3))                      # -0.014 K (mean <= steady, by Jensen)

Knobs: n_orbits / max_orbits (how long to march), steps_per_orbit (temporal resolution), time_safety_factor (default 2.0, the conservative margin on the temporal estimate), convergence_tol_K, energy_tol_K, time_tol_K.

---

5. Generating the figures

The plotting scripts require matplotlib, which is included in the dev extra:

python -m pip install -e ".[dev]"

The figures in the papers are reproduced by scripts under scripts/ and should be run from the repository root:

python scripts/plot_effective_sink.py      # effective sink around the orbit
python scripts/plot_mccalip_correction.py  # the +6.35 K correction vs beta
python scripts/plot_transient.py           # transient temperature ripple
python scripts/plot_edge_on_geometry.py    # the edge-on geometry schematic (paper three)

Each writes a PNG into results/figures/.