{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# SPH Evaluation Kernels\n",
    "\n",
    "This tutorial shows how to use the SPH (Smoothed Particle Hydrodynamics) evaluation machinery in Struphy to reconstruct a smooth field from particle data, and to verify correctness by comparing against known exact solutions.\n",
    "\n",
    "The core idea: given $N$ particles with positions $\\boldsymbol{\\eta}_j \\in [0,1]^3$ and scalar weights $\\rho_j$, the SPH estimate of $\\rho$ at an arbitrary point $\\boldsymbol{\\eta}$ is\n",
    "\n",
    "$$\n",
    "\\rho(\\boldsymbol{\\eta}) \\approx \\sum_{j=0}^{N-1} \\rho_j\\, W_h(\\boldsymbol{\\eta} - \\boldsymbol{\\eta}_j)\\,,\n",
    "$$\n",
    "\n",
    "where $W_h$ is a smoothing kernel with bandwidth $h$. The gradient can be evaluated by differentiating the kernel:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial \\rho}{\\partial \\eta_k}(\\boldsymbol{\\eta}) \\approx \\sum_{j=0}^{N-1} \\rho_j\\,\n",
    "\\frac{\\partial W_h}{\\partial \\eta_k}(\\boldsymbol{\\eta} - \\boldsymbol{\\eta}_j)\\,.\n",
    "$$\n",
    "\n",
    "**What this tutorial covers**\n",
    "\n",
    "1. 1-D density reconstruction — value and first derivative, periodic and non-periodic boundary conditions, three kernel families.\n",
    "2. 2-D density reconstruction — the same workflow on a meshgrid, with a colour-plot comparison against the exact field.\n",
    "\n",
    "The evaluation is triggered via `ParticlesSPH.eval_density`, which internally dispatches to the box-based kernel in `struphy.pic.sph_eval_kernels` (an $\\mathcal{O}(1)$-per-point algorithm that restricts the kernel sum to the 27 neighbouring sorting boxes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from struphy import (\n",
    "    BoundaryParameters,\n",
    "    LoadingParameters,\n",
    "    SortingParameters,\n",
    "    domains,\n",
    "    perturbations,\n",
    ")\n",
    "from struphy.fields_background.equils import ConstantVelocity\n",
    "from struphy.pic.particles import ParticlesSPH"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "## Part 1 — 1-D density reconstruction\n",
    "\n",
    "### Problem setup\n",
    "\n",
    "We want to reconstruct the density field\n",
    "\n",
    "$$\n",
    "\\rho(\\eta_1) = 1.5 + \\cos(2\\pi\\eta_1), \\qquad \\eta_1 \\in [0, 1]\\,,\n",
    "$$\n",
    "\n",
    "from $N$ particles whose weights $\\rho_j$ are set by this exact function.\n",
    "\n",
    "The domain is a 3-D cuboid — the $\\eta_2, \\eta_3$ directions are kept trivial (`l2=r2`, `l3=r3` effectively, all action is in $\\eta_1$).  We use a **tesselation** loading strategy (particles placed on a regular lattice) so the reconstruction error is controlled by the lattice spacing rather than Monte Carlo noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cuboid domain: non-unit physical size, but the SPH evaluation is in logical space [0,1]^3\n",
    "domain_1d = domains.Cuboid(l1=1.0, r1=2.0, l2=10.0, r2=20.0, l3=100.0, r3=200.0)\n",
    "\n",
    "# Background equilibrium: constant density n0 = 1.5\n",
    "background_1d = ConstantVelocity(n=1.5, density_profile=\"constant\")\n",
    "background_1d.domain = domain_1d\n",
    "\n",
    "# Density perturbation: rho = n0 + cos(2 pi eta1)\n",
    "pert_1d = {\"n\": perturbations.ModesCos(ls=(1,), amps=(1.0,))}\n",
    "\n",
    "# Exact field and its eta1-derivative\n",
    "rho_exact   = lambda e1, e2, e3: 1.5 + np.cos(2 * np.pi * e1)\n",
    "drho_exact  = lambda e1, e2, e3: -2 * np.pi * np.sin(2 * np.pi * e1)\n",
    "\n",
    "# Sorting: 24 boxes in eta1, 1 in eta2 / eta3 — purely 1-D\n",
    "boxes_per_dim = (24, 1, 1)\n",
    "sorting_params = SortingParameters(boxes_per_dim=boxes_per_dim)\n",
    "\n",
    "# Kernel bandwidth = one box width in each dimension\n",
    "h1 = 1 / boxes_per_dim[0]\n",
    "h2 = 1 / boxes_per_dim[1]\n",
    "h3 = 1 / boxes_per_dim[2]\n",
    "\n",
    "# Evaluation grid (purely 1-D meshgrid)\n",
    "n_eval = 200\n",
    "eta1_pts = np.linspace(0, 1, n_eval)\n",
    "eta2_pts = np.array([0.0])\n",
    "eta3_pts = np.array([0.0])\n",
    "ee1, ee2, ee3 = np.meshgrid(eta1_pts, eta2_pts, eta3_pts, indexing=\"ij\")\n",
    "\n",
    "print(f\"Domain logical size: [0,1]^3   |   Sorting: {boxes_per_dim}\")\n",
    "print(f\"Kernel bandwidth h1={h1:.4f}, h2={h2:.4f}, h3={h3:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "### 1.1 Periodic boundary condition\n",
    "\n",
    "With a periodic BC the boundary layer receives contributions from particles on the opposite side of the domain. The mirror images are handled automatically inside the distance kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_particles_1d(bc_x, ppb=4, loading=\"tesselation\"):\n",
    "    \"\"\"Construct and initialise a 1-D ParticlesSPH object.\"\"\"\n",
    "    loading_params  = LoadingParameters(ppb=ppb, seed=1607, loading=loading)\n",
    "    boundary_params = BoundaryParameters(bc_sph=(bc_x, \"periodic\", \"periodic\"))\n",
    "    particles = ParticlesSPH(\n",
    "        comm_world=None,\n",
    "        loading_params=loading_params,\n",
    "        boundary_params=boundary_params,\n",
    "        sorting_params=sorting_params,\n",
    "        bufsize=1.0,\n",
    "        domain=domain_1d,\n",
    "        background=background_1d,\n",
    "        perturbations=pert_1d,\n",
    "        n_as_volume_form=True,\n",
    "    )\n",
    "    particles.draw_markers(sort=False)\n",
    "    particles.initialize_weights()\n",
    "    return particles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "particles_periodic = make_particles_1d(ppb=4, bc_x=\"periodic\")\n",
    "\n",
    "kernels_1d = [\"gaussian_1d\", \"linear_1d\", \"trigonometric_1d\"]\n",
    "x_plot = ee1.squeeze()\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4), sharey=True)\n",
    "fig.suptitle(r\"1-D density reconstruction  ($\\rho = 1.5 + \\cos(2\\pi\\eta_1)$, periodic BC)\")\n",
    "\n",
    "for ax, kernel in zip(axes, kernels_1d):\n",
    "    rho_sph = particles_periodic.eval_density(\n",
    "        ee1, ee2, ee3,\n",
    "        h1=h1, h2=h2, h3=h3,\n",
    "        kernel_type=kernel,\n",
    "        derivative=0,\n",
    "    ).squeeze()\n",
    "    rho_ex = rho_exact(x_plot, 0.0, 0.0)\n",
    "    err = np.max(np.abs(rho_sph - rho_ex)) / np.max(np.abs(rho_ex))\n",
    "\n",
    "    ax.plot(x_plot, rho_ex, \"k-\",  lw=2, label=\"exact\")\n",
    "    ax.plot(x_plot, rho_sph, \"r--\", lw=1.5, label=f\"SPH  (err={err:.2e})\")\n",
    "    ax.set_title(kernel)\n",
    "    ax.set_xlabel(r\"$\\eta_1$\")\n",
    "    ax.legend(fontsize=8)\n",
    "    ax.grid(True, alpha=0.3)\n",
    "\n",
    "axes[0].set_ylabel(r\"$\\rho$\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "### 1.2 Derivative evaluation\n",
    "\n",
    "Setting `derivative=1` returns the $\\partial/\\partial\\eta_1$ component of the gradient instead of the field value. The trigonometric kernel is exact for a single cosine mode; the linear kernel requires more particles per box to reach the same accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# More particles per box to resolve the derivative accurately\n",
    "particles_deriv = make_particles_1d(bc_x=\"periodic\", ppb=100)\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4), sharey=True)\n",
    "fig.suptitle(\n",
    "    r\"1-D derivative reconstruction  ($\\partial\\rho/\\partial\\eta_1 = -2\\pi\\sin(2\\pi\\eta_1)$, periodic BC)\"\n",
    ")\n",
    "\n",
    "for ax, kernel in zip(axes, kernels_1d):\n",
    "    drho_sph = particles_deriv.eval_density(\n",
    "        ee1, ee2, ee3,\n",
    "        h1=h1, h2=h2, h3=h3,\n",
    "        kernel_type=kernel,\n",
    "        derivative=1,\n",
    "    ).squeeze()\n",
    "    drho_ex = drho_exact(x_plot, 0.0, 0.0)\n",
    "    err = np.max(np.abs(drho_sph - drho_ex)) / np.max(np.abs(drho_ex))\n",
    "\n",
    "    ax.plot(x_plot, drho_ex, \"k-\",  lw=2, label=\"exact\")\n",
    "    ax.plot(x_plot, drho_sph, \"r--\", lw=1.5, label=f\"SPH  (err={err:.2e})\")\n",
    "    ax.set_title(kernel)\n",
    "    ax.set_xlabel(r\"$\\eta_1$\")\n",
    "    ax.legend(fontsize=8)\n",
    "    ax.grid(True, alpha=0.3)\n",
    "\n",
    "axes[0].set_ylabel(r\"$\\partial\\rho/\\partial\\eta_1$\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "### 1.3 Effect of boundary conditions\n",
    "\n",
    "Struphy supports three SPH boundary conditions:\n",
    "\n",
    "| BC | Description |\n",
    "|---|---|\n",
    "| `\"periodic\"` | Mirror images placed across the periodic boundary. |\n",
    "| `\"mirror\"` | Ghost particles reflected at the wall; enforces zero normal gradient. |\n",
    "| `\"fixed\"` | Ghost particles reflected and negated; enforces zero field value at the wall. |\n",
    "\n",
    "Below we compare all three for the density value.  Note that `mirror` and `fixed` are only meaningful near the boundary; the interior is identical to `periodic`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary_conditions = [\"periodic\", \"mirror\", \"fixed\"]\n",
    "colors = [\"tab:blue\", \"tab:orange\", \"tab:green\"]\n",
    "kernel = \"gaussian_1d\"\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 4))\n",
    "ax.plot(x_plot, rho_exact(x_plot, 0.0, 0.0), \"k-\", lw=1.5, label=\"exact\", zorder=5)\n",
    "\n",
    "for bc, color in zip(boundary_conditions, colors):\n",
    "    p = make_particles_1d(bc_x=bc)\n",
    "    rho_sph = p.eval_density(\n",
    "        ee1, ee2, ee3,\n",
    "        h1=h1, h2=h2, h3=h3,\n",
    "        kernel_type=kernel,\n",
    "        derivative=0,\n",
    "    ).squeeze()\n",
    "    ax.plot(x_plot, rho_sph, \"--\", color=color, lw=2.5, label=f\"bc={bc!r}\")\n",
    "\n",
    "ax.set_title(f\"Boundary condition comparison ({kernel})\")\n",
    "ax.set_xlabel(r\"$\\eta_1$\")\n",
    "ax.set_ylabel(r\"$\\rho$\")\n",
    "ax.legend()\n",
    "ax.grid(True, alpha=0.3)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "### 1.4 Tesselation vs Monte-Carlo loading\n",
    "\n",
    "With **tesselation** loading, particles are placed on a regular lattice (`ppb` per sorting box), giving a smooth, deterministic reconstruction.  With **pseudo-random** loading (Monte Carlo), particle positions are drawn from a uniform distribution, which introduces sampling noise that decays only as $\\mathcal{O}(N^{-1/2})$.\n",
    "\n",
    "Here we compare both strategies at matched particle counts:\n",
    "\n",
    "| Strategy | `ppb` | $N$ (approx.) |\n",
    "|---|---|---|\n",
    "| tesselation | 4 | 96 |\n",
    "| pseudo_random | 40 | 960 |\n",
    "\n",
    "The pseudo-random run uses **ten times** as many particles yet still shows visible fluctuations, illustrating why tesselation is preferred whenever the initial particle positions can be chosen freely."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "ppb_tess = 4\n",
    "ppb_rand = ppb_tess * 10  # ten times more particles\n",
    "\n",
    "particles_tess = make_particles_1d(bc_x=\"periodic\", ppb=ppb_tess, loading=\"tesselation\")\n",
    "particles_rand = make_particles_1d(bc_x=\"periodic\", ppb=ppb_rand, loading=\"pseudo_random\")\n",
    "\n",
    "Np_tess = particles_tess.Np\n",
    "Np_rand = particles_rand.Np\n",
    "print(f\"Tesselation  ppb={ppb_tess}  →  Np={Np_tess}\")\n",
    "print(f\"Pseudo-random  ppb={ppb_rand}  →  Np={Np_rand}\")\n",
    "\n",
    "rho_ex_1d = rho_exact(x_plot, 0.0, 0.0)\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4), sharey=True)\n",
    "fig.suptitle(\n",
    "    r\"Tesselation vs Monte Carlo — 1-D density  ($\\rho = 1.5 + \\cos(2\\pi\\eta_1)$, periodic BC)\"\n",
    ")\n",
    "\n",
    "for ax, kernel in zip(axes, kernels_1d):\n",
    "    rho_tess = particles_tess.eval_density(\n",
    "        ee1, ee2, ee3, h1=h1, h2=h2, h3=h3,\n",
    "        kernel_type=kernel, derivative=0,\n",
    "    ).squeeze()\n",
    "    rho_rand = particles_rand.eval_density(\n",
    "        ee1, ee2, ee3, h1=h1, h2=h2, h3=h3,\n",
    "        kernel_type=kernel, derivative=0,\n",
    "    ).squeeze()\n",
    "\n",
    "    err_tess = np.max(np.abs(rho_tess - rho_ex_1d)) / np.max(np.abs(rho_ex_1d))\n",
    "    err_rand = np.max(np.abs(rho_rand - rho_ex_1d)) / np.max(np.abs(rho_ex_1d))\n",
    "\n",
    "    ax.plot(x_plot, rho_ex_1d, \"k-\", lw=2, label=\"exact\")\n",
    "    ax.plot(x_plot, rho_tess, \"b-\",  lw=1.5,\n",
    "            label=f\"tesselation  Np={Np_tess}  err={err_tess:.2e}\")\n",
    "    ax.plot(x_plot, rho_rand, \"r-\",  lw=1,\n",
    "            label=f\"pseudo_random  Np={Np_rand}  err={err_rand:.2e}\")\n",
    "    ax.set_title(kernel)\n",
    "    ax.set_xlabel(r\"$\\eta_1$\")\n",
    "    ax.legend(fontsize=7)\n",
    "    ax.grid(True, alpha=0.3)\n",
    "\n",
    "axes[0].set_ylabel(r\"$\\rho$\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "## Part 2 — 2-D density reconstruction\n",
    "\n",
    "### Problem setup\n",
    "\n",
    "We extend the test to two dimensions and reconstruct\n",
    "\n",
    "$$\n",
    "\\rho(\\eta_1, \\eta_2) = 1.5 + \\cos(2\\pi\\eta_1)\\cos(2\\pi\\eta_2)\\,,\n",
    "$$\n",
    "\n",
    "together with its partial derivatives\n",
    "\n",
    "$$\n",
    "\\frac{\\partial\\rho}{\\partial\\eta_1} = -2\\pi\\sin(2\\pi\\eta_1)\\cos(2\\pi\\eta_2)\\,,\\qquad\n",
    "\\frac{\\partial\\rho}{\\partial\\eta_2} = -2\\pi\\cos(2\\pi\\eta_1)\\sin(2\\pi\\eta_2)\\,.\n",
    "$$\n",
    "\n",
    "The sorting uses a $12\\times 12\\times 1$ box grid, which matches the `trigonometric_2d`, `gaussian_2d`, and `linear_2d` kernel families."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "domain_2d = domains.Cuboid(l1=1.0, r1=2.0, l2=0.0, r2=2.0, l3=100.0, r3=200.0)\n",
    "\n",
    "background_2d = ConstantVelocity(n=1.5, density_profile=\"constant\")\n",
    "background_2d.domain = domain_2d\n",
    "\n",
    "pert_2d = {\"n\": perturbations.ModesCosCos(ls=(1,), ms=(1,), amps=(1.0,))}\n",
    "\n",
    "rho_2d_exact   = lambda e1, e2, e3: 1.5 + np.cos(2*np.pi*e1) * np.cos(2*np.pi*e2)\n",
    "drho_2d_deta1  = lambda e1, e2, e3: -2*np.pi * np.sin(2*np.pi*e1) * np.cos(2*np.pi*e2)\n",
    "drho_2d_deta2  = lambda e1, e2, e3: -2*np.pi * np.cos(2*np.pi*e1) * np.sin(2*np.pi*e2)\n",
    "\n",
    "boxes_per_dim_2d = (12, 12, 1)\n",
    "h1_2d = 1 / boxes_per_dim_2d[0]\n",
    "h2_2d = 1 / boxes_per_dim_2d[1]\n",
    "h3_2d = 1 / boxes_per_dim_2d[2]\n",
    "\n",
    "# Tesselation loading: 16 particles per box\n",
    "loading_params_2d  = LoadingParameters(ppb=16, loading=\"tesselation\")\n",
    "boundary_params_2d = BoundaryParameters(bc_sph=(\"periodic\", \"periodic\", \"periodic\"))\n",
    "sorting_params_2d  = SortingParameters(boxes_per_dim=boxes_per_dim_2d)\n",
    "\n",
    "particles_2d = ParticlesSPH(\n",
    "    comm_world=None,\n",
    "    loading_params=loading_params_2d,\n",
    "    boundary_params=boundary_params_2d,\n",
    "    sorting_params=sorting_params_2d,\n",
    "    bufsize=1.0,\n",
    "    domain=domain_2d,\n",
    "    background=background_2d,\n",
    "    perturbations=pert_2d,\n",
    "    n_as_volume_form=True,\n",
    ")\n",
    "particles_2d.draw_markers(sort=False)\n",
    "particles_2d.initialize_weights()\n",
    "\n",
    "# 2-D evaluation meshgrid\n",
    "n_eval_2d = 50\n",
    "eta1_2d = np.linspace(0, 1.0, n_eval_2d)\n",
    "eta2_2d = np.linspace(0, 1.0, n_eval_2d)\n",
    "eta3_2d = np.array([0.0])\n",
    "ee1_2d, ee2_2d, ee3_2d = np.meshgrid(eta1_2d, eta2_2d, eta3_2d, indexing=\"ij\")\n",
    "\n",
    "print(f\"Particles: {particles_2d.Np}   |   Boxes: {boxes_per_dim_2d}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "### 2.1 Density value\n",
    "\n",
    "We evaluate the density on the 2-D meshgrid and compare side-by-side with the exact field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel_2d = \"gaussian_2d\"\n",
    "\n",
    "rho_sph_2d = particles_2d.eval_density(\n",
    "    ee1_2d, ee2_2d, ee3_2d,\n",
    "    h1=h1_2d, h2=h2_2d, h3=h3_2d,\n",
    "    kernel_type=kernel_2d,\n",
    "    derivative=0,\n",
    ").squeeze()\n",
    "\n",
    "rho_ex_2d = rho_2d_exact(ee1_2d, ee2_2d, ee3_2d).squeeze()\n",
    "err_2d = np.max(np.abs(rho_sph_2d - rho_ex_2d)) / np.max(np.abs(rho_ex_2d))\n",
    "print(f\"Max relative error: {err_2d:.3e}\")\n",
    "\n",
    "x_plot_2d = ee1_2d.squeeze()\n",
    "y_plot_2d = ee2_2d.squeeze()\n",
    "vmin = rho_ex_2d.min()\n",
    "vmax = rho_ex_2d.max()\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(16, 4))\n",
    "fig.suptitle(f\"2-D density reconstruction  ({kernel_2d}, periodic BC)\")\n",
    "\n",
    "im0 = axes[0].pcolormesh(x_plot_2d, y_plot_2d, rho_ex_2d, vmin=vmin, vmax=vmax, shading=\"auto\")\n",
    "axes[0].set_title(\"Exact\")\n",
    "fig.colorbar(im0, ax=axes[0])\n",
    "\n",
    "im1 = axes[1].pcolormesh(x_plot_2d, y_plot_2d, rho_sph_2d, vmin=vmin, vmax=vmax, shading=\"auto\")\n",
    "axes[1].set_title(\"SPH\")\n",
    "fig.colorbar(im1, ax=axes[1])\n",
    "\n",
    "err_field = np.abs(rho_sph_2d - rho_ex_2d)\n",
    "im2 = axes[2].pcolormesh(x_plot_2d, y_plot_2d, err_field, shading=\"auto\", cmap=\"Reds\")\n",
    "axes[2].set_title(f\"Point-wise error  (max={err_2d:.2e})\")\n",
    "fig.colorbar(im2, ax=axes[2])\n",
    "\n",
    "for ax in axes:\n",
    "    ax.set_xlabel(r\"$\\eta_1$\")\n",
    "    ax.set_ylabel(r\"$\\eta_2$\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17",
   "metadata": {},
   "source": [
    "### 2.2 Partial derivatives\n",
    "\n",
    "`derivative=1` gives $\\partial\\rho/\\partial\\eta_1$ and `derivative=2` gives $\\partial\\rho/\\partial\\eta_2$.  We show both below using the trigonometric kernel, which is spectrally exact for smooth periodic functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# More particles per box for the trigonometric kernel derivative test\n",
    "loading_params_deriv = LoadingParameters(ppb=100, loading=\"tesselation\")\n",
    "particles_2d_deriv = ParticlesSPH(\n",
    "    comm_world=None,\n",
    "    loading_params=loading_params_deriv,\n",
    "    boundary_params=boundary_params_2d,\n",
    "    sorting_params=sorting_params_2d,\n",
    "    bufsize=1.0,\n",
    "    domain=domain_2d,\n",
    "    background=background_2d,\n",
    "    perturbations=pert_2d,\n",
    "    n_as_volume_form=True,\n",
    ")\n",
    "particles_2d_deriv.draw_markers(sort=False)\n",
    "particles_2d_deriv.initialize_weights()\n",
    "\n",
    "kernel_deriv = \"trigonometric_2d\"\n",
    "derivative_info = [\n",
    "    (1, drho_2d_deta1, r\"$\\partial\\rho/\\partial\\eta_1$\"),\n",
    "    (2, drho_2d_deta2, r\"$\\partial\\rho/\\partial\\eta_2$\"),\n",
    "]\n",
    "\n",
    "fig, axes = plt.subplots(2, 3, figsize=(16, 9))\n",
    "fig.suptitle(f\"2-D derivative reconstruction  ({kernel_deriv}, periodic BC)\")\n",
    "\n",
    "for row, (deriv_idx, exact_fn, label) in enumerate(derivative_info):\n",
    "    drho_sph = particles_2d_deriv.eval_density(\n",
    "        ee1_2d, ee2_2d, ee3_2d,\n",
    "        h1=h1_2d, h2=h2_2d, h3=h3_2d,\n",
    "        kernel_type=kernel_deriv,\n",
    "        derivative=deriv_idx,\n",
    "    ).squeeze()\n",
    "    drho_ex = exact_fn(ee1_2d, ee2_2d, ee3_2d).squeeze()\n",
    "    err = np.max(np.abs(drho_sph - drho_ex)) / np.max(np.abs(drho_ex))\n",
    "\n",
    "    vmin_d, vmax_d = drho_ex.min(), drho_ex.max()\n",
    "\n",
    "    im0 = axes[row, 0].pcolormesh(x_plot_2d, y_plot_2d, drho_ex,  vmin=vmin_d, vmax=vmax_d, shading=\"auto\")\n",
    "    axes[row, 0].set_title(f\"Exact {label}\")\n",
    "    fig.colorbar(im0, ax=axes[row, 0])\n",
    "\n",
    "    im1 = axes[row, 1].pcolormesh(x_plot_2d, y_plot_2d, drho_sph, vmin=vmin_d, vmax=vmax_d, shading=\"auto\")\n",
    "    axes[row, 1].set_title(f\"SPH {label}\")\n",
    "    fig.colorbar(im1, ax=axes[row, 1])\n",
    "\n",
    "    err_field = np.abs(drho_sph - drho_ex)\n",
    "    im2 = axes[row, 2].pcolormesh(x_plot_2d, y_plot_2d, err_field, shading=\"auto\", cmap=\"Reds\")\n",
    "    axes[row, 2].set_title(f\"Error  (max={err:.2e})\")\n",
    "    fig.colorbar(im2, ax=axes[row, 2])\n",
    "\n",
    "for ax in axes.flat:\n",
    "    ax.set_xlabel(r\"$\\eta_1$\")\n",
    "    ax.set_ylabel(r\"$\\eta_2$\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19",
   "metadata": {},
   "source": [
    "### 2.3 Kernel comparison in 2-D\n",
    "\n",
    "All three 2-D kernel families are compared for the density value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "kernels_2d = [\"gaussian_2d\", \"linear_2d\", \"trigonometric_2d\"]\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(16, 4))\n",
    "fig.suptitle(r\"Kernel comparison — 2-D density  ($\\rho = 1.5 + \\cos(2\\pi\\eta_1)\\cos(2\\pi\\eta_2)$)\")\n",
    "\n",
    "for ax, kernel in zip(axes, kernels_2d):\n",
    "    rho_sph = particles_2d.eval_density(\n",
    "        ee1_2d, ee2_2d, ee3_2d,\n",
    "        h1=h1_2d, h2=h2_2d, h3=h3_2d,\n",
    "        kernel_type=kernel,\n",
    "        derivative=0,\n",
    "    ).squeeze()\n",
    "    rho_ex = rho_2d_exact(ee1_2d, ee2_2d, ee3_2d).squeeze()\n",
    "    err = np.max(np.abs(rho_sph - rho_ex)) / np.max(np.abs(rho_ex))\n",
    "\n",
    "    im = ax.pcolormesh(\n",
    "        x_plot_2d, y_plot_2d, np.abs(rho_sph - rho_ex),\n",
    "        shading=\"auto\", cmap=\"Reds\",\n",
    "    )\n",
    "    ax.set_title(f\"{kernel}\\nerr={err:.2e}\")\n",
    "    ax.set_xlabel(r\"$\\eta_1$\")\n",
    "    ax.set_ylabel(r\"$\\eta_2$\")\n",
    "    fig.colorbar(im, ax=ax)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "| Task | Key parameter |\n",
    "|---|---|\n",
    "| Field value | `derivative=0` |\n",
    "| $\\partial/\\partial\\eta_1$ | `derivative=1` |\n",
    "| $\\partial/\\partial\\eta_2$ | `derivative=2` |\n",
    "| $\\partial/\\partial\\eta_3$ | `derivative=3` |\n",
    "| Kernel family | `kernel_type` — e.g. `\"gaussian_1d\"`, `\"trigonometric_2d\"`, `\"linear_3d\"` |\n",
    "| Kernel bandwidth | `h1`, `h2`, `h3` — typically set to `1/boxes_per_dim` |\n",
    "\n",
    "The error decreases as `ppb` (particles per box) increases and as the kernel bandwidth matches the inter-particle spacing.  The **trigonometric** kernel achieves spectral accuracy for smooth periodic functions; the **Gaussian** and **linear** kernels are more robust near non-periodic boundaries.\n",
    "\n",
    "For implementation details of the underlying kernel functions see `struphy.pic.sph_eval_kernels` and `struphy.pic.sph_smoothing_kernels`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env (3.12.3)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}