lost_species (lost_species) wrote,
lost_species
lost_species

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

Tags: code, gis, gps, vba
Subscribe

  • This is a recipe

    450 g porridge oats 200 g plain flour 200 g wholemeal flour 2 tsp salt 2 tsp sugar (muscovado) 2 tsp yeast 900 ml warm water 900 ml milk…

  • www.theyworkforyou.com

    FOR THE ATTENTION OF: Mark Hendrick MP Preston Dear Mark Hendrick, In regards to your application for my vote for you during the upcoming election…

  • Hotpoint Repair Services Consumer Review

    So over a month ago my fridge freezer developed an identity crisis and the fridge temperature dropped to a fairly steady -5 deg C. My max/min…

  • Post a new comment

    Error

    default userpic

    Your reply will be screened

    Your IP address will be recorded 

    When you submit the form an invisible reCAPTCHA check will be performed.
    You must follow the Privacy Policy and Google Terms of use.
  • 1 comment