#pragma rtGlobals=1		// Use modern global access method.
#pragma IgorVersion = 4.0
#pragma version = 1.0
#include  "Xray", version>=2.2



// ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
//
// All values are in nm  (that's nanometer)!    no �, no �m, no cm, etc.
//
// ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
//
Static Constant re = 2.81794E-6		// electron radius in (nm)
Static Constant hc = 1.239841757		// keV-nm  (12keV is ~0.1nm)


Menu "multi-layer"
	"update",/Q, reCalcLayersProc("",0,"","")
	"re-set to defaults",SetDefaultsProc("")
	"table of layer definitions",/Q, MakeTableLayers()
End


// definitions:
// names are just a short name to be used in the stack diagram (not really needed)
// a is lattice constant used to get Vc and thickness is n*a  (nm)
// n is number of lattice planes of lattice constant a
//		the thickness of layer j is a[j]*n[j]   (nm)
// noteStr, a wave note to be added to the resultant reflectivity wave
// planes, this contains definition of the layered material, the format of planes is:
//	formula1A:frac1A;formula1B:frac1B		where 1A and 1B are coexisting materials (frac1A+frac1B�1)
//	formula2A:frac2A;formula2B:frac2B		where 2A and 2B are coexisting materials (frac2A+frac2B�1)
//	formula3A:frac3A							where the concentration of 3A is frac3A
//	formula4A									one solid layer of 4A
//	formula5A:frac5A;formula5B:frac5B		where 5A and 5B are coexisting materials (frac5A+frac5B�1)
//	... etc.
//
// an example of layers.planes[] is:
//	SrTiO3										// full STO
//	SrTiO3:0.9;SrRuO3:0.1						// mixed 90%STO and 10%SRO
//	SrTiO3TiO2:0.1;SrRuO3:0.89				// mixed 10%STO and 89% SRO (not sum is not 100%)
//	SrRuO3										// full SRO
//	SrO:0.7;PbO:0.3							// a mixed layer of SrO and PbO
//	PbO:Zr0.2Ti0.8O2							// PZT
//	PbO:Zr0.2Ti0.8O2.5:0.35					// a 35% partially filled layer of PZT
//
Structure ref_layer_structure
	Wave/T planes				// formula of layer, store as a list (not usual formula)
	Wave/T names				// names to use for labels
	Wave a						// lattice constant of each layer
	Wave n						// number of cells in each layer
	String noteStr				// string to use for the wave note
EndStructure
//
Window Table_layers() : Table
	MakeTableLayers()
EndMacro
Function MakeTableLayers() : Table
	String win="TableLayers"
	if (strlen(WinList(win,"","WIN:2" )))
		DoWindow /F $win
		return 0
	endif
	Edit/K=1/W=(591,351,1095,696) :layerDefnition:names,:layerDefnition:planes,:layerDefnition:a,:layerDefnition:n
	DoWindow /C $win
	ModifyTable width(Point)=32,width(:layerDefnition:names)=50,alignment(:layerDefnition:planes)=0,width(:layerDefnition:planes)=312
	ModifyTable width(:layerDefnition:a)=52,width(:layerDefnition:n)=40
End
//
Function totalThickness(layers)			// total thickness of all finite layers
	STRUCT ref_layer_structure &layers
	String wName = UniqueName("thick",1,0)
	Make/N=(numpnts(layers.n))/D $wName
	Wave dt=$wName
	dt = layers.a[p] * layers.n[p]
	dt = numtype(dt[p]) ? 0 : dt[p]
	Variable thickness=sum(dt)
	KIllWaves/Z dt
	return thickness
End

Function EnsureLayersStructExists()	// make sure waves exist to fill the layers structure
	NewDataFolder/O layerDefnition				// put the waves needed for the layer definition here

	Wave planes = :layerDefnition:planes
	if (exists(":layerDefnition:planes")!=1)
		Make/N=(2)/O/T :layerDefnition:planes
	elseif (WaveType(planes))
		Make/N=(2)/O/T :layerDefnition:planes		// force it to be a text wave
	endif
	Wave names = :layerDefnition:names
	if (exists(":layerDefnition:names")!=1)
		Make/N=(2)/O/T :layerDefnition:names
	elseif (WaveType(names))
		Make/N=(2)/O/T :layerDefnition:names		// force it to be a text wave
	endif
	Wave a = :layerDefnition:a
	if (exists(":layerDefnition:a")!=1)
		Make/N=(2)/O/D :layerDefnition:a
	elseif (WaveType(a)!=4)
		Make/N=(2)/O/D :layerDefnition:a			// force it to be a double
	endif
	Wave n = :layerDefnition:n
	if (exists(":layerDefnition:n")!=1)
		Make/N=(2)/O/D :layerDefnition:n
	elseif (WaveType(n)!=4)
		Make/N=(2)/O/D :layerDefnition:n			// force it to be a double
	endif
End


Function reCalcLayersProc(ctrlName,varNum,varStr,varName) : SetVariableControl
	String ctrlName
	Variable varNum
	String varStr
	String varName

	String str = "reCalcLayers"
	str = SelectString(exists(str)<5,str,"reCalcLayers_proto")
	FUNCREF reCalcLayers_proto func = $str
	return func()
End
Function reCalcLayers_proto()
	DoAlert 0, "reCalcLayers() not found"
	return 1
End

Function SetDefaultsProc(ctrlName) : ButtonControl
	String ctrlName

	String str = "SetDefaults"
	str = SelectString(exists(str)<5,str,"SetDefaults_proto")
	FUNCREF SetDefaults_proto SetDefaults = $str
	SetDefaults()
End
Function SetDefaults_proto()
	DoAlert 1, "SetDefaults() not found, nothing done"
	return 1
End



//	*********************************************************************************
//	***************************** start of calculating reflectivity ******************************

Function calcMultiLayerReflectivity(lambda,layers,reflect)
	Variable lambda								// wavelength in vacuum (nm)
	STRUCT ref_layer_structure &layers		// the structure that holds the layers
	Wave reflect								// wave to calculate, assumed to be scaled in 1/nm

	Variable N=numpnts(layers.planes)
	Variable keV=hc/lambda						// photon energy (keV)
	Variable Vc
	Variable/C k0 = 2*PI/lambda				// 2�/lambda
	Variable/C F, one=cmplx(1,0)
	Make/N=(N+1)/O d_nm
	Make/N=(N+1)/C/O k_multi
	d_nm[0] = Inf
	d_nm[1,N] = layers.a[N-p]*layers.n[N-p]
	d_nm = numtype(d_nm) ? 0 : d_nm
	k_multi[0] = k0
	Variable j
	for (j=1;j<=N;j+=1)
		F = getF000(keV,layers.planes[N-j],layers.a[N-j])
		Vc = (layers.a[N-j])^3
		k_multi[j] = k0*(one - 0.5*(F/Vc)*re*lambda^2/PI)
	endfor
	if (!strsearch(GetRTStackInfo(0),"test",0))
		printWave(d_nm)
		print complexVecStr(k_multi)
	endif

	Variable Qmin = leftx(reflect), Qmax = rightx(reflect)-deltax(reflect)
	String xlabel=GetDimLabel(reflect,0,-1)
	SetScale/I x asin(Qmin*lambda/4/PI),asin(Qmax*lambda/4/PI),"rad", reflect	// max angle comes from max Q
	reflect = magsqr(multiLayerAmplitude(k_multi,d_nm,x))	// This is the calculation of reflectivity
	SetScale/I x Qmin,Qmax,xlabel, reflect		// reset scaling from angle back to Q
	KillWaves/Z d_nm,k_multi
End


// k[j] = 2*�*n[j](1+i*kappa[j])/lambda, where n is real index, kappa is imag correction, omega/c =2�/lambda0
//   note, in this routine, k is for the ray, it is NOT the Q vector (=kf-ki)
// d[j] = thickness of each layer (nm),  note that for last layer which is infinite, use d[N]=Inf, use d[0]=anything
// Bragg = angle between incident ray and surface (just like Bragg angle) (radians)
// length of k[] and d[] are N
//
// This routine has been changed to properly account for U (not just A, R, & T)

Static Function/C multiLayerAmplitude(k,d,Bragg)
	Wave/C k						//k-vector of incident wave for every layer k = 2*�*n(1+i*kappa)/lambda (1/nm)
	Wave d							// thickness of each layer (nm)
	Variable Bragg					// angle between incident ray and surface (radians)

//	NVAR once=once
	Variable N=numpnts(d)
	String wName = UniqueName("phi",1,0)
	Make/N=(N)/C/D $wName
	Wave/C/D sine = $wName
	sine[0] = cmplx(sin(PI/2-Bragg),0)
	sine[1,N] = k[p-1]/k[p] * sine[p-1]	// complex Snell's law, note that k1/k2 = n1/n2, since k = 2�n/lambda
	// convert sine[] to phi[]
	Wave/C/D phi = $wName				// note that I am re-using the same wave to be called both sine and phi!
	phi = asin(sine)						// I can do this because I will not need the sin(phi) any more
//	if (once)
//		printf "***  multiLayerAmplitude(k,d,Bragg=%g)\r",Bragg*180/PI
//		phi *= 180/PI
//		print complexVecStr(phi),"�","   angle to the normal"
//		phi *= PI/180
//	endif

	Variable/C q							// q = n2/n1 * cos(phi2) / cos(phi1)
	Variable/C cosine,one=cmplx(1,0),i=cmplx(0,1), eikr
	Variable/C A=one,R=cmplx(0,0),T,U	// amplitudes A=incident, R=reflected, T=transmitted, U=upward
	Variable j
	for (j=N-1;j>0;j-=1)					// calculate for interface between j and j-1
		cosine = cos(phi[j])
		eikr = exp( i*k[j]*d[j]*cosine )
		T =  A*eikr							// A is last one from j+1, propogate backward through layer j to get T
		U = R/eikr							// R is last one from j+1, propogate forward through layer j to get U
		q = k[j]/k[j-1] * cosine/cos(phi[j-1])
		A = 0.5 * ( (one+q)*T + (one-q)*U )	// A for j
		R = U+T-A							// R for j
//		if (once)
//			printf "interface %d%d,  phi={%s, %s},	A=%s		R=%s 		T=%s		U=%s		with q=%s",j-1,j,polarStr(phi[j-1]*180/PI),polarStr(phi[j]*180/PI),polarStr(A),polarStr(R),polarStr(T),polarStr(U),polarStr(q)
//			if (magsqr((T+U)-(A+R))>1e-20)
//				printf ",   (T+U)-(A+R)=%s",polarStr((T+U) - (A+R))
//			endif
//			printf "\r"
//			//	printf "phase factors are 1/exp(ikd cos) = %s and     exp(ikd cos) = %s\r",polarStr(one/eikr),polarStr(eikr)
//		endif
	endfor
	KillWaves/Z phi
	// the transmitted amplitude at the bottom is 1/A
	Variable/C ref = R/A					// the amplitude reflectivity
//	if (once)
//		print  "final amplitude reflectivity =",polarStr(ref)
//	endif
	return ref
End

Static Function/T complexVecStr(wz)
	Wave/C wz
	Variable j
	String str, out

	sprintf out,"%s = {%s",NameOfWave(wz),polarStr(wz[0])
	for (j=1;j<numpnts(wz);j+=1)
		sprintf str,", %s",polarStr(wz[j])
		out += str
	endfor
	out += "}"
	return out
End
//
Static Function/T polarStr(z)
	Variable/C z
	Variable angle = atan2(imag(z),real(z))*180/PI
	String str
	if (abs(angle)<1e-3)
		sprintf str, "%.5g",sqrt(magsqr(z))
	elseif (abs(180-abs(angle))<1e-3)
		sprintf str, "%.5g",-sqrt(magsqr(z))
	else
		sprintf str, "%.5g e^(%.1f�i)",sqrt(magsqr(z)),atan2(imag(z),real(z))*180/PI
	endif
	return str
End


// 'list' is a string like "PbO,0:Zr0.2Ti0.8O2,0.5:1"	 "formula,z/a:formula,z/a:formula,z/a:multiplier"
Static Function/C getF000(keV,list,a)	// structure factor at Q=0
	Variable keV			// photon energy (keV)
	String list				// "formula,z/a:formula,z/a:formula,z/a:multiplier", e.g. "TiO2,0:Sr0,0.5:1;"
	Variable a				// lattice constant (nm)

	Variable/C i=cmplx(0,1)
	Variable/C Fm						// F for a single material, i.e. for F for SrTiO3
	Variable/C F=cmplx(0,0)			// the whole F for this layer is accumulated in this
	Variable N=ItemsInList(list), x
	Variable m, j
	String item							// "formula,z/a,formula,z/a,multiplier", e.g. "TiO2,0:Sr0,0.5:1", 
	String sub							// just the "formula,z/a" part, e.g. "Sr0,0.5"
	String formula						// chemical formula, e.g. "Ti;O,2"
	for (j=0;j<N;j+=1)
		item = StringFromList(j,list)	// this is something like  "TiO2,0:Sr0,0.5:1", one material
		Fm = cmplx(0,0)
		m =0
		do
			sub = StringFromList(m,item,":")		// something like "SrO,0.5"
			formula  = ChemicalFormulaToList(StringFromList(0,sub,","))	// Something like "TiO2"
			if (strlen(formula)>0)
				Fm += F000ofFormula(formula,keV)
			elseif (numtype(str2num(sub))==0)	// the multiplier (which is always last)
				Fm *= cmplx(str2num(sub),0)
				break
			else										// nothing left
				break
			endif
			m += 1
		while(1)
		F += Fm
	endfor
	return F
End
//
Static Function/C F000ofFormula(formula,keV)	// return F of something like "Ti;O,2;", no phases in here, Q=0
	String formula
	Variable keV			// photon energy (keV)
	if (exists("Get_f")!=6)
		Abort "Load the 'X-ray Data' Package, and try again (See the Analysis Menu)"
	endif

//	Variable keV = NumVarOrDefault("root:keV",NaN)
	String item,atom
	Variable m									// number of atoms of one element
	Variable i,x
	Variable/C F=cmplx(0,0)
	FUNCREF Get_f_proto fatom = $"Get_f"
	i = 0
	do
		item = StringFromList(i,formula)		// something like "Ti" or "O,2"
		atom = StringFromList(0,item,",")
		if (strlen(atom)<1)
			break
		endif
		m= str2num(StringFromList(1,item,","))
		m = numtype(m) ? 1 : m				// if no m, default to 1
		F += fatom(atom,0, keV)*m
		i += 1
	while(1)
	return F
End
Function/C Get_f_proto(atom,Q,keV)
	String atom
	Variable Q, keV
	return NaN
End

//	****************************** end of calculating reflectivity ******************************
//	*********************************************************************************



//	*********************************************************************************
//	******************************** start of drawing boxes *********************************

Constant layerBoxheight=15			// height of one layer on DisplayLayers

Function ReDrawAllLayers(layers)
	STRUCT ref_layer_structure &layers	// description of all of the layers

	String win = "DisplayLayers"
	if (strlen(WinList(win,"","WIN:1" )))
		//	DoWindow /F $win		// don't need to bring it to the front
	else
//		Display /W=(4,50,305,540)/K=1
		Display /W=(4,50,305,712)/K=1 
		DoWindow /C $win
	endif

	GetWindow $win gsize					// get V_top, V_bottom, V_left, V_bottom
	SetDrawLayer/W=$win/K UserFront
	Variable xc = 118
	Variable y0=V_bottom-5				// start 5 up from bottom

	Variable N = numpnts(layers.n)			// number of layers in structure
	Variable j,m, i
	String material
	String str1,str2,form1,form2
	String form3,form4
	Variable x
	Variable imaterials, ilayers
	String str = "fixUpLabels"
	str = SelectString(exists(str)<5,str,"fixUpLabels_proto")
	FUNCREF fixUpLabels_proto fixUpLabels = $str

	for (j=0;j<N;j+=1)
		m = layers.n[j]
		if (numtype(m))
			y0 = DrawDots(win,xc,y0,5)	// draw big dots for the thick substrate
		endif
		imaterials = ItemsInList(layers.planes[j])
		ilayers = ItemsInList(StringFromList(0,layers.planes[j]),":")

		if (imaterials==1 && ilayers==1 && m>1)			// full single-layers
			str1 = StringFromList(0,layers.planes[j],":")	// something like "SrTiO3"
			form1 = StringFromList(0,str1,",")
			form1 = fixUpLabels(form1)
			y0 = DrawNsingleLayers(win,m,y0,xc,form1,layers.names[j])
		elseif (imaterials==1 && ilayers==2)					// full bi-layer
			str1 = StringFromList(1,layers.planes[j],":")	// something like "TiO2,0"
			str2 = StringFromList(0,layers.planes[j],":")	// and "SrO,0.5"
			form1 = StringFromList(0,str1,",")
			form2 = StringFromList(0,str2,",")
			form1 = fixUpLabels(form1)
			form2 = fixUpLabels(form2)
			y0 = DrawNbiLayers(win,m,y0,xc,form1,form2,layers.names[j])
		elseif (imaterials==1 && ilayers==1)				// single-layer
			str1 = StringFromList(0,layers.planes[j],":")	// something like "TiO2,0"
			form1 = StringFromList(0,str1,",")
			form1 = fixUpLabels(form1)
			y0 = Draw1Layer(win,y0,xc,form1)
		elseif (imaterials==1 && ilayers==3)				// partial bi-layer
			str1 = StringFromList(1,layers.planes[j],":")	// something like "TiO2,0"
			str2 = StringFromList(0,layers.planes[j],":")	// and "SrO,0.5"
			form1 = StringFromList(0,str1,",")
			form2 = StringFromList(0,str2,",")
			form1 = fixUpLabels(form1)
			form2 = fixUpLabels(form2)
			x = str2num(StringFromList(2,layers.planes[j],":"))	// and "SrO,0.5"
			y0 = DrawMixedLayers(win,y0,xc,x,form1,"")
			y0 = DrawMixedLayers(win,y0,xc,x,form2,"")
		elseif (imaterials==2 && ilayers==3)				// two partial bi-layer
			material = StringFromList(0,layers.planes[j])	// first partial bi-layer
			str1 = StringFromList(1,material,":")			// something like "TiO2,0"
			str2 = StringFromList(0,material,":")			// and "SrO,0.5"
			form1 = StringFromList(0,str1,",")
			form2 = StringFromList(0,str2,",")
			form1 = fixUpLabels(form1)
			form2 = fixUpLabels(form2)
			x = str2num(StringFromList(2,material,":"))	// and "SrO,0.5"

			material = StringFromList(1,layers.planes[j])
			str1 = StringFromList(1,material,":")			// something like "TiO2,0"
			str2 = StringFromList(0,material,":")			// and "SrO,0.5"
			form3 = StringFromList(0,str1,",")
			form4 = StringFromList(0,str2,",")
			form3 = fixUpLabels(form3)
			form4 = fixUpLabels(form4)
			y0 = DrawMixedLayers(win,y0,xc,x,form1,form3)
			y0 = DrawMixedLayers(win,y0,xc,x,form2,form4)

		elseif (imaterials==2 && ilayers==2)				// a mixed single layer
			material = StringFromList(0,layers.planes[j])	// first partial bi-layer
			str1 = StringFromList(0,material,":")			// something like "SrO,0"
			form1 = StringFromList(0,str1,",")
			form1 = fixUpLabels(form1)
			x = str2num(StringFromList(1,material,":"))

			material = StringFromList(1,layers.planes[j])	// first partial bi-layer
			str2 = StringFromList(0,material,":")			// something like "SrO,0"
			form2 = StringFromList(0,str2,",")
			form2 = fixUpLabels(form2)
			y0 = DrawMixedLayers(win,y0,xc,x,form1,form2)
		else
			str1 = layers.planes[j]
			printf "planes[j] = '%s',   imaterials=%d,   ilayers=%d\r",str1,imaterials,ilayers
		endif
	endfor
	return y0
End
Function/T fixUpLabels_proto(str)
	String str
	// str = ReplaceString("Zr0.2Ti0.8O2",str, "Zr,Ti O2")
	// more ReplaceString() commands can go here
	return str
End
Function/T formula2rgb_proto(formula)
	String formula
	return "65535,65535,65535"		// default to white
End


Static Function DrawNsingleLayers(win,N,y0,xc,str,material)
	String win					// name of window
	Variable N					// number of bi-layer
	Variable y0					// bottom of first layer (=451)
	Variable xc					// horizontal center (=118)
	String str					// label for material
	String material				// name of material

	Variable y0start = y0		// save the starting y
	Variable hw = 84			// half width of the box
	Variable i
	Variable useInf=0			// flags that we asked for an infinite number of layers

	if (N<1)
		return y0
	elseif (numtype(N)==1)
		N = 3					// was 3, or 7
		useInf = 1
	elseif (N>1)
		SetDrawEnv/W=$win xcoord= abs,ycoord= abs
		DrawLine/W=$win xc+hw+4,y0,xc+2*hw,y0
	endif

	if (N<5)					// draw all of them		(was 5)
		for (i=0;i<N;i+=1)
			y0 = Draw1Layer(win,y0,xc,str)
		endfor
	else
		for (i=0;i<2;i+=1)
			y0 = Draw1Layer(win,y0,xc,str)	// draw single-layers before the "..."
		endfor
		y0 = DrawDots(win,xc,y0,3)			// put in the "..."
		for (i=0;i<2;i+=1)
			y0 = Draw1Layer(win,y0,xc,str)	// draw single-layers after the "..."
		endfor
	endif
	if (N>1)
		SetDrawEnv/W=$win xcoord= abs,ycoord= abs
		DrawLine/W=$win xc+hw+4,y0,xc+2*hw,y0
	endif
	SetDrawEnv/W=$win xcoord= abs,ycoord= abs,textxjust= 1,textyjust= 1
	DrawText/W=$win xc+1.5*hw+2,(y0+y0start)/2,SelectString(useInf,num2str(N)+" of "+material,material+" substrate")
	return y0
End


Static Function DrawNbiLayers(win,N,y0,xc,str1,str2,material)
	String win					// name of window
	Variable N					// number of bi-layer
	Variable y0					// bottom of first layer (=451)
	Variable xc					// horizontal center (=118)
	String str1, str2			// names for each of the bi-layers
	String material				// name of material

	Variable y0start = y0		// save the starting y
	Variable hw = 84			// half width of the box
	Variable i
	Variable useInf=0			// flags that we asked for an infinite number of layers

	if (N<1)
		return y0
	elseif (numtype(N)==1)
		N = 3					// was 3, or 7
		useInf = 1
	elseif (N>1)
		SetDrawEnv/W=$win xcoord= abs,ycoord= abs
		DrawLine/W=$win xc+hw+4,y0,xc+2*hw,y0
	endif

	if (N<5)					// draw all of them		(was 5)
		for (i=0;i<N;i+=1)
			y0 = Draw1Layer(win,y0,xc,str1)
			y0 = Draw1Layer(win,y0,xc,str2)
		endfor
	else
		for (i=0;i<2;i+=1)
			y0 = Draw1Layer(win,y0,xc,str1)	// draw bi-layers before the "..."
			y0 = Draw1Layer(win,y0,xc,str2)
		endfor
		y0 = DrawDots(win,xc,y0,3)			// put in the "..."
		for (i=0;i<2;i+=1)
			y0 = Draw1Layer(win,y0,xc,str1)	// draw bi-layers after the "..."
			y0 = Draw1Layer(win,y0,xc,str2)
		endfor
	endif
	if (N>1)
		SetDrawEnv/W=$win xcoord= abs,ycoord= abs
		DrawLine/W=$win xc+hw+4,y0,xc+2*hw,y0
	endif
	SetDrawEnv/W=$win xcoord= abs,ycoord= abs,textxjust= 1,textyjust= 1
	DrawText/W=$win xc+1.5*hw+2,(y0+y0start)/2,SelectString(useInf,num2str(N)+" of "+material,material+" substrate")
	return y0
End
//
Static Function DrawMixedLayers(win,y0,xc,f,str1,str2)
	String win					// name of window
	Variable y0					// bottom of first layer
	Variable xc					// x center of box
	Variable f					// f*a + (1-f)*b
	String str1,str2			// name to put in box
	Variable hw = 84			// half width of the box
	String rgb

	String str = "fixUpLabels"
	str = SelectString(exists(str)<5,str,"fixUpLabels_proto")
	FUNCREF fixUpLabels_proto fixUpLabels = $str
	str = "formula2rgb"
	str = SelectString(exists(str)<5,str,"formula2rgb_proto")
	FUNCREF formula2rgb_proto formula2rgb = $str
	String cmd
	//	Variable xm =  xc-hw + (1-f)*2*hw	// a is [xc-hw,xm],   b is [xm,xc+hw]
	Variable xm =  xc-hw + f*2*hw		// a is [xc-hw,xm],   b is [xm,xc+hw]
	if (f>0)
		rgb = formula2rgb(str1)
		sprintf cmd, "SetDrawEnv/W=%s xcoord= abs,ycoord= abs,linethick= 0,fillfgc=(%s)",win,rgb
		Execute cmd
		DrawRect/W=$win xc-hw,y0-layerBoxheight,xm,y0
	endif
	if (f<1)
		rgb = formula2rgb(str2)
		sprintf cmd, "SetDrawEnv/W=%s xcoord= abs,ycoord= abs,linethick= 0,fillfgc=(%s)",win, rgb
		Execute cmd
		DrawRect/W=$win xm,y0-layerBoxheight,xc+hw,y0
	endif
	if (f>0)
		str1 = fixUpLabels(str1)
		SetDrawEnv/W=$win xcoord= abs,ycoord= abs,textxjust= 1,textyjust= 1		// put label in center of region
		DrawText/W=$win (xc-hw+xm)/2,y0-7,str1
	endif
	if (f<1)
		str2 = fixUpLabels(str2)
		SetDrawEnv/W=$win xcoord= abs,ycoord= abs,textxjust= 1,textyjust= 1		// put label in center of region
		DrawText/W=$win (xm+xc+hw)/2,y0-7,str2
	endif
	return (y0-layerBoxheight)
End
//
Static Function Draw1Layer(win,y0,xc,str)
	String win					// name of window
	 Variable y0				// bottom of first layer
	Variable xc					// x center of box
	String str					// name to put in box

	Variable hw = 84			// half width of the box
	String fstr = "formula2rgb"
	fstr = SelectString(exists(fstr)<5,fstr,"formula2rgb_proto")
	FUNCREF formula2rgb_proto formula2rgb = $fstr
	String cmd, rgb = formula2rgb(str)
	sprintf cmd, "SetDrawEnv/W=%s xcoord= abs,ycoord= abs,linethick= 0,fillfgc=(%s)",win,rgb
	Execute cmd
	DrawRect/W=$win xc-hw,y0-layerBoxheight,xc+hw,y0
	SetDrawEnv/W=$win xcoord= abs,ycoord= abs,textxjust= 1,textyjust= 1
	DrawText/W=$win xc,y0-7,str
	return (y0-layerBoxheight)
End
//
Static Function DrawDots(win,xc,y0,nhw)		// use xc=18, y0=490
	String win					// name of window
	Variable xc, y0
	Variable nhw
	Variable hw = nhw<4 ? 2 : 3
	y0-= 2*hw
	SetDrawEnv/W=$win xcoord= abs,ycoord= abs,fillfgc= (0,0,0)
	DrawOval/W=$win xc-hw,y0-hw,xc+hw,y0+hw
	y0 -= nhw*hw
	SetDrawEnv/W=$win xcoord= abs,ycoord= abs,fillfgc= (0,0,0)
	DrawOval/W=$win xc-hw,y0-hw,xc+hw,y0+hw
	y0 -= nhw*hw
	SetDrawEnv/W=$win xcoord= abs,ycoord= abs,fillfgc= (0,0,0)
	DrawOval/W=$win xc-hw,y0-hw,xc+hw,y0+hw
	return y0-nhw/2*hw			// returns next bottom level
End

//	******************************** end of drawing boxesy *********************************
//	*********************************************************************************