# Mathematical model
Here I describe the mathematical model used to estimate the mount axis.

## Celestial sphere and polar axis
We will describe positions on the celestial sphere using normalized 3-D Cartesian vectors

$$
\mathbf{x} = \begin{pmatrix}
x \\
y \\
z
\end{pmatrix},
$$

where $\|\mathbf{x}\| = 1$.
One way to parametrize the celestial sphere is using local alt-azimuthal coordinates.
Using the altitude $A$ and azimuth $\alpha$, we can write

$$
\mathbf{x} = \begin{pmatrix}
\cos(A)\cos(\alpha) \\
-\cos(A)\sin(\alpha) \\
\sin(A)
\end{pmatrix},
$$

where $A = \arcsin(z)$ and $\alpha = \arctan(-y/x)$.

In these coordinates, the $x$ axis points North, the $y$ axis points West, and the $z$ axis points to the
zenith. Azimuth $\alpha$ increases to the East.

The polar axis is given by $A_n = \phi_{obs}$, $\alpha_n = 0$, where $\phi_{obs}$ is the local 
latitude of the observer. In Cartesian coordinates, it is

$$
\mathbf{n} = \begin{pmatrix}
\cos(\phi_{obs}) \\
0 \\
\sin(\phi_{obs})
\end{pmatrix}.
$$

## Apparent motion of the stars
Due to the rotation of the Earth, a star at initial apparent position $\mathbf{x}_0$ is rotated on the
celestial sphere about the polar axis $\mathbf{n}$ with an angular velocity of $\Omega = 2\pi/T$,
where $T\approx 86164.0905\,s$ is the period of one sidereal day.

This motion is described by a rotation matrix 

$$
R(\mathbf{n}, \theta) = \mathbb{1} + \sin(\theta) K_{\mathbf{n}} +(1 - \cos(\theta)) K_{\mathbf{n}}^2.
$$

This is the Rodrigues rotation formula, where the antisymmetric matrix

$$
K_{\mathbf{n}} = \begin{pmatrix}
    0 & -n_z & n_y \\
    n_z & 0 & -n_x \\
    -n_y & n_x & 0
\end{pmatrix}
$$

is the infinitesimal generator of rotations about the axis $\mathbf{n}$.
NB: This formula can be rewritten as $R(\mathbf{n}, \theta) = \exp(\theta K_{\mathbf{n}})$ using a matrix
exponential.

Armed with this, we can finally write the apparent position of a star after a time $t$ as

$$
\mathbf{x}(t) = R(\mathbf{n}, -\Omega t) \mathbf{x}_0,
$$

where the minus sign comes from the fact that the Earth rotates counterclockwise as seen from the North celestial pole, and thus the celestial sphere rotates clockwise.

## Fitting the mount axis

Let's say we have an imperfectly polar-aligned mount that nonetheless rotates at $\Omega=2\pi/T$ (we are neglecting periodic error and similar issues).

If the mount axis is given by the normal vector $\mathbf{m}$, then initially pointing the camera it at a position
$\mathbf{x}_0$ on the celestial sphere, after a time $t$ the camera will point at

$$
\mathbf{x}(t) = R(\mathbf{m}, -\Omega t) \mathbf{x}_0.
$$

Note that the camera simply rotates about an axis different from the polar axis. A perfectly polar-aligned
mount has $\mathbf{m} = \mathbf{n}$.

In the astrotools code, we use plate-solving to obtain a set of measurements of camera pointing positions

$$
\mathbf{x}_1, \mathbf{x}_2, \mathbf{x}_3, \dots, \mathbf{x}_N
$$

and corresponding times $t_{12}, t_{23}, \dots, t_{N-1,N}$ such that $\mathbf{x}_{i+1} = R(\mathbf{m}, -\Omega
t_{i,i+1}) \mathbf{x}_i$.

We want to find the mount axis $\mathbf{m}$ to know how misaligned we are. To do this, we solve
the nonlinear least-squares problem which corresponds to minimizing the sum of squared errors

$$
E(\mathbf{m}) = \frac{1}{2} \sum_{i=1}^{N-1} \| \mathbf{x}_{i+1} - R(\mathbf{m}, -\Omega
t_{i,i+1})\mathbf{x}_i \|^2.
$$

The minimization is done numerically and the mount axis is parametrized by its altitude $A_m$ and azimuth
$\alpha_m$.

## Estimating drift velocity

Once the mount axis is known, we can estimate the velocity with which the camera's field drifts.
First we note that if we are looking at a position $\mathbf{x}(t) = R(\mathbf{n},-\Omega t)\mathbf{x}_0$ on the celestial
sphere, then the apparent motion in the frame of reference of the camera is given by

$$
\mathbf{x}_c(t) = R(\mathbf{m},\Omega t) R(\mathbf{n},-\Omega t)\mathbf{x}_{c,0}.
$$

We approximate the drift by looking at the angle between the camera positions after a time $t$,

$$
\cos\theta(t) = \mathbf{x}_c(t) \cdot \mathbf{x}_{c,0}.
$$

Taylor expanding for small $t$, we arrive at

$$
|\theta(t)| \approx v t, \qquad v = \Omega \|(\mathbf{m} - \mathbf{n})\times \mathbf{x}_{c,0}\|,
$$

where $v$ is the angular drift velocity.

### TODO: Error analysis for drift velocity

## Best rotation to match vectors

Let's assume we have two unit vectors, say the mount axis $\mathbf{m}$ and the polar axis $\mathbf{n}$,
and we would like to find the best rotation $R(\mathbf{q},\theta)$ for a given rotation axis $\mathbf{q}$ (This could be, say the Alt or the Az axis of the mount) to rotate $\mathbf{m}$ onto $\mathbf{n}$.
The motivation for this is to develop a direct feedback method for polar alignment.

We want to solve the optimization problem

$$
\operatorname{min}_{\theta} \frac{1}{2} \| e^{\theta K_\mathbf{q}}\mathbf{m} - \mathbf{n} \|^2.
$$

Taking the derivative of the objective with respect to $\theta$ we obtain,

\begin{align}
\frac{\partial \mathcal{L}}{\partial\theta} &= \frac{1}{2} (K_\mathbf{q} e^{\theta K_\mathbf{q}} \mathbf{m} )^\top( e^{\theta K_\mathbf{q}}\mathbf{m} - \mathbf{n})
    + \frac{1}{2} ( e^{\theta K_\mathbf{q}}\mathbf{m} - \mathbf{n})^\top (K_\mathbf{q} e^{\theta K_\mathbf{q}} \mathbf{m} ) \\
    &= \frac{1}{2} \left( \mathbf{m}^\top K_q^\top e^{-\theta K_\mathbf{q}} e^{\theta K_\mathbf{q}}\mathbf{m}
    - \mathbf{m}^\top K_q^\top e^{-\theta K_\mathbf{q}} \mathbf{n} \right)
    + \frac{1}{2} \left( \mathbf{m}^\top e^{-\theta K_\mathbf{q}} e^{\theta K_\mathbf{q}} K_{\mathbf{q}}\mathbf{m}
    -\mathbf{n}^\top e^{\theta K_\mathbf{q}} K_{\mathbf{q}}\mathbf{m}  \right) \\
    &= \mathbf{m}^\top K_\mathbf{q} e^{-\theta K_\mathbf{q}} \mathbf{n},
\end{align}

where we used that $K_{\mathbf{q}}^\top = - K_{\mathbf{q}}$. We further simplify this by plugging in
the Rodrigues formula,

\begin{align}
    \frac{\partial \mathcal{L}}{\partial\theta} &= \mathbf{m}^\top K_\mathbf{q} e^{-\theta K_\mathbf{q}} \mathbf{n} \\
    &= \mathbf{m}^\top K_\mathbf{q} \left( \mathbb{1} - \sin\theta K_\mathbf{q} + (1-\cos\theta)K_\mathbf{q}^2 \right)\mathbf{n} \\
    &= \mathbf{m}^\top K_\mathbf{q} \mathbf{n} - \sin\theta\, \mathbf{m}^\top K_\mathbf{q}^2 \mathbf{n}
    + \mathbf{m}^\top K_\mathbf{q}^3 \mathbf{n} - \cos\theta\, \mathbf{m}^\top K_\mathbf{q}^3 \mathbf{n} \\
    &= - \sin\theta\, \mathbf{m}^\top K_\mathbf{q}^2 \mathbf{n} + \cos\theta\, \mathbf{m}^\top K_\mathbf{q}\mathbf{n},
\end{align}

where we used that $K_\mathbf{q}^3 = - \|\mathbf{q}\|^2 K_\mathbf{q}$ and that $\|\mathbf{q}\|^2=1$.
To obtain our desired minimizer we finally set the derivative to zero,

\begin{align}
\frac{\partial \mathcal{L}}{\partial\theta} &= 0 \\
\Leftrightarrow  \sin\theta\, \mathbf{m}^\top K_\mathbf{q}^2 \mathbf{n} &= \cos\theta\, \mathbf{m}^\top K_\mathbf{q}\mathbf{n} \\
\Rightarrow \tan\theta &= \frac{\mathbf{m}^\top K_\mathbf{q}\mathbf{n}}{\mathbf{m}^\top K_\mathbf{q}^2 \mathbf{n}}
= \frac{\mathbf{m}\cdot (\mathbf{q} \times \mathbf{n})}{\mathbf{m}\cdot\left[ \mathbf{q}\times(\mathbf{q}\times \mathbf{n}) \right]}.
\end{align}

Out of the two values of $\theta$ that solve this equation, the smaller one is to be preferred.

The using the formulas for $\sin\arctan(x)$ and $\cos\arctan(x)$ we can write down the optimal rotation matrix
as 
$$
R(\mathbf{q},\theta) = \mathbb{1} + \frac{\mathbf{m}^\top K_\mathbf{q} \mathbf{n}}{\sqrt{(\mathbf{m}^\top K_\mathbf{q} \mathbf{n})^2 + (\mathbf{m}^\top K_\mathbf{q}^2 \mathbf{n})^2}} K_\mathbf{q}
+ \left(1 - \frac{\mathbf{m}^\top K_\mathbf{q}^2 \mathbf{n}}{\sqrt{(\mathbf{m}^\top K_\mathbf{q} \mathbf{n})^2 + (\mathbf{m}^\top K_\mathbf{q}^2 \mathbf{n})^2}} \right) K_\mathbf{q}^2.
$$

### Rotation axes corresponding to mount alt and az axes.
We work in the local alt-azimuthal frame established above.
Assuming that the mount is level, its azimuth axis always corresponds simply to
$$
\mathbf{q}_{az} = \begin{pmatrix}0\\0\\1\end{pmatrix}.
$$

The altitude axis is more complicated. It must lie in the $x$-$y$ plane and is perpendicular to the
mount axis $\mathbf{m}$. Projecting the mount axis onto the plane and normalizing we find
$$
\tilde{\mathbf{m}} = \frac{1}{\sqrt{\cos(A)^2}}\begin{pmatrix}\cos(A)\cos(\alpha)\\-\cos(A)\sin(\alpha)\\0\end{pmatrix} = \begin{pmatrix}\cos(\alpha)\\-\sin(\alpha)\\0\end{pmatrix},
$$
where we used that by definition of the reference frame, $0\leq A \leq 90^\circ$.
Rotating this by $90^\circ$ in the plane we finally find
$$
\mathbf{q}_{alt} = \begin{pmatrix}-\sin(\alpha)\\-\cos(\alpha)\\0\end{pmatrix}.
$$
The minus signs make sure that positive rotation angles $\theta$ lead to rotation upwards.
While the azimuth axis is fixed, the altitude axis depends on mount orientation. If the mount is not levelled properly, additional complications arise.

It should be noted that the rotation matrices corresponding to adjustments in Alt and Az do *not* commute. 
In general, optimal adjustment directions must be computed and executed one after the other, alternating
in Alt and Az.


```python

```

a3940e22df308e2405396bdcc1c04e707e9c2b1a
ipynb
Jupyter Notebook
MathematicalModel.ipynb
hronellenfitsch/astrotools
0aaa5a823f754734ee40cf98acac80e9d270a2da
MIT
MathematicalModel.ipynb
hronellenfitsch/astrotools
0aaa5a823f754734ee40cf98acac80e9d270a2da
MIT
MathematicalModel.ipynb
hronellenfitsch/astrotools
0aaa5a823f754734ee40cf98acac80e9d270a2da
MIT
Qwen/Qwen-72B
1. YES
2. YES