Nuclear Physics

Code Examples

Nuclear Properties

import numpy as np

# Physical constants
c = 2.998e8  # m/s
u = 1.66054e-27  # Atomic mass unit (kg)
e = 1.602e-19  # Elementary charge (C)
MeV_to_J = 1.602e-13  # MeV to Joules

# Nuclear masses and binding energy
def atomic_mass_excess(A, Z, mass_excess_MeV):
    """Convert mass excess to atomic mass in u."""
    return A + mass_excess_MeV / 931.494

def binding_energy(A, Z, mass_nucleus):
    """Calculate binding energy from nuclear mass."""
    mp = 1.007276  # Proton mass (u)
    mn = 1.008665  # Neutron mass (u)
    mass_defect = Z * mp + (A - Z) * mn - mass_nucleus
    return mass_defect * 931.494  # MeV

def semi_empirical_mass(A, Z):
    """
    Semi-empirical mass formula (Weizsäcker formula).
    B(A,Z) = a_v A - a_s A^(2/3) - a_c Z(Z-1)/A^(1/3) - a_a (A-2Z)²/A ± δ
    """
    a_v = 15.75
    a_s = 17.8
    a_c = 0.711
    a_a = 23.7
    a_p = 11.18
    
    # Pairing term
    if A % 2 == 0 and Z % 2 == 0:
        delta = a_p / np.sqrt(A)
    elif A % 2 == 1 and Z % 2 == 1:
        delta = -a_p / np.sqrt(A)
    else:
        delta = 0
    
    B = (a_v * A 
         - a_s * A**(2/3) 
         - a_c * Z * (Z - 1) / A**(1/3) 
         - a_a * (A - 2*Z)**2 / A 
         + delta)
    
    return B

def binding_energy_per_nucleon(B, A):
    """Calculate binding energy per nucleon."""
    return B / A

# Example calculations
print("Nuclear binding energy:")
print(f"  Semi-empirical for ⁴⁰Ca (Z=20, A=40):")
B_Ca40 = semi_empirical_mass(40, 20)
print(f"    Total B = {B_Ca40:.1f} MeV")
print(f"    B/A = {B_Ca40/40:.3f} MeV/nucleon")

print(f"\n  Semi-empirical for ⁵⁶Fe (Z=26, A=56):")
B_Fe56 = semi_empirical_mass(56, 26)
print(f"    Total B = {B_Fe56:.1f} MeV")
print(f"    B/A = {B_Fe56/56:.3f} MeV/nucleon")

print(f"\n  Semi-empirical for ²³⁸U (Z=92, A=238):")
B_U238 = semi_empirical_mass(238, 92)
print(f"    Total B = {B_U238:.1f} MeV")
print(f"    B/A = {B_U238/238:.3f} MeV/nucleon")

# Magic numbers (nuclear shell closure)
magic_numbers = [2, 8, 20, 28, 50, 82, 126]
print(f"\nMagic numbers (shell closures): {magic_numbers}")
print("  Nuclei with magic numbers have:")
print("    - Extra binding energy")
print("    - Larger neutron separation energies")
print("    - Spherical ground states")
print("    - Higher first excited states")

Radioactive Decay

import numpy as np

def decay_constant(half_life):
    """λ = ln(2) / t_1/2"""
    return np.log(2) / half_life

def activity(N, half_life):
    """A = λN = (ln2/t_1/2) × N"""
    return decay_constant(half_life) * N

def radioactive_decay(N0, half_life, t):
    """N(t) = N₀ × 2^(-t/t_1/2) = N₀ × exp(-λt)"""
    return N0 * np.exp(-decay_constant(half_life) * t)

def daughter_growth(N0, half_life_parent, half_life_daughter, t):
    """Bateman equations for decay chain."""
    lambda_p = decay_constant(half_life_parent)
    lambda_d = decay_constant(half_life_daughter)
    
    N_parent = N0 * np.exp(-lambda_p * t)
    N_daughter = (lambda_p / (lambda_d - lambda_p)) * N0 * (
        np.exp(-lambda_p * t) - np.exp(-lambda_d * t)
    )
    
    return N_parent, N_daughter

def secular_equilibrium(t, lambda_p, lambda_d):
    """When t >> t_parent, N_daughter/N_parent ≈ λ_p/λ_d."""
    return lambda_p / lambda_d

# Decay constants and half-lives
half_lives = {
    'U-238': 4.468e9,      # years
    'U-235': 7.038e8,      # years
    'Ra-226': 1600,        # years
    'C-14': 5730,          # years
    'Co-60': 5.27,         # years
    'I-131': 8.02,         # days
    'Cs-137': 30.17,       # years
    'Rn-222': 3.82,         # days
}

print("Radioactive decay:")
for isotope, t_half in list(half_lives.items())[:5]:
    lambda_decay = decay_constant(t_half)
    print(f"  {isotope}: t_1/2 = {t_half:.2e}")
    print(f"    λ = {lambda_decay:.2e} s⁻¹")
    print(f"    1/λ = {1/lambda_decay:.2e} s = {1/lambda_decay/3600:.2e} hours")

# Carbon-14 dating
def carbon14_age(C14_ratio, t_half=5730):
    """
    Calculate age from C-14 ratio.
    Assumes initial ratio was 1.25e-12 (living organism)
    """
    lambda_c14 = decay_constant(t_half * 365.25 * 24 * 3600)
    return np.log(C14_ratio) / (-lambda_c14) / (365.25 * 24 * 3600)

# Samples with different C-14 ratios
initial_ratio = 1.25e-12
for ratio in [1.25e-12, 0.625e-12, 0.156e-12, 0.039e-12]:
    age = carbon14_age(ratio / initial_ratio)
    print(f"\n  C-14 ratio = {ratio:.3e}: age = {age:.0f} years")
    print(f"    ~ {age/1000:.1f} thousand years")

# Decay chain (U-238 series)
U238_hl = 4.468e9 * 365.25 * 24 * 3600  # seconds
Pb206_hl = np.inf  # Stable

t_equilibrium = 1e9 * 365.25 * 24 * 3600  # 1 billion years
N_U238 = 1e6  # Atoms
ratio_Pb_U = np.exp(-decay_constant(U238_hl) * t_equilibrium)
print(f"\nU-238 → Pb-206 decay chain:")
print(f"  After 1 billion years: Pb-206/U-238 ≈ {ratio_Pb_U:.2e}")
print(f"  Secular equilibrium reached when t >> 4.5 billion years")

Nuclear Reactions and Q-Values

import numpy as np

def Q_value(m_initial, m_final, c2=931.494):
    """
    Calculate Q-value of nuclear reaction.
    Q = (m_initial - m_final) × c²
    Positive Q: exothermic (releases energy)
    Negative Q: endothermic (absorbs energy)
    """
    return (m_initial - m_final) * c2

def threshold_energy(Q, A_target, A_projectile, c2=931.494):
    """
    Minimum kinetic energy for endothermic reaction.
    E_thresh = |Q| × (A_target + A_projectile) / A_target
    """
    return abs(Q) * (A_target + A_projectile) / A_target

def reaction_rate(n1, n2, sigma, v):
    """R = n₁n₂σv for two-body reactions."""
    return n1 * n2 * sigma * v

def cross_section_beer_lambert(I0, I, x):
    """σ = (1/nx) × ln(I₀/I) from Beer-Lambert law."""
    return np.log(I0 / I) / x

# Masses in atomic mass units (u)
masses_u = {
    'n': 1.008665,
    'p': 1.007825,
    'D': 2.014102,    # Deuterium
    'T': 3.016049,   # Tritium
    'He3': 3.016029,
    'He4': 4.002602,
    'U235': 235.043930,
    'U236': 236.045568,
    'Ba141': 140.914411,
    'Kr92': 91.926156,
    'fission_fragments': 140.9 + 91.9,  # Approximate
}

# Fusion reactions
print("Nuclear reaction Q-values:")
print("\n1. D + T → He-4 + n:")
m_init = masses_u['D'] + masses_u['T']
m_final = masses_u['He4'] + masses_u['n']
Q_DT = Q_value(m_init - m_final)
print(f"   Q = {Q_DT:.1f} MeV (releases energy!)")

print("\n2. D + D → He-3 + n:")
m_init = 2 * masses_u['D']
m_final = masses_u['He3'] + masses_u['n']
Q_DDn = Q_value(m_init - m_final)
print(f"   Q = {Q_DDn:.1f} MeV")

print("\n3. D + D → T + p:")
m_init = 2 * masses_u['D']
m_final = masses_u['T'] + masses_u['p']
Q_DDp = Q_value(m_init - m_final)
print(f"   Q = {Q_DDp:.1f} MeV")

print("\n4. U-235 + n → U-236*:")
m_init = masses_u['U235'] + masses_u['n']
m_final = masses_u['U236']
Q_fission = Q_value(m_init - m_final)
print(f"   Q = {Q_fission:.1f} MeV")

print("\n5. U-236 fission products (typical):")
Q_fission_products = Q_value(masses_u['U236'] - (masses_u['Ba141'] + masses_u['Kr92'] + 3*masses_u['n']))
print(f"   Ba-141 + Kr-92 + 3n: Q = {Q_fission_products:.1f} MeV")
print(f"   Total fission Q ≈ 200 MeV")

# Threshold energy
A_target, A_proj = 12, 1  # C-12 + n
E_thresh = threshold_energy(-5 MeV), A_target, A_proj)
print(f"\nThreshold energy (¹²C + n → ¹³C):")
print(f"   Q = -1.95 MeV (endothermic)")
print(f"   E_thresh = {E_thresh:.2f} MeV")

Nuclear Reactors

import numpy as np

# Neutron transport
def neutron_diffusion(D, Σa, phi):
    """Diffusion equation: D∇²φ - Σaφ = -S"""
    # Simplified 1D slab solution
    L = np.sqrt(D / Σa)  # Diffusion length
    return phi

def four_factor_formula(k_inf, epsilon, p, f):
    """
    Four-factor formula for infinite multiplication.
    k∞ = η × f × ε × p
    η = neutrons per absorption in fuel
    f = thermal utilization factor
    ε = fast fission factor
    p = resonance escape probability
    """
    return eta * f * epsilon * p

def effective_multiplication(k_inf, k_eff):
    """k_eff = k∞ × P_non_leakage"""
    return k_inf * (1 - k_non_leakage)

def reactor_period(T, rho):
    """Reactor period from reactivity."""
    # inhour equation (simplified)
    return T / (1 - 1/k_eff)

# Criticality calculations
k_inf = 1.32  # For U-235 thermal reactor
L = 0.10  # Non-leakage probability for large reactor
k_eff = k_inf * L

print("Nuclear reactor criticality:")
print(f"  k∞ = {k_inf:.2f}")
print(f"  k_eff = {k_eff:.2f}")
print(f"  Reactivity ρ = (k-1)/k = {(k_eff-1)/k_eff*100:.2f}%")

if k_eff > 1:
    print(f"  Supercritical - power increasing!")
elif k_eff < 1:
    print(f"  Subcritical - power decreasing!")