Points on a geodesic line

Question

I am working on a unit sphere. I am interested to place N points on a strait line over the surface of the sphere (geodesic) between two arbitrary points. The coordinate of these points are in spherical coordinate (radians).

How do I compute a set of N equally spaced points along such line. I would like to take the curvature of the sphere into account in my calculation.

I am using python 2.7.9

Answer 1

Let us reason geometrically.

Convert the two given points to Cartesian coordinates.

The angle between the position vectors from the center to P0 and P1 is given by the dot product

cos A = P0.P1

Construct a linear combination of these:

P = (1-t).P0 + t.P1

The angle between P and P0 is given by the dot product with P normalized

cos a = cos kA/N = P.P0/|P| = ((1-t) + t.cos A)/ sqrt((1-t)² + 2.(1-t).t.cos A + t²)

Squaring and rewriting, you obtain a second degree equation in t:

cos²a.(1-t)² + 2.(1-t).t.cos²a.cos A + t².cos²a - (1-t)² - 2.(1-t).t.cos A - t².cos²A = 0

- sin²a.(1-t)² - 2.(1-t).t.sin²a.cos A - t².(cos²A - cos² a) = 0

t²(-sin²a + 2.sin²a.cos A - cos²A + cos²a) + 2.t.sin²a.(1 - cos A) - sin²a = 0

Solve the equation, compute the vector P from its definition and normalize it.

Then revert to spherical coordinates. Varying k between 1 and N-1 will give you the required intermediate points.

---

Alternatively, you can use the Rodrigue's rotation formula around an axis in 3D. The axis is given by the cross-product P0 x P1.

Answer 2

You may consider SLERP - spherical linear interpolation

P = P0*Sin(Omega*(1-t))/Sin(Omega) + P1*Sin(Omega * t)/Sin(Omega)

where Omega is central angle between start and end points (arc of great circle), t is parameter in range [0..1], for i-th point t(i) = i/N