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. ;)
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
ENDFUNCTION
Evtl. kanns ja der Ein oder Andere gebrauchen.
Gruß Julian
Krass! Berücksichtigt das Erdradius und Höhen, oder?
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
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:
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.
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.
a*b or a+b <- if any of these if floating point, the result is, too.