Solution (source code)

= Solution

Let $\mathbf n,\mathbf w,\mathbf u$ be orthonormal <unit vectors> north, west and upward. A direction of <altitude> $\gamma$ and westward <azimuth> $\eta$ is
$$
\mathbf x=\cos\gamma\cos\eta\,\mathbf n+\cos\gamma\sin\eta\,\mathbf w+\sin\gamma\,\mathbf u.
$$
The north <celestial pole> and the upper equatorial meridian direction are
$$
\mathbf p=\cos\theta\,\mathbf n+\sin\theta\,\mathbf u,\qquad
\mathbf e=-\sin\theta\,\mathbf n+\cos\theta\,\mathbf u.
$$
In the <equatorial coordinate system>, the same <unit vector> is $\mathbf x=\sin\delta\,\mathbf p+\cos\delta(\cos h\,\mathbf e+\sin h\,\mathbf w)$. Taking <dot products> with $\mathbf p,\mathbf w,\mathbf e$ respectively gives
$$
\boxed{\begin{aligned}
\sin\delta&=\sin\gamma\sin\theta+\cos\gamma\cos\theta\cos\eta,\\
\cos\delta\sin h&=\sin\eta\cos\gamma,\\
\cos\delta\cos h&=\sin\gamma\cos\theta-\cos\gamma\sin\theta\cos\eta.
\end{aligned}}
$$
This is an <orthogonal transformation> between two <bases>, so it preserves the length of the direction vector. Recover $h$ from both its <sine> and <cosine> to keep the correct quadrant. At the <zenith>, <azimuth> itself is undefined; the vector equations remain meaningful by continuity. Using the more usual eastward <azimuth> changes the sign of the second equation.