10-30-2025, 02:35 PM
Hello,
here is my program for calculating distances on the Earth's surface based on the Vincenty formulas.
I'm going to use it to find the alignment of buildings in Western Europe on the computer. (Three points lie on a straight line if the distance from A to B plus the distance from B to C equals the distance from A to B.)
Less accuracy (plus higher speed) could also be achieved using the Haversine formula.
More accuracy can also be achieved by using the C++ library:
https://geographiclib.sourceforge.io/htm.../code.html
subpage;
https://geographiclib.sourceforge.io/htm...ic.Inverse
In Linux You can install easy GeographicLib via "Synaptic" (Windows and Mac: See there website please).
I need to search out the right name(s) of their functions.
With this lib You can do all possible "geo-calculations".
A very interesting function is concerns the so named "bearing" (p.e. the calculation of the deviation of a ship or aeroplane from a line between to points).
Regards, Rudy
First my Vincenty listing:
And now the little but less accurate Haversine:
here is my program for calculating distances on the Earth's surface based on the Vincenty formulas.
I'm going to use it to find the alignment of buildings in Western Europe on the computer. (Three points lie on a straight line if the distance from A to B plus the distance from B to C equals the distance from A to B.)
Less accuracy (plus higher speed) could also be achieved using the Haversine formula.
More accuracy can also be achieved by using the C++ library:
https://geographiclib.sourceforge.io/htm.../code.html
subpage;
https://geographiclib.sourceforge.io/htm...ic.Inverse
In Linux You can install easy GeographicLib via "Synaptic" (Windows and Mac: See there website please).
I need to search out the right name(s) of their functions.
With this lib You can do all possible "geo-calculations".
A very interesting function is concerns the so named "bearing" (p.e. the calculation of the deviation of a ship or aeroplane from a line between to points).
Regards, Rudy
First my Vincenty listing:
Code: (Select All)
OPTION _EXPLICIT
DECLARE LIBRARY "math"
FUNCTION atan2# (BYVAL y AS DOUBLE, BYVAL x AS DOUBLE) 'arc-tangent of y/x designates quadrant
FUNCTION pow# (BYVAL base AS DOUBLE, BYVAL exponent AS DOUBLE) 'base number raised to exponent
END DECLARE
DIM x$
'Three chapels on a leyline in Belgium (added for my Belgian and Dutch friends leyline investigators)
'(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%20a...apter1.htm
PRINT " For places (1 LeyCentra + 3 chapels) in the village of Retie in Belgium:"
PRINT "Distance between LeyC1 L1 to chapel C7 under Retie is :"; distance#(51.28382, 5.01033, 51.26397, 5.07559)
PRINT "Distance between Chapel C7 to chapel C8 under Retie is :"; distance#(51.26397, 5.07559, 51.25774, 5.09963)
PRINT "Distance between Chapel C8 to chapel C9 under Retie is :"; distance#(51.25774, 5.09963, 51.25193, 5.12125)
PRINT "Total distance LeyC1 L1 to chapel C9 under Retie is :"; distance#(51.28382, 5.01033, 51.25193, 5.12125)
PRINT
PRINT "Distance AB +BC = +- distance AC ---> indicating points on a straight line"
PRINT "The deviation of point B from straight line AB can be calculated from the little sum-difference."
PRINT "Deviation of my GPS itself was 5 meters."
PRINT "Comparision result with http://www.geolifeline.com/distancecalc.php = OK "
PRINT "Comparision can also be made via Google satelite-view"
PRINT "Note that each system has his own deviation, due to different calculation and measuring systems."
PRINT "------------------------------------------------------------------------------------------------"
PRINT
SLEEP
PRINT "Wellknown places and 'towers' in The Netherlands:"
PRINT "Distance Palace (Dam) to the equator with UTM coordinates is.: 5804,225 km."
PRINT "Distance Palace (Dam) to the equator is, according this listing: "; distance#(52.373126, 4.892467, 0.00, 4.892467)
PRINT
PRINT "Distance Amersfoort L-Vrouwen-tower to Palace on the Dam is : "; distance#(52.15517440, 5.38720621, 52.373126, 4.892467)
PRINT "Distance Amersfoort L-Vrouwen-tower to A-dam Wester-tower is : "; distance#(52.15517440, 5.38720621, 52.37453253, 4.88352559)
PRINT "Geoflifeline.com: distance L-Vrouwen-tower to Martini-tower : 142.83216 km"
PRINT " this listing calculate: "; distance#(52.15517440, 5.38720621, 53.21938317, 6.56820053)
PRINT "------------------------------------------------------------------------------------------------"
PRINT
SLEEP
'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 volgens deze berekening..:" distance#(52.373126, 4.892467, 0.00, 4.892467)
'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)
'Print "De afstand Amersfoort Lieve Vrouwentoren 31U tot Amsterdam Westertoren 31U is :" distance#(52.15517440, 5.38720621, 52.37453253, 4.88352559)
'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)
PRINT
PRINT "From Paris,France to the equator and to the Nordpole"
PRINT "Calculated on GeoLifeLine.com 5392.92947 km"
PRINT " following this listing: "; distance#(48.67, 2.33, 0, 2.33)
PRINT "Calculated ony GeoLifeLine 4609.03626 km"
PRINT " following this listing: "; distance#(48.67, 2.33, 90, 2.33)
PRINT "------------------------------------------------------------------------------------------------"
PRINT
SLEEP
PRINT "USA distances between Airports"
PRINT "Distance between Nashville International Airport (BNA), USA"
PRINT " and Los Angeles International Airport (LAX), USA"
PRINT "Rosettacode says between 2886.444 en 2887.259 km"
PRINT "Calculation GeoLifeLine 2884.70058 km"
PRINT "this listing calculates..:"; distance#(36.12, -86.76, 33.94, -118.40)
PRINT "------------------------------------------------------------------------------------------------"
PRINT
SLEEP
PRINT
PRINT "Distance between 'LA' International Airport (LAX), USA"
PRINT " and (JFK Airport), USA: "; distance#(33.95000, -118.40000, 40.63333, -73.78333)
PRINT "------------------------------------------------------------------------------------------------"
PRINT
PRINT
PRINT "A few tests, for fun:"
PRINT "Test how much km in one latitudee (Belgium = 51)"
PRINT "The distance for 1 latitude in Belgium is :"; distance#(51, 1, 51, 2)
PRINT "The distance for 1 latitude on the equator:"; distance#(0, 1, 0, 2)
PRINT "------------------------------------------------------------------------------------------------"
PRINT
PRINT "4 times 90 degrees = length of the equator :"; (4 * distance#(0, 1, 0, 91))
PRINT " WGS 84 voor gps 'says' equator is 40075.0167 km long"
PRINT "------------------------------------------------------------------------------------------------"
PRINT
PRINT "Results of my old Freebasic listing:"
PRINT " The distance for 1 latitude in Belgium is: 70.197 km"
PRINT " The distance for 1 latitude on the equator: 111.319 km"
PRINT " 4 times 90 degrees = length of the equator: 40075.016 km"
SLEEP
END
FUNCTION distance# (Lat1#, Lon1#, Lat2#, Lon2#)
'==============================================
' Geodesic distance in kilometers between two points Lat1# Lon1# and Lat2# Lon2#
' Coordinates expressed in decimal gradians, pe: 52.12345
' Mathematical: Based on Vincenty inverse formula for an ellipsoid (flattened globe).
'====================================================================================
' Sources:
' 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/la...centy.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
'------------------------------------------------------------------------------------
'https://en.wikipedia.org/w/index.php?title=Geodetic_Reference_System_1980
'Flattening = f 0.003 352 810 681 183 637 418
'Reciprocal of flattening = 1 / f = 298.257 222 100 882 711 243
'Semi-minor axis = Polar Radius = b = 6356752.314140347 m
'Aspect ratio = b / a = 0.996 647 189 318 816 363;
'------------------------------------------------------------------------------------
'https://en.wikipedia.org/wiki/World_Geodetic_System
'Ellipsoid reference Semi-major axis a Semi-minor axis b Inverse flattening 1/f
' WGS 84 6378137.0 m 6356752.314245 m 298.257223563
'Constants for WGS-84 ellipsiod:
CONST radian_equator = 6378137 '------semi-major axis a
'CONST radian_poles = 6356752.314245
CONST radian_poles = 6356752.314245 '-semi-minor axis b
'flattening of the elipsoid, (flattening of the earth...)
'Flatt= 0.003352810681183637418
CONST Flatt = 1 / 298.257223563 '= 1 / Inverse flattening
DIM x$
DIM decimals%: decimals% = 3 'round km to 1m
DIM dist#
DIM ToRound#
DIM factor%
'Iterate until Lambda# negligible (10^-12 = 0.006mm)
'Iterate_until# = 0.000000000001
'With the value below, only 30cm deviation is obtained on the total length of the equator
DIM Iterate_until#
Iterate_until# = 0.00000001
DIM number_of_itereations%
DIM L#, U1#, U2#, sinU1#, sinU2#, cosU1#, cosU2#, Lambda#, LambdaP#
DIM sinLambda#, cosLambda#, sinSigma#, cosSigma#, sigma#, sinAlpha#
'mathematical word Duth 'kwadraat' = number ^ 2"
DIM coskwadraatAlpha#, cos2SigmaM#, deltasigma#, C#, UKwadraat#, Upper_A#, Upper_B#
'Radian# = degr# * 3.1415926535897932384626433 / 180 (L# = difference in longitude in radian)
L# = (Lon2# - Lon1#) * _PI / 180
U1# = ATN((1 - Flatt) * TAN(Lat1# * _PI / 180)) 'reduced latitude (*_Pi/180 to convert in radian)
U2# = ATN((1 - Flatt) * TAN(Lat2# * _PI / 180))
sinU1# = SIN(U1#)
cosU1# = COS(U1#)
sinU2# = SIN(U2#)
cosU2# = COS(U2#)
Lambda# = L#
LambdaP# = 2 * _PI
number_of_itereations% = 10
'number_of_itereations% = 6
WHILE (ABS(Lambda# - LambdaP#) > Iterate_until#) AND (number_of_itereations% > 0)
number_of_itereations% = number_of_itereations% - 1
sinLambda# = SIN(Lambda#)
cosLambda# = COS(Lambda#)
sinSigma# = SQR(pow#((cosU2# * sinLambda#), 2) + pow#((cosU1# * sinU2# - sinU1# * cosU2# * cosLambda#), 2))
IF sinSigma# = 0 THEN
PRINT "Coincident points (overlappende punten)"
distance# = 0
x$ = INKEY$
EXIT FUNCTION
END IF
cosSigma# = sinU1# * sinU2# + cosU1# * cosU2# * cosLambda#
sigma# = atan2#(sinSigma#, cosSigma#)
sinAlpha# = cosU1# * cosU2# * sinLambda# / sinSigma#
coskwadraatAlpha# = 1 - pow#(sinAlpha#, 2)
'Avoid division by zero:
IF coskwadraatAlpha# = 0 THEN 'Two points on the equator.
cos2SigmaM# = 0
ELSE
cos2SigmaM# = cosSigma# - 2 * sinU1# * sinU2# / coskwadraatAlpha#
END IF
C# = Flatt / 16 * coskwadraatAlpha# * (4 + Flatt * (4 - 3 * coskwadraatAlpha#))
LambdaP# = Lambda#
WEND
IF number_of_itereations% < 1 THEN
'Later MsgBox
PRINT "Maximum number of iterations has been performed."
PRINT "Causes: the entered Data And Or the limits of Vincenty formula."
x$ = INKEY$
distance# = number_of_itereations%
EXIT FUNCTION
END IF
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
'Round km to 0.001 km = 1 meter
factor% = 10 ^ decimals%
distance# = (INT(.5 + factor% * dist#)) / factor%
END FUNCTION 'distance#-------------------------------------------------------------------------------------
And now the little but less accurate Haversine:
Code: (Select All)
SCREEN _NEWIMAGE(800, 100, 32)
'*** Units: K=kilometers M=miles N=nautical miles
DIM UNIT AS STRING
DIM Distance AS STRING
DIM Result AS DOUBLE
DIM ANSWER AS DOUBLE
'*** Change the To/From Latittude/Logitudes for your run
'*** LAT/LON for Nashville International Airport (BNA)
lat1 = 36.12
Lon1 = -86.67
'*** LAT/LONG for Los Angeles International Airport (LAX)
Lat2 = 33.94
Lon2 = -118.40
'*** Initialize Values
UNIT = "K"
Distance = ""
'Radius = 6378.137
Radius = 6372.8
'*** Calculate distance using Haversine Function
lat1 = (lat1 * _PI / 180)
Lon1 = (Lon1 * _PI / 180)
Lat2 = (Lat2 * _PI / 180)
Lon2 = (Lon2 * _PI / 180)
DLon = Lon1 - Lon2
ANSWER = _ACOS(SIN(lat1) * SIN(Lat2) + COS(lat1) * COS(Lat2) * COS(DLon)) * Radius
'*** Adjust Answer based on Distance Unit (kilometers, miles, nautical miles)
SELECT CASE UNIT
CASE "M"
Result = ANSWER * 0.621371192
Distance = "miles"
CASE "N"
Result = ANSWER * 0.539956803
Distance = "nautical miles"
CASE ELSE
Result = ANSWER
Distance = "kilometers"
END SELECT
'*** Change PRINT statement with your labels for FROM/TO locations
PRINT "The distance from Nashville International to Los Angeles International in "; Distance;
PRINT USING " is: ##,###.##"; Result;
PRINT "."
END