{ "cells": [ { "attachments": {}, "cell_type": "markdown", "id": "mission-introduction", "metadata": {}, "source": [ "# The complete mission, one forward pass at a time\n", "\n", "Follow the launch to separation, then trace the two independent trajectories: upper-stage delivery and booster recovery. Every worked flight result comes from the simulator's canonical Python model.\n", "\n", "This is a five-state planar translation model. Guidance reads the simulated state directly at each solver evaluation, and thrust force is applied immediately. Sensors, state estimation, attitude dynamics and actuator response are outside this model. The Python model computes the complete flight before playback begins.\n", "\n", "**Reading order:** initialize the state; calculate thrust and drag; integrate the five rates until the next event; then follow both branches, playback and results. The two branches are concurrent in physical mission time.\n", "\n", "**Model scope:** educational planar point masses, ideal pointing, spherical nonrotating Earth and a still atmosphere. Internal masses, Isp, drag, guidance and recovery scenario are assumptions. All internal calculations use SI units.\n", "\n", "To reproduce the complete flight, install NumPy 2.4.2 and SciPy 1.17.0 with Python 3.11 or newer, then **Run All**. The source is included once at the end; no external model file is required." ] }, { "attachments": {}, "cell_type": "markdown", "id": "stage1-part", "metadata": {}, "source": [ "## Stage 1: launch to separation\n", "\n", "At each trial time and state, calculate thrust, drag and all five rates. The solver repeats this calculation to advance the flight." ] }, { "attachments": {}, "cell_type": "markdown", "id": "stage1-initial", "metadata": {}, "source": [ "### 01 · Put the rocket on the pad\n", "\n", "Carry just five changing quantities. At time $t$ in seconds, $r$ is distance from Earth's centre (m), $\\theta$ is angle travelled around Earth (rad), $v_r$ is outward speed (m/s), $v_t$ is tangential speed (m/s), and $m$ is the total attached mass (kg). The subscript 0 means the launch state. The model starts at liftoff, at rest on a nonrotating spherical Earth, with full thrust applied immediately.\n", "\n", "**Symbols in this section**\n", "\n", "| Symbol | Meaning | Units |\n", "| --- | --- | --- |\n", "| $t,\\ t_0$ | Elapsed time and launch time; the reference run starts at zero. | s |\n", "| $r,\\ R_E$ | Distance from Earth's centre and the fixed spherical Earth radius. | m |\n", "| $\\theta$ | Signed angle from the launch radius, increasing in the chosen downrange direction. | rad |\n", "| $\\mathbf e_r,\\ \\mathbf e_t$ | Local unit directions: outward from Earth's centre and perpendicular to it toward increasing angle. | dimensionless |\n", "| $v_r,\\ v_t$ | Signed velocity components along $\\mathbf e_r$ and $\\mathbf e_t$; positive means outward and downrange. | m/s |\n", "| $m,\\ m_0,\\ m_u$ | Current attached mass, launch mass, and carried upper assembly mass (upper dry structure, propellant and payload). | kg |\n", "| $h,\\ s$ | Altitude above the model surface and signed downrange arc length on that surface; these are derived outputs. | m |\n", "| $\\mathbf y,\\ (\\cdot)^T,\\ (\\cdot)_0$ | Ordered five-state column vector, transpose, and launch-value subscript; its entries have different units. | mixed, by entry |\n", "\n", "$$\n", "\\mathbf y(t)=\\begin{bmatrix}r&\\theta&v_r&v_t&m\\end{bmatrix}^{\\!T},\\qquad \\mathbf y_0=\\begin{bmatrix}R_E&0&0&0&m_0\\end{bmatrix}^{\\!T}\n", "$$\n", "\n", "$$\n", "m_u=5\\,000+102\\,000+8\\,000=115\\,000\\ \\mathrm{kg}\n", "$$\n", "\n", "$$\n", "m_0=30\\,000+335\\,000+m_u=480\\,000\\ \\mathrm{kg}\n", "$$\n", "\n", "Here $R_E=6\\,378\\,137$ m is Earth radius. The carried upper assembly has mass $m_u$: upper dry mass + upper propellant + payload. Stage one contributes 30,000 kg dry mass, including its captive fairing, and 335,000 kg propellant. These mass allocations are assumptions; the 480,000 kg total matches the published lift-off mass.\n", "\n", "**Why these equations take this form**\n", "\n", "$$\n", "\\begin{aligned}h&=r-R_E,&s&=R_E\\theta\\\\v_r&=\\dot r,&v_t&=r\\dot\\theta\\end{aligned}\n", "$$\n", "\n", "Radius is measured from Earth's centre, not from the ground. Tangential velocity is the local arc rate at radius $r$; ground downrange instead uses the fixed surface radius $R_E$. Consequently $\\dot s=(R_E/r)v_t$, not generally $v_t$.\n", "\n", "$$\n", "\\begin{aligned}\\mathbf e_r&=\\begin{bmatrix}\\cos\\theta\\\\\\sin\\theta\\end{bmatrix},&\\mathbf e_t&=\\begin{bmatrix}-\\sin\\theta\\\\\\cos\\theta\\end{bmatrix}\\\\\\mathbf v&=v_r\\mathbf e_r+v_t\\mathbf e_t\\end{aligned}\n", "$$\n", "\n", "The two entries in each unit vector are components along fixed Cartesian axes in the flight plane. The first axis points from Earth's centre through the pad. The local directions rotate as $\\theta$ changes, even though this Earth does not rotate. A velocity component may be negative; speed is a nonnegative magnitude.\n", "\n", "**Checks and model boundaries**\n", "\n", "The five states describe planar translation and attached mass. There is no attitude, angular rate, out-of-plane motion, pad hold-down, Earth rotation or launch-latitude model. At liftoff the engine force is already applied; engine startup is outside this calculation.\n", "\n", "All worked numbers in this lesson use an 8,000 kg payload and 60,000 kg recovery reserve. The quoted 480,000 kg is the resulting reference-case mass; changing the payload in the executable runner changes the actual initial mass.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Earth constants** — `python/neutron/model.py:16`\n", "\n", "```python\n", "MU = 3.986004418e14 # Earth gravitational parameter [m^3/s^2]\n", "EARTH_RADIUS = 6_378_137.0 # spherical Earth radius [m]\n", "G0 = 9.80665 # Isp reference gravity [m/s^2]\n", "```\n", "\n", "**Vehicle parameters and assumed mass split** — `python/neutron/model.py:27`\n", "\n", "```python\n", "@dataclass(frozen=True)\n", "class Vehicle:\n", " \"\"\"Published geometry/thrust; explicitly assumed internal mass and Isp.\"\"\"\n", " height_m: float = 43.0\n", " diameter_m: float = 7.0\n", " published_liftoff_mass_kg: float = 480_000.0\n", " advertised_leo_payload_kg: float = 13_000.0\n", " booster_engines: int = 9\n", " booster_thrust_n: float = 1_450_000.0 * 4.4482216152605\n", " upper_thrust_n: float = 890_000.0\n", " booster_dry_kg: float = 30_000.0 # assumption, includes captive fairing\n", " booster_propellant_kg: float = 335_000.0\n", " booster_isp_s: float = 330.0\n", " upper_dry_kg: float = 5_000.0\n", " upper_propellant_kg: float = 102_000.0\n", " upper_isp_s: float = 365.0\n", " drag_coefficient: float = 0.35\n", "```\n", "\n", "**Initialize the five-state stack** — `python/neutron/model.py:361`\n", "\n", "```python\n", " area = math.pi * vehicle.diameter_m**2 / 4\n", " upper_wet = vehicle.upper_dry_kg + vehicle.upper_propellant_kg * upper_propellant_scale + payload_kg\n", " initial_mass = vehicle.booster_dry_kg + vehicle.booster_propellant_kg + upper_wet\n", " initial = np.array([EARTH_RADIUS, 0.0, 0.0, 0.0, initial_mass])\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "stage1-pitch", "metadata": {}, "source": [ "### 02 · Point the engine force\n", "\n", "At the current time, linearly interpolate the commanded pitch in degrees, $\\beta_{deg}$, between the knots below. Pitch is measured above the local horizontal: 90° points outward. Convert degrees to radians before evaluating the sine and cosine. The total nine-engine thrust $T$ stays constant through the powered phase.\n", "\n", "**Symbols in this section**\n", "\n", "| Symbol | Meaning | Units |\n", "| --- | --- | --- |\n", "| $\\beta_{deg},\\ \\beta$ | Commanded thrust-direction angle above the local positive horizontal, in degrees and in radians respectively. | °, rad |\n", "| $t_a,\\ t_b;\\ \\beta_a,\\ \\beta_b$ | Times and degree-valued pitches at the two knots bracketing the current time; a and b label endpoints. | s; ° |\n", "| $T$ | Magnitude of the total nine-engine thrust, held constant during the powered ascent. | N |\n", "| $F_r,\\ F_t$ | Signed thrust-force components in the local outward and positive downrange directions; gravity and drag are separate. | N |\n", "| $\\gamma$ | Flight-path angle of velocity above the local horizontal; defined only when speed is nonzero. | rad |\n", "\n", "$$\n", "\\beta_{deg}(t)=\\beta_a+(\\beta_b-\\beta_a)\\frac{t-t_a}{t_b-t_a}\\quad(t_a\\leq t\\leq t_b),\\qquad\\beta=\\frac{\\pi}{180}\\beta_{deg}\n", "$$\n", "\n", "$$\n", "T=1\\,450\\,000\\times4.4482216152605=6\\,449\\,921.342\\ \\mathrm N\n", "$$\n", "\n", "$$\n", "F_r=T\\sin\\beta,\\qquad F_t=T\\cos\\beta\n", "$$\n", "\n", "| Time $t$ (s) | 0 | 12 | 35 | 70 | 115 | 160 |\n", "| --- | --- | --- | --- | --- | --- | --- |\n", "| Pitch $\\beta_{deg}$ (°) | 90 | 90 | 82 | 67 | 52 | 40 |\n", "\n", "$(t_a,\\beta_a)$ and $(t_b,\\beta_b)$ are the adjacent time–pitch knots. The two force components $F_r,F_t$ act along the outward and tangential directions. Outside the knot range the endpoint pitch is held. This is an assumed pitch program with ideal pointing, not a simulated attitude controller.\n", "\n", "**Why these equations take this form**\n", "\n", "$$\n", "\\begin{aligned}\\beta=90^\\circ:&\\quad(F_r,F_t)=(T,0)\\\\\\beta=0^\\circ:&\\quad(F_r,F_t)=(0,T)\\end{aligned}\n", "$$\n", "\n", "These endpoint checks explain the sine on the radial component and cosine on the tangential component. The angle is measured from the horizontal, so the radial force is opposite the angle in the component triangle. Degree labels here describe the directions; trigonometric evaluation in Python uses radians.\n", "\n", "$$\n", "\\gamma=\\operatorname{atan2}(v_r,v_t)\\quad(V>0)\n", "$$\n", "\n", "Pointing the thrust at $\\beta$ does not instantly turn the velocity to the same angle. Forces change velocity over time, so the flight-path angle $\\gamma$ generally differs from $\\beta$. This model commands thrust direction directly; it does not integrate body attitude, gimbal motion or angle-of-attack aerodynamics.\n", "\n", "**Checks and model boundaries**\n", "\n", "The thrust conversion uses $1\\ \\mathrm{lbf}=4.4482216152605\\ \\mathrm N$; lbf means pound-force, not mass. The total is not multiplied by nine again. Constant thrust and specific impulse omit throttle, startup and ambient-pressure variation.\n", "\n", "The pitch knots are an assumed open-loop schedule. Linear interpolation makes pitch continuous at each knot, while its slope can change there. Pitch is held at 40° after the last knot until the mass event stops the powered phase.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Interpolate pitch and resolve thrust** — `python/neutron/model.py:146`\n", "\n", "```python\n", "def ascent_guidance(time, state, vehicle):\n", " \"\"\"An assumed pitch program; angles measured above local horizontal.\"\"\"\n", " pitch = np.interp(time, [0, 12, 35, 70, 115, 160], [90, 90, 82, 67, 52, 40])\n", " beta = math.radians(pitch)\n", " return vehicle.booster_thrust_n * np.array([math.sin(beta), math.cos(beta)])\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "stage1-drag", "metadata": {}, "source": [ "### 03 · Read the air from the current state\n", "\n", "Use the current radius to obtain altitude $h$ (m), then air density $\\rho$ (kg/m³). The atmosphere is still, so its relative speed $V$ is the magnitude of the two velocity components. Reference area $A$ comes from the 7 m diameter $d$; the assumed drag coefficient is $C_D=0.35$.\n", "\n", "**Symbols in this section**\n", "\n", "| Symbol | Meaning | Units |\n", "| --- | --- | --- |\n", "| $h,\\ \\rho,\\ \\rho_0$ | Altitude, atmospheric density at that altitude, and reference sea-level density (1.225 kg/m³). | m; kg/m³ |\n", "| $H$ | Density scale height: 8,500 m, the altitude increase that reduces density by a factor of e above the surface. | m |\n", "| $V=\\lVert\\mathbf v\\rVert$ | Speed relative to the still atmosphere; the magnitude of the two signed velocity components. | m/s |\n", "| $d,\\ A,\\ C_D$ | Vehicle diameter (7 m), circular reference area, and constant assumed drag coefficient (0.35). | m; m²; dimensionless |\n", "| $q_{dyn},\\ D$ | Dynamic pressure and nonnegative drag-force magnitude; pressure becomes force only after multiplying by coefficient and area. | Pa; N |\n", "| $D_r,\\ D_t$ | Signed components of the drag force, opposite the respective components of velocity. | N |\n", "\n", "$$\n", "h=r-R_E,\\qquad \\rho=1.225\\exp\\!\\left[-\\frac{\\max(0,h)}{8500\\ \\mathrm m}\\right]\\ \\mathrm{kg\\,m^{-3}}\n", "$$\n", "\n", "$$\n", "V=\\sqrt{v_r^2+v_t^2},\\qquad A=\\frac{\\pi d^2}{4}=38.485\\ \\mathrm{m^2}\n", "$$\n", "\n", "$$\n", "\\begin{bmatrix}D_r\\\\D_t\\end{bmatrix}=-\\frac12\\rho C_D A V\\begin{bmatrix}v_r\\\\v_t\\end{bmatrix}\n", "$$\n", "\n", "The drag forces $D_r,D_t$ (N) oppose velocity. At launch the speed is zero, so both drag components are zero. Density is a one-layer approximation; wind and weather are absent.\n", "\n", "**Why these equations take this form**\n", "\n", "$$\n", "\\begin{aligned}\\rho&=\\rho_0 e^{-\\max(0,h)/H}\\\\q_{dyn}&=\\frac12\\rho V^2,\\qquad D=q_{dyn}C_DA\\end{aligned}\n", "$$\n", "\n", "The density exponent is dimensionless because altitude and scale height use the same units. The max operation caps density at its sea-level value for a trial state below the surface; the separate impact event terminates physical flight at the surface.\n", "\n", "$$\n", "\\begin{aligned}\\mathbf D&=-D\\frac{\\mathbf v}{V}\\\\&=-\\frac12\\rho C_D A V\\mathbf v\\quad(V>0)\\\\\\mathbf D&=\\mathbf0\\quad(V=0)\\end{aligned}\n", "$$\n", "\n", "Divide velocity by speed to obtain its direction, then reverse it and multiply by drag magnitude. Cancelling one speed factor gives the implementation's expression, which also evaluates safely to zero at launch without division by zero. During descent the radial velocity is negative, so radial drag points outward.\n", "\n", "**Checks and model boundaries**\n", "\n", "Dimensional check: $\\rho V^2$ has units kg/(m·s²) = Pa; multiplying by $A$ gives kg·m/s² = N. $q_{dyn}$ is a pressure, not an acceleration. The solver-component index $q$ in step 5 is unrelated to dynamic pressure.\n", "\n", "This is a one-layer density law with constant drag coefficient and frontal area. It omits lift, winds, atmospheric rotation, Mach-dependent aerodynamics and heating. Fairing opening later does not alter this area or coefficient.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Altitude to density** — `python/neutron/model.py:46`\n", "\n", "```python\n", "def atmosphere(altitude_m):\n", " \"\"\"One-layer density model, not a weather or aerothermal model.\"\"\"\n", " return 1.225 * math.exp(-max(0.0, altitude_m) / 8500.0)\n", "```\n", "\n", "**Velocity-opposing drag components** — `python/neutron/model.py:103`\n", "\n", "```python\n", "def drag_force(state, area_m2, cd):\n", " \"\"\"Components in the local outward/tangential frame; still atmosphere.\"\"\"\n", " radius, _, vr, vt, _ = state\n", " speed = math.hypot(vr, vt)\n", " scale = -0.5 * atmosphere(radius - EARTH_RADIUS) * cd * area_m2 * speed\n", " return np.array([scale * vr, scale * vt])\n", "```\n", "\n", "
\n" ] }, { "attachments": { "polar-frame.svg": { "image/svg+xml": [ "PHN2ZyB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHdpZHRoPSI3NjAiIGhlaWdodD0iMzgwIiB2aWV3Qm94PSIwIDAgNzYwIDM4MCIgcm9sZT0iaW1nIiBhcmlhLWxhYmVsbGVkYnk9InRpdGxlIGRlc2MiPgo8dGl0bGUgaWQ9InRpdGxlIj5BIGxvY2FsIGJhc2lzIHRoYXQgdHVybnMgd2l0aCB0aGUgcm9ja2V0PC90aXRsZT4KPGRlc2MgaWQ9ImRlc2MiPkxlZnQ6IHIgcG9pbnRzIGZyb20gRWFydGjigJlzIGNlbnRyZSB0byB0aGUgcm9ja2V0LCB0aGV0YSBpbmNyZWFzZXMgY291bnRlcmNsb2Nrd2lzZSBmcm9tIGZpeGVkIFgsIHJhZGlhbCBwb2ludHMgb3V0d2FyZCBhbmQgdGFuZ2VudGlhbCBmb2xsb3dzIGluY3JlYXNpbmcgdGhldGEuIFJpZ2h0OiBhdCB0d28gcG9zaXRpb25zLCBib3RoIGJhc2lzIHZlY3RvcnMgdHVybiB0aHJvdWdoIHRoZSBzYW1lIGFuZ2xlIGRlbHRhIHRoZXRhLjwvZGVzYz4KPGRlZnM+PG1hcmtlciBpZD0iYXJyb3ciIHZpZXdCb3g9IjAgMCAxMCAxMCIgcmVmWD0iOCIgcmVmWT0iNSIgbWFya2VyV2lkdGg9IjYiIG1hcmtlckhlaWdodD0iNiIgb3JpZW50PSJhdXRvLXN0YXJ0LXJldmVyc2UiPjxwYXRoIGQ9Ik0wIDAgMTAgNSAwIDEwWiIgZmlsbD0iY29udGV4dC1zdHJva2UiLz48L21hcmtlcj48L2RlZnM+CjxyZWN0IHdpZHRoPSI3NjAiIGhlaWdodD0iMzgwIiByeD0iNiIgZmlsbD0iIzE0MTQxNCIvPgo8ZyBmb250LWZhbWlseT0iQXJpYWwsIHNhbnMtc2VyaWYiIGZvbnQtc2l6ZT0iMTUiIGZpbGw9IiNmNGY0ZjQiPgo8dGV4dCB4PSIyNCIgeT0iMzIiIGZvbnQtc2l6ZT0iMTciIGZvbnQtd2VpZ2h0PSI2MDAiPlBvc2l0aW9uIGFuZCB2ZWxvY2l0eSBjb21wb25lbnRzPC90ZXh0Pgo8dGV4dCB4PSI0MDQiIHk9IjMyIiBmb250LXNpemU9IjE3IiBmb250LXdlaWdodD0iNjAwIj5UaGUgYmFzaXMgY2hhbmdlcyB3aXRoIM64PC90ZXh0Pgo8cGF0aCBkPSJNMzgwIDIyVjM1MCIgc3Ryb2tlPSIjNDA0MDQwIi8+CjxwYXRoIGQ9Ik03MyAyODdIMzM1IE03MyAyODdWNjgiIGZpbGw9Im5vbmUiIHN0cm9rZT0iIzc3NyIgbWFya2VyLWVuZD0idXJsKCNhcnJvdykiLz4KPHRleHQgeD0iMzQwIiB5PSIyOTEiPlg8L3RleHQ+PHRleHQgeD0iNjgiIHk9IjU5Ij5ZPC90ZXh0Pgo8Y2lyY2xlIGN4PSI3MyIgY3k9IjI4NyIgcj0iMzQiIGZpbGw9IiMyNjMzM2YiIHN0cm9rZT0iIzgyOWJhZCIvPgo8dGV4dCB4PSIyNCIgeT0iMzQ0IiBmaWxsPSIjYjBiMGIwIj5FYXJ0aOKAmXMgY2VudHJlIE88L3RleHQ+CjxwYXRoIGQ9Ik03MyAyODdMMjI1IDE2MiIgc3Ryb2tlPSIjZGRkIiBtYXJrZXItZW5kPSJ1cmwoI2Fycm93KSIvPgo8dGV4dCB4PSIxNDIiIHk9IjIwNyI+cjwvdGV4dD4KPHBhdGggZD0iTTExOSAyODdBNDYgNDYgMCAwIDAgMTA5IDI1OCIgZmlsbD0ibm9uZSIgc3Ryb2tlPSIjZGRkIi8+Cjx0ZXh0IHg9IjEyOCIgeT0iMjc0Ij7OuDwvdGV4dD4KPHBhdGggZD0iTTIyNSAxNjJMMzA0IDk3IiBzdHJva2U9IiNmZjkyOWEiIHN0cm9rZS13aWR0aD0iMiIgbWFya2VyLWVuZD0idXJsKCNhcnJvdykiLz4KPHBhdGggZD0iTTIyNSAxNjJMMTYwIDgzIiBzdHJva2U9IiNhN2Q0ZWUiIHN0cm9rZS13aWR0aD0iMiIgbWFya2VyLWVuZD0idXJsKCNhcnJvdykiLz4KPGNpcmNsZSBjeD0iMjI1IiBjeT0iMTYyIiByPSI1IiBmaWxsPSIjZjRmNGY0Ii8+Cjx0ZXh0IHg9IjI4NCIgeT0iNzkiIGZpbGw9IiNmZjkyOWEiPnbhtaMgZeG1ozwvdGV4dD4KPHRleHQgeD0iMTE4IiB5PSI2NyIgZmlsbD0iI2E3ZDRlZSI+duKCnCBl4oKcPC90ZXh0Pgo8dGV4dCB4PSIyMzIiIHk9IjE4NSI+Um9ja2V0PC90ZXh0Pgo8dGV4dCB4PSIyNCIgeT0iMzY3IiBmb250LXNpemU9IjEzIiBmaWxsPSIjYjBiMGIwIj5l4bWjLCBl4oKcOiB1bml0IGRpcmVjdGlvbnMgwrcgduG1oywgduKCnDogc2lnbmVkIGNvbXBvbmVudHM8L3RleHQ+CjxjaXJjbGUgY3g9IjQ1MSIgY3k9IjI4MyIgcj0iNCIgZmlsbD0iI2Y0ZjRmNCIvPgo8cGF0aCBkPSJNNDUxIDI4M0g2MTQgTTQ1MSAyODNMNTY2IDE2OCIgZmlsbD0ibm9uZSIgc3Ryb2tlPSIjNzc3IiBzdHJva2UtZGFzaGFycmF5PSI0IDUiLz4KPHBhdGggZD0iTTYxNCAyODNBMTYzIDE2MyAwIDAgMCA1NjYgMTY4IiBmaWxsPSJub25lIiBzdHJva2U9IiNiMGIwYjAiIHN0cm9rZS1kYXNoYXJyYXk9IjMgNSIvPgo8cGF0aCBkPSJNNTA1IDI4M0E1NCA1NCAwIDAgMCA0ODkgMjQ1IiBmaWxsPSJub25lIiBzdHJva2U9IiNkZGQiIG1hcmtlci1lbmQ9InVybCgjYXJyb3cpIi8+Cjx0ZXh0IHg9IjUxMCIgeT0iMjUyIj7OlM64PC90ZXh0Pjx0ZXh0IHg9IjQyOSIgeT0iMzAzIj5PPC90ZXh0Pgo8ZyBzdHJva2Utd2lkdGg9IjIiIGZpbGw9Im5vbmUiPgo8cGF0aCBkPSJNNjE0IDI4M0g3MDcgTTU2NiAxNjhMNjMyIDEwMiIgc3Ryb2tlPSIjZmY5MjlhIiBtYXJrZXItZW5kPSJ1cmwoI2Fycm93KSIvPgo8cGF0aCBkPSJNNjE0IDI4M1YxOTAgTTU2NiAxNjhMNTAwIDEwMiIgc3Ryb2tlPSIjYTdkNGVlIiBtYXJrZXItZW5kPSJ1cmwoI2Fycm93KSIvPgo8L2c+CjxnIGZpbGw9IiNmNGY0ZjQiPjxjaXJjbGUgY3g9IjYxNCIgY3k9IjI4MyIgcj0iNSIvPjxjaXJjbGUgY3g9IjU2NiIgY3k9IjE2OCIgcj0iNSIvPjwvZz4KPHRleHQgeD0iNjc4IiB5PSIzMDciIGZpbGw9IiNmZjkyOWEiPmXhtaM8L3RleHQ+PHRleHQgeD0iNjIyIiB5PSIyMDciIGZpbGw9IiNhN2Q0ZWUiPmXigpw8L3RleHQ+Cjx0ZXh0IHg9IjYzOCIgeT0iMTAyIiBmaWxsPSIjZmY5MjlhIj5l4bWj4oCyPC90ZXh0Pjx0ZXh0IHg9IjQ3MCIgeT0iOTQiIGZpbGw9IiNhN2Q0ZWUiPmXigpzigLI8L3RleHQ+Cjx0ZXh0IHg9IjQwNSIgeT0iMzQwIiBmaWxsPSIjYjBiMGIwIj5JbmNyZWFzaW5nIM64IHR1cm5zIGXhtaMgdG93YXJkIGXigpwsPC90ZXh0Pgo8dGV4dCB4PSI0MDUiIHk9IjM2MiIgZmlsbD0iI2IwYjBiMCI+YW5kIHR1cm5zIGXigpwgdG93YXJkIOKIkmXhtaMuPC90ZXh0Pgo8L2c+PC9zdmc+Cg==" ] } }, "cell_type": "markdown", "id": "stage1-rates", "metadata": {}, "source": [ "### 04 · Turn forces into five simultaneous rates\n", "\n", "A dot means change per second. Evaluate every right-hand side from the same current state; together they form $\\mathbf f(t,\\mathbf y)=\\dot{\\mathbf y}$. Earth gravity uses $\\mu=3.986004418\\times10^{14}$ m³/s². The specific impulse $I_{sp}=330$ s and reference gravity $g_0=9.80665$ m/s² set the propellant flow.\n", "\n", "**Symbols in this section**\n", "\n", "| Symbol | Meaning | Units |\n", "| --- | --- | --- |\n", "| $\\dot{(\\ )}=d(\\ )/dt$ | A time derivative: rate of change per second in the fixed Earth-centred frame, with local velocity components expressed on the rotating basis. | quantity/s |\n", "| $\\mathbf f(t,\\mathbf y)$ | The right-hand-side function returning all five state rates from the same trial time and state. | mixed, by entry |\n", "| $a_r,\\ a_t$ | Physical acceleration vector components along the current local radial and tangential directions, obtained from net force per unit mass. | m/s² |\n", "| $\\dot v_r,\\ \\dot v_t$ | Rates of the stored velocity components; changing local directions makes these differ from physical acceleration components. | m/s² |\n", "| $\\mu,\\ g(r)$ | Earth's gravitational parameter and local inward gravity magnitude $g(r)=\\mu/r^2$. | m³/s²; m/s² |\n", "| $I_{sp},\\ g_0,\\ c_{eff}$ | Specific impulse (330 s), fixed reference gravity (9.80665 m/s²), and effective exhaust speed, their product. | s; m/s²; m/s |\n", "| $\\ell=r v_t$ | Signed specific angular momentum about Earth's centre; specific means per unit mass. | m²/s |\n", "\n", "$$\n", "\\begin{aligned}\\dot r&=v_r\\\\\\dot\\theta&=\\frac{v_t}{r}\\\\\\dot v_r&=\\frac{v_t^2}{r}-\\frac{\\mu}{r^2}+\\frac{F_r+D_r}{m}\\\\\\dot v_t&=-\\frac{v_rv_t}{r}+\\frac{F_t+D_t}{m}\\\\\\dot m&=-\\frac{\\sqrt{F_r^2+F_t^2}}{I_{sp}g_0}\\end{aligned}\n", "$$\n", "\n", "$$\n", "\\dot v_{r,0}=\\frac{6449921.342}{480000}-9.798285=3.639051\\ \\mathrm{m\\,s^{-2}}\n", "$$\n", "\n", "$$\n", "\\dot m_0=-1993.057383\\ \\mathrm{kg\\,s^{-1}}\n", "$$\n", "\n", "Thrust magnitude $\\sqrt{F_r^2+F_t^2}$ equals $T$ during ascent and zero during a coast. The terms $v_t^2/r$ and $-v_rv_t/r$ account for the rotating local radial/tangential basis. Thrust already represents exhaust momentum: adding a second exhaust or $\\dot m\\,v$ force would count it twice. At launch $\\dot r=\\dot\\theta=\\dot v_t=0$ mathematically.\n", "\n", "**Why these equations take this form**\n", "\n", "$$\n", "\\begin{aligned}\\frac{d\\mathbf e_r}{d\\theta}&=\\mathbf e_t,&\\frac{d\\mathbf e_t}{d\\theta}&=-\\mathbf e_r\\\\\\dot{\\mathbf e}_r&=\\dot\\theta\\mathbf e_t,&\\dot{\\mathbf e}_t&=-\\dot\\theta\\mathbf e_r\\end{aligned}\n", "$$\n", "\n", "Differentiate the sine/cosine unit vectors from step 1. Increasing $\\theta$ turns the outward direction toward the tangent and turns the tangent toward the inward direction. The chain rule contributes $\\dot\\theta=v_t/r$.\n", "\n", "$$\n", "\\begin{aligned}\\mathbf a=\\frac{d\\mathbf v}{dt}&=\\dot v_r\\mathbf e_r+v_r\\dot{\\mathbf e}_r\\\\&\\quad+\\dot v_t\\mathbf e_t+v_t\\dot{\\mathbf e}_t\\\\&=(\\dot v_r-v_t\\dot\\theta)\\mathbf e_r\\\\&\\quad+(\\dot v_t+v_r\\dot\\theta)\\mathbf e_t\\end{aligned}\n", "$$\n", "\n", "Apply the product rule to both velocity components and their directions. The cross terms appear because the directions themselves change; differentiating the two stored numbers alone would miss them.\n", "\n", "$$\n", "\\begin{aligned}a_r&=\\dot v_r-\\frac{v_t^2}{r}\\\\&=-\\frac{\\mu}{r^2}+\\frac{F_r+D_r}{m}\\\\a_t&=\\dot v_t+\\frac{v_rv_t}{r}\\\\&=\\frac{F_t+D_t}{m}\\end{aligned}\n", "$$\n", "\n", "Newton's law applies to physical acceleration. Gravity points entirely inward, so it contributes to $a_r$ only. Solve these identities for the stored component rates to obtain the two velocity equations in the ODE.\n", "\n", "$$\n", "\\begin{gathered}\\boxed{\\dot v_t=a_t-\\frac{v_rv_t}{r}}\\\\a_t=\\frac{F_t+D_t}{m}\\end{gathered}\n", "$$\n", "\n", "Tangential acceleration is not simply $v_rv_t$. The force-driven tangential acceleration is $a_t=(F_t+D_t)/m$; the extra negative term belongs to the derivative of the tangential velocity component. It includes division by radius and is required even on a nonrotating Earth. It is not an extra aerodynamic or engine force.\n", "\n", "$$\n", "\\begin{aligned}\\frac{d\\ell}{dt}=\\frac{d(rv_t)}{dt}&=v_rv_t+r\\dot v_t\\\\&=r a_t\\end{aligned}\n", "$$\n", "\n", "An independent check: with no tangential force (for example gravity-only flight with drag omitted), $a_t=0$ and $rv_t$ stays constant. As an outward-moving vehicle's radius grows, its positive tangential component must fall. During inward motion it grows. Omitting the cross term would violate this angular-momentum identity.\n", "\n", "$$\n", "\\begin{aligned}r&=6.5\\times10^6\\ \\mathrm m\\\\v_r&=1000\\ \\mathrm{m/s},\\quad v_t=2000\\ \\mathrm{m/s}\\\\a_t&=3\\ \\mathrm{m/s^2}\\\\\\frac{v_rv_t}{r}&=0.307692\\ \\mathrm{m/s^2}\\\\\\dot v_t&=3-0.307692\\\\&=2.692308\\ \\mathrm{m/s^2}\\end{aligned}\n", "$$\n", "\n", "This is a chosen local-state example, not a sampled flight result. The units are $(\\mathrm{m/s})(\\mathrm{m/s})/\\mathrm m=\\mathrm{m/s^2}$. With $v_r=-1000$ m/s at the same radius, tangential speed and force, the correction changes sign and $\\dot v_t=3.307692$ m/s². With either velocity component zero, this correction is zero.\n", "\n", "$$\n", "\\begin{aligned}c_{eff}&=I_{sp}g_0=3236.1945\\ \\mathrm{m/s}\\\\\\dot m&=-\\frac{T}{c_{eff}}\\end{aligned}\n", "$$\n", "\n", "Specific impulse is thrust divided by propellant weight flow, so multiplying it by the fixed reference gravity gives effective exhaust speed. Dividing thrust by that speed gives a positive propellant consumption rate; the vehicle mass derivative is its negative. Use $g_0$ here, not altitude-dependent $\\mu/r^2$.\n", "\n", "**Checks and model boundaries**\n", "\n", "Evaluation order follows dependencies, not a sequence of physical changes: read one complete trial state $(r,\\theta,v_r,v_t,m)$ and its time; derive altitude and speed; evaluate density, thrust direction and drag; divide the net forces by that same mass; then form all five rates. The ascent pitch schedule depends only on time, so its evaluation can occur before or after the atmosphere calculation. None of these calculations overwrites the trial state.\n", "\n", "The radial term has the same geometric origin: $\\dot v_r=a_r+v_t^2/r$. In a circular gravity-only trajectory, $v_r=0$ and $v_t^2/r=\\mu/r^2$, so $\\dot v_r=0$ even while physical acceleration is inward. Zero radial velocity derivative does not mean zero acceleration vector.\n", "\n", "Python's rhs names radial and tangential for the thrust-plus-drag acceleration components before adding gravity and the geometric terms. The returned entries [2] and [3] are the complete radial and tangential velocity derivatives. Tiny nonzero horizontal values at a nominal 90° pitch can arise from floating-point cosine evaluation.\n", "\n", "\n", "![Local radial and tangential directions turn as the rocket moves. Static coordinate illustration; not a flight result.](attachment:polar-frame.svg)\n", "\n", "*Local radial and tangential directions turn as the rocket moves. Static coordinate illustration; not a flight result.*\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**The five simultaneous state derivatives** — `python/neutron/model.py:111`\n", "\n", "```python\n", "def rhs(time, state, guidance, area_m2, cd, isp_s):\n", " \"\"\"Polar Newton equations; state = [r, theta, v_r, v_t, mass].\n", "\n", " Pointing is ideal: guidance commands force directly, not attitude/torque.\n", " The changing mass appears in F/m and dm/dt; thrust includes exhaust momentum.\n", " \"\"\"\n", " radius, _, vr, vt, mass = state\n", " force = guidance(time, state)\n", " radial, tangential = (force + drag_force(state, area_m2, cd)) / mass\n", " return [vr, vt / radius,\n", " vt * vt / radius - MU / radius**2 + radial,\n", " -vr * vt / radius + tangential,\n", " -np.linalg.norm(force) / (isp_s * G0)]\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "stage1-integrate", "metadata": {}, "source": [ "### 05 · Advance the state, then repeat\n", "\n", "SciPy RK45 evaluates the rates several times inside a trial time step $\\Delta t$. Intermediate slopes $\\mathbf k_i$ use intermediate trial states, so pitch, density, drag and acceleration are recalculated there. $a_{ij},b_i,\\widehat b_i,c_i$ are the fixed Dormand–Prince coefficients, not vehicle parameters.\n", "\n", "**Symbols in this section**\n", "\n", "| Symbol | Meaning | Units |\n", "| --- | --- | --- |\n", "| $n,\\ i,\\ j,\\ q$ | Accepted-step index; trial-stage index; earlier-stage index; state-component index from 1 to 5. None is time itself. | dimensionless indices |\n", "| $t_n,\\ \\Delta t,\\ \\mathbf y_n$ | Time of the current accepted state, proposed time increment, and current accepted five-state vector. | s; s; state units |\n", "| $\\mathbf k_i$ | Five-rate vector evaluated at trial stage i; its components have each state's units divided by seconds. | state units/s |\n", "| $a_{ij},\\ c_i,\\ b_i,\\ \\widehat b_i$ | Fixed method coefficients: trial-state weights, time fractions, fifth-order weights and embedded fourth-order weights. | dimensionless |\n", "| $\\mathbf y_{n+1}^{(5)},\\ \\mathbf e$ | Proposed fifth-order state and difference between the embedded fifth- and fourth-order estimates; (5) means method order, not exponent. | state units |\n", "| $\\mathrm{atol}_q,\\ \\mathrm{rtol},\\ s_q$ | Component absolute tolerance, dimensionless relative tolerance (10⁻⁷), and combined error scale for that component. | state units; dimensionless; state units |\n", "| $e_q,\\ E$ | One component of the estimated local error, and the root-mean-square of all five scaled errors. | state units; dimensionless |\n", "| $\\tau,\\ \\widetilde{\\mathbf y}_i$ | Dummy time variable inside an integral and the complete trial state passed to rate evaluation i. | s; state units |\n", "| $\\mathbf y^{[i]},\\ \\mathbf y^{(5)}$ | Equivalent web shorthand: $\\mathbf y^{[i]}=\\widetilde{\\mathbf y}_i$ is a trial state; $\\mathbf y^{(5)}=\\mathbf y_{n+1}^{(5)}$ is the fifth-order candidate. | state units |\n", "\n", "$$\n", "\\mathbf k_i=\\mathbf f\\!\\left(t_n+c_i\\Delta t,\\ \\mathbf y_n+\\Delta t\\sum_{j0\\\\v_r(\\Delta t)&\\approx\\dot v_{r,0}\\Delta t\\end{aligned}\n", "$$\n", "\n", "Near launch, outward acceleration first creates velocity, and the time integral of that increasing velocity creates altitude. Euler using only the zero initial velocity would produce no altitude change on its first step. RK45 samples accelerated intermediate states, which is why its first accepted altitude is positive. This describes continuous causation, not permission to overwrite velocity before evaluating the other rates.\n", "\n", "$$\n", "\\mathbf e=\\mathbf y_{n+1}^{(5)}-\\mathbf y_{n+1}^{(4)},\\qquad\\left[\\frac{e_q}{s_q}\\right]=1\n", "$$\n", "\n", "The embedded methods reuse rate evaluations to estimate local error. Scaling each component before combining them prevents metres, radians, metres per second and kilograms from being compared directly. SciPy may use the opposite sign for the difference; the squared norm and acceptance decision are identical.\n", "\n", "$$\n", "\\begin{aligned}s_r\\big|_{r\\approx R_E}&\\approx10^{-3}\\ \\mathrm m\\\\&\\quad+10^{-7}(6\\,378\\,137\\ \\mathrm m)\\\\&\\approx0.639\\ \\mathrm m\\end{aligned}\n", "$$\n", "\n", "The relative tolerance acts on the stored Earth-centred radius, not altitude. Thus a radius absolute tolerance of one millimetre does not imply one-millimetre altitude accuracy. This scale participates in an RMS acceptance test across all five states; it is not a guaranteed per-component or global error bound.\n", "\n", "**Checks and model boundaries**\n", "\n", "The coefficient table uses $j$ for earlier Runge–Kutta stages. The web error-control panel also uses $j$ to select a state component; that role is written $q$ in the download. In the error formulas, $e_j,s_j,\\mathrm{atol}_j$ therefore mean exactly $e_q,s_q,\\mathrm{atol}_q$; neither index denotes a physical variable.\n", "\n", "The full loop is: accepted time and five-state vector → complete trial stage → derived altitude, speed and density plus commanded thrust → drag and net forces → all five rates → further trial stages → candidate state and error estimate → accept the whole vector or retry the whole step. Force evaluation uses a decreasing trial mass on later stages, so rising thrust acceleration is included without a separate mass-first update.\n", "\n", "A rejected trial does not become a new flight state. RK45 retries from the same accepted state with a smaller time increment. Several trial force evaluations can therefore occur before one state is accepted; a function-evaluation count is not a count of accepted steps.\n", "\n", "The 2 s maximum step is a ceiling, and the 2 s output interval is a separate display choice. The solver selects its own first and later steps. Dense output uses a quartic interpolating polynomial inside accepted RK45 steps; it supplies both display samples and states used by the event root search.\n", "\n", "A small numerical error estimate cannot validate the physical assumptions. Check numerical convergence by tightening tolerances and maximum step and comparing event times and separation states. The download's automated checks also compare against the canonical simulator, analytic burn time and separation momentum.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Configure SciPy's adaptive RK45 solver** — `python/neutron/model.py:133`\n", "\n", "```python\n", "def integrate(state, start_s, end_s, guidance, area_m2, cd, isp_s,\n", " events=(), max_step_s=2.0, rtol=1e-7):\n", " \"\"\"Integrate one continuous phase; restart after every discrete change.\"\"\"\n", " solution = solve_ivp(\n", " lambda t, y: rhs(t, y, guidance, area_m2, cd, isp_s),\n", " (start_s, end_s), state, method=\"RK45\", dense_output=True,\n", " events=events, max_step=max_step_s, rtol=rtol,\n", " atol=[1e-3, 1e-11, 1e-5, 1e-5, 1e-4])\n", " if not solution.success:\n", " raise RuntimeError(solution.message)\n", " return solution\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "stage1-cutoff", "metadata": {}, "source": [ "### 06 · Stop burning at the reserve mass\n", "\n", "Keep 60,000 kg of stage-one propellant for later recovery. The booster mass remaining at cutoff, $m_b$, is its dry mass plus this reserve. While the upper assembly is attached, the cutoff stack mass is $m_{cut}$. The solver locates the downward zero crossing of the mass gate $g_{MECO}$ between accepted steps.\n", "\n", "**Symbols in this section**\n", "\n", "| Symbol | Meaning | Units |\n", "| --- | --- | --- |\n", "| $m_b,\\ m_u$ | Booster dry mass plus retained recovery propellant, and the still-attached upper assembly mass. | kg |\n", "| $m_{cut}=m_b+m_u$ | Total attached mass at which the ascent burn stops; retained propellant is still part of this mass. | kg |\n", "| $g_{MECO}(t,\\mathbf y),\\ g(t,\\mathbf y)$ | Event residual: current mass minus cutoff mass; g is the web shorthand for this same gate, not gravitational acceleration. | kg |\n", "| $t_{MECO},\\ t_c$ | Elapsed time of main engine cutoff, located at the downward zero crossing of the mass residual; the two symbols name the same time. | s |\n", "\n", "$$\n", "m_b=30\\,000+60\\,000=90\\,000\\ \\mathrm{kg},\\qquad m_{cut}=m_b+m_u=205\\,000\\ \\mathrm{kg}\n", "$$\n", "\n", "$$\n", "g_{MECO}(t,\\mathbf y)=m-m_{cut}=0\\quad\\text{with decreasing }m\n", "$$\n", "\n", "$$\n", "t_{MECO}=\\frac{m_0-m_{cut}}{T/(I_{sp}g_0)}=137.978968\\ \\mathrm s\n", "$$\n", "\n", "Constant thrust and specific impulse make this burn-time check analytic. The integrated cutoff is at 56.588 km altitude and 40.314 km downrange. A separate downward crossing of $r-R_E=0$ stops integration on ground impact; the powered phase also has a 600 s safety limit. Neither stop is treated as successful MECO.\n", "\n", "**Why these equations take this form**\n", "\n", "$$\n", "\\begin{aligned}m(t)&=m_0-\\frac{T}{I_{sp}g_0}t\\\\m(t_{MECO})&=m_{cut}\\end{aligned}\n", "$$\n", "\n", "Integrate the constant mass derivative during the powered phase and solve for the crossing time. The 275,000 kg burned is initial mass minus cutoff mass; the 60,000 kg reserve is retained, so MECO is not fuel exhaustion. This analytic time works because thrust and specific impulse are constant.\n", "\n", "$$\n", "\\begin{aligned}g_{MECO}>0:&\\quad\\text{keep burning}\\\\g_{MECO}=0:&\\quad\\text{cutoff}\\\\\\dot g_{MECO}&=\\dot m<0\\end{aligned}\n", "$$\n", "\n", "The terminal event requests the decreasing crossing. A sign change between accepted states brackets a root; the solver refines its time using the continuous interpolated solution and returns that root state. The event need not coincide with a 2 s output sample, and the model does not wait for the next display frame.\n", "\n", "**Checks and model boundaries**\n", "\n", "Changing only payload adds the same mass to launch and cutoff, so the analytic burn time is unchanged at fixed thrust and reserve, even though the acceleration and trajectory change. Increasing the reserve reduces the mass burned and shortens the burn.\n", "\n", "The ground residual is altitude in metres and the MECO residual is mass in kilograms; an event function needs a zero crossing, not common units with other events. Ground impact and the 600 s phase limit are failure paths, checked separately before the coast is started.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Set terminal root and crossing direction** — `python/neutron/model.py:126`\n", "\n", "```python\n", "def terminal_event(function, direction=-1):\n", " \"\"\"SciPy locates a zero crossing between accepted RK45 steps.\"\"\"\n", " function.terminal = True\n", " function.direction = direction\n", " return function\n", "```\n", "\n", "**Reserve-mass MECO gate and powered ascent** — `python/neutron/model.py:450`\n", "\n", "```python\n", " coast = lambda t, y: np.zeros(2)\n", " reserve_mass = vehicle.booster_dry_kg + recovery_propellant_kg\n", " meco = terminal_event(lambda t, y: y[4] - upper_wet - reserve_mass)\n", " record(\"liftoff\", 0, \"stack\", \"Nine-engine ascent begins.\")\n", " time, state, impacted, ascent_solution = phase(\"stack\", \"ascent\", initial, 0, 600,\n", " lambda t, y: ascent_guidance(t, y, vehicle), (meco,))\n", " ascent_complete = bool(ascent_solution.t_events[1].size) and not impacted\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "stage1-coast", "metadata": {}, "source": [ "### 07 · Coast for four seconds while the fairing opens\n", "\n", "At MECO, set both thrust components to zero and restart the same ODE from the cutoff state. Gravity and drag still act, while mass is constant. The captive fairing remains on the booster; opening is a prescribed four-second animation and does not change the force model or reference area.\n", "\n", "**Symbols in this section**\n", "\n", "| Symbol | Meaning | Units |\n", "| --- | --- | --- |\n", "| $t_{MECO}=t_c,\\ t_{sep}=t_s$ | Main engine cutoff time and stage-separation time, four seconds later in the reference sequence; each pair is equivalent notation. | s |\n", "| $\\mathbf y(t^-),\\ \\mathbf y(t^+)$ | State immediately before and immediately after a change of phase; minus and plus indicate one-sided limits, not negative and positive times. | state units |\n", "| $f_{open},\\ \\operatorname{clip}(z,0,1)$ | Prescribed fairing-opening fraction, and a clamp that limits its argument z to the interval from zero to one. | dimensionless |\n", "\n", "$$\n", "F_r=F_t=0,\\qquad \\dot m=0,\\qquad \\mathbf y(t_{MECO}^{+})=\\mathbf y(t_{MECO}^{-})\n", "$$\n", "\n", "$$\n", "t_{sep}=t_{MECO}+4\\ \\mathrm s=141.978968\\ \\mathrm s\n", "$$\n", "\n", "$$\n", "f_{open}(t)=\\operatorname{clip}\\!\\left(\\frac{t-t_{MECO}}{4\\ \\mathrm s},0,1\\right)\n", "$$\n", "\n", "The superscripts − and + mean immediately before and after the mode change. The opening fraction $f_{open}$ runs from 0 (closed) to 1 (open). Integration gives 60.581 km altitude at separation. Ground contact remains a terminal event throughout the coast.\n", "\n", "**Why these equations take this form**\n", "\n", "$$\n", "\\begin{aligned}\\dot v_r&=\\frac{v_t^2}{r}-\\frac{\\mu}{r^2}+\\frac{D_r}{m}\\\\\\dot v_t&=-\\frac{v_rv_t}{r}+\\frac{D_t}{m}\\\\\\dot m&=0\\end{aligned}\n", "$$\n", "\n", "Substitute zero thrust into the same equations. Coast means no engine thrust; it does not mean constant velocity or zero acceleration. Gravity, drag and the conversion between physical acceleration and rotating velocity components still apply.\n", "\n", "$$\n", "\\begin{aligned}f_{open}(t_{MECO})&=0\\\\f_{open}(t_{MECO}+2\\ \\mathrm s)&=\\frac12\\\\f_{open}(t_{sep})&=1\\end{aligned}\n", "$$\n", "\n", "The four-second denominator gives a dimensionless fraction. It is prescribed mechanism progress, not an extra ODE state. The fairing stays attached, so opening does not remove mass.\n", "\n", "**Checks and model boundaries**\n", "\n", "At cutoff, position, velocity and mass stay continuous because no impulse or mass drop occurs. Their derivatives may jump when thrust turns off. The solver is restarted at this known change so a single integration phase does not straddle the thrust discontinuity.\n", "\n", "At an event the ordering is explicit: locate and retain the common event time and state, select the next phase's thrust command, then evaluate fresh rates from that unchanged state under the new command. At separation in step 8, apply the mass partition and paired velocity impulse before any subsequent branch integration. Never reuse the old phase's force or rate vector after a mode change.\n", "\n", "The animation is kinematic: it adds no actuator forces, torques, drag-area change, contact or structural dynamics. Separation is attempted only after the full four-second coast without ground contact.\n", "\n", "The function clip limits its first argument to the stated lower and upper bounds.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Restart from MECO for the four-second coast** — `python/neutron/model.py:464`\n", "\n", "```python\n", " if ascent_complete:\n", " record(\"meco\", time, \"stack\", \"Main engine cutoff preserves the requested recovery reserve.\",\n", " recovery_propellant_kg=float(recovery_propellant_kg))\n", " record(\"fairing_opening\", time, \"stack\", \"Captive fairing opens over four seconds (kinematic mechanism).\")\n", " time, state, impacted, _ = phase(\"stack\", \"fairing_opening\", state, time, time + 4, coast)\n", "```\n", "\n", "**Recorded kinematic fairing opening and closure** — `python/neutron/model.py:379`\n", "\n", "```python\n", " def fairing_fraction(time, body):\n", " opening = next((e[\"time_s\"] for e in events if e[\"name\"] == \"fairing_opening\"), None)\n", " if opening is None or body not in (\"stack\", \"booster\"):\n", " return 0.0\n", " if separation_time is None or time <= separation_time:\n", " return float(np.clip((time - opening) / 4, 0, 1))\n", " return float(np.clip(1 - (time - separation_time) / 4, 0, 1))\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "stage1-separate", "metadata": {}, "source": [ "### 08 · Split the mass and preserve momentum\n", "\n", "At the end of the coast, the common position and tangential velocity pass to both bodies. Apply equal and opposite radial impulses with prescribed relative separation speed $\\Delta v=0.5$ m/s. Subscripts $b$ and $u$ denote booster and carried upper assembly; the superscript − denotes the attached state just before separation.\n", "\n", "**Symbols in this section**\n", "\n", "| Symbol | Meaning | Units |\n", "| --- | --- | --- |\n", "| $b,\\ u,\\ (\\cdot)^-$ | Booster subscript, upper-assembly subscript, and the common attached state immediately before separation. | labels |\n", "| $m^-,\\ m_b,\\ m_u$ | Mass before separation and the two positive body masses after it; their sum is unchanged. | kg |\n", "| $\\Delta v=v_{r,u}-v_{r,b}$ | Prescribed relative outward separation speed (0.5 m/s); this is shared between the bodies, not applied in full to each. | m/s |\n", "| $J$ | Magnitude of the instantaneous internal radial impulse: positive outward on the upper assembly, equally negative on the booster. | N·s = kg·m/s |\n", "| $\\Delta K$ | Increase in the pair's translational kinetic energy due to the prescribed separation impulse. | J (joules) |\n", "\n", "$$\n", "m^-=m_b+m_u,\\qquad r_b=r_u=r^-,\\qquad\\theta_b=\\theta_u=\\theta^-,\\qquad v_{t,b}=v_{t,u}=v_t^-\n", "$$\n", "\n", "$$\n", "v_{r,b}=v_r^- -\\Delta v\\frac{m_u}{m^-},\\qquad v_{r,u}=v_r^- +\\Delta v\\frac{m_b}{m^-}\n", "$$\n", "\n", "$$\n", "m_bv_{r,b}+m_uv_{r,u}=m^-v_r^-,\\qquad v_{r,u}-v_{r,b}=\\Delta v\n", "$$\n", "\n", "$$\n", "v_{r,b}=978.963072\\ \\mathrm{m/s},\\qquad v_{r,u}=979.463072\\ \\mathrm{m/s}\n", "$$\n", "\n", "The attached-stack calculation is complete. Carry each separated body's state into its own branch below: upper-stage delivery and booster recovery. Separation geometry and contact dynamics are not modeled.\n", "\n", "**Why these equations take this form**\n", "\n", "$$\n", "\\begin{aligned}m_b(v_{r,b}-v_r^-)&=-J\\\\m_u(v_{r,u}-v_r^-)&=J\\\\\\Delta v&=J\\left(\\frac1{m_u}+\\frac1{m_b}\\right)\\end{aligned}\n", "$$\n", "\n", "Start with equal and opposite impulses. Divide each impulse by that body's mass to obtain its velocity change, then subtract the two velocities to impose the requested relative speed. Over this ideal instantaneous split, the accumulated impulse from continuous external forces is neglected.\n", "\n", "$$\n", "\\begin{aligned}J&=\\Delta v\\frac{m_bm_u}{m_b+m_u}\\\\\\Delta v_{r,b}&=-0.280488\\ \\mathrm{m/s}\\\\\\Delta v_{r,u}&=+0.219512\\ \\mathrm{m/s}\\end{aligned}\n", "$$\n", "\n", "For 90,000 kg of booster and 115,000 kg of upper assembly, the lighter booster receives the larger speed change. Their difference is 0.5 m/s and their mass-weighted sum is zero, giving the momentum-preserving update above. These are increments, not the final radial velocities.\n", "\n", "$$\n", "\\Delta K=\\frac12\\frac{m_bm_u}{m_b+m_u}(\\Delta v)^2\\approx6311\\ \\mathrm J\n", "$$\n", "\n", "Momentum is conserved by the internal impulse, but translational kinetic energy increases. The prescribed impulse implicitly supplies separation energy; this model does not solve a spring, pneumatic system or contact mechanism. Energy conservation of the two bodies' translation alone is therefore not the intended separation check.\n", "\n", "**Checks and model boundaries**\n", "\n", "Both point masses initially occupy the common pre-separation position. The radial impulse changes their radial velocities immediately; position does not jump. No tangential impulse is applied, so both tangential velocities retain the common value.\n", "\n", "There is no mass jettison at this split: the original attached mass is partitioned between two bodies. The captive fairing belongs to booster dry mass. The stage-one calculation ends with two initial separated states, each ready for its own subsequent propagation.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Split mass and apply the relative radial impulse** — `python/neutron/model.py:198`\n", "\n", "```python\n", "def separate(state, booster_mass_kg, relative_speed_m_s=0.5):\n", " \"\"\"Split mass with equal/opposite impulses: conserve linear momentum.\n", "\n", " Both bodies start at the same position. Geometry and contact are not modeled.\n", " The reusable fairing stays within booster dry mass.\n", " \"\"\"\n", " booster, upper = state.copy(), state.copy()\n", " upper[4] = state[4] - booster_mass_kg\n", " booster[4] = booster_mass_kg\n", " booster[2] -= relative_speed_m_s * upper[4] / state[4]\n", " upper[2] += relative_speed_m_s * booster[4] / state[4]\n", " return booster, upper\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "mission-branch-point", "metadata": {}, "source": [ "## One separation. Two trajectories.\n", "\n", "Separation creates two independent initial-value problems at the same mission time. Each branch now carries its own five-number state; $m$ is the mass of that body. Throughout these sections, $h=r-R_E$ means altitude and $\\ell=rv_t$ means specific angular momentum; they are different quantities. Both reuse the atmosphere, drag, ODE and adaptive RK45 loop already established above. Within each solver evaluation, every command and derivative comes from the same trial state: guidance reads it, force limits bound the command, and the ODE returns all five derivatives together. RK45 then constructs its next trial or accepted state. These equations describe dependencies, not sequential in-place updates of radius followed by velocity followed by mass. Phase-specific force commands, reference areas, propellant parameters and event gates use that same ODE. Read either branch forward from separation; the booster lands while the upper stage is still flying." ] }, { "attachments": {}, "cell_type": "markdown", "id": "branch-upper", "metadata": {}, "source": [ "## Upper stage → payload orbit\n", "\n", "Carry the upper state forward through clearance, powered insertion, apogee coast, circularization, deployment and one real unpowered orbit." ] }, { "attachments": {}, "cell_type": "markdown", "id": "upper-clearance", "metadata": {}, "source": [ "### 09 · Coast clear of the booster\n", "\n", "Start with the upper assembly state from separation. Here $t_{sep}$ is the shared stage-separation time. Its clock is still mission elapsed time; it does not restart at zero. For three seconds, feed zero engine force into the same five-state ODE.\n", "\n", "**Symbols in this section.**\n", "\n", "- $t,\\ t_{sep},\\ t_{ign}$ — Mission elapsed time, stage-separation time, and upper-engine ignition time; all share the same clock. **Units:** s.\n", "- $F_r,\\ F_t$ — Engine-force components: positive radially outward and in the direction of increasing angle, respectively. These exclude gravity and drag. **Units:** N.\n", "- $m,\\ \\dot m$ — Mass of the upper stage plus attached payload, and its time derivative. A dot always means differentiation with respect to mission time. **Units:** kg; kg/s.\n", "\n", "$$\n", "t_{ign}=t_{sep}+3\\ \\mathrm s\n", "$$\n", "\n", "The upper engine ignition time is three seconds after the shared separation time.\n", "\n", "$$\n", "F_r=F_t=0,\\qquad\\dot m=0\n", "$$\n", "\n", "Gravity and drag keep changing velocity; mass stays constant during clearance.\n", "\n", "**Computed reference flight.** Separation is at 141.979 s; upper ignition is at 144.979 s. The carried mass is 115,000 kg.\n", "\n", "Use the integrated clearance state as the burn initial condition. The current implementation keeps the 38.485 m² reference area during this coast, then changes it to 12 m² from ignition onward. This timed coast assumes clearance; body contact and plume interaction are not modeled.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Three-second clearance and upper-stage fuel floor** — `python/neutron/model.py:480`\n", "\n", "```python\n", " upper_time, upper, upper_hit, _ = phase(\"upper_stage\", \"separation_coast\", upper, time, time + 3, coast)\n", " fuel_floor = vehicle.upper_dry_kg + payload_kg\n", " empty = terminal_event(lambda t, y: y[4] - fuel_floor)\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "upper-guidance", "metadata": {}, "source": [ "### 10 · Command the upper-stage engine force\n", "\n", "The target altitude $h_{target}$ is 400 km. First command a 220 km reference altitude while building horizontal speed. Let $h_c$ be the commanded altitude and $\\tau$ = 30 s the response time. The requested radial-speed derivative is $a_r^*$. A star marks a requested value. The vacuum-engine thrust is $T_u$ = 890,000 N, with assumed specific impulse $I_{sp,u}$ = 365 s.\n", "\n", "**Symbols in this section.**\n", "\n", "- $R_E,\\ r$ — Fixed spherical Earth radius (6,378,137 m), and vehicle distance from Earth's centre. **Units:** m.\n", "- $h_c,\\ h_{target}$ — Commanded altitude in this burn, and requested final orbital altitude; altitude is measured above the spherical surface. **Units:** m.\n", "- $v_r,\\ v_t$ — Radial speed (positive outward) and tangential speed (positive with increasing angle). In particular, v_t is a linear speed, not an angular rate. **Units:** m/s.\n", "- $\\mu$ — Earth gravitational parameter, 3.986004418 × 10¹⁴; local gravitational acceleration magnitude is μ/r². **Units:** m³/s².\n", "- $m$ — Current upper-stage plus attached-payload mass, including remaining propellant. **Units:** kg.\n", "- $\\tau$ — Chosen radial feedback response time, 30 s; smaller values request faster correction. **Units:** s.\n", "- $a_r^*$ — Requested derivative of radial speed, dot v_r, before command limits; the star denotes a request, not a measured or guaranteed acceleration. **Units:** m/s².\n", "- $T_u,\\ q$ — Upper-engine thrust magnitude (890,000 N), and its signed radial fraction F_r/$T_u$. The clipped fraction lies between −0.5 and 0.95. **Units:** N; dimensionless.\n", "- $F_r,\\ F_t$ — Commanded engine forces in the outward radial and positive tangential directions. **Units:** N.\n", "- $I_{sp,u},\\ g_0$ — Assumed upper-engine specific impulse (365 s), and fixed reference gravity (9.80665 m/s²). g_0 converts specific impulse to effective exhaust speed; it is not local gravity. **Units:** s; m/s².\n", "- $\\dot m$ — Rate of change of this body's mass; negative while propellant is consumed. **Units:** kg/s.\n", "- $e_h$ — Altitude error r − (R_E + $h_c$), used only in the ideal feedback derivation below; distinct from orbital eccentricity e. **Units:** m.\n", "\n", "$$\n", "a_r^*=\\frac{R_E+h_c-r}{\\tau^2}-\\frac{2v_r}{\\tau}\n", "$$\n", "\n", "Altitude error requests acceleration; radial velocity adds damping. Initially $h_c$ = min($h_{target}$, 220 km).\n", "\n", "$$\n", "\\begin{aligned}q=\\operatorname{clip}\\!\\Big[&\\frac{m}{T_u}\\left(a_r^*+\\frac{\\mu}{r^2}-\\frac{v_t^2}{r}\\right),\\\\&-0.5,0.95\\Big]\\end{aligned}\n", "$$\n", "\n", "The radial thrust fraction $q$ compensates gravity and the polar-basis term. The function clip limits its first argument to the stated lower and upper bounds. This command does not compensate drag.\n", "\n", "$$\n", "F_r=T_uq,\\qquad F_t=T_u\\sqrt{1-q^2}\n", "$$\n", "\n", "The two force components give constant total thrust and positive tangential thrust.\n", "\n", "$$\n", "\\dot m=-\\frac{T_u}{I_{sp,u}g_0}\n", "$$\n", "\n", "Use these force components and propellant flow in the established ODE; retain drag with reference area 12 m².\n", "\n", "**Why these equations have this form.**\n", "\n", "$$\n", "\\begin{aligned}e_h&=r-(R_E+h_c)\\\\\\ddot e_h+\\frac{2}{\\tau}\\dot e_h+\\frac{1}{\\tau^2}e_h&=0\\end{aligned}\n", "$$\n", "\n", "For a constant altitude command, ignored drag and an unsaturated engine command, substitution into the ODE gives this critically damped error equation. It explains the feedback gains: altitude error is corrected without oscillation in this ideal limit. Clipping, finite thrust and drag break that ideal equation.\n", "\n", "**Computed reference flight.** At ignition the upper assembly has mass 115,000 kg. The upper burn consumes 248.643 kg/s.\n", "\n", "At each RK45 trial state, read radius, velocities and mass; compute altitude and radial-speed errors; form $a_r^*$; then compute and clip the radial force fraction before constructing both force components. Only after that does the shared ODE combine engine force with drag and gravity to calculate all five state derivatives. The integrator, not the guidance law, updates the state. Constant thrust does not mean constant acceleration: mass and gravity change. Clipping and uncompensated drag can prevent the requested radial acceleration from being achieved. Here $q$ is a scalar thrust fraction, distinct from dynamic pressure or a quaternion.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Requested radial acceleration to fixed-magnitude thrust** — `python/neutron/model.py:153`\n", "\n", "```python\n", "def upper_guidance(time, state, vehicle, target_m):\n", " \"\"\"PD radial acceleration request plus spherical-gravity compensation.\"\"\"\n", " radius, _, vr, vt, mass = state\n", " tau = 30.0\n", " requested = (EARTH_RADIUS + target_m - radius) / tau**2 - 2 * vr / tau\n", " radial_fraction = (requested + MU / radius**2 - vt * vt / radius) * mass / vehicle.upper_thrust_n\n", " beta = math.asin(np.clip(radial_fraction, -0.5, 0.95))\n", " return vehicle.upper_thrust_n * np.array([math.sin(beta), math.cos(beta)])\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "upper-transfer", "metadata": {}, "source": [ "### 11 · Burn until the trajectory reaches the target apogee\n", "\n", "From the current position and velocity, infer the instantaneous unpowered two-body orbit. Its specific energy $\\varepsilon$ and specific angular momentum $\\ell$ determine whether that instantaneous orbit is bound. The dimensionless eccentricity $e$ and semi-latus rectum $p$ (m) locate its low and high points.\n", "\n", "**Symbols in this section.**\n", "\n", "- $r,\\ R_E$ — Current distance from Earth's centre and fixed spherical surface radius. **Units:** m.\n", "- $v_r,\\ v_t$ — Signed outward radial speed and signed tangential speed. **Units:** m/s.\n", "- $\\mu$ — Earth gravitational parameter, 3.986004418 × 10¹⁴. **Units:** m³/s².\n", "- $\\varepsilon$ — Specific orbital mechanical energy: kinetic plus gravitational potential energy per unit vehicle mass, with zero potential at infinity. **Units:** J/kg = m²/s².\n", "- $\\ell$ — Signed specific angular momentum normal to the simulation plane; positive for increasing angle. It is angular momentum per unit mass. **Units:** m²/s.\n", "- $e,\\ p$ — Nonnegative orbital eccentricity, and semi-latus rectum (the orbit's radius 90° from perigee in the two-body ellipse). These are not vehicle mass or thrust parameters. **Units:** dimensionless; m.\n", "- $h_a,\\ h_{target}$ — Apogee altitude of the instantaneous two-body orbit, and requested final orbit altitude. Neither is necessarily the current altitude. **Units:** m.\n", "- $\\uparrow$ — An event function crossing zero from negative to positive. The arrow describes the gate value, not necessarily the vehicle's vertical motion. **Units:** event direction.\n", "\n", "$$\n", "\\varepsilon=\\frac{v_r^2+v_t^2}{2}-\\frac{\\mu}{r},\\qquad\\ell=rv_t\n", "$$\n", "\n", "Energy is in J/kg; angular momentum is in m²/s. These describe the current state, not a prescribed flight path.\n", "\n", "$$\n", "\\begin{aligned}e^2={}&\\left(\\frac{rv_t^2}{\\mu}-1\\right)^2\\\\&+\\left(\\frac{rv_rv_t}{\\mu}\\right)^2\\end{aligned}\n", "$$\n", "\n", "Take the nonnegative square root for eccentricity $e$. This formula avoids cancellation near a circular orbit.\n", "\n", "$$\n", "p=\\frac{\\ell^2}{\\mu},\\qquad h_a=\\frac{p}{1-e}-R_E\n", "$$\n", "\n", "For $\\varepsilon < 0$ and $e < 1$, $h_a$ is the osculating apogee altitude. An unbound state has no finite apogee here.\n", "\n", "$$\n", "h_a-h_{target}=0\\quad(\\uparrow)\n", "$$\n", "\n", "Locate the upward crossing of the target apogee, here $h_{target}$ = 400 km; then cut thrust.\n", "\n", "**Computed reference flight.** Transfer cutoff occurs at 553.236 s, while the vehicle itself is at 219.980 km altitude. Its osculating apogee has reached 400 km.\n", "\n", "The event finder locates the cutoff on the integrated solution. Accept that state, stop the powered phase, and pass the same time, position, velocity and mass into a new coast with zero thrust. Do not move the vehicle to apogee. A downward fuel-floor crossing m − (5,000 kg + payload mass) = 0, ground contact, or the 1,200 s burn limit can end the burn before the target gate.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Osculating energy, angular momentum and apsides** — `python/neutron/model.py:51`\n", "\n", "```python\n", "def orbital_elements(state):\n", " \"\"\"Osculating two-body apsides; a bound ellipse alone is not a safe orbit.\"\"\"\n", " radius, _, vr, vt, _ = map(float, state)\n", " energy = (vr * vr + vt * vt) / 2 - MU / radius\n", " momentum = radius * vt\n", " # Polar eccentricity-vector components avoid cancellation near a circle.\n", " eccentricity = math.hypot(radius * vt * vt / MU - 1, radius * vr * vt / MU)\n", " p = momentum**2 / MU\n", " return {\"perigee_km\": (p / (1 + eccentricity) - EARTH_RADIUS) / 1000,\n", " \"apogee_km\": ((p / (1 - eccentricity) - EARTH_RADIUS) / 1000\n", " if energy < 0 and eccentricity < 1 else None),\n", " \"eccentricity\": eccentricity, \"specific_energy_j_kg\": energy}\n", "```\n", "\n", "**Select transfer target, root gate and powered integration** — `python/neutron/model.py:489`\n", "\n", "```python\n", " transfer = target_altitude_km > TRANSFER_ALTITUDE_KM\n", " def transfer_apogee(t, y):\n", " apogee = orbital_elements(y)[\"apogee_km\"]\n", " return (apogee if apogee is not None else -1e6) - target_altitude_km\n", " transfer_gate = terminal_event(transfer_apogee, 1)\n", " upper_time, upper, upper_hit, solution = phase(\n", " \"upper_stage\", \"transfer_insertion\" if transfer else \"orbit_insertion\", upper, upper_time, upper_time + 1200,\n", " lambda t, y: upper_guidance(t, y, vehicle, min(target_altitude_km, TRANSFER_ALTITUDE_KM) * 1000),\n", " (empty, transfer_gate if transfer else target), 12.0)\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "upper-apogee", "metadata": {}, "source": [ "### 12 · Coast to the real high point\n", "\n", "With the engine off, propagate the cutoff state under the same gravity and drag equations. The vehicle is still climbing until its radial velocity changes sign.\n", "\n", "**Symbols in this section.**\n", "\n", "- $F_r,\\ F_t$ — Outward radial and positive tangential engine forces; both vanish during coast. **Units:** N.\n", "- $m,\\ \\dot m$ — Upper stage plus attached-payload mass, and its derivative. **Units:** kg; kg/s.\n", "- $v_r$ — Outward radial speed, equal to the derivative of altitude. At the high point it changes from climbing to descending. **Units:** m/s.\n", "- $\\downarrow$ — A positive-to-negative zero crossing of the event function; here it distinguishes apogee from a low point. **Units:** event direction.\n", "\n", "$$\n", "F_r=F_t=0,\\qquad\\dot m=0\n", "$$\n", "\n", "The upper stage coasts with reference area 12 m² and retains its cutoff mass.\n", "\n", "$$\n", "v_r=0\\quad(\\downarrow)\n", "$$\n", "\n", "Find the downward zero crossing of radial velocity: this is the propagated apogee event.\n", "\n", "**Computed reference flight.** Apogee occurs at 3270.066 s, after a coast of 2716.830 s.\n", "\n", "Accept the event's integrated time and state before switching guidance. At this event, restart the upper engine from the integrated apogee state. The restart changes force and mass-flow rate, not position, velocity or mass instantaneously. Restart is assumed instantaneous and perfect. An impact or the 6,000 s coast limit prevents this transition.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Coast until radial velocity crosses zero** — `python/neutron/model.py:501`\n", "\n", "```python\n", " apex = terminal_event(lambda t, y: y[2])\n", " upper_time, upper, upper_hit, solution = phase(\"upper_stage\", \"apogee_coast\", upper,\n", " upper_time, upper_time + 6000, coast, (apex,), 12.0, max(sample_step_s, 10.0))\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "upper-circularize", "metadata": {}, "source": [ "### 13 · Make a finite circularization burn\n", "\n", "Now set the altitude command $h_c$ to the final 400 km target and reuse the upper guidance. The burn takes finite time, consumes propellant, and continuously changes the state. Circular speed $v_c$ depends on the current radius.\n", "\n", "**Symbols in this section.**\n", "\n", "- $h_c,\\ h_{target}$ — Commanded altitude for upper guidance and final target altitude, both measured above the spherical surface. **Units:** m.\n", "- $r,\\ \\mu$ — Current geocentric radius and Earth gravitational parameter. **Units:** m; m³/s².\n", "- $v_c,\\ v_t$ — Local two-body circular speed and actual signed tangential speed. Their equality alone does not guarantee a circular orbit. **Units:** m/s.\n", "- $v_r,\\ \\dot v_r$ — Outward radial speed and its time derivative, used in the circular-balance derivation. **Units:** m/s; m/s².\n", "- $F_r,\\ D_r$ — Outward radial engine force and signed radial drag force. Both are zero in the ideal unpowered circular-balance derivation. **Units:** N.\n", "- $\\uparrow$ — The speed-error gate v_t − v_c crosses from negative to positive; v_c itself changes as radius changes. **Units:** event direction.\n", "\n", "$$\n", "h_c=h_{target},\\qquad v_c=\\sqrt{\\frac{\\mu}{r}}\n", "$$\n", "\n", "Use the new altitude command in the same radial request and force calculation.\n", "\n", "$$\n", "v_t-v_c=0\\quad(\\uparrow)\n", "$$\n", "\n", "Cut thrust at the upward crossing of circular speed, unless fuel depletion or ground contact occurs first.\n", "\n", "**Why these equations have this form.**\n", "\n", "$$\n", "\\begin{aligned}v_r&=\\dot v_r=0,\\quad F_r=D_r=0\\\\&\\Longrightarrow\\ \\frac{v_t^2}{r}=\\frac{\\mu}{r^2}\\end{aligned}\n", "$$\n", "\n", "The circular-speed formula comes from unpowered radial balance in the two-body limit. During a finite burn those assumptions do not all hold; this is why crossing circular speed is followed by independent orbit checks.\n", "\n", "**Computed reference flight.** The restart lasts 0.779141 s and ends at 3270.846 s. Tangential speed is 7668.558218 m/s.\n", "\n", "The speed crossing only ends the burn; the next checks decide whether insertion succeeded. For targets at or below 220 km, the first upper burn uses this speed gate directly and skips the transfer coast and restart. The circularization burn has the same 13,000 kg reference fuel floor and a 600 s time limit.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Circular-speed upward-crossing gate** — `python/neutron/model.py:485`\n", "\n", "```python\n", " target = terminal_event(lambda t, y: y[3] - math.sqrt(MU / y[0]), 1)\n", "```\n", "\n", "**Restart the engine at the target-altitude command** — `python/neutron/model.py:508`\n", "\n", "```python\n", " upper_time, upper, upper_hit, solution = phase(\"upper_stage\", \"circularization_burn\", upper,\n", " upper_time, upper_time + 600,\n", " lambda t, y: upper_guidance(t, y, vehicle, target_altitude_km * 1000), (empty, target), 12.0)\n", "```\n", "\n", "**The same finite-thrust upper-stage feedback** — `python/neutron/model.py:153`\n", "\n", "```python\n", "def upper_guidance(time, state, vehicle, target_m):\n", " \"\"\"PD radial acceleration request plus spherical-gravity compensation.\"\"\"\n", " radius, _, vr, vt, mass = state\n", " tau = 30.0\n", " requested = (EARTH_RADIUS + target_m - radius) / tau**2 - 2 * vr / tau\n", " radial_fraction = (requested + MU / radius**2 - vt * vt / radius) * mass / vehicle.upper_thrust_n\n", " beta = math.asin(np.clip(radial_fraction, -0.5, 0.95))\n", " return vehicle.upper_thrust_n * np.array([math.sin(beta), math.cos(beta)])\n", "```\n", "\n", "
\n" ] }, { "attachments": { "orbit-apsides.svg": { "image/svg+xml": [ "PHN2ZyB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHdpZHRoPSI3NjAiIGhlaWdodD0iMzQwIiB2aWV3Qm94PSIwIDAgNzYwIDM0MCIgcm9sZT0iaW1nIiBhcmlhLWxhYmVsbGVkYnk9InRpdGxlIGRlc2MiPgo8dGl0bGUgaWQ9InRpdGxlIj5PcmJpdGFsIHJhZGlpIHN0YXJ0IGF0IEVhcnRo4oCZcyBjZW50cmU8L3RpdGxlPgo8ZGVzYyBpZD0iZGVzYyI+QW4gZWxsaXBzZSB3aXRoIEVhcnRoIGF0IGl0cyBsZWZ0IGZvY3VzLCBwZXJpZ2VlIGF0IHRoZSBsZWZ0IGVuZCBhbmQgYXBvZ2VlIGF0IHRoZSByaWdodCBlbmQuIFRoZSBzZW1pLW1ham9yIGF4aXMgYSBydW5zIGZyb20gdGhlIGVsbGlwc2UgY2VudHJlIHRvIGFwb2dlZS4gUmFkaWkgYXJlIG1lYXN1cmVkIGZyb20gdGhlIGZvY3VzOyBhbHRpdHVkZXMgc3VidHJhY3QgRWFydGjigJlzIHJhZGl1cy48L2Rlc2M+CjxkZWZzPjxtYXJrZXIgaWQ9ImFycm93IiB2aWV3Qm94PSIwIDAgMTAgMTAiIHJlZlg9IjgiIHJlZlk9IjUiIG1hcmtlcldpZHRoPSI2IiBtYXJrZXJIZWlnaHQ9IjYiIG9yaWVudD0iYXV0by1zdGFydC1yZXZlcnNlIj48cGF0aCBkPSJNMCAwIDEwIDUgMCAxMFoiIGZpbGw9ImNvbnRleHQtc3Ryb2tlIi8+PC9tYXJrZXI+PC9kZWZzPgo8cmVjdCB3aWR0aD0iNzYwIiBoZWlnaHQ9IjM0MCIgcng9IjYiIGZpbGw9IiMxNDE0MTQiLz4KPGcgZm9udC1mYW1pbHk9IkFyaWFsLCBzYW5zLXNlcmlmIiBmb250LXNpemU9IjE1IiBmaWxsPSIjZjRmNGY0Ij4KPHRleHQgeD0iMjQiIHk9IjMyIiBmb250LXNpemU9IjE3IiBmb250LXdlaWdodD0iNjAwIj5SYWRpdXMgciBpcyBub3QgYWx0aXR1ZGUgaDwvdGV4dD4KPGVsbGlwc2UgY3g9IjM4MCIgY3k9IjE3NCIgcng9IjI3NiIgcnk9IjExMCIgZmlsbD0ibm9uZSIgc3Ryb2tlPSIjYjBiMGIwIiBzdHJva2Utd2lkdGg9IjEuNSIvPgo8Y2lyY2xlIGN4PSIxMjciIGN5PSIxNzQiIHI9IjE5IiBmaWxsPSIjMjYzMzNmIiBzdHJva2U9IiM4MjliYWQiLz4KPGNpcmNsZSBjeD0iMTI3IiBjeT0iMTc0IiByPSIzIiBmaWxsPSIjZGRkIi8+CjxwYXRoIGQ9Ik0xMjcgMTc0SDEwNiIgc3Ryb2tlPSIjZmY5MjlhIiBzdHJva2Utd2lkdGg9IjIiIG1hcmtlci1lbmQ9InVybCgjYXJyb3cpIi8+CjxwYXRoIGQ9Ik0xMjcgMTc0SDY1MCIgc3Ryb2tlPSIjYTdkNGVlIiBzdHJva2Utd2lkdGg9IjIiIG1hcmtlci1lbmQ9InVybCgjYXJyb3cpIi8+CjxwYXRoIGQ9Ik0xMjcgMTk5VjIyNEgxMDRWMTg3IiBmaWxsPSJub25lIiBzdHJva2U9IiNmZjkyOWEiLz4KPHRleHQgeD0iOTUiIHk9IjI0NSIgZmlsbD0iI2ZmOTI5YSI+cuKCmjwvdGV4dD4KPHRleHQgeD0iNDU1IiB5PSIxNjUiIGZpbGw9IiNhN2Q0ZWUiPnLigpAgPSBhKDEgKyBlKTwvdGV4dD4KPHRleHQgeD0iMjQiIHk9IjEzMSIgZmlsbD0iI2ZmOTI5YSI+UGVyaWdlZTwvdGV4dD48dGV4dCB4PSIyNCIgeT0iMTUxIiBmaWxsPSIjZmY5MjlhIj5y4oKaID0gYSgxIOKIkiBlKTwvdGV4dD4KPHRleHQgeD0iNjU1IiB5PSIxNTEiIGZpbGw9IiNhN2Q0ZWUiPkFwb2dlZTwvdGV4dD4KPGNpcmNsZSBjeD0iMTA0IiBjeT0iMTc0IiByPSI0IiBmaWxsPSIjZmY5MjlhIi8+PGNpcmNsZSBjeD0iNjU2IiBjeT0iMTc0IiByPSI0IiBmaWxsPSIjYTdkNGVlIi8+CjxwYXRoIGQ9Ik0zODAgMTc0VjIxMCBNNjU2IDE4OFYyMTAgTTM4MyAyMDhINjUzIiBmaWxsPSJub25lIiBzdHJva2U9IiNhYWEiIG1hcmtlci1lbmQ9InVybCgjYXJyb3cpIi8+Cjx0ZXh0IHg9IjQ0MCIgeT0iMjMyIj5hID0gc2VtaS1tYWpvciBheGlzPC90ZXh0Pjx0ZXh0IHg9IjI5NSIgeT0iMjUyIiBmaWxsPSIjYjBiMGIwIj5FbGxpcHNlIGNlbnRyZTwvdGV4dD48dGV4dCB4PSIxNDMiIHk9IjE0MiIgZmlsbD0iI2IwYjBiMCI+RWFydGggYXQgZm9jdXM8L3RleHQ+Cjx0ZXh0IHg9IjI0IiB5PSIzMTUiPmjigpogPSBy4oKaIOKIkiBSX0U8L3RleHQ+PHRleHQgeD0iMjYwIiB5PSIzMTUiPmjigpAgPSBy4oKQIOKIkiBSX0U8L3RleHQ+PHRleHQgeD0iNTE1IiB5PSIzMTUiIGZpbGw9IiNiMGIwYjAiPmU6IGVjY2VudHJpY2l0eSAodW5pdGxlc3MpPC90ZXh0Pgo8L2c+PC9zdmc+Cg==" ] } }, "cell_type": "markdown", "id": "upper-qualify", "metadata": {}, "source": [ "### 14 · Test the achieved orbit before releasing anything\n", "\n", "Compute perigee altitude $h_p$ from the cutoff state, using the same orbital quantities as before. The model requires a bound orbit, both apsides near the target, low radial speed, and low eccentricity together.\n", "\n", "**Symbols in this section.**\n", "\n", "- $p,\\ R_E$ — Semi-latus rectum from ℓ²/μ and the fixed spherical Earth radius. **Units:** m.\n", "- $e$ — Nonnegative eccentricity: zero is circular; a bound nondegenerate ellipse has e < 1. **Units:** dimensionless.\n", "- $\\varepsilon$ — Specific orbital mechanical energy; negative is required for a bound two-body orbit. **Units:** J/kg.\n", "- $h_p,\\ h_a,\\ h_{target}$ — Osculating perigee altitude, osculating apogee altitude, and final target altitude. The displayed 2 km tolerance is converted to 2,000 m when comparing SI quantities. **Units:** m.\n", "- $v_r$ — Actual outward radial speed at cutoff; a circular orbit would have zero radial speed. **Units:** m/s.\n", "\n", "$$\n", "h_p=\\frac{p}{1+e}-R_E\n", "$$\n", "\n", "Perigee is the low point of the orbit implied by the current state.\n", "\n", "$$\n", "\\varepsilon<0,\\qquad e\\leq0.001\n", "$$\n", "\n", "The orbit must be bound and nearly circular.\n", "\n", "$$\n", "|h_p-h_{target}|\\leq2\\ \\mathrm{km}\n", "$$\n", "\n", "Perigee must be within the altitude tolerance.\n", "\n", "$$\n", "|h_a-h_{target}|\\leq2\\ \\mathrm{km}\n", "$$\n", "\n", "Apogee must independently be within the same tolerance.\n", "\n", "$$\n", "|v_r|\\leq5\\ \\mathrm{m/s}\n", "$$\n", "\n", "A circular-speed crossing at the wrong altitude or with large radial motion is a failure.\n", "\n", "**Computed reference flight.** Achieved apsides: 399.999868 × 399.999981 km; radial speed 0.00006380 m/s; eccentricity 8.320e-09.\n", "\n", "Success also requires reaching the intended speed gate without impact. Failed insertion is reported as fuel depletion, missed orbit, impact, or timeout; it is never repaired by changing the state. A surviving failed insertion coasts for up to 1,200 s without release.\n", "\n", "\n", "![Perigee radius rₚ and apogee radius rₐ are distances from Earth's centre (m); subtract Earth radius R_E to obtain altitudes hₚ and hₐ (m). Semi-major axis a is half the ellipse's long diameter (m); e is dimensionless eccentricity. Static geometry illustration; not a flight result.](attachment:orbit-apsides.svg)\n", "\n", "*Perigee radius rₚ and apogee radius rₐ are distances from Earth's centre (m); subtract Earth radius R_E to obtain altitudes hₚ and hₐ (m). Semi-major axis a is half the ellipse's long diameter (m); e is dimensionless eccentricity. Static geometry illustration; not a flight result.*\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**All orbit tolerances must pass together** — `python/neutron/model.py:65`\n", "\n", "```python\n", "def orbit_target_error(state, target_altitude_km):\n", " \"\"\"Enter the target set only when BOTH apsides and radial speed agree.\n", "\n", " A proper target orbit is bound, within 2 km at each apsis, |vr| <= 5 m/s,\n", " and e <= 0.001. A negative result satisfies all four constraints.\n", " \"\"\"\n", " orbit = orbital_elements(state)\n", " if orbit[\"apogee_km\"] is None or orbit[\"specific_energy_j_kg\"] >= 0:\n", " return 1e6\n", " return max(abs(orbit[\"perigee_km\"] - target_altitude_km) / 2,\n", " abs(orbit[\"apogee_km\"] - target_altitude_km) / 2,\n", " abs(state[2]) / 5, orbit[\"eccentricity\"] / 0.001) - 1\n", "```\n", "\n", "**Require the speed gate and no impact before declaring success** — `python/neutron/model.py:512`\n", "\n", "```python\n", " elements = orbital_elements(upper)\n", " reached = bool(cutoff_triggered and not upper_hit\n", " and orbit_target_error(upper, target_altitude_km) <= 1e-8)\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "payload-release", "metadata": {}, "source": [ "### 15 · Coast ten seconds, then release the payload\n", "\n", "After verified insertion at the upper-stage cutoff time $t_{cutoff}$, coast for ten seconds before release. Apply the same momentum-conserving separation rule, now with relative radial speed $\\delta v$ = 0.2 m/s. Let $m_p$ be payload mass and $m_s$ the retained upper-stage mass; $m^-$ is their combined mass just before release.\n", "\n", "**Symbols in this section.**\n", "\n", "- $t_{release},\\ t_{cutoff}$ — Payload-release and verified upper-stage cutoff times on the mission clock. **Units:** s.\n", "- $m^-,\\ m_s,\\ m_p$ — Combined mass immediately before release, retained upper-stage mass immediately after release, and payload mass. Superscript − means before the event, not negative mass. **Units:** kg.\n", "- $v_r^-,\\ v_{r,s},\\ v_{r,p}$ — Common outward radial speed immediately before release, and retained-stage and payload radial speeds immediately after release. **Units:** m/s.\n", "- $\\delta v$ — Specified relative outward separation speed v_r,p − v_r,s = 0.2 m/s. It is the difference between the two final speeds, not the increment given to each body. **Units:** m/s.\n", "\n", "$$\n", "t_{release}=t_{cutoff}+10\\ \\mathrm s\n", "$$\n", "\n", "The deployment coast uses zero thrust; impact during this coast blocks release.\n", "\n", "$$\n", "m_s=m^- -m_p,\\qquad m_p=8\\,000\\ \\mathrm{kg}\n", "$$\n", "\n", "Only the retained body loses the payload mass; the total mass is conserved.\n", "\n", "$$\n", "v_{r,s}=v_r^- -\\delta v\\frac{m_p}{m^-}\n", "$$\n", "\n", "The retained upper stage receives the equal-and-opposite recoil.\n", "\n", "$$\n", "v_{r,p}=v_r^- +\\delta v\\frac{m_s}{m^-}\n", "$$\n", "\n", "The payload receives its share of the separation impulse. Position and tangential velocity are unchanged.\n", "\n", "**Why these equations have this form.**\n", "\n", "$$\n", "\\begin{aligned}m_s v_{r,s}+m_p v_{r,p}&=m^-v_r^-\\\\v_{r,p}-v_{r,s}&=\\delta v\\end{aligned}\n", "$$\n", "\n", "Solve conservation of radial linear momentum together with the prescribed relative speed to obtain the two recoil formulas. This is an instantaneous internal impulse at a shared location. It conserves total momentum, but separation supplies a small amount of kinetic energy; no finite mechanism or contact force is modeled.\n", "\n", "**Computed reference flight.** Release occurs at 3280.846 s. Payload mass is 8,000 kg; retained upper-stage mass is 5295.957 kg.\n", "\n", "First finish the ten-second coast and accept its endpoint. At that same mission time, split the accepted mass and apply the equal-and-opposite velocity increments; then start one new integration per body. Reference area is 12 m² for the retained upper stage and 2 m² for the payload. The point masses begin at an identical location; the model assumes successful physical clearance.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Coast, check for impact, then separate the payload** — `python/neutron/model.py:524`\n", "\n", "```python\n", " coast_name, duration = (\"deployment_coast\", 10) if reached else (\"failed_insertion_coast\", 1200)\n", " upper_time, upper, upper_hit, _ = phase(\"upper_stage\", coast_name, upper,\n", " upper_time, upper_time + duration, coast, body_area=12.0)\n", " if upper_hit:\n", " record(\"upper_impact\", upper_time, \"upper_stage\", \"Upper stage contacted the ground during unpowered flight.\")\n", " if reached and not upper_hit:\n", " retained, payload = separate(upper, upper[4] - payload_kg, 0.2)\n", "```\n", "\n", "**The same momentum-conserving split with 0.2 m/s relative speed** — `python/neutron/model.py:198`\n", "\n", "```python\n", "def separate(state, booster_mass_kg, relative_speed_m_s=0.5):\n", " \"\"\"Split mass with equal/opposite impulses: conserve linear momentum.\n", "\n", " Both bodies start at the same position. Geometry and contact are not modeled.\n", " The reusable fairing stays within booster dry mass.\n", " \"\"\"\n", " booster, upper = state.copy(), state.copy()\n", " upper[4] = state[4] - booster_mass_kg\n", " booster[4] = booster_mass_kg\n", " booster[2] -= relative_speed_m_s * upper[4] / state[4]\n", " upper[2] += relative_speed_m_s * booster[4] / state[4]\n", " return booster, upper\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "payload-orbit", "metadata": {}, "source": [ "### 16 · Propagate a whole unpowered revolution\n", "\n", "Use the released payload's specific energy $\\varepsilon_p$ to estimate semi-major axis $a_p$ and orbital period $P$. Then integrate for that duration; the estimated period does not substitute an analytic ellipse for the actual trajectory.\n", "\n", "**Symbols in this section.**\n", "\n", "- $\\mu$ — Earth gravitational parameter, 3.986004418 × 10¹⁴. **Units:** m³/s².\n", "- $\\varepsilon_p,\\ \\ell_p$ — Payload specific energy and signed specific angular momentum immediately after release; the reference values for the drift tests. **Units:** J/kg; m²/s.\n", "- $a_p,\\ P$ — Payload's osculating two-body semi-major axis and corresponding estimated orbital period. The subscript p denotes payload here, whereas $h_p$ denotes perigee. **Units:** m; s.\n", "- $t_{release},\\ t_{req},\\ t_{end}$ — Release time, requested integration endpoint, and actual endpoint returned by the solver, all on the mission clock. **Units:** s.\n", "- $\\theta,\\ \\Delta\\theta$ — Unwrapped Earth-centred angle, positive along the ascent direction, and its total change since release. The angle is not reset at 2π. **Units:** rad.\n", "- $n,\\ t_n$ — Index of a returned solver state, and that state's mission time. These are accepted integration endpoints and the initial state, not the coarser playback frames. **Units:** integer index; s.\n", "- $h(t_n),\\ h_{target}$ — Payload altitude r(t_n) − R_E at each checked state, and target altitude; the 2 km bound means 2,000 m. **Units:** m.\n", "- $\\varepsilon(t_n),\\ \\ell(t_n)$ — Specific energy (v_r² + v_t²)/2 − μ/r and signed specific angular momentum r v_t evaluated at each checked state. **Units:** J/kg; m²/s.\n", "- $\\max_n$ — Largest value over all returned solver states. Relative drift divides by the nonzero release value, making it dimensionless. **Units:** operator.\n", "\n", "$$\n", "a_p=-\\frac{\\mu}{2\\varepsilon_p}\n", "$$\n", "\n", "Evaluate payload energy from its post-release state, including the small separation impulse.\n", "\n", "$$\n", "P=2\\pi\\sqrt{\\frac{a_p^3}{\\mu}}\n", "$$\n", "\n", "This is the requested unpowered propagation duration for both released bodies.\n", "\n", "$$\n", "\\begin{aligned}t_{req}&=t_{release}+P\\\\|t_{end}-t_{req}|&<10^{-3}\\ \\mathrm s\\end{aligned}\n", "$$\n", "\n", "Propagate with zero thrust until requested time $t_{req}$ or ground contact. The actual endpoint $t_{end}$ must reach the requested time within 0.001 s, with no impact.\n", "\n", "$$\n", "\\Delta\\theta\\geq2\\pi-10^{-3}\n", "$$\n", "\n", "The swept angle is $\\Delta\\theta=\\theta(t_{end})-\\theta(t_{release})$; it must also cover a complete revolution within this angular tolerance.\n", "\n", "$$\n", "\\begin{aligned}|h(t_n)-h_{target}|&\\leq2\\ \\mathrm{km}\\\\\\varepsilon(t_n)&<0\\end{aligned}\n", "$$\n", "\n", "Check every returned solver state, indexed by $n$: the initial state, accepted step endpoints and any final event endpoint. All must remain near the target altitude and bound.\n", "\n", "$$\n", "\\max_n\\left|\\frac{\\varepsilon(t_n)-\\varepsilon_p}{\\varepsilon_p}\\right|<10^{-5}\n", "$$\n", "\n", "The maximum relative energy change over these solver states must remain small during the unpowered orbit.\n", "\n", "$$\n", "\\max_n\\left|\\frac{\\ell(t_n)-\\ell_p}{\\ell_p}\\right|<10^{-5}\n", "$$\n", "\n", "The same relative bound applies to specific angular momentum, measured from its release value $\\ell_p$. These are discrete checks, not continuous extrema.\n", "\n", "**Computed reference flight.** The propagated orbit lasts 5553.624 s and sweeps 6.283185307 rad. Altitude spans 399.929456–400.070394 km; energy drift is 6.208e-15.\n", "\n", "The period formula assumes an isolated bound two-body ellipse. The integrated coast retains the still-atmosphere drag model, so energy and angular momentum need only remain within the stated small drift thresholds, rather than being exactly conserved. Earth oblateness, Earth rotation, third bodies and orbital-plane dynamics are absent. The payload-orbit result comes from the integrated solution. This completes the upper branch. Next, return to separation and follow the booster on the shared mission clock.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Calculate period and propagate both unpowered bodies** — `python/neutron/model.py:535`\n", "\n", "```python\n", " energy = orbital_elements(payload)[\"specific_energy_j_kg\"]\n", " semimajor_axis = -MU / (2 * energy)\n", " period = 2 * math.pi * math.sqrt(semimajor_axis**3 / MU)\n", " # Keep long-coast exports compact without changing solver steps.\n", " interval = max(sample_step_s, 10.0)\n", " phase(\"upper_stage\", \"post_release_coast\", retained, upper_time, upper_time + period,\n", " coast, body_area=12.0, sample_interval=interval)\n", " payload_time, _, _, solution = phase(\"payload\", \"free_flight\", payload,\n", " upper_time, upper_time + period, coast, body_area=2.0, sample_interval=interval)\n", " payload_orbit = coast_verification(solution, period, target_altitude_km)\n", "```\n", "\n", "**Check returned solver states, duration and swept angle** — `python/neutron/model.py:79`\n", "\n", "```python\n", "def coast_verification(solution, period_s, target_altitude_km):\n", " \"\"\"Verify an actual propagated revolution; no analytic orbit is substituted.\"\"\"\n", " radius, angle, vr, vt, _ = solution.y\n", " altitude_km = (radius - EARTH_RADIUS) / 1000\n", " energy = (vr**2 + vt**2) / 2 - MU / radius\n", " momentum = radius * vt\n", " energy_change = float(np.max(np.abs((energy - energy[0]) / energy[0])))\n", " momentum_change = float(np.max(np.abs((momentum - momentum[0]) / momentum[0])))\n", " duration = float(solution.t[-1] - solution.t[0])\n", " swept_angle = float(angle[-1] - angle[0])\n", " completed = abs(duration - period_s) < 1e-3 and not solution.t_events[0].size\n", " minimum, maximum = float(altitude_km.min()), float(altitude_km.max())\n", " initial_orbit = orbital_elements(solution.y[:, 0])\n", " success = (completed and np.all(energy < 0) and swept_angle >= 2 * math.pi - 1e-3\n", " and minimum >= target_altitude_km - 2 and maximum <= target_altitude_km + 2\n", " and energy_change < 1e-5 and momentum_change < 1e-5)\n", " return dict(success=bool(success), completed_one_orbit=bool(completed),\n", " period_s=period_s, propagated_duration_s=duration, swept_angle_rad=swept_angle,\n", " perigee_km=initial_orbit[\"perigee_km\"], apogee_km=initial_orbit[\"apogee_km\"],\n", " eccentricity=initial_orbit[\"eccentricity\"],\n", " min_altitude_km=minimum, max_altitude_km=maximum,\n", " relative_energy_change=energy_change, relative_angular_momentum_change=momentum_change)\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "branch-booster", "metadata": {}, "source": [ "## Booster → return and landing\n", "\n", "Carry the booster state forward from the same separation epoch through closure, boostback, descent and feedback landing." ] }, { "attachments": {}, "cell_type": "markdown", "id": "booster-close", "metadata": {}, "source": [ "### 17 · Coast while the upper stage clears the fairing\n", "\n", "Return to the booster state at stage separation, not the end of upper-stage flight. The booster retains its 90,000 kg mass, including 60,000 kg recovery propellant and the fairing. It coasts for four seconds before boostback. The fairing remains open until integrated upper-stage displacement along the release axis reaches 23.2 m, then closes over four seconds.\n", "\n", "**Symbols in this section.**\n", "\n", "- $t,\\ t_{sep},\\ t_{bb},\\ t_{close}$ — Current mission time, shared stage separation, boostback ignition, and clearance-gated fairing closure times. **Units:** s.\n", "- $f_{open}$ — Prescribed fairing opening fraction: 1 means fully open and 0 fully closed. clip limits the result to this interval. **Units:** dimensionless.\n", "- $F_r,\\ F_t$ — Outward radial and positive tangential engine forces; they are zero during the initial separation coast. **Units:** N.\n", "- $m,\\ \\dot m$ — Booster mass including captive fairing and remaining recovery propellant, and its time derivative. **Units:** kg; kg/s.\n", "\n", "$$\n", "t_{bb}=t_{sep}+4\\ \\mathrm s\n", "$$\n", "\n", "Boostback can start at $t_{bb}$, four seconds after separation.\n", "\n", "$$\n", "f_{open}=\\operatorname{clip}\\!\\left(1-\\frac{t-t_{close}}{4\\ \\mathrm s},0,1\\right)\n", "$$\n", "\n", "The clearance root is found from the continuous integrated trajectories. Closure is kinematic: it neither ejects mass nor changes the aerodynamic reference area.\n", "\n", "$$\n", "F_r=F_t=0,\\qquad\\dot m=0\n", "$$\n", "\n", "Advance the booster under gravity and drag with the established 7 m diameter reference area.\n", "\n", "**Computed reference flight.** Boostback begins at 145.979 s. The fairing starts closing at 146.936 s and is closed at 150.936 s. The upper engine started one second earlier, at 144.979 s.\n", "\n", "The branches are physically concurrent. Each starts from its own post-separation state at $t_{sep}$; each force law advances its own state. The final fairing calculation compares both integrated positions at the same mission time. It does not feed back into either force law, even though Python evaluates one branch after the other.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Restart the booster at the shared separation time** — `python/neutron/model.py:549`\n", "\n", "```python\n", " bt, booster, booster_hit, _ = phase(\"booster\", \"clearance_coast\", booster, time, time + 4, coast)\n", "```\n", "\n", "**Initial fairing samples before the clearance-gated hardware pass** — `python/neutron/model.py:379`\n", "\n", "```python\n", " def fairing_fraction(time, body):\n", " opening = next((e[\"time_s\"] for e in events if e[\"name\"] == \"fairing_opening\"), None)\n", " if opening is None or body not in (\"stack\", \"booster\"):\n", " return 0.0\n", " if separation_time is None or time <= separation_time:\n", " return float(np.clip((time - opening) / 4, 0, 1))\n", " return float(np.clip(1 - (time - separation_time) / 4, 0, 1))\n", "```\n", "\n", "**Final dense-state clearance gate, closure and coast pointing** — `python/neutron/model.py:216`\n", "\n", "```python\n", "def apply_kinematic_hardware(trajectories, events, dense_phases):\n", " \"\"\"Prescribed hardware coordinates, sampled from the actual integrated flight.\n", "\n", " Coast pointing uses a bounded cubic pre-alignment before the next ignition.\n", " Powered pointing is the actual force axis. Fairing closure waits for physical\n", " upper-stage clearance. No position, velocity, mass or applied force changes.\n", " \"\"\"\n", " from scipy.optimize import brentq\n", "\n", " def phase_at(body, time):\n", " return next((phase for phase in reversed(dense_phases[body])\n", " if phase[\"start\"] <= time <= phase[\"stop\"]), None)\n", "\n", " def state_xy(body, time):\n", " state = phase_at(body, time)[\"solution\"].sol(time)\n", " return state[0] * np.array([math.cos(state[1]), math.sin(state[1])])\n", "\n", " def insert_samples(body, times):\n", " rows = trajectories[body]\n", " for time in times:\n", " if any(abs(row[\"time_s\"] - time) < 1e-8 for row in rows):\n", " continue\n", " phase = phase_at(body, time)\n", " if phase is None:\n", " continue\n", " state = phase[\"solution\"].sol(time)\n", " radius, angle, vr, vt, mass = map(float, state)\n", " force = phase[\"guidance\"](time, state)\n", " thrust = float(np.linalg.norm(force))\n", " ca, sa = math.cos(angle), math.sin(angle)\n", " fx, fy = float(force[0] * ca - force[1] * sa), float(force[0] * sa + force[1] * ca)\n", " before = max((row for row in rows if row[\"time_s\"] <= time), key=lambda row: row[\"time_s\"])\n", " row = dict(before, time_s=float(time), phase=phase[\"name\"],\n", " altitude_m=radius - EARTH_RADIUS, downrange_m=EARTH_RADIUS * angle,\n", " radial_velocity_m_s=vr, tangential_velocity_m_s=vt,\n", " speed_m_s=math.hypot(vr, vt), mass_kg=mass, thrust_n=thrust,\n", " force_x_n=fx, force_y_n=fy,\n", " dynamic_pressure_pa=0.5 * atmosphere(radius - EARTH_RADIUS) * (vr * vr + vt * vt),\n", " x_m=radius * ca, y_m=radius * sa,\n", " vx_m_s=vr * ca - vt * sa, vy_m_s=vr * sa + vt * ca)\n", " if thrust > 0:\n", " row[\"axis_x\"], row[\"axis_y\"] = fx / thrust, fy / thrust\n", " if body == \"booster\" and any(sample.get(\"leg_open_fraction\", 0) > 0 for sample in rows):\n", " command = next((event[\"time_s\"] for event in events if event[\"name\"] == \"landing_burn\"), math.inf)\n", " row[\"leg_open_fraction\"] = float(np.clip((time - command) / 2, 0, 1))\n", " rows.append(row)\n", " rows.sort(key=lambda row: row[\"time_s\"])\n", "\n", " for body, rows in trajectories.items():\n", " segments = []\n", " index = 0\n", " while index < len(rows):\n", " if rows[index][\"thrust_n\"] > 0:\n", " index += 1\n", " continue\n", " start = index\n", " while index < len(rows) and rows[index][\"thrust_n\"] == 0:\n", " index += 1\n", " if index == len(rows):\n", " break\n", " first, target = rows[start], rows[index]\n", " initial = math.atan2(first[\"axis_y\"], first[\"axis_x\"])\n", " final = math.atan2(target[\"force_y_n\"], target[\"force_x_n\"])\n", " turn = math.atan2(math.sin(final - initial), math.cos(final - initial))\n", " available = target[\"time_s\"] - first[\"time_s\"]\n", " duration = min(available, max(2.0, 1.5 * abs(turn) / COAST_POINTING_RATE_RAD_S))\n", " if duration > 0:\n", " segments.append((target[\"time_s\"] - duration, target[\"time_s\"], initial, turn))\n", " # Include the profile endpoints and 0.1 s kinematic samples even when\n", " # the caller chooses sparse trajectory output. States use RK45 dense\n", " # output here; they are never linearly invented or re-integrated.\n", " insert_samples(body, [float(t) for begin, end, _, _ in segments\n", " for t in np.r_[np.arange(begin, end, 0.1), end]])\n", " for begin, end, initial, turn in segments:\n", " for row in rows:\n", " if row[\"thrust_n\"] == 0 and begin <= row[\"time_s\"] <= end:\n", " u = float(np.clip((row[\"time_s\"] - begin) / (end - begin), 0, 1))\n", " angle = initial + turn * u * u * (3 - 2 * u)\n", " row[\"axis_x\"], row[\"axis_y\"] = math.cos(angle), math.sin(angle)\n", "\n", " opening = next((event[\"time_s\"] for event in events if event[\"name\"] == \"fairing_opening\"), None)\n", " separation = next((event[\"time_s\"] for event in events if event[\"name\"] == \"stage_separation\"), None)\n", " if opening is None or separation is None or not trajectories[\"booster\"] or not trajectories[\"upper_stage\"]:\n", " return\n", " attached = trajectories[\"stack\"][-1]\n", " axis = np.array([attached[\"axis_x\"], attached[\"axis_y\"]])\n", " def clearance(time):\n", " return float((state_xy(\"upper_stage\", time) - state_xy(\"booster\", time)) @ axis) - UPPER_EXIT_CLEARANCE_M\n", " stop = min(dense_phases[body][-1][\"stop\"] for body in (\"booster\", \"upper_stage\"))\n", " times = sorted(set(float(t) for body in (\"booster\", \"upper_stage\")\n", " for phase in dense_phases[body] for t in phase[\"solution\"].t\n", " if separation <= t <= stop))\n", " close = math.inf\n", " for left, right in zip(times, times[1:]):\n", " if clearance(left) <= 0 <= clearance(right):\n", " close = float(brentq(clearance, left, right, xtol=1e-9))\n", " break\n", " events[:] = [event for event in events if event[\"name\"] not in (\"fairing_closing\", \"fairing_closed\")]\n", " if math.isfinite(close):\n", " for name, time, text in [\n", " (\"fairing_closing\", close, \"Captive fairing closes after integrated upper-stage clearance.\"),\n", " (\"fairing_closed\", close + 4, \"Four-second closure completes after upper-stage clearance.\")]:\n", " events.append(dict(name=name, time_s=time, body=\"booster\", description=text,\n", " required_axial_clearance_m=UPPER_EXIT_CLEARANCE_M))\n", " insert_samples(\"booster\", [close, close + 4])\n", " for body in (\"stack\", \"booster\"):\n", " for row in trajectories[body]:\n", " time = row[\"time_s\"]\n", " fraction = min(1.0, max(0.0, (time - opening) / 4))\n", " if time > close:\n", " fraction = max(0.0, 1 - (time - close) / 4)\n", " row[\"fairing_open_fraction\"] = fraction\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "booster-boostback", "metadata": {}, "source": [ "### 18 · Reverse downrange motion, then cut off on predicted miss\n", "\n", "Assume three engines point directly against the positive tangential direction. Let $T_e$ be one engine's thrust, one ninth of the published nine-engine ascent thrust. For cutoff only, estimate remaining fall time $t_f$ under constant local gravity g and the current vertical state. The estimated ground miss $x_{pred}$ is current downrange plus tangential speed times this fall time.\n", "\n", "**Symbols in this section.**\n", "\n", "- $T_e,\\ F_r,\\ F_t$ — One engine's assumed available thrust, outward radial engine force, and signed tangential engine force. Negative F_t points against the positive downrange direction. **Units:** N.\n", "- $r,\\ R_E,\\ h$ — Vehicle geocentric radius, fixed Earth radius, and altitude h = r − R_E. **Units:** m.\n", "- $\\mu,\\ g$ — Earth gravitational parameter and current positive local gravity magnitude μ/r². The predictor holds this g constant during its imagined fall. **Units:** m³/s²; m/s².\n", "- $v_r,\\ v_t$ — Actual radial and tangential speeds at the predictor's starting state; v_r is negative while descending. **Units:** m/s.\n", "- $t_f$ — Nonnegative estimated time remaining until ground contact in the constant-gravity, drag-free predictor; it is a duration, not an absolute mission time. **Units:** s.\n", "- $\\theta$ — Signed, unwrapped angle from the launch-site radius; positive toward the original ascent direction. **Units:** rad.\n", "- $x_{pred},\\ x$ — Predicted signed miss at contact and actual current signed surface-arc distance x = R_E $\\theta$ from the launch site. **Units:** m.\n", "- $\\downarrow$ — Predicted miss crossing zero from positive to negative; this gate ends boostback. **Units:** event direction.\n", "\n", "$$\n", "T_e=716\\,657.927\\ \\mathrm N\n", "$$\n", "\n", "One engine's thrust is the nine-engine ascent thrust divided by nine.\n", "\n", "$$\n", "F_r=0,\\qquad F_t=-3T_e\n", "$$\n", "\n", "Feed these forces into the same ODE with booster Isp = 330 s; this burn also consumes the recovery reserve.\n", "\n", "$$\n", "\\begin{aligned}g&=\\frac{\\mu}{r^2}\\\\t_f&=\\frac{v_r+\\sqrt{v_r^2+2g\\max(0,h)}}{g}\\end{aligned}\n", "$$\n", "\n", "The predictor assumes constant gravity and no drag. It only decides when to stop the burn.\n", "\n", "$$\n", "x_{pred}=R_E\\theta+v_t t_f\n", "$$\n", "\n", "Predicted miss is signed downrange distance from the launch site. This local approximation uses tangential speed as surface-distance rate; exactly, $\\dot x=(R_E/r)v_t$ for $x=R_E\\theta$. It also ignores curvature during the predicted fall.\n", "\n", "$$\n", "x_{pred}=0\\quad(\\downarrow)\n", "$$\n", "\n", "Cut off at the downward zero crossing, or earlier at the dry-mass fuel floor or ground contact.\n", "\n", "**Why these equations have this form.**\n", "\n", "$$\n", "0=\\max(0,h)+v_r t_f-\\tfrac12g t_f^2\n", "$$\n", "\n", "This constant-acceleration height equation is the predictor's model. The nonnegative root of this quadratic gives the stated fall-time expression; the main trajectory continues to use changing spherical gravity, drag and tangential motion.\n", "\n", "**Computed reference flight.** Boostback ends at 194.050 s after 48.071 s. Tangential velocity is then -338.563 m/s.\n", "\n", "The actual path is still integrated with spherical gravity and drag; the predictor never overwrites position. This educational return-to-launch-site scenario differs from Rocket Lab's published downrange sea-recovery baseline. The boostback time limit is 300 s.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Approximate fall time and signed miss for cutoff only** — `python/neutron/model.py:163`\n", "\n", "```python\n", "def predicted_miss(state, target_downrange_m=0.0, contact_height_m=0.0):\n", " \"\"\"Cheap constant-g ballistic predictor used only to stop boostback.\n", "\n", " It does not set a landing position; the integrated trajectory can miss.\n", " Landing feedback subsequently corrects remaining position/velocity errors.\n", " \"\"\"\n", " radius, angle, vr, vt, _ = state\n", " gravity = MU / radius**2\n", " fall_time = (vr + math.sqrt(vr * vr + 2 * gravity * max(0, radius - EARTH_RADIUS - contact_height_m))) / gravity\n", " return EARTH_RADIUS * angle + vt * fall_time - target_downrange_m\n", "```\n", "\n", "**Three-engine reverse thrust and stopping events** — `python/neutron/model.py:550`\n", "\n", "```python\n", " engine = vehicle.booster_thrust_n / vehicle.booster_engines\n", " empty = terminal_event(lambda t, y: y[4] - vehicle.booster_dry_kg)\n", " initial_miss = predicted_miss(booster, recovery_target_downrange_m, contact_height)\n", " correction = -1.0 if initial_miss >= 0 else 1.0\n", " return_target = terminal_event(\n", " lambda t, y: predicted_miss(y, recovery_target_downrange_m, contact_height), correction)\n", " if not booster_hit:\n", " record(\"fairing_closed\", bt, \"booster\", \"Captive fairing closed for recovery.\")\n", " if not booster_hit and booster[4] > vehicle.booster_dry_kg + 1e-3:\n", " record(\"boostback_ignition\", bt, \"booster\",\n", " \"Three-engine range correction targets the fixed recovery ship.\" if ship_recovery\n", " else \"Three-engine boostback targets the launch site (assumed RTLS scenario).\")\n", " bt, booster, booster_hit, solution = phase(\"booster\", \"boostback\", booster, bt, bt + 300,\n", " lambda t, y: np.array([0.0, correction * 3 * engine]), (empty, return_target))\n", " if not booster_hit:\n", " reason = \"predicted zero miss\" if solution.t_events[2].size else \"propellant depletion\" if solution.t_events[1].size else \"time limit\"\n", " record(\"boostback_cutoff\", bt, \"booster\", \"Boostback stops at \" + reason + \".\")\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "booster-gate", "metadata": {}, "source": [ "### 19 · Coast to the landing-burn ignition point\n", "\n", "Coast with zero thrust after boostback. At each state, estimate maximum net upward acceleration $a_{max}$ from three engines, then compute stopping distance $d_{stop}$ using only downward radial speed.\n", "\n", "**Symbols in this section.**\n", "\n", "- $T_e,\\ m$ — One engine's available thrust and current booster mass, including remaining recovery propellant. **Units:** N; kg.\n", "- $\\mu,\\ r$ — Earth gravitational parameter and current geocentric radius; μ/r² is the local gravitational acceleration magnitude. **Units:** m³/s²; m.\n", "- $h,\\ v_r$ — Current altitude r − R_E and outward radial speed; only negative v_r contributes to the stopping-distance estimate. **Units:** m; m/s.\n", "- $a_{max}$ — Estimated upward net acceleration for three engines, floored at 1 m/s². This heuristic estimate is not a guaranteed available deceleration. **Units:** m/s².\n", "- $d_{stop}$ — Estimated vertical stopping distance at constant $a_{max}$; excludes tangential speed and lateral thrust demand. **Units:** m.\n", "- $\\downarrow$ — The altitude-margin gate crossing zero from positive to negative; it triggers landing ignition. **Units:** event direction.\n", "\n", "$$\n", "a_{max}=\\max\\!\\left(1\\ \\mathrm{m/s^2},\\frac{3T_e}{m}-\\frac{\\mu}{r^2}\\right)\n", "$$\n", "\n", "The acceleration estimate uses current mass and gravity, with a 1 m/s² lower bound.\n", "\n", "$$\n", "d_{stop}=\\frac{\\min(v_r,0)^2}{2a_{max}}\n", "$$\n", "\n", "Only descending speed contributes to this approximate stopping-distance gate.\n", "\n", "$$\n", "h-1.5d_{stop}-3\\,000\\ \\mathrm m=0\\quad(\\downarrow)\n", "$$\n", "\n", "Ignite at the downward crossing, leaving a 50% stopping-distance margin plus 3 km.\n", "\n", "**Why these equations have this form.**\n", "\n", "$$\n", "0=\\min(v_r,0)^2-2a_{max}d_{stop}\n", "$$\n", "\n", "The stopping-distance estimate is the familiar squared-speed relation for constant deceleration. Actual available vertical deceleration changes as fuel burns and thrust is redirected laterally; drag and polar geometry also act. The factor 1.5 and 3 km offset are chosen ignition margins, not guarantees of a safe touchdown.\n", "\n", "**Computed reference flight.** The landing burn begins at 368.023 s and 39.640 km altitude. At this state $a_{max}$ = 27.350 m/s² and $d_{stop}$ = 24.426 km.\n", "\n", "Descending through 80 km is recorded at 326.230 s; drag has acted continuously all along. This marker does not switch the atmosphere on. Ground impact remains terminal, and this coast is limited to 1,000 s.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Stopping-distance ignition gate** — `python/neutron/model.py:569`\n", "\n", "```python\n", " def landing_gate(t, y):\n", " height, vr, mass = y[0] - EARTH_RADIUS - contact_height, y[2], y[4]\n", " maximum_deceleration = max(1.0, 3 * engine / mass - MU / y[0]**2)\n", " stopping_distance = min(vr, 0)**2 / (2 * maximum_deceleration)\n", " return height - 1.5 * stopping_distance - 3000\n", "```\n", "\n", "**Coast with ignition gate and nonterminal reentry marker** — `python/neutron/model.py:576`\n", "\n", "```python\n", " gate = terminal_event(landing_gate)\n", " reentry = lambda t, y: y[0] - EARTH_RADIUS - 80_000\n", " reentry.direction = -1\n", " bt, booster, booster_hit, solution = phase(\"booster\", \"coast_and_reentry\", booster, bt, bt + 1000,\n", " coast, (gate, reentry))\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "booster-land", "metadata": {}, "source": [ "### 20 · Use the current descent error to command thrust\n", "\n", "First set a desired downward speed $v_r^*$ from height, clamped to $h_+=\\max(0,h)$. Then request radial acceleration $a_r^*$ to follow that speed. Use the estimated time remaining $\\tau_l$ to drive both downrange $x=R_E\\theta$ and tangential speed toward zero.\n", "\n", "**Symbols in this section.**\n", "\n", "- $r,\\ R_E,\\ h,\\ h_+$ — Geocentric radius, fixed Earth radius, altitude r − R_E, and altitude clamped below at zero. **Units:** m.\n", "- $v_r,\\ v_t,\\ v_r^*$ — Actual outward radial speed, signed tangential speed, and desired radial descent speed. A negative radial value means descent. **Units:** m/s.\n", "- $a_r^*,\\ a_t^*$ — Requested derivatives dot v_r and dot v_t. These are requests for changing the polar velocity components, not the physical acceleration projections before basis correction. **Units:** m/s².\n", "- $\\theta,\\ x$ — Signed angle from launch and signed surface-arc downrange distance x = R_E $\\theta$; x = 0 is the desired landing site. **Units:** rad; m.\n", "- $\\tau_l$ — Estimated remaining landing duration, recomputed from the current state and bounded below by 3 s. **Units:** s.\n", "- $m,\\ \\mu$ — Current booster mass and Earth gravitational parameter. **Units:** kg; m³/s².\n", "- $D_r,\\ D_t$ — Signed drag-force components in the same radial/tangential frame. Subtracting these in the command compensates the drag already included in the ODE. **Units:** N.\n", "- $\\widetilde F_r,\\ \\widetilde F_t$ — Engine-force components that would realize the requested velocity derivatives before any pointing or magnitude limits. A tilde denotes this unconstrained request. **Units:** N.\n", "- $\\mathbf F_c,\\ Q$ — Requested engine-force vector after setting any downward radial component to zero, and its Euclidean magnitude sqrt($\\mathbf F_c$,r² + $\\mathbf F_c$,t²). Bold symbols are vectors. **Units:** N.\n", "- $T_e,\\ \\mathbf F$ — Single-engine available thrust, and final bounded engine-force vector [F_r, F_t] supplied to the dynamics. **Units:** N.\n", "- $\\dot m,\\ g_0$ — Mass-consumption rate and fixed specific-impulse reference gravity (9.80665 m/s²); the assumed booster specific impulse is 330 s. **Units:** kg/s; m/s².\n", "\n", "$$\n", "v_r^*=-\\sqrt{2(6\\ \\mathrm{m/s^2})h_++(1\\ \\mathrm{m/s})^2}\n", "$$\n", "\n", "This soft-descent profile approaches −1 m/s at ground level.\n", "\n", "$$\n", "a_r^*=-\\frac{(6\\ \\mathrm{m/s^2})v_r}{-v_r^*}+\\frac{v_r^*-v_r}{2\\ \\mathrm s}\n", "$$\n", "\n", "The first term follows the changing descent profile; the second corrects radial velocity error.\n", "\n", "$$\n", "\\tau_l=\\max\\!\\left(3\\ \\mathrm s,\\frac{2h_+}{\\max(1\\ \\mathrm{m/s},-v_r)}\\right)\n", "$$\n", "\n", "The time estimate stays finite near touchdown and during shallow descent.\n", "\n", "$$\n", "a_t^*=-\\frac{6x}{\\tau_l^2}-\\frac{4v_t}{\\tau_l}\n", "$$\n", "\n", "The lateral request corrects position and tangential velocity. Using $v_t$ as $\\dot x$ is a near-surface approximation; the trajectory itself retains the exact polar equations.\n", "\n", "$$\n", "\\widetilde F_r=m\\!\\left(a_r^*+\\frac{\\mu}{r^2}-\\frac{v_t^2}{r}\\right)-D_r\n", "$$\n", "\n", "Convert requested radial acceleration into engine force, compensating gravity, polar motion and drag.\n", "\n", "$$\n", "\\widetilde F_t=m\\!\\left(a_t^*+\\frac{v_rv_t}{r}\\right)-D_t\n", "$$\n", "\n", "Convert the tangential request into engine force with the matching polar and drag compensation.\n", "\n", "$$\n", "\\mathbf F_c=\\begin{bmatrix}\\max(0,\\widetilde F_r)\\\\\\widetilde F_t\\end{bmatrix}\n", "$$\n", "\n", "Rectify the radial component so the engine cannot point downward; $\\mathbf F_c$ is this clipped vector.\n", "\n", "$$\n", "Q=\\|\\mathbf F_c\\|\n", "$$\n", "\n", "$Q$ is the magnitude of the rectified force vector, in newtons.\n", "\n", "$$\n", "\\mathbf F=\\mathbf F_c\\frac{\\operatorname{clip}(Q,0.3T_e,3T_e)}{\\max(Q,10^{-9}\\ \\mathrm N)}\n", "$$\n", "\n", "Scale to the assumed envelope from 30% of one engine to three full engines when $Q\\geq10^{-9}$ N. Below that denominator guard, thrust can fall below the minimum; an exactly zero vector remains zero.\n", "\n", "**Why these equations have this form.**\n", "\n", "$$\n", "\\frac{d v_r^*}{dt}=-\\frac{(6\\ \\mathrm{m/s^2})v_r}{-v_r^*}\\qquad(h>0)\n", "$$\n", "\n", "Differentiate the desired speed with respect to height and use dot h = v_r. This gives the first, feed-forward term in $a_r^*$. Adding ($v_r^*$ − v_r)/(2 s) corrects tracking error. The chosen 6 m/s² shapes the descent-speed profile; it is not an additional force or the full engine acceleration.\n", "\n", "$$\n", "\\begin{aligned}\\dot v_t&=-\\frac{v_rv_t}{r}+\\frac{F_t+D_t}{m}\\\\\\dot v_t&=a_t^*\\\\\\Longrightarrow\\quad F_t&=m\\left(a_t^*+\\frac{v_rv_t}{r}\\right)-D_t\\end{aligned}\n", "$$\n", "\n", "Set the desired derivative dot v_t to $a_t^*$ and rearrange the actual ODE. This is why the command contains a positive v_r v_t/r compensation term: it cancels the negative moving-basis term in the velocity derivative. Tangential physical acceleration is (F_t + D_t)/m, not v_r v_t by itself.\n", "\n", "**Computed reference flight.** The landing burn runs from 368.023 s to ground contact at 475.517 s; 3795.595 kg of propellant remains.\n", "\n", "The update order at each RK45 trial state is: read current height, velocities and mass; construct the desired descent speed and landing-time estimate; request the two velocity derivatives; compensate gravity, basis rotation and drag; rectify the radial force; then limit total force. Feed only this final force into the ODE, with mass flow $-\\|\\mathbf F\\|/(330g_0)$. The ODE produces actual derivatives and RK45 advances the state together. Rectification and magnitude limits can change both achieved derivatives, so the requested $a_r^*$ and $a_t^*$ must never be written straight into the state. Pointing and engine switching are ideal: ignition transients, attitude slew and gimbal limits are absent. The landing-burn time limit is 600 s.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Descent feedback, force compensation and engine envelope** — `python/neutron/model.py:175`\n", "\n", "```python\n", "def landing_guidance(time, state, vehicle, target_downrange_m=0.0, contact_height_m=0.0):\n", " \"\"\"Soft-descent velocity law + finite-time lateral position feedback.\n", "\n", " Ideal engine envelope: 0.3 of one engine to three engines. Engine switching,\n", " ignition transients, attitude slew and gimbal limits are not simulated.\n", " \"\"\"\n", " radius, angle, vr, vt, mass = state\n", " height = max(0.0, radius - EARTH_RADIUS - contact_height_m)\n", " desired_vr = -math.sqrt(2 * 6.0 * height + 1.0**2)\n", " radial_request = -6.0 * vr / -desired_vr + (desired_vr - vr) / 2.0\n", " remaining = max(3.0, 2 * height / max(1.0, -vr))\n", " lateral_request = -6 * (EARTH_RADIUS * angle - target_downrange_m) / remaining**2 - 4 * vt / remaining\n", " acceleration = np.array([radial_request + MU / radius**2 - vt * vt / radius,\n", " lateral_request + vr * vt / radius])\n", " area = math.pi * vehicle.diameter_m**2 / 4\n", " force = mass * acceleration - drag_force(state, area, vehicle.drag_coefficient)\n", " # Only upward pointing during the landing burn; bounded total force.\n", " force[0] = max(0.0, force[0])\n", " magnitude = np.linalg.norm(force)\n", " engine = vehicle.booster_thrust_n / vehicle.booster_engines\n", " return force * np.clip(magnitude, 0.3 * engine, 3 * engine) / max(magnitude, 1e-9)\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "booster-contact", "metadata": {}, "source": [ "### 21 · Stop at contact and score what actually happened\n", "\n", "The downward ground crossing ends the flight. Evaluate impact speed $V_{contact}$ and absolute downrange miss $x_{miss}$ from that integrated state, then turn the engine display off. Do not force position or velocity to an ideal landing.\n", "\n", "**Symbols in this section.**\n", "\n", "- $r,\\ R_E$ — Integrated booster distance from Earth's centre and fixed spherical surface radius; contact occurs when they coincide. **Units:** m.\n", "- $v_r,\\ v_t$ — Actual radial and tangential speeds at contact. Neither is reset when the landing is scored. **Units:** m/s.\n", "- $\\theta$ — Signed angle from the launch-site radius at contact. **Units:** rad.\n", "- $V_{contact},\\ x_{miss}$ — Nonnegative total contact speed and absolute surface-arc distance from launch; both must satisfy their own limit. **Units:** m/s; m.\n", "- $\\downarrow$ — Altitude crossing zero from positive to negative; this is a terminal impact/contact event. **Units:** event direction.\n", "\n", "$$\n", "r-R_E=0\\quad(\\downarrow)\n", "$$\n", "\n", "Ground contact is a solver-located event, not a sampled frame chosen by the renderer.\n", "\n", "$$\n", "V_{contact}=\\sqrt{v_r^2+v_t^2}\n", "$$\n", "\n", "Both radial and tangential velocity contribute to contact speed.\n", "\n", "$$\n", "x_{miss}=|R_E\\theta|\n", "$$\n", "\n", "Measure the actual surface-arc distance from the launch site.\n", "\n", "$$\n", "V_{contact}\\leq2\\ \\mathrm{m/s},\\qquad x_{miss}\\leq100\\ \\mathrm m\n", "$$\n", "\n", "Both tests must pass to label the contact a successful landing.\n", "\n", "**Computed reference flight.** Touchdown is at 475.517 s with speed 1.271738 m/s and miss 0.335251 m.\n", "\n", "Every powered recovery phase has a downward mass gate at 30,000 kg booster dry mass. If fuel runs out, the remaining descent is unpowered for up to 1,000 s. Fast or remote contact is an impact; no contact before the phase limits is a recovery timeout. Contact mechanics and motion after impact are outside the model.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Integrate to ground or actual contact within the finite recovery deck** — `python/neutron/model.py:397`\n", "\n", "```python\n", " impact = terminal_event(lambda t, y: y[0] - EARTH_RADIUS)\n", " deck_events = ()\n", " if body == \"booster\" and ship_recovery:\n", " # A fixed flat deck, not an elevated spherical surface everywhere.\n", " # Check its finite along-track footprint at actual downward roots.\n", " def deck_crossing(t, y):\n", " return (y[0] * math.cos(y[1] - target_angle) - EARTH_RADIUS\n", " - deck_height - leg_clearance * leg_fraction(t, body))\n", " deck_crossing.direction = -1\n", " deck_events = (deck_crossing,)\n", " solution = integrate(state, start, end, guidance, body_area, vehicle.drag_coefficient,\n", " isp, (impact, *monitors, *deck_events), max_step_s, rtol)\n", " integrations += solution.nfev\n", " stop = float(solution.t[-1])\n", " hit = bool(solution.t_events[0].size)\n", " if hit and body == \"booster\":\n", " booster_contact_surface = \"ocean\" if ship_recovery else \"ground\"\n", " if deck_events:\n", " for contact_time in solution.t_events[-1]:\n", " contact_state = solution.sol(contact_time)\n", " along_track = contact_state[0] * math.sin(contact_state[1] - target_angle)\n", " if abs(along_track) <= RECOVERY_SHIP_LENGTH_M / 2:\n", " stop, hit = float(contact_time), True\n", " booster_contact_surface = \"deck\"\n", " solution.t_events = [times[times <= stop + 1e-9] for times in solution.t_events]\n", " break\n", " dense_phases[body].append({\"start\": float(start), \"stop\": stop,\n", " \"solution\": solution, \"guidance\": guidance, \"name\": name})\n", "```\n", "\n", "**Score actual contact, preserve its state and cut displayed thrust** — `python/neutron/model.py:591`\n", "\n", "```python\n", " speed = math.hypot(booster[2], booster[3])\n", " miss = abs(EARTH_RADIUS * booster[1] - recovery_target_downrange_m)\n", " success = (speed <= 2.0 and booster_contact_surface == \"deck\" and leg_fraction(bt, \"booster\") == 1\n", " if ship_recovery else speed <= 2.0 and miss <= 100.0)\n", " contact = trajectories[\"booster\"][-1]\n", " if contact[\"thrust_n\"] > 0:\n", " record(\"booster_engine_cutoff\", bt, \"booster\", \"Engine shuts down at ground contact; contact mechanics are outside this model.\")\n", " trajectories[\"booster\"].append(dict(contact, phase=\"touchdown\" if success else \"impact\",\n", " thrust_n=0.0, force_x_n=0.0, force_y_n=0.0))\n", " landing = dict(success=bool(success), touchdown_time_s=bt, speed_m_s=speed,\n", " miss_distance_m=miss, propellant_remaining_kg=max(0.0, booster[4] - vehicle.booster_dry_kg),\n", " contact_surface=booster_contact_surface)\n", " record(\"touchdown\" if success else \"booster_impact\", bt, \"booster\",\n", " (\"Soft landing on the fixed recovery deck.\" if ship_recovery else \"Soft landing within 100 m of launch.\")\n", " if success else \"Surface contact fails the speed or recovery-target criteria.\")\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "playback", "metadata": {}, "source": [ "## 22 · Turn the state into the display\n", "\n", "The Run button calculates the entire mission in a Python worker before playback begins. The two physical branches are complete. To show them together on the shared mission clock, Python evaluates the dense numerical solution at output times. It records regular two-second samples and exact phase endpoints; the apogee coast and long post-release coasts use ten-second samples. Earth-centred drawing coordinates $X,Y$ come directly from the integrated radius and angle.\n", "\n", "**Symbols in this section.**\n", "\n", "- $r,\\ \\theta$ — Integrated geocentric radius and unwrapped angle. $\\theta$ = 0 lies on the launch-site radius; increasing $\\theta$ is the positive tangential direction. **Units:** m; rad.\n", "- $X,\\ Y$ — Cartesian position components in the fixed simulation plane, with Earth at the origin and launch on the positive X axis. These are physical positions before camera projection. **Units:** m.\n", "- $t,\\ t_a,\\ t_b$ — Requested playback mission time and the neighbouring distinct recorded mission times bracketing it, with t_a ≤ t ≤ t_b. **Units:** s.\n", "- $\\lambda$ — Interpolation fraction between those samples, from 0 at the earlier sample to 1 at the later sample. **Units:** dimensionless.\n", "- $z_a,\\ z_b,\\ z_{display}$ — Earlier recorded value, later recorded value, and displayed interpolated value of one numeric field, such as altitude or mass. z is a generic field, not a sixth physical state. **Units:** same units as the selected field.\n", "\n", "$$\n", "X=r\\cos\\theta,\\qquad Y=r\\sin\\theta\n", "$$\n", "\n", "These coordinates are derived for each recorded physical sample.\n", "\n", "$$\n", "\\lambda=\\frac{t-t_a}{t_b-t_a}\n", "$$\n", "\n", "For a display time $t$ between adjacent recorded times $t_a,t_b$, $\\lambda$ is the interpolation fraction.\n", "\n", "$$\n", "z_{display}=(1-\\lambda)z_a+\\lambda z_b\n", "$$\n", "\n", "For each numeric sample field $z$, the browser interpolates the two neighbouring recorded values. The phase label stays with the earlier sample.\n", "\n", "**Computed reference flight.** The verified payload orbit ends at 8834.470 s; the booster has already touched down at 475.517 s.\n", "\n", "Playback speed scales elapsed wall time into display time; it does not change solver accuracy. Interpolated values never feed back into the physical state or event gates. At or beyond a body's last sample, the browser holds its final state. Complete mission success requires verified insertion, payload release, a verified propagated payload orbit, and successful booster landing together. These checks verify the assumed model's results; they do not validate real vehicle performance.\n", "\n", "\n", "
\n", "Associated code · exact implementation excerpts\n", "\n", "**Dense-output samples and Earth-centred display coordinates** — `python/neutron/model.py:393`\n", "\n", "```python\n", " def phase(body, name, state, start, end, guidance, monitors=(), body_area=None, sample_interval=None):\n", " nonlocal integrations, booster_contact_surface\n", " body_area = area if body_area is None else body_area\n", " isp = vehicle.upper_isp_s if body in (\"upper_stage\", \"payload\") else vehicle.booster_isp_s\n", " impact = terminal_event(lambda t, y: y[0] - EARTH_RADIUS)\n", " deck_events = ()\n", " if body == \"booster\" and ship_recovery:\n", " # A fixed flat deck, not an elevated spherical surface everywhere.\n", " # Check its finite along-track footprint at actual downward roots.\n", " def deck_crossing(t, y):\n", " return (y[0] * math.cos(y[1] - target_angle) - EARTH_RADIUS\n", " - deck_height - leg_clearance * leg_fraction(t, body))\n", " deck_crossing.direction = -1\n", " deck_events = (deck_crossing,)\n", " solution = integrate(state, start, end, guidance, body_area, vehicle.drag_coefficient,\n", " isp, (impact, *monitors, *deck_events), max_step_s, rtol)\n", " integrations += solution.nfev\n", " stop = float(solution.t[-1])\n", " hit = bool(solution.t_events[0].size)\n", " if hit and body == \"booster\":\n", " booster_contact_surface = \"ocean\" if ship_recovery else \"ground\"\n", " if deck_events:\n", " for contact_time in solution.t_events[-1]:\n", " contact_state = solution.sol(contact_time)\n", " along_track = contact_state[0] * math.sin(contact_state[1] - target_angle)\n", " if abs(along_track) <= RECOVERY_SHIP_LENGTH_M / 2:\n", " stop, hit = float(contact_time), True\n", " booster_contact_surface = \"deck\"\n", " solution.t_events = [times[times <= stop + 1e-9] for times in solution.t_events]\n", " break\n", " dense_phases[body].append({\"start\": float(start), \"stop\": stop,\n", " \"solution\": solution, \"guidance\": guidance, \"name\": name})\n", " samples = np.unique(np.r_[start, np.arange(start, stop, sample_interval or sample_step_s), stop])\n", " for time, y in zip(samples, solution.sol(samples).T):\n", " radius, angle, vr, vt, mass = y\n", " force = guidance(time, y)\n", " cos_a, sin_a = math.cos(angle), math.sin(angle)\n", " speed = math.hypot(vr, vt)\n", " thrust = float(np.linalg.norm(force))\n", " force_x = float(force[0] * cos_a - force[1] * sin_a)\n", " force_y = float(force[0] * sin_a + force[1] * cos_a)\n", " if thrust > 0:\n", " pointing[body] = (force_x / thrust, force_y / thrust)\n", " trajectories[body].append(dict(\n", " time_s=float(time), phase=name, altitude_m=float(radius - EARTH_RADIUS),\n", " downrange_m=float(EARTH_RADIUS * angle), radial_velocity_m_s=float(vr),\n", " tangential_velocity_m_s=float(vt), speed_m_s=speed, mass_kg=float(mass),\n", " thrust_n=thrust, force_x_n=force_x, force_y_n=force_y,\n", " axis_x=pointing[body][0], axis_y=pointing[body][1],\n", " fairing_open_fraction=fairing_fraction(time, body),\n", " leg_open_fraction=leg_fraction(time, body),\n", " dynamic_pressure_pa=0.5 * atmosphere(radius - EARTH_RADIUS) * speed**2,\n", " x_m=radius * cos_a, y_m=radius * sin_a,\n", " vx_m_s=vr * cos_a - vt * sin_a, vy_m_s=vr * sin_a + vt * cos_a))\n", " final_state = solution.y[:, -1] if stop == float(solution.t[-1]) else solution.sol(stop)\n", " return stop, final_state.copy(), hit, solution\n", "```\n", "\n", "**Browser interpolation and endpoint handling** — `app/trajectory/neutron-types.ts:87`\n", "\n", "```typescript\n", "export function sampleAt(samples: FlightSample[], time: number) {\n", " if (!samples.length || time < samples[0].time_s) return null;\n", " let low = 0;\n", " let high = samples.length - 1;\n", " while (low < high) {\n", " const middle = Math.ceil((low + high) / 2);\n", " if (samples[middle].time_s <= time) low = middle;\n", " else high = middle - 1;\n", " }\n", " const first = samples[low];\n", " const next = samples[Math.min(low + 1, samples.length - 1)];\n", " const span = next.time_s - first.time_s;\n", " const fraction = span > 0 ? Math.min(1, (time - first.time_s) / span) : 0;\n", " const sample = { ...first };\n", " for (const key of Object.keys(first) as (keyof FlightSample)[]) {\n", " if (key !== 'phase')\n", " sample[key] = first[key] + fraction * (next[key] - first[key]);\n", " }\n", " return { sample, index: low, ended: time > samples.at(-1)!.time_s };\n", "}\n", "\n", "```\n", "\n", "
\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "mission-results", "metadata": {}, "source": [ "## The complete calculated mission\n", "\n", "| Quantity | Integrated result |\n", "| --- | ---: |\n", "| Upper-stage cutoff | 3270.845583 s |\n", "| Insertion perigee | 399.999868 km |\n", "| Insertion apogee | 399.999981 km |\n", "| Payload release | 3280.845583 s |\n", "| Verified unpowered orbit | 5553.624179 s |\n", "| Booster touchdown | 475.517241 s |\n", "| Touchdown speed | 1.271738 m/s |\n", "| Landing miss | 0.335251 m |\n", "| Recovery propellant remaining | 3795.595109 kg |\n", "\n", "Mission success requires verified insertion, successful booster landing, payload release, and the propagated payload-orbit verification together." ] }, { "attachments": {}, "cell_type": "markdown", "id": "mission-source-intro", "metadata": {}, "source": [ "## Appendix · execute the exact model\n", "\n", "The next cell is the complete canonical Python source, unmodified and included exactly once. The original command-line guard does not run in a Jupyter namespace without `__file__`. The final cell runs the reference mission and checks its outcomes.\n", "\n", "Model SHA-256: `ec34d96c89016df203ebc80120e57946d56f96a94820f3d2fc81e48fd28d6e66`." ] }, { "cell_type": "code", "execution_count": 1, "id": "canonical-model", "metadata": { "source_id": "canonical-model" }, "outputs": [], "source": [ "\"\"\"Neutron mission, in readable Python: five states, one ODE, explicit events.\n", "\n", "Educational planar translation model; NOT Rocket Lab flight software or a\n", "performance prediction. Ideal pointing, nonrotating spherical Earth, SI units.\n", "SciPy RK45 is the adaptive Dormand–Prince embedded fifth/fourth-order method.\n", "\"\"\"\n", "\n", "from dataclasses import asdict, dataclass\n", "import argparse\n", "import json\n", "import math\n", "\n", "import numpy as np\n", "from scipy.integrate import solve_ivp\n", "\n", "MU = 3.986004418e14 # Earth gravitational parameter [m^3/s^2]\n", "EARTH_RADIUS = 6_378_137.0 # spherical Earth radius [m]\n", "G0 = 9.80665 # Isp reference gravity [m/s^2]\n", "SOURCE = \"https://rocketlabcorp.com/launch/neutron/\"\n", "TRANSFER_ALTITUDE_KM = 220.0 # chosen mission design, not a physical constant\n", "RECOVERY_DECK_HEIGHT_M = 8.0 # original schematic vessel geometry assumption\n", "RECOVERY_SHIP_LENGTH_M = 122.0\n", "RECOVERY_SHIP_WIDTH_M = 50.0 # assumed beam; planar flight has zero cross-track\n", "RECOVERY_LEG_CLEARANCE_M = 1.6588194256045004 # authored deployed foot below engine plane\n", "\n", "\n", "@dataclass(frozen=True)\n", "class Vehicle:\n", " \"\"\"Published geometry/thrust; explicitly assumed internal mass and Isp.\"\"\"\n", " height_m: float = 43.0\n", " diameter_m: float = 7.0\n", " published_liftoff_mass_kg: float = 480_000.0\n", " advertised_leo_payload_kg: float = 13_000.0\n", " booster_engines: int = 9\n", " booster_thrust_n: float = 1_450_000.0 * 4.4482216152605\n", " upper_thrust_n: float = 890_000.0\n", " booster_dry_kg: float = 30_000.0 # assumption, includes captive fairing\n", " booster_propellant_kg: float = 335_000.0\n", " booster_isp_s: float = 330.0\n", " upper_dry_kg: float = 5_000.0\n", " upper_propellant_kg: float = 102_000.0\n", " upper_isp_s: float = 365.0\n", " drag_coefficient: float = 0.35\n", "\n", "\n", "def atmosphere(altitude_m):\n", " \"\"\"One-layer density model, not a weather or aerothermal model.\"\"\"\n", " return 1.225 * math.exp(-max(0.0, altitude_m) / 8500.0)\n", "\n", "\n", "def orbital_elements(state):\n", " \"\"\"Osculating two-body apsides; a bound ellipse alone is not a safe orbit.\"\"\"\n", " radius, _, vr, vt, _ = map(float, state)\n", " energy = (vr * vr + vt * vt) / 2 - MU / radius\n", " momentum = radius * vt\n", " # Polar eccentricity-vector components avoid cancellation near a circle.\n", " eccentricity = math.hypot(radius * vt * vt / MU - 1, radius * vr * vt / MU)\n", " p = momentum**2 / MU\n", " return {\"perigee_km\": (p / (1 + eccentricity) - EARTH_RADIUS) / 1000,\n", " \"apogee_km\": ((p / (1 - eccentricity) - EARTH_RADIUS) / 1000\n", " if energy < 0 and eccentricity < 1 else None),\n", " \"eccentricity\": eccentricity, \"specific_energy_j_kg\": energy}\n", "\n", "\n", "def orbit_target_error(state, target_altitude_km):\n", " \"\"\"Enter the target set only when BOTH apsides and radial speed agree.\n", "\n", " A proper target orbit is bound, within 2 km at each apsis, |vr| <= 5 m/s,\n", " and e <= 0.001. A negative result satisfies all four constraints.\n", " \"\"\"\n", " orbit = orbital_elements(state)\n", " if orbit[\"apogee_km\"] is None or orbit[\"specific_energy_j_kg\"] >= 0:\n", " return 1e6\n", " return max(abs(orbit[\"perigee_km\"] - target_altitude_km) / 2,\n", " abs(orbit[\"apogee_km\"] - target_altitude_km) / 2,\n", " abs(state[2]) / 5, orbit[\"eccentricity\"] / 0.001) - 1\n", "\n", "\n", "def coast_verification(solution, period_s, target_altitude_km):\n", " \"\"\"Verify an actual propagated revolution; no analytic orbit is substituted.\"\"\"\n", " radius, angle, vr, vt, _ = solution.y\n", " altitude_km = (radius - EARTH_RADIUS) / 1000\n", " energy = (vr**2 + vt**2) / 2 - MU / radius\n", " momentum = radius * vt\n", " energy_change = float(np.max(np.abs((energy - energy[0]) / energy[0])))\n", " momentum_change = float(np.max(np.abs((momentum - momentum[0]) / momentum[0])))\n", " duration = float(solution.t[-1] - solution.t[0])\n", " swept_angle = float(angle[-1] - angle[0])\n", " completed = abs(duration - period_s) < 1e-3 and not solution.t_events[0].size\n", " minimum, maximum = float(altitude_km.min()), float(altitude_km.max())\n", " initial_orbit = orbital_elements(solution.y[:, 0])\n", " success = (completed and np.all(energy < 0) and swept_angle >= 2 * math.pi - 1e-3\n", " and minimum >= target_altitude_km - 2 and maximum <= target_altitude_km + 2\n", " and energy_change < 1e-5 and momentum_change < 1e-5)\n", " return dict(success=bool(success), completed_one_orbit=bool(completed),\n", " period_s=period_s, propagated_duration_s=duration, swept_angle_rad=swept_angle,\n", " perigee_km=initial_orbit[\"perigee_km\"], apogee_km=initial_orbit[\"apogee_km\"],\n", " eccentricity=initial_orbit[\"eccentricity\"],\n", " min_altitude_km=minimum, max_altitude_km=maximum,\n", " relative_energy_change=energy_change, relative_angular_momentum_change=momentum_change)\n", "\n", "\n", "def drag_force(state, area_m2, cd):\n", " \"\"\"Components in the local outward/tangential frame; still atmosphere.\"\"\"\n", " radius, _, vr, vt, _ = state\n", " speed = math.hypot(vr, vt)\n", " scale = -0.5 * atmosphere(radius - EARTH_RADIUS) * cd * area_m2 * speed\n", " return np.array([scale * vr, scale * vt])\n", "\n", "\n", "def rhs(time, state, guidance, area_m2, cd, isp_s):\n", " \"\"\"Polar Newton equations; state = [r, theta, v_r, v_t, mass].\n", "\n", " Pointing is ideal: guidance commands force directly, not attitude/torque.\n", " The changing mass appears in F/m and dm/dt; thrust includes exhaust momentum.\n", " \"\"\"\n", " radius, _, vr, vt, mass = state\n", " force = guidance(time, state)\n", " radial, tangential = (force + drag_force(state, area_m2, cd)) / mass\n", " return [vr, vt / radius,\n", " vt * vt / radius - MU / radius**2 + radial,\n", " -vr * vt / radius + tangential,\n", " -np.linalg.norm(force) / (isp_s * G0)]\n", "\n", "\n", "def terminal_event(function, direction=-1):\n", " \"\"\"SciPy locates a zero crossing between accepted RK45 steps.\"\"\"\n", " function.terminal = True\n", " function.direction = direction\n", " return function\n", "\n", "\n", "def integrate(state, start_s, end_s, guidance, area_m2, cd, isp_s,\n", " events=(), max_step_s=2.0, rtol=1e-7):\n", " \"\"\"Integrate one continuous phase; restart after every discrete change.\"\"\"\n", " solution = solve_ivp(\n", " lambda t, y: rhs(t, y, guidance, area_m2, cd, isp_s),\n", " (start_s, end_s), state, method=\"RK45\", dense_output=True,\n", " events=events, max_step=max_step_s, rtol=rtol,\n", " atol=[1e-3, 1e-11, 1e-5, 1e-5, 1e-4])\n", " if not solution.success:\n", " raise RuntimeError(solution.message)\n", " return solution\n", "\n", "\n", "def ascent_guidance(time, state, vehicle):\n", " \"\"\"An assumed pitch program; angles measured above local horizontal.\"\"\"\n", " pitch = np.interp(time, [0, 12, 35, 70, 115, 160], [90, 90, 82, 67, 52, 40])\n", " beta = math.radians(pitch)\n", " return vehicle.booster_thrust_n * np.array([math.sin(beta), math.cos(beta)])\n", "\n", "\n", "def upper_guidance(time, state, vehicle, target_m):\n", " \"\"\"PD radial acceleration request plus spherical-gravity compensation.\"\"\"\n", " radius, _, vr, vt, mass = state\n", " tau = 30.0\n", " requested = (EARTH_RADIUS + target_m - radius) / tau**2 - 2 * vr / tau\n", " radial_fraction = (requested + MU / radius**2 - vt * vt / radius) * mass / vehicle.upper_thrust_n\n", " beta = math.asin(np.clip(radial_fraction, -0.5, 0.95))\n", " return vehicle.upper_thrust_n * np.array([math.sin(beta), math.cos(beta)])\n", "\n", "\n", "def predicted_miss(state, target_downrange_m=0.0, contact_height_m=0.0):\n", " \"\"\"Cheap constant-g ballistic predictor used only to stop boostback.\n", "\n", " It does not set a landing position; the integrated trajectory can miss.\n", " Landing feedback subsequently corrects remaining position/velocity errors.\n", " \"\"\"\n", " radius, angle, vr, vt, _ = state\n", " gravity = MU / radius**2\n", " fall_time = (vr + math.sqrt(vr * vr + 2 * gravity * max(0, radius - EARTH_RADIUS - contact_height_m))) / gravity\n", " return EARTH_RADIUS * angle + vt * fall_time - target_downrange_m\n", "\n", "\n", "def landing_guidance(time, state, vehicle, target_downrange_m=0.0, contact_height_m=0.0):\n", " \"\"\"Soft-descent velocity law + finite-time lateral position feedback.\n", "\n", " Ideal engine envelope: 0.3 of one engine to three engines. Engine switching,\n", " ignition transients, attitude slew and gimbal limits are not simulated.\n", " \"\"\"\n", " radius, angle, vr, vt, mass = state\n", " height = max(0.0, radius - EARTH_RADIUS - contact_height_m)\n", " desired_vr = -math.sqrt(2 * 6.0 * height + 1.0**2)\n", " radial_request = -6.0 * vr / -desired_vr + (desired_vr - vr) / 2.0\n", " remaining = max(3.0, 2 * height / max(1.0, -vr))\n", " lateral_request = -6 * (EARTH_RADIUS * angle - target_downrange_m) / remaining**2 - 4 * vt / remaining\n", " acceleration = np.array([radial_request + MU / radius**2 - vt * vt / radius,\n", " lateral_request + vr * vt / radius])\n", " area = math.pi * vehicle.diameter_m**2 / 4\n", " force = mass * acceleration - drag_force(state, area, vehicle.drag_coefficient)\n", " # Only upward pointing during the landing burn; bounded total force.\n", " force[0] = max(0.0, force[0])\n", " magnitude = np.linalg.norm(force)\n", " engine = vehicle.booster_thrust_n / vehicle.booster_engines\n", " return force * np.clip(magnitude, 0.3 * engine, 3 * engine) / max(magnitude, 1e-9)\n", "\n", "\n", "def separate(state, booster_mass_kg, relative_speed_m_s=0.5):\n", " \"\"\"Split mass with equal/opposite impulses: conserve linear momentum.\n", "\n", " Both bodies start at the same position. Geometry and contact are not modeled.\n", " The reusable fairing stays within booster dry mass.\n", " \"\"\"\n", " booster, upper = state.copy(), state.copy()\n", " upper[4] = state[4] - booster_mass_kg\n", " booster[4] = booster_mass_kg\n", " booster[2] -= relative_speed_m_s * upper[4] / state[4]\n", " upper[2] += relative_speed_m_s * booster[4] / state[4]\n", " return booster, upper\n", "\n", "\n", "COAST_POINTING_RATE_RAD_S = math.radians(60.0)\n", "UPPER_EXIT_CLEARANCE_M = 43.0 - 24.8 + 5.0 # authored enclosure, upper base and margin\n", "\n", "\n", "def apply_kinematic_hardware(trajectories, events, dense_phases):\n", " \"\"\"Prescribed hardware coordinates, sampled from the actual integrated flight.\n", "\n", " Coast pointing uses a bounded cubic pre-alignment before the next ignition.\n", " Powered pointing is the actual force axis. Fairing closure waits for physical\n", " upper-stage clearance. No position, velocity, mass or applied force changes.\n", " \"\"\"\n", " from scipy.optimize import brentq\n", "\n", " def phase_at(body, time):\n", " return next((phase for phase in reversed(dense_phases[body])\n", " if phase[\"start\"] <= time <= phase[\"stop\"]), None)\n", "\n", " def state_xy(body, time):\n", " state = phase_at(body, time)[\"solution\"].sol(time)\n", " return state[0] * np.array([math.cos(state[1]), math.sin(state[1])])\n", "\n", " def insert_samples(body, times):\n", " rows = trajectories[body]\n", " for time in times:\n", " if any(abs(row[\"time_s\"] - time) < 1e-8 for row in rows):\n", " continue\n", " phase = phase_at(body, time)\n", " if phase is None:\n", " continue\n", " state = phase[\"solution\"].sol(time)\n", " radius, angle, vr, vt, mass = map(float, state)\n", " force = phase[\"guidance\"](time, state)\n", " thrust = float(np.linalg.norm(force))\n", " ca, sa = math.cos(angle), math.sin(angle)\n", " fx, fy = float(force[0] * ca - force[1] * sa), float(force[0] * sa + force[1] * ca)\n", " before = max((row for row in rows if row[\"time_s\"] <= time), key=lambda row: row[\"time_s\"])\n", " row = dict(before, time_s=float(time), phase=phase[\"name\"],\n", " altitude_m=radius - EARTH_RADIUS, downrange_m=EARTH_RADIUS * angle,\n", " radial_velocity_m_s=vr, tangential_velocity_m_s=vt,\n", " speed_m_s=math.hypot(vr, vt), mass_kg=mass, thrust_n=thrust,\n", " force_x_n=fx, force_y_n=fy,\n", " dynamic_pressure_pa=0.5 * atmosphere(radius - EARTH_RADIUS) * (vr * vr + vt * vt),\n", " x_m=radius * ca, y_m=radius * sa,\n", " vx_m_s=vr * ca - vt * sa, vy_m_s=vr * sa + vt * ca)\n", " if thrust > 0:\n", " row[\"axis_x\"], row[\"axis_y\"] = fx / thrust, fy / thrust\n", " if body == \"booster\" and any(sample.get(\"leg_open_fraction\", 0) > 0 for sample in rows):\n", " command = next((event[\"time_s\"] for event in events if event[\"name\"] == \"landing_burn\"), math.inf)\n", " row[\"leg_open_fraction\"] = float(np.clip((time - command) / 2, 0, 1))\n", " rows.append(row)\n", " rows.sort(key=lambda row: row[\"time_s\"])\n", "\n", " for body, rows in trajectories.items():\n", " segments = []\n", " index = 0\n", " while index < len(rows):\n", " if rows[index][\"thrust_n\"] > 0:\n", " index += 1\n", " continue\n", " start = index\n", " while index < len(rows) and rows[index][\"thrust_n\"] == 0:\n", " index += 1\n", " if index == len(rows):\n", " break\n", " first, target = rows[start], rows[index]\n", " initial = math.atan2(first[\"axis_y\"], first[\"axis_x\"])\n", " final = math.atan2(target[\"force_y_n\"], target[\"force_x_n\"])\n", " turn = math.atan2(math.sin(final - initial), math.cos(final - initial))\n", " available = target[\"time_s\"] - first[\"time_s\"]\n", " duration = min(available, max(2.0, 1.5 * abs(turn) / COAST_POINTING_RATE_RAD_S))\n", " if duration > 0:\n", " segments.append((target[\"time_s\"] - duration, target[\"time_s\"], initial, turn))\n", " # Include the profile endpoints and 0.1 s kinematic samples even when\n", " # the caller chooses sparse trajectory output. States use RK45 dense\n", " # output here; they are never linearly invented or re-integrated.\n", " insert_samples(body, [float(t) for begin, end, _, _ in segments\n", " for t in np.r_[np.arange(begin, end, 0.1), end]])\n", " for begin, end, initial, turn in segments:\n", " for row in rows:\n", " if row[\"thrust_n\"] == 0 and begin <= row[\"time_s\"] <= end:\n", " u = float(np.clip((row[\"time_s\"] - begin) / (end - begin), 0, 1))\n", " angle = initial + turn * u * u * (3 - 2 * u)\n", " row[\"axis_x\"], row[\"axis_y\"] = math.cos(angle), math.sin(angle)\n", "\n", " opening = next((event[\"time_s\"] for event in events if event[\"name\"] == \"fairing_opening\"), None)\n", " separation = next((event[\"time_s\"] for event in events if event[\"name\"] == \"stage_separation\"), None)\n", " if opening is None or separation is None or not trajectories[\"booster\"] or not trajectories[\"upper_stage\"]:\n", " return\n", " attached = trajectories[\"stack\"][-1]\n", " axis = np.array([attached[\"axis_x\"], attached[\"axis_y\"]])\n", " def clearance(time):\n", " return float((state_xy(\"upper_stage\", time) - state_xy(\"booster\", time)) @ axis) - UPPER_EXIT_CLEARANCE_M\n", " stop = min(dense_phases[body][-1][\"stop\"] for body in (\"booster\", \"upper_stage\"))\n", " times = sorted(set(float(t) for body in (\"booster\", \"upper_stage\")\n", " for phase in dense_phases[body] for t in phase[\"solution\"].t\n", " if separation <= t <= stop))\n", " close = math.inf\n", " for left, right in zip(times, times[1:]):\n", " if clearance(left) <= 0 <= clearance(right):\n", " close = float(brentq(clearance, left, right, xtol=1e-9))\n", " break\n", " events[:] = [event for event in events if event[\"name\"] not in (\"fairing_closing\", \"fairing_closed\")]\n", " if math.isfinite(close):\n", " for name, time, text in [\n", " (\"fairing_closing\", close, \"Captive fairing closes after integrated upper-stage clearance.\"),\n", " (\"fairing_closed\", close + 4, \"Four-second closure completes after upper-stage clearance.\")]:\n", " events.append(dict(name=name, time_s=time, body=\"booster\", description=text,\n", " required_axial_clearance_m=UPPER_EXIT_CLEARANCE_M))\n", " insert_samples(\"booster\", [close, close + 4])\n", " for body in (\"stack\", \"booster\"):\n", " for row in trajectories[body]:\n", " time = row[\"time_s\"]\n", " fraction = min(1.0, max(0.0, (time - opening) / 4))\n", " if time > close:\n", " fraction = max(0.0, 1 - (time - close) / 4)\n", " row[\"fairing_open_fraction\"] = fraction\n", "\n", "\n", "def simulate(payload_kg=8000.0, recovery_propellant_kg=60_000.0,\n", " target_altitude_km=400.0, upper_propellant_scale=1.0,\n", " max_step_s=2.0, rtol=1e-7, sample_step_s=2.0,\n", " recovery_target_downrange_m=0.0):\n", " \"\"\"Run ascent, deployment and independent booster recovery; return JSON data.\n", "\n", " Inputs change actual integrated states. Insufficient fuel, missed orbit and\n", " hard/remote touchdowns remain failures. No endpoint is moved to the target.\n", " \"\"\"\n", " vehicle = Vehicle()\n", " values = [payload_kg, recovery_propellant_kg, target_altitude_km,\n", " upper_propellant_scale, max_step_s, rtol, sample_step_s,\n", " recovery_target_downrange_m]\n", " if not all(math.isfinite(value) for value in values):\n", " raise ValueError(\"All inputs must be finite.\")\n", " if payload_kg <= 0 or not 0 <= recovery_propellant_kg < vehicle.booster_propellant_kg:\n", " raise ValueError(\"Payload must be positive; reserve must be inside stage-one fuel capacity.\")\n", " if not 150 <= target_altitude_km <= 1000 or not 0 < upper_propellant_scale <= 1.5:\n", " raise ValueError(\"Target must be 150–1000 km; upper-stage fuel scale must be (0, 1.5].\")\n", " if min(max_step_s, rtol, sample_step_s) <= 0:\n", " raise ValueError(\"Solver and sampling controls must be positive.\")\n", " if not 0 <= recovery_target_downrange_m <= 2_000_000:\n", " raise ValueError(\"Recovery target must be between 0 and 2,000,000 m downrange.\")\n", "\n", " ship_recovery = recovery_target_downrange_m > 0\n", " target_angle = recovery_target_downrange_m / EARTH_RADIUS\n", " deck_height = RECOVERY_DECK_HEIGHT_M if ship_recovery else 0.0\n", " leg_clearance = RECOVERY_LEG_CLEARANCE_M if ship_recovery else 0.0\n", " contact_height = deck_height + leg_clearance\n", " booster_contact_surface = None\n", "\n", " area = math.pi * vehicle.diameter_m**2 / 4\n", " upper_wet = vehicle.upper_dry_kg + vehicle.upper_propellant_kg * upper_propellant_scale + payload_kg\n", " initial_mass = vehicle.booster_dry_kg + vehicle.booster_propellant_kg + upper_wet\n", " initial = np.array([EARTH_RADIUS, 0.0, 0.0, 0.0, initial_mass])\n", " trajectories = {body: [] for body in (\"stack\", \"booster\", \"upper_stage\", \"payload\")}\n", " dense_phases = {body: [] for body in trajectories}\n", " events = []\n", " integrations = 0\n", " separation_time = None\n", " # Display pointing is the model's ideal thrust axis, not an attitude plant.\n", " # Separated bodies inherit their parent's axis. The explicit kinematic\n", " # hardware pass below pre-aligns coast segments before subsequent ignitions.\n", " pointing = {body: (1.0, 0.0) for body in trajectories}\n", "\n", " def record(name, time, body, description, **details):\n", " events.append(dict(name=name, time_s=float(time), body=body,\n", " description=description, **details))\n", "\n", " def fairing_fraction(time, body):\n", " opening = next((e[\"time_s\"] for e in events if e[\"name\"] == \"fairing_opening\"), None)\n", " if opening is None or body not in (\"stack\", \"booster\"):\n", " return 0.0\n", " if separation_time is None or time <= separation_time:\n", " return float(np.clip((time - opening) / 4, 0, 1))\n", " return float(np.clip(1 - (time - separation_time) / 4, 0, 1))\n", "\n", " def leg_fraction(time, body):\n", " command = next((e[\"time_s\"] for e in events if e[\"name\"] == \"landing_burn\"), None)\n", " if not ship_recovery or body != \"booster\" or command is None:\n", " return 0.0\n", " return float(np.clip((time - command) / 2, 0, 1))\n", "\n", " def phase(body, name, state, start, end, guidance, monitors=(), body_area=None, sample_interval=None):\n", " nonlocal integrations, booster_contact_surface\n", " body_area = area if body_area is None else body_area\n", " isp = vehicle.upper_isp_s if body in (\"upper_stage\", \"payload\") else vehicle.booster_isp_s\n", " impact = terminal_event(lambda t, y: y[0] - EARTH_RADIUS)\n", " deck_events = ()\n", " if body == \"booster\" and ship_recovery:\n", " # A fixed flat deck, not an elevated spherical surface everywhere.\n", " # Check its finite along-track footprint at actual downward roots.\n", " def deck_crossing(t, y):\n", " return (y[0] * math.cos(y[1] - target_angle) - EARTH_RADIUS\n", " - deck_height - leg_clearance * leg_fraction(t, body))\n", " deck_crossing.direction = -1\n", " deck_events = (deck_crossing,)\n", " solution = integrate(state, start, end, guidance, body_area, vehicle.drag_coefficient,\n", " isp, (impact, *monitors, *deck_events), max_step_s, rtol)\n", " integrations += solution.nfev\n", " stop = float(solution.t[-1])\n", " hit = bool(solution.t_events[0].size)\n", " if hit and body == \"booster\":\n", " booster_contact_surface = \"ocean\" if ship_recovery else \"ground\"\n", " if deck_events:\n", " for contact_time in solution.t_events[-1]:\n", " contact_state = solution.sol(contact_time)\n", " along_track = contact_state[0] * math.sin(contact_state[1] - target_angle)\n", " if abs(along_track) <= RECOVERY_SHIP_LENGTH_M / 2:\n", " stop, hit = float(contact_time), True\n", " booster_contact_surface = \"deck\"\n", " solution.t_events = [times[times <= stop + 1e-9] for times in solution.t_events]\n", " break\n", " dense_phases[body].append({\"start\": float(start), \"stop\": stop,\n", " \"solution\": solution, \"guidance\": guidance, \"name\": name})\n", " samples = np.unique(np.r_[start, np.arange(start, stop, sample_interval or sample_step_s), stop])\n", " for time, y in zip(samples, solution.sol(samples).T):\n", " radius, angle, vr, vt, mass = y\n", " force = guidance(time, y)\n", " cos_a, sin_a = math.cos(angle), math.sin(angle)\n", " speed = math.hypot(vr, vt)\n", " thrust = float(np.linalg.norm(force))\n", " force_x = float(force[0] * cos_a - force[1] * sin_a)\n", " force_y = float(force[0] * sin_a + force[1] * cos_a)\n", " if thrust > 0:\n", " pointing[body] = (force_x / thrust, force_y / thrust)\n", " trajectories[body].append(dict(\n", " time_s=float(time), phase=name, altitude_m=float(radius - EARTH_RADIUS),\n", " downrange_m=float(EARTH_RADIUS * angle), radial_velocity_m_s=float(vr),\n", " tangential_velocity_m_s=float(vt), speed_m_s=speed, mass_kg=float(mass),\n", " thrust_n=thrust, force_x_n=force_x, force_y_n=force_y,\n", " axis_x=pointing[body][0], axis_y=pointing[body][1],\n", " fairing_open_fraction=fairing_fraction(time, body),\n", " leg_open_fraction=leg_fraction(time, body),\n", " dynamic_pressure_pa=0.5 * atmosphere(radius - EARTH_RADIUS) * speed**2,\n", " x_m=radius * cos_a, y_m=radius * sin_a,\n", " vx_m_s=vr * cos_a - vt * sin_a, vy_m_s=vr * sin_a + vt * cos_a))\n", " final_state = solution.y[:, -1] if stop == float(solution.t[-1]) else solution.sol(stop)\n", " return stop, final_state.copy(), hit, solution\n", "\n", " coast = lambda t, y: np.zeros(2)\n", " reserve_mass = vehicle.booster_dry_kg + recovery_propellant_kg\n", " meco = terminal_event(lambda t, y: y[4] - upper_wet - reserve_mass)\n", " record(\"liftoff\", 0, \"stack\", \"Nine-engine ascent begins.\")\n", " time, state, impacted, ascent_solution = phase(\"stack\", \"ascent\", initial, 0, 600,\n", " lambda t, y: ascent_guidance(t, y, vehicle), (meco,))\n", " ascent_complete = bool(ascent_solution.t_events[1].size) and not impacted\n", " orbit = {\"success\": False, \"perigee_km\": None, \"apogee_km\": None,\n", " \"eccentricity\": None, \"cutoff_time_s\": None}\n", " landing = {\"success\": False, \"touchdown_time_s\": None, \"speed_m_s\": None,\n", " \"miss_distance_m\": None, \"propellant_remaining_kg\": None}\n", " released = False\n", " payload_orbit = None\n", "\n", " if ascent_complete:\n", " record(\"meco\", time, \"stack\", \"Main engine cutoff preserves the requested recovery reserve.\",\n", " recovery_propellant_kg=float(recovery_propellant_kg))\n", " record(\"fairing_opening\", time, \"stack\", \"Captive fairing opens over four seconds (kinematic mechanism).\")\n", " time, state, impacted, _ = phase(\"stack\", \"fairing_opening\", state, time, time + 4, coast)\n", " if ascent_complete and not impacted:\n", " separation_time = time\n", " booster, upper = separate(state, reserve_mass)\n", " pointing[\"booster\"] = pointing[\"upper_stage\"] = pointing[\"stack\"]\n", " record(\"fairing_open\", time, \"stack\", \"Captive fairing is fully open.\")\n", " record(\"stage_separation\", time, \"stack\", \"Upper stage released with 0.5 m/s relative radial speed.\",\n", " mass_before_kg=float(state[4]), booster_mass_kg=float(booster[4]),\n", " upper_mass_kg=float(upper[4]))\n", " record(\"fairing_closing\", time, \"booster\", \"Fairing remains attached and closes over four seconds.\")\n", "\n", " # UPPER BRANCH: coast clear, ignite, stop at orbit or fuel floor, deploy.\n", " upper_time, upper, upper_hit, _ = phase(\"upper_stage\", \"separation_coast\", upper, time, time + 3, coast)\n", " fuel_floor = vehicle.upper_dry_kg + payload_kg\n", " empty = terminal_event(lambda t, y: y[4] - fuel_floor)\n", " # This monotonic speed crossing cannot skip the narrow near-circular set\n", " # between solver steps; all orbit constraints are checked at its root.\n", " target = terminal_event(lambda t, y: y[3] - math.sqrt(MU / y[0]), 1)\n", " reached = False\n", " if not upper_hit:\n", " record(\"upper_ignition\", upper_time, \"upper_stage\", \"Single vacuum Archimedes engine ignites.\")\n", " transfer = target_altitude_km > TRANSFER_ALTITUDE_KM\n", " def transfer_apogee(t, y):\n", " apogee = orbital_elements(y)[\"apogee_km\"]\n", " return (apogee if apogee is not None else -1e6) - target_altitude_km\n", " transfer_gate = terminal_event(transfer_apogee, 1)\n", " upper_time, upper, upper_hit, solution = phase(\n", " \"upper_stage\", \"transfer_insertion\" if transfer else \"orbit_insertion\", upper, upper_time, upper_time + 1200,\n", " lambda t, y: upper_guidance(t, y, vehicle, min(target_altitude_km, TRANSFER_ALTITUDE_KM) * 1000),\n", " (empty, transfer_gate if transfer else target), 12.0)\n", " cutoff_triggered = bool(solution.t_events[2].size)\n", " if transfer and cutoff_triggered and not upper_hit:\n", " record(\"transfer_cutoff\", upper_time, \"upper_stage\", \"Engine stops when the actual transfer-orbit apogee reaches the requested altitude.\")\n", " apex = terminal_event(lambda t, y: y[2])\n", " upper_time, upper, upper_hit, solution = phase(\"upper_stage\", \"apogee_coast\", upper,\n", " upper_time, upper_time + 6000, coast, (apex,), 12.0, max(sample_step_s, 10.0))\n", " cutoff_triggered = False\n", " if solution.t_events[1].size and not upper_hit:\n", " record(\"apogee\", upper_time, \"upper_stage\", \"Radial velocity crosses zero at the integrated apogee.\")\n", " record(\"circularization_ignition\", upper_time, \"upper_stage\", \"Assumed vacuum-engine restart begins a finite circularization burn.\")\n", " upper_time, upper, upper_hit, solution = phase(\"upper_stage\", \"circularization_burn\", upper,\n", " upper_time, upper_time + 600,\n", " lambda t, y: upper_guidance(t, y, vehicle, target_altitude_km * 1000), (empty, target), 12.0)\n", " cutoff_triggered = bool(solution.t_events[2].size)\n", " elements = orbital_elements(upper)\n", " reached = bool(cutoff_triggered and not upper_hit\n", " and orbit_target_error(upper, target_altitude_km) <= 1e-8)\n", " orbit = dict(elements, success=reached, cutoff_time_s=upper_time,\n", " target_altitude_km=float(target_altitude_km), radial_velocity_m_s=float(upper[2]),\n", " tangential_velocity_m_s=float(upper[3]), circular_speed_m_s=math.sqrt(MU / upper[0]))\n", " if not upper_hit:\n", " 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\"\n", " 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.\")\n", " if upper_hit:\n", " record(\"upper_impact\", upper_time, \"upper_stage\", \"Upper stage contacted the ground.\")\n", " if not upper_hit:\n", " coast_name, duration = (\"deployment_coast\", 10) if reached else (\"failed_insertion_coast\", 1200)\n", " upper_time, upper, upper_hit, _ = phase(\"upper_stage\", coast_name, upper,\n", " upper_time, upper_time + duration, coast, body_area=12.0)\n", " if upper_hit:\n", " record(\"upper_impact\", upper_time, \"upper_stage\", \"Upper stage contacted the ground during unpowered flight.\")\n", " if reached and not upper_hit:\n", " retained, payload = separate(upper, upper[4] - payload_kg, 0.2)\n", " pointing[\"payload\"] = pointing[\"upper_stage\"]\n", " record(\"payload_release\", upper_time, \"payload\", \"Payload released at 0.2 m/s relative radial speed; both bodies propagate independently.\",\n", " payload_mass_kg=float(payload[4]), retained_mass_kg=float(retained[4]), mass_before_kg=float(upper[4]))\n", " released = True\n", " energy = orbital_elements(payload)[\"specific_energy_j_kg\"]\n", " semimajor_axis = -MU / (2 * energy)\n", " period = 2 * math.pi * math.sqrt(semimajor_axis**3 / MU)\n", " # Keep long-coast exports compact without changing solver steps.\n", " interval = max(sample_step_s, 10.0)\n", " phase(\"upper_stage\", \"post_release_coast\", retained, upper_time, upper_time + period,\n", " coast, body_area=12.0, sample_interval=interval)\n", " payload_time, _, _, solution = phase(\"payload\", \"free_flight\", payload,\n", " upper_time, upper_time + period, coast, body_area=2.0, sample_interval=interval)\n", " payload_orbit = coast_verification(solution, period, target_altitude_km)\n", " record(\"payload_orbit_verified\" if payload_orbit[\"success\"] else \"payload_orbit_failed\",\n", " payload_time, \"payload\", \"One full unpowered revolution checked against altitude and conservation limits.\")\n", "\n", " # BOOSTER BRANCH: close fairing, reverse downrange motion, coast, land.\n", " bt, booster, booster_hit, _ = phase(\"booster\", \"clearance_coast\", booster, time, time + 4, coast)\n", " engine = vehicle.booster_thrust_n / vehicle.booster_engines\n", " empty = terminal_event(lambda t, y: y[4] - vehicle.booster_dry_kg)\n", " initial_miss = predicted_miss(booster, recovery_target_downrange_m, contact_height)\n", " correction = -1.0 if initial_miss >= 0 else 1.0\n", " return_target = terminal_event(\n", " lambda t, y: predicted_miss(y, recovery_target_downrange_m, contact_height), correction)\n", " if not booster_hit:\n", " record(\"fairing_closed\", bt, \"booster\", \"Captive fairing closed for recovery.\")\n", " if not booster_hit and booster[4] > vehicle.booster_dry_kg + 1e-3:\n", " record(\"boostback_ignition\", bt, \"booster\",\n", " \"Three-engine range correction targets the fixed recovery ship.\" if ship_recovery\n", " else \"Three-engine boostback targets the launch site (assumed RTLS scenario).\")\n", " bt, booster, booster_hit, solution = phase(\"booster\", \"boostback\", booster, bt, bt + 300,\n", " lambda t, y: np.array([0.0, correction * 3 * engine]), (empty, return_target))\n", " if not booster_hit:\n", " reason = \"predicted zero miss\" if solution.t_events[2].size else \"propellant depletion\" if solution.t_events[1].size else \"time limit\"\n", " record(\"boostback_cutoff\", bt, \"booster\", \"Boostback stops at \" + reason + \".\")\n", " has_fuel = booster[4] > vehicle.booster_dry_kg + 1e-3\n", "\n", " def landing_gate(t, y):\n", " height, vr, mass = y[0] - EARTH_RADIUS - contact_height, y[2], y[4]\n", " maximum_deceleration = max(1.0, 3 * engine / mass - MU / y[0]**2)\n", " stopping_distance = min(vr, 0)**2 / (2 * maximum_deceleration)\n", " return height - 1.5 * stopping_distance - 3000\n", "\n", " if not booster_hit and has_fuel:\n", " gate = terminal_event(landing_gate)\n", " reentry = lambda t, y: y[0] - EARTH_RADIUS - 80_000\n", " reentry.direction = -1\n", " bt, booster, booster_hit, solution = phase(\"booster\", \"coast_and_reentry\", booster, bt, bt + 1000,\n", " coast, (gate, reentry))\n", " if solution.t_events[2].size:\n", " record(\"atmospheric_reentry\", solution.t_events[2][0], \"booster\", \"Descending through 80 km; exponential-atmosphere drag acts continuously.\")\n", " if not booster_hit and solution.t_events[1].size:\n", " record(\"landing_burn\", bt, \"booster\", \"Feedback landing burn begins from the actual descending state.\")\n", " bt, booster, booster_hit, _ = phase(\"booster\", \"landing_burn\", booster, bt, bt + 600,\n", " lambda t, y: landing_guidance(t, y, vehicle, recovery_target_downrange_m, contact_height), (empty,))\n", " if not booster_hit and booster[4] <= vehicle.booster_dry_kg + 1e-3:\n", " record(\"booster_fuel_depleted\", bt, \"booster\", \"Recovery propellant exhausted; remaining flight is ballistic.\")\n", " bt, booster, booster_hit, _ = phase(\"booster\", \"unpowered_descent\", booster, bt, bt + 1000, coast)\n", " if booster_hit:\n", " speed = math.hypot(booster[2], booster[3])\n", " miss = abs(EARTH_RADIUS * booster[1] - recovery_target_downrange_m)\n", " success = (speed <= 2.0 and booster_contact_surface == \"deck\" and leg_fraction(bt, \"booster\") == 1\n", " if ship_recovery else speed <= 2.0 and miss <= 100.0)\n", " contact = trajectories[\"booster\"][-1]\n", " if contact[\"thrust_n\"] > 0:\n", " record(\"booster_engine_cutoff\", bt, \"booster\", \"Engine shuts down at ground contact; contact mechanics are outside this model.\")\n", " trajectories[\"booster\"].append(dict(contact, phase=\"touchdown\" if success else \"impact\",\n", " thrust_n=0.0, force_x_n=0.0, force_y_n=0.0))\n", " landing = dict(success=bool(success), touchdown_time_s=bt, speed_m_s=speed,\n", " miss_distance_m=miss, propellant_remaining_kg=max(0.0, booster[4] - vehicle.booster_dry_kg),\n", " contact_surface=booster_contact_surface)\n", " record(\"touchdown\" if success else \"booster_impact\", bt, \"booster\",\n", " (\"Soft landing on the fixed recovery deck.\" if ship_recovery else \"Soft landing within 100 m of launch.\")\n", " if success else \"Surface contact fails the speed or recovery-target criteria.\")\n", " else:\n", " record(\"recovery_timeout\", bt, \"booster\", \"No ground contact within the recovery time limit.\")\n", " else:\n", " record(\"launch_failure\" if impacted else \"ascent_timeout\", time, \"stack\",\n", " \"Vehicle contacted the ground before separation.\" if impacted else \"Recovery-reserve cutoff was not reached.\")\n", "\n", " apply_kinematic_hardware(trajectories, events, dense_phases)\n", " # Sampling is presentation only; accepted integrator steps locate all events.\n", " events.sort(key=lambda event: event[\"time_s\"])\n", " samples = [sample for path in trajectories.values() for sample in path]\n", " parameters = {\"vehicle\": asdict(vehicle), \"payload_kg\": float(payload_kg),\n", " \"recovery_propellant_kg\": float(recovery_propellant_kg),\n", " \"recovery_target_downrange_m\": float(recovery_target_downrange_m),\n", " \"recovery_target\": {\"kind\": \"ship\" if ship_recovery else \"launch\",\n", " \"downrange_m\": float(recovery_target_downrange_m),\n", " \"deck_height_m\": deck_height, \"leg_clearance_m\": leg_clearance,\n", " \"length_m\": RECOVERY_SHIP_LENGTH_M if ship_recovery else 0.0,\n", " \"width_m\": RECOVERY_SHIP_WIDTH_M if ship_recovery else 0.0,\n", " \"x_m\": EARTH_RADIUS * math.cos(target_angle),\n", " \"y_m\": EARTH_RADIUS * math.sin(target_angle)},\n", " \"target_altitude_km\": float(target_altitude_km),\n", " \"upper_propellant_scale\": float(upper_propellant_scale),\n", " \"source_url\": SOURCE, \"source_checked\": \"2026-09-26\",\n", " \"published_fields\": [\"height_m\", \"diameter_m\", \"published_liftoff_mass_kg\",\n", " \"advertised_leo_payload_kg\", \"booster_engines\", \"booster_thrust_n\", \"upper_thrust_n\"],\n", " \"assumed_fields\": [\"booster_dry_kg\", \"booster_propellant_kg\", \"booster_isp_s\",\n", " \"upper_dry_kg\", \"upper_propellant_kg\", \"upper_isp_s\", \"drag_coefficient\"],\n", " \"initial_mass_kg\": initial_mass, \"solver\": \"SciPy RK45 / Dormand–Prince 5(4)\",\n", " \"rtol\": rtol, \"max_step_s\": max_step_s,\n", " \"landing_limits\": {\"speed_m_s\": 2.0, \"miss_distance_m\": RECOVERY_SHIP_LENGTH_M / 2 if ship_recovery else 100.0},\n", " \"orbit_limits\": {\"apsis_error_km\": 2.0, \"radial_speed_m_s\": 5.0, \"eccentricity\": 0.001},\n", " \"assumptions\": [\"Nonrotating spherical Earth and still exponential atmosphere.\",\n", " \"Planar point mass with ideal instantaneous thrust direction; no attitude dynamics.\",\n", " \"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.\",\n", " \"Assumed mass split, Isp, drag, guidance, reserve, separation impulses and actuator timing.\",\n", " \"Above 220 km: first burn raises transfer apogee, coast reaches apogee, and an assumed vacuum-engine restart circularizes with a finite burn.\",\n", " \"Landing thrust envelope: 30% of one to 100% of three engines, with ideal switching.\",\n", " \"Published 480000 kg is a reference gross mass; actual initial mass varies with payload/fuel.\",\n", " \"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.\"],\n", " \"recovery_scenario\": (\"Fixed downrange sea-recovery target; schematic deck and axial foot contact, no waves or rigid-body contact dynamics.\"\n", " if ship_recovery else \"Educational return to launch site; published baseline uses downrange sea landing.\")}\n", " return {\"model\": \"Neutron planar mission v1\", \"parameters\": parameters,\n", " \"trajectories\": trajectories, \"events\": events,\n", " \"summary\": {\"orbit\": orbit, \"landing\": landing, \"payload_released\": released,\n", " \"payload_orbit\": payload_orbit,\n", " \"mission_success\": bool(orbit[\"success\"] and landing[\"success\"] and released\n", " and payload_orbit and payload_orbit[\"success\"]),\n", " \"duration_s\": max(s[\"time_s\"] for s in samples),\n", " \"max_dynamic_pressure_pa\": max(s[\"dynamic_pressure_pa\"] for s in samples),\n", " \"function_evaluations\": integrations}}\n", "\n", "\n", "if __name__ == \"__main__\" and \"__file__\" in globals():\n", " parser = argparse.ArgumentParser(description=__doc__)\n", " parser.add_argument(\"--json\", action=\"store_true\", help=\"Print complete sampled mission JSON\")\n", " parser.add_argument(\"--payload\", type=float, default=8000.0, help=\"Payload mass in kg\")\n", " args = parser.parse_args()\n", " result = simulate(payload_kg=args.payload)\n", " print(json.dumps(result if args.json else result[\"summary\"], indent=2, allow_nan=False))\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "run-mission", "metadata": { "source_id": "run-mission" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The two branches share a separation epoch and continue independently.\n", "Separation: 141.978968 s\n", "Upper cutoff: 3270.845583 s\n", "Insertion apsides: 399.999868 × 399.999981 km\n", "Payload release: 3280.845583 s\n", "Verified orbit: 5553.624179 s\n", "Payload altitude range: 399.929456–400.070394 km\n", "Booster touchdown: 475.517241 s\n", "Contact speed: 1.271738 m/s; miss: 0.335251 m\n", "Recovery propellant remaining: 3795.595109 kg\n", "Integrated insertion, payload orbit and booster landing checks passed.\n" ] } ], "source": [ "flight = simulate()\n", "summary = flight[\"summary\"]\n", "orbit, landing = summary[\"orbit\"], summary[\"landing\"]\n", "verified = summary[\"payload_orbit\"]\n", "event_times = {event[\"name\"]: event[\"time_s\"] for event in flight[\"events\"]}\n", "assert summary[\"mission_success\"]\n", "assert orbit[\"success\"] and landing[\"success\"] and summary[\"payload_released\"]\n", "assert verified[\"success\"] and verified[\"completed_one_orbit\"]\n", "assert event_times[\"touchdown\"] < event_times[\"transfer_cutoff\"]\n", "print(\"The two branches share a separation epoch and continue independently.\")\n", "print(f\"Separation: {event_times['stage_separation']:.6f} s\")\n", "print(f\"Upper cutoff: {orbit['cutoff_time_s']:.6f} s\")\n", "print(f\"Insertion apsides: {orbit['perigee_km']:.6f} × {orbit['apogee_km']:.6f} km\")\n", "print(f\"Payload release: {event_times['payload_release']:.6f} s\")\n", "print(f\"Verified orbit: {verified['propagated_duration_s']:.6f} s\")\n", "print(f\"Payload altitude range: {verified['min_altitude_km']:.6f}–{verified['max_altitude_km']:.6f} km\")\n", "print(f\"Booster touchdown: {landing['touchdown_time_s']:.6f} s\")\n", "print(f\"Contact speed: {landing['speed_m_s']:.6f} m/s; miss: {landing['miss_distance_m']:.6f} m\")\n", "print(f\"Recovery propellant remaining: {landing['propellant_remaining_kg']:.6f} kg\")\n", "print(\"Integrated insertion, payload orbit and booster landing checks passed.\")\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "mission-sources", "metadata": {}, "source": [ "## Sources\n", "\n", "[Rocket Lab — Neutron](https://rocketlabcorp.com/launch/neutron/): Published vehicle dimensions, thrust and recovery architecture; internal masses and guidance remain simulation assumptions.\n", "\n", "[SciPy RK45](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html): Adaptive Dormand–Prince integration.\n", "\n", "[SciPy solve_ivp](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html): Dense output and terminal-event root detection." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.11" }, "missionForwardPass": { "modelSha256": "ec34d96c89016df203ebc80120e57946d56f96a94820f3d2fc81e48fd28d6e66", "revision": "5ce90245a522", "stage1Revision": "aad264bd8eb4" } }, "nbformat": 4, "nbformat_minor": 5 }