Fixed convolve_density in impedance model

Author: ChristopherMayes

docs/examples/wakefields/resistive_wall_impedance.ipynb

@@ -179,7 +179,8 @@
    "outputs": [],
    "source": [
     "# Compute wakefield arrays (slow due to numerical integration)\n",
-    "zs = np.linspace(0, 200e-6, 30)\n",
+    "# Wake is defined for z < 0 (behind the source particle)\n",
+    "zs = -np.linspace(1e-6, 200e-6, 30)  # negative z values\n",
     "Wz_flat = imp_flat.wakefield(zs, k_max=1e6)\n",
     "Wz_round = imp_round.wakefield(zs, k_max=1e6)"
    ]
@@ -192,9 +193,9 @@
    "outputs": [],
    "source": [
     "fig, ax = plt.subplots(figsize=(8, 4))\n",
-    "ax.plot(zs * 1e6, Wz_flat * 1e-12, label=\"Flat\")\n",
-    "ax.plot(zs * 1e6, Wz_round * 1e-12, \"--\", label=\"Round\")\n",
-    "ax.set_xlabel(r\"$z$ (µm)\")\n",
+    "ax.plot(-zs * 1e6, Wz_flat * 1e-12, label=\"Flat\")  # plot as distance behind source\n",
+    "ax.plot(-zs * 1e6, Wz_round * 1e-12, \"--\", label=\"Round\")\n",
+    "ax.set_xlabel(r\"Distance behind source $|z|$ (µm)\")\n",
     "ax.set_ylabel(r\"$W(z)$ (V/pC/m)\")\n",
     "ax.legend()\n",
     "ax.set_title(\"Wakefield Comparison: Flat vs Round Geometry\")"
@@ -251,11 +252,10 @@
    "outputs": [],
    "source": [
     "# Evaluate pseudomode wakefields (very fast!)\n",
-    "zs_fine = np.linspace(0, 200e-6, 200)\n",
-    "\n",
-    "# Note: pseudomode expects negative z (trailing the source)\n",
-    "Wz_pseudo_flat = wake_flat.pseudomode(-zs_fine)\n",
-    "Wz_pseudo_round = wake_round.pseudomode(-zs_fine)"
+    "# Both methods use z < 0 for trailing particles\n",
+    "zs_fine = -np.linspace(1e-6, 200e-6, 200)  # negative z values\n",
+    "Wz_pseudo_flat = wake_flat.pseudomode(zs_fine)\n",
+    "Wz_pseudo_round = wake_round.pseudomode(zs_fine)"
    ]
   },
   {
@@ -269,17 +269,17 @@
     "fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n",
     "\n",
     "ax = axes[0]\n",
-    "ax.plot(zs * 1e6, Wz_flat * 1e-12, \"o\", label=\"Numerical\", markersize=5)\n",
-    "ax.plot(zs_fine * 1e6, Wz_pseudo_flat * 1e-12, \"-\", label=\"Pseudomode\", alpha=0.8)\n",
-    "ax.set_xlabel(r\"$z$ (µm)\")\n",
+    "ax.plot(-zs * 1e6, Wz_flat * 1e-12, \"o\", label=\"Numerical\", markersize=5)\n",
+    "ax.plot(-zs_fine * 1e6, Wz_pseudo_flat * 1e-12, \"-\", label=\"Pseudomode\", alpha=0.8)\n",
+    "ax.set_xlabel(r\"Distance behind source $|z|$ (µm)\")\n",
     "ax.set_ylabel(r\"$W(z)$ (V/pC/m)\")\n",
     "ax.legend()\n",
     "ax.set_title(\"Flat Geometry\")\n",
     "\n",
     "ax = axes[1]\n",
-    "ax.plot(zs * 1e6, Wz_round * 1e-12, \"o\", label=\"Numerical\", markersize=5)\n",
-    "ax.plot(zs_fine * 1e6, Wz_pseudo_round * 1e-12, \"-\", label=\"Pseudomode\", alpha=0.8)\n",
-    "ax.set_xlabel(r\"$z$ (µm)\")\n",
+    "ax.plot(-zs * 1e6, Wz_round * 1e-12, \"o\", label=\"Numerical\", markersize=5)\n",
+    "ax.plot(-zs_fine * 1e6, Wz_pseudo_round * 1e-12, \"-\", label=\"Pseudomode\", alpha=0.8)\n",
+    "ax.set_xlabel(r\"Distance behind source $|z|$ (µm)\")\n",
     "ax.set_ylabel(r\"$W(z)$ (V/pC/m)\")\n",
     "ax.legend()\n",
     "ax.set_title(\"Round Geometry\")\n",
@@ -307,7 +307,7 @@
    "source": [
     "%%timeit -n 1 -r 3\n",
     "# Numerical integration for 20 points\n",
-    "z_test = np.linspace(0, 100e-6, 20)\n",
+    "z_test = -np.linspace(1e-6, 100e-6, 20)\n",
     "_ = imp_flat.wakefield(z_test, k_max=1e6)"
    ]
   },
@@ -320,8 +320,8 @@
    "source": [
     "%%timeit -n 1000 -r 3\n",
     "# Pseudomode for 1000 points\n",
-    "z_test = np.linspace(0, 100e-6, 1000)\n",
-    "_ = wake_flat.pseudomode(-z_test)"
+    "z_test = -np.linspace(1e-6, 100e-6, 1000)\n",
+    "_ = wake_flat.pseudomode(z_test)"
    ]
   },
   {
@@ -341,6 +341,18 @@
     "- Working with unusual parameter regimes\n",
     "- Need access to the impedance spectrum $Z(k)$"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "592c26c6",
+   "metadata": {},
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b76be5fd",
+   "metadata": {},
+   "source": []
   }
  ],
  "metadata": {

docs/examples/wakefields/resistive_wall_wakefield.ipynb

@@ -18,6 +18,17 @@
     "- Bane, Stupakov, Tu, [EPAC 2006 THPCH073](https://accelconf.web.cern.ch/e06/PAPERS/THPCH073.PDF)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dbaddcb4-a6ff-42ac-b55a-4395b4cb9ef0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%load_ext autoreload\n",
+    "%autoreload 2"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -262,6 +273,75 @@
     "P.drift_to_z()\n",
     "P.wakefield_plot(wake_ref, figsize=(12, 4))"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d867b96b-4047-405a-85ca-008185bb8d11",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "P.drift_to_t()\n",
+    "H, b = P.histogramdd(\"z\", bins=100)\n",
+    "dz = np.diff(b[0])[0]\n",
+    "density = H / dz\n",
+    "n = len(H)\n",
+    "z = np.arange(0, n, 1) * dz\n",
+    "\n",
+    "v = wake.convolve_density(density, dz=dz)\n",
+    "plt.plot(z, v)\n",
+    "\n",
+    "wp = -wake.pseudomode.particle_kicks(P.z, P.weight)\n",
+    "zp = P.z - P.z.min()\n",
+    "plt.scatter(zp, wp)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fdbe6d76-b402-4e5f-ae9c-fa78a3a0456e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pmd_beamphysics.wakefields import ResistiveWallImpedance"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9834d67e-b2cb-4312-b71e-453dcbb4e99a",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2b8e3018-d376-45c1-a77c-671eb1153c22",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "imp = ResistiveWallImpedance(\n",
+    "    radius=0.0025,\n",
+    "    conductivity=6.5e7,  # σ [S/m]\n",
+    "    relaxation_time=27e-15,  # τ [s]\n",
+    "    geometry=\"round\",\n",
+    ")\n",
+    "\n",
+    "wi = imp.convolve_density(density, dz=dz)\n",
+    "plt.plot(z, wi)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c4ab6d5c-9ff5-4a83-b119-0b01710b0134",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.plot(z, wi)\n",
+    "plt.plot(z, v)"
+   ]
   }
  ],
  "metadata": {

pmd_beamphysics/wakefields/resistive_wall.py

@@ -729,10 +729,10 @@ def convolve_density(self, density, dz, offset=0):
         ----------
 
         density: ndarray
-            charge density array
+            charge density array in C / m
 
         dz: float
-            array spacing
+            array spacing in m
 
         offset : float
             Offset coordinates for the center of the grid in [m]. Default: 0