###############################################################
#
# ADIABATIC_FLAME - A freely-propagating, premixed flat flame 
#              
###############################################################

#import :

from Cantera import *
from Cantera.Reactor import *
from Cantera.Func import *
from matplotlib.pylab import *

#################################################################
# Prepare your run
#################################################################
#Parameter values :
	#General
pressure           =   1e5                 # pressure
temperature        =   300.0               # unburned gas temperature
phi_min            =   0.6                 # phi min
phi_max            =   2.                  # phi max
cond 	           =   'HP'                # equilibrate the gas holding 2 thermo properties constant
phi_steps          =   29                  # number of phi's to compute

# Use GRI-Mech 3.0 for the methane/air mixture, and set its initial state
#gas = importPhase('2S_CH4_BFER.cti') 
gas = GRI30()

#Print mole frac and temperature info 
	#on screen
print "******************************************************** "
print "    GRIMECH 3.0 mechanism                                "
print "    Computing Equilibirum for phi from "+str(phi_min)+" to "+str(phi_max)+", T = "+str(temperature)+" K, P = "+str(pressure)+" Pa"
print "    Equilibrate holding " +str(cond)+"  constants        "
print "******************************************************** "
print ""
	#on file
f_all_moles = open('all_mole_fractions_phi.csv', "w")  
writeCSV(f_all_moles,list(gas.speciesNames()))
f_temp_phi  = open('all_temperature_phi.txt','w')

#################
#Stoechiometry and assembling objects:

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[io2] = stoich_O2
x[in2] = stoich_O2*air_N2_O2_ratio

phi      = zeros(phi_steps,'d')
tad      = zeros(phi_steps,'d')

for i in range(phi_steps):
   phi[i]   = phi_min + (phi_max - phi_min)*i/(phi_steps - 1)
   x[ifuel] = phi[i]
   gas.set(T = temperature, P = pressure ,X = x)


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


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

   print "At phi = ", "%10.2f  "%  (phi[i])+", Tadiabatique = ", "%10.2f  "%  (tad[i])+" K"

# On file
   writeCSV(f_all_moles,list(gas.moleFractions()))
   f_temp_phi.write( '%10.3e %10.6f \n' % (phi[i], tad[i]))   

f_all_moles.close()
f_temp_phi.close()

print 'Output written to files    all_mole_fractions_phi.csv    all_temperature_phi.txt'
print ""


#################################################################
# Plot your results
#################################################################
#Plot the adiabatic flame temperature
fig=figure(1)

a=fig.add_subplot(111)
a.plot(phi,tad)
title(r'$T_{adiabatic}$ vs. $\Phi$')
xlabel(r'$\Phi$', fontsize=20)
ylabel("Adiabatic Flame Temperature [K]")
a.xaxis.set_major_locator(MaxNLocator(10)) # this controls the number of tick marks on the axis

fig.text(0.5,0.95,r'Adiabatic Flame Temperature vs equivalence ratio, GRIMECH3.0, P = 1 atm Ti = 300 K',fontsize=12,horizontalalignment='center')
grid()
show()