Equally spaced points on a great circle arc between 2 points

Question

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}}

Thanks in advance.

"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.'s missing motivation, Apr 06 '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 06 '17 at 17:03
This https://math.stackexchange.com/questions/1476391/what-is-the-cross-product-in-spherical-coordinates 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 06 '17 at 23:08
I provide a formula in https://mathematica.stackexchange.com/questions/180422/simplify-closed-formula-for-distance-preserving-great-circle-spherical-coordinat but it's hideously ugly and needs to be simplified. — , Aug 22 '18 at 15:24

J. M.'s missing motivation · Answer 1 · 2017-04-22T02:02:19.023

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]]]

spherical linear interpolation

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.'s missing motivation, Apr 26 '17 at 02:39

score 1 · Answer 2 · answered Apr 07 '17 at 01:02

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}}]

Equally spaced points on a great circle arc between 2 points

2 Answers2

Linked