﻿program rankine_cycle
    implicit none
    
    double precision :: P_boiler, P_cond, T_super, eta_pump, eta_turb
    double precision :: T_sat_h, T_sat_l, h_fg_h, h_fg_l
    double precision :: h1, h2, h2s, h3, h4, h4s
    double precision :: s1, s3, s4s
    double precision :: W_pump, W_turb, Q_in, Q_out, W_net
    double precision :: eta_thermal, BWR, SSC
    double precision :: m_dot, W_net_total, Q_in_total
    double precision :: v_f1  ! specific volume at condenser exit
    integer :: mode  ! 1=simple, 2=reheat
    
    ! Reheat variables
    double precision :: P_reheat, T_reheat
    double precision :: h5, h6, h6s, W_turb2
    
    ! Simplified steam properties (polynomial approximations)
    double precision :: T_sat, h_f, h_g, s_f, s_g, v_f
    
    ! Read inputs
    read(*,*) mode         ! 1=simple Rankine, 2=reheat
    read(*,*) P_boiler     ! Boiler pressure [kPa]
    read(*,*) T_super      ! Superheat temperature [deg-C]
    read(*,*) P_cond       ! Condenser pressure [kPa]
    read(*,*) eta_turb     ! Turbine isentropic efficiency [0-1]
    read(*,*) eta_pump     ! Pump isentropic efficiency [0-1]
    read(*,*) m_dot        ! Mass flow rate [kg/s]
    
    if (mode == 2) then
        read(*,*) P_reheat   ! Reheat pressure [kPa]
        read(*,*) T_reheat   ! Reheat temperature [deg-C]
    end if
    
    ! ============================================
    ! SIMPLIFIED STEAM PROPERTY APPROXIMATIONS
    ! (reasonable for educational/preliminary design)
    ! ============================================
    
    ! Saturation temperature approximation (Antoine-like)
    ! T_sat ≈ f(P) for water
    T_sat_h = 99.974d0 + 28.07d0 * log(P_boiler/101.325d0)  ! Rough approx
    T_sat_l = 99.974d0 + 28.07d0 * log(P_cond/101.325d0)
    
    ! Better approximation for common pressures
    if (P_cond < 20) T_sat_l = 17.5d0 + 13.0d0 * log(P_cond)
    if (P_cond < 5) T_sat_l = 32.88d0   ! ~5 kPa
    if (P_cond < 10) T_sat_l = 45.81d0  ! ~10 kPa
    
    ! State 1: Condenser exit (saturated liquid at P_cond)
    h1 = 4.18d0 * T_sat_l   ! Approximate h_f
    v_f1 = 0.001d0           ! Approximate v_f for liquid water [m3/kg]
    s1 = 4.18d0 * log((T_sat_l + 273.15d0)/273.15d0)  ! Approximate s_f
    
    ! State 2: Pump exit (compressed liquid at P_boiler)
    W_pump = v_f1 * (P_boiler - P_cond)  ! Ideal pump work [kJ/kg] (already in kJ since P in kPa)
    h2s = h1 + W_pump
    h2 = h1 + W_pump / eta_pump
    
    ! State 3: Turbine inlet (superheated steam at P_boiler, T_super)
    ! Approximate superheated steam enthalpy
    ! h ≈ h_g(P) + Cp_steam x (T - T_sat)
    ! h_g ≈ 2675 + 1.0 x (P_boiler/1000 - 0.1) for mid pressures
    h_fg_h = 2257.0d0 - 2.0d0 * (T_sat_h - 100.0d0)  ! Rough h_fg
    h3 = 4.18d0 * T_sat_h + h_fg_h + 2.1d0 * (T_super - T_sat_h)
    
    ! More reasonable approximation
    if (P_boiler >= 500 .and. P_boiler <= 20000) then
        h3 = 2800.0d0 + 2.1d0 * (T_super - T_sat_h)
        if (T_super > 500) h3 = h3 + 0.3d0 * (T_super - 500)
    end if
    
    s3 = 6.5d0 + 0.002d0 * (T_super - 300.0d0)  ! Approximate entropy
    
    ! State 4: Turbine exit (mixture at P_cond)
    ! Isentropic expansion: s4s = s3
    s4s = s3
    ! Quality at state 4s
    h_fg_l = 2392.0d0  ! Approximate h_fg at low pressure
    h4s = h1 + (s4s - s1) / (7.36d0 - s1) * h_fg_l  ! Rough mixture enthalpy
    
    ! More reasonable: assume exit enthalpy based on expansion
    h4s = h3 - (h3 - h1) * 0.85d0  ! Rough isentropic estimate
    
    ! Actual turbine exit
    h4 = h3 - eta_turb * (h3 - h4s)
    
    ! ============================================
    ! CYCLE PERFORMANCE
    ! ============================================
    W_turb = h3 - h4       ! Turbine work [kJ/kg]
    Q_in = h3 - h2         ! Heat input [kJ/kg]
    Q_out = h4 - h1        ! Heat rejected [kJ/kg]
    W_net = W_turb - W_pump / eta_pump  ! Net work [kJ/kg]
    
    if (mode == 2) then
        ! Reheat: second expansion
        h5 = 2800.0d0 + 2.1d0 * (T_reheat - T_sat_h * P_reheat / P_boiler)
        h6s = h5 - (h5 - h1) * 0.8d0
        h6 = h5 - eta_turb * (h5 - h6s)
        
        W_turb2 = h5 - h6
        W_turb = (h3 - h4) + W_turb2
        Q_in = (h3 - h2) + (h5 - h4)  ! Two heat additions
        Q_out = h6 - h1
        W_net = W_turb - W_pump / eta_pump
        h4 = h6  ! For display purposes, state 4 becomes turbine 2 exit
    end if
    
    eta_thermal = W_net / Q_in
    BWR = (W_pump / eta_pump) / W_turb  ! Back work ratio
    SSC = 3600.0d0 / W_net  ! Specific steam consumption [kg/kWh]
    
    ! Total power
    W_net_total = m_dot * W_net  ! kW
    Q_in_total = m_dot * Q_in   ! kW
    
    ! ============================================
    ! OUTPUT
    ! ============================================
    write(*,'(A)') '============================================================'
    if (mode == 1) then
        write(*,'(A)') '   SIMPLE RANKINE CYCLE ANALYSIS'
    else
        write(*,'(A)') '   RANKINE CYCLE WITH REHEAT'
    end if
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- OPERATING CONDITIONS ---'
    write(*,'(A,F12.2,A)')  '  Boiler Pressure         = ', P_boiler, ' kPa'
    write(*,'(A,F12.2,A)')  '  Superheat Temperature   = ', T_super, ' deg-C'
    write(*,'(A,F12.2,A)')  '  Condenser Pressure      = ', P_cond, ' kPa'
    write(*,'(A,F12.4)')    '  Turbine Efficiency      = ', eta_turb
    write(*,'(A,F12.4)')    '  Pump Efficiency         = ', eta_pump
    write(*,'(A,F12.4,A)')  '  Mass Flow Rate          = ', m_dot, ' kg/s'
    if (mode == 2) then
        write(*,'(A,F12.2,A)')  '  Reheat Pressure         = ', P_reheat, ' kPa'
        write(*,'(A,F12.2,A)')  '  Reheat Temperature      = ', T_reheat, ' deg-C'
    end if
    write(*,*)
    write(*,'(A)') '--- STATE POINTS ---'
    write(*,'(A)') '  State   Description              h [kJ/kg]    P [kPa]'
    write(*,'(A)') '  ---'
    write(*,'(A,F12.2,4X,F10.2)') '  1       Condenser exit (sat liq)  ', h1, P_cond
    write(*,'(A,F12.2,4X,F10.2)') '  2       Pump exit (comp liquid)   ', h2, P_boiler
    write(*,'(A,F12.2,4X,F10.2)') '  3       Turbine inlet (superheat) ', h3, P_boiler
    if (mode == 1) then
        write(*,'(A,F12.2,4X,F10.2)') '  4       Turbine exit (mixture)   ', h4, P_cond
    else
        write(*,'(A,F12.2,4X,F10.2)') '  4       HP Turbine exit          ', h3-(h3-h4s)*eta_turb, P_reheat
        write(*,'(A,F12.2,4X,F10.2)') '  5       Reheat exit              ', h5, P_reheat
        write(*,'(A,F12.2,4X,F10.2)') '  6       LP Turbine exit          ', h4, P_cond
    end if
    write(*,*)
    write(*,'(A)') '--- ENERGY ANALYSIS (per kg) ---'
    write(*,'(A,F12.2,A)')  '  Turbine Work (W_t)      = ', W_turb, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  Pump Work (W_p)         = ', W_pump/eta_pump, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Net Work (W_net)        = ', W_net, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Heat Input (Q_in)       = ', Q_in, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Heat Rejected (Q_out)   = ', Q_out, ' kJ/kg'
    write(*,*)
    write(*,'(A)') '--- PERFORMANCE ---'
    write(*,'(A,F12.2,A)')  '  Thermal Efficiency (eta)  = ', eta_thermal*100, ' %'
    write(*,'(A,F12.4,A)')  '  Back Work Ratio (BWR)   = ', BWR*100, ' %'
    write(*,'(A,F12.4,A)')  '  Spec. Steam Consumption = ', SSC, ' kg/kWh'
    write(*,*)
    write(*,'(A)') '--- TOTAL POWER OUTPUT ---'
    write(*,'(A,F12.2,A)')  '  Net Power Output        = ', W_net_total, ' kW'
    write(*,'(A,F12.4,A)')  '  Net Power Output        = ', W_net_total/1000, ' MW'
    write(*,'(A,F12.2,A)')  '  Heat Input Rate         = ', Q_in_total, ' kW'
    write(*,'(A,F12.4,A)')  '  Heat Input Rate         = ', Q_in_total/1000, ' MW'
    write(*,*)
    
    ! Carnot comparison
    write(*,'(A)') '--- CARNOT COMPARISON ---'
    write(*,'(A,F12.2,A)')  '  T_hot (superheat)       = ', T_super + 273.15d0, ' K'
    write(*,'(A,F12.2,A)')  '  T_cold (condenser)      = ', T_sat_l + 273.15d0, ' K'
    write(*,'(A,F12.2,A)')  '  Carnot Efficiency       = ', &
        (1.0d0 - (T_sat_l+273.15d0)/(T_super+273.15d0))*100, ' %'
    write(*,'(A,F12.2,A)')  '  2nd Law Efficiency      = ', &
        eta_thermal / (1.0d0 - (T_sat_l+273.15d0)/(T_super+273.15d0)) * 100, ' %'
    write(*,*)
    write(*,'(A)') '--- NOTES ---'
    write(*,'(A)') '  * Steam properties are approximated for educational use'
    write(*,'(A)') '  * For precise design, use IAPWS-IF97 steam tables'
    write(*,'(A)') '  * Results are indicative of cycle performance trends'
    
end program rankine_cycle
