###############################################################
#
# Constant pressure or Constant-volume reactor, 
#             adiabatic kinetics simulation.
#              
###############################################################

#import :

import sys

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


#################################################################
# Prepare your run
#################################################################
#Mechanisms used for the process
matrice_cas      = ['gri30','VALORANI']
type_plot = ['-','.']

	# Specify the number of time steps
nt        = 2500

	# Specify the time step
dtms	    = 1
dt = dtms * 1.e-6 #s

#Storage
nb_cas           = len(matrice_cas)
tim              = zeros(nt,'d')
nsp_cas = []
Autoignition_cas = []
FinalTemp_cas    = []
temp_cas         = []
dtemp_cas        = []
pres_cas         = []
mfrac_cas        = []

 


fuel_species     = 'CH4'
air_N2_O2_ratio = 3.76

#################
#Enter general parameters
phi        = 1.
Ti         = 1500       # Kelvin
Pi         = 1.e5       # Pascal
cond       = 'HP'     # Const P or Cons V

#Set initial conditions
for i, cas in enumerate(matrice_cas):
  print 'Cas ' + cas +' :\n'

	# Import the mechanism
  cti  = importPhase(cas + '.xml')
		#Number of species
  m                = cti.nSpecies()
  nsp_cas.append(m)
		#Find fuel, nitrogen, and oxygen indices
  ifuel, io2, in2  = cti.speciesIndex([fuel_species, 'O2', 'N2'])
		#Air composition
  stoich_O2        = cti.nAtoms(fuel_species,'C') + 0.25*cti.nAtoms(fuel_species,'H')

	#Initial mass fraction vector
  x_i              = zeros(m,'d')
  x_i[ifuel]       = phi
  x_i[io2]         = stoich_O2
  x_i[in2]         = stoich_O2*air_N2_O2_ratio

	#Set gas state
  cti.set(T = Ti, P = Pi, X = x_i)

	# create some arrays to hold the data
  temp_cas.append(zeros(nt,'d'))
  pres_cas.append(zeros(nt,'d'))
  dtemp_cas.append(zeros(nt-1,'d'))
  mfrac_cas.append(zeros([m,nt],'d'))
#################################################################
# Program starts here
#################################################################
	#Create the batch reactor
  r   = ConstPressureReactor(cti)

	#Set condition
#  if cond == 'HP':
		# Define a wall between the reactor and the environment, and
		# make it flexible, so that the pressure in the reactor is held
		# at the environment pressure.
#    env = Reservoir(air)
#    w   = Wall(r,env)
#    w.set(K = 1.0e6)   # set expansion parameter. dV/dt = KA(P_1 - P_2)
#    w.set(A = 1.0)     # set wall area

	# Now create a reactor network consisting of the single batch reactor
  sim   = ReactorNet([r])

#################
	#Run the simulation
		# Initial simulation time
  time      = 0.0


		#Loop for nt time steps of dt seconds. 
  print 'time [s] ,   T [K] ,   p [Pa] ,   u [J/kg]'
  for n in range(nt):
    time += dt
    sim.advance(time)
    tim[n]          = time
    temp_cas[i][n]  = r.temperature()
    pres_cas[i][n]  = r.pressure()
    for j in range(m):
        mfrac_cas[i][j][n] = r.massFraction(j)
    print '%10.3e %10.3f %10.3f %14.6e' % (tim[n], r.temperature(), 
                                           r.pressure(), r.intEnergy_mass())


#################################################################
# Post processing
#################################################################
	# Get autoignition timing
  Dtmax = [0,0.0]
  for n in range(nt-1):
    dtemp_cas[i][n] = (temp_cas[i][n+1]-temp_cas[i][n])/dt
    if (dtemp_cas[i][n]>Dtmax[1]):
        Dtmax[0] = n
        Dtmax[1] = dtemp_cas[i][n] 
		# Local print
  Autoignition         = (tim[Dtmax[0]]+tim[Dtmax[0] + 1])/2.
  print 'Autoignition time = ' + str(Autoignition)
		# Posterity
  Autoignition_cas.append(Autoignition)
  FinalTemp_cas.append(temp_cas[i][nt-1])


#################################################################
# Save your results if needed
#################################################################
# write output CSV file for importing into Excel
if cond == 'HP':
     csvfile = 'Phi'+ str(phi) +'_P'+str(Pi)+'_T'+str(Ti)+'_HP.csv'
elif cond == 'UV':
     csvfile = 'Phi'+ str(phi) +'_P'+str(Pi)+'_T'+str(Ti)+'_UV.csv'

f = open(csvfile,'w')
for i, cas in enumerate(matrice_cas):
	# Import the mechanism for the correct number of species
  cti  = importPhase(cas + '.xml')
  m    = cti.nSpecies()
  writeCSV(f,[cas])
  writeCSV(f,['Auto ignition time = '] + [Autoignition_cas[i]] )
  writeCSV(f,['Final Temperature = '] + [FinalTemp_cas[i]] )
  writeCSV(f,['Time','Temp','Press.'] + cti.speciesNames())
  for n in range(nt):
    f.write('%10.3e %10.3f %10.3f ' % (tim[n], temp_cas[i][n],
                                           pres_cas[i][n]))
    for j in range(m):
      f.write('%10.6e ' % (mfrac_cas[i][j][n]))
    f.write("\n")
  f.write("\n")
f.close()

print 'output written to '+csvfile


#################################################################
# Plot your results
#################################################################
# create plot
fig=figure(1)

# create first subplot 
a=fig.add_subplot(111)
for i, cas in enumerate(matrice_cas):
  if cas=='1S_CH4_MP1':
    a.plot(tim,temp_cas[i], type_plot[i], color = 'orange', label = cas)
    hold(True)
  else:
    a.plot(tim,temp_cas[i], type_plot[i], label = cas)
    hold(True)
title(r'Temperature evolution vs. Time')
xlabel(r'Time [s]', fontsize=20)
ylabel("Temperature [K]")
a.xaxis.set_major_locator(MaxNLocator(5)) # this controls the number of tick marks on the axis
hold(False)
legend(bbox_to_anchor=(0.9, 0.9,1,1),loc=0)

fig.text(0.5,0.95,r'Autoignition of $CH_{4}$ + Air mixture at $\Phi$ = 1. $T_{i}$ = 1500 K and P = 1 bar',fontsize=22,horizontalalignment='center')
grid()
subplots_adjust(left=0.08, right=0.96, wspace=0.25)

#show()
savefig('Phi'+ str(phi) +'_P'+str(Pi)+'_T'+str(Ti)+'_HP.png', bbox_inches='tight')