Fórmules de Vincenty

De Viquipèdia
Dreceres ràpides: navegació, cerca

Les fórmules de Vincenty formen un algorisme molt eficient per al càlcul de la distància entre dos punts de la superfície d'un el·lipsoide de revolució. Són utilitzades àmpliament en geodèsia per calcular distàncies sobre la superfície de la Terra a causa del fet que requereix un nombre petit d'operacions tot i donar una precisió de 0.5mm (0,000015 "), molt millor que el mètode tradicional de la fórmula del haversine utilitzada anteriorment en trigonometria esfèrica. L'algorisme va ser publicat per Thaddeus Vincenty el 1975.

Utilitza un mètode iteratiu.

Algorisme[modifica | modifica el codi]

El següent algorisme expressa les fórmules d'una forma senzilla de calcular:

a, b = semieixos major i menor del el·lipsoide
φ1, φ2 = latitud geodèsica
L = diferència en longitud
f = aplatament del el·lipsoide (a-b)/a
U1 = atan((1−f).tanφ1) (U is ‘latitud reduïda’) 
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)]} (emprar cos2σm=0 si se està al llarg de l'equador)
}
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λ)
On:
*s és la distancia (mateixes unitats que a, b)
*α1 és l'azimut inicial
*α2 és l'azimut final (en direcció p1→p2)

La fórmula pot no tenir solució per a dos punts gairebé antipodals. Limitar el nombre d'iteracions per evitar aquest cas.

L'el·lipsoide més estès és el WGS84, per al qual:

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

Exemple per una prova[modifica | modifica el codi]

Cas de prova extret de Geoscience Austràlia, usant WGS-84:

Flinders Peak 37 ° 57'03 72.030 "S, 144 ° 25'29 52.440" E
Buninyong 37 ° 39'10 15.610 "S, 143 ° 55'35 38.390" E
Resultat:
s 54.972,271 m
α1 306 ° 52'05 37 "
α2 127 ° 10'25 07 "(≡ 307 · 10'25 07" p1 → p2)

Referències[modifica | modifica el codi]

Direct and Inverse Solutions of Geodesics on the ellipse with application of nested equations . Publicació original de Vincenty.

Enllaços externs[modifica | modifica el codi]

Codi font disponible (llicència LGPL) per Chris Veness