program entropy_exergy
    implicit none
    integer :: i, iostat_val, n_sweep
    double precision :: T_sys, T0, P_sys, P0, mdot, cp_val
    double precision :: h_in, h_out, s_in, s_out, h0, s0
    double precision :: Q_dot_in, Q_dot_out, T_source, T_sink
    double precision :: S_gen, X_dest, psi_in, psi_out
    double precision :: X_in, X_out, eta_ex, X_Q_in, X_Q_out
    double precision :: dH, dS, W_rev, W_act, eta_2nd
    double precision :: T0_sweep, S_gen_sw, X_dest_sw, eta_sw

    ! ── Read inputs ─────────────────────────────────────────────
    read(*,*,iostat=iostat_val) T_sys
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid system temperature input.'
        stop
    end if
    read(*,*,iostat=iostat_val) T0
    read(*,*,iostat=iostat_val) P_sys
    read(*,*,iostat=iostat_val) P0
    read(*,*,iostat=iostat_val) mdot
    read(*,*,iostat=iostat_val) cp_val
    read(*,*,iostat=iostat_val) h_in
    read(*,*,iostat=iostat_val) h_out
    read(*,*,iostat=iostat_val) s_in
    read(*,*,iostat=iostat_val) s_out
    read(*,*,iostat=iostat_val) h0
    read(*,*,iostat=iostat_val) s0
    read(*,*,iostat=iostat_val) Q_dot_in
    read(*,*,iostat=iostat_val) Q_dot_out
    read(*,*,iostat=iostat_val) T_source
    read(*,*,iostat=iostat_val) T_sink
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all inputs.'
        stop
    end if

    ! ── Input validation ────────────────────────────────────────
    if (T0 <= 0.0d0) T0 = 298.15d0
    if (T_source <= 0.0d0) T_source = 1000.0d0
    if (T_sink <= 0.0d0) T_sink = T0
    if (mdot <= 0.0d0) mdot = 1.0d0
    if (cp_val <= 0.0d0) cp_val = 1.005d0

    ! ── Core calculations ───────────────────────────────────────
    ! Entropy generation (general entropy balance for open system)
    S_gen = mdot * (s_out - s_in)
    if (T_source > 0.0d0 .and. Q_dot_in > 0.0d0) then
        S_gen = S_gen - Q_dot_in / T_source
    end if
    if (T_sink > 0.0d0 .and. Q_dot_out > 0.0d0) then
        S_gen = S_gen + Q_dot_out / T_sink
    end if
    if (S_gen < 0.0d0) S_gen = 0.0d0   ! cannot be negative (2nd law)

    ! Exergy destruction (Gouy-Stodola theorem)
    X_dest = T0 * S_gen

    ! Flow exergy (specific)
    psi_in  = (h_in - h0) - T0 * (s_in - s0)
    psi_out = (h_out - h0) - T0 * (s_out - s0)

    ! Exergy rates
    X_in  = mdot * psi_in
    X_out = mdot * psi_out

    ! Heat exergy (Carnot factor)
    X_Q_in = 0.0d0
    if (T_source > 0.0d0) X_Q_in = Q_dot_in * (1.0d0 - T0 / T_source)
    X_Q_out = 0.0d0
    if (T_sink > 0.0d0) X_Q_out = Q_dot_out * (1.0d0 - T0 / T_sink)

    ! Exergy efficiency
    if (abs(X_in) > 1.0d-10) then
        eta_ex = X_out / X_in
    else
        eta_ex = 0.0d0
    end if

    ! Enthalpy and entropy changes
    dH = mdot * (h_out - h_in)
    dS = mdot * (s_out - s_in)

    ! Reversible work and second-law efficiency
    W_rev = X_in - X_out + X_Q_in - X_Q_out
    W_act = dH
    if (abs(W_rev) > 1.0d-10) then
        eta_2nd = W_act / W_rev
    else
        eta_2nd = 0.0d0
    end if

    ! ── Output ──────────────────────────────────────────────────
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   ENTROPY & EXERGY ANALYSIS'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,F12.2,A)')  '  System Temperature        = ', T_sys, ' K'
    write(*,'(A,F12.2,A)')  '  Dead State Temperature T0 = ', T0, ' K'
    write(*,'(A,F12.2,A)')  '  System Pressure           = ', P_sys, ' kPa'
    write(*,'(A,F12.2,A)')  '  Dead State Pressure P0    = ', P0, ' kPa'
    write(*,'(A,F12.4,A)')  '  Mass Flow Rate            = ', mdot, ' kg/s'
    write(*,'(A,F12.4,A)')  '  cp                        = ', cp_val, ' kJ/(kg.K)'
    write(*,'(A,F12.2,A)')  '  h_in                      = ', h_in, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  h_out                     = ', h_out, ' kJ/kg'
    write(*,'(A,F12.6,A)')  '  s_in                      = ', s_in, ' kJ/(kg.K)'
    write(*,'(A,F12.6,A)')  '  s_out                     = ', s_out, ' kJ/(kg.K)'
    write(*,'(A,F12.2,A)')  '  h0 (dead state)           = ', h0, ' kJ/kg'
    write(*,'(A,F12.6,A)')  '  s0 (dead state)           = ', s0, ' kJ/(kg.K)'
    write(*,'(A,F12.4,A)')  '  Q_dot_in                  = ', Q_dot_in, ' kW'
    write(*,'(A,F12.4,A)')  '  Q_dot_out                 = ', Q_dot_out, ' kW'
    write(*,'(A,F12.2,A)')  '  T_source                  = ', T_source, ' K'
    write(*,'(A,F12.2,A)')  '  T_sink                    = ', T_sink, ' K'
    write(*,*)
    write(*,'(A)') '--- ENTROPY RESULTS -----------------------------------------'
    write(*,'(A,ES14.6,A)')  '  Entropy Generation Sgen   = ', S_gen, ' kW/K'
    write(*,'(A,F12.4,A)')   '  Entropy Change dS         = ', dS, ' kW/K'
    write(*,*)
    write(*,'(A)') '--- EXERGY RESULTS ------------------------------------------'
    write(*,'(A,F12.4,A)')   '  Exergy Destruction Xdest  = ', X_dest, ' kW'
    write(*,'(A,F12.4,A)')   '  Flow Exergy In psi_in     = ', psi_in, ' kJ/kg'
    write(*,'(A,F12.4,A)')   '  Flow Exergy Out psi_out   = ', psi_out, ' kJ/kg'
    write(*,'(A,F12.4,A)')   '  Exergy Rate In X_in       = ', X_in, ' kW'
    write(*,'(A,F12.4,A)')   '  Exergy Rate Out X_out     = ', X_out, ' kW'
    write(*,'(A,F12.4,A)')   '  Heat Exergy In X_Q_in     = ', X_Q_in, ' kW'
    write(*,'(A,F12.4,A)')   '  Heat Exergy Out X_Q_out   = ', X_Q_out, ' kW'
    write(*,'(A,F12.4,A)')   '  Irreversibility Rate      = ', X_dest, ' kW'
    write(*,*)
    write(*,'(A)') '--- EFFICIENCY ----------------------------------------------'
    write(*,'(A,F10.6)')     '  Exergy Efficiency eta_ex  = ', eta_ex
    write(*,'(A,F10.2,A)')   '  Exergy Efficiency         = ', eta_ex*100.0d0, ' percent'
    write(*,'(A,F12.4,A)')   '  Reversible Work W_rev     = ', W_rev, ' kW'
    write(*,'(A,F12.4,A)')   '  Actual Energy Change dH   = ', dH, ' kW'
    write(*,'(A,F10.6)')     '  Second Law Efficiency     = ', eta_2nd
    write(*,'(A,F10.2,A)')   '  Second Law Efficiency     = ', eta_2nd*100.0d0, ' percent'
    write(*,*)

    ! ── Sensitivity sweep: T0 from 273 to 323 K ────────────────
    n_sweep = 50
    write(*,'(A)') '--- SENSITIVITY: EXERGY VS DEAD-STATE TEMPERATURE ----------'
    write(*,'(A)') '  T0[K]         Sgen[kW/K]      Xdest[kW]       eta_ex'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 0, n_sweep
        T0_sweep = 273.0d0 + dble(i) * (323.0d0 - 273.0d0) / dble(n_sweep)
        S_gen_sw = mdot * (s_out - s_in)
        if (T_source > 0.0d0 .and. Q_dot_in > 0.0d0) then
            S_gen_sw = S_gen_sw - Q_dot_in / T_source
        end if
        if (T_sink > 0.0d0 .and. Q_dot_out > 0.0d0) then
            S_gen_sw = S_gen_sw + Q_dot_out / T_sink
        end if
        if (S_gen_sw < 0.0d0) S_gen_sw = 0.0d0
        X_dest_sw = T0_sweep * S_gen_sw
        ! recompute exergy with swept T0
        psi_in  = (h_in - h0) - T0_sweep * (s_in - s0)
        psi_out = (h_out - h0) - T0_sweep * (s_out - s0)
        if (abs(mdot * psi_in) > 1.0d-10) then
            eta_sw = (mdot * psi_out) / (mdot * psi_in)
        else
            eta_sw = 0.0d0
        end if
        write(*,'(F10.2,4X,ES14.6,2X,F12.4,4X,F10.6)') &
            T0_sweep, S_gen_sw, X_dest_sw, eta_sw
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Gouy-Stodola theorem: X_dest = T0 * S_gen'
    write(*,'(A)') '  Flow exergy: psi = (h-h0) - T0*(s-s0)'
    write(*,'(A)') '  Heat exergy: X_Q = Q*(1 - T0/T_boundary)'
    write(*,'(A)') '  Exergy efficiency: eta_ex = X_out / X_in'

end program entropy_exergy