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

<p><span class="math-container">$$\cos(\theta_{12}) = \sin(\delta_1)\sin(\delta_2) + \cos(\delta_1)\cos(\delta_2)\cos(\alpha_1 - \alpha_2)$$</span></p>

<p>The problem with that approach is that it loses numerical accuracy very quickly when <span class="math-container">$\theta_{12} \ll 1$</span> radian. Thus the <A HREF="https://en.wikipedia.org/wiki/Haversine_formula" rel="nofollow noreferrer">Haversine formula</A>:
<span class="math-container">$$\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),$$</span>
is often preferable. There is an algorithm that covers all cases, <A HREF="https://en.wikipedia.org/wiki/Vincenty%27s_formulae" rel="nofollow noreferrer">Vincenty's formulae</A>, but it is much more complicated to implement.</p>

<p>You can find a good discussion of practical differences at this <A HREF="http://gis.stackexchange.com/questions/4906/why-is-law-of-cosines-more-preferable-than-haversine-when-calculating-distance-b">Geographic Information Systems Stack Exchange Question</A>.</p>

CC BY-SA 4.0
2019-02-22T15:17:21.273Z
Why would you use Vincenty's formulae for this? Those formulae are for computing geodesics on an ellipsoid by approximating elliptic integrals. But the celestial sphere is a perfect sphere, so no elliptic integrals are required.
PM 2Ring
CC BY-SA 4.0
2019-02-22T15:21:48.210Z
BTW, Vincenty's formulae were developed in 1975 for use on the desktop calculators of the time, modern computers can use superior algorithms. I strongly recommend Charles Karney's [geographiclib](https://geographiclib.sourceforge.io/). FWIW, Dr Karney is a major contributor to several Wikipedia articles on ellipsoids, especially [Geodesics on an ellipsoid](https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid).
PM 2Ring
CC BY-SA 4.0
2016-12-04T22:58:20.823Z
2019-02-24T23:35:18.983Z
2019-02-24T23:35:18.983Z
Sean Lake
<p>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.</p>

<p>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. </p>

<p>Spherical coordinate systems using inclination are diagrammed here:</p>

<p><a href="https://en.wikipedia.org/wiki/Spherical_coordinate_system" rel="nofollow noreferrer">Wikipedia on spherical coordinate systems</a></p>

<p>The equatorial coordinate system is diagrammed here:</p>

<p><a href="https://en.wikipedia.org/wiki/Equatorial_coordinate_system" rel="nofollow noreferrer">Wikipedia on the equatorial coordinate system</a></p>

<p>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.</p>

<p>More info here:</p>

<p><a href="https://en.wikipedia.org/wiki/Right_ascension" rel="nofollow noreferrer">Wikipedia on Right Ascension</a></p>

<p>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. </p>

<p>More here:</p>

<p><a href="https://en.wikipedia.org/wiki/Declination" rel="nofollow noreferrer">Wikipedia on Declination</a></p>

<p>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. </p>

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

<p>LenXY = 1*cos(Declination)</p>

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

<p>X = LenXY*cos(Right Ascension) or</p>

<p>X = cos(Declination)*cos(Right Ascension)</p>

<p>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.</p>

<p>Y = LenXY*sin(Right Ascension) or</p>

<p>Y = cos(declination)*sin(Right Ascension)</p>

<p>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. </p>

<p>Z = 1*sin(Declination) or</p>

<p>Z = sin(Declination)</p>

<p>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. </p>

<p>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:
<a href="http://mpmath.org/" rel="nofollow noreferrer">Link to Python Arbitrary Precision library</a></p>

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

<p><a href="https://stjarnhimlen.se/comp/tutorial.html" rel="nofollow noreferrer">Computing Planetary Positions</a></p>

CC BY-SA 4.0
2019-02-22T03:15:03.457Z
Hi @jeffry7. I believe Sean's answer is most correct except for an "obvious" mistake: the X and Y components cannot be identical. I will try to correct his post and add the equation that should be given (the equation for the dot product).
JohnHoltz
CC BY-SA 4.0
2019-02-24T03:17:15.220Z
Hi @JohnHoltz, that does look better. It might be helpful to add which angle in Sean's answer is which. If, like me, you don't know which symbol means which angle, then using his answer is not straight forward.
jeffry7
CC BY-SA 4.0
2019-02-22T02:28:17.460Z
2019-02-22T02:28:17.460Z
jeffry7
<p>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.</p>

<coordinate>
Angle from 2 point on celestial sphere
pacman
# Angle from 2 point on celestial sphere

(Asked by: username_0 on 2016-12-04)

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.

**Tags:** coordinate


---

## 2 answers

### Answer by username_1 on 2016-12-04. Score: 8

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](https://en.wikipedia.org/wiki/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](https://en.wikipedia.org/wiki/Vincenty%27s_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](http://gis.stackexchange.com/questions/4906/why-is-law-of-cosines-more-preferable-than-haversine-when-calculating-distance-b).

> **Comments:**
> **username_2**: Why would you use Vincenty's formulae for this? Those formulae are for computing geodesics on an ellipsoid by approximating elliptic integrals. But the celestial sphere is a perfect sphere, so no elliptic integrals are required.  
> **username_2**: BTW, Vincenty's formulae were developed in 1975 for use on the desktop calculators of the time, modern computers can use superior algorithms. I strongly recommend Charles Karney's [geographiclib](https://geographiclib.sourceforge.io/). FWIW, Dr Karney is a major contributor to several Wikipedia articles on ellipsoids, especially [Geodesics on an ellipsoid](https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid).  

---

### Answer by username_3 on 2019-02-22. Score: 2

I wanted to make a comment on username_1'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 username_1.

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 username_1'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](https://en.wikipedia.org/wiki/Spherical_coordinate_system)

The equatorial coordinate system is diagrammed here:

[Wikipedia on the equatorial coordinate system](https://en.wikipedia.org/wiki/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](https://en.wikipedia.org/wiki/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](https://en.wikipedia.org/wiki/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](http://mpmath.org/)

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](https://stjarnhimlen.se/comp/tutorial.html)

> **Comments:**
> **username_4**: Hi @username_3. I believe username_1's answer is most correct except for an "obvious" mistake: the X and Y components cannot be identical. I will try to correct his post and add the equation that should be given (the equation for the dot product).  
> **username_3**: Hi @username_4, that does look better. It might be helpful to add which angle in username_1's answer is which. If, like me, you don't know which symbol means which angle, then using his answer is not straight forward.  

---

