_Geo_Diff - Exakte Berechnung der Distanz zweier Punkte auf der Erdoberfläche

  • Titel sagt viel, das Skript den Rest ^^ . Beispiel ist mit dabei.

    Spoiler anzeigen
    [autoit]


    ; Paris Observatory
    $geoParis = 'E-2-20-14-N-48-50-11'

    [/autoit] [autoit][/autoit] [autoit]

    ; US Naval Observatory
    $geoUSA = 'W-77-3-55.5-N-38-55-17'

    [/autoit] [autoit][/autoit] [autoit]

    ConsoleWrite(@LF)
    ConsoleWrite("Distance between the two given locations is " & _Geo_Diff($geoParis, $geoUSA) & " ± 0.05 km!" & @LF)
    ConsoleWrite(@LF)

    [/autoit] [autoit][/autoit] [autoit]

    ; #FUNCTION# ================================================================================
    ; Name ..........: _Geo_Diff
    ; Description ...: Compute the gedesic distance between two Earth surface
    ; ...............: coordinates to an accuracy of about ±50 meters. To get
    ; ...............: this degree of accuracy, this function takes into account
    ; ...............: the spheroidal flattening factor of the Earth rather than
    ; ...............: assuming that the Earth is a perfect sphere.
    ; Syntax ........: _Geo_Diff($geoStart, $geoFinish[, $iDec = 2[, $iUnit = 0]])
    ; Parameters ....: $geoStart - Location Coordinates in GeoData form (see Example)
    ; $geoFinish - Location Coordinates in GeoData form (see Example)
    ; $iDec - [optional] Number of decimal places. Default is 2.
    ; $iUnit - [optional] The unit to use. Default is 0.
    ; 0 = km (± 0.05)
    ; 1 = mi (± 0.03)
    ; 2 = nmi (± 0.03)
    ; Return values .: Distance in the given unit, rounded to the given number of decimal places.
    ; Author ........: minx
    ; Example .......: s. above
    ; ===========================================================================================
    Func _Geo_Diff($geoStart, $geoFinish, $iDec = 2, $iUnit = 0)
    Local $aData[4]
    Local $fTemp, $aTemp

    [/autoit] [autoit][/autoit] [autoit]

    $aTemp = StringSplit($geoStart, "-", 3)
    $aData[0] = ($aTemp[1] * 3600 + $aTemp[2] * 60 + $aTemp[3]) / 3600
    If $aTemp[0] = "E" Then $aData[0] *= -1
    $aData[1] = ($aTemp[5] * 3600 + $aTemp[6] * 60 + $aTemp[7]) / 3600
    If $aTemp[0] = "S" Then $aData[0] *= -1

    [/autoit] [autoit][/autoit] [autoit]

    $aTemp = StringSplit($geoFinish, "-", 3)
    $aData[2] = ($aTemp[1] * 3600 + $aTemp[2] * 60 + $aTemp[3]) / 3600
    If $aTemp[0] = "E" Then $aData[0] *= -1
    $aData[3] = ($aTemp[5] * 3600 + $aTemp[6] * 60 + $aTemp[7]) / 3600
    If $aTemp[0] = "S" Then $aData[0] *= -1

    [/autoit] [autoit][/autoit] [autoit]

    Local $ff = (1 / 298.257)

    [/autoit] [autoit][/autoit] [autoit]

    $UF = 1
    If $iUnit = 1 Then $UF = 1.609344
    If $iUnit = 2 Then $UF = 1.852

    [/autoit] [autoit][/autoit] [autoit]

    ; auxiliary angles
    $F = ($aData[1] + $aData[3]) / 2
    $G = ($aData[1] - $aData[3]) / 2
    $L = ($aData[0] - $aData[2]) / 2

    [/autoit] [autoit][/autoit] [autoit]

    ; (co)sines of auxliliary angles
    $SG = Sin($G * ATan(1) / 45)
    $CG = Cos($G * ATan(1) / 45)
    $SF = Sin($F * ATan(1) / 45)
    $CF = Cos($F * ATan(1) / 45)
    $SL = Sin($L * ATan(1) / 45)
    $CL = Cos($L * ATan(1) / 45)
    $S = ($SG * $CL) ^ 2 + ($CF * $SL) ^ 2
    $C = ($CG * $CL) ^ 2 + ($SF * $SL) ^ 2
    $O = ATan(Sqrt($S / $C))
    $R = Sqrt($S * $C) / $O
    $D = 2 * $O * 6378.14
    $H1 = (3 * $R - 1) / (2 * $C)
    $H2 = (3 * $R + 1) / (2 * $S)

    [/autoit] [autoit][/autoit] [autoit]

    ; distance and angle between start and finish
    ; on the surface of a geodesic sphereoid
    Return Round($D * ((($SF * $CG) ^ 2 * $H1 * $ff + 1) - (($CF * $SG) ^ 2 * $H2 * $ff)) / $UF, $iDec)
    EndFunc

    [/autoit] [autoit][/autoit] [autoit][/autoit]
  • Interessante Lösung.
    Wo stammt dieser Näherungsalgorithmus her und vor allem die Angabe, dass es immer im Bereich ± 0.05 km liegt?
    Für dein Beispiel betrug die Differenz zur exakten Lösung auf dem angegeben Ellipsoid übrigens 6m was ja im Bereich liegt.

    Vielleicht sollte man aber dazusagen: Bei Koordinatenberechnungen auf dem Ellipsoid müssen die angegebenen Koordinaten in einem Bezugssytem ermittelt wurden sein welchem das verwendete Ellipsoid zu Grunde liegt. Sonst sind die ermittelten Werte verfälscht. In Wikipedia z.B. sind die Angaben meist auf WGS84 bezogen (Referenzsystem für GPS-Koordinaten), also muss man auch dessen Werte für die Halbachse und Abplattung verwenden wenn man mit deren Koordinaten rechnen möchte.

    Theoretisch kann man auf dem Ellipsoid bis zu jeder beliebigen Genauigkeit rechnen. Da es sich hier aber um elliptische Integrale handelt welche ein Randwertproblem darstellen, ist die Berechnung ziemlich tricky und nur iterativ zu lösen.
    Hier mal eine Version welche auf etwa 1nm genau rechnet:

    Strecke zweier Punkte auf einem Rotationsellipsoid mit hoher Genauigkeit
    [autoit]

    #Region Hilfsgrößen und Strukturdefinitionen
    Global Const $PI = ACos(-1)
    Global Const $TAG_GEOCOORD = "double B;double L" ; Koordinatenangabe in Dezimalwinkeln (B=Breite/Latitude, L=Länge/Longitude)
    Global Const $TAG_GEOCOORD_HOUR = "int bD;int bM;double bS;int lD;int lM;double lS" ; Koordinateenangabe der Form: Winkel:Minute:Sekunde (b=Breite/Latitude, l=Länge/Longitude, D=Degrees, M=Minuten, S=Sekunden)
    #EndRegion Hilfsgrößen und Strukturdefinitionen

    [/autoit] [autoit][/autoit] [autoit]

    #Region Eingabedaten
    ; Paris Observatory
    $geoParis = 'E-2-20-14-N-48-50-11'
    ConsoleWrite("Paris: " & @TAB & _Coord_PrintHour(_Coords_Degree2Hour(_Coords_GeoData2Degrees($geoParis))) & @CRLF)
    ; US Naval Observatory
    $geoUSA = 'W-77-3-55.5-N-38-55-17'
    ConsoleWrite("USA: " & @TAB & _Coord_PrintHour(_Coords_Degree2Hour(_Coords_GeoData2Degrees($geoUSA))) & @CRLF)
    #EndRegion Eingabedaten

    [/autoit] [autoit][/autoit] [autoit]

    ; Berechne die Distanz der beiden Punkte auf dem Ellipsoid
    Global $a_Ret = _Ellipsoid_DistanceOf2Points(_Coords_GeoData2Degrees($geoParis), _Coords_GeoData2Degrees($geoUSA), 6378137, 1 / 298.257223563)
    ConsoleWrite(StringFormat("Strecke: %15.9f m\n", $a_Ret[0]))
    ConsoleWrite(StringFormat("Azimut im Punkt 1: %15.9f°\n", $a_Ret[1] * 180 / $PI))
    ConsoleWrite(StringFormat("Azimut im Punkt 2: %15.9f°\n", $a_Ret[2] * 180 / $PI))

    [/autoit] [autoit][/autoit] [autoit]

    ; #FUNCTION# ======================================================================================
    ; Name ..........: _Ellipsoid_DistanceOf2Points()
    ; Description ...: Berechnet die Entfernung zweier Punkte auf einem beliebigen Rotationsellipsoid
    ; Syntax ........: _Ellipsoid_DistanceOf2Points($P1, $P2[, $a = 6378137[, $f = 1 / 298.257222101]])
    ; Parameters ....: $P1 - Koordinaten von P1 als DezimalGrad-Angaben in der Struktur $TAG_GEOCOORD
    ; $P2 - Koordinaten von P2 als DezimalGrad-Angaben in der Struktur $TAG_GEOCOORD
    ; $a - [optional] Große Halbachse des Rotationsellipsoiden (weitere: http://de.wikipedia.org/wiki/Referenze…erenzellipsoide) (default:6378137)
    ; $f - [optional] Abplattungsfaktor des Rotationsellipsoiden (weitere: http://de.wikipedia.org/wiki/Referenze…erenzellipsoide) (default:1 / 298.257222101)
    ; $EPS- zu erreichende Genauigkeit
    ; Return values .: Success Array: [Strecke, Azimut1, Azimut2]
    ; Failure
    ; Author ........: AspirinJunkie
    ; Related .......: __Gn, __Mn, __Simp, __Sn, __Tn, CalcIntegral, CalcIntegralAdaptiv, f, IntMonteCarlo, Sign
    ; Remarks .......: Lösungsalgorithmus entsprechend H. Schmidt, RWTH Aachen 2000
    ; =================================================================================================
    Func _Ellipsoid_DistanceOf2Points($P1, $P2, $a = 6378137, $f = 1 / 298.257223563, Const $EPS = 1e-15)
    Local $s_Funktion ; Variable welche die zu integrierende Funktion f(x) als String enthält
    Local $b = $a - $f * $a ; kleine Halbachse des Ellipsoides
    Local $e2 = 1 - (($b / $a) ^ 2) ; Quadrat der 1. numerischen Exzentrizität

    [/autoit] [autoit][/autoit] [autoit]

    Local $B1 = $P1.B * $PI / 180 ; Breite von Punkt 1 in Radiant
    Local $B2 = $P2.B * $PI / 180 ; Breite von Punkt 2 in Radiant
    Local $L1 = $P1.L * $PI / 180 ; Länge von Punkt 1 in Radiant
    Local $L2 = $P2.L * $PI / 180 ; Länge von Punkt 2 in Radiant

    [/autoit] [autoit][/autoit] [autoit]

    ; reduzierte Breiten beta1 und beta2:
    Local $beta1 = (Abs($B1 - ($PI / 2)) < 1E-5) ? $B1 : ATan(Sqrt(1 - $e2) * Tan($B1)) ; zusätzlicher Check ob eine Breite dem Pol entspricht
    Local $beta2 = (Abs($B2 - ($PI / 2)) < 1E-5) ? $B2 : ATan(Sqrt(1 - $e2) * Tan($B2))

    [/autoit] [autoit][/autoit] [autoit]

    ; Ausgangswerte für Restabweichung und Langendifferenz
    Local $deltaL
    Local $deltaLsoll = $L2 - $L1
    Local $deltaLs = $L2 - $L1
    Local $dL = $deltaLsoll
    Local $dL1, $dL2 ; Längendifferenzen für beide Fälle

    [/autoit] [autoit][/autoit] [autoit]

    ; Vorzeichen der Längendifferenz
    Local $y = Sign($deltaLsoll)

    [/autoit] [autoit][/autoit] [autoit]

    ; iterative Berechnung der Intergrationsgrenzen
    Local $tanBetaMax ; tangens der maximalen Breite der geodätischen Linie
    Local $h2, $h ; Clairaut-Konstante
    Local $w1, $w2, $wu, $wo ; Integrationsgrenzen
    Local $deltaL1, $deltaL2 ; Lösung für Fall 1 und Fall 2
    Local $FALL ; Angabe welcher Fall vorherrscht
    Local $d_ItMax = 50, $it = 0 ; Maximale Iterationsanzahl (um Endlosschleifen zu vermeiden)
    Do
    ; tangens der maximalen Breite der geodätischen Linie
    $tanBetaMax = Sqrt(1 - $e2) * (Sqrt(Tan($B1) * Tan($B1) + Tan($B2) * Tan($B2) - 2 * Tan($B1) * Tan($B2) * Cos($deltaLs)) / Sin($deltaLs))

    [/autoit] [autoit][/autoit] [autoit]

    ; Clairaut-Konstante
    $h2 = 1 / (1 + $tanBetaMax * $tanBetaMax)

    [/autoit] [autoit][/autoit] [autoit]

    ; Integrationsgrenzen
    $w1 = ASin(Tan($beta1) / $tanBetaMax)
    $w2 = ASin(Tan($beta2) / $tanBetaMax)

    [/autoit] [autoit][/autoit] [autoit]

    If Abs($dL) < $EPS Or $it >= $d_ItMax Then ExitLoop ; Integrationsgrenzen hinreichend genau bestimmt

    [/autoit] [autoit][/autoit] [autoit]

    $wu = (Abs($w1) < Abs($w2) ? Abs($w1) : Abs($w2)) ; = Min(w1,w2)
    $wo = (Abs($w1) > Abs($w2) ? Abs($w1) : Abs($w2)) ; = Max(w1,w2)

    [/autoit] [autoit][/autoit] [autoit]

    ; Berechnung von wu und wo für Fall 1 (es gibt zwei mögliche Lösungen - daher beide berechnen und mit $deltaLsoll vergleichen)
    $FALL = 1
    $s_Funktion = "sqrt( 1-(" & $e2 * $h2 & " / (" & $h2 & " + " & (1 - $h2) & " * sin(x)*sin(x))))"

    [/autoit] [autoit][/autoit] [autoit]

    $deltaL1 = $y * Abs(_MathInt_Romberg($s_Funktion, $wu, $wo, $EPS)) ; Integralberechnung bis auf 1e-14 Genauigkeit
    ; Addition der Abschnitte symmetrisch zum Äquator
    If ($w1 * $w2) < 0 Then $deltaL1 += $y * 2 * _MathInt_Romberg($s_Funktion, 0, $wu, $EPS)

    [/autoit] [autoit][/autoit] [autoit]

    ; Längendifferenz $dL und Restabweichung DL für Fall I
    $deltaL = $deltaL1
    $dL1 = $deltaLsoll - $deltaL1
    $dL = $dL1

    [/autoit] [autoit][/autoit] [autoit]

    ; Addition der Abschnitte symmetrisch zu Bmax für Fall II
    $deltaL2 = $deltaL1 + $y * 2 * _MathInt_Romberg($s_Funktion, $wo, 0.5 * $PI, $EPS)

    [/autoit] [autoit][/autoit] [autoit]

    ; Längendifferenz $dL und Restabweichung DL für Fall II
    $dL2 = $deltaLsoll - $deltaL2

    [/autoit] [autoit][/autoit] [autoit]

    ; Entscheidung für Fall 1 oder 2
    If (Abs($dL2) < Abs($dL1)) Then ; Fall 2 ist vorhanden - daher dessen Ergebnisse übernehmen
    $FALL = 2
    $deltaL = $deltaL2
    $dL = $dL2
    EndIf

    [/autoit] [autoit][/autoit] [autoit]

    ; Korrektur von $deltaLs zur Berechnung von tanBmax
    $deltaLs = $deltaLs + (1 + ($dL / ($FALL * $deltaL))) * $dL
    $it += 1 ; Iterationszähler inkrementieren
    Until 0
    ; Azimutberechnung
    ; Vorzeichen der Breitendifferenz
    Local $x = Sign($w2 - $w1)
    If $FALL = 2 Then $x = Sign(Sign($B1 + $B2) * $PI - $w2 - $w1)

    [/autoit] [autoit][/autoit] [autoit]

    ; Azimute A1 und A2
    Local $A1 = (1 - $x) * 0.5 * $PI + $x * $y * ASin(Sqrt($h2) / Cos($beta1))
    Local $A2s = (1 - $x) * 0.5 * $PI + $x * $y * ASin(Sqrt($h2) / Cos($beta2))
    If $FALL = 2 Then $A2s = $PI - $A2s
    If $A1 < 0 Then $A1 += 2 * $PI
    If $A2s < 0 Then $A2s += 2 * $PI
    Local $A2 = $A2s + $PI
    If $A2 > $PI Then $A2 -= 2 * $PI

    [/autoit] [autoit][/autoit] [autoit]

    ; Berechnung der Streckenlänge
    ; Integrationsgrenzen v1 und v2 bestimmen
    Local $v1 = ASin(Sin($beta1) / Sqrt(1 - $h2))
    Local $v2 = ASin(Sin($beta2) / Sqrt(1 - $h2))

    [/autoit] [autoit][/autoit] [autoit]

    If (Abs($deltaLsoll - $PI) < 1E-10) Then $FALL = 2 ; check ob deltaLsoll = Pi
    If $FALL = 2 Then $v2 = Sign($B1 + $B2) * $PI - $v2

    [/autoit] [autoit][/autoit] [autoit]

    ; Berechnung der Streckenlänge durch numerische Intergration
    $s_Funktion = "sqrt(1-" & $e2 & "*(1-sin(x)*sin(x)*" & 1 - $h2 & "))"
    Local $S = $a * _MathInt_Romberg($s_Funktion, $v1, $v2, $EPS) ; Strecke zwischen beiden Punkten in m

    [/autoit] [autoit][/autoit] [autoit]

    Local $a_Ret = [$S, $A1, $A2]
    Return $a_Ret
    EndFunc ;==>_Ellipsoid_DistanceOf2Points

    [/autoit] [autoit][/autoit] [autoit][/autoit] [autoit]

    #Region Hilsfunktionen
    ; Signum (Vorzeichen) einer Zahl
    Func Sign($x)
    Return ($x < 0) ? -1 : ($x > 0) ? 1 : 0
    EndFunc ;==>Sign

    [/autoit] [autoit][/autoit] [autoit]

    ; #=============================================
    ; = _Coords_GeoData2Degrees =
    ; =============================================*
    ; Wandelt Koordinaten im Geodata-Format (z.B. E-2-20-14-N-48-50-11) - in normale Dezimalwinkel
    Func _Coords_GeoData2Degrees($s_GeoDataForm)
    ; by AspirinJunkie
    Local $aTemp = StringSplit($s_GeoDataForm, "-", 3)
    Local $t_Ret = DllStructCreate($TAG_GEOCOORD)
    $t_Ret.L = ($aTemp[1] * 3600 + $aTemp[2] * 60 + $aTemp[3]) / 3600 * (($aTemp[0] = "W") ? -1 : 1)
    $t_Ret.B = ($aTemp[5] * 3600 + $aTemp[6] * 60 + $aTemp[7]) / 3600 * (($aTemp[4] = "S") ? -1 : 1)
    Return $t_Ret
    EndFunc ;==>_Coords_GeoData2Degrees

    [/autoit] [autoit][/autoit] [autoit]

    ; #=============================================
    ; = _Coords_Degree2Hour =
    ; =============================================*
    ; Wandelt Koordinaten in Dezimalform in Winkel:Minute:Sekunde-Form (Bogenminuten)
    Func _Coords_Degree2Hour($t_Degree)
    ; by AspirinJunkie
    Local $t_Ret = DllStructCreate($TAG_GEOCOORD_HOUR)

    [/autoit] [autoit][/autoit] [autoit]

    $t_Ret.bD = Int($t_Degree.B)
    Local $temp = $t_Degree.B - $t_Ret.bD
    $t_Ret.bM = Int(60 * $temp)
    $t_Ret.bS = Abs(3600 * $temp - 60 * $t_Ret.bM)
    $t_Ret.bM = Abs($t_Ret.bM)

    [/autoit] [autoit][/autoit] [autoit]

    $t_Ret.lD = Int($t_Degree.L)
    $temp = $t_Degree.L - $t_Ret.lD
    $t_Ret.lM = Int(60 * $temp)
    $t_Ret.lS = Abs(3600 * $temp - 60 * $t_Ret.lM)
    $t_Ret.lM = Abs($t_Ret.lM)
    Return $t_Ret
    EndFunc ;==>_Coords_Degree2Hour

    [/autoit] [autoit][/autoit] [autoit]

    ; Gibt die Koordinaten im Bogenminuten-Format als String zurück
    Func _Coord_PrintHour($t_Hour)
    ; by AspirinJunkie
    Return StringFormat("Lat: %2d° %2d' %2d"".%05d %1s\tLon: %2d° %2d' %2d"".%05d %1s", _
    Abs(Int($t_Hour.bD)), $t_Hour.bM, Int($t_Hour.bS), Int(StringLeft(StringTrimLeft($t_Hour.bS, StringLen(Int($t_Hour.bS)) + 1), 5)), ($t_Hour.bD >= 0) ? "N" : "O", _
    Abs(Int($t_Hour.lD)), $t_Hour.lM, Int($t_Hour.lS), Int(StringLeft(StringTrimLeft($t_Hour.lS, StringLen(Int($t_Hour.lS)) + 1), 5)), ($t_Hour.lD >= 0) ? "E" : "W")
    EndFunc ;==>_Coord_PrintHour
    #EndRegion Hilsfunktionen

    [/autoit] [autoit][/autoit] [autoit][/autoit] [autoit][/autoit] [autoit][/autoit] [autoit]

    #Region Funktionen zur Berechnung einer numerischen Integration
    ; #FUNCTION#=================================================================================
    ; Name...........: _MathInt_Romberg
    ; Description ...: Berechnet das bestimmte Integral eindimensionalen Funktion mit
    ; einer festzulegenden relativen Genauigkeit
    ; Syntax.........: _MathInt_Romberg(Const $a, Const $b, Const $EPS = 0.0001)
    ; Parameters ....: $cb_F - Funktionsvariable oder Funktionsstring einer 1D-Funktion f(x)
    ; $a - Untere Integrationsgrenze
    ; $b - obere Integrationsgrenze
    ; $EPS - zu erreichender relativer Integrationsfehler
    ; Hinweis: != absoluter Fehler, rel. Fehler = Abs(1-(berechneter Integralwert / wahrer Integralwert)
    ; $kMax - Maximale Anzahl an Halbierungsintervallen
    ; Return values .: Integral in den Grenzen $a-$b, @extended = Anzahl benötigter Iterationen
    ; @error = 1: Maximale Iterationsanzahl erreicht
    ; Remarks .......: Romberg-Integrationsalgorithmus
    ; Author ........: AspirinJunkie
    ;============================================================================================
    Func _MathInt_Romberg(Const $cb_F, Const $a, Const $b, Const $EPS = 0.0001, Const $kMax = 30)
    Local $a_R1[$kMax + 1], $a_R2[$kMax + 1]
    Local $h = $b - $a, $n = 1, $sumf, $f

    [/autoit] [autoit][/autoit] [autoit]

    $a_R1[0] = 0.5 * $h * ((IsFunc($cb_F) ? $cb_F($a) : __F($cb_F, StringFormat("%.20f", $a))) + (IsFunc($cb_F) ? $cb_F($b) : __F($cb_F, StringFormat("%.20f", $b)))) ; Erste Approximation

    [/autoit] [autoit][/autoit] [autoit]

    For $k = 1 To $kMax ; Schrittweise halbieren
    $sumf = 0
    For $i = 1 To $n
    $sumf += (IsFunc($cb_F) ? $cb_F($a + ($i - 0.5) * $h) : __F($cb_F, StringFormat("%.20f", $a + ($i - 0.5) * $h)))
    Next
    $a_R2[0] = 0.5 * ($a_R1[0] + $h * $sumf) ; Trapezoid-Formel
    $f = 1e0
    For $j = 1 To $k ; Quadratur-Ordnung erhöhen
    $f *= 4
    $a_R2[$j] = ($f * $a_R2[$j - 1] - $a_R1[$j - 1]) / ($f - 1) ; Neue Approximation
    Next
    If $k > 1 Then ; Check auf Konvergenz
    If (Abs($a_R2[$k] - $a_R1[$k - 1]) <= ($EPS * Abs($a_R2[$k]))) Then ExitLoop
    If (Abs($a_R2[$k]) <= $EPS) And (Abs($a_R2[$k]) <= Abs($a_R2[$k] - $a_R1[$k - 1])) Then ExitLoop
    EndIf
    $h *= 0.5 ; Halbiere Schrittweite
    $n *= 2
    For $j = 0 To $k
    $a_R1[$j] = $a_R2[$j]
    Next
    Next
    If $k >= $kMax Then Return SetError(1, $k, $a_R2[$k - 1])
    Return SetExtended($k, $a_R2[$k])
    EndFunc ;==>_MathInt_Romberg

    [/autoit] [autoit][/autoit] [autoit]

    ; Wrapperfunktion für den Fall, dass die zu integrierende Funktion als String übergeben wird
    Func __F($s_Func, $x)
    $s_Func = StringReplace($s_Func, "exp", "§")
    $s_Func = StringReplace($s_Func, "x", $x)
    $s_Func = StringReplace($s_Func, "§", "exp")
    Return Execute($s_Func)
    EndFunc ;==>__F
    #EndRegion Funktionen zur Berechnung einer numerischen Integration

    [/autoit]


    Ein Hinweis vielleicht noch: Bei deinem Koordinatenparsing negierst du die Koordinatenangaben wenn sie Ostkoordinaten sind. Normal ist es aber eigentlich Westkoordinaten zu negieren.
    Am Ergebnis ändert das aber nichts - eher ein Schönheitsdetail.

    Ein sehr schönes Online-Tool mit dem man sich das ganze berechnen und vor allem anzeigen lassen kann ist hier zu finden: http://geographiclib.sourceforge.net/scripts/geod-google.html (unser Anwendungsfall ist der "Inverse")

    Edit: Hab auch mal die Azimutwinkel im Anfangs und Endpunkt mit hinzugefügt, damit man weiß in welche Richtung man loslaufen muss und in welche man ankommt.
    Edit2: Auf Romberg-Integration umgestellt.
    Edit3: Link zum Online-Tool hinzugefügt

    7 Mal editiert, zuletzt von AspirinJunkie (11. September 2014 um 12:35)

  • Wo stammt dieser Näherungsalgorithmus her

    Laut meinem Gekritzel ist das eine Methode von Andoyer. Das stammt aus meinen Aufzeichnungen ausm Astronomie-Kurs. Wir hatten das damals in Excel umgesetzt, aber ich wollte es nochmal als AutoIt Funktion umschreiben.

    Ein Hinweis vielleicht noch: Bei deinem Koordinatenparsing negierst du die Koordinatenangaben wenn sie Ostkoordinaten sind. Normal ist es aber eigentlich Westkoordinaten zu negieren.
    Am Ergebnis ändert das aber nichts - eher ein Schönheitsdetail.

    Da habe ich genau das Gegenteil gelernt und damit auch in Prüfungen gerechnet ^^

    Theoretisch kann man auf dem Ellipsoid bis zu jeder beliebigen Genauigkeit rechnen.

    Das ist sicher richtig, aber diese Näherung reicht völlig aus, wenn man sie z.B. mit der geläufigen AGPS Triangulation vergleicht.

  • Da habe ich genau das Gegenteil gelernt und damit auch in Prüfungen gerechnet

    Na dann: Glück gehabt. Wenn du nicht gerade mal derartige Prüfungen schreibst, kannst du es ja von nun an richtig herum machen.

    Das ist sicher richtig, aber diese Näherung reicht völlig aus, wenn man sie z.B. mit der geläufigen AGPS Triangulation vergleicht.

    Das wollte ich Dir hiermit ja zum Ausdruck bringen. Wenn deine Koordinaten aus GPS-Messungen stammen müsstest du auch mit dem WGS84-Ellipsoid rechnen anstatt einem beliebigen.
    Für das Beispiel hier hat man schon allein dadurch einen Fehler von 3m der nicht sein müsste.
    Herkömmliche GPS-Messungen sind außerdem deutlich genauer als 50m und eine höhere Rechengenauigkeit für die Distanz daher auch hierbei durchaus wünschenswert.

  • Ja das ist erstmal bisschen verwirrend, aber richtig so.
    Kurz: Die Winkel müssen nicht zusammen 360° ergeben.
    Wenn man sich auf einer Kugel ein Dreieck zeichnet und die Kanten dabei die jeweils kürzesten Verbindungen auf der Kugel sind dann ist die Innenwinkelsumme in diesem Dreieck größer als 180°.
    Die Gesetzmäßigkeiten der planaren Geometrie gelten nicht mehr vollständig.

    Oder ist der WGS84 Geoid so kartoffelig?

    Nein nein es ist ein ganz glattes normales Rotationsellipsoid.
    Kartoffelig ist erst das Geoid.