pro great_circle, lat1,lon1,lat2,lon2,gcd,head 

;---------------------------------------------------------------------
; Purpose:
;   Computes the heading and distance between 2 locations.  
;
;   Some rules for the input locations,  lat1/lon1 and lat2/lon2:
;
;      both can be arrays of the same dimension
;      both can be scalars
;      one can be a scalar and one an array
;
; Input:
;   lat1....initial latitude, degrees N
;   lon1....initial longitude, degrees E
;   lat2....final latitude, degrees N
;   lon2....final longitude, degrees E
;
; Output:
;   head....heading, degrees, north =0, west =270
;   gcd.....great circle distance, km
;
; Source: Mark Hervig
;---------------------------------------------------------------------

;- some constants

   earthr = 6372.0    ; Earth's radius (km)
   pi     = 3.14159   ; yes
   d2r    = pi/180.0  ; converst degrees to radians

;- convert degrees to radians

   lata  = lat1 * d2r  
   lona  = lon1 * d2r
   latb  = lat2 * d2r
   lonb  = lon2 * d2r

;- check the input variable dimensions

   n1 = n_elements(lat1) & m1 = n_elements(lon1)
   n2 = n_elements(lat2) & m2 = n_elements(lon2)

   go = 0
   if (n1 eq n2)            then go = 1
   if (n1 eq 1 and n2 gt 1) then go = 1
   if (n2 eq 1 and n1 gt 1) then go = 1
   if (n1 ne m1)            then go = 0
   if (n2 ne m2)            then go = 0

;- do the calculations 

   if (go eq 1) then begin

     term  = sin(lata)*sin(latb) + cos(lata)*cos(latb)*cos(lonb-lona)

     angd  = acos(term)         ; dist in radians
     gcd   = angd * earthr      ; dist in km

     term  = (sin(latb)-cos(angd)*sin(lata)) / (sin(angd)*cos(lata)) 

     head  = acos(term)/d2r     ; heading in deg

     dlon  = lon1 - lon2
     dlat  = lat1 - lat2

;-   fix western headings

     west = where( (dlon lt 0.0 and abs(dlon) gt 180.0) or (dlon gt 0.0 and abs(dlon) lt 180.0),nk)
     if (nk gt 0) then head(west) = 360.0 - head(west) 

;-   fix due south headdings

     south = where(dlon eq 0.0 and dlat gt 0.0,nk) 
     if (nk gt 0) then head(south) = 180.0          

;-   fix due north headdings

     north = where(dlon eq 0.0 and dlat lt 0.0,nk) 
     if (nk gt 0) then head(north) = 0.0          

   endif

   if (go eq 0) then print,'lat1/lon1 and lat2/lon2 must have equall dimensions'
   if (go eq 0) then print,'OR one set can be a scalar and one set can be an array'
   if (go eq 0) then stop,'in great_circle.pro'

   return
   end