// Peak Measurement Procedures
//	V2.0- 11/16/1999, rewritten to be a separate procedure file.
//	It still places way too many variables and strings in the root: data folder.

#pragma rtGlobals=1		// Use modern global access method.

Menu "Macros"
	"Initialize Most EverythingÉ"
	"Make PeaksÉ"
	"-"
	Submenu "Baseline"
		"Init Baseline Fit..."
		"Add Region To Fit/1"
		"Delete Region From Fit/2"
		"Fit Baseline At RegionsÉ"
		"Show Fit ResidualÉ"
		"Subtract Baseline..."
	End
	Submenu "Area"
		"Init Area Between CursorsÉ"
		"Area Between Cursors/3"
		"Hide Area Tags"
		"List Areas"
	End
	Submenu "Fit Peaks"
		"Init Identify PeaksÉ"
		"Identify Peak With Cursors/4"
		"Delete Pks Between Cursors/5"
		"Auto Identify PeaksÉ/6"
		"-"
		"Fit Peaks And BaselineÉ/7"
		"-"
		"ReportÉ"
		"Show Fitted PeaksÉ"
		"Hide Fitted Peaks"
		"Show Fit ResidualÉ"
	End
	"-"
	"Show Only Data And Base/0"
	"Color Waves"
End


Proc InitPeakMacros()
	Silent 1;PauseUpdate

	String/G g_w,g_wx,g_b,g_keep
	String/G p_w,p_wx,p_b,p_wd,p_wxd,p_bx
	p_w = "peak data wave, including baseline"
	p_wd = "data wave"
	p_wx = "X coordinates for peak data wave"
	p_wxd = "X coordinates for data wave"
	p_b = "baseline"
	p_bx="approximate peak x width at 50%"

	Variable/G g_bx

	Variable/G mp_kind,mp_pnts=8,mp_noise,mp_log,mp_minx,mp_maxx
	String/G p_pts="(2;(4;(8;16;32;64;128;256;512;1024;2048;4096;8192;16384"

	String/G fbar_wr,fbar_fit,fbar_wobase

	Variable/G mpr_pol

	Variable/G mir_what,mir_anno
	String/G mir_pwh="counts;integrate rectangular;integrate trapezoidal"

	Variable/G fp_pol=1,fp_minamp
	
	Variable/G tfp_pol,tfp_extent,tfp_mw,tfp_what
	
	tfp_mw= inf	// no maximum.

	String/G S_funcs
	String s
	s    =" line   2;"
	s += " poly  3 ;"
	s += " poly  4 ;"
	s += " poly  5 ;"
	s += " sin    4;"
	s += " dblexp 4;"
	s += " exp    3;"
	s += " lor    4;"
	s += " gauss  4;"

	S_funcs=s

	Variable/G fpks_extent

	String/G ar_wfit,ar_ow
	Variable/G ar_anno,ar_ex

	Variable/G afp_anno
	
	String/G pfr_title
		
	String/G S_Wave,new,none,calcbase  ///////////////
	new="_New_;"
	none="_None_;"
	calcbase="_Calculated_;"   ////////////////

	// Other global variables and waves
	Make/O/N=5 W_PM
	Variable/G V_tol=1e-10,V_peakX,V_peakP,V_start,V_theEnd,V_startX,V_theEndX

	String/G rb_b
	String/G fpks_b
	String/G afp_b
	String/G rb_ow
	Variable/G afp_pks
	Variable/G g_ptgn
	Variable/G g_atgn
	Variable/G/D V_npnts	// Otherwise Identify Peaks with Cursors fails.
	String/G fpks_weights
	Variable/G fpks_pktype
	Variable/G pfr_order,pfr_bi
	String/G pfr_sort
End

Macro MakePeaks()
	if( exists("g_w") == 0 )
		InitPeakMacros()
	endif
	MakeThosePeaks(,,,,,,)
End

Proc MakeThosePeaks(py,kind,noise,log,minx,maxx,pts)
	String py=g_w
	Variable kind=mp_kind,noise=mp_noise,log=mp_log,pts=mp_pnts,minx=mp_minx,maxx=mp_maxx
	Prompt py,"Output wave name"
	Prompt kind,"peak type",popup,"gaussian;lorentzian"
	Prompt noise,"peak wave noise (enoise val)"
	Prompt log,"x wave increment",popup,"linear;logarithmic"
	Prompt minx,"starting x value"
	Prompt maxx,"ending x value"
	Prompt pts,"points in output waves",popup,p_pts

	Silent 1;PauseUpdate
	mp_kind=kind;mp_noise=noise;mp_log=log;mp_pnts=pts;mp_minx=minx;mp_maxx=maxx
	String ctrX="W_MakeCentersX",ampY="W_MakeAmpsY",widthX="W_MakeWidthsX",px
	Variable kk1,kk3,last=strlen(py)-1,n=2^pts
	if(cmpstr("y",py[last])==0)
		last-=1
	endif
	px=py[0,last]+"X"
	Mk(py,n);Mk(px,n);$py=0
	SetScale x,minx,maxx,$py,$px
	if(log==2)
		$px=minx*(maxx/minx)^(p/n)
	else
		$px=x
	endif
	iterate (numpnts($ctrX))
		if(kind==1)	// gauss
			$py+=$ampY[i]*exp(-((($px[p]-$ctrX[i])/($widthX[i]/2/sqrt(ln(2))))^2))
		else
			if(kind==2) // lor
				kk3= $widthX[i]*$widthX[i]/4
				kk1=$ampY[i]*kk3
				$py+=kk1/(($px[p]-$ctrX[i])^2+kk3)
			endif
		endif
	loop
	$py+=enoise(noise)
	g_w=py
	if(log==2)
		g_wx=px
		AppWv(py,px,"")
	else
		g_wx=calcbase //////////
		AppWv(py,"","")
	endif
EndMacro

Macro InitializeMostEverything()
	if( exists("g_w") == 0 )
		InitPeakMacros()
	endif
	InitializeAlmostEverything(,,)
End

Proc InitializeAlmostEverything(w,wx,wb)
	String w=g_w,wx=g_wx,wb=g_b
	Prompt w,p_w,popup, WaveList("*",";","")
	Prompt wx,p_wx,popup,calcbase+WaveList("*", ";","")   ///////
	Prompt wb,p_b,popup, none+WaveList("*base*",";","")

	Silent 1;PauseUpdate
	g_w=w;g_wx=wx; SBs(wb)
	Redimension/D $w
	SameLen(w,"W_BaseRegion");W_BaseRegion=NaN
	SameLen(w,"W_AreaRegion");W_AreaRegion=NaN
	Mk("W_EstAmpsY",2);Mk("W_EstCentersX",2);Mk("W_EstCentersP",2);Mk("W_EstWidthsX",2)
	Mk("W_EstEdgesP",4);
	Mk("W_AreaNoBase",2);Mk("W_AreaX1",2);Mk("W_AreaX2",2);
	Mk("W_BasePM",2);Mk("W_PeakPM",2)
	if(exists(wx)==1)
		Redimension/D$wx
	endif
	if(exists(wb)==1)
		Redimension/D $wb
	endif
	AppWv(w,wx,"");HideFittedPeaks();HideAreaTags()
End

Macro InitBaseLineFit(w,wx)
	String w=g_w,wx=g_wx
	Prompt w,p_w,popup, WaveList("*",";","")
	Prompt wx,p_wx,popup,calcbase+WaveList("*", ";","") /////////////
	
	Silent 1;PauseUpdate
	g_w=w;g_wx=wx
	SameLen(w,"W_BaseRegion");AppWv("W_BaseRegion",wx,"");Modify mode(W_BaseRegion)=1;W_BaseRegion=NaN;AppCrs(w)
End

Macro AddRegionToFit()
	Silent 1;PauseUpdate
	CheckTwoCursors(g_w);W_BaseRegion[V_start,V_theEnd]=$g_w[p]
End

Macro DeleteRegionFromFit()
	Silent 1;PauseUpdate
	CheckTwoCursors(g_w);W_BaseRegion[V_start,V_theEnd]=NaN
End

Macro FitBaselineAtRegions(w,wx,wr,fit)
	String w=g_w,wx=g_wx,wr=fbar_wr,fit=fbar_fit
	Prompt w,p_w,popup, WaveList("*",";","")
	Prompt wx,p_wx,popup,calcbase+WaveList("*", ";","") //////////////
	Prompt wr,"region wave (baseline weighting wave)",popup, "W_BaseRegion" +";"+none
	Prompt fit,"baseline function",popup,S_funcs
	
	Silent 1;PauseUpdate
	g_w=w;g_wx=wx;fbar_wr=wr;fbar_fit=fit
	ChkLen(w,wx);ChkLen(w,wr)
	String ow="W_BaselineFit",wtmp="W_tmp",opts="";SameLen(w,ow)
	Variable t=0,typ=floor(strsearch(S_funcs,fit,0)/10)
	
	if(exists(wx)==1)
		opts+="/X="+wx
	endif
	if(exists(wr)==1)
		Dup(wr,wtmp)
		$wtmp=(numtype($wtmp)==0)
		opts+="/W="+wtmp
		ar_ex=3
	else
		ar_ex=1
	endif
	t=str2num(fit[7,8])-1
	String command="CurveFit "+fit[1,7]+", $w"+opts
	Execute command
	KillWv(wtmp)
	Mk("W_BaseCoefs",2);W_BaseCoefs={NaN,K0,K1,K2,K3,K4};Redimension/N=(2+t) W_BaseCoefs
	Mk("W_BasePM",2);W_BasePM={1,0,numpnts($w)-1,typ,t};Dup("W_BasePM","W_PM")
	if(exists(wx)==1)
		$ow=PolyMorph(W_BaseCoefs,$wx[p])
	else
		$ow=PolyMorph(W_BaseCoefs,x)
	endif
	AppWv(ow,wx,"")
	SBs(ow)
	ar_wfit=ow
End

Macro SubtractBaseline(w,wx,wb,ow)
	String w=g_w,wx=g_wx,wb=rb_b,ow=rb_ow
	Prompt w,p_w,popup, WaveList("*",";","")
	Prompt wx,p_wx,popup,calcbase+WaveList("*", ";","")  ///////////
	Prompt wb,p_b,popup, "_From Fit Peaks_;"+WaveList("*base*",";","")
	Prompt ow,"output baseline-corrected wave",popup,new+WaveList("*",";","")

	Silent 1;PauseUpdate
	g_w=w;g_wx=wx
	ChkLen(w,wx);ChkLen(w,wb)
	String cw,dcw,tcw="W_tmp0",tdcw="W_PM"
	if(exists(ow)!=1)
		NewWv(w,"NoBase");ow=S_Wave
	else
		SameLen(w,ow)
	endif
	rb_ow=ow
	if(exists(wb)!=1)
		cw="W_BaseCoefs"
		dcw="W_BasePM"
		if(numtype($dcw[1])!=0)
			Abort "Run Fit Peaks first!"
		endif
		Variable terms=$dcw[4],s=$dcw[1],e=$dcw[2],funcs=$dcw[0]
		Dup(cw,tcw);Redimension/N=(2+terms) $cw
		Dup(dcw,tdcw);$tdcw[0]=1
		wb="W_tmp";SameLen(w,wb);$wb=NaN
		if(exists(wx)==1)
			$wb[s,e]=PolyMorph($tcw,$wx[p])
		else
			$wb[s,e]=PolyMorph($tcw,x)
		endif
	endif

	$ow=$w-$wb
	AppWv(ow,wx,"");Modify zero(left)=1
	KillWv(wb);KillWv(tcw)
	g_w=ow
	SBs("_None_")
End

Macro InitAreaBetweenCursors(what,w,wx,wb,anno)
	Variable what=mir_what,anno=mir_anno
	String w=g_w,wx=g_wx,wb=g_b
	Prompt what,"integration type",popup, mir_pwh
	Prompt w,p_wd,popup, WaveList("*",";","")
	Prompt wx,p_wx,popup,calcbase+WaveList("*", ";","") //////////////
	Prompt wb,p_b,popup,none+WaveList("*base*", ";","")
	Prompt anno,"area annotation",popup,none+"_one tag on graph_;_one tag per peak_"
	
	Silent 1;PauseUpdate
	mir_what=what;g_w=w;g_wx=wx;SBs(wb);mir_anno=anno
	ChkLen(w,wx);ChkLen(w,wb);AppWv(w,wx,"");AppWv(wb,wx,"")
	SameLen(w,"W_AreaRegion");AppWv("W_AreaRegion",wx,"");Modify mode(W_AreaRegion)=1;W_AreaRegion=NaN;AppCrs(w)
	HideAreaTags();Mk("W_AreaNoBase",2);Mk("W_AreaX1",2);Mk("W_AreaX2",2)
End
Macro AreaBetweenCursors()
	Silent 1;PauseUpdate
	String w=g_w,wx=g_wx,wb=g_b,text,wv=w
	CheckTwoCursors(w)
	Variable/D s=V_start,en=V_theEnd,wi=0,oi,wlx,wrx
	if (mir_anno!=3)
		W_AreaRegion=NaN;HideAreaTags()
	endif
	W_AreaRegion[s,en]=$w[p]
	iterate (2)
		wlx=pnt2x($wv,s);wrx=pnt2x($wv,en)
		oi=wi
		if(mir_what==1)
			text="count"
			wi=(en-s+1)*mean($wv,wlx,wrx)
		else
		if(mir_what==2)
			text="rectangular area"
			if(exists(wx)==1)
				wi=AreaXYRPt($wx,$wv,s,en)
			else
				wi=(wrx-wlx)*mean($wv,wlx,wrx)
			endif
		else
			text="trapezoidal area"
			if(exists(wx)==1)
				wi=AreaXYTPt($wx,$wv,s,en)
			else
				wi=area($wv,wlx,wrx)
			endif
		endif
		endif
		Print wv+" "+text+" = "+num2str(wi)
		if(exists(wb)!=1)
			break
		endif
		wv=wb
	loop
	if(exists(wb)==1)
		Print w+" - "+wb+" "+text+" = "+num2str(oi-wi)
		wi=oi-wi;text += "\r(less baseline)"
	endif
	if(mir_anno!=1)
		Tag/C/N=$("area"+num2istr(g_atgn))/A=MB $w,(pnt2x($w,s)+pnt2x($w,en))/2, text+"\r= "+num2str(wi)
		if(mir_anno==3)
			g_atgn+=1
		endif
	endif
	InsertPoints 0,1,W_AreaNoBase,W_AreaX1,W_AreaX2;W_AreaNoBase[0]=wi;W_AreaX1[0]=V_startX;W_AreaX2[0]=V_theEndX
	Sort W_AreaX1,W_AreaX1,W_AreaX2,W_AreaNoBase
End

Macro HideAreaTags()
	do
		Tag/K/N=$("area"+num2istr(g_atgn))
		g_atgn-=1
	while (g_atgn>=0)
	g_atgn=0
End

Macro ListAreas()
	DoWindow/F AreaReportTable
	if(V_Flag==0)
		AreaReportTable()
	endif
End

Macro InitIdentifyPeaks(w,wx,wb,pol,box)
	String w=g_w,wx=g_wx,wb=g_b
	Variable box=g_bx,pol=mpr_pol
	Prompt w,p_w,popup, WaveList("*",";","")
	Prompt wx,p_wx,popup,calcbase+WaveList("*", ";","")  ///////////////
	Prompt wb,p_b,popup, none+WaveList("*base*",";","")
	Prompt pol,"peak polarity", popup,"positive;negative"
	Prompt box,p_bx
	
	Silent 1;PauseUpdate
	g_w=w;g_wx=wx;SBs(wb)
	box=abs(box);mpr_pol=pol;g_bx=box

	ChkLen(w,wx);ChkLen(w,wb)
	Mk("W_EstAmpsY",2);Mk("W_EstCentersX",2);Mk("W_EstCentersP",2);Mk("W_EstWidthsX",2)
	Mk("W_EstEdgesP",4)
	DoWindow/F PeakFitGraph
	if(V_Flag==0)
		AppWv(w,wx,"");DoWindow/C PeakFitGraph
	endif
	AppWv(w,wx,"");AppWv(wb,wx,"")
	ShowEdges();AppWv("W_EstAmpsY","W_EstCentersX","");Modify mode(W_EstAmpsY)=8,marker(W_EstAmpsY)=18
	AppCrs(w);HideFittedPeaks()
End

Macro IdentifyPeakWithCursors()
	Silent 1;PauseUpdate
	String w=g_w,wx=g_wx,wb=g_b
	CheckTwoCursors(w)
	Variable s=V_start,en=V_theEnd,npks,ctr,box=g_bx,pol=mpr_pol
	String wtmp="W_tmp",wtmp2="W_tmp2",e0="W_tmp0",e1="W_tmp1"
	SplitEdges(e0,e1);npks=V_npnts
	if(exists(wx)==1)
		box *= sign($wx[en]-$wx[s])
		FindPX(w,wx,box+$wx[s]);box=V_peakP
		FindPX(w,wx,$wx[s]);box=abs(V_peakP-box)
	else
		box*=sign(rightx($w)-leftx($w))
		box=x2pnt($w,box+leftx($w))
	endif
	box=max(1,box);Print "Averaging ",box," points..."
	Duplicate/O/R=[s,en] $w,$wtmp;$wtmp=MyBoxSmooth($w,p+s,box)
	WaveStats/Q $wtmp // Find peak's center pt
	if(pol==1)
		ctr=V_maxLoc
	else
		ctr=V_minLoc
	endif
	ctr=x2pnt($wtmp,ctr)+s
	InsertPoints 0,1,W_EstCentersP,W_EstCenters,W_EstAmpsY,W_EstWidthsX,$e0,$e1
	InsertPoints 0,2,W_EstEdgesP
	W_EstCentersP[0]=ctr
	$e0[0]=s;$e1[0]=en
	W_EstAmpsY[0]=$wtmp[ctr-s]
	KillWv(wtmp);KillWv(wtmp2)
	if(exists(wb)==1)
		W_EstAmpsY[0]-=MyBoxSmooth($wb,ctr,box)
	endif
	if(exists(wx)==1)
		W_EstCenters[0]=$wx[ctr]
	else
		W_EstCenters[0]=pnt2x($w,ctr)
	endif
	W_EstWidthsX[0]=abs(V_theEndX-V_startX)
	MergeEdges(e0,e1);ShowEdges()
End

Macro DeletePksBetweenCursors()
	Silent 1;PauseUpdate
	String w=g_w,wx=g_wx,wb=g_b
	CheckTwoCursors(w)
	Variable s=V_start,en=V_theEnd,npks,ctr,box=g_bx,pol=mpr_pol
	String wtmp="W_tmp",wtmp2="W_tmp2",e0="W_tmp0",e1="W_tmp1"
	SplitEdges(e0,e1);npks=V_npnts
	iterate (npks)
		if((W_EstCentersP[i] >= s)%&(W_EstCentersP[i] <= en))
			W_EstCentersP[i]=NaN;W_EstCenters[i]=NaN;W_EstAmpsY[i]=NaN;W_EstWidthsX[i]=NaN
			$e0[i]=NaN;$e1[i]=NaN		
		endif
	loop
	MergeEdges(e0,e1);ShowEdges()
	if(V_npnts<1)
		DoAlert 0, "no peaks left"
	endif
End

// Find peaks by analyzing peak data's smoothed derivatives (WARNING: sensitive to noise)
Macro AutoIdentifyPeaks()
	if( exists("FindAPeak") == 0 )
		String msg="Please install a shortcut or alias to the \"FindPeaks\" XOP into your Igor Extensions folder."
		msg += " You will find the \"FindPeaks\" XOP in the \"More Extensions:Data Analysis:\" folder."
		Abort msg
	endif
	AutoIdentifyThePeaks(,,,,,,,,)
End

// Find peaks by analyzing peak data's smoothed derivatives (WARNING: sensitive to noise)
Proc AutoIdentifyThePeaks(what,w,wx,wb,mina,pol,box,extent,maxWidth)
	String w=g_w,wx=g_wx,wb=g_b
	Variable mina=fp_minamp
	Variable what=tfp_what,pol=tfp_pol,box=g_bx,extent=tfp_extent,maxWidth=tfp_mw
	Prompt what,"do what with found peaks?",popup, "ADD to existing peak list;OVERWRITE peak list"
	Prompt w,p_w,popup,WaveList("*", ";","")
	Prompt wx,p_wx,popup,calcbase+WaveList("*", ";","")  //////////////
	Prompt wb,p_b,popup,none+WaveList("*base*", ";","")
	Prompt mina,"min pk amplitude, or 0 for auto-level"
	Prompt pol,"peak polarity", popup,"positive;negative" 
	Prompt box,p_bx
	Prompt extent,"extent of Y wave to search",popup,"entire wave;between cursors"
	Prompt maxWidth,"maximum peak x width (INF for no max)"

	Silent 1;PauseUpdate
	tfp_what=what;g_w=w;g_wx=wx;SBs(wb)

	if( what == 1 )	// append
		if( exists("W_EstCentersX") == 0 )	// but waves don't exist!
			what= 2	// overwrite
		endif
	endif

	box=abs(box)
	tfp_pol=pol;g_bx=box;tfp_extent=extent
	String lvlw="W_ShowMinPkLvl"
	ChkLen(w,wx);ChkLen(w,wb)
	Variable s=0,en=numpnts($w)-1
	if(extent==2) // use cursors
		CheckTwoCursors(w)
		s=V_start
		en=V_theEnd
	endif
	if(exists(wx)==1)
		box *= sign($wx[en]-$wx[s])
		FindPX(w,wx,box+$wx[s]);box=V_peakP
		FindPX(w,wx,$wx[s]);box=abs(V_peakP-box)
	else
		box*=sign(rightx($w)-leftx($w))
		box=x2pnt($w,box+leftx($w))
	endif
	box=max(1,box);Print "Averaging ",box," points..."
	if(mina==0)
		Dup(w,lvlw)
		if(exists(wb)==1)
			$lvlw-=$wb
		endif
		Smooth 3, $lvlw
		WaveStats/Q/R=[s,en] $lvlw
		mina=V_sdev*0.5 //Heuristic
		if (pol==2)
			mina *= -1
		endif
		Print "Auto min peak amplitude = ",mina
	endif
	fp_minamp=mina
	SameLen(w,lvlw);$lvlw=NaN
	if(exists(wb)==1)
		$lvlw[s,en]=mina+$wb[p]
	else
		$lvlw[s,en]=mina
	endif
	AppWv(lvlw,wx,"");ResumeUpdate
	PauseUpdate

	if(what==2)
		Mk("W_EstAmpsY",2)
		Mk("W_EstCentersX",2)
		Mk("W_EstCentersP",2)
		Mk("W_EstWidthsX",2)
		Mk("W_EstEdgesP",4)
	endif

	Variable lf=s,rt=en,fl,ce,dX,npks=0
	do
		if( exists(wb) == 1 )
			FindAPeak/B=$wb  mina,pol,box,$w[rt,lf]	//r=>l search
		else
			FindAPeak  mina,pol,box,$w[rt,lf]	//r=>l search
		endif

		if(V_Flag==0)
			// Insert a new Peak
			InsertPoints 0,1,W_EstAmpsY,W_EstCentersX,W_EstCentersP,W_EstWidthsX
			InsertPoints 0,2,W_EstEdgesP
			W_EstCentersP[0]=V_peakP
			fl=floor(V_peakP);ce=ceil(V_peakP)
			if(exists(wx)==1)	// convert f.p. V_peakP f.p. to wx coordinate
				if(fl != ce)
					dX=$wx[ce]-$wx[fl]
					W_EstCentersX[0]=$wx[fl]+(V_peakP-fl)*dX
				else
					W_EstCentersX[0]=$wx[V_peakP]
				endif
			else
				W_EstCentersX[0]=V_peakX
			endif
			rt=fl-floor((box-1)/2)
			npks+=1
			Printf "%d...",npks
		endif
	while ((V_Flag==0) %& ((rt-lf)>1))
	AppWv("W_EstAmpsY","W_EstCentersX","");Modify mode(W_EstAmpsY)=8,marker(W_EstAmpsY)=18;ResumeUpdate
	PauseUpdate
	Printf "\r%d Peaks found...",npks
	if(npks>0)
		Printf "Estimating their sizes..."
		if(exists(wx)==1)
			if(exists(wb)==1)
				EstimatePeakSizes/X=$wx/B=$wb/E=W_EstEdgesP 50,maxWidth,box,npks,W_EstCentersP,$w,W_EstAmpsY,W_EstWidthsX
			else
				EstimatePeakSizes/X=$wx/E=W_EstEdgesP 50,maxWidth,box,npks,W_EstCentersP,$w,W_EstAmpsY,W_EstWidthsX
			endif
		else
			if(exists(wb)==1)
				EstimatePeakSizes/B=$wb/E=W_EstEdgesP 50,maxWidth,box,npks,W_EstCentersP,$w,W_EstAmpsY,W_EstWidthsX
			else
				EstimatePeakSizes/E=W_EstEdgesP 50,maxWidth,box,npks,W_EstCentersP,$w,W_EstAmpsY,W_EstWidthsX
			endif
		endif
	endif
	String e0="W_tmp0",e1="W_tmp1"
	SplitEdges(e0,e1);MergeEdges(e0,e1)
	Print "Done"
	ShowEdges()
End

Proc SplitEdges(e0,e1)
	String e0,e1
	Variable pts=numpnts(W_EstCentersP)
	Make/O/D/N=(pts) $e0,$e1
	$e0=W_EstEdgesP[0+p*2]
	$e1=W_EstEdgesP[1+p*2]
	WaveStats/Q W_EstCentersX // Count non-Nan peaks
End

Proc MergeEdges(e0,e1)
	String e0,e1
	Sort W_EstCentersP,W_EstCentersP,W_EstCentersX,W_EstAmpsY,W_EstWidthsX,$e0,$e1	// NaN's sorted to end of array
	W_EstEdgesP[0,;2]=$e0[p/2]
	W_EstEdgesP[1,;2]=$e1[(p-1)/2]
	KillWaves $e0,$e1
	WaveStats/Q W_EstCentersX
	Redimension/N=(V_npnts+2) W_EstCentersP,W_EstCentersX,W_EstAmpsY,W_EstWidthsX
	Redimension/N=((V_npnts+2)*2) W_EstEdgesP
End

Proc ShowEdges()
	Variable b,x0,x1,p0,p1
	String sx="W_ShowEstWidthsX",sy="W_ShowEstWidthsY"
	WaveStats/Q W_EstAmpsY
	Mk(sx,V_npnts*5+2);Mk(sy,V_npnts*5+2)
	iterate (V_npnts)
		b=i*5;p0=W_EstEdgesP[2*i];p1=W_EstEdgesP[2*i+1]
		if(exists(g_wx)==1)
			x0 =$g_wx[p0]
			x1=$g_wx[p1]
		else
			x0=pnt2x($g_w,p0)
			x1=pnt2x($g_w,p1)
		endif
		$sx[b,b+1]=x0
		$sx[b+2,b+3]=x1
		$sy[b,b+3]=0
		$sy[b+1,b+2]=W_EstAmpsY[i]
		if(exists(g_b)==1)
			$sy[b,b+1]+=$g_b[p0]
			$sy[b+2,b+3]+=$g_b[p1]
		endif
	loop
	AppWv(sy,sx,"")
End

Function/D MyBoxSmooth(w,pp,box)
	Wave/D w
	Variable pp,box // box should be odd
	
	Variable/D result=0
	Variable  hp,top=numpnts(w)-1,terms
	box = 2*trunc(box/2)+1
	pp-=trunc(box/2)
	hp=limit(0,pp+box-1,top)
	pp=limit(0,pp,top)
	terms=hp-pp+1
	do
		result+=w[pp]
		pp+=1
	while (pp<=hp)
	return result/terms
End


Function AreaXYTPt(wx,wy,pa,pb)
	Wave wx,wy	// monotonic wx!
	Variable pa,pb	//pa<pb

	Variable a=0
	do
		a+=(wy[pa]+wy[pa+1])*(wx[pa+1]-wx[pa])/2
		pa+=1
	while(pa<pb)
	return a
End

Function AreaXYRPt(wx,wy,pa,pb)
	Wave wx,wy	// monotonic wx!
	Variable pa,pb	// pa<pb

	Variable a=0
	do
		a+=wy[pa]*(wx[pa+1]-wx[pa])
		pa+=1
	while(pa<pb)
	return a
End

// W_PM[0]=n, number of functions to fit
// W_PM[1]=start p of peak fit
// W_PM[2]=end p of peak fit
// W_PM[3]=tp, type of baseline
// W_PM[4]=nw, baseline terms, could be 0
// W_PM[5]=tp, type of 1st peak function
//...
// w[0]=NaN
// w[1]=K0
// w[2]=1st  coefficient...
Function/D PolyMorph(w,xx)
	Wave/D w
	Variable/D xx
	
	WAVE W_PM=root:W_PM
	Variable/D tp,nw, bs=3,n=W_PM[0],st=2,ii,s,r=w[1] // K0
 	if(n==0)
		return r
	endif
	do	
		tp=W_PM[bs]
		nw=W_PM[bs+1]
		if(nw>0)
			if(tp<4) // line, poly 3 - poly 5 K1*x+K2*x^2...
				ii=nw-1
				s=0
				do
					s=w[ii+st]+xx*s
					ii-=1
				while (ii >= 0)
				r+= s*xx
			else
				if(tp==8) // gauss
					r+=w[st]*exp(-((xx-w[st+1])/w[st+2])^2)
				else
				if(tp==7) // lor
					r+=w[st]/((xx-w[st+1])^2+w[st+2])
				else
				if(tp==6) // exp
					r+= w[st]*exp(-w[st+1]*xx)
				else
				if(tp==5)	// dblexp
					r+= w[st]*exp(-w[st+1]*xx)*exp(-w[st+2]*xx)
				else	// 4, sin
					r+= w[st]*sin(w[st+1]*xx+w[st+2])
				endif
				endif
				endif
				endif
			endif
			st+=nw	
		endif
		bs+=2
		n-=1
	while (n>0)
	return r
End

// This fits peak equations and a baseline equation or wave to the data.
// If the baseline is removed, a constant baseline + peaks will be fit to the data.
Macro FitPeaksAndBaseline(w,wx,extent,wb,pktype,wts)
	String w=g_w,wx=g_wx,wb=fpks_b,wts=fpks_weights
	Variable extent=fpks_extent	,pktype=fpks_pktype
	Prompt w,p_w,popup, WaveList("*",";","")
	Prompt wx,p_wx,popup,calcbase+WaveList("*", ";","")  ///////////////////
	Prompt extent,"extent of peak wave to fit",popup,"entire wave, all peaks;wave & peaks within cursors"
	Prompt wb,p_b,popup, none+S_funcs+"_From Fit Baseline_;"+WaveList("*base*",";","")
	Prompt pktype,"type of peak fit",popup, S_funcs[70,89] // lor 4,gauss 4
	Prompt wts,"weighting wave",popup, none+WaveList("*",";","")
	
	Silent 1;PauseUpdate
	g_w=w;g_wx=wx;SBs(wb);fpks_pktype=pktype
	fpks_extent=extent;fpks_weights=wts
	
	ChkLen(w,wx);ChkLen(w,wb);ChkLen(w,wts)
	ChkLen("W_EstCentersX","W_EstAmpsY");ChkLen("W_EstCentersX","W_EstAmpsY")
	Variable st=0,en=numpnts($w)-1
	if(extent==2)
		CheckTwoCursors(w)
		st=V_start
		en=V_theEnd
	endif

	String tmpw="W_tmp",ow="W_PeakFit"
	Mk("W_PeakFitCoefs",5);Mk("W_PeakPM",5)
	Dup(w,ow);$ow=NaN
	Variable terms=0,pkTyp=0,lim=1,baseTyp = floor(strsearch(S_funcs,wb,0)/10)
	String cmd,opts=" /D="+ow
	WaveStats/Q/R=[st,en] $w
	Variable/D bsgs=V_avg,xdiff=pnt2x($w,en)-pnt2x($w,st)
	if(bsgs==0)
		bsgs=V_tol
	endif
	if(exists(wx)==1)
		opts+=" /X="+wx
		xdiff=$wx[en]-$wx[st]
	endif
	if(exists(wts)==1)
		opts+="/W="+wts
	endif
	AppWv(ow,wx,"")
	// initial baseline guesses
	String hold
	if(baseTyp>=0)// function selected (expect "_poly   1",  "_dblexp 5", etc)
		terms=str2num(wb[7,8])-1 // exclude K0
		cmd="CurveFit/Q/O "+wb[1,7]+", "+w+"[st,en] "+opts
		Execute cmd
		W_PeakFitCoefs={NaN,K0,K1,K2,K3,K4}
		W_PeakPM={1,st,en,baseTyp,terms}
		hold="100000"[0,terms+1]
		afp_b="_From Fit Peaks_"
	else
		if(exists(wb)==1)	// remove base before fit
			Dup(w,tmpw);w=tmpw;$w-=$wb
			terms=0
			W_PeakFitCoefs={NaN,0}
			W_PeakPM={1,st,en,0,0}
			hold="11"
			bsgs=0 // K0 not varied
		else
			afp_b="_From Fit Peaks_"
			if(cmpstr(wb,"_From Fit Baseline_")==0)	// Use base fit to start
				Dup("W_BaseCoefs","W_PeakFitCoefs")
				Dup("W_BasePM","W_PeakPM")
				terms=W_PeakPM[4]
				W_PeakPM[1]=st
				W_PeakPM[2]=en
				hold="100000"[0,terms+1]
			else // _None_
				terms=0
				W_PeakFitCoefs={NaN,bsgs}
				W_PeakPM={1,st,en,0,0}
				hold="10" //  K0 varied
			endif
		endif
	endif
	// prevent singular matrix w/ poly base
	if(W_PeakPM[3]<4)
		W_PeakFitCoefs[1,1+terms]=W_PeakFitCoefs[p]*(W_PeakFitCoefs[p]!=0)+(W_PeakFitCoefs[p]==0)*bsgs/(xdiff^(p-1))
	endif
	WaveStats/Q W_EstCentersX
	Variable npks=V_npnts
	if(npks<1)
		Abort "no peaks in peak center wave "+"W_EstCentersX"
	endif
	Variable/D startOfTerms=2+terms,dst,sz=3,dsz=2,pks=0,kk0
	Redimension/N=(5+npks*2)W_PeakPM
	iterate (npks)
		if((W_EstCentersP[i]>=st) %& (W_EstCentersP[i]<=en))
			dst=3+2*W_PeakPM[0]
			Redimension/N=(startOfTerms+sz) W_PeakFitCoefs
			if(pktype==2)	// gauss
				W_PeakFitCoefs[startOfTerms]=W_EstAmpsY[i]						// K1 
				W_PeakFitCoefs[startOfTerms+1]=W_EstCentersX[i]					// K2
				W_PeakFitCoefs[startOfTerms+2]=W_EstWidthsX[i]/2/sqrt(ln(2)) 	// K3
				pkTyp=8
			else 	// lor
				W_PeakFitCoefs[startOfTerms+2]=W_EstWidthsX[i]^2 / 4 			// K3 
				W_PeakFitCoefs[startOfTerms]=W_EstAmpsY[i]*(W_EstWidthsX[i]^2)/4	// K1
				W_PeakFitCoefs[startOfTerms+1]=W_EstCentersX[i]					// K2
				pkTyp=7
			endif
			W_PeakPM[dst]=pkTyp
			W_PeakPM[dst+1]=3
			hold+= "000"
			startOfTerms +=sz
			dst+=dsz
			W_PeakPM[0]+=1
			pks+=1
		endif
	loop
	if(pks<1)
		Abort "no peaks in between cursors"
	else
		Print "Fitting ",pks," peaks"
	endif
	// Fit peaks using initial guesses for coefficients
	Dup("W_PeakPM","W_PM")
	sprintf cmd, "FuncFit/H=\"%s\" PolyMorph %s, %s[%g,%g]%s", hold, "W_PeakFitCoefs", w, st, en, opts
Print cmd
	ResumeUpdate
	Execute cmd		// FuncFitÉ ("singular matrixÉ" means normalize x values!)
	KillWv(tmpw)
	if(exists(wb)==1)
		$ow[st,en]+=$wb[p]
	endif
	ar_wfit=ow
	ar_ex=extent
End

Macro Report(title,srt,order,basInf)
	String title=pfr_title, srt=pfr_sort,
	Variable order=pfr_order,basInf=pfr_bi
	Prompt title,"report title"
	Prompt srt,"sort report by",popup,"center;amplitude;width;area"
	Prompt order,"sort order",popup,"ascending;descending"
	Prompt basInf,"baseline info",popup,"include;omit"
	
	Silent 1
	pfr_title=title;pfr_sort=srt;pfr_order=order;pfr_bi=basInf
	String center="W_PkCenters",amplitude="W_PkAmps",width="W_PkFWHMs",area="W_PkAreas"

	Variable terms=W_PeakPM[4],btyp=W_PeakPM[3] // baseline terms, not counting K0
	Variable sz=3,dsz=2,st=2+terms,dst=5
	Variable npks=W_PeakPM[0]-1,np1=npks-1,nrows=max(2,npks)
	String text=num2istr(npks)
	Mk(center,nrows);Mk(amplitude,nrows);Mk(width,nrows);Mk(area,nrows)
	if(npks<1)
		Abort "no peaks"
	endif
	DoWindow/F PeakFitGraph
	if(V_Flag==0)
		AppWv(g_w,g_wx,"");AppWv("W_PeakFit",wx,"");ColorWaves();DoWindow/C PeakFitGraph
	else
		NoCrs()
	endif
	
	DoWindow/F PeakReportTable
	if(V_Flag==0)
		PeakReportTable()
	endif

	Variable pkTyp=W_PeakPM[5]
	// Separate the coefficients into separate waves
	$center[0,np1]=W_PeakFitCoefs[st+1+p*sz] // K2 is same
	if(pkTyp==8)	// Gaussian
		$amplitude[0,np1]=W_PeakFitCoefs[st+p*sz]							// K1
		$width[0,np1]=2* sqrt(ln(2))*W_PeakFitCoefs[st+2+p*sz]			// 2*sqrt(ln(2))*K3
		$area[0,np1]=sqrt(Pi)*W_PeakFitCoefs[st+p*sz]*W_PeakFitCoefs[st+2+p*sz]	//  K1*K3*sqrt(Pi)
		text+=" Gaussian peaks"
	else
		if(pkTyp==7) // Lorentzian
			$amplitude[0,np1]=W_PeakFitCoefs[st+p*3] / W_PeakFitCoefs[st+2+p*3]	// K1/K3
			$width[0,np1]=2*sqrt(W_PeakFitCoefs[st+2+p*3])				// 2* sqrt(K3) 
			$area[0,np1]=Pi*W_PeakFitCoefs[st+p*3]/sqrt(W_PeakFitCoefs[st+2+p*3])// Pi*K1/sqrt(K3)
			text+=" Lorentzian peaks"
		else
			Abort "unsupported peak type"
		endif
	endif
	String cmd="Sort "
	if(order==2)
		cmd+="/R  "
	endif
	cmd+="$"+srt+",$center,$amplitude,$width,$area"
	
	Execute cmd		// Sort...
	
	DoWindow/F PeakReportLayout
	if(V_Flag==0)
		PeakReportLayout()
	endif
	if((basInf==1)%&(terms>0)%&(btyp>=0))
		text+="\r\rBaseline function: "+S_funcs[btyp*10+1,btyp*10+7]
		text+="\r\tK0 = "+num2str(W_PeakFitCoefs[1])
		iterate (terms)
			text+="\r\tK"+num2istr(i+1)+"= "+num2str(W_PeakFitCoefs[2+i])
		loop
	endif
	ReplaceText/N=info text
	ReplaceText/N=title "\JC"+title
	Modify rows(PeakReportTable)=nrows
End

Macro ShowFittedPeaks(w,wx,wb,pks,anno)
	String w=g_w,wx=g_wx,wb=afp_b
	Variable pks=afp_pks,anno=afp_anno
	Prompt w,p_w,popup, WaveList("*",";","")
	Prompt wx,p_wx,popup,calcbase+WaveList("*", ";","")   ///////////////
	Prompt wb,p_b,popup, none+"_From Fit Peaks_;"+WaveList("*base*",";","")
	Prompt pks,"peak waves",popup,"off;on"
	Prompt anno,"peak annotation tags",popup, "off;on"

	Silent 1;PauseUpdate
	g_w=w;g_wx=wx;afp_anno=anno;afp_pks=pks;SBs(wb)
	ChkLen(w,wx);ChkLen(w,wb)
	String tdcw="W_PM",tcw="W_tmp",pkw,text,typ
	if ((pks==1)%&(anno==1))
		DoAlert 0,"Both tags and peaks are off!"
	endif
	if(numtype(W_PeakPM[1])!=0)
		Abort "Run FitPeaks first!"
	endif
	Variable s=W_PeakPM[1],en=W_PeakPM[2],npks=W_PeakPM[0]-1,terms,ctr,amp,width,area
	Dup("W_PeakPM",tdcw);Dup("W_PeakFitCoefs",tcw)
	if(cmpstr(wb,"_From Fit Peaks_") == 0)
		terms=W_PeakPM[4]
		Redimension/N=(2+terms+3) $tcw
		Redimension/N=(5+2) $tdcw;$tdcw[0]=2
	else
		$tdcw={2,s,en,0,0,NaN,3}
		$tcw={NaN,0,8,8,8}
		terms = 0
	endif
	Variable sz=3,dsz=2,st=2+terms,cst=2+W_PeakPM[4],dcst=5
	iterate (npks)
		$tcw[st,st+2]=W_PeakFitCoefs[cst+p-st]
		$tdcw[5,6]=W_PeakPM[dcst+p-5]
 		if (pks==2)
			NewWv(w,"_fpk");pkw=S_Wave;$pkw=NaN
			if(exists(wx)==1)
				$pkw[s,en]=PolyMorph($tcw,$wx[p])
			else
				$pkw[s,en]=PolyMorph($tcw,x)
			endif
			if(exists(wb)==1)
				$pkw[s,en]+=$wb[p]
			endif
			AppWv(pkw,wx,"");Modify mode($pkw)=1
		endif
		ctr=$tcw[st+1] // K2
		if(W_PeakPM[5]==8)
			typ="\JCGaussian\r"
			amp=$tcw[st]	//K1
			width=2*sqrt(ln(2))*$tcw[st+2]
			area=sqrt(Pi)*amp*$tcw[st+2]
		else
			typ="\JCLorentzian\r"
			amp=$tcw[st]/$tcw[st+2]	//K1/K3
			width=2*sqrt($tcw[st+2])
			area=Pi*$tcw[st]/sqrt($tcw[st+2])
		endif
		sprintf text, "\JL\Z09amp=%g\rcenter=%g\rwidth=%g\rarea=%g",amp,ctr,width,area
		 if(anno==2)
			FindPX(w,wx,ctr)
			Tag/C/N=$("peak"+num2istr(g_ptgn))/A=MB  $w, V_peakX, typ+text
			g_ptgn+=1
		endif
		dcst+=dsz;cst+=sz
	loop
End

Macro HideFittedPeaks()
	Silent 1;PauseUpdate
	String w,ws=WaveList("*_fpk",";","WIN:"+WinName(0,1))
	Variable pos
	do
		pos=strsearch(ws,";",0)
		if(pos<0)
			break
		endif
		w=ws[0,pos-1]
		ws[0,pos]=""
		RmWv(w);KillWv(w)
	while (1)
	do
		Tag/K/N=$("peak"+num2istr(g_ptgn))
		g_ptgn-=1
	while (g_ptgn>=0)
	g_ptgn=0
	ColorWaves()
End

Proc FindPX(w,wx,xx)
	String w,wx	// wx monotonic
	Variable xx

	if(exists(wx)==1)
		FindLevel/Q $wx, xx 
		V_peakP = x2pnt($wx,V_levelX)
		if(V_Flag==1)
			if((xx<$wx[1])%&($wx[0]<$wx[1]))
				V_peakP=0
			else
				V_peakP=numpnts($w)-1
			endif		
		endif
		V_peakX = pnt2x($w,V_peakP)
	else
		V_peakP = x2pnt($w,xx)
		V_peakX = xx
	endif
End

Proc ChkLen(w1,w2)
	String w1,w2	// Wave names
	if((exists(w1)==1) %& (exists(w2)==1))
		if(numpnts($w1) != numpnts($w2))
			Abort w1+" and "+w2+" have different number of points!"
		endif
	endif
End

Proc Mk(w,n)
	String w
	Variable n
	if(exists(w)==1)
		Redimension/D/N=(n) $w;$w=NaN
	else
		Make/D/N=(n) $w=NaN
	endif	
End

Proc Dup(t,w)
	String t,w
	String n=""
	if (exists(w)==1)
		n=note($w) // keep w's note
	endif
	Duplicate/O $t,$w
	Note/K $w;Note $w,n
End	

Proc SameLen(t,w)
	String t,w
	if(exists(w)!=1)
		Duplicate/O $t,$w
	else
		Redimension/D/N=(numpnts($t)) $w;CopyScales $t,$w
	endif
End	

Proc NewWv(srcw,sfx)
	String srcw,sfx

	Variable ii=1
	S_Wave=srcw[0,17-strlen(sfx)]+sfx
	if(exists(S_Wave)==1)
		String shw=S_Wave
		do
			S_Wave=srcw[0,14-strlen(sfx)]+num2istr(ii)+sfx
			ii+=1
		while ((exists(S_Wave)==1)*(ii<999))
	endif
	if(exists(srcw)==1)
		Duplicate/O $srcw,$S_Wave
	else
		Make/D/O $S_Wave
	endif
	Print "Created wave "+S_Wave
End

Proc RmWv(w)
	String w

	CheckDisplayed/W=$WinName(0,1) $w
	if(V_Flag==1)
		DoWindow/F $WinName(0,1);Remove $w
	endif
End

Proc KillWv(w)
	String w
	if(exists(w)==1)
		CheckDisplayed/A $w
		if(V_Flag==0)
			KillWaves $w
		endif
	endif
End

Proc AppCrs(w)
	String w
	ShowInfo;Cursor/P A,$w,0;Cursor/P B,$w,numpnts($w)-1
End
Proc NoCrs()
	HideInfo;Cursor/K A;Cursor/K B
End
Proc AppWv(w,wx,ax)
	String w,wx,ax
	
	String wn=WinName(0,1)
	if(strlen(wn)==0)
		ax = "Display $w"
	else
		DoWindow/F $wn
		ax="Append"+ax+" $w "
	endif
	if(exists(w)==1)
		CheckDisplayed $w
		if(V_Flag==0)
			if(exists(wx)==1)
				ax+=" vs $wx"
			endif
			Execute ax
			ColorWaves()
		endif
	endif
End
Proc SBs(s)
	String s
	if((exists(s)==1)+ (cmpstr(s+";",none)==0))
		g_b=s;rb_b=s;fpks_b=s;afp_b=s
	else
	if(strsearch(S_funcs,s,0)>=0)
		fpks_b=s
	else
	if(cmpstr(s,"_From Fit Peaks_")==0)
		rb_b=s;afp_b=s
	else
	if(cmpstr(s,"_From Fit Baseline_")==0)
		fpks_b=s
	endif
	endif
	endif
	endif
End
Macro ShowFitResidual(w,wx,wfit,ow,extent,anno)
	String w=g_w,wx=g_wx,wfit=ar_wfit,ow=ar_ow
	Variable anno=ar_anno,extent=ar_ex
	Prompt w,p_wd,popup, WaveList("*",";","")
	Prompt wx,p_wxd,popup,calcbase+WaveList("*", ";","")  ///////////////
	Prompt wfit,"fit wave",popup, WaveList("*Fit",";","")
	Prompt ow,"output residual wave",popup,new+WaveList("*",";","")
	Prompt extent,"extent of residual wave to measure",popup,"entire wave;between cursors;"+"W_BaseRegion"
	Prompt anno,"standard deviation lines annotations",popup,"± 1; ± 2;"+none

	Silent 1;PauseUpdate
	String wf=wfit[2,5]
	if(exists(ow)!=1)
		NewWv(w,wf+"Res");ow=S_Wave
	else
		SameLen(w,ow)
	endif
	g_w=w;g_wx=wx;ar_wfit=wfit;ar_ow=ow;ar_ex=extent;ar_anno=anno
	ChkLen(w,wx);ChkLen(w,wfit)
	Variable s=0,en=numpnts($w)-1
	if(extent==2) // use cursors
		CheckTwoCursors(w)
		s=V_start
		en =V_theEnd
	endif
	if(extent==3)
		if(exists("W_BaseRegion")!=1)
			Abort "W_BaseRegion"+" doesn't exist!"
		endif
		ChkLen(w,"W_BaseRegion")
		$ow=($w[p]-$wfit[p]) / (numtype(W_BaseRegion[p])==0)
	else
		$ow=NaN;$ow[s,en]=$w[p]-$wfit[p]
	endif
	
	WaveStats/R=[s,en] $ow
	String rx="W_Resid"+wf+"X",ry="W_Resid"+wf+"Y"
	if(anno<3)
		Mk(rx,6);Mk(ry,6)
		if(exists(wx)==1)
			$rx[0,;4]=$wx[s]
			$rx[1,;4]=$wx[en]
		else
			$rx[0,;4]=pnt2x($w,s)
			$rx[1,;4]=pnt2x($w,en)
		endif
		$ry[0,;4]=(trunc(p/2)-1)*V_sdev*anno
		$ry[1,;4]=$ry[p-1]
		AppWv(ry,rx,"/R")
		Variable maxv=max(V_max,anno*V_sdev),minv=min(V_min,-anno*V_sdev)
		Modify zero(right)=1,minor(right)=1
		SetAxis right minv-3*(maxv-minv),maxv
	endif
	AppWv(ow,wx,"/R");Modify mode($ow)=3,marker($ow)=8,msize($ow)=1
	Label right "Residual (fit-data)"
End

Macro ShowOnlyDataAndBase()
	Silent 1;PauseUpdate
	String keep=g_w+";"+g_wx+";"
	keep+=g_b+";"+WaveList("*",";","WIN:"+WinName(0,1))
	keep+=g_keep+";";g_keep=""
	String w,wn=WinName(0,1),ws=WaveList("*",";","WIN:"+wn)
	Variable pos,flg=0
	if(strlen(wn)!=0)
		do
			pos=strsearch(ws,";",0)
			if(pos<0)
				break
			endif
			w=ws[0,pos-1]
			ws[0,pos]=""
			if(strsearch(keep,w+";",0)<0)
				RmWv(w);flg=1
			endif
		while (1)
		if(flg)
			ColorWaves()
		endif
	endif
End

Macro ColorWaves() : GraphStyle
	Silent 1;PauseUpdate
	Modify/Z lSmooth=1,lHair=0.5,minor(bottom)=1
	Modify/Z rgb[0]=(65535,0,0)
	Modify/Z rgb[1]=(0,0,65535),rgb[2]=(3009,65535,1882),rgb[3]=(0,0,0)
	Modify/Z rgb[4]=(0,0,65535),rgb[5]=(63953,3661,65535)
	Modify/Z rgb[6]=(37510,1,1),rgb[7]=(27232,40528,22540),rgb[8]=(1531,1314,28456)
	Legend/C/N=default ""
End


Proc CheckTwoCursors(w)
	String w	// the y wave
	Variable/G V_start,V_theEnd,V_startX,V_theEndX // point indices, x vals
	String wa=CsrWave(A),wb=CsrWave(B),t="cursors on target graph "
	if((strlen(wa)==0) + (strlen(wb)==0))
		Abort "expecting two "+t
	endif
	V_start=pcsr(A);V_theEnd=pcsr(B)
	V_startX=hcsr(A);V_theEndX=hcsr(B)
	Variable x0=pnt2x($w,0),x1=pnt2x($w,1),ok
	ok = (numpnts($wa)==numpnts($w)) %& (numpnts($wb)==numpnts($w)) 
	ok = ok %& (x0==pnt2x($wa,0)) %& (x1==pnt2x($wa,1))
	ok = ok %& (x0==pnt2x($wb,0)) %& (x1==pnt2x($wb,1)) 
	if(ok)
		if(V_start>V_theEnd)
			x0=V_Start
			V_start=V_theEnd
			V_theEnd=x0
			x0=V_startX
			V_startX=V_theEndX
			V_theEndX=x0
		endif
		if(V_start==V_theEnd)
			Abort t+"are too close together!"
		endif
	else
		Abort t+"must be on wave "+w+" or one like it"
	endif
End

Window DefinePeaks() : Table
	Edit /W=(25,90,367,257) W_MakeCentersX.y,W_MakeAmpsY.y,W_MakeWidthsX.y
	Modify width(Point)=64
EndMacro

Window PeakReportTable() : Table
	PauseUpdate; Silent 1		// building window...
	Edit /W=(5,42,401,195) W_PkCenters.y,W_PkAmps.y,W_PkFWHMs.y
	Append W_PkAreas.y
	Modify width(Point)=38
	Modify width(W_PkCenters.y)=84
	Modify width(W_PkAmps.y)=84
	Modify width(W_PkFWHMs.y)=84
	Modify width(W_PkAreas.y)=84
EndMacro

Window PeakReportLayout() : Layout
	PauseUpdate; Silent 1		// building window...
	Layout /W=(3,37,462,340) PeakReportTable(115,460,491,552)/O=2 as "Peak Report Layout"
	Textbox /N=title/A=LT/X=39.0909/Y=4.80769 "\JCYour Report Title Here"
	Modify left(title)=246,top(title)=67,width(title)=115,height(title)=14,frame(title)=4
	Textbox /N=info/A=LT/X=38.3636/Y=43.2692 "5 Gaussian peaks\r\rBaseline function: exp    "
	AppendText "\tK0 = 9.3137\r\tK1= 16469\r\tK2= 0.010146"
	Modify left(info)=242,top(info)=347,width(info)=122,height(info)=74,frame(info)=1
	Append PeakFitGraph(63,116,555,325)/O=1
EndMacro

Window PeakEstimates() : Table
	Edit /W=(35,87,576,268) W_EstCentersP.y,W_EstCentersX.y,W_EstAmpsY.y as "Peak Estimates"
	Append W_EstWidthsX.y,W_EstEdgesP.y
	Modify width=98,width(Point)=64
EndMacro

Window AreaReportTable() : Table
	Edit /W=(59,222,420,347) W_AreaX1.y,W_AreaX2.y,W_AreaNoBase.y as "Areas Between Cursors"
	Modify width(Point)=64
	Modify width(W_AreaNoBase.y)=112
EndMacro