"""Neutron mission, in readable Python: five states, one ODE, explicit events. Educational planar translation model; NOT Rocket Lab flight software or a performance prediction. Ideal pointing, nonrotating spherical Earth, SI units. SciPy RK45 is the adaptive Dormand–Prince embedded fifth/fourth-order method. """ from dataclasses import asdict, dataclass import argparse import json import math import numpy as np from scipy.integrate import solve_ivp MU = 3.986004418e14 # Earth gravitational parameter [m^3/s^2] EARTH_RADIUS = 6_378_137.0 # spherical Earth radius [m] G0 = 9.80665 # Isp reference gravity [m/s^2] SOURCE = "https://rocketlabcorp.com/launch/neutron/" TRANSFER_ALTITUDE_KM = 220.0 # chosen mission design, not a physical constant RECOVERY_DECK_HEIGHT_M = 8.0 # original schematic vessel geometry assumption RECOVERY_SHIP_LENGTH_M = 122.0 RECOVERY_SHIP_WIDTH_M = 50.0 # assumed beam; planar flight has zero cross-track RECOVERY_LEG_CLEARANCE_M = 1.6588194256045004 # authored deployed foot below engine plane @dataclass(frozen=True) class Vehicle: """Published geometry/thrust; explicitly assumed internal mass and Isp.""" height_m: float = 43.0 diameter_m: float = 7.0 published_liftoff_mass_kg: float = 480_000.0 advertised_leo_payload_kg: float = 13_000.0 booster_engines: int = 9 booster_thrust_n: float = 1_450_000.0 * 4.4482216152605 upper_thrust_n: float = 890_000.0 booster_dry_kg: float = 30_000.0 # assumption, includes captive fairing booster_propellant_kg: float = 335_000.0 booster_isp_s: float = 330.0 upper_dry_kg: float = 5_000.0 upper_propellant_kg: float = 102_000.0 upper_isp_s: float = 365.0 drag_coefficient: float = 0.35 def atmosphere(altitude_m): """One-layer density model, not a weather or aerothermal model.""" return 1.225 * math.exp(-max(0.0, altitude_m) / 8500.0) def orbital_elements(state): """Osculating two-body apsides; a bound ellipse alone is not a safe orbit.""" radius, _, vr, vt, _ = map(float, state) energy = (vr * vr + vt * vt) / 2 - MU / radius momentum = radius * vt # Polar eccentricity-vector components avoid cancellation near a circle. eccentricity = math.hypot(radius * vt * vt / MU - 1, radius * vr * vt / MU) p = momentum**2 / MU return {"perigee_km": (p / (1 + eccentricity) - EARTH_RADIUS) / 1000, "apogee_km": ((p / (1 - eccentricity) - EARTH_RADIUS) / 1000 if energy < 0 and eccentricity < 1 else None), "eccentricity": eccentricity, "specific_energy_j_kg": energy} def orbit_target_error(state, target_altitude_km): """Enter the target set only when BOTH apsides and radial speed agree. A proper target orbit is bound, within 2 km at each apsis, |vr| <= 5 m/s, and e <= 0.001. A negative result satisfies all four constraints. """ orbit = orbital_elements(state) if orbit["apogee_km"] is None or orbit["specific_energy_j_kg"] >= 0: return 1e6 return max(abs(orbit["perigee_km"] - target_altitude_km) / 2, abs(orbit["apogee_km"] - target_altitude_km) / 2, abs(state[2]) / 5, orbit["eccentricity"] / 0.001) - 1 def coast_verification(solution, period_s, target_altitude_km): """Verify an actual propagated revolution; no analytic orbit is substituted.""" radius, angle, vr, vt, _ = solution.y altitude_km = (radius - EARTH_RADIUS) / 1000 energy = (vr**2 + vt**2) / 2 - MU / radius momentum = radius * vt energy_change = float(np.max(np.abs((energy - energy[0]) / energy[0]))) momentum_change = float(np.max(np.abs((momentum - momentum[0]) / momentum[0]))) duration = float(solution.t[-1] - solution.t[0]) swept_angle = float(angle[-1] - angle[0]) completed = abs(duration - period_s) < 1e-3 and not solution.t_events[0].size minimum, maximum = float(altitude_km.min()), float(altitude_km.max()) initial_orbit = orbital_elements(solution.y[:, 0]) success = (completed and np.all(energy < 0) and swept_angle >= 2 * math.pi - 1e-3 and minimum >= target_altitude_km - 2 and maximum <= target_altitude_km + 2 and energy_change < 1e-5 and momentum_change < 1e-5) return dict(success=bool(success), completed_one_orbit=bool(completed), period_s=period_s, propagated_duration_s=duration, swept_angle_rad=swept_angle, perigee_km=initial_orbit["perigee_km"], apogee_km=initial_orbit["apogee_km"], eccentricity=initial_orbit["eccentricity"], min_altitude_km=minimum, max_altitude_km=maximum, relative_energy_change=energy_change, relative_angular_momentum_change=momentum_change) def drag_force(state, area_m2, cd): """Components in the local outward/tangential frame; still atmosphere.""" radius, _, vr, vt, _ = state speed = math.hypot(vr, vt) scale = -0.5 * atmosphere(radius - EARTH_RADIUS) * cd * area_m2 * speed return np.array([scale * vr, scale * vt]) def rhs(time, state, guidance, area_m2, cd, isp_s): """Polar Newton equations; state = [r, theta, v_r, v_t, mass]. Pointing is ideal: guidance commands force directly, not attitude/torque. The changing mass appears in F/m and dm/dt; thrust includes exhaust momentum. """ radius, _, vr, vt, mass = state force = guidance(time, state) radial, tangential = (force + drag_force(state, area_m2, cd)) / mass return [vr, vt / radius, vt * vt / radius - MU / radius**2 + radial, -vr * vt / radius + tangential, -np.linalg.norm(force) / (isp_s * G0)] def terminal_event(function, direction=-1): """SciPy locates a zero crossing between accepted RK45 steps.""" function.terminal = True function.direction = direction return function def integrate(state, start_s, end_s, guidance, area_m2, cd, isp_s, events=(), max_step_s=2.0, rtol=1e-7): """Integrate one continuous phase; restart after every discrete change.""" solution = solve_ivp( lambda t, y: rhs(t, y, guidance, area_m2, cd, isp_s), (start_s, end_s), state, method="RK45", dense_output=True, events=events, max_step=max_step_s, rtol=rtol, atol=[1e-3, 1e-11, 1e-5, 1e-5, 1e-4]) if not solution.success: raise RuntimeError(solution.message) return solution def ascent_guidance(time, state, vehicle): """An assumed pitch program; angles measured above local horizontal.""" pitch = np.interp(time, [0, 12, 35, 70, 115, 160], [90, 90, 82, 67, 52, 40]) beta = math.radians(pitch) return vehicle.booster_thrust_n * np.array([math.sin(beta), math.cos(beta)]) def upper_guidance(time, state, vehicle, target_m): """PD radial acceleration request plus spherical-gravity compensation.""" radius, _, vr, vt, mass = state tau = 30.0 requested = (EARTH_RADIUS + target_m - radius) / tau**2 - 2 * vr / tau radial_fraction = (requested + MU / radius**2 - vt * vt / radius) * mass / vehicle.upper_thrust_n beta = math.asin(np.clip(radial_fraction, -0.5, 0.95)) return vehicle.upper_thrust_n * np.array([math.sin(beta), math.cos(beta)]) def predicted_miss(state, target_downrange_m=0.0, contact_height_m=0.0): """Cheap constant-g ballistic predictor used only to stop boostback. It does not set a landing position; the integrated trajectory can miss. Landing feedback subsequently corrects remaining position/velocity errors. """ radius, angle, vr, vt, _ = state gravity = MU / radius**2 fall_time = (vr + math.sqrt(vr * vr + 2 * gravity * max(0, radius - EARTH_RADIUS - contact_height_m))) / gravity return EARTH_RADIUS * angle + vt * fall_time - target_downrange_m def landing_guidance(time, state, vehicle, target_downrange_m=0.0, contact_height_m=0.0): """Soft-descent velocity law + finite-time lateral position feedback. Ideal engine envelope: 0.3 of one engine to three engines. Engine switching, ignition transients, attitude slew and gimbal limits are not simulated. """ radius, angle, vr, vt, mass = state height = max(0.0, radius - EARTH_RADIUS - contact_height_m) desired_vr = -math.sqrt(2 * 6.0 * height + 1.0**2) radial_request = -6.0 * vr / -desired_vr + (desired_vr - vr) / 2.0 remaining = max(3.0, 2 * height / max(1.0, -vr)) lateral_request = -6 * (EARTH_RADIUS * angle - target_downrange_m) / remaining**2 - 4 * vt / remaining acceleration = np.array([radial_request + MU / radius**2 - vt * vt / radius, lateral_request + vr * vt / radius]) area = math.pi * vehicle.diameter_m**2 / 4 force = mass * acceleration - drag_force(state, area, vehicle.drag_coefficient) # Only upward pointing during the landing burn; bounded total force. force[0] = max(0.0, force[0]) magnitude = np.linalg.norm(force) engine = vehicle.booster_thrust_n / vehicle.booster_engines return force * np.clip(magnitude, 0.3 * engine, 3 * engine) / max(magnitude, 1e-9) def separate(state, booster_mass_kg, relative_speed_m_s=0.5): """Split mass with equal/opposite impulses: conserve linear momentum. Both bodies start at the same position. Geometry and contact are not modeled. The reusable fairing stays within booster dry mass. """ booster, upper = state.copy(), state.copy() upper[4] = state[4] - booster_mass_kg booster[4] = booster_mass_kg booster[2] -= relative_speed_m_s * upper[4] / state[4] upper[2] += relative_speed_m_s * booster[4] / state[4] return booster, upper COAST_POINTING_RATE_RAD_S = math.radians(60.0) UPPER_EXIT_CLEARANCE_M = 43.0 - 24.8 + 5.0 # authored enclosure, upper base and margin def apply_kinematic_hardware(trajectories, events, dense_phases): """Prescribed hardware coordinates, sampled from the actual integrated flight. Coast pointing uses a bounded cubic pre-alignment before the next ignition. Powered pointing is the actual force axis. Fairing closure waits for physical upper-stage clearance. No position, velocity, mass or applied force changes. """ from scipy.optimize import brentq def phase_at(body, time): return next((phase for phase in reversed(dense_phases[body]) if phase["start"] <= time <= phase["stop"]), None) def state_xy(body, time): state = phase_at(body, time)["solution"].sol(time) return state[0] * np.array([math.cos(state[1]), math.sin(state[1])]) def insert_samples(body, times): rows = trajectories[body] for time in times: if any(abs(row["time_s"] - time) < 1e-8 for row in rows): continue phase = phase_at(body, time) if phase is None: continue state = phase["solution"].sol(time) radius, angle, vr, vt, mass = map(float, state) force = phase["guidance"](time, state) thrust = float(np.linalg.norm(force)) ca, sa = math.cos(angle), math.sin(angle) fx, fy = float(force[0] * ca - force[1] * sa), float(force[0] * sa + force[1] * ca) before = max((row for row in rows if row["time_s"] <= time), key=lambda row: row["time_s"]) row = dict(before, time_s=float(time), phase=phase["name"], altitude_m=radius - EARTH_RADIUS, downrange_m=EARTH_RADIUS * angle, radial_velocity_m_s=vr, tangential_velocity_m_s=vt, speed_m_s=math.hypot(vr, vt), mass_kg=mass, thrust_n=thrust, force_x_n=fx, force_y_n=fy, dynamic_pressure_pa=0.5 * atmosphere(radius - EARTH_RADIUS) * (vr * vr + vt * vt), x_m=radius * ca, y_m=radius * sa, vx_m_s=vr * ca - vt * sa, vy_m_s=vr * sa + vt * ca) if thrust > 0: row["axis_x"], row["axis_y"] = fx / thrust, fy / thrust if body == "booster" and any(sample.get("leg_open_fraction", 0) > 0 for sample in rows): command = next((event["time_s"] for event in events if event["name"] == "landing_burn"), math.inf) row["leg_open_fraction"] = float(np.clip((time - command) / 2, 0, 1)) rows.append(row) rows.sort(key=lambda row: row["time_s"]) for body, rows in trajectories.items(): segments = [] index = 0 while index < len(rows): if rows[index]["thrust_n"] > 0: index += 1 continue start = index while index < len(rows) and rows[index]["thrust_n"] == 0: index += 1 if index == len(rows): break first, target = rows[start], rows[index] initial = math.atan2(first["axis_y"], first["axis_x"]) final = math.atan2(target["force_y_n"], target["force_x_n"]) turn = math.atan2(math.sin(final - initial), math.cos(final - initial)) available = target["time_s"] - first["time_s"] duration = min(available, max(2.0, 1.5 * abs(turn) / COAST_POINTING_RATE_RAD_S)) if duration > 0: segments.append((target["time_s"] - duration, target["time_s"], initial, turn)) # Include the profile endpoints and 0.1 s kinematic samples even when # the caller chooses sparse trajectory output. States use RK45 dense # output here; they are never linearly invented or re-integrated. insert_samples(body, [float(t) for begin, end, _, _ in segments for t in np.r_[np.arange(begin, end, 0.1), end]]) for begin, end, initial, turn in segments: for row in rows: if row["thrust_n"] == 0 and begin <= row["time_s"] <= end: u = float(np.clip((row["time_s"] - begin) / (end - begin), 0, 1)) angle = initial + turn * u * u * (3 - 2 * u) row["axis_x"], row["axis_y"] = math.cos(angle), math.sin(angle) opening = next((event["time_s"] for event in events if event["name"] == "fairing_opening"), None) separation = next((event["time_s"] for event in events if event["name"] == "stage_separation"), None) if opening is None or separation is None or not trajectories["booster"] or not trajectories["upper_stage"]: return attached = trajectories["stack"][-1] axis = np.array([attached["axis_x"], attached["axis_y"]]) def clearance(time): return float((state_xy("upper_stage", time) - state_xy("booster", time)) @ axis) - UPPER_EXIT_CLEARANCE_M stop = min(dense_phases[body][-1]["stop"] for body in ("booster", "upper_stage")) times = sorted(set(float(t) for body in ("booster", "upper_stage") for phase in dense_phases[body] for t in phase["solution"].t if separation <= t <= stop)) close = math.inf for left, right in zip(times, times[1:]): if clearance(left) <= 0 <= clearance(right): close = float(brentq(clearance, left, right, xtol=1e-9)) break events[:] = [event for event in events if event["name"] not in ("fairing_closing", "fairing_closed")] if math.isfinite(close): for name, time, text in [ ("fairing_closing", close, "Captive fairing closes after integrated upper-stage clearance."), ("fairing_closed", close + 4, "Four-second closure completes after upper-stage clearance.")]: events.append(dict(name=name, time_s=time, body="booster", description=text, required_axial_clearance_m=UPPER_EXIT_CLEARANCE_M)) insert_samples("booster", [close, close + 4]) for body in ("stack", "booster"): for row in trajectories[body]: time = row["time_s"] fraction = min(1.0, max(0.0, (time - opening) / 4)) if time > close: fraction = max(0.0, 1 - (time - close) / 4) row["fairing_open_fraction"] = fraction def simulate(payload_kg=8000.0, recovery_propellant_kg=60_000.0, target_altitude_km=400.0, upper_propellant_scale=1.0, max_step_s=2.0, rtol=1e-7, sample_step_s=2.0, recovery_target_downrange_m=0.0): """Run ascent, deployment and independent booster recovery; return JSON data. Inputs change actual integrated states. Insufficient fuel, missed orbit and hard/remote touchdowns remain failures. No endpoint is moved to the target. """ vehicle = Vehicle() values = [payload_kg, recovery_propellant_kg, target_altitude_km, upper_propellant_scale, max_step_s, rtol, sample_step_s, recovery_target_downrange_m] if not all(math.isfinite(value) for value in values): raise ValueError("All inputs must be finite.") if payload_kg <= 0 or not 0 <= recovery_propellant_kg < vehicle.booster_propellant_kg: raise ValueError("Payload must be positive; reserve must be inside stage-one fuel capacity.") if not 150 <= target_altitude_km <= 1000 or not 0 < upper_propellant_scale <= 1.5: raise ValueError("Target must be 150–1000 km; upper-stage fuel scale must be (0, 1.5].") if min(max_step_s, rtol, sample_step_s) <= 0: raise ValueError("Solver and sampling controls must be positive.") if not 0 <= recovery_target_downrange_m <= 2_000_000: raise ValueError("Recovery target must be between 0 and 2,000,000 m downrange.") ship_recovery = recovery_target_downrange_m > 0 target_angle = recovery_target_downrange_m / EARTH_RADIUS deck_height = RECOVERY_DECK_HEIGHT_M if ship_recovery else 0.0 leg_clearance = RECOVERY_LEG_CLEARANCE_M if ship_recovery else 0.0 contact_height = deck_height + leg_clearance booster_contact_surface = None area = math.pi * vehicle.diameter_m**2 / 4 upper_wet = vehicle.upper_dry_kg + vehicle.upper_propellant_kg * upper_propellant_scale + payload_kg initial_mass = vehicle.booster_dry_kg + vehicle.booster_propellant_kg + upper_wet initial = np.array([EARTH_RADIUS, 0.0, 0.0, 0.0, initial_mass]) trajectories = {body: [] for body in ("stack", "booster", "upper_stage", "payload")} dense_phases = {body: [] for body in trajectories} events = [] integrations = 0 separation_time = None # Display pointing is the model's ideal thrust axis, not an attitude plant. # Separated bodies inherit their parent's axis. The explicit kinematic # hardware pass below pre-aligns coast segments before subsequent ignitions. pointing = {body: (1.0, 0.0) for body in trajectories} def record(name, time, body, description, **details): events.append(dict(name=name, time_s=float(time), body=body, description=description, **details)) def fairing_fraction(time, body): opening = next((e["time_s"] for e in events if e["name"] == "fairing_opening"), None) if opening is None or body not in ("stack", "booster"): return 0.0 if separation_time is None or time <= separation_time: return float(np.clip((time - opening) / 4, 0, 1)) return float(np.clip(1 - (time - separation_time) / 4, 0, 1)) def leg_fraction(time, body): command = next((e["time_s"] for e in events if e["name"] == "landing_burn"), None) if not ship_recovery or body != "booster" or command is None: return 0.0 return float(np.clip((time - command) / 2, 0, 1)) def phase(body, name, state, start, end, guidance, monitors=(), body_area=None, sample_interval=None): nonlocal integrations, booster_contact_surface body_area = area if body_area is None else body_area isp = vehicle.upper_isp_s if body in ("upper_stage", "payload") else vehicle.booster_isp_s impact = terminal_event(lambda t, y: y[0] - EARTH_RADIUS) deck_events = () if body == "booster" and ship_recovery: # A fixed flat deck, not an elevated spherical surface everywhere. # Check its finite along-track footprint at actual downward roots. def deck_crossing(t, y): return (y[0] * math.cos(y[1] - target_angle) - EARTH_RADIUS - deck_height - leg_clearance * leg_fraction(t, body)) deck_crossing.direction = -1 deck_events = (deck_crossing,) solution = integrate(state, start, end, guidance, body_area, vehicle.drag_coefficient, isp, (impact, *monitors, *deck_events), max_step_s, rtol) integrations += solution.nfev stop = float(solution.t[-1]) hit = bool(solution.t_events[0].size) if hit and body == "booster": booster_contact_surface = "ocean" if ship_recovery else "ground" if deck_events: for contact_time in solution.t_events[-1]: contact_state = solution.sol(contact_time) along_track = contact_state[0] * math.sin(contact_state[1] - target_angle) if abs(along_track) <= RECOVERY_SHIP_LENGTH_M / 2: stop, hit = float(contact_time), True booster_contact_surface = "deck" solution.t_events = [times[times <= stop + 1e-9] for times in solution.t_events] break dense_phases[body].append({"start": float(start), "stop": stop, "solution": solution, "guidance": guidance, "name": name}) samples = np.unique(np.r_[start, np.arange(start, stop, sample_interval or sample_step_s), stop]) for time, y in zip(samples, solution.sol(samples).T): radius, angle, vr, vt, mass = y force = guidance(time, y) cos_a, sin_a = math.cos(angle), math.sin(angle) speed = math.hypot(vr, vt) thrust = float(np.linalg.norm(force)) force_x = float(force[0] * cos_a - force[1] * sin_a) force_y = float(force[0] * sin_a + force[1] * cos_a) if thrust > 0: pointing[body] = (force_x / thrust, force_y / thrust) trajectories[body].append(dict( time_s=float(time), phase=name, altitude_m=float(radius - EARTH_RADIUS), downrange_m=float(EARTH_RADIUS * angle), radial_velocity_m_s=float(vr), tangential_velocity_m_s=float(vt), speed_m_s=speed, mass_kg=float(mass), thrust_n=thrust, force_x_n=force_x, force_y_n=force_y, axis_x=pointing[body][0], axis_y=pointing[body][1], fairing_open_fraction=fairing_fraction(time, body), leg_open_fraction=leg_fraction(time, body), dynamic_pressure_pa=0.5 * atmosphere(radius - EARTH_RADIUS) * speed**2, x_m=radius * cos_a, y_m=radius * sin_a, vx_m_s=vr * cos_a - vt * sin_a, vy_m_s=vr * sin_a + vt * cos_a)) final_state = solution.y[:, -1] if stop == float(solution.t[-1]) else solution.sol(stop) return stop, final_state.copy(), hit, solution coast = lambda t, y: np.zeros(2) reserve_mass = vehicle.booster_dry_kg + recovery_propellant_kg meco = terminal_event(lambda t, y: y[4] - upper_wet - reserve_mass) record("liftoff", 0, "stack", "Nine-engine ascent begins.") time, state, impacted, ascent_solution = phase("stack", "ascent", initial, 0, 600, lambda t, y: ascent_guidance(t, y, vehicle), (meco,)) ascent_complete = bool(ascent_solution.t_events[1].size) and not impacted orbit = {"success": False, "perigee_km": None, "apogee_km": None, "eccentricity": None, "cutoff_time_s": None} landing = {"success": False, "touchdown_time_s": None, "speed_m_s": None, "miss_distance_m": None, "propellant_remaining_kg": None} released = False payload_orbit = None if ascent_complete: record("meco", time, "stack", "Main engine cutoff preserves the requested recovery reserve.", recovery_propellant_kg=float(recovery_propellant_kg)) record("fairing_opening", time, "stack", "Captive fairing opens over four seconds (kinematic mechanism).") time, state, impacted, _ = phase("stack", "fairing_opening", state, time, time + 4, coast) if ascent_complete and not impacted: separation_time = time booster, upper = separate(state, reserve_mass) pointing["booster"] = pointing["upper_stage"] = pointing["stack"] record("fairing_open", time, "stack", "Captive fairing is fully open.") record("stage_separation", time, "stack", "Upper stage released with 0.5 m/s relative radial speed.", mass_before_kg=float(state[4]), booster_mass_kg=float(booster[4]), upper_mass_kg=float(upper[4])) record("fairing_closing", time, "booster", "Fairing remains attached and closes over four seconds.") # UPPER BRANCH: coast clear, ignite, stop at orbit or fuel floor, deploy. upper_time, upper, upper_hit, _ = phase("upper_stage", "separation_coast", upper, time, time + 3, coast) fuel_floor = vehicle.upper_dry_kg + payload_kg empty = terminal_event(lambda t, y: y[4] - fuel_floor) # This monotonic speed crossing cannot skip the narrow near-circular set # between solver steps; all orbit constraints are checked at its root. target = terminal_event(lambda t, y: y[3] - math.sqrt(MU / y[0]), 1) reached = False if not upper_hit: record("upper_ignition", upper_time, "upper_stage", "Single vacuum Archimedes engine ignites.") transfer = target_altitude_km > TRANSFER_ALTITUDE_KM def transfer_apogee(t, y): apogee = orbital_elements(y)["apogee_km"] return (apogee if apogee is not None else -1e6) - target_altitude_km transfer_gate = terminal_event(transfer_apogee, 1) upper_time, upper, upper_hit, solution = phase( "upper_stage", "transfer_insertion" if transfer else "orbit_insertion", upper, upper_time, upper_time + 1200, lambda t, y: upper_guidance(t, y, vehicle, min(target_altitude_km, TRANSFER_ALTITUDE_KM) * 1000), (empty, transfer_gate if transfer else target), 12.0) cutoff_triggered = bool(solution.t_events[2].size) if transfer and cutoff_triggered and not upper_hit: record("transfer_cutoff", upper_time, "upper_stage", "Engine stops when the actual transfer-orbit apogee reaches the requested altitude.") apex = terminal_event(lambda t, y: y[2]) upper_time, upper, upper_hit, solution = phase("upper_stage", "apogee_coast", upper, upper_time, upper_time + 6000, coast, (apex,), 12.0, max(sample_step_s, 10.0)) cutoff_triggered = False if solution.t_events[1].size and not upper_hit: record("apogee", upper_time, "upper_stage", "Radial velocity crosses zero at the integrated apogee.") record("circularization_ignition", upper_time, "upper_stage", "Assumed vacuum-engine restart begins a finite circularization burn.") upper_time, upper, upper_hit, solution = phase("upper_stage", "circularization_burn", upper, upper_time, upper_time + 600, lambda t, y: upper_guidance(t, y, vehicle, target_altitude_km * 1000), (empty, target), 12.0) cutoff_triggered = bool(solution.t_events[2].size) elements = orbital_elements(upper) reached = bool(cutoff_triggered and not upper_hit and orbit_target_error(upper, target_altitude_km) <= 1e-8) orbit = dict(elements, success=reached, cutoff_time_s=upper_time, target_altitude_km=float(target_altitude_km), radial_velocity_m_s=float(upper[2]), tangential_velocity_m_s=float(upper[3]), circular_speed_m_s=math.sqrt(MU / upper[0])) if not upper_hit: stop_name = "orbit_cutoff" if reached else "upper_fuel_depleted" if upper[4] <= fuel_floor + 1e-4 else "orbit_target_missed" if cutoff_triggered else "insertion_timeout" record(stop_name, upper_time, "upper_stage", "Near-circular target orbit verified at engine cutoff." if reached else "Insertion did not reach the target orbit.") if upper_hit: record("upper_impact", upper_time, "upper_stage", "Upper stage contacted the ground.") if not upper_hit: coast_name, duration = ("deployment_coast", 10) if reached else ("failed_insertion_coast", 1200) upper_time, upper, upper_hit, _ = phase("upper_stage", coast_name, upper, upper_time, upper_time + duration, coast, body_area=12.0) if upper_hit: record("upper_impact", upper_time, "upper_stage", "Upper stage contacted the ground during unpowered flight.") if reached and not upper_hit: retained, payload = separate(upper, upper[4] - payload_kg, 0.2) pointing["payload"] = pointing["upper_stage"] record("payload_release", upper_time, "payload", "Payload released at 0.2 m/s relative radial speed; both bodies propagate independently.", payload_mass_kg=float(payload[4]), retained_mass_kg=float(retained[4]), mass_before_kg=float(upper[4])) released = True energy = orbital_elements(payload)["specific_energy_j_kg"] semimajor_axis = -MU / (2 * energy) period = 2 * math.pi * math.sqrt(semimajor_axis**3 / MU) # Keep long-coast exports compact without changing solver steps. interval = max(sample_step_s, 10.0) phase("upper_stage", "post_release_coast", retained, upper_time, upper_time + period, coast, body_area=12.0, sample_interval=interval) payload_time, _, _, solution = phase("payload", "free_flight", payload, upper_time, upper_time + period, coast, body_area=2.0, sample_interval=interval) payload_orbit = coast_verification(solution, period, target_altitude_km) record("payload_orbit_verified" if payload_orbit["success"] else "payload_orbit_failed", payload_time, "payload", "One full unpowered revolution checked against altitude and conservation limits.") # BOOSTER BRANCH: close fairing, reverse downrange motion, coast, land. bt, booster, booster_hit, _ = phase("booster", "clearance_coast", booster, time, time + 4, coast) engine = vehicle.booster_thrust_n / vehicle.booster_engines empty = terminal_event(lambda t, y: y[4] - vehicle.booster_dry_kg) initial_miss = predicted_miss(booster, recovery_target_downrange_m, contact_height) correction = -1.0 if initial_miss >= 0 else 1.0 return_target = terminal_event( lambda t, y: predicted_miss(y, recovery_target_downrange_m, contact_height), correction) if not booster_hit: record("fairing_closed", bt, "booster", "Captive fairing closed for recovery.") if not booster_hit and booster[4] > vehicle.booster_dry_kg + 1e-3: record("boostback_ignition", bt, "booster", "Three-engine range correction targets the fixed recovery ship." if ship_recovery else "Three-engine boostback targets the launch site (assumed RTLS scenario).") bt, booster, booster_hit, solution = phase("booster", "boostback", booster, bt, bt + 300, lambda t, y: np.array([0.0, correction * 3 * engine]), (empty, return_target)) if not booster_hit: reason = "predicted zero miss" if solution.t_events[2].size else "propellant depletion" if solution.t_events[1].size else "time limit" record("boostback_cutoff", bt, "booster", "Boostback stops at " + reason + ".") has_fuel = booster[4] > vehicle.booster_dry_kg + 1e-3 def landing_gate(t, y): height, vr, mass = y[0] - EARTH_RADIUS - contact_height, y[2], y[4] maximum_deceleration = max(1.0, 3 * engine / mass - MU / y[0]**2) stopping_distance = min(vr, 0)**2 / (2 * maximum_deceleration) return height - 1.5 * stopping_distance - 3000 if not booster_hit and has_fuel: gate = terminal_event(landing_gate) reentry = lambda t, y: y[0] - EARTH_RADIUS - 80_000 reentry.direction = -1 bt, booster, booster_hit, solution = phase("booster", "coast_and_reentry", booster, bt, bt + 1000, coast, (gate, reentry)) if solution.t_events[2].size: record("atmospheric_reentry", solution.t_events[2][0], "booster", "Descending through 80 km; exponential-atmosphere drag acts continuously.") if not booster_hit and solution.t_events[1].size: record("landing_burn", bt, "booster", "Feedback landing burn begins from the actual descending state.") bt, booster, booster_hit, _ = phase("booster", "landing_burn", booster, bt, bt + 600, lambda t, y: landing_guidance(t, y, vehicle, recovery_target_downrange_m, contact_height), (empty,)) if not booster_hit and booster[4] <= vehicle.booster_dry_kg + 1e-3: record("booster_fuel_depleted", bt, "booster", "Recovery propellant exhausted; remaining flight is ballistic.") bt, booster, booster_hit, _ = phase("booster", "unpowered_descent", booster, bt, bt + 1000, coast) if booster_hit: speed = math.hypot(booster[2], booster[3]) miss = abs(EARTH_RADIUS * booster[1] - recovery_target_downrange_m) success = (speed <= 2.0 and booster_contact_surface == "deck" and leg_fraction(bt, "booster") == 1 if ship_recovery else speed <= 2.0 and miss <= 100.0) contact = trajectories["booster"][-1] if contact["thrust_n"] > 0: record("booster_engine_cutoff", bt, "booster", "Engine shuts down at ground contact; contact mechanics are outside this model.") trajectories["booster"].append(dict(contact, phase="touchdown" if success else "impact", thrust_n=0.0, force_x_n=0.0, force_y_n=0.0)) landing = dict(success=bool(success), touchdown_time_s=bt, speed_m_s=speed, miss_distance_m=miss, propellant_remaining_kg=max(0.0, booster[4] - vehicle.booster_dry_kg), contact_surface=booster_contact_surface) record("touchdown" if success else "booster_impact", bt, "booster", ("Soft landing on the fixed recovery deck." if ship_recovery else "Soft landing within 100 m of launch.") if success else "Surface contact fails the speed or recovery-target criteria.") else: record("recovery_timeout", bt, "booster", "No ground contact within the recovery time limit.") else: record("launch_failure" if impacted else "ascent_timeout", time, "stack", "Vehicle contacted the ground before separation." if impacted else "Recovery-reserve cutoff was not reached.") apply_kinematic_hardware(trajectories, events, dense_phases) # Sampling is presentation only; accepted integrator steps locate all events. events.sort(key=lambda event: event["time_s"]) samples = [sample for path in trajectories.values() for sample in path] parameters = {"vehicle": asdict(vehicle), "payload_kg": float(payload_kg), "recovery_propellant_kg": float(recovery_propellant_kg), "recovery_target_downrange_m": float(recovery_target_downrange_m), "recovery_target": {"kind": "ship" if ship_recovery else "launch", "downrange_m": float(recovery_target_downrange_m), "deck_height_m": deck_height, "leg_clearance_m": leg_clearance, "length_m": RECOVERY_SHIP_LENGTH_M if ship_recovery else 0.0, "width_m": RECOVERY_SHIP_WIDTH_M if ship_recovery else 0.0, "x_m": EARTH_RADIUS * math.cos(target_angle), "y_m": EARTH_RADIUS * math.sin(target_angle)}, "target_altitude_km": float(target_altitude_km), "upper_propellant_scale": float(upper_propellant_scale), "source_url": SOURCE, "source_checked": "2026-09-26", "published_fields": ["height_m", "diameter_m", "published_liftoff_mass_kg", "advertised_leo_payload_kg", "booster_engines", "booster_thrust_n", "upper_thrust_n"], "assumed_fields": ["booster_dry_kg", "booster_propellant_kg", "booster_isp_s", "upper_dry_kg", "upper_propellant_kg", "upper_isp_s", "drag_coefficient"], "initial_mass_kg": initial_mass, "solver": "SciPy RK45 / Dormand–Prince 5(4)", "rtol": rtol, "max_step_s": max_step_s, "landing_limits": {"speed_m_s": 2.0, "miss_distance_m": RECOVERY_SHIP_LENGTH_M / 2 if ship_recovery else 100.0}, "orbit_limits": {"apsis_error_km": 2.0, "radial_speed_m_s": 5.0, "eccentricity": 0.001}, "assumptions": ["Nonrotating spherical Earth and still exponential atmosphere.", "Planar point mass with ideal instantaneous thrust direction; no attitude dynamics.", "Powered pointing aligns with thrust. A prescribed cubic coast manoeuvre pre-aligns each body before ignition (60 deg/s limit where coast time permits); no torque dynamics are claimed.", "Assumed mass split, Isp, drag, guidance, reserve, separation impulses and actuator timing.", "Above 220 km: first burn raises transfer apogee, coast reaches apogee, and an assumed vacuum-engine restart circularizes with a finite burn.", "Landing thrust envelope: 30% of one to 100% of three engines, with ideal switching.", "Published 480000 kg is a reference gross mass; actual initial mass varies with payload/fuel.", "Hungry Hippo fairing opens in four seconds and closes only after 23.2 m of integrated axial upper-stage clearance; its mass stays on the booster."], "recovery_scenario": ("Fixed downrange sea-recovery target; schematic deck and axial foot contact, no waves or rigid-body contact dynamics." if ship_recovery else "Educational return to launch site; published baseline uses downrange sea landing.")} return {"model": "Neutron planar mission v1", "parameters": parameters, "trajectories": trajectories, "events": events, "summary": {"orbit": orbit, "landing": landing, "payload_released": released, "payload_orbit": payload_orbit, "mission_success": bool(orbit["success"] and landing["success"] and released and payload_orbit and payload_orbit["success"]), "duration_s": max(s["time_s"] for s in samples), "max_dynamic_pressure_pa": max(s["dynamic_pressure_pa"] for s in samples), "function_evaluations": integrations}} if __name__ == "__main__" and "__file__" in globals(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--json", action="store_true", help="Print complete sampled mission JSON") parser.add_argument("--payload", type=float, default=8000.0, help="Payload mass in kg") args = parser.parse_args() result = simulate(payload_kg=args.payload) print(json.dumps(result if args.json else result["summary"], indent=2, allow_nan=False))