program absorption_refrig
    implicit none
    integer :: system_type, i, iostat_val
    double precision :: T_gen, T_abs, T_cond, T_evap
    double precision :: Q_gen, Q_abs, Q_cond, Q_evap, W_pump
    double precision :: COP, COP_carnot, eta_exergy
    double precision :: x_strong, x_weak, f_circ, m_ref, m_sol
    double precision :: h_ref_vap, h_ref_liq, h_sol_strong, h_sol_weak
    double precision :: P_high, P_low, dT_overlap
    double precision :: T_i, COP_i, f_i
    character(len=40) :: system_name

    read(*,*,iostat=iostat_val) system_type
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid system type input.'
        stop
    end if
    read(*,*,iostat=iostat_val) T_gen
    read(*,*,iostat=iostat_val) T_abs
    read(*,*,iostat=iostat_val) T_cond
    read(*,*,iostat=iostat_val) T_evap
    read(*,*,iostat=iostat_val) m_ref
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all absorption cycle inputs.'
        stop
    end if
    if (T_gen<=0.0d0.or.T_abs<=0.0d0.or.T_cond<=0.0d0.or.T_evap<=0.0d0) then
        write(*,*) 'ERROR: All temperatures must be positive (K).'
        stop
    end if
    if (T_gen <= T_cond) then
        write(*,*) 'ERROR: Generator temp must exceed condenser temp.'
        stop
    end if
    if (T_evap >= T_abs) then
        write(*,*) 'ERROR: Evaporator temp must be below absorber temp.'
        stop
    end if
    if (m_ref <= 0.0d0) m_ref = 1.0d0

    select case(system_type)
    case(1)
        system_name = 'Aqua-Ammonia (NH3-H2O)'
        call aqua_ammonia_cycle(T_gen, T_abs, T_cond, T_evap, m_ref, &
            x_strong, x_weak, f_circ, Q_gen, Q_abs, Q_cond, Q_evap, &
            W_pump, COP, P_high, P_low, &
            h_ref_vap, h_ref_liq, h_sol_strong, h_sol_weak)
    case(2)
        system_name = 'LiBr-H2O (water is refrigerant)'
        call libr_h2o_cycle(T_gen, T_abs, T_cond, T_evap, m_ref, &
            x_strong, x_weak, f_circ, Q_gen, Q_abs, Q_cond, Q_evap, &
            W_pump, COP, P_high, P_low, &
            h_ref_vap, h_ref_liq, h_sol_strong, h_sol_weak)
    case default
        write(*,*) 'ERROR: System type must be 1 (NH3-H2O) or 2 (LiBr-H2O).'
        stop
    end select

    ! Carnot COP for absorption
    ! COP_Carnot = (T_evap/(T_cond-T_evap)) * ((T_gen-T_abs)/T_gen)
    COP_carnot = (T_evap/(T_cond - T_evap)) * ((T_gen - T_abs)/T_gen)
    if (COP_carnot < 0.0d0) COP_carnot = 0.0d0
    eta_exergy = COP / max(COP_carnot, 1.0d-10)

    dT_overlap = T_abs - T_evap   ! temperature lift

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   ABSORPTION REFRIGERATION CYCLE ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')        '  System Type               = ', trim(system_name)
    write(*,'(A,F12.2,A)')  '  Generator Temperature     = ', T_gen, ' K'
    write(*,'(A,F12.2,A)')  '  Absorber Temperature      = ', T_abs, ' K'
    write(*,'(A,F12.2,A)')  '  Condenser Temperature     = ', T_cond, ' K'
    write(*,'(A,F12.2,A)')  '  Evaporator Temperature    = ', T_evap, ' K'
    write(*,'(A,ES12.4,A)') '  Refrigerant Flow          = ', m_ref, ' kg/s'
    write(*,*)
    write(*,'(A)') '--- SOLUTION PROPERTIES -------------------------------------'
    write(*,'(A,F10.4)')    '  Strong Solution x_s       = ', x_strong
    write(*,'(A,F10.4)')    '  Weak Solution x_w         = ', x_weak
    write(*,'(A,F10.4)')    '  Circulation Ratio f       = ', f_circ
    write(*,'(A,ES12.4,A)') '  Solution Flow (strong)    = ', m_ref*f_circ, ' kg/s'
    write(*,'(A,ES12.4,A)') '  High Pressure             = ', P_high, ' Pa'
    write(*,'(A,ES12.4,A)') '  Low Pressure              = ', P_low, ' Pa'
    write(*,*)
    write(*,'(A)') '--- HEAT DUTIES ---------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Generator Heat Qgen       = ', Q_gen, ' W'
    write(*,'(A,ES12.4,A)') '  Absorber Heat Qabs        = ', Q_abs, ' W'
    write(*,'(A,ES12.4,A)') '  Condenser Heat Qcond      = ', Q_cond, ' W'
    write(*,'(A,ES12.4,A)') '  Evaporator Heat Qevap     = ', Q_evap, ' W'
    write(*,'(A,ES12.4,A)') '  Pump Work Wpump           = ', W_pump, ' W'
    write(*,*)
    write(*,'(A)') '--- PERFORMANCE ---------------------------------------------'
    write(*,'(A,F10.4)')    '  COP (cooling)             = ', COP
    write(*,'(A,F10.4)')    '  COP Carnot (absorption)   = ', COP_carnot
    write(*,'(A,F10.4)')    '  Exergetic Efficiency      = ', eta_exergy
    write(*,'(A,F10.2,A)')  '  Temperature Lift          = ', dT_overlap, ' K'
    write(*,*)
    write(*,'(A)') '  Energy Balance Check:'
    write(*,'(A,ES12.4,A)') '    Qgen + Qevap + Wpump   = ', Q_gen+Q_evap+W_pump, ' W'
    write(*,'(A,ES12.4,A)') '    Qcond + Qabs            = ', Q_cond+Q_abs, ' W'
    write(*,*)

    ! COP vs Generator Temperature sweep
    write(*,'(A)') '--- COP VS GENERATOR TEMPERATURE SWEEP ----------------------'
    write(*,'(A)') '  T_gen[K]      COP           f_circ        COP_Carnot'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, 50
        T_i = T_cond + 5.0d0 + (T_gen + 50.0d0 - T_cond - 5.0d0)*dble(i-1)/49.0d0
        select case(system_type)
        case(1)
            call aqua_ammonia_cop(T_i, T_abs, T_cond, T_evap, COP_i, f_i)
        case(2)
            call libr_h2o_cop(T_i, T_abs, T_cond, T_evap, COP_i, f_i)
        end select
        write(*,'(F12.2,2X,F10.4,2X,F10.4,2X,F10.4)') T_i, COP_i, f_i, &
            (T_evap/(T_cond-T_evap))*((T_i-T_abs)/T_i)
    end do
    write(*,*)

    ! COP vs Evaporator Temperature sweep
    write(*,'(A)') '--- COP VS EVAPORATOR TEMPERATURE SWEEP ---------------------'
    write(*,'(A)') '  T_evap[K]     COP           COP_Carnot'
    write(*,'(A)') '  -------------------------------------------'
    do i = 1, 40
        T_i = 240.0d0 + (T_abs - 5.0d0 - 240.0d0)*dble(i-1)/39.0d0
        select case(system_type)
        case(1)
            call aqua_ammonia_cop(T_gen, T_abs, T_cond, T_i, COP_i, f_i)
        case(2)
            call libr_h2o_cop(T_gen, T_abs, T_cond, T_i, COP_i, f_i)
        end select
        write(*,'(F12.2,2X,F10.4,2X,F10.4)') T_i, COP_i, &
            (T_i/max(T_cond-T_i,1.0d0))*((T_gen-T_abs)/T_gen)
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  COP = Qevap / (Qgen + Wpump).'
    write(*,'(A)') '  COP_Carnot = (Te/(Tc-Te))*(Tg-Ta)/Tg.'
    write(*,'(A)') '  f = x_strong / (x_strong - x_weak).'
    write(*,'(A)') '  Simplified property correlations for educational use.'

contains

    subroutine aqua_ammonia_cycle(Tg, Ta, Tc, Te, mref, xs, xw, f, &
            Qg, Qa, Qc, Qe, Wp, cop, Ph, Pl, hrv, hrl, hss, hsw)
        implicit none
        double precision, intent(in) :: Tg, Ta, Tc, Te, mref
        double precision, intent(out) :: xs, xw, f, Qg, Qa, Qc, Qe, Wp, cop
        double precision, intent(out) :: Ph, Pl, hrv, hrl, hss, hsw
        double precision :: h_evap_nh3, Tg_C, Ta_C, Tc_C, Te_C

        Tg_C = Tg - 273.15d0
        Ta_C = Ta - 273.15d0
        Tc_C = Tc - 273.15d0
        Te_C = Te - 273.15d0

        ! Simplified NH3 saturation pressures
        Ph = 1.0d6 * exp(11.0d0 - 3100.0d0/Tc)   ! condenser
        Pl = 1.0d6 * exp(11.0d0 - 3100.0d0/Te)   ! evaporator
        if (Ph < Pl) Ph = Pl * 2.0d0

        ! Simplified concentration correlations
        xs = 0.60d0 - 0.002d0*(Ta_C - 30.0d0)    ! strong solution (rich in NH3)
        xw = 0.30d0 + 0.002d0*(Tg_C - 100.0d0)   ! weak solution (lean)
        if (xs > 0.70d0) xs = 0.70d0
        if (xs < 0.30d0) xs = 0.30d0
        if (xw > xs - 0.02d0) xw = xs - 0.02d0
        if (xw < 0.05d0) xw = 0.05d0

        f = xs / max(xs - xw, 1.0d-10)

        ! Simplified enthalpies (kJ/kg → J/kg)
        h_evap_nh3 = 1370.0d3   ! latent heat NH3
        hrv = 1450.0d3 + 2100.0d0*(Te_C - (-33.0d0))   ! vapor from evaporator
        hrl = 300.0d3 + 4500.0d0*(Tc_C - 25.0d0)       ! liquid from condenser
        hss = -100.0d3 + 3800.0d0*(Ta_C - 25.0d0)       ! strong solution
        hsw = 100.0d3 + 4200.0d0*(Tg_C - 100.0d0)       ! weak solution

        Qe = mref * h_evap_nh3
        Qc = mref * (hrv - hrl)
        Qg = mref * (hrv - hrl) + mref * f * (hsw - hss)
        if (Qg < Qe) Qg = Qe * 1.5d0
        Qa = Qg + Qe - Qc
        Wp = mref * f * 0.001d0 * (Ph - Pl) / 0.80d0   ! pump work
        cop = Qe / max(Qg + Wp, 1.0d-10)
    end subroutine aqua_ammonia_cycle

    subroutine libr_h2o_cycle(Tg, Ta, Tc, Te, mref, xs, xw, f, &
            Qg, Qa, Qc, Qe, Wp, cop, Ph, Pl, hrv, hrl, hss, hsw)
        implicit none
        double precision, intent(in) :: Tg, Ta, Tc, Te, mref
        double precision, intent(out) :: xs, xw, f, Qg, Qa, Qc, Qe, Wp, cop
        double precision, intent(out) :: Ph, Pl, hrv, hrl, hss, hsw
        double precision :: Tg_C, Ta_C, Tc_C, Te_C

        Tg_C = Tg - 273.15d0
        Ta_C = Ta - 273.15d0
        Tc_C = Tc - 273.15d0
        Te_C = Te - 273.15d0

        ! Water saturation pressures (simplified Antoine)
        Ph = 610.78d0 * exp(17.27d0*Tc_C/(Tc_C+237.3d0))
        Pl = 610.78d0 * exp(17.27d0*Te_C/(Te_C+237.3d0))
        if (Ph < Pl) Ph = Pl * 2.0d0

        ! LiBr concentrations (mass fraction LiBr)
        ! Strong = high LiBr (leaves generator), Weak = low LiBr (leaves absorber)
        xw = 0.55d0 + 0.001d0*(Ta_C - 30.0d0)
        xs = 0.62d0 + 0.001d0*(Tg_C - 80.0d0)
        if (xs > 0.70d0) xs = 0.70d0
        if (xs < xw + 0.02d0) xs = xw + 0.02d0
        if (xw < 0.40d0) xw = 0.40d0

        ! f = xs/(xs-xw) for LiBr system (based on LiBr mass balance)
        f = xs / max(xs - xw, 1.0d-10)

        ! Simplified water enthalpies
        hrv = 2500.0d3 + 1900.0d0*(Te_C - 5.0d0)
        hrl = 105.0d3 + 4186.0d0*(Tc_C - 25.0d0)
        hss = 150.0d3 + 2000.0d0*(Tg_C - 80.0d0)
        hsw = 80.0d3 + 2200.0d0*(Ta_C - 30.0d0)

        Qe = mref * (hrv - hrl)
        Qc = mref * (hrv - hrl)
        Qg = mref * hrv + mref*(f-1.0d0)*hss - mref*f*hsw
        if (Qg < 0.0d0) Qg = Qe * 1.4d0
        Qa = Qg + Qe - Qc
        Wp = mref * f * 0.001d0 * (Ph - Pl) / 0.75d0
        cop = Qe / max(Qg + Wp, 1.0d-10)
    end subroutine libr_h2o_cycle

    subroutine aqua_ammonia_cop(Tg, Ta, Tc, Te, cop_out, f_out)
        implicit none
        double precision, intent(in) :: Tg, Ta, Tc, Te
        double precision, intent(out) :: cop_out, f_out
        double precision :: xs, xw, Qg, Qa, Qc, Qe, Wp, Ph, Pl
        double precision :: hrv, hrl, hss, hsw
        call aqua_ammonia_cycle(Tg, Ta, Tc, Te, 1.0d0, xs, xw, f_out, &
            Qg, Qa, Qc, Qe, Wp, cop_out, Ph, Pl, hrv, hrl, hss, hsw)
    end subroutine aqua_ammonia_cop

    subroutine libr_h2o_cop(Tg, Ta, Tc, Te, cop_out, f_out)
        implicit none
        double precision, intent(in) :: Tg, Ta, Tc, Te
        double precision, intent(out) :: cop_out, f_out
        double precision :: xs, xw, Qg, Qa, Qc, Qe, Wp, Ph, Pl
        double precision :: hrv, hrl, hss, hsw
        call libr_h2o_cycle(Tg, Ta, Tc, Te, 1.0d0, xs, xw, f_out, &
            Qg, Qa, Qc, Qe, Wp, cop_out, Ph, Pl, hrv, hrl, hss, hsw)
    end subroutine libr_h2o_cop

end program absorption_refrig