{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# This file is part of nannos\n# License: GPLv3\n%matplotlib notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Photonic crystal slab\n\nMetasurface with holes.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nimport nannos as nn"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will study a benchmark of hole in a dielectric surface similar to\nthose studied in :cite:p:`Fan2002`.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define the lattice\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lattice = nn.Lattice([[1.0, 0], [0, 1.0]], discretization=(2**9, 2**9))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define the layers\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sup = lattice.Layer(\"Superstrate\", epsilon=1)\nms = lattice.Layer(\"Metasurface\", thickness=0.5)\nsub = lattice.Layer(\"Substrate\", epsilon=1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define the pattern and add it to the metasurface layer\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ms.epsilon = lattice.ones() * 12.0\ncirc = lattice.circle(center=(0.5, 0.5), radius=0.2)\nms.epsilon[circ] = 1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualize the permittivity\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure()\ncmap = mpl.colors.ListedColormap([\"#ffe7c2\", \"#232a4e\"])\nbounds = [1, 12]\nnorm = mpl.colors.BoundaryNorm(bounds, cmap.N)\nims = ms.plot(cmap=cmap, ax=plt.gca())\nplt.axis(\"scaled\")\nplt.colorbar(ims[0], ticks=bounds)\nplt.xlabel(\"$x$\")\nplt.ylabel(\"$y$\")\nplt.title(r\"permittitivity $\\varepsilon(x,y)$\")\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define the incident plane wave\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pw = nn.PlaneWave(wavelength=1 / 1.4, angles=(0, 0, 0))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define the simulation\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "stack = [sup, ms, sub]\nsim = nn.Simulation(stack, pw, nh=100)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Compute diffraction efficiencies\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "R, T = sim.diffraction_efficiencies()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Compute diffraction efficiencies per order\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Ri, Ti = sim.diffraction_efficiencies(orders=True)\nnmax = 5\nprint(\"Ri = \", Ri[:nmax])\nprint(\"Ti = \", Ti[:nmax])\nprint(\"R = \", R)\nprint(\"T = \", T)\nprint(\"R+T = \", R + T)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, (axR, axT) = plt.subplots(1, 2, figsize=(4, 2))\nlabels = [f\"({g[0]},{g[1]})\" for g in (sim.harmonics[:, :nmax]).T]\naxR.bar(range(nmax), Ri[:nmax], color=[\"#e69049\"])\naxR.set_xticks(range(nmax))\naxR.set_xticklabels(labels=labels)\naxR.set_xlabel(\"order\")\naxR.set_ylabel(\"reflection $R_{i,j}$\")\naxR.annotate(\n    r\"$R = \\sum_i\\,\\sum_j\\, R_{i,j}=$\" + f\"{sum(Ri[:nmax]):0.4f}\",\n    (0.5, 0.9),\n    xycoords=\"axes fraction\",\n)\n\naxT.bar(range(nmax), Ti[:nmax], color=[\"#4a77ba\"])\naxT.set_xticks(range(nmax))\naxT.set_xticklabels(labels=labels)\naxT.set_xlabel(\"order\")\naxT.set_ylabel(\"transmission $T_{i,j}$\")\naxT.annotate(\n    r\"$T =\\sum_i\\,\\sum_j\\, T_{i,j}=$\" + f\"{sum(Ti[:nmax]):0.4f}\",\n    (0.5, 0.9),\n    xycoords=\"axes fraction\",\n)\nplt.suptitle(f\"Diffraction efficiencies: $R+T=${sum(Ri[:nmax] + Ti[:nmax]):0.4f}\")\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Fig 12 (c) from :cite:p:`Fan2002`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def compute_transmission(fn):\n    pw = nn.PlaneWave(wavelength=1 / fn, angles=(0, 0, 0))\n    sim = nn.Simulation(stack, pw, 100)\n    R, T = sim.diffraction_efficiencies()\n    print(f\"f = {fn} (normalized)\")\n    print(\"T = \", T)\n    return T\n\n\nfreqs_norma = np.linspace(0.25, 0.6, 100)\nfreqs_adapted, transmission = nn.adaptive_sampler(\n    compute_transmission,\n    freqs_norma,\n)\n\n\nplt.figure()\nplt.plot(freqs_adapted, transmission, c=\"#be4c83\")\nplt.xlim(freqs_norma[0], freqs_norma[-1])\nplt.ylim(0, 1)\nplt.xlabel(r\"frequency ($2\\pi c / a$)\")\nplt.ylabel(\"Transmission\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nannos.utils.jupyter\n%nannos_version_table"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 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.10.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}