Neutron.
One forward pass.

Start with five numbers

At time $t$ in seconds, the state $\mathbf y$ is the five numbers carried from one calculation to the next. Subscript $0$ means liftoff: the rocket starts at rest with full thrust applied immediately.

The three terms are stage 1 dry mass, stage 1 propellant, and the carried assembly. That last 115,000 kg includes its structure, fuel and an 8,000 kg payload. These internal masses are model assumptions.

Carry forward: the current time and this five-number state.

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.

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.

MU = 3.986004418e14              # Earth gravitational parameter [m^3/s^2]
EARTH_RADIUS = 6_378_137.0       # spherical Earth radius [m]
G0 = 9.80665                    # Isp reference gravity [m/s^2]

@dataclass(frozen=True)
class Vehicle:
    """Published geometry/thrust; explicitly assumed internal mass and Isp."""
    height_m: float = 43.0
    diameter_m: float = 7.0
    published_liftoff_mass_kg: float = 480_000.0
    advertised_leo_payload_kg: float = 13_000.0
    booster_engines: int = 9
    booster_thrust_n: float = 1_450_000.0 * 4.4482216152605
    upper_thrust_n: float = 890_000.0
    booster_dry_kg: float = 30_000.0   # assumption, includes captive fairing
    booster_propellant_kg: float = 335_000.0
    booster_isp_s: float = 330.0
    upper_dry_kg: float = 5_000.0
    upper_propellant_kg: float = 102_000.0
    upper_isp_s: float = 365.0
    drag_coefficient: float = 0.35

    area = math.pi * vehicle.diameter_m**2 / 4
    upper_wet = vehicle.upper_dry_kg + vehicle.upper_propellant_kg * upper_propellant_scale + payload_kg
    initial_mass = vehicle.booster_dry_kg + vehicle.booster_propellant_kg + upper_wet
    initial = np.array([EARTH_RADIUS, 0.0, 0.0, 0.0, initial_mass])

Point the thrust

The pitch angle $\beta$ is measured above the local horizontal. The model linearly interpolates this assumed schedule and holds the endpoint pitch outside its time range:

$F_r,F_t$ are thrust components (N). The total thrust $T$ stays at 6.4499 MN during ascent. The implementation applies the commanded direction instantly.

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.

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.

def ascent_guidance(time, state, vehicle):
    """An assumed pitch program; angles measured above local horizontal."""
    pitch = np.interp(time, [0, 12, 35, 70, 115, 160], [90, 90, 82, 67, 52, 40])
    beta = math.radians(pitch)
    return vehicle.booster_thrust_n * np.array([math.sin(beta), math.cos(beta)])

Find the air, then the drag

$h$ is altitude (m), $R_E=6\,378\,137\ \mathrm m$ is Earth’s radius, and $V$ is speed through the still atmosphere (m/s).

$\rho$ is air density (kg/m³). The assumed atmosphere has sea-level density 1.225 kg/m³ and an 8,500 m scale height.

$D_r,D_t$ are drag components (N). The assumed drag coefficient is $C_D=0.35$. The reference area is $A=\pi(7\,\mathrm m)^2/4$ = 38.48 m². Negative signs make drag oppose velocity.

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.

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.

def atmosphere(altitude_m):
    """One-layer density model, not a weather or aerothermal model."""
    return 1.225 * math.exp(-max(0.0, altitude_m) / 8500.0)

def drag_force(state, area_m2, cd):
    """Components in the local outward/tangential frame; still atmosphere."""
    radius, _, vr, vt, _ = state
    speed = math.hypot(vr, vt)
    scale = -0.5 * atmosphere(radius - EARTH_RADIUS) * cd * area_m2 * speed
    return np.array([scale * vr, scale * vt])

Calculate all five rates of change

A dot means “change per second.” Insert the thrust and drag just computed into the following equations. Together they form $\dot{\mathbf y}=f(t,\mathbf y)$, the function evaluated by the solver.

Carry forward: five derivatives, all evaluated from the same state. No state component has been advanced yet.

Scroll the diagram sideways to inspect both ends.

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.

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.

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.

def rhs(time, state, guidance, area_m2, cd, isp_s):
    """Polar Newton equations; state = [r, theta, v_r, v_t, mass].

    Pointing is ideal: guidance commands force directly, not attitude/torque.
    The changing mass appears in F/m and dm/dt; thrust includes exhaust momentum.
    """
    radius, _, vr, vt, mass = state
    force = guidance(time, state)
    radial, tangential = (force + drag_force(state, area_m2, cd)) / mass
    return [vr, vt / radius,
            vt * vt / radius - MU / radius**2 + radial,
            -vr * vt / radius + tangential,
            -np.linalg.norm(force) / (isp_s * G0)]

Advance the state, then repeat

The solver uses Dormand–Prince RK45. Within a trial step $\Delta t$ (s), it recalculates forces and derivatives at intermediate times and states, producing slopes $\mathbf k_i$.

$n$ counts accepted steps; $i$ counts solver stages and $j$ earlier stages. $c_i,a_{ij}$ are fixed RK45 coefficients; $\mathbf y^{[i]}$ is a trial state.

$b_i$ are the fifth-order weights. A second set of weights gives a fourth-order estimate. If their scaled difference is small enough, accept $\mathbf y_{n+1}=\mathbf y^{(5)}$ and $t_{n+1}=t_n+\Delta t$. Otherwise shrink the step and retry from $\mathbf y_n$.

Repeat this force-to-state calculation while monitoring the cutoff condition below.

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.

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.

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.

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.

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.

def integrate(state, start_s, end_s, guidance, area_m2, cd, isp_s,
              events=(), max_step_s=2.0, rtol=1e-7):
    """Integrate one continuous phase; restart after every discrete change."""
    solution = solve_ivp(
        lambda t, y: rhs(t, y, guidance, area_m2, cd, isp_s),
        (start_s, end_s), state, method="RK45", dense_output=True,
        events=events, max_step=max_step_s, rtol=rtol,
        atol=[1e-3, 1e-11, 1e-5, 1e-5, 1e-4])
    if not solution.success:
        raise RuntimeError(solution.message)
    return solution

Stop burning at the fuel reserve

Main engine cutoff (MECO) occurs when the stack reaches the carried assembly’s mass plus stage 1 dry mass and its unburned reserve.

$m_{\rm cut}$ is the cutoff threshold. The 60,000 kg reserve is a chosen model input.

$g$ is the event function. SciPy locates its zero crossing within the accepted step and ends powered integration at that time, $t_c$. The arrow $\downarrow$ means a crossing from positive to negative; $\uparrow$ means the reverse.

The integration also stops on downward ground contact, or at 600 s if cutoff is never reached. Those paths do not proceed to a successful separation.

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.

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.

def terminal_event(function, direction=-1):
    """SciPy locates a zero crossing between accepted RK45 steps."""
    function.terminal = True
    function.direction = direction
    return function

    coast = lambda t, y: np.zeros(2)
    reserve_mass = vehicle.booster_dry_kg + recovery_propellant_kg
    meco = terminal_event(lambda t, y: y[4] - upper_wet - reserve_mass)
    record("liftoff", 0, "stack", "Nine-engine ascent begins.")
    time, state, impacted, ascent_solution = phase("stack", "ascent", initial, 0, 600,
                                                  lambda t, y: ascent_guidance(t, y, vehicle), (meco,))
    ascent_complete = bool(ascent_solution.t_events[1].size) and not impacted

Coast for four seconds

$t_s$ is separation time, written as $t_{sep}$ in the branches below. Restart the same five equations from the cutoff state with zero thrust. Gravity and drag still change the velocities; mass remains 205,000 kg. Downward ground contact remains a terminal event.

The captive fairing opens during these four seconds. Its opening is a kinematic display model; it changes neither drag area nor mass in the equations.

Carry forward: the final attached-stack state, ready for an instantaneous mass split.

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.

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.

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.

    if ascent_complete:
        record("meco", time, "stack", "Main engine cutoff preserves the requested recovery reserve.",
               recovery_propellant_kg=float(recovery_propellant_kg))
        record("fairing_opening", time, "stack", "Captive fairing opens over four seconds (kinematic mechanism).")
        time, state, impacted, _ = phase("stack", "fairing_opening", state, time, time + 4, coast)

    def fairing_fraction(time, body):
        opening = next((e["time_s"] for e in events if e["name"] == "fairing_opening"), None)
        if opening is None or body not in ("stack", "booster"):
            return 0.0
        if separation_time is None or time <= separation_time:
            return float(np.clip((time - opening) / 4, 0, 1))
        return float(np.clip(1 - (time - separation_time) / 4, 0, 1))

Separate with momentum conserved

Split the final mass into the booster $m_b$ and the released assembly $m_u$. The fairing stays with stage 1.

The assumed separation speed $\Delta v=0.5\ \mathrm{m/s}$ is the relative outward velocity between the two bodies. Superscript $-$ means just before separation; $+$ means just after.

Both bodies inherit the same $r,\theta,v_t$. The impulses change only their radial velocities. There is no artificial position jump.

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.

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.

def separate(state, booster_mass_kg, relative_speed_m_s=0.5):
    """Split mass with equal/opposite impulses: conserve linear momentum.

    Both bodies start at the same position. Geometry and contact are not modeled.
    The reusable fairing stays within booster dry mass.
    """
    booster, upper = state.copy(), state.copy()
    upper[4] = state[4] - booster_mass_kg
    booster[4] = booster_mass_kg
    booster[2] -= relative_speed_m_s * upper[4] / state[4]
    upper[2] += relative_speed_m_s * booster[4] / state[4]
    return booster, upper