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$ – J. M.'s discontentment Apr 6 '17 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 Apr 6 '17 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 Apr 6 '17 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 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]]]

spherical linear interpolation

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

{ArcTan @@ Most[#], ArcCos[Last[#]/Norm[#]]} & /@ arcDots
| improve this answer | |
  • $\begingroup$ Good, only in my question I have theta, then phi, so your #1 and #2 are the other way around. $\endgroup$ – Morgan Apr 25 '17 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$ – J. M.'s discontentment 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, 

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

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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