###############################################################
#
#   HP EQUILIBRIUM COMPUTATION of the GRI 3.0 mechanism
#              
###############################################################

#import :

from Cantera import *
from Cantera.Reactor import *
from Cantera.Func import *

#################################################################
# Prepare your run
#################################################################
#Parameter values :

	#General
pressure           =   1e5                 # pressure
temperature        =   300.0               # unburned gas temperature
phi        =   1.
cond 	   =   'HP'                  # equilibrate the gas holding 2 thermo properties constant

print ""
print "    Computing Equilibirum at phi = "+str(phi)+", T = "+str(temperature)+" K, P = "+str(pressure)+" Pa"
print "    Equilibrate holding " +str(cond)+"  constants"
print ""
print ""

# Use GRI-Mech 3.0 for the methane/air mixture, and set its initial state
gas = GRI30()

#################
#Stoechiometry :

fuel_species = 'CH4'
nsp = gas.nSpecies()
stoich_O2 = gas.nAtoms(fuel_species,'C') + 0.25*gas.nAtoms(fuel_species,'H')
air_N2_O2_ratio = 3.76
ifuel, io2, in2 = gas.speciesIndex([fuel_species, 'O2', 'N2'])

x = zeros(nsp,'d')
x[ifuel] = phi
x[io2] = stoich_O2
x[in2] = stoich_O2*air_N2_O2_ratio

#################
#Assembling objects :

	#Set gas state to that of the unburned gas
gas.set(T = temperature, P = pressure, X = x)  

print "******************************************************** "
print "    Initial state:"
print "******************************************************** "
print "P =  " , "%10.4e  "%  (gas.pressure())+"    Pa"
print "T =  " , "%10.4e  "%  (gas.temperature())+"    K"
print "V =  " , "%10.4e  "%  (gas.volume_mass())+"    m3/kg"
print "U =  " , "%10.4e  "%  (gas.intEnergy_mass())+"    J/kg"
print "H =  " , "%10.4e  "%  (gas.enthalpy_mass())+"    J/kg"
print "S =  " , "%10.4e  "%  (gas.entropy_mass())+"    J/kg/K"
print ""
print ""


#################################################################
# Program starts here
#################################################################
#Equilibrium:
gas.equilibrate(cond, solver = 1)


#################################################################
# Save and print your results as needed
#################################################################
# On screen    
tad     = gas.temperature()

print ""
print "******************************************************** "
print "    Final state:"
print "    Tadiabatique = "+str(tad)+ " K"
print "******************************************************** "
print "P =  " , "%10.4e  "%  (gas.pressure())+"    Pa"
print "T =  " , "%10.4e  "%  (gas.temperature())+"    K"
print "V =  " , "%10.4e  "%  (gas.volume_mass())+"    m3/kg"
print "U =  " , "%10.4e  "%  (gas.intEnergy_mass())+"    J/kg"
print "H =  " , "%10.4e  "%  (gas.enthalpy_mass())+"    J/kg"
print "S =  " , "%10.4e  "%  (gas.entropy_mass())+"    J/kg/K"
print ""
print ""

#################
# Save the state after the simulation
f_all_moles = open('all_mole_fractions.csv', "w")   # Moles Frac at the end
writeCSV(f_all_moles,list(gas.speciesNames()))
writeCSV(f_all_moles,list(gas.moleFractions()))
f_all_moles.close()

print 'Output written to file    all_mole_fractions.csv'
print ""