Solenoid Modeling¶
We can create solenoid fieldmesh data using an ideal model described in:
Derby, N., & Olbert, S. (2010). Cylindrical magnets and ideal solenoids. American Journal of Physics, 78(3), 229–235. https://doi.org/10.1119/1.3256157
(preprint: https://arxiv.org/abs/0909.3880)
In [1]:
Copied!
from pmd_beamphysics.fields.solenoid import make_solenoid_fieldmesh, fit_ideal_solenoid
from pmd_beamphysics.fields.analysis import solenoid_analysis
from pmd_beamphysics.fields.analysis import check_static_div_equation
from pmd_beamphysics.units import mu_0
from pmd_beamphysics import FieldMesh
import matplotlib.pyplot as plt
import numpy as np
from pmd_beamphysics.fields.solenoid import make_solenoid_fieldmesh, fit_ideal_solenoid
from pmd_beamphysics.fields.analysis import solenoid_analysis
from pmd_beamphysics.fields.analysis import check_static_div_equation
from pmd_beamphysics.units import mu_0
from pmd_beamphysics import FieldMesh
import matplotlib.pyplot as plt
import numpy as np
In [2]:
Copied!
FM = make_solenoid_fieldmesh(
radius=0.05,
L=0.2,
rmax=0.1,
zmin=-0.5,
zmax=0.5,
nr=51,
nz=100,
B0=1,
)
FM = make_solenoid_fieldmesh(
radius=0.05,
L=0.2,
rmax=0.1,
zmin=-0.5,
zmax=0.5,
nr=51,
nz=100,
B0=1,
)
In [3]:
Copied!
FM = make_solenoid_fieldmesh(L=0.1, B0=1, zmin=-0.2, zmax=0.2, rmax=0.01, radius=0.1)
FM.plot()
FM = make_solenoid_fieldmesh(L=0.1, B0=1, zmin=-0.2, zmax=0.2, rmax=0.01, radius=0.1)
FM.plot()
In [4]:
Copied!
FM.plot()
FM.plot()
In [5]:
Copied!
FM.plot(stream=True, mirror="r")
FM.plot(stream=True, mirror="r")
In [6]:
Copied!
FM.plot_onaxis()
FM.plot_onaxis()
This calculates field integrals, and gives the equivalent hard-edge solenoid parameters:
In [7]:
Copied!
z0, Bz0 = FM.axis_values("z", "Bz")
solenoid_analysis(z0, Bz0)
z0, Bz0 = FM.axis_values("z", "Bz")
solenoid_analysis(z0, Bz0)
Out[7]:
{'B0': np.float64(0.9974788851617891), 'int_BL': np.float64(0.19894343802696227), 'int_B2L': np.float64(0.13609559046233594), 'L_hard': np.float64(0.2908139154217565), 'B_hard': np.float64(0.6840918796421487)}
Check the static divergence Maxwell equation
In [8]:
Copied!
check_static_div_equation(FM, rtol=1e-1, plot=True)
check_static_div_equation(FM, rtol=1e-1, plot=True)
Out[8]:
np.True_
Hard edge¶
Making the radius very small approximates a hard-edge model.
Here we expext that $B_z = \mu_0 n I$
In [9]:
Copied!
FM_hard = make_solenoid_fieldmesh(
radius=1e-9,
L=0.2,
rmax=0.1,
zmin=-0.4,
zmax=0.4,
nr=101,
nz=200,
nI=1,
)
FM_hard = make_solenoid_fieldmesh(
radius=1e-9,
L=0.2,
rmax=0.1,
zmin=-0.4,
zmax=0.4,
nr=101,
nz=200,
nI=1,
)
In [10]:
Copied!
# Ckech the on-avis maximum field is mu_0
FM_hard.Bz[0, 0, :].max() == mu_0
# Ckech the on-avis maximum field is mu_0
FM_hard.Bz[0, 0, :].max() == mu_0
Out[10]:
np.True_
In [11]:
Copied!
FM_hard.plot_onaxis()
FM_hard.plot_onaxis()
Check hard-edge analysis again
In [12]:
Copied!
z0, Bz0 = FM_hard.axis_values("z", "Bz")
solenoid_analysis(z0, Bz0)
z0, Bz0 = FM_hard.axis_values("z", "Bz")
solenoid_analysis(z0, Bz0)
Out[12]:
{'B0': np.float64(1.25663706127e-06), 'int_BL': np.float64(2.525903640743712e-07), 'int_B2L': np.float64(3.1741441281553495e-13), 'L_hard': np.float64(0.20100502512562904), 'B_hard': np.float64(1.256637061269991e-06)}
Compare with a real solenoid¶
In [13]:
Copied!
FM2 = FieldMesh("../data/solenoid.h5")
FM2.plot_onaxis()
FM2 = FieldMesh("../data/solenoid.h5")
FM2.plot_onaxis()
In [14]:
Copied!
z0, Bz0 = FM2.axis_values("z", "Bz")
fit = fit_ideal_solenoid(z0, Bz0)
fit
z0, Bz0 = FM2.axis_values("z", "Bz")
fit = fit_ideal_solenoid(z0, Bz0)
fit
Out[14]:
{'B0': np.float64(1.0), 'radius': 0.019104535909029125, 'L': 0.058730560178559606}
In [15]:
Copied!
FM3 = make_solenoid_fieldmesh(
radius=fit["radius"],
L=fit["L"],
rmax=0.1,
zmin=-0.1,
zmax=0.1,
nr=100,
nz=40,
B0=fit["B0"],
)
FM3 = make_solenoid_fieldmesh(
radius=fit["radius"],
L=fit["L"],
rmax=0.1,
zmin=-0.1,
zmax=0.1,
nr=100,
nz=40,
B0=fit["B0"],
)
In [16]:
Copied!
fig, ax = plt.subplots()
z0, Bz0 = FM2.axis_values("z", "Bz")
Bz0 = np.real(Bz0 / Bz0.max())
ax.plot(z0, Bz0, label="Superfish")
z, Bz = FM3.axis_values("z", "Bz")
Bz = np.real(Bz / Bz.max())
ax.plot(z, Bz, label="Ideal")
plt.legend()
ax.set_ylim(0, None)
ax.set_xlabel(r"$z$ (m)")
ax.set_ylabel(r"$B_z$ (T)")
fig, ax = plt.subplots()
z0, Bz0 = FM2.axis_values("z", "Bz")
Bz0 = np.real(Bz0 / Bz0.max())
ax.plot(z0, Bz0, label="Superfish")
z, Bz = FM3.axis_values("z", "Bz")
Bz = np.real(Bz / Bz.max())
ax.plot(z, Bz, label="Ideal")
plt.legend()
ax.set_ylim(0, None)
ax.set_xlabel(r"$z$ (m)")
ax.set_ylabel(r"$B_z$ (T)")
Out[16]:
Text(0, 0.5, '$B_z$ (T)')
Note that the fields are very different off-axis:
In [17]:
Copied!
FM2.plot()
FM3.plot()
FM2.plot()
FM3.plot()
Stylized Plot¶
Here we make a stylized plot, similar to the figure Fig. 1 in Derby & Olbert (2010).
In [18]:
Copied!
FM = make_solenoid_fieldmesh(
radius=0.05,
L=0.2,
rmax=0.5,
zmin=-0.5,
zmax=0.5,
nr=200,
nz=200,
B0=1,
)
FM = make_solenoid_fieldmesh(
radius=0.05,
L=0.2,
rmax=0.5,
zmin=-0.5,
zmax=0.5,
nr=200,
nz=200,
B0=1,
)
In [19]:
Copied!
# Stylize with varying line thickness
linewidth = 0.5 + np.random.normal(scale=0.3, size=(2 * FM.shape[0] - 1, FM.shape[2]))
FM.plot(
stream=True,
mirror="r",
density=8,
linewidth=linewidth,
arrowsize=0,
figsize=(8, 8),
aspect="equal",
cmap="copper",
)
# Stylize with varying line thickness
linewidth = 0.5 + np.random.normal(scale=0.3, size=(2 * FM.shape[0] - 1, FM.shape[2]))
FM.plot(
stream=True,
mirror="r",
density=8,
linewidth=linewidth,
arrowsize=0,
figsize=(8, 8),
aspect="equal",
cmap="copper",
)
/home/runner/miniconda3/envs/beamphysics-dev/lib/python3.13/site-packages/matplotlib/patches.py:3436: RuntimeWarning: invalid value encountered in scalar divide cos_t, sin_t = head_length / head_dist, head_width / head_dist