Soyez le premier à donner votre avis sur cette source.
Snippet vu 35 386 fois - Téléchargée 28 fois
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
25 janv. 2011 à 10:13
16 avril 2008 à 12:31
A(xa,ya), B(xb,yb)
CapLoxo = Atn((xa - xb) / Log(Tan(ya / 2 + pi / 4) / Tan(yb / 2 + pi / 4)))
La distance Loxo ainsi :
DstLoxo = RayonSphere * Abs((ya - yb) / Cos(CapLoxo))
18 janv. 2008 à 14:20
18 janv. 2008 à 13:59
1.609km est la valeur d'un "miles" anglais. En revanche en navigation on utilise le "mile nautique" qui vaut 1.852 km.
29 déc. 2007 à 23:02
Vous n'êtes pas encore membre ?
inscrivez-vous, c'est gratuit et ça prend moins d'une minute !
Les membres obtiennent plus de réponses que les utilisateurs anonymes.
Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.
Le fait d'être membre vous permet d'avoir des options supplémentaires.