x
• Status: Solved
• Priority: Medium
• Security: Public
• Views: 802

VB Solution to return the angle, in degrees, from a given lat, lon to current position.

Hi,

Basically the title says it all. What I need is a routine which will return the angle, in degrees, from a given point (given as Lat,Lon) to my current position (also given as Lat,Lon), using VB6.

Thanks

Dave
0
dmarch100
• 5
• 5
• 3
2 Solutions

Commented:
Here are 2 Functions. 1 for Bearing and another for GeoDistBetween.

'To call Bearing:

Heading = Bearing(LatDest - LatOrig, LongDest - LongOrig)

' ----------------------------------------------------------------------
' Inverse longitudinal tangent function to return angles in degrees.
' This function uses dual arguments and returns an angle in the range
' from 0 to 360 degrees relative to True North.

Public Function Bearing(ArgY, ArgX) As Single
'  Check for special coordinate zero value cases
On Error GoTo Bearing_Err
100   If ArgY = 0 And ArgX >= 0 Then
102      Bearing = 90
104   ElseIf ArgY > 0 And ArgX = 0 Then
106      Bearing = 180
108   ElseIf ArgY = 0 And ArgX < 0 Then
110      Bearing = 270
112   ElseIf ArgY < 0 And ArgX = 0 Then
114      Bearing = 0
Else
116      Bearing = 90 + 45 * Atn(ArgY / ArgX) / Atn1
118      If ArgX < 0 Then
120         Bearing = Bearing + 180
End If
122      If Bearing < 0 Then
124         Bearing = Bearing + 360
End If
End If

'<EhFooter>
Exit Function

Bearing_Err:
Select Case ErrMsgBox("Cyberime2002.modAstro.Bearing")
Case vbAbort
Screen.MousePointer = vbDefault
Exit Function
Case vbRetry
Resume
Case vbIgnore
Resume Next
End Select
'</EhFooter>
End Function

'--------

Private Function GeoDistBetween(Long1 As Single, Lat1 As Single, Long2 As Single, Lat2 As Single) As Single
' V2.0
'
' Compute the gedesic distance between two Earth surface
' coordinates to an accuracy of about ±50 meters.  To get
' this degree of accuracy, this function takes into account
' the spheroidal flattening factor of the Earth rather than
' assuming that the Earth is a perfect sphere.
'
' The arguments (Lng1, Lat1, Lng2, Lat2) are the longitudes
' and latitudes of the two locations expressed in degrees.
On Error GoTo GeoDistBetween_Err

' Positive longitude is west and negative is east.
'
' ===========================================================
' Function written by Jay Tanner - Based on the mathematical
' method of H. Andoyer.
'
' References:
'
' International Earth Rotation Service,
' Annual report for 1996 (Observatorie de Paris, 1997)
'
' Annuaire du Bureau des Longitudes pour 1950 (Paris) pg.145
'
' ==========================================================
' Flattening factor of the earth geoid
Const ff As Double = 1 / 298.257

' Auxiliary working variables
Dim c  As Double
Dim D  As Double
Dim F  As Double
Dim g  As Double
Dim H1 As Double
Dim H2 As Double
Dim L  As Double
Dim O  As Double
Dim r  As Double
Dim s  As Double
Dim W1 As Double
Dim W2 As Double
Dim SG As Double
Dim cG As Double
Dim SF As Double
Dim CF As Double
Dim sl As Double
Dim CL As Double

' Compute auxiliary angles
100   F = (Lat1 + Lat2) / 2
102   g = (Lat1 - Lat2) / 2
104   L = (Long1 - Long2) / 2
106   If g + L = 0 Then
108      MsgBox "Warning." & vbCrLf & "Source and Destin are equal."
Exit Function
End If
' Compute sines and cosines of auxiliary angles
110   SG = Sin(g * Atn(1) / 45)
112   cG = Cos(g * Atn(1) / 45)
114   SF = Sin(F * Atn(1) / 45)
116   CF = Cos(F * Atn(1) / 45)
118   sl = Sin(L * Atn(1) / 45)
120   CL = Cos(L * Atn(1) / 45)

122   W1 = SG * CL: W1 = W1 * W1
124   W2 = CF * sl: W2 = W2 * W2
126   s = W1 + W2

128   W1 = cG * CL: W1 = W1 * W1
130   W2 = SF * sl: W2 = W2 * W2
132   c = W1 + W2

134   O = Atn(Sqr(s / c))
136   r = Sqr(s * c) / O
138   D = 2 * O * 6378.14

140   H1 = (3 * r - 1) / (2 * c)
142   H2 = (3 * r + 1) / (2 * s)

' Compute the angle between the points on a synthetic
' geodesic sphereoid connecting the two points.
144   W1 = SF * cG: W1 = W1 * W1 * H1 * ff + 1
146   W2 = CF * SG: W2 = W2 * W2 * H2 * ff

' Return the distance between the given locations in
' Kilometers
148   GeoDistBetween = D * (W1 - W2)

'<EhFooter>
Exit Function

GeoDistBetween_Err:
Select Case ErrMsgBox("Cyberime2002.frmWorld.GeoDistBetween")
Case vbAbort
Screen.MousePointer = vbDefault
Exit Function
Case vbRetry
Resume
Case vbIgnore
Resume Next
End Select
'</EhFooter>
End Function

Public Function ErrMsgBox(Msg As String) As Integer
'*****WARNING*****
'Do not put error handling in this sub.
'Doing so will clear the Err Object and Erl

100   Msg = "Error Number " & Err.Number & vbCrLf & _
Err.Description & vbCrLf & _
"in " & Msg & vbCrLf & _
"at line " & Erl
102   ErrMsgBox = MsgBox(Msg, vbAbortRetryIgnore + vbCritical, "Error details")

End Function
0

Author Commented:
Hi danaseaman

Thanks for your answer but how does the 'Bearing' function relate to my current position in any way? I need to know the angle in degrees from a given point relative to my current location

Dave
0

Commented:
You call it with 4 parameters. Latitude and Longitude of both source and destination.
It will return the angle relative to the source Lat/Long.

Heading = Bearing(LatDest - LatOrig, LongDest - LongOrig)
0

Author Commented:
Hi danaseaman,

Sorry for not noting the call procedure as clearly shown in your original post.

However, alsthough on the right track it doesn't actually work. If we concentrate on the Bearing function only at line 116 the Atn1 variable is not defined or used and subsequently kicks out a Division by Zero error. With that variable removed from the equation it works but obviously returns an incorrect radial.  My test data is Glasgow airport as the current position and Heathrow as the destination....

Heading = Bearing(51.28 - 55.52, 0.27 - 4.25)

Dave
0

Commented:
Oops, forgot the constant.

Public Const Atn1     As Double = PI / 4

0

Commented:
Public Const PI       As Double = 3.14159265358979
Public Const PI2      As Double = PI * 2
Public Const Atn1     As Double = PI / 4
0

Commented:
Which bearing do you need - rhumb line or great circle?
0

Commented:
Public Enum bTypes
RhumbLine
GreatCircle
End Type

Public Type POSITION
lat As Double
lon As Double
End Type

Const PI As Double = 3.14159265358979
Const e As Double = 0.081813333

Public Function Azimuth(ptFrom As POSITION, ptTo As POSITION, Optional bType As bTypes = RhumbLine) As Double
If bType = GreatCircle Then
Azimuth = AzimuthGC(ptFrom, ptTo)
Else
Azimuth = AzimuthRL(ptFrom, ptTo)
End If
End Function

Public Function Distance(ptFrom As POSITION, ptTo As POSITION, Optional bType As bTypes = RhumbLine) As Double
If bType = GreatCircle Then
Distance = DistanceGC(ptFrom, ptTo)
Else
Distance = DistanceRL(ptFrom, ptTo)
End If
End Function

' Great Circle staff
Private Function AzimuthGC(ptFrom As POSITION, ptTo As POSITION) As Double
AzimuthGC = NormalizeAngle(Deg(Atan2(Sind(ptTo.lon - ptFrom.lon), Tand(ptTo.lat) * Cosd(ptFrom.lat) - Sind(ptFrom.lat) * Cosd(ptTo.lon - ptFrom.lon))))
End Function

Private Function DistanceGC(ptFrom As POSITION, ptTo As POSITION) As Double
DistanceGC = Abs(Arccosd(Sind(ptFrom.lat) * Sind(ptTo.lat) + Cosd(ptFrom.lat) * Cosd(ptTo.lat) * Cosd(ptTo.lon - ptFrom.lon)) * 60)
End Function

' Rhumb Line staff
Private Function AzimuthRL(ptFrom As POSITION, ptTo As POSITION) As Double
Dim dL As Double, dR As Double
dL = ptTo.lon - ptFrom.lon
If dL > 180 Then
dL = dL - 360
ElseIf dL < -180 Then
dL = dL + 360
End If
dR = MCH(ptTo.lat) - MCH(ptFrom.lat)
AzimuthRL = NormalizeAngle(Deg(Atan2(dL * 60, dR)))
End Function

Private Function DistanceRL(ptFrom As POSITION, ptTo As POSITION) As Double
Dim dF As Double, dL As Double, dW As Double, fAver As Double
dF = ptTo.lat - ptFrom.lat
dL = ptTo.lon - ptFrom.lon
fAver = (ptTo.lat + ptFrom.lat) / 2
dW = dL * Cosd(fAver)
DistanceRL = Sqr(dF * dF + dW * dW)
End Function

Private Function Rad(ByVal x As Double) As Double
Rad = x * PI / 180
End Function

Private Function Tand(ByVal x As Double) As Double
End Function

Private Function Deg(ByVal x As Double) As Double
Deg = x * 180 / PI
End Function

'Sin with argument in degrees
Private Function Sind(ByVal x As Double) As Double
End Function

'Cos with argument in degrees
Private Function Cosd(ByVal x As Double) As Double
End Function

Private Function Atan2(ByVal Y As Double, ByVal x As Double) As Double
If x = 0 Then
Atan2 = IIf(Y < 0, 3 * PI / 2, PI / 2)
Else
Atan2 = Atn(Y / x)
If x < 0 Then
If Y < 0 Then Atan2 = Atan2 + PI Else Atan2 = Atan2 - PI
End If
End If
End Function

'ArcCos in degrees
Private Function Arccosd(ByVal x As Double) As Double
On Error Resume Next
If x = 1 Then
Arccosd = 0
ElseIf x = -1 Then
Arccosd = 180
Else
Arccosd = Deg(Atan2(-x, Sqr(1 - x * x))) + 90
End If
On Error GoTo 0
End Function

'ArcSin in degrees
Private Function Arcsind(ByVal x As Double) As Double
If x = 1 Then
Arcsind = 90
ElseIf x = -1 Then
Arcsind = 270
Else
Arcsind = Deg(Atan2(x, Sqr(-x * x + 1)))
End If
End Function

Private Function NormalizeAngle(x As Double) As Double
NormalizeAngle = x
If x < 0 Then NormalizeAngle = x + 360
If x > 360 Then NormalizeAngle = x - 360
End Function

'Meridional parts
Private Function MCH(ByVal x As Double) As Double
Dim a As Double
a = 45 * 60 / Atn(1)
On Error Resume Next
MCH = a * Log(Tan(Rad(45) + Rad(x) / 2) * ((1 - e * Sind(x)) / (1 + e * Sind(x))) ^ (e / 2))
End Function

'=====Using=========
Private Sub Command1_Click()
Dim ptCurrent As POSITION, ptDestination As POSITION
'Glasgow
ptCurrent.lat = 55 + 50 / 60
ptCurrent.lon = -1 * (4 + 15 / 60)
'Heathrow
ptDestination.lat = 51 + 28 / 60 + 11 / 3600
ptDestination.lon = -1 * (0 + 27 / 60 + 5 / 3600)
'Calculation
Dim A_GC As Double, A_RL As Double
Dim D_GC As Double, D_RL As Double
A_GC = Azimuth(ptCurrent, ptDestination, GreatCircle)
A_RL = Azimuth(ptCurrent, ptDestination, RhumbLine)
D_GC = Distance(ptCurrent, ptDestination, GreatCircle)
D_RL = Distance(ptCurrent, ptDestination, RhumbLine)
Dim sMsg As String
sMsg = "From Glasgow to Heathrow:" & vbCrLf & "   Rhumb Line:" & vbCrLf
sMsg = sMsg & "Azimuth: " & Format(A_RL, "000.00 deg.") & vbCrLf & "Distance: " & Format(D_RL, "0000.00 nm.") & vbCrLf
sMsg = sMsg & "    Great Gircle: " & vbCrLf
sMsg = sMsg & "Azimuth: " & Format(A_GC, "000.00 deg.") & vbCrLf & "Distance: " & Format(D_GC, "0000.00 nm.")
MsgBox sMsg, vbInformation
End Sub
0

Commented:
Oops, sorry,
DistanceRL = Sqr(dF * dF + dW * dW) 'is in degrees
'Should be
DistanceRL = 60 * Sqr(dF * dF + dW * dW) 'in nautical miles

0

Commented:
PS
You can check it here: http://www.movable-type.co.uk/scripts/LatLong.html
Note that small difference in Rhumb Line calculations between above page and my code depends on spheroidical Earth calculation which I used(unlike above page used sphere Earth)
0

Author Commented:
Danaseaman and Ark,

Thanks for your efforts so far, however, there are still some problems somewhere and I'm confused.

Danaseaman, your code appears to work given some co-ordinates but is quite obvious it's failing on others given the result.

Ark, your code appears to be giving more accurate results until checking against the online calculator for which you supplied the link.

Example, using the co-ordinates mention above (Glasgow and Heathrow)

Danaseaman's code returns 136 degrees (after reversing the Dest/Cur co-ords)

Ark's code returns 152 degrees

The OnLine calculator returns 209 degrees.

I'll be more than happy to accept a couple of degree difference but the figures shown are a little to wide of the mark.

I will also be happy to split the points between you both should you both come back with a suitable answer

Thank you for your efforts so far

Dave
0

Commented:
One problem in testing is that longitudes you supplied for Glasgow and Heathrow should be negative.

'Glasgow
'55.870503  -4.445725

'Heathrow
'Lat 51.4833 Long -0.4500

Using online calculator with above values I get distance of 540.1 Km and bearing of 149°09&#8242;05&#8243;

Arks code with same values gives 151.62 deg. Rhumb Line or 150.03 deg. Great Circle
Distance looks close too, 299 nm(556 km)

I would go with his Arks solution.

0

Commented:
Thanks for points, glad I could help
Just checked for
'Glasgow
'55.870503  -4.445725

'Heathrow
'Lat 51.4833 Long -0.4500

Great Circle
=================
On-line calculator
A = 150°01&#8242;38&#8243;
D = 554.1 km

My program:
A = 150,027284473195 = 150°01&#8242;38.22&#8243;
D = 298,978134462858 nautical miles (ie. minutes of meridian) Equation 1 nm = 1852 m is average for meridional radius from equator to pole. For 6371 km radius (as given in on-line calculator page formulas) equation is 1 nm = 1853.24877740931 (PI * 6371000 / 180 / 60). In this case D = 554,080862165407 kM. If you need output in kilometers, you can use following function:

'Calculate distance via Great Circle in kilometers
Public Function DistanceGC_KM(ptFrom As POSITION, ptTo As POSITION) As Double
DistanceGC_KM = Rad(Abs(Arccosd(Sind(ptFrom.lat) * Sind(ptTo.lat) + Cosd(ptFrom.lat) * Cosd(ptTo.lat) * Cosd(ptTo.lon - ptFrom.lon)))) * 6371
End Function

'======================
Rhumb Line

On-line calculator:
A = 151°40&#8242;53&#8243;
D = 554.2 km

My program:
A = 151,624973888495 = 151°37&#8242;39.2&#8243;
D = 299,094954464381 = 554,297358690407 kM
Note, that for Rhumb Line Azimuth calculation I used Earth excentricity (more precessious), but for Distance I used approximate formula. More correct is to use already calculated azimuth:
D = dF/Cos(Az)

Private Function DistanceRL(ptFrom As POSITION, ptTo As POSITION) As Double
Dim dF As Double, dL As Double, Az As Double
dF = ptTo.lat - ptFrom.lat
dL = ptTo.lon - ptFrom.lon
If dL > 180 Then
dL = dL - 360
ElseIf dL < -180 Then
dL = dL + 360
End If
dR = MCH(ptTo.lat) - MCH(ptFrom.lat)
Az = Atan2(dL * 60, dR)
DistanceRL = dF / Cos(Az) * 60
End Function
0
Question has a verified solution.

Are you are experiencing a similar issue? Get a personalized answer when you ask a related question.

Have a better answer? Share it in a comment.