Spherical Kinematics

Here we detail the algorithm for computing the propagation of the position of an asteroid given observations. This method uses the laws of kinematic motion in spherical coordinates.

Deprecated since version 2022.2.1: This method has not been implemented nor have the calculations below been proven to be correct and work for the purposes of being a PropagationEngine.

Definitions

We are using spherical coordinates, namely the coordinates \((r, \theta, \phi)\) for the radial distance, polar angle, and azimuthal angle respectively. In relation to astrometric coordinates RA and DEC, \((\alpha, \delta)\):

\[\alpha = \phi \qquad \delta = \frac{\pi}{2} - \theta\]

Observations taken by the Opihi telescoped and processed by OpihiExarata are represented as \((\alpha_n, \delta_n, t_n)\). These correspond to the temporal spherical coordinates \((r=1, \phi_n, \delta_n, t_n)\). For \(t\) is the absolute time of the observation; UNIX time or Julian date time works best.

Moreover, in spherical coordinates, the position, velocity, and acceleration vectors are given as: (See Keplerian Ellipses Chapter 2 Reed 2019.)

\[\begin{split}\mathbf{r} &= r \hat{\mathbf{r}} \\ \mathbf{v} &= \dot{r} \hat{\mathbf{r}} + r \dot\theta \hat{\boldsymbol\theta } + r \dot\phi \sin\theta \hat{\boldsymbol\phi} \\ \mathbf{a} &= \left(\ddot{r} - r\dot\theta^2 - r\dot\phi^2\sin^2\theta \right) \hat{\mathbf{r}} \\ &\quad + \left( r\ddot\theta + 2\dot{r}\dot\theta - r\dot\phi^2\sin\theta\cos\theta \right) \hat{\boldsymbol\theta } \\ &\quad + \left( r\ddot\phi\sin\theta + 2\dot{r}\dot\phi\sin\theta + 2 r\dot\theta\dot\phi\cos\theta \right) \hat{\boldsymbol\phi}\end{split}\]

Where the basis vectors of the spherical coordinates are:

\[\begin{split}\hat{\mathbf r} &= \sin\theta \cos\phi \hat{\mathbf x} + \sin\theta \sin\phi \hat{\mathbf y} + \cos\theta \hat{\mathbf z} \\ \hat{\boldsymbol\theta} &= \cos\theta \cos\phi \hat{\mathbf x} + \cos\theta \sin\phi \hat{\mathbf y} - \sin\theta \hat{\mathbf z} \\ \hat{\boldsymbol\phi} &= - \sin\phi \hat{\mathbf x} + \cos\phi \hat{\mathbf y} + 0 \hat{\mathbf z}\end{split}\]

However, for the purposes of finding by propagation, we only care about the location of the asteroid on the sky as determined by the celestial coordinates. The distance from the origin of the coordinate system does not change. Thus:

\[r = 1 \qquad \dot{r} = \ddot{r} = 0\]

And thus the kinematic vectors are:

\[\begin{split}\mathbf{r} &= \hat{\mathbf{r}} \\ \mathbf{v} &= \dot\theta \hat{\boldsymbol\theta } + \dot\phi \sin\theta \hat{\boldsymbol\phi} \\ \mathbf{a} &= \left(-\dot\theta^2 - \dot\phi^2\sin^2\theta \right) \hat{\mathbf{r}} + \left(\ddot\theta - \dot\phi^2\sin\theta\cos\theta \right) \hat{\boldsymbol\theta } + \left(\ddot\phi\sin\theta + 2 \dot\theta\dot\phi\cos\theta \right) \hat{\boldsymbol\phi}\end{split}\]

We can convert these vectors to Cartesian coordinates using the following matrix transformation. The spherical to Cartesian transformation matrix \(\mathbf{R}\) can be derived from the defining angles of the spherical coordinate system \(\theta\) and \(\phi\). Where \(\mathbf{u}_\text{cart} = \mathbf{R} \mathbf{u}_\text{sph}\).:

\[\begin{split}\mathbf{R} = \begin{bmatrix} \sin\theta\cos\phi & \cos\theta\cos\phi & -\sin\phi \\ \sin\theta\sin\phi & \cos\theta\sin\phi & \cos\phi \\ \cos\theta & -\sin\theta & 0 \end{bmatrix}\end{split}\]

Thus, in Cartesian coordinates:

\[\begin{split}\mathbf{r} &= \mathbf{R} \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} \\ \mathbf{v} &= \mathbf{R} \begin{bmatrix} 0 \\ \dot\theta \\ \dot\phi \sin\theta \end{bmatrix} \\ \mathbf{a} &= \mathbf{R} \begin{bmatrix} -\dot\theta^2 - \dot\phi^2\sin^2\theta \\ \ddot\theta - \dot\phi^2\sin\theta\cos\theta \\ \ddot\phi\sin\theta + 2 \dot\theta\dot\phi\cos\theta \end{bmatrix}\end{split}\]

Generally, as kinematics are defined based on the time derivatives, the general equation of motion for constant acceleration can be used to determine the movement of the position vector representing the asteroid:

\[\mathbf{r} \triangleq \mathbf{r} \qquad \mathbf{v} \triangleq \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} \qquad \mathbf{a} \triangleq \frac{\mathrm{d}^2\mathbf{r}}{\mathrm{d}t^2} \qquad \mathbf{j} \triangleq \frac{\mathrm{d}^3\mathbf{r}}{\mathrm{d}t^3} = 0\]

\[\mathbf{r} = \mathbf{r}_0 + \mathbf{v}_0 \tau + \frac{1}{2} \mathbf{a} \tau^2\]

Where \(\tau\) is the time interval between the defining time of observation to the current time.

Deriving Rates

Multiple observations from Opihi provides multiple sightings of an asteroid at many different points in the sky, providing multiple RA and DEC coordinates, \(\alpha_n\) and \(\delta_n\) at time \(t_n\). We have a total of \(N\) RA-DEC observations. (The propagation calculation will need to be redone for a new observation set \(N' = N + 1\).)

We can convert this to spherical coordinates with \(\phi_n = \alpha_n\) and \(\theta_n = \frac{\pi}{2} - \delta_n\).

These multiple observations allows for the determination of the rates of change of spherical coordinates for the asteroid, namely: (For the time difference \(t_\Delta = t_{n+1} - t_n\).)

\[ \begin{align}\begin{aligned}\dot\theta_p = \frac{\theta_{n+1} - \theta_{n}}{t_{n+1} - t_n}\\\dot\phi_p = \frac{\phi_{n+1} - \phi_{n}}{t_{n+1} - t_n}\\t'_p = \frac{1}{2} \left( t_{n+1} + t_n \right)\end{aligned}\end{align} \]

and,

\[ \begin{align}\begin{aligned}\ddot\theta_q = \frac{\dot\theta_{p+1} - \dot\theta_{p}}{t'_{p+1} - t'_p}\\\ddot\phi_q = \frac{\dot\phi_{p+1} - \dot\phi_{p}}{t'_{p+1} - t'_p}\end{aligned}\end{align} \]

The first order rates changes over time. As such, it is required that two observations be reserved as special observations which the first order rates are calculated and to established the spherical coordinate system itself. Although it does not need to be the first two observations, it is often connivent to use them. As such, using the first two observations \(n=0\) and \(n=1\), we have:

\[\begin{split}\theta &= \theta_0 \\ \phi &= \phi_0 \\ \dot\theta &= \dot\theta_0 = \frac{\theta_1 - \theta_0}{t_1 - t_0} \\ \dot\phi &= \dot\phi_0 = \frac{\phi_1 - \phi_0}{t_1 - t_0} \\\end{split}\]

Because we assume constant acceleration (\(\mathbf{j} = 0\)), the second differential values are assumed to be constant and thus an average is more representational of the value. (A mean or median is valid.)

\[ \begin{align}\begin{aligned}\ddot\theta = \frac{1}{Q} \sum_q^Q \ddot\theta_q \approx \mathrm{median} (\ddot\theta_q)\\\ddot\phi = \frac{1}{Q} \sum_q^Q \ddot\phi_q \approx \mathrm{median} (\ddot\phi_q)\end{aligned}\end{align} \]

In the case for \(N=2\), then the total number of derived angular first order rates is \(P=1\). As such the second order rates cannot be calculated and \(Q=0\) (the cardinality of the arrays are zero). By default, for this special case:

\[\#(\ddot\theta_q) = \#(\ddot\phi_q) = 0 \implies Q = 0 \longrightarrow \ddot\theta = 0 \quad \ddot\phi = 0\]

Spherical Motion

With the 0th, 1st, and 2nd order rates calculated from the set of \(N\) observations, the kinematic vectors can be calculated. The special observations defining the coordinate system and the velocities also define the initial vectors from which kinematics shall be applied to. The acceleration vector, being constant means \(\mathbf{a}_0 = \mathbf{a}\). Namely, these vectors are, in Cartesian coordinates,

\[\begin{split}\mathbf{r_0} &= \begin{bmatrix} \sin\theta\cos\phi & \cos\theta\cos\phi & -\sin\phi \\ \sin\theta\sin\phi & \cos\theta\sin\phi & \cos\phi \\ \cos\theta & -\sin\theta & 0 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} \\ \mathbf{v_0} &= \begin{bmatrix} \sin\theta\cos\phi & \cos\theta\cos\phi & -\sin\phi \\ \sin\theta\sin\phi & \cos\theta\sin\phi & \cos\phi \\ \cos\theta & -\sin\theta & 0 \end{bmatrix} \begin{bmatrix} 0 \\ \dot\theta \\ \dot\phi \sin\theta \end{bmatrix} \\ \mathbf{a} &= \begin{bmatrix} \sin\theta\cos\phi & \cos\theta\cos\phi & -\sin\phi \\ \sin\theta\sin\phi & \cos\theta\sin\phi & \cos\phi \\ \cos\theta & -\sin\theta & 0 \end{bmatrix} \begin{bmatrix} -\dot\theta^2 - \dot\phi^2\sin^2\theta \\ \ddot\theta - \dot\phi^2\sin\theta\cos\theta \\ \ddot\phi\sin\theta + 2 \dot\theta\dot\phi\cos\theta \end{bmatrix}\end{split}\]

All three of these vectors are constant in future time. The position at a set of future observations at time(s) \(t^+_i\) can be calculated using the kinematic equation; the time intervals \(\tau_i\) being \(\tau_i = t^+_i - t_0\):

\[\mathbf{r}^+_i = \mathbf{r}_0 + \mathbf{v}_0 \left(t^+_i - t_0\right) + \frac{1}{2} \mathbf{a} \left(t^+_i - t_0\right)^2\]

Celestial Sphere

These new future position vectors \(\mathbf{r}^+_i\) are in Cartesian coordinates. The calculations should be done in Cartesian, provided the conversion earlier.

Each position vector can be represented as:

\[\begin{split}\mathbf{r}^+_i = X_i \hat{\mathbf{x}} + Y_i \hat{\mathbf{y}} + Z_i \hat{\mathbf{z}} = \begin{bmatrix} X_i \\ Y_i \\ Z_i \end{bmatrix}\end{split}\]

These Cartesian coordinate position vectors, centered on the origin, represents where the asteroid is on the celestial sphere in the future at an observation time of \(t^+_i\). From these Cartesian coordinates, we can extract their location in spherical coordinates,

\[\begin{split}r^+_i &= \sqrt{X_i^2 + Y_i^2 + Z_i^2} \\ \theta^+_i &= \arccos\left(\frac{Z_i}{r^+_i}\right) = \arccos\left(\frac{Z_i}{\sqrt{X_i^2 + Y_i^2 + Z_i^2}}\right) \\ \phi^+_i &= \arctan\!2(Y_i, X_i) \simeq \arctan\left(\frac{Y_i}{X_i}\right)\end{split}\]

Note

In order to properly handle the quadrant issue, the 2-argument arctangent is required. Moreover, if the 2-argument arctangent function returns in a range \(-\pi \leq \angle \leq \pi\), it can be converted to the usual range of \(0 \leq \phi \leq 2\pi\) with: \(\phi = \angle \mod 2\pi\) or \(\phi = \angle \mod 360^\circ\)

These spherical coordinate locations can then be converted into future RA and DEC temporal coordinates \((\alpha^+_i, \delta^+_i, t^+_i)\):

\[\begin{split}\alpha^+_i &= \phi^+_i \\ \delta^+_i &= \frac{\pi}{2} - \theta^+_i \\ t^+_i &= t^+_i\end{split}\]

Lemmas

Derivation of Vector Equation of Motion

Newton’s second law and constant acceleration stipulates:

\[\mathbf{F} = m \mathbf{a} = m \ddot{\mathbf{r}} \qquad \dot{\mathbf{F}} = 0\]

This thus provides the differential equation of motion (For constant \(\mathbf{F}\).)

\[\ddot{\mathbf{r}} = \frac{\mathrm{d}^2\mathbf{r}}{\mathrm{d}t^2} = \frac{\mathbf{F}}{m}\]

We define based on the laws of integrations (and in essence the fundamental theorem of calculus):

\[\dot{\mathbf{f}} \triangleq \frac{\mathrm{d}\mathbf{f}}{\mathrm{d}t} \Longleftrightarrow \int \dot{\mathbf{f}} \mathrm{d} t = \mathbf{f} + \mathbf{C}\]

\[\int \frac{\mathrm{d}\mathbf{f}}{\mathrm{d}t} \mathrm{d} t = \mathbf{f}\]

We can solve the differential equation of motion:

\[\begin{split}\ddot{\mathbf{r}} &= \frac{\mathbf{F}}{m} \\ \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} \right) &= \frac{\mathbf{F}}{m} \\ \int \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} \right) \mathrm{d}t &= \int \frac{\mathbf{F}}{m} \mathrm{d}t = \frac{\mathbf{F}}{m} \int 1 \mathrm{d}t = \frac{\mathbf{F}}{m} t + \mathbf{C_1} \\ \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} &= \frac{\mathbf{F}}{m} t + \mathbf{C_1} \\ \int \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} \mathrm{d}t &= \int \frac{\mathbf{F}}{m} t + \mathbf{C_1} \mathrm{d}t = \frac{\mathbf{F}}{m} \int t \mathrm{d}t + \int \mathbf{C_1} \mathrm{d}t = \frac{\mathbf{F}}{m} \frac{1}{2} t^2 + \mathbf{C_1} t + \mathbf{C_2} \\ \mathbf{r} &= \frac{\mathbf{F}}{m} \frac{1}{2} t^2 + \mathbf{C_1} t + \mathbf{C_2}\end{split}\]

For the initial conditions:

\[\begin{split}t = 0 &\implies \mathbf{r} = \mathbf{C_2} = \mathbf{r_0} \\ t = 0 &\implies \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} = \mathbf{C_1} = \mathbf{v_0} \\ \dot{\mathbf{F}} = 0 &\implies \frac{\mathrm{d}^2\mathbf{r}}{\mathrm{d}t^2} = \frac{\mathbf{F}}{m} = \mathbf{a_0} = \mathbf{a}\end{split}\]

Thus, the total valid solution is:

\[\mathbf{r} = \mathbf{r_0} + \mathbf{v_0} t + \frac{1}{2} \mathbf{a} t^2\]