﻿program refrigeration_cycle
    implicit none
    
    double precision :: T_evap, T_cond, eta_comp
    double precision :: h1, h2, h2s, h3, h4
    double precision :: P_evap, P_cond_p
    double precision :: W_comp, Q_L, Q_H, COP_ref, COP_hp
    double precision :: m_dot, Q_L_total, W_comp_total
    double precision :: COP_carnot_ref, COP_carnot_hp
    double precision :: T_evap_K, T_cond_K
    integer :: refrig  ! 1=R134a, 2=R410A, 3=R22, 4=R717(ammonia)
    character(len=30) :: refrig_name
    
    ! Read inputs
    read(*,*) refrig
    read(*,*) T_evap      ! Evaporator temperature [deg-C]
    read(*,*) T_cond      ! Condenser temperature [deg-C]
    read(*,*) eta_comp    ! Compressor isentropic efficiency [0-1]
    read(*,*) m_dot       ! Mass flow rate [kg/s]
    
    T_evap_K = T_evap + 273.15d0
    T_cond_K = T_cond + 273.15d0
    
    ! ============================================
    ! REFRIGERANT PROPERTIES (simplified approx)
    ! ============================================
    select case(refrig)
    case(1)
        refrig_name = 'R-134a'
        ! Approximate saturation properties
        ! h_g at evaporator (sat vapor) ≈ 400 + 0.4*T
        h1 = 398.0d0 + 0.5d0 * T_evap       ! Sat vapor at T_evap [kJ/kg]
        h3 = 255.0d0 - 1.2d0 * (T_cond - 40)  ! Sat liquid at T_cond [kJ/kg]
        h2s = h1 + (42.0d0 + 0.8d0*(T_cond - T_evap))  ! Approx isentropic compression
        P_evap = 292.0d0 + 8.0d0 * T_evap    ! Rough P_sat [kPa]
        P_cond_p = 292.0d0 + 8.0d0 * T_cond
        
    case(2)
        refrig_name = 'R-410A'
        h1 = 420.0d0 + 0.4d0 * T_evap
        h3 = 260.0d0 - 1.3d0 * (T_cond - 40)
        h2s = h1 + (48.0d0 + 0.9d0*(T_cond - T_evap))
        P_evap = 800.0d0 + 20.0d0 * T_evap
        P_cond_p = 800.0d0 + 20.0d0 * T_cond
        
    case(3)
        refrig_name = 'R-22'
        h1 = 405.0d0 + 0.45d0 * T_evap
        h3 = 250.0d0 - 1.1d0 * (T_cond - 40)
        h2s = h1 + (40.0d0 + 0.75d0*(T_cond - T_evap))
        P_evap = 498.0d0 + 12.0d0 * T_evap
        P_cond_p = 498.0d0 + 12.0d0 * T_cond
        
    case(4)
        refrig_name = 'R-717 (Ammonia)'
        h1 = 1450.0d0 + 1.5d0 * T_evap
        h3 = 370.0d0 - 3.5d0 * (T_cond - 40)
        h2s = h1 + (180.0d0 + 3.0d0*(T_cond - T_evap))
        P_evap = 430.0d0 + 15.0d0 * T_evap
        P_cond_p = 430.0d0 + 15.0d0 * T_cond
        
    case default
        refrig_name = 'R-134a'
        h1 = 398.0d0 + 0.5d0 * T_evap
        h3 = 255.0d0 - 1.2d0 * (T_cond - 40)
        h2s = h1 + (42.0d0 + 0.8d0*(T_cond - T_evap))
        P_evap = 292.0d0 + 8.0d0 * T_evap
        P_cond_p = 292.0d0 + 8.0d0 * T_cond
    end select
    
    ! State 2: Compressor exit (actual)
    h2 = h1 + (h2s - h1) / eta_comp
    
    ! State 4: Throttling valve exit (h4 = h3, isenthalpic)
    h4 = h3
    
    ! ============================================
    ! CYCLE CALCULATIONS
    ! ============================================
    W_comp = h2 - h1          ! Compressor work [kJ/kg]
    Q_L = h1 - h4             ! Refrigeration effect [kJ/kg]
    Q_H = h2 - h3             ! Heat rejected [kJ/kg]
    
    COP_ref = Q_L / W_comp    ! COP (refrigeration)
    COP_hp = Q_H / W_comp     ! COP (heat pump)
    
    ! Carnot COPs
    COP_carnot_ref = T_evap_K / (T_cond_K - T_evap_K)
    COP_carnot_hp = T_cond_K / (T_cond_K - T_evap_K)
    
    ! Total capacities
    Q_L_total = m_dot * Q_L    ! kW
    W_comp_total = m_dot * W_comp  ! kW
    
    ! ============================================
    ! OUTPUT
    ! ============================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   VAPOR COMPRESSION REFRIGERATION CYCLE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- OPERATING CONDITIONS ---'
    write(*,'(A,A)')        '  Refrigerant             = ', trim(refrig_name)
    write(*,'(A,F12.2,A)')  '  Evaporator Temp (T_L)   = ', T_evap, ' deg-C'
    write(*,'(A,F12.2,A)')  '  Condenser Temp (T_H)    = ', T_cond, ' deg-C'
    write(*,'(A,F12.2,A)')  '  Evaporator Pressure     ≈ ', P_evap, ' kPa'
    write(*,'(A,F12.2,A)')  '  Condenser Pressure      ≈ ', P_cond_p, ' kPa'
    write(*,'(A,F12.4)')    '  Compressor Efficiency   = ', eta_comp
    write(*,'(A,F12.4,A)')  '  Mass Flow Rate          = ', m_dot, ' kg/s'
    write(*,*)
    write(*,'(A)') '--- STATE POINTS ---'
    write(*,'(A)') '  State   Description              h [kJ/kg]    P [kPa]'
    write(*,'(A)') '  ---'
    write(*,'(A,F12.2,4X,F10.2)') '  1       Compressor inlet (sat vap)', h1, P_evap
    write(*,'(A,F12.2,4X,F10.2)') '  2       Compressor exit (superheat)', h2, P_cond_p
    write(*,'(A,F12.2,4X,F10.2)') '  2s      Isentropic comp exit      ', h2s, P_cond_p
    write(*,'(A,F12.2,4X,F10.2)') '  3       Condenser exit (sat liq)  ', h3, P_cond_p
    write(*,'(A,F12.2,4X,F10.2)') '  4       Throttle exit (mixture)   ', h4, P_evap
    write(*,*)
    write(*,'(A)') '--- ENERGY ANALYSIS (per kg) ---'
    write(*,'(A,F12.2,A)')  '  Compressor Work (W)     = ', W_comp, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Refrigeration (Q_L)     = ', Q_L, ' kJ/kg'
    write(*,'(A,F12.2,A)')  '  Heat Rejected (Q_H)     = ', Q_H, ' kJ/kg'
    write(*,*)
    write(*,'(A)') '--- PERFORMANCE ---'
    write(*,'(A,F12.4)')    '  COP (refrigeration)     = ', COP_ref
    write(*,'(A,F12.4)')    '  COP (heat pump)         = ', COP_hp
    write(*,'(A,F12.4)')    '  Carnot COP (refrig)     = ', COP_carnot_ref
    write(*,'(A,F12.4)')    '  Carnot COP (heat pump)  = ', COP_carnot_hp
    write(*,'(A,F12.2,A)')  '  2nd Law Efficiency      = ', COP_ref/COP_carnot_ref*100, ' %'
    write(*,*)
    write(*,'(A)') '--- TOTAL CAPACITY ---'
    write(*,'(A,F12.2,A)')  '  Cooling Capacity        = ', Q_L_total, ' kW'
    write(*,'(A,F12.2,A)')  '  Cooling Capacity        = ', Q_L_total*3412.14d0/1000, ' BTU/hr (x1000)'
    write(*,'(A,F12.2,A)')  '  Cooling Capacity        = ', Q_L_total/3.517d0, ' tons'
    write(*,'(A,F12.2,A)')  '  Compressor Power        = ', W_comp_total, ' kW'
    write(*,'(A,F12.2,A)')  '  Heat Rejection          = ', m_dot * Q_H, ' kW'
    write(*,*)
    
    ! P-h diagram data
    write(*,'(A)') '--- P-h DIAGRAM DATA ---'
    write(*,'(A)') '  State    h [kJ/kg]   P [kPa]   Process'
    write(*,'(A)') '  ---'
    write(*,'(A,F10.2,4X,F10.2,4X,A)') '  1      ', h1, P_evap, 'Start (sat vapor)'
    write(*,'(A,F10.2,4X,F10.2,4X,A)') '  2      ', h2, P_cond_p, 'Compression'
    write(*,'(A,F10.2,4X,F10.2,4X,A)') '  3      ', h3, P_cond_p, 'Condensation'
    write(*,'(A,F10.2,4X,F10.2,4X,A)') '  4      ', h4, P_evap, 'Throttling'
    write(*,'(A,F10.2,4X,F10.2,4X,A)') '  1      ', h1, P_evap, 'Evaporation (back to 1)'
    write(*,*)
    write(*,'(A)') '--- NOTES ---'
    write(*,'(A)') '  * Refrigerant properties are approximated'
    write(*,'(A)') '  * For precise design, use REFPROP or CoolProp libraries'
    write(*,'(A)') '  * COP_HP = COP_REF + 1 (energy balance)'
    
end program refrigeration_cycle
