Dispersion_relation_chi_2_and_interpolations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving phase and frequency matching conditions for $\\chi^{(2)}$  computationally\n",
    "\n",
    "# 1D interpolation for solution space\n",
    "\n",
    "Refractive index from:\n",
    "\n",
    "http://refractiveindex.info/?shelf=main&book=LiNbO3&page=Zelmon-o"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notes\n",
    "\n",
    "There are three frequencies, $\\omega_1,\\omega_2,\\omega_3$. By the frequency matching condition, we impose $\\omega_3 = -\\omega_1 - \\omega_2$. \n",
    "\n",
    "The user must be watchful of the range for which the index of refraction equation is valid. In the case blow this corresponds to $\\omega \\in \\pm [6,60]$ in units of $10^{13} Hz$. \n",
    "\n",
    "In the code below, `fixed_omega1 is the value` dedicated to a fixed value of $\\omega_1$. There is a single degree of freedom remaining since $\\omega_3$ depends on $\\omega_1,\\omega_2$. We initialize a list of tuples `oms` with the values that will correspond to the initial values of $\\omega_1$ and $\\omega_2$. We must make sure that all three values are inside the admissible range. When the values are in the 'corners' of the ranges, Newton's method convergence slows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with 3-wave mixing interactions. \n",
    "\n",
    "The conditions we need are:\n",
    "\\begin{align}\n",
    "&T \\left|\\omega_1+\\omega_2+\\omega_3 \\right| \\le \\epsilon_1\\\\\n",
    "&\\left| L\\left(k_1+k_2+k_3\\right) \\right|= \\le \\epsilon_2 \\implies \\left| \\frac{L}{c} \\left( n(\\omega_1)\\omega_1+n(\\omega_2)\\omega_2+n(\\omega_3)\\omega_3 \\right)) \\right| \\le \\epsilon_2.\n",
    "\\end{align}\n",
    "In practice we use small cutoffs $\\epsilon_1,\\epsilon_2$ instead of zeros. These are dimensionless numbers. The $T$ is the duration of the simulation or decay time of the system modes (whichever is smaller), and $\\epsilon_1$ maps to the fractional error in the result. $L$ is the legnth of the cavity, and $\\epsilon_2$ is chosen so that $L \\Delta k / 2$  is sufficiently small so that that $\\text{sinc}^2(L \\Delta_k /2)$ is small ($L \\Delta k / 2$ should be on the order of magnitude of 10).\n",
    "\n",
    "One way to solve a set of equations is to use Newton's method a few times. We do this here to simulateneously solve the phase and frequency matching conditions for $\\chi^{(2)}$ nonlinearity.\n",
    "The main idea is to trace along the solutions, incrementing $\\omega_1$ at each step and finding the new $\\omega_2$ and $\\omega_3$ at each step.\n",
    "\n",
    "Another point of consideration when considering the cutoffs $\\epsilon_1$ and $\\epsilon_2$ is the number of modes being dropped. If this number if very large the approximation is not guaranteed to work well.\n",
    "\n",
    "Later I also interpolate the dispersion relations in the appropriate region. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "import numpy as np\n",
    "import scipy.constants\n",
    "from sympy.utilities.autowrap import ufuncify\n",
    "import time\n",
    "from scipy import interpolate\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from sympy import init_printing\n",
    "init_printing() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_arr(arr):\n",
    "    fig = plt.figure(figsize=(15,15))\n",
    "    ax = fig.add_subplot(111)\n",
    "    cax = ax.matshow(np.asmatrix(arr), interpolation='nearest')\n",
    "    fig.colorbar(cax)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## from https://www.andreas-jung.com/contents/a-python-decorator-for-measuring-the-execution-time-of-methods\n",
    "\n",
    "def timeit(method):\n",
    "    def timed(*args, **kw):\n",
    "        ts = time.time()\n",
    "        result = method(*args, **kw)\n",
    "        te = time.time()\n",
    "        print '%r %2.2f sec' % \\\n",
    "              (method.__name__, te-ts)\n",
    "        return result\n",
    "    return timed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lambd,omega,omega1,omega2,omega3,omega4 = sp.symbols('lambda omega omega_1 omega_2 omega_3 omega_4')\n",
    "l2 = lambd **2\n",
    "\n",
    "def n_symb(pol='o'):\n",
    "    s = 1.\n",
    "    if pol == 'o':\n",
    "        s += 2.6734 * l2 / (l2 - 0.01764)\n",
    "        s += 1.2290 * l2 / (l2 - 0.05914)\n",
    "        s += 12.614 * l2 / (l2 - 474.6)\n",
    "    else:\n",
    "        s += 2.9804 * l2 / (l2 - 0.02047)\n",
    "        s += 0.5981 * l2 / (l2 - 0.0666)\n",
    "        s += 8.9543 * l2 / (l2 - 416.08)\n",
    "    return sp.sqrt(s)\n",
    "\n",
    "def k_symb(symbol=omega,pol='o'):\n",
    "    '''k is accurate for omega inputs between 6-60.'''\n",
    "    return ((n_symb(pol=pol) * symbol )\n",
    "                .subs(lambd,scipy.constants.c / (symbol*1e7))) ## / scipy.constants.c\n",
    "\n",
    "def make_newton_update_expression(dispersion_difference_function,symbols):\n",
    "    expression_diff = [sp.diff(dispersion_difference_function,sym) for sym in symbols]\n",
    "    return [- dispersion_difference_function / der for der in expression_diff]\n",
    "\n",
    "dispersion_difference_function = k_symb()\n",
    "update_expression = make_newton_update_expression(dispersion_difference_function.expand(),symbols=[omega])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/UAAAAyBAMAAADvkOZbAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAuyLvRGYQdpmJVM0y\nq91NqyaOAAAACXBIWXMAAA7EAAAOxAGVKw4bAAARN0lEQVR4Ae1bfYxcVRU/O7szuzOzXwXTEMPH\nQGMCkthCocZWZRJJqGLsxIqVEO0IiYgEdgMJn5oOHyGttOnWpmFXwI5p8INIdlGByFK6f5YY7URE\nili7oIlgA7Sl3UpLWX/nd+5983Y+dme6s5sN3Zv2vnffPfec3/2deW/fe/f9REpKfHy+nB4MHCzJ\nvEhn2ZH5A6cLA+efLhOdn2cZA78sOzJ/4HRhYOPpMtH5eZYxMFp2pJ4DkfULRmDfnpPYhmtHxDXl\nz3sHBvoiA1vyIhfed7bI97c9JrB4POVs7+xPi1w08EmR57K789aNrpktREisGkcDt7+zYWCgPDpm\nIzJw9wgrs0ne0f5jYYWumSgkyxiLPwKyUG68+0EwN/AXMQgxJa3ItzFHG9ERnvo6wMWG6jAuN71F\n4kdwtLNbkvnYm+KashZ3zoULJXFUYgXZnpbn5daU3JCJXGq2bSJ/l0hOBjNyxfgHYt2SWFbuvwFH\n6FYrIiRW9aqBo4B50kXfkcVBg4DZyM3ptiFWZtM8Pt4nrGDVwBKAI1nG2O/lYUbISUc+drksShkE\nklbk25hTG4npCE99HeASfXUYl5v+TWQpjt7eLatFHhXXlPtFOmRdXj6Sjqw0L4kukfiQvCyy12yf\nELlYWo9L16j84dWUWPctm8bK/U//CN2yIkJiVa8auC0lUmD02N79WREHAbORs6SpmxVtJPn4bfjx\naNXQUgRHslg15WSBBklkJDLUsUQ6cgaBpDmCFaFOgDbCEa6nHnjxdD3WZbZfFNmPS/113XKOyGDW\nmiL4QT0mX0rHTkjzqHR82FSQxJicENkutN0icqu0XCY9Q1KAS+uW1rEy9404QLdaESHjq1sNHMe/\nERd9XRYHaKuzaTmMFiva4JeMwkp3GlcCcCSL1aqUuY93S2y0a0jiBw0CSSvybRNQG+EIT30d2JLZ\nOozLTXdlNffxZLecFFmVt6aatfSiwjU/OYbcdxWk5XjskMg6s121ULapEa75mgJ2O+L1cGNLQC8R\nEqsG0MAoV/noxdzrbJqWoI8VbWY+9ySL1ZkKC6XlUDae6emW1g/RuArsKmlFvjkB2ghHFKnn6Fqq\nZC1Gk9q8l5WvJ7tj/wM6zTeaWpq0urVP686xnpy0HJMVOO/Taisd49dm0BF7W+SVNRvFumf6vHcI\nNb4WDYzS7aMXc68WXS888qBVtJGOa+4asUoHNa4EP0xHFhg78dJdaQbYf3SbNOO81xsqYHakKcGc\nAyegNsEIR33N6FbWbFnFsOUDifUmuyO4ZduDVKPJouf1RW9ldX9V3+I+af9AL/i7MmqL38AJ9ER+\n9pDIQlmctu6Zzr0htPhApYFxyc+Iix7knhY9z0hznhVtpCkbOWIVhjawFHNPslDFTmRkOSPEx5/F\nBV+SRw2CkaYE2xw4AbXxIzz1NcO7p2bLKoadBfzZRO5x3u9B7tHUEjuodfwyrd+RxTnNfbI3sT2v\ntpL463vs2ZpGd9Oodc947omQ8RWVBhbpEQdOgtzTomdM2paxoo2anxtUutOgUsy9kRW/LDaelQey\n6v6pr72flufk18cNgiPN8a0GmIDa+BGOeu2prSyszay61XqR3yCf/pqPppbWIW6WZkXacv6yuuaH\n20fUVr4tre+jR5qXoWo7PJvXfMZHVAYWuRe88g9SMfe06BqV6BFWtFHzzSlfaatBJZR7UbJQfSSy\nNo29aEHWLZPonS99aBAcaY5vGIA5s3EjHPXaU1sZrc2sqlVrTmJ9mk/cSQ3mBU0WfUS+QOS8vMhd\nvJlL4NeLZ4Lv0RaP+XvykZQ0He/Qv2e419Nu8sDRDa0CehXhLxgf/hlY5I0AnD/vbTadQxI9yoo2\n8rreq7BqKDabsyIkWawOIPcZRGlOSbuS1oFewBSSpgQbQk7AbGyEp75mgO2Fmk0rG/5KIp/et2//\nF3oXAHJW0EypYQ/84uKF3CdycgEelaOYgcjOhNrecRAzynSNIfe4l2kbc90znXtF+ANiBRIGlhh+\nDy66zz0R9uK5OnqEFW30PnVzlpVOo3El+GGSLFZ4DcLzHn+O5BX8TxYIIUbSlGDjmxMwGxvhqa8Z\nXGumZtOKhu05aU2BSb7beVjf7WoTp3VO5BLc9Kdkt8hteH2Ct2SrU+361Axb/ISbUth2Hm7DHcKQ\ndc/4eb9a+MIMcVEYWNpxD+Ci+9wbwnb8vV/CijaSEzlglQ5uXAlyT7JY7cHfew2Ac1rOiF4qe3D+\nA6aR5gjGHDgB2ghHuB4dWWPpSNVoWMXsMwP9/0RXF/KYj31e2OwZwsU1J3KNJI5J5O2BTaPyX7lh\nRO7N3qy/NNhenZU1Es3IYG8kh5+Edc947omQ8fHOyQInlFSCK/69p4W8JjfkraLNhYIssMKABpYg\n9ySLVUdfbDnu3aXlWYn2tj4bewZP/QqTpBX55gRoIxzheurA1lmHbSXTXePjuNOIn3cyE3vkRyPC\npt5vLsrjcW/9fWlJ4pX5qNx05z9gtf5BuFDb9r39aZGV6z+BpZ7HHxLrbnv5hF7iGl3olhURMj4e\nPSxwZDnCEdymK3ZmxCAoQoAlYlS0iWzYmxVWDcVXBEeyWEn/phF9OJI/bQBf/evTDgJJK/JtE6AN\nR1hPPejODxlzsS3Unt/9eDOwRqe3I6u1W8/i7nx1GjCwEa9huITFe544/lJrmfvf8iR6Dencrec+\nwlElj7e4fjVND8znXlmYXpnzuY/xRLfcd+lim813PvfTy7uOnpO57wawRWmFB4B9Wlvue+zdph6Y\nz72yML0yF3MfTWNOPb0SGcE2rg2X+8VcbNP2fO7JwvSquZh7vA9CzvNYVcM2mdWWnfeLudim7fnc\nk4XpVXMx91eKYJ0iJYMFzC3J+U285j86PPyv4eGd2oP3MnOtKKz48PCLK4aHdQJdcw3f+PjBOYoQ\nuHr1pXtcYt9ZAowr8d+f9122mqYH5v/eKwvTK3PxvMdqS14ukkSbLqz9lPOz8z5YbJvP/fSybqPn\nYu7P1jW23+JVjj7PLQzl3q1n6ZFTP+/7z8rS5UxXp85s/D4sIMxCqRfhbDC3tX/5ya0ZrA/oh3Wj\nJAHnPdaN3HqWHjnl3MdT+KRgNkq9zBYx/US/w5iFUifCWWHuuydGbv+szh3Zti83uISFdSOuZ5GV\nU859py3Rzzy1dTIbAvSmLojPQqkT4awxx6nvwldSOP8rlbLc+yU+vxVT3FEIhp9QzjvpTONz1gmF\nKjIq8biHPmr1TNpGMZ6ahxRn7DatHjVxpjhTI+vRvQpvzZzYzmSB9O7kf2a/2zao/83vc4Mm1jNU\n51aO0PRuAUKLTttShH5iUyAkT24OjjGvvAOa+pgrKgNLmTNJ3yTMcebbR6TalxvXh6jhLvV02LPt\njiw+K8LDVYFCMODuphEr6G8mlDaV3lGJxz30caRJ2yjG02XEkOKM3fgwAboNauJy0pFHrpe5gc55\nLOV2/MaJ7Ux0R+8m/6PgzT5486aDWb/HbRWEpncLEBps2pYitIlNiZA8caxnzNHn0NTDnFMG7ihn\njpK+MHOwMewuDDcgt9YvN3AbyCU+bm3d7358cmdCsEDhRq+dvRbDX/KeUOkdPijbK9xDJ0eatE0l\nZXTXllKhDC3YbVo91cQlVJVmQjn2mPfyGh9kPapHuaV3BnUKu1Z8DRWU/7g9d3WrgpB6N1ahCVdC\naLAZaFKE5MlZqJ4O36055Z3hqYc5TrAic+tUBEkhIZkLYXez5qanIOEvN8JdpftOsuaX+vSZsE8V\nd6tSaqn6taD0uz2f+y0qvcOtxXbhHno5soNmBdZwF1KcsZtaPZPDUXHG77nY49yXbc7Bd2JZHOWW\n3hnUfQj2qaeLA+I5t+9yXwUhfosoAUKDTdsShG5i+vXVpAjJk1k4xow+h0bqYc7o0xeypcxR0hdm\nLkiWj6Pb5lF8NFdb6XJLfG6r7vBJVq8JwQKFmx6M5rRG8blXFZkp8bwIjyNLco8BVznFGbvRHsxQ\nDmeKM5dChKxWTmJ4Hp3cqncn/+MHcdIXyv23vAuX+yoIS3NPXG4OhsMhNKHc1AgdTzpWFX8oZ3ok\nuq2LuWLuMbKEOYggw8wFyQoHw3ffG8PtSfZ7cqqnmyhjUMWdSceoDtPJLErLl+Wr5sfnniqyFfp1\nu9eT4VqHvxeqb3OaOEMHB84C3Ujd204OR8WZY5Y95r+k9tIQ29I7g9rAaAq5dwhb+lpGbLDLfRWE\nBs5Ue4ZQYZuQcCJCO2g/suoIxVR0Onm3Z/RpOwKVUl3MGX1uIaaEOUj6qCYUYy7AbpO2uvWIjIbb\nk+z7JT63NXfbnBCM6rBoGsN7etuff3Kh+fG5p4qMSjynJ0P3NqRf9W14tbQY4+gunkEHZXqU6KpW\nz+RwqjjDEwmuqByomwrFywFtS+8WlANvlKfFIZQ/PvnzrDlwua+C0MCFEBJXJYQGeyqETkWnXpwm\n0HR03wCY2DGpjzmjryJzqusLM+eTNYG09mP25caEg5Ubi3Oqp5soXYNWwIRg1K/hXheXsjzuPw+a\nC597qsioxPMiPKfVM31b06jLvaoMzEJFCChb0yaHoyrNcu962F1SeTmg38I7g9rAPuTeIcTt+7gb\n63JfHaGCU9GbcazRKyK0g8z9JAjxl1mv9GrhNIFGXw4x8G69bubOxTDmtZQ5iCDDzNGmDNch+3JD\nY1cpsRffRTkwUuma3zokeke5Nk39ml8dNEehhTZTkakSz4vwvFZvcwrGbVBrEN29aJiFOkZpXtY1\nCm1MlKo0S6HrYbevblaE774w8ZqPTnjXoByYGEHuJyKU4mpldYQKrohQo1dC6A/qlWkShFQtmoXX\nBCp9V+IzdvwR8NOpnTmlrzJzS7Nh5mhThut3bWkfcoptl1vic1u604f6AwD/TSrc3Opg0Y8/7y9V\n6R0O78/6Pb4OeJ13ACrGczN4AzZmoY4h7Go63qmauOYUVWk8q8LvEYqRbA/3eIN57HJL72jszzL3\nNwlyX4bQv7n2uEoQ4uV0/Agrh1CjV0LoDk6FMGZKQHhxe0Yfbrd38wHBpoG6FubcBJmIicxdoCLI\nMHNBsoIAurP54gnNSRpNTk/ntnTXUxB9bl/7FSrccrY6WPThZqAXm44MDu8M9qjVW6H6tmYV4xmz\nKoxzFuqYWj3K4ZwqjcxqT7WyAFCy6OSW3tHYifMMAZ7at+/9V8sQutwHuEoQElwIocKuiNAdnAqh\n6fzUi9sz+laK3CORvtC0Er1sVMPFubkJaiJKmBtXEWSYuSBZoRDg6roJzUkaUerp8BxiW7rbk8MZ\nbdKx5m7xq4OBEzcDPSmaUlTicU/7dST+HdA/wfiF8spVVJyxGx47D1MOh/NezrAUsgfjK5bVJrbT\ndzsPm3eT/zElqg0sQxg67yshJLgQQka3OegEightivZXSXuqF4xxc9A9o69pBFNPpkKDamEux4tG\nJeYuURFkmLkgWaEQuHN7a0JzsoZK1opLfXQ3iPgUglG/5lcHAyd+Bler9I5KPO55rR6lbRGK8TiD\nouKMUr6oavUohzPFmTGrIasViu3gnVt6N/mfy/0hKUPoc18FIcGFEBIXbUsQmlBuaoSm87M5QJMI\nrEpf/HMPDL52eXhWtTDnlIGaiBLmKOkLMxckKxxDunQFv7bCJb53wtI1Ku5MOqb6tWB10PvzM6CK\njEo87nUWTKtn0jaK8biMSNGbyfQW5eGDWj1q4qg4M6Ece3yAki3FdvDOLb0zqA2ULeOvlCH0ua+C\n0PRuRYSccEWEPDg1QlMCcg7KmDIB5V3svBdal/aGJ1MLc0ZfJeZM0ldkjjaWrHCQJj4zh480cN/P\noIEuG+3KP9832u80/c0Kc/HRaaKcbHgsNVnvnOi7fk6gKAMxK8y19JXFnT/wcWbg/yMSvcAgjX+/\nAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\omega \\sqrt{1.0 + \\frac{2402.73209483501}{- 0.01764 \\omega^{2} + 898.755178736817} + \\frac{1104.57011466755}{- 0.05914 \\omega^{2} + 898.755178736817} + \\frac{11336.8978245862}{- 474.6 \\omega^{2} + 898.755178736817}}$$"
      ],
      "text/plain": [
       "       _______________________________________________________________________\n",
       "      ╱                2402.73209483501                     1104.57011466755  \n",
       "ω⋅   ╱  1.0 + ────────────────────────────────── + ───────────────────────────\n",
       "    ╱                       2                                    2            \n",
       "  ╲╱          - - -0.01764⋅ω  + 898.755178736817   - - -0.05914⋅ω  + 898.75517\n",
       "\n",
       "______________________________________________\n",
       "                    11336.8978245862          \n",
       "─────── + ─────────────────────────────────── \n",
       "                         2                    \n",
       "8736817   - -474.6⋅(-1)⋅ω  + 898.755178736817 "
      ]
     },
     "execution_count": 160,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dispersion_difference_function.expand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "k_o_vs_omega = ufuncify([omega],k_symb())\n",
    "k_e_vs_omega = ufuncify([omega],k_symb(pol='e'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#####  Different granularities for the range of omegas\n",
    "\n",
    "omegas1 = np.asarray([6.+ i*5 for i in range(int(1e1))])\n",
    "omegas2 = np.asarray([6.+ i*5e-1 for i in range(int(1e2))])\n",
    "omegas3 = np.asarray([6.+ i*5e-2 for i in range(int(1e3))])\n",
    "omegas4 = np.asarray([6.+ i*5e-3 for i in range(int(1e4))])\n",
    "omegas5 = np.asarray([6.+ i*5e-4 for i in range(int(1e5))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x16159fdd0>"
      ]
     },
     "execution_count": 173,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEKCAYAAAAxXHOuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FOX2wPHvoUpvCl5AqoqUJHSkSVAuoDQvSFMxECA0\nAeWnCFiAKypgQ0F6UVAE5CI9gEhCFUJLpVkQUYrUUIQQkvf3x0xwCRuIYZPdbM7nefJkd3Z257wz\ns3P2vNPEGINSSinlKtncHYBSSinvoolFKaWUS2liUUop5VKaWJRSSrmUJhallFIupYlFKaWUS6U5\nsYhIIxHZ78pgMoqIlBWRRBFJU/tFZLiITHd1XMmmcVhEHk/hNY+Z98nnpYisFpFu7o7LGTvOCvbj\nKSLyurtjSpIR65Q3EJEAEdns8PyiiJRzX0TOiUgTETnq8DxaRB5zZ0yO0vt7miMVARwGehpjNjgO\nN8ZsASqnV2AZIFUn8IhIE+BLY8wDN95ozHvpFlUqeOC8vzEvjTFPuTOQO3CMs587A0nO3evUnYhI\nWeAwkMMYk+jmcByXYwF3BnIHjnFWc2cgyaX39zTTdYWJSHZXjPNPJkkqk5C3c/F8dce0xQWfkbYJ\nu3HeuUjS9yDFeejqNqa1R8FF087Uy8vd8d9NV1jyUu+wiPyfiESIyDkR+VpEcjm83lpE9tqvbRER\nH4fXXhORn0Tkgl0yPu3wWoA9/kcichoY6SSWkSLyjYjME5HzQIBYhtmfe0pEFohI4RTa0l1E9tnT\n/0lEguzheYHVQEm75L4gIvfb05vn8P62dtxnRWSDiDySmvkiIsVEZIU9/IyIbEwWWo0U3uds3g8T\nkRj7c2Y5zvvbEZGCIjJXRP60P+d1h9dumfcikk1EPrDn6U9Aq2SfFyIigQ7v3ywi79vz5mcRaXmn\n+e7YRhEZKiLHgdkiEiUirRzGyWHH4ZdC214VkWMi8ruI9MDhB4KIzBGR/95pOdxp3t5hvT5sxx8B\nXLLn3Wt2PBdEZL+INLXHdck6dScikstefkdE5LhYXYK57deGish2+btbs589z3MBSfPkvB17vRTW\njwoi8r2InLbXqS9FpKDD9B+x15Fz9me3SbZMJovIKhG5CPiLSFERWS4isSKyHaiYrD2O3ZtzRGSS\niKy0Y/xBRMo7jDtBRH6zP2uniDRyeC35NmSYiFwWkSIO49S023TLRltE7hGRz+3lFQ3USfb6ja5t\nEaljTz/WXgYf2MOTupV7i8gf9t//OXyGSArbNIf3BorIEeB7Ecltz//T9vzeISL32eM7fk9FRN4Q\nkV9F5ITdjoLJPvcFe535U0RG3H4tA4wxt/3DKn8fdzK8CfBbsvG2AyWAwsA+IMh+rQZwEqiN9Yun\nmz1+Tvv1DkAJ+3FH4JLD8wAgHuiPlQhzO4llJBAHtLGf5wYGA9uAfwE5gSnAfPv1skACkM1+/iRQ\nzn7cGLgMVHfWTofpzbUfP2zH+ziQHXgV+BGry+BO8+VdYLLdruxAw1TOT2fzPhIoaY+7BfjvnZat\n/d65wLdAXnu+HAR6pDTvgb52LEnT2pBsXoYAgQ7vjwMC7eXeF/jDYdp3mu/x9jzKaU/7FWCBw/vb\nAREptKslcByryzAP8JUdZwX79TlJ8ygVy8HpvOXO6/VhYI/93tz2uvIbf6/bZYDyrlynUrG8PwaW\nAoWAfMAy4B37NQFCgbeAB4GzgG+y74w4fJaz9aMi8ARWN3sx+/M+ssfPYbfjNftxU+AC8JDDMjkH\nPOrwPV5g/90DVAV+BzY5xJB8mZ4CatnxfIn9nbdff9aeX9mAl+31I1cK25B7gJVAH4f3fwR8ksJ8\nHYuVfAsBpYAobv2OPm4/3gY8Zz/OC9R1mMeJWOvqPUA14E+H991pm5YIfG6/NzcQZC/f3PayrQHk\nd/I9DQQO2Z+RF/gff6+LSZ87DcgF+AJXgUq3Xc9cnFi6OjwfB0y2H08GRid7/wGgcQrT3OuwgAOA\nX+8Q40ggNNmwfUBTh+f/Aq7ZK9VNicXJ530LDHTWTicbgTe4eWMnWCv/Y6mYL6PtaVVMYb6n9D5n\n8763w/MngR9TsWyzYX2ZKjkMCwI2pDTvge9x2IgB/+b2ieWQw7h57HGLp3K+X8XeSDssw1j+/nJ8\nA7ySwmfNAt51eP4QKSeWOy0Hp/P2Tuu1/d4Ah9cqAiewN7zpsU6lYplfwk5m9vP6wC8Oz8sCZ7C+\nP0OTDb/pO+Ns/XAyvXbAbvtxY+BYstfnA285LJPPk62f17ATjz3sHW5OLInJlun0ZMtq321iOwv4\nOMz/5NuQTsAWh1iOA7VT+KyfgX87PO9Nyokl1J5esWSfkbQRd2zvOGCG/Tg127SyDq/3wPoh5OMk\nXsfv6Xqgr8NrDzv53H85vL4D6HS75e7qPsyTDo//AvLbj8sC/2eXiWdF5BxQGuuXHHaZldSdcA7r\nl8m9Dp91lDtLPk5Z4NukaWItlHisX3k3EZEn7bL5jD39J5NN/3ZKAkeSnhhrzh/F+tWSJKX58j7W\nCrnOLm9fS/bZKb3Pmd8dHh+x47qTe7F+Of6W7L2OsSefryWTDTvC7Z1IemCMuYK1kcwPqZrvp4wx\n8Q7vPw5sBTqISCF7/K9SmK6zOFPaP3Cn5ZDSvL3tep38vcaYn4GXgFHASRGZLyL3pxB7WtepFNnd\nIHmB3Q7fi2CsyiJpWkewNjplsRLnndy0fohIcbG65n63u5S+5O9l+q/k43P79e0+rIot+fy/nRMO\nj2+aLyLyilhdr0nbmYLcfjuzDKgs1oELzYHzxphdKUy35D+IsydQCThgd085dicbJ5/juL7daZvm\n+N65wFpggb08xjnrxiPZ+mY/zpHsc//R+pZRO8eOYpXbRe2/IsaY/MaYhSJSBpgO9LeHFwFiuHkj\nYFIxjeTj/AY8mWya+eyN0w1i9R8vBsYD99nTD3aY/p2mfQxrgTt6gJsXsPOAjblkjHnFGFMRaAsM\nEbvPPQ0ecHhc1o7rTk5jrZiO8ZcF/nAMM9l7jjuZ1j+WivnubNpgfVm6YXWZbku+PO8Qp9NlmYrl\nkNK8TXG9TqkNxpgFxpjG/D3fxjkJKc3r1B2cxtooVHWIubAxplDSCPZGrj5WZfqBY+gpfGby4e9i\n/equaowpDDzP38v0GDfPS7C6A1Na304B15O9p0wKcdyWiDTG6lJ8xmE7c4HbrG/GmDhgEdb69jww\nj5Qlb1uK3wtjzM/GmGeNMfdhrf+LRSRPUqjc2t6k9S0127QbbTDGJBhj3jbGVAUaAK2BF1KIPfk2\nIJ6bk8k/ktrEksveEZT090+POJgB9BWRugAikk9EnhKRfFj9vInAabF2bvbA6lu8W9OAd+3EhYjc\nJyJtHV5PWqFy2X+njTGJIvIk1q+TJCeBYuKwAzKZRUArEWkq1s7kV7C6cH64U4Ai0kpEknZGXsT6\nEiWksn3JDRCRUiJSFBiB1S+dNJ1EcXIMvbEOG10EvCMi+e1fZi9z+y/QImCQPa0iWP3laXGn+Z6S\npUBNYBBWkrldnN1FpLJYB2G8ldKIqVgOKc3b263XzqbzsL2e5MLqariCte47iz1N65Q9nZSWt7Fj\nnuCwE7eUiDS3H99rvx4IdAda28sFrI18Isl2njtRAKu77aKIlMLamCfZAfwl1kECOUTEH2tj97Wz\nD7LXzyXAKBHJIyJVsLrf0iI/1sbyjFgHMLxlx3on87DmRRtu/734BhguIoVFpDTwYkojishz9rwG\nq2vXcPN68Kbd3qpY3VlJ61tqt2lJ0/EXkWpiHYxxCav9zrYvXwMvi0g5EcmP1d24wPx9WPk/Ppoy\ntYllFdYvnSv2/5FOxknxl70xZjdWn+Mku4Q7hL2CGGP2Ax9i7Yw8gdUNtiWVcd3OJ1il7DoRicXa\n6VU3ebzGmEtYG6lv7Ni62O9Liv0g1oz/xS5Bb+q6MMYcwvo1Mwnry9cKa//QdcfppOAhYL1YR8Bs\nBT4zxmxKxfucmQ+sA37C2kH6DoCIPID1yywqhfcNwlqmvwCbsM7ZmXOb6czAKq8jgF1YO/oc3Snu\nVM33FN9szFV7muWxNjopjbcGmIB1cMEhrF/gKbndcoAU5u3t1mvHtjrIjbWT9xTWr8T7gOFOYk/z\nOpWK5f2a3Y7tdlfVOqw+dbA2XN8aY9YaY84CvYAZIlLE7sZ8B9hqfw/qOvtwrP1VtYDzwAoc1g+7\nW7MN8BRW9TQJ6GaM+fE27RqIlQCOA7PtP0ep/Z6stf8OYe3v+ItUdLEbY7ZhbfT3GGNuN/5orIri\nMLCGW3/0OMbZEogRkQtYB1N0tqujJBuxltF3wHhjTNK6m6ptmoP7sXoFYrF6gUKwuiaTjzsbK2lu\nwuoS/gvru5nS595xnou9M+YfE5FZWL82ThpjfJO99n9Y/db32isoIjIc65fQdWCwMWZdmiasbiEp\nnMRqv/YcUMUY4zFnmd8tEXkTawens7Le1dNKcd56Im9c3u4mIt8DXxljkic1V0+nLNYPvJzG/Seh\n3pU7nnl/G3OAiSTLzHYZ+G8cdgaJSGWsIywqY+3cXC8iD5m0ZjWVasaYlHZuZ0p2d1RP4Dl3x+KJ\nvG15u5uI1ME6TLftncZ11SQzaDrpKs077411WZFzTl76mJv7VcE65HCBMea6MeZXrO6ElEpp9c9l\niQQtIr2wuhtWGWO2ZtBks8S8VbcSkc+xugoHG2MuZ9BkvWJ9u5uK5Rb2jqSjxpgokZsSbylu3vH4\nBzcfYqjugjGmgrtjyAjGmJnAzAyeZpaYt+pWxpjuGTy9I1iHV2d6Lkss9uFyI7C6wZRSSmVRrqxY\nKgLlgAixypXSwB776JE/uPn489LcfOz6DSLiFaWgUkplNGOMR+yjudsTJMX+wxgTbYy53xhTwRhT\nHutkrhrGmD+B5UBn+/jx8ljXIQpL6UNvd6mAzP43cuRIt8eg7dO2afu878+T3M3VjedjHUf9sFhX\nDO2RbBTD30lnH9ZJX/uwrhbc33janFBKKeUSae4KM8Y8e4fXKyR7/h7g0TczUkopdfcy3Y2+Mjt/\nf393h5CuvLl93tw20PYp10nzmffpRUS0l0wppf4hEcF4yc57pZRS6iaaWJRSSrmUJhallFIupYlF\nKaWUS2liUUop5VKaWJRSSrmUJhallFIupYlFKaWUS2liUUop5VKaWJRSKp0ZA7Nnw6BB7o4kY7j0\nDpJKKaVuFhsL/fpBZCQsXHDjou9eTSsWpZRKJzt2QI0aUCx/HHtbv0nVD5LfXcQ7aWJRSikXS0yE\nceOgTRuYGRTGxG21yLk/Et59192hZQi9urFSSrnQ8ePwwguQePkK3/qOpODSuTBhAnTuDJJ+3WB6\ndWOllPJCwcFQsyZ0fWAL60/7UfDcEWvnSpcu6ZpUPI3uvFdKqbsUFwfDh8PqRZfY2WAEpdcshkmT\noH17d4fmFlqxKKXUXTh0CBo0gHw7NhCTw5fS+WMhOjrLJhXQxKKUUmliDMydCy3rxzI3Tx/++1sA\n2SdPgi++gKJFbxo30SQyccdE+q/q76ZoM5Z2hSml1D904QL07w/5NgVzIFcfclVpCauioVChW8b9\n48If9FjWgwtxF5j3n3luiDbjacWilFL/wM6d0NTvLP3DujMlW39yzZsN06c7TSoLoxdSc3pNGpdp\nzJbALTxU7CE3RJzxtGJRSqlUSEyEDz+E6DFL2ZJ9AHmeaw/vRUH+/LeMe+7KOV4MfpHdx3azsutK\n6pSq44aI3SfNFYuIzBKRkyIS6TBsvIjsF5FwEfmfiBR0eG24iPxov978bgNXSqmMcuIEdHniFNXH\nd2Vm0VfJs2wBTJzoNKl8/8v3+E31o+g9RdnTZ0+WSypwd11hc4AWyYatA6oaY6oDPwLDAUSkCtAJ\nqAw8CUwWyUIHdSulMq21awyjKi9k1k4fnnihNDljIqBx41vGu3r9KkPWDiFgaQAz2sxg4lMTyZsz\nrxsidr80d4UZY7aISNlkw9Y7PN0OdLAftwUWGGOuA7+KyI9AXWBHWqevlFLp6do1GDv4OHU/78+H\nJQ6Sb+EyqFfP6bjhJ8J5fsnzVL6vMhF9IyiWt1gGR+tZ0nPnfSCw2n5cCjjq8Nof9jCllPI4B/Yb\nxjw8l0Gz/XisX1XyHdzrNKkkJCYwbss4ms9rzrBGw1j0zKIsn1QgnXbei8jrQLwx5uu0vH/UqFE3\nHvv7++Pv7++awJRS6jaMgfnjjlL8rT4MKnGMQj+sRWrWcDruL+d+IWBpADmy5WBX0C7KFCqTobGG\nhoYSGhqaodNMrbu6CKXdFbbCGOPrMKw70Bt43BgTZw8bBhhjzDj7+RpgpDHmlq4wvQilUsodzpw2\nLPr3DDpHvU5C/0Hc9+EwyJnzlvGMMczYM4PXN7zO8EbDeenRl8gm7j9zw5MuQnm3FYvgcNcaEWkJ\nvAo8lpRUbMuBr0TkY6wusAeBsLuctlJKucTWLw9jevWiVeGL5A8LIVfNak7HO3bxGL2W9+Lk5ZNs\n7L6RKvdVyeBIM4e7Odx4PrANeFhEfhORHsBEID/wnYjsEZHJAMaYfcAiYB/Wfpf+WpYopdzt2tVE\nljabSOWAOhR/vgVlft+WYlJZEL2AGtNqUKdkHbb33K5J5Tb0fixKqSzp5+BDXOjYk3tyJXL/6tkU\nebSS0/HO/HWGAasHEHEygrlPz/XY81I8qSvM/R2DSimVgcz1BLZ3eJ/CrRpw+amOPPLnphSTyuof\nV+M71Zd/5f8Xe4Ky5smOaaGXdFFKZRnntsRwqm0g2eLzcW5NGI2aV3A63sW4i/zfuv9j3c/r+PI/\nX9K0fNMMjjRz04pFKeX94uP5uccYEpv4E1krEL8/1/NgCkll85HN+E31IyExgch+kZpU0kArFqWU\nV7sWFs6pNj04GluCY3N388xzzs83uXr9Km9seIP5UfOZ1noabSq1yeBIvYcmFqWUd4qL4/RLY8g2\ncxqLfd7nuegXuPc+5/u29xzfQ7dvu1H5XuuSLPfluy+Dg/UumliUUl7H7AjjXPtAdp6qyOkx4Qwa\nWhJnl729nnid9za/x8SwiXzc4mOe9XkWvT7u3dPEopTyHleu8NfQkcTNmMv4+ycQEN6ZJ6s4TxQH\nTh/ghW9foPA9hdnTZw+lC5bO4GC9l+68V0p5hy1buPyQH+tnHeGTwEhGH+xCZSdJJdEk8sn2T2g0\nuxE9qvdg7fNrNam4mFYsSqnM7dIlrg8dweW5i3kl9yQ6L2/PqGbORz1y/gg9lvXg6vWr/NDzhyxz\nq+CMphWLUirz2rCBuEd8Wf31eV56Ipqxh9rTzElSMcbwefjn1J5Rm+YVm7O5x2ZNKulIKxalVOYT\nG4t5dSgXF62mr5lKy4mtmN0Npzvoj188Tp+VfTgSe4T13dbjd79fxsebxWjFopTKXIKDuV7Fh1Ur\nDZ0qR/NOeCteeOHWpGKMYX7UfKpPq45fCT929t6pSSWDaMWilMoczp6FIUO4tHoj3a/NpvawZqx6\nFbJnv3XUPy//Sd+VfTl45iCrnl1F7ZK1Mz7eLEwrFqWU51u6lMRqPqzdVgD/olGM2NCMYcOcJ5Vv\nYr7Bd4ovlYpVYnfQbk0qbqAVi1LKc506BYMG8dfmXbyQsICSLRuzeRzkyXPrqKf/Om1d3v5EBEu7\nLOXR0o9mfLwK0IpFKeWJjIGFCzE+PoT8WBqfhAh6z23Mp586Tyrf7v8W3ym+PFDwAfb22atJxc20\nYlFKeZbjx6F/f65GHiQw3zKulatH2FooVuzWUc9eOcug4EHs+GMH33T8hoZlGmZ8vOoWWrEopTyD\nMTB3LsbPjx8uVqVi7F5avFWPb75xnlRWHlqJzxQfiuUpRnifcE0qHkQrFqWU+x09Cn36cO3IMQaU\nWcv+qzXYshPKl7911PNXz/PSmpfYdGQT89vPp0m5Jhkfr7otrViUUu5jDEyfDjVrEpmvPmVP7qRC\nhxps3Og8qaz5aQ0+U3zIlzMfkf0iNal4KK1YlFLucfgw9OrF9fMXGV47hJXR1Vi5FmrVunXUC3EX\neGXdK6z7eR1z2s2hWYUULgamPIJWLEqpjJWYCBMnQp06/FSxBQ/9uY2rD1Zj927nSeX7X77Hd4ov\nAJH9IjWpZAJasSilMs6hQ9CzJ4nXExnXeiuTVlVi9mxo0eLWUS9du8TQ74ay4tAKZraZSYsHnYyk\nPFKaKxYRmSUiJ0Uk0mFYERFZJyIHRWStiBRyeG24iPwoIvtFpPndBq6UykQSEuCDD6BBA36v35Ea\nFzex62IlIiKcJ5WNv27Ed4ovV65fIapflCaVTEaMMWl7o0gj4BIw1xjjaw8bB5wxxowXkdeAIsaY\nYSJSBfgKqAOUBtYDDxknExcRZ4OVUplVTAwEBpKYNx9Ta81k5BcV+PBD6ObkasSXr11mxPcjWLx/\nMdNaT6P1w63dE3MmJCIYYzzivspprliMMVuAc8kGtwO+sB9/ATxtP24LLDDGXDfG/Ar8CNRN67SV\nUplAfDyMGQNNmnCydSCNr67nf3srsHs3Tq9GvOnIJqpPq87Zq2eJ6helSSUTc/U+luLGmJMAxpgT\nIlLcHl4K+MFhvD/sYUopbxQeDj16YEqU4MuX9/Dyx2V46y148UXIluzn7OVrlxn+/XD+t/9/TH5q\nMu0eaeeemJXLpPfO+zT1aY0aNerGY39/f/z9/V0UjlIqXcXFWVXKtGmcG/E+XYNf4MxSYcsWeOSR\nW0cPORxCz+U9aVy2MVH9oiiap2jGx5xJhYaGEhoa6u4wnErzPhYAESkLrHDYx7If8DfGnBSR+4EQ\nY0xlERkGGGPMOHu8NcBIY8wOJ5+p+1iUyozCwiAwECpWZGmLKQSNKsmLL8Lw4ZAz582jXoy7yGvr\nX2PFoRVMbTWVVg+3ck/MXsQr9rHYxP5Lshzobj8OAJY5DO8iIrlEpDzwIBB2l9NWSnmCK1dg6FBo\n25aLg9+gyz1LGT6xJKtXw1tv3ZpU1v+yHp8pPsRdjyOqX5QmFS+U5q4wEZkP+APFROQ3YCQwFvhG\nRAKBI0AnAGPMPhFZBOwD4oH+WpYo5QW2bLGqlBo1+P7jSAJeLU7HjjDn81svbx97NZZXv3uVNT+t\nYXqb6bR8sKVbQlbp7666wtKDdoUplQlcugQjRsDixVx5fxIvb27PmjUwZw40bXrr6Gt+WkPQiiCe\nfPBJ3m/+PgVzF8z4mL2cJ3WF6Zn3Sql/ZsMG6NULGjVix6xonhtYlMaNISICChW6edRzV84xZN0Q\nQg6HMLvdbL0cSxah1wpTSqVObCz06QMBAcR/NJHhpebydGBRPvjAqlSSJ5UVB1fgM8WHvDnyEtUv\nSpNKFqIVi1LqzoKDraTSsiXRC6J5rn8hype3qpTixW8e9eyVswxeM5htR7cx7z/zaFreSd+Y8mpa\nsSilUnb2LHTvDv37kzBjNmMrTKfp04UYMgS+/fbWpLL0wFKqTa5G0XuKEtk3UpNKFqUVi1LKuaVL\nYcAAaN+eX5ZF0a1ffnLnhl27oGzZm0c9/ddpBgYPZPex3SzquIhGZRq5J2blEbRiUUrd7NQp6NoV\nXn0V8/UCplWbSL0n8tOpE6xff2tSWbxvMT5TfCiZvyThfcM1qSitWJRSNmNg0SIYPBi6deP4mFkE\nvpiXU6dg0yaoXPnm0f+8/CcDVg8g6mQUSzotof4D9d0Tt/I4WrEopeD4cejQAUaPxixdxtc136d6\ng7zUqwc//HBzUjHG8HXU1/hO8aVC4Qrs7bNXk4q6iVYsSmVlxsC8efDKKxAUxKlPv6bfS7nZvx9W\nrYLatW8e/cSlE/Rb1Y9DZw6xvOty6pbSu1+oW2nFolRWdfQotGoFH30Ea9eypOYYfOvkpmJF2L37\n5qRijGFexDx8p/hS5d4q7Anao0lFpUgrFqWyGmNgxgx4/XUYNIizQcMYOCQnO3fCkiVQP1mv1m+x\nv9FnZR+OXzxO8HPB1CpZyz1xq0xDKxalspLDh6FZM5g5E0JCWFXzTXxr5eTee617czkmlUSTyOSd\nk6k1vRaNHmjEzt47NamoVNGKRamsIDERPvsMRo+GoUOJ7TmEIUNzsGEDfPklJL+X3qEzh+i1vBfx\nifFs7L6RKvdVcUvYKnPSikUpb3foEDRpAgsWwNatrK85FN+aOciZEyIjb04q1xOvM37reBrMakCH\nyh3Y0mOLJhX1j2nFopS3SkiAjz+GsWPhrbe4FDCAocOzs2KF1RPWosXNo0eciKDn8p4UvqcwO3vv\npHyR8u6JW2V6mliU8kYxMdYNuPLlg7AwNv1egR41oXFjiIqCwoX/HjXuehxjNo1h2u5pjG02lh7V\neyDiEbf1UJmUJhalvEl8PIwbBxMmwDvvcOX53rz+ZjYWLICpU6Ft25tH3/77dgKXBfJwsYcJ7xtO\nyQIl3RO38iqaWJTyFuHh0KMHlCgBe/aw43gZAmpB9epWlVKs2N+jXr52mTc2vMGCmAV80vITOlbp\nqFWKchndea9UZhcXB2++Cc2bw0svEbc0mOFTytCuHbz9trXP3jGpfP/L9/hM8eH0ldNE94umU9VO\nmlSUS2nFolRmFhZm7UupWBHCw9lzoiQBdeDBB62bcJUo8feo56+e59V1r7L257VMbT2Vpx56yn1x\nK6+mFYtSmdGVKzB0qLXT5I03iP9mKaNnlKRlS2vwkiU3J5XlB5dTbXI1cmbPSXT/aE0qKl1pxaJU\nZrNli1Wl1KgBkZFE/1mcgPrW3Rz37oVSpf4e9c/LfzIoeBC7j+/mq/Zf0aRcE/fFrbIMrViUyiwu\nX4ZBg6BTJxg7lutfLWTs7OI0bQr9+sHq1X8nFWMMX0V+hc8UH8oUKkNE3whNKirDpEvFIiIvAz2B\nRCAK6AHkAxYCZYFfgU7GmNj0mL5SXmfDBujVCxo1guhoDp4qSkAjyJsXdu6EcuX+HvVo7FH6rerH\nb7G/sbLrSuqUquO2sFXW5PKKRURKAgOBmsYYX6zk1RUYBqw3xlQCNgDDXT1tpbxObCz06QMBATBx\nIglz5vLR50Vp2BCef966VXBSUkk0iUzbNY2a02tSt1RddgXt0qSi3CK99rFkB/KJSCKQB/gDK5Ek\n1eJfAKFs15++AAAdP0lEQVRYyUYp5UxwsJVUWra0qpQThQh8DHLkgB07rAPBkvx09id6r+jNX/F/\nERoQStXiVd0Xt8ryXF6xGGOOAR8Cv2EllFhjzHqghDHmpD3OCaC4q6etlFc4exa6d4f+/WH2bBKm\nTOfDmYVo2BC6dIGQkL+TyvXE63yw7QMenfkobR9uy7bAbZpUlNu5vGIRkcJAO6x9KbHANyLyHGCS\njZr8+Q2jRo268djf3x//5Nf0VspbLV0KAwZA+/YQFcXBP/KnWKVEnIig14peFMhVgB29dlCxaMWU\nP1d5ndDQUEJDQ90dhlNiTIrb97R9oMgzQAtjTG/7eTfgUeBxwN8Yc1JE7gdCjDGVnbzfuDompTze\nqVPWEV+7dllVSoPGTJgA770HI0dauSab3b9wJf4Kb296m5l7ZupFI9UNIoIxxiNWhPQ43Pg34FER\nuUestf0JYB+wHOhujxMALEuHaSuVuRgDCxeCjw+ULg0RERws3pjGjWHFCqtKGTjw76Sy8deN+E31\n46ezPxHZL5LAGoGaVJTHcXlXmDEmTEQWA3uBePv/dKAAsEhEAoEjQCdXT1upTOX4casUOXAAli0j\noXa9G1XKqFHWLpakhHL+6nmGfjeU4J+CmfTkJNo90s6toSt1Oy7vCrtb2hWmvJ4xMG8evPIKBAXB\nm29y8Nfc9OgBuXLBrFk370tZsn8JA4MH0q5SO9574j0K3VPIfbErj+VJXWF6SRelMtLRo9YhxMeO\nwdq1JPjWuHGTx9GjrTPok6qUYxeP8eLqF9l3ah8LOiygcdnG7o1dqVTSS7oolRGMgenToWZNqF8f\nwsI4kKcGjRrBypXWRYqTdtAnmkSm756O31Q/qt5XlfC+4ZpUVKaiFYtS6e3wYetyLBcvQkgICZWr\npVilHDx9kKCVQcRdj2PDCxvwKeHj3tiVSgOtWJRKL4mJMHEi1KkDLVrAtm0cyFGNRo1g1aqbq5T4\nhHje3fwuDWc3pP0j7dkauFWTisq0tGJRKj0cOgQ9e1rJZetWEh6sxEcfWbejT16l7PxjJ71W9KJk\ngZLsCtpFucLl3Bq6UndLE4tSrpSQwI1+rrfeggEDOPBjdno0gnvusaqUChWsUS9fu8ybIW8yP2o+\nHzb/kGd9ntVzUpRX0K4wpVwlJgYaNLBujBIWRsKAQbz/UXYaNYJu3eD77/9OKmt/Wku1KdU49dcp\novpF8Zzvc5pUlNfQikWpuxUfD+PHw4QJMGYM9O7NgUPZ6NEI8uSx7pdSvrw16pm/zvDy2pfZdGQT\nU1tPpeWDLd0bu1LpQCsWpe5GeDjUrQubN8Pu3ST06sP7H2ajcWN44QXrfinly1t3dJwfNZ+qk6tS\nLE8xovtHa1JRXksrFqXSIi7Oqk6mTYP334cXXuDAQaF7Q+uujmFhf1cpv8X+duOOjsu6LKNe6Xru\njV2pdKYVi1L/VFgY1KoFkZEQHs715wJ4/wOhcWPrRo9JVUpCYgITd0yk5rSaNCjdgN1BuzWpqCxB\nKxalUuvKFesa9nPnWvtTOncmOkYIfBoKFLj53vMxf8bQa0UvcmbLyZbALTxy7yNuDV2pjKQVi1Kp\nsWUL+PnBkSMQGUl8hy68PUZo2hR69/773vNx1+MYGTIS/y/8CfALILR7qCYVleVoxaLU7Vy+DMOH\nw+LFMGkStG/P3r3QozmUKgV79sADD1ijbv1tK71W9KJSsUqE9wmnVMFS7o1dKTfRxKJUSjZssK7x\n1agRREcTl68ob78BM2bABx/A88+DiHWvlGHrh7H84HI+ffJTOlTuoOekqCxNu8KUSi421rq0fUCA\nda2vuXPZ8WNRatSAffsgIsI64REMC6MXUuWzKgjCvgH7eKbKM5pUVJanFYtSjoKDraTSsiVER/NX\nzkK89Qp89RV88gl07GhVKYfPHab/6v4cjT3K4k6LafBAA3dHrpTH0IpFKYCzZ6F7d+t+wLNnw/Tp\nbIoohJ+fdU+uyEjo1AmuJ8Yzfut46syoQ5OyTdjTZ48mFaWS0YpFqaVLrevXt28PUVFcIj/DXrQG\nf/YZtLNvL7/99+30WdmHEvlKsKPXDioWrXj7z1Uqi9LEorKuU6dg0CDYtQsWLIDGjVm/3jp8uGlT\niIqCIkUg9mosI74fwZIDS/iw+Yd0rdZV96ModRvaFaayHmNg4ULw8YHSpSEigljfxvTubd1CZepU\nqzescGHD4n2LqTK5CtcSrhHTP0Yvba9UKmjForKW48etbq8DB2DZMqhXj5UrrRtvtWljVSkFC8KR\n80cYsHoAv5z7hQUdFug955X6B7RiUVmDMdalWPz8oEoV2LuXMw/Wo1s3GDzYemnyZMib/zofbvuQ\nWtNr8WjpRwnvG65JRal/KF0qFhEpBMwEqgGJQCBwCFgIlAV+BToZY2LTY/pK3eToUesQ4mPHYO1a\nqFGD//0PBg6Ezp2tI77y5bNuERy0MoiieYryQ88feKjYQ+6OXKlMKb0qlk+A1caYyoAfcAAYBqw3\nxlQCNgDD02naSlmMgenToWZNqF8fwsI4WbIGHTvC669bV2n5+GNIyHGBwcGDafN1G4Y8OoT13dZr\nUlHqLrg8sYhIQaCxMWYOgDHmul2ZtAO+sEf7Anja1dNW6obDh6FZM5g5E0JCMG+8yfzFufD1hQcf\ntO7P1aABLD2wlKqTq3Lx2kVi+sfQza+b7pxX6i6lR1dYeeC0iMzBqlZ2AS8BJYwxJwGMMSdEpHg6\nTFtldYmJ1skno0fD0KEwZAh/nMxBv3bw66+wahXUrg1HY48ycOlADpw+wLz/zMO/nL+7I1fKa6RH\nYskB1AQGGGN2icjHWN1gJtl4yZ/fMGrUqBuP/f398ff3d32UyvscOmQdL5yYCFu3Yh6uxJw5MGyY\ndUL94sWQPUcCn2yfxNub3mZg3YEsfGYhuXPkdnfkSv1joaGhhIaGujsMp8SYFLfvaftAkRLAD8aY\nCvbzRliJpSLgb4w5KSL3AyH2Ppjk7zeujkl5uYQEa2fJ2LHw1lswYABHfs9O795w5gzMmQO+vrDn\n+B6CVgSRP1d+praeqvdJUV5FRDDGeEQ/rsv3sdjdXUdF5GF70BNADLAc6G4PCwCWuXraKguKibF2\nlqxeDWFhJL44iMnTslO7Njz+OOzYARUeucSQtUN48qsnebHui4QEhGhSUSodpdcJkoOAr0QkJ/AL\n0APIDiwSkUDgCNApnaatsoL4eBg/3rpF8Jgx0Ls3P/2SjV6PQ1wcbNoElSvDioMreDH4RfzL+RPd\nL5r78t3n7siV8nou7wq7W9oVpu4oPBx69IASJWD6dBJKleHTT+Gdd+CNN6zzU05c/oNBawYRdTKK\nqa2n8nj5x90dtVLpyqu7wpRKN3Fx8Oab0Ly5dbp8cDAxF8vQsKF1dZbt22HgoASm7J6E31Q/qt5X\nlch+kZpUlMpgeq0wlTmEhUFgIFSsCOHhXLu3JGPftm7waPeEEflnOA1m9yFX9lxs6rGJKvdVcXfU\nSmVJmliUZ7tyBUaOtC7mNWECdO5M2E6hZwsoVw727oUixS/z2vpRfBHxBe8+8S6BNQLJJlqMK+Uu\n+u1TnmvLFuuikUeOQGQkf7XtwiuvCm3bwogRsHw5RF5ZTbUp1Th+6ThR/aLoVbOXJhWl3EwrFuV5\nLl+G4cOtMxonTYL27QkJsbq76tWzLm1//Z7jdF48mN3HdzOt9TSaV2zu7qiVUjb9aac8y4YN1g24\nzp+H6GjOP96eoCAICLB6wuZ9mcjiX6fgO9WXB4s+SHS/aE0qSnkYrViUZ4iNta7ttXq1dQvHVq1Y\nvty6FEubNhAdDUeuRNFwdhDZJBshASFUK17N3VErpZzQxKLcLzjYul9Ky5YQHc2fcYUY1AX27IGv\nvoI6Df7ivxv/y6y9sxjTdAy9a/XW/ShKeTD9dir3OXsWune3ypLZszHTpjNveSF8fKBsWYiIgKul\n1lJtcjWOxB4hql8UfWr30aSilIfTikW5x9Kl1r3n27eHqCh+O5ufvq3gjz+sS9s/8MhJeq5+mR9+\n/4HJT03myYeedHfESqlU0p9+KmOdOgVdu8Krr8KCBSR+MpHJc/NTqxY0bAhhOxPZKzPwmeLDAwUf\nIKZ/jCYVpTIZrVhUxjAGFi2yLsXSrRvMmsXBo3np1cS6fcqmTZBYLIYnvuzD9cTrrH9hPb4lfN0d\ntVIqDbRiUenv+HHo0MG6q+OyZcS/+z5jP81Lo0bQuTOs23CFr46/gf8X/jzr8yxbA7dqUlEqE9PE\notKPMdalWPz8oEoV2LuXvbnqUbcuhITAzp3wyFPrqT7dl4NnDhLRN4L+dfqTPVt2d0eulLoL2hWm\n0sfRo9YhxMeOwdq1XK1cg9GjYNYseP99eLLDKf7vuyFsPrKZSU9NovXDrd0dsVLKRbRiUa5lDEyf\nDjVrQv36EBbGlss18PODn36CiAhDgu9sfKZWo0S+EkT3j9akopSX0YpFuc7hw9CrF1y4ACEhXCxb\njeFD4NtvrUt+VW58gC4r+/BX/F+seW4NNf5Vw90RK6XSgVYs6u4lJlo3RqlTB1q0gB9+IPhoNapV\ng6tXYVf4VcILj6TR7EY8U/kZtvfcrklFKS+mFYu6O4cOQc+eVnLZupXTxSrxcg/YutXan5K9Ygj+\nC/tS9b6qhPcNp3TB0u6OWCmVzrRiUWmTkAAffAANGkDHjpiNm1gYXgkfH7j3XgjZcZovL3UnYGkA\n45uNZ0nnJZpUlMoitGJR/1xMjHWb4Hz5ICyMY/dUoF8Ha+f8kiWGQ3nmUveLoXSt1pWY/jEUyF3A\n3RErpTKQViwq9eLj4Z13wN8fAgMx361n5oYKVK8O1avDwvWHeP3QE3wa9imrnl3FhJYTNKkolQVp\nxaJSJzwcevSAEiVg925+ji9D73/DpUsQvC6OVbHj8J/3KW889gYv1n2RHNl01VIqq0q3ikVEsonI\nHhFZbj8vIiLrROSgiKwVkULpNW3lQnFx8Oab0Lw5DB5MwspgPvymDPXqQatWMG7BJp7fWp1dx3ax\np88eXnr0JU0qSmVx6dkVNhjY5/B8GLDeGFMJ2AAMT8dpK1cIC4NatSAyEsLDiarVnfoNhJUrYc3G\ns+x/qBfdlj7LO4+/w7IuyyhTqIy7I1ZKeYB0SSwiUhp4CpjpMLgd8IX9+Avg6fSYtnKBK1es2wS3\nbQtvvEHcwqWMnFaSxx+HXr0MPT/5ijZrqpInRx72DdhH+8rtERF3R62U8hDp1WfxMfAq4NjdVcIY\ncxLAGHNCRIqn07TV3diyxTriq0YNiIxkx+Hi9KwNFSvCtxt/YvTufpzadoplXZZRt1Rdd0erlPJA\nLk8sItIKOGmMCRcR/9uMalJ6YdSoUTce+/v74+9/u49RLnH5MgwfDosXw6RJXG7RnjffhK+/hg8n\nXOPwvz7g6ZUfMazRMN2PopQHCA0NJTQ01N1hOCXGpLh9T9sHirwLPA9cB/IABYBvgdqAvzHmpIjc\nD4QYYyo7eb9xdUzqDjZssK7x1agRTJjA+j1FCQqy7ujYZegWhm7uQ7nC5fjsqc8oV7icu6NVSjkh\nIhhjPKJP2uWJ5aYPF2kC/J8xpq2IjAfOGGPGichrQBFjzDAn79HEklEuXLBuEbx6NUydyrkGrXjl\nFfjuO3h/0jk2yDBW/riSCS0m8EyVZ3Q/ilIezJMSS0aeIDkW+LeIHASesJ8rdwkOhmrVrMvcR0fz\n7bVWVKsGue8xjPrfAl7+sSrZs2Unpn8MHat21KSilEq1dK1Y0kIrlnR29iwMGQIbN8KMGZyo1oyB\nA60jiv878TBz/uzHsYvHmNZ6GvUfqO/uaJVSqZRVKxblbkuXgo8PFCiAiYzi89+b4ecHFR6Kp9u0\ncQwIr0PTck3ZHbRbk4pSKs300J6s4NQpGDQIdu2Cr7/m1zKPEdQBTp+G8fO38+GhIEr+XpKw3mFU\nKFLB3dEqpTI5rVi8mTGwcKFVpZQuTcKeCD7Z+xi1a0PDx2OpO7o/w8PbM6LxCIKfC9akopRyCa1Y\nvNXx4zBgABw4AMuWsa9APXo2h5y5DG8uWMz4yJdoTWti+sdQJE8Rd0erlPIiWrF4G2Ng7lzw84Mq\nVbi2Yy//XVuPJk2g9fO/UqBPG6YfGsXCZxYyrc00TSpKKZfTisWbHD0KffrAsWOwdi1h8TXo2QDK\nlLtOn7kT+DhqLEN8h/Bt5yXkyp7L3dEqpbyUVizewBiYPh1q1oT69bkcEsaQeTVo2xY6DwnjWOs6\nbD+1hu29tjOi8QhNKkqpdKUVS2Z3+LB1OZYLFyAkhO9PVqN3LajT6AJtJr/BpF8W8UHzD3jO5zk9\nyVEplSG0YsmsEhNh4kSoUwdatODc6h/o+XE1uvcwdP3vErZWr4LJ8Rcx/WN43vd5TSpKqQyjFUtm\ndOgQ9OxpJZetW1kSU4mB1aFZh9/wfXcgS2IPMb/DfB4r+5i7I1VKZUFasWQmCQnwwQfQoAF07MiJ\nRZt45vVKDH/9Oh3en8Cq0jV59IHahPcJ16SilHIbrVgyi5gY6wZc+fJhdoTx+aYKvFYDWgftJl+z\nIKLjC7Gt5zYeLvawuyNVSmVxehFKTxcfD+PHw4QJMGYMh5v1JqhvNk7FXqTawLf47uR8xjcbzwt+\nL+h+FKWyML0IpUqd8HCoWxc2byYhbDcTrvShTr1slHp8OWe7ViVHgXPE9I8hoHqAJhWllMfQisUT\nxcXBO+/A1KkwfjwxtQPo2UuQgn9QoPNAfv0rmmmtp9G0fFN3R6qU8hBasaiUhYVBrVoQEcG1sHBG\nH+lOk8cTKddpIj/9uzr1K/gQ2S9Sk4pSymPpzntPceUKjBxpXedrwgR2lOtMz9ZCsarhlB4ZxPF8\nedjcejOP3PuIuyNVSqnb0orFE2zZYl008sgRLv8QyZCwLrTt+Bfl+rzCgTotGNSgLyEBIZpUlFKZ\nglYs7nT5MgwfDosXw6RJrC/YnqAnoOy/V5H75QEUKd2YqOZRFM9X3N2RKqVUqmlicZcNG6xrfDVq\nxPkt0QwZU5S1245R4aXBHGUvs1vPpFmFZu6OUiml/jHtCstoFy5Yl7YPCICJE/lfu7lUeawQvxSb\nTFygH02qVSKqX5QmFaVUpqUVS0YKDraSSsuWnFgfzYARhdjzRxRFXwniesHsbGwdStXiVd0dpVJK\n3RVNLBnh7FkYMgQ2bsTMms3cY814xf8KD/ceweVHZzK88Rh61exFNtECUimV+bl8SyYipUVkg4jE\niEiUiAyyhxcRkXUiclBE1opIIVdP2yMtXQo+PlCgAEdXR/HUR814+6v13DPEh9K+PxPRN4KgWkGa\nVJRSXsPlZ96LyP3A/caYcBHJD+wG2gE9gDPGmPEi8hpQxBgzzMn7vePM+1OnYNAg2LWLxBmzmHHw\nMUaMOcUDvYdwvuBmJreazFMPPeXuKJVSXsKrz7w3xpwwxoTbjy8B+4HSWMnlC3u0L4CnXT1tj2AM\nLFxoVSmlS/PLtxE88d/GjFv7OQyoRrNHSxDTP0aTilLKa6XrPhYRKQdUB7YDJYwxJ8FKPiLifSdn\nHD8OAwbAgQMkLFnGpJ31GPX0IYp170vR4hdY3CaYmv+q6e4olVIqXaVbYrG7wRYDg40xl0Qkef9W\niv1do0aNuvHY398ff3//9AjRdYyBefPglVcgKIiDo76me3/heIUx0GsCA5q8zsB6A8mRTY+VUEq5\nRmhoKKGhoe4Ow6l0ubqxiOQAVgLBxphP7GH7AX9jzEl7P0yIMaayk/dmrn0sR49ahxAfO8b1GXP4\ncEMN3v1yK3k7B1GrQgU+e2oSZQuXdXeUSikv59X7WGyzgX1JScW2HOhuPw4AlqXTtDOGMTBjBtSs\nCfXrEz07jNoDyzPxcF/ydOvExPajWdF1uSYVpVSWkx5HhTUENgFRWN1dBhgBhAGLgAeAI0AnY8x5\nJ+/3/Irl8GHo3RtiY4mfPod3llflo+DFZG/1Ep1rtGVss/cofE9hd0eplMpCPKli0Rt9/ROJifDZ\nZzB6NAwdyu4mQ3hu8B+cqz+AwuUOM/s/02lYpqG7o1RKZUGelFh0b3JqHToEPXtCYiLXQrby1tcP\n8tkbE6HVGF5r8jJDGy4hV/Zc7o5SKaXcThPLnSQkwMcfw9ix8NZb7Kk/gM5993OmYQN8O+Xl8/Y/\n8FCxh9wdpVJKeQxNLLcTEwOBgZAvH/Fbwxj9dSk+Hv1fsreczAdPvUuvmj31UixKKZWMJhZn4uNh\n/HiYMAHGjCHq0d480z+MY7Xb0KjDg8zuEE6pgqXcHaVSSnkkTSzJxcaCvz+UKMH1Hbt5++uijB81\nhNz+C5n5n0/oVLUjIh6xf0wppTySHhWWnDGwYQP7//U4/3n1O4749qFV1ceY1v4jiuUt5r64lFLq\nNvSoMA+WkCi8u7cG78zoQb5HQ/m2yzRaPtTC3WEppVSmoYklmRPnY3nvvC9d/tOBSf+JJn+u/O4O\nSSmlMhXtCnPiyPkjeikWpVSm4kldYZpYlFLKC3hSYtGTMJRSSrmUJhallFIupYlFKaWUS2liUUop\n5VKaWJRSSrmUJhallFIupYlFKaWUS2liUUop5VKaWJRSSrmUJhallFIupYlFKaWUS2liUUop5VIZ\nmlhEpKWIHBCRQyLyWkZOWymlVMbIsMQiItmASUALoCrQVUQeyajpe4rQ0FB3h5CuvLl93tw20PYp\n18nIiqUu8KMx5ogxJh5YALTLwOl7BG9fub25fd7cNtD2KdfJyMRSCjjq8Px3e5hSSikvojvvlVJK\nuVSG3UFSRB4FRhljWtrPhwHGGDMu2Xh6+0illEoDT7mDZEYmluzAQeAJ4DgQBnQ1xuzPkACUUkpl\niBwZNSFjTIKIvAisw+qCm6VJRSmlvE+GVSxKKaWyBo/Zee+NJ0+KyCwROSkikQ7DiojIOhE5KCJr\nRaSQO2NMKxEpLSIbRCRGRKJEZJA93Fval1tEdojIXrt9I+3hXtE+sM4tE5E9IrLcfu5NbftVRCLs\n5RdmD/Om9hUSkW9EZL/9HaznSe3ziMTixSdPzsFqk6NhwHpjTCVgAzA8w6NyjevAEGNMVaA+MMBe\nZl7RPmNMHNDUGFMDqA48KSJ18ZL22QYD+xyee1PbEgF/Y0wNY0xde5g3te8TYLUxpjLgBxzAk9pn\njHH7H/AoEOzwfBjwmrvjclHbygKRDs8PACXsx/cDB9wdo4vauRRo5o3tA/ICu4A63tI+oDTwHeAP\nLLeHeUXb7PgPA8WSDfOK9gEFgZ+dDPeY9nlExULWOnmyuDHmJIAx5gRQ3M3x3DURKYf1q3471ort\nFe2zu4r2AieA74wxO/Ge9n0MvAo47mT1lraB1a7vRGSniPSyh3lL+8oDp0Vkjt2VOV1E8uJB7fOU\nxJKVZeqjJ0QkP7AYGGyMucSt7cm07TPGJBqrK6w0UFdEquIF7RORVsBJY0w4cLvzHjJd2xw0NMbU\nBJ7C6qZtjBcsO1sOoCbwmd3Gy1i9PB7TPk9JLH8AZRyel7aHeaOTIlICQETuB/50czxpJiI5sJLK\nPGPMMnuw17QviTHmAhAKtMQ72tcQaCsivwBfA4+LyDzghBe0DQBjzHH7/ymsbtq6eMeyA6tH56gx\nZpf9/H9YicZj2ucpiWUn8KCIlBWRXEAXYLmbY3IV4eZfhcuB7vbjAGBZ8jdkIrOBfcaYTxyGeUX7\nROTepKNqRCQP8G9gP17QPmPMCGNMGWNMBazv2gZjTDdgBZm8bQAikteupBGRfEBzIAovWHYAdnfX\nURF52B70BBCDB7XPY85jEZGWWEc6JJ08OdbNId01EZmPtXO0GHASGIn16+kb4AHgCNDJGHPeXTGm\nlYg0BDZhfWGN/TcC64oKi8j87fMBvsBaH7MBC40x74hIUbygfUlEpAnwf8aYtt7SNhEpD3yLtU7m\nAL4yxoz1lvYBiIgfMBPICfwC9ACy4yHt85jEopRSyjt4SleYUkopL6GJRSmllEtpYlFKKeVSmliU\nUkq5lCYWpZRSLqWJRSmllEtpYlFKKeVSmliUUkq51P8DouL4bh+ZMOYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x169bb9a50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(omegas3,k_o_vs_omega(omegas3))\n",
    "plt.plot(omegas3,k_e_vs_omega(omegas3))\n",
    "plt.plot(omegas3,omegas3*2.3)\n",
    "plt.title('Linear relationship, ordinary dispersion, extraordinary dispersion')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "######"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dispersion_difference_function = k_symb(omega1,pol='o')+k_symb(omega2,pol='o')+k_symb(omega3,pol='e')\n",
    "\n",
    "symbols=[omega1,omega2]\n",
    "update_expression = make_newton_update_expression(dispersion_difference_function.expand(),\n",
    "                                                    symbols=[omega1,omega2,omega3])\n",
    "\n",
    "update_expression = [ex.subs(omega3,-omega1-omega2) for ex in update_expression]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "newton_update_func_0 = ufuncify( symbols , update_expression[0])\n",
    "newton_update_func_1 = ufuncify( symbols , update_expression[1])\n",
    "newton_update_func_2 = ufuncify( symbols , update_expression[2])\n",
    "\n",
    "omegas_dict = {}\n",
    "\n",
    "ind_range = range(2,5)\n",
    "for i in ind_range:\n",
    "    omegas_dict[i] = np.asarray([6.+ j*5*pow(10,-i) for j in range(1+pow(10,i))])\n",
    "\n",
    "def trace_omegas(omegas,arg2 = 24.):\n",
    "    fixed_omega_1_init=omegas[0]\n",
    "    omegas_delta = omegas[1]-omegas[0]\n",
    "    oms = [(fixed_omega_1_init, arg2)]\n",
    "    converged_points = []\n",
    "    for i in range(int(1e5)):\n",
    "        arg2 += newton_update_func_1(*(oms[-1]))\n",
    "        oms.append((omegas[0],arg2))\n",
    "        if sum([ (oms[-1][i] - oms[-2][i])**2 for i in range(2)]) < omegas_delta * 1e-15:\n",
    "            break\n",
    "    for arg1 in omegas:\n",
    "        for i in range(50):\n",
    "            for j in range(10):\n",
    "                arg2 += newton_update_func_1(*(oms[-1]))\n",
    "                oms.append((arg1,arg2))\n",
    "            if sum([ (oms[-1][k] - oms[-2][k])**2 for k in range(2)]) < omegas_delta * 1e-15:\n",
    "                break\n",
    "        converged_points.append(arg2)\n",
    "    return np.asarray(converged_points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "i = 2\n",
      "'timed_trace_omegas' 0.50 sec\n",
      "i = 3\n",
      "'timed_trace_omegas' 0.52 sec\n",
      "i = 4\n",
      "'timed_trace_omegas' 4.87 sec\n",
      "'interpolate_omegas' 0.01 sec\n"
     ]
    }
   ],
   "source": [
    "@timeit\n",
    "def timed_trace_omegas(*args,**kw):\n",
    "    return trace_omegas(*args, **kw)\n",
    "\n",
    "@timeit\n",
    "def interpolate_omegas(omegas_few,omegas_many,omegas_few_values):\n",
    "    f = interpolate.UnivariateSpline(omegas_few, omegas_few_values,k=2,s=1e-10)\n",
    "    return f(omegas_many)\n",
    "\n",
    "omegas_2_converged = {}\n",
    "for i in ind_range:\n",
    "    print \"i = \" + str(i)\n",
    "    omegas_2_converged[i] = timed_trace_omegas(omegas_dict[i])\n",
    "\n",
    "omega_values_interpolated = interpolate_omegas(omegas_dict[3],\n",
    "    omegas_dict[4],omegas_2_converged[3] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def trace_omegas_BENCHMARK(omegas,arg2 = 24.):\n",
    "    '''\n",
    "    Trace through triples of frequencies constrained by \n",
    "    :math:`omega1+omega2+omega3 = 0`.\n",
    "    \n",
    "    At each step of incrementing omega1, find the new omega2 \n",
    "    and omega3 by applying Newton's method to omega2 \n",
    "    and the constraint to obtain omega3.\n",
    "    \n",
    "    At each step apply 200 Newton's steps to converge.\n",
    "    \n",
    "    Args:\n",
    "        omegas (list of floats):\n",
    "            The frequencies in the range desired\n",
    "        fixed_omega_3 (optional [float]): initial value of omega3.\n",
    "        \n",
    "    Returns:\n",
    "        oms (list of floats): A history of the Newton's steps\n",
    "        converged_points (list of tuples): \n",
    "    '''\n",
    "    fixed_omega1=omegas[0]\n",
    "    omegas_delta = omegas[1]-omegas[0]\n",
    "    oms = [(fixed_omega1, arg2)]  ## k is accurate for omega inputs between 6-60.\n",
    "    converged_points = []\n",
    "    for i in range(int(1e5)):\n",
    "        arg2 += newton_update_func_1(*(oms[-1]))\n",
    "        oms.append((omegas[0],arg2))\n",
    "        if sum([ (oms[-1][i] - oms[-2][i])**2 for i in range(2)]) < omegas_delta * 1e-10 :\n",
    "            break\n",
    "    for arg1 in omegas:\n",
    "        for i in range(200):\n",
    "            arg2 += newton_update_func_1(*(oms[-1]))\n",
    "            oms.append((arg1,arg2))\n",
    "        converged_points.append((arg1,arg2))\n",
    "    return oms,converged_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def trace_omegas(omegas,arg2 = 24.):\n",
    "    '''\n",
    "    Trace through triples of frequencies constrained by \n",
    "    :math:`omega1+omega2+omega3 = 0`.\n",
    "    \n",
    "    At each step of incrementing omega1, find the new omega2 \n",
    "    and omega3 by applying Newton's method to omega2 \n",
    "    and the constraint to obtain omega3.\n",
    "    \n",
    "    \n",
    "    Trace through the omegas, when finding the next\n",
    "    step with Newton's method occasionally check the error \n",
    "    and allow terminating early.\n",
    "    \n",
    "    Args:\n",
    "        omegas (list of floats):\n",
    "            The frequencies in the range desired\n",
    "        fixed_omega_3 (optional [float]): initial value of omega3.\n",
    "        \n",
    "    Returns:\n",
    "        oms (list of floats): A history of the Newton's steps\n",
    "        converged_points (list of tuples): \n",
    "    '''\n",
    "    fixed_omega1=omegas[0]\n",
    "    omegas_delta = omegas[1]-omegas[0]\n",
    "    oms = [(fixed_omega1, arg2)]  ## k is accurate for omega inputs between 6-60.\n",
    "    converged_points = []\n",
    "    for i in range(int(1e5)):\n",
    "        arg2 += newton_update_func_1(*(oms[-1]))\n",
    "        oms.append((omegas[0],arg2))\n",
    "        if sum([ (oms[-1][i] - oms[-2][i])**2 for i in range(2)]) < omegas_delta * 1e-10 :\n",
    "            break\n",
    "    for arg1 in omegas:\n",
    "        for i in range(50):\n",
    "            for j in range(10):\n",
    "                arg2 += newton_update_func_1(*(oms[-1]))\n",
    "                oms.append((arg1,arg2))\n",
    "            if sum([ (oms[-1][k] - oms[-2][k])**2 for k in range(2)]) < omegas_delta * 1e-10:\n",
    "                break\n",
    "        converged_points.append((arg1,arg2))\n",
    "    return oms,converged_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "converged_dict = {}\n",
    "converged_dict_BENCHMARK = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time to converge:  1.592868\n",
      "time to converge:  0.80961\n",
      "time to converge:  6.283256\n"
     ]
    }
   ],
   "source": [
    "for ind,omegas in zip([3,4,5],[omegas3,omegas4,omegas5]):\n",
    "    t0 = time.clock()\n",
    "    oms,converged_points = trace_omegas(omegas)\n",
    "    t1 = time.clock()\n",
    "    print \"time to converge: \", t1-t0\n",
    "    k_os = k_o_func(omegas)\n",
    "    k_es = k_e_func(omegas)\n",
    "    converged_omegas = [(arg1,arg2,30-arg1-arg2) for arg1,arg2 in converged_points]\n",
    "    converged_dict[ind] = omegas,converged_omegas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time to converge:  0.445815\n",
      "time to converge:  2.811241\n",
      "time to converge:  28.655134\n"
     ]
    }
   ],
   "source": [
    "for ind,omegas in zip([3,4,5],[omegas3,omegas4,omegas5]):\n",
    "    t0 = time.clock()\n",
    "    oms,converged_points = trace_omegas_BENCHMARK(omegas)\n",
    "    t1 = time.clock()\n",
    "    print \"time to converge: \", t1-t0\n",
    "    k_os = k_o_func(omegas)\n",
    "    k_es = k_e_func(omegas)\n",
    "    converged_omegas_BENCHMARK = [(arg1,arg2,30-arg1-arg2) for arg1,arg2 in converged_points]\n",
    "    converged_dict_BENCHMARK[ind] = omegas,converged_omegas_BENCHMARK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIkAAAJPCAYAAAD404TLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Wm4ZQV95/vf2tOZai6KQuYZZJbJAYVKEMRZ1Dad7qRj\n7CFD981wTaKJaUUTk2B3hh7SGTrTjTe5GqM4tAMIWDgigigzIlDMVUDNp86wp3Vf7FNFoQVS1Kmz\nz/D5PM9+1rDX2vtfPOcN32ettYuyLAMAAADAwlbp9wAAAAAA9J9IBAAAAIBIBAAAAIBIBAAAAEBE\nIgAAAAAiEgEAAACQaYpERVEsLYriY0VR3FkUxe1FUby4KIrlRVFcVRTF3UVRXFkUxdLp+C4AAAAA\npt90XUn035J8rizLFyY5PcldSd6d5OqyLE9Icm2S35ym7wIAAABgmhVlWe7bBxTFkiQ3l2V5zA/s\nvyvJBWVZbiiK4qAka8uyPHGfvgwAAACA/WI6riQ6KsmTRVH8bVEU3y6K4i+LohhOsrosyw1JUpbl\n+iQHTsN3AQAAALAfTEckqiU5M8mflmV5ZpId6d1q9oOXKO3bJUsAAAAA7De1afiMh5M8VJbljVPb\nH08vEm0oimL1brebPb6nk4uiEI8AAAAApllZlsXeHL/PkWgqAj1UFMXxZVl+L8mFSW6fer09yeVJ\nfibJp57lM/Z1DJjVLrvsslx22WX9HgP2K3/nLAT+zlkI/J2zEPg7ZyEoir3qQ0mm50qiJPmlJP9Q\nFEU9yX1JfjZJNck/FUXxjiQPJHnbNH0XAAAAANNsWiJRWZbfTXLOHt565XR8PgAAAAD713Q8uBr4\nEdasWdPvEWC/83fOQuDvnIXA3zkLgb9z2LOi388DKoqi7PcMAAAAAPNJURR7/eBqVxIBAAAAIBIB\nAAAAIBIBAAAAEJEIAAAAgIhEAAAAAEQkAgAAACAiEQAAAAARiQAAAACISAQAAABARCIAAAAAIhIB\nAAAAEJEIAAAAgIhEAAAAAEQkAgAAACAiEQAAAAARiQAAAACISAQAAABARCIAAAAAIhIBAAAAs9zH\nvnRXXviT78m/vOwf+z3KvCYSAQAAALPOzjBUPfrkvO1NL81jj9yTk486tN9jzWtFWZb9HaAoyn7P\nAAAAAPTfx750V977lx/O9775yXQ3Ppqlp1+Ut7z2J/JH//F1WbpooN/jzSlFUaQsy2Kvzul3oBGJ\nAAAAYOEShvYPkQgAAACY9YSh/U8kAgAAAGYlYWhmiUQAAADArCEM9Y9IBAAAAPSVMDQ7iEQAAADA\njBOGZh+RCAAAAJgRwtDsJhIBAAAA+40wNHeIRAAAAMC0EobmJpEIAAAA2GfC0NwnEgEAAADPizA0\nv4hEAAAAwHMmDM1fIhEAAADwrIShhUEkAgAAAH6IMLTwiEQAAABAEmFooROJAAAAYAEThthJJAIA\nAIAFRhhiT0QiAAAAWACEIX4UkQgAAADmKWGIvSESAQAAwDwiDPF8iUQAAAAwxwlDTAeRCAAAAOYg\nYYjpJhIBAADAHCEMsT+JRAAAADCLCUPMFJEIAAAAZhlhiH4QiQAAAGAWEIboN5EIAAAA+kQYYjYR\niQAAAGAGCUPMViIRAAAA7GfCEHOBSAQAAAD7gTDEXNO3SFQUxbokW5N0k7TKsjy3KIrlST6a5Igk\n65K8rSzLrXs4VyQCAABg1hGGmMv6GYnuS3JWWZabd9t3eZKNZVl+qCiKdyVZXpblu/dwrkgEAADA\nrCAMMV/0MxLdn+Tssiw37rbvriQXlGW5oSiKg5KsLcvyxD2cKxIBAADQN8IQ81G/ryTakqST5C/K\nsvyroig2l2W5fLdjNpVluWIP54pEAAAAzChhiPnu+USi2jR993llWT5WFMWqJFcVRXF3kh8sP0oQ\nAAAAfbOnMPT2n/uAMARTpiUSlWX52NTyiaIoPpnk3CQbiqJYvdvtZo8/0/mXXXbZrvU1a9ZkzZo1\n0zEWAAAAC5wwxEKxdu3arF27dp8+Y59vNyuKYjhJpSzL0aIoRpJcleT9SS5Msqksy8s9uBoAAICZ\n4lYy6NMziYqiOCrJFendTlZL8g9lWf5BURQrkvxTksOSPJDkbWVZbtnD+SIRAAAA+0QYgqfr24Or\n94VIBAAAwPMhDMEzE4kAAACY14QheG5EIgAAAOYdYQj2nkgEAADAvCAMwb4RiQAAAJizhCGYPiIR\nAAAAc4owBPuHSAQAAMCsJwzB/icSAQAAMCsJQzCzRCIAAABmDWEI+kckAgAAoK+EIZgdRCIAAABm\nnDAEs49IBAAAwIwQhmB2E4kAAADYb4QhmDtEIgAAAKaVMARzk0gEAADAPhOGYO4TiQAAAHhehCGY\nX0QiAAAAnjNhCOYvkQgAAIBnJQzBwiASAQAA8EOEIVh4RCIAAACSCEOw0IlEAAAAC5gwBOwkEgEA\nACwwwhCwJyIRAADAAiAMAT+KSAQAADBPCUPA3hCJAAAA5hFhCHi+RCIAAIA5ThgCpoNIBAAAMAcJ\nQ8B0E4kAAADmCGEI2J9EIgAAgFlMGAJmikgEAAAwywhDQD+IRAAAALOAMAT0m0gEAADQJ8IQMJuI\nRAAAADNIGAJmK5EIAABgPxOGgLlAJAIAANgPhCFgrhGJAAAApokwBMxlIhEAAMA+EIaA+UIkAgAA\n2EvCEDAfiUQAAADPgTAEzHciEQAAwDMQhoCFRCQCAADYjTAELFQiEQAAsOAJQwAiEQAAsEAJQwBP\nJxIBAAALhjAE8MxEIgAAYF4ThgCeG5EIAACYd4QhgL0nEgEAAPOCMASwb0QiAABgzhKGAKaPSAQA\nAMwpwhDA/iESAQAAs54wBLD/iUQAAMCsJAwBzCyRCAAAmDWEIYD+EYkAAIC+EoYAZgeRCAAAmHHC\nEMDs09dIVBRFJcmNSR4uy/INRVEsT/LRJEckWZfkbWVZbt3DeSIRAADMMcIQwOzW70j0q0nOSrJk\nKhJdnmRjWZYfKoriXUmWl2X57j2cJxIBAMAcIAwBzB19i0RFURya5G+TfDDJ/z0Vie5KckFZlhuK\nojgoydqyLE/cw7kiEQAAzFLCEMDc9HwiUW2avvuPk/x6kqW77VtdluWGJCnLcn1RFAdO03cBAAD7\n0Z7C0Nt/7gPCEMA8t8+RqCiK1ybZUJbld4qiWPMsh7pcCAAAZilhCIDpuJLovCRvKIriNUmGkiwu\niuLDSdYXRbF6t9vNHn+mD7jssst2ra9ZsyZr1qyZhrEAAIBnIwwBzB9r167N2rVr9+kzpu3B1UlS\nFMUFSd459UyiD6X34OrLPbgaAABmB88YAlgY+vlMoj35gyT/VBTFO5I8kORt+/G7AACAZ+CKIQCe\ni2m9kuh5DeBKIgAAmHauGAJY2J7PlUQiEQAAzBPCEAA7iUQAALDACEMA7IlIBAAAC4AwBMCPIhIB\nAMA8JQwBsDdEIgAAmEeEIQCeL5EIAADmOGEIgOkgEgEAwBwkDAEw3UQiAACYI4QhAPYnkQgAAGYx\nYQiAmSISAQDALCMMAdAPIhEAAMwCwhAA/SYSAQBAnwhDAMwmIhEAAMwgYQiA2UokAgCA/UwYAmAu\nEIkAAGA/EIYAmGtEIgAAmCbCEABzmUgEAAD7QBgCYL4QiQAAYC8JQwDMRyIRAAA8B8IQAPOdSAQA\nAM9AGAJgIRGJAABgN8IQAAuVSAQAwIInDAGASAQAwAIlDAHA04lEAAAsGMIQADwzkQgAgHlNGAKA\n50YkAgBg3hGGAGDviUQAAMwLwhAA7BuRCACAOUsYAoDpIxIBADCnCEMAsH+IRAAAzHrCEADsfyIR\nAACzkjAEADNLJAIAYNYQhgCgf0QiAAD6ShgCgNlBJAIAYMYJQwAw+4hEAADMCGEIAGY3kQgAgP1G\nGAKAuUMkAgBgWglDADA3iUQAAOwzYQgA5j6RCACA50UYAoD5RSQCAOA5E4YAYP4SiQAAeFbCEAAs\nDCIRAAA/RBgCgIVHJAIAIIkwBAALnUgEALCACUMAwE4iEQDAAiMMAQB7IhIBACwAwhAA8KOIRAAA\n85QwBADsDZEIAGAeEYYAgOdLJAIAmOOEIQBgOohEAABzkDAEAEw3kQgAYI4QhgCA/UkkAgCYxYQh\nAGCmiEQAALOMMAQA9ENfIlFRFANJvpykkaSW5J/Lsnx/URTLk3w0yRFJ1iV5W1mWW/dwvkgEAMwr\nwhAA0G99u5KoKIrhsizHiqKoJvlakl9K8pYkG8uy/FBRFO9Ksrwsy3fv4VyRCACY84QhAGA26fvt\nZkVRDKd3VdEvJPlwkgvKstxQFMVBSdaWZXniHs4RiQCAOUkYAgBmq35eSVRJclOSY5L8aVmWv1kU\nxeayLJfvdsymsixX7OFckQgAmDOEIQBgLng+kag2HV9clmU3yYuKoliS5IqiKE5O8oPlRwkCAOak\nPYWht//cB4QhAGBemZZItFNZltuKolib5JIkG4qiWL3b7WaPP9N5l1122a71NWvWZM2aNdM5FgDA\nXhOGAIC5ZO3atVm7du0+fcZ0/LrZAUlaZVluLYpiKMmVSf4gyQVJNpVlebkHVwMAc4FbyQCA+aIv\nzyQqiuLUJP9PksrU66NlWX6wKIoVSf4pyWFJHkjytrIst+zhfJEIAOgbYQgAmI/6/utmz4dIBADM\nNGEIAJjvRCIAgGcgDAEAC4lIBACwG2EIAFioRCIAYMEThgAARCIAYIEShgAAnk4kAgAWDGEIAOCZ\niUQAwLwmDAEAPDciEQAw7whDAAB7TyQCAOYFYQgAYN+IRADAnCUMAQBMH5EIAJhThCEAgP1DJAIA\nZj1hCABg/xOJAIBZSRgCAJhZIhEAMGsIQwAA/SMSAQB9JQwBAMwOIhEAMOOEIQCA2UckAgBmhDAE\nADC7iUQAwH4jDAEAzB0iEQAwrYQhAIC5SSQCAPaZMAQAMPeJRADA8yIMAQDMLyIRAPCcCUMAAPOX\nSAQAPCthCABgYRCJAIAfIgwBACw8IhEAkEQYAgBY6EQiAFjAhCEAAHYSiQBggRGGAADYE5EIABYA\nYQgAgB9FJAKAeUoYAgBgb4hEADCPCEMAADxfIhEAzHHCEAAA00EkAoA5SBgCAGC6iUQAMEcIQwAA\n7E8iEQDMYsIQAAAzRSQCgFlGGAIAoB9EIgCYBYQhAAD6TSQCgD4RhgAAmE1EIgCYQcIQAACzlUgE\nAPuZMAQAwFwgEgHAfiAMAQAw14hEADBNhCEAAOYykQgA9oEwBADAfCESAcBeEoYAAJiPRCIAeA6E\nIQAA5juRCACegTAEAMBCIhIBwG6EIQAAFiqRCIAFTxgCAACRCIAFShgCAICnE4kAWDCEIQAAeGYi\nEQDzmjAEAADPjUgEwLwjDAEAwN4TiQCYF4QhAADYNyIRAHOWMAQAANNHJAJgThGGAABg/+hLJCqK\n4tAkf59kdZJukv9dluV/L4pieZKPJjkiybokbyvLcusezheJABYQYQgAAPa/fkWig5IcVJbld4qi\nWJTkpiRvTPKzSTaWZfmhoijelWR5WZbv3sP5IhHAPCcMAQDAzJoVt5sVRfHJJP9z6nVBWZYbpkLS\n2rIsT9zD8SIRwDwkDAEAQP/0PRIVRXFkkrVJTknyUFmWy3d7b1NZliv2cI5IBDBPCEMAADA7PJ9I\nVJvGL1+U5J+T/HJZlqNFUfxg+VGCAOahPYWht//cB4QhAACYY6YlEhVFUUsvEH24LMtPTe3eUBTF\n6t1uN3v8mc6/7LLLdq2vWbMma9asmY6xANhPhCEAAJhd1q5dm7Vr1+7TZ0zL7WZFUfx9kifLsvy/\nd9t3eZJNZVle7sHVAHOfW8kAAGDu6Nevm52X5MtJbk3vlrIyyW8luSHJPyU5LMkDSd5WluWWPZwv\nEgHMUsIQAADMTX1/cPXzIRIBzC7CEAAAzH0iEQDPizAEAADzi0gEwHMmDAEAwPwlEgHwrIQhAABY\nGEQiAH6IMAQAAAuPSARAEmEIAAAWOpEIYAEThgAAgJ1EIoAFRhgCAAD2RCQCWACEIQAA4EcRiQDm\nKWEIAADYGyIRwDwiDAEAAM+XSAQwxwlDAADAdBCJAOYgYQgAAJhuIhHAHCEMAQAA+5NIBDCLCUMA\nAMBMEYkAZhlhCAAA6AeRCGAWEIYAAIB+E4kA+kQYAgAAZhORCGAGCUMAAMBsJRIB7GfCEAAAMBeI\nRAD7gTAEAADMNSIRwDQRhgAAgLlMJALYB8IQAAAwX4hEAHtJGAIAAOYjkQjgORCGAACA+U4kAngG\nwhAAALCQiEQAuxGGAACAhUokAhY8YQgAAEAkAhYoYQgAAODpRCJgwRCGAAAAnplIBMxrwhAAAMBz\nIxIB844wBAAAsPdEImBeEIYAAAD2jUgEzFnCEAAAwPQRiZiX2s3xPHz/V3P/uq/miS3rsmnH+mye\n3JjNza3Z2tmRZtlOq+ykWXbSSjetspNuklqK1IpKaqnsWo5UBrK4OpzF9ZEsri/J4sbSLB1amRWL\nD8lBq0/JwYedmyUrjkylWu33P3tBEIYAAAD2D5GIOWtix6bc/p2P5pZ7r8ydm2/PXZOP5YHueDZ0\nu9nYTpZUkwMrRVZUallS1LK0GMjS6lCWVkcyWB1Io1JPo9pIvdLIQGUgRVGk3W33XmUr7W47rW4r\no+2xjHbGs707kdGymdFuK6NlO9vSzaZuN1s6SbdMlteSFUUlKyrVHFAM5OD60hw2eFCOWHZcjjro\nzBx97IVZdfDpYtLzIAwBAADsfyIRc0K308ltN304a2/523xzy235TntL7m11c2A1ObrayAn1lTlx\n8TE5YfVZOfqw83Lk8RdmcGTFjM23beP9efTBG/PY47dl/aZ7sn70oTy045E83NqUR7tjeazbzoZO\nmU6ZrK4lh1TqObK6OMcMviDHL39hTjjigrzwlEuzaNkhMzbzbCcMAQAAzCyRiFlr/QM35Ipr3pOr\nN96Qr7W2pVsmZ9eHc87w0TnvyFfnZS//5TkXVTatvzP33nN1vv/wN3LPptvz/fFHsq6zPQ90mnms\nnSyrJodVqzmyOpIThw7JKSvPyBknvDHHn3ZpqrVGv8ff74QhAACA/hGJmFUeumdt/vGa38gVW27O\nbe12zq0P5MIlJ+W1L/qFnHbOO+b1rVrt5njuveNzufPeK3PXEzfnztF1+V57S+7rtLOtmxxeK3Jc\ndSgnDByUk5efktOPuSSnnv3TaQws6vfo+0QYAgAAmB1EIvquOTmaj3383+YvH/pkvt1q5vzGcN5y\n0Cvz1tf+SZasPKrf480Km9bfme9+5x9z68Nfzu3b7sndrY35fqeZJzvJUbUiJ9UW57ThI3LmIefn\n3DN/NqsPP6vfIz8rYQgAAGD2EYnom/UP3JAPXfET+Ycd67K6WsnPrnxZ/u2b/14Y2gtbnrg3N37r\nr3LTg9fku6Pfzx3trbmn1c2iSnJCrZFT6qvyohWn55wT35RTzvyp1BpDfZtVGAIAAJjdRCJm3P13\nXZkPfvZn8k/jG3LxwJK88+z35KU/9hv9Hmve6LSbufuWf84Nt38sN2+8ObdNrs9dncls6SQn1qs5\no74y5yw/Iy89+W055ayf3q/POhKGAAAA5g6RiBn1gT89Pf918y25dPCA/PYlf53jTnlDv0daMDY8\neFO+9q0/zzcfuS7fGX8ot3Ymsr2bvLBWyxn1A3LOyjNz3mn/Oiee/hP79OwnYQgAAGBuEomYUTd9\n7c+y+qCTc+gx5/d7FJI8ev/1+fq3/jzffOyr+fbEw7m9M5mxbnJSrZYzG6tz3kHnZc1LfimHHH3e\ns36OMAQAADD3iUTA0zx0z9p89ca/zDfWfy03Nh/Lra1WllaSF9VG8uJFx+Xlx74pL335r+TT33xM\nGAIAAJhHRCLgWXXazdx8/V/my3d8JJ/fcFtuGt+WLZUyxePJgd3FecnKU/OBn/jNnHrm6/o9KgAA\nAPtAJAKe1Z5uJXvTK8/Pq47/Tr69YW1umHgot7SbGSiSM2pDefHIMfmx496al13wzjQGFvV7fAAA\nAJ4jkQj4IXv7jKFup5M7b/7/8qXv/m2+tvnmfKu1JY+1y5xar+Wlg4fmgsMuyprz35Vlq47pw78G\nAACA50IkApJM/8OnNzx4U6792h/my49dl280N+TuVidH1ys5t74yr1j10lz00l/NYcetmf5/CAAA\nAM+LSMSM+r133pQjjhvMv/75k/s9CpnZXyXbsW19vnzdh7L2/v+Tb4w/kJtbzSyvJOfUl+S8pafm\nwtN/Nqee/fZUqtVp/V4AAACeG5GIGfWbP3tj/tffnZr/8FO35L98+Jx+j7MgzZafq++0m7nxa/8z\n197xj/na9jtzY3sszTI5uzaY8xefmItP+Zmc9bJfTLXWmLGZAAAAFjKRiBn3D39+e371Px6Zc1/4\n/Vxx8ymp1105sr/NljD0o9z9nY/nqhv/NNdtvinXt7ZltJucUx/M+YtPyEUn/XTOefn/JRoBAADs\nJyIRfXHbTRvy1leUGag384kbRnLMCSv7PdK8M1fC0LP5/u3/J1dd/ydZu+nGfKO1Ndu7ydn1gbxi\n0fG5+KSfyrmv+BXRCAAAYJqIRPTN6PbJvOn0e3LLA0fldz54b37u3af1e6Q5bz6EoWdz3x2fy5XX\n/3G+tPFb+WZra7Z0k7PqA3nFouNy8Yk/lZdc8KuiEQAAwPPUt0hUFMVfJ3ldkg1lWZ42tW95ko8m\nOSLJuiRvK8ty6x7OFYnmkf/8czfmf/zvU/LKs+7IP3zltAwM1vo90pwy38PQs7n/ritz5df/KGs3\n3pDr21uyqdOLRuePHJdLTv5pVxoBAADshX5GopcnGU3y97tFosuTbCzL8kNFUbwryfKyLN+9h3NF\nonnma1c/mHe8oZFq0c2ffbyZCy45st8jzWoLOQw9m3V3fzFXff2Pcu2T1+frrS3ZXiYvrg/lx5ac\nklef+Ys55ayf9utpAAAAz6Cvt5sVRXFEks/sFonuSnJBWZYbiqI4KMnasixP3MN5ItE8ND7Wyr9Z\nc0uu+tYpectFt+QvPnumh1rvRhjae9+79Yp87vo/zpc235yvt0ZTLZKX1Rfnx5afnde9/F056sRX\n9XtEAACAWWO2RaJNZVmu2O39p23vtl8kmseu+PDd+fV/vyL1ajv//R/Hc9Ebj+73SH0jDE2fbqeT\nW771N/nczX+eL227Ize0JrKimpxXX5ELD3xFXrPmP2f14Wf1e0wAAIC+me2RaGNZlj/0s1ci0fw3\nOdHO2y/8bj779VPy42fckb/6wvE5YPVIv8eaEcLQzOi0m/n62g/l83f+Q748dm9ubrZyeK3IKxqr\nc/Fhr8rFF74vS1Ye1e8xAQAAZsxsi0R3Jlmz2+1mXyrL8oV7OK983/vet2t7zZo1WbNmzbTMxOzy\n1S8+kF/+iW4e3roqv/Dzd+WyPz273yPtF8JQ/03s2JQvXfu7uereK3LdxEO5q9XJifVqLhg8LBcd\n/cb8+IXvzeDID13YCAAAMGetXbs2a9eu3bX9/ve/v6+R6Mj0ItGpU9uXJ9lUluXlHlzN7n7/127K\n//iT47J8cHt+60Nb8q9/8eR+j7TPhKHZbdvG+3PVNe/PFx++Kl+Z3JAHO92cVW/kwkUn5rVn/Hxe\n9JL/4CHYAADAvNLPXzf7xyRrkqxMsiHJ+5J8MsnHkhyW5IEkbyvLcssezhWJFqDR7ZP5pbfckiuu\nPjXHH/Bwfu+vq7nw9XPrdiBhaO569P7r89nrfidXPvHVfLW1LWWZvLyxJBetfElev+a9OeTo8/o9\nIgAAwD7p6+1mz5dItLA98uDW/Kc33Zdrbj45Zx95T37/75bmxRcc2u+xnpEwNP90O53ceuPf5TM3\n/c9cs/3O3NiazCHVIhcMHJRXH/G6XHTRBzKy5KB+jwkAALBXRCLmrFu+tT6/9pNP5hv3Hp8zD703\n7/1fg7PmyiJhaGGZHNuaa695fz7//Y/nusmHc0+rmzMa9Vw4fGxee+q/zbmv+BW3pgEAALOeSMSc\nd9tNG/KbP7M+191+Yk4+8MH82oc6ecvPnDjjcwhD7PTEI7fks9e+L1dtuC5fbm7OeJmcVx/JK1ec\nnde//Ddz1Imv6veIAAAAP0QkYt5Yd8/m/Nq/uj/X3HRSDh7emLf+9GP57f/+otTr++8KDmGI5+LO\nb380n/nWH+XqLbfk+tZEVlWLXNBYlVcdcnFefdEHsmTl7LgCDgAAWNhEIuadLZvG876fvy2f+dRR\nGWsP5KKX3533/tmROe6kA6bl84Uh9kW7OZ4vX/O7+ez3PpK14w/kzlYnp9ZruXD4mLzhtH/v1jQA\nAKBvRCLmrU6nm7/6L7flr/9wJHc8eXhOXvVA3vqObfmV3zl9r68uEobYXzY//r187ur/nM89enWu\na27KRJmc31icSw44L2+88Hez+vCz+j0iAACwQIhELAh33vJ4fv+XH8x1Xz0m491Gzjv9e3nn5Svy\n8ouOeMZzhCH64dZv/X0+9a0/zhe335Ebm80cXa/kxwcPy+uP/8mseeV7U2sM9XtEAABgnhKJWFA6\nnW7+3/91R/7uDyu56YFjc9Dg5rzsZQ/llz54cM58ycHCELPK+PYncuVV78ln1306a5uP5/FOmZc2\nhvKq5WflTa/4bQ/ABgAAppVIxIK1betE/ui3bssVHxnKXZuOSlH5btqVj2TJ6Vtz6b94vTDErHPf\nHZ/LJ7/6e7lqy7fz9eZ4VleL/FhjdV571Btz8UW/k6HFq/o9IgAAMIeJRCxIP3jF0LKTXp1Tijdn\n7O4X5a5Nh+eQwU150RkP51/9p8V53b88NtVqpd8jw9M0J0dz3dW/k8987yP50uTDua/VzdmNRi5e\ncnLecPav5NRz/k2/RwQAAOYYkYgF47neSrbpybH82QfvzNWfGMqtDx2ZeqWTU49al9f8ZCfveOeJ\nWbJ0sI//Ctiz9Q/ckE9d+9584cmv5yvN7RkskjWNFXn1wa/May/6vSxbdUy/RwQAAGY5kYh5bV+f\nMdTpdPORv7gz//xXzXz39sOyvrk0xy5Zn1NP35BLf3ZRLv03x7vKiFmn2+nkm1/+w3z61r/ONWP3\n5bZWO6dZoV5jAAAgAElEQVTWa7lw5NhcesYv5qyX/WIq1b37hT8AAGD+E4mYd/bnw6fvvu2J/M1/\nXZfrr1mUux49LK2ymhMOfDgvevGWvPHtK/LKNxwlGjHrbH78e/nsF9+Tzz12Ta5tbk6SrKkvy+sO\nfmVed/EfuMoIAABIIhIxT/TrV8m++Kn78pE/25hbb1qRezcenG6KHLP80Zx8+qZc/LahvPXtJ2Rg\nsLbfvh/2VrfTybe/8We54uY/zdU7vp9bW+2cXq/l4kUn5tKzfymnnfMOVxkBAMACJRIxZ822n6vv\ndLpZ+4V1+eTfbMx3v7k4964/OFs7Qzlq8YYcc/STOfP8dt78M4fllLNWz/hs8Ew2rb8zn/7ib+Vz\n67+Utc2taRTJjzVW5rWHviqvvfj3s3jF4f0eEQAAmCEiEXPKbAtDP8ot31qfj/7Fw/nON+pZ98AB\neWDHgRmuNHP4sg057sTNeelF1bzlHUflkMOX9ntUSLfTyQ1f+ZN86pa/zFVj9+bOVidnNRp51eKT\ncum578zJZ/9Uv0cEAAD2I5GIWW+uhaFn02p1cuUn7ssXP74lt904nAcePTAPTy7PyupYDl72ZI44\naktOe0mZi958YM49/xDPN6Kvnnjklnzq6t/K5x//Sq5rbstIkfz4wKq87ojX5zWv+oMMLV7V7xEB\nAIBpJBIxK82nMPSjjG6fzOc+dl++9oXR3P3dgTz0yMo8MrYy3bKSQ4Y35tAXbMzxp47nnDVDuejS\nw1x1RF902s18fe2H8unb/y5fHL8/97S6OafRyCVLT8ulL3l3TjjjLf0eEQAA2EciEbPGQgpDz8W3\nr380X/jYY7nlG8m6e5fmsc0rsqG1NIsrEzloeHMOWr0lRx4/ntNeUs+Pv+EFOfkMzzpi5mx48KZ8\n8pr35HNPfC1fbo1meaXIhY0D8/qj35yLL/7dDI6s6PeIAADAXhKJ6CthaO9MTrTz1S8+lK98YWPu\nurnIw+sWZf2mFVk/uSzVdHPQ0OasXrElBx++I8ecVOb0l43kglcfmtUvWNzv0ZnHOu1mvnzNB/Pp\nO/8+X5x4MA+0u3lxfSCXLDszb37Fb+fok17T7xEBAIDnQCRixglD06/T6ebm6x/LdZ99PLff2M7D\n9w/miSeW5Ikdy/JEe1EWFc2sGtyaA5ZtzQsOHc1RL+zm9BcP5+WXHJwjjl7e7/GZZx6572u54kvv\nzec3Xp+vNseyulrkosFD8qbj/1V+7KLLUmsM9XtEAABgD0QiZtThb/ilPHTdh4WhGdRqdXLDdY/k\nm1/amLu+085D9w7k8ccX58nty7KhvTi1dHNAfTTLR7Zn5fLtWX3oRA4/rsyJZwzl7FesyomnrvIA\nbZ63dnM811z5n/Opez6SL04+mie7Zc6vL8rrDjw/b7ro97PqkNP6PSIAADBFJGJGXX/7o3nhESuF\noVmi0+nmrlufyLeueyJ3f2ciD36/yPpHBrNp06JsGlucja1F6aSSlbUdWT60LQcs354DXzCegw7v\n5ojjannhixbnzPNWu52N5+x7t16Rj3/1g/n8tltzU7OZk+q1XDxybN5y1q/kjBf/u1Sq1X6PCAAA\nC5ZIBDyrdfdszjevW587btqRdXcn6x8eyObNw9m6YyRbJxdlU3cojXSyvDaWZYOjWbp4R1YcMJ4D\nD2nl8GMrOebkwZx+zgF54RmuSOLptm96MJ/5wm/kM49elWsnN6deJBc2Dsgbjnh9XnPJ5RlavKrf\nIwIAwIIiEgH7ZOfVSLfcsDH33DaWh75f5vFHatn05HC2bB/J1vGRbG6PZDK1LC0msrg+liVDO7Jk\n0XiWLp/IAavbOeiw5NBj6jn25EU55awDcsjhS/v9z2KGdTud3PCVP8knvvvnuWr8vnx/6uHXr1l+\nVt5ywWU58oSL+j0iAADMeyIRMCOe3LAjt3778dx7+2geuGcyjz1Y5sn1tWzdNJCt24azfXw421vD\n2dIdSlJmWWUiixtTQWnxeJaumMyKVe2sXJ284PBaDj1mMMeetDTHn7wyixa7fXG+efT+6/OJa9+T\nz278Rr7aHM9h1UpeOXhoLj3p7Tn/wvekWmv0e0QAAJh3RCJgVul0unn0oW257cYnc99dO/Lg95vZ\n8FCRjY/Xsn3rQLaPDmbHxFBGW0MZ7QxktGxkMO0srk5mUW08I4PjWTQ8kcVLJrNsZSsrDuxm9aGV\nHHxEI0edsCjHnrwshx6x1K1vc8jk2NZcfdVv55P3fixXNzdkWzdZ01iS165ekzdd/AdZcdAL+z0i\nAADMCyIRMKe1Wp3cd/em3HP7ljz4/bE8uq6VJx4rs+nxSrZubmT79sGMjg1mx+RQRttD2dYdSDuV\nLCqaGak0M1ybyPDARIYHJzMyMplFS1tZsqyT5avKrDyoyEGHNXLw4UM58tglOeaFKzI0XO/3P3nB\nu+Pb/5iPX/+hfH77HflOs5VT67VcsvjEvOWcd+a0c9/e7/EAAGDOEomABWfLpvHc973NeeCe7Xn0\ngYk8/kgrGzeU2fxkke1bahndVs+OsYGMTQxmrDmQsfZgdnQa2ZF6BtLJokozw9XJDNcnMjw4keGh\nZoZHmhlZ1MniZZ0sWV5m2aoiK1ZVs/rQgRx02FAOPWJxDj96WQYGa/3+588rW564N5+68jfy2ceu\nybXNrRkpkgsHDswbj7o0r3rV72VwZEW/RwQAgDlDJAJ4jlqtTh68b0vWfW9bHl63I+sfbOaJxzrZ\n/GQyuqWS0e21jO2oZ2yskYnJRsabA5loNzLeGchYWc9EahlIO8NFK0PVZgZrzQzWJzM8OJmhwVaG\nR1oZWdzK4mXdLF2RLF9VycrVtaw4sJ4DXzCU1QcP5+AjFmfp8kG3y+1Bp93M19d+KFfc9te5cvzB\nPNjp5mX1wbx6xTl564//bg495vx+jwgAALOaSAQwQ1qtTh5Zty0P3rctjz00licea+bJ9a1sebLM\n1o1FRrdOhaaxesbHpiJTs5GJTiOT3XomuvWMp5YyyVDaGay0eq9qMwP1ZgbqrQwOtDI41MrQcDvD\nI52MLO5m0dIyS5YXWbqykhWr6lm5upEDXzCUgw4bycGHLZm3t9A98L1r84nr3pfPb7ox32hN5Kha\nJZcMHZ23nPELOeflv5xKtdrvEQEAYFYRiQDmmK2bJ/LIA9uy/tEdeeLRiWx8vJnNT7SzdVM32zcn\no9uKjG2vZmxHNRPjtUxM1DPRbGSyWc9Eu5HJTj0T3UYmylomUks93QwWrQwWrQxU2mlUW2nUWmnU\n22nUWxlotDMw2M7gYCeDw+0MjZQZXlRmZEmyaEmRJcurWXZALctW1rPiwIGsOnA4qw4ezgEHjsya\nK57Gtz+Rz33hXfnkA5/O1c2NSZJXNlbmTUe8Ia+55PIMLV7V5wkBAKD/RCKABazT6ebxx0bz2MOj\n2fDwWDY9OZmtG9vZurmd0S1ltm8ps2M0GR8tMj5WzcRYNZMT1UxO1jPZrKXZrqfZqqfZqaXZqWey\nrGWyW8tkammlkoF0MlC0M1C006i0M1Bt9SJUvRehBhrtDAy0MjDUycBgJ4ND3QwOlxkcSoYXJ8OL\ni4wsqmTxsmoWL61l6cp6li5vZPnKwaxcPZQDDhzZ6+c8dTudfOur/y0f/86f5Qvj9+X+djcvrQ/m\ndSvPzVt+/PdyyNHn7af/2gAAMLuJRADsF5MT7Ty+fjRPrB/Lxg0T2fTEZLY82crWze1s29TN6LYy\nO7YnY6NFJsYqaU70AlSzWU2zVU2zWUurXUurM7Xs1tLsTi3LaiZTTTPVVFNmIJ3Ui04aUzGqXmmn\nXm2nUW2nVmtPXRXVSb3RzsBAJwOD3V6UGi5TjKzLxlUfz/2rrs+9ww/n4M7inDlxQs5ddGmOP+YN\nWbJiIEuXT72WDWbpisHU625VAwBg/hGJAJizOp1utm+dzBPrd2TLxsls2TSZzRsns31LO9u3tDO6\ntZsdo92MbU/GdyQTY70gNTlRZHKymubEzhhVTbNVz2TZzI6jvpJtx16d7Ud+KenUM/D9S1K56w3p\n3rcmrc5QWqmmnm7qU2Fq56tWdFKrdFKvtnvLWju1aif1eif1Wie1eieNRif1RjeNge5UqEoGhsoM\nDCWDw8nQSJHhRZWMLK5mZHE1i5bUMrK4liXLGlm6fCCLlzeyeMlAFi8dEKoAAJh2IhEA7EG308k3\nv/yH+edb/iJXjq/LA+1uzmsM5TXLX5KLznl/hkZOyZbNE9m6aTLbt7ayfWsrO7Z1MrqtnfGxbsZ3\nlBkfLTMxlkyMJ5PjRSYnirSaRZqT1TQnK2m1qmk1q2m1q2m3q2l3qml1ar1lt5Z2t5pWt5pWWUur\nrKSdStqpppVKqummnm5q6aZWdHuRateyF6yqlW5q1fbUspNarZNatdtb1rup1zup1cvUG93U62Xq\nA900GkljsExjIBkYSgYGiwwMJ4NDlQyNVDK8qJrB4UoWLa5neFEti5bUM7KoniXLBrJoScOVVgAA\nc5hIBADPwbq7v5h/Xvu+fHbLt3NDczIn1Ku5ZOS4vPWsX8kZL/53M/praZ1ONztGm9m6eSKj25oZ\n3drK6GgzO7a1M7ajnbHtnYyPdTI+2s3ERDcTY2Umxso0J5LJ8aTZTJpTwao1WUmrVaTVrKTdrqTd\nqvSCVbuSVruaTrcXrzqdStrdatrdqX1lJe2yOvXqBaxWKmmlmsoPBKxqylSnAla16PZelan1ytR6\nZWq92gta1Wo31drO7TK1em+73uimVksvcjWSWqNMo5HUG0l9IGkMJo1GkYHhIgMDlTQGiwwNVzMw\nWMnQSDWDw9UMDdcyvKiWoZFaRhbVM7yonpFFjYwsbghcAMCCJhIBwF4a3fJIPvP5X8+nHv58rm1u\nyUCRXDRwYC495l/k4os/mIHhpf0esW86nW7Gx1q7Atb4WDs7trUyMd7JjtHecnKsm/HxTibHu5kY\n66Y5WaY50YtYzWZv2WomrVbSmizSbiXtVi9mtVtF2q1KOu1K2p0inXalF7CmQlanU0mnu/NV7S3L\nqfWpsNUpK7te7RRpp5JOeqGrkm5qKVNLN9WdkavorVeLbipFObXs7opclZ3LSjeVSplqpdNbVnfu\nL3vBq9JNpVqmVitTrZap1rupTm3X6km1VqZW763Xar3oVav3Xo2BIvVG0QthjUrqjSIDg5XUB4oM\nDFUyMFBNY7CSoeFaGgOVDI/UMzBYzcBQNcOL6hkaqmdwuJah4fqs+dVBAGD2EYkAYB90O5187Uu/\nn0/c9te5cvzBPNzp5hWN4bz+gJfnzRddngMPO6PfI/Ic7QxcO0abvcA12s74WDujo81M7OikOdnN\nxHgnzYluJie6aU7ufJVpNadek0mr1Vv24lYvdnVaU7GrXaTTLtJuJ5120Qtc7SLdqejV7VTS7RZp\ndyrpdnvRq9st0u1W0u729nXLXgTrToWublnstj4VwNLb106RbnrLztT6zhBWTblrfVf0mlqv7FwW\nz7CsdFMUZaq7LStFmUp15/vlrlBWVHpRrLJzWZ0KaLUylUovjlWrZaq1pFpNqvUy1WovlFVrmQpo\nSb1e9LYbSa3Wi2a1+lPLxkAvmtXqRQYGqqk3KmkMVNIY6F1JVm/0olm9UcngUC2NgWqGhuppDFYz\nNFx3FRkARCQCgGl1/11X5mPXvTef3fKd3Nhs5qR6La9edELeeu6v5bRz397v8VjgdoawibF2xsdb\nGRttZXK8k/GxdibG22lOdjM+1k6r2U27VaY50U2z2U27WabZ7KbTTpqT3bTbZdrNpN3uxbFOO2k3\nk06nF8Y67aTdTrqdpN0q0un0oli3U0zFsd6+bqdIp9uLZN1uevGr+1Qs63aLdMteTCvLYiqS7bYs\nn9rulFPHlJV089Ry5/5eJNtt/9R2J0U66V1dVU0vkFV3RbQyRcpeANu5Xey+7KYokkp6kawodju+\n2G3f1HZvWe7aX5mKa0WlfNqyUpla3+1VTEW4aqUX14rKzu2kUu2mUkkq1UxFuExt7xbfapmKcru9\nKsVT6/WkWi1Sq/VCW2//VIirV6aWRSq13na90XtVa0VqjSIDjWrqjWqqtSKNRjXVepFarZKBgdqu\n9Xq9mlq9klr9qXVXtgHMLiIRAOwn2zc9mE9/4dfzqUeuzLXNrVlUSS5qrM6bjv2JXHTJB9MYWNTv\nEYEprVYnkxPtTIy1MznZzthoK+12N63JbprNTlrNbjrtbiabnbSb5dR2b9lq9UJap1Om3Sp7261e\nPOu0e/va7ae2e8teZOt2s2u7005vuzO1vjOmdXvr3U7S6RQpu0U6O/d1i5Td9GLb1Hu9cLbzvadi\n21PvVVKWmVpO7SuLXZGt3LX+VIgryyJlpvZl6rwUu5bdne8/bT1Tme2H17tT68XOEJc8bVlkZ5ib\n2i7KFMmucLcr1qWcCnXlbsc9PebtPHf3iFf8QNRLsfv+7Ip6RZHdQt9u70+Fu12BryhTVDK13gt3\nu46rZCoKTh0zFfp2xb6p9Z3vV6s7I2B2ey+7tqu1p9aLSi/2pVKmWimmPr+3v1JNKkWRylQQ7H12\nMTVfUqkUqdZ629Xq1PrOYypJrVoklaRWq+zav/OcWq2Sokhq9Uoqld53VIre/kq1t6zWe8fXqpUU\n1d7nVCpF6vXqD23X6r3z6vVq7xwBEfpGJAKAGdBpN/OVa38/n7j9b3LVxENZ3ylzfmMkb1i9Jm9+\n1X/JioNe2O8RAWZUp9NNt1um3eqm1eo8fTk5td7uptMq0+500251p6JbN81WJ912eu+3y13LTrsX\n6brtMu1Ob7vbmYpznalX++nbvfjWW3a7SafVuwqu053atyvUZVfU63afepW7gl1vvRf2evGud0wx\nFeV2D3l52rKcOqbTnQpyU9vlVPDrhbtMRbzsCn47/z+uWxZJnjqnnFrPVOxLdt+/+zm77d/DMb3t\npwJfWT59PVNh8IePza7P+aF9P7DsZufn5Gnf1ds7FfimgmGe47Ioev+vWMluy+Kpz9oZF4vdvu1p\n5+x+7O6hMWVSPD1c7h4mk0xFx94sO78jO4/buX+373tqudvnF+XTt3ctf+D4qdme9t7TPmdq1srU\nSLuFz6fOeWq7Ujz9mKIonr6d3vb/396dh0dZHnof/90zw0zCGtlrQEQWWYQimwoqEQSysa8Ba3tO\nT5fTHoUCeqTYgtWjiIqg1vP29fi2PZUQCDtZ2AlUUBCQgoIGQREim+xLmGFm7vePDDS12gokPJnJ\n93Ndc2We+1n4eV33NYYfz3NPSUlq/nodo787x2UUKUhLvVyR7F8Zd5Xe/rtzSh3n/vtxlykpSV3G\nXBm/UpS6Lmcx6pKUoL6Dmv2zjyWIkggAAEd88mGOstdPUe6ZHXo/cEkdvFWUVrOdht//lJq3TXc6\nHgAAjipdIoYiJWHwUlihcMnPcNgqdMn+3XYw9Nf3l8+31upSsGQ8HFTJ/ki5aMO68jMUspFjInf5\nXd6OFI1he7lMLDlfMiXjl48J2yvXkTUKR7atSgrEkrv+rBQ2ClsbKfwUKff+9n04VGr8q8eVel0u\nJyWVKhf/+pJVpIzUlbLx8njJtr0yXlIq6sr7r+4rKRIj1wxH6i+rvztGl/98mch/v/nrcaXOkUpn\nKrXv8vFXysjI8fZvy80r50vSlXO+UpRGxnulHtKMOV3Kbb7GEkoiAAAcdvJooRbkj9fiI2u1LnBe\nN7uNUuKbaEi7n+qepAlyuVlQFwAAAOWPkggAgAok4D+n5flPaOHebK3wH1VQUm9vHQ1qOkipyc8r\nrlptpyMCAAAgRlESAQBQQYVDIW3d+Lqyt72i/OJ9+iwY1r3eeA2oe5+G9H1B9RLbOx0RAAAAMYSS\nCACAKLG/cI3mrp2knFPbtCUQ0B1VPEqr0UbDuz2pVncOczoeAAAAohwlEQAAUejM8U+1MH+CFh9a\nqbWBs6rrMkqJa6TBbX+o+3s9yTpGAAAAuGqURAAARLlgoFgrl0/SwsLZWu4/rGIr9fLdpEGN+6lf\n6ouKr1HP6YgAAACIApREAADEmO3vvqG5701X3vk9+iQYUrcqcepX5x4N6z1VDZt0dToeAAAAKihK\nIgAAYtjBveuVvWaSck68p02X/GrlcSu1eiuNuPsJte38kNPxAAAAUIFQEgEAUEmcO1WkJXnjtaho\nmVYHTquWyyjZd7MGt/qeHujzlNwer9MRAQAA4CBKIgAAKqFQMKC1KyZrwUd/0jL/FzodturlraWB\nicnqn/qSqickOh0RAAAANxglEQAA0M73/lfZm6Yp79xH+igY0l1VfEqv3UXDe01V4m3dnY4HAACA\nG6BClkTGmGRJMyS5JL1prX3+K/spiQAAKCeH929W9sontPT4O9p46aKae9xKrdZCw7uMU4e7f+R0\nPAAAAJSTClcSGWNckgol9ZL0haT3JI201n5U6hhKIgAAboDzZw4rJ+8xLTqYq1X+k6rukpJ9N2tI\nq4dZxwgAACDGVMSS6G5Jk621KZHtJyTZ0ncTURIBAHDjhYIBrV/1tLJ3/V7L/F/oTNjqQW+CBjVO\nV/+06YqvUc/piAAAALgOFbEkGiKpr7X2x5HthyR1tdY+WuoYSiIAABy2Y/MfNGfT88o9v0d7gyF1\n98ZrQN3uGtr3JdVLbO90PAAAAFwlSiIAAHDdDuwp0Ny1k7Tk5BZtCQTUvopH6bXaacT9v1HztulO\nxwMAAMC3cC0lkae8wkQUSbql1HajyNjfmDJlypX3SUlJSkpKKudYAADgmzRukaTxLTZovKRTx/Zq\nYf54LTq8WtMW9lOjpS4lx9+qEXc+os7dH5HL7XY6LgAAACQVFBSooKDguq5R3ncSuSV9rJKFqw9J\n2iwpw1q7u9Qx3EkEAEAU8F84reXLJ2r+3mytCHwpt6S+vvoa0iJDfZKfk8cb73REAAAARFS4x80k\nyRiTLGmmJJekN621U7+yn5IIAIAoEw6F9E7Bi5q787+1rPhzHQ1b9fTW1KCb+2hg2gxVT0h0OiIA\nAEClViFLon8agJIIAICo99H72Zqz8Wnlnt2tXcGg7qri04A6d2l47xfUsElXp+MBAABUOpREAADA\ncYf3b9bclY9r6fFNeufSRbXyuJVWo7VG3vOkWncc4XQ8AACASoGSCAAAVCjnThVpcd44LSxaptWB\nM6rrMkqJa6xh7X6k7g9MZOFrAACAckJJBAAAKqxgoFirl/9K8wrfUr7/iIKS+njraMhtQ5Wc/Lx8\nVWs5HREAACBmUBIBAICoEA6FtHXj68re9opyL+zTwVBYPbzVNLBhLw1Oma6Ees2cjggAABDVKIkA\nAEBU2rcrT3PW/VpLz+zQXwKX1MnrVb+EjhrZ8zk1bpHkdDwAAICoQ0kEAACi3vFDH2jesgladGy9\n3g4Uq5nHpdRqLTWi6wR9964fOh0PAAAgKlASAQCAmFJ89phy8x/T/M+XaJX/pKq7pBRfooa2/r56\n9J4st8frdEQAAIAKiZIIAADErFAwoPWrnlb2rt9rmf8LnQlb9fbdpCG39Fd62nTFVavtdEQAAIAK\ng5IIAABUGjs2/0FZm55X7vlCfRoM635vVQ1s8ICGps5k4WsAAFDpURIBAIBK6bOPV2rO2klacnq7\ntgcuqbPXq/43ddbIB6cp8bbuTscDAAC44SiJAABApXf80AfKzh+vRV/+WRsCxWrpcSu9RmtldPu1\nWt05zOl4AAAANwQlEQAAQCnnzxzW4pyxWlCUr9WBM6rvMkqLb6IRdz6iLveOkcvtdjoiAABAuaAk\nAgAA+AbBQLGW5z+heZ/M1jL/MbklJfsaaFjLh9Sr79PyeOOdjggAAFBmKIkAAAC+hXAopHcKXtSc\nHa8r/+IBnQhb9fLW0pDG/dQ/bbria9RzOiIAAMB1oSQCAAC4Bh9ueUtZ7zyrpecLtTcY0r3eeA2q\n10NDU6ardsPWTscDAAC4apREAAAA1+nAngJlrZmoJae2aVsgoI5er/ondNTIns+pcYskp+MBAAB8\nK5REAAAAZejE4d3Kzv+FFh1br7cDxWrucSu9WktldH9SbTqOcjoeAADAN6IkAgAAKCfFZ49pce5Y\nLTiQq1WB06rjMkqNv0XD2/277kmawDelAQCACoWSCAAA4AYIBoq1cvkkzS/M1DL/EVlJfX31NbRF\nhvokP8c3pQEAAMdREgEAANxg4VBIm/88Q3O2v6bc4v06Frbq5a2pIY3S1D9tuqrVbOh0RAAAUAlR\nEgEAADhs97Y5ytr4jJae2609wZC6e+M1sO59Gpbykup85w6n4wEAgEqCkggAAKACObh3veasnqjF\nJ7doayCgO71V1K9WB43sOVVNWvZ0Oh4AAIhhlEQAAAAV1MmjhZqfP04Lj6zV+sAFNfO4lFqtpTLu\nnqh2XR52Oh4AAIgxlEQAAABRoPjsMS3Nm6AFny/RqsApJbiMUuMaaVi7H6v7AxP5pjQAAHDdKIkA\nAACiTCgY0OplTyr74/9Vvv+IQpL6+upqWLOR6p3ynLy+6k5HBAAAUYiSCAAAIIqFQyFt2fCqsra9\novziz3Q4bNXTW1ODE5M1MH0m35QGAAC+NUoiAACAGPLx9vnK2vCUlp7dpY8j35Q2qN79Gpbysmo3\nbO10PAAAUIFREgEAAMSoon0blLXq8SvflNbJ69WAmzpr5IPTlHhbd6fjAQCACoaSCAAAoBI4cXi3\nsvN/oYXH1uvtQLFaedzqX6OtRt4/RS3bDXI6HgAAqAAoiQAAACqZc6eKtDhvnOYX5WuN/6wSPS6l\nxd+mkV0mqGO3nzgdDwAAOISSCAAAoBIL+M8pP3eCsvdla4X/hKq5pNS4Rhpxx491b89fyuV2Ox0R\nAADcIJREAAAAkCSFggGtW/mU5uz+vfIuHtIlScm+uhrWfJT6pkyVxxvvdEQAAFCOKIkAAADwd8Kh\nkLZseFVZ22Yqt3i/joWtHvTW0pDG/dQ/bbria9RzOiIAAChjlEQAAAD4p3Zty1TWxme15NxH2hsM\n6aiw+vkAAB74SURBVH5vVQ1u0FNDUmcooV4zp+MBAIAyQEkEAACAq7K/cI2y1jyhxae36y+BS+rq\n9WlQnbs0os90Nbilk9PxAADANaIkAgAAwDU7emC75q4Yr0VfbtC7l/xq6/Gof612yujxjG5rk+p0\nPAAAcBUoiQAAAFAmzp74XPNzx2jBoZUqCJzXrR6X0qu2UMY9v1S7Lg87HQ8AAPwTlEQAAAAocxfP\nn1Bu3gRlf7ZQKwOnVNtllBLXWBkdfq677h8vl9vtdEQAAPAVlEQAAAAoV8FAsVYv/5XmFv5Jy/xH\nJUkpvgYafvvD6pX8jNwer8MJAQCAREkEAACAGygcCumdghc1Z8fryrt4QKfCVr19N2norYOUnjpd\nvqq1nI4IAEClRUkEAAAAx+zY/AdlvTtVSy/s0f5gWA94q2vwzX00OPVl1ah9i9PxAACoVCiJAAAA\nUCHs25Wn2eue1JLTO/VhMKh7qsRpQN1uGpH8suoltnc6HgAAMY+SCAAAABXO4f2blbVighad2KT3\nAgF18FZR/1odNLLnVDVp2dPpeAAAxCRKIgAAAFRoJ48Wan7+OC04slZ/DlxQc49b/aq3Uka3X6l1\nxxFOxwMAIGZQEgEAACBqFJ89psW5YzXvQI5WB86ovssoLf5WZXT6hTp1+5lcbrfTEQEAiFqURAAA\nAIhKAf85rcyfqOy9WVrm/1JeSSlx39GI1v+iHr0ny+3xOh0RAICoQkkEAACAqBcOhbR+9TOa+8Eb\nyvMX6UJY6uuro6G3DVVK2ovy+qo7HREAgAqPkggAAAAxZ+uG/9acLdOVc2GfDoXD6umtqaGN0jQw\nfabia9RzOh4AABUSJREAAABi2sfb5ytzwxQtObtbnwRDut9bVUMaPqihaTNUs05Tp+MBAFBhUBIB\nAACg0thfuEaz1/ynFp/arh2Xgrrb69Ogut01Ivll1Uts73Q8AAAcRUkEAACASunI51s1Z8U4LTy+\nSZsDft3praIBCZ2U8eDzatTsfqfjAQBww1ESAQAAoNI7ebRQ8/J+oflH12pDoFi3e9zqX6OtMno8\nrRZ39Hc6HgAANwQlEQAAAFDKuVNFWpQ7VvO/WK41/rNq7HGpX9XmGnXPJLXr8rDT8QAAKDc3vCQy\nxgyVNEVSa0ldrLXbSu2bKOlfJQUljbHWrviGa1ASAQAAoNz5L5xWTt44zftsoVb4T6q2yyg9volG\n3vmIutw7Ri632+mIAACUGSdKotslhSX9TtKEyyWRMaa1pExJXSQ1krRKUouva4MoiQAAAHCjBQPF\nWrl8kuYWvqV8/zF5JaXEfUcj2/xQPXpPoTACAEQ9xx43M8aslTS+VEn0hCRrrX0+sp0vaYq1dtPX\nnEtJBAAAAMeEQyGtX/2M5n7whnIvFskvqa+3roa3GKW+KVPl8cY7HREAgKtWkUqiVyW9Y63NjGz/\nj6Q8a+2CrzmXkggAAAAVQjgU0taNr2v21peVU/yZvgxb9fYmaFiTAUpPm664arWdjggAwLdyLSWR\n61tcdKUxZkep187Iz37XHhUAAACoeFxut7rc94imj92nwolh/bnvn9QsrqF+s+9Pqv9yHQ2YVkN/\nfGuIzp0qcjoqAABlzvPPDrDW9r6G6xZJalxqu1Fk7GtNmTLlyvukpCQlJSVdwx8JAAAAlK22nR/S\ns50f0rOS9u3KU2bBJP324BL9/LUF6l4lToPr99CwlJdVu2Frp6MCACq5goICFRQUXNc1yvJxswnW\n2q2R7TaSZkm6S1KipJVi4WoAAADEiKJ9G5S16nEtOrlF2wIBdfZ6NbB2V2X0eUkNm3R1Oh4AAI58\nu9lASa9KqivplKTt1tqUyL6Jkn4o6ZKkMdbaFd9wDUoiAAAARK3jhz7QnLyxWvjlBr1z6aLaejwa\nUKu9RvWcqltvv5ab8gEAuH6OLVx9PSiJAAAAECvOnvhc83Ie1YLDq1QQOK9mHpf6V2+tjG6/UuuO\nI5yOBwCoRCiJAAAAgAqi+OwxLckdp+wDS7Q6cEYNXEb9qt6mkZ3Hq1P3f3c6HgAgxlESAQAAABVQ\nwH9Oy/Of0NxPsrTMf1zVXVJqXCNltP+ZuiU9Lpfb7XREAECMoSQCAAAAKrhQMKC1KyZr7u4/KM9/\nWFZSiq+Bht/+sHolPyO3x+t0RABADKAkAgAAAKJIOBTSpvUvafb23yr34gGdDlv19dXWsKZDlJo+\nXV5fdacjAgCiFCURAAAAEMW2v/uGMjdNU27xPhUFw+rlq6mhjdI0MH2m4mvUczoeACCKUBIBAAAA\nMeLj7fOVuWGKlpzdrU+CId3nraohDXpqSOoMJdRr5nQ8AEAFR0kEAAAAxKD9hWuUteYJLTr1vnZc\nCupur08D63TTyJQZqpfY3ul4AIAKiJIIAAAAiHFHPt+qOSvGaeHxTdoc8KuDt4oGJNypUQ++oEbN\n7nc6HgCggqAkAgAAACqRk0cLNS/vF1pwdK3eDhSrpcetATXaKqPH02pxR3+n4wEAHERJBAAAAFRS\n588c1qKcMZpXlK81/rNq7HGpX9Xmyrh7otp3/YHT8QAANxglEQAAAAD5L5xWbv4EZX86Xyv8J1Xb\nZZQe30SjOo1Tp24/k8vtdjoiAKCcURIBAAAA+BvBQLFWLJuouXtmaZn/S8UZKdWXqIz2P1X3ByZS\nGAFAjKIkAgAAAPCNwqGQ1q2coqxdbyr34iFZSSm+BhrR6vvq2fdpuT1epyMCAMoIJREAAACAbyUc\nCumdgheVteO3yrl4QOfDUrKvjkY0z1Df1GnyeOOdjggAuA6URAAAAACuydYN/63MLS8q58KnOha2\n6u1N0PBbByk97WX5qtZyOh4A4CpREgEAAAC4bh9ueUuz3vkvLT1fqM+DYfX01dDQxBQNTJ+pajUb\nOh0PAPAtUBIBAAAAKFN7PliizHVPavHZXdoTDOk+b1UNbfighqbNUM06TZ2OBwD4BpREAAAAAMrN\n/sI1mr3mP7Xo1HZ9EAzqnipxGlzvXo1IfUW1G7Z2Oh4AoBRKIgAAAAA3xOH9mzV7xXgtPLFZWwMB\ndfJ6Nah2V43qO0MNbunkdDwAqPQoiQAAAADccMcPfaA5eWO14Mu39U7Ar3ZVPBqUcKdG9Zqmxi2S\nnI4HAJUSJREAAAAAR505/qmycx7V/CNr9OfABbX0uDWgRluNSvovNW+b7nQ8AKg0KIkAAAAAVBjn\nzxzWopwxmleUr9X+s2ricalftZZ6qPuv1KbjKKfjAUBMoyQCAAAAUCFdPH9CObnjNHf/Iq0MnFZ9\nl1F61aYa1XmCOnX/d6fjAUDMoSQCAAAAUOEFA8XKz3tMcz/J0jL/cdVwGaXHNVJGh//QXfePl8vt\ndjoiAEQ9SiIAAAAAUSUUDGjN8l9pzkd/VK7/iNyS0uK+o5FtfqgevadQGAHANaIkAgAAABC1wqGQ\n3l7zrLJ2/k65F4vkl5Tsq6vhLUarT/Jz8njjnY4IAFGDkggAAABATAiHQtqy4VXN3jZTOcX7dTJs\n1cd3k4Y3HarU9Ony+qo7HREAKjRKIgAAAAAx6S+b3lTmpueVc2GvvgiF1dNbU8Map2tA2gzF16jn\ndDwAqHAoiQAAAADEvI/ez1bmhqe05Nxu7Q2GleStpiHf6a2h/V5T9YREp+MBQIVASQQAAACgUvn0\no+XKXDtRi0/v1K5gUPd64zW4Xg8NS5upm+q3dDoeADiGkggAAABApVW0b4Nmr3xMi05t0fuBS+rq\n9WlQnbs1su8M1W/cwel4AHBDURIBAAAAgKSjB7Yra/lYLTz+rjYH/OrgraKBCZ00qveLSrytu9Px\nAKDcURIBAAAAwFecPFqo7NwxWnBsnd4OFKuNx6MBtdpp1APPqWmrvk7HA4ByQUkEAAAAAP/AuVNF\nmp/zqOZ/sVxrA+fVzONS/+qtNfrep3R7hyFOxwOAMkNJBAAAAADfUvHZY1qcO1bZB3K02n9GN7td\nSq96m0bf/YS+e9cPnY4HANeFkggAAAAArkHAf055OeM099N5WuE/qZtcRunxTZTRcYw6d39ELrfb\n6YgAcFUoiQAAAADgOgUDxVq5fJLmFr6lfP8x+SSlxSVqZLuf6N6ev6QwAhAVKIkAAAAAoAyFQyGt\nWzlFWbveVO7FQwpLSvHV14jbv69eyc/I7fE6HREAvhYlEQAAAACUk3AopE3rX9Ls7a8p5+JBnQ1b\nJfvqaHjzkUpJfUEeb7zTEQHgCkoiAAAAALhBtm38nWa/94KWXNino2GrPt4EDb91kNLTXpavai2n\n4wGo5CiJAAAAAMABu7Zl6q0NT2vJ+UIdDIbVy1dTQxulaWD6TMXXqOd0PACVECURAAAAADiscOdC\nzVr/ay0+u1t7gyE94K2mYTcna1D6TFVPSHQ6HoBKgpIIAAAAACqQTz9arsy1E7Xo9E7tDgZ1nzde\nQxv00rD0V1SzTlOn4wGIYZREAAAAAFBBHdy7XpmrHtPCk9u041JQ3bxxGlLvfg1PnaHaDVs7HQ9A\njKEkAgAAAIAocHj/Zs1eMV4LT2zW1kBAXbxeDa5zjzJSXlG9xPZOxwMQAyiJAAAAACDKHCvaoaz8\nsVp4fKM2Bfzq6PVq4E2dldH7Jd3c9G6n4wGIUpREAAAAABDFTh4t1NzcR7Xg6DptuHRRd3g8GpjQ\nQaN7vaDGLZKcjgcgilASAQAAAECMOHP8U83LHav5h1dpfeCCbve4NaDmHRqd9Kxua5PqdDwAFRwl\nEQAAAADEoPNnDmvh0kc0r2iZ1gTO6TaPSwOqt9Goe6fo9g5DnI4HoAKiJAIAAACAGHfx/AktXjpG\ncw8s0Wr/GSW6XRpQrYVGd3tSbTs/5HQ8ABUEJREAAAAAVCL+C6eVmz9Bc/fN04rAKdVzGfWvepsy\nujymjt1+4nQ8AA6iJAIAAACASioYKNbyvMeV9clsLfMfV02XUXr8LRrdcaw6d39ELrfb6YgAbiBK\nIgAAAACAQsGAVi97Ulkf/1F5/qPySUqLS9So7/5c3ZIepzACKoEbXhIZY6ZJ6ifJL2mvpH+x1p6J\n7Jso6V8lBSWNsdau+IZrUBIBAAAAQDkJh0Jat+o3yvrwDeVePCQrKc3XUBlt/009ek+hMAJilBMl\n0YOS1lhrw8aYqZKstXaiMaaNpFmSukhqJGmVpBZf1wZREgEAAADAjREOhfROwYvK3PGaci8e1EUr\npfjqaWTLh/VgyrNye7xORwRQRhx93MwYM1DSEGvt94wxT6ikMHo+si9f0hRr7aavOY+SCAAAAABu\nsHAopK0bX1fm1ulaWrxfp8JWyb7aGtE8QympL8jjjXc6IoDr4HRJtETSbGvtbGPMq5LesdZmRvb9\nj6Q8a+2CrzmPkggAAAAAHLb93TeUuWmallzYq6Nhqz7eBA1rOljpqdPlq1rL6XgArlK5lETGmJWS\nGpQekmQlTbLWLo0cM0lSR2vtkMg2JREAAAAARKld2zL11oanteR8oQ4Gw+rlq6lhjdM1sN+riqtW\n2+l4AL4FR+4kMsb8QNKPJPW01vojY1993GyZpMnf9LjZ5MmTr2wnJSUpKSnpujIBAAAAAMpG4c6F\nmrX+11p8drf2BUN6wFtdwxKTNajfq6pWs6HT8QBEFBQUqKCg4Mr2U089dcMXrk6W9JKk+621x0uN\nX164+i5JiZJWioWrAQAAACCq7duVp8yCSVp85gN9FAzqPm9VDW34oIamzVDNOk2djgegFCe+3WyP\nJK+kywXRu9ban0X2TZT0Q0mXJI2x1q74hmtQEgEAAABAlDmwp0CZqx/XolPva2cwqG5V4jSkfg8N\nT3tFN9Vv6XQ8oNJzdOHqa0VJBAAAAADR7fD+zcpcPk6LTr6nrYGAunp9Glynm0amzFC9xPZOxwMq\nJUoiAAAAAICjjhXtUFb+WC04vlGbA3519Ho16KYuGtV3uho26ep0PKDSoCQCAAAAAFQYJ48Wam7u\no5p/dJ02XrqoOzweDUq4U6N6TVPjFklOxwNiGiURAAAAAKBCOnP8U83LHat5h1fpz4ELut3j1sCa\n7TS651Q1bdXX6XhAzKEkAgAAAABUeOfPHNbCpY8ou2iZ1gbO6TaPWwOqt9boHr9Ry3aDnI4HxARK\nIgAAAABAVLl4/oQWLx2juQeWaLX/jBp5XOpftYVGd3tSbTs/5HQ8IGpREgEAAAAAopb/wmnl5k/Q\n3H3ztCJwSvVdRv2rNtOoux5Xh7t/5HQ8IKpQEgEAAAAAYkIwUKzleY8r65PZyvcfV4LLqF98E43q\nNE6duv1MLrfb6YhAhUZJBAAAAACIOaFgQKvyf6k5hX9Snv+o4oyUFtdIo9r/h+5JmkBhBHwNSiIA\nAAAAQEwLh0Jat3KKsna9qdyLhyRJaXHf0ci2P1KPB39NYQREUBIBAAAAACqNcCikjQXTNHvH68op\nPii/pFRffY28/fvqlfyM3B6v0xEBx1ASAQAAAAAqpXAopC0bXtXsbTO1tHi/Toetkn11NLJ5hvqm\nTpPHG+90ROCGoiQCAAAAAEDSto2/0+z3XtDSC/t0LGzVx3eTRjQdqtT06fL6qjsdDyh3lEQAAAAA\nAHzFh1ve0qyNz2jx+T36IhxWb28tjbh1sNLTXpavai2n4wHlgpIIAAAAAIB/4KP3szXr7clafP5j\nfR4Mq5evpoY1TtfAfq8qrlptp+MBZYaSCAAAAACAb2nPB0uUue5JLTy7S/uCIfX0VtewRqkamD5T\n1Wo2dDoecF0oiQAAAAAAuAb7duUps2CSFp3ZqY+DIfXwVtXQ7/TR0H6vqXpCotPxgKtGSQQAAAAA\nwHXaX7hGmasf16LTf9GuYFD3euM1tEEvDUt/RTXrNHU6HvCtUBIBAAAAAFCGivZtUObKCVpwcot2\nBoPqXiVOQ+onaXj6a0qo18zpeMA3oiQCAAAAAKCcHN6/WZnLx2nhyff0/qWA7qri05C692pk2quq\n3bC10/GAv0FJBAAAAADADXD0wHZlLR+j+cff1ZZAQF28Xg2uc48yUl5RvcT2TscDKIkAAAAAALjR\njh/6QHPyxmrel29rU8Cvjl6vBtfuqlF9Z6jBLZ2cjodKipIIAAAAAAAHnTxaqLm5j2r+0XXaeOmi\n2lepokEJnTS6z8u6uendTsdDJUJJBAAAAABABXHq2F7NyxujeUfWaEOgWG08Hg1K6KDRvV5Q4xZJ\nTsdDjKMkAgAAAACgAjp74nPNzx2jeYdWaF3gglp5PBpYq51GPfCcmrbq63Q8xCBKIgAAAAAAKrjz\nZw5r4dJHlF20TGsD59TM49bAGm00OulZNW+b7nQ8xAhKIgAAAAAAokjx2WNalDNG2QdztNp/Vrd6\nXBpQvbVG3/uUbu8wxOl4iGKURAAAAAAARKmL508oJ3ec5u5fpJWB07rZ5dKA6i31UPdfqU3HUU7H\nQ5ShJAIAAAAAIAYE/OeUk/MLzd03TysCp9TAZdSvajM9dM8kte/6A6fjIQpQEgEAAAAAEGOCgWLl\n5z2mrE9ma7n/hGq7jPpVbarRXR5Xx24/cToeKihKIgAAAAAAYlgwUKwVyyZqTuEs5fm/VE2XUXr8\nLXqo03h16vYzudxupyOigqAkAgAAAACgkggFA1q97EllffxH5fmPKt5I6XGNNarDf+iu+8dTGFVy\nlEQAAAAAAFRC4VBI61ZOUeaHbyjXf0QeSelxicpo/1N1f2AihVElREkEAAAAAEAlFw6F9PaaZzV7\n5//R0otfyEpK8zVUxh0/Vo8Hf01hVElQEgEAAAAAgCvCoZDeKXhRmTteU07xQfklpfnqa2SrH6hn\n36fl9nidjohyQkkEAAAAAAC+VjgU0pYNr+qtrS8r5+IBnQ9bpcTV1YiW31OflKkURjGGkggAAAAA\nAHwrW95+TZlbpmtp8Wc6GbZK9tVWRvNR6ps6TR5vvNPxcJ0oiQAAAAAAwFXb/u4bytw0TUsu7NWx\nsFUf300a0XSoUtOny+ur7nQ8XANKIgAAAAAAcF0+3PKWZm18Rosv7NEXobB6e2tpxK2DlZ72snxV\nazkdD98SJREAAAAAACgzH72frVlvT9bi8x/r82BYvXw1Naxxugb2e1Vx1Wo7HQ//ACURAAAAAAAo\nF4U7F2rW+l9r8dnd2hcMqae3uoY1StXA9JmqVrOh0/HwFZREAAAAAACg3O3blafMgkladGanPg6G\nlOStpqE399WQ9FdUPSHR6XgQJREAAAAAALjB9heuUebqx7Xo9F+0KxjUvd54DW3QS8P7/VY1at/i\ndLxKi5IIAAAAAAA4pmjfBs1aOU4LT27TzmBQ3avEaUj9JA1Pf00J9Zo5Ha9SoSQCAAAAAAAVwuH9\nm5W5fJwWnNys9wOXdLfXp6H17tOI1FdUu2Frp+PFPEoiAAAAAABQ4Rw9sF1Zy8do/vF3tSUQUBev\nV4Pr3KOMlFdUL7G90/FiEiURAAAAAACo0I4f+kCz88ZowZcbtPmSXx2reDWodleN6jtDDW7p5HS8\nmEFJBAAAAAAAosbJo4Wam/uo5h9dp42XLqp9lSoanNBZo/pM181N73Y6XlSjJAIAAAAAAFHp1LG9\nmpc3RvOOrNHbgWK19Xg0KKGDRvd6QY1bJDkdL+pQEgEAAAAAgKh39sTnmpfzqOYdXqn1gQtq5fFo\nYK12eqjXi2rSsqfT8aICJREAAAAAAIgp588c1sKljyi7aJm61Wyt//zpZqcjRQVKIgAAAAAAAFxT\nSeQqrzAAAAAAAACIHpREAAAAAAAAuL6SyBjzG2PMX4wx7xtjlhljGpbaN9EYs8cYs9sY0+f6owIA\nAAAAAKC8XO+dRNOstd+11t4pKVfSZEkyxrSRNFxSa0kpkl43xlzVc3BALCkoKHA6AlDumOeoDJjn\nqAyY56gMmOfA17uukshae67UZjVJ4cj7/pKyrLVBa+1nkvZI6no9fxYQzfifECoD5jkqA+Y5KgPm\nOSoD5jnw9TzXewFjzDOSHpZ0StIDkeFESe+UOqwoMgYAAAAAAIAK6J/eSWSMWWmM2VHqtTPys58k\nWWuftNbeImmWpEfKOzAAAAAAAADKnrHWls2FjGksKdda294Y84Qka619PrJvmaTJ1tpNX3Ne2QQA\nAAAAAADAFdbaq1of+roeNzPGNLfWfhLZHCjpo8j7JZJmGWNeVsljZs0lbf66a1xtYAAAAAAAAJS9\n612TaKoxpqVKFqzeL+mnkmSt3WWMmStpl6RLkn5my+qWJQAAAAAAAJS5MnvcDAAAAAAAANHrny5c\nXZ6MMcnGmI+MMYXGmP90MgtQVowxbxpjjhhjdpQau8kYs8IY87ExZrkxppaTGYHrZYxpZIxZY4z5\nMPKFBo9GxpnriBnGGJ8xZpMx5v3IPJ8cGWeeI6YYY1zGmG3GmCWRbeY4Yo4x5jNjzF8in+mbI2PM\ndcQUY0wtY0y2MWZ35Pf0u652njtWEhljXJJek9RXUltJGcaYVk7lAcrQ71Uyr0t7QtIqa+3tktZI\nmnjDUwFlKyhpnLW2raR7JP088hnOXEfMsNb6JT1grb1TUgdJKcaYrmKeI/aMUckyEZcxxxGLwpKS\nrLV3Wmu7RsaY64g1MyXlWWtbS/quStaNvqp57uSdRF0l7bHW7rfWXpKUJWmAg3mAMmGtfVvSya8M\nD5D0x8j7P6pkoXcgallrD1trt0fen5O0W1IjMdcRY6y1FyJvfSpZy9GKeY4YYoxpJClV0v+UGmaO\nIxYZ/f3ff5nriBnGmJqS7rPW/l6SrLVBa+1pXeU8d7IkSpR0oNT2wcgYEIvqW2uPSCV/uZZU3+E8\nQJkxxtyqkrss3pXUgLmOWBJ5DOd9SYclrbTWvifmOWLLy5IeU0kBehlzHLHISlppjHnPGPNvkTHm\nOmJJU0lfGmN+H3mE+P8aY6rqKue5o2sSAZUYK8YjJhhjqkuaJ2lM5I6ir85t5jqimrU2HHncrJGk\nrsaYtmKeI0YYY9IkHYncGWr+waHMccSC7tbajiq5c+7nxpj7xOc5YotHUkdJv43M9fMqedTsqua5\nkyVRkaRbSm03iowBseiIMaaBJBljGko66nAe4LoZYzwqKYj+ZK1dHBlmriMmWWvPSCqQlCzmOWJH\nd0n9jTH7JM2W1NMY8ydJh5njiDXW2kORn8ckLVLJ8id8niOWHJR0wFq7JbI9XyWl0VXNcydLovck\nNTfGNDHGeCWNlLTEwTxAWTL623+RWyLpB5H335e0+KsnAFHo/0naZa2dWWqMuY6YYYype/kbQIwx\n8ZJ6q2T9LeY5YoK19pfW2lustbep5HfxNdba70laKuY4Yogxpmrk7mcZY6pJ6iNpp/g8RwyJPFJ2\nwBjTMjLUS9KHusp5bqx17o46Y0yySlbfdkl601o71bEwQBkxxmRKSpJUR9IRSZNV8q8V2ZIaS9ov\nabi19pRTGYHrZYzpLmm9Sn7BspHXLyVtljRXzHXEAGNMO5Us8OiKvOZYa//LGFNbzHPEGGNMD0nj\nrbX9meOINcaYppIWquT3FY+kWdbaqcx1xBpjzHdV8kUEVSTtk/Qvkty6innuaEkEAAAAAACAioGF\nqwEAAAAAAEBJBAAAAAAAAEoiAAAAAAAAiJIIAAAAAAAAoiQCAAAAAACAKIkAAAAAAAAgSiIAAAAA\nAACIkggAAAAAAACS/j/9o0oPjsSv4gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x16164dd90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(20,10))\n",
    "for ind,(omegas,converged_omegas) in converged_dict.iteritems():\n",
    "    #plt.figure(ind)\n",
    "    plt.plot(omegas[::pow(10,ind-3)],converged_omegas[::pow(10,ind-3)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# omegas_converged = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# for i in range(3,6):\n",
    "#     omegas,converged_omegas = converged_dict[i]\n",
    "#     omegas_converged[i] =  [tup[1] for tup in converged_omegas]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolations "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: Using UnivariateSpline with k=2 and a small value of `s` seems to offer a good balance of speed and precision."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@timeit\n",
    "def interpolate_omegas(omegas_few,omegas_many,omegas_few_values):\n",
    "    '''\n",
    "    A method for interpolation from omegas_few to omegas_many.\n",
    "    \n",
    "    Args:\n",
    "        omegas_few (list of floats):\n",
    "            points to use to interpolate\n",
    "        omegas_many (list of floats):\n",
    "            points to interpolate to. Should have omegas_few included.\n",
    "        omegas_few_values (list of floats):\n",
    "            The values of the omegas_few that have already been computed\n",
    "            \n",
    "    Returns:\n",
    "        A list of length len(omegas_many) interpolating the values of omegas_few.\n",
    "    '''\n",
    "    f = interpolate.UnivariateSpline(omegas_few, omegas_few_values,k=2,s=1e-10)\n",
    "    return f(omegas_many)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "'interpolate_omegas' 0.03 sec\n"
     ]
    }
   ],
   "source": [
    "omega_values_interpolated = interpolate_omegas(np.asarray(omegas3),np.asarray(omegas5[1:99900]),np.asarray(omegas_converged[3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x128c23410>"
      ]
     },
     "execution_count": 159,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEKCAYAAAACS67iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8lfXd//HXJ4skJCQhQNhhyRYRFbfEhYM66qrS1lm7\nbLV12wXeLda7vXs7av1ZF3UhVq0K2ioqRNyCIEOWshIIJMyQEAIZ398f58L7GBPIOMmVc5338/HI\ng5zrOtd1fb7XObzP93yvEXPOISIi0S/O7wJERCQyFOgiIgGhQBcRCQgFuohIQCjQRUQCQoEuIhIQ\nCvRWZmZLzewkv+toDDNba2antNG2pprZf7XFttobM5toZq9H+rmNWFervL5mVmtmAyK9Xmk6BXor\nc86NdM7Nbcxz2zJQW0r/iZvPOTfNOXdmpJ/rI13M0k4o0APEzKwNN9fu/xObWbzfNdTVHmuKgLZ8\n38kBKNBbWXiv28wmmdlzZvaEme0ysyVmNsab9yTQF5jpzbvZm36Mmb1vZjvMbKGZjQtb9xwz+4OZ\nvWdmu4H+3rS7zOxjMys1s5fMLDNsmXO9YaDtZjbbzIY2UPdRZvaBt92NZvZXM0vw5r1D6D/xYq/W\ni73p3/Jq3OHVdGjY+g43s0+9mqYDyQ1sN8lbfnjYtC5mVmFmXRqxnbVmdquZLQLKzSzOzG4zsw1e\nrcvN7GTvuV8b9jGzcWZWGPa43uXqqbmTmT1pZiXe9n8dNu8Kr8b/NbOtwCRv2rthzxlvZiu89vzN\nzPLN7Oqw5cOfW2tmPzKzVd5r+EDYvAFm9raZbfVqedrMOtVXc536x5rZpvAOgZl929uHB3wv1LOu\nOftrb6D+oWY2y8y2efv04rB5Z5vZ597+LjSzGw9Wu9ThnNNPK/4Aa4FTvN8nARXAGYQC8S7gwzrP\nPTnscU9gK3CG9/hU73G293gOsA4YSujDOcGbVggMA1KAF4CnvOcPBsqBU4B44BbgCyChnlrHAGO9\nOvsCnwPXh9VWC/QPe3w4UAwc6S3zfW99id7POuB6b7sXAvuA/2pgnz0K/D7s8U+Bfx9sO2FtWODt\nuw5emwuAHG9+3/11A1PDawDGAQVh+6re5eqp90ngJSAVyAVWAld5864Aqrw2xHk1XQHM9eZ3AUqB\n87z51wN7gavDlp9bZ7/PANKBPkAJMN6bN9B7jyQA2UA+8L/1vRfracMXwKlhj/8J3NKE98KAsPfk\n1WHzwtua6u3Ty711HQZsAYZ684uA47zfM4DRfv//jbYf9dDb3nvOuTdc6F37FDCqzvzwr6/fA15z\nzr0B4Jx7G5gPnB32nH8451Y452qdc9XetKecc8udc3uA3wIXe72vS4BXnXOznXM1wP8QCv3j6hbp\nnFvgnPvEhRQADxMKvIZqvRZ4yDk331vmKULBdIz3k+Ccu985V+OcexGYd4B99CxwWdjjicAzjdjO\nfvc554qcc3uBGiAJGGlmCc65Aufc2gNse79GLWdmccB3gNudcxXOufXAXwh90Oy30Tn3oPca7a2z\nirOApc65V7z59xP6wDqQPzrnypxzhYQCdDSAc261c+5t51y1c24bcA/ffM0aMp3QfsbM0gm9x6Z7\n623Me6ExvgWsdc496a1rEfAisL+Xvg8YYWbpzrlS59xnzdhGTFOgt73NYb9XAMleKNQnF7jE+2q9\n3cx2AMcD3cOeU1jPcuHT1hPqIXch1Gtdv3+G96FSCPSquwIzO8TMZnpfxXcCU7x1NCQXuKlOrb29\nbfYENtZ5/vq6KwgzB0jxvurnEurJvdyI7ey3IayNq4FfAJOBYjObZmbh+69eDSzXo56ndiHUIy6o\n07bwfVrfa7Rfz3rmb6jviWHCA78CSAMws25m9qw3TLQTeJoDv2bhpgHfNrNE4ALgU+8DoznvhYbk\nAsfUee0mAjne/AuBCcB6b+jmmIZWJPVToLcvdQ80FgJPOuc6ez9Zzrl059yfD7AMhL6K75dL6Cv/\nVkJfaXPreW59AfL/gOXAQOdcJvBrDnzwqxCYUqfWNOfcc8Amvvmh0behFTnnagl95Z9IqKf+qnNu\ndyO289Uq6qxvunPuxLC2/7f3725CwwD79TjIcnfXU+5WQvs3fL/m8vUPsAMdQN7E118vCH1ANcdd\nhIY/Rniv2fdo5AFL59xyQh9EZxPa59PCZjflvVB3n9btfOTXee06Oed+5tXwqXPufKAr8Aqh94A0\ngQLdf+H/MTYD4acCPg2c4x00izOzZO/AXU8O7HvewadU4E7gea83/k9ggpmdbGYJFjrwWgl8WM86\n0oFdzrkKCx04/Umd+XVrfQT4sZmNBTCzjt5Bro7e+qvN7Ofedi8gNCZ7IM8SGsqYyNfD5UDb+QYz\nG+y1N4nQV/o9hEIP4DPgbDPL8nrtNzRyua+EffhMMbM07xvFLwkNpzXGa4SGdc41s3gz+xn/12Nt\nqnRCx0jKzKwXoWMkTTGN0D44EXi+znoP9F4I9xlwgZmlmNkg4Jqwea8Cg83se977INHMjvTeq4kW\nOue+kzccWEZo2EuaQIHe+g52el/4/LuB33pfR290zm0gdLDsV4QOHq0Hbub/XreG1v0U8AShHnkS\nXlA551YR6rU94K1vAnBO2Nh7+PpuBr5rZruAv+ONp4aZDDzp1XqRc+5TQuPbD5jZdmAVoQNiOOeq\nCH2NvwrYRmjM9MUD7hTnPiHU2+sB/CdseoPbaWCfdCC0X7d4+6MrcEfYflpM6IDt63XaeKDl6rqe\n0NDHGmAu8LRzbuqB2hfWnv3748+EevtDCR0nqTvW3lD7wt0JHAHsBGbyzX18sPfidOAk4G3n3Paw\n6Qd7L4Sv9x5C31g2Ezro/PRXT3KuHBgPXEponxYR2sdJ3lO+D6z1hnV+iDemL41noY7bAZ5g9hih\ngxnFzrlR3rQ/AecQetOtJnREf1cr1yqNYGZzCB0UfdzvWqTpvIPXG4CJzrl3/K5HoktjeuhTCZ1m\nF24WoXG60YROd2qo5yIiB+ENqWWYWQdC49MAH/lZk0Sngwa6c+49YEedaW95Y4cQeuM19yCORF67\nv4JTvuFYQt90SwgNg51Xz+mNIgd10CEXAO9Az8z9Qy515s0Apjvnpn1zSRERaSstOihqoUucqxTm\nIiL+q/d+DI1hZlcSOmf1gHcHNDMNAYiININzrkk3PmtsD90IO1/azM4kdI7ruY0Z63Pt4B4HrfUz\nadIk32tQ+9Q2tS94P81x0EA3s2nAB4QuCCgws6uAvxK63PhNM1tgZg82a+siIhIxBx1ycc7Vd3J/\noy6aEBGRtqMrRVsoLy/P7xJaVZDbF+S2gdoXixp12mKLNmDmWnsbIiJBY2a4VjooKiIi7ZwCXUQk\nIBToIiIBoUAXEQkIBbqISEAo0EVEAiLwgb5u8w4ybzzB7zJERFpd4AO9Z3YnStM+Yddu3V5aRIIt\n8IGelBhP/J4cFq/d5HcpIiKtKvCBDpBc3YulBRv9LkNEpFXFRKBnWC9WFinQRSTYYiLQuyT1Yu02\nBbqIBFtMBHrP9F5s3KVAF5Fgi4lAz83qRUmlAl1Egi0mAv2Q7r0orVWgi0iwxUSgj+jTi4r4Ir/L\nEBFpVTER6KMH9qQqeSO1tfpDGyISXDER6N07p0FtEuuKd/hdiohIq4mJQAdIquzForUaRxeR4IqZ\nQE+r7cWyQgW6iARXzAR6VkIvvixRoItIcMVMoOek9GL9DgW6iARXzAR636zeFJVt8LsMEZFWEzOB\nPrRHX7ZWFfhdhohIq4mZQB/dL5ddcev8LkNEpNXETKAfPTSXqtQCXVwkIoEVM4HePSsdq05hWcEW\nv0sREWkVMRPoAMl7+zFv1Tq/yxARaRUxFeiZ5LKkcL3fZYiItIqYCvTuyf1YVbzO7zJERFpFTAV6\nbmYuBbvUQxeRYIqpQB+a04+SvQp0EQmmmAr0Ubm5lNo6v8sQEWkVMRXoRw/pR2XyepzTuegiEjwH\nDXQze8zMis1scdi0LDObZWYrzewNM8to3TIjo3+PTMD0hy5EJJAa00OfCpxRZ9rtwFvOuSHAbOCO\nSBfWGswguTKXj1ZoHF1Eguegge6cew+o26U9D3jC+/0J4PwI19VqMunP/NVr/C5DRCTimjuG3s05\nVwzgnNsMdItcSa2rd+ogPt/0pd9liIhEXEKE1nPAo4yTJ0/+6ve8vDzy8vIitNmmG5x9CPOLPvVt\n+yIi9cnPzyc/P79F67DGnPFhZrnATOfcKO/xciDPOVdsZt2BOc65YQ0s69rTWSV/fvEtprw7hZ33\nzvG7FBGRBpkZzjlryjKNHXIx72e/GcCV3u9XAK80ZaN+On7YIZQlfuF3GSIiEdeY0xanAR8Ag82s\nwMyuAu4GTjezlcCp3uOocNTg3tR22MbW0gq/SxERiaiDjqE75yY2MOu0CNfSJhIT4kmq6M/cpau5\n4PhD/S5HRCRiYupK0f2y3CA++VJnuohIsMRkoPdMHsTSjQp0EQmWmAz0wdmHsKZUB0ZFJFhiMtAP\n6zOITfvUQxeRYInJQD9xxCHsSlAPXUSCJSYD/eihfajtsI2NW8r9LkVEJGJiMtATE+JJqRjCmwtX\n+F2KiEjExGSgA+TYcN5ftczvMkREIiZmA31Q5nCWbFagi0hwxGygj+kzjHW7FegiEhwxG+jjhg9n\nW5wCXUSCI2YDPW/UQKpTNrKjbI/fpYiIRETMBnpqciIdKgby5sKVfpciIhIRMRvoAF3ccN5budzv\nMkREIiKmA31gxnAWFGocXUSCIaYD/ajcEXy5a4nfZYiIRERMB/pZo0ezNWGR32WIiERETAf6uFED\nqemwlbWbdvpdiohIi8V0oCfEx5FeMYpXPlIvXUSiX0wHOkBu0mjeWfWZ32WIiLRYzAf6YT0OY+kW\nBbqIRL+YD/RTho1mY40CXUSiX8wH+jlHj2RPx5XsrtzndykiIi0S84HeNTOVpIp+/Gee/tiFiES3\nmA90gF5xY3ht4Xy/yxARaREFOjAmZywfb5jndxkiIi2iQAfOGjWWdVUf+12GiEiLKNCBi04YzZ6O\nK9hWqnuji0j0UqADGR2TSd09nOffW+h3KSIizaZA9wzocDSvL9Gwi4hELwW655g+Y1m45RO/yxAR\naTYFuue8I8dSZOqhi0j0UqB7zjxyCDWJO1m0erPfpYiINIsC3ZMQH0e3yhP4x5x3/S5FRKRZFOhh\njux6Im9/OdfvMkREmqVFgW5mvzSzpWa22MyeMbOkSBXmh28fcRJf7FWgi0h0anagm1lP4OfAGOfc\nKCABuDRShfnhsnFjqExdw9pNO/wuRUSkyVo65BIPdDSzBCAVKGp5Sf5JTU4ka/fRTH37fb9LERFp\nsmYHunOuCPgLUABsBHY6596KVGF+OSzrJN5YoWEXEYk+LRlyyQTOA3KBnkCamU2MVGF+OXfUSSzb\nrUAXkeiT0IJlTwPWOOe2A5jZv4DjgGl1nzh58uSvfs/LyyMvL68Fm21dV552NDfO/5z1xaXk5mT4\nXY6IxIj8/Hzy8/NbtA5zzjVvQbOxwGPAUcBeYCowzzn3tzrPc83dhl+6/HI8Pzjsp9x95fl+lyIi\nMcrMcM5ZU5ZpyRj6J8ALwEJgEWDAw81dX3tybLfxzFw+y+8yRESapNk99EZvIAp76C+8u5jLXr6A\nqr986XcpIhKj2rSHHmQXHH8otfG7mbNotd+liIg0mgK9HnFxRr+a8TwyW8MuIhI9FOgNGD9wPPmF\nCnQRiR4K9Ab8fMLpbOowh9LyvX6XIiLSKAr0BgzP7UZ65QjunznH71JERBpFgX4AJ3Y9n+kLX/G7\nDBGRRlGgH8D1489jBa9QXVPrdykiIgelQD+AM44cTEJ1Jk/Pnu93KSIiB6VAP4gxqefxyLsv+12G\niMhBKdAP4gcnns/88peorY2uq11FJPYo0A/iqtPHUhO/h+ffXex3KSIiB6RAP4i4OOOo5Eu5581n\n/S5FROSAFOiNcOP4y5hfOV3DLiLSrinQG+HC40cRX5vKI69/6HcpIiINUqA3QlyccULmZTw4d7rf\npYiINEiB3ki/PnciS9x0yvfs87sUEZF6KdAb6ZTRA8nYO4JJ02b4XYqISL0U6E1w6ZBreHLpo36X\nISJSLwV6E0z57oVs6zCPD5cV+F2KiMg3KNCboHOnFEa4y/jV81P9LkVE5BsU6E10x5nX8G754+yr\nqvG7FBGRr1GgN9HEkw8npboHv3tmpt+liIh8jQK9Ga4cdgMPfXaf32WIiHyNAr0Z/vvyiyhL/ILn\n3vnM71JERL6iQG+G1ORETut0Hb999X6/SxER+YoCvZn+etUP+TLhJZauLfa7FBERQIHebIN7ZzO8\n9jKuffxev0sREQEU6C3y0Pdv4+Oqh1ldtN3vUkREFOgtccLIXAbVnM81j+iMFxHxnwK9hf526R3M\n3fM3CreU+l2KiMQ4BXoLnX7EIHL3ncXVDz3gdykiEuMU6BHw/yb+hrd336uxdBHxlTnXun8n08xc\na2+jPRh+64/omNiJeVP+7HcpIhIAZoZzzpqyjHroETLtR5P5tOZxPli23u9SRCRGKdAjZPTAHhyf\ndB2XT/2t36WISIxSoEfQsz+/hbVxs3hm9gK/SxGRGNSiQDezDDN73syWm9nnZnZ0pAqLRr27pnN5\n7z/w4xnXUV1T63c5IhJjWtpDvw/4t3NuGHAYsLzlJUW3R356Nc4cP3zwH36XIiIxptmBbmadgBOd\nc1MBnHPVzrldEassSiXEx/H3cx7kiQ13sGaTTmMUkbbTkh56f2CrmU01swVm9rCZpUSqsGj23VPG\nMMxdxIR77vC7FBGJIS0J9ARgDPA359wYoAK4PSJVBcBrN03hC17jL/+a43cpIhIjElqw7Aag0Dk3\n33v8AnBbfU+cPHnyV7/n5eWRl5fXgs1Gh9ycTH5z2N+5/YNr+G7eYrp3TvO7JBFpx/Lz88nPz2/R\nOlp0paiZvQNc65xbZWaTgFTn3G11nhMTV4o2ZOBNV9IxKY3Ff9S9XkSk8fy4UvR64Bkz+4zQWS53\ntXB9gfPWzffyefUrGnoRkVane7m0gSnPvc6k+T9k2Q0LGdw72+9yRCQKNKeHrkBvI0f95mYKdq9k\n019mEBfXpNdIRGKQbs7Vjs359V3sdiVc9Bf9dSMRaR0K9DaSlpLEjCue5eWtdzF11jy/yxGRAFKg\nt6FTDh/AbSP+zrVvXsjiNZv9LkdEAkZj6D7ImzyZBaVvsnHKbNJTO/hdjoi0QzooGiWqa2rpc9OF\nZHbowud/fFgHSUXkG3RQNEokxMcx79dPsrbqQy7801/9LkdEAkKB7pPeXdN5+5pXmbn9T9z02At+\nlyMiAaBA99HxI/rxzIRXuWfVT3lg5ly/yxGRKNeSm3NJBHxn3GjWlzzLDe9dTI/Mt7jwxEP9LklE\nopR66O3ArRefys8G3cclM8/k35+s8LscEYlSCvR24r5rL+Wa3Ls494XTeGvBF36XIyJRSIHejjx8\n3RVM7DmZM589lblL1vhdjohEGQV6O/PkL37ARTl3cMqTJ/P6vFV+lyMiUUSB3g5Nv/knXN73d0x4\nPo/p+Qv9LkdEooTOcmmnHv/5NWQ/nsXE/5zB9t0v8tMJJ/pdkoi0c+qht2N/vvoCphwxjZ/NvZA7\n/jHD73JEpJ3TvVyiwJNvz+OqN87nW51v5qVbf6F7v4jEAN2cK8A+WLaeUx89h4GJxzPvv+4npUOi\n3yWJSCtSoAfchi27OHzKpdRSzce3TGdQr85+lyQirUR3Wwy43l07UfinGQxIG8mwe47kmdkL/C5J\nRNoR9dCj1E2PPc89X/yUK3rezdTrr/G7HBGJMA25xJhXP17Ohc9dSG78sbz7q/vJyerod0kiEiEa\ncokx3zp6GOt+/Qm1VNHn90fw9Nuf+l2SiPhIPfSAuP7h6Tyw+nrOzLiRV269hcSEeL9LEpEW0JBL\njPvg8wLOevhyLA5evnIqeYf197skEWkmDbnEuONG9GXL/7zNCd0mcMqzR3H+3feyr6rG77JEpI2o\nhx5Qsz5dxSVPXUuN7eWZ7zzGuceM8LskEWkCDbnI11TX1PL9ex/huS2/4bikn/DSjbfTNTPV77JE\npBE05CJfkxAfx7M3/YiPr15I4e6V9PjDcG6Z+iK1tfqAFQki9dBjyP++NJtfvXs9HV13nrzsfiaM\nHe53SSLSAA25yEFVVFYx8d4HmbHzD4x0lzHtJ79hZL9ufpclInVoyEUOKjU5kZdvv4HlP1tGfLwx\n6qFhnHznnWzaXuZ3aSLSQuqhx7h3Fq/hiid+S2H8bC7K+Q2P/uRa0lOT/C5LJOZpyEWa7ZnZC/nF\nzDvYkbCci7rfzsM/vppOHTv4XZZIzFKgS4s99saH3PH6H9ieuIhvd7uVh390LVnpKX6XJRJzfAl0\nM4sD5gMbnHPn1jNfgR6FnnhzPre99ge2JH3CmRm/5MEf/JDcnAy/yxKJGX4F+i+BI4BOCvTgee6d\nRdw+879Zn/g6o7mS+797PSeM7Od3WSKB1+ZnuZhZb+Bs4NGWrEfar++MO4y1/zONj65YREJcAic9\nfQR9b7qEx974GH1Oi7QvLeqhm9nzwBQgA7hJPfTg27i1jOsefZzXtt5LUnVXLu7/Y/58xaW6pYBI\nhLXpkIuZTQDOcs79zMzyCAX6OfU8z02aNOmrx3l5eeTl5TVrm9J+VFXX8McX3uDBjx+ipMP7jLLv\n8vvzf8Q5R+smYCLNkZ+fT35+/leP77zzzjYN9LuA7wHVQAqQDvzLOXd5neephx5wH3xewK3PPcqH\nlY+SVjWIiwf8gN9PvIAe2Wl+lyYStXw7bdHMxqEhl5i3e08Vk6fP4Kkl/6C4w7sMqDqPnxx7Ob84\n/2QS4nVRskhTKNCl3Vi6rpjfPvcsr29+gqqEbRyd+j1uO/P7nHvsML9LE4kKurBI2qV/vrOYP73x\nFAurp5FYlc1xmRdzy9mXcNZRQ/wuTaTdUqBLu1ZdU8vD/3mfh99/niU1L5BU3YUTsi7h5rMv5owj\nFO4i4RToEjWqqmt56LX3eeSD5/ncvUBiVTZHpJ3HVcedw5Xjj9KYu8Q8BbpEpeqaWh5/80P+8cFM\nFpTPpCpxG4PcBC469Bx+ee7pdMno6HeJIm1OgS6BMOezNdz/+kzyN89gZ+o8ulacyAk9zuDqceM5\ne+wQ4uKa9B4XiUoKdAmcdZt38td/z+K1FbP40s0CYCCnc9aQ8Vx31qkc0quLzxWKtA4FugRaba3j\njfmreCx/Fu8WzaIkdS6pew7h0I6nc/aIcVxz+vH06pLud5kiEaFAl5hSVrGPx2d9xL8WvcWine9Q\nmvopHSuGMyz1JM4aPo5rTjuR3JxMv8sUaRYFusS0neWVPDX7E17+7B0Wbn+HHR0/JmXPIAZ3OIlx\n/Y/nkuOO4djhfTQGL1FBgS4SpqxiH8/M+ZQZi95h4dYPKUn8CFw83auOZXTXYzj70GO4dNwRZHfS\nnSKl/VGgixxAba3jvaXr+OeHH/Hu2o/4svIjKjouJXXPUPonHsPYXmM5a/QYvnX0MFI6JPhdrsQ4\nBbpIE23fVclzcxfyxucfsmjLfIrcAvYlF9Jx90hyk8ZwRM8xnDFqDOceM5L0VP3RbGk7CnSRCNi0\nvYyXPljEW8sWsKh4ARtrF7C345ekVAyhV9wYRnQZxXEDR3L2kYcyol83TEPy0goU6CKtZOvOPbzy\n0RLmrFjA4uIlFFQuYVfyEqhNJHPvoeSmjGRU95GcNPRQJowdQfcsnT4pLaNAF2lDtbWOhauLeH3B\nUj5as5RlW5ewqWYpezouJ35fVzpXjaRPyhCGdRvC2AFDOWXUEIbndtNZNtIoCnSRdmDvvhreXbqG\nOcs+57OClXyxcwWbq1ZS3mElLq6GjnuG0D1hCIMyh3JY7yEcN3gIJ40cRFZ6st+lSzuiQBdp51YW\nbuXtRSv4ePVKlpWspGD3SnbEr6AqdT3xe3JIqx5ITuIA+nUayNDuAzii/0BOGDGA/t2zMA3WxxQF\nukiUqtxXxScrC/ho1RoWF6zhi22r2Vixhu1uNXtSVmPEkVI5kGwbQK/UgQzsPIBhPfozun9fxg7p\nQ9dMnUsfNAp0kQCqrXWs2rCd95etZsHaNazcsobC8tVsrVpHWXwBVSmFWFU6yXv7kklfcpL70jej\nL4O69eXQPn05cnAuw/p0Iz5O95iPJgp0kRhUU1vL8oItfPpFAUsLC1hVUkBBaQHFlQWUugL2dCjA\nJZWSUNGb1JreZMb1omtyD3qm9yS3c08G9+jJ8L49GT2gp66abUcU6CJSr22le1i4egNLCwv5YvMm\n1m8roqi8iK2VRZTWFlGRUER1chFW24GkvT3pWNuTrISedEvpSY+0nvTJ6kG/LjkM6tGNoX260b97\nFvHxGtNvTQp0EWm22lrH2k07WbS2iBUbi/iyuIjCnUUUVxSxvaqIstoSKuNKqEoqwSXuJq6yKx2q\nupFKNzrFdyO7Qw7dOnajV2Y3+mZ3Y0BODgNzunFIr6507qQzeJpKgS4ibaKsYi8rCrewamMJa4pL\nWL+tmKLSEkp2l7B9bwm7qkvYHVfMvoQSalNKoLoD8fs6k1STTbLLJi0um06JnclKzqZLajY5nbLp\nmZlNr86d6ZeTzcAe2fTpmhnT4/4KdBFpd2prHcU7y1i9aRvrS7ZTsHUbRTu2UbxrG1t3b2f7nm2U\nVm2jvGYbe9jO3vhtVCdtg6QybF8GCfuySarpTDKZpMRlkBafSXpSBhnJmWSlZNAlLZOu6RnkZGbS\nPTOD3tmZ9O6aQU5WWlR/ICjQRSQwKiqrWV+yg4It29mwbTvFO0sp2bWTreWlbK/YSWllKaV7d1Je\nvZOK2lIq3U72xZVSFb+T2qRSSNiD7UsnvjqDhOpMklwGHehEclwaKfHpdExIo2NSGp2S0umUnEZG\nSjpZHdPonJZOdnoaXTqlkZOZTk5WGt2z0klOSmzT9ivQRUQ8FZXVbNiyiw1bd1JcWkpx6U62l5ex\nY3c5OyrKKKssZ9feMsr3lbO7qpyKmjIqa8vZ68rYRznVceVUx5dRm1AOSWVQG49VpxFXnU5CTRoJ\nriNjsyZ+9OivAAAEgklEQVQw+3e/a5X6mxPouumziARSanICg/t0ZnCfzi1eV02NY0fZXkp2lrOl\ntIytZWVsL99NTmZGBCqNHPXQRUTaoeb00KP3iIGIiHyNAl1EJCAU6CIiAaFAFxEJCAW6iEhAKNBF\nRAJCgS4iEhAKdBGRgGh2oJtZbzObbWafm9kSM7s+koWJiEjTtKSHXg3c6JwbARwLXGdmQyNTVvTI\nz8/3u4RWFeT2BbltoPbFomYHunNus3PuM+/3cmA50CtShUWLoL+pgty+ILcN1L5YFJExdDPrB4wG\nPo7E+kREpOlaHOhmlga8ANzg9dRFRMQHLbrbopklAK8C/3HO3dfAc3SrRRGRZmjTP3BhZk8CW51z\nNzZ7JSIiEhHNDnQzOx6YCywBnPfzK+fc65ErT0REGqvV/8CFiIi0jVa7UtTMzjSzFWa2ysxua63t\ntCUze8zMis1scdi0LDObZWYrzewNM2tff5OqkRq6UCxA7etgZh+b2UKvfZO86YFoH4CZxZnZAjOb\n4T0OUtvWmdki7/X7xJsWpPZlmNnzZrbc+z94dHPa1yqBbmZxwAPAGcAI4LKAXHQ0lVCbwt0OvOWc\nGwLMBu5o86oio6ELxQLRPufcXuBk59zhhE6xPcvMxhKQ9nluAJaFPQ5S22qBPOfc4c65sd60ILXv\nPuDfzrlhwGHACprTPudcxH+AYwid+bL/8e3Aba2xrbb+AXKBxWGPVwA53u/dgRV+1xihdr4MnBbE\n9gGpwHzgqKC0D+gNvAnkATO8aYFom1f/WiC7zrRAtA/oBKyuZ3qT29daQy69gMKwxxsI7lWk3Zxz\nxRC6ehbo5nM9LRZ2odhHhN5QgWifNySxENgMvOmcm0dw2ncPcAuhkxP2C0rbINSuN81snpn9wJsW\nlPb1B7aa2VRvyOxhM0ulGe3T3RYjL6qPMtdzoVjd9kRt+5xztS405NIbGGtmIwhA+8xsAlDsQrfi\nONB5y1HXtjDHO+fGAGcTGg48kQC8dp4EYAzwN6+NuwmNajS5fa0V6BuBvmGPe3vTgqjYzHIAzKw7\nUOJzPc3mXSj2AvCUc+4Vb3Jg2refc24XkA+cSTDadzxwrpmtAZ4FTjGzp4DNAWgbAM65Td6/WwgN\nB44lGK8dhEYwCp1z873HLxIK+Ca3r7UCfR4wyMxyzSwJuBSY0UrbamvG13tBM4Arvd+vAF6pu0AU\neRxY5r5+1W8g2mdmXfafJWBmKcDphG4oF/Xtc879yjnX1zk3gND/tdnOue8DM4nytgGYWar3zREz\n6wiMJ3T9S9S/dgDesEqhmQ32Jp0KfE4z2tdq56Gb2ZmEjtzGAY855+5ulQ21ITObRuigUzZQDEwi\n1Ft4HugDrAcucc7t9KvG5mroQjHgE+CfRH/7DgWeIPR+jAOec85NMbPOBKB9+5nZOOAm59y5QWmb\nmfUHXiL0nkwAnnHO3R2U9gGY2WHAo0AisAa4Coinie3ThUUiIgGhg6IiIgGhQBcRCQgFuohIQCjQ\nRUQCQoEuIhIQCnQRkYBQoIuIBIQCXUQkIP4/uYpTGfmA5V0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12caf7a10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(omegas5[1:99900],omega_values_interpolated)\n",
    "plt.plot(omegas3,omegas_converged[3])\n",
    "plt.title('interpolated versus original values')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total abs error 0.520267025544\n",
      "max error 7.03401227735e-05\n"
     ]
    }
   ],
   "source": [
    "print 'total abs error', sum(abs(omega_values_interpolated - np.asarray(omegas_converged[5][1:99900])))\n",
    "print 'max error', max(abs(omega_values_interpolated - np.asarray(omegas_converged[5][1:99900])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Experimenting with 2D interpolation (could be used for four-wave mixing)\n",
    "\n",
    "We try to use 2-D cubic interpolation from a grid of length 1001 to a grid of length 10001. Time is 3.67 sec."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "omegas_dict = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for i in range(3,6):\n",
    "    omegas_dict[i] = [6.+ j*5*pow(10,-i) for j in range(1+pow(10,i))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  6.        6.00005   6.0001  ...,  10.9999   10.99995  11.     ]\n"
     ]
    }
   ],
   "source": [
    "print np.asarray(omegas_dict[5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  6.      6.005   6.01  ...,  10.99   10.995  11.   ]\n"
     ]
    }
   ],
   "source": [
    "print np.asarray(omegas_dict[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = omegas_dict[3]\n",
    "y = omegas_dict[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACYAAAAPBAMAAAB6sKBjAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAVO8Qq5l2zWaJMt0i\nu0SCRuA9AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAmUlEQVQYGWNgAAHmyM4FEALMZWBUZmBgF2De\nAyHAXJOQzwwMTQwM0yEEmMvABhS7wsAgbwAmwFww8YWB4b0AmICLMX8FiiWDiAS4GM9PBob1c0DE\nBIQYUMn6OSACIYZNLwPQeHmQHfICcL0MdxkY+g3ABEIM6OZwkJvDGRBi3ALMWgxgAizG6vTdg4F5\nWgowDEAEmAsJCBQSABtTQ/E9ZhpaAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$1001$$"
      ],
      "text/plain": [
       "1001"
      ]
     },
     "execution_count": 169,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "z = np.zeros((len(x),len(y)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for i in xrange(len(x)):\n",
    "    for j in xrange(len(y)):\n",
    "        z[i,j] = x[i]+y[j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = interpolate.interp2d(x,y,z,kind='q')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_new = omegas_dict[4]\n",
    "y_new = omegas_dict[4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAADAAAAAPBAMAAABQAFHMAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAVO8Qq5l2zWaJMt0i\nu0SCRuA9AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAmUlEQVQYGWNgAAHmyM4FSARYjIFRmYGBXYB5\nDxIBFjMJ+czA0MTAMB2JAIsxsAElrjAwyBsgCLAYmPjCwPBeAEHAJZi/AiWS4UQCXILnJwPD+jlw\nYgJCAqh4/Rw4gZDAaRQD0F55kOVQAm4Uw10Ghn4DBIGQAHowHORBKIGQ4BZg1mJAEGAJVqfvHgzM\n01KAgQgnwGKQkMQkAUHBWQq0DGvVAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$10001$$"
      ],
      "text/plain": [
       "10001"
      ]
     },
     "execution_count": 167,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(x_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@timeit\n",
    "def get_values(x_new,y_new):\n",
    "    return f(x_new,y_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "'get_values' 3.67 sec\n"
     ]
    }
   ],
   "source": [
    "z = get_values(x_new,y_new)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}