{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# This file is part of gyptis\n# License: MIT\n%matplotlib notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# A scattering simulation in 2D\n\nIn this example we will study the diffraction of a plane wave by a nanorod\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We first need to import the Python packages\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport gyptis as gy\n\nplt.ion()\nplt.close(\"all\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define the wavelength of operation. This can be choosen arbitrarily, but\nit is best for numerical stability not to use too big/small numbers\n(e.g. if we work in optics, it is better to assume the\nunits are in microns amd use ``wavelength = 0.8`` rather than considering\nmeters and ``wavelength = 800e-9``).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "wavelength = 0.8"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Silicon permittivity at this wavelength:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "epsilon_rod = 13.646 - 1j * 0.048"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now define the geometry using the class :class:`~gyptis.BoxPML` a\nsubclass of :class:`~gyptis.Geometry`.\nThis is a rectangular box in 2D (with the argument `dim=2`) centered at the\norigin by default and surrounded by Cartesian Perfectly Matched Layers.\nThe important arguments are its size `box_size` and  the witdh of the\nPMLs along $x$ and $y$.\nAnd optinal argument ``Rcalc`` can be used\nto build a circular boundary containing the scatterer(s) in order to compute\ncross sections.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "geom = gy.BoxPML(\n    dim=2,\n    box_size=(5 * wavelength, 5 * wavelength),\n    pml_width=(wavelength, wavelength),\n    Rcalc=2.4 * wavelength,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we add a circular rod.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "radius = wavelength * 0.5\nrod = geom.add_circle(0, 0, 0, radius)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We use the boolean operation :meth:`~gyptis.BoxPML.fragment` to substract\nthe rod from the box and get the remaining entities:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rod, *box = geom.fragment(rod, geom.box)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Add physical domains:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "geom.add_physical(box, \"box\")\ngeom.add_physical(rod, \"rod\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And set the mesh sizes. A good practice is to have a mesh size that\nis smaller than the wavelength in the media to resolve the field\nso ``size = wavelength / (n*pmesh)``, with ``n`` the refractive index.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pmesh = 10\ngeom.set_pml_mesh_size(wavelength / pmesh)\ngeom.set_size(\"box\", wavelength / pmesh)\ngeom.set_size(\"rod\", wavelength / (pmesh * np.real(epsilon_rod) ** 0.5))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now build and mesh the geometry. The method\n:meth:`~gyptis.BoxPML.build` takes an ``interactive`` boolean argument\nto open and wisualize the geometry in ``gmsh`` (usefull for debugging).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "geom.build()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. attention::\n      A geometry object cannot be modified after the method\n      :meth:`~gyptis.BoxPML.build` has been called: you should create a new object,\n      define the geometry and set mesh parameters before building.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualize the mesh:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(figsize=(2, 2))\ngeom.plot_subdomains(ax=ax)\ngeom.plot_mesh(ax=ax, color=\"red\", lw=0.2)\nplt.axis(\"equal\")\nplt.xlabel(\"$x$ (\u03bcm)\")\nplt.ylabel(\"$y$ (\u03bcm)\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualize the subdomains:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(figsize=(3, 2.3))\nout = geom.plot_subdomains(markers=True, ax=ax)\nplt.axis(\"scaled\")\nplt.xlabel(\"$x$ (\u03bcm)\")\nplt.ylabel(\"$y$ (\u03bcm)\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We define the incident plane wave and plot it. The angle is in radian and\n``theta=0`` corresponds to a wave travelling from the bottom.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pw = gy.PlaneWave(wavelength=wavelength, angle=0, dim=2, domain=geom.mesh, degree=2)\n\nfig, ax = plt.subplots(1, 2, figsize=(4, 1.8))\nfig, ax, maps, cb = pw.plot(ax=ax)\nplt.axis(\"scaled\")\nfor a in ax[:2]:\n    geom.plot_subdomains(ax=a)\n    a.set_xlabel(\"$x$ (\u03bcm)\")\n    a.set_ylabel(\"$y$ (\u03bcm)\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Initialize the simulation. By default, materials properties (``epsilon`` and\n``mu`` arguments are ``None``, and they will be initialized to unity in all\nsubdomains). The values for the PMLs are constructed automatically.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "epsilon = dict(box=1, rod=epsilon_rod)\n\nsimulation = gy.Scattering(\n    geom,\n    epsilon=epsilon,\n    source=pw,\n    degree=2,\n    polarization=\"TE\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We are now ready to solve the problem with the FEM:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "simulation.solve()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The attribute ``simulation.solution`` is a dictionary with keys\n``diffracted`` and ``total`` containing the diffracted and total fields.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(simulation.solution)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lets plot the results. By default, we visualize the real part of the total field:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(1, figsize=(2.5, 2))\nsimulation.plot_field(ax=ax)\ngeom_lines = geom.plot_subdomains()\nplt.xlabel(r\"$x$ (\u03bcm)\")\nplt.ylabel(r\"$y$ (\u03bcm)\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "But we can switch to plot the module of the diffracted field:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(1, figsize=(2.5, 2))\nsimulation.plot_field(ax=ax, field=\"diffracted\", type=\"module\", cmap=\"inferno\")\ngeom_lines = geom.plot_subdomains(color=\"white\")\nplt.xlabel(r\"$x$ (\u03bcm)\")\nplt.ylabel(r\"$y$ (\u03bcm)\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Compute cross sections (in 2D those are rather scattering width *etc.*,\nor cross sections per unit length) and check energy conservation (optical theorem).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cs = simulation.get_cross_sections()\nprint(f'scattering width = {cs[\"scattering\"]:0.4f}\u03bcm')\nprint(f'absorption width = {cs[\"absorption\"]:0.4f}\u03bcm')\nprint(f'extinction width = {cs[\"extinction\"]:0.4f}\u03bcm')\nbalance = (cs[\"scattering\"] + cs[\"absorption\"]) / cs[\"extinction\"]\nprint(f\"    balance      = {balance}\")\nassert np.allclose(cs[\"extinction\"], cs[\"scattering\"] + cs[\"absorption\"], rtol=1e-3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lets recalculate for TM polarization\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "simulationTM = gy.Scattering(\n    geom,\n    epsilon=epsilon,\n    source=pw,\n    degree=2,\n    polarization=\"TM\",\n)\n\nsimulationTM.solve()\nfig, ax = plt.subplots(1, figsize=(2.5, 2))\nsimulationTM.plot_field(ax=ax, field=\"diffracted\", type=\"module\", cmap=\"inferno\")\ngeom_lines = geom.plot_subdomains(color=\"white\")\nplt.xlabel(r\"$x$ (\u03bcm)\")\nplt.ylabel(r\"$y$ (\u03bcm)\")\nplt.tight_layout()\n\n\ncs = simulationTM.get_cross_sections()\nprint(f'scattering width = {cs[\"scattering\"]:0.4f}\u03bcm')\nprint(f'absorption width = {cs[\"absorption\"]:0.4f}\u03bcm')\nprint(f'extinction width = {cs[\"extinction\"]:0.4f}\u03bcm')\nbalance = (cs[\"scattering\"] + cs[\"absorption\"]) / cs[\"extinction\"]\nprint(f\"    balance      = {balance}\")\nassert np.allclose(cs[\"extinction\"], cs[\"scattering\"] + cs[\"absorption\"], rtol=1e-3)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import gyptis.utils.jupyter\n%gyptis_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.11.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}