DONNE LA DISTANCE (EN M) ENTRE DEUX POINTS DE COORDONNÉES GÉOGRAPHIQUES CONNUES

Attribute VB_Name = "modDistance_2006_05_29_B"
Option Explicit

Public Const PI      As Double = 3.14159265358979     
Public Const PI_DEMI As Double = 1.5707963267949
Public Const EPSILON As Double = 0.00001 ' 1E-05       

Type Pt_GEOGRAPHIQUE
    Lon      As Double ' radian
    Lat      As Double ' radian
    LonDegDD As Double ' decimal degree
    LatDegDD As Double ' decimal degree
End Type



Function dtR(deg As Double) As Double
    dtR = deg * (PI / 180)
End Function

Function arctan2(y As Double, x As Double) As Double
    If Abs(x) < 0.000000000000001 Then   '1E-15
        ' donc x = 0
        arctan2 = IIf(y < 0, -PI_DEMI, PI_DEMI)
    Else
        arctan2 = Arctan(y / x) + IIf(x > 0#, 0#, IIf(y < 0#, -PI, PI))
    End If
End Function
'---------------------------------------------------------------------------------------
' Procedure : distVincenty
' DateTime  : 29-05-2006  20:26
' Author    : CADRATURE
' Purpose   : Donne la distance (en m) entre deux points de coordonnées géographiques connues
'           :        en faisant usage de l'ellispoide WGS-84
'           :        inspiré d'un script javascript
'           :        http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
'           :        http://www.movable-type.co.uk/
'---------------------------------------------------------------------------------------
Public Function distVincenty(p1 As Pt_GEOGRAPHIQUE, p2 As Pt_GEOGRAPHIQUE) As Double

    On Error GoTo Error_distVincenty

    Dim sMot As String
    Dim sTemp As String
    Dim sFile As String
    Dim cP As Integer
    Dim i As Integer
    Dim k As Integer
    Dim fA As Byte
    Dim fB As Byte
    Dim aT As Double
    Dim bT As Double
    Dim f As Double
    Dim A As Double
    Dim B As Double
    Dim C As Double
    Dim cos2SigmaM As Double
    Dim cosSigma As Double
    Dim cosSqAlpha As Double
    Dim deg As Double
    Dim deltaSigma As Double
    Dim iterLimit As Double
    Dim L As Double
    Dim lambda As Double
    Dim s As Double
    Dim sigma As Double
    Dim sinAlpha As Double
    Dim sinLambda As Double
    Dim sinSigma As Double
    Dim sinU2 As Double
    Dim U1 As Double
    Dim U2 As Double
    Dim uSq As Double
    Dim sinU1 As Double
    Dim cosU1 As Double
    Dim cosU2 As Double
    Dim lambdaP As Double
    Dim cosLambda As Double
    Dim rA As Double
    Dim rB As Double
    Dim rC As Double
    Dim rD As Double
    Dim rE As Double
    Dim dTempA As Double
    Dim dTempB As Double
    Dim dTempC As Double
    Dim oldEcart As Double



    aT = 6378137
    bT = 6356752.3142
    f = 1 / 298.257223563 ' WGS-84 ellipsoid
    L = (1000000 * p2.Lon - 1000000 * p1.Lon) / 1000000
    U1 = Arctan((1 - f) * Tan(p1.Lat))
    U2 = Arctan((1 - f) * Tan(p2.Lat))
    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)
        oldEcart = Abs(lambda - lambdaP)


        DoEvents
        iterLimit = iterLimit - 1
        sinLambda = sin(lambda)
        cosLambda = cos(lambda)
        sinSigma = Sqr((cosU2 * sinLambda) * (cosU2 * sinLambda) + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda))
        If (sinSigma = 0) Then
            distVincenty = 0 ' co-incident points
            Exit Function
        End If
        cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
        sigma = arctan2(sinSigma, cosSigma)
        DoEvents
        sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
        cosSqAlpha = 1 - sinAlpha * sinAlpha

        DoEvents

        If (cosSqAlpha = 0) Then
            distVincenty = Abs(aT * L) ' two points on equator
        End If

        cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha

        DoEvents
        C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
        lambdaP = lambda





        dTempA = (-1 + 2 * cos2SigmaM * cos2SigmaM)
        dTempB = (cos2SigmaM + (C * cosSigma) * dTempA)
        dTempC = (1 - C) * (f * sinAlpha)
        lambda = L + dTempC * (sigma + C * sinSigma * dTempB)




        DoEvents

        rA = (1 - C) * (f * sinAlpha)
        rB = (sigma + C * sinSigma)
        rC = (cos2SigmaM + (C * cosSigma) * (-1 + 2 * cos2SigmaM * cos2SigmaM))
        lambdaP = L + rA * rB * rC

        Debug.Print iterLimit & " " & Abs(lambda - lambdaP) * 100000#
        DoEvents
    Wend
    If (iterLimit = 0) And ((Abs(lambda - lambdaP) - oldEcart) / oldEcart) > 0.01 Then
        distVincenty = -1  ' formula failed to converge
        Debug.Assert False
        Exit Function
    End If
    uSq = cosSqAlpha * (aT * aT - bT * bT) / (bT * bT)


    'A-------------
    rA = (-768 + uSq * (320 - 175 * uSq))
    rB = (4096 + uSq * rA)
    A = 1 + uSq / 16384 * rB



    'B-------------
    B = 0
    rA = (-128 + uSq * (74 - 47 * uSq))
    rB = (256 + uSq * rA)
    B = uSq / 1024 * rB


    'deltaSigma
    rA = (-3 + 4 * cos2SigmaM * cos2SigmaM)
    rB = (-3 + 4 * sinSigma * sinSigma)
    rC = ((6 * cos2SigmaM) * rA)
    rD = (-1 + 2 * cos2SigmaM * cos2SigmaM)
    rE = (4 * (cosSigma * rD - B) / rC)
    deltaSigma = B * sinSigma * (cos2SigmaM + B / rE)

    '----------
    s = bT * A * (sigma - deltaSigma)

    distVincenty = dblRoundOff(s, 3) ' round to 1mm precision


Sortir_distVincenty:
    On Error GoTo 0
    Exit Function

Error_distVincenty:


End Function

Sub Main()
    Dim dDistance As Double
    Dim p2 As Pt_GEOGRAPHIQUE
    Dim p1 As Pt_GEOGRAPHIQUE
    
    p1.LatDegDD = 52.874
    p1.LonDegDD = 4.389
    
    p2.LatDegDD = 45.001
    p2.LonDegDD = 15.716
    
    
    p1.Lat = dtR(p1.LatDegDD)
    p1.Lon = dtR(p1.LonDegDD)
    
    p2.Lat = dtR(p2.LatDegDD)
    p2.Lon = dtR(p2.LonDegDD)
    
    
    dDistance = distVincenty(p1, p2)
    Debug.Assert False
End Sub