PureBasic Forums - English

Posted: **Fri Oct 14, 2011 4:00 am**

To calculate distances between two points on the Earth using latitude and longitude?

http://www.movable-type.co.uk/scripts/latlong.html

I'm sure it has to be fairly simple using PB.

Posted: **Fri Oct 14, 2011 12:49 pm**

Hi,

no, it's not difficult. But...

I think PB has a bug in the ATan2(x, y) function.
In my opinion are the operands x and y switched!

Here the code for haversine:

Code: Select all

Procedure.f haversine(lat1.f, lon1.f, lat2.f, lon2.f)
  
  dLat.f = Radian(lat2 - lat1)
  dLon.f = Radian(lon2 - lon1)
  lat1 = Radian(lat1)
  lat2 = Radian(lat2)
  
  a.f = (Sin(dLat / 2) * Sin(dLat / 2)) + ((Cos(lat1) * Cos(lat2)) * (Sin(dLon / 2) * Sin(dLon / 2)))
    
  ProcedureReturn 2 * (ATan2(Sqr(a), Sqr(1 - a)))
  
EndProcedure

Debug haversine(50.0, 5.0, 58.0, 3.0) * 6372.797560856

But it gives wrong Results.

This one is working:

Code: Select all

Procedure.f haversine(lat1.f, lon1.f, lat2.f, lon2.f)
  
  dLat.f = Radian(lat2 - lat1)
  dLon.f = Radian(lon2 - lon1)
  lat1 = Radian(lat1)
  lat2 = Radian(lat2)
  
  a.f = (Sin(dLat / 2) * Sin(dLat / 2)) + ((Cos(lat1) * Cos(lat2)) * (Sin(dLon / 2) * Sin(dLon / 2)))
  
  ProcedureReturn 2 * (ATan2(Sqr(1 - a), Sqr(a)))
  
EndProcedure

Debug haversine(50.0, 5.0, 58.0, 3.0) * 6372.797560856

Bernd

P.S.: I opened a bug thread for that.
http://www.purebasic.fr/english/viewtop ... =4&t=47855

Posted: **Fri Oct 14, 2011 5:20 pm**

An optimized version

Code: Select all

; <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
; Project name : Haversine Calculations
; File Name : Haversine.pb
; File version: 1.0.1
; Programmation : OK
; Programmed by : Guimauve
; Date : 14-10-2011
; Last Update : 14-10-2011
; PureBasic code : 4.60
; Plateform : Windows, Linux, MacOS X
; <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

#Earth_Radius = 6372.8 ; km
#Moon_Radius = 1737.4 ; km

Procedure.d DMS_To_DD(Direction.s, Degrees.d, Minutes.d, Seconds.d)
  
  If Direction = "N" Or Direction = "E"
    AngleDD.d = Degrees + Minutes / 60 + Seconds / 3600
  ElseIf Direction = "S" Or Direction = "W"
    AngleDD = -1 * (Degrees + Minutes / 60 + Seconds / 3600)
  EndIf
  
  ProcedureReturn AngleDD
EndProcedure 

Procedure.d Haversine(PlanetRadius.d, Lat1.d, Lon1.d, Lat2.d, Lon2.d)

  Sin_dLat.d = Sin(Radian(Lat2 - Lat1) / 2)
  Sin_dLon.d = Sin(Radian(Lon2 - Lon1) / 2)
  a.d = Sin_dLat * Sin_dLat + Sin_dLon * Sin_dLon * Cos(Radian(Lat1)) * Cos(Radian(Lat2))
  
  ProcedureReturn PlanetRadius * 2 * ATan2(Sqr(1 - a), Sqr(a))
EndProcedure 

Lat1.d = DMS_To_DD("N", 50, 3, 59)
Lon1.d = DMS_To_DD("W", 5, 42, 53)

Lat2.d = DMS_To_DD("N", 58, 38, 38)
Lon2.d = DMS_To_DD("W", 3, 4, 12)

Debug Haversine(#Earth_Radius, Lat1, Lon1, Lat2, Lon2)
Debug Haversine(#Moon_Radius, Lat1, Lon1, Lat2, Lon2)

; <<<<<<<<<<<<<<<<<<<<<<<
; <<<<< END OF FILE <<<<<
; <<<<<<<<<<<<<<<<<<<<<<<

This version can be used to calculate distance on any planet or moon.

Edit :

DMS_To_DD() conversion function added.

Best regards
Guimauve

Posted: **Fri Oct 14, 2011 6:30 pm**

Thanks so much folks! I was sure someone had worked with it here.

Posted: **Sat Oct 15, 2011 2:40 pm**

The other function on the page (Spherical Law) also seems to work pretty well

Code: Select all

Procedure.d SphericalLaw(PlanetRadius.d, Lat1.d, Lon1.d, Lat2.d, Lon2.d)
  ProcedureReturn ACos(Sin(Radian(Lat1)) * Sin(Radian(Lat2)) + Cos(Radian(Lat1)) * Cos(Radian(Lat2)) * Cos(Radian(Lon2 - Lon1))) * PlanetRadius
EndProcedure

Posted: **Sun Oct 16, 2011 12:42 am**

You could do that SphericalLaw procedure as a macro.

Code: Select all

Macro SphericalLaw(PlanetRadius, Lat1, Lon1, Lat2, Lon2)
  ACos(Sin(Radian(Lat1)) * Sin(Radian(Lat2)) + Cos(Radian(Lat1)) * Cos(Radian(Lat2)) * Cos(Radian(Lon2 - Lon1))) * PlanetRadius
EndMacro

Posted: **Sun Oct 16, 2011 1:00 am**

Say Sey, what would be the benefit of using a macro?

I'm tracking mileage in Michigan, latitude 42.7 North. The spherical law of cosines seems to be less accurate than Haversine the further you get from the equator.

Posted: **Sun Oct 16, 2011 6:40 am**

RichAlgeni wrote:The spherical law of cosines seems to be less accurate than Haversine the further you get from the equator.

I didn't check for other places.

I also tried to convert the Haversine formula to asm (only checked the coordinates from the example by Guimauve).
It's a bit faster and has the potential to be a bit more accurate since the fpu works with 80 bits and doesn't have to convert the temporary values to 64 bits to store them in a variable but the difference is not that big.

Code: Select all

Procedure.d HaverSine(PlanetRadius.d, Lat1.d, Lon1.d, Lat2.d, Lon2.d)
  EnableASM
  FLD Lat2
  FLD Lat1
  FLD Lon2
  FLD Lon1
  !fsubp st1, st0
  !fld qword [HaversineDegRad]
  !fmul st1, st0
  !fmul st2, st0
  !fmulp st3, st0
  !fld st2
  !fsub st0, st2
  !fld dword [HaversineHalf]
  !fmul st1, st0
  !fmulp st2, st0
  !fsin
  !fmul st0, st0
  !fstp st4
  !fsin
  !fmul st0, st0
  !fstp st4
  !fcos
  !fmulp st3, st0
  !fcos
  !fmulp st2, st0
  !faddp st1, st0
  ; st0 = a
  !fld1
  !fsub st0, st1
  !fstp st2
  !fsqrt
  !fstp st2
  !fsqrt
  !fpatan
  !fadd st0, st0
  FLD PlanetRadius
  !fmulp st1, st0
  DisableASM
  ProcedureReturn
  !HaversineHalf:
  !dd 0x3f000000
  !HaversineDegRad:
  !dd 0xa2529d39, 0x3f91df46
EndProcedure

To be complete also SphericalLaw converted to asm

Code: Select all

Procedure.d SphericalLaw(PlanetRadius.d, Lat1.d, Lon1.d, Lat2.d, Lon2.d)
  EnableASM
  !fld qword [SphericalLawDegRad]
  FLD Lon2
  FLD Lon1
  !fsubp st1, st0
  !fmul st0, st1
  !fcos
  FLD Lat2
  !fmul st0, st2
  !fsincos
  FLD Lat1
  !fmul st0, st4
  !ffree st4
  !fsincos
  !fmulp st2, st0
  !fmulp st2, st0
  !fmulp st2, st0
  !faddp st1, st0
  !fld1
  !fld st1
  !fmul st0, st0
  !fsubp st1, st0
  !fsqrt
  !fstp st2
  !fpatan
  FLD PlanetRadius
  !fmulp st1, st0
  DisableASM
  ProcedureReturn
  !SphericalLawDegRad:
  !dd 0xa2529d39, 0x3f91df46
EndProcedure

Personally I don't see a big difference in results.

Posted: **Sun Oct 16, 2011 3:33 pm**

Just the right place to jump in.

Did anybody try to implement Vincenty formula?

I tried to convert this JavaScript but failed:

Code: Select all

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
/* Vincenty Inverse Solution of Geodesics on the Ellipsoid (c) Chris Veness 2002-2011             */
/*                                                                                                */
/* 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                                             */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

/**
 * Calculates geodetic distance between two points specified by latitude/longitude using 
 * Vincenty inverse formula for ellipsoids
 *
 * @param   {Number} lat1, lon1: first point in decimal degrees
 * @param   {Number} lat2, lon2: second point in decimal degrees
 * @returns (Number} distance in metres between points
 */
function distVincenty(lat1, lon1, lat2, lon2) {
  var a = 6378137, b = 6356752.314245,  f = 1/298.257223563;  // WGS-84 ellipsoid params
  var L = (lon2-lon1).toRad();
  var U1 = Math.atan((1-f) * Math.tan(lat1.toRad()));
  var U2 = Math.atan((1-f) * Math.tan(lat2.toRad()));
  var sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
  var sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
  
  var lambda = L, lambdaP, iterLimit = 100;
  do {
    var sinLambda = Math.sin(lambda), cosLambda = Math.cos(lambda);
    var sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) + 
      (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
    if (sinSigma==0) return 0;  // co-incident points
    var cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
    var sigma = Math.atan2(sinSigma, cosSigma);
    var sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
    var cosSqAlpha = 1 - sinAlpha*sinAlpha;
    var cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
    if (isNaN(cos2SigmaM)) cos2SigmaM = 0;  // equatorial line: cosSqAlpha=0 (§6)
    var 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)));
  } while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0);

  if (iterLimit==0) return NaN  // formula failed to converge

  var uSq = cosSqAlpha * (a*a - b*b) / (b*b);
  var A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
  var B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
  var deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
    B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
  var s = b*A*(sigma-deltaSigma);
  
  s = s.toFixed(3); // round to 1mm precision
  return s;
  
  // note: to return initial/final bearings in addition to distance, use something like:
  var fwdAz = Math.atan2(cosU2*sinLambda,  cosU1*sinU2-sinU1*cosU2*cosLambda);
  var revAz = Math.atan2(cosU1*sinLambda, -sinU1*cosU2+cosU1*sinU2*cosLambda);
  return { distance: s, initialBearing: fwdAz.toDeg(), finalBearing: revAz.toDeg() };
}

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

My attempt:

Code: Select all

Procedure.d distVincenty(lat1.d, lon1.d, lat2.d, lon2.d)
  Protected a.d=6378.137, b.d= 6356.752, f.d=(a-b)/a;1/298.257223563
  Protected rlat1.d=Radian(lat1)
  Protected rlat2.d=Radian(lat2)
  Protected rlon1.d=Radian(lon1)
  Protected rlon2.d=Radian(lon2)
  Protected L.d=Radian(lon2-lon1)
  Protected U1.d= ATan(1-f)*Tan(rlat1)
  Protected U2.d= ATan(1-f)*Tan(rlat2)
  Protected sinU1.d=Sin(U1)
  Protected cosU1.d=Cos(U1)
  Protected sinU2.d=Sin(U2)
  Protected cosU2.d=Cos(U2)
  
  Protected lambda.d=L, lambdaP.d, iterLimit.i=100
  
  Repeat
    Protected sinLambda.d=Sin(lambda)
    Protected cosLambda.d=Cos(lambda)
    Protected sinSigma.d=Sqr((cosU2*sinLambda) * (cosU2*sinLambda) + (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda))
    If sinSigma=0
      ProcedureReturn 0
    EndIf
    Protected cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda
    Protected sigma.d=ATan2(cosSigma, sinSigma)
    Protected sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
    Protected cosSqAlpha = 1 - sinAlpha*sinAlpha
    Protected cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha
    If (IsNAN(cos2SigmaM))
      cos2SigmaM=0
    EndIf
    Protected C.d=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)))
    iterLimit-1
  Until Abs(lambda-lambdaP)<1e-12 Or iterLimit<1
  If iterLimit=0
    ProcedureReturn NaN()
  EndIf
  Protected uSq.d=cosSqAlpha * (a*a - b*b) / (b*b)
  Protected AA.d=1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))
  Protected BB.d=uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)))
  Protected deltaSigma.d=BB*sinSigma*(cos2SigmaM+BB/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-BB/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)))

 Protected s.d = b*AA*(sigma-deltaSigma)
  
  ProcedureReturn s
EndProcedure

I noticed, that I get an N-S error, the E-W part seems to be ok.
Where is my error?

Posted: **Sun Oct 16, 2011 4:38 pm**

And the error is .... here:

Code: Select all

Protected U1.d= ATan(1-f)*Tan(rlat1)
Protected U2.d= ATan(1-f)*Tan(rlat2)

Corected:

Code: Select all

Protected U1.d= ATan((1-f)*Tan(rlat1))
Protected U2.d= ATan((1-f)*Tan(rlat2))

Bernd

Posted: **Sun Oct 16, 2011 5:38 pm**

To make life easier:
http://www.purebasic.fr/english/viewtop ... 12&t=47881

Bernd

Posted: **Sun Oct 16, 2011 6:26 pm**

Since I have now so many input,
I combined all together:

Code: Select all

#Version = "1.01"

#MiddleEarthRadius = 6371.000785

;GRS 80 WGS84
#ERadius = 6378.1370
#PRadius = 6356.752314

Enumeration
  #CalcButton = 100
  #Lat1Gadget
  #Long1Gadget
  #Lat2Gadget
  #Long2Gadget
  #HaversineGadget
  #SphericalGadget
  #VincentyGadget
EndEnumeration



Procedure.f Haversine(lat1.f, lon1.f, lat2.f, lon2.f)

  dLat.f = Radian(lat2 - lat1)
  dLon.f = Radian(lon2 - lon1)
  lat1 = Radian(lat1)
  lat2 = Radian(lat2)

  a.f = (Sin(dLat / 2) * Sin(dLat / 2)) + ((Cos(lat1) * Cos(lat2)) * (Sin(dLon / 2) * Sin(dLon / 2)))

  ProcedureReturn 2 * (ATan2(Sqr(1 - a), Sqr(a))) * #MiddleEarthRadius

EndProcedure



Procedure.d SphericalLaw(Lat1.d, Lon1.d, Lat2.d, Lon2.d)
  ProcedureReturn ACos(Sin(Radian(Lat1)) * Sin(Radian(Lat2)) + Cos(Radian(Lat1)) * Cos(Radian(Lat2)) * Cos(Radian(Lon2 - Lon1))) * #MiddleEarthRadius
EndProcedure



Procedure.d Vincenty(lat1.d, lon1.d, lat2.d, lon2.d)
  Protected a.d=#ERadius, b.d=#PRadius, f.d=(a-b)/a
  Protected rlat1.d=Radian(lat1)
  Protected rlat2.d=Radian(lat2)
  Protected rlon1.d=Radian(lon1)
  Protected rlon2.d=Radian(lon2)
  Protected L.d=Radian(lon2-lon1)
  Protected U1.d= ATan((1-f)*Tan(rlat1))
  Protected U2.d= ATan((1-f)*Tan(rlat2))
  Protected sinU1.d=Sin(U1)
  Protected cosU1.d=Cos(U1)
  Protected sinU2.d=Sin(U2)
  Protected cosU2.d=Cos(U2)

  Protected lambda.d=L, lambdaP.d, iterLimit.i=100

  Repeat
    Protected sinLambda.d=Sin(lambda)
    Protected cosLambda.d=Cos(lambda)
    Protected sinSigma.d=Sqr((cosU2*sinLambda) * (cosU2*sinLambda) + (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda))
    If sinSigma=0
      ProcedureReturn 0
    EndIf
    Protected cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda
    Protected sigma.d=ATan2(cosSigma, sinSigma)
    Protected sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
    Protected cosSqAlpha = 1 - sinAlpha*sinAlpha
    Protected cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha
    If (IsNAN(cos2SigmaM))
      cos2SigmaM=0
    EndIf
    Protected C.d=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)))
    iterLimit-1
  Until Abs(lambda-lambdaP)<1e-12 Or iterLimit<1
  If iterLimit=0
    ProcedureReturn NaN()
  EndIf
  Protected uSq.d=cosSqAlpha * (a*a - b*b) / (b*b)
  Protected AA.d=1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))
  Protected BB.d=uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)))
  Protected deltaSigma.d=BB*sinSigma*(cos2SigmaM+BB/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-BB/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)))

Protected s.d = b*AA*(sigma-deltaSigma)

  ProcedureReturn s
EndProcedure



Procedure.f ConvertGeoToDeg(Value$)
 
  Protected Result.f = 0.0
  Protected Direction$
  Protected Direction = 1
  Protected Pos
 
  ReplaceString(Value$, ",", ".", #PB_String_InPlace)
 
  If FindString(Value$, "°", 1) = 0
    ProcedureReturn ValF(Value$)
  EndIf
 
  ReplaceString(Value$, "''", #DQUOTE$ + " ", #PB_String_InPlace)
 
  Direction$ = Left(Value$, 1)
  If Direction$ > "A"
    Value$ = Mid(Value$, 2)
  Else
    Direction$ = Right(Value$, 1)
  EndIf
 
  If Direction$ = "S" Or Direction$ = "W"
    Direction = -1
  EndIf
 
  Pos = FindString(Value$, "°", 1)
  If Pos
    Result = Val(Value$)
    Value$ = Mid(Value$, Pos + 1)
    Pos = FindString(Value$, "'", 1)
    If Pos
      Result + (Val(Value$) / 60)
      Value$ = Mid(Value$, Pos + 1)
      Pos = FindString(Value$, #DQUOTE$, 1)
      If Pos
        Result + (ValF(Value$) / 3600)
      EndIf
    EndIf
  EndIf
 
  ProcedureReturn Result * Direction
 
EndProcedure


Lat1.f = 0
Long1.f = 0
Lat2.f = 0
Long2.f = 0


OpenWindow(0, 0, 0, 290, 260, "Geographic distance calculator V" + #VERSION, #PB_Window_SystemMenu|#PB_Window_ScreenCentered)


TextGadget(1, 60, 10, 50, 20, "Latitude")
TextGadget(2, 170, 10, 50, 20, "Longitude")

TextGadget(3, 10, 40, 50, 20, "Point 1")
StringGadget(#Lat1Gadget, 60, 40, 100, 20, "")
StringGadget(#Long1Gadget, 170, 40, 100, 20, "")

TextGadget(4, 10, 70, 50, 20, "Point 2")
StringGadget(#Lat2Gadget, 60, 70, 100, 20, "")
StringGadget(#Long2Gadget, 170, 70, 100, 20, "")

ButtonGadget(#CalcButton, 10, 110, 260, 30, "Calculate")

TextGadget(5, 10, 160, 100, 20, "Haversine:")
StringGadget(#HaversineGadget, 120, 160, 80, 20, "")
TextGadget(6, 210, 160, 20, 20, "km")

TextGadget(7, 10, 190, 100, 20, "Spherical law:")
StringGadget(#SphericalGadget, 120, 190, 80, 20, "")
TextGadget(8, 210, 190, 20, 20, "km")

TextGadget(9, 10, 220, 100, 20, "Vincenty:")
StringGadget(#VincentyGadget, 120, 220, 80, 20, "")
TextGadget(10, 210, 220, 20, 20, "km")

SetActiveGadget(#Lat1Gadget)

Exit = #False
Repeat
 
  Event = WaitWindowEvent()
 
  Select Event
    Case #PB_Event_Gadget
      If EventGadget() = #CalcButton
        Lat1 = ConvertGeoToDeg(GetGadgetText(#Lat1Gadget))
        Long1 = ConvertGeoToDeg(GetGadgetText(#Long1Gadget))
        Lat2 = ConvertGeoToDeg(GetGadgetText(#Lat2Gadget))
        Long2 = ConvertGeoToDeg(GetGadgetText(#Long2Gadget))
       
        SetGadgetText(#HaversineGadget, StrF(Haversine(Lat1, Long1, Lat2, Long2)))
        SetGadgetText(#SphericalGadget, StrF(SphericalLaw(Lat1, Long1, Lat2, Long2)))
        SetGadgetText(#VincentyGadget, StrF(Vincenty(Lat1, Long1, Lat2, Long2)))
      EndIf
    Case #PB_Event_CloseWindow
      Exit = #True
  EndSelect
 
Until Exit

Have fun,

Bernd

P.S.: I don't know why Vincenty gives a bit other results.

Posted: **Sun Oct 16, 2011 6:33 pm**

If you look at the earth from the side, it's not a perfect circle but more like an ellipse.
It's a bit flattened at the poles. Haversine and SphericalLaw ignore this so they are a bit less accurate.

Posted: **Sun Oct 16, 2011 6:41 pm**

Hi,

V1.01 is born.

Small correction of the middle earth radius,
and constants for Vincenty added.
Activate the first input gadget.
Now you can use , or . as decimal separator.

Posted: **Sun Oct 16, 2011 6:58 pm**

infratec wrote:And the error is .... here:
Code: Select all
Protected U1.d= ATan(1-f)*Tan(rlat1)
Protected U2.d= ATan(1-f)*Tan(rlat2)
Corected:
Code: Select all
Protected U1.d= ATan((1-f)*Tan(rlat1))
Protected U2.d= ATan((1-f)*Tan(rlat2))
Bernd

Thanks for opening my eyes!

And here is the Procedure with Vincenty's modification:

Code: Select all

Procedure.d distVincenty2(lat1.d, lon1.d, lat2.d, lon2.d)
  Protected a.d=6378.137, b.d= 6356.752, f.d=(a-b)/a;1/298.257223563
  Protected rlat1.d=Radian(lat1)
  Protected rlat2.d=Radian(lat2)
  Protected rlon1.d=Radian(lon1)
  Protected rlon2.d=Radian(lon2)
  Protected L.d=Radian(lon2-lon1)
  Protected U1.d= ATan((1-f)*Tan(rlat1))
  Protected U2.d= ATan((1-f)*Tan(rlat2))
  Protected sinU1.d=Sin(U1)
  Protected cosU1.d=Cos(U1)
  Protected sinU2.d=Sin(U2)
  Protected cosU2.d=Cos(U2)
  
  Protected lambda.d=L, lambdaP.d, iterLimit.i=100
  
  Repeat
    Protected sinLambda.d=Sin(lambda)
    Protected cosLambda.d=Cos(lambda)
    Protected sinSigma.d=Sqr((cosU2*sinLambda)*(cosU2*sinLambda)+(cosU1*sinU2-sinU1*cosU2*cosLambda)*(cosU1*sinU2-sinU1*cosU2*cosLambda))
    If sinSigma=0
      ProcedureReturn 0
    EndIf
    Protected cosSigma.d=sinU1*sinU2+cosU1*cosU2*cosLambda
    Protected sigma.d=ATan2(cosSigma, sinSigma)
    Protected sinAlpha.d=cosU1*cosU2*sinLambda/sinSigma
    Protected cosSqAlpha.d=1-sinAlpha*sinAlpha
    Protected cos2SigmaM.d=cosSigma-2*sinU1*sinU2/cosSqAlpha
    If (IsNAN(cos2SigmaM))
      cos2SigmaM=0
    EndIf
    Protected C.d=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)))
    iterLimit-1
  Until Abs(lambda-lambdaP)<1e-12 Or iterLimit<1
  If iterLimit=0
    ProcedureReturn NaN()
  EndIf
  Protected uSq.d=cosSqAlpha*(a*a-b*b)/(b*b)
  Protected k1.d=(Sqr(1+uSq)-1)/(Sqr(1+uSq)+1) ; ---- Vincenty's modification ----
  Protected AA.d=(1+1/4*k1*k1)/(1-k1) ; ---- Vincenty's modification ----
  Protected BB.d=k1*(1-3/8*k1*k1) ; ---- Vincenty's modification ----
  Protected deltaSigma.d=BB*sinSigma*(cos2SigmaM+BB/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM) - BB/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)))
  Protected s.d=b*AA*(sigma-deltaSigma)
  
  
  ProcedureReturn s
EndProcedure

PureBasic Forums - English

Has anyone done the 'Haversine' formula?

Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?

Re: Has anyone done the 'Haversine' formula?