From 94149d03ea96c1c6290ed7649ca49e4b279b465b Mon Sep 17 00:00:00 2001
From: Richard West <rwest@mit.edu>
Date: Fri, 12 Jun 2009 14:43:52 -0400
Subject: [PATCH] initial commit.

hopefully I have most of the important files
---
 _rp_svg.html             |   2 +
 _rp_svg2.html            |   2 +
 canteraSimulation 2.py   | 171 +++++++++++++
 canteraSimulation.py     | 506 +++++++++++++++++++++++++++++++++++++++
 canteraSimulationHX.py   | 459 +++++++++++++++++++++++++++++++++++
 mixmaster.py             |   4 +
 prettyDot.py             |  42 ++++
 templates/rxnlist.html   |  39 +++
 templates/stylesheet.css |  57 +++++
 toMixMaster.py           | 107 +++++++++
 toMolFile.py             | 110 +++++++++
 11 files changed, 1499 insertions(+)
 create mode 100644 _rp_svg.html
 create mode 100644 _rp_svg2.html
 create mode 100644 canteraSimulation 2.py
 create mode 100644 canteraSimulation.py
 create mode 100644 canteraSimulationHX.py
 create mode 100755 mixmaster.py
 create mode 100644 prettyDot.py
 create mode 100644 templates/rxnlist.html
 create mode 100644 templates/stylesheet.css
 create mode 100644 toMixMaster.py
 create mode 100644 toMolFile.py

diff --git a/_rp_svg.html b/_rp_svg.html
new file mode 100644
index 00000000..206a3b7a
--- /dev/null
+++ b/_rp_svg.html
@@ -0,0 +1,2 @@
+<object  data="rxnpath.svg" name="rxnpath" height=100% width=100%
+type="image/svg+xml" >
diff --git a/_rp_svg2.html b/_rp_svg2.html
new file mode 100644
index 00000000..9f6a7622
--- /dev/null
+++ b/_rp_svg2.html
@@ -0,0 +1,2 @@
+<object  data="rxnpath2.dot.svg" name="rxnpath" height=200% width=200%
+type="image/svg+xml" >
diff --git a/canteraSimulation 2.py b/canteraSimulation 2.py
new file mode 100644
index 00000000..84484221
--- /dev/null
+++ b/canteraSimulation 2.py	
@@ -0,0 +1,171 @@
+# this version updated by RHW in May 2008
+"""
+Batch reactor
+
+        +------------------------+
+        |                        |
+        |                        |
+        |  reactor      (T,V)    |
+        |                        |
+        |                        |
+        +------------------------+
+
+"""
+
+######## SETTINGS
+
+saveCSVever=False  # whether to save the [T,P,massfractions] profiles as CSV files (for MixMaster)
+saveCSV=3 # saveCSV is False or an integer saying how many steps to do between saves (1=save every step)
+
+if True:
+    ctifile='chem' # the .cti extension is added later
+    clearfigs=True
+    format='-'
+
+
+####### END OF SETTINGS
+
+
+import sys
+from Cantera import *
+from Cantera.Reactor import *
+from Cantera.Func import *
+from math import *
+import random
+import pylab
+import scipy
+from pylab import *
+from scipy import *
+from scipy.optimize import leastsq
+
+import time
+def prettyTime(seconds):
+	"""returns a string of the form HH:MM:SS when given a number of seconds"""
+	return time.strftime("%H:%M:%S",time.gmtime(seconds))
+	
+
+# unit conversion factors
+mm = 0.001
+cm = 0.01
+hours = 60*60
+atm=101325
+
+# reset()  # reset Cantera (clear all cached data)
+
+
+# create a gas object
+gas = importPhase(ctifile+'.cti', 'chem')
+
+io2= gas.speciesIndex('O2(1)')
+assert io2==0, 'Was expecting O2(1) to be first species'
+maxtime=4*hours # seconds
+
+iTarget=gas.speciesIndex('C16H26O2(119)')
+
+fuelmassfracs="n-undecane(2):0.05, n-tridecane(3):0.19, SPC(4):0.11, n-hexadecane(5):0.25, SPC(6):0.12, n-nonadecane(7):0.18, n-heneicosane(8):0.10"
+gas.set(Y=fuelmassfracs)
+
+# Temperature
+temperature= 135 +273 # Kelvin
+#temperature=500
+gas.set(T=temperature)
+
+
+#The density of petroleum diesel is about 0.85 kg/l (7.09 lbs/gallon) en.wikipedia.org/wiki/Diesel
+density= 0.85*1000 # [kg/m3]
+gas.set(Rho=density)
+
+
+
+#  Oxygen concentration:
+# (doi: 10.1021/ie071481z ): At 25C and 1 atm of oxygen, O2 is sparsely 
+# soluble in octane and decane and its liquid-phase mole fractions were 
+# reported to be 2.05E-3 and 2.18E-3 respectively. Using this
+# data and the partial pressure of O2 used under the given
+# experimental conditions, the liquid-phase concentration of O2
+# was calculated by assuming Henry's Law behavior. 
+#  
+#  using van't Hoff eqn with C=1700K gives k(140C) = 4.896*k(298K)
+#  ie. concentration for same partial pressuer is ~5 times higher than at 298K
+#  so for Po2 = 0.2 atm at 140C we have mole fraction of roughly 2E-3.
+#
+# not sure what to do for diesel, so use the same!
+molefracO2=0.002
+
+
+
+def setOneMoleFrac(gas,index,molefrac):
+	molefracs=gas.moleFractions()
+	molefracs*=(1-molefrac)/(1-molefracs[index])
+	molefracs[index]=molefrac
+	gas.set(X=molefracs)
+	
+setOneMoleFrac(gas,io2,molefracO2)
+
+# create a reactor
+reactor = Reactor(gas, volume = 1.0)
+# create a network
+sim = ReactorNet([reactor])  
+           
+
+
+if saveCSV:
+    filehandle = open("MassFractionProfile.csv",'w')
+    output=['Time (s)','T','P']
+    output.extend(gas.speciesNames())
+    output.extend(['foo'])
+    writeCSV(filehandle, output)
+    
+
+if sim.time()>0.0:
+    print "warning: t>0"
+sim.setInitialTime(0.0)
+
+
+def fuelmassfraction(gas):
+	"""gives the mass fraction of unreacted fuel 
+	   (components 1-7, because O2 is 0)"""
+	fraction=0
+	for i in range(1,8):
+		fraction+=gas.massFraction(i)
+	return fraction
+	
+	
+initialFuelMassFraction=fuelmassfraction(gas)
+def conversion(gas):
+	return 1-fuelmassfraction(gas)/initialFuelMassFraction
+	
+steps=0
+
+# observed deposit formation is 8E-4 grams per gram of fuel
+finalConversion=8E-4
+
+while conversion(gas)<finalConversion and sim.time()<maxtime:
+    print "time: %.3gs (%s)\t frac used: %.2g"%(sim.time(),prettyTime(sim.time()),conversion(gas))
+    
+    sim.step(maxtime)
+    steps+=1
+
+    setOneMoleFrac(gas,io2,molefracO2)
+    
+
+    
+    if saveCSV and steps%saveCSV:
+        output=[sim.time(), gas.temperature(), gas.pressure()]
+        #output.extend(gas.massFractions())
+        output.extend([gas.massFraction(i) for i in range(gas.nSpecies())])
+        output.extend([0.0])
+        writeCSV(filehandle, output)
+
+if saveCSV and not filehandle.closed: 
+    filehandle.flush()
+    filehandle.close()     
+    
+simtime=sim.time()
+conversion=conversion(gas)
+
+print "(%d steps)\t"%steps,
+
+print "time (s)\t%.5g\t(%s)\tfrac used\t%.4g"%(simtime,prettyTime(simtime),conversion) 
+
+
diff --git a/canteraSimulation.py b/canteraSimulation.py
new file mode 100644
index 00000000..37021a9c
--- /dev/null
+++ b/canteraSimulation.py
@@ -0,0 +1,506 @@
+# this version updated by RHW in May 2008
+"""
+Batch reactor
+
+        +------------------------+
+        |                        |
+        |                        |
+        |  reactor      (T,V)    |
+        |                        |
+        |                        |
+        +------------------------+
+
+"""
+
+######## SETTINGS
+
+saveCSVever=False  # whether to save the [T,P,massfractions] profiles as CSV files (for MixMaster)
+saveCSV=3 # saveCSV is False or an integer saying how many steps to do between saves (1=save every step)
+
+if True:
+    ctifile='chem' # the .cti extension is added later
+    clearfigs=True
+    format='-'
+
+
+molefracPeroxide=1e-6
+
+
+####### END OF SETTINGS
+
+
+import sys
+from Cantera import *
+from Cantera.Reactor import *
+from Cantera.Func import *
+from math import *
+import random
+import pylab
+import scipy
+from pylab import *
+from scipy import *
+from scipy.optimize import leastsq
+
+import time
+def prettyTime(seconds):
+	"""returns a string of the form HH:MM:SS when given a number of seconds"""
+	return time.strftime("%H:%M:%S",time.gmtime(seconds))
+	
+
+# unit conversion factors
+mm = 0.001
+cm = 0.01
+hours = 60*60
+atm=101325
+
+reset()  # reset Cantera (clear all cached data)
+
+
+
+for molefracPeroxide in [10**i for i in range(-9,-1)]:
+	# create a gas object
+	gas = importPhase(ctifile+'.cti', 'chem')
+	
+	io2= gas.speciesIndex('O2(1)')
+	assert io2==0, 'Was expecting O2(1) to be first species'
+	maxtime=4*hours # seconds
+	
+	iPeroxide=gas.speciesIndex('peroxydcbz(9)')
+	
+	fuelmassfracs="n-undecane(2):0.05, n-tridecane(3):0.19, SPC(4):0.11, n-hexadecane(5):0.25, SPC(6):0.12, n-nonadecane(7):0.18, n-heneicosane(8):0.10"
+	gas.set(Y=fuelmassfracs)
+	
+	# Temperature
+	temperature= 135 +273 # Kelvin
+	#temperature=500
+	gas.set(T=temperature)
+	
+	
+	#The density of petroleum diesel is about 0.85 kg/l (7.09 lbs/gallon) en.wikipedia.org/wiki/Diesel
+	density= 0.85*1000 # [kg/m3]
+	gas.set(Rho=density)
+	
+	
+	
+	#  Oxygen concentration:
+	# (doi: 10.1021/ie071481z ): At 25C and 1 atm of oxygen, O2 is sparsely 
+	# soluble in octane and decane and its liquid-phase mole fractions were 
+	# reported to be 2.05E-3 and 2.18E-3 respectively. Using this
+	# data and the partial pressure of O2 used under the given
+	# experimental conditions, the liquid-phase concentration of O2
+	# was calculated by assuming Henry's Law behavior. 
+	#  
+	#  using van't Hoff eqn with C=1700K gives k(140C) = 4.896*k(298K)
+	#  ie. concentration for same partial pressuer is ~5 times higher than at 298K
+	#  so for Po2 = 0.2 atm at 140C we have mole fraction of roughly 2E-3.
+	#
+	# not sure what to do for diesel, so use the same!
+	molefracO2=0.002
+	
+	
+	
+	def setOneMoleFrac(gas,index,molefrac):
+		molefracs=gas.moleFractions()
+		molefracs*=(1-molefrac)/(1-molefracs[index])
+		molefracs[index]=molefrac
+		gas.set(X=molefracs)
+		
+	setOneMoleFrac(gas,io2,molefracO2)
+	
+	
+	
+	setOneMoleFrac(gas,iPeroxide,molefracPeroxide)
+#	print "initial mole fraction of Peroxide:",molefracPeroxide
+	
+	# create a reactor
+	reactor = Reactor(gas, volume = 1.0)
+	# create a network
+	sim = ReactorNet([reactor])  
+	           
+	
+	
+	if saveCSV:
+	    filehandle = open("MassFractionProfile.csv",'w')
+	    output=['Time (s)','T','P']
+	    output.extend(gas.speciesNames())
+	    output.extend(['foo'])
+	    writeCSV(filehandle, output)
+	    
+	
+	if sim.time()>0.0:
+	    print "warning: t>0"
+	sim.setInitialTime(0.0)
+	
+	
+	def fuelmassfraction(gas):
+		"""gives the mass fraction of unreacted fuel 
+		   (components 1-7, because O2 is 0)"""
+		fraction=0
+		for i in range(1,8):
+			fraction+=gas.massFraction(i)
+		return fraction
+		
+		
+	initialFuelMassFraction=fuelmassfraction(gas)
+	def conversion(gas):
+		return 1-fuelmassfraction(gas)/initialFuelMassFraction
+		
+	steps=0
+	
+	# observed deposit formation is 8E-4 grams per gram of fuel
+	finalConversion=8E-4
+	
+	while conversion(gas)<finalConversion and sim.time()<maxtime:
+	    #print "time: %.3gs (%s)\t frac used: %.2g"%(sim.time(),prettyTime(sim.time()),conversion(gas))
+	    
+	    sim.step(maxtime)
+	    steps+=1
+	
+	    setOneMoleFrac(gas,io2,molefracO2)
+	    
+	
+	    
+	    if saveCSV and steps%saveCSV:
+	        output=[sim.time(), gas.temperature(), gas.pressure()]
+	        #output.extend(gas.massFractions())
+	        output.extend([gas.massFraction(i) for i in range(gas.nSpecies())])
+	        output.extend([0.0])
+	        writeCSV(filehandle, output)
+	
+	if saveCSV and not filehandle.closed: 
+	    filehandle.flush()
+	    filehandle.close()     
+	    
+	simtime=sim.time()
+	conversion=conversion(gas)
+	
+	print "(%d steps)\t"%steps,
+	#print "time (s)\t%.5g\t(%s)\tfrac used\t%.4g"%(simtime,prettyTime(simtime),conversion) 
+	print "molefracPeroxide\t%g\ttime(s)\t%g"%(molefracPeroxide,simtime)
+
+
+
+assert False, "stop"
+
+# temperature=1500
+# pressure=1.0e5
+
+f=figure(1)
+title("Effect of pressure")
+Ts=range(901,2500,50)
+Ps=[1e5*10**(i) for i in arange(-1,2.1,1)] # was ,0.25)]
+Ps.insert(0,0) # Pratsinis
+Ps.insert(0,-1) # Raghavan
+if clearfigs: clf()
+meshTemps=[]
+meshPressures=[]
+meshLogKs=[]
+for P in Ps:
+    temps=[]
+    logKs=[]
+    
+    if P==-1:
+        if not clearfigs: continue
+        label="Raghavan (2001)"
+        for T in [646, 2044]:
+            temps.append(T)
+            logKs.append(RaghavanLogK(T))    
+    elif P==0:
+        if not clearfigs: continue
+        label="Pratsinis (1990)"
+        for T in [973, 1273]:
+            temps.append(T)
+            logKs.append(PratsinisLogK(T))
+    else:
+        for T in Ts:
+            label="%.2g bar"%(P/1e5)
+            X="TiCl4:0.01, O2:0.05, Ar:0.94"
+            s=Simulation(T,P,X)
+            print "T (K)\t%d\tP (bar)\t%.2g\t "%(T,P/1e5),
+            logK=s.getLogK()
+            if logK>-Inf:
+                temps.append(T)
+                logKs.append(logK)
+    
+        meshTemps.append(temps)
+        meshPressures.append([P for T in temps]) 
+        meshLogKs.append(logKs)
+    temps=array(temps)
+    logKs=array(logKs)
+    plot( 1000.0/temps, logKs, format, label=label, linewidth=2 )
+    
+ylabel("ln k")
+xlabel("1000/T")
+if clearfigs:
+    legend(loc=0)
+    drawPratsinisBox()
+    drawRaghavanBox()
+
+
+### 3D plot
+import matplotlib.axes3d as p3
+f=figure(2)
+ax=p3.Axes3D(f)
+show()
+
+
+##RAGHAVAN'S
+#temps=array([  1046.82525236,  1073.85699043,
+#        1151.80924153,  1257.77744668,  1298.76346761,  1406.39240553,
+#        1439.28607796,  1514.65377135,  1604.81167741,  1645.52873364,
+#        1715.11423198,  1740.23047543,  1806.84217932,  2043.4652303 ])
+#x=1000.0/temps
+#y=log(array([   551205.35156326,
+#         429845.99298875,   736792.54475485,   964331.41922609,
+#         887997.53785585,  1198889.6613466 ,   942090.5684987 ,
+#        1436794.95733164,  1706284.57763703,  1810558.77650401,
+#        1458560.42080549,  1526135.31697436,  1695978.55319514,
+#        2384852.8012917 ]))/log(10.0)
+## z given by Raghavan's arrhenius rate       
+#z=array([RaghavanLogK(T) for T in temps])
+
+
+        
+# z calculated from Raghavan's measurement (assuming lorentzian, with his Eact)
+def RaghavanLogKexpt(Tp,NNo,tc):
+    Eact=108681 # J/mol
+   # Eact= 88800 # J/mol (pratsinis)
+    molargasconstant=8.314472 # J/mol/K
+    # NNo=NNo-0.1;
+    f=0.2
+    if NNo<f: return 15.0
+    if NNo>=1.0: return -15.0
+    numerator= -1.0*( log((NNo-f)/(1-f)) )
+    denominator= tc * sqrt( pi * molargasconstant * Tp / Eact)
+    #print "Tp=%g\tNNo=%g\tTc=%g"%(Tp,NNo,tc)
+    k=numerator/denominator
+    #print "k=",numerator,'/',denominator,"\t",k
+    #print "log(k)=",log(k)
+    return log(k)
+    
+temps=array([  646.,   757.,  1047.,  1074.,  1152.,  1258.,  1299.,  1407.,
+        1440.,  1516.,  1607.,  1648.,  1718.,  1743.,  1809.,  2044.])       
+measuredNNos=array([ 1.  ,  0.97,  1.02,  0.9 ,  0.96,  1.04,  0.94,  0.49,  0.44,
+        0.6 ,  0.33,  0.35,  0.16,  0.24,  0.17,  0.17])
+reactiontimes=array([ 0.467,  0.166,  0.058,  0.063,  0.04 ,  0.032,  0.032,  0.027,
+        0.023,  0.024,  0.019,  0.019,  0.018,  0.018,  0.013,  0.01 ])
+pressures=array([  165163.21252451,   252971.75364402,   551205.35156326,
+         429845.99298875,   736792.54475485,   964331.41922609,
+         887997.53785585,  1198889.6613466 ,   942090.5684987 ,
+        1436794.95733164,  1706284.57763703,  1810558.77650401,
+        1458560.42080549,  1526135.31697436,  1695978.55319514,
+        2384852.8012917 ])
+x=1000.0/temps 
+y=log(pressures)/log(10.0)
+z=array([RaghavanLogKexpt( temps[i], measuredNNos[i], reactiontimes[i]) for i in range(temps.size)])
+
+RaghavanMeasuredInverseTs=x
+RaghavanMeasuredPressures=y
+RaghavanMeasuredLogKs=z
+
+ax.plot3D(x,y,z,'o-',label="Raghavan (2001)", linewidth=2)
+save('lineXrag.out',x)
+save('lineYrag.out',y)
+save('lineZrag.out',z)
+
+temps=arange(973.0,1274,50)
+x=1000.0/temps
+y=log(array([OneAtm for T in temps]))/log(10.0)
+z=array([PratsinisLogK(T) for T in temps])
+ax.plot3D(x,y,z,label="Pratsinis (1990)", linewidth=2)
+
+PratData=array([[ 0.883, 1.1112 ],[ 0.8981, 1.0121 ],[ 0.9149, 0.9306 ],[ 0.9317, 0.6962 ],[ 0.9499, 0.4882 ],[ 0.9677, 0.1884 ],[ 0.9832, -0.2687 ],[ 1.0073, -0.2841 ]])
+x=PratData[:,0]
+y=log(array([OneAtm for t in x]))/log(10.0)
+z=PratData[:,1]
+ax.plot3D(x,y,z,'o-',label="Pratsinis (1990) 1:5", linewidth=2)
+save('lineXprat.out',x)
+save('lineYprat.out',y)
+save('lineZprat.out',z)
+
+PratsinisMeasuredInverseTs=x
+PratsinisMeasuredPressures=y
+PratsinisMeasuredLogKs=z
+
+
+def plotRaghavanMeasurements():
+    temps=array([  646.,   757.,  1047.,  1074.,  1152.,  1258.,  1299.,  1407.,
+            1440.,  1516.,  1607.,  1648.,  1718.,  1743.,  1809.,  2044.])       
+    measuredNNos=array([ 1.  ,  0.97,  1.02,  0.9 ,  0.96,  1.04,  0.94,  0.49,  0.44,
+            0.6 ,  0.33,  0.35,  0.16,  0.24,  0.17,  0.17])
+    reactiontimes=array([ 0.467,  0.166,  0.058,  0.063,  0.04 ,  0.032,  0.032,  0.027,
+            0.023,  0.024,  0.019,  0.019,  0.018,  0.018,  0.013,  0.01 ])
+    pressures=array([  165163.21252451,   252971.75364402,   551205.35156326,
+             429845.99298875,   736792.54475485,   964331.41922609,
+             887997.53785585,  1198889.6613466 ,   942090.5684987 ,
+            1436794.95733164,  1706284.57763703,  1810558.77650401,
+            1458560.42080549,  1526135.31697436,  1695978.55319514,
+            2384852.8012917 ])
+            
+    x=1000.0/temps 
+    for i in range(temps.size):
+        temp=temps[i]
+        measuredNNo=measuredNNos[i]
+        pressure=pressures[i]
+        reactiontime=reactiontimes[i]
+        errorInNNo=0.1
+        x=1000.0/temp
+        y=RaghavanLogKexpt( temp, measuredNNo, reactiontime )
+        ymin =RaghavanLogKexpt( temp, measuredNNo+errorInNNo, reactiontime )
+        ymax =RaghavanLogKexpt( temp, measuredNNo-errorInNNo, reactiontime )
+        ybar=array([ymin,ymax])
+        xbar=array([x, x])
+        if y>9:
+            plot( [x], [y], 'b^', label="Raghavan", linewidth=1 )
+        elif y<-5:
+            plot( [x], [y], 'bv', label="Raghavan", linewidth=1 )
+        else:
+            plot( [x], [y], 'ro', label="Raghavan", linewidth=1 )
+        
+        plot( xbar, ybar, 'b-_', label="Raghavan", linewidth=1 )
+        print "Tp\t%dK\tNNo\t%g\tt\t%g\tk_eff\t%.4g\tmin\t%.4g\tmax\t%.4g"%(temp,measuredNNo, reactiontime, y,ymin,ymax)
+    
+def plotPratsinisMeasurements():
+    plot( PratsinisMeasuredInverseTs, PratsinisMeasuredLogKs, 'go', label="Pratsinis", linewidth=2 )
+
+
+X=1000.0/array(meshTemps)
+Y=log(array(meshPressures))/log(10.0)
+Z=array(meshLogKs)
+surf1=ax.plot_surface(X,Y,Z)
+frame1=ax.plot_wireframe(X,Y,Z)
+
+Zprat=array([[PratsinisLogK(1000.0/i) for i in j] for j in X])
+surf2=ax.plot_surface(X,Y,Zprat)
+frame2=ax.plot_wireframe(X,Y,Zprat)
+frame1._wrapped.set_color('r')
+
+Zrag=array([[RaghavanLogK(1000.0/i) for i in j] for j in X])
+surf3=ax.plot_surface(X,Y,Zrag)
+frame3=ax.plot_wireframe(X,Y,Zrag)
+frame2._wrapped.set_color('g')
+
+save('meshX.out', X)
+save('meshY.out', Y)
+save('meshZme.out',Z)
+save('meshZprat.out',Zprat)
+save('meshZrag.out',Zrag)
+
+
+for surf in [surf1, surf2, surf3]:
+    surf._wrapped.set_alpha(0.2)
+    surf._wrapped.set_linewidth(0)
+    
+for frame in [frame1, frame2, frame3]:
+    frame._wrapped.set_alpha(0.5)
+    
+ax.set_xlabel('1000/T')
+ax.set_ylabel('log_{10} P')
+ax.set_zlabel('ln k')
+
+
+## now add things
+figure(1)
+if clearfigs:
+    plotRaghavanMeasurements()
+    plot( PratsinisMeasuredInverseTs, PratsinisMeasuredLogKs, 'go', label="Pratsinis", linewidth=2 )
+
+#### O2 
+figure(3)
+if clearfigs: clf()
+P=1.0e5
+title("Effect of excess O2 (no argon)")
+o2xs=[-1,0,1,4,9]
+for o2x in o2xs:
+    temps=[]
+    logKs=[]
+    if o2x==-1:
+        if not clearfigs: continue
+        label="Raghavan (2001)"
+        for T in [646, 2044]:
+            temps.append(T)
+            logKs.append(RaghavanLogK(T)) 
+    elif o2x==0:
+        if not clearfigs: continue
+        label="Pratsinis (1990)"
+        for T in [973, 1273]:
+            temps.append(T)
+            logKs.append(PratsinisLogK(T))
+    else:
+        for T in Ts:
+            #label="O2/TiCl4=%.2g"%(o2x)
+            ar=0.0
+            ticl=(1.0-ar)/(1+o2x)
+            o2=ticl*o2x
+            X="TiCl4:%f, O2:%f, Ar:%f"%(ticl,o2,ar)
+            label=X
+            s=Simulation(T,P,X)
+            print "T\t%dK\tO2/TiCl4\t%.2g\t"%(T,o2x),
+            logK=s.getLogK()
+            if logK>-Inf:
+                temps.append(T)
+                logKs.append(logK)
+    temps=array(temps)
+    logKs=array(logKs)
+    plot( 1000.0/temps, logKs, format, label=label, linewidth=2 )
+    
+
+
+ylabel("ln k")
+xlabel("1000/T")
+if clearfigs:
+    legend(loc=0)
+    drawPratsinisBox()
+    drawRaghavanBox()
+    plotRaghavanMeasurements()
+    plot( PratsinisMeasuredInverseTs, PratsinisMeasuredLogKs, 'go', label="Pratsinis", linewidth=2 )
+
+figure(4)
+if clearfigs: clf()
+title("Effect of argon fraction (4:1 O2:TiCl4)")
+P=1.0e5
+argons=[-2,-1,0,0.9,0.99]
+o2x=4
+for argon in argons:
+    temps=[]
+    logKs=[]
+    if argon==-2:
+        if not clearfigs: continue
+        label="Raghavan (2001)"
+        for T in [646, 2044]:
+            temps.append(T)
+            logKs.append(RaghavanLogK(T)) 
+    elif argon==-1:
+        if not clearfigs: continue
+        label="Pratsinis (1990)"
+        for T in [973, 1273]:
+            temps.append(T)
+            logKs.append(PratsinisLogK(T))
+    else:
+        for T in Ts:
+            #label="O2/TiCl4=%.2g"%(o2x)
+            ar=argon
+            ticl=(1.0-ar)/(1+o2x)
+            o2=ticl*o2x
+            X="TiCl4:%f, O2:%f, Ar:%f"%(ticl,o2,ar)
+            label=X
+            s=Simulation(T,P,X)
+            print "T (K)\t%d\tArgon\t%.2g\t"%(T,argon),
+            logK=s.getLogK()
+            if logK>-Inf:
+                temps.append(T)
+                logKs.append(logK)
+    temps=array(temps)
+    logKs=array(logKs)
+    plot( 1000.0/temps, logKs, format, label=label, linewidth=2 )
+
+
+
+ylabel("ln k")
+xlabel("1000/T")
+if clearfigs:
+    legend(loc=0)
+    drawPratsinisBox()
+    drawRaghavanBox()
+    plotRaghavanMeasurements()
+    plot( PratsinisMeasuredInverseTs, PratsinisMeasuredLogKs, 'go', label="Pratsinis", linewidth=2 )
diff --git a/canteraSimulationHX.py b/canteraSimulationHX.py
new file mode 100644
index 00000000..833c728f
--- /dev/null
+++ b/canteraSimulationHX.py
@@ -0,0 +1,459 @@
+# this version updated by RHW in May 2008
+"""
+Batch reactor
+
+        +------------------------+
+        |                        |
+        |                        |
+        |  reactor      (T,V)    |
+        |                        |
+        |                        |
+        +------------------------+
+
+"""
+
+######## SETTINGS
+
+saveCSVever=False  # whether to save the [T,P,massfractions] profiles as CSV files (for MixMaster)
+saveCSV=3 # saveCSV is False or an integer saying how many steps to do between saves (1=save every step)
+
+if True:
+    ctifile='chem' # the .cti extension is added later
+    clearfigs=True
+    format='-'
+
+
+####### END OF SETTINGS
+
+
+import sys
+from Cantera import *
+from Cantera.Reactor import *
+from Cantera.Func import *
+from math import *
+import random
+import pylab
+import scipy
+from pylab import *
+from scipy import *
+from scipy.optimize import leastsq
+
+import time
+def prettyTime(seconds):
+	"""returns a string of the form HH:MM:SS when given a number of seconds"""
+	return time.strftime("%H:%M:%S",time.gmtime(seconds))
+	
+
+# unit conversion factors
+mm = 0.001
+cm = 0.01
+hours = 60*60
+atm=101325
+
+# reset()  # reset Cantera (clear all cached data)
+
+
+# create a gas object
+gas = importPhase(ctifile+'.cti', 'chem')
+
+maxtime=4*hours # seconds
+
+iTarget=gas.speciesIndex('HXD13(1)')
+
+molefracs="HXD13(1):6.829e-4, CH4(2):0.104, H2(3):1.56e-2, N2:0.8797"
+gas.set(X=molefracs)
+
+
+# Temperature
+temperature= 1350# Kelvin
+
+gas.set(T=temperature)
+gas.set(P=1*atm)
+
+
+# create a reactor
+reactor = Reactor(gas, volume = 1.0)
+# create a network
+sim = ReactorNet([reactor])  
+           
+if saveCSV:
+    filehandle = open("MassFractionProfile.csv",'w')
+    output=['Time (s)','T','P']
+    output.extend(gas.speciesNames())
+    output.extend(['foo'])
+    writeCSV(filehandle, output)
+    
+
+if sim.time()>0.0:
+    print "warning: t>0"
+sim.setInitialTime(0.0)
+
+
+def conversion(gas):
+	"""gives the conversion of iTarget"""
+	fraction=0
+	for i in range(1,8):
+		fraction+=gas.massFraction(i)
+	return fraction
+	
+	
+initialMoleFraction=gas.moleFraction(iTarget)
+def conversion(gas):
+	""" gives the conversion of iTarget"""
+	return 1-gas.moleFraction(iTarget)/initialMoleFraction
+	
+steps=0
+
+
+finalConversion=0.9
+
+while conversion(gas)<finalConversion and sim.time()<maxtime:
+    print "time: %.3gs (%s)\t frac used: %.2g"%(sim.time(),prettyTime(sim.time()),conversion(gas))
+    
+    sim.step(maxtime)
+    steps+=1
+    
+    if saveCSV and steps%saveCSV:
+        output=[sim.time(), gas.temperature(), gas.pressure()]
+        #output.extend(gas.massFractions())
+        output.extend([gas.massFraction(i) for i in range(gas.nSpecies())])
+        output.extend([0.0])
+        writeCSV(filehandle, output)
+
+if saveCSV and not filehandle.closed: 
+    filehandle.flush()
+    filehandle.close()     
+    
+simtime=sim.time()
+conversion=conversion(gas)
+
+print "(%d steps)\t"%steps,
+
+print "time (s)\t%.5g\t(%s)\tfrac used\t%.4g"%(simtime,prettyTime(simtime),conversion) 
+
+
+
+assert False, "stop"
+
+# temperature=1500
+# pressure=1.0e5
+
+f=figure(1)
+title("Effect of pressure")
+Ts=range(901,2500,50)
+Ps=[1e5*10**(i) for i in arange(-1,2.1,1)] # was ,0.25)]
+Ps.insert(0,0) # Pratsinis
+Ps.insert(0,-1) # Raghavan
+if clearfigs: clf()
+meshTemps=[]
+meshPressures=[]
+meshLogKs=[]
+for P in Ps:
+    temps=[]
+    logKs=[]
+    
+    if P==-1:
+        if not clearfigs: continue
+        label="Raghavan (2001)"
+        for T in [646, 2044]:
+            temps.append(T)
+            logKs.append(RaghavanLogK(T))    
+    elif P==0:
+        if not clearfigs: continue
+        label="Pratsinis (1990)"
+        for T in [973, 1273]:
+            temps.append(T)
+            logKs.append(PratsinisLogK(T))
+    else:
+        for T in Ts:
+            label="%.2g bar"%(P/1e5)
+            X="TiCl4:0.01, O2:0.05, Ar:0.94"
+            s=Simulation(T,P,X)
+            print "T (K)\t%d\tP (bar)\t%.2g\t "%(T,P/1e5),
+            logK=s.getLogK()
+            if logK>-Inf:
+                temps.append(T)
+                logKs.append(logK)
+    
+        meshTemps.append(temps)
+        meshPressures.append([P for T in temps]) 
+        meshLogKs.append(logKs)
+    temps=array(temps)
+    logKs=array(logKs)
+    plot( 1000.0/temps, logKs, format, label=label, linewidth=2 )
+    
+ylabel("ln k")
+xlabel("1000/T")
+if clearfigs:
+    legend(loc=0)
+    drawPratsinisBox()
+    drawRaghavanBox()
+
+
+### 3D plot
+import matplotlib.axes3d as p3
+f=figure(2)
+ax=p3.Axes3D(f)
+show()
+
+
+##RAGHAVAN'S
+#temps=array([  1046.82525236,  1073.85699043,
+#        1151.80924153,  1257.77744668,  1298.76346761,  1406.39240553,
+#        1439.28607796,  1514.65377135,  1604.81167741,  1645.52873364,
+#        1715.11423198,  1740.23047543,  1806.84217932,  2043.4652303 ])
+#x=1000.0/temps
+#y=log(array([   551205.35156326,
+#         429845.99298875,   736792.54475485,   964331.41922609,
+#         887997.53785585,  1198889.6613466 ,   942090.5684987 ,
+#        1436794.95733164,  1706284.57763703,  1810558.77650401,
+#        1458560.42080549,  1526135.31697436,  1695978.55319514,
+#        2384852.8012917 ]))/log(10.0)
+## z given by Raghavan's arrhenius rate       
+#z=array([RaghavanLogK(T) for T in temps])
+
+
+        
+# z calculated from Raghavan's measurement (assuming lorentzian, with his Eact)
+def RaghavanLogKexpt(Tp,NNo,tc):
+    Eact=108681 # J/mol
+   # Eact= 88800 # J/mol (pratsinis)
+    molargasconstant=8.314472 # J/mol/K
+    # NNo=NNo-0.1;
+    f=0.2
+    if NNo<f: return 15.0
+    if NNo>=1.0: return -15.0
+    numerator= -1.0*( log((NNo-f)/(1-f)) )
+    denominator= tc * sqrt( pi * molargasconstant * Tp / Eact)
+    #print "Tp=%g\tNNo=%g\tTc=%g"%(Tp,NNo,tc)
+    k=numerator/denominator
+    #print "k=",numerator,'/',denominator,"\t",k
+    #print "log(k)=",log(k)
+    return log(k)
+    
+temps=array([  646.,   757.,  1047.,  1074.,  1152.,  1258.,  1299.,  1407.,
+        1440.,  1516.,  1607.,  1648.,  1718.,  1743.,  1809.,  2044.])       
+measuredNNos=array([ 1.  ,  0.97,  1.02,  0.9 ,  0.96,  1.04,  0.94,  0.49,  0.44,
+        0.6 ,  0.33,  0.35,  0.16,  0.24,  0.17,  0.17])
+reactiontimes=array([ 0.467,  0.166,  0.058,  0.063,  0.04 ,  0.032,  0.032,  0.027,
+        0.023,  0.024,  0.019,  0.019,  0.018,  0.018,  0.013,  0.01 ])
+pressures=array([  165163.21252451,   252971.75364402,   551205.35156326,
+         429845.99298875,   736792.54475485,   964331.41922609,
+         887997.53785585,  1198889.6613466 ,   942090.5684987 ,
+        1436794.95733164,  1706284.57763703,  1810558.77650401,
+        1458560.42080549,  1526135.31697436,  1695978.55319514,
+        2384852.8012917 ])
+x=1000.0/temps 
+y=log(pressures)/log(10.0)
+z=array([RaghavanLogKexpt( temps[i], measuredNNos[i], reactiontimes[i]) for i in range(temps.size)])
+
+RaghavanMeasuredInverseTs=x
+RaghavanMeasuredPressures=y
+RaghavanMeasuredLogKs=z
+
+ax.plot3D(x,y,z,'o-',label="Raghavan (2001)", linewidth=2)
+save('lineXrag.out',x)
+save('lineYrag.out',y)
+save('lineZrag.out',z)
+
+temps=arange(973.0,1274,50)
+x=1000.0/temps
+y=log(array([OneAtm for T in temps]))/log(10.0)
+z=array([PratsinisLogK(T) for T in temps])
+ax.plot3D(x,y,z,label="Pratsinis (1990)", linewidth=2)
+
+PratData=array([[ 0.883, 1.1112 ],[ 0.8981, 1.0121 ],[ 0.9149, 0.9306 ],[ 0.9317, 0.6962 ],[ 0.9499, 0.4882 ],[ 0.9677, 0.1884 ],[ 0.9832, -0.2687 ],[ 1.0073, -0.2841 ]])
+x=PratData[:,0]
+y=log(array([OneAtm for t in x]))/log(10.0)
+z=PratData[:,1]
+ax.plot3D(x,y,z,'o-',label="Pratsinis (1990) 1:5", linewidth=2)
+save('lineXprat.out',x)
+save('lineYprat.out',y)
+save('lineZprat.out',z)
+
+PratsinisMeasuredInverseTs=x
+PratsinisMeasuredPressures=y
+PratsinisMeasuredLogKs=z
+
+
+def plotRaghavanMeasurements():
+    temps=array([  646.,   757.,  1047.,  1074.,  1152.,  1258.,  1299.,  1407.,
+            1440.,  1516.,  1607.,  1648.,  1718.,  1743.,  1809.,  2044.])       
+    measuredNNos=array([ 1.  ,  0.97,  1.02,  0.9 ,  0.96,  1.04,  0.94,  0.49,  0.44,
+            0.6 ,  0.33,  0.35,  0.16,  0.24,  0.17,  0.17])
+    reactiontimes=array([ 0.467,  0.166,  0.058,  0.063,  0.04 ,  0.032,  0.032,  0.027,
+            0.023,  0.024,  0.019,  0.019,  0.018,  0.018,  0.013,  0.01 ])
+    pressures=array([  165163.21252451,   252971.75364402,   551205.35156326,
+             429845.99298875,   736792.54475485,   964331.41922609,
+             887997.53785585,  1198889.6613466 ,   942090.5684987 ,
+            1436794.95733164,  1706284.57763703,  1810558.77650401,
+            1458560.42080549,  1526135.31697436,  1695978.55319514,
+            2384852.8012917 ])
+            
+    x=1000.0/temps 
+    for i in range(temps.size):
+        temp=temps[i]
+        measuredNNo=measuredNNos[i]
+        pressure=pressures[i]
+        reactiontime=reactiontimes[i]
+        errorInNNo=0.1
+        x=1000.0/temp
+        y=RaghavanLogKexpt( temp, measuredNNo, reactiontime )
+        ymin =RaghavanLogKexpt( temp, measuredNNo+errorInNNo, reactiontime )
+        ymax =RaghavanLogKexpt( temp, measuredNNo-errorInNNo, reactiontime )
+        ybar=array([ymin,ymax])
+        xbar=array([x, x])
+        if y>9:
+            plot( [x], [y], 'b^', label="Raghavan", linewidth=1 )
+        elif y<-5:
+            plot( [x], [y], 'bv', label="Raghavan", linewidth=1 )
+        else:
+            plot( [x], [y], 'ro', label="Raghavan", linewidth=1 )
+        
+        plot( xbar, ybar, 'b-_', label="Raghavan", linewidth=1 )
+        print "Tp\t%dK\tNNo\t%g\tt\t%g\tk_eff\t%.4g\tmin\t%.4g\tmax\t%.4g"%(temp,measuredNNo, reactiontime, y,ymin,ymax)
+    
+def plotPratsinisMeasurements():
+    plot( PratsinisMeasuredInverseTs, PratsinisMeasuredLogKs, 'go', label="Pratsinis", linewidth=2 )
+
+
+X=1000.0/array(meshTemps)
+Y=log(array(meshPressures))/log(10.0)
+Z=array(meshLogKs)
+surf1=ax.plot_surface(X,Y,Z)
+frame1=ax.plot_wireframe(X,Y,Z)
+
+Zprat=array([[PratsinisLogK(1000.0/i) for i in j] for j in X])
+surf2=ax.plot_surface(X,Y,Zprat)
+frame2=ax.plot_wireframe(X,Y,Zprat)
+frame1._wrapped.set_color('r')
+
+Zrag=array([[RaghavanLogK(1000.0/i) for i in j] for j in X])
+surf3=ax.plot_surface(X,Y,Zrag)
+frame3=ax.plot_wireframe(X,Y,Zrag)
+frame2._wrapped.set_color('g')
+
+save('meshX.out', X)
+save('meshY.out', Y)
+save('meshZme.out',Z)
+save('meshZprat.out',Zprat)
+save('meshZrag.out',Zrag)
+
+
+for surf in [surf1, surf2, surf3]:
+    surf._wrapped.set_alpha(0.2)
+    surf._wrapped.set_linewidth(0)
+    
+for frame in [frame1, frame2, frame3]:
+    frame._wrapped.set_alpha(0.5)
+    
+ax.set_xlabel('1000/T')
+ax.set_ylabel('log_{10} P')
+ax.set_zlabel('ln k')
+
+
+## now add things
+figure(1)
+if clearfigs:
+    plotRaghavanMeasurements()
+    plot( PratsinisMeasuredInverseTs, PratsinisMeasuredLogKs, 'go', label="Pratsinis", linewidth=2 )
+
+#### O2 
+figure(3)
+if clearfigs: clf()
+P=1.0e5
+title("Effect of excess O2 (no argon)")
+o2xs=[-1,0,1,4,9]
+for o2x in o2xs:
+    temps=[]
+    logKs=[]
+    if o2x==-1:
+        if not clearfigs: continue
+        label="Raghavan (2001)"
+        for T in [646, 2044]:
+            temps.append(T)
+            logKs.append(RaghavanLogK(T)) 
+    elif o2x==0:
+        if not clearfigs: continue
+        label="Pratsinis (1990)"
+        for T in [973, 1273]:
+            temps.append(T)
+            logKs.append(PratsinisLogK(T))
+    else:
+        for T in Ts:
+            #label="O2/TiCl4=%.2g"%(o2x)
+            ar=0.0
+            ticl=(1.0-ar)/(1+o2x)
+            o2=ticl*o2x
+            X="TiCl4:%f, O2:%f, Ar:%f"%(ticl,o2,ar)
+            label=X
+            s=Simulation(T,P,X)
+            print "T\t%dK\tO2/TiCl4\t%.2g\t"%(T,o2x),
+            logK=s.getLogK()
+            if logK>-Inf:
+                temps.append(T)
+                logKs.append(logK)
+    temps=array(temps)
+    logKs=array(logKs)
+    plot( 1000.0/temps, logKs, format, label=label, linewidth=2 )
+    
+
+
+ylabel("ln k")
+xlabel("1000/T")
+if clearfigs:
+    legend(loc=0)
+    drawPratsinisBox()
+    drawRaghavanBox()
+    plotRaghavanMeasurements()
+    plot( PratsinisMeasuredInverseTs, PratsinisMeasuredLogKs, 'go', label="Pratsinis", linewidth=2 )
+
+figure(4)
+if clearfigs: clf()
+title("Effect of argon fraction (4:1 O2:TiCl4)")
+P=1.0e5
+argons=[-2,-1,0,0.9,0.99]
+o2x=4
+for argon in argons:
+    temps=[]
+    logKs=[]
+    if argon==-2:
+        if not clearfigs: continue
+        label="Raghavan (2001)"
+        for T in [646, 2044]:
+            temps.append(T)
+            logKs.append(RaghavanLogK(T)) 
+    elif argon==-1:
+        if not clearfigs: continue
+        label="Pratsinis (1990)"
+        for T in [973, 1273]:
+            temps.append(T)
+            logKs.append(PratsinisLogK(T))
+    else:
+        for T in Ts:
+            #label="O2/TiCl4=%.2g"%(o2x)
+            ar=argon
+            ticl=(1.0-ar)/(1+o2x)
+            o2=ticl*o2x
+            X="TiCl4:%f, O2:%f, Ar:%f"%(ticl,o2,ar)
+            label=X
+            s=Simulation(T,P,X)
+            print "T (K)\t%d\tArgon\t%.2g\t"%(T,argon),
+            logK=s.getLogK()
+            if logK>-Inf:
+                temps.append(T)
+                logKs.append(logK)
+    temps=array(temps)
+    logKs=array(logKs)
+    plot( 1000.0/temps, logKs, format, label=label, linewidth=2 )
+
+
+
+ylabel("ln k")
+xlabel("1000/T")
+if clearfigs:
+    legend(loc=0)
+    drawPratsinisBox()
+    drawRaghavanBox()
+    plotRaghavanMeasurements()
+    plot( PratsinisMeasuredInverseTs, PratsinisMeasuredLogKs, 'go', label="Pratsinis", linewidth=2 )
diff --git a/mixmaster.py b/mixmaster.py
new file mode 100755
index 00000000..f7133632
--- /dev/null
+++ b/mixmaster.py
@@ -0,0 +1,4 @@
+#! python
+from MixMaster import MixMaster
+o =  MixMaster()
+
diff --git a/prettyDot.py b/prettyDot.py
new file mode 100644
index 00000000..7fce9241
--- /dev/null
+++ b/prettyDot.py
@@ -0,0 +1,42 @@
+# this version updated by RHW in November 2008
+"""
+make the dot file prettier
+"""
+import os, sys, re
+
+StripLineLabels=True
+if StripLineLabels:
+	print "stripping edge (line) labels"
+
+infile=file('rxnpath.dot')
+outfile=file('rxnpath2.dot','w')
+
+# replace this:
+#  s10 [ fontname="Helvetica", label="C11H23J"];
+# with this:
+#  s10 [ shapefile="mols/C11H23J.png" label="" width="1" height="1" imagescale=true fixedsize=true color="white" ];
+ 
+reSize=re.compile('size=\"5,6\"\;page=\"5,6\"')
+reNode=re.compile('(?P<node>s\d+)\ \[\ fontname=\"Helvetica\",\ label=\"(?P<label>[^\"]*)\"\]\;')
+
+for line in infile:
+	(line,changedSize)=reSize.subn('size="12,12";page="12,12"',line)
+	match=reNode.search(line)
+	if match:
+		if os.path.isfile("pics/%s.png"%match.group('label')):
+			line='%s [ image="pics/%s.png" label="" width="0.5" height="0.5" imagescale=false fixedsize=false color="none" ];\n'%(match.group('node'),match.group('label'))
+
+	# rankdir="LR" to make graph go left>right instead of top>bottom
+
+	if StripLineLabels:
+		line=re.sub('label=\"\s*[\d.]+\"','label=""',line)
+		
+	# change colours
+	line=re.sub('color="0.7,\ (.*?),\ 0.9"',r'color="1.0, \1, 0.7*\1"',line)
+		
+	outfile.write(line)
+
+outfile.close()
+infile.close()
+        
+print "now try:\n dot -O -Tsvg rxnpath2.dot"
\ No newline at end of file
diff --git a/templates/rxnlist.html b/templates/rxnlist.html
new file mode 100644
index 00000000..96e8a507
--- /dev/null
+++ b/templates/rxnlist.html
@@ -0,0 +1,39 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN">
+<html lang="en">
+<head>
+    <title>{{ title }}</title>
+	<link rel="stylesheet" type="text/css" href="templates/stylesheet.css" media="all">
+</head>
+<body>
+	<h1>{{ title }}</h1>
+	<h2>{{ reactionList|length }} reactions</h2>
+
+<table class="reactionList">
+{% for reaction in reactionList %}
+ <tr class="reactionRow">
+	<td class="reactionNumber"> {{ reaction._num }} </td>
+  <td class="reactionSide">
+   {% for species in reaction._r %}
+    {% if reaction._r[species]!=1 %}
+     {{ reaction._r[species] }}
+    {% endif %}
+    <img src="pics/{{ species }}.png" class="speciesPic" >
+    {% if not loop.last %} + {% endif %}
+   {% endfor %}
+  </td>
+  <td>=</td>
+  <td class="reactionSide">
+   {% for species in reaction._p %}
+    {% if reaction._p[species]!=1 %}
+     {{ reaction._p[species] }}
+    {% endif %}
+    <img src="pics/{{ species }}.png" class="speciesPic" >
+    {% if not loop.last %} + {% endif %}
+   {% endfor %}
+  </td>
+ </tr>
+{% endfor %}
+</table>
+
+</body>
+</html>
\ No newline at end of file
diff --git a/templates/stylesheet.css b/templates/stylesheet.css
new file mode 100644
index 00000000..82f55b5e
--- /dev/null
+++ b/templates/stylesheet.css
@@ -0,0 +1,57 @@
+body {
+background-color: #ffffff;
+color: #111;
+font-size: 12px;
+font-family: Verdana, Arial, SunSans-Regular, Sans-Serif;
+}
+
+table.reactionList {
+	border-collapse: collapse;
+	font-size: 16px;
+}
+
+
+tr.reactionRow td{
+	border-bottom:1px solid #111;
+}
+
+.residuals {
+	font-size: 10px;
+	float:left;
+}
+
+table.groupList {
+	clear:both;
+	font-size: 10px;
+	border: 1px solid #111;
+	width: 100%;
+}
+tr.groupRow td {
+	border-bottom-width: 0px;
+}
+
+tr.unmatchedAtomsRow td {
+	border-bottom-width: 0px;
+}
+
+table.multimatchedAtomsList {
+	font-size: 10px;
+	border: 1px solid #111;
+	width: 100%;
+	background: #fb8;
+}
+
+tr.multimatchedAtomsRow td {
+	border-bottom: 1px solid #111;
+}
+
+img.speciesPic {
+	vertical-align: middle;
+}
+td.speciesPic {
+	text-align: center;
+}
+td.reactionSide{
+	vertical-align: middle;
+	text-align: center;
+}
\ No newline at end of file
diff --git a/toMixMaster.py b/toMixMaster.py
new file mode 100644
index 00000000..ebb0a31f
--- /dev/null
+++ b/toMixMaster.py
@@ -0,0 +1,107 @@
+# this version updated by RHW in January 2009
+"""
+load the chem.inp file from RMG and get it ready for mixmaster
+"""
+import os, sys, shutil
+#from Cantera import *
+#from Cantera.Reactor import *o
+#from Cantera.Func import *
+#from math import *
+#import random
+#import pylab
+#import scipy
+#from pylab import *	
+#from scipy import *
+#from scipy.optimize import leastsq
+
+#RMGworkingDir="/Users/rwest/XCodeProjects/RMG/ExampleWorkFolder/"
+RMGworkingDir="/Volumes/Rainier/RMG3/WorkingFolders/dieselWithOOH/"
+#RMGworkingDir="/Users/rwest/XCodeProjects/RMG2-sounceforge/WorkingDirectory"
+#RMGworkingDir="/Users/rwest/XCodeProjects/RMG/WorkingFolders/v3_from_CVS/diesel_surrogate_.5mod"
+
+RMGworkingDir="/Users/rwest/Documents/AcaPostdoc/ENI/from Amrit/diesel surrogate runs/GP tol=0.1 time=4hrs"
+
+# convert the chemkin file from RMG into a cantera file chem.cti
+from Cantera import ck2cti
+infile='chem.inp'
+oldpath=os.path.join(RMGworkingDir,'chemkin',infile)
+newpath=os.path.join(os.getcwd(),infile)
+print "copying %s to %s"%(oldpath,newpath)
+shutil.copy2(oldpath, newpath) # copy it to the current folder
+
+thermodb=''
+trandb=''
+nm='chem'
+ck2cti.ck2cti(infile = infile, thermodb = thermodb,  trandb = trandb, idtag = nm, debug=0, validate=1)
+
+
+# convert the Final_Model.txt into approprita CSV file
+
+print "NB  ForMixMaster.csv is wrong. MixMaster wants MASS fractions and we are giving it MOLE fractions" 
+
+temperature=273+150
+pressure=208*101325
+print " using these settings:\n Temperature: %f K \t Pressure: %f Pa\n"%(temperature,pressure)
+
+# load file
+resultFile=file(os.path.join(RMGworkingDir,'Final_Model.txt'))
+# search for "Mole Fraction Profile Output"
+line=resultFile.next()
+while (line.find('Mole Fraction Profile Output')<0):
+    line=resultFile.next()
+# add "T \t P" to the  following line
+titles=resultFile.next()
+print "Species:",titles
+output=titles.strip()+"\tT\tP\tnothing\n"
+# add the temperature IN KELVIN and pressure IN PASCAL to all the following nonblank lines
+line=resultFile.next()
+while (line.strip()):
+    output += line.strip() + "\t%f\t%f\t0\n"%(temperature,pressure)
+    line=resultFile.next()
+# turn whitespaces into commas
+# save the output
+outputFile=file('ForMixMaster.csv','w')
+outputFile.write(output.replace('\t',','))
+outputFile.close()
+
+
+### make pretty table of species
+import ctml_writer
+from ctml_writer import *
+# these lists store top-level entries
+ctml_writer._elements = []
+ctml_writer._species = []
+ctml_writer._speciesnames = []
+ctml_writer._phases = []
+ctml_writer._reactions = []
+ctml_writer._atw = {}
+ctml_writer._enames = {}
+ctml_writer._valsp = ''
+ctml_writer._valrxn = ''
+ctml_writer._valexport = ''
+ctml_writer._valfmt = ''
+
+execfile('chem.cti')
+
+for rxn in ctml_writer._reactions:
+	for side in [rxn._r,rxn._p]:
+		for sp in side:
+			num=side[sp]
+			if num!=1: print num,
+			print sp
+		
+import jinja2
+env = jinja2.Environment(loader = jinja2.FileSystemLoader('templates'))
+template = env.get_template('rxnlist.html')
+outstring=template.render(title=infile, reactionList=ctml_writer._reactions)
+outfile=file('ReactionList'+'.html','w')
+outfile.write(outstring)
+outfile.close()
+
+print outstring
+
+
+# load mixmaster
+from MixMaster import MixMaster
+#o=MixMaster()
+#o.loadmech('','chem.cti')
\ No newline at end of file
diff --git a/toMolFile.py b/toMolFile.py
new file mode 100644
index 00000000..0f49d9eb
--- /dev/null
+++ b/toMolFile.py
@@ -0,0 +1,110 @@
+# this version updated by RHW in November 2008
+"""
+load the chem.inp file from RMG and get it ready for mixmaster
+"""
+import os, sys, shutil, re
+#from Cantera import *
+#from Cantera.Reactor import *o
+#from Cantera.Func import *
+#from math import *
+#import random
+#import pylab
+#import scipy
+#from pylab import *
+#from scipy import *
+#from scipy.optimize import leastsq
+
+import openbabel, pybel
+# please cite:
+# Pybel: a Python wrapper for the OpenBabel cheminformatics toolkit
+# Noel M O'Boyle, Chris Morley and Geoffrey R Hutchison
+# Chemistry Central Journal 2008, 2:5
+# doi:10.1186/1752-153X-2-5
+
+
+
+picfolder='pics'
+molfolder='mols'
+for path in [picfolder,molfolder]:
+    os.path.isdir(path) or os.mkdir(path)
+
+
+periodicTableByNumber={ 1: 'H',  2: 'He',  3: 'Li',  4: 'Be',  5: 'B',  6: 'C',  7: 'N',  8: 'O',  9: 'F',  10: 'Ne',  11: 'Na',  12: 'Mg',  13: 'Al',  14: 'Si',  15: 'P',  16: 'S',  17: 'Cl',  18: 'Ar',  19: 'K',  20: 'Ca',  21: 'Sc',  22: 'Ti',  23: 'V',  24: 'Cr',  25: 'Mn',  26: 'Fe',  27: 'Co',  28: 'Ni',  29: 'Cu',  30: 'Zn',  31: 'Ga',  32: 'Ge',  33: 'As',  34: 'Se',  35: 'Br',  36: 'Kr',  37: 'Rb',  38: 'Sr',  39: 'Y',  40: 'Zr',  41: 'Nb',  42: 'Mo',  43: 'Tc',  44: 'Ru',  45: 'Rh',  46: 'Pd',  47: 'Ag',  48: 'Cd',  49: 'In',  50: 'Sn',  51: 'Sb',  52: 'Te',  53: 'I',  54: 'Xe',  55: 'Cs',  56: 'Ba',  57: 'La',  58: 'Ce',  59: 'Pr',  60: 'Nd',  61: 'Pm',  62: 'Sm',  63: 'Eu',  64: 'Gd',  65: 'Tb',  66: 'Dy',  67: 'Ho',  68: 'Er',  69: 'Tm',  70: 'Yb',  71: 'Lu',  72: 'Hf',  73: 'Ta',  74: 'W',  75: 'Re',  76: 'Os',  77: 'Ir',  78: 'Pt',  79: 'Au',  80: 'Hg',  81: 'Tl',  82: 'Pb',  83: 'Bi',  84: 'Po',  85: 'At',  86: 'Rn',  87: 'Fr',  88: 'Ra',  89: 'Ac',  90: 'Th',  91: 'Pa',  92: 'U',  93: 'Np',  94: 'Pu',  95: 'Am',  96: 'Cm',  97: 'Bk',  98: 'Cf',  99: 'Es',  100: 'Fm',  101: 'Md',  102: 'No',  103: 'Lr',  104: 'Rf',  105: 'Db',  106: 'Sg',  107: 'Bh',  108: 'Hs',  109: 'Mt',  110: 'Ds',  111: 'Rg',  112: 'Uub',  113: 'Uut',  114: 'Uuq',  115: 'Uup',  116: 'Uuh',  117: 'Uus',  118: 'Uuo'}
+periodicTableBySymbol=dict([(val, key) for key, val in periodicTableByNumber.items()])   
+
+OBMolBondTypes={'S':1, 'D':2, 'T':3, 'B':5 }
+
+RMGworkingDir="/Users/rwest/XCodeProjects/RMG/ExampleWorkFolder/"
+RMGworkingDir="/Volumes/Rainier/RMG3/WorkingFolders/dieselWithOOH/"
+#RMGworkingDir="/Users/rwest/XCodeProjects/RMG2-sounceforge/WorkingDirectory"
+#RMGworkingDir="/Users/rwest/XCodeProjects/RMG/WorkingFolders/v3_from_CVS/diesel_surrogate_.5mod"
+RMGworkingDir="/Users/rwest/Documents/AcaPostdoc/ENI/from Amrit/diesel surrogate runs/GP tol=0.1 time=4hrs"
+
+# copy the RMG dictionary file
+infile='RMG_Dictionary.txt'
+oldpath=os.path.join(RMGworkingDir,infile)
+newpath=os.path.join(os.getcwd(),infile)
+print "copying %s to %s"%(oldpath,newpath)
+shutil.copy2(oldpath, newpath) # copy it to the current folder
+
+# load file
+RMGfile=file(newpath)
+
+
+for i in range(1,30000):
+    print 'Molecule', i,'\t',
+    name=''
+    try:
+        while name=='':
+            name=RMGfile.next().strip()
+    except StopIteration:
+        print 'No more molecules'
+        break
+    print name
+    graph=[]
+    line=RMGfile.next()
+    while line.strip():
+        graph.append(line)
+        line=RMGfile.next()
+    # now have 'name' and 'graph'
+    
+    mol = openbabel.OBMol()
+    #print 'Should print 0 (atoms)'
+    #print mol.NumAtoms()
+
+    re_bond=re.compile('\{(?P<atomnum>\d+),(?P<bondtype>[SDTB])\}')
+    for line in graph:
+        #print 'line:',line.strip()
+        if len(line.split())>3:
+            (number, element, radical, bonds)=line.split(None,3)
+        else:
+            (number, element, radical )=line.split(None)
+        a = mol.NewAtom()
+        a.SetAtomicNum(periodicTableBySymbol[element])  # 6 for a carbon atom
+        #a.SetVector(0.0, 1.0, 2.0) # coordinates
+        if int(radical[0]): # the [0] is so we take the first character of the string, in case it's something like "2T"
+            a.SetSpinMultiplicity(int(radical[0])+1)
+            # note that for non-radicals it's 0, but single radicals are 2, double radicals are 3...
+            # http://openbabel.org/wiki/Radicals_and_SMILES_extensions#How_OpenBabel_does_it
+        for bond in bonds.split():
+            #print bond
+            matchobject=re_bond.match(bond)
+            if matchobject:
+                fromAtom=int(number)
+                toAtom=int(matchobject.group('atomnum'))
+                bondType=matchobject.group('bondtype')
+                if toAtom>fromAtom:
+                    continue # because toAtom hasn't been placed yet!
+                # print "%s bond from %d to %d"%(bondType,fromAtom,toAtom)
+                mol.AddBond(fromAtom,toAtom,OBMolBondTypes[bondType])
+            else:
+                raise "couldn't figure out this bond: %s"%bond
+    pymol=pybel.Molecule(mol)
+    print pymol.write().strip(), 
+    chemkinformula=pymol.formula+'J'*(pymol.spin-1)
+    print chemkinformula
+    if pymol.OBMol.NumHvyAtoms()>1:
+        pymol.removeh()
+    pymol.draw(filename=os.path.join('pics',name+'.png'), update=True, show=False)
+    pymol.write(format='mol',filename=os.path.join('mols',name+'.mol'),overwrite=True)
+