Wiki Home

Thadeus Vincenty Algorithm


Namespace: WIN_COM_API
Vincenty's formulae is used to calculate the distance between two points on the surface of an spheroid like The Earth.

Here's a VFP implementation of that, courtesy of long-time friend of the Wiki, Dr. F. N. Kautzmann.

Function distVincenty(lat1, lon1, lat2, lon2)
* Dr. F. N. Kautzmann, III
* Program date: June 2001
* Converted from MATLAB/public sources
* Free for non-commercial use
  SET DECIMALS TO 18
  PI = 4 * Atan(1)
  a = 6378137
  b = 6356752.3142
  f = 1 / 298.257223563
  L = ((lon2 - lon1) * (4 * Atan(1) / 180))
  U1 = Atan((1 - f) * Tan((lat1* (4 * Atan(1) / 180))))
  U2 = Atan((1 - f) *Tan((lat2 * (4 * Atan(1) / 180))))
  sinU1 = Sin(U1)
  cosU1 = Cos(U1)
  sinU2 = Sin(U2)
  cosU2 = Cos(U2)
  **********  For Antipodal Points ****
  *** Correct for antipodal errors at exact poles by adjusting 0.6 millimeters ****
  if abs(pi/2-abs(lat1)) < 1e-12
    lat1 = sign(lat1)*(pi/2-(1e-12))
     endif
     if abs(pi/2-abs(lat2)) < 1e-12
    lat2 = sign(lat2)*(pi/2-(1e-12))
  endif
******************  Vars *************************************
  Lambda = L
  lambdaP = 2 * (4 * atan(1)) && '2*PI
  iterLimit = 50
  ******** Main *******
  DO While Abs(Lambda - lambdaP) > 0.000000000001 And iterLimit > 0
    iterLimit = iterLimit - 1
    sinLambda = Sin(Lambda)
    cosLambda = Cos(Lambda)
    sinSigma = Sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 *   cosLambda))
    If (sinSigma = 0) Then
        s = 0 && coincident points
    EndIf
    cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
    sigma = atn2(sinSigma, cosSigma)
    sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
            cosSqAlpha = 1 - sinAlpha * sinAlpha
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha
            If isdigit(STR(cos2SigmaM))=.T.
            Else
                cos2SigmaM = 0  && equatorial line: cosSqAlpha=0
            EndIf
            C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
            lambdaP = Lambda
            Lambda = L + (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)))
     Enddo
        If (iterLimit = 0)
            s = 0  && formula failed to converge
            answer='Failed to Converge'
            RETURN answer
        EndIf
        uSq = cosSqAlpha * (a * a - b * b) / (b * b)
        VinA = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
        VinB = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
        deltaSigma = VinB * sinSigma * (cos2SigmaM + VinB / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - VinB / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)))
        s = b * VinA * (sigma - deltaSigma)
        ** Return **
        answer=s
      RETURN answer
ENDFUNC
****** sample data
lat1 = 29.7730050010
lon1 = -95.7272334343
lat2 = 29.77282732323
lon2 = -95.7272276567
Answer w/o bearing is:  19.7201133002100620000 millimeters distance between lat/lon pairs
Note: Be certain enter the sign of the longitudes as applicable
( Topic last updated: 2009.11.15 04:48:03 PM )