function calpha,t1,t2,p1,p2

t1=double(t1)
t2=double(t2)
p1=double(p1)
p2=double(p2)

;R1 = (cos(p1)*sin(t1), sin(p1)*sin(t1), cos(t1))
;R2 = (cos(p2)*sin(t2), sin(p2)*sin(t2), cos(t2))
R12= sqrt(2*(1-sin(t1)*sin(t2)*cos(p1-p2) - cos(t1)*cos(t2)) )

cos_alpha = (sin(t1)+sin(t2))*sin((p2-p1)/2) / R12       

return,cos_alpha

end

n=100

t1=randomn(seed,n)*10
t2=randomn(seed,n)*10
p1=randomn(seed,n)*10
p2=randomn(seed,n)*10

ca=calpha(t1,t2,p1,p2)
print,'FIRST TEST'
print,'cos(alpha)=',ca
print,'alpha=',acos(ca)/!dtor

print,'----------------------------------------------------------------------'
pmm,ca

pause

print,'SECOND TEST'
;; second test--centroid 1 at (1,0,0), centroid 2 at random points
;         selected to have 0.8 <x < 0.9

n=1000

p1 = 0
t1 = !pi/2

p2=randomn(seed,n)
t2=randomn(seed,n)*6
x=cos(p2)*sin(t2)
w=where(x lt 0.9 and x gt 0.8)


p2=p2(w)
t2=t2(w)


ca=calpha(t1,t2,p1,p2)
a=acos(ca)
ww=sort(a)
a=a(ww)

print,'alpha (degrees)=',a/!dtor
pmm,a

end