Practical drawing of geodesics I want to use a computer to draw geodesics on a known parameterized surface of revolution, starting from a known point and at a known angle to the meridian.
What would be the easiest way of doing this? 
Some methods I've considered:


*

*Using Clairaut's relation ($r \cos \theta = C$), and using
small increments, but I'm afraid that this is not a general solution, as this will cause the geodesic to
"get stuck" on parallels whenever the angle crosses zero.

*Use the relation that for a surface given by:
$$\left(\varphi(v)\cos u, \varphi(v)\sin u, \psi(v)\right)$$ The geodesics are given by: $$\frac{du}{dv}=\frac{\sqrt{\varphi'^2+\psi'^2}}{\varphi\sqrt{\varphi^2-c^2}}\ \longrightarrow\ u(v_1) = u_0 + \int_{v_0}^{v_1} \frac{\sqrt{\varphi'^2+\psi'^2}}{\varphi\sqrt{\varphi^2-c^2}} dv$$ This presents the problem of finding the correct value of $c$ given the starting angle, and as user levap pointed out, has the same "getting stuck" problem, since $u(v)$ is no longer a function.

*Various numerical methods which work for any triangulated surface, but I feel it would be a shame to use these when I know the exact form of the surface.


Are there any easier methods or modifications to the above?
 A: Assume that we are given a curve $\gamma \colon J \rightarrow \mathbb{R}^3$ of the form $\gamma(v) = (\varphi(v), 0, \psi(v))$ such that $\frac{d \gamma}{dv}(v) \neq 0$ and $\varphi(v) > 0$ for all $v \in J$. The surface of revolution obtained by revolving $\gamma$ around the $z$-axis is 
$$ S = \{ (\varphi(v) \cos u, \varphi(v) \sin u, \psi(v)) \, | \, v \in J, u \in [0,2\pi] \}. $$
Since $S$ cannot be covered with a single coordinate chart, it is more useful to consider $S$ as the image of the map $X \colon \mathbb{R} \times J \rightarrow \mathbb{R}^3$ given by
$$ X(u,v) = (\varphi(v) \cos(u), \varphi(v) \sin u, \psi (v)). $$ 
The map $X$ is a local diffeomorphism and the pullback of the first fundamental form from $S$ to $\mathbb{R} \times J$ is given by
$$ g(u,v) = \begin{pmatrix} \varphi(v)^2 & 0 \\ 0 & \varphi'(v)^2 + \psi'(v)^2 \end{pmatrix} = \begin{pmatrix} E(v) & 0 \\ 0 & F(v) \end{pmatrix}. $$
Geodesics on $\mathbb{R} \times J$ with the pullback metric will map under $X$ to geodesics of $S$. Using $g$, we can calculate the Christoffel symbols and write explicitly the geodesic equation $\nabla_{\dot{\alpha}(t)}{\dot{\alpha}} = 0$ for a curve $\alpha(t) = (u(t), v(t))$ in $\mathbb{R} \times J$. Written shortly, the equations are
$$ \ddot{u} + \frac{E'}{E} \dot{u} \dot{v} = 0, \\
\ddot{v} - \frac{E'}{2G} \dot{u}^2 + \frac{G'}{2G} \dot{v}^2 = 0. $$
More explicitly, we have
$$ \ddot{u}(t) + \frac{2 \varphi(v(t)) \varphi'(v(t))}{\varphi(v(t))^2} \dot{u}(t) \dot{v}(t) = 0, \\
\ddot{v}(t) - \frac{\varphi(v(t)) \varphi'(v(t))}{\varphi'(v(t))^2 + \psi'(v(t))^2} \dot{u}(t)^2 + \frac{\varphi'(v(t)) \varphi''(v(t)) + \psi'(v(t)) \psi''(v(t))}{\varphi'(v(t))^2 + \psi'(v(t))^2} \dot{v}(t)^2 = 0. \label{equation1}\tag{1}$$
This is a system of (coupled) non-linear second order differential equations for $(u(t), v(t))$ which can be solved numerically given initial conditions at some (arbitrary) time $t = t_0$ of the form
$$ \alpha(t_0) = (u_0, v_0), \,\,\, \dot{\alpha}(t_0) = (u_0', v_0') $$
and you can use them to find and draw the geodesic $X(\alpha(t))$. The geodesic will be of constant speed (with respect to the metric $g$) and won't be generally a graph (so $\alpha(t)$ won't be of the form $\alpha(t) = (t, v(t))$ or $\alpha(t) = (u(t), t)$).
The description above doesn't really use the fact that we are looking for geodesics on a surface of revolution that has a special symmetry. The special symmetry gives a conserved quality that must be constant along geodesics. This is called Clairaut's relation:
$$ \varphi(v(t)) \cos \theta(t) = C $$ 
where $\theta(t)$ is the angle the geodesic $\alpha(t)$ makes with the parallel $v = v(t)$ and $C$ is some constant. In fact, the converse almost holds in the sense that if $\alpha$ is a regular curve that doesn't coincide with a parallel on any interval and satisfies Clairaut's relation then $\alpha$ is a reparametrization of a geodesic.
If we assume that $\alpha$ has the form of a graph $\alpha(u) = (u, v(u))$ (note that this is different than the equation you present in which $\alpha(v) = (v(u), u)$) and plug $\alpha$ into Clairaut's relation, we obtain the following equation for $v$:
$$  \frac{ dv}{du} = \pm \frac{1}{C} \frac{\varphi(v(u)) \sqrt{\varphi^2(v(u)) - C^2}}{\sqrt{ \varphi'(v(u))^2 + \psi'(v(u))^2}} = F(v). \label{eq2}\tag{2} $$ 
Thus, instead of having two second order equations, we have reduced the problem to a single first order equation for $v$. Again, given initial conditions $v(u_0) = v_0$ with $\varphi(v_0) \neq \pm C$, the equation has a unique solution which can found numerically. However, this equation has the drawback that it cannot really used to find geodesics which, at $(u_0, v_0)$, are tangent to a parallel. The reason is that the problem is not well-posted - if $\alpha$ is tangent to a parallel at $(u_0, v_0)$ then $\cos (\theta) = 1$ and hence, $\varphi(v_0) = C$ and $F(v_0) = 0$. The right hand side is not a Lipschitz function and you don't have uniqueness of solutions. Indeed, the solution $v(u) \equiv v_0$ is a solution of ($\ref{eq2}$) but it is not necessarily a geodesic - it will be a geodesic if and only if $\varphi'(v_0) = 0$. If it is not a geodesic, then equation ($\ref{eq2}$) will have two solutions.
The advantage of seeking a parametrization of the form $(u, v(u))$ lies in the fact that any geodesic which is not a meridian can be parametrized in this way. A parametrization of the form $(u(v), v)$ is possible only for geodesics that doesn't become tangent to a parallel at some point. If we write the equation for $u(v)$ (which is the inverse function of $v(u)$), we get (more or less, up to a constant that seems to have been lost) the equation you bring in your second point. In some sense, this improves the situation because the equation you get for $u(v)$ is just a direct integral which can be evaluated numerically. However, it suffers from a similar problem as equation ($\ref{eq2}$) because if the geodesic becomes parallel to a tangent at some point, the integrand blows up and the derivative $\frac{du}{dv}$ approaches $\pm \infty$. This is not surprising, because one can't expect such a parametrization to hold near a point of tangency to a parallel.
Taking all this into account, I can offer a few different methods for drawing geodesics:


*

*Just use equations ($\ref{equation1}$). They look more complicated but unless you will have performance issues, solving them numerically will get you want without worrying about whether the geodesic becomes tangent to a parallel or not.

*Use equation ($\ref{eq2}$) and the symmetry involved. More explicitly, let us assume you want to find a geodesic that is not a meridian, emanating forward from $(u_0, v_0)$, whose angle with the parallel $v = v_0$ is $\cos \theta_0 \neq 1$. Set $C = \varphi(v_0) \cos \theta_0$ and solve equation ($\ref{eq2}$) numerically to get (some approximation of) $v(u)$. There are three possible cases:


*

*The geodesic won't become tangent to a parallel at any point.

*The geodesic won't become tangent to a parallel but will become asymptotic to a parallel. This means that the solution will satisfy $\lim_{u \to \infty} v(u) = v_1$ and $F(v(u)) \neq 0$ for all $u > u_0$. This is possible if and only if the parallel itself is a geodesic, which can be detected by checking whether $\varphi'(v_1) = 0$.

*The geodesic will become tangent to a parallel after some finite number of windings around the surface. This means that the solution will satisfy $\lim_{u \to u_1} v(u) = v_1$ and $F(v_1) = 0$ (so $v'(u_1) = 0$). To continue the solution beyond $u_1$, use the fact that the geodesic bounces back from the parallel tracing its reflection as Andrew D. Hwang notes. More formally, the solution will satisfy $v(u) = v(2u_1 - u)$ for $u_1 < u < u_1 + (u_1 - u_0)$. Then repeat the process.


*The method described in the previous item does not allow you to find a geodesic that starts as tangent to a parallel. In order to do that, you can combine equations ($\ref{equation1}$) and ($\ref{eq2}$). Namely, assume that you want to find a geodesic $\alpha(t) = (u(t), v(t))$ satisfying $\alpha(t_0) = (u_0, v_0)$ and $\dot{\alpha}(t_0) = (1,0)$. Linearizing equations ($\ref{equation1}$), we have
$$ u(t) \approx u_0 + (t - t_0), \\
v(t) \approx v_0 + \frac{1}{2}\frac{\varphi(v_0) \varphi'(v_0)}{\varphi'(v_0)^2 + \psi'(v_0)^2} (t - t_0)^2 $$
Use this to find $\alpha(t_1)$ for some small $t_1 > t_0$ and then use the method of the previous item to continue solving for the geodesic.
A: Triangulation need not any more be adopted once you know Clairaut's Law as operative and geodesic trajectory is found by quadrature for the surfaces of revolution.
EDIT2:
The parametrization you gave is not full, valid for all meridians.
Wlog you could consider 
$$ X(u,v) = (v \cos u, v \sin u, \psi (v) ) $$
$C$ is determined by initial condition:  $ C_{initial}=  v_0 \sin \theta_{0}. $
It is beneficial with respect to easily crossing  the parallels in numerical computations to work on arc length basis (primes) using Liouville's theorem for zero geodesic curvature lines aka  geodesics.
$$ \theta^{'}(s) = -\sin \theta \sin \phi / v $$
where the meridian has slope: 
$$ \dfrac{dv}{dz} = \tan \phi $$
and polar coordinate rotation:
$$ v\,u^{'} = \sin \theta $$ 
EDIT1:
Say you want to practically compute and construct a great circle on a sphere. Using Runge-Kutta_4 numerically integrate the set of equations to define the circle in space and supply boundary conditions.
$$ \theta ^{'} = -\sin \theta \sin \phi / v ;\;  \phi^{'} = 1/ (a \cos \phi ) ;\;   v^{'} = \sin \phi \cos \theta;\;  z^{'} = \cos \phi \cos \theta\; ; u{'}=  \sin \theta / v ;  $$
(The following is more advanced, after some experience gained in geodesic tracing:)
Geodesics appear to cross all parallels before the one corresponding to
$$ r_{min} = C, \theta = \pi/2 $$
There are however, three ways geodesics move during turn-around at above extreme parallel: 
The returning behavior types are commonly traced for:


*

*Positive/Negative Gauss curvatures: Geodesics are returning after meeting parallel circle $ r_{min}< C_{initial} $ at turn-around.


Non-returning types are for Negative Gauss curvature:


*Geodesics shoot through after meeting $ r_{min}= C$ if $ r_{min} > C_{initial} $

*Geodesics are asymptotic (never reach)  $ r_{min}= C$ if $ r_{min} = C_{initial} $
The latter are demonstrated in WolframAlpha Demo:
ONE_sheet_Hyperboloid_geodesics
