program diffusion_reaction
    implicit none
    integer :: geom_type, reaction_order, i, n_points, iostat_val
    double precision :: Rpellet, Lslab, De, k, Cs, rhoCat, targetPhi
    double precision :: Lc, phi, eta, wp, rateIntrinsic, rateObs, phi_i, eta_slab_i, eta_cyl_i, eta_sph_i, wp_i
    character(len=80) :: geomName, regime

    read(*,*,iostat=iostat_val) geom_type
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid geometry input.'
        stop
    end if
    read(*,*,iostat=iostat_val) reaction_order
    read(*,*,iostat=iostat_val) Rpellet
    read(*,*,iostat=iostat_val) Lslab
    read(*,*,iostat=iostat_val) De
    read(*,*,iostat=iostat_val) k
    read(*,*,iostat=iostat_val) Cs
    read(*,*,iostat=iostat_val) rhoCat
    read(*,*,iostat=iostat_val) targetPhi
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all diffusion-reaction inputs.'
        stop
    end if
    if (De <= 0.0d0 .or. k < 0.0d0 .or. Cs <= 0.0d0) then
        write(*,*) 'ERROR: De and Cs must be positive; k must be non-negative.'
        stop
    end if

    select case(geom_type)
    case(1)
        if (Lslab <= 0.0d0) then
            write(*,*) 'ERROR: Slab half-thickness must be positive.'
            stop
        end if
        geomName = 'Slab / Flat Pellet'
        Lc = Lslab
    case(2)
        if (Rpellet <= 0.0d0) then
            write(*,*) 'ERROR: Cylinder radius must be positive.'
            stop
        end if
        geomName = 'Long Cylinder'
        Lc = Rpellet
    case(3)
        if (Rpellet <= 0.0d0) then
            write(*,*) 'ERROR: Sphere radius must be positive.'
            stop
        end if
        geomName = 'Sphere'
        Lc = Rpellet
    case default
        write(*,*) 'ERROR: Geometry must be 1 slab, 2 cylinder, or 3 sphere.'
        stop
    end select

    if (reaction_order == 1) then
        phi = Lc*sqrt(k/De)
        rateIntrinsic = k*Cs
    else
        phi = Lc*sqrt(max(k,0.0d0)/(De*Cs))
        rateIntrinsic = k
    end if
    eta = effectiveness(geom_type, phi)
    rateObs = eta*rateIntrinsic
    wp = rateObs*Lc**2/(Cs*De)

    if (phi < 0.3d0 .and. wp < 0.3d0) then
        regime = 'Kinetic-control likely; internal diffusion small'
    else if (phi < 3.0d0 .or. wp < 1.0d0) then
        regime = 'Mixed kinetic/diffusion regime'
    else
        regime = 'Strong internal diffusion limitation likely'
    end if

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   DIFFUSION WITH CHEMICAL REACTION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A,A)')        '  Geometry Name            = ', trim(geomName)
    write(*,'(A,I8)')       '  Reaction Order Model     = ', reaction_order
    write(*,'(A,ES12.4,A)') '  Characteristic Length    = ', Lc, ' m'
    write(*,'(A,ES12.4,A)') '  Effective Diffusivity    = ', De, ' m2/s'
    write(*,'(A,ES12.4)')   '  Rate Constant k          = ', k
    write(*,'(A,ES12.4,A)') '  Surface Concentration Cs = ', Cs, ' mol/m3'
    write(*,*)
    write(*,'(A)') '--- REACTION-DIFFUSION RESULTS -----------------------------'
    write(*,'(A,ES12.4)')   '  Thiele Modulus (phi)     = ', phi
    write(*,'(A,ES12.4)')   '  Effectiveness Factor (eta) = ', eta
    write(*,'(A,ES12.4)')   '  Intrinsic Surface Rate   = ', rateIntrinsic
    write(*,'(A,ES12.4)')   '  Observable Rate          = ', rateObs
    write(*,'(A,ES12.4)')   '  Weisz-Prater Criterion   = ', wp
    write(*,'(A,A)')        '  Regime Assessment        = ', trim(regime)
    write(*,*)
    write(*,'(A)') '--- EFFECTIVENESS PROFILE -----------------------------------'
    write(*,'(A)') '  phi          eta_slab      eta_cylinder  eta_sphere    WP_if_sphere'
    write(*,'(A)') '  -------------------------------------------------------------------'
    n_points = 80
    do i=0,n_points-1
        phi_i = 0.01d0 * (1000.0d0)**(dble(i)/dble(n_points-1))
        eta_slab_i = effectiveness(1, phi_i)
        eta_cyl_i = effectiveness(2, phi_i)
        eta_sph_i = effectiveness(3, phi_i)
        wp_i = eta_sph_i*phi_i**2
        write(*,'(F10.5,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
            phi_i, eta_slab_i, eta_cyl_i, eta_sph_i, wp_i
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  First-order Thiele modulus: phi = Lc*sqrt(k/De).'
    write(*,'(A)') '  Slab eta = tanh(phi)/phi.'
    write(*,'(A)') '  Sphere eta = 3/phi^2*(phi*coth(phi)-1).'
    write(*,'(A)') '  Cylinder eta approximated by interpolation between slab and sphere behavior.'
    write(*,'(A)') '  Observable Weisz-Prater: C_WP = r_obs*Lc^2/(Cs*De).'

contains
    double precision function effectiveness(g, ph)
        implicit none
        integer, intent(in) :: g
        double precision, intent(in) :: ph
        double precision :: cothv, slab, sph
        if (ph < 1.0d-8) then
            effectiveness = 1.0d0
            return
        end if
        slab = tanh(ph)/ph
        if (ph > 50.0d0) then
            cothv = 1.0d0
        else
            cothv = cosh(ph)/sinh(ph)
        end if
        sph = 3.0d0/(ph*ph)*(ph*cothv - 1.0d0)
        if (g == 1) then
            effectiveness = slab
        else if (g == 3) then
            effectiveness = sph
        else
            effectiveness = 0.5d0*(slab+sph)
        end if
        if (effectiveness > 1.0d0) effectiveness = 1.0d0
        if (effectiveness < 0.0d0) effectiveness = 0.0d0
    end function effectiveness
end program diffusion_reaction