Regenerative Rankine Cycle

program rankine_regenerative
    implicit none
    integer :: fwh_type, n_extractions, i, iostat_val
    double precision :: P_boiler, T_boiler, P_cond, P_ext1, P_ext2
    double precision :: eta_turbine, eta_pump
    double precision :: h1, s1, h2a, h2s, h3a, h3s, h4, h5, h6a, h6s
    double precision :: h_f_cond, h_fg_cond, s_f_cond, s_fg_cond, T_sat_cond
    double precision :: h_f_ext1, h_fg_ext1, s_f_ext1, s_fg_ext1, T_sat_ext1
    double precision :: h_f_ext2, h_fg_ext2, s_f_ext2, s_fg_ext2, T_sat_ext2
    double precision :: h_f_boil, s_f_boil, T_sat_boil, v_f_cond, v_f_ext1
    double precision :: x2, x3, y1, y2, Wt, Wp, Wnet, Qin, Qout
    double precision :: eta_th, bwr, m_dot, P_output, SSC
    double precision :: h_fw, h_mix, v_f_ext2
    double precision :: P_i, eta_i, y_i
    double precision, parameter :: Ru = 8.314462d0
    character(len=40) :: fwh_name

    read(*,*,iostat=iostat_val) P_boiler
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid boiler pressure input.'
        stop
    end if
    read(*,*,iostat=iostat_val) T_boiler
    read(*,*,iostat=iostat_val) P_cond
    read(*,*,iostat=iostat_val) n_extractions
    read(*,*,iostat=iostat_val) P_ext1
    read(*,*,iostat=iostat_val) P_ext2
    read(*,*,iostat=iostat_val) fwh_type
    read(*,*,iostat=iostat_val) eta_turbine
    read(*,*,iostat=iostat_val) eta_pump
    read(*,*,iostat=iostat_val) m_dot
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all Rankine inputs.'
        stop
    end if
    if (P_boiler <= 0.0d0 .or. P_cond <= 0.0d0) then
        write(*,*) 'ERROR: Pressures must be positive.'
        stop
    end if
    if (P_boiler <= P_cond) then
        write(*,*) 'ERROR: Boiler pressure must exceed condenser pressure.'
        stop
    end if
    if (eta_turbine <= 0.0d0 .or. eta_turbine > 1.0d0) eta_turbine = 0.88d0
    if (eta_pump <= 0.0d0 .or. eta_pump > 1.0d0) eta_pump = 0.85d0
    if (n_extractions < 1) n_extractions = 1
    if (n_extractions > 2) n_extractions = 2
    if (P_ext1 <= P_cond .or. P_ext1 >= P_boiler) P_ext1 = sqrt(P_boiler*P_cond)
    if (n_extractions == 2) then
        if (P_ext2 <= P_cond .or. P_ext2 >= P_ext1) &
            P_ext2 = sqrt(P_cond*P_ext1)
    end if
    if (m_dot <= 0.0d0) m_dot = 1.0d0

    if (fwh_type == 1) then
        fwh_name = 'Open feedwater heater (direct contact)'
    else
        fwh_name = 'Closed feedwater heater (shell-tube)'
        fwh_type = 2
    end if

    ! Compute saturation properties at key pressures
    call sat_props(P_cond, T_sat_cond, h_f_cond, h_fg_cond, s_f_cond, s_fg_cond, v_f_cond)
    call sat_props(P_ext1, T_sat_ext1, h_f_ext1, h_fg_ext1, s_f_ext1, s_fg_ext1, v_f_ext1)
    if (n_extractions == 2) then
        call sat_props(P_ext2, T_sat_ext2, h_f_ext2, h_fg_ext2, s_f_ext2, s_fg_ext2, v_f_ext2)
    end if
    call sat_props(P_boiler, T_sat_boil, h_f_boil, h_fg_cond, s_f_boil, s_fg_cond, v_f_cond)

    ! State 1: Turbine inlet (superheated steam at P_boiler, T_boiler)
    call superheat_props(P_boiler, T_boiler, T_sat_boil, h_f_boil, s_f_boil, h1, s1)

    ! State 2: After first extraction (isentropic expansion to P_ext1)
    call sat_props(P_ext1, T_sat_ext1, h_f_ext1, h_fg_ext1, s_f_ext1, s_fg_ext1, v_f_ext1)
    x2 = (s1 - s_f_ext1) / max(s_fg_ext1, 1.0d-10)
    if (x2 > 1.0d0) then
        ! Still superheated at extraction
        call superheat_from_entropy(P_ext1, s1, T_sat_ext1, h_f_ext1, s_f_ext1, s_fg_ext1, h2s)
    else
        h2s = h_f_ext1 + x2 * h_fg_ext1
    end if
    h2a = h1 - eta_turbine * (h1 - h2s)

    ! State 3: Turbine exit (isentropic expansion to P_cond)
    call sat_props(P_cond, T_sat_cond, h_f_cond, h_fg_cond, s_f_cond, s_fg_cond, v_f_cond)
    x3 = (s1 - s_f_cond) / max(s_fg_cond, 1.0d-10)
    if (x3 > 1.0d0) then
        call superheat_from_entropy(P_cond, s1, T_sat_cond, h_f_cond, s_f_cond, s_fg_cond, h3s)
    else
        h3s = h_f_cond + x3 * h_fg_cond
    end if
    h3a = h1 - eta_turbine * (h1 - h3s)

    ! State 4: Condenser exit (saturated liquid at P_cond)
    h4 = h_f_cond

    ! Pump work and extraction fraction calculation
    if (n_extractions == 1) then
        if (fwh_type == 1) then
            ! OPEN FWH: 1 extraction
            ! Pump 1: P_cond → P_ext1
            call sat_props(P_cond, T_sat_cond, h_f_cond, h_fg_cond, s_f_cond, s_fg_cond, v_f_cond)
            h5 = h4 + v_f_cond*(P_ext1 - P_cond)/eta_pump   ! after pump 1
            ! FWH energy balance: y*h2a + (1-y)*h5 = h_f_ext1
            h_fw = h_f_ext1
            y1 = (h_fw - h5) / max(h2a - h5, 1.0d-10)
            if (y1 < 0.0d0) y1 = 0.0d0
            if (y1 > 1.0d0) y1 = 1.0d0
            ! Pump 2: P_ext1 → P_boiler
            h6a = h_fw + v_f_ext1*(P_boiler - P_ext1)/eta_pump
            ! Turbine work per unit mass at boiler
            Wt = (h1 - h2a) + (1.0d0 - y1)*(h2a - h3a)
            ! Pump work
            Wp = (1.0d0 - y1)*(h5 - h4) + (h6a - h_fw)
            Qin = h1 - h6a
        else
            ! CLOSED FWH: 1 extraction
            h5 = h4 + v_f_cond*(P_boiler - P_cond)/eta_pump
            ! Closed FWH: extracted steam heats feedwater
            ! y*h2a + h5 = y*h_f_ext1 + h_fw (drain to trap)
            ! Assume drain cascaded back: h_fw = h5 + y*(h2a - h_f_ext1)
            y1 = (h_f_ext1 - h5) / max((h2a - h_f_ext1) + (h_f_ext1 - h5), 1.0d-10)
            if (y1 < 0.0d0) y1 = 0.0d0
            if (y1 > 1.0d0) y1 = 1.0d0
            h6a = h5 + y1*(h2a - h_f_ext1)
            Wt = (h1 - h2a) + (1.0d0 - y1)*(h2a - h3a)
            Wp = (h5 - h4)
            Qin = h1 - h6a
        end if
        y2 = 0.0d0
    else
        ! 2 extractions — open FWH simplified
        call sat_props(P_ext2, T_sat_ext2, h_f_ext2, h_fg_ext2, s_f_ext2, s_fg_ext2, v_f_ext2)
        ! Intermediate extraction state
        ! For simplicity: second extraction between ext1 and condenser
        ! Recompute h at P_ext2 from isentropic expansion
        x2 = (s1 - s_f_ext2) / max(s_fg_ext2, 1.0d-10)
        if (x2 > 1.0d0) then
            call superheat_from_entropy(P_ext2, s1, T_sat_ext2, h_f_ext2, &
                                        s_f_ext2, s_fg_ext2, h_mix)
        else
            h_mix = h_f_ext2 + x2 * h_fg_ext2
        end if
        h_mix = h1 - eta_turbine*(h1 - h_mix)  ! actual h at ext2

        ! Pump 1: P_cond → P_ext2
        h5 = h4 + v_f_cond*(P_ext2 - P_cond)/eta_pump
        ! Open FWH 1 at P_ext2
        y2 = (h_f_ext2 - h5) / max(h_mix - h5, 1.0d-10)
        if (y2 < 0.0d0) y2 = 0.0d0
        if (y2 > 0.5d0) y2 = 0.5d0

        ! Pump 2: P_ext2 → P_ext1
        h_fw = h_f_ext2 + v_f_ext2*(P_ext1 - P_ext2)/eta_pump
        ! Open FWH 2 at P_ext1
        y1 = (h_f_ext1 - h_fw) / max(h2a - h_fw, 1.0d-10)
        if (y1 < 0.0d0) y1 = 0.0d0
        if (y1 > 0.5d0) y1 = 0.5d0

        ! Pump 3: P_ext1 → P_boiler
        h6a = h_f_ext1 + v_f_ext1*(P_boiler - P_ext1)/eta_pump

        Wt = (h1 - h2a) + (1.0d0-y1)*(h2a - h_mix) + (1.0d0-y1-y2)*(h_mix - h3a)
        Wp = (1.0d0-y1-y2)*(h5-h4) + (1.0d0-y1)*(h_fw-h_f_ext2) + (h6a-h_f_ext1)
        Qin = h1 - h6a
    end if

    Wnet = Wt - Wp
    Qout = (1.0d0 - y1 - y2) * (h3a - h4)
    eta_th = Wnet / max(Qin, 1.0d-10)
    bwr = Wp / max(Wt, 1.0d-10)
    P_output = m_dot * Wnet / 1000.0d0    ! kW
    SSC = 3600.0d0 / max(Wnet, 1.0d-10)   ! kg/kWh

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   RANKINE CYCLE — REGENERATIVE ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Boiler Pressure           = ', P_boiler, ' Pa'
    write(*,'(A,F12.2,A)')  '  Boiler Temperature        = ', T_boiler, ' K'
    write(*,'(A,ES12.4,A)') '  Condenser Pressure        = ', P_cond, ' Pa'
    write(*,'(A,I8)')       '  Number of Extractions     = ', n_extractions
    write(*,'(A,ES12.4,A)') '  Extraction 1 Pressure     = ', P_ext1, ' Pa'
    if (n_extractions == 2) &
    write(*,'(A,ES12.4,A)') '  Extraction 2 Pressure     = ', P_ext2, ' Pa'
    write(*,'(A,A)')        '  FWH Type                  = ', trim(fwh_name)
    write(*,'(A,F10.4)')    '  Turbine Isentropic Eff    = ', eta_turbine
    write(*,'(A,F10.4)')    '  Pump Isentropic Eff       = ', eta_pump
    write(*,'(A,ES12.4,A)') '  Mass Flow Rate            = ', m_dot, ' kg/s'
    write(*,*)
    write(*,'(A)') '--- SATURATION DATA -----------------------------------------'
    write(*,'(A,F12.2,A)')  '  T_sat (boiler)            = ', T_sat_boil, ' K'
    write(*,'(A,F12.2,A)')  '  T_sat (ext1)              = ', T_sat_ext1, ' K'
    if (n_extractions == 2) &
    write(*,'(A,F12.2,A)')  '  T_sat (ext2)              = ', T_sat_ext2, ' K'
    write(*,'(A,F12.2,A)')  '  T_sat (condenser)         = ', T_sat_cond, ' K'
    write(*,*)
    write(*,'(A)') '--- STATE POINTS (per kg at boiler) -------------------------'
    write(*,'(A)') '  Point   h[kJ/kg]      s[kJ/kgK]     T[K]          P[Pa]'
    write(*,'(A)') '  ----------------------------------------------------------------'
    write(*,'(A,F12.2,2X,F10.4,2X,F12.2,2X,ES12.4)') '  1 turb in  ', h1/1000, s1/1000, T_boiler, P_boiler
    write(*,'(A,F12.2,2X,A,2X,F12.2,2X,ES12.4)')     '  2 ext1     ', h2a/1000, '  ---     ', T_sat_ext1, P_ext1
    write(*,'(A,F12.2,2X,A,2X,F12.2,2X,ES12.4)')     '  3 turb out ', h3a/1000, '  ---     ', T_sat_cond, P_cond
    write(*,'(A,F12.2,2X,A,2X,F12.2,2X,ES12.4)')     '  4 cond out ', h4/1000,  '  ---     ', T_sat_cond, P_cond
    write(*,'(A,F12.2,2X,A,2X,A,2X,ES12.4)')         '  6 boil in  ', h6a/1000, '  ---     ', '  ---       ', P_boiler
    write(*,*)
    write(*,'(A)') '--- EXTRACTION FRACTIONS ------------------------------------'
    write(*,'(A,F10.5)')    '  y1 (extraction 1)         = ', y1
    if (n_extractions == 2) &
    write(*,'(A,F10.5)')    '  y2 (extraction 2)         = ', y2
    write(*,'(A,F10.5)')    '  Condenser fraction        = ', 1.0d0-y1-y2
    write(*,*)
    write(*,'(A)') '--- ENERGY PER UNIT MASS ------------------------------------'
    write(*,'(A,F12.2,A)')  '  Turbine Work Wt           = ', Wt/1000, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Pump Work Wp              = ', Wp/1000, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Net Work Wnet             = ', Wnet/1000, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Heat Input Qin            = ', Qin/1000, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Heat Rejected Qout        = ', Qout/1000, ' kJ/kg'
    write(*,*)
    write(*,'(A)') '--- CYCLE PERFORMANCE ---------------------------------------'
    write(*,'(A,F10.4)')    '  Thermal Efficiency        = ', eta_th
    write(*,'(A,F10.2,A)')  '  Thermal Efficiency        = ', eta_th*100.0d0, ' percent'
    write(*,'(A,F10.4)')    '  Back-Work Ratio           = ', bwr
    write(*,'(A,F12.2,A)')  '  Power Output              = ', P_output, ' kW'
    write(*,'(A,F10.4,A)')  '  Specific Steam Consumption= ', SSC, ' kg/kWh'
    write(*,*)

    ! Efficiency vs extraction pressure sweep
    write(*,'(A)') '--- EFFICIENCY VS EXTRACTION PRESSURE SWEEP -----------------'
    write(*,'(A)') '  P_ext1[kPa]   y1            eta_th        Wnet[kJ/kg]'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, 40
        P_i = P_cond + (P_boiler - P_cond) * dble(i) / 41.0d0
        call compute_single_extraction(P_boiler, T_boiler, P_cond, P_i, &
            eta_turbine, eta_pump, h1, s1, h_f_cond, v_f_cond, eta_i, y_i)
        write(*,'(F12.2,2X,F10.5,2X,F10.5,2X,F12.2)') P_i/1000.0d0, y_i, eta_i, &
            eta_i * (h1 - h_f_cond - v_f_cond*(P_boiler-P_cond)/eta_pump) / 1000.0d0
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Simplified water saturation correlations (curve fits).'
    write(*,'(A)') '  Isentropic turbine: h_out = h_in - eta_t*(h_in - h_out_s).'
    write(*,'(A)') '  Pump work: w_p = v_f*(P_out - P_in)/eta_p.'
    write(*,'(A)') '  Open FWH: y*h_ext + (1-y)*h_fw = h_f_sat.'

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)
        ! Curve-fit approximations for water (valid ~5 kPa to 22 MPa)
        Tsat = 373.15d0 + 42.0d0*logP - 0.8d0*logP**2   ! K
        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  ! J/kg
        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   ! J/(kg K)
        if (sf < 0.3d3) sf = 0.3d3 + 0.5d3*P_MPa
        sfg = hfg / Tsat
        vf = 0.001d0 * (1.0d0 + 0.0002d0*P_MPa)   ! m3/kg
    end subroutine sat_props

    subroutine superheat_props(P_pa, T, Tsat, hf_boil, sf_boil, h_out, s_out)
        implicit none
        double precision, intent(in) :: P_pa, T, Tsat, hf_boil, sf_boil
        double precision, intent(out) :: h_out, s_out
        double precision :: P_MPa, hfg_b, sfg_b, hg, sg, cp_steam, dT
        P_MPa = P_pa / 1.0d6
        hfg_b = (2258.0d0 - 180.0d0*log(P_MPa) - 40.0d0*log(P_MPa)**2)*1000.0d0
        if (hfg_b < 200.0d3) hfg_b = 200.0d3
        hg = hf_boil + hfg_b
        sfg_b = hfg_b / Tsat
        sg = sf_boil + sfg_b
        cp_steam = 2100.0d0   ! J/(kg K) approximate
        dT = T - Tsat
        if (dT < 0.0d0) dT = 0.0d0
        h_out = hg + cp_steam * dT
        s_out = sg + cp_steam * log(max(T/Tsat, 1.0d0))
    end subroutine superheat_props

    subroutine superheat_from_entropy(P_pa, s_target, Tsat, hf, sf, sfg, h_out)
        implicit none
        double precision, intent(in) :: P_pa, s_target, Tsat, hf, sf, sfg
        double precision, intent(out) :: h_out
        double precision :: hfg, hg, sg, cp_steam, T_est, dT
        hfg = sfg * Tsat
        hg = hf + hfg
        sg = sf + sfg
        cp_steam = 2100.0d0
        ! s = sg + cp ln(T/Tsat) => T = Tsat * exp((s-sg)/cp)
        T_est = Tsat * exp((s_target - sg)/cp_steam)
        dT = T_est - Tsat
        if (dT < 0.0d0) dT = 0.0d0
        h_out = hg + cp_steam * dT
    end subroutine superheat_from_entropy

    subroutine compute_single_extraction(Pb, Tb, Pc, Pe, etat, etap, &
                                          h1in, s1in, hfc, vfc, eta_out, y_out)
        implicit none
        double precision, intent(in) :: Pb, Tb, Pc, Pe, etat, etap, h1in, s1in, hfc, vfc
        double precision, intent(out) :: eta_out, y_out
        double precision :: Tse, hfe, hfge, sfe, sfge, vfe, xe, h2si, h2ai
        double precision :: h3si, h3ai, x3i, h5i, h6i, Wti, Wpi, Qini
        double precision :: Tsc, hfc_local, hfgc, sfc, sfgc, vfc2

        call sat_props(Pe, Tse, hfe, hfge, sfe, sfge, vfe)
        call sat_props(Pc, Tsc, hfc_local, hfgc, sfc, sfgc, vfc2)

        xe = (s1in - sfe)/max(sfge,1.0d-10)
        if (xe > 1.0d0) then
            call superheat_from_entropy(Pe, s1in, Tse, hfe, sfe, sfge, h2si)
        else
            h2si = hfe + xe*hfge
        end if
        h2ai = h1in - etat*(h1in - h2si)

        x3i = (s1in - sfc)/max(sfgc,1.0d-10)
        if (x3i > 1.0d0) then
            call superheat_from_entropy(Pc, s1in, Tsc, hfc_local, sfc, sfgc, h3si)
        else
            h3si = hfc_local + x3i*hfgc
        end if
        h3ai = h1in - etat*(h1in - h3si)

        h5i = hfc_local + vfc2*(Pe-Pc)/etap
        y_out = (hfe - h5i)/max(h2ai - h5i, 1.0d-10)
        if (y_out < 0.0d0) y_out = 0.0d0
        if (y_out > 1.0d0) y_out = 1.0d0
        h6i = hfe + vfe*(Pb-Pe)/etap
        Wti = (h1in-h2ai)+(1.0d0-y_out)*(h2ai-h3ai)
        Wpi = (1.0d0-y_out)*(h5i-hfc_local)+(h6i-hfe)
        Qini = h1in - h6i
        eta_out = (Wti-Wpi)/max(Qini,1.0d-10)
    end subroutine compute_single_extraction

end program rankine_regenerative