(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 '================================================================================= ' Copyright lost_species 2008 LGPL http://www.fsf.org/licensing/licenses/lgpl.html '================================================================================= 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 L = toRad(lon2 - lon1) 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