From c8792374dde74701f67582ca781cacf967781663 Mon Sep 17 00:00:00 2001 From: Christopher Mayes <31023527+ChristopherMayes@users.noreply.github.com> Date: Mon, 23 Dec 2024 11:21:42 -0800 Subject: [PATCH] Add field fitting, optional B0 to set nI --- docs/examples/fields/solenoid_modeling.ipynb | 154 +++++++++++-------- pmd_beamphysics/fields/analysis.py | 58 +++++++ pmd_beamphysics/fields/solenoid.py | 104 ++++++++++++- 3 files changed, 248 insertions(+), 68 deletions(-) diff --git a/docs/examples/fields/solenoid_modeling.ipynb b/docs/examples/fields/solenoid_modeling.ipynb index 0df8c06..4ea8648 100644 --- a/docs/examples/fields/solenoid_modeling.ipynb +++ b/docs/examples/fields/solenoid_modeling.ipynb @@ -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" ] }, { @@ -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", @@ -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", ")" ] }, @@ -65,6 +74,33 @@ "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, @@ -72,7 +108,7 @@ "metadata": {}, "outputs": [], "source": [ - "check_static_div_equation(FM, rtol=1e-2, plot=True)" + "check_static_div_equation(FM, rtol=1e-1, plot=True)" ] }, { @@ -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", @@ -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": [ @@ -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()" ] } ], diff --git a/pmd_beamphysics/fields/analysis.py b/pmd_beamphysics/fields/analysis.py index d3569b4..d4e5a8b 100644 --- a/pmd_beamphysics/fields/analysis.py +++ b/pmd_beamphysics/fields/analysis.py @@ -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 diff --git a/pmd_beamphysics/fields/solenoid.py b/pmd_beamphysics/fields/solenoid.py index 09b7eb7..3063a5b 100644 --- a/pmd_beamphysics/fields/solenoid.py +++ b/pmd_beamphysics/fields/solenoid.py @@ -3,6 +3,8 @@ from pmd_beamphysics.units import mu_0 from pmd_beamphysics import FieldMesh +from scipy.optimize import curve_fit + from scipy.integrate import quad @@ -266,7 +268,8 @@ def make_solenoid_fieldmesh( nz: int = 40, radius: float = None, L: float = None, - nI: float = 1.0, + nI: float = None, + B0: float = None, ): """ Generates a 2D cylindrically symmetric ideal solenoid FieldMesh. @@ -299,9 +302,16 @@ def make_solenoid_fieldmesh( The inner radius of the solenoid (in meters). Default is 0.1. L : float, optional The half length of the solenoid (in meters). Default is 0.2. + nI : float, optional The product of the number of turns per unit length (n) and current (I) in amperes. - This determines the current density of the solenoid. Default is 1.0. + This determines the current density of the solenoid. + default: None + + B0 : float, optional + Peak on-axis magnetic field in T + This can be used instead of nI to scale the field. + default: None Returns ------- @@ -314,6 +324,15 @@ def make_solenoid_fieldmesh( >>> print(field_mesh) """ + if nI is None and B0 is not None: + nI = B0 * np.hypot(radius, L / 2) / (mu_0 * L / 2) + elif nI is None and B0 is None: + B0 = 1 + elif nI is not None and B0 is None: + pass + else: + raise ValueError("Cannot set nI and B0. Choose one") + # Form coordinate mesh rs = np.linspace(rmin, rmax, nr) zs = np.linspace(zmin, zmax, nz) @@ -353,3 +372,84 @@ def make_solenoid_fieldmesh( data = dict(attrs=attrs, components=components) return FieldMesh(data=data) + + +def _fit_ideal_solenoid_model(z: np.ndarray, radius: float, L: float) -> np.ndarray: + """ + Computes the normalized magnetic field on the solenoid's axis. + + This function calculates the on-axis magnetic field of an ideal solenoid + normalized to a maximum value of 1. + + Parameters + ---------- + z : np.ndarray + Array of axial positions (z) where the magnetic field is evaluated. + radius : float + Radius of the solenoid (meters). + L : float + Total length of the solenoid (meters). + + Returns + ------- + np.ndarray + Normalized magnetic field values at the specified z positions. + """ + a = radius + b = L / 2 + + A = np.hypot(a, b) / (2 * b) + term1 = (z + b) / np.hypot(z + b, a) + term2 = (z - b) / np.hypot(z - b, a) + return A * (term1 - term2) + + +def fit_ideal_solenoid( + z0: np.ndarray, Bz0: np.ndarray, radius0: float = 0.1, L0: float = 0.1 +) -> dict: + """ + Fits the ideal solenoid model to experimental or simulated data. + + This function fits the normalized on-axis magnetic field of an ideal solenoid + to the provided data using a non-linear least squares optimization. It returns + the fitted solenoid parameters: radius, length, and the normalization factor (B0). + + Parameters + ---------- + z0 : np.ndarray + Array of axial positions (z) corresponding to the measured magnetic field data. + Bz0 : np.ndarray + Array of measured magnetic field values (Bz) at the positions in `z0`. + radius0 : float, optional + Initial guess for the radius of the solenoid (meters). Default is 0.1 m. + L0 : float, optional + Initial guess for the total length of the solenoid (meters). Default is 0.1 m. + + Returns + ------- + dict + Dictionary containing the fitted parameters: + - 'B0' : float + Normalization factor for the magnetic field (peak field at z = 0). + - 'radius' : float + Fitted solenoid radius (meters). + - 'L' : float + Fitted solenoid length (meters). + """ + + B0 = Bz0.max() + Bz_data = Bz0 / B0 + z_data = z0 - z0[np.argmax(Bz_data)] # Shift z so maximum is at z=0 + + # Fit the data to the model + popt, pcov = curve_fit(_fit_ideal_solenoid_model, z_data, Bz_data, p0=[radius0, L0]) + + # Extract fitted parameters + radius_fit, L_fit = popt + + d = {} + d["B0"] = B0 + d["radius"] = float(abs(radius_fit)) + d["L"] = float(abs(L_fit)) + + return d