program combined_cycle
    implicit none
    integer :: i, iostat_val
    double precision :: T1, T3, rp, gamma_g, cp_g, eta_comp, eta_turb_g
    double precision :: P_boiler, T_steam, P_cond, eta_turb_s, eta_pump
    double precision :: T_pinch_min, m_gas, m_steam
    double precision :: T2s, T2a, T4s, T4a, T_exhaust
    double precision :: W_comp, W_turb_g, Wnet_brayton, Qin_brayton, eta_brayton
    double precision :: h1_s, h2s_s, h2a_s, h3_s, h4_s, s1_s
    double precision :: Wt_s, Wp_s, Wnet_rankine, Qin_rankine, eta_rankine
    double precision :: Q_HRSG, T_pinch_actual, T_stack
    double precision :: Wnet_total, Qin_total, eta_combined, heat_rate
    double precision :: T_sat_boil, hf_boil, hfg_boil, sf_boil, sfg_boil, vf_boil
    double precision :: T_sat_cond, hf_cond, hfg_cond, sf_cond, sfg_cond, vf_cond
    double precision :: hg_boil, sg_boil, cp_steam_sh, dT_sh
    double precision :: rp_i, eta_b_i, eta_c_i, T4a_i
    double precision :: x4s, h4s_s_val, h4a_s_val
    double precision, parameter :: Ru = 8.314462d0

    read(*,*,iostat=iostat_val) T1
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid compressor inlet temperature.'
        stop
    end if
    read(*,*,iostat=iostat_val) T3
    read(*,*,iostat=iostat_val) rp
    read(*,*,iostat=iostat_val) gamma_g
    read(*,*,iostat=iostat_val) cp_g
    read(*,*,iostat=iostat_val) eta_comp
    read(*,*,iostat=iostat_val) eta_turb_g
    read(*,*,iostat=iostat_val) P_boiler
    read(*,*,iostat=iostat_val) T_steam
    read(*,*,iostat=iostat_val) P_cond
    read(*,*,iostat=iostat_val) eta_turb_s
    read(*,*,iostat=iostat_val) eta_pump
    read(*,*,iostat=iostat_val) T_pinch_min
    read(*,*,iostat=iostat_val) m_gas
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all combined cycle inputs.'
        stop
    end if
    if (T1<=0.0d0.or.T3<=0.0d0.or.rp<=1.0d0) then
        write(*,*) 'ERROR: T1, T3 positive and rp > 1 required.'
        stop
    end if
    if (gamma_g<=1.0d0) gamma_g = 1.4d0
    if (cp_g<=0.0d0) cp_g = 1005.0d0
    if (eta_comp<=0.0d0.or.eta_comp>1.0d0) eta_comp = 0.86d0
    if (eta_turb_g<=0.0d0.or.eta_turb_g>1.0d0) eta_turb_g = 0.90d0
    if (eta_turb_s<=0.0d0.or.eta_turb_s>1.0d0) eta_turb_s = 0.88d0
    if (eta_pump<=0.0d0.or.eta_pump>1.0d0) eta_pump = 0.85d0
    if (T_pinch_min<=0.0d0) T_pinch_min = 15.0d0
    if (m_gas<=0.0d0) m_gas = 100.0d0
    if (P_boiler<=0.0d0) P_boiler = 8.0d6
    if (P_cond<=0.0d0) P_cond = 10.0d3

    ! ===================== BRAYTON (GAS TURBINE) =====================
    ! 1→2: Isentropic compression
    T2s = T1 * rp**((gamma_g-1.0d0)/gamma_g)
    T2a = T1 + (T2s - T1) / eta_comp

    ! 3→4: Isentropic expansion
    T4s = T3 / rp**((gamma_g-1.0d0)/gamma_g)
    T4a = T3 - eta_turb_g * (T3 - T4s)

    W_comp = cp_g * (T2a - T1)          ! J/kg gas
    W_turb_g = cp_g * (T3 - T4a)        ! J/kg gas
    Wnet_brayton = W_turb_g - W_comp
    Qin_brayton = cp_g * (T3 - T2a)
    eta_brayton = Wnet_brayton / max(Qin_brayton, 1.0d-10)

    T_exhaust = T4a   ! gas turbine exhaust temperature

    ! ===================== HRSG =====================
    ! Saturation properties at steam boiler and condenser pressures
    call sat_props(P_boiler, T_sat_boil, hf_boil, hfg_boil, sf_boil, sfg_boil, vf_boil)
    call sat_props(P_cond, T_sat_cond, hf_cond, hfg_cond, sf_cond, sfg_cond, vf_cond)

    ! Superheated steam state at turbine inlet
    hg_boil = hf_boil + hfg_boil
    sg_boil = sf_boil + sfg_boil
    cp_steam_sh = 2100.0d0   ! J/(kg K) approximate
    dT_sh = T_steam - T_sat_boil
    if (dT_sh < 0.0d0) dT_sh = 0.0d0
    h1_s = hg_boil + cp_steam_sh * dT_sh
    s1_s = sg_boil + cp_steam_sh * log(max(T_steam/T_sat_boil, 1.0d0))

    ! Feedwater state (saturated liquid at condenser pressure, pumped)
    h4_s = hf_cond
    h2a_s = hf_cond + vf_cond * (P_boiler - P_cond) / eta_pump

    ! HRSG energy balance: m_gas * cp_g * (T_exhaust - T_stack) = m_steam * (h1_s - h2a_s)
    ! Pinch point constraint: T_stack >= T_sat_boil + T_pinch_min  (simplified)
    T_pinch_actual = T_exhaust - T_sat_boil
    T_stack = T_sat_boil + T_pinch_min

    if (T_stack > T_exhaust - 10.0d0) then
        T_stack = T_exhaust - 10.0d0
    end if
    if (T_stack < T1 + 20.0d0) T_stack = T1 + 20.0d0

    Q_HRSG = m_gas * cp_g * (T_exhaust - T_stack)
    Qin_rankine = h1_s - h2a_s
    m_steam = Q_HRSG / max(Qin_rankine, 1.0d-10)

    ! ===================== RANKINE (STEAM TURBINE) =====================
    ! Isentropic expansion to P_cond
    x4s = (s1_s - sf_cond) / max(sfg_cond, 1.0d-10)
    if (x4s > 1.0d0) then
        h4s_s_val = hf_cond + hfg_cond + cp_steam_sh * &
                    (T_sat_cond*exp((s1_s-sg_boil)/cp_steam_sh) - T_sat_cond)
    else
        h4s_s_val = hf_cond + x4s * hfg_cond
    end if
    h4a_s_val = h1_s - eta_turb_s * (h1_s - h4s_s_val)

    Wt_s = h1_s - h4a_s_val                         ! J/kg steam
    Wp_s = h2a_s - hf_cond                           ! J/kg steam
    Wnet_rankine = Wt_s - Wp_s
    eta_rankine = Wnet_rankine / max(Qin_rankine, 1.0d-10)

    ! ===================== COMBINED =====================
    Wnet_total = m_gas * Wnet_brayton + m_steam * Wnet_rankine   ! W total
    Qin_total = m_gas * Qin_brayton
    eta_combined = Wnet_total / max(Qin_total, 1.0d-10)
    heat_rate = 3600.0d0 / max(eta_combined, 1.0d-10)   ! kJ/kWh

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   COMBINED CYCLE (BRAYTON + RANKINE) ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- BRAYTON (GAS TURBINE) -----------------------------------'
    write(*,'(A,F12.2,A)')  '  T1 (compressor inlet)     = ', T1, ' K'
    write(*,'(A,F12.2,A)')  '  T2s (isentropic)          = ', T2s, ' K'
    write(*,'(A,F12.2,A)')  '  T2a (actual)              = ', T2a, ' K'
    write(*,'(A,F12.2,A)')  '  T3 (turbine inlet)        = ', T3, ' K'
    write(*,'(A,F12.2,A)')  '  T4s (isentropic)          = ', T4s, ' K'
    write(*,'(A,F12.2,A)')  '  T4a (exhaust)             = ', T4a, ' K'
    write(*,'(A,F12.2)')    '  Pressure Ratio rp         = ', rp
    write(*,'(A,F12.2,A)')  '  Compressor Work           = ', W_comp/1000, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Gas Turbine Work          = ', W_turb_g/1000, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Brayton Net Work          = ', Wnet_brayton/1000, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Brayton Heat Input        = ', Qin_brayton/1000, ' kJ/kg'
    write(*,'(A,F10.4)')    '  Brayton Efficiency        = ', eta_brayton
    write(*,*)
    write(*,'(A)') '--- HRSG (HEAT RECOVERY) ------------------------------------'
    write(*,'(A,F12.2,A)')  '  Gas Exhaust Temperature   = ', T_exhaust, ' K'
    write(*,'(A,F12.2,A)')  '  Stack Temperature         = ', T_stack, ' K'
    write(*,'(A,F12.2,A)')  '  Pinch Temperature Delta   = ', T_pinch_actual, ' K'
    write(*,'(A,F12.2,A)')  '  HRSG Heat Recovery        = ', Q_HRSG/1.0d6, ' MW'
    write(*,'(A,ES12.4,A)') '  Steam Mass Flow           = ', m_steam, ' kg/s'
    write(*,'(A,F12.4)')    '  Steam-to-Gas Mass Ratio   = ', m_steam/m_gas
    write(*,*)
    write(*,'(A)') '--- RANKINE (STEAM TURBINE) ---------------------------------'
    write(*,'(A,ES12.4,A)') '  Boiler Pressure           = ', P_boiler, ' Pa'
    write(*,'(A,F12.2,A)')  '  Steam Temperature         = ', T_steam, ' K'
    write(*,'(A,ES12.4,A)') '  Condenser Pressure        = ', P_cond, ' Pa'
    write(*,'(A,F12.2,A)')  '  Steam Turbine Work        = ', Wt_s/1000, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Pump Work                 = ', Wp_s/1000, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Rankine Net Work          = ', Wnet_rankine/1000, ' kJ/kg'
    write(*,'(A,F10.4)')    '  Rankine Efficiency        = ', eta_rankine
    write(*,*)
    write(*,'(A)') '--- COMBINED CYCLE PERFORMANCE ------------------------------'
    write(*,'(A,F12.2,A)')  '  Brayton Power             = ', m_gas*Wnet_brayton/1.0d6, ' MW'
    write(*,'(A,F12.2,A)')  '  Rankine Power             = ', m_steam*Wnet_rankine/1.0d6, ' MW'
    write(*,'(A,F12.2,A)')  '  Total Power               = ', Wnet_total/1.0d6, ' MW'
    write(*,'(A,F12.2,A)')  '  Total Heat Input          = ', Qin_total/1.0d6, ' MW'
    write(*,'(A,F10.4)')    '  Combined Efficiency       = ', eta_combined
    write(*,'(A,F10.2,A)')  '  Combined Efficiency       = ', eta_combined*100, ' percent'
    write(*,'(A,F10.2,A)')  '  Heat Rate                 = ', heat_rate, ' kJ/kWh'
    write(*,'(A,F10.4)')    '  eta_check (1-(1-eB)(1-eR))= ', &
        1.0d0-(1.0d0-eta_brayton)*(1.0d0-eta_rankine*Q_HRSG/(m_gas*Qin_brayton))
    write(*,*)

    ! Efficiency vs pressure ratio sweep
    write(*,'(A)') '--- EFFICIENCY VS PRESSURE RATIO SWEEP ----------------------'
    write(*,'(A)') '  rp            eta_Brayton   eta_Combined  T_exhaust[K]'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, 50
        rp_i = 2.0d0 + 38.0d0*dble(i-1)/49.0d0
        T2s = T1 * rp_i**((gamma_g-1.0d0)/gamma_g)
        T2a = T1 + (T2s - T1)/eta_comp
        T4s = T3 / rp_i**((gamma_g-1.0d0)/gamma_g)
        T4a_i = T3 - eta_turb_g*(T3 - T4s)
        eta_b_i = (cp_g*(T3-T4a_i) - cp_g*(T2a-T1)) / max(cp_g*(T3-T2a), 1.0d-10)
        ! Simplified combined
        if (T4a_i > T_stack) then
            eta_c_i = eta_b_i + (1.0d0-eta_b_i)*eta_rankine* &
                      (T4a_i-T_stack)/(T4a_i-T1)
        else
            eta_c_i = eta_b_i
        end if
        write(*,'(F10.2,2X,F10.5,2X,F10.5,2X,F12.2)') rp_i, eta_b_i, eta_c_i, T4a_i
    end do
    write(*,*)

    ! T-Q diagram for HRSG
    write(*,'(A)') '--- HRSG TQ DIAGRAM DATA ------------------------------------'
    write(*,'(A)') '  Q_frac        T_gas[K]      T_steam[K]'
    write(*,'(A)') '  -------------------------------------------'
    do i = 0, 50
        call hrsg_tq(dble(i)/50.0d0, T_exhaust, T_stack, &
                     T_sat_boil, T_steam, hf_cond, h2a_s, h1_s, hf_boil, hg_boil)
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Brayton: T2s=T1*rp^((g-1)/g); T4s=T3/rp^((g-1)/g).'
    write(*,'(A)') '  eta_combined ~ 1-(1-eta_B)(1-eta_R * Q_HRSG/Qin_B).'
    write(*,'(A)') '  HRSG: m_gas*cp*(T_exh-T_stack)=m_steam*(h1-h2a).'

contains

    subroutine sat_props(P_pa, Tsat, hf, hfg, sf, sfg, vf)
        implicit none
        double precision, intent(in) :: P_pa
        double precision, intent(out) :: Tsat, hf, hfg, sf, sfg, vf
        double precision :: P_MPa, logP
        P_MPa = P_pa / 1.0d6
        if (P_MPa < 0.001d0) P_MPa = 0.001d0
        logP = log(P_MPa)
        Tsat = 373.15d0 + 42.0d0*logP - 0.8d0*logP**2
        if (Tsat < 290.0d0) Tsat = 290.0d0 + 15.0d0*P_MPa
        if (Tsat > 647.0d0) Tsat = 647.0d0
        hf = (417.0d0 + 390.0d0*logP + 35.0d0*logP**2) * 1000.0d0
        if (hf < 100.0d3) hf = 100.0d3 + 200.0d3*P_MPa
        if (hf > 2100.0d3) hf = 2100.0d3
        hfg = (2258.0d0 - 180.0d0*logP - 40.0d0*logP**2) * 1000.0d0
        if (hfg < 200.0d3) hfg = 200.0d3
        if (hfg > 2500.0d3) hfg = 2500.0d3
        sf = (1.303d0 + 0.82d0*logP + 0.06d0*logP**2) * 1000.0d0
        if (sf < 0.3d3) sf = 0.3d3 + 0.5d3*P_MPa
        sfg = hfg / Tsat
        vf = 0.001d0 * (1.0d0 + 0.0002d0*P_MPa)
    end subroutine sat_props

    subroutine hrsg_tq(frac, T_exh, T_stk, T_sat, T_sh, hfc, h2a, h1, hfb, hgb)
        implicit none
        double precision, intent(in) :: frac, T_exh, T_stk, T_sat, T_sh
        double precision, intent(in) :: hfc, h2a, h1, hfb, hgb
        double precision :: Q_tot, Q_i, T_gas, T_stm
        double precision :: Q_eco, Q_evap, Q_sh, frac_eco, frac_evap
        Q_tot = h1 - h2a
        Q_eco = hfb - h2a
        Q_evap = hgb - hfb
        Q_sh = h1 - hgb
        frac_eco = Q_eco / max(Q_tot, 1.0d-10)
        frac_evap = Q_evap / max(Q_tot, 1.0d-10)
        Q_i = frac * Q_tot
        T_gas = T_exh - (T_exh - T_stk)*frac
        if (frac < frac_eco) then
            T_stm = h2a/max(hfb,1.0d-10)*T_sat * frac/max(frac_eco,1.0d-10)
            T_stm = T_sat * frac/max(frac_eco,1.0d-10) + 273.15d0*(1.0d0-frac/max(frac_eco,1.0d-10))
            if (T_stm < 300.0d0) T_stm = 300.0d0
        else if (frac < frac_eco + frac_evap) then
            T_stm = T_sat
        else
            T_stm = T_sat + (T_sh - T_sat)*(frac - frac_eco - frac_evap)/max(1.0d0-frac_eco-frac_evap,1.0d-10)
        end if
        write(*,'(F10.4,2X,F12.2,2X,F12.2)') frac, T_gas, T_stm
    end subroutine hrsg_tq

end program combined_cycle