1. The Problem: The Invisible Power and Peril of Magnetic Fields
From the life-saving images of Magnetic Resonance Imaging (MRI) to the colossal energies of the Large Hadron Collider, powerful electromagnets are cornerstone technologies of modern science and medicine. They operate on a simple principle—running an electrical current through a coil of wire to generate a magnetic field. Yet, this field, while immensely useful, is invisible, non-intuitive, and can extend far beyond the physical boundaries of the device. This presents a fundamental challenge for the engineers and physicists who design them: how can we confidently predict a magnet’s performance and ensure its safety before a single wire is wound?
A simple textbook formula for an “ideal” solenoid can estimate the field deep within its core, but it falls apart in the real world. Real magnets are of finite length and have windings with physical thickness, creating complex “fringing fields” at their ends. It is in these details that performance is defined and safety is determined.
2. The Motivation: Why Simulation is Essential
Relying on rough estimates or building physical prototypes for high-power magnets is both prohibitively expensive and dangerously impractical. The motivation for accurate simulation is therefore driven by three critical needs:
- Performance Prediction: An engineer designing a 7 Tesla magnet for a research application needs to know if their design will actually achieve 7 T, and over what volume the field will be uniform. An underperforming magnet can render an entire experiment useless, while an overpowered one might be unnecessarily costly. Simulation allows designers to iterate through dozens of configurations—adjusting the coil’s length, thickness, and current—to meet precise specifications.
- Safety and Regulatory Compliance: This is arguably the most important driver for simulation. Stray magnetic fields can pose serious risks:
- High-Field Occupational Hazards: Extremely strong fields (e.g., above the 2 Tesla head exposure limit cited in French law) can induce vertigo or nausea in personnel. More critically, they can exert powerful forces on ferromagnetic tools, creating a lethal projectile risk (the “missile effect”), and can instantly disable or destroy critical medical implants like pacemakers. Simulation is essential to map out this high-risk zone.
- Low-Field Public Hazards: Even weak stray fields can be dangerous. The 5 Gauss (0.5 mT) limit is established to protect the public and those with sensitive medical devices from interference. A simulation must be able to predict how far this “5 Gauss line” extends into surrounding rooms, hallways, or even adjacent floors, defining the mandatory public exclusion zone.
- Cost Optimization: The materials for high-power magnets, particularly superconducting wire, are extraordinarily expensive. The operational costs, especially for cooling cryogenics, are also substantial. Simulation provides the tool to optimize a design, achieving the target magnetic field with the minimum amount of material and the most efficient geometry, saving potentially millions of dollars in construction and operational costs.
To address these challenges, we turn to numerical simulation. While complex commercial Finite Element Method (FEM) packages are the industry gold standard, a powerful and highly instructive simulation can be built from first principles. This article presents a Python-based tool that does just that. By implementing the fundamental Biot-Savart law, we will build a magnet model from the ground up, providing deep physical insight while delivering the accurate field maps needed for a comprehensive safety and performance analysis.
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import ellipk, ellipe
# --- 1. DEFINE PARAMETERS ---
# Coil Geometry
I_total = 3950000.0 # Total Ampere-turns
R_in = 0.15 # Inner radius (m)
R_out = 0.20 # Outer radius (m)
L = 0.6 # Length of the coil (m)
# Discretization (Resolution of the simulation)
N_radial = 40
N_axial = 50
# Calculation Domain (Area of interest)
z_pts = np.linspace(0, 6.0, 200) # z from 0 to 6.0m
r_pts = np.linspace(0, 6.0, 200) # r from 0 to 6.0m
# --- 2. SETUP AND CALCULATION ---
print("Starting simulation... This may take several minutes.")
N_filaments = N_radial * N_axial
i_filament = I_total / N_filaments
Z, R_grid = np.meshgrid(z_pts, r_pts)
r_filaments = np.linspace(R_in, R_out, N_radial)
z_filaments = np.linspace(-L/2, L/2, N_axial)
Bz_total = np.zeros_like(Z)
Br_total = np.zeros_like(Z)
# Main calculation loop (numerical integration)
for r_c in r_filaments:
for z_c in z_filaments:
z_rel = Z - z_c
m = (4 * r_c * R_grid) / ((r_c + R_grid)**2 + z_rel**2)
m = np.clip(m, 0, 1)
K = ellipk(m)
E = ellipe(m)
mu_0 = 4 * np.pi * 1e-7
prefactor = (mu_0 * i_filament) / (2 * np.pi)
denom_br = (r_c - R_grid)**2 + z_rel**2
denom_br[denom_br == 0] = 1e-12
sqrt_term = np.sqrt((r_c + R_grid)**2 + z_rel**2)
Bz_coil = (prefactor / sqrt_term) * (K + (r_c**2 - R_grid**2 - z_rel**2) / denom_br * E)
Br_coil = np.zeros_like(Z)
non_axis_mask = R_grid > 1e-9
term1 = (r_c**2 + R_grid[non_axis_mask]**2 + z_rel[non_axis_mask]**2)
term2 = denom_br[non_axis_mask]
Br_coil[non_axis_mask] = (prefactor * z_rel[non_axis_mask] / (R_grid[non_axis_mask] * sqrt_term[non_axis_mask])) * \
(-K[non_axis_mask] + (term1 / term2) * E[non_axis_mask])
Bz_total += Bz_coil
Br_total += Br_coil
print("Calculation complete.")
# --- 3. DATA ANALYSIS AND INTERPOLATION ---
B_mag = np.sqrt(Bz_total**2 + Br_total**2)
critical_5_gauss_T = 0.0005
critical_2_T = 2.0
# Find 5 Gauss positions
B_axial = B_mag[0, :]
z_5_gauss_m = np.interp(critical_5_gauss_T, B_axial[::-1], z_pts[::-1])
z0_index = np.argmin(np.abs(z_pts))
B_radial = B_mag[:, z0_index]
r_5_gauss_m = np.interp(critical_5_gauss_T, B_radial[::-1], r_pts[::-1])
# Find 2 Tesla positions
z_peak_idx = np.argmax(B_axial)
z_2T_m = np.interp(critical_2_T, B_axial[z_peak_idx:][::-1], z_pts[z_peak_idx:][::-1])
r_out_idx = np.argmin(np.abs(r_pts - R_out))
r_2T_m = np.interp(critical_2_T, B_radial[r_out_idx:][::-1], r_pts[r_out_idx:][::-1])
print(f"Definitive 5 Gauss limit found at z = +/- {z_5_gauss_m:.2f} m and r = +/- {r_5_gauss_m:.2f} m")
print(f"Definitive 2 Tesla limit found at z = +/- {z_2T_m:.2f} m and r = +/- {r_2T_m:.2f} m")
# --- 4. PLOTTING AND VISUALIZATION ---
# Reconstruct full plots for symmetric visualization
full_z_pts = np.concatenate((-z_pts[::-1], z_pts[1:]))
full_B_axial = np.concatenate((B_axial[::-1], B_axial[1:]))
full_r_pts = np.concatenate((-r_pts[::-1], r_pts[1:]))
full_B_radial = np.concatenate((B_radial[::-1], B_radial[1:]))
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 14))
# Plot 1: Axial Field
ax1.plot(full_z_pts, full_B_axial, label='Field on Axis')
ax1.axhline(y=critical_2_T, color='red', linestyle='-.', lw=2, label='2 T Head Limit')
ax1.axhline(y=critical_5_gauss_T, color='darkorange', linestyle='--', lw=2, label='5 Gauss Entrée Limit')
ax1.set_title('Fully Annotated Magnetic Field Along the Central Axis')
ax1.set_xlabel('Distance along z-axis (m)')
ax1.set_ylabel('Magnetic Field Strength |B| (Tesla)')
ax1.grid(True, linestyle='--', alpha=0.7)
ax1.legend(loc='upper right')
ax1.annotate(f'5 Gauss limit:\nz = +/- {z_5_gauss_m:.2f} m',
xy=(z_5_gauss_m, critical_5_gauss_T), xytext=(z_5_gauss_m - 2, 2.5),
arrowprops=dict(facecolor='black', shrink=0.05),
bbox=dict(boxstyle="round,pad=0.3", fc="wheat", ec="black", lw=1, alpha=0.9))
ax1.annotate(f'2 Tesla limit:\nz = +/- {z_2T_m:.2f} m',
xy=(z_2T_m, critical_2_T), xytext=(z_2T_m + 1, 3.5),
arrowprops=dict(facecolor='red', shrink=0.05, width=1.5, headwidth=8),
bbox=dict(boxstyle="round,pad=0.3", fc="mistyrose", ec="red", lw=1, alpha=0.9))
# Plot 2: Radial Field
ax2.plot(full_r_pts, full_B_radial, label='Field in Mid-Plane')
ax2.axhline(y=critical_2_T, color='red', linestyle='-.', lw=2, label='2 T Head Limit')
ax2.axhline(y=critical_5_gauss_T, color='darkorange', linestyle='--', lw=2, label='5 Gauss Entrée Limit')
ax2.set_title('Fully Annotated Magnetic Field in the Mid-Plane (z=0)')
ax2.set_xlabel('Distance from center, r (m)')
ax2.set_ylabel('Magnetic Field Strength |B| (Tesla)')
ax2.grid(True, linestyle='--', alpha=0.7)
ax2.legend(loc='upper right')
ax2.annotate(f'5 Gauss limit:\nr = +/- {r_5_gauss_m:.2f} m',
xy=(r_5_gauss_m, critical_5_gauss_T), xytext=(r_5_gauss_m - 2, 2.5),
arrowprops=dict(facecolor='black', shrink=0.05),
bbox=dict(boxstyle="round,pad=0.3", fc="wheat", ec="black", lw=1, alpha=0.9))
ax2.annotate(f'2 Tesla limit:\nr = +/- {r_2T_m:.2f} m',
xy=(r_2T_m, critical_2_T), xytext=(r_2T_m + 1.5, 3.5),
arrowprops=dict(facecolor='red', shrink=0.05, width=1.5, headwidth=8),
bbox=dict(boxstyle="round,pad=0.3", fc="mistyrose", ec="red", lw=1, alpha=0.9))
plt.tight_layout()
plt.show()