#pragma rtGlobals=3		// Use modern global access method and strict wave access.

// Modified by Ken Aikin from code that Steve Brown had for doing 1-D chemical model of plume.
// Modified code now does dispersion in X and Y with chemistry.
// To run, choose Macros->Run_Model2()
// Prompts for number of downwind rows, and cross-wind columns, as well as timestep for each iteration, and the dispersion fraction.
// This dispersion fraction is used to disperse the plume by colvolving that wave with the concentrations for each species.
// Dispersion is calculated with the wave {dispersion fraction, 1-2*(dispersion fraction), dispersion fraction} i.e. {0.05, 0.9, 0.05}
// The sum of that wave should always be 1 to conserve mass.
// KA 20151102

// Mods:
// Clean up indentation, and replace many spaces with tabs in DispersionMechanismCompile(). Delete empty "else" statements.
// Add function removeDuplicatesInTextWave() to cleanup Names wave with possible duplicates.
// Remove unused variables in RunDispersionModel()	// T and Prs
// KA 20160510

// Add columns to table after running Mechanism_Setup.
// Need to change Macro RateConst() to be run sometime (when T and/or Prs change) and change macro to function.
// Add new wave BackgroundConc to table. Use this instead of having to add variables to code for passing in values.
// Add the menu items to "Macros" menu for running Chem/Dispersion model.
// KA 20160511

// Fix bugs in DispersionMechanismCompile() that caused out of bounds error for Names wave while looking for matching string in reactants and products waves.
// Put in break statements as more graceful way to exit loops.
// KA 20160512

Menu "Macros"
	"---"
	"Dispersion Mechanism_Setup"
	"Dispersion Mechanism Compile"
	"Run Dispersion Model"
	"DoManyModelRuns"
	"transect"
	"---"
end


// Macro to call RunDispersionModel() repeatedly for doing lots of runs. Store runs in data folders using SaveModelRun() in Results.ipf
// Specific case for varying N2O5_LossRate,MonoBkg,AlkBkg,AldBkg,Dispersion,NO_emissions,O3_background in model runs.
// Make it more general later?
macro DoManyModelRuns()
        variable refNum, n, num
        string runName

        Print "Select a comma-delimited file containing model run parameters"
        NewDataFolder/O/S root:MultipleRuns
        Loadwave/J/A/W/O                // prompt for comma-delimited file of parameters for model runs.

        num = numpnts(N2O5_LossRate)
        SetDataFolder root:

        do
                // set up parameters loaded from comma-delimited file above
                BackgroundConc[5] = root:MultipleRuns:MonoBkg[n]                        // CHANGE THESE to row number of background value
                BackgroundConc[10] = root:MultipleRuns:AlkBkg[n]                         // CHANGE THESE
                BackgroundConc[12] = root:MultipleRuns:AldBkg[n]                         // CHANGE THESE
                BackgroundConc[1] = root:MultipleRuns:O3_background[n]          // CHANGE THESE
                RateCoefficients[5] = root:MultipleRuns:N2O5_LossRate[n]*(1-root:MultipleRuns:ClNO2_)
                RateCoefficients[6] = root:MultipleRuns:N2O5_LossRate[n]*(root:MultipleRuns:ClNO2_)
                runName = "run_"+num2str(root:MultipleRuns:N2O5_loss[n])+"_"+num2str(root:MultipleRuns:Mon[n])+"_"+num2str(root:MultipleRuns:Disp[n])+"_"+num2str(root:MultipleRuns:NO_emissions[n])+"_"+num2str(root:MultipleRuns:O3_background[n])
                //sprintf runName, "run_%03d", n  // do simple miethod for now. Make n+last# if running multiple times in one experiment
                PRINT runName           // debug
                DispersionMechanismCompile("ChemKinetic")
                RunDispersionModel(14, root:MultipleRuns:Dispersion[n], 60, root:MultipleRuns:NO_emissions[n], .1)       // RunDispersionModel(TotalTime, DiffPerSec, TimeStep, NOconc, NOwidth) // make sure these are right for the defaults
                ExtractSpeciesFromLayer("All")
                transect(720)
                Display NOy_transect
                append NO_transect, NO2_transect, NO3_transect, N2O5_transect, HNO3_transect, ClNO2_transect, MonNO3_transect, AlkNO3_transect, Pan_transect, orgN_transect
                Make/O/N=17 Areas
                InitializeAlmostEverything("NO_transect","_calculated_","_none_")
                InitAreaBetweenCursors1(3,"NO_transect","_calculated_","_None_",1)
                AreaBetweenCursors1()
                Areas[0] = W_AreaNoBase[0]
                 InitAreaBetweenCursors1(3,"NO2_transect","_calculated_","_None_",1)
                AreaBetweenCursors1()
                Areas[2] = W_AreaNoBase[0]
                 InitAreaBetweenCursors1(3,"NO3_transect","_calculated_","_None_",1)
                AreaBetweenCursors1()
                Areas[3] = W_AreaNoBase[0]
                 InitAreaBetweenCursors1(3,"N2O5_transect","_calculated_","_None_",1)
                AreaBetweenCursors1()
                Areas[4] = W_AreaNoBase[0]
                 InitAreaBetweenCursors1(3,"HNO3_transect","_calculated_","_None_",1)
                AreaBetweenCursors1()
                Areas[8] = W_AreaNoBase[0]
                 InitAreaBetweenCursors1(3,"ClNO2_transect","_calculated_","_None_",1)
                AreaBetweenCursors1()
                Areas[9] = W_AreaNoBase[0]
                 InitAreaBetweenCursors1(3,"MonNO3_transect","_calculated_","_None_",1)
                AreaBetweenCursors1()
                Areas[6] = W_AreaNoBase[0]
                 InitAreaBetweenCursors1(3,"AlkNO3_transect","_calculated_","_None_",1)
                AreaBetweenCursors1()
                Areas[11] = W_AreaNoBase[0]
                 InitAreaBetweenCursors1(3,"Pan_transect","_calculated_","_None_",1)
                AreaBetweenCursors1()
                Areas[14] = W_AreaNoBase[0]
                InitAreaBetweenCursors1(3,"orgN_transect","_calculated_","_None_",1)
                AreaBetweenCursors1()
                Areas[15] = W_AreaNoBase[0]
                InitAreaBetweenCursors1(3,"NOy_transect","_calculated_","_None_",1)
                AreaBetweenCursors1()
                Areas[16] = W_AreaNoBase[0]
                SaveModelRun(runName)
                n+=1
        while(n<num)
end



// This is the main macro to run after doing the initial setup to compile the chemical formulas.
// Then just choose number of downwind and cross-wind rows and columns, and the time step.
Macro RunDispersionModel(TotalTime, DiffPerSec, TimeStep, NOconc, NOwidth)	// O3bkg
	Variable TotalTime=2, DiffPerSec=5e-4, TimeStep=60, NOconc=40,NOwidth=1	//O3bkg=50, // unused: Prs=900, T=285
	Prompt TotalTime, "Model run time (hrs)"
	//Prompt NumCols, "Number of rows for plume matrix (cross-wind, 0.1km)"
	Prompt DiffPerSec, "Dispersion Fraction (1e-4 m2/s)"
	Prompt TimeStep, "Time step (sec)"
//	Prompt O3bkg, "O3 background (ppbv)"
	Prompt NOconc, "NO plume conc (ppbv)"
	Prompt NOwidth, "NO plume width (km)"
		
	Silent 1; PauseUpdate
	
	Variable NumRows = Totaltime*60*60/Timestep	//Calculate number of rows needed for total time and timestep
	Variable DiffFraction = DiffperSec*Timestep		//Convert diffusion from per second to per timestep
	Variable NumCols=201								//If you change this, you have to change the plume width/box calculation in the function below
//	Variable/G NDensity								//Number density from rate constant setup	//KA THIS IS GLOBAL SO DON'T NEED TO PASS INTO DORUN_MODEL2
			
	
	If (!exists("CurrChemKineticFct"))
		Abort "Must compile ODE first!"
	Endif
	
	Print "\r\rDoing model run with ",NumRows, " rows and ",NumCols, " columns, and time step of ",TImeStep, " with dispersion of ",DiffPerSec
	Print "Initial conditions: NO conc: ", NOconc, " NO plume width: ", NOwidth	// print all other non-zero bkg values also?
	
	DoRun_Model2(NumRows, NumCols, TImeStep, DiffFraction, NOconc, NOwidth)	// O3bkg, NDensity
	
	Print "Refreshing image matrices..."
	RefreshImageMatrices()
	Print "Done"
end

function DoRun_Model2(NumRows, NumCols, TimeStep, DiffFraction, NOconc, NOwidth)		// O3bkg, NDensity
	Variable NumRows, NumCols, TimeStep, DiffFraction, NOconc, NOwidth					// O3bkg, NDensity

	Variable NumSpecies, n, m, zLayer, starttime, fract
	Variable NOconcPL	// variable for initial plume concentration.
	String cmd, SpeciesList, CurrSpecies
	wave/T Names		// the waves listed in 'Names' get the returned results from the chem reactions after calling IntegrateODE
	wave InitialConcentrations, BackgroundConc, CalTime

	
	SVAR CurrChemKineticFct
	NVAR NDensity				//Number density from rate constant setup. WHEN DOES THAT FUNCTION GET RUN??
		
	starttime = DateTime
	SpeciesList = ""
	NumSpecies = numpnts(Names)
	fract = NumRows/10	// define interval for messages within loop to see progress.

	// initial concentrations.

//	O3concBK = O3bkg*NDensity*1e-9	//Convert O3 background to number density
	NOconcPL = NOconc*NDensity*1e-9	//Convert NO concentration to number density
	
	if (DiffFraction>1)
		Abort "Error in DiffFraction"
	endif
	
	// Z is in order of Names wave
	Make/D/O/N=(NumRows, NumCols, NumSpecies) resultsMatrix=0	// concentrations.  Row=downwind distance. Column=crosswind dist. Layer=species.
	Make/D/O/N=(NumCols) DownWindSlice
//	Make/D/O/N=3 sm_coefs={0.05, 0.9, 0.05}	// wave that gets convolved with concentration to do downwind dispersion. Default
	Make/D/O/N=3 sm_coefs={DiffFraction, 1-2*DiffFraction, DiffFraction}	// wave that gets convolved with concentration to do downwind dispersion. Larger middle number is less dispersion. Wave needs to add up to 1.
	CalTime[1] = TimeStep
		
	for(n=1; n<NumSpecies; n+=1)
		resultsMatrix[0][][n] = BackgroundConc[n]	*Ndensity*1e-9		// background level everywhere crosswind in first row.
	endfor

// init plume matrix with background values for NO and one cell for NO plume emission, as well as background values for all other species.
	Variable p1 = (NumCols-1)/2-round(NOwidth/2*10)
	Variable p2 = (NumCols-1)/2+round(NOwidth/2*10)	
	// ASSUME NO IS ALWAYS IN ROW 0. 
	resultsMatrix[0][p1,p2][0] = NOconcPL	//NO initial plume concentration in middle of crosswind dist. All these layer numbers need to be updated if the order of "NAMES" list changes
	resultsMatrix[0][0,(p1-1)][0] = 0		//NO background, set to 0
	resultsMatrix[0][(p2+1),200][0] = 0		//NO background, set to 0

	
	for (n=0; n<NumSpecies; n+=1)					// Run once outside of loop to set up speciesList for the IntegrateODE command.
		SpeciesList += Names[n] + ","
	endfor
	SpeciesList = RemoveEnding(SpeciesList, ",")				// trim off trailing comma
	
	for (n=1; n<NumRows; n+=1)								// row (downwind distance)
		if (mod(n,fract)==0)									// n is a multiple of fract (i.e. print out every 10 percent of completion of loop.)
			print "Loop is ",100*n/NumRows, "% complete"
		endif
		for (m=0; m<NumCols; m+=1)							// col (crosswind distance)
			InitialConcentrations = resultsMatrix[n-1][m][p]	// Get current values into InitialConcentrations wave
			for (zLayer=0; zLayer<NumSpecies; zLayer+=1)	// transfer from InitConc to first wave element of NO, NO2, etc.
				wave sp = $Names[zLayer]						// set wave sp to the name in text wave Names
				sp[0] = InitialConcentrations[zLayer]			// put conc from InitialConcentrations into NO, NO2, etc in preparation for running chemistry.
			endfor

			// Do chemistry
			cmd= "IntegrateODE/U=1000/X=CalTime "+CurrChemKineticFct+" KK, {"+ SpeciesList + "}"	// chem mechanism. Puts results in waves NO, O3, NO2, NO3, N2O5, RONO2, HNO3. Make sure they are same length as CalTime
			Execute cmd
			
			// Store results in resultsMatrix
			for (zLayer=0; zLayer<NumSpecies; zLayer+=1)
				wave sp = $Names[zLayer]					// set wave sp to the name in text wave Names
				resultsMatrix[n][m][zLayer] = sp[1]		// loop thru all the species waves and store the new conc after running the chemistry.
			endfor
			
		endfor	// end of a crosswind slice (one row is complete)
		
		// Do dispersion for whole row of matrix for each species
		for (zLayer=0; zLayer<NumSpecies; zLayer+=1)
			DownWindSlice = resultsMatrix[n][p][zLayer]
			FilterFIR/E=0/COEF=sm_coefs DownWindSlice		// this is a custom smooth operation which convolves the wave with the coeficients sm_coefs.
			// put dispersed results (DownWindSlice) back into resultsMatrix
			resultsMatrix[n][][zLayer] = DownWindSlice[q]	// overwrite the chemistry-only data that was in this row with the chem&disp data.

		endfor
	endfor
	resultsMatrix[][][] *= 1e9/NDensity
	print "Calculation time", DateTime - starttime,"seconds"
end

// Run this to refresh all matrices that have already been extracted from the resultsMatrix after running the model.
function RefreshImageMatrices()
	//variable timestep
	variable n, len, numRows, numCols
	string NamesList, species, speciesMat
	wave/T Names
	wave resultsMatrix
	NVAR NDensity
	
	numRows = DimSize(resultsMatrix, 0)
	numCols = DimSize(resultsMatrix, 1)

	//variable timedownwind = numrows*timestep

	NamesList = TextWave2List(Names)
	len = numpnts(Names)
	for(n=0; n<len; n+=1)
		species = StringFromList(n, NamesList)
		speciesMat = species+"_mat"
		if (exists(speciesMat)==1)
			wave w = $speciesMat
			Redimension/N=(numrows, numCols) w
			w = resultsMatrix[p][q][n]
			//w *= 1e9/NDensity					//Convert to ppbv units
		endif
	endfor
end


Macro DispersionMechanism_Setup()
//  This macro will setup the input waves needed to generate the reaction mechanism
//  After running this macro, enter the reaction mechanism
//  Reactant and Product names need to be less than 11 characters in length
//  Do not use Reactant and Product as species names.
//
//  Enter the Rate coefficients for the reactions in the RateCoefficients wave
//
//  The time base is set for 60 sec.
//              This can be changed in the command window as desired later:  i.e.  caltime = 1e-3*exp(x/5) (unequally spaced points)
//
	Make/O/T/N=100 Reactant1, Reactant2, Product1, Product2, Product3, Names
	Reactant1 = ""
	Reactant2 = ""
	Product1 = ""
	Product2 = ""
	Product3 = ""
	Names = ""
	Make/O/N=100 RateCoefficients, BackgroundConc
	RateCoefficients=0
	BackgroundConc=0
	Make/O/N=2 CalTime=p*60	// initialize CalTime to 60 sec time step.
	
	
	Edit/W=(5.25,42.5,1122,236) reactant1,reactant2,product1,product2, product3, RateCoefficients		// make new table of input waves.
	DoWindow/K ReactionMechanism							// kill existing table if it exists
	DoWindow/C/T ReactionMechanism,"Reaction Mechanism"	// new table
	ModifyTable  size=11,width=90
	ModifyTable width(ratecoefficients)=100
	Variable /G dilution = 0
	IncludeChemKineticFuncs()									// load subroutines for reaction mechanism from experiment folder
end

Macro DispersionMechanismCompile(fctname)
	String fctname = "ChemKinetic"
	Prompt fctname, "Select Kinetic Function", popup FunctionList("ChemKinetic*",";","KIND:2")
	
	Variable/G NDensity = 1.865e19 //CHANGE WITH T AND P
	String/G CurrChemKineticFct = fctname
	Variable Rnumber, SpeciesNumber, SpeciesDiff
	Variable flag,n,m, nn, position, position2
	Variable na,nb,nc,nd, ne
	Variable TotalLen, breakpos
	String temp, temp1, temp2, temp3, temp4, temp5
	String  TotalTemp, totaltemp1, totaltemp2, totaltemp3, totaltemp4,  totaltemp5
	String dilutiontemp
	String sa,sb,sc,sd, se
	String s0,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10
	String Plus, Equal
	String MyList, tempNote, tempNote2
	String/G dydtstr = "" // string to be passed to ChemKineticsFunction
	Silent 1; PauseUpdate
	Plus = "+"
	Equal = "="

	s0 = ""
	s1 = " "
	s2 = S1 + " "
	s3 = s2 + " "
	s4 = s3 + " "
	s5 = s4 + " "
	s6 = s5 + " "
	s7 = s6 + " "
	s8 = s7 + " "
	s9 = s8 + " "
	s10 = s9 + " "

	n=0
	//printf "\r\t\t\t\tReaction Mechanism \t\t\t\t\t\t\t\t\t\t Rate Coefficient\r"
	
	Do	// loop 0-100 over reactants and products to find number of reactions. Break on empty row for each of the 5 waves.

		temp = Reactant1[n] + Reactant2[n] + Product1[n] + Product2[n] + Product3[n]
		if (stringmatch(temp,"")) //Empty line in rxn. mechanism means jump to the end of the loop
			Rnumber = n  // This is the number of reactions
			break

		else					//There is at least one reactant or product
        		na = 10 - strlen(Reactant1[n])
        		temp = "s" + num2str(na)
        		sa =($temp) + Reactant1[n] + plus
        		if(strlen(Reactant2[n]) == 0)
        			sa =($temp) + Reactant1[n]
        		endif
         		if (strlen(Reactant1[n]) == 0)
         			sa = "                   "
        		endif


		      nb = 10 - strlen(Reactant2[n])
        		temp = "s" + num2str(nb)
        		sb =($temp) + Reactant2[n] + equal
         		if (strlen(Reactant2[n]) == 0)
         			sb = "                     " + equal

        		endif

        		nc = 10 - strlen(Product1[n])
        		temp = "s" + num2str(nc)
        		sc =($temp) + Product1[n]
        		if (strlen(Product1[n]) == 0)
         			sc = "                     "
        		endif

        		nd = 10 - strlen(Product2[n])
        		temp = "s" + num2str(nd)
        		sd = plus + ($temp) + Product2[n]
         		if (strlen(Product2[n]) == 0)
         			sd = "                  "
        		endif

 			ne = 10 - strlen(Product3[n])
        		temp = "s" + num2str(ne)
        		se = plus + ($temp) + Product3[n]
         		if (strlen(Product3[n]) == 0)
         			sd = "                  "
        		endif
 			
			//print "\t\t\t\t" + sa  + sb + sc + sd+ se + "\t\t\t\t\t" + num2str(ratecoefficients[n])
		endif

  		n += 1
 	while (n < 100)
 //	print "\r\r"
 //	Print " ODEs in Compound Format"


//  Transfer Rate Coefficient Data to KK wave
	Make/O/N=(Rnumber) KK
//   n = 0
//   Do
//   	kk[n] = RateCoefficients[n]
//		n +=1
//		while(n < Rnumber)
	KK = RateCoefficients

// Capture species names and write to Names[n] wave

	n = 0
	Make/O/N=600/T Names = ""
	Do 
		nn = n*5
		Names[nn] = Reactant1[n]
		Names[nn+1] = Reactant2[n]
		Names[nn+2] = Product1[n]
		Names[nn+3] = Product2[n]
		Names[nn+4] = Product3[n]

		n += 1
	while (n < rnumber)

//  Find number of species in Reaction Names List Remove Blanks from Names List

	n = FindLastNonEmtpyRowInTextwave(Names)
	if ((n!=-1) && (n<numpnts(Names)-1))
		DeletePoints n+1, 600, Names	// delete points from last non-blank row to end
	endif
	Do									// loop until there are no blank rows found
		FindValue/TEXT="" Names		// finds first blank entry in Names
		if (V_value==-1)				// No blank cells found.
			SpeciesNumber = numpnts(Names)
		else
			DeletePoints V_value, 1, Names
		endif
	While(V_value!=-1)


//  Remove Duplicate Names from Names List
	removeDuplicatesInTextWave("Names")
	SpeciesNumber = numpnts(Names)		// redefine in case there were duplicates, and Names is shorter now.
	                

	Make/O/N=(SpeciesNumber) InitialConcentrations

//  Make Output waves for Species with wave Names from Names[n]

	n =0
	Do
		Make/O/D/N=2 $Names[n] = 0
		n +=1
	while (n < SpeciesNumber)

	DoWindow/F ReactionMechanism
	if (V_flag==0)
		Abort "Table not found."	// Make option to create it?
	endif

//	Append more columns to table if they aren't present already.
	string wlist = WaveList("*", ";", "WIN:")
	if (WhichListItem("KK", wlist)==-1)
		AppendToTable KK
		ModifyTable size(KK)=14,rgb(KK)=(65535,0,0), width(KK) = 100
	endif
	if (WhichListItem("Names", wlist)==-1)
		AppendToTable Names
		ModifyTable size(Names)=14, width(Names)=100
	endif
	if (WhichListItem("BackgroundConc", wlist)==-1)
		AppendToTable BackgroundConc
		ModifyTable size(BackgroundConc)=14, rgb(BackgroundConc)=(65535,0,0), width(BackgroundConc)=150
	endif
	if (WhichListItem("InitialConcentrations", wlist)==-1)
		AppendToTable InitialConcentrations
		ModifyTable size(initialconcentrations)=14, rgb(InitialConcentrations)=(65535,0,0), width(initialconcentrations)=150
	endif
	if (WhichListItem("CalTime", wlist)==-1)
		AppendToTable CalTime
		ModifyTable size(CalTime)=14, rgb(CalTime)=(65535,0,0)
	endif


// Print  ODEs  for Reaction Mechanism in Chemical Format

	n = 0
	Do		//while (n < SpeciesNumber)
		m=0
		Totaltemp = ""
		Totaltemp1 = ""
		Totaltemp2 = ""
		Totaltemp3 = ""
		Totaltemp4 = ""
		Totaltemp5 = ""

		Do // Find Reactant 1 in Names
			temp1 = ""
			If(stringmatch(Names[n], Reactant1[m]))
				temp1 = "-" + num2str(ratecoefficients[m]) + " * [" + Reactant1[m] + "] * ["+Reactant2[m] + "]"

				if(strlen(Reactant1[m]) == 0)
					temp1 = "-" + num2str(ratecoefficients[m]) + " * " +  "["+Reactant2[m] + "]"
				endif

				if(strlen(Reactant2[m]) == 0)
					temp1 = "-" + num2str(ratecoefficients[m]) + " * [" + Reactant1[m] + "]"
				endif

			endif
			Totaltemp1 += temp1

			m += 1
		While(m <rnumber)

		m = 0
		Do // Find Reactant 2 in Names
			temp2 = ""
			If(stringmatch(Names[n], Reactant2[m]))
				temp2 = "-" + num2str(ratecoefficients[m]) + " * [" + Reactant1[m] + "] * ["+Reactant2[m] + "]"
	
				if(strlen(Reactant1[m]) == 0)
					temp2 = "-" + num2str(ratecoefficients[m]) + " * " +  "["+Reactant2[m] + "]"
				endif
				
				if(strlen(Reactant2[m]) == 0)
					temp2 = "-" + num2str(ratecoefficients[m]) + " * [" + Reactant1[m] + "]"
				endif
	
			endif
		
			Totaltemp2 += temp2
	
			m += 1
		While(m <rnumber)


		m = 0
		Do	// Find Product 1 in Names
			temp3 = ""
			If(stringmatch(Names[n], Product1[m]))
				temp3 = "+" + num2str(ratecoefficients[m]) + " * [" + Reactant1[m] + "] * ["+Reactant2[m] + "]"

				if(strlen(Reactant1[m]) == 0)
					temp3 = "+" + num2str(ratecoefficients[m]) + " * " +  "["+Reactant2[m] + "]"
				endif
				
				if(strlen(Reactant2[m]) == 0)
					temp3 = "+" + num2str(ratecoefficients[m]) + " * [" + Reactant1[m] + "]"
				endif

			endif
			Totaltemp3 += temp3
	
			m += 1
		While(m <rnumber)

		m = 0
		Do		// Find Product 2 in Names
			temp4 = ""
			If(stringmatch(Names[n], Product2[m]))
				temp4 = "+" + num2str(ratecoefficients[m]) + " * [" + Reactant1[m] + "] * ["+Reactant2[m] + "]"
		
				if(strlen(Reactant1[m]) == 0)
					temp4 = "+" + num2str(ratecoefficients[m]) + " * " +  "["+Reactant2[m] + "]"
				endif
				if(strlen(Reactant2[m]) == 0)
					temp4 = "+" + num2str(ratecoefficients[m]) + " * [" + Reactant1[m] + "]"
				endif
			endif
			Totaltemp4 += temp4

			m += 1
		While(m <rnumber)

		m = 0
		Do		// Find Product 3 in Names
			temp5 = ""
			If(stringmatch(Names[n], Product3[m]))
				temp5 = "+" + num2str(ratecoefficients[m]) + " * [" + Reactant1[m] + "] * ["+Reactant2[m] + "]"
		
				if(strlen(Reactant1[m]) == 0)
					temp5 = "+" + num2str(ratecoefficients[m]) + " * " +  "["+Reactant2[m] + "]"
				endif
				if(strlen(Reactant2[m]) == 0)
					temp5 = "+" + num2str(ratecoefficients[m]) + " * [" + Reactant1[m] + "]"
				endif
			endif
			Totaltemp5 += temp5

			m += 1
		While(m <rnumber)


		Totaltemp = totaltemp1 + totaltemp2 + totaltemp3 + totaltemp4 + totaltemp5 // combine all instances where the name was found

		if(stringmatch(totaltemp, ""))

		else
			//Print "d["+Names[n] + "]/dt = " + Totaltemp
		endif


		n +=1
	while (n < SpeciesNumber)

// Print  ODEs  for Reaction Mechanism in ChemKinetic (integrateODE) Format

//*******************************************
//  Make Notebook to contain ODEs to be Cut and Pasted

//                                Mylist = WinList("ReactionCopy",";","WIN:16")
//                                if (strlen(MyList) == 0)                
//                                        else                    
//                                        DoWindow/K ReactionCopy                                                         
//                                endif
//                                        NewNotebook/N=ReactionCopy/F=1/V=1/K=0/W=(5,42,505,337) 
//                                        Notebook ReactionCopy defaultTab=36, statusWidth=238, backRGB=(65535,0,0), pageMargins={72,72,72,72}
//                                        Notebook ReactionCopy showRuler=1, rulerUnits=1, updating={1, 60}
//                                        Notebook ReactionCopy newRuler=Normal, justification=0, margins={0,0,468}, spacing={0,0,0}, tabs={}, rulerDefaults={"Times",20,0,(0,0,0)}
//                                        
//                                        Notebook ReactionCopy ruler=Normal, text = " Cut and Paste these ODEs into ChemKinetic Function\r\r"
//***********************************************
	n = 0
	Do		// while (n < SpeciesNumber)
		m=0
		Totaltemp = ""
		Totaltemp1 = ""
		Totaltemp2 = ""
		Totaltemp3 = ""
		Totaltemp4 = ""
		Totaltemp5 = ""
		Do
			temp1 = ""
			If(stringmatch(Names[n], Reactant1[m]))
				nn = 0
				position = 0
				Do
					if(stringmatch(Names[nn],Reactant1[m]))
						position = nn
						break	//nn = 500
					endif
					nn +=1
//				while(nn < rnumber*4)
				while(nn < SpeciesNumber)	// KA try this to see if it fixes compile bug
		                        
				nn = 0
				position2 = 0
				Do
					if(stringmatch(Names[nn],Reactant2[m]))
						position2 = nn
						break	//nn = 500
					endif
					nn +=1
//				while(nn < rnumber*4)
				while(nn < SpeciesNumber)	// KA try this to see if it fixes compile bug

				temp1 = "-" + "pw[" + num2str(m) + "]*" + "yw["+ num2str(position) + "]*" + "yw["+ num2str(position2) + "]"
		
				if(strlen(Reactant1[m]) == 0)
					temp1 = "-" + "pw[" +   num2str(m) + "]*" + "yw["+num2str(position2) + "]"
				endif
				if(strlen(Reactant2[m]) == 0)
					temp1 = "-" + "pw[" +  num2str(m)+ " ]*"  + "yw[" + num2str(position) + "]"
				endif
			endif
			Totaltemp1 = Totaltemp1 + temp1

			m += 1
		While(m <rnumber)

		m = 0
		Do
			temp2 = ""
			If(stringmatch(Names[n], Reactant2[m]))

				nn = 0
				position = 0
				Do
					if(stringmatch(Names[nn],Reactant1[m]))
						position = nn
						break	// nn = 500
					endif
					nn +=1
//				while(nn < rnumber*4)
				while(nn < SpeciesNumber)	// KA try this to see if it fixes compile bug

				nn = 0
				position2 = 0
				Do
					if(stringmatch(Names[nn],Reactant2[m]))
						position2 = nn
						break	// nn = 500
					endif
					nn +=1
//				while(nn < rnumber*4)
				while(nn < SpeciesNumber)	// KA try this to see if it fixes compile bug


				temp2 = "-" + "pw[" + num2str(m) + "]*" + "yw["+ num2str(position) + "]*" + "yw["+ num2str(position2) + "]"

				if(strlen(Reactant1[m]) == 0)
					temp2 = "-" + "pw[" +   num2str(m) + "]*" + "yw["+num2str(position2) + "]"
				endif
				if(strlen(Reactant2[m]) == 0)
					temp2 = "-" + "pw[" +  num2str(m)+ " ]*"  + "yw[" + num2str(position) + "]"
				endif

			endif
			Totaltemp2 = Totaltemp2 + temp2

			m += 1
		While(m <rnumber)


		m = 0
		Do
			temp3 = ""
			If(stringmatch(Names[n], Product1[m]))

				nn = 0
				position = 0
				Do
					if(stringmatch(Names[nn],Reactant1[m]))
						position = nn
						break	// nn = 500
					endif
					nn +=1
//				while(nn < rnumber*4)
				while(nn < SpeciesNumber)	// KA try this to see if it fixes compile bug

				nn = 0
				position2 = 0
				Do
					if(stringmatch(Names[nn],Reactant2[m]))
						position2 = nn
						break	// nn = 500
					endif
					nn +=1
//				while(nn < rnumber*4)
				while(nn < SpeciesNumber)	// KA try this to see if it fixes compile bug


				temp3 = "+" + "pw[" + num2str(m) + "]*" + "yw["+ num2str(position) + "]*" + "yw["+ num2str(position2) + "]"

				if(strlen(Reactant1[m]) == 0)
					temp3 = "+" + "pw[" +   num2str(m) + "]*" + "yw["+num2str(position2) + "]"
				endif
				if(strlen(Reactant2[m]) == 0)
					temp3 = "+" + "pw[" +  num2str(m)+ " ]*"  + "yw[" + num2str(position) + "]"
				endif

			endif
			Totaltemp3 = Totaltemp3 + temp3

			m += 1
		While(m <rnumber)

		m = 0
		Do
			temp4 = ""
			If(stringmatch(Names[n], Product2[m]))

				nn = 0
				position = 0
				Do
					if(stringmatch(Names[nn],Reactant1[m]))
						position = nn
						break	// nn = 500
					endif
					nn +=1
//				while(nn < rnumber*4)
				while(nn < SpeciesNumber)	// KA try this to see if it fixes compile bug

				nn = 0
				position2 = 0
				Do
					if(stringmatch(Names[nn],Reactant2[m]))
						position2 = nn
						break	// nn = 500
					endif
					nn +=1
//				while(nn < rnumber*4)
				while(nn < SpeciesNumber)	// KA try this to see if it fixes compile bug


				temp4 = "+" + "pw[" + num2str(m) + "]*" + "yw["+ num2str(position) + "]*" + "yw["+ num2str(position2) + "]"

				if(strlen(Reactant1[m]) == 0)
					temp4 = "+" + "pw[" +   num2str(m) + "]*" + "yw["+num2str(position2) + "]"
				endif
				if(strlen(Reactant2[m]) == 0)
					temp4 = "+" + "pw[" +  num2str(m)+ " ]*"  + "yw[" + num2str(position) + "]"
				endif

			endif
			Totaltemp4 = Totaltemp4 + temp4

			m += 1
		While(m <rnumber)

		m = 0

		Do
			temp5 = ""
			If(stringmatch(Names[n], Product3[m]))

				nn = 0
				position = 0
				Do
					if(stringmatch(Names[nn],Reactant1[m]))
						position = nn
						break	//nn = 500
					endif
					 nn +=1
//				while(nn < rnumber*4)
				while(nn < SpeciesNumber)	// KA try this to see if it fixes compile bug

				nn = 0
				position2 = 0
				Do
					if(stringmatch(Names[nn],Reactant2[m]))
						position2 = nn
						break	// nn = 500
					endif
					nn +=1
//				while(nn < rnumber*4)
				while(nn < SpeciesNumber)	// KA try this to see if it fixes compile bug


				temp5 = "+" + "pw[" + num2str(m) + "]*" + "yw["+ num2str(position) + "]*" + "yw["+ num2str(position2) + "]"

				if(strlen(Reactant1[m]) == 0)
					temp5 = "+" + "pw[" +   num2str(m) + "]*" + "yw["+num2str(position2) + "]"
				endif
				if(strlen(Reactant2[m]) == 0)
					temp5 = "+" + "pw[" +  num2str(m)+ " ]*"  + "yw[" + num2str(position) + "]"
				endif

			endif
			Totaltemp5 = Totaltemp5 + temp5

			m += 1
		While(m <rnumber)


		Totaltemp = totaltemp1 + totaltemp2 + totaltemp3 + totaltemp4 + totaltemp5

// include dilution
		if (exists("dilution"))
			dilutiontemp = "-d*yw["+ num2str(n) + "]"
			Totaltemp += dilutiontemp
		endif

		TotalLen = StrLen(Totaltemp)        

		if (TotalLen>0)
                	
			if (TotalLen>390) // string too long for one line in function
				breakpos = StrSearch(Totaltemp, "+", 300) // find first "+" past 300 chars
				if (breakpos<0 || breakpos> 390)
					breakpos = StrSearch(Totaltemp, "-", 300)) // Check if there's maybe a "-"
				endif
				if (breakpos>300)
					//Print "dydt["+ num2str(n) + "]=" + Totaltemp[0,breakpos-1]
					//Print "dydt["+ num2str(n) + "]+=" + Totaltemp[breakpos, StrLen(Totaltemp)-1]
					dydtstr += "dydt["+ num2str(n) + "]="+Totaltemp[0,breakpos-1] + "\r" 
					dydtstr += "dydt["+ num2str(n) + "]+="+ Totaltemp[breakpos, StrLen(Totaltemp)-1] + "\r" 
				else
					//print Totaltemp
					abort "something is strange with this string"
				endif
					
			else
				//Print "dydt["+ num2str(n) + "]=" + Totaltemp
				dydtstr += "dydt["+ num2str(n) + "]=" + Totaltemp+ "\r" 
			endif

		endif

		n +=1
	while (n < SpeciesNumber)

	RemoveChemKineticFuncs()
	Execute/Q /P "EditChemKineticFuncs()"
	IncludeChemKineticFuncs()

end

// Keep for now to refer to original setup for model without dispersion.
//Macro Run_Model2_ORIG()	// ORIGINAL CODE
//
//	Silent 1
//	
//	If (!exists("CurrChemKineticFct"))
//		Abort "Must compile ODE first!"
//	Endif
//	String cmd, SpeciesList, CurrSpecies
//	Variable n, speciesnumber
//	Variable npnts
//	SetDataFolder root:
//	Silent 1
//	PauseUpdate
//	Variable starttime = DateTime
//	n=0
//	SpeciesList = ""
//	SpeciesNumber = numpnts(Names)
//	Do
//		if(n == 0)
//			SpeciesList += Names[n]
//		else
//			SpeciesList = SpeciesList + "," + Names[n]
//		endif
//		n += 1
//	while (n < numpnts(Names))
//	
//	n=0
//	Do
//		CurrSpecies = Names[n]
//		print $CurrSpecies
//		
//		$CurrSpecies[0] = InitialConcentrations[n] // restart: initialize to original values
//		n += 1
//	while(n < SpeciesNumber)
//
//	cmd= "IntegrateODE/U=1000/X=CalTime "+CurrChemKineticFct+" KK, {"+ SpeciesList + "}"
//	//Print cmd
//	Execute cmd
//	//print "Calculation time", DateTime - starttime,"seconds"
//
//end

Macro Transect(p3)
	variable p3
	
	make/N=201/D/O NO_transect, NO2_transect, O3_transect, NO3_transect, N2O5_transect, HNO3_transect, ClNO2_transect, MonNO3_transect, AlkNO3_transect, PAN_transect, Ace_transect, NOy_transect, OrgN_transect
	SetScale/I x -10,10,"", NO_transect, NO2_transect, O3_transect, NO3_transect, N2O5_transect, HNO3_transect, ClNO2_transect, MonNO3_transect, AlkNO3_transect, PAN_transect, Ace_transect, NOy_transect, OrgN_transect
	
	NO_transect = NO_mat[p3][p]
	NO2_transect = NO2_mat[p3][p]
	O3_transect = O3_mat[p3][p]
	HNO3_Transect = HNO3_mat[p3][p]
	N2O5_Transect = N2O5_mat[p3][p]
	NO3_Transect = NO3_mat[p3][p]
	ClNO2_Transect = ClNO2_mat[p3][p]
	MonNO3_transect = MonNO3_mat[p3][p]
	AlkNO3_transect = AlkNO3_mat[p3][p]
	PAN_transect = PAN_mat[p3][p]
	Ace_transect = Ace_mat[p3][p]
	NOy_transect = NO_transect + NO2_transect + HNO3_transect + 2*N2O5_transect + NO3_transect + ClNO2_transect + MonNO3_transect + AlkNO3_transect + PAN_transect
	OrgN_transect = MonNO3_transect + AlkNO3_transect
	//print p3
End Macro

function removeDuplicatesInTextWave(textwaveStr)
	string textwaveStr
	
	variable n, len, start
	string species
	wave/T tw = $textwaveStr
	
	len = numpnts(tw)
	
	for (n=0; n<len; n+=1)
		start = 0
		species = tw[n]
		do
			FindValue/TEXT=species/TXOP=2/S=(start) tw
			if ((V_value!=-1) && (V_value!=n))
				DeletePoints V_value, 1, tw
				len-=1
			else
				start = v_value+1	// found species, but it may have been itself (n==V_value), so increment to see if there is another one.
			endif
		while((v_value!=-1) && (start<len))
	endfor
end

// Return the last row of a text wave that contains text. Returns -1 if whole wave is empty.
function FindLastNonEmtpyRowInTextwave(twave)
	wave/T twave
	
	variable, n, len, foundIt
	len = numpnts(twave)
	for(n=len-1; n>=0; n-=1)
		if (cmpstr(twave[n], "")!=0)
			foundIt = 1
			break
		endif
	endfor
	
	if (foundIt)
		return n
	else
		return -1
	endif
end