{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "shape-0",
   "metadata": {
    "tags": [
     "readme"
    ]
   },
   "source": [
    "# Resolving a polygon boundary\n",
    "\n",
    "A concave polygon with slanted edges makes rasterization error visible. All lengths\n",
    "are in nm. `delta` limits the raster-cell size, not the size of the final cuboids:\n",
    "merging adjacent occupied cells into larger cuboids preserves the rasterized geometry\n",
    "exactly. Refining `delta` improves the staircase approximation along oblique edges.\n",
    "The previous axis-aligned L-shape could be represented exactly by two large cuboids;\n",
    "a small cuboid count alone does not imply a coarse approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "shape-1",
   "metadata": {
    "tags": [
     "readme"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from microcubed import cuboidize\n",
    "\n",
    "polygon = np.array([(0, 0), (120, 15), (95, 65), (55, 45), (35, 115), (-15, 80)])\n",
    "thickness = 20\n",
    "magnetization = [0, 0, 8e5]\n",
    "delta = 0.5\n",
    "shape = cuboidize(polygon, t=thickness, delta=delta, mag=magnetization)\n",
    "field = shape.Bfield([50, 50, -60])\n",
    "print(f\"Raster spacing: {delta} nm; merged cuboids: {len(shape)}\")\n",
    "print(\"Field at (50, 50, -60) nm (T):\", field.ravel())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "shape-2",
   "metadata": {},
   "source": [
    "## Compare boundary resolution\n",
    "\n",
    "The orange dashed curve is the input polygon. Blue rectangles are projected\n",
    "cuboids, including their internal boundaries. Compare the staircase error at\n",
    "10 nm, 2 nm, and 0.5 nm; large rectangles in the interior are an exact compression\n",
    "of occupied cells, not a loss of resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "shape-3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.patches import Polygon\n",
    "\n",
    "closed_polygon = np.vstack([polygon, polygon[0]])\n",
    "fig, axes = plt.subplots(1, 3, figsize=(12, 4.5), layout=\"constrained\")\n",
    "for ax, spacing in zip(axes, (10, 2, delta)):\n",
    "    model = shape if spacing == delta else cuboidize(polygon, t=thickness, delta=spacing, mag=magnetization)\n",
    "    for cuboid in model:\n",
    "        ax.add_patch(\n",
    "            Polygon(cuboid.union_boundary(\"xy\")[0].T, facecolor=\"#dbeafe\", edgecolor=\"#2563eb\", linewidth=0.35)\n",
    "        )\n",
    "    ax.plot(*closed_polygon.T, color=\"#c2410c\", linestyle=\"--\", linewidth=1.5, label=\"Input polygon\")\n",
    "    ax.set(\n",
    "        xlim=(-25, 130),\n",
    "        ylim=(-10, 125),\n",
    "        aspect=\"equal\",\n",
    "        xlabel=\"x (nm)\",\n",
    "        ylabel=\"y (nm)\",\n",
    "        title=f\"delta = {spacing:g} nm; {len(model)} cuboids\",\n",
    "    )\n",
    "    ax.legend(fontsize=8, loc=\"lower right\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "shape-4",
   "metadata": {},
   "source": [
    "## Check field convergence\n",
    "\n",
    "Evaluate identical exterior points for every resolution. The reference uses\n",
    "`delta = 0.25 nm`; it is a finer raster approximation, not an exact polygon solution.\n",
    "The relative error below is the RMS vector-field difference divided by the reference\n",
    "RMS vector-field magnitude. Raster errors need not decrease monotonically at every\n",
    "step. Area error and field error measure different properties; repeat this check\n",
    "at the observation distances needed for your application."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "shape-5",
   "metadata": {},
   "outputs": [],
   "source": [
    "probe_x = np.linspace(-25, 130, 81)\n",
    "probes = np.vstack([probe_x, np.full_like(probe_x, 45), np.full_like(probe_x, -60)])\n",
    "reference = cuboidize(polygon, t=thickness, delta=0.25, mag=magnetization).Bfield(probes)\n",
    "reference_norm = np.linalg.norm(reference)\n",
    "polygon_area = (\n",
    "    abs(np.dot(polygon[:, 0], np.roll(polygon[:, 1], -1)) - np.dot(polygon[:, 1], np.roll(polygon[:, 0], -1))) / 2\n",
    ")\n",
    "errors = []\n",
    "print(\"delta (nm) | cuboids | area error (%) | relative field error (%)\")\n",
    "for spacing in (10, 5, 2, 1, delta):\n",
    "    model = shape if spacing == delta else cuboidize(polygon, t=thickness, delta=spacing, mag=magnetization)\n",
    "    values = model.Bfield(probes)\n",
    "    area = sum(np.prod(m.size.ravel()[:2]) for m in model)\n",
    "    error = np.linalg.norm(values - reference) / reference_norm\n",
    "    errors.append(error)\n",
    "    print(\n",
    "        f\"{spacing:10g} | {len(model):7d} | {100 * abs(area - polygon_area) / polygon_area:14.5f} | {100 * error:24.5f}\"\n",
    "    )\n",
    "assert np.isfinite(field).all() and np.isfinite(reference).all()\n",
    "assert np.isfinite(shape.dBfield(probes)).all()\n",
    "assert errors[-1] < errors[0], \"The refined geometry should improve this exterior-field comparison.\"\n",
    "fig, ax = plt.subplots(figsize=(6, 3.5), layout=\"constrained\")\n",
    "ax.loglog((10, 5, 2, 1, delta), np.array(errors) * 100, \"o-\")\n",
    "ax.set(xlabel=\"Raster spacing (nm)\", ylabel=\"Relative field error (%)\", title=\"Compared with the 0.25 nm raster\")\n",
    "ax.grid(which=\"both\", alpha=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "shape-6",
   "metadata": {},
   "source": [
    "## Refined field map and material boundary\n",
    "\n",
    "The field is sampled on a 201 × 201 grid at `z = -60 nm`, below the 20 nm-thick\n",
    "structure. This sampling resolution is independent of the source raster spacing.\n",
    "Thin black lines show the cuboid decomposition; the orange outline marks the input\n",
    "polygon. The cyan `union_boundary(\"xy\")` follows the actual\n",
    "rasterized material, retaining the notch. It removes internal cuboid edges and\n",
    "returns separate closed loops for disconnected pieces and holes. All outlines\n",
    "are XY projections onto the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "shape-7",
   "metadata": {
    "tags": [
     "readme"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = shape.plot_2d(x=(-40, 145, 201), y=(-25, 140, 201), z=-60, component=\"z\")\n",
    "for index, cuboid in enumerate(shape):\n",
    "    ax.plot(\n",
    "        *cuboid.union_boundary(\"xy\")[0],\n",
    "        color=\"black\",\n",
    "        linewidth=0.3,\n",
    "        alpha=0.35,\n",
    "        label=\"Projected cuboids\" if index == 0 else None,\n",
    "    )\n",
    "ax.plot(*np.vstack([polygon, polygon[0]]).T, color=\"#ff9500\", linewidth=1.5, label=\"Input polygon\")\n",
    "for index, boundary in enumerate(shape.union_boundary(\"xy\")):\n",
    "    ax.plot(*boundary, color=\"#00ffff\", linewidth=1.5, label=\"Union boundary\" if index == 0 else None)\n",
    "ax.set(xlabel=\"x (nm)\", ylabel=\"y (nm)\", title=\"Bz (T) at z = -60 nm; delta = 0.5 nm\", aspect=\"equal\")\n",
    "ax.legend(loc=\"upper right\", fontsize=7, facecolor=\"#555555\", labelcolor=\"white\", framealpha=0.95)\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}