trace a geodesic between two points on a sphere

Question

I'd like to draw a geodesic onto a array of pixels. I already have a way to draw a pixel at a given longitude latitude, and I am given the longitude and latitude of each end point. What is the parametric equation for tracing a geodesic between two given points?

Some have suggested a way to find the distance between two points. I'm not looking for a distance, rather some sort of parameterization I can iterate through which will generate successive points along the great circle joining the two given points.

Answer 1

I'd transform the points to Cartesian coordinates, consider them as vectors, form their cross product to determine an axis of rotation, and rotate one into the other about this axis (using the axis-angle representation of rotations). If necessary, you can transform the resulting Cartesian coordinates back to longitudes and latitudes.

Answer 2

Let's say the coordinates of the pixels are $(\theta_i, \phi_i)$. $\theta$ is measured from the equator, so it varies between $\pm \pi/2$ and $\phi$ is measured in the horizontal plane, from $0$ to $2\pi$. This is like the standard spherical coordinate system shifted by $\pi/2$. Then the unit vectors in Cartesian coordinates for a point are described by: $$\begin{align}x_i&=\cos\theta_i\cos\phi_i\\y_i&=\cos\theta_i\sin\phi_i\\z_i&=\sin\theta_i\end{align}$$ As mentioned in another answer, we need to find a rotation axis. For that, we use the cross product: $$\vec n=\vec P_1\times\vec P_2$$ By components, this becomes: $$\begin{align}n_x&=y_1z_2-z_1y_2\\n_y&=z_1x_2-x_1z_2\\n_z&=x_1y_2-y_1x_2\end{align}$$ Note that $\vec n$ doess not necessarily have length $1$, so you will need to use $\hat n=\frac{\vec n}{|\vec n|}$. Also, $\vec n$ might be $0$. Then you have two points either identical or on diametrically opposite positions. In that case there are infinitely many rotation axes you can choose, in the perpendicular plane to the diameter. So for example, if the two points are the North and South poles, then any meridian would work. In that case, I would calculate the angle between the diameter and each of the $x, y, z$ axes, choose the largest one of them (or the first of the largest, if they are equal), and use the Gram-Schmidt orthogonalization process to find $\hat n$. The reason I would choose the largest is that it will introduce less numerical errors.

Note that $|\vec n|=|\vec P_1||\vec P_2|\sin\alpha=\sin\alpha$., where $\alpha$ is the angle between the vectors. If $|\vec n|=0$, you need to check if the points are identical $\alpha =0$ or opposite $\alpha =\pi$.

With the given information, we can now proceed to find the intermediate points. for that, we create an array of $\alpha_i$ from $0$ to $\alpha$, and we rotate the vector $\vec P_1$ towards vector $\vec P_2$ in the plane perpendicular to $\hat n$ using all $\alpha_i$ rotation angles. You can use the formula on wikipedia for the rotation matrix: $$R(\hat n, \alpha_i)=\begin{pmatrix}\cos\alpha_i+u_x^2(1-\cos\alpha_i)&u_xu_y(1-\cos\alpha_i)-u_z\sin\alpha_i&u_xu_z(1-\cos\alpha_i)+u_y\sin\alpha_i\\u_yu_x(1-\cos\alpha_i)+u_z\sin\alpha_i&\cos\alpha_i+u_y^2(1-\cos\alpha_i)&u_yu_z(1-\cos\alpha_i)-u_x\sin\alpha_i\\u_xu_z(1-\cos\alpha_i)-u_y\sin\alpha_i&u_zu_y(1-\cos\alpha_i)+u_x\sin\alpha_i&\cos\alpha_i+u_z^2(1-\cos\alpha_i)\end{pmatrix}$$

I've used the components of $\hat n$: $$u_j=\frac{n_j}{|\vec n|}$$Then the intermediate points are $$\vec P_i=R(\hat n,\alpha_i)\vec P_1$$This is a simple matrix multiplication for each point. You can now calculate back the $\theta_i$ and $\phi_i$ coordinates, if you want to use those.