The great circle path may be found using
spherical trigonometry; this is the spherical version of the
inverse geodetic problem. If a navigator begins at
P1 = (φ1,λ1) and plans to travel the great circle to a point at point
P2 = (φ2,λ2) (see Fig. 1, φ is the latitude, positive northward, and λ is the longitude, positive eastward), the initial and final courses α1 and α2 are given by
formulas for solving a spherical triangle :\begin{align} \tan\alpha_1&=\frac{\cos\phi_2\sin\lambda_{12}}{ \cos\phi_1\sin\phi_2-\sin\phi_1\cos\phi_2\cos\lambda_{12}},\\ \tan\alpha_2&=\frac{\cos\phi_1\sin\lambda_{12}}{-\cos\phi_2\sin\phi_1+\sin\phi_2\cos\phi_1\cos\lambda_{12}},\\ \end{align} where λ12 = λ2 − λ1 and the quadrants of α1,α2 are determined by the signs of the numerator and denominator in the tangent formulas (e.g., using the
atan2 function). The
central angle between the two points, σ12, is given by :\tan\sigma_{12}=\frac{\sqrt{(\cos\phi_1\sin\phi_2-\sin\phi_1\cos\phi_2\cos\lambda_{12})^2 + (\cos\phi_2\sin\lambda_{12})^2}}{\sin\phi_1\sin\phi_2+\cos\phi_1\cos\phi_2\cos\lambda_{12}}.{{refn|group=note|A simpler formula is : \cos\sigma_{12}=\sin\phi_1\sin\phi_2+\cos\phi_1\cos\phi_2\cos\lambda_{12}; however, this is numerically less accurate if σ12 small.}}{{refn|group=note|These equations for α1,α2,σ12 are suitable for implementation on modern calculators and computers. For hand computations with logarithms,
Delambre's analogies were usually used: : \begin{align} \cos\tfrac12(\alpha_2+\alpha_1) \sin\tfrac12\sigma_{12} &= \sin\tfrac12(\phi_2-\phi_1) \cos\tfrac12\lambda_{12},\\ \sin\tfrac12(\alpha_2+\alpha_1) \sin\tfrac12\sigma_{12} &= \cos\tfrac12(\phi_2+\phi_1) \sin\tfrac12\lambda_{12},\\ \cos\tfrac12(\alpha_2-\alpha_1) \cos\tfrac12\sigma_{12} &= \cos\tfrac12(\phi_2-\phi_1) \cos\tfrac12\lambda_{12},\\ \sin\tfrac12(\alpha_2-\alpha_1) \cos\tfrac12\sigma_{12} &= \sin\tfrac12(\phi_2+\phi_1) \sin\tfrac12\lambda_{12}. \end{align} McCaw yields for the angle between the great circles through that point to the North on one hand and to on the other hand ::\cos\theta_{t,N} = \cos\theta_{s,t}\cos\theta_{s,N}+\sin\theta_{s,t}\sin\theta_{s,N}\cos p. ::\sin\varphi_t = \cos\theta_{s,t}\sin\varphi_s +\sin\theta_{s,t}\cos\varphi_s\cos p. The
sine formula yields ::\frac{\sin p}{\sin \theta_{t,N}} = \frac{\sin(\lambda_t-\lambda_s)}{\sin\theta_{s,t}}. Solving this for and insertion in the previous formula gives an expression for the tangent of the
position angle, ::\sin\varphi_t = \cos\theta_{s,t}\sin\varphi_s +\frac{\sin(\lambda_t-\lambda_s)}{\sin p}\cos\varphi_t\cos\varphi_s\cos p; ::\tan p = \frac{\sin(\lambda_t-\lambda_s)\cos\varphi_t\cos\varphi_s}{\sin\varphi_t-\cos\theta_{s,t}\sin\varphi_s}.
Further details Because the brief derivation gives an angle between 0 and which does not reveal the sign (west or east of north ?), a more explicit derivation is desirable which yields separately the sine and the cosine of such that use of the
correct branch of the inverse tangent allows to produce an angle in the full range . The computation starts from a construction of the great circle between and . It lies in the plane that contains the sphere center, and and is constructed rotating by the angle around an axis . The axis is perpendicular to the plane of the great circle and computed by the normalized vector
cross product of the two positions: ::\mathbf{\omega} = \frac{1}{R^2\sin \theta_{s,t}}\mathbf{s}\times \mathbf{t} = \frac{1}{\sin \theta_{s,t}}\left(\begin{array}{c} \cos\varphi_s\sin\lambda_s\sin\varphi_t -\sin\varphi_s\cos\varphi_t\sin\lambda_t \\ \sin\varphi_s\cos\lambda_t\cos\varphi_t -\cos\varphi_s\sin\varphi_t\cos\lambda_s \\ \cos\varphi_s\cos\varphi_t\sin(\lambda_t-\lambda_s) \end{array}\right). A right-handed tilted coordinate system with the center at the center of the sphere is given by the following three axes: the axis , the axis ::\mathbf{s}_\perp = \omega \times \frac{1}{R}\mathbf{s} = \frac{1}{\sin\theta_{s,t}} \left(\begin{array}{c} \cos\varphi_t\cos\lambda_t(\sin^2\varphi_s+\cos^2\varphi_s\sin^2\lambda_s)-\cos\lambda_s(\sin\varphi_s\cos\varphi_s\sin\varphi_t+\cos^2\varphi_s\sin\lambda_s\cos\varphi_t\sin\lambda_t)\\ \cos\varphi_t\sin\lambda_t(\sin^2\varphi_s+\cos^2\varphi_s\cos^2\lambda_s)-\sin\lambda_s(\sin\varphi_s\cos\varphi_s\sin\varphi_t+\cos^2\varphi_s\cos\lambda_s\cos\varphi_t\cos\lambda_t)\\ \cos\varphi_s[\cos\varphi_s\sin\varphi_t-\sin\varphi_s\cos\varphi_t\cos(\lambda_t-\lambda_s)] \end{array}\right) and the axis . A position along the great circle is ::\mathbf{s}(\theta) = \cos\theta \mathbf{s}+\sin\theta \mathbf{s}_\perp,\quad 0\le\theta\le 2\pi. The compass direction is given by inserting the two vectors and and computing the gradient of the vector with respect to at . ::\frac{\partial}{\partial\theta}\mathbf{s}_{\mid \theta=0}=\mathbf{s}_\perp. The angle is given by splitting this direction along two orthogonal directions in the plane tangential to the sphere at the point . The two directions are given by the partial derivatives of with respect to and with respect to , normalized to unit length: ::\mathbf{u}_N = \left( \begin{array}{c} -\sin\varphi_s\cos\lambda_s\\ -\sin\varphi_s\sin\lambda_s\\ \cos\varphi_s \end{array}\right); ::\mathbf{u}_E = \left(\begin{array}{c} -\sin\lambda_s\\ \cos\lambda_s\\ 0 \end{array} \right); ::\mathbf{u}_N\cdot \mathbf{s} = \mathbf{u}_E\cdot \mathbf{u}_N =0 points north and points east at the position . The position angle projects into these two directions, ::\mathbf{s}_\perp = \cos p \,\mathbf{u}_N+\sin p\, \mathbf{u}_E, where the positive sign means the positive position angles are defined to be north over east. The values of the cosine and sine of are computed by multiplying this equation on both sides with the two unit vectors, ::\cos p = \mathbf{s}_\perp \cdot \mathbf{u}_N =\frac{1}{\sin\theta_{s,t}}[\cos\varphi_s\sin\varphi_t - \sin\varphi_s\cos\varphi_t\cos(\lambda_t-\lambda_s)]; ::\sin p = \mathbf{s}_\perp \cdot \mathbf{u}_E =\frac{1}{\sin\theta_{s,t}}[\cos\varphi_t\sin(\lambda_t-\lambda_s)]. Instead of inserting the convoluted expression of , the evaluation may employ that the
triple product is invariant under a circular shift of the arguments: ::\cos p = (\mathbf{\omega}\times \frac{1}{R}\mathbf{s})\cdot \mathbf{u}_N = \omega\cdot(\frac{1}{R}\mathbf{s}\times \mathbf{u}_N). If
atan2 is used to compute the value, one can reduce both expressions by division through and multiplication by , because these values are always positive and that operation does not change signs; then effectively ::\tan p = \frac{\sin(\lambda_t-\lambda_s)}{\cos\varphi_s\tan\varphi_t -\sin\varphi_s\cos(\lambda_t-\lambda_s)}. ==Finding way-points==