1. Add these functions to your Visual Basic project.

' Converts decimal degrees to radians.
' Divide radians by DDtoRAD to convert back to decimal degrees.
Public Const DDtoRAD = 3.14159265358979 / 180#
Public a As Double, e2 As Double

' Change the library path if necessary. The geoddist function provides
' access to the pe_geodesic_distance function in the Projection Engine
' library, pe90.dll. The function inputs are:
'
'      a = semimajor axis of a spheroid
'     e2 = eccentricity, e2 = f*(2-f) where f = flattening
'   lam1 = longitude of the first point in radians
'   phi1 = latitude of the first point in radians
'   lam2 = longitude of the second point in radians
'   phi2 = latitude of the second point in radians
'
' The function returns:
'
'   distance = geodesic distance between the two points in meters
'   az12 = geodetic azimuth from point 1 to point 2 in radians
'   az21 = geodetic azimuth from point 2 to point 1 in radians
'
' The returned distance is a geodesic distance. A great circle distance
' is on a sphere. The geodesic is on a spheroid/ellipsoid. If you wish
' the great circle distance, make sure that the geographic coordinate
' system uses a sphere. Ex: esriSRGeoCS_Authalicsphere (R = 6371000 m)


Private Declare Function geoddist Lib "c:\arcgis\arcexe90\bin\pe90.dll" _
  Alias "pe_geodesic_distance" (ByVal a As Double, ByVal e2 As Double, _
  ByVal lam1 As Double, ByVal phi1 As Double, _
  ByVal lam2 As Double, ByVal phi2 As Double, _
  distance As Double, azi12 As Double, azi21 As Double) As Long
  
' Change the library path if necessary. The geodcoord function provides
' access to the pe_geodesic_coordinate function in the Projection Engine
' library, pe82.dll. The function inputs are:
'
'      a = semimajor axis of a spheroid
'     e2 = eccentricity, e2 = f*(2-f) where f = flattening
'   lam1 = longitude of the first point in radians
'   phi1 = latitude of the first point in radians
'   distance = geodesic distance between the two points in meters
'   az12 = geodetic azimuth from point 1 to point 2 in radians
'
' The function returns:
'
'   lam2 = longitude of the second point in radians
'   phi2 = latitude of the second point in radians
'
' The returned coordinate is based on a geodesic distance. A great circle
' is on a sphere. The geodesic is on a spheroid/ellipsoid. If you wish
' to use a great circle distance, make sure that the geographic coordinate
' system uses a sphere. Ex: esriSRGeoCS_Authalicsphere (R = 6371000 m)

Private Declare Function geodcoord Lib "c:\arcgis\arcexe90\bin\pe90.dll" _
  Alias "pe_geodesic_coordinate" (ByVal a As Double, ByVal e2 As Double, _
  ByVal lam1 As Double, ByVal phi1 As Double, _
  ByVal distance As Double, ByVal azi12 As Double, _
  lam2 As Double, phi2 As Double) As Long

Public Sub main()

' By default, these functions calculate the geodesic distance or coordinate.
' A geodesic distance is calculated on a spheroid. A great circle distance
' is calculated on a sphere. If you want the great circle distance, set
' the a parameter to the radius of the sphere and set e2 to zero.

   geodesic_distance
   geodesic_coordinate
   
End Sub

Public Sub geodesic_distance(ByVal lon1 As Double, ByVal lat1 As Double, _
  ByVal lon2 As Double, ByVal lat2 As Double, _
  dist As Double, a12 As Double, a21 As Double)

  ' Initialization of variables needed by the geoddist function.
  Dim lam1 As Double, phi1 As Double
  Dim lam2 As Double, phi2 As Double
  Dim distance As Double
  Dim azi12 As Double, azi21 As Double

  ' Defining the geographic coordinate system.
  ' SetSpheroid retrieves the spheroid parameters. Change the
  ' geographic coordinate system in the SetSpheroid subroutine. 
  SetSpheroid a, e2
  
  ' Definition of two input points and conversion to radians. 
  lam1 = lon1 * DDtoRAD
  phi1 = lat1 * DDtoRAD
  lam2 = lon2 * DDtoRAD
  phi2 = lat2 * DDtoRAD

  ' The function is called. 
  geoddist a, e2, lam1, phi1, lam2, phi2, distance, azi12, azi21
  dist = distance
  a12 = azi12 / DDtoRAD
  a21 = azi21 / DDtoRAD
  
  ' Print or store the results. The azimuths are converted back to
  ' decimal degrees.
'  Debug.Print distance
'  Debug.Print azi12 / DDtoRAD
'  Debug.Print azi21 / DDtoRAD 

End Sub

Public Sub geodesic_coordinate(ByVal lon1c As Double, ByVal lat1c As Double, _
     ByVal distc As Double, ByVal azi12c As Double, _
     lon2c As Double, lat2c As Double)

  ' Initialization of variables needed by the geodcoord function. 
  Dim lam1 As Double, phi1 As Double
  Dim lam2 As Double, phi2 As Double
  Dim distance As Double
  Dim a12 As Double
  
  ' Defining the geographic coordinate system.
  ' SetSpheroid retrieves the spheroid parameters. Change the
  ' geographic coordinate system in the SetSpheroid subroutine. 
  SetSpheroid a, e2
  
  ' Convert the input point and azimuth to radians. 
  lam1 = lon1c * DDtoRAD
  phi1 = lat1c * DDtoRAD
  a12 = azi12c * DDtoRAD
  
  ' Assign the distance 
  distance = distc

  ' The function is called. 
  geodcoord a, e2, lam1, phi1, distance, a12, lon2c, lat2c
  
  ' Print or store the results. The second point is converted back to
  ' decimal degrees. 
  lon2c = lon2c / DDtoRAD
  lat2c = lat2c / DDtoRAD

End Sub

Public Sub SetSpheroid(a As Double, e2 As Double)
  Dim f As Double
  
  Dim pSpatialReferenceEnv As SpatialReferenceEnvironment
  Set pSpatialReferenceEnv = New SpatialReferenceEnvironment
  
  ' Defining the geographic coordinate system.
  ' Change esriSRGeoCSType enumeration to the appropriate geographic
  ' coordinate system. 
  Dim pGCS As IGeographicCoordinateSystem
  Dim pSpatialReference As ISpatialReference
  Set pGCS = pSpatialReferenceEnv.CreateGeographicCoordinateSystem(esriSRGeoCS_NAD1983)
  
  ' This section retrieves the spheroid parameters from the geographic
  ' coordinate system. 
  Dim pDatum As IDatum
  Dim pSpheroid As ISpheroid
  Set pDatum = pGCS.Datum
  Set pSpheroid = pDatum.Spheroid
  a = pSpheroid.SemiMajorAxis
  f = pSpheroid.Flattening
  e2 = f * (2# - f)

End Sub