program conduction_shape_factors
    implicit none
    
    ! Pi constant
    real(8), parameter :: pi = 3.14159265358979d0
    
    ! Inputs
    integer :: config_idx
    real(8) :: L_val, z_val, r_val, r1_val, r2_val, w_val, D_val, D1_val, D2_val, e_val, s_val, Di_val, Do_val, ri_val, ro_val
    real(8) :: k_val
    real(8) :: T1_val, T2_val
    
    ! Outputs
    real(8) :: S_factor, R_th, Q_rate
    integer :: err_flag
    character(len=100) :: err_msg
    character(len=50) :: config_name
    
    ! Inline acosh function helper
    ! acosh(x) = log(x + sqrt(x^2 - 1.0)) for x >= 1.0
    real(8) :: x_temp
    
    err_flag = 0
    err_msg = ""
    S_factor = 0.0d0
    
    ! Read configuration index
    read *, config_idx
    
    ! Read global parameters: k, T1, T2
    read *, k_val
    read *, T1_val
    read *, T2_val
    
    select case (config_idx)
    case (1)
        ! 1. Buried horizontal cylinder in semi-infinite medium
        config_name = "Buried Horizontal Cylinder"
        read *, L_val  ! Length [m]
        read *, z_val  ! Depth [m]
        read *, r_val  ! Radius [m]
        
        if (L_val <= 0.0d0 .or. z_val <= 0.0d0 .or. r_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Dimensions must be positive values."
        else if (z_val < r_val) then
            err_flag = 1
            err_msg = "Burial depth (z) must be greater than or equal to cylinder radius (r)."
        else
            x_temp = z_val / r_val
            S_factor = 2.0d0 * pi * L_val / log(x_temp + sqrt(x_temp**2 - 1.0d0))
        end if
        
    case (2)
        ! 2. Buried vertical cylinder in semi-infinite medium
        config_name = "Buried Vertical Cylinder"
        read *, L_val  ! Length [m]
        read *, D_val  ! Diameter [m]
        
        if (L_val <= 0.0d0 .or. D_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Dimensions must be positive values."
        else if (L_val < 2.0d0 * D_val) then
            err_flag = 1
            err_msg = "Formula valid for L >= 2D (long cylinder/well)."
        else
            S_factor = 2.0d0 * pi * L_val / log(4.0d0 * L_val / D_val)
        end if
        
    case (3)
        ! 3. Two parallel buried cylinders in infinite medium
        config_name = "Two Parallel Cylinders"
        read *, L_val   ! Length [m]
        read *, w_val   ! Spacing between centers [m]
        read *, r1_val  ! Radius 1 [m]
        read *, r2_val  ! Radius 2 [m]
        
        if (L_val <= 0.0d0 .or. w_val <= 0.0d0 .or. r1_val <= 0.0d0 .or. r2_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Dimensions must be positive values."
        else if (w_val <= r1_val + r2_val) then
            err_flag = 1
            err_msg = "Centerline spacing (w) must be greater than r1 + r2 to prevent overlap."
        else
            x_temp = (w_val**2 - r1_val**2 - r2_val**2) / (2.0d0 * r1_val * r2_val)
            if (x_temp < 1.0d0) then
                err_flag = 1
                err_msg = "Invalid parameters for acosh calculation."
            else
                S_factor = 2.0d0 * pi * L_val / log(x_temp + sqrt(x_temp**2 - 1.0d0))
            end if
        end if
        
    case (4)
        ! 4. Isothermal sphere in semi-infinite medium
        config_name = "Buried Isothermal Sphere"
        read *, r_val  ! Radius [m]
        read *, z_val  ! Depth [m]
        
        if (r_val <= 0.0d0 .or. z_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Dimensions must be positive values."
        else if (z_val < r_val) then
            err_flag = 1
            err_msg = "Burial depth (z) must be greater than or equal to sphere radius (r)."
        else
            S_factor = 4.0d0 * pi * r_val / (1.0d0 - r_val / (2.0d0 * z_val))
        end if
        
    case (5)
        ! 5. Cylinder centered in square solid
        config_name = "Cylinder in Square Solid"
        read *, L_val  ! Length [m]
        read *, w_val  ! Width of square side [m]
        read *, D_val  ! Cylinder diameter [m]
        
        if (L_val <= 0.0d0 .or. w_val <= 0.0d0 .or. D_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Dimensions must be positive values."
        else if (w_val < D_val) then
            err_flag = 1
            err_msg = "Square side width (w) must be greater than cylinder diameter (D)."
        else
            S_factor = 2.0d0 * pi * L_val / log(1.08d0 * w_val / D_val)
        end if
        
    case (6)
        ! 6. Eccentric cylinders
        config_name = "Eccentric Concentric Cylinders"
        read *, L_val   ! Length [m]
        read *, D1_val  ! Inner cylinder outer diameter [m]
        read *, D2_val  ! Outer cylinder inner diameter [m]
        read *, e_val   ! Eccentricity (offset) [m]
        
        if (L_val <= 0.0d0 .or. D1_val <= 0.0d0 .or. D2_val <= 0.0d0 .or. e_val < 0.0d0) then
            err_flag = 1
            err_msg = "Dimensions must be positive, eccentricity non-negative."
        else if (D2_val <= D1_val + 2.0d0 * e_val) then
            err_flag = 1
            err_msg = "Outer diameter D2 must be larger than D1 + 2*e to prevent wall collision."
        else
            x_temp = (D1_val**2 + D2_val**2 - 4.0d0 * e_val**2) / (2.0d0 * D1_val * D2_val)
            if (x_temp < 1.0d0) then
                err_flag = 1
                err_msg = "Invalid parameters for acosh calculation."
            else
                S_factor = 2.0d0 * pi * L_val / log(x_temp + sqrt(x_temp**2 - 1.0d0))
            end if
        end if
        
    case (7)
        ! 7. Edge of two adjoining walls
        config_name = "Edge of Two Adjoining Walls"
        read *, D_val  ! Length of edge [m]
        
        if (D_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Edge length (D) must be positive."
        else
            S_factor = 0.54d0 * D_val
        end if
        
    case (8)
        ! 8. Corner of three walls
        config_name = "Corner of Three Walls"
        read *, L_val  ! Wall thickness [m]
        
        if (L_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Wall thickness (L) must be positive."
        else
            S_factor = 0.15d0 * L_val
        end if
        
    case (9)
        ! 9. Disk on semi-infinite surface
        config_name = "Isothermal Disk on Surface"
        read *, r_val  ! Disk radius [m]
        
        if (r_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Radius (r) must be positive."
        else
            S_factor = 4.0d0 * r_val
        end if
        
    case (10)
        ! 10. Hemisphere on surface
        config_name = "Hemisphere on Surface"
        read *, r_val  ! Radius [m]
        
        if (r_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Radius (r) must be positive."
        else
            S_factor = 2.0d0 * pi * r_val
        end if
        
    case (11)
        ! 11. Row of cylinders in semi-infinite solid
        config_name = "Row of Cylinders (Floor Heating)"
        read *, L_val  ! Length [m]
        read *, s_val  ! Center-to-center spacing [m]
        read *, D_val  ! Cylinder diameter [m]
        read *, z_val  ! Burial depth [m]
        
        if (L_val <= 0.0d0 .or. s_val <= 0.0d0 .or. D_val <= 0.0d0 .or. z_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Dimensions must be positive values."
        else if (s_val <= D_val) then
            err_flag = 1
            err_msg = "Spacing (s) must be greater than cylinder diameter (D)."
        else if (z_val < D_val / 2.0d0) then
            err_flag = 1
            err_msg = "Burial depth (z) must be at least half the diameter."
        else
            x_temp = (s_val / (pi * D_val)) * sinh(2.0d0 * pi * z_val / s_val)
            if (x_temp <= 0.0d0) then
                err_flag = 1
                err_msg = "Logarithm argument is non-positive."
            else
                S_factor = 2.0d0 * pi * L_val / log(x_temp)
            end if
        end if
        
    case (12)
        ! 12. Concentric cylinders
        config_name = "Concentric Cylinders (1D Radial)"
        read *, L_val   ! Length [m]
        read *, Di_val  ! Inner diameter [m]
        read *, Do_val  ! Outer diameter [m]
        
        if (L_val <= 0.0d0 .or. Di_val <= 0.0d0 .or. Do_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Dimensions must be positive values."
        else if (Do_val <= Di_val) then
            err_flag = 1
            err_msg = "Outer diameter (Do) must be greater than inner diameter (Di)."
        else
            S_factor = 2.0d0 * pi * L_val / log(Do_val / Di_val)
        end if
        
    case (13)
        ! 13. Concentric spheres
        config_name = "Concentric Spheres (1D Radial)"
        read *, ri_val  ! Inner radius [m]
        read *, ro_val  ! Outer radius [m]
        
        if (ri_val <= 0.0d0 .or. ro_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Dimensions must be positive values."
        else if (ro_val <= ri_val) then
            err_flag = 1
            err_msg = "Outer radius (ro) must be greater than inner radius (ri)."
        else
            S_factor = 4.0d0 * pi * ri_val * ro_val / (ro_val - ri_val)
        end if
        
    case default
        err_flag = 1
        err_msg = "Invalid configuration index chosen."
    end select
    
    ! Calculate thermal resistance and heat flow rate if valid
    if (err_flag == 0) then
        if (k_val <= 0.0d0) then
            err_flag = 1
            err_msg = "Thermal conductivity (k) must be positive."
        else
            R_th = 1.0d0 / (S_factor * k_val)
            Q_rate = S_factor * k_val * (T1_val - T2_val)
        end if
    end if
    
    ! Write output
    print *, '=================================================='
    print *, '  2D CONDUCTION SHAPE FACTOR CALCULATION REPORT'
    print *, '=================================================='
    print *, ''
    print *, 'CONFIGURATION DETAILS:'
    print *, '--------------------------------------------------'
    print *, '  Selected Profile: ', trim(config_name)
    print *, '  Index Code:       ', config_idx
    print *, ''
    
    if (err_flag /= 0) then
        print *, '❌ ERROR DETECTED:'
        print *, '  ', trim(err_msg)
        print *, '=================================================='
        stop
    end if
    
    print *, 'THERMO-PHYSICAL INPUTS:'
    print *, '--------------------------------------------------'
    print '(A,F10.4,A)', '  Conductivity of Solid (k):        ', k_val, ' W/(m.K)'
    print '(A,F10.2,A)', '  Boundary Temperature 1 (T1):      ', T1_val, ' C'
    print '(A,F10.2,A)', '  Boundary Temperature 2 (T2):      ', T2_val, ' C'
    print '(A,F10.2,A)', '  Temperature Difference (dT):      ', T1_val - T2_val, ' C'
    print *, ''
    
    print *, 'CALCULATED THERMAL METRICS:'
    print *, '--------------------------------------------------'
    print '(A,F12.6,A)', '  Conduction Shape Factor (S):      ', S_factor, ' m'
    print '(A,F12.6,A)', '  Conduction Thermal Resistance (R):', R_th, ' C/W'
    print '(A,F12.2,A)', '  Conduction Heat Flow Rate (Q):    ', Q_rate, ' W'
    print *, '=================================================='
    
end program conduction_shape_factors