{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Evaluate block to execute LaTeX definitions) \n",
    "$ \\newcommand{\\ybar}{\\overline{y}} $\n",
    "$ \\renewcommand{\\d}[2]{\\frac{d #1}{d #2}} $\n",
    "$ \\newcommand{\\dd}[2]{\\frac{d^2 #1}{d #2^2}} $\n",
    "$ \\newcommand{\\pd}[2]{\\frac{\\partial #1}{\\partial #2}} $\n",
    "$ \\newcommand{\\pdd}[2]{\\frac{\\partial^2 #1}{\\partial #2^2}} $\n",
    "$ \\renewcommand{\\b}{\\beta} $\n",
    "$ \\newcommand{\\m}{\\mu} $\n",
    "$ \\renewcommand{\\v}[1]{\\mathbf{#1}} $\n",
    "$ \\newcommand{\\N}{\\mathcal{N}} $\n",
    "$ \\renewcommand{\\l}{\\lambda} $\n",
    "$ \\newcommand{\\ol}[1]{\\overline{#1}} $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating the Dispersion Relation for Zonostrophic Instability\n",
    "\n",
    "This is the supplementary material for the paper [Dynamics of zonal flows: failure of wave-kinetic theory, and new geometrical optics approximations](https://www.cambridge.org/core/journals/journal-of-plasma-physics/article/dynamics-of-zonal-flows-failure-of-wave-kinetic-theory-and-new-geometrical-optics-approximations/5B51BB7D026E19E8F21801568ED6EA75).  The article may also be found on arXiv [here](https://arxiv.org/abs/1604.06904).\n",
    "\n",
    "The dispersion relation for zonostrophic instability is shown in Figure 2 in the paper.  We show how to derive the general form of the dispersion relation for zonostrophic instability, and then how to compute it numerically for a given set of parameters.  Python code is included.\n",
    "\n",
    "In each case, we reduce the general form of the dispersion relation to a simpler form by using a fluctuation spectrum that consists of a thin ring in wavevector space $\\sim \\delta(k - k_f)$ centered at $k_f$.  This reduction was first considered for CE2 by Srinivasan and Young (2012); more details can be found in Parker (Ph.D. Thesis, 2014).\n",
    "\n",
    "If this document is being read as a PDF or webpage, an interactive ipython notebook can be found [here](https://github.com/jeffbparker/ZonalFlowWaveKinetic), or at a [direct link](https://github.com/jeffbparker/ZonalFlowWaveKinetic/raw/master/Supplemental_Material.ipynb) (right click, save as).  In the ipython notebook one can follow along the mathematics and run the Python code inside the notebook using the [Jupyter](http://jupyter.org/) software.  To view it non-interactively within a browser without the Jupyter software, see  [here](http://nbviewer.jupyter.org/github/jeffbparker/ZonalFlowWaveKinetic/blob/master/Supplemental_Material.ipynb).\n",
    "\n",
    "## Coordinate Convention\n",
    "This document uses the geophysical convention for coordinates.  To use the standard plasma coordinates, let $(x,y,\\ybar,\\b,U) \\mapsto (-y, x, \\ol{x}, \\kappa, -U)$.\n",
    "\n",
    "## General Form for Zonostrophic Instability\n",
    "### Asymptotic WKE\n",
    "We start from Eq. (3.2).  There is a homogeneous equilibrium, independent of $\\ybar$, at $\\N_H(k_x,k_y) = F(k_x,k_y) / 2\\b\\m$, $U=0$.  Using the form\n",
    "  \\begin{equation}\n",
    "    \\N = \\N_H + e^{iq\\ybar} e^{\\l t} \\N_1(k_x, k_y),\n",
    "  \\end{equation}\n",
    "  \\begin{equation}\n",
    "    U = e^{iq\\ybar} e^{\\l t} U_1,\n",
    "  \\end{equation}\n",
    "we linearize Eq. (3.2) for perturbations about the equilibrium.  We obtain\n",
    "  \\begin{equation}\n",
    "    \\l \\N_1 = iqk_x U_1 \\pd{\\N_H}{k_y} - \\frac{2i\\b q k_x k_y}{\\ol{k}^4} \\N_1 - 2\\m \\N_1\n",
    "  \\end{equation}\n",
    "  \\begin{equation}\n",
    "    (\\l + \\m) U_1 = iq \\int \\frac{d\\v{k}}{(2\\pi)^2} \\frac{k_x k_y}{\\ol{k}^4} \\b \\N_1.\n",
    "  \\end{equation}\n",
    "  \n",
    "The first equation is solved for $\\N_1$ in terms of $U_1$ and then substituted into the second equation.  This procedure yields a single nonlinear equation for the eigenvalues $\\l$ corresponding to zonostrophic instability.  One finds\n",
    "  \\begin{equation}\n",
    "    \\l + \\m = -q^2 \\int \\frac{d\\v{k}}{(2\\pi)^2} \\frac{k_x^2 k_y}{(\\l + 2\\m)\\ol{k}^4 + 2i\\b q k_x k_y} \\b \\pd{\\N_H}{k_y}\n",
    "  \\end{equation}\n",
    "  \n",
    "When $q$ is large, $\\l \\sim q$.\n",
    "\n",
    "### Zonostrophic Instability in CE2-GO\n",
    "We follow the same procedure, starting from Eq. (4.1).  After linearizing about the homogeneous equilibrium $W_H(k_x,k_y) = F(k_x,k_y) / 2\\m$, we find\n",
    "  \\begin{equation}\n",
    "    \\l W_1 = iqk_x U_1 \\pd{}{k_y} \\left[ \\left( 1 - \\frac{q^2}{\\ol{k}^2} \\right) W_H \\right] - \\frac{2i\\b q k_x k_y}{\\ol{k}^4} W_1 - 2\\m W_1,\n",
    "  \\end{equation}\n",
    "  \\begin{equation}\n",
    "    (\\l + \\m)U_1 = iq \\int \\frac{d\\v{k}}{(2\\pi)^2} \\frac{k_x k_y}{\\ol{k}^4} W_1.\n",
    "  \\end{equation}\n",
    "  \n",
    "Solving for $W_1$ and substituting into the second equation yields\n",
    "  \\begin{equation}\n",
    "    \\l + \\m = -q^2 \\int \\frac{d\\v{k}}{(2\\pi)^2} \\frac{k_x^2 k_y}{(\\l + 2\\m)\\ol{k}^4 + 2i\\b q k_x k_y} \\pd{}{k_y} \\left[ \\left( 1 - \\frac{q^2}{\\ol{k}^2} \\right) W_H \\right]\n",
    "  \\end{equation}\n",
    "\n",
    "### Zonostrophic Instability in the WKE\n",
    "We follow the same procedure, starting from Eq. (4.3).  We linearize about the homogeneous equilibrium $\\N_H(k_x,k_y) = F(k_x,k_y) / 2\\m \\b$, and find\n",
    "\n",
    "  \\begin{equation}\n",
    "    \\l \\N_1 = iqk_x U_1 \\left( 1 - \\frac{q^2}{\\ol{k}^2} \\right) \\pd{\\N_H}{k_y} - \\frac{2i\\b q k_x k_y}{\\ol{k}^4} \\N_1 - \\frac{q^2 F}{\\b^2} U_1 - 2 \\m \\N_1,\n",
    "  \\end{equation}\n",
    "  \\begin{equation}\n",
    "    (\\l + \\m) U_1 = iq \\int \\frac{d\\v{k}}{(2\\pi)^2} \\frac{k_x k_y}{\\ol{k}^4} (\\b \\N_1 - U_1'' \\N_H).\n",
    "  \\end{equation}\n",
    "  \n",
    "Under the assumption that $F$ and $\\N_H$ are even in $k_x$ and in $k_y$, the integral over $U_1'' \\N_H$ vanishes.  Solving for $\\N_1$ gives\n",
    "  \\begin{equation}\n",
    "    \\N_1 = \\frac{iqk_x U_1 (1 - q^2/\\ol{k}^2)}{\\l + 2\\m + 2i\\b q k_x k_y / \\ol{k}^4} \\pd{\\N_H}{k_y} - \\frac{q^2 F U_1}{\\b^2 (\\l + 2\\m + 2i\\b q k_x k_y / \\ol{k}^4)}\n",
    "  \\end{equation}\n",
    "Plugging this into the zonal flow equation yields the dispersion relation,\n",
    "  \\begin{equation}\n",
    "    \\l + \\m = -q^2 \\int \\frac{d\\v{k}}{(2\\pi)^2} \\left[ \\frac{\\b k_x^2 k_y (1 - q^2/\\ol{k}^2)}{(\\l+2\\m) \\ol{k}^4 + 2i\\b q k_x k_y} \\pd{\\N_H}{k_y} + \\frac{qF}{\\b} \\frac{i k_x k_y}{(\\l + 2\\m) \\ol{k}^4 + 2i\\b q k_x k_y} \\right]\n",
    "  \\end{equation}\n",
    "\n",
    "The last term here has no analog in the other models and its presence is a consequence of the invalidity of this model --- wave action is not conserved during this instability.\n",
    "\n",
    "### Zonostrophic Instability in CE2\n",
    "The same procedure can be applied to Eq. (2.3) with a little more algebraic complexity.  See Srinivasan and Young (2012) or Parker (PhD Thesis, Section 3.2) for the derivation.  Using Eqs. (3.25) and (3.26) of Parker (PhD Thesis), and taking $\\ol{q}^2 = q^2$ (corresponding to $L_d^{-2} = 0$ for the zonal flows) and the limit of zero viscosity $\\nu=0$, the dispersion relation is\n",
    "  \\begin{equation}\n",
    "    \\l + \\m = -q \\int \\frac{d\\v{k}}{(2\\pi)^2} \\frac{k_x^2 k_y}{(\\l + 2\\m) \\ol{h}^2_1 \\ol{h}^2_2 + 2i\\b qk_xk_y } \\left[ \\left( 1 - \\frac{q^2}{\\ol{h}_1^2} \\right) W_H(k_x, k_y + \\tfrac{1}{2}q) - \\left(1 - \\frac{q^2}{\\ol{h}_2^2} \\right) W_H(k_x, k_y - \\tfrac{1}{2}q) \\right],\n",
    "  \\end{equation}\n",
    "\n",
    "\n",
    "where $\\ol{h}_{1,2}^2 = k_x^2 + (k_y \\pm \\tfrac{1}{2} q)^2 + L_d^{-2}$.  After a Taylor expansion for small $q$, this equation reduces to the CE2-GO dispersion relation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reduction for Thin-ring, isotropic forcing\n",
    "A thin-ring forcing isotropic in wavevector space is a convenient simplification.  We use\n",
    "  \\begin{equation}\n",
    "    F(k) = 4\\pi \\varepsilon k_f \\delta(k-k_f).\n",
    "  \\end{equation}\n",
    "If $L_d^{-2} = 0$, then $\\varepsilon$ is equal to the energy input into the system by the forcing.  With this forcing, the dispersion relations above can be simplified by converting the integrals to polar coordinates.  One is left with only one integral, the angle integral, that must be computed numerically.\n",
    "\n",
    "For some of these computations we will integrate by parts, and we need the relation\n",
    "\n",
    "  \\begin{equation}\n",
    "    \\pd{}{k_y} \\frac{k_y}{(\\l + 2\\m) \\ol{k}^4 + 2i\\b qk_x k_y} = \\frac{(\\l + 2\\m) \\ol{k}^2 (\\ol{k}^2 - 4k_y^2)}{[(\\l + 2\\m) \\ol{k}^4 + 2i\\b qk_x k_y]^2}\n",
    "  \\end{equation}\n",
    "\n",
    "### Asymptotic WKE\n",
    "First we apply integration by parts in $k_y$ for the dispersion relation, then substitute $\\N_H = F/2\\m \\b$.  This yields\n",
    "\n",
    "  \\begin{equation}\n",
    "    \\l + \\m = q^2 \\int \\frac{d\\v{k}}{(2\\pi)^2} \\frac{F}{2\\m} \\frac{(\\l + 2\\m) \\ol{k}^2 k_x^2 (\\ol{k}^2 - 4k_y^2)} {[(\\l + 2\\m) \\ol{k}^4 + 2i\\b qk_x k_y]^2}\n",
    "  \\end{equation}\n",
    "\n",
    "For isotropic $F$, we use polar coordinates with $k_x = k\\sin \\phi$, $k_y = -k \\cos \\phi$.  The dispersion relation becomes\n",
    "\n",
    "  \\begin{equation}\n",
    "    \\l + \\m = q^2 \\int_0^\\infty dk\\, \\frac{F(k)}{4\\pi \\m} k^2 (\\l + 2\\m)(1+m) \\times \\int_0^{2\\pi} \\frac{d\\phi}{2\\pi} \\frac{\\sin^2 \\phi (1 + m - 4\\cos^2 \\phi)}{ [(\\l + 2\\m)k^2 (1+m)^2 - 2i\\b q \\cos\\phi \\sin\\phi]^2}\n",
    "  \\end{equation}\n",
    "\n",
    "where $m = (k L_d)^{-2}$.  Substituting the thin-ring forcing $F = 4\\pi \\varepsilon k_f \\delta(k-k_f)$, we obtain the final form of the dispersion relation,\n",
    "\n",
    "  \\begin{equation}\n",
    "    \\l + \\m = q^2 \\frac{\\varepsilon}{\\m} k_f^4 (\\l + 2\\m) (1 + m_f) \\int_0^{2\\pi} \\frac{d\\phi}{2\\pi} \\frac{(1 + m_f - 4\\cos^2\\phi) \\sin^2 \\phi}{[(\\l + 2\\m)k_f^2 (1 + m_f)^2 - 2i\\b q \\cos\\phi \\sin \\phi]^2}\n",
    "  \\end{equation}\n",
    "\n",
    "with $m_f = (k_f L_d)^{-2}$.\n",
    "\n",
    "Now, this is a nonlinear equation that we can solve for $\\l$.  For most choices of the spectrum W_H, an unstable perturbation with $\\text{Re } \\l > 0$, has an eigenvalue $\\l$ that is purely real.  It can be seen that if $\\l$ is real, the imaginary part of the $\\phi$ integral vanishes, so that one can work only with the real part.  Working only with real $\\l$ rather than complex $\\l$ has the advantage that robust one-dimensional root finders can be used.  For the integrand above, we have\n",
    "    \\begin{equation}\n",
    "    \\frac{A}{(B + iC)^2} = \\frac{A(B^2 - C^2)}{(B^2 + C^2)^2} + \\text{imag. part}\n",
    "    \\end{equation}\n",
    "where $A$, $B$, and $C$ are assumed real.  If $\\l$ is real, the imaginary part vanishes after integration, so we drop this term.  We do indeed find real solutions for $\\l$.\n",
    "\n",
    "Below is Python code to solve the dispersion for zonostrophic instability within the Asymptotic WKE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Some python setup\n",
    "from __future__ import division\n",
    "import numpy as np\n",
    "import scipy.optimize\n",
    "import scipy.integrate\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "matplotlib.style.use('ggplot')\n",
    "plt.rcParams['figure.figsize'] = (8.0, 8.0)\n",
    "font = {'family' : 'normal',\n",
    "        'weight' : 'bold',\n",
    "        'size'   : 14}\n",
    "plt.rc('font', **font)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x19bca860>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiAAAAIACAYAAAC7JdcqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XtAVGXCBvDnPcKAwCAoEmJ5AUULQ8mUwluhW263tShS\ndtvspqulhZdvt7yjaK1572a2rVpKkpvlbbVaa0UNUFNJvIOaQoIIwiByGc77/cHnfJkXGJgzZ2Z4\nfn/NMHPOPPNm+HjOed8jpJQSRERERHak6B2AiIiImh4WECIiIrI7FhAiIiKyOxYQIiIisjsWECIi\nIrI7FhAiIiKyOxYQIiIisjs3vQPY0vz585Geng4AiIqKwrhx4276/p07d2L9+vXIzc2FwWBAt27d\nEB8fj6CgIHvEJSIiarJc5gjId999Zykf9bFt2zYsXrwYp06dgr+/P6SUSE9Px9SpU1FSUlKvfWRl\nZTU0LtUTx1h7HGP74Dhrj2OsPVuOsUsUkHPnzmH58uUICwtDy5Yt63y/2WzG6tWrAQD33HMPlixZ\ngvnz58PT0xMlJSVYt25dvT6Xf9i1xzHWHsfYPjjO2uMYa48F5FdUVcWSJUugKArGjh0LRan7K2Vn\nZ8NkMgGoPVUDAP7+/ggLCwMA7N+/X7vARERE5PwFJCUlBSdOnMCLL76I1q1b12ubCxcuWB77+vpa\nHrdo0QIAUFhYaNuQREREdBWnvgg1JycHX331Ffr164c+ffo0en913ZcvKyvrqsNPcXFxjf5MujmO\nsfY4xvbBcdYex1h7cXFxSElJsTwPDw9HeHh4g/bl1AXk559/hqqqSEtLQ0ZGBgCgsrISALB79278\n+c9/xtKlS9G8efOrtmvVqpXlcWlp6TWPAwICrvt51xvovLy8xn8RuiGj0Wg5XUba4BjbB8dZexxj\n7QUHB9us6Dn9KRgAqK6uRmVlpaV8ALXXhlRVVUFKicTERCQkJCA5ORkAEBoaCh8fHwBAWloaAKCo\nqAjHjh0DAPTo0cPO34CIiKhpceojIPfddx/uu+++q3728ssvo7Cw8Kp1QPLz81FYWIji4mIAgJub\nG4YNG4Zly5YhPT0dY8aMgclkQkVFBXx9fTFkyBB7fxUiIqImxakLyM0IIW76s0GDBsHT0xMbNmxA\nbm4u3N3dERUVhfj4ePj5+dkzKhERUZMjZF1XXtJN8RoQbfGcrvY4xvbBcdYex1h7wcHBNtuXS1wD\nQkRERM6FBYSIiIjsjgWEiIiI7I4FhIiIiOyOBYSIiIjsjgWEiIiI7I4FhIiIiOyOBYSIiIjsjgWE\niIiI7I4FhIiIiOyOBYSIiMjG5OkTekdweCwgRERENqR+ux7q+29Cll/SO4pDc9m74RIREdmbuuVf\nkNu3Qpk4B8LLW+84Do0FhIiIyAbUjWsg07+vLR/+rfSO4/BYQIiIiBpBSgm5fjXk3l1QJsyGaOGv\ndySnwAJCRETUQFJKyH+tgMz6EcrE2RDGFnpHchosIERERA0gpYRM+QfksSwo42dB+PjqHcmpsIAQ\nERFZSaoqZPKHkKdPQBk/E8LLR+9ITocFhIiIyApSVSE/fQ/ylzNQEhIhmnvpHckpsYAQERHVk1Rr\nIJcvhrxwHsqr0yE8m+sdyWmxgBAREdWDrKmB/HgBpKkEythpEB4eekdyaiwgREREdZDmaqjL5gFV\nlVDGTIFwN+gdyemxgBAREd2ErK6GuvQtAIAy+g0Id3edE7kGFhAiIqIbkFWVUN9/EzB4QHlpPIQb\ny4etsIAQERFdh6yshPruLAhjC4jnEyCaNdM7kkthASEiIvoNWXEZ6juzIFoGQAwfC6GwfNgaCwgR\nEdGvyMvlUBfPgGhzG8SfRkMoit6RXBJHlYiI6P/I8jKoC6ZC3NqR5UNjPAJCREQEQJaVQl0wDSIs\nHCLuBQgh9I7k0lhAiIioyZOlF2uPfHTrCfHEn1k+7IAFhIiImjRZUgx13mSIntEQj8WzfNgJCwgR\nETVZsvgC1PmTIaLug/LI03rHaVJYQIiIqEmSF85DnTcJYsBgKA8+oXecJocFhIiImhx5/lztaZdB\nj0EZ9JjecZokFhAiImpSZH4e1PlTIH4fC+W+h/SO02SxgBARUZMhfzkDdf5UiMeGQen3gN5xmjQW\nECIiahJk7unadT5in4Vy7/16x2nyWECIiMjlyZ9zapdXj3sBSu/+eschsIAQEZGLk6eOQ12cCOWP\noyB6Rusdh/6P0xeQbdu24ZtvvkF+fj4qKipgNBrRsWNHDBkyBF27dr3hdp9//jnWrl173deSk5Oh\ncP1/IiKnJ7OPQH03CcqzYyC699Y7Dv2K0xeQI0eO4MKFC2jdujVUVUVubi727duHrKwsLFiwAAEB\nATfd3tfXF7fccovluRCCq+AREbkAeSwL6gdvQnn+NYhuPfWOQ7/h9AVkxIgRcHP7/6+xbds2LF26\nFFVVVcjJyamzgERGRmL06NFaxyQiIjuSRzKhfjgXyovjIe7ooXccug6nLyBubm44fPgwPv30U1RW\nViIvLw8AYDAYEBISUuf26enp2LVrF7y9vdGxY0cMHToUHTp00Dg1ERFpRWbtg/qP+VBG/hWiSze9\n49ANOH0BAYCysjKcOHHC8rxFixYYN25cnUc/FEWBn58fmjVrZjl1c/DgQcyaNYslhIjICcnM3VCX\nL4Yy+nWITnfoHYduQkgppd4hbKWkpARffPEFtmzZglatWmHmzJlo1arVdd977tw5GI1GeHt7AwAy\nMzORlJQEAIiJicHIkSOv2SYrKwtZWVmW53FxcTCZTBp8E7rCYDCgqqpK7xgujWNsHxxn7cl9aSj9\n4O/wnpgEt0636x3HJRmNRqSkpFieh4eHIzw8vEH7cqkCAgCXL1/G8OHDAQCPP/44hg4dWu9tX3jh\nBZSVlSEiIgKTJk2q1zZXTvmQNoxGI0uexjjG9sFx1pbcswPys2UQY6ZCtA/VO47LCg4Ottm+nHqu\naVVVFf7zn/9c9a+KPXv2WB5XVFQAABITE5GQkIDk5GTLa5s2bUJxcbHleWZmJsrKygAArVu31jo6\nERHZiJr2PdTPPoLPG3NZPpyIU18DYjab8eGHH+Ljjz9GUFAQqqurkZ+fD6D24tT+/WtXu8vPz0dh\nYeFVhWPz5s1YuXIlAgIC4OHhgdzcXACAp6cnHn74Yft/GSIispq68z+QX34CZVwimrUPBXiUyWk4\ndQFxd3dHv379cPz4cZw/fx5msxn+/v4ICwvDY489ds0smF+v7/HEE08gLS0NZ86cQUFBAQIDA9Gl\nSxfExsaiTZs29v4qRERkJXX7FsiNKVDGz4IIulXvOGQll7sGxN54DYi2eN5cexxj++A425a6bSPk\n119CGTcTIrD2H40cY+3Z8hoQpz4CQkRETY/69ZeQ322CMiEJIuCWujcgh8QCQkRETkPd/Dnkzv9A\nmTgboiUnDDgzFhAiInJ4UkrIDZ9B7k6tLR9+LfWORI3EAkJERA5NSgn55aeQBzKgTEyC8PXXOxLZ\nAAsIERE5LCkl5Np/Qh4+AGV8EoTRV+9IZCMsIERE5JCklJCfLYPMPlI71dbbqHcksiEWECIicjhS\nVSFXfQB59mTtVFsvb70jkY2xgBARkUORag3kyncgC36BkjADwtNL70ikARYQIiJyGLKmBvKfCyEv\nFkF5dTqEh6fekUgjLCBEROQQpNkM+Y/5kOWXoIyZCuHhoXck0hALCBER6U6aq6F+OBcwm6G8MgnC\n3aB3JNKYoncAIiJq2mR1FdT35gAAlNGvs3w0ETwCQkREupFVlVDfnQ3h5Q3xwjgIN/611FTwvzQR\nEelCVlZAXTITwr8VxPBXIZo10zsS2RELCBER2Z2sKIe6OBEisA3En1+BUFg+mhoWECIisitZfgnq\n4hkQbTtA/PEvEAovR2yKWECIiMhu5CUT1AXTIEK7Qgx9CUIIvSORTlhAiIjILqSpFOqCKRC3d4d4\n8jmWjyaOBYSIiDQnS4uhzp8K0b03xJA/sXwQCwgREWlLXrwAdd4UiF79IB4dyvJBAFhAiIhIQ7Lo\nPNR5kyH6DILy0FN6xyEHwgJCRESakIX5UOdPgbjvISgPDNE7DjkYFhAiIrI5WfBLbfl4YAiUmEf0\njkMOiAWEiIhsSp47W3vB6SNxUPoP1jsOOSgWECIishmZ9zPUBVMhhvwJSp9BeschB8YCQkRENiHP\nnoS6cDrEk89Buec+veOQg2MBISKiRpOns2uXVx86AkqvvnrHISfAAkJERI0iTx6DumQmlD+Nhrjr\nXr3jkJNgASEiogaTJw5BfW8OlGfHQnTvpXccciIsIERE1CDy6EGoS9+C8sI4iPBIveOQk2EBISIi\nq8lD+6F+NA/KSxMgbu+udxxyQiwgRERkFXlwL9SPF0L5y98gwsL1jkNOigWEiIjqTR7IgLpiCZSX\nJ0GEdtU7DjkxFhAiIqoXuXcX1FXvQxk7FaJDZ73jkJNjASEiojqpGdshU/4B5bUZEO1C9I5DLoAF\nhIiIbkrdtQ3yi5VQEhIh2rbXOw65CBYQIiK6ITX1a8j1yVDGz4Joc6vecciFsIAQEdF1qd9vhvz3\nWigTkiBuCdY7DrkYpy8g27ZtwzfffIP8/HxUVFTAaDSiY8eOGDJkCLp2vfkV2jt37sT69euRm5sL\ng8GAbt26IT4+HkFBQXZKT0TkmNRv10N+ux7KhNkQrfk7kWxP0TtAYx05cgQXLlxA69at0bZtW5hM\nJuzbtw9JSUkoLCy84Xbbtm3D4sWLcerUKfj7+0NKifT0dEydOhUlJSV2/AZERI5F3foF5LaNUCay\nfJB2nP4IyIgRI+Dm9v9fY9u2bVi6dCmqqqqQk5ODgICAa7Yxm81YvXo1AOCee+5BQkICiouL8dpr\nr6GkpATr1q3D8OHD7fUViIgchrpxDWTa91AmzoHwb6V3HHJhTl9A3NzccPjwYXz66aeorKxEXl4e\nAMBgMCAk5PpTxbKzs2EymQAAUVFRAAB/f3+EhYUhMzMT+/fvt094IiIHIaWEXL8acu+u2iMfLfz1\njkQuzukLCACUlZXhxIkTluctWrTAuHHjrnv0AwAuXLhgeezr63vVdgBueuqGiMjVSCkhv1gJ+dOe\n2gtOff30jkRNgNNfAwIAvXr1wpo1a/Dhhx9i8ODBKCkpweLFi68qGvUhpdQoIRGRY5JSQqb8A/LQ\nfpYPsiuXOAJyRYsWLTB06FBs2bIFFy5cwDfffIOhQ4de875Wrf7/vGZpaek1j2905CQrKwtZWVmW\n53FxcTAajbaKT9dhMBg4xhrjGNuHI46zVFVcXr4YNSePwXvaQijePnpHahRHHGNXlJKSYnkcHh6O\n8PCG3ZDQqQtIVVUVUlNT0a9fPxgMBgDAnj17LK9XVFQAABITE1FcXIzevXtj2LBhCA0NhY+PD8rK\nypCWlobo6GgUFRXh2LFjAIAePXpc9/OuN9BXriUhbRiNRo6xxjjG9uFo4yxVFfLT9yDzfoby6nRc\nUiXgQPkawtHG2BUZjUbExcXZZF9OXUDMZjM+/PBDfPzxxwgKCkJ1dTXy8/MB1F6c2r9/fwBAfn4+\nCgsLUVxcbHlt2LBhWLZsGdLT0zFmzBiYTCZUVFTA19cXQ4YM0e07ERFpTao1kMuXQF4oqL23i2dz\nvSNRE+TUBcTd3R39+vXD8ePHcf78eZjNZstslscee+yaWTBCCMvjQYMGwdPTExs2bEBubi7c3d0R\nFRWF+Ph4+PnxHCgRuSZZUwP58QJIUwmUsdMgPDz0jkRNlJC88rJRrkz7JW3wkKr2OMb24QjjLM1m\nqMveBqoqoIx6HcLgWuXDEcbY1QUH225Jfqc+AkJERPUjq6uhLn0LAKCMngTh7q5zImrqWECIiFyc\nrKqE+v6bgMEDykvjIdxYPkh/LCBERC5MVlZCfXcWhI8vxAvjIJo10zsSEQAWECIilyUrLkN9ZxaE\nfwDEc2MhFJYPchwsIERELkheLoe6eAZE0K0Qz4xm+SCHwwJCRORiZHkZ1IXTIdqHQgwbCaG4xF03\nyMWwgBARuRBZVgp1wTSIzndAPP3iVesfETkSFhAiIhchTSVQ50+BCL8LIvZZlg9yaCwgREQuQJYU\nQ503GaJnNMRj8Swf5PBYQIiInJwsvgB1/mSIqAFQHrn2DuBEjogFhIjIickL56HOmwTR/0Eog2P1\njkNUbywgREROSp4/V3vaZdCjUAb9Qe84RFZhASEickIyP6/2tMvgJ6Hc/5DecYisxgJCRORk5C9n\na2e7PDYMSr8H9I5D1CAsIERETkTmnq5d5+OJP0OJjtE7DlGDsYAQETkJ+XNO7fLqTz0PJWqA3nGI\nGoUFhIjICchTx6EuToTyx1EQPaP1jkPUaCwgREQOTmYfgfpuEpQ/vwLRI0rvOEQ2wQJCROTA5LEs\nqB+8CeW51yDu7Kl3HCKbYQEhInJQ8kgm1A/nQnlxPMQdPfSOQ2RTLCBERA5IZu2D+o/5UEb+FaJL\nN73jENkcCwgRkYORmbuhLl8MZfTrEJ3u0DsOkSZYQIiIHIjclwb1k3ehvDIZIqSL3nGINMMCQkTk\nIOSeHVCTP4Ty6nSI9qF6xyHSFAsIEZEDUNO+h1z7TyivzYC4raPecYg0xwJCRKQzded/IL/8BErC\nTIi27fSOQ2QXLCBERDpSt2+B3JgCZfwsiKBb9Y5DZDcsIEREOlG3bYTcug7KhFkQgcF6xyGyKxYQ\nIiIdqF9/CfndJigTZ0ME3KJ3HCK7YwEhIrIzdfPnkDu/rS0fLVvrHYdIFywgRER2IqWE3PAZ5O7U\n2vLh10rvSES6YQEhIrIDKSXkl59C7k+HMjEJwtdf70hEumIBISLSmJQSFZ9+AJm5B8qEJAhjC70j\nEemOBYSISENSSsjPlsF86njtbBdvo96RiBwCCwgRkUakqkKu+gDy7En4TH4bZTVS70hEDkPROwAR\nkSuSag3kyiWQeT9DSZgB4eWjdyQih8IjIERENiZraiD/uRDyYhGU16ZDeHjqHYnI4bCAEBHZkDSb\nIf8xH7L8EpQxUyE8PPSOROSQGlVAKioqAACenmz3RETSXA31w7mA2QzllUkQ7ga9IxE5LKsKyMGD\nB5GRkYGjR4/i7NmzMJvNtTtxc8Ott96KsLAw9O7dG3feeacmYX9rw4YN2LdvH/Ly8mAymeDr64vO\nnTvjySefRLt2N76j5Oeff461a9de97Xk5GQoCi+NISLryOoqqB+8BSjNoIx+HcLNXe9IRA6tzgJi\nNpvx7bffYsOGDSgsLISPjw86duyIAQMGwMfHB1JKXLp0CQUFBdi1axe+/vprBAQE4JFHHsHvfvc7\nuLlpd5Zny5YtKCwsRFBQEJo3b468vDykp6fjwIEDmDdvHgICAm66va+vL2655f/vwSCEgBBCs7xE\n5JpkVSXUd2dDNPeCeHE8hIa/94hcRZ3/l4wdOxbV1dUYMGAAoqOjERISctP35+TkYNeuXVi3bh02\nbtyId99912ZhfysmJgb9+vVDYGAgAGDjxo345JNPUFFRgYyMDDz00EM33T4yMhKjR4/WLB8RuT5Z\nWQH1nVkQLfwhnnsNolkzvSMROYU6C8gf/vAH3H///TAY6ncuMyQkBCEhIXj66aexbdu2Rge8mdjY\n2KueR0REWB7X58hLeno6du3aBW9vb3Ts2BFDhw5Fhw4dbB2TiFyUrCiHujgRonUbiGdfgVBYPojq\nq86LHR588MF6l49fc3d3x4MPPtigUA21ceNGALWnVu69996bvldRFPj5+SEwMBAXL17Evn37MHny\nZJw6dcoOSYnI2cnyS1AXTodocxvEs2NYPoisJKSUTr80n9lsxgcffIDU1FR4eXnh9ddfR1hY2A3f\nf+7cORiNRnh7ewMAMjMzkZSUBKD2tM7IkSOvu11WVhaysrIsz+Pi4mAymWz4Tei3DAYDqqqq9I7h\n0jjG1lPLTLg0eyKadb4DzYePqde1Yxxn7XGMtWc0GpGSkmJ5Hh4ejvDw8Abtq15XSh04cADt27eH\nn5+f5Wdms/m6pzl++eUXZGZm2u3oh8lkwty5c3H06FG0bNkSf/vb39C+ffubbhMUFHTV84iICPj4\n+KCsrAyFhYU33O56A80Coi2j0cgx1hjH2DrSVAp1wRSIrhGQscNRVlZWr+04ztrjGGvPaDQiLi7O\nJvuq13zT2bNn4+DBg5bnJpMJf/zjH6/62RUnTpzAxx9/bJNwdTl79izeeOMNHD16FB06dEBSUtI1\n5SMxMREJCQlITk62/GzTpk0oLi62PM/MzLT8EmndurVdshOR85GlF6HOmwTRrSfEU89z1hxRIzj1\nXLG3334bBQUFAICamhrMmzfP8trAgQMRExOD/Px8FBYWXlU4Nm/ejJUrVyIgIAAeHh7Izc0FULug\n2sMPP2zfL0FETkFeLII6bzJEr74Qjw5j+SBqJKcuIFcWQgOAM2fOXPVaZGTkVc9//cviiSeeQFpa\nGs6cOYOCggIEBgaiS5cuiI2NRZs2bbQNTURORxYV1paP6BgoD9vm8DNRU+fUBeSdd96p8z3XW4dk\n4MCBGDhwoBaRiMjFyMJ8qPOnQNz3eygPPK53HCKX4dQFhIhIS7Lgl9ry8bshUAY+onccIpdS7wKS\nn5+PEydOAADKy8sBALm5udfciC4/P9+G8YiI9CHPnYU6fyrEw3FQBgzWOw6Ry6l3AUlJSblq7i8A\nu812ISKyJ5n3M9QFUyGG/AlKn0F6xyFySfUqIKNGjdI6BxGRQ5BnT9aucPrkcCj33K93HCKXVa8C\nct9991m101/PTiEichbydDbUxTMgho6A0quv3nGIXFq9FiJLT0+v9w4rKiowe/bsBgciItKDPHkM\n6qLpUP44iuWDyA7qVUAWLVqEjIyMOt9XVlaGxMTEq+6XQkTk6OSJQ1CXzITy7FiIu25+I0siso16\nFZBOnTph4cKFNy0hFy9exPTp05GTk3PDm7kRETkaefQg1PfmQHk+AaJ7L73jEDUZ9Sogb7zxBjp3\n7nzDElJQUICpU6ciLy8Pr776KmJiYmwelIjI1uSh/VA/eBPKSxMgut2ldxyiJqVeBcTT0/OGJeTs\n2bOYNm0aiouL8T//8z+4914eviQixycP7oX60Twoo16HuL273nGImpx6FRAA8PDwwBtvvIGwsDBL\nCcnJycG0adNQUVGBSZMmoUePHlpmJSKyCXkgA+rHC6G8PAkiLFzvOERNkpBSSms2qKqqwpw5c3D0\n6FG4u7vDYDBg0qRJ6NChg0YRHVteXp7eEVya0WiEyWTSO4ZLa2pjLPfugrrqfShjpkJ07Gy3z21q\n46wHjrH2goODbbaveq0DcmUJ9itiY2Px/vvvo7S0FM888wzMZvM17+nUqZPNQhIR2YKasR1yzUdQ\nXpsO0S5U7zhETVq9CsikSZNu+Nr17jYLAGvWrGlYIiIiDai7tkF+sRJKQiLErR30jkPU5HEpdiJy\neWrq15Drk6GMnwnR5ja94xARNFqKnYjIUajfb4b891ooE5IgbrHd+Wsiapx6HwHp2rUrunTpgttv\nvx3t27fXOhcRUaOp366H/HY9lAmzIVoH6R2HiH6lXgWkVatWyMjIwK5duwAAXl5eCAsLQ9euXdG1\na1d06tQJ7u7umgYlIrKGuvULyP9ugTJxNkSrQL3jENFv1HsabnV1NU6cOIGjR4/iyJEjOHr0KMrL\nywEAbm5uCAkJsRSSLl26wMfHR9PgjoLTcLXFaXXac8UxVjeugUz7Hsr4WRD+rfSOA8A1x9nRcIy1\nZ8tpuFavA/JrP//8s6WMHDlyBIWFhQAARVGQnJxss5COjAVEW/yFoj1XGmMpJeT61ZB7d0EZNxPC\nr6XekSxcaZwdFcdYe3ZfB+RG2rVrh9tuuw1hYWHo1KkTdu3ahWPHjkFVVVvlIyKqFykl5BcrIX/a\nU3vBqa+f3pGI6CasLiBVVVU4fvw4jhw5giNHjuD48eO4fPkyfHx8EBYWhmHDhqFr165aZCUiui4p\nJWTKPyCPHawtHz6+ekciojrUq4BkZGRYTrXk5ORAVVUEBwejS5cuiI6ORpcuXWx6WIaIqL6kqkIm\nfwh5+gSUcbMgvJvG9WdEzq5eBWTevHlo1qwZoqOjERsbi7CwsCZzkSkROS6pqpCfvgeZ93PtCqfN\nvfSORET1VK8C0rVrV2RnZyM1NRWHDx9Gly5dLDNebrvtNgghtM5JRHQVqdZALl8CeSEfymszIDyb\n6x2JiKxQ71kwZrMZOTk5lms/jh07BpPJhObNm6Nz586W6bedO3eGh4eH1rkdBmfBaItXtWvPGcdY\n1tRAfrwA0lQC5eVJEB6eekeqkzOOs7PhGGvPYabhnj17FkePHrVMw83Pz4eiKOjQoQPmzJljs5CO\njAVEW/yFoj1nG2NpNkNd9jZQVQFl1OsQBuf4B4+zjbMz4hhrz2Gm4d56660IDg5Ghw4d0K5dO6Sn\np+PYsWPIycmxVT4iIgtZXQ116VsAAGX0JAiuwEzktKwuIBUVFTh27JhlVszx48dRWVkJADAYDAgP\nD+c0XCKyOVlVCfX9NwGDB5SXxkO4sXwQObN6FZC0tDTLtR+nT5+2LDTm6+uLiIgIywWpISEhUBRF\n08BE1PTIykqo786C8PGFeGEcRLNmekciokaqVwFZsGABACAoKAj9+/e3FI42bdpoGo6ISFZchvrO\nLAj/AIjnxkIoLB9ErqBeBWT8+PHo2rUrfH25uiAR2Y+8XA518QyIoFshnhnN8kHkQupVQHr37q11\nDiKiq8jyMqgLp0O0D4UYNhKCp3eJXEqjZsEQEWlBlpVCXTANovMdEE+/yMUOiVwQCwgRORRZehHq\ngqkQ4XdBxD7L8kHkolhAiMhhyJJiqPMmQ/SMhngsnuWDyIWxgBCRQ5DFF6DOnwwRNQDKI0P1jkNE\nGmMBISLdyQvnoc6bBNH/QSiDY/WOQ0R2YPVl5RUVFSgsLLzh64WFhZaVUYmI6iLPn4M693WImIdZ\nPoiaEKuPgKxYsQLZ2dn4+9//ft3X586di86dO+PFF19sdLi6bNiwAfv27UNeXh5MJhN8fX3RuXNn\nPPnkk2h20uW2AAAgAElEQVTXrt1Nt925cyfWr1+P3NxcGAwGdOvWDfHx8QgKCtI8NxHVkvl5tadd\nBj8J5f6H9I5DRHZk9RGQzMzMm64L0qtXLxw4cKBRoepry5YtyMrKgoeHBwIDA1FUVIT09HRMmTLl\npkdptm3bhsWLF+PUqVPw9/eHlBLp6emYOnUqSkpK7JKdqKmTv5yB+vYkiEeGsnwQNUFWF5CLFy/C\n39//hq/7+fmhqKioUaHqKyYmBkuWLMGiRYuwYMECPPPMMwBqTxNlZGRcdxuz2YzVq1cDAO655x4s\nWbIE8+fPh6enJ0pKSrBu3Tq7ZCdqymTuaajzpkA8/gyUfg/oHYeIdGB1AfH19cXZs2dv+PrZs2fh\n7e3dqFD1FRsbi8DAQMvziIgIy2M3t+ufXcrOzobJZAIAREVFAQD8/f0RFhYGANi/f79WcYkIgPw5\np3adj6eegxIdo3ccItKJ1QUkMjIS3377LbKzs6957cSJE/j2228RGRlpk3DW2rhxI4DaknTvvfde\n9z0XLlywPP71vW1atGgBADc9dUNEjSNPHYe6cBqUYSOhRA3QOw4R6cjqi1Dj4uKwb98+TJ48GZGR\nkbj11lsBAGfOnMH+/fvRokULPP300zYPejNmsxkffPABUlNT4eXlhYkTJ8JoNFq1Dyllne/JyspC\nVlaW5XlcXJzVn0PWMRgMHGON2WuMzceycGnJTHiPnAj3u/to/nmOhn+Wtccxto+UlBTL4/DwcISH\nhzdoP1YXED8/P8yZMwerVq3C7t27sXfvXgBA8+bN0a9fP8THx8PPz69BYRrCZDJh7ty5OHr0KFq2\nbIm//e1vaN++/Q3f36pVK8vj0tLSax4HBATccNvrDfSV0zmkDaPRyDHWmD3GWB7LgvrBm1Ceew0V\nXSJQ0QT/m/LPsvY4xtozGo2Ii4uzyb4atBCZn58fXn75ZUgpLX9x+/r62n3Z5LNnz+Ktt95CQUEB\nOnTogL/+9a9o2bLlVe9JTExEcXExevfujWHDhiE0NBQ+Pj4oKytDWloaoqOjUVRUhGPHjgEAevTo\nYdfvQOTq5OEDUD+cC+WlCRB38P8vIqpVZwG5ck3ElSMDN7pG4tfXVvz6/Vp6++23UVBQAACoqanB\nvHnzLK8NHDgQMTExyM/PR2FhIYqLiwHUXpw6bNgwLFu2DOnp6RgzZgxMJhMqKirg6+uLIUOGaJ6b\nqKmQWfug/mM+lL/8DaJLN73jEJEDqbOAvPzyywCAVatWwc3NzfK8LmvWrGlcsnowm82Wx2fOnLnq\ntd9eCPvrozODBg2Cp6cnNmzYgNzcXLi7uyMqKsrup4+IXJnM3A11+WIoo1+H6HSH3nGIyMEIWcfV\nl99//z0AYMCAARBCWJ7X5b777mtkNOeQl5endwSXxnO62tNijOWPP0D99D0or0yGCOli0307K/5Z\n1h7HWHvBwcE221edR0B+WySaSrEgooZRd++A/OxDKK9Oh2gfqnccInJQVq8DcujQoatmj/xWaWkp\nDh061KhQROSc1LTvIdcsg5Iwg+WDiG7K6gIyY8YMZGZm3vD1gwcPYsaMGY0KRUTOR935LeS/lkNJ\nmAlxa0e94xCRg2vQNNybMZvNUBSrew0ROTF1+xbIjSlQxs+CCLpV7zhE5ATqVUDKy8tRXl5ueW4y\nma47HbesrAw7duy4Zi0OInJd6raNkF9/CWVCEkRgG73jEJGTqFcB2bRpE9auXWt5vnz5cixfvvyG\n74+Pj290MCJyfOrX6yC/21xbPgJu0TsOETmRehWQ7t27w9PTE1JKrFq1Cn369EHHjlef4xVCwMPD\nA6GhoQgJCdEkLBE5DnXz55A7/wNl4myIlq31jkNETqZeBSQsLMxyu/rKykpERUWhXbt2mgYjIsck\npYTc8Bnk7lQoE5Mg/FrVvRER0W9YfRHqU089ZXlcWlpqWQo9MDDwqtvbE5HrkVJCrvsEMnN3bfnw\n9dc7EhE5qQbNgjly5AhWrFiBnJycq37eqVMnPPPMM+jatatNwhGR45BSQq79J+ThA1DGJ0EY+Q8O\nImo4qwvIkSNHMHPmTHh6euKRRx5B27ZtAQC5ubnYvn07Zs6ciSlTprCEELkQKSXkZ8sgs4/UTrX1\nNuodiYicnNUFZM2aNWjdujVmzpwJo/HqX0KPP/44Jk+ejDVr1mDatGk2C0lE+pGqCrnqfcizp6CM\nmwnh5a13JCJyAVavGHbixAkMHDjwmvIBAD4+PoiJiUF2drZNwhGRvqRaA7liCeQvZ2qXV2f5ICIb\nsfoIiKIoqK6uvuHrZrMZQohGhSIi/cmaGsh/LoS8WFR7YzkPT70jEZELsfoISNeuXbF161bk5+df\n81pBQQG2bt2K22+/3SbhiEgf0myGXPY2pKkUypipLB9EZHNWHwGJj4/H1KlTMW7cOPTs2RPBwcEA\ngLy8POzduxdubm5cCZXIiUlzNdSlc4EaM5RXJkG4G/SOREQuSEgppbUb5ebmIjk5GQcOHEBVVRUA\nwGAwoEePHhg6dKhlZkxTkJeXp3cEl2Y0GmEymfSO4dJ+Pcayugrq+28Cbm5QRkyEcHPXOZ3r4J9l\n7XGMtXfloIMtNGgdkLZt22LChAlQVRWlpaUAAF9fX94Fl8iJyapKqO/OhvDyhnhhHISbzW+WTURk\nYXVjOHTokKV0KIoCPz8/+Pn5WcpHaWkpDh06ZNuURKQpWVkBdXEihG8LiBfHs3wQkeasLiAzZsxA\nZmbmDV8/ePAgZsyY0ahQRGQ/8nI51EXTIQICIZ57FaJZM70jEVETYPN/5pjNZp6KIXISsvwSyt6Z\nCdGmHcQf/wLB/3eJyE7qVUDKy8tRXl5ueW4ymVBYWHjN+8rKyrBjxw60bNnSdgmJSBPykgnqgmkw\n3H4n5BPDuX4PEdlVvQrIpk2bsHbtWsvz5cuXY/ny5Td8P6fhEjk2aSqFumAKxO3d0fzZMSgrK9M7\nEhE1MfUqIN27d4enpyeklFi1ahX69OmDjh07XvUeIQQ8PDwQGhqKkJAQTcISUePJ0mKo86dCdO8N\nMeRPPPJBRLqoVwEJCwtDWFgYAKCyshJRUVFo166dpsGIyPbkxQtQ502B6N0f4pGnWT6ISDdWX4T6\n1FNPAaid7bJ3716cP38eABAYGIi77roL3bp1s21CIrIJWXQe6rzJEH0GQXnoKb3jEFETZ3UBqaio\nwMKFC7Fv3z4AgLd37d0xd+/ejU2bNqFHjx5ISEiApyfvHUHkKGRhfm35uP9hKA8M0TsOEZH1BWTl\nypXYt28fYmNj8fvf/x5GoxFA7cyYzZs344svvsDKlSsxYsQIm4clIuvJgl+gzp8C8cAQKDGP6B2H\niAhAAxYi++GHHzBw4EDExcVZygdQuwb/008/jZiYGPzwww82DUlEDSPPnYX69iSIh55k+SAih2J1\nAZFSokOHDjd8/WavEZH9yLyfa0+7/CEeSv/BeschIrqK1QUkMjISe/fuveHrP/74IyIjIxsVioga\nR549WTvVNnY4lD6D9I5DRHQNq68BiY2NxaJFi/Dmm29i8ODBCAoKAgD88ssv2LJlC4qKivDMM8+g\npKTkqu1atGhhm8REdFPydDbUxTOgDBsBcXdfveMQEV2XkFJKazZ4+umnG/RBa9asadB2ji4vL0/v\nCC7NaDTCZDLpHcNpyJyjUN+ZBeWZlyEi76nXNhxj++A4a49jrL3g4GCb7atBR0C4eBGR45EnDkF9\nbw6U4WMhInrpHYeI6KasLiBxcXFa5CCiRpBHf4K69O9QXhgHEc5rsIjI8VldQIjIschD+6EuexvK\niIkQt3fXOw4RUb2wgBA5MfnTXqj/XAhl1OsQYeF6xyEiqrcGFZCsrCx89913yM/Px6VLl/Db61iF\nEJg/f75NAhLR9cn96VBXvgPl5UkQoV31jkNEZBWrC8j69euxatUqGAwGBAcH6z699vDhw/jqq6+Q\nnZ2N0tJSALU3zHvyySdvut3nn3+OtWvXXve15ORkKIrVS6QQ2Y3cuwvqqvehjJ0K0aGz3nGIiKxm\ndQHZsGEDunbtir/+9a/w8vLSIpNVTp48iQMHDqBNmzaWAmINX19f3HLLLZbnQgjO8iGHpqb/F/Lz\nj6G8NgOiXYjecYiIGsTqAlJVVYW+ffs6RPkAgP79+2PQoEFQVRXPPvus1dtHRkZi9OjRGiQjsj11\n138gv/gESkIiRNv2eschImowqwtIt27dcPr0aS2yNIiPjw8AoKKiokHbp6enY9euXfD29kbHjh0x\ndOhQ3s+GHJKa+jXk+mQo42dBtLlV7zhERI1i9YUOL7zwAg4dOoQvv/zymuXWnY2iKPDz80NgYCAu\nXryIffv2YfLkyTh16pTe0Yiuon63GXLjGigTklg+iMglWH0EpGXLlrj//vuxatUqJCcnw83N7Zpr\nJoQQ+OSTT2wWUgv9+vXDQw89BG9vbwBAZmYmkpKSUF1dja1bt2LkyJHXbJOVlYWsrCzL87i4OBiN\nRrtlbooMBkOTH+OKzWtR9c2X8Jm+CM0C29h8/xxj++A4a49jbB8pKSmWx+Hh4QgPb9gSAFYXkOTk\nZHz55Zdo2bIlQkNDHeZaEGtduYneFREREfDx8UFZWRkKCwuvu831Bpr3HdBWU7+3g7rlX5Dbt0IZ\nn4Ty5j6ABmPR1MfYXjjO2uMYa89oNNpsRXSrC8i3336Lnj17YsKECU4zVTUxMRHFxcXo3bs3hg0b\nBgDYtGkToqOj4e/vD6D2CEhZWRkAoHXr1rplJbpC3fgZZPp/oUycA+HfSu84REQ2ZXUBqampQWRk\npMOUj4yMjGtO92zevBnbt29H586dMWbMGOTn56OwsBDFxcVXvWflypUICAiAh4cHcnNzAQCenp54\n+OGH7fodiH5NSgn51SrIH3+AMmE2RAt/vSMREdmc1S3irrvuwqFDh7TI0iDl5eUoKChAQUGB5WeX\nLl1Cfn7+VYUDwFXXqjzxxBOIiIhATU0NCgoKEBgYiH79+uHNN99E27Zt7Zaf6NeklJD/WgF5IAPK\nRJYPInJdQv52HfU65ObmYtGiRQgNDUVMTAwCAgKuezRE7xVS7SUvL0/vCC6tKZ3TlVJCrvkI8vgh\nKAkzIHx87fK5TWmM9cRx1h7HWHvBwcE225fVp2DGjRsHADh9+jS2bdt2w/etWbOm4amImhipqpDJ\nSyFPZ0MZPxPCy0fvSEREmrK6gMTGxnKpciIbkqoK+cm7kOfO1q5w2tw5Z5YREVnD6gJiq+k3RARI\ntQZy+WLIokIor06H8GyudyQiIruw+iLU/fv3w8rLRojoOqTZDPnRfMiLRVDGTGX5IKImxeojIHPm\nzIGfnx/69u2Lfv368b4pRA0gzdVQl80DqiqhjJkC4W7QOxIRkV1ZPQtmz5492L59O3788UdUV1fj\ntttuQ//+/dG3b1+0bNlSq5wOi7NgtOWKV7XL6mqoS98ChIAy4n8g3N11zeOKY+yIOM7a4xhrz5az\nYKwuIFdcvnwZP/zwA3bs2IFDhw5BCIHw8HD0798fUVFR8PDwsFlIR8YCoi1X+4Uiqyqhvj8H8PCE\n8uIECDerD0LanKuNsaPiOGuPY6w9hyggv1ZUVIQdO3YgNTUVP//8Mzw8PNC7d28MGDAAd955py1y\nOiwWEG250i8UWVkJ9d1ZEMYWEM8nQDRrpnckAK41xo6M46w9jrH2dF0H5HpUVYXZbIbZbAZQe0fC\nn376CampqWjfvj1eeeUVtGvXzhYfReSUZMVlqEtmQrQKhBg+BkJxjPJBRKSXBheQ8vJy7Nq1C6mp\nqTh69CiaNWuGyMhIxMfH46677oIQAnv27MGKFSvw/vvvY86cObbMTeQ0ZPklqItnQAS3g/jTaAgH\nuY8SEZGerC4gGRkZSE1Nxb59+1BdXY3Q0FAMHz4cffv2hY/P1as39u7dGyaTCR999JHNAhM5E3mp\nDOqi6RDtO0EMG8HyQUT0f6wuIPPmzUPLli3x0EMPYcCAAXXeuK19+/bo169fgwMSOStZVgp1wVSI\nsDsh4p7nCsJERL9idQGZNGkS7rzzznr/Mu3UqRM6depkdTAiZyZLL0KdPwXizrshnvgzywcR0W9Y\nfTzYYDBg69atV/1s586dePXVV/HSSy9h+fLlUFXVZgGJnI28WAT17UkQkfeyfBAR3YDVBWTNmjU4\nfPiw5XleXh7effddKIqCkJAQ/Pvf/8a///1vm4Ykchay+EJt+YgaAOUP8SwfREQ3YHUBOXPmDDp3\n7mx5vn37dhgMBiQlJeH1119H//798d1339k0JJEzkBcKoM59HaLf76A8zJs2EhHdjNUF5PLly/D2\n9rY8P3DgACIiIuDlVXsL8a5du6KgoMB2CYmcgDx/DurcNyBiHoHy4BN6xyEicnhWFxB/f3+cPXsW\nQO0KqCdPnkRERITl9cuXL6OZg6zwSGQPMj8P6ttvQAyOhTLoMb3jEBE5BatnwfTu3RtbtmxBdXU1\nTpw4AXd3d/Tq1cvy+qlTpxAYGGjTkESOSv5yBur8qRB/iIfS93d6xyEichpWF5C4uDiUlJQgNTUV\nXl5eePnll9GiRQsAtaujpqenY/DgwTYPSuRoZO5pqAumQcQ+C+Xe+/WOQ0TkVGxyM7orVFVFRUUF\nDAYD3BzgLp/2wJvRactRby4lf86GujgR4ukXofRy7oX2HHWMXQ3HWXscY+3Z8mZ0Nl0XWlEUeHl5\nobCw0Ja7JXIo8uRxqAunQ4kf6fTlg4hILzY7TFFaWmq5Od2JEyewZs0aW+2ayGHI7CNQ302C8uwY\niO699Y5DROS0GlVAqqqqLDen++mnn1BTU4OgoCA8+uijtspH5DDksYNQP3gLyvOvQXTrqXccIiKn\nZnUBkVIiMzMTqamp2L17NyoqKgAAMTExePTRR216fojIUcjDB6AuexvKSxMgbu+udxwiIqdX7wKS\nk5OD7du344cffsDFixcRFBSEhx9+GJ06dcJbb72FHj16sHyQS5IHf4T68QIoI/8K0aWb3nGIiFxC\nvQpIQkIC8vLy0LJlS/Tt2xd9+vRBSEgIAODcuXOaBiTSkzywG+qKxVBGvwHR6Xa94xARuYx6FZC8\nvDwEBgYiPj4ed999N9zd3bXORaQ7+eMPUD99D8qYKRAdw/SOQ0TkUupVQEaOHIkdO3Zg0aJF8PDw\nQK9evdCnTx90785z4eSa1N07ID/7EMqr0yHah+odh4jI5dSrgMTExCAmJgZFRUVITU3Fjh07kJqa\nCh8fH9xxxx0AwNuOk8tQ076DXLsCSkIixK0d9I5DROSSGrwS6unTp5GamoqdO3eiqKgIvr6+iIyM\nxN13342IiAh4enraOqtD4kqo2rL3yobqzm8hv/y0tnwEt7Pb5+qJq0faB8dZexxj7dlyskmjl2KX\nUiIrKwvbt29HRkYGLl++DHd3d3z66ae2yujQWEC0Zc9fKOp/t0BuToGSMBMiqK1dPtMR8Je2fXCc\ntccx1p4tC0ijV0IVQqBbt27o1q0bXnrpJezevRs7duywRTYiu1H/sxHymy+hjE+CCGyjdxwiIpdX\nZwEpLy+Hl5dXvXbm7u6O6OhoREdHW70tkV7Ur9dBfv9vKBNnQ7QK1DsOEVGTUOfN6EaNGoXVq1ej\noKCg3js9f/48Vq1ahdGjRzcqHJHW1E0pkP/dCmUCywcRkT3VeQRk1KhRSElJwVdffYXQ0FBEREQg\nJCQEt9xyC7y9vQEAZWVlKCgoQE5ODjIzM5GdnY22bdti1KhRmn8BooaQUkJuSIbcs7P2yIdfS70j\nERE1KXUWkHvuuQe9e/fG3r178d1332HDhg0wm83Xfa+7uzt69OiB2NhY3HXXXZyaSw5JSgm57hPI\nzN1QJiRB+PrpHYmIqMmp10WoiqKgV69e6NWrF6qrq5GTk4Pc3FyUlZUBqL3yuG3btggJCYGbW6Ov\nayXSjJQS8vOPIY9k1l5wavTVOxIRUZNkdVtwd3dHly5d0KVLFy3yEGlGqirkZ8sgTx6rLR/ePnpH\nIiJqspz+cMXhw4fx1VdfITs7G6WlpQCAp556Ck8++WSd2+7cuRPr169Hbm4uDAYDunXrhvj4eAQF\nBWkdm+xMqirkqvchc0/XLjLm5a13JCKiJq3OWTCO7uTJkzhw4ACMRqNV223btg2LFy/GqVOn4O/v\nDykl0tPTMXXqVJSUlGiUlvQg1RrIFUsgz52F8tp0lg8iIgfg9AWkf//+WLFiBWbPnl3vbcxmM1av\nXg2g9iLbJUuWYP78+fD09ERJSQnWrVunVVyyM1lTA/mPhZBF56GMnQbhyXVpiIgcgdMXEB8fHxgM\nBqu2yc7OtizXGxUVBQDw9/dHWFjtLdf3799v25CkC2k2Qy57G/JSKZQxUyA8msb9iYiInIHTXwPS\nEBcuXLA89vX9/1kQLVq0AAAUFhbaPRPZlqyuhvrhXECtgfLyJAh360oqERFpS5MCYjKZrL4mwxHU\ndV++rKwsZGVlWZ7HxcU55fd0JgaDweoxllVVuPReEpoZDPAamwjh5q5ROtfQkDEm63Gctccxto+U\nlBTL4/DwcISHhzdoPzYrIBUVFfjxxx+RlpaGAwcOYMWKFbbatc21atXK8vjKzJlfPw4ICLjudtcb\naN55UVvW3t1SVlZCfS8JwtsI8VwCyi5XAKjQLqAL4B1E7YPjrD2OsfaMRiPi4uJssq9GFZDy8nLs\n2bMH6enpyMzMhKqq6N69OxISEmwSzlYSExNRXFyM3r17Y9iwYQgNDYWPjw/KysqQlpaG6OhoFBUV\n4dixYwCAHj166JyYGkJWVkBdMhPCPwBi+FiIZs30jkRERDdgdQEpKyvD7t27kZaWhoMHD6JZs2bo\n0aMH/vKXv6Bnz57w9LTvhX4ZGRn45JNPrvrZ5s2bsX37dnTu3BljxoxBfn4+CgsLUVxcDABwc3PD\nsGHDsGzZMqSnp2PMmDEwmUyoqKiAr68vhgwZYtfvQI0nL5dDXZwIcUswxJ9fhlBYPoiIHFm9CkhV\nVRX++9//Ii0tDYcOHYLBYEDPnj0xZswYZGRkYPDgwZYZJPZWXl5+zZ16L126hEuXLl1zKuXX96YZ\nNGgQPD09sWHDBuTm5sLd3R1RUVGIj4+Hnx/vDeJMZHkZ1EUzIG7rCBH/FwjF6Sd3ERG5PCHruvIS\nwN///nccOXIEvXr1QlRUFCIiIiz3fKmursaiRYtw55134sEHH9Q8sKPJy8vTO4JLq+ucrrxkgrpg\nGkSn2yGefpE3QGwAnje3D46z9jjG2gsODrbZvup1BMTLywtLly6Fu/u1swnc3d0xbtw4vPPOOzh+\n/DhGjBhh9bocRA0hTSVQ50+FCO8BETuc5YOIyInU61j16NGjr1s+LDtRFIwdOxbNmzfHpEmTrjkl\nQmRrsrQY6tuTICJ6sXwQETmhehUQpZ7n1F944QX07NkTiYmJjQpFdDPy4gWocydB9OoL5fE/sXwQ\nETkhmy9ENnToUM4iIc3IovNQ502G6PsAlN/H6h2HiIgaSJOVUO09FZeaBlmYX1s+Yh6B8rs/6B2H\niIgaoUneC4acjyzIgzpvCsTgJ6Dc/7DecYiIqJFYQMjhyV/OQl0wFeKRp6H0b3pTvYmIXBELCDm0\nmjMnoc6fDPH4M1CiB+odh4iIbIQFhByWPHMSZYtnQDz5HJSoAXrHISIiG2IBIYckT5+AujgRXs+/\nisrwnnrHISIiG+NNM8jhyJyjUBfNgPLMaBjuuU/vOEREpAEeASGHIo8fgvr+HCjPvQpx5916xyEi\nIo2wgJDDkEd/grr071BeHAdxR6TecYiISEMsIOQQ5KF9UD+aD2XERIiuEXrHISIijbGAkO7kT3ug\n/nMRlFGvQ3S+Q+84RERkBywgpCu5Pw3qynehvDwJIrSr3nGIiMhOWEBIN3LvTqirl0J5dRpE+056\nxyEiIjtiASFdqOn/hfz8YyivzYC4raPecYiIyM5YQMju1F3/gVz3CZSEmRBt2+kdh4iIdMACQnal\npn4NueEzKONmQbS5Ve84RESkExYQshv1u02QW76AMmEWRGCw3nGIiEhHLCBkF+o3X0Fu2whlQhJE\n6yC94xARkc5YQEhz6r//BbnjaygTZkO0aq13HCIicgAsIKQpdeNnkOn/hTJxNoRfK73jEBGRg2AB\nIU1IKSG/WgW5L622fPj66x2JiIgcCAsI2ZyUEvJfyyGz9tde82FsoXckIiJyMCwgZFNSSsg1H0Ge\nOFw728XbqHckIiJyQCwgZDNSVSGTl0KezoYyLhHCy0fvSERE5KBYQMgmpFoD+cl7kOdyoSQkQjT3\n0jsSERE5MBYQajRZUwO5fDFkcWHtjeU8m+sdiYiIHBwLCDWKNJshP14AeckEZcxUCA8PvSMREZET\nYAGhBpPmaqjL3gaqq6G8MhnC3aB3JCIichIsINQgsroa6tK3ACGgjHodwt1d70hEROREWEDIarKq\nEur7cyA8mkO8OB7CjX+MiIjIOvybg6wiKyugvpsEYfSDeP41iGbN9I5EREROiAWE6k1WlENdMhMi\nIAji2VcgFJYPIiJqGBYQqhdZfgnq4hkQwe0g/jQaQlH0jkRERE6MBYTqJC+VQV04DaJjZ4ihI1g+\niIio0VyigOzcuRPr169Hbm4uDAYDunXrhvj4eAQFBd1wm88//xxr16697mvJyclQ+JcsAECaSqEu\nnArR5U6Ip56HEELvSERE5AKcvoBs27YNS5cuBQAEBgairKwM6enpOHLkCObOnYsWLW5+J1ZfX1/c\ncsstludCCP4l+39k6UWo86dARNwN8fifOS5ERGQzTl1AzGYzVq9eDQC45557kJCQgOLiYrz22mso\nKSnBunXrMHz48JvuIzIyEqNHj7ZDWuciLxbVlo+7+0A8Oozlg4iIbMqpC0h2djZMJhMAICoqCgDg\n7++PsLAwZGZmYv/+/XXuIz09Hbt27YK3tzc6duyIoUOHokOHDlrGdniyqLC2fNx7P5SH4/SOQ0RE\nLlwYTlkAABLPSURBVMipL3S4cOGC5bGvr6/l8ZXTLoWFhTfdXlEU+Pn5ITAwEBcvXsS+ffswefJk\nnDp1SpO8zkBeKID69hsQ/R5g+SAiIs049RGQG5FS1vmefv364aGHHoK3tzcAIDMzE0lJSaiursbW\nrVsxcuRIrWM6HHn+HNR5kyF+9wcoAx/VOw4REbkwpy4grVq1sjwuLS295nFAQMANt/3tDJmIiAj4\n+PigrKzshkdOsrKykJWVZXkeFxcHo9HYoOyOpibvDMrmTUbzx/8Ij0GP6R3HwmAwuMwYOyqOsX1w\nnLXHMbaPlJQUy+Pw8HCEh4c3aD9OXUBCQ0MtpSEtLQ3R0dEoKirCsWPHAAA9evQAACQmJqK4uBi9\ne/fGsGHDAACbNm1CdHQ0/P39AdQeASkrKwMAtG7d+rqfd72BvnINijOTv5yBOn8qxB/iURV1P6oc\n6DsZjUaXGGNHxjG2D46z9jjG2jMajYiLs83peacuIG5ubhg2bBiWLVuG9PR0jBkzBiaTCRUVFfD1\n9cWQIUMAAPn5+SgsLERxcfH/tnf30VHV+R3HP/c2T0IyhCSEBEF5aEANReAcknAW5IBUW2gtEqQh\neg5oj7WHGo/4sLUuKkSxWxHWQusDrKCwh2i0EggC6pFtgbhJLAukG57cKIqJkg1gEmADSe7tH2nm\nBCEwmczcm8m8X3/defjdfOfL9yQf7r0z4127fft2bdiwQUlJSYqOjlZ1dbUkKSYmRjNnznTl9bjB\n/va4rFeWyJgzX2bWVLfLAQCEiZAOIJI0ffp0xcTEqLi4WNXV1YqMjFRmZqZyc3MVHx9/yXM7vpV0\n9uzZKi0t1YkTJ1RbW6vk5GSNGjVK2dnZSk1NdfpluML+pkrWvy2VkfOgzAmT3S4HABBGDNuXKzbR\nqZqaGrdL8Iv91ReyVufLvG+hjPET3S6nUxxSDT567Az6HHz0OPgGDRoUsH2F/BEQdJ39+8OyXn1R\n5vxHZNw6we1yAABhiAASZuxjv5P1+r/KfGCRjNHj3S4HABCmCCBhxD58UNbal2U++ISMm291uxwA\nQBgjgIQJ+3e/lbXuFzL/4Z9kjBztdjkAgDBHAAkD9sFyWW+vlvmPP5Mx4ia3ywEAgADS29m//UzW\nr16TmfesjGFpbpcDAIAkAkivZn2+R/a7v5T56BIZN4xwuxwAALwIIL2U9Ztfy/7Pt2U+ulTG4KFu\nlwMAwCUIIL2QtfcT2Vs2yXz8eRmpQ9wuBwCAyxBAehnrv3bI3vGezMdfkJFyvdvlAABwRQSQXsT6\ntFj2J1tkPvGijAEpbpcDAECnCCC9hPXRZtn/vUPmky/KSEx2uxwAAK6KANILWB8Wyv7Nr9uOfCQk\nuV0OAADXRAAJYbZty95aIHtficwnlsmIT3C7JAAAfEIACVG2bcvevEF2xf+0hQ9PvNslAQDgMwJI\nCLJtW3bhOtnH/rctfMR63C4JAIAuIYCEGNuyZL+zRvbx38t87AUZfWPdLgkAgC4jgIQQ27Jk/+pV\n2TXftH3CaZ++bpcEAIBfCCAhwrZaZb+1Wvapk23f7RLTx+2SAADwGwEkBNitrbLXvSK78QeZjzwn\nIzrG7ZIAAOgWAkgPZ7e0yP7lCtkX/ijz4cUyoqLdLgkAgG4jgPRgdnOzrDUvSbYtc+HPZERGul0S\nAAABQQDpoezmi7Je+7kUGSXzwcdlRBA+AAC9BwGkB7IvXJD1Hy/IiPXIeGCRjAj+mQAAvQt/2XoY\nu+mPslY/LyNxgIwFj8gw/8TtkgAACDgCSA9inz8na9VSGYNukHHfQhmm6XZJAAAEBQGkh7DPnZX1\nynMyhqXJyPl7wgcAoFcjgPQAdmODrF88I+OmMTLueUCGYbhdEgAAQUUAcZndcEbWymdl3JohY9Z9\nhA8AQFgggLjI/uGUrBXPyMi4TcZf/S3hAwAQNgggLrFP/UHWysUyJt0h8y+z3S4HAABHEUBcYP/h\ne1krFsu4/a9l/vnfuF0OAACOI4A4zD5ZI2vlMzL+Ilvm1BlulwMAgCsIIA6yvzvRdsHpXfNkTr7D\n7XIAAHANAcQh9rdfyXplqYzs+TInTnW7HAAAXEUAcYD9dVXbJ5zmPChzwmS3ywEAwHUEkCCzvzwq\n699fkHnfQhnjJ7pdDgAAPQIBJIjsLw7Jeu1fZC54RMaYCW6XAwBAj0EACRL7SIWsNctl/t1jMtLH\nuV0OAAA9Sq8IICUlJdq6dauqq6sVFRWl0aNHKzc3VykpKUFZdy125X5Zb66U+dBPZYz6s27tCwCA\n3ijkv3J1165dWrVqlY4fP67+/fvLtm2VlZXp2WefVX19fcDXXYt98PO28LHwnwkfAAB0IqQDSEtL\nizZt2iRJysrK0urVq7Vy5UrFxMSovr5emzdvDui6a7E+3yPr7VUy856R8ae3+PeiAAAIAyEdQKqq\nqtTY2ChJyszMlCT1799fI0eOlCQdOHAgoOuuxtrzsezCN2U+9ryMYSO7vB4AgHAS0teAnDp1yrvt\n8Xi82/369ZMk1dXVBXRdZ6xPtsj+tFjmEy/KGDioS2sBAAhHIR1AOmPbdlDWVVZWqrKy0nt77ty5\niiguUHP5HsUtXSUzaaBfPxedi4qKUlxcnNtl9Gr02Bn0OfjosTMKCwu92+np6UpPT/drPyEdQBIT\nE73bDQ0Nl20nJSUFdN2VGn3h0EGZP/25zkX3kf7/tA4CJy4uznu6DMFBj51Bn4OPHgdfXFyc5s6d\nG5B9hfQ1ICNGjFBsbKwkqbS0VJJ0+vRpHTt2TJI0duxYSVJ+fr4WLVqkgoKCLq3zhfn48zLiPNd+\nIgAA8ArpIyARERGaN2+e1q5dq7KyMuXl5amxsVFNTU3yeDyaNWuWJOnkyZOqq6vTmTNnurTOF0ZE\nZFBeGwAAvVlIBxBJmj59umJiYlRcXKzq6mpFRkYqMzNTubm5io+Pv+S5hmH4tQ4AAASWYft7xSYk\nSTU1NW6X0KtxTjf46LEz6HPw0ePgGzQocO/0DOlrQAAAQGgigAAAAMcRQAAAgOMIIAAAwHEEEAAA\n4DgCCAAAcBwBBAAAOI4AAgAAHEcAAQAAjiOAAAAAxxFAAACA4wggAADAcQQQAADgOAIIAABwHAEE\nAAA4jgACAAAcRwABAACOI4AAAADHEUAAAIDjCCAAAMBxBBAAAOA4AggAAHAcAQQAADiOAAIAABxH\nAAEAAI4jgAAAAMcRQAAAgOMIIAAAwHEEEAAA4DgCCAAAcBwBBAAAOI4AAgAAHEcAAQAAjiOAAAAA\nxxFAAACA4wggAADAcQQQAADgOAIIAABwHAEEAAA4LsLtArqjtbVVH3zwgXbv3q1Tp06pX79+yszM\nVE5OjmJiYq66dsmSJTp8+PBl9990001aunRpsEoGAAAK8QDy6quvau/evTJNUykpKaqtrdWOHTv0\n9ddf67nnnrvqWsMwJEkDBw6Ux+Px3j9kyJCg1gwAAEI4gHz11Vfau3evJOn+++/XHXfcoX379uml\nl17SoUOHVF5eroyMjGvuJzs7W1OmTAl2uQAAoIOQvQZk//793u32oDF+/HhFRkZKkg4cOODTft56\n6y3de++9ysvL05o1a1RfXx/4YgEAwCVCNoCcOnXKu92vXz9JbadV4uLiLnu8M1FRUUpMTJTH41Ft\nba0+/fRTLV68WBcvXgxO0QAAQFIPPAXzzjvvaPPmzVd9zrWu7/DF/PnzNXjwYEVEtLWgoKBARUVF\nqq2tVXl5uSZNmnTZmsrKSlVWVnpvz507V4MGDep2Lbi69lCJ4KHHzqDPwUePg6+wsNC7nZ6ervT0\ndL/20+MCyPDhw696TYZhGIqPj1diYqL3vvr6esXHx8u2bTU2NkrSJY9fydChQy+5PWnSJBUVFUmS\n6urqrrjmx40uLCzU3Llzr/pz0D30OPjosTPoc/DR4+ALZI97XADJyMjw6eLRsWPH6t1335UklZWV\n6c4779S+ffvU3NzsfVySTp8+rfz8fBmGodzcXE2YMEENDQ0qLS3VlClTFB0dLUkqKSnx7nvAgAGB\nflkAAKCDHhdAfDV8+HD95Cc/UUlJidavX6+dO3fq5MmTkqRbbrnFG2JaW1v13XffSZLOnz8vSbpw\n4YLefPNNrV+/XikpKbpw4YL3mpHBgwcrMzPThVcEAED4CNkAIkkPP/ywUlNTtXv3btXW1srj8Sgr\nK0s5OTmXPbf9cz8kyePxaPbs2aqoqND333+vixcv6vrrr1dGRobuuusu73Uh1+LveS/4jh4HHz12\nBn0OPnocfIHssWHbth2wvQEAAPggZN+GCwAAQhcBBAAAOI4AAgAAHEcAAQAAjgvpd8EES0lJibZu\n3arq6mpFRUVp9OjRys3NVUpKSlDWhSN/evXee+/p/fffv+JjBQUFMk3ydLvDhw9ry5YtqqqqUkND\ngyTpnnvu0Zw5c665ljn2nb99ZpZ9V1xcrP3796umpkaNjY3yeDxKS0vTnDlzdMMNN1x1LbPsG397\n3N05JoD8yK5du/TGG29IkpKTk3X27FmVlZXpyJEjWr58ufd7ZwK1Lhx1t1cej0cDBw703jYM45K3\nWaPt26IPHjyo1NRU7x9GXzDHXeNvn9sxy9e2c+dO1dXVKSUlRdddd51qampUVlamgwcPasWKFUpK\nSrriOmbZd/72uJ2/c0wA6aClpUWbNm2SJGVlZWnRokU6c+aMHn30UdXX12vz5s1asGBBwNaFo0D0\naty4cVq4cKED1Yau2267TdOnT5dlWZo/f75Pa5jjrvOnzx0xy9c2bdo0TZ48WcnJyZKkbdu2aePG\njWpqalJ5eblmzJhx2RpmuWv86XFH/s4xAaSDqqoq73fJtH8aav/+/TVy5EhVVFTowIEDAV0XjgLR\nq7KyMn322Wfq27evhg0bppycnMu+2yfcxcbGSpKampp8XsMcd50/fe6IWb627OzsS26PGTPGu93Z\nh0Yyy13jT4878neOOdHYQfvHsUtth5TatR+q6+xL6vxdF4662yvTNBUfH6/k5GT98MMP2r9/vxYv\nXqzjx48Hpd5wwhw7i1n2z7Zt2yS1zejEiROv+BxmuXt86XG77swxR0B84O+HxfIhs77zpVeTJ0/W\njBkz1LdvX0lSRUWFli1bpubmZn300Ud66KGHgl1mWGKOA49Z7rqWlha9/vrr2rNnj/r06aMnn3xS\ncXFxXdoHs3x1Xe1xd+eYANJBYmKid7vjBWXt251diOPvunDUnV79+Mr1MWPGKDY2VmfPnuV/NAHA\nHDuHWe6axsZGLV++XEePHlVCQoKeeuop3XjjjZ0+n1nuuq72WOr+HHMKpoMRI0Z4z+mWlpZKkk6f\nPq1jx45JksaOHStJys/P16JFi1RQUNCldfC/x5L04Ycf6syZM97bFRUVOnv2rCRpwIABjtTfmzDH\nzmCWu+fbb7/V008/raNHj2ro0KFatmzZZX8YmeXu8afHUvfnmCMgHURERGjevHlau3atysrKlJeX\np8bGRjU1Ncnj8WjWrFmSpJMnT6qurs7beF/Xwf8eS9L27du1YcMGJSUlKTo6WtXV1ZKkmJgYzZw5\n05XX01OVl5dr48aNl9y3fft27d69W2lpacrLy2OOA8CfPrc/h1n2zcsvv6za2lpJUmtrq1asWOF9\n7Pbbb9e0adOY5W7yp8dS9+eYAPIj06dPV0xMjIqLi1VdXa3IyEhlZmYqNzdX8fHxlzy34/ucu7Iu\n3Pnb49mzZ6u0tFQnTpxQbW2tkpOTNWrUKGVnZys1NdXpl9GjnT9/3vsLpd25c+d07ty5yw4/M8f+\n87fPzLLvWlpavNsnTpy45LFx48ZdcptZ9o+/Pe7uHBs2V+UAAACHcQ0IAABwHAEEAAA4jgACAAAc\nRwABAACOI4AAAADHEUAAAIDjCCAAAMBxBBAAAOA4AggAAHAcAQQAADiOAAIAABzHl9EBCBmnT59W\nYWGhBg4cqObmZo0ZM0aHDh3S7Nmz3S4NQBdxBARASGhoaNCSJUs0adIk3X333crIyFB+fj7fHguE\nKAIIgJCwYcMGDRkyRKNHj5YkxcbGqrW1VSNHjnS5MgD+IIAA6PEaGxtVUlKirKws731HjhxRQkKC\nEhMTXawMgL8IIAB6vC+++EKWZenmm2/23nfkyBGlpaW5WBWA7iCAAOjxmpubFRUVpaSkJO99R48e\n5fQLEMIIIAB6vLS0NJmmqZaWFknSxx9/rG+++YYAAoQw3oYLoMdLSEjQggULtG7dOiUkJOj8+fOK\niIjQ8OHD3S4NgJ8IIABCwtSpUzV16lRJUlFRkYYOHaqICH6FAaGKUzAAQs6XX37J6RcgxBFAAISc\nqqoq3gEDhDiOXwIIGQcOHFBRUZHq6uq0Y8cO9enTR2PHjnW7LAB+MGzbtt0uAgAAhBdOwQAAAMcR\nQAAAgOMIIAAAwHEEEAAA4DgCCAAAcBwBBAAAOI4AAgAAHPd/IfQLtRnad7oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x19cb2c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def lambda_of_params(beta, mu, epsilon, Ldinvsq, kf, q):\n",
    "    \"\"\"Given parameters, calculate the eigenvalue lambda for the asymptotic WKE.\"\"\"\n",
    "    \n",
    "    funres = lambda lamb: disp_relation_residual(beta, mu, epsilon, Ldinvsq, kf, q, lamb)\n",
    "    domainmin = -2*mu\n",
    "    domainmax = 100\n",
    "    lamb = scipy.optimize.brentq(funres, domainmin, domainmax) #1D root finder\n",
    "    return lamb\n",
    "\n",
    "def disp_relation_residual(beta, mu, epsilon, Ldinvsq, kf, q, lamb):\n",
    "    \"\"\"Compute the residual for the dispersion relation.  lamb must be real.\"\"\"\n",
    "    mf = Ldinvsq / (kf*kf)\n",
    "    lhs = lamb + mu\n",
    "    rhs = q*q*epsilon/mu * kf**4 * (lamb + 2*mu) * (1+mf) * polarint(beta, mu, epsilon, Ldinvsq, kf, q, lamb)\n",
    "    res = rhs - lhs\n",
    "    return res\n",
    "\n",
    "def polarint(beta, mu, epsilon, Ldinvsq, kf, q, lamb):\n",
    "    \"\"\"Carry out the polar integral\"\"\"\n",
    "    fun = lambda phi: polarint_integrand(beta, mu, epsilon, Ldinvsq, kf, q, lamb, phi)\n",
    "    out, err = scipy.integrate.quad(fun, 0, 2*np.pi, epsabs=2e-6, epsrel=2e-6)\n",
    "    out = out / (2*np.pi)\n",
    "    return out\n",
    "\n",
    "def polarint_integrand(beta, mu, epsilon, Ldinvsq, kf, q, lamb, phi):\n",
    "    \"\"\" As explained above, integrand is A/(B+iC)^2.  Keep only the real part.\n",
    "    So use A(B^2 - C^2) / (B^2 + C^2)^2.\"\"\"\n",
    "    c = np.cos(phi)\n",
    "    s = np.sin(phi)\n",
    "    mf = Ldinvsq / (kf*kf)\n",
    "    A = (1 + mf - 4*c*c) * s * s\n",
    "    B = (lamb + 2*mu) * kf**2 * (1+mf)**2\n",
    "    C = -2*beta*q*c*s\n",
    "    \n",
    "    out = A*(B*B - C*C)  / (B*B + C*C)**2\n",
    "    return out\n",
    "\n",
    "# Now, let's go and calculate the dispersion relation lambda(q)\n",
    "beta = 1\n",
    "mu = 0.02\n",
    "epsilon = 1\n",
    "Ldinvsq = 1\n",
    "kf = 1\n",
    "\n",
    "q1 = np.logspace(-4, 0, 150)\n",
    "q2 = np.linspace(1.02, 2.2, 150)\n",
    "q = np.concatenate([q1, q2])\n",
    "lamb = np.zeros_like(q)\n",
    "\n",
    "for j in range(len(q)):\n",
    "    lamb[j] = lambda_of_params(beta, mu, epsilon, Ldinvsq, kf, q[j])\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.plot(q, lamb)\n",
    "plt.xlabel(r'$q$')\n",
    "plt.ylabel(r'$\\lambda$ (Asymptotic WKE)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CE2-GO\n",
    "The only difference for the general dispersion relation in CE2-GO as compared to the asymptotic WKE is the presence of a $(1 - q^2 / \\ol{k}^2)$ term.  For the same isotropic thin-ring forcing, the dispersion relation reduces to the same result except that in CE2-GO there is an additional factor of \n",
    "  \\begin{equation}\n",
    "    1 - \\frac{q^2}{k_f^2 (1+m_f)}\n",
    "  \\end{equation}\n",
    "in front of the polar integral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x19e1e518>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh4AAAIACAYAAADNH2WhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VOWhN/Dfc5KZTJaZ7CGEnUhAElmFsG8Gsa0WqpAr\n3Ndqb7erFisuVV+Lrdxbe62KVfStlF6r5VNQqqWACwqENZCAyCJhDyAkgYQsJJOEyTbP+0ckhZKQ\nSTLnPHNmft/Pxw+T5JyT3zxNw49znvMcIaWUICIiIjKApjoAERERBQ4WDyIiIjIMiwcREREZhsWD\niIiIDMPiQURERIZh8SAiIiLDsHgQERGRYYJVB+isI0eOYM2aNcjPz0dVVRUAYM6cOZg9e7bHx1i8\neDFyc3MBAOnp6Xjsscd0yUpERETNTHvG4/Tp0zhw4ADsdnun9t+8eXNL6eiovLy8Tu1HnuMY649j\nbAyOs/44xvrz5hibtnhMmjQJ7777Ll544YUO73vhwgW88847SElJQUxMTIf35w+5/jjG+uMYG4Pj\nrD+Osf5YPABERETAarV2eD+3240lS5ZA0zQ88sgj0DTTDgEREZHpBNzfuqtWrcLJkyfxox/9CPHx\n8arjEBERBRTTTi7tjFOnTmHNmjWYOHEixo8f7/F+eXl515xmyszM1CMeXYVjrD+OsTE4zvrjGOsv\nMzMTq1atavk4NTUVqampnTqWMPvTaV0uF+6//34A7d/VsmXLFvzhD3+AxWJpucRSV1cHANA0DRaL\nBUuXLkVoaGi737eoqMgL6aktdrsdTqdTdQy/xjE2BsdZfxxj/SUlJXntWH59xmPRokWoqKjA6NGj\nMXfu3JbPNzQ0XLet2+1GfX09TN7DiIiIfJppi8fu3buxfPnyaz73ySefYNu2bRgwYADmz5+P4uJi\nlJaWoqKiAgAwZcoUTJky5Zp9Hn74YZSWlnIdDyIiIgOYtnjU1taipKTkms/V1NSgpqYGcXFx13xe\nCNHu8TzZhoiIiLrG9HM8VOEcD33xmq3+OMbG4Djrj2OsP2/O8Qi422mJiIhIHRYPIiIiMgyLBxER\nERmGxYOIiIgMw+JBREREhmHxICIiIsOweBAREZFhWDyIiIjIMCweREREZBgWDyIiIjIMiwcREREZ\nhsWDiIiIDMPiQURERIZh8SAiIiLDsHgQERGRYVg8iIiIyDAsHkRERGQYFg8iIiIyDIsHERERGYbF\ng4iIiAzD4kFERESGYfEgIiIiw7B4EBERkWFYPIiIiMgwLB5ERERkGBYPIiIiMgyLBxERERmGxYOI\niIgMw+JBREREhmHxICIiIsOweBAREZFhWDyIiIjIMCweREREZBgWDyIiIjIMiwcREREZhsWDiIiI\nDMPiQURERIZh8SAiIiLDsHgQERGRYVg8iIiIyDAsHkRERGQYFg8iIiIyDIsHERERGYbFg4iIiAzD\n4kFERESGCVYdoLOOHDmCNWvWID8/H1VVVQCAOXPmYPbs2Tfcb926ddi3bx+KiorgdDrhcDgwYMAA\nzJ49G7179zYiOhERUcAy7RmP06dP48CBA7Db7R3ab/369cjLy0NISAgSEhJQXl6O3NxcLFy4EKWl\npTqlJSIiIsDEZzwmTZqEjIwMuN1u3H///R7vN23aNEycOBEJCQkAgI8++gjLly+Hy+XC7t278e1v\nf1uvyERERAHPtMUjIiICAOByuTq03z333HPNx0OGDGl5HRxs2uEgIiIyhYD/m/ajjz4CADgcDowd\nO1ZxGiJ9SCmBynKg6Cxk0Vmg6ByqL9fAHRULJHSHiO8O9EuBCI9QHZWI/FzAFo/Gxka89dZb2L59\nO8LCwvDkk0+2OV8kLy8PeXl5LR9nZmZ2eG4JdYzVauUYd4KUEvJSOZoKzsB97jSaCs6gqeAMZOHX\nQFAwgnr2hdazL4IGDIYlJg715wvgLi6E+8BuNJ46hqDeybCMGAPrxOnQYuJVvx2/wJ9l/XGMjbFq\n1aqW16mpqUhNTe3UcQKyeDidTrz00ks4duwYYmJi8PTTT6NPnz5tbt/aADudTr1jBjS73c4x7gDZ\n2Ai5/TPIj/8GNDUASb0hkno3/zliHERSbwh7JCSApm/+C7Hb0XDVGGv1dZDHD6FuXy5cT/wHROpw\niOkzIfqlqHpbfoE/y/rjGOvPbrcjMzPTK8fy6+KxaNEiVFRUYPTo0Zg7dy4AoKCgAC+++CJKSkrQ\nt29fPPXUU4iJiVGclKhzpJTAV1/A/bc/A1Ex0OYvBHr3hxCiw8cS1hAgbSRE2kjIe74Pmb0J7rf+\nB6LfQIg5P4CITdDhHRBRoDFt8di9ezeWL19+zec++eQTbNu2DQMGDMD8+fNRXFyM0tJSVFRUtGzz\n8ssvo6SkBADQ1NSEV155peVrt912G6ZNm2bMGyDqInnuNNx/exuoKIM2+wfAkFs7VThaI8IiIKbP\nhJx0B+Rnf4f7vxdA3HYXxLfmQAQFeeV7EFFgMm3xqK2tbSkQV9TU1KCmpgZxcXHXfP7qX8aNjY0t\nr8+dO3fNdsOHD9chKZF3yUtlkP/4K+TBPRB33QsxcQaETndkiZAQiO/OhRyfAfdflkDmPQvtx09A\nxMS1vzMRUSuElFKqDmFGRUVFqiP4NV6zvZ6sc0F+/g/ITesgJmRAfHsORFjn70Lp6BhLtxty/YeQ\nm9ZB+8HPIdJGdvp7BxL+LOuPY6y/pKQkrx3LtGc8iAKFlBIyZwvk6uUQN90M7dlXIOITDc8hNA3i\n23MgbxoM9x9/B/G9+6CNzzA8BxGZG4sHkQ+TTU2Q7/0R8sRhaD/9BUTyINWRIFJSoT3xG7h//2u4\nqyoh7rjba3NLiMj/mfZZLUT+Tta54P7DbyGLi6A99aJPlI4rRGJPaE+9CJm7BXLNX1XHISITYfEg\n8kGy6hLcLz8LERYB7ZHnIELDVEe6joiOhfb4byD3ZsO9YY3qOERkEiweRD5GXiiE+39+AZE2EuIH\nP4cItqiO1CZhd0B7dBHkxjVw78xSHYeITIDFg8iHyJNH4H7pGYhvzYY2c54p5k6I2Hhojz4P+eE7\nkHn7VMchIh/H4kHkI+TenXD/vxeg/eDn0CberjpOh4juvaD95Bdw/+9iyNJi1XGIyIexeBD5APfG\nNXC/twzao7827foYYmAaxB33wP2H/4FsqFcdh4h8FIsHkULS7Yb7/T9Bbvsc2tMvQvROVh2pS8T0\nmRAJ3SH/+ha4NiERtYbFg0gRKSXk8jchz+Y33y7rBw9hE0JA3D8fMv8o5O5tquMQkQ9i8SBSRH70\nPuS509DmPwcR3vmlz32NsIVC+48FkO//CbKyov0diCigsHgQKeDeuQkyeyO0+QshbKGq43id6DcA\nYuIMuJe/yUsuRHQNFg8ig8nD+yA/eAfaz38FERmtOo5uxF3/BpSVQO7arDoKEfkQFg8iA8mC03D/\naTG0/3waonsv1XF0JYItzZdcPvgzZNUl1XGIyEeweBAZRJaXwv36f0HM/QlESqrqOIYQvfpBjJkC\nuXq56ihE5CNYPIgMIGtr4H79eYjb7oQ2aqLqOIYSd82F/Gov5OkTqqMQkQ9g8SDSmWxsgPut/4FI\nSYW4/Xuq4xhOhIZBfO8+uFcuhXS7VcchIsVYPIh0JKWE/MsbgDUE4t4fm+LZK3oQY6cCAGQOJ5oS\nBToWDyIdybUrIC8UQvvxkxBakOo4yghNgzb3J5Crl0PW16mOQ0QKsXgQ6cSdswUydyu0n/0SIiRE\ndRzlRL8UoP9AyC2fqo5CRAqxeBDpQF4ohHz/T9AeegbCEaU6js/QvjsP8rO/Q7ouq45CRIqweBB5\nmWyoh/uPv4OYOQ+iZz/VcXyK6NEHYtBQyE3rVEchIkVYPIi8TH7wDhDfHWLyt1RH8UnirnshN66F\nrK1WHYWIFGDxIPIiuS8H8sBuaN//WcDewdIekdgDYuhoyA1rVEchIgVYPIi8RJZdhHv5m9B+/IRf\nPW1WD+LbcyC3fMq5HkQBiMWDyAtkUxPcf3oZYvosiORBquP4PJHQHUhJg8zeqDoKERmMxYPIC+Ta\nlc2LhM0IvJVJO0ub8T3IDWsgm5pURyEiA7F4EHWRPHIAMnsjtB8ugND4fylPif4Dgeg4yL3ZqqMQ\nkYH4W5KoC2RVBdxvvwrtPx6FcESrjmM62ozvQX62GlJK1VGIyCAsHkSdJN1uuN/+PcTYaRCDh6mO\nY05DRgH1LuDYV6qTEJFBWDyIOklu+Afgugzx3Xmqo5iW0DSI2+6Ce/PHqqMQkUFYPIg6QZ4vgFz/\nIbQfPQ4RHKw6jqmJ9CnA0YOQl8pVRyEiA7B4EHWQdLvhfvd1iLvmQsR1Ux3H9ERoGMTI8by1lihA\nsHgQdZDc8gkgBMSUb6uO4jfE5Dsgt30G6eattUT+jsWDqANkWQnkupXQvj+ft856kehzE2CPBPL2\nqY5CRDrjb04iD0kp4V7+JkTGTIjuPVXH8Tti8h1wb12vOgYR6YzFg8hDctdmoPISxIy7VUfxS2LU\nRODEYcjyi6qjEJGOWDyIPCCrKiA/+DO0Bx7hXSw6EbZQiFvHQ+ZuVR2FiHTE4kHkAbnijxDjMyD6\nJKuO4tfEmKmQuzZzJVMiP8biQdQOuS8H8txpiLvuVR3F/yUPAurrgHOnVSchIp2weBDdgKythnvF\nUmj3/wzCGqI6jt8TmgaRPgUyd4vqKESkExYPohuQH7wDMWw0REqa6igBQ4yZDJm7jWt6EPkpFg+i\nNsgjByDzvoS4+37VUQKK6N4LiIoBjvLBcUT+iMWDqBWysQHuv74Fbe5PIULDVMcJOGLMZMicLapj\nEJEOWDyIWiE3rQMSukMMS1cdJSCJUZMgD+RC1tWpjkJEXmbaBQmOHDmCNWvWID8/H1VVVQCAOXPm\nYPbs2e3um52djbVr16KwsBBWqxVpaWmYN28eEhMT9Y5NJiAvlTc/efap36mOErBEZDTQOxk4vA8Y\nPkZ1HCLyItOe8Th9+jQOHDgAu93eof2ysrLw+uuv48yZM4iOjoaUErm5uXjuuedQWVmpU1oyE/nh\nuxATbodI7KE6SkATI8ZC7tulOgYReZlpi8ekSZPw7rvv4oUXXvB4n8bGRqxYsQIAMGbMGCxZsgSL\nFy+GzWZDZWUlVq9erVdcMgl58gjk0YMQ35mjOkrAE8PGQB78ArKxUXUUIvIi0xaPiIgIWK3WDu2T\nn58Pp9MJAEhPb752Hx0djZSUFADA/v37vRuSTEW6m+Be+UeIe+6HsHFCqWoiOhZI6A4c590tRP7E\ntMWjM8rKylpeOxyOlteRkZEAgNLSUsMzke+QOzYAVitE+mTVUegbYsRYyC95uYXIn5h2cqk3tfdc\niLy8POTl5bV8nJmZ2eG5JdQxVqvV0DF2VzvhXLsSEU+/iOCrSqk/M3qMO6NpQgaqn/85In76BIQW\npDpOp5hhnM2OY2yMVatWtbxOTU1Fampqp44TUMUjNja25fWVO2Gufh0XF9fqfq0N8JVLNqQPu91u\n6Bi7VywFhqXjclwiECD/2xo9xp0SEQkZbofzwBcQNw1WnaZTTDHOJscx1p/dbkdmZqZXjuXXl1oW\nLVqEBQsWYOXKlQCA5ORkREREAABycnIAAOXl5Th+/DgAYNiwYWqCklKy4DTkFzsgZv0f1VGoFWI4\nL7cQ+RPTnvHYvXs3li9ffs3nPvnkE2zbtg0DBgzA/PnzUVxcjNLSUlRUVAAAgoODMXfuXCxbtgy5\nubmYP38+nE4nXC4XHA4HZs2apeKtkEJSSrhXLoO4ay5ERGBcYjEbMWIs3G/+BnLOf0AIoToOEXWR\naYtHbW0tSkpKrvlcTU0NampqrrtkcvUvq4yMDNhsNqxbtw6FhYWwWCxIT0/HvHnzEBUVZUh28h3y\nix1AbQ3E5Bmqo1BbevYF3G7gQgHQvZfqNETURUK2N7OSWlVUVKQ6gl8z4pqtrK+De+GD0H74OERK\n5yZJmZmZrou7//IGkNgT2u3mOytppnE2K46x/pKSkrx2LL+e40F0I3LTOqBvSkCWDrMRaSMh875U\nHYOIvIDFgwKSdFZBfr4a2vfuUx2FPHHzUCD/GGSdS3USIuoiFg8KSPLj9yFGTeTzWExChIYBfZKB\no1zFlMjsWDwo4MiLFyBzt0Dcea/qKNQBIm0k5KG9qmMQURexeFDAkauXQ9x2F4SDdzGZibhlBOSh\nve2uNExEvo3FgwKKPHMC8ngexHTz3R0R8Hr0BRobgGLeUUZkZiweFDCklHB/8A7Ed++FCLGpjkMd\nJISASB3Byy1EJsfiQYHj0F6gsgJi/HTVSaiTxC28rZbI7Fg8KCBIdxPcH74L7Z7vQwSZ8ymnBGDQ\nEODkEcjGRtVJiKiTWDwoIMhdm4HQcGBouuoo1AUi3A7EJQJnTqiOQkSdxOJBfk/W1UGuWQFt9gN8\nyJgfEINugTzG9TyIzIrFg/yezFoH9EuBSB6kOgp5gRjI4kFkZiwe5NdkjZNLo/ublFTg1HHIhgbV\nSYioE1g8yK/Jz/8BMXwsl0b3IyIsAkjsAZw+rjoKEXUCiwf5LemshNy6HuI7maqjkJeJgWmQx3m5\nhciMWDzIb8n1f29+EFxsguoo5GVi4C2QfGAckSmxeJBfkpfKIbM3QnxnjuoopIcBqcCZE5AN9aqT\nEFEHsXiQX5KffgAxdhpEVKzqKKQDERoGJPUGTh1THYWIOojFg/yOLL8ImbsV4lv3qI5COhIpabzc\nQmRCLB7kd+THqyAm3c7H3vs5kZIKmX9EdQwi6iAWD/Ir8uIFyC93Qtz+PdVRSG/9BzbP83A3qU5C\nRB3A4kF+RX70PsTU70BEOFRHIZ2JCAfgiAaKzqqOQkQdwOJBfkNeKIA8uAciY6bqKGQQkTwI8uRR\n1TGIqANYPMhvyHXvQWR8FyIsXHUUMkryQOAUiweRmbB4kF+QhV9DHj0IcdtdqqOQgUT/QZD5vKWW\nyExYPMgvyI/eh7h9FoQtVHUUMlJSL8BZCemsUp2EiDzE4kGmJ8+fgzz2FcTkb6mOQgYTWhDQbwAv\ntxCZCIsHmZ785G8Qt93Fsx0BqvlyC4sHkVmweJCpyZIiyENfQky7U3UUUkQkD4Tk0ulEpsHiQaYm\nP/kAYuq3m5/dQYGp/0DgzEnIJi4kRmQGLB5kWrK0GHJ/LsRt31UdhRQSYRFATBxQcEZ1FCLyAIsH\nmZZc/yHEpBkQ4RGqo5BiInkQL7cQmQSLB5mSrCiD/CIbYjpXKSUAfW4Cvj6pOgUReYDFg0xJfvZ3\niPG3QdgjVUchHyD6JEOezVcdg4g8wOJBpiMrKyB3beYTaOmfevYFigshG+pVJyGidrB4kOnIz/8B\nkT4ZIjJadRTyEcJiBRKSgIKvVUchonaweJCpSGcV5I4NEHfcrToK+RhebiEyBxYPMhW5cS3EreMh\nYuJVRyFf0zuZE0yJTIDFg0xDXq6F3PYpxB33qI5CPkj0uQny7CnVMYioHSweZBpy23qIwcMh4hNV\nRyFf1LMfcP4sZGOD6iREdAMsHmQKsqEecsNanu2gNomQECAuESg6qzoKEd0AiweZgty1GejVD6JX\nP9VRyIeJ3smQX3OCKZEvY/EgnyfdTZCf/R3at3i2g9rRhxNMiXwdiwf5vi93AfZIYECq6iTk4zjB\nlMj3sXiQT5NSwv3ph9DuuAdCCNVxyNf16gcUfg3Z2Kg6CRG1IVh1gK7Izs7G2rVrUVhYCKvVirS0\nNMybNw+JiW3f9VBXV4e//e1v2LNnD8rLy6FpGuLj4zFu3DjMmjULmsYu5lOO7Aca6oEho1QnIRMQ\ntlAgJh44f665hBCRzzFt8cjKysLSpUsBAAkJCaiurkZubi6OHj2Kl156CZGRrT88bNmyZdi+fTsA\noGfPnnC5XDh37hzef/99aJqGWbNmGfYeqH3uTz+EuOMeCBZC8pDo1Q+y8AwnIhP5KFP+Nm9sbMSK\nFSsAAGPGjMGSJUuwePFi2Gw2VFZWYvXq1W3ue+TIEQDA0KFD8corr+C1116DzWYDAJSWluofnjwm\nTx8HSoogRk9SHYXMJKk3UMhbaol8lSmLR35+PpxOJwAgPT0dABAdHY2UlBQAwP79+9vc9+abbwYA\nHDhwAI8//jh+/vOfw+VyYcCAATzb4WPc6z+EmD4LIti0J+ZIAdGjD2QhHxZH5KtMWTzKyspaXjsc\njpbXVy6v3OjMxU9/+lNMnDgRAFBQUIDS0lIEBwejd+/esNvtOiWmjmoqOgucOAwx8XbVUchsevTm\nImJEPsyv/ikppWx3m48//hjbt2/HgAED8Itf/AI1NTX41a9+hU2bNkFKiZ/+9KfX7ZOXl4e8vLyW\njzMzM1lSdOZa9SeE3D4ToXF8GJxerFarX/4cy/CbUFldhYggDSIsXHUcvx1nX8IxNsaqVataXqem\npiI1tXNLHJiyeMTGxra8rqqquu51XFxcq/vV19fj/fffB9B8icbhcMDhcGDw4MHYtWsXvvrqq1b3\na22Ar1zqIe+TVRWQO7MgFv0/NHKcdWO32/335zixJ5zHD0MkD1KdxL/H2UdwjPVnt9uRmZnplWOZ\n8lJLcnIyIiIiAAA5OTkAgPLychw/fhwAMGzYMADAokWLsGDBAqxcuRJA8620brcbQPM8EaB5ourZ\ns82nZUNCQox7E9QmufkTWMZOhXBEqY5CJsV5HkS+y5RnPIKDgzF37lwsW7YMubm5mD9/PpxOJ1wu\nFxwOR8sk0eLiYpSWlqKiogJAc2O7+eabceTIEezatQv5+fmor6/HpUuXAABTpkxR9ZboG7KuDnLr\neoQ8/zpqVYch8+I8DyKfZcriAQAZGRmw2WxYt24dCgsLYbFYkJ6ejnnz5iEq6tp/KV+94uUvfvEL\n/OMf/8CePXtQVlaG4OBgJCcn4/bbb2fx8AFy1yYgeRCCknoDPHVKnSSS+sD91V7VMYioFUJ6MiOT\nrlNUVKQ6gt+R7ia4Fz4E7YGfwzEinddsdebP18XlpTK4Fz2KoMXLVUfx63H2FRxj/SUlJXntWKac\n40F+av9uINwO3HSz6iRkdpExQFMTZNUl1UmI6F+weJDPcH++GtqM7/FhcNRlQojmeR6cYErkc1g8\nyCfIk0eAygpg+BjVUchPiB59IDnBlMjnsHiQT3Bv+AfE9JkQWpDqKOQvkvrwjAeRD2LxIOVkSRFw\nPA9ifIbqKORHeMaDyDd1qXi4XC64XC5vZaEAJTeshZg0AyLEpjoK+ZNv5njwxj0i39KhdTwOHTqE\n3bt349ixYygoKEBjY2PzQYKD0bNnT6SkpGD06NG45ZZbdAlL/kfWOCF3b4X2/Juqo5CfEeF2ICQU\nKC8FYvnMHyJf0W7xaGxsxMaNG7Fu3TqUlpYiIiIC/fr1w+TJkxEREQEpJWpqalBSUoKdO3fi888/\nR1xcHO68805Mnz4dwXykOd2A3PY5xNDREFExqqOQP0rsARQXsHgQ+ZB2W8EjjzyChoYGTJ48GePG\njUP//v1vuP2pU6ewc+dOrF69Gh999BHefJP/kqXWycZGyM0fQ/vZs6qjkJ8S3XtCXiiEGDxcdRQi\n+ka7xWPmzJmYOnUqrFarRwfs378/+vfvj3/7t39DVlZWlwOS/5Jf7gTiEyF6J6uOQv6qWw/gQoHq\nFER0lXYnl86YMcPj0nE1i8WCGTNmdCoU+T8pJeTGtdAyvqs6Cvkxkdh8xoOIfEenJmA0NDTgwoUL\nuHz5Mmw2G7p37w6LxeLtbOTPTh0DqquAoaNUJyF/1r0ncJ5nPIh8SYeKx/Hjx/Hhhx/i0KFDLXe0\nAM13taSlpWH27NkYMGCA10OS/5Eb1kBMu5MLhpG+ouOA2mpIVy2ELUx1GiJCB4rH+vXr8e677wIA\nBg0ahD59+sBms8HlcuHrr7/GwYMHcfDgQTzwwAO8xEI3JMtKII8ehPbAfNVRyM8JTQMSkoDiIqDP\nTarjEBE8LB7Hjx/Hn//8ZwwaNAgPP/wwEhISrtumpKQEf/jDH/DnP/8Z/fr1Q0pKitfDkn+QWR9D\njJvGf4GSIUT3npDnCyBYPIh8gkcrl65duxaJiYn45S9/2WrpAICEhAQ888wz6NatG9auXevVkOQ/\npOsyZPZGiGl3qo5CgYJ3thD5FI+Kx7FjxzB58uR2J5BarVZMnjwZx44d80o48j9y5yZg4C0Qcd1U\nR6FA0a07UHJedQoi+oZHxaO2thbR0dEeHTA6Ohq1tbVdCkX+SbrdkJs+4i20ZCgR3x2SxYPIZ3hU\nPKKiolBY6Nm98IWFhYiKiupSKPJTefsAmw246WbVSSiQdEsCSs7zYXFEPsKj4jF06FBs2rQJJSUl\nN9yupKQEmzZtwtChQ70SjvyLO+uj5ltohVAdhQKIiHAAQjSvG0NEynlUPO6++2643W4sXLgQO3bs\nuGYND6D5QXI7duzAc889Bykl7r77bl3CknnJ4iLg65MQoyaqjkKBKKF78y21RKSckB6efzx69CgW\nL16MyspKWK1WJCUltazjUVRUhPr6ejgcDjz++OMYNGiQ3rmVKyriL7GOcL//JyDYAu2e+z3a3m63\nw+l06pwqsAXSGLuXvQKkDoM27jbDv3cgjbMqHGP9JSUlee1YHi8gNmjQILz66qvYsGED9u7di4KC\nArhcLthsNvTt2xcjRozA9OnTERER4bVw5B+k6zLkrs3QFr6qOgoFKt7ZQuQzOrRkenh4OGbNmoVZ\ns2bplYf8kMzZDKSkQsS2vgYMke4SugMHv1Cdgojg4RwPos6SUkJmfQyNC4aRQiIhqXmeEREp5/EZ\nj/3798Nms7XM33C5XHj77bev2y4+Ph5z5szxXkIyt6MHm+8oGHiL6iQUyBK6Axebb6nlXVVEanlU\nPPLy8vDb3/4Wjz/+eMvnGhoasHXrVlgsFmjaP0+c1NXVYfDgwUhNTfV+WjIdd9bHEFO/w1/2pFa4\nvfnP2up/viYiJTwqHps3b0bv3r0xevTo67729NNPIy0treXjJ554Aps3b2bxIMiyEuBEHsQPF6iO\nQgFOCAHEdQMuXmDxIFLM42e13HrrrR4dcPTo0XxWCwEA5JZPIcZOg7CFqo5CBMQnAqXFqlMQBTyP\nikd5efmOBm3rAAAgAElEQVR1T6W1WCwYO3bsdcujx8XFoaKiwnsJyZRkfV3zU2infkt1FCIAgIjr\nBnmRxYNINY8utQQHB6OhoeGaz9lsNjz66KPXbdvY2HjNnA8KTHL3NqDvAIgE7y06Q9QlcYlAwWnV\nKYgCnkcNISEhASdPnvTogCdPnrzu7AgFluZbaD+CNu07qqMQtRDx3SAvXlAdgyjgeVQ8hg8fjp07\nd+L8+Ruv/FdUVITs7GyMGDHCK+HIpPKPAHV1wODhqpMQ/VMc53gQ+QKPisedd96J0NBQ/PrXv0ZO\nTg7cbvc1X3e73di5cyeef/55hIWF4c47uVhUIJNZH0NM/TYEL7mRL4lNACpKIZuaVCchCmgezfFw\nOBx4+umn8dJLL+HVV19t8yFxUVFReOqpp+BwOPTOTT5KXiqDzNsH7f88pDoK0TWExQI4ooCK0uZb\na4lICY9XLk1OTsbixYvx+eefY+/evSgsLMTly5dbHhI3cuRITJ8+HeHh4XrmJR8nt34GMXoSRBh/\nDsgHXVnLg8WDSJkOPSQuLCyMD4mjNsnGBsjtn0F7/L9VRyFqlYhNgCwrAdfRJVKnyxfhm5qacPjw\nYdTW1nojD5mY/CIbSOoN0b2X6ihErYtNAMouqk5BFNC6XDycTieef/55nDp1yht5yMR4Cy35vJh4\noKxEdQqigMbbDsgr5JkTQNUlYMgo1VGI2iRiEyDLecaDSCUWD/IKuXU9xKQZEFqQ6ihEbYtN4BkP\nIsVYPKjLZG015Jc7ISZMVx2F6MZi4oBLZZBuruVBpEqH7mppjcPhwBtvvHHdw+IocMhdmyFSR0A4\n+DNAvk1YrEC4Hai8BETHqo5DFJC6fMZD0zTEx8fDYrF4Iw+ZjJSy+TLLFD6FlkyCE0yJlOryGQ+V\nsrOzsXbtWhQWFsJqtSItLQ3z5s1DYmLiDferrq7GBx98gC+++AIVFRUIDQ1Fr169cP/996Nv377G\nhPcXxw81/zkgVW0OIg+1rOVx082qoxAFJI+Lx7lz57B69WoUFBTAbrdj/PjxmDp1KoS4dime7du3\n44033sD777/v9bBXy8rKwtKlSwE0Pz23uroaubm5OHr0KF566SVERka2ul91dTWeeeYZlJSUQNM0\nJCUlQdM0nDp1CsXFxSweHSS3fAox+VvX/RwQ+azYeIB3thAp41HxOH/+PJ599lk0NTWhV69eKCoq\nwtKlS7F582Y8/vjjhs/vaGxsxIoVKwAAY8aMwYIFC1BRUYFHH30UlZWVWL16NR544IFW9125ciVK\nSkoQExODX/3qVy1nR6SUaGhoMOot+AVZWQF5eB+0+x5WHYXIczHxQNFZ1SmIApZHxeO9996DzWbD\nokWLWv6i3rZtG95++208++yzePbZZ5GUlKRr0Kvl5+fD6XQCANLT0wEA0dHRSElJwcGDB7F///42\n983JyQEAdOvWDa+++iqKiooQHx+PO+64A7fffrv+4f2I3LEBYuR4PpeFTEXEJsD91V7VMYgClkeT\nS0+cOIE77rjjmrkTkyZNwm9+8xtomoaFCxfi5MmTuoX8V2VlZS2vr34S7pXLK6Wlpa3uV1VVherq\nagDAkSNHUF5ejqioKBQWFuJ///d/8dlnn+mY2r9IdxPkts8gJnNSKZlMdFzzE2qJSAmPiofT6Wz1\nckqPHj3wX//1X4iNjcWiRYtueKbBCFLKG369qemf9+7b7Xa88cYbeO2115CSkgIAWL9+va75/MpX\nXwKR0RB9klUnIeqY6Digoqz97YhIFx5daklISMDZs61fE42KisKvf/1rvPjii/jd736HYcOGeTVg\na2Jj/3n/fVVV1XWv4+LiWt3P4XAgODgYjY2N6N69O0JCQgAA/fr1w/Hjx3HxYusTzvLy8pCXl9fy\ncWZmJux2e5ffh5lVZ3+OkBmzEKLTOFit1oAfY70F6hjLiAhUNtQhwhIMYQvV/fsF6jgbiWNsjFWr\nVrW8Tk1NRWpq5+5m9Kh4DB48GLt27cJ9992HoKDrl8QOCwvDL3/5S7z66qvYu1f/a6fJycmIiIhA\ndXU1cnJyMG7cOJSXl+P48eMA0FJ+Fi1ahIqKCowePRpz585FUFAQBg8ejIMHD+L8+fOoq6uDxWLB\nmTNnAKDNeSqtDfCVOSaBSF68APeJw3D/8AnU6zQOdrs9oMfYCAE9xlGxcJ77GiKxh+7fKqDH2SAc\nY/3Z7XZkZmZ65VgeFY8pU6agsrIS+fn5LZcl/pXFYsETTzyBv/zlL/j666+9Eq4twcHBmDt3LpYt\nW4bc3FzMnz8fTqcTLpcLDocDs2bNAgAUFxejtLQUFRUVLfvee++9OHz4MJxOJ372s58hJCSk5UzH\n7Nmzdc3tL+T2zyDGTIX45owRkelcmedhQPEgomt5VDySk5Px2GOPtbudpmlt3sbqbRkZGbDZbFi3\nbh0KCwthsViQnp6OefPmXTcf5eo1JpKTk/H888/jvffew4kTJ+B2u5GamorZs2dj8ODBhmQ3M9nQ\nALljI7Rf/FZ1FKJOE9GxkBVl4OozRMYTsr0ZmQAOHDiAPn36XPMXemNjI4KDr+8t58+fx8GDBzFj\nxgzvJvUxRUVFqiMo4d69DXL75wh6/L91/T48daq/QB5j94fvArZQaN/xzqnjGwnkcTYKx1h/3lwy\nw6O7Wl544QUcOnSo5WOn04l///d/v+ZzV5w8eRJvv/221wKSb5FbP4U2+Q7VMYi6JjoWuMQ7W4hU\n6PJD4ihwyMKzQHERMGyM6ihEXSKi4yB5Sy2REiwe5DG59VOICdMhWrnERmQq0bFcRIxIERYP8oh0\nXYbM3Qoxyb/n7lCA4CJiRMqweJBH5J7twIDBEDHxqqMQdZ09EqitgWyoV52EKOB4fM68uLi45Xks\ntbW1AIDCwkLYbLbrtiP/I7d9Bu2ue1XHIPIKoWlAVEzzWY+E7qrjEAUUj4vHqlWrrlkuFQDvXgkQ\n8uwpoLICSBuhOgqR90RGN/9cs3gQGcqj4vHggw/qnYN8mNz+OcSEDAjt+uXyiUwrMhqoLFedgijg\neLxkOgUmWVcHuXsbtOdeUx2FyKtEZAxkZQVXLyUyGCeX0g3JL3YAyYMgYjmplPwMz3gQKeFR8aio\nqMCjjz6K995774bbvffee1iwYAEqKyu9Eo7Uk9s/gzbpdtUxiLwvMhq4VNH+dkTkVR4Vj08//RTV\n1dWYOXPmDbebOXMmnE4nPvnkE6+EI7Vk4VmgtAS4ZZTqKEReJ6KaL7UQkbE8Kh779u3DuHHjEBoa\nesPtQkNDMX78eOzdu9cr4Ugtuf0ziPEZEEGcVEp+KDKGl1qIFPCoeFy4cAG9e/f26IC9evXiWh5+\nQDbUQ+ZugZiQoToKkT6ivrmdlogM5fHkUimlV7cj3yb37gR63wQRn6g6CpE+IiKBy7WQjQ2qkxAF\nFI+KR0JCQsuqpe3Jz89HQkJCl0KRepxUSv5OaBpgdwBVl1RHIQooHhWPESNGIDs7G4WFhTfcrrCw\nEDt27MDIkSO9Eo7UkBcKgAuFwNDRqqMQ6SsyhpdbiAzmUfG46667EBoaiueffx47duxAU1PTNV9v\namrCjh07sGjRIoSGhuLOO+/UJSwZQ27/HGLsNIhgi+ooRPriWh5EhvNo5VKHw4FnnnkGL730EpYs\nWYKlS5ciKSkJNpsNLpcLRUVFqK+vR0xMDJ566ik4HA69c5NOZEMD5K7N0J56UXUUIt2JyGjIS1y9\nlMhIHj8krn///njllVewYcMG7N27FwUFBbh8+TJCQ0PRt29f3HrrrZg+fTrCwsL0zEs6k/tzgKTe\nEN2SVEch0h8vtRAZzuPiAQBhYWGYOXNmuwuJkXnJ7Z9DTJqhOgaRMSKjgHNnVKcgCijtzvGora3t\n9MG7si8ZT5YUAedOQwwfqzoKkSGEPQrSybtaiIzUbvF48MEHsWLFCpSUlHh80IsXL+Kvf/0rHnro\noS6FI2PJ7Rsgxk6FsHBSKQUIeyTg5LOliIzU7qWWBx98EKtWrcKaNWuQnJyMIUOGoH///ujWrRvC\nw8MBANXV1SgpKcGpU6dw8OBB5Ofno0ePHnjwwQd1fwPkHbKxAXLnJmhP/EZ1FCLjOCKBKhYPIiO1\nWzzGjBmD0aNHY+/evdi8eTPWrVuHxsbGVre1WCwYNmwY7rnnHowYMQJCcK64aRzYA3RLgujeS3US\nIuPYo3jGg8hgHk0u1TQNo0aNwqhRo9DQ0IBTp06hsLAQ1dXVAAC73Y4ePXqgf//+CA7u0HxV8hHu\n7Z9BTOSkUgowYeFAfR1kQwMvMRIZpMMtwWKxYODAgRg4cKAeeUgBWVoMnDkJ8dD/VR2FyFBCiOZl\n052VQEyc6jhEAcHjh8SR/5I7NkCkT4awhqiOQmQ8TjAlMhSLR4CTTU2Q2RshJvKBcBSg7FEAb6kl\nMgyLR6D76gsgNgGiZ1/VSYiUEI5ISN7ZQmQYFo8A5972Gc92UGDjpRYiQ7F4BDBZUQbkH4W4dYLq\nKETq8FILkaF0KR5Op1OPw5KXyZ2bIG4dDxFiUx2FSB0uIkZkKK8tuuFyufDll18iJycHBw4cwLvv\nvuutQ5MOpNsNmb0R2o+fUB2FSClhj4Sbl1qIDNOl4lFbW4svvvgCubm5OHjwINxuN4YOHYoFCxZ4\nKx/p5UQeYA0B+g5QnYRILa5eSmSoDheP6upq7NmzBzk5OTh06BCCgoIwbNgw/Od//idGjhwJm42n\n7c1A7tgAMT6Dy9oTOSI5x4PIQB4Vj/r6emzduhU5OTk4fPgwrFYrRo4cifnz52P37t244447kJKS\nondW8hJZWw15YA+0zB+pjkKknr15joeUkkWcyAAeFY/f//73OHr0KEaNGoUnn3wSQ4YMaXkmy8iR\nI/Haa6/hlltuwYwZfNaHGcjd24DBQyHsDtVRiJQT1hBAE0B9HcCJ1kS686h4hIWFYenSpbC08hAl\ni8WCxx57DG+88QZOnDiBn/zkJ7BarV4PSt4jd2yENuvfVccg8h1hdqDayeJBZACPbqd96KGHWi0d\nLQfRNDzyyCMIDQ3Fs88+i5KSEq8FJO+S504DVZeAwcNURyHyHRF2oKZKdQqigOBR8dA0z5b7+OEP\nf4iRI0di0aJFXQpF+pHZGyHGTYPQglRHIfId4Xagplp1CqKA4PUFxO699168/PLL3j4seYFsaIDM\n3QoxPkN1FCLfEmGHrObCh0RG0GXlUt5S65vk/hygZ1+I+ETVUYh8igh38FILkUH4rJYAIndsgJgw\nXXUMIt8T8c3kUiLSHYtHgJBlJcDX+RDDx6iOQuR7wiOAGhYPIiN47VktKmRnZ2Pt2rUoLCyE1WpF\nWloa5s2bh8REzy4lLF68GLm5uQCA9PR0PPbYY3rGVUpmb4IYPbF5zQIiula4Ayg4ozoFUUAw7RmP\nrKwsvP766zhz5gyio6MhpURubi6ee+45VFa2/9yFzZs3t5QOfyfd7uYn0Y7nZRai1ghOLiUyjCmL\nR2NjI1asWAEAGDNmDJYsWYLFixfDZrOhsrISq1evvuH+Fy5cwDvvvIOUlBTExMQYEVmtoweAsHCI\nPsmqkxD5pnA7L7UQGcSUxSM/Px9OZ/MvifT0dABAdHR0y/Ni9u/f3+a+brcbS5YsaVn0zNM1SsxM\n7tjISaVEN8LJpUSGMeXfumVlZS2vHY5/Pm8kMjISAFBaWtrmvqtWrcLJkyfxox/9CPHx8fqF9BGy\nxgl56EuI9CmqoxD5rnAHz3gQGcTUk0v/lZTyhl8/deoU1qxZg4kTJ2L8+PEeHzcvLw95eXktH2dm\nZsJut3c6p5HqsjegcXg6whO7q47SIVar1TRjbFYc43+SYWGovFyDiPAwr6/qy3HWH8fYGKtWrWp5\nnZqaitTU1E4dx5TFIzY2tuV1VVXVda/j4uJa3e/s2bNwu93IycnB7t27AQB1dXUAgD179uD73/8+\nli5ditDQ0Gv2a22Ar1zq8WVSSrg3fgRtzg9MkfdqdrvddJnNhmP8L2yhcJYUQ4R79y8wjrP+OMb6\ns9vtyMzM9MqxTHmpJTk5GREREQCAnJwcAEB5eTmOHz8OABg2rPkBaIsWLcKCBQuwcuXKa/ZvaGhA\nXV1dS+kAmud+1NfXt3vWxFTO5gOXa4BBQ1QnIfJ94ZznQWQEU57xCA4Oxty5c7Fs2TLk5uZi/vz5\ncDqdcLlccDgcmDVrFgCguLgYpaWlqKioAABMmTIFU6ZMueZYDz/8MEpLS/1yHQ+5YyPE+AyIAJhA\nS9RlEQ6gugrolqQ6CZFfM2XxAICMjAzYbDasW7cOhYWFsFgsSE9Px7x58xAVFXXNtkKIdo/nyTZm\nIhvqIfdsh7bw96qjEJkDb6klMoRpiwcATJgwARMmTGjz62+++Wa7x/BkGzOS+3OB3v0hYv3/zh0i\nb7iyiJh//ROEyPfwHLyfktnNl1mIyEOh4cDlWtUpiPwei4cfkuWlwOkTEMP4QDgij4WGA5erVacg\n8nssHn5I5myGuHU8RAgfCEfksTCe8SAyAouHn5FSNj+JdtxtqqMQmUtoGFBbozoFkd9j8fA3+UcA\nTQD9B6pOQmQuoeGQl1k8iPTG4uFn5M4siHEZfnd7MJHeBC+1EBmCxcOPyDoX5N6dEGOnqI5CZD68\n1EJkCBYPPyK/3AUkD4KIim1/YyK6Fm+nJTIEi4cfkdkboY2bpjoGkTmFhQO1vJ2WSG8sHn5ClhYD\nhWeAoemqoxCZU2gYz3gQGYDFw0/InVkQoyZBWCyqoxCZk8UKSAnZUK86CZFfY/HwA9Lthty5iUuk\nE3WBEOKbRcQ4wZRITywe/uBEHmALBXr3V52EyNxCw4BaXm4h0hOLhx+48kA4rt1B1EWhPONBpDcW\nD5OTrlrI/bsh0ierjkJkfrzUQqQ7Fg+Tk19kAwPTIBxRqqMQmR/vbCHSHYuHycnsTdDG84FwRN4g\nQsMhuXopka5YPExMFhcBJUVA2q2qoxD5B87xINIdi4eJyZ2bINInQwQHq45C5B/4vBYi3bF4mJR0\nN0Hu2gwxjpdZiLyGT6gl0h2Lh1kdOQg4oiB69lWdhMh/8FILke5YPEyqee0Onu0g8iYRGsbJpUQ6\nY/EwIVlbDXnoS4jRk1RHIfIvvNRCpDsWDxOSu7dDDB4GEW5XHYXIv4SG8VILkc5YPEyID4Qj0kmI\nDairU52CyK+xeJiMLDoLVJQCqcNURyHyP1YbUO9SnYLIr7F4mIzcuQlizFQILUh1FCL/Y7MBdSwe\nRHpi8TAR2dQEmbOFd7MQ6cVqA1wuSClVJyHyWyweZpL3JRCbAJHYU3USIr8kgoOBIA1obFAdhchv\nsXiYiNyZxZVKifQWEgq4eLmFSC8sHiYha6ohD++HuHWC6ihE/i0khBNMiXTE4mEScs+VtTsiVEch\n8m8840GkKxYPk5C7siDGTVMdg8j/hdiAusuqUxD5LRYPE5AXCoDSYiB1hOooRP4vhLfUEumJxcME\n5K4tEOmTIYK4dgeR7lg8iHTF4uHjpNsNmbOZl1mIDCJCbJAsHkS6YfHwdce+AsIjIHr2U52EKDDw\njAeRrlg8fBwnlRIZLITPayHSE4uHD5Ouy5D7d0OMnqQ6ClHgCLHxdloiHbF4+DD55U5gwGAIR7Tq\nKESBg5daiHTF4uHD5M4saLzMQmQsFg8iXbF4+ChZVgIUngGGjFYdhSiwsHgQ6YrFw0fJXZshbp0A\nYbGojkIUWFg8iHTF4uGDpJTNxWMsL7MQGY3reBDpi8XDF+UfBTQB9EtRnYQo8ISE8nZaIh0Fqw7Q\nFdnZ2Vi7di0KCwthtVqRlpaGefPmITExsc191q1bh3379qGoqAhOpxMOhwMDBgzA7Nmz0bt3bwPT\nt03uyoIYOw1CCNVRiAJPSAhvpyXSkWnPeGRlZeH111/HmTNnEB0dDSklcnNz8dxzz6GysrLN/dav\nX4+8vDyEhIQgISEB5eXlyM3NxcKFC1FaWmrgO2idrK+D/CIbYswU1VGIAlNIKOd4EOnIlMWjsbER\nK1asAACMGTMGS5YsweLFi2Gz2VBZWYnVq1e3ue+0adOwZMkSvPbaa3j11Vdx3333AQBcLhd2795t\nSP4bkQd2A32SIWLiVUchCkwhISweRDoyZfHIz8+H0+kEAKSnpwMAoqOjkZLSPCdi//79be57zz33\nICEhoeXjIUOGtLwODlZ/5YmTSokU4xkPIl2ZsniUlZW1vHY4HC2vIyMjAaBDl0w++uijluOMHTvW\nSwk7R1ZWAPlHIEaozUEU0EJCgHoXpJSqkxD5JfX/xPeijvyiaGxsxFtvvYXt27cjLCwMTz75JOx2\ne6vb5uXlIS8vr+XjzMzMNrftCtfWT+EeNRFhcbzMYrVadRlj+ieOcdsuQcAeaoOwWLt8LI6z/jjG\nxli1alXL69TUVKSmpnbqOKYsHrGxsS2vq6qqrnsdFxd3w/2dTideeuklHDt2DDExMXj66afRp0+f\nNrdvbYCvXOrxFikl3Js/gTb3J14/thnZ7XaOg844xjcQYoOz9CJEhKP9bdvBcdYfx1h/drsdmZmZ\nXjmWKYtHcnIyIiIiUF1djZycHIwbNw7l5eU4fvw4AGDYsGEAgEWLFqGiogKjR4/G3LlzAQAFBQV4\n8cUXUVJSgr59++Kpp55CTEyMsvfS4twpwHUZGNC5BklEXhRiA+rqgAjVQYj8jymLR3BwMObOnYtl\ny5YhNzcX8+fPh9PphMvlgsPhwKxZswAAxcXFKC0tRUVFRcu+L7/8MkpKSgAATU1NeOWVV1q+dttt\nt2HaNDUTO+XOLIixUyE0U067IfIvITag7rLqFER+yZTFAwAyMjJgs9mwbt06FBYWwmKxID09HfPm\nzUNUVNQ12169EFdjY2PL63Pnzl2z3fDhw/UN3QbZ2Ai5exu0p19U8v2J6F9cOeNBRF5n2uIBABMm\nTMCECRPa/Pqbb7553efeeOMNPSN1zqG9QEJ3iIQk1UmICAAsVqChXnUKIr/E8/o+wL0rC2Ic1+4g\n8hkWC4sHkU5YPBST1VXAkQMQt7Z95oaIDMYzHkS6YfFQTO7ZDpE2EiKM0+eJfIbFAtnQoDoFkV9i\n8VCs+W4WXmYh8iXCEsIzHkQ6YfFQSJ4/B1SUAYOHqY5CRFfjHA8i3bB4KCR3ZkGkT4YIClIdhYiu\nxjkeRLph8VBEupsgc7bwbhYiX2SxApzjQaQLFg9Vjh4EHFEQPdp+RgwRKcIzHkS6YfFQ5MoS6UTk\ngzjHg0g3LB4KSFct5ME9EKMnqY5CRK3hGQ8i3bB4KCD37gJS0iAcUe1vTETG4xwPIt2weCggd2VB\n49odRL7LYgUa+JA4Ij2weBhMlpUABWeAIaNURyGitlgsPONBpBMWD4PJnC0Qt46HsFhURyGiNgiL\nFZJzPIh0weJhICkl5K7NXCKdyNdxjgeRblg8jHT6OCAl0H+g6iREdCOc40GkGxYPAzWf7ZgKIYTq\nKER0I5zjQaQbFg+DyIYGyC+2Q4yZojoKEbWH63gQ6YbFwyhf7QGS+kDEdVOdhIjawzkeRLph8TCI\n+5vLLERkAjzjQaQbFg8DSGcVcOwQxMjxqqMQkSf4rBYi3bB4GEDu2QZxy60QoWGqoxCRJ6y81EKk\nFxYPA/BJtEQmE8zbaYn0wuKhM1l0FrhUDtw8VHUUIvJUUBAgAdnUpDoJkd9h8dCZzNkMkT4ZIihI\ndRQi8pAQgvM8iHTC4qEj6W6CzNkKMY5LpBOZDud5EOmCxUNPxw4BdgdEjz6qkxBRR3GeB5EuWDx0\nxEmlRCbGZdOJdMHioRPpugx5YDfE6EmqoxBRZ3ARMSJdsHjoRH65C7jpZghHtOooRNQZXDadSBcs\nHjqROZuhcVIpkXlZOMeDSA8sHjqQ5ReBs6eAoaNVRyGizuIcDyJdsHjoQOZsgRg5DsJiVR2FiDqL\nczyIdMHi4WVSSkg+iZbI9ITFCsniQeR1LB7eduYk0NQIJN+sOgkRdQXPeBDpgsXDy+SuLIix05qX\nXCYi8+IcDyJdsHh4kWxsgNyzHWLMFNVRiKirrCE840GkAxYPbzq0F+jeEyI+UXUSIuqqYD4kjkgP\nLB5e5N7ZfJmFiPwA53gQ6YLFw0tkdRVw9CDEyPGqoxCRN3COB5EuWDy8RO7ZAZE2EiIsXHUUIvIG\nK894EOmBxcNLrtzNQkR+IpjFg0gPLB5eIC8UAOUXgcHDVEchIm+xhQKXa1WnIPI7LB5eIHdthhg9\nCSIoSHUUIvISERkNWXVJdQwiv8Pi0UXS3QS5Mwti3G2qoxCRN0VGA5XlqlMQ+Z1g1QG6Ijs7G2vX\nrkVhYSGsVivS0tIwb948JCbeeB2Nzu7XqsMHgMhoiJ59O/cmiMg3RcYAl1g8iLzNtGc8srKy8Prr\nr+PMmTOIjo6GlBK5ubl47rnnUFlZ6fX92iJ3boIYn9GVt0JEvigsHGhshKxzqU5C5FdMWTwaGxux\nYsUKAMCYMWOwZMkSLF68GDabDZWVlVi9erVX92uLrHFCHvoSYvSkrr0hIvI5QghebiHSgSmLR35+\nPpxOJwAgPT0dABAdHY2UlBQAwP79+726X1vk7m0QaSMgwiM6/iaIyPdFxQCXKlSnIPIrpiweZWVl\nLa8dDkfL68jISABAaWmpV/dri8zexEmlRP4sMgayksWDyJtMPbn0X0kpddkvLy8PeXl5LR9nZmYi\nrOIiqp2VsKdPgNB4G623Wa1W2O121TH8Gse4fbXx3RDkqkFIF8aJ46w/jrExVq1a1fI6NTUVqamp\nnTqOKYtHbGxsy+uqqqrrXsfFxXl1v9YGuHr9amDMVFTXcIEhPdjt9pbLYqQPjnH73KHhaCg+j/ou\njBPHWX8cY/3Z7XZkZmZ65VimvNSSnJyMiIjmeRU5OTkAgPLychw/fhwAMGxY8wqiixYtwoIFC7By\n5coO7ecJmbMFYvIML7wbIvJZUTEAL7UQeZUpz3gEBwdj7ty5WLZsGXJzczF//nw4nU64XC44HA7M\nmi2DNs8AAAyoSURBVDULAFBcXIzS0lJUVFR0aD+PDLoFIiZej7dHRD5CRMbAzbtaiLzKlMUDADIy\nMmCz2bBu3ToUFhbCYrEgPT0d8+bNQ1RU1DXbCiE6td+NaLd/z2vvhYh8VFQ0z3gQeZmQnZ2RGeCK\niopUR/BrvGarP45x+6SzEu6FDyHo93/t9DE4zvrjGOsvKSnJa8cy5RwPIiJDhNuB+jrIirL2tyUi\nj7B4EBG1QWgaxIy74V7+Zqdv1yeia7F4EBHdgPjOHKCiFPLDdyBdl1XHITI9004uJSIyggi2QJu/\nEPLDv8D99I8gBg8D+iRDxHcHEhKBiEggPALCYlUdlcgUWDyIiNohYuIhfvw4ZHkp5JH9QMHXcJ84\nDJScB2qcQG01ANH8RNsQGxBsASwWINiCalsomoQAgoIBLQgiKAgIDgYcUUBIKFBdCbjdzfsLAEI0\n/wdx1Wtc9TkAQmv+s2Wbq15DANpV27a2zTXfo6033eYXPPrUDb/Q2rG7cIy6kBC46+pusG1bh27t\nCx08RitfEBNvb/7fmVrF4kFE5CEREwcxPuO6z0spgYb65gJSXwc0NAKN9UBDA0KsFrgrKwF3E9DU\nCNnkBpoamm/TdbmA7r0ALQiABCS++fOq/658/urXkIBbXrst/mUb6f5mW3fztlcf++r9W9PmdJZW\nvtDW3BdvHKPNQ1y7fZPVAtTVe/792vx0W9t28D1yPtANsXgQEXWREAKwhjT/9y8sdjtcV93q2eY/\nnKnTwng7ralwcikREREZhsWDiIiIDMPiQURERIZh8SAiIiLDsHgQERGRYVg8iIiIyDAsHkRERGQY\nFg8iIiIyDIsHERHR/2/v/mOiruM4jr+OgSHgeR4EYrOQhj+SMawFtKSmsdpsa/4e0h/ZH63WctM/\naq211Vj+QUb/tDVbm226pZFLSwNtrT8wmkdzIJsKGaYROugA5ZDd4ld/0H0HKnr3vbsPd93z8df3\nvt/v5/zce+9xL78/7gtjCB4AAMAYggcAADCG4AEAAIwheAAAAGMIHgAAwBiCBwAAMIbgAQAAjCF4\nAAAAYwgeAADAGIIHAAAwhuABAACMIXgAAABjCB4AAMAYggcAADCG4AEAAIwheAAAAGMIHgAAwBiC\nBwAAMIbgAQAAjCF4AAAAYwgeAADAGIIHAAAwhuABAACMIXgAAABjCB4AAMAYggcAADCG4AEAAIwh\neAAAAGMIHgAAwJjk2Z6AHWNjY/rmm2/U2Niovr4+zZ8/X6WlpaqsrFRqauqM4/x+vw4dOqSOjg79\n/fff8vv9crvdeuyxx7RhwwY5nU6DnwIAgMQTl0c8Pv30Ux0+fFher1c5OTkaHBxUQ0ODampq7jrO\n5/OpoaFBV65c0YIFC5Senq6enh7V19frgw8+MDR7AAASV9wd8fjjjz/0888/S5JefvllPfvsszpz\n5ow+/PBDnT9/Xs3NzSopKbnj2JSUFL344ouqqKhQWlqaxsfH9fHHH+vXX3/VlStXdPnyZeXl5Rn8\nNAAAJJa4O+LR0tJiLQcCxqOPPqqUlBRJUmtr64xjXS6XXnjhBaWlpUmSkpKSVFhYaG0PvAcAAIiO\nuAsefX191vL8+fMlSQ6HQ/Pmzbtt+734/X79+OOPkqRHHnlEDzzwQARnCgAAbhUzp1oOHTqkI0eO\n3HWf9957L2L/Xn9/v2pqatTV1aXFixdr586dM+577tw5nTt3znq9detWLVq0KGJzwZ0FwiSihxqb\nQZ2jjxpHX11dnbW8cuVKrVy50tb7xEzwyM/P19NPPz3jdofDIZfLpczMTGvdjRs35HK5NDExIZ/P\nJ0nTts/k0qVLqqmp0fXr17V8+XK99dZbSk9Pn3H/WwtcV1enrVu3BvOxYBM1jj5qbAZ1jj5qHH2R\nrHHMBI+SkpIZLwqdqri4WF999ZUkyePx6LnnntOZM2c0MjJibZcmj2hUV1fL4XCoqqpKjz/+uCSp\nublZn3zyif755x+Vl5frtddeU3JyzJQBAID/tbj7xs3Pz9eTTz6ppqYmffHFFzpx4oR6enokTV6n\nEQgvY2NjunbtmiRpeHhYkjQwMKDa2lpJkxeWXrt2bdrpm1deeYW7WgAAiKK4Cx6S9MYbbyg3N1eN\njY3q7e2V0+lUWVmZKisrb9vX4XBYy6Ojo9by+Pi4fv/992n7BQLKvdg9r4XgUePoo8ZmUOfoo8bR\nF8kaOyYmJiYi9m4AAAB3EXe30wIAgPhF8AAAAMYQPAAAgDEEDwAAYExc3tUSLU1NTfruu+/U3d2t\nOXPmqLCwUFVVVVq4cGFUxiUiO7X6+uuvdfjw4TtuO3jwoJKSyM8BFy5c0LfffqvOzk4NDg5KkrZs\n2aLNmzffcyx9HDy7daaXg3fs2DG1tLTo6tWr8vl8cjqdKigo0ObNm/Xggw/edSy9HBy7NQ63jwke\n//npp5/02WefSZKys7M1NDQkj8ej9vZ27dmzx3ouTKTGJaJwa+V0OpWTk2O9djgc026XxuTTm8+e\nPavc3FzrCzEY9HFo7NY5gF6+txMnTsjr9WrhwoWaO3eurl69Ko/Ho7Nnz6q2tlZZWVl3HEcvB89u\njQPs9jHBQ5O/7/Hll19KksrKyrRr1y4NDAxo586dunHjho4cOaLt27dHbFwiikStVq1apddff93A\nbOPXU089pYqKCo2Pj+ull14Kagx9HDo7dZ6KXr63tWvXqry8XNnZ2ZKk48eP68CBA/L7/Wpubta6\ndetuG0Mvh8ZOjaey28cED0mdnZ3Ws15KS0slSQsWLNDSpUvV1tam1tbWiI5LRJGolcfj0S+//KL0\n9HQtWbJElZWV/NLsLTIyMiRNPnk5WPRx6OzUeSp6+d42bdo07XVRUZG1PNNjLujl0Nip8VR2+5gT\nipL6+vqsZafTaS0HDsl5vd6IjktE4dYqKSlJLpdL2dnZun79ulpaWvTuu+/q8uXLUZlvIqGPzaKX\n7Tl+/LikyR594okn7rgPvRyeYGocEE4fc8TjLuz+qCs/Bhu8YGpVXl6udevWWU8Qbmtr0+7duzUy\nMqKTJ0/q1VdfjfY0ExJ9HHn0cuhGR0e1d+9enTp1SmlpaXrzzTc1b968kN6DXr67UGscbh8TPCRl\nZmZay1MvFAssz3SBjd1xiSicWt16JXpRUZEyMjI0NDTE/2AigD42h14Ojc/n0549e9TR0SG32623\n335bDz300Iz708uhC7XGUvh9zKkWSQ8//LB1zvb06dOSpP7+fv3222+SpOLiYklSdXW1du3apYMH\nD4Y0DvZrLEnff/+9BgYGrNdtbW0aGhqSJN1///1G5v9/Qh+bQS+H56+//tI777yjjo4O5eXlaffu\n3bd9IdLL4bFTYyn8PuaIhyYvotm2bZs+//xzeTwe7dixQz6fT36/X06nU+vXr5ck9fT0yOv1WgUP\ndhzs11iS6uvrtX//fmVlZem+++5Td3e3JCk1NVXPP//8rHyeWNXc3KwDBw5MW1dfX6/GxkYVFBRo\nx44d9HEE2KlzYB96OTgfffSRent7JUljY2Oqra21tj3zzDNau3YtvRwmOzWWwu9jgsd/KioqlJqa\nqmPHjqm7u1spKSkqLS1VVVWVXC7XtH2n3qccyrhEZ7fGGzdu1OnTp9XV1aXe3l5lZ2dr2bJl2rRp\nk3Jzc01/jJg2PDxs/SEJuHnzpm7evHnbYWb62D67daaXgzc6Omotd3V1Tdu2atWqaa/pZXvs1jjc\nPnZMcNUNAAAwhGs8AACAMQQPAABgDMEDAAAYQ/AAAADGEDwAAIAxBA8AAGAMwQMAABhD8AAAAMYQ\nPAAAgDEEDwAAYAzBAwAAGMND4gDEvP7+ftXV1SknJ0cjIyMqKirS+fPntXHjxtmeGoAQccQDQEwb\nHBzU+++/r9WrV2vDhg0qKSlRdXU1T3MF4hTBA0BM279/vxYvXqzCwkJJUkZGhsbGxrR06dJZnhkA\nOwgeAGKWz+dTU1OTysrKrHXt7e1yu93KzMycxZkBsIvgASBmXbx4UePj41qxYoW1rr29XQUFBbM4\nKwDhIHgAiFkjIyOaM2eOsrKyrHUdHR2cZgHiGMEDQMwqKChQUlKSRkdHJUk//PCD/vzzT4IHEMe4\nnRZAzHK73dq+fbv27dsnt9ut4eFhJScnKz8/f7anBsAmggeAmLZmzRqtWbNGknT06FHl5eUpOZk/\nXUC84lQLgLhx6dIlTrMAcY7gASBudHZ2ckcLEOc4Xgkg5rW2turo0aPyer1qaGhQWlqaiouLZ3ta\nAGxwTExMTMz2JAAAQGLgVAsAADCG4AEAAIwheAAAAGMIHgAAwBiCBwAAMIbgAQAAjCF4AAAAY/4F\nQhB54HSyOSkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x19c87278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def CE2GO_lambda_of_params(beta, mu, epsilon, Ldinvsq, kf, q):\n",
    "    \"\"\"Given parameters, calculate the eigenvalue lambda.\"\"\"\n",
    "    \n",
    "    funres = lambda lamb: CE2GO_disp_relation_residual(beta, mu, epsilon, Ldinvsq, kf, q, lamb)\n",
    "    domainmin = -2*mu\n",
    "    domainmax = 100\n",
    "    lamb = scipy.optimize.brentq(funres, domainmin, domainmax)  # 1D root finder\n",
    "    return lamb\n",
    "\n",
    "def CE2GO_disp_relation_residual(beta, mu, epsilon, Ldinvsq, kf, q, lamb):\n",
    "    \"\"\"Compute the residual for the dispersion relation.  lamb must be real.\"\"\"\n",
    "    mf = Ldinvsq / (kf*kf)\n",
    "    lhs = lamb + mu\n",
    "    rhs = q*q*epsilon/mu * kf**4 * (lamb + 2*mu) * (1+mf) * (1-q*q/(kf*kf*(1+mf))) * CE2GO_polarint(beta, mu, epsilon, Ldinvsq, kf, q, lamb)\n",
    "    res = rhs - lhs\n",
    "    return res\n",
    "\n",
    "def CE2GO_polarint(beta, mu, epsilon, Ldinvsq, kf, q, lamb):\n",
    "    \"\"\"Carry out the polar integral.  Same as for asymptotic WKE\"\"\"\n",
    "    fun = lambda phi: CE2GO_polarint_integrand(beta, mu, epsilon, Ldinvsq, kf, q, lamb, phi)\n",
    "    out, err = scipy.integrate.quad(fun, 0, 2*np.pi, epsabs=2e-6, epsrel=2e-6)\n",
    "    out = out / (2*np.pi)\n",
    "    return out\n",
    "\n",
    "def CE2GO_polarint_integrand(beta, mu, epsilon, Ldinvsq, kf, q, lamb, phi):\n",
    "    \"\"\" As explained above, integrand is A/(B+iC)^2.  Keep only the real part.\n",
    "    So use A(B^2 - C^2) / (B^2 + C^2)^2.\"\"\"\n",
    "    c = np.cos(phi)\n",
    "    s = np.sin(phi)\n",
    "    mf = Ldinvsq / (kf*kf)\n",
    "    A = (1 + mf - 4*c*c) * s * s\n",
    "    B = (lamb + 2*mu) * kf**2 * (1+mf)**2\n",
    "    C = -2*beta*q*c*s\n",
    "    \n",
    "    out = A*(B*B - C*C)  / (B*B + C*C)**2\n",
    "    return out\n",
    "\n",
    "# Now, let's go and calculate the dispersion relation lambda(q)\n",
    "beta = 1\n",
    "mu = 0.02\n",
    "epsilon = 1\n",
    "Ldinvsq = 1\n",
    "kf = 1\n",
    "\n",
    "q1 = np.logspace(-4, 0, 150)\n",
    "q2 = np.linspace(1.02, 2.2, 150)\n",
    "q = np.concatenate([q1, q2])\n",
    "lamb = np.zeros_like(q)\n",
    "\n",
    "for j in range(len(q)):\n",
    "    lamb[j] = CE2GO_lambda_of_params(beta, mu, epsilon, Ldinvsq, kf, q[j])\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.plot(q, lamb)\n",
    "plt.xlabel(r'$q$')\n",
    "plt.ylabel(r'$\\lambda$ (CE2-GO)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### WKE\n",
    "For the WKE, the dispersion relation is calculated similarly (details and code are omitted here).  However, a brief note is warranted on how the dispersion relation in Figure 2 is obtained.  Because the WKE is invalid for calculating the dispersion relation, it is not wholly surprising to see strange behavior.  When one calculates the dispersion relation about an equilibrium that balances forcing and dissipation, one finds that for some values of $q$, the eigenvalue $\\l$ becomes complex, even though the correct answer is that $\\l$ is real.  The source of this can be traced to a linearization of the $F/(\\b - U'')$ term, which gives $F U_1'' / \\b^2$.\n",
    "\n",
    "The WKE dispersion relation in Figure 2 is obtained by neglecting the $F U_1'' / \\b^2$ term, which is justifiable through an alternate procedure.  There are two ways of obtaining the equilibrium incoherent spectrum that we linearize about.  The first way to realize the equilibrium is the route we have been using, which is a balance between external forcing and linear dissipation.  An alternate route is to remove forcing and dissipation, in which case *any* homogeneous spectrum is an equilibrium.  This alternate procedure yields effectively the same dispersion relation, although within the WKE linearization it has the effect of removing the $FU_1''/\\b^2$ term in the linearization.  This procedure yields a real $\\l$ and is the one shown in Figure 2.\n",
    "\n",
    "The main point here is solely to demonstrate quantitatively that the WKE is not correct, which is to be expected because one had to assume that $U$ varied slowly in time in order to derive the WKE.\n",
    "\n",
    "\n",
    "### CE2\n",
    "The CE2 dispersion relation is obtained in a similar way (details and code omitted).  Since there is no $k_y$ derivative, one does not use an integration by parts but rather a shift of the integration variable, after which the isotropic form for $W_H$ can be inserted.  For details, see Srinivsan and Young (2012) or Parker (Ph.D. thesis, section 3.2.3)."
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}