# Geodesics on a Sphere: An Exploration with Maxima/wxMaxima

The sphere is a comparatively simple surface but its general geodesics, aside from the lines of longitude and the equator, are more complex. 
Nevertheless, the closed equations of the general sphere geodesic can be derived although we shall not do so here.  Our purpose is to utilize
the dereived equations to explore, using Maxima/wxMaxima, some of the geometric properties of the sphere geodesic curve.

Differential geometry (DG) is a vast subject that cannot be expounded in this simple article.  The reader is assumed to have a basic knowledge
of DG as this article attempts only to demonstrate the use of Maxima/wxMaxima to eliminate the sheer drudgery of very complex DG computations.
Also, graphical plots and narrated videos are given using Geomview/StageTools.  These graphics are essential in conveying the underlying DG
concepts.  After all, geometry, differential or otherwise, has always been primarily visual.

A familiarity with the Maxima syntax is also assumed.

The goal, perhaps, of this article, and my other articles on DG, is to elevate one beyond the tedious mathematical formalism of DG and to provide
a more visceral insight.

The actual Maxima/wxMaxima script for this article is located here.  Please download to reproduce these results on your own machine.

## The parametric sphere and intrinsic properties thereof

The first step is to initialize the Maxima/wxMaxima environment:

kill(all)$ load(vect)$ load(draw)$
set_draw_defaults(head_angle=10, head_length=0.02, line_width=2, proportional_axes=xyz, nticks=200, xyplane=0, surface_hide=true,
   point_size=3, point_type=filled_circle, point_size=1)$

A standard parametric equation for the sphere is now defined and also the first and second partial derivatives.  The variable names follow the
convention stated in my previous article of the Monkey Saddle.

In this case of the sphere, the parameter \(u\) is colatitude (θ), 0-π, and \(v\) is longitude (φ), 0-2π

f: [sin(u)*cos(v), sin(u)*sin(v), cos(u)];
fu:diff(f,u);
fv:diff(f,v);
fuu:diff(fu,u);
fuv:diff(fu,v);
fvu:diff(fv,u);
fvv:diff(fv,v);

The sphere metric and inverse metric is now determined, as these are intrinsic geometric properties and will be needed later:

E:trigsimp(fu.fu)$
F:trigsimp(fu.fv)$
G:trigsimp(fv.fv)$
g:trigsimp(matrix([fu.fu, fu.fv],[fv.fu, fv.fv]));
ginv:trigsimp(invert(g));

A check to see if the matrix product gives the unit matrix:

trigsimp(g.ginv);

We now derive the unit surface normal, \( n \), of this spherical surface:

tmp:trigsimp(express(fu~fv));
n:trigsimp(tmp/sqrt(tmp.tmp));

The Christoffel symbols of the sphere are also determined.  These are also intrinsic geometric properties that will be utilized later:

load(ctensor)$
init_ctensor()$

dim:2$ ct_coords:[u,v]$
lg:g$
ug:ginv$

christof(mcs);

We need to re-express the Christoffel symbols as the Maxima result differs from what is used in standard texts on DG: 

Γ111:mcs[1,1,1];
Γ112:mcs[1,2,1];
Γ121:mcs[2,1,1];
Γ122:mcs[2,2,1];
Γ211:mcs[1,1,2];
Γ212:mcs[1,2,2];
Γ221:mcs[2,1,2];
Γ222:mcs[2,2,2];

## The general geodesics of the sphere

The geodesics on the sphere can be derived exactly.  The method, which gives the colatitude, u, and longitude, v, in terms of a parameter, t, is described in this linked file:  Cf. sphere.geodesics.pdf

The parametric equations are as follows:

\[ u = acos(\sqrt{(1-m^2)} \cdot \sin(t/R)) \]

\[ v = \atan(m \cdot \tan(t/R)) + v0 \]

In this example, we let \( m = -1/2 \), \( v0 = 0 \), and \( R = 1 \).

Also, to simplify the equations a bit, we limit \( t \in [-\pi/2, \pi/2] \) which puts a nice constraint on \( \cos(t) \).

/* We will limit t in our plots to [-π/2, π/2] to simplify the equations*/
assume(cos(t)>=0);
m:-1/2$ 
c:trigsimp(ev(f, u=acos(sqrt(1-m^2)*sin(t)), v=atan(m*tan(t))));
ct:trigsimp(diff(c,t));

Now comes a plot of the general geodesic.  Note that our limits on the parameter, \( t \), lead to a "semicircle" but that is all that we require for this demonstration:

wxdraw3d(color=gray50, parametric_surface(f[1], f[2], f[3], u,0,π,v,0,2*π),
    color=blue, parametric(c[1], c[2], c[3], t, -π/2, π/2)
)$

Check if geodesic is unit speed, i.e. if the tangent length is constant and equal to \( 1 \):

trigreduce(trigsimp(sqrt(ct.ct)));

Find the unit normal (determined above) in terms of the geodesic curve parameter \( t \):

nt:trigsimp(ev(n, u=acos(sqrt(1-m^2)*sin(t)), v=atan(m*tan(t))));

Check if unit length:

trigsimp(nt.nt);

### The local basis vectors, fu and fv, are now expressed in terms of the geodesic curve parameter t:

fut:trigsimp(ev(fu, u=acos(sqrt(1-m^2)*sin(t)), v=atan(m*tan(t))));
fvt:trigsimp(ev(fv, u=acos(sqrt(1-m^2)*sin(t)), v=atan(m*tan(t))));

This local basis (shown in green) is now plotted at several points along the geodesic curve:

wxdraw3d(color=gray50, parametric_surface(f[1], f[2], f[3], u, 0, π, v, 0, 2*π), color=red, parametric(c[1], c[2], c[3], t, -π/2, π/2),
    color=dark_green, vector(ev(c, t=0), ev(fut, t=0)),
    vector(ev(c, t=0), ev(fvt, t=0)),
    vector(ev(c, t=1/2), ev(fut, t=1/2)),
    vector(ev(c, t=1/2), ev(fvt, t=1/2)),
    vector(ev(c, t=1), ev(fut, t=1)),
    vector(ev(c, t=1), ev(fvt, t=1)),
    vector(ev(c, t=π/2-0.01), ev(fut, t=π/2-0.01)),
    vector(ev(c, t=π/2-0.01), ev(fvt, t=π/2-0.01))
)$

### Express the tangent to the geodesic curve in terms of the above basis vectors

Find the tangent to geodesic curve \( c \) in terms of un-normalized basis vectors via vector projection,
i.e. \( ct = a \cdot fut + b \cdot fvt \)

The derivatives of the coeffients, \( a, b \), are also determined for later use.

a:trigsimp((ct.fut)/(fut.fut));
b:trigsimp((ct.fvt)/(fvt.fvt));
at:trigsimp(diff(a,t));
bt:trigsimp(diff(b,t));

Check if this basis expression equals the geodesic curve tangent (derivative) that was determined above:

trigsimp(trigreduce(trigsimp((a*fut+b*fvt))));
ct;

OK.  They are equal.

## The Covariant Derivative of the Geodesic Tangent Vector and General Vectors

### The covariant derivative of the geodesic tangent

We now come to the heart of the exploration.  The covariant derivative of the geodesic tangent should be \( 0 \) everywhere because geodesics
experience no curvature within the surface (i..e. geodesic curvature).  The geodesic is the generalization of the Euclidian straight line. 

First, we need to express the \( u,v \) curve parameters in terms of the geodesic curve parameter \( t \):

ut:trigsimp(diff(acos(sqrt(1-m^2)*sin(t)),t));
vt:trigsimp(diff(atan(m*tan(t)),t));

We will now calculate the covariant derivative of the geodesic tangent using the Christoffel symbols:

tmp:(trigsimp((at+Γ111*a*ut + Γ112*a*vt + Γ121*b*ut + Γ122*b*vt)*fut + (bt+Γ211*a*ut + Γ212*a*vt + Γ221*b*ut + Γ222*b*vt)*fvt)); 

covd:trigsimp(ev(tmp, u=acos(sqrt(1-m^2)*sin(t)), v=atan(m*tan(t))));

Indeed!  The covariant derivative of the curve tangent is the zero vector everywhere on the geodesic.  Once again, the power of Maxima/wxMaxima
to deal with very complex expessions is quite evident.

The covariant derivative can also be determined by subtracing the projection of the directional derivative onto the unit normal from the
directional derivate itself.  This method, in contrast to the above, is an example of extrinsic geometry, but it serves as a useful check:

dirdiv:trigsimp(diff(ct, t));

covdn:trigsimp(dirdiv-(dirdiv.nt)*nt);

As expected, the result is the same, but the intrinsic method involving Christoffel symbols is preferred.

Now, both the geodesic basis vectors (green) and the geodesic tangent vector (red) are plotted.

wxdraw3d(color=gray50, parametric_surface(f[1], f[2], f[3], u, 0, π, v, 0, 2*π), color=blue, parametric(c[1], c[2], c[3], t, -π/2, π/2),
    color=dark_green, vector(ev(c, t=1/2), ev(fut, t=1/2)),
    vector(ev(c, t=1/2), ev(fvt, t=1/2)),
    color=red, vector(ev(c, t=1/2), ev(ct, t=1/2))
)$

### The covariant derivative of a general vector along the geodesic curve

The geodesic curve tangent is only one of an infinity of vectors in the tangent plane at a given curve point.  We now derive the covariant derivative
of a general geodesic curve vector: \( a \cdot fut + b \cdot fvt \).  To make things simple, we use \( a=1, b=1 \) as the vector coefficients.  The coefficient
derivatives are also determined:

a:1;
b:1;
at:trigsimp(diff(a,t));
bt:trigsimp(diff(b,t));

The vector, \( V \), is now expressed:

V:trigsimp(a*fut+b*fvt);

Using the intrinsic method, with Christoffel symbols, the covariant derivative of vector \( V \) along the geodesic curve is now calculated:

tmp:(trigsimp((at+Γ111*a*ut + Γ112*a*vt + Γ121*b*ut + Γ122*b*vt)*fut + (bt+Γ211*a*ut + Γ212*a*vt + Γ221*b*ut + Γ222*b*vt)*fvt));

covd:trigsimp(ev(tmp, u=acos(sqrt(1-m^2)*sin(t)), v=atan(m*tan(t))));

As expected, the covariant derivative of this general tangent vector is not eveywhere equal to \( 0 \).

We now plot the general vector, \( V \), in blue, and it's covariant derivative (in black) at several points along the geodesic curve: 

wxdraw3d(color=gray50, parametric_surface(f[1], f[2], f[3], u,0,π,v,0,2*π), color=red, parametric(c[1], c[2], c[3], t, -π/2, π/2),
    color=blue, vector(ev(c, t=1), ev(V, t=1)),
    color=black, vector(ev(c, t=1), ev(covd, t=1)),
    color=blue, vector(ev(c, t=1/2), ev(V, t=1/2)),
    color=black, vector(ev(c, t=1/2), ev(covd, t=1/2)),
    color=blue, vector(ev(c, t=0), ev(V, t=0)),
    color=black, vector(ev(c, t=0), ev(covd, t=0)),
    color=blue, vector(ev(c, t=-1/2), ev(V, t=-1/2)),
    color=black, vector(ev(c, t=-1/2), ev(covd, t=-1/2))
)$

As a final check, both the general vector, \( V \), and the covariant derivative should be in the tangent plane, i.e. both should be
perpendiclar to the unit normal at the same point (their dot products should be \( 0 \)).

Is covd perpendicular to the unit normal?

trigsimp(sqrt(covd.nt));

Is \( V \) also in tangent plane (i.e. perpendicular to the normal)?

trigsimp(sqrt(V.nt));

## Final Comments

Again, this article hopefully demonstrates how the use of the Maxima/WxMaxima computer algebra software can eliminate the tediousness of DG and allow a visual insight to be easily developed.

Fortunately, exact geodesic equations are possible for the sphere and this allows geodesic properties to be readily explored.  For more complex surfaces, however, this will not necessarily be the case.

dirdiv2:trigsimp(diff(V, t));

covdn2:trigsimp(dirdiv2-(dirdiv2.nt)*nt);

covd;