Calculate distance on Earth surface via Vincenty formula
Posted: Mon Dec 23, 2024 6:00 pm
Dear Pb programmers,
I am busy learning PB, and to learn the geometric functions of PB I rewrote an old Freebasic listing for PB. I don't want to tell you that it wasn't simple. OK, you can say I'm a bit crazy, but to learn to swim you have to jump into the water. So...
To be honest, it took some sweating to get it working properly. I also realize now that the next step in my PB learning process will be to get the hang of all the screen output.
I am publishing my old Freebasic listing and its output.
Below that the Purebasic listing (and you can see the output on your screen...).
I have added as many comments as possible to make everything clear.
For my ley line research friends I will expand it later with routines like: "Find out from a number of geographical coordinates where three or more points of them lie on a straight line on the earth's surface. (This is of course the projection of those points on a flat surface, but this will take us too far off here.)
Curious about the replies of the experienced PB programmers and possible tips to improve the program even more.
My old freebasic listing:
Output from the freebasic listing:
I am busy learning PB, and to learn the geometric functions of PB I rewrote an old Freebasic listing for PB. I don't want to tell you that it wasn't simple. OK, you can say I'm a bit crazy, but to learn to swim you have to jump into the water. So...
To be honest, it took some sweating to get it working properly. I also realize now that the next step in my PB learning process will be to get the hang of all the screen output.
I am publishing my old Freebasic listing and its output.
Below that the Purebasic listing (and you can see the output on your screen...).
I have added as many comments as possible to make everything clear.
For my ley line research friends I will expand it later with routines like: "Find out from a number of geographical coordinates where three or more points of them lie on a straight line on the earth's surface. (This is of course the projection of those points on a flat surface, but this will take us too far off here.)
Curious about the replies of the experienced PB programmers and possible tips to improve the program even more.
My old freebasic listing:
Code: Select all
Declare Function distance (lat1 As Double, lon1 As Double, lat2 As Double, lon2 As Double) as Double
#define radian(x) ( (x) * 3.1415926535897932384626433 / 180 )
'Freebasic:
'#define allows to declare text-based preprocessor macros.
'Once the compiler has seen a #define,
'it will start replacing further occurrences of identifier with body, body may be empty.
'The expansion is done recursively,
'until there is nothing more to expand and the compiler can continue analyzing the resulting code.
'#undef can be used to make the compiler forget about a #define.
declare FUNCTION Round (x AS DOUBLE, Digit_Count AS INTEGER = 0) AS DOUBLE
Dim x As String
Print "De afstand tussen LC1 en K7 onder Retie is:" ; distance(51.28382, 5.01033, 51.26397, 5.07559) ; " kilometer"
Print "De afstand tussen LC1 en K8 onder Retie is:" ; distance(51.28382, 5.01033, 51.25774, 5.09963) ; " kilometer"
Print "De afstand tussen LC1 en K9 onder Retie is:" ; distance(51.28382, 5.01033, 51.25193, 5.12125) ; " kilometer"
Print
Print "De afstand tussen Kapel K7 en Kapel K8 onder Retie is:" ; distance(51.26397, 5.07559, 51.25774, 5.09963) ; " kilometer"
Print "De afstand tussen Kapel K8 en Kapel K9 onder Retie is:" ; distance(51.25774, 5.09963, 51.25193, 5.12125) ; " kilometer"
Print "De afstand tussen Kapel K7 en Kapel K9 onder Retie is:" ; distance(51.26397, 5.07559, 51.25193, 5.12125) ; " kilometer"
Print
Print
Print "De afstand tussen het Paleis op de Dam tot de evenaar is volgens de UTM coordinaten 5804,225 km."
Print "De afstand tussen het Paleis op de Dam tot de evenaar is .........................:" ; distance(52.373126, 4.892467, 0.00, 4.892467) ; " kilometer"
Print
Print "De afstand Amersfoort Lieve Vrouwentoren 31U tot het Paleis op de Dam 31U is :" ; distance(52.15517440, 5.38720621, 52.373126, 4.892467) ; " kilometer"
Print "De afstand Amersfoort Lieve Vrouwentoren 31U tot Amsterdam Westertoren 31U is :" ; distance(52.15517440, 5.38720621, 52.37453253, 4.88352559) ; " kilometer"
Print "Volgens Geoflifeline.com is de afstand tussen Lieve Vrouwetoren en Martinitoren : 142.83216 kilometer"
Print "De afstand Amersfoort Lieve Vrouwentoren 31U tot Groningen Martinitoren 32U is :" ; distance(52.15517440, 5.38720621, 53.21938317, 6.56820053) ; " kilometer"
Print
Print "De afstand van Parijs tot de evenaar is de volgens GeoLifeLine.com 5392.92947 km"
Print "De afstand van Parijs tot de evenaar is volgens deze berekening .:" ; distance(48.67, 2.33, 0, 2.33) ; " kilometer"
Print "De afstand van Parijs tot de noordpool is volgen GeoLifeLine 4609.03626 km"
Print "De afstand van Parijs tot de noordpool is volgens deze berekening:" ; distance(48.67, 2.33, 90, 2.33) ; " kilometer"
Print
Print "De afstand tussen Nashville International Airport (BNA), USA"
Print " en Los Angeles International Airport (LAX), USA"
Print "volgens Rosettacode tussen 2886.444 en 2887.259 km"
print "volgens GeoLifeLine ....... 2884.70058 km"
Print "volgens deze berekening...:" ; distance(36.12, -86.76, 33.94, -118.40) ; " kilometer"
Print
Print
Print "De afstand tussen Los Angeles International Airport (LAX), USA"
Print " (JFK), USA " & distance(33.95000, -118.40000, 40.63333,-73.78333) ; " kilometer"
Print "De afstand tussen Los Angeles International Airport (LAX), USA"
Print " en punt D is " & distance(33.95000, -118.40000, 34.50000,-116.50000) ; " kilometer"
Print
Print "Test hoeveel kilometer in een breedtegraad op latitude 51"
Print "De afstand voor 1 breedtegraad in Belgie is :" ; distance(51, 1, 51, 2) ; " kilometer"
Print "De afstand op de evenaar voor 1 breedtegraad :" ; distance(0, 1, 0, 2) ; " kilometer"
Print "4 maal 90 graden = lengte van de evenaar :" ; (4 * distance(0, 1, 0, 91)) ; " kilometer"
Print " Volgens WGS 84 voor gps is de evenaar 40075.0167 kilometer"
Input x
Function distance (lat1 As Double, lon1 As Double, lat2 As Double, lon2 As Double) as Double
' Geodetische afstand in kilometer tussen twee punten lat1 lon1 en lat2 lon2
' Coordinaten uitgedrukt in decimale gradianen.
' Wiskundig: Vincenty inverse formule voor een ellipsoide (afgeplatte aardbol).
'=================================================================================
' Code has been ported by lost_species from www.aliencoffee.co.uk to VBA from javascript published at:
' http://www.movable-type.co.uk/scripts/latlong-vincenty.html
' * from: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the
' * Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975
' * http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
'Constanten voor de WGS-84 ellipsiode:
Dim radian_equator As Double = 6378137
Dim radian_poles As Double = 6356752.314245
Dim f As Double = 1 / 298.257223563 'Flattening of the elipsoid, afplatting van de aarde
'Itereer tot λ verwaarloosbaar is (10^-12 ≈ 0.006mm)
'Dim ItereerTot As Double = 0.000000000001
'Met onderstaande waarde wordt slechts 30cm afwijking bekomen op de totale lengte van de evenaar
Dim ItereerTot As Double = 0.00000001
Dim L As Double
Dim U1 As Double
Dim U2 As Double
Dim sinU1 As Double
Dim sinU2 As Double
Dim cosU1 As Double
Dim cosU2 As Double
Dim Lambda As Double
Dim LambdaP As Double
Dim aantalIteraties As Integer
Dim sinLambda As Double
Dim cosLambda As Double
Dim sinSigma As Double
Dim cosSigma As Double
Dim sigma As Double
Dim sinAlpha As Double
Dim coskwadraatAlpha As Double
Dim cos2SigmaM As Double
Dim C As Double
Dim uKwadraat As Double
Dim upper_A As Double
Dim upper_B As Double
Dim deltaSigma As Double
Dim dist As Double
L = radian(lon2 - lon1) '= het longitude verschil
U1 = Atn((1 - f) * Tan(radian(lat1))) 'reduced latitude
U2 = Atn((1 - f) * Tan(radian(lat2)))
sinU1 = Sin(U1)
cosU1 = Cos(U1)
sinU2 = Sin(U2)
cosU2 = Cos(U2)
Lambda = L
LambdaP = 2 * 3.141592653589793238
'aantalIteraties = 10
aantalIteraties = 6
While (Abs(Lambda - LambdaP) > ItereerTot) And (aantalIteraties > 0)
aantalIteraties = aantalIteraties - 1
sinLambda = Sin(Lambda)
cosLambda = Cos(Lambda)
sinSigma = Sqr(((cosU2 * sinLambda) ^ 2) + ((cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ^ 2))
If sinSigma = 0 Then 'overlappende punten (Coïncident points)
Return 0
End If
cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
sigma = ATan2(sinSigma, cosSigma)
sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
coskwadraatAlpha = 1 - (sinAlpha ^ 2)
'Vermijd deling door nul:
If coskwadraatAlpha = 0 Then 'Twee punten op de evenaar.
cos2SigmaM = 0
Else
cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / coskwadraatAlpha
End If
C = f / 16 * coskwadraatAlpha * (4 + f * (4 - 3 * coskwadraatAlpha))
LambdaP = Lambda
Lambda = L + (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * (cos2SigmaM ^ 2))))
Wend
If aantalIteraties < 1 Then
Print "Maximum aantal iteraties is uitgevoerd. Oorzaken: de ingegeven data en of de grenzen van Vincenty formule."
Return 0
End If
uKwadraat = coskwadraatAlpha * (radian_equator ^ 2 - radian_poles ^ 2) / (radian_poles ^ 2)
upper_A = 1 + uKwadraat / 16384 * (4096 + uKwadraat * (-768 + uKwadraat * (320 - 175 * uKwadraat)))
upper_B = uKwadraat / 1024 * (256 + uKwadraat * (-128 + uKwadraat * (74 - 47 * uKwadraat)))
deltaSigma = upper_B * sinSigma * (cos2SigmaM + upper_B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM ^ 2) - upper_B / 6 * cos2SigmaM * (-3 + 4 * sinSigma ^ 2) * (-3 + 4 * cos2SigmaM ^ 2)))
'Voor het resultaat in meters:
'dist = (radian_poles * upper_A * (sigma - deltaSigma))
'Gedeeld door 1000 voor kilometers:
dist = (radian_poles * upper_A * (sigma - deltaSigma)) / 1000
dist=Round(dist,3)
Return dist
End Function
FUNCTION Round (x AS DOUBLE, Digit_Count AS INTEGER = 0) AS DOUBLE
CONST MaxRoundPlaces = 16
DIM AS DOUBLE ReturnAnswer, Dec_Place
DIM AS INTEGER N
IF (Digit_Count >= 0) AND (Digit_Count < MaxRoundPlaces) THEN
N = Digit_Count
ELSE
N = 0
END IF
Dec_Place = 10 ^ N
ReturnAnswer = INT((ABS (x) * Dec_Place) + 0.5)
Round = SGN(x) * ReturnAnswer / Dec_Place
END FUNCTION
My Purebasic listing:For compleetness I show the results of my old freebasic listing.
Sorry it is in Dutch.
Feel free to translate it in Your language via Google translate.
(The PB listing outputs in Englisch)
Belgium: (kapel = chapel, see also PB listing)
De afstand tussen LC1 en K7 onder Retie is: 5.061 kilometer
De afstand tussen LC1 en K8 onder Retie is: 6.874 kilometer
De afstand tussen LC1 en K9 onder Retie is: 8.516 kilometer
De afstand tussen Kapel K7 en Kapel K8 onder Retie is: 1.816 kilometer
De afstand tussen Kapel K8 en Kapel K9 onder Retie is: 1.642 kilometer
De afstand tussen Kapel K7 en Kapel K9 onder Retie is: 3.457 kilometer
The Netherlands:
De afstand tussen het Paleis op de Dam tot de evenaar is volgens de UTM coordinaten 5804,225 km.
De afstand tussen het Paleis op de Dam tot de evenaar is .........................: 5804.862 kilometer
De afstand Amersfoort Lieve Vrouwentoren 31U tot het Paleis op de Dam 31U is : 41.582 kilometer
De afstand Amersfoort Lieve Vrouwentoren 31U tot Amsterdam Westertoren 31U is : 42.169 kilometer
Volgens Geoflifeline.com is de afstand tussen Lieve Vrouwetoren en Martinitoren : 142.83216 kilometer
De afstand Amersfoort Lieve Vrouwentoren 31U tot Groningen Martinitoren 32U is : 142.832 kilometer
France:
De afstand van Parijs tot de evenaar is de volgens GeoLifeLine.com 5392.92947 km
De afstand van Parijs tot de evenaar is volgens deze berekening .: 5392.929 kilometer
De afstand van Parijs tot de noordpool is volgen GeoLifeLine 4609.03626 km
De afstand van Parijs tot de noordpool is volgens deze berekening: 4609.036 kilometer
USA:
De afstand tussen Nashville International Airport (BNA), USA
en Los Angeles International Airport (LAX), USA
volgens Rosettacode tussen 2886.444 en 2887.259 km
volgens GeoLifeLine ....... 2884.70058 km
volgens deze berekening...: 2884.701 kilometer
De afstand tussen Los Angeles International Airport (LAX), USA
(JFK), USA 3981.601 kilometer
De afstand tussen Los Angeles International Airport (LAX), USA
en punt D is 185.389 kilometer
Some other tests (see Englisch output of the Purebasis listing.)
Test hoeveel kilometer in een breedtegraad op latitude 51
De afstand voor 1 breedtegraad in Belgie is : 70.197 kilometer
De afstand op de evenaar voor 1 breedtegraad : 111.319 kilometer
4 maal 90 graden = lengte van de evenaar : 40075.016 kilometer
Volgens WGS 84 voor gps is de evenaar 40075.0167 kilometer
Code: Select all
Procedure.d distance (lat1.d, lon1.d, lat2.d, lon2.d)
; Geodetische afstand in kilometer tussen twee punten lat1 lon1 en lat2 lon2
; Coordinaten uitgedrukt in decimale gradianen.
; Wiskundig: Vincenty inverse formule voor een ellipsoide (afgeplatte aardbol).
;=================================================================================
; My source:
; Code has been ported by lost_species from www.aliencoffee.co.uk to VBA from javascript published at:
; http://www.movable-type.co.uk/scripts/latlong-vincenty.html
; * from: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the
; * Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975
; * http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
;Constants for WGS-84 ellipsiod:
Define.d radian_equator = 6378137
Define.d radian_poles = 6356752.314245
Define.d flatt = 1 / 298.257223563 ;Flattening of the elipsoid, (afplatting van de aarde)
;#Define.d Radian(x) = ( (x) * 3.1415926535897932384626433 / 180 )
;Info about Freebasic's #Define:
;#define allows to declare text-based preprocessor macros.
;Once the compiler has seen a #define,
;it will start replacing further occurrences of identifier with body, body may be empty.
;The expansion is done recursively,
;until there is nothing more to expand and the compiler can continue analyzing the resulting code.
;#undef can be used to make the compiler forget about a #define.
;#define maakt het mogelijk om tekstgebaseerde preprocessor-macro's te declareren.
;Zodra de compiler een #define heeft gezien,
;zal deze het verder aantreffen van de identifier vervangen door 'body', body mag leeg zijn.
;De uitbreiding gebeurt recursief,
;tot er niets meer is om uit te breiden en de compiler de resulterende code kan blijven analyseren.
;#undef kan worden gebruikt om de compiler een #define te laten vergeten.
;Iterate until lambda negligible (10^-12 = 0.006mm)
;Dim Iterate_until As Double = 0.000000000001
;With the value below, only 30cm deviation is obtained on the total length of the equator
Define.d Iterate_until = 0.00000001
Define.l number_of_iterations
Define.d L, U1, U2 , sinU1 , sinU2 , cosU1 , cosU2 , Lambda , LambdaP
Define.d sinLambda
Define.d cosLambda
Define.d sinSigma
Define.d cosSigma
Define.d sigma
Define.d sinAlpha
;mathematical word Duth 'kwadraat' = number ^ 2"
Define.d coskwadraatAlpha , cos2SigmaM , C , uKwadraat , upper_A ,upper_B
Define.d deltaSigma , dist
;Freebasic:
;#Define.d Radian(x) = ( (x) * 3.1415926535897932384626433 / 180 )
;L = Radian(lon2 - lon1) ;= difference in longitude
;For Purebasic (less digits otherwise PB-error...) and function changed to:
L= (lon2 -lon1) * 3.1415926535897932384 / 180
;L = Radian(lon2 - lon1) ;= difference in longitude
U1 = ATan((1 - flatt) * Tan(Radian(lat1))) ;reduced latitude
U2 = ATan((1 - flatt) * Tan(Radian(lat2)))
sinU1 = Sin(U1)
cosU1 = Cos(U1)
sinU2 = Sin(U2)
cosU2 = Cos(U2)
Lambda = L
LambdaP = 2 * 3.141592653589793238
;number_of_iterations = 10
number_of_iterations = 6
While (Abs(Lambda - LambdaP) > Iterate_until) And (number_of_iterations > 0)
number_of_iterations = number_of_iterations - 1
sinLambda = Sin(Lambda)
cosLambda = Cos(Lambda)
;=============================================================================================
;Freebasic:
;sinSigma = Sqr(((cosU2 * sinLambda) ^ 2) + ((cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ^ 2))
;----------------------------------------------------------------------------------------------
;Purebasic: Pow((cosU2 * sinLambda), 2)+Pow((cosU1 * sinU2 - sinU1 * cosU2 * cosLambda), 2)
;square root: Result.f(.d) = Sqr(Number.f(.d))
;=============================================================================================
sinSigma = Sqr(Pow((cosU2 * sinLambda), 2)+Pow((cosU1 * sinU2 - sinU1 * cosU2 * cosLambda), 2))
If sinSigma = 0 ;Coincident points (overlappende punten)
ProcedureReturn 0
EndIf
cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
;===================================================================================
;Freebasic: result = ATan2( y, x )
;y Vertical And x Horizontal component of the ratio.
;Return Value : The angle whose tangent is y/x, in radians, in the range [-Pi..Pi].
;-----------------------------------------------------------------------------------
;Purebasic; result.f(.d) = ATan2(x.f(.d), y.f(.d))
; This function calculates the value ATan(y/x)
;RUDY: NOTE THAT THE PARAMETERS y,x ARE IN EACH OTHER'S PLACE IN BOTH LANGUAGES!!!
;===================================================================================
;Freebasic: sigma = ATan2(sinSigma, cosSigma)
;For Purebasic changed to:
sigma = ATan2(cosSigma, sinSigma)
sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
coskwadraatAlpha = 1 - Pow(sinAlpha ,2)
;Avoide division bij zero:
If coskwadraatAlpha = 0 ;Two points on the equator.
cos2SigmaM = 0
Else
cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / coskwadraatAlpha
EndIf
C = f / 16 * coskwadraatAlpha * (4 + f * (4 - 3 * coskwadraatAlpha))
LambdaP = Lambda
Wend
If number_of_iterations < 1
OpenConsole ()
;PrintN( "Maximum aantal iteraties is uitgevoerd.
;PrintN ("Oorzaken: de ingegeven Data en of de grenzen van Vincenty formule.")
PrintN( "Maximum number of iterations has been performed.")
PrintN( "Causes: the entered Data And Or the limits of Vincenty formula.")
Input ()
CloseConsole ()
ProcedureReturn 0
EndIf
uKwadraat = coskwadraatAlpha * (Pow(radian_equator,2) - Pow(radian_poles,2)) / (Pow(radian_poles,2))
upper_A = 1 + uKwadraat / 16384 * (4096 + uKwadraat * (-768 + uKwadraat * (320 - 175 * uKwadraat)))
upper_B = uKwadraat / 1024 * (256 + uKwadraat * (-128 + uKwadraat * (74 - 47 * uKwadraat)))
deltaSigma = upper_B * sinSigma * (cos2SigmaM + upper_B / 4 * (cosSigma * (-1 + 2 * Pow(cos2SigmaM, 2)) - upper_B / 6 * cos2SigmaM * (-3 + 4 * Pow(sinSigma , 2)) * (-3 + 4 * Pow(cos2SigmaM,2))))
;Result in meter:
;dist = (radian_poles * upper_A * (sigma - deltaSigma))
;Devision by 1000 for km:
dist = (radian_poles * upper_A * (sigma - deltaSigma)) / 1000
;dist=Round(dist,3) freebasic
ProcedureReturn dist
EndProcedure
;Main
;Define.d Radian(x) = ( (x) * 3.1415926535897932384626433 / 180 )
;Define.s x
OpenConsole()
;Three chapels on a leyline in Belgium (added for my Belgian and Dutch friends leyline investigators)
;For them I did not delete the Dutch texts...
;(purposes:
;Later calculation of the deviation of the K8 from the straight 'projectionline' K7-K9) on a flat surface.
;Idem for French villages on straight lines conc. the "Eleusis Alaisia" book from the Xavier Guichard(+) FR
; http://www.ancient-wisdom.com/eleuse%20alaise%20chapter1.htm
;PrintN( "Three places (chapels) in the village of Retie in Belgium:")
;PrintN( "Distance between LC1 and K7 under Retie is :"+ StrD(distance(51.28382, 5.01033, 51.26397, 5.07559),3))
;PrintN( "Distance between LC1 and K8 under Retie is :" +StrD(distance(51.28382, 5.01033, 51.25774, 5.09963),3))
;PrintN( "Distance between LC1 and K9 under Retie is :" +StrD(distance(51.28382, 5.01033, 51.25193, 5.12125),3))
;PrintN( "De afstand tussen LC1 en K7 onder Retie is :"+ StrD(distance(51.28382, 5.01033, 51.26397, 5.07559),3))
;PrintN( "De afstand tussen LC1 en K8 onder Retie is :" +StrD(distance(51.28382, 5.01033, 51.25774, 5.09963),3))
;PrintN( "De afstand tussen LC1 en K9 onder Retie is :" +StrD(distance(51.28382, 5.01033, 51.25193, 5.12125),3))
;PrintN(" ")
;PrintN( "Distance between chapel K7 and K8 under Retie is :" +StrD(distance(51.26397, 5.07559, 51.25774, 5.09963),3))
;PrintN( "Distance between chapel K8 and K9 under Retie is :" +StrD(distance(51.25774, 5.09963, 51.25193, 5.12125),3))
;PrintN( "Distance between chapel K7 and K9 under Retie is :" +StrD(distance(51.26397, 5.07559, 51.25193, 5.12125),3))
;PrintN( "De afstand tussen Kapel K7 en Kapel K8 onder Retie is :" +StrD(distance(51.26397, 5.07559, 51.25774, 5.09963),3))
;PrintN( "De afstand tussen Kapel K8 en Kapel K9 onder Retie is :" +StrD(distance(51.25774, 5.09963, 51.25193, 5.12125),3))
;PrintN( "De afstand tussen Kapel K7 en Kapel K9 onder Retie is :" +StrD(distance(51.26397, 5.07559, 51.25193, 5.12125),3)) PrintN(" ")
PrintN(" ")
PrintN( "Wellknown places and 'towers' in The Netherlands:")
PrintN( "Distance Palace (Dam) to the equator with UTM coordinates is.: 5804,225 km.")
PrintN( "Distance Palace (Dam) to the equator is, according this listing: " +StrD(distance(52.373126, 4.892467, 0.00, 4.892467),3))
PrintN(" ")
PrintN( "Distance Amersfoort L-Vrouwen-tower to Palace on the Dam is : " +StrD(distance(52.15517440, 5.38720621, 52.373126, 4.892467),3))
PrintN( "Distance Amersfoort L-Vrouwen-tower to A-dam Wester-tower is : " +StrD(distance(52.15517440, 5.38720621, 52.37453253, 4.88352559),3))
PrintN( "Geoflifeline.com: distance L-Vrouwen-tower to Martini-tower : 142.83216 km")
PrintN( " this listing calculate: " +StrD(distance(52.15517440, 5.38720621, 53.21938317, 6.56820053),3))
;PrintN( "De afstand tussen het Paleis op de Dam tot de evenaar is volgens de UTM coordinaten 5804,225 km.")
;PrintN( "De afstand tussen het Paleis op de Dam tot de evenaar is volgens deze berekening..:" +StrD(distance(52.373126, 4.892467, 0.00, 4.892467),3))
;PrintN(" ")
;PrintN( "De afstand Amersfoort Lieve Vrouwentoren 31U tot het Paleis op de Dam 31U is :" +StrD(distance(52.15517440, 5.38720621, 52.373126, 4.892467),3))
;PrintN( "De afstand Amersfoort Lieve Vrouwentoren 31U tot Amsterdam Westertoren 31U is :" +StrD(distance(52.15517440, 5.38720621, 52.37453253, 4.88352559),3))
;PrintN( "Volgens Geoflifeline.com is de afstand tussen Lieve Vrouwetoren en Martinitoren : 142.83216 kilometer")
;PrintN( "De afstand Amersfoort Lieve Vrouwentoren 31U tot Groningen Martinitoren 32U is :" +StrD(distance(52.15517440, 5.38720621, 53.21938317, 6.56820053),3))
PrintN(" ")
PrintN( "From Paris,France to the equator and to the Nordpole")
PrintN( "Calculated on GeoLifeLine.com 5392.92947 km")
PrintN( " following this listing: " +StrD(distance(48.67, 2.33, 0, 2.33),3))
PrintN( "Calculated ony GeoLifeLine 4609.03626 km")
PrintN( " following this listing: " +StrD(distance(48.67, 2.33, 90, 2.33),3))
PrintN(" ")
PrintN( "USA distances between Airports")
PrintN( "Distance between Nashville International Airport (BNA), USA")
PrintN( " and Los Angeles International Airport (LAX), USA")
PrintN( "Rosettacode says between 2886.444 en 2887.259 km")
PrintN( "Calculation GeoLifeLine 2884.70058 km")
PrintN( "this listing calculates..:" +StrD(distance(36.12, -86.76, 33.94, -118.40),3))
PrintN(" ")
PrintN(" ")
PrintN( "Distance between 'LA' International Airport (LAX), USA")
PrintN( " and (JFK Airport), USA: " + StrD(distance(33.95000, -118.40000, 40.63333,-73.78333),3))
PrintN(" ")
PrintN(" ")
PrintN(" ")
PrintN( "A few tests, for fun:")
PrintN( "Test how much km in one latitudee (Belgium = 51)")
PrintN( "The distance for 1 latitude in Belgium is :" +StrD(distance(51, 1, 51, 2),3))
PrintN( "The same ont het equator for 1 latitude :" +StrD(distance(0, 1, 0, 2),3))
PrintN(" ")
PrintN( "4 times 90 degrees = length of the equator :" +StrD(4 * distance(0, 1, 0, 91),3))
PrintN( " WGS 84 voor gps 'says' equator is 40075.0167 km long")
PrintN(" ")
PrintN(" ")
PrintN("Results of my old Freebasic listing:")
PrintN(" The distance for 1 latitude in Belgium is : 70.197 kilometer")
PrintN("The same ont het equator for 1 latitude 111.319 kilometer")
PrintN(" 4 times 90 degrees = length of the equator: 40075.016 kilometer")
Input()
CloseConsole()