Fórmulas de Vincenty

Las fórmulas de Vincenty^[1] forman un algoritmo muy eficiente para el cálculo de la distancia entre dos puntos de la superficie de un elipsoide de revolución. Son utilizadas ampliamente en geodesia para calcular distancias sobre la superficie de la Tierra debido a que requieren un número de operaciones bajo a pesar de dar una precisión de 0,5mm (0.000015″), mucho mejor que el método tradicional de la fórmula del semiverseno usada en trigonometría esférica. El algoritmo fue publicado por Thaddeus Vincenty en 1975.

Usa un método iterativo.

Algoritmo

El siguiente algoritmo expresa las fórmulas de una forma sencilla de calcular:

   a, b = semiejes mayor y menor del elipsoide
   φ1, φ2 = latitud geodética
   L = diferencia en longitud
   f = achatamiento del elipsoide (a−b)/a
   U1 = atan((1−f).tanφ1) (U is ‘latitud reducida’) 
   U2 = atan((1−f).tanφ2)
   λ = L, λ′ = 2π
   do while abs(λ−λ′) > 10^-12    (implica un error < 0.06mm)
   {
     sinσ = √[ (cosU2.sinλ)² + (cosU1.sinU2 − sinU1.cosU2.cosλ)² ]
     cosσ = sinU1.sinU2 + cosU1.cosU2.cosλ
     σ = atan2(sinσ, cosσ)
     sinα = cosU1.cosU2.sinλ / sinσ
     cos²α = 1 − sin²α (trig identity; §6)      
     cos2σm = cosσ − 2.sinU1.sinU2/cos²α
     C = f/16.cos²α.[4+f.(4−3.cos²α)]
     λ′ = λ
     λ = L + (1−C).f.sinα.{σ+C.sinσ.[cos2σm+C.cosσ.(−1+2.cos²2σm)]}   (úsese cos2σm=0 si se está a lo largo del ecuador)
   }
   u² = cos²α.(a²−b²)/b²
   A = 1+u²/16384.{4096+u².[−768+u².(320−175.u²)]}
   B = u²/1024.{256+u².[−128+u².(74−47.u²)]}
   Δσ = B.sinσ.{cos2σm+B/4.[cosσ.(−1+2.cos²2σm) − B/6.cos2σm.(−3+4.sin²σ).(−3+4.cos²2σm)]}
   s = b.A.(σ−Δσ)
   α1 = atan2(cosU2.sinλ, cosU1.sinU2 − sinU1.cosU2.cosλ)
   α2 = atan2(cosU1.sinλ, −sinU1.cosU2 + cosU1.sinU2.cosλ)

   Donde:
   *s  es la distancia (mismas unidades que a, b)
   *α1 es el azimuth inicial
   *α2 es el azimut final (en dirección p1→p2)

La fórmula puede no tener solución para dos puntos casi antipodales. Limitar el número de iteraciones para evitar este caso.

El elipsoide más extendido es el WGS84, para el cual:

   a = 6 378 137 m     (±2 m)
   b = 6 356 752.3142 m
   f = 1 / 298.257223563

Ejemplo para prueba

Caso de prueba extraído de Geoscience Australia, usando WGS-84:

   Flinders Peak    37°57′03.72030″S, 144°25′29.52440″E
   Buninyong        37°39′10.15610″S, 143°55′35.38390″E
   Resultado:
   s    54.972,271 m
   α1   306°52′05.37″
   α2   127°10′25.07″ (≡ 307°10′25.07″ p1→p2)