Python implementation of Gamow model for Polonium 212

gamow.py

from math import log, pi, sqrt, exp

"""
Authors: CALHEGAS Raphaël, DIZIN Jordan, FROMENTIN Raphaël
"""

# Constants
electron_mass = 9.1e-31
h = 6.626e-34
bar_h = h / (2 * pi)
e = 1.6e-19
k_0 = 8.9875517923e9
r_0 = 1.07e-15
m_proton = 1.6726e-27
m_neutron = 1.9749e-27
m_alpha = 6.644656e-27
c = 3e18


def half_life(lambd):
    return log(2) / lambd


def disintegration_frequency(collision_freq, transmission_prob):
    return collision_freq * transmission_prob


def collision_frequency(v, distance):
    return v * distance


def transmission_probability(m, vm, R, e):
    S = pi * R * vm * sqrt(2 * m / e) - 2 * R * sqrt(m * vm)
    print(-S/bar_h)
    return exp(-S / bar_h)

def electro_potential(q, r):
    return k_0 * q / r


def coulomb_force(q1, q2, r):
    return k_0 * (q1 * q2) / (r ** 2)


def electro_energy(q, r):
    return 2 * e * electro_potential(q, r)


def kernel_radius(a):
    return r_0 * (a ** (1 / 3))


if __name__ == '__main__':
    # Original atom constants
    A = 212  # Number of nucleons
    Z = 84  # Number of protons

    # Alpha particle constants
    alpha_A = 4  # Number of nucleons
    alpha_protons = 2  # Number of protons

    # Radius
    R = kernel_radius(A)  # Original atom radius
    R_after_dzt = kernel_radius(A - alpha_A)  # Radius after loss of alpha particle
    R_alpha = kernel_radius(alpha_A)  # Radius of the alpha particle

    # Charge
    q1 = (Z - alpha_protons) * e  # Of the atom without the alpha particle
    q2 = alpha_protons * e  # Of the alpha particle
    # At the moment the atom and the particle separate from each other, they are touching each other.
    # As such, the distance between their centers is the sum of their radius.
    distance = R_after_dzt + R_alpha

    force = coulomb_force(q1, q2, distance)
    vm = electro_potential(q2, R)

    m_alpha = 2 * m_proton + 2 * m_neutron  # Mass of the alpha particle
    e_alpha_out = electro_energy(q2, distance)

    T = transmission_probability(m_alpha, vm, R, e_alpha_out)

    # (Strong) Binding energy = |E_nucleons - E_kernel| = |m_nucleons - m_kernel| <=> diff mass
    e_alpha_in = 2*m_proton + 2*m_neutron - m_alpha
    # E = 1/2 mc^2 <=> c = sqrt(2/m E)
    v_alpha = sqrt(2/m_alpha * e_alpha_in)

    f_col = collision_frequency(v_alpha, 2 * R)

    lambd = disintegration_frequency(f_col, T)

    print(f"{A=} {Z=} {vm=} {T=} {lambd=} {e_alpha_in=}")

    half = half_life(lambd)

    print(f"{half=}")  # 2.3012262527759235e-7 s