Fórmulas de Vincenty
Las fórmulas de Vincenty 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 requiere 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)
Referencias
Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations. Publicación original de Vincenty.
Enlaces externos
Código fuente disponible (licencia LGPL) por Chris Veness