Define a function to plot a geodesic arc between 2 points on a unit sphere

Questions
,
1

I need to define a function to plot a geodesic arc between two points on a unit sphere. Can someone help me with this? Thanks.

2

How about GeoGraphics:

GeoGraphics[GeoPath[{{30, -20}, {33, 35}}], GeoBackground -> None, 
 GeoRange -> "World", GeoProjection -> "Orthographic", 
 GeoGridLines -> Automatic]
3

I had seen this example in the geodesic documentation and decided not to use it because the plots I need have to be in 3D. I tried using Graphics3D but couldn’t make it work with GeoPath. Sorry I wasn’t more explicit in my post.

4

Some years ago I made this utility for precisely the same goal:

geodesicArc[v1_List, v2_List, center_List, n_Integer] :=
  
  Module[{\[Alpha]0 = VectorAngle[v1 - center, v2 - center], \[Beta], 
    f},
   f[\[Beta]_] = (
    EuclideanDistance[v1, center] Sec[
      1/2 (\[Alpha]0 - 2 \[Beta])] Sin[\[Beta]])/
    EuclideanDistance[v1, v2];
   Line[Table[
     center + 
      Norm[v1 - center]*
       Normalize[t*(v1 - center) + (1 - t) (v2 - center)], {t, 
      Map[f, Subdivide[\[Alpha]0, n]]}]]];
Graphics3D[{Sphere[], 
  geodesicArc[FromSphericalCoordinates[{1, .5, .1}], 
   FromSphericalCoordinates[{1, 2, -.1}], {0, 0, 0}, 15]}]

The parameter n is the number of segments in the arc.

5

Thank you for this. Would you have an example of calling this function? Thanks.

6

Yes, at the bottom of the code you can find an example.

7

Yes, thank you! The routine works perfectly.

8

That is an interesting question.
I pursued another more geometric approach. The idea is to turn the sphere in a manner that both points lie in the x,y plane, define a vector circling in this plane and turn it back :

Turning matrices

dx[a_] := ( {
   {1, 0, 0},
   {0, Cos[a], -Sin[a]},
   {0, Sin[a], Cos[a]}
  } )
dy[a_] := ( {
   {Cos[a], 0, -Sin[a]},
   {0, 1, 0},
   {Sin[a], 0, Cos[a]}
  } )
dz[a_] := ( {
   {Cos[a], -Sin[a], 0},
   {Sin[a], Cos[a], 0},
   {0, 0, 1}
  } )

factor to convert degrees into arc-lengths

f = Pi/180.;

Coordinates of the two points, converted to polar coordintates

c1 = {30, -20};
c2 = {33, 35};

cc1 = Function[{x, y}, f {x, 90 - y}] @@ c1;
cc2 = Function[{x, y}, f {x, 90 - y}] @@ c2;

Function to convert polar coordinates to a vector in R3

xx[{a_, b_}] := {Cos[a] Sin[b], Sin[a] Sin[b], Cos[b]}

Now do the turning job

x1 = xx[cc1];
x2 = xx[cc2];
x1a = dy[cc1[[2]] - Pi/2].dz[-cc1[[1]]].x1;
x2a = dy[cc1[[2]] - Pi/2].dz[-cc1[[1]]].x2;
alpha = ArcTan[x2a[[2]], x2a[[3]]];
x1b = dx[-alpha].x1a // Chop;
x2b = dx[-alpha].x2a // Chop;

And define the geodesic

geo[a_] := dz[cc1[[1]]].dy[Pi/2 - cc1[[2]]].dx[alpha].dz[a].x1b

Show the whole thing

pl1 = Show[
  ParametricPlot3D[
   geo[x],
   {x, 0, Pi},
   PlotStyle -> {Blue, Thick}
   ],
  Graphics3D[{Red, PointSize[.02], Point /@ {x1, x2}}],
  Graphics3D[{Opacity[.3], Sphere[]}]
  ]

Here is the function describing the geodesic:

geo[t]

It should be easy to comprise the steps into a function returning the expression of the geodesic.

Unfortunately I cannot compare geo[t] with the results of Gianluca, because his code does not run on my machine (too old version of Mma)

9

Yet another approach. Using RotationMatrix and ParametricPlot3D

greatCircleArc[{\[Theta]1_, \[Phi]1_}, {\[Theta]2_, \[Phi]2_}] := 
 Module[
  	{v1, v2, normal, \[Alpha]},
  	v1 = FromSphericalCoordinates@{1, \[Theta]1, \[Phi]1};
  	v2 = FromSphericalCoordinates@{1, \[Theta]2, \[Phi]2};
  	normal = Cross[v1, v2];
  	\[Alpha] = VectorAngle[v1, v2];
  	ParametricPlot3D[RotationMatrix[t, normal].v1, {t, 0, \[Alpha]}, PlotStyle -> Black]
 ]

Show[
 	Graphics3D@Sphere[],
 	greatCircleArc[{0.5, 0.1}, {2, -0.1}]
]