function arclength,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
return,meanradius*arcangle(loc1,loc2,meanradius)
end