/* 
PROPANE
Débit de gaz passant dans un injecteur en fonction de
la pression d'alimentation
et Puissance thermique
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main (int argc, char *argv[])
{
const double pa=100000.0,  // pression atmosphérique (Pa)
		ta=288.0,           // température du gaz avant injecteur (K)
		gamma=1.14,	        // cp / cv
		r=189.0,            // cp - cv (J/kg/K)
		s=1.0e-6,           // m2
		PC = 46100;         // PCI (kJ/kg)
double cp,     // chaleur massique à pression constante
        p,      // pression 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 */
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))) ;

for(p=pa; p<5.01*pa; p=p+0.002*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;
	printf(" %.4lf  %.4lf   \n", 
		(p-pa)/pa, P ) ;
}

return 0 ;
}