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

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