program participating_media_radiation_p1
    implicit none
    integer :: iostat_val
    double precision :: L_gap_m, T1_hot_K, T2_cold_K, a_abs_m, sig_scat_m, eps1_emiss, eps2_emiss
    double precision :: beta_ext_m, albedo_omega, tau_0, k_rosseland_WmK, T_mean_K
    double precision :: qr_p1_Wm2, qr_p1_kWm2, qr_transparent_kWm2, qr_opt_thick_kWm2
    double precision :: G0_Wm2, GL_Wm2, denom_p1
    double precision, parameter :: SIGMA = 5.670374d-8 ! W/(m2.K4)

    ! Read inputs
    read(*,*,iostat=iostat_val) L_gap_m       ! Medium Gap Thickness [m] (e.g. 0.50)
    read(*,*,iostat=iostat_val) T1_hot_K      ! Hot Wall Temp [K] (e.g. 1400.0)
    read(*,*,iostat=iostat_val) T2_cold_K     ! Cold Wall Temp [K] (e.g. 600.0)
    read(*,*,iostat=iostat_val) a_abs_m       ! Absorption Coefficient [1/m] (e.g. 2.0)
    read(*,*,iostat=iostat_val) sig_scat_m    ! Scattering Coefficient [1/m] (e.g. 1.0)
    read(*,*,iostat=iostat_val) eps1_emiss    ! Hot Surface Emissivity (e.g. 0.85)
    read(*,*,iostat=iostat_val) eps2_emiss    ! Cold Surface Emissivity (e.g. 0.80)

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

    if (L_gap_m <= 0.0d0 .or. T1_hot_K <= 0.0d0 .or. T2_cold_K <= 0.0d0 .or. a_abs_m < 0.0d0) then
        write(*,*) 'ERROR: Geometry, temperatures, and absorption must be positive.'
        stop
    end if

    beta_ext_m = a_abs_m + sig_scat_m
    albedo_omega = sig_scat_m / max(1.0d-5, beta_ext_m)
    tau_0 = beta_ext_m * L_gap_m

    T_mean_K = 0.5d0 * (T1_hot_K + T2_cold_K)
    k_rosseland_WmK = (16.0d0 * SIGMA * (T_mean_K**3)) / (3.0d0 * max(0.01d0, beta_ext_m))

    ! Optically Transparent Baseline (No medium)
    qr_transparent_kWm2 = (SIGMA * (T1_hot_K**4 - T2_cold_K**4) / &
                          (1.0d0/eps1_emiss + 1.0d0/eps2_emiss - 1.0d0)) / 1000.0d0

    ! P1 Approximation for 1D Planar Slab
    denom_p1 = (1.0d0/eps1_emiss - 0.5d0) + (1.0d0/eps2_emiss - 0.5d0) + 0.75d0 * tau_0
    qr_p1_Wm2 = (SIGMA * (T1_hot_K**4 - T2_cold_K**4)) / max(0.1d0, denom_p1)
    qr_p1_kWm2 = qr_p1_Wm2 / 1000.0d0

    ! Pure Rosseland Diffusion Limit
    qr_opt_thick_kWm2 = (k_rosseland_WmK * (T1_hot_K - T2_cold_K) / L_gap_m) / 1000.0d0

    G0_Wm2 = 4.0d0 * SIGMA * (T1_hot_K**4) - 2.0d0 * (2.0d0 - eps1_emiss)/eps1_emiss * qr_p1_Wm2
    GL_Wm2 = 4.0d0 * SIGMA * (T2_cold_K**4) + 2.0d0 * (2.0d0 - eps2_emiss)/eps2_emiss * qr_p1_Wm2

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — PARTICIPATING MEDIA RADIATION (P1 MODEL)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Gap Thickness / Extinction= ', L_gap_m, ' m / ', beta_ext_m, ' 1/m'
    write(*,'(A,F10.3,A,F10.3)')   'Optical Thickness tau_0 / w= ', tau_0, ' / ', albedo_omega
    write(*,'(A,F10.1,A,F10.1,A)') 'Hot Wall T1 / Cold Wall T2= ', T1_hot_K, ' K / ', T2_cold_K, ' K'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'P1 RADIATIVE HEAT FLUX qr = ', qr_p1_kWm2, ' kW/m2'
    write(*,'(A,F10.2,A)')  'Transparent Medium Flux   = ', qr_transparent_kWm2, ' kW/m2'
    write(*,'(A,F10.2,A)')  'Rosseland Diffusion Flux  = ', qr_opt_thick_kWm2, ' kW/m2'
    write(*,'(A,F10.3,A)')  'Rosseland Conductivity k  = ', k_rosseland_WmK, ' W/(m.K)'
    write(*,'(A,F10.2,A,F10.2,A)')'Incident Rad G(0) / G(L)  = ', G0_Wm2/1000.0d0, ' kW/m2 / ', GL_Wm2/1000.0d0, ' kW/m2'
    write(*,'(A)') '============================================================'

end program participating_media_radiation_p1