How to calculate 3D position and update it

Question

I'm trying to calculate the new position of an object moving over the surface of a 3D sphere after some time Δt, from X1 to X2 always onto the surface of the sphere. The measurements/data I have of such object's displacement and sphere are:

• Translational velocity (cm/seg)
• Rotational velocity (deg/seg)
• Direction of displacement (deg)
• Radius of the sphere (cm) Here is an example of what I mean.

Then the idea is to update for each Δt the position of the object on the sphere. I suspect some linear algebra is needed but pretty lost at the moment. Any help would be greatly appreciated!

This is how the data frame looks like

time	translational velocity	rotational velocity	direction
0.01	2.36	6.45	0.78
0.02	1.12	1.19	0.62
0.03	1.67	4.45	1.51
0.04	1.25	5.39	1.67

Answer 1

Let's say you have a sphere of radius $r$ centered at origin, and point $\vec{p} = (p_x, p_y, p_z)$ somewhere on the sphere so $\lVert \vec{p} \rVert = \sqrt{ \vec{p} \cdot \vec{p} } = \sqrt{p_x^2 + p_y^2 + p_z^2} = r$.

Let's say that at time $t = 0$ the point starts traveling on the sphere in direction $\vec{v}$ with velocity $\lVert\vec{v}\rVert$. Because $\vec{v}$ is tangential to the sphere, $\vec{v} \perp \vec{p}$. If you pick a random direction $\vec{w}$ in 3D, you can take just its tangential component, $$\vec{v} = \vec{w} - \frac{\vec{p} \cdot \vec{w}}{\vec{p} \cdot \vec{p}} \vec{p}$$ essentially performing a single step of the Gram-Schmidt process, making $\vec{v}$ the part of $\vec{w}$ that is orthogonal to $\vec{p}$ and thus planar in the sphere tangent plane at $\vec{p}$.

If the point travels in a local straight line, it will travel along a great circle (a circle with its center at the center of the sphere, starting point on the perimeter, and the initial direction in the plane of the circle). The simplest way to express this is by rotation around axis $\hat{a}$ perpendicular to the great circle, $$\hat{a} = \frac{\vec{v} \times \vec{p}}{\left\lVert \vec{v} \times \vec{p} \right\rVert} \tag{1}\label{G1}$$ If $v = \lVert\vec{v}\rVert$ is the (2D) velocity, the angular velocity $\omega$ is $$\omega = \frac{v}{2 \pi r} = \frac{v}{2 \pi \lVert \vec{p} \rVert} \tag{2}\label{G2}$$ and we can describe the trajectory of $\vec{p}$ as a function of time $t$ via $$\vec{p}(t) = \mathbf{R}(\hat{a}, \omega t) \vec{p}_0 \tag{3a}\label{G3a}$$ where $\vec{p}_0 = \vec{p}(0)$ is the location at time $t = 0$, and $$\begin{aligned} \mathbf{R}(\hat{a}, \omega t) & = \left[ \begin{matrix} c + a_x^2 u & a_x a_y u - a_z s & a_x a_z u + a_y s \\ a_y a_x u + a_z s & c + a_y^2 u & a_y a_z u - a_x s \\ a_z a_x u - a_y s & a_z a_y u + a_x s & c + a_z^2 u \\ \end{matrix} \right] \\ c & = \cos(\omega t) \\ s & = \sin(\omega t) \\ u & = 1 - \cos(\omega t) = 1 - c \\ \end{aligned} \tag{3b}\label{G3b}$$ Alternatively (and equivalently), you can use Rodrigues' rotation formula, $$\begin{aligned} \vec{p}(t) &= \vec{p}_0 c + \left(\hat{a} \times \vec{p}_0\right) s + \hat{a} \left( \hat{a} \cdot \vec{p}_0 \right) (1 - c) \\ c &= \cos(\omega t) \\ s &= \sin(\omega t) \\ \end{aligned} \tag{3c}\label{G3c}$$ In both cases, we just rotate $\vec{p}_0$ around unit axis vector $\hat{a}$ by angle $\omega t$.

To find out the instant direction, the same applies to $\vec{v}$: just replace $\vec{p}_0$ by $\vec{v}_0$ in the formula you wish to use, and you get the exact local direction vector $\vec{v}(t)$ as a function of time $t$.

(By notation $\vec{p}(t)$, I am referring to a vector-valued function in a single variable $t$. While the latter looks simpler, if you need both $\vec{p}(t)$ and $\vec{v}(t)$, the matrix form may require fewer basic arithmetic operations.)

---

Let's say you know the starting point $\vec{p}_0$, and the destination point $\vec{p}_1$, both on the sphere ($\lVert \vec{p}_0 \rVert = \lVert \vec{p}_1 \rVert = r$).

This time, the unit axis vector $\hat{a}$ is $$\hat{a} = \frac{\vec{p}_0 \times \vec{p}_1}{\lVert \vec{p}_0 \times \vec{p}_1\rVert} \tag{4}\label{G4}$$
The total travel distance $L$ (measured along the surface of the sphere) is $$L = 2 \pi r \arccos\left(\theta_\Delta\right) \tag{4a}\label{G4a}$$ where $\theta_\Delta$ is the angular separation, $$\theta_\Delta = \frac{\vec{p}_0 \cdot \vec{p}_1}{\lVert\vec{p}_0\rVert \, \lVert\vec{p}_1\rVert} \tag{4b}\label{G4b}$$ and therefore the arrival time $t$ with (sphere surface) velocity $v$ is $$t = \frac{L}{v} \tag{4c}\label{G4c}$$ and the angular velocity $\omega$ is $$\omega = \frac{\theta_\Delta}{t} \tag{4d}\label{G4d}$$ Alternatively, if you know the arrival time $t$, then the (sphere surface) velocity $v$ is $$v = \frac{L}{t} \tag{4e}\label{G4e}$$

---

Treating the current position as an $r$-length Cartesian vector yields a simple simulation. However, it also means that one can easily adjust the altitude of the point also. For example, to move $\vec{p}$ to altitude $h$ just scales it to length $r + h$, $$\vec{p}^\prime = \frac{r + h}{\lVert \vec{p} \rVert} \vec{p} \tag{5}\label{G5}$$ and to change the altitude of $\vec{p}$ by $h$ is just $$\vec{p}^\prime = \frac{\lVert \vec{p} \rVert + h}{\lVert \vec{p} \rVert} \vec{p} \tag{6}\label{G6}$$ where the$\,^\prime$ is there just to denote that it is the point after the altitude adjustment.

Furthermore, you can use 3D squared distance between two points to determine collisions et cetera. The velocity $v$ is similarly measured along the sphere surface, as opposed to along the elliptical/non-spherical curve if altitude changes would be taken into account; just as we do with air travel on Earth.

Finally, if $z$ axis is towards the north pole and $x$ axis towards zero meridian on the equator, then $$\begin{aligned} \text{latitude} &= \displaystyle \arcsin\left(\frac{p_z}{\lVert\vec{p}\rVert}\right) \\ \text{longitude} &= \displaystyle \operatorname{atan2}\left(\frac{p_y}{\lVert\vec{p}\rVert}, ~ \frac{p_x}{\lVert\vec{p}\rVert}\right) \\ \end{aligned} \tag{7}\label{G7}$$ where $\operatorname{atan2}(y, x)$ is the two-argument version of arcus tangent provided by most programming languages; similar to $\arctan(y/x)$, but takes the signs of $y$ and $x$ into account, so that the result covers all four quadrants (full 360°).

---

Conversion between a 3D direction vector $\vec{v}$ and a single direction angle $\varphi$ – compass reading, with 0° towards North, 180° towards South, and 90° West – is possible, if we handle directions at north pole (where every direction is south) and at south pole (where every direction is north) specially: there, the direction angle is towards that longitude on the equator.

If $\hat{N}$ is a 3D unit vector towards local North, and $\hat{W}$ is a 3D unit vector towards local West, then heading $\varphi$ corresponds to 3D unit direction $\hat{d}$, $$\hat{d} = \hat{N}\cos\varphi + \hat{W}\sin\varphi \tag{8}\label{G8}$$ If you are at 3D point $\vec{p} = (p_x, p_y, p_z)$ and not at a pole ($p_x^2 + p_y^2 \ne 0$), then $$\hat{N} = \frac{1}{\sqrt{(p_x^2 + p_y^2)(p_x^2 + p_y^2 + p_z^2)}} \left[ \begin{matrix} - p_x p_z \\ - p_y p_z \\ p_x^2 + p_y^2 \\ \end{matrix} \right], \quad \hat{W} = \frac{1}{\sqrt{p_x^2 + p_y^2 + p_z^2}} \left[ \begin{matrix} p_y \\ -p_x \\ p_z \\ \end{matrix} \right] \tag{9a}\label{G9a}$$ but at the poles ($p_x = p_y = 0$), $$\hat{N} = \left[ \begin{matrix} 1 \\ 0 \\ 0 \end{matrix} \right], \quad \hat{W} = \left[ \begin{matrix} 0 \\ 1 \\ 0 \end{matrix} \right] \tag{9b}\label{G9b}$$

If you have direction vector $\vec{d}$ at point $\vec{p}$, calculate the north and west unit vectors $\hat{N}$ and $\hat{W}$, then ensure the direction vector is tangential to the sphere, $$\vec{d}^\prime = \vec{d} - \displaystyle\frac{\vec{p} \cdot \vec{d}}{\vec{p} \cdot \vec{p}} \vec{p}$$ and the heading angle $\varphi$ is then $$\varphi = \operatorname{atan2}\left( \hat{W} \cdot \vec{d}^\prime , ~ \hat{N} \cdot \vec{d}^\prime \right) \tag{G9c}\label{G9c}$$

Answer 2

You need to solve Euler-Lagrange equation of motion along the geodesic curve in your metric space: \begin{equation} \frac{d^2x^{\mu}}{dt^2} + \Gamma^{\mu}_{\alpha\beta}\frac{dx^{\alpha}}{dt}\frac{dx^{\beta}}{dt} = 0, \end{equation} where t is the time, $x^{\mu}$ are coordinate components, $\Gamma^{\mu}_{\alpha\beta}$ is the Christoffel symbols of the metric. The Cristoffel symbols $\Gamma^{\mu}_{\alpha\beta}$ are given by the following formula \begin{equation} \Gamma^{\mu}_{\alpha\beta} = \frac{1}{2}g^{\mu\nu}(\frac{\partial g_{\nu\alpha}}{\partial x^{\beta}} + \frac{\partial g_{\nu\beta}}{\partial x^{\alpha}} - \frac{\partial g_{\alpha\beta}}{\partial x^{\nu}}), \end{equation} where $g_{\alpha\beta}$ is the metric tensor. $g_{\alpha\beta}$ and $g^{\beta\gamma}$ are linked by the following equation \begin{equation} g_{\alpha\beta}g^{\beta\gamma} = \delta_{\alpha}^{\,\gamma} = \delta^{\gamma}_{\,\alpha}, \end{equation} where $\delta^{\gamma}_{\,\alpha}$ is the Kronecker symbol: \begin{equation} \delta^{\gamma}_{\,\alpha} = \begin{cases} 1,\text{ if }\gamma=\alpha\\ 0,\text{ otherwise} \end{cases}. \end{equation} So the matrix $g^{\alpha\beta}$ is inverse of the matrix $g_{\alpha\beta}$. The metric tensor is the tensor, which links the length $ds$ of differential minor arc with the coordinate differentials $dx^{\alpha}$: \begin{equation} ds^2 = g_{\mu\nu}dx^{\mu}dx^{\nu}. \end{equation} In the case of the sphere \begin{equation} g = \begin{pmatrix} 1 & 0\\ 0 & \sin^2\theta, \end{pmatrix}, \end{equation} so $ds^2$ is given by the following formula \begin{equation} ds^2 = d\theta^2 + \sin^2\theta d\phi, \end{equation} where $\theta$ is a polar angle, $\phi$ is axial angle. $\theta$ and $\phi$ are the components of the vector $x^{\mu}$ in this case. The differential $dx^{mu}$ can be written in the following form $dx^{\mu} = \frac{dx^{\mu}}{dt}dt$, where $dt$ is the time differential and $\frac{dx^{\mu}}{dt}$ is the generalized speed. Let us consider the Lagrange functional of the following form \begin{equation} \mathcal{L}[x^{\mu}] = \int^{s_1}_0 ds = \int^{s_1}_{0}\sqrt{g_{\mu\nu}dx^{\mu}dx^{\nu}} = \int^{t}_{0}\sqrt{g_{\mu\nu}\frac{dx^{\mu}}{dt}dt\frac{dx^{\nu}}{dt}dt} = \int^{t}_{0}\sqrt{g_{\mu\nu}\frac{dx^{\mu}}{dt}\frac{dx^{\nu}}{dt}}dt. \end{equation} We can use the calculus of variations and find a function $dx^{\mu}(t)$ that minimizes our functional. To do this, you need to write down the Euler-Lagrange equation. In the case of this functional, it has the form \begin{equation} \frac{d^2x^{\mu}}{dt^2} + \Gamma^{\mu}_{\alpha\beta}\frac{dx^{\alpha}}{dt}\frac{dx^{\beta}}{dt} = 0. \end{equation} Since the functional $\mathcal{L}[x^{\mu}] $ has the meaning of length, the solution to the previous equation is a curve of shortest length connecting two points in a given metric space. Such curves are called geodesic curves. The non-zero Christoffel symbols in the case of sphere metric are (see the following link http://einsteinrelativelyeasy.com/index.php/general-relativity/34-christoffel-symbol-exercise-calculation-in-polar-coordinates-part-ii) \begin{equation} \Gamma^{\phi}_{\phi\theta}=\Gamma^{\phi}_{\theta\phi} = \frac{\cos\theta}{\sin\theta}, \\ \Gamma^{\theta}_{\phi\phi} = -\sin\theta\cos\theta, \end{equation} so the system of motion equation has the following form \begin{equation} \begin{cases} \frac{d^2\theta}{dt^2} - \sin\theta\cos\theta\left(\frac{d\phi}{dt}\right)^2 = 0,\\ \frac{d^2\phi}{dt^2} + 2\frac{\cos\theta}{\sin\theta}\frac{d\theta}{dt}\frac{d\phi}{dt} = 0. \end{cases} \end{equation} Let us introduce $\theta$-speed $u_{\theta}$ and $\phi$-speed $u_{\phi}$: \begin{equation} u_{\theta} = \frac{d\theta}{dt},\\ u_{\phi} = \frac{d\phi}{dt}. \end{equation} In this terms the system of motion equations has the following form: \begin{equation} \begin{cases} \frac{du_{\theta}}{dt} - \sin\theta\cos\theta u^2_{\phi} = 0,\\ \frac{du_{\phi}}{dt} + 2\frac{\cos\theta}{\sin\theta}u_{\theta}u_{\phi} = 0. \end{cases} \end{equation}

The useful information about motion along the sphere geodesics can be found here Geodesics of the Unit Sphere using Christoffel symbols.

The simplest way to simulate motion along geodesic curve is to rotate sphere in a such way that the motion will be at a constant $\phi^{\prime}$, for example, $\phi^{\prime} = 0$. In this case only $\theta^{\prime}$ changes with the time. Then calculate $\theta^{\prime}$ as $\theta^{\prime} = \omega t$. Where $\omega$ is a cyclic frequency. Then rotate sphere back into initial basis. There $\theta^{\prime}$ and $\phi^{\prime}$ are coordinates of the point object in the rotated basis.