"""First-principles reconstruction of NASA/TP-2000-210387 Figure 2. This is a transparent perfect-gas cycle calculation, not flight-test data. The paper's Figure 2 case does not tabulate all reconstruction closures. See methods JSON for explicit choices (ambient-matched nozzle, unit pre-ionizer and cooling-expansion pressure ratios, and calorically perfect air constants). """ from __future__ import annotations import csv, json, math from pathlib import Path OUT=Path(__file__).resolve().parent GAMMA=1.4 R_AIR=287.05 CP_AIR=1004.675 T_AMBIENT=220.0 FUEL_HEAT=45e6 ETA_COMBUSTION=1.0 DIFFUSER_PI=0.70 BURNER_PI=0.95 NOZZLE_PI=0.98 ETA_GENERATOR=0.90 ETA_ACCELERATOR=0.90 PREIONIZER_FRACTION=0.05 T_LIMIT=3000.0 def point(mach: float, extraction: float, enforce_accelerator_cap: bool=False) -> dict: """Solve one idealized NASA-reference station chain, per kg inlet air.""" if mach <= 0 or not (0 <= extraction < 1.0): raise ValueError("Mach must be positive and extraction must be a fraction below one") z=PREIONIZER_FRACTION; g=GAMMA; cp=CP_AIR t0a=T_AMBIENT*(1+(g-1)*mach**2/2) t01=t0a/(1-z*extraction) t02=t01 t03=t02*(1-extraction) if t03 >= T_LIMIT: return {"mach":mach,"extraction":extraction,"status":"burner-cap-infeasible","t0a_K":t0a,"t03_K":t03} def fuel_ratio(t04: float) -> float: den=ETA_COMBUSTION*FUEL_HEAT-cp*t04 if den <= 0 or t04 < t03: raise ValueError("No positive burner fuel solution") return cp*(t04-t03)/den def t06_for(t04: float) -> tuple[float,float]: f=fuel_ratio(t04) t06=t04 + ((1-z)*extraction/(1-z*extraction))*t0a/(1+f) return t06,f t04=T_LIMIT if enforce_accelerator_cap: hi=T_LIMIT; lo=t03 if t06_for(lo)[0] > T_LIMIT: return {"mach":mach,"extraction":extraction,"status":"accelerator-cap-infeasible","t0a_K":t0a,"t03_K":t03} for _ in range(80): mid=(lo+hi)/2 if t06_for(mid)[0] <= T_LIMIT: lo=mid else: hi=mid t04=lo t06,f=t06_for(t04) # Source equations (36), (53): component total-pressure ratios. pig=(1-extraction/ETA_GENERATOR)**(g/(g-1)) added=((1-z)*extraction/(1-z*extraction))*t0a/((1+f)*t04) pia=(1+ETA_ACCELERATOR*added)**(g/(g-1)) # Explicit reconstruction closures where Figure 2 does not publish a value. pi_preionizer=1.0 pi_cooling=1.0 pi_total=pi_preionizer*DIFFUSER_PI*pig*BURNER_PI*pi_cooling*pia*NOZZLE_PI pe_over_pa=1.0 # Explicit ambient-matched nozzle closure for this reconstruction. exit_pressure_factor=(1+(g-1)*mach**2/2)*(pi_total/pe_over_pa)**((g-1)/g) p0e_over_pa=(1+(g-1)*mach**2/2)**(g/(g-1))*pi_total if exit_pressure_factor <= 1: return {"mach":mach,"extraction":extraction,"status":"no-positive-exit-mach-at-ambient-pressure","t0a_K":t0a,"t03_K":t03,"t04_K":t04,"t06_K":t06,"fuel_air_ratio":f,"pi_total":pi_total,"p0e_over_pa":p0e_over_pa} me=math.sqrt(2/(g-1)*(exit_pressure_factor-1)) ue=me*math.sqrt(g*R_AIR*t06/(1+(g-1)*me**2/2)) ua=mach*math.sqrt(g*R_AIR*T_AMBIENT) momentum=(1+f)*ue-ua pressure=0.0 # Fully-expanded closure pe=pa, so NASA Eq. (64)'s pressure term is zero here. generator_J_per_kg_air=extraction*cp*t02 return { "mach":mach,"extraction":extraction,"status":"calculated", "t0a_K":t0a,"t01_K":t01,"t02_K":t02,"t03_K":t03,"t04_K":t04,"t05_K":t04,"t06_K":t06, "fuel_air_ratio":f,"pi_generator":pig,"pi_accelerator":pia,"pi_total":pi_total,"p0e_over_pa":p0e_over_pa,"exit_pressure_factor":exit_pressure_factor,"exit_static_pressure_over_ambient":pe_over_pa, "exit_mach":me,"exit_speed_m_s":ue,"freestream_speed_m_s":ua, "momentum_specific_thrust_N_per_kg_s":momentum, "pressure_specific_thrust_N_per_kg_s":pressure, "specific_thrust_N_per_kg_s":momentum+pressure, "generator_extraction_J_per_kg_air":generator_J_per_kg_air, "preionizer_load_J_per_kg_air":z*generator_J_per_kg_air, "accelerator_addition_J_per_kg_air":(1-z)*generator_J_per_kg_air, "accelerator_limit_exceeded":t06>T_LIMIT+1e-8, "accelerator_cap_enforced":enforce_accelerator_cap, } def digitized_reference() -> dict[float,list[tuple[float,float]]]: # Rounded visual readings; never presented as tabulated NASA model outputs. return { 0.0:[(1,.46),(1.25,.51),(1.5,.67),(1.75,.81),(2,.92),(2.25,1.01),(2.5,1.06),(2.75,1.09),(3,1.10),(3.25,1.10),(3.5,1.09),(3.75,1.06),(4,1.03),(4.25,.99),(4.5,.94),(4.75,.89),(5,.84),(5.25,.80),(5.5,.75),(5.75,.64),(6,.53),(6.25,.49),(6.5,.42),(6.75,.35),(7,.27),(7.25,.20),(7.5,.12),(7.75,.04),(7.9,0)], .25:[(1.75,.16),(2,.37),(2.25,.55),(2.5,.69),(2.75,.80),(3,.88),(3.25,.93),(3.5,.96),(3.75,.97),(4,.97),(4.25,.95),(4.5,.94),(4.75,.91),(5,.87),(5.25,.83),(5.5,.79),(5.75,.73),(6,.58),(6.25,.55),(6.5,.53),(6.75,.53),(7,.50),(7.25,.44),(7.5,.39),(7.75,.34),(8,.28),(8.25,.22),(8.5,.15),(8.75,.09),(9,.04),(9.1,0)], .5:[(2.9,0),(3,.20),(3.25,.38),(3.5,.51),(3.75,.59),(4,.66),(4.25,.72),(4.5,.75),(4.75,.77),(5,.78),(5.25,.77),(5.5,.76),(5.75,.74),(6,.71),(6.25,.71),(6.5,.69),(6.75,.66),(7,.62),(7.25,.60),(7.5,.56),(7.75,.53),(8,.49),(8.25,.45),(8.5,.41),(8.75,.38),(9,.34),(9.25,.30),(9.5,.25),(9.75,.21),(10,.17),(10.25,.13),(10.5,.10),(10.75,.05),(11,.02),(11.1,0)] } def interp(series, x): if xseries[-1][0]: return None for (x0,y0),(x1,y1) in zip(series,series[1:]): if x0<=x<=x1: return y0+(y1-y0)*(x-x0)/(x1-x0) return series[-1][1] def main(): rows=[] machs=[round(1+i*.05,2) for i in range(221)] for capped in (False,True): mode="figure2-burner-cap-case" if not capped else "dual-cap-sensitivity" for extraction in (0,.25,.5): for mach in machs: r=point(mach,extraction,capped); r["model_case"]=mode; rows.append(r) fields=sorted(set().union(*(r.keys() for r in rows))) with (OUT/"nasa-mhd-cycle-equation-reproduction.csv").open("w",newline="",encoding="utf-8") as f: w=csv.DictWriter(f,fieldnames=fields);w.writeheader();w.writerows(rows) (OUT/"nasa-mhd-cycle-model-results.json").write_text(json.dumps(rows,separators=(",",":"))+"\n",encoding="utf-8") residuals=[] refs=digitized_reference() for extraction,series in refs.items(): for mach,y in series: r=point(mach,extraction,False) if r["status"]=="calculated": residuals.append((extraction,mach,r["specific_thrust_N_per_kg_s"]/1000-y)) rms=math.sqrt(sum(e*e for _,_,e in residuals)/len(residuals)) audit={"points_compared":len(residuals),"rms_residual_kN_s_per_kg":rms,"max_abs_residual_kN_s_per_kg":max(abs(e) for _,_,e in residuals),"largest_residuals":sorted(residuals,key=lambda t:abs(t[2]),reverse=True)[:10]} audit["status"]="Historical approximate-project-reading diagnostics only; not validated residuals against NASA numerical output or total reading uncertainty." (OUT/"nasa-mhd-cycle-overlay-audit.json").write_text(json.dumps(audit,indent=2)+"\n",encoding="utf-8") svg_lines=[] colors={0.0:"#f2098c",0.25:"#168dcc",0.5:"#b27800"} W,H=1200,720; L,R,T,B=100,330,70,590 X=lambda x:L+(x-1)/(12-1)*(W-L-R) Y=lambda y:B-y/1.2*(B-T) svg_lines.extend([f'','Equation-based cycle reconstruction compared with visual readings from NASA Figure 2','For extraction ratios zero, one quarter and one half, solid curves are calculated from the simplified NASA equations. Dashed curves are a separate approximate visual digitization of the source graph. They are not measurements. All lines are conditional model results.','']) for j in range(7): y=1.2*j/6; yy=Y(y); svg_lines.append(f'{y:.1f}') for m in range(1,13): xx=X(m);svg_lines.append(f'{m}') svg_lines.extend([f'',f'',f'Flight Mach number',f'Specific thrust F / air mass flow (kN·s/kg)']) for e,series in refs.items(): col=colors[e] digit=" ".join(f'{X(m):.1f},{Y(v):.1f}' for m,v in series) svg_lines.append(f'') calc=[] for row in rows: if row['model_case']=='figure2-burner-cap-case' and row['extraction']==e and row['status']=='calculated': calc.append((row['mach'],row['specific_thrust_N_per_kg_s']/1000)) pts=" ".join(f'{X(m):.1f},{Y(v):.1f}' for m,v in calc if 0<=v<=1.2) svg_lines.append(f'') svg_lines.extend(['Generator extraction']) for idx,e in enumerate((0,.25,.5)): yy=130+idx*62; col=colors[e] svg_lines.extend([f'{e:g}: equation model',f'historical graph reading']) svg_lines.append(f'Dashed curves are approximate readings,under review; not NASA numerical data.Solid curves add stated gas/nozzleclosure assumptions; see methods.') svg_lines.append('') (OUT/"nasa-mhd-cycle-equation-comparison.svg").write_text("\n".join(svg_lines)+"\n",encoding="utf-8") checks=[] for e in (0,.25,.5): for m in (2,4,6,8): r=point(m,e,False) if r['status']=='calculated': burner_lhs=(1+r['fuel_air_ratio'])*CP_AIR*r['t04_K'] burner_rhs=CP_AIR*r['t03_K']+ETA_COMBUSTION*r['fuel_air_ratio']*FUEL_HEAT acc_lhs=(1+r['fuel_air_ratio'])*CP_AIR*(r['t06_K']-r['t05_K']) acc_rhs=(1-PREIONIZER_FRACTION)*e*CP_AIR*r['t02_K'] assert math.isclose(burner_lhs,burner_rhs,rel_tol=1e-12,abs_tol=1e-7) assert math.isclose(acc_lhs,acc_rhs,rel_tol=1e-12,abs_tol=1e-7) assert math.isclose(r['exit_static_pressure_over_ambient'],1.0,abs_tol=0) assert math.isclose(r['specific_thrust_N_per_kg_s'],(1+r['fuel_air_ratio'])*r['exit_speed_m_s']-r['freestream_speed_m_s'],rel_tol=1e-12) checks.append((m,e)) for m in (2,4,6): zero=point(m,0,False); assert zero['pi_generator']==1 and zero['pi_accelerator']==1 and zero['accelerator_addition_J_per_kg_air']==0 for m in (2,4,6,8): limited=point(m,.25,True) if limited['status']=='calculated': assert limited['t04_K']<=T_LIMIT+1e-6 and limited['t06_K']<=T_LIMIT+1e-6 figure2_rows=[r for r in rows if r['model_case']=='figure2-burner-cap-case' and r['status']=='calculated'] assert all(r['accelerator_limit_exceeded'] for r in figure2_rows if r['extraction']>0) assert all(not r['accelerator_limit_exceeded'] for r in figure2_rows if r['extraction']==0) dual_cap_rows=[r for r in rows if r['model_case']=='dual-cap-sensitivity' and r['status']=='calculated'] assert all(not r['accelerator_limit_exceeded'] for r in dual_cap_rows) point_m6=point(6.0,.25,False) plot_m6=next(y for x,y in digitized_reference()[.25] if x==6.0) pointwise_diagnostic={ "mach":6.0, "generator_enthalpy_extraction_ratio":0.25, "digitized_source_plot_specific_thrust_kN_s_per_kg":plot_m6, "estimated_visual_reading_interval_kN_s_per_kg":[plot_m6-.05,plot_m6+.05], "reconstruction_specific_thrust_kN_s_per_kg":point_m6["specific_thrust_N_per_kg_s"]/1000, "reconstruction_minus_visual_read_kN_s_per_kg":point_m6["specific_thrust_N_per_kg_s"]/1000-plot_m6, "reconstruction_pressure_thrust_N_s_per_kg":point_m6["pressure_specific_thrust_N_per_kg_s"], "reconstruction_temperature_after_accelerator_K":point_m6["t06_K"], "reported_general_accelerator_temperature_limit_K":T_LIMIT, "generator_extraction_J_per_kg_air":point_m6["generator_extraction_J_per_kg_air"], "preionizer_allocation_J_per_kg_air":point_m6["preionizer_load_J_per_kg_air"], "accelerator_addition_J_per_kg_air":point_m6["accelerator_addition_J_per_kg_air"], "energy_split_residual_J_per_kg_air":point_m6["generator_extraction_J_per_kg_air"]-point_m6["preionizer_load_J_per_kg_air"]-point_m6["accelerator_addition_J_per_kg_air"], "interpretation":"The old 0.58 reading and residual are retained as audit history only; the mismatch inference is withdrawn. Current raster observations are manually read, not NASA-tabulated data. Full uncertainty is incomplete; agreement or disagreement is not validated." } with (OUT/"figure2-direct-raster-trace.csv").open(encoding="utf8",newline="") as trace_file: trace=list(csv.DictReader(trace_file)) obs=next(r for r in trace if r["nominal_mach_for_column_selection"]=="6.0" and r["identity_status"]=="clear") actual_mach=float(obs["mach_from_pixel"]) model_at_pixel=point(actual_mach,.25,False)["specific_thrust_N_per_kg_s"]/1000 pointwise_diagnostic["legacy_comparison_status"]="withdrawn_as_evidence_of_source_mismatch" pointwise_diagnostic["current_raster_observation"]={ "mach":actual_mach,"ordinate_kN_s_per_kg":float(obs["ordinate_from_pixel_kN_s_per_kg"]), "same_column_ink_band_not_total_uncertainty":[float(obs["ordinate_band_from_pixel_min"]),float(obs["ordinate_band_from_pixel_max"])], "model_at_same_mach_kN_s_per_kg":model_at_pixel, "model_minus_reading_kN_s_per_kg":model_at_pixel-float(obs["ordinate_from_pixel_kN_s_per_kg"]), "status":"manual raster observation; total calibration and selection uncertainty incomplete; no validated agreement or source error"} assert abs(pointwise_diagnostic["energy_split_residual_J_per_kg_air"]) < 1e-6 assert pointwise_diagnostic["reconstruction_temperature_after_accelerator_K"] > T_LIMIT methods={ "validation":{"station_cases_energy_and_thrust_checks_passed":len(checks),"zero-extraction_limiting_cases_passed":3,"dual_temperature_cap_checks_passed":True,"burner_cap_mode_positive_extraction_exceeds_accelerator_cap":True,"dual_cap_mode_has_no_accelerator_cap_exceedance":True,"figure_overlay_is_comparison_only":True}, "title":"Equation-based reconstruction of NASA/TP-2000-210387 Figure 2 reference cycle", "source":{"record":"https://ntrs.nasa.gov/citations/20000095934","pdf":"https://ntrs.nasa.gov/api/citations/20000095934/downloads/20000095934.pdf","equations":"Printed pp. 5–18 (PDF pp. 11–24); Figure 2 printed p. 21 (PDF p. 27)","copyright":"NASA NTRS record says U.S. Government work; public use permitted."}, "status":"Preliminary conditional reconstruction from the simplified NASA equations; not author code, test data, or a vehicle-performance prediction. Independent reviewers checked its algebra under the stated closures. It is compared with separate approximate visual readings without fitting; the historical readings are under review and do not establish a mismatch with NASA numerical results. This is not a validated reproduction.", "figure2_reported_inputs":{"freestream_static_temperature_K":T_AMBIENT,"fuel_heating_value_MJ_per_kg":45,"combustion_efficiency":ETA_COMBUSTION,"diffuser_total_pressure_ratio":DIFFUSER_PI,"burner_total_pressure_ratio":BURNER_PI,"nozzle_total_pressure_ratio":NOZZLE_PI,"generator_isentropic_efficiency":ETA_GENERATOR,"accelerator_isentropic_efficiency":ETA_ACCELERATOR,"preionizer_power_fraction":PREIONIZER_FRACTION,"generator_enthalpy_extraction_ratios":[0,.25,.5],"burner_exit_stagnation_temperature_limit_K":T_LIMIT}, "gas_constants_added_for_reconstruction":{"gamma":GAMMA,"R_J_kg_K":R_AIR,"cp_J_kg_K":CP_AIR,"reason":"NASA Figure 2 text says it uses air properties but does not tabulate these exact reconstruction constants. They are explicit perfect-air choices here."}, "closures_added_for_this_reconstruction":{"preionizer_total_pressure_ratio":1.0,"cooling_expansion_total_pressure_ratio":1.0,"nozzle_exit_pressure":"ambient-matched (pe=pa); the source schematic supports an ambient-pressure exit as a generic cycle boundary, but applying it at every Mach point is an ideal perfectly expanded nozzle choice, not a published fixed-nozzle match. It forces the Eq. 64 pressure-thrust term to zero."}, "two_temperature_limit_modes":{"figure2-burner-cap-case":"Sets burner exit T0,4=3000 K as stated for the Figure 2 calculation; independently flags the report's additional 3000 K accelerator-exit limit. By the model equation T0,6=T0,4+positive accelerator heat addition, every successful positive-extraction case exceeds that separate limit. The no-extraction case does not. These are conditional model consequences, not tested operation.","dual-cap-sensitivity":"For sensitivity only, lowers selected burner-exit T0,4 until both report constraints T0,4<=3000 K and T0,6<=3000 K are satisfied. This is not the Figure 2 curve."}, "pointwise_diagnostic_mach6_extraction025":pointwise_diagnostic, "equation_checks":{"freestream":"T0,a=Ta[1+(gamma-1)Ma^2/2]","preionizer_and_diffuser":"T0,1=T0,2=T0,a/(1-zeta*epsilon_g)","generator":"T0,3=T0,2(1-epsilon_g); pi_g=(1-epsilon_g/eta_g)^(gamma/(gamma-1))","burner":"(1+f)cpT0,4=cpT0,3+eta_c*f*q_f","accelerator_power":"epsilon_a(1+f)cpT0,5=(1-zeta)epsilon_g*cpT0,2","accelerator_pressure":"NASA Eq. 53","nozzle_and_thrust":"NASA Eqs. 57, 62, 64 with p_e=p_a closure"}, "digitized_curve_overlay":"nasa-mhd-cycle-overlay-audit.json summarizes unfit model-versus-source-redraw residuals (95 visual readings, RMS 0.0401 and maximum absolute 0.1261 kN s/kg). These are historical differences against old project readings under review; the old allowances are not validated total uncertainty. Do not infer NASA source disagreement from these residuals. The overlay is audit history, not validation.", "source_equation_cross_reference_review":{"source_locator":"Final NASA/TP-2000-210387: Figure 2 discussion, printed p. 20 / PDF p. 26; Eqs. 36 and 53, PDF pp. 18 and 22. Earlier NASA preprint 20000032786: Eqs. 26 and 30, PDF p. 8; Figure 2 discussion, PDF p. 11.","reported_reference":"The final report's Figure 2 discussion points to Eqs. 25 and 28 for generator/accelerator pressure ratios, but the final report's equation sequence identifies 25 and 28 as diffuser relations.","equation_sequence_observed":"By equation content, final-report Eq. 36 is the generator total-pressure ratio and Eq. 53 is the accelerator total-pressure ratio. The 1999 preprint uses earlier numbering: Eq. 26 is the generator pressure-ratio relation and Eq. 30 is the accelerator pressure-ratio relation; its Figure 2 discussion explicitly points to 26 and 30 for those component ratios.","reconstruction_choice":"Use the final report's component equations by content (Eqs. 36 and 53), with the 1999 preprint used only to document the earlier numbering and cross-reference history.","status":"The numbering comparison supports a likely stale or miscited reference in the final report, but its cause is not proven. This does not validate the Figure 2 curve, the project's digitization or reconstruction, nor flight performance."}, "limitations":["No variable cp/gamma, real-gas chemistry, detailed ionization, mass capture/bleed, heat transfer or vehicle drag.","No flight-proven generator-accelerator loop; HVEPS and MAPX results are separate evidence.","Specific thrust is normalized per inlet-air mass flow and cannot be converted to total vehicle force without a matched captured-air flow and installation boundary."] } (OUT/"nasa-mhd-cycle-equation-methods.json").write_text(json.dumps(methods,indent=2,ensure_ascii=False)+"\n",encoding="utf-8") print(json.dumps({"rows":len(rows),"overlay":audit,"figure2_case_infeasible":sum(r['status']!='calculated' for r in rows if r['model_case']=='figure2-burner-cap-case'),"figure2_positive_extraction_cap_conflict_assertion_passed":True,"dual_cap_assertion_passed":True},indent=2)) if __name__=="__main__":main()