;FUNCTION GCDIST
;  Purpose
;    Calculate the angular distance between two points on a sphere
;  Input
;    lat, lon  - One location. May be an array, but both lat and lon
;                need to be the same size
;    lat0,lon0 - The other location. May be an array, but if array,
;                lat and lon need to be either scalar or arrays of
;                the same size
;    /deg      - Set if input angles are in degrees, default is radians
;                Output is in same units
;  Return
;    Distance between points. If location is a scalar and the other is
;    an array, returns an array of same shape, with the distance from
;    the scalar point to each of the array points. If both are arrays,
;    returns an array of distance from one array point to corresponding
;    point in other array.
function gcdist,lat,lon,lat0,lon0,deg=deg
  if keyword_set(deg) then begin
    fac=180d/!dpi
  end else begin
    fac=1d
  end
  right=!dpi/2d
  a=right-(lat/fac);
  b=right-(lat0/fac);
  C=(lon-lon0)/fac;
  dist=acos(cos(a)*cos(b)+sin(a)*sin(b)*cos(C))
  return,dist*fac
end