# 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!

• Possible duplicate of A time parameterization of geodesics on the sphere Commented Jul 6, 2018 at 22:32
• I’m voting to close this question because OP stopped interacting with the users trying to help. Commented Jul 24, 2023 at 11:10
• @Lucas: So radius of sphere $a=r$ ? Commented Apr 8 at 3:53

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

• Am I mistaken, or is your equation a formula for the plane $ax+by+cz=0$? If so I think that is not what the OP wanted. Commented Jul 8, 2018 at 6:38
• Yes, or you can think of it as a dot product g.p = 0, where g is on the geodesic and p is the pole of the geodesic. The OP said he got nowhere with this approach, but I don't think he definitely wanted something else. Commented Jul 8, 2018 at 7:36

You want a relation between longitude ( in this case spherical polar angle $$\theta$$) the latitude ( in this case co-latitude $$\phi$$), and inclination of great circle $$\alpha$$ to the plane of the equator.

We can find this from standard trigonometrical relations:

Ellipse in polar coordinates $$\frac{x^{2}}{a^{2}} +\frac{y^{2}}{b^{2}} =1, (x,y)=r(\cos \theta, \sin \theta), \frac{\cos^{2} \theta}{a^{2}} +\frac{\sin^{2}\theta}{b^{2}} = \frac{1}{r^{2}} \tag 1$$ Plug in inclination and co-latitude $$\sin \alpha=\frac{b}{a},~\cos \phi=\frac{r}{a} \tag2$$ Simplify to find the required geodesic /great circle equation

$${\cos^{2} \theta}+\frac{\sin^{2}\theta}{\sin^{2}\alpha} =\sec^{2}\phi \tag3$$

which further simplifies to:

$$\sin \theta= \tan \alpha \tan \phi \tag4$$

For cardinal points $$(A,B)$$:

$$\theta=0,\phi==0;~~\theta=\pi/2,\sin \alpha=\cos \phi,~\alpha+\phi_{max}=\pi/2 ; \tag 5$$

When the given ellipse polar projection makes initial angle by rotation.. with principal axis $$\theta_{0},$$

$$\sin (\theta-\theta_{0}) = \tan \alpha \tan \phi \tag6$$

The sphere and the geodesic (great circle in red) can be 3D pictured in a CAS. For the particular case shown the boundary conditions are:

$$r_i=a,~\theta_0=0,~\alpha=30^{\circ}.$$