pro genpeak_positions,abc,lambda,ttheta,unique=unique
;
a=abc.a
b=abc.b
c=abc.c
;
; Convert angles to radians
alpha=!dtor*abc.al
beta=!dtor*abc.be
gamma=!dtor*abc.ga
tthetamin=!dtor*abc.tthmin
tthetamax=!dtor*abc.tthmax
;
; Initialization
ttheta=fltarr(4)
;;res=dialog_message(["Is there a better way to generate maxhkl?","e.g. max(a,b,c)/dmin"],/question)
maxhkl=15
;
; A little bit of processing to get the metric tensor and the peak positions
;
g=[[a*a,a*b*cos(gamma),a*c*cos(beta)],[a*b*cos(gamma),b*b,b*c*cos(alpha)],[c*a*cos(beta),c*b*cos(alpha),c*c]]
inv_g=invert(g,/double)
;
dmin=lambda/(2.0*sin(tthetamax/2.0))
dmax=lambda/(2.0*sin(tthetamin/2.0))
;
for h=-maxhkl,maxhkl do begin
	for k=-maxhkl,maxhkl do begin
		for l=-maxhkl,maxhkl do begin
			dstar=sqrt([h,k,l]#inv_g#[h,k,l])
			if(dstar gt 0.0) then begin
				d=1.0/dstar
				if (d gt dmin and d lt dmax) then begin
					ttheta=[[ttheta],[h,k,l,2.0*asin(lambda/(2.0*d))*!radeg]]; "visible ttheta"
				endif
			endif
		endfor
	endfor
endfor
;
; Unsorted peak positions
nrefl=n_elements(ttheta)/4-1
if (nrefl gt 0) then begin
	ttheta=ttheta[*,1:nrefl]
;
; Sorted peak positions
	ttheta=ttheta[*,sort(ttheta[3,*])]
	ttheta[3,*]=round(1000.0*ttheta[3,*])*0.001
	if (unique) then ttheta=ttheta[*,uniq(ttheta[3,*])]
endif else begin
	ttheta=0.0
endelse
;
end
