44 "cell_type": "markdown",
55 "id": "f49966c6",
66 "metadata": {},
7- "source": [
8- "# scattnlay vs miepython Near-Field Validation\n",
9- "\n",
10- "This notebook validates `miepython` near-field electric and magnetic fields against\n",
11- "precomputed references from **scattnlay**:\n",
12- "\n",
13- "- Project: [scattnlay](https://github.com/ovidiopr/scattnlay)\n",
14- "- Implementation: C++ (with Python bindings / CLI)\n",
15- "- Scope: electromagnetic scattering and near fields for multilayered spheres\n",
16- "\n",
17- "Reference arrays are generated once by:\n",
18- "\n",
19- "`python generate_scattnlay_reference_fields.py`\n",
20- "\n",
21- "and stored in `docs/data/` as `.npy` files.\n"
22- ]
7+ "source": "# scattnlay vs miepython Near-Field Validation\n\nThis notebook validates `miepython` near-field electric and magnetic fields against\nprecomputed references from **scattnlay**:\n\n- Project: [scattnlay](https://github.com/ovidiopr/scattnlay)\n- Implementation: C++ (with Python bindings / CLI)\n- Scope: electromagnetic scattering and near fields for multilayered spheres\n\nReference arrays are generated by:\n\n`python docs/generate_scattnlay_reference_fields.py`\n\nand stored in `docs/data/` as `.npy` files.\n"
238 },
249 {
2510 "cell_type": "markdown",
5136 "id": "6fc2d052",
5237 "metadata": {},
5338 "outputs": [],
54- "source": [
55- "from pathlib import Path\n",
56- "\n",
57- "import matplotlib.pyplot as plt\n",
58- "import numpy as np\n",
59- "\n",
60- "from miepython.field import eh_near_cartesian\n",
61- "\n",
62- "\n",
63- "def find_data_dir(default=Path(\"data\")):\n",
64- " candidates = [default, Path(\"docs\") / default]\n",
65- " for candidate in candidates:\n",
66- " if candidate.exists():\n",
67- " return candidate.resolve()\n",
68- " raise FileNotFoundError(\n",
69- " \"Could not find reference data directory. \"\n",
70- " \"Run: python docs/generate_scattnlay_reference_fields.py\"\n",
71- " )\n",
72- "\n",
73- "\n",
74- "def load_scattnlay_case_from_npy(case_key, data_dir=None):\n",
75- " base = find_data_dir() if data_dir is None else Path(data_dir).resolve()\n",
76- " paths = {\n",
77- " \"X\": base / f\"scattnlay_{case_key}_X.npy\",\n",
78- " \"Z\": base / f\"scattnlay_{case_key}_Z.npy\",\n",
79- " \"E\": base / f\"scattnlay_{case_key}_E.npy\",\n",
80- " \"H\": base / f\"scattnlay_{case_key}_H.npy\",\n",
81- " }\n",
82- "\n",
83- " missing = [str(path) for path in paths.values() if not path.exists()]\n",
84- " if missing:\n",
85- " msg = \"Missing precomputed scattnlay files:\\n\" + \"\\n\".join(missing)\n",
86- " msg += \"\\nRun: python docs/generate_scattnlay_reference_fields.py\"\n",
87- " raise FileNotFoundError(msg)\n",
88- "\n",
89- " X = np.load(paths[\"X\"])\n",
90- " Z = np.load(paths[\"Z\"])\n",
91- " E = np.load(paths[\"E\"])\n",
92- " H = np.load(paths[\"H\"])\n",
93- " return X, Z, E, H\n",
94- "\n",
95- "\n",
96- "def summarize_errors(label, Eabs_s, Eabs_m, Habs_s, Habs_m):\n",
97- " E_rel_err = np.abs(Eabs_m - Eabs_s) / np.maximum(Eabs_s, 1e-12)\n",
98- " H_rel_err = np.abs(Habs_m - Habs_s) / np.maximum(Habs_s, 1e-12)\n",
99- "\n",
100- " print(f\"[{label}] scattnlay max |E| = {np.max(Eabs_s):.6g}\")\n",
101- " print(f\"[{label}] miepython max |E| = {np.max(Eabs_m):.6g}\")\n",
102- " print(f\"[{label}] scattnlay max |H| (eta-scaled) = {np.max(Habs_s):.6g}\")\n",
103- " print(f\"[{label}] miepython max |H| = {np.max(Habs_m):.6g}\")\n",
104- " print(f\"[{label}] E relative error median/max = {np.median(E_rel_err):.3e} / {np.max(E_rel_err):.3e}\")\n",
105- " print(f\"[{label}] H relative error median/max = {np.median(H_rel_err):.3e} / {np.max(H_rel_err):.3e}\")\n",
106- "\n",
107- " return E_rel_err, H_rel_err\n",
108- "\n",
109- "\n",
110- "def plot_comparison(title, Xg, Zg, radius, Eabs_s, Eabs_m, E_rel_err, Habs_s, Habs_m, H_rel_err):\n",
111- " fig, axes = plt.subplots(2, 3, figsize=(13, 8), constrained_layout=True)\n",
112- "\n",
113- " plots = [\n",
114- " (Eabs_s, \"scattnlay |E|\", \"turbo\"),\n",
115- " (Eabs_m, \"miepython |E|\", \"turbo\"),\n",
116- " (E_rel_err, \"|E| relative error\", \"viridis\"),\n",
117- " (Habs_s, \"scattnlay |H| (eta-scaled)\", \"magma\"),\n",
118- " (Habs_m, \"miepython |H|\", \"magma\"),\n",
119- " (H_rel_err, \"|H| relative error\", \"viridis\"),\n",
120- " ]\n",
121- "\n",
122- " for ax, (arr, ttl, cmap) in zip(axes.flat, plots):\n",
123- " im = ax.imshow(\n",
124- " arr.T,\n",
125- " origin=\"lower\",\n",
126- " extent=(Xg.min(), Xg.max(), Zg.min(), Zg.max()),\n",
127- " cmap=cmap,\n",
128- " aspect=\"equal\",\n",
129- " )\n",
130- " circle_color = \"w\" if cmap in (\"magma\", \"viridis\") else \"k\"\n",
131- " ax.add_patch(plt.Circle((0.0, 0.0), radius, fill=False, color=circle_color, lw=1.1))\n",
132- " ax.set_title(ttl)\n",
133- " ax.set_xlabel(\"x\")\n",
134- " ax.set_ylabel(\"z\")\n",
135- " fig.colorbar(im, ax=ax, fraction=0.046, pad=0.02)\n",
136- "\n",
137- " fig.suptitle(title, y=1.02)\n"
138- ]
39+ "source": "from pathlib import Path\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfrom miepython.field import eh_near_cartesian\n\n\ndef find_data_dir(default=Path(\"data\")):\n candidates = [default, Path(\"docs\") / default]\n for candidate in candidates:\n if candidate.exists():\n return candidate.resolve()\n raise FileNotFoundError(\n \"Could not find reference data directory. \"\n \"Run: python docs/generate_scattnlay_reference_fields.py\"\n )\n\n\ndef load_scattnlay_case_from_npy(case_key, data_dir=None):\n base = find_data_dir() if data_dir is None else Path(data_dir).resolve()\n paths = {\n \"X\": base / f\"scattnlay_{case_key}_X.npy\",\n \"Z\": base / f\"scattnlay_{case_key}_Z.npy\",\n \"E\": base / f\"scattnlay_{case_key}_E.npy\",\n \"H\": base / f\"scattnlay_{case_key}_H.npy\",\n }\n\n missing = [str(path) for path in paths.values() if not path.exists()]\n if missing:\n msg = \"Missing precomputed scattnlay files:\\n\" + \"\\n\".join(missing)\n msg += \"\\nRun: python docs/generate_scattnlay_reference_fields.py\"\n raise FileNotFoundError(msg)\n\n X = np.load(paths[\"X\"])\n Z = np.load(paths[\"Z\"])\n E = np.load(paths[\"E\"])\n H = np.load(paths[\"H\"])\n return X, Z, E, H\n\n\ndef summarize_errors(label, Eabs_s, Eabs_m, Habs_s, Habs_m):\n E_rel_err = np.abs(Eabs_m - Eabs_s) / np.maximum(Eabs_s, 1e-12)\n H_rel_err = np.abs(Habs_m - Habs_s) / np.maximum(Habs_s, 1e-12)\n\n print(f\"[{label}] scattnlay max |E| = {np.max(Eabs_s):.6g}\")\n print(f\"[{label}] miepython max |E| = {np.max(Eabs_m):.6g}\")\n print(f\"[{label}] scattnlay max |H| (eta-scaled) = {np.max(Habs_s):.6g}\")\n print(f\"[{label}] miepython max |H| = {np.max(Habs_m):.6g}\")\n print(f\"[{label}] E relative error median/max = {np.median(E_rel_err):.3e} / {np.max(E_rel_err):.3e}\")\n print(f\"[{label}] H relative error median/max = {np.median(H_rel_err):.3e} / {np.max(H_rel_err):.3e}\")\n\n return E_rel_err, H_rel_err\n\n\ndef plot_comparison(title, Xg, Zg, radius, Eabs_s, Eabs_m, E_rel_err, Habs_s, Habs_m, H_rel_err):\n fig, axes = plt.subplots(2, 3, figsize=(13, 8), constrained_layout=True)\n\n plots = [\n (Eabs_s, \"scattnlay |E|\", \"turbo\"),\n (Eabs_m, \"miepython |E|\", \"turbo\"),\n (E_rel_err, \"|E| relative error\", \"viridis\"),\n (Habs_s, \"scattnlay |H| (eta-scaled)\", \"magma\"),\n (Habs_m, \"miepython |H|\", \"magma\"),\n (H_rel_err, \"|H| relative error\", \"viridis\"),\n ]\n\n for ax, (arr, ttl, cmap) in zip(axes.flat, plots):\n im = ax.imshow(\n arr.T,\n origin=\"lower\",\n extent=(Xg.min(), Xg.max(), Zg.min(), Zg.max()),\n cmap=cmap,\n aspect=\"equal\",\n )\n circle_color = \"w\" if cmap in (\"magma\", \"viridis\") else \"k\"\n ax.add_patch(plt.Circle((0.0, 0.0), radius, fill=False, color=circle_color, lw=1.1))\n ax.set_title(ttl)\n ax.set_xlabel(\"x\")\n ax.set_ylabel(\"z\")\n fig.colorbar(im, ax=ax, fraction=0.046, pad=0.02)\n\n fig.suptitle(title, y=1.02)\n"
13940 },
14041 {
14142 "cell_type": "markdown",
333234 },
334235 "nbformat": 4,
335236 "nbformat_minor": 5
336- }
237+ }
0 commit comments