<p>The position angle P of a body ($\alpha_1, \delta_1$) with respect to another body ($\alpha_2, \delta_2$) can be calculated from
$$tan(P)={sin(\Delta\alpha)\over cos(\delta_2)tan(\delta_1)-sin(\delta_2)cos(\Delta\alpha)}$$
where $\Delta\alpha = \alpha_1-\alpha_2$. If the denominator is negative, the position angle lies in the range of 90 to 270 degrees.</p>

<p>Reference: Jean Meeus, Astronomical Algorithms, Second Edition, </p>

CC BY-SA 3.0
2018-03-01T05:17:38.613Z
IDL procedure [POSANG](https://idlastro.gsfc.nasa.gov/ftp/pro/astro/posang.pro) uses a similar formula, citing Smart.
Mike G
CC BY-SA 3.0
2018-03-04T00:45:59.770Z
Could anyone be bothered to check whether my prescription yield the same formula? I suspect it does.
ProfRob
CC BY-SA 3.0
2018-03-01T03:17:37.227Z
2018-03-01T03:17:37.227Z
JohnHoltz
<p>First, it should be noted that position angle is not just defined by "two positions". The starting point is distinct from the ending point, resulting in a difference of <span class="math-container">$180^\circ$</span>, depending on which one is first. </p>

<p>The previous answers were good, I just want to offer a different perspective/derivation. One way to define the position angle is that it's the angle from North counter-clockwise to the direction in question as measured on an <A HREF="https://en.wikipedia.org/wiki/Orthographic_projection_in_cartography" rel="nofollow noreferrer">orthographic projection</A> that has your starting point as the origin (assuming the astronomy standard of "North up, East left"). </p>

<p>The algebra behind the words begins by defining the unit vector/coordinates that define the starting point:
<span class="math-container">$$\hat{n}_0 = \left[\begin{array}{c}
  \cos\alpha_0 \cos\delta_0 \\
  \sin\alpha_0 \cos\delta_0 \\
  \sin\delta_0
\end{array}\right],$$</span>
with your ending position taking the same form with <span class="math-container">$0\rightarrow 1$</span>. Orthographic projection on any vector is simply the process of projecting off the component of the vector in the direction of <span class="math-container">$\hat{n}_0$</span>. The standard formula for this is
<span class="math-container">$$\vec{v}_\perp = \vec{v} - \hat{n}_0 (\vec{v}\cdot\hat{n}_0).\tag1$$</span></p>

<p>In principle, you could apply Equation (1) numerically with <span class="math-container">$\vec{v}$</span> as the North pole, and then as <span class="math-container">$\hat{n}_1$</span>, then the dot product between the normalized results will give you the position angle. </p>

<p>Doing it that way is probably a mistake, though. See, most positions in astronomy, <span class="math-container">$\hat{n}_0$</span> and <span class="math-container">$\hat{n}_1$</span>, will be separated by a small angle, so Equation (1) will suffer badly from <A HREF="https://en.wikipedia.org/wiki/Loss_of_Significance" rel="nofollow noreferrer">loss of significance</A>. That is why it's a good idea to continue the derivation to derive a formula that doesn't have this problem.</p>

<p>Examining the structure of <span class="math-container">$\hat{n}_0$</span> is worth doing because it will simplify the algebra tremendously. It is what you get if you start with the <span class="math-container">$x$</span>-direction unit vector, <span class="math-container">$\hat{x}$</span>, rotate around the <span class="math-container">$y$</span>-axis by <span class="math-container">$\delta_0$</span>, and then rotate around the <span class="math-container">$z$</span>-axis by <span class="math-container">$\alpha_0$</span>. Viewed this way, finding the two vectors that are perpendicular to <span class="math-container">$\hat{n}_0$</span> that we need is pretty straightforward - just apply the same rotation matrices to <span class="math-container">$\hat{y}$</span> and <span class="math-container">$\hat{z}$</span>. In other words, you can read them off from the columns of the rotation matrix
<span class="math-container">\begin{align}
 R &amp;= \left[\begin{array}{ccc} 
  \cos\alpha_0 &amp; -\sin\alpha_0 &amp; 0 \\
  \sin\alpha_0 &amp; \cos\alpha_0 &amp; 0 \\
  0 &amp; 0 &amp; 1 
\end{array}\right] \left[\begin{array}{ccc} 
  \cos\delta_0 &amp; 0 &amp; -\sin\delta_0 \\
  0 &amp; 1 &amp; 0 \\
  \sin\alpha_0 &amp; 0 &amp; \cos\alpha_0
\end{array}\right] \\
 &amp; = \left[\begin{array}{ccc} 
  \cos\alpha_0 \cos\delta_0 &amp; -\sin\alpha_0 &amp; -\cos\alpha_0 \sin\delta_0 \\
  \sin\alpha_0 \cos\delta_0 &amp; \cos\alpha_0 &amp; -\sin\alpha_0 \sin\delta_0 \\
  \sin\delta_0 &amp; 0 &amp; \cos\delta_0
\end{array}\right].\tag2
\end{align}</span>
I chose how to apply the signs to the sine functions in the two rotation matrices so that the first column would match <span class="math-container">$\hat{n}_0$</span>. </p>

<p>Let's call the second column of (2) <span class="math-container">$\hat{E}'$</span>, and the third column <span class="math-container">$\hat{N}'$</span>. Notice that if we rotate a standard set of <span class="math-container">$x$</span>-<span class="math-container">$y$</span> axes by 90 degrees counter-clockwise, then the <span class="math-container">$x$</span>-axis corresponds to North and the <span class="math-container">$y$</span> to East. Thus, we can use the standard formula for the 2-dimensional component of a vector and its polar angle if we identify <span class="math-container">$\hat{n}_1\cdot\hat{N}'$</span> as the <span class="math-container">$x$</span>-component and <span class="math-container">$\hat{n}_1\cdot\hat{E}'$</span> as the <span class="math-container">$y$</span>. That formula is
<span class="math-container">\begin{align} 
 P &amp;= \operatorname{atan2}\left(y,\,x\right)\\
   &amp; = \operatorname{atan2}\left(\sin\delta_1\cos\delta_0 - \sin\delta_0\cos\delta_1 \cos(\alpha_1-\alpha_0),\,\cos\delta_0\sin(\alpha_1-\alpha_0)\right). \tag3
\end{align}</span></p>

<p>Because <span class="math-container">$\cos\delta &gt; 0$</span> for all declinations, you could make the formula look more like the standard one from textbooks. My own instinct is to avoid using <span class="math-container">$\tan\delta$</span> because it diverges near the poles. This formula will work for all <span class="math-container">$\alpha$</span> and <span class="math-container">$\delta$</span>, as long as the two points are distinct. The only bit that remains is to fiddle with the <span class="math-container">$x$</span>-like argument to make it behave well, numerically, when the points are near each other. To do so, use <span class="math-container">$\sin\delta_1\cos\delta_0 = \sin(\delta_1-\delta_0) + \sin\delta_0\cos\delta_1$</span> and <span class="math-container">$1 - \cos(\alpha_1-\alpha_0) = 2\sin^2\left(\frac{\alpha_1 - \alpha_0}{2}\right)$</span> to get
<span class="math-container">$$P  = \operatorname{atan2}\left(\sin(\delta_1-\delta_0) + 2\sin\delta_0\cos\delta_1 \sin^2\left(\frac{\alpha_1-\alpha_0}{2}\right),\,\cos\delta_0\sin(\alpha_1-\alpha_0)\right). \tag4$$</span></p>

<p>In principle, you'd want to investigate when it's best to use (3) or (4) numerically. In practice, I suspect that (4) will do better than (3), in terms of numerical accuracy, in the vast majority of cases that astronomers care about.</p>

CC BY-SA 4.0
2019-08-01T10:09:24.720Z
2019-08-01T10:09:24.720Z
Sean Lake
<p>Assuming you mean the angle between the meridian line through A and the great circle that goes through points A and B, then it goes something like this.</p>

<p>Define vectors from the origin to A and B assuming they lie on a unit sphere, such that $x_A= \cos \delta_A \cos \phi_A$, $ y_A = \cos \delta_A \sin \phi_A$ and $z_A= \sin \delta_A$, and similar for B. Here $\phi$ refers to right Ascension and $\delta$ is declination.</p>

<p>The great circles in question define planes that pass through the origin. A normal to the plane defined by OAB is given by the vector product $\vec{n_1}=\vec{A}\times \vec{B}$. Similarly a great circle passing through the O, A and the NCP $(0, 0, 1)$ has a normal of $\vec{n_2}= \vec{A}\times (0,0,1) = (y_A, -x_A, 0)$.</p>

<p>The angle you are looking for <em>is</em> the angle between these two normal vectors, which can be found from scalar product in the usual way.
$$\cos \theta = \frac{\vec{n_1}\cdot \vec{n_2}}{|n_1||n_2|}=\frac{\cos \phi_{A}(y_{A}z_{B} - y_{B}z_A) + \sin \phi_{A}(x_Az_B - x_Bz_A)}{|\vec{A}\times \vec{B}|}$$</p>

CC BY-SA 3.0
2018-02-28T05:06:20.120Z
The second plane should pass through the NCP, the origin, and A.
Mike G
CC BY-SA 3.0
2018-02-28T08:33:30.207Z
@MikeG I think that's it.
ProfRob
CC BY-SA 3.0
2018-02-28T11:35:05.600Z
if this doesn't get an upvote I dunno what does!  :)
Fattie
CC BY-SA 3.0
2018-02-27T23:22:38.093Z
2018-02-28T08:33:27.067Z
2018-02-28T08:33:27.067Z
ProfRob
<p>Let's say I have the (RA,Dec) for two positions on the sky. I want to measure the position angle between these two positions. In other words, if I had the Cartesian coordinate vectors for these two positions (say [x,y,z]) then I think the position angle PA would simply be the dot product of the two vectors, so basically cos(PA). However, RA and Dec are spherical coordinates (corresponding to $\phi$ and $\theta$ respectively), not Cartesian coordinates. How do I get cos(PA) between the two positions on the sky if I only have their (RA,Dec) instead of [x,y,z]? </p>

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

<blockquote>
  <p>Position angle, usually abbreviated PA, is [a] convention for measuring angles on the sky in astronomy. The International Astronomical Union defines it as the angle measured relative to the north celestial pole (NCP), turning positive into the direction of the right ascension. In the standard (non-flipped) images this is a counterclockwise measure relative to the axis into the direction of positive declination.</p>
</blockquote>

CC BY-SA 3.0
2018-02-27T12:37:29.600Z
One way would be to use your dot product! Convert each position to a unit vector using $ \mathbf{n} = cos(\phi) cos(\theta) \mathbf{\hat{x}} \ + \ sin(\phi) cos(\theta) \mathbf{\hat{y}} \ + \ sin(\theta) \mathbf{\hat{z}}$ then the cosine of the angle between the two will be their inner product. And just fyi, Smart's [Spherical Astronomy](https://archive.org/details/SphericalAstronomy) is archived!
uhoh
CC BY-SA 3.0
2018-02-27T16:41:42.897Z
Astronomers conventionally use $\alpha$ for right ascension and $\delta$ for declination. Also you seem to be describing [angular separation](https://en.wikipedia.org/wiki/Angular_distance); [position angle](https://en.wikipedia.org/wiki/Position_angle) is something else.
Mike G
CC BY-SA 3.0
2018-02-27T18:17:31.390Z
http://www.gyes.eu/calculator/calculator_page1.htm gives the forumla !
Fattie
CC BY-SA 3.0
2018-02-27T18:21:53.263Z
@MikeG Sorry for the confusion, I am talking about position angle, not angular separation/distance. Position angle is the angle of a vector pointing from one (RA,Dec) to another (RA,Dec), usually measured East of North. Angular separation is much easier to calculate.
quantumflash
CC BY-SA 3.0
2018-02-27T18:22:14.347Z
@Fattie see my previous response, sorry for the confusion. I'm talking about position angle, not angular separation.
quantumflash
CC BY-SA 3.0
2018-02-27T18:25:49.327Z
@uhoh Thanks a lot, do you know if there is a more recent version/explanation of this? I can't find your exact formula in Chapter 2 of that book. Also, using something like astropy in Python, I can compute the Position Angle from one (RA,Dec) to another (RA,Dec) directly. After that, would cos(PA) be equivalent to the cos(angle) I get from your approach?
quantumflash
CC BY-SA 3.0
2018-02-27T18:31:43.857Z
@quantumflash I think it's best to wait for someone to post a proper answer with the optimal equations. If you have found them in Smart, you can post an answer to your own question, that's perfectly fine in Stack Exchange. You can also search this site to see if there are other answers with the equations. The problem with spherical trig is that some forms of th equations, while correct mathematically, break down when you are using them to calculate because of zeros in the denominator, so it's best to use a recommended formulation.
uhoh
CC BY-SA 3.0
2018-02-27T18:33:25.353Z
@RobJeffries Unfortunately I can't edit my comment to change "convert to" to "express as", but a normal vector can certainly *point* to a coordinate on the celestial sphere.
uhoh
CC BY-SA 3.0
2018-02-27T18:42:24.363Z
OP, I appreciate the link to wikipedia's ideas on what astronomers apparently (joke) mean by "Position Angle".  it is totally psychotic - astronomers are nuts  :)   Not that, apart from that issue as I believe R Jeffries pointed out, you last sentence does not parse - RA-Dec is completely unrelated to XYZ.  (Unless you mean "on a unit sphere"? or some such.)
Fattie
CC BY-SA 3.0
2018-02-27T18:44:52.913Z
I *think* I might know what you want: between two quaternions, q' = q-1q2 is a rotation rotate you from q1 to q2  .. is that it ???
Fattie
CC BY-SA 3.0
2018-02-27T18:47:29.820Z
Mind you, it's possible you want the **angle of the projection of the second one, away from the projection of North (?), onto the normal plane of the first one**.  Could that be it?!  (the "eyepiece" thing in the wiki article freaks me out!)
Fattie
CC BY-SA 3.0
2018-02-27T18:49:47.977Z
@RobJeffries I think there is a relation between [RA,Dec] and xyz -- I think the point of RA,Dec is to *project* the xyz onto the 2D celestial sky. I guess another way to put it is: can I just do dot product of the two (RA,Dec) vectors, and would that give me a position angle (East of North)?
quantumflash
CC BY-SA 3.0
2018-02-27T18:51:17.610Z
@Fattie yes exactly, (RA,Dec) are 2D coordinates on the celestial sky *after* you project 3D xyz onto the celestial sphere. In principle if you had something like redshift or distance, you could combine that with (RA,Dec) to get a 3D system where the redshift or distance gives you the radial component of the spherical coordinate system.
quantumflash
CC BY-SA 3.0
2018-02-27T18:53:53.167Z
@quantumflash, however the celestial sphere is *not a fixed distance away*.
Fattie
CC BY-SA 3.0
2018-02-27T18:54:29.707Z
thus, "I think there is a relation between [RA,Dec] and xyz -- I think the point of RA,Dec is to project the xyz onto the 2D celestial sky"  either that makes unfortunately absolutely *no* sense, or, there is some shared convention astronomers use, to which I am not privy.
Fattie
CC BY-SA 3.0
2018-02-27T19:23:28.850Z
Do you mean what is the angle between the meridian and a great circle that joins points A and B? Voting to close as unclear until this is tidied up.
ProfRob
CC BY-SA 3.0
2018-02-27T21:13:32.987Z
FWIW @RobJeffries, notice the comment by a folk under the answer - I'm thinking too it might be the angle between two great circles, for the three relevant points.
Fattie
<coordinate><declination><astrometry><right-ascension><alignment>
Measuring misalignment between two positions on sky
quantumflash
# Measuring misalignment between two positions on sky

(Asked by: username_0 on 2018-02-27)

Let's say I have the (RA,Dec) for two positions on the sky. I want to measure the position angle between these two positions. In other words, if I had the Cartesian coordinate vectors for these two positions (say [x,y,z]) then I think the position angle PA would simply be the dot product of the two vectors, so basically cos(PA). However, RA and Dec are spherical coordinates (corresponding to $\phi$ and $\theta$ respectively), not Cartesian coordinates. How do I get cos(PA) between the two positions on the sky if I only have their (RA,Dec) instead of [x,y,z]?

<https://en.wikipedia.org/wiki/Position_angle>

> Position angle, usually abbreviated PA, is [a] convention for measuring angles on the sky in astronomy. The International Astronomical Union defines it as the angle measured relative to the north celestial pole (NCP), turning positive into the direction of the right ascension. In the standard (non-flipped) images this is a counterclockwise measure relative to the axis into the direction of positive declination.

**Tags:** coordinate, declination, astrometry, right-ascension, alignment
> **Comments:**
> **username_1**: One way would be to use your dot product! Convert each position to a unit vector using $ \mathbf{n} = cos(\phi) cos(\theta) \mathbf{\hat{x}} \ + \ sin(\phi) cos(\theta) \mathbf{\hat{y}} \ + \ sin(\theta) \mathbf{\hat{z}}$ then the cosine of the angle between the two will be their inner product. And just fyi, Smart's [Spherical Astronomy](https://archive.org/details/SphericalAstronomy) is archived!  
> **username_2**: Astronomers conventionally use $\alpha$ for right ascension and $\delta$ for declination. Also you seem to be describing [angular separation](https://en.wikipedia.org/wiki/Angular_distance); [position angle](https://en.wikipedia.org/wiki/Position_angle) is something else.  
> **username_3**: http://www.gyes.eu/calculator/calculator_page1.htm gives the forumla !  
> **username_0**: @MikeG Sorry for the confusion, I am talking about position angle, not angular separation/distance. Position angle is the angle of a vector pointing from one (RA,Dec) to another (RA,Dec), usually measured East of North. Angular separation is much easier to calculate.  
> **username_0**: @Fattie see my previous response, sorry for the confusion. I'm talking about position angle, not angular separation.  
> **username_0**: @uhoh Thanks a lot, do you know if there is a more recent version/explanation of this? I can't find your exact formula in Chapter 2 of that book. Also, using something like astropy in Python, I can compute the Position Angle from one (RA,Dec) to another (RA,Dec) directly. After that, would cos(PA) be equivalent to the cos(angle) I get from your approach?  
> **username_1**: @quantumflash I think it's best to wait for someone to post a proper answer with the optimal equations. If you have found them in Smart, you can post an answer to your own question, that's perfectly fine in Stack Exchange. You can also search this site to see if there are other answers with the equations. The problem with spherical trig is that some forms of th equations, while correct mathematically, break down when you are using them to calculate because of zeros in the denominator, so it's best to use a recommended formulation.  
> **username_1**: @RobJeffries Unfortunately I can't edit my comment to change "convert to" to "express as", but a normal vector can certainly *point* to a coordinate on the celestial sphere.  
> **username_3**: OP, I appreciate the link to wikipedia's ideas on what astronomers apparently (joke) mean by "Position Angle".  it is totally psychotic - astronomers are nuts  :)   Not that, apart from that issue as I believe R Jeffries pointed out, you last sentence does not parse - RA-Dec is completely unrelated to XYZ.  (Unless you mean "on a unit sphere"? or some such.)  
> **username_3**: I *think* I might know what you want: between two quaternions, q' = q-1q2 is a rotation rotate you from q1 to q2  .. is that it ???  
> **username_3**: Mind you, it's possible you want the **angle of the projection of the second one, away from the projection of North (?), onto the normal plane of the first one**.  Could that be it?!  (the "eyepiece" thing in the wiki article freaks me out!)  
> **username_0**: @RobJeffries I think there is a relation between [RA,Dec] and xyz -- I think the point of RA,Dec is to *project* the xyz onto the 2D celestial sky. I guess another way to put it is: can I just do dot product of the two (RA,Dec) vectors, and would that give me a position angle (East of North)?  
> **username_0**: @Fattie yes exactly, (RA,Dec) are 2D coordinates on the celestial sky *after* you project 3D xyz onto the celestial sphere. In principle if you had something like redshift or distance, you could combine that with (RA,Dec) to get a 3D system where the redshift or distance gives you the radial component of the spherical coordinate system.  
> **username_3**: @quantumflash, however the celestial sphere is *not a fixed distance away*.  
> **username_3**: thus, "I think there is a relation between [RA,Dec] and xyz -- I think the point of RA,Dec is to project the xyz onto the 2D celestial sky"  either that makes unfortunately absolutely *no* sense, or, there is some shared convention astronomers use, to which I am not privy.  
> **username_4**: Do you mean what is the angle between the meridian and a great circle that joins points A and B? Voting to close as unclear until this is tidied up.  
> **username_3**: FWIW @RobJeffries, notice the comment by a folk under the answer - I'm thinking too it might be the angle between two great circles, for the three relevant points.  

---

## 3 answers

### Answer by username_4 on 2018-02-27. Score: 7 | Accepted answer

Assuming you mean the angle between the meridian line through A and the great circle that goes through points A and B, then it goes something like this.

Define vectors from the origin to A and B assuming they lie on a unit sphere, such that $x\_A= \cos \delta\_A \cos \phi\_A$, $ y\_A = \cos \delta\_A \sin \phi\_A$ and $z\_A= \sin \delta\_A$, and similar for B. Here $\phi$ refers to right Ascension and $\delta$ is declination.

The great circles in question define planes that pass through the origin. A normal to the plane defined by OAB is given by the vector product $\vec{n\_1}=\vec{A}\times \vec{B}$. Similarly a great circle passing through the O, A and the NCP $(0, 0, 1)$ has a normal of $\vec{n\_2}= \vec{A}\times (0,0,1) = (y\_A, -x\_A, 0)$.

The angle you are looking for *is* the angle between these two normal vectors, which can be found from scalar product in the usual way.
$$\cos \theta = \frac{\vec{n\_1}\cdot \vec{n\_2}}{|n\_1||n\_2|}=\frac{\cos \phi\_{A}(y\_{A}z\_{B} - y\_{B}z\_A) + \sin \phi\_{A}(x\_Az\_B - x\_Bz\_A)}{|\vec{A}\times \vec{B}|}$$

> **Comments:**
> **username_2**: The second plane should pass through the NCP, the origin, and A.  
> **username_4**: @MikeG I think that's it.  
> **username_3**: if this doesn't get an upvote I dunno what does!  :)  

---

### Answer by username_5 on 2018-03-01. Score: 4

The position angle P of a body ($\alpha\_1, \delta\_1$) with respect to another body ($\alpha\_2, \delta\_2$) can be calculated from
$$tan(P)={sin(\Delta\alpha)\over cos(\delta\_2)tan(\delta\_1)-sin(\delta\_2)cos(\Delta\alpha)}$$
where $\Delta\alpha = \alpha\_1-\alpha\_2$. If the denominator is negative, the position angle lies in the range of 90 to 270 degrees.

Reference: Jean Meeus, Astronomical Algorithms, Second Edition,

> **Comments:**
> **username_2**: IDL procedure [POSANG](https://idlastro.gsfc.nasa.gov/ftp/pro/astro/posang.pro) uses a similar formula, citing Smart.  
> **username_4**: Could anyone be bothered to check whether my prescription yield the same formula? I suspect it does.  

---

### Answer by username_6 on 2019-08-01. Score: 1

First, it should be noted that position angle is not just defined by "two positions". The starting point is distinct from the ending point, resulting in a difference of $180^\circ$, depending on which one is first.

The previous answers were good, I just want to offer a different perspective/derivation. One way to define the position angle is that it's the angle from North counter-clockwise to the direction in question as measured on an [orthographic projection](https://en.wikipedia.org/wiki/Orthographic_projection_in_cartography) that has your starting point as the origin (assuming the astronomy standard of "North up, East left").

The algebra behind the words begins by defining the unit vector/coordinates that define the starting point:
$$\hat{n}\_0 = \left[\begin{array}{c}
\cos\alpha\_0 \cos\delta\_0 \\
\sin\alpha\_0 \cos\delta\_0 \\
\sin\delta\_0
\end{array}\right],$$
with your ending position taking the same form with $0\rightarrow 1$. Orthographic projection on any vector is simply the process of projecting off the component of the vector in the direction of $\hat{n}\_0$. The standard formula for this is
$$\vec{v}\_\perp = \vec{v} - \hat{n}\_0 (\vec{v}\cdot\hat{n}\_0).\tag1$$

In principle, you could apply Equation (1) numerically with $\vec{v}$ as the North pole, and then as $\hat{n}\_1$, then the dot product between the normalized results will give you the position angle.

Doing it that way is probably a mistake, though. See, most positions in astronomy, $\hat{n}\_0$ and $\hat{n}\_1$, will be separated by a small angle, so Equation (1) will suffer badly from [loss of significance](https://en.wikipedia.org/wiki/Loss_of_Significance). That is why it's a good idea to continue the derivation to derive a formula that doesn't have this problem.

Examining the structure of $\hat{n}\_0$ is worth doing because it will simplify the algebra tremendously. It is what you get if you start with the $x$-direction unit vector, $\hat{x}$, rotate around the $y$-axis by $\delta\_0$, and then rotate around the $z$-axis by $\alpha\_0$. Viewed this way, finding the two vectors that are perpendicular to $\hat{n}\_0$ that we need is pretty straightforward - just apply the same rotation matrices to $\hat{y}$ and $\hat{z}$. In other words, you can read them off from the columns of the rotation matrix
\begin{align}
R &= \left[\begin{array}{ccc}
\cos\alpha\_0 & -\sin\alpha\_0 & 0 \\
\sin\alpha\_0 & \cos\alpha\_0 & 0 \\
0 & 0 & 1
\end{array}\right] \left[\begin{array}{ccc}
\cos\delta\_0 & 0 & -\sin\delta\_0 \\
0 & 1 & 0 \\
\sin\alpha\_0 & 0 & \cos\alpha\_0
\end{array}\right] \\
& = \left[\begin{array}{ccc}
\cos\alpha\_0 \cos\delta\_0 & -\sin\alpha\_0 & -\cos\alpha\_0 \sin\delta\_0 \\
\sin\alpha\_0 \cos\delta\_0 & \cos\alpha\_0 & -\sin\alpha\_0 \sin\delta\_0 \\
\sin\delta\_0 & 0 & \cos\delta\_0
\end{array}\right].\tag2
\end{align}
I chose how to apply the signs to the sine functions in the two rotation matrices so that the first column would match $\hat{n}\_0$.

Let's call the second column of (2) $\hat{E}'$, and the third column $\hat{N}'$. Notice that if we rotate a standard set of $x$-$y$ axes by 90 degrees counter-clockwise, then the $x$-axis corresponds to North and the $y$ to East. Thus, we can use the standard formula for the 2-dimensional component of a vector and its polar angle if we identify $\hat{n}\_1\cdot\hat{N}'$ as the $x$-component and $\hat{n}\_1\cdot\hat{E}'$ as the $y$. That formula is
\begin{align}
P &= \operatorname{atan2}\left(y,\,x\right)\\
& = \operatorname{atan2}\left(\sin\delta\_1\cos\delta\_0 - \sin\delta\_0\cos\delta\_1 \cos(\alpha\_1-\alpha\_0),\,\cos\delta\_0\sin(\alpha\_1-\alpha\_0)\right). \tag3
\end{align}

Because $\cos\delta > 0$ for all declinations, you could make the formula look more like the standard one from textbooks. My own instinct is to avoid using $\tan\delta$ because it diverges near the poles. This formula will work for all $\alpha$ and $\delta$, as long as the two points are distinct. The only bit that remains is to fiddle with the $x$-like argument to make it behave well, numerically, when the points are near each other. To do so, use $\sin\delta\_1\cos\delta\_0 = \sin(\delta\_1-\delta\_0) + \sin\delta\_0\cos\delta\_1$ and $1 - \cos(\alpha\_1-\alpha\_0) = 2\sin^2\left(\frac{\alpha\_1 - \alpha\_0}{2}\right)$ to get
$$P = \operatorname{atan2}\left(\sin(\delta\_1-\delta\_0) + 2\sin\delta\_0\cos\delta\_1 \sin^2\left(\frac{\alpha\_1-\alpha\_0}{2}\right),\,\cos\delta\_0\sin(\alpha\_1-\alpha\_0)\right). \tag4$$

In principle, you'd want to investigate when it's best to use (3) or (4) numerically. In practice, I suspect that (4) will do better than (3), in terms of numerical accuracy, in the vast majority of cases that astronomers care about.


---

