Great Circle on Spherical Earth
The shortest path between two points on the (spherical) Earth lies on the great circle going through these two points. The great circle is  as the name suggests  the largest possible circle on the sphere. All great circles of a given sphere have both the same radius and the same center as the sphere. The great circle's plane passes through the center of the sphere and divides it into 2 equal hemispheres.
Great Circle between Two Points
The equation of the great circle is computed most easily in cartesian coordinates. So, if the points were given in spherical coordinates we would transform them to cartesian coordinates first.
Let
p_{1} = (x_{1}, y_{1}, z_{1})^{T} 
(1a) 
p_{2} = (x_{2}, y_{2}, z_{2})^{T} 
(1b) 
be the vectors corresponding to the two points P_{1} and P_{2}, respectively. The great circle is located in the plane OP_{1}P_{2} spanned by these two vectors.
If we find two orthonormal vectors u and v in this plane then the equation of the great circle will be
c = r (u cos ω + v sin ω) 
(2) 
0° ≤ ω < 360° 
The obvious choice for the vector u is the normalized vector p_{1}:
u = p_{1} / R. 
(3) 
The vector product
w = p_{1} x p_{2} 
(4) 
is orthogonal to the plane OP_{1}P_{2}. The unit vector
v = u x w / w. 
(5) 
is in the plane OP_{1}P_{2} and is orthogonal to u.
Another way to compute the vector v is by using the GramSchmidt orthonormalization process. The angle δ between p_{1} and p_{2} can be computed using the scalar product:
p_{1} · p_{2} = R^{2} cos δ. 
(6) 
The projection of p_{2} to u is:
p_{u} = R cos(δ) u. 
(7) 
The vector
p_{v} = p_{2}  p_{u} 
(8) 
is orthogonal to p_{1} and is in the plane OP_{1}P_{2}. Hence the normalized vector v is our searched vector:
v = p_{v} / p_{v} 
(9) 
The above calculations break if the vectors p_{1} and p_{2} are parallel, because then they do not specify a plane. This happens if the corresponding points are either identical or antipodal. In this case there are infinite many great circles through P_{1} and P_{2}. Their equations are still given by
c = r (u cos ω + v sin ω) 
(10) 
where v is an arbitrary unit vector orthogonal to u.
Shortest Distance between Two Points
The shortest path between two points on the earth surface is the circular arc on the great circle between these points. So, the shortest distance d is
d = R δ 
(11) 
where R is the earth radius and δ is the angle between the two points. From the equation (6) it follows:
cos δ = p_{1} · p_{2} / R^{2} = (x_{1}x_{2} + y_{1}y_{2} + z_{1}z_{2}) / R^{2}. 
(12) 
By substituting spherical coordinates and using the identity:
cos(λ_{1}  λ_{2}) = cosλ_{1} cosλ_{2} + sinλ_{1} sinλ_{2} 
(13) 
we get:
cos δ = cos(λ_{1}  λ_{2}) cos φ_{1} cos φ_{2} + sin φ_{1} sin φ_{2} 
(14) 
Finally, the shortest distance d is:
d = R δ = R acos[cos(λ_{1}  λ_{2}) cos φ_{1} cos φ_{2} + sin φ_{1} sin φ_{2}]. 
(15) 
Heading (Bearing, Azimuth)
Heading (also called bearing or azimuth) is the angle α between a great circle and a meridian. Heading plays an important role in navigation. When we travel from the point P_{1} to the point P_{2} on the shortest path, then the heading at point P_{1} specifies the starting direction of our journey. The heading changes during the voyage, of course.
Heading as angle between planes
Heading α is the angle between the great circle plane OP_{1}P_{2} and the meridian plane OP_{1}P_{2}, or, equivalently between the respective normal vectors of these two planes. The vector orthogonal to the meridian plane is
m = p_{1} x (0, 0, 1)^{T} = (y_{1}, x_{1}, 0)^{T} 
(16) 
The vector orthogonal to the path plane is w = (w_{x}, w_{y}, w_{z})^{T} = p_{1} x p_{1}.
The dot product of the normalized vectors w and m equals the cosine of the heading angle α:
cosα = w/w · m/m 
(17) 
When we can take arccos() we get the angle α, 0 ≤ α ≤ 180°. The usual convention is that a positive angle denotes travelling eastwards while a negative angle denotes travelling westwards. Under the assumption that we travel along the shorter arc from P_{1} to P_{2} we can determine the direction by considering the vector product w. When it points upwards with respect to the positive zaxis the journey goes eastwards, when it points downwards the journey goes westwards. So, our formula for the heading angle is finally:
α = acos[cos(α)] if w_{z} ≥ 0 
(18a) 
α = acos[cos(α)] if w_{z} < 0 
(18a) 
Simplification of the geometry
To simplify the following calculations I will make two assumptions without loss of generality:
 The radius of the earth is 1. This is justified because only angles are considered.
 The point P_{1} has the longitude 0 and the point P_{2} has the longitude dλ = λ_{2}  λ_{1}. I just start counting the longitude at the meridian passing through P_{1}. The vector p_{1} is hence (x_{1}, 0, z_{1})^{T}.
Under the simplifying assumptions we get:
m = (0, x_{1}, 0)^{T} 
(19) 
and
w_{x} = z_{1} y_{2} 
(20a) 
w_{y} = x_{1} z_{2} + z_{1} x_{2} 
(20b) 
w_{z} = x_{1} y_{2} 
(20c) 
and
cosα = w_{y} / w 
(21) 
Heading as angle in a spherical triangle
The heading is the angle at the vertex P_{1} of the spherical triangle NP_{1}P_{2}, where N is the north pole. The sides of this triangle are arcs of the respective great circles and their lengths are equal to the respective angles at the center of the Earth.
Now we can apply the cosine rule formula for the side P_{2}N that is opposite to the angle α (see the equation (10) in Weisstein, Eric W. Spherical Trigonometry in MathWorld):
cosΦ_{2} = cosΦ_{1} cosδ + sinΦ_{1} sinδ cosα 
(22) 
Solving for cosα and taking into account that
Φ_{1} = 90°  φ_{1} 
(23a) 
Φ_{2} = 90°  φ_{2} 
(23b) 
we get:
cosα = 

(24) 
Tangens formula
Let us have a look at the geometry of the vector w:
In the equation (21) we computed in fact α using the cosine formula in the triangle OBC. Now, we will compute it using the tangens formula:
tanα =  w_{xz} / w_{y} 
(25) 
Here, w_{xz} is the projection of the vector w onto the XZplane:
w_{xz}^{2} = w_{x}^{2} + w_{z}^{2} = z_{1}^{2} y_{2}^{2} + x_{1}^{2} y_{2}^{2} = y_{2}^{2}(x_{1}^{2} + z_{1}^{2}) = y_{2}^{2} 
(26) 
The last equality follows from the fact that p_{1} = 1 and y_{1} = 0. Substituting (26) and (20b) into (25) we get:
tanα = 

(27) 
Finally, using spherical coordinates we obtain:
tanα = 

(28) 
Even if this formula looks more complicated than the formula (24) it is often recommended because it avoids the calculation of the angle δ between the two points P_{1} and P_{2}. However, significant care must be taken by applying atan() to the equation (28). A simpler way is to use atan2() instead of atan():
dl = lon2  lon1 y2 = sin(dl)*cos(lat2) wy = cos(lat1)*sin(lat2)  sin(lat1)*cos(lat2)*cos(dl) alpha = atan2(y2, wy)
References
 Wikipedia: Greatcircle distance
 Chris Veness: Calculate distance, bearing and more between Latitude/Longitude points
 Weisstein, Eric W. "Great Circle." From MathWorldA Wolfram Web Resource
 Weisstein, Eric W. "Spherical Trigonometry" From MathWorldA Wolfram Web Resource
 PROJ.4  Geodesic Calculations
 Ed Williams: Aviation Formulary
 James R. Clynch: Paths Between Points on Earth
 Mathforum: Latitude & Longitude
 Mathforum: Bearing Between Two Points