Let the radius of Earth (assumed perfectly spherical) be $r_E$. Attach a coordinate frame to Earth, with its center at the center of Earth, and its $z$ axis passing through the two poles of Earth. And let the $x$ axis pass through the equator at Greenwich longitude. In this coordinate frame the coordinates of points on the surface of Earth are given by
$ P = r_E ( \cos \theta \cos \phi, \cos \theta \sin \phi, \sin \theta ) $
where $\theta$ is the Latitude angle (positive in the northern hemisphere, and negative in the southern hemisphere), and $\phi$ is the Longitude (positive east of Greenwich and negative west of it).
Now suppose we have two points $A$ and $B$ on the surface of Earth whose Latitudes and Longitudes are known, then we can compute their Cartesian coordinates from the above formula. The normal vector to the great circle connecting $A$ with $B$ is along the unit vector
$ U = \dfrac{A \times B}{\| A \times B \|} $
Orthogonalizing the vectors $A$ and $B$, by taking $A$ as a reference, then the orthogonal vector to $A$ lying on this great circle is given by
$ C = U \times A $
Now we can express the parametric equation of the great circle as
$ P(t) = A \cos t + C \sin t $
Note that $P(0) = A $
The tangent vector to the great circle $P(t)$ at $A$ is given by $P'(0)$
$P'(t) = - A \sin t + C \cos t $
Therefore
$P'(0) = C $
To find our directions at point $A$ to point in the direction of $C$, we have to build a local coordinate system as follows
The local North unit vector is given by
$ N = \dfrac{ \partial (\cos \theta \cos \phi, \cos \theta \sin \phi, \sin \theta )}{\partial \theta } = ( -\sin \theta \cos \phi , -\sin \theta \sin \phi, \cos \theta ) $
The local East unit vector is given by
$ E = (- \sin \phi , \cos \phi , 0 ) $
And the radial unit vector is
$ V = E \times N = ( \cos \theta \cos \phi, \cos \theta \sin \phi, \sin \theta) $
Define the rotation matrix corresponding to this local frame as R
$ R = [E , N , V ] $
Then the local unit direction in the basis $R$ is
$ C' = \dfrac{1}{r_E} R^T C $
The $z$ component of $C'$ is zero. And the first and second coordinates of $C'$ determine the components of this unit vector along the East and North directions respectively. The angle $\psi$ from the North (the bearing) is
$ \psi = ATAN2( C'_y , C'_x) $
Example: Let $A$ represent Ottawa, Canada. The latitude is $45.424721^\circ$, and the longitude is $-75.695^\circ $
And let $B$ represent Mecca, Saudi Arabia. The latitude is $21.426388^\circ$ and the longitude is $39.82555^\circ$
With these values the bearing is calculated as
$ \psi = 57.1667^\circ $ (Clock wise from local North).
The following program summarizes these calculations outlined above.
Public Sub mecca()
' rE is kept untouched
Dim a(3), b(3), u(3), c(3), cp(3) As Double
Dim n(3), e(3), v(3) As Double
Dim r(3, 3), rt(3, 3) As Double
p = WorksheetFunction.Pi()
th1 = 45.424721 * p / 180
phi1 = -75.695 * p / 180
th2 = 21.426388 * p / 180
phi2 = 39.82555 * p / 180
' calculate A and B
a(1) = Cos(th1) * Cos(phi1)
a(2) = Cos(th1) * Sin(phi1)
a(3) = Sin(th1)
b(1) = Cos(th2) * Cos(phi2)
b(2) = Cos(th2) * Sin(phi2)
b(3) = Sin(th2)
Call cross(a, b, u)
Call normalize1(u)
Call cross(u, a, c)
For i = 1 To 3
v(i) = a(i)
Next i
e(1) = -Sin(phi1)
e(2) = Cos(phi1)
e(3) = 0
n(1) = -Sin(th1) * Cos(phi1)
n(2) = -Sin(th1) * Sin(phi1)
n(3) = Cos(th1)
For i = 1 To 3
r(i, 1) = e(i)
r(i, 2) = n(i)
r(i, 3) = v(i)
Next i
Call transpose(3, 3, r, rt)
Call multiplyv(3, 3, rt, c, cp)
psi = WorksheetFunction.Atan2(cp(2), cp(1))
MsgBox (psi * 180 / p)
End Sub