import subprocess
from numarray import *
from SDToolbox import *
import re

i = 0
npoints = 10
ind_time_therm = zeros(npoints, Float)
ind_time = zeros(npoints, Float)
ind_len_therm = zeros(npoints, Float)
ind_len = zeros(npoints, Float)
exo_time_therm = zeros(npoints, Float)
exo_len_therm = zeros(npoints, Float)
cj_speed = zeros(npoints, Float)
phi  = zeros(npoints, Float)
theta_effective  = zeros(npoints, Float)
Ts = zeros(npoints, Float)


P1 = 100000 #Starting pressure in Pa

T1 = 300 #Temperature in k
#mechanism file (can change this to xml once this is run once)
mech = 'h2air_highT.cti'
xmlmech = 'h2air_highT.xml' #mechanism file (must be xml for znd)
output_frequency = 3; #Plot every nth time step, in this case 3rd
output_file_name = 'znd_phi_series'

gas1 = importPhase(mech)
gas = importPhase(mech)


while (i < npoints): 
        phi[i] = 0.7 + 0.6/npoints*(i);

        q = 'H2:' + str(2*phi[i]) + ',O2:1.,N2:3.76' # mole fractions
        print '%s: Phi = %s' % (i+1,phi[i])
        print q 
        
	# here the CJ speed is found for the above ICs
	[cj_speed[i],R2] = CJspeed(P1, T1, q, mech, 0);
	U1 = cj_speed[i];  #shock speed

	#Now make an input file and use it to call the C++ executable for        # for a znd simulation

	#Note that this simulation automatically stops when the mach # becomes greater than 0.995
 
	#Note: More relevant calucations for screen output and output files resume below this section
	###########################################################
        ##!!!! Do not change the below unless you know what you are doing!
        file = open('znd.inp', 'w')
        file.write('#Initial_Temperature_(K)\n%.4e\n' % T1);
        file.write('#Initial_Pressure_(Pa)\n%.4e\n' % P1);
        file.write('#Shock_Velocity\n%.4e\n' % U1);
        file.write('#Mole_Fractions_(string)\n%s\n'% q);
        file.write('#Mechanism_File_-_xml_format_(string)\n%s\n' % xmlmech);
        file.write('#Writing_output_frequency_flag\n%s\n' % output_frequency);
        file.write('#print_progress(1or0)\n%s\n' % 0); # use 1 to print znd progress to screen
        file.write('#Output_file_title(string)\n%s\n' % 'znd');
        file.write('#Species_Number\n%s\n' % 0);
        file.write('#Species\n%s\n' % -1);
        file.close()
        p=subprocess.Popen("./znd", stdout=subprocess.PIPE, stdin=subprocess.PIPE)
        p.stdin.write('znd.inp\n')
        p.stdin.close()
        p.wait()
 
        file = open('znd.txt','r')
        x = file.read()
        p=re.compile("Time to peak Thermicity")
        y = p.search(x).span()
        x = x[y[1]+2:];
        p= re.compile('[01234567890]+.[0-9e-]*')
        x = p.findall(x)
	
	ind_time_therm[i] = float(x[0]); 
	ind_time[i]  = float(x[1]);
	ind_len_therm[i]  = float(x[2]);
	ind_len[i]  = float(x[3]);
	exo_time_therm[i] = float(x[4]);
	exo_len_therm[i] = float(x[5]);
	file.close()
	###########################################################
	

	# Approximate effective activation energy for ZND Detonation
        gas = PostShock_fr(U1,P1,T1, q, mech)
        Ts[i] = gas.temperature()

	U1new = U1*1.02
	gas = PostShock_fr(U1*1.02,P1,T1, q, mech)
    	Ta = gas.temperature();

        ###########################################################
        ##!!!! Do not change the below unless you know what you are doing!
        file = open('znd.inp', 'w')
        file.write('#Initial_Temperature_(K)\n%.4e\n' % T1);
        file.write('#Initial_Pressure_(Pa)\n%.4e\n' % P1);
        file.write('#Shock_Velocity\n%.4e\n' % U1new);
        file.write('#Mole_Fractions_(string)\n%s\n' % q);
        file.write('#Mechanism_File_-_xml_format_(string)\n%s\n' % xmlmech);
        file.write('#Writing_output_frequency_flag\n%s\n' % output_frequency);
        file.write('#print_progress(1or0)\n%s\n' % 0); # use 1 to print znd progress to screen
        file.write('#Output_file_title(string)\n%s\n' % 'znd');
        file.write('#Species_Number\n%s\n' % 0);
        file.write('#Species\n%s\n' % -1);
        file.close()
        p=subprocess.Popen("./znd", stdout=subprocess.PIPE, stdin=subprocess.PIPE)
        p.stdin.write('znd.inp\n')
        p.stdin.close()
        p.wait()

	file = open('znd.txt','r')
        x = file.read()
        p=re.compile("Time to peak Thermicity")
        y = p.search(x).span()
        x = x[y[1]+2:];
        p= re.compile('[01234567890]+.[0-9e-]*')
        x = p.findall(x)

	ind_len_therm_a  = float(x[2]);
        file.close()
        ###########################################################


        U1new = U1*0.98
        gas = PostShock_fr(U1*0.98,P1,T1, q, mech)
	Tb = gas.temperature();

        ###########################################################
        ##!!!! Do not change the below unless you know what you are doing!
        file = open('znd.inp', 'w')
        file.write('#Initial_Temperature_(K)\n%.4e\n' % T1);
        file.write('#Initial_Pressure_(Pa)\n%.4e\n' % P1);
        file.write('#Shock_Velocity\n%.4e\n' % U1new);
        file.write('#Mole_Fractions_(string)\n%s\n' % q);
        file.write('#Mechanism_File_-_xml_format_(string)\n%s\n' % xmlmech);
        file.write('#Writing_output_frequency_flag\n%s\n' % output_frequency);
        file.write('#print_progress(1or0)\n%s\n' % 0); # use 1 to print znd progress to screen
        file.write('#Output_file_title(string)\n%s\n' % 'znd');
        file.write('#Species_Number\n%s\n' % 0);
        file.write('#Species\n%s\n' % -1);
        file.close()
        p=subprocess.Popen("./znd", stdout=subprocess.PIPE, stdin=subprocess.PIPE)
        p.stdin.write('znd.inp\n')
        p.stdin.close()
        p.wait()

        file = open('znd.txt','r')
        x = file.read()
        p=re.compile("Time to peak Thermicity")
        y = p.search(x).span()
        x = x[y[1]+2:];
        p= re.compile('[01234567890]+.[0-9e-]*')
        x = p.findall(x)
 
	ind_len_therm_b  = float(x[2]);
        file.close()
        ###########################################################

        theta_effective[i] = 1./Ts[i]*((log(ind_len_therm_a)-log(ind_len_therm_b))/((1./Ta)-(1./Tb))); # Approximate effective activation energy for ZND Detonation

	i = i + 1





print "\n\nEffective Activation Energy"

i = 0
while i < npoints:
	print "phi = %s,\n theta_effective = %s" % (phi[i], theta_effective[i])
        print "Induction time based on thermicity = %s" % (ind_time_therm[i])
        print "Induction length based on thermicity = %s" % (ind_len_therm[i])

      
        i = i+1


##################################################################################################
# CREATE OUTPUT TEXT FILE
##################################################################################################
d = datetime.date.today();
fn = 'znd_phi_series.plt';
outfile = open(fn, 'w');
outfile.write('# ZND DETONATION\n');
outfile.write('# CALCULATION RUN ON %s\n\n' % d);
outfile.write('# INITIAL CONDITIONS\n');
outfile.write('# TEMPERATURE (K) %0.2f\n' % T1);
outfile.write('# PRESSURE (Pa) %0.2f\n' % P1 );
outfile.write('# INITIAL SPECIES ARE: 2*Phi*H2 + O2 + 3.76*N2\n'); 
outfile.write('# THE OUTPUT DATA COLUMNS ARE:\n');
outfile.write('Variables = "Equivalence Ratio state 1","CJ Speed", "Induction time", "Exothermic time", "Effective Activation Energy "\n');
for i in range(npoints):
     outfile.write('%.4E \t %.4E \t %.4E \t %.4E \t %.4E \n'% (phi[i],cj_speed[i], ind_time_therm[i], exo_time_therm[i], theta_effective[i]));