#pragma rtGlobals=1		// Use modern global access method.
#pragma IgorVersion = 4.0
#pragma version = 3.1



Function FindDepositionProfile(slope,s1,s2,tth,p1,p2,optimizeType)
	String slope				// name of data wave with slope (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=long mirror, 0.06=short) (m)
	Variable tth				// twice the grazing angle on mirror (usually 0.006) (rad)
	Variable p1,p2			// points over which to do optimization and to consider when computing deposition
	Variable optimizeType	// type of optimization: 1=slope at 0,  2=total area,  3=max thickness  (usually 1)

	Variable length
	Variable partialRange	// flags whether to consider whole mirror or only a partial range for deposition (is it necessary to fix edges?)
	partialRange = (p2>p1) || (p1<0) || (p2<0) || (numtype(p1)==2) || (numtype(p2)==2)

	Variable prompt_slope,prompt_s1, prompt_s2, prompt_tth, prompt_optType
	prompt_slope = (exists(slope)!=1)
	prompt_s1 =  !( s1>=1)
	prompt_s2 =  !(s2>=.01 )
	prompt_tth =  !( tth>=1e-5 && tth<.02 )
	prompt_optType =  !( optimizeType>=1 && optimizeType<=3 )
	if (prompt_slope || prompt_s1 || prompt_s2 || prompt_tth || prompt_optType)
		if (prompt_slope)
			slope = StringFromList(0, WaveList("*Rad",";",""))
		endif
		Prompt slope, "name of wave with slope data",popup,WaveList("*Rad",";","")
		DoPrompt "choose profile to deposit upon", slope
		if (V_flag)
			DoAlert 0, "no data wave chosen, cannot compute deposition"
			return 1
		endif
		if (exists(slope)!=1)
			return 1										// nothing to be done
		endif
		Wave rad=$slope
		length = NumberByKey("length", note(rad),"=")
		s1 = (s1>=1) ? s1 : NumVarOrDefault("root:Packages:KBmirror:last_s1", 60)	// distance from Source to center of mirror (m)
		if (numtype(s2) || s2<.001)
			s2 = SelectNumber(stripeType(length)-1 ,0.100,0.130,0.060)// if bad s2, choose s2 to match stripe
		endif
		tth = (tth>=1e-5 && tth<.02) ? tth :  .006
		optimizeType = prompt_optType ? 1 : optimizeType
		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)", popup, s2list
		Prompt tth, "twice glancing angle on mirror (rad)"
		Prompt optimizeType, "object of minimization" popup, "slope at 0;total area;max thickness"
		Prompt partialRange, "deposit over the whole mirror" popup, "Yes;No"
		DoPrompt "make a deposition profile", s1,s2,tth,partialRange, optimizeType
		partialRange -= 1									// turn value into true/false flag
		if (V_flag)
			return 1
		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)"
			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)
		DoAlert 0, "cannot use s1="+num2str(s1)+"  or  s2="+num2str(s2)
		return 1
	endif
	NVAR last_s1 = root:Packages:KBmirror:last_s1  ;  last_s1 = s1
	NVAR last_s2 = root:Packages:KBmirror:last_s2  ;  last_s2 = s2

	Wave rad=$slope
	length = NumberByKey("length", note(rad),"=")
	String depositName		// name of wave with what to deposit
	String idealHeightName	// name of wave with ideal height

	Wave rad=$slope
	Wave radX = $StringByKey("Xdata", note(rad),"=")
	String smm = num2istr(round(length*1000))
	depositName = "deposit"+smm+"mm"
	idealHeightName = "idealHeight"+smm+"mm"

	if (numtype(p1) && numtype(p2) && (p1<p2))		// bad p1, p2 passed
		if (stripeType(length)==2)
			NVAR range0 = root:Packages:KBmirror:rangeS0	// fit range (in points) for short mirror
			NVAR range1 = root:Packages:KBmirror:rangeS1
		else
			NVAR range0 = root:Packages:KBmirror:rangeL0	// fit range (in points) for long mirror
			NVAR range1 = root:Packages:KBmirror:rangeL1
		endif
		p1=range0
		p2=range1
	endif
	if (optimizetype==2)
		Abort "optimizing by total area does not work yet"
	endif

	// set p1 and p2 range of interest for the deposition
	if (!partialRange)
		p1 = 0										// select the whole mirror
		p2 = numpnts(rad)-1
	else												//	put up graph to select important range
		PutUpGraph_GetRange(p1,p2,rad,radX,"deposit_range")
		// call to p1 and p2 is by reference, they are set by PutUpGraph_GetRange
	endif

	Wave height=$HeightFromLTP(rad)				// make the measured height profile using angles from LTP
	Duplicate/O height $idealHeightName
	Wave idealHeight=$idealHeightName
	String xunits = WaveUnits(radX,-1)
	Variable Xscaling = LengthScalng2meters(xunits)
	idealHeight = MirrorIdealHeight(s1,s2,tth,radX[p]*Xscaling)

	WaveStats/Q height
	Variable yoffset = 4*(V_max - V_min)			// yoffset big enough that: ideal - rotated(measuredHeight) > 0
	//	tilt the mirror to get the minimum amount of deposition

	// find the optimum tilt
	Duplicate /O radX optimize_xx__, optimize_temp__,  temp_rotate__X, temp_rotate__Y
	Duplicate /O height optimize_height__			// current measured height from LTP
	Duplicate /O idealHeight optimize_idealHeight__
	Make/N=1/O tiltWave_={yoffset,p1,p2,optimizeType}
	Optimize /A=0 /L=-0.001 /H=0.001 /T=1e-12 /I=100 /Q  Deposited_Optimize, tiltWave_
	Variable bestTilt = V_minloc

	// optimize_temp__ is now the properly tilted profile, but at too big an offset
	WaveStats /Q /R=[p1, p2] optimize_temp__	// V_min is new offset (for region of interest)
	optimize_temp__ -= V_min
	// optimize_temp__ is now the profile of what needs to be deposited
	Variable maxThicknes = V_max - V_min

	Xscaling = LengthScalng2meters(WaveUnits(radX,-1))
//	Variable xend = round(1000*abs(radX[inf] - radX[0])/2)/1000
//	Make/N=(2000*xend+1)/O $depositName
	// switch to 0.1 mm spacing, from 1mm spacing
	Variable scale10s = ceil(log(49/abs(radX[inf]-radX[0])))			// chooses scale to give at least 50 points
	scale10s = 10^scale10s
	Variable xend = round(scale10s*abs(radX[inf] - radX[0])/2)/scale10s
	Make/N=(2*scale10s*xend+1)/O $depositName

	Wave deposit=$depositName
	String newNote=note(radX)
	newNote = ReplaceStringByKey("Ydata",newNote,NameOfWave(rad),"=")
	newNote = ReplaceStringByKey("Xdata",newNote,NameOfWave(radX),"=")
	newNote = ReplaceStringByKey("length",newNote,num2str(length),"=")
	newNote = ReplaceStringByKey("s1",newNote,num2str(s1),"=")
	newNote = ReplaceStringByKey("s2",newNote,num2str(s2),"=")
	newNote = ReplaceStringByKey("p1",newNote,num2str(p1),"=")
	newNote = ReplaceStringByKey("p2",newNote,num2str(p2),"=")
	newNote = ReplaceStringByKey("x1",newNote,num2str(radX[p1]),"=")
	newNote = ReplaceStringByKey("x2",newNote,num2str(radX[p2]),"=")
	newNote = ReplaceStringByKey("optimizeType",newNote,num2str(optimizeType),"=")
	newNote = ReplaceStringByKey("maxThicknes",newNote,num2str(maxThicknes),"=")
	Note /K deposit
	Note deposit newNote
	Note /K idealHeight
	Note idealHeight newNote
	SetScale/I x (-xend),(xend),"m", deposit
	SetScale d 0,0,"m", deposit
	deposit = interp(x,radX,optimize_temp__)
	printf "optimized tilt = %.2f (µrad),   over range of interest is [%.1f, %.1f] mm\r",bestTilt*1.e6,1e3*radX[p1],1e3*radX[p2]
	printf "       deposited area of %.3g (m^2),  and a maximum thickness of %d (nm)\r",area(deposit,-xend,xend), maxThicknes*1e9
	KillWaves/Z optimize_xx__, optimize_idealHeight__, optimize_height__
	KillWaves/Z optimize_temp__, temp_rotate__X, temp_rotate__Y, tiltWave_

	DoAlert 1, "run check on resulting profile"
	if (V_flag==1)
		CheckDepositionProfile(radX,height,deposit)
	endif
	return 0
End

Function Deposited_Optimize(ww,tilt) : FitFunc
	//	this optimization function does one of three types of optimization
	//	slope at 0,  total area,  or max thickness
	Wave ww			// ww[0] = y offset to ensure ideal - rotated(measuredHeight) > 0
						// ww[1], ww[2] are points over which to do optimization
						// ww[3], type of optimization used  0=slope at 0,  1=total area,  2=max thickness
	Variable tilt
	// it is important that ww[0] is set so big that ideal is always > rotated(measuredHeight)
	// then afterwards excess offset is subtracted off

	Wave measuredHeight = optimize_height__	// y values of measuredHeight (at the xx)
	Wave ideal = optimize_idealHeight__		// y values of calculated perfect shape (at the xx)
	Wave xx = optimize_xx__					// x values for measuredHeight and ideal
	Wave temp=optimize_temp__				// measuredHeight tilted to optimal angle, + real big offset
	Wave rotate_X=temp_rotate__X				// work space used to properly rotate measuredHeight
	Wave rotate_Y=temp_rotate__Y

	Variable ct = cos(tilt)
	Variable st = sin(tilt)

	//	rotate 'measuredHeight',   rotation from (x0,y0) to (x,y)
	//	x = x0*cos(tilt) - y0*sin(tilt)
	//	y = x0*sin(tilt) + y0*cos(tilt)
	// make temp a rotated version of 'measuredHeight' (= optimize_height__)
	rotate_Y = xx[p]*st + measuredHeight[p]*ct	// rotate x and y values into new waves
	rotate_X = xx[p]*ct - measuredHeight[p]*st
	temp = interp(xx[p],rotate_X,rotate_Y)// interpolate rotated back to the original x grid

	// temp becomes difference between idel and rotated ww[0] big enough so temp[i] > 0
	temp = ideal[p] - temp[p] + ww[0]		// difference between ideal and rotated (+ and offset)

	// now choose the statistic that you want to minimize, and return that
	switch(ww[3])
		case 1:							// minimize slope at middle (=0)
			// to make the slope of the deposition zero at the center, so depositon rate is unchanged at middle
			Variable i = round(BinarySearchInterp(xx,0))				// middle of waves
			return abs((temp[i+2]-temp[i-2]) / (xx[i+2]-xx[i-2]))// return average slope in the middle
		case 2:							// minimize total area
			Variable x1 = round(BinarySearchInterp(xx,ww[1]))		// left edge of range
			Variable x2 = round(BinarySearchInterp(xx,ww[2]))		// right edge of range
			return faverageXY(xx,temp,x1,x2)							// does not work?
		case 3:							// minimize maximum thickness deposited
			WaveStats/Q/R=[ww[1],ww[2]] temp
			return V_max
	endswitch
End

Function/T Display_BothDepositProfiles(longName,shortName)
	String longName, shortName
	String list
	if (strlen(longName)<1 || strlen(shortName)<1)
		if (strlen(longName)<1)
			list = CheckWaveListForValues(WaveList("*Deposit*",";",""),"length",.082,.092)
			longName = StringFromList(0, list)
		endif
		if (strlen(longName)<1)
			longName = StringFromList(0, WaveList("*Deposit*",";",""))
		endif
		if (strlen(shortName)<1)
			list = CheckWaveListForValues(WaveList("*Deposit*",";",""),"length",.033,.041)
			shortName = StringFromList(0, list)
		endif
		Prompt longName, "pick name of long mirror deposition profile", popup, "none;"+WaveList("*Deposit*",";","")
		Prompt shortName, "pick name of short mirror deposition profile", popup, "none;"+WaveList("*Deposit*",";","")
		DoPrompt "deposition profiles to display", longName,shortName
		if (V_flag)
			return ""
		endif
	endif

	Variable i
	list = "slope at middle;total area;max thickness"
	Display
	DoWindow /C $UniqueName("GraphDeposit",6,0)
	SetWindow kwTopWin note="deposition"
	if (!stringmatch(longName,"none"))
		AppendToGraph $longName
		i = NumberByKey("optimizeType",note($longName),"=")
		TextBox/C/N=text0/F=0 longName+" minimizes "+StringFromList(i-1, list)
	endif
	if (!stringmatch(shortName,"none"))
		AppendToGraph $shortName
		ModifyGraph/Z rgb($shortName)=(1,4,52428)
		i = NumberByKey("optimizeType",note($shortName),"=")
		TextBox/C/N=text0/F=0 shortName+" minimizes "+StringFromList(i-1, list)
	endif
	ModifyGraph/Z gfSize=14, lowTrip(bottom)=0.1, lowTrip(left)=0.001
	ModifyGraph/Z lSize=2, tick=2, mirror=1, minor=1, standoff=0
	Label left "depth to deposit  (\\U)"
	Label bottom "position along mirror  (\\U)"
	ModifyGraph zero(left)=2

	Variable x1,x2
	if (!stringmatch(longName,"none"))
		x1 = NumberByKey("x1",note($longName),"=")
		x2 = NumberByKey("x2",note($longName),"=")
		SetDrawLayer UserFront
		SetDrawEnv xcoord= bottom,ycoord= prel,linefgc= (65535,0,0),dash= 1
		DrawLine x1,0,x1,1
		SetDrawEnv xcoord= bottom,ycoord= prel,linefgc= (65535,0,0),dash= 1
		DrawLine x2,0,x2,1
	endif
	if (!stringmatch(shortName,"none"))
		x1 = NumberByKey("x1",note($shortName),"=")
		x2 = NumberByKey("x2",note($shortName),"=")
		SetDrawLayer UserFront
		SetDrawEnv xcoord= bottom,ycoord= prel,linefgc= (1,4,52428),dash= 1
		DrawLine x1,0,x1,1
		SetDrawEnv xcoord= bottom,ycoord= prel,linefgc= (1,4,52428),dash= 1
		DrawLine x2,0,x2,1
	endif
	return WinName(0,1)
End




Function CheckDepositionProfile(radX,height,deposit)
	Wave radX				// x-positions of measurement
	Wave height				// height of mirror surface determined from slopes (m)
	Wave deposit			// what we want to deposit

	if (!(WaveExists(radX) && WaveExists(height)&& WaveExists(deposit)))
		String radXname, heightName,depositName
		Prompt radXname,"name of slope Xdata", popup, WaveList("*RadX",";","")
		Prompt heightName,"name of height (measured) wave", popup, WaveList("height*",";","")
		Prompt depositName,"name deposition profile", popup, WaveList("deposit*",";","")
		DoPrompt "pick waves to check deposition", radXname, heightName,depositName
		if (V_flag)
			return 1
		endif
		Wave radX = $radXname
		Wave height = $heightName
		Wave deposit = $depositName
	endif
	if (!(WaveExists(radX) && WaveExists(height)&& WaveExists(deposit)))
		DoAlert 0, "Needed waves not found, returning"
		return 1
	endif

	String spotError, sCorrelate = ""
	Variable s1, s2

	Variable length=NumberByKey("length",note(radX),"=")
	String mm = num2istr(round(length*1000))

	String sCorrected = "corrected"+mm+"mm"
	String sCorrected_deriv = "corrected"+mm+"mm_deriv"
	Duplicate/O radX $sCorrected
	Wave corrected = $sCorrected
	SetScale d 0,0,"m", corrected
	corrected = height + deposit(radX[p])
	DifferentiateXY(radX, corrected, sCorrected_deriv)
	Wave corrected_deriv = $sCorrected_deriv
	if (stringmatch(WaveUnits(radX,1),"m")&&stringmatch(WaveUnits(corrected,1),"m"))
		SetScale d 0,0,"rad", $sCorrected_deriv
	endif
	Note/K corrected_deriv
	Note corrected_deriv, note(deposit)
	s2 = NumberByKey("s2",note(deposit),"=")
	s1 = NumberByKey("s1",note(deposit),"=")
	String windowID = "check"+mm+"mm"
	Fit_a_Stripe(sCorrected_deriv,NameOfWave(radX),s1,s2,3,windowID)
	if (exists("M_Covar"))
		Wave M_Covar=M_Covar
		sCorrelate="correlate"+mm+"mm"
		Duplicate/O M_Covar, $sCorrelate				// make the correlation matrix
		Wave corelate=$sCorrelate
		corelate = M_Covar[p][q]/sqrt(M_Covar[p][p]*M_Covar[q][q])
		KillWaves M_Covar
	endif
	// make the spot error
	spotError = "spotSize"+mm+"mmCheck"
	String sres = "Res_"+sCorrected_deriv
	SetScale d 0,0,"rad", $sres
	Duplicate/O $sres $spotError
	Wave res = $sres
	Wave spot = $spotError
	spot = spotErrorCircle(radX[p],res[p],s2)
	SetScale d 0,0,"m", $spotError

	DisplayOneFit("fit_"+sCorrected_deriv,sres,spotError,sCorrelate,windowID)
	SetWindow kwTopWin note="check"
End





Function CallMakeDepositionForMetrology(deposit)
	String deposit
	if (strlen(deposit)<1)
		Prompt deposit "deposition profile to write", popup,WaveList("deposit*mm",";","")
		DoPrompt "pick deposition profile",deposit
		if (V_flag)
			return 1
		endif
	endif
	if (exists(deposit)==1)
		MakeDepositionForMetrology($deposit)
		printf "saved file %s%s\r", getdatafolder(1),deposit
	else
		Abort "deposition profile named '"+deposit+"' does not exist"
	endif
End
Function MakeDepositionForMetrology(deposit)
	Wave deposit
	if (LengthScalng2meters(WaveUnits(deposit,-1))!=1)
		Abort "position along mirrors must be in meters, not mm"
	endif

	String list = note(deposit)
	String wname=NameOfWave(deposit)+"_nm"
	Duplicate/O $NameOfWave(deposit) $wname
	Wave nm=$wname
	nm *= 1e9
	Variable xlo = 1000 * DimOffset(deposit, 0)
	Variable xhi = 1000* (DimDelta(deposit, 0) * (numpnts(deposit)-1) + DimOffset(deposit, 0))
	SetScale/I x (xlo),(xhi),"mm", nm
	SetScale d 0,0,"nm", nm

	Variable flip = NumberByKey("flip", list, "=")
	if (flip)
		Variable N=numpnts(nm)
		Variable temp
		SetScale/I x (-xhi),(-xlo),"mm", nm				// flip the curve nm
		Duplicate/O nm, flip_wave_temp__
		nm = flip_wave_temp__[N-1-p]
		KillWaves flip_wave_temp__
	endif

	String fileName=StringByKey("fileName",list,"=")		// original slope data file
	Variable i = strsearch(fileName,".",0)
	fileName =  fileName[0,i]+"nm"
	Variable refNum
	Open /C="R*ch" /M="File for Deposition Profile" /T="TEXT" refNum as fileName
	if (strlen(S_fileName)<1)
		Abort
	endif
	fprintf refNum,"# thickness to deposit on a mirror blank to make it elliptical\r"
	Variable seconds=NumberByKey("MODTIME",WaveInfo(deposit, 0))
	fprintf refNum,"Date_Computed=%s %s\r", Secs2Date(seconds,2), Secs2Time(seconds,1)
	fprintf refNum,"Date_Measured=%s\r"  , StringByKey("dateTaken",list,"=")
	fprintf refNum,"Slope_File_Name=%s\r"  , StringByKey("fileName",list,"=")
	fprintf refNum,"Igor_location=%s\r"  , StringByKey("path",list,"=")
	if (flip)
		fprintf refNum,"need to flip for x-rays (only)\r"
	endif
	fprintf refNum,"region_lo_mm=%.2f\r"  , NumberByKey("x1",list,"=")*1e3
	fprintf refNum,"region_hi_mm=%.2f\r"  , NumberByKey("x2",list,"=")*1e3
	fprintf refNum,"source_dist_m=%.3f\r"  , NumberByKey("s1",list,"=")
	fprintf refNum,"object_dist_m=%.3f\r"  , NumberByKey("s2",list,"=")
	i = NumberByKey("optimizeType",list,"=")
	fprintf refNum,"optimization_type=%s\r"  , StringFromList(i-1, "slope at middle;total area;max thickness")
	fprintf refNum,"\rposition_mm\theight_nm\r"
	Close refNum
	Duplicate/O nm, wave_nm_temp__
	wave_nm_temp__ = x
	Save /A=2/J wave_nm_temp__,nm as S_fileName
	KillWaves wave_nm_temp__
End