#pragma rtGlobals=1		// Use modern global access method.
#include <Scatter Plot Matrix 2>
#include <All Gizmo Procedures>
#include <Autosize Images>
#include <Wave Arithmetic Panel>
#include <Multi-peak fitting 2.0>

// written by Matt Sfeir, Brookhaven National Laboratory
// msfeir@bnl.gov
// v.3 February 11 2011


Menu "Load Waves"
	"Load Helios TA Data...", ProcessTA()
End

Menu "TA Analysis"
	"Load Helios TA Data...", ProcessTA()
	"Extract Kinetics", ExtractKin()
	"Extract Spectra", ExtractSpec()
End

Function ProcessTA()  // CMD - . aborts the procedure
	Variable refNum
	String message = "Select a file" 
	String outputPath, fileHandle, wName, wName2
	String fileFilters = "Helios TA Data Files (*.csv):.csv;" 
	fileFilters += "All Files:.*;" 
	Open /D /R /F=fileFilters /M=message refNum 
	outputPath = S_fileName
	printf "%s\r" outputPath
	fileHandle = ParseFilePath(5, outputPath, ":", 0,0)
	fileHandle = ParseFilePath(3, outputPath, ":", 0,0)
	fileHandle = fileHandle
	printf "%s\r" fileHandle
	LoadWave/J/M/U={0,2,0,2}/D/A=$fileHandle/K=0 outputPath
	Wave w =outputPath
	// WaveExists  1 valid 0 invalid
	wName = fileHandle +"0"
	Variable b=0
	if (WaveExists ($wName) == 1)
		wName2 = fileHandle + num2str(b+1)
		if (WaveExists ($wName2) == 1)
			wName = wName2
			b += 1
		else
			wName = fileHandle + num2str(b)
		endif		
	endif
	Duplicate $wName w
	String RPName = "RP_" + wName
	String CPName = "CP_" + wName
	Wave RP =RPName
	Wave CP =CPName
	Duplicate $RPName RP
	Duplicate $CPName CP

//	Wave w,RP,CP	
	Variable m, n, p, q, r
	
// remove NaNs and Infs.  Replace with zeros.

	m = Dimsize(w,0) // numer of rows
	printf "%f\r" m
	n = Dimsize(w,1) // number of columns	
	
	Make/N=(m,n) /D/O Zmat
	Zmat= (numtype(w) == 1) ? 0  : w
	Zmat= (numtype(Zmat) == 2) ? 0  : Zmat
	
	Make/N=(m) /D/O ContX
	ContX= (numtype(RP) == 1) ? 0  : RP
	ContX= (numtype(ContX) == 2) ? 0  : ContX

// remove extra rows - often contain bad data

	q = Dimsize(ContX,0)
	
	if (q  > 512)
		Redimension /N=512 ContX
		Redimension /N=(512,-1) Zmat
	endif
	
	m = Dimsize(Zmat,0) // numer of rows
	n = Dimsize(Zmat,1) // number of columns
	p = m*n
			
	Make/N=(n) /D/O ContY
	ContY= (numtype(CP) == 1) ? 0  : CP
	ContY= (numtype(ContY) == 2) ? 0  : ContY

	Make/N=(p) /D/O tempXwave
	Make/N=(p) /D/O tempYwave
	Make/N=(p) /D/O tempZwave
	
	Make/N=(m+1) /D/O ImageX
	Make/N=(n+1) /D/O ImageY
	
	Make/N=(M,N,3) /O ParametricWave

// create new processed row, column positions and z matrix
// can be used to generate contour plot
	
	Variable i,j
	
	i=0
	do
		j=0
		do	
			tempXwave[j+m*i]=ContX[j]
			tempYwave[j+m*i]=ContY[floor((j+m*i)/m)]
			tempZwave[j+m*i]=Zmat[j][i]
			j +=1
		while (j <=m)
	i +=1	
	while (i <= n)
	
	Concatenate/O {tempXwave, tempYwave, tempZwave}, wScatterMatrix
	KillWaves tempXwave, tempYwave, tempZwave
	
// Format for Image Plot
// Create x and y wave for image.  Points define start and end of pixel positions, not center 
// so need to generate N+1 points from N.
// Add point at [-1] at same distance from [0] as [1]
// and shift by [0] - [1] distance OR pick midpoint between consecutive points
// assuming smallest increment at start of wave
	
	Variable g, h
	
	ImageX[0]= ContX[0]-((ContX[1]-ContX[0])/2)
//	g = 0
//	do
//		ImageX[g+1]=ContX[g]+((ContX[1]-ContX[0])/2)
//		g += 1
//	while (g <= m)
	g=1
	do
		ImageX[g] = (ContX[g] + ContX[g-1])/2
		g += 1
	while (g < m)
	ImageX[m] = ContX[m-1] + ((ContX[m-1] - ContX[m-2])/2)


	ImageY[0]= ContY[0]-((ContY[1]-ContY[0])/2)
//	h=0
//	do
//		ImageY[h+1]=ContY[h]+((ContY[1]-ContY[0])/2)
//		h += 1
//	while (h <= n)
	h=1
	do
		ImageY[h] = (ContY[h] + ContY[h-1])/2
		h+= 1
	while (h < n)
	ImageY[n] = ContY[n-1] + ((ContY[n-1] - ContY[n-2])/2)



// generate Parametic data.
// used for uneven spacings in surface plots

	Variable k,l
	
	k=0
	do
		l=0
		do
			ParametricWave[k][l][0]=ContX[k]
			ParametricWave[k][l][1]=ContY[l] 
			ParametricWave[k][l][2]=Zmat[k][l]
			l +=1
		while (l <= n)
	k +=1
	while (k <= m)

// clean up wave names
	
//	String ZmatName = NameOfWave(w) + "_zmat"
	String ZmatName = wName + "_zmat"
	Duplicate/O Zmat, $ZmatName
	KillWaves Zmat
	
//	String ContXName = NameOfWave(w) + "_xc"
	String ContXName = wName + "_xc"
	Duplicate/O ContX, $ContXName
	KillWaves ContX
	
//	String ContYName = NameOfWave(w) + "_yc"
	String ContYName = wName + "_yc"
	Duplicate/O ContY, $ContYName
	KillWaves ContY
	
//	String ImageXName = NameOfWave(w) + "_xim"
	String ImageXName = wName + "_xim"
	Duplicate/O ImageX, $ImageXName
	KillWaves ImageX
	
//	String ImageYName = NameOfWave(w) + "_yim"
	String ImageYName = wName + "_yim"
	Duplicate/O ImageY, $ImageYName
	KillWaves ImageY
	
//	String ScatMatName = NameOfWave(w) + "_sca"
	String ScatMatName = wName + "_sca"
	Duplicate/O wScatterMatrix, $ScatMatName
	KillWaves wScatterMatrix
	
//	String ParaName = NameOfWave(w) + "_par"
	String ParaName = wName + "_par"
	Duplicate/O ParametricWave, $ParaName
	KillWaves ParametricWave, w, RP, CP
	
	Display;AppendImage $ZmatName vs {$ImageXName,$ImageYName};DelayUpdate
	ModifyImage $ZmatName ctab= {*,*,Rainbow,0}
	ShowInfo
	
	GetWindow kwTopWin wtitle
	String imTitle = S_Value
	imTitle = RemoveEnding(imTitle, ":")
	Cursor/I /W=$imTitle A $ZmatName 1,1
	
	return 0
	
End

Function ExtractKin()
//	WAVE/Z w = CsrWaveRef(A)
//	if (!WaveExists(w))
//		return NaN
//	endif
//	String zmatName = NameOfWave(w)

	String info= ImageInfo("","",0)
	String zmatName = StringByKey("ZWAVE",info)
//	String xwName = StringByKey("XWAVE",info)
	String xwName = RemoveEnding(zmatName, "zmat")
	xwName = xwName + "xc"
	Variable xpos, lambda
	xpos = pcsr(A)
	lambda = hcsr(A)
	Prompt zmatName, "Enter matrix wave name: "
	Prompt xpos, "Enter cursor postion for desired wavelength: "
	Prompt xwName, "Enter wavelength wave name: "
	DoPrompt "Generate Kinetics", zmatName, xpos, xwName
	if (V_Flag)
		return -1
		endif
//	Print zmatName; Print xwName
//	Variable tName = NumberByKey("RECREATION",CsrInfo(A))
	KineticsPlot(zmatName,xpos,xwName)
	return 0
End

Function ExtractSpec()
//	WAVE/Z w = CsrWaveRef(A)
//	if (!WaveExists(w))
//		return NaN
//	endif
//	String zmatName = NameOfWave(w)

	String info= ImageInfo("","",0)
	String zmatName = StringByKey("ZWAVE",info)
	String ywName = RemoveEnding(zmatName, "zmat")
	ywName = ywName + "yc"
	Variable ypos, stime	
	ypos = qcsr(A)
	stime = vcsr(A)
	Prompt zmatName, "Enter matrix wave name: "
	Prompt ypos, "Enter cursor postion for desired spectra: "
	Prompt ywName, "Enter time wave name: "
	DoPrompt "Generate Spectra", zmatName, ypos, ywName
	if (V_Flag)
		return -1
		endif
//	Print zmatName; Print ywName
	SpectraPlot(zmatName,ypos,ywName)
	return 0
End

Function KineticsPlot(KMat, lambdaPos, lambdaWave)
	String KMat, lambdaWave
	Variable lambdaPos
	
	String tWave = RemoveEnding(KMat, "zmat") + "yc"
	Duplicate $KMat, kM
	Duplicate $lambdaWave, lM
	
	Variable m,n
	m = Dimsize(kM,0) // numer of rows
	n = Dimsize(kM,1) // number of columns

	String kWaveName
	Variable l = lM[lambdaPos]
	String lName = num2str(l)
//	printf "%s\r" lName
//	Print l; Print KMat
//	printf "%s\r" KMat
	kWaveName = RemoveEnding(KMat, "zmat") + "k" +lName
	// Start of added code
	kWaveName = CleanupName(kWaveName, 1)
//	Print kWaveName
	// End of added code
	Make/D/O/N=(n) $kWaveName
	Wave w = $kWaveName
	Variable i,j
	i=lambdaPos
	j=0
	do
		w[j] = kM[i][j]
		j+=1
	while (j<=n)
//	Duplicate $kWaveName w
//	KillWaves tempKWave
	KillWaves kM, lM
	
	Display w vs $tWave

return 0
end

Function SpectraPlot(SMat, timePos, timeWave)
	String SMat, timeWave
	Variable timePos
	
	String sWave = RemoveEnding(SMat, "zmat") + "xc"
	
	Wave sM = $SMat
	Wave tM = $timeWave
	
	Variable m,n
	m = Dimsize(sM,0) // numer of rows
	n = Dimsize(sM,1) // number of columns
	
	String tWaveName
	Variable tme = tM[timePos]
	String tmeName = num2str(tme)
	tWaveName = RemoveEnding(SMat, "zmat") + "s" +tmeName
	tWaveName = CleanupName(tWaveName, 1)
		
	Make/D/O/N=(m) $tWavename
	Wave v = $tWavename
	
	Variable i,j
	i=timePos
	j=0
	do
		v[j] = sM[j][i]
		j+=1
	while (j<=m)
	
	Display v vs $sWave
	
return 0
end

// Misc Fit Functions
// Keep t0 and pulsewidth identical for all decays
// G are rates

Function TwoExpGaussInf(w,x) : FitFunc
	Wave w
	Variable x

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ f(x) = y0+(amp1/G1)*ExpGauss(x-t0,G1,pw)+(amp2/G2)*ExpGauss(x-t0,G2,pw)+(Ainf/1E-9)*ExpGauss(x-t0,1E-9,pw)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 8
	//CurveFitDialog/ w[0] = y0
	//CurveFitDialog/ w[1] = amp1
	//CurveFitDialog/ w[2] = t0
	//CurveFitDialog/ w[3] = pw
	//CurveFitDialog/ w[4] = G1
	//CurveFitDialog/ w[5] = amp2
	//CurveFitDialog/ w[6] = G2
	//CurveFitDialog/ w[7] = Ainf

	return w[0]+(w[1]/w[4])*ExpGauss(x-w[2],w[4],w[3])+(w[5]/w[6])*ExpGauss(x-w[2],w[6],w[3])+(w[7]/1E-9)*ExpGauss(x-w[2],1E-9,w[3])
End

Function ThreeExpGaussMatt(w,x) : FitFunc
	Wave w
	Variable x

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ f(x) = y0+(amp1/G1)*ExpGauss(x-t0,G1,pw)+(amp2/G2)*ExpGauss(x-t0,G2,pw)+(amp3/G3)*ExpGauss(x-t0,G3,pw)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 9
	//CurveFitDialog/ w[0] = y0
	//CurveFitDialog/ w[1] = amp1
	//CurveFitDialog/ w[2] = t0
	//CurveFitDialog/ w[3] = pw
	//CurveFitDialog/ w[4] = G1
	//CurveFitDialog/ w[5] = amp2
	//CurveFitDialog/ w[6] = G2
	//CurveFitDialog/ w[7] = amp3
	//CurveFitDialog/ w[8] = G3

	return w[0]+(w[1]/w[4])*ExpGauss(x-w[2],w[4],w[3])+(w[5]/w[6])*ExpGauss(x-w[2],w[6],w[3])+(w[7]/w[8])*ExpGauss(x-w[2],w[8],w[3])
End
Function ExpGaussDampedOsc(w,t) : FitFunc
	Wave w
	Variable t

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ f(t) = (A/G1)*ExpGauss(t-t0,G1,sig)+ (B/G2)*ExpGauss(t-t0,G2,sig)+(C/G3)*ExpGauss(t-t0,G3,sig)*(D+sin(w*t+phi))
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ t
	//CurveFitDialog/ Coefficients 11
	//CurveFitDialog/ w[0] = A
	//CurveFitDialog/ w[1] = t0
	//CurveFitDialog/ w[2] = sig
	//CurveFitDialog/ w[3] = G1
	//CurveFitDialog/ w[4] = B
	//CurveFitDialog/ w[5] = G2
	//CurveFitDialog/ w[6] = C
	//CurveFitDialog/ w[7] = G3
	//CurveFitDialog/ w[8] = D
	//CurveFitDialog/ w[9] = w
	//CurveFitDialog/ w[10] = phi

	return (w[0]/w[3])*ExpGauss(t-w[1],w[3],w[2])+ (w[4]/w[5])*ExpGauss(t-w[1],w[5],w[2])+(w[6]/w[7])*ExpGauss(t-w[1],w[7],w[2])*(w[8]+sin(w[9]*t+w[10]))
End

Function ExpGaussDampedOsc2(w,t) : FitFunc
	Wave w
	Variable t

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ f(t) = (A/G1)*ExpGauss(t-t0,G1,sig)+ (B/G2)*ExpGauss(t-t0,G2,sig)+(C/G3)*ExpGauss(t-t0,G3,sig)*(D+sin(w*t+phi))+(E/G4)*ExpGauss(t-t0,G4,sig)*(F+sin(w2*t+phi2))
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ t
	//CurveFitDialog/ Coefficients 16
	//CurveFitDialog/ w[0] = A
	//CurveFitDialog/ w[1] = t0
	//CurveFitDialog/ w[2] = sig
	//CurveFitDialog/ w[3] = G1
	//CurveFitDialog/ w[4] = B
	//CurveFitDialog/ w[5] = G2
	//CurveFitDialog/ w[6] = C
	//CurveFitDialog/ w[7] = G3
	//CurveFitDialog/ w[8] = D
	//CurveFitDialog/ w[9] = w
	//CurveFitDialog/ w[10] = phi
	//CurveFitDialog/ w[11] = E
	//CurveFitDialog/ w[12] = G4
	//CurveFitDialog/ w[13] = F
	//CurveFitDialog/ w[14] = w2
	//CurveFitDialog/ w[15] = phi2

	return (w[0]/w[3])*ExpGauss(t-w[1],w[3],w[2])+ (w[4]/w[5])*ExpGauss(t-w[1],w[5],w[2])+(w[6]/w[7])*ExpGauss(t-w[1],w[7],w[2])*(w[8]+sin(w[9]*t+w[10]))+(w[11]/w[12])*ExpGauss(t-w[1],w[12],w[2])*(w[13]+sin(w[14]*t+w[15]))
End

Function ExpGaussDampedOsc3(w,t) : FitFunc
	Wave w
	Variable t

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ f(t) = (A/G1)*ExpGauss(t-t0,G1,sig)+ (C/G3)*ExpGauss(t-t0,G3,sig)*(D+sin(w*t+phi))
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ t
	//CurveFitDialog/ Coefficients 9
	//CurveFitDialog/ w[0] = A
	//CurveFitDialog/ w[1] = t0
	//CurveFitDialog/ w[2] = sig
	//CurveFitDialog/ w[3] = G1
	//CurveFitDialog/ w[4] = C
	//CurveFitDialog/ w[5] = G3
	//CurveFitDialog/ w[6] = D
	//CurveFitDialog/ w[7] = w
	//CurveFitDialog/ w[8] = phi

	return (w[0]/w[3])*ExpGauss(t-w[1],w[3],w[2])+ (w[4]/w[5])*ExpGauss(t-w[1],w[5],w[2])*(w[6]+sin(w[7]*t+w[8]))
End