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