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)