;-------------------------------------------------------------
;+
; 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