"""Finite-resource, real-nitrogen steady-flow opportunity screen.

Compare one expansion with reheated two-stage expansion using the same nitrogen
batch, prepared hot-gas stream and total gas-pressure-loss target. This does not
close startup, hardware mass, materials, thermal stress or a vehicle duty cycle.
"""
import os
os.environ['OPENBLAS_NUM_THREADS']='1';os.environ['MPLCONFIGDIR']='/tmp/cryopulse-mpl'
import json,math
from pathlib import Path
from CoolProp.CoolProp import PropsSI
import numpy as np
R=Path(__file__).resolve().parents[1]
P0=101325.;TA=293.15;PP=1.02e6;P1=1e6;PF=P0+2000.;M=2.;MH=30.;CPH=1100.;RGH=287.05;TH=700.;DUTY=20.;ET=.7;EG=.9;EM=.9;EP=.7;EA=.65;CONTROL=100.
H0=PropsSI('H','P',P0,'Q',0,'Nitrogen');S0=PropsSI('S','P',P0,'Q',0,'Nitrogen');RHO0=PropsSI('D','P',P0,'Q',0,'Nitrogen');HA=PropsSI('H','P',P0,'T',TA,'Nitrogen');SA=PropsSI('S','P',P0,'T',TA,'Nitrogen');EL=H0-HA-TA*(S0-SA)
def ex(h,p=P0):return h-HA-TA*(PropsSI('S','P',p,'H',h,'Nitrogen')-SA)
def expander(p,h,out):
 si=PropsSI('S','P',p,'H',h,'Nitrogen');hs=PropsSI('H','P',out,'S',si,'Nitrogen');w=ET*(h-hs);ho=h-w;so=PropsSI('S','P',out,'H',ho,'Nitrogen');assert so>=si-1e-8;return ho,w,so-si

def ua(p,hin,hout,hotin,hotout):
 hs=np.linspace(hin,hout,161);cold=np.array([PropsSI('T','P',p,'H',h,'Nitrogen') for h in hs]);hot=hotout+(hs-hin)/(hout-hin)*(hotin-hotout);d=hot-cold
 if min(d)<=0:return None,float(min(d))
 a,b=d[:-1],d[1:];lm=np.divide(b-a,np.log(b/a),out=(a+b)/2,where=abs(b-a)>1e-9);return float(sum(M*np.diff(hs)/DUTY/lm)),float(min(d))

def cycle(tin=300.,stages=1,dp_hot=500.,pmid=330000.):
 pump=(PP-P0)/(RHO0*EP);hp=H0+pump;h1=PropsSI('H','P',P1,'T',tin,'Nitrogen');q1=M*(h1-hp);qr=0.;sgen=[]
 if stages==1:ho,w,sg=expander(P1,h1,PF);wt=M*w;sgen=[M*sg]
 else:
  hm,w,sg=expander(P1,h1,pmid);preheat=pmid-10000.;hr=PropsSI('H','P',preheat,'T',tin,'Nitrogen');qr=M*(hr-hm);ho,w2,sg2=expander(preheat,hr,PF);wt=M*(w+w2);sgen=[M*sg,M*sg2]
 q=q1+qr;C=MH*CPH;fan_shaft=MH*CPH*TH*((1+dp_hot/P0)**(RGH/CPH)-1)/EA;fan_dc=fan_shaft/EM;hotin=TH+fan_shaft/C;hotout=hotin-q/C;regen=wt*EG;pump_dc=M*pump/EM;control=CONTROL*DUTY;net=regen-pump_dc-fan_dc-control
 # Assign converter/motor loss heat to an ambient sink. Fan shaft work enters
 # the finite hot stream and is included in its outlet enthalpy, not discarded.
 rejected=wt*(1-EG)+(pump_dc-M*pump)+(fan_dc-fan_shaft)+control;enthalpy_out_minus_in=M*(ho-H0)+C*(hotout-TH);energy=enthalpy_out_minus_in-(pump_dc+fan_dc+control-regen-rejected)
 hot_ex=lambda T:C*((T-TA)-TA*math.log(T/TA));resource=EL*M+hot_ex(TH)-hot_ex(hotout);outgas=ex(ho)*M;destroyed=resource-outgas-net;total_entropy_generation=M*(PropsSI('S','P',P0,'H',ho,'Nitrogen')-S0)+C*math.log(hotout/TH)+rejected/TA;second_law_residual=destroyed-TA*total_entropy_generation;assert abs(energy)<.001 and destroyed>=0 and total_entropy_generation>=0 and abs(second_law_residual)<.001
 if stages==2:
  uar,pinchr=ua(preheat,hm,hr,hotin,hotin-qr/C);uab,pinchb=ua(P1,hp,h1,hotin-qr/C,hotout)
 else:uar=0.;pinchr=None;uab,pinchb=ua(P1,hp,h1,hotin,hotout)
 return {'physical_measurements':0,'stages':stages,'inlet_temperature_K':tin,'interstage_expander_outlet_Pa':pmid if stages==2 else None,'hot_gas_total_pressure_loss_Pa_target':dp_hot,'steady_nitrogen_kg':M,'prepared_nitrogen_inventory_kg':6.,'remaining_N2_before_unknown_startup_kg':4.,'hot_gas_kg':MH,'hot_gas_inlet_K':TH,'hot_gas_outlet_K':hotout,'heat_into_nitrogen_J':q,'turbine_shaft_J':wt,'generator_DC_J':regen,'liquid_pump_DC_J':pump_dc,'hot_gas_blower_DC_J':fan_dc,'control_DC_J':control,'net_steady_DC_J':net,'gas_exit_temperature_K':PropsSI('T','P',P0,'H',ho,'Nitrogen'),'boiler_required_UA_W_K':uab,'reheater_required_UA_W_K':uar,'minimum_boiler_pinch_K':pinchb,'minimum_reheat_pinch_K':pinchr,'resource_availability_used_J':resource,'retained_exit_N2_availability_J':outgas,'exergy_destroyed_with_ambient_loss_rejection_J':destroyed,'net_DC_over_resource_availability_used':net/resource,'turbine_entropy_generation_J_K':sum(sgen),'total_entropy_generation_J_K':total_entropy_generation,'exergy_entropy_identity_residual_J':second_law_residual,'complete_steady_first_law_residual_J':energy,'full_6kg_liquid_preparation_reversible_floor_J':6*EL,'startup_and_hardware_credit':'None. Startup, metal transients, materials and installed hardware mass are not closed; no vehicle energy credit.'}
def main():
 rows=[cycle(t,s,dp) for t in [200.,250.,300.,320.,350.] for s in [1,2] for dp in [100.,200.,500.,1000.,2000.]];one=cycle();two=cycle(stages=2);out={'physical_measurements':0,'classification':'STEADY_COMPONENT_OPPORTUNITY_NOT_CLOSED_HARDWARE_OR_VEHICLE','selected_single_stage':one,'selected_reheated_two_stage':two,'net_steady_DC_gain_J':two['net_steady_DC_J']-one['net_steady_DC_J'],'gain_fraction':two['net_steady_DC_J']/one['net_steady_DC_J']-1,'cases':rows,'assumptions':{'nitrogen_real_fluid':'CoolProp Nitrogen','turbine_isentropic_efficiency':ET,'generator_efficiency':EG,'pump_hydraulic_efficiency':EP,'motor_efficiency':EM,'hot_gas_model':'Ideal generic dry heat-source gas, Cp=1100 J/kgK, R=287.05 J/kgK; no combustion chemistry credited.','hot_gas_total_mass_kg':MH,'duty_s':DUTY,'heat_source_note':'Finite external hot stream. Engine waste heat is a future integration option, requiring measured heat, pressure loss and fuel accounting.'},'limits':['Steady-flow control volume only. Do not compare these net values directly with the P3.6 startup-inclusive 76 kJ account.','The hot gas is an explicit prepared resource; it is not ambient free heat or a new energy source.','Both layouts share the same total pressure-loss target; the extra core would have to meet that target in actual hardware.','Cold preparation, unspent liquid and cold outlet availability remain distinct.','Heat-transfer area, pressure duty, phase distribution, material gradients, startup inventory and total installed mass remain unqualified.','Driving-model energy credit is zero.']};(R/'output/reheated-cycle.json').write_text(json.dumps(out,indent=2)+'\n');print(json.dumps({k:v for k,v in out.items() if k not in ['cases','limits','assumptions']},indent=2))
if __name__=='__main__':main()