Angle
subtended at the eye by neighboring stars W. Kahan
September 17, 2004 1:01 pm
Spherical polar coordinates {R, H, A} convert to Cartesian (X, Y, Z) :
R = radius, H = angle of Elevation above horizon, A = angle of Azimuth from compass North.
R > 0 , -
π
/2
≤
H
≤
π
/2 , -
π
≤
A
≤
π
.
X = distance North, Y = distance West, Z = distance Up ;
X = R · cos H · cos A , Y = R · cos H · sin A , Z = R · sin H .
Neighboring stars: Put one at (X, Y, Z) ‹—› { R , H , A } for any arbitrary R > 0 ;
X = R·cos H · cos A , Y = R·cos H · sin A , Z = R·sin H ;
another at (X+x, Y+y, Z+z) ‹—› { R , H+h , A+a } ;
X+x = R·cos(H+h)·cos(A+a) , Y+y = R·cos(H+h)·sin(A+a) , Z+z = R·sin(H+h) .
The angle v subtended at the eye by the stars satisfies (2·R·sin v/2)
2
= x
2
+ y
2
+ z
2
, so
(2·sin v/2)
2
= (cos(H+h)·cos(A+a) – cos H · cos A )
2
+ ( cos(H+h)·sin(A+a) – cos H · sin A )
2
+ ( sin(H+h) – sin H )
2
= 2 – 2·sin(H+h)·sin H - 2·cos(H+h)·cos(H)·cos(a)
= 2·(1 – cos(h)) + 2·(1 – cos a)·cos(H+h)·cos(H)
= 4·sin
2
(h/2) + 4·sin
2
(a/2)·cos(H+h)·cos(H) .
Conclusion: Of the two azimuths A and A+a only their difference a mod 2
π
matters; then
v = arccos( sin(H+h)·sin(H) + cos(H+h)·cos(H)·cos(a) )
= 2·arcsin
√
(
sin
2
(h/2) + sin
2
(a/2)·cos(H+h)·cos(H)
)
.
The first formula malfunctions at small subtended angles. The second is numerically fine at small
angles v but loses almost half the precision carried if v
≈
π
though no cancellation occurs.
Example: a = 179
.
999˚ , h = 52˚ , H = –26˚ ; carrying 10 sig. dec., this v = 180˚ instead of 179
.
999101˚ .
We can do better. Since -
π
/2
≤
H
≤
π
/2 and -
π
/2
≤
H+h
≤
π
/2 ,
0
≤
cos(H+h)·cos(H) = ( cos(2H+h) + cos(h) )/2 = cos
2
(H+h/2) – sin
2
(h/2)
= cos
2
(h/2) - sin
2
(H+h/2) .
Therefore
tan
2
(v/2) = sin
2
(v/2)/(1 – sin
2
(v/2))
=
(
sin
2
(h/2) + sin
2
(a/2)·cos(H+h)·cos(H)
)/(
cos
2
(h/2) – sin
2
(a/2)·cos(H+h)·cos(H)
)
=
(
sin
2
(h/2)·cos
2
(a/2) + sin
2
(a/2)·cos
2
(H+h/2)
)/(
cos
2
(h/2)·cos
2
(a/2) + sin
2
(a/2)·sin
2
(H+h/2)
) .
Hence follows a numerically accurate and efficient formula:
v = 2·arctan
√
((
th·(1 + ta + TH) + ta
)/(
1 + TH·(1 + ta·(1+th))
))
wherein ta = tan
2
(a/2) , th = tan
2
(h/2) and TH = tan
2
(H+h/2) . Only if h = 0 and H =
±π
/2
does TH =
∞
. Therefore presubstitute 0·
∞
= 0 and
∞
/
∞
= 1/TH to handle
all
special cases.
( Actually, infinite floating-point values of tan(…) can’t arise unless angles are in degrees.)
The same formula works with astronomical Declination instead of Elevation and Right
Ascension instead of Azimuth. A similar formula would work for distance over the surface of
the Earth, using Latitude instead of Elevation and Longitude instead of Azimuth, if the
Earth were perfectly spherical.
This document was created with FrameMaker404