#pragma rtGlobals=1		// Use modern global access method.
#pragma ModuleName=StereoGraphicProjection
#pragma version = 2.4
#include "LatticeSym", version>=3.1


Menu "Analysis"
	"Stereographic Projection", MakeStereo(NaN,NaN,NaN,NaN,NaN)
End


Function MakeStereo(Hz,Kz,Lz,hklmax,fntSize,[Qmax])
	Variable Hz,Kz,Lz			// hkl of the pole
	Variable hklmax
	Variable fntSize
	Variable Qmax				// maximum Q (1/nm)

	Variable Qpassed = !ParamIsDefault(Qmax)		// a Qmax was passed
	Qmax = Qpassed ? Qmax : Inf						// if not passed, do not limit the Q (it is limited by hklmax)
	if (!(fntSize>0) || !(hklmax>0) || numtype(Hz+Kz+Lz))
		Hz = numtype(Hz) ? 3 : Hz
		Kz = numtype(Kz) ? -1 : Kz
		Lz = numtype(Lz) ? 7 : Lz
		fntSize = fntSize>0 ? fntSize : 9
		hklmax = hklmax>0 ? hklmax : 6
		String hklStr
		sprintf hklStr,"%g %g %g",Hz,Kz,Lz
		Prompt hklStr,"(hkl) of the pole"
		Prompt fntSize,"size of font for hkl"
		Prompt hklmax,"maximum hkl to accept"
		Prompt Qmax,"maximum Q to show (a redundant limit) (1/nm)"
		DoPrompt "hkl max & font size",hklStr,hklmax,Qmax,fntSize
		if (V_flag)
			return 1
		endif
		sscanf hklStr, "%g %g %g",Hz,Kz,Lz
		Qmax = numtype(Qmax) ? Inf : Qmax
		if (numtype(Qmax)!=1)
			printf "MakeStereo(%g,%g,%g,  %g,%g, Qmax=%g)\r",Hz,Kz,Lz,hklmax,fntSize,Qmax
		else
			printf "MakeStereo(%g,%g,%g,  %g,%g)\r",Hz,Kz,Lz,hklmax,fntSize
		endif
	endif
	if (!(fntSize>0) || !(hklmax>0) || numtype(Hz+Kz+Lz))
		return 1
	endif
	if (setLattice())
		return 1
	endif

	STRUCT crystalStructure xtal
	if (FillCrystalStructDefault(xtal))				//fill the lattice structure with current values
		DoAlert 0, "No Lattice, please set one"
		return 1
	endif

	Variable maxWaveLen=9000
	Make/N=(maxWaveLen)/O longs_MakeStereo,lats_MakeStereo
	Make/N=(maxWaveLen,3)/W/O hkls_MakeStereo
	Wave hkls=hkls_MakeStereo
	Wave longs=longs_MakeStereo, lats=lats_MakeStereo
	longs=NaN
	lats=NaN
	Make/N=3/O/D qhat_MakeStereo
	Wave qhat = qhat_MakeStereo
	Variable qmag										// magnitude of qvector

	Variable hmax, kmax, lmax						// get separate limits on h,k, & l based on hklmax and Qmax
	qmag = sqrt(xtal.as0^2 + xtal.as1^2 + xtal.as2^2)
	hmax = min(ceil(Qmax/qmag),hklmax)
	qmag = sqrt(xtal.bs0^2 + xtal.bs1^2 + xtal.bs2^2)
	kmax = min(ceil(Qmax/qmag),hklmax)
	qmag = sqrt(xtal.cs0^2 + xtal.cs1^2 + xtal.cs2^2)
	lmax = min(ceil(Qmax/qmag),hklmax)
	Variable H,K,L
	Variable N=0, lat, long
	String str
	// Variable timer = startMSTimer
	for (L=0;L<=lmax;L=signedAdvance(L))
		for (K=0;K<=kmax;K=signedAdvance(K))
			for (H=0;H<=hmax;H=signedAdvance(H))
				if (H==0 && K==0 && L==0)		// skip the (000)
					continue
				endif
				if (magsqr(Fstruct(xtal,h,k,l)) < 0.001)
					continue
				endif
				qhat[0] = H*xtal.as0+K*xtal.bs0+L*xtal.cs0		// qmag*qhat =  recip x hkl
				qhat[1] = H*xtal.as1+K*xtal.bs1+L*xtal.cs1
				qhat[2] = H*xtal.as2+K*xtal.bs2+L*xtal.cs2
				qmag = normalize(qhat)
				if (qmag > Qmax)
					continue
				endif
				hkl2LongLat(Hz,Kz,Lz,qhat,long,lat,xtal)
				if (lat>=-0.01 && N<maxWaveLen)
					longs[N] = long
					lats[N] = lat
					hkls[N][0]=H  ;  hkls[N][1]=K  ;  hkls[N][2]=L
					N += 1
				endif
			endfor
		endfor
	endfor
	// printf "loop execution took %.2f sec\r",stopMSTimer(timer)*1e-6
	Redimension/N=(N) longs,lats
	lats = abs(lats[p]) < 1e-12 ? 0 : lats[p]		// This is needed because of a bug in Project
	Project/M=1 longs, lats
	Wave W_XProjection=W_XProjection, W_YProjection=W_YProjection

	// remove duplicates
	Variable i, j
	for (i=N-1;i>=0;i-=1)
		for (j=0;j<i;j+=1)
			if (abs(W_XProjection[i]-W_XProjection[j])+abs(W_YProjection[i]-W_YProjection[j]) < 0.01)
				W_XProjection[i] = NaN
				W_YProjection[i] = NaN
				break
			endif
		endfor
	endfor

	Display/W=(358.5,41,848.25,518.75)/K=1	// draw the graph
	NewFreeAxis/O/L emptyLeft
	NewFreeAxis/O/B emptyBottom
	ModifyGraph gfMult=60,width={Aspect,1},tick=3,noLabel=2,standoff=0,axThick=0,freePos=0,freePos=0
	SetAxis emptyLeft -2,2
	SetAxis emptyBottom -2,2
//	sprintf str, "\Zr142#%d   %s\r\\Zr141(%d %d %d)\r90\\F'Symbol'°\\F]0",xtal.SpaceGroup,getSymString(xtal.SpaceGroup),Hz,Kz,Lz
//	sprintf str, "\Zr142#%d   %s\r\\Zr141(%s)\r90\\F'Symbol'°\\F]0",xtal.SpaceGroup,getSymString(xtal.SpaceGroup),hkl2str(Hz,Kz,Lz)
	if (strlen(xtal.desc))
		sprintf str, "\\JR\\Zr200%s\r\\Zr071#%d   %s\r\\Zr141(%s)\r90\\F'Symbol'°\\F]0",xtal.desc,xtal.SpaceGroup,getSymString(xtal.SpaceGroup),hkl2str(Hz,Kz,Lz)
	else
		sprintf str, "\\JR\\Zr142#%d   %s\r\\Zr141(%s)\r90\\F'Symbol'°\\F]0",xtal.SpaceGroup,getSymString(xtal.SpaceGroup),hkl2str(Hz,Kz,Lz)
	endif
	TextBox/N=pole/F=0/A=RT/X=0/Y=0 str
	Button zoomIn,pos={5,5},size={30,20},proc=ZoomStereoButtonProc,title="+"
	Button zoomOut,pos={5,25},size={30,20},proc=ZoomStereoButtonProc,title="Ñ"

	String list = ReplaceNumberByKey("stereoRadius","",90,"=")
	list = ReplaceNumberByKey("Hz",list,Hz,"=")
	list = ReplaceNumberByKey("Kz",list,Kz,"=")
	list = ReplaceNumberByKey("Lz",list,Lz,"=")
	SetWindow kwTopWin, userdata(stereoRadius)=list

	SetDrawLayer/K UserFront
	if (1)
		SetDrawEnv xcoord= emptyBottom,ycoord= emptyLeft,fillpat= 0
		SetDrawEnv linefgc= (43520,43520,43520)
		DrawOval -0.01,-0.01,0.01,0.01
	else
		SetDrawEnv xcoord= emptyBottom,ycoord= prel
		DrawLine 0,0.49,0,0.51						// draw a small cross at the center
		SetDrawEnv xcoord= prel,ycoord= emptyLeft
		DrawLine 0.49,0,0.51,0
	endif

	SetDrawEnv fillpat= 0
	DrawOval 0,0,1,1								// draw big circle
	Variable low
	for (i=0;i<N;i+=1)								// draw all of the hkl on the graph
		if (numtype(W_XProjection[i]+W_YProjection[i]))
			continue
		endif
		hklstr = hkl2str(hkls[i][0],hkls[i][1],hkls[i][2])
		hklstr = ReplaceString(",",hklstr,"")
		low = lowOrder(hkls[i][0],hkls[i][1],hkls[i][2])	// flag, a low order reflection, show in larger font & bold
		SetDrawEnv xcoord=emptyBottom,ycoord=emptyLeft,textxjust=1,textyjust=1,fstyle=low,fsize=(fntsize+4*low)
		DrawText W_XProjection[i],W_YProjection[i],hklstr
	endfor

	KillWaves/Z W_YProjection,W_XProjection, W_cross
	KillWaves/Z xhat_hkl2LongLat,yhat_hkl2LongLat,zhat_hkl2LongLat
	KillWaves/Z lats_MakeStereo, longs_MakeStereo, hkls_MakeStereo, qhat_MakeStereo
	KillWaves/Z FindPerpVector_hat, W_cross
End
//
Function signedAdvance(i)
	Variable i
	i = round(i)
	return (i<0) ? abs(i) : -(i+1)
End
//
Function ZoomStereoButtonProc(ctrlName) : ButtonControl
	String ctrlName

	String list = GetUserData("","","stereoRadius")
	Variable Hz = NumberByKey("Hz",list,"=")
	Variable Kz = NumberByKey("Kz",list,"=")
	Variable Lz = NumberByKey("Lz",list,"=")
	Variable stereoRadius = NumberByKey("stereoRadius",list,"=")
	stereoRadius = numtype(stereoRadius) ? 90 : stereoRadius

	if (stringmatch(ctrlName,"zoomIn"))
		stereoRadius -= 5
	elseif (stringmatch(ctrlName,"zoomOut"))
		stereoRadius += 5
	else
		return 1
	endif
	stereoRadius = limit(stereoRadius,5,90)
	Variable d = 2*sin(stereoRadius*PI/180) / (1+cos(stereoRadius*PI/180))
	SetAxis emptyLeft -d,d
	SetAxis emptyBottom -d,d
	String str
	STRUCT crystalStructure xtal
	if (FillCrystalStructDefault(xtal))					//fill the lattice structure with current values
		sprintf str, "\\Zr200(%s)\r%g\\F'Symbol'°\\F]0",hkl2str(Hz,Kz,Lz),stereoRadius
	else
		if (strlen(xtal.desc))
			sprintf str, "\\JR\\Zr200%s\r\\Zr071#%d   %s\r\\Zr141(%s)\r%g\\F'Symbol'°\\F]0",xtal.desc,xtal.SpaceGroup,getSymString(xtal.SpaceGroup),hkl2str(Hz,Kz,Lz),stereoRadius
		else
			sprintf str, "\\JR\\Zr142#%d   %s\r\\Zr141(%s)\r%g\\F'Symbol'°\\F]0",xtal.SpaceGroup,getSymString(xtal.SpaceGroup),hkl2str(Hz,Kz,Lz),stereoRadius
		endif
	endif
	TextBox/C/N=pole str
	list = ReplaceNumberByKey("stereoRadius",list,stereoRadius,"=")
	SetWindow kwTopWin, userdata(stereoRadius)=list
	return 0
End
//
Static Function hkl2LongLat(Hz,Kz,Lz,qhat,long,lat,xtal)
	Wave qhat					// unit vector in direction of qvector
	Variable &long, &lat
	Variable Hz,Kz,Lz			// hkl of the pole
	STRUCT crystalStructure &xtal

	Make/N=3/O/D xhat_hkl2LongLat,yhat_hkl2LongLat,zhat_hkl2LongLat
	Wave xhat=xhat_hkl2LongLat, yhat=yhat_hkl2LongLat, zhat=zhat_hkl2LongLat
	zhat[0] = Hz*xtal.as0+Kz*xtal.bs0+Lz*xtal.cs0	// zhat =  recip x {Hz,Kz,Lz}
	zhat[1] = Hz*xtal.as1+Kz*xtal.bs1+Lz*xtal.cs1
	zhat[2] = Hz*xtal.as2+Kz*xtal.bs2+Lz*xtal.cs2
	FindPerpVector(zhat,xhat)

	Cross zhat, xhat
	Wave W_Cross=W_Cross
	yhat = W_Cross
	normalize(xhat)
	normalize(zhat)
	normalize(yhat)

	lat = 90 - acos(min(MatrixDot(zhat,qhat),1))*180/PI
	long = atan2(MatrixDot(yhat,qhat),MatrixDot(xhat,qhat))*180/PI
	return 0
End
//
Function FindPerpVector(ref,perp)	// find a vector perpendicular to ref
	Wave ref							// reference vector, find a vec perp to this
	Wave perp							// resulting perpendicular vector, not normalized
	Make/N=3/O/D FindPerpVector_hat
	Wave hat=FindPerpVector_hat

	hat = abs(ref)
	WaveStats/M=1/Q hat
	hat = 0
	hat[V_maxLoc] = 1
	Variable dot = MatrixDot(ref,hat)
	perp = ref - dot*hat
	if (norm(perp)<1e-9)				// ref is exacly along x, y, or z
		perp = 0
		perp[V_minloc]=1
	endif
	Cross ref, perp
	Wave W_cross=W_cross
	perp = W_cross
//	KillWaves/Z FindPerpVector_hat, W_cross
	return 0
End

Static Function lowOrder(h,k,l)		// returns true if (hkl) is a low order reflection
	Variable h,k,l
	h = abs(h)
	k = abs(k)
	l = abs(l)

	Variable zeros = (!h + !k + !l)
	if (zeros >= 2)
		return 1						// (100) type 
	elseif (zeros==1 && (h==k || h==l || k==l))
		return 1						// (110) type
	elseif (h==k && h==l)
		return 1						// (111) type
	endif
	return 0							// not a low order reflection
End


//Function ALLOW_All(h,k,l)				// allow all hkl, always returns 1
//	Variable h,k,l
//	return 1
//End
////
//Static Function ALLOW_FC(h,k,l)				// for face-centered, hkl must be all even or all odd
//	Variable h,k,l
//	return !mod(h+k,2) && !mod(k+l,2)
//End
////
//Static Function ALLOW_BC(h,k,l)				// for body-centered, !mod(round(h+k+l),2)
//	Variable h,k,l
//	return !mod(round(h+k+l),2)
//End
////
//Static Function ALLOW_CC(h,k,l)				// for C-centered, !mod(round(h+k),2)
//	Variable h,k,l
//	return !mod(round(h+k),2)
//End
////
//Static Function ALLOW_AC(h,k,l)				// for A-centered, !mod(round(k+l),2)
//	Variable h,k,l
//	return !mod(round(k+l),2)
//End
////
//Static Function ALLOW_RHOM_HEX(h,k,l)		// for rhombohedral hexagonal, allowed are -H+K+L=3n or H-K+L=3n
//	Variable h,k,l
//	return !mod(-h+k+l,3) || !mod(h-k+l,3)
//End
////
//Static Function ALLOW_HEXAGONAL(h,k,l)		// for hexagonal, // forbidden are: H+2K=3N with L odd
//	Variable h,k,l
//	return mod(h+2*k,3) || !mod(l,2)
//End



Function InitStereoGraphicPackage()
	InitLatticeSymPackage()							// used to initialize this package
End