Great Circle Formula in VB
Paul Hazen
I am trying to use a version of the great circle formula in a vb 6.0 app.
Here is the formula:
cos D = (sin p sin a) + (cos p cos a cos |dl|)
p and a are the latitudes of P and A
|dl| is the absolute value of the difference in longitude between P and A
Here is an example of the formula in action:
1. Convert to decimal degrees:
Paris: 48.87° N, 02.33° E
Austin: 30.27° N, 97.74° W
2. cos D = [sin(48.87) * sin(30.27)] + [cos(48.87) * cos(30.27) *
cos(|-97.74 - 2.33|)]
3. cos D = [0.753 * 0.504] + [0.658 * 0.864 * -0.174] = 0.281
4. D = cos-1 (0.281) = 73.68°
5. D = 73.68° * 111.23 kilometers/degree = 8,195.44 kilometers (5,092.03
miles)
Anyone have any ideas how I can get this coded into VB? I have tried several
different methods, none work. I get hung up on step four. The cos-1.
Any help would be appreciated!
Paul Hazen
MCP
cooch
is my old QB function for acos:
FUNCTION Acos# (value#)
absval# = ABS(value#)
ang# = 0
delta# = 1
GOTO ang
delta: delta# = delta# / 10
ang: ang# = ang# + delta#
IF ang# > pi# THEN
actval# = 0
ELSE
actval# = COS(ang#)
END IF
IF actval# > absval# THEN GOTO ang
IF delta# < 0.00001 THEN GOTO GotAng
ang# = ang# - delta#
GOTO delta
GotAng:
IF value# < 0 THEN ang# = pi# - ang#
Acos# = ang#
END FUNCTION
Hope it helps
Cooch (a newcomer to this forum)
In article <#yl3khrTAHA.265@cppssbbsa05>,
"Paul Hazen" <pha...@ka.net> wrote:
> .... I get hung up on step four .. cos-1.
Sent via Deja.com http://www.deja.com/
Before you buy.
Ken James
has formulae for lots of functions that ought to be in the standard VB
library, but aren't.
Inverse Cosine is defined there as Arccos(X) =Atn((-X/Sqr(-X*X+1))+2*Atn(1),
which is correct, but needs protection from values of X very close to 1.0
which may result in Divide By Zero errors.
Your Great Circle formula may be great for tennis balls, but the Earth is
not spherical! You might be horrified at the real formula's complexity.
Ken James
To add insult to injury, that 2*Atn(1) is Pi/2, and Pi ought to be a
predefined constant too!
Paul Hazen <pha...@ka.net> wrote in message
news:#yl3khrTAHA.265@cppssbbsa05...
Lenny Abbey
that it makes your code "spaghetti code" in some sort of de facto fashion.
Lenny
"cooch" <couc...@onthenet.com.au> wrote in message
news:8utdj6$j5t$1...@nnrp1.deja.com...
Roger Hamlett
> is my old QB function for acos:
>
> FUNCTION Acos# (value#)
> absval# = ABS(value#)
> ang# = 0
> delta# = 1
> GOTO ang
> delta: delta# = delta# / 10
> ang: ang# = ang# + delta#
> IF ang# > pi# THEN
> actval# = 0
> ELSE
> actval# = COS(ang#)
> END IF
> IF actval# > absval# THEN GOTO ang
> IF delta# < 0.00001 THEN GOTO GotAng
> ang# = ang# - delta#
> GOTO delta
> GotAng:
> IF value# < 0 THEN ang# = pi# - ang#
> Acos# = ang#
> END FUNCTION
Sledgehammer... :-)
The reason that ACos is not coded, is that it is not needed. You have the
Atn function.
So the statement:
IIf(x<>0,Atn(Sqr(1#-x*x)/x)),pi/2)
will return the aCos of 'x'.
Best Wishes
Roger Hamlett
"Ken James" <bristec_...@compuserve.com> wrote in message
news:8uthfm$202$1...@newsg3.svr.pol.co.uk...
Sent: Wednesday, November 15, 2000 8:18 AM
Subject: Re: Great Circle Formula in VB
> Appendix D of Microsoft's VB6 Language Reference (Derived Math
Functions)
> has formulae for lots of functions that ought to be in the standard VB
> library, but aren't.
>
> Inverse Cosine is defined there as Arccos(X)
>=Atn((-X/Sqr(-X*X+1))+2*Atn(1),
> which is correct, but needs protection from values of X very close to
>1.0
> which may result in Divide By Zero errors.
Much better, since this works for all the circle segments (what I posted
would only work for the first segment). :-)
I was posting quickly, and did the formula 'out of my head' (which is
where I normally am, some people would say... :-)
Does the language evaluate 2*Atn(1) quicker than using a constant for this
value?. Anyone know?.
so:
Const Piby2 = 1.570796327
IIf((Abs(x) < 0.9999), Atn(-x / Sqr(-x * x + 1)) + Piby2, 0)
Would seem to be a 'safe', and quick answer, or (if 2*Atn(1) is faster
than a constant, the same but using this instead of the constant.
> Your Great Circle formula may be great for tennis balls, but the Earth
is
> not spherical! You might be horrified at the real formula's complexity.
Yes, but the error is not great over normal distances. Depends what you
actually want to do...
> Ken James
>
> To add insult to injury, that 2*Atn(1) is Pi/2, and Pi ought to be a
> predefined constant too!
The same thought as me. I am suspicious that the internal fpu, has some
things defined internally as constants, so 'Atn(1)', may simply return
Pi/2, without performing any maths at all. This would be faster than
transferring a float to the fpu, and even doubling this, may be the
fastest solution. It would involve doing a quick 'benchmark', performing a
few thousand conversions to find out.
Just gone and done this:-)
Can confirm that using the 'constant', is _fractionally_ faster. The
difference is tiny, and I suspect that 'Atn(1)', is actually returned
faster than the constant, but by the time the mutiplication is done, the
'edge' swings the other way. The difference is between 0.202, and 0.23
seconds, performing 100000 conversions (on my machine here), so it isn't
exactly 'large'.
Best Wishes
Bertie Wooster
> Eeek!.
> Sledgehammer... :-)
> The reason that ACos is not coded, is that it is not needed. You have the
> Atn function.
> So the statement:
> IIf(x<>0,Atn(Sqr(1#-x*x)/x)),pi/2)
> will return the aCos of 'x'.
Actually it will blow up because IIF evaluates both expressions anyway. You
need IF/THEN/ELSE.
Bertie
Roger Hamlett
"Bertie Wooster" <ber...@thedrones.club.uk> wrote in message
news:eHOR6ov...@cppssbbsa02.microsoft.com...
Best Wishes
Norm Cook
Option Explicit
Private PI#
Private Sub Form_Load()
Dim PLa#, ALa# 'Latitudes
Dim PLo#, ALo# 'Longitudes
Dim DLo# 'Delta Longitude
Dim DistKM# 'Kilometers
Dim DistMi# 'Miles
Dim DistDeg#, CosD# 'intermediate variables
PI = 4 * Atn(1)
PLo = 2.33: ALo = -97.74
DLo = Rads(Abs(PLo - ALo))
PLa = Rads(48.87)
ALa = Rads(30.27)
CosD = Sin(PLa) * Sin(ALa) + Cos(PLa) * Cos(ALa) * Cos(DLo)
DistDeg = Deg(ACos(CosD))
DistKM = Kilometers(DistDeg)
DistMi = Miles(DistKM)
Debug.Print DistKM, DistMi
End Sub
Private Function Rads(ByVal Deg#) As Double
Rads = Deg * PI / 180
End Function
Private Function Deg(ByVal Rads#) As Double
Deg = Rads * 180 / PI
End Function
Private Function Kilometers(ByVal Deg#) As Double
Kilometers = Deg * 111.23
End Function
Private Function Miles(ByVal km#) As Double
Miles = km / 1.609
End Function
Private Function ACos(value As Double) As Double
If Abs(value) <> 1 Then
ACos = 1.5707963267949 - Atn(value / Sqr(1 - value * value))
Else
ACos = 3.14159265358979 * Sgn(value)
End If
End Function
Ken James
avoid DivideByZero type errors!
2. Paris to Austin, Tx is not a "normal" distance, so I would expect the
error to be several kiliometres.
3. Anyone who likes complicated equations cold start by looking at
www.wgs84.com
Ken James
Roger Hamlett <ro...@ttelmah.demon.co.uk> wrote in message
news:PtuQ5.1248$q05....@news11-gui.server.ntli.net...
Tim Arheit
<ro...@ttelmah.demon.co.uk> wrote:
>
>"Ken James" <bristec_...@compuserve.com> wrote in message
>news:8uthfm$202$1...@newsg3.svr.pol.co.uk...
>Sent: Wednesday, November 15, 2000 8:18 AM
>Subject: Re: Great Circle Formula in VB
>
>
>> Appendix D of Microsoft's VB6 Language Reference (Derived Math
>Functions)
>> has formulae for lots of functions that ought to be in the standard VB
>> library, but aren't.
>>
>> Inverse Cosine is defined there as Arccos(X)
>>=Atn((-X/Sqr(-X*X+1))+2*Atn(1),
>> which is correct, but needs protection from values of X very close to
>>1.0
>> which may result in Divide By Zero errors.
>Much better, since this works for all the circle segments (what I posted
>would only work for the first segment). :-)
>I was posting quickly, and did the formula 'out of my head' (which is
>where I normally am, some people would say... :-)
>Does the language evaluate 2*Atn(1) quicker than using a constant for this
>value?. Anyone know?.
A constant would be quicker than a multiplication and a function call
(regardless of what the function does). Just be sure to carry it out
to the maximum number of digits the variables you are using support
(ie. single vs double)
-Tim
Roger Hamlett
"Norm Cook" <norm...@arn.net> wrote in message
news:e9p7ZswTAHA.278@cppssbbsa03...
context of VB, since it evaluates both states whatever the condition in
VB.
However your 'formula' is wrong on the 'else' state, since it gives
Acos(1) as PI (should be Pi/2)...
Also you have used the formula I posted, which is wrong outside the first
quadrant.
Private Function ACos(value As Double) As Double
If Abs(value) <> 1 Then
ACos = 1.5707963267949 + Atn(-value / Sqr(1 - value * value))
Else
ACos = 1.5707963267949 * (1-value)
End If
End Function
I think this is 'right', over the entire -1 to 1 range... :-)
Fun Huh!.
Matt Sargent
Paste the following into a bas module:
Option Explicit
Enum gsType
gsMiles = 1
gsKilometers = 2
End Enum
'The Polar radius is: 3949.99 miles, 6356.8927 kilometers
'The Equatorial radius is: 3963.34 miles, 6378.1602 kilometers
Enum Radius
Miles = 3956.665 ' Average of Polar and Equatorial
Kilometers = 6367.52645 ' Average of Polar and Equatorial
End Enum
Private Function ArcCos(X As Double) As Double
ArcCos = Atn(-X / Sqr(-X * X + 1)) + 2 * Atn(1)
End Function
Public Function GreatCircle(Lat1 As Double, Lon1 As Double, _
Lat2 As Double, Lon2 As Double, Mode As gsType) As Long
' This is the simple Great Circle calculation, which assumes a
' spherical earth, not accounting for the equatorial bulge. An
' average of the polar and equatorial radius of the earth is used.
' Mins, and secs must be converted to degree decimals:
' (Mins * 60 + Secs) / 3600
Dim dWork As Double
' Pre-calc radians for speed.
Lat1 = Rad(Lat1)
Lat2 = Rad(Lat2)
dWork = ArcCos( _
(Sin(Lat1) * Sin(Lat2)) + _
(Cos(Lat1) * Cos(Lat2) * Cos(Rad(Lon2 - Lon1))) _
)
If Mode = gsMiles Then
GreatCircle = dWork * Radius.Miles
ElseIf Mode = gsKilometers Then
GreatCircle = dWork * Radius.Kilometers
End If
End Function
Private Function Rad(X As Double) As Double
' Radian = X / 180 * Pi
' Pi = 4 * Atn(1)
Rad = X / 45 * Atn(1)
End Function
On Tue, 14 Nov 2000 22:59:08 -0500, "Paul Hazen" <pha...@ka.net>
wrote:
>Hello all,
>I am trying to use a version of the great circle formula in a vb 6.0 app.
>Here is the formula:
>cos D = (sin p sin a) + (cos p cos a cos |dl|)
>p and a are the latitudes of P and A
>|dl| is the absolute value of the difference in longitude between P and A
>
>Here is an example of the formula in action:
>1. Convert to decimal degrees:
>Paris: 48.87° N, 02.33° E
>Austin: 30.27° N, 97.74° W
>
>2. cos D = [sin(48.87) * sin(30.27)] + [cos(48.87) * cos(30.27) *
>cos(|-97.74 - 2.33|)]
>3. cos D = [0.753 * 0.504] + [0.658 * 0.864 * -0.174] = 0.281
>4. D = cos-1 (0.281) = 73.68°
>5. D = 73.68° * 111.23 kilometers/degree = 8,195.44 kilometers (5,092.03
>miles)
>
>Anyone have any ideas how I can get this coded into VB? I have tried several
>different methods, none work. I get hung up on step four. The cos-1.
>Any help would be appreciated!
>
>Paul Hazen
>MCP
>
>
Jeff Johnson
"Lenny Abbey" <LAb...@mindspring.com> wrote in message
news:8uti9q$nf4$1...@slb6.atl.mindspring.net...
> Haven't you heard? They jump all over you here if you use GOTO. It seems
> that it makes your code "spaghetti code" in some sort of de facto fashion.
I'm a vocal defender of GoTo when it makes sense, but that was just plain
scary....
Colin Young
fair, we were warned it was QuickBasic code...
Colin
"Jeff Johnson" <pawp...@com.geocities> wrote in message
news:euEMf#3TAHA.278@cppssbbsa03...
Lenny Abbey
code much easier to understand.
Lenny
"Colin Young" <x...@nospam.com> wrote in message
news:#ebekp9TAHA.193@cppssbbsa05...
Bertie Wooster
> I always use type declaration characters in my variable names. It makes
the
> code much easier to understand.
You are one sick puppy. Get help immediately.
Bertie
Lenny Abbey
Lenny
"Bertie Wooster" <ber...@thedrones.club.uk> wrote in message
news:OvmjACIUAHA.243@cppssbbsa03...