program linde_hampson
    implicit none
    integer :: gas_type, iostat_val, i, n_sweep
    double precision :: T_inlet, P_high, P_low, T_precool, eta_comp, mdot
    double precision :: Tc_g, Pc_g, Tb_g, hfg_liq, cp_g, M_g, R_spec
    double precision :: a_vdw, b_vdw, R_gas
    double precision :: mu_JT, T_inv, dT_JT, T_after_JT
    double precision :: h1, h2, hf, y_liq, dh_dep
    double precision :: W_comp, W_per_liq, COP_liq, COP_Carnot, FOM
    double precision :: P_sw, y_sw, Wl_sw, mu_sw, dh_sw, h2_sw
    character(len=30) :: gas_name

    R_gas = 8.314462d0

    read(*,*,iostat=iostat_val) gas_type
    if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid gas type.'; stop; end if
    read(*,*,iostat=iostat_val) T_inlet
    read(*,*,iostat=iostat_val) P_high
    read(*,*,iostat=iostat_val) P_low
    read(*,*,iostat=iostat_val) T_precool
    read(*,*,iostat=iostat_val) eta_comp
    read(*,*,iostat=iostat_val) mdot
    if (iostat_val /= 0) then; write(*,*) 'ERROR: Failed to read all inputs.'; stop; end if

    if (T_inlet <= 0.0d0) T_inlet = 300.0d0
    if (P_high <= 0.0d0) P_high = 20.0d0
    if (P_low <= 0.0d0) P_low = 0.1d0
    if (P_high <= P_low) then; write(*,*) 'ERROR: P_high must exceed P_low.'; stop; end if
    if (T_precool <= 0.0d0) T_precool = T_inlet
    if (eta_comp <= 0.0d0 .or. eta_comp > 1.0d0) eta_comp = 0.80d0
    if (mdot <= 0.0d0) mdot = 1.0d0

    select case(gas_type)
    case(1)
        gas_name='Nitrogen (N2)'; Tc_g=126.2d0; Pc_g=3.39d0
        Tb_g=77.4d0; hfg_liq=198.8d0; cp_g=1.04d0; M_g=28.014d0
    case(2)
        gas_name='Oxygen (O2)'; Tc_g=154.6d0; Pc_g=5.04d0
        Tb_g=90.2d0; hfg_liq=213.1d0; cp_g=0.918d0; M_g=32.0d0
    case(3)
        gas_name='Air'; Tc_g=132.5d0; Pc_g=3.77d0
        Tb_g=78.8d0; hfg_liq=201.0d0; cp_g=1.005d0; M_g=28.97d0
    case(4)
        gas_name='Methane (CH4)'; Tc_g=190.6d0; Pc_g=4.60d0
        Tb_g=111.7d0; hfg_liq=510.0d0; cp_g=2.22d0; M_g=16.04d0
    case default
        gas_name='Argon (Ar)'; Tc_g=150.7d0; Pc_g=4.86d0
        Tb_g=87.3d0; hfg_liq=161.1d0; cp_g=0.520d0; M_g=39.95d0
        gas_type=5
    end select

    R_spec = R_gas / M_g * 1000.0d0    ! J/(kg.K) -> convert to kJ/(kg.K) later
    ! a_vdw in Pa.m^6/mol^2, b_vdw in m^3/mol
    a_vdw = 27.0d0 * R_gas**2 * Tc_g**2 / (64.0d0 * Pc_g * 1.0d6)
    b_vdw = R_gas * Tc_g / (8.0d0 * Pc_g * 1.0d6)

    ! JT coefficient (molar basis, K/Pa)
    ! mu_JT_mol = (1/cp_mol)*(2a/(RT) - b)
    ! cp_mol = cp_g * M_g / 1000  [kJ/(mol.K)] -> J/(mol.K) = cp_g*M_g
    mu_JT = (1.0d0 / (cp_g * M_g)) * (2.0d0*a_vdw/(R_gas*T_precool) - b_vdw)
    ! mu_JT is K/Pa (per mole basis but cancels), convert effect to K
    ! For mass basis: same formula gives K/Pa for the gas

    T_inv = 2.0d0 * a_vdw / (R_gas * b_vdw)

    ! Temperature drop through JT valve (no HX, single pass)
    dT_JT = mu_JT * (P_high - P_low) * 1.0d6
    T_after_JT = T_precool + dT_JT  ! dT_JT is negative for cooling

    ! ── Linde cycle with perfect counter-flow HX ──────────────
    ! Enthalpy departure for vdW gas: dh = P*(b - 2a/(RT)) per mole
    ! Per kg: dh_dep = (P_high*1e6)*(b_vdw - 2*a_vdw/(R_gas*T_precool)) / (M_g/1000)
    ! This is in J/kg, convert to kJ/kg
    dh_dep = (P_high*1.0d6) * (b_vdw - 2.0d0*a_vdw/(R_gas*T_precool)) &
             / (M_g * 1.0d-3) / 1000.0d0   ! kJ/kg

    ! h1: returning gas at (T_inlet, P_low) — approximately ideal
    h1 = cp_g * T_inlet   ! kJ/kg (ref at 0 K)

    ! h2: gas before JT at (T_precool, P_high) — with departure
    h2 = cp_g * T_precool + dh_dep   ! kJ/kg

    ! hf: saturated liquid at (Tb, P_low)
    hf = cp_g * Tb_g - hfg_liq   ! kJ/kg

    ! Liquid yield
    if (abs(h1 - hf) > 1.0d-10) then
        y_liq = (h1 - h2) / (h1 - hf)
    else
        y_liq = 0.0d0
    end if
    if (y_liq < 0.0d0) y_liq = 0.0d0
    if (y_liq > 1.0d0) y_liq = 1.0d0

    ! Compressor work (isothermal, per kg total flow)
    ! W = R_spec * T_inlet * ln(P_high/P_low) / eta_comp [kJ/kg]
    W_comp = (R_gas / (M_g*1.0d-3)) * T_inlet * log(P_high/P_low) / (eta_comp * 1000.0d0)
    ! kJ/kg of total gas compressed

    ! Work per kg of liquid produced
    if (y_liq > 1.0d-10) then
        W_per_liq = W_comp / y_liq
    else
        W_per_liq = 1.0d10
    end if

    ! COP of liquefaction
    if (W_per_liq > 1.0d-10 .and. W_per_liq < 1.0d9) then
        COP_liq = hfg_liq / W_per_liq
    else
        COP_liq = 0.0d0
    end if

    ! Carnot COP for refrigeration from Tb to T_inlet
    COP_Carnot = Tb_g / max(T_inlet - Tb_g, 1.0d0)

    ! Figure of merit
    if (COP_Carnot > 1.0d-10) then
        FOM = COP_liq / COP_Carnot
    else
        FOM = 0.0d0
    end if

    ! ── Output ──────────────────────────────────────────────────
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   LINDE-HAMPSON GAS LIQUEFACTION CYCLE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')        '  Gas                       = ', trim(gas_name)
    write(*,'(A,F12.2,A)')  '  Inlet Temperature         = ', T_inlet, ' K'
    write(*,'(A,F12.4,A)')  '  High Pressure P_high      = ', P_high, ' MPa'
    write(*,'(A,F12.4,A)')  '  Low Pressure P_low        = ', P_low, ' MPa'
    write(*,'(A,F12.2,A)')  '  Pre-cooling Temperature   = ', T_precool, ' K'
    write(*,'(A,F10.4)')    '  Compressor Efficiency     = ', eta_comp
    write(*,'(A,F12.4,A)')  '  Make-up Flow Rate         = ', mdot, ' kg/s'
    write(*,*)
    write(*,'(A)') '--- GAS PROPERTIES ------------------------------------------'
    write(*,'(A,F12.2,A)')  '  Critical Temperature Tc   = ', Tc_g, ' K'
    write(*,'(A,F12.4,A)')  '  Critical Pressure Pc      = ', Pc_g, ' MPa'
    write(*,'(A,F12.2,A)')  '  Boiling Point Tb          = ', Tb_g, ' K'
    write(*,'(A,F12.2,A)')  '  hfg (liquid)              = ', hfg_liq, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  cp                        = ', cp_g, ' kJ/(kg.K)'
    write(*,'(A,F12.2,A)')  '  Molar Mass                = ', M_g, ' g/mol'
    write(*,*)
    write(*,'(A)') '--- JOULE-THOMSON ANALYSIS ----------------------------------'
    write(*,'(A,ES14.6,A)') '  mu_JT                     = ', mu_JT, ' K/Pa'
    write(*,'(A,F12.4,A)')  '  mu_JT                     = ', mu_JT*1.0d6, ' K/MPa'
    write(*,'(A,F12.2,A)')  '  Inversion Temperature     = ', T_inv, ' K'
    write(*,'(A,F12.4,A)')  '  Delta_T (JT, no HX)      = ', dT_JT, ' K'
    write(*,'(A,F12.2,A)')  '  T after JT (no HX)       = ', T_after_JT, ' K'
    if (T_precool < T_inv) then
    write(*,'(A)')          '  T_precool < T_inv => JT cooling ACTIVE'
    else
    write(*,'(A)')          '  WARNING: T_precool > T_inv — JT heating!'
    end if
    write(*,*)
    write(*,'(A)') '--- LINDE CYCLE (with counter-flow HX) ----------------------'
    write(*,'(A,F12.4,A)')  '  h1 (return gas, T_in)     = ', h1, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  h2 (before JT, T_pre)     = ', h2, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  hf (sat liquid, Tb)       = ', hf, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  Enthalpy departure        = ', dh_dep, ' kJ/kg'
    write(*,'(A,F10.6)')    '  Liquid Yield y            = ', y_liq
    write(*,'(A,F10.2,A)')  '  Liquid Yield              = ', y_liq*100.0d0, ' percent'
    write(*,'(A,F12.4,A)')  '  Liquid production rate    = ', mdot*y_liq, ' kg/s'
    write(*,*)
    write(*,'(A)') '--- ENERGY & PERFORMANCE ------------------------------------'
    write(*,'(A,F12.4,A)')  '  Compressor Work (per kg)  = ', W_comp, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  Work per kg liquid        = ', W_per_liq, ' kJ/kg_liq'
    write(*,'(A,F12.2,A)')  '  Total compressor power    = ', mdot*W_comp, ' kW'
    write(*,'(A,F10.6)')    '  COP liquefaction          = ', COP_liq
    write(*,'(A,F10.4)')    '  COP Carnot (Tb to Tin)    = ', COP_Carnot
    write(*,'(A,F10.6)')    '  Figure of Merit FOM       = ', FOM
    write(*,*)

    ! ── Sensitivity sweep: P_high ──────────────────────────────
    n_sweep = 40
    write(*,'(A)') '--- SENSITIVITY: YIELD VS HIGH PRESSURE --------------------'
    write(*,'(A)') '  P_high[MPa]   y[%]          W_liq[kJ/kg]'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, n_sweep
        P_sw = 2.0d0 + dble(i-1)*(40.0d0 - 2.0d0)/dble(n_sweep-1)
        mu_sw = (1.0d0/(cp_g*M_g))*(2.0d0*a_vdw/(R_gas*T_precool)-b_vdw)
        dh_sw = (P_sw*1.0d6)*(b_vdw-2.0d0*a_vdw/(R_gas*T_precool))/(M_g*1.0d-3)/1000.0d0
        h2_sw = cp_g*T_precool + dh_sw
        if (abs(h1-hf) > 1.0d-10) then
            y_sw = (h1 - h2_sw)/(h1 - hf)
        else
            y_sw = 0.0d0
        end if
        if (y_sw < 0.0d0) y_sw = 0.0d0
        if (y_sw > 1.0d0) y_sw = 1.0d0
        Wl_sw = (R_gas/(M_g*1.0d-3))*T_inlet*log(P_sw/P_low)/(eta_comp*1000.0d0)
        if (y_sw > 1.0d-10) then
            Wl_sw = Wl_sw / y_sw
        else
            Wl_sw = 99999.0d0
        end if
        write(*,'(F10.4,4X,F10.4,4X,F12.4)') P_sw, y_sw*100.0d0, Wl_sw
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Van der Waals EOS for enthalpy departure.'
    write(*,'(A)') '  mu_JT = (1/cp)(2a/(RT)-b) [K/Pa].'
    write(*,'(A)') '  Yield: y = (h1-h2)/(h1-hf) with counter-flow HX.'
    write(*,'(A)') '  Isothermal compressor: W = RT*ln(P2/P1)/eta.'
    write(*,'(A)') '  COP_liq = hfg / W_per_kg_liquid.'

end program linde_hampson
