Equally spaced points on a great circle arc between 2 points

I have 2 points on a sphere, given by $(\theta, \phi)$. How can I find $n$ equally spaced points on the smaller arc of the great circle between them, without converting the coordinates to another system?

I am sure that the smaller arc does not cross the singularity line $\phi=0,2\pi$.

coordinates = {{θ1, ϕ1}, {θ2, ϕ2}}


• "without converting the coordinates to another system" - so you don't allow conversion to Cartesian coordinates, even if it can make things a bit easier? Anyway, you might want to see this. – J. M. is away Apr 6 '17 at 16:59
• No, because I have done some after-processing of the points in the case when they lie on the arc phi = 0, 2*pi. If I convert them, and then again, I will lose this processing. – Morgan Apr 6 '17 at 17:03
• This math.stackexchange.com/questions/1476391/… gives you the cross product in spherical. That will be normal to the plane holding the great circle through your two points and seems to be a large step towards what you want. – Bill Apr 6 '17 at 23:08
• I provide a formula in mathematica.stackexchange.com/questions/180422/… but it's hideously ugly and needs to be simplified. – barrycarter Aug 22 '18 at 15:24

I remain mystified by the OP's insistence against a Cartesian conversion, since it eases the task of subdivision in this case. Nevertheless, here is a simple demonstration of "slerping" (previously shown here):

sphericalDistance[{u1_, v1_}, {u2_, v2_}] :=
InverseHaversine[Haversine[v1 - v2] + Sin[v1] Sin[v2] Haversine[u1 - u2]]

With[{n = 10, p1 = N[{π/6, π/4}], p2 = N[{-π/2, π/2}]},
Block[{d, h},
d = sphericalDistance[p1, p2]; h = Range[0, 1, 1/(n - 1)];
arcDots = (Sin[Transpose[{1 - h, h}] d]/Sin[d]).
(Append[Sin[#2] Through[{Cos, Sin}[#1]], Cos[#2]] & @@@ {p1, p2});
Graphics3D[{Sphere[], Red, Sphere[arcDots, 0.02]}, Boxed -> False]]]


If the longitude and colatitude angles are needed, here is some code for obtaining them:

{ArcTan @@ Most[#], ArcCos[Last[#]/Norm[#]]} & /@ arcDots

• Good, only in my question I have theta, then phi, so your #1 and #2 are the other way around. – Morgan Apr 25 '17 at 19:17
• @Morgan, I'll leave the conversion to your preferred convention up to you; I'm accustomed to putting the longitude before the colatitude. – J. M. is away Apr 26 '17 at 2:39

A first thought:

f[s_, t_] := {Sin[s] Cos[t], Sin[s] Sin[t], Cos[s]}
fun[{a_, b_}, {c_, d_}, n_] :=
Module[{p1 = f[a, b], p2 = f[c, d], ax, angle, pts},
ax = Cross[p1, p2];
angle = VectorAngle[p1, p2];
pts = Table[RotationMatrix[j angle/n, ax].p1, {j, 0, n, 1}];
Graphics3D[{Opacity[0.5], Sphere[], Opacity[1], Red,
PointSize[0.02], Point[{p1, p2}], PointSize[0.01], Black,
Point[pts]}]]


For example:

Manipulate[
fun[u, v,
n], {u, {0, 0}, {2 Pi, 2 Pi}}, {{v, {1, 1}}, {0, 0}, {2 Pi,
2 Pi}}, {n, {5, 10, 20}}]