program multisurface_enclosure
    implicit none
    integer :: n, i, j
    real(8), allocatable :: A(:), eps(:), BC_val(:)
    integer, allocatable :: BC_type(:)
    real(8), allocatable :: F(:, :)
    real(8), allocatable :: AM(:, :), CM(:), J_vec(:)
    real(8), allocatable :: G(:), q_net(:), Q_net_val(:), Eb(:), T_C(:)
    real(8) :: sigma, q_sum
    logical :: error
    
    sigma = 5.670374d-8
    
    ! Read N
    read(*, *) n
    
    allocate(A(n), eps(n), BC_type(n), BC_val(n))
    allocate(F(n, n))
    allocate(AM(n, n), CM(n), J_vec(n))
    allocate(G(n), q_net(n), Q_net_val(n), Eb(n), T_C(n))
    
    ! Read surface properties
    do i = 1, n
        read(*, *) A(i), eps(i), BC_type(i), BC_val(i)
    end do
    
    ! Read view factor matrix
    do i = 1, n
        read(*, *) F(i, :)
    end do
    
    ! Construct matrix equations
    do i = 1, n
        if (BC_type(i) == 1) then
            ! Temperature imposed (BC_val is in C, convert to K)
            T_C(i) = BC_val(i)
            Eb(i) = sigma * (T_C(i) + 273.15d0)**4
            
            do j = 1, n
                if (i == j) then
                    AM(i, j) = 1.0d0 - (1.0d0 - eps(i)) * F(i, j)
                else
                    AM(i, j) = - (1.0d0 - eps(i)) * F(i, j)
                end if
            end do
            CM(i) = eps(i) * Eb(i)
        else
            ! Net flux density imposed (BC_val is in W/m2)
            q_net(i) = BC_val(i)
            T_C(i) = 0.0d0 ! Will be calculated later
            Eb(i) = 0.0d0 ! Will be calculated later
            
            do j = 1, n
                if (i == j) then
                    AM(i, j) = 1.0d0 - F(i, j)
                else
                    AM(i, j) = - F(i, j)
                end if
            end do
            CM(i) = q_net(i)
        end if
    end do
    
    ! Solve linear system
    call solve_system(n, AM, CM, J_vec, error)
    
    if (error) then
        print *, "ERROR: Linear system is singular or poorly conditioned."
        stop
    end if
    
    ! Post processing
    q_sum = 0.0d0
    do i = 1, n
        ! Irradiation G(i) = sum( F(i, j) * J(j) )
        G(i) = 0.0d0
        do j = 1, n
            G(i) = G(i) + F(i, j) * J_vec(j)
        end do
        
        if (BC_type(i) == 1) then
            ! Temp was imposed, calculate flux
            ! q_net = J - G
            q_net(i) = J_vec(i) - G(i)
            Q_net_val(i) = A(i) * q_net(i)
        else
            ! Flux was imposed, calculate temperature
            ! q_net = J - G (should match input flux, let's use input flux for stability)
            if (eps(i) > 1.0d-6) then
                Eb(i) = J_vec(i) + ((1.0d0 - eps(i)) / eps(i)) * q_net(i)
                if (Eb(i) > 0.0d0) then
                    T_C(i) = (Eb(i) / sigma)**0.25d0 - 273.15d0
                else
                    T_C(i) = -273.15d0
                end if
            else
                Eb(i) = 0.0d0
                T_C(i) = -273.15d0 ! Absolute zero if perfect reflector
            end if
            Q_net_val(i) = A(i) * q_net(i)
        end if
        q_sum = q_sum + Q_net_val(i)
    end do
    
    ! Output results in key-value format
    print *, "N=", n
    do i = 1, n
        print '(A,I0,A,F20.6)', "SURF_", i, "_AREA=", A(i)
        print '(A,I0,A,F20.6)', "SURF_", i, "_EPS=", eps(i)
        print '(A,I0,A,I0)', "SURF_", i, "_BCTYPE=", BC_type(i)
        print '(A,I0,A,F20.6)', "SURF_", i, "_T=", T_C(i)
        print '(A,I0,A,F20.6)', "SURF_", i, "_QNET_DEN=", q_net(i)
        print '(A,I0,A,F20.6)', "SURF_", i, "_QNET=", Q_net_val(i)
        print '(A,I0,A,F20.6)', "SURF_", i, "_J=", J_vec(i)
        print '(A,I0,A,F20.6)', "SURF_", i, "_G=", G(i)
        print '(A,I0,A,F20.6)', "SURF_", i, "_EB=", Eb(i)
    end do
    print '(A,F20.8)', "Q_SUM=", q_sum
    
contains

    subroutine solve_system(n, A_mat, C_vec, X_vec, err)
        integer, intent(in) :: n
        real(8), intent(inout) :: A_mat(n, n)
        real(8), intent(inout) :: C_vec(n)
        real(8), intent(out) :: X_vec(n)
        logical, intent(out) :: err
        integer :: i, j, k, max_row
        real(8) :: factor, pivot, temp_val
        real(8) :: temp_row(n)
        
        err = .false.
        X_vec = 0.0d0
        
        do i = 1, n
            ! Find pivot with partial pivoting
            max_row = i
            do k = i + 1, n
                if (abs(A_mat(k, i)) > abs(A_mat(max_row, i))) then
                    max_row = k
                end if
            end do
            
            ! Swap rows if necessary
            if (max_row /= i) then
                temp_row = A_mat(i, :)
                A_mat(i, :) = A_mat(max_row, :)
                A_mat(max_row, :) = temp_row
                
                temp_val = C_vec(i)
                C_vec(i) = C_vec(max_row)
                C_vec(max_row) = temp_val
            end if
            
            pivot = A_mat(i, i)
            if (abs(pivot) < 1.0d-12) then
                err = .true.
                return
            end if
            
            ! Normalize row i
            C_vec(i) = C_vec(i) / pivot
            A_mat(i, :) = A_mat(i, :) / pivot
            
            ! Eliminate other rows
            do j = 1, n
                if (j /= i) then
                    factor = A_mat(j, i)
                    A_mat(j, :) = A_mat(j, :) - factor * A_mat(i, :)
                    C_vec(j) = C_vec(j) - factor * C_vec(i)
                end if
            end do
        end do
        
        X_vec = C_vec
    end subroutine solve_system

end program multisurface_enclosure
