program view_factors
    implicit none

    ! Inputs
    integer :: geom_type
    double precision :: dim1, dim2, dim3

    ! Constants
    double precision, parameter :: PI = 3.141592653589793d0

    ! Calculated properties
    double precision :: F12, F21
    double precision :: X, Y, S, term1, term2, term3, term4, term5, ln_term
    double precision :: r1, r2, A1, A2
    integer :: iostat_val

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

    read(*,*,iostat=iostat_val) dim1
    read(*,*,iostat=iostat_val) dim2
    read(*,*,iostat=iostat_val) dim3
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read dimension parameters.'
        stop
    end if

    ! Basic validations
    if (dim1 <= 0.0d0 .or. dim2 <= 0.0d0 .or. dim3 <= 0.0d0) then
        write(*,*) 'ERROR: All geometric dimensions must be strictly positive.'
        stop
    end if

    F12 = 0.0d0
    F21 = 0.0d0

    select case (geom_type)
    case (1)
        ! Aligned parallel rectangles: dim1 = width (W), dim2 = height (H), dim3 = separation (L)
        X = dim1 / dim3
        Y = dim2 / dim3

        term1 = log(((1.0d0 + X**2) * (1.0d0 + Y**2)) / (1.0d0 + X**2 + Y**2)) / 2.0d0
        term2 = X * sqrt(1.0d0 + Y**2) * atan(X / sqrt(1.0d0 + Y**2))
        term3 = Y * sqrt(1.0d0 + X**2) * atan(Y / sqrt(1.0d0 + X**2))
        term4 = X * atan(X)
        term5 = Y * atan(Y)

        F12 = (2.0d0 / (PI * X * Y)) * (term1 + term2 + term3 - term4 - term5)
        F21 = F12  ! Symmetric area

    case (2)
        ! Perpendicular rectangles with common edge: dim1 = height (H), dim2 = length (L), dim3 = common edge (W)
        X = dim1 / dim3
        Y = dim2 / dim3

        term1 = X * atan(1.0d0/X)
        term2 = Y * atan(1.0d0/Y)
        term3 = sqrt(X**2 + Y**2) * atan(1.0d0 / sqrt(X**2 + Y**2))
        
        ln_term = ((1.0d0 + X**2)*(1.0d0 + Y**2)) / (1.0d0 + X**2 + Y**2) * &
                  ( (X**2 * (1.0d0 + X**2 + Y**2)) / ((1.0d0 + X**2)*(X**2 + Y**2)) )**(X**2) * &
                  ( (Y**2 * (1.0d0 + X**2 + Y**2)) / ((1.0d0 + Y**2)*(X**2 + Y**2)) )**(Y**2)
        term4 = log(ln_term) / 4.0d0

        F12 = (1.0d0 / (PI * X)) * (term1 + term2 - term3 + term4)
        
        ! Reciprocity: A1 = H*W, A2 = L*W => F21 = F12 * A1 / A2 = F12 * H / L
        F21 = F12 * dim1 / dim2

    case (3)
        ! Coaxial parallel disks: dim1 = radius 1 (R1), dim2 = radius 2 (R2), dim3 = separation (L)
        r1 = dim1 / dim3
        r2 = dim2 / dim3
        S = 1.0d0 + (1.0d0 + r2**2) / (r1**2)

        F12 = 0.5d0 * (S - sqrt(S**2 - 4.0d0 * (dim2/dim1)**2))
        
        ! Reciprocity: A1 = PI*R1^2, A2 = PI*R2^2 => F21 = F12 * R1^2 / R2^2
        F21 = F12 * (dim1 / dim2)**2

    case (4)
        ! Concentric cylinders (infinite): dim1 = inner radius (r1), dim2 = outer radius (r2), dim3 = length (ignored)
        if (dim1 >= dim2) then
            write(*,*) 'ERROR: Inner radius r1 must be less than outer radius r2.'
            stop
        end if
        F12 = 1.0d0
        F21 = dim1 / dim2

    case (5)
        ! Concentric spheres: dim1 = inner radius (r1), dim2 = outer radius (r2), dim3 = ignored
        if (dim1 >= dim2) then
            write(*,*) 'ERROR: Inner radius r1 must be less than outer radius r2.'
            stop
        end if
        F12 = 1.0d0
        F21 = (dim1 / dim2)**2

    case default
        write(*,*) 'ERROR: Invalid geometry code.'
        stop
    end select

    ! Output results
    write(*,'(A)') '============================================================'
    write(*,'(A)') '            RADIATION VIEW FACTOR SOLVER'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A,I2)')         '  Geometry Code           = ', geom_type
    write(*,'(A,F12.4)')    '  Dimension 1 (d1)        = ', dim1
    write(*,'(A,F12.4)')    '  Dimension 2 (d2)        = ', dim2
    write(*,'(A,F12.4)')    '  Dimension 3 (d3)        = ', dim3
    write(*,*)
    write(*,'(A)') '--- CALCULATED VIEW FACTORS --------------------------------'
    write(*,'(A,F12.6)')    '  View Factor F1->2       = ', F12
    write(*,'(A,F12.6)')    '  View Factor F2->1       = ', F21
    write(*,*)
    write(*,'(A)') '--- RECIPROCITY VERIFICATION -------------------------------'
    if (geom_type == 1) then
        write(*,'(A)') '  Symmetric geometry: A1 = A2, F12 = F21.'
    else if (geom_type == 2) then
        write(*,'(A,F12.6)') '  Area ratio A1/A2 (H/L)  = ', dim1 / dim2
        write(*,'(A,F12.6)') '  F21 calculated          = ', F21
    else if (geom_type == 3) then
        write(*,'(A,F12.6)') '  Area ratio A1/A2        = ', (dim1 / dim2)**2
        write(*,'(A,F12.6)') '  F21 calculated          = ', F21
    else if (geom_type == 4) then
        write(*,'(A,F12.6)') '  Area ratio A1/A2 (r1/r2)= ', dim1 / dim2
        write(*,'(A,F12.6)') '  F21 calculated          = ', F21
    else if (geom_type == 5) then
        write(*,'(A,F12.6)') '  Area ratio A1/A2        = ', (dim1 / dim2)**2
        write(*,'(A,F12.6)') '  F21 calculated          = ', F21
    end if

end program view_factors