Fòrmules de Vincenty

Les fòrmules de Vincenty formen un algoritme molt eficient per al càlcul de la distància entre dos punts de la superfície d'un elipsoide de revolució. Són utilisades àmpliament en geodèsia per a calcular distàncies sobre la superfície de la Terra degut a que requerixen un número d'operacions baix a pesar de donar una precisió de 0,5mm (0.000015″), molt millor que el método tradicional de la fòrmula del semiverseno usada en trigonometria esfèrica. L'algoritme va ser publicat per Thaddeus Vincenty en 1975.

Usa un método iterativo.

Algoritme

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

   a, b = semiejes major i menor de l'elipsoide
   φ1, φ2 = latitut geodética
   L = diferència en llongitut
   f = achatamiento de l'elipsoide (a−b)/a U1
   = nuguen((1−f).tanφ1) (O is ‘latitut reduïda’) 
   U2 = nuguen((1−f).tanφ2)
   λ = L, λ′ = 2π
   do while abs(λ−λ′) > 10-12    (implica un error < 0.06mm)
   {
     senseσ = √[ (cosU2.senseλ)² + (cosU1.sinU2 − sinU1.cosU2.cosλ)² ]
     cosσ = sinU1.sinU2 + cosU1.cosU2.cosλ
     σ = nuguen2(senseσ, cosσ)
     senseα = cosU1.cosU2.senseλ / senseσ
     cos²α = 1 − sense²α (trig identity; §6)      
     cos2σm = cosσ − 2.sinU1.sinU2/cos²α
     C = f/16.cos²α.[4+f.(4−3.cos²α)]
     λ′ = λ
     λ = L + (1−C).f.senseα.{σ+C.senseσ.[cos2σm+C.cosσ.(−1+2.cos²2σm)]}   (use's cos2σm=0 si s'està a lo llarc de l'equador)
   }
   o² = cos²α.(a²−b²)/b²
   A = 1+o²/16384.{4096+o².[−768+o².(320−175.o²)]}
   B = o²/1024.{256+o².[−128+o².(74−47.o²)]}
   Δσ = B.senseσ.{cos2σm+B/4.[cosσ.(−1+2.cos²2σm) − B/6.cos2σm.(−3+4.sense²σ).(−3+4.cos²2σm)]}
   s = b.A.(σ−Δσ)
   α1 = nuguen2(cosU2.senseλ, cosU1.sinU2 − sinU1.cosU2.cosλ)
   α2 = nuguen2(cosU1.senseλ, −sinU1.cosU2 + cosU1.sinU2.cosλ)

   A on:
   s  és la distància (mateixes unitats que a, b)
   *α1 és el azimuth inicial
   *α2 és el azimut final (en direcció p1→p2)

La fòrmula pot no tindre solució per a dos punts casi antipodales. Llimitar el número de iteraciones per a evitar este cas.

L'elipsoide 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

Eixemple per a prova

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

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