Week8_IB-2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 879,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enter the wavelength: 2.405\n"
     ]
    }
   ],
   "source": [
    "l = float(input(\"Enter the wavelength: \"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 880,
   "metadata": {},
   "outputs": [],
   "source": [
    "nxsq = 2.72186 + 0.01668 / (l**2 - 0.01546) - 0.00983 * l**2\n",
    "nysq = 2.78252 + 0.01767 / (l**2 - 0.01562) - 0.0111 * l**2\n",
    "nzsq = 2.83314 + 0.01834 / (l**2 - 0.01584) - 0.01179 * l**2\n",
    "\n",
    "nx2wsq = 2.72186 + 0.01668 / ((l/2)**2 - 0.01546) - 0.00983 * (l/2)**2\n",
    "ny2wsq = 2.78252 + 0.01767 / ((l/2)**2 - 0.01562) - 0.0111 * (l/2)**2\n",
    "nz2wsq = 2.83314 + 0.01834 / ((l/2)**2 - 0.01584) - 0.01179 * (l/2)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 881,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.6678945679329873 2.7213805609515305 2.7681258549224537\n",
      "2.719305640233575 2.7788226372930414 2.8289152702723195\n"
     ]
    }
   ],
   "source": [
    "print(nxsq, nysq, nzsq)\n",
    "print(nx2wsq, ny2wsq, nz2wsq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 882,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 883,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 391,
   "metadata": {},
   "outputs": [],
   "source": [
    "def objective(x):\n",
    "    phi = x[0] \n",
    "    # phi is a scalar and not a vector, so no need for the square barakets; I would also simplify the notation and \n",
    "    # call the independent variable phi from the beginning,to avoid any confusion\n",
    "\n",
    "    kx = np.sin(1/2*np.pi) * np.cos(phi)\n",
    "    ky = np.sin(1/2*np.pi) * np.sin (phi)\n",
    "    kz = np.cos(1/2*np.pi)\n",
    "    \n",
    "    a1 = nxsq**(-1)\n",
    "    b1 = nysq**(-1)\n",
    "    c1 = nzsq**(-1)\n",
    "\n",
    "    a2 = nx2wsq**(-1)\n",
    "    b2 = ny2wsq**(-1)\n",
    "    c2 = nz2wsq**(-1)\n",
    "\n",
    "    B1 = (-(kx**2)*(b1+c1) - (ky**2)*(a1+c1) - (kz**2)*(a1+b1))\n",
    "    C1 = (kx**2)*b1*c1 + (ky**2)*a1*c1 + (kz**2)*a1*b1\n",
    "\n",
    "    B2 = (-(kx**2)*(b2+c2) - (ky**2)*(a2+c2) - (kz**2)*(a2+b2))\n",
    "    C2 = (kx**2)*b2*c2 + (ky**2)*a2*c2 + (kz**2)*a2*b2\n",
    "\n",
    "    eq1 = np.sqrt(2)/np.sqrt(-B1 - np.sqrt(B1**2 - 4*C1))\n",
    "    eq2 = np.sqrt(2)/np.sqrt(-B2 + np.sqrt(B2**2 - 4*C2))\n",
    "    \n",
    "    # this needs to be absolute value of the difference\n",
    "    return abs(eq1 - eq2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 392,
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds_phi = (0*np.pi/180, 90*np.pi/180)\n",
    "bounds = [bounds_phi]\n",
    "phi = 20*np.pi/180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 393,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     fun: 7.99948995933164e-11\n",
      "     jac: array([-0.01158281])\n",
      " message: 'Optimization terminated successfully.'\n",
      "    nfev: 28\n",
      "     nit: 6\n",
      "    njev: 6\n",
      "  status: 0\n",
      " success: True\n",
      "       x: array([0.32217152])\n"
     ]
    }
   ],
   "source": [
    "result = minimize(objective, phi, method = 'SLSQP', bounds = bounds)\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 394,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[18.45906857]\n"
     ]
    }
   ],
   "source": [
    "print(result.x *180/np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 895,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd5xU5dn/8c+1jWXZpS+9lwWkqIAIgortUdEEbIndKBF7ijXJkzxJTH9+iTFWsGONqEQsoDE2FBVZOgjSpLelLcsCW6/fHzPrs+ICC+zMmfJ9v17z2pk5Z+Zch2G/c/Y+97lvc3dERCR5pARdgIiIRJeCX0QkySj4RUSSjIJfRCTJKPhFRJKMgl9EJMko+CWhmdlTZvb7ALbbxMwmmVmRmc02s37RrkFkfxT8IpHxd6AQaAo8BvzLzOr0983MepnZe2ZWaGbLzOy8fZafZmaLzWy3mb1vZh2rLetmZjPN7CszO7su65LYp+AXqWNm1hD4HvAbdy9z9wcJ/a6dVofbSAMmAW8Q+nIZAzxrZnnh5c2BicCvwsvzgRervcVvgR8BxwO/rqu6JD4o+CXqzKyrmW0zs/7hx23MbIuZDTezi8xs5j7r32Zmr+7nvYab2Voz+0X4PVaa2WX7WbeJmb1hZgVmtj18v1215R+Y2e/MbFq4iebf4QCtWj7YzD4xsx1mNtfMhu9nF/sBhe6+otpzM4BBtfoHqp2eQBvg7+5e4e7vAdOAK8LLzwcWuvtL7r4X+A1wtJn1DC9PBdLCP9PrsC6JAwp+iTp3Xw7cBTxnZlnAk8BT7v4B8BrQ2cx6VXvJ5cAzB3jLVkBzoC1wFfCImfWoYb2U8LY6Ah2APcAD+6xzKXA10ALIAG4HMLO2wJvA7wkdQd8OvGJmufupZ/U+z60OP/8tZvZQ+Mukptu8/eyz7ee5PuH7vYG5VQvcvRhYHn4e4G7gPmB6+L4kEQW/BMLdHwWWEgqe1sB/h58vIdQkcTmAmfUGOhFq0jiQX7l7ibt/SCigv1fDNre6+yvuvtvdi4A/ACfvs9qT7r7E3fcAE4Bjws9fDkx298nuXunu7xBqPhmxn3oGVA9w4OYD/Fvc6O6N93Pb30nhxcBm4A4zSzez/wrvS1Z4eTahcwzVFQI54W1+4e5Hu3sHd5+0v9okMSn4JUiPEjpCvT8c+FXGA5eamRFqupjg7iVm1sHMdlXdqq2/PXxEW2UVoWaQbzCzLDMbZ2arzGwnMBVobGap1VbbWO3+bkIBCqG/Ei7aJ8yHEfrS2tcGYF71AAceDj9fJ9y9DBgFnBOu+TZCX1Rrw6vsAhru87KGQFFd1SDxS8EvgTCzbOBe4HHgN2bWtGqZu38GlAInEmp6eSb8/Gp3z666VXu7JmbWoNrjDsD6GjZ7G9ADON7dGwInVZVTi5LXAM/sczTewN3/XMO684EW4ROwVXoAn9f0xmY2tvoX2j63hfsryN3nufvJ7t7M3c8EulTbxkLg6GrbaAB0DT8vSU7BL0H5BzDT3X9IqGlm7D7LnybU/l7u7h/X4v1+a2YZZnYicC7wUg3r5BBq198R/qI5lN4szwLfMbMzzSzVzDLDJ5bb7buiu+8E3gF+YWb1zOxcQqH8Xk1v7O7XV/9C2+fWu6bXAJhZv3AdWWZ2O6G/Pp4KL/4X0MfMLjCzTOB/CP0VsvgQ9lkSlIJfos7MRgJnAdeHn7oV6L9Pb5xnCDUDHeikbpWNwHZCR/nPAdfvJ+DuBeoDW4DPgLdqW7O7rwFGAr8ACgj9BXAH+/8duhUYCGwjdEL4++5eWdvt1dIVhJqPNhPqKnpGVZOZuxcAFxA6j7GdULfNi+t4+xKnTBOxSCwys/qEAq2/uy89wHrDgWfd/VtH3iJSMx3xS6y6AZhxoNAXkcOTdvBVRKLLzFYSOuE6KuBSRBKSmnpERJKMmnpERJJMXDT1NG/e3Dt16hR0GSIicWXmzJlb3P1bw4rERfB36tSJ/Pz8oMsQEYkrZraqpufV1CMikmQU/CIiSUbBLyKSZBT8IiJJRsEvIpJkFPwiIklGwS8ikmQSOvgnzVnHs5/V2I1VRCSmLdlUxD3vLKGgqOTgKx+ihA7+KfM38sTHXwVdhojIIZu3tpD73l3KntKKOn/vhA7+vFY5rNxazN6yuv+HExGJpM1FewFonpNR5++d2MHfMptKhxUFxQdfWUQkhhQUlZBdL42sjLofWSfBgz8HCLWViYjEk4KiEnJz6kXkvRM6+Ds1a0B6qin4RSTuFBSVkJut4D9kGWkpdG7eQMEvInGnYJeO+A9bXssclmzaFXQZIiKHRE09RyCvZQ6rt+1md2l50KWIiNTK3rIKivaWK/gPV9UJ3mWbddQvIvGh6qItBf9hymuZDcCXG9XOLyLxoWCXgv+IdGzWgIy0FJbqiF9E4sTXR/zq1XN4UlOMbrnZOuIXkbhRFfwtdMR/+Hq0ymGpunSKSJzYXFSCGTRtUPfDNUCSBH/3ltmsL9zLzr1lQZciInJQBUUlNGuQQVpqZCI6KYK/R7hnz1L15xeROFBQVELzCLXvQ5IEv8bsEZF4EsmrdiFJgr9t4/rUT09V8ItIXNgSwat2IUmCPyXFyGuZreAXkZjn7hQUldAiJzNi20iK4AforjF7RCQO7NxTTmlFZfwe8ZvZT81soZktMLMXzCzTzDqb2XQzW2pmL5pZZPor7aNHyxwKikrYXlwajc2JiByWgl2hmbfiMvjNrC3wI2Cgu/cBUoGLgb8Af3f37sB2YHSkaqiue3joBjX3iEgs2xzhq3Yh8k09aUB9M0sDsoANwKnAy+Hl44FREa4BCF3EBQp+EYltkR6gDSIY/O6+DvgrsJpQ4BcCM4Ed7l41RvJaoG1NrzezMWaWb2b5BQUFR1xPq4aZ5NRLUzu/iMS0uA5+M2sCjAQ6A22ABsDZNazqNb3e3R9x94HuPjA3N7cu6iGvVQ5f6ohfRGJYQVEJGWkpNMys+0nWq0Syqed04Ct3L3D3MmAicALQONz0A9AOWB/BGr4hr2U2SzcV4V7jd42ISOCq5to1s4htI5LBvxoYbGZZFtqD04AvgPeBC8PrXAVMimAN35DXMoftu8u+HutaRCTWRPqqXYhsG/90QidxZwHzw9t6BLgLuNXMlgHNgMcjVcO+8jRmj4jEuEjOtVslco1IgLv/Gvj1Pk+vAAZFcrv7UxX8X24sYmi35kGUICJyQAVFJfTv2CSi20iaK3cBmmdn0CQrnaWbdYJXRGJPWUUl23aXRmwClipJFfxmRl7LHM3GJSIxaVtxKe6R7coJSRb8EGruWbppl3r2iEjMifRcu1WSMPizKSopZ0Ph3qBLERH5hmhcvAVJGfwaukFEYpOCP0LUpVNEYtXmolBLRCSnXYQkDP4mDTLIzamnoRtEJOYUFJXQMDONzPTUiG4n6YIf/m/oBhGRWBKNq3YhaYM/NBtXZaV69ohI7IjGVbuQxMG/p6yCdTv2BF2KiMjXQsEfubl2qyRt8AO6kEtEYkpoknUd8UfE19MwaugGEYkRxSXlFJdWqKknUhpmptOmUSZLdMQvIjFiy67oXLULSRr8AN3DJ3hFRGJBtC7egiQO/h6tclhWsIsK9ewRkRig4I+C7i2yKS2vZNXW4qBLERH5emZABX8EacweEYklm3eWkJpiNMnKiPi2kjb4v+7Zo3Z+EYkBBUUlNGuQQWpK5CZZr5K0wZ+VkUb7pvU1Zo+IxIRoDdcASRz8AD1a5mjMHhGJCdEargGSPPi7t8xhRUExpeWVQZciIkmuoKgkKn34IcmDv0fLHMornZXq2SMiAaqsdLbsKqFFQwV/xFWd4NWYPSISpB17yiivdB3xR0PX3GxSDLXzi0ig/u/irciPzAlJHvyZ6al0atZAPXtEJFDRvGoXkjz4IXQhl+bfFZEgFewKzbWr4I+SvJbZrNxazN6yiqBLEZEktXmnjvijKq9VDpUOywt01C8iwSgoKqF+eioNMiI7yXoVBb/G7BGRgFVdtWsW+eEaQMFPp2YNSEsxjdkjIoGJ5lW7oOAnIy2FLrkNNBuXiAQmmlftgoIfCM/Gpfl3RSQg0RygDRT8QGjohjXb9lBcUh50KSKSZErKK9ixu4wWCv7oygsP3bBss9r5RSS6tu4qBaLXlRMU/MD/9ezRFbwiEm3RvmoXIhz8ZtbYzF42s8VmtsjMhphZUzN7x8yWhn82iWQNtdGxWQMy0lI0Zo+IRF3CBT/wD+Atd+8JHA0sAn4GvOvu3YF3w48DlZpidMvNVpdOEYm6peEm5uaJ0KvHzBoCJwGPA7h7qbvvAEYC48OrjQdGRaqGQ5HXMlsXcYlIVH2ybAt/f2cJgzo3pXWj6IzMCZE94u8CFABPmtlsM3vMzBoALd19A0D4Z4uaXmxmY8ws38zyCwoKIlhmSF6rHDYU7mXn3rKIb0tEZOH6QsY8M5NOzbN49IqBUbtqFyIb/GlAf+Bhdz8WKOYQmnXc/RF3H+juA3NzcyNV49fyWoRO8KqdX0Qibc223fzgyRnkZKYx/ppBNMpKj+r2Ixn8a4G17j49/PhlQl8Em8ysNUD45+YI1lBrPVpVjdmjdn4RiZytu0q48onPKS2v5OlrBtG6Uf2o1xCx4Hf3jcAaM+sRfuo04AvgNeCq8HNXAZMiVcOhaNu4PvXTUzUNo4hETHFJOdc8NYP1O/bw+FUD6R7uSh5taRF+/1uA58wsA1gBXE3oy2aCmY0GVgMXRbiGWklJMbq3zGaphm4QkQgoq6jkpudnMX9dIeOuGMjATk0DqyWiwe/uc4CBNSw6LZLbPVx5LXP44MvNuHtUT7SISGIrr6jk1glz+eDLAv50fl/OOKploPXoyt1qjuvUhC27Slm4fmfQpYhIgqgK/dfnrudnZ/fkkkEdgi5JwV/dGUe1IjXFmDx/Q9CliEgCKK+o5LaX5vLa3PXcdVZPrj+5a9AlAQr+b2jaIIMhXZoxef4G3D3ockQkjlVUOre9NJdJc9Zz51k9uGF4bIQ+KPi/5ey+rVi5dTeLNugkr4gcnopK57YJc5g0Zz13nNmDG4d3C7qkb1Dw7+PM3q1IMZiyQM09InLoqkL/1XDo33RKbIU+KPi/pXl2PY7v3Iw31dwjIoeootK5/aW5MR36oOCv0Yi+rVhRUPz1qHkiIgdTUenc8dJc/jV7Hbf/V17Mhj4o+Gt0Zp9WmMGb89TcIyIHV1Hp3PHyXCbOXsdtZ+Rx86ndgy7pgBT8NWiRk8lxnZqqnV9EDqqi0rnz5XlMnLWOW8/I45bTYjv0QcG/XyP6tGLJpl0s0xAOIrIfFZXOXa/M45VZa/np6Xn8KA5CHxT8+3VWn9YATJm/MeBKRCQWVYX+yzPX8pPTu/Pj0+Mj9OEgY/WYWSZwLnAi0AbYAywA3nT3hZEvLzitGmUyoGMT3py/IS7+dBOR6CmrqOSul+cxcfY6fnJ6d35yel7QJR2S/R7xm9lvgGnAEGA6MA6YAJQDfw5PlN4vGkUGZUTf1izeWMSKAvXuEZGQ3aXljHk6/+sTufEW+nDgI/4Z7v6b/Sy7x8xaAMGPNhRBZ/Vpxe/e+IIpCzbGdNcsEYmO7cWlXDN+BnPX7OCP5/Xl0uPjMwL3e8Tv7m8e6IXuvtnd8+u+pNjRtnF9jmnfWL17RIR1O/Zw4dhPWLh+Jw9dNiBuQx9qMR6/mb0O7HsJayGQD4xz972RKCxWnNO3NX+YvIjVW3fToVlW0OWISAC+3FjEVU98TnFpOc9cM4jjuzQLuqQjUptePSuAXcCj4dtOYBOQF36c0M7q0wqAyTrqF0lKM1Zu46Kxn1DpzoTrhsR96EPtZuA61t1Pqvb4dTOb6u4nmVlC9+wBaN80i37tGjFl/oaYGUtbRKLjnS82cfPzs2jbuD7jrxlE+6aJ8Vd/bY74c83s68as8P3m4YelEakqxpzdpzVz1xaydvvuoEsRkSh5ccZqrnsmn56tcnjp+iEJE/pQu+C/DfjYzN43sw+Aj4A7zKwBMD6SxcWKEX1DzT26mEsk8bk7D76/jLtemc+w7rk8f+1gmmXXC7qsOnXQph53n2xm3YGegAGLq53QvTeSxcWKjs0a0LtNQyYv2MC1J3UJuhwRiZDKSufuN77gqU9Wct6xbfnLBf3ISEu8AQ4OukdmlgXcAdzs7nOA9mZ2bsQrizEj+rZm9uodrN+xJ+hSRCQCSsoruOWfs3nqk5X8cFhn/nbR0QkZ+lC7pp4nCbXlDwk/Xgv8PmIVxaizw7173lqg5h6RRFO4u4yrn5zBm/M28IsRPfnluUeRkmJBlxUxtQn+ru7+v0AZgLvvIdTkk1S65GbTs1UOk+erW6dIIlm5pZjzHp7GjJXbuOd7RzPmpMTvvVeb4C81s/qEL+Iys65ASUSrilEj+rYmf9V2NhYm9DVrIklj+oqtjHpoGtuKS3l29PGc379d0CVFRW2C/9fAW4Ta9p8D3gXujGhVMaqqd8/bC9XcIxLvXp65lssfn07TBhm8euPQhLgwq7Zq06vnHTObBQwm1MTzY3ffEvHKYlC3Fjnktcxm8vwNXHVCp6DLEZHDUFnp/PXfX/LQB8sZ2q0ZD106gEZZ6UGXFVX7DX4z67/PU1WN2x3MrIO7z4pcWbHr7D6tue+9pWwu2kuLnMygyxGRQ7CntIJbJ8xhyoKNXDKoA3eP7E16amL23DmQAx3x/y38MxMYCMwldMTfj9D4/MMiW1psGtG3Nf94dylvL9zEFYM7Bl2OiNTS5p17+eHT+cxfV8gvz+nF6GGdMUu6firAgYdlPsXdTwFWAf3dfaC7DwCOBZZFq8BYk9cymy65DZii3j0icWPh+kJGPjiNZZt38egVA/nhiV2SNvShdid3e7r7/KoH7r4AOCZyJcU2M+Ocvq35bMVWtu5Kys5NInHlnS82cdHYTwF4+foTOP2olgFXFLzaBP8iM3vMzIab2clm9iiwKNKFxbKz+7Sm0uHthZuCLkVE9sPdeXTqCsY8k0+3FtlMumkoR7VpGHRZMaE2wX81sBD4MfAT4Ivwc0mrV+scOjXL0sxcIjGqtLySX/xrPn+YvIiz+7TixTFDaNFQnTGq1KY7517g7+GbEGruGdG3NeOmrmB7cSlNGmQEXZKIhG0u2suNz84if9V2bjqlK7ed0SOhh184HPs94jez183sO2b2rQ6uZtbFzO42s2siW17sGtG3NRWVzjtfqLlHJFbMXbOD794/jQXrC7n/kmO548yeCv0aHKip51rgRGCxmc0ws8lm9p6ZfQWMA2a6+xMH24CZpZrZbDN7I/y4s5lNN7OlZvaimcXl4XLvNg1p37Q+b6p3j0hMeGXmWi4a9ympKcYrN5zAd45uE3RJMWu/TT3uvpHQ0Ax3mlknoDWwB1ji7ocyFdWPCZ0Mrjqr8hfg7+7+TzMbC4wGHj700oNlZozo05rHP/6KbcWlNFVzj0ggyisq+cPkRTw5bSVDujTjwcv66/fxIGp1yZq7r3T3T919zqGEvpm1A84BHgs/NuBU4OXwKuOBUYdWcuy4YEA7yiudZz5dFXQpIklpW3EpVz7xOU9OW8nVQzvx9OhBCv1aiPS1yvcS+quhMvy4GbDD3cvDj9cCbWt6oZmNMbN8M8svKCiIcJmHJ69lDqf3asFTn3zF7tLyg79AROrMF+t38t0HPiZ/1Xb+etHR/Po7yTn8wuGI2L9SeJauze4+s/rTNazqNb3e3R8JXy08MDc3NyI11oUbhndl++4yJsxYE3QpIknj9bnrOf/haZRXOBOuG8KFA5JjOOW6UpupF83M8syshx3aNc5Dge+a2Urgn4SaeO4FGptZ1bmFdsD6Q6w5pgzo2JTjOjXh0Y++oqyi8uAvEJHDVlHp/HnKYm55YTa92zTitVuGckz7xkGXFXcOGPxmdhShi7cmhm+Lws8dlLv/3N3buXsn4GLgPXe/DHgfuDC82lXApMOsPWbcMLwr63bs4fW5cf0dJhLTCneXcc1TMxj74XIuGdSBF64drBFyD9PBjvifBH7h7n3cvTdwF/DcEW7zLuBWM1tGqM3/8SN8v8Cd0qMFPVrmMPbD5VRW1thyJSJHYMmmIkY9NI1Plm/hD+f14U/n903YidCj4UAXcN1A6MRrQzO70syuBBoDzczsxkPZiLt/4O7nhu+vcPdB7t7N3S9y97gf6czMuGF4V5Zs2sV7izcHXY5IQnlj3npGPTiNor3lPH/tYC47XsOhH6kDfWUOANKB46rdBgIZ4WVSzbn9WtOuSX0e+mAZ7jrqFzlS5RWV/HHyIm5+fjY9W+Xw5o+GcVynpkGXlRAOdAHXD81sOfAPd18GoaEagHPdfXS0CowXaakpjDmpC/8zaSEzVm5nUGf9BxU5XFt3lXDz87P5dMVWrhjckV+de5SadurQwf4lbwU+NLMHzOx+YCrw08iXFZ8uGtCepg0yGPvh8qBLEYlbc9fs4Dv3f8zM1dv5fxf243ej+ij069gB/zXdfRKh8XoWEOrdM9zdX41GYfGofkYqV5/QifcWb2bRhp1BlyMSd16csZqLxn6KmTHxhhO4aGD7oEtKSAf9Gg2fjB0bviXtlIu1deWQTjTISGWcjvpFaq2kvIKfT5zPXa/M5/guTXnjlmH0adso6LISlv5+qmONstK59PgOvD5vA2u2HcpYdiLJaUPhHr437jNe+Hw1Nw7vylNXD9IcFxGm4I+A0cO6kGLw6Ecrgi5FJKZ9unwr5973Mcs2FTH28v7ceVZPUjV+fsQp+COgVaNMzj+2HS/OWMMWTcgu8i3uzmMfreDyx6fTOCudSTcP5aw+rYMuK2ko+CNkzMldKK2o5KlpK4MuRSSmFJeUc8sLs/n9m4s4vVcLXr1pKN1a5ARdVlJR8EdI19xszjyqFU9/upJdJRqyWQRgecEuRj04jcnzN3DnWT0Ye/kAcjK/NburRJiCP4KuH96VnXvLeWH66qBLEQncWws2MPKBaWwtLuWZ0cdz4/BuHNqAv1JXFPwRdEz7xpzQtRmPfbyCkvKKoMsRCUR5RSV/nrKY65+dRdcW2bxxyzCGdmsedFlJTcEfYTcM78qmnSW8Ontd0KWIRN2WXSVc8fjnjP1wOZcd34EJ1w2mTeP6QZeV9BT8ETasW3P6tG3IuA9XUKEhmyWJzFq9nXPv+5hZ4aEX/nBeX+qlpQZdlqDgjzgz44aTu7FiSzH/Xrgx6HJEIs7deeazVXx/3KekpxkTb9TQC7FGwR8FZ/VpRadmWTz84XIN2SwJbU9pBbe9NJdfvbqAYd2a88bNJ9K7jYZeiDUK/ihITTHGnNSVeWsL+WT51qDLEYmIVVuLOf/hT/jX7HX89PQ8Hr/qOBplqatmLFLwR8n5/duSm1NPQzZLQnp/8Wa+c//HrN+xhyd+cBw/Pr07KRp6IWYp+KMkMz2V0cM689HSLcxdsyPockTqRGWlc/+7S7lm/AzaNcnijVuGcUqPFkGXJQeh4I+iy47vQJOsdP737cVq65e4V7S3jOufncnf3lnCqGPaMvHGE2jfNCvosqQWFPxRlJOZzi2ndmfasq1MXbol6HJEDlvV0AvvLt7Mr79zFPd872gy09VVM14o+KPsssEdaN+0Pn+esphK9euXOPTOF5sY9cA0duwu49nRx3P10M4aeiHOKPijrF5aKrf/Vw8WbdjJpLm6mlfiR2Wlc+9/lnDt0/l0at6A124ZxpCuzYIuSw6Dgj8A3+nXhj5tG/LXt5ewt0xj+Ejs27m3jDHP5HPvf5ZyQf92vHT9ENpq6IW4peAPQEqK8bOzerFuxx6e/WxV0OWIHNCyzUWMemAaH3xZwN0je/PXi/qpPT/OKfgDMqx7c07s3pwH3l9G4Z6yoMsRqdHbCzcy8oFp7NxbxvPXDubKIZ3Unp8AFPwB+tnZPSncU8bDH+iiLoktlZXO3/79Jdc9M5NuLXN4/ZZhDOrcNOiypI4o+APUu00jRh3TlienfcX6HXuCLkcECPXPv/bpfO5/bxnfG9iOF8cMpnUjtecnEgV/wG49Iw93uPc/S4IuRYQ123ZzwcOf8MGSAn43sjd/uUDt+YlIwR+w9k2zuGJIR16euZYlm4qCLkeS2IyV2xj54DQ2Fu7l6WsGcYXa8xOWgj8G3HxKNxrUS+MvUxYHXYokqZfy13Dpo5/RuH46r940VFMjJjgFfwxo0iCDG4Z35d3Fm5m+QsM2S/RUVDp/mryIO16ex6DOTfnXjUPpkpsddFkSYQr+GHHN0M60apjJn6ZoADeJjl0l5Vz3zEzGTV3B5YM78NTVgzR+fpJQ8MeIzPRUbj0jjzlrdjBlgaZolMhau303Fz78Ce9/uZm7R/bm96P6kp6qOEgW+qRjyAUD2pHXMpv/9/aXlFVUBl2OJKiZq7Yz6sFprNuxhyd/cBxXDukUdEkSZRELfjNrb2bvm9kiM1toZj8OP9/UzN4xs6Xhn00iVUO8SU0x7jyzJ19tKeafM9YEXY4koImz1nLJI5+RXS+Nf904lJPycoMuSQIQySP+cuA2d+8FDAZuMrOjgJ8B77p7d+Dd8GMJO61XCwZ1aso//rOU4pLyoMuRBFFZ6fzvW4u5dcJcBnRswqs3DaVbC53ETVYRC3533+Dus8L3i4BFQFtgJDA+vNp4YFSkaohHZsbPRvRky64SHv1oRdDlSALYW1bBjc/N4qEPlnPJoA48PXoQjbMygi5LAhSVNn4z6wQcC0wHWrr7Bgh9OQA1TtBpZmPMLN/M8gsKCqJRZszo36EJZ/dpxSNTV1BQVBJ0ORLHdpWUc/WTM3hr4UZ+eU4v/nheH53ElcgHv5llA68AP3H3nbV9nbs/4u4D3X1gbm7ytUPecWYPSsorue/dpUGXInFqe3Eplz36GZ+v3Ma93z+GH57YRVfiChDh4DezdEKh/5y7Tww/vcnMWoeXtwY2R7KGeNUlN5tLBrXnhc9X89WW4qDLkTizsXAv3xv3KYs2FjHu8gGMOrZt0CVJDIlkrx4DHgcWufs91Ra9BlwVvn8VMClSNcS7H53WnfTUFO5/T0f9UnsrtxRz4dhP2FC4l/FXD+L0o1oGXZLEmEge8Q8FrgBONbM54dsI4M/AGWa2FDgj/Fhq0CInk+8f1xx3YFYAAAxNSURBVJ7X5qxnQ6GGbZaDW7RhJxeO/ZTiknKev/Z4zYkrNYpkr56P3d3cvZ+7HxO+TXb3re5+mrt3D//cFqkaEsHoYZ1x4MlpK4MuRWLczFXb+f64T0lLMV66fgj92jUOuiSJUTq9H+PaN81iRN/WPD99NTv3aopGqdnUJQVc/th0mjbI4KXrh9CtRU7QJUkMU/DHgetO6sKuknKen7466FIkBk2ev4HR42fQsVkWE64fQvumWUGXJDFOwR8H+rRtxNBuzXhy2leUlmsMH/k/E2as4ebnZ9GvXWNevG4ILXIygy5J4oCCP06MOakrm3aWMGnOuqBLkRjx6NQV3PnKPIZ1z+WZ0YNoVF9DKkvtKPjjxEndm9OzVQ6PTF1BZaXG609m7s7f/v0lf5i8iHP6tuaxKweSlZEWdFkSRxT8ccLMuO7kLizdvIsPluiat2Q2/pOV3P/eMr4/sD33XXIsGWn6NZZDo/8xceTcfm1o0yiTcR9q8LZk9fHSLfzuzUWc3qslfzq/L6kpGoJBDp2CP46kp6ZwzbDOTP9qG3PW7Ai6HImyr7YUc+NzM+mWm829Fx9DikJfDpOCP85cPKgDOZlpPDJ1edClSBTt3FvGD8fPIC01hceuGkh2PbXpy+FT8MeZ7HppXD64I28t2MhKDd6WFCoqnVuen82qrbt56LL+6qcvR0zBH4euPqETaSkpPPax2vqTwZ+nLOLDJQXcPbIPg7to7B05cgr+ONSiYSbnHduWl/LXsnWXJmpJZC/lr+HRj77iqiEdufT4DkGXIwlCwR+nrj2pMyXllTz96aqgS5EImblqG//9rwUM7daMX517VNDlSAJR8Mepbi1yOL1XS57+dCV7SiuCLkfq2Lode7jumZm0aZzJg5f2J03TJUod0v+mOHbdyV3YvruMl2auCboUqUO7S8u5dnw+JWWVPHbVQE2MLnVOwR/HBnZsQv8OjXnso68or9DgbYmgstK5/aW5LN64k/suPVbDK0tEKPjjmJkx5qSurN62m7cWbgy6HKkD9723lMnzN/Lzs3txSo8WQZcjCUrBH+fOOKolnZs34JGpK3DX4G3xbMr8Ddz7n6VcOKAdPzyxc9DlSAJT8Me51BTj2hO7MG9tIZ+t0CyW8WrBukJunTCX/h0a84fz+mCm4RgkchT8CeD8/m1pnp3BOA3jEJe27iphzNP5NM5KZ+wVA6iXlhp0SZLgFPwJIDM9lauGdOKDLwv4cmNR0OXIIXB3fj5xPlt2lfLolQM1g5ZEhYI/QVwxpCP101N5ZKqGcYgnE2et499fbOL2M/Po07ZR0OVIklDwJ4jGWRl8/7j2TJqzjg2Fe4IuR2ph/Y49/Oa1hQzq1JTRw7oEXY4kEQV/Ahk9rDMO3P/esqBLkYOorHTufHkeFe789aKjNaGKRJWCP4G0b5rFD07oxPPTV/PsZxrDJ5Y9O30VHy/bwn+f04sOzTTMskSXZnNIMD8/uycrtxTzP5MW0KZxJqf2bBl0SbKPr7YU88fJizg5L5dLB2nETYk+HfEnmLTUFO6/9Fh6t2nEzc/PZv7awqBLkmrKKyq5dcIc6qWl8pcL+qm/vgRCwZ+AsjLSePwHA2mSlcE142ewdvvuoEuSsHFTVzB79Q7uHtmbVo3UdVOCoeBPUC1yMnny6uPYW1bB1U/OoHBPWdAlJb0v1u/k3v8s4Zy+rfnu0W2CLkeSmII/geW1zGHcFQNYubWYG56dSWm5RvAMSkl5BbdOmEOj+hn8bpSGZJBgKfgT3Aldm/OXC/rxyfKt/GziPA3kFpB//GcpizcW8ZcL+tK0gcbXl2CpV08SOL9/O9Zu38M97yyhfZMsfnpGXtAlJZWZq7Yz9sPlfH9ge07rpV5WEjwFf5K45dRurN62m3+8u5R2Tepz0cD2QZeUFHaXlnPbhDm0blSfX57bK+hyRAAFf9IwM/50fl82Fu7l5xPn07pRfYZ1bx50WQnvz1MWs3Lrbl64djA5melBlyMCqI0/qaSnpvDQ5f3p1iKbG56dyeKNO4MuKaF9vHQLT3+6imuGdmZI12ZBlyPytUCC38zOMrMvzWyZmf0siBqSVcPMdJ74wXHUz0jlmidnsGnn3qBLSkiFe8q44+W5dM1twJ1n9Qi6HJFviHrwm1kq8CBwNnAUcImZHRXtOpJZm8b1eeIHx1G4p4xrnppBcUl50CUlnN++vpDNRSXc871jyEzXxCoSW4I44h8ELHP3Fe5eCvwTGBlAHUmtT9tGPHBZfxZvLOKm52dRXqE+/nXlrQUbmThrHTed0o2j2zcOuhyRbwki+NsCa6o9Xht+7hvMbIyZ5ZtZfkFBQdSKSyan9GjB70b2oVtuNim6oKjO5GSmcXqvltxyaregSxGpURC9empKmG9dVeTujwCPAAwcOFBXHUXIpcdrdMi6NrRbc4Z2U48piV1BHPGvBap3Im8HrA+gDhGRpBRE8M8AuptZZzPLAC4GXgugDhGRpBT1ph53Lzezm4G3gVTgCXdfGO06RESSVSBX7rr7ZGByENsWEUl2unJXRCTJKPhFRJKMgl9EJMko+EVEkozFw4xMZlYArAq6jjrUHNgSdBFRlGz7C9rnZBHr+9zR3XP3fTIugj/RmFm+uw8Muo5oSbb9Be1zsojXfVZTj4hIklHwi4gkGQV/MB4JuoAoS7b9Be1zsojLfVYbv4hIktERv4hIklHwi4gkGQV/hJjZE2a22cwW7Ge5mdl94Qnn55lZ/2jXWNdqsc/DzazQzOaEb/8T7Rrrkpm1N7P3zWyRmS00sx/XsE5Cfc613OdE+5wzzexzM5sb3uff1rBOPTN7Mfw5TzezTtGv9BC4u24RuAEnAf2BBftZPgKYQmhGssHA9KBrjsI+DwfeCLrOOtzf1kD/8P0cYAlwVCJ/zrXc50T7nA3IDt9PB6YDg/dZ50ZgbPj+xcCLQdd9oJuO+CPE3acC2w6wykjgaQ/5DGhsZq2jU11k1GKfE4q7b3D3WeH7RcAivj1/dEJ9zrXc54QS/ux2hR+mh2/79ooZCYwP338ZOM0sdieyVvAHp1aTziegIeE/maeYWe+gi6kr4T/tjyV0NFhdwn7OB9hnSLDP2cxSzWwOsBl4x933+zm7ezlQCDSLbpW1p+APTq0mnU8wswiNHXI0cD/wasD11AkzywZeAX7i7jv3XVzDS+L+cz7IPifc5+zuFe5+DKE5wgeZWZ99Vomrz1nBH5ykm3Te3XdW/cnsoVnY0s2secBlHREzSycUgM+5+8QaVkm4z/lg+5yIn3MVd98BfACctc+irz9nM0sDGhHDzZ4K/uC8BlwZ7vUxGCh09w1BFxVJZtaqqt3TzAYR+v+3NdiqDl94Xx4HFrn7PftZLaE+59rscwJ+zrlm1jh8vz5wOrB4n9VeA64K378QeM/DZ3pjUSBz7iYDM3uBUO+G5ma2Fvg1oZNCuPtYQnMOjwCWAbuBq4OptO7UYp8vBG4ws3JgD3BxLP9y1MJQ4Apgfrj9F+AXQAdI2M+5NvucaJ9za2C8maUS+hKb4O5vmNndQL67v0boy/AZM1tG6Ej/4uDKPTgN2SAikmTU1CMikmQU/CIiSUbBLyKSZBT8IiJJRsEvIpJkFPyS9MzsPTObHL4waX/r1DezD8Nd+o50e7lm9taRvo/I4VLwS9Jz91OBEuCcA6x2DTDR3SvqYHsFwAYzG3qk7yVyOBT8IiFTgMsOsPwyYBJ8Pd78G1ULzOwBM/tB+P5KM/ujmX1qZvlm1t/M3jaz5WZ2fbX3e/Ug2xOJGAW/SMjFwH+ZWcN9F5hZBtDF3VfW8r3WuPsQ4CPgKUJXsg4G7q62Tj5w4pEULHK4FPyS9MysL6FBtZ4HLqhhlebAjkN4y9fCP+cTmnilKNy8s7dqzBdCw/u2OcySRY6Igl8EfgL8HXiWmptf9gCZ+zxXfRjefU8Kl4R/Vla7X/W4anyszPD7ikSdgl+SmpnlAicTmipvGtDRzL5xJO7u24FUM6se/n3CPX0ygCHAofb2yQNqnJtYJNIU/JLsrgMec/ey8OMXqHlkxX8Dw6o93g68BXwKvAv83MwaHMJ2TwHePPRyRY6cRucUqQUzOxa41d2vMLPhwO3ufu4RvN9UYGT4rwmRqNIRv0gtuPts4P26uoALuEehL0HREb+ISJLREb+ISJJR8IuIJBkFv4hIklHwi4gkGQW/iEiS+f/upQcWkh5STQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = [0.7,0.8,0.9,1,1.1,1.2,1.3,1.9,2,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3,3.1,3.15,3.2]\n",
    "y = [90,57.36,40.52,28.68,18.46,6.37,0,-np.inf,0,8.1,14.82,19.99,24.62,29.04,33.37,37.74,42.21,46.9,51.92,57.48,60.56,90]\n",
    "plt.plot(x,y)\n",
    "plt.xlabel('λ (μm)')\n",
    "plt.ylabel('Φ (deg)')\n",
    "plt.title('xy-plane θ = 90°')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 601,
   "metadata": {},
   "outputs": [],
   "source": [
    "def objective(x):\n",
    "    theta = x[0] \n",
    "    \n",
    "    kx = np.sin(theta) * np.cos(1/2*np.pi)\n",
    "    ky = np.sin(theta) * np.sin(1/2*np.pi)\n",
    "    kz = np.cos(theta)\n",
    "    \n",
    "    a1 = nxsq**(-1)\n",
    "    b1 = nysq**(-1)\n",
    "    c1 = nzsq**(-1)\n",
    "\n",
    "    a2 = nx2wsq**(-1)\n",
    "    b2 = ny2wsq**(-1)\n",
    "    c2 = nz2wsq**(-1)\n",
    "\n",
    "    B1 = (-(kx**2)*(b1+c1) - (ky**2)*(a1+c1) - (kz**2)*(a1+b1))\n",
    "    C1 = (kx**2)*b1*c1 + (ky**2)*a1*c1 + (kz**2)*a1*b1\n",
    "\n",
    "    B2 = (-(kx**2)*(b2+c2) - (ky**2)*(a2+c2) - (kz**2)*(a2+b2))\n",
    "    C2 = (kx**2)*b2*c2 + (ky**2)*a2*c2 + (kz**2)*a2*b2\n",
    "\n",
    "    eq1 = np.sqrt(2)/np.sqrt(-B1 - np.sqrt(B1**2 - 4*C1))\n",
    "    eq2 = np.sqrt(2)/np.sqrt(-B2 + np.sqrt(B2**2 - 4*C2))\n",
    "    \n",
    "    return abs(eq1 - eq2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 602,
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds_theta = (0*np.pi/180, 90*np.pi/180)\n",
    "bounds = [bounds_theta]\n",
    "theta = 20*np.pi/180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 603,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     fun: 0.0016750781558880767\n",
      "     jac: array([0.])\n",
      " message: 'Optimization terminated successfully.'\n",
      "    nfev: 15\n",
      "     nit: 5\n",
      "    njev: 5\n",
      "  status: 0\n",
      " success: True\n",
      "       x: array([0.])\n"
     ]
    }
   ],
   "source": [
    "result = minimize(objective, phi, method = 'SLSQP', bounds = bounds)\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 604,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.]\n"
     ]
    }
   ],
   "source": [
    "print(result.x *180/np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 896,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxcdb3/8dcnS9Mm6ZalC13IwtYWCi2ltHVDQVpBRBQV2cqi4HJVLooXuVcFXHFBf14VBREqyOJFlEVWWZWthJZSukE3utDSNOmSJk2zfX5/zAnEmKZpkjNnZs77+XjMI7OczPkMQ99z8pnv+X7N3RERkfjIiroAERFJLgW/iEjMKPhFRGJGwS8iEjMKfhGRmFHwi4jEjIJfMpaZXWlmt0Zdh0iqUfCLhMjMnjGz80N8/mFmNs/MtgSXKzs9XmZmT5hZg5ktN7MTOjw20MzuNrNNZvbfYdUoqUfBL9LPzGyAmX3LzN4AZgG/DT4Ajg9hdz8D8oEyYDpwTqcPmtuBhUAx8N/AXWZWGjx2DlAV/O5MM6sMoT5JQQp+STozu8zM/tzpvv81s5+b2Uwz29Xh0mhma/fyPGVm5mZ2kZm9GRy5frWb/f6fmW02sx1m9rSZTerw2M1m9isz+5uZ1ZnZCx2D0MwOM7NHzazWzFaY2Se7eYk3AmcCnwCeAb4C3Ao8aGYn9uS/0X44BfiRuze4+9pg3xcENR8CTAW+7e673f3PwGLg48HvZne4ZAU/JQYU/BKFW4E5ZjYMwMxygE8Bt7j7c+5e6O6FwHDgeRJHrd15P3AwcCJwecd2RicPBtuNABYAf+z0+KeBq4L9rgS+F9RXADwK3Bb87qeBX3f84GhnZmOAs4EL3X1+cPdud78O+APwta4KM7MzzWx7N5fx3bx+63T98OD6JGC1u9d1eHxRcD/ALcAMYC1Q5e6vdbMPySAKfkk6d98EPE3iiBhgDrDV3V/qtOkvgHoSLYruXOXu9e6+GLiJRDB3td/fu3udu+8BrgSONLOhHTa5293nu3sLiQ+Fo4L7Pwysdfeb3L3F3RcAfwZO72I3BwY/X+7isYUk2ipd1Xabuw/r5rJuL6/9IRIfdoPN7CASR/v5wWOFwI5O2+8ABgf7rHf3k919hLt/ay/PLxlIwS9RmUfiyJjg5y0dHzSzi4HjgDPdvS24r2MLqOMR8PoO198ADui8MzPLNrMfmtkqM9tJ4igXoKTDZps7XG8gEZyQCPNjOx6BA2cBo7p4Xe0BfUQXj03usN/+8mVgN/A6cA+Jv442BI/tAoZ02n4IUIfEmoJfovJXYLKZHU7iiPrttouZvQf4DnCqu799xNreAgouHY+Ax3W4Ph54s4v9nQmcCpwADOWdI2/rYtvO1gNPdToCL3T3z3fe0N03kGgJ/c7Mjgzuzgu+cD0P+FFXOzCzszp9sHW+dNnqcfdadz/L3Ue5+yQS/6bbW0xLgAozG9zhV44M7pcYU/BLJNy9EbiLREjObw9yMxsH3Amcux8952+aWX7Qcz8/+P3OBgN7gBoSrZDv70e59wOHmNk5ZpYbXI4xswl72f584E/AfSRG9fwCuBCY7e5/7+oX3P2PnT7YOl+6bPWYWaWZFQd/0XwIuAj4bvCcr5FoOX07GLp5Gom/Ov7c1XNJfCj4JUrzSLREOrZ5jifRQrmrw9Huvo5QnyLxZexjwE/c/ZEutvkDiTbQRmApiS+NeyT4cvRE4AwSf01sBq4B8vayfZO7X+3u44HngM+5+7vd/cme7nM/HE1ipE4d8APgLHfv+N/rDGAasA34IXC6u1eHUIekEdNCLBKVoH2xHBjl7jt78ftlwBogN/hCVkR6QEf8EgkzywIuBe7oTeiLSO/lRF2AxE8wLv4tEq2XORGXIxI7avWIiMSMWj0iIjGTFq2ekpISLysri7oMEZG08tJLL21199LO96dF8JeVlVFVVRV1GSIiaSWYIfbfqNUjIhIzCn4RkZhR8IuIxIyCX0QkZhT8IiIxo+AXEYkZBb+ISMxkdPD/deFGbn2+y2GsIiIpbcXmOq599DWq6/b0+3NndPA/+OombnpmTdRliIjst8Ubd/CLx16noan/ZxzP6OAvLylkXW0DrW2aiE5E0kttfeJIv6hgQL8/d0YHf0VJAc2tzsZtu6MuRURkv9TUNzEgO4vCvP6fWSejg7+spACA1Vt3RVyJiMj+qd3VRFHBAMys3587o4O/PAj+tVvrI65ERGT/1NY3hdLmgQwP/pLCAQzOy2GNgl9E0kxNfRPFhQr+/WZmlJcWsFrBLyJppra+ieH5Cv5eKSsu0BG/iKQdtXr6oLykgI3bd7OnpTXqUkREemRPSyu79rRQrODvnYrSAtxhXU1D1KWIiPRIbX0TAEXq8fdO+dtDOtXuEZH0ULMrEfw64u+l9rH86vOLSLp4+4i/IC+U58/44B8yMJeSwgEayy8iaeOd4NcRf6+Vl2hIp4ikj5p6tXr6rLxEQzpFJH3U1u8hO8sYOig3lOePRfCXlRRQXbeHusbmqEsREdmnxMlbuWRl9f88PRCT4K8IvuB9Q0M6RSQNhHnyFsQk+MtLCgEN6RSR9KDg7wcHFudjBmuqFfwikvpq6psoDmkoJ4Qc/Gb2n2a2xMxeNbPbzWygmZWb2Qtm9rqZ3Wlm4X2sBQbmZnPA0EGs0bz8IpIG0vaI38zGAF8Gprn74UA2cAZwDfAzdz8Y2AZcGFYNHZWXFLBGPX4RSXEtrW1sb2hOz+AP5ACDzCwHyAc2AR8A7goenwd8NOQagCD4q3fhrvV3RSR1bWtIjD4May5+CDH43X0j8BNgHYnA3wG8BGx39/Zl4zcAY7r6fTO7yMyqzKyqurq6z/WUlxSws7Hl7TPiRERSUdhn7UK4rZ7hwKlAOXAAUAB8qItNuzwEd/fr3X2au08rLS3tcz3lmrNHRNJATf0eIE2DHzgBWOPu1e7eDNwNzAKGBa0fgLHAmyHW8DYFv4ikg9q3p2tIz1E964AZZpZviWXijweWAk8ApwfbzAXuCbGGt40dPoicLFPwi0hKS+tWj7u/QOJL3AXA4mBf1wP/BVxqZiuBYuDGsGroKCc7i/HF+Qp+EUlp7XPxD88PZ54eSIy6CY27fxv4dqe7VwPTw9zv3pRr/V0RSXG19U0My88lJzu8hkwsztxtV15SwNqaetraNKRTRFJTbX0TRfnhntcar+AvLaCxuY3NOxujLkVEpEs19XtC7e9D3IJfI3tEJMWFPV0DxDT4NUuniKSq2vqmUM/ahZgF/8jBAxmUm631d0UkJbW1OdtCnqcHYhb8WVlGmZZhFJEUtWN3M61tTlGIJ29BzIIfEqtxKfhFJBWFvch6u9gFf1lJPutqG2hubYu6FBGRf7GtIfyzdiGGwV9eUkhrm7Nh2+6oSxER+RftZ+0q+PvZO0M6tRqXiKSWtydo06ie/lXRPqRT6++KSIqpTcKUzBDD4B9eMIChg3L1Ba+IpJya+iYK83LIy8kOdT+xC354Z84eEZFUkoyzdiGmwV9RUsAatXpEJMUo+ENUXlLAmzsa2d3UGnUpIiJvq9nVFPoYfohp8JcFX/Cq3SMiqURH/CFqH9KpOXtEJFW4eyL4Qx7KCTEPfs3SKSKpYteeFppa29TqCUtBXg4jh+RpSKeIpIx3FlkPd4I2iGnwA5Rp/V0RSSHJmqANYhz8FaUF6vGLSMqoTdI8PRDj4C8vKaCmvokdDc1RlyIi0qHVo+APTXlJIQBrNKRTRFJATZImaINYB38+oFk6RSQ11NbvIS8ni0G54c7TAzEO/nFF+WQZrNnaEHUpIiLU1CfO2jWz0PcV2+DPy8lm7PB8jewRkZSQrJO3IMbBD4kveNXqEZFUsK2+KSlj+EHBz5rqetw96lJEJObaWz3JEPvgr29qpXrXnqhLEZGYS9YEbaDgB9Dc/CISqcbmVhqaWhX8yfDOwusKfhGJTjKna4CYB/8BwwYxIDtLwS8ikUrmdA0Q8+DPzjIOLNaQThGJVk194nvGZJy1CzEPfmgf0qngF5HoJHNKZlDwU15awBs1DbS2aUiniEQjmRO0QcjBb2bDzOwuM1tuZsvMbKaZFZnZo2b2evBzeJg17Et5cQFNrW28uX13lGWISIzV1DeRm20MGZiTlP2FfcT//4CH3P0w4EhgGXA58Ji7Hww8FtyOjEb2iEjUanc1MTw/OfP0QIjBb2ZDgPcCNwK4e5O7bwdOBeYFm80DPhpWDT1RXqrgF5Fo1STx5C0I94i/AqgGbjKzhWb2OzMrAEa6+yaA4OeIrn7ZzC4ysyozq6qurg6tyNLCPArzchT8IhKZ2vo9SRvRA+EGfw4wFbjO3acA9exHW8fdr3f3ae4+rbS0NKwaMTPKSvJZreAXkYjUJnGCNgg3+DcAG9z9heD2XSQ+CN4ys9EAwc8tIdbQI+UlhVp/V0Qik8wJ2iDE4Hf3zcB6Mzs0uOt4YClwLzA3uG8ucE9YNfRUeUkBG7Y1sKelNepSRCRmmlraqGtsSWqPP+yxQ18C/mhmA4DVwPkkPmz+ZGYXAuuAT4Rcwz5VlBTQ5rC+toGDRgyOuhwRiZFtDckdww8hB7+7vwxM6+Kh48Pc7/4qC4Z0rq6uV/CLSFLV7EruBG2gM3eBxElcAGtr1OcXkeRK9lm7oOAHYGh+LsUFAzSkU0SSrjaCVo+CP1BeUsBqLcgiIkn20tpaAEoHZ8ZwzrRSplk6RSTJHlv2FvOee4Ozjh3PsHwd8SddeUkBW+r2UL+nJepSRCQGNmxr4NI/LWLSAUP45ocnJnXfCv5AhSZrE5EkaWpp4z9uW0hbm/Prs6YyMDc7qftX8Ac0WZuIJMsPH1zOy+u386PTJ3NgMKowmRT8gQOLgiGdCn4RCdFDr27m98+s4bxZZXzoiNGR1KDgDwwakM0BQwfqiF9EQrOupoHL7lrEkWOHcsVJEyKrQ8HfQXlpgWbpFJFQ7Glp5Yu3LcCAX545lQE50cWvgr+DxFj+Xbhr/V0R6T/uzpX3LmXxxh389JNHMa4oP9J6FPwdTB47jJ2NLby+ZVfUpYhIBrnl+Te4ff46vnBcJR+cODLqcroP/mBx9F+Z2StmVm1m68zsATP7opkNTVaRyTKrshiAZ1dujbgSEckUz67cylX3LeWECSP42omH7vsXkmCvwW9mDwKfAR4G5gCjgYnA/wADgXvM7CPJKDJZxg7PZ3xRPs+uqom6FBHJAG/U1POF2xZQWVrAzz51FFlZyVlMfV+6m5b5HHfvfOi7C1gQXH5qZiWhVRaRWZXFPLB4E61tTnaKvEkikn7qGpv5zLwqAH537jEMHpgbcUXv2OsRfxeh36tt0s3MymJ2Nraw5M0dUZciImmqtc255I6XWb21nl+fNZXxxdF+mdvZPr/cNbM6M9vZ6bLezP5iZhXJKDKZZrb3+dXuEZFe+ukjK3hs+RauPGUisypTrzHSk1E91wKXAWOAscDXgBuAO4Dfh1daNEYMHsjBIwoV/CLSK/e8vJFfP7mKM48dz9kzDoy6nC71JPjnuPtv3b3O3Xe6+/XASe5+JzA85PoiMauymBfX1NLU0hZ1KSKSRhat387X73qF6eVFXHnKJMxS83vCngR/m5l90syygssnOzyWkWc6zawsYXdzK4s2bI+6FBFJE1t2NnLRLVWUDs7jurOiPTN3X3pS2VnAOcAW4K3g+tlmNgj4jxBri8yMiiLM4NmVaveIyL41Nrfy2Vteoq6xhRvOnUZxYfJW0+qNfQa/u69291PcvcTdS4PrK919t7v/MxlFJtuw/AEcfsBQnl2VcYOWRKSfuTtX3L2YReu3c+0nj2LC6CFRl7RPPRnVc4iZPWZmrwa3J5vZ/4RfWrRmVRazcN12dje1Rl2KiKSwG/6xmrsXbuTSDx7CnMNHRV1Oj/Sk1XMD8A2gGcDdXwHOCLOoVDCzspim1jZeemNb1KWISIp6YvkWfvDgck4+YjRf+sBBUZfTYz0J/nx3n9/pvoxfmPaYsiJyskztHhHp0sotdXz59oVMGDWEH39icsqO4OlKT4J/q5lVEozgMbPTgU2hVpUCCvJyOGrcMI3nF5F/s3XXHi64uYq83CxumDuN/AHdzX6TenoS/F8EfgscZmYbgUuAz4daVYqYVVnMKxu2s7OxOepSRCRFNDS1cOG8KrbUNXLDudMYM2xQ1CXtt56O6jkBKAUOc/d3u/va0CtLATMrS2hzmL+6NupSRCQFtLY5X779ZV7ZsJ1fnDGFKePT8xzWvf59YmaX7uV+ANz92pBqShlTxg8jLyeLZ1fVcEIKLJ4gItFxd666bwl/X/YWV31kEidOSo8RPF3prjE1OPh5KHAMcG9w+xTg6TCLShUDc7OZVjZcX/CKCDf8YzV/eO4NLnpvBXNnlUVdTp90Ny3zVe5+FVACTHX3r7r7V4GjSUzWFguzKktYvrmOml17oi5FRCJy36I3+f4Dyzl58mgun3NY1OX0WU++3B0PNHW43QSUhVJNCmqfpvl59flFYmn+mlq++qdFHFM2nJ9+4siUWUWrL3oyBukWYL6Z/YXEkM7TgHmhVpVCJo8ZSmFeDs+u2srJk0dHXY6IJNHKLXV89g9VjC0axA3nTmNgbnbUJfWLfQa/u38vWH/3PcFd57v7wnDLSh052VlMLy/iOY3nF4mVLXWNnHfTi+RmG/POn86w/AFRl9RvuhvVU+juuwDcvX2d3b1uk8lmVRbz+PItbNqxm9FD02/Mrojsn/o9LVx4cxU1u5q48+IZjCtKraUT+6q7Hv89ZvZTM3uvmRW032lmFWZ2oZk9DMzZ1w7MLNvMFprZ/cHtcjN7wcxeN7M7zSzlP0bb+/w66hfJfC2tbXzp9oUseXMHvzxzCpPHDou6pH7X3aie44HHgIuBJcFauzXArcAoYK6739WDfXwFWNbh9jXAz9z9YGAbcGFvi0+WCaOGMCw/V9M3iGQ4d+db9y7h8eVbuPrUwzl+Qmaev9Ntj9/dHwAe6O2Tm9lY4GTge8Clljj76wPAmcEm84Arget6u49kyMoyZlYU89yqGtw9rSZjEpGeu+6pVdz2wjo+f1xlyq6X2x/CXhvs58DXgfbFa4uB7e7ePrvnBhKLuP8bM7vIzKrMrKq6ujrkMvdtVmUxG7fvZl1tQ9SliEgI7nl5Iz96aAUfOfIALjvx0KjLCVVowW9mHwa2uPtLHe/uYtMu1+119+vdfZq7TystLQ2lxv0xs7IEgGe0HKNIxnlm5Va+9n+LOLa8iB9/YnJGjNXvTphH/O8CPmJma4E7SLR4fg4MM7P2FtNY4M0Qa+g3laUFjBicp+kbRDLM0jd3cvEtL1FeUsD1504jLyczxup3p9vgN7N8M/u4mf2nmV1iZh8zsx6Na3L3b7j7WHcvI7Fi1+PufhbwBHB6sNlc4J4+1J80Zsa7Dip5u88vIulvfW0D5900n8EDc5h3wXSGDsqNuqSk2GvwBwuuPA1MBy4lMWnbHGChmZ3Wh33+F4kveleS6Pnf2IfnSqqZlcXU1Dfx2lsZf+qCSMbbVt/E3Jvm09jcyrwLpsfqHJ3uRvVcBUx393ozO9HdvwNgZiXAM8BferoTd38SeDK4vprEh0namRWM53921VYOHTV4H1uLSKpqbG7lM3+oYsO23dxywXQOGRmvf8/dtXqMd0bjdOxttIZXTmobOzyf8UX5Gs8vksZa25wv3b6QBeu28fNPHcWxFcVRl5R03R3xXwk8HczTM8LMriDxZexxwDfCLy01zaos5m+LN9Ha5mRn+Df/IpnG3fnWPa/y6NK3uPKUiZx0RDwnXuzuzN0/Ae8HXgGuBfYAjwPHuPvdySkv9cysLKausYUlb+6IuhQR2U+/emIlf3xhHZ97XyXnvas86nIis68zd3cBPZmWITZmvt3nr8nIOTxEMtWfqtbzk0de47QpY/j67Mw+QWtfwj5zN+OMGDyQg0cUqs8vkkaeWLGFb9y9mPccXMI1H8/8E7T2RcHfC7Mqi3lxTS1NLW373lhEIrVo/Xa+cOsCDhs1mOvOPpoBOYo9/RfohZmVJexubmXRhu1RlyIi3Vi7tZ4Lbn6R4sIB3HT+MRTm9WTRwcyn4O+FGRVFmMGzmrdHJGVt3bWHuTfNp82deRdMZ8TggVGXlDIU/L0wLH8Akw4YwjOat0ckJTU0tXDBzS/y1s5GbjzvGCpLC6MuKaUo+HtpVmUJC9dtY3dTbM9nE0lJrW3Ol29fyKsbd/C/n57K1PHDoy4p5Sj4e2lmZTHNrU7VG7VRlyIiHXzn/qX8fdkWvn3KJD44MTNX0OorBX8vHVNWRE6WaVinSAr5/T/XcPOza7nw3eXMnVUWdTkpS8HfS4V5ORw5bpiCXyRFPLJkM9/521JmTxrJFSdNiLqclKbg74NZlcUs3rCdnY3NUZciEmuL1m/ny3csZPKYofz8U1M0j9Y+KPj7YGZlMW0OVWvV5xeJyvraBi6cV0VJYR6/m3sMgwZk/gpafaXg74Op44czIDuL51cr+EWisGN3Mxfc/CJ7Wlq56bxjKB2cF3VJaUHB3wcDc7M5avwwnl+tPr9IsjW1tPH5W19ibU09vz37aA6O2WIqfaHg76MZFcW8unGH+vwiSeTuXPGXxTy7qoYffGwysw4qibqktKLg76MZFUXq84sk2S8fX8ldL23gK8cfzOlHj426nLSj4O8j9flFkuuvCzfy00df42NTxnDJCQdHXU5aUvD30cDcbI4apz6/SDI8v7qGr9/1CseWF/GDjx+BmYZt9oaCvx/MqChSn18kZKuqd3HxLS8xrmgQ158zjbwcDdvsLQV/P5hRofH8ImGqrW/i/JteJCfLuOm86QzNz426pLSm4O8HU9TnFwlNa5vzlTsWsnlHIzfMncb44vyoS0p7Wo6mHwwaoD6/SFh+8djr/OP1rXz/tCM0xXI/0RF/P1GfX6T/PbliC794/HU+PnUsn54+LupyMoaCv5+ozy/SvzZsa+CSO1/m0JGD+e5HD9cInn6k4O8n6vOL9J89La184Y8LaG11rjv7aE281s/U4+8n6vOL9J/v3L+UVzbs4LfnHE15SUHU5WQcHfH3I/X5Rfru7gUbuPX5dVz83gpmTxoVdTkZScHfj9TnF+mb5Zt3csVfFjO9vIjLZh8adTkZS8Hfj9r7/C+ozy+y3+oam/n8rQsYPDCXX545hZxsxVNY1OPvR+rzi/SOu3PZ/73CutoGbv/sDEYMHhh1SRlNH6n9bEZFEYs37qBOfX6RHrvxn2t4aMlmLp9zGNPLi6IuJ+Mp+PvZO33+bVGXIpIW5q+p5QcPLmfOpFF85j3lUZcTC6EFv5mNM7MnzGyZmS0xs68E9xeZ2aNm9nrwM6POwX5nPL/aPSL7sqWukS/etoDxRfn86BOTdZJWkoR5xN8CfNXdJwAzgC+a2UTgcuAxdz8YeCy4nTHU5xfpmZbWNr5020LqGpu57uypDBmoGTeTJbTgd/dN7r4guF4HLAPGAKcC84LN5gEfDauGqKjPL7JvP35kBS+sqeX7px3BYaOGRF1OrCSlx29mZcAU4AVgpLtvgsSHAzBiL79zkZlVmVlVdXV1MsrsN8eqzy/Sradeq+a3T63mrGPH87GpWjM32UIPfjMrBP4MXOLuO3v6e+5+vbtPc/dppaWl4RUYgqnjh5ObbWr3iHShubWNq+9bQkVJAd86ZWLU5cRSqMFvZrkkQv+P7n53cPdbZjY6eHw0sCXMGqKgPr/I3t0+fx2rquu54qQJWj4xImGO6jHgRmCZu1/b4aF7gbnB9bnAPWHVEKUZFcXq84t0smN3Mz979DVmVRZz/IQuu7ySBGEe8b8LOAf4gJm9HFxOAn4IfNDMXgc+GNzOOBrPL/LvfvXESrbvbua/T56goZsRCm3KBnf/J7C3d/b4sPabKjr2+d9/mI5sRN6oqefmZ9byiaPHMumAoVGXE2s6czck6vOL/KtrHlpOTrbx1RM162bUFPwhUp9fJOHFtbU8sHgzn3tfJSOHaAK2qCn4Q6Q+vwi0tTnfvX8po4YM5LPvqYi6HEHBH6q3+/xr1O6R+Lp30Zss2rCDy2YfqrVzU4SCP0Tv9Pm1MIvE0+6mVq55aDmHjxnCaVPGRF2OBBT8IZtRUcyr6vNLTN34z9Vs2tHI/5w8kawsDd9MFQr+kM2oKKa1zal6Q31+iZctdY38+slVzJ40khkVxVGXIx0o+EOmeXskrq595DWaW9u4/EMToi5FOlHwh0x9fomjZZt2cmfVes6dWUZ5SUHU5UgnCv4kUJ9f4sTd+d7fljF0UC5f/sDBUZcjXVDwJ4H6/BInT66o5p8rt/KV4w9maL5W1UpFCv4kUJ9f4qK5tY3v/m0pFSUFnD3jwKjLkb1Q8CeB+vwSF3cEc+1/46QJ5GYrXlKV3pkkUZ9fMt2O3c387O+vM7OimBM0135KU/AnybHl6vNLZvv1EyvZ1tCkufbTgII/SaYeOEx9fslY62oauOmZtZw+dSyHj9Fc+6lOwZ8k+QNyOHKs+vySmX740DKys4yvzdZc++lAwZ9E6vNLJpq/JjHX/ueP01z76ULBn0Tt4/lfXKujfskMbW3O1fcvYfRQzbWfThT8STStbDgFA7J5dOlbUZci0i/+vGADr27cyeUfOkxz7acRBX8SDczN5rjDRvDo0rdobfOoyxHpk/o9Lfzo4RVMGT+Mjxx5QNTlyH5Q8CfZ7Emj2LqriQXrNKxT0ttvnlpFdd0evvnhiRq+mWYU/En2/kNLGZCdxcOvbo66FJFe27CtgeufXs2pRx3A1PHDoy5H9pOCP8kGD8xl1kHFPLRkM+5q90h6uuahFZjBf805LOpSpBcU/BGYM2kUG7btZummnVGXIrLfXnqjlvsWvclF76nggGGDoi5HekHBH4ETJo4ky+DhJRrdI+klMXxzGSOH5HHx+yqjLkd6ScEfgZLCPKYdWMQjS9Tnl/Ryz6KNLFq/na/PPoyCvJyoy5FeUvBH5MRJI1m+uY61W+ujLkWkRxqaWrjmwRVMHjuU06aMiboc6QMFf0RmTxoFwMM66pc0cf3Tq9m8s5FvfngiWVkavpnOFPwRGVeUz8TRQxT8khY27djNb55axQXoNYMAAAbfSURBVMmTR3NMWVHU5UgfKfgjNOfwUSxYt50tOxujLkWkWz9+aAVtDpdr+GZGUPBHqL3d84jm7pEU9vL67dy9cCOfeXc544ryoy5H+oGCP0KHjCykrDhf7R5JWe7O1fctoaQwjy+8/6Coy5F+ouCPkJkxe9IonltVw44GzdEvqee+VzaxYN12Lpt9CIUavpkxFPwRO3HSKFranMdXqN0jqaWxuZVrHlzOxNFDOP3ocVGXI/0okuA3szlmtsLMVprZ5VHUkCqmjBvGiMF5PPyqgl9Sy+/+sZqN23fzzQ9PJFvDNzNK0oPfzLKBXwEfAiYCnzazicmuI1VkZRknThrJU69V09jcGnU5IgC8tbORXz+5itmTRjKzsjjqcqSfRXHEPx1Y6e6r3b0JuAM4NYI6UsbsSaPY3dzK069VR12KCAA/eXgFza1tXHHShKhLkRBEEfxjgPUdbm8I7vsXZnaRmVWZWVV1dWYH4oyKYt5/aCkDc7V0nUTP3SkrKeALxx3EgcUFUZcjIYjia/qumoX/NjG9u18PXA8wbdq0jJ64Pjc7i5vOnx51GSJAYrTZFzV0M6NFccS/Aeg4RGAs8GYEdYiIxFIUwf8icLCZlZvZAOAM4N4I6hARiaWkt3rcvcXM/gN4GMgGfu/uS5Jdh4hIXEVyKp67PwA8EMW+RUTiTmfuiojEjIJfRCRmFPwiIjGj4BcRiRlzT/1zo8ysGngj6jr6UQmwNeoikihurxf0muMi1V/zge5e2vnOtAj+TGNmVe4+Leo6kiVurxf0muMiXV+zWj0iIjGj4BcRiRkFfzSuj7qAJIvb6wW95rhIy9esHr+ISMzoiF9EJGYU/CIiMaPgD4mZ/d7MtpjZq3t53MzsF8GC86+Y2dRk19jfevCajzOzHWb2cnD5VrJr7E9mNs7MnjCzZWa2xMy+0sU2GfU+9/A1Z9r7PNDM5pvZouA1X9XFNnlmdmfwPr9gZmXJr3Q/uLsuIVyA9wJTgVf38vhJwIMkViSbAbwQdc1JeM3HAfdHXWc/vt7RwNTg+mDgNWBiJr/PPXzNmfY+G1AYXM8FXgBmdNrmC8BvgutnAHdGXXd3Fx3xh8TdnwZqu9nkVOAPnvA8MMzMRienunD04DVnFHff5O4Lgut1wDL+ff3ojHqfe/iaM0rw3u0KbuYGl86jYk4F5gXX7wKON7OulplNCQr+6PRo0fkMNDP4k/lBM5sUdTH9JfjTfgqJo8GOMvZ97uY1Q4a9z2aWbWYvA1uAR919r++zu7cAO4Di5FbZcwr+6PRo0fkMs4DE3CFHAv8L/DXievqFmRUCfwYucfednR/u4lfS/n3ex2vOuPfZ3Vvd/SgSa4RPN7PDO22SVu+zgj86sVt03t13tv/J7IlV2HLNrCTisvrEzHJJBOAf3f3uLjbJuPd5X685E9/ndu6+HXgSmNPpobffZzPLAYaSwm1PBX907gXODUZ9zAB2uPumqIsKk5mNau97mtl0Ev//1URbVe8Fr+VGYJm7X7uXzTLqfe7Ja87A97nUzIYF1wcBJwDLO212LzA3uH468LgH3/SmokjW3I0DM7udxOiGEjPbAHybxJdCuPtvSKw5fBKwEmgAzo+m0v7Tg9d8OvB5M2sBdgNnpPI/jh54F3AOsDjo/wJcAYyHjH2fe/KaM+19Hg3MM7NsEh9if3L3+83saqDK3e8l8WF4i5mtJHGkf0Z05e6bpmwQEYkZtXpERGJGwS8iEjMKfhGRmFHwi4jEjIJfRCRmFPwSe2b2uJk9EJyYtLdtBpnZU8GQvr7ur9TMHurr84j0loJfYs/dPwDsAU7uZrMLgLvdvbUf9lcNbDKzd/X1uUR6Q8EvkvAgcFY3j58F3ANvzzd/f/sDZvZLMzsvuL7WzL5vZs+ZWZWZTTWzh81slZl9rsPz/XUf+xMJjYJfJOEM4EQzG9L5ATMbAFS4+9oePtd6d58J/AO4mcSZrDOAqztsUwW8py8Fi/SWgl9iz8yOIDGp1m3Ax7vYpATYvh9PeW/wczGJhVfqgvZOY/ucLySm9z2glyWL9ImCXwQuAX4G3ErX7ZfdwMBO93Wchrfzl8J7gp9tHa63326fH2tg8LwiSafgl1gzs1LgfSSWynsGONDM/uVI3N23Adlm1jH8Dw9G+gwAZgL7O9rnEKDLtYlFwqbgl7i7GPiduzcHt2+n65kVHwHe3eH2NuAh4DngMeAbZlawH/t9P/C3/S9XpO80O6dID5jZFOBSdz/HzI4DvubuH+7D8z0NnBr8NSGSVDriF+kBd18IPNFfJ3AB1yr0JSo64hcRiRkd8YuIxIyCX0QkZhT8IiIxo+AXEYkZBb+ISMz8f4qnX+72I7cFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = [0.7,0.8,0.9,1,1.05,1.1,2.4,2.5,2.6,2.7,2.8,2.9,3,3.1,3.15,3.2]\n",
    "y = [90,52.41,30.93,8.74,0,-np.inf,0,9.74,19.6,26.78,33.19,39.4,45.72,52.47,56.14,90]\n",
    "plt.plot(x,y)\n",
    "plt.xlabel('λ (μm)')\n",
    "plt.ylabel('θ (deg)')\n",
    "plt.title('yz-plane Φ = 90°')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 884,
   "metadata": {},
   "outputs": [],
   "source": [
    "def objective(x):\n",
    "    theta = x[0] \n",
    "    \n",
    "    kx = np.sin(theta) * np.cos(0)\n",
    "    ky = np.sin(theta) * np.sin(0)\n",
    "    kz = np.cos(theta)\n",
    "    \n",
    "    a1 = nxsq**(-1)\n",
    "    b1 = nysq**(-1)\n",
    "    c1 = nzsq**(-1)\n",
    "\n",
    "    a2 = nx2wsq**(-1)\n",
    "    b2 = ny2wsq**(-1)\n",
    "    c2 = nz2wsq**(-1)\n",
    "\n",
    "    B1 = (-(kx**2)*(b1+c1) - (ky**2)*(a1+c1) - (kz**2)*(a1+b1))\n",
    "    C1 = (kx**2)*b1*c1 + (ky**2)*a1*c1 + (kz**2)*a1*b1\n",
    "\n",
    "    B2 = (-(kx**2)*(b2+c2) - (ky**2)*(a2+c2) - (kz**2)*(a2+b2))\n",
    "    C2 = (kx**2)*b2*c2 + (ky**2)*a2*c2 + (kz**2)*a2*b2\n",
    "\n",
    "    eq1 = np.sqrt(2)/np.sqrt(-B1 - np.sqrt(B1**2 - 4*C1))\n",
    "    eq2 = np.sqrt(2)/np.sqrt(-B2 + np.sqrt(B2**2 - 4*C2))\n",
    "    \n",
    "    return abs(eq1 - eq2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 885,
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds_theta = (0*np.pi/180, 40*np.pi/180)\n",
    "bounds = [bounds_theta]\n",
    "theta = 20*np.pi/180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 886,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     fun: 0.0006290130823924933\n",
      "     jac: array([0.])\n",
      " message: 'Optimization terminated successfully.'\n",
      "    nfev: 12\n",
      "     nit: 4\n",
      "    njev: 4\n",
      "  status: 0\n",
      " success: True\n",
      "       x: array([2.77555756e-17])\n"
     ]
    }
   ],
   "source": [
    "result = minimize(objective, phi, method = 'SLSQP', bounds = bounds)\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 887,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.59027734e-15]\n"
     ]
    }
   ],
   "source": [
    "print(result.x *180/np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 897,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd5xU9b3/8ddnK1vouyBFerOB4CIKFmxYYn4msUUNMcaaRBMjJhrvvbmmXcu1xcSoWGO9iS0aBRsKSBEEkSYs7NL77tLZZevn98cMuCILA+zMmdl5Px+PeezMOWfnvEHmM1+/53u+X3N3REQkeaQEHUBERGJLhV9EJMmo8IuIJBkVfhGRJKPCLyKSZFT4RUSSjAq/JAUzu9PMXgg6h0g8UOEXiREzm2xmV0Xx/c3M7jGzsvDjXjOzevuvMbM1ZvahmbWNVg6Jfyr8IlFkZhlm9lszWw4MBR4PfwGcEYXTXQd8BxgA9AfOB64P58gGrgX6AM8AN0fh/JIgVPglUGbW08w2mtmg8OuOZlZqZsPN7EQz217vsdPMljXwPt3MzM3sunCrdq2ZjdrHeV8xs3VmtsXMJprZUfX2PWtmj5jZO2a2zcymmVnPevv7mdkH4dyFZnbJPv6ITwGXAxcDk4FfAC8AY81sxIH8XUXgSuB+d1/l7quB+4EfhfelAKl7PCRJqfBLoNy9GLgNeDHcKn0GeNbdx7v7VHfPdfdcoDXwKfDyft7yNKA3MAK43czObOC4seHj2gGfAy/usf8y4Hfh8xYBfwIwsxzgA+Cl8O9eBvyt/hfHLmbWCfgBcLW7Tw9vrnD3R4HngFv3FszMLjezzft4dGngz3QUMLve69nhbbj7duDvQCHwY+ChBt5DkoAKvwTO3Z8AFgPTgA7Af+zlsIeBHQ3sq+937r7D3ecS+hK5rIFzPu3u29y9ErgTGGBmLesd8rq7T3f3GkJfCseGt58PLHP3Z9y9xt0/B14DLtrLabqGf36xl32zgG4NZHvJ3Vvt47GigT97LrCl3ustQO6ufn53/4u7H+buw919QwPvIUlAhV/ixRPA0cBfwsV4NzO7HhgOXO7udeFt9buA6reAV9Z7vhzouOeJzCzVzO42s2Iz2wosC+/Kq3fYunrPywkVVQgV8yH1W+DAFcBhe/kz7SrQx+xlX/96520s24EW9V63ALa7ZmKUPajwS+DMLJdQ18NTwJ1m1qbevpOBPwAXuPvu1uyuLqDwo34L+PB6z7sAa/ZyysuBC4AzgZZ81fK2vRy7p5XAhD1a4Lnu/pM9D3T3VYS6hJ40swHhzZnhkT0/Au7d2wnM7Io9vtj2fDTU1TOf0IXdXQaEt4l8jQq/xIM/AzPd/RrgHeAxADM7HPgH8EN3XxThe/2XmWWH+9yvCv/+npoDlUAZkA38zwFkfRvoY2YjzSw9/BhsZkc0cPxVwD+BfxMa1fMwcDVwtrt/uLdfcPcX9/hi2/PRUFfPc8AtZtbJzDoCo4BnD+DPJklChV8CZWYXAOcAN4Q33QIMMrMrgDMIdaG8Wq+1u78W7ARCF2PHAfe5+/t7OeY5Qt1Aq4EvCV00joi7byN04fj7hP5vYh1wD5DZwPFV7v57d+8CTAVucPeT3H18pOc8AI8T+oKZC8wj9CX6eBTOIwnO1P0nTYGZdQOWAunhC7Ii0gC1+EVEkowKv4hIklFXj4hIklGLX0QkyaQFHSASeXl53q1bt6BjiIgklJkzZ5a6e/6e2xOi8Hfr1o0ZM2YEHUNEJKGEZ4X9BnX1iIgkmagVfjM73Mw+NrMFZjbfzH4R3n6nma02sy/Cj/OilUFERL4pml09NcAod//czJoDM83sg/C+B939viieW0REGhC1wu/ua4G14efbzGwB0Cla5xMRkcjEpI8/fDv9QELzrQPcaGZzzOxpM2vdwO9cZ2YzzGxGSUlJLGKKiCSFqBf+8JS7rwE3u/tW4FGgJ6GFLdYSWh7uG9x9tLsXuHtBfv43RiOJiMhBimrhN7N0QkX/RXd/HcDd17t7bXhBjSeA46OZQUREvi5qffzh5d6eAha4+wP1tncI9/8DfJfQ9LEija62ztm+s4atO6tDj4rw84pqtu6sYfvOGtrkZtC5VRadW2fRqXUW2RkJcWuLyCGJ5r/yYcBIYK6Z7Vpz9A7gMjM7FnBCS89dH8UM0gS4O8vKylmzuSJctKvZtrNmdwHftW1XYd+1b1vlgc/O3CYng86tw18ErbLo3Dp795dC59bZ5Gbqi0ESXzRH9Uxi70vZjYnWOaXpKK+qYWpxGR8XbmB8YQmrNlV84xgzaJ6ZRousdFo0S6dFVhpd2mTTIiud5s3SwtvSadHs68e0aBZ6np2ZysYdVazaVMGqTeWs2lTB6s0VrNpUwcJ12xi3YAOVNXVfO2fLrPR6XwzZXz0PfzG0zEqP1V+RyEFT80XigrtTXLKD8YUbmLCohGlLNlJVW0d2RipDe+Zxw6k96d0uN1TAw8U8JyONlJRIlsltWPsWzWjfohnHdf3m4DJ3p3R71R5fCqHnS0p2MHFRKRXVtV/7nTY5Gfz89F788MRuh5xNJFpU+CUwu1r14wtLGL9oAys3hlr1vdrl8sMTuzK8bzsGd29NZlpqIPnMjPzmmeQ3z2Rgl71/Mez6P4ZdXwoTF5Vy57+/ZOy8dfzvRQPo0jY7gOQi+5YQ8/EXFBS4JmlLfO7OktIdoUJfuIFpSzdSVVNHVnoqw3q15dS+7RjeJ5/D2yRusXR3Xpmxij+8/SW17tx+bj9+MKSrWv8SCDOb6e4F39iuwi/RVFFVy9QlpeFiX8KKjeUA9MzPYXjfdgzvm8/x3dsE1qqPljWbK7jttTl8sriUE3q04X8vGpDQX2iSmFT4JWbWbqlg7Nx1jF9UwqdLyna36of2bMvwvvkM79suKYqgu/OPz1byx3cWUOfOb847giuO76LWv8SMCr9E3c7qWh6bUMyj44uprKmjR34Ow/t81apvlt60WvWRWr25gtvDrf+hPdtyz4X9k+KLT4Knwi9R4+58uGADv397Pis3VnB+/w6MGtGX7nk5QUeLG+7Oy9NX8qd3vgTgjm8dweXHdyF0n6NIdDRU+DWqRw7J0tId/O7f8xlfWELvdrm8dM0QhvbKCzpW3DEzLh/ShVP65HHba3P4jzfmMXbuOu6+8Bg6t1brX2JLLX45KOVVNTzycRFPTFxKRloKN5/ZmyuHdiM9VYu67Y+78+K0Fdw1ZgFmxh3nHcFlxx+u1r80OrX4pVG4O2PnreOPb3/Jmi07+d7ATtx+bj/atWgWdLSEYWb84ISunNonn9tem8Mdb8xl7Ly13H1hfzq1ygo6niQBtfglYkUbtnHnW18yqaiUfoc15w/fOZrB3doEHSuh1dU5L05bzl1jF5Jixn9+6wguHazWvzQOXdyVg7a9soaHxy3m6UlLyc5I5daz+3L58V1IU7dOo1m5sZxfvTqbT5ds5JQ++dxz4TF0aKnWvxwaFX45YO7OW7PX8D9jFrB+ayWXFhzOr87pS15uZtDRmqS6OueFacu5a8xC0lKM/zr/SC4u6KzWvxw09fHLAVm4biu/fXM+05du5JhOLXnsB8ftdb4aaTwpKcYPT+zGqX3y+dWrc/j1a3MYM28td31PrX9pXGrxy9dsqajmoQ8X8dzU5TRvlsavz+7HpYMPJ1V3m8ZUXZ3z3NRl3PNuIVW1deRkpJKdkUZWRipZ6alkZ6SSlRH6udft6V9tzw7vywofm52RStvcDC06kwTU4pd9qqtzXp+1mrvHLqBsRxVXDOnCqLP60jonI+hoSSklxfjRsO4M79uO1z5fxbadNZRX1VBRXUdFVQ3lVbVs21nDhq2VlFfXUFFVS3lVLRXVtUTSlsvJSOWmM3pz1bBuTW6eJNk/tfiFeau38Ns35/H5is0M7NKKP1xwNEd3ahl0LDkI7k5lTR3lVbWhL4rwF0LoS6Fm9/P356/nwwXr6dY2m99++0hO79c+6OgSBbq4K3v1l3GLeeDDRbTNyeC2c/px4aDOmkQsSUxYVMLv/j2fJSU7GN43n/86/0h65ucGHUsakQq/fMPjE4q5a+xCLji2I7+/4GgtG5iEqmrqeG7qMv784WJ21tRy1bDu3HR6L5o307+FpqChwq+B2Elq101D5/fvwAOXHKuin6Qy0lK45uQefHTrcL47sBOjJy7htPsm8MqMldTVxX+jUA6OCn8S+tes1fznv+ZxRr92PHjpsRqxI+Q3z+Teiwbw5s+GcXibLH716hy+++gUvli5OehoEgUq/Enm/fnrGPXKbE7o3pZHrhikSdXkawYc3orXbhjKA5cMYM3mCr7zyGRufWU2G7btDDqaNCJ96pPIpMWl3PjSLI7u1JInrixI2oVRZN9SUozvDerMx7cO54ZTe/LmF6s5/b4JjJ5YTFVNXdDxpBGo8CeJmcs3ce1zM+iRn8PfrxpMbqZu4ZB9y81M4/Zz+/H+L09lSPc2/M+YhZzz0EQ+LtwQdDQ5RCr8SWD+mi1c9cx02rfI5Lmrj6dVtm7Kksh1z8vhqR8N5pmrBgNw1TOf8eNnP2Np6Y6Ak8nBUuFv4opLtvPDp6aTm5nGC9cMoV1zzZsvB+e0vu149+ZTuOO8fkxfupERD07grrEL2F5ZE3Q0OUAq/E3Yqk3l/ODJaZjBC9cM0RJ/csgy0lK47pSefHTrqVxwbCcen7CE0+4bz2szV2n4ZwJR4W+iNmzdyRVPTmNHZQ3P/XgIPXRHpjSids2bcd/FA3jjp0Pp2CqLUa/M5vtPfMraLRVBR5MIqPA3QZvLqxj51HRKtlXy7I+P58iOLYKOJE3UwC6teeMnQ7n3wv7MX72F8/78CR8v1MXfeKfC38Rsr6zhyqens7RsB0/8sIBBmkNfoiwlxbhk8OH8+6aTOKxlFlc9+xl3jVlAda2GfsYrFf4mZGd1LVc/+xnz1mzlkcsHMaxXXtCRJIn0yM/ljZ8O5YohXXh84hIufXwqqzer6yceqfA3EVU1dfzkhZlMX7aRBy4ZwFlHappdib1m6an86bvH8NfLB7Jo/XbO+/MnfPDl+qBjyR6iVvjN7HAz+9jMFpjZfDP7RXh7GzP7wMwWh3+qL+IQ1dY5v/znF3xcWMKfvnMMFxzbKehIkuTO79+Rt286icPbZHHtczP4w9tf6q7fOBLNFn8NMMrdjwBOAH5mZkcCtwPj3L03MC78Wg5SXZ3zm9fn8M6ctdxxXj8uH9Il6EgiAHTLy+G1nwzlyhO78tSkpVz8+FRWbiwPOpYQxcLv7mvd/fPw823AAqATcAHw9/Bhfwe+E60MTZ2788d3FvDPGav4+em9uO6UnkFHEvmazLRUfnfB0Tx6xSCWlGznvIc/4d1564KOlfRi0sdvZt2AgcA0oL27r4XQlwPQroHfuc7MZpjZjJKSkljETDgPfbiYpycv5UdDu/HLs/oEHUekQece04F3bjqZHnk53PDCTO58az6VNbVBx0paUS/8ZpYLvAbc7O5bI/09dx/t7gXuXpCfnx+9gAnqyU+W8Odxi7n4uM789vwjMdOc+hLfurTN5pUbhvLjYd15dsoyLnp0KsvLNN9PEKJa+M0snVDRf9HdXw9vXm9mHcL7OwC62+MAvTx9BX98ZwHfOqYDd1/YX2vkSsLISEvht98+ktEjj2N52Q7Of3gS78xZG3SspBPNUT0GPAUscPcH6u16C7gy/PxK4M1oZWiK3pq9hjvemMvwvvlaPUsS1oijDmPML06mZ7tcfvbS5/zXv+axs1pdP7ESzRb/MGAkcLqZfRF+nAfcDZxlZouBs8KvJQJTiku55R9fMLhbGx694jgy0nQbhiSuzq2z+ef1J3Ltyd15/tPlfO9vUzTVc4yYe/zPqFdQUOAzZswIOkag6uqc8x7+hIrqWt6+6SSaN9Pi6NJ0jFuwnlGvzKa6po67LuzP/xvQMehITYKZzXT3gj23q8mYIMbMW8vCddu45aw+KvrS5JxxRHvG/Pxk+nVowc9fnsVvXp+rrp8oUuFPADW1dTzwwSL6tM/l/P5qCUnT1LFVFv933QnccGpPXp6+gu88MpklJduDjtUkqfAngDdmrWZJyQ5uOauvLuZKk5aemsLt5/bjmasGs2FbJRc9NpV5q7cEHavJUeGPc1U1dfx53GKO6dSSs4/SxGuSHE7r247XfzKUrPRULnviU2Yu3xR0pCZFhT/O/WPGSlZtqmDUiD66SUuSSre8HP5x/Qm0zclg5FPTmFpcFnSkJkOFP47trK7lrx8tZnC31pzaR3cvS/LZNeSzU6ssfvTMdK3u1UhU+OPYC58uZ/3WSm4d0VetfUla7Vo04x/Xn0ivdrlc9/wMxs7Vnb6HSoU/Tm2vrOFv44s5uXceQ3q0DTqOSKDa5GTw0rUncEynlvzspc95/fNVQUdKaCr8ceqZSUvZuKOKUSP6Bh1FJC60zErn+auHMKR7W0a9MpsXpy0POlLCUuGPQ1vKqxn9yRLOPKI9xx7eKug4InEjJzONZ64azPA++fzHG/N48pMlQUdKSCr8cWj0J8Vsr6xh1AjNsS+yp2bpqTw+soBzjz6MP76zgIfHLSYRpp6JJyr8caZ0eyXPTF7G+f07ckSHFkHHEYlLGWkp/OWygXxvYCce+GAR97xbqOJ/ANKCDiBf9+j4YnZW13Lzmb2DjiIS19JSU7jv4gFkZaTy2IRiKqpq+O9vH6X1KSKgwh9H1m6p4PlPl3PhoM70zM8NOo5I3EtJMf74naPJSk/lyUlLKa+q5e4L+2tqk/1Q4Y8jf/moCHfn52eotS8SKTPjP751BNmZaTw8bjEV1bU8eOmxpKeqJ7shKvxxYkVZOf/8bCWXHd+Fw9tkBx1HJKGYGbec1YfsjFTuHrswdNf75YNolp4adLS4pK/EOPHQuEWkphg3nt4r6CgiCeuGU3vyhwuO4sMFG7j2uRmUV9UEHSkuqfDHgaIN2/jXrNVcObQb7Vs0CzqOSEIbeWI37rt4AJOLSrny6els21kddKS4o8IfBx78YDFZ6anccGrPoKOINAkXHdeZhy8byKwVm7niyWlsLq8KOlJcUeEP2LzVW3hn7lquPqk7bXIygo4j0mSc378jj488joXrtvH90Z9Ssq0y6EhxQ4U/YA98sIiWWelcfXKPoKOINDlnHNGep68czPKyci59fCprt1QEHSkuqPAHaObyTXy0cAPXndKDlllaQF0kGk7qncdzVx9PybZKLn5sKis3lgcdKXAq/AG6//1C8nIzuGpYt6CjiDRpg7u14cVrh7BtZw1X//0zdlQm92gfFf6ATCkqZUpxGT8d3ovsDN1OIRJt/Tu34pHLB1G0YTu/enV2Us/to8IfAHfnvvcL6dCyGZcP6RJ0HJGkcVLvPG4/tx9j5q7j0QnFQccJjAp/AD4u3MDnKzZz0+m9dWehSIxde3IPvj2gI//7XiHjC5NzDV8V/hirq3Pue28RXdpkc3FB56DjiCQdM+OeC4+hb/vm/PzlWSwv2xF0pJhT4Y+xd+ev48u1W7n5zN6aREokINkZaYweWYCZcf3zM5NuagdVnhiqrXMe+GARvdvlcsGxnYKOI5LUurTN5i+XDWTR+m386tU5SXWxV4U/hv41azVFG7Zzy1l9NF+4SBw4pU8+vz6nH+/MWcvoicmzfq8Kf4xU19bx0LhFHNWxBWcfdVjQcUQk7PpTevCtYzpwz7sL+WRxSdBxYiJqhd/MnjazDWY2r962O81stZl9EX6cF63zx5t/zljJyo0V3Dqir5aGE4kjZsa9F/WnT/vm3PjSLFaUNf07e6PZ4n8WOGcv2x9092PDjzFRPH/c2Fldy1/GFXFc19YM75sfdBwR2UNOZhqPjzwOd+e655v+PP5RK/zuPhHYGK33TyQvTlvBuq07GTWiD2Zq7YvEo65tc3j4soEUrt/Gba/NbdIXe4Po47/RzOaEu4JaN3SQmV1nZjPMbEZJSeL2u+2orOFvHxcxrFdbhvbMCzqOiOzD8L7tuHVEX/49ew1PfrI06DhRE+vC/yjQEzgWWAvc39CB7j7a3QvcvSA/P3G7R56dsoyyHVWMGtE36CgiEoGfDu/JuUcfxl1jFzBpcWnQcaIipoXf3de7e6271wFPAMfH8vyxtqWimscnFHNGv3YM6tLg/9yISBwxM+67eAC92uVy08ufN8lpnGNa+M2sQ72X3wXmNXRsU/DkJ0vYurOGW0b0CTqKiByA0MXeAmrqnOufn0lFVW3QkRpVNIdzvgxMBfqa2Sozuxq418zmmtkc4DTgl9E6f9DKtlfy9KSlfOuYDhzVsWXQcUTkAHXPy+Hh7w9kwbqt/Ob1pnVnb9Qmgnf3y/ay+alonS/evD1nLTuqarnpjF5BRxGRg3Rav3aMOqsP972/iKM7teSaJrJEqu7cjZLxhRvo1jabfoe1CDqKiByCnw7vxdlHteeusQuZUtQ0Lvbus/Cb2Ylm9kh4+GWJma0wszFm9jMzU/9FA3ZW1zJ1SRnD+7YLOoqIHKKUFOP+S46le14ON748i1WbEv9ib4OF38zGAtcA7xG6A7cDcCTwn0Az4E0z+3+xCJlopi3dyM7qOk7VXboiTUJuZhqjRx5HdU0dN7wwk53ViX2xd18t/pHufrW7v+Xua9y9xt23u/vn7n6/uw8HpsQoZ0IZX7iBzLQUTuzRNugoItJIeuTn8tD3j2Xe6q3c8Xpi39nbYOF39/12ZkVyTDKaUFjCCT3aallFkSbmjCPa88sz+/D6rNU8O2VZ0HEO2n4v7prZNjPbusdjpZm9YWZN4xJ3I1pRVs6S0h2ajE2kibrp9F6cdWR7/vjOAj5dUhZ0nIMSyaieB4BfAZ2AzsCthO66/T/g6ehFS0zjF4UWb9aFXZGmKSXFeOCSAXRtm83PXvyc1Zsrgo50wCIp/Oe4++Puvs3dt7r7aOA8d/8HoHkI9jC+sISubbPpnpcTdBQRiZLmzdIZPbKAypo6bng+8S72RlL468zsEjNLCT8uqbcvca9uRMHO6lqmFJcyvI+6eUSaul7tcnnw0mOZu3oLj00oDjrOAYmk8F8BjAQ2AOvDz39gZlnAjVHMlnCmh4dxqptHJDmcdWR7zjnqMJ6YuITS7ZVBx4nYfgu/uy9x92+7e56754efF7l7hbtPikXIRPFx4QYy0lI4QcM4RZLGr87py86aOv76UVHQUSIWyaiePmY2btfauWbW38z+M/rREs+uYZxZGRrGKZIseubncklBZ16ctjxhpnCOpKvnCeA3QDWAu88Bvh/NUIlo9zBO9e+LJJ1fnNGHFDPuf78w6CgRiaTwZ7v79D22Ne2ViA/CrmGcp/VT/75IsjmsZTN+fFJ33py9hi/XbA06zn5FUvhLzawn4RE8ZnYRoWUTpR4N4xRJbjec2pMWzdK5972FQUfZr0gK/8+Ax4F+ZrYauBn4SVRTJRgN4xSRllnp/HR4T8YXljC1OL7v6I10VM+ZQD7Qz91PcvdlUU+WQDSMU0QArhzajQ4tm3H3uwvjehK3BlfgMrNbGtgOgLs/EKVMCWd8YYmGcYoIzdJT+eWZffj1a3N4b/46zjm6w/5/KQD7avE3Dz8KCHXtdAo/biA0L7+EjV+0QcM4RQSA7w3qRK92udz7XiE1tXVBx9mrfU3L/Dt3/x2QBwxy91HuPgo4jtBkbQKs3FjOkhIN4xSRkLTUFH59dl+WlOzglZmrgo6zV5Fc3O0CVNV7XQV0i0qaBDS+cNdsnCr8IhJy1pHtOa5rax76cBEVVfE3gVskhf95YLqZ3Wlm/w1MA/4e3ViJY3xhCV3aaBiniHzFzLjtnH6s31rJM1OWBh3nGyIZ1fMn4CpgE7AZuMrd74p2sEQQGsZZxvC++bsveouIABzfvQ1n9GvHo+OL2Vxetf9fiKF9Lbaeu+t5eJ3dP4cfs/Z2TDL6bNlGKqpr1c0jInv1q3P6sr2yhr+Nj69pm/fV4n/TzO43s1PMbHc/hpn1MLOrzew94JzoR4xfHy8MDeM8sUde0FFEJA71O6wF3xvYmWenLGNNHK3Uta9RPWcA44DrgfnhtXbLgBeAw4Ar3f3V2MSMT+MXbWBI9zYaxikiDfrlWb3B4aEPFwUdZbcGb+ACcPcxwJgYZUkou4Zx/mBI16CjiEgc69w6m5EnduWZyUu59uQe9G7fPOhIEY3qkb3QME4RidTPTutFTkYa974XH9M2q/AfJA3jFJFItcnJ4PpTe/DBl+uZuXxj0HFU+A+GhnGKyIH68UndyW+eyd1jg5/AbZ+F38yyzexCM/ulmd1sZt8zs+xYhYtXGsYpIgcqOyONX5zRm8+WbeKjhRsCzbKvcfwXAROB44FbCE3Ydg4wy8y+u783NrOnzWzDrrV6w9vamNkHZrY4/LP1of8RYm/XbJwaxikiB+LSwYfTPS+He98tpLYuuFb/vlr8vwNOdffbgFJ3/4O7XwcMA+6O4L2f5Zvj/G8Hxrl7b0JDRW8/8MjBG1+oYZwicuDSU1MYNaIPheu38cas1YHl2FfhN2DXnKL1v5oimnHI3ScCe17FuICv5vn5O/CdSN4rnqzcWE5xyQ4tuiIiB+W8ozvQv3NLHvxgETurg5nAbV+F/05gopn9HmhnZneY2d+AycBvDvJ87d19LUD4Z4PV08yuM7MZZjajpKTkIE/X+MYvCmVR/76IHIyUlNAEbqs3V/DCp8uDydDQDnf/J3AaMAd4AKgEPgIGu/vr0Q7m7qPdvcDdC/Lz46fITigs4fA2WfTQME4ROUjDeuVxcu88/vpxEVt3Vsf8/Psc1ePu2939VXd/wN3vDz/fcQjnW29mHQDCP4O9tH2Aamrr+HRJGSf31jBOETk0t53Tj83l1YyesCTm5471OP63gCvDz68E3ozx+Q/J7FVb2F5Zw0m9NJpHRA7N0Z1a8u0BHXlq0lI2bN0Z03NHrfCb2cvAVKCvma0ys6sJjQY6y8wWA2cR2eiguDG5qBQzOFGLqotIIxh1Vh+qa+v487jFMT3vPidpOxTuflkDu86I1jmjbXJRKUd1bEHrnIygo4hIE9AtL4fLh3ThxWkruObkHjGbAkZTNkSovKqGWSs2M6ynunlEpPHcdHpvMtNSuO6tWU8AAA0BSURBVO/92E3gpsIfoc+WbaKqto6h6t8XkUaU3zyTa07qzjtz1jJn1eaYnFOFP0JTikrJSE1hcLeEnGVCROLYtaf0oE1OBve8uzAm51Phj9CkolIGdmlFdkbULouISJJq3iydG0/rxeSiMj5ZHP0bVlX4I7BxRxVfrt2qYZwiEjVXnNCFzq2zuHvsQuqiPIGbCn8EphaX4Y7690UkajLTUhk1og/z12zl7blro3ouFf4ITC4uJTczjQGdWwYdRUSasAsGdKLfYc25//1Cqmrq9v8LB0mFPwJTiko5oUcb0lL11yUi0ZOSYtxwak+Wl5WzaP226J0nau/cRKzaVM6ysnKGavy+iMRAXm4mQFSnbFbh348pRWUAnNRbhV9Eoi8jLVSWK9XVE5xJRaXkN8+kd7vcoKOISBLI3F341eIPhLszpbiUYT3bahpmEYmJzPRw4a9Wiz8Qheu3Ubq9SsM4RSRmMtNCa3mrqycgk8P9+8NU+EUkRtTVE7ApRaV0z8uhU6usoKOISJLI1MXd4FSHl1kc2lOLrohI7GSmh7t61Mcfe3NWbWZHVa3m5xGRmFJXT4AmLS4LLbOoFr+IxFBaipFiaMqGIEwuDi2z2CpbyyyKSOyYGRlpKerjj7XQMoubNJpHRAKRmZaqwh9r05dupLrWtb6uiAQiMy1FffyxNqW4LLzMYpugo4hIEspMT9GonlibtLiUQV1bkZWRGnQUEUlC6uqJsV3LLKqbR0SCoq6eGJtSXArAME3DLCIBydSontiaXFRG88w0+nfSMosiEozMtFT18cfSlOJShvRoq2UWRSQwGerqiZ2VG8tZXlbOsF66W1dEgqOunhja3b+vG7dEJECZ6amasiFWJhWVaZlFEQmcWvwxUlfnTCnSMosiErxoD+dMi9o774OZLQO2AbVAjbsXBJGjvsL12yjbUaVuHhEJXLRH9QRS+MNOc/fSAM//NZOL1L8vIvEhM11dPTExpbiMHnk5dNQyiyISsMy0FKpq66ir86i8f1CF34H3zWymmV0XUIbdqmvrmLakjKEaxikicSAzLTRPWFVtdFr9QXX1DHP3NWbWDvjAzBa6+8T6B4S/EK4D6NKlS1TDzF4ZWmZR8/OISDzI2LX8YnUdzdIbf7LIQFr87r4m/HMD8AZw/F6OGe3uBe5ekJ+fH9U8k4pKtcyiiMSN3evu1kZnZE/MC7+Z5ZhZ813PgRHAvFjnqG9KURlHd2ypZRZFJC5k1mvxR0MQLf72wCQzmw1MB95x93cDyAHAjsoaPtcyiyISRzLD3TvRGtkT8z5+d18CDIj1eRsyfdlGaupc8/OISNzY3eKP0k1cST+cc0pRKRmpKRR01TKLIhIfvir8TaerJ65MKirjuK6ttcyiiMSNXcM5m1Iff9wo217JgrVb1c0jInElM11dPVEzpbgM0DQNIhJf1NUTRVOKS2memcYxWmZRROKICn8UTSrSMosiEn++6uNXV0+jWrmxnJUbKzhJ/fsiEmd2tfijNVdP0hZ+TcMsIvFKo3qiZFJRKe2aZ9JLyyyKSJz5alSPCn+jqatzphaXMaxXnpZZFJG4k5Gq4ZyNbuE6LbMoIvErJcXISI3eKlxJWfinFO/q39eFXRGJT5lpKerjb0yTi0rpkZdDh5ZaZlFE4lNo3V119TSKqpo6pi3dqG4eEYlr6uppRLNXbaa8qlbdPCIS1zLTU1X4G8ukxeFlFnuoxS8i8SszLYUqdfU0jinFpRzTqSUts9ODjiIi0qDMNHX1NIodlTXMWrGZoT3V2heR+JaZlqpRPY1h+tLQMosn6cKuiMQ5jeppJJOLSslIS6GgW+ugo4iI7JO6ehrJpKJSCrq2plm6llkUkfiWmaZRPYesdHslC9dt0/h9EUkIoRa/unoOya5lFof21Ph9EYl/GZqy4dBNKSqleTMtsygiiUF9/I1gUlEpJ2iZRRFJEKE7d9XVc9BWlJWzalOFhnGKSMII3blbh7s3+nsnReGfrGmYRSTBZKalUOdQU6fCf1AmFZXSvkUmPfO1zKKIJIbd6+5GoZ+/yRf+3css9tQyiyKSOHavu1vd+P38Tb7wL1i3lY07qhiq/n0RSSCZadFbcL3JF/4pRaHx++rfF5FEoq6eQzCpqJQe+VpmUUQSS8buFn8T6eoxs3PMrNDMiszs9midp6qmjulLN2oYp4gknN1dPVG4ezfmhd/MUoFHgHOBI4HLzOzIaJzri5Wbqaiu1fz7IpJwmlpXz/FAkbsvcfcq4P+AC6JxoklFpaQYnNhD/fsiklh2j+ppIl09nYCV9V6vCm/7GjO7zsxmmNmMkpKSgztRq2ZcdFxnLbMoIgknLzeT8445jNbZGY3+3haN24H3eUKzi4Gz3f2a8OuRwPHuflNDv1NQUOAzZsyIVUQRkSbBzGa6e8Ge24No8a8CDq/3ujOwJoAcIiJJKYjC/xnQ28y6m1kG8H3grQByiIgkpbRYn9Dda8zsRuA9IBV42t3nxzqHiEiyinnhB3D3McCYIM4tIpLsmvyduyIi8nUq/CIiSUaFX0Qkyajwi4gkmZjfwHUwzKwEWB6DU+UBpTE4z6FIhIygnI0pETKCcjamxsrY1d3z99yYEIU/Vsxsxt7ucosniZARlLMxJUJGUM7GFO2M6uoREUkyKvwiIklGhf/rRgcdIAKJkBGUszElQkZQzsYU1Yzq4xcRSTJq8YuIJBkVfhGRJJN0hd/MnjazDWY2r4H9ZmYPhxeCn2Nmg+Iw4xXhbHPMbIqZDYh1xnCOfeasd9xgM6s1s4tilW2P8+83p5kNN7MvzGy+mU2IZb56Gfb3372lmf3bzGaHc14VQMbDzexjM1sQzvCLvRwTD5+hSHIG+jmKJGO9Yxv3M+TuSfUATgEGAfMa2H8eMBYw4ARgWhxmHAq0Dj8/N4iMkeQMH5MKfERoNtaL4jEn0Ar4EugSft0uTnPeAdwTfp4PbAQyYpyxAzAo/Lw5sAg4co9j4uEzFEnOQD9HkWQM72v0z1DStfjdfSKhD0xDLgCe85BPgVZm1iE26UL2l9Hdp7j7pvDLTwmtYhZzEfxdAtwEvAZsiH6ivYsg5+XA6+6+Inx8IFkjyOlAczMzIDd8bE0ssu0O4L7W3T8PP98GLOCba2bHw2dovzmD/hxF+HcJUfgMJV3hj0BEi8HHkasJta7ijpl1Ar4LPBZ0lv3oA7Q2s/FmNtPMfhh0oAb8FTiC0FKlc4FfuHtdUGHMrBswEJi2x664+gztI2d9gX6OGsoYrc9QIAuxxDnby7a4HPNqZqcR+gd7UtBZGvAQcJu714YaqXErDTgOOAPIAqaa2afuvijYWN9wNvAFcDrQE/jAzD5x962xDmJmuYRaoTfv5fxx8xnaT85dxwT6OdpPxqh8hlT4vykhFoM3s/7Ak8C57l4WdJ4GFAD/F/4HmwecZ2Y17v6vYGN9wyqg1N13ADvMbCIwgFCfazy5CrjbQx2/RWa2FOgHTI9lCDNLJ1SoXnT31/dySFx8hiLIGfjnKIKMUfkMqavnm94CfhgemXACsMXd1wYdqj4z6wK8DoyMw1bpbu7e3d27uXs34FXgp3FY9AHeBE42szQzywaGEOpvjTcrCP1fCWbWHugLLIllgPD1haeABe7+QAOHBf4ZiiRn0J+jSDJG6zOUdC1+M3sZGA7kmdkq4L+BdAB3f4zQlfPzgCKgnFArK94y/hZoC/wt3BKo8QBmG4wgZ1zYX053X2Bm7wJzgDrgSXff5xDVIHICfwCeNbO5hLpTbnP3WE8vPAwYCcw1sy/C2+4AutTLGfhnKMKcQX+OIskYFZqyQUQkyairR0Qkyajwi4gkGRV+EZEko8IvIpJkVPhFRJKMCr8kPTP7yMzGhG+maeiYLDObYGapjXC+/PDwUZFAqPBL0nP304FK4Fv7OOzHhCZyq22E85UAa81s2KG+l8jBUOEXCRkLXLGP/VcQusN319z9b+/aYWZ/NbMfhZ8vM7P/MbOpZjbDzAaZ2XtmVmxmN9R7v3/t53wiUaPCLxLyfWCEmbXYc4eZZQA93H1ZhO+10t1PBD4BngUuIjQv/e/rHTMDOPlQAoscLBV+SXpmdgzQEngJuHAvh+QBmw/gLd8K/5xLaHGPbeHunZ1m1iq8bwPQ8SAjixwSFX4RuBl4EHiBvXe/VADN9thWf47cPS8KV4Z/1tV7vuv1rvmxmoXfVyTmVPglqZlZPnAq8A93nwx0NbOvtcTDqzSlmln94n90eKRPBnAioeXxDkQfIOYTwYmACr/I9YRm46wOv36ZUH//nt7n6wt1bALeBaYC44DfmFnOAZz3NOCdA48rcug0O6dIBMxsIHCLu480s+HAre5+/iG830TggnprvorEjFr8IhFw91nAx411AxfwgIq+BEUtfhGRJKMWv4hIklHhFxFJMir8IiJJRoVfRCTJqPCLiCSZ/w/Fq79XR5zzzgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = [0.95,1.05,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2,2.1,2.2,2.3,2.4,2.41]\n",
    "y = [0,9.82,14.35,19.34,22.08,23.62,24.36,24.52,24.2,22.47,22.36,20.85,18.92,16.45,13.21,8.39,0]\n",
    "plt.plot(x,y)\n",
    "plt.xlabel('λ (μm)')\n",
    "plt.ylabel('θ (deg)')\n",
    "plt.title('xz-plane Φ = 0°')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}