Aerospace starter - Nozzle branches

MODEL "Aerospace starter - Nozzle branches"
DESCRIPTION "Both branches of the isentropic area-Mach relation, with a sonic endpoint and input contracts."
VERSION "1.0.0"
TAGS "aerospace", "starter", "propulsion", "roots"

strict

input AreaRatio = 2.5
  VALIDATE BETWEEN 1 AND 10 MESSAGE "Area/throat area must be at least one; this starter accepts 1 to 10."
input Gamma = 1.4
  VALIDATE BETWEEN 1.1 AND 1.67 MESSAGE "Use a constant heat-capacity ratio from 1.1 to 1.67."
input StagnationTemperature_K = 900
  VALIDATE BETWEEN 100 AND 3000 MESSAGE "Use a positive absolute temperature from 100 to 3000 K."
input StagnationPressure_Pa = 500000
  VALIDATE BETWEEN 1000 AND 100000000 MESSAGE "Use absolute pressure from 1000 to 100,000,000 Pa."
GasConstant_J_kg_K = 287.05

AreaMach = (mach) => (1 / mach) * ((2 / (Gamma + 1)) * (1 + (Gamma - 1) * mach^2 / 2))^((Gamma + 1) / (2 * (Gamma - 1)))
# The endpoint M=1 belongs to both branches. Invalid source defaults stay errors.
## Subsonic mathematical branch; pressure boundary conditions still select actual flow.
output MachSubsonic = AreaRatio < 1 THEN #NUM! ELSE (AreaRatio == 1 THEN 1 ELSE ROOT((mach) => AreaMach(mach) - AreaRatio, 0.01, 1, 1e-10, 200))
## Supersonic mathematical branch; this equation alone does not rule out shocks.
output MachSupersonic = AreaRatio < 1 THEN #NUM! ELSE (AreaRatio == 1 THEN 1 ELSE ROOT((mach) => AreaMach(mach) - AreaRatio, 1, 10, 1e-10, 200))
output SupersonicTemperature_K = StagnationTemperature_K / (1 + (Gamma - 1) * MachSupersonic^2 / 2)
output SupersonicPressure_Pa = StagnationPressure_Pa * (SupersonicTemperature_K / StagnationTemperature_K)^(Gamma / (Gamma - 1))
output SupersonicVelocity_m_s = MachSupersonic * SQRT(Gamma * GasConstant_J_kg_K * SupersonicTemperature_K)
output SubsonicResidual = AreaMach(MachSubsonic) - AreaRatio
output SupersonicResidual = AreaMach(MachSupersonic) - AreaRatio
output Assumptions = "Steady, quasi-one-dimensional, adiabatic reversible perfect-gas flow; constant gamma and R. No boundary layer, heat addition, shocks, chemistry, or back-pressure selection. Numerical units are K, Pa, J/(kg K), m/s."
output ScenarioGuide = "Save Baseline at area ratio 2.5: Mach values 0.2395428431 and 2.4427648455. Set AreaRatio to 1: both become 1. A write below 1 should show the input message. Restore Baseline, save, export .gridoc, and reopen."

END MODEL