#pragma rtGlobals=1		// Use modern global access method.
#pragma version=3.00
#pragma IgorVersion=6.10
//#pragma ModuleName=BaselineSpline
//#pragma hide=1

#include <Readback ModifyStr>

// written by tony withers, tony.withers@uwo.ca

// 3.00 4/11/16
// added option for Akima style spline

// 2.60 3/30/16
// added option to supress node addition
// this is a toggle that swiches on the /NI option in GraphWaveEdit
// updated a lot of other code
// changed some of the defaults (now saves nodes by default), so behaviour is not the same as previous version

// 2.50
// log mode to save node positions to file
// can load node positions from a file

// 2.40
// Added textbox reminder of how to toggle graph mode
// renamed procedures for consistency
// 2.30 5/16/11
// bug fix - code was choking on non-existant wave in first call to ACW_SplineClear() 
// added panel to change fitting options
// 2.10
// call ACW_SplineClear() from ACW_SplineSubtract()

// 2.0  8/6/09
// incorporated suggestions from JJ Weimer that improve functionality and clarity of code
// and compatibility with 6.10
// ... and then wrecked his nice programming with my clumsy code.

// 1.43b
// allow refit function to switch interpolate flag when number of nodes is small to allow linear fit

// 1.42 6/8/08
// set y axis to manual range to avoid crazy rescaling if spline shoots out of range.
// show subtracted spectrum whilst fitting

// 1.41 10/9/07
// for XY data we now interpolate directly into XY baseline
// improved guesses for initial node positions for unevenly spaced XY data

// 1.40   10/9/07
// fixed bug that occurred with baseline subtraction of XY data

// 1.30  7/24/07
// handles XY data
// baseline is waveform, but data wave need not be (though it should be reasonably closely spaced).

// 1.20 7/3/07
// uses new graphwaveedit flag. now requires version 6.0.1.


// note: editing these values won't change the action of the procedure
Static StrConstant thePackage="BaselineSpline"
Static StrConstant thePackageFolder="root:Packages:BaselineSpline"
Static StrConstant theProcedureFile = "BaselineSpline.ipf"
Static Constant thePackageVersion = 2.6
Static Constant hasHelp = 0

// uncomment the next line to use JJ Weimer style menu definitions:
//#define Weimerized

#ifdef Weimerized
	Menu "Misc"
		Submenu "Packages"
			"Spline Baseline",/Q, ACW_SplineInitFit()
		end
	end
	
	Menu "Macros"
		ACW_MenuSplineToggleFitting(),/Q, ACW_SplineToggleFitting()
		ACW_MenuSplineToggleInsert(), /Q, ACW_SplineToggleInsert()
		ACW_MenuSplineLoadNodes(),/Q, ACW_SplineLoadNodes()
		ACW_MenuSplineSubtract(),/Q, ACW_SplineSubtract()
		ACW_MenuSplineClear(),/Q, ACW_SplineClear()
		ACW_MenuSplineSettings(),/Q, ACW_SplineSettings() 
	end
#else
	Menu "Analysis"
		Submenu "Packages"
			"Spline Baseline",/Q, ACW_SplineInitFit()
		end
	end
	
	Menu "Macros"	
		Submenu "Spline Baseline"
			"Initalize...",/Q, ACW_SplineInitFit()
			"Settings...",/Q, ACW_SplineSettings() 		
			// following Macro menus appear when editing is possible	
			ACW_MenuSplineToggleFitting(), /Q, ACW_SplineToggleFitting()
			ACW_MenuSplineToggleInsert(), /Q, ACW_SplineToggleInsert()
			ACW_MenuSplineLoadNodes(),/Q, ACW_SplineLoadNodes()
			ACW_MenuSplineSubtract(), /Q , ACW_SplineSubtract()	
			ACW_MenuSplineClear(),/Q, ACW_SplineClear()
		end
	end
#endif

// Macro menus appear after fit is initialized
// avoid allowing dynamic menus to create the package folder before we even initialize.
function /s ACW_MenuSplineSettings()
	if (ACW_SplineCheckDFR()==0)
		return ""
	endif
	return "Spline Settings..."
end


function /s ACW_MenuSplineLoadNodes()
	if (ACW_SplineCheckDFR()==0)
		return ""
	endif
	DFREF dfr = ACW_SplineGetDFREF()
	NVAR V_edit= dfr:V_edit
	if (V_edit>-1)
		return "Load saved nodes..."
	endif
	return ""
end

Function/S ACW_MenuSplineToggleFitting()	
	if (ACW_SplineCheckDFR()==0)
		return ""
	endif
	DFREF dfr = ACW_SplineGetDFREF()
	NVAR V_edit= dfr:V_edit

	switch(V_edit)
		case -1:
			return ""
			break
		case 0:
			return "Adjust Nodes/1"
			break
		case 1:
			return "Adjust Nodes!"+num2char(18)+"/1"
			break
	endswitch
end

Function/S ACW_MenuSplineToggleInsert()	
	if (ACW_SplineCheckDFR()==0)
		return ""
	endif
	DFREF dfr = ACW_SplineGetDFREF()
	NVAR /SDFR=dfr V_NoInsert, V_edit
	if (V_edit==-1)
		return ""
	endif
	if(V_NoInsert)
		return "Allow New Nodes/2"
	endif
	return "Allow New Nodes!"+num2char(18)+"/2"
end

// Macro menu appears when spline curve is attached to graph
Function/S ACW_MenuSplineSubtract()
	if (ACW_SplineCheckDFR()==0)
		return ""
	endif
	DFREF dfr = ACW_SplineGetDFREF()
	NVAR V_edit= dfr:V_edit
	
	if (V_edit >= 0) // V_edit is set to -1 when we're done with fitting
		return "Subtract Baseline"
	else
		return ""
	endif
end

// clear spline line menu
Function/S ACW_MenuSplineClear()	
	if (ACW_SplineCheckDFR()==0)
		return ""
	endif
	DFREF dfr = ACW_SplineGetDFREF()
	NVAR V_edit= dfr:V_edit

	if (V_edit >= 0)
		return "Clear Spline"
	else
		return ""
	endif
end


// initialize spline baseline package
function ACW_SplineInitFit()
	
	if (strlen(WinList("*",";","WIN:1"))==0)  // no graphs 
		doalert 0, "Spline baseline requires a trace plotted in a graph window."
		return 0
	endif
	
	// make sure the top graph is visible
	dowindow /F $WinName(0,1)
	
	// clear any package detritus from the last used spline graph
	ACW_SplineClear() 
	
	// initialize (or reinitialize) the package data folder
	ACW_SplineResetPackageFolder()
	DFREF dfr = ACW_SplineGetDFREF()
		
	// define the globals
	svar /SDFR=dfr S_Data,S_Xwave,S_GraphName
	nvar /SDFR=dfr V_edit
	wave SplineFitGlobals=dfr:SplineFitGlobals
	wave W_nodesX=dfr:W_nodesX
	wave W_nodesY=dfr:W_nodesY
	wave W_spline_dependency=dfr:W_spline_dependency
	wave W_base=dfr:W_base
	
	// define some local variables	
	variable i, nthRange, nn=SplineFitGlobals[%nodes]-3
	nn=max(1, nn) // just making sure
	string s_traceName
	variable p_low, p_high, v1, v2
	variable pmin, pmax, delP
	string s_info
	variable trace_offset=0
	

	prompt s_traceName, "Data Wave", popup, TraceNameList("",";",1+4) // what does bit 2 do?
	prompt nn, "Number of Nodes", popup, "4;5;6;7;8;9;10;"
	DoPrompt "", s_traceName, nn
	if (V_Flag)
		ACW_SplineClear()
		return 0
	endif
	
	SplineFitGlobals[%nodes]=nn+3
	if (stringmatch(s_traceName, "_none_"))
		doalert 0, "No spectra in top graph window. Inititalisation failed."
		ACW_SplineClear()
		return 0
	endif
	
	S_Data=s_traceName
	
	S_GraphName=WinName(0, 1)
	string newName="SF_"+S_GraphName
	if (checkname(newName, 6)!=0) 
		newName=uniquename(newName, 6, 0)
	endif
	RenameWindow  $S_GraphName $newName
	S_GraphName=WinName(0, 1)
	
	wave data=TraceNametoWaveRef("",s_tracename)
	duplicate /O data dfr:W_base,  dfr:W_sub
	wave W_base=dfr:W_base
	wave w_sub= dfr:W_sub
	
	variable resetNodesX=1
	if(dimsize(W_nodesX, 0)==SplineFitGlobals[%nodes] && (SplineFitGlobals[%options]&2^1) )
		resetNodesX=0
	endif
	redimension/N=(SplineFitGlobals[%nodes]) W_nodesX, W_nodesY
	
	SplineFitGlobals[%isXY]=0
	// find out whether data is XY
	S_Xwave=StringByKey("XWAVE", TraceInfo("",s_tracename,0))
	if (strlen(S_Xwave))
		wave Xwave=XWaveRefFromTrace("", s_tracename )
		S_Xwave[0]=StringByKey("XWAVEDF", TraceInfo("",s_tracename,0))
		SplineFitGlobals[%isXY]=1
		
		// check for monotonic X wave
		Differentiate Xwave /D=dfr:W_diff
		wave W_diff=dfr:W_diff
		Wavestats /Q W_diff
		killwaves W_diff
		if (V_max>0 && V_min<0)
			doalert 0, "X wave is not monatonic!"			
			ACW_SplineClear()
			return 0
		endif
	
	endif
		
	// manually scale y axis
	getaxis /Q left
	setAxis left, V_min, V_max
	
	// identify trace to be corrected
	if (SplineFitGlobals[%options]&2^0)
		ACW_SplineHighlightData(1)
	endif
	
	// determine where to place nodes on the bottom axis	
	if(SplineFitGlobals[%fitWithinGraph]) // keep nodes within window 
		GetAxis/Q bottom
		
		if (SplineFitGlobals[%isXY])
			findlevel /Q/P Xwave, V_max
			if (V_flag)
				v1=numpnts(Xwave)
			else
				v1=V_LevelX
			endif
			findlevel /Q/P Xwave, V_min
			if (V_flag)
				v2=0
			else
				v2=V_LevelX
			endif
			
			pmax=max(v1,v2)
			pmin=min(v1,v2)
		else
			pmax=floor(max( x2pnt(data, V_max), x2pnt(data,V_min) ) ) -1
			pmin=ceil (min( x2pnt(data,V_max), x2pnt(data,V_min) ) ) +1
		endif
	else // spead nodes across entire range of wave and autoscale bottom axis
		pmin=0
		pmax=numpnts(data)-1

		// autoscale bottom axis before adding spline nodes
		// baseline is subtracted from entire wave, so show the whole range.
		getaxis /Q  bottom
		if (v_min>v_max)
			setaxis /r/a bottom
		else
			setaxis /A bottom
		endif		
	endif
	
	delP=(pmax-pmin)/(SplineFitGlobals[%nodes]-1)
	
	// initialize node positions. set nodes to follow roughly the data	
	if (SplineFitGlobals[%isXY])
		nthRange=(Xwave[pmax]-Xwave[pmin])/(SplineFitGlobals[%nodes]-1)
		if (resetNodesX)
			W_nodesX=Xwave[pmin]+p*nthRange
		endif
	else
		nthRange=(pnt2x(data, pmax)-pnt2x(data, pmin))/(SplineFitGlobals[%nodes]-1)
		if (resetNodesX)
			W_nodesX=pnt2x(data, pmin)+p*nthRange
		endif
	endif
		
	for (i=0;i<(SplineFitGlobals[%nodes]);i+=1)
		if (SplineFitGlobals[%isXY])
			// necessary for unevenly spaced data
			findlevel /Q/P Xwave, W_nodesX[i]-nthRange/24
			if (V_flag)
				p_low=pmin
			else
				p_low= V_LevelX
			endif
			
			findlevel /Q/P Xwave, W_nodesX[i]+nthRange/24
			if (V_flag)
				p_high=pmax
			else
				p_high= V_LevelX
			endif		
		else					
			p_low=  max(pmin, pmin+i*delP - delP/24 ) 
			p_high= min(pmax, pmin+i*delP + delP/24 )
		endif
		wavestats /Q/M=1/R=[p_low , p_high] /Z  data
		W_nodesY[i]=V_avg  // this handles a few NANs	
	endfor

	
	// if data are offset, then apply offsets to baseline and nodes
	s_info=traceinfo("",s_tracename,0)
	trace_offset=GetNumFromModifyStr(s_info,"offset","{",1)
	
	
	AppendToGraph W_spline_dependency // need to have this on graph to force update while in drawing mode
	ModifyGraph hideTrace(W_spline_dependency)=1
	
	appendtograph W_nodesY vs W_nodesX
	ModifyGraph mode(W_nodesY)=3,marker(W_nodesY)=11,rgb(W_nodesY)=(0,39168,0)
	ModifyGraph offset(W_nodesY)={0,trace_offset}
	
	ACW_SplineUpdate(W_nodesX)
	
	if (SplineFitGlobals[%isXY])
		appendtograph W_base vs Xwave
		if (SplineFitGlobals[%showSub])
			appendtograph W_sub vs Xwave
		endif
	else
		appendtograph W_base
		if (SplineFitGlobals[%showSub])
			appendtograph W_sub
		endif
	endif
	
	ModifyGraph offset(W_base)={0,trace_offset}
	ModifyGraph rgb(W_base)=(0,0,52224)
	modifygraph live(W_base)=1
	
	if (SplineFitGlobals[%showSub])
		ModifyGraph rgb(W_sub)=(0,0,0)
		ModifyGraph zero(left)=1
		W_sub=data-W_base
	endif
	
	// set a dependency to trigger interpolation AND update graph when we adjust a node position
	W_spline_dependency:=ACW_SplineUpdate(root:Packages:SplineFit:W_nodesY)	
	// won't work if we change package folder
	
	ReorderTraces $"W_nodesY",{$"W_base"}
	V_edit = 0
	ACW_SplineToggleFitting() // turn editing on	
	
	if(SplineFitGlobals[%options]&2^2) //verbose mode	
		print "Spline fit initialized with data wave "+GetWavesDataFolder(data, 2 )
	endif
	
	ACW_SplineEnableControls() // in case settings panel is open
	
	return 0	
end

// initialize (or reinitialize) the package data folder
Function /DF ACW_SplineResetPackageFolder()
	
	NewDataFolder /O root:Packages
	NewDataFolder /O root:Packages:SplineFit	
	// Create a data folder reference variable
	DFREF dfr = root:Packages:SplineFit
	
	string /G dfr:S_GraphName=""
	string /G dfr:S_Data=""
	string /G dfr:S_Xwave=""
	variable /G dfr:V_edit=0  // global var to store edit status. 
	variable /G dfr:V_NoInsert=0 // allow node insertion
	wave /Z W_nodesX=dfr:W_nodesX
	if (WaveExists(W_nodesX)==0)	
		make /O/n=(1) dfr:W_nodesX, dfr:W_nodesY
	endif
	make /O/n=1  dfr:W_spline_dependency
	make/O/n=1  dfr:W_base
	
	// create a wave to store some global variables	
	wave /Z SplineFitGlobals=dfr:SplineFitGlobals
	if (waveexists(w_globals)==0) // retain these settings when we reinitialize
		make /o/n=10 dfr:SplineFitGlobals=nan 
		wave SplineFitGlobals=dfr:SplineFitGlobals
		setdimlabel  0, 0, null, SplineFitGlobals
		setdimlabel  0, 1, flagE, SplineFitGlobals // flag for interpolate2
		setdimlabel  0, 2, flagF, SplineFitGlobals // flag for interpolate2
		setdimlabel  0, 3, flagJ, SplineFitGlobals // flag for interpolate2
		setdimlabel  0, 4, flagT, SplineFitGlobals // flag for interpolate2
		setdimlabel  0, 5, isXY, SplineFitGlobals
		setdimlabel  0, 6, showSub, SplineFitGlobals
		setdimlabel  0, 7, nodes, SplineFitGlobals
		setdimlabel  0, 8, fitWithinGraph, SplineFitGlobals
		setdimlabel  0, 9, options, SplineFitGlobals 
		// options
		// bit 0: highlight data 
		// bit 1: remember nodes
		// bit 2: verbose mode
		// bit 3: save nodes
		// bit 4: use Akima spline
		SplineFitGlobals[%flagE]=2
		SplineFitGlobals[%flagF]=1
		SplineFitGlobals[%flagJ]=1
		SplineFitGlobals[%flagT]=2
		SplineFitGlobals[%showSub]=1
		SplineFitGlobals[%nodes]=5
		SplineFitGlobals[%fitWithinGraph]=1 // place nodes within range of graph
		SplineFitGlobals[%options]=2^0+2^2+2^3
	endif
	
	return dfr
end

Function/DF ACW_SplineGetDFREF()
	DFREF dfr = root:Packages:SplineFit
	if (DataFolderRefStatus(dfr) != 1)	
		DFREF dfr = ACW_SplineResetPackageFolder()
	endif
	return dfr
End

Function ACW_SplineCheckDFR()
	DFREF dfr = root:Packages:SplineFit
	if (DataFolderRefStatus(dfr) != 1)	
		return 0
	endif
	return 1
End

// dependency function to define the spline curve
function ACW_SplineUpdate(w)
	wave w
	
	DFREF dfr = ACW_SplineGetDFREF()
	
	wave g=dfr:SplineFitGlobals
	wave W_nodesX=dfr:W_nodesX
	wave W_nodesY=dfr:W_nodesY
	wave W_base=dfr:W_base
	
	SVAR S_Xwave=dfr:S_Xwave
	
	// check for small number of points
	variable SplineType=g[%flagT]
	if (numpnts(W_nodesX)<4)
		SplineType=1
	endif
	
	if ((g[%options]&2^4) && numpnts(W_nodesX)<5)
		doalert 0, "Akima spline requires at least 5 nodes\r - reverting to standard cubic spline"
		g[%options]=g[%options]%^2^4
		ACW_SplineEnableControls()
	endif
	
	if (g[%isXY])	
		wave w_x=$S_Xwave
		if (g[%options]&2^4)
			if (ACW_SplineCalcIota(W_nodesX, W_nodesY))
				wave w_iota=dfr:w_iota
				
				MultiThread W_base=ACW_SplineAkima(w_x[p], W_nodesX, W_nodesY, w_iota)
			endif
		else
			Interpolate2 /E=(g[%flagE]) /F=(g[%flagF]) /J=(g[%flagJ]) /T=(SplineType) /I=3  /X=w_x /Y=W_base W_nodesX, W_nodesY
		endif
	else
		if (g[%options]&2^4)
			if (ACW_SplineCalcIota(W_nodesX, W_nodesY))
				wave w_iota=dfr:w_iota	
				MultiThread W_base=ACW_SplineAkima(x, W_nodesX, W_nodesY, w_iota)
			endif
		else
			Interpolate2 /E=(g[%flagE]) /F=(g[%flagF]) /J=(g[%flagJ]) /T=(SplineType) /I=3  /Y=W_base  W_nodesX, W_nodesY
		endif
	endif
	
	if (g[%showSub])
		SVAR s_data=dfr:S_Data
		wave data=TraceNametoWaveRef("",S_Data) 
		wave w_sub= dfr:w_sub
		w_sub=data-w_base
	endif
	
	return nan	
end

// function to subtract the spline curve
function ACW_SplineSubtract()
	
	DFREF dfr = ACW_SplineGetDFREF()
	
	SVAR s_data=dfr:S_Data
	NVAR V_edit=dfr:V_edit
	wave SplineFitGlobals=dfr:SplineFitGlobals
	wave W_base=dfr:W_base
	wave data=TraceNametoWaveRef("",S_Data)
	wave W_nodesX=dfr:W_nodesX
	wave W_nodesY=dfr:W_nodesY
	string strSubtracted, strBaseline, strNodes
	strSubtracted=CleanupName( nameofwave(data)+"_SpBLCorr",0)
	strBaseline=CleanupName( nameofwave(data)+"_SpBL",0)
	strNodes=CleanupName( nameofwave(data)+"_SpNodes",0)

	if (exists(strSubtracted))
		doalert 1, strSubtracted+" exists. Overwrite?"
		if(V_flag==2)
			print "Spline baseline subtraction aborted"
			return 0
		endif
	endif
	duplicate /o data $strSubtracted
	wave nob=$strSubtracted
	
	if (exists(strBaseline))
		doalert 1, strBaseline+" exists. Overwrite?"
		if(V_flag==2)
			print "Spline baseline subtraction aborted"
			return 0
		endif
	endif
	duplicate /o W_base $strBaseline
	
	if (SplineFitGlobals[%options]&2^3)
		if (exists(strNodes))
			doalert 1, strNodes+" exists. Overwrite?"
			if(V_flag==2)
				print "Spline baseline subtraction aborted"
				return 0
			endif
		endif
		duplicate /o W_nodesX $strNodes
		wave nd=$strNodes
		redimension /N=(-1, 2) nd
		nd[][0]=W_nodesX[p]
		nd[][1]=W_nodesY[p]
		if(SplineFitGlobals[%options]&2^2) //verbose mode
			print "Spline fit saved node positions in wave "+GetWavesDataFolder(nd, 2 )	
		endif
	endif
	
	
	wave bl=$strBaseline
	
	nob=data-bl
		
	if(SplineFitGlobals[%options]&2^2) //verbose mode	
		print "Spline fit created baseline wave "+GetWavesDataFolder(bl, 2 )
		print "Spline fit created baseline subtracted wave "+GetWavesDataFolder(nob, 2 )
		if(SplineFitGlobals[%options]&2^4)
			print "Spline fit used Akima type spline"
		endif
	endif
	
	removefromgraph /Z $strSubtracted, $strBaseline, W_base,  W_nodesY, w_sub
		
	ACW_SplineClear()	
	
	// find out whether data is XY	
	string s_Xwave=StringByKey("XWAVE", TraceInfo("",NameOfWave(data),0))
	
	if (strlen(s_Xwave))  // XY data
		wave Xwave=XWaveRefFromTrace("", NameOfWave(data) )	
		appendtograph nob, bl vs Xwave
	else
		appendtograph nob, bl
	endif
	
	// would be nice to put some of this in the history...
	
	ModifyGraph rgb($nameofwave(bl))=(0,65280,0)	
	ModifyGraph rgb($nameofwave(nob))=(0,0,0)	
	V_edit=-1
	BuildMenu "macros"	
	
	return 1
end

// function to toggle between normal graph and draw mode
function ACW_SplineToggleFitting()
	DFREF dfr = ACW_SplineGetDFREF()

	NVAR /SDFR=dfr V_edit, V_NoInsert
	SVAR /SDFR=dfr S_GraphName
	
	variable isMac=abs(cmpstr(IgorInfo(2)[0], "W"))
	string msg
	
	if (V_edit)
		GraphNormal /W=$S_GraphName
		ModifyGraph /W=$S_GraphName mode(W_nodesY)=2
		V_edit=0		
	else
		ModifyGraph /W=$S_GraphName mode(W_nodesY)=3
		if (V_NoInsert)
			GraphWaveEdit /W=$S_GraphName /NI/T=1 /M $"W_nodesY"
		else
			GraphWaveEdit /W=$S_GraphName /T=1 /M $"W_nodesY"
		endif
		V_edit=1
	endif
	
	sprintf msg, "\Z14\\K(32768,65280,0)%s-1 to %s node editing", selectstring(isMac, "CTRL", "CMD"), selectstring(V_edit, "resume", "pause")
	TextBox/C/N=SplineFitTextBox/F=0/B=1/A=LT/X=5.00/Y=5.00 msg
	
	buildmenu "macros"
	return 1
end


function ACW_SplineToggleInsert()
	DFREF dfr = ACW_SplineGetDFREF()

	NVAR /SDFR=dfr V_edit, V_NoInsert
	SVAR /SDFR=dfr S_GraphName
	V_NoInsert=1-V_NoInsert
	if (V_edit)
		if (V_NoInsert)
			GraphWaveEdit /W=$S_GraphName /NI/T=1 /M $"W_nodesY"
		else
			GraphWaveEdit /W=$S_GraphName /T=1 /M $"W_nodesY"
		endif
	endif
	buildmenu "macros"
	return 1
end


// function to clear the spline curve
function ACW_SplineClear()
	DFREF dfr = ACW_SplineGetDFREF()

	NVAR V_edit=dfr:V_edit
	SVAR /SDFR=dfr S_GraphName, S_Data	
	wave W_spline_dependency=dfr:W_spline_dependency
	wave SplineFitGlobals=dfr:SplineFitGlobals
	
	V_edit = -1
	W_spline_dependency=NaN
	
	if (strlen(S_GraphName)&&WinType(S_GraphName))
		GraphNormal /W=$S_GraphName
		RemoveFromGraph /W=$S_GraphName/Z W_spline_dependency,W_base,W_nodesY, w_sub
		TextBox /W=$S_GraphName/K/N=SplineFitTextBox
		
		// undo highlighting of data		
		if (SplineFitGlobals[%options]&2^0)
			ACW_SplineHighlightData(0)
		endif
		
		// remove prepended SF_ from graph name
		if (WinType(S_GraphName))
			string newName=S_GraphName[3,inf]
			if (checkname(newName, 6)!=0) // in case another window stole the name
				newName=uniquename(newName, 6, 0)
			endif
			RenameWindow  $S_GraphName $newName
		endif
	endif
	
	S_Data=""	
	S_GraphName=""	
	BuildMenu "macros"	
	return 1
end


Function ACW_SplineCheckboxes(ctrlName,checked) : CheckBoxControl
	String ctrlName
	Variable checked
	
	DFREF dfr = ACW_SplineGetDFREF()
	wave SplineFitGlobals=dfr:SplineFitGlobals
	SVAR /SDFR=dfr S_Data, S_GraphName
	
	strswitch(ctrlName)
		case "checkGraphWidth":
			SplineFitGlobals[%fitWithinGraph]=checked
			return 1
		case "checkShowSub":
			SplineFitGlobals[%showSub]=checked
			if (SVAR_exists(S_GraphName))
				dowindow $S_GraphName
				if (V_flag==1)
					if (checked)
						checkdisplayed /W=$S_GraphName w_sub
						if (V_flag<1) 
							wave w_sub= root:Packages:SplineFit:W_sub
							appendtograph /W=$S_GraphName W_sub
							ModifyGraph /W=$S_GraphName rgb(W_sub)=(0,0,0)
							ModifyGraph /W=$S_GraphName zero(left)=1
						endif
					else
						dowindow $S_GraphName
						if (V_flag==1)
							removefromgraph /W=$S_GraphName /Z w_sub
						endif
					endif
				endif
			endif
			return 1
		case "checkHighlightData":
			SplineFitGlobals[%options]=SplineFitGlobals[%options]%^2^0 // toggle 1st bit		
			ACW_SplineHighlightData(checked) // update graph if neccessary
			return 1		
		case "checkRememberNodes":
			SplineFitGlobals[%options]=SplineFitGlobals[%options]%^2^1
			return 1
		case "checkVerbose":
			SplineFitGlobals[%options]=SplineFitGlobals[%options]%^2^2
			return 1		
		case "checkSaveNodes":
			SplineFitGlobals[%options]=SplineFitGlobals[%options]%^2^3
			return 1	
		case "checkAkima":
			SplineFitGlobals[%options]=SplineFitGlobals[%options]%^2^4
			
			ACW_SplineEnableControls()
			NVAR V_edit=dfr:V_edit
			if(V_edit==1)
				make /free/n=1 w
				ACW_SplineUpdate(w)
			endif
			
			
			return 1	
	endswitch
	return 0
End

Function ACW_SplineButtons(ctrlName)
	String ctrlName
	DisplayHelpTopic "interpolate2"
End


Function ACW_SplinePopMenu(ctrlName,popNum,popStr)
	String ctrlName
	Variable popNum
	String popStr
	
	DFREF dfr = ACW_SplineGetDFREF()
	wave SplineFitGlobals=dfr:SplineFitGlobals
	
	strswitch(ctrlName)
		case "popE":
			SplineFitGlobals[%flagE]=popNum
			return 1
		case "popJ":
			SplineFitGlobals[%flagJ]=popNum-1
			return 1
		case "popT":
			SplineFitGlobals[%flagT]=popNum
			return 1	
	endswitch
	return 0
End

//Create panel to edit settings
Function ACW_SplineSettings() 
	DFREF dfr = ACW_SplineGetDFREF()
	
	wave SplineFitGlobals=dfr:SplineFitGlobals
	DoWindow/K SplineFitPanel
	NewPanel /N=SplineFitPanel /K=1/W=(549,75,854,350) as "Spline Fit Parameters"
	ModifyPanel fixedSize=1
	SetDrawLayer ProgBack
	
	GroupBox groupFlags,pos={20,19},size={247,51},title="Interpolate2 flags"
	
	PopupMenu popE,pos={31,42},size={46,21}, title="/E", proc=ACW_SplinePopMenu
	PopupMenu popE,value="1;2", mode=(SplineFitGlobals[%flagE]), popvalue=num2str(SplineFitGlobals[%flagE])
	PopupMenu popE,disable= 2*(SplineFitGlobals[%options]&2^4)
	SetVariable setvarFlagF,pos={90,44},size={45,16},bodyWidth=30,title="/F"
	SetVariable setvarFlagF,value= SplineFitGlobals[%flagF], limits={0,inf,0.1}	
	SetVariable setvarFlagF,disable= 2*(SplineFitGlobals[%options]&2^4)
	PopupMenu popJ,pos={150,42},size={44,21},title="/J", proc=ACW_SplinePopMenu
	PopupMenu popJ,value="0;1;2", mode=(SplineFitGlobals[%flagJ]+1), popvalue=num2str(SplineFitGlobals[%flagJ])
	PopupMenu popJ,disable= 2*(SplineFitGlobals[%options]&2^4)
	PopupMenu popT,pos={210,42},size={46,21},title="/T", proc=ACW_SplinePopMenu
	PopupMenu popT,value="1;2;3", mode=(SplineFitGlobals[%flagT]), popvalue=num2str(SplineFitGlobals[%flagT])
	PopupMenu popT,disable= 2*(SplineFitGlobals[%options]&2^4)
	
	CheckBox checkAkima,pos={33.00,79.00},size={73.00,16.00},title="Akima Spline"
	CheckBox checkAkima,value= (SplineFitGlobals[%options]&2^4), proc=ACW_SplineCheckboxes
	
	Button buttonHelp,pos={150,79},size={125,20},proc=ACW_SplineButtons,title="Interpolate2 help"
	
	variable x0=38, x1=153, y0=140
	GroupBox groupSettings,pos={20.00,118.00},size={247.00,141.00}
	SetVariable setvarNodes,pos={x0,y0},size={63,16},bodyWidth=30,title="nodes"
	SetVariable setvarNodes,value= SplineFitGlobals[%nodes], limits={5,inf,1}
	CheckBox checkVerbose,pos={x0,y0+30},size={88,14},proc=ACW_SplineCheckboxes,title="Write to history"
	CheckBox checkVerbose,value= SplineFitGlobals[%options]&2^2
	CheckBox checkShowSub,pos={x0,y0+60},size={100,14},proc=ACW_SplineCheckboxes,title="Show subtraction"
	CheckBox checkShowSub,value= SplineFitGlobals[%showSub]	
	CheckBox checkRememberNodes,pos={x0,y0+90},size={140,14},proc=ACW_SplineCheckboxes,title="Remember node positions"
	CheckBox checkRememberNodes,value= SplineFitGlobals[%options]&2^1	
	
	CheckBox checkSaveNodes,pos={x1,y0},size={140,14},proc=ACW_SplineCheckboxes,title="Save node positions"
	CheckBox checkSaveNodes,value= SplineFitGlobals[%options]&2^3	
	CheckBox checkHighlightData,pos={x1,y0+30},size={119,14},proc=ACW_SplineCheckboxes,title="Highlight data wave"
	CheckBox checkHighlightData,value= (SplineFitGlobals[%options]&2^0)	
	CheckBox checkGraphWidth,pos={x1,y0+60},size={89,14},proc=ACW_SplineCheckboxes,title="Fit within graph"
	CheckBox checkGraphWidth,value= 1
	
	return 1
End

function ACW_SplineEnableControls()
	DFREF dfr = ACW_SplineGetDFREF()	
	wave SplineFitGlobals=dfr:SplineFitGlobals
	
	PopupMenu /Z popE,win=SplineFitPanel, disable= 2*(SplineFitGlobals[%options]&2^4)/16
	PopupMenu /Z popJ,win=SplineFitPanel, disable= 2*(SplineFitGlobals[%options]&2^4)/16
	PopupMenu /Z popT,win=SplineFitPanel, disable= 2*(SplineFitGlobals[%options]&2^4)/16
	SetVariable /Z setvarFlagF,win=SplineFitPanel,disable= 2*(SplineFitGlobals[%options]&2^4)/16
	// make sure we don't get out of sync:
	CheckBox /Z checkAkima,win=SplineFitPanel,value= (SplineFitGlobals[%options]&2^4)
end


// function to highlight selected data wave in graph window
// increases line or marker size by 1 to highlight
Function ACW_SplineHighlightData(On)
	variable On // 1 to highlight, 0 to undo
	
	DFREF dfr = ACW_SplineGetDFREF()
	
	SVAR /SDFR=dfr S_Data, S_GraphName
	
	if (!(strlen(S_GraphName)))
		return 0
	endif	
	if (!(strlen(S_Data)))
		return 0
	endif
	
	dowindow $S_GraphName
	if (V_flag!=1)
		return 0
	endif
	
	string S_trace=TraceInfo(S_GraphName, S_Data, 0)
	variable mode=GetNumFromModifyStr(S_trace,"mode","",0)
	variable size
	if (mode==3) //data is plotted with markers
		size=GetNumFromModifyStr(S_trace,"msize","",0)
		if (size==0) // don't mess with auto sized markers
			return 0
		endif
		if (On)
			ModifyGraph/Z /W=$S_GraphName msize($S_Data)=size+1
		else
			ModifyGraph/Z /W=$S_GraphName msize($S_Data)=size-1
		endif
	else
		size=GetNumFromModifyStr(S_trace,"lSize","",0)
		if (On)
			ModifyGraph/Z /W=$S_GraphName lSize($S_Data)=size+1
		else
			ModifyGraph/Z /W=$S_GraphName lSize($S_Data)=size-1
		endif
	endif	
	
	return 1
End

Function ACW_SplineLoadNodes()
	DFREF dfr = ACW_SplineGetDFREF()
	
	wave SplineFitGlobals=dfr:SplineFitGlobals
		
	string str_NodeWave
	variable LoadY
	
	prompt str_NodeWave, "Select saved nodes", popup, wavelist("*",";","DIMS:2") 
	prompt LoadY, " ", popup, "Load X positions only;Load X and Y positions;"
	DoPrompt "", str_NodeWave, LoadY
	if (V_Flag)
		return 0
	endif
	
	wave w_nodes=$str_NodeWave
	loadY-=1
	
	wave W_nodesX=dfr:W_nodesX
	wave W_nodesY=dfr:W_nodesY
	
	redimension /N=(dimsize(w_nodes, 0)), W_nodesX
	redimension /N=(dimsize(w_nodes, 0)), W_nodesY
	
	W_nodesX=w_nodes[p][0]
	
	SplineFitGlobals[%nodes]=numpnts(W_nodesX)	
	
	if (LoadY)
		W_nodesY=w_nodes[p][1]
		return 1
	endif

	SVAR s_data=dfr:S_Data
	wave /Z data=TraceNametoWaveRef("",S_Data) 
	if (waveexists(data)==0)
		return 0
	endif
		
	// figure out Y values from data
	// average data over range 5 points either side of X position
	variable i, p_low, p_high
	
	for (i=0;i<(SplineFitGlobals[%nodes]);i+=1)
		if (SplineFitGlobals[%isXY])
			// necessary for unevenly spaced data
			SVAR S_Xwave=dfr:S_Xwave
			wave xwave=$S_Xwave
			
			findlevel /Q/P Xwave, W_nodesX[i]
			if (V_flag)
				return 0
			else
				p_low= V_LevelX-5
				p_high= V_LevelX+5
			endif					
		else					
			p_low= x2pnt(data, W_nodesX[i])-5
			p_high= x2pnt(data, W_nodesX[i])+5
		endif
		wavestats /Q/M=1/R=[p_low , p_high] /Z  data
		W_nodesY[i]=V_avg  // this handles a few NANs	
	endfor
	return 1
end



// **********     Akima spline code     *************

// this procedure contains code posted to IgorExchange by Michael Bongard
// a couple of edits by tony.withers@uwo.ca to allow extrapolation beyond end nodes
 
// Akima.ipf: Routines to implement Akima-spline fitting, based on
// H. Akima, Journ. ACM, Vol 17, No 4, 1970 p 589-602
// M. Bongard, 11/17/09

// First, call calcIota to generate interpolation information; 
// then you can interpolate using Akima's spline method with the akima() function.
 
Function ACW_SplineCalcIota(knotX, knotY)
	WAVE knotX // knot X locations
	WAVE knotY // knot Y locations
 	
 	// removed some sanity checks
 	DFREF dfr = ACW_SplineGetDFREF()
 	
	Variable numKnots = numpnts(knotX)
 
	if(numKnots < 5)
		Print "calcIota: ERROR -- Akima spline algorithm requires at least 5 knots. Aborting..."
		return -1
	endif
 
	// Make intermediate ai, bi, mi arrays
	Make/D/FREE/N=(numKnots + 4)  kX, kY
 
	Variable i, j
	for(i = 2, j = 0; j < numKnots; i += 1, j += 1)
		kX[i] = knotX[j]
		kY[i] = knotY[j]
	endfor
 
	// Handle end-point extrapolation
	Make/D/FREE/N=5 endX, endY
	// RHS: end points are last three in dataset
	Variable endStartPt
	// RHS
	endStartPt = numPnts(kX) - 5
	endX = kX[p + endStartPt]
	endY = kY[p + endStartPt]
 
	ACW_SplineCalcEndPoints(endX, endY)
 
	kX[numpnts(kX)-2] = endX[3]
	kX[numpnts(kX)-1] = endX[4]
	kY[numpnts(kX)-2] = endY[3]
	kY[numpnts(kX)-1] = endY[4]
 
	// LHS: end points are first three in dataset, but reversed in ordering
	// (i.e. point 3 in Akima's notation == index 0)
	endX = 0
	endY = 0
 
	for(i = 0, j = 2; i < 3; i += 1, j -= 1)
		endX[j] = knotX[i]
		endY[j] = knotY[i]
	endfor
 
	ACW_SplineCalcEndPoints(endX, endY)
 
	kX[1] = endX[3]
	kX[0] = endX[4]
	kY[1] = endY[3]
	kY[0] = endY[4]
 
	// kX, kY are now properly populated, along with all necessary extrapolated endpoints
	// computed as specified in Akima1970
 
	Make/D/FREE/N=(numKnots + 4 - 1)  mK
	mK = (kY[p + 1] - kY[p]) / (kX[p + 1] - kX[p])	
 
	Make/O/N=(numKnots) dfr:w_Iota /wave=w_iota
 
	Variable denom, m1,m2,m3,m4
	for(i = 2, j = 0; j < numKnots; i += 1, j += 1)
		m1 = mK[i - 2]
		m2 = mK[i - 1]
		m3 = mK[i]
		m4 = mK[i + 1]
 
		denom = abs(m4 - m3) + abs(m2 - m1)
		if(denom == 0)
			w_Iota[j] = 0.5 * (m2 + m3)
			continue
		endif
 
		w_Iota[j] = ( abs(m4 - m3)*m2 + abs(m2 - m1)*m3 ) / denom
	endfor
	return 1
End
 
// Given: 5-point knot wave knotX, knotY, with i=[0,2] representing the last three
// knot locations from data, compute end knots i=[3,4] appropriately.
ThreadSafe Function ACW_SplineCalcEndPoints(kX, kY)
	WAVE kX, kY // knot X,Y coordinate locations, respectively
 
	// Sanity checks
	if( (numPnts(kX) != numpnts(kY)) || numpnts(kX) != 5)
		Print "calcEndPoints: ERROR -- must have 5 points in knot wave! Aborting..."
		return -1
	endif
 
	// First, compute X locations of knots, according to relations in eq. 8 of Akima1970:
	kX[3] = kX[1] + kX[2] - kX[0]
	kX[4] = 2*kX[2] - kX[0]
 
	// Now all kX are known, so let's set up the line segment slope waves
	// ai, bi, mi of eq. (12)-(14).
	Make/N=4/FREE ai, bi, mi
 
	ai = kX[p+1] - kX[p]
	// ai is now determined completely
 
	bi = kY[p + 1] - kY[p]
	mi = bi/ai
	// bi, mi determined on i=[0,1]	
 
	// Determine remainder of quantities by applying solutions of eq (9)
	kY[3] = (2*mi[1] - mi[0])*(kX[3] - kX[2]) + kY[2]
	mi[2] = (kY[3] - kY[2]) / (kX[3] - kX[2])
 
	kY[4] = (2*mi[2] - mi[1])*(kX[4] - kX[3]) + kY[3]
	mi[3] = (kY[4] - kY[3]) / (kX[4] - kX[3])
End
 
ThreadSafe Function ACW_SplineAkima(x, knotX, knotY, knotIota)
	Variable x
	WAVE knotX, knotY, knotIota
 
	Variable i1, i2
 
	// Find where x
	Variable done = 0
	i2 = -1
	Variable numKnots = numpnts(knotX)
	do
		i2 += 1
 
		if(knotX[i2] == x)
			return knotY[i2]
		endif
 
		done = (knotX[i2] > x) ? 1 : 0
 
	while(!done && (i2 < numKnots-1)) // edited this to allow extrapolation beyond end knot
	i1 = i2 - 1
	
	// this edit to allow extrapolation beyond start knot
	if (i1<0)
		i1+=1
		i2+=1
	endif
 
	Variable x1, x2, y1, y2, iota1, iota2
	x1 = knotX[i1]
	y1 = knotY[i1]
	iota1 = knotIota[i1]
 
	x2 = knotX[i2]
	y2 = knotY[i2]
	iota2 = knotIota[i2]
 
	Variable p0, p1, p2, p3, tmp
	p0 = y1
	p1 = iota1
	p2 = ( 3 * (y2 - y1)/(x2 - x1) - 2*iota1 - iota2 ) / (x2 - x1)
	p3 = (iota1 + iota2 - 2 * (y2 - y1)/(x2 - x1)) / (x2 - x1)^2
 
	tmp = x - x1
 
	return p0 + p1 * tmp + p2 * tmp^2 + p3*tmp^3
End