program transient_conduction_1d
    implicit none

    ! Variable declarations
    integer :: geom_type, n_points, i, n_terms
    real(8) :: dim_char, k_mat, rho_mat, cp_mat, h_conv
    real(8) :: temp_init, temp_inf
    real(8) :: alpha_diff, biot, fourier, time_val
    real(8) :: x_ratio
    real(8) :: zeta1, coeff_c1, theta_star, theta_center
    real(8) :: temp_result, q_ratio
    real(8) :: pi, fo_val, x_pos, temp_x
    real(8) :: zeta_n, delta_fo
    character(len=20) :: geom_name

    pi = 3.14159265358979d0

    ! Read inputs
    read *, geom_type    ! 1=Plane Wall, 2=Long Cylinder, 3=Sphere
    read *, dim_char     ! Half-thickness L (wall) or radius r0
    read *, k_mat        ! Thermal conductivity [W/m.K]
    read *, rho_mat      ! Density [kg/m3]
    read *, cp_mat       ! Specific heat [J/kg.K]
    read *, h_conv       ! Convection coefficient [W/m2.K]
    read *, temp_init    ! Initial temperature [deg-C]
    read *, temp_inf     ! Ambient temperature [deg-C]
    read *, x_ratio      ! Position ratio x/L or r/r0 (0=center, 1=surface)
    read *, time_val     ! Time [s]

    select case (geom_type)
        case (1)
            geom_name = "PLANE WALL"
        case (2)
            geom_name = "LONG CYLINDER"
        case (3)
            geom_name = "SPHERE"
        case default
            print *, 'ERROR: Invalid geometry type'
            stop
    end select

    ! Thermal diffusivity
    alpha_diff = k_mat / (rho_mat * cp_mat)

    ! Biot number
    biot = h_conv * dim_char / k_mat

    ! Fourier number
    if (time_val > 0.0d0) then
        fourier = alpha_diff * time_val / (dim_char**2)
    else
        fourier = 0.0d0
    end if

    ! =============================================
    ! DISPLAY HEADER
    ! =============================================
    print *, '================================================='
    print *, '   TRANSIENT CONDUCTION - ONE-TERM APPROXIMATION'
    print *, '   (Heisler Chart Method)'
    print *, '================================================='
    print *, ''
    print *, '   Geometry: ', trim(geom_name)
    print *, ''

    print *, '================================================='
    print *, '          INPUT PARAMETERS'
    print *, '================================================='
    print *, ''

    select case (geom_type)
        case (1)
            print '(A,F12.6,A)', '   Half-thickness (L):        ', dim_char, ' m'
            print '(A,F12.4,A)', '                              ', dim_char*1000.0d0, ' mm'
        case (2)
            print '(A,F12.6,A)', '   Radius (r0):               ', dim_char, ' m'
            print '(A,F12.4,A)', '                              ', dim_char*1000.0d0, ' mm'
        case (3)
            print '(A,F12.6,A)', '   Radius (r0):               ', dim_char, ' m'
            print '(A,F12.4,A)', '                              ', dim_char*1000.0d0, ' mm'
    end select

    print '(A,F12.4,A)', '   Thermal Conductivity (k):  ', k_mat, ' W/m.K'
    print '(A,F12.2,A)', '   Density (rho):             ', rho_mat, ' kg/m3'
    print '(A,F12.2,A)', '   Specific Heat (cp):        ', cp_mat, ' J/kg.K'
    print '(A,F12.2,A)', '   Convection Coeff. (h):     ', h_conv, ' W/m2.K'
    print '(A,F12.2,A)', '   Initial Temp. (Ti):        ', temp_init, ' deg-C'
    print '(A,F12.2,A)', '   Ambient Temp. (T_inf):     ', temp_inf, ' deg-C'
    print '(A,F12.4)',   '   Position ratio (x/L or r/r0):', x_ratio
    print '(A,F12.4,A)', '   Time:                      ', time_val, ' s'
    print *, ''

    ! =============================================
    ! DIMENSIONLESS NUMBERS
    ! =============================================
    print *, '================================================='
    print *, '       DIMENSIONLESS PARAMETERS'
    print *, '================================================='
    print *, ''
    print '(A,ES14.6,A)', '   Thermal Diffusivity (a):   ', alpha_diff, ' m2/s'
    print '(A,F14.6)',    '   Biot Number (Bi):          ', biot
    print '(A,F14.6)',    '   Fourier Number (Fo):       ', fourier
    print *, ''

    ! Validity check
    if (fourier < 0.2d0 .and. time_val > 0.0d0) then
        print *, '   WARNING: Fo < 0.2'
        print *, '   -------------------------------------------'
        print *, '   The one-term approximation may not be'
        print *, '   accurate. Consider using the full series'
        print *, '   solution or numerical methods.'
        print *, '   Accuracy improves significantly for Fo > 0.2'
    else if (time_val > 0.0d0) then
        print *, '   STATUS: Fo >= 0.2'
        print *, '   -------------------------------------------'
        print *, '   One-term approximation IS VALID.'
        print *, '   Error is typically less than 2%.'
    end if
    print *, ''

    if (biot < 0.1d0) then
        print *, '   NOTE: Bi < 0.1'
        print *, '   The lumped system analysis could also be'
        print *, '   used for this problem (simpler method).'
        print *, ''
    end if

    ! =============================================
    ! FIND FIRST EIGENVALUE
    ! =============================================
    call find_eigenvalue(geom_type, biot, zeta1)
    call compute_coefficient(geom_type, zeta1, coeff_c1)

    print *, '================================================='
    print *, '       EIGENVALUE ANALYSIS'
    print *, '================================================='
    print *, ''

    select case (geom_type)
        case (1)
            print *, '   Transcendental Equation:'
            print *, '   zeta * tan(zeta) = Bi'
        case (2)
            print *, '   Transcendental Equation:'
            print *, '   zeta * J1(zeta) / J0(zeta) = Bi'
        case (3)
            print *, '   Transcendental Equation:'
            print *, '   1 - zeta * cot(zeta) = Bi'
    end select

    print *, ''
    print '(A,F14.8)',    '   First eigenvalue (zeta1):  ', zeta1
    print '(A,F14.8)',    '   Coefficient (C1):          ', coeff_c1
    print *, ''

    ! =============================================
    ! TEMPERATURE AT SPECIFIED POINT AND TIME
    ! =============================================
    if (time_val > 0.0d0) then
        ! Center temperature
        theta_center = coeff_c1 * exp(-(zeta1**2) * fourier)

        ! Temperature at position
        call compute_theta_position(geom_type, zeta1, coeff_c1, fourier, x_ratio, theta_star)

        temp_result = temp_inf + (temp_init - temp_inf) * theta_star

        print *, '================================================='
        print *, '       TEMPERATURE RESULTS'
        print *, '================================================='
        print *, ''
        print *, '   One-term approximation:'
        print *, '   theta* = C1 * exp(-zeta1^2 * Fo) * f(x)'
        print *, ''
        print '(A,F12.6)',    '   theta* (center):           ', theta_center
        print '(A,F12.6)',    '   theta* (at x/L or r/r0):   ', theta_star
        print '(A,F12.2,A)',  '   T (center):                ', &
            temp_inf + (temp_init - temp_inf) * theta_center, ' deg-C'
        print '(A,F12.2,A)',  '   T (at position):           ', temp_result, ' deg-C'
        print '(A,F12.2,A)',  '   T (surface, x=1):          ', &
            temp_inf + (temp_init - temp_inf) * theta_center * &
            spatial_factor(geom_type, zeta1, 1.0d0), ' deg-C'
        print *, ''

        ! Energy transfer
        call compute_energy_ratio(geom_type, zeta1, coeff_c1, fourier, q_ratio)

        print *, '================================================='
        print *, '       ENERGY TRANSFER'
        print *, '================================================='
        print *, ''
        print '(A,F12.6)',    '   Q/Qmax:                    ', q_ratio
        print '(A,F12.2,A)',  '   Percentage:                ', q_ratio * 100.0d0, ' %'
        print *, ''
        print *, '   Q/Qmax represents the fraction of maximum'
        print *, '   possible energy transfer that has occurred'
        print *, '   by time t.'
        print *, ''
    end if

    ! =============================================
    ! TEMPERATURE PROFILE AT TIME t
    ! =============================================
    if (time_val > 0.0d0) then
        print *, '================================================='
        print *, '   TEMPERATURE PROFILE AT t =', time_val, ' s'
        print *, '================================================='
        print *, ''

        select case (geom_type)
            case (1)
                print *, '   x/L       |  T (deg-C)   |  theta*'
            case default
                print *, '   r/r0      |  T (deg-C)   |  theta*'
        end select
        print *, '   -------------------------------------------'

        n_points = 20
        do i = 0, n_points
            x_pos = dble(i) / dble(n_points)
            call compute_theta_position(geom_type, zeta1, coeff_c1, fourier, x_pos, theta_star)
            temp_x = temp_inf + (temp_init - temp_inf) * theta_star
            write(*, '(3X,F8.4,5X,A,2X,F10.3,5X,A,2X,F10.6)') &
                x_pos, '|', temp_x, '|', theta_star
        end do
        print *, ''
    end if

    ! =============================================
    ! TEMPERATURE EVOLUTION AT CENTER
    ! =============================================
    print *, '================================================='
    print *, '   TEMPERATURE EVOLUTION AT CENTER'
    print *, '================================================='
    print *, ''
    print *, '   Fo         |  t (s)       |  T_center   |  theta*'
    print *, '   --------------------------------------------------------'

    n_points = 40
    delta_fo = 5.0d0 / dble(n_points)

    do i = 0, n_points
        fo_val = dble(i) * delta_fo
        if (fo_val < 0.001d0) fo_val = 0.001d0

        theta_center = coeff_c1 * exp(-(zeta1**2) * fo_val)
        temp_x = temp_inf + (temp_init - temp_inf) * theta_center

        write(*, '(3X,F8.4,5X,A,2X,F12.2,3X,A,2X,F10.3,3X,A,2X,F10.6)') &
            fo_val, '|', fo_val * dim_char**2 / alpha_diff, &
            '|', temp_x, '|', theta_center
    end do
    print *, ''

    ! =============================================
    ! SURFACE TEMPERATURE EVOLUTION
    ! =============================================
    print *, '================================================='
    print *, '   SURFACE TEMPERATURE EVOLUTION'
    print *, '================================================='
    print *, ''
    print *, '   Fo         |  t (s)       |  T_surface  |  theta*_s'
    print *, '   --------------------------------------------------------'

    do i = 0, n_points
        fo_val = dble(i) * delta_fo
        if (fo_val < 0.001d0) fo_val = 0.001d0

        theta_center = coeff_c1 * exp(-(zeta1**2) * fo_val)
        theta_star = theta_center * spatial_factor(geom_type, zeta1, 1.0d0)
        temp_x = temp_inf + (temp_init - temp_inf) * theta_star

        write(*, '(3X,F8.4,5X,A,2X,F12.2,3X,A,2X,F10.3,3X,A,2X,F10.6)') &
            fo_val, '|', fo_val * dim_char**2 / alpha_diff, &
            '|', temp_x, '|', theta_star
    end do
    print *, ''

    ! =============================================
    ! RECOMMENDATIONS
    ! =============================================
    print *, '================================================='
    print *, '       DESIGN RECOMMENDATIONS'
    print *, '================================================='
    print *, ''

    if (biot < 0.1d0) then
        print *, '   LOW Bi (< 0.1): Internal resistance negligible'
        print *, '   - Lumped system analysis is also applicable'
        print *, '   - Temperature is nearly uniform inside'
    else if (biot < 1.0d0) then
        print *, '   MODERATE Bi (0.1 - 1.0):'
        print *, '   - Both internal and external resistances matter'
        print *, '   - One-term approximation works well for Fo > 0.2'
    else if (biot < 10.0d0) then
        print *, '   HIGH Bi (1 - 10):'
        print *, '   - Internal resistance dominates'
        print *, '   - Significant temperature gradients inside'
    else
        print *, '   VERY HIGH Bi (> 10):'
        print *, '   - Surface temperature approaches T_inf quickly'
        print *, '   - Large internal temperature gradients'
        print *, '   - Consider if surface boundary approaches'
        print *, '     a constant temperature condition'
    end if
    print *, ''

    print *, '================================================='
    print *, '          CALCULATION COMPLETE'
    print *, '================================================='
    print *, ''

contains

    ! =========================================================
    ! BESSEL FUNCTION J0 (Taylor series, 20 terms)
    ! =========================================================
    function bessel_j0(x) result(j0)
        implicit none
        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, 20 terms)
    ! =========================================================
    function bessel_j1(x) result(j1)
        implicit none
        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)
        implicit none
        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)
            case (3) ! Sphere: sin(zeta * r/r0) / (zeta * r/r0)
                if (x_r < 1.0d-10) then
                    sf = 1.0d0  ! limit as r->0
                else
                    sf = sin(zeta * x_r) / (zeta * x_r)
                end if
        end select
    end function spatial_factor

    ! =========================================================
    ! FIND FIRST EIGENVALUE via Newton-Raphson
    ! =========================================================
    subroutine find_eigenvalue(gtype, bi, zeta)
        implicit none
        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

            case (3) ! 1 - zeta * cot(zeta) = Bi
                if (bi < 0.01d0) then
                    zeta = sqrt(3.0d0 * bi)
                else if (bi < 1.0d0) then
                    zeta = sqrt(3.0d0 * bi)
                else if (bi < 10.0d0) then
                    zeta = 2.0d0 + 0.05d0 * bi
                else
                    zeta = pi - 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)

                case (3) ! f = 1 - zeta * cos(zeta)/sin(zeta) - Bi
                    if (abs(sin(zeta)) < 1.0d-12) then
                        zeta = zeta - 0.01d0
                        cycle
                    end if
                    f_val = 1.0d0 - zeta * cos(zeta) / sin(zeta) - bi
                    df_val = -cos(zeta)/sin(zeta) + zeta/(sin(zeta)**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
                case (3)
                    if (zeta >= pi) zeta = pi - 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)
        implicit none
        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)
            case (3) ! C1 = 4*(sin(zeta) - zeta*cos(zeta)) / (2*zeta - sin(2*zeta))
                c1 = 4.0d0 * (sin(zeta) - zeta * cos(zeta)) / &
                     (2.0d0 * zeta - sin(2.0d0 * zeta))
        end select
    end subroutine compute_coefficient

    ! =========================================================
    ! COMPUTE THETA* AT A POSITION
    ! =========================================================
    subroutine compute_theta_position(gtype, zeta, c1, fo, x_r, theta)
        implicit none
        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)
        implicit none
        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
            case (3) ! Q/Qmax = 1 - 3*theta0*(sin(zeta)-zeta*cos(zeta))/zeta^3
                qr = 1.0d0 - 3.0d0 * theta0 * &
                     (sin(zeta) - zeta * cos(zeta)) / (zeta**3)
        end select

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

end program transient_conduction_1d