Great Circle Distance

Can anyone suggest a more accurate set of functions to calculate great circle distance?

Here’s what I am using, but the result doesn’t seem to be very accurate.

input data: lat:|N |33 |45 |43.00 long:|W |84 |27 |29.00 |

getDistance( long As String, lat As String ) As Double

dim d1, m1, s1 As Integer
dim d2, m2, s2 As Integer
Dim decLat, decLong, dist As Double

dim a, b, c As String

a = mid( long, 3, 4 )
b = mid( long, 8, 2 )
c = mid( long, 12, 7 )

d1 = val( a )
m1 = val( b )
s1 = val( c )

a = mid( lat, 3, 4 )
b = mid( lat, 8, 2 )
c = mid( lat, 12, 7 )

d2 = val( a )
m2 = val( b )
s2 = val( c )

decLong = decimalDeg( d1, m1, s1 )
if Left( long, 1 ) = “W” Then decLong = decLong * -1
decLat = decimalDeg( d2, m2, s2 )
if Left( lat, 1 ) = “S” Then decLat = decLat * -1
dist = greatCircleDistance( decLong, decLat )
Return dist

========================

getGreatCircleDistance( long As Double, lat As Double ) As Double

Const earthRadius = 3958.75 // miles
Const pi = 3.1415926535897932384626433832795
Const deg2rad = 0.01745329252
Dim x1, y1, x2, y2, a, b, c, distance As Double

Const myLat = 34.04527778
Const myLong = -118.54694444

x1 = myLong * deg2rad
y1 = mylat * deg2rad
x2 = long * deg2rad
y2 = lat * deg2rad

a = sin(( y2-y1 ) / 2.0 )^2
b = sin(( x2-x1 ) / 2.0 )^2
c = sqrt( a + cos( y2 ) * cos( y1 ) * b )

distance = 2 * asin( c ) * earthRadius

Return distance

====================

If you assume Earth is spherical, then you first need to find the solid angle (as measured from the Earth’s centre) between the two points on its surface. So:

  1. Convert the lat-long for the two points to 3D cartesian coordinates, so you have V1.x, V1.y, V1.z and V2.x, V2.y, V2.z, where V1 and V2 indicate vectors or points and .x, .y and .z are their direction cosines,
  2. The Dot Product of two vectors gives the solid angle between them, i.e. V1 dot V2 = acos(angle),
  3. Then get the great circle distance from the angle you’ve just calculated and Earth’s radius.

If you want to take into account the slight flattening of the Earth (about 1 in 300), it gets more involved.

Regards - Richard.

This is the formula that I use for great circle distance.

Function GCdistance(LatitudeLocation1 as double, LongitudeLocation1 as double, LatitudeLocation2 as double, LongitudeLocation2 as double) as double
  'Latitude and Longitude are decimal values
  DistPart1=Sin(Radians(LatitudeLocation1))*Sin(Radians(LatitudeLocation2))+_
    Cos(Radians(LatitudeLocation1))*Cos(Radians(LatitudeLocation2))*_
    Cos(Abs(Radians(LongitudeLocation2-LongitudeLocation1)))

  Miles= Atan(Sqrt((1-DistPart1^2)/DistPart1^2)) * 3959.871
  km= Atan(Sqrt((1-DistPart1^2)/DistPart1^2))* 6372.795
  return miles 'or km
End Function

Function Radians(d as double) as double
  const pi = 3.1415926535897932384626433832795
  return pi/180*d
End Function 

Thanks.

Robert:

Would you be kind enough to also show a function for converting from h.m.s to decimal?

I am using

decDegree( d As integer, m As integer, s As Integer ) as double

Dim tSeconds As Integer
Dim fractionalPart As Double
Dim result As Double

tSeconds = ( m * 60 ) + s

fractionalPart = tSeconds / 3600
result = d + fractionalPart

Return result

Since I don’t know the exact string format of your input, I can’t help with that without more info. This assumes that you’ve already extracted the degrees minutes and seconds from the input string, and converted them to numeric values:

Function DegMinSec2Decimal(d As Double, m As Double, s As Double, dir As String) as Double 'Convert degrees (or hours)/minutes/seconds to decimal degrees (or hours) 'West or South values are treated as negative. 'Parameters d,m,s are degrees, minutes and seconds respectively. 'Parameter dir should contain one of N,E,S,W. Any other value will give a zero result static dSign() As Integer = Array(0,1,1,-1,-1) return (d+m/60+s/3600)*dSign(instr("NESW",dir)) End Function

Thank you all!

Here’s the final code:

======
Function radians(d as double) As Double
const pi = 3.1415926535897932384626433832795
return pi/180*d
End Function

======
Function decimalDeg(dir As string, d As integer, m As integer, s As Double) As Double
Dim tSeconds As Integer
Dim fractionalPart As Double
Dim result As Double

tSeconds = ( m * 60 ) + s
fractionalPart = tSeconds / 3600 result = d + fractionalPart
if dir = “S” Then result = result * -1
if dir = “W” Then result = result * -1
Return result
End Function

======
Function GCdistance(Lat as double, Long as double) As double
Const myLat = 34.04527778
Const myLong = -118.54694444
Dim DistPart1, Miles, km As Double

DistPart1=Sin(Radians(myLat))*Sin(Radians(Lat))+_ Cos(Radians(myLat))Cos(Radians(Lat))_ Cos(Abs(Radians(Long-myLong)))

Miles= Atan(Sqrt((1-DistPart1^2)/DistPart1^2)) * 3959.871
km= Atan(Sqrt((1-DistPart1^2)/DistPart1^2))* 6372.795
Return miles 'or km
End Function