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.

  • $\begingroup$ "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. $\endgroup$ Commented Apr 6, 2017 at 16:59
  • $\begingroup$ 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. $\endgroup$
    – Morgan
    Commented Apr 6, 2017 at 17:03
  • $\begingroup$ 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. $\endgroup$
    – Bill
    Commented Apr 6, 2017 at 23:08
  • $\begingroup$ I provide a formula in mathematica.stackexchange.com/questions/180422/… but it's hideously ugly and needs to be simplified. $\endgroup$
    – user1722
    Commented Aug 22, 2018 at 15:24

2 Answers 2


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
  • $\begingroup$ Good, only in my question I have theta, then phi, so your #1 and #2 are the other way around. $\endgroup$
    – Morgan
    Commented Apr 25, 2017 at 19:17
  • $\begingroup$ @Morgan, I'll leave the conversion to your preferred convention up to you; I'm accustomed to putting the longitude before the colatitude. $\endgroup$ Commented Apr 26, 2017 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, 

For example:

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

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.