Attribute VB_Name = "mod_DME"
Option Explicit
Private Const Pi As Double = 3.14159265358979

Public Sub logthis(str As String)

    Dim fp As Integer
    
    fp = FreeFile()

    Open "goat.log" For Append As #fp
    Print #fp, str
    Close #fp
    
End Sub


'Based on vb code example from
'http://msdn.microsoft.com/en-us/library/6x627e5f%28v=vs.90%29.aspx
'returns the tokenth string delimited by one or delimiters
Public Function getToken(ByVal s1 As String, ByVal token As Integer, ByVal delimiter As String) As String
    On Error GoTo gettokenError
    Dim i As Long
    Dim sa() As String 'string array
    Dim lastnonempty As Integer
    Dim result As String

    sa = Split(s1, delimiter)
    lastnonempty = -1
    
    For i = 0 To UBound(sa)
        If (sa(i) <> "") Then
            lastnonempty = lastnonempty + 1
            sa(lastnonempty) = sa(i)
        End If
    Next

    ReDim Preserve sa(lastnonempty + 1)
    
    result = sa(token - 1)
    
    getToken = result
    Exit Function
    
gettokenError:
    getToken = ""
End Function




'Intermediate Point on Great Circle
Public Function IPGC(Latitude1 As Double, Longitude1 As Double, Latitude2 As Double, Longitude2 As Double, nm As Double, Fraction As Double) As Double()

    Dim D As Double
    Dim lat As Double
    Dim lon As Double
    Dim x As Double
    Dim y As Double
    Dim Z As Double
    Dim a As Double
    Dim b As Double
    Dim lat1 As Double
    Dim lon1 As Double
    Dim lat2 As Double
    Dim lon2 As Double

    
    D = (nm * Pi) / (180 * 60)
    lat1 = Latitude1 * (Pi / 180)
    lon1 = Longitude1 * (Pi / 180)
    lat2 = Latitude2 * (Pi / 180)
    lon2 = Longitude2 * (Pi / 180)

    a = Sin((1 - Fraction) * D) / Sin(D)
    b = Sin(Fraction * D) / Sin(D)
    x = a * Cos(lat1) * Cos(lon1) + b * Cos(lat2) * Cos(lon2)
    y = a * Cos(lat1) * Sin(lon1) + b * Cos(lat2) * Sin(lon2)
    Z = a * Sin(lat1) + b * Sin(lat2)
    
    lat = Atan2(Z, SqRT(x * x + y * y)) * (180 / Pi)
    lon = Atan2(y, x) * (180 / Pi)

    Dim rar(2) As Double
    
    rar(0) = lat
    rar(1) = lon
    
    IPGC = rar
    
End Function







'The LAT4RD and LON4RD functions NO LONGER expect +N, +W, and -S, -E
Public Function Lat4RD(Latitude As Double, Longitude As Double, nm As Double, radial As Double) As Double

    Dim D As Double
    Dim lat As Double
    Dim x As Double
    Dim y As Double
    Dim Z As Double
    Dim lon As Double
    Dim lat1 As Double
    Dim lon1 As Double
    Dim tc As Double
    
    
    
    D = (nm * Pi) / (180 * 60)
    lat1 = Latitude * (Pi / 180)
    lon1 = Longitude * (Pi / 180)
    

    tc = radial * (Pi / 180)

    lat = ASin(Sin(lat1) * Cos(D) + Cos(lat1) * Sin(D) * Cos(tc))
    x = lon1 - ASin(Sin(tc) * Sin(D) / Cos(lat)) + Pi
    y = 2 * Pi
    Z = fmod(x, y)
    lon = Z - Pi
    
    Lat4RD = lat * (180 / Pi)



End Function


'The LAT4RD and LON4RD functions NO LONGER expect +N, +W, and -S, -E
Public Function Lon4RD(Latitude As Double, Longitude As Double, nm As Double, radial As Double) As Double

    Dim D As Double
    Dim lat As Double
    Dim x As Double
    Dim y As Double
    Dim Z As Double
    Dim lon As Double
    Dim lat1 As Double
    Dim lon1 As Double
    Dim tc As Double
    
    D = (nm * Pi) / (180 * 60)
    lat1 = Latitude * (Pi / 180)
    lon1 = -Longitude * (Pi / 180)

    
    tc = radial * (Pi / 180)

    lat = ASin(Sin(lat1) * Cos(D) + Cos(lat1) * Sin(D) * Cos(tc))
    x = lon1 - ASin(Sin(tc) * Sin(D) / Cos(lat)) + Pi
    y = 2 * Pi
    Z = fmod(x, y)
    lon = Z - Pi
    Lon4RD = -lon * (180 / Pi)

End Function


Function calculate_distance(lat_1 As Double, lon_1 As Double, lat_2 As Double, lon_2 As Double) As Double

    Dim Cosd As Double
    Dim D As Double
    
    Dim lat1 As Double
    Dim long1 As Double
    Dim lat2 As Double
    Dim long2 As Double
    
    On Local Error Resume Next
    
    lat1 = lat_1
    long1 = lon_1
    lat2 = lat_2
    long2 = lon_2
    
    'Convert to radians
    lat1 = Radians(lat1)
    long1 = Radians(long1)
    lat2 = Radians(lat2)
    long2 = Radians(long2)
    
    Cosd = Sin(lat1) * Sin(lat2) + Cos(lat1) * Cos(lat2) * Cos(long1 - long2)
    
    D = ACos(Cosd) * ((180 * 60) / Pi)
    

    'Return the distance in nautical miles
    calculate_distance = D
    
End Function

Function GreatCircleRoute(lat_1 As Double, lon_1 As Double, lat_2 As Double, lon_2 As Double) As Double

    Dim lat1, lon1, lat2, lon2, nm, D, tc As Double
    

    

    nm = calculate_distance(lat_1, lon_1, lat_2, lon_2)
    
    If (nm = 0) Then
        GreatCircleRoute = -1
        Exit Function
    End If
        

    D = (3.1415617 / (180 * 60)) * nm


    lat1 = lat_1 * 3.1415617 / 180
    lon1 = -lon_1 * 3.1415617 / 180
    lat2 = lat_2 * 3.1415617 / 180
    lon2 = -lon_2 * 3.1415617 / 180


    If (Sin(lon2 - lon1) < 0) Then
        tc = ACos((Sin(lat2) - Sin(lat1) * Cos(D)) / (Sin(D) * Cos(lat1)))
    Else
        tc = 2 * 3.1415617 - ACos((Sin(lat2) - Sin(lat1) * Cos(D)) / (Sin(D) * Cos(lat1)))
    End If

    tc = tc * (180 / 3.1415617)
    
    GreatCircleRoute = tc



End Function

'Function ACos(r As Double) As Double
'    If r = -1 Then
'        ACos = 4 * Atn(1)
'    Else
'        ACos = Atn(-r / Sqr(-r * r + 1)) + 2 * Atn(1)
'    End If
'End Function



Function Radians(x As Double) As Double
    Radians = x * (Pi / 180)
End Function

'This is an improvment over VB6's Atn function, which returns the value in
'radians. This one returns the value in Decimal:

Function ArcTan(Value As Double) As Double

ArcTan = Atn(Value) * 180 / Pi

End Function

'
'ASin, ACos, ACot, ASec, ACsc - Missing inverse trig functions
' arc sine
' error if value is outside the range [-1,1]

Function ASin(Value As Double) As Double
    If Abs(Value) <> 1 Then
        ASin = Atn(Value / Sqr(1 - Value * Value))
    Else
        ASin = 1.5707963267949 * Sgn(Value)
    End If
End Function

' arc cosine
' error if NUMBER is outside the range [-1,1]

Function ACos(ByVal Number As Double) As Double
    If (Number < -1) Or (Number > 1) Then
        ACos = 3.1415
        Exit Function
    End If
    If Abs(Number) <> 1 Then
        ACos = 1.5707963267949 - Atn(Number / Sqr(1 - Number * Number))
    ElseIf Number = -1 Then
        ACos = 3.14159265358979
    ElseIf Number = 1 Then '--> Acos=0 (implicit)
    End If
    
End Function

' arc cotangent
' error if NUMBER is zero

Function ACot(Value As Double) As Double
    ACot = Atn(1 / Value)
End Function

' arc secant
' error if value is inside the range [-1,1]

Function ASec(Value As Double) As Double
    ' NOTE: the following lines can be replaced by a single call
    '            ASec = ACos(1 / value)
    If Abs(Value) <> 1 Then
        ASec = 1.5707963267949 - Atn((1 / Value) / Sqr(1 - 1 / (Value * Value)))
    Else
        ASec = 3.14159265358979 * Sgn(Value)
    End If
End Function

' arc cosecant
' error if value is inside the range [-1,1]

Function ACsc(Value As Double) As Double
    ' NOTE: the following lines can be replaced by a single call
    '            ACsc = ASin(1 / value)
    If Abs(Value) <> 1 Then
        ACsc = Atn((1 / Value) / Sqr(1 - 1 / (Value * Value)))
    Else
        ACsc = 1.5707963267949 * Sgn(Value)
    End If
End Function

Public Function fmod(x As Double, y As Double) As Double

    Dim i As Long
    Dim Z As Double


    i = Int(x / y)

    Z = x - (y * i)


    fmod = Z


End Function

Public Function SqRT(x As Double) As Double

    SqRT = Exp(Log(x) / 2)

End Function



Public Function PDP(altitude As Double, target As Double, distance As Double, ias As Double) As Double

    Dim T As Integer
    Dim flag As Integer
    Dim alt As Double
    Dim drop As Double
    Dim miles2go As Double
    Dim tas As Double
    Dim dist As Double
    Dim despro As Double
    Dim vv As Double
    Dim maxvv As Double
    Dim minvv As Double
    Dim ittg As Double
    Dim lowest_dist As Double
    Dim lowest_drop As Double
    Dim optimal As Double
    Dim Index As Double
    Dim x As Double
    Dim bm As Double
    Dim t1 As Double
    Dim t2 As Double
    Dim i As Integer
    Dim j As Integer
    Dim k As Integer
    Dim s As String
    
    'initialization
    tas = (ias + ((ias / 50) * (altitude / 1000)))
    
    minvv = (altitude - target) / (60 * (distance / ias))
    maxvv = (altitude - target) / (60 * (distance / tas))

    'estimate initial time to go, in minutes
    ittg = 60 * (distance / ias)
    
    'iterate through each possible descent rate (by 1)
    For vv = minvv To maxvv
        alt = altitude
        dist = distance
        'Iterate through time (by 1 minute)
        For T = 1 To ittg
       
            tas = ias + ((ias / 50) * (alt / 1000))
            If (alt < 10000) Then tas = 250 + ((250 / 50) * (alt / 1000))  'adjusting for 10/250 speed restriction
            alt = alt - vv  'feet per minute
            dist = dist - (tas / 60)
            drop = alt - target

            If ((drop < 0) Or (dist < 0)) Then
                Exit For
            End If
            
            If (flag = 0) Then
                Index = dist + drop
                lowest_dist = dist
                lowest_drop = drop
                optimal = vv
                flag = 1
            Else
                x = dist + drop
                If (x < Index) Then
                    optimal = vv
                    Index = x
                    lowest_dist = dist
                    lowest_drop = drop
                End If
            End If
        Next T
    Next vv

   
    PDP = optimal


End Function





 
Public Function Atan2(ByVal y As Double, ByVal x As Double) As Double
 
    If y > 0 Then
      If x >= y Then
        Atan2 = Atn(y / x)
      ElseIf x <= -y Then
        Atan2 = Atn(y / x) + Pi
      Else
        Atan2 = Pi / 2 - Atn(x / y)
      End If
    Else
      If x >= -y Then
        Atan2 = Atn(y / x)
      ElseIf x <= y Then
        Atan2 = Atn(y / x) - Pi
      Else
        Atan2 = -Atn(x / y) - Pi / 2
      End If
    End If
 
End Function


