':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
':::                                                                         :::
':::  This routine calculates the distance between two points (given the     :::
':::  latitude/longitude of those points). It is being used to calculate     :::
':::  the distance between two locations using  GeoDataSource(TM) products   :::
':::                                                                         :::
':::  Definitions:                                                           :::
':::    South latitudes are negative, east longitudes are positive           :::
':::                                                                         :::
':::  Passed to function:                                                    :::
':::    lat1, lon1 = Latitude and Longitude of point 1 (in decimal degrees)  :::
':::    lat2, lon2 = Latitude and Longitude of point 2 (in decimal degrees)  :::
':::    unit = the unit you desire for results                               :::
':::           where: 'M' is statute miles (default)                         :::
':::                  'K' is kilometers                                      :::
':::                  'N' is nautical miles                                  :::
':::                                                                         :::
':::  Worldwide cities and other features databases with latitude longitude  :::
':::  are available at https://www.geodatasource.com	                     :::
':::                                                                         :::
':::  For enquiries, please contact sales@geodatasource.com                  :::
':::                                                                         :::
':::  Official Web site: https://www.geodatasource.com                       :::
':::                                                                         :::
':::             GeoDataSource.com (C) All Rights Reserved 2022              :::
':::                                                                         :::
':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

const pi = 3.14159265358979323846

Function distance(lat1, lon1, lat2, lon2, unit)
  Dim theta, dist
  If lat1 = lat2 And lon1 = lon2 Then
    distance = 0
  Else
    theta = lon1 - lon2
    dist = sin(deg2rad(lat1)) * sin(deg2rad(lat2)) + cos(deg2rad(lat1)) * cos(deg2rad(lat2)) * cos(deg2rad(theta))
    response.write "dist = " & dist & "<br>"
    dist = acos(dist)
    dist = rad2deg(dist)
    response.write "dist = " & dist & "<br>"
    distance = dist * 60 * 1.1515
    Select Case ucase(unit)
    Case "K"
      distance = distance * 1.609344
    Case "N"
      distance = distance * 0.8684
    End Select
  End If
End Function


'::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
':::  This function get the arccos function using arctan function   :::
'::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
Function acos(rad)
  If Abs(rad) <> 1 Then
    acos = pi/2 - Atn(rad / Sqr(1 - rad * rad))
  ElseIf rad = -1 Then
    acos = pi
  End If
End function


'::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
':::  This function converts decimal degrees to radians             :::
'::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
Function deg2rad(Deg)
  deg2rad = cdbl(Deg * pi / 180)
End Function

'::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
':::  This function converts radians to decimal degrees             :::
'::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
Function rad2deg(Rad)
  rad2deg = cdbl(Rad * 180 / pi)
End Function

msgbox(distance(32.9697, -96.80322, 29.46786, -98.53506, "M") & " Miles<br>")
msgbox(distance(32.9697, -96.80322, 29.46786, -98.53506, "K") & " Kilometers<br>")
msgbox(distance(32.9697, -96.80322, 29.46786, -98.53506, "N") & " Nautical Miles<br>")