#pragma rtGlobals=1		// Use modern global access method.
#include <Wave Loading>
// This procedure analyzes multiple data sets with 2D gaussian fits
Function Analyze2DGauss(wlow,whigh,num)
	Variable wlow
	Variable whigh
	Variable num // number of points per data set 
	
	CloseAllGraphs()
	KillWaves/A/Z
	String theItem
	String wvnm
	String scnnm
	String clnm
	String dly
	String gaussfit
	String slopenm = "slope"
	String gfit_coefs
	Variable T
	Variable j //for organizing data sets
	Variable index=0
	Variable index1=0
	Make/N=0 Delay, Gausscor //Makes T delay wave and Gausscor wave; we will insert another point for every 2D IR spectrum
	Silent 1; PauseUpdate
	NewPath/O dataFolder //dialog box to get the data folder path
	String files=ListFilesOfType("dataFolder", 1, ".dat")
	String callist=ListFilesofType("dataFolder",1,".calib")
	String calibfile=StringFromList(0, callist)
	LoadWave/G/D/O/N=calib/P=datafolder calibfile
	Axes(calib0)
	do
		theItem= StringFromList(index, files)
		if (strlen(theItem) == 0)
			break								// ran out of items
		endif
		//Print(theItem)
		LoadWave/G/M/D/A=wave/P=datafolder theItem
		sscanf theItem, "%4s%d", wvnm, T
		wvnm += num2str(T)
		rename wave0 $wvnm
		if(index==0)
			Chop($wvnm,w1,w3,whigh,wlow)
		endif
		if(index>0)
			Chop2($wvnm)
		endif		
		
		Display2D(wvnm)  //This displays a contour plot for each loaded spectrum.  Tab out if the data set is large to avoid crashes.
		InsertPoints index, 1, Delay, Gausscor
	
	///////////////////////////////////////////////////////	
		//Gaussianfit($wvnm, w1, w3)
		//scnnm= wvnm+"gaussfit"
		//rename g_fit gaussfit
		//AppendMatrixContour $gaussfit vs {w1,w3}
		//ModifyContour $gaussfit labels=0, rgbLines=(0,0,0), autoLevels={*,0,11}
		//Wave coefs=goefs
		//Gausscor[index]=coefs[6]
		
		Gausscor[index] =Gaussianfit($wvnm,w1,w3)
		gaussfit=wvnm+"gaussfit"
		   rename g_fit $gaussfit
		   AppendMatrixContour $gaussfit vs {w1,w3}
		   ModifyContour $gaussfit labels=0, rgbLines=(0,0,0), autoLevels={*,0,11}
		
		Delay[index] = T
		//Display2D(wvnm)  //This displays a contour plot for each loaded spectrum.  Tab out if the data set is large to avoid crashes.
		index += 1
	while (1)		// loop until break above

//For multiple data sets with "num" points per data set:
		for (j=0; j<=(numpnts(Gausscor)/num)-1; j+=1)//Sort the corr waves in to individual data sets  (num is input for # of points per data set)
			if (j<((numpnts(Gausscor)/num)-1))
				scnnm ="corr"+num2str(j+1)
				dly = "delay"+ num2str(j+1)
				Duplicate/R=[0+num*j,(num+(j*num)-1)] Gausscor, $scnnm		
				Duplicate/R=[0+num*j,(num+(j*num)-1)] delay, $dly // in order to sort properly, we must sort each peice of Gausscor with its corresponding time delay peice
				sort $dly $dly, $scnnm //sorts the data by delaytime
				killwaves $dly
			Else
				scnnm ="corr"+num2str(j+1)
				dly = "delay"+ num2str(j+1)
				Duplicate/R=[0+num*j,(num+(j*num)-1)] Gausscor, $scnnm
				Duplicate/R=[0+num*j,(num+(j*num)-1)] delay, delaysort // in order to sort properly, we must sort each cls peice of cls with its corresponding time delay peice
				sort delaysort delaysort, $scnnm //sorts the data by delaytime
				Endif
			endfor
	
//	Duplicate/R=[0,num-1] delay delaysort
//	killwaves delay
//	sort delaysort delaysort, $scnnm//cannot sort original delay wave because we need to preserve its ordering to sort the slope data at each frequency
		
	for (j=0; j<=(numpnts(Gausscor)/num)-1; j+=1)//Sort the corr waves in to individual data sets  (num is input for # of points per data set)
			scnnm ="corr"+num2str(j+1)
			if (j==0)
				Display $scnnm vs delaysort
			else 
				appendtograph $scnnm vs delaysort
			Endif
			endfor

	ModifyGraph mode=3,marker=19,msize=3,rgb($scnnm)=(0,0,0)
	killwaves Delay //we no longer need the unsorted 2D data and delay wave now that all data is sorted
	Execute "TileWindows/O=1/C"
End

Function ftircls(coef, x): FitFunc
Wave coef				//c0 = offset, c1 = absorbance amplitude, c2 = center freq, c3 = T2*, c4 = Delta magnitude, C5 = rel amp 1, c6 = tau1, c7 = rel amp 2, c8 = tau2, c9 = rel amp static, c10 = lifetime
Variable x

Variable i
Variable S1
S1 = 0
for(i=0; i<8192; i+=2)  //Doing a brute force FT
	S1 += Cos(2*pi*x*.00003*i)*Cos(2*pi*coef[2]*0.00003*i)*exp(-(i/coef[3]))*exp(-coef[4]*coef[4]*((coef[5]/(coef[5]+coef[7]+coef[9]))*(coef[6]*i + coef[6]*coef[6]*(exp(-i/coef[6])-1)) + (coef[7]/(coef[5]+coef[7]+coef[9]))*(coef[8]*i + coef[8]*coef[8]*(exp(-i/coef[8])-1)) + (coef[9]/(coef[5]+coef[7]+coef[9]))*i*i*0.5)) * exp(-i/(2*coef[10]))
endfor

S1 = coef[0] + coef[1]*S1

return S1

end

Function Axes(n)
wave n
Make/N=512/D/O w1
Make/N=1024/D/O w3
w1=n[5]+x/(3e-5*1024*n[6])
w3=10e6*(1/(n[0]*x^2+n[1]*x+n[2])-1/(n[4]))
end


Function CLS(data, w1, w3, n, lim4)
	Wave data
	Wave w1
	Wave w3
	Wave n
	Wave lim4
	variable  i
	Variable l= Dimsize(data,0)
	Variable x0=w1[0]
	Variable dx=w1[1]-w1[0]
	Variable w=Dimsize(data,1)

	Make/O/N=(w)  slice
	Make/O/N=(l) centline
	SetScale/p x  x0, dx,"", centline
	Make/O/N=(l) diagamp
	SetScale/p x  x0, dx,"", diagamp
	WaveStats/Q data
	for(i=0;i<l;i+=1)
		slice = data[i][x]
		//Smooth 40,slice  //This smooths the w3 slice that is fit to the CLS.  Tab out to disable smoothing.
		FindPeak/Q/N/M=(V_min/4)/B=15 slice 
	    //centline[i] =w3[V_peakloc]
	    centline[i]=10e6*(1/(n[0]*(V_peakloc+lim4)^2+n[1]*(V_peakloc+lim4)+n[2])-1/(n[4]))
	    diagamp[i] =V_PeakVal
	endfor
	KillWaves slice
end

Function Chop(fulldata, w1,w3,maximum,minimum)
wave fulldata
wave w1
wave w3
variable maximum
variable minimum
Variable l= Dimsize(fulldata,0)
Variable w=Dimsize(fulldata,1)
variable a=0
variable b=0
variable c=0
variable d=0
variable restw1=l-a
variable restw3=w-c
do
	a+=1
while(w1[a]<maximum)
do
	b+=1
while(w1[b]<minimum)
do
	c+=1
while(w3[c]<maximum)
do
	d+=1
while(w3[d]<minimum)
deletepoints/M=0 a,restw1,fulldata
deletepoints/M=1 c,restw3,fulldata
deletepoints/M=0 0,b,fulldata
deletepoints/M=1 0,d,fulldata
deletepoints/M=0 a,restw1,w1
deletepoints/M=0 c,restw3,w3
deletepoints/M=0 0,b,w1
deletepoints/M=0 0,d,w3
Make/N=1/D/O lim1; lim1=a
Make/N=1/D/O lim2; lim2=b
Make/N=1/D/O lim3; lim3=c
Make/N=1/D/O lim4; lim4=d
End

Function Chop2(n)
wave n
Variable l= Dimsize(n,0)
Variable w=Dimsize(n,1)
wave lim1
wave lim2
wave lim3
wave lim4
variable restw1=l-lim1[0]
variable restw3=w-lim3[0]
deletepoints/M=0 lim1[0],restw1,n
deletepoints/M=1 lim3[0],restw3,n
deletepoints/M=0 0,lim2[0],n
deletepoints/M=1 0,lim4[0],n
End

Function Display2D(corr)
	String corr
	
	String corr_cls
	String fit_corr_cls
	Display;AppendMatrixContour $corr vs {w1,w3}
	SetAxis/A
	ModifyContour $corr labels=0
	//corr_cls = corr+"cls"
	//AppendToGraph $corr_cls
	//ModifyGraph mode($corr_cls)=3,marker($corr_cls)=8,msize($corr_cls)=3;DelayUpdate
	//ModifyGraph rgb($corr_cls)=(0,0,65535)
	//fit_corr_cls = "fit_" + corr_cls
	//AppendToGraph $fit_corr_cls
End	

Function CloseAllGraphs()
String name
do
name = WinName(0,1) // name of the front graph
if (strlen(name) == 0)
break // all done
endif
DoWindow/K $name // Close the graph
while(1)
End

ss
Function Gaussianfit(data, w1, w3)
	Wave data
	Wave w1
	Wave w3
	duplicate/O data g_fit
	//string gaussfit
	Make/O/D gcoefs={0, -0.05, 2049.3, 3.5, 2047.2, 3.5, 0.3}
	Make/O/T/N=6 T_Constraints 
	T_Constraints[0] = {"K1< 0","K2 > 0","K3 > 0","K4 > 0","K5 > 0","K6 > 0","K6<1"}
	//CurveFit/NTHR=0/G Gauss2D, kwCWave=gcoefs, data /X=w1 /Y=w3 /D=g_fit/C=T_Constraints//Original code, fits to entire data, works fine
	CurveFit/G/NTHR=0 Gauss2D kwCWave=gcoefs,  data[33,41][127,152] /X=w1 /Y=w3 /C=T_Constraints//New line, should fit only a small region of the spectrum 
	Return gcoefs[6]
End