#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.

// by tony.withers@uwo.ca
// Version 1.0, 31/1/17

// Autocorrelation anlysis, 
// following Salje et al. (2000) EJM 12:503-519

// Plot a baseline-corrected, x-scaled data wave and select a region 
// for analysis by dragging mouse to make a marquee. Right click
// within marquee, select "Autocorrelation Analysis" and click on the
// name of the data wave.
// If the calculation is successful, a plot of fitted Gaussian widths
// (parameter K2) as a function of offset is created.
// Drag the red bars to select a region within which to fit a parabolic
// function.
// The intercept of the fit provides an estimate of deltaCorr.

Menu "GraphMarquee", dynamic
	submenu "Autocorrelation Analysis"
		TraceNameList("", ";",1), AutoCorrelationAnalysis()
	end
End

// use marquee to select region for analysis
// might be slow if you have densely sampled data
function AutoCorrelationAnalysis()
	GetLastUserMenuInfo // S_value will be trace name
	string tracename=S_value
	wave w=TraceNameToWaveRef("", tracename )
	
	// assume wave is correctly scaled
	string StrXAxisName=stringbykey("XAXIS", TraceInfo("", tracename, 0 ))
	getmarquee /K $StrXAxisName
	
	duplicate /free /R=(V_left,V_right) w w_data
	
	// initialize (or reinitialize) the package data folder
	NewDataFolder /O root:Packages
	NewDataFolder /O root:Packages:CorrelationAnalysis
	
	// Create a data folder reference variable
	DFREF dfr = root:Packages:CorrelationAnalysis
	
	make /o /n=10001 dfr:Correlation /WAVE=correlation
	setscale /i x, -500, 500,Correlation
	correlation = correlationStep(w_data, pnt2x(correlation, p)) 
	
	make /o /n=1 dfr:intercept=nan, dfr:interceptSD=nan
	
	// do many fits over different ranges and extract K2 fit parameter values
	FitToGetK2(correlation)	
	makeK2fitPlot()
	
end

function correlationStep(w, step)
	wave w
	variable step
	variable i, lim=numpnts(w), v_corr=0
	variable var_x
	for(i=0;i<lim;i+=1)
		var_x=pnt2x(w, i)+step
		if (var_x<max(rightx(w), leftx(w)) && var_x>min(rightx(w), leftx(w)))
			//v_corr+=dimdelta(w,0)*w[i]*w(pnt2x(w, i)+step)
			v_corr+=w[i]*w(pnt2x(w, i)+step)
		endif
	endfor
	return v_corr
end

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

	return w[0]*exp(-(x/(w[1]*sqrt(2)/2.354))^2)
End

function FitToGetK2(w)
	wave w // the correlation wave with 10001 points
	
	DFREF dfr = root:Packages:CorrelationAnalysis
	
	make /n=500 /o dfr:K2fit /WAVE=K2fit, dfr:K2fit_sd /WAVE=K2fit_sd
	K2fit=nan
	K2fit_sd=nan
	setscale /P x, 0, dimdelta(w, 0), K2fit, K2fit_sd
	variable startpoint=3
	
	variable i ,lowpnt=5000-startpoint, highpnt=5000+startpoint
	make /free /n=2 w_coef
	variable HeightGuess=w[5000]
	variable V_FitError=0
	for (i=startpoint;i<500;i+=1)
		w_coef={HeightGuess,50}
		V_FitError=0
		FuncFit /Q CorrelationGaussFit W_coef w[lowpnt,highpnt] /D 
		if(V_FitError==0 && abs(W_coef[1])<1e3)
			K2fit[i]=abs(W_coef[1])
			wave w_sigma=w_sigma
			K2fit_sd[i]=W_sigma[1]
		endif		
		lowpnt--; highpnt++	
	endfor	
end

function makeK2fitPlot()
	DFREF dfr = root:Packages:CorrelationAnalysis
	
	wave intercept=dfr:intercept, interceptSD=dfr:interceptSD, correlation=dfr:correlation
	wave K2fit=dfr:K2fit, K2fit_sd=dfr:K2fit_sd
	
	killwindow /Z K2FitGraph
	Display /K=1/W=(35,44,576,485)/N=K2FitGraph K2fit as "K2 fit"
	ModifyGraph rgb(K2fit)=(0,0,65535), standoff=0
	ErrorBars K2fit SHADE= {0,0,(0,0,0,0),(0,0,0,0)},wave=(K2fit_sd,K2fit_sd)

	doupdate
	getaxis /Q/W=K2FitGraph left
	variable midY=(V_max-V_min)/2
	Cursor/F/S=2/N=1/H=2/C=(65535,0,0) C  K2fit 5, midY
	Cursor/F/S=2/N=1/H=2/C=(65535,0,0) D  K2fit 10, midY
	AppendToGraph intercept
	ModifyGraph mode(intercept)=3,marker(intercept)=18,useMrkStrokeRGB(intercept)=1
	ErrorBars intercept Y,wave=(interceptSD,interceptSD)
	
	updateK2fitArea(5, 10)
	
	ModifyGraph rgb(intercept)=(3,52428,1), rgb(fit_K2fit)=(3,52428,1)
	
	SetWindow K2FitGraph, hook(K2Hook) = K2WindowHook
	
	Label left "K2 (cm\\S-1\\M)"
	Label bottom "Offset (cm\\S-1\\M)"
	
	Display/W=(0.5,0.5,1,0.9)/PG=($"",$"",PR,$"")/HOST=#  correlation
	ModifyGraph rgb(correlation)=(0,0,0)
	SetActiveSubwindow ##	
end

Function K2WindowHook(s)
	STRUCT WMWinHookStruct &s
	
	Variable hookResult = 0	
	switch(s.eventCode)
		case 7:
			if (stringmatch(s.cursorname, "C") || stringmatch(s.cursorname, "D"))
				 variable x1=hcsr(C, "K2FitGraph")
				 variable x2=hcsr(D, "K2FitGraph")
				 updateK2fitArea(x1, x2)	 
			endif
			hookResult = 1
			break					
	endswitch

	return hookResult	// If non-zero, we handled event and Igor will ignore it.
End

function updateK2fitArea(x1, x2)
	variable x1, x2
	
	DFREF dfr = root:Packages:CorrelationAnalysis
	wave K2fit=dfr:K2fit,K2fit_sd=dfr:K2fit_sd 
	wave intercept=dfr:intercept
	wave interceptSD=dfr:interceptSD
	wave /T textresults=dfr:textresults
	
	// redraw shded area
	SetDrawLayer /W=K2FitGraph ProgBack
	DrawAction /W=K2FitGraph getgroup=gK2, delete, begininsert
	SetDrawEnv /W=K2FitGraph gstart,gname=gK2
	SetDrawEnv /W=K2FitGraph xcoord= bottom,ycoord= prel,linethick= 0
	SetDrawEnv /W=K2FitGraph fillfgc= (65535,0,65535,5000)
	DrawRect /W=K2FitGraph x1,0,x2,1
	SetDrawEnv /W=K2FitGraph gstop,gname=gK2
	DrawAction /W=K2FitGraph endinsert
	
	// don't let vertical axis autoscale
	// poly fit can go to large y values
	getaxis /W=K2FitGraph/Q left 
	setaxis /W=K2FitGraph left, v_min, V_max 
	
	// update poly fit
	variable V_FitError=0
	CurveFit /Q/X=1 poly 3, K2fit(x1,x2) /W=K2fit_sd /I=1 /D 
	wave w_coef=w_coef, W_sigma=W_sigma
	intercept=W_coef[0]
	interceptSD=W_sigma[0]
	
	// update textbox
	string cmd
	sprintf cmd, "Range: %g-%g\rIntercept = %g ± %g", hcsr(C), hcsr(D), W_coef[0], W_sigma[0]
	TextBox /W=K2FitGraph /C/N=text0/F=0/A=MC cmd
	
end