#pragma rtGlobals=1		// Use modern global access method.
#pragma IgorVersion = 4.0
#pragma version = 3.0
#include  "KB_mirrors multi 3.0", version>=3.0
#include  "KB_mirrors deposit 3.0", version>=3.0
#include  "KB_mirrors rays 3.0", version>=3.0
#include  "KB_mirrors util 3.0", version>=3.0
#include <Differentiate XY>

// #include "KB_mirrors", version>=3.0
// for packages menu
//		"K-B Mirror analysis",Execute/P "INSERTINCLUDE  \"KB_mirrors 3.0\", version>=3.0;Execute/P "COMPILEPROCEDURES ";Execute/P "KB_mirror_InitPackage()"
//		help = {"Load procedures for analyzing K-B mirrors"}

Menu "K-B mirrors"
	"Fit an LTP Mirror Profile", Call_Fit_a_Stripe("",NaN,NaN,"")
	"Show height from LTP", DisplayHeightProfile("","")
	"Load LTP slope file", Load_APS_LTPslopeFile("")
	Submenu "Multiple"
		"Fit Multiple Stripes", Fit_multiple_stripes("",3)
		"Load Multiple Stripes", Load_multiple_Stripes("","")
		"Save Multiple Stripes for Metrology", Save_multiple_Depositions("")
	End
	"-"
	"Make A Deposition Profile", FindDepositionProfile("",NaN,NaN,NaN,NaN,NaN,NaN)
	"Save A Deposition Profile", CallMakeDepositionForMetrology("")
	"Display Both Deposit Profiles", Display_BothDepositProfiles("","")
	"-"
	"Show the NoteBook", DoWindow /F KB_Mirror_Results
	"Put Plots in NoteBook", AppendPlotsToNotebook("","","")
	"-"
	"Show Spot Size", CallShowSpotSize()
	"Display the actual rays", MakeRaysPlot()
	"Make Perfect Ellipse", MakePerfectEllipse(0,0,0,0)
	"Make Perfect Circle", MakePerfectCircle(0,0,0)
	"-"
	Submenu "Debugging"
		"Debugging ON", root:Packages:KBmirror:debugLevel=10
		"Debugging OFF", root:Packages:KBmirror:debugLevel=0
		"Debugging graph", Graph_debugging()
	End
End



Function Load_APS_LTPslopeFile(fileName)
	String filename								// if null string, you are prompted
	// loads LTP files from APS metrology

	KB_mirror_InitPackage()					// ensure package was inited
	String line, dateTaken, comment, wnote="flip=0"
	Variable refNum, V_flag
	if (strlen(PathList("raw",";",""))<2)		// path 'raw' does not exist, create it
		NewPath raw
	endif

	Open /Z=2/R/P=raw/M="pick a slope data file (.SLP)" refNum as fileName
	if (V_flag)
		Abort "unable to open data file"
	endif
	String fullFileName = S_fileName
	// go down to first line starting with "LTPII"
	Variable i=0
	do
		FReadLine refNum, line
	while (strsearch(line,"LTPII",0)!=0)
	FReadLine refNum, line ; FReadLine refNum, line		// skip 2 more lines
	FReadLine refNum, dateTaken
	FReadLine refNum, comment

	// check to see if this file has had the curvature removed
	FStatus refNum								// get V_filePos, current file position
	i = 0										// check to see if curvature was removed from this file
	do											// look for the measured x-units
		FReadLine refNum, line
		i += 1
	while (i<40 && !stringmatch(line,"I2S * PS1*"))
	if (i<40)									// found the 'PS1'
		DoAlert 1, "this file had curvature removed, continue?"
		if (V_flag!=1)
			Close refNum
			return -1
		endif
	endif
	FSetPos refNum, V_filePos					// return to previous position in file

	// now read in the measured x-units
	String xunits
	i = 0
	do											// look for the measured x-units
		FReadLine refNum, line
		xunits = StringByKey("X unit", line,":",";")
		i += 1
	while (i<40 && strlen(xunits)<1)
	xunits = ZapControlCodes(xunits)
	if (strlen(xunits)<1)
		xunits = "mm"
	endif
	Variable Xscaling=LengthScalng2meters(xunits)
	Close refNum

	dateTaken = ZapControlCodes(dateTaken)
	comment = ZapControlCodes(comment)
	wnote += ";dateTaken="+dateTaken
	if (strlen(comment)>2)
		wnote += ";comment="+comment
	endif

	i=FindFirstDataLine(fullFileName)			// first data line in file
	if (i<0)
		Abort "unable to re-open data file '"+S_fileName+"'"
	endif
	LoadWave /A/B="C=1,F=0,N=Temp_radX;  C=1,F=0,N=Temp_rad;"/J/L={0,i,0,0,2}/O fullFileName

	// trim off trailing NaN's, which come from extra blank lines
	Wave wav = $StringFromList(0, S_waveNames)
	Variable N = FindFirstNaN(wav)
	if (N>0)
		i = 0
		do
			line = StringFromList(i, S_waveNames)
			Redimension/N=(N) $line
			i += 1
		while (i<=ItemsInList("S_waveNames"))
	elseif (N==0)						// nothing read in
		i = 0							// kill the waves
		do
			line = StringFromList(i, S_waveNames)
			KillWaves/Z $line
			i += 1
		while (strlen(line)>0)
		Abort "nothing read in, check that there are no lines after data ends"
	endif

	print "    from directory  '"+S_Path+"'"
	wnote += ";fileName="+S_fileName+";path="+S_path
	wnote += ";"

	Wave Temp_rad=Temp_rad
	Wave Temp_radX=Temp_radX
	Temp_radX *= Xscaling						// put x scaling into m not mm

	Variable x0,x1,xmid, y0
	x0 = Temp_radX[0]							// center the mirror
	x1 = Temp_radX[inf]
	xmid = (x0+x1)/2
	Temp_radX -= xmid
	y0 = faverage(Temp_rad,-inf,inf)			// make slope symmetric about zero
	Temp_rad -= y0
	SetScale d 0,0,"rad", Temp_rad
	SetScale d 0,0,"m", Temp_radX
	String name,nameX

	Variable length=abs(x1-x0)
	Variable lengthmm = round(1000*length)
	wnote += "length="+num2str(length)+";"
	name ="s"+num2istr(lengthmm)+"mmRad"
	nameX = name+"X"

	wnote += "Xdata="+nameX+";Ydata="+name+";"
	Note Temp_rad, wnote
	Note Temp_radX, wnote
	Rename Temp_rad $name
	Rename Temp_radX $nameX
End

Function FindFirstNaN(wav)
	Wave wav
	Variable N=numpnts(wav)
	Variable i=-1
	do
		i += 1
	while(i<N && !numtype(wav[i]))
	return (numtype(wav[i])) ? i : -1
End




Function/T DisplayHeightProfile(Srad,SradX)
	String Srad						// slopes from LTP
	String SradX					// x positions of slopes

	if (exists(Srad)!=1 || exists(SradX)!=1)
		Prompt Srad, "choose wave with slope data", popup, ListOfSlopeWaves()
		DoPrompt "select data from LTP", Srad
		if (V_flag)
			Abort "user aborted"
		endif
		SradX = StringByKey("Xdata", note($Srad),"=")
		if (exists(SradX)!=1)
			Prompt Srad, "choose wave with slope data", popup, WaveList("*",";","")
			Prompt SradX, "choose X wave that goes with slope", popup, WaveList("*X",";","")
			DoPrompt "re-select data from LTP", Srad,SradX
			if (V_flag)
				Abort "user aborted"
			endif
		endif
	endif

	Variable left,top=150						// position of graph to create
	Variable width=425							// graph width

	Wave rad = $Srad
	Wave radX = $SradX
	String waveNote=note(radX)
	Variable length = NumberByKey("length",waveNote,"=")
	Wave height=$HeightFromLTP(rad)			// make the measured height profile using angles from LTP

	// now make the window
	left=0
	if (stripeType(length)==2)						// a standard short mirror
		String info=StringByKey("SCREEN1", IgorInfo(0))
		Variable right = str2num(StringFromList(2, info[strsearch( Info,",RECT=",0)+6,inf] , ","))-10
		left = (10+2*width > right) ? (right - (10+width)) : width+10	// try to fit graphs side by side
	endif
	Display /W=(5+left,top,5+width+left,top+290) height vs radX
	PauseUpdate
	ModifyGraph mode=0, lSize=2, rgb=(0,3,65280)
	ModifyGraph tick=2, minor=1
	ModifyGraph mirror=1, lowTrip=0.1, standoff=0
	ModifyGraph lblPos(left)=40,lblPos(left)=40
	ModifyGraph wbRGB=(57000,62000,65535)
	Label bottom "distance  (\\U)"
	DoWindow /C $UniqueName("GraphHeight",6,0)
	// ellipse, ww={xoffset,angOffset,s2,tth,s1}
	// circle,  ww={angOffset,roc}
	String title = StringByKey("title", waveNote,"=")
	String num
	num = "\\JR"+title+"\\Z10"
	if (strlen(title)>0)
		num += "\r"
	endif
	num += "file = '"+StringByKey("fileName", waveNote,"=")+"'"
	Textbox/N=text2/F=0/S=3/A=MT/Y=4/B=1 num
	AppendText/N=text2 StringByKey("dateTaken", waveNote,"=")
	num = StringByKey("comment", waveNote,"=")
	if (strlen(num)>0)
		AppendText/N=text2 num
	endif
	AppendText/N=text2 StringByKey("path", waveNote,"=")
	SetWindow kwTopWin note="height"
	return WinName(0,1)
End










Function Call_Fit_a_Stripe(data,s1,fit_s2,windowID)
	String data				// name of slope data wave to fit
	Variable s1				// distance from center of Source to center of mirror (m)
	Variable fit_s2
	String windowID

	Variable s2											// distance from ideal center of mirror to focus
	Variable length										// length of mirror (m)
	String dataX=""										// name of data x-wave
	if (exists(data)==1)
		dataX = StringByKey("Xdata", note($data),"=")
	endif
	KB_mirror_InitPackage()							// ensure package was inited

	Variable needPrompt=0
	needPrompt = needPrompt || (s1<1) || numtype(s1)
	needPrompt = needPrompt || (fit_s2<1) || (fit_s2>4) || numtype(fit_s2)
	needPrompt = needPrompt || (exists(dataX)!=1)
	if (needPrompt)
		Prompt data, "choose wave with slope data", popup, ListOfSlopeWaves()
		DoPrompt "mirror slope data to fit", data
		if (V_flag)
			Abort "user aborted"
		endif
		dataX = StringByKey("Xdata", note($data),"=")
		if (exists(dataX)!=1)
			Prompt data, "choose wave with slope data", popup, WaveList("*",";","")
			Prompt dataX, "choose X wave that goes with slope", popup, WaveList("*X",";","")
			DoPrompt "mirror slope data to fit", data,dataX
			if (V_flag)
				Abort "user aborted"
			endif
		endif
		length = NumberByKey("length",note($data),"=")
		if (numtype(length))
			Wave xdataX=$dataX
			length=abs(xdataX[numpnts(xdataX)-1] - xdataX[0])
		endif
		if (numtype(s1) || (s1<1))
			s1 = numtype(s1) ? NumVarOrDefault("root:Packages:KBmirror:last_s1", 60) : s1
			s1 = SelectNumber(stripeType(length)-1,s1,64.5,36.5)	// for standard long and short mirrors
		endif
		if (numtype(s2) || (s2<=0))							// s2 is starting point of fit
			s2 = numtype(s2) ? NumVarOrDefault("root:Packages:KBmirror:last_s2", 0.1) : s2
			s2 = SelectNumber(stripeType(length)-1,s2,0.130,0.060)	// for standard long and short mirrors
		endif
		fit_s2 = (fit_s2>=1 && fit_s2 <=4) ? fit_s2 : NumVarOrDefault("root:Packages:KBmirror:fit_s2", 3)
		Prompt fit_s2, "type of fit", popup, "vary s2;Fix s2;Fix s2 & xoffset;circle"
		Prompt windowID, "identifier for window"

		String s1list=MakePromptString(s1,"64.5 (source on 34);36.5 (L5 on 34);other") ; s1=0
		String s2list=MakePromptString(s2,"0.13 (long mirror);0.060 (short mirror);other") ; s2=0
		Prompt s1, "distance from mirror to source (m)", popup, s1list
		Prompt s2, "distance from mirror to focus (m),  fixed value or start for fitting", popup, s2list
		DoPrompt "fit a mirror stripe", fit_s2,windowID,s1,s2
		if (V_flag)
			Abort "user aborted"
		endif
		s1 = str2num(StringFromList(s1-1,s1list))
		s2 = str2num(StringFromList(s2-1,s2list))
		if (numtype(s1) || numtype(s2))
			Prompt s1, "distance from mirror to source (m)"
			Prompt s2, "distance from mirror to focus (m),  fixed value or start for fitting"
			s1 = str2num(StringFromList(0,s1list))		// reset to original
			s2 = str2num(StringFromList(0,s2list))
			DoPrompt "make a deposition profile", s1,s2		// let user enter anything
			if (V_flag)
				return 1
			endif
		endif
	endif
	if (numtype(s1) || numtype(s2) || s1<1 || s2<.001)
		Abort "cannot use s1="+num2str(s1)+"  or  s2="+num2str(s2)
	endif
	NVAR last_s1 = root:Packages:KBmirror:last_s1
	NVAR last_s2 = root:Packages:KBmirror:last_s2
	NVAR lastfit_s2 = root:Packages:KBmirror:fit_s2
	last_s1 = s1
	last_s2 = s2
	lastfit_s2 = fit_s2

	String title,spotError
	spotError = "spotSize"+num2istr(round(length*1000))+"mm"

	s2 = Fit_a_Stripe(data,dataX,s1,s2,fit_s2,windowID)

	if (s2<0)
		Abort "fit failed"
	endif

	String Scorelate = ""
	if (exists("M_Covar"))
		Wave M_Covar=M_Covar
		Scorelate = "correlate_"+num2istr(round(length*1000))+"mm"
		Duplicate/O M_Covar, $Scorelate					// make the correlation matrix
		Wave corelate=$Scorelate
		corelate = M_Covar[p][q]/sqrt(M_Covar[p][p]*M_Covar[q][q])
		KillWaves M_Covar
	endif

	String sww = "ww"+data
	String Sres = "Res_"+data
	String fitWave = "fit_"+data

	// first make the spot error
	Duplicate/O $Sres $spotError
	Wave res = $Sres
	Wave spot = $spotError
	Wave radX = $dataX
	spot = spotErrorCircle(radX[p],res[p],s2)
	SetScale d 0,0,"m", $spotError
	// done making spot error  

	DisplayOneFit(fitWave,Sres,spotError,Scorelate,windowID)		//  Display results of the fit
	SetWindow kwTopWin note="fitLTP"
End

Function Fit_a_Stripe(data,dataX,s1,s2,fit_s2,windowID)
	String data, dataX
	Variable s1				// distance from center of Source to center of mirror (m)
	Variable s2				// distance from ideal center of mirror to focus (m)
	Variable fit_s2			// 1=vary s2,   2=Fix s2,   3=Fix s2 & xoffset,   4=circle
	String windowID

	if (s1<1)
		Abort "s1 is too small, s1="+num2str(s1)+" m"
	endif
	String sww,Sres,title,fitWave
	Variable roc										// radius of curvature (used with circle only)
	Variable iflip
	Variable xoffset0									// x offet of the mirror
	Variable angOffset0 = 1e-4
	Variable tth0 = .006
	Variable width=425									// graph width
	String windowName
	Wave rad = $data
	Wave radX = $dataX
	Variable length=abs(radX[0]-radX[inf])			// length of mirror (m)
	sprintf windowName, "Graph_%s_%dmm",windowID,round(length*1000)
	sprintf title, "%d mm Mirror",round(length*1000)
	windowName = UniqueName(windowName,6,0)
	xoffset0 = (fit_s2==3) ? 0 : -5e-3 + 18e-3		// xoffset = 0 when not fitting it
	switch(stripeType(length))	
		case 1:											// a standard long mirror
			NVAR range0 = root:Packages:KBmirror:rangeL0	// fit range (in points) for long mirror
			NVAR range1 = root:Packages:KBmirror:rangeL1
			angOffset0 = 4e-5
			break
		case 2:											// a standard short mirror
			NVAR range0 = root:Packages:KBmirror:rangeS0	// fit range (in points) for short mirror
			NVAR range1 = root:Packages:KBmirror:rangeS1
			String info=StringByKey("SCREEN1", IgorInfo(0))
			Variable right = str2num(StringFromList(2, info[strsearch( Info,",RECT=",0)+6,inf] , ","))-10
			angOffset0 = 1e-4
			break
	endswitch

	sww = "ww"+data
	Sres = "Res_"+data
	fitWave = "fit_"+data
	if (LengthScalng2meters(WaveUnits(radX,-1))!=1)
		Abort "position along mirrors must be in meters, not mm"
	endif
	Variable n = numpnts(rad)
	range0 = (numtype(range0)) ? n/7 : range0
	range1 = (numtype(range1)) ? 6*n/7 : range1

	printf "For a %s",title
	if (fit_s2==4)
		printf ",    fit to a circle\r"
	else
		printf ",    the distance to the source %g m\r", s1
	endif
	printf "  data in '%s',    x-values in '%s', and     parameters in '%s\r",data,dataX,sww

	Variable r0=range0
	Variable r1=range1
	PutUpGraph_GetRange(r0,r1,rad,radX,windowName)
	range0 = r0
	range1 = r1

	String HoldStr
	if (fit_s2==1)
		HoldStr="00001"	// fix s1, vary {xofset, angle offset, s2, tth}
	endif
	if (fit_s2==2)
		HoldStr="00101"	// fix s1, s2, vary {xofset, angle offset, tth}
	endif
	if (fit_s2==3)
		HoldStr="10101"	// fix s1, s2, xoffset, vary {angle offset, tth}
	endif

	KillWaves/Z $Sres
	KillWaves/Z M_Covar
	printf "  fitting points in range [%d,%d]\r",range0,range1

	if ((fit_s2==1) || (fit_s2==2) || (fit_s2==3))
		Variable chi0, chi1
		String epsWave = "eps"+data
		Make/O/D/N=5 $sww, $epsWave
		Wave ww=$sww
		Wave eps = $epsWave
		ww={xoffset0,angOffset0,s2,tth0,s1}
		eps = {0.001, 0.0001,  0.01,  0.0001,  1}
		FuncFit /Q/H=HoldStr xyFit ww rad[range0,range1] /D /R /X=radX /E=$epsWave
		chi0 = V_chisq

		FlipWave(rad)				// flip waves and try again to find best fit
		iflip = FlipWave(radX)
		ww={xoffset0,angOffset0,s2,tth0,s1}
		FuncFit /Q/H=HoldStr xyFit ww rad[range0,range1] /D /R /X=radX /E=$epsWave
		chi1 = V_chisq
		if (chi0<=chi1)				// un-flip to original if that was better
			FlipWave(rad)
			iflip = FlipWave(radX)
		endif
		ww={xoffset0,angOffset0,s2,tth0,s1}				// and for the final fit
		FuncFit /H=HoldStr/M=2 xyFit ww rad[range0,range1] /D /R /X=radX /E=$epsWave
		chi0 = V_chisq
		s1 = ww[4]
		s2 = ww[2]
	else
		Make/O/D $sww={angOffset0,s1}						// to fit a circle  {angle offset, radius}
		Wave ww=$sww
		FuncFit circleFit ww rad[range0,range1] /D /R /X=radX
		roc = ww[1]
	endif
	DoWindow /K $windowName							// Kill the window, all down with it

	if (cmpstr(GetDataFolder(0) ,"root"))
		title = "\\f01"+GetDataFolder(0)+"\\f00,  "+title
	endif
	if (iflip)
		title += ",   \\f02Flipped\\f00 "
	endif
	String newNote=note(rad)
	newNote = ReplaceNumberByKey("fitType", newNote, fit_s2, "=" )
	newNote = ReplaceStringByKey("Ydata", newNote, data,"=")
	newNote = ReplaceStringByKey("Xdata", newNote, dataX,"=")
	newNote = ReplaceStringByKey("title", newNote, title,"=")
	newNote = ReplaceNumberByKey("length", newNote, length, "=" )
	Wave W_sigma=W_sigma
	if ((fit_s2==1) || (fit_s2==2) || (fit_s2==3))
		newNote = ReplaceNumberByKey("s1", newNote, s1, "=" )
		newNote = ReplaceNumberByKey("s2", newNote, s2, "=" )
		if (fit_s2==1)
			newNote = ReplaceNumberByKey("s2_sig", newNote, W_sigma[2], "=" )
		endif
		newNote = ReplaceNumberByKey("xoffset", newNote, ww[0], "=" )
		if (W_sigma[0]!=0)
			newNote = ReplaceNumberByKey("xoffset_sig", newNote, W_sigma[0], "=" )
		endif
		newNote = ReplaceNumberByKey("angOffset", newNote, ww[1], "=" )
		newNote = ReplaceNumberByKey("angOffset_sig", newNote, W_sigma[1], "=" )
		newNote = ReplaceNumberByKey("tth", newNote, ww[3], "=" )
		newNote = ReplaceNumberByKey("tth_sig", newNote, W_sigma[3], "=" )
	endif
	if (fit_s2==4)
		newNote = ReplaceNumberByKey("angOffset", newNote, ww[0], "=" )
		newNote = ReplaceNumberByKey("angOffset_sig", newNote, W_sigma[0], "=" )
		newNote = ReplaceNumberByKey("roc", newNote, ww[1], "=" )
		newNote = ReplaceNumberByKey("roc_sig", newNote, W_sigma[1], "=" )
		newNote = ReplaceNumberByKey("s1", newNote, s1, "=" )
		newNote = ReplaceNumberByKey("s2", newNote, s2, "=" )
	endif
	newNote = ReplaceNumberByKey("range0", newNote, range0, "=" )
	newNote = ReplaceNumberByKey("range1", newNote, range1, "=" )

	Note/K $fitWave
	Note $fitWave newNote
	Note/K $Sres
	Note $Sres newNote
	if (exists("M_Covar"))
		Note/K M_Covar
		Note M_Covar newNote
	endif

	CopyScales /I $data, $Sres
	return s2
End

Function DisplayOneFit(fitWave,Sres,spotError,corelate,windowID)
	String fitWave							// result fitted wave
	String Sres								// residual error from the fit
	String spotError						// holds name of wave with spot errors
	String corelate
	String windowID						// name of graph window to create

	Variable left,top						// position of graph to create
	Variable width=425						// graph width
	String windowName

	String waveNote=note($fitWave)
	String data=StringByKey("Ydata", waveNote,"=")	// slope data used in fit
	String dataX=StringByKey("Xdata", waveNote,"=")	// x-values of original slope data
	Variable fitType=NumberByKey("fitType", waveNote,"=")
	Variable length=NumberByKey("length", waveNote,"=")
	Variable s1=NumberByKey("s1", waveNote,"=")	// distance from Source to center of mirror (m)
	Variable s2	=NumberByKey("s2", waveNote,"=")	// distance from ideal center of mirror to focus (m)
	if (numtype(fitType) || numtype(s1) || numtype(s2))
		Abort "wave note from the fitted wave does not have proper info"
	endif
	if (exists(data)!=1 || exists(dataX)!=1)
		Abort "Slope data or X values cannot be found from wave note"
	endif
	if (strlen(windowID)>0)
		windowID = "_"+windowID
	endif

	sprintf windowName, "Graph%s_%dmm",windowID,round(length*1000)
	windowName = UniqueName(windowName,6,0)

	// first make the spot error
	Wave spot = $spotError
	Variable rms_error = goodRMS($Sres)				// rms error (rad)
	Variable spotFWHM = FWHM_of_errorWave(spot)	// spot FWHM (m)
	WaveStats/Q $spotError
	spotError = spotError+"2"
	Make/O  $spotError={V_min, V_max}
	SetScale/I x Nan,Nan,"", $spotError
	SetScale d 0,0,"m", $spotError
	// done making spot error  

	// now make the window
	left=0
	switch(stripeType(length))	
		case 1:											// a standard long mirror
			top=42										// positions graph vertically
			break
		case 2:											// a standard short mirror
			top=42										// positions graph vertically
			String info=StringByKey("SCREEN1", IgorInfo(0))
			Variable right = str2num(StringFromList(2, info[strsearch( Info,",RECT=",0)+6,inf] , ","))-10
			left = (10+2*width > right) ? (right - (10+width)) : width+10	// try to fit graphs side by side
			break
		default:
		top=150										// positions graph vertically
	endswitch
	Display /W=(5+left,top,5+width+left,top+290) $Sres vs $dataX
	PauseUpdate
	AppendToGraph/L=main_left $data vs $dataX
	AppendToGraph/L=main_left $fitWave
	AppendToGraph/R $spotError
	ModifyGraph mode[1]=3, marker[1]=8, marker[1]=8
	ModifyGraph lSize=2, opaque=1, rgb[1]=(0,43520,65280)
	ModifyGraph tick=2, minor=1
	ModifyGraph mirror(bottom)=1,mirror(main_left)=1
	ModifyGraph lowTrip=0.1, standoff=0
	ModifyGraph lblMargin(right)=12					//	ModifyGraph gfSize=14
	ModifyGraph lblPos(left)=40,lblPos(main_left)=40
	ModifyGraph freePos(main_left)=0
	ModifyGraph axisEnab(left)={0.8,1}
	ModifyGraph axisEnab(main_left)={0,0.75}
	ModifyGraph axisEnab(right)={0.8,1}

	if (stringmatch(data,"corrected*_deriv"))			// a check plot
		ModifyGraph wbRGB=(65535,62000,57000)
	endif
	Label bottom "distance  (\\U)"
	Cursor/P A $data (NumberByKey("range0",waveNote,"="))
	Cursor/P B $data (NumberByKey("range1",waveNote,"="))
	DoWindow /C $windowName

	printf "total rms error = %.3f  µrad,   spot size = %.3f µm\r" rms_error*1e6, spotFWHM*1.e6
	String num, num0, num1, num2="", num3
	Wave W_sigma=W_sigma

	// ellipse, ww={xoffset,angOffset,s2,tth,s1}
	// circle,  ww={angOffset,roc}
	Variable roc,roc_sig
	Variable angOffset, angOffset_sig
	Variable s1_sig, s2_sig
	Variable tth, tth_sig, xoffset, xoffset_sig

	roc = NumberByKey("roc",waveNote,"=")				// and for the circle
	roc_sig = NumberByKey("roc_sig",waveNote,"=")

	angOffset = NumberByKey("angOffset",waveNote,"=")	// for both circle and ellipse
	angOffset_sig = NumberByKey("angOffset_sig",waveNote,"=")

	s1 = NumberByKey("s1",waveNote,"=")					// for ellipse
	s1_sig = NumberByKey("s1_sig",waveNote,"=")
	s2 = NumberByKey("s2",waveNote,"=")
	s2_sig = NumberByKey("s2_sig",waveNote,"=")
	tth = NumberByKey("tth",waveNote,"=")
	tth_sig = NumberByKey("tth_sig",waveNote,"=")
	xoffset = NumberByKey("xoffset",waveNote,"=")
	xoffset_sig = NumberByKey("xoffset_sig",waveNote,"=")

	if (fitType==1)
		sprintf num "center-focus = %.2f ± %.2f mm" (s2-xoffset)*1000, sqrt(s2_sig^2-xoffset_sig^2)*1000
		Textbox/N=text0/F=0/S=3/A=LC/B=1 num
		sprintf num, "s2 = %.2f ± %.2f mm\rxoff = %.2f ± %.2f mm" s2*1000,s2_sig*1000,xoffset*1000,xoffset_sig*1000
		AppendText/N=text0 num
	endif
	if (fitType==2)
		sprintf num "center-focus = %.2f ± %.2f mm" (s2-xoffset)*1000, xoffset_sig*1000
		Textbox/N=text0/F=0/S=3/A=LC/B=1 num
		sprintf num, "s2 = %.1f mm,  xoff = %.2f ± %.2f mm" s2*1000,xoffset*1000,xoffset_sig*1000
		AppendText/N=text0 num
	endif
	if (fitType==3)
		sprintf num, "(center-focus) = s2 = %.1f mm" s2*1000
		Textbox/N=text0/F=0/S=3/A=LC/B=1 num
	endif
	if (fitType==4)					// for fitting a circle
		sprintf num "radius of curvature = %.2f ± %.2f m" ,roc, roc_sig
		Textbox/N=text0/F=0/S=3/A=LC/B=1 num
		sprintf num, "tilt = %.1f ± %.1f µrad\\Z]0" angOffset*1e6,angOffset_sig*1e6
		AppendText/N=text0 num
		sprintf num, "mirror angle = %.3f mrad\\Z]0" mirrorAngleCircle(s1,s2,roc)*1e3
		AppendText/N=text0 num
	else
		sprintf num, "mirror angle = %.2f ± %.2f mrad" tth*500, tth_sig*500
		AppendText/N=text0 num
		sprintf num, "tilt = %.2f ± %.2f mrad,  s1=%gm\\Z]0" angOffset*1e3,angOffset_sig*1000,s1
		AppendText/N=text0 num
	endif
	sprintf num, "slope %.3f µrad rms\rspot %.3f µm FWHM" rms_error*1e6, spotFWHM*1e6
	Textbox/N=text1/F=0/S=3/A=LT/B=1 num

	String title = StringByKey("title", waveNote,"=")
	num = "\\JR"+title+"\\Z10"
	if (strlen(title)>0)
		num += "\r"
	endif
	num += "file = '"+StringByKey("fileName", waveNote,"=")+"'"
	Textbox/N=text2/F=0/S=3/A=RB/Y=3/B=1 num
	AppendText/N=text2 StringByKey("dateTaken", note($data),"=")
	num = StringByKey("comment", note($data),"=")
	if (strlen(num)>0)
		AppendText/N=text2 num
	endif
	AppendText/N=text2 StringByKey("path", note($data),"=")

	if (exists(corelate)==1)
		TextBox/N=text3/F=0/S=3/B=1/T={54,99,150,191}/A=RC/X=4/Y=-16 "	cross-correlations"
		Wave cor=$corelate
		if (fitType==1)
			AppendText/N=text3 "\\Z10\tx offset\tang offset\t   s2\ttheta"
			sprintf num0, "x offset\t%7.4f\t%7.4f\t%7.4f\t%7.4f\r", 		cor[0][0], cor[0][1], cor[0][2], cor[0][3]
			sprintf num1, "ang offset\t%7.4f\t%7.4f\t%7.4f\t%7.4f\r", 	cor[1][0], cor[1][1], cor[1][2], cor[1][3]
			sprintf num2, "s2\t%7.4f\t%7.4f\t%7.4f\t%7.4f\r", 			cor[2][0], cor[2][1], cor[2][2], cor[2][3]
			sprintf num3, "theta\t%7.4f\t%7.4f\t%7.4f\t%7.4f", 			cor[3][0], cor[3][1], cor[3][2], cor[3][3]
		elseif (fitType==3)
			num0 = ""
			AppendText/N=text3 "\\Z10\tang offset\ttheta"
			sprintf num1, "ang offset\t%7.4f\t%7.4f\r", 	cor[1][1], cor[1][3]
			sprintf num3, "theta\t%7.4f\t%7.4f", 			cor[3][1], cor[3][3]
		else
			AppendText/N=text3 "\\Z10\t\tang offset\ttheta"
			sprintf num0, "x offset\t%7.4f\t%7.4f\t%7.4f\r", 		cor[0][0], cor[0][1], cor[0][3]
			sprintf num1, "ang offset\t%7.4f\t%7.4f\t%7.4f\r", 	cor[1][0], cor[1][1], cor[1][3]
			sprintf num3, "theta\t%7.4f\t%7.4f\t%7.4f", 			cor[3][0], cor[3][1], cor[3][3]
		endif
		AppendText/N=text3 num0+num1+num2+num3
	endif
End



Function stripeType(length)			// identify standard type mirrors for setting defaults
	Variable length										// length of mirror (m)
	if (abs(length-.087)<.005)						// a standard long mirror
		return 1
	elseif (abs(length-.037)<.004)					// a standard short mirror
		return 2
	endif
	return 0
End




Function circleFit(ww,xx) : FitFunc
	Wave ww
		//	ww[0] = angle offset (radians)
		//	ww[1] = radius, (radius/2 = focal length) (m)
	Variable xx								// x position  (m)
	Variable r2 = ww[1]^2					// radius^2 of mirror (m)
	return ww[0] + xx / sqrt(r2 - xx^2)	// slope for the circle + offset  (radians)
End

Function xyFit(ww,xx) : FitFunc
	Wave ww
		//	ww[0] = xoffset (m)  (increasing offset moves mirror toward focus)
		//	ww[1] = angle offset (radians)
		//	ww[2] = s2 (distance from center of mirror to focus) (m)
		//	ww[3] = tth (twice grazing angle of center of mirror) (radians)
		//	ww[4] = s1 (distance from center of mirror to source) (m)
	Variable xx												// input position (along mirror, entered in m)
	Variable w = ww[0]+xx									// full position along mirror (m)
	return ellipseAngle(ww[4],ww[2],ww[3],w)+ ww[1]	// the angle of an ellipse (angle in radians)
End

Function ellipseAngle(s1,s2,tth,w)
	// returns angle to ellipse (rad)
	Variable s1				// distance from center of mirror to source (usually 36.5) (m)
	Variable s2				// distance from center of mirror to focus (0.13 for long mirror) (m)
	Variable tth				// twice the grazing angle on mirror (usually 0.006) (rad)
	Variable w				// position along mirror (not x-axis) (m)
							// center of mirror is determined knowing tth
	Variable F2 = (s1^2+s2^2-2*s1*s2*cos(PI-tth))/4	// F is dist. from cent of ellipse to foci (law of cosines)
	Variable F = sqrt(F2)
	Variable a2= (s1+s2)^2 / 4			// a is semi-major axis (s1+s2)=2a
	Variable b2 = a2 - F2					// b is semi-minor axis a^2= b^2+F^2
	Variable x0								// x-coordinate of mirror center (along x-axis, not mirror)
	Variable alpha							// angle opposite side s1
	x0 = F - s2*cos(alpha)					// F = x0 + s2 cos(alpha)
	alpha = asin( s1/(2*F) * sin(PI-tth) )// law of sines
	Variable mu = alpha - tth/2				// mu is angle between x-axis and middle of mirror

	Variable cosine = cos(mu)
	Variable xp = x0 + w*cosine			// x-coordinate (along x-axis, not mirror)
	Variable b = sqrt(b2)

	Variable xs2 = xp^2
	if (xs2>a2)								// requested x value past end of ellipse
		return NaN
	endif
	Variable yp								// dy/ dw
	yp = cosine*(  (b*cosine*xp) / (a2*sqrt(1 - xs2/a2)) - sin(mu)  )
//	yp = b*xp/a2/sqrt(1-x2/a2) - tan(mu)			// old, same to order of cos(mu)^2
	return atan(yp)						// return angle not slope
End



Function MirrorIdealHeight(s1,s2,tth,w)	// returns height above the tangent to ellipse (m)
// ellipse translated so center of mirror at (0,0) and rotated so center spot is horizontal
// this routine uses rotation matricies, not a linear approximation
	Variable s1				// distance from center of mirror to source (usually 36.5) (m)
	Variable s2				// distance from center of mirror to focus (0.13 for long mirror) (m)
	Variable tth				// twice the grazing angle on mirror (usually 0.006) (rad)
	Variable w				// position along mirror (not x-axis) (m)

	// center of mirror is determined knowing tth, independent of anything else
	Variable F2 = (s1^2+s2^2-2*s1*s2*cos(PI-tth))/4	// F is dist. from cent of ellipse to foci (law of cosines)
	Variable F = sqrt(F2)
	Variable a2= (s1+s2)^2 / 4			// a is semi-major axis (s1+s2)=2a
	Variable b2 = a2 - F2					// b is semi-minor axis a^2= b^2+F^2

	Variable x0,y0							// coordinates of mirror center
	Variable alpha							// angle opposite side s1
	alpha = asin( s1/(2*F) * sin(PI-tth) )// law of sines
	Variable mu = alpha - tth/2				// mu is angle between x-axis and middle of mirror, will need later
	x0 = F - s2*cos(alpha)					// F = x0 + s2 cos(alpha)
	y0 = -sqrt((1-x0^2/a2) * b2)		// y0 is y-coord along y-axis of mirror center, take lower root

	Variable x1 = w*cos(mu)				// change from dist along w to along x-axis
	Variable x2 = (x1+x0)^2
	if (x2>a2)								// requested x value past end of ellipse
		return NaN
	endif
	Variable y1 = -sqrt((1-x2/a2) * b2) - y0	// y-coord along y-axis of mirror center, take lower root
	Variable y = y1*cos(-mu) + x1*sin(-mu)	// rotate the y1 into final frame
	return y
End

Function MakePerfectEllipse(s1,s2,length,tth)
	Variable s1					// distance from center of mirror to source (m)
	Variable s2					// distance from center of mirror to focus (m)
	Variable length				// length of mirror (m)
	Variable tth					// twice grazing angle (rad)

	Variable s10=s1, s20=s2, l0=length, tth0=tth
	s1 = ((s1<.1) || numtype(s1)) ? 36.5 : s1			// default to 36.5 m = s1
	s2 = ((s2<.001) || numtype(s2)) ? 0.13 : s2		// default to 0.13 m = s2
	if (abs(s2-.13)<.001)
		length = ((length<1) || numtype(length)) ? .1 : length
	elseif (abs(s2-.06)<.001)
		length = ((length<1) || numtype(length)) ? .04 : length
	endif
	tth = ((tth<1e-5) || numtype(tth)) ? .006 : tth

	if (s10!=s1 || s20!=s2 || l0!=length || tth0!=tth)
		Prompt s1, "distance from source to mirror center, s1 (m)"
		Prompt s2, "distance from mirror center to focus, s2 (m)"
		Prompt length, "mirror length (m)"
		Prompt tth, "twice glancing angle on mirror (rad)"
		DoPrompt "ellipse parameters", s1,s2, length, tth
		if (V_flag)
			return 1
		endif
	endif

	String height, slope
	Variable mm = round(length*1000)
	height = "ellipseHeight"+num2istr(mm)+"mm"
	slope = "ellipseSlope"+num2istr(mm)+"mm"
	if (s10!=s1 || s20!=s2 || l0!=length || tth0!=tth)
		printf "„MakePerfectEllipse(s1=%g, s2=%g, length=%g, tth=%g, height='%s', slope='%s')\r",s1,s2,length,tth,height,slope
	endif

	PauseUpdate
	Make/N=1001/O/D $height, $slope
	SetScale d 0,0,"m", $height
	SetScale d 0,0,"slope", $slope
	SetScale/I x (-length/2),(length/2),"m", $height, $slope
	Wave hh = $height
	Wave sl = $slope
	hh = MirrorIdealHeight(s1,s2,tth,x)
	sl = tan(ellipseAngle(s1,s2,tth,x))		// slope not angle
	String wnote=""
	sprintf wnote, "s1=%g;s2=%g;tth=%g;fl=%g;length=%g", s1,s2,tth,1/(1/s1+1/s2),length

	Note/K $height
	Note/K $slope
	Note $height wnote
	Note $slope wnote
	return 0
End


Function MakePerfectCircle(fl,length,tth)
	Variable fl					// focal length mirror at tth (m)
	Variable length				// length of mirror (m)
	Variable tth					// twice grazing angle (rad)

	Variable fl0=fl, l0=length, tth0=tth
	fl = ((fl<.1) || numtype(fl)) ? 1.0 : fl		// default to 1m focal length
	if (abs(fl-.13)<.001)
		length = ((length<.001) || numtype(length)) ? .1 : length
	elseif (abs(fl-.06)<.001)
		length = ((length<.001) || numtype(length)) ? .04 : length
	endif
	tth = ((tth<1e-5) || numtype(tth)) ? .006 : tth
	length = ((length<.001) || numtype(length)) ? .2: length

	if (fl0!=fl || l0!=length || tth0!=tth)
		Prompt fl, "focal length (m)"
		Prompt length, "mirror length (m)"
		Prompt tth, "twice glancing angle on mirror (rad)"
		DoPrompt "circle parameters", fl, length, tth
		if (V_flag)
			return 1
		endif
	endif

	String height, slope
	Variable mm = round(length*1000)
	height = "circleHeight"+num2istr(mm)+"mm"
	slope = "circleSlope"+num2istr(mm)+"mm"
	if (fl0!=fl || l0!=length || tth0!=tth)
		printf "„MakePerfectCircle(fl=%g, length=%g, tth=%g, height='%s', slope='%s')\r",fl,length,tth,height,slope
	endif

	PauseUpdate
	Make/N=1001/O/D $height, $slope
	SetScale d 0,0,"m", $height
	SetScale d 0,0,"slope", $slope
	SetScale/I x (-length/2),(length/2),"m", $height, $slope
	Wave hh = $height
	Wave sl = $slope

	Variable radius=2*fl/sin(tth/2)	// radius of curvature (m)
	Variable r2 = radius^2
	hh = radius - sqrt(r2 - x*x)
	sl = x / sqrt(r2 - x*x)

	String wnote=""
	sprintf wnote, "radius=%g;fl=%g;", radius,fl
	Note/K $height
	Note/K $slope
	Note $height wnote
	Note $slope wnote
	return 0
End



Function/T HeightFromLTP(rad)
	// give a wave "rad" containing angles (in radians) from LTP, make a new wave with heights,
	// and returns name of heights as results.  Both rad and height use radX for x-values
	Wave rad						// measured LTP angles
	String radName = NameOfWave(rad)
	String radXname = radName+"X"
	Wave radX = $radXname			// x position where angles measured

	String xunits = WaveUnits(radX,-1)
	String yunits = WaveUnits(rad,-1)
	Variable Xscaling = LengthScalng2meters(xunits)
	Variable Yscaling = AngleScalng2radians(yunits)
	Variable N = numpnts(rad)
	String command
	sprintf command, "Interpolate/T=2/N=%d/E=1/Y=rad_CS %s /X=%s", 4*N,radName, radXname
	Execute command
	Wave rad_CS=rad_CS
	rad_CS *= Yscaling
	Variable dx, x0					// used for x-axis rescaling
	dx = DimDelta(rad_CS, 0)*Xscaling
	x0 = DimOffset(rad_CS, 0)*Xscaling

	SetScale/P x (x0),(dx),"m", rad_CS
	SetScale d 0,0,"rad", rad_CS
	// rad_CS is now the measured angle.  Its x-scale is in meters, and y-scale in radians

	rad_CS = tan(rad_CS[p])		// convert angle to slope
	Integrate/T rad_CS				// integrate rad_CS to get heights
	SetScale d 0,0,"m", rad_CS
	String heightName = radName[0,strsearch(radName,"Rad",0)-1]
	if (stringmatch(heightName, "s*" ))		// it is of standard type, e.g. s83mmRad
		heightName = "height"+heightName[1,inf]		// remove leading 's'
	else
		heightName = "height_"+heightName
	endif
//	print "making height wave =",heightName	// name of new wave that will contain heights
	Make/N=(N)/O $heightName
	Wave height = $heightName
	SetScale d 0,0,"m", height
	height = rad_CS(radX[p]*Xscaling)
	WaveStats/Q height
	height -= V_min
	KillWaves/Z rad_CS
	return heightName
End

Function/T HeightFromScaledSlope(rad)
	// with a wave conaining angles (in radians) make a new wave with heights
	Wave rad						// measured LTP angles
	String radName = NameOfWave(rad)

	String xunits = WaveUnits(rad,0)
	String yunits = WaveUnits(rad,-1)
	Variable Xscaling = LengthScalng2meters(xunits)
	Variable Yscaling = AngleScalng2radians(yunits)
	Variable N = numpnts(rad)
	String heightName = radName[0,strsearch(radName,"Rad",0)-1]
	if (stringmatch(heightName, "s*" ))		// it is of standard type, e.g. s83mmRad
		heightName = "height"+heightName[1,inf]		// remove leading 's'
	else
		heightName = "height_"+heightName
	endif
//	print "making height wave =",heightName	// name of new wave that will contain heights

	Duplicate/O $radName $heightName
	Wave height = $heightName
	height *= Yscaling
	Variable dx, x0					// used for x-axis rescaling
	dx = DimDelta($heightName, 0)*Xscaling
	x0 = DimOffset($heightName, 0)*Xscaling
	SetScale/P x (x0),(dx),"m", $heightName
	SetScale d 0,0,"rad", $heightName
	// height is now the measured angle scaled for meters in x and radians in y value

	height = tan(height[p])			// convert angle to slope
	Integrate/T height				// integrate height to get heights
	SetScale d 0,0,"m", $heightName

	WaveStats/Q height
	height -= V_min
	return heightName
End



Static Function FindFirstDataLine(fileNameStr)	// find first line of data in data file
	String fileNameStr							// returns negative for an error

	Variable refNum, V_flag
	Open /R/Z=1 refNum  as fileNameStr
	if (V_flag)
		return -1
	endif

	String line, dataType=""
	Variable bad=-1
	Variable  i=0												// line counter
	do
		i += 1
		FReadLine refNum, line
		if (stringmatch(line,"*Data type:*"))					// found the data type line
			sscanf line, "Data type:%s", dataType
			bad = (strsearch(dataType,"SLOPE",0)<0) ? -2 : 0	// = 0 if SLOPE
		endif
	while((!stringmatch(line,"*Data starts*")) && strlen(line)>0)	// check for start or EOF
	Close refNum

	if (bad==-2)												// alert for wrong type of data
		DoAlert 0, "an LTP file, but data type is '"+dataType+"', not SLOPE"
		print "an LTP file, but data type is '"+dataType+"', not SLOPE"
	endif
	bad = stringmatch(line,"*Data starts*") ? bad : -1			// error unless start of data
	return (bad  ? bad : i)										// return line number or error
End



Function spotErrorCircle(xx,angError,s2)
	Variable xx			// postion on mirror (m)
	Variable angError	// angular error  (rad)
	Variable s2			// distance from center of mirror to focus (m)
	angError = (angError==0) ? NaN : angError	// exactly zero indicates unfitted region
	return (s2-xx) * angError * 2				// ray deviation is twice error angle
End



Function mirrorAngleCircle(s1,s2,roc)	// compute mirror angle for a cylinder from roc
	Variable s1, s2					// distance from source & focus (m)
	Variable roc					// radious of curvature (m)
	// r/2 * sin(theta) = f,    and 1/f = 1/s1 + 1/s2
	Variable sine_theta = (2/roc) / (1/s1 + 1/s2)
	return asin(sine_theta)
End



Function FlipWave(wav)
	Wave wav
	Variable n=numpnts(wav)
	Duplicate wav, flip_wave_temp__
	wav = -flip_wave_temp__[n-1-p]
	KillWaves flip_wave_temp__
	String noteStr=note(wav)
	Variable iflip = !NumberByKey("flip",noteStr,"=")
	noteStr = ReplaceNumberByKey("flip", noteStr, iflip,"=")
	Note /K wav
	Note wav,noteStr
	return iflip
End



Function /T ListOfSlopeWaves()
	String big = WaveList("*",";","")
	String item,list=""
	Variable i=0
	do
		item = StringFromList(i, big)
		if (stringmatch(item, "*X" ))
		elseif (stringmatch(item, "*Height" ))
		elseif (stringmatch(item, "ww*" ))
		elseif (stringmatch(item, "*Covar" ))
		elseif (stringmatch(item, "eps*" ))
		elseif (stringmatch(item, "correlate*" ))
		elseif (stringmatch(item, "W_sigma*" ))
		elseif (stringmatch(item, "*_nm" ))
		elseif (stringmatch(item, "Res_*" ))
		elseif (stringmatch(item, "W_ParamConfidenceInterval*" ))
		elseif (stringmatch(item, "spotSize*" ))
		else
			list += item+";"
		endif
		i += 1
	while (strlen(item)>0)
	return list

//	String slope
//	Prompt slope, "choose wave with slope data", popup, list
//	DoPrompt "pick a wave", slope
//	return slope
End