;FUNCTION GCAZ
;  Purpose
;    Calculate the azimuth (compass bearing) from one point to another
;  Input
;    lat0,lon0  - Reference location. May be an array, but both lat and lon
;                need to be the same size
;    lat ,lon  - 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
;    Angle between great circle from reference point to other point, and due north
;    from the reference point, at the reference point.
;      North =   0deg =     0 rad
;      East  =  90deg =  pi/2 rad
;      South = 180deg =  pi   rad
;      West  = 270deg = 3pi/2 rad
;
function gcaz,_lat0,_lon0,_lat,_lon,deg=deg
  if keyword_set(deg) then begin
    lat=_lat/180d*!dpi
    lon=_lon/180d*!dpi
    lat0=_lat0/180d*!dpi
    lon0=_lon0/180d*!dpi
  end else begin
    lat=_lat
    lon=_lon
    lat0=_lat0
    lon0=_lon0
  end
  dlon=lon-lon0
  y=sin(dlon)*cos(lat)
  x=cos(lat0)*sin(lat)-sin(lat0)*cos(lat)*cos(dlon)
  az=atan(y,x)
  if keyword_set(deg) then az*=(180d/!dpi)
  return,az
end