01-simple_geometries.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defining geometries for MEEP/MPB\n",
    "\n",
    "Simple shapes can be directly defined in PyMEEP. Later, we will see how to import arbitrary shapes made with more specialized tools.\n",
    "For visualization, we will call the Python packages [Matplotlib](https://matplotlib.org/) and [Ipyvolume](https://ipyvolume.readthedocs.io/en/latest/).\n",
    "\n",
    "References :\n",
    "\n",
    "1. https://mpb.readthedocs.io/en/latest/Python_Tutorial/\n",
    "2. https://mpb.readthedocs.io/en/latest/Python_User_Interface/\n",
    "2. http://www.simpetuscloud.com/projects.html#mpb_waveguide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import meep and mpb (from meep)\n",
    "import meep as mp\n",
    "from meep import mpb\n",
    "\n",
    "# arrays\n",
    "import numpy as np\n",
    "\n",
    "# plotting\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cross-section of a strip waveguide\n",
    "\n",
    "As a (useful) example, let's characterize a strip waveguide (Figure from [3])\n",
    "\n",
    "![alt text](http://www.simpetus.com/files/SOI_waveguide_bands.png)\n",
    "\n",
    "The left image is the geometry we want to simulate. Since we have translational invariance along $x$, the propagation direction, we will only consider the $yz$ cross-section. The right image contains the information we want : dispersion relation (allowed modes), and field profiles.\n",
    "\n",
    "## 1. Computational domain\n",
    "\n",
    "Since the computer has finite memory, our geometry will live inside a \"super\" geometry. For MPB, this \"supercell\" always has periodic boundary conditions. But if it is big enough, we don't expect it to have an effect on localized modes, which is what we are trying to solve for a waveguide.\n",
    "\n",
    "### 1.1 A note on geometry\n",
    "\n",
    "By default, everything is 3D. For a waveguide, we can use the translational invariance along the propagation direction (say, $x$) to do a 2D simulation. Then, the only thing that matters is the geometries' intersection with the yz plane (x=0). We can set this when we define the size of the computational cell, controlled by the geometry_lattice argument, an object of the meep.Lattice class: we can set some of the dimensions to have a size 0, which reduces the dimensionality of the system. [1]\n",
    "\n",
    "### 1.2 A note on units\n",
    "\n",
    "Since Maxwell's equations with linear permittivity are [scale-invariant](https://mpb.readthedocs.io/en/latest/Python_Tutorial/#a-few-words-on-units), a simulation in MPB with units of $x$ outputs frequencies equivalent to $f=x/\\lambda$. By default, dimensions are in microns. For instance, to look for the wavelength response at $1.5\\mu$m, we look at frequency $1/1.5 = 0.6667$ (the physical frequency outside of MPB would be $f=c/\\lambda$).\n",
    "\n",
    "Let's define the supercell :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc_y = 2  # supercell width (um)\n",
    "sc_z = 2  # supercell height (um)\n",
    "resolution = 64  # pixels/um\n",
    "geometry_lattice = mp.Lattice(size=mp.Vector3(0, sc_y, sc_z)) # Computational domain object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we used variables to parametrize the supercell. The resolution should be a multiple of 2 for better results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters for the waveguide\n",
    "w = 0.5  # Si width (um)\n",
    "h = 0.25  # Si height (um)\n",
    "\n",
    "# Materials\n",
    "Si = mp.Medium(index=3.45)\n",
    "SiO2 = mp.Medium(index=1.45)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Geometry\n",
    "\n",
    "The geometry that will be passed to the solver is a list of meep objects. Objects typically have a shape and index. Later objects in the list will be overlaid \"on top\" of earlier objects, effectively replacing them (useful to draw holes). See https://meep.readthedocs.io/en/latest/Python_User_Interface/#geometricobject for common geometries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the 2d blocks for the strip and substrate\n",
    "geometry = [mp.Block(size=mp.Vector3(mp.inf, mp.inf, 0.5 * (sc_z - h)),\n",
    "                     center=mp.Vector3(z=0.25 * (sc_z + h)), material=SiO2),\n",
    "            mp.Block(size=mp.Vector3(mp.inf, w, h), material=Si)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 - Visualizing the geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While we could directly convert the geometry object above into a visualization, MPB and MEEP offer an easy way to display the dielectric distribution once the simulation has been initialized. The function we are interested in is `ModeSolver.init_params(p, reset_fields)` [2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Working in 3 dimensions.\n",
      "Grid size is 1 x 128 x 128.\n",
      "Solving for 1 bands at a time.\n",
      "Creating Maxwell data...\n",
      "Mesh size is 3.\n",
      "Lattice vectors:\n",
      "     (1, 0, 0)\n",
      "     (0, 2, 0)\n",
      "     (0, 0, 2)\n",
      "Cell volume = 4\n",
      "Reciprocal lattice vectors (/ 2 pi):\n",
      "     (1, -0, 0)\n",
      "     (-0, 0.5, -0)\n",
      "     (0, -0, 0.5)\n",
      "Geometric objects:\n",
      "     block, center = (0,0,0.5625)\n",
      "          size (1e+20,1e+20,0.875)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "     block, center = (0,0,0)\n",
      "          size (1e+20,0.5,0.25)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "Geometric object tree has depth 2 and 8 object nodes (vs. 2 actual objects)\n",
      "Initializing epsilon function...\n",
      "Allocating fields...\n",
      "Solving for band polarization: .\n",
      "Initializing fields to random numbers...\n"
     ]
    }
   ],
   "source": [
    "ms = mpb.ModeSolver(\n",
    "    geometry_lattice=geometry_lattice,\n",
    "    geometry=geometry,\n",
    "    resolution=resolution)\n",
    "ms.init_params(mp.NO_PARITY, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's ask the modesolver for the permittivity distribution it will simulate and plot it :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epsilon: 1-11.9025, mean 1.82332, harm. mean 1.3534, 47.7539% > 1, 7.55167% \"fill\"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAEKCAYAAABewe3GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAYh0lEQVR4nO3dfYxc1X3G8e+zjh0KRQViXsxbcCorKm6KQ1cGSlWgwZGx0jqkqWTUEkJTWSSxSqKqqtVIpC+qRNOmlSgEd5M4GCkFtSqEbbrEOFYKoSmJX2r8EqB2HCpcr2xjHFxKBNnk1z/mLJ2MZ3bveu+ce+/k+UhXc+895945dwZ+PnvmvCgiMDOz/huqugBmZj8pHHDNzDJxwDUzy8QB18wsEwdcM7NMHHDNzDKpNOBKWi/psKTdPdIl6S5J+yTtlHR5W9pySc+ltLX5Sm1mdnKqruHeByyfIv0GYFHaVgP3AkiaA9yT0i8FbpJ0aV9LamY2S5UG3Ih4Anhpiiwrgfuj5SngDEkLgKXAvojYHxGvAw+mvGZmtfWmqgswjQuAF9qOD6Rz3c5f0e0GklbTqh1z6qmn/uLb3vY2hoaqrtibDaYXXniBl156SbO5h6SZDH/dGBFT/ZVcK3UPuN2+uJji/IknI0aAEYB3vOMdMTo6yimnnFJeCc3sDcuXZ49983O/4WzUPeAeAC5qO74QOAjM63F+Sm9+85tZuHBhqQU0s/83d+7cUu4jFaskN20umLr/bT0KfCD1VrgSeDkixoEtwCJJCyXNA1alvGY2AIaGhgptTVNpDVfSA8C1wHxJB4BPAnMBImIdMAasAPYBrwK3prQJSWuAjcAcYH1E7Mn+AGbWF0VruE1TacCNiJumSQ/goz3SxmgFZDMbIJIccM3McnHANTPLxAHXzCwTB1wzswwkNbIHQhEOuGZWO67hmpll4oBrZpaJA66ZWSYOuGZmGfhHMzOzjFzDNTPLxAHXzCwTB1wzsww8eY2ZWUYOuGZmmbiXgplZJq7hmpll4DZcM7OMBjXgVtpQImm5pOck7ZO0tkv6H0jakbbdkn4o6ayU9rykXSlta/7Sm1m/TNZyp9uaprIarqQ5wD3AMlrLoW+RNBoR357MExF/Cfxlyv9rwMcj4qW221wXES9mLLaZZTCoP5pV+VRLgX0RsT8iXgceBFZOkf8m4IEsJTOzyhSt3TaxhltlwL0AeKHt+EA6dwJJpwLLgX9qOx3AY5K2SVrdt1KaWXaDGnCr/NGs26cVPfL+GvBvHc0JV0fEQUnnAJskPRsRT5zwJq1gvBrg4osvnm2ZzSyDJgbTIqqs4R4ALmo7vhA42CPvKjqaEyLiYHo9DDxMq4niBBExEhHDETF89tlnz7rQZtZ/g1rDrTLgbgEWSVooaR6toDramUnSzwDXAI+0nTtN0umT+8C7gd1ZSm1mfTeoAbeyJoWImJC0BtgIzAHWR8QeSbel9HUp643AYxHxv22Xnws8nD7wNwF/HxFfyVd6M+sXT0DeJxExBox1nFvXcXwfcF/Huf3AZX0unplVpIm11yI80szMascB18wsEwdcM7MMmvqDWBEOuGZWO4MacAfzp0Aza7ShoaFC23QknSLpW5KelrRH0p90ySNJd6VJtHZKurwvD4VruGZWQyXWcF8DfjUiXpE0F3hS0qMR8VRbnhuARWm7Arg3vZbOAdfMaqXMNtyICOCVdDg3bZ1TCKwE7k95n5J0hqQFETFeSiHauEnBzGpnBiPN5kva2radMJGVpDmSdgCHgU0R8c2OLIUn0pot13DNrHZmUMN9MSKGp8oQET8Elkg6g9YI1Z+PiPapAGYykdasOOCaWe30Y2hvRHxP0r/Smuq1PeDOZCKtWXGTgpnVSpkTkEs6O9VskfRTwPXAsx3ZRoEPpN4KVwIv96P9FlzDNbMaKrGXwgJgg1pLeg0B/xARX+6YJGsMWAHsA14Fbi3rzTs54JpZ7ZTYS2En8M4u59e17Qfw0VLecBoOuGZWO4M60swB18xqxwHXzCwDT0BuZpaRa7hmZpk44JqZZTKoAbfShhJJyyU9l6ZFW9sl/VpJL0vakbY7il5rZs1U5sCHuqmshps6It8DLKM1tG6LpNGI+HZH1q9HxHtO8loza6AmBtMiqqzhLgX2RcT+iHgdeJDWNGn9vtbMaq6sCcjrpsoSF50S7ao0W/ujkhbP8FokrZ6cuu3IkSNllNvM+mxQmxSqDLhFpkTbDrw1Ii4D/hb40gyubZ2MGImI4YgYPvvss0+6sGaWxyC34VYZcKedEi0ijkfEK2l/DJgraX6Ra82suRxwy7cFWCRpoaR5wCpa06S9QdJ5Sp+qpKW0ynu0yLVm1lyDGnAr66UQEROS1gAbgTnA+ojY0zFt2vuBD0uaAL4PrEoz+3S9tpIHMbPSNfEHsSIqHfiQmgnGOs61T5t2N3B30WvNrPmaWnstwiPNzKx2HHDNzDJxwDUzy8QB18wsEwdcM7MMPAG5mVlGruGamWXigGtmlokDrplZBh74YGaWkQOumVkm7qVgZpaBmxTMzDJywDUzy8QB18wsEwdcM7MMPLTXzCwj13DNzDIZ1IBbab1d0nJJz0naJ2ltl/TfkrQzbd+QdFlb2vOSdknaIWlr3pKbWT95EcmSSZoD3AMso7Xs+RZJoxHx7bZs3wWuiYhjkm4ARoAr2tKvi4gXsxXazLJoYjAtosomhaXAvojYDyDpQWAl8EbAjYhvtOV/CrgwawnNLLum1l6LqLJJ4QLghbbjA+lcLx8CHm07DuAxSdskre51kaTVkrZK2nrkyJFZFdjM8hgaGiq0TUfSRZK+JukZSXsk3d4lz7WSXk7Nkzsk3dGXh6LaGm63f8Kia0bpOloB95fbTl8dEQclnQNskvRsRDxxwg0jRmg1RTA8PNz1/mZWLyXWcCeA34+I7ZJOB7ZJ2tTRdAnw9Yh4T1lv2kuVNdwDwEVtxxcCBzszSfoF4HPAyog4Onk+Ig6m18PAw7SaKMxsAJT1o1lEjEfE9rT/P8AzTP2XdF9VGXC3AIskLZQ0D1gFjLZnkHQx8BBwc0T8Z9v509K/Vkg6DXg3sDtbyc2sb4oG2xRw5082GaZtqubFS4B3At/sknyVpKclPSppcV8ejAqbFCJiQtIaYCMwB1gfEXsk3ZbS1wF3AG8BPpM+3ImIGAbOBR5O594E/H1EfKWCxzCzPphBk8KLKSZMd7+fBv4J+FhEHO9I3g68NSJekbQC+BKwaCblLarSgQ8RMQaMdZxb17b/u8DvdrluP3BZ53kzGwxlDu2VNJdWsP1iRDzUmd4egCNiTNJnJM3vR5dTjzQzs9op60cztW70eeCZiPjrHnnOAw5FREhaSqup9Wi3vLPlgGtmtVJyP9yrgZuBXZJ2pHN/BFwMb/xF/X7gw5ImgO8DqyKiLz2aHHDNrHbKCrgR8STdu6C257kbuLuUN5yGA66Z1c6gjjRzwDWz2nHANTPLwBOQm5ll5BquWfKFL3yBxx9/vOpilO6aa64B4NZbb624JOaAa2aWiQOuWfL444+zYcOGqovRN67hVs8B18wsg0GegNwB18xqx70UzMwycQ3XzCwTB1wzswzchmtmlpEDrplZJv7RzMwsk0Gt4Vb6z4ik5ZKek7RP0tou6ZJ0V0rfKenyoteaWTPNcBHJRqks4EqaA9wD3ABcCtwk6dKObDfQWsxtEbAauHcG15pZQw1qwC3UpCBpM/DptOjj5LmRiOi5JHEBS4F9aUFIJD0IrAS+3ZZnJXB/Wu7iKUlnSFoAXFLgWjNrqDoHU0lvBn6DVhx6I4ZGxJ9Od23RGu5C4A8lfbLt3LRLE0/jAuCFtuMD6VyRPEWuBUDS6sk1648cOTLLIptZDjWv4T5Cq4I3Afxv2zatoj+afQ94F3CXpH8GfvskCtmp26fVuXBbrzxFrm2djBgBRgCGh4f7sjCcmZWnAROQXxgRy0/mwqJPpYiYiIiP0Frf/UngnJN5wzYHgIvaji8EDhbMU+RaM2uomtdwvyHpHSdzYdGAu25yJyLuAz4IPHYyb9hmC7BI0kJJ84BVwGhHnlHgA6m3wpXAyxExXvBaM2uomgfcXwa2pV5SOyXtkrSzyIWFmhQi4u86jrcBvzPzcv7YPSYkrQE2AnOA9RGxR9JtKX0dMAasAPYBrwK3TnXtbMpjZvVR5x/NaPWOOimVDnxIvR7GOs6116YD+GjRa81sMNQ54EbEf53stR5pZma10tQ+tkU44JpZ7dS8l8JJc8A1s9pxDdfMLBMHXDOzDNyGa2aWkQOumVkm/tHMzCwT13AHwGuvvcZ3v/tdTjnllKqL0mhDQ0O8733vq7oYpZusVY2Pj1dckub6wQ9+MOt7uA13QPzoRz/i+PHjHD9+vOqiNNrtt99edRH66vDhw1UXobEmJiZKuY8DrplZJoMacAezZdrMGq2s2cIkXSTpa5KekbRH0gl/nqXZCLuunVg213DNrFZKnoB8Avj9iNgu6XRa0ypuioj25bja1068gtbaiVeUVYB2ruGaWe2UVcONiPGI2J72/wd4hhOX43pj7cSIeAqYXDuxdK7hmlntzKANd76krW3HI2lZrW73vAR4J/DNjqReaySW3l3FAdfMamcGAffFiJh2QVtJP01rebCPRURnN6XCayTOlgOumdVOmb0UJM2lFWy/GBEPdcmSbY1Et+GaWa0Ubb8t2EtBwOeBZyLir3tk67V2YukqCbiSzpK0SdLe9Hpmlzw9u3NI+mNJ/y1pR9pW5H0CM+unoaGhQlsBVwM3A7/aHi8k3Ta5fiKtpbr201o78bPAR/ryUFTXpLAW2BwRd0pam47/sCPPdN05/iYi/ipjmc0sk7KaFCLiSbq30bbn6bl2YtmqalJYCWxI+xuA93ZmKNidw8wGUM2XST9pVQXccyfbSNLrOVNl7tGdY00aFbK+W5NE27WrJW2VtPXYsWOzL7mZ9VWZbbh107eAK+mrknZ32VbO8D7dunPcC/wssIRWX7lP97o+IkYiYjgihs88s2dcNrMaGdSA27c23Ii4vleapEOSFkTEeBrR0XV6pl7dOSLiUFuezwJfLq/kZla1QZ2AvKqnGgVuSfu3AI90ZpiqO0fHsLsbgd19KqeZZeYmhfLdCSyTtBdYlo6RdL6ksZSna3eOlPYpSbsk7QSuAz6eufxm1keDGnAr6RYWEUeBd3U5fxBYkfZ7dueIiJv7WkAzq1QTg2kRHtprZrXjgGtmlokDrplZBip3AvJaccA1s9pxDdfMLBMHXDOzTBxwzcwyaGof2yIccM2sdvyjmZlZJq7hmpll4oBrZpaB23DNzDJywDUzy8QB18wsE/dSMDPLwG24ZmYZOeCamWXigGtmlsmgBtxKWqYlnSVpk6S96bXr+uWSnk9rl+2QtHWm15tZMw3qmmZV/RS4FtgcEYuAzem4l+siYklEDJ/k9WbWIJMTkBfZmqaqEq8ENqT9DcB7M19vZjXmGm65zo2IcYD0ek6PfAE8JmmbpNUncT2SVkvaKmnrsWPHSiq+mfXToAbcvv1oJumrwHldkj4xg9tcHREHJZ0DbJL0bEQ8MZNyRMQIMAKwePHimMm1ZlaNJgbTIvoWcCPi+l5pkg5JWhAR45IWAId73ONgej0s6WFgKfAEUOh6M2ueptZei6iqSWEUuCXt3wI80plB0mmSTp/cB94N7C56vZk1l380K9edwDJJe4Fl6RhJ50saS3nOBZ6U9DTwLeBfIuIrU11vZoPBbbglioijwLu6nD8IrEj7+4HLZnK9mQ2GJgbTIjzSzMxqpam11yKa1whiZgOvrCYFSeslHZa0u0f6tZJeTqNZd0i6o/SHaeMarpnVTok13PuAu4H7p8jz9Yh4T1lvOBUHXDOrnbJ6IETEE5IuKeVmJXCTgpnVStHmhFQLnj85kjRtq6e7fxdXSXpa0qOSFpf8OD/GNVwzq50ZNCm82DGx1UxtB94aEa9IWgF8CVg0i/tNyTVcM6udXP1wI+J4RLyS9seAuZLmz/rGPbiGa2a1k6tbmKTzgEMREZKW0qqEHu3X+zngmlntlBVwJT0AXEurrfcA8ElgLkBErAPeD3xY0gTwfWBVRPRtkisHXDOrlckJyMsQETdNk343rW5jWTjgmlntDOpIMwdcM6sdB1wzs0wccM3MMhjkyWsccM2sdpo4uXgRDrhmVjuu4ZqZZeKAa2aWwSC34VbSUCLpLEmbJO1Nr2d2yfP2tkmBd0g6LuljKe2PJf13W9qK/E9hZv0yqGuaVdUyvRbYHBGLgM3p+MdExHMRsSQilgC/CLwKPNyW5W8m09OkE2Y2IBxwy7US2JD2NwDvnSb/u4DvRMR/9bVUZlYLXia9XOdGxDhAej1nmvyrgAc6zq2RtDOtWXRCk4SZNdMMJyBvlL4FXElflbS7y7ZyhveZB/w68I9tp+8FfhZYAowDn57i+tWTs8EfO3bsJJ7EzHIb1IDbt14KEXF9rzRJhyQtiIhxSQuAw1Pc6gZge0Qcarv3G/uSPgt8eYpyjAAjAIsXL+7btGtmVp4mBtMiqmpSGAVuSfu3AI9MkfcmOpoTUpCedCPQdQlkM2umQa3hVhVw7wSWSdoLLEvHSDpf0hs9DiSdmtIf6rj+U5J2SdoJXAd8PE+xzSyHQQ24lQx8iIijtHoedJ4/CKxoO34VeEuXfDf3tYBmVpkyJyCvG480M7PaaWLttQgHXDOrHQdcM7NMHHDNzDJo6g9iRTjgmlnt+EczM7NMXMM1M8vEAdfMLAO34ZqZZeSAa2aWiQOumVkGHtprZpaRa7hmZpk44JqZZeKAa2aWiQOumVkG7odrZpaReymYmWXiGq6ZWSaDGnArqbdL+k1JeyT9SNLwFPmWS3pO0j5Ja9vOnyVpk6S96fXMPCU3s34ruoBkkaAsab2kw5K6ruytlrtSjNkp6fLSH6hNVQ0lu4H3AU/0yiBpDnAPcANwKXCTpEtT8lpgc0QsAjanYzMbECWu2nsfsHyK9BuARWlbDdw768JPoZKAGxHPRMRz02RbCuyLiP0R8TrwILAypa0ENqT9DcB7+1NSM6vC0NBQoW06EfEE8NIUWVYC90fLU8AZkhaU9BgnqHMb7gXAC23HB4Ar0v65ETEOEBHjks7pdRNJq2n9ywXw2pIlS7r+adFw84EXqy5Enwzqsw3qc719tjfYtm3bRknzC2Y/RdLWtuORiBiZwdt1izMXAOMzuEdhfQu4kr4KnNcl6RMR8UiRW3Q5FzMtR/rwR1KZtkZEzzbjphrU54LBfbZBfq7Z3iMipmoCKFspcaaovgXciLh+lrc4AFzUdnwhcDDtH5K0INVuFwCHZ/leZvaTaao4U7o69y7eAiyStFDSPGAVMJrSRoFb0v4tQJEas5lZp1HgA6m3wpXAy5PNlf1QVbewGyUdAK4C/kXSxnT+fEljABExAawBNgLPAP8QEXvSLe4ElknaCyxLx0XMpG2nSQb1uWBwn83PlYGkB4B/B94u6YCkD0m6TdJtKcsYsB/YB3wW+EhfyxPRt+YKMzNrU+cmBTOzgeKAa2aWyUAH3NkOIa6rokObJT0vaZekHWV01+mX6T7/3MMvy1Tg2a6V9HL6jnZIuqOKcs5U3YbMNkZEDOwG/Bytjtj/Cgz3yDMH+A7wNmAe8DRwadVln+a5PgWsTftrgb/oke95YH7V5Z3mWab9/IEVwKO0+kxeCXyz6nKX+GzXAl+uuqwn8Wy/AlwO7O6R3sjvrN/bQNdwY/ZDiOtqkIY2F/n8sw6/LFET/9sqJGo2ZLYpBjrgFtRraF+d/djQZqDX0OYAHpO0LQ1xrqMin38TvyMoXu6rJD0t6VFJi/MUre+a+p31VZ3nUiikLkOIyzbVc83gNldHxME018QmSc+mmkmdFPn8a/kdFVCk3NuBt0bEK5JWAF+iNXNV0zX1O+urxgfc6O8Q4spM9VySCg1tjoiD6fWwpIdp/Ylbt4Bb5POv5XdUwLTljojjbftjkj4jaX5ENH1im6Z+Z33lJoWphxDX1bRDmyWdJun0yX3g3bTmIa6bIp9/1uGXJZr22SSdpzSxq6SltP6fPJq9pOVr6nfWX1X/atfPDbiR1r+0rwGHgI3p/PnAWFu+FcB/0vpF+RNVl7vAc72F1sTre9PrWZ3PReuX8afTtqfOz9Xt8wduA25L+6I1Gf13gF306HFSx63As61J38/TwFPAL1Vd5oLP9QCtKQx/kP4f+9CgfGf93Dy018wsEzcpmJll4oBrZpaJA66ZWSYOuGZmmTjgmpll4oBrZpaJA66ZWSYOuFY7kv5M0u1tx38u6feqLJNZGTzwwWpH0iXAQxFxuaQhWiPqlkbEIAx5tZ9gjZ+8xgZPRDwv6aikdwLnAv/hYGuDwAHX6upzwAdpTVG5vtqimJXDTQpWS2l2rV3AXGBRRPyw4iKZzZpruFZLEfG6pK8B33OwtUHhgGu1lH4suxL4zarLYlYWdwuz2pF0KbAP2BwRe6suj1lZ3IZrZpaJa7hmZpk44JqZZeKAa2aWiQOumVkmDrhmZpn8H25tnLQ0vt/4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "\n",
    "n = np.sqrt(ms.get_epsilon())\n",
    "\n",
    "pos = ax.imshow(n.T, cmap='gray_r', interpolation='none', extent=[-sc_y/2,sc_y/2,-sc_z/2,sc_z/2] )\n",
    "cbar = fig.colorbar(pos, ax=ax)\n",
    "cbar.set_label('n')\n",
    "ax.set_xlabel('y')\n",
    "ax.set_ylabel('z')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.5D photonic crystal\n",
    "\n",
    "It is just as easy to define different elements and to extend to 3D. The periodic boundary conditions of MPB make a cell containing a single element automatically a photonic crystal. To simulate defect modes, one can proceed as above i.e. have a larger computational domain with the bounded mode decaying before reaching the boundary.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Working in 3 dimensions.\n",
      "Grid size is 16 x 16 x 64.\n",
      "Solving for 1 bands at a time.\n",
      "Creating Maxwell data...\n",
      "Mesh size is 3.\n",
      "Lattice vectors:\n",
      "     (0.5, 0, 0)\n",
      "     (0, 0.5, 0)\n",
      "     (0, 0, 2)\n",
      "Cell volume = 0.5\n",
      "Reciprocal lattice vectors (/ 2 pi):\n",
      "     (2, -0, 0)\n",
      "     (-0, 2, -0)\n",
      "     (0, -0, 0.5)\n",
      "Geometric objects:\n",
      "     cylinder, center = (0,0,0)\n",
      "          radius 0.2, height 0.875, axis (0, 0, 1)\n",
      "Geometric object tree has depth 1 and 1 object nodes (vs. 1 actual objects)\n",
      "Initializing epsilon function...\n",
      "Allocating fields...\n",
      "Solving for band polarization: .\n",
      "Initializing fields to random numbers...\n"
     ]
    }
   ],
   "source": [
    "# Supercell\n",
    "sc_x = 0.5 # um\n",
    "sc_y = 0.5 # um\n",
    "sc_z = 2 # um\n",
    "resolution = 32  # pixels/um\n",
    "geometry_lattice = mp.Lattice(size=mp.Vector3(sc_x, sc_y, sc_z)) # Computational domain object\n",
    "\n",
    "# Objects\n",
    "geometry = [mp.Cylinder(radius=0.2, height=0.5 * (sc_z - h), material=Si)] # 200nm radius cylinder\n",
    "\n",
    "# ModeSolver\n",
    "ms = mpb.ModeSolver(\n",
    "    geometry_lattice=geometry_lattice,\n",
    "    geometry=geometry,\n",
    "    resolution=resolution)\n",
    "ms.default_material = SiO2 # Default fill for background\n",
    "ms.init_params(mp.NO_PARITY, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usually (say on your own machine), you can use Mayavi to easily visualize in 3D with `sim.plot3D`. The graphical backend of the cluster doesn't allow for Mayavi, so we can look at cuts :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epsilon: 2.1025-11.9025, mean 4.25806, harm. mean 2.62689, 100% > 1, 21.9955% \"fill\"\n",
      "(16, 16, 64)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAEWCAYAAAAq1S8mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAd2UlEQVR4nO3df5QdZZ3n8fcnIYASNED4ESAxjJtVwAHBXmCWOIqCJ+HIBnfcPUEWGQ9uDsjPXWZWRmbFmfHMeHQPuwMJcDKCgIMwLj8k44YfkcXBiDDpsCEYAhIjLDEZYgg/lUED3/2jns6Wl3u7771d96m+3Z/XOff0vVVPVX27aT795KmqpxQRmJlZ702quwAzs4nCgWtmlokD18wsEweumVkmDlwzs0wcuGZmmThwbdQkHS/pKUmvSjq17nrMxioHrg1L0vWSvjxCsz8HFkfE1Ij4To66hiNphqRlkjZLCkmzR2g/W9L9kn4l6QlJJ+ap1CYaB65V4V3Aum42lLRLxbUAvAncDfxBm+1vBv4PsA9wKXCrpH17UJdNcA7cMUzSuyVtl3R0+nygpG2SPizp30la3dD+YklNe5iS9pb0jdTre2GonaQ/lLSyoW1I+heSFgGnA/8lDRf8fZP9/hT4HeDvU5vdUp3LUu0bJP3HUvsvSbpV0t9Kehn4w4b9HZj2M/T6laSOboeMiOci4ipg1UhtJf1L4Gjgsoh4LSJuAx4jhXWq929L7Wenn88u6fP3JX1Z0oNDPyNJ+0i6SdLLklaN1MO2icOBO4ZFxE+BzwM3SXo78A3g+oj4PrAMOETSoaVN/gPwzRa7+ybwduBwYD/gv7dx/KXATcBX03DBKU3avBv4v8Apqc3rFD3GTcCBwCeBv5T00dJmC4BbgWlp/+X9bU77mRoRU4E7gFsAJM2V9OIwr7kjfU9NHA5sjIhXSsseTcvbtRA4AzgIeDfwI4r/VnsD64HLuqjLxqFe/HPOKhQRfyPpFOBhIIB/k5a/LunvKEL2UkmHA7OB7zbuQ9IMYD6wT0S8kBb/Qy/qlTQTmAt8PCL+GVgj6esUgXRfavaj0ljva8Ps6/PAe9P+iIiVFCFdpanASw3LXqIIz3Z9I/1xRNJdwGER8b30+X8Cf1FFodb/3MPtD38DvA+4MvUgh9wAfEqSKALt2w3rh8wEtpfCtpcOTMcq9xif4bcD7NmRdiJpPnAhcGpEtAzlCrwKvKNh2TuAV5q0beW50vvXmnye2l1pNt44cMc4SVOB/wFcC3xJ0t5D6yLiIeDXwAeBT9F6OOFZYG9JzXqHv6QYahg63gEN6zudTm5zOtaepWWzgJ+3u09J76H4Y/LvI+LZ0vIPNozvNr4+2GGtUJzs+52Geo/k/58E/K2fD9D48zFrmwN37PtrYHVEfBb4X8A1DetvBBYDO9I/ud8iIrYAdwFXSdpL0hRJv59WPwocLun9knYHvtSw+XMUJ8XakgLyQeCvJO0u6QjgLBrGaluR9A7gTuBPG7+fiPhBeXy3yesHpf3sDuyWPu6WPjer9yfAGuCyVO8ngCOA21KTNcDvS5ol6Z3An7T5ozB7CwfuGCZpATAPODst+s/A0ZJOLzX7JsVwQ6ve7ZAzgN8ATwBbgYtgZ+D8OfA94CmgMbSvBQ5LJ6Xavcb2NIrx5M0UJ70ui4gVbW57NPAe4PJy77XNbcteoxgugOJ73jksIekaSeU/XAuBAeAF4CvAJyPiFwCp7r8D1gKraTJGbtYueQLy/ibpbRQBenREPFV3PWbWmnu4/e8cYJXD1mzs82VhfUzS04AAz19g1gc8pGBmlomHFMzMMhmXQwrTp0+P2bNn112G2bj19NNPs23bNo1mHx3OkXFPRMwbzfHGgnEZuLNnz2ZwcLDuMszGrYGBgdyHnJ77gL0wLgPXzPpDcVf6yMbLuSYHrpnVZtKk9k4jvfHGGz2uJA8HrpnVpt0e7njhwDWzWkhy4JqZ5eLANTPLxIFrZpaJA9fMLANJbV+lMF44cM2sNu7hmpll4sA1M8vEgWtmlokD18wsA580MzPLyD1cM7NMHLhmZpk4cM3MMpiIk9fUOmItaZ6kJyVtkHRJk/WnS1qbXg9KOrKOOs2sN4ZCd6TXCPvYXdI/SnpU0jpJf9akTcsskfS0pMckrZHU00fF1NbDlTQZWAKcBGwCVklaFhGPl5r9DPhQRLwgaT6wFDg2f7Vm1gsVXaXwOvCRiHhV0hRgpaS7IuKhUpuRsuSEiNhWRTHDqXNI4RhgQ0RsBJB0C7AA2Bm4EfFgqf1DwMFZKzSznqpiSCGK5++8mj5OSa9oaDMmsqTOIYWDgGdLnzelZa2cBdzVaqWkRZIGJQ3+4he/qKhEM+uVdocTUihPH/r/O70WNexrsqQ1wFZgRUQ8PMyhG7MkgHslrW7cb9Xq7OE2+9PW9Elxkk6g+CHNbbWziFhK8c8EBgYGxscT58zGuQ56uNsiouWjgiPiDeD9kqYBd0h6X0T8uMnxmmXJ8RGxWdJ+wApJT0TEA+1/F+2rs4e7CZhZ+nwwsLmxkaQjgK8DCyLi+Uy1mVkGVZw0K4uIF4HvA/OaHKtplkTE5vR1K3AHxXBnT9QZuKuAOZIOkbQrsBBYVm4gaRZwO3BGRPykhhrNrIcmTZrU1ms4kvZNPVskvQ04EXiioU3TLJG0h6Q9h94DHwPe0jOuSm1DChGxQ9J5wD3AZOC6iFgn6ey0/hrgi8A+wFXpr9yO4f5ZYWb9o8LrcGcAN6QrnyYB346I77aZJftTDEFAkYffioi7qyiqGRUn+MaXgYGBGBzs6eV0fWXlypVdbXfllVd2vE23P/cXX3yxq+2mTZvW8TYDA939zT7//PM73mbu3JanHfrawMAAg4ODo0rL3XbbLWbMmNFW22eeeWb1eOhs+U4zM6vNRLvTzIFrZrVx4JqZZeLANTPLwBOQm5ll5B6umVkmDlwzs0wcuGZmGUzECcgduGZWGweumVkmvkrBzCwT93DNzDLwGK6ZWUYOXBuzzjvvvK62u/rqq7va7s033+xqu5y2b9/e8TYbN27s6li33nprx9ucc845XR1r8eLFXW3Xbxy4ZmaZ+KSZmVkGHsM1M8vIgWtmlslEC9yJNYBiZmNKFU/tlbS7pH+U9KikdZL+rEkbSbpC0gZJayUdXVo3T9KTad0lPfg2d3IP18xqU1EP93XgIxHxqqQpwEpJd0XEQ6U284E56XUscDVwbHrw5BLgJGATsErSsoh4vIrCGjlwzawWVU1AHsWTcF9NH6ekV+PTcRcAN6a2D0maJmkGMBvYEBEbU023pLY9CVwPKZhZbToYUpguabD0WtSwn8mS1gBbgRUR8XDDoQ4Cni193pSWtVreE+7hmlltOhhS2DbcY9Ij4g3g/ZKmAXdIel9E/Lh8qGabDbO8Jxy4Zlabqq9SiIgXJX0fmAeUA3cTMLP0+WBgM7Bri+U94SEFM6tFu8MJbVylsG/q2SLpbcCJwBMNzZYBn05XKxwHvBQRW4BVwBxJh0jaFViY2vaEe7hmVpuKergzgBvSFQeTgG9HxHclnQ0QEdcAy4GTgQ3Ar4DPpHU7JJ0H3ANMBq6LiHVVFNWMA7cm3UxEs2TJkh5UYu3qZjKf3P/N+m3Sm4quUlgLHNVk+TWl9wGc22L75RSB3HMOXDOrzUS708yBa2a18OQ1ZmYZOXDNzDJx4JqZZeIJyM3MMvAYrplZRhMtcGvtz480D6Wk90r6kaTXJf1RHTWaWe9UcadZP6mth9vmPJTbgQuAU2so0cx6bDyFaTvq7OEeQ5qHMiJ+DQzNQ7lTRGyNiFXAb+oo0Mx6yz3cfJrNQ3lsTbWYWWZVTUDeT+oM3ErnoUwTEi8CmDVrVre7MbOMxlPvtR11/nlpNT9lVyJiaUQMRMTAvvvuO+rizKz3PKSQz855KIGfU8xD+aka6+nKypUru9ruqquuqrgSG0+6/f1YuHBhx9vMnTu3q2NVYTyFaTtqC9xW81CW57CUdAAwCLwDeFPSRcBhEfFyXXWbWTXGW++1HbXe+NBsHsqGOSz/iWKowczGIQeumVkmvkrBzCwDDymYmWVUReBKmgncCBwAvAksjYi/bmjzx8Dp6eMuwKHAvhGxXdLTwCvAG8CO4R7HPloOXDOrTUU93B3AxRHxiKQ9gdWSVpSnCYiIrwFfS8c8BfhPEbG9tI8TImJbFcUMx4FrZrWpInDT4863pPevSFpPcSfr4y02OQ24edQH7sLEGrE2szFj6Nbedl7AdEmDpdeiFvucTfEE34dbrH87MA+4rbQ4gHslrW6136q4h2tmtemgh7ttpLFVSVMpgvSiYa7VPwX4YcNwwvERsVnSfsAKSU9ExAPtFtYJ93DNrDZV3doraQpF2N4UEbcP03QhDcMJEbE5fd0K3EExk2FPOHDNrDZVBK6KBtcC6yPi8mHavRP4EHBnadke6UQbkvYAPgb8uIJvrSkPKZhZbSq6SuF44AzgMUlr0rIvALPgt+5e/QRwb0T8srTt/sAdqY5dgG9FxN1VFNWMA9fMalHVjQ8RsZLm0702trseuL5h2UbgyFEX0SYH7ihdccUVXW0X0fXUvzYBdPv70c3vY52zhfnWXjOzTHxrr5lZJg5cM7MMPHmNmVlGDlwzs0x80szMLBP3cM3MMvAYrplZRg5cM7NMHLhmZpk4cM3MMhiagHwiceCaWW3cw7WOrF69uu4SzHbqt99HB66ZWSYOXDOzTBy4ZmYZ+MYHM7OMJtpVChPruzWzMaWih0jOlHS/pPWS1km6sEmbD0t6SdKa9Ppiad08SU9K2iDpkh58mzu5h2tmtaloSGEHcHFEPJKewLta0oqIeLyh3Q8i4uMNx58MLAFOAjYBqyQta7JtJdzDNbNatNu7HSmUI2JLRDyS3r8CrAcOarOMY4ANEbExIn4N3AIsGMW3NSwHrpnVporAbdjfbOAo4OEmq39P0qOS7pJ0eFp2EPBsqc0m2g/rjnlIwcxq08FJs+mSBkufl0bE0nIDSVOB24CLIuLlhu0fAd4VEa9KOhn4DjCH5o9X79kjtR24ZlabDnqv2yJiYJj9TKEI25si4vbG9eUAjojlkq6SNJ2iRzuz1PRgYHO7RXWq1iGFkc4OqnBFWr9W0tF11Glm1atqDFdFg2uB9RFxeYs2B6R2SDqGIvueB1YBcyQdImlXYCGwrMJv87fU1sNt8+zgfIpu/xzgWODq9NXMxoGKrlI4HjgDeEzSmrTsC8AsgIi4BvgkcI6kHcBrwMKICGCHpPOAe4DJwHURsa6KopoZMXBTMTdFxAsVH3vn2cF0nKGzg+XAXQDcmH4wD0maJmlGRGypuBYzq0EVgRsRK2k+FltusxhY3GLdcmB5u8eTtBvwB8BsShkaEX8+0rbt9HAPoOh9PgJcB9yTAnC0mp0dbOy9tjqD+JbAlbQIWAQwa9asCsprzwsvVP13yKx7/fb72Ke39t4JvASsBl7vZMMRAzci/lTSfwU+BnwGWCzp28C1EfHTLood0s7ZwbbPIKYzlksBBgYGenaW0cyq0ccTkB8cEfO62bCt7zb1aP8pvXYAewG3SvpqNwdN2jk7mPUMopnlVfV1uJk8KOl3u9lwxMCVdIGk1cBXgR8CvxsR5wAfoBjH6FY7ZweXAZ9OVyscB7zk8Vuz8aNPA3cuxe3DT6arpx6TtLadDdsZw50O/NuIeKa8MCLelPTxFtuMKCKanh2UdHZafw3FQPbJwAbgVxRDGmY2TozBMG3H/G43bGcM94vDrFvf7YHT9m85O5iCduh9AOeO5hhmNnb1Y+A2dj474TvNzKwWY3S4oKccuGZWmz69SqFrDlwzq417uGZmmThwzcwy8BiumVlGDlwzs0x80szMLBP3cK0je+21V1fb9dusTtYfuv19rIPHcM3MMnLgmpll4sA1M8vEgWtmlkEfT0DetYn13ZrZmFLRU3tnSrpf0npJ6yRd2KTN6Wnu2rWSHpR0ZGnd02lO2zWSBnvwbe7kHq6Z1aaiIYUdwMUR8YikPSkmB1/R8ATwnwEfiogXJM2neBxX+RmKJ0TEtiqKGY4D18xqU9FTe7eQHiwbEa9IWk/xsNnHS20eLG3yEMXjurLzkIKZ1aaDIYXpkgZLr0Ut9jcbOAp4eJjDngXcVfocwL2SVrfab1XcwzWzWnR448O2iBgYYX9TgduAiyLi5RZtTqAI3LmlxcdHxGZJ+wErJD0REQ+0W1gn3MM1s9pMmjSprddIJE2hCNubIuL2Fm2OAL4OLIiI54eWR8Tm9HUrcAdwTAXfWlMOXDOrTUVXKQi4FlgfEZe3aDMLuB04IyJ+Ulq+RzrRhqQ9gI8BP67o23sLDymYWW0qukrheOAM4DFJa9KyLwCzYOeDab8I7ANclY65Iw1R7A/ckZbtAnwrIu6uoqhmHLij9IEPfKCr7TZu3FhxJWbd/z7WoarJayJiJTDsjiLis8BnmyzfCBz51i16w4FrZrXxrb1mZplMtFt7HbhmVgvPh2tmlpED18wsEweumVkmDlwzs0wcuGZmGUzECcgduGZWG/dwzcwyceCamWXiwDUzy2Ai3vhQy4i1pL0lrZD0VPq6V4t210naKqln06WZWX2qmg+3X9TVw70EuC8iviLpkvT5803aXQ8sBm7MWFtHLrjggq62u/XWWzveJiK6Opb1n257ft3+PtbFPdw8FgA3pPc3AKc2a5Qec7E9V1FmllcVE5D3k7p6uPunJ20SEVvSs4TMbAIZb2Hajp4FrqTvAQc0WXVpj463CFgEMGvWrF4cwswq5sCtSESc2GqdpOckzUi92xnA1gqOtxRYCjAwMODBTrM+MNECt64x3GXAmen9mcCdNdVhZjWq4ioFSTMl3S9pvaR1ki5s0kaSrpC0QdJaSUeX1s2T9GRad0kPvs2d6grcrwAnSXoKOCl9RtKBkpYPNZJ0M/Aj4D2SNkk6q5Zqzaxy7Z4wa6MXvAO4OCIOBY4DzpV0WEOb+cCc9FoEXJ1qmAwsSesPA05rsm1lajlplp4J/9EmyzcDJ5c+n5azLjPLq6KHSG4Bhk7CvyJpPXAQ8Hip2QLgxiiurXxI0rQ0nDkb2JAeJomkW1Lb8raV8Z1mZlabDgJ3uqTB0uel6bxN4/5mA0cBDzesOgh4tvR5U1rWbPmx7RbVKQeumdWmg8DdFhEDI+xrKnAbcFFEvNy4uskmMczynnDgmlltqrpKQdIUirC9KSJub9JkEzCz9PlgYDOwa4vlPTF+blI2s74yNAF5BVcpCLgWWB8Rl7dotgz4dLpa4TjgpTT2uwqYI+kQSbsCC1PbnnAP18xqU1EP93jgDOAxSWvSsi8AswAi4hpgOcUJ+Q3Ar4DPpHU7JJ0H3ANMBq6LiHVVFNWMA3eU5s6d29V2n/vc5zreZsmSJV0dy/pPN78f0P3vY10qukphJc3HYsttAji3xbrlFIHccw5cM6vNRLvTzIFrZrXw5DVmZhmNp8nF2+HANbPauIdrZpaJA9fMLAOP4ZqZZeTANTPLxIFrZpaJr1IwM8vAY7hmZhk5cM3MMnHgmpll4sC1LBYvXpztWFdffXVX27355psVV9LfujnBc84553R1rJy/H3Vy4JqZZTA0AflE4sA1s9q4h2tmlokD18wsEweumVkGVd74IOk64OPA1oh4X5P1fwycnj7uAhwK7BsR2yU9DbwCvAHsGOlx7KMxsUaszWxMqeKpvcn1wLxWKyPiaxHx/oh4P/AnwD9ExPZSkxPS+p6FLbiHa2Y1qqqHGxEPSJrdZvPTgJsrOXCH3MM1s9oMDSuM9AKmSxosvRZ1eby3U/SEbystDuBeSau73W+73MM1s1p0OIa7raJ/7p8C/LBhOOH4iNgsaT9ghaQnIuKBCo71Fu7hmlltOujhVmUhDcMJEbE5fd0K3AEcU+UByxy4ZlabnIEr6Z3Ah4A7S8v2kLTn0HvgY8CPKzlgEx5SMLPaVHVrr6SbgQ9TjPVuAi4DpgBExDWp2SeAeyPil6VN9wfuSKG+C/CtiLi7kqKacOCaWS2q7L1GxGlttLme4vKx8rKNwJGVFNEGB24f6XYGqYULF3a13ZVXXtnxNoODg10d68UXX+xqu2nTpnW8zcBAd+dezj///I63mTt3blfHmih8p5mZWSYTLXBrOWkmaW9JKyQ9lb7u1aTNTEn3S1ovaZ2kC+uo1cx6p4arFGpV11UKlwD3RcQc4L70udEO4OKIOBQ4DjhX0mEZazSzHnPg5rEAuCG9vwE4tbFBRGyJiEfS+1eA9cBB2So0s54amoC8orkU+kJdY7j7R8QWKII13eHRUrpH+ijg4WHaLAIWAcyaNauyQs2sd8ZT77UdPQtcSd8DDmiy6tIO9zOV4r7niyLi5VbtImIpsBRgYGAgOjmGmdXDgVuRiDix1TpJz0makXq3M4CtLdpNoQjbmyLi9h6VamY1mWiBW9fgyDLgzPT+TEq32g1R8V/iWmB9RFyesTYzy6DdE2bjKZTrCtyvACdJego4KX1G0oGSlqc2xwNnAB+RtCa9Tq6nXDPrBZ80yyAingc+2mT5ZuDk9H4lMH7+tJnZW4yn3ms7fKeZmdXGgWtmlsF4G59thwN3Auh2AhVPvGK95sA1M8vEgWtmlsHQrb0TiQPXzGrjHq6ZWSYTLXAnVn/ezMaUqu40k3SdpK2Smj4AUtKHJb1Uuonqi6V18yQ9KWmDpGZTxVbGPVwzq02FPdzrgcXAjcO0+UFEfLzh+JOBJRR3vG4CVklaFhGPV1VYmXu4ZlaLKudSiIgHgO1dlHEMsCEiNkbEr4FbKObr7gkHrpnVpoO5FKZLGiy9FnVxuN+T9KikuyQdnpYdBDxbarOJHj7owEMKZlabDoYUtkVEd49bLjwCvCsiXk2TYH0HmEPz+Vp6Np+2e7hmVptc0zNGxMsR8Wp6vxyYImk6RY92ZqnpwcDmUR+wBfdwzawWOedSkHQA8FxEhKRjKDqbzwMvAnMkHQL8HFgIfKpXdThwzaw2VQWupJuBD1OM9W4CLgOmAETENcAngXMk7QBeAxZGRAA7JJ0H3ANMBq6LiHWVFNWEA9fMalPVrb0RcdoI6xdTXDbWbN1yYHmzdVVTEfLji6RfAM9UvNvpwLaK99lL/VRvP9UK/VVvr2p9V0TsO5odSLqbor52bIuIeaM53lgwLgO3FyQNjvIsaVb9VG8/1Qr9VW8/1ToR+CoFM7NMHLhmZpk4cNu3tO4COtRP9fZTrdBf9fZTreOex3DNzDJxD9fMLBMHrplZJg7cFiTtLWmFpKfS172atJkp6X5J6yWtk3RhHbWmWkasN7UbdqLmHtc47ETPKlyR1q+VdHTuGku1jFTreyX9SNLrkv6ojhob6hmp3tPTz3StpAclHVlHnROdA7e1S4D7ImIOcF/63GgHcHFEHAocB5wr6bCMNZa1Uy8UEzVnv4C8NNHzfOAw4LQmP6v5FDM4zQEWAVdnLTJps9btwAXAf8tc3lu0We/PgA9FxBHAX+CTabVw4La2ALghvb8BOLWxQURsiYhH0vtXgPX0cC7NEYxYL4xqoubRamei5wXAjVF4CJgmaUbuQmmj1ojYGhGrgN/UUF+jdup9MCJeSB8fopgVyzJz4La2f0RsgSJYgf2GayxpNnAU8HDPK2uuo3pr0M5Ez1kngx7GWKmjXZ3WexZwV08rsqYm9OQ1kr4HHNBk1aUd7mcqcBtwUUS8XEVtLY5TSb01aWei56yTQQ9jrNTRrrbrlXQCReDO7WlF1tSEDtyIOLHVOknPSZoREVvSP2u3tmg3hSJsb4qI23tUKlBNvTVqZ6LnrJNBD2Os1NGutuqVdATwdWB+RDyfqTYr8ZBCa8uAM9P7M4E7GxuomMzzWmB9RFyesbZmRqy3ZqtIEz1L2pVioudlDW2WAZ9OVyscB7w0NEySWTu1jiUj1itpFnA7cEZE/KSGGg0gIvxq8gL2oTjb/1T6undafiCwPL2fS/FPt7XAmvQ6eazWmz7fDGyhONmzCTgrY40nAz8BfgpcmpadDZyd3ovibPtPgceAgRr/+49U6wHp5/cyxVMDNgHvGMP1fh14ofR7OlhXrRP55Vt7zcwy8ZCCmVkmDlwzs0wcuGZmmThwzcwyceCamWXiwDUzy8SBa2aWiQPXaiPpX6X5WXeXtEeaU/h9dddl1iu+8cFqJenLwO7A24BNEfFXNZdk1jMOXKtVuvd/FfDPwL+OiDdqLsmsZzykYHXbG5gK7EnR0zUbt9zDtVpJWkbxhIJDgBkRcV7NJZn1zISeD9fqJenTwI6I+FZ6LteDkj4SEf+77trMesE9XDOzTDyGa2aWiQPXzCwTB66ZWSYOXDOzTBy4ZmaZOHDNzDJx4JqZZfL/ACrbrbbDtx4CAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "\n",
    "n = np.sqrt(ms.get_epsilon())\n",
    "print(np.shape(n))\n",
    "nxy = n[:,:,int(sc_z*resolution/2)]\n",
    "\n",
    "pos = ax.imshow(nxy, cmap='gray_r', interpolation='none', extent=[-sc_x/2,sc_x/2,-sc_y/2,sc_y/2] )\n",
    "cbar = fig.colorbar(pos, ax=ax)\n",
    "cbar.set_label('n')\n",
    "ax.set_title('xy cut for z={0}um'.format(sc_z/2))\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epsilon: 2.1025-11.9025, mean 4.25806, harm. mean 2.62689, 100% > 1, 21.9955% \"fill\"\n",
      "(16, 16, 64)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAEWCAYAAAAq1S8mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAX/UlEQVR4nO3df5BlZX3n8fcHRA0/BHFWYAAdYmZVNBrJiEbRYAIWQ5kdTEgWcBWzulMYjVpraqXKnxtTlZjastwsGjJRAiQIcVVwTEYR2WTRRS0GFxFEZSBYjEOYBfkhjlGHfPePe5q6Nrenb8/c+5zp7ver6lTfc85zz/neZvj00895zulUFZKk6dun7wIkabkwcCWpEQNXkhoxcCWpEQNXkhoxcCWpEQNXU5PkxUluTfJQktP6rkfqm4Gr3ZLkwiR/NE+zPwTOq6oDq+qKFnXNJ8lZSb6b5IdJrkhyaN81afkwcDVNTwVu3p03JnnMhGshybOAvwBeDRwG7AA+POnzSHMxcBehJE9L8v0kx3XrK5Pck+TEJL+d5PpZ7d+WZGQPM8mhSf4qybYk9820S/LaJF+a1baS/EKS9cCrgP/SDRd8ZsRxbwN+HvhM1+ZxXZ0bu9q3JPlPQ+3fm+QTSf4myYPAa2cdb2V3nJllR5KF3ib5KuAzVXVNVT0EvAv4zSQHdee4I8lJs2r6m+71qu7z/26SO7vv1TlJnp/kxiT3JzlvgfVomZl4L0LTV1W3JXk7cEmSXwb+Criwqv4xyeOAv0jyzKq6pXvLfwDm+vX/r4GHgGd1X180xvk3JHkRsLWq3jlHm6cluQN4fVV9ASDJpQx6vCuBZwBXJbm9qq7u3rYO+G3gNcDjZh1vG3DgzHqSS+g6DElOAP5uFyW/oqq+1H3Ga4eOeVuSnwD/Frh+rjfP8gJgNfBSYCPwOeAkYD/g/yb5n1X1v8c8lpYZA3eRqqq/TPIbwFeBAv5dt/3HSf6WQci+o/s1ehUjAinJEcBa4ElVdV+3eSphkeRo4AQG4fcvwA1JPsLg1/uZwP3y0Fjvj3ZxrLczCOwTALowPWSMMg4EHpi17QHgoHE/B/C+rv7PJ/khcGlVbe/q+iLwPKb0PdTi55DC4vaXwLOB/1FVPx7afhFwVpIwCLSPz9o/42jg+0NhO00ru3P9YGjbd4Ejh9bvnO8gSdYCbwFOq6o5Q3kODwFPmLXtCcAPRrSdy91Dr380Yv1ApDkYuItUkgOBDwIfBd47fLW9qr4C/AR4CXAWg2GDUe4EDk0yqnf4Q2D/ofMdPmv/QsdPt3XnGu5NPgX43rjHTPJ0Bj9Mfqeq7hza/pJZ47uzl5d0TW8Gnjv0vp9nMHTxnW7Tz3xmYPZnlvaIgbt4/Xfg+qp6PfD3wPmz9l8MnAfs7H7lfpSqugv4LPDhJE9Msl+Sl3a7vw48K8kvJXk88N5Zb7+bwUWxsXQBeS3wx0ken+Q5wOuAS8Z5f5InAJ8G3jn781TVF7upZ3MtX+yaXgL8RhfQBzCYtvapoV73DcAZ3fdhDXD6uJ9PGoeBuwglWQecApzTbfrPwHFJXjXU7K8ZDDfM1bud8Wrgp8C3gO3AWwGq6jsMAukLwK3A7ND+KHBsd3V+3Dm2ZzIYT94GXA68p6quGvO9xwFPBz4w3Hsd870AVNXNDL5nlzD4rAcBvzfU5F3A04D7gP8KfGwhx5fmEx9AvjQl+TkGoXJcVd3adz2S7OEuZW8ArjNspb2H08KWoG7+awCfXyDtRRxSkKRGHFKQpEaW5JDC/vvvX4ccMs6NR5J2x/3338+OHTuyJ8dY4LMwrqyqU/bkfHuDJRm4hxxyCOvXr++7DGnJ2rBhQ+tTrmh9wmlYkoEraXEY3H0+v6VyrcnAldSbffYZ7zLSww8/POVK2jBwJfVm3B7uUmHgSupFEgNXkloxcCWpEQNXkhoxcCWpgSRjz1JYKgxcSb2xhytJjRi4ktSIgStJjRi4ktSAF80kqSF7uJLUiIErSY0YuJLUwHJ8eE2vI9ZJTkny7SRbkpw7Yv+rktzYLdcmeW4fdUqajpnQnW9ZKnrr4SbZF/gQcDKwFbguycaq+uZQs38CfrWq7kuyFtgAvKB9tZKmwVkK7RwPbKmq2wGSXAasAx4J3Kq6dqj9V4CjmlYoaaqWUu91HH3+eDkSuHNofWu3bS6vAz47184k65NsTrJ5x44dEypR0rSMO5ywlEK5zx7uqO/iyL8Ul+RlDAL3hLkOVlUbGAw5sHLlyqXxF+ekJW4phek4+gzcrcDRQ+tHAdtmN0ryHOAjwNqqurdRbZIaMHDbuQ5YneQY4HvAGcBZww2SPAX4FPDqqvpO+xIlTZMXzRqpqp1J3gRcCewLXFBVNyc5p9t/PvBu4EnAh7ufhDurak1fNUuanKU2PjuOXm98qKpNwKZZ284fev164PWt65LUhoErSY0YuJLUiIErSY0st8BdXpcIJe01Zh5APs4yxrEuSLI9yU1z7D8xyQNJbuiWdw/tuyPJN7rtmyf4ER/FHq6k3kywh3shcB5w8S7afLGqXjHHvpdV1T2TKmYuBq6k3kwqcKvqmiSrJnKwKXJIQVJvFvAshRUzz0rplvW7cbpfSfL1JJ9N8qyh7QV8Psn1u3ncsdnDldSLBd74cM8e3vT0NeCpVfVQklOBK4DV3b4XV9W2JE8Grkryraq6Zg/ONSd7uJJ60+ppYVX1YFU91L3eBOyXZEW3vq37uh24nMGjY6fCwJXUm0nNUphPksPTJXeS4xlk371JDkhyULf9AODlwMiZDpPgkIKk3kzqolmSS4ETGYz1bgXeA+wHjzwu4HTgDUl2Aj8CzqiqSnIYcHlXx2OAj1XV5yZS1AgGrqReTPLhNVV15jz7z2MwbWz29tuBZn8r0cCV1JvldqeZgSupNwauJDXiA8glqQEfQC5JDRm4ktSIgStJjRi4ktSIgStJDcw8gHw5MXAl9cYeriQ1YuBKUiMGriQ14I0PktSQgStJjThLQZIasYcrSQ04hitJDRm4ktSIgStJjXjRTJIacAxXkhpaboHba38+ySlJvp1kS5JzR+x/RpIvJ/lxkj/oo0ZJ0zPTy51vWSp66+Em2Rf4EHAysBW4LsnGqvrmULPvA28GTuuhRElTtpTCdBx99nCPB7ZU1e1V9RPgMmDdcIOq2l5V1wE/7aNASdNlD7edI4E7h9a3Ai/oqRZJjfkA8rZG/diq3T5Ysh5YD3DwwQfv7mEkNbSUeq/j6PPHy1bg6KH1o4Btu3uwqtpQVWuqas3++++/x8VJmr7lNqTQZ+BeB6xOckySxwJnABt7rEdSY5MK3CQXJNme5KY59p+Y5IEkN3TLu4f27XK21CT1NqRQVTuTvAm4EtgXuKCqbk5yTrf//CSHA5uBJwD/muStwLFV9WBfdUuajAn3Xi8EzgMu3kWbL1bVK2bVMM5sqYnp9caHqtoEbJq17fyh1//MYKhB0hI0qcCtqmuSrNqNtz4yW6qrZ2a21FQCd3ldIpS0V9lnn33GWibkV5J8Pclnkzyr2zZqttSRkzrhbN7aK6kXCxxSWJFk89D6hqrasIDTfQ14alU9lORU4ApgNROeLTUfA1dSbxYQuPdU1ZrdPc/wdZ+q2pTkw0lWMOHZUvMxcCX1ptWUr+4C/N1VVUmOZzCcei9wP91sKeB7DGZLnTWtOgxcSb2ZVOAmuRQ4kcHQw1bgPcB+8MiF+NOBNyTZCfwIOKOqChg5W2oiRY1g4ErqxSRv7a2qM+fZfx6DaWOj9j1qttS0GLiSerOU7iIbh4ErqTcGriQ1YuBKUiMGriQ1sNSeBDYOA1dSb3wAuSQ1Yg9XkhoxcCWpAcdwJakhA1eSGvGimSQ1Yg9XkhpwDFeSGjJwJakRA1eSGjFwJamBST6AfLEwcCX1xh6uJDVi4EpSIwauJDVi4EpSA974IEkNOUtBkhqxhytJjRi4ktSAY7iS1JCBK0mNeNFMkhpZbj3cXn+8JDklybeTbEly7oj9SfJn3f4bkxzXR52SJm9mDHecZYxjXZBke5Kb5mn3/CQPJzl9aNsdSb6R5IYkmyfw0ebUW+Am2Rf4ELAWOBY4M8mxs5qtBVZ3y3rgz5sWKWmqJhW4wIXAKfOca1/g/cCVI3a/rKp+qarWLPhDLMC8QwpJ3gRcUlX3TfjcxwNbqur27jyXAeuAbw61WQdcXFUFfCXJIUmOqKq7JlyLpB5Makihqq5JsmqeZr8PfBJ4/p6cK8njgN8CVjGUoVX1h/O9d5we7uHAdUk+3g0BTGrQ5UjgzqH1rd22hbYBIMn6JJuTbN6xY8eESpQ0TQvo4a6Y+f+7W9Yv8DxHAq8Ezh+xu4DPJ7l+zON+mkFncCfww6FlXvP2cKvqnUneBbwc+F3gvCQfBz5aVbeNc5I5jAru2o02M3VuADYArFy5cmQbSXuPBT6A/J49/HX/g8Dbq+rhEX3GF1fVtiRPBq5K8q2qumYXxzqqqnY5fDGXsT5t9yv9P3fLTuCJwCeS/OnunLSzFTh6aP0oYNtutJG0SE1wDHc+a4DLktwBnA58OMlpAFW1rfu6HbicwXDnrlyb5Bd3p4h5AzfJm5NcD/wp8H+AX6yqNwC/zGAcY3ddB6xOckySxwJnABtntdkIvKabrfBC4AHHb6Wlo1XgVtUxVbWqqlYBnwB+r6quSHJAkoO6Wg5g8Jv8Lmc6ACcA13czrG7sZjjcOE4d48zDXQH8ZlV9d9YH+NckrxjnJKNU1c7ugtyVwL7ABVV1c5Jzuv3nA5uAU4EtwA4GQxqSlohJXRJKcilwIoOx3q3Ae4D94JEsmcthwOVdHY8BPlZVn5vndGt3t85xxnDfvYt9t+zuibv3b2IQqsPbzh96XcAb9+QckvZeE5ylcOYC2r526PXtwHMXeK7vzt9qNO80k9SLCY7PLhoGrqTe+CwFSWrEHq4kNWLgSlIDjuFKUkMGriQ14kUzSWrEHq4kNeAYriQ1ZOBKUiMGriQ1YuBKUgMLfAD5kmDgSuqNPVxJasTAlaRGDFxJasTAlaQGvPFBkhpyloIkNWIPV5IaMXAlqQHHcCWpIQNXkhrxopkkNeCQgiQ1ZOBKUiMGriQ1stwCd3mNWEvaq8yM4863jHGcC5JsT3LTPO2en+ThJKcPbTslybeTbEly7gQ+1pwMXEm9mHkA+TjLGC4ETpnnfPsC7weunLXtQ8Ba4FjgzCTH7u5nmo+BK6k3k+rhVtU1wPfnafb7wCeB7UPbjge2VNXtVfUT4DJg3W5+nHkZuJJ6M6nAHeM8RwKvBM6ftetI4M6h9a3dtqnwopmk3iwgTFck2Ty0vqGqNizgVB8E3l5VD88656gCagHHXRADV1IvFth7vaeq1uzB6dYAl3XnWwGcmmQngx7t0UPtjgK27cF5dqmXwE1yKPC3wCrgDuB3quq+Ee0uAF4BbK+qZ7esUdL0tbq1t6qOmXmd5ELg76rqiiSPAVYnOQb4HnAGcNa06uhrDPdc4OqqWg1c3a2PciHzXHmUtHhNcFrYpcCXgacn2ZrkdUnOSXLOrt5XVTuBNzGYuXAL8PGqunkCH22kvoYU1gEndq8vAv4RePvsRlV1TZJVrYqS1NakbnyoqjMX0Pa1s9Y3AZsmUsg8+grcw6rqLoCquivJk3uqQ1JPfHjNBCX5AnD4iF3vmNL51gPrAQ4++OBpnELShBm4E1JVJ821L8ndSY7oerdH8LMTkXf3fBuADQArV66c2rQOSZOz3AK3r4tmG4Gzu9dnA5/uqQ5JPZrgrb2LQl+f5E+Ak5PcCpzcrZNkZZJHBq9HXXnspVpJEzfuDIWl1Avu5aJZVd0L/PqI7duAU4fWx77yKGnxWUphOg7vNJPUGwNXkhoxcCWpEQNXkhqYeQD5cmLgSuqNPVxJasTAlaRGDFxJamCp3dQwDgNXUm+8aCZJjdjDlaRGDFxJasAxXElqyMCVpEYMXElqxFkKktSAY7iS1JCBK0mNGLiS1IiBK0mNGLiS1IAPIJekhuzhSlIjyy1wl1d/XtJeZWYu7nzLGMe5IMn2JDfNsX9dkhuT3JBkc5IThvbdkeQbM/sm+PEexR6upF5M+MaHC4HzgIvn2H81sLGqKslzgI8Dzxja/7KqumdSxczFwJXUm0ldNKuqa5Ks2sX+h4ZWDwBqIideIIcUJPVmAUMKK7qhgJll/W6c65VJvgX8PfAfh3YV8Pkk1+/OcRfCHq6k3ixgSOGeqlqzJ+eqqsuBy5O8FHgfcFK368VVtS3Jk4Grknyrqq7Zk3PNxR6upF6M27ud9EyGLkyflmRFt76t+7oduBw4fqInHGLgSupNq8BN8gvpDpTkOOCxwL1JDkhyULf9AODlwMiZDpPgkIKk3kyq95rkUuBEBmO9W4H3APsBVNX5wG8Br0nyU+BHwL/vZiwcxmCYAQZ5+LGq+txEihrBwJXUmwnOUjhznv3vB94/YvvtwHMnUsQYDFxJvfAB5JLUkIErSY0st8DtZZZCkkOTXJXk1u7rE0e0OTrJPyS5JcnNSd7SR62SpqePaWF96mta2LnA1VW1msE9zueOaLMTeFtVPRN4IfDGJMc2rFHSlBm4bawDLupeXwScNrtBVd1VVV/rXv8AuAU4slmFkqZq5gHk4yxLRV9juIdV1V0wCNbulro5dQ+leB7w1V20WQ+sBzj44IMnVqik6VlKvddxTC1wk3wBOHzErncs8DgHAp8E3lpVD87Vrqo2ABsAVq5c2cuTgCQtjIE7IVV10lz7ktyd5Iiud3sEsH2OdvsxCNtLqupTUypVUk+WW+D2NTiyETi7e3028OnZDbr7nj8K3FJVH2hYm6QG+np4TZ/6Ctw/AU5OcitwcrdOkpVJNnVtXgy8Gvi17k9f3JDk1H7KlTQNXjRroKruBX59xPZtwKnd6y8BS+dHm6RHWUq913F4p5mk3hi4ktTAUhufHYeBK6k3Bq4kNWLgSlIDM7f2LicGrqTe2MOVpEYMXElqxMCVpEYMXElqwHm4ktSQsxQkqRF7uJLUiIErSQ04hitJDRm4ktTIcrtolqql9/cWk/w/4LsTPuwK4J4JH3OaFlO9i6lWWFz1TqvWp1bVv9mTAyT5HIP6xnFPVZ2yJ+fbGyzJwJ2GJJurak3fdYxrMdW7mGqFxVXvYqp1OVhe/XlJ6pGBK0mNGLjj29B3AQu0mOpdTLXC4qp3MdW65DmGK0mN2MOVpEYMXElqxMCdQ5JDk1yV5Nbu6xNHtDk6yT8kuSXJzUne0ketXS3z1tu1uyDJ9iQ39VDjKUm+nWRLknNH7E+SP+v235jkuNY1DtUyX63PSPLlJD9O8gd91DirnvnqfVX3Pb0xybVJnttHncudgTu3c4Grq2o1cHW3PttO4G1V9UzghcAbkxzbsMZh49QLcCHQfAJ5kn2BDwFrgWOBM0d8r9YCq7tlPfDnTYvsjFnr94E3A/+tcXmPMma9/wT8alU9B3gfXkzrhYE7t3XARd3ri4DTZjeoqruq6mvd6x8AtwBHNqvwZ81bL0BVXcMgLFo7HthSVbdX1U+AyxjUPGwdcHENfAU4JMkRrQtljFqrantVXQf8tIf6Zhun3mur6r5u9SvAUY1rFAburhxWVXfBIFiBJ++qcZJVwPOAr069stEWVG8PjgTuHFrfyqN/OI3TpoW9pY5xLbTe1wGfnWpFGmlZP7wmyReAw0fsescCj3Mg8EngrVX14CRqm+M8E6m3J6MeCzV7TuI4bVrYW+oY19j1JnkZg8A9YaoVaaRlHbhVddJc+5LcneSIqrqr+7V2+xzt9mMQtpdU1aemVCowmXp7tBU4emj9KGDbbrRpYW+pY1xj1ZvkOcBHgLVVdW+j2jTEIYW5bQTO7l6fDXx6doMMHub5UeCWqvpAw9pGmbfenl0HrE5yTJLHAmcwqHnYRuA13WyFFwIPzAyTNDZOrXuTeetN8hTgU8Crq+o7PdQogKpyGbEAT2Jwtf/W7uuh3faVwKbu9QkMfnW7EbihW07dW+vt1i8F7mJwsWcr8LqGNZ4KfAe4DXhHt+0c4JzudRhcbb8N+Aawpsf//vPVenj3/XsQuL97/YS9uN6PAPcN/Tvd3Fety3nx1l5JasQhBUlqxMCVpEYMXElqxMCVpEYMXElqxMCVpEYMXElqxMBVb5I8v3s+6+OTHNA9U/jZfdclTYs3PqhXSf4IeDzwc8DWqvrjnkuSpsbAVa+6e/+vA/4FeFFVPdxzSdLUOKSgvh0KHAgcxKCnKy1Z9nDVqyQbGfyFgmOAI6rqTT2XJE3Nsn4ervqV5DXAzqr6WPd3ua5N8mtV9b/6rk2aBnu4ktSIY7iS1IiBK0mNGLiS1IiBK0mNGLiS1IiBK0mNGLiS1Mj/B89eJ1G5rYo/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "\n",
    "n = np.sqrt(ms.get_epsilon())\n",
    "print(np.shape(n))\n",
    "nxy = n[:,:,0]\n",
    "\n",
    "pos = ax.imshow(nxy, cmap='gray_r', interpolation='none', extent=[-sc_x/2,sc_x/2,-sc_y/2,sc_y/2] )\n",
    "cbar = fig.colorbar(pos, ax=ax)\n",
    "cbar.set_label('n')\n",
    "ax.set_title('xy cut for z={0}um'.format(0))\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or use ipyvolume (only works on Jupyter Notebook, not Jupyter Lab, see [this](https://github.com/maartenbreddels/ipyvolume/issues/295)) :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(16, 16, 64)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9a17cc7656ab46e1b1336d666edf5d5c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(VBox(children=(HBox(children=(Label(value='levels:'), FloatSlider(value=1.0, max=1.0, step=0.00…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import ipyvolume as ipv\n",
    "ipv.quickvolshow(n, level=[1.45,3.45], extent=[[0, sc_x],[0, sc_y],[0, sc_z]])\n",
    "ipv.pylab.xlim(0,sc_x)\n",
    "ipv.pylab.ylim(0,sc_y)\n",
    "ipv.pylab.zlim(0,sc_z)\n",
    "ipv.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mp [~/home/ELE559/anaconda3/envs/mp/]",
   "language": "python",
   "name": "conda_mp"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}