L'algorisme de Vincenty és l'estàndard d'or per a la geodèsia perquè modela la Terra com un el·lipsoide oblat (WGS-84), a diferència de mètodes més simples com Haversine que assumeixen una esfera.
Aquesta precisió requereix resoldre un problema complex: la distància geodèsica depèn de la longitud sobre l'el·lipsoide, la qual al seu torn depèn de la distància. Aquesta dependència circular es resol mitjançant un procés iteratiu.
Abans de començar, definim les constants de l'el·lipsoide de referència (habitualment WGS-84):
Per simplificar els càlculs, projectem la latitud geodèsica ($\phi$) sobre una esfera auxiliar. Això s'anomena Latitud Reduïda ($U$).
Obtenim $U_1$ (punt inici) i $U_2$ (punt final). També precalculem els seus sinus i cosinus, ja que els farem servir constantment:
Aquest és el punt de partida de la iteració:
Entrem al cor de l'algorisme. Repetim els següents passos fins que $\lambda$ deixi de canviar (convergeixi).
Amb el valor actual de $\lambda$, calculem la distància angular entre els punts sobre l'esfera auxiliar.
Primer, component sinus (utilitzant identitats vectorials per precisió):
$$ \sin \sigma = \sqrt{(\cos U_2 \sin \lambda)^2 + (\cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos \lambda)^2} $$Segon, component cosinus:
$$ \cos \sigma = \sin U_1 \sin U_2 + \cos U_1 \cos U_2 \cos \lambda $$Finalment, l'angle $\sigma$ (en radians):
$$ \sigma = \arctan2(\sin \sigma, \cos \sigma) $$Determina la direcció del camí en el punt on creua l'equador. Si $\sin \sigma = 0$ (punts coincidents), evitem la divisió per zero.
$$ \sin \alpha = \frac{\cos U_1 \cos U_2 \sin \lambda}{\sin \sigma} $$Necessitem el cosinus quadrat per al següent pas (identitat trigonomètrica):
$$ \cos^2 \alpha = 1 - \sin^2 \alpha $$Determina on es troba el centre de la línia respecte a l'equador. Això és crucial perquè la curvatura canvia segons la latitud.
$$ \cos(2\sigma_m) = \cos \sigma - \frac{2 \sin U_1 \sin U_2}{\cos^2 \alpha} $$Aquest coeficient quantifica l'efecte de l'aplatament $f$ en funció de la direcció del camí ($\cos^2 \alpha$).
$$ C = \frac{f}{16} \cos^2 \alpha \cdot [4 + f (4 - 3 \cos^2 \alpha)] $$Recalculem $\lambda$ afegint un terme de correcció basat en tot el que hem calculat anteriorment.
$$ \lambda_{nou} = L + (1 - C) \cdot f \cdot \sin \alpha \cdot \left\{ \sigma + C \sin \sigma \left[ \cos(2\sigma_m) + C \cos \sigma (-1 + 2 \cos^2(2\sigma_m)) \right] \right\} $$Després de calcular $\lambda_{nou}$, comprovem:
$$ |\lambda_{nou} - \lambda_{antic}| < 10^{-12} $$Si la diferència és menor que aquest llindar (aprox. 0.06mm), aturem el bucle. Si no, actualitzem $\lambda_{antic} = \lambda_{nou}$ i tornem al Pas A.
Un cop hem sortit del bucle (tenim el $\lambda$ correcte), calculem la distància final en metres.
Definim dos nous paràmetres derivats de les dimensions de l'el·lipsoide ($u^2$) i constants de correcció ($A$ i $B$):
$$ u^2 = \cos^2 \alpha \cdot \frac{a^2 - b^2}{b^2} $$ $$ A = 1 + \frac{u^2}{16384} (4096 + u^2 (-768 + u^2 (320 - 175 u^2))) $$ $$ B = \frac{u^2}{1024} (256 + u^2 (-128 + u^2 (74 - 47 u^2))) $$Calculem el terme delta sigma ($\Delta \sigma$):
$$ \Delta \sigma = B \sin \sigma \left\{ \cos(2\sigma_m) + \frac{B}{4} \left[ \cos \sigma (-1 + 2 \cos^2(2\sigma_m)) - \frac{B}{6} \cos(2\sigma_m) (-3 + 4 \sin^2 \sigma) (-3 + 4 \cos^2(2\sigma_m)) \right] \right\} $$La distància geodèsica $s$ és:
$$ s = b \cdot A \cdot (\sigma - \Delta \sigma) $$