#!/usr/bin/env python
"""

  Constant-pressure, adiabatic kinetics simulation.

"""
import sys
from Cantera import *
from Cantera.Reactor import *
from Cantera.Func import *
from Cantera import rxnpath
import numpy as np
import matplotlib.pyplot as plt 
from pylab import *

# Set intial pressures and temperature of the reactor
temperature = 1200    # Kelvin
pressure    = 1.01325e5     # Pa
print "Computing auto-ignition delay for T="+str(temperature)+" K, P="+str(pressure)+" Pa" 

# specify the cti file and the gas used
rxnmech    =  'ethane_reduced.cti'                 # reaction mechanism file
mix        =  'gas'                 # gas mixture model

# specify the initial composition of your reactor
comp       =  'aC2H6:0.741, H2O:0.259'  # composition of the mixture

# Numerical parameters
npoints = 10000                            # Number of iteraions
timestep = 5.e-5                             # Timestep

# Initialize the simulation
gri3 = importPhase(rxnmech,mix)
gri3.set(T = temperature, P = pressure ,Y = comp)
r   = Reactor(gri3)
sim = ReactorNet([r])
time = 0.0
tim = zeros(npoints,'d')
tempe = zeros(npoints,'d')
press = zeros(npoints,'d')
H2 = zeros(npoints,'d')
CH4 = zeros(npoints,'d')
C2H2 = zeros(npoints,'d')
C2H4 = zeros(npoints,'d')
C2H6 = zeros(npoints,'d')
C3H6 = zeros(npoints,'d')
C4H8 = zeros(npoints,'d')
C4H6 = zeros(npoints,'d')
C4H10 = zeros(npoints,'d')
H = zeros(npoints,'d')
CH3 = zeros(npoints,'d')
C2H3 = zeros(npoints,'d')
C2H5 = zeros(npoints,'d')
C3H7 = zeros(npoints,'d')
a1C4H9 = zeros(npoints,'d')
a2C4H9 = zeros(npoints,'d')
C44 = zeros(npoints,'d')
C3H5 = zeros(npoints,'d')
H2O = zeros(npoints,'d')

# Print the state before the simulation
print r

# Start the simulation
fichier1 = open("temprature.txt", "w")
fichier2 = open("mole_fractions.txt", "w")
fichier3 = open("pressure.txt", "w")
for n in range(npoints):
    time += timestep 
    sim.advance(time)
    tim[n]   = time
    tempe[n] = r.temperature()
    press[n] = r.pressure()/1.01325e5
    H2[n]    = r.moleFraction('aH2')
    CH4[n]   = r.moleFraction('aCH4') 
    C2H2[n]  = r.moleFraction('aC2H2') 
    C2H4[n]  = r.moleFraction('aC2H4')
    C2H6[n]  = r.moleFraction('aC2H6')
    C3H6[n]  = r.moleFraction('aC3H6')
    C4H8[n]  = r.moleFraction('a1C4H8') 
    C4H6[n]  = r.moleFraction('a1,3C4H6') 
    C4H10[n] = r.moleFraction('anC4H10') 
    H[n]     = r.moleFraction('aH.') 
    CH3[n]   = r.moleFraction('aCH3.') 
    C2H3[n]  = r.moleFraction('aC2H3.') 
    C2H5[n]  = r.moleFraction('aC2H5.') 
    C3H7[n]  = r.moleFraction('a1C3H7.') 
    a1C4H9[n]= r.moleFraction('a1C4H9.') 
    a2C4H9[n]= r.moleFraction('a2C4H9.') 
    C44[n]   = r.moleFraction('a1C4-4.') 
    C3H5[n]  = r.moleFraction('aC3H5.') 
    H2O[n]   = r.moleFraction('H2O') 
    fichier1.write( '%10.3e %10.6f \n' % (sim.time(),r.temperature()))   
    fichier3.write( '%10.3e %10.6f \n' % (sim.time(),r.pressure()/1.01325e5))   
    fichier2.write( '%10.3e %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f \n' % (sim.time(),r.moleFraction('aH2'),r.moleFraction('aCH4'),r.moleFraction('aC2H2'),r.moleFraction('aC2H4'),r.moleFraction('aC2H6'),r.moleFraction('aC3H6'),r.moleFraction('a1C4H8'),r.moleFraction('a1,3C4H6'),r.moleFraction('anC4H10'),r.moleFraction('aH.'),r.moleFraction('aCH3.'),r.moleFraction('aC2H3.'),r.moleFraction('aC2H5.'),r.moleFraction('a1C3H7.'),r.moleFraction('a1C4H9.'),r.moleFraction('a2C4H9.'),r.moleFraction('a1C4-4.'),r.moleFraction('aC3H5.'),r.moleFraction('H2O')))   
fichier1.close()
fichier2.close()

# Print the state after the simulation
print r

#Calcul de la derive de la temperature 
Temperature = np.array(tempe)
Time        = np.array(tim)
H2mol       = np.array(H2)
CH4mol      = np.array(CH4)
C2H2mol     = np.array(C2H2)
C2H4mol     = np.array(C2H4)
C2H6mol     = np.array(C2H6)
C3H6mol     = np.array(C3H6)
C4H8mol     = np.array(C4H8)
C4H6mol     = np.array(C4H6)
C4H10mol    = np.array(C4H10)
Hmol        = np.array(H)
CH3mol      = np.array(CH3)
C2H3mol     = np.array(C2H3)
C2H5mol     = np.array(C2H5)
C3H7mol     = np.array(C3H7)
a1C4H9mol   = np.array(a1C4H9)
a2C4H9mol   = np.array(a2C4H9)
C44mol      = np.array(C44)
C3H5mol     = np.array(C3H5)
H2Omol      = np.array(H2O)


dTdt           = np.diff(Temperature)
dTdt_max       = dTdt.max()
dTdt_max_index = np.where(dTdt==dTdt_max)[0]
#print "Temps d'autoallumage:"+str(t_AI) 

dH2dt     = np.diff(H2mol)
dCH4dt    = np.diff(CH4mol)
dC2H2dt   = np.diff(C2H2mol)
dC2H4dt   = np.diff(C2H4mol)
dC2H6dt   = np.diff(C2H6mol)
dC3H6dt   = np.diff(C3H6mol)
dC4H8dt   = np.diff(C4H8mol)
dC4H6dt   = np.diff(C4H6mol)
dC4H10dt  = np.diff(C4H10mol)
dHdt      = np.diff(Hmol)
dCH3dt    = np.diff(CH3mol)
dC2H3dt   = np.diff(C2H3mol)
dC2H5dt   = np.diff(C2H5mol)
dC3H7dt   = np.diff(C3H7mol)
da1C4H9dt = np.diff(a1C4H9mol)
da2C4H9dt = np.diff(a2C4H9mol)
dC44dt    = np.diff(C44mol)
dC3H5dt   = np.diff(C3H5mol)
dH2Odt    = np.diff(H2Omol)

dH2dt_max=dH2dt.max()
dH2dt_max_ind = np.where(dH2dt==dH2dt_max)[0]
H2_ind = H2mol[dH2dt_max_ind]
H2_max = H2mol.max()
print ""
print "H2_max:"+str(H2_max)
print "dH2dt_max:"+str(dH2dt_max)
print "H2:"+str(H2_ind)

dCH4dt_max    = dCH4dt.max() 
dCH4dt_max_ind = np.where(dCH4dt==dCH4dt_max)[0]
CH4_ind = CH4mol[dCH4dt_max_ind]
CH4_max = CH4mol.max()
print ""
print "CH4_max:"+str(CH4_max)
print "dCH4dt_max:"+str(dCH4dt_max)
print "CH4:"+str(CH4_ind)

dC2H2dt_max   = dC2H2dt.max()
dC2H2dt_max_ind = np.where(dC2H2dt==dC2H2dt_max)[0]
C2H2_ind = C2H2mol[dC2H2dt_max_ind]
C2H2_max = C2H2mol.max()
print ""
print "C2H2_max:"+str(C2H2_max)
print "dC2H2dt_max:"+str(dC2H2dt_max)
print "C2H2:"+str(C2H2_ind)

dC2H4dt_max   = dC2H4dt.max()
dC2H4dt_max_ind = np.where(dC2H4dt==dC2H4dt_max)[0]
C2H4_ind = C2H4mol[dC2H4dt_max_ind]
C2H4_max = C2H4mol.max()
print ""
print "C2H4_max:"+str(C2H4_max)
print "dC2H4dt_max:"+str(dC2H4dt_max)
print "C2H4:"+str(C2H4_ind)

dC2H6dt_max   = dC2H6dt.max()
dC2H6dt_max_ind = np.where(dC2H6dt==dC2H6dt_max)[0]
C2H6_ind = C2H6mol[dC2H6dt_max_ind]
C2H6_max = C2H6mol.max()
print ""
print "C2H6_max:"+str(C2H6_max)
print "dC2H6dt_max:"+str(dC2H6dt_max)
print "C2H6:"+str(C2H6_ind)

dC3H6dt_max   = dC3H6dt.max()
dC3H6dt_max_ind = np.where(dC3H6dt==dC3H6dt_max)[0]
C3H6_ind = C3H6mol[dC3H6dt_max_ind]
C3H6_max = C3H6mol.max()
print ""
print "C3H6_max:"+str(C3H6_max)
print "dC3H6dt_max:"+str(dC3H6dt_max)
print "C3H6:"+str(C3H6_ind)

dC4H8dt_max   = dC4H8dt.max()
dC4H8dt_max_ind = np.where(dC4H8dt==dC4H8dt_max)[0]
C4H8_ind = C4H8mol[dC4H8dt_max_ind]
C4H8_max = C4H8mol.max()
print ""
print "C4H8_max:"+str(C4H8_max)
print "dC4H8dt_max:"+str(dC4H8dt_max)
print "C4H8:"+str(C4H8_ind)

dC4H6dt_max   = dC4H6dt.max()
dC4H6dt_max_ind = np.where(dC4H6dt==dC4H6dt_max)[0]
C4H6_ind = C4H6mol[dC4H6dt_max_ind]
C4H6_max = C4H6mol.max()
print ""
print "C4H6_max:"+str(C4H6_max)
print "dC4H6dt_max:"+str(dC4H6dt_max)
print "C4H6:"+str(C4H6_ind)

dC4H10dt_max  = dC4H10dt.max()
dC4H10dt_max_ind = np.where(dC4H10dt==dC4H10dt_max)[0]
C4H10_ind = C4H10mol[dC4H10dt_max_ind]
C4H10_max = C4H10mol.max()
print ""
print "C4H10_max:"+str(C4H10_max)
print "dC4H10dt_max:"+str(dC4H10dt_max)
print "C4H10:"+str(C4H10_ind)

dHdt_max      = dHdt.max()
dHdt_max_ind = np.where(dHdt==dHdt_max)[0]
H_ind = Hmol[dHdt_max_ind]
H_max = Hmol.max()
print ""
print "H_max:"+str(H_max)
print "dHdt_max:"+str(dHdt_max)
print "H:"+str(H_ind)

dCH3dt_max    = dCH3dt.max()
dCH3dt_max_ind = np.where(dCH3dt==dCH3dt_max)[0]
CH3_ind = CH3mol[dCH3dt_max_ind]
CH3_max = CH3mol.max()
print ""
print "CH3_max:"+str(CH3_max)
print "dCH3dt_max:"+str(dCH3dt_max)
print "CH3:"+str(CH3_ind)

dC2H3dt_max   = dC2H3dt.max() 
dC2H3dt_max_ind = np.where(dC2H3dt==dC2H3dt_max)[0]
C2H3_ind = C2H3mol[dC2H3dt_max_ind]
C2H3_max = C2H3mol.max()
print ""
print "C2H3_max:"+str(C2H3_max)
print "dC2H3dt_max:"+str(dC2H3dt_max)
print "C2H3:"+str(C2H3_ind)

dC2H5dt_max   = dC2H5dt.max()
dC2H5dt_max_ind = np.where(dC2H5dt==dC2H5dt_max)[0]
C2H5_ind = C2H5mol[dC2H5dt_max_ind]
C2H5_max = C2H5mol.max()
print ""
print "C2H5_max:"+str(C2H5_max)
print "dC2H5dt_max:"+str(dC2H5dt_max)
print "C2H5:"+str(C2H5_ind)

dC3H7dt_max   = dC3H7dt.max()
dC3H7dt_max_ind = np.where(dC3H7dt==dC3H7dt_max)[0]
C3H7_ind = C3H7mol[dC3H7dt_max_ind]
C3H7_max = C3H7mol.max()
print ""
print "C3H7_max:"+str(C3H7_max)
print "dC3H7dt_max:"+str(dC3H7dt_max)
print "C3H7:"+str(C3H7_ind)

da1C4H9dt_max = da1C4H9dt.max()
da1C4H9dt_max_ind = np.where(da1C4H9dt==da1C4H9dt_max)[0]
a1C4H9_ind = a1C4H9mol[da1C4H9dt_max_ind]
a1C4H9_max = a1C4H9mol.max()
print ""
print "a1C4H9_max:"+str(a1C4H9_max)
print "da1C4H9dt_max:"+str(da1C4H9dt_max)
print "a1C4H9:"+str(a1C4H9_ind)

da2C4H9dt_max = da2C4H9dt.max()
da2C4H9dt_max_ind = np.where(da2C4H9dt==da2C4H9dt_max)[0]
a2C4H9_ind = a2C4H9mol[da2C4H9dt_max_ind]
a2C4H9_max = a2C4H9mol.max()
print ""
print "a2C4H9_max:"+str(a2C4H9_max)
print "da2C4H9dt_max:"+str(da2C4H9dt_max)
print "a2C4H9:"+str(a2C4H9_ind)

dC44dt_max    = dC44dt.max()
dC44dt_max_ind = np.where(dC44dt==dC44dt_max)[0]
C44_ind = C44mol[dC44dt_max_ind]
C44_max = C44mol.max()
print ""
print "C44_max:"+str(C44_max)
print "dC44dt_max:"+str(dC44dt_max)
print "C44:"+str(C44_ind)

dC3H5dt_max   = dC3H5dt.max()
dC3H5dt_max_ind = np.where(dC3H5dt==dC3H5dt_max)[0]
C3H5_ind = C3H5mol[dC3H5dt_max_ind]
C3H5_max = C3H5mol.max()
print ""
print "C3H5_max:"+str(C3H5_max)
print "dC3H5dt_max:"+str(dC3H5dt_max)
print "C3H5:"+str(C3H5_ind)

dH2Odt_max    = dH2Odt.max()
dH2Odt_max_ind = np.where(dH2Odt==dH2Odt_max)[0]
H2O_ind = H2Omol[dH2Odt_max_ind]
H2O_max = H2Omol.max()
print ""
print "H2O_max:"+str(H2O_max)
print "dH2Odt_max:"+str(dH2Odt_max)
print "H2O:"+str(H2O_ind)


print """Calculation performed, output writen in xmgrace.txt"""

# Plot xmgrace.txt if option '- plot' is used
args = sys.argv
if len(args) > 1 and args[1] == '-plot':
    tmp=np.loadtxt('mole_fractions.txt')
    plt.plot(tmp[:,0],tmp[:,1],'k',label='X(H)') 
    plt.plot(tmp[:,0],tmp[:,2],'r',label='X(CH4)') 
    plt.plot(tmp[:,0],tmp[:,3],'g',label='X(H2O)') 
    plt.legend(loc='lower right',prop={'size':12})
    plt.show()