function arcangle,loc1,loc2,a

if n_params() eq 2 then a=0

if a eq 0 then begin
   loc1xyz=latlon2xyz(loc1,0,0)
   loc2xyz=latlon2xyz(loc2,0,0)
   meanradius=(sqrt(loc1xyz[0]^2+loc1xyz[1]^2+loc1xyz[2]^2)+sqrt(loc2xyz[0]^2+loc2xyz[1]^2+loc2xyz[2]^2))/2
endif else meanradius=a
loc1xyz=latlon2xyz(loc1,0,meanradius)
loc2xyz=latlon2xyz(loc2,0,meanradius)
chord=loc1xyz-loc2xyz
return,2*asin(sqrt(chord[0]^2+chord[1]^2+chord[2]^2)/(2*sqrt(loc1xyz[0]^2+loc1xyz[1]^2+loc1xyz[2]^2)))
end