program transient_diffusion_solid
    implicit none
    integer :: geom_type, i, n_points, n_time, iostat_val
    double precision :: L, R, Aref, D, hm, C0, Cs, time, rho_s
    double precision :: Lc, Volume, Area_s, Fo, Bi, Ccenter, Csurf, Cavg, Mabs, Mfrac
    double precision :: x, xstar, theta, Cx, relpos, tprof, Foprof, fracprof, Mprof
    character(len=80) :: geom_name
    double precision, parameter :: PI = 3.1415926535897932384626433832795d0

    read(*,*,iostat=iostat_val) geom_type
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid geometry type input.'
        stop
    end if
    read(*,*,iostat=iostat_val) L
    read(*,*,iostat=iostat_val) R
    read(*,*,iostat=iostat_val) Aref
    read(*,*,iostat=iostat_val) D
    read(*,*,iostat=iostat_val) hm
    read(*,*,iostat=iostat_val) C0
    read(*,*,iostat=iostat_val) Cs
    read(*,*,iostat=iostat_val) time
    read(*,*,iostat=iostat_val) rho_s
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all transient diffusion parameters.'
        stop
    end if

    if (D <= 0.0d0) then
        write(*,*) 'ERROR: Diffusivity D must be positive.'
        stop
    end if
    if (time < 0.0d0) then
        write(*,*) 'ERROR: Time cannot be negative.'
        stop
    end if
    if (hm < 0.0d0) then
        write(*,*) 'ERROR: Surface mass transfer coefficient cannot be negative.'
        stop
    end if

    select case (geom_type)
    case (1)
        if (L <= 0.0d0) then
            write(*,*) 'ERROR: Slab half-thickness L must be positive.'
            stop
        end if
        if (Aref <= 0.0d0) Aref = 1.0d0
        geom_name = 'Plane Wall / Slab'
        Lc = L
        Volume = 2.0d0 * L * Aref
        Area_s = 2.0d0 * Aref
    case (2)
        if (R <= 0.0d0) then
            write(*,*) 'ERROR: Cylinder radius R must be positive.'
            stop
        end if
        if (Aref <= 0.0d0) Aref = 1.0d0
        geom_name = 'Long Cylinder'
        Lc = R
        Volume = PI * R**2 * Aref
        Area_s = 2.0d0 * PI * R * Aref
    case (3)
        if (R <= 0.0d0) then
            write(*,*) 'ERROR: Sphere radius R must be positive.'
            stop
        end if
        geom_name = 'Sphere'
        Lc = R
        Volume = 4.0d0/3.0d0 * PI * R**3
        Area_s = 4.0d0 * PI * R**2
    case default
        write(*,*) 'ERROR: Invalid geometry. Use 1 slab, 2 cylinder, or 3 sphere.'
        stop
    end select

    Fo = D * time / (Lc**2)
    if (hm <= 0.0d0) then
        Bi = 0.0d0
    else
        Bi = hm * Lc / D
    end if

    Ccenter = concentration_at(geom_type, 0.0d0, Fo, Bi, C0, Cs)
    Csurf = concentration_at(geom_type, 1.0d0, Fo, Bi, C0, Cs)
    Cavg = average_concentration(geom_type, Fo, Bi, C0, Cs)
    if (abs(Cs - C0) > 1.0d-30) then
        Mfrac = (Cavg - C0) / (Cs - C0)
    else
        Mfrac = 0.0d0
    end if
    if (Mfrac < 0.0d0 .and. Cs > C0) Mfrac = 0.0d0
    if (Mfrac > 1.0d0 .and. Cs > C0) Mfrac = 1.0d0
    Mabs = (Cavg - C0) * Volume

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   TRANSIENT DIFFUSION IN SOLIDS ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A,A)')        '  Geometry Name           = ', trim(geom_name)
    write(*,'(A,I2)')       '  Geometry Type Code      = ', geom_type
    write(*,'(A,ES12.4,A)') '  Characteristic Length   = ', Lc, ' m'
    write(*,'(A,ES12.4,A)') '  Volume                  = ', Volume, ' m3'
    write(*,'(A,ES12.4,A)') '  Exposed Area            = ', Area_s, ' m2'
    write(*,*)
    write(*,'(A)') '--- DIFFUSION INPUTS ----------------------------------------'
    write(*,'(A,ES12.4,A)') '  Diffusivity (D)         = ', D, ' m2/s'
    write(*,'(A,ES12.4,A)') '  Surface h_m             = ', hm, ' m/s'
    write(*,'(A,ES12.4,A)') '  Initial Concentration   = ', C0, ' units'
    write(*,'(A,ES12.4,A)') '  Surface Equilibrium Cs  = ', Cs, ' units'
    write(*,'(A,ES12.4,A)') '  Elapsed Time            = ', time, ' s'
    write(*,'(A,ES12.4,A)') '  Solid Density           = ', rho_s, ' kg/m3'
    write(*,*)
    write(*,'(A)') '--- DIMENSIONLESS GROUPS ------------------------------------'
    write(*,'(A,ES12.4)')   '  Fourier Number (Fo)     = ', Fo
    write(*,'(A,ES12.4)')   '  Biot Number (Bi_m)      = ', Bi
    write(*,*)
    write(*,'(A)') '--- TRANSIENT DIFFUSION RESULTS -----------------------------'
    write(*,'(A,ES12.4,A)') '  Center Concentration    = ', Ccenter, ' units'
    write(*,'(A,ES12.4,A)') '  Surface Concentration   = ', Csurf, ' units'
    write(*,'(A,ES12.4,A)') '  Average Concentration   = ', Cavg, ' units'
    write(*,'(A,ES12.4,A)') '  Total Absorbed Mass     = ', Mabs, ' kg-equivalent'
    write(*,'(A,ES12.4)')   '  Mass Fraction Progress  = ', Mfrac
    if (rho_s > 0.0d0) write(*,'(A,ES12.4,A)') '  Avg Concentration ppmw  = ', Cavg / rho_s * 1.0d6, ' ppmw'
    write(*,*)

    write(*,'(A)') '--- CONCENTRATION PROFILE -----------------------------------'
    write(*,'(A)') '  x [m]       x*            theta         C(x,t)        relative_pos'
    write(*,'(A)') '  ----------------------------------------------------------'
    n_points = 61
    do i = 0, n_points-1
        xstar = dble(i) / dble(n_points-1)
        x = xstar * Lc
        theta = theta_at(geom_type, xstar, Fo, Bi)
        Cx = Cs + theta * (C0 - Cs)
        relpos = xstar
        write(*,'(F10.6,2X,F10.6,2X,ES12.4,2X,ES12.4,2X,F10.6)') x, xstar, theta, Cx, relpos
    end do
    write(*,*)

    write(*,'(A)') '--- MASS UPTAKE PROFILE -------------------------------------'
    write(*,'(A)') '  t [s]       Fo            fraction      M_abs [kg]'
    write(*,'(A)') '  ----------------------------------------------------------'
    n_time = 60
    do i = 0, n_time-1
        tprof = time * dble(i) / dble(n_time-1)
        if (time <= 0.0d0) tprof = 0.0d0
        Foprof = D * tprof / (Lc**2)
        if (abs(Cs-C0) > 1.0d-30) then
            fracprof = (average_concentration(geom_type, Foprof, Bi, C0, Cs) - C0) / (Cs - C0)
        else
            fracprof = 0.0d0
        end if
        Mprof = fracprof * (Cs - C0) * Volume
        write(*,'(F10.3,2X,ES12.4,2X,ES12.4,2X,ES12.4)') tprof, Foprof, fracprof, Mprof
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  1D Fickian diffusion in solids with Fourier-series eigenfunction expansions.'
    write(*,'(A)') '  Slab finite Bi: lambda tan(lambda) = Bi; Phi = cos(lambda x*).'
    write(*,'(A)') '  Cylinder and sphere use standard radial Fourier solutions for prescribed surface concentration when Bi is large.'
    write(*,'(A)') '  For finite Bi in cylinder/sphere, a one-term effective approximation is used for surface-resistance correction.'

contains
    double precision function theta_at(g, xs, Fo_in, Bi_in)
        implicit none
        integer, intent(in) :: g
        double precision, intent(in) :: xs, Fo_in, Bi_in
        integer :: n, maxn
        double precision :: lam, coef, sumv, root
        double precision, parameter :: PIc = 3.1415926535897932384626433832795d0
        if (Fo_in <= 0.0d0) then
            theta_at = 1.0d0
            return
        end if
        if (Bi_in <= 1.0d-12) then
            theta_at = 1.0d0
            return
        end if
        maxn = 160
        sumv = 0.0d0
        if (g == 1) then
            do n = 1, maxn
                call slab_root(n, Bi_in, root)
                coef = 4.0d0 * sin(root) / (2.0d0*root + sin(2.0d0*root))
                sumv = sumv + coef * cos(root*xs) * exp(-root*root*Fo_in)
            end do
        else if (g == 2) then
            ! Large-Bi radial cylinder approximation based on J0 zeros; finite Bi adjusted by Bi/(Bi+1).
            do n = 1, maxn
                lam = besselj0_zero_approx(n)
                coef = 2.0d0 / (lam * besselj1(lam))
                sumv = sumv + coef * besselj0(lam*xs) * exp(-lam*lam*Fo_in)
            end do
            if (Bi_in < 100.0d0) sumv = 1.0d0 - (1.0d0 - sumv) * Bi_in / (Bi_in + 1.0d0)
        else
            ! Sphere, prescribed surface concentration series.
            if (xs <= 1.0d-10) then
                do n = 1, maxn
                    lam = dble(n) * PIc
                    coef = 2.0d0 * (-1.0d0)**(n+1)
                    sumv = sumv + coef * exp(-lam*lam*Fo_in)
                end do
            else
                do n = 1, maxn
                    lam = dble(n) * PIc
                    coef = 2.0d0 * (-1.0d0)**(n+1)
                    sumv = sumv + coef * sin(lam*xs)/(lam*xs) * exp(-lam*lam*Fo_in)
                end do
            end if
            if (Bi_in < 100.0d0) sumv = 1.0d0 - (1.0d0 - sumv) * Bi_in / (Bi_in + 1.0d0)
        end if
        theta_at = max(-0.05d0, min(1.05d0, sumv))
    end function theta_at

    double precision function concentration_at(g, xs, Fo_in, Bi_in, C0_in, Cs_in)
        implicit none
        integer, intent(in) :: g
        double precision, intent(in) :: xs, Fo_in, Bi_in, C0_in, Cs_in
        concentration_at = Cs_in + theta_at(g, xs, Fo_in, Bi_in) * (C0_in - Cs_in)
    end function concentration_at

    double precision function average_concentration(g, Fo_in, Bi_in, C0_in, Cs_in)
        implicit none
        integer, intent(in) :: g
        double precision, intent(in) :: Fo_in, Bi_in
        double precision, intent(in) :: C0_in, Cs_in
        integer :: j, N
        double precision :: xs, w, integ, theta_avg
        N = 800
        integ = 0.0d0
        do j = 0, N
            xs = dble(j)/dble(N)
            if (j == 0 .or. j == N) then
                w = 0.5d0
            else
                w = 1.0d0
            end if
            if (g == 1) then
                integ = integ + w * theta_at(g, xs, Fo_in, Bi_in)
            else if (g == 2) then
                integ = integ + w * theta_at(g, xs, Fo_in, Bi_in) * 2.0d0 * xs
            else
                integ = integ + w * theta_at(g, xs, Fo_in, Bi_in) * 3.0d0 * xs**2
            end if
        end do
        theta_avg = integ / dble(N)
        average_concentration = Cs_in + theta_avg * (C0_in - Cs_in)
    end function average_concentration

    subroutine slab_root(n, Bi_in, root)
        implicit none
        integer, intent(in) :: n
        double precision, intent(in) :: Bi_in
        double precision, intent(out) :: root
        integer :: it
        double precision :: a, b, c, fa, fc, eps
        double precision, parameter :: PIc = 3.1415926535897932384626433832795d0
        eps = 1.0d-10
        a = (dble(n)-1.0d0)*PIc + eps
        b = (dble(n)-0.5d0)*PIc - eps
        if (Bi_in > 1.0d8) then
            root = (dble(n)-0.5d0)*PIc
            return
        end if
        fa = a*tan(a) - Bi_in
        do it = 1, 80
            c = 0.5d0*(a+b)
            fc = c*tan(c) - Bi_in
            if (fa*fc <= 0.0d0) then
                b = c
            else
                a = c
                fa = fc
            end if
        end do
        root = 0.5d0*(a+b)
    end subroutine slab_root

    double precision function besselj0(x)
        implicit none
        double precision, intent(in) :: x
        integer :: k
        double precision :: term, sumv
        if (abs(x) < 12.0d0) then
            term = 1.0d0; sumv = 1.0d0
            do k = 1, 80
                term = term * (-(x*x)/4.0d0) / (dble(k)*dble(k))
                sumv = sumv + term
                if (abs(term) < 1.0d-14) exit
            end do
            besselj0 = sumv
        else
            besselj0 = sqrt(2.0d0/(3.141592653589793d0*x)) * cos(x - 0.25d0*3.141592653589793d0)
        end if
    end function besselj0

    double precision function besselj1(x)
        implicit none
        double precision, intent(in) :: x
        integer :: k
        double precision :: term, sumv
        if (abs(x) < 12.0d0) then
            term = x/2.0d0; sumv = term
            do k = 1, 80
                term = term * (-(x*x)/4.0d0) / (dble(k)*dble(k+1))
                sumv = sumv + term
                if (abs(term) < 1.0d-14) exit
            end do
            besselj1 = sumv
        else
            besselj1 = sqrt(2.0d0/(3.141592653589793d0*x)) * cos(x - 0.75d0*3.141592653589793d0)
        end if
    end function besselj1

    double precision function besselj0_zero_approx(n)
        implicit none
        integer, intent(in) :: n
        double precision, parameter :: PIc = 3.1415926535897932384626433832795d0
        ! McMahon approximation; adequate for profile visualization and series convergence.
        besselj0_zero_approx = (dble(n)-0.25d0)*PIc + 1.0d0/(8.0d0*(dble(n)-0.25d0)*PIc)
    end function besselj0_zero_approx
end program transient_diffusion_solid
