program Annular_Fin_Efficiency
    implicit none

    ! ----------------------------------------------------------------
    !  ANNULAR (CIRCULAR) FIN EFFICIENCY CALCULATOR
    !
    !  Correct formula from Incropera et al., "Fundamentals of Heat
    !  and Mass Transfer", 7th ed., Table 3.5, Case 4:
    !
    !            2 * r1         K1(m*r1) * I1(m*r2c) - I1(m*r1) * K1(m*r2c)
    !  eta_f = ─────────── * ──────────────────────────────────────────────────
    !          m*(r2c²-r1²)    I0(m*r1) * K1(m*r2c) + K0(m*r1) * I1(m*r2c)
    !
    !  where r2c = r2 + t/2  (corrected radius for convective tip)
    !        m   = sqrt(2h / k*t)
    !
    !  Modified Bessel functions computed via Abramowitz & Stegun §9.8
    ! ----------------------------------------------------------------

    real(8) :: r1, t, r2, r2c, k, h, T0, Tinf
    real(8) :: L, m, mR1, mR2c
    real(8) :: I0_r1, I1_r1, K0_r1, K1_r1
    real(8) :: I0_r2c, I1_r2c, K0_r2c, K1_r2c
    real(8) :: num, den, eta
    real(8) :: Asurf, Qmax, Q_fin, theta_base, epsilon_fin, Biot_fin
    real(8) :: C1, C2, r_i, mr_i, I0_ri, K0_ri, theta_r, T_r
    real(8), parameter :: PI = 3.14159265358979323846d0
    integer :: i, n_points

    write(*,*) '================================================================'
    write(*,*) '      ANNULAR (CIRCULAR) FIN EFFICIENCY CALCULATOR'
    write(*,*) '================================================================'
    write(*,*)
    write(*,*) 'Reading input parameters...'

    ! Inputs
    read(*,*) r1     ! Inner radius [m]
    read(*,*) t      ! Fin thickness [m]
    read(*,*) r2c    ! Corrected tip radius r2c = r2 + t/2 [m]
    read(*,*) k      ! Thermal conductivity [W/m.K]
    read(*,*) h      ! Convection coefficient [W/m2.K]
    read(*,*) T0     ! Base temperature [deg-C]
    read(*,*) Tinf   ! Ambient temperature [deg-C]

    r2 = r2c - t / 2.0d0
    L  = r2 - r1

    write(*,*)
    write(*,*) '================================================================'
    write(*,*) '                    INPUT PARAMETERS'
    write(*,*) '================================================================'
    write(*,*)
    write(*,'(A,F10.4,A)') ' Inner radius          (r1):    ', r1*1000.0d0,  ' mm'
    write(*,'(A,F10.4,A)') ' Outer radius          (r2):    ', r2*1000.0d0,  ' mm'
    write(*,'(A,F10.4,A)') ' Corrected tip radius (r2c):    ', r2c*1000.0d0, ' mm'
    write(*,'(A,F10.4,A)') ' Fin thickness          (t):    ', t*1000.0d0,   ' mm'
    write(*,'(A,F10.4,A)') ' Fin radial length (L=r2-r1):   ', L*1000.0d0,   ' mm'
    write(*,'(A,F10.3,A)') ' Thermal Conductivity   (k):    ', k,            ' W/m.K'
    write(*,'(A,F10.3,A)') ' Convection Coefficient (h):    ', h,            ' W/m2.K'
    write(*,'(A,F10.2,A)') ' Base Temperature       (T0):   ', T0,           ' deg-C'
    write(*,'(A,F10.2,A)') ' Ambient Temperature  (Tinf):   ', Tinf,         ' deg-C'
    write(*,*)

    ! Fin parameter
    m    = sqrt(2.0d0 * h / (k * t))
    mR1  = m * r1
    mR2c = m * r2c

    write(*,*) '================================================================'
    write(*,*) '                   FIN PARAMETERS'
    write(*,*) '================================================================'
    write(*,*)
    write(*,'(A,F12.4,A)') ' Fin parameter        (m):      ', m,        ' 1/m'
    write(*,'(A,F12.6  )') ' m * r1:                        ', mR1
    write(*,'(A,F12.6  )') ' m * r2c (corrected):           ', mR2c
    write(*,'(A,F12.6  )') ' Radius ratio (r2/r1):          ', r2 / r1
    write(*,'(A,F12.6  )') ' Corrected ratio (r2c/r1):      ', r2c / r1
    write(*,*)

    ! Compute Bessel functions
    call bessel_I0(mR1,  I0_r1)
    call bessel_I1(mR1,  I1_r1)
    call bessel_K0(mR1,  K0_r1)
    call bessel_K1(mR1,  K1_r1)
    call bessel_I0(mR2c, I0_r2c)
    call bessel_I1(mR2c, I1_r2c)
    call bessel_K0(mR2c, K0_r2c)
    call bessel_K1(mR2c, K1_r2c)

    write(*,*) '================================================================'
    write(*,*) '              BESSEL FUNCTION VALUES'
    write(*,*) '================================================================'
    write(*,*)
    write(*,'(A,F14.8)') ' I0(m*r1)  = ', I0_r1
    write(*,'(A,F14.8)') ' I1(m*r1)  = ', I1_r1
    write(*,'(A,F14.8)') ' K0(m*r1)  = ', K0_r1
    write(*,'(A,F14.8)') ' K1(m*r1)  = ', K1_r1
    write(*,*)
    write(*,'(A,F14.8)') ' I0(m*r2c) = ', I0_r2c
    write(*,'(A,F14.8)') ' I1(m*r2c) = ', I1_r2c
    write(*,'(A,F14.8)') ' K0(m*r2c) = ', K0_r2c
    write(*,'(A,F14.8)') ' K1(m*r2c) = ', K1_r2c
    write(*,*)

    ! ----------------------------------------------------------------
    ! Fin efficiency (Incropera Table 3.5, Case 4):
    !
    !           2 r1          K1(mr1)*I1(mr2c) - I1(mr1)*K1(mr2c)
    ! eta_f = ───────── × ──────────────────────────────────────────
    !         m(r2c²-r1²)   I0(mr1)*K1(mr2c) + K0(mr1)*I1(mr2c)
    ! ----------------------------------------------------------------
    num = K1_r1 * I1_r2c - I1_r1 * K1_r2c
    den = I0_r1 * K1_r2c + K0_r1 * I1_r2c

    if (abs(den) < 1.0d-15) then
        write(*,*) 'ERROR: Denominator near zero — check geometry inputs.'
        stop
    end if

    eta = (2.0d0 * r1 / (m * (r2c**2 - r1**2))) * (num / den)
    if (eta > 1.0d0) eta = 1.0d0
    if (eta < 0.0d0) eta = 0.0d0

    ! Performance
    theta_base  = T0 - Tinf
    Asurf       = 2.0d0 * PI * (r2c**2 - r1**2)
    Qmax        = h * Asurf * theta_base
    Q_fin       = eta * Qmax
    epsilon_fin = Q_fin / max(h * 2.0d0 * PI * r1 * t * theta_base, 1.0d-30)
    Biot_fin    = h * t / (2.0d0 * k)

    write(*,*) '================================================================'
    write(*,*) '                  PERFORMANCE RESULTS'
    write(*,*) '================================================================'
    write(*,*)
    write(*,'(A,F12.4,A)') ' Fin efficiency         (eta):  ', eta * 100.0d0,     ' %'
    write(*,'(A,F12.6,A)') ' Fin corrected area (A_fin):    ', Asurf * 1.0d4,     ' cm2'
    write(*,'(A,F12.4,A)') ' Max heat transfer   (Q_max):   ', Qmax,              ' W'
    write(*,'(A,F12.4,A)') ' Actual heat transfer (Q_fin):  ', Q_fin,             ' W'
    write(*,'(A,F12.4  )') ' Fin effectiveness (epsilon):   ', epsilon_fin
    write(*,'(A,F12.6  )') ' Biot number (Bi = h*t/2k):    ', Biot_fin
    write(*,*)

    ! Validity
    write(*,*) '================================================================'
    write(*,*) '                 VALIDITY ASSESSMENT'
    write(*,*) '================================================================'
    write(*,*)
    if (Biot_fin < 0.1d0) then
        write(*,'(A,F8.5,A)') ' [OK]    Bi = ', Biot_fin, ' < 0.1  →  1-D fin assumption VALID.'
    else
        write(*,'(A,F8.5,A)') ' [WARN]  Bi = ', Biot_fin, ' >= 0.1 →  1-D introduces error; use 2-D analysis.'
    end if
    write(*,*)
    if (eta > 0.9d0) then
        write(*,*) ' [EXCELLENT]  eta > 90% — highly efficient fin design.'
    else if (eta >= 0.7d0) then
        write(*,*) ' [GOOD]       eta 70–90% — acceptable performance.'
    else if (eta >= 0.5d0) then
        write(*,*) ' [MARGINAL]   eta 50–70% — consider increasing k or reducing L.'
    else
        write(*,*) ' [POOR]       eta < 50%  — fin is under-performing. Redesign recommended.'
    end if
    write(*,*)
    if (epsilon_fin < 2.0d0) then
        write(*,*) ' [WARNING]    Effectiveness < 2 — fin usage may not be justified.'
    else if (epsilon_fin > 10.0d0) then
        write(*,*) ' [EXCELLENT]  Effectiveness > 10 — fin greatly enhances heat transfer.'
    else
        write(*,*) ' [GOOD]       Effectiveness in range 2–10 — fin is beneficial.'
    end if
    write(*,*)

    ! Temperature profile
    ! theta(r)/theta_b = C1*I0(mr) + C2*K0(mr)
    ! BC1: theta(r1) = theta_b  =>  C1*I0(mr1) + C2*K0(mr1) = 1
    ! BC2: dT/dr|r2c = 0        =>  C1*I1(mr2c) - C2*K1(mr2c) = 0
    ! => C2 = C1 * I1(mr2c) / K1(mr2c)
    ! => C1 = K1(mr2c) / [I0(mr1)*K1(mr2c) + I1(mr2c)*K0(mr1)]
    C1 = K1_r2c / (I0_r1 * K1_r2c + I1_r2c * K0_r1)
    C2 = C1 * I1_r2c / K1_r2c

    write(*,*) '================================================================'
    write(*,*) '          RADIAL TEMPERATURE DISTRIBUTION'
    write(*,*) '================================================================'
    write(*,*)
    write(*,*) '  r [mm]    |  T [deg-C]    |  theta/theta_b'
    write(*,*) ' ---------------------------------------------------'

    n_points = 20
    do i = 0, n_points
        r_i     = r1 + (real(i, 8) / real(n_points, 8)) * (r2 - r1)
        mr_i    = m * r_i
        call bessel_I0(mr_i, I0_ri)
        call bessel_K0(mr_i, K0_ri)
        theta_r = C1 * I0_ri + C2 * K0_ri
        T_r     = Tinf + theta_r * theta_base
        write(*,'(2X,F10.3,2X,A,2X,F10.4,4X,A,2X,F8.5)') &
            r_i*1000.0d0, '|', T_r, '|', theta_r
    end do
    write(*,*)

    ! Material comparison
    write(*,*) '================================================================'
    write(*,*) '                MATERIAL COMPARISON'
    write(*,*) '================================================================'
    write(*,*)
    write(*,*) '  (Same geometry, h, T0, Tinf — material conductivity varies)'
    write(*,*)
    call material_eta('Aluminum (6061-T6) ', 200.0d0, m, k, r1, r2c)
    call material_eta('Copper (pure)      ', 385.0d0, m, k, r1, r2c)
    call material_eta('Steel (low-carbon) ',  50.0d0, m, k, r1, r2c)
    call material_eta('Brass (Cu65/Zn35)  ', 110.0d0, m, k, r1, r2c)
    call material_eta('Cast Iron          ',  52.0d0, m, k, r1, r2c)
    write(*,*)

    write(*,*) '================================================================'
    write(*,*) '               CALCULATION COMPLETE'
    write(*,*) '================================================================'

contains

    ! Abramowitz & Stegun §9.8.1–9.8.8 polynomial approximations

    subroutine bessel_I0(x, res)
        implicit none
        real(8), intent(in)  :: x
        real(8), intent(out) :: res
        real(8) :: tx, p
        if (x <= 3.75d0) then
            tx = (x / 3.75d0)**2
            res = 1.0d0 + tx*(3.5156229d0 + tx*(3.0899424d0 + tx*(1.2067492d0 &
                + tx*(0.2659732d0 + tx*(0.0360768d0 + tx*0.0045813d0)))))
        else
            tx = 3.75d0 / x
            p = 0.39894228d0 + tx*(0.01328592d0 + tx*(0.00225319d0 &
                + tx*(-0.00157565d0 + tx*(0.00916281d0 + tx*(-0.02057706d0 &
                + tx*(0.02635537d0 + tx*(-0.01647633d0 + tx*0.00392377d0)))))))
            res = (exp(x) / sqrt(x)) * p
        end if
    end subroutine bessel_I0

    subroutine bessel_I1(x, res)
        implicit none
        real(8), intent(in)  :: x
        real(8), intent(out) :: res
        real(8) :: tx, p
        if (x <= 3.75d0) then
            tx = (x / 3.75d0)**2
            res = x * (0.5d0 + tx*(0.87890594d0 + tx*(0.51498869d0 &
                + tx*(0.15084934d0 + tx*(0.02658733d0 + tx*(0.00301532d0 &
                + tx*0.00032411d0))))))
        else
            tx = 3.75d0 / x
            p = 0.39894228d0 + tx*(-0.03988024d0 + tx*(-0.00362018d0 &
                + tx*(0.00163801d0 + tx*(-0.01031555d0 + tx*(0.02282967d0 &
                + tx*(-0.02895312d0 + tx*(0.01787654d0 - tx*0.00420059d0)))))))
            res = (exp(x) / sqrt(x)) * p
        end if
    end subroutine bessel_I1

    subroutine bessel_K0(x, res)
        implicit none
        real(8), intent(in)  :: x
        real(8), intent(out) :: res
        real(8) :: tx, p, I0_val
        if (x <= 2.0d0) then
            tx = (x / 2.0d0)**2
            call bessel_I0(x, I0_val)
            p = -0.57721566d0 + tx*(0.42278420d0 + tx*(0.23069756d0 &
                + tx*(0.03488590d0 + tx*(0.00262698d0 + tx*(0.00010750d0 &
                + tx*0.0000074d0)))))
            res = -log(x / 2.0d0) * I0_val + p
        else
            tx = 2.0d0 / x
            p = 1.25331414d0 + tx*(-0.07832358d0 + tx*(0.02189568d0 &
                + tx*(-0.01062446d0 + tx*(0.00587872d0 + tx*(-0.00251540d0 &
                + tx*0.00053208d0)))))
            res = (exp(-x) / sqrt(x)) * p
        end if
    end subroutine bessel_K0

    subroutine bessel_K1(x, res)
        implicit none
        real(8), intent(in)  :: x
        real(8), intent(out) :: res
        real(8) :: tx, p, I1_val
        if (x <= 2.0d0) then
            tx = (x / 2.0d0)**2
            call bessel_I1(x, I1_val)
            p = 1.0d0 + tx*(0.15443144d0 + tx*(-0.67278579d0 &
                + tx*(-0.18156897d0 + tx*(-0.01919402d0 + tx*(-0.00110404d0 &
                - tx*0.00004686d0)))))
            res = log(x / 2.0d0) * I1_val + p / x
        else
            tx = 2.0d0 / x
            p = 1.25331414d0 + tx*(0.23498619d0 + tx*(-0.03655620d0 &
                + tx*(0.01504268d0 + tx*(-0.00780353d0 + tx*(0.00325614d0 &
                - tx*0.00068245d0)))))
            res = (exp(-x) / sqrt(x)) * p
        end if
    end subroutine bessel_K1

    subroutine material_eta(name, k_mat, m_ref, k_ref, r1_in, r2c_in)
        implicit none
        character(len=*), intent(in) :: name
        real(8), intent(in) :: k_mat, m_ref, k_ref, r1_in, r2c_in
        real(8) :: m_mat, mR1m, mR2cm
        real(8) :: I0r1, I1r1, K0r1, K1r1
        real(8) :: I0r2c, I1r2c, K0r2c, K1r2c
        real(8) :: num_m, den_m, eta_m
        ! Scale m: m = sqrt(2h/kt), so m_mat = m_ref * sqrt(k_ref / k_mat)
        m_mat  = m_ref * sqrt(k_ref / k_mat)
        mR1m   = m_mat * r1_in
        mR2cm  = m_mat * r2c_in
        call bessel_I0(mR1m,  I0r1)
        call bessel_I1(mR1m,  I1r1)
        call bessel_K0(mR1m,  K0r1)
        call bessel_K1(mR1m,  K1r1)
        call bessel_I0(mR2cm, I0r2c)
        call bessel_I1(mR2cm, I1r2c)
        call bessel_K0(mR2cm, K0r2c)
        call bessel_K1(mR2cm, K1r2c)
        num_m = K1r1 * I1r2c - I1r1 * K1r2c
        den_m = I0r1 * K1r2c + K0r1 * I1r2c
        if (abs(den_m) < 1.0d-15) then
            write(*,'(3X,A20,A)') name, ':  near-singular — skipped'
            return
        end if
        eta_m = (2.0d0 * r1_in / (m_mat * (r2c_in**2 - r1_in**2))) * (num_m / den_m)
        if (eta_m > 1.0d0) eta_m = 1.0d0
        if (eta_m < 0.0d0) eta_m = 0.0d0
        write(*,'(3X,A20,A,F7.2,A,F7.3,A)') name, ':  k = ', k_mat, &
            ' W/m.K   →   eta = ', eta_m * 100.0d0, ' %'
    end subroutine material_eta

end program Annular_Fin_Efficiency