;-------------------------------------------------------------
;+
; NAME:
;       SPHDIST
; PURPOSE:
;       Angular distance between points on a sphere.
; CALLING SEQUENCE:
;       d = sphdist(long1, lat1, long2, lat2)
; INPUTS:
;       long1 = longitude of point 1, scalar or vector
;       lat1 = latitude of point 1, scalar or vector
;       long2 = longitude of point 2, scalar or vector
;       lat2 = latitude of point 2, scalar or vector
;
; OPTIONAL KEYWORD INPUT PARAMETERS:
;       /DEGREES - means angles are in degrees, else radians.
; OUTPUTS:
;       d = angular distance between points (in radians unless /DEGREES
;           is set.)
; PROCEDURES CALLED:
;       RECPOL, POLREC
; NOTES:
;       (1) The procedure GCIRC is similar to SPHDIST(), but may be more 
;           suitable for astronomical applications.
;
;       (2) If long1,lat1 are scalars, and long2,lat2 are vectors, then
;           SPHDIST returns a vector giving the distance of each element of 
;           long2,lat2 to long1,lat1.   Similarly, if long1,lat1 are vectors,
;           and long2, lat2 are scalars, then SPHDIST returns a vector giving
;           giving the distance of each element of long1,lat1 to to long2,lat2. 
;           If both long1,lat1 and long2,lat2 are vectors then SPHDIST returns
;           vector giving the distance of each element of long1,lat1 to the 
;           corresponding element of long2, lat2.   If the input vectors are 
;           not of equal length, then excess elements of the longer ones will 
;           be ignored.
; MODIFICATION HISTORY:
;       R. Sterner, 5 Feb, 1991
;       R. Sterner, 26 Feb, 1991 --- Renamed from sphere_dist.pro
;
; Copyright (C) 1991, Johns Hopkins University/Applied Physics Laboratory
; This software may be used, copied, or redistributed as long as it is not
; sold and this copyright notice is reproduced on each copy made.  This
; routine is provided as is without any express or implied warranties
; whatsoever.  Other limitations apply as described in the file disclaimer.txt.
;	Converted to IDL V5.0   W. Landsman   September 1997
;-
;-------------------------------------------------------------
 
	function sphdist, long1, lat1, long2, lat2, $
	  help=hlp, degrees=degrees
 
	if (n_params(0) lt 4) or keyword_set(hlp) then begin
	  print,' Angular distance between points on a sphere.'
	  print,' d = sphdist(long1, lat1, long2, lat2)'
	  print,'   long1 = longitude of point 1.         in'
	  print,'   lat1 = latitude of point 1.           in'
	  print,'   long2 = longitude of point 2.         in'
	  print,'   lat2 = latitude of point 2.           in'
	  print,'   d = angular distance between points.  out'
	  print,' Keywords:'
	  print,'   /DEGREES means angles are in degrees, else radians.'
	  print,' Notes: points 1 and 2 may be arrays.'
	  return, -1
	endif
 
	cf = 1.0
	if keyword_set(degrees) then cf = !radeg
 
	;--- Convert both points to rectangular coordinates. ---
	polrec, 1.0, lat1/cf, rxy, z1
	polrec, rxy, long1/cf, x1, y1
	polrec, 1.0, lat2/cf, rxy, z2
	polrec, rxy, long2/cf, x2, y2
 
	;--- Compute vector dot product for both points. ---
	cs = x1*x2 + y1*y2 + z1*z2
 
	;--- Compute the vector cross product for both points. ---
	xc = y1*z2 - z1*y2
	yc = z1*x2 - x1*z2
	zc = x1*y2 - y1*x2
	sn = sqrt(xc*xc + yc*yc + zc*zc)
 
	;--- Convert to polar.  ------
	recpol, cs, sn, r, a
	return, cf*a
 
	end