Angle from 2 point on celestial sphere

Question

I'd like to find the angle in radians between two equatorial coordinate i.e from RA dec of point 1 and RA and dec of point 2.

Answer 1

There are multiple formulae for how to calculate this angle. The simplest is to construct the unit vectors: $$\begin{align} \hat{n}_i & = \left[\begin{array}{c} \cos \delta_i \cos \alpha_i \\ \cos \delta_i \sin \alpha_i \\ \sin \delta_i \end{array}\right] \end{align}$$ where $\delta$ is the declination and $\alpha$ is the right ascension, then take the dot product of the vectors for the two positions, and finally use $\hat{n}_1 \cdot \hat{n}_2 = \cos \theta_{12}$. The dot product results in the following equation:

$$\cos(\theta_{12}) = \sin(\delta_1)\sin(\delta_2) + \cos(\delta_1)\cos(\delta_2)\cos(\alpha_1 - \alpha_2)$$

The problem with that approach is that it loses numerical accuracy very quickly when $\theta_{12} \ll 1$ radian. Thus the Haversine formula: $$\theta_{12} = 2\operatorname{arcsin} \left(\sqrt{\sin^2\left(\frac{\delta_2 - \delta_1}{2}\right) + \cos \delta_1 \cos \delta_2 \sin^2 \left(\frac{\alpha_2 - \alpha_1}{2}\right)}\right),$$ is often preferable. There is an algorithm that covers all cases, Vincenty's formulae, but it is much more complicated to implement.

You can find a good discussion of practical differences at this Geographic Information Systems Stack Exchange Question.

Answer 2

I wanted to make a comment on Sean's answer, but I do not have the reputation to do that. I guess I will take a stab at an answer with apologies to Sean.

When converting from spherical coordinates to cartesian, you have to pay attention to the assumptions of the spherical coordinate system being used. In this case the assumption used in the equatorial coordinate system about the meaning of the Declination angle is different then the assumption used in the first equation in Sean's answer. We know this because the value for x is (cos δi * sin αi). This is correct for inclination, or the case where the angle is between the vector and the Z axis. For equatorial coordinates the angle is between the vector and the XY plane.

Spherical coordinate systems using inclination are diagrammed here:

Wikipedia on spherical coordinate systems

The equatorial coordinate system is diagrammed here:

Wikipedia on the equatorial coordinate system

So to calculate the Cartesian vectors for a given equatorial point from the usual astronomical notation to degrees or radians. Right Ascension is given in hours, minutes and seconds. An hour is 15 degrees. (The are 24 hours in a complete revolution and 360/24 is 15.) A minute is 1/60 of an hour. 15/60 gives 1/4 so every minute is 1/4 of a degree. Finally, a second is 1/60 of a minute. Multiplying the fraction of a minute in a second (1/60) by the number of degrees in a minute (1/4) gives the number degrees in a second. (1/(60*4)) = 1/240. Multiply the hours by 15, then the minutes by 1/4 and finally the seconds by 1/240, then add the three numbers together.

More info here:

Wikipedia on Right Ascension

There is a similar bit of work for Declination. Declination is in Degrees, minutes, seconds. We don't have to convert from hours to degrees since the declination does not use the idea of hours. Minutes and seconds convert as with Right Ascension.

More here:

Wikipedia on Declination

After this conversion you can transform the coordinates from spherical to Cartesian. In order to do this, we need to know the radius, or distance from the origin of the coordinate system and the object. Equatorial points don't usually come with this information. We can use 1 because the distance doesn't change the angle between the two points. It also makes the math easier.

To get X, we start by calculating the projection of the vector on the XY plane.

LenXY = 1*cos(Declination)

We next need to calculate the projection of the vector on the XY plane onto the X axis.

X = LenXY*cos(Right Ascension) or

X = cos(Declination)*cos(Right Ascension)

The process for Y is similar but of course we want the projection on the Y axis. This will be sin(Right Ascension) instead of cos.

Y = LenXY*sin(Right Ascension) or

Y = cos(declination)*sin(Right Ascension)

Finally Z is the projection of the original vector on the Z axis. Since Declination is relative to the Z axis, although measured with respect to the XY plane, we can calculate the projection directly.

Z = 1*sin(Declination) or

Z = sin(Declination)

I believe the Haversine would work as it is meant for work with longitude and latitude. Latitude makes use of a declination angle the same way Declination in equatorial coordinates does.

Since the issue appears to be the limit of the amount of a precision typical floating point numbers in computers have, I would suggest looking at math libraries that support arbitrary precision. A library for Python is here: Link to Python Arbitrary Precision library

Finally, here is a link to a description on how to calculate the position of the moon and planets given their orbital elements.

Computing Planetary Positions