Calculate distance on Earth surface via Vincenty formula

Rudy M · Post by **Rudy M** » 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:

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

Output from the freebasic 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

My Purebasic listing:

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()

infratec · Post by **infratec** » Mon Dec 23, 2024 7:27 pm

You should not print something in your distance procedure.
Else you can not use it in a window program.
Simply return -1, -2 ... in such error cases.

And inside of Procedures you should use Protected instead of Define.

And ... always use EnableExplicit as your first command in any program.
Then you can avoid trouble with typos or double defined names.

Here are some other routines to compare:

viewtopic.php?t=47853

Rudy M · Post by **Rudy M** » Tue Dec 24, 2024 9:25 am

@Infratec:
Thanks for the "good programming" tips. I'll try to adjust the listing.
EnableExplicit... indeed, I normally always use it but now I forgot; Shame on my...
I also change the listing to avoid "printing" within the procedure.
Protected... not yet used in PB, just looked in the PB manual and you're right!!!

Thanks for the forum link on the same subject, very instructive!
I know Haversine (that's what I started with back then), but switched to Vincenty because of the higher accuracy.
(I think the online program Geofline.com also works with Vincenty's given the very similar numerical values.)
Greetings,
Rudy M

SMaag · Post by **SMaag** » Tue Dec 24, 2024 3:05 pm

Some comments on your code:

The Purebasic transformation looks well done.

But there is a line with undefined variables
C = f / 16 * coskwadraatAlpha * (4 + f * (4 - 3 * coskwadraatAlpha))

C and f are not defined

You should always use. EnableExplicite
at the beginning of a code. So the compiler checks undefined variables

Some more PureBasic specific things:

Use of define in Procedures:

Define in Procedure is same or nearly same as Protected
For temporary variables in Procedures better use Protected instead of Define.
I guess protected is the newer command and indicates it is a temporary variable created on the Stack.
But in actual Prebasic there is no difference between Define and Protected in a Procedure.
Outside of a Procedure you can use only Define and Global. See Purebasic Help of Define, Global and Protected

For PI you can use directly the integrated PI constand #PI
Some of us define 2*Pi as
#PI2 = 2*#PI

For number_of_iterations better use .i, not .l

Purebasic do not have an integrated command for x². You can use the Power Operator or some of us do, a Square Macro

Code: Select all

Macro Square(X)
 ((X) * (X))
EndMacro

Define var = 4
Define result

result = Square(var)

debug result

Rudy M · Post by **Rudy M** » Tue Dec 24, 2024 5:44 pm

Thank You all very much for the programming hints!!!

I changed the listing with all Your hints,
(Enable Explicit, avoid printing in a proc. and the hint from user Smaag conc for using Purebasics #PI variable).

Not yet implemented: The very usable hint from Smaag to use a macro to calculate ^2
To be honest: Till to day I did not tested out the PB macro possibility's.

While controlling the calculated results I remember me the possibility to calculate the deviation from an inner point B between points A and B.
But that's for later. (Two possibility's: Simpel calculation... or 100% trigonometric calculation...)

So, thanks to You I can present the updated listing.

Rudy M

Code: Select all

EnableExplicit

Procedure.d distance (Lat1.d, Lon1.d, Lat2.d, Lon2.d)

; Geodesic distance in kilometers between two points lat1 lon1 and Lat2 Lon2
; Coordinates expressed in decimal gradians.
; Mathematical: 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/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:
    Protected.d radian_equator = 6378137
    Protected.d radian_poles = 6356752.314245
    Protected.d flatt = 1 / 298.257223563    ;Flattening of the elipsoid, (afplatting van de aarde)
       
    ;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
    Protected.d Iterate_until = 0.00000001   
    
    Protected.l number_of_iterations   
    Protected.d f, L, U1, U2 , sinU1 , sinU2 , cosU1 , cosU2 , Lambda , LambdaP
    Protected.d sinLambda 
    Protected.d cosLambda 
    Protected.d sinSigma 
    Protected.d cosSigma 
    Protected.d sigma
    Protected.d sinAlpha 
    ;mathematical word Duth 'kwadraat' = number ^ 2"
    Protected.d coskwadraatAlpha , cos2SigmaM , C , uKwadraat , upper_A ,upper_B
    Protected.d deltaSigma , calculated
    
    L= (Lon2 -lon1) * #PI / 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)       
        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 
        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 number of iterations has been performed.") 
      PrintN( "Causes: the entered Data And Or the limits of Vincenty formula.")
      Input ()
      CloseConsole ()  
      ProcedureReturn number_of_iterations
    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:
    calculated = (radian_poles * upper_A * (sigma - deltaSigma)) / 1000
       
    ProcedureReturn calculated   

EndProcedure 

OpenConsole()
  PrintN (" ")
  PrintN (" ")
 ;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%20alaise%20chapter1.htm
  PrintN( "         For places (1 LeyCentra + 3 chapels) in the village of Retie in Belgium:")
  PrintN( "Distance between LeyC1  L1 to chapel C7 under Retie is :"+ StrD(distance(51.28382, 5.01033, 51.26397, 5.07559),3))    
  PrintN( "Distance between Chapel C7 to chapel C8 under Retie is :" +StrD(distance(51.26397, 5.07559, 51.25774, 5.09963),3))
  PrintN( "Distance between Chapel C8 to chapel C9 under Retie is :" +StrD(distance(51.25774, 5.09963, 51.25193, 5.12125),3))
  PrintN( "Total distance   LeyC1  L1 to chapel C9 under Retie is :" +StrD(distance(51.28382, 5.01033, 51.25193, 5.12125),3))
  PrintN (" ")
  PrintN( "Distance A tot B = +- distance AB etc...---> indicating a straight line")
  PrintN( "The deviation of point B from straight line AB can be calculated from the little sum-difference.")
  PrintN( "Deviation of my GPS itself was 5 meters.")
  PrintN( "Comparision result with   http://www.geolifeline.com/distancecalc.php  = OK ")
  PrintN( "Comparision can also be made via Google satelite-view")
  PrintN( "Note that each system has his own deviation, due to different calculation and measuring systems.") 
  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 on the 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 km")
  PrintN("The same ont het equator for 1 latitude        111.319 km")
  PrintN(" 4 times 90 degrees = length of the equator: 40075.016 km")
 
  Input()
  
CloseConsole()

Rudy M · Post by **Rudy M** » Tue Dec 24, 2024 5:53 pm

It was possible to edit my previous post and replace the wrong listing.

Now its working again.
26/12/2025
Rudy M

PureBasic Forums - English

Calculate distance on Earth surface via Vincenty formula

Calculate distance on Earth surface via Vincenty formula

Re: Calculate distance on Earth surface via Vincenty formula

Re: Calculate distance on Earth surface via Vincenty formula

Re: Calculate distance on Earth surface via Vincenty formula

Re: Calculate distance on Earth surface via Vincenty formula

Re: Calculate distance on Earth surface via Vincenty formula