## July 30th, 2008

### (no subject)

If you happen to need the visual basic code to calculate the distance between 2 co-ordinates on the earths surface taking into account that the earth is not a perfect sphere but ignoring hills. Then you are in luck!

```Option Explicit

Private Const PI = 3.14159265358979
Private Const EPSILON As Double = 1e-12

Public Sub test()
MsgBox distVincenty(52.874, 4.389, 45.001, 15.716)
End Sub

Private Function toRad(ByVal degrees As Double) As Double
toRad = degrees * (PI / 180)
End Function

Private Function Atan2(ByVal X As Double, ByVal Y As Double) As Double
' code nicked from:
' http://en.wikibooks.org/wiki/Programming:Visual_Basic_Classic/Simple_Arithmetic#Trigonometrical_Functions
' If you re-use this watch out the x and y have been reversed.
If Y > 0 Then
If X >= Y Then
Atan2 = Atn(Y / X)
ElseIf X <= -Y Then
Atan2 = Atn(Y / X) + PI
Else
Atan2 = PI / 2 - Atn(X / Y)
End If
Else
If X >= -Y Then
Atan2 = Atn(Y / X)
ElseIf X <= Y Then
Atan2 = Atn(Y / X) - PI
Else
Atan2 = -Atn(X / Y) - PI / 2
End If
End If
End Function

Public Function distVincenty(ByVal lat1 As Double, ByVal lon1 As Double, ByVal lat2 As Double, ByVal lon2 As Double) As Double
'=================================================================================
' Calculate geodesic distance (in m) between two points specified by latitude/longitude (in numeric degrees)
' using Vincenty inverse formula for ellipsoids
'=================================================================================
' Code has been ported by lost_species from www.aliencoffee.co.uk to VBA from javascript published at:
' http://www.movable-type.co.uk/scripts/latlong-vincenty.html
' * from: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the
' *       Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975
' *       http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
'=================================================================================
'=================================================================================

Dim low_a As Double
Dim low_b As Double
Dim f As Double
Dim L As Double
Dim U1 As Double
Dim U2 As Double
Dim sinU1 As Double
Dim sinU2 As Double
Dim cosU1 As Double
Dim cosU2 As Double
Dim lambda As Double
Dim lambdaP As Double
Dim iterLimit As Integer
Dim sinLambda As Double
Dim cosLambda As Double
Dim sinSigma As Double
Dim cosSigma As Double
Dim sigma As Double

Dim sinAlpha As Double
Dim cosSqAlpha As Double
Dim cos2SigmaM As Double
Dim C As Double
Dim uSq As Double
Dim upper_A As Double
Dim upper_B As Double
Dim deltaSigma As Double
Dim s As Double

low_a = 6378137
low_b = 6356752.3142
f = 1 / 298.257223563 'WGS-84 ellipsiod

U1 = Atn((1 - f) * Tan(toRad(lat1)))
U2 = Atn((1 - f) * Tan(toRad(lat2)))
sinU1 = Sin(U1)
cosU1 = Cos(U1)
sinU2 = Sin(U2)
cosU2 = Cos(U2)

lambda = L
lambdaP = 2 * PI
iterLimit = 20

While (Abs(lambda - lambdaP) > EPSILON) And (iterLimit > 0)
iterLimit = iterLimit - 1

sinLambda = Sin(lambda)
cosLambda = Cos(lambda)
sinSigma = Sqr(((cosU2 * sinLambda) ^ 2) + ((cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ^ 2))
If sinSigma = 0 Then
distVincenty = 0  'co-incident points
Exit Function
End If
cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
sigma = Atan2(cosSigma, sinSigma)
sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
cosSqAlpha = 1 - sinAlpha * sinAlpha

If cosSqAlpha = 0 Then 'check for a divide by zero
cos2SigmaM = 0 '2 points on the equator
Else
cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha
End If

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 ^ 2))))
Wend

If iterLimit < 1 Then
MsgBox "iteration limit has been reached, something didn't work."
Exit Function
End If

uSq = cosSqAlpha * (low_a ^ 2 - low_b ^ 2) / (low_b ^ 2)
upper_A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
upper_B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
deltaSigma = upper_B * sinSigma * (cos2SigmaM + upper_B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM ^ 2) _
- upper_B / 6 * cos2SigmaM * (-3 + 4 * sinSigma ^ 2) * (-3 + 4 * cos2SigmaM ^ 2)))

s = low_b * upper_A * (sigma - deltaSigma)
distVincenty = s

End Function
```

Useful resources