{ "cells": [ { "attachments": {}, "cell_type": "markdown", "id": "stage1-introduction", "metadata": {}, "source": [ "# Stage one, from ignition to separation\n", "\n", "A complete forward pass through the simulator: initialize one state, turn it into forces and rates, advance it, stop at the fuel reserve, coast, then separate. Every worked number is computed from the same Python physics functions used by the simulator.\n", "\n", "**Scope.** Educational planar point-mass model with ideal pointing, spherical nonrotating Earth and a still atmosphere. Vehicle internal masses, specific impulse, drag coefficient and pitch program are assumptions. This notebook ends at separation. All internal calculations use SI units.\n", "\n", "Read the eight steps in order; executable Python follows in the appendix. Use **Run All** in Jupyter after installing NumPy 2.4.2 and SciPy 1.17.0 (Python 3.11 or newer)." ] }, { "attachments": {}, "cell_type": "markdown", "id": "stage1-initial", "metadata": {}, "source": [ "## 1 · 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": [ "## 2 · 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": [ "## 3 · 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": [ "## 4 · 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": [ "## 5 · 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": [ "## 6 · 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": [ "## 7 · 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", "\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": [ "## 8 · 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", "This is the end of the stage-one ascent forward pass. The upper assembly appears only as carried mass and the other side of the separation impulse. Recovery and upper-stage flight are outside this notebook. 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": "stage1-results", "metadata": {}, "source": [ "## The computed path\n", "\n", "Downrange is arc length $x=R_E\\theta$. These are sampled integrated states.\n", "\n", "| Time (s) | Altitude (km) | Downrange (km) | Radial speed (m/s) | Tangential speed (m/s) | Mass (kg) |\n", "| ---: | ---: | ---: | ---: | ---: | ---: |\n", "| 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 480000.000 |\n", "| 40.000 | 3.484 | 0.330 | 187.300 | 36.113 | 400277.705 |\n", "| 80.000 | 15.845 | 5.384 | 442.228 | 255.053 | 320555.409 |\n", "| 120.000 | 40.286 | 24.144 | 802.734 | 742.602 | 240833.114 |\n", "| 137.979 | 56.588 | 40.314 | 1017.239 | 1084.841 | 205000.000 |\n", "| 141.979 | 60.581 | 44.612 | 979.244 | 1083.905 | 205000.000 |" ] }, { "attachments": {}, "cell_type": "markdown", "id": "stage1-appendix", "metadata": {}, "source": [ "## Appendix · run the forward pass\n", "\n", "The next cell copies each required canonical model definition exactly once. The following cell supplies only stage-one phase orchestration and trace output. There are no upper-flight or recovery definitions. Canonical model SHA-256: `ec34d96c89016df203ebc80120e57946d56f96a94820f3d2fc81e48fd28d6e66`.\n", "\n", "The final cell runs the reference case and checks analytic burn time plus mass and momentum at separation." ] }, { "cell_type": "code", "execution_count": 1, "id": "stage1-canonical", "metadata": { "source_id": "stage1-canonical" }, "outputs": [], "source": [ "MODEL_SHA256 = \"ec34d96c89016df203ebc80120e57946d56f96a94820f3d2fc81e48fd28d6e66\"\n", "\n", "from dataclasses import asdict, dataclass\n", "\n", "\n", "import argparse\n", "\n", "\n", "import json\n", "\n", "\n", "import math\n", "\n", "\n", "import numpy as np\n", "\n", "\n", "from scipy.integrate import solve_ivp\n", "\n", "\n", "MU = 3.986004418e14 # Earth gravitational parameter [m^3/s^2]\n", "\n", "\n", "EARTH_RADIUS = 6_378_137.0 # spherical Earth radius [m]\n", "\n", "\n", "G0 = 9.80665 # Isp reference gravity [m/s^2]\n", "\n", "\n", "SOURCE = \"https://rocketlabcorp.com/launch/neutron/\"\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 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 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" ] }, { "cell_type": "code", "execution_count": 2, "id": "stage1-runner", "metadata": { "source_id": "stage1-runner" }, "outputs": [], "source": [ "def stage1_sample(time, state, phase):\n", " \"\"\"Display units are derived from the five-state solution, never reintegrated.\"\"\"\n", " radius, angle, vr, vt, mass = map(float, state)\n", " return dict(time_s=float(time), phase=phase,\n", " altitude_m=radius - EARTH_RADIUS, downrange_m=EARTH_RADIUS * angle,\n", " radial_velocity_m_s=vr, tangential_velocity_m_s=vt, mass_kg=mass)\n", "\n", "\n", "def run_stage1(payload_kg=8000.0, recovery_propellant_kg=60_000.0,\n", " max_step_s=2.0, rtol=1e-7, sample_step_s=2.0):\n", " \"\"\"The canonical stack's launch, powered ascent, four-second coast and split.\"\"\"\n", " vehicle = Vehicle()\n", " if not all(math.isfinite(value) for value in\n", " (payload_kg, recovery_propellant_kg, max_step_s, rtol, sample_step_s)):\n", " raise ValueError(\"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 min(max_step_s, rtol, sample_step_s) <= 0:\n", " raise ValueError(\"Solver and sampling controls must be positive.\")\n", "\n", " area = math.pi * vehicle.diameter_m**2 / 4\n", " upper_wet = vehicle.upper_dry_kg + vehicle.upper_propellant_kg + payload_kg\n", " initial_mass = vehicle.booster_dry_kg + vehicle.booster_propellant_kg + upper_wet\n", " reserve_mass = vehicle.booster_dry_kg + recovery_propellant_kg\n", " initial = np.array([EARTH_RADIUS, 0.0, 0.0, 0.0, initial_mass])\n", " powered = lambda t, y: ascent_guidance(t, y, vehicle)\n", " coast = lambda t, y: np.zeros(2)\n", " impact = terminal_event(lambda t, y: y[0] - EARTH_RADIUS)\n", " meco = terminal_event(lambda t, y: y[4] - upper_wet - reserve_mass)\n", "\n", " ascent = integrate(initial, 0, 600, powered, area, vehicle.drag_coefficient,\n", " vehicle.booster_isp_s, (impact, meco), max_step_s, rtol)\n", " if ascent.t_events[0].size or not ascent.t_events[1].size:\n", " raise RuntimeError(\"Stage one did not reach the reserve-mass MECO gate.\")\n", " meco_time = float(ascent.t[-1])\n", " meco_state = ascent.y[:, -1].copy()\n", " opening = integrate(meco_state, meco_time, meco_time + 4, coast, area,\n", " vehicle.drag_coefficient, vehicle.booster_isp_s,\n", " (impact,), max_step_s, rtol)\n", " if opening.t_events[0].size:\n", " raise RuntimeError(\"Ground contact interrupted the fairing-opening coast.\")\n", " separation_time = float(opening.t[-1])\n", " separation_state = opening.y[:, -1].copy()\n", " booster, upper = separate(separation_state, reserve_mass)\n", "\n", " samples = []\n", " for solution, name in ((ascent, \"ascent\"), (opening, \"fairing_opening\")):\n", " start, stop = float(solution.t[0]), float(solution.t[-1])\n", " times = np.unique(np.r_[start, np.arange(start, stop, sample_step_s), stop])\n", " for time, state in zip(times, solution.sol(times).T):\n", " sample = stage1_sample(time, state, name)\n", " sample[\"fairing_open_fraction\"] = (0.0 if name == \"ascent\" else\n", " float(np.clip((time - meco_time) / 4, 0, 1)))\n", " samples.append(sample)\n", "\n", " first_state = ascent.y[:, 1].copy()\n", " first_step = dict(stage1_sample(ascent.t[1], first_state, \"ascent\"), state=first_state.tolist())\n", " rates = np.array(rhs(0, initial, powered, area, vehicle.drag_coefficient, vehicle.booster_isp_s))\n", " milestones = []\n", " for identifier, label, time, state, phase in (\n", " (\"launch\", \"Launch\", 0.0, initial, \"ascent\"),\n", " (\"first-step\", \"First accepted RK45 step\", ascent.t[1], first_state, \"ascent\"),\n", " (\"meco\", \"Main engine cutoff\", meco_time, meco_state, \"ascent\"),\n", " (\"separation\", \"Just before separation\", separation_time, separation_state, \"fairing_opening\"),\n", " ):\n", " milestones.append(dict(stage1_sample(time, state, phase), id=identifier,\n", " label=label, state=state.tolist()))\n", "\n", " return dict(\n", " modelSha256=MODEL_SHA256,\n", " initialValues=dict(payload_kg=payload_kg, recovery_propellant_kg=recovery_propellant_kg,\n", " upper_wet_mass_kg=upper_wet, initial_mass_kg=initial_mass,\n", " reserve_mass_kg=reserve_mass, cutoff_mass_kg=reserve_mass + upper_wet,\n", " earth_radius_m=EARTH_RADIUS, mu_m3_s2=MU, g0_m_s2=G0,\n", " booster_thrust_n=vehicle.booster_thrust_n, booster_isp_s=vehicle.booster_isp_s,\n", " booster_dry_kg=vehicle.booster_dry_kg,\n", " booster_propellant_kg=vehicle.booster_propellant_kg,\n", " upper_dry_kg=vehicle.upper_dry_kg, upper_propellant_kg=vehicle.upper_propellant_kg,\n", " diameter_m=vehicle.diameter_m, drag_coefficient=vehicle.drag_coefficient,\n", " state=initial.tolist()),\n", " firstEvaluation=dict(thrust_n=float(np.linalg.norm(powered(0, initial))),\n", " density_kg_m3=atmosphere(0), area_m2=area,\n", " gravity_m_s2=MU / EARTH_RADIUS**2,\n", " radial_acceleration_m_s2=float(rates[2]),\n", " mass_flow_kg_s=float(rates[4]), derivative=rates.tolist()),\n", " firstAcceptedStep=first_step,\n", " mecoTime=meco_time, separationTime=separation_time,\n", " milestones=milestones, samples=samples,\n", " mecoState=meco_state.tolist(), separationState=separation_state.tolist(),\n", " boosterState=booster.tolist(), upperState=upper.tolist(),\n", " solver=dict(method=\"RK45\", rtol=rtol, atol=[1e-3, 1e-11, 1e-5, 1e-5, 1e-4],\n", " max_step_s=max_step_s, sample_step_s=sample_step_s,\n", " accepted_powered_steps=len(ascent.t) - 1,\n", " accepted_coast_steps=len(opening.t) - 1,\n", " function_evaluations=ascent.nfev + opening.nfev))\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "stage1-run", "metadata": { "source_id": "stage1-run" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Stage one: launch → MECO → four-second coast → separation\n", "MECO: 137.978968 s; separation: 141.978968 s\n", "Time (s) Altitude (km) Downrange (km) vr (m/s) vt (m/s) Mass (kg)\n", " 0.000000 0.000000 0.000000 0.000000 0.000000 480000.000\n", " 0.036069 0.000002 0.000000 0.131295 0.000000 479928.111\n", "137.978968 56.587922 40.314212 1017.238595 1084.841261 205000.000\n", "141.978968 60.580828 44.612168 979.243560 1083.905092 205000.000\n", "Analytical burn time, mass and momentum checks passed.\n" ] } ], "source": [ "stage1 = run_stage1()\n", "print(\"Stage one: launch → MECO → four-second coast → separation\")\n", "print(f\"MECO: {stage1['mecoTime']:.6f} s; separation: {stage1['separationTime']:.6f} s\")\n", "print(\"Time (s) Altitude (km) Downrange (km) vr (m/s) vt (m/s) Mass (kg)\")\n", "for row in stage1[\"milestones\"]:\n", " print(f\"{row['time_s']:10.6f} {row['altitude_m']/1000:14.6f} {row['downrange_m']/1000:17.6f} \"\n", " f\"{row['radial_velocity_m_s']:11.6f} {row['tangential_velocity_m_s']:11.6f} {row['mass_kg']:12.3f}\")\n", "before, booster, upper = (np.array(stage1[key]) for key in (\"separationState\", \"boosterState\", \"upperState\"))\n", "assert abs(booster[4] + upper[4] - before[4]) < 1e-8\n", "np.testing.assert_allclose(booster[4] * booster[2:4] + upper[4] * upper[2:4], before[4] * before[2:4], rtol=1e-14)\n", "assert abs(upper[2] - booster[2] - 0.5) < 1e-10\n", "analytic_meco = ((stage1[\"initialValues\"][\"initial_mass_kg\"] - stage1[\"initialValues\"][\"cutoff_mass_kg\"])\n", " / -stage1[\"firstEvaluation\"][\"mass_flow_kg_s\"])\n", "assert abs(stage1[\"mecoTime\"] - analytic_meco) < 1e-8\n", "print(\"Analytical burn time, mass and momentum checks passed.\")\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "stage1-sources", "metadata": {}, "source": [ "## Sources\n", "\n", "[Rocket Lab — Neutron](https://rocketlabcorp.com/launch/neutron/): published dimensions, thrust and lift-off mass.\n", "\n", "[SciPy RK45](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html): Dormand–Prince integration.\n", "\n", "[SciPy solve_ivp](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html): tolerances, dense output and terminal roots.\n", "\n", "[SciPy RK implementation](https://github.com/scipy/scipy/blob/v1.17.0/scipy/integrate/_ivp/rk.py): fixed tableau and acceptance rule." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.11" }, "stage1": { "modelSha256": "ec34d96c89016df203ebc80120e57946d56f96a94820f3d2fc81e48fd28d6e66", "revision": "aad264bd8eb4" } }, "nbformat": 4, "nbformat_minor": 5 }