{
  "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# Elliptical holes\n\nConvergence checks.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport nannos as nn"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will study the convergence on a benchmark case from\n:cite:p:`Schuster2007`.\nFirst we define the main function that performs the simulation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def array_ellipse(nh, formulation, psi):\n    wl = 500 + 1e-6  # avoid Wood-Rayleigh anomaly\n    d = 1000\n    Nx = 2**9\n    Ny = 2**9\n    lattice = nn.Lattice(([d, 0], [0, d]), discretization=(Nx, Ny))\n    pw = nn.PlaneWave(wavelength=wl, angles=(0, 0, psi))\n    epsgrid = lattice.ones() * (1.75 + 1.5j) ** 2\n    ell = lattice.ellipse((0.5 * d, 0.5 * d), (1000 / 2, 500 / 2), rotate=45)\n    epsgrid[ell] = 1\n\n    sup = lattice.Layer(\"Superstrate\", epsilon=1)\n    sub = lattice.Layer(\"Substrate\", epsilon=1.5**2)\n    st = lattice.Layer(\"Structured\", thickness=50)\n    st.epsilon = epsgrid\n\n    sim = nn.Simulation([sup, st, sub], pw, nh, formulation=formulation)\n    order = (0, 0)\n    R, T = sim.diffraction_efficiencies(orders=True)\n    r = sim.get_order(R, order)\n    return R, T, r, sim\n\n\n#\n# sim = array_ellipse(100, \"original\")\n# lay = sim.get_layer(\"Structured\")\n# lay.plot()\n# plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Perform the simulation for different formulations and number\nof retained harmonics:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "NH = [100, 200, 300, 400, 500, 600]\nformulations = [\"original\", \"tangent\"]\n\n\ndef run_convergence(psi):\n    nhs = {f: [] for f in formulations}\n    rs = {f: [] for f in formulations}\n\n    for nh in NH:\n        print(\"============================\")\n        print(\"number of harmonics = \", nh)\n        print(\"============================\")\n        for formulation in formulations:\n            Ri, Ti, r, sim = array_ellipse(nh, formulation=formulation, psi=psi)\n            R = np.sum(Ri)\n            T = np.sum(Ti)\n            print(\"formulation = \", formulation)\n            print(\"nh0 = \", nh)\n            print(\"nh = \", sim.nh)\n            print(\"r = \", r)\n            print(\"R = \", R)\n            print(\"T = \", T)\n            print(\"R+T = \", R + T)\n            print(\"-----------------\")\n            nhs[formulation].append(sim.nh)\n            rs[formulation].append(r)\n\n    return nhs, rs"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the results:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "markers = {\"original\": \"^\", \"tangent\": \"o\"}\ncolors = {\n    \"original\": \"#d4b533\",\n    \"tangent\": \"#4cb7c6\",\n}\n\nplt.ion()\n\nfor psi in [45, -45]:\n    nhs, rs = run_convergence(psi)\n    plt.figure(figsize=(2, 2))\n\n    for formulation in formulations:\n        plt.plot(\n            nhs[formulation],\n            rs[formulation],\n            \"-\",\n            color=colors[formulation],\n            marker=markers[formulation],\n            label=formulation,\n        )\n        plt.pause(0.1)\n    plt.legend()\n    plt.xlabel(\"number of Fourier harmonics $n_h$\")\n    plt.ylabel(\"$R_{0,0}$\")\n    t = \"\" if psi == 45 else \"-\"\n    plt.title(rf\"$\\psi = {t}45\\degree$\")\n    plt.ylim(0.16, 0.2)\n    plt.tight_layout()\n    plt.show()"
      ]
    },
    {
      "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
}