I am currently working on a sub-problem of a larger problem which involves being able to efficiently align an equatorial telescope mount with greater precision than the conventional amateur method, with the added advantage of being able to do so without needing to see the polar point of the celestial equator (basically trying to implement this from scratch for myself for precision astrophotography purposes$^\dagger$). While I have fleshed out the theoretical framework of my solution, it needs to actually be implemented in code, and not merely handwaved through.
First, some background:
The standard system astronomers use to assign coordinates to points on the celestial sphere essentially projects the latitude and longitude lines on Earth into the sky, on so-called lines of Declination (Dec) and Right Ascension (RA), respectively. Zero degrees of Declination corresponds to the celestial equator, whereas zero Right Ascension corresponds to an "arbitrary" great circle similar to the arbitrariness of the prime meridian (no, they are not the same with one projected to the other). What we get is essentially a spherical coordinate system with all but two points represented uniquely by a pair of angles $(\theta, \phi)$; note, however, that the pole points of Earth's rotational axis don't have unique Right Ascension coordinates due to gimbal lock.
On the other hand, it is occasionally useful to express the position of a point on the celestial sphere in terms of its altitude angle above one's local horizon (Alt) and its azimuth angle (Az) with respect to true north, again from one's own perspective. This is, as before, a spherical coordinate system with a pair of angles $(\alpha, \omega)$, but this time the angles are measured with respect to different axes. Note that here, lockout happens at the zenith / nadir, which does not have a unique azimuth coordinate.
Question:
Suppose I have several points on the celestial sphere whose position I have expressed in both (RA, Dec) coordinates as well as (Alt, Az) coordinates (the latter relative to a specific time & location). How can I use this information to compute a conversion from the former coordinate system to the latter when those systems differ only in where the point $(0,0)$ is said to be located?
My intuition had me give change-of-basis a stab, or maybe I could find a transformation matrix $A$ whose values I can find by chugging out $A\mathbf{x}_1 = \mathbf{y}_1$ and $A\mathbf{x}_2 = \mathbf{y}_2$ explicitly & solving for unknowns, but this doesn't feel right because the "origins" (where both reference angles are $0$) change, whereas a linear transformation would preserve the origin. So it feels like a spherical version of an affine-but-not-linear transformation, and I'm a bit stumped about how to deal with it. In $\mathbb{R}^2$, when we translate the origin, say to $(1,0)$, then we may translate all other points $(x,y)$ by simply adding $1$ to the $x$-coordinate. However, I am suspicious that this "just rotate & translate everything by the same amount" will work on a sphere without some sort of scaling factors.
For what it's worth, I'm also able to express both gimbal-lockout points and both origins of either coordinate system in terms of the other.
Thanks for reading & for any feedback!
$^\dagger$ The idea is to first get my mount's RA axis roughly aligned. Next, I'll pan the telescope along this mount-axis only (keeping the other mount-axis locked), taking $3$ different pictures along the way. Plate-solving these pictures yields $3$ points whose (RA, Dec) coordinates are known exactly (if the RA axis of the mount is perfectly aligned with the rotational axis of the Earth, these points will have equal Declination; otherwise, Declination will drift). These $3$ points uniquely determine a plane; linear algebra can be used to determine the pole point of the mount's alignment by finding the vector orthogonal to that plane whose tail is at the center of the Earth. Annoyingly, adjustment of the mount's alignment is carried out exclusively along the Altitude-Azimuth axes. As such, I need to be able to convert the (RA, Dec) coordinates of my mount's pole point into (Alt, Az) so that I know how much correction needs to be applied in those directions to get the mount's computed pole point to agree with that of Earth's rotation (the latter of which has Alt-Az coordinates $(\text{Lat}, 0)$, where $\text{Lat}$ is your current latitude on Earth).
Best Answer
Let the standard astronomical system be
$$ X(\theta,\phi, R) = R \left( \cos \theta \cos \phi, \ \cos \theta \ \sin \phi, \sin \theta\right)$$
with $\theta=0$ equatorial plane, $\theta=\pm \pi/2$ the polar directions and $\phi=0$ the Greenwich meridian.
In order to find the map to a second system with polar direction $\{\theta_0,\phi_0\}$ you have to construct the rotation matrix in this plane spanned by the two polar vectors with an angle of rotation parameter $\theta_0$.
So the single task is to find the generator of such a rotation. The usual method is to take the cross product and norm it to unity
$$ \vec {n'} = (0,0,1) \times X(\theta_0,\phi_0,1)$$ But its much easier, the rotation plane is the plane of constant $\phi_0$ and the unit vector of rotation is simply the unit vector in the equatorial plane at $\phi=\phi_0\pm\pi/2$
Then the group of rotations in this plane is
$$\text{MatrixExp}\left[t \left(\left( \begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \\ \end{array} \right) \sin \left(\phi _0\right)+\left( \begin{array}{ccc} 0 & 0 & 1 \\ 0 & 0 & 0 \\ -1 & 0 & 0 \\ \end{array} \right) \cos \left(\phi _0\right)\right)\right] = \left( \begin{array}{ccc} \cos ^2\left(\frac{t}{2}\right)-\sin ^2\left(\frac{t}{2}\right) \cos \left(2 \phi _0\right) & \sin ^2\left(\frac{t}{2}\right) \sin \left(2 \phi _0\right) & \sin (t) \cos \left(\phi _0\right) \\ \sin ^2\left(\frac{t}{2}\right) \sin \left(2 \phi _0\right) & \sin ^2\left(\frac{t}{2}\right) \cos \left(2 \phi _0\right)+\cos ^2\left(\frac{t}{2}\right) & -\sin (t) \sin \left(\phi _0\right) \\ \sin (t) \left(-\cos \left(\phi _0\right)\right) & \sin (t) \sin \left(\phi _0\right) & \cos (t) \\ \end{array} \right)$$
Setting $t=\pi/2 - \theta_0$ for t=0, if th polar axes coincide, one finally has the compact formula
$$X' = \left( \begin{array}{c} \cos \left(\phi _0\right) \left(\sin \left(\theta _0\right) \cos (\theta ) \cos \left(\phi +\phi _0\right)+\sin (\theta ) \cos \left(\theta _0\right)\right)+\cos (\theta ) \sin \left(\phi _0\right) \sin \left(\phi +\phi _0\right) \\ \cos (\theta ) \sin \left(\phi +\phi _0\right) \cos \left(\phi _0\right)-\sin \left(\phi _0\right) \left(\sin \left(\theta _0\right) \cos (\theta ) \cos \left(\phi +\phi _0\right)+\sin (\theta ) \cos \left(\theta _0\right)\right) \\ \sin (\theta ) \sin \left(\theta _0\right)-\cos (\theta ) \cos \left(\theta _0\right) \cos \left(\phi +\phi _0\right) \\ \end{array} \right)$$
Disclaimer: This is a result produced by Mathematica GPT, so better control all of the results. For Mathematica users