Has anyone done the 'Haversine' formula?

Just starting out? Need help? Post your questions and find answers here.
User avatar
RichAlgeni
Addict
Addict
Posts: 914
Joined: Wed Sep 22, 2010 1:50 am
Location: Bradenton, FL

Has anyone done the 'Haversine' formula?

Post by RichAlgeni »

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

Re: Has anyone done the 'Haversine' formula?

Post by infratec »

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
User avatar
Guimauve
Enthusiast
Enthusiast
Posts: 742
Joined: Wed Oct 22, 2003 2:51 am
Location: Canada

Re: Has anyone done the 'Haversine' formula?

Post by Guimauve »

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
User avatar
RichAlgeni
Addict
Addict
Posts: 914
Joined: Wed Sep 22, 2010 1:50 am
Location: Bradenton, FL

Re: Has anyone done the 'Haversine' formula?

Post by RichAlgeni »

Thanks so much folks! I was sure someone had worked with it here.
wilbert
PureBasic Expert
PureBasic Expert
Posts: 3870
Joined: Sun Aug 08, 2004 5:21 am
Location: Netherlands

Re: Has anyone done the 'Haversine' formula?

Post by wilbert »

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
Seymour Clufley
Addict
Addict
Posts: 1233
Joined: Wed Feb 28, 2007 9:13 am
Location: London

Re: Has anyone done the 'Haversine' formula?

Post by Seymour Clufley »

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
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."
User avatar
RichAlgeni
Addict
Addict
Posts: 914
Joined: Wed Sep 22, 2010 1:50 am
Location: Bradenton, FL

Re: Has anyone done the 'Haversine' formula?

Post by RichAlgeni »

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

Re: Has anyone done the 'Haversine' formula?

Post by wilbert »

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.
User avatar
Lord
Addict
Addict
Posts: 847
Joined: Tue May 26, 2009 2:11 pm

Re: Has anyone done the 'Haversine' formula?

Post by Lord »

Just the right place to jump in. :D

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?
Image
infratec
Always Here
Always Here
Posts: 6817
Joined: Sun Sep 07, 2008 12:45 pm
Location: Germany

Re: Has anyone done the 'Haversine' formula?

Post by infratec »

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
infratec
Always Here
Always Here
Posts: 6817
Joined: Sun Sep 07, 2008 12:45 pm
Location: Germany

Re: Has anyone done the 'Haversine' formula?

Post by infratec »

infratec
Always Here
Always Here
Posts: 6817
Joined: Sun Sep 07, 2008 12:45 pm
Location: Germany

Re: Has anyone done the 'Haversine' formula?

Post by infratec »

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.
Last edited by infratec on Mon Oct 17, 2011 7:24 am, edited 6 times in total.
wilbert
PureBasic Expert
PureBasic Expert
Posts: 3870
Joined: Sun Aug 08, 2004 5:21 am
Location: Netherlands

Re: Has anyone done the 'Haversine' formula?

Post by wilbert »

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

Re: Has anyone done the 'Haversine' formula?

Post by infratec »

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.
User avatar
Lord
Addict
Addict
Posts: 847
Joined: Tue May 26, 2009 2:11 pm

Re: Has anyone done the 'Haversine' formula?

Post by Lord »

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
Image
Post Reply