Appendix D — Python Toolkit Reference

Throughout this textbook, you have built a Nuclear Data Analysis Toolkit — a collection of Python modules that implement the physics from each chapter as working, testable code. This appendix serves as a complete reference: what each module does, how the modules fit together, common API patterns, and example workflows that combine multiple modules to answer real nuclear physics questions.

The toolkit is not a substitute for professional nuclear data codes such as TALYS, EMPIRE, GEANT4, or MESA. It is a learning toolkit — every line of code corresponds to physics you derived in the text, and every plot you generate reinforces the physical intuition behind the calculation. By the time you complete the capstone project (Chapter 34), the toolkit will be a portfolio-quality demonstration of your ability to translate nuclear physics into quantitative computation.

D.1 Installation and Setup

D.1.1 Requirements

The toolkit requires Python 3.9 or later and the following packages:

numpy>=1.24.0
scipy>=1.10.0
matplotlib>=3.7.0
pandas>=2.0.0

Optional packages (used in specific modules):

periodictable>=1.7.0    # Nuclear data lookups (Ch. 2)
jupyter>=1.0.0          # Interactive notebooks
plotly>=5.15.0          # Interactive 3D plots (Ch. 23, 25)

D.1.2 Installation

# Clone or download the textbook repository
git clone https://github.com/datafield-dev/nuclear-physics-textbook.git
cd nuclear-physics-textbook

# Create a virtual environment (recommended)
python -m venv .venv
source .venv/bin/activate  # Linux/macOS
.venv\Scripts\activate     # Windows

# Install dependencies
pip install -r requirements.txt

D.1.3 Directory Structure

Each chapter's code lives in its own code/ subdirectory:

nuclear-physics/
├── part-01-foundations/
│   ├── chapter-01-discovery-of-the-nucleus/
│   │   └── code/
│   │       ├── binding_energy_curve.py
│   │       ├── nuclide_chart.py
│   │       └── project-checkpoint.md
│   ├── chapter-02-nuclear-properties/
│   │   └── code/
│   │       ├── nuclear_data_access.py
│   │       └── project-checkpoint.md
│   ...
├── requirements.txt
└── appendices/
    └── d-python-toolkit.md  (this file)

To import modules from other chapters, either add the relevant code/ directory to your Python path or copy the needed file into your working directory. Each module is designed to be self-contained — it imports only from the standard library and the packages listed in requirements.txt.

D.2 Module Reference

The following table lists every module in the toolkit, organized by chapter. Each entry gives the module name, the chapter where it is introduced, a brief description, and the key functions or classes it provides.

Part I: Foundations

Part II: Nuclear Structure

Part III: Radioactive Decay

Part IV: Nuclear Reactions

Part V: Nuclear Astrophysics

Part VI: Applications

Part VIII: Capstone

D.3 Common API Patterns

The toolkit follows consistent design patterns across all modules.

D.3.1 Data Loading

All modules that need nuclear data use nuclear_data_access.py as the data layer:

from nuclear_data_access import get_binding_energy, get_mass, get_half_life

# Get binding energy in MeV
B = get_binding_energy(Z=82, A=208)  # 208Pb → 1636.430 MeV

# Get atomic mass in u
m = get_mass(Z=92, A=235)  # 235U → 235.043930 u

# Get half-life as (value, unit) tuple
t_half = get_half_life(Z=27, A=60)  # 60Co → (5.2714, 'yr')

D.3.2 Plotting Convention

All plotting functions follow a consistent interface:

def plot_something(data, ax=None, save_path=None, **kwargs):
    """
    Parameters
    ----------
    data : array-like or specific input
        The physics input.
    ax : matplotlib.axes.Axes, optional
        If provided, plot on this axis. If None, create a new figure.
    save_path : str, optional
        If provided, save figure to this path (PDF recommended).
    **kwargs : dict
        Passed to the underlying matplotlib call.

    Returns
    -------
    fig, ax : matplotlib Figure and Axes
    """

This pattern allows both standalone use (plot_something(data) pops up a figure) and embedding in multi-panel figures:

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
plot_binding_energy_curve(ax=axes[0])
plot_semf_residuals(ax=axes[1])
plot_shell_energies(Z=50, ax=axes[2])
fig.tight_layout()
plt.savefig('nuclear_structure_overview.pdf')

D.3.3 Physical Constants

All modules import physical constants from a shared dictionary to ensure consistency:

# Every module uses these values (from Appendix C)
CONSTANTS = {
    'hbar_c': 197.3269804,       # MeV fm
    'e2': 1.4399764,             # MeV fm (= e^2 / 4 pi epsilon_0)
    'mp': 938.27208816,          # MeV/c^2
    'mn': 939.56542052,          # MeV/c^2
    'me': 0.51099895,            # MeV/c^2
    'u_to_MeV': 931.49410242,    # MeV/c^2 per u
    'alpha': 1/137.035999084,    # fine structure constant
    'r0': 1.2,                   # fm (nuclear radius parameter)
    'mu_N': 3.15245125844e-14,   # MeV/T (nuclear magneton)
}

D.3.4 Error Handling

Modules validate physical inputs and raise informative exceptions:

# This will raise ValueError with a physics-based message
from reaction_kinematics import compute_q_value
Q = compute_q_value(projectile=(1, 1), target=(8, 16),
                    ejectile=(0, 1), residual=(9, 16))
# ValueError: Reaction 16O(p,n)16F: Q = -16.19 MeV.
# Threshold energy = 17.20 MeV (lab). Specify E_lab >= threshold.

D.4 Example Workflows

D.4.1 Complete Analysis of an Alpha Decay

This workflow combines three modules to analyze the alpha decay of $^{212}$Po:

import numpy as np
import matplotlib.pyplot as plt

# Step 1: Get nuclear data
from nuclear_data_access import get_binding_energy, get_mass

B_parent = get_binding_energy(Z=84, A=212)    # 212Po
B_daughter = get_binding_energy(Z=82, A=208)  # 208Pb
B_alpha = get_binding_energy(Z=2, A=4)        # 4He

Q_alpha = B_daughter + B_alpha - B_parent
print(f"Q(alpha) = {Q_alpha:.3f} MeV")       # 8.954 MeV

# Step 2: Compute tunneling probability and half-life
from alpha_tunneling import gamow_half_life

t_half = gamow_half_life(Z_parent=84, A_parent=212, Q_alpha=Q_alpha)
print(f"Predicted t_1/2 = {t_half:.2e} s")    # Compare to 0.299 μs

# Step 3: Plot the Coulomb barrier and tunneling wavefunction
from alpha_tunneling import plot_barrier

fig, ax = plot_barrier(Z_daughter=82, A_daughter=208,
                       Q_alpha=Q_alpha, show_wkb=True)
ax.set_title(r'$^{212}$Po $\alpha$ Decay: Coulomb Barrier')
plt.savefig('po212_alpha_barrier.pdf')

D.4.2 Stellar Burning: From pp Chain to Silicon Burning

import numpy as np
import matplotlib.pyplot as plt
from stellar_burning import energy_generation_rate

# Compute energy generation for hydrogen burning
T_range = np.logspace(6.5, 8.0, 200)  # 3 × 10^6 to 10^8 K
rho = 150.0  # g/cm^3 (solar core-like)
X_H, X_He = 0.70, 0.28  # mass fractions

eps_pp = [energy_generation_rate('pp', T, rho, X_H=X_H) for T in T_range]
eps_cno = [energy_generation_rate('cno', T, rho, X_H=X_H, X_CNO=0.02)
           for T in T_range]

fig, ax = plt.subplots(figsize=(8, 6))
ax.loglog(T_range, eps_pp, 'b-', label='pp chain')
ax.loglog(T_range, eps_cno, 'r-', label='CNO cycle')
ax.axvline(1.57e7, color='gold', ls='--', label=r'$T_\odot$ core')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel(r'$\varepsilon$ (erg g$^{-1}$ s$^{-1}$)')
ax.set_title('Hydrogen Burning: pp vs CNO')
ax.legend()
plt.savefig('pp_vs_cno.pdf')

D.4.3 Neutron Star Mass-Radius Curve

from tov_solver import solve_tov, plot_mass_radius
import matplotlib.pyplot as plt

# Solve TOV for multiple equations of state
eos_names = ['SLy4', 'APR', 'Polytrope_gamma2']
fig, ax = plt.subplots(figsize=(8, 6))

for eos in eos_names:
    M_array, R_array = solve_tov(eos=eos, central_density_range=(1e14, 3e15),
                                  n_points=200)
    ax.plot(R_array, M_array, label=eos)

# Add observational constraints
ax.axhspan(1.97, 2.18, alpha=0.2, color='gray',
           label=r'PSR J0740+6620 ($2.08 \pm 0.07\; M_\odot$)')

ax.set_xlabel('Radius (km)')
ax.set_ylabel(r'Mass ($M_\odot$)')
ax.set_xlim(8, 16)
ax.set_ylim(0, 2.5)
ax.set_title('Neutron Star Mass-Radius Relation')
ax.legend(fontsize=9)
plt.savefig('neutron_star_mr.pdf')

D.4.4 Capstone: Full System Analysis

The capstone module provides a one-command analysis of any nucleus:

from capstone_analysis import NuclearSystem

# Analyze 56Fe — the peak of the B/A curve
fe56 = NuclearSystem(Z=26, A=56)

# Generate complete report
fe56.report(output_dir='fe56_analysis/')
# Creates:
#   fe56_analysis/binding_energy.pdf     — SEMF comparison
#   fe56_analysis/shell_structure.pdf    — Level diagram, J^pi prediction
#   fe56_analysis/decay_properties.pdf   — Stability analysis
#   fe56_analysis/reaction_summary.pdf   — Key reaction Q-values
#   fe56_analysis/astro_context.pdf      — Nucleosynthesis site
#   fe56_analysis/report.md              — Written summary

# Access individual properties
print(f"B/A = {fe56.binding_energy_per_nucleon:.4f} MeV")
print(f"SEMF prediction: {fe56.semf_binding_energy:.1f} MeV")
print(f"SEMF residual: {fe56.semf_residual:.1f} MeV")
print(f"Shell model J^pi: {fe56.predicted_spin_parity}")
print(f"Astrophysical origin: {fe56.nucleosynthesis_site}")

D.5 Module Dependency Graph

The modules build on each other in a logical progression that mirrors the textbook structure. The core dependency chain is:

nuclear_data_access.py  (Ch. 2 — foundation data layer)
    ├── binding_energy_curve.py  (Ch. 1)
    ├── nuclide_chart.py  (Ch. 1)
    ├── semf_fit.py  (Ch. 4)
    ├── shell_model.py  (Ch. 6)
    ├── alpha_tunneling.py  (Ch. 13)
    │       └── uses angular_momentum.py  (Ch. 5) for barrier calculations
    ├── beta_spectrum.py  (Ch. 14)
    ├── decay_chains.py  (Ch. 12)
    ├── reaction_kinematics.py  (Ch. 17)
    ├── breit_wigner.py  (Ch. 18)
    ├── fission_yields.py  (Ch. 20)
    ├── fusion_rates.py  (Ch. 21)
    │       └── uses gamow_peak.py  (Ch. 21) for rate integrals
    ├── stellar_burning.py  (Ch. 22)
    │       └── uses fusion_rates.py  (Ch. 21) for reaction rates
    ├── r_process_path.py  (Ch. 23)
    │       └── uses nuclide_chart.py  (Ch. 1) for visualization
    ├── bbn_network.py  (Ch. 24)
    │       └── uses fusion_rates.py  (Ch. 21) for rate inputs
    ├── tov_solver.py  (Ch. 25) — self-contained ODE solver
    ├── reactor_physics.py  (Ch. 26)
    └── capstone_analysis.py  (Ch. 34 — imports ALL above)

Standalone modules (no nuclear data dependency): - angular_momentum.py / clebsch_gordan.py (Ch. 5) — pure angular momentum algebra - weisskopf.py / transition_rates.py (Ch. 9) — needs only $A$ and multipole order - stopping_power.py / bragg_peak.py (Ch. 16) — needs only projectile and target material - dose_calculation.py (Ch. 27) — needs only source and geometry parameters - accelerator_physics.py (Ch. 30) — needs only beam and magnet parameters

D.6 Troubleshooting

D.6.1 Common Issues

"Module not found" when importing across chapters

Each module lives in its chapter's code/ directory. To import from another chapter:

import sys
sys.path.insert(0, '/path/to/nuclear-physics/part-01-foundations/'
                   'chapter-02-nuclear-properties/code')
from nuclear_data_access import get_binding_energy

Or set the PYTHONPATH environment variable to include all code/ directories.

AME2020 data file not found

The nuclear_data_access.py module expects a data/ directory containing ame2020.csv (or equivalent). If the data file is missing, the module falls back to a built-in subset of ~300 commonly used nuclei. To use the full AME2020 dataset, download it from the AME website and place it in the data/ directory as described in Chapter 2.

Matplotlib plots not appearing

If running in a non-interactive environment (script, SSH session), add at the top of your script:

import matplotlib
matplotlib.use('Agg')  # Non-interactive backend
import matplotlib.pyplot as plt

Plots will be saved to files but will not display interactively. In Jupyter notebooks, use %matplotlib inline or %matplotlib widget.

scipy.special functions unavailable

Some angular momentum functions (wigner_3j, wigner_6j, clebsch_gordan) were added to scipy.special in SciPy 1.11.0 (June 2023). If your SciPy version is older:

pip install --upgrade scipy>=1.11.0

If upgrading is not possible, the angular_momentum.py module includes pure-Python fallback implementations that are slower but produce identical results.

Numerical instabilities in the TOV solver

The TOV equation becomes stiff near the neutron star surface where the pressure drops rapidly. If solve_tov fails or produces unphysical results:

# Use tighter tolerances and a stiff solver
M, R = solve_tov(eos='SLy4',
                 central_density=5e14,
                 method='Radau',        # stiff-capable
                 rtol=1e-10,
                 atol=1e-12,
                 max_step=10.0)         # meters

Slow performance in the BBN network solver

The bbn_network.py module solves a stiff ODE system with widely separated timescales. If the solver is slow or fails to converge, ensure you are using a stiff integrator:

from bbn_network import solve_bbn
abundances = solve_bbn(eta=6.1e-10, method='BDF', rtol=1e-8)

The default BDF (backward differentiation formula) method handles the stiffness well. Avoid RK45 for this problem.

D.6.2 Validation

Every module includes a __main__ block that runs self-validation when executed directly:

python binding_energy_curve.py    # Plots B/A curve, compares to AME2020
python semf_fit.py                # Fits SEMF, prints parameters
python alpha_tunneling.py         # Computes 10 known alpha half-lives
python tov_solver.py              # Generates M-R curve for SLy4 EOS

Expected output and reference values are documented in each module's docstring and in the corresponding chapter's project-checkpoint.md file.

D.6.3 Unit Conventions

The toolkit uses the following unit conventions throughout:

All functions document their expected input and output units in their docstrings.

D.7 Extending the Toolkit

The toolkit is designed to be extended. Each module follows a consistent structure:

"""
module_name.py — Brief description
Chapter N: Chapter Title

Physics: One-paragraph summary of the physics implemented.
References: Krane Ch. X, Wong Ch. Y, or original paper.
"""

import numpy as np
import matplotlib.pyplot as plt
# ... other imports

# Physical constants (from Appendix C)
HBAR_C = 197.3269804  # MeV fm
# ... etc.

def main_physics_function(input_params):
    """Compute the main physical quantity.

    Parameters
    ----------
    input_params : type
        Description with units.

    Returns
    -------
    result : type
        Description with units.

    Examples
    --------
    >>> main_physics_function(82, 208)
    1636.430
    """
    # Implementation
    pass

def plot_results(data, ax=None, save_path=None):
    """Standard plotting function."""
    pass

if __name__ == '__main__':
    # Self-validation: compare to known results
    print("Running validation...")
    result = main_physics_function(test_input)
    expected = known_value
    assert abs(result - expected) < tolerance, f"FAIL: got {result}, expected {expected}"
    print(f"PASS: {result:.3f} (expected {expected:.3f})")

    # Generate standard plot
    plot_results(result, save_path='validation_plot.pdf')
    print("Validation complete.")

To add a new module: 1. Create it in the appropriate chapter's code/ directory. 2. Follow the template above (docstring, constants, functions, __main__ validation). 3. Use nuclear_data_access.py for any nuclear data lookups. 4. Use the standard plotting convention (accept ax parameter). 5. Document units explicitly in all docstrings. 6. Update the project-checkpoint.md in the same directory.

The complete toolkit source code is available in each chapter's code/ directory. For the latest version, corrections, and additional examples, see the textbook repository. Contributions and bug reports are welcome via the issue tracker.