There is a common formula for computing the initial bearing for a course from lat1, lon1 to lat2, lon2. (And yes, I'm aware that the course will change when traveling).
atan2( sin(lon2-lon1) * cos(lat2), cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon2-lon1))
The formulary where I took this from mostly uses a spherical model for simplicity.
My questions about this formula are:
What is the appropriate formula for an ellipsoid model such as WGS84?
If it is the same formula, does this hold generally? I.e. is it the same formula for any ellipsoid model because we deformed the coordinate system, not the coordinates?
Somewhat like Vincentry's formula, which uses a reduced latitude u = atan((1 - f) * tan(lat)). Also note that I'm looking for a global solution, including for the poles.