Skip to content

Commit

Permalink
Add field fitting, optional B0 to set nI
Browse files Browse the repository at this point in the history
  • Loading branch information
@ChristopherMayes
ChristopherMayes committed Dec 23, 2024
1 parent 6cdaf19 commit c879237
Show file tree
Hide file tree
Showing 3 changed files with 248 additions and 68 deletions.
154 changes: 88 additions & 66 deletions docs/examples/fields/solenoid_modeling.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,14 @@
"id": "0fe9a7b5-2538-4fb1-93f5-c3dd9e4b1023",
"metadata": {},
"source": [
"# Solenoid Modeling"
"# Solenoid Modeling\n",
"\n",
"We can create solenoid fieldmesh data using an ideal model described in:\n",
"\n",
"Derby, N., & Olbert, S. (2010). Cylindrical magnets and ideal solenoids.\n",
"American Journal of Physics, 78(3), 229–235. https://doi.org/10.1119/1.3256157\n",
"\n",
"(preprint: https://arxiv.org/abs/0909.3880)\n"
]
},
{
Expand All @@ -15,7 +22,9 @@
"metadata": {},
"outputs": [],
"source": [
"from pmd_beamphysics.fields.solenoid import make_solenoid_fieldmesh\n",
"from pmd_beamphysics.fields.solenoid import make_solenoid_fieldmesh, fit_ideal_solenoid\n",
"\n",
"from pmd_beamphysics.fields.analysis import solenoid_analysis\n",
"\n",
"from pmd_beamphysics.fields.analysis import check_static_div_equation\n",
"from pmd_beamphysics.units import mu_0\n",
Expand All @@ -37,11 +46,11 @@
" radius=0.05,\n",
" L=0.2,\n",
" rmax=0.1,\n",
" zmin=-0.4,\n",
" zmax=0.4,\n",
" nr=101,\n",
" nz=200,\n",
" nI=1,\n",
" zmin=-0.5,\n",
" zmax=0.5,\n",
" nr=51,\n",
" nz=100,\n",
" B0=1,\n",
")"
]
},
Expand All @@ -65,14 +74,41 @@
"FM.plot_onaxis()"
]
},
{
"cell_type": "markdown",
"id": "6b27bc8a-e07f-4f5c-980f-b07ee4cd64b7",
"metadata": {},
"source": [
"This calculates field integrals, and gives the equivalent hard-edge solenoid parameters:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d7c3bb9f-47ed-4b0d-86f0-c786e7ffec9f",
"metadata": {},
"outputs": [],
"source": [
"z0, Bz0 = FM.axis_values(\"z\", \"Bz\")\n",
"solenoid_analysis(z0, Bz0)"
]
},
{
"cell_type": "markdown",
"id": "ed57f328-93f2-449f-84c6-5300dc554d83",
"metadata": {},
"source": [
"Check the static divergence Maxwell equation"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8a317af4-d705-4481-bc3a-0ea9d767cbb2",
"metadata": {},
"outputs": [],
"source": [
"check_static_div_equation(FM, rtol=1e-2, plot=True)"
"check_static_div_equation(FM, rtol=1e-1, plot=True)"
]
},
{
Expand Down Expand Up @@ -126,6 +162,25 @@
"FM_hard.plot_onaxis()"
]
},
{
"cell_type": "markdown",
"id": "0d8fd0e8-1ed4-4c68-8a10-d56b0def6838",
"metadata": {},
"source": [
"Check hard-edge analysis again"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9abcf3c2-fb43-4d09-9b55-23499a9b3505",
"metadata": {},
"outputs": [],
"source": [
"z0, Bz0 = FM_hard.axis_values(\"z\", \"Bz\")\n",
"solenoid_analysis(z0, Bz0)"
]
},
{
"cell_type": "markdown",
"id": "bfb0168b-cf80-47d4-a77d-23391391d112",
Expand All @@ -148,26 +203,39 @@
{
"cell_type": "code",
"execution_count": null,
"id": "c8ce66f4-1b14-45f5-b0d7-6e402868cf05",
"id": "891bc4fb-497b-4848-9bad-b7321cb751e1",
"metadata": {},
"outputs": [],
"source": [
"z0, Bz0 = FM2.axis_values(\"z\", \"Bz\")\n",
"\n",
"fit = fit_ideal_solenoid(z0, Bz0)\n",
"fit"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8a74f0b7-86cd-4235-96cc-27a05f5c7f5c",
"metadata": {},
"outputs": [],
"source": [
"FM3 = make_solenoid_fieldmesh(\n",
" radius=0.0191,\n",
" L=0.02936 * 2,\n",
" radius=fit[\"radius\"],\n",
" L=fit[\"L\"],\n",
" rmax=0.1,\n",
" zmin=-0.1,\n",
" zmax=0.1,\n",
" nr=100,\n",
" nz=40,\n",
" nI=1,\n",
" B0=fit[\"B0\"],\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1313ec9a-2338-46ae-a6b0-03da50510c7a",
"id": "d40a160f-09e3-41a0-8eb3-66b0e3cca10d",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -182,72 +250,26 @@
"plt.legend()\n",
"ax.set_ylim(0, None)\n",
"ax.set_xlabel(r\"$z$ (m)\")\n",
"ax.set_ylabel(r\"B_z (T)\")"
"ax.set_ylabel(r\"$B_z$ (T)\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "78ca976c-b608-461d-95ce-d7b147bba952",
"cell_type": "markdown",
"id": "4f19b4b3-2d7b-4f47-a238-1f2c227c32ac",
"metadata": {},
"outputs": [],
"source": [
"FM2.plot()\n",
"FM3.plot()"
"Note that the fields are very different off-axis:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e61c96b2-2ad9-4920-bb80-6a4d54549ab3",
"id": "f4e59992-b296-4482-851e-25014c1ce9ec",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.optimize import curve_fit\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"# Define the Bz model from the given formula\n",
"def Bz_model(z, a, b):\n",
" \"\"\"\n",
" Normalized ideal on-axis model\n",
" \"\"\"\n",
" A = np.hypot(a, b) / (2 * b)\n",
" term1 = (z + b) / np.hypot(z + b, a)\n",
" term2 = (z - b) / np.hypot(z - b, a)\n",
" return A * (term1 - term2)\n",
"\n",
"\n",
"z_data = z0\n",
"Bz_data = Bz0\n",
"\n",
"# Normalize the data\n",
"Bz_data = Bz_data / np.max(Bz_data) # Normalize Bz to a maximum of 1\n",
"z_data = z_data - z_data[np.argmax(Bz_data)] # Shift z so maximum is at z=0\n",
"\n",
"# Fit the data to the model\n",
"popt, pcov = curve_fit(Bz_model, z_data, Bz_data, p0=[0.1, 0.1])\n",
"\n",
"# Extract fitted parameters\n",
"a_fit, b_fit = popt\n",
"a_fit = abs(a_fit)\n",
"b_fit = abs(b_fit)\n",
"print(f\"Fitted parameters: a = {a_fit}, b = {b_fit}\")\n",
"\n",
"# Plot the data and the fitted curve\n",
"z_fit = np.linspace(np.min(z_data), np.max(z_data), 500)\n",
"Bz_fit = Bz_model(z_fit, a_fit, b_fit)\n",
"\n",
"plt.figure(figsize=(8, 6))\n",
"plt.scatter(z_data, Bz_data, label=\"Normalized Data\", color=\"blue\", alpha=0.7)\n",
"plt.plot(z_fit, Bz_fit, label=\"Fitted Model\", color=\"red\", linewidth=2)\n",
"plt.xlabel(\"z (normalized)\")\n",
"plt.ylabel(\"Bz (normalized)\")\n",
"plt.title(\"Fitting Bz Model to Data\")\n",
"plt.legend()\n",
"plt.grid()\n",
"plt.show()"
"FM2.plot()\n",
"FM3.plot()"
]
}
],
Expand Down
58 changes: 58 additions & 0 deletions pmd_beamphysics/fields/analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1343,3 +1343,61 @@ def plot_maxwell_curl_equations_cartesian(FM, ix=None, iy=None, plot_diff=True):

fig.suptitle(rf"Fields evaluated at $x=${x[ix]:0.6f}, $y=${y[iy]:0.6f} meters.")
plt.tight_layout()


def solenoid_analysis(z0, Bz0):
r"""
This computes the hard edge equivalent magnetic field and length
from on-axis solenoid field data using two integrals:
$B_1 = \int B/B_0 dz$
$B_2 = \int (B/B_0)^2 dz$
These give:
$L_\text{hard} = B_1$^2/B_2$
$B_\text{hard} = B_0 * B_2 / B_1$
Parameters
----------
z0: array-like
On-axis z points in meters
Bz0: array-like
On-axis Bz points in T
Returns
-------
dict with:
'B0': float
maximum field in T
'int_BL': float
\int B(z) dz
'int_B2L': float
\int [B(z)]^2 dz
'L_hard': float
hard-edge length in meters
'B_hard':
hard-edge field in T
"""
B0 = Bz0.max()

Bz0 = Bz0 / B0
BL = np.trapezoid(Bz0, z0)
B2L = np.trapezoid(Bz0**2, z0)
d = {}
d["B0"] = B0
d["int_BL"] = BL * B0
d["int_B2L"] = B2L * B0**2
d["L_hard"] = BL**2 / B2L
d["B_hard"] = B2L / BL * B0

return d
Loading

0 comments on commit c879237

Please sign in to comment.