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

  EQUILIBRATE, adiabatic kinetics simulation.

"""
import sys

from Cantera import *
import numpy as np
import matplotlib.pyplot as plt

# -----------------------------------------------------------------------
#  1_ THERMO PROPERTIES
# -----------------------------------------------------------------------

# use GRI-Mech 3.0 for the methane/air mixture, and set its initial state
gas = importPhase('NAJM.cti')
nsp = gas.nSpecies()

# gaseous fuel species
fuel_species =  'CH4'
 
# air composition
air_N2_O2_ratio =  3.76

# find fuel, nitrogen, and oxygen indices
ifuel, io2, in2 = gas.speciesIndex([fuel_species, 'O2', 'N2'])

# Stoechiometric coefficient for O2
stoich_O2 = gas.nAtoms(fuel_species,'C') + 0.25*gas.nAtoms(fuel_species,'H')
#         accounts for formation of CO2     accounts for formation of H2O (4 H's per O2)

# Specify the initial composition of your reactor

print ""
print "-------------------------------------------------------- "
print "    THERMO PROPERTIES: "
print "-------------------------------------------------------- "
print ""
phi_min    = input('Enter phi min : ')
phi_min    = float(phi_min)
phi_max    = input('Enter phi max : ')
phi_max    = float(phi_max)
phi_steps  = input('Enter phi steps : ')
phi_steps  = int(phi_steps)
print ""
print "Computing equilibirum from phi_min="+str(phi_min)+" to phi_max="+str(phi_max)
print "                          for"+str(phi_steps)+"intervals"
print " "

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

# Specify intial pressures and temperature of the reactor

temperature = input('Enter temperature in kelvin : ')
temperature = float(temperature)       # Kelvin
print ""
pressure = input('Enter pressure in bar : ')
pressure = float(pressure)*1e5         # Pascal
print ""   
print "Computing equilibrium at T="+str(temperature)+" K, P="+str(pressure)+" Pa"
print " "


# Specify the equilibrium  conditions: Pression or Volume constant

print "-------------------------------------------------------- "
print "    Equilibirum conditions: "
print "-------------------------------------------------------- "
print ""
print "Eequilibrate the gas holding 2 thermo properties constant"
print "['TP','TV','HP','SP','UV']"
print ""
print " For a constant volume equilibrium, enter :      UV "
print " For a constant pressure equilibrium, enter :    HP "

cond  = raw_input('Specify the equilibrium condition : ')
cond  = str(cond)

print ""
print "Equilibrate holding" +str(cond)+"  constants"
print ""


# Output files

f_temp_phi   = open('temperature_phi-'+str(temperature)+'-'+str(pressure)+'_'+str(cond)+'.txt','w')

f_all_moles = open('all_mole_fractions_phi-'+str(temperature)+'-'+str(pressure)+'_'+str(cond)+'.csv', "w")  
writeCSV(f_all_moles,list(gas.speciesNames()))


# Initial reactor state, depending to the inputs: phi, pressure, temperature

for i in range(phi_steps):
   phi[i]   = phi_min + (phi_max - phi_min)*i/(phi_steps - 1)
   x        = zeros(nsp,'d')
   x[ifuel] = phi[i]
   x[io2]   = stoich_O2
   x[in2]   = stoich_O2*air_N2_O2_ratio
   print "phi = "+str(phi[i])
   gas.set(T = temperature, P = pressure ,X = x)


# -----------------------------------------------------------------------
#  2_ EQUILIBRIUM 
# -----------------------------------------------------------------------

   gas.equilibrate(cond, solver = 1)

# Balance record    
   tad[i]  = gas.temperature()
   names   = gas.speciesNames()

   print "At phi ="+str(phi[i])+", Tadiabatique = "+str(tad[i])+ " K"
   print " "
   f_temp_phi.write( '%10.3e %10.6f \n' % (phi[i], tad[i]))

   writeCSV(f_all_moles,list(gas.moleFractions()))    

f_temp_phi.close()
f_all_moles.close()

print 'Output written to file    all_mole_fractions_phi-'+str(temperature)+'-'+str(pressure)+'_'+str(cond)+'.csv'
print '                          temperature_phi-'+str(temperature)+'-'+str(pressure)+'_'+str(cond)+'.txt'

# Plot xmgrace.txt if option '- plot' is used
args = sys.argv
if len(args) > 1 and args[1] == '-plot':
    tmp=np.loadtxt('temperature_phi-'+str(temperature)+'-'+str(pressure)+'_'+str(cond)+'.txt')
    plt.plot(tmp[:,0],tmp[:,1],'k',label='Adiabatic flame temperature (K), '+str(cond)+' fixed')
    plt.xlabel('$\Phi$');
    plt.ylabel('Temperature (K)')
#    plt.legend(loc='upper right',prop={'size':12})
    plt.legend(bbox_to_anchor=(0.9, 0.9,1,1),loc=8)
    plt.grid()
    plt.savefig('plot_equil_TP-'+str(temperature)+'-'+str(pressure)+'_'+str(cond)+'.png', bbox_inches='tight')

import os
print 'Directory: '+os.getcwd()
print ""