program cyclone_particle_separator_leith_licht
    implicit none
    integer :: iostat_val, cyclone_model
    double precision :: Q_gas_m3h, T_gas_C, P_gas_bar, D_cyclone_m, rho_p_kgm3, dp_mean_um
    double precision :: rho_g, mu_g, Q_gas_m3s, a_inlet_m, b_inlet_m, De_exit_m
    double precision :: Vin_ms, Ne_turns, d50_um, d50_m, eta_grade_pct, NH_heads, dp_cyclone_Pa, dp_cyclone_mbar
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) cyclone_model   ! 1=Stairmand High-Efficiency, 2=Swift High-Efficiency, 3=Lapple General
    read(*,*,iostat=iostat_val) Q_gas_m3h       ! Gas Flow Rate [m3/h] (e.g. 3600.0)
    read(*,*,iostat=iostat_val) T_gas_C         ! Gas Temperature [deg C] (e.g. 150.0)
    read(*,*,iostat=iostat_val) P_gas_bar       ! Operating Pressure [bar abs] (e.g. 1.05)
    read(*,*,iostat=iostat_val) D_cyclone_m     ! Cyclone Barrel Diameter Dc [m] (e.g. 0.80)
    read(*,*,iostat=iostat_val) rho_p_kgm3      ! Dust Particle Density [kg/m3] (e.g. 2200.0)
    read(*,*,iostat=iostat_val) dp_mean_um      ! Mean Dust Particle Size [microns] (e.g. 10.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for cyclone calculation.'
        stop
    end if

    if (Q_gas_m3h <= 0.0d0 .or. D_cyclone_m <= 0.0d0 .or. rho_p_kgm3 <= 0.0d0) then
        write(*,*) 'ERROR: Flow rate, diameter and particle density must be positive.'
        stop
    end if

    ! Gas properties
    rho_g = (P_gas_bar * 1.0d5 * 0.02896d0) / (8.314d0 * (T_gas_C + 273.15d0))
    mu_g = 1.82d-5 * (((T_gas_C + 273.15d0) / 293.15d0)**0.7d0)
    Q_gas_m3s = Q_gas_m3h / 3600.0d0

    ! Stairmand High-Efficiency Proportions
    if (cyclone_model == 1) then ! Stairmand
        a_inlet_m = 0.50d0 * D_cyclone_m
        b_inlet_m = 0.20d0 * D_cyclone_m
        De_exit_m = 0.50d0 * D_cyclone_m
        Ne_turns = 6.0d0
        NH_heads = 6.4d0
    else if (cyclone_model == 2) then ! Swift
        a_inlet_m = 0.44d0 * D_cyclone_m
        b_inlet_m = 0.21d0 * D_cyclone_m
        De_exit_m = 0.40d0 * D_cyclone_m
        Ne_turns = 6.0d0
        NH_heads = 9.2d0
    else ! Lapple
        a_inlet_m = 0.50d0 * D_cyclone_m
        b_inlet_m = 0.25d0 * D_cyclone_m
        De_exit_m = 0.50d0 * D_cyclone_m
        Ne_turns = 5.0d0
        NH_heads = 8.0d0
    end if

    ! Inlet Velocity [m/s]
    Vin_ms = Q_gas_m3s / max(1.0d-5, a_inlet_m * b_inlet_m)

    ! Lapple / Leith-Licht d50 Cut Diameter [m]
    d50_m = sqrt((9.0d0 * mu_g * b_inlet_m) / (2.0d0 * PI * Ne_turns * Vin_ms * max(1.0d0, rho_p_kgm3 - rho_g)))
    d50_um = d50_m * 1.0d6

    ! Fractional Collection Efficiency for mean particle size
    eta_grade_pct = (1.0d0 / (1.0d0 + (d50_um / max(0.01d0, dp_mean_um))**2)) * 100.0d0

    ! Cyclone Pressure Drop [Pa & mbar]
    dp_cyclone_Pa = 0.5d0 * rho_g * (Vin_ms**2) * NH_heads
    dp_cyclone_mbar = dp_cyclone_Pa / 100.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — CYCLONE DUST SEPARATOR (LEITH-LICHT)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.2,A)') 'Gas Flow Rate / Barrel Diam  = ', Q_gas_m3h, ' m3/h / ', D_cyclone_m, ' m'
    write(*,'(A,F10.2,A,F10.1,A)') 'Inlet Velocity / Gas Temp    = ', Vin_ms, ' m/s / ', T_gas_C, ' deg C'
    write(*,'(A,F10.1,A,F10.1,A)') 'Particle Size / Solid Density= ', dp_mean_um, ' microns / ', rho_p_kgm3, ' kg/m3'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'CUT-OFF DIAMETER (d50)       = ', d50_um, ' microns'
    write(*,'(A,F10.2,A)')  'GRADE COLLECTION EFFICIENCY  = ', eta_grade_pct, ' %'
    write(*,'(A,F10.1,A,F10.2,A)') 'CYCLONE PRESSURE DROP (ΔP)   = ', dp_cyclone_mbar, ' mbar (', dp_cyclone_Pa/1000.0d0, ' kPa)'
    write(*,'(A,F10.2)')    'Inlet Velocity Heads (NH)    = ', NH_heads
    write(*,'(A)') '============================================================'

end program cyclone_particle_separator_leith_licht