program gray_radiation
    implicit none

    ! Inputs
    integer :: mode
    integer :: geom_type
    double precision :: T1_C, T2_C
    double precision :: A1, A2
    double precision :: eps1, eps2
    double precision :: F12

    ! Shields parameters
    integer :: n_shields
    double precision, dimension(5) :: eps_s
    double precision, dimension(5) :: Ts_K, Ts_C

    ! Constants
    double precision, parameter :: SIGMA = 5.670374d-8

    ! Calculated properties
    double precision :: T1_K, T2_K
    double precision :: Eb1, Eb2
    double precision :: R_s1, R_sp, R_s2, R_total
    double precision :: J1, J2, Q12
    double precision :: Q_no_shields, Q_with_shields, R_shields, reduction
    double precision :: temp_term
    integer :: i, iostat_val

    ! Read inputs
    read(*,*,iostat=iostat_val) mode
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid calculation mode.'
        stop
    end if

    if (mode == 1) then
        ! Mode 1: Two-surface enclosure
        read(*,*,iostat=iostat_val) geom_type
        read(*,*,iostat=iostat_val) T1_C
        read(*,*,iostat=iostat_val) T2_C
        read(*,*,iostat=iostat_val) A1
        read(*,*,iostat=iostat_val) A2
        read(*,*,iostat=iostat_val) eps1
        read(*,*,iostat=iostat_val) eps2
        read(*,*,iostat=iostat_val) F12
    else if (mode == 2) then
        ! Mode 2: Radiation shields
        read(*,*,iostat=iostat_val) T1_C
        read(*,*,iostat=iostat_val) T2_C
        read(*,*,iostat=iostat_val) A1
        read(*,*,iostat=iostat_val) eps1
        read(*,*,iostat=iostat_val) eps2
        read(*,*,iostat=iostat_val) n_shields
        do i = 1, 5
            read(*,*,iostat=iostat_val) eps_s(i)
        end do
    else
        write(*,*) 'ERROR: Invalid calculation mode.'
        stop
    end if

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

    ! Basic validation
    T1_K = T1_C + 273.15d0
    T2_K = T2_C + 273.15d0
    if (T1_K <= 0.0d0 .or. T2_K <= 0.0d0) then
        write(*,*) 'ERROR: Temperatures must be above absolute zero.'
        stop
    end if

    if (A1 <= 0.0d0) then
        write(*,*) 'ERROR: Area 1 must be positive.'
        stop
    end if

    if (eps1 <= 0.0d0 .or. eps1 > 1.0d0 .or. eps2 <= 0.0d0 .or. eps2 > 1.0d0) then
        write(*,*) 'ERROR: Surface emissivities must be in range ]0.0, 1.0].'
        stop
    end if

    Eb1 = SIGMA * T1_K**4
    Eb2 = SIGMA * T2_K**4

    ! Select Mode for logic execution
    if (mode == 1) then
        ! Mode 1: Two-surface gray enclosures
        select case (geom_type)
        case (1)
            ! Parallel plates
            if (F12 <= 0.0d0 .or. F12 > 1.0d0) then
                write(*,*) 'ERROR: View factor must be between 0.0 and 1.0.'
                stop
            end if
        case (2, 3)
            ! Concentric cylinders/spheres => F12 = 1.0
            F12 = 1.0d0
            if (A1 >= A2) then
                write(*,*) 'ERROR: Inner area A1 must be less than outer area A2.'
                stop
            end if
        case default
            write(*,*) 'ERROR: Invalid geometry type selected.'
            stop
        end select

        ! Resistor network resistances
        R_s1 = (1.0d0 - eps1) / (eps1 * A1)
        R_sp = 1.0d0 / (A1 * F12)
        R_s2 = (1.0d0 - eps2) / (eps2 * A2)
        R_total = R_s1 + R_sp + R_s2

        ! Heat exchange
        Q12 = (Eb1 - Eb2) / R_total

        ! Radiosities
        J1 = Eb1 - Q12 * R_s1
        J2 = Eb2 + Q12 * R_s2

        ! Output Results
        write(*,'(A)') '============================================================'
        write(*,'(A)') '      TWO-SURFACE GRAY ENCLOSURE THERMAL SOLVER'
        write(*,'(A)') '============================================================'
        write(*,*)
        write(*,'(A,F12.2,A)')  '  T1 (Blackbody Node Eb1) = ', T1_C, ' deg-C'
        write(*,'(A,F12.2,A)')  '  T2 (Blackbody Node Eb2) = ', T2_C, ' deg-C'
        write(*,'(A,ES12.4,A)') '  Eb1                     = ', Eb1, ' W/m2'
        write(*,'(A,ES12.4,A)') '  Eb2                     = ', Eb2, ' W/m2'
        write(*,*)
        write(*,'(A)') '--- RESISTOR NETWORK VALUES --------------------------------'
        write(*,'(A,ES12.4,A)') '  Surface Resistance R1   = ', R_s1, ' 1/m2'
        write(*,'(A,ES12.4,A)') '  Space Resistance R12    = ', R_sp, ' 1/m2'
        write(*,'(A,ES12.4,A)') '  Surface Resistance R2   = ', R_s2, ' 1/m2'
        write(*,'(A,ES12.4,A)') '  Total Resistance R_tot  = ', R_total, ' 1/m2'
        write(*,*)
        write(*,'(A)') '--- RADIOSITY NODES AND FLUX -------------------------------'
        write(*,'(A,ES12.4,A)') '  Radiosity Node J1       = ', J1, ' W/m2'
        write(*,'(A,ES12.4,A)') '  Radiosity Node J2       = ', J2, ' W/m2'
        write(*,'(A,ES12.4,A)') '  Net Heat Transfer Rate  = ', Q12, ' W'
        write(*,*)

    else if (mode == 2) then
        ! Mode 2: Multi-shield Radiation Screens
        if (n_shields < 1 .or. n_shields > 5) then
            write(*,*) 'ERROR: Number of shields must be between 1 and 5.'
            stop
        end if

        ! Calculate resistance without shields
        R_total = (1.0d0/eps1 + 1.0d0/eps2 - 1.0d0) / A1
        Q_no_shields = (Eb1 - Eb2) / R_total

        ! Calculate resistance of shields
        R_shields = 0.0d0
        do i = 1, n_shields
            if (eps_s(i) <= 0.0d0 .or. eps_s(i) > 1.0d0) then
                write(*,*) 'ERROR: Shield emissivity must be in range ]0, 1].'
                stop
            end if
            R_shields = R_shields + (2.0d0 / eps_s(i) - 1.0d0) / A1
        end do

        Q_with_shields = (Eb1 - Eb2) / (R_total + R_shields)
        reduction = 100.0d0 * (1.0d0 - Q_with_shields / Q_no_shields)

        ! Calculate shield temperatures sequentially
        ! Q = A1 * sigma * (Ts_{i-1}^4 - Ts_i^4) / (1/eps_s(i-1) + 1/eps_s(i) - 1)
        do i = 1, n_shields
            if (i == 1) then
                temp_term = T1_K**4 - (Q_with_shields / (A1 * SIGMA)) * (1.0d0/eps1 + 1.0d0/eps_s(1) - 1.0d0)
            else
                temp_term = Ts_K(i-1)**4 - (Q_with_shields / (A1 * SIGMA)) * (1.0d0/eps_s(i-1) + 1.0d0/eps_s(i) - 1.0d0)
            end if

            if (temp_term < 0.0d0) then
                Ts_K(i) = 0.0d0
            else
                Ts_K(i) = temp_term**0.25d0
            end if
            Ts_C(i) = Ts_K(i) - 273.15d0
        end do

        ! Output Results
        write(*,'(A)') '============================================================'
        write(*,'(A)') '        RADIATION SHIELDS SYSTEM SOLVER'
        write(*,'(A)') '============================================================'
        write(*,*)
        write(*,'(A,F12.2,A)')  '  T1 (Hot Surface)        = ', T1_C, ' deg-C'
        write(*,'(A,F12.2,A)')  '  T2 (Cold Surface)       = ', T2_C, ' deg-C'
        write(*,'(A,I2)')         '  Number of Shields       = ', n_shields
        write(*,'(A,ES12.4,A)') '  Q without Shields       = ', Q_no_shields, ' W'
        write(*,'(A,ES12.4,A)') '  Q with Shields          = ', Q_with_shields, ' W'
        write(*,'(A,F12.2,A)')  '  Heat Transfer Reduction = ', reduction, ' %'
        write(*,*)
        write(*,'(A)') '--- SHIELD TEMPERATURE DISTRIBUTION ------------------------'
        write(*,'(A)') '  Shield #      Emissivity     Temp [K]     Temp [deg-C]'
        write(*,'(A)') '  ----------------------------------------------------------'
        do i = 1, n_shields
            write(*,'(I5,10X,F6.3,9X,F8.2,7X,F8.2)') i, eps_s(i), Ts_K(i), Ts_C(i)
        end do
        write(*,*)
    end if

end program gray_radiation
