{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d178512",
   "metadata": {},
   "source": [
    "# Building a PhotonicDevice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6e4c910c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Successfully imported lumapi\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import imodulator\n",
    "import shapely\n",
    "import openbandparams as obp\n",
    "from imodulator.ElectroOpticalModel import InGaAsPElectroOpticalModel\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1daf272",
   "metadata": {},
   "source": [
    "Building a `PhotonicDevice` is where you, as a designer, will be spending the majority of your time. It is a very tedius and cumbersome task, as you have to get all the material parameters right as well as the geometry. Luckily, with the `Imodulator` you only have to do this once. While using this software, we have found it easier to to wrap the creation of the `PhotonicDevice` in a python class, but feel free to do as you see fit. \n",
    "\n",
    "For educational purposes, we will build the class `InP_EOPM` in a very inneficient way, but bear with us. Our goal is to simulate the following:\n",
    "\n",
    "<div style=\"text-align: center;\">\n",
    "  <img src=\"support/parametrization.png\" alt=\"Alt text\" style=\"max-width: 80%; height: auto;\">\n",
    "</div>\n",
    "\n",
    "For the optical refractive indices and dielectric constants of the InGaAsP alloys, we will use `openbandparams`. \n",
    "\n",
    "Before we go to the geometry, we shall define all the mesh settings. For the full workflow of designing a modulator, we must handle 4 meshes: optical mode solver, poisson-drift-diffusion solver, RF solver and electro-optic calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "58b35ae5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def tand_fitted_bcb(x):\n",
    "    \"\"\"\n",
    "    Fitted to results from https://link.springer.com/article/10.1007/s10762-009-9552-0\n",
    "    \n",
    "    x must be in GHz\n",
    "    \"\"\"\n",
    "    out =  0.0093839 - 0.01790336 * np.exp(-0.04773444 * (x - -4.64170761))\n",
    "\n",
    "    if isinstance(x, (list, np.ndarray)):\n",
    "        x = np.asarray(x)\n",
    "    \n",
    "        out[np.where(out<0.001)] = 0.001\n",
    "    else:\n",
    "        if out < 0.001:\n",
    "            out = 0.001\n",
    "    return out\n",
    "\n",
    "class InP_EOPM:\n",
    "    def __init__(\n",
    "            self,\n",
    "            **kwargs\n",
    "    ):\n",
    "        \n",
    "        self.e = 1.60e-19 # electron charge in C\n",
    "        self.e0 = 8.85e-12 # vacuum permittivity in F/m\n",
    "        \n",
    "        self.w_sig_metal = 5 # Width of signal metal in um\n",
    "        self.metal_sep = 10 # Separation between signal and ground metals in um\n",
    "        self.h_metal = 4 # Height of metals in um\n",
    "        self.w_gnd_metal = 10\n",
    "        \n",
    "        self.w_wg = 1\n",
    "        self.h_n = 0.4\n",
    "        self.h_wg1 = 0.5\n",
    "        self.h_wg2 = 0.3\n",
    "        self.h_p1 = 1\n",
    "        self.h_p2 = 0.2\n",
    "\n",
    "        self.h_box = 4\n",
    "\n",
    "        self.w_window = 100\n",
    "        self.h_bottom = 30\n",
    "        self.h_top = 30\n",
    "\n",
    "        for kwarg, value in kwargs.items():\n",
    "            if hasattr(self, kwarg):\n",
    "                setattr(self, kwarg, value)\n",
    "    \n",
    "    def _make_meshes(self):\n",
    "        # optical mesh\n",
    "        self.optical_mesh_settings = {\n",
    "            'substrate': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'background': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'box': {'resolution': 0.3, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'sig_metal': {'resolution': 10, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'n_metal_left': {'resolution': 10, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'n_metal_right': {'resolution': 10, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'bcb': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'n': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'wg1': {'resolution': 0.1, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'wg2': {'resolution': 0.1, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'p1': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'p2': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "        }\n",
    "\n",
    "        # RF mesh\n",
    "        self.rf_mesh_settings = {\n",
    "            'substrate': {'resolution': 5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'background': {'resolution': 5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'box': {'resolution': 3, 'SizeMax': 3, 'distance': 0.1},\n",
    "            'sig_metal': {'resolution': 3, 'SizeMax': 3, 'distance': 0.1},\n",
    "            'n_metal_left': {'resolution': 3, 'SizeMax': 3, 'distance': 0.1},\n",
    "            'n_metal_right': {'resolution': 3, 'SizeMax': 3, 'distance': 0.1},\n",
    "            'bcb': {'resolution': 5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'n': {'resolution': 5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'wg1': {'resolution': 1, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'wg2': {'resolution': 1, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'p1': {'resolution': 5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'p2': {'resolution': 5, 'SizeMax': 5, 'distance': 0.1},\n",
    "        }\n",
    "\n",
    "        # eo mesh\n",
    "        self.eo_mesh_settings = {\n",
    "            'substrate': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'background': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'box': {'resolution': 0.3, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'sig_metal': {'resolution': 10, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'n_metal_left': {'resolution': 10, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'n_metal_right': {'resolution': 10, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'bcb': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'n': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'wg1': {'resolution': 0.1, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'wg2': {'resolution': 0.1, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'p1': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'p2': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "        }\n",
    "\n",
    "        self.charge_mesh_settings = {\n",
    "            'substrate': {'resolution': 0.5},\n",
    "            'background': {'resolution': 0.5},\n",
    "            'sig_metal': {'resolution': 0.01},\n",
    "            'n_metal_left': {'resolution': 0.01},\n",
    "            'n_metal_right': {'resolution': 0.01},\n",
    "            'n': {'resolution': 0.003},\n",
    "            'wg1': {'resolution': 0.002},\n",
    "            'wg2': {'resolution': 0.002},\n",
    "            'p1': {'resolution': 0.003},\n",
    "            'p2': {'resolution': 0.003},\n",
    "        }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "233623f2",
   "metadata": {},
   "source": [
    "Now that we have all the mesh settings, we shall create the geometry which is handled by `shapely`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c22db1f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _create_polygons(self):\n",
    "    #We will now set the RF properties of the metals and the BCB\n",
    "    freq = np.linspace(0.1,100, 100) #GHz. This will be the simulation frequency\n",
    "            \n",
    "    eps_rf_metal = 1 - 1j*6e7/(2*np.pi*freq*1e9 * self.e0)\n",
    "    eps_rf_metal = np.asarray([freq, eps_rf_metal])\n",
    "    \n",
    "    bcb_eps_real = 2.65*np.ones(100)\n",
    "    bcb_eps_imag = bcb_eps_real * tand_fitted_bcb(freq)\n",
    "\n",
    "    bcb_eps = bcb_eps_real - 1j*bcb_eps_imag\n",
    "    bcb_eps = np.asarray([freq, bcb_eps])\n",
    "    \n",
    "    #Now we create the PhotoPolygons\n",
    "    self.substrate = imodulator.SemiconductorPolygon(\n",
    "        shapely.box(\n",
    "            -self.w_window/2,\n",
    "            -self.h_box - self.h_bottom,\n",
    "            self.w_window/2,\n",
    "            -self.h_box\n",
    "        ),\n",
    "        rf_eps = 11.7,\n",
    "        name = 'substrate',\n",
    "        optical_material=3**2,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['substrate'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['substrate'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['substrate'],\n",
    "    )\n",
    "\n",
    "    self.background = imodulator.InsulatorPolygon(\n",
    "        shapely.box(\n",
    "            -self.w_window/2,\n",
    "            -self.h_box - self.h_bottom,\n",
    "            self.w_window/2,\n",
    "            self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2 + self.h_metal + self.h_top\n",
    "        ),\n",
    "        rf_eps = 1,\n",
    "        optical_material=1,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['background'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['background'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['background'],\n",
    "        name = 'background'\n",
    "    )\n",
    "\n",
    "    self.box = imodulator.InsulatorPolygon(\n",
    "        shapely.box(\n",
    "            -self.w_window/2,\n",
    "            -self.h_box,\n",
    "            self.w_window/2,\n",
    "            0\n",
    "        ),\n",
    "        rf_eps = 3.9 - 1j*3.9*0.001,\n",
    "        optical_material=1.44**2,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['box'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['box'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['box'],\n",
    "        name = 'box'\n",
    "    )\n",
    "\n",
    "    n_obp_material = obp.GaInPAs(T=300, As = 0, a = obp.InP.a())\n",
    "    self.n = imodulator.SemiconductorPolygon(\n",
    "        shapely.box(\n",
    "            -self.w_sig_metal/2 - self.metal_sep - self.w_gnd_metal,\n",
    "            0,\n",
    "            self.w_sig_metal/2 + self.metal_sep + self.w_gnd_metal,\n",
    "            self.h_n\n",
    "        ),\n",
    "        rf_eps = n_obp_material.dielectric(T=300),\n",
    "        optical_material=n_obp_material.refractive_index(T=300)**2,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['n'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['n'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['n'],\n",
    "        charge_mesh_settings=self.charge_mesh_settings['n'],\n",
    "        name = 'n',\n",
    "        electro_optic_module=InGaAsPElectroOpticalModel,\n",
    "        electro_optic_module_kwargs={\n",
    "            'y': 0,\n",
    "            'T': 300,\n",
    "            'BF_model': 'vinchant'\n",
    "        },\n",
    "        charge_transport_simulator_kwargs={\n",
    "            'sol_obp_material': n_obp_material,\n",
    "            'sol_Nd': 1e18\n",
    "        }\n",
    "    )\n",
    "\n",
    "    wg1_obp_material = obp.GaInPAs(T=300, As = 0.53, a = obp.InP.a())\n",
    "    self.wg1 = imodulator.SemiconductorPolygon(\n",
    "        shapely.box(\n",
    "            -self.w_wg/2,\n",
    "            self.h_n,\n",
    "            self.w_wg/2,\n",
    "            self.h_n + self.h_wg1\n",
    "        ),\n",
    "        rf_eps = wg1_obp_material.dielectric(T=300),\n",
    "        optical_material=wg1_obp_material.refractive_index(T=300)**2,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['wg1'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['wg1'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['wg1'],\n",
    "        charge_mesh_settings=self.charge_mesh_settings['wg1'],\n",
    "        name = 'wg1',\n",
    "        electro_optic_module=InGaAsPElectroOpticalModel,\n",
    "        electro_optic_module_kwargs={\n",
    "            'y': 0.53,\n",
    "            'T': 300,\n",
    "            'BF_model': 'vinchant'\n",
    "        },\n",
    "        charge_transport_simulator_kwargs={\n",
    "            'sol_obp_material': wg1_obp_material,\n",
    "            'sol_Nd': 1e16\n",
    "        }\n",
    "    )\n",
    "\n",
    "    wg2_obp_material = obp.GaInPAs(T=300, As = 0, a = obp.InP.a())\n",
    "    self.wg2 = imodulator.SemiconductorPolygon(\n",
    "        shapely.box(\n",
    "            -self.w_wg/2,\n",
    "            self.h_n + self.h_wg1,\n",
    "            self.w_wg/2,\n",
    "            self.h_n + self.h_wg1 + self.h_wg2\n",
    "        ),\n",
    "        rf_eps = wg2_obp_material.dielectric(T=300),\n",
    "        optical_material=wg2_obp_material.refractive_index(T=300)**2,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['wg2'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['wg2'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['wg2'],\n",
    "        charge_mesh_settings=self.charge_mesh_settings['wg2'],\n",
    "        name = 'wg2',\n",
    "        electro_optic_module=InGaAsPElectroOpticalModel,\n",
    "        electro_optic_module_kwargs={\n",
    "            'y': 0,\n",
    "            'T': 300,\n",
    "            'BF_model': 'vinchant'\n",
    "        },\n",
    "        charge_transport_simulator_kwargs={\n",
    "            'sol_obp_material': wg2_obp_material,\n",
    "            'sol_Nd': 1e16\n",
    "        }\n",
    "    )\n",
    "\n",
    "    p1_obp_material = obp.GaInPAs(T=300, As = 0, a = obp.InP.a())\n",
    "    self.p1 = imodulator.SemiconductorPolygon(\n",
    "        shapely.box(\n",
    "            -self.w_wg/2,\n",
    "            self.h_n + self.h_wg1 + self.h_wg2,\n",
    "            self.w_wg/2,\n",
    "            self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1\n",
    "        ),\n",
    "        rf_eps = p1_obp_material.dielectric(T=300),\n",
    "        optical_material=p1_obp_material.refractive_index(T=300)**2,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['p1'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['p1'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['p1'],\n",
    "        charge_mesh_settings=self.charge_mesh_settings['p1'],\n",
    "        name = 'p1',\n",
    "        electro_optic_module=InGaAsPElectroOpticalModel,\n",
    "        electro_optic_module_kwargs={\n",
    "            'y': 0,\n",
    "            'T': 300,\n",
    "            'BF_model': 'vinchant'\n",
    "        },\n",
    "        charge_transport_simulator_kwargs={\n",
    "            'sol_obp_material': wg2_obp_material,\n",
    "            'sol_Na': 1e17\n",
    "        }\n",
    "    )\n",
    "\n",
    "    p2_obp_material = obp.GaInAs(T=300, a = obp.InP.a())\n",
    "    self.p2 = imodulator.SemiconductorPolygon(\n",
    "        shapely.box(\n",
    "            -self.w_wg/2,\n",
    "            self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1,\n",
    "            self.w_wg/2,\n",
    "            self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2\n",
    "        ),\n",
    "        rf_eps = p2_obp_material.dielectric(T=300),\n",
    "        optical_material=p2_obp_material.refractive_index(T=300)**2,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['p2'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['p2'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['p2'],\n",
    "        charge_mesh_settings=self.charge_mesh_settings['p2'],\n",
    "        name = 'p2',\n",
    "        electro_optic_module=InGaAsPElectroOpticalModel,\n",
    "        electro_optic_module_kwargs={\n",
    "            'y': 0,\n",
    "            'T': 300,\n",
    "            'BF_model': 'vinchant'\n",
    "        },\n",
    "        charge_transport_simulator_kwargs={\n",
    "            'sol_obp_material': p2_obp_material,\n",
    "            'sol_Na': 1e19\n",
    "        }\n",
    "    )\n",
    "\n",
    "    self.bcb_far_left = imodulator.InsulatorPolygon(\n",
    "        shapely.box(\n",
    "            -self.w_window/2,\n",
    "            0,\n",
    "            -self.w_sig_metal/2 - self.metal_sep - self.w_gnd_metal,\n",
    "            self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2\n",
    "        ),\n",
    "        rf_eps = bcb_eps,\n",
    "        optical_material=1.56**2,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['bcb'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['bcb'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['bcb'],\n",
    "        name = 'bcb_far_left'\n",
    "    )\n",
    "\n",
    "    self.bcb_far_right = imodulator.InsulatorPolygon(\n",
    "        shapely.box(\n",
    "            self.metal_sep + self.w_gnd_metal + self.w_sig_metal/2,\n",
    "            0,\n",
    "            self.w_window/2,\n",
    "            self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2\n",
    "        ),\n",
    "        rf_eps = bcb_eps,\n",
    "        optical_material=1.56**2,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['bcb'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['bcb'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['bcb'],\n",
    "        name = 'bcb_far_right'\n",
    "    )\n",
    "\n",
    "    self.bcb_left = imodulator.InsulatorPolygon(\n",
    "        shapely.box(\n",
    "            -self.w_sig_metal/2 - self.metal_sep,\n",
    "            self.h_n,\n",
    "            -self.w_wg/2,\n",
    "            self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2\n",
    "        ),\n",
    "        rf_eps = bcb_eps,\n",
    "        optical_material=1.56**2,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['bcb'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['bcb'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['bcb'],\n",
    "        name = 'bcb_left'\n",
    "    )\n",
    "\n",
    "    self.bcb_right = imodulator.InsulatorPolygon(\n",
    "        shapely.box(\n",
    "            self.w_wg/2,\n",
    "            self.h_n,\n",
    "            self.w_sig_metal/2 + self.metal_sep,\n",
    "            self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2\n",
    "        ),\n",
    "        rf_eps = bcb_eps,\n",
    "        optical_material=1.56**2,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['bcb'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['bcb'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['bcb'],\n",
    "        name = 'bcb_right'\n",
    "    )\n",
    "\n",
    "    self.sig_metal = imodulator.MetalPolygon(\n",
    "        shapely.box(\n",
    "            -self.w_sig_metal/2,\n",
    "            self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2,\n",
    "            self.w_sig_metal/2,\n",
    "            self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2 + self.h_metal\n",
    "        ),\n",
    "        rf_eps = eps_rf_metal,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['sig_metal'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['sig_metal'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['sig_metal'],\n",
    "        name = 'sig_metal',\n",
    "        calculate_current=True,\n",
    "        d_buffer_current=min(self.w_sig_metal/20, self.h_metal/20, 0.05)\n",
    "    )\n",
    "\n",
    "    self.n_metal_left = imodulator.MetalPolygon(\n",
    "        shapely.box(\n",
    "            -self.w_sig_metal/2 - self.metal_sep - self.w_gnd_metal,\n",
    "            self.h_n,\n",
    "            -self.w_sig_metal/2 - self.metal_sep,\n",
    "            self.h_n + self.h_metal\n",
    "        ),\n",
    "        rf_eps = eps_rf_metal,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['n_metal_left'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['n_metal_left'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['n_metal_left'],\n",
    "        name = 'n_metal_left',\n",
    "        calculate_current=False\n",
    "    )\n",
    "\n",
    "    self.n_metal_right = imodulator.MetalPolygon(\n",
    "        shapely.box(\n",
    "            self.w_sig_metal/2 + self.metal_sep,\n",
    "            self.h_n,\n",
    "            self.w_sig_metal/2 + self.metal_sep + self.w_gnd_metal,\n",
    "            self.h_n + self.h_metal\n",
    "        ),\n",
    "        rf_eps = eps_rf_metal,\n",
    "        eo_mesh_settings=self.eo_mesh_settings['n_metal_right'],\n",
    "        rf_mesh_settings=self.rf_mesh_settings['n_metal_right'],\n",
    "        optical_mesh_settings=self.optical_mesh_settings['n_metal_right'],\n",
    "        name = 'n_metal_right',\n",
    "        calculate_current=False\n",
    "    )\n",
    "\n",
    "\n",
    "InP_EOPM._create_polygons = _create_polygons\n",
    "\n",
    "eopm = InP_EOPM()\n",
    "eopm._make_meshes()\n",
    "eopm._create_polygons()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "765fc3ad",
   "metadata": {},
   "source": [
    "Now that we have created all the polygons, we can initialize our device. A very important aspect to consider is hierarchy. When meshing, the algorithms must now which boundaries to prioritize over others, and of there are overlapping polygons, then we need to know which to keep and which to exclude. The way we do this is through a list. The top elements have priority over bottom elements. We haven't quite figured out what is the proper rule to put polygons, but our reccomendations are:\n",
    "\n",
    "- Always put metal polygons on top. This will ensure that RF calculations have the metalic boundaries well defined.\n",
    "- Put the semiconductor structures in a way that makes sense. That is, in our case, we have a vertical structure, so we will place all the polygons from top to bottom.\n",
    "- Prioritize semiconductor over insulator polygons.\n",
    "- It is mandatory to always include a background polygon. If you don't include it we will create it for you internally. But might as well create it yourself, so that you can control how far away the simulation boundaries are.\n",
    "- Always put the background with least priority.\n",
    "- When doing RF simulations, you can flag `MetalPolygon` with `calculate_current`. By doing so, we will be including in the mesh a line with which we will integrate the magnetic field for current calculations. Any line that is created is put as top priority.\n",
    "- Always make sure you have no empty polygons.\n",
    "\n",
    "It is possible that you get some error due to this hierarchy. Sadly the error messages stemming from `gmsh` will not be very clear, so make sure you try various sensible configurations until you have something that is meshed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5bfaef70",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _initialize_device(self):\n",
    "    photo_polygons = [\n",
    "        self.sig_metal,\n",
    "        self.n_metal_left,\n",
    "        self.n_metal_right,\n",
    "        self.p2,\n",
    "        self.p1,\n",
    "        self.wg2,\n",
    "        self.wg1,\n",
    "        self.n,\n",
    "        self.box,\n",
    "        self.bcb_left,\n",
    "        self.bcb_right,\n",
    "        self.bcb_far_left,\n",
    "        self.bcb_far_right,\n",
    "        self.substrate,\n",
    "        self.background\n",
    "    ]\n",
    "    \n",
    "    #Just in case there are empty polygons\n",
    "    idxs_to_remove = []\n",
    "    for i, poly in enumerate(photo_polygons):\n",
    "        if np.isclose(poly.polygon.bounds[1], poly.polygon.bounds[3]):\n",
    "            idxs_to_remove.append(i)\n",
    "    for i in idxs_to_remove[::-1]:\n",
    "        del photo_polygons[i]\n",
    "    self.device = imodulator.PhotonicDevice(\n",
    "        photo_polygons\n",
    "    )\n",
    "\n",
    "InP_EOPM._initialize_device = _initialize_device\n",
    "eopm = InP_EOPM()\n",
    "eopm._make_meshes()\n",
    "eopm._create_polygons()\n",
    "eopm._initialize_device()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4547dc7d",
   "metadata": {},
   "source": [
    "And we're **DONE**. Now we can inspect our device"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "494c1855",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAJNCAYAAACsgOMnAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAO2VJREFUeJzt3Qmc1XW9P/43IDPDriyCXBZxSTPXxAwtA+WKZZZX46/lAup1Cy2Fm8ItcU3ccs21nwuVBi5p5Vbm7hVMcUlcSEuNK4KawADDLDDn//h+uzMxbALfmTmcc57Px+PrnOV7zvf9Zcb5zOt8lm+bXC6XCwAAgAzaZnkxAABAQrAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsACg4Tz31VBx00EHRt2/faNOmTdx3331Nns/lcjFx4sTYYostokOHDjF8+PB466238lYvQCkQLAAoOEuWLIlddtklrr322tU+f8kll8TVV18dN9xwQzz33HPRqVOnGDFiRFRXV7d6rQClok0u+VgHAApU0mNx7733xsEHH5zeT5q1pCdj3Lhx8V//9V/pYwsXLozevXvHbbfdFocffnieKwYoTpvku4CNTX19fcyZMye6dOmSNlYApSj543zRokXpH+ht2xZW5/Y777wTc+fOTYc/NejWrVvsueeeMW3atDUGi5qamnRbsT345JNPokePHtoDoGTl1qM9ECxWkoSK/v3757sMgI3C7Nmzo1+/flFIklCRSHooVpTcb3hudSZNmhTnnntui9cHUKztgWCxkqSnouEfr2vXrvkuByAvKisr0w9ZGn4nloIJEybE2LFjG+8nw6cGDBgQP39pfHTsUp7X2gDypWpRTRy920Xr1B4IFitp6O5OQoVgAZS6QhwC1KdPn/TrvHnz0lWhGiT3d9111zW+rry8PN1WloSKjl0qWqhagOJpDwpr4CwAfIpBgwal4eLRRx9t0gOTrA41ZMiQvNYGUMz0WABQcBYvXhxvv/12kwnbL7/8cnTv3j0dvnTaaafFBRdcENtuu20aNM4666x04mHDylEAND/BAoCC88ILL8SwYcMa7zfMjRg1alS6pOwZZ5yRXuvihBNOiAULFsSXvvSlePjhh6OiwpAmgJYiWABQcIYOHZougbi2scDnnXdeugHQOsyxAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzDbJ/hYkcrlcVFVV5bsMgCY6duwYbdq0yXcZAJQAwaKZQsWXvvSlePbZZ/NdCkATe31xSDzz7P8IFwC0OEOhmkHSUyFUABujZ6dPi0Xz5ue7DABKgB6LZjbjzKnRuWvnfJcBlLiq2urY7dxD09v11XX5LgeAEiBYNLMkVHTZrFu+ywBKXLuasnyXAECJMRQKAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAAAonWBx/fXXx8477xxdu3ZNtyFDhsRDDz3U+Hx1dXWMGTMmevToEZ07d45DDz005s2bl9eaAQCgVBRMsOjXr19cdNFFMWPGjHjhhRdi3333jW9+85vx2muvpc+ffvrp8bvf/S7uuuuuePLJJ2POnDlxyCGH5LtsAAAoCQUTLA466KD42te+Fttuu2185jOfiR//+Mdpz8T06dNj4cKFcfPNN8fll1+eBo7dd989br311nj22WfT59empqYmKisrm2wAFLbly5fHWWedFYMGDYoOHTrE1ltvHeeff37kcrl8lwZQtAomWKzcYEyZMiWWLFmSDolKejHq6upi+PDhjftsv/32MWDAgJg2bdpa32vSpEnRrVu3xq1///6tcAYAtKSLL744HUL705/+NN544430/iWXXBLXXHNNvksDKFqbRAF59dVX0yCRzKdIeivuvffe2GGHHeLll1+OsrKy2HTTTZvs37t375g7d+5a33PChAkxduzYxvtJj4VwAVDYkh7rZLjsgQcemN7fcsst41e/+lX86U9/yndpAEWroILFdtttl4aIZOjT3XffHaNGjUrnU2RRXl6ebgAUj7322ituuumm+Mtf/pIOn33llVfimWeeSYfMrm1obLI1MDQWoIiDRdIrsc0226S3k3kUzz//fFx11VVx2GGHRW1tbSxYsKBJr0WyKlSfPn3yWDEA+TB+/Pg0GCTDYtu1a5cOoU3m5h1xxBFrHRp77rnntmqdAMWkIOdYNKivr08/XUpCRvv27ePRRx9tfG7WrFnx97//PR06BUBpufPOO+P222+PO+64I1588cWYPHlyXHbZZenXtQ2NTXrEG7bZs2e3as0Aha5geiySX/hf/epX0wnZixYtShuLJ554In7/+9+nk66PO+64dK5E9+7d0+tcnHrqqWmo+OIXv5jv0gFoZT/4wQ/SXovDDz88vb/TTjvFe++9l/ZKJMNoV8fQWIASCRYffvhhHH300fHBBx+kQSK5WF4SKv793/89ff6KK66Itm3bphfGS3oxRowYEdddd12+ywYgD6qqqtI2YUXJkKikpxuAEg8WyXUq1qaioiKuvfbadAOgtCXXPkrmVCS93J/73OfipZdeSiduH3vssfkuDaBoFUywAIB1lVyvIrlA3ne/+920x7tv375x4oknxsSJE/NdGkDREiwAKDpdunSJK6+8Mt0AaB0FvSoUAACwcRAsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADLbJPtbAACs6n8XfBQ1y+ryXQYlbvPOm0aXio75LqMkCBYAQLOqz9XH5Y/dHY/95aV8lwJRvkn7OOuAo+Lz/bfNdylFT7AAAJrV9U//Lh6b9VK0WdYmNu3QKd/lUMLqli+LxXXVccHDv4wLv3FcbN97QL5LKmqCBQDQrMOfHpg5PeKWiNzsXMyPxfkuiRLXaZuKWHJEdUx+7g8x6Rv/me9yiprJ2wBAs1lcszQimVYxO9+VwD8tebs6/ZlMfzZpUXosAIAW8737vxUdO5fnuwxKUN3SZfGTr07NdxklRbAAAFpMEio6de6Q7zIoQbXtrEjW2gyFAgAAMhMsAACA0gkWkyZNij322CO6dOkSm2++eRx88MExa9asJvtUV1fHmDFjokePHtG5c+c49NBDY968eXmrGQAASkXBBIsnn3wyDQ3Tp0+PRx55JOrq6mL//fePJUuWNO5z+umnx+9+97u466670v3nzJkThxxySF7rBgCAUlAwk7cffvjhJvdvu+22tOdixowZsc8++8TChQvj5ptvjjvuuCP23XffdJ9bb701PvvZz6Zh5Itf/GKeKgcAgOJXMD0WK0uCRKJ79+7p1yRgJL0Yw4cPb9xn++23jwEDBsS0adPW+D41NTVRWVnZZAMAAEogWNTX18dpp50We++9d+y4447pY3Pnzo2ysrLYdNNNm+zbu3fv9Lm1zd3o1q1b49a/f/8Wrx8AAIpNQQaLZK7FzJkzY8qUKZnfa8KECWnvR8M2e7ZLhQIUg/fffz+OPPLIdEGPDh06xE477RQvvPBCvssCKFoFM8eiwSmnnBL3339/PPXUU9GvX7/Gx/v06RO1tbWxYMGCJr0WyapQyXNrUl5enm4AFI/58+envdrDhg2Lhx56KHr16hVvvfVWbLbZZvkuDaBoFUywyOVyceqpp8a9994bTzzxRAwaNKjJ87vvvnu0b98+Hn300XSZ2USyHO3f//73GDJkSJ6qBiAfLr744nRoa7KIR4OV243VzblLtgbm3AEU6VCoZPjTL3/5y3TVp+RaFsm8iWRbunRp+nwyP+K4446LsWPHxuOPP55O5j7mmGPSUGFFKIDS8tvf/jYGDx4cI0eOTFcQ3G233eJnP/vZWl9jzh1AiQSL66+/Pp0DMXTo0Nhiiy0at6lTpzbuc8UVV8TXv/71tMciWYI2GQL161//Oq91A9D6/va3v6Xtxrbbbhu///3v4+STT47vfe97MXny5DW+xpw7gBIaCvVpKioq4tprr003AEpXsnpg0mNx4YUXpveTHotk0Y8bbrghRo0atdrXmHMHUCI9FgCwrpIe7R122KHJY8kFU5N5dwC0DMECgKKTrAiVLOCxor/85S8xcODAvNUEUOwECwCKzumnnx7Tp09Ph0K9/fbb6cIfN910U7oQCAAtQ7AAoOjsscce6fLkv/rVr2LHHXeM888/P6688so44ogj8l0aQNEqmMnbALA+klUCkw2A1qHHAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACA0goWTz31VBx00EHRt2/faNOmTdx3331Nns/lcjFx4sTYYostokOHDjF8+PB466238lYvAACUioIKFkuWLIlddtklrr322tU+f8kll8TVV18dN9xwQzz33HPRqVOnGDFiRFRXV7d6rQAAUEo2iQLy1a9+Nd1WJ+mtuPLKK+NHP/pRfPOb30wf+/nPfx69e/dOezYOP/zwVq4WAABKR0H1WKzNO++8E3Pnzk2HPzXo1q1b7LnnnjFt2rQ1vq6mpiYqKyubbAAAQIkGiyRUJJIeihUl9xueW51JkyalAaRh69+/f4vXCgAAxaZogsWGmjBhQixcuLBxmz17dr5LAgCAglM0waJPnz7p13nz5jV5PLnf8NzqlJeXR9euXZtsAABAiQaLQYMGpQHi0UcfbXwsmS+RrA41ZMiQvNYGQH5ddNFF6TLlp512Wr5LAShaBbUq1OLFi+Ptt99uMmH75Zdfju7du8eAAQPSBuOCCy6IbbfdNg0aZ511VnrNi4MPPjivdQOQP88//3zceOONsfPOO+e7FICiVlA9Fi+88ELstttu6ZYYO3Zseju5KF7ijDPOiFNPPTVOOOGE2GOPPdIg8vDDD0dFRUWeKwcgH5J24Igjjoif/exnsdlmm+W7HICiVlA9FkOHDk2vV7EmSTf3eeedl24AMGbMmDjwwAPTpciTHu21SZYfT7YGlh8HKOJgAQDrasqUKfHiiy+mQ6HWRbL8+LnnntvidQEUq4IaCgUA6yJZOvz73/9+3H777es8HNby4wDZ6LEAoOjMmDEjPvzww/j85z/f+Njy5cvjqaeeip/+9KfpkKd27dqtsvx4sgGwYQQLAIrOfvvtF6+++mqTx4455pjYfvvt48wzz1wlVACQnWABQNHp0qVL7Ljjjk0e69SpU/To0WOVxwFoHuZYAAAAmemxAKAkPPHEE/kuAaCo6bEAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyMwF8iAP/lY3O96ufm+DXlsW7WNIp92ivG1Zs9cFALChBAtoZe8tez++fNl3Ijc7t8Hv0Xtgj3jxjPuiTZs2zVobAMCGMhQKWtkrVbMyhYrEvPf+EUsWVTZbTQAAWemxgDx64sAfRscOHdd5/z+UvxI/uv3O9HZu+bIWrAwAYP0IFpBHSajo3KHTOu9f3r68ResBANhQhkIBAACZCRYAAEBmggUAAJCZYAEAAGQmWAAAAJkJFgAAQGaCBQAAkJlgAQAAZCZYAAAAmbnyNhuNXC4XS2urW/24Hcoqok2bNq1+3ELiewMAfBrBgo3mD9dv/uS78cLfXm31Y++x1Y5x37jr/QG7EX5vBm+1U/xm3HW+NwBQAAyFYqOQfBqejz9cE8//bWYsrlqcl2MXgnx+b5LjLqpempdjAwAt3GPxzjvvxNNPPx3vvfdeVFVVRa9evWK33XaLIUOGREVFxfq+Hazi4VFXRaeKDi1+nKV1NbH/Laemt+uXL2/x4xWD+8ZNjY7lLf+9qa6tjm9c9v+lt+vrcy1+PFqONgOgdKxzsLj99tvjqquuihdeeCF69+4dffv2jQ4dOsQnn3wSf/3rX9MG4ogjjogzzzwzBg4c2LJVU9SSUNGpQ8cWP07bTXTYra8kVHSpaPnvzSZtDX0qdNoMgNKzTsEi+XSprKwsRo8eHffcc0/079+/yfM1NTUxbdq0mDJlSgwePDiuu+66GDlyZEvVDMBGTJsBUJrWKVhcdNFFMWLEiDU+X15eHkOHDk23H//4x/Huu+82Z40AFBBtBkBpWqdgsbYGYmU9evRINwBKkzYDoDRt8HKzH374YbrV19c3eXznnXdujroAKCLaDIDit97BYsaMGTFq1Kh444030vXtE8ka88nt5Otyq+sA8H+0GQClY72DxbHHHhuf+cxn4uabb05X+nDhKgDWRJsBUDrWO1j87W9/S1f52GabbVqmIgCKhjYDoHSs90L+++23X7zyyistUw0ARUWbAVA61rvH4v/9v/+XjpedOXNm7LjjjtG+ffsmz3/jG9+IUlZVWx3tasryXUbBqapZmt/j11ZH+1aqoaa2tvH2J9WLoyb+NZl15THoDbcbLF5W3Xh7SU3ys7a0JL43ZS5mWHDftwbaDIDSsd7BIrmo0f/8z//EQw89tMpzpToRr+EPwMRu5x6a11qK7d+ztY6z+8RvRT4ccP+kDX7t5887LErhezP4h99slWMWs9b6vq2ONgOgdKz3x4CnnnpqHHnkkfHBBx+kywauuJVqA1FVVZXvEopKVf2/PpVvSdXL/tVzwLr3HrSG6rqaVjlOqfjHgo/zdmxtBkDpWO8ei3/84x9x+umnp6t7sKp7fnB+9Oi1Wb7LKDjzFy+K/zhvQno716bpOvet4bLD74ge3bq32vFe2+TF+N/l7yQfJTc+trSqKn7ziztW2fcbxxwYPXr+8wJii3KL4tHaJ2J4u33jyM2+k3zk2+K1Lli8IEZdfGx6Ox+fe1/3nXGxeZdN83DkwraganH85y8uTm/X5jGoaTMASsd6B4tDDjkkHn/88dh6661bpqIC17GsIjpXdMp3GQWnti6/n1x2LK+Irh06ttrxhsSXItp/qcljC3Lz4zfxz2Bx4VGj4r9/MTm9fXjfkdGnW5/G/b4bJ0ZrWlaX356dTuUdoksH/0+tr+W51g/oq6PNACgd6x0skvXIJ0yYEM8880zstNNOq0zE+973vtec9UFJ6tn5X5/QW/efQqbNACgdG7QqVOfOnePJJ59MtxUlfwBpJABooM0AKB3rHSzeeeedlqkEgKKjzQAoHRaHh41Q5QorjeVzqVAAgBbrsTj22H+uDrMmt9xyy/q+JbBSgPivyT9rvF1da+lVCpc2A6B0rHewmD9/fpP7dXV16RVVFyxYEPvuu29z1gasdKVuKDTaDIDSsd7B4t57713lseRCRyeffLLlBKEFtG9fke8SYINpMwBKx3oHi9Vp27ZtjB07NoYOHRpnnHFGc7wlsMInvFBMCq3NqF5Sm9bMuqmpqotYoaO1an515GpzJT3ENV/Lhuf7+PlWu3SFa2TVRNTX1Kf/P7N+1uffrFmCReKvf/1rLFu2rLnejhK2oGpRtG9T1uLHmV9V2Xh70dIFUV7eLvJpUfXCxtsTv3F2nPfbc9PbuVx9zF+8IG91LVzyr7oWLFkQ5e1a/g+s5DgNTF4vToXUZhy5y6R8l1DQrv6Pe/JdAkT8JOKdmBuHjDs735UUtfUOFsmnTCtKGv0PPvggHnjggRg1alRz1kYJWfGPx2/fPrHVj3/mnWufYNraGkJF4rtXnxIbi2NvaN2rfidq8nzlb7LRZgCUjvUOFi+99FKT+0n3cK9eveInP/nJp67+0VquvfbauPTSS2Pu3Lmxyy67xDXXXBNf+MIX8l0WsAGq66yKVcgKoc34NKf+8j+jvGPL96IWk9rltXHnrN/G3AUfRslaEhHX/9/t70ZEx1Y+frJq+XX/d/vkiOgUpSf5zPL/RhPvN2ibOHHrIfmuqCBVVdXG6MOntkywePzxx2NjNnXq1PQTshtuuCH23HPPuPLKK2PEiBExa9as2HzzzVv8+AsWL452GqD1tmDJ4sbbE750WGxa1rVVjluzvC5q6yM6tO+yEYxBzUXNsrqoq6+JLpu0j7r6ZVEXy6NTRec81/bPumrraqNTu7JWqWVxbVVc8PQv0tv19cvS4XGsn4VVizeK4WQbe5uxLjp1LouKTuX5LqPAlMdxgw+Ltxa8E3XLS3OeWPXCmngonkhvHz5439hs89ZNFlULqmPydX9Mbx+/457RuUeHKEVvL/ooXl84L47dcUhUWAxlg9SvxxDoZptjsbG4/PLL4/jjj49jjjkmvZ8EjKTLPVkrffz48avsX1NTk24NKiv/Ne5+g4bxXPWvISxsmEnPrFsqpjScfnfDR25sqKrqpfkugRLUvt0msUOPbaNULW5f1Rgs9thiu9ji37q16vEXdloSk+OfwWLvnltG996t84HdxmZYn9L9GcyHdYogBxxwQEyfPv1T91u0aFFcfPHF6VCkfKitrY0ZM2bE8OHDm3S7J/enTZu22tdMmjQpunXr1rj179+/FSsGaHnVC+tb9XiF0mYA0LzWqcdi5MiRceihh6Z/eB900EExePDg6Nu3b1RUVKQXP3r99dfjmWeeiQcffDAOPPDAdH5DPnz88cexfPny6N27d5PHk/tvvvnmal8zYcKEJpMLkx6L9Q0XKw4N+fm+18Wmm/VY79pLXi4XS5fVRHV1TXRts9lGMCyJfEr6AKuXLY3q+qXRrUO+h4IVpvnVC+PYP56U3u7drV+rHrtQ2gxoLYsWVEXHitYdJFI5P5lk8U8W16O1rNNP+XHHHRdHHnlk3HXXXekchptuuikWLvznEpRJg7/DDjuk8xief/75+OxnPxuFpLy8PN2aS/fyXtGrommwYT10yXcBUBwqyv81nru1g1kxtxmwIcOkz/n25LzWUlNdmvNcaH3rHJ+TP76ThiLZEkkjsXTp0ujRo0e0b98+NgY9e/aMdu3axbx585o8ntzv06dP3uoCKDWF0GZAqagWLGglG9wv1zAnYWNSVlYWu+++ezz66KNx8MEHp4/V19en9085ZeO5FgBAqdkY2wxoSSv2FH770iGxac/WXZGoamFtTD7tf9Lbnbta1YzWUXSrQiXzJZKLLiVjepNrVyTLzS5ZsqRxlSgAgNbUtUfH6NGnc6ses7yiuvG2aWq0lnVfmLZAHHbYYXHZZZfFxIkTY9ddd42XX345Hn744VUmdANQvJIV//bYY4/o0qVLeg2jpBc7uZ4RAC2n6IJFIhn29N5776XXp3juuefSC+UBUDqefPLJGDNmTLrs7SOPPBJ1dXWx//77pz3YALSMohsKBQBJT/WKbrvttrTnIrnW0T777LPa1zTHBVMBStl691gk8xeeeuqplqkGgKKysbQZDcvddu/efY37uGAqQCsHi+SXc3Il62233TYuvPDCeP/99zOWAECx2hjajGR1wNNOOy323nvv2HHHHde4X3LB1KTehm327NmtWidAyQWL++67L20YTj755PTCR1tuuWV89atfjbvvvjsdwwoAG1Obkcy1mDlzZkyZMuVTr73RtWvXJhsALTx5u1evXumyrq+88ko6OXqbbbaJo446Kvr27Runn356vPXWWxvytgAUoXy2GcliHvfff388/vjj0a9fvxY7DgAZV4X64IMP0tU2ki254vXXvva1ePXVV2OHHXaIK664ovmqBKDgtWabkcvl0lBx7733xmOPPRaDBg1q1vcHoBmCRdJ1fc8998TXv/71GDhwYNx1113p2NU5c+bE5MmT449//GPceeedcd55563vWwNQZPLVZiTDn375y1/GHXfckV7LYu7cuem2dOnSZj0OABmWm91iiy3SiXDf/va3409/+lN6EbqVDRs2LDbddNP1fWsAiky+2ozrr78+/Tp06NAmj996660xevToZj0WABsYLJLu6pEjR0ZFRcUa90kaiHfeeWd93xqAIpOvNiMZCgXARh4skgl3ALAutBkApSPT5G0AAICEYAEAAGQmWAAAAJkJFgAAQGaCBQAAkJlgAQAAZCZYAAAArX8dC9ZuQc2CaFddlu8ygBK3oHpB4+1cuFgcAC1PsGjmK7we+djxea0FYGVVVVX5LgGAEmAoFECRq6tblO8SACgBeiyaQZs2bRpvP/itSdFzsx55rQfgH0sq46t3/Fd6u3efbvkuB4ASIFg0s16bdovePbrnuwygxG1S0Xa1H34AQEsxFAoAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyMyqUAAALaiqsi4Wd6hu1WMuXVjbeHuF6/hCixIsAACaWW6Fv+Z/cdpTea2ltrour8endBgKBQBQxGpr6vNdAiVCjwUAQDNb8cKU37vs29G9Z+dWPf7ihUvisu/fkd7u2q1Dqx6b0iVYAAC0oM16dY1eW3Rr1WOWdyhbbciBlmQoFAAAkJlgAQAAZCZYAAAAmQkWAABAZoIFAACQmWABAABkZrlZAIAWtLiyKioqWvdPriULqxpvr3ARcGhRggUAQDPLrfDX/E9O/UVea6mprsvr8SkdhkIBABSx2prl+S6BEqHHAgCgma14teszzh0VPXt0bdXjL1pYFef/8Ob0dteuFa16bEqXYAEA0IK69ewcPXu3brAoq2i72pADLclQKAAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMw2yf4WAFC8lixYGstql+e7DArMkoVLG28vrqyKyvLW/ZNr0aJ/HX/hgupo26Fdqx6f4rG0qq74gsWPf/zjeOCBB+Lll1+OsrKyWLBgwSr7/P3vf4+TTz45Hn/88ejcuXOMGjUqJk2aFJtsUjCnCcBG5prRt+a7BArcpAn5/Rn6/nfvy+vxKR0F8xd3bW1tjBw5MoYMGRI333zzKs8vX748DjzwwOjTp088++yz8cEHH8TRRx8d7du3jwsvvDAvNQMAQKkomGBx7rnnpl9vu+221T7/hz/8IV5//fX44x//GL17945dd901zj///DjzzDPjnHPOSXs5VqempibdGlRWVrbQGQBQiEbd840o69Q+32VQaHK5qK1eFstrlkeHzcqiTZs2rXz4XNQtXRZRsyw6dytv9eNTPGqX1MWNh9xfXMHi00ybNi122mmnNFQ0GDFiRDo06rXXXovddtttta9Lhko1hBYAWFmPf+sU5Z1X/+EUQLGrWVxbeqtCzZ07t0moSDTcT55bkwkTJsTChQsbt9mzZ7d4rQAAUGzyGizGjx+fds2tbXvzzTdbtIby8vLo2rVrkw0AAFg/eR0KNW7cuBg9evRa99lqq63W6b2SSdt/+tOfmjw2b968xucAAIAiDRa9evVKt+aQrBaVLEn74Ycfxuabb54+9sgjj6Q9EDvssEOzHAMAACjwydvJNSo++eST9GuytGxyPYvENttsk16zYv/9908DxFFHHRWXXHJJOq/iRz/6UYwZMyYd7gQAALScggkWEydOjMmTJzfeb1jlKbkY3tChQ6Ndu3Zx//33p6tAJb0XnTp1Si+Qd9555+WxagAAKA0FEyyS61es6RoWDQYOHBgPPvhgq9UEAAAU2XKzALCya6+9NrbccsuoqKiIPffcc5VFPgBoPoIFAEVp6tSpMXbs2Dj77LPjxRdfjF122SW9cGqyyAcAJTwUCgDWx+WXXx7HH398HHPMMen9G264IR544IG45ZZb0usorasln1THsprlLVgpwMarZkndOu8rWABQdGpra2PGjBkxYcKExsfatm0bw4cPj2nTpq32NTU1NenWoLKyMv1649fubYWKAQqfoVAAFJ2PP/44XZq8d+/eTR5P7ifLka/OpEmTolu3bo1b//79W6lagOKgxwIAItLejWROxoo9Fkm4OOfpo6Oic1leawPIl+rFtXHOl3++TvsKFgAUnZ49e6bXN5o3b16Tx5P7ffr0We1rkoupru6Cqr027xYdurjQKlCalnb81xDRT2MoFABFp6ysLHbfffd49NFHGx+rr69P7ycXUQWg+emxAKAoJcOaRo0aFYMHD44vfOELceWVV8aSJUsaV4kCoHkJFgAUpcMOOyw++uijmDhxYjphe9ddd42HH354lQndADQPwQKAonXKKaekGwAtzxwLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAoDSCxbvvvhvHHXdcDBo0KDp06BBbb711nH322VFbW9tkvz//+c/x5S9/OSoqKqJ///5xySWX5K1mAAAoJZtEAXjzzTejvr4+brzxxthmm21i5syZcfzxx8eSJUvisssuS/eprKyM/fffP4YPHx433HBDvPrqq3HsscfGpptuGieccEK+TwEAAIpaQQSLAw44IN0abLXVVjFr1qy4/vrrG4PF7bffnvZg3HLLLVFWVhaf+9zn4uWXX47LL79csAAAgBZWEEOhVmfhwoXRvXv3xvvTpk2LffbZJw0VDUaMGJEGkPnz56/xfWpqatLejhU3AACgBILF22+/Hddcc02ceOKJjY/NnTs3evfu3WS/hvvJc2syadKk6NatW+OWzM0AAAAKKFiMHz8+2rRps9YtmV+xovfffz8dFjVy5Mh0nkVWEyZMSHs/GrbZs2dnfk8AACg1eZ1jMW7cuBg9evRa90nmUzSYM2dODBs2LPbaa6+46aabmuzXp0+fmDdvXpPHGu4nz61JeXl5ugEAAAUaLHr16pVu6yLpqUhCxe677x633nprtG3btLNlyJAh8cMf/jDq6uqiffv26WOPPPJIbLfddrHZZpu1SP0AAEABzbFIQsXQoUNjwIAB6SpQH330UTpvYsW5E9/5znfSidvJ9S5ee+21mDp1alx11VUxduzYvNYOAACloCCWm016HpIJ28nWr1+/Js/lcrn0azLx+g9/+EOMGTMm7dXo2bNnTJw40VKzAADQCgoiWCTzMD5tLkZi5513jqeffrpVagIAAApsKBQAALBxEywAAIDMBAsAACAzwQKAovLuu++mKwQOGjQoOnToEFtvvXWcffbZUVtbm+/SAIpaQUzeBoB19eabb0Z9fX3ceOONsc0228TMmTPj+OOPjyVLlqRLlgPQMgQLAIrKAQcckG4Nttpqq5g1a1Zcf/31aw0WNTU16dagsrKyxWsFKCaGQgFQ9BYuXBjdu3df6z6TJk1Kr4nUsPXv37/V6gMoBoIFAEUtubjqNddcEyeeeOJa95swYUIaQBq22bNnt1qNAMVAsACgIIwfPz7atGmz1i2ZX7Gi999/Px0WNXLkyHSexdqUl5dH165dm2wArDtzLAAoCOPGjYvRo0evdZ9kPkWDOXPmxLBhw2KvvfaKm266qRUqBChtggUABaFXr17pti6SnookVOy+++5x6623Rtu2OugBWppgAUBRSULF0KFDY+DAgekqUB999FHjc3369MlrbQDFTLAAoKg88sgj6YTtZOvXr1+T53K5XN7qAih2+oYBKCrJPIwkQKxuA6DlCBYAAEBmggUAAJCZYAEAAGQmWAAAAJkJFgAAQGaCBQAAkJlgAQAAZCZYAAAAmQkWAABAZoIFAACQmWABAABkJlgAAACZCRYAAEBmggUAAJCZYAEAAGQmWAAAAJkJFgAAQGaCBQAAkJlgAQAAZCZYAAAAmQkWAABAZoIFAACQmWABAABkJlgAAACZCRYAAEBmggUAAJCZYAEAAGQmWAAAAJkJFgAAQGaCBQAAkJlgAQAAZCZYAAAAmQkWAABAZoIFAACQmWABAABkJlgAAACZCRYAAEBmggUAAJCZYAEAAGQmWAAAAJkJFgAAQGaCBQAAkJlgAQAAlE6w+MY3vhEDBgyIioqK2GKLLeKoo46KOXPmNNnnz3/+c3z5y19O9+nfv39ccskleasXAABKScEEi2HDhsWdd94Zs2bNinvuuSf++te/xre+9a3G5ysrK2P//fePgQMHxowZM+LSSy+Nc845J2666aa81g0AAKVgkygQp59+euPtJDyMHz8+Dj744Kirq4v27dvH7bffHrW1tXHLLbdEWVlZfO5zn4uXX345Lr/88jjhhBPW+L41NTXptmJAAQAAirTHYkWffPJJGiT22muvNFQkpk2bFvvss08aKhqMGDEi7eGYP3/+Gt9r0qRJ0a1bt8YtGUIFAAAUcbA488wzo1OnTtGjR4/4+9//Hr/5zW8an5s7d2707t27yf4N95Pn1mTChAmxcOHCxm327NkteAYAAFCc8hoskuFMbdq0Wev25ptvNu7/gx/8IF566aX4wx/+EO3atYujjz46crlcphrKy8uja9euTTYAAKCA5liMGzcuRo8evdZ9ttpqq8bbPXv2TLfPfOYz8dnPfjYdtjR9+vQYMmRI9OnTJ+bNm9fktQ33k+cAAIAiDRa9evVKtw1RX1+ffm2YeJ2Eix/+8IeNk7kTjzzySGy33Xax2WabNWPVAABAQc6xeO655+KnP/1pusrTe++9F4899lh8+9vfjq233joNFInvfOc76cTt4447Ll577bWYOnVqXHXVVTF27Nh8lw8AAEWvIIJFx44d49e//nXst99+aQ9EEh523nnnePLJJ9M5EolkRadk7sU777wTu+++ezrMauLEiWtdahYAACih61jstNNOaS/Fp0nCxtNPP90qNQEAAAXWYwEAAGzcBAsAACAzwQIAAMhMsAAAADITLAAoWsm1jnbddddo06ZNumQ5AC1HsACgaJ1xxhnRt2/ffJcBUBIKYrlZAFhfDz30UHp9o3vuuSe9vaEW/aMq6mqWNWttAIWienHtOu8rWABQdObNmxfHH3983HfffelFVtd12FSyNaisrEy/jt/r5harE6CYGAoFQFHJ5XIxevToOOmkk2Lw4MHr/LpJkyZFt27dGrf+/fu3aJ0AxUaPBQAFYfz48XHxxRevdZ833ngjHf60aNGimDBhwnq9f7L/2LFjm/RYJOHipF9/M8o7td/gugEKWc2SurjhkN+s076CBQAFYdy4cWlPxNpstdVW8dhjj8W0adOivLy8yXNJ78URRxwRkydPXu1rk/1Xfk2ie68uUdG5LGP1AIXJHAsAik6vXr3S7dNcffXVccEFFzTenzNnTowYMSKmTp0ae+65ZwtXCVC6BAsAisqAAQOa3O/cuXP6deutt45+/frlqSqA4mfyNgAAkJkeCwCK2pZbbpmuFAVAy9JjAQAAZCZYAAAAmRkK1cw+Wboo2i2xLCGQX59ULWq8XZ+rz2stAJQGwaKZjbj9zHyXANDEJ/Pnx8B8FwFA0TMUqhn07NkzunTuku8yAFarXcdVL/oGAM1Nj0UzaNu2bSxYuCDmzno3li+pznc5AOnwp6SnIgkVO+41ON/lAFACBItmDBd9P7tVvssAaGT4EwCtyVAoAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDMBAsAACAzwQIAAMhMsAAAADITLAAAgMw2yf4WxSWXy6VfKysr810KQN40/A5s+J1YihrOvWZJbb5LAcibht+B69IeCBYrWbRoUfq1f//++S4FYKP4nditW7co5fbgJwfcme9SAAqiPWiTK+WPo1ajvr4+5syZE126dIk2bdpEIXyqmISg2bNnR9euXaNYlcp5JpxrcSq0c02ahqQR6du3b7RtW5qjZvPdHhTaz0xLKPV/g1I//4R/g8j7v8H6tAd6LFaS/IP169cvCk3yg1YK/8OVynkmnGtxKqRzLdWeio2tPSikn5mWUur/BqV+/gn/BpHXf4N1bQ9K82MoAACgWQkWAABAZoJFgSsvL4+zzz47/VrMSuU8E861OJXSudI8/Mz4Nyj180/4N4iC+jcweRsAAMhMjwUAAJCZYAEAAGQmWAAAAJkJFgAAQGaCRRGoqamJXXfdNb0y7Msvv9zkuT//+c/x5S9/OSoqKtKrNl5yySVRSN5999047rjjYtCgQdGhQ4fYeuut05URamtri+o8V3TttdfGlltumZ7LnnvuGX/605+ikE2aNCn22GOP9OrFm2++eRx88MExa9asJvtUV1fHmDFjokePHtG5c+c49NBDY968eVHoLrroovT/y9NOO63oz5X8/74vZuvaFhSbYmsPmrvtKCUXraY92RgJFkXgjDPOSC+zvrpLwO+///4xcODAmDFjRlx66aVxzjnnxE033RSF4s0334z6+vq48cYb47XXXosrrrgibrjhhvjv//7vojrPBlOnTo2xY8emDeaLL74Yu+yyS4wYMSI+/PDDKFRPPvlk+of09OnT45FHHom6urr0+7VkyZLGfU4//fT43e9+F3fddVe6/5w5c+KQQw6JQvb888+nP7c777xzk8eL8VzJ/+/7YrcubUGxKcb2oLnbjlLx/Brak41SstwshevBBx/Mbb/99rnXXnstWTY499JLLzU+d9111+U222yzXE1NTeNjZ555Zm677bbLFbJLLrkkN2jQoKI8zy984Qu5MWPGNN5fvnx5rm/fvrlJkyblta7m9OGHH6Y/q08++WR6f8GCBbn27dvn7rrrrsZ93njjjXSfadOm5QrRokWLcttuu23ukUceyX3lK1/Jff/73y/ac2Xj+H1filZuC4pNKbQHWdqOUrFoDe3JxkqPRQFLhk8cf/zx8Ytf/CI6duy4yvPTpk2LffbZJ8rKyhofSz7tSLoS58+fH4Vq4cKF0b1796I7z6RLP+lxGT58eONjbdu2Te8n51gsku9fouF7mJxz8knUiue9/fbbx4ABAwr2vJNP2Q488MAm51Ss58rG8fu+FK3cFhSTUmkPsrQdpWLMGtqTjZVgUaCS6xqOHj06TjrppBg8ePBq95k7d2707t27yWMN95PnCtHbb78d11xzTZx44olFd54ff/xxLF++fLXnUkjnsTbJUIZkfOjee+8dO+64Y/pYcm5JKNx0002L4rynTJmSDltIxgevrNjOlY3n932pWV1bUExKoT3I2naUgilraU82VoLFRmb8+PHp5Jy1bclY0+QX6qJFi2LChAlRzOe5ovfffz8OOOCAGDlyZPrJHYX5ycvMmTPTX5bFaPbs2fH9738/br/99nSyJaxNqfy+XxttAeui2NuOYmpPNsl3ATQ1bty49JOptdlqq63iscceS7tDy8vLmzyXfJp1xBFHxOTJk6NPnz6rrDbTcD95rhDOs0EywXXYsGGx1157rTIpe2M+z/XRs2fPaNeu3WrPpZDOY01OOeWUuP/+++Opp56Kfv36NT6enFvS7b9gwYImn+QX4nknQxeSiZWf//znGx9LPnVMzvmnP/1p/P73vy+ac2Xj+n1fqJqzLSgmxd4eNEfbUexmfEp7kqwQl/yMbHTyPcmDDfPee+/lXn311cbt97//fTqp6e67787Nnj27yaTm2traxtdNmDCh4CY1/+///m86cenwww/PLVu2bJXni+U8GybrnXLKKU0m6/3bv/1bQU/Wq6+vTycgJpMO//KXv6zyfMOE5uRnt8Gbb75ZkBOaKysrm/x/mWyDBw/OHXnkkentYjpXNq7f96Xg09qCYlOM7UFzth3FrvJT2pONlWBRJN55551VVglJ/ojp3bt37qijjsrNnDkzN2XKlFzHjh1zN954Y66QGpJtttkmt99++6W3P/jgg8atmM6zQVJ7eXl57rbbbsu9/vrruRNOOCG36aab5ubOnZsrVCeffHKuW7duuSeeeKLJ96+qqqpxn5NOOik3YMCA3GOPPZZ74YUXckOGDEm3YrDyKh7FfK7k7/d9sVuXtqDYFGN70NxtR6n5SgGsCiVYFHlD88orr+S+9KUvpb+ckk86LrroolwhufXWW9PzWt1WTOe5omuuuSb9w7OsrCz9xGr69Om5Qram71/yvW2wdOnS3He/+9205ykJhf/xH/9RNH8wrNwQFPO50jpKMVisa1tQbIqtPWjutqPUfKUAgkWb5D/5Ho4FAAAUNqtCAQAAmQkWAABAZoIFAACQmWABAABkJlgAAACZCRYAAEBmggUAAJCZYAEAAGQmWEAe3HzzzbH//vu3+HEefvjh2HXXXaO+vr7FjwXA+tMeUEwEC2hl1dXVcdZZZ8XZZ5/d4sc64IADon379nH77be3+LEAWD/aA4qNYAGt7O67746uXbvG3nvv3SrHGz16dFx99dWtciwA1p32gGIjWMAG+uijj6JPnz5x4YUXNj727LPPRllZWTz66KNrfN2UKVPioIMOavLY0KFD47TTTmvy2MEHH5w2Ag223HLLuOCCC+Loo4+Ozp07x8CBA+O3v/1tWsc3v/nN9LGdd945XnjhhSbvkxwreeyvf/1rM5w1ACvTHsA/CRawgXr16hW33HJLnHPOOekv6kWLFsVRRx0Vp5xySuy3335rfN0zzzwTgwcP3qBjXnHFFeknWy+99FIceOCB6fGShuXII4+MF198Mbbeeuv0fi6Xa3zNgAEDonfv3vH0009v0DEBWDvtAfyTYAEZfO1rX4vjjz8+jjjiiDjppJOiU6dOMWnSpDXuv2DBgli4cGH07dt3g4934oknxrbbbhsTJ06MysrK2GOPPWLkyJHxmc98Js4888x44403Yt68eU1elxzvvffe26BjAvDptAcgWEBml112WSxbtizuuuuudFJceXn5GvddunRp+rWiomKDjpV0bTdIPnVK7LTTTqs89uGHHzZ5XYcOHaKqqmqDjgnAutEeUOoEC8goGas6Z86cdAm/d999d6379ujRI9q0aRPz589v8njbtm2bdFcn6urqVnl9sqJHg+R91vTYyssJfvLJJ2lXPQAtR3tAqRMsIIPa2tp0POthhx0W559/fvznf/7nKp8OrSiZyLfDDjvE66+/3uTx5Jf8Bx980Hh/+fLlMXPmzGZbzjBp7HbbbbdmeT8AVqU9AMECMvnhD3+YjpFNlu9LxrMm41qPPfbYtb5mxIgR6YS9Fe27777xwAMPpNubb74ZJ598cjr+tjlMnz497Y4fMmRIs7wfAKvSHoBgARvsiSeeiCuvvDJ+8YtfpOuQJ93Xye1ktY3rr79+ja877rjj4sEHH0wboAZJ4zNq1Kh0BY+vfOUrsdVWW8WwYcOapc5f/epX6WTCjh07Nsv7AdCU9gD+qU1u5YF8QItLVu34/Oc/HxMmTGjR43z88cex3XbbpcsfDho0qEWPBcD60x5QTPRYQB5ceuml6QWMWloyefC6667TiABspLQHFBM9FgAAQGZ6LAAAgMwECwAAIDPBAgAAyEywAAAAMhMsAACAzAQLAAAgM8ECAADITLAAAAAyEywAAIDI6v8HqpCFUYvsGPoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "gs = fig.add_gridspec(1,2)\n",
    "\n",
    "ax1 = fig.add_subplot(gs[0,0])\n",
    "ax2 = fig.add_subplot(gs[0,1])\n",
    "\n",
    "for ax in [ax1, ax2]:\n",
    "    eopm.device.plot_polygons(\n",
    "        color_polygon=\"black\",\n",
    "        color_line=\"green\",\n",
    "        color_junctions=\"blue\",\n",
    "        fill_polygons=True,\n",
    "        fig=fig,\n",
    "        ax=ax,\n",
    "    )\n",
    "\n",
    "ax2.set_xlim(-5,5)\n",
    "ax2.set_ylim(-5,10)\n",
    "\n",
    "ax1.set_xlabel('x (um)')\n",
    "ax1.set_ylabel('y (um)')\n",
    "ax2.set_xlabel('x (um)')\n",
    "ax2.set_ylabel('y (um)')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88230354",
   "metadata": {},
   "source": [
    "For the moment we cannot check how the physical properties look like. But we will be able to do once we start working with the different simulators.\n",
    "\n",
    "Because we deal with shapely, it is possible to have arbitrary shapes with each polygon as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd8ce97f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAJNCAYAAACsgOMnAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWh5JREFUeJzt3Qd4VFX+xvE3fZJAQhVEijQRpQqIgCIIUsS2KmtBBUVEBF3AVeC/u7rFFVddRRFBXEVXcUHsIqgIiiKhF+kqAiJVBJKQZGZS5v+cg4mEZsJNcjMz38/zXHOn5f4GJCe/OfeeNyIQCAQEAAAAAA5EOnkxAAAAABg0FgAAAAAco7EAAAAA4BiNBQAAAADHaCwAAAAAOEZjAQAAAMAxGgsAAAAAjtFYAAAAAHCMxgIAAACAYzQWAAAAAByjsQAABJ0vvvhCV1xxhWrVqqWIiAi9++67hR4PBAJ68MEHdfrppys+Pl7du3fXt99+61q9ABAOaCwAAEEnIyNDLVu21IQJE477+GOPPaZnnnlGkyZN0uLFi5WYmKiePXvK6/WWea0AEC4iAuZjHQAAgpSZsXjnnXd09dVX29tmWDMzGffdd5/++Mc/2vtSU1NVo0YNvfzyy7rhhhtcrhgAQlO02wWUN3l5edq5c6cqVqxoBysACEfml/P09HT7C3pkZHBNbm/ZskW7d++2pz/lS05OVvv27ZWSknLCxsLn89ntyPFg//79qlq1KuMBgLAVKMZ4QGNxFNNU1KlTx+0yAKBc2L59u2rXrq1gYpoKw8xQHMnczn/seMaOHau//e1vpV4fAITqeEBjcRQzU5H/h5eUlOR2OQDgirS0NPshS/7PxHAwZswYjRw5suC2OX2qbt26+s9XE5VQId7V2gDALZmHsnRHpyFFGg9oLI6SP91tmgoaCwDhLhhPAapZs6b9umfPHrsqVD5zu1WrVid8XVxcnN2OZpqKhIoJpVQtAITOeBBcJ84CAPAb6tevb5uLuXPnFpqBMatDdejQwdXaACCUMWMBAAg6hw4d0nfffVfogu1Vq1apSpUq9vSl4cOH6+GHH1bjxo1to/GXv/zFXniYv3IUAKDk0VgAAILOsmXL1LVr14Lb+ddG9O/f3y4p+8ADD9isizvvvFMHDx7UhRdeqI8++kgej8fFqgEgtNFYAACCTpcuXewSiCc7F/jvf/+73QAAZYNrLAAAAAA4RmMBAAAAwDEaCwAAAACO0VgAAAAAcIzGAgAAAIBjNBYAAAAAHKOxAAAAAOAYjQUAAAAAx2gsAAAAADhGYwEAAADAMRoLAAAAAI7RWAAAAABwjMYCAAAAgGM0FgAAAAAci3b+LWAEAgFlZma6XQYAFJKQkKCIiAi3ywAAhAEaixJqKi688EItXLjQ7VIAoJAOHTvoqwVf0VwAAEodp0KVADNTQVMBoDxKWZiiA+k/u10GACAMMGNRwqateUEVEiu6XQaAMOfN9Oq6ZgPsfnbA73Y5AIAwQGNRwkxTkVwx2e0yAIS52Mg4t0sAAIQZToUCAAAA4BiNBQAAAADHaCwAAAAAOEZjAQAAAMAxGgsAAAAAjtFYAAAAAHCMxgIAAACAYzQWAAAAAByjsQAAAADgGI0FAAAAAMdoLAAAAAA4RmMBAAAAwDEaCwAAAACO0VgAAAAAcIzGAgAAAIBjNBYAAAAAHKOxAAAAAOAYjQUAAAAAx2gsAAAAADhGYwEAAADAMRoLAAAAAI7RWAAAAABwjMYCAAAAgGM0FgAAAAAco7EAAAAA4BiNBQAAAADHaCwAAAAAOBbt/FsAAACgpHizvXro04f17b5v7e2oiGj1b9NPlze9zO3SgJNixgIAAKAceWXFa9r00zfK+zygvPEBZe/K1pRlr+r7n7e4XRpwUjQWAAAA5cTyHSs1e9Mn0mZJn0naJyW871GuP1dPLnhGvhy/2yUCJ0RjAQAAUA6kedP07MKJUqYU90Fswf2ZO7yK+zJWP6bu0KsrprpaI3AyNBYAAAAuCwQCem7RZB3IPKj4TzzyHfTr9FpJGj6mp33c96Vf2iLN3Dhbq3Z+7Xa5wHHRWAAAALhs3ub5WvTDEkWujVDWKq8iIyN0z6he6tS1ibr1aioFJI+ZxciSnvlqgtJ9h9wuGTgGjQUAAICL9qTv1X+WTpEOSlGzo+x9193UVo2b1rT7/YdcrBo1K8q736/4Tz3an3VAkxa9YGc5gPKExgIAAMAluXl5GvfVs8ryZSnhQ4+yM3PUuEl1XdOvfcFz4uNj7exFRGSEspZ7FbEuQl9tS9EXWxa4WjtwNBoLAAAAl7y77n1t2LtR0UuilfmtV3FxUbpnVG9FRRX+Fa3Juafrmuvb2P2YWVFSmvT84hf106F9LlUOHIvGAgAAwAUml+J/q6dLe6S8T/Psff0HX6TTa1c67vOvu6W9GjaqJn9Gjp3dyPRl6umvJigvcPi1gNuCprGYOHGiWrRooaSkJLt16NBBs2fPLnjc6/Vq6NChqlq1qipUqKBrr71We/bscbVmAACA4zF5FCaXIsefa3Mq8nLydF67uurep/kJXxMdHaV7RvdWbGyUMjd5Fb08Smv3rNP762eWae1A0DcWtWvX1qOPPqrly5dr2bJluuSSS3TVVVdp3bp19vERI0bogw8+0IwZMzR//nzt3LlT11xzjdtlAwAAHOPVla/bXIq4L2JtTkVSUpyG/LGHIiIiTvq6M+pW1q2DOtn9wCcBaa/02spp2npgWxlVDoRAY3HFFVfosssuU+PGjXXWWWfpn//8p52ZWLRokVJTU/Xiiy/qySeftA1HmzZtNGXKFC1cuNA+fjI+n09paWmFNgBAcMvNzdVf/vIX1a9fX/Hx8WrYsKH+8Y9/sIoOygWTQzFzwyybS2HzKSTdNbK7KlVOKNLre1zZUq3a1FZudp4SPvAox5+jp74cL38uqdxwV9A0FkcPGNOmTVNGRoY9JcrMYmRnZ6t79+4Fzzn77LNVt25dpaSknPR7jR07VsnJyQVbnTp1yuAdAABK07/+9S97Cu2zzz6rDRs22NuPPfaYxo8f73ZpCHMmf8LkUMj7Sy5FQOrW82y169iwyN/DzGrc/ceeqlgxTpnbvYpdEKttB3/Q1JXTSrV2IKQaizVr1thZiri4ON1111165513dM4552j37t2KjY1VpUqFL3aqUaOGfexkxowZY2c88rft27eX8rsAAJQ2M2NtTpft06ePzjzzTF133XXq0aOHlixZ4nZpCGNmxszkT5gcivg5HptLYfIp+g/pUuzvVblqogYPv8TuZ8/3S9uk99d/qDW715ZC5UAINhZNmjTRqlWrtHjxYg0ZMkT9+/fX+vXrHX1P06TkXxCevwEAglvHjh01d+5cffPNN/b26tWrtWDBAvXu3fuEr+HUWJQ2kzth8idMDoXNo/glXTs+IfaUvl/7ixqrS/cmMmf4eT6IU8Ab0NMLJuiQP6PEawdCrrEwsxKNGjWy11CYU5hatmypp59+WjVr1pTf79fBgwcLPd+sCmUeAwCEl9GjR+uGG26wp8XGxMSodevWGj58uPr163fC13BqLEqTyZswuRMmf8LmUEi65vrzbD6FE7cN7aLqp1WQd59P8XM92pf5s14wxwFcEFSNxdHy8vLsJ0ym0TADh/l0Kt+mTZv0ww8/2GswAADh5Y033tDUqVP1+uuva8WKFXrllVf0xBNP2K8nwqmxKC0mZ8LkTZjcCZM/YXIoGjSqputuucDx905IjNM9D5jVpKSspV5FbJDmb1mgBVsWlkjtQHFEK0iYH/hmCttckJ2enm4Hi88//1wff/yx/WRp4MCBGjlypKpUqWJPZ7rnnntsU3HBBc7/0QIAgsv9999fMGthNG/eXNu2bbOzEuY02hOdGms2oKSZax9M3oTJnTD5EyaH4t5RvWwuRUlo2qK2rurbWu++sVIxH8bIXztbkxa/oKY1zlbVhColcgwgpGYs9u7dq1tvvdVeZ9GtWzctXbrUNhWXXnqpffypp57S5ZdfboPxOnfubE+Bevvtt90uGwDggszMTEVGFh7ioqKi7Ew3UJZMvsRrK/9n8yZs7oSkW+7opDPqlewv/L/v31Fn1q8i/6FsJczy6JAvw64+RSo3ylLQzFiYnIqT8Xg8mjBhgt0AAOHNZB+ZvCMzy33uuedq5cqVNuvo9ttvd7s0hBGTK2HyJUzOhMmbyMz2quV5tdXzqpYlfqyYmCjdO7q3Rg37nzI3eBW1IkqrI9bow40f6Yqml5X48YCgnrEAAKCoTF6FWWL27rvvVtOmTfXHP/5RgwcPtiF5QFl5feV0my9hciZM3kSFCrG6+/7fTtc+VXXqV1W/gR0P3/g4IO2T/rt8qrYf/LFUjgccjcYCABByKlasqHHjxtnrKrKysrR582Y9/PDDdnVBoCys2b1O762fafMlbM6EpMHDu6lK1QqletzeV7dW81ZnKNd/OJXbBAg/+eUzys7NKdXjAgaNBQAAQAkyORJPL3jW5krYfImAdHG3s3RB58alfuzIyAgNvb+HEhNjlbnNq9ivYrTlwFZNW/1GqR8boLEAAAAoQS8secnmSZhcCZMvUb16om4b1rXMjl+1ekUNuvfw8bI/z5G2S2+vfU/r92wssxoQnmgsAAAASsiCrQs1//svbZ6EzZWIkO4Z1VOJiWW7lHGnrk10UZdGCuQFFP9+nAK+gMZ9NV6Z/swyrQPhhcYCAACgBPycuV+TFr0gpcvmSRhXXtfa5ky4YeC93VS1WqKyfvLJM8+jvYd+0n+WvuxKLQgPNBYAAAAOmbyIZ756zuZHmBwJkydRr34VXd+/g2s1JVaI0zC7CpXkXeyVNknzNn+ulG2LXasJoY3GAgAAwKFZGz/S6l1f2/wIkyMRExOpP4zurZhYdyPDmrWuo8t/dzg3I87MohySnls0WfszD7haF0ITjQUAAIADJifivyum2twImx8hqd/tHW2uRHlw4+2dVLdeZfnSspXwkUfp3nQ9mzJJAbNcFVCCaCwAAABOkcmHMDkRfn+2zY0w+RHNW9ZS79+1VnlhZk1MKnd0dKQy13oVtSpKK3as1OxNn7hdGkIMjQUAAMApMvkQJifC5EWY3AiTHzH0gZ42T6I8qdewmm4ccIHdj/jIXGkuvbz8Ve1I3el2aQghNBYAAACnwORCvLPuPZsTYfMiJA26p4vNkSiPLr+ujc5tfrpyfLlKmOmR3+/XUwvGKyePVG6UDBoLAACAYjJ5ECYXIs8bULxJ184L6MIujdTpkrNVXtlU7gd6Kj4hRplbvIpNidF3P2/WG1+/5XZpCBE0FgAAAMX04tKXbS6EyYfI2utT1aoJGnjPJSrvqtdI0h1Du9j9nM9ypB3Sm2ve1qafvnG7NIQAGgsAAIBiMDkQczd/bnMhbD6EpKEP9FCFih4Fg4u6n60OFzVQXu7h2ZY8X8CeEpWVffi9AKeKxgIAAKCITP6DyYEweRA2F8Jcu/C7lmreuq6CRUREhAb9obsqV4lX1m6fPJ/HaXf6Hr207BW3S0OQo7EAAAAoApP7YPIfTA6EyYMwuRB16lXSjQM7KdhUTPJo6B972H3vQp/0rTTn27lasn2Z26UhiNFYAAAAFMFH38yx+Q8mB8LkQZhciHtH9Vasy+nap6pl23rqfWVzux83M0bKkCakPK+DWalul4YgRWMBAADwG0zew5Rl/5X2/5IDIemG/u11ZqPqCmb9Bl2kM2ony5earYSPPUrNStUEUrlximgsAAAATsLkPJiLm03ug0nXNjkQ5zSraXMhgl1c3OFU7qioCGV+7VXk15Fa+uNye1oUUFw0FgAAACcx4+u3bd5DTMrh/AeTAzFsVC9FRYXGr1ENzjpN19/S3u5HzY6QDsheyL0rbbfbpSHIhMa/CAAAgFJg8h1mrHnb5j3kmtwHSQPvvtjmQYSSK69vqyZNayjbeziV2+v3adyC8crNy3W7NAQRGgsAAIDjMLkO5hSoPF/e4byH3IAuuLCBOl/aVKHGzL7cM7qXPJ5oZW72KmZxtDbt+1Zvrn3H7dIQRGgsAAAAjsNcrG3yHUzOg8l7qFw5XncO72ZzIEJRjdOTdfvdne1+7qe50i5p+uo39e2+79wuDUGCxgIAAOAoJs/hk28/tfkONufBpGvf30MVk+IVyrr0PFfndzjzcCr3+x7l+fM0bsGz8uUc/jMATobGAgAA4Agmx8HkOSjzl3wHyeY9mNyHUGdmY+4c0V3JlTzK2uVV3Pw47UjbqZeXv+p2aQgCNBYAAAC/MPkNz6U8b/McbLp2arbNeTB5D+EiuVKC7r7vUrvvW+CTNkuzN32i5TtWul0ayjkaCwAAgF98+t08Lflxmc1zMLkOJt/B5DyYvIdwcl77+urR5xy775kZa2dvnl04UWneNLdLQzlGYwEAACDZ3IYXl75scxxsnoOk3998vs15CEe33HmxTq+VJO8BvxLmeHQg66CeWzSZVG6cEI0FAAAIeyavweQ2mPwGk+Ng8hyaND1NV93QTuHKEx+je0b1UmRkhDJXehW5NkKLfliieZvnu10ayikaCwAAEPbeWvuuzW0w+Q0mx8HkOdwzqnfIpGufqsZNa+q6m9ra/ahZUdJB6T9Lp2hP+l63S0M5FN7/WgAAQNj7bt9mm9dgchtsfoOk24Z0Vo1ayW6XVi5c06+9GjepruzMHCV86FGWL0vjvnpWuXl5bpeGcobGAgAAhC2Tz2DStXP9uYdzG3IDatfhTHXtda7bpZWvVO5R5gL2KGV+61X00mht2LtR76573+3SUM7QWAAAgLD1yvLXbE6DyWswuQ0mv2HwiO4hm659qk6vXUn9Bx9ecjdvTp60R/rf6un6/uctbpeGcoTGAgAAhKUVO1Zp1qaPbU6DzWuQNGSkCYdLcLu0cql7n+Y6r11d5eXk2dmdHH+une3x5fjdLg3lBI0FAAAIOyaPYfzC56SsX3IaJF162Tlqc0EDt0srt8wszpA/9lBSUpyydngV90Wstqf+qFdXvu52aSgnaCwAAED4pWsvmmxzGeI/9ticBpPXcOvgi90urdyrVDlBd43sbvd9X/qlrdLMDbO0eufXbpeGcoDGAgAAhJXPvp9v8xgi1kQoa5XX5jSYvAaT24Df1q5jQ3XrebYUkDwfxEpe6ZmFzyndd8jt0uAyGgsAABA2TP7CC0umSKlS9Owoe5/JaTB5DSi6/nd3UY2aFeX92a/4OR79nLlfzy/+D6ncYY7GAgAAhAWTu2DyF0wOg03XzsxRo7Oq25wGFE98fKyd5YmIjFDWcq8i1kdowdaF+mLLArdLg4toLAAAQFh4b/37Nn/B5DCYPAaTy3DvaNK1T1WTc0/XNdefZ/djPoyS0qTnF7+onw7tc7s0uIR/SQAAIOSZvIXXV023+Qs2h8GcznPnRTafAafuulsuUING1eTPOJzKnenL1DMLJygvQCp3OKKxAAAAIc3kLJi8BZO7kGDStXPybB5D98ubu11a0IuOjtK9o3opNjZKmZu8il4erTW71+n99R+6XRpcQGMBAABC2msrX7d5CyZ3IXOHVxWT4nTXfZeSrl1CzqhXRbfc0cnuBz7Jk/aaP/P/aeuBH9wuDWWMxgIAAIQsk6/wwYZZNm/B5i6YdO0R3VW5SqLbpYWUnle1VMvzais3O08JH5hU7hw99eUzys7Ndrs0lCEaCwAAEJJMroLJVzA5CzZvISBd0uNstevU0O3SQo6Z/bn7/h6qUCFWmdu9il0Qq20Hf9DUldPcLg1liMYCAACEpMmL/2PzFUzOgslbqFGjogbc3cXtskJWlaoVNHh4N7ufPd8v/WBW4pppr7lAeKCxAAAAIWf+9wv05daFNl/B5iz8kq4dnxDrdmkh7YLOjXVxt7NkcvLi349TwBvQM19NUIY/0+3SUAZoLAAAQEgxOQomBdrkKth8BUm/+/15atLsdLdLCwu3D+uq6tUTlbXPJ89cj37K2KfJS150uyyUARoLAAAQMkx+gslRMHkKCbM8Nl+hQaOq6nvrBW6XFjYSEuN0z6ieMotueZd6FbHRzCB9aZO5EdpoLAAAQMgwK0CZc/pNnkLmRq9iYqN0z6jeNm8BZadpi9q68rrWdj/mwxgpXZq06AV7zQtCF40FAAAICSY34dUVr9scBZunIOmWOzqqdr0qbpcWlq4f0FH16leRPz3bzh4d8mXoma+eI5U7hNFYAACAoGfyEkxugslPMDkKJk/B5Cr0uqqV26WFrZiYKP1hdG/FxEQqc4NXUSuitHrX15q18SO3S0MpobEAAABBb+qq6TY3weQnmBwFk6dgchVI13ZXnfpV1e/2jnY/4hNJ+6T/rpiq7Qd/dLs0lAIaCwAAENTW7l6v99Z9YHMTsr84nK595/BuNlcB7uv9u9Zq3rKWcny5djbJ78/WUwvGKzs3x+3SUMJoLAAAQNAy+QhPf/WsAr7A4dyEPNkchQ6dG7tdGn4RGRmhoQ/0VGJirDK3eRXzVYy+379F01a/4XZpCNfGYuzYsWrXrp0qVqyo0047TVdffbU2bdpU6Dler1dDhw5V1apVVaFCBV177bXas2ePazUDAIDS9cKSl2xOgudTj81NMPkJtw3r6nZZOErV6hU16J7Dqec5n+dIP0rvrHtP6/dsdLs0hGNjMX/+fNs0LFq0SHPmzFF2drZ69OihjIyMgueMGDFCH3zwgWbMmGGfv3PnTl1zzTWu1g0AAErHV1tT9Pn3X9icBJuXECENG2U+GY9zuzQcR6dLztaFXRopkHd4dinPG7CzTZmkcoeMaAWJjz4qvILAyy+/bGculi9frs6dOys1NVUvvviiXn/9dV1yySX2OVOmTFHTpk1tM3LBBQTjAAAQKkwewsRFk20+gslJ8CtbV17bSue0qO12aTiJO+7tpg1rdurnvZnyzIvTnt579eLSl3VPp7vdLg3hNGNxNNNIGFWqHF6b2jQYZhaje/fuBc85++yzVbduXaWkpJzw+/h8PqWlpRXaAABAOU/X/uo5m4uQMNtjcxJMXoLJTUD5llghTkMf6GH3vYt90iZp7ubPlbJtsdulIVwbi7y8PA0fPlydOnVSs2bN7H27d+9WbGysKlWqVOi5NWrUsI+d7NqN5OTkgq1OnTqlXj8AADh1szZ+bPMQTC5C5nqvzUm41+QlxAbNiRhhrXnrurr8dy3tfpxJ5T4kPbdosg5kHXS7NIRjY2GutVi7dq2mTZvm+HuNGTPGzn7kb9u3by+RGgEA7tqxY4duvvlmu6BHfHy8mjdvrmXLlrldFhwy+Qf/XfGazUOwuQiSbrqtg+rWr+p2aSiGGwd2Up16leRLy1bCRx6le9M1fuFEBQIBt0tDODUWw4YN08yZM/XZZ5+pdu1fz6OsWbOm/H6/Dh4s3O2aVaHMYycSFxenpKSkQhsAILgdOHDAzmrHxMRo9uzZWr9+vf7973+rcuXKbpcGB0zugck/MDkIJg/B5CKYfITLrjnP7dJQTLGx0bp3VG9FR0cqc61XUauitGLHSn30zRy3S0M4NBamgzVNxTvvvKN58+apfv36hR5v06aNHUDmzp1bcJ9ZjvaHH35Qhw4dXKgYAOCWf/3rX/bUVrOIx/nnn2/HDLOSYMOGDU/4Gq65K/+mfz3D5h+YHASTh2ByEUw+gslJQPA5s1F13dC/vd2PMGv07JemLPuvdqTtdLs0hHpjYU5/eu211+yqTybLwlw3YbasrCz7uLk+YuDAgRo5cqSdzTAXc9922222qWBFKAAIL++//77atm2rvn372hUEW7durRdeeOGkr+Gau/Jtw96Nenvtuzb/wOYgmBWGhnWx+QgIXpdf10bnNKt5RCq3X+O+HK+cPFK5g1HQNBYTJ06010B06dJFp59+esE2ffr0guc89dRTuvzyy20wnlmC1pwC9fbbb7taNwCg7H3//fd23GjcuLE+/vhjDRkyRPfee69eeeWVE76Ga+7KL5NzMG7Bs8orSNcOqNPFjXRht7PdLg0ORUVFatioXopPiFHmFq9iUmL07c+bNeNrfn8LRkGzfEJRLubxeDyaMGGC3QAA4cusHmhmLB555BF728xYmEU/Jk2apP79+5/wmjuzofx5cdkr2nNorzxz45S116eqVRN0x72HM6sQ/KrXSNLAuy/Ws098qtzPcqT60oyIt3XeGa3UpPpZbpeHUJyxAACgqMyM9jnnnFPoPhOYaq67Q3BZ9MMSzf3uM5t3YHMPzOnR9/dQhYoet0tDCep8aVNdcGED5eUGFP9BnPJ8eXaWKivb63ZpKAYaCwBAyDErQpkFPI70zTffqF69eq7VhOIzuQbPpTxvcw5s3oGkPle3UPPz6rpdGkpYRESE7hzeXZUrxytrt09xn8VpV/puezE3ggeNBQAg5IwYMUKLFi2yp0J99913duGPyZMn24VAEDynQJtcgzRvuuI/8ti8A5N7cNMdF7pdGkpJxSSPnY0yfCk+6Vvpk28/1ZLt5M8ECxoLAEDIadeunV2e/H//+5+aNWumf/zjHxo3bpz69evndmkooo+/mWNzDUy+QZbJOYiOtLkHJv8Aoatl23rqfWVzux/3YayUKU1IeV4Hs1LdLg1FQGMBAAhJZpXANWvWyOv1asOGDRo0aJDbJaGITI7BlOWv2lwDm29gkppvbW9zDxD6+g26SGfUTpbvoN+mcqdmpdpT4kjlLv9oLAAAQLlh8gtMjoEJLMxP1256bk1d3reN26WhjMTFReve0b0VFRWhzK+9ilwTqSU/LtOn381zuzT8BhoLAABQbpj8ApNjYPIMTK5BfHyM7hndy+YdIHw0OOs0/f7m8+1+1KwI6YD04tKX7QXdKL/4VwoAAMqFb376VjPWvC3t0OE8A0kDh15scw4Qfq66oZ2aND1N2d5cJcz0yOv32SVoc/Ny3S4NJ0BjAQAAXGfyCp5aMN7mF9gcg9yALuhU3+YbIDyZWap7RvWWxxOtzM1exSyO1qafvtFba991uzScAI0FAABw3cvL/2tPc/F8HmdzDEyewZ0jutt8A4SvGrWSdduQznY/99NcaZc0ffWb+m7fZrdLw3HQWAAAAFct/XG5Pv7mU+k7ybvwcLr23X+8VBWT4t0uDeVA117nql2HM39J5fYo159rZ7d8OYf/X0H5QWMBAABcY/IJnl04yeYVxM2Mtff1uqKZWrU70+3SUE6YWavBI7oruZJHWTu9ivs8zi5J/Mry19wuDUehsQAAAK4wuQQmn8DkFNh07YN+m19w86CL3C4N5UxypQQNGdnd7vu+8knfS7M2fawVO1a5XRqOQGMBAABc8el3n9l8ApNTkPW11+YW3DO6t+I8MW6XhnKozQUNdOll59h9zwexUpY0fuFzSvOmu10afkFjAQAAypy5UPvFpVNsPoHNKZDUt187NTzrNLdLQzl26+CLdXqtJHkP+BX/sUcHsg5q4qLJpHKXEzQWAACgTJkcApNHYHIJTD6BySkweQVX33g4EA04EY8JTBzVS5GREcpa5VXE2gil/LBYn30/3+3SQGMBAADK2ttr37N5BCaXwOQTmJwCk1dAujaKonHTmrruprZ2P3pWlJQqvbBkivYc2ut2aWGPf8EAAKDMmPyBaatn2DwCm0sgacBdnW1eAVBU1/Rrr0ZnVVd2Zo6d9cryZenpBROUm5fndmlhjcYCAACUCZM7YPIHTA6BySMwuQTtLjhTl/Q+1+3SEGTM7Na95kL/uChlfutV9NJord+7Qe+tf9/t0sIajQUAACgTJnfA5A/EzY+zeQQml2DwSNK1cWpOr11J/e88vDRx4NM8aY/0+qrp+v7nLW6XFrZoLAAAQKkzeQMmd8DkD/gWHE5MNrkEJp8AOFXdL2+u89rVVW52nuLf9yjnl1Ruf67f7dLCEo0FAAAoVSZnYPzCiTZ3wOYPSLq09zk2lwBwwsx23XXfpaqYFKesHV7FfRGr7ak/6tUVr7tdWliisQAAAKXG5AuYnIEDWQcU/4nH5g+YHIJb77rY7dIQIipXSdSQEb+kcn/pl7ZKH2yYpdU7v3a7tLBDYwEAAErN599/YXMGTN5A1kqvzR8wOQQmjwAoKe06NdQlPc6WAr/MinmlZxY+p0O+Q26XFlZoLAAAQKkwuQKTl7xkcwZs3oCka29sY3MIgJI24O4uqlGjorw/+xU/x6OfM/fr+cX/cbussEJjAQAASpzJEzC5AiZfIOFDj80bMLkDJn8AKA3xCbF2NizCpHIv9ypifYS+3LpQ879f4HZpYYPGAgAAlLj31n9gcwVMvkDmN16bN2B+6YuOPjxzAZSGJs1O1+9+f57djzGzZGmysxY/Zexzu7SwQGMBAABK1Pf7t+r1VdNsroDNF5B0650XqVadym6XhjDQ99YL1KBRNfkP5ShhlkeZvkw989UE5QVI5S5tNBYAAKDEmPyAcQvG2zyBhA88Nl/A5Axcenlzt0tDmDCzYmZ2LCY2SpkbvYpeHq01u9fZlaJQumgsAABAiTH5AT8c3K64L2OV+aPX5guYnAHStVGWaterolvu6Gj3A5/kST8d/n9z64Ef3C4tpNFYAACAErF615rDnwpvlXxfHE4+vmtEN5szAJS1Xle1UsvzattZswSbyp1jZ9Oyc7PdLi1k0VgAAADHTF6AOY/d5AfYHIGA1LXH2Tq/UyO3S0OYMrNkd9/fQxUqxCpzu1exC2K19cA2TV013e3SQhaNBQAAcMysvGNyA0x+gMkRMHkCt93dxe2yEOaqVK2gO4d3s/vZZhbtB+m9dR9o7e71bpcWkmgsAACAI19sWWDzAkxugM0PiIzQsFE9ba4A4LYOnRvr4m5nySwKFf9BnAK+gJ7+6lll+DPdLi3k0FgAAIBTZvIBJi36j80LsLkBkq7u21pnN6vldmlAgduGdVX16onK+sknz6ce+//tCyYVHiWKxgIAAJwSkwtgrqvI9GfavACTG1C/YVX1vbWD26UBhSQmxtlZNLM4mXepVxEbpc+//0JfbU1xu7SQQmMBAABOycwNs2w+QPSyaJsXYHID7h3dWzExpGuj/DmnRW1deW0rux/zYYyULk1cNNleG4SSQWMBAACKbduBH/Tqiv/ZfACbEyDZ3ACTHwCUV9cP6Kh69avIn56thNkeHfJlaPxXE0nlLiE0FgAAoFhMDsBTJg/An23zAUxOQMvWZ6jnlYc/DQbKq5jY6F9m1SKVud6rqJWRWrVrtWZt/Njt0kICjQUAACgWkwNg8gBivzqcD2ByAu5+oKciI0nXRvlXt35V3XTb4euAIj6OkPZJ/13xmrYf/NHt0oIejQUAACgys/6/yQEweQDZ8w+na5ucAJMXAASLy645T81a1FKOL1cJMz3y+3+ZhcvNcbu0oEZjAQAAisSs+2/W/zc5ADYPIE/qfEljmxMABBMzuzb0gZ5KTIxV5lavYr6K0ff7t2j61zPcLi2o0VgAAIAiMev+m/X/PXM9Ng/A5ALcfs8lbpcFnJJqp1XUHcMOp8PnfJ4j/Si9vfZdbdi70e3SghaNBQAA+E1mvX+z7r9Z/9+7xGvzAA5/4hvndmnAKbuw29nqdHEjBfIOz8Ll+QIat+BZZWVnuV1aUKKxAAAAJ7U/c78mLXrBrvtv1/+XdMW1rXRuy9pulwY4dse9l6hq1QRl7fHJMzdOew7t1X+Wvux2WUGJxgIAAPxGuvZEpfsO2XX/zfr/JgfghgEd3S4NKBEVKno09P4edt+72Cd9I8397jMt+mGJ26UFHRoLAABwQrM3fWLX+Tfr/Zt1/836/zYHIDba7dKAEtP8vLrqc3ULux9nZuUypOdSnteBrINulxZUaCwAAMBxmXX9X1n+ql3n3673L9n1/00OABBqbrrjQtWpW0m+1GzFz/YozZuuZxdOVCAQcLu0oEFjAQAAjmHW8zfr+pv1/c06/2a9f7Puv1n/HwhFsb+kckdFRyprrVdRq6O0fMdKffzNHLdLCxo0FgAA4BhmPX+zrn/swhi7zr9Z79+sAkW6NkLZmY2q68Zb29v9iI/MygXSlOWvakfaTrdLCwo0FgAAoBCzjr9Zz9+s65/92eEkYrPev1n3Hwh1l/dto6bn1lSO93Aqt8/n07gvxysnj1Tu30JjAQAAjloF6jm7nv/hdO2AOnVuaNf7B8JBVFSkho3qpfj4GGV+71VMSoy+/XmzZm6Y7XZp5R5LOgAAgALmOtVD/gxpvuy6/saSlC26+crnFI7O79hA947upXC044f9+uf/vau0NK/Ckd+fa79mz8uWmkoZ5t8FTorGAgAAFIiKjFSfs3tp2tszCu7Lzs6TzBaGvpz3Tdg2FutW/6if9h5yu4xyISonUj3PutTtMso9GgsAAFDIZU166e0r35O/k18Kz35CypT0H7eLKB8qn5ukA93SFLbipS7NOqtaIsss/xYaCwAAUEiSp6L6t7lZH2z4MCzX8DehaP5Uv9tllBv+qGypilTJk6y46DiFmwpxFfT7lte5XUZQCKrG4osvvtDjjz+u5cuXa9euXXrnnXd09dVXFzxufvg99NBDeuGFF3Tw4EF16tRJEydOVOPGjV2tGwCAYGNOhzJbOHp47qNalrrC7TLKnWEdh6htbXJMECKrQmVkZKhly5aaMGHCcR9/7LHH9Mwzz2jSpElavHixEhMT1bNnT3m94XnREQAAAFBWgmrGonfv3nY7HjNbMW7cOP35z3/WVVddZe/773//qxo1aujdd9/VDTfcUMbVAgAAAOEjqGYsTmbLli3avXu3unfvXnBfcnKy2rdvr5SUlBO+zoSepKWlFdoAAAAAhGljYZoKw8xQHMnczn/seMaOHWsbkPytTp06pV4rAAAAEGpCprE4VWPGjFFqamrBtn37drdLAgAAAIJOyDQWNWvWtF/37NlT6H5zO/+x44mLi1NSUlKhDQAAAECYNhb169e3DcTcuXML7jPXS5jVoTp06OBqbQAAdz366KOKiIjQ8OHD3S4FAEJWUK0KdejQIX333XeFLthetWqVqlSporp169oB4+GHH7a5FabR+Mtf/qJatWoVyroAAISXpUuX6vnnn1eLFi3cLgUAQlpQzVgsW7ZMrVu3tpsxcuRIu//ggw/a2w888IDuuece3XnnnWrXrp1tRD766CN5PB6XKwcAuMGMA/369bPBqZUrV3a7HAAIaUE1Y9GlSxebV3EiZpr773//u90AABg6dKj69OljlyI3M9onY5YfN1s+lh8HgBBuLAAAKKpp06ZpxYoV9lSoojDLj//tb38r9boAIFQF1alQAAAUhVk6/A9/+IOmTp1a5NNhWX4cAJxhxgIoZ2YvnaNlq1fpj/3vUXwM1wcBp2L58uXau3evzjvvvIL7cnNz9cUXX+jZZ5+1pzxFRUUds/y42QAAp4bGAihHsrK9emTQk9IOKTYyWqNuH+F2SUBQ6tatm9asWVPovttuu01nn322Ro0adUxTAQBwjsYCKEdmfvmRbSqML95bpFG3u10REJwqVqyoZs2aFbovMTFRVatWPeZ+AEDJ4BoLoByZ+dbHBftpK9O0Y/9OV+sBAAAoKmYsgHLCl+PTlgXbfr0jS3r93Td1/+33ulkWEDI+//xzt0sAgJDGjAVQTsxO+VSBnQEpQqrdoY6978uZKW6XBQAAUCQ0FkA58f6M2fZr5XOrqM/AK+z+gRUHtTt1r8uVAQAA/DYaC6Ac8Of6tXnBFrt//mXtdX7X8xWZECllSNM/eMvt8gAAAH4TjQVQDsxZ+rnytufZ06B6Xddb0THRqtvxTPvYZ+8vcLs8AACA30RjAZQD78740H5NblJJVWpWsfvdru1mv/68bL9+St/nan0AAAC/hcYCcFl2bra+mf+d3W97WbuC+zt076gIT4SULr3x4TsuVggAAPDbaCwAl322coHytuXZ/cv69im4PyYuRnU61LX7c9/7wrX6AAAAioLGAnDZ22+8b78mNU5StTOqFXqs6+8Onw7109J92p95wJX6AAAAioLGAnBRTl6uNn7+rd0/r3fbYx6/sNeFioiLkFKlN2e/60KFAAAARUNjAbhowZqFyt2aa/cv+/2vp0Hli/XE6oz2te3+nHfnl3l9AAAARUVjAbjoTXMaVECq0KCCatStcdznXHxVF/t19+I9Ss1KK+MKAQAAiobGAnBJbl6u1s3baPdb9TzvhM+7uM/FioiJkA5Ib33yQRlWCAAAUHQ0FoBLUjYsVc7mnBOeBpUvLsGjmm1Pt/sfvzu3zOoDAAAoDhoLwCUzpr9rT4NKqJeoMxqecdLndr7qYvt1Z8oupXsPlVGFAAAARUdjAbggL5CnNfPW2/2WPVr95vO7XnmJFG1iuKV35x1O6QYAAChPaCwAFyz9ZoWyv822+32uP/FpUPniK8Srxnk17f7st+eUen0AAADFRWMBuOCN6e9IeVJ87XjVOetwuvZvufCqi+zX7Qt3KNOfWcoVAgAAFA+NBVDGAoGAVs1da/ebdW9R5NddckU3KcrEcEvvz59dihUCAAAUH40FUMZWfv+1/Jv8RT4NKl+FShVUveVpdn/mW5+UWn0AAACngsYCKGPT3nhbypU8NT0685z6xXptxys62a8/fLVdWdneUqoQAACg+GgsgDI+DWrFnNV2/5zuzRQREVGs1196TQ/7rzawO6BZC5i1AAAA5QeNBVCG1mxbL98Gn92/7PrLiv36ipUrqmrzanb//Te5zgIAAJQfZmV8oNx8mu/NPPxLd1nyJMQVe+bgVE2b8ZZkwrajpaWfLtOyucsKHqt5Rk1dfF0XRUb92u//uOlHfTV7gfLy8grui4+Jt1+3LNgmX45PcdFxpV53OPzdAAAAZ2gsUC6YX1zv7vNHrV16ODSuLJ3b7hxN/PCJMvkF9tsVmw/v5EgfPz3rmMcrVquott3aFdz+1+CxSvs+9bjfK7AzoG17t+usWo1C9u+mWbumeu7Df9NcAAAQBGgsUC6YT8Pd+MXVWLd0vVLTs1QpKaHUj3XtwCs1Of0V5WX/OgNh5G7KlTKkPXv3FLo/K+1wXkXEWRGKrFD4zMXTm9ZUg5rFu/g72P5u1i7doPTUdCVVSnLl+AAAoBQbiy1btujLL7/Utm3blJmZqerVq6t169bq0KGDPB5Pcb8dcIz7356gxMTDp/uUJr/Xp0euGmL3c3MDKgs3dL/Wbkfr2v4K5XxvzpE6vv733aSB19wstz363xE2Bby0+bx+PXDTE3Y/N+fEfy4o/xgzACB8FLmxmDp1qp5++mktW7ZMNWrUUK1atRQfH6/9+/dr8+bNdoDo16+fRo0apXr16pVu1QhppqlIrFj6v7zGRHN6TXGZpiKxYmKpHyc6OqbUj4HSxZgBAOGnSI2F+XQpNjZWAwYM0FtvvaU6deoUetzn8yklJUXTpk1T27Zt9dxzz6lv376lVTMAoBxjzACA8FSkxuLRRx9Vz549T/h4XFycunTpYrd//vOf2rp1a0nWCAAIIowZABCeitRYnGyAOFrVqlXtBgAIT4wZABCeTnlVqL1799rtyPX1jRYtWpREXQCAEMKYAQChr9iNxfLly9W/f39t2LDBrm9vmDXmzb75mpubWxp1AgCCEGMGAISPYjcWt99+u8466yy9+OKLdqUPgqsAACfCmAEA4aPYjcX3339vV/lo1Kh0034BAMGPMQMAwkfhKN8i6Natm1avXl061QAAQgpjBgCEj2LPWPznP/+x58uuXbtWzZo1U0xM4SCrK6+8UuHMm+lVbGSc22UEnaxMr+t/b1mxUa4dP/BL8rc/y6/M9MyCc9Dzcg9f6Jrty1FWhjd8/25ceu/BzO2/t3yMGQAQPordWJhQo6+++kqzZ88+5rFwvRAv/4JE47pmA1ytJdT+PMvqODe0ulnlwbt/fdtuR5t673S7hePfze/bDCyTY4aysvp7Ox7GDAAIH8U+Feqee+7RzTffrF27dtllA4/cwnWAyMzMdLuEkJLjyyqT42R7/WVynFCSmV02f2Z+X3aZHCdcHDpw0LVjM2YAQPgo9ozFzz//rBEjRtjVPXCsCXP/pKpVCHsqrtT96Rrc7SG7H6Gy/3R1zKdPqNppyXLLqx9N1rbvN+vApP3HPNb7+d66sNuFcsuhnw9pbPuxdj8QUfZ/Ny++8WdVrVmlzI8b7A4eOKQBVz1o9wO/nFLnBsYMAAgfxW4srrnmGn322Wdq2LBh6VQU5OLi45VQIcHtMoKO35vj6vHjEzxKdPHv7a7rhittX6oemHS3vT34kWv1/P+9Zfdbt22lxMRE12rLy3Lvl1LDk+BRPP+mis3nd/ffVD7GDAAIH8VuLMx65GPGjNGCBQvUvHnzYy7Eu/fee0uyPiAsVaxeoWCfdf8RzBgzACB8nNKqUBUqVND8+fPtdiTzCxCDBAAgH2MGAISPYjcWW7ZsKZ1KAAAhhzEDAMJHsVeFAlD6MtKyysVSoQAAAKU2Y3H77bef9PGXXnqpuN8SwFENxMT73ijYN6F5QLBizACA8FHsxuLAgQOFbmdnZ9tE1YMHD+qSSy4pydoAmMYik0wHBC/GDAAIH8VuLN55551j7jNBR0OGDGE5QaAUxHgKr6IDBBPGDAAIH8VuLI4nMjJSI0eOVJcuXfTAAw+UxLcE8IscH+nECC3BNmZ4s3yKjOKSxHCSa37uZv56+63XFjk6xdWtZcOdHv/bTXsOf5+8gOQ/fGquN9NbojVCQfEzsEwbC2Pz5s3KySkfgUwIbodS0xQVVfqf0h86mPbr/v40eaLd/cUh48Chgv1hD96mZ/8+5ZdbEUrfl+5q8nbB/sFMxUfElfox01MzCva5dj00BdOYcfsFg90uAS6b9t8lCmeZa73SWumxR550uxSUc8VuLMynTEd3w7t27dKHH36o/v37l2RtCCMB/frb47hb/1Tmx3/s8v9TefJrUyE9felTKi/GDp5U5sf0eov+SQnKH8YMAAgfxW4sVq5cecyUdvXq1fXvf//7N1f/KCsTJkzQ448/rt27d6tly5YaP368zj//fLfLAnAKvNnB8ak2gnfM+C1//WS0PIket8tAGVu0baneXfvhqX8DM/H69C/7JgeygsqWmWx+5pf9P0hKPPVvVSuxuu5re7MqxCaUVHUIIlkZPg2//InSaSw+++wzlWfTp0+3n5BNmjRJ7du317hx49SzZ09t2rRJp512WqkfP23/IcVGxJb6cUJN+v5fT30ZMH24EqtXLvVjmk9O/Vk+5XizlVg50bVzYI+sJzvLr1xvtuKT45XtzZbf71N85Qpys7TDdWUrJ8uvhGRTS+kXk7E/XVN+f3imxsxlHdz/62lrKJrUg4fKRRZKeR8ziiI+PkHxCTQW4eaiszpq/f5NqhSfrLOqNyr2670HvXpT79n9PvWvUoVKZdtZeNO8ek9v2f0bm/VSUrVTO350ZLRaVD9L8dGlfxosyqlAZNlfY1FePPnkkxo0aJBuu+02e9s0GGbK3ayVPnr06GOe7/P57JYvLa34v8AcOWiPvPzRU64dh718/Ti3S0A5cs9NY90uIehlZh5xFSqAIomJitGwTnee8uvTEw8VNBbNK7VUlSpVVJYORaYXNBbn1zhXVWuW7fERnorUgvTq1UuLFv32igjp6en617/+ZU9FcoPf79fy5cvVvXv3QtPu5nZKSspxXzN27FglJycXbHXq1CnDigGg9KUWY0WPkhAsYwYAoGQVacaib9++uvbaa+0v3ldccYXatm2rWrVqyePx2PCj9evXa8GCBZo1a5b69Oljr29ww759+5Sbm6saNWoUut/c3rhx43FfM2bMmEIXF5oZi+I2F0eeGjLmjYdVqSqfChRbICCf16csb7YSKiW7floSXBaQfFle5WV5lVRGp1+FmvQDqfrH9WPsfnL10j8NNBjHDKCsZKQdUmxk2WYSHUo78nTIMj00wliRGouBAwfq5ptv1owZM+w1DJMnT1Zqaqp9zAz455xzjr2OYenSpWratKmCSVxcnN1KSsWqlVW5jAdxADhaZMyv13qVdWMWymMGcCqnSY//g7vLtPp82a4eH+GjyNdYmF++zUBhNsMMEllZWapatapiYspHMnC1atUUFRWlPXsOB7rkM7dr1qzpWl0AEG6CYcwAwoXP53e7BISJU754O/+ahPIkNjZWbdq00dy5c3X11Vfb+/Ly8uztYcOGuV0eAISt8jhmAKXpyJnCG1+5XRUrJ5Xp8TMOZGpq/8l2v0IlVjVD2Qi5VaHM9RImdMmc02uyK8xysxkZGQWrRAEAAJSlStWSVKl62V5/6Yl173RIhK+iL0wbJK6//no98cQTevDBB9WqVSutWrVKH3300TEXdAMAQpdZ8a9du3aqWLGizTAys9gmzwgAUHpCrrEwzGlP27Zts/kUixcvtkF5AIDwMX/+fA0dOtQueztnzhxlZ2erR48edgYbAFA6Qu5UKAAAzEz1kV5++WU7c2Gyjjp37nzc15REYCoAhLNiz1iY6xe++OKL0qkGABBSysuYkb/c7cnSjwlMBYAybizMD2eTZN24cWM98sgj2rFjh8MSAAChqjyMGWZ1wOHDh6tTp05q1qzZCZ9nAlNNvfnb9u3by7ROAAi7xuLdd9+1A8OQIUNs8NGZZ56p3r17680337TnsAIAUJ7GDHOtxdq1azVt2rTfzN5ISkoqtAEASvni7erVq9tlXVevXm0vjm7UqJFuueUW1apVSyNGjNC33357Kt8WABCC3BwzzGIeM2fO1GeffabatWuX2nEAAA5Xhdq1a5ddbcNsJvH6sssu05o1a3TOOefoqaeeKrkqAQBBryzHjEAgYJuKd955R/PmzVP9+vVL9PsDAEqgsTBT12+99ZYuv/xy1atXTzNmzLDnru7cuVOvvPKKPv30U73xxhv6+9//XtxvDQAIMW6NGeb0p9dee02vv/66zbLYvXu33bKyskr0OAAAB8vNnn766fZCuBtvvFFLliyxIXRH69q1qypVqlTcbw0ACDFujRkTJ060X7t06VLo/ilTpmjAgAEleiwAwCk2Fma6um/fvvJ4PCd8jhkgtmzZUtxvDQAIMW6NGeZUKABAOW8szAV3AAAUBWMGAIQPRxdvAwAAAIBBYwEAAADAMRoLAAAAAI7RWAAAAABwjMYCAAAAgGM0FgAAAAAco7EAAAAAUPY5Fji5jANpio2Mc7sMAGEu/UB6wT5hcQCAskBjUQKOHLT/ed0YV2sBgKNlZWa6XQIAIAxwKhQAhLhDBw+6XQIAIAwwY1ECIiIiCvbv/fetSk5KdrUeAEhPy9S4+160+1WSK7hdDgAgDNBYlLDK1ZJVtUoVt8sAEOY8ntjjfvgBAEBp4VQoAAAAAI7RWAAAAABwjMYCAAAAgGM0FgAAAAAco7EAAAAA4BirQgEAAJSizINZiok6VMbHzCjYPyLHFyhVNBYAAAAlLHDEb/Ov9Hve1Vp83mxXj4/wwalQAAAAISzXm+N2CQgTzFgAAACUsCODKf/5yL2qnJRcpsdPS0vX6P8bZ/eTk+LL9NgIXzQWAAAApahSlWRVq1KlTI8ZExdTsB+hX5scoDRxKhQAAAAAx2gsAAAAADhGYwEAAADAMRoLAAAAAI7RWAAAAABwjMYCAAAAgGMsNwsAAFCK0tIOKTYqpsyPme/XDHCgdNFYAAAAlLBA4Ndf58c88JSrtfi82a4eH+GDU6EAAABCmM9HY4GywYwFAABACYuI+DXt+t7Hb1al5KQyPX56Wqae+uPLdj8x2VOmx0b4orEAAAAoRUnVklWpauUyPWaUJ/a4TQ5QmjgVCgAAAIBjNBYAAAAAHKOxAAAAAOAYjQUAAAAAx2gsAAAAADhGYwEAAADAMRoLAAAAAI7RWAAAAABwjMYCAAAAgGM0FgAAAAAco7EAAAAA4BiNBQAAAADHop1/CwAAQtehAxnK8ee4XQaC8P+bgv3UTMVExZTt8dOyCvbTDmQqJvZQmR4foSMrwxd6jcU///lPffjhh1q1apViY2N18ODBY57zww8/aMiQIfrss89UoUIF9e/fX2PHjlV0dNC8TQBAOfOPK//ldgkIcuOGT3H1+A/1f97V4yN8BM1v3H6/X3379lWHDh304osvHvN4bm6u+vTpo5o1a2rhwoXatWuXbr31VsXExOiRRx5xpWYAAAAgXARNY/G3v/3Nfn355ZeP+/gnn3yi9evX69NPP1WNGjXUqlUr/eMf/9CoUaP017/+1c5yHI/P57NbvrS0tFJ6BwCAYPTES3crPiHO7TIQZAIByef1K9ubraQK8YqIiCjjAiSvN1v+7BwlVfKU/fERMrIyfRoxcGJoNRa/JSUlRc2bN7dNRb6ePXvaU6PWrVun1q1bH/d15lSp/KYFAICjVatSRQmJHrfLAABXZGZ4w29VqN27dxdqKoz82+axExkzZoxSU1MLtu3bt5d6rQAAAECocbWxGD16tJ2aO9m2cePGUq0hLi5OSUlJhTYAAAAAxePqqVD33XefBgwYcNLnNGjQoEjfy1y0vWTJkkL37dmzp+AxAAAAACHaWFSvXt1uJcGsFmWWpN27d69OO+00e9+cOXPsDMQ555xTIscAAAAAEOQXb5uMiv3799uvZmlZk2dhNGrUyGZW9OjRwzYQt9xyix577DF7XcWf//xnDR061J7uBAAAAKD0BE1j8eCDD+qVV14puJ2/ypMJw+vSpYuioqI0c+ZMuwqUmb1ITEy0AXl///vfXawaAAAACA9B01iY/IoTZVjkq1evnmbNmlVmNQEAAAAIseVmAQA42oQJE3TmmWfK4/Goffv2xyzyAQAoOTQWAICQNH36dI0cOVIPPfSQVqxYoZYtW9rgVLPIBwAgjE+FAgCgOJ588kkNGjRIt912m709adIkffjhh3rppZdsjlJRpR3MULY/pxQrBYDyKyvTV+Tn0lgAAEKO3+/X8uXLNWbMmIL7IiMj1b17d6WkpBz3NT6fz2750tLS7Nd7BzxVBhUDQPDjVCgAQMjZt2+fXZq8Ro0ahe43t81y5MczduxYJScnF2x16tQpo2oBIDQwYwEAgGRnN8w1GUfOWJjm4u9v3CxPQqyrtQGAW7yZfj34+9eK9FwaCwBAyKlWrZrNN9qzZ0+h+83tmjVrHvc1Jkz1eIGqydUrK74CQasAwlPWoaJfY8GpUACAkBMbG6s2bdpo7ty5Bffl5eXZ2yZEFQBQ8pixAACEJHNaU//+/dW2bVudf/75GjdunDIyMgpWiQIAlCwaCwBASLr++uv1008/6cEHH7QXbLdq1UofffTRMRd0AwBKBo0FACBkDRs2zG4AgNLHNRYAAAAAHKOxAAAAAOAYjQUAAAAAx2gsAAAAADhGYwEAAADAMRoLAAAAAI7RWAAAAABwjMYCAAAAgGM0FgAAAAAco7EAAAAA4BiNBQAAAADHaCwAAAAAOEZjAQAAAMAxGgsAAAAAjtFYAAAAAHCMxgIAAACAYzQWAAAAAByjsQAAAADgGI0FAAAAAMdoLAAAAAA4RmMBAAAAwDEaCwAAAACO0VgAAAAAcIzGAgAAAIBjNBYAAAAAHKOxAAAAAOAYjQUAAAAAx2gsAAAAADhGYwEAAADAMRoLAAAAAI7RWAAAAABwjMYCAAAAgGM0FgAAAAAco7EAAAAA4BiNBQAAAADHaCwAAAAAOEZjAQAAAMAxGgsAAAAAjtFYAAAAAHCMxgIAAACAYzQWAAAAAByjsQAAAADgGI0FAAAAAMdoLAAAAAA4RmMBAAAAIDwai61bt2rgwIGqX7++4uPj1bBhQz300EPy+/2Fnvf111/roosuksfjUZ06dfTYY4+5VjMAAAAQTqIVBDZu3Ki8vDw9//zzatSokdauXatBgwYpIyNDTzzxhH1OWlqaevTooe7du2vSpElas2aNbr/9dlWqVEl33nmn228BAAAACGlB0Vj06tXLbvkaNGigTZs2aeLEiQWNxdSpU+0MxksvvaTY2Fide+65WrVqlZ588kkaCwAAAKCUBcWpUMeTmpqqKlWqFNxOSUlR586dbVORr2fPnrYBOXDgwAm/j8/ns7MdR24AAAAAwqCx+O677zR+/HgNHjy44L7du3erRo0ahZ6Xf9s8diJjx45VcnJywWauzQAAAAAQRI3F6NGjFRERcdLNXF9xpB07dtjTovr27Wuvs3BqzJgxdvYjf9u+fbvj7wkAAACEG1evsbjvvvs0YMCAkz7HXE+Rb+fOneratas6duyoyZMnF3pezZo1tWfPnkL35d82j51IXFyc3QAAAAAEaWNRvXp1uxWFmakwTUWbNm00ZcoURUYWnmzp0KGD/vSnPyk7O1sxMTH2vjlz5qhJkyaqXLlyqdQPAAAAIIiusTBNRZcuXVS3bl27CtRPP/1kr5s48tqJm266yV64bfIu1q1bp+nTp+vpp5/WyJEjXa0dAAAACAdBsdysmXkwF2ybrXbt2oUeCwQC9qu58PqTTz7R0KFD7axGtWrV9OCDD7LULAAAAFAGgqKxMNdh/Na1GEaLFi305ZdflklNAAAAAILsVCgAAAAA5RuNBQAAAADHaCwAAAAAOEZjAQAIKVu3brUrBNavX1/x8fFq2LChHnroIfn9frdLA4CQFhQXbwMAUFQbN25UXl6enn/+eTVq1Ehr167VoEGDlJGRYZcsBwCUDhoLAEBI6dWrl93yNWjQQJs2bdLEiRNP2lj4fD675UtLSyv1WgEglHAqFAAg5KWmpqpKlSonfc7YsWNtJlL+VqdOnTKrDwBCAY0FACCkmXDV8ePHa/DgwSd93pgxY2wDkr9t3769zGoEgFBAYwEACAqjR49WRETESTdzfcWRduzYYU+L6tu3r73O4mTi4uKUlJRUaAMAFB3XWAAAgsJ9992nAQMGnPQ55nqKfDt37lTXrl3VsWNHTZ48uQwqBIDwRmMBAAgK1atXt1tRmJkK01S0adNGU6ZMUWQkE/QAUNpoLAAAIcU0FV26dFG9evXsKlA//fRTwWM1a9Z0tTYACGU0FgCAkDJnzhx7wbbZateuXeixQCDgWl0AEOqYGwYAhBRzHYZpII63AQBKD40FAAAAAMdoLAAAAAA4RmMBAAAAwDEaCwAAAACO0VgAAAAAcIzGAgAAAIBjNBYAAAAAHKOxAAAAAOAYjQUAAAAAx2gsAAAAADhGYwEAAADAMRoLAAAAAI7RWAAAAABwjMYCAAAAgGM0FgAAAAAco7EAAAAA4BiNBQAAAADHaCwAAAAAOEZjAQAAAMAxGgsAAAAAjtFYAAAAAHCMxgIAAACAYzQWAAAAAByjsQAAAADgGI0FAAAAAMdoLAAAAAA4RmMBAAAAwDEaCwAAAACO0VgAAAAAcIzGAgAAAIBjNBYAAAAAHKOxAAAAAOAYjQUAAAAAx2gsAAAAADhGYwEAAADAMRoLAAAAAI7RWAAAAABwjMYCAAAAgGM0FgAAAAAco7EAAAAA4BiNBQAAAADHaCwAAAAAOEZjAQAAACB8Gosrr7xSdevWlcfj0emnn65bbrlFO3fuLPScr7/+WhdddJF9Tp06dfTYY4+5Vi8AAAAQToKmsejataveeOMNbdq0SW+99ZY2b96s6667ruDxtLQ09ejRQ/Xq1dPy5cv1+OOP669//asmT57sat0AAABAOIhWkBgxYkTBvmkeRo8erauvvlrZ2dmKiYnR1KlT5ff79dJLLyk2NlbnnnuuVq1apSeffFJ33nnnCb+vz+ez25ENCgAAAIAQnbE40v79+20j0bFjR9tUGCkpKercubNtKvL17NnTznAcOHDghN9r7NixSk5OLtjMKVQAAAAAQrixGDVqlBITE1W1alX98MMPeu+99woe2717t2rUqFHo+fm3zWMnMmbMGKWmphZs27dvL8V3AAAAAIQmVxsLczpTRETESbeNGzcWPP/+++/XypUr9cknnygqKkq33nqrAoGAoxri4uKUlJRUaAMAAAAQRNdY3HfffRowYMBJn9OgQYOC/WrVqtntrLPOUtOmTe1pS4sWLVKHDh1Us2ZN7dmzp9Br82+bxwAAAACEaGNRvXp1u52KvLw8+zX/wmvTXPzpT38quJjbmDNnjpo0aaLKlSuXYNUAAAAAgvIai8WLF+vZZ5+1qzxt27ZN8+bN04033qiGDRvahsK46aab7IXbAwcO1Lp16zR9+nQ9/fTTGjlypNvlAwAAACEvKBqLhIQEvf322+rWrZudgTDNQ4sWLTR//nx7jYRhVnQy115s2bJFbdq0sadZPfjggyddahYAAABAGOVYNG/e3M5S/BbTbHz55ZdlUhMAAACAIJuxAAAAAFC+0VgAAAAAcIzGAgAAAIBjNBYAAAAAHKOxAACELJN11KpVK0VERNglywEApYfGAgAQsh544AHVqlXL7TIAICwExXKzAAAU1+zZs22+0VtvvWX3T1X6gUxl+3NKtDYACBbeDH+Rn0tjAQAIOXv27NGgQYP07rvv2pDVop42ZbZ8aWlp9uufr3mx1OoEgFDCqVAAgJASCAQ0YMAA3XXXXWrbtm2RXzd27FglJycXbHXq1CnVOgEg1DBjAQAICqNHj9a//vWvkz5nw4YN9vSn9PR0jRkzpljf3zx/5MiRhWYsTHPxzxcHy5MQd8p1A0Aw82b69KeBzxfpuTQWAICgcN9999mZiJNp0KCB5s2bp5SUFMXFFW4GzOxFv3799Morrxz3teb5R7/GqFy9iuITPQ6rB4DglJXhLfJzaSwAAEGhevXqdvstzzzzjB5++OGC2zt37lTPnj01ffp0tW/fvpSrBIDwRWMBAAgpdevWLXS7QoUK9mvDhg1Vu3Ztl6oCgNDHxdsAAAAAHGPGAgAQ0s4880y7UhQAoHQxYwEAAADAMRoLAAAAAI5xKlQJy0jLVGxkrNtlAAhzh1IzC/bz8jgNCABQ+mgsStgT97zodgkAUMjBAwfcLgEAEAY4FaoEVKtWrWA5QwAob+I9CW6XAAAIA8xYlIDIyEilpqbqm80blOnNcLscALCnP5mZCtNUdDj/QrfLAQCEARqLEmwuzm58rttlAAAAAK7gVCgAAAAAjtFYAAAAAHCMxgIAAACAYzQWAAAAAByjsQAAAADgGI0FAAAAAMdoLAAAAAA4RmMBAAAAwDEaCwAAAACO0VgAAAAAcIzGAgAAAIBjNBYAAAAAHKOxAAAAAOAYjQUAAAAAx2gsAAAAADhGYwEAAADAMRoLAAAAAI7RWAAAAABwjMYCAAAAgGM0FgAAAAAci3b+LUJLIBCwX9PS0twuBQBck/8zMP9nYjjKf+9ZmT63SwEA1+T/DCzKeEBjcZT09HT7tU6dOm6XAgDl4mdicnKywnk8uO+GJ90uBQCCYjyICITzx1HHkZeXp507d6pixYqKiIhQMHyqaJqg7du3KykpSaEqXN6nwXsNTcH2Xs3QYAaRWrVqKTIyPM+adXs8CLb/Z0pDuP8ZhPv7N/gzkOt/BsUZD5ixOIr5A6tdu7aCjfkfLRz+wYXL+zR4r6EpmN5ruM5UlLfxIJj+nykt4f5nEO7v3+DPQK7+GRR1PAjPj6EAAAAAlCgaCwAAAACO0VgEubi4OD300EP2aygLl/dp8F5DUzi9V5QM/p/hzyDc37/Bn4GC6s+Ai7cBAAAAOMaMBQAAAADHaCwAAAAAOEZjAQAAAMAxGgsAAAAAjtFYhACfz6dWrVrZZNhVq1YVeuzrr7/WRRddJI/HY1MbH3vsMQWTrVu3auDAgapfv77i4+PVsGFDuzKC3+8Pqfd5pAkTJujMM8+076V9+/ZasmSJgtnYsWPVrl07m1582mmn6eqrr9amTZsKPcfr9Wro0KGqWrWqKlSooGuvvVZ79uxRsHv00Uftv8vhw4eH/HuF+z/vQ1lRx4JQE2rjQUmPHeHk0eOMJ+URjUUIeOCBB2zM+vEi4Hv06KF69epp+fLlevzxx/XXv/5VkydPVrDYuHGj8vLy9Pzzz2vdunV66qmnNGnSJP3f//1fSL3PfNOnT9fIkSPtgLlixQq1bNlSPXv21N69exWs5s+fb3+RXrRokebMmaPs7Gz795WRkVHwnBEjRuiDDz7QjBkz7PN37typa665RsFs6dKl9v/bFi1aFLo/FN8r3P95H+qKMhaEmlAcD0p67AgXS08wnpRLZrlZBK9Zs2YFzj777MC6devMssGBlStXFjz23HPPBSpXrhzw+XwF940aNSrQpEmTQDB77LHHAvXr1w/J93n++ecHhg4dWnA7Nzc3UKtWrcDYsWNdrask7d271/6/On/+fHv74MGDgZiYmMCMGTMKnrNhwwb7nJSUlEAwSk9PDzRu3DgwZ86cwMUXXxz4wx/+ELLvFeXj5304OnosCDXhMB44GTvCRfoJxpPyihmLIGZOnxg0aJBeffVVJSQkHPN4SkqKOnfurNjY2IL7zKcdZirxwIEDClapqamqUqVKyL1PM6VvZly6d+9ecF9kZKS9bd5jqDB/f0b+36F5z+aTqCPf99lnn626desG7fs2n7L16dOn0HsK1feK8vHzPhwdPRaEknAZD5yMHeFi6AnGk/KKxiJImVzDAQMG6K677lLbtm2P+5zdu3erRo0ahe7Lv20eC0bfffedxo8fr8GDB4fc+9y3b59yc3OP+16C6X2cjDmVwZwf2qlTJzVr1szeZ96baQorVaoUEu972rRp9rQFc37w0ULtvaL8/LwPN8cbC0JJOIwHTseOcDDtJONJeUVjUc6MHj3aXpxzss2ca2p+oKanp2vMmDEK5fd5pB07dqhXr17q27ev/eQOwfnJy9q1a+0Py1C0fft2/eEPf9DUqVPtxZbAyYTLz/uTYSxAUYT62BFK40m02wWgsPvuu89+MnUyDRo00Lx58+x0aFxcXKHHzKdZ/fr10yuvvKKaNWses9pM/m3zWDC8z3zmAteuXbuqY8eOx1yUXZ7fZ3FUq1ZNUVFRx30vwfQ+TmTYsGGaOXOmvvjiC9WuXbvgfvPezLT/wYMHC32SH4zv25y6YC6sPO+88wruM586mvf87LPP6uOPPw6Z94ry9fM+WJXkWBBKQn08KImxI9Qt/43xxKwQZ/4fKXfcvsgDp2bbtm2BNWvWFGwff/yxvajpzTffDGzfvr3QRc1+v7/gdWPGjAm6i5p//PFHe+HSDTfcEMjJyTnm8VB5n/kX6w0bNqzQxXpnnHFGUF+sl5eXZy9ANBcdfvPNN8c8nn9Bs/l/N9/GjRuD8oLmtLS0Qv8uzda2bdvAzTffbPdD6b2ifP28Dwe/NRaEmlAcD0py7Ah1ab8xnpRXNBYhYsuWLcesEmJ+ialRo0bglltuCaxduzYwbdq0QEJCQuD5558PBNNA0qhRo0C3bt3s/q5duwq2UHqf+UztcXFxgZdffjmwfv36wJ133hmoVKlSYPfu3YFgNWTIkEBycnLg888/L/T3l5mZWfCcu+66K1C3bt3AvHnzAsuWLQt06NDBbqHg6FU8Qvm9wr2f96GuKGNBqAnF8aCkx45wc3EQrApFYxHiA83q1asDF154of3hZD7pePTRRwPBZMqUKfZ9HW8Lpfd5pPHjx9tfPGNjY+0nVosWLQoEsxP9/Zm/23xZWVmBu+++2848mabwd7/7Xcj8wnD0QBDK7xVlIxwbi6KOBaEm1MaDkh47ws3FQdBYRJj/uH06FgAAAIDgxqpQAAAAAByjsQAAAADgGI0FAAAAAMdoLAAAAAA4RmMBAAAAwDEaCwAAAACO0VgAAAAAcIzGAgAAAIBjNBaAC1588UX16NGj1I/z0UcfqVWrVsrLyyv1YwEAio/xAKGExgIoY16vV3/5y1/00EMPlfqxevXqpZiYGE2dOrXUjwUAKB7GA4QaGgugjL355ptKSkpSp06dyuR4AwYM0DPPPFMmxwIAFB3jAUINjQVwin766SfVrFlTjzzySMF9CxcuVGxsrObOnXvC102bNk1XXHFFofu6dOmi4cOHF7rv6quvtoNAvjPPPFMPP/ywbr31VlWoUEH16tXT+++/b+u46qqr7H0tWrTQsmXLCn0fcyxz3+bNm0vgXQMAjsZ4ABxGYwGcourVq+ull17SX//6V/uDOj09XbfccouGDRumbt26nfB1CxYsUNu2bU/pmE899ZT9ZGvlypXq06ePPZ4ZWG6++WatWLFCDRs2tLcDgUDBa+rWrasaNWroyy+/PKVjAgBOjvEAOIzGAnDgsssu06BBg9SvXz/dddddSkxM1NixY0/4/IMHDyo1NVW1atU65eMNHjxYjRs31oMPPqi0tDS1a9dOffv21VlnnaVRo0Zpw4YN2rNnT6HXmeNt27btlI4JAPhtjAcAjQXg2BNPPKGcnBzNmDHDXhQXFxd3wudmZWXZrx6P55SOZaa285lPnYzmzZsfc9/evXsLvS4+Pl6ZmZmndEwAQNEwHiDc0VgADplzVXfu3GmX8Nu6detJn1u1alVFRETowIEDhe6PjIwsNF1tZGdnH/N6s6JHPvN9TnTf0csJ7t+/307VAwBKD+MBwh2NBeCA3++357Nef/31+sc//qE77rjjmE+HjmQu5DvnnHO0fv36QvebH/K7du0quJ2bm6u1a9eW2HKGZrBr3bp1iXw/AMCxGA8AGgvAkT/96U/2HFmzfJ85n9Wc13r77bef9DU9e/a0F+wd6ZJLLtGHH35ot40bN2rIkCH2/NuSsGjRIjsd36FDhxL5fgCAYzEeADQWwCn7/PPPNW7cOL366qt2HXIzfW32zWobEydOPOHrBg4cqFmzZtkBKJ8ZfPr3729X8Lj44ovVoEEDde3atUTq/N///mcvJkxISCiR7wcAKIzxADgsInD0iXwASp1ZteO8887TmDFjSvU4+/btU5MmTezyh/Xr1y/VYwEAio/xAKGEGQvABY8//rgNMCpt5uLB5557jkEEAMopxgOEEmYsAAAAADjGjAUAAAAAx2gsAAAAADhGYwEAAADAMRoLAAAAAI7RWAAAAABwjMYCAAAAgGM0FgAAAAAco7EAAAAA4BiNBQAAAAA59f+CjswZhTUKFQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mThe Kernel crashed while executing code in the current cell or a previous cell. \n",
      "\u001b[1;31mPlease review the code in the cell(s) to identify a possible cause of the failure. \n",
      "\u001b[1;31mClick <a href='https://aka.ms/vscodeJupyterKernelCrash'>here</a> for more info. \n",
      "\u001b[1;31mView Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "class InP_EOPM_tree:\n",
    "    def __init__(\n",
    "            self,\n",
    "            **kwargs\n",
    "    ):\n",
    "        \n",
    "        self.e = 1.60e-19 # electron charge in C\n",
    "        self.e0 = 8.85e-12 # vacuum permittivity in F/m\n",
    "        \n",
    "        self.w_sig_metal = 5 # Width of signal metal in um\n",
    "        self.metal_sep = 10 # Separation between signal and ground metals in um\n",
    "        self.h_metal = 4 # Height of metals in um\n",
    "        self.w_gnd_metal = 10\n",
    "        \n",
    "        self.w_wg = 1\n",
    "        self.h_n = 0.4\n",
    "        self.h_wg1 = 0.5\n",
    "        self.h_wg2 = 0.3\n",
    "        self.h_p1 = 1\n",
    "        self.h_p2 = 0.2\n",
    "\n",
    "        self.h_box = 4\n",
    "\n",
    "        self.w_window = 100\n",
    "        self.h_bottom = 30\n",
    "        self.h_top = 30\n",
    "\n",
    "        self.base_tree = 10\n",
    "        self.height_tree = 50\n",
    "        self.width_tree = 20\n",
    "\n",
    "        for kwarg, value in kwargs.items():\n",
    "            if hasattr(self, kwarg):\n",
    "                setattr(self, kwarg, value)\n",
    "    \n",
    "    def _make_meshes(self):\n",
    "        # optical mesh\n",
    "        self.optical_mesh_settings = {\n",
    "            'substrate': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'background': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'box': {'resolution': 0.3, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'sig_metal': {'resolution': 10, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'n_metal_left': {'resolution': 10, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'n_metal_right': {'resolution': 10, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'bcb': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'n': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'wg1': {'resolution': 0.1, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'wg2': {'resolution': 0.1, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'p1': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'p2': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "        }\n",
    "\n",
    "        # RF mesh\n",
    "        self.rf_mesh_settings = {\n",
    "            'substrate': {'resolution': 5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'background': {'resolution': 5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'box': {'resolution': 3, 'SizeMax': 3, 'distance': 0.1},\n",
    "            'sig_metal': {'resolution': 3, 'SizeMax': 3, 'distance': 0.1},\n",
    "            'n_metal_left': {'resolution': 3, 'SizeMax': 3, 'distance': 0.1},\n",
    "            'n_metal_right': {'resolution': 3, 'SizeMax': 3, 'distance': 0.1},\n",
    "            'bcb': {'resolution': 5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'n': {'resolution': 5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'wg1': {'resolution': 1, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'wg2': {'resolution': 1, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'p1': {'resolution': 5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'p2': {'resolution': 5, 'SizeMax': 5, 'distance': 0.1},\n",
    "        }\n",
    "\n",
    "        # eo mesh\n",
    "        self.eo_mesh_settings = {\n",
    "            'substrate': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'background': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'box': {'resolution': 0.3, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'sig_metal': {'resolution': 10, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'n_metal_left': {'resolution': 10, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'n_metal_right': {'resolution': 10, 'SizeMax': 0.2, 'distance': 0.1},\n",
    "            'bcb': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'n': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'wg1': {'resolution': 0.1, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'wg2': {'resolution': 0.1, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'p1': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "            'p2': {'resolution': 0.5, 'SizeMax': 5, 'distance': 0.1},\n",
    "        }\n",
    "\n",
    "        self.charge_mesh_settings = {\n",
    "            'substrate': {'resolution': 0.5},\n",
    "            'background': {'resolution': 0.5},\n",
    "            'sig_metal': {'resolution': 0.01},\n",
    "            'n_metal_left': {'resolution': 0.01},\n",
    "            'n_metal_right': {'resolution': 0.01},\n",
    "            'n': {'resolution': 0.003},\n",
    "            'wg1': {'resolution': 0.002},\n",
    "            'wg2': {'resolution': 0.002},\n",
    "            'p1': {'resolution': 0.003},\n",
    "            'p2': {'resolution': 0.003},\n",
    "        }\n",
    "\n",
    "    def _create_polygons(self):\n",
    "        #We will now set the RF properties of the metals and the BCB\n",
    "        freq = np.linspace(0.1,100, 100) #GHz. This will be the simulation frequency\n",
    "            \n",
    "        eps_rf_metal = 1 - 1j*6e7/(2*np.pi*freq*1e9 * self.e0)\n",
    "\n",
    "        bcb_eps_real = 2.65*np.ones(100)\n",
    "        bcb_eps_imag = bcb_eps_real * tand_fitted_bcb(freq)\n",
    "\n",
    "        bcb_eps = bcb_eps_real - 1j*bcb_eps_imag\n",
    "        \n",
    "        #Now we create the PhotoPolygons\n",
    "        self.substrate = imodulator.SemiconductorPolygon(\n",
    "            shapely.box(\n",
    "                -self.w_window/2,\n",
    "                -self.h_box - self.h_bottom,\n",
    "                self.w_window/2,\n",
    "                -self.h_box\n",
    "            ),\n",
    "            rf_eps = 11.7,\n",
    "            name = 'substrate',\n",
    "            optical_material=3**2,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['substrate'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['substrate'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['substrate'],\n",
    "        )\n",
    "\n",
    "        self.background = imodulator.InsulatorPolygon(\n",
    "            shapely.box(\n",
    "                -self.w_window/2,\n",
    "                -self.h_box - self.h_bottom,\n",
    "                self.w_window/2,\n",
    "                self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2 + self.h_metal + self.h_top\n",
    "            ),\n",
    "            rf_eps = 1,\n",
    "            optical_material=1,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['background'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['background'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['background'],\n",
    "            name = 'background'\n",
    "        )\n",
    "\n",
    "        self.box = imodulator.InsulatorPolygon(\n",
    "            shapely.box(\n",
    "                -self.w_window/2,\n",
    "                -self.h_box,\n",
    "                self.w_window/2,\n",
    "                0\n",
    "            ),\n",
    "            rf_eps = 3.9 - 1j*3.9*0.001,\n",
    "            optical_material=1.44**2,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['box'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['box'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['box'],\n",
    "            name = 'box'\n",
    "        )\n",
    "\n",
    "        n_obp_material = obp.GaInPAs(T=300, As = 0, a = obp.InP.a())\n",
    "        self.n = imodulator.SemiconductorPolygon(\n",
    "            shapely.box(\n",
    "                -self.w_sig_metal/2 - self.metal_sep - self.w_gnd_metal,\n",
    "                0,\n",
    "                self.w_sig_metal/2 + self.metal_sep + self.w_gnd_metal,\n",
    "                self.h_n\n",
    "            ),\n",
    "            rf_eps = n_obp_material.dielectric(T=300),\n",
    "            optical_material=n_obp_material.refractive_index(T=300)**2,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['n'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['n'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['n'],\n",
    "            charge_mesh_settings=self.charge_mesh_settings['n'],\n",
    "            name = 'n',\n",
    "            electro_optic_module=InGaAsPElectroOpticalModel,\n",
    "            electro_optic_module_kwargs={\n",
    "                'y': 0,\n",
    "                'T': 300,\n",
    "                'BF_model': 'vinchant'\n",
    "            },\n",
    "            charge_transport_simulator_kwargs={\n",
    "                'sol_obp_material': n_obp_material,\n",
    "                'sol_Nd': 1e18\n",
    "            }\n",
    "        )\n",
    "\n",
    "        wg1_obp_material = obp.GaInPAs(T=300, As = 0.53, a = obp.InP.a())\n",
    "        self.wg1 = imodulator.SemiconductorPolygon(\n",
    "            shapely.box(\n",
    "                -self.w_wg/2,\n",
    "                self.h_n,\n",
    "                self.w_wg/2,\n",
    "                self.h_n + self.h_wg1\n",
    "            ),\n",
    "            rf_eps = wg1_obp_material.dielectric(T=300),\n",
    "            optical_material=wg1_obp_material.refractive_index(T=300)**2,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['wg1'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['wg1'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['wg1'],\n",
    "            charge_mesh_settings=self.charge_mesh_settings['wg1'],\n",
    "            name = 'wg1',\n",
    "            electro_optic_module=InGaAsPElectroOpticalModel,\n",
    "            electro_optic_module_kwargs={\n",
    "                'y': 0.53,\n",
    "                'T': 300,\n",
    "                'BF_model': 'vinchant'\n",
    "            },\n",
    "            charge_transport_simulator_kwargs={\n",
    "                'sol_obp_material': wg1_obp_material,\n",
    "                'sol_Nd': 1e16\n",
    "            }\n",
    "        )\n",
    "\n",
    "        wg2_obp_material = obp.GaInPAs(T=300, As = 0, a = obp.InP.a())\n",
    "        self.wg2 = imodulator.SemiconductorPolygon(\n",
    "            shapely.box(\n",
    "                -self.w_wg/2,\n",
    "                self.h_n + self.h_wg1,\n",
    "                self.w_wg/2,\n",
    "                self.h_n + self.h_wg1 + self.h_wg2\n",
    "            ),\n",
    "            rf_eps = wg2_obp_material.dielectric(T=300),\n",
    "            optical_material=wg2_obp_material.refractive_index(T=300)**2,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['wg2'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['wg2'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['wg2'],\n",
    "            charge_mesh_settings=self.charge_mesh_settings['wg2'],\n",
    "            name = 'wg2',\n",
    "            electro_optic_module=InGaAsPElectroOpticalModel,\n",
    "            electro_optic_module_kwargs={\n",
    "                'y': 0,\n",
    "                'T': 300,\n",
    "                'BF_model': 'vinchant'\n",
    "            },\n",
    "            charge_transport_simulator_kwargs={\n",
    "                'sol_obp_material': wg2_obp_material,\n",
    "                'sol_Nd': 1e16\n",
    "            }\n",
    "        )\n",
    "\n",
    "        p1_obp_material = obp.GaInPAs(T=300, As = 0, a = obp.InP.a())\n",
    "        self.p1 = imodulator.SemiconductorPolygon(\n",
    "            shapely.box(\n",
    "                -self.w_wg/2,\n",
    "                self.h_n + self.h_wg1 + self.h_wg2,\n",
    "                self.w_wg/2,\n",
    "                self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1\n",
    "            ),\n",
    "            rf_eps = p1_obp_material.dielectric(T=300),\n",
    "            optical_material=p1_obp_material.refractive_index(T=300)**2,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['p1'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['p1'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['p1'],\n",
    "            charge_mesh_settings=self.charge_mesh_settings['p1'],\n",
    "            name = 'p1',\n",
    "            electro_optic_module=InGaAsPElectroOpticalModel,\n",
    "            electro_optic_module_kwargs={\n",
    "                'y': 0,\n",
    "                'T': 300,\n",
    "                'BF_model': 'vinchant'\n",
    "            },\n",
    "            charge_transport_simulator_kwargs={\n",
    "                'sol_obp_material': wg2_obp_material,\n",
    "                'sol_Na': 1e17\n",
    "            }\n",
    "        )\n",
    "\n",
    "        p2_obp_material = obp.GaInAs(T=300, a = obp.InP.a())\n",
    "        self.p2 = imodulator.SemiconductorPolygon(\n",
    "            shapely.box(\n",
    "                -self.w_wg/2,\n",
    "                self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1,\n",
    "                self.w_wg/2,\n",
    "                self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2\n",
    "            ),\n",
    "            rf_eps = p2_obp_material.dielectric(T=300),\n",
    "            optical_material=p2_obp_material.refractive_index(T=300)**2,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['p2'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['p2'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['p2'],\n",
    "            charge_mesh_settings=self.charge_mesh_settings['p2'],\n",
    "            name = 'p2',\n",
    "            electro_optic_module=InGaAsPElectroOpticalModel,\n",
    "            electro_optic_module_kwargs={\n",
    "                'y': 0,\n",
    "                'T': 300,\n",
    "                'BF_model': 'vinchant'\n",
    "            },\n",
    "            charge_transport_simulator_kwargs={\n",
    "                'sol_obp_material': p2_obp_material,\n",
    "                'sol_Na': 1e19\n",
    "            }\n",
    "        )\n",
    "\n",
    "        self.bcb_far_left = imodulator.InsulatorPolygon(\n",
    "            shapely.box(\n",
    "                -self.w_window/2,\n",
    "                0,\n",
    "                -self.w_sig_metal/2 - self.metal_sep - self.w_gnd_metal,\n",
    "                self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2\n",
    "            ),\n",
    "            rf_eps = bcb_eps,\n",
    "            optical_material=1.56**2,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['bcb'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['bcb'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['bcb'],\n",
    "            name = 'bcb_far_left'\n",
    "        )\n",
    "\n",
    "        self.bcb_far_right = imodulator.InsulatorPolygon(\n",
    "            shapely.box(\n",
    "                self.metal_sep + self.w_gnd_metal + self.w_sig_metal/2,\n",
    "                0,\n",
    "                self.w_window/2,\n",
    "                self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2\n",
    "            ),\n",
    "            rf_eps = bcb_eps,\n",
    "            optical_material=1.56**2,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['bcb'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['bcb'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['bcb'],\n",
    "            name = 'bcb_far_right'\n",
    "        )\n",
    "\n",
    "        self.bcb_left = imodulator.InsulatorPolygon(\n",
    "            shapely.box(\n",
    "                -self.w_sig_metal/2 - self.metal_sep,\n",
    "                self.h_n,\n",
    "                -self.w_wg/2,\n",
    "                self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2\n",
    "            ),\n",
    "            rf_eps = bcb_eps,\n",
    "            optical_material=1.56**2,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['bcb'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['bcb'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['bcb'],\n",
    "            name = 'bcb_left'\n",
    "        )\n",
    "\n",
    "        self.bcb_right = imodulator.InsulatorPolygon(\n",
    "            shapely.box(\n",
    "                self.w_wg/2,\n",
    "                self.h_n,\n",
    "                self.w_sig_metal/2 + self.metal_sep,\n",
    "                self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2\n",
    "            ),\n",
    "            rf_eps = bcb_eps,\n",
    "            optical_material=1.56**2,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['bcb'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['bcb'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['bcb'],\n",
    "            name = 'bcb_right'\n",
    "        )\n",
    "\n",
    "\n",
    "        self.sig_metal = imodulator.MetalPolygon(\n",
    "            shapely.Polygon(\n",
    "                [\n",
    "                    (-self.base_tree/2,0+self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2),\n",
    "                    (-self.base_tree/2, self.base_tree+self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2),\n",
    "                    (-self.width_tree/2, self.base_tree+self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2),\n",
    "                    (0, self.height_tree + self.base_tree+self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2),\n",
    "                    (self.width_tree/2, self.base_tree+self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2),\n",
    "                    (self.base_tree/2, self.base_tree+self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2),\n",
    "                    (self.base_tree/2,0+self.h_n + self.h_wg1 + self.h_wg2 + self.h_p1 + self.h_p2)\n",
    "                ]\n",
    "            ),\n",
    "            rf_eps = eps_rf_metal,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['sig_metal'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['sig_metal'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['sig_metal'],\n",
    "            name = 'sig_metal',\n",
    "            calculate_current=True,\n",
    "            d_buffer_current=min(self.w_sig_metal/20, self.h_metal/20, 0.05)\n",
    "        )\n",
    "\n",
    "        self.n_metal_left = imodulator.MetalPolygon(\n",
    "            shapely.box(\n",
    "                -self.w_sig_metal/2 - self.metal_sep - self.w_gnd_metal,\n",
    "                self.h_n,\n",
    "                -self.w_sig_metal/2 - self.metal_sep,\n",
    "                self.h_n + self.h_metal\n",
    "            ),\n",
    "            rf_eps = eps_rf_metal,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['n_metal_left'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['n_metal_left'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['n_metal_left'],\n",
    "            name = 'n_metal_left',\n",
    "            calculate_current=False\n",
    "        )\n",
    "\n",
    "        self.n_metal_right = imodulator.MetalPolygon(\n",
    "            shapely.box(\n",
    "                self.w_sig_metal/2 + self.metal_sep,\n",
    "                self.h_n,\n",
    "                self.w_sig_metal/2 + self.metal_sep + self.w_gnd_metal,\n",
    "                self.h_n + self.h_metal\n",
    "            ),\n",
    "            rf_eps = eps_rf_metal,\n",
    "            eo_mesh_settings=self.eo_mesh_settings['n_metal_right'],\n",
    "            rf_mesh_settings=self.rf_mesh_settings['n_metal_right'],\n",
    "            optical_mesh_settings=self.optical_mesh_settings['n_metal_right'],\n",
    "            name = 'n_metal_right',\n",
    "            calculate_current=False\n",
    "        )\n",
    "\n",
    "\n",
    "    def _initialize_device(self):\n",
    "        photo_polygons = [\n",
    "            self.sig_metal,\n",
    "            self.n_metal_left,\n",
    "            self.n_metal_right,\n",
    "            self.p2,\n",
    "            self.p1,\n",
    "            self.wg2,\n",
    "            self.wg1,\n",
    "            self.n,\n",
    "            self.box,\n",
    "            self.bcb_left,\n",
    "            self.bcb_right,\n",
    "            self.bcb_far_left,\n",
    "            self.bcb_far_right,\n",
    "            self.substrate,\n",
    "            self.background\n",
    "        ]\n",
    "        \n",
    "        #Just in case there are empty polygons\n",
    "        idxs_to_remove = []\n",
    "        for i, poly in enumerate(photo_polygons):\n",
    "            if np.isclose(poly.polygon.bounds[1], poly.polygon.bounds[3]):\n",
    "                idxs_to_remove.append(i)\n",
    "        for i in idxs_to_remove[::-1]:\n",
    "            del photo_polygons[i]\n",
    "        self.device = imodulator.PhotonicDevice(\n",
    "            photo_polygons\n",
    "        )\n",
    "\n",
    "\n",
    "eopm = InP_EOPM_tree(\n",
    "    height_tree = 4,\n",
    "    base_tree = 2,\n",
    "    width_tree = 6\n",
    ")\n",
    "eopm._make_meshes()\n",
    "eopm._create_polygons()\n",
    "eopm._initialize_device()\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "gs = fig.add_gridspec(1,2)\n",
    "\n",
    "ax1 = fig.add_subplot(gs[0,0])\n",
    "ax2 = fig.add_subplot(gs[0,1])\n",
    "\n",
    "for ax in [ax1, ax2]:\n",
    "    eopm.device.plot_polygons(\n",
    "        color_polygon=\"black\",\n",
    "        color_line=\"green\",\n",
    "        color_junctions=\"blue\",\n",
    "        fill_polygons=True,\n",
    "        fig=fig,\n",
    "        ax=ax,\n",
    "    )\n",
    "\n",
    "ax2.set_xlim(-5,5)\n",
    "ax2.set_ylim(-5,10)\n",
    "\n",
    "ax1.set_xlabel('x (um)')\n",
    "ax1.set_ylabel('y (um)')\n",
    "ax2.set_xlabel('x (um)')\n",
    "ax2.set_ylabel('y (um)')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "imodulator_venv",
   "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
