It is currently Fri Sep 20, 2019 11:45 pm

All times are UTC + 1 hour




Post new topic Reply to topic  [ 27 posts ]  Go to page 1, 2  Next
Author Message
 Post subject: Has anyone done the 'Haversine' formula?
PostPosted: Fri Oct 14, 2011 4:00 am 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Wed Sep 22, 2010 1:50 am
Posts: 796
Location: Bradenton, FL
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.


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Fri Oct 14, 2011 12:49 pm 
Offline
Addict
Addict

Joined: Sun Sep 07, 2008 12:45 pm
Posts: 4310
Location: Germany
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:
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:
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.
viewtopic.php?f=4&t=47855


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Fri Oct 14, 2011 5:20 pm 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Wed Oct 22, 2003 2:51 am
Posts: 743
Location: Canada
An optimized version

Code:
; <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
; 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


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Fri Oct 14, 2011 6:30 pm 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Wed Sep 22, 2010 1:50 am
Posts: 796
Location: Bradenton, FL
Thanks so much folks! I was sure someone had worked with it here.


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Sat Oct 15, 2011 2:40 pm 
Offline
PureBasic Expert
PureBasic Expert

Joined: Sun Aug 08, 2004 5:21 am
Posts: 3443
Location: Netherlands
The other function on the page (Spherical Law) also seems to work pretty well
Code:
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


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Sun Oct 16, 2011 12:42 am 
Offline
Addict
Addict

Joined: Wed Feb 28, 2007 9:13 am
Posts: 1037
Location: London
You could do that SphericalLaw procedure as a macro.

Code:
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

_________________
JACK WEBB: "Coding in C is like sculpting a statue using only sandpaper. You can do it, but the result wouldn't be any better. So why bother? Just use the right tools and get the job done."


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Sun Oct 16, 2011 1:00 am 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Wed Sep 22, 2010 1:50 am
Posts: 796
Location: Bradenton, FL
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.


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Sun Oct 16, 2011 6:40 am 
Offline
PureBasic Expert
PureBasic Expert

Joined: Sun Aug 08, 2004 5:21 am
Posts: 3443
Location: Netherlands
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:
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:
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.


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Sun Oct 16, 2011 3:33 pm 
Offline
Enthusiast
Enthusiast

Joined: Tue May 26, 2009 2:11 pm
Posts: 585
Just the right place to jump in. :D

Did anybody try to implement Vincenty formula?

I tried to convert this JavaScript but failed:

Code:
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
/* 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:
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?

_________________
Image


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Sun Oct 16, 2011 4:38 pm 
Offline
Addict
Addict

Joined: Sun Sep 07, 2008 12:45 pm
Posts: 4310
Location: Germany
And the error is .... here:
Code:
Protected U1.d= ATan(1-f)*Tan(rlat1)
Protected U2.d= ATan(1-f)*Tan(rlat2)

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

Bernd


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Sun Oct 16, 2011 5:38 pm 
Offline
Addict
Addict

Joined: Sun Sep 07, 2008 12:45 pm
Posts: 4310
Location: Germany
To make life easier:
viewtopic.php?f=12&t=47881

Bernd


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Sun Oct 16, 2011 6:26 pm 
Offline
Addict
Addict

Joined: Sun Sep 07, 2008 12:45 pm
Posts: 4310
Location: Germany
Since I have now so many input,
I combined all together:
Code:
#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.


Last edited by infratec on Mon Oct 17, 2011 7:24 am, edited 6 times in total.

Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Sun Oct 16, 2011 6:33 pm 
Offline
PureBasic Expert
PureBasic Expert

Joined: Sun Aug 08, 2004 5:21 am
Posts: 3443
Location: Netherlands
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.


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Sun Oct 16, 2011 6:41 pm 
Offline
Addict
Addict

Joined: Sun Sep 07, 2008 12:45 pm
Posts: 4310
Location: Germany
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.


Top
 Profile  
Reply with quote  
 Post subject: Re: Has anyone done the 'Haversine' formula?
PostPosted: Sun Oct 16, 2011 6:58 pm 
Offline
Enthusiast
Enthusiast

Joined: Tue May 26, 2009 2:11 pm
Posts: 585
infratec wrote:
And the error is .... here:
Code:
Protected U1.d= ATan(1-f)*Tan(rlat1)
Protected U2.d= ATan(1-f)*Tan(rlat2)

Corected:
Code:
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:
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

_________________
Image


Top
 Profile  
Reply with quote  
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 27 posts ]  Go to page 1, 2  Next

All times are UTC + 1 hour


Who is online

Users browsing this forum: No registered users and 11 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum

Search for:
Jump to:  

 


Powered by phpBB © 2008 phpBB Group
subSilver+ theme by Canver Software, sponsor Sanal Modifiye