Finding the equation of a geodesic passing through two given points I am trying to find the equation of a geodesic (otherwise known as a great circle or great circle arc) on the surface of a sphere of given radius $a$ through two points on the sphere. I am given the coordinats of the two points in spherical coordinates, $(r,\theta_1,\phi_1)$ and $(r,\theta_2,\phi_2)$. The question is, how do I do this?
Upon googling, the only remotely useful result I found was from a Wolfram Mathworld site, where they derived the following equation:
$$ a\cos (u) \sin (v) \sin (c_2) + a\sin(u)\sin(v)\cos(c_2)-\frac{a\cos(v)}{\sqrt{(a/c_1)^2-1}} = 0  $$
where $u:=\theta$ and $v:=\frac{1}{2}\pi-\phi$, and $c_1$, $c_2$ are arbitrary constants (of integration?!). Then I was able to derive quite easily
$$ \cos(\theta)\cos(\phi)\sin(c_2) + \sin(\phi)\cos(\phi)\cos(c_2) -\frac{\sin(\phi)}{\sqrt{(a/c_1)^2-1}} = 0 $$
but this gets me nowhere. Firstly, I do not understand its derivation at all (the part regarding the partial derivates $P$, $Q$ and $R$ defeats me). Secondly, I have no idea why the constants $c_1$, $c_2$ come into play and thirdly, I do not know how even to use this formula to derive anything I want since the coordinates of the given points are not there in the first place.
I have tried to derive a formula on my own by using vectors, trying to describe the geodesic as the intersection of a plane $a_1x+a_2y+a_3z=0$ and the sphere, etc. but I have got nowhere. Any help is appreciated, whether it is in helping to explain Wolfram's formula or deriving a new one. Thanks in advance!
 A: Not sure what your $\theta, \phi$ are. I use $\phi$ for latitude, $\lambda$ for longitude. With the usual convention about axes,
$$x = \cos\phi\cos\lambda,\quad y=\cos\phi\sin\lambda,\quad z = \sin\phi.$$
Your idea $ax + by + cz = 0$ for a geodesic looks OK. Since the geodesic goes through $(\phi_1,\lambda_1)$ and $(\phi_2,\lambda_2)$ we get (up to multiplication by a constant)
$$\eqalign{a &=  \cos\phi_1\sin\phi_2\sin\lambda_1 - \cos\phi_2\sin\phi_1\sin\lambda_2,\cr
  b &= -\cos\phi_1\sin\phi_2\cos\lambda_1 + \cos\phi_2\sin\phi_1\cos\lambda_2,\cr
  c &= \cos\phi_1\cos\phi_2\sin(\lambda_2 - \lambda_1)}.$$
The equation $ax + by + cz = 0$ in terms of $(\phi,\lambda)$ is then
$$[\cos\phi_1\sin\phi_2\sin(\lambda - \lambda_1)
 -\cos\phi_2\sin\phi_1\sin(\lambda - \lambda_2)]\cos\phi
 -\cos\phi_1\cos\phi_2\sin(\lambda_2 - \lambda_1)\sin\phi = 0,$$
which can be re-arranged as
$$\tan\phi = {{\cos\phi_1\sin\phi_2\sin(\lambda - \lambda_1) - \cos\phi_2\sin\phi_1\sin(\lambda - \lambda_2)}\over{\cos\phi_1\cos\phi_2\sin(\lambda_2 - \lambda_1)}}.$$
The denominator is non-zero provided the geodesic doesn't pass though the poles. The last equation can be expressed as
$$\tan\phi = k\cos(\lambda - \lambda_0)$$
where $k$ and $\lambda_0$ are constants. This I think is equivalent to the formula given at:
http://www.damtp.cam.ac.uk/user/reh10/lectures/nst-mmii-handout2.pdf
