#################################################################""
# Puissance thermique passant dans un injecteur de 1 mm2 
# en fonction de la pression d'alimentation 
# Didier Descamps, février 2023
# usage : python3 puissance.py x > tutu.dat 
# où x=b pour le méthane B, h pour le méthane H, et p pour le propane 
#################################################################"

import sys
from math import pi, sqrt 

if(len(sys.argv)!=2):
	sys.exit("il faut entrer le gaz, b, h ou p") 
else:
	gaz=(sys.argv[1])

s=1.0e-6         #  surface de l'injecteur (m3)
pa=100000.0      #  pression atmosphérique (Pa)
ta=288.0         #  température du gaz avant injecteur (K)

if gaz=="b" :    # méthane B
    gamma=1.3	 #  cp/cv 
    r=450.0      #  cp-cv (J/kg/K) 
    PC = 40200.0 #  PCI (kJ/kg) 
elif gaz=="h" :  # méthane H
    gamma=1.3 
    r=505.0   
    PC = 49600.0   
elif gaz=="p" :  # propane
    gamma=1.14	 
    r=189.0  
    PC = 46100.0   
else : 
    sys.exit(" gaz non répertorié !")

# autres variables : 
# cp,       chaleur massique à pression constante
# p,        pression absolue d'alimentation (Pa)
# pc,       pression d'alimentation critique sub/sonique (consante)
# pchoc,    pression à l'amont du choc
# pspec,    pc/p , en sonique (constante sans dim.)
# v,        vitesse au col subsonique
# vc,       vitesse au col sonique (constante)
# tc,       température col sonique (constante)
# t,        température col subsonique
# rho,      masse volumique (kg/m3)
# qm,       débit masse (kg/s)
# P         puissance thermique (kW)

# quelques constantes (pour un gaz donné) 
cp = gamma * r / ( gamma - 1.0 ) 
pc = pa * pow((0.5*(gamma+1.0)),(gamma/(gamma-1.0)))
tc = 2.0 * ta /(gamma+1.0) 
vc = sqrt(2.0 * cp * (ta-tc)) 
pspec = pow((tc/ta),(gamma/(gamma-1.0))) 

# /*****************************
# print("  cp ", cp ) 
# print("  pc ", pc ) 
# print("  tc ", tc ) 
# print("  vc ", vc ) 
# print("pspec ", pspec ) 

p=pa 
while (p<=5.0*pa) : 
    if (p<pc) :
        t = ta * pow((pa/p),((gamma-1.0)/gamma)) 
        v= sqrt(2.0 * cp * (ta-t)) 
        rho = pa / r / t 
        qm = rho * s * v
    else :
        pchoc = p * pspec 
        rho = pchoc / r / tc 
        qm = rho * s * vc 
    P = PC * qm
    print("{0:3.3f}".format((p-pa)/pa), "  {0:3.4f}".format(P)) 
    p=p+0.002*pa