program gas_radiation
    implicit none

    ! Inputs
    double precision :: Tg_C, Ts_C
    double precision :: pc, pw, L, Pt

    ! Internal temperatures
    double precision :: Tg_K, Ts_K
    double precision :: tg, ts

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

    ! Emissivity variables
    double precision :: L_cm, pL_c, pL_w, pL_sum
    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
    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, temp_term

    integer :: i, j, iostat_val

    ! 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) Ts_C
    read(*,*,iostat=iostat_val) pc
    read(*,*,iostat=iostat_val) pw
    read(*,*,iostat=iostat_val) L
    read(*,*,iostat=iostat_val) Pt

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

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

    if (pc < 0.0d0 .or. pw < 0.0d0 .or. Pt <= 0.0d0) then
        write(*,*) 'ERROR: Gas pressures must be positive.'
        stop
    end if

    if ((pc + pw) > Pt) then
        write(*,*) 'ERROR: Sum of partial pressures cannot exceed total pressure.'
        stop
    end if

    if (L <= 0.0d0) then
        write(*,*) 'ERROR: Mean beam length must be positive.'
        stop
    end if

    ! Units conversion: path length to cm
    L_cm = L * 100.0d0
    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)

    ! =======================================================
    ! 1. EMISSIVITY CALCULATIONS (at gas temperature Tg)
    ! =======================================================

    ! Base Emissivity CO2
    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)

    ! Base Emissivity H2O
    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
    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 = 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 / (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 - overlap_eps
    if (eps_total < 0.0d0) eps_total = 0.0d0

    ! =======================================================
    ! 2. ABSORPTIVITY CALCULATIONS (for surface at Ts)
    ! =======================================================
    ts = Ts_K / 1000.0d0
    x_s = log(ts)

    ! CO2 Absorptivity (evaluated at Ts and scaled pL)
    pL_c_s = pc * L_cm * (Ts_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)
    alpha_c = Cc * eps_0c_s * (Tg_K / Ts_K)**0.65d0

    ! H2O Absorptivity
    pL_w_s = pw * L_cm * (Ts_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)
    alpha_w = Cw * eps_0w_s * (Tg_K / Ts_K)**0.45d0

    ! Overlap correction for absorptivity
    overlap_alpha = 0.0d0
    pL_sum_s = (pc + pw) * L_cm * (Ts_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 / (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

    ! =======================================================
    ! OUTPUT DATA
    ! =======================================================
    write(*,*) '============================================================'
    write(*,*) '       LECKNER MIXTURE GAS EMISSION ENGINE'
    write(*,*) '============================================================'
    write(*,*)
    write(*,'(A,F12.2,A)')  '  Gas Temperature (Tg)    = ', Tg_C, ' deg-C'
    write(*,'(A,F12.2,A)')  '  Wall Temperature (Ts)   = ', Ts_C, ' deg-C'
    write(*,'(A,F12.4,A)')  '  Path Length (L)         = ', L, ' m'
    write(*,'(A,F12.4,A)')  '  Total Pressure (Pt)     = ', Pt, ' bar'
    write(*,'(A,F12.4,A)')  '  CO2 Partial Pressure    = ', pc, ' bar'
    write(*,'(A,F12.4,A)')  '  H2O Partial Pressure    = ', pw, ' bar'
    write(*,*)
    write(*,*) '--- GAS EMISSIVITY PROPERTIES ------------------------------'
    write(*,'(A,F12.6)')    '  CO2 Base Emissivity     = ', eps_0c
    write(*,'(A,F12.6)')    '  CO2 Pressure Coeff Cc   = ', Cc
    write(*,'(A,F12.6)')    '  CO2 Corrected Emiss     = ', eps_c
    write(*,'(A,F12.6)')    '  H2O Base Emissivity     = ', eps_0w
    write(*,'(A,F12.6)')    '  H2O Pressure Coeff Cw   = ', Cw
    write(*,'(A,F12.6)')    '  H2O Corrected Emiss     = ', eps_w
    write(*,'(A,F12.6)')    '  Mixture Overlap Corr    = ', overlap_eps
    write(*,'(A,F12.6)')    '  Total Gas Emissivity    = ', eps_total
    write(*,*)
    write(*,*) '--- GAS ABSORPTIVITY PROPERTIES ----------------------------'
    write(*,'(A,F12.6)')    '  CO2 Absorptivity        = ', alpha_c
    write(*,'(A,F12.6)')    '  H2O Absorptivity        = ', alpha_w
    write(*,'(A,F12.6)')    '  Mixture Overlap Corr    = ', overlap_alpha
    write(*,'(A,F12.6)')    '  Total Gas Absorptivity  = ', alpha_total
    write(*,*)

end program gas_radiation