program carnot_heatpump
    implicit none
    integer :: cycle_type, iostat_val, i, n_sweep
    double precision :: T_hot, T_cold, Q_desired, real_frac
    double precision :: eta_carnot, COP_ref, COP_hp
    double precision :: W_cycle, Q_hot, Q_cold
    double precision :: COP_ref_real, COP_hp_real, eta_real
    double precision :: W_real, Q_hot_real, Q_cold_real
    double precision :: T_sw, eta_sw, COP_ref_sw, COP_hp_sw
    character(len=40) :: cycle_name

    ! ── Read inputs ─────────────────────────────────────────────
    read(*,*,iostat=iostat_val) T_hot
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid T_hot input.'
        stop
    end if
    read(*,*,iostat=iostat_val) T_cold
    read(*,*,iostat=iostat_val) Q_desired
    read(*,*,iostat=iostat_val) cycle_type
    read(*,*,iostat=iostat_val) real_frac
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all inputs.'
        stop
    end if

    ! ── Input validation ────────────────────────────────────────
    if (T_hot <= 0.0d0) then
        write(*,*) 'ERROR: T_hot must be positive (K).'
        stop
    end if
    if (T_cold <= 0.0d0) then
        write(*,*) 'ERROR: T_cold must be positive (K).'
        stop
    end if
    if (T_hot <= T_cold) then
        write(*,*) 'ERROR: T_hot must exceed T_cold.'
        stop
    end if
    if (Q_desired <= 0.0d0) Q_desired = 100.0d0
    if (cycle_type < 1 .or. cycle_type > 3) cycle_type = 1
    if (real_frac <= 0.0d0 .or. real_frac > 1.0d0) real_frac = 0.50d0

    if (cycle_type == 1) then
        cycle_name = 'Carnot Heat Engine'
    else if (cycle_type == 2) then
        cycle_name = 'Refrigerator (Reversed Carnot)'
    else
        cycle_name = 'Heat Pump (Reversed Carnot)'
    end if

    ! ── Carnot limits ───────────────────────────────────────────
    eta_carnot = 1.0d0 - T_cold / T_hot
    COP_ref = T_cold / (T_hot - T_cold)
    COP_hp  = T_hot  / (T_hot - T_cold)

    ! ── Energy balance based on cycle type ──────────────────────
    if (cycle_type == 1) then
        ! Heat engine: Q_desired = Q_hot (heat input)
        Q_hot = Q_desired
        W_cycle = Q_hot * eta_carnot
        Q_cold = Q_hot - W_cycle
    else if (cycle_type == 2) then
        ! Refrigerator: Q_desired = Q_cold (cooling load)
        Q_cold = Q_desired
        W_cycle = Q_cold / max(COP_ref, 1.0d-10)
        Q_hot = Q_cold + W_cycle
    else
        ! Heat pump: Q_desired = Q_hot (heating load)
        Q_hot = Q_desired
        W_cycle = Q_hot / max(COP_hp, 1.0d-10)
        Q_cold = Q_hot - W_cycle
    end if

    ! ── Real cycle estimates ────────────────────────────────────
    eta_real = eta_carnot * real_frac
    COP_ref_real = COP_ref * real_frac
    COP_hp_real  = COP_hp  * real_frac

    if (cycle_type == 1) then
        W_real = Q_desired * eta_real
        Q_hot_real = Q_desired
        Q_cold_real = Q_hot_real - W_real
    else if (cycle_type == 2) then
        W_real = Q_desired / max(COP_ref_real, 1.0d-10)
        Q_cold_real = Q_desired
        Q_hot_real = Q_cold_real + W_real
    else
        W_real = Q_desired / max(COP_hp_real, 1.0d-10)
        Q_hot_real = Q_desired
        Q_cold_real = Q_hot_real - W_real
    end if

    ! ── Output ──────────────────────────────────────────────────
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   CARNOT, REVERSED CARNOT & HEAT PUMP'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,F12.2,A)')  '  Hot Reservoir T_hot       = ', T_hot, ' K'
    write(*,'(A,F12.2,A)')  '  Hot Reservoir T_hot       = ', T_hot-273.15d0, ' C'
    write(*,'(A,F12.2,A)')  '  Cold Reservoir T_cold     = ', T_cold, ' K'
    write(*,'(A,F12.2,A)')  '  Cold Reservoir T_cold     = ', T_cold-273.15d0, ' C'
    write(*,'(A,F12.2,A)')  '  Q_desired                 = ', Q_desired, ' kW'
    write(*,'(A,A)')        '  Cycle Type                = ', trim(cycle_name)
    write(*,'(A,F10.4)')    '  Real Cycle Factor         = ', real_frac
    write(*,*)
    write(*,'(A)') '--- CARNOT LIMITS -------------------------------------------'
    write(*,'(A,F10.6)')    '  Carnot Efficiency         = ', eta_carnot
    write(*,'(A,F10.2,A)')  '  Carnot Efficiency         = ', eta_carnot*100.0d0, ' percent'
    write(*,'(A,F10.4)')    '  COP Refrigeration (max)   = ', COP_ref
    write(*,'(A,F10.4)')    '  COP Heat Pump (max)       = ', COP_hp
    write(*,'(A)')          '  COP_hp = COP_ref + 1      (verified)'
    write(*,*)
    write(*,'(A)') '--- ENERGY BALANCE (IDEAL) ----------------------------------'
    write(*,'(A,F12.4,A)')  '  Q_hot                     = ', Q_hot, ' kW'
    write(*,'(A,F12.4,A)')  '  Q_cold                    = ', Q_cold, ' kW'
    write(*,'(A,F12.4,A)')  '  W_cycle                   = ', W_cycle, ' kW'
    write(*,*)
    write(*,'(A)') '--- REAL CYCLE ESTIMATES ------------------------------------'
    write(*,'(A,F10.4)')    '  Real Efficiency           = ', eta_real
    write(*,'(A,F10.2,A)')  '  Real Efficiency           = ', eta_real*100.0d0, ' percent'
    write(*,'(A,F10.4)')    '  Real COP Refrigeration    = ', COP_ref_real
    write(*,'(A,F10.4)')    '  Real COP Heat Pump        = ', COP_hp_real
    write(*,'(A,F12.4,A)')  '  Real W_cycle              = ', W_real, ' kW'
    write(*,'(A,F12.4,A)')  '  Real Q_hot                = ', Q_hot_real, ' kW'
    write(*,'(A,F12.4,A)')  '  Real Q_cold               = ', Q_cold_real, ' kW'
    write(*,*)

    ! ── Sensitivity sweep: T_cold from 200 K to T_hot-5 ────────
    n_sweep = 40
    write(*,'(A)') '--- SENSITIVITY: PERFORMANCE VS T_COLD ----------------------'
    write(*,'(A)') '  T_cold[K]     eta_Carnot    COP_ref       COP_hp'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, n_sweep
        T_sw = 200.0d0 + dble(i-1) * (T_hot - 5.0d0 - 200.0d0) / dble(n_sweep - 1)
        if (T_sw >= T_hot) T_sw = T_hot - 1.0d0
        eta_sw = 1.0d0 - T_sw / T_hot
        COP_ref_sw = T_sw / max(T_hot - T_sw, 1.0d-10)
        COP_hp_sw  = T_hot / max(T_hot - T_sw, 1.0d-10)
        write(*,'(F10.2,4X,F10.6,4X,F10.4,4X,F10.4)') T_sw, eta_sw, COP_ref_sw, COP_hp_sw
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Carnot efficiency: eta = 1 - T_cold/T_hot'
    write(*,'(A)') '  COP_refrigeration = T_cold / (T_hot - T_cold)'
    write(*,'(A)') '  COP_heat_pump = T_hot / (T_hot - T_cold)'
    write(*,'(A)') '  COP_hp = COP_ref + 1'
    write(*,'(A)') '  Real cycle: multiply ideal COP by factor (0.4-0.6 typical)'
    write(*,'(A)') '  Energy balance: Q_hot = Q_cold + W'

end program carnot_heatpump
