Distance & Course Calculations

IP076

Somewhat of a Newbie
Local time
Today, 13:02
Joined
Sep 12, 2004
Messages
24
I've been using the following code to determine distance between two lat long cooridinates.

Code:
Function CalculateDistance(lat1, lon1, lat2, lon2) As Double
        Dim A As Double
        Dim B As Double
        Dim C As Double
        
        Dim X As Double
        
        A = Cos(lat1 * PI / 180) * Cos(lat2 * PI / 180) * Cos(lon1 * PI / 180) * Cos(lon2 * PI / 180)
        B = Cos(lat1 * PI / 180) * Sin(lon1 * PI / 180) * Cos(lat2 * PI / 180) * Sin(lon2 * PI / 180)
        C = Sin(lat1 * PI / 180) * Sin(lat2 * PI / 180)
        If (A + B + C) >= 1 Or (A + B + C) <= -1 Then
            CalculateDistance = 0
        Else
            X = (A + B + C)
            
            CalculateDistance = (Atn(-X / Sqr(-X * X + 1)) + 2 * Atn(1)) * RadiusEarth
        End If
End Function

Wish I could take credit for it...but I can't...although, I do not remember what site I got it from. This is by far the best one I've found, as it takes into account great circle routings for determining the distance.

I have also defined the following Constants:
Code:
Public Const PI = 3.14159265358979
Public Const RadiusEarth = 3437.74677 ' Nautical Miles
Public Const CircumEarth = (RadiusEarth * 2 * PI)
Public Const NMPerDegreeLat = CircumEarth / 360

This will return your distance value in Nautical Miles.

I'm wondering if anyone has any advice on how to modify this function to also obtain an initial course. I have done it, but the course I come up with, doesn't apply the great circle method... I think if I new what side of the triangle the A, B, and C above represented, I could figure it out....however, my small amount of trig knowledge has escaped my mind as of late, just hoping maybe someone can help out a little bit.

Thanks!
 
I am working on a solution!!! - I know this one for sure -

~Begin Confident Voice~ I am a surveyor ~End Confident Voice"

How do you want it returned? Azimuth in decimal degree's?, Bearing? N42E
 
Try this

Private Function CalcAz(lat1 As Double, lon1 As Double, lat2 As Double, lon2 As Double, d As Double) As Double
Dim X As Double
Dim phi As Double

X = Cos((Sin(lat2) - Sin(lat1) * Cos(d)) / (Sin(d) * Cos(lat1)))

If Sin(lon2 - lon1) < 0 Then phi = X

If Sin(lon2 - lon1) > 0 Then phi = 2 * PI - X

CalcAz = phi

End Function
 
IP076 said:
I'm wondering if anyone has any advice on how to modify this function to also obtain an initial course.

Bit drunk now but I'll have a look at Bowditch tomorrow to refresh myself with nautical mathematics. :)
 
Thanks for the info....I'll have to try out a couple of the above.

As for return...just a simple Decimal degrees would be perfect, such as 081.2 or 236.5.

Hey Surjer, whats the "d" in your function?
 
Last edited:
the d is the distance they are apart..

After a bit of thought I need to take the aCos instead of Cos
 
Surjer said:
the d is the distance they are apart..

After a bit of thought I need to take the aCos instead of Cos


Ahh..yes...Distance....

aCos is a standard VB function? I get a sub or function not defined compile error if I change all of your Cos to aCos.

Here's what I've done:

Code:
Public Function CalcAz(lat1 As Double, lon1 As Double, lat2 As Double, lon2 As Double) As Double
Dim X As Double
Dim phi As Double
Dim d As Double

d = CalculateDistance(lat1, lon1, lat2, lon2)

X = aCos((Sin(lat2) - Sin(lat1) * aCos(d)) / (Sin(d) * aCos(lat1)))

If Sin(lon2 - lon1) < 0 Then phi = X

If Sin(lon2 - lon1) > 0 Then phi = 2 * PI - X

CalcAz = phi
End Function

I changed the location of the d....not sure if it makes any difference, I wasnt sure if you could reference another function in the set for this function.

About the return value, I'm looking for a course (Degrees), as used for Aviation or Maritime purposes.
 
Last edited:
well it should be this but there is no freaking aCos in VBA.. Hmm what to do what to do... Maybe derive it from aTan but its way to late to be thinking this hard!

I'll check back monday and see then...

Code:
Public Function CalcAz(lat1 As Double, lon1 As Double, lat2 As Double, lon2 As Double) As Double
Dim X As Double
Dim phi As Double
Dim d As Double

d = CalculateDistance(lat1, lon1, lat2, lon2)

X = aCos((Sin(lat2) - Sin(lat1) * Cos(d)) / (Sin(d) * Cos(lat1)))

If Sin(lon2 - lon1) < 0 Then phi = X

If Sin(lon2 - lon1) > 0 Then phi = 2 * PI - X

CalcAz = phi
End Function
 
Last edited:
Code:
Public Function aCOS(ByVal nValue As Double, Optional fRadians As Boolean = True) As Double
        
        aCOS = -Atn(nValue / Sqr(1 - nValue * nValue)) + PI / 2
        If fRadians = False Then aCOS = aCOS * (PI / 180)
End Function

I found that code....but, I guess I'm having some trig problems with your supplied code...brain fart really...that whole value after the aCos (in your code) is going to have to be in Radians, correct? I believe that is true for the Atn function to work properly.

Really starting to think I should have paid more attention in trig!!
:D

Why couldn't bill just include Acos? Arghhhh....
 
Last edited:
This is in radians... The second part of your code you just posted will change it to degrees... VB does all its math in radians...
 
Code:
Public Function CalcAz(lat1 As Double, lon1 As Double, lat2 As Double, lon2 As Double) As Double
'This Function takes into account great circle routings

Dim x As Double
Dim phi As Double
Dim d As Double
Dim Lat1Radians As Double
Dim Lon1Radians As Double
Dim Lat2Radians As Double
Dim Lon2Radians As Double
Dim DistanceRadians As Double
Dim InitialCourseRadians As Double
Dim Course As Double

Lat1Radians = DegreesRadians(lat1)
Lat2Radians = DegreesRadians(lat2)
Lon1Radians = DegreesRadians(lon1)
Lon2Radians = DegreesRadians(lon2)

d = CalculateDistance(lat1, lon1, lat2, lon2) / NMPerDegreeLat
DistanceRadians = DegreesRadians(d)

x = (Sin(Lat2Radians) - Sin(Lat1Radians) * Cos(DistanceRadians)) / (Cos(Lat1Radians) * Sin(DistanceRadians))
    
InitialCourseRadians = aCOS(x, True)
phi = RadiansDegrees(InitialCourseRadians)

If lon2 > lon1 Then
    Course = 0 + Abs(phi)
Else
    Course = 360 - Abs(phi)
End If

CalcAz = Course
End Function

Got the above to work in all except one situation....It doesnt work if you cross the International Date Line....lets say sailing from Hawaii to Japan, it says I shoudl start out going North East. I'm fiddling with it trying to find a solution.
 
Okay, Here is my whole module
Code:
Option Compare Database
Option Explicit

Public Const PI = 3.14159265358979
Public Const RadiusEarth = 3437.74677 ' Nautical Miles
Public Const CircumEarth = (RadiusEarth * 2 * PI)
Public Const NMPerDegreeLat = CircumEarth / 360

Function CalculateDistance(lat1, lon1, lat2, lon2) As Double
        Dim A As Double
        Dim B As Double
        Dim C As Double
        
        Dim x As Double
        
        A = Cos(lat1 * PI / 180) * Cos(lat2 * PI / 180) * Cos(lon1 * PI / 180) * Cos(lon2 * PI / 180)
        B = Cos(lat1 * PI / 180) * Sin(lon1 * PI / 180) * Cos(lat2 * PI / 180) * Sin(lon2 * PI / 180)
        C = Sin(lat1 * PI / 180) * Sin(lat2 * PI / 180)
        If (A + B + C) >= 1 Or (A + B + C) <= -1 Then
            CalculateDistance = 0
        Else
            x = (A + B + C)
            
            CalculateDistance = (Atn(-x / Sqr(-x * x + 1)) + 2 * Atn(1)) * RadiusEarth
        End If
End Function

Public Function ConvertToDecimal(hdg As Boolean, deg As Double, min As Double, sec As Double)

    Dim degdec As Double
    Dim mindec As Double
    Dim secdec As Double
    
    secdec = sec / 60
    mindec = (min + secdec) / 60
    degdec = deg + mindec
    
    If hdg = True Then
        ConvertToDecimal = 0 - degdec
    Else
        ConvertToDecimal = degdec
    End If
    
End Function

Public Function NMPerDegreeLong(MiddleLat As Double)
    NMPerDegreeLong = (Cos(MiddleLat * (PI / 180)) * CircumEarth) / 360
End Function

Public Function DegreesRadians(Degrees As Double)

        DegreesRadians = Degrees * (PI / 180)
End Function

Public Function RadiansDegrees(Radians As Double)
        RadiansDegrees = Radians * (180 / PI)
End Function
Public Function aCOS(ByVal nValue As Double, Optional fRadians As Boolean = True) As Double
        
        aCOS = -Atn(nValue / Sqr(1 - nValue * nValue)) + PI / 2
        If fRadians = False Then aCOS = aCOS * (PI / 180)
End Function


Public Function WhatIsCourse(lat1 As Double, lon1 As Double, lat2 As Double, lon2 As Double) As String
'This Function takes into account great circle routings

Dim x As Double
Dim phi As Double
Dim d As Double
Dim Lat1Radians As Double
Dim Lon1Radians As Double
Dim Lat2Radians As Double
Dim Lon2Radians As Double
Dim DistanceRadians As Double
Dim InitialCourseRadians As Double
Dim Course As String
Dim DeltaLongitude As Double
Dim DeltaLatitude As Double

Lat1Radians = DegreesRadians(lat1)
Lat2Radians = DegreesRadians(lat2)
Lon1Radians = DegreesRadians(lon1)
Lon2Radians = DegreesRadians(lon2)

'-----------Determine Change (Delta) in Longitude------------

    DeltaLongitude = DeltaLong(lon1, lon2)

'------------------------------------------------------------
'-----------Determine Change (Delta) in Latitude-------------
    DeltaLatitude = DeltaLat(lat1, lat2)
'------------------------------------------------------------
'----------Get Distance from CalculateDistance Function------
d = CalculateDistance(lat1, lon1, lat2, lon2) / NMPerDegreeLat
'And Convert it to radians
DistanceRadians = DegreesRadians(d)

'Determine appropriate angle in radians
x = (Sin(Lat2Radians) - Sin(Lat1Radians) * Cos(DistanceRadians)) / (Cos(Lat1Radians) * Sin(DistanceRadians))
InitialCourseRadians = aCOS(x, True)
'convert angle to degrees
phi = RadiansDegrees(InitialCourseRadians)

'ensure angle matches direction of travel
If lon2 > lon1 Then
        Course = 0 + Abs(phi)
Else
        Course = 360 - Abs(phi)
End If

'allow for travel across date line
If DeltaLongitude > 180 Then
    If lon2 > lon1 Then
        Course = 360 - Abs(phi)
    Else
        Course = 0 + Abs(phi)
    End If
End If

'----Properly Format Output for 3 Digits before Decimal------
If Course >= 0 And Course < 10 Then
        Course = "00" & Format(Course, "standard") & " degrees"
    ElseIf Course >= 10 And Course < 100 Then
        Course = "0" & Format(Course, "standard") & " degrees"
    Else
        Course = Format(Course, "standard") & " degrees"
    End If
    
WhatIsCourse = Course
End Function

Public Function DeltaLong(lon1 As Double, lon2 As Double) As Double
If lon1 < 0 And lon2 < 0 Then
    DeltaLong = Abs(Abs(lon1) - Abs(lon2))
Else
    DeltaLong = Abs(lon1) + Abs(lon2)
End If
End Function

Public Function DeltaLat(lat1 As Double, lat2 As Double) As Double

    DeltaLat = (90 - lat1) + (90 - lat2)
    
End Function

Works Awesome....There are a few extra items in there, that I may or may not need in the future. THere is only one fault. If I try to move along a line of latitude it gives me an error with the aCos function. If I try to move along a line of longitude, it always returns 0. I was working on some IF and ELSE statements but was having a problem figuring them out. Thought maybe you could figure it out, well, I've had enough of this for a while....time for beer. Thanks for the help!
 

Users who are viewing this thread

Back
Top Bottom