program rankine_orc
    implicit none
    integer :: fluid_type, iostat_val, i, n_sweep, iter
    double precision :: T_source, T_sink, T_superheat, eta_t, eta_p, mdot, T_pinch
    double precision :: Tc, Pc, Tb, M_fluid, cp_liq, cp_vap, hfg_Tb, rho_liq
    double precision :: T_evap, T_cond, P_evap, P_cond
    double precision :: h1, s1, h2s, h2a, h3, h4, v_f
    double precision :: T2s, hfg_evap, hfg_cond
    double precision :: Wt, Wp, Wnet, Qin, Qout, eta_th, eta_Carnot, eta_rel
    double precision :: P_output, SSC, bwr
    double precision :: Ts_sw, eta_sw, Wnet_sw
    double precision :: Te_sw, Tc_sw, Pe_sw, Pcd_sw
    double precision :: h1s, s1s, h2ss, h2as, h3s, h4s, vfs
    double precision :: hfg_e_sw, hfg_c_sw, T2ss, Qins
    character(len=40) :: fluid_name

    read(*,*,iostat=iostat_val) fluid_type
    if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid fluid type.'; stop; end if
    read(*,*,iostat=iostat_val) T_source
    read(*,*,iostat=iostat_val) T_sink
    read(*,*,iostat=iostat_val) T_superheat
    read(*,*,iostat=iostat_val) eta_t
    read(*,*,iostat=iostat_val) eta_p
    read(*,*,iostat=iostat_val) mdot
    read(*,*,iostat=iostat_val) T_pinch
    if (iostat_val /= 0) then; write(*,*) 'ERROR: Failed to read all inputs.'; stop; end if

    if (T_source <= 0.0d0) T_source = 423.15d0
    if (T_sink <= 0.0d0) T_sink = 303.15d0
    if (T_superheat < 0.0d0) T_superheat = 5.0d0
    if (eta_t <= 0.0d0 .or. eta_t > 1.0d0) eta_t = 0.82d0
    if (eta_p <= 0.0d0 .or. eta_p > 1.0d0) eta_p = 0.75d0
    if (mdot <= 0.0d0) mdot = 10.0d0
    if (T_pinch <= 0.0d0) T_pinch = 10.0d0
    if (T_source <= T_sink + 2.0d0*T_pinch) then
        write(*,*) 'ERROR: T_source must exceed T_sink + 2*T_pinch.'
        stop
    end if

    ! ── Fluid properties ───────────────────────────────────────
    select case(fluid_type)
    case(1)
        fluid_name='R245fa (1,1,1,3,3-Pentafluoropropane)'
        Tc=427.2d0; Pc=3.651d0; Tb=288.3d0; M_fluid=134.05d0
        cp_liq=1360.0d0; cp_vap=900.0d0; hfg_Tb=196000.0d0; rho_liq=1320.0d0
    case(2)
        fluid_name='R134a (1,1,1,2-Tetrafluoroethane)'
        Tc=374.2d0; Pc=4.059d0; Tb=247.1d0; M_fluid=102.03d0
        cp_liq=1430.0d0; cp_vap=850.0d0; hfg_Tb=217000.0d0; rho_liq=1210.0d0
    case(3)
        fluid_name='n-Pentane (C5H12)'
        Tc=469.7d0; Pc=3.370d0; Tb=309.2d0; M_fluid=72.15d0
        cp_liq=2320.0d0; cp_vap=1670.0d0; hfg_Tb=358000.0d0; rho_liq=626.0d0
    case(4)
        fluid_name='Toluene (C7H8)'
        Tc=591.8d0; Pc=4.109d0; Tb=383.8d0; M_fluid=92.14d0
        cp_liq=1710.0d0; cp_vap=1130.0d0; hfg_Tb=363000.0d0; rho_liq=867.0d0
    case default
        fluid_name='R1233zd(E) (trans-1-Chloro-3,3,3-trifluoropropene)'
        Tc=438.8d0; Pc=3.624d0; Tb=291.4d0; M_fluid=130.50d0
        cp_liq=1250.0d0; cp_vap=820.0d0; hfg_Tb=195000.0d0; rho_liq=1290.0d0
        fluid_type=5
    end select

    ! ── Evaporator / condenser conditions ──────────────────────
    T_evap = T_source - T_pinch
    T_cond = T_sink + T_pinch
    if (T_evap >= Tc) T_evap = Tc - 5.0d0
    if (T_cond >= T_evap) then
        write(*,*) 'ERROR: T_evap must exceed T_cond.'
        stop
    end if

    ! Pressure from Clausius-Clapeyron approximation
    ! ln(P/Pc) ~ (hfg*M/(R*1000)) * (1/Tc - 1/T) simplified
    ! Use: P = Pc * exp( -5.0*(Tb/T - 1.0)*(Tc/Tb) )  (rough fit)
    P_evap = Pc * exp(-5.0d0 * (Tb/T_evap - 1.0d0) * Tc/Tb)
    if (P_evap > 0.90d0*Pc) P_evap = 0.90d0*Pc
    if (P_evap < 0.01d0) P_evap = 0.01d0
    P_cond = Pc * exp(-5.0d0 * (Tb/T_cond - 1.0d0) * Tc/Tb)
    if (P_cond > P_evap) P_cond = P_evap * 0.1d0
    if (P_cond < 0.001d0) P_cond = 0.001d0

    ! Latent heat adjusted with temperature (Trouton-Watson)
    hfg_evap = hfg_Tb * ((Tc - T_evap)/(Tc - Tb))**0.38d0
    hfg_cond = hfg_Tb * ((Tc - T_cond)/(Tc - Tb))**0.38d0
    if (hfg_evap < 0.0d0) hfg_evap = hfg_Tb*0.5d0
    if (hfg_cond < 0.0d0) hfg_cond = hfg_Tb*0.5d0

    ! ── State points (J/kg, J/(kg.K)) ─────────────────────────
    ! Reference: saturated liquid at Tb => h=0, s=0

    ! State 1: superheated vapor at P_evap
    h1 = cp_liq*(T_evap - Tb) + hfg_evap + cp_vap*T_superheat
    s1 = cp_liq*log(T_evap/Tb) + hfg_evap/T_evap + cp_vap*log((T_evap+T_superheat)/T_evap)

    ! State 2s: isentropic expansion to P_cond
    ! For dry ORC fluid, exit is superheated
    ! s2s = s1 => cp_liq*ln(T_cond/Tb) + hfg_cond/T_cond + cp_vap*ln(T2s/T_cond) = s1
    ! => T2s = T_cond * exp( (s1 - cp_liq*ln(T_cond/Tb) - hfg_cond/T_cond) / cp_vap )
    T2s = T_cond * exp((s1 - cp_liq*log(T_cond/Tb) - hfg_cond/T_cond) / cp_vap)
    if (T2s < T_cond) T2s = T_cond
    h2s = cp_liq*(T_cond - Tb) + hfg_cond + cp_vap*(T2s - T_cond)

    ! State 2a: actual turbine exit
    h2a = h1 - eta_t * (h1 - h2s)

    ! State 3: saturated liquid at condenser
    h3 = cp_liq * (T_cond - Tb)

    ! State 4: compressed liquid after pump
    v_f = 1.0d0 / rho_liq   ! m3/kg
    h4 = h3 + v_f * (P_evap - P_cond) * 1.0d6 / eta_p  ! Pa*m3/kg = J/kg

    ! ── Energy analysis (kJ/kg) ────────────────────────────────
    Wt = (h1 - h2a) / 1000.0d0
    Wp = (h4 - h3) / 1000.0d0
    Wnet = Wt - Wp
    Qin = (h1 - h4) / 1000.0d0
    Qout = (h2a - h3) / 1000.0d0
    eta_th = Wnet / max(Qin, 1.0d-10)
    eta_Carnot = 1.0d0 - T_sink / T_source
    eta_rel = eta_th / max(eta_Carnot, 1.0d-10)
    P_output = mdot * Wnet
    bwr = Wp / max(Wt, 1.0d-10)
    SSC = 3600.0d0 / max(Wnet, 1.0d-10)

    ! ── Output ──────────────────────────────────────────────────
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   ORGANIC RANKINE CYCLE (ORC)'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')        '  Working Fluid             = ', trim(fluid_name)
    write(*,'(A,F12.2,A)')  '  Heat Source Temperature   = ', T_source, ' K'
    write(*,'(A,F12.2,A)')  '  Heat Sink Temperature     = ', T_sink, ' K'
    write(*,'(A,F12.2,A)')  '  Superheat                 = ', T_superheat, ' K'
    write(*,'(A,F10.4)')    '  Turbine Isentropic Eff    = ', eta_t
    write(*,'(A,F10.4)')    '  Pump Isentropic Eff       = ', eta_p
    write(*,'(A,F12.4,A)')  '  Mass Flow Rate            = ', mdot, ' kg/s'
    write(*,'(A,F12.2,A)')  '  Pinch Point DeltaT        = ', T_pinch, ' K'
    write(*,*)
    write(*,'(A)') '--- FLUID PROPERTIES ----------------------------------------'
    write(*,'(A,F12.2,A)')  '  Critical Temperature Tc   = ', Tc, ' K'
    write(*,'(A,F12.4,A)')  '  Critical Pressure Pc      = ', Pc, ' MPa'
    write(*,'(A,F12.2,A)')  '  Normal Boiling Point Tb   = ', Tb, ' K'
    write(*,'(A,F12.2,A)')  '  Molar Mass M              = ', M_fluid, ' g/mol'
    write(*,'(A,F12.2,A)')  '  cp_liquid                 = ', cp_liq, ' J/(kg.K)'
    write(*,'(A,F12.2,A)')  '  cp_vapor                  = ', cp_vap, ' J/(kg.K)'
    write(*,'(A,F12.2,A)')  '  hfg at Tb                 = ', hfg_Tb/1000.0d0, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Liquid density            = ', rho_liq, ' kg/m3'
    write(*,*)
    write(*,'(A)') '--- CYCLE PRESSURES & TEMPERATURES --------------------------'
    write(*,'(A,F12.2,A)')  '  T_evaporator              = ', T_evap, ' K'
    write(*,'(A,F12.2,A)')  '  T_condenser               = ', T_cond, ' K'
    write(*,'(A,F12.6,A)')  '  P_evaporator              = ', P_evap, ' MPa'
    write(*,'(A,F12.6,A)')  '  P_condenser               = ', P_cond, ' MPa'
    write(*,'(A,F12.2,A)')  '  hfg at T_evap             = ', hfg_evap/1000.0d0, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  hfg at T_cond             = ', hfg_cond/1000.0d0, ' kJ/kg'
    write(*,*)
    write(*,'(A)') '--- STATE POINTS (kJ/kg, relative to sat liq at Tb) --------'
    write(*,'(A,F12.2,A)')  '  h1 (turb inlet)           = ', h1/1000.0d0, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  h2s (turb exit, isen)     = ', h2s/1000.0d0, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  h2a (turb exit, actual)   = ', h2a/1000.0d0, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  h3 (cond exit, sat liq)   = ', h3/1000.0d0, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  h4 (pump exit)            = ', h4/1000.0d0, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  T2s (turb exit temp)      = ', T2s, ' K'
    write(*,'(A,F12.2,A)')  '  T1 (turb inlet temp)      = ', T_evap+T_superheat, ' K'
    write(*,*)
    write(*,'(A)') '--- ENERGY PER UNIT MASS ------------------------------------'
    write(*,'(A,F12.4,A)')  '  Turbine Work Wt           = ', Wt, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  Pump Work Wp              = ', Wp, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  Net Work Wnet             = ', Wnet, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  Heat Input Qin            = ', Qin, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  Heat Rejected Qout        = ', Qout, ' kJ/kg'
    write(*,*)
    write(*,'(A)') '--- CYCLE PERFORMANCE ---------------------------------------'
    write(*,'(A,F10.6)')    '  Thermal Efficiency        = ', eta_th
    write(*,'(A,F10.2,A)')  '  Thermal Efficiency        = ', eta_th*100.0d0, ' percent'
    write(*,'(A,F10.6)')    '  Carnot Efficiency         = ', eta_Carnot
    write(*,'(A,F10.2,A)')  '  Relative Efficiency       = ', eta_rel*100.0d0, ' percent of Carnot'
    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 Consump    = ', SSC, ' kg/kWh'
    write(*,*)

    ! ── Sensitivity sweep: T_source ────────────────────────────
    n_sweep = 40
    write(*,'(A)') '--- SENSITIVITY: EFFICIENCY VS SOURCE TEMPERATURE ----------'
    write(*,'(A)') '  T_source[K]   eta_th        Wnet[kJ/kg]   P_out[kW]'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, n_sweep
        Ts_sw = (T_sink+30.0d0) + dble(i-1)*(T_source+80.0d0 - (T_sink+30.0d0))/dble(n_sweep-1)
        Te_sw = Ts_sw - T_pinch
        if (Te_sw >= Tc) Te_sw = Tc - 5.0d0
        Tc_sw = T_cond
        Pe_sw = Pc * exp(-5.0d0*(Tb/Te_sw - 1.0d0)*Tc/Tb)
        if (Pe_sw > 0.90d0*Pc) Pe_sw = 0.90d0*Pc
        if (Pe_sw < 0.01d0) Pe_sw = 0.01d0
        Pcd_sw = P_cond
        hfg_e_sw = hfg_Tb * ((Tc-Te_sw)/(Tc-Tb))**0.38d0
        if (hfg_e_sw < 0.0d0) hfg_e_sw = hfg_Tb*0.3d0
        hfg_c_sw = hfg_cond
        h1s = cp_liq*(Te_sw-Tb)+hfg_e_sw+cp_vap*T_superheat
        s1s = cp_liq*log(Te_sw/Tb)+hfg_e_sw/Te_sw+cp_vap*log((Te_sw+T_superheat)/Te_sw)
        T2ss = Tc_sw*exp((s1s-cp_liq*log(Tc_sw/Tb)-hfg_c_sw/Tc_sw)/cp_vap)
        if (T2ss < Tc_sw) T2ss = Tc_sw
        h2ss = cp_liq*(Tc_sw-Tb)+hfg_c_sw+cp_vap*(T2ss-Tc_sw)
        h2as = h1s - eta_t*(h1s-h2ss)
        h3s = cp_liq*(Tc_sw-Tb)
        vfs = 1.0d0/rho_liq
        h4s = h3s + vfs*(Pe_sw-Pcd_sw)*1.0d6/eta_p
        Qins = (h1s-h4s)/1000.0d0
        Wnet_sw = ((h1s-h2as)-(h4s-h3s))/1000.0d0
        if (Qins > 0.0d0) then
            eta_sw = Wnet_sw / Qins
        else
            eta_sw = 0.0d0
        end if
        write(*,'(F10.2,4X,F10.6,4X,F12.4,4X,F12.2)') Ts_sw, eta_sw, Wnet_sw, mdot*Wnet_sw
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Pressure: Clausius-Clapeyron approx with empirical fit.'
    write(*,'(A)') '  Latent heat: Watson correlation hfg ~ hfg_Tb*((Tc-T)/(Tc-Tb))^0.38.'
    write(*,'(A)') '  Isentropic expansion assuming dry fluid (superheated exit).'
    write(*,'(A)') '  Pump work: w_p = v_f*(P_evap-P_cond)/eta_p.'
    write(*,'(A)') '  Simplified constant cp model (educational accuracy).'

end program rankine_orc
