//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////   
// data analysis software for SPELEEM - author: A.Locatelli 2004/04/07
// requires Igor 4.08 or later										
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////         

#pragma rtGlobals=1		// Use modern global access method.
#include  <All IP Procedures>
#include <Percentile and Box Plot>

Macro SpeLeem()
	// initialization of all global variables and strings
	variable/G scanNpnts, scanFrom, scanTo
	variable/G layer,energy,frameno
	variable/G lastp0=44, lastp1=473, lastq0=134, lastq1=563
	variable/G lastcursAp,lastcursAq, lastcursBp,lastcursBq
	variable/G xmcdupper, xmcdlower
	variable/G imageRange=0.3
	variable/G helRescaling=1
	variable/G lastWidth=5
	variable/G keg
	variable/G kefeg
	variable/G xzerog=231
	variable/G posfeg
	variable/G conversiong=28.63
	variable/G pa_thresh_low=128
	variable/G pa_thresh_high=128
	variable/G pa_nlevels=1
	String/G imagePathStr="C:data"
	String/G lastImagePathStr=""
	String/G imageTypeStr="png"
	String/G separatorStr="-"
	String/G scanNameStr=""
	String/G scanNumberStr=""
	String/G fileNumberStr=""
	String/G lastImage=""
	String/G fileNameStr=""
	String/G cmdStrg=""
	// creates speleem panel
	Dowindow/k mainw
	NewPanel/W=(0,10,600,80) as "Speleem image analysis"
	Dowindow/c mainw
	button browseImagePathButton, proc=browsePath ,pos={20,10} , size={100,20}, title="browse path" 
	button loadImageButton, proc=loadBWImage ,pos={130,10} ,  size={100,20}, title="load image" 
	button loadImageScanButton, proc=loadImageScan ,pos={240,10},  size={100,20}, title="GUI data"
	popupmenu imageTypePopup, pos={350,10}, size={100,20}, proc=setImageType, mode=1,value="PNG;tiff;tif;t16;jpeg;jpg;dat",title="image type: "
	button changeImagePathButton, proc=changePath ,pos={20,40} , size={100,20}, title="change path"
	button loadXMCDImageButton, proc=loadXMCDImage ,pos={130,40} ,  size={100,20}, title="load XMCD" 
	button loadSequenceButton, proc=loadSequence ,pos={240,40} ,  size={100,20}, title="load sequence" 
	button makeMovieButton, proc=sequenceToMovie ,pos={490,10} ,  size={100,20}, title="make movie" 
	popupmenu separatorPopup, pos={350,40}, size={100,20}, proc=setSeparatorType, mode=1,value="-;#",title="frame no. char: "
	button xmcdHrmButton, proc=loadXMCD_HRM ,pos={490,40} ,  size={100,20}, title="XMCD HRM" 
End

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// path functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////   

Function browsePath(ctrlname): buttoncontrol  //160203
	string ctrlname
	pathInfo imagePath
	SVAR imagePathStr=root:imagePathStr
	if (V_flag==1)   // path exists already
		newPath/M="please specify a path" imagePath imagePathStr	
	endif
	//enter path
	newPath/M="please specify a path"/O imagePath
	updatePath(1)
end

Function changePath(ctrlname): buttoncontrol //160203
	string ctrlname
	newPath/O imagePath getPath()
      	updatePath(1)
end

Function/S getPath()  //160203
	SVAR imagePathStr=root:imagePathStr
	String s=imagePathStr
	Prompt s, "path to experiment: "
	DoPrompt "please specify a path", s
	return s
End

Function updatePath(verbose) //160203	
	variable verbose
	SVAR imagePathStr=root:imagePathStr
	SVAR scanNameStr=root:scanNameStr
	SVAR scanNumberStr=root:scanNumberStr
	pathInfo imagePath
	if (V_flag==1) //ok
		imagePathStr=S_Path
		scanNameStr=returnFolderRootNameFromPath(imagePathStr)
		scanNumberStr=returnFolderNumberFromPath(imagePathStr)
		if (verbose==1)
		print "image path has been set to: ",imagePathStr,"scan folder: ", scanNameStr, "scan number: ", scanNumberStr
		endif
	else
		print "image path has not been set"
	endif
End

Function incrementPath()
	SVAR imagePathStr=root:imagePathStr
	SVAR scanNameStr=root:scanNameStr
	SVAR scanNumberStr=root:scanNumberStr
	pathInfo imagePath
	if (V_flag==1) //ok
		imagePathStr=S_Path
		variable len=strlen(imagePathstr)
		scanNumberStr=returnFolderNumberFromPath(imagePathStr)
		variable numlen=strlen(scannumberStr)
		String newNumberStr=num2ScanNumberStr(str2num(scanNumberStr)+1)
		String newPathStr = imagePathStr[0,len-numlen-2]+newNumberStr+":"
		newPath/O/Q/Z imagePath newPathStr
		if (V_flag==0)
			// general variables update
			imagePathStr=newPathStr
			scanNameStr=returnFolderRootNameFromPath(imagePathStr)
			scanNumberStr=returnFolderNumberFromPath(imagePathStr)
			return 0
		else
			return 1
		endif
	else 
		return 1
	endif
end

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// basic functionality
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Function/S OpenImage(filename,darkfield)
	String fileName
	variable darkfield
	SVAR imageTypeStr=root:imageTypeStr
	SVAR imagePathStr=root:imagePathStr
	SVAR scanNameStr=root:scanNameStr
	SVAR lastImage=root:lastImage
	String prevType = imageTypeStr
	// check if file exists
	Variable refNum
	if (stringMatch(filename,"")==1)
		fileName="image"
		Open/D/P=imagePath/R/M="Open image"/T="????" refNum as fileName
	else
		Open/P=imagePath/R/T="????"/Z=1 refNum as fileName
		imageTypeStr = returnImageType(S_fileName)
	endif
	// read image if exists	
	if (stringMatch(S_fileName,"")!=1)
		if (stringMatch(returnImageType(S_fileName),"tiff")==1||stringMatch(returnImageType(S_fileName),"tif")==1)
			ImageLoad/T=tiff/O/N=tmpimage/Q/P=imagePath S_fileName
			rename tmpimage0 tmpImage
		elseif (stringMatch(returnImageType(S_fileName),"t16")==1)
			ImageLoad/T=tiff/O/N=tmpimage/Q/P=imagePath S_fileName
			rename tmpimage0 tmpImage
		elseif ((stringMatch(returnImageType(S_fileName),"jpeg")==1)||(stringMatch(returnImageType(S_fileName),"JPG")==1))
			ImageLoad/T=JPG/O/N=tmpimage/Q/P=imagePath S_fileName	
		elseif ((stringMatch(returnImageType(S_fileName),"PNG")==1)||(stringMatch(returnImageType(S_fileName),"png")==1))
			ImageLoad/T=PNG/O/N=tmpimage/Q/P=imagePath S_fileName
		elseif ((stringMatch(returnImageType(S_fileName),"DAT")==1)||(stringMatch(returnImageType(S_fileName),"dat")==1))	
			LoadDatImage(S_fileName)
			S_path = S_Filename
		endif
		// close refNum
		// update globals
		imagePathStr=S_path
		updatePath(0)
		// give name
		lastImage = shortenWaveName(returnFileNameFromFullPath(S_filename))
		lastImage = removeForbiddenChar(lastImage)
		variable np=DimSize(tmpimage,0)
		variable nq=DimSize(tmpimage,1)
		make/n=(np,nq)/o $lastImage
		// 16 bit operation for png - RGB to gray for jpegs
		if ((stringMatch(returnImageType(S_fileName),"PNG")==1)||(stringMatch(returnImageType(S_fileName),"png")==1))
  	   		PNGimageTo2DWave($lastImage,tmpImage)
  	   	else
  	   		if  ((stringMatch(returnImageType(S_fileName),"jpeg")==1)||(stringMatch(returnImageType(S_fileName),"jpg")==1))
  	   			JPGimageTo2DWave($lastImage,tmpImage)
  	   		else
				duplicate/o tmpimage $lastimage
			endif
		endif
		redimension/U $lastImage
		// subtracting darkfiled
      		if (darkfield == 1)
      			darkFieldStack($lastimage,0)
      		endif
		killwaves/Z tmpimage
		print "opened",returnFileNameFromFullPath(S_filename),"image rows: ",np,"image columns: ",nq
		return lastimage
	else
		imageTypeStr = prevType
		return ""
	endif
end
	
Function loadBWImage(ctrlname): buttoncontrol	//160203
	string ctrlname
	SVAR imagePathStr=root:imagePathStr
	SVAR scanNameStr=root:scanNameStr
	SVAR lastImage=root:lastImage
	String lastStr
	variable darkfield=1 //2
       	prompt darkfield, "subtact darkfield", popup, "no;yes"
       	doPrompt "open image", darkfield
        if ((V_flag ==0 )&& (darkfield == 2))
		lastStr = OpenImage("",1)
	else
		lastStr = OpenImage("",0)
	endif
	if (stringMatch(lastStr,"")!=1)	//add controls
			NewImage/S=2 $lastStr
			PauseUpdate
			showinfo
			ControlBar 62
			ModifyGraph width=250*dimsize($lastimage,0)/dimsize($lastimage,1),height=250
			Execute "addBWcontrols()"
	endif	
End

Function loadImageScan(ctrlname): buttoncontrol	 //160203
	string ctrlname
	string tmpStr,newName
	variable j=0
	SVAR imagePathStr=root:imagePathStr
	SVAR imagetypeStr=root:imageTypeStr
	SVAR scanNameStr=root:scanNameStr
	SVAR scanNumberStr=root:scanNumberStr
	SVAR lastImage=root:lastImage
	SVAR separatorStr = root:separatorStr
	NVAR scanFrom=root:scanFrom
	NVAR scanTo=root:scanTo
	NVAR energy=root:energy
	NVAR layer=root:layer
	// load all waves in trace.txt
	// load info ImageNum STV     MOBJ Mesh BeamCurrent Box1 Box2 Box3 (kin energy scan)
	//             ImageNum Energy --------- Mesh BeamCurrent Box1 Box2 Box3 (photon energy scan)
	variable lastindex = strsearch(scanNameStr,"_",Inf,1)
	String nameStr=scanNameStr[0,lastindex-1]
	string scanType = scanNameStr[lastindex+1,strlen(scanNameStr)]
	print "scan type = ",scanType
	// find rootname
	print "file root name = ",nameStr
	string scaninfoStr=nameStr+"_trace.txt"
	print "load auxiliary waves from file: ", scaninfoStr
	Loadwave/A/D/G/O/Q/P=imagePath/W scaninfoStr
	lastindex = strsearch(nameStr,"_",Inf,1)
	scanNumberStr=nameStr[lastIndex+1,strlen(nameStr)]
	tmpStr = "ImageNum"+nameStr[lastIndex,strlen(nameStr)]
	print tmpStr,nameStr,scanNumberStr
	variable scanNpnts=numpnts($tmpStr)
	// rename waves
	if (exists(tmpStr)==1)
		newName = nameStr + "_NUM" 
		duplicate/o $tmpStr  $newName
		killwaves $tmpStr
	endif
	tmpStr = "STV_"+scanNumberStr
	if ( exists(tmpStr)==1 && (stringmatch(scanType,"XPS")==1) || (stringmatch(scanType,"IV")==1))
		//print "XPS scan type"
		newName =nameStr + "_STV" 
		duplicate/o $tmpStr  $newName
		scanFrom=returnWaveElement($newname,0)
		scanTo=returnWaveElement($newname,scanNpnts-1)
		killwaves $tmpStr
	endif
	tmpStr = "Energy_"+scanNumberStr
	if (exists(tmpStr)==1 && stringmatch(scanType,"XAS")==1)
		newName = NameStr + "_NRG" 
		duplicate/o $tmpStr  $newName
		scanFrom=returnWaveElement($newname,0)
		scanTo=returnWaveElement($newname,scanNpnts-1)
		killwaves $tmpStr
	endif
	tmpStr = "BeamCurrent_"+scanNumberStr
	if (exists(tmpStr)==1)
		newName = nameStr + "_BCR" 
		duplicate/o $tmpStr  $newName
		killwaves $tmpStr
	endif
	tmpStr = "Mesh_"+scanNumberStr
	if (exists(tmpStr)==1)
		newName = nameStr + "_MSH" 
		duplicate/o $tmpStr  $newName
		killwaves $tmpStr
	endif
	tmpStr = "MOBJ_"+scanNumberStr
	if (exists(tmpStr)==1 && stringmatch(scanType,"FOCUS")==1)
		newName = nameStr + "_OBJ" 
		duplicate/o $tmpStr  $newName
		killwaves $tmpStr
	endif
	// darkfield
    	variable darkfield=2
       	prompt darkfield, "subtact darkfield", popup, "no;yes"
       	doPrompt "open image", darkfield	
        if ((V_flag ==0 )&& (darkfield == 2))
		darkfield=1
	else
		darkfield=0
	endif
	//load first image in order to dimension the array
	print "number of images: ",scanNpnts,"start energy: ",scanFrom,"end energy: ",scanTo
	tmpStr = nameStr +"_"+num2ScanNumberStr(j)+"."+imageTypeStr 
	tmpStr = OpenImage(tmpStr,darkfield)
	if  (stringMatch(tmpStr,"")!=1)
    	// initialise the stack
   	variable np=DimSize($tmpStr,0)
   	variable nq=DimSize($tmpStr,1)      
      	make/n=(np,nq,scanNpnts)/O/W $nameStr
  	imageToStack($nameStr,$tmpStr,0)
    	j+=1
      // load the rest		
      do
       	tmpStr=nameStr +"_"+num2ScanNumberStr(j)+"."+imageTypeStr //scanNameStr +"_"+num2ScanNumberStr(j)+"."+imageTypeStr 
       	tmpStr = OpenImage(tmpStr,darkfield)
       	imageToStack($nameStr,$tmpStr,j)
       	killwaves $tmpStr
       	j+=1
       while (j<scanNpnts && (stringMatch(tmpStr,"")!=1) )
	// show
	if (stringMatch(tmpStr,"")!=1)
		NewImage/S=2 $nameStr
		ModifyImage $nameStr plane=0
		PauseUpdate
		showinfo
		ControlBar 81
		ModifyGraph width=250*np/nq,height=250
		//add controls
		String commandStr="addScanControls()"
		if ((stringMatch(num2Str(np),"NaN")!=1)&&(stringMatch(num2Str(nq),"NaN")!=1)) 
			Execute commandStr
			layer=0
       			energy=scanFrom
		endif
	endif
	endif
End

Function LoadSequence(ctrlname): buttoncontrol	 //160203
	string ctrlname
	// globals
	SVAR scanNameStr = root:scanNameStr
	SVAR imageTypeStr = root:imageTypeStr
	NVAR frameno,layer
	// prompt scan parameters
	String rootStr,fromStr,toStr,afterStr
	String fileNameStr=""
	variable step=1
	variable ref
	//variable lastindex = strlen(scanNameStr)-strsearch(scanNameStr,"_",3)
	variable lastindex = strsearch(scanNameStr,"_",Inf,1)
	string scanType = scanNameStr[lastindex+1,strlen(scanNameStr)]
	print "scan type = ",scanType
	rootStr=scanNameStr[0,lastindex-1]+"_"
	Prompt rootStr, "Enter Root-Name: "
	Prompt fromStr, "from: "
	Prompt toStr,"to:"
	Prompt afterStr,"enter string after number"
	Prompt step,"step:"
	DoPrompt "Load sequence", rootStr,fromStr,toStr,afterStr,step
	variable inputok = V_flag
	
	//  if user wants to continue
	if (inputok==0)
		// darkfield
    		variable darkfield=2
       		prompt darkfield, "subtact darkfield", popup, "no;yes"
       		doPrompt "open image", darkfield	
        		if ((V_flag ==0 )&& (darkfield == 2))
			darkfield=1
		else
			darkfield=0
		endif
		// avoids running low on VM	
		variable from=str2num(fromStr)
		variable to=str2num(toStr)
		variable memorylimit=2000
		if ((to-from)/step > memorylimit)
			to = from+memorylimit*step
		endif
		variable npoints = round((to-from)/step)+1
	 	variable len=strlen(fromStr)
	 	variable i=from
	 	variable index=0
	 	String sequenceNameStr=removeForbiddenChar(rootStr)
	 	String sequenceInfoStr=sequenceNameStr+"_INFO"
	 	make/n=(1,1,1)/O/U/W $sequenceNameStr
	 	make/n=(npoints)/o $sequenceInfoStr
	 	variable np=NaN,nq=NaN
		do
			// filename
			filenameStr=rootStr+returnSequenceNumber(i,len)+afterStr+"."+imageTypeStr 
       			filenameStr = OpenImage(filenameStr,darkfield)
       			// if exisits
       			If (stringMatch(filenameStr,"")!=1&&exists(filenameStr)==1)
       				// first loaded image sets the array dimension
       				if (index==0)
					np=DimSize($filenameStr,0)
					nq=DimSize($filenameStr,1)
					redimension/n=(np,nq,npoints)/U/W $sequenceNameStr
					Wavestats/Q $filenameStr
				endif
				redimension/U/W $sequenceNameStr
				imageToStack($sequenceNameStr,$filenameStr,index)
				entryToWave($sequenceInfoStr,i,index)
       				killwaves/Z $filenameStr
       				index+=1
       			endif
       			i+=step
       		while (i<=to )
       		// if successful, redimension and show
       		if (index > 0)
       			redimension/n=(np,nq,index) $sequenceNameStr
       			redimension/n=(index) $sequenceInfoStr
			NewImage/S=2 $sequenceNameStr
			ModifyImage $sequenceNameStr plane=0
			PauseUpdate
			showinfo
			ControlBar 81
			ModifyGraph width=250*np/nq,height=250
			// controls
			String commandStr="addSequenceControls()"
			if ((stringMatch(num2Str(np),"NaN")!=1)&&(stringMatch(num2Str(nq),"NaN")!=1)) 
				Execute commandStr
			endif
		endif
	endif
End

Function loadXMCDImage(ctrlname): buttoncontrol	
	string ctrlname
	string tmpStr
	variable j=0
	SVAR imagePathStr=root:imagePathStr
	SVAR imagetypeStr=root:imageTypeStr
	SVAR scanNameStr=root:scanNameStr
	SVAR scanNumberStr=root:scanNumberStr
	SVAR separatorStr=root:separatorStr
	SVAR lastImage=root:lastImage
	NVAR scanFrom=root:scanFrom
	NVAR scanTo=root:scanTo
	NVAR energy=root:energy
	NVAR layer=root:layer
	NVAR helRescaling=root:helRescaling
	// darkfield
    	variable darkfield=2
       	prompt darkfield, "subtact darkfield", popup, "no;yes"
       	doPrompt "open image", darkfield	
        if ((V_flag ==0 )&& (darkfield == 2))
		darkfield=1
	else
		darkfield=0
	endif
	// open beamline info text file in order to understand whether phase is negative or positive
	tmpStr = scanNameStr +"_beam"+separatorStr+num2ScanNumberStr(j)+"."+"txt"
	variable phaseSign = getPhaseSign(tmpStr)  // is 000 negative? returns sign
	print tmpStr, "has helicity: ", phasesign
	// open image 000
	tmpStr = scanNameStr+separatorStr+num2ScanNumberStr(j)+"."+imageTypeStr
	print "load file: "+tmpstr 
	tmpStr = OpenImage(tmpStr,darkfield)	
//removeBackGround($tmpStr) //background removal on raw data
	// initialise the stack  
	variable np=DimSize($tmpStr,0)
	variable nq=DimSize($tmpStr,1)     
      make/n=(np,nq,2)/o/d $scanNameStr
      j+=1
      String tmpStr1 = scanNameStr+separatorStr+num2ScanNumberStr(j)+"."+imageTypeStr
      print "load file: "+tmpstr1 
      tmpStr1 = OpenImage(tmpStr1,darkfield)
 //removeBackGround($tmpStr1) //background removal on raw data
	if (phaseSign == -1) // 000 is negative phase
       	imageToStack($scanNameStr,$tmpStr,0)
       	imageToStack($scanNameStr,$tmpStr1,1)
    	else
       	imageToStack($scanNameStr,$tmpStr,1)
       	imageToStack($scanNameStr,$tmpStr1,0)
       endif
       // rescaling for plus helcity
       variable coef = helRescaling
       prompt coef, "rescaling factor for helicity minus"
       doPrompt "XMCD input", coef
       if (V_flag==0)
       		MultiplyStackLayer($scanNameStr,coef,0)
       		helRescaling = coef
       endif
       killwaves/Z $tmpStr, $tmpStr1
       // sho raw data
	if (stringMatch(tmpStr,"")!=1 || stringMatch(tmpStr1,"")!=1)
		NewImage/S=2 $scanNameStr
		ModifyImage $scanNameStr plane=0
		layer=0
		if ((stringMatch(num2Str(np),"NaN")!=1)&&(stringMatch(num2Str(nq),"NaN")!=1))	
			execute "addRawDataControls()"
			ModifyGraph width=250*np/nq,height=250
		endif
		Execute "AutoSizeImage(1,1)" 
	endif
	// xmcd subtraction
	if (stringMatch(tmpStr,"")!=1 && stringMatch(tmpStr1,"")!=1 && (phaseSign != 0) )
		String xmcdStr=scanNameStr +"_"+"XMCD"
		make/n=(np,nq,1)/O/D $xmcdStr
		calcXMCD($scanNameStr,$xmcdStr,-1)
		// display
 		NewImage/S=2 $xmcdStr
		if ((stringMatch(num2Str(np),"NaN")!=1)&&(stringMatch(num2Str(nq),"NaN")!=1)) 
			Execute "addXMCDControls()"
			ModifyGraph width=250*np/nq,height=250
			TextBox/C/N=text0/A=LT xmcdStr
			DoUpdate
			ModifyGraph width=0,height=0
			ModifyImage $xmcdstr ctab= {-0.3,0.3,Grays,0}
		endif
	endif
	incrementPath()
end

Function loadXMCD_HRM(ctrlname): buttoncontrol	
	string ctrlname
	string tmpStr
	variable j=0
	SVAR imagePathStr=root:imagePathStr
	SVAR imagetypeStr=root:imageTypeStr
	SVAR scanNameStr=root:scanNameStr
	SVAR scanNumberStr=root:scanNumberStr
	SVAR separatorStr=root:separatorStr
	SVAR lastImage=root:lastImage
	NVAR scanFrom=root:scanFrom
	NVAR scanTo=root:scanTo
	NVAR energy=root:energy
	NVAR layer=root:layer
	NVAR helRescaling=root:helRescaling
	// darkfield
    	variable darkfield=2
       	prompt darkfield, "subtact darkfield", popup, "no;yes"
       	doPrompt "open image", darkfield	
        if ((V_flag ==0 )&&(darkfield == 2))
		darkfield=1
	else
		darkfield=0
	endif
	String filenameStr=""
	String beamlineStr=""
	variable refnum
	variable ok=0
	do
		beamlineStr = scanNameStr +"_beam"+separatorStr+num2ScanNumberStr(j)+".txt"
		open/R/P=imagePath/t="TEXT"/Z=1 refnum as beamlineStr
		if (V_Flag==0)
			ok+=1
			variable phaseSign = getPhaseSign(beamlineStr) 
			close refnum
			if (phaseSign>0) 
				print "phase is positive"
			else
				print "phase is negative"
			endif
			filenameStr=scanNameStr +separatorStr+num2ScanNumberStr(j)+"."+imageTypeStr 	
       			filenameStr = OpenImage(filenameStr,darkfield)
       			if (j==0)
       				// initialise the stack  
				variable np=DimSize($filenameStr,0)
				variable nq=DimSize($filenameStr,1)     
       				make/n=(np,nq,2)/o/d $scanNameStr
       				make/n=(np,nq)/o/d plus,minus
       			endif
       		else
       			filenameStr=scanNameStr +separatorStr+num2ScanNumberStr(j)+"."+imageTypeStr 	
       			filenameStr = OpenImage(filenameStr,darkfield)
		endif
		if (stringMatch(filenameStr,"")!=1&&exists(filenameStr)==1)
			if (phaseSign>0)
       				plus=addImages($filenameStr,plus)
       			else
       				minus=addImages($filenameStr,minus)
       			endif
			killwaves/Z $filenameStr
		endif
		j+=1
	while (strlen(filenameStr)>0)
	imageToStack($scanNameStr,plus,1)
       	imageToStack($scanNameStr,minus,0)
       	 killwaves/Z plus, minus
       	// rescaling for plus helcity
   	variable coef = helRescaling
	prompt coef, "rescaling factor for helicity minus"
	doPrompt "XMCD input", coef
	if (V_flag==0)
       		MultiplyStackLayer($scanNameStr,coef,0)
       		helRescaling = coef
	endif
       	// sho raw data
	if (ok > 1)
		NewImage/S=2 $scanNameStr
		ModifyImage $scanNameStr plane=0
		layer=0
		if ((stringMatch(num2Str(np),"NaN")!=1)&&(stringMatch(num2Str(nq),"NaN")!=1))	
			execute "addRawDataControls()"
			ModifyGraph width=250*np/nq,height=250
		endif
		Execute "AutoSizeImage(1,1)" 
		String xmcdStr=scanNameStr +"_"+"XMCD"
		make/n=(np,nq,1)/O/D $xmcdStr
		calcXMCD($scanNameStr,$xmcdStr,-1)
		// display
 		NewImage/S=2 $xmcdStr
		if ((stringMatch(num2Str(np),"NaN")!=1)&&(stringMatch(num2Str(nq),"NaN")!=1)) 
			Execute "addXMCDControls()"
			ModifyGraph width=250*np/nq,height=250
			TextBox/C/N=text0/A=LT xmcdStr
			DoUpdate
			ModifyGraph width=0,height=0
			ModifyImage $xmcdstr ctab= {-0.3,0.3,Grays,0}
		endif
	endif
end
	
Function sequenceToMovie(ctrlname): buttoncontrol	 //160203
	string ctrlname
	SVAR scanNameStr = root:scanNameStr
	SVAR imageTypeStr = root:imageTypeStr
	NVAR frameno,layer
  	NVAR lastp0=root:lastp0
	NVAR lastp1=root:lastp1
	NVAR lastq0=root:lastq0
	NVAR lastq1=root:lastq1
	String rootStr,fromStr,toStr
	String fileNameStr=""
	variable step=1
	variable ref
	variable V_init=0
	variable redimensionFlag=1
	rootStr=scanNameStr
	variable deltaT=10
	Prompt rootStr, "Enter Root-Name: "
	Prompt fromStr, "from: "
	Prompt toStr,"to:"
	Prompt step,"step:"
	Prompt redimensionFlag,  "redimension image",popup "no;yes"
	Prompt deltaT, "increment (s)"
	DoPrompt "Load sequence", rootStr,fromStr,toStr,step,redimensionFlag,deltaT
	variable t=0
	variable inputok = V_flag
	if (inputok==0)
		// darkfield
    		variable darkfield=2
       		prompt darkfield, "subtact darkfield", popup, "no;yes"
       		doPrompt "open image", darkfield	
        		if ((V_flag ==0 )&& (darkfield == 2))
			darkfield=1
		else
			darkfield=0
		endif
		variable from=str2num(fromStr)
		variable to=str2num(toStr)
	 	variable len=strlen(fromStr)
	 	variable i=from
	 	variable index=0
	 	variable np=-1
	 	variable nq=-1
		Variable nframes = 2
		Prompt nframes, "frames per seconds"
		DoPrompt "Set speed (frames per second)", nframes
		String textboxStr=""
		String movieStr = rootStr
		make/n=(1,1,1)/o $rootStr=0
		do
			filenameStr=rootStr+returnSequenceNumber(i,len)+"."+imageTypeStr 
       			filenameStr = OpenImage(filenameStr,darkfield)
       			If (stringMatch(filenameStr,"")!=1&&exists(filenameStr)==1)
       				np=DimSize($filenameStr,0)
				nq=DimSize($filenameStr,1)
       				make/n=(np,nq,1)/O/W $rootStr
  				imageToStack($rootStr,$filenameStr,0)
  				if (redimensionFlag == 2)
  					redimensionImage(rootStr,lastp0,lastp1, lastq0, lastq1)
  				endif		
				if (V_init == 0)   			
					// create window
					String winStr = "movie_sequence_utility"
					Dowindow $winStr
					if (V_flag==0)
						Display
						DoWindow/C $winStr
						appendImage/L/T $rootStr
						ModifyGraph nticks=0
						ModifyGraph standoff=0
						// ModifyGraph width=0,height=0
						ModifyGraph margin=5
						SetAxis/A/R left
						TextBox/C/N=text0/A=LT ""
						execute "autosizeimage(1,1)"
						doUpdate
					else
						DoWindow/F $winStr
					endif
					V_init = 1
					NewMovie /F=(nframes)/L/O/Z as movieStr
				endif
				t=(i)*deltaT
       				textboxStr ="frame # "+ num2Str(i)+" time: "+ num2Str(t)+" s"
				TextBox/C/N=text0 textboxStr 
				doUpdate
       				AddMovieFrame
       				killwaves/Z $filenameStr
       				index+=1
       			endif
       			i+=step
       		while ((i<=to))
       		closeMovie
	endif
End

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// button functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Function normalise(ctrlname): buttoncontrol //160203
	String ctrlname
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	String normWaveStr
	Prompt normWaveStr, "enter norm wave",popup,WaveList("n*",";","")
	DoPrompt "Choose normalisation wave", normWaveStr 
	variable doNormalise=V_Flag
	variable npNormalise = -1
	variable nqNormalise = -1
	variable normRms=1
	variable np=DimSize($tmpStr,0)
	variable nq=DimSize($tmpStr,1)
	if (doNormalise==0)
		npNormalise = dimSize($normWaveStr,0)
		nqNormalise = dimSize($normWaveStr,1)
		redimension/D $normWaveStr
		Wavestats/Q $normWaveStr
		normRms=V_rms
		variable coef=2
		if (np==npNormalise&&nq==nqNormalise)
			backGroundCorrection($tmpStr,$normWaveStr,coef,normRms)
		endif
	endif
end

Function normaliseSequence(ctrlname): buttoncontrol //160203
	String ctrlname
	NVAR frameno = root:frameno
	NVAR layer = root:layer
	NVAR energy = root:energy
	String sourceNameStr=returnListFirstElement(imageNameList("",""),";")
	sourceNameStr=removeQuotes(sourceNameStr)
	variable pp=dimsize($sourceNameStr,0)
	variable qq=dimsize($sourceNameStr,1)	
	variable n=dimsize($sourceNameStr,2)
	make/n=(pp,qq)/o tmp
	// norm wave
	String normWaveStr
	Prompt normWaveStr, "enter norm wave",popup,WaveList("n*",";","")
	DoPrompt "Choose normalisation wave", normWaveStr 
	Wavestats/Q $normWaveStr
	variable normRms=V_rms
	variable coef=V_rms
	variable npp=dimsize($normWaveStr,0)
	variable nqq=dimsize($normWaveStr,1)
	if ((npp==pp)&&(nqq==qq))
	variable i=0
		do
		setGlobals(sourceNameStr,i)
		layer=i
		extract2DWaveFromStack($sourceNameStr,i,tmp)		
		redimension/D tmp
		backGroundCorrection(tmp,$normWaveStr,coef,normRms)
		redimension/U tmp
		//removeBackGround(tmp)
		insert2DWaveinStack($sourceNameStr,i,tmp)
		
		i+=1
	while(i<n)
	else
		Print "***** normalisation wave dimension does not match target! ******************"
	endif
	killWaves/z tmp
end

Function shoLayer(ctrlnamesel,varNum,varStr,varName): setVariableControl   //160203
	String ctrlnamesel
	Variable varNum
	String varStr
	String varName
	// globals
	NVAR frameno = root:frameno
	NVAR layer = root:layer
	NVAR energy = root:energy
	String sourceNameStr=returnListFirstElement(imageNameList("",""),";")
	sourceNameStr=removeQuotes(sourceNameStr)
	// auto scaling image range
	variable pp=dimsize($sourceNameStr,0)
	variable qq=dimsize($sourceNameStr,1)
	make/n=(pp,qq)/o tmp
	extract2DWaveFromStack($sourceNameStr,layer,tmp)
	waveStats/Q tmp
	//ModifyImage $sourceNameStr ctab= {V_min,V_max,,0}; delayUpdate;
	killWaves/z tmp
	// update globals
	variable n=dimsize($sourceNameStr,2)
	setGlobals(sourceNameStr,varNum)
	if ((varNum<0)||(varNum>=n)) 
		if (varNum<0)
			varNum=0
			layer=0	
		endif
		if (varNum>=n)
			varNum=n-1
			layer=n-1
		endif
	else
		modifyImage $sourcenameStr plane=varNum
	endif
End

Function setGlobals(nameStr,x) // sets appropriate value of layer,energy,frameNo
	String nameStr
	variable x
	NVAR frameno = root:frameno
	NVAR layer = root:layer
	NVAR energy = root:energy
	SVAR fileNameStr = root:fileNameStr
	String testStr=nameStr+"_INFO"
	String xwaveStr=""
	if (exists(testStr)==1)
		xwaveStr=testStr
		frameno=returnWaveElement($xwaveStr,x)
		fileNameStr = nameStr+num2ScanNumberStr(frameno)
	endif
	testStr=nameStr+"_NRG"
	if (exists(testStr)==1)
		xwaveStr=testStr
		energy=returnWaveElement($xwaveStr,x)
		frameno=x
		fileNameStr = nameStr+"_"+num2ScanNumberStr(x)
	endif
	testStr=nameStr+"_OBJ"
	if (exists(testStr)==1)
		wave w = $testStr
		xwaveStr=testStr
		energy=returnWaveElement($xwaveStr,x)
		frameno=x
		fileNameStr = nameStr+"_"+num2ScanNumberStr(x)
	endif
	testStr=nameStr+"_STV"	
	if (exists(testStr)==1)
		wave w = $testStr
		if (w[0]!=w[numpnts(w)-1])
			xwaveStr=testStr
			energy=returnWaveElement($xwaveStr,x)
    			frameno=x
			fileNameStr = nameStr+"_"+num2ScanNumberStr(x)
		endif
	endif
end

Function setImageType(ctrlName,popNum,popStr) : PopupMenuControl //160203
	String ctrlName
	Variable popNum
	String popStr
	SVAR imageTypeStr = root:imageTypeStr
      imageTypeStr=popStr
      print "image type has been set to: ",imageTypeStr
End

Function setSeparatorType(ctrlName,popNum,popStr) : PopupMenuControl //160203
	String ctrlName
	Variable popNum
	String popStr
	SVAR separatorStr = root:separatorStr
      separatorStr=popStr
      print "frame no. character is set to: ",separatorStr
End

Function shoMovie(ctrlnamesho): buttoncontrol //160203
	String ctrlnamesho
	NVAR frameno = root:frameno
	NVAR layer = root:layer
	NVAR energy = root:energy
	String sourceNameStr=returnListFirstElement(imageNameList("",""),";")
	sourceNameStr=removeQuotes(sourceNameStr)
	variable pp=dimsize($sourceNameStr,0)
	variable qq=dimsize($sourceNameStr,1)	
	variable n=dimsize($sourceNameStr,2)
	make/n=(pp,qq)/o tmp
	String testStr,xwaveStr
	variable i=0
	
	layer=i
		extract2DWaveFromStack($sourceNameStr,i,tmp)
		waveStats/Q tmp
		//ModifyImage $sourceNameStr ctab= {V_min,V_max,,0}; delayUpdate;
		
	do
		setGlobals(sourceNameStr,i)
		layer=i
		extract2DWaveFromStack($sourceNameStr,i,tmp)
		waveStats/Q tmp
		//ModifyImage $sourceNameStr ctab= {V_min,V_max,,0}; delayUpdate;
		modifyImage $sourceNameStr plane=i
		doUpdate
		i+=1
	while(i<n)
	killWaves/z tmp
End

Function makeMovie(ctrlnamesho): buttoncontrol //160203
	String ctrlnamesho
	NVAR frameno = root:frameno
	NVAR layer = root:layer
	NVAR energy = root:energy
	String sourceNameStr=returnListFirstElement(imageNameList("",""),";")
	sourceNameStr=removeQuotes(sourceNameStr)
	ModifyGraph nticks=0,minor=0
	variable n=dimsize($sourceNameStr,2)
	String testStr,xwaveStr,textboxStr,tmpStr,twaveStr,timeWaveStr
	variable i=0
	variable nframes=2
	Prompt nframes, "movie speed: "
	DoPrompt "Set speed (frames per second)", nframes
	if (V_flag==0)
		//TextBox/C/N=text0/A=LB""
		//TextBox/C/N=text1/A=RB ""
		NewMovie /F=(nframes)/L/O/Z as sourceNameStr
		do
			layer=i
			testStr = sourceNameStr+"_INFO"
			twaveStr = sourceNameStr+"_TEMP"
			timewaveStr = sourceNameStr+"_TIME"
			xwaveStr=""
			if (exists(twaveStr)==1)
				xwaveStr=testStr
				frameno=returnWaveElement($testStr,i)
				textboxStr ="image # "+ num2str(returnWaveElement($xwaveStr,i)) 
				//tmpStr="T = "+num2str(returnWaveElement($twaveStr,frameno)) +" C"
				sprintf tmpStr "%*.*f\r" 4, 2, returnWaveElement($twaveStr,frameno)
				tmpStr = "T= "+ tmpStr + "C"
				TextBox/C/N=text0 textboxStr
				TextBox/C/N=text1 tmpStr
			endif
			if (exists(timeWaveStr)==1)
				xwaveStr=testStr
				frameno=returnWaveElement($xwaveStr,i)
				//textboxStr ="frame # "+ num2str(returnWaveElement($xwaveStr,i)) 
				sprintf tmpStr "%*.*f" 4, 1, returnWaveElement($timeWaveStr,i)
				tmpStr=tmpstr +" s"
				//TextBox/C/N=text0 textboxStr
				TextBox/C/N=text1 tmpStr
			endif
			testStr=sourceNameStr+"_STV"	
			if (exists(testStr)==1)
				xwaveStr=testStr
				energy=returnWaveElement($xwaveStr,i)
				textboxStr ="frame # "+ num2str(i) 
				tmpStr="STV = "+num2str(energy) +" eV"
				TextBox/C/N=text0 textboxStr
				TextBox/C/N=text1 tmpStr
			endif
			testStr=sourceNameStr+"_NRG"
			if (exists(testStr)==1)
				xwaveStr=testStr
				energy=returnWaveElement($xwaveStr,i)
				textboxStr ="frame # "+ num2str(i) 
				tmpStr="hv = "+num2str(energy) +" eV"
				TextBox/C/N=text0 textboxStr
				TextBox/C/N=text1 tmpStr
			endif
			modifyImage $sourceNameStr plane=i
			doUpdate
			AddMovieFrame
			i+=1
		while(i<n)
		TextBox/K/N=text0
		TextBox/K/N=text1
		closeMovie
	endif
End

Function shoPixelScan(ctrlnamepix): buttoncontrol //160203
	String ctrlnamepix
	variable i
	Silent 1; PauseUpdate
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	variable n=dimsize($tmpStr,2)
	String pixelStr=tmpStr
	if (stringMatch(num2Str(pcsr(A)),"NaN")!=1)
		pixelStr=pixelStr+"_p"+num2Str(pcsr(A))+"_q"+num2Str(qcsr(A))
		make/N=(n)/O $pixelStr
		extract1DWaveFromStack($tmpStr,pcsr(A),qcsr(A),$pixelStr)
		//normalisation
		normaliseScan(tmpStr,pixelStr)
		String graphNameStr="SCAN_"+tmpStr
		stackToFrontAndDisplay(graphNameStr,tmpStr,pixelStr)
		Execute "addPixelScanControls()"
	endif
End

Function normaliseXAS(): graphmarquee 
	//String tmpStr=returnListFirstElement(WaveList("*ROI", ";","WIN:"),";")
	//
	//
       String tmpStr
	Prompt tmpStr, "enter wave",popup WaveList("*",";","WIN:")
	doPrompt "choose wave" tmpStr
	if (V_flag == 0)
		tmpStr= removeQuotes(tmpStr)
		String meshStr=tmpStr[0,strlen(tmpStr)-4]+"MSH"
		String darkStr=tmpStr[0,strlen(tmpStr)-4]+"DARK"
		String targetStr = tmpStr+"N"
		String nrgStr = tmpStr[0,strlen(tmpStr)-4]+"NRG"
		duplicate/o $tmpStr $targetStr
		wave w = $targetStr
		wave d = $darkStr
		wave m = $meshStr
		w-=d
		w/= (m-0.0127)
		display $targetStr vs $nrgStr
	endif
       print tmpStr
end

Function shoRoiScan(ctrlnameroi): buttoncontrol	//160203
	String ctrlnameroi
	Silent 1; PauseUpdate
	variable i=0
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	variable p=dimsize($tmpStr,0)
	variable q=dimsize($tmpStr,1)
	variable r=dimsize($tmpStr,2)
	if (r==0)
		r=1
		redimension/N=(p,q,r) $tmpStr 
	endif
	String roiStr=tmpStr+"_ROI"
	roiStr = giveRoiName(roiStr)
	make/N=(r)/O/U $roiStr
	make/n=(p,q)/o myLayer
	Silent, 1
	PauseUpdate
	do
		extract2DWaveFromStack($tmpStr,i,myLayer)
		redimension/D mylayer
		ImageStats/R=M_RoiMask myLayer
		entryToWave($roiStr,V_avg,i)
		i+=1
	while(i<=r)
	killwaves myLayer
	//normalisation
	//normaliseScan(tmpStr,roiStr)
	String graphNameStr="SCAN_"+tmpStr
	stackToFrontAndDisplay(graphNameStr,tmpStr,roiStr)
	Execute "addROIScanControls()"
End

Function normaliseScan(rootnameStr,targetStr)
	String rootnameStr,targetStr
	String testStr=rootnameStr+"_STV"
	String normStr=""
	if (exists(testStr)==1)
		normStr=rootnameStr+"_BCR"
		normaliseWave($targetStr,$normStr)
	endif
	testStr=rootnameStr+"_NRG"	
	if (exists(testStr)==1)
		normStr=rootnameStr+"_MSH"
		normaliseWave($targetStr,$normStr)
	endif
end

Function/S giveRoiName(defaultName)
	String defaultName
	String previousName=defaultName
	Prompt defaultName, "Please give ROI scan a name: "		
	DoPrompt "ROI Scan name", defaultName
	if (V_flag==0)
		return defaultName
	else
		return previousName
	endif
End

Function shoRoi(ctrlname): buttoncontrol //160203	
	string ctrlname
	WMCreateImageROIPanel();
End

Function InvertImage(ctrlname): buttoncontrol
	string ctrlname	
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	wave w=$tmpStr 
	Wavestats/Q $tmpStr  // [a,b]
	redimension/D $tmpStr
	w=-1*w // [-b,-a]
	w+=(V_max+V_min)  // [a,b]
End

Function saveJpeg(ctrlname): buttoncontrol	//160203
	string ctrlname
	NVAR layer=root:layer
	NVAR frameno=root:frameno
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr=removeQuotes(tmpStr)
	if (dimsize($tmpStr,2)==0)
		layer = 0
	endif
	make/n=(dimSize($tmpStr,0),dimSize($tmpStr,1))/o/u thislayer    //=$tmpStr
	extract2DWaveFromStack($tmpStr,layer,thislayer)
	if (dimsize($tmpStr,2)>1)
		tmpStr=tmpStr+"_"+num2ScanNumberStr(frameno)
	endif
	PathInfo savePath
	if ( V_Flag==0) // if path does not exist ask 
		newPath/M="please specify a path for Saving data: " savePath
	endif
	variable refnum
	Open/D/P=savePath  refnum as tmpStr+".jpg"  // get path
	if (strlen(S_filename)>0)
		String newPathStr=returnPathFromFullPath(S_filename)
		newPath/O/Q SavePath newPathStr //update path
		ImageSave/O/P=savePath/T="JPEG" thislayer as S_filename //save
		print "data saved in: "+newPathStr+" as "+S_filename
	else
		print "save operation cancelled"
	endif
	killwaves thislayer
	close/A
End

Function saveTiff(ctrlname): buttoncontrol	//160203
	string ctrlname
	NVAR layer=root:layer
	NVAR frameno=root:frameno
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr=removeQuotes(tmpStr)
	wave w=$tmpStr
	if (dimsize($tmpStr,2)==0)
		layer = 0
	endif
	make/n=(dimSize($tmpStr,0),dimSize($tmpStr,1))/o thislayer
	extract2DWaveFromStack($tmpStr,layer,thisLayer)
	if (dimsize($tmpStr,2)>1)
		tmpStr=tmpStr+"_"+num2ScanNumberStr(frameno)
	endif
	PathInfo savePath
	if ( V_Flag==0) // if path does not exist 
		newPath/M="please specify a path for Saving data: " savePath
	endif
	variable refnum
	Open/D/P=savePath  refnum as tmpStr+".tif"  // get path
	if (strlen(S_filename)>0)
		redimension/U/W thisLayer
		String newPathStr=returnPathFromFullPath(S_filename)
		print S_filename
		newPath/O/Q SavePath newPathStr //update path
		ImageSave/D=16/O/T="TIFF"/P=savePath  thislayer as S_filename  //save
		print "data saved in: "+newPathStr+" as "+S_filename
	else
		print "save operation cancelled"
	endif
	close/a
	killwaves thislayer
End

Function scanToJpeg(ctrlnamesho): buttoncontrol //160203
	String ctrlnamesho
	NVAR frameno = root:frameno
	NVAR layer = root:layer
	NVAR energy = root:energy
	SVAR filenameStr = root:filenameStr
	String sourceNameStr=returnListFirstElement(imageNameList("",""),";")
	sourceNameStr=removeQuotes(sourceNameStr)
	variable n=dimsize($sourceNameStr,2)
	String testStr,xwaveStr
	wave thislayer
	make/n=(dimsize($sourceNameStr,0),dimsize($sourceNameStr,1)) thislayer
	PathInfo savePath
	if ( V_Flag==0)
		newPath/M="please specify a path for Saving data: " savePath
	endif
	variable i=0
	String temp,tempi
	do
		setGlobals(sourceNameStr,i)
		layer=i
		modifyImage $sourceNameStr plane=i
		doUpdate
		extract2DWaveFromStack($sourceNameStr,i,thislayer)
		PathInfo savePath
		if (i<10)
			temp="00"+num2str(i)
		else
			if (i<100)
				temp="0"+num2str(i)
			else
				temp=num2str(i)
			endif
		endif	
		tempi=sourceNameStr+"_"+temp+".jpg"
		print "saving "+sourceNameStr+" as "+tempi+" in "+S_path
		ImageSave/O/P=savePath/T="JPEG" thislayer as tempi
		//print "saving "+fileNameStr+" as "+fileNameStr+".tif"+" in "+S_path
		//ImageSave/D=16/O/T="TIFF"/P=savePath thislayer as fileNameStr+".tif"
		i+=1
	while(i<n)
	killwaves thislayer
End

Function PQaverage(ctrlname): buttoncontrol  //160203	
	string ctrlname
	Silent 1; pauseUpdate
	String sep = ";"
	String listStr=WaveList("*_p*", sep, "WIN:")
	String tmpStr=returnListFirstElement(listStr,sep)
	tmpStr=removeQuotes(tmpStr)
	variable i=strsearch(tmpStr, "_p",0)
	String sourceNameStr=tmpStr[0,i-1]
	String graphStr="SCAN_"+sourceNameStr
	if (bringToFront(graphStr)==1) // the graph exists
		if (strlen(listStr)>0)         //at least one pq wave exists
      			String avgStr=sourceNameStr+"_avg"
       		make/n=(dimsize($sourceNameStr,2))/o $avgStr=0
       		wave wavg=$avgStr
       		wave wtmp=$tmpStr
       		i=0
       		do
       			tmpStr=returnListFirstElement(listStr,sep)
       			wavg+=wtmp
       			listStr=listStr[strlen(tmpStr)+1,strlen(listStr)]
       			i+=1       			
       		while (strlen(listStr)>0)
       		wavg/=i
       		stackToFrontAndDisplay(graphStr,sourceNameStr,avgStr)
       	endif
       endif
 end
 
Function ImageSubset() : GraphMarquee //160203
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	NVAR lastp0=root:lastp0
	NVAR lastp1=root:lastp1
	NVAR lastq0=root:lastq0
	NVAR lastq1=root:lastq1
	GetMarquee/Z/K left, top
	if (V_Flag == 0)
		Print "There is no marquee"
	else
		lastp0=round(V_top)
		lastp1=round(V_bottom)
		lastq0=round(V_left)
		lastq1=round(V_right)
		redimensionImage(tmpStr, lastp0, lastp1, lastq0, lastq1)
	endif
	Execute "AutoSizeImage(1,1)"
End

Function previousSubset() : GraphMarquee  //160203
	NVAR lastp0=root:lastp0
	NVAR lastp1=root:lastp1
	NVAR lastq0=root:lastq0
	NVAR lastq1=root:lastq1
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	redimensionImage(tmpStr,lastp0,lastp1, lastq0, lastq1)
	Execute "AutoSizeImage(1,1)"
End

Function prevSubset(ctrlname): buttoncontrol   //160203
	string ctrlname
	NVAR lastp0=root:lastp0
	NVAR lastp1=root:lastp1
	NVAR lastq0=root:lastq0
	NVAR lastq1=root:lastq1
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	redimensionImage(tmpStr,lastp0,lastp1, lastq0, lastq1)
	Execute "AutoSizeImage(1,1)"
	//ModifyGraph nticks=0; ModifyImage $tmpstr ctab= {*,*,YellowHot,0}
End

 Function redimensionImage(imageStr,row1,row2,col1,col2)  //160203
 	String imageStr
 	variable row1,row2,col1,col2
 	variable nRows,nCols, nLayers
 	wave imagewave=$imageStr
 	nCols=dimSize(imagewave,0)
 	nRows=dimSize(imagewave,1)
 	nLayers=dimSize(imagewave,2)
 	if (row1>nRows)
		row1=nRows
	endif
	if (row2>nRows)
		row2=nRows
	endif
	if (col1>nCols)
		col1=nCols
	endif
	if (row2>nCols)
		col2=nCols
	endif
	variable c1,c2,r1,r2
	c1=min(col1,col2)
	c2=max(col1,col2)
	r1=min(row1,row2)
	r2=max(row1,row2)
	DeletePoints/M=0   (c2+1),(nCols-1), imagewave
 	DeletePoints/M=0   0,(c1-1), imagewave
 	DeletePoints/M=1   (r2+1),(nRows-1), imagewave
 	DeletePoints/M=1   0,(r1-1), imagewave
 	if (nLayers <= 1)
 		redimension/N=(c2-c1,r2-r1) imagewave
 	else
 		redimension/N=(c2-c1,r2-r1,nLayers) imagewave
 	endif
 	
 	redimension/N=(c2-c1,r2-r1,nLayers) imagewave
 end
 
Function despikeScan(ctrlname): buttoncontrol   //160203
	string ctrlname
	Silent 1; pauseUpdate
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	despike(tmpStr) 	
end

Function despike(imageStr)   //160203
 	String imageStr
 	wave w=$imageStr
 	WaveStats/Q w
 	Variable lowerThreshold=V_min,upperThreshold=V_max
	Prompt lowerThreshold, "Enter Lower Threshold: "
	Prompt upperThreshold, "Enter Upper Threshold: "
	DoPrompt "Set threshold", lowerThreshold, upperThreshold
	if (V_flag==0)
 		w[0,2][0,2]=upperThreshold
 		w[0,2][2,4]=lowerThreshold
 		w=w*(w<upperThreshold)+upperThreshold*(w>=upperThreshold)
 		w=w*(w>lowerThreshold)+lowerThreshold*(w<=lowerThreshold)
 		String tmpStr=returnListFirstElement(imageNameList("",""),";")
		tmpStr= removeQuotes(tmpStr)
		ModifyImage $tmpStr ctab= {lowerThreshold,upperThreshold,,0}
 	endif
 end
 
Function rescaleScan(ctrlname): buttoncontrol  //160203	
	string ctrlname
	Silent 1; pauseUpdate
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	rescale(tmpStr)
end

Function/S rescale(imageStr)  //160203
 	String imageStr	
 	wave w=$imageStr
 	WaveStats/Q w
 	Variable newMin=0,newMax=255
	Prompt newMin, "Enter new min: "
	Prompt newMax, "Enter new max: "		
	DoPrompt "Set Range", newMin, newMax
	if (V_flag==0)
 		redimension/D w
 		w=(w-V_min)/(V_max-V_min)*(newMax-newMin)+newMin
 	endif
 	//n=dimsize(w,2)
 	//w[0,10][0,10][]=newMin
 	//w[0,10][5,10][]=newMax
 end

Function shoProfile(ctrlname): buttoncontrol	  //160203
	string ctrlname
	// flag
	variable isok=0
	// get wavename
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	//  globals
	NVAR lastcursAp = root:lastcursAp
	NVAR lastcursBp = root:lastcursBp
	NVAR lastcursAq = root:lastcursAq
	NVAR lastcursBq = root:lastcursBq
	NVAR layer=root:layer
	NVAR frameno=root:frameno
	// cursors
	variable p1,p2,q1,q2
	String proStr=""
	if (stringMatch(num2Str(pcsr(A)),"NaN")!=1 && stringMatch(num2Str(pcsr(B)),"NaN")!=1)
 		p1=pcsr(A)
		p2=pcsr(B)
		q1=qcsr(A)
		q2=qcsr(B)
		// store curs pos
		lastcursAp=pcsr(A)
		lastcursBp=pcsr(B)
		lastcursAq=qcsr(A)
		lastcursBq=qcsr(B)
	else
		p1= lastcursAp
		p2 = lastcursBp
		q1= lastcursAq
		q2 = lastcursBq
	endif
	// profile width
	setGlobals(tmpStr,layer)
	variable halfwidth=getProfileHalfWidth()
	if (halfwidth >= 0)
		isok=getProfile(tmpStr,halfwidth,p1,p2,q1,q2,layer,frameno)
	endif
	return isok
end

Function getProfile(stackStr,halfwidth,p1,p2,q1,q2,layer,i)
	String stackStr
	variable halfwidth,p1,p2,q1,q2,layer,i
	String proStr
	variable isok=0
	if (halfwidth>=0  && stringMatch(num2Str(p1),"NaN")!=1 && stringMatch(num2Str(p2),"NaN")!=1)
		isok=1
		Silent 1; pauseUpdate
		variable np=dimsize($stackStr,0)
		variable nq=dimsize($stackStr,1)
		make/n=(np,nq)/o mylayer=0
		extract2DWaveFromStack($stackStr,layer,mylayer)
		variable npixels=round(sqrt((p2-p1)^2+(q2-q1)^2 ))+1
		if (dimSize($stackStr,2)<=1)
			proStr=stackStr+"_FPR"
		else
			proStr=stackStr+"_FPR_"+num2ScanNumberStr(i)
		endif
		//  create profile wave
		make/n=(npixels)/o/d $proStr 
		wave wpro=$proStr
		wpro=0
		make/n=(npixels)/o pw,qw,pfloor,pceil,qfloor,qceil,c1,c2,c3,c4
		variable theta=-atan2((q2-q1),(p2-p1))
		variable j=-abs(halfwidth)
		do
			pw[]=p1+(p2-p1)/(npixels-1)*p-j*cos(theta-Pi/2)	
			pfloor[] = floor(pw[p])	
			pceil[] =pfloor[p]+1
			qw[]=q1+(q2-q1)/(npixels-1)*p+j*sin(theta-Pi/2)
			qfloor=floor(qw[p])	
			qceil=qfloor[p]+1
			c1[]=(1-abs(pw[p]-pfloor[p]))*(1-abs(qw[p]-qfloor[p]))
			c2[]=(1-abs(pw[p]-pceil[p]))*(1-abs(qw[p]-qfloor[p]))
			c3[]=(1-abs(pw[p]-pfloor[p]))*(1-abs(qw[p]-qceil[p]))
			c4[]=(1-abs(pw[p]-pceil[p]))*(1-abs(qw[p]-qceil[p]))
			wpro[]=wpro[p]+(c1[p]*myLayer[pfloor[p]][qfloor[p]]+c2[p]*myLayer[pceil[p]][qfloor[p]]+c3[p]*myLayer[pfloor[p]][qceil[p]]+c4[p]*myLayer[pceil[p]][qceil[p]])
			// wpro[] = wpro[p]+Interp2D (myLayer, pw[p],qw[p])
			j+=1
		while(j<halfwidth+1)
		wpro/=(halfwidth*2+1)
		killwaves myLayer, pw,qw,pfloor,pceil,qfloor,qceil,c1,c2,c3,c4
		//  draw
		// drawing environment
		SetDrawLayer UserFront
		setDrawEnv xcoord= top,ycoord= left,linefgc= (65280,0,0)		
		DrawLine p1,q1,p2,q2
		// display
		String graphNameStr="LINESCAN_"+stackStr
		bringToFrontAndDisplay(graphNameStr,proStr)
		execute "addProfileControls()"
	endif
	return isok
end

Function getProfileHalfWidth()  //160203
	NVAR lastWidth=root:lastWidth
	variable width=lastWidth
	prompt width, "width",popup "1;3;5;7;9;11;13;15;17;19;21;23;25;27;29;31;33;35;37;39;41;43;45;47;49"
	DoPrompt "Line profile width", width
	width-=1
	if (V_flag==0)
		lastWidth = width+1
	endif
	return width
end

Function shoFramehistogram(ctrlname): buttoncontrol	  //160203
	string ctrlname
	Silent 1; pauseUpdate
	NVAR layer=root:layer
	String imageStr=returnListFirstElement(imageNameList("",""),";")
	imageStr= removeQuotes(imageStr)
	calcFrameHistogram(imageStr)
 end
 
 Function calcFrameHistogram(waveStr)
 	String waveStr
 	NVAR layer=root:layer
 	wave M_ROIMask=root:M_ROIMask
	variable np=dimsize($waveStr,0)
	variable nq=dimsize($waveStr,1)
	variable mp=dimsize(M_ROIMask,0)
	variable mq=dimsize(M_ROIMask,1)
	// graph name 
 	String graphNameStr="HISTOGRAM_"+waveStr
 	// histogram name
 	String tmpStr
	if ((dimSize($waveStr,2)<=1))
		tmpStr=waveStr+"_HST"
	else
		tmpStr=waveStr+"_"+num2ScanNumberStr(layer)+"_HST"
	endif
	// copy the stack layer into "mylayer"
	make/n=(np,nq)/o myLayer
	extract2DWaveFromStack($waveStr,layer,mylayer)
	Wavestats myLayer	
	// mask. if compatible
	if (np==mp && nq==mq)
		myLayer[][]=(M_ROIMask[p][q]==0)*myLayer[p][q]+(M_ROIMask[p][q]==1)*(V_min-1)
	endif
	// 255 channels in hystogram
	make/n=255/o $tmpStr
	redimension/N=(np*nq) myLayer
	Histogram /B={V_min,(V_max-V_min)/254,255} myLayer, $tmpStr
	//
	variable i=0
	do
		if(stringMatch(graphNameStr[i],".")==1)
			graphNameStr=graphNameStr[0,i-1]+"_"+graphNameStr[i+1,strlen(graphNameStr)-1]
		endif
		i+=1
	while(i<strlen(graphNameStr))
	bringToFrontAndDisplay(graphNameStr,tmpStr)
	execute "addHistoScanControls()"
	Killwaves mylayer
 end
 
  Function binning_2x2(ctrlname): buttoncontrol	  //160203
		string ctrlname
		String tmpStr=returnListFirstElement(imageNameList("",""),";")
		tmpStr= removeQuotes(tmpStr)
		variable pp=round(dimSize($tmpStr,0)/2)
		variable qq=round(dimSize($tmpStr,1)/2)
		variable zz=dimSize($tmpStr,2)
		wave w=$tmpStr
		make/n=(pp,qq,zz)/o tmp
		tmp[][][]=w[p*2][q*2][r]+w[p*2+1][q*2][r]+w[p*2][q*2+1][r]+w[p*2+1][q*2+1][r]
		duplicate/o tmp $tmpStr
		killwaves/z tmp
end
 
 Function toUB(ctrlname): buttoncontrol	  //160203
	string ctrlname
		String tmpStr=returnListFirstElement(imageNameList("",""),";")
		tmpStr= removeQuotes(tmpStr)
		rescale(tmpStr)
		redimension/B/U $tmpStr
		WMCreateImageParticlePanel();
end

Function doPA(ctrlname): buttoncontrol	  //160203
	string ctrlname
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	// performs check on dimensions
	String infoStr=WaveInfo($tmpStr,0)
	variable type=NumberByKey("NUMTYPE", infoStr)	
	if (type != 72)
		rescale(tmpStr)
		redimension/B/U $tmpStr
	endif
	//
	NVAR pa_thresh_low=root:pa_thresh_low
	NVAR pa_thresh_high=root:pa_thresh_high
	NVAR pa_nlevels=root:pa_nlevels
	variable nlevels=pa_nlevels
	variable thrs12=pa_thresh_low
	variable thrs13=pa_thresh_high
	prompt thrs12, "lower threshold"
	prompt thrs13, "upper threshold"
	prompt nlevels, "levels",popup "two;three;"
	DoPrompt "threshold", thrs12,thrs13,nlevels
	if (V_flag==0)
		if (nlevels==1)
			thrs13=thrs12
		endif
		// 
		variable grey,dark,bright
		variable levelgrey,leveldark,levelbright
		variable particlearea=30
		//
		redimension/B/U $tmpStr
 		variable np=dimsize($tmpStr,0)
 		variable nq=dimsize($tmpstr,1)
 		make/n=(np,nq)/o tmplayer=0
		make/n=(np,nq)/o tmpmask=0
		duplicate/o root:M_ROIMask neutralmask
		imagestats neutralmask
		variable roiarea=np*nq*(1-V_avg)
 		// analyse dark
 		duplicate/o $tmpStr tmplayer,tmpmask
 		ImageThreshold/o/Q/T=(thrs12)/M=0/R=M_ROIMask tmpmask 
		ImageAnalyzeParticles/E/M=3/Q/A=(particlearea)/O/W stats, tmpmask
		// calculates the avg value within mask
		duplicate/o root:M_Particle tmpmask
		tmpmask-=16		
		imageStats/M=1/R=tmpmask tmplayer
		doupdate
		leveldark=V_avg
		// total area
		dark=sum(root:W_ImageObjArea, 0, numpnts(root:W_ImageObjArea) )/roiarea
		// start creating mask for grey
		duplicate/o M_ROIMask neutralmask
		neutralmask*=(208)
		neutralmask+=(tmpmask[p][q]-48)
		// analyse bright
		//  copy layer into tmp
		duplicate/o $tmpStr tmplayer,tmpmask
		// create and append overlay
		//wave wroi=M_ROiMask
		//wroi=0
		//wover[][]+=wroi[p][q]
		//CheckDisplayed $overlayYStr
		// if( V_Flag==0 )
				//AppendImage M_particle
				//ModifyImage ''#1  explicit=1, eval={64,65000,128,128}	
		//endif
		// bright
		ImageThreshold/O/I/Q/T=(thrs13)/M=0/R=M_ROIMask tmpmask 
		ImageAnalyzeParticles/E/M=3/Q/A=(particlearea)/O stats, tmpmask
		duplicate/o root:M_Particle tmpmask
		tmpmask-=16
		imageStats/M=1/R=tmpmask tmplayer
		levelbright=V_avg
		// total area
		bright=sum(root:W_ImageObjArea, 0, numpnts(root:W_ImageObjArea) )/roiarea
		// neutral
		String textStr=""
		if (nlevels==2)
			neutralmask+=(tmpmask[p][q]-48)
			imageStats/M=1/R=neutralmask tmplayer
			levelgrey=V_avg
			ImageAnalyzeParticles/E/W/Q/M=3/A=(particlearea) stats, neutralmask
			// total area
			grey=sum(root:W_ImageObjArea, 0, numpnts(root:W_ImageObjArea) )/roiarea
 			print "three level *************************************************"
 			print "dark","grey","bright","leveldark","levelgrey","levelbright"
			print dark,grey,bright,leveldark,levelgrey,levelbright
			textStr = "dark: "+num2Str(dark)+" at "+num2Str(leveldark)+"\r"+ "grey: "+num2Str(grey)+" at "+num2Str(levelgrey)+"\r"+"bright: "+num2Str(bright)+" at "+num2Str(levelbright)
 		else
 			print "two level *************************************************"
 			print "dark","bright","leveldark","levelbright"
			print (dark+(1-bright))/2,(bright+(1-dark))/2,NaN,leveldark,levelbright,NaN
			textStr = "dark: "+num2Str( (dark+(1-bright))/2)+" at "+num2Str(leveldark)+"\r"+ "bright: "+num2Str((bright+(1-dark))/2)+" at "+num2Str(levelbright)
 		endif
 		// overlay
 		String overlayXStr=tmpStr+"_OVRX"
 		String overlayYStr=tmpStr+"_OVRY"
		duplicate/o  root:W_BoundaryY $overlayYStr
		duplicate/o  root:W_BoundaryX $overlayXStr
		CheckDisplayed $overlayYStr
		if( V_Flag==0 )
				AppendTograph/L/T $overlayYStr vs $overlayXStr
		endif
 		// textbox
 		textBox/C/N=text0/A=RB textstr
		pa_thresh_low=thrs12
		pa_thresh_high=thrs13
		pa_nlevels=nlevels
		killwaves/z tmplayer,tmpmask,neutralmask
	endif	
end

Function Imagelayout(ctrlname): buttoncontrol	
	String ctrlname
	NVAR xmcdupper=root:xmcdupper
	NVAR xmcdlower=root:xmcdlower
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	ImageStats/R=M_RoiMask $tmpStr
	variable pp=dimsize($tmpStr,0)
	variable qq=dimsize($tmpStr,1)
	variable rpp=dimsize(M_ROiMask,0)
	variable rqq=dimsize(M_ROiMask,1)
	variable mini=min(V_min,0)
	if ((pp==rpp) &&(qq==rqq))
		MCPmask($tmpStr,M_ROiMask,mini)
	endif
end

Function Stacklayout(): GraphMarquee	
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	ImageStats/R=M_RoiMask $tmpStr
	variable pp=dimsize($tmpStr,0)
	variable qq=dimsize($tmpStr,1)
	variable zz=dimsize($tmpStr,2)
	variable rpp=dimsize(M_ROiMask,0)
	variable rqq=dimsize(M_ROiMask,1)
	variable mini=min(V_min,0)
	if ((pp==rpp) &&(qq==rqq))
		variable i=0
		make/n=(pp,qq)/o tmp
		for (i=0; i<zz; i+=1)
			setGlobals(tmpStr,i)
			extract2DWaveFromStack($tmpStr,i,tmp)			
			MCPmask(tmp,M_ROiMask,mini)
			insert2DWaveinStack($tmpStr,i,tmp)
		endfor
		killwaves/z tmp
	endif
end

Function SubtractFrameZeroFromStack(): GraphMarquee	
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	variable pp=dimsize($tmpStr,0)
	variable qq=dimsize($tmpStr,1)
	variable zz=dimsize($tmpStr,2)
	variable rpp=dimsize(M_ROiMask,0)
	variable rqq=dimsize(M_ROiMask,1)
	if ((pp==rpp) &&(qq==rqq))
		variable i=0
		make/n=(pp,qq)/o tmp
		extract2DWaveFromStack($tmpStr,0,tmp)			
		for (i=0; i<zz; i+=1)
			setGlobals(tmpStr,i)
			subtractImageFromStack($tmpStr,tmp,i)
		endfor
		killwaves/z tmp
	endif
end

Function AverageValue() : GraphMarquee 
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	GetMarquee/Z/K left, top
	if (V_Flag == 0)
		Print "There is no marquee"
	else
		variable p0=round(V_top)
		variable p1=round(V_bottom)
		variable q0=round(V_left)
		variable q1=round(V_right)
		variable np=dimsize($tmpStr,0)
		variable nq=dimsize($tmpStr,1)
		variable avg
		duplicate/o $tmpStr mymask
		mymask=1
		redimension/B/U mymask
		mymask[q0,q1][p0,p1]=0 
		ImageStats/R=mymask $tmpStr
		avg=V_avg
		print "the avg value of: ",tmpStr," is:",avg
		killwaves mymask
	endif
end

Proc resizeImageWindow(ctrlname): buttoncontrol	
	String ctrlname
	AutoSizeImage(,1)
End

Function killseries(ctrlname): buttoncontrol	
	String ctrlname
	String stackStr=returnListFirstElement(imageNameList("",""),";")
	String wavelistStr = WaveList("*", ";","WIN:") 
	// kill image stack
	stackStr= removeQuotes(stackStr)
	String graphStr=WinName(0,1)
	print "killing...",graphStr," and ",stackStr
	DoWindow/k $graphStr
	killwaves/z $stackStr
	// kill traces
	variable n=ItemsInList(wavelistStr) 
	variable m=0
	String wavenameStr = ""
	do
		wavenameStr=StringFromList(m, wavelistStr) 
		wave w = $wavenameStr
		killwaves w
		print "killing....",wavenameStr
		m+=1
	while(m < n)
End

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// parameters functions
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Function getPhaseSign(fileNameStr) //-1,+1, 0 err
	String fileNameStr
	Variable refNum	
	Open/P=imagePath/R/Z=1 refNum as fileNameStr
	if (stringMatch(S_fileName,"")!=1)
		String str=""
		FReadLine /T=(num2char(13)) refNum, str		// Read all into string variable
		return sign(str2num(getBcsParameter(str,"Undulator 1 phase")))
	else
		return 0
	endif
	Close refNum
End

Function getPhase_old(fileNameStr) // obsolete, use with data older than 2002/04
	String fileNameStr
	Variable refNum	
	Open/R/Z=1 refNum as fileNameStr
	//print "opening file: ",returnImageName(S_fileName), "of type: ", returnImageType(S_fileName)
	if (stringMatch(S_fileName,"")!=1)
		String str=""
		FReadLine /T=(num2char(13)) refNum, str		// Read all into string variable
		return str2num(getBcsParameter(str,"phase1"))
	else
		return NaN
	endif
	Close refNum
End

Function getTemperature(fileNameStr) //
	String fileNameStr
	Variable refNum	
	Open/P=imagePath/R/Z=1 refNum as fileNameStr
	//print "opening file: ",returnImageName(S_fileName), "of type: ", returnImageType(S_fileName)
	if (stringMatch(S_fileName,"")!=1)
		String str=""
		FReadLine /T=(num2char(13)) refNum, str		// Read all into string variable
		return str2num(getLEEMParameter(str,"Sample Temp."))
	else
		return NaN
	endif
	Close refNum
End

Function xmcdAvgTemperature()
	SVAR imagePathStr=root:imagePathStr
	SVAR imagetypeStr=root:imageTypeStr
	SVAR scanNameStr=root:scanNameStr
	SVAR scanNumberStr=root:scanNumberStr
	SVAR separatorStr=root:separatorStr
	variable temperature0,temperature1
	String fileNameStr0,fileNameStr1
	// store variables 
	String oldimagePathStr=imagePathStr
	String oldscanNameStr=scanNameStr
	String oldscanNumberStr=scanNumberStr
	// _TEMP wave	
	variable len=strlen(scanNameStr)
	String  tempStr=scanNameStr[0,len-4]+"_TEMP"
	//String  infoStr=scanNameStr[0,len-4]+"_INFO"
	Make/n=999/o $tempStr=NaN  //, $infoStr=p
	wave tempw=$tempStr
	variable imax=0
	do
		fileNameStr0 = scanNameStr +"_leem"+separatorStr+num2ScanNumberStr(0)+"."+"txt"
		fileNameStr1 = scanNameStr +"_leem"+separatorStr+num2ScanNumberStr(1)+"."+"txt"
		temperature0 = getTemperature(fileNameStr0)  
		temperature1 = getTemperature(fileNameStr1)  
		imax=Str2num(scanNumberStr)
		tempw[imax]=(temperature0+temperature1)/2
	while(incrementPath()==0)
	redimension/n=(imax) tempw  //,$infoStr
	// reset gen variables
	imagePathStr=oldimagePathStr
	scanNameStr=oldscanNameStr
	scanNumberStr=oldscanNumberStr
end

 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// image fft functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Function imageFFT(ctrlname): buttoncontrol	
	String ctrlname
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	// image dimensions
	wave w=$tmpStr
	variable np=dimSize(w,0)
	variable nq=dimSize(w,1)
	// introducing subset waves for fft - correlation
	//variable mp=calcPower2(np)
	//variable mq=calcPower2(nq)
	variable mp=(round(np/2)*2)
	variable mq=(round(np/2)*2)
	print mp,mq
	String fftStr=tmpStr+"_FT"
	wave wfft=$fftStr
	make/n=(mp,mq)/o $fftStr, hann
	variable mr=round(np/2-mp/2)
	variable mc=round(nq/2-mq/2)  // +nq/12 burn in MCP
	wfft=w[mr+p][mc+q]
	wavestats/Q wfft
	wfft-=V_avg
	// windowing (for making nice FFTs)
	hann=sin(pi*p/mp)*sin(pi*q/mq)
	wfft*=hann
	wave w=$fftStr
	FFT/OUT=3/Z/DEST=tfft $fftStr
	SetScale/P x 0,1,"", tfft;DelayUpdate
	SetScale/P y 0,1,"", tfft
	redimension/R tfft
	duplicate/o tfft $fftStr
	 //w=sqrt(real(w)^2+imag(w)^2)
	 killwaves/z hann,tfft
	String graphNameStr=removeForbiddenChar("FFT_"+tmpStr)
	BringImageToFront(graphNameStr,fftstr)
	Execute "addBWControls()"
	//SetScale/I y -0.5,0.5,"1/pixel", wfft
	//SetScale/I x 0,0.5,"1/pixel", wfft
end

Function interpolateimage(): graphMarquee
	String ctrlname
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	// zero binning
	variable bin=3
	// image dimensions
	wave w=$tmpStr
	variable np=dimSize(w,0)
	variable nq=dimSize(w,1)
	variable nlayer=dimSize(w,2)
	variable maxpq=max(np,nq)
	variable flag=1
	prompt flag, "hanning centre",popup "zero;f/2"
	prompt bin, "zero binning",popup "two;four;eight"
	DoPrompt "fft interpolate", flag,bin
	bin=2^bin
	PauseUpdate
	if (nlayer < 2 && V_flag ==0)	
		// introducing subset waves for fft - correlation
		variable mp=2*calcPower2(maxpq)
		String fftStr=tmpStr+"_FT"
		make/n=(mp,mp)/o $fftStr=0
		wave wfft=$fftStr
		wavestats/Q w
		variable mr=round(mp/2-np/2)
		variable mc=round(mp/2-nq/2)  
		wfft[mr,mr+np][mc,mc+nq]=(w[p-mr][q-mc]-V_avg)*sin(pi*(p-mr)/np)*sin(pi*(q-mc)/nq)
		FFT wfft
		np=dimSize(wfft,0)
		nq=dimSize(wfft,1)
		variable newp=(np-1)*bin+1
		variable newq=nq*bin
		String ifftStr=tmpStr+"_IFT"
		make/n=(newp,newq)/o/c $ifftStr=cmplx(0,0)
		wave/c wi=$ifftStr
		mr=0
		mc=round(newq/2-nq/2)
		if (flag == 1)
			wi[mr,mr+np][mc,mc+nq]=cmplx(real(wfft[p-mr][q-mc])*cos(0.5*pi*(p-mr)/np)*sin(pi*(q-mc)/nq),imag(wfft[p-mr][q-mc]*cos(0.5*pi*(p-mr)/np)*sin(pi*(q-mc)/nq)))
		endif
		if (flag == 2)
			wi[mr,mr+np][mc,mc+nq]=cmplx(real(wfft[p-mr][q-mc])*sin(pi*(p-mr)/np)*sin(pi*(q-mc)/nq),imag(wfft[p-mr][q-mc]*sin(pi*(p-mr)/np)*sin(pi*(q-mc)/nq)))
		endif
		IFFT wi
		String graphNameStr=removeForbiddenChar("IFT_"+tmpStr)
		BringImageToFront(graphNameStr,ifftstr)
		Execute "addBWControls()"
		SetScale/P x 0,1,"", wi; DelayUpdate
		SetScale/P y 0,1,"", wi
		Execute "autosizeImage(0.166,1)"
		doUpdate
	endif
end
 
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// xmcd functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Function shoLoop(ctrlname): buttoncontrol	
	String ctrlname
	NVAR layer = root:layer
	String sourceStr=returnListFirstElement(imageNameList("",""),";")
	sourceStr= removeQuotes(sourceStr)
	variable i=0
	do
		modifyImage $sourceStr plane=0
		layer = 0
		doUpdate
		modifyImage $sourceStr plane=1
		layer = 1
		doUpdate
		i+=1
	while(i<15)
End

Function rexmcd(ctrlname): buttoncontrol	
	String ctrlname
	String stackStr=returnListFirstElement(imageNameList("",""),";")
	stackStr= removeQuotes(stackStr)
	String xmcdStr=stackStr +"_"+"XMCD"
	variable np=dimSize($stackStr,0)
	variable nq=dimSize($stackStr,1)
	make/n=(np,nq,1)/o/u $xmcdStr
	calcXMCD($stackStr,$xmcdStr,-1) 
	// display
 	NewImage/S=2 $xmcdStr
	ModifyImage $xmcdStr plane=0
	if ((stringMatch(num2Str(np),"NaN")!=1)&&(stringMatch(num2Str(nq),"NaN")!=1)) 
		Execute "addXMCDControls()"
		ModifyGraph width=250*np/nq,height=250
		TextBox/C/N=text0/A=LT xmcdStr
		DoUpdate
		ModifyGraph width=0,height=0
		ModifyImage $xmcdstr ctab= {-0.3,0.3,Grays,0}
	endif	
end

Function NormaliseXMCDbackground() : GraphMarquee //160203
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	GetMarquee/Z/K left, top
	if (V_Flag == 0)
		Print "There is no marquee"
	else
		variable p0=round(V_top)
		variable p1=round(V_bottom)
		variable q0=round(V_left)
		variable q1=round(V_right)
		variable np=dimsize($tmpStr,0)
		variable nq=dimsize($tmpStr,1)
		variable avg0,avg1
		make/n=(np,nq)/O mystack0,myStack1
		extract2DWaveFromStack($tmpStr,0,mystack0)
		extract2DWaveFromStack($tmpStr,1,mystack1)
		ImageStats/G={p0,p1,q0,q1}/M=1 myStack0
		avg0=V_avg
		ImageStats/G={p0,p1,q0,q1}/M=1 myStack1
		avg1=V_avg
		print avg0, avg1
		MultiplyStackLayer($tmpStr,avg1/avg0,0)
		killwaves mystack0,mystack1
		print "rescaling factor:", avg1/avg0
	endif
End

Function XMCDlayout(ctrlname): buttoncontrol	
	String ctrlname
	NVAR xmcdupper=root:xmcdupper
	NVAR xmcdlower=root:xmcdlower
	NVAR imageRange=root:imageRange
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	variable  ir=imageRange
	prompt ir, "image range"
	DoPrompt "input image range", ir
	if (V_Flag==0)
		// average value to zero
		variable avg2zeroFlag=2
		prompt avg2zeroFlag, "subtract average value?",popup "yes;no"
		DoPrompt "xmcd", avg2zeroFlag
		if (avg2zeroFlag==1 && V_Flag==0)
			ImageStats/R=M_RoiMask $tmpStr
			subtractOffset($tmpStr,(V_avg))
		endif
		MCPmask($tmpStr,M_ROiMask,0)
		xmcdlower=-ir
		xmcdupper=+ir
		imageRange=ir
		ModifyImage $tmpStr ctab= {xmcdlower,xmcdupper,,0}
	endif
end

Function bkgRemove(ctrlname): buttoncontrol	
	String ctrlname
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	redimension/D $tmpStr
	Imagestats/R=M_RoiMask $tmpStr
	duplicate/o M_RoiMask tmpMask
	tmpMask-=1
	ImageRemoveBackground /R=tmpMask/O/P=3 $tmpStr
	killwaves/z tmpMask
	subtractOffset($tmpStr,(-V_avg))
End

Function MCPmask(w1,w2,k)
	wave w1
	wave w2
	variable k
	w1=(1-w2)*w1+abs(w2)*k
end

Function despikeXMCD(ctrlname): buttoncontrol	
	String ctrlname
	NVAR xmcdupper=root:xmcdupper
	NVAR xmcdlower=root:xmcdlower
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	wave w=$tmpStr
 	Variable lowerThreshold=xmcdlower,upperThreshold=xmcdupper
	Prompt lowerThreshold, "Enter Lower Threshold: "
	Prompt upperThreshold, "Enter Upper Threshold: "
	DoPrompt "Set threshold", lowerThreshold, upperThreshold
	if (V_flag==0)
 		w=w*(w<upperThreshold)+upperThreshold*(w>=upperThreshold)
 		w=w*(w>lowerThreshold)+lowerThreshold*(w<=lowerThreshold)
 		w[0,2][0,4]=lowerThreshold
 		w[0,2][2,4]=upperThreshold
 		ModifyImage $tmpStr ctab= {lowerThreshold,upperThreshold,,0}
 	endif
End

Function save_XMCD_Jpeg(ctrlname): buttoncontrol	
	String ctrlname
	NVAR xmcdupper=root:xmcdupper
	NVAR xmcdlower=root:xmcdlower
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	despikeXMCD("")
	saveJpeg("")
End

Function makeAll(ctrlname): buttoncontrol	
	String ctrlname
	NVAR xmcdupper=root:xmcdupper
	NVAR xmcdlower=root:xmcdlower
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	wave w=$tmpStr
	bkgRemove("")
	XMCDlayout("")
	prevSubset("")
	despikeXMCD("")
	saveJpeg("")
End

Function removeBackGround(w)
	wave w
	redimension/D w
	if (exists("M_RoiMask")==1)
		ImageStats/R=M_RoiMask w
		print V_avg
		duplicate/o M_RoiMask tmpMask
		tmpMask-=1
		ImageRemoveBackground /R=tmpMask/O/P=3 w
		killwaves/z tmpMask
		subtractOffset(w,(-V_avg))
	endif
End

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// correlation functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Function/C autoCalcStackCorrelation(stackStr,layer1,layer2,filterFlag)
	String stackStr
	variable layer1,layer2,filterFlag
	// initialises pdrift, qdrift ....
	variable pdrift=0
	variable qdrift=0
	variable this_pdrift,this_qdrift
	variable pmin=0,qmin=0,vmin=+inf
	// stack dimensions
	wave wstack=$stackStr
	variable np=dimSize(wstack,0)
	variable nq=dimSize(wstack,1)
	// introducing subset waves for fft - correlation
	variable mp=calcPower2(np/2)
	variable mq=calcPower2(nq)
	make/n=(mp,mq)/o/d w0,w1,hann
	variable mr=round(np/2-mp/2)
	variable mc=round(nq/2-mq/2+nq/12)
	w0=wstack[mr+p][mc+q][layer1]
	w1=wstack[mr+p][mc+q][layer2]
	wavestats/Q w0
	w0-=V_avg
	wavestats/Q w1
	w1-=V_avg
	//sharpening
	if (filterFlag == 2) 
		imageFilter findEdges w0 //ImageFilter sharpen w0
		ImageFilter findEdges w1 //ImageFilter sharpen w1
	endif
	if (filterFlag == 3) 
		ImageFilter sharpen w0
		ImageFilter sharpen w1
	endif
	// windowing
	hann=sin(pi*p/mp)*sin(pi*q/mq)
	w0*=hann
	w1*=hann
	// duplicate/o w1 wdebug
	// FFT
	FFT w0
	FFT w1
	// 2D correlation
	wave/C wc0=w0
      	wave/C wc1=w1
      	wc0=conj(wc0)
      	wc0*=wc1
	IFFT wc0
	// duplicate/o wc0 wdebug
	ImageStats wc0
      	// print V_MaxrowLoc,mp,V_MaxColLoc,mq
      	// calc drift for p,q
      	if (V_MaxColLoc >= mq/2)
       	qdrift = V_MaxColLoc -mq 
      	else
      		qdrift=V_MaxColLoc     
 	endif
 	if (V_MaxRowLoc >= mp/2)
       	pdrift = V_MaxRowLoc-mp
	else      
      		pdrift=V_MaxRowLoc     
      endif
      killwaves/z w0,w1,wc0,wc1,hann
	return cmplx(pdrift,qdrift)
end

Function/C autoImageCorrelation(imageStr0,imageStr1,filterFlag)
	String imageStr0,imageStr1
	variable filterFlag
	// initialises pdrift, qdrift ....
	variable pdrift=0
	variable qdrift=0
	variable this_pdrift,this_qdrift
	variable pmin=0,qmin=0,vmin=+inf
	// stack dimensions
	wave wimage0=$imageStr0
	wave wimage1=$imageStr1
	variable np=dimSize(wimage1,0)
	variable nq=dimSize(wimage1,1)
	// introducing subset waves for fft - correlation
	variable mp=calcPower2(np/2)
	variable mq=calcPower2(nq)
	make/n=(mp,mq)/o/d w0,w1,hann
	variable mr=round(np/2-mp/2)
	variable mc=round(nq/2-mq/2+nq/12)  // +nq/12 burn in MCP
	w0=wimage0[mr+p][mc+q]
	w1=wimage1[mr+p][mc+q]
	wavestats/Q w0
	w0-=V_avg
	wavestats/Q w1
	w1-=V_avg
	//sharpening
	if (filterFlag == 2) 
		imageFilter findEdges w0 //ImageFilter sharpen w0
		ImageFilter findEdges w1 //ImageFilter sharpen w1
	endif
	if (filterFlag == 3) 
		ImageFilter sharpen w0
		ImageFilter sharpen w1
	endif
	// windowing
	// windowing (for making nice FFTs)
	hann=sin(pi*p/mp)*sin(pi*q/mq)
	w0*=hann
	w1*=hann
	// FFT
	FFT w0
	FFT w1
	// 2D correlation
	wave/C wc0=w0
      	wave/C wc1=w1
      	wc0=conj(wc0)
      	wc0*=wc1
	IFFT wc0
	ImageStats wc0
      	// print V_MaxrowLoc,mp,V_MaxColLoc,mq
      	// calc drift for p,q
      	if (V_MaxColLoc >= mq/2)
       	qdrift = V_MaxColLoc -mq 
      	else
      		qdrift=V_MaxColLoc     
 	endif
 	if (V_MaxRowLoc >= mp/2)
       	pdrift = V_MaxRowLoc-mp
	else      
      		pdrift=V_MaxRowLoc     
      endif
      killwaves/z w0,w1,wc0,wc1,hann
	return cmplx(pdrift,qdrift)
end

Function/C marqueeCorrelation(imageStr0,imageStr1,filterFlag)
	String imageStr0,imageStr1
	variable filterFlag
	// initialises pdrift, qdrift ....
	variable pdrift=0
	variable qdrift=0
	variable this_pdrift,this_qdrift
	variable pmin=0,qmin=0,vmin=+inf
	// image dimensions - assumes input has matching dimensions
	wave wimage0=$imageStr0
	wave wimage1=$imageStr1
	variable np=dimSize(wimage1,0)
	variable nq=dimSize(wimage1,1)
	wavestats/Q wimage0
	variable V_avg0=V_avg
	wavestats/Q wimage1
	variable V_avg1=V_avg
	// introducing subset waves for fft - correlation
	variable mp=2*calcPower2(np)
	variable mq=2*calcPower2(nq)
	make/n=(mp,mq)/o w0=0,w1=0
	variable mr=round(mp/2-np/2)
	variable mc=round(mq/2-nq/2)  
	w0[mr,mr+np][mc,mc+nq]=(wimage0[p-mr][q-mc]-V_avg0)*sin(pi*(p-mr)/np)*sin(pi*(q-mc)/nq)
	w1[mr,mr+np][mc,mc+nq]=(wimage1[p-mr][q-mc]-V_avg1)*sin(pi*(p-mr)/np)*sin(pi*(q-mc)/nq)
	// filtering
	if (filterFlag == 2) 
		imageFilter findEdges w0 //ImageFilter sharpen w0
		ImageFilter findEdges w1 //ImageFilter sharpen w1
	endif
	if (filterFlag == 3) 
		ImageFilter sharpen w0
		ImageFilter sharpen w1
	endif
	doupdate
	// FFT
	FFT w0
	FFT w1
	// 2D correlation
	wave/C wc0=w0
      	wave/C wc1=w1
      	wc0=conj(wc0)
      	wc0*=wc1
	IFFT wc0
	ImageStats wc0
      	// calc drift for p,q
      	if (V_MaxColLoc >= mq/2)
       	qdrift = V_MaxColLoc -mq 
      	else
      		qdrift=V_MaxColLoc     
 	endif
 	if (V_MaxRowLoc >= mp/2)
       	pdrift = V_MaxRowLoc-mp
	else      
      		pdrift=V_MaxRowLoc     
      	endif
     	killwaves/z w0,w1,wc0,wc1
	return cmplx(pdrift,qdrift)
end

function calcPower2(x)
	variable x
	variable p2=2
	do
		p2=p2*2	
	while (p2 < x)
	return p2/2
end

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// drift correction functions and co-addition
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Function driftCorrection(ctrlname): buttoncontrol	
	String ctrlname
	String stackStr=returnListFirstElement(imageNameList("",""),";")
	stackStr= removeQuotes(stackStr)
	variable filterFlag=2
	prompt filterFlag, "perform image filtering when calculating image correlation?",popup "none;find edges;sharpen"
	DoPrompt "drift correction", filterFlag
	print "calculating correlation for: ", stackStr
	variable/C pq=autoCalcStackCorrelation(stackStr,0,1,filterFlag)
	variable pdrift=real(pq)
	variable qdrift=imag(pq)
	Prompt pdrift, "suggested column drift correction"
	Prompt qdrift, "suggested row drift correction"
	DoPrompt "drift correction", pdrift, qdrift
	if (V_flag == 1)
		pdrift=0
		qdrift=0
	endif
	String xmcdStr=stackStr +"_"+"XMCD"
	wave wstack=$stackStr
	variable np=dimSize(wstack,0)
	variable nq=dimSize(wstack,1)
	make/n=(np,nq,1)/o/u $xmcdStr
	calcXMCDdrift(wstack,$xmcdStr,-1,pdrift,qdrift)
End

Function shift2DriftCorrect(ctrlname): buttoncontrol
	String ctrlname
	String stackStr=returnListFirstElement(imageNameList("",""),";")
	stackStr= removeQuotes(stackStr)
	// stack dimensions
	wave wstack=$stackStr
	variable np=dimSize(wstack,0)
	variable nq=dimSize(wstack,1)
	variable nlayers=dimSize(wstack,2)
	// drift
	variable/C pqdrift
	variable pp=0
	variable qq=0
	variable ilayer
	variable niter=5
	variable jiter=0
	variable sumppqq=0
	// coadd
	make/n=(np,nq)/o/d wtmp1,wtmp2
	wtmp1=0
	addShiftedImage(wtmp1,wstack,0,0,0)
	// input correlation alghoritm
	variable filterFlag=0
	prompt filterFlag, "perform image filtering when calculating image correlation?",popup "none;find edges;sharpen"
	DoPrompt "drift correction", filterFlag
	if (V_Flag == 0)
	// loop
	do
		ilayer=1 
		sumppqq=0
		do
			// extract layer
			extract2DWaveFromStack(wstack,ilayer,wtmp2)  
			// calc corr between layer zero and i
			pqdrift=autoImageCorrelation("wtmp1","wtmp2",filterFlag)
			pp=real(pqdrift)
			qq=imag(pqdrift)
			print "iteration ",jiter,"shift by:",stackStr,"layer",ilayer,"drift correction:",pp,qq
			// update tmp1
			addShiftedImage(wtmp1,wstack,ilayer,pp,qq)
			// shift
			ShiftStackLayer(wstack,ilayer,pp,qq)
			sumppqq=sumppqq+abs(pp)+abs(qq)
			ilayer+=1
		while(ilayer<=nlayers)
		jiter+=1
	while (sumppqq>1)
	endif
	killwaves/z wtmp1,wtmp2
	doUpdate
End

Function marqueeDriftCorrection(): GraphMarquee 
	String ctrlname
	GetMarquee/Z/K left, top
	if (V_Flag == 0)
		Print "There is no marquee"
	else
		// input correlation alghoritm
		variable filterFlag=0
		prompt filterFlag, "perform image filtering when calculating image correlation?",popup "none;find edges;sharpen"
		DoPrompt "drift correction", filterFlag
		if (V_Flag == 0)
			// marquee position
			variable q0=round(V_top)
			variable q1=round(V_bottom)
			variable p0=round(V_left)
			variable p1=round(V_right)
			// stack dimensions  
			String stackStr=returnListFirstElement(imageNameList("",""),";")
			stackStr= removeQuotes(stackStr)
			wave wstack=$stackStr
			variable np=dimSize(wstack,0)
			variable nq=dimSize(wstack,1)
			variable nlayers=dimSize(wstack,2)
			make/n=(p1-p0,q1-q0)/o tmpw0,tmpwi
			// layer 0 to tmpw
			extractMarqueeFromStack(wstack,0,tmpw0,p0,q0)
			// drift loop
			variable ilayer=1
			variable/c pqdrift
			variable pp,qq
			variable sumppqq,jiter=0
			do
			ilayer=1 
			sumppqq=0
			extractMarqueeFromStack(wstack,0,tmpw0,p0,q0)
				do
					extractMarqueeFromStack(wstack,ilayer,tmpwi,p0,q0)
					// calc corr between layer zero and i
					pqdrift=marqueeCorrelation("tmpw0","tmpwi",filterFlag)
					pp=real(pqdrift)
					qq=imag(pqdrift)
					// shift
					print "iteration:",jiter,stackStr," - layer: ",ilayer," drift correction:",pp,qq
					ShiftStackLayer(wstack,ilayer,pp,qq)
					// coadd to tmp
					tmpw0[][]+=wstack[p0+p][q0+q][ilayer]
					sumppqq=sumppqq+abs(pp)+abs(qq)
					ilayer+=1
				while (ilayer < nlayers)
			jiter+=1
			while (sumppqq>1 && jiter < 5)	
			killwaves/z tmpwi,tmpw0
		endif
	endif
end

Function driftcorrect_and_savetiff(): GraphMarquee 
	String ctrlname
	SVAR imagetypeStr=root:imageTypeStr
	SVAR imagePathStr=root:imagePathStr
	GetMarquee/Z/K left, top
	// marquee position
	variable q0=round(V_top)
	variable q1=round(V_bottom)
	variable p0=round(V_left)
	variable p1=round(V_right)
	// stack dimensions  
	if (V_Flag == 0)
		Print "There is no marquee"
	else
		// wave name
		String nameStr=returnListFirstElement(imageNameList("",""),";")
		nameStr= removeQuotes(nameStr)
		// input correlation alghoritm
		variable filterFlag=0
		prompt filterFlag, "perform image filtering when calculating image correlation?",popup "none;find edges;sharpen"
		DoPrompt "drift correction", filterFlag
		String tmpStr,outStr
		if (V_Flag == 0)
			// reference image
			make/n=(p1-p0,q1-q0)/o tmpw0,tmpwi
			extractMarqueeFromStack($nameStr,0,tmpw0,p0,q0)
			newimage/S=2 tmpw0
			nameStr = nameStr[0,strlen(nameStr)-4]
			variable i=0
			variable j=0
			variable/C pqdrift
			variable pp,qq
			do 
				// open image
				if (i<10)
					tmpStr = namestr+"00"+num2str(i)+"."+imageTypeStr 
				else
					if (i<100)
						tmpStr = namestr+"0"+num2str(i)+"."+imageTypeStr 
					else
						tmpStr = namestr+num2str(i)+"."+imageTypeStr 
					endif
				endif
				tmpStr = OpenImage(tmpStr,2)
				// drift correct
				print "performing drift correction"
				j=0
				do
					
					extractMarqueeFromStack($tmpStr,0,tmpwi,p0,q0)
					pqdrift=marqueeCorrelation("tmpw0","tmpwi",filterFlag)
					pp=real(pqdrift)
					qq=imag(pqdrift)
					ShiftStackLayer($tmpStr,0,pp,qq)
					print "iteration ",j," drift correction:",pp,qq
					j+=1
				while (pp!=0 && qq!=0 && j<5)
				// save image
				redimension/U/W $tmpStr
				outStr=tmpStr+".tif"
				ImageSave/D=16/O/T="TIFF"/P=imagePath  $tmpStr as outStr  //save
				print "data saved in: ",imagePathStr," as ", outStr
				// update loop
				i+=1
			while (stringMatch(tmpStr,"")!=1)		
		endif
	endif
end


Function manualDrift(ctrlname): buttoncontrol
	String ctrlname
	String stackStr=returnListFirstElement(imageNameList("",""),";")
	stackStr= removeQuotes(stackStr)
	NVAR layer=root:layer
	variable p=0
	prompt p, "drift correction for rows: "
	variable q=0
	prompt q, "drift correction for columns: "
	DoPrompt "manual drift correction", p,q
	if (V_flag == 0)
		wave wstack=$stackStr
		ShiftStackLayer(wstack,layer,p,q)
		doUpdate
	endif
end

Function smartCoadd(ctrlname): buttoncontrol
	String ctrlname
	String stackStr=returnListFirstElement(imageNameList("",""),";")
	stackStr= removeQuotes(stackStr)
	// stack dimensions
	wave wstack=$stackStr
	variable np=dimSize(wstack,0)
	variable nq=dimSize(wstack,1)
	variable nlayers=dimSize(wstack,2)
	// coadded wave
	String coaddStr=stackStr+"_CO"
	make/n=(np,nq)/o $coaddStr
	wave wcoadd=$coaddStr
	wcoadd=0
	redimension/d wcoadd
	duplicate/o wcoadd tmpw
	// correlation ??
	variable corrFlag=1
	prompt corrFlag, "calculate correlation and perform drift correction",popup "no;yes"
	DoPrompt "co-add", corrFlag
	if (V_Flag==0)
		variable filterFlag=0
		if (corrFlag == 2)
		// input correlation alghoritm
			prompt filterFlag, "perform image filtering when calculating image correlation?",popup "none;find edges;sharpen"
			DoPrompt "drift correction", filterFlag
		endif
		// drift
		variable/C pqdrift
		variable pp=0
		variable qq=0
		// loop
		variable ilayer=0
		addShiftedImage(wcoadd,wstack,0,0,0)
		do
			// extract layer
			extract2DWaveFromStack(wstack,ilayer+1,tmpw)     
			// calc corr with coadd
			if (corrFlag == 2)
				wcoadd/=(ilayer+1)
				pqdrift=autoImageCorrelation(coaddStr,"tmpw",filterFlag)
				// extract layer from stack and coadd
				pp=real(pqdrift)
				qq=imag(pqdrift)
				print "co-add:",stackStr,"layer",ilayer,"drift correction:",pp,qq
				wcoadd*=(ilayer+1)
			else 
				print "co-add:",stackStr,"layer",ilayer
			endif
			addShiftedImage(wcoadd,wstack,ilayer+1,pp,qq)
			ilayer+=1
		while(ilayer<nlayers-1)
		wcoadd/=nlayers
		killwaves tmpw
		String graphNameStr=removeForbiddenChar("COADD_"+stackStr)
		BringImageToFront(graphNameStr,coaddStr)
		Execute "addBWControls()"
	endif
End

Function defectCoadd(): GraphMarquee 
	String ctrlname
	GetMarquee/Z/K left, top
	if (V_Flag == 0)
		Print "There is no marquee"
	else
		// weigh on max or min ???
		variable minMaxFlag=0
		prompt minMaxFlag, "calculate centre of mass of object",popup "min dark;max bright"
		DoPrompt "defect co-add", minMaxFlag
		// marquee position
		variable q0=round(V_top)
		variable q1=round(V_bottom)
		variable p0=round(V_left)
		variable p1=round(V_right)
		// stack dimensions  
		String stackStr=returnListFirstElement(imageNameList("",""),";")
		stackStr= removeQuotes(stackStr)
		wave wstack=$stackStr
		variable np=dimSize(wstack,0)
		variable nq=dimSize(wstack,1)
		variable nlayers=dimSize(wstack,2)
		make/n=(p1-p0,q1-q0)/o tmpw
		// layer 0 to tmpw0
		extractMarqueeFromStack(wstack,0,tmpw,p0,q0)
		// calc centre of object
		variable/C w0_centre= weightImage(tmpw,minMaxFlag)
		variable pp0 = real(w0_centre)
		variable qq0 = imag(w0_centre)
		variable ppi=0
		variable qqi=0
		variable/C wi
		// coadd wave
		String coaddStr=stackStr+"_CO"
		make/n=(np,nq)/o/d $coaddStr
		wave wcoadd=$coaddStr
		addShiftedImage(wcoadd,wstack,0,0,0)
		// coadd loop
		variable ilayer=1
		do
			extractMarqueeFromStack(wstack,ilayer,tmpw,p0,q0)
			wi=weightImage(tmpw,minMaxFlag)
			ppi=round(real(wi)-pp0)
			qqi=round(imag(wi)-qq0)
			print "defect co-add:",stackStr,"layer",ilayer,"drift correction:",ppi,qqi
			addShiftedImage(wcoadd,wstack,ilayer,ppi,qqi)
			ilayer+=1
		while (ilayer < nlayers)
		wcoadd/=nlayers
		killwaves tmpw
		String graphNameStr=removeForbiddenChar("COADD_"+stackStr)
		BringImageToFront(graphNameStr,coaddStr)
		Execute "addBWControls()"
	endif
end

Function/C weightImage(w,Flag)
	wave w
	variable Flag
	ImageStats w
	Redimension/D w
	if (Flag == 1)  // rescaling to [0,1]
		w=(V_min-w)/(V_max-V_min)
	else
		if (flag == 2)
			w=(w-V_min)/(V_max-V_min)
		endif
	endif
	// weigthed sum
	duplicate/o w wpq
	ImageStats w
	variable wtot = V_avg
	wpq=w[p][q]*p
	ImageStats wpq
	variable pp=V_avg/wtot
	duplicate/o w wpq
	wpq=w[p][q]*q
	ImageStats wpq
	variable qq=V_avg/wtot
	killwaves/z wpq
	return cmplx(pp,qq)
end

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// profile functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

function powerSpectrum(ctrlname): buttoncontrol	
	String ctrlname
	String tmpStr=returnListFirstElement(WaveList("*", ";","WIN:") ,";")
	tmpStr= removeQuotes(tmpStr)
	wave w=$tmpStr
	waveStats/Q w
	variable n=numpnts(w)
	n=calcPower2(n)*2
	String fftStr = "FFT_"+tmpStr
	duplicate/o w $fftStr
	wave/c wfft=$fftStr
	wfft[]=cmplx((w[p]-V_avg)*(sin(Pi*p/n))^2,0)
	redimension/n=(n)/d wfft
	FFT wfft
	wfft[] =cmplx(sqrt(real(wfft[p])^2+imag(wfft[p])^2),0)
	redimension/R wfft
	SetScale/I x 0,0.5,"1/pixel", wfft
	String graphStr = "powerFFT_"+tmpStr
	BringToFrontAndDisplay(graphStr,fftStr)
	Execute "addFFTControls()"
end

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// dispersive plane functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Function Dispersive(ctrlname): buttoncontrol	
	String ctrlname
	pauseUpdate
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	// get profile
	variable isok = shoProfile("")
	if (isok)
		DoUpdate;
		Execute "addDispersiveControls()"
		disableDispersiveButtons()
	endif		
end

Function/S stackLayerToDispersiveKE(stackStr,halfwidth,i)
	String stackStr
	variable halfwidth,i
	//  globals
	NVAR lastcursAp = root:lastcursAp
	NVAR lastcursBp = root:lastcursBp
	NVAR lastcursAq = root:lastcursAq
	NVAR lastcursBq = root:lastcursBq
	NVAR layer=root:layer
	NVAR frameno=root:frameno
	NVAR keg=root:keg
	NVAR xzerog=root:xzerog
	NVAR conversiong=root:conversiong
	// cursors
	variable p1,p2,q1,q2
	p1= lastcursAp
	p2 = lastcursBp
	q1= lastcursAq
	q2 = lastcursBq
	setGlobals(stackStr,i)
	// get profile
	if (halfwidth >= 0)
		getProfile(stackStr,halfwidth,p1,p2,q1,q2,i,frameno)
	endif
	String dispersiveStr=stackStr+"_FPR_"+num2ScanNumberStr(frameno)
	wave dispY=$dispersiveStr
	// to  KE
	variable m=numpnts(dispY)
	String dispersiveXStr=dispersiveStr+"_X"
	make/n=(m)/o $dispersiveXStr
	wave wx=$dispersiveXStr
	wx= toKE(p,keg,xzerog,conversiong)
	// removeFromGraph $dispersiveStr
	appendToGraph $dispersiveStr vs wx
	Label left "\\F'Arial'\\Z14intensity (a.u.)";DelayUpdate
	Label bottom "\\F'Arial'\\Z14relative kinetic energy (eV)";DelayUpdate
	SetAxis/A bottom
	ModifyGraph rgb($dispersiveStr)=(abs(enoise(65280)),abs(enoise(65280)),abs(enoise(65280)))
	// cursors
	cursorsOnWave(dispersiveStr,25,numpnts($dispersiveStr)-25)
	// return wavename
	return dispersiveStr
end
	
Macro placeCRS()
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	Cursor/I/P B, $tmpStr, 209,436
	Cursor/I/P A, $tmpStr, 457,135
end

function cursorsOnWave(wStr,p1,p2)
	String wStr
	variable p1,p2
	Cursor/P A, $wStr, p1	// cursor A 
	Cursor/P B, $wStr, p2    // cursor B 
end

function dispersiveScaleCorrection(pix,zero)
	variable pix,zero
	// pix = x axis (in pixels) of disp plane profile from csra=(457,135) to crsB=(209,436)
	// zero = position in pixels of electron beam
	variable diff=(pix-zero)+190
	return -0.51029+(-0.0078623*diff+4.1059e-05*diff^2-7.213e-08*diff^3)
end

function toKE(pix,ke,zero,eV2pix)
	variable pix,ke,zero,eV2pix
	// pix = x axis (in pixels) of disp plane profile from csra=(457,135) to crsB=(209,436)
	// zero = position in pixels of electron beam
	// ke = start voltage of e-beam
	// eV2pix = conversion factor 1ev to pixels
	return ke+(pix-zero)/eV2pix+dispersiveScaleCorrection(pix,zero) 
end

function toBE(pix,ke,kefe,zero,posfe,eV2pix)
	variable pix,ke,kefe,zero,posfe,eV2pix
	// pix = x axis (in pixels) of disp plane profile from csra=(457,135) to crsB=(209,436)
	// zero = position in pixels of electron beam
	// ke = start voltage of e-beam
	// kefe = start voltage of e-beam for fermi edge disp plane image
	// pixel pos of fermi edge
	// eV2pix = conversion factor 1ev to pixels
	return  toKe(posfe,kefe,zero,eV2pix)-toKE(pix,ke,zero,eV2pix)
end

function normCoef(p,zero)
	variable p,zero
	variable ref=190
	variable delta=zero-ref
	return (3035.1+5.9521*(p-delta)-0.082288*(p-delta)^2+0.00034568*(p-delta)^3-4.4797e-07*(p-delta)^4)/2982.62
end

function normDP(tmpStr,zero)
	String tmpStr
	variable zero
	wave w=$tmpStr
	w[]=w[p]/normcoef(p,zero)
end

macro DP_2_KE(tmpStr,ke,zero,eV2pix)
	string tmpStr
	variable ke=83
	variable zero=210
	variable ev2pix=29.6
	Prompt tmpStr, "enter wave",popup WaveList("*",";","WIN:")
	String tmpYStr=tmpStr+"_Y"
	duplicate/o $tmpStr $tmpYStr
	variable offset=$tmpStr[0]
	$tmpYstr-=offset
	normDP(tmpYStr,zero)
	String tmpXStr=tmpStr+"_X"
	duplicate/o $tmpStr $tmpXStr
	variable n=numpnts($tmpXStr)-1
	// to account for wrong integration direction
	// otheriwse, the command would be: $tmpXstr[]=toKE(p,ke,zero,ev2pix)
	$tmpXstr[]=toKE(n-p,ke,zero,ev2pix)
	display $tmpYStr vs $tmpXstr
	Label left "Intensity (a.u.)"
	Label bottom "Kinetic Energy (eV)"
end

Function Shift_EnergyScale(ctrlname): buttoncontrol	
	String ctrlname
	variable ev=0
	String wavenameStr
	Prompt eV, "input energy shift (eV):"
	Prompt wavenameStr,"enter wave:", popup WaveList("*X",";","WIN:")
	doPrompt "Shift binding energy", ev, wavenameStr
	if (V_flag == 0)
		wave wx=$wavenameStr
		wx+=eV
		SetAxis/A
	endif
end

Function ShiftAll_EnergyScale(ctrlname): buttoncontrol	
	String ctrlname
	variable ev=0
	Prompt eV, "input energy shift (eV):"
   	doPrompt "energy shift", ev
	String wavelistStr = WaveList("*X", ";","WIN:") 
	if (V_flag == 0)
		variable n=ItemsInList(wavelistStr) 
		variable m=0
		String wavenameStr = ""
		variable x1,x2
		do
			wavenameStr=StringFromList(m, wavelistStr) 
			wave wx = $wavenameStr
			wx+=eV
			m+=1
		while(m < n)
	setAxis/A
	endif
end

Function nanOutside(ctrlname): buttoncontrol	
	String ctrlname
	String waveStr
	Prompt waveStr, "enter wave",popup WaveList("*_FPR",";","WIN:")
	doPrompt "Shirley background removal" waveStr
	if (stringMatch(num2Str(pcsr(A)),"NaN")==1)
		Cursor/P A, $waveStr,  0
	endif
	if (stringMatch(num2Str(pcsr(B)),"NaN")==1)
		Cursor/P B, $waveStr, (numpnts($waveStr)-1)
	endif
	if (V_flag == 0)
		variable pa=pcsr(A)
		variable pb=pcsr(B)
		variable p1=min(pa,pb)-1
		variable p2=max(pa,pb)+1
		print "setting to NaN ...", waveStr, "outside", p1, p2
		wave wy = $waveStr
		variable np=numpnts(wy)
		wy[0,p1]=NaN
		wy[p2,np-1]=NaN
	endif
end

// in  order to fiind xzerog (if e-gun pos NA): make/n=400/o test= toBE(pixfe1,kefe1,kefe2,p,posfe2,28.63)
Function pixel2Binding(ctrlname): buttoncontrol	
	// ONLY the first FPR wave in the graph is considered
	String ctrlname
	NVAR keg=root:keg
	NVAR kefeg=root:kefeg
	NVAR xzerog=root:xzerog
	NVAR posfeg=root:posfeg
	NVAR conversiong=root:conversiong
	variable ke=keg
	variable kefe=kefeg
	variable xzero=xzerog
	variable posfe=posfeg
	variable conversion=conversiong	
	variable normFlag = 1
	prompt ke, "start voltage of DP image (eV):"
	prompt kefe, "start voltage of Fermi edge image (eV):"
	prompt posfe, "position of fermi edge (pixels);"
	prompt xzero, "postion of e-beam (pixels)"
	prompt conversion, "conversion factor (pixels per eV):"
	prompt normFlag, "normalise", popup "no; yes;"
	doPrompt "converting pixels to eV", ke,kefe,posfe,xzero,conversion,normFlag
	if (V_flag == 0)
		Silent, 1; pauseUpdate
		String wavelistStr = WaveList("*_FPR", ";","WIN:") 
		String wavenameStr=StringFromList(0, wavelistStr) 
		variable m=numpnts($wavenameStr)
		if (normFlag == 2)
			normDP(wavenameStr,m/2)
		endif
		String tmpXStr=wavenameStr+"_X"
		make/n=(m)/o $tmpXStr
		wave wx=$tmpXStr
		wx= toBE(p,ke,kefe,xzero,posfe,conversion)
		print "creating X wave for:",wavenameStr," Start x = ",wx[0]," end x = ",wx[m-1]
		removeFromGraph $wavenameStr
		appendToGraph $wavenameStr vs $tmpXstr
		Label left "\\F'Arial'\\Z16intensity (a.u.)";DelayUpdate
		Label bottom "\\F'Arial'\\Z16binding energy (eV)";DelayUpdate
		SetAxis/A/R bottom
		enableDispersiveOperations()
		variable deltaCrs=25
		cursorsOnWave(wavenameStr,deltaCrs,m-deltaCrs)
		keg=ke
		kefeg=kefe
		xzerog=xzero
		posfeg=posfe
		conversiong=conversion
	endif
end

Function pixel2kinetic(ctrlname): buttoncontrol	
	// ONLY the FPR wave
	String ctrlname
	NVAR keg=root:keg
	NVAR xzerog=root:xzerog
	NVAR conversiong=root:conversiong
	variable ke=keg
	variable xzero=xzerog
	variable conversion=conversiong
	variable normFlag = 1	
	prompt ke, "start voltage of DP image (eV):"
	prompt xzero, "postion of e-beam (pixels)"
	prompt conversion, "conversion factor (pixels per eV):"
	prompt normFlag, "normalise", popup "no; yes;"
	doPrompt "converting pixels to eV", ke,xzero,conversion,normFlag
	if (V_flag == 0)
		Silent, 1; pauseUpdate
		String wavelistStr = WaveList("*_FPR", ";","WIN:") 
		String wavenameStr=StringFromList(0, wavelistStr) 
		variable m=numpnts($wavenameStr)
		if (normFlag == 2)
			normDP(wavenameStr,m/2)
		endif
		String tmpXStr=wavenameStr+"_X"
		make/n=(m)/o $tmpXStr
		wave wx=$tmpXStr
		wx= toKE(p,ke,xzero,conversion)
		print "creating X wave for:",wavenameStr," Start x = ",wx[0]," end x = ",wx[m-1]
		removeFromGraph $wavenameStr
		appendToGraph $waveNameStr vs $tmpXstr
		Label left "\\F'Arial'\\Z16intensity (a.u.)";DelayUpdate
		Label bottom "\\F'Arial'\\Z16kinetic energy (eV)";DelayUpdate
		SetAxis/A bottom
		enableDispersiveOperations()
		variable deltaCrs=25
		cursorsOnWave(wavenameStr,deltaCrs,m-deltaCrs)
		keg=ke
		xzerog=xzero
		conversiong=conversion
	endif
end
			
Function Shirley_equidistant(x,w,n)
	variable x
	wave w
	variable n
	return Area(w,x,n)
end

Function Shirley(x,wx,wy,x2)
	variable x
	wave wx, wy
	variable x2
	return areaXY(wx, wy, x, x2 )
end

function removeShirleyBackGround(ctrlname): buttoncontrol	
	String ctrlname
	String waveStr
	Prompt waveStr, "enter wave",popup WaveList("*_FPR",";","WIN:")
	doPrompt "Shirley background removal" waveStr
	if (V_flag == 0)
		makeShirelyWave(WaveStr)
	endif
end
	
function makeShirelyWave(waveStr)
	String waveStr
	variable v1=NaN,v2=NaN
	if (stringMatch(num2Str(pcsr(A)),"NaN")!=1 && stringMatch(num2Str(pcsr(B)),"NaN")!=1 )
		String shirleyYStr = waveStr+"_SY"
		String shirleyXStr = waveStr+"_SX"
		String shirleyBkgStr = waveStr+"_SB"
		String waveXStr = waveStr+"_X"
		variable pxa=min(pcsr(A),pcsr(B))
		variable pxb=max(pcsr(A),pcsr(B))
		pxb=min(pxb,numpnts($waveStr))
		variable m=abs(pxb-pxa)
		make/n=(m)/o $shirleyXStr,$ShirleyYStr,$shirleyBkgStr
		Silent, 1; pauseUpdate
		wave syw = $shirleyYStr
		wave sxw = $shirleyXStr
		wave w = $waveStr
		wave wx=$waveXStr
		wave wb = $shirleyBkgStr
		// shirely X wave
       		sxw[]=wx[p+pxa]
       		variable x=wx[pxa]								// x=pnt2x($waveStr,pxa)
		variable n=wx[pxb]								// n=pnt2x($waveStr,pxb)
		// shirley Y wave
		syw[] = Shirley(wx[p+pxa],wx,w,wx[pxa+m-1])		// syw[] = Shirley_equidistant(pnt2x(w,p+pxa),w,n)
       		v1=mean(w,pxa,pxa+10) 						// v1=mean(w,pnt2x(w,pxa),pnt2x(w,pxa+10))
       		v2=mean(w,pxb-10,pxb) 						// v2=mean(w,pnt2x(w,pxb-10),pnt2x(w,pxb))
       		variable scalingfactor=(v1-v2)/syw[0]
       		syw*=scalingfactor
       		syw+=v2
       		wb[]=syw[p]
       		syw[] = w[p+pxa]-syw[p]
       		String resultStr=""
       		if (x<n)
       			resultStr=" the peak area is: "+num2Str(AreaXY($shirleyXStr,$shirleyYStr,x,n))+"\r the bkg area is: "+num2Str(AreaXY($shirleyXStr,$shirleybkgStr,x,n))+"\r the bkg avg value is: "+num2Str(faverageXY($shirleyXStr, $shirleybkgStr, n, x ))
       		else
       			resultStr=" the peak area is: "+num2Str(AreaXY($shirleyXStr,$shirleyYStr,n,x))+"\r the bkg area is: "+num2Str(AreaXY($shirleyXStr,$shirleybkgStr,n,x))+"\r the bkg avg value is: "+num2Str(faverageXY($shirleyXStr, $shirleybkgStr, n, x ))
       		endif
       		CheckDisplayed $shirleyYStr
       		if (V_flag==0)
       				appendToGraph $shirleyYStr vs $shirleyXStr
       				appendToGraph $shirleyBkgStr vs $shirleyXStr	
       		endif
       		textBox/C/N=text0/A=RB resultStr
       		return faverageXY($shirleyXStr, $shirleybkgStr, n, x )//v1+v2)/2
       endif
  end  
  
 Function calcAreaShirley(xStr,yStr,left,right)
 	String xStr, yStr
 	variable left,right
 	return AreaXY($xStr,$yStr,min(left,right),max(left,right))
 end
  
 Macro ShirleyBackGround(waveYStr,waveXStr, npntsbkg)
  	String waveYStr
	Prompt waveYStr, "enter Y wave",popup WaveList("*",";","WIN:")
	String waveXStr
	Prompt waveXStr, "enter X wave",popup WaveList("*",";","WIN:")
	variable npntsbkg=1
	Prompt npntsbkg, "number of bkground points"
	if (stringMatch(num2Str(pcsr(A)),"NaN")!=1 && stringMatch(num2Str(pcsr(B)),"NaN")!=1)
		String shirleyXStr=waveXStr+"_SH"
		String shirleyYStr=waveYStr+"_SH"
		String bkgStr=waveYStr+"_BK"
		variable pxa=min(pcsr(A),pcsr(B))
		variable pxb=max(pcsr(A),pcsr(B))
		pxb=min(pxb,numpnts($waveYStr))
		variable m=abs(pxb-pxa)
		make/n=(m)/o $shirleyXStr,$ShirleyYStr,$bkgStr
		Silent, 1; pauseUpdate
		// shirely X wave
       		$shirleyXStr[]=$waveXStr[p+pxa]
       		variable x=$waveXStr[pxa]								
		variable n=$waveXStr[pxb]								
		// shirley Y wave
		$shirleyYStr[] = Shirley($waveXStr[p+pxa],$waveXStr,$waveYStr,$waveXStr[pxa+m-1])		
		variable v1,v2
       		v1=mean($waveYStr,pxa,pxa+npntsbkg) 						
       		v2=mean($waveYStr,pxb-npntsbkg,pxb) 						
       		variable scalingfactor=(v1-v2)/$shirleyYStr[0]
       		$shirleyYStr*=scalingfactor
       		$shirleyYStr+=v2
       		$bkgstr[]=$shirleyYStr[p]
       		$shirleyYStr[] = $waveYStr[p+pxa]-$shirleyYStr[p]
       		String resultStr=""
       		if (x<n)
       			resultStr=" the peak area is: "+num2Str(AreaXY($shirleyXStr,$shirleyYStr,x,n))+" \r the bkg avg value is: "+num2Str(v1) 
       		else
       			resultStr=" the peak area is: "+num2Str(AreaXY($shirleyXStr,$shirleyYStr,n,x))+" \r the bkg avg value is: "+num2Str(v1) 
       		endif
       		print waveYStr, resultStr
       		CheckDisplayed $shirleyYStr
       		if (V_flag==0)
       			appendToGraph $shirleyYStr vs $shirleyXStr
       			appendToGraph $bkgStr vs $shirleyXStr
       		endif
       		textBox/C/N=text0/A=LT resultStr
       endif
  end   
  
  Function fermiEdge_T0_Pix(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) = 
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 6
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 6
	//CurveFitDialog/ w[0] = k1
	//CurveFitDialog/ w[1] = k2
	//CurveFitDialog/ w[2] = k3
	//CurveFitDialog/ w[3] = x0
	//CurveFitDialog/ w[4] = sigma
	//CurveFitDialog/ w[5] = amp
	variable poli=0
	variable x0=w[3]
	variable sigma=w[4]/(2*sqrt(ln(2))*sqrt(2))
	Variable nsteps = 50
	Variable integral = 0
	Variable start = -4*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable y
	Variable normG = 1/sqrt(2*Pi*sigma^2)
	variable k=w[5]
	do
		y=x-start
		poli=0 //w[1]*(y-x0)+w[2]*(y-x0)^2
             	integral+=normG*exp(-(start)^2/(2*sigma^2))*(w[0]+poli/(exp(y-x0)+1)+k/(1+exp((y-x0)/(1/1))))*delta                                  
		start+=delta 	
	while (start <= stop)
	return integral
	
	
end

Function fermiEdge_BE(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) = 
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 6
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 6
	//CurveFitDialog/ w[0] = k1
	//CurveFitDialog/ w[1] = k2
	//CurveFitDialog/ w[2] = k3
	//CurveFitDialog/ w[3] = x0
	//CurveFitDialog/ w[4] = sigma
	//CurveFitDialog/ w[5] = amp
	variable poli=0
	variable x0=w[3]
	variable sigma=w[4]/(2*sqrt(ln(2))*sqrt(2))
	Variable nsteps = 50
	Variable integral = 0
	Variable start = -4*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable y
	Variable normG = 1/sqrt(2*Pi*sigma^2)
	variable k=w[5]
	do
		y=x-start
		poli=w[1]*(-y+x0)+w[2]*(y-x0)^2
             	integral+=normG*exp(-(start)^2/(2*sigma^2))*(w[0]+poli/(exp(-y+x0)+0.001)+k/(1+exp((y-x0)/0.025)))*delta                                  
		start+=delta 	
	while (start <= stop)
	return integral
end

function fit_FermiEdge(ctrlname): buttoncontrol	
	String ctrlname
	String wavenameStr
	Prompt wavenameStr, "enter wave",popup WaveList("*_FPR",";","WIN:")
	doPrompt "target wave" wavenameStr
	if (stringMatch(num2Str(pcsr(A)),"NaN")!=1 && stringMatch(num2Str(pcsr(B)),"NaN")!=1 && V_flag == 0)
		make/n=6/o W_coef=0
		W_coef=0
		W_coef[4]=9
		wavestats/Q/R=(pnt2x($wavenameStr,pcsr(A)),pnt2x($wavenameStr,pcsr(B)) ) $wavenameStr
		W_coef[0]=V_min
		W_coef[5]=(V_max-V_min)*0.5
		duplicate/o $wavenameStr tmp
		differentiate tmp
		wavestats/Q/R=(pnt2x($wavenameStr,pcsr(A)),pnt2x($wavenameStr,pcsr(B)) ) tmp
		W_coef[3]=V_minloc
		killwaves tmp
		FuncFit/H="010000"/Q fermiEdge_T0_Pix W_coef $wavenameStr[pcsr(A),pcsr(B)] /D
		String msgStr = "Fermi edge @ "+num2str(W_coef[3])+"\r"+ "resolution: "+num2Str(W_coef[4]/28.63) 
		TextBox/C/N=text1/A=MC msgStr
	endif
end

function fit_FermiEdge_Pix(ctrlname): buttoncontrol	
	String ctrlname
	String wavenameStr
	Prompt wavenameStr, "enter wave",popup WaveList("*_FPR",";","WIN:")
	doPrompt "target wave" wavenameStr
	if (stringMatch(num2Str(pcsr(A)),"NaN")!=1 && stringMatch(num2Str(pcsr(B)),"NaN")!=1 && V_flag == 0)
		make/n=6/o W_coef=0
		W_coef=0
		W_coef[4]=3
		wavestats/Q/R=(pnt2x($wavenameStr,pcsr(A)),pnt2x($wavenameStr,pcsr(B)) ) $wavenameStr
		W_coef[0]=V_min
		W_coef[5]=(V_max-V_min)*0.5
		duplicate/o $wavenameStr tmp
		differentiate tmp
		wavestats/Q/R=(pnt2x($wavenameStr,pcsr(A)),pnt2x($wavenameStr,pcsr(B)) ) tmp
		W_coef[3]=V_minloc
		killwaves tmp
		FuncFit/H="010000"/Q fermiEdge_T0_Pix W_coef $wavenameStr[pcsr(A),pcsr(B)] /D
		String msgStr = "Fermi edge @ "+num2str(W_coef[3])+"\r"+ "resolution: "+num2Str(W_coef[4]) 
		TextBox/C/N=text1/A=MC msgStr
	endif
end

function fit_FermiEdge_BE()
	String wavenameStr
	Prompt wavenameStr, "enter y wave",popup WaveList("*_FPR",";","WIN:")
	String xwavenameStr
	Prompt xwavenameStr, "enter x wave",popup WaveList("*_X",";","WIN:")
	doPrompt "target wave" wavenameStr, xwavenameStr
	if (stringMatch(num2Str(pcsr(A)),"NaN")!=1 && stringMatch(num2Str(pcsr(B)),"NaN")!=1 && V_flag == 0)
		make/n=6/o W_coef=0
		W_coef=0
		W_coef[4]=0.3
		wavestats/Q/R=(pnt2x($wavenameStr,pcsr(A)),pnt2x($wavenameStr,pcsr(B)) ) $wavenameStr
		W_coef[0]=V_min
		W_coef[5]=(V_max-V_min)*0.5
		duplicate/o $wavenameStr tmp
		differentiate tmp
		wavestats/Q/R=(pnt2x($wavenameStr,pcsr(A)),pnt2x($wavenameStr,pcsr(B)) ) tmp
		W_coef[3]=0
		killwaves tmp
		FuncFit/H="010000"/Q fermiEdge_BE W_coef $wavenameStr[pcsr(A),pcsr(B)]  /X=$xwavenameStr /D
		String msgStr = "Fermi edge @ "+num2str(W_coef[3])+"\r"+ "resolution: "+num2Str(W_coef[4]) 
		TextBox/C/N=text1/A=MC msgStr
	endif
end

function myFitFunction (ctrlname): buttoncontrol	
	String ctrlname
	SVAR cmdStrg=root:cmdStrg
	String cmdStr=cmdStrg
	Prompt cmdStr, "enter command for macro"
	doPrompt "my fit command" cmdStr
	if (V_flag == 0)
		cmdStrg=cmdStr
		Execute cmdStr
	endif
end	

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// wave functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Function loadDatImage(fileName)
	String fileName
	killwaves/z tmpImage0, tmpImage
	String GBstr = "GBLoadWave/T={80,80}/W=1/B=1/O/N=tmpImage/P=imagePath/Q"+" "+"\""+fileName+"\"" 
	Execute GBstr  //  load GB
	wave w0=tmpImage0
	variable np=w0[20] 
	variable nq=w0[21]
	DeletePoints 0,196, w0
	redimension/n=(np,nq) w0
	rename tmpImage0 tmpImage
	flipimage(tmpImage)
end

Function flipImage(w)
	wave w
	duplicate/o w tmp
	variable nq = dimsize(w,1)
	w[][]=tmp[p][nq-1-q]
	killwaves/z tmp
End


Function imageToStack(wstack,wlayer,ilayer)
	wave wstack, wlayer
	variable ilayer
	imagetransform/P=(ilayer)/D=wlayer setplane wstack
End

Function PNGimageToStack(wstack,wlayer,ilayer)
	wave wstack, wlayer
	variable ilayer
	imagetransform/P=(ilayer)/D=wlayer setplane wstack
	//wstack[][][ilayer]=256*wlayer[p][q][1]+wlayer[p][q][2]
End

Function PNGimageTo2DWave(wave2D,wpng)
	wave wave2D, wpng
	wave2D[][]=256*wpng[p][q][1]+wpng[p][q][2]
End	

Function  JPGimageTo2DWave(wave2D,wpng)
	wave wave2D, wpng
	wave2D[][]=wpng[p][q][2]
End		

Function entryToWave(w,entry,i)
	wave w
	variable entry,i
	w[i]=entry
End	

Function MultiplyStackLayer(wstack,coef,ilayer)
	wave wstack
	variable coef,ilayer
	redimension/d wstack
	wstack[][][ilayer]=wstack[p][q][ilayer]*coef
End

Function extract1DWaveFromStack(wstack,p1,q1,wdest)
	wave wstack
	variable p1,q1
	wave wdest
	variable n=dimsize(wstack,2)
	wdest[0,n-1]=wstack[p1][q1][p]
End

Function extract2DWaveFromStack(wstack,i,wdest)
	wave wstack
	variable i
	wave wdest
	wdest[][]=wstack[p][q][i]
End

Function insert2DWaveinStack(wstack,i,win)
	wave wstack
	variable i
	wave win
	wstack[][][i]=win[p][q]
End

Function subtractImageFromStack(wstack,w,layer)
	wave wstack
	wave w
	variable layer
	wstack[][][layer]-=w[p][q]
end

Function darkFieldStack(wstack,layer)
	wave wstack
	variable layer
	variable np=dimsize(wstack,0)
	variable nq=dimsize(wstack,1)
	make/n=(np,nq)/o/b/u tmpROI
	make/n=(np,nq)/o/w tmp
	tmpROI=1
	variable rp=round(np/8)
	variable rq=round(np/8)
	tmpROI[0,rp][0,rq]=0
	tmpROI[np-1-rp,np-1][0,rq]=0
	tmpROI[0,rp][nq-1-rq,nq-1]=0
	tmpROI[np-1-rp,np-1][nq-1-rq,nq-1]=0
	extract2DWaveFromStack(wstack,layer,tmp)
	ImageStats/R=tmpROI/M=1 tmp
	tmp-=(V_min) //-(V_avg-V_min)) //avoids overflow //debug print V_avg, V_min
	insert2DWaveinStack(wstack,layer,tmp)
	killwaves/z tmp,tmpROI 
end

Function calcDarkfield(wStack,i)
	wave wstack
	variable i
	variable darkfield
	make/n=(dimSize(wstack,0),dimSize(wstack,1)), tmp
	extract2DWaveFromStack(wstack,i,tmp)
	duplicate/o tmp mymask
	mymask=1
	redimension/B/U mymask
	mymask[200,250][200,250]=0 
	ImageStats/R=mymask tmp
	darkfield=V_avg
	killwaves mymask, tmp
	return darkfield
end

Function returnDarkfield(wStack,layer)
	wave wstack
	variable layer
	variable np=dimsize(wstack,0)
	variable nq=dimsize(wstack,1)
	make/n=(np,nq)/o/b/u tmpROI
	make/n=(np,nq)/o tmp
	tmpROI=1
	variable rp=round(np/8)
	variable rq=round(np/8)
	tmpROI[0,rp][0,rq]=0
	tmpROI[np-1-rp,np-1][0,rq]=0
	tmpROI[0,rp][nq-1-rq,nq-1]=0
	tmpROI[np-1-rp,np-1][nq-1-rq,nq-1]=0
	extract2DWaveFromStack(wstack,layer,tmp)
	ImageStats/R=tmpROI/M=1 tmp
	killwaves/z tmp,tmpROI 
	return V_avg
end

Function addImages(w1,w2)
	wave w1,w2
	return w1+w2
End	


Function extractMarqueeFromStack(wstack,ilayer,wdest,p0,q0)
	wave wstack
	variable ilayer
	wave wdest
	variable p0,q0
	wdest[][]=wstack[p+p0][q+q0][ilayer]
end

Function returnWaveElement(w,i)
	wave w
	variable i
	return w[i]
End	

Function normaliseWave(w,wnorm)
	wave w
	wave wnorm
	w/=wnorm
End

Function divideWave(w,d)
	wave w
	variable d
	w/=d
End

Function subtractOffset(w,offset)
	wave w
	variable offset
	w-=offset
End

Function backGroundCorrection(win,wbkg,coef,rmsbkg)
	wave win
	wave wbkg
	variable coef,rmsbkg
	redimension/D win
	win=win/(wbkg+10)*rmsbkg
end

Function calcXMCD(wstack,wxmcd,signxmcd)
	wave wstack
	wave wxmcd
	variable signxmcd
	if (signxmcd<0)
		wxmcd[][][0]=(wstack[p][q][0]-wstack[p][q][1])/(wstack[p][q][0]+wstack[p][q][1])
	else
		wxmcd[][][0]=(-wstack[p][q][0]+wstack[p][q][1])/(wstack[p][q][0]+wstack[p][q][1])
	endif
End	

Function calcXMCDdrift(wstack,wxmcd,signxmcd,pdrift,qdrift)
	wave wstack
	wave wxmcd
	variable signxmcd
	variable pdrift,qdrift
	if (signxmcd<0)
		wxmcd[][][0]=(wstack[p][q][0]-wstack[p+pdrift][q+qdrift][1])/(wstack[p][q][0]+wstack[p+pdrift][q+qdrift][1])
	else
		wxmcd[][][0]=(-wstack[p][q][0]+wstack[p+pdrift][q+qdrift][1])/(wstack[p][q][0]+wstack[p+pdrift][q+qdrift][1])
	endif
End	

Function addShiftedImage(wtarget,wstack,layer,pdrift,qdrift)
	wave wtarget
	wave wstack
	variable layer,pdrift,qdrift
	wtarget[][]+=wstack[p+pdrift][q+qdrift][layer]
End	

Function ShiftStackLayer(wstack,ilayer,pdrift,qdrift)
	wave wstack
	variable ilayer
	variable pdrift,qdrift
	variable np=dimsize(wstack,0)
	variable nq=dimsize(wstack,1)
	make/n=(np,nq)/o tmpw
	tmpw=wstack[p][q][ilayer]
	wstack[][][ilayer]=tmpw[p+pdrift][q+qdrift]
End	

function writeCol(targetStr,inStr,colno)
	String targetStr,inStr
	variable colno
	wave wtarget=$targetStr
	wave win=$inStr
	wtarget[][colno]=win[p]
end

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// plot functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Function stackToFrontAndDisplay(graphNameStr,sourceNameStr,waveNameStr)
	String graphNameStr
	String sourceNameStr
	String waveNameStr
	String xlabel=""
	String testStr=sourceNameStr+"_NUM"
	String xwaveStr=""
	if (exists(testStr)==1)
		xlabel="\\F'Arial'\\Z14Frame number"
		xwaveStr=testStr
	endif
	testStr=sourceNameStr+"_STV"	
	if (exists(testStr)==1)
		xlabel="\\F'Arial'\\Z14Kinetic energy (eV)"
		xwaveStr=testStr
	endif
	testStr=sourceNameStr+"_NRG"
	if (exists(testStr)==1)
		xlabel="\\F'Arial'\\Z14Photon energy (eV)"
		xwaveStr=testStr
	endif
	testStr=sourceNameStr+"_OBJ"
	if (exists(testStr)==1)
		xlabel="\\F'Arial'\\Z14Photon energy (eV)"
		xwaveStr=testStr
	endif
	Dowindow $graphNameStr
	if (V_flag==0)
		Display
		DoWindow/C $graphNameStr
		appendToGraph $waveNameStr vs $xwaveStr
		Label bottom xlabel
		ModifyGraph mirror=1
		Label left "\\F'Arial'\\Z14Intensity (A.U.)";DelayUpdate
		Legend/C/N=text0/A=LT
		ModifyGraph rgb($waveNameStr)=(65280,16384,16384)
	else
		DoWindow/F $graphNameStr
		String empty=WaveList(waveNameStr, ";", "WIN:"+graphNameStr)
		if (strlen(empty)==0)
			appendtoGraph $waveNameStr vs $xwaveStr
			ModifyGraph rgb($waveNameStr)=(abs(enoise(65280)),abs(enoise(65280)),abs(enoise(65280)))
		endif
	endif
End

Function bringToFront(graphNameStr)
	String graphNameStr
	Dowindow $graphNameStr
	if (V_flag==0)
	else
		DoWindow/F $graphNameStr
	endif
	return V_flag
End

Function bringImageToFront(graphNameStr,waveNameStr)
	String graphNameStr,waveNameStr
	Dowindow $graphNameStr
	if (V_flag==0)
		Display
		DoWindow/C $graphNameStr
		appendImage/L/T $waveNameStr
		SetAxis/A/R left
		ModifyGraph nticks=30,minor=0,sep=1
		Execute "AutoSizeImage(1,1)"
	else
		DoWindow/F $graphNameStr
	endif
	return V_flag
End
 
Function bringToFrontAndDisplay(graphNameStr,waveNameStr)
	String graphNameStr
	String waveNameStr
	graphNameStr= removeForbiddenChar(graphNameStr)
	Dowindow $graphNameStr
	if (V_flag==0)
		Display
		DoWindow/C $graphNameStr
		appendToGraph $waveNameStr
		Legend/C/N=text0/A=LT
		ModifyGraph rgb($waveNameStr)=(abs(enoise(65280)),abs(enoise(65280)),abs(enoise(65280)))
		ModifyGraph mirror=1
	else
		DoWindow/F $graphNameStr
		String empty=WaveList(waveNameStr, ";", "WIN:"+graphNameStr)
		if (strlen(empty)==0)
			appendToGraph $waveNameStr
			ModifyGraph rgb($waveNameStr)=(abs(enoise(65280)),abs(enoise(65280)),abs(enoise(65280)))
		endif
	endif
End
 
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// all string functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Function/S shortenWaveName(str)  // removes file type from wavename - truncates at 20 char
	String str
	variable i
	if (strsearch(str, ".tiff",0)>=0)
		str=str[0,strsearch(str, ".tiff",0)-1]
	endif
	if (strsearch(str, ".tif",0)>=0)
		str=str[0,strsearch(str, ".tif",0)-1]
	endif
	if (strsearch(str, ".t16",0)>=0)
		str=str[0,strsearch(str, ".t16",0)-1]
	endif
	if (strsearch(str, ".T16",0)>=0)
		str=str[0,strsearch(str, ".T16",0)-1]
	endif
	if (strsearch(str, ".PNG",0)>=0)
		str=str[0,strsearch(str, ".PNG",0)-1]
	endif
	if (strsearch(str, ".png",0)>=0)
		str=str[0,strsearch(str, ".png",0)-1]
	endif
	if (strsearch(str, ".JPG",0)>=0)
		str=str[0,strsearch(str, ".JPG",0)-1]
	endif
	if (strsearch(str, ".jpeg",0)>=0)
		str=str[0,strsearch(str, ".jpeg",0)-1]
	endif
	if (strsearch(str, ".jpg",0)>=0)
		str=str[0,strsearch(str, ".jpg",0)-1]
	endif
	if (strsearch(str, ".dat",0)>=0)
		str=str[0,strsearch(str, ".dat",0)-1]
	endif
	if (strsearch(str, ".DAT",0)>=0)
		str=str[0,strsearch(str, ".DAT",0)-1]
	endif
	if (strlen(str)>20)
		str=str[0,20]
	endif
	return str
end

Function/S removeQuotes(str)
	String str
	// debug print "folderRootNameFrom Path input:", str
	// get rid of last ' 
	variable n=strlen(str)
	if (stringMatch(str[n-1],"'")==1)
		str=str[0,n-2]
	endif
	// get rid of first ' 
	n=strlen(str)
	if (stringMatch(str[0],"'")==1)
		str=str[1, n-1]
	endif
	// debug print "folderRootNameFrom Path output:", str
	return str
end

Function/S removeForbiddenChar(str)
	String str 
	variable n=strlen(str)
	variable i
	for(i=0;i<n;i+=1)
		if(stringMatch(str[i],"#")==1)
			str=str[0,i-1]+"_"+str[i+1,n-1]
		endif
		if(stringMatch(str[i],"-")==1)
			str=str[0,i-1]+"_"+str[i+1,n-1]
		endif
	endfor
	return str
end
 
Function/S returnFolderRootNameFromPath(str)
	String str
	// debug print "folderRootNameFrom Path input:", str
	//variable lastindex = strlen(str)-strsearch(str,"_",3)
	//str=str[0,lastindex-3]
	// get rid of last ' 
	variable n=strlen(str)
	if (stringMatch(str[n-1],"'")==1)
		str=str[0,n-2]
	endif
	// get rid of first ' 
	n=strlen(str)
	if (stringMatch(str[0],"'")==1)
		str=str[1, n-1]
	endif
	// get rid of last :
	n=strlen(str)
	if (stringMatch(str[n-1],":")==1)
		str=str[0,n-2]
	endif
	// if *-NNN.* get rid of it
	n=strlen(str)
	variable j=n-1
	do
		j-=1;			
	while(j>=0 && stringMatch(str[j],".")==0)
	if (j>0 && j-5>=0 && stringMatch(str[j-4],"-")==1 && isNumber(str[j-3,j-1])==1)
		str=str[0,j-5]
	endif
	// extract folder name
	n=strlen(str)
	j=n-1
	do
		j-=1;			
	while(j>=0 && stringMatch(str[j],":")==0)
	str=str[j+1,n-1]
	// debug print "folderRootNameFrom Path output:", str
	return str
end

Function isNumber(str)
	String str
	variable n = strlen(str)
	do
		n-=1
	while (0<=n && stringMatch(num2Str(str2num(str[n])),"NaN")==0)
	if (n>=0)
		return 0
	else
		if (stringMatch(str,"")==1)
			return 0 
		else
			return 1
		endif
	endif
end

Function/S returnFolderNumberFromPath(str)
	String str	
      // debug print "returnFolderNumberFromPath Path input:", str
	variable lastindex = strlen(str)-strsearch(str,"_",3)
	str=str[0,lastindex-3]
	print "string", str
	// get rid of last ' 
	variable n=strlen(str)
	if (stringMatch(str[n-1],"'")==1)
		str=str[0,n-2]
	endif
	// get rid of first ' 
	n=strlen(str)
	if (stringMatch(str[0],"'")==1)
		str=str[1, n-1]
	endif
	// get rid of last :
	n=strlen(str)
	if (stringMatch(str[n-1],":")==1)
		str=str[0,n-2]
	endif	
	// if *-NNN.* get rid of it
	n=strlen(str)
	variable j=n-1
	do
		j-=1;			
	while(j>=0 && stringMatch(str[j],".")==0)
	if (j>0 && j-5>=0 && stringMatch(str[j-4],"-")==1 && isNumber(str[j-3,j-1])==1)
		str=str[0,j-5]
	endif
	// look for last *_NNN and checks whether is number
	n=strlen(str)
	j=n-1
	do
		j-=1;			
	while(j>=0 && stringMatch(str[j],"_")==0)
	str=str[j+1,n-1]
	if (isNumber(str)==0)
		str=""
	endif
	// debug print "returnFolderNumberFromPath Path output:", str
	return str
end

Function/S returnFileNumberFromFileName(str)
	String str
	// if *-NNN.* get it
	variable n=strlen(str)
	variable j=n-1
	do
		j-=1;			
	while(j>=0 && stringMatch(str[j],".")==0)
	if (j>0 && j-4>=0 && stringMatch(str[j-4],"-")==1 && isNumber(str[j-3,j-1])==1)
		str=str[j-3,j-1]
	else
		str=""
	endif
	return str
end

Function/S returnPathFromFullPath(str)
	String str
	variable n=strlen(str)
	variable j=n-1
	do
		j-=1;			
	while(j>=0 && stringMatch(str[j],":")==0)
	return str[0,j]
end

Function/S returnFileNameFromFullPath(str)
	String str
	variable n=strlen(str)
	variable j=n-1
	do
		j-=1;			
	while(j>=0 && stringMatch(str[j],":")==0)
	return str[j+1,n-1]
end

Function/S returnListFirstElement(str,separator)
	String str
	String separator
	variable n=strlen(str)
	variable j=0
	do
		j+=1;			
	while(j<n && stringMatch(str[j],separator)==0)
	str=str[0,j-1]
	return str
end

Function/S num2ScanNumberStr(p)
	variable p
	if (p<10) 
		return "00"+num2Str(p)
	else
		if (p<100)
			return "0"+num2Str(p)
		else
			return num2Str(p)
		endif
	endif
end
 
Function/S returnImageName(str)
	String str
	variable n=strlen(str)
	variable j=n-1
	do
		j-=1;			
	while(j>=0 && stringMatch(str[j],":")==0)
	str=str[j+1,n-1]
	return str
end

Function/S returnImageType(str)
	String str
	variable n=strlen(str)
	variable j=n-1
	do
		j-=1;			
	while(j>=0 && stringMatch(str[j],".")==0)
	str=str[j+1,n-1]
	return str
end

Function/S getBcsParameter(inputStr,paramStr)
	String inputStr,paramStr
	variable is_tcp=strsearch(inputStr,"-----------",0)
	variable place = 2
	if (is_tcp>1) 
		place =3
	endif	
	variable i1=strsearch(inputStr,paramStr,0)
	String str=inputStr[i1,strlen(inputStr)]
	variable i
	for(i=0;i<place;i+=1)
		i1=strsearch(str,",",0)	
		str=str[i1+1,strlen(str)]	
	endfor
	i1=strsearch(str,",",0)
	str=str[0,i1-1]
	return str	
End 

Function/S getLEEMParameter(inputStr,paramStr)
	String inputStr,paramStr
	variable i1=strsearch(inputStr,paramStr,0)
	String str=inputStr[i1,strlen(inputStr)]
	variable i
	i1=strsearch(str,":",0)	
	str=str[i1+1,strlen(str)]	
	i1=strsearch(str,"\n",0)
	str=str[0,i1-1]
	return str	
End 

Function/S returnSequenceNumber(i,len)
	variable i,len
	variable ilen=strlen(num2Str(i))
	if (ilen<=len)
		if (ilen<len)
			variable j=len-ilen
			String out=""
			variable k=0
			do
				out+="0"
				k+=1
			while (k<j)
			return out+ num2str(i)
		else
			return num2str(i)
		endif
	else
		return ""
	endif
End

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// layout procedures
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Proc addBWcontrols()
	       //add controls
		PauseUpdate
		showinfo
		ControlBar 101
		//first line
		button resizeButton, proc=resizeImageWindow, pos={1,1} , size={78,18}, title="resize"	
		button saveFrameButton, proc=saveJpeg ,pos={81,1} , size={78,18}, title="save to jpeg"
		button saveTiffFrameButton, proc=saveTiff ,pos={161,1} , size={78,18}, title="save to tiff"
		button killButton, proc=killSeries ,pos={241,1} , size={78,18}, title="kill"	
		// second line
		button rescaleButton, proc=rescaleScan ,pos={1,21} , size={78,18}, title="rescale"
		button invertButton, proc=InvertImage ,pos={81,21} , size={78,18}, title="invert"
		button despikeButton, proc=despikeScan ,pos={161,21} , size={78,18}, title="threshold"
		button normButton, proc=normalise ,pos={241,21} , size={78,18}, title="normalise"  
		//third line
		button showRoi, proc=shoRoi ,pos={1,41} , size={78,18}, title="ROI tools" 
		button bkgButton, proc=bkgRemove ,pos={81,41} , size={78,18}, title="remove bkg" 
		button layoutButton, proc=Imagelayout ,pos={161,41} , size={78,18}, title="layout" 
		button profileButton, proc=shoProfile ,pos={241,41} , size={78,18}, title="profile"   
		//fourth line
		button dispersiveButton, proc=dispersive ,pos={1,61} , size={78,18}, title="dispersive"
		button frameHistogramButton, proc=shoFramehistogram ,pos={81,61} , size={78,18}, title="histo"
		button fftButton, proc=imageFFT ,pos={161,61} , size={78,18}, title="FFT"
		button prevSubset, proc=prevSubset, pos={241,61} , size={78,18}, title="prev subset"	
		//fifth line
		button redimensionButton, proc=toUB, pos={1,81} , size={78,18}, title=" red /U/B"	
		button paButton, proc=doPA ,pos={81,81} , size={78,18}, title="particle"
		button binningButton, proc=binning_2x2, pos={161,81} , size={78,18}, title="bin 2x2"
		
end

Proc addScanControls()  //160203
	       //add controls
		PauseUpdate
		showinfo
		ControlBar 121
		//first line
		button resizeButton, proc=resizeImageWindow, pos={1,1} , size={78,18}, title="resize"	
		SetVariable setLayer, proc=shoLayer, pos={81,1} , size={78,18}, frame=1, font="Helvetica",title="layer" , value=layer
		ValDisplay keValDisplay, pos={161,1} , size={78,18}, frame=1, font="Helvetica",title="eV:" , value=energy
		button movieButton, proc=shoMovie ,pos={241,1} , size={78,18}, title="scan movie"
		// second line
		button saveFrameButton, proc=saveJpeg ,pos={1,21} , size={78,18}, title="layer to jpeg"
		button saveTiffFrameButton, proc=saveTiff ,pos={81,21} , size={78,18}, title="layer to tiff" 
		button saveScanButton, proc=scanToJpeg ,pos={161,21} , size={78,18}, title="scan to jpeg"   
		button killButton, proc=killSeries ,pos={241,21} , size={78,18}, title="kill scan"	
		// third line
		button showRoi, proc=shoRoi ,pos={1,41} , size={78,18}, title="ROI tools" 
		button roiIntegrateButton, proc=shoRoiScan ,pos={81,41} , size={78,18}, title="ROI scan"
		button pixelIntegrateButton, proc=shoPixelScan ,pos={161,41} , size={78,18}, title="pixel scan"  
		button profileButton, proc=shoProfile ,pos={241,41} , size={78,18}, title="profile" 
		//forth line
		button frameHistogramButton, proc=shoFramehistogram ,pos={1,61} , size={78,18}, title="layer histo"
		button despikeButton, proc=despikeScan ,pos={81,61} , size={78,18}, title="threshold"  
		button rescaleButton, proc=rescaleScan ,pos={161,61} , size={78,18}, title="rescale"
		button invertButton, proc=invertImage ,pos={241,61} , size={78,18}, title="invert"
		// fifth line
		button correlarionButton, proc=smartCoadd ,pos={1,81} , size={78,18}, title="co-add"	
		button driftButton, proc=shift2DriftCorrect ,pos={81,81} , size={78,18}, title="drift corr"	
		button manualDriftButton, proc=manualDrift ,pos={161,81} , size={78,18}, title="drift layer"
		button makeMovieButton, proc=makeMovie ,pos={241,81} , size={78,18}, title="make movie"   	
		// sixth line
		button userButton1, proc=userMacro1 ,pos={1,101} , size={78,18}, title="usermacro 1"	
		button userButton2, proc=userMacro2 ,pos={81,101} , size={78,18}, title="usermacro 2"
		button normaliseButton, proc=normaliseSequence, pos={161,101} , size={78,18}, title="normalise"
		button binningButton, proc=binning_2x2, pos={241,101} , size={78,18}, title="bin 2x2"
end

Proc addSequenceControls()  //160203
                 //add controls
		PauseUpdate
		showinfo
		ControlBar 121
		// first line
		button resizeButton, proc=resizeImageWindow, pos={1,1} , size={78,18}, title="resize"	
		SetVariable shoLayerVar, proc=shoLayer, pos={81,1} , size={78,18}, frame=1, font="Helvetica",title="layer" , value=layer
		ValDisplay frameValDisplay, pos={161,1} , size={78,18}, frame=1, font="Helvetica",title="#" ,value=frameno 
 		button movieButton, proc=shoMovie ,pos={241,1} , size={78,18}, title="scan movie"
 		// second line
		button saveFrameButton, proc=saveJpeg ,pos={1,21} , size={78,18}, title="layer to jpeg"
		button saveTiffFrameButton, proc=saveTiff ,pos={81,21} , size={78,18}, title="frame to tiff" 
		button saveScanButton, proc=scanToJpeg ,pos={161,21} , size={78,18}, title="scan to jpeg"   
		button killButton, proc=killSeries ,pos={241,21} , size={78,18}, title="kill"	
		//third line
		button showRoi, proc=shoRoi ,pos={1,41} , size={78,18}, title="ROI tools" 
		button roiIntegrateButton, proc=shoRoiScan ,pos={81,41} , size={78,18}, title="ROI scan"
		button pixelIntegrateButton, proc=shoPixelScan ,pos={161,41} , size={78,18}, title="pixel scan"  
		button profileButton, proc=shoProfile ,pos={241,41} , size={78,18}, title="profile"
		//forth line
		button frameHistogramButton, proc=shoFramehistogram ,pos={1,61} , size={78,18}, title="layer histo"
		button despikeButton, proc=despikeScan ,pos={81,61} , size={78,18}, title="threshold"  
		button rescaleButton, proc=rescaleScan ,pos={161,61} , size={78,18}, title="rescale"
		button invertButton, proc=invertImage ,pos={241,61} , size={78,18}, title="invert"
		//fifth line	
		button correlarionButton, proc=smartCoadd ,pos={1,81} , size={78,18}, title="co-add"
		button driftButton, proc=shift2DriftCorrect ,pos={81,81} , size={78,18}, title="drift corr"		
		button manualDriftButton, proc=manualDrift ,pos={161,81} , size={78,18}, title="drift layer"	
		button makeMovieButton, proc=makeMovie ,pos={241,81} , size={78,18}, title="make movie"   	
		// sixth line
		button userButton1, proc=userMacro1 ,pos={1,101} , size={78,18}, title="usermacro 1"
		button userButton2, proc=userMacro2 ,pos={81,101} , size={78,18}, title="usermacro 2"
		button normaliseButton, proc=normaliseSequence, pos={161,101} , size={78,18}, title="normalise"
		button binningButton, proc=binning_2x2, pos={241,101} , size={78,18}, title="bin 2x2"
end

Proc addXMCDControls()  //160203	
		//add controls
		PauseUpdate
		showinfo
		ControlBar 81
		// first line
		button resizeButton, proc=resizeImageWindow, pos={1,1} , size={78,18}, title="resize"
		button saveFrameButton, proc=save_XMCD_Jpeg ,pos={81,1} , size={78,18}, title="frame to jpeg" 
		button saveTiffFrameButton, proc=saveTiff ,pos={161,1} , size={78,18}, title="frame to tiff" 		
		button killButton, proc=killSeries ,pos={241,1} , size={78,18}, title="kill"	
		// second line 
		button showRoi, proc=shoRoi ,pos={1,21} , size={78,18}, title="ROI tools" 	
		button bkgButton, proc=bkgRemove ,pos={81,21} , size={78,18}, title="remove bkg" 
		button layoutButton, proc=XMCDlayout ,pos={161,21} , size={78,18}, title="XMCDlayout"
		button prevSubset, proc=prevSubset, pos={241,21} , size={78,18}, title="prev subset" 
		// third
		button despikeButton, proc=despikeXMCD ,pos={1,41} , size={78,18}, title="threshold" 
		button rescaleButton, proc=rescaleScan ,pos={81,41} , size={78,18}, title="rescale"
		button invertButton, proc=InvertImage ,pos={161,41} , size={78,18}, title="invert"
		button makeAllFrameButton, proc=makeAll ,pos={241,41} , size={78,18}, title="all" 
		// fourth line
		button profileButton, proc=shoProfile ,pos={1,61} , size={78,18}, title="profile"   
		button frameHistogramButton, proc=shoFramehistogram ,pos={81,61} , size={78,18}, title="histo"		
End

Proc addRawDataControls()
	        //add controls
		PauseUpdate
		showinfo
		ControlBar 61
		//first line
		button loopButton, proc=shoLoop ,pos={1,1} , size={78,18}, title="loop"
		SetVariable setLayer, proc=shoLayer, pos={81,1} , size={78,18}, frame=1, font="Helvetica",title="layer" , value=layer
		button saveFrameButton, proc=saveJpeg ,pos={161,1} , size={78,18}, title="layer to jpeg"  
		button killButton, proc=killSeries ,pos={241,1} , size={78,18}, title="kill"	
		// second line
		button despikeButton, proc=despikeScan ,pos={1,21} , size={78,18}, title="threshold"  
		button rescaleButton, proc=rescaleScan ,pos={81,21} , size={78,18}, title="rescale"
		button driftButton, proc=driftCorrection ,pos={161,21} , size={78,18}, title="drift correct"
		button rexmcdButton, proc=rexmcd ,pos={241,21} , size={78,18}, title="re-xmcd" 	
		button manualDriftButton, proc=manualDrift ,pos={321,21} , size={78,18}, title="drift layer"
		// third line
		button resizeButton, proc=resizeImageWindow, pos={1,41} , size={78,18}, title="resize"
		button profileButton, proc=shoProfile ,pos={81,41} , size={78,18}, title="profile"   
		button frameHistogramButton, proc=shoFramehistogram ,pos={161,41} , size={78,18}, title="histo"	
		button prevSubset, proc=prevSubset, pos={241,41} , size={78,18}, title="prev subset"	
		button normaliseButton, proc=normaliseSequence, pos={321,41} , size={78,18}, title="normalise"	
End

Proc addDispersiveControls()
	ControlBar 41
	button shiftButton, proc=Shift_EnergyScale ,pos={1,1} , size={79,18}, title="Shift wave"  
	button shiftAllButton, proc=ShiftAll_EnergyScale ,pos={81,1} , size={79,18}, title="Shift All"  
	button pixel2bindingButton, proc=pixel2Binding ,pos={161,1} , size={79,18}, title="2 binding" 
	button pixel2kineticButton, proc=pixel2kinetic ,pos={241,1} , size={79,18}, title="2 kinetic" 
	button shirleyButton, proc=removeShirleyBackGround ,pos={321,1} , size={79,18}, title="Shirley"
	button NaNButton, proc=nanOutside ,pos={401,1} , size={78,18}, title="NaN outside"	
	button killButton, proc=killSeries ,pos={481,1} , size={78,18}, title="kill"	
	button fermiButton, proc=fit_FermiEdge ,pos={1,21} , size={79,18}, title="Fermi Edge" 
	button myfitButton, proc=myFitFunction, pos={81,21} , size={79,18}, title="fit" 
End

Proc addProfileControls()
	ControlBar 21
	button fftButton, proc=powerSpectrum ,pos={561,1} , size={79,18}, title="FFT"
	button killButton, proc=killSeries ,pos={481,1} , size={78,18}, title="kill"	
End

Proc addPixelScanControls()
	ControlBar 21
	button averageButton, proc=PQaverage ,pos={1,1} , size={178,18}, title="pixel average"  
	button killButton, proc=killSeries ,pos={401,1} , size={78,18}, title="kill"	
End

Proc addRoiScanControls()
	ControlBar 21 
	button killButton, proc=killSeries ,pos={401,1} , size={78,18}, title="kill"	
End

Proc addHistoScanControls()
	ControlBar 21 
	button fitHistoButton, proc=fit_Histogram,pos={321,1} , size={78,18}, title="fit"	
	button killButton, proc=killSeries ,pos={401,1} , size={78,18}, title="kill"	
End

Proc addFFTControls()
	ControlBar 21
	button killButton, proc=killSeries ,pos={401,1} , size={78,18}, title="kill"	
end

function disableDispersiveButtons()
	button shiftButton disable=2
	button shiftAllButton disable=2
	button shirleyButton disable=2
	button fermiButton disable=0
	button fftButton disable=1
end

function enableDispersiveOperations()
	button shiftButton disable=0
	button shiftAllButton disable=0
	button shirleyButton disable=0
	button fermiButton disable=2
end

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// macro menu
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// t16 XMCD
macro XMCD()
	 manual_XMCD_16bit()
end

function manual_XMCD_16bit()
	//general variables
	SVAR imagePathStr=root:imagePathStr
	SVAR imagetypeStr=root:imageTypeStr
	SVAR scanNameStr=root:scanNameStr
	SVAR scanNumberStr=root:scanNumberStr
	SVAR lastImage=root:lastImage
	NVAR scanFrom=root:scanFrom
	NVAR scanTo=root:scanTo
	NVAR energy=root:energy
	NVAR layer=root:layer
	NVAR helrescaling=root:helrescaling
	//other variables
	string tmpStr, tmpStr1, stackStr
	variable j=0
	// darkfield
    	variable darkfield=2
       	prompt darkfield, "subtact darkfield", popup, "no;yes"
       	doPrompt "open image", darkfield	
        if ((V_flag ==0 )&& (darkfield == 2))
		darkfield=1
	else
		darkfield=0
	endif
	// requires input from the user
	Variable refNum
	String fileName="image"
	Open/D/P=imagePath/R/M="Open helicity minus"/T="????" refNum as fileName
	print "will open image (hel minus): ",returnImageName(S_fileName), "of type: ", returnImageType(S_fileName)
	if (stringMatch(S_fileName,"")!=1)
		// initialise the stack   	
		tmpStr = OpenImage(S_filename,darkfield)
		variable np=DimSize($tmpStr,0)
		variable nq=DimSize($tmpStr,1)    
      		print "the XMCD couple will be named: ",tmpStr
      		// first image is always negative and therefore put in layer 0
      		stackStr=tmpStr+"MP"
      		make/n=(np,nq,2)/o/u $stackStr
      		imageToStack($stackStr,$tmpStr,0)
      		killwaves/z $tmpStr
      		// open plus helicity
		Open/D/P=imagePath/R/M="Open helicity plus"/T="????" refNum as fileName
		print "will open image (hel plus): ",returnImageName(S_fileName), "of type: ", returnImageType(S_fileName)
		if (stringMatch(S_fileName,"")!=1)
			tmpStr1 = OpenImage(S_filename,darkfield)
			imageToStack($stackStr,$tmpStr1,1)
			killwaves/z $tmpStr1
			// cosmetics
			NewImage/S=2 $stackStr
			ModifyImage $stackStr plane=0
			layer=1
			doUpdate
			// commands
			if ((stringMatch(num2Str(np),"NaN")!=1)&&(stringMatch(num2Str(nq),"NaN")!=1))	
				Execute "addRawDataControls()"
				ModifyGraph width=250*np/nq,height=250
			endif
		// rescaling for plus helcity
       		variable coef = helRescaling
       		prompt coef, "rescaling factor for helicity minus"
       		doPrompt "XMCD input", coef
       		if (V_flag==0)
       			MultiplyStackLayer($stackStr,coef,0)
       			helRescaling = coef
       		endif
			// xmcd subtraction
			String xmcdStr=stackStr +"_"+"XMCD"
			make/n=(np,nq,1)/o $xmcdStr
			redimension/D $xmcdStr,$stackStr
			calcXMCD($stackStr,$xmcdStr,-1)
			// display
 			NewImage/S=2 $xmcdStr
			ModifyImage $xmcdStr plane=0
			layer=0
			if ((stringMatch(num2Str(np),"NaN")!=1)&&(stringMatch(num2Str(nq),"NaN")!=1)) 
				//add controls
				Execute "addXMCDControls()"
				ModifyGraph width=250*np/nq,height=250
				TextBox/C/N=text0/A=LT xmcdStr
				DoUpdate
				ModifyGraph width=0,height=0
				ModifyImage $xmcdstr ctab= {-0.3,0.3,Grays,0}
			endif			
		endif
	endif
end

// channelplatemask
Macro make_MCP_ROI()
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	String ctrlname
	tmpStr= removeQuotes(tmpStr)
	variable pp=dimsize($tmpStr,0)
	variable qq=dimsize($tmpStr,1)
	SetDrawEnv arrow=1,linefgc= (65280,21760,0),xcoord= top,ycoord= left;DelayUpdate
	SetDrawEnv fillpat= 0
	variable x0,y0,x1,y1
	variable d
	if (qq > 700)
		d=320
		x0=90
		y0=35
	else
		if (pp > 500)
			d=420
			x0=125
			y0=50
		else
			if (pp > 250)
				d=210
				x0=62
				y0=25
			else
				if (pp > 128)
					d=105
					x0=30
					y0=12
				endif
			endif
		endif
	endif
	x1=x0+d
	y1=y0+d
	DrawOval x0,y0,x1,y1	
end

Macro mask_MCPspot(width)
       variable width = 15
     	if (stringMatch(num2Str(pcsr(A)),"NaN")!=1)
     		variable x=pcsr(A)
     		variable y=qcsr(A)
     		print x,y
     		variable x1=round(x-width)
     		variable x2=round(x+width)
		variable y1=round(y-width)
		variable y2=round(y+width)
		M_ROIMask[x1,x2][y1,y2]=1
	endif
end

Macro loadXMCDsequence()
	 XMCDSequence()
end

function XMCDSequence() 
	string ctrlname
	SVAR scanNameStr = root:scanNameStr
	SVAR imageTypeStr = root:imageTypeStr
	NVAR frameno,layer
	String rootStr,fromStr,toStr,PathStr
	String fileNameStr=""
	variable step=1
	variable ref
	rootStr="2003_06_17_"
	fromStr="014"
	toStr="060"
	PathStr="F:data_2003_06_17:2003_06_17_"
	String tempStr
	Prompt rootStr, "Enter Root-Name: "
	Prompt fromStr, "from: "
	Prompt toStr,"to:"
	Prompt step,"step:"
	Prompt pathStr,"Path"
	DoPrompt "Load sequence", rootStr,fromStr,toStr,step,PathStr
	variable inputok = V_flag
	String normWaveStr
	Prompt normWaveStr, "enter norm wave",popup,WaveList("*",";","")
	DoPrompt "Choose normalisation wave", normWaveStr 
	variable doNormalise=V_Flag
	variable npNormalise = -1
	variable nqNormalise = -1
	variable maxbkg=1
	if (doNormalise==0)
		npNormalise = dimSize($normWaveStr,0)
		nqNormalise = dimSize($normWaveStr,1)
		redimension/D $normWaveStr
		Wavestats/Q $normWaveStr
		maxbkg=V_rms
	endif
	if (inputok==0)
		print "inputok"
		variable from=str2num(fromStr)
		variable to=str2num(toStr)
	 	variable len=strlen(fromStr)
	 	variable i=from
	 	variable index=0
	 	String sequenceNameStr=rootStr
	 	make/n=(1,1,1)/o/W/U $sequenceNameStr
	 	String sequenceInfoStr=rootStr+"_INFO"
	 	make/n=((to-from)/step+1)/o $sequenceInfoStr
	 	String tempInfoStr=rootStr+"_TEMP"
	 	make/n=((to-from)/step+1)/o $tempInfoStr
	 	variable coef
	 	if (stringMatch(imageTypeStr,"png")==1||stringMatch(imageTypeStr,"PNG")==1||stringMatch(imageTypeStr,"t16")==1||stringMatch(imageTypeStr,"tif")==1)
	 		coef=65536
	 	else
	 		coef=255
	 	endif
	 	variable np=-2,nq=-2
		do
			filenameStr=rootStr+returnSequenceNumber(i,len)+"_XMCD."+imageTypeStr 
			print "loading:"+filenameStr
       		filenameStr = OpenImage(filenameStr,0)
       		If (stringMatch(filenameStr,"")!=1)
       			if (index==0)
					np=DimSize($filenameStr,0)
					nq=DimSize($filenameStr,1)
					redimension/n=(np,nq,(to-from)/step+1) $sequenceNameStr
					Wavestats/Q $filenameStr
					coef=0.5*coef/V_max
				endif
				if (np==npNormalise&&nq==nqNormalise)
					backGroundCorrection($filenameStr,$normWaveStr,coef,maxbkg)
				endif
				redimension/W/U $sequenceNameStr
				imageToStack($sequenceNameStr,$filenameStr,index)
				entryToWave($sequenceInfoStr,i,index)
				tempStr=PathStr+returnSequenceNumber(i,len)+":2003_06_17_"+returnSequenceNumber(i,len)+"_LEEM-000.txt"
				print getTemperature(tempStr)
				entryToWave($tempInfoStr,getTemperature(tempStr),index)
				entryToWave($sequenceInfoStr,i,index)
       			killwaves/Z $filenameStr
       			index+=1
       		endif
       		i+=step
       	while ((i<=to))
       	if (index > 0)
       		redimension/n=(np,nq,index-1) $sequenceNameStr
       		redimension/n=(index-1) $sequenceInfoStr
			NewImage/S=2 $sequenceNameStr
			ModifyImage $sequenceNameStr plane=0
			PauseUpdate
			showinfo
			ControlBar 81
			ModifyGraph width=250*np/nq,height=250
			// controls
			String commandStr="addSequenceControls()"
			if ((stringMatch(num2Str(np),"NaN")!=1)&&(stringMatch(num2Str(nq),"NaN")!=1)) 
				Execute commandStr
			endif
		endif
	endif
End

Proc WannaKnowPhase()
	variable refnum
	Open/D/P=imagePath/R/D/M="wanna know phase" refNum 
	print S_filename, "\n"
	print , "has phase equal to: ",getPhase_old(S_filename)
End
 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// user macros
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Function userMacro1(ctrlname): buttoncontrol //160203
	String ctrlname
	Execute "DrawCristallographicAxes(,,,)" 	// Execute is used when you need to
									 	// run a macro from inside a function
End

Function userMacro2(ctrlname): buttoncontrol //160203
	String ctrlname
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	// tmpStr is the name of the image stack wave
	tmpStr= removeQuotes(tmpStr)
	variable selection=5
	prompt selection, "layer operations on the image stack",popup "displersive plane;profile;histogram;particle;LEED series; FFT series"
	DoPrompt "defect co-add", selection
	if (V_flag == 0) 
		if (selection == 1)
			// extracts a dispersive plane line profile for each layer in the stack; fits the profile
			// save the fit result, W_coef, in a 2D wave called FIT_tmpStr 
			dispersiveSeries(tmpStr)
		endif
		if (selection == 2)
			// plot the line profile for each layer in the stack 
			profileSeries(tmpStr)
		endif
		if (selection == 3)
			// calculates and fit the image histogram for each layer in the stack.
			histogramSeries(tmpstr)
		endif
		if (selection == 4)
			// calculates and fit the image histogram for each layer in layer... and much more!
			particleSeries(tmpstr)
		endif
		if (selection == 5)
			LEEDseries(tmpStr)
		endif
		if (selection == 6)
			FFTseries(tmpStr)
		endif
		
	endif
End

macro rename_kappa(name)
	String name="n_layers"
	String newname1 = "G_width_"+name
	duplicate/o  resultG $newname1
	String newname2 = "kappa_width_norm"+name
	duplicate/o k_norm $newname2
	//String newname3 = "kappa_parallel_"+name
	//duplicate/o  kappap $newname3
	//String newname4= "kappa_norm_"+name
	//duplicate/o  kappan $newname4
	//AppendToGraph $newname3 vs $newname4
	AppendToGraph $newname1 vs $newname2
	ModifyGraph rgb($newname1)=(0,12800,52224), marker($newname1)=8;DelayUpdate
	ModifyGraph mode($newname1)=3
	Label bottom "\\F'Arial'\\Z18k\\Bnorm \\M\\F'Arial'\\Z18(Å\\S-1\\M\\F'Arial'\\Z18)"
	Label left "\\F'Arial'\\Z18\\F'Symbol'\\F'Arial'Gaussian width (Å\\S-1\\M\\F'Arial'\\Z18)"
	ModifyGraph axOffset=0
	ModifyGraph margin=80,width=226.772,height=226.772
	SetAxis left 0,0.8 ;DelayUpdate
	SetAxis bottom 0,4
	Make/n=2/o lcoef={0,0.03}
	CurveFit/H="10"/NTHR=0 line kwCWave= lcoef,  $newname1 /X=$newname2 /D
end

Function Gau_2D(w,x,y) : FitFunc
	Wave w
	Variable x
	Variable y

	//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,y) = C+abs(A)*exp(-((x-x0)^2+(y-y0)^2)/w)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 2
	//CurveFitDialog/ x
	//CurveFitDialog/ y
	//CurveFitDialog/ Coefficients 5
	//CurveFitDialog/ w[0] = C
	//CurveFitDialog/ w[1] = A
	//CurveFitDialog/ w[2] = x0
	//CurveFitDialog/ w[3] = y0
	//CurveFitDialog/ w[4] = w

	return w[0]+abs(w[1])*exp(-((x-w[2])^2+(y-w[3])^2)/w[4])
End

Function Lor_2D(w,x,y) : FitFunc
	Wave w
	Variable x
	Variable y

	//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,y) = w1+abs(w2)/((x-w3)^2+(y-w4)^2+w5)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 2
	//CurveFitDialog/ x
	//CurveFitDialog/ y
	//CurveFitDialog/ Coefficients 7
	//CurveFitDialog/ w[0] = w1
	//CurveFitDialog/ w[1] = w2
	//CurveFitDialog/ w[2] = w3
	//CurveFitDialog/ w[3] = w4
	//CurveFitDialog/ w[4] = w5

	return w[0]+abs(w[1])/((x-w[2])^2+(y-w[3])^2+abs(w[4]))
End

Function GL_2D(w,x,y) : FitFunc
	Wave w
	Variable x
	Variable y

	//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,y) = w1+abs(w2)/((x-w3)^2+(y-w4)^2+w5)+abs(w5/((x-w2)^2+(y-w3)^2+w6))
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 2
	//CurveFitDialog/ x
	//CurveFitDialog/ y
	//CurveFitDialog/ Coefficients 7
	//CurveFitDialog/ w[0] = w1
	//CurveFitDialog/ w[1] = w2
	//CurveFitDialog/ w[2] = w3
	//CurveFitDialog/ w[3] = w4
	//CurveFitDialog/ w[4] = w5
	//CurveFitDialog/ w[5] = w6
	//CurveFitDialog/ w[6] = w7

	return w[0]+abs(w[1])*abs(w[5])*exp(-((x-w[2])^2+(y-w[3])^2)/w[4])+abs(w[5])/((x-w[2])^2+(y-w[3])^2+abs(w[6]))
End


Function LEEDseries(tmpStr)
	String tmpStr
	tmpStr=returnListFirstElement(imageNameList("",""),";")
	// tmpStr is the name of the image stack wave
	Wave W_coef = root:W_coef
	NVAR energy = root:energy
	NVAR layer = root:layer
	tmpStr= removeQuotes(tmpStr)
	variable np=dimSize($tmpStr,0)
	variable nq=dimSize($tmpStr,1)
	variable nr=dimSize($tmpStr,2)
	variable i=0

	// get cursor position
	variable pA = pcsr(A)  
	variable qA = qcsr(A)  
	variable L60eV=NaN	
	
	make/n=(np,nq)/o tmplayer=0
	extract2DWaveFromStack($tmpStr,0,tmplayer)
	make/n=(np) wtmp
	wtmp=tmplayer[p][qA]
	CurveFit lor wtmp
	
	variable fit_type=2 // 1 Gaussian, 2 lorenzian, 3 both
	Prompt fit_type, "enter peak shape",popup "Gaussian; Lorentzian; Gaussian + Lorentzian"
	variable fit_all=1
	Prompt fit_all, "fit all",popup "no; yes"
	variable verse
	Prompt verse, "start from?", popup "first; last"
	variable time_step=10
	Prompt time_step, "time step (s)"	
	variable delta=65 //round(sqrt(W_coef[3])*4.5)
	Prompt delta, "fit size (pixels)"
	String name="000"
	Prompt name, "wave name suffix"
	variable stw=1
	Prompt stw, "subtract transfer width", popup "yes; no"
	variable twidth=5
	Prompt twidth, "transfer width"
	DoPrompt "fit LEED scan", fit_type, fit_all, verse, time_step, delta, name, stw, twidth
	variable WF=1
	
	String ResultG ="Gau_width"+name
	String ResultL = "Lor_width"+name
	String ampG = "Gau_amp"+name
	String ampL = "Lor_amp"+name
	String chisq = "$chisq"+name
	String k_norm = "k_norm"+name
	String index = "index"+name
	
	make/n=(nr)/o $ResultG, $resultL, $ampG, $ampL, $index, $k_norm, $chisq, transferw
	wave GW=$ResultG
	wave LW=$ResultL
	wave GA=$ampG
	wave LA=$ampL
	wave IND=$index
	wave KN=$k_norm
	wave CHI=$chisq
	
	GW=NaN
	LW=NaN
	GA=NaN
	LA=NaN
	IND=NaN
	KN=NaN
	CHI=NaN
	transferw=NaN
	
	if (fit_type == 1)
	       W_coef={89.684,549.94,pA,qA,217.88}  // #100 7939330
	elseif (fit_type == 2)
			W_coef={60,2860,pA,qA,10}  // #100 chi=640000
			W_coef={51.44,25525,pA,qA,13.418}
 	elseif (fit_type == 3)
		//W_coef={67.681,0.0039498,235.23,248.44,28.154,63132,108.02}	// #100 343154
		W_coef={56.99,0.015648,pA,qA,18.705,6579.5,20} // graphite 20090211_003
	endif
	
	if (fit_all == 1)
		if (fit_type == 1)
			 W_coef={63.042,2000,pA,qA,9.9233}
		elseif (fit_type == 2)
			 W_coef={55.761,570.1,pA,qA,30}
			 //W_coef={51.44,25525,pA,qA,13.418}
			 W_coef={66.612,35819,pA,qA,5}
		elseif  (fit_type == 3)
			//W_coef={62.704,572.12,pA,qA,152.9,0.23169,-36229}
			//W_coef={56.99,0.015648,pA,qA,18.705,6579.5,20} 
			  W_coef={58.54,0.0022236,pA,qA,302.82,1.0572e+05,37.42}
		endif//
	endif
	
	make/n=1000/o pippo
	make/n=4 poppo={-0.018357,0.069301,23.078,13.963}
	
	display pippo

	variable first,last,step
	if (verse == 1)
		i=0
		first=0
		last=nr
		step=1 
	elseif (verse == 2)
		i=nr
		last=nr
		first=0
		step=-1
	endif
	
	do
		layer = i
		setGlobals(tmpstr,layer)
		extract2DWaveFromStack($tmpStr,i,tmplayer)
		
		if (i==first)
			NewImage/S=2 tmplayer
			ModifyImage tmplayer ctab= {*,*,Terrain,0}
			ShowInfo
			Cursor/I/A=1 A tmplayer pA-delta,qA-delta
			Cursor/I/A=1 B tmplayer pA+delta,qA+delta
		endif
			
			doUpdate
			
			//if ( (energy == 31 ) || (energy == 53) ||(energy == 75)  || (fit_all == 2) )  // single layer
			//if ( (energy == 34) || (energy == 43) ||(energy == 55)  || (energy == 65) || (fit_all == 2) )  // double layer
			//if  ( (energy == 22) ||  (energy == 27)  ||  (energy == 35)  ||  (energy == 41) || (energy == 55) ||  (energy == 82)  || (energy == 66) || (fit_all == 2)) // double layer
			//if ( (energy == 20)   || (energy == 40)  ||   (energy == 50) ||  (energy == 78)  || (fit_all == 2)) 
			if  ( (energy >= 20) &&  (energy <= 58) || (fit_all == 2))
			 			
				if (fit_type == 1)
					FuncFitMD/H="00000" Gau_2D W_coef  tmplayer[pcsr(A),pcsr(B)][qcsr(A),qcsr(B)] /D  
					GW[i]=sqrt(abs(W_coef[4])/2)*2.35/2    //gaussian half width
			       	LW[i]=NaN
			       	GA[i]=W_coef[1]
			       	LA[i]=NaN
					CHI[i]=V_chisq
					KN[i]=energy
					IND[i]=i
				
				elseif (fit_type == 2) // Lorenzian
					FuncFitMD/H="00000" Lor_2D W_coef  tmplayer[pcsr(A),pcsr(B)][qcsr(A),qcsr(B)] /D
					GW[i]=NaN     //gaussian half width
			      		LW[i]=sqrt(abs(W_coef[4]))
					GA[i]=NaN
					LA[i]=W_coef[1]
					CHI[i]=V_chisq
					KN[i]=energy
					IND[i]=i
				
				elseif (fit_type == 3) // Gaussian + Lorenzian
					FuncFitMD/H="0000000" GL_2D W_coef  tmplayer[pcsr(A),pcsr(B)][qcsr(A),qcsr(B)] /D
					GW[i]=sqrt(abs(W_coef[4])/2)*2.35/2    //gaussian half width
			      		LW[i]=sqrt(abs(W_coef[6]))
					GA[i]=W_coef[1]
					LA[i]=W_coef[5]
					CHI[i]=V_chisq
					KN[i]=energy
					IND[i]=i
					
					pippo=GL_2D(W_coef,p+W_coef[2],W_coef[3])
					Differentiate pippo
					Differentiate pippo
					Wavestats pippo
					CurveFit gauss  kwCWave=poppo, pippo[V_maxloc-5,V_maxloc+5] /D 
					LW[i] = poppo[2]
					GW[i] = poppo[2]
				endif	
		
		endif

		i=i+step

	while  (step*(i-last)<=0) 

	if (stw==1)
		twidth/=2 //pixels
		duplicate/o KN transferw
		transferw=twidth*sqrt(KN/60)   // KN still in eV!
		GW=sqrt(GW^2-transferw^2)
		LW=sqrt(LW^2-transferw^2)
	endif

	// k space calibration
	variable kspcali=0.0031*2*Pi  //1/2.464* 1/ 0.866025 = 0.468629  and  0.469/npixel
	GW*=kspcali 
	LW*= kspcali
	variable eV2J = 1.602E-19
	variable hbar = 1.101E-34
	variable m = 9.1e-31
	KN = sqrt( (KN-WF) * eV2J)
	KN = KN* (sqrt(2*m)/hbar) * 1E-10  // in A-1
	//converts IND to time
	IND*=time_step
	
	if (fit_all == 2)
		if (fit_type == 1)
			display GW vs IND
			Label bottom "\\F'Arial'\\Z18time (s)"
			Label left "\\F'Arial'\\Z18half width (2\\F'Symbol'p \\F'Arial'Å\\S-1\\M\\F'Arial'\\Z18)"
			ModifyGraph width=283.465,height=226.772
			ModifyGraph margin(left)=57,margin(bottom)=57,margin(top)=28,margin(right)=28
			SetAxis left 0,1;DelayUpdate
			SetAxis bottom 0,4
			display GA vs IND
			ModifyGraph mode=3
			ModifyGraph marker(GW)=19
			ModifyGraph width=283.465,height=226.772
			ModifyGraph margin(left)=57,margin(bottom)=57,margin(top)=28,margin(right)=28
		elseif (fit_type == 2)
			display LW vs IND
			ModifyGraph mode=3
			ModifyGraph marker(LW)=19
			ModifyGraph rgb(LW)=(0,12800,52224)
			Label bottom "\\F'Arial'\\Z18time (s)"
			Label left "\\F'Arial'\\Z18half width (2\\F'Symbol'p \\F'Arial'Å\\S-1\\M\\F'Arial'\\Z18)"
			ModifyGraph width=283.465,height=226.772
			ModifyGraph margin(left)=57,margin(bottom)=57,margin(top)=28,margin(right)=28
			SetAxis left 0,1;DelayUpdate
			SetAxis bottom 0,4
			display LA vs IND
			ModifyGraph mode=3
			ModifyGraph marker(LA)=19	
			ModifyGraph rgb(LA)=(0,12800,52224)
			Label bottom "\\F'Arial'\\Z18time (s)"
			ModifyGraph width=283.465,height=226.772
			ModifyGraph margin=57
		elseif (fit_type == 3)
			display GW vs IND
			appendToGraph LW vs IND
			ModifyGraph mode=3
			ModifyGraph marker(GW)=19
			ModifyGraph marker(LW)=19
			SetAxis left 0,1;DelayUpdate
			SetAxis bottom 0,4
			Label bottom "\\F'Arial'\\Z18time (s)"
			Label left "\\F'Arial'\\Z18half width (2\\F'Symbol'p \\F'Arial'Å\\S-1\\M\\F'Arial'\\Z18)"
			ModifyGraph rgb(LW)=(0,12800,52224)
			ModifyGraph width=283.465,height=226.772
			ModifyGraph margin(left)=57,margin(bottom)=57,margin(top)=28,margin(right)=28
			display GA vs IND
			appendToGraph/R LA vs IND
			ModifyGraph mode=3
			ModifyGraph marker(GA)=19
			ModifyGraph marker(LA)=19			
			ModifyGraph rgb(LA)=(0,12800,52224)
			Label bottom "\\F'Arial'\\Z18time (s)"
			ModifyGraph width=283.465,height=226.772
			ModifyGraph margin(left)=57,margin(bottom)=57,margin(top)=28,margin(right)=28
		endif
		display CHI vs IND
		Label left "\\F'Arial'\\Z18\\F'Symbol'c\\F'Arial'\\S2"
		Label bottom "\\F'Arial'\\Z18time (s)"
		ModifyGraph width=283.465,height=226.772
		ModifyGraph margin(left)=57,margin(bottom)=57,margin(top)=28,margin(right)=28
	else
		display GW vs KN
		appendToGraph LW vs KN
		Label bottom "\\F'Arial'\\Z18k\\Bnorm"
		Label left "\\F'Arial'\\Z18half width (2\\F'Symbol'p \\F'Arial'Å\\S-1\\M\\F'Arial'\\Z18)"
		ModifyGraph rgb(LW)=(0,12800,52224)
		ModifyGraph width=283.465,height=226.772
		ModifyGraph mode=3
		ModifyGraph marker=19
		ModifyGraph margin(left)=57,margin(bottom)=57,margin(top)=28,margin(right)=57
		SetAxis left 0,1
		SetAxis bottom 0,4.5
		
		Make/n=2/o W_line={0,0.03}
		CurveFit/H="10" line kwCWave=W_line,  GW /X=KN /D 
		variable thetaG=0.5*asin(W_line[1])*180/PI
		Make/n=2/o W_line={0,0.03}
		CurveFit/H="10" line kwCWave=W_line,  LW /X=KN /D 
		variable thetaL=0.5*asin(W_line[1])*180/PI
		String text=""
		if (fit_type == 1)
			text="\\F'Symbol'\\Z14Dq\\F'arial'\\Z14(G) = "+num2str(thetaG) 
		elseif (fit_type == 2)
			text=" \\F'Symbol'\\Z14Dq\\F'arial'\\Z14(L) = "+num2str(thetaL)
		elseif (fit_type == 3)
			text="\\F'Symbol'\\Z14Dq\\F'arial'\\Z14(G) = "+num2str(thetaG) +" \\F'Symbol'\\Z14Dq\\F'arial'\\Z14(L) = "+num2str(thetaL)
		endif
		TextBox/C/N=text1/F=0/A=LT text

		variable K60eV=sqrt((60-WF)*eV2J)*(sqrt(2*m)/hbar)*1E-10
		L60eV=sqrt((W_line[1]*K60eV)^2+(twidth/2*kspcali)^2)
		
		display GA vs KN
		appendToGraph/R LA vs KN
		ModifyGraph mode=3
		ModifyGraph marker=19
		ModifyGraph rgb(LA)=(0,12800,52224)
		Label bottom "\\F'Arial'\\Z18k\\Bnorm"
		ModifyGraph width=283.465,height=226.772
		ModifyGraph margin(left)=57,margin(bottom)=57,margin(top)=28,margin(right)=57
		display $chisq vs KN
		Label left "\\F'Arial'\\Z18\\F'Symbol'c\\F'Arial'\\S2"
		Label bottom "\\F'Arial'\\Z18k\\Bnorm"
		ModifyGraph mode=3
		ModifyGraph marker(CHI)=19
		ModifyGraph width=283.465,height=226.772
		ModifyGraph margin(left)=57,margin(bottom)=57,margin(top)=28,margin(right)=28	
	endif
	
	print "Lorentzian width at 60 eV: ",L60eV
	
End

function E2k(x)
	variable x
	variable eV2J = 1.602E-19
	variable hbar = 1.101E-34
	variable m = 9.1e-31
	variable kn
	kn = sqrt( (x-2) * eV2J)
	return kn* (sqrt(2*m)/hbar) * 1E-10  // in A-1
end

// analysis of series of images
Function histogramSeries(tmpstr)
	String tmpstr // series name
	NVAR layer = root:layer
	NVAR frameno = root:frameno
	wave W_fit=root:W_fit
	variable n=dimsize($tmpStr,2)
	variable i=0
	String hstStr
	String fitStr="FIT_"+tmpStr
	make/n=(20,n)/o $fitStr
	wave wfit=$fitStr
	wfit=NaN
	do
		layer = i
		setGlobals(tmpstr,layer)
		CalcFrameHistogram(tmpStr)
		// fit
		hstStr = tmpstr+"_"+num2ScanNumberStr(layer)+"_HST"
		//if (i < 13) //20020420
			Fit_Historgram_2Level(hstStr,W_coef)
		//else
			//Fit_Historgram_3Level(hstStr,W_coef)
		//endif
		if (W_fit[0]<2e6) // chisquare
			wfit[][i]=W_fit[p]
		else
			wfit[][i]=NaN
		endif
		wfit[19][i]=frameno
		doUpdate
		// kill window and waves therein
 		killSeries("")
 		doUpdate
		i+=1
	while (i<n)
end

// analysis of series of images
Function particleSeries(tmpstr)
	String tmpstr // series name
	NVAR layer = root:layer
	NVAR frameno = root:frameno
	wave W_fit=root:W_fit
	variable n=dimsize($tmpStr,2)
	variable i=0
	String hstStr
	String fitStr="FIT_"+tmpStr
	make/n=(20,n)/o $fitStr
	wave wfit=$fitStr
	wfit=NaN
	variable thrs12=NaN
	variable thrs13=NaN
	// make the tmp wave for image particle analysis
	variable np=dimsize($tmpStr,0)
	variable nq=dimsize($tmpStr,1)
	duplicate/o root:M_ROIMask neutralmask
	imagestats neutralmask
	variable roiarea=np*nq*(1-V_avg)
	variable totArea=0
	variable particlearea=10
	duplicate/o M_ROIMask myROI
	variable j
	do
		layer = i
		setGlobals(tmpstr,layer)
		CalcFrameHistogram(tmpStr)
		// fit
		hstStr = tmpstr+"_"+num2ScanNumberStr(layer)+"_HST"
		//Fit_Historgram_WY(hstStr)
		Fit_Historgram_3Level(hstStr,W_coef)
		//saves coeff
		wfit[][i]=W_fit[p]
		wfit[19][i]=frameno
		doUpdate
		thrs12=0.6*wfit[2]+0.4*wfit[1]
		thrs13=0.6*wfit[2]+0.4*wfit[3]
		j=0
		do
		// analyse dark
			make/n=(np,nq)/o tmplayer=0
			make/n=(np,nq)/o tmpmask=0
			//  copy layer into tmp
			extract2DWaveFromStack($tmpStr,i,tmplayer)
			extract2DWaveFromStack($tmpStr,i,tmpmask)
			redimension/b/u tmpmask
			// dark
			ImageThreshold/o/Q/T=(thrs12)/M=0/R=M_ROIMask tmpmask 
			ImageAnalyzeParticles/E/M=3/Q/A=(particlearea)/o stats, tmpmask
			// calculates the avg value within mask
			duplicate/o root:M_Particle tmpmask
			tmpmask-=16		
			imageStats/M=1/R=tmpmask tmplayer
			doupdate
			wfit[11][i]=V_avg
			// total area
			wfit[14][i]=sum(root:W_ImageObjArea, 0, numpnts(root:W_ImageObjArea) )/roiarea
			// start creating mask for grey
			duplicate/o M_ROIMask neutralmask
			neutralmask*=(208)
			neutralmask+=(tmpmask[p][q]-48)
		// analyse bright
			//  copy layer into tmp
			extract2DWaveFromStack($tmpStr,i,tmplayer)
			extract2DWaveFromStack($tmpStr,i,tmpmask)
			redimension/b/u tmplayer
			// bright
			ImageThreshold/O/I/Q/T=(thrs13)/M=0/R=M_ROIMask tmpmask 
			ImageAnalyzeParticles/E/M=3/Q/A=(particlearea)/o stats, tmpmask
			duplicate/o root:M_Particle tmpmask
			tmpmask-=16
			imageStats/M=1/R=tmpmask tmplayer
			wfit[13][i]=V_avg
			// total area
			wfit[16][i]=sum(root:W_ImageObjArea, 0, numpnts(root:W_ImageObjArea) )/roiarea
			// neutral
			neutralmask+=(tmpmask[p][q]-48)
			imageStats/M=1/R=neutralmask tmplayer
			wfit[12][i]=V_avg
			ImageAnalyzeParticles/E/W/Q/M=3/A=(particlearea) stats, neutralmask
			// total area
			wfit[15][i]=sum(root:W_ImageObjArea, 0, numpnts(root:W_ImageObjArea) )/roiarea
 			duplicate/o neutralmask goophy
 			killwaves/z tmplayer,tmpmask,neutralmask
 			thrs12=0.6*wfit[12]+0.4*wfit[11]
			thrs13=0.6*wfit[12]+0.4*wfit[13]
 			doUpdate
 			j+=1
 		while (j<5)
 		// kill window and waves therein
 		killSeries("")
		i+=1
	while (i<n)
	//killwaves/z tmplayer,tmpmask,neutralmask
end

Function profileSeries(ctrlname): buttoncontrol //160203
	String ctrlname
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	//  globals
	NVAR lastcursAp = root:lastcursAp
	NVAR lastcursBp = root:lastcursBp
	NVAR lastcursAq = root:lastcursAq
	NVAR lastcursBq = root:lastcursBq
	NVAR layer=root:layer
	NVAR frameno=root:frameno
	NVAR energy=root:energy
	WAVE W_coef=root:W_coef
	// cursors
	variable p1,p2,q1,q2
	if (stringMatch(num2Str(pcsr(A)),"NaN")!=1 && stringMatch(num2Str(pcsr(B)),"NaN")!=1)
 		p1=pcsr(A)
		p2=pcsr(B)
		q1=qcsr(A)
		q2=qcsr(B)
		// store curs pos
		lastcursAp=pcsr(A)
		lastcursBp=pcsr(B)
		lastcursAq=qcsr(A)
		lastcursBq=qcsr(B)
	else
		p1= lastcursAp
		p2 = lastcursBp
		q1= lastcursAq
		q2 = lastcursBq
	endif
	// profile width
	variable halfwidth=getProfileHalfWidth()
	// loop
	variable n=dimsize($tmpStr,2)
	variable i=0
	String testStr,xwaveStr,proStr
	variable np
	String waterfallStr="WF_"+tmpStr
	make/n=(1,1)/o $waterfallStr
	wave w=$waterfallStr
	variable de=1
	make/n=(n) result=NaN
	wave r=result
	do
		testStr=tmpStr+"_INFO"
		xwaveStr=""
		if (exists(testStr)==1)
			xwaveStr=testStr
			frameno=returnWaveElement($xwaveStr,i)
		endif
		testStr=tmpStr+"_STV"	
		if (exists(testStr)==1)
			xwaveStr=testStr
			energy=returnWaveElement($xwaveStr,i)
			variable energy0=returnWaveElement($xwaveStr,0)
			variable energy1=returnWaveElement($xwaveStr,1)
			de= energy1-energy0
		endif
		testStr=tmpStr+"_NRG"
		if (exists(testStr)==1)
			xwaveStr=testStr
			energy=returnWaveElement($xwaveStr,i)
		endif
		testStr=tmpStr+"_OBJ"
		if (exists(testStr)==1)
			xwaveStr=testStr
			energy=returnWaveElement($xwaveStr,i)
		endif
		layer=i
		getProfile(tmpStr,halfwidth,p1,p2,q1,q2,i,i)
		proStr = tmpStr+"_FPR_"+num2ScanNumberStr(i)
		waveStats $proStr
		if (i<16)
			W_coef={75.157,V_max*100,V_maxloc,119.38}
			//CurveFit lor kwCWave=W_coef,  $proStr /D
			redimension/n=5 W_coef
		else
			W_coef[2]=V_maxloc
			//FuncFit fit_MDC W_coef $proStr /D
			r[i] = w_coef[2] 
		endif
		
		
		
		//W_coef[1]=V_avg*100
		
		//CurveFit lor kwCWave=W_coef,  $proStr /D
		//redimension/n=5 W_coef
		
		//divideWave($proStr,V_max)
		if (i==0)
			np = numpnts($proStr)
			print "i eq 0:", proStr, np, n
			redimension/n=(np,n) w
		endif
		writecol(waterfallStr,proStr,i)
		
		i+=1
	while(i<n)
	// cosmetics
	NewImage/S=2 $waterfallStr
	//ModifyImage $waterfallStr ctab= {*,*,Yellow,0}
	//SetScale/P x 0,45.977,"", $waterfallStr
	//SetScale/P y 0,10,"", $waterfallStr
	SetScale/P y energy0,de,"", $waterfallStr
	SetAxis/A left
	Label left "\\F'Arial'\\Z16energy (eV)";DelayUpdate; 
	Label top "\\F'Arial'\\Z16Profile (pixels)"
	ModifyGraph margin=38
	ModifyGraph width=283.465,height=283.465
	// controls
	execute "addBWcontrols()"
End

Function fit_MDC(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) = w1+abs(w2)/((x-w3)^2+w4)+w5*x
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 2
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 7
	//CurveFitDialog/ w[0] = w1
	//CurveFitDialog/ w[1] = w2
	//CurveFitDialog/ w[2] = w3
	//CurveFitDialog/ w[3] = w4
	//CurveFitDialog/ w[4] = w5

	return w[0]+abs(w[1])/((x-w[2])^2+abs(w[3]))+w[4]*x

end

Function dispersiveSeries(ctrlname): buttoncontrol //160203
	String ctrlname
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	//  globals
	NVAR lastcursAp = root:lastcursAp
	NVAR lastcursBp = root:lastcursBp
	NVAR lastcursAq = root:lastcursAq
	NVAR lastcursBq = root:lastcursBq
	NVAR layer=root:layer
	NVAR frameno=root:frameno
	NVAR energy=root:energy
	// cursors
	variable p1,p2,q1,q2
	if (stringMatch(num2Str(pcsr(A)),"NaN")!=1 && stringMatch(num2Str(pcsr(B)),"NaN")!=1)
 		p1=pcsr(A)
		p2=pcsr(B)
		q1=qcsr(A)
		q2=qcsr(B)
		// store curs pos
		lastcursAp=pcsr(A)
		lastcursBp=pcsr(B)
		lastcursAq=qcsr(A)
		lastcursBq=qcsr(B)
	else
		p1= lastcursAp
		p2 = lastcursBp
		q1= lastcursAq
		q2 = lastcursBq
	endif
	// profile width
	variable halfwidth=getProfileHalfWidth()
	// variables
	String xStr,sxStr,syStr,dispStr
	String testStr,xwaveStr
	variable bkg,peakarea
	variable n=dimsize($tmpStr,2)
	variable m=1
	variable i=0
	// waterfall imageplot
	String waterfallStr="WF_"+tmpStr
	make/n=(1,n)/o $waterfallStr
	wave w=$waterfallStr
	// fit parameter array
	String fitStr = "FIT_"+tmpStr
	make/n=(20,n)/o $fitStr
	wave wfit=$fitStr
	// other parameters
	String statStr="R_"+tmpStr
	make/n=(2,n)/o $statStr
	wave wresult=$statStr 
	variable avg0, avg
	do
		testStr=tmpStr+"_INFO"
		xwaveStr=""
		if (exists(testStr)==1)
			xwaveStr=testStr
			frameno=returnWaveElement($xwaveStr,i)
		endif
		testStr=tmpStr+"_STV"	
		if (exists(testStr)==1)
			xwaveStr=testStr
			energy=returnWaveElement($xwaveStr,i)
		endif
		testStr=tmpStr+"_NRG"
		if (exists(testStr)==1)
			xwaveStr=testStr
			energy=returnWaveElement($xwaveStr,i)
		endif
		testStr=tmpStr+"_OBJ"
		if (exists(testStr)==1)
			xwaveStr=testStr
			energy=returnWaveElement($xwaveStr,i)
		endif
		// select the i-th layer of the sequence
		layer=i
		// from the i-th layer from the stack the line profile is extracted and converted to KE 
		dispStr = stackLayerToDispersiveKE(tmpStr,halfwidth,i)
		if (i==0)
			m=numpnts($dispStr)
			redimension/n=(numpnts($dispStr),n) w
			avg0=mean($dispStr,m-50,m-1)
		endif
		avg=mean($dispStr,m-50,m-1)
		duplicate/o $dispstr tmp
		tmp=avg0/avg*tmp
		duplicate/o tmp $dispStr
		// Shirley wave names
		sxStr=dispStr+"_SX" //energy 
		syStr=dispStr+"_SY" //intensity
		xStr=dispStr+"_X"    
		// remove Shirley bkg from wave.  Calculates bkg avg value and peak integral
		bkg=makeShirelyWave(dispStr)
		peakarea=calcAreaShirley(sxStr,syStr,pcsr(A),pcsr(B))
		// store data in the waterfall 2D-wave
		
		writecol(waterfallStr,dispStr,i)
		wresult[0][i]=peakarea
		wresult[1][i]=bkg
		// fit the dispersive plane line profile
		if (i==0)
			//testFitFunction($syStr,$sxStr)
			//doublet_KE($syStr,$sxStr,1)
			//single_KE($syStr,$sxStr,1)
		else
			//testFitFunction($syStr,$sxStr)
			//doublet_KE($syStr,$sxStr,1)
			//single_KE($syStr,$sxStr,1)
			//Au4f_XY_KE(syStr)
		endif
		
		// store the fit result
		writecol(fitStr,"W_coef",i)
		// update plots
		doUpdate
		i+=1
	while(i<n)
	// cosmetics on waterfall
	NewImage/S=2 $waterfallStr
	ModifyImage $waterfallStr ctab= {*,*,YellowHot,0}
	SetScale/P x 0,1,"", $waterfallStr
	SetScale/P y 0,1 ,"", $waterfallStr
	Label left "\\F'Arial'\\Z16scan number";DelayUpdate; 
	Label top "\\F'Arial'\\Z16pixel"
	ModifyGraph margin=38
	Execute "addBWcontrols()"
End

function getFrequency(w, fov,delta,ptest,qtest)
	wave w
	variable fov
	variable delta
	variable ptest
	variable qtest

	Wave W_coef = root:W_coef
	Wave W_sigma= root:W_sigma
	imageStats w
	variable pmax = V_maxRowLoc
	variable qmax = V_maxColLoc
	if (ptest !=  0)
		pmax = ptest
	endif
	
	if (qtest != 0)
		qmax = qtest
	endif	
	variable np=dimSize(w,0)
	variable nq=dimSize(w,1)
	print pmax,qmax
	//W_coef={3.8431e+05,-3.0581e+06,5.5974,344.13,1.5318}
	//W_coef[2]=pmax
	//W_coef[3]=qmax
	//FuncFitMD/H="00110" Gau_2D W_coef  w[pmax-6,pmax+6][qmax-6,qmax+6] /D
	//FuncFitMD Gau_2D W_coef  w[pmax-4,pmax+4][qmax-4,qmax+4] /D 
	//pmax=W_coef[2]/(np*2) //*256/np
	//qmax=(nq-W_coef[3])/nq //*128/nq
	
	  //W_coef={64841,3.6581e+06,6.9213,0.77347,295.96,0.79846,-0.14126}
	  
	 //W_coef={24786,8.8213e+05,34.614,17.256,297.82,7.0834,-0.95392}
	  //W_coef[2]=pmax
	 //W_coef[4]=qmax
	W_coef={21211,48798,68.745,24.078,4.0489,3.2188,0.34135}  //061
	// W_coef={24786,8.8213e+05,34.614,17.256,297.82,7.0834,-0.95392} //060
	 
	variable pmin= max(0,pmax-delta)
	variable qmin= max(0,qmax-delta)
	CurveFit Gauss2D  kwCWave=W_coef,  w[pmin,pmax+delta][qmin,qmax+delta] /D 
	pmax=W_coef[2]/(np*2) //*256/np
	//qmax=(nq-W_coef[4])/nq //*128/nq
	
	qmax=(W_coef[4])/nq
	variable d=sqrt(pmax^2+qmax^2) 	
	variable error1=sqrt((W_sigma[2]/(2*np)*2*pmax/d)^2+(W_sigma[4]/(nq)*2*qmax/d)^2)
	variable error2 =(fov*1000/435) *error1/d
	print "error = ", error2
	return 1/sqrt(pmax^2+qmax^2)*fov*1000/435	
end

Function getFFT(w,fov,delta,ptest,qtest)
	wave w
	variable fov
	variable delta
	variable ptest
	variable qtest
	variable np=dimSize(w,0)
	variable nq=dimSize(w,1)
	variable mp=(round(np/2)*2)
	variable mq=(round(np/2)*2)
	make/n=(mp,mq)/o w1, hann
	//variable mr=round(np/2-mp/2)
	//variable mc=round(nq/2-mq/2)  // +nq/12 burn in MCP
	//w1=w[mr+p][mc+q]
	w1=w[p][q]
	wavestats/Q w1
	w1-=V_avg
	// windowing (for making nice FFTs)
	hann=sin(pi*p/mp)*sin(pi*q/mq)
	w1*=hann
	FFT/OUT=3/Z/DEST=tfft w1
	redimension/R tfft
	duplicate/o tfft w1
	killwaves/z hann,tfft
	SetScale/P x 0,1,"", w1;DelayUpdate
	SetScale/P y 0,1,"", w1
	return getFrequency(w1,fov,delta,ptest,qtest)
end

function getArea(w)
	wave w
	variable thres=2
	Imagestats w
	return 1
end

Function FFTseries(tmpStr)
	String tmpStr
	tmpStr=returnListFirstElement(imageNameList("",""),";")
	// tmpStr is the name of the image stack wave
	Wave W_coef = root:W_coef
	NVAR energy = root:energy
	NVAR layer = root:layer
	tmpStr= removeQuotes(tmpStr)
	
	variable np=dimSize($tmpStr,0)
	variable nq=dimSize($tmpStr,1)
	variable nr=dimSize($tmpStr,2)
	variable i

	make/n=(np,nq)/o tmplayer=0
	extract2DWaveFromStack($tmpStr,0,tmplayer)
	String ResultP ="Period_"+tmpStr
	String index = "index"+tmpStr
	
	make/n=(nr)/o $ResultP, $index
	wave PW=$ResultP
	wave IND=$index
	
	variable startfrom=30
	Prompt startfrom, "start from frame"
	variable fov=4
	Prompt fov, "field of view (um)"
	variable delta=30
	Prompt delta, "fit interval (pixels)"
	variable ptest=68
	Prompt ptest, "p test"
	variable qtest=4
	Prompt qtest, "q test"
	DoPrompt "period 2D FFT", startfrom, fov, delta, ptest, qtest

	PW=NaN
	IND=NaN
	
	i=0
	do
		layer = i
		setGlobals(tmpstr,layer)
		
		if (i >= startfrom)
		extract2DWaveFromStack($tmpStr,i,tmplayer)
		PW[i] = getFFT(tmplayer,fov,delta,ptest,qtest)
		if (PW[i]<10)
			PW[i]=NaN
		endif
		//print GW[i]
		
		//GW[i]=getFrequency(tmplayer, 2)
		Else
		PW[i] = NaN
		endif
		doUpdate
			
	
		i=i+1
	while  (i<=nr) 
	
End

Macro Style3(Axisbottom,Axisleft)
	String Axisbottom="Temperature (C)"
	String Axisleft="Stripe period (nm)"
	ModifyGraph mirror(left)=2,mirror(bottom)=1,nticks(bottom)=10;DelayUpdate
	Label left "\\F'Arial'\\Z18"+axisleft+"";DelayUpdate
	Label bottom "\\F'Arial'\\Z18"+Axisbottom+""
	ModifyGraph margin(left)=62,margin(bottom)=62,margin(top)=14,margin(right)=14;DelayUpdate
	ModifyGraph width=250,height=250	
	ModifyGraph fSize=14
End




/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// fitting functions and macros
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// basic bricks: doniac-sunijic
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// Binding energy
Function DS_BE(alpha,gammaL,pos,Ampli,x) 
	Variable alpha,gammaL,pos,Ampli
	Variable x
	alpha=abs(alpha)
	return abs(Ampli) * exp (gammln(1-alpha))*cos( 0.5*Pi*alpha+(1-alpha)*atan(2*(pos-x)/gammaL))/((x-pos)^2+ 4*gammaL^2)^(0.5*(1-alpha))
End

// kinetic energy
Function DS_KE(alpha,gammaL,pos,Ampli,x) 
	Variable alpha,gammaL,pos,Ampli
	Variable x
	alpha=abs(alpha)
	//return abs(Ampli) * exp (gammln(1-alpha,1e-7))*cos( 0.5*Pi*alpha-(1-alpha)*atan(2*(pos-x)/gammaL))/((x-pos)^2+ 4*gammaL^2)^(0.5*(1-alpha))
	return abs(Ampli) * exp (gammln(1-alpha,1e-7))*cos( 0.5*Pi*alpha-(1-alpha)*atan(2*(pos-x)/gammaL))/((x-pos)^2+ 4*gammaL^2)^(0.5*(1-alpha))
End

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// multipeak 
Function multipeakCDS_BE(w,x) : FitFunc
	Wave w
	Variable x
	Variable npeaks = (numpnts(w)-6)/2
	Variable i=0
	Variable j=0
	Variable sigma = w[3]/(2*sqrt(2*ln(2)))  
	Variable nsteps = 50
	Variable integral = 0
	Variable start = -4*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable sumDS,poli,y,posi
	Variable normG = 1/sqrt(2*Pi*sigma^2)
	do
		y=x-start
		sumDS=0
		for(i=0;i<npeaks;i+=1)
			posi=w[6]
			if (i>0)
				posi=w[6]+w[6+2*i]
			endif
			sumDS+=DS_BE(w[5],w[4],posi,w[7+2*i],y)
		endfor
             poli=w[0]+w[1]*y+w[2]*y^2
             integral+=normG*exp(-(start)^2/(2*sigma^2))*(sumDS+poli)*delta
		start+=delta 	
	while (start <= stop)
	return integral
end

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// one peak
Function onepeakCDS_KE(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) = sum(i) sum (DS (x(i)-y) * exp () * dy)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 9
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 8
	//CurveFitDialog/ w[0] = k0
	//CurveFitDialog/ w[1] = k1
	//CurveFitDialog/ w[2] = k2
	//CurveFitDialog/ w[3] = gaussfwhm
	//CurveFitDialog/ w[4] = gamma
	//CurveFitDialog/ w[5] = alpha
	//CurveFitDialog/ w[6] = pos1
	//CurveFitDialog/ w[7] = A1
	
	Variable i=0
	Variable j=0
	Variable sigma = w[3]/(2*sqrt(2*ln(2)))  
	Variable nsteps = 50
	Variable integral = 0
	Variable start = -4*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable ds1,poli,y
	Variable normG = 1/sqrt(2*Pi*sigma^2)
	do
		y=x-start
		ds1=DS_KE(w[5],w[4],w[6],w[7],y)
             poli=w[0]+w[1]*y+w[2]*y^2
             integral+=normG*exp(-(start)^2/(2*sigma^2))*(ds1+poli)*delta
		start+=delta 	
	while (start <= stop)
	return integral
end

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// two peaks; same Gauss,Lor,alpha
Function twopeakCDS_BE(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) = sum(i) sum (DS (x(i)-y) * exp () * dy)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 10
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 10
	//CurveFitDialog/ w[0] = k0
	//CurveFitDialog/ w[1] = k1
	//CurveFitDialog/ w[2] = k2
	//CurveFitDialog/ w[3] = gaussfwhm
	//CurveFitDialog/ w[4] = gamma
	//CurveFitDialog/ w[5] = alpha
	//CurveFitDialog/ w[6] = pos1
	//CurveFitDialog/ w[7] = A1
	//CurveFitDialog/ w[8] = deltapos2
	//CurveFitDialog/ w[9] = A2
	
	Variable i=0
	Variable j=0
	Variable sigma = w[3]/(2*sqrt(2*ln(2)))  
	Variable nsteps = 80
	Variable integral = 0
	Variable start = -4*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable ds1,ds2,poli,y
	Variable normG = 1/sqrt(2*Pi*sigma^2)
	Variable pos2 = w[6]+w[8]
	do
		y=x-start
		ds1=DS_BE(w[5],w[4],w[6],w[7],y)
             ds2=DS_BE(w[5],w[4],pos2,w[9],y)
             poli=w[0]+w[1]*y+w[2]*y^2
             integral+=normG*exp(-(start)^2/(2*sigma^2))*(ds1+ds2+poli)*delta
		start+=delta 	
	while (start <= stop)
	return integral
end

// two peaks; same Gauss,Lor,alpha
Function   twopeakCDS_KE(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) = sum(i) sum (DS (x(i)-y) * exp () * dy)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 10
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 11
	//CurveFitDialog/ w[0] = k0
	//CurveFitDialog/ w[1] = k1
	//CurveFitDialog/ w[2] = k2
	//CurveFitDialog/ w[3] = gaussfwhm
	//CurveFitDialog/ w[4] = gamma
	//CurveFitDialog/ w[5] = alpha
	//CurveFitDialog/ w[6] = pos1
	//CurveFitDialog/ w[7] = A1
	//CurveFitDialog/ w[8] = deltapos2
	//CurveFitDialog/ w[9] = A2
	//CurveFitDialog/ w[9] = gmma2
	Variable i=0
	Variable j=0
	Variable sigma = w[3]/2.35 //     (2*sqrt(2*ln(2)))  
	// s2/2 = g2/4*ln2
	Variable nsteps = 100
	Variable integral = 0
	Variable start = -3*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable ds1,ds2,poli,y
	Variable normG = 1/sqrt(2*Pi*sigma^2)
	Variable pos2 = w[6]+w[8]
	do
		y=x-start
		ds1=DS_KE(w[5],w[4],w[6],w[7],y)
             	ds2=DS_KE(w[5],w[11],pos2,w[9],y)
             	poli=0 //w[0]+w[1]*y+w[2]*y^2
             	integral+=normG*exp(-(start)^2/(2*sigma^2))*(ds1+ds2+poli)*delta
		start+=delta 	
	while (start <= stop)
	integral+=w[0]+w[1]*x+w[2]*x^2
	return integral
end

threepeakCDS_KE
Function   threepeakCDS_KE(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) = sum(i) sum (DS (x(i)-y) * exp () * dy)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 12
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 12
	//CurveFitDialog/ w[0] = k0
	//CurveFitDialog/ w[1] = k1
	//CurveFitDialog/ w[2] = k2
	//CurveFitDialog/ w[3] = gaussfwhm
	//CurveFitDialog/ w[4] = gamma
	//CurveFitDialog/ w[5] = alpha
	//CurveFitDialog/ w[6] = pos1
	//CurveFitDialog/ w[7] = A1
	//CurveFitDialog/ w[8] = deltapos2
	//CurveFitDialog/ w[9] = A2
	//CurveFitDialog/ w[10] = deltapos3
	//CurveFitDialog/ w[10] = A3
	Variable i=0
	Variable j=0
	Variable sigma = w[3]/2.35 //     (2*sqrt(2*ln(2)))  
	// s2/2 = g2/4*ln2
	Variable nsteps = 100
	Variable integral = 0
	Variable start = -3*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable ds1,ds2,ds3,poli,y
	Variable normG = 1/sqrt(2*Pi*sigma^2)
	Variable pos2 = w[6]+w[8]
	Variable pos3 = w[6]+w[10]
	do
		y=x-start
		ds1=DS_KE(w[5],w[4],w[6],w[7],y)
             	ds2=DS_KE(w[5],w[4],pos2,w[9],y)
             	ds3=DS_KE(w[5],w[4],pos3,w[11],y)
             	poli=0 //w[0]+w[1]*y+w[2]*y^2
             	integral+=normG*exp(-(start)^2/(2*sigma^2))*(ds1+ds2+ds3+poli)*delta
		start+=delta 	
	while (start <= stop)
	integral+=w[0]+w[1]*x+w[2]*x^2
	return integral
end

// two peaks; same Gauss,Lor,alpha
Function   fourpeakCDS_KE(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) = sum(i) sum (DS (x(i)-y) * exp () * dy)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 12
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 14
	//CurveFitDialog/ w[0] = k0
	//CurveFitDialog/ w[1] = k1
	//CurveFitDialog/ w[2] = k2
	//CurveFitDialog/ w[3] = gaussfwhm
	//CurveFitDialog/ w[4] = gamma
	//CurveFitDialog/ w[5] = alpha
	//CurveFitDialog/ w[6] = pos1
	//CurveFitDialog/ w[7] = A1
	//CurveFitDialog/ w[8] = deltapos2
	//CurveFitDialog/ w[9] = A2
	//CurveFitDialog/ w[10] = deltapos3
	//CurveFitDialog/ w[11] = A3
	//CurveFitDialog/ w[12] = deltapos4
	//CurveFitDialog/ w[13] = A4
	Variable i=0
	Variable j=0
	Variable sigma = w[3]/2.35 //     (2*sqrt(2*ln(2)))  
	// s2/2 = g2/4*ln2
	Variable nsteps = 100
	Variable integral = 0
	Variable start = -3*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable ds1,ds2,ds3,ds4,poli,y
	Variable normG = 1/sqrt(2*Pi*sigma^2)
	Variable pos2 = w[6]+w[8]
	Variable pos3 = w[6]+w[10]
	Variable pos4 = w[6]+w[12]
	do
		y=x-start
		ds1=DS_KE(w[5],w[4],w[6],w[7],y)
             	ds2=DS_KE(w[5],w[4],pos2,w[9],y)
             	ds3=DS_KE(w[5],w[4],pos3,w[11],y)
             	ds4=DS_KE(w[5],w[4],pos4,w[13],y)
             	poli=0 //w[0]+w[1]*y+w[2]*y^2
             	integral+=normG*exp(-(start)^2/(2*sigma^2))*(ds1+ds2+ds3+ds4+poli)*delta
		start+=delta 	
	while (start <= stop)
	integral+=w[0]+w[1]*x+w[2]*x^2
	return integral
end
	
	
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Macro fit_singlepeak_KE(waveYStr,hold)
	String waveYStr
	Prompt waveYStr, "enter wave",popup,WaveList("*_FPR_SY",";","WIN:")
	variable hold=1
	Prompt hold, "start coefficient:", popup "guess; use existing; use existing and hold L G alpha"
	variable n=strlen(waveYStr)-1
	String waveXstr=waveYStr[0,n-1]+"X"	
	single_KE($waveYStr,$waveXStr,hold)
end

Function single_KE(waveY,waveX,hold)
	wave waveY,waveX
	variable hold
	wavestats/Q waveY
	variable np=numpnts(waveY)-1
	variable theta = (waveY[0]-waveY[np])/(waveX[0]-waveX[np])
	variable peakpos=waveX[V_maxloc]
	if (numpnts(waveY) > -1 && numpnts(waveX)>-1 )
		if (hold == 1) 
			make/n=8/o W_coef=0
			W_coef[0]=V_min
			W_coef[1]=0
			W_coef[2]=0
			W_coef[3]=1.6 //1
			W_coef[4]=0.35
			W_coef[5]=0.1
			W_coef[6]=peakpos
			W_coef[7]=V_max
			FuncFit/H="00111100"/Q onepeakCDS_KE W_coef waveY /X=waveX /D
			FuncFit/H="00101000" onepeakCDS_KE W_coef waveY /X=waveX /D 
		endif
		if (hold == 2)
			FuncFit/H="00101000"/Q onepeakCDS_KE W_coef waveY /X=waveX /D 
		endif	
		if (hold == 3)
			FuncFit/H="00011100"/Q  onepeakCDS_KE W_coef  waveY /X=waveX /D /R
		endif
		ModifyGraph mode(waveY)=3;DelayUpdate
		ModifyGraph marker(waveY)=41;DelayUpdate
		ModifyGraph rgb(waveY)=(0,34816,52224)
		
	endif
end 

Macro fit_doublet_KE(waveYStr,hold)
	String waveYStr
	Prompt waveYStr, "enter wave",popup,WaveList("*_FPR_SY",";","WIN:")
	variable hold=1
	Prompt hold, "start coefficient:", popup "guess; use existing; use existing and hold L G alpha"
	variable n=strlen(waveYStr)-1
	String waveXstr=waveYStr[0,n-1]+"X"	
	doublet_KE($waveYStr,$waveXStr,hold)
end

Function doublet_KE(waveY,waveX,hold)
	wave waveY,waveX
	variable hold
	wavestats waveY
	variable np=numpnts(waveY)-1
	variable theta = (waveY[0]-waveY[np])/(waveX[0]-waveX[np])
	variable peakpos=waveX(V_maxloc)
	if (numpnts(waveY) > -1 && numpnts(waveX)>-1 )
		if (hold == 1) 
			make/n=11/o W_coef=0
			W_coef[0]=V_min
			W_coef[1]=0
			W_coef[2]=0
			W_coef[3]=0.25//1
			W_coef[4]=1.6  //3054
			W_coef[5]=0.1  //4487
			W_coef[6]=peakpos
			W_coef[7]=V_max
			W_coef[8]=2.02 //-2.1
			W_coef[9]=V_max*0.1 //*0.8
			W_coef[10]=1.6//*0.8
			print V_min, peakpos
			//FuncFit/H="0011110000"/Q twopeakCDS_KE W_coef waveY /X=waveX /D
			FuncFit/H="00111100001" twopeakCDS_KE W_coef waveY /X=waveX /D
			FuncFit/H="00110000000" twopeakCDS_KE W_coef waveY /X=waveX /D  
		endif
		if (hold == 2)
			FuncFit/H="00011100000" twopeakCDS_KE W_coef waveY /X=waveX /D
		endif	
		if (hold == 3)
			FuncFit/H="0001110000"/Q twopeakCDS_KE W_coef  waveY /X=waveX /D /R
		endif
		//ModifyGraph mode(waveY)=3;DelayUpdate
		//ModifyGraph marker(waveY)=41;DelayUpdate
		//ModifyGraph rgb(waveY)=(0,34816,52224)
		
	endif
end 

Macro fit_Ti_KE(waveYStr,npeaks,hold)
	String waveYStr
	Prompt waveYStr, "enter wave",popup,WaveList("*_Y",";","WIN:")
	variable npeaks=3
	Prompt npeaks, "number of peaks:", popup "3; 4"
	variable hold=1
	Prompt hold, "start coefficient:", popup "guess; use existing; use existing and hold L G alpha"
	variable n=strlen(waveYStr)-2
	String waveXstr=waveYStr[0,n]+"X"	
	print waveXstr, waveYStr
	npeaks+=2
	if (npeaks == 3)
		threepeaks_KE($waveYStr,$waveXStr,hold)
	else 
		fivepeaks_KE($waveYStr,$waveXStr,hold)
	endif
	if (duplicate_W_coef()==0)
			String c1Str=waveYStr+"_C1"
			String c2Str=waveYStr+"_C2"
			String c3Str=waveYStr+"_C3"
			String c4Str=waveYStr+"_C4"
			String c5Str=waveYStr+"_C5"
			String dataStr =waveYStr+"_data"
			String fitStr =waveYStr+"_fit"
			String bkgStr =waveYStr+"_bkg"
			duplicate/o $waveYStr $c1Str, $c2Str, $c3Str, $c4Str, $c5Str,$dataStr, $fitStr, $bkgStr
			duplicate/o W_coef W_tmp
			if (npeaks == 3)
				$fitStr= threepeakCDS_KE(W_tmp,$waveXStr)
			else if (npeaks == 4)
				$fitStr= fourpeakCDS_KE(W_tmp,$waveXStr)
			else
				$fitStr= fivepeakCDS_KE(W_tmp,$waveXStr)
			endif
			W_tmp[9]=0
			W_tmp[11]=0
			$c1Str= threepeakCDS_KE(W_tmp,$waveXStr)
			duplicate/o W_coef W_tmp
			W_tmp[7]=0
			W_tmp[11]=0
			$c2Str= threepeakCDS_KE(W_tmp,$waveXStr)
			duplicate/o W_coef W_tmp
			W_tmp[9]=0
			W_tmp[7]=0
			$c3Str= threepeakCDS_KE(W_tmp,$waveXStr)
			W_tmp[11]=0
			$bkgStr= threepeakCDS_KE(W_tmp,$waveXStr)
			$c1Str-=$bkgStr
			$c2Str-=$bkgStr
			$c3Str-=$bkgStr  
			$fitStr-=$bkgStr
			$dataStr-=$bkgStr
			display $dataStr vs $waveXStr
			append  $c1Str $c2Str $c3Str $fitStr vs $waveXStr
			if (npeaks < 5)
				$c4Str= fourpeakCDS_KE(W_tmp,$waveXStr)
				append  $c4Str vs $waveXStr
				modifyGraph lsize($c4Str)=2;DelayUpdate
				ModifyGraph rgb($c4Str)=(52224,0,41728)
			endif
			
			if (npeaks == 5)
				$c5Str= fivepeakCDS_KE(W_tmp,$waveXStr)
				append  $c5Str vs $waveXStr
				modifyGraph lsize($c5Str)=2;DelayUpdate
				ModifyGraph rgb($c5Str)=(52224,0,41728)
			endif
			
			ModifyGraph rgb($fitStr)=(0,0,0)
			ModifyGraph mode($dataStr)=3,marker($dataStr)=8;DelayUpdate
			ModifyGraph rgb($dataStr)=(0,0,65280)
			ModifyGraph rgb($c2Str)=(0,43520,65280)
			ModifyGraph rgb($c3Str)=(0,52224,0)
			ModifyGraph lsize($c1Str)=2,lsize($c2Str)=2,lsize($c3Str)=2
			Label left "\\F'Arial'\\Z16intensity (a.u.)";DelayUpdate
			Label bottom "\\F'Arial'\\Z16kinetic energy (eV)"
			ModifyGraph nticks(left)=0
			ModifyGraph fSize(bottom)=12,font(bottom)="Arial"
		endif
end

Function threepeaks_KE(waveY,waveX,hold)
	wave waveY,waveX
	variable hold
	wavestats waveY
	variable np=numpnts(waveY)-1
	variable theta = (waveY[0]-waveY[np])/(waveX[0]-waveX[np])
	variable peakpos=waveX(V_maxloc)
	if (numpnts(waveY) > -1 && numpnts(waveX)>-1 )
		if (hold == 1) 
			make/n=12/o W_coef=0
			W_coef[0]=V_min
			W_coef[1]=0
			W_coef[2]=0
			W_coef[3]=0.25//1
			W_coef[4]=1.5  //3054
			W_coef[5]=0.08  //4487
			W_coef[6]=peakpos
			W_coef[7]=V_max
			W_coef[8]=1.1 //0.921 //-2.1
			W_coef[9]=V_max*0.2//*0.8
			W_coef[10]=2.35585//*0.8
			W_coef[11]=V_max*0.5//*0.8
			print V_min, peakpos
			//FuncFit/H="0011110000"/Q twopeakCDS_KE W_coef waveY /X=waveX /D
			FuncFit/H="011100001010" threepeakCDS_KE W_coef waveY /X=waveX /D
			FuncFit/H="011100000000" threepeakCDS_KE W_coef waveY /X=waveX /D
			//FuncFit/H="000100000000" threepeakCDS_KE W_coef waveY /X=waveX /D  
		endif
		if (hold == 2)
			FuncFit/H="011100001010" threepeakCDS_KE W_coef waveY /X=waveX /D
			
		endif	
		if (hold == 3)
			FuncFit/H="0001110000"/Q twopeakCDS_KE W_coef  waveY /X=waveX /D /R
		endif
		//ModifyGraph mode(waveY)=3;DelayUpdate
		//ModifyGraph marker(waveY)=41;DelayUpdate
		//ModifyGraph rgb(waveY)=(0,34816,52224)
		
	endif
end 

Function fourpeaks_KE(waveY,waveX,hold)
	wave waveY,waveX
	variable hold
	wavestats waveY
	variable np=numpnts(waveY)-1
	variable theta = (waveY[0]-waveY[np])/(waveX[0]-waveX[np])
	variable peakpos=waveX(V_maxloc)
	if (numpnts(waveY) > -1 && numpnts(waveX)>-1 )
		if (hold == 1) 
			make/n=14/o W_coef=0
			//W_coef[0]=V_min
			//W_coef[1]=0
			//W_coef[2]=0
			//W_coef[3]=0.25//1
			//W_coef[4]=1.07  //3054
			//W_coef[5]=0.073  //4487
			//W_coef[6]=123.636
			//W_coef[7]=V_max*0.3
			//W_coef[8]=0.866 //0.921 //-2.1
			//W_coef[9]=V_max*0.3//*0.8
			//W_coef[10]=2.10//*0.8
			//W_coef[11]=V_max*0.6//*0.8
			//W_coef[12]=3.7//*0.8
			//W_coef[13]=V_max*0.2//*0.8
			
			print V_min, peakpos
			//FuncFit/H="0011110000"/Q twopeakCDS_KE W_coef waveY /X=waveX /D
			  W_coef={1014.1,0,0,0.25,1.331,0.10409,123.46,14239,1.2532,14019,2.5197,26460,3.6465,15212}
			FuncFit/H="01111100101010" fourpeakCDS_KE W_coef waveY /X=waveX /D
			//FuncFit/H="01110000000000" fourpeakCDS_KE W_coef waveY /X=waveX /D
			//FuncFit/H="011100000000" threepeakCDS_KE W_coef waveY /X=waveX /D
			//FuncFit/H="000100000000" threepeakCDS_KE W_coef waveY /X=waveX /D  
		endif
		if (hold == 2)
			//FuncFit/H="011100001010" threepeakCDS_KE W_coef waveY /X=waveX /D
			    W_coef={1014.1,0,0,0.25,1.331,0.10409,123.46,14239,1.2532,14019,2.5197,26460,3.6465,15212}
			    FuncFit/H="01110000000000" fourpeakCDS_KE W_coef waveY /X=waveX /D
		endif	
		if (hold == 3)
			//FuncFit/H="0001110000"/Q twopeakCDS_KE W_coef  waveY /X=waveX /D /R
		endif
		//ModifyGraph mode(waveY)=3;DelayUpdate
		//ModifyGraph marker(waveY)=41;DelayUpdate
		//ModifyGraph rgb(waveY)=(0,34816,52224)
		
	endif
end 

Function fivepeaks_KE(waveY,waveX,hold)
	wave waveY,waveX
	variable hold
	wavestats waveY
	variable np=numpnts(waveY)-1
	variable theta = (waveY[0]-waveY[np])/(waveX[0]-waveX[np])
	variable peakpos=waveX(V_maxloc)
	if (numpnts(waveY) > -1 && numpnts(waveX)>-1 )
		if (hold == 1) 
			make/n=14/o W_coef=0
			//W_coef[0]=V_min
			//W_coef[1]=0
			//W_coef[2]=0
			//W_coef[3]=0.25//1
			//W_coef[4]=1.07  //3054
			//W_coef[5]=0.073  //4487
			//W_coef[6]=123.636
			//W_coef[7]=V_max*0.3
			//W_coef[8]=0.866 //0.921 //-2.1
			//W_coef[9]=V_max*0.3//*0.8
			//W_coef[10]=2.10//*0.8
			//W_coef[11]=V_max*0.6//*0.8
			//W_coef[12]=3.7//*0.8
			//W_coef[13]=V_max*0.2//*0.8
			
			print V_min, peakpos
			//FuncFit/H="0011110000"/Q twopeakCDS_KE W_coef waveY /X=waveX /D
			  W_coef={1014.1,0,0,0.25,1.331,0.10409,123.46,14239,1.2532,14019,2.5197,26460,3.6465,15212}
			FuncFit/H="01111100101010" fourpeakCDS_KE W_coef waveY /X=waveX /D
			//FuncFit/H="01110000000000" fourpeakCDS_KE W_coef waveY /X=waveX /D
			//FuncFit/H="011100000000" threepeakCDS_KE W_coef waveY /X=waveX /D
			//FuncFit/H="000100000000" threepeakCDS_KE W_coef waveY /X=waveX /D  
		endif
		if (hold == 2)
			//FuncFit/H="011100001010" threepeakCDS_KE W_coef waveY /X=waveX /D
			    W_coef={1014.1,0,0,0.25,1.331,0.10409,123.46,14239,1.2532,14019,2.2,26460,2.9,7000,3.6465,7212}
			    FuncFit/H="0111000000000000" fivepeakCDS_KE W_coef waveY /X=waveX /D
		endif	
		if (hold == 3)
			//FuncFit/H="0001110000"/Q twopeakCDS_KE W_coef  waveY /X=waveX /D /R
		endif
		//ModifyGraph mode(waveY)=3;DelayUpdate
		//ModifyGraph marker(waveY)=41;DelayUpdate
		//ModifyGraph rgb(waveY)=(0,34816,52224)
		
	endif
end 


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// three peaks; same Gauss, Lor, alpha
Function threepeakCDS_BE(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) = sum(i) sum (DS (x(i)-y) * exp () * dy)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 12
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 12
	//CurveFitDialog/ w[0] = k0
	//CurveFitDialog/ w[1] = k1
	//CurveFitDialog/ w[2] = k2
	//CurveFitDialog/ w[3] = gaussfwhm
	//CurveFitDialog/ w[4] = gamma
	//CurveFitDialog/ w[5] = alpha
	//CurveFitDialog/ w[6] = pos1
	//CurveFitDialog/ w[7] = A1
	//CurveFitDialog/ w[8] = deltapos2
	//CurveFitDialog/ w[9] = A2
	//CurveFitDialog/ w[10] = deltapos3
	//CurveFitDialog/ w[11] = A3
	
	Variable i=0
	Variable j=0
	Variable sigma = w[3]/(2*sqrt(2*ln(2)))  
	Variable nsteps = 80
	Variable integral = 0
	Variable start = -4*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable ds1,ds2,ds3,poli,y
	Variable normG = 1/sqrt(2*Pi*sigma^2)
	Variable pos2 = w[6]+w[8]
	Variable pos3 = w[6]+w[10]
	do
		y=x-start
		ds1=DS_BE(w[5],w[4],w[6],w[7],y)
             ds2=DS_BE(w[5],w[4],pos2,w[9],y)
             ds3=DS_BE(w[5],w[4],pos3,w[11],y)
             poli=w[0]+w[1]*y+w[2]*y^2
             integral+=normG*exp(-(start)^2/(2*sigma^2))*(ds1+ds2+ds3+poli)*delta
		start+=delta 	
	while (start <= stop)
	return integral
end

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// four peaks;  same gauss, lor, alpha
Function fourpeakCDS_BE(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) = sum(i) sum (DS (x(i)-y) * exp () * dy)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 14
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 14
	//CurveFitDialog/ w[0] = k0
	//CurveFitDialog/ w[1] = k1
	//CurveFitDialog/ w[2] = k2
	//CurveFitDialog/ w[3] = gaussfwhm
	//CurveFitDialog/ w[4] = gamma
	//CurveFitDialog/ w[5] = alpha
	//CurveFitDialog/ w[6] = pos1
	//CurveFitDialog/ w[7] = A1
	//CurveFitDialog/ w[8] = pos2
	//CurveFitDialog/ w[9] = A2
	//CurveFitDialog/ w[10] = pos3
	//CurveFitDialog/ w[11] = A3
	//CurveFitDialog/ w[12] = pos4
	//CurveFitDialog/ w[13] = A4
	
	Variable i=0
	Variable j=0
	Variable sigma = w[3]/(2*sqrt(2*ln(2)))  
	Variable nsteps = 60
	Variable integral = 0
	Variable start = -4*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable ds1,ds2,ds3,ds4,poli,y
	Variable normG = 1/sqrt(2*Pi*sigma^2)

	do
		y=x-start
		ds1=DS_BE(w[5],w[4],w[6],w[7],y)
             ds2=DS_BE(w[5],w[4],w[8],w[9],y)
             ds3=DS_BE(w[5],w[4],w[10],w[11],y)
             ds4=DS_BE(w[5],w[4],w[12],w[13],y)
             poli=w[0]+w[1]*y+w[2]*y^2
             integral+=normG*exp(-(start)^2/(2*sigma^2))*(ds1+ds2+ds3+ds4+poli)*delta
		start+=delta 	
	while (start <= stop)
	return integral
end

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// four peaks;  same gauss, lor, alpha
Function fivepeakCDS_BE(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) = sum(i) sum (DS (x(i)-y) * exp () * dy)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 14
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 14
	//CurveFitDialog/ w[0] = k0
	//CurveFitDialog/ w[1] = k1
	//CurveFitDialog/ w[2] = k2
	//CurveFitDialog/ w[3] = gaussfwhm
	//CurveFitDialog/ w[4] = gamma
	//CurveFitDialog/ w[5] = alpha
	//CurveFitDialog/ w[6] = pos1
	//CurveFitDialog/ w[7] = A1
	//CurveFitDialog/ w[8] = pos2
	//CurveFitDialog/ w[9] = A2
	//CurveFitDialog/ w[10] = pos3
	//CurveFitDialog/ w[11] = A3
	//CurveFitDialog/ w[12] = pos4
	//CurveFitDialog/ w[13] = A4
	
	Variable i=0
	Variable j=0
	Variable sigma = w[3]/(2*sqrt(2*ln(2)))  
	Variable nsteps = 60
	Variable integral = 0
	Variable start = -4*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable ds1,ds2,ds3,ds4,ds5,poli,y
	Variable normG = 1/sqrt(2*Pi*sigma^2)

	do
		y=x-start
		ds1=DS_BE(w[5],w[4],w[6],w[7],y)
             ds2=DS_BE(w[5],w[4],w[8],w[9],y)
             ds3=DS_BE(w[5],w[4],w[10],w[11],y)
             ds4=DS_BE(w[5],w[4],w[12],w[13],y)
             ds5=DS_BE(w[5],w[4],w[14],w[15],y)
             poli=w[0]+w[1]*y+w[2]*y^2
             integral+=normG*exp(-(start)^2/(2*sigma^2))*(ds1+ds2+ds3+ds4+ds5+poli)*delta
		start+=delta 	
	while (start <= stop)
	return integral
end






// forur peaks; different gauss, lor, alpha; delta peak
Function fourpeakCDS_Rel(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) = sum(i) sum (DS (x(i)-y) * exp () * dy)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 14
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 14
	//CurveFitDialog/ w[0] = k0
	//CurveFitDialog/ w[1] = k1
	//CurveFitDialog/ w[2] = k2
	//CurveFitDialog/ w[3] = gaussfwhm
	//CurveFitDialog/ w[4] = gamma
	//CurveFitDialog/ w[5] = alpha
	//CurveFitDialog/ w[6] = pos1
	//CurveFitDialog/ w[7] = A1
	//CurveFitDialog/ w[8] = del2
	//CurveFitDialog/ w[9] = A2
	//CurveFitDialog/ w[10] = del3
	//CurveFitDialog/ w[11] = A3
	//CurveFitDialog/ w[12] = pos4
	//CurveFitDialog/ w[13] =A4
	
	Variable i=0
	Variable j=0
	Variable sigma = w[3]/(2*sqrt(2*ln(2)))  
	Variable nsteps = 50
	Variable integral = 0
	Variable start = -4*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable ds1,ds2,ds3,ds4,poli,y
	Variable normG = 1/sqrt(2*Pi*sigma^2)

	do
		y=x-start
		ds1=DS_BE(w[5],w[4],w[6],w[7],y)
             ds2=DS_BE(w[5],w[4],w[6]+w[8],w[9],y)
             ds3=DS_BE(w[5],w[4],w[6]+w[10],w[11],y)
             ds4=DS_BE(w[5],w[4],w[6]+w[12],w[13],y)
             poli=w[0]+w[1]*y+w[2]*y^2
             integral+=normG*exp(-(start)^2/(2*sigma^2))*(ds1+ds2+ds3+ds4+poli)*delta
		start+=delta 	
	while (start <= stop)
	return integral
end

// fit doublet, two componentes(Bulk, surface); different alpha, lor for each peak 
Function fourpeakCDS_doublet(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) = sum(i) sum (DS (x(i)-y) * exp () * dy)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 14

	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 14
	//CurveFitDialog/ w[0] = k0
	//CurveFitDialog/ w[1] = k1
	//CurveFitDialog/ w[2] = k2
	//CurveFitDialog/ w[3] = gaussfwhm
	//CurveFitDialog/ w[4] = gamma1
	//CurveFitDialog/ w[5] = alpha1
	//CurveFitDialog/ w[6] = gamma2
	//CurveFitDialog/ w[7] = alpha2
	//CurveFitDialog/ w[8] = scls
	//CurveFitDialog/ w[9] = pos1
	//CurveFitDialog/ w[10] = pos2
	//CurveFitDialog/ w[11] = A1
	//CurveFitDialog/ w[12] = A2
	//CurveFitDialog/ w[13] = A3
	
	// preparing gaussian convolution
	Variable i=0
	Variable j=0
	Variable sigma = w[3]/(2*sqrt(2*ln(2)))  
	Variable nsteps = 50
	Variable integral = 0
	Variable start = -4*sigma
	Variable stop = -start
	Variable delta = 2* abs (start) / nsteps
	Variable ds1,ds2,ds3,ds4,poli,y
	Variable normG = 1/sqrt(2*Pi*sigma^2)
        // convolution loop
	do
		y=x-start
		// DS(alpha,gammaL,pos,Ampli)
		ds1=DS_BE(w[5],w[4],w[9],w[11],y) 				// pos 1 bulk
             	ds2=DS_BE(w[7],w[6],w[9]+w[8],w[11]*w[13],y) 		// pos 1 surface
             	ds3=DS_BE(w[5],w[4],w[10],w[12],y)				// pos 2 bulk
             	ds4=DS_BE(w[7],w[6],w[10]+w[8],w[12]*w[13],y)		// pos 2 surface
             	poli=w[0]+w[1]*y+w[2]*y^2
             	integral+=normG*exp(-(start)^2/(2*sigma^2))*(ds1+ds2+ds3+ds4+poli)*delta
		start+=delta 	
	while (start <= stop)
	return integral
end

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// fit Au 4f doublet
Function Au4f_XY_KE(waveYStr)
	String waveYStr
		String waveXStr
	variable n=strlen(waveYStr)-1
	wavestats/Q $waveYStr
	variable np=numpnts($waveYStr)-1
	//variable theta = ($waveYStr[0]-$waveYStr[np])/($waveXStr[0]-$waveXStr[np])
	//variable peakpos=$waveXStr[V_maxloc]
	 Make/n=10/o W_coef={1064.8,-34.664,0.25,0.25485,0.54055,0.034444,77.281,9869.8,-3.5592,6885.3}
	if (exists(waveYStr) ==1 && exists(waveXStr)==1)
		//if (hold == 1) 
		//	make/n=10/o W_coef=0
		//	W_coef[0]=V_min
		//	W_coef[1]=0
		//	W_coef[2]=0
		//	W_coef[3]=0.25
		//	W_coef[4]=0.8
		//	W_coef[5]=0.0
		//	W_coef[6]=peakpos
		//	W_coef[7]=V_max
		//	W_coef[8]=-3.6
		//	W_coef[9]=V_max*0.8
			FuncFit/H="0011110000"/Q twopeakCDS_KE W_coef $waveYStr /X=$waveXStr /D
			FuncFit/H="0001100000" twopeakCDS_KE W_coef $waveYStr /X=$waveXStr /D /R
		//endif
	endif
end

// fit Au 4f doublet
Macro fit_Au4f_XY_KE(waveYStr,waveXStr,hold)
	String waveYStr
	Prompt waveYStr, "enter Y wave",popup,WaveList("*",";","WIN:")
	String waveXStr
	Prompt waveXStr, "enter X wave",popup,WaveList("*",";","WIN:")
	variable hold=1
	Prompt hold, "start coefficient:", popup "guess; use existing; use existing and hold L G alpha"
	variable n=strlen(waveYStr)-1
	//String waveXstr=waveYStr[0,n-1]+"X"	
	wavestats/Q $waveYStr
	variable np=numpnts($waveYStr)-1
	variable theta = ($waveYStr[0]-$waveYStr[np])/($waveXStr[0]-$waveXStr[np])
	variable peakpos=$waveXStr[V_maxloc]
	if (exists(waveYStr) ==1 && exists(waveXStr)==1)
		if (hold == 1) 
			make/n=13/o W_coef=0
			W_coef[0]=V_min
			W_coef[1]=0
			W_coef[2]=0
			
			W_coef[3]=0.35
			W_coef[4]=0.32
			W_coef[5]=0.01
			
			W_coef[6]=96 //77.5 //76.8 dp4//peakpos
			W_coef[7]=V_max
			W_coef[8]=3.7  //0.41224 //3.6
			W_coef[9]=V_max*0.8
			
			//W_coef[10]=-3.761
			//W_coef[11]=0.702
			//W_coef[12]=0.4
			//FuncFit/H="0011110000"/Q twopeakCDS_KE W_coef $waveYStr /X=$waveXStr /D
			//FuncFit/H="0001000000" twopeakCDS_KE W_coef $waveYStr /X=$waveXStr /D /R
			FuncFit/H="0000110000000" twopeakCDS_KE W_coef $waveYStr[pcsr(A),pcsr(B)] /X=$waveXStr /D /R
		endif
		if (hold == 2)
			//FuncFit/H="0000000000"/Q twopeakCDS_KE W_coef  $waveYStr /X=$waveXStr /D /R
			FuncFit/H="0001111010111" fourpeakCDS_KE W_coef $waveYStr[pcsr(A),pcsr(B)] /X=$waveXStr /D /R
			//FuncFit/H="0001000000000" fourpeakCDS_KE W_coef $waveYStr[pcsr(A),pcsr(B)] /X=$waveXStr /D /R
		endif	
		if (hold == 3)
			FuncFit/H="0001110000"/Q twopeakCDS_KE W_coef  $waveYStr /X=$waveXStr /D /R
		endif
		ModifyGraph mode($waveYStr)=3;DelayUpdate
		ModifyGraph marker($waveYStr)=41;DelayUpdate
		ModifyGraph rgb($waveYStr)=(0,34816,52224)
		String tmpStr = "Res_"+waveYStr
		ModifyGraph rgb($tmpStr)=(0,0,0)
		tmpStr = "fit_"+waveYStr
		ModifyGraph rgb($tmpStr)=(0,0,0)
		if (duplicate_W_coef()==0)
			String c1Str=waveYStr+"_C1"
			String c2Str=waveYStr+"_C2"
			String bkgStr=waveYStr+"_bkg"
			String dataStr =waveYStr+"_data"
			String fitStr =waveYStr+"_fit"
			duplicate/o $waveYStr $c1Str, $c2Str, $bkgStr, $dataStr, $fitStr
			duplicate/o W_coef W_tmp
			$fitStr= fourpeakCDS_KE(W_tmp,$waveXStr)
			W_tmp[9]=0
			$c1Str= fourpeakCDS_KE(W_tmp,$waveXStr)
			duplicate/o W_coef W_tmp
			W_tmp[7]=0
			$c2Str= fourpeakCDS_KE(W_tmp,$waveXStr)
			W_tmp[9]=0
			$bkgStr= fourpeakCDS_KE(W_tmp,$waveXStr) 
			$c1Str-=$bkgStr
			$c2Str-=$bkgStr
			$fitStr-=$bkgStr
			$dataStr-=$bkgStr
			display $dataStr vs $waveXStr
			append  $c1Str $c2Str $fitStr vs $waveXStr
			ModifyGraph rgb($fitStr)=(0,0,0)
			ModifyGraph mode($dataStr)=3,marker($dataStr)=8;DelayUpdate
			ModifyGraph rgb($dataStr)=(0,0,65280)
			ModifyGraph rgb($c2Str)=(0,43520,65280)
			ModifyGraph lsize($c1Str)=2,lsize($c2Str)=2
			Label left "\\F'Arial'\\Z16intensity (a.u.)";DelayUpdate
			Label bottom "\\F'Arial'\\Z16kinetic energy (eV)"
			ModifyGraph nticks(left)=0
			ModifyGraph fSize(bottom)=12,font(bottom)="Arial"
		endif
	endif
end

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// fit Au 4f doublet
Macro fit_Au4f72_XY(waveYStr,waveXStr,hold)
	String waveYStr
	Prompt waveYStr, "enter Y wave",popup,WaveList("*",";","WIN:")
	String waveXStr
	Prompt waveXStr, "enter X wave",popup,WaveList("*",";","WIN:")
	variable hold=1
	Prompt hold, "start coefficient:", popup "guess; use existing; use existing and hold L G alpha"
	//variable n=strlen(waveYStr)-1
	//String waveXstr=waveYStr[0,n-1]+"X"	
	wavestats/Q $waveYStr
	variable np=numpnts($waveYStr)-1
	variable theta = ($waveYStr[0]-$waveYStr[np])/($waveXStr[0]-$waveXStr[np])
	variable peakpos=$waveXStr[V_maxloc]
	if (exists(waveYStr) ==1 && exists(waveXStr)==1)
		if (hold == 1) 
			make/n=10/o W_coef=0
			W_coef[0]=V_min
			W_coef[1]=0
			W_coef[2]=0
			W_coef[3]=0.22
			W_coef[4]=0.330
			W_coef[5]=0.072
			W_coef[6]=peakpos
			W_coef[7]=V_max
			W_coef[8]=-3.6
			W_coef[9]=V_max*0.8
			FuncFit/H="0111100000"/Q twopeakCDS_KE W_coef $waveYStr /X=$waveXStr /D
			FuncFit/H="0001000000" twopeakCDS_KE W_coef $waveYStr /X=$waveXStr /D /R
			//FuncFit/H="0010000010" twopeakCDS_KE W_coef $waveYStr /X=$waveXStr /D /R
		endif
		if (hold == 2)
			FuncFit/H="0000000000"/Q twopeakCDS_KE W_coef  $waveYStr /X=$waveXStr /D /R
		endif	
		if (hold == 3)
			FuncFit/H="0001110000"/Q twopeakCDS_KE W_coef  $waveYStr /X=$waveXStr /D /R
		endif
		ModifyGraph mode($waveYStr)=3;DelayUpdate
		ModifyGraph marker($waveYStr)=41;DelayUpdate
		ModifyGraph rgb($waveYStr)=(0,34816,52224)
		String tmpStr = "Res_"+waveYStr
		ModifyGraph rgb($tmpStr)=(0,0,0)
		tmpStr = "fit_"+waveYStr
		ModifyGraph rgb($tmpStr)=(0,0,0)
		duplicate_W_coef()
	endif
end


Macro fit_Au4f_2c_XY(waveYStr,hold)
	String waveYStr
	Prompt waveYStr, "enter wave",popup,WaveList("*_FPR_SY",";","WIN:")
	variable hold=1
	Prompt hold, "hold Gaussian and lorenzian width (1)", popup "no;yes"
	variable n=strlen(waveYStr)-1
	String waveXstr=waveYStr[0,n-1]+"X"	
	wavestats/Q $waveYStr
	variable np=numpnts($waveYStr)-1
	variable theta = ($waveYStr[0]-$waveYStr[np])/($waveXStr[0]-$waveXStr[np])
	variable peakpos=$waveXStr[V_maxloc]
	if (exists(waveYStr) ==1 && exists(waveXStr)==1)
		if (hold == 2) 
			FuncFit/H="00001111000000" fourpeakCDS_doublet W_coef $waveYStr /X=$waveXStr /D /R
		else 
			make/n=14/o W_coef=0
			W_coef[0]=V_min
			W_coef[1]=theta
			W_coef[2]=0
			W_coef[3]=0.26		// guess for instumental resolution
			W_coef[4]=0.23  		// lor 1
			W_coef[5]=0.004 		// alpha 1
			W_coef[6]=0.23		// lor 2
			W_coef[7]=0.004		// apha 2
			W_coef[8]=-0.19		// SCLS
			W_coef[9]=peakpos+0.05
			W_coef[10]=W_coef[9]+3.65
			W_coef[11]=V_max
			W_coef[12]=V_max*0.8
			W_coef[13]=0.1
			FuncFit/H="00111111111000"/Q fourpeakCDS_doublet W_coef $waveYStr /X=$waveXStr /D
			FuncFit/H="00011111100000" fourpeakCDS_doublet W_coef $waveYStr /X=$waveXStr /D
			FuncFit/H="00011010000000" fourpeakCDS_doublet W_coef $waveYStr /X=$waveXStr /D /R
		endif
		duplicate_W_coef()
	endif	
end

function duplicate_W_coef()
	String coefStr
	Prompt coefStr, "enter name for coef wave:"
	doPrompt "duplicate W_coef:", coefStr
	if (V_flag == 0)
		duplicate/o W_coef $coefStr
		return 0
	else
		return 1
	endif
end

// fit tungsten 4f; need shirley bkg removal before use
Macro fit_W4f_bulk_and_surface_XY(waveYStr,hold)
	String waveYStr
	Prompt waveYStr, "enter wave",popup,WaveList("*_FPR_SY",";","WIN:")
	variable hold=1
	Prompt hold, "hold Gaussian and lorenzian width (1)", popup "no;yes"
	variable n=strlen(waveYStr)-1
	String waveXstr=waveYStr[0,n-1]+"X"	
	wavestats/Q $waveYStr
	variable np=numpnts($waveYStr)-1
	variable theta = ($waveYStr[0]-$waveYStr[np])/($waveXStr[0]-$waveXStr[np])
	variable peakpos=$waveXStr[V_maxloc]
	if (exists(waveYStr) ==1 && exists(waveXStr)==1)
		if (hold == 2) 
			FuncFit/H="00001111000000" fourpeakCDS_doublet W_coef $waveYStr /X=$waveXStr /D /R
		else 
			make/n=14/o W_coef=0
			W_coef[0]=V_min
			W_coef[1]=theta
			W_coef[2]=0
			W_coef[3]=0.25		// guess for instumental resolution
			W_coef[4]=0.060  		// lor 1
			W_coef[5]=0.035 		// alpha 1
			W_coef[6]=0.084		// lor 2
			W_coef[7]=0.063		// apha 2
			W_coef[8]=-0.321		// SCLS
			W_coef[9]=peakpos+0.05
			W_coef[10]=W_coef[9]+2.2
			W_coef[11]=V_max
			W_coef[12]=V_max*0.8
			W_coef[13]=0.07
			FuncFit/H="00111111111000"/Q fourpeakCDS_doublet W_coef $waveYStr /X=$waveXStr /D
			FuncFit/H="00001111100000" fourpeakCDS_doublet W_coef $waveYStr /X=$waveXStr /D
			FuncFit/H="00001111000000" fourpeakCDS_doublet W_coef $waveYStr /X=$waveXStr /D /R
		endif
		duplicate_W_coef()
	endif	
end

Macro fit_W4f_oxide_XY(waveYStr,hold)
	String waveYStr
	Prompt waveYStr, "enter wave",popup,WaveList("*_FPR_SY",";","WIN:")
	variable hold=1
	Prompt hold, "hold Gaussian and lorenzian width (1)", popup "no;yes"
	variable n=strlen(waveYStr)-1
	String waveXstr=waveYStr[0,n-1]+"X"	
	wavestats/Q $waveYStr
	variable np=numpnts($waveYStr)-1
	variable theta = ($waveYStr[0]-$waveYStr[np])/($waveXStr[0]-$waveXStr[np])
	variable peakpos=$waveXStr[V_maxloc]
	if (exists(waveYStr) ==1 && exists(waveXStr)==1)
		if (hold == 2) 
			FuncFit/H="00001010000000" fourpeakCDS_doublet W_coef $waveYStr /X=$waveXStr /D /R
		else 
			make/n=14/o W_coef=0
			W_coef[0]=V_min
			W_coef[1]=theta
			W_coef[2]=0
			W_coef[3]=0.25		// guess for instumental resolution
			W_coef[4]=0.060  		// lor 1
			W_coef[5]=0.035 		// alpha 1
			W_coef[6]=0.084		// lor 2
			W_coef[7]=0.063		// apha 2
			W_coef[8]=0.71	// SCLS
			W_coef[9]=peakpos+0.05
			W_coef[10]=W_coef[9]+2.2
			W_coef[11]=V_max
			W_coef[12]=V_max*0.8
			W_coef[13]=0.5
			FuncFit/H="00111111111000"/Q fourpeakCDS_doublet W_coef $waveYStr /X=$waveXStr /D
			FuncFit/H="00001111100000" fourpeakCDS_doublet W_coef $waveYStr /X=$waveXStr /D
			FuncFit/H="00000101000000" fourpeakCDS_doublet W_coef $waveYStr /X=$waveXStr /D /R
		endif
		duplicate_W_coef()
	endif	
end



// fit Wf4
function fit_W4f_bulk_and_surface(wx,wy)
	wave wx,wy
	wavestats/Q wy
			//make/n=14/o W_coef=0
			//W_coef[0]=2810
			//W_coef[1]=-182
			//W_coef[2]=3.5
			//W_coef[3]=0.287		// guess for instumental resolution
			//W_coef[4]=0.060  		// lor 1
			//W_coef[5]=0.035 		// alpha 1
			//W_coef[6]=0.084		// lor 2
			//W_coef[7]=0.063		// apha 2
			//W_coef[8]=-0.321		// SCLS
			//W_coef[9]=wx[V_maxloc]
			//W_coef[10]=W_coef[9]+2.2
			//W_coef[11]=V_max
			//W_coef[12]=V_max*0.66
			//W_coef[13]=0.08
	FuncFit/H="00001111000000" fourpeakCDS_doublet W_coef wy  /X=wx /D /R
end

function twoLayerUptake(w,x)
	//w0 = x0
	//w1=b1
	//w2=b2
	wave  w
	variable x
	if (x<w[0])
	return x*w[1]
	else
	return x*w[1]-(x-w[0])*w[2]
	endif
end

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// histogram functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

function fit_histogram(ctrlname): buttoncontrol	
	String ctrlname
	String waveStr
	Prompt waveStr, "enter wave",popup WaveList("*_HST",";","WIN:")
	variable levels=2
	Prompt levels, "levels", popup "2;3"
	doPrompt "fit histogram", waveStr, levels
	duplicate/o $waveStr fit_histogr
	wave W_coef = root:W_coef
	levels+=1
	if (levels==2)
		redimension/n=6/d W_coef
		Fit_Historgram_2Level(waveStr,W_coef)
	else
		redimension/n=8/d W_coef
		Fit_Historgram_3Level(waveStr,W_coef)
	endif
end

function Fit_Historgram_YW(hstStr) //obsolete
	String hstStr
	Wavestats/Q $hstStr
	wave w=$hstStr
	// init guess	
	Make/n=8/o W_coef={804.67,60,7,70,128,10,44.185,160}
	// fit
	FuncFit/H="00100000"/Q threeLevelStripeSimulation W_coef $hstStr /D   
	// save parameters
		Make/n=10/o W_fit
		W_fit[0]=W_coef[1]*(100-W_coef[6])/100
		W_fit[1]=100-W_coef[1]
		W_fit[2]=W_coef[1]*(W_coef[6])/100
		W_fit[3]=W_coef[3]
		W_fit[4]=W_coef[4]
		W_fit[5]=W_coef[7]
		W_fit[6]=V_chisq	 
	// printout
	print hstStr, ": coef 1 (dark)", W_coef[1]*(100-W_coef[6])/100," at: ", W_coef[3]
	print hstStr, ": 100- coef 1 (neutral)",100-W_coef[1]," at: ", W_coef[4]
	print hstStr, ": coef 1 (bright)",W_coef[1]*(W_coef[6])/100," at: ", W_coef[7]
	print hstStr, ": chisquare",V_chisq	
End

Function twoGaussHistogram(w,x)
	//w0=amplitude (scaling coef)
	//w1=perc_level1 (percentage of histogram level 1)
	//w2=pos_level1
	//w3=pos_level2
	//w4=sigma (gaussian width)
	wave w
	variable x
	variable i1=0,i2=0,iu=0,ilev=0
	variable p1=w[2]
	variable p2=w[3]
	variable s=w[4]
	i1=1/sqrt(2*pi*s^2)*exp(-(x-p1)^2/(2*s^2))  //first gaussian
	i2=1/sqrt(2*pi*s^2)*exp(-(x-p2)^2/(2*s^2)) //second gaussian
	variable a1,a2,a3
	a1=abs(w[1])
	a2=1-w[1]
	return ((a1)*i1+(a2)*i2)*w[0]  
end

Function twoLevelSimulation(w,x)
	//w0=amplitude (scaling coef)
	//w1=perc_level1 (percentage of histogram level 1)
	//w2=perc_boxcar (lateral res/stripe)
	//w3=pos_level1
	//w4=pos_level2
	//w5=sigma (gaussian width)
	wave w
	variable x
	variable i1=0,i2=0,iu=0,ilev=0
	variable p1=w[3]
	variable p2=w[4]
	variable s=w[5]
	i1=1/sqrt(2*pi*s^2)*exp(-(x-p1)^2/(2*s^2))  //first gaussian
	i2=1/sqrt(2*pi*s^2)*exp(-(x-p2)^2/(2*s^2)) //second gaussian
	variable a1,a2,a3
	a1=abs(w[1])
	a2=1-w[1]
	a3=abs(w[2])
	variable pa,pb,level
	if ((a1>a3)&&(a2>a3))
		// normal distibution
		if ((p1<=x)&&(x<=p2))
			iu=1/abs(p2-p1)*abs(gammp(0.5,((p2-x)/s)^2/2)/2+gammp(0.5,((p1-x)/s)^2/2)/2)
		elseif  ((x<p1)&&(x<p2))
			iu=1/abs(p2-p1)*abs(gammp(0.5,((p2-x)/s)^2/2)/2-gammp(0.5,((p1-x)/s)^2/2)/2)
		elseif ((p1<x)&&(p2<x))
        			iu=1/abs(p2-p1)*abs(gammp(0.5,((x-p1)/s)^2/2)/2-gammp(0.5,((x-p2)/s)^2/2)/2)
		endif
		return ((a1-a3)*i1+(a2-a3)*i2+2*a3*iu)*w[0]  
	elseif (a3>=a1)
		level=p2+sign(p1-p2)*abs(p2-p1)*(a1)/a3
		pa=min(p2,level)
		pb=max(p2,level)
		if ((pa<=x)&&(x<=pb))
			iu=1/abs(pb-pa)*abs(gammp(0.5,((pb-x)/s)^2/2)/2+gammp(0.5,((pa-x)/s)^2/2)/2)
		elseif  ((x<pa)&&(x<pb))
			iu=1/abs(pb-pa)*abs(gammp(0.5,((pb-x)/s)^2/2)/2-gammp(0.5,((pa-x)/s)^2/2)/2)
		elseif ((pa<x)&&(pb<x))
        			iu=1/abs(pb-pa)*abs(gammp(0.5,((x-pa)/s)^2/2)/2-gammp(0.5,((x-pb)/s)^2/2)/2)
		endif
		//return ((a2-a3)*i2+(a1+a3)*iu)*w[0]  //uniform distrib
		ilev=1/sqrt(2*pi*s^2)*exp(-(x-level)^2/(2*s^2)) // true boxcar
		return ((a2-a3)*i2+(2*a1)*iu+ilev*(a3-a1))*w[0] 
	elseif (a3>=a2)
		level=p1+sign(p2-p1)*abs(p1-p2)*(a2)/a3
		pa=min(p1,level)
		pb=max(p1,level)
		if ((pa<=x)&&(x<=pb))
			iu=1/abs(pb-pa)*abs(gammp(0.5,((pb-x)/s)^2/2)/2+gammp(0.5,((pa-x)/s)^2/2)/2)
		elseif  ((x<pa)&&(x<pb))
			iu=1/abs(pb-pa)*abs(gammp(0.5,((pb-x)/s)^2/2)/2-gammp(0.5,((pa-x)/s)^2/2)/2)
		elseif ((pa<x)&&(pb<x))
        			iu=1/abs(pb-pa)*abs(gammp(0.5,((x-pa)/s)^2/2)/2-gammp(0.5,((x-pb)/s)^2/2)/2)
		endif
		//return ((a1-a3)*i1+(a2+a3)*iu)*w[0]  //uniform distrib
		ilev=1/sqrt(2*pi*s^2)*exp(-(x-level)^2/(2*s^2)) //true boxcar
		return ((a1-a3)*i1+(2*a2)*iu+ilev*(a3-a2))*w[0] 
		
	endif
End

Function threeGaussHistogram(w,x)
	// weighted sum of three two level hystogram  andone one level
	//w0=scaling
	//w1=pos_level1
	//w2=pos_level2
	//w3=pos_level3
	//w4=perc grey
	//w5=perc dark
	//w6=sigma (gaussian width)
	wave w
	variable x
	wave tmp
	make/n=6/o tmp
	// grey and bright
	tmp[0]=(1-w[5])/2
	tmp[1]=w[4]
	tmp[2]=w[2]
	tmp[3]=w[3]
	tmp[4]=w[6]
	variable itot=twoGaussHistogram(tmp,x)
	// grey and dark
	tmp[0]=(w[5])/2
	tmp[1]=1-w[4]
	tmp[2]=w[1]
	tmp[3]=w[2]
	itot+=twoGaussHistogram(tmp,x)
	// bright and dark
	tmp[0]=(1-w[4])/2
	tmp[1]=w[5]
	tmp[2]=w[1]
	tmp[3]=w[3]
	itot+=twoGaussHistogram(tmp,x)
	// grey
	itot+=0.5*w[4]/sqrt(2*pi*w[6]^2)*exp(-(x-w[2])^2/(2*w[6]^2))
	killwaves tmp
	return w[0]*itot
end

Function threeLevelSimulation(w,x)
	// weighted sum of three two level hystogram  andone one level
	//w0=scaling
	//w1=pos_level1
	//w2=pos_level2
	//w3=pos_level3
	//w4=perc grey
	//w5=perc dark
	//w6=perc_boxcar12 (lateral res/stripe or fov width*100)
	//w7=sigma (gaussian width)
	// assume p1<p2<p3
	wave w
	variable x
	wave tmp
	make/n=6/o tmp
	// grey and bright
	tmp[0]=(1-w[5])/2
	tmp[1]=w[4]
	tmp[2]=w[6]
	tmp[3]=w[2]
	tmp[4]=w[3]
	tmp[5]=w[7]
	variable itot=twoLevelSimulation(tmp,x)
	// grey and dark
	tmp[0]=(w[5])/2
	tmp[1]=1-w[4]
	tmp[3]=w[1]
	tmp[4]=w[2]
	itot+=twoLevelSimulation(tmp,x)
	// bright and dark
	tmp[0]=(1-w[4])/2
	tmp[1]=w[5]
	tmp[3]=w[1]
	tmp[4]=w[3]
	itot+=twoLevelSimulation(tmp,x)
	// grey
	itot+=0.5*w[4]/sqrt(2*pi*w[7]^2)*exp(-(x-w[2])^2/(2*w[7]^2))
	killwaves tmp
	return w[0]*itot
End

function Fit_Historgram_3Level(hstStr,W_coef)
	String hstStr
	wave W_coef
	Wavestats/Q $hstStr
	wave w=$hstStr
	make/o W_coef={1.303e+05,82.81,V_maxloc,168.69,0.52604,0.48132,0.21,9.8608} //1706
	//W_coef={1.3159e+05,110.07,109.39,150.75,0.12328,0.18373,0.36,14.534}
	//W_coef={1.3111e+05,88.592,122.94,149.71,0.4282,0.31848,0.22,9.8503} //030621_002
	// W_coef={1.3103e+05,111.89,145.8,179.65,0.41137,0.31409,0.22,10.769}//030620_120
	//W_coef={1.3086e+05,V_maxloc-60,V_maxloc-40,V_maxloc,0.23,0.07084,0.22,11.689}//030620_98
	//W_coef={1.3134e+05,130.49,164.09,188.15,0.38519,0.31776,0.22,11.689}//030620_102
	// W_coef={1.3177e+05,113.41,125.11,155.16,0.4,0.21372,0.2,15.692}// 105?
	//  W_coef={1.3137e+05,85.555,128.55,151.55,0.44154,0.35412,0.44137,10.66}  // 030620_101 & 103
	// W_coef={1.0466e+05,94.092,124.75,177.87,0.51345,0.81438,0.3,13.631} //20020420
	FuncFit/H="00000010" threeLevelSimulation W_coef $hstStr /D   
	
	//W_coef={1.0466e+05,94.092,124.75,177.87,0.51345,0.81438,13.631} //20020420
	//FuncFit/H="0000000" threeGausshistogram W_coef $hstStr /D   
	//w0=scaling
	//w1=pos_level1
	//w2=pos_level2
	//w3=pos_level3
	//w4=perc grey
	//w5=perc dark
	//w6=perc_boxcar12 (lateral res/stripe width*100)
	//w7=sigma (gaussian width)
	Make/n=20/o W_fit
	W_fit[0]=V_chisq
	W_fit[4]=(1-W_coef[4])*W_coef[5]
	W_fit[5]=W_coef[4]
	W_fit[6]=(1-W_coef[4])*(1-W_coef[5])
	W_fit[1]=W_coef[1]
	W_fit[2]=W_coef[2]
	W_fit[3]=W_coef[3]	 
	// printout
	print hstStr, ": coef 1 (dark)", W_fit[4]*100," at: ", W_coef[1]
	print hstStr, ": 100- coef 1 (neutral)",W_fit[5]*100," at: ", W_coef[2]
	print hstStr, ": coef 1 (bright)",W_fit[6]*100," at: ", W_coef[3]
	print hstStr, "sum",W_fit[4]+W_fit[5]+W_fit[6]
	print hstStr, ": chisquare",V_chisq
	print "dark","neutral","bright","level dark", "level neutral", "level bright"
	print W_fit[4]*100,	W_fit[5]*100, W_fit[6]*100, W_coef[1],  W_coef[2], W_coef[3]
End

function Fit_Historgram_2Level(hstStr,W_coef)
	String hstStr
	wave W_coef
	Wavestats/Q $hstStr
	wave w=$hstStr
	//W_coef={1.2558e+05,0.5,0.15,70,170,12.413} //20030620_030-065
	W_coef={1.2558e+05,0.5,0.22,100,140,12.413} //20020420
	// make/o W_coef={1.3105e+05,0.3,0.16,V_maxloc-30,V_maxloc,10.616} //20030620_085-095
	//w1=perc_level1 (percentage of histogram level 1)
	//w2=perc_boxcar (lateral res/stripe width*100)
	//w3=pos_level1
	//w4=pos_level2
	//w5=sigma (gaussian width)
	FuncFit/H="001000" twoLevelSimulation W_coef $hstStr /D   
	
	//W_coef={1.2558e+05,0.5,100,140,12.413} 
	//FuncFit/H="00000" twoGaussHistogram W_coef $hstStr /D   
	//w7=sigma (gaussian width)
	Make/n=20/o W_fit
	W_fit[0]=V_chisq
	W_fit[4]=W_coef[1]
	W_fit[5]=(1-W_coef[1]) //W_coef[1]
	W_fit[6]=NaN //(1-W_coef[1])
	W_fit[1]=W_coef[3] //3
	//W_fit[2]=W_coef[4] //128 //W_coef[3]
	W_fit[3]=NaN
		
	W_fit[2]=W_coef[4] //
	// printout
	print hstStr, ": (dark)", NaN," at: ", 0
	print hstStr, ": coef 1(neutral)",W_fit[5]*100," at: ", W_coef[3]
	print hstStr, ": 1 -coef 1 (bright)",W_fit[6]*100," at: ", W_coef[2]  //4
	print hstStr, ": chisquare",V_chisq
	print "dark","neutral","bright","level dark", "level neutral", "level bright"
	print W_fit[4]*100,W_fit[5]*100, "NaN",W_coef[3],W_coef[2],"NaN"
End

Function threeLevelStripeSimulation(w,x)
	//w0=amplitude (scaling coef)
	//w1=perc_ferromag (percentage of FM phase)
	//w2=resol   (lat res of the microscope  (30 to 40 nm) / stripe periodicty * 100)
	//w3=pos_dark (dark level)
	//w4=pos_para (neutal level)
	//w5=sigma (gaussian width)
	//w6=howmuch_bri (how miuch bright %)
	//w7=pos_bri (bright level)
	// init coef
	//Make/n=8/o W_coef={804.67,70,7,70,128,20,44.185,170}  // 06/17
	//W_coef={804.67,50,7,90,128,10,44.185,150} // 06/21
	//W_coef={804.67,80,5,90,128,10,44.185,150} // 06/21
	Wave w
	Variable x
	variable sigma=2
	variable intensity1=0
	variable intensity2=0
	variable level=0
	if ((w[1]>=w[2])&&(100-w[1])>=w[2])
	//lower and middle peak
       	intensity1=(w[1]-w[2])*1/sqrt(2*pi*w[5]^2)*exp(-(x-w[3])^2/(2*w[5]^2))+(100-w[1]-w[2])*1/sqrt(2*pi*w[5]^2)*exp(-(x-w[4])^2/(2*w[5]^2))
       	if ((w[3]<=x)&&(x<=w[4]))
       		intensity1+=2*w[2]/abs(w[3]-w[4])*(gammp(0.5,((w[4]-x)/w[5])^2/2)/2+gammp(0.5,((x-w[3])/w[5])^2/2)/2)
       	elseif  ((x<w[3])&&(x<w[4]))
             intensity1+=2*w[2]/abs(w[3]-w[4])*abs(gammp(0.5,((w[4]-x)/w[5])^2/2)/2-gammp(0.5,((w[3]-x)/w[5])^2/2)/2)
         elseif   ((w[3]<x)&&(w[4]<x))
             intensity1+=2*w[2]/abs(w[3]-w[4])*abs(gammp(0.5,((x-w[3])/w[5])^2/2)/2-gammp(0.5,((x-w[4])/w[5])^2/2)/2)
       	endif
       	//higher and middle peak
        intensity2=(w[1]-w[2])*1/sqrt(2*pi*w[5]^2)*exp(-(x-w[7])^2/(2*w[5]^2))+(100-w[1]-w[2])*1/sqrt(2*pi*w[5]^2)*exp(-(x-w[4])^2/(2*w[5]^2))
       	if ((w[4]<=x)&&(x<=w[7]))
       	intensity2+=2*w[2]/abs(w[7]-w[4])*(gammp(0.5,((w[7]-x)/w[5])^2/2)/2+gammp(0.5,((x-w[4])/w[5])^2/2)/2)
       	elseif  ((w[4]>x)&&(x<w[7]))
             intensity2+=2*w[2]/abs(w[7]-w[4])*abs(gammp(0.5,((w[7]-x)/w[5])^2/2)/2-gammp(0.5,((w[4]-x)/w[5])^2/2)/2)
             elseif   ((w[7]<x)&&(x>w[4]))
             intensity2+=2*w[2]/abs(w[7]-w[4])*abs(gammp(0.5,((x-w[4])/w[5])^2/2)/2-gammp(0.5,((x-w[7])/w[5])^2/2)/2)
       endif
       elseif ((w[1]<w[2])&&(100-w[1])>w[2])
             level=(w[1]*w[3]+w[4]*(w[2]-w[1]))/w[2]
       	intensity1=(w[2]-w[1])*1/sqrt(2*pi*w[5]^2)*exp(-(x-level)^2/(2*w[5]^2))+(100-(w[2]-w[1])-2*w[1])*1/sqrt(2*pi*w[5]^2)*exp(-(x-w[4])^2/(2*w[5]^2))
       	if ((level<=x)&&(x<=w[4]))
       	intensity1+=(2*w[1])/abs(w[4]-level)*(gammp(0.5,((w[4]-x)/w[5])^2/2)/2+gammp(0.5,((x-level)/w[5])^2/2)/2)
       	elseif  ((level>x)&&(x<w[4]))
             intensity1+=(2*w[1])/abs(w[4]-level)*abs(gammp(0.5,((w[4]-x)/w[5])^2/2)/2-gammp(0.5,((level-x)/w[5])^2/2)/2)
             elseif   ((level<x)&&(x>w[4]))
             intensity1+=(2*w[1])/abs(w[4]-level)*abs(gammp(0.5,((x-level)/w[5])^2/2)/2-gammp(0.5,((x-w[4])/w[5])^2/2)/2)
       	endif
       	//higher and middle peak
       	level=(w[1]*w[7]+w[4]*(w[2]-w[1]))/w[2]
       	intensity2=(w[2]-w[1])*1/sqrt(2*pi*w[5]^2)*exp(-(level-x)^2/(2*w[5]^2))+(100-(w[2]-w[1])-2*w[1])*1/sqrt(2*pi*w[5]^2)*exp(-(w[4]-x)^2/(2*w[5]^2))
       	if ((level>=x)&&(x>=w[4]))
       	intensity2+=(2*w[1])/abs(w[4]-level)*(gammp(0.5,((w[4]-x)/w[5])^2/2)/2+gammp(0.5,((level-x)/w[5])^2/2)/2)
       	elseif  ((level>x)&&(x<w[4]))
             intensity2+=(2*w[1])/abs(w[4]-level)*abs(gammp(0.5,((level-x)/w[5])^2/2)/2-gammp(0.5,((w[4]-x)/w[5])^2/2)/2)
             elseif   ((level<x)&&(x>w[4]))
             intensity2+=(2*w[1])/abs(w[4]-level)*abs(gammp(0.5,((x-w[4])/w[5])^2/2)/2-gammp(0.5,((level-x)/w[5])^2/2)/2)
       	endif
        elseif (((100-w[1])<w[2])&&(w[1]>w[2]))
       	level=((100-w[1])*w[4]+w[3]*(w[2]-(100-w[1])))/w[2]
       	intensity1=(100-(w[2]-(100-w[1]))-2*(100-w[1]))*1/sqrt(2*pi*w[5]^2)*exp(-(w[3]-x)^2/(2*w[5]^2))+(w[2]-(100-w[1]))*1/sqrt(2*pi*w[5]^2)*exp(-(level-x)^2/(2*w[5]^2))
       	if ((w[3]<=x)&&(x<=level))
       	intensity1+=2*(100-w[1])/abs(level-w[3])*(gammp(0.5,((x-w[3])/w[5])^2/2)/2+gammp(0.5,((level-x)/w[5])^2/2)/2)
       	elseif  ((level>x)&&(x<w[3]))
             intensity1+=2*(100-w[1])/abs(level-w[3])*abs(gammp(0.5,((level-x)/w[5])^2/2)/2-gammp(0.5,((w[3]-x)/w[5])^2/2)/2)
             elseif   ((level<x)&&(x>w[3]))
             intensity1+=2*(100-w[1])/abs(level-w[3])*abs(gammp(0.5,((x-w[3])/w[5])^2/2)/2-gammp(0.5,((x-level)/w[5])^2/2)/2)
       	endif
       	//higher and middle peak
       	level=((100-w[1])*w[4]+w[7]*(w[2]-(100-w[1])))/w[2]
       	intensity2=(100-(w[2]-(100-w[1]))-2*(100-w[1]))*1/sqrt(2*pi*w[5]^2)*exp(-(w[7]-x)^2/(2*w[5]^2))+(w[2]-(100-w[1]))*1/sqrt(2*pi*w[5]^2)*exp(-(level-x)^2/(2*w[5]^2))       	
       	if ((w[7]>=x)&&(x>=level))
       	intensity2+=2*(100-w[1])/abs(level-w[7])*(gammp(0.5,((x-level)/w[5])^2/2)/2+gammp(0.5,((w[7]-x)/w[5])^2/2)/2)
       	elseif  ((level>x)&&(x<w[7]))
             intensity2+=2*(100-w[1])/abs(level-w[7])*abs(gammp(0.5,((w[7]-x)/w[5])^2/2)/2-gammp(0.5,((level-x)/w[5])^2/2)/2)
             elseif   ((level<x)&&(x>w[7]))
             intensity2+=2*(100-w[1])/abs(level-w[7])*abs(gammp(0.5,((x-level)/w[5])^2/2)/2-gammp(0.5,((x-w[7])/w[5])^2/2)/2)
       	endif
       endif
      //weighting 
      return w[0]*(intensity2*w[6]+ intensity1*(100-w[6]))/100    
End

// MnAs functions

Macro prepMnAs()
	// parameters for 3LEVEL
	Make/o W_coef={1.303e+05,82.81,125.33,168.69,0.52604,0.48132,0.3,9.8608} //170603
	//W_coef={77321,88.455,119.83,150.08,0.44831,0.27373,0.3,11.096}//210603
	 W_coef={1.204e+05,99.76,127.76,169.77,0.61861,0.76224,0.3,14.144} //200402
end

Macro MnAs(tmpStr)
	// parameters for 3LEVEL
	String tmpStr
	Prompt tmpStr, "enter wave",popup WaveList("FIT_*",";","")
	variable m=dimsize($tmpStr,1)
	// historgram fit result
	make/n=(m)/o dark; dark[]=$tmpStr[4][p]
	make/n=(m)/o neutral; neutral[]=$tmpStr[5][p]
	make/n=(m)/o bright; bright[]=$tmpStr[6][p]
	make/n=(m)/o darkpos; darkpos[]=$tmpStr[1][p]
	make/n=(m)/o neutralpos; neutralpos[]=$tmpStr[2][p]
	make/n=(m)/o brightpos; brightpos[]=$tmpStr[3][p]
	make/n=(m)/o chisq; chisq[]=$tmpStr[0][p]
	make/n=(m)/o tmp=128-neutralpos
	make/n=(m)/o temperature=NaN
	brightpos+=tmp
	darkpos+=tmp
	neutralpos+=tmp
	String infoStr=tmpStr[4,strlen(tmpStr)-1]+"_TEMP"
	temperature[] =p //$infoStr[$tmpStr[19][p]]
	// particle analysis result
	make/n=(m)/o darkp; darkp[]=$tmpStr[14][p]
	make/n=(m)/o neutralp; neutralp[]=$tmpStr[15][p]
	make/n=(m)/o brightp; brightp[]=$tmpStr[16][p]
	make/n=(m)/o darkposp; darkposp[]=$tmpStr[11][p]
	make/n=(m)/o neutralposp; neutralposp[]=$tmpStr[12][p]
	make/n=(m)/o brightposp; brightposp[]=$tmpStr[13][p]
	tmp=128-neutralposp
	brightposp+=tmp
	darkposp+=tmp
	neutralposp+=tmp
	// data layout and graph 
	Display dark,neutral,bright vs temperature
	append darkp,neutralp,brightp vs temperature
	Label left "\\Z14\\F'Arial'area fraction (%)";DelayUpdate
	Label bottom "\\Z14\\F'Arial'Temperature"
	ModifyGraph mirror=1
	ModifyGraph mode=3,marker=19,rgb(dark)=(0,0,52224);DelayUpdate
	ModifyGraph rgb(neutral)=(39168,39168,39168),rgb(bright)=(65280,43520,0)
	ModifyGraph rgb(darkp)=(0,0,13056),rgb(neutralp)=(32768,32768,65280);DelayUpdate
	ModifyGraph rgb(brightp)=(52224,17408,0)
	Legend/C/N=text0/A=LT
	ModifyGraph width=283,height=283
	SetAxis left 0,1
	// 
	Display darkpos,neutralpos,brightpos vs temperature
	append darkposp,neutralposp,brightposp vs temperature
	ModifyGraph mode=3,marker=19,rgb(darkpos)=(0,0,52224);DelayUpdate
	ModifyGraph rgb(neutralpos)=(39168,39168,39168),rgb(brightpos)=(65280,43520,0)
	ModifyGraph rgb(darkposp)=(0,0,13056),rgb(neutralposp)=(32768,32768,65280);DelayUpdate
	ModifyGraph rgb(brightposp)=(52224,17408,0)
	ModifyGraph mirror=1
	Label left "\\Z14\\F'Arial'area fraction (%)";DelayUpdate
	Label bottom "\\Z14\\F'Arial'Temperature (C)"
	Label left "\\Z14\\F'Arial'histogram level (0-255)"
	Legend/C/N=text0/A=LT
	ModifyGraph width=283,height=283
	SetAxis left 0,255 
	// 
	display chisq vs temperature
	
	killwaves/z tmp
end

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// useful macros
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
			
Macro dataFitLayout1()
	ModifyGraph margin=40; ModifyGraph margin(top)=20,margin(right)=20
	ModifyGraph margin(left)=30
	ModifyGraph width=240,height=240
	label left "\F'Arial'\Z16intensity (a.u.)"
	Label bottom "\\F'Arial'\\Z16binding energy (eV)"
	ModifyGraph grid(bottom)=1,mirror=1,nticks(left)=0,nticks(bottom)=10
end

Macro dataFitLayout2()
	Label bottom "\\F'Arial'\\Z16Time (min)"
	ModifyGraph margin(left)=71,margin(bottom)=42,margin(top)=28,margin(right)=37;DelayUpdate
	ModifyGraph width=255,height=170
	label left "\F'Arial'\Z16peak intensity (a.u.)"
	Label bottom "\\F'Arial'\\Z16time (min)"
	ModifyGraph grid(bottom)=1,mirror(bottom)=1,nticks(left)=5,nticks(bottom)=10
end

Macro LEEDLayout()
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	ModifyGraph nticks=0
	ModifyGraph margin=1
	ModifyImage $tmpStr ctab= {*,*,YellowHot,0}
end

function convert_T_mv(x)
	variable x
	return 25+5.4115+67.558*x-1.9296*x^2+0.096065*x^3-0.0016197*x^4+273.16
end	
	
Macro mv_to_Kelvin(mv)
	variable mv
	print convert_T_mv(mv) 
end	

Macro DrawCristallographicAxes(rotangle,len,labelx,labely)
variable rotangle=-19.6486 // 2003 03 08
variable len=250
string labelx="[1,-1,0]"
string labely="[0,0,1]"
variable posx=pcsr(A)
variable posy=qcsr(A)
variable offset=-10
SetDrawEnv arrow=1,linefgc= (65280,21760,0),xcoord= top,ycoord= left;DelayUpdate
DrawLine posx-len*cos(rotangle*Pi/180),posy+len*sin(rotangle*Pi/180),posx+len*cos(rotangle*Pi/180),posy-len*sin(rotangle*Pi/180)
SetDrawEnv arrow=1,linefgc= (65280,21760,0),xcoord= top,ycoord= left;DelayUpdate
DrawLine posx+len*sin(rotangle*Pi/180),posy+len*cos(rotangle*Pi/180),posx-len*sin(rotangle*Pi/180),posy-len*cos(rotangle*Pi/180)
SetDrawEnv xcoord= top,ycoord= left,textrgb=(65280,21760,0);DelayUpdate
DrawText posx+len*cos(rotangle*Pi/180),posy-len*sin(rotangle*Pi/180)+offset,labelx
SetDrawEnv xcoord= top,ycoord= left,textrgb=(65280,21760,0);DelayUpdate
DrawText posx-len*sin(rotangle*Pi/180),posy-len*cos(rotangle*Pi/180)+offset,labely
End

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// stupid stuff
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

proc PartialPressure(ppo,pph)
	variable ppo,pph
	print ppo*4.4/1.67*1013/760,pph*4.4/1.67*1013/760
end

function frontSize()
	String tmpStr=returnListFirstElement(imageNameList("",""),";")
	tmpStr= removeQuotes(tmpStr)
	variable np=dimSize($tmpStr,0)
	variable nq=dimSize($tmpStr,1)
	waveStats/Q $tmpStr
	variable v_thres=(V_min+V_max)*0.5
	String sizeStr="SZ_"+tmpStr
	make/n=(nq)/o $sizeStr=NaN
	wave wtmp=$tmpStr
	wave wsz=$sizeStr
	variable p1=0
	variable p2=np-1
	variable pmed
	variable q=0
	do
		p1=0
		do
			p1+=1
		while(wtmp[p1][q]<V_thres)
		p2=np-1
		do
			p2-=1
		while(wtmp[p2][q]<V_thres)
		pmed=(p1+p2)*0.5
		// DimOffset(waveName, dim) + ScaledDimPos *DimDelta(waveName,dim)
		wsz[q]=abs(p1-pmed)+abs(p2-pmed)
		//wsz[q]=2*abs(p1-186)  //+abs(p2-pmed)
		q+=1
	while(q<nq)
end

Function MnAs_peridicity()
	variable p1=1/0.0274
	variable p2=1/0.0215 //45
	variable l=sqrt(p1^2+p2^2)
	variable theta = asin(p1/l)
	variable p=p1*cos(Pi/2-theta)
	print p/435*5
	return p
end