Radiation in Participating Media (Non-Grey)

program participating_media
    implicit none

    ! Inputs
    double precision :: Tg_C, Tw_C
    double precision :: Pt, x_co2, x_h2o
    double precision :: V, Aw, eps_w

    ! Calculated parameters
    double precision :: Le, L_cm
    double precision :: pc, pw, pL_c, pL_w, pL_sum
    double precision :: Tg_K, Tw_K
    double precision :: tg, tw

    ! Leckner matrices
    double precision :: c(0:3, 0:2), d(0:3, 0:2)

    ! Emissivity variables
    double precision :: x_g, y_c, y_w
    double precision :: eps_0c, eps_0w
    double precision :: Pec, Pew, Cc, Cw
    double precision :: eps_c, eps_w_corr
    double precision :: zeta, overlap_eps, eps_total

    ! Absorptivity variables
    double precision :: pL_c_s, pL_w_s, pL_sum_s
    double precision :: x_s, y_c_s, y_w_s
    double precision :: eps_0c_s, eps_0w_s
    double precision :: alpha_c, alpha_w, overlap_alpha, alpha_total
    double precision :: Pec_s, Pew_s, Cc_s, Cw_s

    ! Heat Transfer variables
    double precision :: sigma, q_net, Q_total
    double precision :: temp_term
    integer :: i, j, iostat_val

    sigma = 5.670374d-8

    ! Initialize Leckner matrices
    ! CO2 Emissivity Matrix (c_ij)
    c(0,0) = -3.5430d0; c(0,1) = 1.6370d0;  c(0,2) = -0.5050d0
    c(1,0) = 0.3800d0;  c(1,1) = -0.2760d0; c(1,2) = 0.1060d0
    c(2,0) = -0.1180d0; c(2,1) = 0.0540d0;  c(2,2) = -0.0160d0
    c(3,0) = 0.0170d0;  c(3,1) = -0.0070d0; c(3,2) = 0.0020d0

    ! H2O Emissivity Matrix (d_ij)
    d(0,0) = -2.2118d0;   d(0,1) = -1.1987d0;  d(0,2) = 0.035596d0
    d(1,0) = 0.85667d0;   d(1,1) = 0.93048d0;  d(1,2) = -0.14391d0
    d(2,0) = -0.10833d0;  d(2,1) = -0.17156d0; d(2,2) = 0.045915d0
    d(3,0) = -0.016333d0; d(3,1) = 0.01988d0;  d(3,2) = -0.005830d0

    ! Read inputs
    read(*,*,iostat=iostat_val) Tg_C
    read(*,*,iostat=iostat_val) Tw_C
    read(*,*,iostat=iostat_val) Pt
    read(*,*,iostat=iostat_val) x_co2
    read(*,*,iostat=iostat_val) x_h2o
    read(*,*,iostat=iostat_val) V
    read(*,*,iostat=iostat_val) Aw
    read(*,*,iostat=iostat_val) eps_w

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read participating media parameters.'
        stop
    end if

    ! Validations
    Tg_K = Tg_C + 273.15d0
    Tw_K = Tw_C + 273.15d0
    if (Tg_K <= 0.0d0 .or. Tw_K <= 0.0d0) then
        write(*,*) 'ERROR: Temperatures must be above absolute zero.'
        stop
    end if

    if (Pt <= 0.0d0) then
        write(*,*) 'ERROR: Total pressure must be positive.'
        stop
    end if

    if (x_co2 < 0.0d0 .or. x_h2o < 0.0d0) then
        write(*,*) 'ERROR: Molar fractions must be non-negative.'
        stop
    end if

    if ((x_co2 + x_h2o) > 1.0d0) then
        write(*,*) 'ERROR: Sum of molar fractions cannot exceed 1.0.'
        stop
    end if

    if (V <= 0.0d0 .or. Aw <= 0.0d0) then
        write(*,*) 'ERROR: Volume and Area must be positive.'
        stop
    end if

    if (eps_w <= 0.0d0 .or. eps_w > 1.0d0) then
        write(*,*) 'ERROR: Wall emissivity must be in the range (0, 1].'
        stop
    end if

    ! Step 1: Calculate Mean Beam Length (Le)
    Le = 3.6d0 * V / Aw
    L_cm = Le * 100.0d0

    ! Step 2: Partial Pressures
    pc = x_co2 * Pt
    pw = x_h2o * Pt
    pL_c = pc * L_cm
    pL_w = pw * L_cm
    pL_sum = (pc + pw) * L_cm

    tg = Tg_K / 1000.0d0
    x_g = log(tg)

    ! =======================================================
    ! Step 3: GAS EMISSIVITY CALCULATIONS (at Tg)
    ! =======================================================
    ! CO2 Base Emissivity
    y_c = log(max(1.0001d0, pL_c))
    temp_term = 0.0d0
    do i = 0, 3
        do j = 0, 2
            temp_term = temp_term + c(i, j) * (x_g**i) * (y_c**j)
        end do
    end do
    eps_0c = exp(temp_term)

    ! H2O Base Emissivity
    y_w = log(max(1.0001d0, pL_w))
    temp_term = 0.0d0
    do i = 0, 3
        do j = 0, 2
            temp_term = temp_term + d(i, j) * (x_g**i) * (y_w**j)
        end do
    end do
    eps_0w = exp(temp_term)

    ! Pressure corrections at Tg
    Pec = Pt + 0.28d0 * pc
    Cc = 1.0d0 + ((Pec - 1.0d0) / (Pec - 1.0d0 + 0.28d0)) * 0.15d0 * &
         (1.0d0 - tanh(log10(max(1.0001d0, pL_c))))

    Pew = Pt + 0.85d0 * pw
    Cw = 1.0d0 + ((Pew - 1.0d0) / (Pew - 1.0d0 + 0.5d0)) * 0.35d0 * &
         (1.0d0 - tanh(log10(max(1.0001d0, pL_w))))

    eps_c = eps_0c * Cc
    eps_w_corr = eps_0w * Cw

    ! Overlap correction (zeta)
    overlap_eps = 0.0d0
    if ((pc + pw) > 1.0d-8) then
        zeta = pw / (pc + pw)
        if (pL_sum >= 1.0d0) then
            overlap_eps = (zeta / (10.7d0 + 101.0d0 * zeta) - 0.0089d0 * zeta**10.4d0) * &
                          (log10(pL_sum))**2.76d0
        end if
    end if

    eps_total = eps_c + eps_w_corr - overlap_eps
    if (eps_total < 0.0d0) eps_total = 0.0d0

    ! =======================================================
    ! Step 4: GAS ABSORPTIVITY CALCULATIONS (at Tw)
    ! =======================================================
    tw = Tw_K / 1000.0d0
    x_s = log(tw)

    ! CO2 Absorptivity (evaluated at Tw and scaled path length)
    pL_c_s = pc * L_cm * (Tw_K / Tg_K)
    y_c_s = log(max(1.0001d0, pL_c_s))
    temp_term = 0.0d0
    do i = 0, 3
        do j = 0, 2
            temp_term = temp_term + c(i, j) * (x_s**i) * (y_c_s**j)
        end do
    end do
    eps_0c_s = exp(temp_term)

    ! CO2 pressure correction for absorptivity (evaluated at scaled path length)
    Pec_s = Pt + 0.28d0 * pc
    Cc_s = 1.0d0 + ((Pec_s - 1.0d0) / (Pec_s - 1.0d0 + 0.28d0)) * 0.15d0 * &
         (1.0d0 - tanh(log10(max(1.0001d0, pL_c_s))))

    alpha_c = Cc_s * eps_0c_s * (Tg_K / Tw_K)**0.65d0

    ! H2O Absorptivity
    pL_w_s = pw * L_cm * (Tw_K / Tg_K)
    y_w_s = log(max(1.0001d0, pL_w_s))
    temp_term = 0.0d0
    do i = 0, 3
        do j = 0, 2
            temp_term = temp_term + d(i, j) * (x_s**i) * (y_w_s**j)
        end do
    end do
    eps_0w_s = exp(temp_term)

    ! H2O pressure correction for absorptivity (evaluated at scaled path length)
    Pew_s = Pt + 0.85d0 * pw
    Cw_s = 1.0d0 + ((Pew_s - 1.0d0) / (Pew_s - 1.0d0 + 0.5d0)) * 0.35d0 * &
         (1.0d0 - tanh(log10(max(1.0001d0, pL_w_s))))

    alpha_w = Cw_s * eps_0w_s * (Tg_K / Tw_K)**0.45d0

    ! Overlap correction for absorptivity
    overlap_alpha = 0.0d0
    pL_sum_s = (pc + pw) * L_cm * (Tw_K / Tg_K)
    if ((pc + pw) > 1.0d-8) then
        zeta = pw / (pc + pw)
        if (pL_sum_s >= 1.0d0) then
            overlap_alpha = (zeta / (10.7d0 + 101.0d0 * zeta) - 0.0089d0 * zeta**10.4d0) * &
                            (log10(pL_sum_s))**2.76d0
        end if
    end if

    alpha_total = alpha_c + alpha_w - overlap_alpha
    if (alpha_total < 0.0d0) alpha_total = 0.0d0

    ! =======================================================
    ! Step 5: NET HEAT EXCHANGE (Hottel grey wall formulation)
    ! =======================================================
    temp_term = (1.0d0 / eps_w) + alpha_total - 1.0d0
    if (temp_term > 1.0d-6) then
        q_net = (sigma * (eps_total * Tg_K**4 - alpha_total * Tw_K**4)) / temp_term
    else
        q_net = 0.0d0
    end if
    Q_total = q_net * Aw

    ! Output results in key-value format
    print *, "TG_C=", Tg_C
    print *, "TW_C=", Tw_C
    print *, "PT=", Pt
    print *, "X_CO2=", x_co2
    print *, "X_H2O=", x_h2o
    print *, "V=", V
    print *, "AW=", Aw
    print *, "EPS_W=", eps_w
    print *, "LE=", Le
    print *, "PC=", pc
    print *, "PW=", pw
    print *, "EPS_C=", eps_c
    print *, "EPS_W_CORR=", eps_w_corr
    print *, "OVERLAP_EPS=", overlap_eps
    print *, "EPS_G=", eps_total
    print *, "ALPHA_C=", alpha_c
    print *, "ALPHA_W=", alpha_w
    print *, "OVERLAP_ALPHA=", overlap_alpha
    print *, "ALPHA_G=", alpha_total
    print *, "QNET_W=", q_net
    print *, "Q_TOTAL=", Q_total

end program participating_media