' Reference: NASA/TM—2005-213658 - " Parametric (On-Design) Cycle Analysis for a
' Separate-Exhaust Turbofan Engine With Interstage Turbine Burner" by
' Liew, K. H.; Urip, E.; Yang, S. L.; Siow, Y. K.; Marek, C. J 
Sub itb()
  Dim rc, rb, rt, a0, v0, eta_r, tau_alp_b, tau_alp_itb
  Dim tau_c, tau_b, tau_itb
  Dim T19_T0, fc_mc, ff_mf, denu, nume, v19, v9
  '
  Style = vbYes
  Style1 = vbYes
  ' Msg = "Do you want to jump to the next alpha?"
  Msg1 = "Division by zero ==> Check zero inputs for gama's (eqn 1,2,3)"
  Msg2 = "Imaginary sound speed ==> Check negative inputs for gamac, gc or T0 (eqn 4)"
  Msg3 = "Division by zero ==> gamac cannot = 1 (eqn 7)"
  Msg4 = "Division by zero ==> Check zero inputs for cpc, T0, ecH, ecL & ef (eqn 10, _
    11, 12, 13 & 14)"
  Msg5 = "Division by zero ==> Check zero inputs for tau_r & tau_c (eqn 16)"
  Msg6 = "Division by zero ==> tau_cH, tau_cL & tau_f cannot = 1 (eqn 17, 18, 19)"
  Msg7 = "Invalid (negative) tau_tH ==> Refer to equation 21"
  Msg8 = "Imaginary M9 ==> Pt9/P9 < 1 (eqn 31 & 34)"
  Msg9 = "Imaginary M19 ==> Pt19/P19 < 1 (eqn 36 & 37)"
  Msg10 = "Division by zero => denominator of 'fb' eqn cannot = 0 (eqn 20)"
  Msg11 = "Division by zero => denominator of 'fitb' eqn cannot = 0 (eqn 22)"
  wdata = 1
  fbi = 0
  fitbi = 0
  'warn = 0
  '
  'If gamac = 0# Or gamat = 0# Or gamaitb = 0# Then
    'response = MsgBox(Msg1, style, Title)
    'Msg1 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", Mo=" & M0 & ", 
      alpha=" & alpha
    'response = MsgBox(Msg1, Style, Title)
    'If singlept <> 1 Then jump = MsgBox(Msgsheet, Style1, Title)
    'GoTo 10
  'End If
  '
  'added on 5/13/04 for controling ITB switch -1
  itbinit = 1
  If Sheets("Input").ComboBox2.Value = "OFF" Then itbinit = 0
    
  'added on 5/13/04 for controling ITB switch -2
  'Constants
  Call enthalpy(1, T_o, 0, 0, Cpc, gamac)
  Call enthalpy(1, Tt4, 0, 0, Cpt, gamat)
  rc = (gamac - 1#) * Cpc / gamac           '(1)
  rt = (gamat - 1#) * Cpt / gamat           '(2)
  '
  If gamac < 0# Or rc < 0# Or T_o < 0# Or gc < 0# Then
    response = MsgBox(Msg2, Style, Title)
    'Msg2 = "Occurs at " & "PIc,PIf,Mo,ByPass = " & pi_cH * pi_cL & "," & _
      pi_f & "," & M0 & "," & alpha
    'Msg2 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", _
      Mo=" & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
     'response = MsgBox(Msg2, Style, Title)
    'warn = 1
    'If singlept <> 1 Or last = 0 Then jump = MsgBox(Msgsheet, Style1, Title)
    GoTo 10
  End If
  a0 = (gamac * rc * T_o * gc) ^ 0.5        '(4)
  v0 = a0 * M0                              '(5)
  '
  If gamac = 1# Then
    response = MsgBox(Msg3, Style, Title)
    'Msg3 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", _
      Mo=" & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
    'response = MsgBox(Msg3, Style, Title)
    'warn = 1
    'If singlept <> 1 Then jump = MsgBox(Msgsheet, Style1, Title)
     GoTo 10
  End If
  '
  tau_r = (gamac - 1#) * (0.5 * M0 ^ 2) + 1# '(6)
  pi_r = tau_r ^ (gamac / (gamac - 1))       '(7)
  '
  If M0 <= 1# Then                           '(8)
    eta_r = 1#
  End If
  
  If M0 > 1 And M0 <= 5 Then
    eta_r = 1# - 0.075 * (M0 - 1#) ^ 1.35
  End If
  If M0 > 5 Then
    eta_r = 800# / (935# + M0 ^ 4)
  End If
  pi_d = pi_dmax * eta_r                     '(9)
  '
  If pi_cL < pi_f Then pi_cL = pi_f
  pi_cH = pi_c / pi_cL
  '
  If Cpc = 0# Or T_o = 0# Or gamac = 0# Or e_cH = 0# Or e_cL = 0# Or e_f = 0# Then
      response = MsgBox(Msg4, Style, Title)
      'Msg4 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", _ 
        Mo=" & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
      'response = MsgBox(Msg4, Style, Title)
      'If singlept <> 1 Then jump = MsgBox(Msgsheet, Style1, Title)
      'warn = 1
      GoTo 10
  End If
  24
  tau_alp_b = Cpt * Tt4 / (Cpc * T_o)              '(10)
  tau_cH = pi_cH ^ ((gamac - 1#) / (gamac * e_cH)) '(12)
  tau_cL = pi_cL ^ ((gamac - 1#) / (gamac * e_cL)) '(13)
  tau_f = pi_f ^ ((gamac - 1#) / (gamac * e_f))    '(14)
  tau_c = tau_cH * tau_cL                          '(15)
  '
  If tau_r = 0# Or tau_c = 0# Then
      response = MsgBox(Msg5, Style, Title)
      'Msg5 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", Mo=" e_cH
        & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
      'response = MsgBox(Msg5, Style, Title)
      'warn = 1
      'If singlept <> 1 Then jump = MsgBox(Msgsheet, Style1, Title)
      GoTo 10
  End If
  tau_b = Tt4 / (T_o * tau_r * tau_c)             '(16)
  '
  If tau_cH = 1# Or tau_cL = 1# Or tau_f = 1# Then
      response = MsgBox(Msg6, Style, Title)
      'Msg6 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", _
        Mo=" & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
      'response = MsgBox(Msg6, Style, Title)
      'warn = 1
      'If singlept <> 1 Then jump = MsgBox(Msgsheet, Style1, Title)
    GoTo 10
  End If
  eta_cH = ((tau_cH ^ e_cH) - 1#) / (tau_cH - 1#) '(17)
  eta_cL = ((tau_cL ^ e_cL) - 1#) / (tau_cL - 1#) '(18)
  eta_f = ((tau_f ^ e_f) - 1#) / (tau_f - 1#)     '(19)
  '
  'add MSH on 7/30/04 -1
  'Component: MAIN BURNER
  h_0 = Cpc * T_o
  temp_fb = eta_b * hPRb / h_0 - tau_alp_b
  If temp_fb = 0# Then
      response = MsgBox(Msg10, Style, Title)
      'Msg7 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", _ 
        Mo=" & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
      'response = MsgBox(Msg7, Style, Title)
      'If singlept <> 1 Then jump = MsgBox(Msgsheet, Style1, Title)
      'warn = 1
    GoTo 10
  End If
  'Fuel model 1: CSH
  If Sheets("Input").ComboBox4.Value = "CSH" Then
      fb = (tau_alp_b - tau_r * tau_c) / temp_fb       '(20a)
  Else
    'Fuel model 2: MSH
      loop3 = 0
      Tt3 = T_o * tau_r * tau_c                        'eqn(20b)
      Call enthalpy(0, Tt3, 0, h_t3, 0, 0)
      f_4i = (tau_alp_b - tau_r * tau_c) / temp_fb     'eqn(20c),initial value for f_4i
  3
    If loop3 > 20 Then
      GoTo 10
    End If
    Call enthalpy(0, Tt4, f_4i, h_t4, 0, 0)
    temp_fb = eta_b * hPRb - h_t4
    If temp_fb = 0# Then
      response = MsgBox(Msg10, Style, Title)
      'Msg7 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", _ 
        Mo=" & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
      'If singlept <> 1 Then jump = MsgBox(Msg_jump, Style1, Title)
      'warn = 1
      GoTo 10
    End If
    fb = (h_t4 - h_t3) / temp_fb                     'eqn(20d)
    If Abs(fb - f_4i) > 0.0001 Then
      f_4i = fb
      loop3 = loop3 + 1
      GoTo 3
    End If
  End If

  'lkh1 - added cp=cp(T,f) and gama=gama(T,f) on 11/15/04
  Err_fb = Abs(fbi - fb)
  If Err_fb > 0.0001 Then
    fbi = fb
    Call enthalpy(1, Tt4, fb, 0, Cpt, gamat)
    GoTo 24
  End If
  'lkh2
  'add MSH on 7/30/04 -2
  '
  tau_tH = 1# - (tau_cH - 1#) * tau_cL * tau_r / ((1# + fb) * tau_alp_b * effm_tH) '5/11/04 '(21)
  If tau_tH < 0# Then
    response = MsgBox(Msg7, Style, Title)
    'Msg7 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", _ 
      Mo=" & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
    'response = MsgBox(Msg7, Style, Title)
    Cells(22, 4).Select
    Msg7 = "The Tt4 lower limit (as indicated) is too low to converge, increase _
      Tt4 lower limit and try again!"
    response = MsgBox(Msg7, Style, Title)
    'warn = 1
    GoTo 10
  End If
  If itbinit = 0 Then
    Tt4p5 = Tt4 * tau_tH
    Call enthalpy(1, Tt4p5, 0, 0, Cpitb, gamaitb)
    tau_alp_itb = Cpitb * Tt4p5 / (Cpc * T_o)              '(11)
    fitb = 0#
    tau_itb = 1#
    GoTo 5
  End If
  Call enthalpy(1, Tt4p5, 0, 0, Cpitb, gamaitb)

  25
  tau_alp_itb = Cpitb * Tt4p5 / (Cpc * T_o)
  If tau_alp_itb <= (tau_alp_b * tau_tH) Then
    itbinit = 0
    Tt4p5 = Tt4 * tau_tH
    Call enthalpy(1, Tt4p5, 0, 0, Cpitb, gamaitb)
    fitb = 0#
    tau_itb = 1#
    tau_alp_itb = Cpitb * Tt4p5 / (Cpc * T_o)
    GoTo 5
  End If
  '
  If itbinit = 1 Then
    'add MSH on 8/1/04 for ITB -1
    'Component: INTERSTAGE-TURBINE BURNER
    temp_fitb = eta_itb * hPRitb / h_0 - tau_alp_itb
    If temp_fitb = 0# Then
      response = MsgBox(Msg11, Style, Title)
      'Msg7 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", _
        Mo=" & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
      'response = MsgBox(Msg7, Style, Title)
      'If singlept <> 1 Then jump = MsgBox(Msg_jump, Style1, Title)   
      'warn = 1
      GoTo 10
    End If
    'Fuel model 1: CSH
    If Sheets("Input").ComboBox4.Value = "CSH" Then
      fitb = (tau_alp_itb - tau_alp_b * tau_tH) / (temp_fitb) * (1# + fb) '9/3/04   '(22a)
    Else
      'Fuel model 2: MSH
      loop4 = 0
      Tt4p4 = Tt4 * tau_tH  '(22b)
      Call enthalpy(0, Tt4p4, fb, h_t4p4, 0, 0)
      f_4p5i = (1# + fb) * (tau_alp_itb - tau_alp_b * tau_tH) / (temp_fitb)   '(22c),initial value for f_4p5i
  4
      If loop4 > 20 Then
        loop4 = 21
        GoTo 10
      End If
      Call enthalpy(0, Tt4p5, f_4p5i, h_t4p5, 0, 0)
      temp_fitb = eta_itb * hPRitb - h_t4p5
    If temp_fitb = 0# Then
      response = MsgBox(Msg11, Style, Title)
      'Msg7 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", _
        Mo=" & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
      'response = MsgBox(Msg7, Style, Title)
      'warn = 1
      'If singlept <> 1 Then jump = MsgBox(Msg_jump, Style1, Title)
      GoTo 10
    End If
    fitb = (1# + fb) * (h_t4p5 - h_t4p4) / (temp_fitb) '(22d)
    If Abs(fb + fitb - f_4p5i) > 0.0001 Then
      f_4p5i = fb + fitb
      loop4 = loop4 + 1
      GoTo 4
    End If
  End If
  'lkh1 - added cp=cp(T,f) and gama=gama(T,f) on 11/15/04
  Err_fitb = Abs(fitbi - fitb)
  If Err_fitb > 0.0001 Then
    fitbi = fitb
    Call enthalpy(1, Tt4p5, fitb, 0, Cpitb, gamaitb)
    GoTo 25
  End If
  'lkh2
  'add MSH on 8/1/04 for ITB -2
  tau_itb = Tt4p5 / (T_o * tau_r * tau_c * tau_b * tau_tH) '(24)
  End If
  5
  ritb = (gamaitb - 1#) * Cpitb / gamaitb                  '(3)
  f = fb + fitb                                            '(23)
  tau_tL = 1# - (tau_cL - 1# + alpha * (tau_f - 1#)) * tau_r / ((1# + fb + fitb) e_cH
    * tau_alp_itb * effm_tL)                                     '(26)
  pi_tH = tau_tH ^ (gamat / ((gamat - 1#) * e_tH))         '(27)
  pi_tL = tau_tL ^ (gamaitb / ((gamaitb - 1#) * e_tL))     '(28) 
  eta_tH = (1# - tau_tH) / (1# - tau_tH ^ (1# / e_tH))     '(29) 
  eta_tL = (1# - tau_tL) / (1# - tau_tL ^ (1# / e_tL))     '(30)
  '
  Pt9_P9 = P0_P9 * pi_r * pi_d * pi_cL * pi_cH * pi_b * pi_tH * pi_itb * pi_tL * pi_n   '(31)
  tempM9 = Pt9_P9 ^ ((gamaitb - 1#) / gamaitb)
  If tempM9 < 1# Then
    response = MsgBox(Msg8, Style, Title)
    'Msg8 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", _
      Mo=" & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
    'response = MsgBox(Msg8, Style, Title)
    'warn = 1
    'If singlept <> 1 Then jump = MsgBox(Msgsheet, Style1, Title)
    GoTo 10
  End If
  Tt9_T0 = tau_r * tau_c * tau_b * tau_tH * tau_itb * tau_tL '(32)
  T9_T0 = Tt9_T0 / Pt9_P9 ^ ((gamaitb - 1#) / gamaitb)       '(33)
  '
  M9 = ((tempM9 - 1#) * 2# / (gamaitb - 1#)) ^ 0.5           '(34)'5/5/04
  M8 = M9
  If M9 >= 1 Then M8 = 1
  v9_a0 = M9 * (gamaitb * ritb * T9_T0 / (gamac * rc)) ^ 0.5 '(35)
  v9 = v9_a0 * a0
  v9_v0 = v9 / v0
  '
  Pt19_P19 = P0_P19 * pi_r * pi_d * pi_f * pi_fn             '(36)
  tempM19 = Pt19_P19 ^ ((gamac - 1#) / gamac)
  If gamac < 1# Or tempM19 < 1# Then
    response = MsgBox(Msg9, Style, Title)
    'Msg8 = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", _  
      Mo=" & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
    'response = MsgBox(Msg9, Style, Title)
    'warn = 1
    'If singlept <> 1 then jump = MsgBox(Msgsheet, Style1, Title)
    GoTo 10
  End If
  M19 = ((tempM19 - 1#) * 2# / (gamac - 1#)) ^ 0.5            '(37)
  M18 = M19
  If M19 >= 1 Then M18 = 1
  Tt19_T0 = tau_r * tau_f                                     '(38)
  T19_T0 = Tt19_T0 / Pt19_P19 ^ ((gamac - 1#) / gamac)        '(39)
  v19_a0 = M19 * T19_T0 ^ 0.5                                 '(40)
  v19 = v19_a0 * a0
  v19_v0 = v19 / v0
  '
  fc_mc = (1# + f) * v9_a0 - M0 + (1# + f) * (ritb / rc) * (T9_T0 / v9_a0) * _ 
    (1# - P0_P9) * (1# / gamac)  '(41)
  ff_mf = v19_a0 - M0 + (T19_T0 / v19_a0) * (1# - P0_P19) / gamac '(42)
  F_mo = (a0 / gc) * (fc_mc + alpha * ff_mf) / (1# + alpha) '(43)
  '
  unitchange = 1#
  If Sheets("Input").ComboBox1.Value = "SI" Then
    unitchange = 1000000#
  Else
    unitchange = 3600#
  End If
  s = f / ((1# + alpha) * F_mo) * unitchange             '(44)

  If F_mo < 0 Or s < 0 Then wdata = 0
    ''Efficiencies
    corenume = (1# + f) * v9 * v9 - v0 * v0              'corrected on 6/16/04
    fanume = alpha * (v19 * v19 - v0 * v0)
    nume = 0.5 / gc * (corenume + fanume)                '(45)
    denu = fb * hPRb * eta_b + fitb * hPRitb * eta_itb   '(46)
    '
    eta_thermal = nume / denu * 100#                     '(47)
    eta_P = (1# + alpha) * F_mo * v0 / nume * 100#       '(48)
    eta_O = (1# + alpha) * F_mo * v0 / denu * 100#       '(49)
    '
    'Mach = 0.3 at HP turbine exit
    M4p4 = 0.3
    Tt4p4 = tau_r * tau_c * tau_b * tau_tH * T_o
    T4p4 = Tt4p4 / (1# + 0.5 * (gamat - 1#) * M4p4 * M4p4)
    v4p4 = M4p4 * (gamat * rt * gc * T4p4) ^ 0.5
    dTitb = Tt4p4 * (tau_itb - 1#)
    GoTo 11
  10
    'If warn = 1 Then                                   'added on 11/30/04
    Msg = "Occurs at " & "PIc=" & pi_cH * pi_cL & ", PIf=" & pi_f & ", _ 
      Mo=" & M0 & ", alpha=" & alpha & ", Tt4=" & Tt4
    response = MsgBox(Msg, Style, Title)
    If last = 1 Then
      jump = 1
      GoTo 11
    End If
    If singlept <> 1 Then jump = MsgBox(Msgsheet, Style1, Title)
  End If
  11
End Sub