program throttling_jt
    implicit none
    integer :: gas_type, iostat_val, i, n_sweep
    double precision :: T_inlet, P_inlet, P_outlet
    double precision :: Tc, Pc, omega, cp_gas, M_gas
    double precision :: Tc_c, Pc_c, omega_c, cp_c, M_c
    double precision :: R_gas, a_vdw, b_vdw
    double precision :: mu_JT, T_outlet, dT, T_inv
    double precision :: T_boil_est, P_H2O
    double precision :: T_sw, mu_sw
    character(len=40) :: gas_name, cooling_str, liq_str

    R_gas = 8.314462d0   ! J/(mol.K)

    ! ── Read inputs ─────────────────────────────────────────────
    read(*,*,iostat=iostat_val) gas_type
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid gas type input.'
        stop
    end if
    read(*,*,iostat=iostat_val) T_inlet
    read(*,*,iostat=iostat_val) P_inlet
    read(*,*,iostat=iostat_val) P_outlet
    read(*,*,iostat=iostat_val) Tc_c
    read(*,*,iostat=iostat_val) Pc_c
    read(*,*,iostat=iostat_val) omega_c
    read(*,*,iostat=iostat_val) cp_c
    read(*,*,iostat=iostat_val) M_c
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all inputs.'
        stop
    end if

    ! ── Input validation ────────────────────────────────────────
    if (T_inlet <= 0.0d0) T_inlet = 300.0d0
    if (P_inlet <= 0.0d0) P_inlet = 10.0d0
    if (P_outlet <= 0.0d0) P_outlet = 0.1d0
    if (P_outlet >= P_inlet) then
        write(*,*) 'ERROR: P_outlet must be less than P_inlet.'
        stop
    end if

    ! ── Gas properties ──────────────────────────────────────────
    if (gas_type == 1) then
        gas_name = 'Nitrogen (N2)'
        Tc = 126.2d0; Pc = 3.39d0; omega = 0.037d0
        cp_gas = 29.1d0; M_gas = 28.014d0
    else if (gas_type == 2) then
        gas_name = 'Oxygen (O2)'
        Tc = 154.6d0; Pc = 5.04d0; omega = 0.022d0
        cp_gas = 29.4d0; M_gas = 32.0d0
    else if (gas_type == 3) then
        gas_name = 'Carbon Dioxide (CO2)'
        Tc = 304.2d0; Pc = 7.38d0; omega = 0.224d0
        cp_gas = 37.1d0; M_gas = 44.01d0
    else if (gas_type == 4) then
        gas_name = 'Methane (CH4)'
        Tc = 190.6d0; Pc = 4.60d0; omega = 0.012d0
        cp_gas = 35.7d0; M_gas = 16.04d0
    else if (gas_type == 5) then
        gas_name = 'Hydrogen (H2)'
        Tc = 33.2d0; Pc = 1.30d0; omega = -0.217d0
        cp_gas = 28.8d0; M_gas = 2.016d0
    else if (gas_type == 6) then
        gas_name = 'Air'
        Tc = 132.5d0; Pc = 3.77d0; omega = 0.035d0
        cp_gas = 29.1d0; M_gas = 28.97d0
    else
        gas_name = 'Custom Gas'
        Tc = Tc_c; Pc = Pc_c; omega = omega_c
        cp_gas = cp_c; M_gas = M_c
        if (Tc <= 0.0d0) Tc = 200.0d0
        if (Pc <= 0.0d0) Pc = 5.0d0
        if (cp_gas <= 0.0d0) cp_gas = 30.0d0
        if (M_gas <= 0.0d0) M_gas = 28.0d0
    end if

    ! ── Van der Waals constants ─────────────────────────────────
    ! Pc in MPa -> Pa for calculation
    a_vdw = 27.0d0 * R_gas**2 * Tc**2 / (64.0d0 * Pc * 1.0d6)
    b_vdw = R_gas * Tc / (8.0d0 * Pc * 1.0d6)

    ! ── Joule-Thomson coefficient ───────────────────────────────
    ! mu_JT = (1/cp) * (2a/(RT) - b) in K/Pa
    mu_JT = (1.0d0 / cp_gas) * (2.0d0 * a_vdw / (R_gas * T_inlet) - b_vdw)
    ! Convert to K/MPa
    ! mu_JT is already in K/Pa, multiply by 1e6 for K/MPa
    ! Actually: dT = mu_JT [K/Pa] * dP [Pa]
    ! P_inlet, P_outlet in MPa, so dP in Pa = (P_outlet - P_inlet)*1e6

    ! Temperature change
    dT = mu_JT * (P_outlet - P_inlet) * 1.0d6
    T_outlet = T_inlet + dT

    ! Inversion temperature (van der Waals)
    T_inv = 2.0d0 * a_vdw / (R_gas * b_vdw)

    ! Boiling point estimate (rough: T_boil ~ 0.6 * Tc for most gases)
    T_boil_est = 0.6d0 * Tc

    ! Cooling or heating?
    if (mu_JT > 0.0d0) then
        cooling_str = 'COOLING (mu_JT > 0)'
    else if (mu_JT < 0.0d0) then
        cooling_str = 'HEATING (mu_JT < 0)'
    else
        cooling_str = 'NO CHANGE (ideal gas)'
    end if

    ! Liquefaction feasibility
    if (T_outlet < T_boil_est .and. T_outlet > 0.0d0) then
        liq_str = 'POSSIBLE — T_out < T_boil'
    else
        liq_str = 'NOT FEASIBLE at these conditions'
    end if

    ! ── Output ──────────────────────────────────────────────────
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   THROTTLING VALVE & JOULE-THOMSON EFFECT'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')        '  Gas                       = ', trim(gas_name)
    write(*,'(A,F12.2,A)')  '  Inlet Temperature T_in    = ', T_inlet, ' K'
    write(*,'(A,F12.4,A)')  '  Inlet Pressure P_in       = ', P_inlet, ' MPa'
    write(*,'(A,F12.4,A)')  '  Outlet Pressure P_out     = ', P_outlet, ' MPa'
    write(*,*)
    write(*,'(A)') '--- GAS PROPERTIES ------------------------------------------'
    write(*,'(A,F12.2,A)')  '  Critical Temperature Tc   = ', Tc, ' K'
    write(*,'(A,F12.4,A)')  '  Critical Pressure Pc      = ', Pc, ' MPa'
    write(*,'(A,F12.4)')    '  Acentric Factor omega     = ', omega
    write(*,'(A,F12.4,A)')  '  cp (molar)                = ', cp_gas, ' J/(mol.K)'
    write(*,'(A,F12.4,A)')  '  Molar Mass M              = ', M_gas, ' g/mol'
    write(*,*)
    write(*,'(A)') '--- VAN DER WAALS CONSTANTS ---------------------------------'
    write(*,'(A,ES14.6,A)') '  a                         = ', a_vdw, ' Pa.m^6/mol^2'
    write(*,'(A,ES14.6,A)') '  b                         = ', b_vdw, ' m^3/mol'
    write(*,*)
    write(*,'(A)') '--- JOULE-THOMSON COEFFICIENT -------------------------------'
    write(*,'(A,ES14.6,A)') '  mu_JT                     = ', mu_JT, ' K/Pa'
    write(*,'(A,F12.6,A)')  '  mu_JT                     = ', mu_JT*1.0d6, ' K/MPa'
    write(*,'(A,A)')        '  Effect                    = ', trim(cooling_str)
    write(*,*)
    write(*,'(A)') '--- TEMPERATURE CHANGE --------------------------------------'
    write(*,'(A,F12.4,A)')  '  Delta_T                   = ', dT, ' K'
    write(*,'(A,F12.2,A)')  '  T_outlet                  = ', T_outlet, ' K'
    write(*,'(A,F12.2,A)')  '  T_outlet                  = ', T_outlet-273.15d0, ' C'
    write(*,*)
    write(*,'(A)') '--- INVERSION TEMPERATURE -----------------------------------'
    write(*,'(A,F12.2,A)')  '  T_inversion (vdW max)     = ', T_inv, ' K'
    write(*,'(A,F12.2,A)')  '  T_inversion               = ', T_inv-273.15d0, ' C'
    if (T_inlet < T_inv) then
    write(*,'(A)')          '  T_in < T_inv => COOLING region'
    else
    write(*,'(A)')          '  T_in > T_inv => HEATING region'
    end if
    write(*,*)
    write(*,'(A)') '--- LIQUEFACTION FEASIBILITY --------------------------------'
    write(*,'(A,F12.2,A)')  '  T_boiling (est.)          = ', T_boil_est, ' K'
    write(*,'(A,A)')        '  Liquefaction              = ', trim(liq_str)
    write(*,*)
    write(*,'(A)') '--- IDEAL GAS COMPARISON ------------------------------------'
    write(*,'(A)')          '  Ideal gas: mu_JT = 0, no temperature change.'
    write(*,'(A,F12.4,A)')  '  Real gas Delta_T          = ', dT, ' K'
    write(*,*)

    ! ── Sensitivity sweep: T from 50 to 800 K ─────────────────
    n_sweep = 40
    write(*,'(A)') '--- SENSITIVITY: MU_JT VS INLET TEMPERATURE ----------------'
    write(*,'(A)') '  T_in[K]       mu_JT[K/MPa]    T_out[K]'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, n_sweep
        T_sw = 50.0d0 + dble(i-1) * (800.0d0 - 50.0d0) / dble(n_sweep - 1)
        mu_sw = (1.0d0/cp_gas) * (2.0d0*a_vdw/(R_gas*T_sw) - b_vdw)
        write(*,'(F10.2,4X,F14.6,4X,F12.2)') T_sw, mu_sw*1.0d6, &
            T_sw + mu_sw*(P_outlet-P_inlet)*1.0d6
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Van der Waals: a=27R^2Tc^2/(64Pc), b=RTc/(8Pc)'
    write(*,'(A)') '  mu_JT = (1/cp)*(2a/(RT) - b)  [K/Pa]'
    write(*,'(A)') '  T_inversion = 2a/(Rb)  (maximum inversion temp)'
    write(*,'(A)') '  Isenthalpic process: h_in = h_out'
    write(*,'(A)') '  T_out = T_in + mu_JT * (P_out - P_in)'
    write(*,'(A)') '  Ideal gas: mu_JT = 0 (no intermolecular forces)'

end program throttling_jt
