Has anyone done the 'Haversine' formula?
- RichAlgeni
- Addict
- Posts: 914
- Joined: Wed Sep 22, 2010 1:50 am
- Location: Bradenton, FL
Has anyone done the 'Haversine' formula?
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.
http://www.movable-type.co.uk/scripts/latlong.html
I'm sure it has to be fairly simple using PB.
Re: Has anyone done the 'Haversine' formula?
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:But it gives wrong Results.
This one is working:
Bernd
P.S.: I opened a bug thread for that.
http://www.purebasic.fr/english/viewtop ... =4&t=47855
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
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
P.S.: I opened a bug thread for that.
http://www.purebasic.fr/english/viewtop ... =4&t=47855
Re: Has anyone done the 'Haversine' formula?
An optimized version
This version can be used to calculate distance on any planet or moon.
Edit :
DMS_To_DD() conversion function added.
Best regards
Guimauve
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 <<<<<
; <<<<<<<<<<<<<<<<<<<<<<<
Edit :
DMS_To_DD() conversion function added.
Best regards
Guimauve
- RichAlgeni
- Addict
- Posts: 914
- Joined: Wed Sep 22, 2010 1:50 am
- Location: Bradenton, FL
Re: Has anyone done the 'Haversine' formula?
Thanks so much folks! I was sure someone had worked with it here.
Re: Has anyone done the 'Haversine' formula?
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
-
- Addict
- Posts: 1233
- Joined: Wed Feb 28, 2007 9:13 am
- Location: London
Re: Has anyone done the 'Haversine' formula?
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."
- RichAlgeni
- Addict
- Posts: 914
- Joined: Wed Sep 22, 2010 1:50 am
- Location: Bradenton, FL
Re: Has anyone done the 'Haversine' formula?
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.
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.
Re: Has anyone done the 'Haversine' formula?
I didn't check for other places.RichAlgeni wrote:The spherical law of cosines seems to be less accurate than Haversine the further you get from the equator.
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
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
Re: Has anyone done the 'Haversine' formula?
Just the right place to jump in.
Did anybody try to implement Vincenty formula?
I tried to convert this JavaScript but failed:
My attempt:
I noticed, that I get an N-S error, the E-W part seems to be ok.
Where is my error?
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() };
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
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
Where is my error?
Re: Has anyone done the 'Haversine' formula?
And the error is .... here:
Corected:
Bernd
Code: Select all
Protected U1.d= ATan(1-f)*Tan(rlat1)
Protected U2.d= ATan(1-f)*Tan(rlat2)
Code: Select all
Protected U1.d= ATan((1-f)*Tan(rlat1))
Protected U2.d= ATan((1-f)*Tan(rlat2))
Re: Has anyone done the 'Haversine' formula?
Since I have now so many input,
I combined all together:
Have fun,
Bernd
P.S.: I don't know why Vincenty gives a bit other results.
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
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.
Re: Has anyone done the 'Haversine' formula?
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.
It's a bit flattened at the poles. Haversine and SphericalLaw ignore this so they are a bit less accurate.
Re: Has anyone done the 'Haversine' formula?
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.
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.
Re: Has anyone done the 'Haversine' formula?
Thanks for opening my eyes!infratec wrote:And the error is .... here:Corected:Code: Select all
Protected U1.d= ATan(1-f)*Tan(rlat1) Protected U2.d= ATan(1-f)*Tan(rlat2)
BerndCode: Select all
Protected U1.d= ATan((1-f)*Tan(rlat1)) Protected U2.d= ATan((1-f)*Tan(rlat2))
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