surfaces of revolution [Maple worksheet]

Like the sphere, surfaces of revolution have a grid formed from the vertical meridians (like lines of longitude) and the horizontal parallels (circles, like circles of latitude).
These always intersect at right angles, so they form an orthogonal grid on the surface.

When the horizontal radius of a parallel circle from the vertical axis of symmetry is a local extremum (blue circles), this parallel is an obvious geodesic or autoparallel curve. If you imagine a tiny car on the surface with its steering wheel locked to move straight ahead, it will remain on such a curve. [Click here for a terrific video of such a trip on a two holed torus---skip down to the motorcycle on the blacktop on the orange surface below the two holed torus video.]
The same is true of the meridians. The nontrivial geodesics follow neither a meridian or extremal parallel.

We parametrize such a surface by an arclength radial coordinate r along the meridians and an azimuthal angle φ around the vertical z axis along the parallels measured from the x axis in the usual counterclockwise sense in the z = 0 plane.

The ring torus is obtained by rotating a circle perpendicular to its plane about an axis outside the circle. Generalizing the sphere terminology, it has an inner and outer equator, a North Polar Circle and a South Polar Circle. Apart from the inner equator itself, geodesics are divided into 2 classes, those which pass through the inner equator, and those which do not. Like the great circles on the sphere, it also has closed geodesics which retrace their own path once they return to it, but they are a discrete spectrum.

The right image is a unit ring torus with  (a,b) = (2,1), so the origin of  coordinates (r,φ) = (0,0) is at (x,y,z) = (3,0,0).
Note red outer equator, green inner equator, blue prime meridian, while the North Polar Circle and the South Polar Circle are not indicated.

If we consider the problem of free motion constrained to the torus surface, we replace the centrifugal potential of the Kepler problem using R(r) instead of r in the formula for the kinetic energy: U(r)  = Ɩ2/R(r)2, while removing the gravitational potential. [See below.]

For any surface of revolution we can easily evaluate the conserved energy and angular momentum, noting that the arclength differential obeys the infinitesimal Pythagorean theorem to combine the orthogonal differentials dr and R(r) dφ just like in ordinary polar coordinates in the plane where R(r) = r. We can set the mass of the moving particle to m = 1, since it is irrelevant here (so we are using the energy and angular momentum per unit mass below).

The centrifugal potential term left in the kinetic energy after using conservation of angular momentum to describe the radial motion alone acts like the potential for the fictitious centrifugal force due to the rotational motion around the axis of symmetry.

If we choose unit speed for an arclength parametrization, this corresponds to energy E = 1/2, but then the centrifugal potential depends on the angular momentum.
However, it is more useful to choose unit angular momentum l = 1 for nonradial geodesics, and let the energy levels parametrize the radial motion in a fixed centrifugal potential that does not change.

Here is a graph of the potential for the unit ring torus, setting the angular momentum parameter to Ɩ = 1.

one period of periodic potential, sufficient for bound geos (E ≤ 1/2)

many periods of periodic potential needed for unbound geos (E > 1/2)
shown are the energy levels for those geos making
exactly m windings of the torus during 1 azimuthal revolution

Bound orbits (like green energy levels) have energy less than 1/2 and do not pass through the donut hole, or exactly energy 1/2 for the inner equator, with the outer equator at the minimum of the potential.

Unbound orbits (blue energy levels) pass through the donut hole an infinite number of times with a periodic potential needed to describe their unbounded r values.

The closed geodesics where the geo returns to its starting point after m radial oscillations and n azimuthal revolutions around the symmetry axis are interesting.
Let [m,n; p] characterize a closed geo with  m radial oscillations and n azimuthal revolutions around the symmetry axis, p =1 for unbound, p = 0 for bound.
Here is the (m,n; p) = [7,1;1] case shown above.

The  [7,1;1] closed geo.
On the right: red, blue, green curves to show x, y, z as functions of the arclength t.
The end of the first blue cycle shows the first azimuthal revolution and gives the period T,
which is the total arclength in the arclength parametrization used in the numerical solver,
with initial data posed at the origin (r,φ) = (0,0) at (x,y,z) = (3,0,0).

This is trial and error at which I stopped short of higher precision,
so we are not exactly at the origin of coordinates at (3,0,0).
Its initial angle from the vertical is about 6.82.

To use the torus-demo Maple worksheet, one only has to choose an initial angle of a geodesic with the vertical direction at the origin of coordinates on the outer equator at the prime meridian. The numerical differential equation solver then finds the path parametrized by the arclength t (like time) for a given interval of that parameter. It then gives you the above two plots and calculates the endpoint position and total length. Solving for one revolution (the end of the first cycle of the blue curve), one sees where one first recrosses the prime meridian. One can thus decrease the angle to azimuthally contract the geodesic or increase the angle to azimuthally stretch the geodesic to try to get it to move closer and closer to the starting point, by intelligent trial and error of course.

This shows the initial data plane with the geodesic for angle 15 degrees. This makes about 2.5 oscillations about the outer equator before returning to the prime meridian just below the inner equator.

Play this game. Click here if you have Maple.

There are many interesting things one can investigate about geodesics. One of these is the focusing of geodesics by positive curvature which occurs in the outer region between the polar circles. The largest curvature affects are along the outer equator, where focusing of the geodesics causes them to fall back to the equator and cross, leading to shorter paths connecting points on the outer equator after that which take a short cut by moving away from that equator. This effect is responsible for gravitational lensing in general relativity.

The spray of geodesics emanating from the origin of coordinates on the left cross on the outer equator at the right. One can easily calculate this minimum crossing distance.

Kepler problem

The nonrelativistic and relativistic Kepler problems have instead the radial motion energy equation with extra potentials for the gravitational force.

This follows easily from Gravitation (Misner, Thorne and Wheeler: Wiki, Amazon) chapter 25 for timelike geodesics in the Schwarzschild spacetime, p.656. The details are discussed in the dr bob differential geometry textbook, section 8.3.