program rotary_regenerator_ljungstrom_aph
    implicit none
    integer :: iostat_val, fuel_type
    double precision :: m_dot_gas_kgs, Tg_in_C, m_dot_air_kgs, Ta_in_C
    double precision :: D_rotor_m, L_depth_m, N_rpm, matrix_mass_tonnes
    double precision :: cp_g, cp_a, cp_mat, C_g, C_a, C_min, C_max, C_star, C_r, C_r_star
    double precision :: h_gas, h_air, A_gas_m2, A_air_m2, NTU_g, NTU_a, NTU_0, eps_cf, eps_reg
    double precision :: Ta_out_C, Tg_out_C, Q_duty_MW, acid_dew_point_C
    character(len=32) :: corrosion_status
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) fuel_type       ! 1=Coal Boiler Flue Gas, 2=Heavy Fuel Oil, 3=Natural Gas Boiler, 4=Biomass
    read(*,*,iostat=iostat_val) m_dot_gas_kgs   ! Flue Gas Mass Flow [kg/s] (e.g. 180.0)
    read(*,*,iostat=iostat_val) Tg_in_C         ! Flue Gas Inlet Temp [deg C] (e.g. 340.0)
    read(*,*,iostat=iostat_val) m_dot_air_kgs   ! Combustion Air Mass Flow [kg/s] (e.g. 165.0)
    read(*,*,iostat=iostat_val) Ta_in_C         ! Combustion Air Inlet Temp [deg C] (e.g. 25.0)
    read(*,*,iostat=iostat_val) D_rotor_m       ! Rotor Outer Diameter [m] (e.g. 7.50)
    read(*,*,iostat=iostat_val) L_depth_m       ! Rotor Matrix Depth [m] (e.g. 1.80)
    read(*,*,iostat=iostat_val) N_rpm           ! Rotor Rotational Speed [rpm] (e.g. 1.50)
    read(*,*,iostat=iostat_val) matrix_mass_tonnes ! Rotor Matrix Mass [tonnes] (e.g. 45.0)

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

    if (m_dot_gas_kgs <= 0.0d0 .or. m_dot_air_kgs <= 0.0d0 .or. D_rotor_m <= 0.0d0 .or. N_rpm <= 0.0d0) then
        write(*,*) 'ERROR: Flow rates, rotor diameter, and rotational speed must be positive.'
        stop
    end if

    ! Flue gas properties and acid dew points
    if (fuel_type == 1) then ! Coal
        cp_g = 1080.0d0; cp_a = 1015.0d0; acid_dew_point_C = 135.0d0
    else if (fuel_type == 2) then ! Heavy Fuel Oil
        cp_g = 1090.0d0; cp_a = 1015.0d0; acid_dew_point_C = 140.0d0
    else if (fuel_type == 3) then ! Natural Gas
        cp_g = 1060.0d0; cp_a = 1015.0d0; acid_dew_point_C = 60.0d0
    else ! Biomass
        cp_g = 1075.0d0; cp_a = 1015.0d0; acid_dew_point_C = 120.0d0
    end if

    cp_mat = 480.0d0 ! Carbon / Corten weathering steel matrix

    C_g = m_dot_gas_kgs * cp_g
    C_a = m_dot_air_kgs * cp_a
    C_min = min(C_g, C_a)
    C_max = max(C_g, C_a)
    C_star = C_min / C_max

    C_r = (matrix_mass_tonnes * 1000.0d0) * cp_mat * (N_rpm / 60.0d0)
    C_r_star = C_r / C_min

    ! Matrix heat transfer area estimation (Compact corrugated elements: ~800 m2/m3)
    A_gas_m2 = 0.50d0 * (PI * (D_rotor_m**2) / 4.0d0) * L_depth_m * 850.0d0
    A_air_m2 = 0.50d0 * (PI * (D_rotor_m**2) / 4.0d0) * L_depth_m * 850.0d0

    h_gas = 65.0d0  ! W/(m2.K)
    h_air = 60.0d0  ! W/(m2.K)

    NTU_g = (h_gas * A_gas_m2) / C_g
    NTU_a = (h_air * A_air_m2) / C_a
    NTU_0 = 1.0d0 / ((1.0d0 / max(0.1d0, NTU_g)) + (1.0d0 / max(0.1d0, NTU_a)))

    ! Kays & London Regenerator Effectiveness Correlation
    if (C_star == 1.0d0) then
        eps_cf = NTU_0 / (1.0d0 + NTU_0)
    else
        eps_cf = (1.0d0 - exp(-NTU_0 * (1.0d0 - C_star))) / (1.0d0 - C_star * exp(-NTU_0 * (1.0d0 - C_star)))
    end if

    ! Matrix thermal capacity correction factor
    eps_reg = eps_cf * (1.0d0 - (1.0d0 / (9.0d0 * max(0.5d0, C_r_star)**1.93d0)))
    eps_reg = min(0.92d0, max(0.10d0, eps_reg))

    ! Outlet temperatures
    Ta_out_C = Ta_in_C + eps_reg * (C_min / C_a) * (Tg_in_C - Ta_in_C)
    Tg_out_C = Tg_in_C - eps_reg * (C_min / C_g) * (Tg_in_C - Ta_in_C)

    Q_duty_MW = (C_a * (Ta_out_C - Ta_in_C)) / 1.0d6

    if (Tg_out_C < acid_dew_point_C) then
        corrosion_status = 'WARNING: ACID DEW POINT CORROSION'
    else
        corrosion_status = 'SAFE (ABOVE ACID DEW POINT)'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — ROTARY REGENERATOR (LJUNGSTRÖM APH)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Rotor Diameter / Depth       = ', D_rotor_m, ' m / ', L_depth_m, ' m'
    write(*,'(A,F10.2,A,F10.2,A)') 'Rotational Speed / Mass      = ', N_rpm, ' rpm / ', matrix_mass_tonnes, ' tonnes'
    write(*,'(A,F10.2,A,F10.1,A)') 'Matrix Capacity Ratio (Cr*)  = ', C_r_star, ' (Effectiveness = ', eps_reg*100.0d0, '%)'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'HEAT RECOVERY DUTY (Q)       = ', Q_duty_MW, ' MWth'
    write(*,'(A,F10.1,A)')  'Combustion Air Preheated Tout= ', Ta_out_C, ' deg C'
    write(*,'(A,F10.1,A)')  'Flue Gas Cooled Exit Temp    = ', Tg_out_C, ' deg C'
    write(*,'(A,F10.1,A,A)')'Acid Dew Point Margin        = ', Tg_out_C - acid_dew_point_C, ' K | ', trim(corrosion_status)
    write(*,'(A)') '============================================================'

end program rotary_regenerator_ljungstrom_aph