# PureBasic Forum

 It is currently Fri Mar 22, 2019 1:33 am

 All times are UTC + 1 hour

 Page 1 of 2 [ 27 posts ] Go to page 1, 2  Next
 Print view Previous topic | Next topic
Author Message
 Post subject: Has anyone done the 'Haversine' formula?Posted: Fri Oct 14, 2011 4:00 am
 Enthusiast

Joined: Wed Sep 22, 2010 1:50 am
Posts: 786
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

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Fri Oct 14, 2011 12:49 pm
 Addict

Joined: Sun Sep 07, 2008 12:45 pm
Posts: 4115
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

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Fri Oct 14, 2011 5:20 pm
 Enthusiast

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

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Fri Oct 14, 2011 6:30 pm
 Enthusiast

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

Top

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Sat Oct 15, 2011 2:40 pm
 PureBasic Expert

Joined: Sun Aug 08, 2004 5:21 am
Posts: 3319
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

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Sun Oct 16, 2011 12:42 am
 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

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Sun Oct 16, 2011 1:00 am
 Enthusiast

Joined: Wed Sep 22, 2010 1:50 am
Posts: 786
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

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Sun Oct 16, 2011 6:40 am
 PureBasic Expert

Joined: Sun Aug 08, 2004 5:21 am
Posts: 3319
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

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Sun Oct 16, 2011 3:33 pm
 Enthusiast

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

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?

_________________

Top

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Sun Oct 16, 2011 4:38 pm
 Addict

Joined: Sun Sep 07, 2008 12:45 pm
Posts: 4115
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

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Sun Oct 16, 2011 5:38 pm
 Addict

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

Bernd

Top

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Sun Oct 16, 2011 6:26 pm
 Addict

Joined: Sun Sep 07, 2008 12:45 pm
Posts: 4115
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

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Sun Oct 16, 2011 6:33 pm
 PureBasic Expert

Joined: Sun Aug 08, 2004 5:21 am
Posts: 3319
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

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Sun Oct 16, 2011 6:41 pm
 Addict

Joined: Sun Sep 07, 2008 12:45 pm
Posts: 4115
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

 Post subject: Re: Has anyone done the 'Haversine' formula?Posted: Sun Oct 16, 2011 6:58 pm
 Enthusiast

Joined: Tue May 26, 2009 2:11 pm
Posts: 536
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

_________________

Top

 Display posts from previous: All posts1 day7 days2 weeks1 month3 months6 months1 year Sort by AuthorPost timeSubject AscendingDescending
 Page 1 of 2 [ 27 posts ] Go to page 1, 2  Next

 All times are UTC + 1 hour

#### Who is online

Users browsing this forum: Blue, falsam and 15 guests

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

Search for:
 Jump to:  Select a forum ------------------ PureBasic    Coding Questions    Game Programming    3D Programming    Assembly Programming    The PureBasic Editor    The PureBasic Form Designer    General Discussion    Feature Requests and Wishlists    Tricks 'n' Tips Bug Reports    Bugs - Windows    Bugs - Linux    Bugs - Mac OSX    Bugs - Documentation OS Specific    AmigaOS    Linux    Windows    Mac OSX Miscellaneous    Announcement    Off Topic Showcase    Applications - Feedback and Discussion    PureFORM & JaPBe    TailBite

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