Geodesic arc length

From Fig. 3 and Eqn. (17) we see that:

$$\cos\alpha=a(1-e^2\cos^2\theta)^{1/2}\frac{\mathrm{d}\theta}{\mathrm{d}s}=\pm(\cos^2\theta-\sin^2\alpha_0)^{1/2}/\cos\theta$$ (18)

and thus

$$\frac{\mathrm{d}s}{\mathrm{d}\theta}=\pm a\frac{\cos\theta(1-e^2\cos^2\theta)^{1/2}}{(\cos^2\alpha_0-\sin^2\theta)^{1/2}}$$ (19)

with the sign being that of $\cos\alpha$. We now substitute $\sin\theta=\sin\sigma\cos\alpha_0$, where $a\sigma$ is the arc-length along the great circle on the auxiliary sphere, measured from where it crosses the equator in a Northerly direction, obtaining:

$$\frac{\mathrm{d}s}{\mathrm{d}\sigma}=a(1-e^2\cos^2\theta)^{1/2}=b(1+u^2\sin^2\sigma)^{1/2}$$ (20)

where $u^2\equiv e^2\cos^2\alpha_0/(1-e^2)$. By expanding in a power series in $u^2$ and integrating term by term[3,5], we obtain $s$ in the form:

$$s/b = \sigma(1+u^2/4-3u^4/64+5u^6/256-175u^8/16384+\ldots)$$

$$-\sin 2\sigma (u^2/8)(1-u^2/4+15u^4/128-35u^6/512+\ldots)$$

$$-\sin 4\sigma (u^4/256)(1-3u^2/4+35u^4/64-\ldots)$$

$$-\sin 6\sigma (u^6/3072)(1-5u^2/4+\ldots)$$

$$-\sin 8\sigma(5u^8/131072)(1-\ldots)$$

$$- \ldots$$ (21)

The distance between two points $s(\sigma_2, \alpha_0)-s(\sigma_1, \alpha_0)$, on a geodesic arc is best obtained, after differencing Eqn. (21), by using the identity $\sin( 2n\sigma_2) -\sin( 2n\sigma_1)=2\cos (2n\sigma_m) \sin n(\sigma_2-\sigma_1)$, where $\sigma_m=(\sigma_1+\sigma_2)/2$. This avoids excessive loss of significant digits when the two points are close together.

Vincenty[5] has rearranged a subset of the resulting equations into nested forms more suitable for computation:

$$\tan\sigma_1 = \tan\phi_1/\cos\alpha_1$$

$$\sin\alpha_0 = \cos\phi_1\sin\alpha_1$$

$$u^2 = e^2 \cos^2\alpha_0/(1-e^2)$$

$$A = 1 +\frac{u^2}{16384}(4096+u^2(-768+u^2(320-175u^2)))$$

$$B = \frac{u^2}{1024}(256+u^2(-128+u^2(74-47u^2)))$$

$$\sigma_m = \sigma_1 +\sigma/2$$

$$\Delta\sigma = B\sin\sigma (\cos 2\sigma_m+\frac{B}{4}(\cos\sigma (-1+2\cos^2 2\sigma_m)$$

$$\quad -\frac{B}{6}\cos 2\sigma_m (-3 +4\sin^2\sigma)(-3+4\cos^2 2\sigma_m))))$$

$$s = bA(\sigma- \Delta\sigma)$$ (22)

In equations (22) the origins of $\sigma$ and $s$ have been shifted from the equator to the initial point ($1$), where the reduced latitude is $\phi_1$ and the azimuth of the geodesic is $\alpha_1$.