Multidimensional Transient Conduction

program multidimensional_transient
    implicit none

    ! Variable declarations
    integer :: geom_type, i
    real(8) :: k_mat, rho_mat, cp_mat, alpha_diff
    real(8) :: h1, h2, h3
    real(8) :: temp_init, temp_inf
    real(8) :: dim1, dim2, dim3
    real(8) :: pos1, pos2, pos3
    real(8) :: time_val
    
    ! 1D results for components
    real(8) :: bi1, bi2, bi3
    real(8) :: fo1, fo2, fo3
    real(8) :: zeta1_1, zeta1_2, zeta1_3
    real(8) :: c1_1, c1_2, c1_3
    real(8) :: theta1, theta2, theta3
    real(8) :: q_ratio1, q_ratio2, q_ratio3
    
    ! Combined results
    real(8) :: theta_combined, temp_combined
    real(8) :: q_ratio_combined
    
    real(8), parameter :: pi = 3.141592653589793d0

    ! Read inputs
    read(*,*) geom_type      ! 1=Short Cylinder, 2=Rectangular Bar, 3=Rectangular Block, 4=Semi-infinite Cylinder, 5=Semi-infinite Rectangular Bar
    read(*,*) k_mat          ! Thermal conductivity [W/m.K]
    read(*,*) rho_mat        ! Density [kg/m3]
    read(*,*) cp_mat         ! Specific heat [J/kg.K]
    read(*,*) h1, h2, h3     ! Convection coefficients for directions 1, 2, 3
    read(*,*) temp_init      ! Initial temperature [deg-C]
    read(*,*) temp_inf       ! Ambient fluid temperature [deg-C]
    read(*,*) dim1, dim2, dim3 ! Specific geometry dimensions (m)
    read(*,*) pos1, pos2, pos3 ! Specific coordinates (m)
    read(*,*) time_val       ! Time elapsed [s]

    ! Safety limits
    if (k_mat <= 0.0d0) k_mat = 1.0d0
    if (rho_mat <= 0.0d0) rho_mat = 1000.0d0
    if (cp_mat <= 0.0d0) cp_mat = 1000.0d0
    if (time_val <= 0.0d0) time_val = 1.0d-5
    
    ! Thermal diffusivity
    alpha_diff = k_mat / (rho_mat * cp_mat)

    ! Initialize 1D component variables
    bi1 = 0.0d0; bi2 = 0.0d0; bi3 = 0.0d0
    fo1 = 0.0d0; fo2 = 0.0d0; fo3 = 0.0d0
    zeta1_1 = 0.0d0; zeta1_2 = 0.0d0; zeta1_3 = 0.0d0
    c1_1 = 0.0d0; c1_2 = 0.0d0; c1_3 = 0.0d0
    theta1 = 1.0d0; theta2 = 1.0d0; theta3 = 1.0d0
    q_ratio1 = 0.0d0; q_ratio2 = 0.0d0; q_ratio3 = 0.0d0
    
    ! Compute 1D components based on combined geometry choice
    select case (geom_type)
        case (1) ! Short Cylinder = Plane Wall (dim1) x Infinite Cylinder (dim2)
            ! Component 1: Plane Wall
            if (dim1 > 0.0d0) then
                bi1 = h1 * dim1 / k_mat
                fo1 = alpha_diff * time_val / (dim1**2)
                call find_eigenvalue(1, bi1, zeta1_1)
                call compute_coefficient(1, zeta1_1, c1_1)
                call compute_theta_position(1, zeta1_1, c1_1, fo1, pos1/dim1, theta1)
                call compute_energy_ratio(1, zeta1_1, c1_1, fo1, q_ratio1)
            end if
            
            ! Component 2: Infinite Cylinder
            if (dim2 > 0.0d0) then
                bi2 = h2 * dim2 / k_mat
                fo2 = alpha_diff * time_val / (dim2**2)
                call find_eigenvalue(2, bi2, zeta1_2)
                call compute_coefficient(2, zeta1_2, c1_2)
                call compute_theta_position(2, zeta1_2, c1_2, fo2, pos2/dim2, theta2)
                call compute_energy_ratio(2, zeta1_2, c1_2, fo2, q_ratio2)
            end if
            
            ! Combine
            theta_combined = theta1 * theta2
            q_ratio_combined = 1.0d0 - (1.0d0 - q_ratio1) * (1.0d0 - q_ratio2)

        case (2) ! Rectangular Bar = Wall 1 (dim1) x Wall 2 (dim2)
            ! Component 1: Wall 1
            if (dim1 > 0.0d0) then
                bi1 = h1 * dim1 / k_mat
                fo1 = alpha_diff * time_val / (dim1**2)
                call find_eigenvalue(1, bi1, zeta1_1)
                call compute_coefficient(1, zeta1_1, c1_1)
                call compute_theta_position(1, zeta1_1, c1_1, fo1, pos1/dim1, theta1)
                call compute_energy_ratio(1, zeta1_1, c1_1, fo1, q_ratio1)
            end if
            
            ! Component 2: Wall 2
            if (dim2 > 0.0d0) then
                bi2 = h2 * dim2 / k_mat
                fo2 = alpha_diff * time_val / (dim2**2)
                call find_eigenvalue(1, bi2, zeta1_2)
                call compute_coefficient(1, zeta1_2, c1_2)
                call compute_theta_position(1, zeta1_2, c1_2, fo2, pos2/dim2, theta2)
                call compute_energy_ratio(1, zeta1_2, c1_2, fo2, q_ratio2)
            end if
            
            ! Combine
            theta_combined = theta1 * theta2
            q_ratio_combined = 1.0d0 - (1.0d0 - q_ratio1) * (1.0d0 - q_ratio2)

        case (3) ! Rectangular Block = Wall 1 (dim1) x Wall 2 (dim2) x Wall 3 (dim3)
            ! Component 1: Wall 1
            if (dim1 > 0.0d0) then
                bi1 = h1 * dim1 / k_mat
                fo1 = alpha_diff * time_val / (dim1**2)
                call find_eigenvalue(1, bi1, zeta1_1)
                call compute_coefficient(1, zeta1_1, c1_1)
                call compute_theta_position(1, zeta1_1, c1_1, fo1, pos1/dim1, theta1)
                call compute_energy_ratio(1, zeta1_1, c1_1, fo1, q_ratio1)
            end if
            
            ! Component 2: Wall 2
            if (dim2 > 0.0d0) then
                bi2 = h2 * dim2 / k_mat
                fo2 = alpha_diff * time_val / (dim2**2)
                call find_eigenvalue(1, bi2, zeta1_2)
                call compute_coefficient(1, zeta1_2, c1_2)
                call compute_theta_position(1, zeta1_2, c1_2, fo2, pos2/dim2, theta2)
                call compute_energy_ratio(1, zeta1_2, c1_2, fo2, q_ratio2)
            end if
            
            ! Component 3: Wall 3
            if (dim3 > 0.0d0) then
                bi3 = h3 * dim3 / k_mat
                fo3 = alpha_diff * time_val / (dim3**2)
                call find_eigenvalue(1, bi3, zeta1_3)
                call compute_coefficient(1, zeta1_3, c1_3)
                call compute_theta_position(1, zeta1_3, c1_3, fo3, pos3/dim3, theta3)
                call compute_energy_ratio(1, zeta1_3, c1_3, fo3, q_ratio3)
            end if
            
            ! Combine
            theta_combined = theta1 * theta2 * theta3
            q_ratio_combined = 1.0d0 - (1.0d0 - q_ratio1) * (1.0d0 - q_ratio2) * (1.0d0 - q_ratio3)

        case (4) ! Semi-infinite Cylinder = Semi-infinite Solid (pos1) x Infinite Cylinder (dim2)
            ! Component 1: Semi-infinite Solid (evaluates at depth pos1)
            call compute_semi_infinite(pos1, h1, k_mat, alpha_diff, time_val, theta1)
            q_ratio1 = -1.0d0 ! Not applicable
            
            ! Component 2: Infinite Cylinder (dim2)
            if (dim2 > 0.0d0) then
                bi2 = h2 * dim2 / k_mat
                fo2 = alpha_diff * time_val / (dim2**2)
                call find_eigenvalue(2, bi2, zeta1_2)
                call compute_coefficient(2, zeta1_2, c1_2)
                call compute_theta_position(2, zeta1_2, c1_2, fo2, pos2/dim2, theta2)
                call compute_energy_ratio(2, zeta1_2, c1_2, fo2, q_ratio2)
            end if
            
            ! Combine
            theta_combined = theta1 * theta2
            q_ratio_combined = -1.0d0 ! Infinite Volume: N/A

        case (5) ! Semi-infinite Rectangular Bar = Semi-infinite Solid (pos1) x Plane Wall (dim2)
            ! Component 1: Semi-infinite Solid (evaluates at depth pos1)
            call compute_semi_infinite(pos1, h1, k_mat, alpha_diff, time_val, theta1)
            q_ratio1 = -1.0d0 ! Not applicable
            
            ! Component 2: Plane Wall (dim2)
            if (dim2 > 0.0d0) then
                bi2 = h2 * dim2 / k_mat
                fo2 = alpha_diff * time_val / (dim2**2)
                call find_eigenvalue(1, bi2, zeta1_2)
                call compute_coefficient(1, zeta1_2, c1_2)
                call compute_theta_position(1, zeta1_2, c1_2, fo2, pos2/dim2, theta2)
                call compute_energy_ratio(1, zeta1_2, c1_2, fo2, q_ratio2)
            end if
            
            ! Combine
            theta_combined = theta1 * theta2
            q_ratio_combined = -1.0d0 ! Infinite Volume: N/A

        case default
            print *, "ERROR: Invalid geometry type selection"
            stop
    end select
    
    ! Resulting Temperature
    temp_combined = temp_inf + (temp_init - temp_inf) * theta_combined

    ! ==================================================
    ! DISPLAY REPORT
    ! ==================================================
    print *, "=================================================="
    print *, "   MULTIDIMENSIONAL TRANSIENT CONDUCTION REPORT"
    print *, "               (Product Solution)"
    print *, "=================================================="
    print *, ""
    
    print *, "1. THERMOPHYSICAL PROPERTIES"
    print *, "--------------------------------------------------"
    print '(A, F10.4, A)', "  Conductivity (k):        ", k_mat, " W/m.K"
    print '(A, F10.2, A)', "  Density (rho):           ", rho_mat, " kg/m3"
    print '(A, F10.2, A)', "  Specific Heat (cp):      ", cp_mat, " J/kg.K"
    print '(A, ES12.4, A)',"  Thermal Diffusivity (a):  ", alpha_diff, " m2/s"
    print '(A, F10.2, A)', "  Initial Temp (Ti):       ", temp_init, " C"
    print '(A, F10.2, A)', "  Ambient Fluid Temp (Tinf):", temp_inf, " C"
    print *, ""
    
    print *, "2. GEOMETRY CONFIGURATION"
    print *, "--------------------------------------------------"
    select case (geom_type)
        case (1)
            print *, "  Type: Short Cylinder (Plane Wall x Infinite Cylinder)"
            print '(A, F10.4, A)', "  Wall Half-Thickness (L1):", dim1, " m"
            print '(A, F10.4, A)', "  Cylinder Radius (r0):     ", dim2, " m"
            print '(A, F10.4, A)', "  Axial Coordinate (x):    ", pos1, " m"
            print '(A, F10.4, A)', "  Radial Coordinate (r):   ", pos2, " m"
            print '(A, F10.2, A)', "  Wall Convection (h1):    ", h1, " W/m2.K"
            print '(A, F10.2, A)', "  Cylinder Convection (h2):", h2, " W/m2.K"
        case (2)
            print *, "  Type: Rectangular Bar (Wall 1 x Wall 2)"
            print '(A, F10.4, A)', "  Wall 1 Half-Thickness(L1):", dim1, " m"
            print '(A, F10.4, A)', "  Wall 2 Half-Thickness(L2):", dim2, " m"
            print '(A, F10.4, A)', "  Coordinate x1:           ", pos1, " m"
            print '(A, F10.4, A)', "  Coordinate x2:           ", pos2, " m"
            print '(A, F10.2, A)', "  Convection Face 1 (h1):  ", h1, " W/m2.K"
            print '(A, F10.2, A)', "  Convection Face 2 (h2):  ", h2, " W/m2.K"
        case (3)
            print *, "  Type: Rectangular Block (Wall 1 x Wall 2 x Wall 3)"
            print '(A, F10.4, A)', "  Wall 1 Half-Thickness(L1):", dim1, " m"
            print '(A, F10.4, A)', "  Wall 2 Half-Thickness(L2):", dim2, " m"
            print '(A, F10.4, A)', "  Wall 3 Half-Thickness(L3):", dim3, " m"
            print '(A, F10.4, A)', "  Coordinate x1:           ", pos1, " m"
            print '(A, F10.4, A)', "  Coordinate x2:           ", pos2, " m"
            print '(A, F10.4, A)', "  Coordinate x3:           ", pos3, " m"
            print '(A, F10.2, A)', "  Convection Face 1 (h1):  ", h1, " W/m2.K"
            print '(A, F10.2, A)', "  Convection Face 2 (h2):  ", h2, " W/m2.K"
            print '(A, F10.2, A)', "  Convection Face 3 (h3):  ", h3, " W/m2.K"
        case (4)
            print *, "  Type: Semi-Infinite Cylinder (Semi-Infinite Solid x Infinite Cylinder)"
            print '(A, F10.4, A)', "  Cylinder Radius (r0):     ", dim2, " m"
            print '(A, F10.4, A)', "  Semi-Infinite Depth (x): ", pos1, " m"
            print '(A, F10.4, A)', "  Radial Coordinate (r):   ", pos2, " m"
            print '(A, F10.2, A)', "  End Convection (h1):     ", h1, " W/m2.K"
            print '(A, F10.2, A)', "  Cylinder Convection (h2):", h2, " W/m2.K"
        case (5)
            print *, "  Type: Semi-Infinite Rectangular Bar (Semi-Infinite Solid x Plane Wall)"
            print '(A, F10.4, A)', "  Wall Half-Thickness (L1):", dim2, " m"
            print '(A, F10.4, A)', "  Semi-Infinite Depth (x): ", pos1, " m"
            print '(A, F10.4, A)', "  Wall Coordinate (y):     ", pos2, " m"
            print '(A, F10.2, A)', "  End Convection (h1):     ", h1, " W/m2.K"
            print '(A, F10.2, A)', "  Wall Convection (h2):    ", h2, " W/m2.K"
    end select
    print '(A, F10.2, A)', "  Elapsed Time (t):        ", time_val, " s"
    print *, ""

    print *, "3. 1D SUB-PROBLEM ANALYSIS"
    print *, "--------------------------------------------------"
    select case (geom_type)
        case (1)
            print *, "  [Sub-Problem 1: Plane Wall]"
            call print_1d_metrics(bi1, fo1, zeta1_1, c1_1, theta1)
            print *, ""
            print *, "  [Sub-Problem 2: Infinite Cylinder]"
            call print_1d_metrics(bi2, fo2, zeta1_2, c1_2, theta2)
        case (2)
            print *, "  [Sub-Problem 1: Plane Wall 1]"
            call print_1d_metrics(bi1, fo1, zeta1_1, c1_1, theta1)
            print *, ""
            print *, "  [Sub-Problem 2: Plane Wall 2]"
            call print_1d_metrics(bi2, fo2, zeta1_2, c1_2, theta2)
        case (3)
            print *, "  [Sub-Problem 1: Plane Wall 1]"
            call print_1d_metrics(bi1, fo1, zeta1_1, c1_1, theta1)
            print *, ""
            print *, "  [Sub-Problem 2: Plane Wall 2]"
            call print_1d_metrics(bi2, fo2, zeta1_2, c1_2, theta2)
            print *, ""
            print *, "  [Sub-Problem 3: Plane Wall 3]"
            call print_1d_metrics(bi3, fo3, zeta1_3, c1_3, theta3)
        case (4)
            print *, "  [Sub-Problem 1: Semi-Infinite Solid]"
            print '(A, F12.6)',    "    Dimensionless Temp (theta1):  ", theta1
            print *, ""
            print *, "  [Sub-Problem 2: Infinite Cylinder]"
            call print_1d_metrics(bi2, fo2, zeta1_2, c1_2, theta2)
        case (5)
            print *, "  [Sub-Problem 1: Semi-Infinite Solid]"
            print '(A, F12.6)',    "    Dimensionless Temp (theta1):  ", theta1
            print *, ""
            print *, "  [Sub-Problem 2: Plane Wall]"
            call print_1d_metrics(bi2, fo2, zeta1_2, c1_2, theta2)
    end select
    print *, ""

    print *, "4. COMBINED RESULTS (PRODUCT SOLUTION)"
    print *, "--------------------------------------------------"
    print '(A, F12.6)',    "  Dimensionless Temp (theta*):    ", theta_combined
    print '(A, F12.2, A)',  "  Combined Temperature T(pos, t): ", temp_combined, " C"
    
    if (q_ratio_combined >= 0.0d0) then
        print '(A, F12.6)',    "  Heat Transfer Ratio Q/Qmax:     ", q_ratio_combined
        print '(A, F12.2, A)',  "  Percentage Q/Qmax:              ", q_ratio_combined * 100.0d0, " %"
    else
        print *, "  Heat Transfer Ratio Q/Qmax:      N/A (Infinite Volume)"
    end if
    print *, ""
    
    print *, "=================================================="
    print *, "              CALCULATION COMPLETE"
    print *, "=================================================="

contains

    ! =========================================================
    ! PRINTS 1D SUB-PROBLEM METRICS WITH ACCURACY CHECKS
    ! =========================================================
    subroutine print_1d_metrics(bi, fo, zeta, c, theta)
        real(8), intent(in) :: bi, fo, zeta, c, theta
        print '(A, F12.4)', "    Biot Number (Bi):             ", bi
        print '(A, F12.4)', "    Fourier Number (Fo):          ", fo
        print '(A, F12.6)', "    First Eigenvalue (zeta1):     ", zeta
        print '(A, F12.6)', "    Coefficient (C1):             ", c
        print '(A, F12.6)', "    Dimensionless Temp (theta):   ", theta
        
        if (fo < 0.2d0) then
            print *, "    WARNING: Fo < 0.2. One-term approximation may be inaccurate."
        else
            print *, "    STATUS: Fo >= 0.2. One-term approximation is valid."
        end if
    end subroutine print_1d_metrics

    ! =========================================================
    ! STABLE EXP(U^2) * ERFC(U) TO AVOID FLOATING OVERFLOW
    ! =========================================================
    function exp_erfc_stable(u) result(val)
        real(8), intent(in) :: u
        real(8) :: val, u2
        
        if (u > 10.0d0) then
            u2 = u * u
            ! Asymptotic series expansion
            val = 1.0d0 / (u * sqrt(pi)) * (1.0d0 - 0.5d0 / u2 + 0.75d0 / (u2 * u2) - 1.875d0 / (u2 * u2 * u2))
        else
            val = exp(u * u) * erfc(u)
        end if
    end function exp_erfc_stable

    ! =========================================================
    ! SOLVES THE SEMI-INFINITE CONVECTION SUB-PROBLEM
    ! =========================================================
    subroutine compute_semi_infinite(x, h, k, alpha, t, theta)
        real(8), intent(in) :: x, h, k, alpha, t
        real(8), intent(out) :: theta
        real(8) :: arg1, arg2
        
        arg1 = x / (2.0d0 * sqrt(alpha * t))
        arg2 = h * sqrt(alpha * t) / k
        
        ! Math equation: theta = erf(arg1) + exp(h*x/k + h^2*alpha*t/k^2) * erfc(arg1 + arg2)
        ! Stably computed as: erf(arg1) + exp(-arg1^2) * exp_erfc_stable(arg1 + arg2)
        theta = erf(arg1) + exp(-arg1**2) * exp_erfc_stable(arg1 + arg2)
        
        if (theta < 0.0d0) theta = 0.0d0
        if (theta > 1.0d0) theta = 1.0d0
    end subroutine compute_semi_infinite

    ! =========================================================
    ! BESSEL FUNCTION J0 (Taylor series, 25 terms)
    ! =========================================================
    function bessel_j0(x) result(j0)
        real(8), intent(in) :: x
        real(8) :: j0, term, x2
        integer :: m

        j0 = 1.0d0
        term = 1.0d0
        x2 = (x / 2.0d0)**2

        do m = 1, 25
            term = -term * x2 / (dble(m)**2)
            j0 = j0 + term
            if (abs(term) < 1.0d-15) exit
        end do
    end function bessel_j0

    ! =========================================================
    ! BESSEL FUNCTION J1 (Taylor series, 25 terms)
    ! =========================================================
    function bessel_j1(x) result(j1)
        real(8), intent(in) :: x
        real(8) :: j1, term, x2
        integer :: m

        j1 = x / 2.0d0
        term = x / 2.0d0
        x2 = (x / 2.0d0)**2

        do m = 1, 25
            term = -term * x2 / (dble(m) * dble(m + 1))
            j1 = j1 + term
            if (abs(term) < 1.0d-15) exit
        end do
    end function bessel_j1

    ! =========================================================
    ! SPATIAL FACTOR at position x_r (x/L or r/r0)
    ! =========================================================
    function spatial_factor(gtype, zeta, x_r) result(sf)
        integer, intent(in) :: gtype
        real(8), intent(in) :: zeta, x_r
        real(8) :: sf

        select case (gtype)
            case (1) ! Plane wall: cos(zeta * x/L)
                sf = cos(zeta * x_r)
            case (2) ! Cylinder: J0(zeta * r/r0)
                sf = bessel_j0(zeta * x_r)
        end select
    end function spatial_factor

    ! =========================================================
    ! FIND FIRST EIGENVALUE via Newton-Raphson
    ! =========================================================
    subroutine find_eigenvalue(gtype, bi, zeta)
        integer, intent(in) :: gtype
        real(8), intent(in) :: bi
        real(8), intent(out) :: zeta
        real(8) :: f_val, df_val, correction
        real(8) :: j0v, j1v
        integer :: iter

        ! Initial guess
        select case (gtype)
            case (1) ! zeta * tan(zeta) = Bi
                if (bi < 0.01d0) then
                    zeta = sqrt(bi)
                else if (bi < 1.0d0) then
                    zeta = sqrt(bi) * (1.0d0 - bi/6.0d0)
                else if (bi < 10.0d0) then
                    zeta = 1.0d0 + 0.07d0 * bi
                else
                    zeta = pi/2.0d0 - 0.1d0
                end if

            case (2) ! zeta * J1(zeta) / J0(zeta) = Bi
                if (bi < 0.01d0) then
                    zeta = sqrt(2.0d0 * bi)
                else if (bi < 1.0d0) then
                    zeta = sqrt(2.0d0 * bi)
                else if (bi < 10.0d0) then
                    zeta = 1.5d0 + 0.05d0 * bi
                else
                    zeta = 2.4048d0 - 0.1d0
                end if
        end select

        ! Newton-Raphson iterations
        do iter = 1, 200
            select case (gtype)
                case (1) ! f = zeta * tan(zeta) - Bi
                    if (abs(cos(zeta)) < 1.0d-12) then
                        zeta = zeta - 0.01d0
                        cycle
                    end if
                    f_val = zeta * tan(zeta) - bi
                    df_val = tan(zeta) + zeta / (cos(zeta)**2)

                case (2) ! f = zeta * J1(zeta) / J0(zeta) - Bi
                    j0v = bessel_j0(zeta)
                    j1v = bessel_j1(zeta)
                    if (abs(j0v) < 1.0d-12) then
                        zeta = zeta - 0.01d0
                        cycle
                    end if
                    f_val = zeta * j1v / j0v - bi
                    ! Derivative using Bessel recurrences
                    df_val = j1v/j0v + zeta * (j0v * (j0v/zeta - j1v) - &
                             j1v * (-j1v)) / (j0v**2)
            end select

            if (abs(df_val) < 1.0d-15) exit

            correction = f_val / df_val

            ! Damping for stability
            if (abs(correction) > 0.5d0) then
                correction = sign(0.5d0, correction)
            end if

            zeta = zeta - correction

            ! Keep zeta positive and in first root range
            if (zeta <= 0.0d0) zeta = 0.01d0

            select case (gtype)
                case (1)
                    if (zeta >= pi/2.0d0) zeta = pi/2.0d0 - 0.001d0
                case (2)
                    if (zeta >= 2.4048d0) zeta = 2.4048d0 - 0.001d0
            end select

            if (abs(correction) < 1.0d-12) exit
        end do
    end subroutine find_eigenvalue

    ! =========================================================
    ! COMPUTE C1 COEFFICIENT
    ! =========================================================
    subroutine compute_coefficient(gtype, zeta, c1)
        integer, intent(in) :: gtype
        real(8), intent(in) :: zeta
        real(8), intent(out) :: c1

        select case (gtype)
            case (1) ! C1 = 4*sin(zeta) / (2*zeta + sin(2*zeta))
                c1 = 4.0d0 * sin(zeta) / (2.0d0 * zeta + sin(2.0d0 * zeta))
            case (2) ! C1 = 2 * J1(zeta) / (zeta * (J0^2 + J1^2))
                c1 = 2.0d0 / zeta * bessel_j1(zeta) / &
                     (bessel_j0(zeta)**2 + bessel_j1(zeta)**2)
        end select
    end subroutine compute_coefficient

    ! =========================================================
    ! COMPUTE THETA* AT A POSITION
    ! =========================================================
    subroutine compute_theta_position(gtype, zeta, c1, fo, x_r, theta)
        integer, intent(in) :: gtype
        real(8), intent(in) :: zeta, c1, fo, x_r
        real(8), intent(out) :: theta

        theta = c1 * exp(-(zeta**2) * fo) * spatial_factor(gtype, zeta, x_r)

        ! Clamp to valid range
        if (theta < 0.0d0) theta = 0.0d0
        if (theta > 1.0d0) theta = 1.0d0
    end subroutine compute_theta_position

    ! =========================================================
    ! COMPUTE ENERGY RATIO Q/Qmax
    ! =========================================================
    subroutine compute_energy_ratio(gtype, zeta, c1, fo, qr)
        integer, intent(in) :: gtype
        real(8), intent(in) :: zeta, c1, fo
        real(8), intent(out) :: qr
        real(8) :: theta0

        theta0 = c1 * exp(-(zeta**2) * fo)

        select case (gtype)
            case (1) ! Q/Qmax = 1 - theta0*sin(zeta)/zeta
                qr = 1.0d0 - theta0 * sin(zeta) / zeta
            case (2) ! Q/Qmax = 1 - 2*theta0*J1(zeta)/zeta
                qr = 1.0d0 - 2.0d0 * theta0 * bessel_j1(zeta) / zeta
        end select

        if (qr < 0.0d0) qr = 0.0d0
        if (qr > 1.0d0) qr = 1.0d0
    end subroutine compute_energy_ratio

end program multidimensional_transient