Geografische Distanz nach Thaddeus Vincenty

[Julian]

Ich bastel gerade ein GPS Programm und da ist es immer gut die Distanz zwischen zwei Punkten zu wissen. Ich habe den Algorithmus von Thaddeus Vincenty mal nach GLBasic portiert. Ein Testergebnis von München nach Berlin war bei mir 9mm größer als das vom JavaScript (oder meinem alten Python Script). Aber ich denke bei ~500km ist eine Abweichung von 9mm vollkommen i.O.

Code: (glbasic) [Select]

FUNCTION distance#: lat1, lon1, lat2, lon2
//source: http://www.movable-type.co.uk/scripts/latlong-vincenty.html
    //WGS-84 ellipsiod
    LOCAL a1 = 6378137.0, b = 6356752.3142, f = 1/298.257223563
   
    LOCAL L = (lon2-lon1);
    LOCAL U1 = ATAN((1-f)*TAN(lat1),1)
    LOCAL U2 = ATAN((1-f)*TAN(lat2),1)
    LOCAL sinU1 = SIN(U1)
    LOCAL sinU2 = SIN(U2)
    LOCAL cosU1 = COS(U1)
    LOCAL cosU2 = COS(U2)
    LOCAL lambda = L
    LOCAL iterLimit = 100
    LOCAL lambdaP = 0
   
    LOCAL cosSqAlpha
    LOCAL sinSigma
    LOCAL cos2SigmaM
    LOCAL cosSigma
    LOCAL sigma
   
    WHILE (ABS(lambda-lambdaP) > 1e-12 AND iterLimit > 0)
    LOCAL sinLambda = SIN(lambda)
        LOCAL cosLambda = COS(lambda)
        sinSigma = SQR((cosU2 * sinLambda) * (cosU2 * sinLambda) + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda))
        IF sinSigma = 0 THEN RETURN 0
        cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
        sigma = ATAN(sinSigma, cosSigma)
        LOCAL sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
        cosSqAlpha = 1 - sinAlpha * sinAlpha
        cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha
        IF NOT IsNumeric(cos2SigmaM) THEN cos2SigmaM = 0
        LOCAL 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*cos2SigmaM)))
    WEND
    IF iterLimit = 0 THEN RETURN -1
    LOCAL uSq = cosSqAlpha * (a1*a1 -     b*b) / (b*b)
    LOCAL A2 = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))
    LOCAL B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)))
    LOCAL deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)))
    LOCAL s = b*A2*(deg2rad(sigma)-deltaSigma)
    RETURN s
ENDFUNCTION

IMPORT "C" double __stdcall atof(char*);

FUNCTION IsNumeric%: string$

INLINE
double t ;
t = atof( string_Str ) ;
if (t == (double) 0.0) // retval if atof() can't convert
return (DGNat) FALSE ;
return (DGNat) TRUE ;
ENDINLINE

ENDFUNCTION

FUNCTION deg2rad#: deg
LOCAL pi = 3.14159265358979
RETURN deg*pi/180
ENDFUNCTIONEvtl. kanns ja der Ein oder Andere gebrauchen.

Gruß Julian

[Kitty Hello]

Krass! Berücksichtigt das Erdradius und Höhen, oder?

[Julian]

Die Genauigkeit der Formel liegt bei 0.5mm und orientiert sich am WGS 84 Elipsoid welcher die Erde beschreibt. Allerdings rechnet GLBasic nicht so genau (der gleiche Code in Python ist genauer). Wie gesagt ich hatte knapp 1cm Abweichung bei ~500km Distanz (das sind 2x10^-6% Abweichung). Das reicht für meine Zwecke vollkommen.

Gruß Julian

[Slydog]

Die Erde ist nicht flach? ha
Faszinierende Funktion!

Könnte Ihr Genauigkeit Probleme durch harte Werte verursacht werden, in eine Ganzzahl umgewandelt und nicht eine Gleitkommazahl wie kann erwarten gewesen wäre? Dies würde unerwünschte Rundungen.

ZB:

Code: (glbasic) [Select]

LOCAL B = uSq/1024 * (256 + USQ * (-128 + USQ * (* 74-47 USQ)))
Die 1024 würde ein Integer sein, und Sie würden das ".0" hinzufügen, um eine float sein.

Code: (glbasic) [Select]

LOCAL B = uSq/1024.0 * (256,0 + USQ * (-128,0 + USQ * (74,0 bis 47,0 * USQ)))
Aber das kann nicht dein Fall sein, wenn die restlichen Formel noch Schwimmer ist.
Ich bin nicht sicher ohne Prüfung.

[Original ENGLISH]
The Earth isn't flat? ha
Fascinating function!

Could your accuracy issues be caused by hard values being converted to an integer and not a floating point number as may have been expected?  This would cause undesired rounding.

Eg:
LOCAL B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)))

The 1024 would be an integer, and you would add the ".0" to be a float.
LOCAL B = uSq/1024.0 * (256.0+uSq*(-128.0+uSq*(74.0-47.0*uSq)))

But this may not be your case if the remaining formula is still float.
I'm not sure without testing.

[Kitty Hello]

a*b or a+b <- if any of these if floating point, the result is, too.