Attribute VB_Name = "Arrow"

Public Et As Double
Public Mf As Double
Public Kf As Double
Public ArrowStraightness As Double 'mm
Public NockTravel As Double 'mm
Public RestMisalignment As Double 'mm
Public NockWeight As Double
Public InsertWeight As Double
Public FlightWeight As Double

Public TubeWeight As Double 'gr
Public TubeLength As Double 'inch assume to start at front of arrow

'    Function Arrow_Fc(StaticSpine As Double, ComponentWeights As Range, ComponentPositions As ) As Double
'        'Fc =   k/(Ys.L)  x  Wt / sum(Widi)
'        Dim k As Double
'        Dim W As Double
'        Dim Wt As Double
'        Dim L As Double
'        k = 5000000
'        W = 0
'        Wt = 0
'        L = 0
'        Dim i As Integer
'        For i = 1 To ComponentWeights.Count
'            Wt = Wt + ComponentWeights(i)
'            W = W + ComponentWeights(i) * ComponentPositions(i)
'            If ComponentPositions(i) > L Then L = ComponentPositions(i)
'        Next
'        Arrow_Fc = k / (StaticSpine * L) * Wt / W
'    End Function

'    Function Arrow_FcPerc(Fc As Double, DrawWeight As Double) As Double
'        Arrow_FcPerc = Fc / DrawWeight * 100
'    End Function
'
'    Function Arrow_FOC(ComponentWeights As Range, ComponentPositions As Range) As Double
'    '  FOC  = (  sum(Widi)  / (Wt.L)  - 0.5 ) x 100
'    Dim Widi As Double
'    Dim Wt As Double
'    Dim L As Double
'        Widi = 0
'        Wt = 0
'        L = 0
'        Dim i As Integer
'        For i = 1 To ComponentWeights.Count
'            Wt = Wt + ComponentWeights(i)
'            Widi = Widi + ComponentWeights(i) * ComponentPositions(i)
'            If ComponentPositions(i) > L Then L = ComponentPositions(i)
'        Next
'        Arrow_FOC = (Widi / (Wt * L) - 0.5) * 100
'    End Function
'
'    Function Arrow_FOC2(SumWiDi As Double, Wt As Double, L As Double) As Double
'    '  FOC  = (  sum(Widi)  / (Wt.L)  - 0.5 ) x 100
'        Arrow_FOC2 = (SumWiDi / (Wt * L) - 0.5) * 100
'    End Function
'
'    'Trad bows
'    Function Arrow_Frequency(PileWeight As Double, ShaftLength As Double, Spine As Double, ShaftGPI As Double) As Double
'        'Fn = 6135788 ( 0.5596 + 3/Wp)/L^2/(S.Ws)^0.5
'        'Wp is the pile weight of the arrow, L the shaft length, S the s-spine, and Ws the shaft's weight per inch (as GPI).
'        Arrow_Frequency = 6135788 * (0.5596 + 3 / PileWeight) / ShaftLength ^ 2 / (Spine * ShaftGPI) ^ 0.5
'    End Function
'
'    Private Function Y2M(yards As Double) As Double
'        Y2M = yards * 0.9144
'    End Function
'
'    Private Function I2M(inches As Double) As Double
'        I2M = inches * 2.54 / 100
'    End Function
'
'    Private Function Gr2Kg(grain As Double) As Double
'        Gr2Kg = grain * 0.065 / 1000
'    End Function
'
'    'Trad bows
'    Function Arrow_FrequencyForBow(DrawLength As Double, BraceHeight As Double, DrawForce As Double, ArrowWeight As Double) As Double
'        a = DrawForce / 2.2 * 10 / Gr2Kg(ArrowWeight) * 0.58
'        v = (2 * a * I2M(DrawLength - BraceHeight)) ^ 0.5
'        t = (2 * I2M(DrawLength - BraceHeight) / a) ^ 0.5 + I2M(BraceHeight) / 2 / v
'        Arrow_FrequencyForBow = 1 / t
'    End Function
'
'    't = 2 * ps / v
'    Function Arrow_FrequencyForBow2(DrawLength As Double, BraceHeight As Double, ArrowSpeed As Double, ArrowWeight As Double) As Double
'        t = 2 * I2M(DrawLength - BraceHeight - 1.25) / (ArrowSpeed / 3 * 0.9144) '+ I2M(BraceHeight) / 2 / (ArrowSpeed / 3 * 0.9144)
'        Arrow_FrequencyForBow2 = 1 / t
'    End Function
'
''    Function Arrow_PenetrationDepth(ArrowWeight As Double, ArrowSpeed As Double, BroadHead_Ro As Double, BoneBreakage As Double) As Double
''        va = (ArrowSpeed ^ 2 - 450240 * BoneBreakage / ArrowWeight) ^ 0.5
''        va = va - 1500 / (va - 18)
''        Arrow_PenetrationDepth = ArrowWeight * va / 430 / BroadHead_Ro
''    End Function
'
''Function Arrow_PenetrationDepth2(m As Double, v As Double, Ro As Double, EkLoss As Double) As Double
''    'using v^2 = u^2 -2as, ie integrating over distance (1mm)
''    Dim s As Double
''    Dim vv As Double
''    Dim ss As Double
''    Dim a As Double
''    Dim k As Double
''    k = 0.001
''
''    s = 0.001  ' 1mm
''    If Ro < 2 Then Exit Function
''
''    v = m * v ^ 2 / 450240 - EkLoss
''    If v < 0 Then Exit Function
''    v = (v * 450240 / m) ^ 0.5
''
''    v = v * 12 * 2.54 / 100 'm/s
''    vv = v ^ 2 'm.m / s.s
''    ss = 0 ' m.m
''    Ro = Ro / 2.2 * 9.8 'Newton
''    m = m * 0.065 / 1000 'kg
''    While vv > 0.1
''        a = (Ro + k * Ro * vv) / m
''        vv = vv - 2 * a * s
''        ss = ss + s
''    Wend
''    If ss >= 0 Then Arrow_PenetrationDepth2 = ss * 100 / 2.54 'as inhces
''End Function
'
'Public Function Arrow_GetLaunchAngle(Distance As Double, Velocity As Double, Height As Double) As Double
'    s = Distance * 0.9144
'    v = Velocity * 12 * 2.54 / 100
'    H = Height * 0.9144
'    a = 0.01
'    n = 0
'    Delta = 0.01
'    dd = 0
'    Do
'        vx = Cos(a) * v
'        vy = Sin(a) * v
'        t = s / vx
'        y = vy * t - 9.8 / 2 * t ^ 2
'        If Abs(y - H) < 0.0001 Then
'            Arrow_GetLaunchAngle = a
'            Exit Function
'        ElseIf y > H Then
'            If dd > 0 Then Delta = Delta / 10
'            dd = -1
'            a = a - Delta
'        Else
'            If dd < 0 Then Delta = Delta / 10
'            dd = 1
'            a = a + Delta
'        End If
'        n = n + 1
'    Loop Until n > 1000
'    Arrow_GetLaunchAngle = a
'End Function
'
'Function Arrow_RangeError(Distance As Double, Velocity As Double, Height As Double, DistanceErr As Double) As Double
'Dim a As Double
'Dim b As Double
'Dim vx As Double
'Dim vy As Double
'Dim t As Double
'Dim g As Double
'Dim dy As Double
'    a = Arrow_GetLaunchAngle(Distance, Velocity, Height)
'    b = Atn(Height / (Distance + DistanceErr)) - Atn(Height / Distance)
'    vx = Cos(a + b) * (Velocity / 3)
'    vy = Sin(a + b) * (Velocity / 3)
'    t = (Distance + DistanceErr) / vx
'    g = -9.8 / 0.9144
'    dy = vy * t + g / 2 * t ^ 2 - Height
'    Arrow_RangeError = dy * 36
'End Function
'
'Public Sub Initialize()
'    Et = 177.6723
'    Mf = 397.5365
'    Kf = 0.000000818
'    ArrowStraightness = 0.006 * 25.4 'mm
'    NockTravel = 3 'mm
'    RestMisalignment = 0.5
'    NockWeight = 10
'    InsertWeight = 15
'    FlightWeight = 40
'End Sub
'
'Public Function Arrow_SetVMKmodel(Weight1 As Double, Speed1 As Double, Weight2 As Double, Speed2 As Double, Weight3 As Double, Speed3 As Double, selection As Integer) As Double
'    If Speed1 * Weight1 > 0 And Speed2 * Weight2 > 0 And Speed3 * Weight3 > 0 Then
'        'uses guass elimination to calculate unknowns in Et=(m+Mf+Kf.m^3)v^2/450240
'        Dim t(7, 3)
'        t(1, 1) = Weight1 ^ 3 * Speed1 ^ 2
'        t(1, 2) = Speed1 ^ 2
'        t(1, 3) = Speed1 ^ 2 * Weight1
'        t(2, 1) = Weight2 ^ 3 * Speed2 ^ 2
'        t(2, 2) = Speed2 ^ 2
'        t(2, 3) = Speed2 ^ 2 * Weight2
'        t(3, 1) = Weight3 ^ 3 * Speed3 ^ 2
'        t(3, 2) = Speed3 ^ 2
'        t(3, 3) = Speed3 ^ 2 * Weight3
'        t(4, 1) = t(1, 1) - t(2, 1)
'        t(4, 2) = t(1, 2) - t(2, 2)
'        t(4, 3) = t(1, 3) - t(2, 3)
'        t(5, 1) = t(3, 1) - t(1, 1)
'        t(5, 2) = t(3, 2) - t(1, 2)
'        t(5, 3) = t(3, 3) - t(1, 3)
'        t(6, 2) = t(5, 2) / t(5, 1) * t(4, 1)
'        t(6, 3) = t(5, 3) / t(5, 1) * t(4, 1)
'        t(7, 2) = t(6, 2) - t(4, 2)
'        t(7, 3) = t(6, 3) - t(4, 3)
'        Mf = -t(7, 3) / t(7, 2)
'        Kf = (-t(6, 3) - Mf * t(6, 2)) / t(4, 1)
'        Et = (-t(1, 3) - t(1, 2) * Mf - t(1, 1) * Kf) / -450240
'    ElseIf Speed1 * Weight1 > 0 And Speed2 * Weight2 > 0 Then
'        Mf = (Weight2 * Speed2 ^ 2 - Weight1 * Speed1 ^ 2) / (Speed1 ^ 2 - Speed2 ^ 2)
'        Kf = 0
'        Et = Weight1 * Speed1 ^ 2 / 450240
'    ElseIf Speed1 * Weight1 > 0 Then
'        Mf = 0
'        Kf = 0
'        Et = Weight1 * Speed1 ^ 2 / 450240
'    Else
'        Mf = 0
'        Kf = 0
'        Et = 0
'    End If
'    Select Case selection
'    Case 1: Arrow_SetVMKmodel = Mf
'    Case 2: Arrow_SetVMKmodel = Kf
'    Case 3: Arrow_SetVMKmodel = Et
'    End Select
'End Function
'
'Public Function Arrow_VMK(Et, Mf, Kf, Optional speed, Optional Weight) As Double
'    If IsMissing(speed) Then
'        Arrow_VMK = (Et * 450240 / (Weight + Mf + Kf * Weight ^ 3)) ^ 0.5
'    ElseIf IsMissing(Weight) Then
'        'Wt + Kf * Wt ^ 3 = Et * 450240 / speed ^ 2 - Mf
'        Dim R  As Double
'        R = Et * 450240 / speed ^ 2 - Mf
'        Dim H As Double, L As Double, m As Double, dm As Double
'        H = 1500
'        L = 100
'        Do
'            m = (H + L) / 2
'            dm = (m + Kf * m ^ 3) - R
'            If dm > 0 Then H = m
'            If dm < 0 Then L = m
'        Loop Until Abs(dm) < 0.01
'        Arrow_VMK = Int(m * 100) / 100
'    End If
'End Function
'
'Public Function Arrow_Get_GPI(PointWeight As Double, TotalWeight As Double, ShaftLength As Double) As Double
'    If FlightWeight = 0 Then Initialize
'    Arrow_Get_GPI = (TotalWeight - PointWeight - InsertWeight - NockWeight - FlightWeight - TubeWeight) / ShaftLength
'End Function

'Public Function Arrow_Deflection(PointWeight As Double, gpi As Double, L As Double, Spine As Double, BH As Double, DL As Double, Optional fps As Double = 0) As Double
'    Dim dt As Double, Ws As Double, Wp As Double, Wt As Double
'    Dim Fd As Double, Fg As Double, Fa As Double
'    Dim d As Double, eF As Double, v As Double, m As Double, t As Double, ct As Double, R As Double
'    Dim ps As Double, a As Double
'    Dim TubeLen As Double
'    Dim WeightTube As Double
'
'    If FlightWeight = 0 Then Initialize
'
'
'    'APPROXIMATION: Assume constant acceleration
'    dt = 0.01 / 1000 'milliseconds - time slice
'
'    InsertWeight = Range("InsertWeight")
'    FlightWeight = Range("FlightsWeight")
'    InsertWeight = Range("InsertWeight")
'    TubeWeight = Range("WeightTubeLength") * Range("WeightTubeGPI")
'    TubeLen = Range("WeightTubeLength")
'
'
'    Wt = PointWeight + InsertWeight + NockWeight + FlightWeight + L * gpi + WeightTube
'    Wp = PointWeight + InsertWeight
'    Ws = Wt - NockWeight - FlightWeight - Wp - WeightTube
'
'    'vmk model, works in imperic
'    If fps = 0 Then fps = (Et * 450240 / (Wt + Mf + Kf * Wt ^ 3)) ^ 0.5
'
'    'change the rest to metric units
'    Wt = Wt * 0.065 / 1000 'kg
'    Wp = Wp * 0.065 / 1000 'pile weight
'    Ws = Ws * 0.065 / 1000 'shaft weight
'    DL = DL * 2.54 / 100 'm
'    BH = BH * 2.54 / 100 'm
'    L = L * 2.54 / 100 'm
'    v = fps * 12 * 25.4 / 1000 'm/s - arrow speed
'
'    'the max G force the arrow experience as axial pressure F=ma
'    ps = DL - BH - (1.25 * 2.54 / 100)
'    Dim Li As Double
'    Li = L - (TubeLength * 2.54 / 100) / 2
'    Fg = (Wp + Ws / 2 + (TubeWeight * 0.065 / 1000) * (Li / L)) * v ^ 2 / (2 * ps)
'
'    'the mass that resists the buckling is 2/3 the shaft weight:
'    'the shaft and tube weight is factored as the integration of a parabole approx. the bend of the arrow
'    Li = L - (TubeLength * 2.54 / 100)
'    m = Ws * 2 / 3 + (TubeWeight * 0.065 / 1000) * (4 / 3 * (Li / L) ^ 3 - 2 * (Li / L) ^ 2 + 2 / 3)
'
'    'get the time the arrow spent accelerating from v = at and vv=2as
'    t = 2 * ps / v
'    If t > 1 Then Exit Function
'
'    Dim s0 As Double
'    s0 = ArrowStraightness * L + NockTravel + RestMisalignment
'    s0 = s0 / 1000 'm
'
'    Dim EI48 As Double 'Pascal, or N/m^2
'    EI48 = (10 * 0.88) * (28 * 2.54 / 100) ^ 3 / (Spine / 1000 * 2.54 / 100) 'D = FL^3/(EI48)
'    Debug.Print EI48 / Spine
'    d = s0
'    v = 0
'    ct = 0
'
'    Do
'        'fraction of axial force available for bending
'        R = d / ((L / 2) ^ 2 - d ^ 2) ^ 0.5
'
'        'the force the arrow is resisiting the bending
'        Fa = (d - s0) * EI48 / L ^ 3 'initial bending is without force
'
'        eF = Fg * R - Fa
'
'        a = eF / m
'        d = d + v * dt + a * dt ^ 2 / 2
'        v = v + a * dt
'        ct = ct + dt
'    Loop While ct < t And d < L / 2
'
'    Arrow_Deflection = d * 1000 'mm
'End Function

Function Arrow_PenetrationDepth3( _
    ByVal Ro As Double, _
    ByVal m As Double, _
    ByVal v As Double, _
    ByVal dia As Double, _
    ByVal BHArea As Double, _
    ByVal ds As Double, _
    ByVal KELoss1 As Double, _
    ByVal BodyDepth As Double, _
    ByVal KELoss2 As Double)
    
    Dim s As Double, d As Double, PI As Double, KE As Double, vv As Double
    Arrow_PenetrationDepth3 = "err"
    PI = 3.141592654
        
    s = 0 'initialize distance penetrated
    If Ro < 2 Then Exit Function 'too little Ro, model will crash
    If ds <= 0 Then ds = 0.001 '1 mm steps
        
    'limits the speed to about 140 fps
    v = v - ((450 / v) ^ 6 - 1)
    If v < 0 Then Exit Function
    
    'SI units
    Ro = Ro / 2.2 * 9.8 'lbs to Newton
    m = m * 0.065 / 1000 'grains to kilogram
    v = v * 12 * 2.54 / 100 'fps to meter per second
    dia = dia * 2.54 / 100 'inches to meter
    circ = PI * dia 'shaft circumference in meter
    BHArea = BHArea / 1550.0031 'inches squared to meter squared
    BodyDepth = BodyDepth * 2.54 / 100 'inches to m
    KELoss1 = KELoss1 * 1.359431604   'foot pounds to J
    KELoss2 = KELoss2 * 1.359431604   'foot pounds to J
    
    'constants
    pt = 1400#  'denisty of soft tissue kg/m3
    pb = 1055#  'density of blood kg/m3
    u = 0.0035  'viscosity of blood
    
    If KELoss1 > 0 Then 'modeling bone penetration as loss of energy
    ' this is valid as penetration into bone is related to KE
        KE = m * v ^ 2 / 2
        KE = KE - KELoss1
        If KE < 0 Then Exit Function  'if we don't get through the bone, drop the simulation
        v = (KE * 2 / m) ^ 0.5 'get the speed exiting at back of bone
    End If
    
    vv = v ^ 2  'meter per second squared - we will loop using vv
    s = 0
    L = ds
    Dim ShaftArea As Double
    Do While vv > 1 And s < 10
        If s > BodyDepth + 0.75 Then            '      [       ]  ------------->
            L = 0
            Ro = 0
            BHArea = 0
            ShaftArea = 0
            Arrow_PenetrationDepth3 = 999
            Exit Function
        ElseIf s > BodyDepth And s > 0.75 Then  '      [   ----]------->
            L = (0.75 - (s - BodyDepth))
            ShaftArea = L * circ
            BHArea = 0
            Ro = 0
        ElseIf s > BodyDepth Then               '    --[-------]--->
            L = BodyDepth
            ShaftArea = L * circ
            BHArea = 0
            Ro = 0
            If KELoss2 > 0 Then  ''exit rib breakage energy loss
                KE = m * vv / 2 - KELoss2
                If KE < 0 Then                  '    --[-------]
                    'if we don't get through the bone, drop the simulation
                    'StoppedByRib = True
                    Exit Do
                End If
                vv = (KE * 2 / m)
                KELoss2 = 0
            End If
        ElseIf s > 0.75 Then                    '      [  ------------>       ]
            L = 0.75
            ShaftArea = L * circ
        Else                                    ' ------[----> ]
            L = s
            ShaftArea = L * circ
        End If
        
        v = vv ^ 0.5
        FrontalArea = (dia / 2) ^ 2 * PI
        Fd = 0.5 * (Ro / 200) * FrontalArea * pt * vv
        Fs = 0.02915 * ((BHArea + ShaftArea) / ((pb * v * (L + ds) / u) ^ 0.2)) * pb * vv
        Ft = Ro + Fd + Fs
        'd = (Ro + (Ro * a + aa) * vv) / m
        d = Ft / m
        vv = vv - 2 * d * ds
        s = s + ds
    Loop
    If s >= 10 Then
        Arrow_PenetrationDepth3 = 999
    Else
        Arrow_PenetrationDepth3 = s * 100 / 2.54 'and return the result in inches
    End If
End Function

Function Arrow_ReqFPS( _
    ByVal Ro As Double, _
    ByVal m As Double, _
    ByVal RD As Double, _
    ByVal dia As Double, _
    ByVal a As Double, _
    ByVal ds As Double, _
    ByVal KELoss1 As Double, _
    ByVal BodyWidth As Double, _
    ByVal EKLoss2 As Double, _
    ByVal ReqDepth As Double)
    
    Dim v As Double
    v = 50
    Dim R
    R = Arrow_PenetrationDepth(Ro, m, v, dia, a, ds, KELoss1, BodyWidth, EKLoss2)
    While v < 350 And (R = "err" Or R <= RD)
        v = v + 1
        R = Arrow_PenetrationDepth(Ro, m, v, dia, a, ds, KELoss1, BodyWidth, EKLoss2)
    Wend
    If v < 350 Then Arrow_ReqFPS = v Else Arrow_ReqFPS = "n.a."
End Function


Function Arrow_Ro(BladeCount As Double, CD As Double, L As Double, Ac As Double, At As Double, dia As Double, coc As Integer)
    'CD - cutting diameter
    'L - length of the blade's cutting edge
    'At - tip angle
    'Ac - cutting edge angle
    
    At = At / 180 * 3.14
    Ac = Ac / 180 * 3.14
    
    Dim x As Double 'the base length of the blade
    Dim H As Double 'blade height
    H = (CD - dia) / 2
    x = (L ^ 2 - H ^ 2) ^ 0.5
    
    Dim Atc As Double
    If coc = 1 Then Atc = Tan(Ac / 2) Else Atc = 1

    Dim St As Double 'frontal surfact area
    St = (dia / 2) ^ 2 * 3.14
    
    At = Tan(At / 2)
    Ac = Tan(Ac / 2)
    
    Arrow_Ro = 44 * (St * At * Atc) ^ 0.5 + 44 * BladeCount * H ^ 2 / (L + x) * Ac
End Function