Computes rigorous great circle arc distances.
Input position can either be either radians, sexagesimal RA, Dec or
degrees. All computations are double precision.
GCIRC, U, RA1, DC1, RA2, DC2, DIS
U -- integer = 0,1, or 2: Describes units of inputs and output:
0: everything radians
1: RAx in decimal hours, DCx in decimal
degrees, DIS in arc seconds
2: RAx and DCx in degrees, DIS in arc seconds
RA1 -- Right ascension or longitude of point 1
DC1 -- Declination or latitude of point 1
RA2 -- Right ascension or longitude of point 2
DC2 -- Declination or latitude of point 2
DIS -- Angular distance on the sky between points 1 and 2
See U above for units; double precision
"Haversine formula" see
(1) If RA1,DC1 are scalars, and RA2,DC2 are vectors, then DIS is a
vector giving the distance of each element of RA2,DC2 to RA1,DC1.
Similarly, if RA1,DC1 are vectors, and RA2, DC2 are scalars, then DIS
is a vector giving the distance of each element of RA1, DC1 to
RA2, DC2. If both RA1,DC1 and RA2,DC2 are vectors then DIS is a
vector giving the distance of each element of RA1,DC1 to the
corresponding element of RA2,DC2. If the input vectors are not the
same length, then excess elements of the longer ones will be ignored.
(2) The function SPHDIST provides an alternate method of computing
a spherical distance.
(3) The haversine formula can give rounding errors for antipodal
Written in Fortran by R. Hill -- SASC Technologies -- January 3, 1986
Translated from FORTRAN to IDL, RSH, STX, 2/6/87
Vector arguments allowed W. Landsman April 1989
Prints result if last argument not given. RSH, RSTX, 3 Apr. 1998
Remove ISARRAY(), V5.1 version W. Landsman August 2000
Added option U=2 W. Landsman October 2006
Use double precision for U=0 as advertised R. McMahon/W.L. April 2007
Use havesine formula, which has less roundoff error in the
milliarcsecond regime W.L. Mar 2009