; It was found on 26 March 2005 that this algorithm is incapable of
; calculating azimuths to the south of 'startpoint'

function geodesicazimuth,startpoint,endpoint,nosouth=ns
f=1/298.257222101
e=sqrt(2*f-f^2)
a=6378137.0
phi=(startpoint[0]+endpoint[0])/2
n=a/sqrt(1-e^2*(sin(phi))^2)
m=a*(1-e^2)/(1-e^2*(sin(phi))^2)^(3/2)
eta=e*cos(phi)/sqrt(1-e^2)
v=sqrt(1+eta^2)
t=tan(phi)
dphi=endpoint[0]-startpoint[0]
dlambda=endpoint[1]-startpoint[1]
ssinalpha=n*dlambda*cos(phi)*(1-(dlambda*sin(phi))^2/24+(1+eta^2-9*eta^2*t^2)*dphi^2/(24*v^4))
scosalpha=m*dphi*cos(dlambda/2)*(1+(1-2*eta^2)*(dlambda*cos(phi))^2/24+(eta^2-(eta*t)^2)*dphi^2/(8*v^4))
dalpha=dlambda*sin(phi)*(1+(1+eta^2)*(dlambda*cos(phi))^2/12+(3+8*eta^2)*dphi^2/(24*v^4))
; The next line is from the original version of this function.
; alphabar=atan(ssinalpha/scosalpha)
; And the next line is the modification of 26 March 2005:
alphabar=atan(ssinalpha,scosalpha)
azimuth=alphabar-dalpha/2
; The next line was added on 29 March 2005 in order to provide quasi-backward compatibility via keyword nosouth
if keyword_set(ns) then if azimuth gt !dpi/2 then azimuth=azimuth-!dpi else if azimuth lt -!dpi/2 then azimuth=azimuth+!dpi
return,azimuth
end