"""Reproduce the bounded post-combustor energy/reinjection example.
Standard library only. This is a synthetic teaching model, not an engine
simulation or measurement. Run from any directory with Python 3.
"""
from __future__ import annotations
import csv
import json
import math
from pathlib import Path
from xml.etree import ElementTree
ROOT = Path(__file__).resolve().parent
MASS_FLOW_KG_S = 1.0
BASELINE_EXIT_SPEED_M_S = 2_000.0
INCOMING_SPEED_M_S = 1_500.0
GENERATOR_EFFICIENCY = 0.60
AUXILIARY_LOAD_KW = 12.0
ACCELERATOR_SHARE = 0.70
ACCELERATOR_EFFICIENCY = 0.60
def calculate(extraction_fraction: float) -> dict[str, float]:
if not 0 <= extraction_fraction <= 0.8:
raise ValueError("extraction_fraction must be between 0 and 0.8")
initial_kw = 0.5 * MASS_FLOW_KG_S * BASELINE_EXIT_SPEED_M_S**2 / 1000
removed_kw = extraction_fraction * initial_kw
generated_kw = GENERATOR_EFFICIENCY * removed_kw
generator_loss_kw = removed_kw - generated_kw
served_aux_kw = min(generated_kw, AUXILIARY_LOAD_KW)
unmet_aux_kw = AUXILIARY_LOAD_KW - served_aux_kw
bus_kw = generated_kw - served_aux_kw
accelerator_input_kw = ACCELERATOR_SHARE * bus_kw
returned_kw = ACCELERATOR_EFFICIENCY * accelerator_input_kw
accelerator_loss_kw = accelerator_input_kw - returned_kw
export_kw = bus_kw - accelerator_input_kw
final_stream_kw = initial_kw - removed_kw + returned_kw
final_speed = math.sqrt(2 * final_stream_kw * 1000 / MASS_FLOW_KG_S)
baseline_thrust_n = MASS_FLOW_KG_S * (
BASELINE_EXIT_SPEED_M_S - INCOMING_SPEED_M_S
)
final_thrust_n = MASS_FLOW_KG_S * (final_speed - INCOMING_SPEED_M_S)
accounted_kw = (
final_stream_kw + generator_loss_kw + served_aux_kw
+ accelerator_loss_kw + export_kw
)
return {
"extraction_fraction": extraction_fraction,
"initial_stream_kinetic_power_kW": initial_kw,
"removed_stream_power_kW": removed_kw,
"generator_electric_output_kW": generated_kw,
"generator_loss_kW": generator_loss_kw,
"auxiliary_load_served_kW": served_aux_kw,
"auxiliary_load_unmet_kW": unmet_aux_kw,
"accelerator_electrical_input_kW": accelerator_input_kw,
"accelerator_power_returned_to_stream_kW": returned_kw,
"accelerator_loss_kW": accelerator_loss_kw,
"electrical_export_kW": export_kw,
"final_stream_kinetic_power_kW": final_stream_kw,
"final_exit_speed_m_s": final_speed,
"baseline_momentum_thrust_N": baseline_thrust_n,
"final_momentum_thrust_N": final_thrust_n,
"momentum_thrust_change_N": final_thrust_n - baseline_thrust_n,
"energy_accounted_kW": accounted_kw,
}
def make_figure(rows: list[dict[str, float]], target: Path) -> None:
width, height = 1040, 740
left, right, top, middle, bottom = 100, 990, 120, 365, 610
x = lambda f: left + (right - left) * f / 0.8
thrust_min = math.floor(min(0, min(row["momentum_thrust_change_N"] for row in rows)) / 100) * 100
y_thrust = lambda value: top + (middle - top) * (0 - value) / (-thrust_min)
y_export = lambda value: middle + (bottom - middle) * (300 - value) / 300
thrust_points = " ".join(
f"{x(row['extraction_fraction']):.2f},{y_thrust(row['momentum_thrust_change_N']):.2f}"
for row in rows
)
export_points = " ".join(
f"{x(row['extraction_fraction']):.2f},{y_export(row['electrical_export_kW']):.2f}"
for row in rows
)
grid = []
for value in range(0, int(thrust_min) - 1, -200):
yy = y_thrust(value)
grid.append(f'{value}')
for value in (0, 100, 200, 300):
yy = y_export(value)
grid.append(f'{value}')
for frac in (0, .2, .4, .6, .8):
xx = x(frac)
grid.append(f'{frac:.0%}')
svg = f''''''
target.write_text(svg, encoding="utf-8")
def make_mobile_figure(rows: list[dict[str, float]], target: Path) -> None:
width, height = 460, 820
left, right, top, middle, bottom = 74, 430, 150, 400, 650
x = lambda f: left + (right - left) * f / 0.8
thrust_min = math.floor(min(0, min(row["momentum_thrust_change_N"] for row in rows)) / 100) * 100
y_thrust = lambda value: top + (middle - top) * (0 - value) / (-thrust_min)
y_export = lambda value: middle + (bottom - middle) * (300 - value) / 300
thrust_points = " ".join(
f"{x(row['extraction_fraction']):.2f},{y_thrust(row['momentum_thrust_change_N']):.2f}"
for row in rows
)
export_points = " ".join(
f"{x(row['extraction_fraction']):.2f},{y_export(row['electrical_export_kW']):.2f}"
for row in rows
)
elements = []
for value in range(0, int(thrust_min) - 1, -200):
yy = y_thrust(value)
elements.append(f'{value}')
for value in (0, 100, 200, 300):
yy = y_export(value)
elements.append(f'{value}')
for frac in (0, .2, .4, .6, .8):
xx = x(frac)
elements.append(f'{frac:.0%}')
svg = f''''''
target.write_text(svg, encoding="utf-8")
def main() -> None:
rows = [calculate(i / 100) for i in range(81)]
with (ROOT / "coupled-ledger.csv").open("w", newline="", encoding="utf-8") as f:
writer = csv.DictWriter(f, fieldnames=list(rows[0]))
writer.writeheader()
writer.writerows(rows)
make_figure(rows, ROOT / "coupled-ledger.svg")
make_mobile_figure(rows, ROOT / "coupled-ledger-mobile.svg")
checks = 0
for row in rows:
assert abs(row["energy_accounted_kW"] - row["initial_stream_kinetic_power_kW"]) < 1e-9
assert row["generator_loss_kW"] >= -1e-12
assert row["accelerator_loss_kW"] >= -1e-12
assert row["auxiliary_load_unmet_kW"] >= -1e-12
assert row["electrical_export_kW"] >= -1e-12
assert row["final_momentum_thrust_N"] == (
MASS_FLOW_KG_S * (row["final_exit_speed_m_s"] - INCOMING_SPEED_M_S)
)
checks += 6
assert rows[0]["momentum_thrust_change_N"] == 0
assert rows[0]["electrical_export_kW"] == 0
assert rows[0]["final_exit_speed_m_s"] == BASELINE_EXIT_SPEED_M_S
checks += 3
# Independent reference arithmetic for the published slider defaults.
ref = calculate(.20)
assert abs(ref["generator_electric_output_kW"] - 240) < 1e-10
assert abs(ref["accelerator_power_returned_to_stream_kW"] - 95.76) < 1e-10
assert abs(ref["electrical_export_kW"] - 68.4) < 1e-10
assert abs(ref["final_exit_speed_m_s"] - math.sqrt(2 * 1695.76 * 1000)) < 1e-10
assert abs(ref["momentum_thrust_change_N"] - (ref["final_exit_speed_m_s"] - 2000)) < 1e-10
checks += 5
assert all(rows[i + 1]["electrical_export_kW"] >= rows[i]["electrical_export_kW"] for i in range(80))
assert all(rows[i + 1]["momentum_thrust_change_N"] <= rows[i]["momentum_thrust_change_N"] for i in range(80))
checks += 2
# Independent identity: with fixed mass flow and incoming momentum, the
# momentum-contribution change equals mdot times the exit-speed change.
for row in rows:
independent_delta = MASS_FLOW_KG_S * (
row["final_exit_speed_m_s"] - BASELINE_EXIT_SPEED_M_S
)
assert abs(row["momentum_thrust_change_N"] - independent_delta) < 1e-9
checks += 1
# Inspect the actual generated SVG polylines, not only the plot mapping.
for filename, top, bottom in (
("coupled-ledger.svg", 120, 365),
("coupled-ledger-mobile.svg", 150, 400),
):
root = ElementTree.parse(ROOT / filename).getroot()
thrust_line = root.findall("{http://www.w3.org/2000/svg}polyline")[0]
points = [tuple(map(float, point.split(","))) for point in thrust_line.attrib["points"].split()]
assert len(points) == len(rows)
assert min(y for _, y in points) >= top - 1e-9
assert max(y for _, y in points) <= bottom + 1e-9
checks += 3
methods = {
"classification": "synthetic analytical teaching model; no measurements; no uncertainty distribution",
"boundary": "After-combustor baseline nozzle-exit station; compare stream kinetic-energy and momentum-flux changes at fixed mass flow, incoming speed, and pressure-thrust contribution. This is a steady, instantaneous vehicle-frame teaching convention, not a vehicle-dynamics model.",
"velocity_convention": "V_infinity and V_e are positive axial flow-speed magnitudes at the inlet and nozzle-exit stations, expressed in one instantaneous vehicle-fixed frame treated as steady for this calculation.",
"thrust_sign_convention": "F_momentum = mass_flow * (V_e - V_infinity) is reported as the forward-thrust contribution on the vehicle. A negative delta means this contribution decreases relative to the baseline. It is not the signed force on the gas.",
"mass_flow_assumption": "The same selected mass flow is imposed at both comparison stations and in both cases. Fuel mass addition, bleed, capture/spillage changes, and flow-rate changes are not modeled.",
"pressure_term_assumption": "The nozzle pressure-thrust term is set to zero for the displayed absolute example and held unchanged when calculating its delta. This is an imposed comparison assumption, not a consequence of energy closure.",
"fixed_reference_case": {
"mass_flow_kg_s": MASS_FLOW_KG_S,
"baseline_exit_speed_m_s": BASELINE_EXIT_SPEED_M_S,
"incoming_speed_m_s": INCOMING_SPEED_M_S,
"generator_efficiency_assumption": GENERATOR_EFFICIENCY,
"auxiliary_load_kW": AUXILIARY_LOAD_KW,
"accelerator_bus_share": ACCELERATOR_SHARE,
"accelerator_efficiency_assumption": ACCELERATOR_EFFICIENCY,
},
"equations": {
"initial_stream_kinetic_power_kW": "0.5 * mass_flow_kg_s * baseline_exit_speed_m_s^2 / 1000",
"generator_electric_output_kW": "generator_efficiency * extraction_fraction * initial_stream_kinetic_power_kW",
"generator_loss_kW": "removed_stream_power_kW - generator_electric_output_kW",
"served_auxiliary_load_kW": "min(generator_electric_output_kW, auxiliary_load_kW)",
"accelerator_power_returned_to_stream_kW": "accelerator_efficiency * accelerator_bus_share * max(0, generator_electric_output_kW - served_auxiliary_load_kW); an assumed lumped electric-to-directed-stream-kinetic-power conversion, not a modeled electromagnetic force field",
"final_exit_speed_m_s": "sqrt(2 * (initial_stream_kinetic_power_kW - removed_stream_power_kW + accelerator_power_returned_to_stream_kW) * 1000 / mass_flow_kg_s)",
"momentum_thrust_N": "forward-thrust contribution on the vehicle = mass_flow_kg_s * (positive downstream exit-speed magnitude - positive downstream inlet-speed magnitude); the displayed absolute example sets pressure thrust to zero, and its change imposes an unchanged pressure term",
"energy_identity": "final_stream_kinetic_power + generator_loss + served_auxiliary_load + accelerator_loss + electrical_export = initial_stream_kinetic_power",
},
"sweep": {"extraction_fraction_min": 0, "extraction_fraction_max": .8, "step": .01, "rows": len(rows)},
"limitations": [
"Not a complete vehicle, inlet, combustor, nozzle, MHD, or engine model; it is a defined same-stream accounting example.",
"No measurements support slider values or selected efficiencies; they are editable teaching assumptions.",
"Fixed mass flow, incoming speed, reference frame, baseline nozzle exit station, and pressure-thrust contribution. No added fuel mass, bleed, capture/spillage change, or changed mass flow is represented.",
"Velocity values are positive axial speed magnitudes in a steady, instantaneous vehicle-frame convention. This is not a vehicle-dynamics model; a negative delta means the forward momentum-thrust contribution on the vehicle is lower than baseline.",
"Pressure-thrust cancellation is an imposed assumption. Actual exit pressure, exit area, nozzle rematching/choking, and coupled pressure change are not calculated.",
"No inlet capture or spillage, pressure drag, shock/viscous/thermal losses, conductivity, magnet mass, cooling, combustion matching, installation mass, aircraft drag, range, or fuel use.",
"Generator and accelerator conversion losses, served auxiliaries, and electrical export are treated as sinks outside the tracked stream; dissipated heat and later interaction with the gas/vehicle are not followed.",
"Does not represent the distinct NASA generator-before-combustor bypass cycle and does not evidence a tested integrated loop.",
"Electrical export is energy left at the bus, not a demonstrated vehicle benefit; thrust change is only the defined momentum-flux difference.",
"Accelerator power return is a lumped assumption that electrical input becomes directed kinetic power in this same stream. Electrode geometry, conductivity, current closure, electromagnetic body force, field generation/self-excitation, and coupling are not modeled or demonstrated.",
],
"verification": {"passed": True, "checks": checks, "method": "81-point sweep; dimensional energy closure; nonnegative loss/load/export checks; independent fixed-inlet momentum-delta identity; generated desktop/mobile SVG thrust-curve bounds; reference and endpoint/monotonicity checks."},
}
(ROOT / "coupled-methods.json").write_text(json.dumps(methods, indent=2) + "\n", encoding="utf-8")
print(json.dumps({"rows": len(rows), "checks_passed": checks, "reference": ref}, indent=2))
if __name__ == "__main__":
main()