function geodesic,startpoint,azimuth,distance

s=distance
alpha0=azimuth
phi0=startpoint[0]     ;latitude
lambda0=startpoint[1]  ;longitude
f=1/298.257222101
e=sqrt(2*f-f^2)
a=6378137.0
n=a/sqrt(1-e^2*(sin(phi0))^2)
m=a*(1-e^2)/(1-e^2*(sin(phi0))^2)^(3/2)
convergence=1

;first guess

endpoint=startpoint+[s*sin(alpha0)/(n*cos(phi0)),s*cos(alpha0)/m]

while (convergence gt 0.0000017453292) do begin

oldpoint=endpoint
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)^1.5
eta=sqrt((e*cos(phi))^2/(1-e^2))
v=sqrt(1+eta^2)
t=tan(phi)
change=endpoint-startpoint
;alphachange=change[1]*sin(phi)*(1+(1+eta^2)*(change[1]*cos(phi))^2/12+(3+8*eta^2)*change[0]^2/(24*v^4))
alpha=alpha0+change[1]*sin(phi)*(1+(1+eta^2)*(change[1]*cos(phi))^2/12+(3+8*eta^2)*change[0]^2/(24*v^4))/2
endpoint[0]=phi0+s*cos(alpha)*(1.-(1.-2.*eta^2)*(change[1]*cos(phi))^2/24.-eta^2*(1-t^2)*change[0]^2/(8.*v^4))/(m*cos(change[1]/2))
endpoint[1]=lambda0+s*sin(alpha)*(1+(change[1]*sin(phi))^2/24-(1+eta^2-9*eta^2*t^2)*change[0]^2/(24*v^4))/(n*cos(phi))
convergence=mag(endpoint-oldpoint)

endwhile

return,endpoint


end