Public Const Pi = 3.14159265358979

Public Function Radians(x)
    Radians = Pi * x / 180#
End Function

Private Function arccos(x)
    arccos = Atn(-x / Sqr(-x * x + 1)) + Pi / 2
End Function

Private Function arcsin(x)
    arcsin = Pi / 2# - arccos(x)
End Function

Public Function Bearing(Lat1, Lon1, Lat2, Lon2)
' Lat1, Lon1 - lat and lon for position 1
' Lat2, Lon2 - lat and lon for position 2
' Returns initial bearing (degrees) of great circle route from
' position 1 to position 2
' Assumes input is North = + , East = +
    If (Lon1 = Lon2) Then
        If (Lat1 < Lat2) Then
            Bearing = 0#
        Else
            Bearing = 180#
        End If
    Else
        dd = Posdist(Lat1, -Lon1, Lat2, -Lon2)
        arad = arccos((Sin(Radians(Lat2)) - Sin(Radians(Lat1)) * Cos(Radians(dd / 60))) / _
             (Sin(Radians(dd / 60)) * Cos(Radians(Lat1))))
        Bearing = arad * 180 / Pi
        If (Sin(Radians(Lon2 - Lon1)) < 0) Then
            Bearing = 360 - Bearing
        End If
    End If
End Function

       
Public Function Posdist(Lat1, Lon1, Lat2, Lon2)
' Lat1, Lon1 - lat and lon for position 1
' Lat2, Lon2 - lat and lon for position 2
' Returns distance between 2 positions in nautical miles
    If (Lat1 = Lat2 And Lon1 = Lon2) Then
       Posdist = 0
    Else
       Rlat1 = Radians(Lat1)
       Rlat2 = Radians(Lat2)
       Rlon = Radians(Lon2 - Lon1)
       Posdist = 60 * (180 / Pi) * arccos(Sin(Rlat1) * Sin(Rlat2) + _
              Cos(Rlat1) * Cos(Rlat2) * Cos(Rlon))
    End If
End Function


Public Function NewPosLat(Lat1, Lon1, Bearing, Distance)
' Computes latitude of a new position at a particular initial
' bearing and distance from an initial position (lat1, lon1)
' following a great circle route
' north = + east = +
' reference for formulas - std mathematical tables pg 176
' Latitude will not be correct if looped around pole in north-south direction
' nor across the equator - possibly this could be corrected
' Lat1, Lon1 - original position in decimal degrees
' Initial bearing to new point in true degrees
' Distance to new point in nautical miles
' Returns latitude in decimal degrees for new point
       If (Bearing = 0 Or Bearing = 180 Or Bearing = 360) Then
           If (Bearing = 180) Then
               NewPosLat = Lat1 - Distance / 60
           Else
               NewPosLat = Lat1 + Distance / 60
           End If
       Else
           a = Radians(90 - Lat1)
           Rlon1 = Radians(Lon1)
           b = Pi * Distance / (60 * 180)
           CAngle = Radians(Bearing)
           APlusB = 2# * Atn(Cos((a - b) / 2#) / (Cos((a + b) / 2#) * Tan(CAngle / 2#)))
           AMinusB = 2# * Atn(Sin((a - b) / 2#) / (Sin((a + b) / 2#) * Tan(CAngle / 2#)))
           BAngle = (APlusB - AMinusB) / 2#
           Aangle = (AMinusB + APlusB) / 2#
           If (Lat1 < 0) Then
               NewPosLat = (180# * arcsin(Sin(b) * Sin(CAngle) / Sin(BAngle)) / Pi) - 90#
           Else
               NewPosLat = Abs((180# * arcsin(Sin(b) * Sin(CAngle) / Sin(BAngle)) / Pi) - 90#)
           End If
       End If
       If (NewPosLat > 90) Then
           NewPosLat = NewPosLat - 90
       Else
           If (NewPosLat < -90) Then
               NewPosLat = -(NewPosLat + 180)
           End If
       End If
End Function


Public Function NewPosLon(Lat1, Lon1, Bearing, Distance)
' Computes longitude of a new position at a particular initial
' bearing and distance from an initial position (lat1, lon1)
' following a great circle route
' north = + east = +
' reference for formulas - std mathematical tables pg 176
' Longitude will not be correct if taken over either pole
' Lat1, Lon1 - original position in decimal degrees
' Initial bearing to new point in true degrees
' Distance to new point in nautical miles
' Returns longitude in decimal degrees for new point
       If (Bearing = 0) Bearing=360
       a = Radians(90 - Lat1)
       Rlon1 = Radians(Lon1)
       b = Pi * Distance / (60 * 180)
       CAngle = Radians(Bearing)
       APlusB = 2# * Atn(Cos((a - b) / 2#) / (Cos((a + b) / 2#) * Tan(CAngle / 2#)))
       AMinusB = 2# * Atn(Sin((a - b) / 2#) / (Sin((a + b) / 2#) * Tan(CAngle / 2#)))
       BAngle = (APlusB - AMinusB) / 2#
       Aangle = (AMinusB + APlusB) / 2#
       NewPosLon = Lon1 + BAngle * 180 / Pi
       If (NewPosLon > 180) Then
          NewPosLon = NewPosLon - 360
       Else
          If (NewPosLon < -180) Then
             NewPosLon = NewPosLon + 360
          End If
       End If
End Function
    

Public Function RetDist(Height, RadPerReticle, Reticles)
' Height in meters
' RaderReticle = Radians per reticle mark
' Reticles = number of reticles below horizon
' RetDist = distance in nautical miles
    x = Sqr(2 * 6366 * Height / 1000 + (Height / 1000) ^ 2)
    If (Reticles = 0) Then
        RetDist = x / 1.852
    Else
        Angle = Atn(x / 6366)
        RetDist = (1 / 1.852) * ((6366 + Height / 1000) * Sin(Angle + Reticles * RadPerReticle) - _
                   Sqr(6366 ^ 2 - ((6366 + Height / 1000) * Cos(Angle + Reticles * RadPerReticle)) ^ 2))
    End If
End Function

Public Function RetDist7x50(Height, Reticles)
' Height in meters
' RadPerReticle = Radians per reticle mark
' Reticles = number of reticles below horizon
' RetDist = distance in nautical miles
    RadPerReticle = 0.00497
    RetDist7x50 = RetDist(Height, RadPerReticle, Reticles)
End Function

Public Function RetDistBE(Height, Reticles)
' Height in meters
' RadPerReticle = Radians per reticle mark
' Reticles = number of reticles below horizon
' RetDist = distance in nautical miles
    RadPerReticle = 0.00136
    RetDistBE = RetDist(Height, RadPerReticle, Reticles)
End Function


Public Function MinSecToDeg(MinSec)
' Converts position field from  degrees.minutes seconds
' specified as deg.mmss (e.g., 36 degrees, 58 minutes and 34 seconds is
' 36.5834 )to decimal degrees
  xMinSec = Abs(MinSec)
  MinSecToDeg = Sign(MinSec) * (Int(xMinSec) + Int((xMinSec - Int(xMinSec)) * 100 + 0.00001) _
            / 60 + 100 * (xMinSec * 100 - Int(xMinSec * 100 + 0.00001)) / 3600)
End Function

Public Function DegToMinSec(DecDeg)
' Converts position field from decimal degrees to
' degrees.minutes seconds (deg.mmss)
  xDecDeg = Abs(DecDeg)
  Ideg = Int(xDecDeg)
  Decim = xDecDeg - Ideg
  imin = Int(Decim * 60)
  isec = Int((Decim * 60 - imin) * 60 + 0.5)
  DegToMinSec = Sign(DecDeg) * (Ideg + imin / 100 + isec / 10000)
End Function

Public Function DegToMinTen(DecDeg)
' Converts position field from decimal degrees to degrees.minutes tenths
' specified as deg.mmtt
  xDecDeg = Abs(DecDeg)
  Ideg = Int(xDecDeg)
  Decim = xDecDeg - Ideg
  imin = Decim * 0.6
  DegToMinTen = Sign(DecDeg) * (Ideg + imin)
End Function

Public Function MinTenToDeg(MinTen)
' Converts position field from  degrees.minutes hundreds of minutes
' specified as deg.mmtt (e.g., 36 degrees, 58.77 minutes is
' 36.5877 )to decimal degrees
  xMinTen = Abs(MinTen)
  MinTenToDeg = Sign(MinTen) * (Int(xMinTen) + Int((xMinTen - Int(xMinTen)) * 100 + 0.00001) / 60 + _
             100 * (xMinTen * 100 - Int(xMinTen * 100 + 0.00001)) / 6000)

End Function

Public Function ClinoDist(Altitude, Angle)
' Computes ground distance from altitude and angle from horizon
' Distance units are the same as altitude units
  ClinoDist = Altitude * Tan(Radians(90 - Angle))
End Function

Public Function ClinoArcDist(Altitude, Angle)
' Computes ground distance from altitude and angle from horizon
' Altitude should be specified in meters
' Output distance units are nautical miles
  EarthRadius = 3440
  Height = Altitude / 1852
  Alpha = ((2 * Height * EarthRadius + Height ^ 2) ^ 0.5) / EarthRadius
  If (Radians(Angle) < Alpha) Then
      ClinoArcDist = "Beyond horizon"
  Else
      ClinoArcDist = EarthRadius * (Radians(Angle) - arccos((EarthRadius + Height) * Cos(Radians(Angle)) / EarthRadius))
  End If
End Function

Public Function DistRet(Height, RadPerReticle, Distance)
' Height in meters
' RaderReticle = Radians per reticle mark
' Distance = distance in nautical miles
' DistRet = number of reticles from the horizon
   x1 = Distance / (6366 / 1.852)
   Angle = arccos(6366 / (6366 + Height / 1000))
   x2 = (x1 - Angle) ^ 2 / (2 * x1)
   DistRet = x2 / RadPerReticle
End Function

Public Function SASB(b, a, c)
' Triangle with sides a,b,c and angles A,B,C
'     Angle A is between b and c opposite a
'     Angle B is between c and a opposite b
'     Angle C is between a and b opposite c
' SASB - returns angle B
SideA = SASa(b, a, c)
SASB = 180 * arccos((c ^ 2 - b ^ 2 + SideA ^ 2) / (2 * c * SideA)) / Pi
End Function

Public Function SASC(b, a, c)
' Triangle with sides a,b,c and angles A,B,C
'     Angle A is between b and c opposite a
'     Angle B is between c and a opposite b
'     Angle C is between a and b opposite c
' SASC - returns angle C
SideA = SASa(b, a, c)
SASC = 180 * arccos((b ^ 2 - c ^ 2 + SideA ^ 2) / (2 * SideA * b)) / Pi
End Function


Public Function SASa(b, a, c)
' Triangle with sides a,b,c and angles A,B,C
'     Angle A is between b and c opposite a
'     Angle B is between c and a opposite b
'     Angle C is between a and b opposite c
' SASa - returns length of a
SASa = (b ^ 2 + c ^ 2 - 2 * b * c * Cos(Radians(a))) ^ 0.5
End Function

Public Function ASAb(b, a, c)
' Triangle with sides a,b,c and angles A,B,C
'     Angle A is between b and c opposite a
'     Angle B is between c and a opposite b
'     Angle C is between a and b opposite c
' ASAb - returns length of b
AngleA = 180 - b - c
ASAb = a * Sin(Radians(b)) / Sin(Radians(AngleA))
End Function

Public Function ASAc(b, a, c)
' Triangle with sides a,b,c and angles A,B,C
'     Angle A is between b and c opposite a
'     Angle B is between c and a opposite b
'     Angle C is between a and b opposite c
' ASAc - returns length of c
AngleA = 180 - b - c
ASAc = a * Sin(Radians(c)) / Sin(Radians(AngleA))
End Function

Public Function SSSA(a, b, c)
' Triangle with sides a,b,c and angles A,B,C
'     Angle A is between b and c opposite a
'     Angle B is between c and a opposite b
'     Angle C is between a and b opposite c
' SSSA - returns angle A
SSSA = 180 * arccos((b ^ 2 + c ^ 2 - a ^ 2) / (2 * b * c)) / Pi
End Function

Public Function SSSB(a, b, c)
' Triangle with sides a,b,c and angles A,B,C
'     Angle A is between b and c opposite a
'     Angle B is between c and a opposite b
'     Angle C is between a and b opposite c
' SSSB - returns angle B
SSSB = 180 * arccos((a ^ 2 + c ^ 2 - b ^ 2) / (2 * a * c)) / Pi
End Function

Public Function SSSC(a, b, c)
' Triangle with sides a,b,c and angles A,B,C
'     Angle A is between b and c opposite a
'     Angle B is between c and a opposite b
'     Angle C is between a and b opposite c
' SSSC - returns angle C
SSSC = 180 * arccos((a ^ 2 + b ^ 2 - c ^ 2) / (2 * a * b)) / Pi
End Function

Public Function SSAC(a, c, AngleA)
' Triangle with sides a,b,c and angles A,B,C
'     Angle A is between b and c opposite a
'     Angle B is between c and a opposite b
'     Angle C is between a and b opposite c
' SSAC - returns angle C
SSAC = 180 * arcsin(c * Sin(Radians(AngleA)) / a) / Pi
End Function

Public Function SSAB(a, c, AngleA)
' Triangle with sides a,b,c and angles A,B,C
'     Angle A is between b and c opposite a
'     Angle B is between c and a opposite b
'     Angle C is between a and b opposite c
' SSAB - returns angle B
SSAB = 180 - AngleA - SSAC(a, c, AngleA)
End Function

Public Function SSASb(a, c, AngleA)
' Triangle with sides a,b,c and angles A,B,C
'     Angle A is between b and c opposite a
'     Angle B is between c and a opposite b
'     Angle C is between a and b opposite c
' SSASB - returns side b
SSASb = a * Sin(Radians(SSAB(a, c, AngleA))) / Sin(Radians(AngleA))
End Function

Function ClosestDistance(x1, y1, x2, y2, xp, yp)
' Computes the closest distance between a point (xp,yp) and the line
' segment defined by the points (x1,y1) and (x2,y2).
  Slope = (y2 - y1) / (x2 - x1)
  Int1 = y1 - x1 * Slope
  Int2 = yp + xp / Slope
  x0 = (Int2 - Int1) * Slope / (Slope ^ 2 + 1)
  y0 = Int2 - x0 / Slope
  If (y1 <> y2) Then
     If ((y0 - y1) * Sign(y2 - y1) < 0) Then
         ClosestDistance = ((yp - y1) ^ 2 + (xp - x1) ^ 2) ^ 0.5
     ElseIf ((y0 - y2) * Sign(y2 - y1) > 0) Then
         ClosestDistance = ((yp - y2) ^ 2 + (xp - x2) ^ 2) ^ 0.5
     Else
         ClosestDistance = ((yp - y0) ^ 2 + (xp - x0) ^ 2) ^ 0.5
     End If
  Else
     If ((x0 - x1) * Sign(x2 - x1) < 0) Then
         ClosestDistance = ((yp - y1) ^ 2 + (xp - x1) ^ 2) ^ 0.5
     ElseIf ((x0 - x2) * Sign(x2 - x1) > 0) Then
         ClosestDistance = ((yp - y2) ^ 2 + (xp - x2) ^ 2) ^ 0.5
     Else
         ClosestDistance = ((yp - y0) ^ 2 + (xp - x0) ^ 2) ^ 0.5
     End If
  End If
End Function
  
Function Sign(x)
' Sign function - if x<0 = -1, else 1
   If (x < 0) Then
     Sign = -1
   Else
    Sign = 1
   End If
End Function