Bethe-Bloch

bethebloch.jl

module BethebBloch
import Unitful: Quantity, 𝐋, 𝐌, 𝐓, @u_str
export ∂ᵨₓE, ∂ₓE, calcrange, energyloss

uLength = Quantity{<:Number,𝐋}
uEnergy = Quantity{<:Number,𝐋^2 * 𝐌 * 𝐓^(-2)}
uDensity = Quantity{<:Number,𝐌 * 𝐋^(-3)}
uStoppingPower = Quantity{<:Number,𝐋^4 * 𝐓^(-2)}

mₚ = 938.272u"MeV/c^2"
mₙ = 939.565u"MeV/c^2"

# Default step of energy ΔE
ΔE₀ = 0.5u"MeV"
# Default step of range Δx
Δx₀ = 1.0u"mm"
# limitation of Bethe-Bloch eq.
Tₗᵢₘ = 2.0u"MeV"

"""
    ∂ᵨₓE(T::Quantity,z,a,Z,A)
    ∂ᵨₓE(T::Number,z,a,Z,A)

Calculate the stopping power `dE/d(ρx)` of charged particles using Bethe-Bloch Formula.

# Arguments
- `T`: Kinetic energy of incident particle.
- `z`: charge of incident particle in units of e.
- `a`: Atomic weight of incident particle.
- `Z`: Atomic number of absorbing material.
- `A`: Atomic weight of absorbing material.

# Examples
```julia-repl
julia> using Unitful

julia> ∂ᵨₓE(100u"MeV", 1, 1, 13, 27)
5.6938828026926505 cm² MeV g⁻¹

julia> ∂ᵨₓE(400, 2, 4, 13, 27)
22.788468359490793 cm² MeV g⁻¹
```
"""
function ∂ᵨₓE(T::uEnergy, z, a, Z, A)
    # Check T > limitation of Bethe-Bloch Tₗᵢₘ
    if T < Tₗᵢₘ
        return NaN
    end

    # Constant part: 2πNₐrₑ²mₑc²
    rₑ = 2.817e-13u"cm"
    mₑ = 0.510998u"MeV/c^2"
    mₑc² = mₑ * 1u"c^2"
    Nₐ = 6.02214076e23u"g^(-1)"
    Const = 2π * Nₐ * rₑ^2 * mₑc²

    # Calculate Mass
    n = a - z
    m = z * mₚ + n * mₙ
    mc² = m * 1u"c^2"
    
    # E = T + mc² = γmc²
    γ = T / mc² + 1
    β = √(1 - 1 / γ^2)

    # Wₘₐₓ: maximum energy transfer in a single collision
    η = β * γ
    s = mₑ / m
    Wₘₐₓ = 2mₑc² * η^2 / (1 + 2s * √(1 + η^2) + s^2)

    # I: mean excitation potential
    if Z < 13
        I = (12Z + 7) * 1e-6u"MeV"
        else
        I = (9.76 + 58.8 * Z^(-1.19)) * Z * 1e-6u"MeV"
    end
    return Const * Z / A * z^2 / β^2 * (
        log((2mₑc² * γ^2 * β^2 * Wₘₐₓ) / I^2)
        - 2β^2
)
end

function ∂ᵨₓE(T::Number, z, a, Z, A)
    return ∂ᵨₓE(T * 1u"MeV", z, a, Z, A)
end

# Calculate dE/dx
"""
    ∂ₓE(stoppingpower::Quantity, ρ::Quantity)
    ∂ₓE(stoppingpower::Quantity, ρ::Number)
    ∂ₓE(stoppingpower::Number, ρ::Quantity)
    ∂ₓE(stoppingpower::Number, ρ::Number)

Calculate the mean energy loss per distance travelled `dE/dx` of charged particles using Bethe-Bloch Formula.

# Arguments
- `stoppingpower`: stoppingpower `dE/d(ρx)`.
- `ρ`: Density of absorbing material.

# Examples
```julia-repl
julia> using Unitful

julia> ∂ₓE(2.259e+1u"MeV*cm^2/g", 2.7u"g/cm^3")
60.993 MeV cm⁻¹

julia> ∂ₓE(2.259e+1u"MeV*cm^2/g", 2.7)
60.993 MeV cm⁻¹

julia> ∂ₓE(2.259e+1, 2.7u"g/cm^3")
60.993 MeV cm⁻¹

julia> ∂ₓE(2.259e+1, 2.7)
60.993 MeV cm⁻¹
```
"""
function ∂ₓE(stoppingpower::uStoppingPower, ρ::uDensity)
    return stoppingpower * ρ
end

function ∂ₓE(stoppingpower::uStoppingPower, ρ::Number)
    return stoppingpower * ρ * 1u"g/cm^3"
end

function ∂ₓE(stoppingpower::Number, ρ::uDensity)
    return stoppingpower * 1u"MeV*cm^2/g" * ρ
end

function ∂ₓE(stoppingpower::Number, ρ::Number)
    return stoppingpower * 1u"MeV*cm^2/g" * ρ * 1u"g/cm^3"
end

"""
    ∂ₓE(T,z,a,Z,A,ρ)

Calculate the mean energy loss per distance travelled `dE/dx` of charged particles using Bethe-Bloch Formula.

# Arguments
- `T`: Kinetic energy of incident particle.
- `z`: charge of incident particle in units of e.
- `a`: Atomic weight of incident particle.
- `Z`: Atomic number of absorbing material.
- `A`: Atomic weight of absorbing material.
- `ρ`: Density of absorbing material.

# Examples
```julia-repl
julia> using Unitful

julia> ∂ₓE(100u"MeV", 1, 1, 13, 27, 2.7u"g/cm^3")
15.373483567270156 MeV cm⁻¹

julia> ∂ₓE(100u"MeV", 1, 1, 13, 27, 2.7)
15.373483567270156 MeV cm⁻¹

julia> ∂ₓE(100, 1, 1, 13, 27, 2.7)
15.373483567270156 MeV cm⁻¹

julia> ∂ₓE(400u"MeV", 2, 4, 13, 27, 2.7u"g/cm^3")
61.528864570625146 MeV cm⁻¹

julia> ∂ₓE(400, 2, 4, 13, 27, 2.7u"g/cm^3")
61.528864570625146 MeV cm⁻¹
```
"""
function ∂ₓE(T, z, a, Z, A, ρ)
    stoppingpower = ∂ᵨₓE(T, z, a, Z, A)
    return ∂ₓE(stoppingpower, ρ)
end

"""
    calcrange(T::Quantity, z, a, Z, A, ρ; ΔE::Quantity)
    calcrange(T::Number, z, a, Z, A, ρ; ΔE::Quantity)

Calculate the mean range `∫(dx/dE)dE` of charged particles.

# Arguments
- `T`: Kinetic energy of incident particle.
- `z`: charge of incident particle in units of e.
- `a`: Atomic weight of incident particle.
- `Z`: Atomic number of absorbing material.
- `A`: Atomic weight of absorbing material.
- `ρ`: Density of absorbing material.
- `ΔE`: Step of energy. Default `ΔE=$ΔE₀`.

# Examples
```julia-repl
julia> using Unitful

julia> calcrange(100u"MeV", 1, 1, 13, 27, 2.7u"g/cm^3")
0.036842499703887754 m

julia> calcrange(100u"MeV", 1, 1, 13, 27, 2.7u"g/cm^3") |> u"cm"
3.6842499703887754 cm

julia> calcrange(400u"MeV", 2, 4, 13, 27, 2.7u"g/cm^3") |> u"cm"
3.6856078193822186 cm

julia> calcrange(400u"MeV", 2, 4, 13, 27, 2.7) |> u"cm"
3.6856078193822186 cm

julia> calcrange(400, 2, 4, 13, 27, 2.7) |> u"cm"
3.6856078193822186 cm

julia> calcrange(400u"MeV", 2, 4, 13, 27, 2.7u"g/cm^3", ΔE=0.1u"MeV") |> u"cm"
3.6856088872357953 cm
```
"""
function calcrange(T::uEnergy, z, a, Z, A, ρ; ΔE::uEnergy=ΔE₀)
    r = 0u"mm"
    while T - Tₗᵢₘ > ΔE
        # ∫(T→Tₗᵢₘ+ϵ) (dx/dE)dE
        ∂ₓE₀ = ∂ₓE(T, z, a, Z, A, ρ)
        ∂ₓE₁ = ∂ₓE(T - ΔE, z, a, Z, A, ρ)
        r += ΔE * (1 / ∂ₓE₀ + 1 / ∂ₓE₁) / 2
        T -= ΔE
    end
    # Tₗᵢₘ + ΔE > T > Tₗᵢₘ
    # ∫(Tₗᵢₘ+ϵ→Tₗᵢₘ) (dx/dE)dE
    ∂ₓE₀ = ∂ₓE(T, z, a, Z, A, ρ)
    ∂ₓE₁ = ∂ₓE(Tₗᵢₘ, z, a, Z, A, ρ)
    r += (T - Tₗᵢₘ) * (1 / ∂ₓE₀ + 1 / ∂ₓE₁) / 2
    return r
end

function calcrange(T::Number, z, a, Z, A, ρ; ΔE::uEnergy=ΔE₀)
    return calcrange(T * 1u"MeV", z, a, Z, A, ρ, ΔE=ΔE)
end

"""
    energyloss(T::Quantity, z, a, Z, A, ρ, depth::Quantity; Δx::Quantity)
    energyloss(T::Quantity, z, a, Z, A, ρ, depth::Number; Δx::Quantity)
    energyloss(T::Number, z, a, Z, A, ρ, depth::Quantity; Δx::Quantity)
    energyloss(T::Number, z, a, Z, A, ρ, depth::Number; Δx::Quantity)

Calculate the total energy loss in absorbing materials `∫(dE/dx)dx`.

# Arguments
- `T`: Kinetic energy of incident particle.
- `z`: charge of incident particle in units of e.
- `a`: Atomic weight of incident particle.
- `Z`: Atomic number of absorbing material.
- `A`: Atomic weight of absorbing material.
- `ρ`: Density of absorbing material.
- `depth`: Length of absorbing material.
- `Δx`: Step of range. Default `Δx=$Δx₀`.

# Examples
```julia-repl
julia> using Unitful

julia> energyloss(100u"MeV", 1, 1, 13, 27, 2.7u"g/cm^3", 2u"cm")
35.63318223527754 MeV

julia> energyloss(100u"MeV", 1, 1, 13, 27, 2.7u"g/cm^3", 5u"cm")
Stop at x = 3.687475418620985 cm
99.99999999999999 MeV

julia> energyloss(400u"MeV", 2, 4, 13, 27, 2.7u"g/cm^3", 2u"cm")
142.63206659101576 MeV

julia> energyloss(400u"MeV", 2, 4, 13, 27, 2.7u"g/cm^3", 5u"cm")
Stop at x = 3.68609556130568 cm
399.99999999999983 MeV

julia> energyloss(400, 2, 4, 13, 27, 2.7, 5)
Stop at x = 3.68609556130568 cm
399.99999999999983 MeV

julia> energyloss(400u"MeV", 2, 4, 13, 27, 2.7u"g/cm^3", 2u"cm", Δx=0.1u"mm")
142.63033203631417 MeV
```
"""
function energyloss(T::uEnergy, z, a, Z, A, ρ, depth::uLength; Δx::uLength=Δx₀)
    Eloss = 0u"MeV"
    x = depth
    while x > Δx
        r = calcrange(T, z, a, Z, A, ρ)
        if r > Δx
            ∂ₓE₀ = ∂ₓE(T, z, a, Z, A, ρ)
            ∂ₓE₁ = ∂ₓE(T - ∂ₓE₀ * Δx, z, a, Z, A, ρ)
            ΔT = Δx * (∂ₓE₀ + ∂ₓE₁) / 2
            Eloss += ΔT
            T -= ΔT
            x -= Δx
        else
            # stop
            x -= r
            println("Stop at x = $(depth - x |> u"cm")")
            Eloss += T
            return Eloss |> u"MeV"
        end
    end
    # x < Δx
    r = calcrange(T, z, a, Z, A, ρ)
    if r > x
        ∂ₓE₀ = ∂ₓE(T, z, a, Z, A, ρ)
        ∂ₓE₁ = ∂ₓE(T - ∂ₓE₀ * x, z, a, Z, A, ρ)
        ΔT = x * (∂ₓE₀ + ∂ₓE₁) / 2
        Eloss += ΔT
    else
        # stop
        x -= r
        println("Stop at x = $(depth - x |> u"cm")")
        Eloss += T
    end
    return Eloss |> u"MeV"
end

function energyloss(T::uEnergy, z, a, Z, A, ρ, depth::Number; Δx::uLength=Δx₀)
    energyloss(T, z, a, Z, A, ρ, depth * 1u"cm", Δx=Δx)
end

function energyloss(T::Number, z, a, Z, A, ρ, depth::uLength; Δx::uLength=Δx₀)
    energyloss(T * 1u"MeV", z, a, Z, A, ρ, depth, Δx=Δx)
end

function energyloss(T::Number, z, a, Z, A, ρ, depth::Number; Δx::uLength=Δx₀)
    energyloss(T * 1u"MeV", z, a, Z, A, ρ, depth * 1u"cm", Δx=Δx)
end

end