The functions below do about the same thing. One just takes longer and is
just a bit more accurate.

I converted the long function from VB code by Jay Tanner who
based his code on the mathematical method of H. Andoyer.

The short function is all over the web and I'm not sure which one I used.

Oh, ya, if you've got questions, don't ask me! I'm not the mathmatician.

#COMPILE EXE
#DEBUG ERROR ON
#DIM ALL

%LongCompute = 1

MACRO Pi = 3.14159265358979##

#IF %LongCompute
FUNCTION Geo_Dist_Between (BYVAL Lat1 AS EXTENDED, _
                           BYVAL Lng1 AS EXTENDED, _
                           BYVAL Lat2 AS EXTENDED, _
                           BYVAL Lng2 AS EXTENDED, _
                       OPT BYVAL UoM   AS LONG      ) AS EXTENDED

  DIM C   AS LOCAL EXTENDED
  DIM CF  AS LOCAL EXTENDED
  DIM CL  AS LOCAL EXTENDED
  DIM CG  AS LOCAL EXTENDED
  DIM D   AS LOCAL EXTENDED
  DIM F   AS LOCAL EXTENDED
  DIM G   AS LOCAL EXTENDED
  DIM H1  AS LOCAL EXTENDED
  DIM H2  AS LOCAL EXTENDED
  DIM L   AS LOCAL EXTENDED
  DIM O   AS LOCAL EXTENDED
  DIM R   AS LOCAL EXTENDED
  DIM S   AS LOCAL EXTENDED
  DIM SF  AS LOCAL EXTENDED
  DIM SG  AS LOCAL EXTENDED
  DIM SL  AS LOCAL EXTENDED
  DIM W1  AS LOCAL EXTENDED
  DIM W2  AS LOCAL EXTENDED
  DIM T1  AS LOCAL EXTENDED
  DIM T2  AS LOCAL EXTENDED

  F = (Lat1 + Lat2) / 2                 ' Compute auxiliary angles
  G = (Lat1 - Lat2) / 2
  L = (Lng1 - Lng2) / 2

  T1 = Pi / 180##                       ' Compute sines and cosines of auxiliary angles
  SG  = SIN(G * T1) : CG = COS(G * T1)
  SF  = SIN(F * T1) : CF = COS(F * T1)
  SL  = SIN(L * T1) : CL = COS(L * T1)

  W1 = SG * CL : W1 = W1 * W1
  W2 = CF * SL : W2 = W2 * W2
  S  = W1 + W2

  W1 = CG * CL: W1 = W1 * W1
  W2 = SF * SL: W2 = W2 * W2
  C  = W1 + W2

  O = ATN(SQR(S / C))
  R = SQR(S * C) / O
  D = 2 * O

  SELECT CASE UoM
    CASE 1    : D = D * 3963.10##   ' Miles
    CASE 2    : D = D * 3443.90##   ' Nautical Miles
    CASE ELSE : D = D * 6378.14##   ' Kilometers
  END SELECT

  H1 = (3 * R - 1) / (2 * C)
  H2 = (3 * R + 1) / (2 * S)


  T1 = 0.0033528131778969144                ' Compute the angle between the points on a
  W1 = SF * CG : W1 = W1 * W1 * H1 * T1 + 1 ' synthetic sphere connecting the two points.

  W2 = CF * SG : W2 = W2 * W2 * H2 * T1

  FUNCTION = D * (W1 - W2)                  ' Return the distance between the given locations

END FUNCTION

#ELSE

FUNCTION Geo_Dist_Between (BYVAL Lat1 AS EXTENDED, _
                           BYVAL Lng1 AS EXTENDED, _
                           BYVAL Lat2 AS EXTENDED, _
                           BYVAL Lng2 AS EXTENDED, _
                       OPT BYVAL UoM   AS LONG      ) AS EXTENDED

  DIM Tmp AS LOCAL EXTENDED

  Lng1 = Lng1 / 57.29577951308232##
  Lat1 = Lat1 / 57.29577951308232##
  Lng2 = Lng2 / 57.29577951308232##
  Lat2 = Lat2 / 57.29577951308232##

  Tmp  = (SIN(Lat1) * SIN(Lat2))
  Tmp  = Tmp + (COS(Lat1) * COS(Lat2) * COS(Lng1 - Lng2))
  Tmp  = (Pi / 2) - ATN(Tmp / SQR(1-Tmp*Tmp))  'ArcCos

  SELECT CASE UoM
    CASE 1    : Tmp = Tmp * 3963.10##   ' Miles
    CASE 2    : Tmp = Tmp * 3443.90##   ' Nautical Miles
    CASE ELSE : Tmp = Tmp * 6378.14##   ' Kilometers
  END SELECT

  FUNCTION = Tmp

END FUNCTION
#ENDIF
'
'-----------------------------------------
'
FUNCTION PBMAIN ()

  DIM Cty1  AS LOCAL STRING
  DIM D     AS LOCAL LONG
  DIM Dist  AS LOCAL EXTENDED
  DIM Lat1  AS LOCAL EXTENDED
  DIM Lat2  AS LOCAL EXTENDED
  DIM Lng1  AS LOCAL EXTENDED
  DIM Lng2  AS LOCAL EXTENDED
  DIM Miles AS LOCAL STRING
  DIM Txt   AS LOCAL STRING
  DIM UoM   AS LOCAL LONG

  UoM  = 1   ' 0 = Km, 1 = Miles, 2 = Nautical Miles

  SELECT CASE UoM
    CASE 1    : Miles = " Miles "
    CASE 2    : Miles = " Naut. Miles "
    CASE ELSE : Miles = " Kilometers "
  END SELECT

  Lat1 = VAL(READ$(1))
  Lng1 = VAL(READ$(2))
  Cty1 = READ$(3)

  FOR D = 4 TO DATACOUNT STEP 3
    Lat2 = VAL(READ$(D  ))
    Lng2 = VAL(READ$(D+1))
    Dist = Geo_Dist_Between(Lat1,Lng1,Lat2,Lng2,UoM)
    Txt  = Txt & FORMAT$(ROUND(Dist,2),"00#.#0") & Miles & Cty1 & " -> " & READ$(D+2) & $CRLF
  NEXT

  Txt = RTRIM$(Txt)
#IF %DEF(%PB_CC32)
  FOR UoM = 1 TO PARSECOUNT(Txt,$CRLF)
    PRINT PARSE$(Txt,$CRLF,UoM)
  NEXT
  INPUT FLUSH
  WAITKEY$
#ELSE
  MSGBOX Txt
#ENDIF

  '    Latitude   Longitude    City                  Zip   Rand McNally distance
  DATA 39.929597, -91.374150,  "Quincy, IL"      : ' 62301
  DATA 38.631551, -90.193000,  "St. Louis, MO"   : ' 63101 Should be 138 miles
  DATA 39.800950, -89.649990,  "Springfield, IL" : ' 62701 Should be 118 miles
  DATA 41.886456, -87.623250,  "Chicago, IL"     : ' 60601 Should be 316 miles

END FUNCTION