{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Tutorial: Rapid Teukolsky mode amplitudes"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Waveforms are constructed from simulated trajectories by computing Teukolsky mode amplitudes $A_{\\ell m k n}$ at the trajectory points $(p, e, x_I)$, conditioning on the MBH spin $a$. As the $A_{\\ell m k n}$ are expensive to compute directly, in FEW we perform interpolation over pre-computed datasets in order to obtain mode amplitudes on a timescale of milliseconds.\n",
                "\n",
                "In this short tutorial, we will review the tools provided in FEW for generating waveform mode amplitudes."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "### Spline interpolation"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "FEW provides spline interpolation routines that are valid over a fixed domain of Kerr eccentric equatorial orbits, i.e. for varying $(a, p, e)$ and $|x_I|=1$. As mode amplitudes do not vary strongly with respect to $a$, we perform bicubic interpolation with respect to $(p,e)$ and linearly interpolate with respect to $a$. Mode indices $(\\ell, m, k, n)$ in the set $\\ell \\in [2, 10]$, $m \\in [-\\ell, \\ell]$, $k = 0$ (equatorial) and $n \\in [-55, 55]$ are supported. A notable restriction is that currently only a single input (i.e., a float) for $a$ is supported.\n",
                "\n",
                "For ease of waveform construction, we interpolate mode amplitudes that have been projected from the spin-weighted spheroidal harmonic basis to the spin-weighted _spherical_ harmonic basis (see e.g. [2102.02713](https://arxiv.org/abs/2102.02713), [2310.19706](https://arxiv.org/abs/2310.19706) and references therein). This basis admits the symmetry \n",
                "\n",
                "$$ A_{\\ell, -m, -k, -n} = (-1)^{\\ell + k} A_{\\ell m k n}^*,$$\n",
                "\n",
                "which we exploit in waveform generation in order to reduce the number of mode amplitude evaluations required. \n",
                "\n",
                "Once amplitude classes are instantiated, the necessary data files for interpolation are downloaded automatically. These can be reasonably large (a few GB), so ensure you are on a reasonably fast internet connection for your first run."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 1,
            "metadata": {},
            "outputs": [],
            "source": [
                "from few.amplitude.ampinterp2d import AmpInterpKerrEccEq\n",
                "import numpy as np\n",
                "\n",
                "kerr_amp_spline = AmpInterpKerrEccEq()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Specific modes can be selected by providing a list of tuples of $(\\ell, m, k, n)$ values as the `specific_modes` kwarg.\n",
                "\n",
                "Notice this returns a dictionary with keys as the mode tuple and values as the mode values."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 2,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "text/plain": [
                            "({(2, 2, 0, 0): array([0.21239204-0.06664877j]),\n",
                            "  (7, -3, 0, 1): array([1.38197053e-06-1.77646921e-06j])},\n",
                            " array([0.21239204-0.06664877j]))"
                        ]
                    },
                    "execution_count": 2,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "# Get mode amplitudes at an example point\n",
                "a = 0.87\n",
                "p = 7.0\n",
                "e = 0.4\n",
                "xI = 1.0\n",
                "specific_modes = [(2,2,0,0),(7,-3,0,1)]\n",
                "specific_amplitudes = kerr_amp_spline(a, p, e, xI, specific_modes=specific_modes)\n",
                "specific_amplitudes, specific_amplitudes[(2,2,0,0)]"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "If you call `kerr_amp_spline` without the key `specific_modes`, it will output all the modes:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 3,
            "metadata": {},
            "outputs": [],
            "source": [
                "all_amplitudes = kerr_amp_spline(a, p, e, xI)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "One can easily retrieve a specific mode via the use of `special_index_map`."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 4,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "text/plain": [
                            "True"
                        ]
                    },
                    "execution_count": 4,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "# Extracting the (2,2,0) mode\n",
                "mode_location = kerr_amp_spline.special_index_map[(2,2,0,0)]  # or special_index_map_arr[2,2,0,0]\n",
                "np.allclose(\n",
                "    specific_amplitudes[(2,2,0,0)], \n",
                "    all_amplitudes[0,mode_location]\n",
                "    )"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "The arguments of `kerr_amp_spline` can also be arrays (with the exception of the spin, which must always be a float).\n",
                "\n",
                "The data grids cannot generally be rectilinear in $(a, p, e)$ space, (for example to avoid computing the modes inside the separatrix). Instead, we start from a data points on a rectilinear grids $(u, w, y, z)$ all with ranges $[0, 1]$, and then map into space via the `apex_of_uwyz` utility:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 5,
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "Total trajectory points: 10201\n"
                    ]
                }
            ],
            "source": [
                "import matplotlib.pyplot as plt\n",
                "from few.utils.mappings.kerrecceq import kerrecceq_backward_map, z_of_a, y_of_x\n",
                "\n",
                "# A rectilinear grid\n",
                "uv = np.linspace(1e-5, 1 - 1e-5, 101)\n",
                "wv = np.linspace(1e-5, 1 - 1e-5, 101)\n",
                "\n",
                "a = 0.84\n",
                "z = z_of_a(a)\n",
                "y = y_of_x(1)\n",
                "\n",
                "# For a given rectilinear grid, compute the associated (a,p,e,x) grid\n",
                "u_grid, w_grid = np.asarray(np.meshgrid(uv, wv, indexing=\"ij\")).reshape(2, -1)\n",
                "a_grid, p_grid, e_grid, xI_grid = kerrecceq_backward_map(\n",
                "    u_grid, w_grid, np.ones_like(u_grid) * y, np.ones_like(u_grid) * z, kind=\"amplitude\"\n",
                ")\n",
                "\n",
                "print(\"Total trajectory points:\", w_grid.shape[0])"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 6,
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "Indices of interest: [1165 3495]\n",
                        "True\n",
                        "False\n"
                    ]
                }
            ],
            "source": [
                "# All the mode amplitudes for the given (a,p,e,x) data points.\n",
                "teuk_all_amps = kerr_amp_spline(a, p_grid, e_grid, xI_grid).reshape(\n",
                "    uv.size, wv.size, -1\n",
                ")\n",
                "\n",
                "# Specific mode amplitudes for the given (a,p,e,x) data points.\n",
                "specific_modes = [(2, 2, 0, 0), (7, -3, 0, 1)]\n",
                "teuk_specific_amps = kerr_amp_spline(a, p_grid, e_grid, xI_grid, specific_modes=specific_modes)\n",
                "\n",
                "# We can find the index to these modes to check\n",
                "inds = np.array([kerr_amp_spline.special_index_map[lmn] for lmn in specific_modes])\n",
                "print(\"Indices of interest:\", inds)\n",
                "\n",
                "# Making sure they are the same\n",
                "print(\n",
                "    np.allclose(\n",
                "        teuk_specific_amps[(2, 2, 0, 0)].reshape(uv.size, wv.size),\n",
                "        teuk_all_amps[:, :, inds[0]],\n",
                "    )\n",
                ")\n",
                "\n",
                "# to check -m modes we need to take the conjugate\n",
                "print(\n",
                "    np.allclose(\n",
                "        teuk_specific_amps[(7, -3, 0, 1)].reshape(uv.size, wv.size),\n",
                "        np.conj(teuk_all_amps[:, :, inds[1]]),\n",
                "    )\n",
                ")\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 7,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAF4CAYAAAC/wIoGAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWONJREFUeJzt3XtUFHeeN/53A9IYTeN446KIxMRLwkQd2ERgTLxEfEjGyeRJVmbcR6KBWTl4GUM0I+M8ER03ZDIJIRchulFZE+OyXnKZsyTKTqJinOSJBPc40USjTkBtJDg7giaCQv3+4EeHpruhq7su36p6v87pc+yyqutbXV3feveHb1XbJEmSQERERERkEiF6N4CIiIiISEkMuERERERkKgy4RERERGQqDLhEREREZCoMuERERERkKgy4RERERGQqDLhEREREZCoMuERERERkKgy4RERERGQqDLhEREREZCoMuERERESkioMHD2LOnDmIjY2FzWbD22+/3ecyra2tWL16NeLj42G32zFmzBhs2bJF1nrDAmwvEREREVGvrl69iokTJ2LhwoV4+OGH/Vpm7ty5uHjxIjZv3oxbb70VjY2NuHHjhqz1MuASERERkSoyMjKQkZHh9/zvv/8+Dhw4gDNnzmDw4MEAgNGjR8ter+UCbkdHBy5cuICbb74ZNptN7+YQkQlJkoSWlhbExsYiJMScI8HYlxLpK5h+5tq1a2hrawtq3T2Pe7vdDrvdHvBrdnn33XeRnJyMZ599Fq+//joGDBiAn/70p/jd736H/v37+/06lgu4Fy5cQFxcnN7NICILqK+vx8iRI/VuhirYlxKJQW4/c+3aNYwaNQDffNMR8DoHDhyIK1euuE1bs2YNCgsLA37NLmfOnMGhQ4cQERGBt956C01NTcjLy8Pf/vY3WeNwLRdwb775ZgDAtKFZCAsJ92uZtrGxfr9++MkLAbUr0PX1Rom2+KJUGwOl5rYRBetGRxv2N21z9Tdm1LVt9fX1cDgcOreGyFwevu2JPucJtJ9pa2vDN990YP8nwzFwoPy/vly5ImHa3Y0ex74S1Vug869DNpsN27dvR2RkJACguLgYjzzyCDZs2OB3FddyAberpB4WEu53wO0Ii/D79f19TaXW50v4F+cABdriixJtDIYS7zOR2sz8p/uubXM4HAy4REHKiFns9lzOOS7QfmbgQBsG3hzIEKrOyq9ax35MTAxGjBjhCrcAMGHCBEiShHPnzuG2227z63UsF3DlaBsv70+L4V+c03ydetC7jUq8z0RERHrpGWjpe2lpadi5cyeuXLmCgQMHAgBOnjyJkJAQWUMxGHB90DvEBUPNAGjk94WIiEgPVg60V65cwVdffeV6fvbsWRw9ehSDBw/GqFGjUFBQgPPnz2Pbtm0AgHnz5uF3v/sdFi5ciLVr16KpqQkrV67EY489xovMghVIiBOlemv2cMvqLRERGYGVQ213R44cwfTp013P8/PzAQCPPvooysvL4XQ6UVdX5/r/gQMHoqqqCkuXLkVycjKGDBmCuXPnYv369bLWy4Dbg14hToTw2BvR20dERKQnBlrvpk2bBkmSfP5/eXm5x7Tx48ejqqoqqPVaNuC2jY1F2FdN3z8PIsCJUlUUpR1qMfv2ERGRcTDQis2yARcQpyopSjt8EaF9DLdERKQ3hlrjsHTAVUKwwUup8KhWABQh3BIREemBgda4GHCDwKqiNvg+ExGRFhhozYMBV0es3hIREemLodacGHADZPaqoijh1uzvMxERacssgXb3qecRGfma3s0QFgNuAES55y2gTgAUJdwSEREpwSyhlvzHgEvCYvWWiIgCxVBrbQy4MrF6qw2GW//2Bd8nIqJODLTUHQOuDAy3pIZg3veeyzLwEpGVMNSSLwy4fmJw0I4V3mu1vkx0f10rvI9EZD0MteQPBlw/KBUUWL21Nq3f4671MegSkdEx1JJcDLgaETlAitQ2s4UxEd5bBl0iMhoGWgoWA24fRAsFSrdHhADWRbT3Ohgiva9dGHSJSGQMtaQkBtxeiDY0gcRnhH3dNn4kQy4RCYGhltTCgOuFqCd/Vm/FJdJ76Q+GXCLSC0MtaYEBtwczh8juRG2X0Rj5fWTIJSKtMNQq6z3nBjQ3N+vdDKFZNuCGn7yAsJBwVdehZPgxcxAx6rYZOdx2YcglIrUw1JKeLBtwrUykYGbEcCXS+6cEhlwiUgpDLYmCAVclolZvzRbOtGbW948hl4iCwWBLomHAFZyZQ4eRts2swbY7hlwikoOhlkTGgKsCUcOQqO0SHd83IqJODLVkFAy4CuPQBP8YoVIo2numBVZxiagnhloyohC9G0DWY4QAZcVwS8ZWWlqKhIQEREREICkpCdXV1b3Ov337dkycOBE33XQTYmJisHDhQly6dEmj1pIRZMQsZrglw2LAVRCrt+Zg9ffL6ttvRBUVFVi+fDlWr16N2tpaTJ06FRkZGairq/M6/6FDh5CVlYXs7Gx8/vnn2LlzJz799FPk5ORo3HISTVeoZbAlo+MQBZMTLayIXL0V7b0i8ldxcTGys7NdAbWkpAR79+5FWVkZioqKPOb/+OOPMXr0aCxbtgwAkJCQgEWLFuHZZ5/VtN0kBoZZMiNWcBUiavWW/MNw647vh3G0tbWhpqYG6enpbtPT09Nx+PBhr8ukpqbi3LlzqKyshCRJuHjxInbt2oUHHnjA53paW1vR3Nzs9iBjY6WWzIwVXAWIGgZEa5eowV2094lIjqamJrS3tyMqKsptelRUFBoaGrwuk5qaiu3btyMzMxPXrl3DjRs38NOf/hQvv/yyz/UUFRVh7dq1iradtMdAa3zvOTfo3QRDYAVXMEqFQNFCG8MtkbpsNpvbc0mSPKZ1OX78OJYtW4annnoKNTU1eP/993H27Fnk5ub6fP2CggJcvnzZ9aivr1e0/aQuVmtJLwcPHsScOXMQGxsLm82Gt99+2+9lP/roI4SFhWHSpEmy18sKbpA4NMGYGGzJLIYOHYrQ0FCPam1jY6NHVbdLUVER0tLSsHLlSgDAnXfeiQEDBmDq1KlYv349YmJiPJax2+2w2+3KbwCphoGWRHD16lVMnDgRCxcuxMMPP+z3cpcvX0ZWVhZmzpyJixcvyl4vA24QRA1JorVLtOAu2vsjKt4T1xjCw8ORlJSEqqoqPPTQQ67pVVVVePDBB70u8+233yIszL37Dw0NBdBZ+SVjY7AlLfQch+/rS3BGRgYyMjJkv/6iRYswb948hIaGyqr6dmHAFQSHJmhDtPeHSAn5+fmYP38+kpOTkZKSgk2bNqGurs415KCgoADnz5/Htm3bAABz5szBL3/5S5SVlWH27NlwOp1Yvnw57rrrLsTGxuq5KRQghlqS6+2WiYiQ+sle7tqV6wD2IS4uzm36mjVrUFhYqEjbtm7ditOnT+ONN97A+vXrA3oNBtwAMSgZD/cZmVVmZiYuXbqEdevWwel0IjExEZWVlYiPjwcAOJ1Ot3viLliwAC0tLXjllVfwxBNPYNCgQZgxYwZ+//vf67UJFCAGW9JLfX09HA6H67lSQ5hOnTqFVatWobq62uMvTXIw4AZA6aDE6q36RHtviJSWl5eHvLw8r/9XXl7uMW3p0qVYunSpyq0iNTDUkggcDodbwFVCe3s75s2bh7Vr12Ls2LFBvRYDrkyiBiVR2yUCvjdEZAYMtmR2LS0tOHLkCGpra7FkyRIAQEdHByRJQlhYGPbt24cZM2b49VoMuDoTqcqpJFG2i+GWiIyOwZaswuFw4NixY27TSktL8cEHH2DXrl1ISEjw+7UYcGUQNSyJ1i6GWyKi4DDUkllcuXIFX331lev52bNncfToUQwePBijRo1yuwg2JCQEiYmJbssPHz4cERERHtP7woDrJzXCkhJBkCHOO74vRGREDLZkNkeOHMH06dNdz/Pz8wEAjz76KMrLyz0uglUKA64fGJb8J0L1lvuLiIyGwZb8YcSf6Z02bVqv99f2dhFsd4WFhQHdfowBtw9qhSVWb9XB94SIjITBlkgdDLi9EDncikjv7WK4JSKjYLAlUhcDrg+ihyXR2sdwS0TUNwZbIm0w4HqhZlji0ATl8f0gIpEx1BJpjwG3GwalwOhZveU+IyJRMdgS6ceyAbdtbCw6wiI0XacZq7cMt+al97ATIqNisCXSn2UDrhEx0H2P7wURiYbBlkgcDLgaMWM1TK9tYrglIpEw2BKJhwHXIBjqOvF9ICJRMNgSiYsBVwOs3iqD4VY7ZvzMEimFwZZIfAy4KuOFZcoQ7T0gIuthsCUyjhC9G0C9Y7Dje0BE+mO4Jb2959ygdxMMRfeAW1paioSEBERERCApKQnV1dW9zr99+3ZMnDgRN910E2JiYrBw4UJcunRJo9bKY8Y/82q9TQy32jPj55YoUBkxixluiQxI14BbUVGB5cuXY/Xq1aitrcXUqVORkZGBuro6r/MfOnQIWVlZyM7Oxueff46dO3fi008/RU5OjsYt14bVw53Vt5+I9MNgS2Rsugbc4uJiZGdnIycnBxMmTEBJSQni4uJQVlbmdf6PP/4Yo0ePxrJly5CQkIAf//jHWLRoEY4cOeJzHa2trWhubnZ7aMGMVTAtt4nhloj0wGBLZA66Bdy2tjbU1NQgPT3dbXp6ejoOHz7sdZnU1FScO3cOlZWVkCQJFy9exK5du/DAAw/4XE9RUREiIyNdj7i4OEW3Qy2iBTwzBnbyxP1MVsVgS2QuugXcpqYmtLe3Iyoqym16VFQUGhoavC6TmpqK7du3IzMzE+Hh4YiOjsagQYPw8ssv+1xPQUEBLl++7HrU19cruh3eBBsSRAu3WrP69hORdhhsicxJ94vMbDab23NJkjymdTl+/DiWLVuGp556CjU1NXj//fdx9uxZ5Obm+nx9u90Oh8Ph9lCTGStgHJpgDWb87BL1hsGWyLx0C7hDhw5FaGioR7W2sbHRo6rbpaioCGlpaVi5ciXuvPNOzJ49G6WlpdiyZQucTqcWzVadaAGP4ZbIGOTckWbBggWw2WwejzvuuEPDFuuHVVsi89Mt4IaHhyMpKQlVVVVu06uqqpCamup1mW+//RYhIe5NDg0NBdBZ+dUbhyYEzsrbLgJWb41N7h1pXnzxRTidTtejvr4egwcPxj/+4z9q3HJtMdgSWYeuQxTy8/Px2muvYcuWLThx4gQef/xx1NXVuYYcFBQUICsryzX/nDlzsGfPHpSVleHMmTP46KOPsGzZMtx1112IjY3VazNMS6vQw3BLFBy5d6SJjIxEdHS063HkyBH8z//8DxYuXKhxy7XDYEtkLbr+VG9mZiYuXbqEdevWwel0IjExEZWVlYiPjwcAOJ1OtwrEggUL0NLSgldeeQVPPPEEBg0ahBkzZuD3v/+9XpvgwuptYKy63SJh9dbYuu5Is2rVKrfpvd2RpqfNmzfjvvvuc/W93rS2tqK1tdX1XKtbLgaLwZbImnQNuACQl5eHvLw8r/9XXl7uMW3p0qVYunSpyq2Sx4wBwYzbRJ64n40vkDvSdOd0OvHee+/hzTff7HW+oqIirF27Nqi2aonBlsyEP9Mrn+53USDxqpgcmkBkPHLuSNNdeXk5Bg0ahJ/97Ge9zqfHLRcDxXBLRLpXcI2OFbDAMNzqj59dcwjkjjRdJEnCli1bMH/+fISHh/c6r91uh91uD7q9alI72Hrrt3gcEYmJAVdnogU9LTpr0bbZinhSNo/ud6R56KGHXNOrqqrw4IMP9rrsgQMH8NVXXyE7O1vtZqpKrWDrT1/Vcx4eW0RiYMANAi8sk8+K20yktvz8fMyfPx/JyclISUnBpk2bPO5Ic/78eWzbts1tuc2bN+Puu+9GYmKiHs1WhBrhNph+qmtZBl0ifXEMboDM2HmZcZvIE/ez+WRmZqKkpATr1q3DpEmTcPDgwV7vSAMAly9fxu7duw1bvVXjnrZt40cq9iWcX+aJOh08eBBz5sxBbGwsbDYb3n777V7n37NnD2bNmoVhw4bB4XAgJSUFe/fulb1eVnADoERAEK3z49AEa2C4NS+5d6SJjIzEt99+q3Kr1KFGsFVD2/iRPObI8q5evYqJEydi4cKFePjhh/uc/+DBg5g1axaefvppDBo0CFu3bsWcOXPwySefYPLkyX6vlwGXNMFwS0TBEm04gr+vz5BLVpaRkYGMjAy/5y8pKXF7/vTTT+Odd97BH//4RwZcNbF6K59e23t5jH9XfEeebu17JhPgSZaMzChVWyKr6PljL2rdaaWjowMtLS0YPHiwrOUYcGUwY7hVmx7b62+w7T6/2UMuwy0Z1czpRQCA3m9iJo/W/RKruCSiDy+ORdgV+YH0xtVWAPsQFxfnNn3NmjUoLCxUpnHdPP/887h69Srmzp0razkGXD+ZtXMy03bJDbY9lzVryDXTPiZr6Qq3SrJakYFILfX19XA4HK7nalRvd+zYgcLCQrzzzjsYPny4rGUZcP2gVEAQrWM1y9CEYIItEYmnZ7A1ax9MZGQOh8Mt4CqtoqIC2dnZ2LlzJ+677z7Zy/M2YX1g9SswWpxILo+xKxpuzRiU+fklo1Gjagsw3JJxvefcoHcTNLdjxw4sWLAAb775Jh544IGAXoMV3F4oGQ5E61zVDD5ahVvqHcMtGY2Zwy3H4ZJVXblyBV999ZXr+dmzZ3H06FEMHjwYo0aN8vghmh07diArKwsvvvgipkyZ4voZ8v79+yMyMtLv9bKC64OZw62a1N5Wpau2ZsUTKRmJGj/a0MVK/S+RiI4cOYLJkye7bvGVn5+PyZMn46mnngLg+UM0GzduxI0bN7B48WLExMS4Hr/61a9krZcV3B6sEAyMuo0Mtv4x6v4la5o5vQjh8B1Eg/k8M9wS6W/atGmQJMnn//f8IZr9+/crsl7LBtzwkxfQcfstnf82+J/r5TDitjLYEpmTWkMSAPH6XiLSlmUDLsBKl5IYbsXAzzQRwy0RcQyuqkTrZNUKPwy3YmC4JSIi6sSAaxFGCz8Mt/IYbf8SqUW0wkIXHqNE2mLAVYmonazS1NhOhlt5eOIko1PqAjOr9LtE1DcGXBWI1skaZWgCbwEmH8MtGZVatwUjIgIYcClAaoRbkofhloysqw/hz/ASkRoYcBUmWidrhBAkSriNPN2qdxP8ZoT9SqQV0fpdIqVY8Wd6lcKAa2JGGJogSrg1EoZbsgIzfc7NtC1ERsGAqyArVBEYbolIScH2KVbod4lIPgZchYjWyapRMWC41R8rQUTfE63fJSJxMOCS5kQMt0YYf8twS2Q8PG6J9MGAqwDRqggiV29FDLdGwJMkWYk/n3fR+l0iEgsDrskw3MonevWW4ZbImHjsEumHATdIrCL4R9RwKzqeIMnMAv18s98lor4w4AZBtE5W1Ootw21gGG5JjtLSUiQkJCAiIgJJSUmorq7udf7W1lasXr0a8fHxsNvtGDNmDLZs2aJRazuJ1ocSkXmE6d0AEpcVwq2owxMYbkmOiooKLF++HKWlpUhLS8PGjRuRkZGB48ePY9SoUV6XmTt3Li5evIjNmzfj1ltvRWNjI27cuKFxy+UzSijmMUykLwbcAInWySrdmVoh3IqKJ0aSq7i4GNnZ2cjJyQEAlJSUYO/evSgrK0NRUZHH/O+//z4OHDiAM2fOYPDgwQCA0aNH97qO1tZWtLZ+/4Wwubk54PZmxCwGfPQx/PwTkRI4RCEAZg+3SjBCuBWxeiviviSxtbW1oaamBunp6W7T09PTcfjwYa/LvPvuu0hOTsazzz6LESNGYOzYsVixYgW+++47n+spKipCZGSk6xEXF6fodvhDtL7XFx7HRPpjBVcmo3SwwQh2G40QbkXEkyIFoqmpCe3t7YiKinKbHhUVhYaGBq/LnDlzBocOHUJERATeeustNDU1IS8vD3/72998jsMtKChAfn6+63lzc3PAIdcK/ShRsN5zbtC7CYbGgGtwIg5NMALRqrcMtxQsm83m9lySJI9pXTo6OmCz2bB9+3ZERkYC6Bzm8Mgjj2DDhg3o37+/xzJ2ux12u35fXo3SN/FYJhIDhyjIYJQOVk+s3srHEyIFY+jQoQgNDfWo1jY2NnpUdbvExMRgxIgRrnALABMmTIAkSTh3Tr/Po9GPBaO3n8hMGHD9JGK4Fa16a5RwK1L1lidEClZ4eDiSkpJQVVXlNr2qqgqpqalel0lLS8OFCxdw5coV17STJ08iJCQEI0eK19eJ2P8SkdgYcP0gYufKcGt8DLeklPz8fLz22mvYsmULTpw4gccffxx1dXXIzc0F0Dl+NisryzX/vHnzMGTIECxcuBDHjx/HwYMHsXLlSjz22GNehydQ33g8E4mFY3D7IGK4FY2Rwq0o1VueDElJmZmZuHTpEtatWwen04nExERUVlYiPj4eAOB0OlFXV+eaf+DAgaiqqsLSpUuRnJyMIUOGYO7cuVi/fr1em0BEpCgG3F6IGm5Fqt4aKdyKguGW1JCXl4e8vDyv/1deXu4xbfz48R7DGkQkaj/cHY9pIvFwiIIPRuhUlWClcCtC9ZYnQiLveGwQkZIYcHtoGz9S6HCr5ElA5O00I57AieQxQh/F45qodwcPHsScOXMQGxsLm82Gt99+u89lDhw4gKSkJEREROCWW27Bq6++Knu9DLjdiN6ZitSRsnorj0j7joiUweOaqG9Xr17FxIkT8corr/g1/9mzZ3H//fdj6tSpqK2txW9+8xssW7YMu3fvlrVey47BbRsbi46wCL2boRsrDU3QG0+CRN+bOb1I7yYQkYYyMjKQkZHh9/yvvvoqRo0ahZKSEgCd9+g+cuQInnvuOTz88MN+vw4ruAYhytAEI4ZbPau3DLdEgeFf1IjE1tzc7PZobVXmXPvnP/8Z6enpbtNmz56NI0eO4Pr1636/jmUruCQfw608PAES+cdox4rR2kvG855zg+rrqLswFCH95f8lu+O7awCAuLg4t+lr1qxBYWFh0O1qaGjw+BXGqKgo3LhxA01NTYiJifHrdRhwDUCE6q0Rw62eeAIkCpzo1VsiAurr6+FwOFzP7XblcoLNZnN7LkmS1+m9YcAVnAjh1qj0qt4y3BKZF49vok4Oh8Mt4ColOjoaDQ0NbtMaGxsRFhaGIUOG+P06HINLfWL11n88+RGZF49vIvWlpKR4/AjNvn37kJycjH79+vn9OroH3NLSUiQkJCAiIgJJSUmorq7udf7W1lasXr0a8fHxsNvtGDNmDLZs2aJRa7UlQvXWqOFWj+otT35EwbPaX5qIzO7KlSs4evQojh49CqDzNmBHjx51/Xx4QUEBsrKyXPPn5ubi66+/Rn5+Pk6cOIEtW7Zg8+bNWLFihaz16jpEoaKiAsuXL0dpaSnS0tKwceNGZGRk4Pjx4xg1apTXZebOnYuLFy9i8+bNuPXWW9HY2IgbN25o3HJjYbhVH8MtUWCMcuwYpZ1Eojly5AimT5/uep6fnw8AePTRR1FeXg6n0+kKuwCQkJCAyspKPP7449iwYQNiY2Px0ksvybpFGKBzwC0uLkZ2djZycnIAACUlJdi7dy/KyspQVOR5r8T3338fBw4cwJkzZzB48GAAwOjRo7Vssmb07kyNGm71oPe+IiJ18RgnCty0adNcF4l5U15e7jHt3nvvxWeffRbUenUbotDW1oaamhqPe52lp6fj8OHDXpd59913kZycjGeffRYjRozA2LFjsWLFCnz33Xc+19Pa2upxrzbRiTA0wai0rt7yxEekHKv1V0SkHt0quE1NTWhvb/d6r7OeV891OXPmDA4dOoSIiAi89dZbaGpqQl5eHv72t7/5HIdbVFSEtWvXKt5+I7Da0AStMdwSmR+PcyJj0v0iM2/3OvN1n7OOjg7YbDZs374dd911F+6//34UFxejvLzcZxW3oKAAly9fdj3q6+sV3wYl6d2ZGjncalm91Xs/EZH6eJwTGZduFdyhQ4ciNDTU673OelZ1u8TExGDEiBGIjIx0TZswYQIkScK5c+dw2223eSxjt9sVvfmwUQRSvWW49Q9PekSBmTnd89oKIiI16FbBDQ8PR1JSkse9zqqqqpCamup1mbS0NFy4cAFXrlxxTTt58iRCQkIwcqTxx24pFZw4jk09DLdEyul+PInWb/FYJzI2XYco5Ofn47XXXsOWLVtw4sQJPP7446irq0Nubi4Az3ujzZs3D0OGDMHChQtx/PhxHDx4ECtXrsRjjz2G/v3767UZitC7M2X1tm967yMi0gaPddLTe84NejfBFHS9TVhmZiYuXbqEdevWwel0IjExEZWVlYiPjwcAj3ujDRw4EFVVVVi6dCmSk5MxZMgQzJ07F+vXr9drE4TDoQlERERkdboGXADIy8tDXl6e1//zdm+08ePHewxrMDo9hyYYOdxqiRUdIvWINDyBxzqROeh+FwWiQHFoAhEpicc6kXkw4OqM1dvAMNwSERGRL7oPUaDgWS3caoXhlkg9/h5f3fsqNb/Y8ngnMhdWcHWkV4dq9HCrRfWWJzsymtLSUiQkJCAiIgJJSUmorq72Oe/+/fths9k8Hl988YWGLe7U2xd0rfoqHu9E5sOAqxPe8zYwDLdEnioqKrB8+XKsXr0atbW1mDp1KjIyMtzuQuPNl19+CafT6Xp4+7EcvRj9izgR6YsB18A4NEF5DLdkRMXFxcjOzkZOTg4mTJiAkpISxMXFoaysrNflhg8fjujoaNcjNDRUoxb3Tst+isc8kTkx4OqAQxMCo3b1lic6MqK2tjbU1NQgPT3dbXp6ejoOHz7c67KTJ09GTEwMZs6ciQ8//LDXeVtbW9Hc3Oz2MDoe80TmxYCrMb2GJjDc9o4nOjKqpqYmtLe3Iyoqym16VFQUGhoavC4TExODTZs2Yffu3dizZw/GjRuHmTNn4uDBgz7XU1RUhMjISNcjLi5OVjtnTi/ymOatHzN6X0VEYuBdFAzIauNu1cZwS2Zgs9ncnkuS5DGty7hx4zBu3DjX85SUFNTX1+O5557DPffc43WZgoIC5Ofnu543NzfLDrldRDjmRGgDEamHFVwNcWhCYHhrICLfhg4ditDQUI9qbWNjo0dVtzdTpkzBqVOnfP6/3W6Hw+FweyiNd00gq3vPuUHvJpgGA67BcGiCcniSIzMIDw9HUlKSx0+YV1VVITU11e/Xqa2tRUxMjNLNIyLSBYcoaESJMGW1oQkMt0T+yc/Px/z585GcnIyUlBRs2rQJdXV1yM3NBdA5vOD8+fPYtm0bAKCkpASjR4/GHXfcgba2NrzxxhvYvXs3du/erVmbe/ZnrN4SkZIYcDXAoQli4QmOzCYzMxOXLl3CunXr4HQ6kZiYiMrKSsTHxwMAnE6n2z1x29rasGLFCpw/fx79+/fHHXfcgf/8z//E/fffr9cmaILHPpF1MOAaBIcmKIMnODKrvLw85OXlef2/8vJyt+dPPvkknnzySQ1a5R9/+istfuSFiMyDY3BVpkegYrj1juGWyLp4/BNZCwOuAVht7K0aeHIjEoMexyKPfyLrYcBVkR4XlrF664knNyKxdO/XjN5nEZGYGHBVwqEJ8jHcEpHS2AcQWRMvMhOYlYYm8AISMSnxGWTAoGAF2j/ws0dkXQy4KuDQBDHw5NY3Lb5E+bMO7itzy4hZDHj5HLDfIiK1BDRE4Z/+6Z+wadMmnDx5Uun2GB6HJsjHoQnqaxs/0utDFCK3jYKnx/5kH0AkjtLSUiQkJCAiIgJJSUmorq7udf7t27dj4sSJuOmmmxATE4OFCxfi0qVLstYZUMAdOHAgiouLMX78eMTGxuIXv/gFXn31VXzxxReBvBz1IOdkwHDriSc2z8BoRGbYBvLEfUnk3XvODXo3QRUVFRVYvnw5Vq9ejdraWkydOhUZGRluPz7T3aFDh5CVlYXs7Gx8/vnn2LlzJz799FPk5OTIWm9AAXfjxo344osvcOHCBRQXFyMyMhIvvvgi7rjjDkv/ljl/jld/Vg23VgiDZtzG6upq/J//83+QkpKC8+fPAwBef/11HDp0SOeWqSOQ4zOQL8FW7QeIRFRcXIzs7Gzk5ORgwoQJKCkpQVxcHMrKyrzO//HHH2P06NFYtmwZEhIS8OMf/xiLFi3CkSNHZK03qLso3HzzzfjBD36AH/zgBxg0aBDCwsIQHR0dzEsaFocmyKd09dZqJzWzhT25jL79u3fvxuzZs9G/f3/U1taitbXzeGhpacHTTz+tc+vUp1b/ZbV+gEgvzc3Nbo+uPqy7trY21NTUID093W16eno6Dh8+7PV1U1NTce7cOVRWVkKSJFy8eBG7du3CAw88IKt9AQXcX//615gyZQqGDh2K3/72t2hra0NBQQEuXryI2traQF6SwKEJwbDKSc3ooU4tRnxf1q9fj1dffRX/+q//in79+rmmp6am4rPPPtOxZURkBeH14bDXyX+E14cDAOLi4hAZGel6FBUVeayjqakJ7e3tiIqKcpseFRWFhoYGr+1KTU3F9u3bkZmZifDwcERHR2PQoEF4+eWXZW1fQHdR+MMf/oBhw4ZhzZo1ePDBBzFhwoRAXsY0ODRBHoZbeaz02VBC9/dL5M/Gl19+iXvuucdjusPhwN///nftGyQguX2FyPubyGzq6+vhcDhcz+1234U3m83m9lySJI9pXY4fP45ly5bhqaeewuzZs+F0OrFy5Urk5uZi8+bNfrcvoIBbW1uLAwcOYP/+/Xj++ecRGhqKe++9F9OmTcO0adMsFXg5NEFfZj2hMdQqo+t9FPFzEhMTg6+++gqjR492m37o0CHccsst+jRKZV37g30YkfE5HA63gOvN0KFDERoa6lGtbWxs9KjqdikqKkJaWhpWrlwJALjzzjsxYMAATJ06FevXr/f7Wq+AhihMnDgRy5Ytw549e/DNN99g7969uOmmm7Bs2TIkJiYG8pKWxqEJgRExtATLaH9mNwoRhzAsWrQIv/rVr/DJJ5/AZrPhwoUL2L59O1asWIG8vDy9m2c4ZuwPiIwuPDwcSUlJqKqqcpteVVWF1NRUr8t8++23CAlxj6ehoaEAOiu//gr4hx5qa2uxf/9+7N+/H9XV1WhubsakSZMwffr0QF/ScNihysNw651IocsKRKnqPvnkk7h8+TKmT5+Oa9eu4Z577oHdbseKFSuwZMkSXdtmNHrvSyLyLT8/H/Pnz0dycjJSUlKwadMm1NXVITc3FwBQUFCA8+fPY9u2bQCAOXPm4Je//CXKyspcQxSWL1+Ou+66C7GxsX6vN6CA+4Mf/ABXrlzBxIkTMW3aNPzyl7/EPffc02ep2kyU6lCtVL1VillOZgy2+hIh6P7Lv/wLVq9ejePHj6OjowO33347Bg4cqFt7RMKf7yYyh8zMTFy6dAnr1q2D0+lEYmIiKisrER8fDwBwOp1u98RdsGABWlpa8Morr+CJJ57AoEGDMGPGDPz+97+Xtd6AAu7rr79uuUDbHcOtfDxZfY/BVix6B92bbroJycnJuqxbS+FfnEPb+JGK92Vm+cJLZGZ5eXk+h16Vl5d7TFu6dCmWLl0a1DoDCrg/+clPglqpkbEzlY9DEzox2IpN76BrVhkxiwF+9olIYwGPwbUiJU98VqneMtwy2BoNg66yut5POceBv/0G9xER+cKA6yeGW30Z8UTGYGtsDLpEpJX3nBv0boLpMOD6gSe4wChVvTXa+89gay4MumLi/iCi3gR0H1wrUboTtUr11orhVrT7rJKyuG+Dp1SfZqR+gYj0wQpuLxhuA2PVcEvmx2qu+njHFe3o2W/xGCK1MeB6wQNPf0bZBwy21iRi0C0tLcUf/vAHOJ1O3HHHHSgpKcHUqVP7XO6jjz7Cvffei8TERBw9elT9hgZJpPdcFEbsh5RqMz8P5AsDbjdqHiis3vrPCB2WEU8opLy28SOF+LxWVFRg+fLlKC0tRVpaGjZu3IiMjAwcP34co0aN8rnc5cuXkZWVhZkzZ+LixYsatpjkYH/jm9z3RoTjlbRh2YAbfvICwkLCNVkXw63/RO98eKKhnkSo5hYXFyM7Oxs5OTkAgJKSEuzduxdlZWUoKiryudyiRYswb948hIaG4u2339aotb711YeI3j8Ei/2L+vx9j83+WbMCywZcUpYVxs3x5EO90aua29bWhpqaGqxatcptenp6Og4fPuxzua1bt+L06dN44403sH79+j7X09raitbW74/z5uZmWe008pd3NbA/EZs/+4chWGwMuCqzSvVWCaJ2FjwRkb+6Pishx89ots6mpia0t7cjKirKbXpUVBQaGhq8LnPq1CmsWrUK1dXVCAvz7zRQVFSEtWvXBt3eQInaP/iL/Yj59LVPjf6ZNToGXEEYOdyaeWgCT0oUiLaxsUCjtuu02WxuzyVJ8pgGAO3t7Zg3bx7Wrl2LsWPH+v36BQUFyM/Pdz1vbm5GXFxc4A02OfYdxACsLwZcFVmhgzNruLXCviNzGDp0KEJDQz2qtY2NjR5VXQBoaWnBkSNHUFtbiyVLlgAAOjo6IEkSwsLCsG/fPsyYMcNjObvdDrtdvS/ivfUlIvYR3rDfIDn4eVEXA65KrDA0geGWSH/h4eFISkpCVVUVHnroIdf0qqoqPPjggx7zOxwOHDt2zG1aaWkpPvjgA+zatQsJCQmKtW3mdN8XuJkF+wsiMTHg6syo4VYJooVbnqjIqPLz8zF//nwkJycjJSUFmzZtQl1dHXJzcwF0Di84f/48tm3bhpCQECQmJrotP3z4cERERHhMV0ow/Zxo/QTAvoLICBhwVWCFzi/Y6q1oJy0r7DMyr8zMTFy6dAnr1q2D0+lEYmIiKisrER8fDwBwOp2oq6vTuZW+GeEuLOwjSC1/+rBA7yaYEgOuwjg0oW8ihVuetMgs8vLykJeX5/X/ysvLe122sLAQhYWFyjfKBNhHEBkTA65OrBpuRWLUE5fanx0z7WMyNj2/DBu1fyCiTgy4CmKH2DcRqrdG2U96fQnytV4GX7ICo/QPRNQ7BlyFcGhC3xhufTPCZ6JnGxl4SQm+Pkda9xei9g1EFJgQvRtQWlqKhIQEREREICkpCdXV1X4t99FHHyEsLAyTJk1St4F+YLjtG8Otp8tj7K6HERm9/URAZ78gWt9ARMHTtYJbUVGB5cuXo7S0FGlpadi4cSMyMjJw/PhxjBo1yudyly9fRlZWFmbOnImLFy9q2GJrMnq4FeXkZeYg2H3bWNml7gL53GvRZ4jSLxCROnSt4BYXFyM7Oxs5OTmYMGECSkpKEBcXh7Kysl6XW7RoEebNm4eUlBSNWuqbFaq3wbB6uLVildOK20yB0ePLECu2RNagW8Bta2tDTU0N0tPT3aanp6fj8OHDPpfbunUrTp8+jTVr1vi1ntbWVjQ3N7s9lGKFcBvMCcjK4ZYBrxPDLsmlZr/BYEtkHboNUWhqakJ7e7vHb6VHRUV5/KZ6l1OnTmHVqlWorq5GWJh/TS8qKsLatWuDbm9PVugojRpu9do3DHG963p/OISBtGaF/pqI3Ol+kZnNZnN7LkmSxzQAaG9vx7x587B27VqMHTvW79cvKCjA5cuXXY/6+vqg2yyXEYMPw63/WKGUh++XNcycXqR3EwAw3BJZlW4V3KFDhyI0NNSjWtvY2OhR1QWAlpYWHDlyBLW1tViyZAkAoKOjA5IkISwsDPv27cOMGTM8lrPb7bDblT2Zmn1oglErbFqeyIy4X0XDiq61edvvSn85Zrglsi7dAm54eDiSkpJQVVWFhx56yDW9qqoKDz74oMf8DocDx44dc5tWWlqKDz74ALt27UJCQoLqbQbYYfZFr+qtVvuFwVZ5DLrmpdfxwn6ajOJPHxbo3QTT0nWIQn5+Pl577TVs2bIFJ06cwOOPP466ujrk5uYC6BxekJWV1dnQkBAkJia6PYYPH46IiAgkJiZiwIABqrdXbqdpxDBktKEJWl0RzT+rq4/vsbUp1X8w3BKJR+5vHrS2tmL16tWIj4+H3W7HmDFjsGXLFlnr1PU+uJmZmbh06RLWrVsHp9OJxMREVFZWIj4+HgDgdDpRV1enZxNdGG57p1e4VZsR96PRsaJLgWK4JRJPIL95MHfuXFy8eBGbN2/GrbfeisbGRty4cUPWem2SJElKbIBRNDc3IzIyEvcNz0FYSLhfyzDc9s6M4daI+9CMjBpyb9y4hoOH1uHy5ctwOBx6N0cVXX2pr22cOb3I53GkxvhbhlsyomCGKPR1DPa13JjfPI3QiAjZ622/dg2nn/4N6uvr3dbr65qnu+++Gz/60Y/cfuNgwoQJ+NnPfoaiIs+LUd9//338/Oc/x5kzZzB48GDZ7eui+10URMdOs3dmC7f8M7lYuD+sgeGWSHs310u4+esAHvWdddG4uDhERka6Ht7CaiC/efDuu+8iOTkZzz77LEaMGIGxY8dixYoV+O6772Rtn65DFEQWaIdpxJNxoFUyrcOt2sGWxMVhC+QLwy2RPrxVcHsK5DcPzpw5g0OHDiEiIgJvvfUWmpqakJeXh7/97W+yxuEy4HrBcCsetU5iRtxnVnZ5jN0wn1lSH8MtkX4cDoffQyP8/c0DoPMWsDabDdu3b0dkZCQAoLi4GI888gg2bNiA/v37+7VODlHoxmq/UW6UcbcMt9Qdhy0YW89+J9C+xEp9NZFRyf3NAwCIiYnBiBEjXOEW6ByzK0kSzp3zv7+wbMBtGxvrCrRKBFujnXCtHG4ZkMyB+1Bsau4fhlsiY+j+mwfdVVVVITU11esyaWlpuHDhAq5cueKadvLkSYSEhGDkSP+PfcsGXCVZ6URr5HDLYGs+3J/GFkh/wnBLZCxyfvMAAObNm4chQ4Zg4cKFOH78OA4ePIiVK1fiscce83t4AsAxuEEz4glW9IvK1KrakjnxAjTrYLglMh65v3kwcOBAVFVVYenSpUhOTsaQIUMwd+5crF+/XtZ6GXAtxmrhlsHWOngBmjhmTi8CvBx7wewfkcKtkv0KP7NkBXl5ecjLy/P6f+Xl5R7Txo8f7zGsQS4G3CAYLTxZKdwabd+QMhhyjUNOv6JHuNWqD/FnPfxME8nHMbgBskqAYrglo7Hq/pfzW++HDh1CWloahgwZgv79+2P8+PF44YUXNGyt/7T6Se6eD5F4a5+I7SR5gvkVM+obK7gBMGKnInIFQKkTmBH3C6nDapVcub/1PmDAACxZsgR33nknBgwYgEOHDmHRokUYMGAA/vmf/1mHLdCWmfoKb9tipc8+kS8MuDIZsWMUeWiC1cJtS7z3G1sH4+avJcVf0wysdPFZcXExsrOzkZOTAwAoKSnB3r17UVZW5vXnMydPnozJkye7no8ePRp79uxBdXW1UAGXY/IDw9BLxIArixE7R1HDrdmDrRpBVu66GHw7mb2a2/Vb76tWrXKb3ttvvfdUW1uLw4cP93qVcmtrK1pbv38fm5ubZbWz+z7wp38xex+htZ7vg5mPCSKAAddvRuwkGW61o2Wg9VfPNlk58Jo55AbyW+9dRo4ciW+++QY3btxAYWGhqwLsTVFREdauXetXm0Q4NkVog8gYeMnsGHD9YKWO0gjhVoT9IWKg7Uv3Nlsx7Jo55ALyfuu9S3V1Na5cuYKPP/4Yq1atwq233opf/OIXXuctKChAfn6+63lzczPi4uKCb7gXwfQTIvQPRsTAS2bDgNsHo3aWgXRODLe9M2Ko9cWqYdeMITeQ33rvkpCQAAD44Q9/iIsXL6KwsNBnwLXb7bDbgz/++upnAu0njNpXi4qBl4yOtwnrhVE7TIZb5bTE21wPszL79vVk1OPal0B+690bSZLcxtgqSc1wxNtlaYO3JiOjYQXXCyMfwGYMt1rvDyuFve6sVNU1WyU3Pz8f8+fPR3JyMlJSUrBp0yaP33o/f/48tm3bBgDYsGEDRo0ahfHjxwPovC/uc889h6VLl+q2DYD8vsLIfbWRdX/fzXQckbkw4HZj9M5SxI7GKOHWqqHWl673w8xB10whV+5vvXd0dKCgoABnz55FWFgYxowZg2eeeQaLFi1StZ29fZlmuDUmhl0SlWUDbnOCHaHh5ukgRbxjghEuFGGw7Z3Zg66ZQq6c33pfunSp7tXaQDHYiovjdkkkHINrYVYOt1YbdxosM79fDEz687e/4L4yFo7ZJT1ZtoJrJiKNuxV5SIJZA5qWWuJtpqzmmqmSq7eMmMXAA2MUf10GJePiMAbSAwOuwZkl3DLYGodZhy0w5CqjZz/Q9Z766nf86TcYbs2DYbfTnz4s0LsJpseAa2AMt71jsFWXWau5FBylj2eGW/Ni2CU1MeAaFMOtbwy22jFbNZdVXG311Xcw3FoHwy4pjReZGRDDrXdmvhBKdGZ63xmqxMD9YF28OI2UwAquBYgWblm1NSczVXNZyVWWtz6ot/6D4YYAVnUpOKzgGowoB7ko4ZZVW/GYZX8wZAUvkP6K7zt5w6ouycWAayCiDE0QIdwy2IrNLPuGJ1R1+OpD+H5TX7qCLj8r1BcOUTAII4dbpYMtGQPvskBEauo6t4jyl00SCyu4BsBw24nh1njMUGlnpSg4PfsiVm9Jaazokjes4AqO4ZbB1gyMXs3lRWfqEuGHXoz8+aROvCiNumPAFZjVwy2Drbkw5JK3vkSUL8K9LW/kz61VcfgCMeCaiNLhllVbUprRQy75T4tgoVU/4W09/BwbA4OudTHgCkruwWiWcCt6sG0d1abYa9nrwhV7LSMxcshlFbd3M6cXATL6gUD7DBH6ie5tMOrn2UoYdK2HF5kJiOFWf62j2rw+jLYOUYm0r+XixSy+9XxvuvdNwfzqYRdRL1rsapeIbSN3vCBNH6WlpUhISEBERASSkpJQXV3t13IfffQRwsLCMGnSJNnrZMAVjBXDrd4nBtFCpkhtUZORwwBPkMGT8x7q3UfIwbBrDHoG3T99WKDLevVSUVGB5cuXY/Xq1aitrcXUqVORkZGBurq6Xpe7fPkysrKyMHPmzIDWy4ArEKuFWz1PAkYKkEZqq1wMAdbkb79h9KDIsCs+/nCE+oqLi5GdnY2cnBxMmDABJSUliIuLQ1lZWa/LLVq0CPPmzUNKSkpA62XAFYTRwm2wHYLWHb5ZqqJGb783Rj3584Tov0C+PBv1c+ELg674GHTlaW5udnu0tnrmmLa2NtTU1CA9Pd1tenp6Og4fPuzztbdu3YrTp09jzZo1AbePF5kJwIjhNlB6BFsz6r5dZrhYzagXnvGiM+96e0/86T/MHAS7ts2In3ersMoFaY6zrQgLk3+s3bjR+b7ExcW5TV+zZg0KCwvdpjU1NaG9vR1RUVFu06OiotDQ0OD19U+dOoVVq1ahuroaYWGBx1QGXJ3pfQCZMdyaNdT60rW9Rg+6Rg255F0gX8TNHGx7YtAVH7/A9q6+vh4Oh8P13G73nQ9sNvdjW5Ikj2kA0N7ejnnz5mHt2rUYO3ZsUO1jwNWR3j/kYKZwa7VQ641Zgq7R8CTYu+79TG99iJXCbXcMumKzSjU3EA6Hwy3gejN06FCEhoZ6VGsbGxs9qroA0NLSgiNHjqC2thZLliwBAHR0dECSJISFhWHfvn2YMWOGX+1jwNWJkcItg62xtI5qM2zIZRXXmqwabrtj0BUbg25gwsPDkZSUhKqqKjz00EOu6VVVVXjwwQc95nc4HDh27JjbtNLSUnzwwQfYtWsXEhIS/F43A64OGG6Dw1DbNyNXc40YclnFJaUw6IqNQVe+/Px8zJ8/H8nJyUhJScGmTZtQV1eH3NxcAEBBQQHOnz+Pbdu2ISQkBImJiW7LDx8+HBERER7T+8K7KGjMCuFWrauFzXgHAbUZ9f0yYkVP76uv5dxIfc+ePZg1axaGDRsGh8OBlJQU7N27V/E2+TM8wYj7Wgu864LYeMcF/2VmZqKkpATr1q3DpEmTcPDgQVRWViI+Ph4A4HQ6+7wnbiAYcDVklXCrNAbb4Bj1/ePJ3X9yb6R+8OBBzJo1C5WVlaipqcH06dMxZ84c1NbWatpu7uO+MeiKjUHXP3l5efjrX/+K1tZW1NTU4J577nH9X3l5Ofbv3+9z2cLCQhw9elT2OjlEQSN6hlstLiZTK9iScow8Ntco9Bqq0P1G6gBQUlKCvXv3oqysDEVFRR7zl5SUuD1/+umn8c477+CPf/wjJk+e7HUdra2tbve5bG5u9tmenn2Xtz6FoU0eIw7dsRIOXRAPK7gaMEq4DfSbqNInKqNWHI3AaO8tQ1DfAr2RencdHR1oaWnB4MGDfc5TVFSEyMhI16PnPTC7+NPfcb8GhtVc8bGaKw4GXJUZKdwGQsnO1mjhy8iM9D4b7YSu9QkukBup9/T888/j6tWrmDt3rs95CgoKcPnyZdejvr7e7f+Tc4rdnnf1PzzhK89ox4TVcNiCGDhEQUVmDrdKB1vSHocsqKc5wQ4c0nad/t5IvacdO3agsLAQ77zzDoYPH+5zPrvd3uuN3P3FcKYM3m1BfBy2oC/dK7giXvmrBIbbvrFiqz+jvP8MRb7JvZF6dxUVFcjOzsZ//Md/4L777lOzmaQSHhvi636e/dOHBTq2xFp0DbhGvfK3L2YNt0qN/2KwFYtR9gVP5N51v5F6d1VVVUhNTfW53I4dO7BgwQK8+eabeOCBBxRtk6/hCdyH6uDYXPFx2IL2dB2ioMWVv1ozQrjVu2pL4uFwBWOTcyN1oDPcZmVl4cUXX8SUKVNc1d/+/fsjMjIyqLaEf3HOa18USB8SaH9h1c8y77RA9D3dAm7Xlb+rVq1ym670lb9ybm0TLIZb3xhsxWeEkMsTuHeZmZm4dOkS1q1bB6fTicTExF5vpL5x40bcuHEDixcvxuLFi13TH330UZSXl2vdfDdK9BXdX0P0z7TSODaXqJNuAVerK3+Lioqwdu3aoNrqDzOGWwZb62HINa68vDzk5eV5/b+eobW3m6orpXt/409folZfYdWwy+OErE73i8yCvfK3oqKi1yt/+7q1TbAiT7fqFm7bxo8UNtxynK1xcb9RoAK9Wlyrz5zV+iWOyyUr062Cq8SVvzt37uzzyl+lbm3jTaCduVLh1l96hFsRjR75jSKv89dzwxR5HZGJXslldUpcPfum3voTvfqKrvWK/BlXCocskFXpFnC7X/n70EMPuaZXVVXhwQcf9Lncjh078Nhjj2HHjh2KX/krB4ckeBIl2CoVZOW8vhlDL0MuqUmE/sJqQZfHC1mJrndREOnKX3/pWbUFxAy3ep6o1A6z/urZDrMEXtFDLompr/5HhHDbnVWCLkMuWYmuAddIV/4G80skIodbowVbUQJtX7q30+hhV+SQyxO2WHreIsxb/yJauO1O5M+6UnjMkFXo/lO9ol3521OwP7HHcBs8o4RaX7rab/SgS2QFDLlE5qB7wBWRUr8brWW41WpIglbB1uih1hsjB12RT/o8WYujbfxIV19ktOptd1YYssDjhszOsgHXcbYVYWHq3UJFyzslmKVqa8ZQ641Rg67IIZf0lZxT3Oc8Rgm33Zn9M887LJCZWTbgqsls4VbNE5NVQq03Rgy6op7wWY0SmxHDbRdRP/NK4vFDZsSAqyCzjbdlsNWGEYMuUW/M9gMDDLmkhCOv5evdBEthwFWIqONtRaraMtT2bvTIbwwRckU92fMErS9f1y4YuXrbHcflEhmL7j/VawZmCrdq/JTl6JHfMNz6ySjvk1lCCylL7sWuRmT2z77Zqu9kXazgBkHUIQmiVG2NEtZEwyELgWMFSn/d+x+zhkFR/4qhFB5HZAas4AbITOFW6aotK7bKEP09NGt4IfKH2T//rOSS0bGCK5NSwRbQf0iC0qGWlCf6uFwRK1msPonB7AEQEPPzT0SdGHD9pGSwBZQNt3reIYHBVn2ih1wiILDxt331H/zc64tfFsnIGHB7oXSoBfQfksBga0wih1xWsahLV3/UWz8jp+/oPi8///pgyCWjYsD9/6kRZnsyQ9WWwVY/Iodc0fCkLKZg+g+RL75kyCUSj2UvMgs/eQHhX5xzPdSmZ7hV4iIyXjgmBlH3gRXGW1JwlPrs8hjQBy86o2CUlpYiISEBERERSEpKQnV1tc959+zZg1mzZmHYsGFwOBxISUnB3r17Za/TsgFXK23jR+oWbhlszYn7wz88IevDW5+j9GdW1H6JIZfIU0VFBZYvX47Vq1ejtrYWU6dORUZGBurq6rzOf/DgQcyaNQuVlZWoqanB9OnTMWfOHNTW1spaLwOuivwNtv6E25Z4m+xwGwxRTyDUScR9Y/aTO/XOV/+k5meVx4H2GHJJruLiYmRnZyMnJwcTJkxASUkJ4uLiUFZW5nX+kpISPPnkk/iHf/gH3HbbbXj66adx22234Y9//KOs9TLgqsDIVVsGW+PgfiLR9Ox7tPiM8jgg0kdzc7Pbo7XV8+e629raUFNTg/T0dLfp6enpOHz4sF/r6ejoQEtLCwYPHiyrfbzITGF6Bttg8CRBShDtYhteHKO+OQ88D4yL9JiuZZ8i2gWYoh0HSuNxZQ7hJy8gLET+5zSkozNvxMXFuU1fs2YNCgsL3aY1NTWhvb0dUVFRbtOjoqLQ0NDg1/qef/55XL16FXPnzpXVTgZcheh5+69gK7YimRX9hSqvW9UwXpXX1ZtoJ3bST2lpKf7whz/A6XTijjvuQElJCaZOnep1XqfTiSeeeAI1NTU4deoUli1bhpKSEm0brDDRjgWGXDK7+vp6OBwO13O73Xe+sdncM40kSR7TvNmxYwcKCwvxzjvvYPjw4bLax4CrACNWbfUOtmoFWX/XZ6bAyxM7dV3EUVpairS0NGzcuBEZGRk4fvw4Ro0a5TF/a2srhg0bhtWrV+OFF15QvD169S+iHQtEXY68lq93ExTncDjcAq43Q4cORWhoqEe1trGx0aOq21NFRQWys7Oxc+dO3HfffbLbxzG4QdB7rG0g9BpjOyv6C7eH3kRrT7D0/sIiMitcFCP3Io7Ro0fjxRdfRFZWFiIjPYcXyNUSbxPm4iqRjgVR3hO1WOHYouCEh4cjKSkJVVVVbtOrqqqQmprqc7kdO3ZgwYIFePPNN/HAAw8EtG5WcANgxOEIWnf6RguNXe01U2WXrKHrIo5Vq1a5TZdzEYc/Wltb3S4iaW5u9jqfSAFTBGb/iwaHKlBf8vPzMX/+fCQnJyMlJQWbNm1CXV0dcnNzAQAFBQU4f/48tm3bBqAz3GZlZeHFF1/ElClTXNXf/v37y/pCzoArg7/BFvAv3Jop2Bot0Ppi5KAr0p9nzX5SF4kSF3H4o6ioCGvXrlXs9dQk0rFAZHWZmZm4dOkS1q1bB6fTicTERFRWViI+Ph5A5zUB3e+Ju3HjRty4cQOLFy/G4sWLXdMfffRRlJeX+71eBtw+yAm1gDhVWy2CrVlCrTdGDrrkyQpVpkAv4vBXQUEB8vO/H0fY3NyMuLg4NCd83+eJVL0VKeSa/QufFY4vCk5eXh7y8vK8/l/P0Lp//35F1smA64XcUAtYJ9iaOdR6Y7Sgy5O69QRzEYccdrvd51XSoo41Fel4MDuGXBKNZQNu29hYdIRFKPJaIgxHUDPYWi3UejMr+gvDhFyylu4XcTz00EOu6VVVVXjwwQd1bBl1xy98RNqybMANlr8VW8D/cCtSsGWo9WSUai6rVtYj9yIOADh69CgA4MqVK/jmm29w9OhRhIeH4/bbb5e17pa4zv5NpOEJ3fF40A6ruCQSBlyZ1Ai2gPxwq8bJhKHWP0ao5opyUhepamXmk6/cizgAYPLkya5/19TU4M0330R8fDz++te/Kto2b/2K6MePWkQ6HojMjgHXD3JCLWCsYMtQGxgjhFyyFjkXcQCdF6Gpqbe+Reu/hojyhc8KzPxFkoyFAdcHuaEWYLC1GoZcIqAtrg23jGxxm+Zv/2LFY4hVXCJtWDbgNifYERouP8R6Y5Rgy1CrPJFP0KJUrXhCtxa5/YxWx5Aox4MVsIpLIrBswFWCEYItQ636RA655I4nXnUF2t9Y7Rjilz4i9YXo3QCjaYm3uR7+aB3VJivcjh75jSLhdlb0Fwy3GhL1vRb1ynainrQ4hng8aEdOAYhIDazg+iGQA1WPiq2oIcsqrFaFIgKAUbFNAOyK9D88hohIKQy4XgTzzVPrYMtQS30RYewh/yRL5M4KxwSHBHU68lp+3zOR4iwbcFvibAiNUOZPKKzWUhdWoMiKlOyP1D6GRPjCR0Tqs2zADZYevzrGUGsMDLliY1WJRGGFKi6RXhhwZQgk1ALBBVuGWlICq1ZkJKzimge/UJJeGHB7EWigBVittTpWcT2xWmVO06NOAuindzOIiNxYNuC2xbUhpL+yd0mzQqh9xPGZ3/Puav6Rii0RH0MuUeB4/BBRMCwbcJViplArJ7wG+3pWD79ERIA1/rLBYQqkBwbcABg51CodYgPVsx0MvOrjuEOi7/F4IDI3Btw+GPknc0UJs/7oaqvZgi7/zComVpSMgccPEQXKsgF3VGwTwgbYVXt9rQOtkcJsb7pvh9nCLhEREWnDsgFXDVqFWrOE2b484vjMFCGXVajvWWG8oRVxyBH1hX81Ia0x4AZBi0BrlTDri1mHLuiF4w5JaT+7+b8BuN+RxixfTrXCL35EymPA9ROrs/riCZPIWIxwzPILH5F5WTbgTo86iYiB+t6cnGFWHiOcMH3hMAWiwPDYIaJAWDbg6oGBNnhGDrlEVsPjlYj0woCrEoZZ9fCkSURERL1hwA0SgywRERF5c+S1fL2bYFkMuL1geBWXEau4HEtIRFbGW4WRlkL6nsWcfnbzf+MRx2e9PojMSKlf5zO6lnib3k0gIrKE0tJSJCQkICIiAklJSaiuru51/gMHDiApKQkRERG45ZZb8Oqrr8pep2UDLhkfv4QYU+uoNr2bQBricUpkbRUVFVi+fDlWr16N2tpaTJ06FRkZGairq/M6/9mzZ3H//fdj6tSpqK2txW9+8xssW7YMu3fvlrVeBlwiIiIiUkVxcTGys7ORk5ODCRMmoKSkBHFxcSgrK/M6/6uvvopRo0ahpKQEEyZMQE5ODh577DE899xzstZruTG4ktQ5/ufKlQ6dW0JKuHblut5NkOXG1Va9mwAA6Pjumq7rb78mxvHX3qbOeMD2ts73t6u/MSM5fWmwx6max43ex0IXUY4Jtal1zImqublZ9dcOtJ+5IbUBAXzsbkhtbuvvYrfbYbfb3aa1tbWhpqYGq1atcpuenp6Ow4cPe339P//5z0hPT3ebNnv2bGzevBnXr19Hv37+/YaB5QJuS0sLAGDa3RyHaA779G6ATEZrLwWjpaUFkZGRejdDFfL60mA/9zxuyJgiX1+t+jrk9jPh4eGIjo7G/oZtAa9z4MCBiIuLc5u2Zs0aFBYWuk1rampCe3s7oqKi3KZHRUWhoaHB62s3NDR4nf/GjRtoampCTEyMX220XMCNjY1FfX09br75Zths5rzIpLm5GXFxcaivr4fD4dC7OYrj9hmf2bdRkiS0tLQgNjZW76aoxkx9qdk/j9w+Y/O1fYH2MxERETh79iza2gK/HkKSJI/jvmf1true83pbvq/5vU3vjeUCbkhICEaOHKl3MzThcDhMebB34fYZn5m30ayV2y5m7EvN/HkEuH1G5237Au1nIiIiEBERoUSzejV06FCEhoZ6VGsbGxs9qrRdoqOjvc4fFhaGIUOG+L1uXmRGRERERIoLDw9HUlISqqqq3KZXVVUhNTXV6zIpKSke8+/btw/Jycl+j78FGHCJiIiISCX5+fl47bXXsGXLFpw4cQKPP/446urqkJubCwAoKChAVlaWa/7c3Fx8/fXXyM/Px4kTJ7BlyxZs3rwZK1askLVeyw1RsAK73Y41a9b0Oh7GyLh9xmeFbSTjMPvnkdtnbEbfvszMTFy6dAnr1q2D0+lEYmIiKisrER8fDwBwOp1u98RNSEhAZWUlHn/8cWzYsAGxsbF46aWX8PDDD8tar00y831siIiIiMhyOESBiIiIiEyFAZeIiIiITIUBl4iIiIhMhQGXiIiIiEyFAdegSktLkZCQgIiICCQlJaG6utrnvHv27MGsWbMwbNgwOBwOpKSkYO/evRq2Vj4529fdRx99hLCwMEyaNEndBgZJ7va1trZi9erViI+Ph91ux5gxY7BlyxaNWhsYudu4fft2TJw4ETfddBNiYmKwcOFCXLp0SaPWklkVFhbCZrO5PaKjo3td5sCBA0hKSkJERARuueUWvPrqqxq1Vr7Ro0d7bJ/NZsPixYu9zr9//36v83/xxRcat9y7gwcPYs6cOYiNjYXNZsPbb7/t9v+SJKGwsBCxsbHo378/pk2bhs8//7zP1929ezduv/122O123H777XjrrbdU2oLe9bZ9169fx69//Wv88Ic/xIABAxAbG4usrCxcuHCh19csLy/3uk+vXbum8taIjQHXgCoqKrB8+XKsXr0atbW1mDp1KjIyMtxus9HdwYMHMWvWLFRWVqKmpgbTp0/HnDlzUFtbq3HL/SN3+7pcvnwZWVlZmDlzpkYtDUwg2zd37lz86U9/wubNm/Hll19ix44dGD9+vIatlkfuNh46dAhZWVnIzs7G559/jp07d+LTTz9FTk6Oxi0nM7rjjjvgdDpdj2PHjvmc9+zZs7j//vsxdepU1NbW4je/+Q2WLVuG3bt3a9hi/3366adu29Z1g/x//Md/7HW5L7/80m252267TYvm9unq1auYOHEiXnnlFa///+yzz6K4uBivvPIKPv30U0RHR2PWrFloaWnx+Zp//vOfkZmZifnz5+O///u/MX/+fMydOxeffPKJWpvhU2/b9+233+Kzzz7D//2//xefffYZ9uzZg5MnT+KnP/1pn6/rcDjc9qfT6dTkl8qEJpHh3HXXXVJubq7btPHjx0urVq3y+zVuv/12ae3atUo3TRGBbl9mZqb029/+VlqzZo00ceJEFVsYHLnb995770mRkZHSpUuXtGieIuRu4x/+8AfplltucZv20ksvSSNHjlStjWQNcvuDJ598Uho/frzbtEWLFklTpkxRuGXq+NWvfiWNGTNG6ujo8Pr/H374oQRA+p//+R9tGxYAANJbb73let7R0SFFR0dLzzzzjGvatWvXpMjISOnVV1/1+Tpz586V/tf/+l9u02bPni39/Oc/V7zNcvTcPm/+3//7fxIA6euvv/Y5z9atW6XIyEhlG2cCrOAaTFtbG2pqapCenu42PT09HYcPH/brNTo6OtDS0oLBgwer0cSgBLp9W7duxenTp7FmzRq1mxiUQLbv3XffRXJyMp599lmMGDECY8eOxYoVK/Ddd99p0WTZAtnG1NRUnDt3DpWVlZAkCRcvXsSuXbvwwAMPaNFkMrlTp04hNjYWCQkJ+PnPf44zZ874nPfPf/6zx2d39uzZOHLkCK5fv652U4PS1taGN954A4899hhsNluv806ePBkxMTGYOXMmPvzwQ41aGJyzZ8+ioaHBbf/Y7Xbce++9vZ4ffO1Tf8+Zerp8+TJsNhsGDRrU63xXrlxBfHw8Ro4ciZ/85CfC/oVWSwy4BtPU1IT29nZERUW5TY+KikJDQ4Nfr/H888/j6tWrmDt3rhpNDEog23fq1CmsWrUK27dvR1iY2D/OF8j2nTlzBocOHcJf/vIXvPXWWygpKcGuXbt8jrHTWyDbmJqaiu3btyMzMxPh4eGIjo7GoEGD8PLLL2vRZDKxu+++G9u2bcPevXvxr//6r2hoaEBqaqrP8d0NDQ1eP7s3btxAU1OTFk0O2Ntvv42///3vWLBggc95YmJisGnTJuzevRt79uzBuHHjMHPmTBw8eFC7hgaoq/+Qe/7ztU/9PWfq5dq1a1i1ahXmzZsHh8Phc77x48ejvLwc7777Lnbs2IGIiAikpaXh1KlTGrZWPGKnAfKp57dzSZL6/MYOADt27EBhYSHeeecdDB8+XK3mBc3f7Wtvb8e8efOwdu1ajB07VqvmBU3O/uvo6IDNZsP27dsRGRkJACguLsYjjzyCDRs2oH///qq3NxBytvH48eNYtmwZnnrqKcyePRtOpxMrV65Ebm4uNm/erEVzyaQyMjJc//7hD3+IlJQUjBkzBv/2b/+G/Px8r8t4++x6my6azZs3IyMjA7GxsT7nGTduHMaNG+d6npKSgvr6ejz33HO45557tGhm0AI5/wV6ztTL9evX8fOf/xwdHR0oLS3tdd4pU6ZgypQprudpaWn40Y9+hJdffhkvvfSS2k0VFgOuwQwdOhShoaEe3zwbGxs9vqH2VFFRgezsbOzcuRP33Xefms0MmNzta2lpwZEjR1BbW4slS5YA6AyEkiQhLCwM+/btw4wZMzRpuz8C2X8xMTEYMWKEK9wCwIQJEyBJEs6dOyfMxSFdAtnGoqIipKWlYeXKlQCAO++8EwMGDMDUqVOxfv16xMTEqN5usoYBAwbghz/8oc/qVnR0tNfPblhYGIYMGaJFEwPy9ddf47/+67+wZ88e2ctOmTIFb7zxhgqtUlbX3S8aGhrc+oS+zn++9mlf50y9XL9+HXPnzsXZs2fxwQcf9Fq99SYkJAT/8A//YPkKLocoGEx4eDiSkpJcV8p2qaqqQmpqqs/lduzYgQULFuDNN98Uelyj3O1zOBw4duwYjh496nrk5uZi3LhxOHr0KO6++26tmu6XQPZfWloaLly4gCtXrrimnTx5EiEhIRg5cqSq7Q1EINv47bffIiTEvTsKDQ0F8H31jEgJra2tOHHihM8vTSkpKR6f3X379iE5ORn9+vXTookB2bp1K4YPHx5Q/15bW2uIL5EJCQmIjo522z9tbW04cOBAr+c/X/u0t2X00hVuT506hf/6r/8K6EuVJEk4evSoIfapqvS5to2C8e///u9Sv379pM2bN0vHjx+Xli9fLg0YMED661//KkmSJK1atUqaP3++a/4333xTCgsLkzZs2CA5nU7X4+9//7tem9ArudvXk+h3UZC7fS0tLdLIkSOlRx55RPr888+lAwcOSLfddpuUk5Oj1yb0Se42bt26VQoLC5NKS0ul06dPS4cOHZKSk5Olu+66S69NIJN44oknpP3790tnzpyRPv74Y+knP/mJdPPNN/v8LJ45c0a66aabpMcff1w6fvy4tHnzZqlfv37Srl279NqEPrW3t0ujRo2Sfv3rX3v8X8/te+GFF6S33npLOnnypPSXv/xFWrVqlQRA2r17t5ZN9qmlpUWqra2VamtrJQBScXGxVFtb67qLwDPPPCNFRkZKe/bskY4dOyb94he/kGJiYqTm5mbXa8yfP9/tji0fffSRFBoaKj3zzDPSiRMnpGeeeUYKCwuTPv74Y6G27/r169JPf/pTaeTIkdLRo0fdztetra0+t6+wsFB6//33pdOnT0u1tbXSwoULpbCwMOmTTz7RfPtEwoBrUBs2bJDi4+Ol8PBw6Uc/+pF04MAB1/89+uij0r333ut6fu+990oAPB6PPvqo9g33k5zt60n0gCtJ8rfvxIkT0n333Sf1799fGjlypJSfny99++23GrdaHrnb+NJLL0m333671L9/fykmJkb6p3/6J+ncuXMat5rMJjMzU4qJiZH69esnxcbGSv/7f/9v6fPPP3f9v7fP4v79+6XJkydL4eHh0ujRo6WysjKNWy3P3r17JQDSl19+6fF/Pbfv97//vTRmzBgpIiJC+sEPfiD9+Mc/lv7zP/9Tw9b2rus2Zr7OVx0dHdKaNWuk6OhoyW63S/fcc4907Ngxt9e49957Pc5vO3fulMaNGyf169dPGj9+vG6BvrftO3v2rNf/AyB9+OGHrtfouX3Lly+XRo0aJYWHh0vDhg2T0tPTpcOHD2u/cYKxSRL//kdERERE5sExuERERERkKgy4RERERGQqDLhEREREZCoMuERERERkKgy4RERERGQqDLhEREREZCoMuERERERkKgy4RERERGQqDLhEREREZCoMuERERERkKgy4RERERGQqDLhEPowePRolJSVu0yZNmoTCwkJd2kNEFKxp06ZhyZIlWLJkCQYNGoQhQ4bgt7/9LSRJ0rtpRIpiwCUiIrKQf/u3f0NYWBg++eQTvPTSS3jhhRfw2muv6d0sIkWF6d0AIiIi0k5cXBxeeOEF2Gw2jBs3DseOHcMLL7yAX/7yl3o3jUgxrOASERFZyJQpU2Cz2VzPU1JScOrUKbS3t+vYKiJlMeAS+RASEuIxLu369es6tYaIiIj8xYBL5MOwYcPgdDpdz5ubm3H27FkdW0REFLyPP/7Y4/ltt92G0NBQnVpEpDwGXCIfZsyYgddffx3V1dX4y1/+gkcffZQnACIyvPr6euTn5+PLL7/Ejh078PLLL+NXv/qV3s0iUhQvMiPyoaCgAGfOnMFPfvITREZG4ne/+x0ruERkeFlZWfjuu+9w1113ITQ0FEuXLsU///M/690sIkXZJN78joiIyBKmTZuGSZMmedzjm8hsOESBiIiIiEyFAZeIiIiITIVDFIiIiIjIVFjBJSIiIiJTYcAlIiIiIlNhwCUiIiIiU2HAJSIiIiJTYcAlIiIiIlNhwCUiIiIiU2HAJSIiIiJTYcAlIiIiIlP5/wBYeE1bukOvXAAAAABJRU5ErkJggg==",
                        "text/plain": [
                            "<Figure size 800x400 with 3 Axes>"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                }
            ],
            "source": [
                "# Look at the contours of the (2,2,0,0) mode using teuk_all_amps\n",
                "plt.figure(figsize=(8, 4))\n",
                "plt.subplot(121)\n",
                "cb = plt.contourf(\n",
                "    u_grid.reshape(uv.size, wv.size),\n",
                "    w_grid.reshape(uv.size, wv.size),\n",
                "    np.abs(teuk_all_amps[:, :, kerr_amp_spline.special_index_map[(2, 2, 0, 0)]]),\n",
                ")\n",
                "plt.xlabel(\"u\")\n",
                "plt.ylabel(\"w\")\n",
                "plt.subplot(122)\n",
                "cb = plt.contourf(\n",
                "    p_grid.reshape(uv.size, wv.size),\n",
                "    e_grid.reshape(uv.size, wv.size),\n",
                "    np.abs(teuk_all_amps[:, :, kerr_amp_spline.special_index_map[(2, 2, 0, 0)]]),\n",
                ")\n",
                "plt.colorbar(cb)\n",
                "plt.xlabel(\"p\")\n",
                "plt.ylabel(\"e\")\n",
                "plt.show()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 8,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAGCCAYAAAALybB5AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWfxJREFUeJzt3X1UVPedP/D3wMhAVDA+ISgisfEpJOrCWoFgjCYYkpq2m1ZSuxItZOWgGEI0lbW/iNYt0SaEmATUjUpt1HJ8StJTjGHXKqg1Gwm2RowxagOYQYqbAD4EFO7vD5eJwzzP3Od5v86Zc5zLvXO/V7jf+76f+d57DYIgCCAiIiIi0okApRtARERERCQmBlwiIiIi0hUGXCIiIiLSFQZcIiIiItIVBlwiIiIi0hUGXCIiIiLSFQZcIiIiItIVBlwiIiIi0hUGXCIiIiLSFQZcIiIiItIVBlwiIiIilaqqqsLs2bMRGRkJg8GAd999V9L1FRQUwGAwWL2GDRsm6TqlwIBLREREpFLXrl3DxIkT8eabb8q2zvvuuw9ms9nyOnXqlGzrFotR6QYQERERkX2pqalITU11+PPOzk786le/wvbt2/HNN98gNjYWa9euxfTp071ep9Fo1GTV9k6s4BIRERFp1IIFC3D06FH84Q9/wN/+9jf89Kc/xWOPPYZz5855/Znnzp1DZGQkYmJi8PTTT+PChQsitlgeBkEQBKUbQURERETOGQwG7Nu3Dz/60Y8AAOfPn8e9996LxsZGREZGWuZ75JFHMGXKFPzmN7/xeB379+/H9evXMWbMGFy+fBlr1qzBZ599htOnT2PQoEFibYrkWMElIiIi0qBPPvkEgiBgzJgx6Nevn+V1+PBhnD9/HgDw97//3eaisd6vxYsXWz4zNTUVTz31FO6//3488sgj+NOf/gQA+N3vfqfINnqLY3CJiIiINKi7uxuBgYGoqalBYGCg1c/69esHABg+fDjOnDnj9HPuvvtuhz/r27cv7r//fp+GPCiBAZeIiIhIgyZPnoyuri40NzcjOTnZ7jx9+vTBuHHjvF5HR0cHzpw54/Dz1YoBl4iIiEilrl69ii+++MLy/uLFizh58iQGDhyIMWPG4Oc//znS09Px6quvYvLkyWhpacHBgwdx//334/HHH/d4fUuXLsXs2bMxcuRINDc3Y82aNWhra8Mzzzwj5mZJjheZEREREanUoUOH8PDDD9tMf+aZZ1BWVoabN29izZo12LZtGy5duoRBgwYhISEBq1atwv333+/x+p5++mlUVVWhpaUFQ4YMwdSpU/HrX/8aEyZMEGNzZMOAS0RERESiKywsxN69e/HZZ58hJCQEiYmJWLt2LcaOHet0ue3bt2PdunU4d+4cwsLC8Nhjj+GVV17x6C4OvIsCEREREYnu8OHDWLRoEY4fP47KykrcunULKSkpuHbtmsNljhw5gvT0dGRkZOD06dPYtWsXPv74Y2RmZnq0blZwiYiIiEhy//jHPzB06FAcPnwY06ZNszvPK6+8gtLSUsttzgDgjTfewLp169DQ0OD2uvzuIrPu7m589dVX6N+/PwwGg9LNISIdEgQB7e3tiIyMRECAPr8oY19KpCxf+plvv/0WnZ2dPq27935vMplgMpmcLtfa2goAGDhwoMN5EhMTsWLFClRUVCA1NRXNzc3YvXs3nnjiCY8b6VcaGhoEAHzxxRdfkr8aGhqU7vIkw76UL77U8fK0n7lx44YweEiAT+vs16+fzbSVK1c6XW93d7cwe/Zs4cEHH3TZxl27dgn9+vUTjEajAEB48sknhc7OTo+20+8quP379wcATA9fAGNAkKzrvjEhwuU8IXVmUT7H13WIsR53uNsWV+RoqzNibQfpw63uThy6vNXS3+hRz7Y1NDQgNDRU4dYQ6ctT4150OY+3/UxnZyda/tGNA8eHoW8/z79huna1G7OmNtns+66qt4sXL8bf/vY3HDlyxOl8dXV1WLJkCV566SXMmjULZrMZy5YtQ1ZWFjZv3ux2O/0u4PaU1I0BQbIHXKMx2OnPQz69BLjRJlef43J5N7fb1/W4tQ4Rfgc3Yocr/ocs998SaYOev7rv2bbQ0FAGXCIfpY5YYvXek2OKt/1M334B6Nff+yFUnuz7OTk5eP/991FVVYURI0Y4nbewsBBJSUlYtmwZAOCBBx5A3759kZycjDVr1iAiwr2CltK5gDx0I3a4T8uHfHpJlvW4w922EBER6UnvQKtXgiAgJycH+/btw6FDhxATE+NymevXr8NotI6nPY8hFjy4LwIDrkzkCIz+iP+vRESkdv4SaHtbtGgRduzYgffeew/9+/dHU1MTACAsLAwhISEAgPz8fFy6dAnbtm0DAMyePRvPPvssSktLLUMUcnNzMWXKFERGRrq9bgZclXCnmqmn6q0Y1NJOVqKJiKg3fw21dyotLQUATJ8+3Wr61q1bMX/+fACA2WxGfX295Wfz589He3s73nzzTbzwwgsYMGAAZsyYgbVr13q0bgZcGbgKYv4YkHzdZrWEWyIiIoCB1h53hhSUlZXZTMvJyUFOTo5P62bA9RN6q94SEREpiYFW3RhwFeaPwZPVWyIi0iKGWu1gwJWYGsKYnkK0FtpIRET6wECrXQy4CtJT8HSXP443JiIibWCg1Q8GXJ3TU4hWYxsZ2ImItI2hVp8YcCXkLJD5YzDyZZvVGG6JiEh79BJo93y2DmFhG5Vuhmox4KqcL8FOT9VbIiIib+kl1JL7GHAVwOqtZxjAiYjIUwy1/o0BVyJihDJWb9XdPn88USEiUisGWroTA67M/DEU+eM2ExGR9BhqyREGXFItNVdviYhIGQy15A4GXBl5UsnUy/AEb6u3DLdERNSDoZY8xYArAYYz/eOwCyIi6TDQkq8YcHWG1VsiItIihloSEwOuTOQankBERKQVDLUkFQZcHWH1Vh4cnkBE5D2GWpIDA64MWL11jz9vOxGRnjHUimt/43q0tbUp3QxVY8CVmFzVPj1Ub7VAz9tGRCQmhlpSEgMuqQKrt0RE2sdQS2rBgCuBkE8veRXYvA15rN4SEZGSGGxJbRhwJcKg5z6tVG/5OyUi+g5DLakZA67GqSl0qaktREQkPoZa0goGXJWQuoqp1iqpr+36emyQ1fu7z3b69HlERGSNoZa0KEDpBpA+eFO99SXcfj02yCbc9kyXAqvTpHYlJSWIiYlBcHAw4uLiUF1d7XT+7du3Y+LEibjrrrsQERGBBQsW4MqVKzK1lrQgdcQShlvSLAZcFdDDxWVykirEEmlVeXk5cnNzsWLFCtTW1iI5ORmpqamor6+3O/+RI0eQnp6OjIwMnD59Grt27cLHH3+MzMxMmVtOatMTahlsSes4RIF8Jmf1Volwy+otqV1RUREyMjIsAbW4uBgHDhxAaWkpCgsLbeY/fvw4Ro0ahSVLboeYmJgYLFy4EOvWrZO13aQODLOkR6zgapQ/hi53wy0rvORPOjs7UVNTg5SUFKvpKSkpOHbsmN1lEhMT0djYiIqKCgiCgMuXL2P37t144oknHK6no6MDbW1tVi/SNlZqSc9YwVWYP15c5k2blAqt/ngiQdrS0tKCrq4uhIeHW00PDw9HU1OT3WUSExOxfft2pKWl4dtvv8WtW7fw5JNP4o033nC4nsLCQqxatUrUtpP8GGi1b3/jeqWboAms4GqQmkKXHG1hRZbINYPBYPVeEASbaT3q6uqwZMkSvPTSS6ipqcEHH3yAixcvIisry+Hn5+fno7W11fJqaGgQtf0kLVZryd+wgqtjeqjeMtwSOTd48GAEBgbaVGubm5ttqro9CgsLkZSUhGXLlgEAHnjgAfTt2xfJyclYs2YNIiIibJYxmUwwmUzibwBJhoGW/BkruAryJoD6U/VW6XCrpv9rIkeCgoIQFxeHyspKq+mVlZVITEy0u8z169cREGDd/QcGBgK4XfklbWO1logVXJKRXBVlPuyB/E1eXh7mzZuH+Ph4JCQkYNOmTaivr7cMOcjPz8elS5ewbds2AMDs2bPx7LPPorS0FLNmzYLZbEZubi6mTJmCyMhIJTeFvMRAS2SNAVenpA6TrN4SqUdaWhquXLmC1atXw2w2IzY2FhUVFYiOjgYAmM1mq3vizp8/H+3t7XjzzTfxwgsvYMCAAZgxYwbWrl2r1CaQlxhsiexjwPVCT3j0JQRpfXiCpzzZXqXDLZEWZWdnIzs72+7PysrKbKbl5OQgJydH4laRFBhqSSsKCwuxd+9efPbZZwgJCUFiYiLWrl2LsWPHurX80aNH8dBDDyE2NhYnT570aN0cg+sDNV7EBWi7equGcKvlEwki0i+OrSWtOXz4MBYtWoTjx4+jsrISt27dQkpKCq5du+Zy2dbWVqSnp2PmzJlerZsVXA8pFWq1HLrUeiJARKQFDLWkVR988IHV+61bt2Lo0KGoqanBtGnTnC67cOFCzJ07F4GBgXj33Xc9XjcDrgKkDHys3vpGyycSRKQfDLWkdr2fZujOrQRbW1sBAAMHDnQ639atW3H+/Hm88847WLNmjVftY8D1ACuRnnP3/0yscMs7KBCRljHYklz+8M33YbrVx+PlOq7eBLAPUVFRVtNXrlyJgoICh8sJgoC8vDw8+OCDiI2NdTjfuXPnsHz5clRXV8No9D6mMuC6Sclwq5aqolraIRW9bx8RqReDLblDTY/pbWhoQGhoqOW9q+rt4sWL8be//Q1HjhxxOE9XVxfmzp2LVatWYcyYMT61jwFXZloenuApuau3RERaw2BLWhUaGmoVcJ3JycnB+++/j6qqKowYMcLhfO3t7Thx4gRqa2uxePFiAEB3dzcEQYDRaMSHH36IGTNmuLVOBlw3OApqclT8WFUkItIfBlvyB4IgICcnB/v27cOhQ4cQExPjdP7Q0FCcOnXKalpJSQkOHjyI3bt3u1z+Tgy45BZPg7YS1Vtfxt/yRIKI5MBgS/5k0aJF2LFjB9577z30798fTU1NAICwsDCEhIQAsH7SYkBAgM343KFDhyI4ONjpuF17GHBdEPNrf38anuAODk0gIn/AUEv+qrS0FAAwffp0q+lbt27F/PnzAdg+aVEsDLhOKB0a1VJVlKp6qxZq+X8mIn1hsCV/JwiCy3nsPWnxTgUFBU7vzuAIA64DrkKamkKR1gIlIH71lrcHIyK1YLAlUh4Drh1SBEZPP1MtAZrVWyIi9zDYEqkHA24v7gQ0hiLfcOwtEekJgy2R+jDg/h+tXgCmpoqpUm3xdngCT1SIyBcMtkTq5bcB98aECBiNwUo3wy61BC8p2sHqLRFpHYMtkfr5bcD1llrCp9qoqZLsDv4eichTDLZE2hGgdANKSkoQExOD4OBgxMXFobq62un827dvx8SJE3HXXXchIiICCxYswJUrV2Rqree0OjxBK9Vb3j2BiOTAcEtK29+4XukmaIqiFdzy8nLk5uaipKQESUlJ2LhxI1JTU1FXV4eRI0fazH/kyBGkp6fjtddew+zZs3Hp0iVkZWUhMzMT+/btk7y9fDSvfazeis/Z/6kW2k+kFwy2RNqkaMAtKipCRkYGMjMzAQDFxcU4cOAASktLUVhYaDP/8ePHMWrUKCxZcrvDiYmJwcKFC7Fu3TpZ260GrN7qiye/z97zMvASiY/BlkjbFBui0NnZiZqaGqSkpFhNT0lJwbFjx+wuk5iYiMbGRlRUVEAQBFy+fBm7d+/GE0884XA9HR0daGtrs3rJRWuVTa3xZniC2sLgjdjhPv+diPEZRHRb6oglDLdEOqBYwG1paUFXVxfCw8OtpoeHh6OpqcnuMomJidi+fTvS0tIQFBSEYcOGYcCAAXjjjTccrqewsBBhYWGWV1RUlFft5fAE+9wJVqze2pIilDLoEnmPwZZIXxS/yMxgMFi9FwTBZlqPuro6LFmyBC+99BJqamrwwQcf4OLFi8jKynL4+fn5+WhtbbW8GhoaRG2/I7y4TFpard7KEUIZconcx2BLpE+KjcEdPHgwAgMDbaq1zc3NNlXdHoWFhUhKSsKyZcsAAA888AD69u2L5ORkrFmzBhERETbLmEwmmEwmn9qqhmCkRqzeuk/u0NmzPv7tEjnGYEukX4pVcIOCghAXF4fKykqr6ZWVlUhMTLS7zPXr1xEQYN3kwMBAALcrv1qmhiCihjZIRcltU7Kiymqu//Dklovz58+HwWCwed13330ytlg5rNoS6Z+iQxTy8vLw9ttvY8uWLThz5gyef/551NfXW4Yc5OfnIz093TL/7NmzsXfvXpSWluLChQs4evQolixZgilTpiAyMlKSNnoTjLQ6PMETSlZvtXTvWzX8vtTQBpJWzy0XV6xYgdraWiQnJyM1NRX19fV253/99ddhNpstr4aGBgwcOBA//elPZW65vBhsifyHorcJS0tLw5UrV7B69WqYzWbExsaioqIC0dHRAACz2WzVQc+fPx/t7e1488038cILL2DAgAGYMWMG1q5dq9Qm6Aart+JSW6i8ETtc179jf+fpLRd7Lrrt8e677+Lrr7/GggULZGuz3BhsifyL4o/qzc7ORnZ2tt2flZWV2UzLyclBTk6OxK26Ta5AoLXgwbG3zqkt3PZgyNWnnlsuLl++3Gq6s1su9rZ582Y88sgjluKCPR0dHejo6LC8l/OWi75gsCXyT4oHXL3xh+EJSvJ0eILcgU7tvyOGXP3x5paLdzKbzdi/fz927NjhdL7CwkKsWrXKp7bKicGW9ISP6fWc4rcJUyt/CgFib6u/Vm/VHm5J3zy55eKdysrKMGDAAPzoRz9yOp9St1z0xvTHOGyNyN+xgmuHt4HPm4Dj7rrUEp6UbIeaq7dq+f24g1VcffHmlos9BEHAli1bMG/ePAQFOT8xFeOWi1JjsCWiHqzg+jlWb32npXDbQ4ttJvu8ueVij8OHD+OLL75ARkaGlE2U3PTH1jLcEpEVVnB7YWXLMS2FIrl+j1r6PyH9ysvLw7x58xAfH4+EhARs2rTJ5paLly5dwrZt26yW27x5M77//e8jNjZWiWaLgsGWiOxhwL2DL6HIH4YnuCJl9VaN977Vyu/FEQ5V0A9Pb7kIAK2trdizZw9ef/11JZrsMwZbInLGbwNuSJ0ZxgD/+zr9Tp6EGy2FOTlCm5b+P5xhyNUPT2+5GBYWhuvXr0vcKmkw3BKRK34bcMWkl7DjC3+q3vL3TaQMBlsichcvMlOI3oYnqIXU1Ug9/j70uE2kPwy3ROQJVnB9pNVwoNfhCVLi/wOR/CwPbPBg/+OwGyJiBVfFtBKo1DI8QcqDmlZ+F97S+/aRNvFpZETkLVZwfeBtKFC6usDqLRGpWc9whJD/e89+iIg8xQou+YTVW/3wl+0kdeNYWyJr+xvXK90ETWLA9ZLUYYBhQx34eyCSz53hVulvuohI2zhEQWZKd9piDk9QS/VWKgy3RPJIHbHE4f7m6X6odB9LROrACq4XGHzURYoDmr/+jv11u0k50x9by787IhIdA64KSdXZa6V6S0T+geNtifStsLAQ//zP/4z+/ftj6NCh+NGPfoSzZ8+6XO7w4cOIi4tDcHAw7rnnHmzYsMHjdTPgesiX8Mmvztyj9MVl/l5N8vftJyIicRw+fBiLFi3C8ePHUVlZiVu3biElJQXXrl1zuMzFixfx+OOPIzk5GbW1tfj3f/93LFmyBHv27PFo3RyD6wEe+PWPv2Mi5fWcuHL8LZG2ffDBB1bvt27diqFDh6KmpgbTpk2zu8yGDRswcuRIFBcXAwDGjx+PEydO4JVXXsFTTz3l9rpZwVUZfx+eoHT1loiIiFxra2uzenV0dLhcprW1FQAwcOBAh/P85S9/QUpKitW0WbNm4cSJE7h586bb7WMF102+Bk+GMfVj9fY7N2KH82+WiEjHDjbei8C7TB4v13X9dpCNioqymr5y5UoUFBQ4XE4QBOTl5eHBBx9EbGysw/mampoQHh5uNS08PBy3bt1CS0sLIiIi3GonA66KqCFgaeXiMrHDlxr+74n8SeqIJQD3OyLNamhoQGhoqOW9yeQ8LC9evBh/+9vfcOTIEZefbTAYrN4LgmB3ujMMuG7QevjRSiVOqXvfav33S6RHHH9LpG6hoaFWAdeZnJwcvP/++6iqqsKIESOczjts2DA0NTVZTWtubobRaMSgQYPcbh/H4LogRvhhxysu/n/Kg8GfpOTo74v7N9F3tP6YXkEQsHjxYuzduxcHDx5ETEyMy2USEhJQWVlpNe3DDz9EfHw8+vTp4/a6GXBVQg1hQisXl4lJDf/vREREerRo0SK888472LFjB/r374+mpiY0NTXhxo0blnny8/ORnp5ueZ+VlYUvv/wSeXl5OHPmDLZs2YLNmzdj6dKlHq2bAdcJPYQfvVVDxNwePfx+ifSI+yaRPpSWlqK1tRXTp09HRESE5VVeXm6Zx2w2o76+3vI+JiYGFRUVOHToECZNmoRf//rXWL9+vUe3CAM4BtchsTpYdwKZGjpzrVxcRvLi3RRIC/g3SqROPReHOVNWVmYz7aGHHsInn3zi07pZwbVDDYFTDFrp9N0dnsDqLREREbmDFdxeGHxs6a16K/Xv2N3/L6XGHROplVZOyolI/Rhw7yB28NHL8AQp6SXkeXMS0LOMXv4PiJTAUExE9vjtEIUbEyJwI3a41UtP9Nbpi7U9Yv+evx4b5HOFW4zPkJLe9g29KikpQUxMDIKDgxEXF4fq6mqn83d0dGDFihWIjo6GyWTC6NGjsWXLFpla6xj/3ohIDKzg+jl/ujWYFOFW7M9jNZe8UV5ejtzcXJSUlCApKQkbN25Eamoq6urqMHLkSLvLzJkzB5cvX8bmzZvxve99D83Nzbh165bMLXeO33IQkbcYcCWileEJWqC2arSUoZ8hl7xRVFSEjIwMZGZmAgCKi4tx4MABlJaWorCw0Gb+Dz74AIcPH8aFCxcwcOBAAMCoUaPkbDIRkaT8doiCnokVCNX8tbmnxDqZkOP/RE//7yS9zs5O1NTUICUlxWp6SkoKjh07ZneZ999/H/Hx8Vi3bh2GDx+OMWPGYOnSpVY3X++to6MDbW1tVi9vpY5Y4nIed/YDtZ38EpF6sIIrAa1Ub7VwcZmaDmByBk9WcsldLS0t6OrqQnh4uNX08PBwm+e597hw4QKOHDmC4OBg7Nu3Dy0tLcjOzsb//u//OhyHW1hYiFWrVonSZnt9T8inl1TRLxKRPrCCS3bpqYro60FT7ReByYHBQ/0MBoPVe0EQbKb16O7uhsFgwPbt2zFlyhQ8/vjjKCoqQllZmcMqbn5+PlpbWy2vhoYG0beBiL6zv3G90k3QNAZcnXG34ukv1Vsxwq1S/D1Uk3sGDx6MwMBAm2ptc3OzTVW3R0REBIYPH46wsDDLtPHjx0MQBDQ2NtpdxmQyITQ01OolFQ5PICJfMeCKTCvDE8g1NQRMNbSB1C0oKAhxcXGorKy0ml5ZWYnExES7yyQlJeGrr77C1atXLdM+//xzBAQEYMSIEZK2l4hIDgy4ZEMNoUrp6q0a/g+I3JWXl4e3334bW7ZswZkzZ/D888+jvr4eWVlZAG4PL0hPT7fMP3fuXAwaNAgLFixAXV0dqqqqsGzZMvziF79ASEiIUptBRCQaXmQmMymrt3oanqAktYVbXnBGrqSlpeHKlStYvXo1zGYzYmNjUVFRgejoaACA2WxGfX29Zf5+/fqhsrISOTk5iI+Px6BBgzBnzhysWbNGqU2w9EtK7H+u+kQOhyDSHgZcEemhE1RbuPOWtyFeL9tP/ic7OxvZ2dl2f1ZWVmYzbdy4cTbDGpTi6R0U5H6yYc98eujjifwFhyjohBaqt+7y9SCihW30lBqCtx7/X8k/eft4dj0+1p1IrxhwZeTvHaPav2ZXQ4gkImti9xti9MP+3pcTaQEDrkj08NWVGgKeUtVbNWw7kb+T8vZgYldfGXKJ1I0BVya8uEy91VuthFuttJPIU1L3TVJ9PkMukXox4IqA1VtxKFG9VcN2E/mT1BFLZF2fVsMzEfmGAddPsBO2pcVwq8U2E92pd1+khwIBEakPA64M1DA8QUlyDE9ggCfSNnsnb737Dk/7O7n6BfY/JLb9jeuVboLmMeD6SAsB0xU1VAXl/n9UwzZrEQ/kpBVy/61y3yBSFwZciamh09P7xWWebh/DLZG+qaHfJSJlMeD6QKyq49djg7wKXWKsXw1hz5ft8Mdwq4dtIOqhpzCqp20h0jrFA25JSQliYmIQHByMuLg4VFdXO52/o6MDK1asQHR0NEwmE0aPHo0tW7bI1FppSRFc2OESkR5o4XaIRKQeRiVXXl5ejtzcXJSUlCApKQkbN25Eamoq6urqMHLkSLvLzJkzB5cvX8bmzZvxve99D83Nzbh165bMLXevs3Wno1WyGif1ut0ZnsDqLRGJtW+qIdzeiB2ui2sziLRO0YBbVFSEjIwMZGZmAgCKi4tx4MABlJaWorCw0Gb+Dz74AIcPH8aFCxcwcOBAAMCoUaPkbLJqsAP1DMMtkfqEfHrJYShV88NhiEj9FBui0NnZiZqaGqSkpFhNT0lJwbFjx+wu8/777yM+Ph7r1q3D8OHDMWbMGCxduhQ3btxwuJ6Ojg60tbVZvXylleqtni8uU0OlRkkM7ORPxOpz5aKmthD5K8UquC0tLejq6kJ4eLjV9PDwcDQ1Ndld5sKFCzhy5AiCg4Oxb98+tLS0IDs7G//7v//rcBxuYWEhVq1aJXr7faV0QFF6/YB8VWg1bCsRSYeBkoh6U/wiM4PBYPVeEASbaT26u7thMBiwfft2TJkyBY8//jiKiopQVlbmsIqbn5+P1tZWy6uhocGn9qqhkqCGNijJk21juCVSLz33U0SkLMUC7uDBgxEYGGhTrW1ubrap6vaIiIjA8OHDERYWZpk2fvx4CIKAxsZGu8uYTCaEhoZavZSmdOjS+sVlRKRN0x9bazPN1/5IrSFZre0i8heKBdygoCDExcWhsrLSanplZSUSExPtLpOUlISvvvoKV69etUz7/PPPERAQgBEjRkjaXkD6UMZg6Bqrt0T6d2df6O99HhF5R9EhCnl5eXj77bexZcsWnDlzBs8//zzq6+uRlZUF4PbwgvT0dMv8c+fOxaBBg7BgwQLU1dWhqqoKy5Ytwy9+8QuEhIQotRlWXAUwOUKXszaoIfTJccBSw3YSkbRYJSU92t+4Xukm6IKiATctLQ3FxcVYvXo1Jk2ahKqqKlRUVCA6OhoAYDabUV9fb5m/X79+qKysxDfffIP4+Hj8/Oc/x+zZs7F+vfR/DKwiuEfKuyfwYGZNqRDP3wOJgX0qkX+oqqrC7NmzERkZCYPBgHfffdflMmI81EvR++ACQHZ2NrKzs+3+rKyszGbauHHjbIY1SE2sJ+j4Gkh8vbjMX6qa/rKdRP5MCydafOgDEXDt2jVMnDgRCxYswFNPPeXWMmI81EvxgEv6IeUYYncPZgy3RNrgzj7NcEikfampqUhNTXV7frEe6qX4bcLUTq7qLZ/aQ0T+xl6/6E5fqIXqLZHe9X6IVkdHhyif681DvexhBdcJsaoHYlQV9TA8gdVbIiIidbj2ZSgCgoM9Xq77228BAFFRUVbTV65ciYKCAp/b5c1DvexhBdcBT8IYqwmsQBMpraSkBDExMQgODkZcXByqq6sdznvo0CEYDAab12effSZji/0Djw+kVw0NDVYP0srPzxflc715qJc9DLh2iDnuy52qohjhUO3VW2+xekvkWnl5OXJzc7FixQrU1tYiOTkZqampVnehsefs2bMwm82W17333itTi11z1A8zMBKpQ++HaJlMJlE+15uHetnDgNuLp+FWjs5WDxdaSLkNUobb9tHddl9EalJUVISMjAxkZmZi/PjxKC4uRlRUFEpLS50uN3ToUAwbNszyCgwMlKnFRET2ifVQLwbc/xPy6SXRQ5gaqopytEGq4QlKVWrcCbIMuqQWnZ2dqKmpQUpKitX0lJQUHDt2zOmykydPRkREBGbOnIk///nPTuft6OiwuahEbK76ElZvibTn6tWrOHnyJE6ePAkAuHjxIk6ePGn5hkmqh3r5bcANqTNbQq3UFz/5wteLy9RAK9Vbb0Kr0kFXDSdRpKyWlhZ0dXUhPDzcanp4eDiamprsLhMREYFNmzZhz5492Lt3L8aOHYuZM2eiqqrK4XoKCwsRFhZmefW+wMSV6Y+ttXrPv10i/3DixAlMnjwZkydPBnD7KbaTJ0/GSy+9BEC6h3rxLgoScbfzlvLiLFZv3cdqLGmdwWCwei8Igs20HmPHjsXYsWMt7xMSEtDQ0IBXXnkF06ZNs7tMfn4+8vLyLO/b2to8Drnu0sOwLCK6bfr06RAEweHPpXqol99WcH2llqqpWtrhiNqrt2JVYBmQSSmDBw9GYGCgTbW2ubnZpqrrzNSpU3Hu3DmHPzeZTDYXlchJ7X2dM1puO8lrf6NnVUpyjAFXAmJVTv21iuHOwUCscCsmhlxSQlBQEOLi4myqHZWVlUhMTHT7c2praxERESF284iIFMEhCl5Qy9m40rcG0/K9b6UKo+2ju9H/PM8bSV55eXmYN28e4uPjkZCQgE2bNqG+vh5ZWVkAbg8vuHTpErZt2wYAKC4uxqhRo3Dfffehs7MT77zzDvbs2YM9e/Yotg3O+hNP+tzefZ+W+yki8h4Drsg8CZZ673i9qUDLUb1lpZX0Ji0tDVeuXMHq1athNpsRGxuLiooKREdHA7C9iKOzsxNLly7FpUuXEBISgvvuuw9/+tOf8Pjjjyu1CRbefnPlqF/oma73/paIrDHgekiu6q0vwxN4dbJjcoRbVnFJCdnZ2cjOzrb7s94Xcbz44ot48cUXZWiVfWL3Ue583tdjgxhyifwIj8IecBVu5azeKj1MwlX71Vi9ZeWWSHvE7HeJyH8w4KoQq7f2MdwS0Z087RP03D8SkTUGXDepqYqg9uqtN6TcJiXCLQM1kfu86VMYVonIGQZcNygdKN2llg5fitubebttDJpE2tG775Ci71VLP0lE0mLAFYGnHaazaoU/3vtWqhMIhlsifWJIJSJXGHBdUFv11lF75Orwpbi4zBUezIi0bfpjay3/9nV/FqM/YJ9CpH9eBdyf//zn2LRpEz7//HOx26Mqcj1Rqwert+JRQ/VWDW0g0iq1FReISFu8Crj9+vVDUVERxo0bh8jISPzsZz/Dhg0b8Nlnn4ndPsVIFW59uUBL7R2+Wqq3DJakFdXV1fjXf/1XJCQk4NKl2/vP73//exw5ckThlsnHkz6RlVfSs/2N65Vugq54FXA3btyIzz77DF999RWKiooQFhaG119/Hffdd58unmWu9iDZm1qGJ3hKiv9nhlvp+eM3DVLYs2cPZs2ahZCQENTW1qKjowMA0N7ejt/85jcKt05+/LsiIjH5NAa3f//+uPvuu3H33XdjwIABMBqNGDZsmFhtU4S7oUvs6q2rzl1roVsMnv4fM9ySlqxZswYbNmzAf/7nf6JPnz6W6YmJifjkk08UbJny7PV3rN4SkSe8Cri//OUvMXXqVAwePBi/+tWv0NnZifz8fFy+fBm1tbVit1E2UoZbqaileutp9UVN9xUmUsLZs2cxbdo0m+mhoaH45ptv5G+Qn2EfQ6RvRm8W+u1vf4shQ4Zg5cqV+OEPf4jx48eL3S7J3ZgQAaMxWOlmAGD1Vgxqrd62j+5G//O8WQnZioiIwBdffIFRo0ZZTT9y5AjuueceZRqlUgyjROQpr468tbW1WLFiBf7nf/4H06ZNw7Bhw5CWlobS0lKcOXNG7DaqircdrRRP/1JLp6909Vat4ZbImYULF+K5557DRx99BIPBgK+++grbt2/H0qVLkZ2drXTzZMXxt0QkNq8quBMnTsTEiROxZMkSAMBf//pXFBcXY8mSJeju7kZXV5eojVQLtQRKuUkRzsXCcEta9eKLL6K1tRUPP/wwvv32W0ybNg0mkwlLly7F4sWLlW6eJHr3ofb6lt4nwP7a7xKRb7wKuMDtKu6hQ4dw6NAhVFdXo62tDZMmTcLDDz8sZvt0z9vhCVrt9Dn2lug7//Ef/4EVK1agrq4O3d3dmDBhAvr166d0s4iINM+rgHv33Xfj6tWrmDhxIqZPn45nn30W06ZNQ2hoqNjtUw1fgpeaK6C+UvKrRVZvSQ/uuusuxMfHK90MVeIJLxF5y6uA+/vf/173gfZOUnWyWqjeihnOxazeMtzepueTJ/I/vKCWiMTiVcD9wQ9+IHY7VMvXMKnnAMILQ4jIlemPrXX6c/YjRCQF3r/ICSW/HmP11jFWb4n0T46+jkMgiPSLAdcBqTs+rVctxGw/wy2R/9Hzt1tEntrfuF7pJugOA64dYoVbbztwNVRvxSTWuDqGWyLtctV/cfwtEYmJAbcXOUKkFqq3clVXtBra/ZUW/nZJ+9gvEJGvGHDvIGanKnZAVFOH70nIYfWWiIiI5Oa3Afeb7wXh67HWLzl4e2swOamtestwS+RaSUkJYmJiEBwcjLi4OFRXV7u13NGjR2E0GjFp0iRpG2gHvxEgIqn4bcDVEn+v3hKRc+Xl5cjNzcWKFStQW1uL5ORkpKamor6+3ulyra2tSE9Px8yZM2VqqX139hVq6u/EwiBPJD8GXAmIfXGZHrF6SySeoqIiZGRkIDMzE+PHj0dxcTGioqJQWlrqdLmFCxdi7ty5SEhIkKmlt/EOCkQkNQZcGXlzFi93NcPZgUfu6i3DLZFrnZ2dqKmpQUpKitX0lJQUHDt2zOFyW7duxfnz57Fy5Uq31tPR0YG2tjarlzv0WJElIvVjwBUZq7eu8YBHJJ6WlhZ0dXUhPDzcanp4eDiamprsLnPu3DksX74c27dvh9Ho3gMtCwsLERYWZnlFRUX53Pbe2DcQkVgYcGXC6q1nWL11jV/z0p0MBoPVe0EQbKYBQFdXF+bOnYtVq1ZhzJgxbn9+fn4+WltbLa+Ghga3l7XXl/X0Kf50ck/kj6qqqjB79mxERkbCYDDg3XffdTr/3r178eijj2LIkCEIDQ1FQkICDhw44PF6GXBVQKxAqIVQyAoNkbgGDx6MwMBAm2ptc3OzTVUXANrb23HixAksXrwYRqMRRqMRq1evxl//+lcYjUYcPHjQ7npMJhNCQ0OtXkRErly7dg0TJ07Em2++6db8VVVVePTRR1FRUYGamho8/PDDmD17Nmpraz1ar3vfTZFbHFXU5Kze9oTc/ueVOXdxFtZ5YZm28UpwdQoKCkJcXBwqKyvx4x//2DK9srISP/zhD23mDw0NxalTp6ymlZSU4ODBg9i9ezdiYmJEa9v0x9YCbu73Spz88lsQIumlpqYiNTXV7fmLi4ut3v/mN7/Be++9hz/+8Y+YPHmy25/DgKswtXydL9bwBF8x3BJ5Li8vD/PmzUN8fDwSEhKwadMm1NfXIysrC8Dt4QWXLl3Ctm3bEBAQgNjYWKvlhw4diuDgYJvpatG7X1DqBJ6IvtP7QlOTyQSTyST6erq7u9He3o6BAwd6tBwDrkiUqt6qKRCKUb0lIs+lpaXhypUrWL16NcxmM2JjY1FRUYHo6GgAgNlsdnlPXCnY2+/v7CvdOcG318e1j+5myCXd2N+4XpH19rsQgECT5/tRV8ftZXpfaLpy5UoUFBSI0TQrr776Kq5du4Y5c+Z4tBwDroJ8rd6KFW5ZvSXSvuzsbGRnZ9v9WVlZmdNlCwoKJDkw+UoPfQKH9pBeNTQ0WI3Fl6J6u3PnThQUFOC9997D0KFDPVqWAVcEYo7j0uo4VV+rt2rbHiKSTu/Q5803PKziEilL6otNy8vLkZGRgV27duGRRx7xeHn2DhKS6szdVRj0pNNXS/WWiPTJ2wIAT3qJ/NfOnTsxf/587NixA0888YRXn6F4wC0pKUFMTAyCg4MRFxeH6upqt5Y7evQojEYjJk2aJG0DXRDzwQ5aHafK6q38ePU3aYW31Vki0oerV6/i5MmTOHnyJADg4sWLOHnypOW6gPz8fKSnp1vm37lzJ9LT0/Hqq69i6tSpaGpqQlNTE1pbWz1ar6IBt7y8HLm5uVixYgVqa2uRnJyM1NRUlxdDtLa2Ij09HTNnzpSppZ5TqnorFrmqt1JsT7+YVrsv8h6r+SS2nhNjrZ7YE5F7Tpw4gcmTJ1tu8ZWXl4fJkyfjpZdeAmB7EezGjRtx69YtLFq0CBEREZbXc88959F6FR2DW1RUhIyMDGRmZgK4fe+zAwcOoLS0FIWFhQ6XW7hwIebOnYvAwECXT8SQktzVW7HDoBhVQDXdOcFViO0X04qrF8Nkag0ReYPVWyJ9mT59OgRBcPjz3hfBHjp0SJT1KlbB7ezsRE1NDVJSUqymp6Sk4NixYw6X27p1K86fP4+VK1dK3USvKVntUstFF3IGdk8qtHJXctXy+yBSA38bWsNvPoiUo1gFt6WlBV1dXTaPkgwPD7d55GSPc+fOYfny5aiurobR6F7TOzo60NHRYXnf+8bE3tJz9dbdTtmX25yJsT3ehlVWcomU46x/8bRf4J0UiMgRxXsGg8Fg9V4QBJtpANDV1YW5c+di1apVGDNmjNufX1hYiLCwMMur942JxSbFGbsav7JTemiCr5VYjsklUg7H3xKR1BQLuIMHD0ZgYKBNtba5udmmqgsA7e3tOHHiBBYvXgyj0Qij0YjVq1fjr3/9K4xGIw4ePGh3Pfn5+WhtbbW8GhoafG671u+coPStwXwN7P4eTv3ta17SptQRS5RuglPcj4j0TbEhCkFBQYiLi0NlZSV+/OMfW6ZXVlbihz/8oc38oaGhOHXqlNW0kpISHDx4ELt370ZMTIzd9Yj9bGS5w6EnYVCur+p8qd6qKdxyqIJ7OI6QvOHpECY1flNFRNql6F0U8vLyMG/ePMTHxyMhIQGbNm1CfX09srKyANyuvl66dAnbtm1DQEAAYmNjrZYfOnQogoODbaarkVqqt86oPcj4e+WWSEvU1r/JTe39KanD/sb1SjdBtxQNuGlpabhy5QpWr14Ns9mM2NhYVFRUIDo6GoDtvdGUpubqrbt8/VpOqeotwy2Rvvh7ACYiaSkacAEgOzsb2dnZdn/W+95ovRUUFKCgoED8RtnhSzD0pnrraRj0dXiCrwFdqoMVgy2RPtx9thMhn16y2x9yeAIRiU3xuyhogatwq5WvoqSs3rrizQFM6+FWijHRvDCGiIjINQZcH7kKt3JUb30ldfVWzdUZrYdoqWnl5I1ITbjfECmPAdcFuStm3oRBdyqFSlZvvcHgSaRfd54Uq/kEmIi0iwHXCV+HJvjDnROkqN4y3BKRlDjUh0j/FL/ITI2U6vykqmQoVb1luCWiO8n9TRAR+S9WcHtxNwyqpXrry4VMartzAsOtc6w6kZbZ6284PIGIpOK3FdwBX3TCaPQuHHoTbl3RWvVW7KEJDLfqwwtlSCxSnZzJ9fRGT3C/IVIH9fUOOuUsEHobbqWu3so1NEGP4VaNB14iNVDbdQhEpE88CntILUMT3CFV1UTM7dFjuCXyd9MfW+tyHg5PICIpMeB6QCuP43VFLdVbhlsifVJzlZZj2Yn8AwOum7wNhVJ19Ep9BS7WUAuGW8/IfVDmOELtKSkpQUxMDIKDgxEXF4fq6mqH8x45cgRJSUkYNGgQQkJCMG7cOLz22muSts8f7qDA/YY8sb9xvdJN0DW/vcjME1J1WlJWbx0FIl+qt2quyhD5s/LycuTm5qKkpARJSUnYuHEjUlNTUVdXh5EjR9rM37dvXyxevBgPPPAA+vbtiyNHjmDhwoXo27cv/u3f/k2ydrIPISK5sILrgrvh1tPqrS/hVo0XMLF6+x01/n5I34qKipCRkYHMzEyMHz8excXFiIqKQmlpqd35J0+ejJ/97Ge47777MGrUKPzrv/4rZs2a5bTq64ve/SjH3xKR1HgkdkKrXzfJXb1luCVSTmdnJ2pqapCSkmI1PSUlBceOHXPrM2pra3Hs2DE89NBDDufp6OhAW1ub1Yu+o9XjBZFeMeA64ElnJWf1VkpSj5FjuPUOx9+SMy0tLejq6kJ4eLjV9PDwcDQ1NTlddsSIETCZTIiPj8eiRYuQmZnpcN7CwkKEhYVZXlFRUW61T00XdampLUQkLQZcO3wNt1Jy9fW3L9VbR8QI6wy3RNIyGAxW7wVBsJnWW3V1NU6cOIENGzaguLgYO3fudDhvfn4+WltbLa+Ghga323ZnP6nWE3wi0hdeZPZ/xKxaqa16q/TQBH8Ktxx/S3IbPHgwAgMDbaq1zc3NNlXd3mJiYgAA999/Py5fvoyCggL87Gc/szuvyWSCyWTyup28wIyI5OS3R+OQOjNCPr1keXlD7qEJ3lZviUi/goKCEBcXh8rKSqvplZWVSExMdPtzBEFAR0eH2M2TjJpOJjmsh0h9WMH1ktru6ajWC8vUXr29ejFM6SY4xfG35I68vDzMmzcP8fHxSEhIwKZNm1BfX4+srCwAt4cXXLp0Cdu2bQMAvPXWWxg5ciTGjRsH4PZ9cV955RXk5OQotg1ERGJiwBWRktVbb/lzuCXSi7S0NFy5cgWrV6+G2WxGbGwsKioqEB0dDQAwm82or6+3zN/d3Y38/HxcvHgRRqMRo0ePxssvv4yFCxeK1qb20d3ofz4AIZ9esvQzSo6/5TdcRP6FAdcLnlZvpe7UpbiwzFf+GG7V9JUp+Z/s7GxkZ2fb/VlZWZnV+5ycHEmrtf50IRm/9SBSJx6RPaTEU760Vr31x3ArBQ5PID24ETvcaZ/SL6aVfQYRiY4B1wNSPwTBG2JXbxluicgXE597zfJvVydpd/YX7DuISEwMuCond/WWiEgMrvouOQOtVN+G8FsP8tb+xvVKN0H3GHDdpKfqLYcmiI/jb4mcE7M/5P5GRK6wl3CDUuHWm05cioqCXsOtmm8RxvG3pGda7C+IyHslJSWIiYlBcHAw4uLiUF1d7XT+7du3Y+LEibjrrrsQERGBBQsW4MqVKx6tkwHXBbV+le9tAPI0rOs13BKR9PR84qTnbSMSU3l5OXJzc7FixQrU1tYiOTkZqampVrcuvNORI0eQnp6OjIwMnD59Grt27cLHH3+MzMxMj9bLgOuEs3Crp+qtL3eAYLglIkdc3UFBDrz/LZGyioqKkJGRgczMTIwfPx7FxcWIiopCaWmp3fmPHz+OUaNGYcmSJYiJicGDDz6IhQsX4sSJEx6tlwHXASXDrStiV28d8ad7WfpC7PGAHJ5AenNnX8KTYiJ9aGtrs3rZe9R3Z2cnampqkJKSYjU9JSUFx44ds/u5iYmJaGxsREVFBQRBwOXLl7F792488cQTHrWPD3roxVUIlKMa4SwwyXVhGYcmEJE3/OHEmCeFpAcDvuiE0eh5gebWrds5JCoqymr6ypUrUVBQYDWtpaUFXV1dCA8Pt5oeHh6OpqYmu5+fmJiI7du3Iy0tDd9++y1u3bqFJ598Em+88YZH7fTbgHtjQgT6f3Z7wLJY42y11rH789AENV9gRqQH9k7Gxeg3eAcFInVoaGhAaGio5b3JZHI4r8FgsHovCILNtB51dXVYsmQJXnrpJcyaNQtmsxnLli1DVlYWNm/e7Hb7/DbgAp4HW71Vbx1xFdS1Hm7FxOEJROrF8bdE0gkNDbUKuPYMHjwYgYGBNtXa5uZmm6puj8LCQiQlJWHZsmUAgAceeAB9+/ZFcnIy1qxZg4iICLfax1NhN7kKt0pVb+UemsBwS0TuUvoCMynwpJDIfUFBQYiLi0NlZaXV9MrKSiQmJtpd5vr16wgIsI6ngYGBAG5Xft3FgOsGucKtN9VbbzDcEpFcevoV9h9E/ikvLw9vv/02tmzZgjNnzuD5559HfX09srKyAAD5+flIT0+3zD979mzs3bsXpaWluHDhAo4ePYolS5ZgypQpiIyMdHu9fj1EwR1qrdwC3lVv9VhR8ZRY4285PIHIPo6TJaIeaWlpuHLlClavXg2z2YzY2FhUVFQgOjoaAGA2m63uiTt//ny0t7fjzTffxAsvvIABAwZgxowZWLt2rUfrZcB1QO4g6Gn1VqvjblNGnnX4sw/rx4qyDiJSxsTnXgNGS/f5ngRnKU4YeVJI5J3s7GxkZ2fb/VlZWZnNtJycHOTk5Pi0TgZcO9wNtxya4B5nobb3fP4acnkxDOnB3Wc70T6ahxUiUh6/R7rD12ODZA+33hBzaIJawq238ytF61/BshJFUnDUz3D8LdF39jeuV7oJfsFvT7W/+V4QAk3eDUMQM9yKXb0V656+gG8HJbUGVd7/lkh6d5/txNdjgzR3b3BneFJIpC3aLkMpQOkO25tO1pvqrZLhVq3huAcvLiOyT+n+sQeH/BCR31Zw1UCOC8t8eRSvp9QeTIlIu5QcFsSTQiLtYQXXA3oZmiDFuFsthFs1Dk9gpYn0zt3+RI37JxFpFwOum5T+6k3MoQnO6DXcioUXlxGpG08aiQjgEAW3iB1u1To0wdNw60/Blog8o3RRQCw8KSTSJm2XoyTWPrpb1k5arKEJDLe2xPj6kxeXkZqVlJQgJiYGwcHBiIuLQ3V1tcN59+7di0cffRRDhgxBaGgoEhIScODAAVHa0f98gKUPEuv2YFr/5oSI5MdewwGpgq2nHbVYocSfwy2R3pWXlyM3NxcrVqxAbW0tkpOTkZqaavX4yztVVVXh0UcfRUVFBWpqavDwww9j9uzZqK2tlbnl6saTQiLt4hCFXqSs2Co1NIHhVn04TpDEVFRUhIyMDGRmZgIAiouLceDAAZSWlqKwsNBm/uLiYqv3v/nNb/Dee+/hj3/8IyZPnixHkyXB/YqIerCC+3/kHo5wJ6mHJohFq+FWjcMT5MZKlH51dnaipqYGKSkpVtNTUlJw7Ngxtz6ju7sb7e3tGDhwoMN5Ojo60NbWZvWyhyGTiNRA20dtH1y9p9sSauUItmIOTVBi3K1Ww60aMQCQmFpaWtDV1YXw8HCr6eHh4WhqanLrM1599VVcu3YNc+bMcThPYWEhwsLCLK+oqCiH8+rhAjOeFBJpm98GXDmJPTTBXWKE25SRZzUdblm95YHaXxgMBqv3giDYTLNn586dKCgoQHl5OYYOHepwvvz8fLS2tlpeDQ0NVj+f+NxrVu89OYF2tp+6u//xxJGI7sQxuBLT8rhbpYLth/VjFVmvHHgQJrENHjwYgYGBNtXa5uZmm6pub+Xl5cjIyMCuXbvwyCOPOJ3XZDLBZDI5/HlPn3O7P9J2BZcnhUTap3hpSi23ttECfwi3YmL1lgdqfxAUFIS4uDhUVlZaTa+srERiYqLD5Xbu3In58+djx44deOKJJ6RuJhEB2N+4Xukm+A1Fj956v7WNmNVbMS4q86dwq0as3pJU8vLy8Pbbb2PLli04c+YMnn/+edTX1yMrKwvA7eEF6enplvl37tyJ9PR0vPrqq5g6dSqamprQ1NSE1lZx7lsrNzH3LZ4UEumDokMU9HxrG6nH3Xp6URnDredYvSWtSEtLw5UrV7B69WqYzWbExsaioqIC0dHRAACz2WxVONi4cSNu3bqFRYsWYdGiRZbpzzzzDMrKynxuDx/wQERKUyzg9tzaZvny5VbTpbi1TUdHh+W9o1vbyEWscbf26CHcijX+VozhCWJi9Zaklp2djezsbLs/6x1aDx06JFk72kd3o59kn05E5B7FTo/VeGsbsYhZdfB13K2Wwq2aaL1yxOotKcGb/cbXk1EOTyAiexQ/iit9axuxKTHuluHWGqu3ROrQc5tBvfUxRKR+ig1RUMutbcSkpnDrLjUdeNRyezBWb4mUJ/d+yP2GSF8UO5Lr7dY2SlxU5uvtwNQUbsXC6i2ROsjRv3D/IiJHFL2LQl5eHubNm4f4+HgkJCRg06ZNNre2uXTpErZt2wbgu1vbvP7665Zb2wBASEgIwsKUCzaehltXeldvxb4dmB6DrVhYvSXyjVh3UJAT9xsi/VE04Krt1jZy8XRogj3ejrtVa7gVY3iCr9VbscMtq0vkbzy9g4KjfdadfZH7FxE5o/ijetVyaxtvqWncrVbDrRg4NIFVKFIPLfU13G+I9Enb38cqjOFWHGq4uIxDE4iIiPRD8QquFrkKQwy38mL1lkhZ96x/Vdb1ibWP8cSQSL+0XbZSgDfh1lP+FG71Vr3l0ATyZ55cYObL+FsiIlfYk3jA23DrSfXWn8KtGNR2YRmRv9NKv8MTQ5Lb/sb1SjfBr3CIght8uaKX4dYxX6u3HJrAgzSpR//zARBi7P9MzG9qOASIiNzB8pULDLf6xaEJROLx9amKgHzfqHDfIdI/VnDt8KSTZbj1jtLVWw5NINIeVm+JyF1+e5TvdyEA/c/bf7lLjHBrD8Otc2oLt6zeEt3mbh+k5PAi7jtE8ispKUFMTAyCg4MRFxeH6upqt5Y7evQojEYjJk2a5PE6/Tbg+kqscNu7eutNuE0ZeVYz4VZvGG6JbvP1Eb1yfKvCfYdIfuXl5cjNzcWKFStQW1uL5ORkpKamWj2p1p7W1lakp6dj5syZXq2XAdcLagu3WqKn6i3DLZF8ODyBSJuKioqQkZGBzMxMjB8/HsXFxYiKikJpaanT5RYuXIi5c+ciISHBq/Uy4HqI4dZ7egq3RKQtPDkkEldbW5vVq6Ojw2aezs5O1NTUICUlxWp6SkoKjh075vCzt27divPnz2PlypVet48XmbnJWfWA4VZ6agu3rN4Suaf3ia29fVmOB+gQkbWQOjOMAbYXubtyq/v2/hgVFWU1feXKlSgoKLCa1tLSgq6uLoSHh1tNDw8PR1NTk93PP3fuHJYvX47q6moYjd7HVAZcF1x1rAy37vGlestwy3BL6qfmfon7D5H4GhoaEBoaanlvMpkczmswGKzeC4JgMw0Aurq6MHfuXKxatQpjxozxqX0MuA64E2K0Gm7n3f3d1wK//zrR589zRclwKzaGWyJ5sXpLpE6hoaFWAdeewYMHIzAw0KZa29zcbFPVBYD29nacOHECtbW1WLx4MQCgu7sbgiDAaDTiww8/xIwZM9xqHwPu//GkE3UVONQYbu8MtfamSxV0xXyCkTd4URmRukk9Np77EJFygoKCEBcXh8rKSvz4xz+2TK+srMQPf/hDm/lDQ0Nx6tQpq2klJSU4ePAgdu/ejZgYB49LtMNvA+6ALzphNHresSoRbqUItlqgpqEJrCIROebuLcLk/kaG4ZZIeXl5eZg3bx7i4+ORkJCATZs2ob6+HllZWQCA/Px8XLp0Cdu2bUNAQABiY2Otlh86dCiCg4NtprvCy8rdFPLpJU2F23l3H1M83Opl3K1S4ZYHZ/KEJzdSN5vNmDt3LsaOHYuAgADk5uZ6vd4HNrzh9bKu8MSSSPvS0tJQXFyM1atXY9KkSaiqqkJFRQWio6MB3O6PXN0T1xt+W8F1l7shw9vH74o9JEHpUNuD4dY3DLfkiZ4bqZeUlCApKQkbN25Eamoq6urqMHLkSJv5Ozo6MGTIEKxYsQKvvfaaZO1y1Q9IOTyB+xCRemRnZyM7O9vuz8rKypwuW1BQYHN3BnewguuAOxXbHmoIt75UbMUef8tw6xsemMlTnt5IfdSoUXj99deRnp6OsDDfhw240195um/7sv9xHyK12d+4Xukm+B2/reD23P/tRuxwrztDd4YkAPKEW7VguPUND8zkqZ4bqS9fvtxquqsbqRMR6ZnfBtweagq3SgRbMau3DLe+Ybglb3hzI3VvdHR0WD2pqK2tzafPc7bPsnpLRL7iEAUv6CXcionh1jc8KJOv3L2RurcKCwsRFhZmefV+ilFvStwikPsREfXw+wquJ3wJtoA44VbMYCtW9dbbA5magi3AcEva5OmN1L2Vn5+PvLw8y/u2tjaXIbdH731dquotEVEPVnDdcCN2uM9VW4Zbawy3tzHckq/uvJH6nSorK5GYKN4QJJPJZHly0Z1PMOob7dtQBbFwXyKiO7GC64S9UNtDziEJYg9HYLi1xieUkdZ5ciP1HidPngQAXL16Ff/4xz9w8uRJBAUFYcKECV61Yd7dx3zuW7zdF7k/EVFvDLi9OAu1gP1gC/hXuOV4W9/wYExiS0tLw5UrV7B69WqYzWbExsa6vJH65MmTLf+uqanBjh07EB0djb///e8erXvGiHOYd/fHDn/uyfAEb3B/IiJ7/Dbg3pgQAaMx2O35PQm2gO/hVoqLyJQMt2oKtgDDLemPpzdSFwRBsrZ4009w7C0RiclvA667HAVbQLr72zLcWmO4JdIHR/syhyYQkdgYcB3wNNgC6qzaAr6HW6WGJDDYEmmbrye3znC/IiJn/DbgfvO9IASaHIdYezwNtoD/hltWbW/jQZj07ukBH8GdG/KIWb3lfkVErvhtwPWEo2AL+B5upXxggxbDLYMtkTb19DdSP+CB+xYRuYMB1wkpgy3Aqm1verhDAsADMBHgXn/AC8uISCoMuHdwFmh7OAu2gLYvJNNDsAVYtSVSG7H2ce5jROQuvw24V+/pRkCw60DbQ83BFtBmuGWwJfJfnu6v3M9Iq/Y3rle6CX7JbwOuu8QKtgCrtj30EmwBHnSJ7tTTn9zZP9jb3xluiUhqDLh3cBVm78Rg6zkGWyLyFPc1IvKG3wbcvtFtCLyrw+PllA62gG/hlsHWNzzYEtm6s08Ss3rL/Y2IvOW3AdcTnoRaQD/BlheP3caDLJFr7vQxDLdEJBcGXDs8DbQ99DIUgdXa23iAJfKeL30B9z0i8pXfBtwZI87B1K+PaJ+ntnDLYOs9HlyJvOOsD3F3n+b+R0Ri8NuAKxaxgy1DrTJ4UCUST+++geGWiOTGgOslMYMtQ60yeDAlEoezfojhloiUwIDrIX8Mtgy1ROTIH775PoDv+hdv+gvul0QkNgZcD/3+60SvQ66WLhgTM9Qq/bx5HjyJlOFq3+e+SURSYcD1Qu+g2hN4fQ2w9ngaatUQaAFWaYn80Z39CMMtESmJAVcEYgZbVmk9xwMlkXIONt6LG5et+yCGWyJSGgOuwrQYaBlmiai3nv7FWf/A/ZeI5MKAKyOGWc/wYEikH9yfiUhODLgi8ybE3smbQKv1MMsDH5F2XfsyFGGXHFdvuX8TkRL8NuAebLwXgXeZFFm3v1ZleaAj0i+GWyJb+xvXK90Ev+W3AVdsvjw84U5ihFglAiwPZET+qd+FANz9d+s+h/0BESnNbwPutS9DERAcLOs6tVqB5cGKiNzF/oKIeispKcFvf/tbmM1m3HfffSguLkZycrLD+Q8fPoy8vDycPn0akZGRePHFF5GVleXROv024Pa7EIBAk7j3fvWFnKGVByAiEsuALzoBYzD7FSKyq7y8HLm5uSgpKUFSUhI2btyI1NRU1NXVYeTIkTbzX7x4EY8//jieffZZvPPOOzh69Ciys7MxZMgQPPXUU26v128D7oAvOmE0qifguoMHECJSm5A6M4wBQUo3g4hUqqioCBkZGcjMzAQAFBcX48CBAygtLUVhYaHN/Bs2bMDIkSNRXFwMABg/fjxOnDiBV155hQHXGUEQAAB9Pv1Sc53yLaUbQERuudV9+xuZnv5Gj3q2rWdbichWW1ub5J/tbT9zS+gEur1cDrbbZjKZYDJZX7zf2dmJmpoaLF++3Gp6SkoKjh07Zvfz//KXvyAlJcVq2qxZs7B582bcvHkTffr0caudfhdw29vbAQCHLm9VuCVEpHft7e0ICxPnAlS1YV9K5FpY2EbJ1+FpPxMUFIRhw4bhUJP3+26/fv0QFRVlNW3lypUoKCiwmtbS0oKuri6Eh4dbTQ8PD0dTU5Pdz25qarI7/61bt9DS0oKIiAi32uh3ATcyMhINDQ3o378/DAaD0s2RRFtbG6KiotDQ0IDQ0FClmyM6bp/26X0bBUFAe3s7IiMjlW6KZPTUl+r975Hbp22Ots/bfiY4OBgXL15EZ6f3374IgmCz3/eu3t6p97z2lnc1v73pzvhdwA0ICMCIESOUboYsQkNDdbmz9+D2aZ+et1GvldseeuxL9fz3CHD7tM7e9nnbzwQHByNYhjtJDR48GIGBgTbV2ubmZpsqbY9hw4bZnd9oNGLQoEFur1tbV1kRERERkSYEBQUhLi4OlZWVVtMrKyuRmJhod5mEhASb+T/88EPEx8e7Pf4WYMAlIiIiIonk5eXh7bffxpYtW3DmzBk8//zzqK+vt9zXNj8/H+np6Zb5s7Ky8OWXXyIvLw9nzpzBli1bsHnzZixdutSj9frdEAV/YDKZsHLlSqfjYbSM26d9/rCNpB16/3vk9mmb1rcvLS0NV65cwerVq2E2mxEbG4uKigpER0cDAMxmM+rr6y3zx8TEoKKiAs8//zzeeustREZGYv369R7dIgwADIKe72NDRERERH6HQxSIiIiISFcYcImIiIhIVxhwiYiIiEhXGHCJiIiISFcYcDWqpKQEMTExCA4ORlxcHKqrqx3Ou3fvXjz66KMYMmQIQkNDkZCQgAMHDsjYWs95sn13Onr0KIxGIyZNmiRtA33k6fZ1dHRgxYoViI6OhslkwujRo7FlyxaZWusdT7dx+/btmDhxIu666y5ERERgwYIFuHLlikytJb0qKCiAwWCweg0bNszpMocPH0ZcXByCg4Nxzz33YMOGDTK11nOjRo2y2T6DwYBFixbZnf/QoUN25//ss89kbrl9VVVVmD17NiIjI2EwGPDuu+9a/VwQBBQUFCAyMhIhISGYPn06Tp8+7fJz9+zZgwkTJsBkMmHChAnYt2+fRFvgnLPtu3nzJn75y1/i/vvvR9++fREZGYn09HR89dVXTj+zrKzM7u/022+/lXhr1I0BV4PKy8uRm5uLFStWoLa2FsnJyUhNTbW6zcadqqqq8Oijj6KiogI1NTV4+OGHMXv2bNTW1srccvd4un09WltbkZ6ejpkzZ8rUUu94s31z5szBf//3f2Pz5s04e/Ysdu7ciXHjxsnYas94uo1HjhxBeno6MjIycPr0aezatQsff/wxMjMzZW456dF9990Hs9lseZ06dcrhvBcvXsTjjz+O5ORk1NbW4t///d+xZMkS7NmzR8YWu+/jjz+22raeG+T/9Kc/dbrc2bNnrZa799575WiuS9euXcPEiRPx5ptv2v35unXrUFRUhDfffBMff/wxhg0bhkcffRTt7e0OP/Mvf/kL0tLSMG/ePPz1r3/FvHnzMGfOHHz00UdSbYZDzrbv+vXr+OSTT/D//t//wyeffIK9e/fi888/x5NPPunyc0NDQ61+n2azWZYnlamaQJozZcoUISsry2rauHHjhOXLl7v9GRMmTBBWrVoldtNE4e32paWlCb/61a+ElStXChMnTpSwhb7xdPv2798vhIWFCVeuXJGjeaLwdBt/+9vfCvfcc4/VtPXr1wsjRoyQrI3kHzztD1588UVh3LhxVtMWLlwoTJ06VeSWSeO5554TRo8eLXR3d9v9+Z///GcBgPD111/L2zAvABD27dtned/d3S0MGzZMePnlly3Tvv32WyEsLEzYsGGDw8+ZM2eO8Nhjj1lNmzVrlvD000+L3mZP9N4+e/7nf/5HACB8+eWXDufZunWrEBYWJm7jdIAVXI3p7OxETU0NUlJSrKanpKTg2LFjbn1Gd3c32tvbMXDgQCma6BNvt2/r1q04f/48Vq5cKXUTfeLN9r3//vuIj4/HunXrMHz4cIwZMwZLly7FjRs35Giyx7zZxsTERDQ2NqKiogKCIODy5cvYvXs3nnjiCTmaTDp37tw5REZGIiYmBk8//TQuXLjgcN6//OUvNn+7s2bNwokTJ3Dz5k2pm+qTzs5OvPPOO/jFL34Bg8HgdN7JkycjIiICM2fOxJ///GeZWuibixcvoqmpyer3YzKZ8NBDDzk9Pjj6nbp7zFRSa2srDAYDBgwY4HS+q1evIjo6GiNGjMAPfvAD1X5DKycGXI1paWlBV1cXwsPDraaHh4ejqanJrc949dVXce3aNcyZM0eKJvrEm+07d+4cli9fju3bt8NoVPfD+bzZvgsXLuDIkSP49NNPsW/fPhQXF2P37t0Ox9gpzZttTExMxPbt25GWloagoCAMGzYMAwYMwBtvvCFHk0nHvv/972Pbtm04cOAA/vM//xNNTU1ITEx0OL67qanJ7t/urVu30NLSIkeTvfbuu+/im2++wfz58x3OExERgU2bNmHPnj3Yu3cvxo4di5kzZ6Kqqkq+hnqpp//w9Pjn6Hfq7jFTKd9++y2WL1+OuXPnIjQ01OF848aNQ1lZGd5//33s3LkTwcHBSEpKwrlz52RsrfqoOw2QQ73PzgVBcHnGDgA7d+5EQUEB3nvvPQwdOlSq5vnM3e3r6urC3LlzsWrVKowZM0au5vnMk99fd3c3DAYDtm/fjrCwMABAUVERfvKTn+Ctt95CSEiI5O31hifbWFdXhyVLluCll17CrFmzYDabsWzZMmRlZWHz5s1yNJd0KjU11fLv+++/HwkJCRg9ejR+97vfIS8vz+4y9v527U1Xm82bNyM1NRWRkZEO5xk7dizGjh1reZ+QkICGhga88sormDZtmhzN9Jk3xz9vj5lKuXnzJp5++ml0d3ejpKTE6bxTp07F1KlTLe+TkpLwT//0T3jjjTewfv16qZuqWgy4GjN48GAEBgbanHk2NzfbnKH2Vl5ejoyMDOzatQuPPPKIlM30mqfb197ejhMnTqC2thaLFy8GcDsQCoIAo9GIDz/8EDNmzJCl7e7w5vcXERGB4cOHW8ItAIwfPx6CIKCxsVE1F4f08GYbCwsLkZSUhGXLlgEAHnjgAfTt2xfJyclYs2YNIiIiJG83+Ye+ffvi/vvvd1jdGjZsmN2/XaPRiEGDBsnRRK98+eWX+K//+i/s3bvX42WnTp2Kd955R4JWiavn7hdNTU1WfYKr45+j36mrY6ZSbt68iTlz5uDixYs4ePCg0+qtPQEBAfjnf/5nv6/gcoiCxgQFBSEuLs5ypWyPyspKJCYmOlxu586dmD9/Pnbs2KHqcY2ebl9oaChOnTqFkydPWl5ZWVkYO3YsTp48ie9///tyNd0t3vz+kpKS8NVXX+Hq1auWaZ9//jkCAgIwYsQISdvrDW+28fr16wgIsO6OAgMDAXxXPSMSQ0dHB86cOePwpCkhIcHmb/fDDz9EfHw8+vTpI0cTvbJ161YMHTrUq/69trZWEyeRMTExGDZsmNXvp7OzE4cPH3Z6/HP0O3W2jFJ6wu25c+fwX//1X16dVAmCgJMnT2ridyopZa5tI1/84Q9/EPr06SNs3rxZqKurE3Jzc4W+ffsKf//73wVBEITly5cL8+bNs8y/Y8cOwWg0Cm+99ZZgNpstr2+++UapTXDK0+3rTe13UfB0+9rb24URI0YIP/nJT4TTp08Lhw8fFu69914hMzNTqU1wydNt3Lp1q2A0GoWSkhLh/PnzwpEjR4T4+HhhypQpSm0C6cQLL7wgHDp0SLhw4YJw/Phx4Qc/+IHQv39/h3+LFy5cEO666y7h+eefF+rq6oTNmzcLffr0EXbv3q3UJrjU1dUljBw5UvjlL39p87Pe2/faa68J+/btEz7//HPh008/FZYvXy4AEPbs2SNnkx1qb28XamtrhdraWgGAUFRUJNTW1lruIvDyyy8LYWFhwt69e4VTp04JP/vZz4SIiAihra3N8hnz5s2zumPL0aNHhcDAQOHll18Wzpw5I7z88suC0WgUjh8/rqrtu3nzpvDkk08KI0aMEE6ePGl1vO7o6HC4fQUFBcIHH3wgnD9/XqitrRUWLFggGI1G4aOPPpJ9+9SEAVej3nrrLSE6OloICgoS/umf/kk4fPiw5WfPPPOM8NBDD1neP/TQQwIAm9czzzwjf8Pd5Mn29ab2gCsInm/fmTNnhEceeUQICQkRRowYIeTl5QnXr1+XudWe8XQb169fL0yYMEEICQkRIiIihJ///OdCY2OjzK0mvUlLSxMiIiKEPn36CJGRkcK//Mu/CKdPn7b83N7f4qFDh4TJkycLQUFBwqhRo4TS0lKZW+2ZAwcOCACEs2fP2vys9/atXbtWGD16tBAcHCzcfffdwoMPPij86U9/krG1zvXcxszR8aq7u1tYuXKlMGzYMMFkMgnTpk0TTp06ZfUZDz30kM3xbdeuXcLYsWOFPn36COPGjVMs0DvbvosXL9r9GQDhz3/+s+Uzem9fbm6uMHLkSCEoKEgYMmSIkJKSIhw7dkz+jVMZgyDw+z8iIiIi0g+OwSUiIiIiXWHAJSIiIiJdYcAlIiIiIl1hwCUiIiIiXWHAJSIiIiJdYcAlIiIiIl1hwCUiIiIiXWHAJSIiIiJdYcAlIiIiIl1hwCUiIiIiXWHAJXJg1KhRKC4utpo2adIkFBQUKNIeIiJfTZ8+HYsXL8bixYsxYMAADBo0CL/61a8gCILSTSMSFQMuERGRH/nd734Ho9GIjz76COvXr8drr72Gt99+W+lmEYnKqHQDiIiISD5RUVF47bXXYDAYMHbsWJw6dQqvvfYann32WaWbRiQaVnCJiIj8yNSpU2EwGCzvExIScO7cOXR1dSnYKiJxMeASORAQEGAzLu3mzZsKtYaIiIjcxYBL5MCQIUNgNpst79va2nDx4kUFW0RE5Lvjx4/bvL/33nsRGBioUIuIxMeAS+TAjBkz8Pvf/x7V1dX49NNP8cwzz/AAQESa19DQgLy8PJw9exY7d+7EG2+8geeee07pZhGJiheZETmQn5+PCxcu4Ac/+AHCwsLw61//mhVcItK89PR03LhxA1OmTEFgYCBycnLwb//2b0o3i0hUBoE3vyMiIvIL06dPx6RJk2zu8U2kNxyiQERERES6woBLRERERLrCIQpEREREpCus4BIRERGRrjDgEhEREZGuMOASERERka4w4BIRERGRrjDgEhEREZGuMOASERERka4w4BIRERGRrjDgEhEREZGuMOASERERka78f0EkiZVcgb+pAAAAAElFTkSuQmCC",
                        "text/plain": [
                            "<Figure size 800x400 with 3 Axes>"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                }
            ],
            "source": [
                "# Look at the contours of the (7,-3, 1) mode using teuk_specific_amps\n",
                "plt.figure(figsize=(8, 4))\n",
                "plt.subplot(121)\n",
                "cb = plt.contourf(\n",
                "    u_grid.reshape(uv.size, wv.size),\n",
                "    w_grid.reshape(uv.size, wv.size),\n",
                "    np.abs(teuk_specific_amps[(7, -3, 0, 1)].reshape(uv.size, wv.size)),\n",
                ")\n",
                "plt.xlabel(\"u\")\n",
                "plt.ylabel(\"w\")\n",
                "plt.subplot(122)\n",
                "cb = plt.contourf(\n",
                "    p_grid.reshape(uv.size, wv.size),\n",
                "    e_grid.reshape(uv.size, wv.size),\n",
                "    np.abs(teuk_specific_amps[(7, -3, 0, 1)].reshape(uv.size, wv.size)),\n",
                ")\n",
                "plt.colorbar(cb)\n",
                "plt.xlabel(\"p\")\n",
                "plt.ylabel(\"e\")\n",
                "plt.show()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "### ROMAN amplitude generation"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "ROMAN uses a reduced order model representing the mode amplitude space. It then trains a neural network on this reduced model. The neural network is evaluated and the resulting matrix is transformed from the reduced basis back to the full mode space.  \n",
                "\n",
                "This approach for amplitude interpolation was the primary method used for the original Schwarzschild eccentric waveforms. It is retained for compatibility (both backwards and, perhaps, future). The interface is very similar to that of the spline interpolation classes."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 9,
            "metadata": {},
            "outputs": [],
            "source": [
                "from few.amplitude.romannet import RomanAmplitude\n",
                "\n",
                "# initialize ROMAN class\n",
                "amp = RomanAmplitude(buffer_length=5000)  # buffer_length creates memory buffers"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 10,
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "Total trajectory points: 2500\n"
                    ]
                },
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAG2CAYAAABViX0rAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAO7ZJREFUeJzt3Xt4VNW9//FPLiRDoYlyCwGTEC1CIFohaSFEvBQMxR5vPS1B2qCnUImgElNbQ4EmcJQoWoi3RGJRDlUxT4u3tumB8LNIMFSPMemx0iIqNTFOhGBNAE3GJPv3B4cpwyRhZpKZ2TPzfj3PPI+zZu09a7mU/eG79+wdZhiGIQAAABMJ9/cAAAAAzkRAAQAApkNAAQAApkNAAQAApkNAAQAApkNAAQAApkNAAQAApkNAAQAApkNAAQAApkNAAQAApuP3gFJaWqrk5GRZLBalpaWpurq6174333yzwsLCnF6TJ0/24YgBAAge7hyHd+/e3eNx+O9//7tDv88++0zLli1TfHy8LBaLUlJSVFlZ6da4/BpQKioqlJeXp5UrV6qurk4zZ87U3Llz1dDQ0GP/hx56SFar1f5qbGzUsGHD9P3vf9/HIwcAIPC5exw+5cCBAw7H4/Hjx9s/s9lsuuqqq/SPf/xDv/3tb3XgwAE98cQTGjt2rFtjC/PnwwKnTZumqVOnqqyszN6WkpKi66+/XsXFxWfd/sUXX9R3v/tdHTp0SElJSd4cKgAAQcfd4/Du3bt15ZVX6p///KfOOeecHvf5+OOP64EHHtDf//53DRo0yOOxRXq8ZT/ZbDbV1taqoKDAoT0rK0s1NTUu7WPz5s2aPXt2n+Gko6NDHR0d9vfd3d369NNPNXz4cIWFhXk2eABASDAMQ8eOHdOYMWMUHu69kw7t7e2y2Wz93o9hGE7HtujoaEVHRzv17c9xeMqUKWpvb9ekSZO0atUqXXnllfbPXn75ZWVkZGjZsmV66aWXNHLkSC1YsEB33323IiIiXJ6L3wJKS0uLurq6FBcX59AeFxen5ubms25vtVr1xz/+Uc8++2yf/YqLi7VmzZp+jRUAENoaGxt13nnneWXf7e3tSkgcopYj3f3e19ChQ3X8+HGHtsLCQhUVFTn19eQ4HB8fr/LycqWlpamjo0O//vWvNWvWLO3evVuXXXaZJOmDDz7QK6+8oh/84AeqrKzUwYMHtWzZMnV2duoXv/iFy3PxW0A55cyk11P668mWLVt0zjnn6Prrr++z34oVK5Sfn29/39raqsTERE27vECRkRaPxgwAnhq83+rvIcANnd027f7kKX31q1/12nfYbDa1HOnWjj+P1pChnldpThzv1pzpzWpsbFRMTIy9vafqyencOQ5PmDBBEyZMsL/PyMhQY2OjHnzwQXtA6e7u1qhRo1ReXq6IiAilpaXp448/1gMPPBAYAWXEiBGKiIhwSmmHDx92SnNnMgxDTz75pHJychQVFdVn395KW5GRFgIKAJ8Z/Nemk/8Q3vefWTAnX1wSMGRouIZ+tf+nkWJiYhwCSm/6cxw+3fTp0/X000/b38fHx2vQoEEOp3NSUlLU3Nwsm8121uP2KX77FU9UVJTS0tJUVVXl0F5VVaUZM2b0ue2rr76q9957T4sWLfLmEAFgQNjDCWAi/TkOn66urk7x8fH295mZmXrvvffU3f2vU1bvvvuu4uPjXQ4nkp9P8eTn5ysnJ0fp6enKyMhQeXm5GhoalJubK+nk6ZmmpiZt3brVYbvNmzdr2rRpSk1N9cewAcBlhBOYmbvH4ZKSEo0bN06TJ0+WzWbT008/re3bt2v79u32fd5666165JFHtHz5ct1+++06ePCg1q1bpzvuuMOtsfk1oGRnZ+vo0aNau3atrFarUlNTVVlZaf9VjtVqdfotdmtrq7Zv366HHnrIH0MGAJcQTBAI3D0O22w23XXXXWpqatLgwYM1efJk/eEPf9DVV19t75OQkKCdO3fqzjvv1MUXX6yxY8dq+fLluvvuu90am1/vg+IPbW1tio2NVeasIq5BAeAVhJPg0dlt0y7rJrW2trp0XYcnTh2X9v51TL+uQTl+rFuXpn7s1bH6kt9vdQ8AwYRwAgwMAgoADBDCCTBw/H4fFAAIdAQTYOARUADAQwQTwHs4xQMAHiCcAN5FBQUA3EAwAXyDCgoAuIhwAvgOFRQAOAuCCeB7VFAAoA+EE8A/qKAAQA8IJoB/EVAA4DQEE8AcCCgAQh6hBDAfAgqAkEUwAcyLgAIg5BBMAPMjoAAIGQQTIHAQUAAEPYIJEHgIKACCFsEECFwEFABBh2ACBD4CCoCgQCgBggsBBUBAI5gAwYmAAiAgEUyA4EZAARAwCCVA6CCgADA9ggkQeggoAEyLYAKELgIKAFMhlACQCCgATIJgAuB0BBQAfkMoAdAbAgoAnyOYADgbAgoAnyGYAHAVAQWA1xFMALiLgALAKwglAPqDgAJgQBFMAAwEAgqAfiOUABhoBBQAHiGUAPAmAgoAlxFKAPgKAQVAnwglAPyBgAKgRwQTAP5EQAFgRygBYBYEFCDEEUoAmBEBBQhBhBIAZkdAAUIEoQRAICGgAEGMUAIgUIX7ewAABs7gvzY5vADgbEpLS5WcnCyLxaK0tDRVV1e7tN1rr72myMhIXXLJJU6flZSUaMKECRo8eLASEhJ05513qr293a1xUUEBAhxBBICnKioqlJeXp9LSUmVmZmrTpk2aO3eu9u/fr8TExF63a21t1cKFCzVr1ix98sknDp8988wzKigo0JNPPqkZM2bo3Xff1c033yxJ2rhxo8tjo4ICBCCqJAAGwoYNG7Ro0SItXrxYKSkpKikpUUJCgsrKyvrcbsmSJVqwYIEyMjKcPtu3b58yMzO1YMECjRs3TllZWbrxxhv15ptvujU2AgoQIAglAFzR1tbm8Oro6Oixn81mU21trbKyshzas7KyVFNT0+v+n3rqKb3//vsqLCzs8fNLL71UtbW1euONNyRJH3zwgSorK/Wd73zHrXlwigcwMcIIEDqe+2yaojsHebx9x/EvJb2ghIQEh/bCwkIVFRU59W9paVFXV5fi4uIc2uPi4tTc3Nzjdxw8eFAFBQWqrq5WZGTPEWL+/Pk6cuSILr30UhmGoc7OTt16660qKChwaz4EFMBkCCUA+qOxsVExMTH299HR0X32DwsLc3hvGIZTmyR1dXVpwYIFWrNmjS688MJe97d7927de++9Ki0t1bRp0/Tee+9p+fLlio+P1+rVq12eBwEF8DMCCYCBFBMT4xBQejNixAhFREQ4VUsOHz7sVFWRpGPHjunNN99UXV2dbrvtNklSd3e3DMNQZGSkdu7cqW9961tavXq1cnJytHjxYknSRRddpBMnTuiWW27RypUrFR7u2tUlBBTADwglAPwtKipKaWlpqqqq0g033GBvr6qq0nXXXefUPyYmRm+//bZDW2lpqV555RX99re/VXJysiTp888/dwohERERMgxDhmG4PD4CCuADBBIAZpSfn6+cnBylp6crIyND5eXlamhoUG5uriRpxYoVampq0tatWxUeHq7U1FSH7UeNGiWLxeLQfs0112jDhg2aMmWK/RTP6tWrde211yoiIsLlsfn9Vzzu3iCmo6NDK1euVFJSkqKjo3XBBRfoySef9NFoAdfxqxsAZpedna2SkhKtXbtWl1xyifbs2aPKykolJSVJkqxWqxoaGtza56pVq/STn/xEq1at0qRJk7Ro0SLNmTNHmzZtcms/YYY79ZYBVlFRoZycHIcbxPzqV7/q8wYx1113nT755BPdc889+trXvqbDhw+rs7NTM2bMcOk729raFBsbq8xZRYqMtAzkdADCCBBkOrtt2mXdpNbWVpeu6/DEqePSbXtvUPTQ/v2K59FLX/DqWH3Jr6d4Tr9BjHTy1rg7duxQWVmZiouLnfr/93//t1599VV98MEHGjZsmCRp3Lhxvhwy4IRQAgADz2+neDy5QczLL7+s9PR0rV+/XmPHjtWFF16ou+66S1988UWv39PR0eF00xqgvzh9AwDe5bcKiic3iPnggw+0d+9eWSwWvfDCC2ppadHSpUv16aef9nodSnFxsdasWTPg40doIYgAgG/5/SJZV28QI538vXVYWJieeeYZffOb39TVV1+tDRs2aMuWLb1WUVasWKHW1lb7q7GxccDngOBElQQA/MdvFRR3bxAjSfHx8Ro7dqxiY2PtbSkpKTIMQx999JHGjx/vtE10dPRZ76IHSFRJAMBM/FZBOf0GMaerqqrq9Rc5mZmZ+vjjj3X8+HF727vvvqvw8HCdd955Xh0vghNVEgAwJ7+e4snPz9evfvUrPfnkk/rb3/6mO++80+kGMQsXLrT3X7BggYYPH67/+I//0P79+7Vnzx799Kc/1Y9+9CMNHjzYX9NAgCGUAID5+fVnxtnZ2Tp69KjWrl0rq9Wq1NTUPm8QM3ToUFVVVen2229Xenq6hg8frnnz5umee+7x1xQQAAgiABB4/H6r+6VLl2rp0qU9frZlyxantokTJzqdFgLORCgBgMDm94ACDBRCCQAEDwIKAhqhBACCEwEFAYdQAgDBj4AC0yOQAEDoIaDAlAglABDaCCgwDUIJAOAUAgr8ilACAOgJAQU+RygBAJwNAQVeRyABALiLgAKvIJQAAPqDgIIBQygBAAwUAgr6hVACAPAGAgrcRigBAHgbAQUuIZQAAHyJgIJeEUoAAP5CQIEdgQQAYBYElBBHKAEAmBEBJQQRSgAAZkdACQEEEgBAoCGgBClCCQAgkBFQggihBAAQLAgoAY5QAgAIRgSUAEMgAQCEAgJKACCUAABCDQHFpAglAIBQRkAxCQIJAAD/QkDxI0IJAAA9I6D4GKEEAICzI6D4AKEEAAD3EFC8hFACAIDnCCgDiFACAMDAIKD0E6EEAICBR0BxE4EEAADvC/f3AALB4L822V8AAAST0tJSJScny2KxKC0tTdXV1S5t99prrykyMlKXXHKJ02fbt2/XpEmTFB0drUmTJumFF15we1wElF4QSgAAwa6iokJ5eXlauXKl6urqNHPmTM2dO1cNDQ19btfa2qqFCxdq1qxZTp/t27dP2dnZysnJ0V/+8hfl5ORo3rx5ev31190aW5hhGIZbWwS4trY2xcbGKnNWkSIjLQ6fEUYAAKfr7LZpl3WTWltbFRMT45XvOHVcum3vDYoeOsjj/XQc/1KPXvqCW2OdNm2apk6dqrKyMntbSkqKrr/+ehUXF/e63fz58zV+/HhFREToxRdfVH19vf2z7OxstbW16Y9//KO97dvf/rbOPfdcbdu2zeX5hHwFhUoJACCYtLW1Obw6Ojp67Gez2VRbW6usrCyH9qysLNXU1PS6/6eeekrvv/++CgsLe/x83759TvucM2dOn/vsScheJDt4v1WR4VH+HgYAAJKkVz4ar4ivRHu8fdfnJ4NIQkKCQ3thYaGKioqc+re0tKirq0txcXEO7XFxcWpubu7xOw4ePKiCggJVV1crMrLnCNHc3OzWPnsTsgEFAIBg1NjY6HCKJzq679ATFhbm8N4wDKc2Serq6tKCBQu0Zs0aXXjhhQOyz74QUAAACCIxMTEuXYMyYsQIRUREOFU2Dh8+7FQBkaRjx47pzTffVF1dnW677TZJUnd3twzDUGRkpHbu3KlvfetbGj16tMv77EvIX4MCAEAoioqKUlpamqqqqhzaq6qqNGPGDKf+MTExevvtt1VfX29/5ebmasKECaqvr9e0adMkSRkZGU773LlzZ4/77AsVFAAAQlR+fr5ycnKUnp6ujIwMlZeXq6GhQbm5uZKkFStWqKmpSVu3blV4eLhSU1Mdth81apQsFotD+/Lly3XZZZfp/vvv13XXXaeXXnpJu3bt0t69e90aGwEFAIAQlZ2draNHj2rt2rWyWq1KTU1VZWWlkpKSJElWq/Ws90Q504wZM/Tcc89p1apVWr16tS644AJVVFTYKyyuCtn7oMyOX8KveAAAffLlfVAmPfezfv+KZ//89V4dqy9xDQoAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdvweU0tJSJScny2KxKC0tTdXV1b323b17t8LCwpxef//73304YgAA4G1+DSgVFRXKy8vTypUrVVdXp5kzZ2ru3LlnvSnMgQMHZLVa7a/x48f7aMQAAMAX/BpQNmzYoEWLFmnx4sVKSUlRSUmJEhISVFZW1ud2o0aN0ujRo+2viIgIH40YAAD4gt8Cis1mU21trbKyshzas7KyVFNT0+e2U6ZMUXx8vGbNmqU//elPffbt6OhQW1ubwwsAAJib3wJKS0uLurq6nB6/HBcX5/SY5lPi4+NVXl6u7du36/nnn9eECRM0a9Ys7dmzp9fvKS4uVmxsrP2VkJAwoPMAAAADz+8PCwwLC3N4bxiGU9spEyZM0IQJE+zvMzIy1NjYqAcffFCXXXZZj9usWLFC+fn59vdtbW2EFAAATM5vFZQRI0YoIiLCqVpy+PBhp6pKX6ZPn66DBw/2+nl0dLRiYmIcXgAAwNz8FlCioqKUlpamqqoqh/aqqirNmDHD5f3U1dUpPj5+oIcHAAD8yK+nePLz85WTk6P09HRlZGSovLxcDQ0Nys3NlXTy9ExTU5O2bt0qSSopKdG4ceM0efJk2Ww2Pf3009q+fbu2b9/uz2kAAIAB5teAkp2draNHj2rt2rWyWq1KTU1VZWWlkpKSJElWq9Xhnig2m0133XWXmpqaNHjwYE2ePFl/+MMfdPXVV/trCgAAwAvCDMMw/D0IX2pra1NsbKxmxy9RZHiUv4cDADCxzm6bdlk3qbW11WvXMJ46Lk167meK+Eq0x/vp+rxD++ev9+pYfcnvt7oHAAA4EwEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYTqS/BwAAgCe+SB3r9e/o7GyXrF7/GvSAgAIA8ApfBAgELwIKAIQ4ggTMiIACAAGGQIFQQEABAC8jUADuI6AAQC8IFoD/EFAABCXCBRDYCCgATImAAYQ2AgoAryFkAPAUAQXAWRE0APgat7oHQtAXqWPdegEIXqWlpUpOTpbFYlFaWpqqq6t77bt3715lZmZq+PDhGjx4sCZOnKiNGzc69HniiSc0c+ZMnXvuuTr33HM1e/ZsvfHGG26PiwoKEAQIEQA8UVFRoby8PJWWliozM1ObNm3S3LlztX//fiUmJjr1HzJkiG677TZdfPHFGjJkiPbu3aslS5ZoyJAhuuWWWyRJu3fv1o033qgZM2bIYrFo/fr1ysrK0jvvvKOxY13/syrMMAxjwGYaANra2hQbG6vZ8UsUGR7l7+EAvSJ0AP7X2dmu1/5fkVpbWxUTE+OV7zh1XJr03M8U8ZVoj/fT9XmH9s9f79ZYp02bpqlTp6qsrMzelpKSouuvv17FxcUu7eO73/2uhgwZol//+tc9j6urS+eee64effRRLVy40KV9SlRQAJ8idAAwC5vNptraWhUUFDi0Z2VlqaamxqV91NXVqaamRvfcc0+vfT7//HN9+eWXGjZsmFvjI6AA/UToAGAmbW1tDu+jo6MVHe1cmWlpaVFXV5fi4uIc2uPi4tTc3Nznd5x33nk6cuSIOjs7VVRUpMWLF/fat6CgQGPHjtXs2bPdmAUBBegRoQOAr534MEbhFovH23e3t0uSEhISHNoLCwtVVFTU63ZhYWEO7w3DcGo7U3V1tY4fP64///nPKigo0Ne+9jXdeOONTv3Wr1+vbdu2affu3bK4OTcCCkIGoQNAKGhsbHS4BqWn6okkjRgxQhEREU7VksOHDztVVc6UnJwsSbrooov0ySefqKioyCmgPPjgg1q3bp127dqliy++2O15EFAQ8AgeAPAvMTExLl0kGxUVpbS0NFVVVemGG26wt1dVVem6665z+fsMw1BHR4dD2wMPPKB77rlHO3bsUHp6uuuDPw0BBaZE6AAA78vPz1dOTo7S09OVkZGh8vJyNTQ0KDc3V5K0YsUKNTU1aevWrZKkxx57TImJiZo4caKkk/dFefDBB3X77bfb97l+/XqtXr1azz77rMaNG2ev0AwdOlRDhw51eWwEFPgc4QMAzCE7O1tHjx7V2rVrZbValZqaqsrKSiUlJUmSrFarGhoa7P27u7u1YsUKHTp0SJGRkbrgggt03333acmSJfY+paWlstls+t73vufwXWe7FuZM3AcFA4bgASDY+PI+KEn339Pvi2Q/vHuVV8fqS1RQ4DICCADAVwgokET4AACYCwElBBA+AACBhoAS4AgfAIBgREAxOQIIACAUEVD8jAACAIAzAoqXEUAAAHAfAWUAEEIAABhYBBQXEUIAAPAdAsppCCEAAJhDyAaULybFKzLS81sKAwAA7wn39wAAAADOREABAACmQ0ABAACmQ0ABAACmQ0ABAACm4/eAUlpaquTkZFksFqWlpam6utql7V577TVFRkbqkksu8e4AAQCAz/k1oFRUVCgvL08rV65UXV2dZs6cqblz56qhoaHP7VpbW7Vw4ULNmjXLRyMFAAC+5HFAqa6u1g9/+ENlZGSoqalJkvTrX/9ae/fudXkfGzZs0KJFi7R48WKlpKSopKRECQkJKisr63O7JUuWaMGCBcrIyPB0+AAAwMQ8Cijbt2/XnDlzNHjwYNXV1amjo0OSdOzYMa1bt86lfdhsNtXW1iorK8uhPSsrSzU1Nb1u99RTT+n9999XYWGhS9/T0dGhtrY2hxcAADA3jwLKPffco8cff1xPPPGEBg0aZG+fMWOG3nrrLZf20dLSoq6uLsXFxTm0x8XFqbm5ucdtDh48qIKCAj3zzDOKjHTtJrjFxcWKjY21vxISElzaDgAA+I9HAeXAgQO67LLLnNpjYmL02WefubWvsLAwh/eGYTi1SVJXV5cWLFigNWvW6MILL3R5/ytWrFBra6v91djY6Nb4AACA73n0LJ74+Hi99957GjdunEP73r17df7557u0jxEjRigiIsKpWnL48GGnqop08vTRm2++qbq6Ot12222SpO7ubhmGocjISO3cuVPf+ta3nLaLjo5WdHS0izMDAABm4FEFZcmSJVq+fLlef/11hYWF6eOPP9Yzzzyju+66S0uXLnVpH1FRUUpLS1NVVZVDe1VVlWbMmOHUPyYmRm+//bbq6+vtr9zcXE2YMEH19fWaNm2aJ1MBAAAm5FEF5Wc/+5laW1t15ZVXqr29XZdddpmio6N111132asbrsjPz1dOTo7S09OVkZGh8vJyNTQ0KDc3V9LJ0zNNTU3aunWrwsPDlZqa6rD9qFGjZLFYnNoBAEBg8yigSNK9996rlStXav/+/eru7takSZM0dOhQt/aRnZ2to0ePau3atbJarUpNTVVlZaWSkpIkSVar9az3RAEAAMEnzDAMw9+D8KW2tjbFxsYqc1aRIiMt/h4OAC/754Qofw8BAayro13vbPq5WltbFRMT45XvOHVcSrr/HoVbPD8udbe368O7V3l1rL7kcQUFQOjgIA/A1wgogIkRDACEKgIK4CLCAgD4DgEFQYMAAQDBg4ACnyNIAADOhoACtxAuAAC+QEAJMQQMAEAgIKAEKIIGACCYEVBMgLABAIAjAooXEDgAAOgfAoobCB4AAPhGyAcUQgcAAOYTsgHls69FKSKacAIAgBmF+3sAAAAAZyKgAAAA0yGgAAAQwkpLS5WcnCyLxaK0tDRVV1f32vf555/XVVddpZEjRyomJkYZGRnasWNHr/2fe+45hYWF6frrr3d7XAQUAABCVEVFhfLy8rRy5UrV1dVp5syZmjt3rhoaGnrsv2fPHl111VWqrKxUbW2trrzySl1zzTWqq6tz6vvhhx/qrrvu0syZMz0aGwEFAIAQtWHDBi1atEiLFy9WSkqKSkpKlJCQoLKysh77l5SU6Gc/+5m+8Y1vaPz48Vq3bp3Gjx+v3/3udw79urq69IMf/EBr1qzR+eef79HYCCgAAASRtrY2h1dHR0eP/Ww2m2pra5WVleXQnpWVpZqaGpe+q7u7W8eOHdOwYcMc2teuXauRI0dq0aJFnk1CIfwzYwAAzGToB+GKiPa8btDVcXLbhIQEh/bCwkIVFRU59W9paVFXV5fi4uIc2uPi4tTc3OzSd/7yl7/UiRMnNG/ePHvba6+9ps2bN6u+vt69CZyBgAIAQBBpbGxUTEyM/X10dHSf/cPCwhzeG4bh1NaTbdu2qaioSC+99JJGjRolSTp27Jh++MMf6oknntCIESM8GP2/EFAAAAgiMTExDgGlNyNGjFBERIRTteTw4cNOVZUzVVRUaNGiRfrNb36j2bNn29vff/99/eMf/9A111xjb+vu7pYkRUZG6sCBA7rgggtcmgfXoAAAEIKioqKUlpamqqoqh/aqqirNmDGj1+22bdumm2++Wc8++6y+853vOHw2ceJEvf3226qvr7e/rr32Wl155ZWqr693Ov3UFyooAACEqPz8fOXk5Cg9PV0ZGRkqLy9XQ0ODcnNzJUkrVqxQU1OTtm7dKulkOFm4cKEeeughTZ8+3V59GTx4sGJjY2WxWJSamurwHeecc44kObWfDQEFAIAQlZ2draNHj2rt2rWyWq1KTU1VZWWlkpKSJElWq9XhniibNm1SZ2enli1bpmXLltnbb7rpJm3ZsmVAx0ZAAQAghC1dulRLly7t8bMzQ8fu3bvd3r+nwYVrUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOnwLB706dgF3f4eQsj46vv8fQEATiGgeBkHeLgqUP9bIVgB8IaQDSjHz+9WuCUwDwiAmZghWBGSgOATsgEFQPDwVUgiCAG+Q0ABABcNZBAi7AB9I6AAgB8MVNgh6CBYEVAAIIAN9OktAg/MgoACALDjNBbMgoACAPAKwg76g4ACADC9gQg7hJzAQkABAIQET0JOd7v/7/MTqoiTAADAdAgoAADAdPweUEpLS5WcnCyLxaK0tDRVV1f32nfv3r3KzMzU8OHDNXjwYE2cOFEbN2704WgBAIAv+PUalIqKCuXl5am0tFSZmZnatGmT5s6dq/379ysxMdGp/5AhQ3Tbbbfp4osv1pAhQ7R3714tWbJEQ4YM0S233OKHGQAAAG/wawVlw4YNWrRokRYvXqyUlBSVlJQoISFBZWVlPfafMmWKbrzxRk2ePFnjxo3TD3/4Q82ZM6fPqgsAAAg8fgsoNptNtbW1ysrKcmjPyspSTU2NS/uoq6tTTU2NLr/8cm8MEQAA+InfTvG0tLSoq6tLcXFxDu1xcXFqbm7uc9vzzjtPR44cUWdnp4qKirR48eJe+3Z0dKijo8P+vq2trX8DBwAAXuf3i2TDwsIc3huG4dR2purqar355pt6/PHHVVJSom3btvXat7i4WLGxsfZXQkLCgIwbAAB4j98qKCNGjFBERIRTteTw4cNOVZUzJScnS5IuuugiffLJJyoqKtKNN97YY98VK1YoPz/f/r6trY2QAgCAyfmtghIVFaW0tDRVVVU5tFdVVWnGjBku78cwDIdTOGeKjo5WTEyMwwsAAJibX39mnJ+fr5ycHKWnpysjI0Pl5eVqaGhQbm6upJPVj6amJm3dulWS9NhjjykxMVETJ06UdPK+KA8++KBuv/12v80BAAAMPL8GlOzsbB09elRr166V1WpVamqqKisrlZSUJEmyWq1qaGiw9+/u7taKFSt06NAhRUZG6oILLtB9992nJUuW+GsKAADAC8IMwzD8PQhfamtrU2xsrJLuv0fhFou/hwMAMLHu9nZ9ePcqtba2eu0SgVPHpclL1iki2vPjUldHu97Z9HOvjtWX/P4rHgAAgDMRUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkQUAAACGGlpaVKTk6WxWJRWlqaqqure+1rtVq1YMECTZgwQeHh4crLy+ux32effaZly5YpPj5eFotFKSkpqqysdGtcfr2TrD8NSWpTxFd6f4ZPIDp+KNbfQwAABJCKigrl5eWptLRUmZmZ2rRpk+bOnav9+/crMTHRqX9HR4dGjhyplStXauPGjT3u02az6aqrrtKoUaP029/+Vuedd54aGxv11a9+1a2xhWxACUZDk1t9+n0EIgAIbBs2bNCiRYu0ePFiSVJJSYl27NihsrIyFRcXO/UfN26cHnroIUnSk08+2eM+n3zySX366aeqqanRoEGDJMn+CBt3EFDgsYEMRIQdABgYbW1tDu+jo6MVHR3t1M9ms6m2tlYFBQUO7VlZWaqpqfH4+19++WVlZGRo2bJleumllzRy5EgtWLBAd999tyIiIlzeDwEFpjAQYYeQAyCQnfOeTZGRnl8a2tlpkyQlJCQ4tBcWFqqoqMipf0tLi7q6uhQXF+fQHhcXp+bmZo/H8cEHH+iVV17RD37wA1VWVurgwYNatmyZOjs79Ytf/MLl/RBQEDT6G3IIOACCQWNjo8PDAnuqnpwuLCzM4b1hGE5t7uju7taoUaNUXl6uiIgIpaWl6eOPP9YDDzxAQAE8QRUHQDCIiYlx6WnGI0aMUEREhFO15PDhw05VFXfEx8dr0KBBDqdzUlJS1NzcLJvNpqioKJf2Q0ABBhBVHACBIioqSmlpaaqqqtINN9xgb6+qqtJ1113n8X4zMzP17LPPqru7W+HhJ09Zvfvuu4qPj3c5nEgEFMBU+hNwCDcA3JWfn6+cnBylp6crIyND5eXlamhoUG5uriRpxYoVampq0tatW+3b1NfXS5KOHz+uI0eOqL6+XlFRUZo0aZIk6dZbb9Ujjzyi5cuX6/bbb9fBgwe1bt063XHHHW6NjYACBAnCDQB3ZWdn6+jRo1q7dq2sVqtSU1NVWVlp/1mw1WpVQ0ODwzZTpkyx/3Ntba2effZZJSUl6R//+Iekkxfp7ty5U3feeacuvvhijR07VsuXL9fdd9/t1tgIKAA8CjeEGiA4LF26VEuXLu3xsy1btji1GYZx1n1mZGToz3/+c7/GRUAB4BF3Qg1hBoC7CCgAvM6VMEOIAXA6AgoAU3C1IkOQAUIDAQVAQDlbkCHAAMGBgAIgqHA6CQgOBBQAIYcQA5gfAQUAesCpJMC/CCgA4AECDOBdBBQA8AICDNA/BBQA8AMCDNA3AgoAmBAX8iLUEVAAIEARYhDMCCgAEMQ4lYRARUABgBBGgIFZEVAAAL3iNBL8hYACAOgXqjDwBgIKAMCrqMLAEwQUAIDfUYXBmQgoAADTI8CEHgIKACDg9RVgCC+BKWQDyrfOO6jooYO0s2GCv4cCAPAiqi+BKWQDyilZiQdc6keQAYDg1FeAaftbtA9HgtOFfEBxlatBRiLMAECwGJLU5u8hhCwCihdQlQEAoH8IKH7kSpAhxAAAQhEBxeSoxgAAQhEBJUhQjQEABBMCSgghxAAAAgUBBQ7OFmIIMAAAXyCgwC0EGACALxBQMKAIMACAgUBAgU/1FWAILwCAUwgoMA3CCwDglHB/D6C0tFTJycmyWCxKS0tTdXV1r32ff/55XXXVVRo5cqRiYmKUkZGhHTt2+HC08JesxAO9vgAAwcevFZSKigrl5eWptLRUmZmZ2rRpk+bOnav9+/crMTHRqf+ePXt01VVXad26dTrnnHP01FNP6ZprrtHrr7+uKVOm+GEGMAMqLwAQfMIMwzD89eXTpk3T1KlTVVZWZm9LSUnR9ddfr+LiYpf2MXnyZGVnZ+sXv/iFS/3b2toUGxur2/beoOihgzwaN4ID4QXA2XR93qH989ertbVVMTExXvmOU8elzFlFioy0eLyfzs52vfb/irw6Vl/yWwXFZrOptrZWBQUFDu1ZWVmqqalxaR/d3d06duyYhg0b1mufjo4OdXR02N+3tfFkSpzUW+WF4AIA/ue3gNLS0qKuri7FxcU5tMfFxam5udmlffzyl7/UiRMnNG/evF77FBcXa82aNf0aK0ILp4wAwP/8/iuesLAwh/eGYTi19WTbtm0qKirSSy+9pFGjRvXab8WKFcrPz7e/b2trU0JCgucDRkij6gIAvuG3gDJixAhFREQ4VUsOHz7sVFU5U0VFhRYtWqTf/OY3mj17dp99o6OjFR0d3e/xAn0huADAwPJbQImKilJaWpqqqqp0ww032Nurqqp03XXX9brdtm3b9KMf/Ujbtm3Td77zHV8MFfAYwQUAPOPXUzz5+fnKyclRenq6MjIyVF5eroaGBuXm5ko6eXqmqalJW7dulXQynCxcuFAPPfSQpk+fbq++DB48WLGxsX6bB+AuggsA9M2vN2rLzs5WSUmJ1q5dq0suuUR79uxRZWWlkpKSJElWq1UNDQ32/ps2bVJnZ6eWLVum+Ph4+2v58uX+mgIwoLgZHQBfc+eGqZL06quvKi0tTRaLReeff74ef/xxpz4lJSWaMGGCBg8erISEBN15551qb293a1x+v0h26dKlWrp0aY+fbdmyxeH97t27vT8gwISouADwBndvmHro0CFdffXV+vGPf6ynn35ar732mpYuXaqRI0fq3//93yVJzzzzjAoKCvTkk09qxowZevfdd3XzzTdLkjZu3Ojy2PweUAB4juACoD82bNigRYsWafHixZJOVj527NihsrKyHm+Y+vjjjysxMVElJSWSTt5c9c0339SDDz5oDyj79u1TZmamFixYIEkaN26cbrzxRr3xxhtujY2AAgQhggsQus68IWlvv2b15Iap+/btU1ZWlkPbnDlztHnzZn355ZcaNGiQLr30Uj399NN644039M1vflMffPCBKisrddNNN7k1DwIKEEIILoB5Dd5vVWR4lMfbd3bbJMnpXl+FhYUqKipy6u/JDVObm5t77N/Z2amWlhbFx8dr/vz5OnLkiC699FIZhqHOzk7deuutTkHobAgoAAguQBBpbGx0eBbP2e4F5u4NU3vqf3r77t27de+996q0tFTTpk3Te++9p+XLlys+Pl6rV692eR4EFAC9IrgAgScmJsalhwV6csPU0aNH99g/MjJSw4cPlyStXr1aOTk59utaLrroIp04cUK33HKLVq5cqfBw135ATEAB4LaegguhBQgsntwwNSMjQ7/73e8c2nbu3Kn09HQNGjRIkvT55587hZCIiAgZhmGvtriCgAJgQFBtAQKPuzdMzc3N1aOPPqr8/Hz9+Mc/1r59+7R582Zt27bNvs9rrrlGGzZs0JQpU+yneFavXq1rr71WERERLo+NgALAq6i2AOaVnZ2to0ePau3atbJarUpNTe3zhqnJycmqrKzUnXfeqccee0xjxozRww8/bP+JsSStWrVKYWFhWrVqlZqamjRy5Ehdc801uvfee90aW5jhTr0lCLS1tSk2Nla37b1B0UMH+Xs4AP4PoQVm1PV5h/bPX6/W1laXruvwxKnj0uz4Jf3+Fc8u6yavjtWXqKAAMAUqLQBOR0ABYFqEFiB0EVAABJQzQwuBBQhOBBQAAY0qCxCcCCgAgg5VFiDwhWxAmX/O6xr61X/dSObX/5zhx9EA8CYCCxB4QjagnCnn3J6f3EhwAYIPp4UA8yOgnEVPwYXQAgQfqiyAuRBQPHBmaCGwAMHn9MBCWAF8j4AyAKiyAMGN6grgewQULyG0AMGLwAJ4HwHFhzg1BAQnTgcBA4+A4kdUWYDgQ3UFGBgEFJOhygIEFwIL4BkCiskRWIDgQmABXENACTAEFiC4cP0K0DMCSoAjsADBg+oK8C8ElCBDYAGCB9UVhDICSpAjsADBgeoKQg0BJcScHlgIK0DgorqCYEdACWFUV4DgQHUFwYiAAjsCCxAcqK4gGBBQ0CtOBwGBj+oKAhUBBS6hugIEBwILAgUBBR6hugIEB04HwawIKOg3qitAcKC6AjMhoGDAUV0BggPVFfgTAQVeRVgBggPVFfgaAQU+w6kgIHhQXYG3EVDgN1RXgOBAdQXeQECBKRBWgOBBdQUDgYAC0yGsAMGDsAJPEVBgaoQVIHhwKgjuIKAgYHCRLRBcqK6gLwQUBCyqK0DwIKzgTAQUBAXCChA8OBUEiYCCIERYAYIL1ZXQREBBUCOsAMGFsBI6CCgIGYQVILgQVoJbuL8HAPhDzrk1Tr8KAhC4shIPOF27gsBGBQUhjaoKEFyoqgQPAgrwfwgrQHAhrAQ2v5/iKS0tVXJysiwWi9LS0lRdXd1rX6vVqgULFmjChAkKDw9XXl6e7waKkMIpICC4cAqod+4chyXp1VdfVVpamiwWi84//3w9/vjjTn22b9+uSZMmKTo6WpMmTdILL7zg9rj8GlAqKiqUl5enlStXqq6uTjNnztTcuXPV0NDQY/+Ojg6NHDlSK1eu1Ne//nUfjxah6FRQIawAweFUUCGsnOTucfjQoUO6+uqrNXPmTNXV1ennP/+57rjjDm3fvt3eZ9++fcrOzlZOTo7+8pe/KCcnR/PmzdPrr7/u1tjCDMMw+jW7fpg2bZqmTp2qsrIye1tKSoquv/56FRcX97ntFVdcoUsuuUQlJSVufWdbW5tiY2O1969jNPSrfi8gIQBx+gcIPr2dAur6vEP7569Xa2urYmJivPLdp45Ls+OXKDI8yuP9dHbbtMu6ya2xunscvvvuu/Xyyy/rb3/7m70tNzdXf/nLX7Rv3z5JUnZ2ttra2vTHP/7R3ufb3/62zj33XG3bts3l+fjtGhSbzaba2loVFBQ4tGdlZammZuD+ttrR0aGOjg77+9bWVknSiePdA/YdCC03RO6VJD332TQ/jwTAQLl82F/1ykfjndq7Pj95/PDF3+U7DZvUj0NTp2GTdDLwnC46OlrR0dFO/T05Du/bt09ZWVkObXPmzNHmzZv15ZdfatCgQdq3b5/uvPNOpz7uFhT8FlBaWlrU1dWluLg4h/a4uDg1NzcP2PcUFxdrzZo1Tu1zpg/cdyBUuX9OFUBgOnr0qGJjY72y76ioKI0ePVq7m5/q976GDh2qhIQEh7bCwkIVFRU59fXkONzc3Nxj/87OTrW0tCg+Pr7XPu4e2/3+K56wsDCH94ZhOLX1x4oVK5Sfn29//9lnnykpKUkNDQ1e+4/NF9ra2pSQkKDGxkavlR19IRjmEQxzkIJjHsEwB4l5mElra6sSExM1bNgwr32HxWLRoUOHZLPZ+r2vno6hPVVPTufucbin/me2D8Sx3W8BZcSIEYqIiHBKVIcPH3ZKXv3RW2krNjY2YP+HOV1MTAzzMIlgmIMUHPMIhjlIzMNMwsO9e82ixWKRxWLx6necyZPj8OjRo3vsHxkZqeHDh/fZx91ju9+uEo2KilJaWpqqqqoc2quqqjRjBhchAgDgTZ4chzMyMpz679y5U+np6Ro0aFCffdw9tvv1FE9+fr5ycnKUnp6ujIwMlZeXq6GhQbm5uZJOnp5pamrS1q1b7dvU19dLko4fP64jR46ovr5eUVFRmjRpkj+mAABAwHL3OJybm6tHH31U+fn5+vGPf6x9+/Zp8+bNDr/OWb58uS677DLdf//9uu666/TSSy9p165d2rt3r3uDM/zsscceM5KSkoyoqChj6tSpxquvvmr/7KabbjIuv/xyh/6SnF5JSUkuf197e7tRWFhotLe3D9AM/IN5mEcwzMEwgmMewTAHw2AeZhIMczgbd4/Du3fvNqZMmWJERUUZ48aNM8rKypz2+Zvf/MaYMGGCMWjQIGPixInG9u3b3R6XX++DAgAA0BPuVAYAAEyHgAIAAEyHgAIAAEyHgAIAAEwnoAPKnj17dM0112jMmDEKCwvTiy++6PC5YRgqKirSmDFjNHjwYF1xxRV65513zrrfgXhMtDu8MY8tW7YoLCzM6dXe3u63eTz//POaM2eORowYobCwMPtPxs/Gl+vhjTmYbS2+/PJL3X333brooos0ZMgQjRkzRgsXLtTHH3981v2aZS08nYPZ1kKSioqKNHHiRA0ZMkTnnnuuZs+e7dJTX83255Qn8/D1epxtDqdbsmSJwsLCXHp+jK/XIlQEdEA5ceKEvv71r+vRRx/t8fP169drw4YNevTRR/U///M/Gj16tK666iodO3as130O1GOi3eGNeUgn795otVodXt68U+HZ5nHixAllZmbqvvvuc3mfvl4Pb8xBMtdafP7553rrrbe0evVqvfXWW3r++ef17rvv6tprr+1zn2ZaC0/nIJlrLSTpwgsv1KOPPqq3335be/fu1bhx45SVlaUjR470uk8z/jnlyTwk367H2eZwyosvvqjXX39dY8aMOes+/bEWIcPtHyablCTjhRdesL/v7u42Ro8ebdx33332tvb2diM2NtZ4/PHHe93PvHnzjG9/+9sObXPmzDHmz58/4GPuyUDN46mnnjJiY2O9ONK+nTmP0x06dMiQZNTV1Z11P/5cj4Gag5nX4pQ33njDkGR8+OGHvfYx61qc4socAmEtWltbDUnGrl27eu1jtj+neuLKPPy5Hr3N4aOPPjLGjh1r/PWvfzWSkpKMjRs39rkff69FMAvoCkpfDh06pObmZofHQkdHR+vyyy/v9THSUu+Pku5rG2/ydB7SybvtJiUl6bzzztO//du/qa6uztvDHXBmWw9PmX0tWltbFRYWpnPOOafXPmZfC1fmIJl7LWw2m8rLyxUbG6uvf/3rvfYz+1q4Og/JXOvR3d2tnJwc/fSnP9XkyZNd2sbsaxHIgjagnHpQkbuPfB6ox0QPFE/nMXHiRG3ZskUvv/yytm3bJovFoszMTB08eNCr4x1oZlsPT5h9Ldrb21VQUKAFCxb0+UA3M6+Fq3Mw61r8/ve/19ChQ2WxWLRx40ZVVVVpxIgRvfY361q4Ow+zrcf999+vyMhI3XHHHS5vY9a1CAZ+fRaPL3jyyOeBeEz0QHN3TNOnT9f06dPt7zMzMzV16lQ98sgjevjhh702Tm8w43q4w8xr8eWXX2r+/Pnq7u5WaWnpWfubcS3cmYNZ1+LKK69UfX29Wlpa9MQTT9ivYRg1alSv25hxLdydh5nWo7a2Vg899JDeeustt/89mnEtgkHQVlBGjx4tSW4/8nmgHhM9UDydx5nCw8P1jW98w+9/U3SX2dZjIJhlLb788kvNmzdPhw4dUlVVVZ+VB8mca+HuHM5klrUYMmSIvva1r2n69OnavHmzIiMjtXnz5l77m3EtJPfncSZ/rkd1dbUOHz6sxMRERUZGKjIyUh9++KF+8pOfaNy4cb1uZ9a1CAZBG1CSk5M1evRoh0c+22w2vfrqq30+8nmgHhM9UDydx5kMw1B9fb3i4+O9MUyvMdt6DAQzrMWpA/vBgwe1a9cuDR8+/KzbmG0tPJnDmcywFj0xDEMdHR29fm62tejN2ebRU39/rUdOTo7+93//V/X19fbXmDFj9NOf/lQ7duzodbtAWYtAFNCneI4fP6733nvP/v7QoUOqr6/XsGHDlJiYqLy8PK1bt07jx4/X+PHjtW7dOn3lK1/RggUL7NssXLhQY8eOVXFxsaQBfEy0n+exZs0aTZ8+XePHj1dbW5sefvhh1dfX67HHHvPbPD799FM1NDTY71Vx4MABSSf/BnKqUuTv9fDGHMy2FmPGjNH3vvc9vfXWW/r973+vrq4u+98Ahw0bpqioqB7nYaa18HQOZluL4cOH695779W1116r+Ph4HT16VKWlpfroo4/0/e9/376Nv9fCW/Pw9Xqc7f/vM0PuoEGDNHr0aE2YMKHXOfhjLUKGX347NED+9Kc/GZKcXjfddJNhGCd/oltYWGiMHj3aiI6ONi677DLj7bffdtjH5Zdfbu9/ykA8Jtrf88jLyzMSExONqKgoY+TIkUZWVpZRU1Pj13k89dRTPX5eWFjY6zwMw7fr4Y05mG0tTv1EuqfXn/70p17nYRjmWQtP52C2tfjiiy+MG264wRgzZowRFRVlxMfHG9dee63xxhtvOOzD32vhrXn4ej3O9v/3mXr6mbEZ1iJUhBmGYfQz4wAAAAyooL0GBQAABC4CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMJ2AflgggP674oorlJqaKkl6+umnFRERoVtvvVX/+Z//qbCwMD+PDkCoooICQP/1X/+lyMhIvf7663r44Ye1ceNG/epXv/L3sACEMB4WCIS4K664QocPH9Y777xjr5gUFBTo5Zdf1v79+/08OgChigoKAE2fPt3hdE5GRoYOHjyorq4uP44KQCgjoAAAANMhoADQn//8Z6f348ePV0REhJ9GBCDUEVAAqLGxUfn5+Tpw4IC2bdumRx55RMuXL/f3sACEMH5mDEALFy7UF198oW9+85uKiIjQ7bffrltuucXfwwIQwggoADRo0CCVlJSorKzM30MBAEmc4gEAACZEQAEAAKbDjdoAAIDpUEEBAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACm8/8B2S5Pu05O+TQAAAAASUVORK5CYII=",
                        "text/plain": [
                            "<Figure size 640x480 with 2 Axes>"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                }
            ],
            "source": [
                "p = np.linspace(10.0, 14.0)\n",
                "e = np.linspace(0.1, 0.7)\n",
                "\n",
                "p_all, e_all = np.asarray([temp.ravel() for temp in np.meshgrid(p, e)])\n",
                "\n",
                "print(\"Total trajectory points:\", p_all.shape[0])\n",
                "teuk_modes = amp(0.0, p_all, e_all, np.ones_like(p_all))  # a, p, e, xI\n",
                "\n",
                "# look at the contours of the (2,2,0) mode\n",
                "cb = plt.contourf(\n",
                "    p,\n",
                "    e,\n",
                "    np.abs(teuk_modes[:, amp.special_index_map[(2, 2, 0, 0)]].reshape(len(p), len(e))),\n",
                ")\n",
                "plt.colorbar(cb)\n",
                "plt.xlabel(\"p\")\n",
                "plt.ylabel(\"e\")\n",
                "plt.show()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": []
        }
    ],
    "metadata": {
        "kernelspec": {
            "display_name": "few-env",
            "language": "python",
            "name": "python3"
        },
        "language_info": {
            "codemirror_mode": {
                "name": "ipython",
                "version": 3
            },
            "file_extension": ".py",
            "mimetype": "text/x-python",
            "name": "python",
            "nbconvert_exporter": "python",
            "pygments_lexer": "ipython3",
            "version": "3.12.11"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 2
}