{
    "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, n)$ in the set $\\ell \\in [2, 10]$, $m \\in [-\\ell, \\ell]$, $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 -n} = (-1)^\\ell A_{\\ell m 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": 4,
            "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, 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): array([0.21239204-0.06664877j]),\n",
                            "  (7, -3, 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),(7,-3,1)]\n",
                "specific_amplitudes = kerr_amp_spline(a, p, e, xI, specific_modes=specific_modes)\n",
                "specific_amplitudes, specific_amplitudes[(2,2,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": 5,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "text/plain": [
                            "True"
                        ]
                    },
                    "execution_count": 5,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "# Extracting the (2,2,0) mode\n",
                "mode_location = kerr_amp_spline.special_index_map[(2,2,0)]\n",
                "np.allclose(\n",
                "    specific_amplitudes[(2,2,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": 6,
            "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 apex_of_uwyz, 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 = apex_of_uwyz(\n",
                "    u_grid, w_grid, np.ones_like(u_grid) * y, np.ones_like(u_grid) * z\n",
                ")\n",
                "\n",
                "print(\"Total trajectory points:\", w_grid.shape[0])"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 7,
            "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), (7, -3, 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)].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, 1)].reshape(uv.size, wv.size),\n",
                "        np.conj(teuk_all_amps[:, :, inds[1]]),\n",
                "    )\n",
                ")\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 8,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAF4CAYAAAC/wIoGAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAATjBJREFUeJzt3Qt8VOWdN/B/LiQBIUGMXI1GUECLBAtLCugKNcgii7K+VooWIgKuCl1K1i5E5CZIrGjEdYNUMOKNkkotuoWiNoUqBckS5F3qCshFEhEC2RYCURJI5v38H98ZZ5KZZC7n8lx+38/n1M4wJ3POzJzn/OY/z/OcOI/H4yEAAAAAAE3Eu70BAAAAAABWQsAFAAAAAK0g4AIAAACAVhBwAQAAAEArCLgAAAAAoBUEXAAAAADQCgIuAAAAAGgFARcAAAAAtIKACwAAAABaQcAFAAAAAK0g4AIAAACALT788EMaO3Ysde/eneLi4mjDhg2trlNXV0dz586lq666ipKTkykzM5OKi4sjet7EGLYZAAAAACCk2tpaysrKogceeIDuuusuCsc999xDVVVV9PLLL9M111xDx48fp8bGRooEAi4AAAAA2GL06NFiCdfmzZvpT3/6Ex0+fJg6deok7uMKbqSMC7j8DeCrr76iDh06iFI5AIDVPB4PnT17VvwkFx+vX08wtKMAarcz58+fp/r6+qift+lxz90IeLHCu+++S4MGDaKnn36aXn/9dbrkkkvojjvuoMWLF1Pbtm3D/jvGBVxulDMyMtzeDAAwQGVlJV1xxRWkG7SjAOq2M+fPn6crr7yETp2K7Cd/r/bt29O5c+cC7luwYAEtXLiQrMCV223btlFKSgr99re/perqanrkkUfof//3f+mVV14J++8YF3C54sCGp0+ixr7NS95JB75qcf363t3Dfq7W/la4InlOK1i13QCmuthYT1urX/O1N7rx7hefWFNTU93eHADt/J9r/9W2dqa+vl6E2607O1P79pH9AnPunIeGZ59sduxbVb31/kLEFeI333yT0tLSxH2FhYV0991304oVK8Ku4hoXcL1l9cT4JGpMTGn273x/S4KtE0prfysc9X2vcPRNStr3JZEF2w0A37U3uu4Xn+AQcAFiN7rb9KjzQ7TtTPv2cdS+Q6RdqBptP/a7detGPXr08IVbdt1114muEV9++SVde+21Yf0d4wKuo0HRgnALAAAAegda+M6wYcPorbfeEt0guDsEO3DggOhnHElXDGMDLv/snxhhKNU9cFoRygEAACCQyYH23LlzdPDgQd/tI0eO0J49e8QMCVdeeSXl5+fTsWPH6LXXXhP/fu+994oBZZMnT6ZFixaJPrg///nPxTRjGGSmAd3DNAAAgM5MDrX+du3aRSNGjPDdzsvLE//Nzc2lNWvWiDluKyoqfP/OVdsPPviAfvrTn4rZFC677DIxL+6SJUsoEgi4NgROFSuhKm4zAACALBBogxs+fLjoPxsKh9ym+vbtK0JuLBBwJQx4qN4CAADIDYFWbgi4knEj3MoU7gEAAGSEQKsWBFyLISwCAADoAaFWXQi4ElVVUb0FAABwj0qB9jefP0tpaavd3gxpIeAi5AEAABhLpVAL4UPAlSQoo3oLAABgPwRaMyDgGjqjAcItAACYAqHWPAi4EtA1RAMAALgBgRYQcF2uiKJrgvvCfQ/wugEAyAuhFvwh4LYAlVV9WPFehvobCL4AAO5AqIVQjA+4boYTVG/1eG39n8uU1xcAwC0ItRAO4wOuFRBq5CBDxR1hFwDAWgi0EA2jAy6qt+qTIdS2tm06vu4AAHZCqIVYGR1wVQ1OoNb7g6ALANA6hFqwEgJujKIJLajemhFsg227Lu8DAIAVEGqj8/vjRVRTU+P2ZkjN2ICbdOArovgktzcDNA+1TaGaCwCmQ6gFJxgbcN0KVKjemhlsm0I1FwBMglALTkPAjQECin10DrdeCLkAoDOEWnATAq6DUL1tnQnB1h9CLgDoBKEWZIGAqzGVgpNpwdYfQi4AqAyhFmSEgBtl0Io0kJgc4FqD1wYhFwDUg2ALMkPA1ZQKYQnBNhBCLgDIDqEWVIGA6wAEuebwmgSHkAsAskGoBRXFu70BKpI9gMi8fRzgEG4BrFdUVESZmZmUkpJC2dnZVFZW1uLjly9fTn369KG2bdtSRkYGzZo1i86fP+/Y9oL8oda7AKgIFVw/dgQvhLnv4LUID6q4EKmSkhLKy8ujlStXinDL4XXUqFG0f/9+6ty5c7PHr127lubMmUPFxcU0dOhQOnDgAN1///0UFxdHhYWFruwDyAGBFnSBgKsZWYMRwm1kEHIhEhxKp02bRpMnTxa3Oehu3LhRBFgOsk1t376dhg0bRvfee6+4zZXfCRMm0M6dOx3fdnAfQi3oCF0UIoTQERl0SQCwV319PZWXl1NOTo7vvvj4eHF7x44dQdfhqi2v4+3GcPjwYdq0aRPdfvvtQR9fV1cnrnvvv4D60AVBTb8/XuT2JigBFVyNuifIFr4RbAHsV11dTQ0NDdSlS5eA+/n2vn37gq7DlVte76abbiKPx0MXL16khx56iB577LGgjy8oKKBFixbZsv3gLARaMAUquGALhNvY4TUEu2zdupWWLl1KK1asoN27d9Pbb78tujQsXrw46OPz8/PpzJkzvqWystLxbYboYcAYuOnDDz+ksWPHUvfu3UU//w0bNoS97p///GdKTEykAQMGRPy8qODaVCE1uXqLYAbgnPT0dEpISKCqqqqA+/l2165dg64zb948mjhxIk2dOlXcvuGGG6i2tpYefPBBmjt3ruji4C85OVksoBYEWpABty1ZWVn0wAMP0F133RX2eqdPn6ZJkybRrbfe2qx9CwcCLlgGwRbAeUlJSTRw4EAqLS2lcePGifsaGxvF7RkzZgRd5+uvv24WYjkkM+6yAOpCqAUn1DTph9/Sl+DRo0eLJVLcbYq7U3HbFEnV1wsB14ZgZmL1FuHWHphNAcLBU4Tl5ubSoEGDaPDgwWKaMK6aeGdV4CpIjx49RF9axj8X8swLN954o5hW7ODBg6Kqy/d7gy6oBcEWIrXhbBaleNpEtM75cxeI6H0xd7a/BQsW0MKFCy3btldeeUUMfn3jjTdoyZIlUf0NBFyIGcItgLvGjx9Pp06dovnz59OJEydEf7XNmzf7Bp5VVFQEVGwff/xx0ReO/3vs2DG6/PLLRbh98sknXdwLiAaCLbihsrKSUlNTfbet7ML0+eefi+kNP/roI9H/NloIuGGStYrm9nYh3ALIgbsjhOqSwIPK/PFJgysuvIB6EGrBbampqQEB1yo8Iwx3S+BZW3r37h3T30LAtZhJgc+kfQUAcBuCLeju7NmztGvXLvrkk098X9h5TAGPDeAv5u+//z798Ic/DOtvGR9wVQ5pblZvVX7dAABUgmALpkhNTaW9e/cG3MfTGf7xj3+k9evX09VXXx323zI+4FoZJE0JfabsJwCAWxBqQRfnzp0TA1m9jhw5Qnv27KFOnTrRlVdeKebZ5rEAr732mhgr0K9fv4D1O3fuTCkpKc3ubw0CrqLcqt4i3AIA2AfBFnS7TO+uXbtoxIgRAbO+MJ75Zc2aNXT8+HExENZqRgdcK8Oak8EP4RYAQC8ItqCr4cOHtzi/NofclvD0Y9FMQWZ0wIXwIdwCAFgPwRbAHsYG3Pre3cPaeben4ZJhexBuAQCshWALYC9jA66VdA6AOu8bAICTEGoBnIOAqxCnq7cItwAAsUOwBXAeAm6MgVLXEKjrfqlGti4yABA+BFsA9yDgKsLJoINwCwAQPQRbAPch4EIAhFsAgOgg2ALIAwFXgTDoVPUW4RYAIHIItgDyQcANwbS+jwi38jHtMwigGgRbAHkh4EoOIQcAQC4ItuAGFS/T6yYE3CjpVPHUaV90gS82APJBsAVQBwKu4RBuAQBahmALoJ54tzegqKiIMjMzKSUlhbKzs6msrKzFxy9fvpz69OlDbdu2pYyMDJo1axadP3/e0eqZLoPLEG4BAFqGcAugJlcruCUlJZSXl0crV64U4ZbD66hRo2j//v3UuXPnZo9fu3YtzZkzh4qLi2no0KF04MABuv/++ykuLo4KCwtd2QdVIdzKC90TANyHYAugNlcDLofSadOm0eTJk8VtDrobN24UAZaDbFPbt2+nYcOG0b333ituc+V3woQJtHPnzpDPUVdXJxavmpoaMj3kINzKC+EWwF0ItgB6cK2LQn19PZWXl1NOTs53GxMfL27v2LEj6DpcteV1vN0YDh8+TJs2baLbb7895PMUFBRQWlqab+FuDTFtN8IhAICWwRbhFkAfrlVwq6urqaGhgbp06RJwP9/et29f0HW4csvr3XTTTeTxeOjixYv00EMP0WOPPRbyefLz80U3CP8KbkshV4YKGqq3ZpLhswdgGoRaAD25PsgsElu3bqWlS5fSihUraPfu3fT222+LLg2LFy8OuU5ycjKlpqYGLKZCuAUA+A7CLYC+XAu46enplJCQQFVVVQH38+2uXbsGXWfevHk0ceJEmjp1Kt1www30T//0TyLwcjeExsZGLQKiXVU8hFu5oXoLTs5IM3z4cDE4t+kyZswYMgG6IwDoz7WAm5SURAMHDqTS0lLffRxS+faQIUOCrvP111+Lfrr+OCQz7rIAwSHcyg3hFqyakWbBggXi162srCwxI83JkyeDPp5//Tp+/Lhv+ctf/iLa0h/96EekMwRbAHO4OosCN8i5ubk0aNAgGjx4sJgmrLa21jerwqRJk6hHjx6iQsvGjh0rZl648cYbRYXi4MGDoqrL93uDrspBw+3nBwA1RTojTadOnQJur1u3jtq1a6dtwEWoBTCPqwF3/PjxdOrUKZo/fz6dOHGCBgwYQJs3b/YNPKuoqAio2D7++OPiZzT+77Fjx+jyyy8X4fbJJ5+0fVtVrYKqut2mwJcasGpGGh5QG+6MNE29/PLL9OMf/5guueQSraZbZAi3oIPfHy9yexOU4/qlemfMmCGWUIPK/CUmJoqf4HjRjR1BB+FWbgi34NaMNP64ry53UeCQGwr/irZo0SJSCYItgNmUmkUBwodwCwDh4GDLg3a5m1goXB0+c+aMb6msrCRZoZ8tAEhRwZWFm9U0VPLMg/cc3JyRxovHPHD/2yeeeKLFx/F0i7zIzu5gG6xwgGMZQE4IuBpWQ1XbXtPghAh2zUgzbty4gBlpQnX/8nrrrbdE39qf/OQnpDK7gm04bWnTx+D4BpADAq7LrG4MEW7lhpMfyDAjjX/3BA7Fl112GanK6nAbaxvqXR/HOoC70AcXwCE44YGdM9I888wzYkYano1mz549zWak4flu/e3fv5+2bdtGU6ZMIRXZ0dfWygIBig0A3/rwww/FjFfdu3cXM2Ft2LCBWsLzdI8cOVLMlMVXn+VrI7z33nsUKVRwWwkeKjVSKm2raRBuQaYZaVifPn2UvUCOzMG26d/FsQ+mq62tFRefeeCBB+iuu+4KKxBzwOUr1Xbs2JFeeeUVEZB37twproMQLgRcF1nZ8CHcygsnOACzgm3T50AbACYbPXq0WMLFXaz8cdB955136D//8z8RcMFeZ3qFN5o67dB3E8MDAJgWbgF0VdPkYi92zrTCg2bPnj3b7AqMrUHAdakBVLF6G26w9X+86SEXlRuA2Nw64tuBcUkKh1tUcUFGW6p6U+K5yM7rF2v5nP4+ZWRkBNzPF+BauHAh2YHHF5w7d47uueeeiNYzPuCq3ug40VBHGmybrmtqyFX9swUgS7i1Eiq3oBoZL9NbWVkpBoB52VW9Xbt2rbiKIndR6Ny5c0TrGh9w3aBK8Ikl2JoeclV5jwFUCLZWHU8ItwDW4HDrH3DtwBehmTp1qpivOycnJ+L1jQ64qocQuxprq4KtqVT/XAHoVrVlCLcA6vjVr34lZl3gkDtmzJio/obRAdeNxlD2SgTCbWwQbgGip3O4RT9cMNW5c+fo4MGDvttHjhwRc3XzoLErr7yS8vPz6dixY/Taa6/5uiXwhWuef/55ys7OphMnToj727ZtS2lpaWE/r7EXekg68JXbmyAdO8OtCcEZJy+A6GdI0DncAphs165dYnov7xRffOVF/v98YRrGF6Hhi9F4vfTSS3Tx4kWaPn06devWzbfMnDkzoudFBVdBVjfYJoRPuyHcAkSHg22STccWwi2A+4YPH97iRWXWrFnT6oVpomFsBVfV7gkIt/JBuAWIjl1VW4ZwC2A2BFyDIdzGDuEWQD4ItwCAgOsQ2aq3CLexQ7gFgHChvQBwFgKuIlQPt7rNg4uTFYD97Vs0xxmqtwDAEHAdaBxlCkOo3Or1fgKoClcpAwA7IeAqwKpGG+E2dgi3AGbNDQ4AasI0YYZAuI0Ngi0AADjt98eL3N4EZaGCK3kwsqIqgXAbG4RbALmPO9mrt2hDAJyHgCtxI6lLuFV5gBlOTAByk63dBgA5IOBqHI5kCLcqc/v9AzABAioA2AEBV9NGH+E2Ngi3APJTIRyjLQFwBwKuQg2liuFWxe4JOCEBAACoDQFXwpAUS9iWKdyqCOEWQI3jUKeiBABYDwFXIwi3sUG4BXCW7sec7vsHIDMEXE2qtzKGW5W6J+BEBKAOVG8BoDUIuBI1lrJsh2kQbkEHRUVFlJmZSSkpKZSdnU1lZWUtPv706dM0ffp06tatGyUnJ1Pv3r1p06ZN5CSd2zy0KwDuwpXMNIDqbXRwAgJdlJSUUF5eHq1cuVKE2+XLl9OoUaNo//791Llz52aPr6+vp5EjR4p/W79+PfXo0YOOHj1KHTt2dGR7bx1REPVxqXMoBgDrIOBKEpp06pqgAoRb0ElhYSFNmzaNJk+eLG5z0N24cSMVFxfTnDlzmj2e7//rX/9K27dvpzZt2oj7uPobSl1dnVi8ampqbNkPXaB9ASvgMr2xQRcFhSHcRgcnH9AJV2PLy8spJyfHd198fLy4vWPHjqDrvPvuuzRkyBDRRaFLly7Ur18/Wrp0KTU0NAR9fEFBAaWlpfmWjIwMcgOqtwAQLuMDrpUNppPVW5nDrczdExBuQTfV1dUimHJQ9ce3T5w4EXSdw4cPi64JvB73u503bx49++yztGTJkqCPz8/PpzNnzviWyspKW/ZFB2hjAOSALgpgDJx4AL7V2Ngo+t++9NJLlJCQQAMHDqRjx47RsmXLaMGCBc0ez4PQeHETqrcAEAkEXJeheusMhFvQVXp6ugipVVVVAffz7a5duwZdh2dO4L63vJ7XddddJyq+3OUhKSmJ3KD6car69gPoxPguCqo1bDKHW1nhpAM64zDKFdjS0tKACi3f5n62wQwbNowOHjwoHud14MABEXzdCrctQfUWACJldMB1u9F0+/lNqN4i3IIJeIqwVatW0auvvkqfffYZPfzww1RbW+ubVWHSpEmiH60X/zvPojBz5kwRbHnGBR5kxoPOIDpoawDkgi4KCkH1NjI44YApxo8fT6dOnaL58+eLbgYDBgygzZs3+waeVVRUiJkVvHgWhPfee49mzZpF/fv3F/PgctidPXs2yUa3QgAAOMPYgFvfu7tlOx9NkIq00ZY93MpUvUWwBRPNmDFDLMFs3bq12X3cfeHjjz8mWah83Kq87QC6MrqLgioQbsOHEw0AOAltDoCcEHAVqN5CeHCiAZDf6G7h9/NFWwmgvg8//JDGjh1L3bt3p7i4ONqwYUOr6/CvTt///vfF9ITXXHMNrVmzJuLnRcCVHKq34UG4BVCDTqEV7Q5A63jAa1ZWFhUVhXfp4SNHjtCYMWNoxIgRtGfPHvrZz35GU6dOFeMGImFsH1wVGnfZw60scJIB0I9OQRggUr8/Hl4YVMHo0aPFEq6VK1fS1VdfLa6u6J2je9u2bfTcc8/RqFGjwv47qODGwPRgJUP11vT3AEAHKh7HKm4zgFVqamoClro66/LAjh07KCcnJ+A+DrZ8fyRQwXUQqrfWwgkGANyAtgd0UPFVOsW3TYloncZvzvumGvTHl/heuHChJdvFUx16pzj04tscpL/55htq27ZtWH8HAVdCKoRbt6u3OMEA6AvdEwDkVllZSampqb7bPBhMNgi4DgUsnRpshFsAMBXaHwAS4dY/4Fqpa9euVFVVFXAf3+bnC7d6y9AHVzIqVG/dhJMLALgF7Q+A/fgiNKWlpQH3ffDBB+L+SLgecHnaiMzMTEpJSaHs7GwqKytr8fGnT58W10vv1q2bKIn37t2bNm3aRDJD9dYaOLkA6Kfpca1TewkAROfOnRPTffHinQaM/z9fQpzl5+fTpEmTfI9/6KGH6PDhw/Rv//ZvtG/fPlqxYgX9+te/FpcWV6aLQklJCeXl5YkpITjcLl++XIyU279/P3Xu3LnZ4+vr62nkyJHi39avXy+un3706FHq2LGjFkEL1dvQEG4BwE1ogwCis2vXLjGnrRfnPpabmysu4HD8+HFf2GU8RdjGjRtFoH3++efpiiuuoNWrV0c0RZjrAbewsJCmTZtGkydPFrc56PJOFRcX05w5c5o9nu//61//Stu3b6c2bdqI+7j6K7NwqxEqhFs3qrc4qQCA29AOAURv+PDh5PF4Qv57sKuU8TqffPJJDM/qYhcFrsaWl5cHzHUWHx8vboea6+zdd98VfTC4iwJPGdGvXz9aunQpNTQ0hHwenput6XxtEDmEWwCwG7onAIBVXAu41dXVIpgGm+uM50ALhvtkcNcEXo/73c6bN09c6WLJkiUhn6egoIDS0tJ8S9O52+wMXTpVb52GcAugn1tHFJBq0BYBqMn1QWaRaGxsFP1vX3rpJRo4cCCNHz+e5s6dK7o2hMKdl8+cOeNbeO42kLt6ixMKAMgAbRGAulzrg5uenk4JCQlB5zrjOdCC4ZkTuO8tr+fF1yjmii93eUhKSmq2Ds+04MYExKjeRgcnFAAzj3d0TwD41u+PF7m9CVpwrYLLYZSrsP5znXGFlm+Hmuts2LBhdPDgQfE4rwMHDojgGyzcWr7NFocvFcKtk9VbhFsAkAXaIwC1udpFgaeKWLVqFb366qv02Wef0cMPP0y1tbW+WRV4XjTuYuDF/86zKMycOVMEW55xgQeZ8aAzsB7CLQCYCO0RgPpcnSaM+9CeOnWK5s+fL7oZDBgwgDZv3uwbeMbzovHMCl48QOy9994Tc6P1799fzIPLYXf27Nkkk3B+alOheusUnEwAAABAm4DLZsyYIZZgtm7d2uw+7r7w8ccfk9NMC2FOVW9Ne10BoDmZ+t+iTQLQg1KzKKhAh+otwi0AmAhtEoA+EHDBFTiRAJgNbQAAaN1FQaeGGNXb8ODEBgDhtJtN20s72ye0SwB6QQUXHIWTCIA9ioqKKDMzk1JSUig7O5vKyspavPZ7XFxcwMLrycTJYgDaJQD9IOA6yPTqLU4iAPYoKSkR0y4uWLCAdu/eTVlZWTRq1Cg6efJkyHVSU1Pp+PHjvuXo0aNkIrRLAHpCwHWwe4LMEG4B1FVYWEjTpk0Tc4hff/314vLl7dq1o+Li4pDrcNWWrxrpXbzTM9rl1hEF2hQDAEB+CLgOMbnBRrgFsA9fpry8vJxycnJ89/H84Xx7x44dIdc7d+4cXXXVVWJ+8TvvvJM+/fTTkI+tq6ujmpqagMWSbQ9SGEDXBACwAgKuBVC9DQ0nEAB7VVdXU0NDQ7MKLN/mC+gE06dPH1Hdfeedd+iNN94Qlz8fOnQoffll8OO1oKCA0tLSfAuHYtXbBRm2AaCp3x8vcnsTtIGAa3j1FuEWwDx8wRy+FDpfPfKWW26ht99+my6//HL65S9/GfTxfMn0M2fO+JbKykrj2koAUAumCYsxoKlevbULwi2AM9LT0ykhIYGqqqoC7ufb3Lc2HG3atKEbb7yRDh48GPTfk5OTxaILtE8A+kMF12YmVm9x8gBwTlJSEg0cOJBKS0t993GXA77NldpwcBeHvXv3Urdu3Uh3aJ8AzIAKrqEQbgH0wVOE5ebm0qBBg2jw4MG0fPlyqq2tFbMqMO6O0KNHD9GXlj3xxBP0gx/8gK655ho6ffo0LVu2TEwTNnXqVMe2uemvXzIXAwBAPQi4NnZPMK3BRrgFcMf48ePp1KlTNH/+fDGwjPvWbt682TfwrKKiQsys4PW3v/1NTCvGj7300ktFBXj79u1iijGdv5SjjQIwBwKugeyo3uLEAeCuGTNmiCWYrVu3Btx+7rnnxCILJ4oBaKMAzII+uDaRtXqLcAsAbkJ7AQBOQMANArMnhA8nKwCIlNPtJ9opAPMg4NrApOotAICO7SUAqA0BNwqo3n4LVREAkP3LOdopADNhkJkh1Qirq7c4aaj/hQrvIThldLfpRC4UBvAZBzAXAq4BDSLCrZzc/iUgkufHew52fNZkLQgAuOH3x4vc3gStRNVFgScNf+WVV+jQoUNkmpZCgQmNNYJO5J+XUIsu+6HavoC7nPq8oK0CkEdRURFlZmZSSkoKZWdnU1lZWYuP54vV9OnTh9q2bUsZGRk0a9YsOn/+vP0VXL40JF8RZ8qUKeLqOLfccgsNHz5c/Pfaa6+N5k+CAtVbnDBaZnLQC7Xv+MxA089DpMdJNG0YPncA8igpKRFXW1y5cqUItxxeR40aRfv376fOnTs3e/zatWtpzpw5VFxcTEOHDqUDBw7Q/fffT3FxcVRYWGhvBXf16tXiCSsrK+npp5+m9u3b07PPPkt9+/alK64w8yQvY/UW4dY+qGKGB69Tyz766CP6yU9+QkOGDKFjx46J+15//XXatm0b6U7GNhMArMehlK+cyJcO56slctBt166dCLDB8FUVhw0bRvfee6+o+t522200YcKEVqu+ls6iwJd4vOyyy8R/O3bsSImJiXT55ZeTqloLcaaenBFumwc1iB5ey2/95je/EVUM/gnuk08+obq6b7+QnjlzhpYuXer25ikL7RWA/WpqagIWb/vVVH19PZWXl1NOTo7vPr5sON/esWNH0HW4asvreAPt4cOHadOmTXT77bfb30XhscceE5d+5Eb5uuuuE10TuJz893//9yLsgj7VW1NPFiYHL7dfa1M+c0uWLBGVDB7TsG7dOt/9XLngfwMAsFNSZRIlpCRFtE7D+UbxX+4X62/BggW0cOHCZo+vrq6mhoYG6tKlS8D9fHvfvn1Bn4Mrt7zeTTfdRB6Phy5evEgPPfSQyJ62B9ynnnpKVGp5h+666y7q3bs3mUy2n9oQbqODUCsHUwIv9z/jokBTaWlpdPr0adKR3f1vdf2sAMimsrKSUlNTfbeTk63LQVxA5V+xVqxYIfrsHjx4kGbOnEmLFy+mefPm2RtwuXL7pz/9SWwE973lQWfegWa8qBh40T3BTHhf5adr4O3atatouLmPmT/uf9uzZ0/SmR1FAV0+FwAqSE1NDQi4oaSnp1NCQgJVVVUF3M+3uQ0MhkPsxIkTaerUqeL2DTfcQLW1tfTggw/S3LlzRRcH2/rgZmVl0b/8y7/Q22+/TadOnRJ9IzjkTp8+XXRZMAmqt+pB/0+16fL+8aALrkrs3LlTjA7+6quv6M0336RHH32UHn74Ybc3DwAgZpwNBw4cSKWlpb77GhsbxW0eXBvM119/3SzEckhm3GXB1gouPwFXcbmCywtXHLiTcf/+/UUlVzeqnEgRbvV4HyG691S1zy2PW+CG/tZbbxUNOndX4J/5OOD+9Kc/JdNF0p6p9t4DmCQvL49yc3Np0KBBNHjwYDFNGFdkeVYFxuMQeMpZnn6WjR07Vsy8cOONN/q6KHBVl+/3Bl3bAm6nTp3o3LlzopLLgZYrETfffLOYScEkslVvraDbiQKh1hyqdWXgqi3/3Pbzn/9cNODcpvIUOjztoo6imQMXANQ3fvx48Wv//Pnz6cSJEzRgwADavHmzb+BZRUVFQMX28ccfF+0j/5enT+QxXxxun3zyyYieN6qA+8Ybb4hAG07/CxXIfiJ0qnqrw+vAcBIF/8+B7J9r/gmPg62ubh3xbVXGrsKA7O8vABDNmDFDLMFwTwB/POUsT2LASyyiCrhjxowhU5gSlnQ4SZjyXoE53Rh0guMTAJwUVcAFubonxFq9Vf2kjxMnhAthVw3html4D0EXvz9e5PYmaMf4gKt6A2nl5XhVg2ALsUDYVRveMwBoifEBN5oAJVP11sSTBIItmNpfV3U6tZ0AIDejA67qJzPTuiYg2ILdUNVVo13DewMArTE64KrMpHCLYAtuQFUXdOdm24rjCuxmbMBNOvAVUXySkT+xqdKwINiCDFSp6hYVFdGyZcvEPJM8R/kLL7wgJlVvzbp162jChAl055130oYNG0h2Mr8HblGxrbRqm/F5gFCMDbgqNxixVG9VaAxkfu3BbLJWdUtKSsTVglauXCmu/MNXCho1ahTt37+fOnfuHHK9L774Qlw5jec1B3mhTbTutZHt2AX7IOCCNNCIgypkC7p8WUu+oqT30pccdDdu3EjFxcXiksDBNDQ00H333UeLFi2ijz76iE6fPm3rNobz61drX95leb3tgjZQntdY98+aCRBwFeueoGP1Fo06qEqG7gv19fVUXl5O+fn5vvv4spc5OTm0Y8eOkOs98cQToro7ZcoUEXBbUldXJxavmpoai7beXGj31H9/ZD2nwrcQcBVqeBBuAeTlVlW3urpaVGO913X34tv79u0Lus62bdvo5Zdfpj179oT1HAUFBaLS6yZZ27Bwoa3TT2vvqeqfWdUh4CpUvY2WjAcZGnvQFX+2L148T3SSpHT27FmaOHEirVq1itLT08Nah6vD3MfXv4KbkZHR6nq3jigw9qI1aOMAAdhdCLiKiLbxl/EAQsMPYB0OqQkJCVRVVRVwP9/u2rVrs8cfOnRIDC4bO3as777Gxkbx38TERDEwrVevXgHrJCcni8UtMrZjwaBtg0jg82IvBNwm8IGzD15bAOslJSXRwIEDqbS0lMaNG+cLrHx7xowZzR7ft29f2rt3b8B9jz/+uKjsPv/882FVZuE7aNcgVqVbvus/D9ZBwFWge4Lq1VucAADsxd0HcnNzadCgQWLuW54mrLa21jerwqRJk6hHjx6iL21KSgr169cvYP2OHTuK/za9X4Y2VJZ2zB/aNAD5IeBKDuEWAFozfvx4OnXqFM2fP19c6GHAgAG0efNm38CziooKMbOCrFTof4u2DEAtCLgakiHc4mQA4CzujhCsSwLbunVri+uuWbOGZIS2DACihYDbSkPmZvcEFaoaup8Q3Hj/VX3fAXSiUzsGYCIEXM24XfFQ8aTgdh/rSLYH4Rd0E+oz7UZbpmL7BQDBIeBKKpog42a4VeXEIFuYtWr7EXwBzGjDACA8Uow6KCoqoszMTDG6Nzs7m8rKysJab926dRQXF+ebGkeX7gkIt9bg96/poitT9hPUI/vnkdsvWdswAFC4gltSUiKmuFm5cqUItzy9zahRo8Rk43yd9FB4ovJHH32Ubr75Zke3F74j40lB9pOpW68DKrygEie+sMvYfgGARhXcwsJCmjZtmpiv8frrrxdBt127dlRcXBxyHb7u+n333Seujd6zZ0/SiSrVW1lODqhchgevE8jIjS9eqNgCmMHVgFtfX0/l5eWUk5Pz3QbFx4vbO3bsCLneE088Iaq7U6ZMafU56urqxHXT/ReZuydEysRwi6AWO7yGYGKb5nbbBQCGdFGorq4W1VjvZORefHvfvn1B19m2bRu9/PLLtGfPnrCeg6/cw5VeHasZJoVbBDFnXlt0ZQCr3DqigGSBYAtgHte7KESCr5U+ceJEWrVqFaWnp4e1Tn5+Pp05c8a3VFZWkoxkDxZu/ayHKqOzUNkFN9nxpR3hFsBMrlZwOaQmJCRQVVVVwP18u2vXrs0ef+jQITG4bOzYsb77GhsbxX8TExPFwLRevXoFrJOcnCyWSKhwcneyeuv0CUKF198EqOyCnZz4TCHcguxKt+S7vQnacjXgJiUl0cCBA6m0tNQ31RcHVr4d7JKTffv2pb179wbc9/jjj4vK7vPPP08ZGRlKNoAyd01w8vVBsJUXwi6odEzL0K4DgOFdFHiKMO5y8Oqrr9Jnn31GDz/8MNXW1opZFdikSZNENwPG8+T269cvYOnYsSN16NBB/H8OzLrTMdziJ3G14P0Cmds2hFsA+UR6vYPTp0/T9OnTqVu3buJX+N69e9OmTZvUmgd3/PjxdOrUKZo/fz6dOHGCBgwYQJs3b/YNPKuoqBAzKzjF6RO3jNUwJ04QCEjq876HMn6GQW52fWYQbgHkUxLh9Q54hq2RI0eKf1u/fj316NGDjh49KgqaSgVcxt0RgnVJYFu3bm1x3TVr1hjTGOow+TmCrX7QfQFkoFp7DqCymiZTrrY03sn/egeMg+7GjRvF9Q7mzJnT7PF8/1//+lfavn07tWnTRtzH1V/luiiYLJIwoHq4xc/aZsD7DNGItX1DuAWIXIdKD3U4GuFS6RHr8pintLQ038JTslp1vYN3332XhgwZIroo8K/53AV16dKlYlpZ5Sq4bqjv3d3VnZet0mXXCQJhx0zovgBOQbgFcF5lZSWlpqb6boeq3kZzvYPDhw/TH//4R3HFWu53e/DgQXrkkUfowoULtGDBgrC30diAq1IYs7t6i3ALdkHQhaaafhaibd8QbAHck5qaGhBwrcSzaXH/25deeklMJcuzbR07doyWLVuGgCs7mbom2HGSQLCFphB0zYN2AADSI7zeAeOZE7jvLa/ndd1114mJCLjLQ7gzZqEPrsGsDrfofwmtwWcErIDqLYAakvyud+Dlvd4B97MNZtiwYaJbgvdCXuzAgQMi+EYyHSwC7v/n1ElXhuqt1ZfdRWiBSOEzA9G2cQi3AGrJi+B6B4z/nWdRmDlzpgi2POMCDzLjQWeRQBcFSdkZbq2EkAKxQNcFPd06ooDIhrYB4RZAPeMjvN4Bz9Dw3nvv0axZs6h///5iHlwOu7Nnz47oeRFwHRTuSVyFcItgC1ZC0DVDLO+vTOHWyvYPn3kwwYwIr3fA3Rc+/vjjmJ4TAdehsOZ2I4ZwCyrgz5bbxwo4I5Iv8k6HWyfbuHCeC8eEnkq3fPezPFgPfXAlY0f11qqTA/pNghPwObP/Wu9vv/02DRo0SFz68pJLLhE/Gb7++uskI6cuHe6/yKbp9sm8rQCyQAVX864JVoZbACeh24J913rv1KkTzZ07l/r27StGJf/ud78TAz74sbye7nRqz4LtC44ZAARcrRo6O8Ktaq/P2avibH8OvlwhOAfdFloX6bXehw8fHnCbB3DwCOdt27bZFnD938Nwv8yja1V0EHoBEHC1rd7qHm6dCLLRPDfCrz1QzQ3Ne613/2l2WrvWuz+PxyMui8nV3l/84hdBH1NXVycWr5qaGtvbDvz6ZK2mrwOOJdAdAq4EZAu3sp0Q3AyzsW4rAq+1UM215lrv7MyZM2L6HQ6ufMWgFStW0MiRI4M+tqCggBYtWkRO0f0LugwQeEF3RgdcuxtANxoMHcKtSoG2NQi81kM11xodOnSgPXv20Llz58RVhbgPb8+ePZt1X2BcHeZ/96/g8lyVMpKhDVMRAi/oxuiAq1v1VtVwq1OgbQ0Cr3VQzY3+Wu/ebgzXXHON+P88iwJfYYgrtcECbnJyslicaPOibccQbK2FwAuqM3aasJqr3a/eyhJu3ZhuhoOedzEZXofYYKqk6K71Hgyv49/P1krhhqNo2jF8BpyBqclANajgaiDWcOsUhLjwXx9UdiNjejWXuw/k5uaKuW0HDx4spglreq137m/LFVrG/+XH9urVS4TaTZs2iXlwX3zxRVIJwpb7r7vJxx3IDQFX8eqtCuEWwTZyCLuRM7lvbqTXeufw+8gjj9CXX35Jbdu2FfPhvvHGG+Lv2Kmldi+StgzBVh4IuyArBFyLOXmAy9xXDaHW+tcSQTc8plZzI7nW+5IlS8Rih9HdphON6UV2QbiVF8IuyAQB1wV2XLFMlpMDgq19EHTDZ2rIlUE0X7zDXQfhVh0Iu+A2BFyDuibYeXJAsHUOui+Ex+QuCzJp7fVHuNUfwm5zpVu+uzAL2AMBV0EyhVsEW3ehqts6VHOdFaqtieXLPcKtPhB2wSkIuIpVb2UJtwi2ckHQbRlCrpzCac8QbvWFsAt2QsB1iBvhFsHWPAi6oaHLgnoQbs2B4xOshoBrAScOSLfDLYKtWhB0Q0M1Vw6ttWkIt2ZCVResYuyVzFSq3roZbnGVLbXhvQsO4ckZ3oASaRuI9wcYrp4GsUDAjVFr3zBVDbcItvrAexkcTppywvsCwSDoQqQQcCXmZrgF/SDoNocTplxtG94PaA2CLoQLAVfi6m0krDjgEYDMgPc4EE6WcsD7AJFA9wVoDQaZScrJ67Ij8JgHg9ACYQS3fZp+0Y/2EuMAoeD4hWAQcCWs3iLcglP4/UfI/Q5mWHCHDFdZxHGgPgRd8IeAawMVwi2CLXihmhsIITd6t44oIGohZARr32QZO9DS+jg21IKpxoAh4EbBrgMG4RbchGrudxByo+N0f0in2rJgz4NjRQ2o6poLAVeS6q0T4Vb2YFt3Zb3lfzO5Isnyv6kzVHO/g5Brf3sYTVsmSzvmvx04XuQnU9At3ZLv9iYYAbMoRMiOg8OkcMshNtTi5POBGp8Xt2GEtn3tXKSvrcyzvHi3Tdbtg+9g5gV3FBUVUWZmJqWkpFB2djaVlZWFtd66desoLi6Oxo0bF/FzIuBayO5pwaKtdrjZ6MoaLBF6W4eT9bdwMnSfSp9FhF01IOg6p6SkhPLy8mjBggW0e/duysrKolGjRtHJkydbXO+LL76gRx99lG6++eaonhcBV5HqrSo/5akcGlXdbjvhJP0tnAjdqd6qHhRV334TIOjar7CwkKZNm0aTJ0+m66+/nlauXEnt2rWj4uLikOs0NDTQfffdR4sWLaKePXtG9bwIuC5Wb3UJtzoGQx33KVo4SX8LJ0Fn6fSZQ1VXfgi6kampqQlY6uqCFwDr6+upvLyccnJyfPfFx8eL2zt27Aj595944gnq3LkzTZkyhaKFQWYWVG9lCrdONaAmBT//fTV50BpmWcDAs3DxaxSqXQynTdM5CGIgp9xkGoxmt9QjdZSYGNmxdvHit69LRkZGwP3c/WDhwoXNHl9dXS2qsV26dAm4n2/v27cv6HNs27aNXn75ZdqzZ09E29YUAq4LVA23JoXaUEwPuwi5CLl2DqLVOdg2haArN5OCbjQqKyspNTXVdzs52aL5rM+epYkTJ9KqVasoPT09pr+FgOtC9ValcItQ2/prY1rQRchFyLWjXTMp3PpD0JUbgm5wHG79A24oHFITEhKoqqoq4H6+3bVr12aPP3TokBhcNnbsWN99jY2N4r+JiYm0f/9+6tWrF4UDfXAlrGbIEG7R/zR8Jr5WpoYRf+ivZx18ntDXXXbooxudpKQkGjhwIJWWlgYEVr49ZMiQZo/v27cv7d27V3RP8C533HEHjRgxQvz/pl0jWoIKroPVW6vDrV3BFqJjWkUXlVxUciMRqm1DqAuEiq7cUNGNHE8RlpubS4MGDaLBgwfT8uXLqba2VsyqwCZNmkQ9evSggoICMU9uv379Atbv2LGj+G/T+1uDCm6UdAu3JlYh7WLSa4mqkzyV3EgmUuf+bTy35KWXXioWHtEc7sTrdvW/hUA4tuSGim74xo8fT8888wzNnz+fBgwYICqxmzdv9g08q6iooOPHj5PVUMFtgZPf0NwKt6YEMTeYVNFFNVeOidR5fkkOt1wh4YnUub8aT7XT1NatW2nChAk0dOhQEYh/8Ytf0G233UaffvqpqKQ4VQBAgGsdji25oaIbnhkzZoglGG6PWrJmzRqKBiq4ElRv3Qi3JlUZ3WbK62xyWHG7khPpROpvvvkmPfLII6Kawn3eVq9e7esXFwzPcdl03sumBk0tDHmiD/b6mPx5iRSqufJzuw2A5hBwQ7Dq25hs4RbB1h2mvO4mn4TdOsFFO5G6v6+//pouXLhAnTp1Cvrv3DcuLS3Nt4Qz0KOlts/kz0ksEHTlhm4LckHAtbF6K1O4NSVgyc6E98DkE7AbJ7eWJlI/ceJEWH9j9uzZ1L1794CQ7C8/P5/OnDnjW3gOzHDhhG89k48x1YNu6ZZ8x7fHVOiDG4QTfWkiuRZ7rEwIVSoxoW+uyf0G+di+ZP95UsVTTz1F69atE/3guD9uMDyJu2UTuSOcWQKzLcgP/XPdJUUFV8bRv3ZWb50Kt6jayg3vjb5qrnauahnpROr+eGQzB9z333+f+vfvb/OWItzaAa+p/PArhqEB1zv6l69jvHv3bsrKyhKjf0+ePBn08d7Rv1u2bBH9y7gvGI/+PXbsmCXbY8U3LSumxoml0UKwVYfO7xVOvHJOpO719NNP0+LFi8V0PTw/pWXbs+9LXxuIE7sz0DdXfuifa2AXBf/Rv4xH/27cuFGM/p0zZ07Q0b/+ePTvb37zG9GY82TBbldvreh3G2u4BfXw+6ZjlwWTuyrIOpE642nBeE7KtWvXil/PvH1127dvLxY7RNOuRdue6XgshQPHG4AkAdc7+pcHMNg1+pent+HFK9j0Nk6xM9wi2KoPIRdimUj91KlTIrRyWOXpv5pOpM5tq9eLL74o2t+777474O/wL2kLFy4kN1nRlvn/DR2PqZagby6ABAG3pdG/+/bts2T0L1csFi1aFFP3BCuqtwi3EA6EXHBiIvUvvvjClm3wb0P927xw2ja72jFTwy6OOTCd631wrRj9+9vf/jbk6N9wp7dRMdzq3H/TZLq+r+gjaIZoxiA49XnX9dgKBcccmMzVCq4Vo3//8Ic/tDj618rpbWQLt7LIvOKUZX/riy8vt+xvqU7Hai6qSuZpqX1zqx0zYao+LxxzYKpEWUb/jhs3LmD0b6if2ryjf5988kl67733LBn9G2v11slw62awtTLIRvIcJodeHUMu6C+ctk+GL+mmBF30ywUTuT6Lggqjf+2aEkz2cOtEoI1mO0wLvLqFXFSU9OU/RVhLbZwM4VbnYywUHHtgkkTTR//GUr2NpWtCJOHWqZOBLIG2NSYGXt1OwDjRmku2cGtiNRfHHpjA9YDr5uhfuy6fp1K4VSXUhrMPugddhFyQ1aCphWH9oiVruNX5OAsGxx6YQIqA64bUI3VEiSmWV29VCLc6hNrW9kvXsKvbyRcnWv1420CVR/CbUM3FsQe6MzbgxkLFcKtrqDWxqqtbyAUzqFC9Ne1YQ8h11q7VeW5vglGUngfXDrHOnBBruLV6nkYOeqaFWxP2X8WwEIrKlT7Q//Oq8raHA8cf6AoB14HqbSTh1iq6Brto6fh66HTixUlWff5jGnR7P3U61oLh90u39wwAATeC6q0K4VbHIGcl3V4b3U+8oJZgbaEun1Fd9qMlCLmgE/TBjZEM4Va30GY3nfvnqgz9AfWjWyjUvU8uw3EIukAFN4bqbazh1or+tgi30dPltdMpRKCCpD7d30Orx0kAgD0QcGPomhBruI0FuiJYQ5fXESdckI3un0md90/3LylgBgTcKIWq3todbnUJZLLR4TXV5YSLkyuA+3AcguqM74NrZdcEO8OtDgFMdvwaq94vV5c+gugHqB5uFyMNRa21a7Ifj7ocb6HgOASVGR1wVQi3CLbOwgA0gNi11NZF0qapcHVChFwAORnbRSHpwFeOPh/CrVpUfu3RVQF0PK5k7p6lyzEXCo5FUJGxAbc1VlZvI238ZG7ITaLye6D7CRfU++xZdTzJelzqfswh5EIsioqKKDMzk1JSUig7O5vKyspCPnbVqlV0880306WXXiqWnJycFh8fCgKujeE2mulkZG28TYX3w104qerxXll9HMlaBEDIBWiupKSE8vLyaMGCBbR7927KysqiUaNG0cmTJ4M8mmjr1q00YcIE2rJlC+3YsYMyMjLotttuo2PHjlEkEHDDEG241aHBBnVDru4nW5DDoKmFrh0/qh6bACYpLCykadOm0eTJk+n666+nlStXUrt27ai4uDjo499880165JFHaMCAAdS3b19avXo1NTY2UmlpaUTPi4Ab5Zy3VoVbBFs1qPoe6RByUTVS9/PmxHEj27GpwzHXEhyPwGpqagKWuro6Cqa+vp7Ky8tFNwOv+Ph4cZurs+H4+uuv6cKFC9SpUyeKhNGzKETbNcHKcAvq0GEaMQCrpR2qC9pOOtm+yXZsYmYFaGrX6jyScbB9Ynxkn9P4xm8zDncb8MfdDxYuXNjs8dXV1dTQ0EBdunQJuJ9v79u3L6znnD17NnXv3j0gJIcDAbcFuoXbkV3D+zBF4oMTfckksp1ITTnZ4oQa3iCOZcuW0YkTJ0QftxdeeIEGDx4c9LGffvopzZ8/X1RWjh49Ss899xz97Gc/I5XJdmzqcNy1BMek2SorKyk1NdV3Ozk5eFfOWD311FO0bt060S+XB6hFAgE3RPXWrnDrVLC1I8yG+zy6h17ZTqQA3kEc3LeNRygvX75cDOLYv38/de7cOehPfj179qQf/ehHNGvWrJifn9tG/zbQrV+ncGwCOCM1NTUg4IaSnp5OCQkJVFVVFXA/3+7atWuL6z7zzDMi4P7hD3+g/v37R7yN6IMbQb9bmcMtB03/xU2ybAfo1S8Qff+sG8Txd3/3d6La++Mf/9i2yotbZOr6pcNx1xIck9CapKQkGjhwYMAAMe+AsSFDhoRc7+mnn6bFixfT5s2badCgQRQNVHCDiKR6G04DZleDq0KA9N9GnSq7KlaKdP/J1FTeQRz5+flRD+JoDQ8g8R9EwoNKVAiYMtD9uENXBWgN/7qUm5srgip3m+JfmGpra8UXcjZp0iTq0aMHFRQUiNu/+MUvRBeqtWvXirlzudsVa9++vVjCZXzAjaVrghvhVoVQ29q26xJ0VQy5qsPJlGwZxNEaPvEsWrQoomkU3YRjE0Ae48ePp1OnTonQymGVp//iyqy3zaqoqBBfyr1efPFF8cX97rvvDmsgWyhGB1xVwq3KoVb3oKvaiVT3ahLYg6vDXIXxr+D6j6L2tocyVW9lOjZ1P+7wxRNaM2PGDLEEwwPI/H3xxRdkBWMDbn3v7q3ufLTh1opGXrdQG2ofdQi5AG6KZRBHuLifrm59dcFaCLkgGwwyC1G9dSvcmjY4S4f9lalqZcLAFwxssWYQh9Xvh4zHgUzbpPpxB6AaBNwI+pDZGW51CHqxUH3fZTqRgnm4+8CqVavo1Vdfpc8++4wefvjhZoM4/Aehcf+2PXv2iIX/P1/jnf//wYMHSTc4Np2DL58gE2O7KETa79aucKt6sLOSTn1zwV74OTS2QRxfffUV3XjjjQHzTfJyyy23NOsPFysT58o2tS8ugEyMDrhWhFsEW+up2jdXpkEtrcGJ1uxBHDz1jsfjseyzFKwdbKmNc/rLrErHpurw5RNkYWwXhZqr3Qm3pndFCBdeIwD9j18Tj3P0xQVwhrEB1+lwi2AbORVfL5X6+6l+okV/P/f9n2v/Nebj1qnjXKVjU3U4NkEGCLg2h1sE29jgtQOQe7pFK45X045z1b9cAqjA+IAbSbjlYBtuuEWwtY5qryMqRWBSVy+VLmiDY9M5qOJ+Z9fq7y6SAs4xOuBGGm51DWQqwGtqD9UrSTiJuutsRpylxyiOcwCwirEB179htircomprL5VeW1SKAEDnL5cAsjM24FoZbhFsnYPX2Xo40UK06jPqLT8u7T7G8eXTOfiFBdyEgGtBuAUIBidSAAAAdxgfcIOF23AGk6Fq6x687uAPVSL9jkdTjnETfj3B8QluiTf5p7VQ4bYlCLZywHsA4L4ru1eTivDrCoD+jL5UbzTh1kl3p+5u9THra77vyLaA3pcIxaV7AQBAJwi4YYRbO4NtOCE2mvVNCL78vjh1LXsAcJYpxze+XALYw/iA62TVNtYwG+1zmRB2wWzcz6/DUY/bm2GcEV0OEFEbUpEqv67oAMcnuMHogGt31dbJQGtq2FWhyoMTKQAAgLMSzR4ckWx5sJUl1La2fboFXYgdfioF2ajwBRYA5GRswA0m2nAre6gNtc26hFycBAHc1bQN1KVtAQB1GTtNWKxTf3GD7l1UpfK2qwbTEoGugrUjaFsig/lwAaxnfAU3kmCrY6ONLgsAEK1xHf5vyDqJTr8SAYB64k0e/RtuuFW9UhsO1fcPF34AVIj0ZOexjV9WAPRlfAVX18AHAOA2VHHBdLtW57m9CcYytoJrcrU2FFP3G8zqCwgAAPpDwP3/TA62/vAa2Ac/hwIAADgj0eTBEe07IN/rBNOFAQAAAEPCg2ZQxQUAAACrFBUVUWZmJqWkpFB2djaVlZW1+Pi33nqL+vbtKx5/ww030KZNmyJ+TgRcAAAAALBFSUkJ5eXl0YIFC2j37t2UlZVFo0aNopMnTwZ9/Pbt22nChAk0ZcoU+uSTT2jcuHFi+ctf/hLR8yLgAgAAgO0wlZ+ZCgsLadq0aTR58mS6/vrraeXKldSuXTsqLi4O+vjnn3+e/uEf/oF+/vOf03XXXUeLFy+m73//+/Qf//EfET2vcX1wPR6P+O+5c41ub4rU/iFuF204m0WquVhbRzJr/OY8ya7hvLrHRkP9t8e32xrqzwe0Nya3o+fPXZD6uJblmFT5uFPxGHVKTU2N7X872nbmoqeeqDGKdYLsV3Jysliaqq+vp/LycsrPz/fdFx8fTzk5ObRjx46gz8H3c8XXH1d8N2zYENG2Ghdwz549K/47PBsj2lv3PqlHxW0GndubtLQ0MrsdteKYxHENakp7fa507UxSUhJ17dqVtp54Larna9++PWVkZATcx90PFi5c2Oyx1dXV1NDQQF26dAm4n2/v2xf8Ii4nTpwI+ni+PxLGBdzu3btTZWUldejQgeLi9Py5hL9Z8YeP9zM1NZV0g/1Tn+77yBUVPulwe6Mj3dpR3T+P2D899y/adiYlJYWOHDkiqqvR4OdtetwHq966zbiAy6XxK664gkzAB4KOB7sX9k99Ou+jjpVb3dtRnT+PDPun3/5F286kpKSIxW7p6emUkJBAVVVVAffzba4iB8P3R/L4UDDIDAAAAAAsx90hBg4cSKWlpb77Ghsbxe0hQ4YEXYfv9388++CDD0I+PhTjKrgAAAAA4AweMJabm0uDBg2iwYMH0/Lly6m2tlbMqsAmTZpEPXr0oIKCAnF75syZdMstt9Czzz5LY8aMoXXr1tGuXbvopZdeiuh5EXA1xH1huMO3jH1irID9U58J+wjq0P3ziP1Tm+r7N378eDp16hTNnz9fDBQbMGAAbd682TeQrKKiQnR78ho6dCitXbuWHn/8cXrsscfo2muvFTMo9OvXL6LnjfPoOo8NAAAAABgJfXABAAAAQCsIuAAAAACgFQRcAAAAANAKAi4AAAAAaAUBV1FFRUWUmZkpJmrOzs6msrKykI9dtWoV3XzzzXTppZeKha8B3dLjVds/fzydCF9hZdy4caTT/p0+fZqmT59O3bp1EyNpe/fuTZs2bSKd9pGnjunTpw+1bdtWXLVn1qxZdP78ece2F/TElw/lNsF/6du3b4vrvPXWW+Ix/Nm94YYbpD7W+Bhrun+8cHsRzJo1a5o91okJ/8P14Ycf0tixY8XVuXjbePS8Px4Xz6PxuS3ktoLPZ59//rlt5xQn9+/ChQs0e/Zs8Zm75JJLxGN4Cq2vvvrK8s+4CRBwFVRSUiLmleNpQ3bv3k1ZWVk0atQoOnnyZNDHb926lSZMmEBbtmyhHTt2iPBw22230bFjx0iH/fP64osv6NFHHxVhXmaR7h9fTnHkyJFi/9avX0/79+8XX1p43kBd9pGnhJkzZ454/GeffUYvv/yy+Bs8RQxArL73ve/R8ePHfcu2bdtCPnb79u2ivZwyZQp98skn4ssyL3/5y19IRv/1X/8VsG88IT770Y9+FHIdvhqW/zpHjx4lWfD8qNxecCAN5umnn6Z///d/p5UrV9LOnTtFEOS2paUvw9GeU5zev6+//lps37x588R/3377bdHe33HHHZZ+xo3B04SBWgYPHuyZPn2673ZDQ4One/funoKCgrDWv3jxoqdDhw6eV1991aPL/vE+DR061LN69WpPbm6u58477/TIKtL9e/HFFz09e/b01NfXe1QR6T7yY3/4wx8G3JeXl+cZNmyY7dsKeluwYIEnKysr7Mffc889njFjxgTcl52d7fnnf/5njwpmzpzp6dWrl6exsTHov7/yyiuetLQ0jwo4ovz2t7/13eZ96tq1q2fZsmW++06fPu1JTk72/OpXv7LtnOnU/gVTVlYmHnf06FHLPuOmQAVXMVzNKy8vFz/LePEEyXybq7Ph4G+J/FNIp06dSJf9e+KJJ6hz586i6iKzaPbv3XffFZco5J8ceWJsnux66dKl1NDQQLrsI0/szet4fzY8fPiw+Fn49ttvd2y7QV/8Ezb/3NuzZ0+67777xMTyofBn1P+zy7jaF2776vax98Ybb9ADDzwgfqYO5dy5c3TVVVeJX/PuvPNO+vTTT0kFR44cERcK8H9/0tLSRJeDUO+PFedMN505c0a8lx07drTsM24KBFzFVFdXi2DjvQKIF9/mAz8c3MeHD4Smjbiq+8c/xfBP2vyzveyi2T8Oe9w1gdfj0Mc/X/ElDJcsWUK67OO9994rvqTcdNNN1KZNG+rVqxcNHz4cXRQgZhx+uN8pXznpxRdfFCGJuzGdPXs26OP5MxpL++om7s/J/fXvv//+kI/hfu7FxcX0zjvviDDc2NgovmB++eWXJDvvexDJ+2PFOdMt3O2Cz9fcZYa7lVj1GTcFLtVrmKeeekoMxOJ+uTINLIgWH8ATJ04U4TY9PZ10xCcgrk7zdbgTEhJo4MCBov/0smXLRJ8yHfDnkavSK1asEI31wYMHxfXIFy9eLAI9QLRGjx7t+//9+/cXny+uXv7617+W/hefSPEXfd5fLmCEwr8G8eLF4fa6666jX/7yl+J4Aznwr6z33HOPGFTHobUlJn3GI4GAqxgOcRxyqqqqAu7n2127dm1x3WeeeUYE3D/84Q/iINBh/w4dOiQGX/GoVP9AyBITE0UHfa4Gqvz+8Whhrmryel58QuLqA//8lpSURDKJZh85xPIXlalTp4rbPIqYB2M8+OCDNHfu3IDrlAPEgn/q5VlI+EtUMPwZjaZ9dRsPFOO2nQcmRYLblhtvvDHk6yET73vA7we3i158e8CAAZafM90Ot/ye/vGPf2yxehvNZ9wUOGsohsMMV/BKS0sDAh3f9v9WHmzkKX87558wBg0aRLrsH0+FsnfvXtqzZ49v4RGnI0aMEP+f+5ip/v4NGzZMNFTe4M4OHDggGnjZwm20+8j9wpuGWG+g/3YsBoA1uP8pfzH2D0j++DPq/9llPDNBS+2rDF555RXxS8+YMWMiWo9/vuc2NNTrIZOrr75ahFL/96empkbMphDq/Yn2nOl2uOU+tfyF5bLLLrP8M24Mt0e5QeTWrVsnRo2uWbPG8z//8z+eBx980NOxY0fPiRMnxL9PnDjRM2fOHN/jn3rqKU9SUpJn/fr1nuPHj/uWs2fPenTYv6Zkn0Uh0v2rqKgQs17MmDHDs3//fs/vfvc7T+fOnT1Llizx6LKPPAqY95FHQh8+fNjz/vvvi5HgPKIdIBb/+q//6tm6davnyJEjnj//+c+enJwcT3p6uufkyZNBP4v8mMTERM8zzzzj+eyzz8Rns02bNp69e/d6ZMWzAlx55ZWe2bNnN/u3pvu3aNEiz3vvvec5dOiQp7y83PPjH//Yk5KS4vn00089MuDz0ieffCIWjiiFhYXi/3tnEeDzGbcl77zzjue///u/RVt/9dVXe7755hvf3+AZWV544YWw2yNZ9o9nyrnjjjs8V1xxhWfPnj0B5+u6urqQ+9faZ9xUCLiK4g83N2gcXHkKlI8//tj3b7fccosIeV5XXXWVOJCaLtxw67B/qgXcaPZv+/btYqoibqR5yrAnn3xSTI2myz5euHDBs3DhQhFq+WSbkZHheeSRRzx/+9vfXNp60MX48eM93bp1E5/DHj16iNsHDx5s8Xj79a9/7endu7dY53vf+55n48aNHplxYOU2nb8AN9V0/372s5/5jssuXbp4br/9ds/u3bs9stiyZUvQ85V3H3iqsHnz5olt5/bw1ltvbbbffM5ren5rqT2SZf84oAb7N154vVD719pn3FRx/D9uV5EBAAAAAKyCPrgAAAAAoBUEXAAAAADQCgIuAAAAAGgFARcAAAAAtIKACwAAAABaQcAFAAAAAK0g4AIAAACAVhBwAQAAAEArCLgAAAAAoBUEXAAAAADQCgIuAAAAAGgFARcghMzMTFq+fHnAfQMGDKCFCxe6tk0AALEYPnw4zZgxQyxpaWmUnp5O8+bNI4/H4/amAVgKARcAAMAgr776KiUmJlJZWRk9//zzVFhYSKtXr3Z7swAslWjtnwMAAACZZWRk0HPPPUdxcXHUp08f2rt3r7g9bdo0tzcNwDKo4AIAABjkBz/4gQi3XkOGDKHPP/+cGhoaXN0uACsh4AKEEB8f36xf2oULF1zbHgAAAAgPAi5ACJdffjkdP37cd7umpoaOHDni6jYBAMRq586dAbc//vhjuvbaaykhIcG1bQKwGgIuQAg//OEP6fXXX6ePPvpI9FHLzc3FCQAAlFdRUUF5eXm0f/9++tWvfkUvvPACzZw50+3NArAUBpkBhJCfny8qtv/4j/8optNZvHgxKrgAoLxJkybRN998Q4MHDxZf2jncPvjgg25vFoCl4jyY/A4AAMCYeXB5Pu+mc3wD6AZdFAAAAABAKwi4AAAAAKAVdFEAAAAAAK2gggsAAAAAWkHABQAAAACtIOACAAAAgFYQcAEAAABAKwi4AAAAAKAVBFwAAAAA0AoCLgAAAABoBQEXAAAAAEgn/w+IpPz8+fStVAAAAABJRU5ErkJggg==",
                        "text/plain": [
                            "<Figure size 800x400 with 3 Axes>"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                }
            ],
            "source": [
                "# Look at the contours of the (2,2,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)]]),\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)]]),\n",
                ")\n",
                "plt.colorbar(cb)\n",
                "plt.xlabel(\"p\")\n",
                "plt.ylabel(\"e\")\n",
                "plt.show()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 9,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAGCCAYAAAALybB5AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAATzZJREFUeJzt3Qt4VOW1//FFEkJAIIgBgsARRQUpCjxQKOCFtigqpXLaKqIFpIpHJS2V6l8QuShKvEbEolTk4qlFqB7RHgWsTYvKIZCayCkqoFyUoFzbQrhJJJn/s17PpJlkJpnLntm37+d5dnGme2f2QOad36xZ77sbBQKBgAAAAAAekWb3CQAAAABWIuACAADAUwi4AAAA8BQCLgAAADyFgAsAAABPIeACAADAUwi4AAAA8BQCLgAAADyFgAsAAABPIeACAADAUwi4AAAADvXuu+/K8OHD5cwzz5RGjRrJa6+9ltTHmzlzpnmcmlu3bt3EbQi4AAAADnXs2DHp2bOnzJs3L2WP+a1vfUv27NlTva1du1bcJsPuEwAAAEB4V111ldkiOXnypEydOlVeeuklOXTokPTo0UMeeeQRGTx4cNyPmZGRIbm5ueJmVHABAABcKi8vT4qKimTZsmXyt7/9Ta699lq58sor5dNPP437Z+qx2hJxzjnnyI033ii7du0St2kUCAQCdp8EAAAA6qf9sCtWrJARI0aY2xo8NYTqnxpIg4YMGSL9+vWT2bNnx/wYq1atkqNHj0rXrl1Ne8L9998vX3zxhXz44YfSokULcQtaFAAAAFxo06ZNUllZKeeff36dtoUzzjjD/PeWLVvkggsuqPfn3HPPPfLwww+b/67ZDnHRRRdJ//795ayzzpLf//73cvPNN4tbEHABAABcSCut6enpUlJSYv6sqXnz5uZPrfBu3ry53p8TDMPhtGrVygTobdu2iZsQcAEAAFyod+/epoK7f/9+ueSSS8Luk5mZmdAyXxqit2/fLqNHjxY3IeACAAA4lAbMmtXTnTt3ysaNG6V169amsqqTwMaMGSNPPPGECbwHDhyQwsJC014wbNiwmB/vrrvuMuvualvCl19+KTNmzDDV4VGjRombMMkMAADAodasWSPf/e5369w/duxYWbJkiXz99dfy4IMPyn/+53+ayWA5OTnyne98x0wOu/DCC2N+vOuvv95cXOLvf/+7tGnTRi6++GJ56KGHpEuXLuImBFwAAABYLj8/X1599VUz0a1p06YycOBAs0avrtBQnzlz5sizzz5rVofQwP6Tn/zE/KysrKyoH5t1cAEAAGC5d955RyZMmCDr16+Xt99+21Sbr7jiCnN1tkiWLl0qkydPNq0ROjlu4cKFsnz5crn33ntjemwquAAAAEg67Q9u27atCb6XXnppxAtXaLDVPuKgX/3qV7Jhw4aYLhnsu0lmVVVVpmlaFyvWBZMBwGpaNzhy5IhZeD0tzXtflDGOAu4eZ7766iupqKiI+3Frv+6bNGlitoYcPnzY/KkT5CLRNoYXX3xRiouLzcUqduzYIStXrox9FYeAz5SVlWnFmo2NjS3pm443XsQ4ysbm3nHmxIkTgZw2aXE/XvPmzevcN2PGjAYft7KyMjBs2LDAoEGDGtz3qaeeCjRu3DiQkZFhfv5tt90WiJXvKrjBy8wNbjdOMtIy6933RPf2de5r+vGeqPaLJNzx8fwcKx4vWsk4r1Q/ByCVTlVVyJp9i111WctYBJ9XWVmZtGzZ0u7TATznx93+X9LGmYqKCjl4oEreWp8rpzWPrfJ77GiVDP3O3jqv/Wiqt9qLq5f7bajNQFeN0EsMP/PMM+YqarpE2sSJE2XWrFkybdq0qM/VdwE3WFbXcNtQwM3ICJ2t1/TDL0RqHXOiR4eo/xLDHV/f41mhoefYkFieX7Ik+hwAu3j16/vg89I3OAIukLirOv4i7ve9eMeZ05qnSfMW8bVQxfra177aN954wyw/1rFjx3r31RCr7Qi33HKLua1LnemktFtvvVWmTp0adTuG3dkFNYKk1Uygdtg5AQDgd7UDrVcFAgH5+c9/LitWrDCV2bPPPrvBY44fP14nxAYvQxzLuggE3BSFRbch3AIAYA2/BNpwbQm67Nfrr79uWin27t1r7s/Ozjbr4iq9CluHDh3MOrdKr6JWUFBgrsoWbFHQqq7eHwy60SDgOiAgO7F66xReeR4AAH/xa6itSS/WoAYPHhxy/+LFi+Wmm24y/60Xc6hZsb3vvvtM24X+qVdm06upabjVq6nFgoAbgZ8rmH5+7gAAxINAW1c0LQXaulBTRkaGuciDbokg4HowCCZS9XTqcwIAwEkItM5GwPVge0K8nHQuAAA4CYHWXQi4UQQ9+kABAPAfQq17EXA9Jt4wTvUWAOB3bgq0/7XlUcnO/o3dp+FYBNw4RRsI3dKe4ERUzgEAyeamUIvoEXA9FLKo3gIAUD8CrT8QcH2OcAsA8DpCrf8QcOOoiNKeAACAcxFoQcD1iHjaE5wcst3UGgIAsB+hFjURcAEAgCsRahEJATdJUtme4LXqLQAAkRBqEQ0Crg85PdzSngAACCLQIh4EXJcjDAIAvIZQi0QRcJNQ/XTy6gl2Pz4AAOEQamElAi4chYo0APgHoTY+q3bPlfLycrtPw9EIuC4OWrGeI9VbAIDdCLVIBQJuChEwAQB+RKhFqhFwLb6CWaoqv16s3rqhag4AiA6hFnYi4AIAAEsQauEUBNwUsbOCSvUWAJAshFo4EQHXQk5tTwAAwGoEWzgZAdfj3FC9BQC4A6EWbkHATUFQtDJkpqp6+8+umSG3T99akbTHoiINAM5FqIUbpdl9AkieeIN17XAb6T4AzjFv3jzp3LmzZGVlSf/+/aW4uLje/efMmSNdu3aVpk2bSqdOneTOO++Ur776KmXnC+eH2uAGuBEVXIt4oQppR4j1wt8bYLfly5fLpEmTZP78+SbcangdOnSobN26Vdq2bVtn/6VLl8rkyZNl0aJFMnDgQPnkk0/kpptukkaNGklBQYEtzwHOQKCFVxBwk8yu9oRYHzeacKv7JLNVAUB8NJSOHz9exo0bZ25r0H3zzTdNgNUgW9u6detk0KBBcsMNN5jbWvkdNWqUbNiwIeXnDvsRauFFtCjAtvYDqrdA4ioqKqSkpESGDBlSfV9aWpq5XVRUFPYYrdrqMcE2hh07dsjKlSvl6quvDrv/yZMnzXXva25wP1oQ3GnV7rl2n4IrUMH1oFiqt/TWAu528OBBqayslHbt2oXcr7e3bNkS9hit3OpxF198sQQCATl16pTcdtttcu+994bdPz8/X+6///6knD9Si0ALv6CCa0FgTEUl0mvVTq89H8BN1qxZI7Nnz5ZnnnlGSktL5dVXXzUtDbNmzQq7/5QpU+Tw4cPVW1lZWcrPGfFjwhj8iAqux9agpXoL+EtOTo6kp6fLvn37Qu7X27m5uWGPmTZtmowePVpuueUWc/vCCy+UY8eOya233ipTp041LQ41NWnSxGxwFwIt/IyA69NqJ+EW8IbMzEzp06ePFBYWyogRI8x9VVVV5nZeXl7YY44fP14nxGpIVtqyAPci1ALfIOB66Kv2aKu3dodbJ/2dAV6gS4SNHTtW+vbtK/369TPLhGlFNriqwpgxY6RDhw6ml1YNHz7crLzQu3dvs6zYtm3bTFVX7w8GXbgLwRYIRcBF1FgiDHCmkSNHyoEDB2T69Omyd+9e6dWrl6xevbp64tmuXbtCKrb33XefWfNW//ziiy+kTZs2Jtw+9NBDNj4LxINgC4RHwE1SH61V/bfRVjup3gL+pu0IkVoSdFJZTRkZGTJjxgyzwX0ItXCL/Px8M4lVV3TRqybqEoWPPPKIuYpiNJYtW2bW6L7mmmvktddei+mxWUUBAAAXYCUEuM0777wjEyZMkPXr18vbb78tX3/9tVxxxRWmhaohn332mdx1111yySWXxPXYVHB9VI2kegsA7kOohVutXr065PaSJUvM5cP1QjOXXnppxON0be8bb7zRrL/93nvvyaFDh2J+bAJuEjixPcHucAsAiB6hFk5WXutqhtEuJajraKvWrVvXu98DDzxggvDNN99sAm48CLhICaq3ANAwgi1SdZneZYf6S5NTjWM65uTRr0VkhXTq1Cnkfu3nnzlzZr3H6vKFv/zlL2XQoEHSo0ePiPutXbtWFi5cKBs3bpREEHAdymnVW1ZQAIDkIdjCTcrKyqRly5bVt6Op3mov7ocffmgCbCRHjhwxF6FZsGCBuYhNInwfcDUgRgqTdlyJzIuo3gJAeARbuFHLli1DAm5DdIWXN954Q959913p2LFjxP22b99uJpfpsoU1K7/B1V+2bt0qXbp0ieoxfR1wEwmwTght9N4CgDsRbOEHgUBAfv7zn8uKFSvMcoVnn312vft369ZNNm3aFHKfrtetld2nnnqqTmtEfXwdcJPBiqqvE8KzVe0JTnsuAGAXQi38ZsKECbJ06VJ5/fXXpUWLFuZCNCo7O9usi1v7SotZWVl1+nNbtWpl/qyvbzccAq6HUb0FAPsRbOFXzz77rPlz8ODBIfcvXrxYbrrpprBXWrQKAdelVUk39Ac74e8JAOxCsIXfBQKBBvepfaXF2nTt3HgQcD3ankD1FgDsQbAF7OfbgHuie3vXPvlUVm/j7b+legvAbwi2gHO4NeMlnRtaACKhegsAqUOwBZyHgOsgXql6euV5AEB9CLaAcxFwPVZZtrJ6y9XLAKAugi3cfplePyDgxoirntWP6i0AryLYAu5h/cJjMZo3b5507tzZLO7bv39/KS4urnf/OXPmSNeuXc0CwXpFizvvvFO++uorS8/JjrAaTTBMZfUWAPCvYEu4BdzF1gru8uXLZdKkSTJ//nwTbjW8Dh061FxruG3btnX216thTJ48WRYtWiQDBw6UTz75xCwU3KhRIykoKLDlOXhVPO0JbqzeNvShwY3PCYB1CLaAO9kacDWUjh8/XsaNG2dua9B98803TYDVIFvbunXrZNCgQXLDDTeY21r5HTVqlGzYsEH8juptcirz4Y4h9ALeR7AF3M22FoWKigopKSmRIUOG/Otk0tLM7aKiorDHaNVWjwm2MezYsUNWrlwpV199dcTHOXnypJSXl4dsVoeVRFsarGhPsJuTQ5/+3Vn592f1zwPgHLQjAN5gWwX34MGDUllZKe3atQu5X29v2bIl7DFaudXjLr74YnP5t1OnTsltt90m9957b8THyc/Pl/vvv9/TIc7q6q0XVk9IRQANPoZTfy8ARI9QC3iL7ZPMYqHXK549e7Y888wzUlpaKq+++qppaZg1a1bEY6ZMmSKHDx+u3srKylJ6zn7gpIBnR3WVii7gXlRsAW+yrYKbk5Mj6enpsm/fvpD79XZubm7YY6ZNmyajR4+WW265xdy+8MIL5dixY3LrrbfK1KlTTYtDbU2aNDFbtMEs1UHFC+0JTmH33xMVXcBdCLaAd9lWwc3MzJQ+ffpIYWFh9X1VVVXm9oABA8Iec/z48TohVkOy0pYFP4Yqu9sTnBDmnFZBddK5wD9iWXJx8ODBZvWZ2tuwYcPED6jaAt5n6yoKukTY2LFjpW/fvtKvXz+zTJhWZIOrKowZM0Y6dOhg+mjV8OHDzcoLvXv3NgP4tm3bTFVX7w8G3WRxQpCDe4Kknhu/M3Dqkova3qUTfYP+/ve/S8+ePeXaa68VLyPUAv5ha8AdOXKkHDhwQKZPny579+6VXr16yerVq6snnu3atSukYnvfffeZKoP++cUXX0ibNm1MuH3ooYfEjyHOz9VbJ4fbIFoW4NQlF1u3bh1ye9myZdKsWTPPBlyCLeA/tl+qNy8vz2yRJpXVlJGRITNmzDCbFxB8vBtua6Kai1QsuagTaqNdcrG2hQsXyvXXXy+nnXZaxOUWdQuqvdyikxFu4QWrds+1+xRcx/aAC3cGOTsCm5P/PhpCyIWTllysSXt1P/zwQxNyU7HcYqoQbAF/c9UyYVZq+vEeVwcuu9sTUs3N4dZLzwHeo8FWV6TReRBeWG5x8JWPmA2Av1HBtSmAuLmal+pzJxgC1i65GKSTerX/9oEHHqh3v/qWW3QSgi0A8XsF1+lSObnMybwWbr32fGC/eJZcDHr55ZdNb+1Pf/pTcTOqtgBqI+B6vNrq5vYEr4ZBrz4v2EeXCFuwYIG88MILsnnzZrn99tvrLLlYcxJazfaEESNGyBlnnCFuRbAFEA4tCnHya0hJVdj3+t8vk85g55KLStfIXbt2rfzxj38UNyLYAqgPAdeB/bdOXvs2FbweboMIubBryUXVtWtX264AmSjCLYCGEHARtVSEMb+EWwCxI9gCiBY9uHAMP4ZbPz5nINnhlm9GAFDBTfFA6db2hGS/Yfg56NGqAERxwQYfjxEAYkcFNw5+DmPJwN8ngHC4GhnAZXrjRQXXJexc+zaZ1UXC7Teo4gKh7Qg1Xw+MEwBiRQXXp4OnE1ZP8OPfO4D6MZEMgBUIuCmUSP+tXagqpo4T//0Bt4VbxiwAihYFF/Di2reEOQA1e20jjQmMFQDiQQU3xkqAnwbbZFVC/PR3GCv+buDHqi2/9wCsRsB1CL8M8H55ngAaRr8t4G35+fny7W9/W1q0aCFt27aVESNGmMuEN+Tll1+Wbt26SVZWllx44YWycuXKmB+bgOvwaqhd7QnJqN4SbqPD3xP8Ljj+xPpaoP8WcJZ33nlHJkyYIOvXr5e3335bvv76a7niiivk2LFjEY9Zt26djBo1Sm6++Wb54IMPTCjW7cMPP4zpsX3fg0uYAAAAsN7q1atDbi9ZssRUcktKSuTSSy8Ne8xTTz0lV155pdx9993m9qxZs0w4/vWvfy3z58+P+rGp4DogDPshZPvhOVqJvy8AgFOVl5eHbCdPnozquMOHD5s/W7duHXGfoqIiGTJkSMh9Q4cONffHwvcVXCfzSnsCYQ1A2KuUMTYAtvnz7vMkvVmTmI6pPP5NkO3UqVPI/TNmzJCZM2fWe2xVVZX88pe/lEGDBkmPHj0i7rd3715p165dyH16W++Pha8DbqqCl5/7wgi38ePqZvAz+m8B5yorK5OWLVtW327SpOGgrL242ke7du1aSQVfB1w/BUA7J5cBQLRjH2MQ8C+rds8VJ2rZsmVIwG1IXl6evPHGG/Luu+9Kx44d6903NzdX9u3bF3Kf3tb7Y0EPbgR2D7JWtyfYgeotAAD+FQgETLhdsWKF/PnPf5azzz67wWMGDBgghYWFIffpJDO9PxYE3CSHNbuDsl3nSri1Bn+PAAC3mjBhgrz44ouydOlSsxau9tHqduLEiep9xowZI1OmTKm+PXHiRLP6whNPPCFbtmwxvb3vv/++Ccqx8G3APdG9vfiFEy7NCwDRov8W8IZnn33WrJwwePBgad++ffW2fPny6n127dole/bsqb49cOBAE4ife+456dmzp7zyyivy2muv1TsxLRx6cB04iLu9PYGqo7WYbAYAcGuLQkPWrFlT575rr73WbInwbQUXoawKUIRbAABgNyq4YXipWua19oRYq9tee/6Al3lp7AVgLwKuTRPMnNSe4PTqbSJ/J8Fj3R50aVOAX9B/C8AKtCjAseFWw6lVgd/KnwU40bx586Rz586SlZUl/fv3l+Li4nr3P3TokJnhrBM+dJH2888/X1auXJmy8wWAZCLgelg0VUunVj+SFUYJufAinZE8adIkc7nM0tJSM/NYr92+f//+sPtXVFTI5ZdfLp999pmZobx161ZZsGCBdOjgrB56PpgCiBctCg7i1oHcyuptKv4O9DHc3rIA1FRQUCDjx4+XcePGmdvz58+XN998UxYtWiSTJ0+us7/e/49//EPWrVsnjRs3Nvdp9TdVrur4CxELxg2nfkAHYD8quEkST/+tG7kt3NrxWFbx0u8NrKPV2JKSEhkyZEj1fWlpaeZ2UVFR2GP+8Ic/mKsCaYtCu3btzPqSs2fPlsrKyrD7nzx5UsrLy0M2q3+Xa4+ZbnyNAn64TK9bEHA9Girc1J5g19eQvIHCCw4ePGiCqQbVmvS2XjEonB07dpjWBD1O+26nTZtmrhr04IMPht0/Pz9fsrOzq7dOnToldM6Rxh63jrcAnIeA65DQ58awZcWbkd3P2+7HB+xQVVUlbdu2NVcK6tOnj4wcOVKmTp1qWhvC0cto6tWIgltZWVnKzxkAYkEPboqlokLhluqtU8IlPblws5ycHElPT5d9+/aF3K+3c3Nzwx6jKydo760eF3TBBReYiq+2PGRmhr42dZUF3VI1DkYzNjhhDAPgXFRwk8APA2+iQd0p4dap5xMJX+GiNg2jWoUtLCwMqdDqbe2zDWfQoEGybds2s1/QJ598YoJv7XALAG5EwPVhuEo0gHst3AJup0uE6TJfL7zwgmzevFluv/12OXbsWPWqCmPGjDFtBkH6/+sqChMnTjTBVldc0ElmOunMLnx4A2AlWhRSOOA6pT3BTk4Ot7QqwK20h/bAgQMyffp002bQq1cvWb16dfXEs127dpmVFYJ0kthbb70ld955p1x00UVm/VsNu/fcc48t568fuoPjox3tCQ2NzX74Vg7wGgIuUhbSnRxuAbfLy8szWzhr1qypc5+2L6xfv178KpaxLLgvQRdwD1oUfIZVItx/nnyVCz9Jxrcq8b6G9Dhef4A7EHAtDn+x/gwrA1Wyv16Pd2B3Q2h08/kCSG1AJeQCzkfATREnDIh2VG8JiwDqXKY3gmT237pxTgWA+BFwI2DwSvzvg3ALwAnjbbJ+Pu8TgHMRcG3klvYEP4Zbt58/AHeHZ/jbqt1z7T4F1yPgWsjJM2ydfG4A/KN2IEzm2ET4BPyLgPt/GGTD82P11g3Pw82/U0C0r7fa30w59YM6r0fAeQi4NnFLe4KXQiEA/0h16CTkAs5CwI1hoHJreEvVrGO3/v0AsE8ygiFhEwABN8khkoHW3QjtgLs+qNs55jLeA85he8CdN2+edO7cWbKysqR///5SXFxc7/6HDh2SCRMmSPv27aVJkyZy/vnny8qVK8VNUtmeQPUWAAD4ja0Bd/ny5TJp0iSZMWOGlJaWSs+ePWXo0KGyf//+sPtXVFTI5ZdfLp999pm88sorsnXrVlmwYIF06MCnZjsRbgHEI/gBPJoJZm6poDrhHACIZNj54AUFBTJ+/HgZN26cuT1//nx58803ZdGiRTJ58uQ6++v9//jHP2TdunXSuHFjc59Wf1PBTyEulgHaT38vAADAHWyr4Go1tqSkRIYMGfKvk0lLM7eLiorCHvOHP/xBBgwYYFoU2rVrJz169JDZs2dLZWVlxMc5efKklJeXh2z1fY1v5afvcD/LDe0JcEeIp1IEv4hmLOP1AMARAffgwYMmmGpQrUlv7927N+wxO3bsMK0Jepz23U6bNk2eeOIJefDBByM+Tn5+vmRnZ1dvnTp1svy5eClIUr0FkCyDr3wkKaHUaeHWaecD+JHtk8xiUVVVJW3btpXnnntO+vTpIyNHjpSpU6ea1oZIpkyZIocPH67eysrKYnrMYIjTaqmT1pu1G+EWAAA4lW09uDk5OZKeni779u0LuV9v5+bmhj1GV07Q3ls9LuiCCy4wFV9tecjMrBu6dKUF3RJRM+Q6RTLaE6g6AEi1RCeYMW7Ba1btnmv3KXiCbRVcDaNahS0sLAyp0Opt7bMNZ9CgQbJt2zazX9Ann3xigm+4cGunZPff2skrzyMWfnzOgBO4tQWM4A34uEVBlwjTZb5eeOEF2bx5s9x+++1y7Nix6lUVxowZY1oMgvT/11UUJk6caIKtrrigk8x00pldg6cTB99kVm8JegCsYMXYSYgE4MiAqz20jz/+uEyfPl169eolGzdulNWrV1dPPNu1a5fs2bOnen+dIPbWW2/JX//6V7nooovkF7/4hQm74ZYUs2KwJMwBAADE791335Xhw4fLmWeeKY0aNZLXXnutwWN0BSydY3XWWWeZNlNdElaXinXNOrgqLy/PbOGsWbOmzn3avrB+/XpxMruXB4sV1VsAqZZo9dUN1Vs9Ryd+ywekkn4zrxfy+tnPfiY/+tGPojrmuuuuM3OyFi5cKOeee64pdtZsT3VFwLVTLAOPkyaY1YfBFIBXJpgxngHud9VVV5ktWvpN/jvvvGOWhm3dunXcF/Vy1TJhsB7VWwBu44bqLeBl5bUuoKUtBVbRi3r17dtXHn30UenQoYOcf/75ctddd8mJEydi+jm+ruAmGuiirS44uT0hGoRbAAC85djnLSUtKyumY6q++sr8WfuiWTNmzJCZM2dacl5auV27dq1kZWXJihUrzIXB7rjjDvn73/8uixcvjvrnUMH10BV1Yv06jypIbAj6cLJ58+aZr/H0TaF///5SXFwccd8lS5aYyR41Nz3ODRi3APuVlZWFXESr5opXidJeWx2Tfve730m/fv3k6quvloKCArPiVixVXN8G3KYf/2t1htoIMv/C3wXgfMuXLzfLLmoVpbS01EzoGDp0qOzfvz/iMS1btjQTN4Lb559/Lk7ilf5bAjm8qGXLliFbohfUqkmvbaCtCdnZ2SEX9QoEArJ79+6of45vA26qOLU9we5B90iXqogbgNhodWP8+PFmDfHu3buby5c3a9as3mV1tEKiV40MbsHlGZNl8JWPNBhe3TKZF0Dy6EW9vvzySzl69Gj1fXrtg7S0NOnYsWPUP4eA6xHJqHZYXb2NNsQSdIHo6WXKS0pKZMiQIdX36RuB3i4qKop4nL556BqT2kt3zTXXyEcffRRxX51AUntSSaIfsOMZX+z+YA4gdjrW6HUOdFM7d+40/63XOlDa3qAX9gq64YYb5IwzzjAf2D/++GOzju7dd99tlhlr2rRp1I9LwI1z0AwXKN0y+EZznskIt/EcQ9AF6qcTMCorK+tUYPX23r17wx7TtWtXU919/fXX5cUXXzQ9bwMHDoz49V9+fr75ujC41Z5gAsAaq3bPFa95//33pXfv3mZT2k6l/60X+VLaIhUMu6p58+by9ttvy6FDh8xqCjfeeKO5UMTcubH93bCKgkvU99Wdk3vVrAio+jNabOezGGAVvWCObkEabrXH7Te/+Y3MmjWrzv5aYdE3pSCt4CYz5Lq5gAAg1ODBg03/bH2TXmvr1q2bCbmJIOAmkVsnaFl13lRfgeTLycmR9PR0c9WfmvS29tZGo3Hjxqaism3btrD/v04gsXISCQAkG2WxGK+sE4lbqgupOM9ktBYQloHwMjMzpU+fPlJYWFh9n7Yc6O2aVdr6aIvDpk2bzOxlu3h5gplb3h8AL6GCGwcntQRYfS6JVm+TGURpVQDC0/aBsWPHmn41XTdyzpw55vrvOklD6QQOXXZHe2nVAw88IN/5znfMNd61z+2xxx4zy4Tdcsst4vaAWHsM83JwBhAZAddHy4Mlu4pAlRWwx8iRI+XAgQNm0oZOLOvVq5e5nntw4plO4NCVFYL++c9/mmXFdN/TTz/dVIDXrVtnlhize3yM90N7pJ+p9xNyAf/xfcDlqyNrQnmqwi1VXCC8vLw8s4WzZs2akNtPPvmk2fwydhFyAf8hKVgQ8OwKybFUOpJ5jlRuASRLQ2NXtOO2Wyf9AogPAdcn7QnJOmc7wi2BGvCWeMe5WMctQi7gHwTcGAdeJ00ws7t6S9AEYKVYxlfCKoD6+Drg1gx+bhssrQza8Tx3u8Ot3Y8PwL4P54mM124b6wHEx9cB1wrJ7G1NRXsCgz2AZBt85SN2nwIAnyHgejwwJiOAUz0FYNcY6bQxFrDKqt1z7T4FT4kr4Oqi4YsXL5bt27eLW53oHvsVe5zSf2vVebixNQGAPyTz2zFCMuB9afFeGlKviHPeeedJp06d5Kc//ak8//zz8umnn4qXuH3dRKvfIJwWbp12PkCs3nvvPTN+6iV1v/jimw+uv/3tb2Xt2rXiF8FxNpoP7gRTAEkNuBpmP/nkEykrK5NHH31UmjdvLk888YR069ZNOnbsKH5RO0C6aXmwWM+VMAlY67/+679k6NCh0rRpU/nggw/k5MmT5v7Dhw/L7Nmz7T49xyHcAkhZD65e4vGMM84wf7Zq1UoyMjKkTZs24jZuGjijbU/gCm2Asz344IMyf/58WbBggTRu3Lj6/kGDBklpaan4GeMXAFsC7r333isDBw404Xby5Mny1VdfmT/1uuZaiYCzUb0F7Ld161a59NJL69yfnZ0thw4dsuWc/MRNhQ0AscuI4xh5+OGHTaV2xowZ8qMf/UjOP/988bqGKqdOak+wsvpBuAWSIzc3V7Zt2yadO3cOuV/7b8855xzbzsuJCKMAUlLB1Srt1KlTpbi42Hyd1qFDB7nhhhvkueeeM725fmDHV2hWrJ4QyxsF4RZInvHjx8vEiRNlw4YN0qhRI/nyyy/ld7/7ndx1111y++23ixfVHn/CTTCjPQGAbRXcnj17mu0Xv/iFuf2///u/8uSTT8qECROkqqpKKisrxS2oDLibhvAW21nOGe6jbV06Xn7/+9+X48ePm3aFJk2amID785//3O7TcwzGaAApC7iBQMBUcdesWWM2/UqtvLxcLrroIrnsssvEC+xaIiyZ7QlUb73FKesyIz5atdVvwu6++27TqnD06FHp3r27WZUGAGBDwG3durUZjLWKq4FWv2q75JJLzEoKfpDM5cGcEGYIt0Dq6LriGmxRF9VbACkNuC+++KIJtC1bthQ/DJ5uqZRZVb0FAKsMvvKRqPel/xaArQF32LBhlp0AGkb1FoBXuaWAAMBdmJ2TIKcsD2ZF9ZZwCyBVGhrv+NYJfrJq91y7T8FzfBtwD52bGdfAy1doABC7hgIrYysAK/k24Lrl67NEH5/qLQC3oXoLIFEEXI+3J0SDcAsAALyEgOthVEG8sTayE79ZALzwe8wYCXgXATcGqe4Ri2bwp3oLQM2bN086d+4sWVlZ0r9/f3Mp9WgsW7bMXHRixIgRYpea45gXQ6fdQR7wIwKuA6p0yXhML75JAAhv+fLlMmnSJJkxY4aUlpaai/AMHTpU9u/fX+9xn332mbk0sK5r7tdvQwB4EwE3zk/dTgiQVG8BqIKCAnNFyXHjxpmros2fP1+aNWsmixYtinhMZWWl3HjjjXL//ffLOeeck7Rzc8JYCcB/CLge/EormjcUwi3gDRUVFVJSUiJDhgypvi8tLc3cLioqinjcAw88IG3btpWbb765wcc4efKklJeXh2xeCMNUkgHvIuBGiTUaATjRwYMHTTW2Xbt2Iffr7b1794Y9Zu3atbJw4UJZsGBBVI+Rn58v2dnZ1VunTp0SDq7BD/GMrQCSgYDr0uXBEnlToHoL+NeRI0dk9OjRJtzm5OREdcyUKVPk8OHD1VtZWVlUxw2+8pEEzxaA27377rsyfPhwOfPMM82E1tdee63e/V999VW5/PLLpU2bNtKyZUsZMGCAvPXWWzE/bkYC5wyXtifA3ZiRjZo0pKanp8u+fftC7tfbubm5dfbfvn27mVymbzhBVVXffOjNyMiQrVu3SpcuXUKOadKkidkAIFbHjh0zE19/9rOfyY9+9KOoArEG3NmzZ0urVq1k8eLFZrzasGGD9O7dO+rHJeBaGCKCldEW25NbGKd6m3r06sGpMjMzpU+fPlJYWFi91JcGVr2dl5dXZ/9u3brJpk2bQu677777TGX3qaeeirn9wCp8OIdfrdo9V7zsqquuMlu05syZE3Jbg+7rr78u//3f/03AtTrExBooNURGE3KtDk0NvUEQbgFv0iXCxo4dK3379pV+/fqZNwitmuiqCmrMmDHSoUMH00ur6+T26NEj5Hitkqja9ydjTIpn3Ks9diW7iACgfrUnmibzWx79wK4fwFu3bh3TcQRch1UZGqocMyEDQG0jR46UAwcOyPTp083Esl69esnq1aurJ57t2rXLrKzgNNGMZ+E+mEdbRHACWorgVM13pEl6k9heR5Unv9m/9jc9ugb3zJkzJRkef/xxOXr0qFx33XUxHUfAtYjd1VGqt4C/aTtCuJYEtWbNmnqPXbJkSZLO6puKbe3xqXboizR+MW4BzlRWVmYmgAUlq3q7dOlSs1a3tijosoaxIODa9Imbnk4AiIxwCzhXy5YtQwJuMuhlxG+55RZ5+eWXQ9b5jpY7vuNxSP9tsqsM8bYnUL31D77uhNskq62LcQ3wrpdeesnMIdA/hw0bFtfPcETAnTdvnnTu3NlMfujfv78UFxdHne51TbXgzGE7OH2Qdfr5uQHVdiC1/beMW4B3HD16VDZu3Gg2tXPnTvPfOjcguM62ToSt2Zagt5944gmTCXVegW66BrerAu7y5cvNDGBtUC4tLTVrpQ0dOlT2799f73G6juNdd90ll1xyifhBvNVbKzU/+3DYzU5umWgCIDyWBwO87f333zfLewWX+NLMp/+tk2LVnj17qsOueu655+TUqVMyYcIEad++ffU2ceJEd/XgFhQUyPjx46uXs5k/f768+eabsmjRIpk8eXLYY/SylDfeeKNpPH7vvffk0KFD4vaKYDK+eraqCtJQiA3+/0d3ZlvyeAD8+22Il6q3tBQBIoMHD5ZAIBD1JNeGJsVGy9byV0VFhZSUlIQ0D+tSNnq7qKgo4nEPPPCAmU138803J/X8oum/tVOyz8kJFVr8C2+W8AJ+jwGkgq0V3IMHD5pqbHCtxiC9vWXLlrDHrF27VhYuXFjdy9GQkydPmi3S4sS1qwqxDL6RKg1Wf20ez9q3iVZB4gm2egxVXADx8lL1FoC9XNXAqFeyGD16tCxYsMBcfz0aeuWe7Ozs6s2uy1Ba3Z6QrOptolVbr1V8mWAGxGfwlY80+KHdzm/GeG0D3mZrBVdDanp6uuzbty/kfr2dm5tbZ//t27ebyWXDhw8PuYSbysjIkK1bt0qXLl1CjtHZedrQXLOCG2vIDTcIp6rSkMrqrdfCKQD7pCK8uumKZgBSy9aRITMzU/r06SOFhYUhgVVvDxgwoM7+3bp1k02bNlUvN6HbD3/4Q/nud79r/jtccNWrawQXJI5lYeJUXRLX6uqtE8ItQRmA39sT6DVGNFbtnmv3KXiW7asoaHV17Nix0rdvX+nXr5/MmTNHjh07Vr2qgq6F1qFDB9NqoOvk9ujRI+T4Vq1amT9r35/MASnegTjWr8RSFbIJpM7HmyUAAC4KuCNHjpQDBw6Y9dB0Id9evXrJ6tWrqyee6dpourICklO9dXO45atJwD30A37ND2o1xzKvVW8B2M/2gKvy8vLMFk5D66HVXj/N6qpqrH1ksYQuu6tybg63qcAkFMB6qfpmCoC/UQJL0QBsRXuCldXbZIdbwjMAALCLrwNuPOvfeuGrNMKnu9hd6QeSyQtjam28ZgH7+TrgOnUATHb1FgDsZPeVIWk/AryPgBvDINxQkIzUf+ukwZTqbXSc9G8GAABi49uA22pbhSMnQMT6+LFUb70UbllBAXCf2uNbot8+MQ4AiITRwQX9WZG+zvNruPUTevngxSXCvMwvzxNwOgJulANSqi7ukAyEWwBOYXf/LQB/IOCmeP1bqyaXRRu4Cbfu/FACuFnPiU/afQoAfI6ACzgUX3XCy+xa/YUPsIA/EHCTOMEsloGU6m30mFgCuI/dE3hTgQ+liMWq3XPtPgVP831SCDcgxbo8WCKPlSxeDbcAwps3b5507txZsrKypH///lJcXBxx31dffVX69u0rrVq1ktNOO0169eolv/3tb1N6vgCQTL4PuE5gdfWWcBs/vr6EGy1fvlwmTZokM2bMkNLSUunZs6cMHTpU9u/fH3b/1q1by9SpU6WoqEj+9re/ybhx48z21ltvJeX8gh/ug+MaF6cBkGy+DrhWrZ4Q7itzghISwVediEVBQYGMHz/ehNTu3bvL/PnzpVmzZrJo0aKw+w8ePFj+/d//XS644ALp0qWLTJw4US666CJZu3Ztys8dAJLB1wE3lf1hsQQWqrcAolVRUSElJSUyZMiQ6vvS0tLMba3QNiQQCEhhYaFs3bpVLr300rD7nDx5UsrLy0O2aDnpw34yz4UPpYCzEHBtXqMx2mBNuE3+BDMnvRED0Tp48KBUVlZKu3btQu7X23v37o143OHDh6V58+aSmZkpw4YNk6efflouv/zysPvm5+dLdnZ29dapUyfLnwcAWMm3Abfpx3ss+TlWtyfEG7D9EG79gkoQUqFFixayceNG+etf/yoPPfSQ6eFds2ZN2H2nTJliAnFwKysri/uDPP23AFIhIyWP4lLJXj3ByuotAH/KycmR9PR02bdvX8j9ejs3NzficdrGcO6555r/1lUUNm/ebCq12p9bW5MmTczmpG/HnLRcIB9KAedxzghhI6esz0j1FkCstMWgT58+po82qKqqytweMGBA1D9Hj9FeWwDwAiq4Se6/TXb11k/h1g/9t1SCEA9tLxg7dqxZ27Zfv34yZ84cOXbsmFlVQY0ZM0Y6dOhgKrRK/9R9dQUFDbUrV6406+A+++yzSfmddkIRwSmvcQCpQcB1yPJgqZ7cBsA7Ro4cKQcOHJDp06ebiWXacrB69erqiWe7du0yLQlBGn7vuOMO2b17tzRt2lS6desmL774ovk5yUS7FYBUIeDawA/V26M7s+0+BcBX8vLyzBZO7cljDz74oNmSoefEJ6vHLi0COKF6m0x86wI4k+97cK0cfGtXb6Md+OKp3jo53CYD7QmA+16zwddWfWOcjmV+G88AJJ/vA25QzQHYCV+jOeEcACCZagZbQi78ZNXuuXafgucRcFNcUaxdMaZ6C0X1Fl5X+0N7IuNYrONvsr6l4XULOBcB1+b2hFirt4Rbb7YnAAAA6/g64KZ68kM01VtaE5y9oDuA6EX6oM8HdcBf5s2bJ507d5asrCzp37+/FBcX17u/LnXYtWtXs8qLXhr8zjvvlK+++iqmxyQ5xNF/21DgStbXVm55U3DLCgpOqd7yNSe8Sj/Ue3UJRF63QHSWL19u1uqeMWOGlJaWSs+ePWXo0KGyf//+sPsvXbpUJk+ebPbXKywuXLjQ/Ix7771XYuHbgHuie/uUBqVEq7duCbcAUN/YluqxzCkfZAG/KigokPHjx5sLz3Tv3l3mz58vzZo1k0WLFoXdf926dTJo0CC54YYbTNX3iiuukFGjRjVY9a3NtwEX7kB7AuAutFkB3ldeXh6yRbrMd0VFhZSUlMiQIUOq79OLzujtoqKisMcMHDjQHBMMtDt27DBXW7z66qtjOkffX+gh1q/P4mlPoHrrPE6p6vA1J7wqWa8xJ3zo5XULL2i1rUIyMmJ7PZ069c3rWvtia9J2gpkzZ9bZ/+DBg1JZWVl9VcUgvb1ly5awj6GVWz3u4osvlkAgIKdOnZLbbrst5hYF3wfcRCsPyQ5KhFsAXsF4BnhDWVmZtGzZsvp2kyZNLPvZeuXF2bNnyzPPPGMmpG3btk0mTpwos2bNkmnTpkX9cwi4Kf5U7/WVE6ycYOaESk0yUQWC19k9wcwp39QAXtOyZcuQgBtJTk6OpKeny759+0Lu19u5ublhj9EQO3r0aLnlllvM7QsvvFCOHTsmt956q0ydOtW0OETD2wnCAYNvIkuRUe1IDt70gNTx0gd4xQdTIHqZmZnSp08fKSwsrL6vqqrK3B4wYEDYY44fP14nxGpIVtqyEC0quAksDxZrUIqleku49TbeJOFFPSc+KdIl/LcvjGmAP02aNEnGjh0rffv2lX79+pk1brUiq6sqqDFjxkiHDh0kPz/f3B4+fLhZeaF3797VLQpa1dX7g0E3GgRci0QzuczrlQ0rJas9geot4F5eb1sCvGjkyJFy4MABmT59uuzdu1d69eolq1evrp54tmvXrpCK7X333SeNGjUyf37xxRfSpk0bE24feuihmB7XtwH30LmZEv3ngNRKpNJxxb9tDXv/H3d1lWRzywUeAHhfMj7M8s0LrLBq91zxm7y8PLNFmlRWU0ZGhlmVQbdE8HE4SrG2JzS0NJjVrQkabCOF2+D/7xZer97yJgm/oT0BQKoRcC1oFbAzsDQUbGvvCwDJph8m9UO9l9qw+GAKuAsBNwlSVb11UmClPSE6vEnCy5wSaJ3ybQ0A+xBwE2xPsCqwpCrcOikU+7E9AUB8mGAGIBa+HzGSXXGItnrrpZCK8Kjewo+i/fDu5G+BeO0C7uP7gBurWCeXWV29dWK4teqNieotAACwgm+XCYv34g6xfKqPpnobbbh1YrBF9KgAwW/s6sflAy0ARQXXIZflrQ/hNn682QGppUWBVFwGPVX4cAq4k28D7tFzYq8uxDK5zKrqrdPDrdPbE5yAN0ikwrx586Rz586SlZVlLm9ZXFwccd8FCxbIJZdcIqeffrrZhgwZUu/+dq9/6+XxAUByMGpYNKjWV71NpDUB8aN6C79Yvny5ud67XvmntLRUevbsKUOHDpX9+/dHvHLQqFGj5C9/+YsUFRVJp06d5IorrjCXxXQzq1/zfDgF3IuAmwRWfT3n9OqtVbxcneENEqlQUFAg48ePl3Hjxkn37t1l/vz50qxZM1m0aFHY/X/3u9/JHXfcYa4J361bN3n++eelqqpKCgsLEz4XPlgCcALvJguL1deekIzqrRvCrZOX9eFNFn5RUVEhJSUlps0gKC0tzdzW6mw0jh8/Ll9//bW0bt067P9/8uRJKS8vD9lqG3zlI4674EMi+HAKuBsB14HVWzeEW6tQvQUSc/DgQamsrJR27dqF3K+39+7dG9XPuOeee+TMM88MCck15efnS3Z2dvWmLQ0NjX2xtGHV92E52jGCD7UAavJuukhQpEE12dVbt4RbqreANzz88MOybNkyWbFihZmgFs6UKVPk8OHD1VtZWVnE8c5LKygAVlu1e67dp+Abvl4H1+rA1NDA7oRw+8ddXcUpqN4CicvJyZH09HTZt29fyP16Ozc3t95jH3/8cRNw//SnP8lFF10Ucb8mTZqYzS94/QLu54iE4cTlbaIRa/XWC5Vbp3NC9ZY3R6RSZmam9OnTJ2SCWHDC2IABAyIe9+ijj8qsWbNk9erV0rdvX3EzJ7zuATiL7QHXicvb1KwsRrv2bSLVW7eFWyvaE7xcvQVSTcdQ/fD/wgsvyObNm+X222+XY8eOmVUV1JgxY0ybQdAjjzwi06ZNM6ssaHFBe3V1O3r0qDiJHeMEH1ABb7A9ZThpeRur1K7esuatf6o4vDnCDiNHjjTtBtOnTzdj48aNG01lNjjxbNeuXbJnz57q/Z999lmz+sJPfvITad++ffWmPyPRQKrjH2MeAF/34AaXt6lZWUjG8ja6BYVb3ibRyWU1q7deb01wavXWCeEWsFNeXp7ZIn3zVdNnn30mXhlTrHzt8wEV8A5bK7hOWd4mkYHTTa0JTppg5kW8OQIA4Ay2tyg4YXkbqyaXxVK9dVvlVlG9jYxwC4Qf5+Id61Ldf8trGPAWW1sU3LK8TbwDX6TqrRvDrRWYWAb4k455yfgGyQkfbgE4k62Jw2nL20RaPSGe3lunTbJI9M3FqRd2cMIbHJUfoO7Y56YP8ryGAe/JcMLyNmPHjjVBtV+/fjJnzpw6y9t06NDB9NIGl7fRmcJLly6tXt5GNW/e3GxOH/jcNOhbidYEwPvMCgoe+OAMwP0ynLC8zYEDB0xo1bCqS9zUXt5GV1YIt7xNTbqO7syZM5MawBKt3ro13PImBCCVohmPnfABF4Bz2R5wnbi8TUMDZ30rJzgx3Nq9egLVW8Dbek58UqSLOz/M8zoGvIlZP1EOfImsnOCmwd7q6i3hFvC+WNf/BoBkI+BGObksHKdPLHNC9daLCLdAeFaMgalsT+C1jFRatXuu3afgKwTcGKu3DV3YoTaqt96r3gKIbcwL90Gb3n4AyUTAjUGsE8vsDreJVG8Jt+FR8QHsXeea6i2AaPg+4MbbnuD0cOs1hFvA2dzUh8trGfA+XwfcSBWHcO0J9VVvnchr1Vu78YYIJJ8Xxw4A9vDtaNJ8R+hTp3rr3HDrhOotAPvHGSvGAj6sAv7g24Ab79JgNau3Tg23Xlo5wQnhljdEoGE6Hjph/APgPPPmzTNXn83KypL+/ftLcXFxVMctW7ZMGjVqJCNGjIj5MQm4UYh25QS3D+5Oq94SbgH/aGj8oHoLuNPy5ctl0qRJ5oqzpaWl0rNnTxk6dKjs37+/3uP0wl533XWXXHLJJXE9LgE3wsDZUO+tE9e8TaR6S7itizdDAAASU1BQIOPHj5dx48ZJ9+7dZf78+dKsWTNZtGhRxGMqKyvlxhtvlPvvv1/OOeecuB6XgGtRoHF79TYRXpwYQrgFkvOh2671b3lNA9YpLy8P2U6ePBl2v4qKCikpKZEhQ4ZU35eWlmZuFxUVRfz5DzzwgLRt21ZuvvnmuM8xQ3zOiuqtU8KtHdVbL04q440Q8F57AoBQTT/eIxlpsV286lTVN6/FTp06hdyv7QczZ86ss//BgwdNNbZdu3Yh9+vtLVu2hH2MtWvXysKFC2Xjxo2SCF8H3JqDZn2hhtaE1LH7jYxwC8THKR/0w+F1DVirrKxMWrZsWX27SZMmlvzcI0eOyOjRo2XBggWSk5OT0M/ydcANJ9zKCVYN6qNPX2f+/O0/B4oT0HcbijdBwHkSHRd4XQPW03BbM+BGoiE1PT1d9u3bF3K/3s7Nza2z//bt283ksuHDh1ffV1X1TZExIyNDtm7dKl26dInqHH0bcFttqxDJyIo4AAbbE6xoTQgG22QFXTuWBSPcAlDnzH0i6m+2wn2o9mIPP4BvZGZmSp8+faSwsLB6qS8NrHo7Ly9PauvWrZts2rQp5L777rvPVHafeuqpOq0R9fFtwA0nXO9tomqHW6t5oe+WcAsgHKq3gPtNmjRJxo4dK3379pV+/frJnDlz5NixY2ZVBTVmzBjp0KGD5Ofnm3Vye/ToEXJ8q1atzJ+172+I7z86x9t721D1VoNtssNtvAi3/8IbIPy4kPpHH30kP/7xj83+uoi6vuEkg5cuOgMkYtXuueJXI0eOlMcff1ymT58uvXr1MpPHVq9eXT3xbNeuXbJnzx7LH9fXFdxow0084TYV4nnzINz+C+EWXltIXdeX1HCrgVUXUtd+NV1qp7bjx4+btSWvvfZaufPOO22bYFbfmEL1FvCOvLy8sC0Jas2aNfUeu2TJkrge0/cV3EjtCTWrtzURbq1BuAXsW0j929/+tjz22GNy/fXXWzb72UkrtvD6BpDm5/Xfolk5IZZlwWIJt4lMMCPcJoY3P3hJvAup243JZQCSydctCrUlUr11ar+tItz+C+EWXhPPQuqx0qsU1bxSkV65KJkSGSd4jQNQvv8IXd/KCTWrt04Jt7FWb50SbvUNi3ALuJPObs7Ozq7eolmqp+ZYlar2BF7jAIJ8H3BrC1Zvo21NiCfcxtue4OZwayd90+OND14V60Lq8ZgyZYocPny4etOrGCVydcdkTi4DAPF7wI123dtI1VvCbcPsfrMi2MJPC6kHBRdSHzBggCWPoRPRglcuivYKRqnGax1ATb7twT3RvX2dJx+uemtluI1XqsIt/baA9xdSD05M+/jjj6v/+4svvjBrUzZv3lzOPffcmB5bx8j6rs5YezxKRvWW1zqA2nwbcKOdWGaleKq3bgy3dgdbxRse/LaQ+oEDB8xC6nv37jWLqddeSF1XVgj68ssvpXfv3tW3dRF23S677LIG16RM1kovhFsAVvJ9wK0tWdVbwm1q8GYHv4plIXW9glkgEEjRmQFA6vm6B7e+6q3drQmE29gRbgFnizTeUL0FYDXfVnAPnZsp6bXua2g2cLzhNtbqbSzhlmDLmxxgp/rGxWQuD8brHkB9fF3BjWViGeE2PMIt4F8/7vb/6ox1kcYwK6u3vO4BNMS3FdxYJpYRbusi2ALY88OzUv6YvPYBRMO3AffoOVXV5etor1jm5nDrlaotb26As8bRaMYpq3tvATdZtXuu3afgS74NuEHJaE2IJdwSbKNDsAWc57SzyuM+ltYEAMnk+4BbH8KtM6osvKkBzvS9jp+GHdOSMbmMcQBALDL8XHlIb3YyYvXWKeGWYAvAzcKNR7GOLYwFAGLl24BbkxV9t9GG22RWbb0QbHkjA9zj+lYbzGI88VzIJlqMCQDi4fuAa0XfrdXh1q6qLcEWQKIamlwWyzjDuAAgXr4OuE4LtwRbAG4VzRhHuAWQKhn+nhzROO5wa3fVlmALwMlqj1GEWwCp5NuAW5td4dZPwZY3LcDbrFg9gXECgBUIuDGKJtxaHWwJtQD8UL1lvABgFQJulNVbgm30eJMC/EPHxvrGPsItADsQcKNgRbhNZbAl1AJIteAYF+v4xdgBIBkIuA4JtoRaAH5rTWAMAZAsBNz/C7LBNgWrJpClItimOtTyZgQg2vGOcAvATgTc/2NFsE12qCXQAnCa4JhYc1wj3AKwGwE3SokGWzeEWt50ACSqvjGLMQZAqhBwHRZqUxVoeaMBkKhlh/qbsTCasY4xB360avdcu0/Btwi4FvfWxhpqUxFoeWMBkArhxjPGHwB28G3A/fPu8yS9WZMG97My0CY7zPJGAiDV42hwDCTcAnAS3wbcZLcdJCvM8oYBwGlqj3eMUwDs5tuAe+zzlpKWleW4MMsbAwA3jaOda4x/jF8Awpk3b5489thjsnfvXunZs6c8/fTT0q9fP4nk5ZdflmnTpslnn30m5513njzyyCNy9dVXSyx8G3Cb70iT9Cap75flDQCAl8bRIMY2AOEsX75cJk2aJPPnz5f+/fvLnDlzZOjQobJ161Zp27Ztnf3XrVsno0aNkvz8fPnBD34gS5culREjRkhpaan06NFDouXbgNtqW4VkZCR29TDFoA7Az+No0y1/t/s0ADhYQUGBjB8/XsaNG2dua9B98803ZdGiRTJ58uQ6+z/11FNy5ZVXyt13321uz5o1S95++2359a9/bY6Nlu8CbiAQMH82/vBzyUjLTPjnnbLgnAB4y6mqipDxxsvj6CkLxlHAq8rLy5P+s+MdZ04FKkSq4jgmzPNq0qSJ2WqrqKiQkpISmTJlSvV9aWlpMmTIECkqKgr7GHq/Vnxr0orva6+9FtO5+i7gHjlyxPy5Zt9iu08FgA/Gm+zs6K5w6CaMo0B0srN/47hxJjMzU3Jzc2XN3vhev82bN5dOnTqF3DdjxgyZOXNmnX0PHjwolZWV0q5du5D79faWLVvC/nzt0w23v94fC98F3DPPPFPKysqkRYsW0qhRI/Ei/WSlv3z6PFu2bClew/NzP68/R62o6JuOjjde5LVx1Ou/jzw/bz6/eMeZrKws2blzp6muxkMft/brPlz11m6+C7haGu/YsaP4gb4QvPhiD+L5uZ+Xn6MXK7deH0e9/PuoeH7ee37xjjNZWVlmS7acnBxJT0+Xffv2hdyvt7WKHI7eH8v+kSQ+ywoAAAAI0w7Rp08fKSwsrL6vqqrK3B4wYICEo/fX3F/pJLNI+0fiuwouAAAAUkMnjI0dO1b69u1r1r7VZcKOHTtWvarCmDFjpEOHDmZZMDVx4kS57LLL5IknnpBhw4bJsmXL5P3335fnnnsupscl4HqQ9sJow7cTe2KswPNzPz88R7iH138feX7u5vbnN3LkSDlw4IBMnz7dTBTr1auXrF69unoi2a5du0zbU9DAgQPN2rf33Xef3HvvveZCD7qCQixr4KpGAa+uYwMAAABfogcXAAAAnkLABQAAgKcQcAEAAOApBFwAAAB4CgHXpebNmyedO3c2CzX3799fiouLI+67YMECueSSS+T00083m14Dur793fb8atLlRPQKKyNGjBAvPb9Dhw7JhAkTpH379mYm7fnnny8rV64ULz1HXTqma9eu0rRpU3PVnjvvvFO++uqrlJ0vvEkvH6pjQs2tW7du9R7z8ssvm330d/fCCy909GtNX2O1n59uOl6Es2TJkjr7pmLB/2i9++67Mnz4cHN1Lj03nT1fk86L19n4OhbqWKHvZ59++mnS3lNS+fy+/vprueeee8zv3GmnnWb20SW0vvzyS8t/x/2AgOtCy5cvN+vK6bIhpaWl0rNnTxk6dKjs378/7P5r1qyRUaNGyV/+8hcpKioy4eGKK66QL774Qrzw/II+++wzueuuu0yYd7JYn59eTvHyyy83z++VV16RrVu3mg8tum6gV56jLgkzefJks//mzZtl4cKF5mfoEjFAor71rW/Jnj17qre1a9dG3HfdunVmvLz55pvlgw8+MB+Wdfvwww/Fif7617+GPDddEF9de+21EY/Rq2HVPObzzz8Xp9D1UXW80EAazqOPPipz586V+fPny4YNG0wQ1LGlvg/D8b6npPr5HT9+3JzftGnTzJ+vvvqqGe9/+MMfWvo77hu6TBjcpV+/foEJEyZU366srAyceeaZgfz8/KiOP3XqVKBFixaBF154IeCV56fPaeDAgYHnn38+MHbs2MA111wTcKpYn9+zzz4bOOeccwIVFRUBt4j1Oeq+3/ve90LumzRpUmDQoEFJP1d424wZMwI9e/aMev/rrrsuMGzYsJD7+vfvH/iP//iPgBtMnDgx0KVLl0BVVVXY/3/x4sWB7OzsgBtoRFmxYkX1bX1Oubm5gccee6z6vkOHDgWaNGkSeOmll5L2npmq5xdOcXGx2e/zzz+37HfcL6jguoxW80pKSszXMkG6QLLe1upsNPRTon4V0rp1a/HK83vggQekbdu2puriZPE8vz/84Q/mEoX6laMujK2LXc+ePVsqKyvFK89RF/bWY4JfG+7YscN8LXz11Ven7LzhXfoVtn7de84558iNN95oFpaPRH9Ha/7uKq32RTu+2v3ae/HFF+VnP/uZ+Zo6kqNHj8pZZ51lvs275ppr5KOPPhI32Llzp7lQQM1/n+zsbNNyEOnfx4r3TDsdPnzY/Fu2atXKst9xvyDguszBgwdNsAleASRIb+sLPxra46MvhNqDuFufn34Vo19p69f2ThfP89Owp60JepyGPv36Si9h+OCDD4pXnuMNN9xgPqRcfPHF0rhxY+nSpYsMHjyYFgUkTMOP9p3qlZOeffZZE5K0jenIkSNh99ff0UTGVztpP6f26990000R99E+90WLFsnrr79uwnBVVZX5gLl7925xuuC/QSz/Pla8Z9pF2y70/VpbZrStxKrfcb/gUr0+8/DDD5uJWNqX66SJBfHSF/Do0aNNuM3JyREv0jcgrU7rdbjT09OlT58+pn/6scceMz1lXqC/j1qVfuaZZ8xgvW3bNnM98lmzZplAD8Trqquuqv7viy66yPx+afXy97//veO/8YmVftDX56sFjEj02yDdgjTcXnDBBfKb3/zGvN7gDPot63XXXWcm1WlorY+ffsdjQcB1GQ1xGnL27dsXcr/ezs3NrffYxx9/3ATcP/3pT+ZF4IXnt337djP5Smel1gyEKiMjwzToazXQzf9+OltYq5p6XJC+IWn1Qb9+y8zMFCeJ5zlqiNUPKrfccou5rbOIdTLGrbfeKlOnTg25TjmQCP2qV1ch0Q9R4ejvaDzjq910opiO7ToxKRY6tvTu3Tvi34eTBP8N9N9Dx8Ugvd2rVy/L3zPtDrf6b/rnP/+53uptPL/jfsG7hstomNEKXmFhYUig09s1P5WHm3mqn871K4y+ffuKV56fLoWyadMm2bhxY/WmM06/+93vmv/WHjO3//sNGjTIDFTB4K4++eQTM8A7LdzG+xy1L7x2iA0G+m/mYgDW0P5T/WBcMyDVpL+jNX93la5MUN/46gSLFy823/QMGzYspuP063sdQyP9fTjJ2WefbUJpzX+f8vJys5pCpH+feN8z7Q632lOrH1jOOOMMy3/HfcPuWW6I3bJly8ys0SVLlgQ+/vjjwK233hpo1apVYO/eveb/Hz16dGDy5MnV+z/88MOBzMzMwCuvvBLYs2dP9XbkyJGAF55fbU5fRSHW57dr1y6z6kVeXl5g69atgTfeeCPQtm3bwIMPPhjwynPUWcD6HHUm9I4dOwJ//OMfzUxwndEOJOJXv/pVYM2aNYGdO3cG/ud//icwZMiQQE5OTmD//v1hfxd1n4yMjMDjjz8e2Lx5s/ndbNy4cWDTpk0Bp9JVAf7t3/4tcM8999T5/2o/v/vvvz/w1ltvBbZv3x4oKSkJXH/99YGsrKzARx99FHACfV/64IMPzKYRpaCgwPx3cBUBfT/TseT1118P/O1vfzNj/dlnnx04ceJE9c/QFVmefvrpqMcjpzw/XSnnhz/8YaBjx46BjRs3hrxfnzx5MuLza+h33K8IuC6lv9w6oGlw1SVQ1q9fX/3/XXbZZSbkBZ111lnmhVR704HbC8/PbQE3nue3bt06s1SRDtK6ZNhDDz1klkbzynP8+uuvAzNnzjShVt9sO3XqFLjjjjsC//znP206e3jFyJEjA+3btze/hx06dDC3t23bVu/r7fe//33g/PPPN8d861vfCrz55psBJ9PAqmO6fgCurfbz++Uvf1n9umzXrl3g6quvDpSWlgac4i9/+UvY96vgc9ClwqZNm2bOXcfD73//+3Wet77n1X5/q288csrz04Aa7v/TTY+L9Pwa+h33q0b6P3ZXkQEAAACr0IMLAAAATyHgAgAAwFMIuAAAAPAUAi4AAAA8hYALAAAATyHgAgAAwFMIuAAAAPAUAi4AAAA8hYALAAAATyHgAgAAwFMIuEAEnTt3ljlz5oTc16tXL5k5c6Zt5wQAiRg8eLDk5eWZLTs7W3JycmTatGkSCATsPjXAUgRcAAB85IUXXpCMjAwpLi6Wp556SgoKCuT555+3+7QAS2VY++MAAICTderUSZ588klp1KiRdO3aVTZt2mRujx8/3u5TAyxDBRcAAB/5zne+Y8Jt0IABA+TTTz+VyspKW88LsBIBF4ggLS2tTl/a119/bdv5AACA6BBwgQjatGkje/bsqb5dXl4uO3futPWcACBRGzZsCLm9fv16Oe+88yQ9Pd22cwKsRsAFIvje974nv/3tb+W9994zPWpjx47lDQCA6+3atUsmTZokW7dulZdeekmefvppmThxot2nBViKSWZABFOmTDEV2x/84AdmOZ1Zs2ZRwQXgemPGjJETJ05Iv379zId2Dbe33nqr3acFWKpRgMXvAADwzTq4up537TW+Aa+hRQEAAACeQsAFAACAp9CiAAAAAE+hggsAAABPIeACAADAUwi4AAAA8BQCLgAAADyFgAsAAABPIeACAADAUwi4AAAA8BQCLgAAADyFgAsAAADxkv8PlTMyTSsdEiUAAAAASUVORK5CYII=",
                        "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, 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, 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": 5,
            "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": 6,
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "Total trajectory points: 2500\n"
                    ]
                },
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAG2CAYAAACkgiamAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAANuBJREFUeJzt3QtUVWXex/E/FwEv4SVGUVLRrIxMbUBQy6yJYs3rZLamiVomLFdjV82iTK0EtSbSGoamGCkns3RmZHrL6q1ebKKsHFEmzJV5IXVS8ALilCKYYHDe9TzznjMcbp5z4Jyz9z7fz1o72PvsvdnbI52f/+eyg2w2m00AAAAMItjfFwAAANAc4QQAABgK4QQAABgK4QQAABgK4QQAABgK4QQAABgK4QQAABgK4QQAABgK4QQAABgK4QQAABiKIcJJXl6exMbGSkREhCQlJUlJSUm7+15zzTUSFBTUapkyZYpPrxkAACvIc+MzePXq1a0+f9VxLe3evVumTp0qvXv3lp49e8q4ceOkvLzcPOGkoKBAMjIyJCsrS7Zt2yZjxoyRlJQUOXbsWJv7v/XWW3L06FHH8vXXX0tISIj86le/8vm1AwBgZgVufgYrkZGRTp/DBw8edHp9//79ctVVV8nIkSNl48aN8tVXX8miRYvaDDHtCfL3g/9USlOJ6sUXX9TrTU1NMnjwYJkzZ44sWLDgnMfn5uZKZmam/gNS6QwAAHjnM1hVTh588EE5ceJEu+e87bbbpFu3brJmzRrxVKj4UUNDg5SWlsrChQsd24KDgyU5OVmKi4tdOscrr7yi/yDaCyb19fV6sVN/8N99952cf/75uhwFAEB71L/fT506JYMGDdKfT95y5swZ/ZnYFdcb1OKzLTw8XC9d9RlcW1srQ4cO1Z+nP/3pT+Xpp5+Wyy67TL+mtr3//vvy6KOP6grMl19+KcOGDdM/Y9q0aW7diN8cPnxYVW1smzdvdto+b948W2Ji4jmP37p1qz5efW1PVlaW3oeFhYWFhcXTpaKiwuYtP/zwgy3qJ8Fdcp29evVqtU19DnbVZ7Da97XXXrN9+eWXto0bN9p+8Ytf2CIjIx1/PkePHtXn7NGjhy0nJ0fvl52dbQsKCtL7u8qvlZPOUlWTyy+/XBITE9vdR6U11Z5md/LkSRkyZIgkTV4goaGut38BQGd133XU35cAN/3Y1CAbq16V8847z2s/Q1Uwjlc3yYYt0dKzl+fVmbraJkkZXykVFRW6X4hdW1UTT02YMEEvdhMnTpRLL71UXnrpJXnyySd15US56aab5KGHHtLfjx07VjZv3iz5+fkyefJkl36OX8NJVFSU7sxaVVXltF2tR0dHd3hsXV2drFu3TpYuXdrhfu2Vs1QwIZwA8JXuXx8WCQ7z92XAQ77oBqCCSa/zOt90pIJJ83Dijc9gO9W35IorrpB9+/Y5zhkaGipxcXFO+6kAs2nTJnOM1gkLC5P4+HgpKipybFOpS603T2ZteeONN3RfkjvuuMMHVwoAnQwmgMGEdeIz2K6xsVF27NghAwcOdJxTdbAtKytz2u+bb77R/VRc5fdmHdXkkp6eLgkJCbp5Ro2+UVWRmTNn6tfT0tIkJiZGsrOzWzXpqM41qmMrABgRoQRGl+HmZ7BqrRg/fryMGDFCj9h59tln9VDiX//6145zzps3T1JTU+Xqq6+Wa6+9VgoLC+V//ud/9LBi04QTdQPV1dV6OHBlZaVum1I3MmDAAP26mrSlZQ9plchUeejDDz/001UDQMcIJjCDVDc/g7///nuZNWuW3rdv37668qL6kzRvxrn55pt1/xIVaB544AG55JJL5M0339Rzn5hmnhNfq6mp0TPWXXndYvqcAPAKgom1OsR+dPQlPZjClX4cnflc2vT1oE71Oak91SRXjTri1Wv1Fb/PEAsAVkIwATrP7806AGAFhBKg61A5AYBOIpgAXYvKCQB4iFACeAeVEwDwAMEE8B4qJwDgBkIJ4H1UTgDARQQTwDeonADAORBKAN+icgIAHSCYAL5H5QQA2kAoAfyHcAIA/49AAhgD4QRAwCOUAMZCOAEQsAglgDERTgAEHEIJYGyEEwABg1ACmAPhBIDlEUoAcyGcALAsQglgToQTAJZDKAHMjXACwBIIJIB1EE4AmBqhBLAewgkA0yGQANZGOAFgGoQSIDAQTgAYHqEECCyEEwCGRCABAhfhBIChEEoAEE4AGAKhBIAd4QSA3xBIALSFcALA5wglADpCOAHgM4QSAK4gnADwKgIJAHcRTgB4BaEEgKcIJwC6FKEEQGcRTgB0GoEEQFcinADwCIEEgLcQTgC4jEACwBcIJwA6RCAB4GuEEwCtEEgA+BPhBIADoQSAERBOgABHIAFgNIQTIAARSAAYGeEECBAEEgBmQTgBLIowAsCsgv19AQC6NpDYFwBwRV5ensTGxkpERIQkJSVJSUmJS8etW7dOgoKCZNq0aU7ba2trZfbs2XLBBRdI9+7dJS4uTvLz88UdVE4AkyOIAPBUQUGBZGRk6PCggklubq6kpKRIWVmZ9O/fv93jDhw4II888ohMmjSp1WvqfB9//LGsXbtWh54PP/xQ7rvvPhk0aJBMnTrVpeuicgKYEBUSAF0hJydHZs2aJTNnznRUOHr06CGrVq1q95jGxkaZPn26LFmyRIYPH97q9c2bN0t6erpcc801OpzcddddMmbMGJcrMgrhBDAJAgkAV9TU1Dgt9fX1be7X0NAgpaWlkpyc7NgWHBys14uLi9s9/9KlS3VV5c4772zz9YkTJ8q7774rhw8fFpvNJp988ol88803csMNN4iraNYBDIwgAgSOdSeSJPzHbh4fX197VkTWy+DBg522Z2VlyeLFi1vtf/z4cV0FGTBggNN2tb5nz542f8amTZvklVdeke3bt7d7HS+88IKulqg+J6GhoTrwrFy5Uq6++mqX74VwAhgMgQRAZ1RUVEhkZKRjPTw8XLrCqVOnZMaMGTpoREVFdRhOtmzZoqsnQ4cOlc8++0zuv/9+3eekeZWmI4QTwM8IIwC6UmRkpFM4aY8KGCEhIVJVVeW0Xa1HR0e32n///v26I+yNN97o2NbU1KS/qgqJ6kSrAshjjz0m69evlylTpujXRo8erSstzz33HOEEMDICCQB/CwsLk/j4eCkqKnIMB1ZhQ62rocAtjRw5Unbs2OG07YknntAVleeff143J505c0bOnj2rm3KaUyHIHmRcQTgBfIAwAsCIMjIy9MiahIQESUxM1EOJ6+rq9OgdJS0tTWJiYiQ7O1vPgzJq1Cin4/v06aO/2rerwDN58mSZN2+enuNENet8+umn8vrrr+uRQaYZrePu5C8nTpzQbVcDBw7U7WgXX3yxfPDBBz67XsBVjK4BYHSpqam6uSUzM1PGjh2rm18KCwsdnWTLy8vl6NGjbp1TTc42btw4PdxYDU9+5pln5De/+Y3cc8895qicuDv5ixr2dP311+vX/vu//1unuYMHDzqSG+BvBBEAZjN79uw2m3GUjRs3dnjs6tWrW21T/VVeffXVTl1TqFEmf1FUSHn//ff15C8LFixotb/a/t133+kJXrp1+/dwK1V1AfyJQAIAXctvzTqeTP6ihiVNmDBBN+uokpNq43r66af1OO32qMlnWk5IA3QWTTYA4D1+q5x4MvnLP//5Tz1fv2rHUv1M9u3bp+frVz2D1SQzbVGdeNQUu0BnEEIAwHdMNVpHDUNS/U1efvllPSxJDYFS0+M+++yz7YaThQsX6n4tdqpy0nL2PKAtBBIACLBw4u7kL4oaoaP6mqjj7C699FKprKzUzURqCFNLakRPV82OB2sjjABAgPc5aT75i5198hfVr6QtV155pW7KaT6Ri3qYkAotbQUT4FzoOwIAxuPXeU5Uc4uao/+1116T3bt3y7333ttq8hfVLGOnXlejdebOnatDiRrZozrEqg6ygLthhEACAMYU6u/JX6qrq/XkL6ppRk0A03Lyl+ZT4Kq+Ihs2bJCHHnpIz9Wv5jlRQWX+/Pl+vAsYHSEEAMwl1GyTv6gmH/W0Q6AjBBIAMC+/hxOgqxBIAMAaCCcwNQIJAFgP4QSmQhgBAOsjnMDwCCQAEFgIJzAkAgkABC7CCQyDQAIAUAgn8CsCCQCgJcIJfI5AAgDoCOEEXkcYAQC4g3ACryCQAAA8RThBlyGQAAC6AuEEnUIgAQB0NcIJ3EYgAQB4E+EELiGQAAB8hXCCNhFGAAD+QjiBA4EEAGAEhJMARyABABgN4SQAEUgAAEZGOAkAhBEAgJkQTiyKQAIAMCvCiYUQSAAAVkA4MTkCCQDAaggnJkMYAQBYHeHEBAgkAIBAQjgxIMIIACCQEU4MgkACAMC/EU78iEACAEBrhBMfI5AAANAxwokPEEgAAHAd4cRLCCQAAHiGcNKFCCQAAHQe4aQTCCMAAHQ9wombCCQAAHhXsJfPb5lAYl8AALCSvLw8iY2NlYiICElKSpKSkhKXjlu3bp0EBQXJtGnTnLbbbDbJzMyUgQMHSvfu3SU5OVn27t3r1jURTtpBIAEAWF1BQYFkZGRIVlaWbNu2TcaMGSMpKSly7NixDo87cOCAPPLIIzJp0qRWry1fvlx+//vfS35+vmzdulV69uypz3nmzBmXr4tw0gyBBAAQSHJycmTWrFkyc+ZMiYuL04GiR48esmrVqnaPaWxslOnTp8uSJUtk+PDhraomubm58sQTT8hNN90ko0ePltdff12OHDkib7/9tsvXFfDhhEACALCSmpoap6W+vr7N/RoaGqS0tFQ3u9gFBwfr9eLi4nbPv3TpUunfv7/ceeedrV779ttvpbKy0umcvXv31s1FHZ2zpYDtENt911EJDQ7z92UAAKB9fOgiCekR7vHxjaf/HUIGDx7stF012SxevLjV/sePH9dVkAEDBjhtV+t79uxp82ds2rRJXnnlFdm+fXubr6tgYj9Hy3PaX3NFwIYTAACsqKKiQiIjIx3r4eGeB57mTp06JTNmzJCVK1dKVFSUeBPhBAAAC4mMjHQKJ+1RASMkJESqqqqctqv16OjoVvvv379fd4S98cYbHduampr019DQUCkrK3Mcp86hRus0P+fYsWNdvoeA73MCAEAgCgsLk/j4eCkqKnIKG2p9woQJrfYfOXKk7NixQzfp2JepU6fKtddeq79XzUnDhg3TAaX5OVW/FzVqp61ztofKCQAAASojI0PS09MlISFBEhMT9Uiburo6PXpHSUtLk5iYGMnOztbzoIwaNcrp+D59+uivzbc/+OCD8tRTT8lFF12kw8qiRYtk0KBBreZD6QjhBACAAJWamirV1dV60jTVYVU1vRQWFjo6tJaXl+sRPO549NFHdcC566675MSJE3LVVVfpc6pw46ogmxqUHEBUeUkNa0oeeDejdQAAHfqxqUE+OvqSnDx50qV+HJ35XIpb92inR+vsum25V6/VV+hzAgAADIVwAgAADIVwAgAADIVwAgAADIVwAgAADIVwAgAADMUQ4SQvL09iY2P1GGj15MKSkpJ29129erUEBQU5Le6MnQYAAMbm93BSUFCgZ6hTT03ctm2bjBkzRlJSUuTYsWPtHqPGbx89etSxHDx40KfXDAAALBxOcnJyZNasWXqq3Li4OMnPz5cePXrIqlWr2j1GVUvU3P32peWjmQEAgHn5NZw0NDRIaWmpJCcn/+eCgoP1enFxcbvH1dbWytChQ/VDhm666SbZuXNnu/vW19fr2feaLwAAwLj8Gk6OHz8ujY2NrSofal3N8d+WSy65RFdV3nnnHVm7dq1+guLEiRPl0KFDbe6vHlakpgW2LyrQAAAA4/J7s4671COX1VMS1cOJJk+eLG+99Zb85Cc/kZdeeqnN/RcuXKifM2BfKioqfH7NAABAzPFU4qioKAkJCZGqqiqn7Wpd9SVxRbdu3eSKK66Qffv2tfl6eHi4XgAAgDn4tXISFhYm8fHxUlRU5NimmmnUuqqQuEI1C+3YsUMGDhzoxSsFAAABUTlR1DDi9PR0SUhIkMTERMnNzZW6ujo9ekdRTTgxMTG674iydOlSGT9+vIwYMUJOnDghzz77rB5K/Otf/9rPdwIAACwRTlJTU6W6uloyMzN1J1jVl6SwsNDRSba8vFyP4LH7/vvv9dBjtW/fvn115WXz5s16GDIAADC/IJvNZpMAooYSq1E7yQPvltDgMH9fDgDAwH5sapCPjr6kB1SoCUC9+bkUt+5RCenheR/JxtP1suu25V69Vl8x3WgdAABgbYQTAABgKIQTAABgKIQTAABgKIQTAABgKIQTAABgKIQTAABgKIQTAABgKIQTAABgKIQTAABgKIQTAABgKIQTAABgKIQTAABgKIQTAABgKKH+vgAAADzxw6gYr/+MH388I3LU6z8GLRBOAACmDQ+wJsIJAAQ4QgSMhnACACZCkEAgIJwAgJcRKAD3EE4AoB2ECsA/CCcALIlgAZgX4QSAIREugMBFOAHgNQQMAJ4gnAA4J0IGAF9i+nogQMOGOwsA68rLy5PY2FiJiIiQpKQkKSkpaXfft956SxISEqRPnz7Ss2dPGTt2rKxZs8bx+tmzZ2X+/Ply+eWX69cHDRokaWlpcuTIEbeuicoJYAEECACeKCgokIyMDMnPz9fBJDc3V1JSUqSsrEz69+/fav9+/frJ448/LiNHjpSwsDB57733ZObMmXpfddzp06dl27ZtsmjRIhkzZox8//33MnfuXJk6dap88cUXLl9XkM1ms0kAqampkd69e0vywLslNDjM35cDtImwARiDerbO34sWy8mTJyUyMtKrn0tx6x6VkB7hHp+n8XS97LptuVvXqgLJuHHj5MUXX9TrTU1NMnjwYJkzZ44sWLDApXP89Kc/lSlTpsiTTz7Z5uv/+Mc/JDExUQ4ePChDhgxx6ZxUTgAfInQAMIqGhgYpLS2VhQsXOrYFBwdLcnKyFBcXn/N4Vdv4+OOPdZVl2bJl7e6nwlJQUJBuCnIV4QToJAIHACOpqalxWg8PD9dLS8ePH5fGxkYZMGCA03a1vmfPng7DRkxMjNTX10tISIj84Q9/kOuvv77Nfc+cOaP7oNx+++1uVZ4IJ0AbCBwAfK3uYKQER0R4fHzTmTP6q2qWaS4rK0sWL14sXeW8886T7du3S21trRQVFek+K8OHD5drrrnGaT/VOfbWW2/VFZYVK1a49TMIJwgYBA4AgaCiosKpStFW1USJiorSlY+qqiqn7Wo9Ojq63fOrpp8RI0bo79Vond27d0t2drZTOLEHE9XPRDX9uNtfh3ACUyNwAIAzFQRcCQNqtE18fLyufkybNs3RIVatz549W1yljlFNPC2Dyd69e+WTTz6R888/X9xFOIEhEToAwPtUk0x6erqeu0SNqFFDievq6vTwYEXNUaL6l6jKiKK+qn0vvPBCHUg++OADPc+JvdlGBZNbbrlFDydWw4xVn5bKykrHMGQViFxBOIFPEToAwDhSU1OlurpaMjMzdYhQzTSFhYWOTrLl5eW6GcdOBZf77rtPDh06JN27d9fznaxdu1afRzl8+LC8++67+nt1ruZUFaVlv5T2MM8JugzBA4DV+HKek6HLnup0h9iD85/w6rX6CpUTuITgAQDwFcIJNMIHAMAoCCcBgOABADATwonJETwAAFZDODE4wgcAINAQTvyM8AEAgDPCiZcRPgAAcA/hpAsQQAAA6DqEExcRQAAA8A3CSTMEEAAA/C9gw8kPcQMlNNTzaYIBAIB3/OdpPgAAAAZAOAEAAIZCOAEAAIZCOAEAAIZCOAEAAIZiiHCSl5cnsbGxEhERIUlJSVJSUuLScevWrZOgoCCZNm2a168RAAAESDgpKCiQjIwMycrKkm3btsmYMWMkJSVFjh071uFxBw4ckEceeUQmTZrks2sFAAAGDieff/653HHHHTJhwgQ5fPiw3rZmzRrZtGmTW+fJycmRWbNmycyZMyUuLk7y8/OlR48esmrVqnaPaWxslOnTp8uSJUtk+PDhnt4CAACwSjh58803dXWje/fu8uWXX0p9fb3efvLkSXn66addPk9DQ4OUlpZKcnLyfy4oOFivFxcXt3vc0qVLpX///nLnnXee82eoa6upqXFaAACAxcLJU089pSscK1eulG7dujm2X3nllbppxlXHjx/XVZABAwY4bVfrlZWVbR6jKjOvvPKK/tmuyM7Olt69ezuWwYMHu3x9AADAJOGkrKxMrr766lbb1Yf/iRMnxFtOnTolM2bM0MEkKirKpWMWLlyoKzr2paKiwmvXBwAA/PRsnejoaNm3b58eYdOyquFOHxAVMEJCQqSqqsppu1pXP6Ol/fv3646wN954o2NbU1OT/hoaGqpD04UXXuh0THh4uF4AAICFKyeqA+vcuXNl69ateijvkSNH5E9/+pMePXPvvfe6fJ6wsDCJj4+XoqIip7Ch1lVH25ZGjhwpO3bskO3btzuWqVOnyrXXXqu/p8kGAIAArZwsWLBAh4jrrrtOTp8+rZt4VHVChZM5c+a4dS41jDg9PV0SEhIkMTFRcnNzpa6uTo/eUdLS0iQmJkb3HVHzoIwaNcrp+D59+uivLbcDAIAACieqWvL444/LvHnzdPNObW2tHgbcq1cvt8+Vmpoq1dXVkpmZqTvBjh07VgoLCx2dZMvLy/UIHgAAEBg8CifNm2VUKOms2bNn66UtGzdu7PDY1atXd/rnA7Cm7y8J8/clwOQa65tE/tPzAGYIJwACAx/yAHyJcAIYGKEAQCAinAAuIigAgG8QTmAZhAcAsAbCCXyOEAEA6AjhBC4jVAAAfIFwEmAIGAAAoyOcmBQhAwBgVYQTAyBoAADwH4QTLyBsAADgOcKJGwgdAAB4X8CHEwIHAADGErDh5MSIMAkJJ5gAAGA0wf6+AAAAgOYIJwAAwFAIJwAABLC8vDyJjY2ViIgISUpKkpKSknb3XblypUyaNEn69u2rl+Tk5A73v+eeeyQoKEhyc3PduibCCQAAAaqgoEAyMjIkKytLtm3bJmPGjJGUlBQ5duxYm/tv3LhRbr/9dvnkk0+kuLhYBg8eLDfccIMcPny41b7r16+XLVu2yKBBg9y+LsIJAAABKicnR2bNmiUzZ86UuLg4yc/Plx49esiqVava3P9Pf/qT3HfffTJ27FgZOXKk/PGPf5SmpiYpKipy2k+FlTlz5uj9u3Xr5vZ1EU4AALCQmpoap6W+vr7N/RoaGqS0tFQ3zdgFBwfrdVUVccXp06fl7Nmz0q9fP8c2FVZmzJgh8+bNk8suu8yjewjYocQAABhJr38GS0i45zWDxvp/H6uaWppTTTaLFy9utf/x48elsbFRBgwY4LRdre/Zs8elnzl//nzdbNM84CxbtkxCQ0PlgQce8PBOCCcAAFhKRUWFREZGOtbDw8O98nOeeeYZWbdune6HojrTKqoS8/zzz+v+K6ojrKdo1gEAwEIiIyOdlvbCSVRUlISEhEhVVZXTdrUeHR3d4c947rnndDj58MMPZfTo0Y7tn3/+ue5MO2TIEF09UcvBgwfl4Ycf1iOCXEU4AQAgAIWFhUl8fLxTZ1Z759YJEya0e9zy5cvlySeflMLCQklISHB6TfU1+eqrr2T79u2ORTX7qP4nGzZscPnaaNYBACBAZWRkSHp6ug4ZiYmJej6Suro6PXpHSUtLk5iYGMnOznb0J8nMzJQ///nPuhJSWVmpt/fq1Usv559/vl6aU6N1VCXmkksucfm6CCcAAASo1NRUqa6u1oFDBQ01RFhVROydZMvLy/UIHrsVK1boUT633HKLS51uPUU4AQAggM2ePVsvbVGdXZs7cOCA2+f35Bj6nAAAAEMhnAAAAEMhnAAAAEMhnAAAAEMhnAAAAEMhnAAAAEMhnAAAAEMhnAAAAEMhnAAAAEMhnAAAAEMhnAAAAEPh2Tro0KkLm/x9CQHjvP38WwEAFMKJl/HhDqv/XSFUAehqARtOaoc3SXCEOT8MACMxQqgiIAHWErDhBIB1+CogEYIA3yCcAIAfQhBBB2gf4QQATBx0CDmwIsIJAJhYVzdpEXZgBIQTAIADTVcwAsIJAMArCDrwFOEEABAQQYeAYx6EEwBAQPAk4DSd8f88PoGIGAkAAAyFcAIAAAzFEOEkLy9PYmNjJSIiQpKSkqSkpKTdfd966y1JSEiQPn36SM+ePWXs2LGyZs0an14vAACwcDgpKCiQjIwMycrKkm3btsmYMWMkJSVFjh071ub+/fr1k8cff1yKi4vlq6++kpkzZ+plw4YNPr92AABgwXCSk5Mjs2bN0gEjLi5O8vPzpUePHrJq1ao297/mmmvk5ptvlksvvVQuvPBCmTt3rowePVo2bdrk82sHAAAWCycNDQ1SWloqycnJ/7mg4GC9rioj52Kz2aSoqEjKysrk6quv9vLVAgAAyw8lPn78uDQ2NsqAAQOctqv1PXv2tHvcyZMnJSYmRurr6yUkJET+8Ic/yPXXX9/mvmoftdjV1NR04R0AAICuZsp5Ts477zzZvn271NbW6sqJ6rMyfPhw3eTTUnZ2tixZssQv1wkAAEwWTqKionTlo6qqymm7Wo+Ojm73ONX0M2LECP29Gq2ze/duHULaCicLFy7U4aV55WTw4MFdeh8AAMAifU7CwsIkPj5eVz/smpqa9PqECRNcPo86pnnTTXPh4eESGRnptAAAAOPye7OOqmqkp6fruUsSExMlNzdX6urq9OgdJS0tTfcvUZURRX1V+6qROiqQfPDBB3qekxUrVvj5TgAAgCXCSWpqqlRXV0tmZqZUVlbqZprCwkJHJ9ny8nLdjGOngst9990nhw4dku7du8vIkSNl7dq1+jwAAMD8gmxqPG4AUX1OevfuLUOXPSXBERH+vhwAgIE1nTkjB+c/oUeJeqtbgP1z6bK7n5aQcM8/lxrrz8jOlx7z6rUGzCRsAAAAzRFOAACAoRBOAACAoRBOAACAoRBOAACAoRBOAACAoRBOAACAoRBOAAAIYHl5eRIbGysRERGSlJQkJSUl7e67c+dO+eUvf6n3DwoK0rO6t+Xw4cNyxx13yPnnn68nTL388svliy++MM8Msf7Sc2iNhPRo+3k8ZlX7bW9/XwIAwEQKCgr0Y2Ty8/N1MFFhIyUlRcrKyqR///6t9j99+rQMHz5cfvWrX8lDDz3U5jm///57ufLKK+Xaa6+V//3f/5Wf/OQnsnfvXunbt6/L1xWw4cSKeg076dOfRxgCAHPLycmRWbNmOZ5np0LK+++/L6tWrZIFCxa02n/cuHF6Udp6XVm2bJkMHjxYXn31Vce2YcOGuXVdhBMYIgwRdACg66bDby48PFwvLTU0NEhpaaksXLjQsU09yy45OVmKi4vFU++++66uvqjqyqeffqof3queiadCkKsIJ7BM0CHgADCzPvsaJDTU866gP/7YoL+qqkVzWVlZsnjx4lb7Hz9+XBobGx0P2rVT63v27PH4Ov75z3/KihUrdHPRY489Jv/4xz/kgQcekLCwMElPT3fpHIQTWEZnAw7hBoAVVFRUOD34r62qiTc1NTVJQkKCPP3003r9iiuukK+//lo3GRFOADdRvQFgBZGRkS49lTgqKkpCQkKkqqrKabtaj46O9vjnDxw4UOLi4py2XXrppfLmm2+6fA7CCdCFqN4AMIuwsDCJj4+XoqIimTZtmqPqodZnz57t8XnVSB012qe5b775RoYOHeryOQgngEXCDcEGgLtUvxDV1KKaYRITE/VQ4rq6OsfonbS0NN2hNTs729GJdteuXY7v1Xwm27dvl169esmIESP0djXEeOLEibpZ59Zbb9Xzprz88st6cRXhBLAIgg0Ad6Wmpkp1dbVkZmZKZWWljB07VgoLCx2dZMvLy/UIHrsjR47oPiR2zz33nF4mT54sGzdu1NvUUOP169frUUBLly7Vw4hV6Jk+fbrL1xVks9lsEmBDrHr37i1x6x6VkB6+7SQEWAmBBoGg6cwZOTj/CTl58qRL/Tg687l05XWLJTQ0wuPz/PjjGfl70WKvXquvUDkB4PVKDUEGgDsIJwAMEWQIMADsCCcADIEAA8COcALANAgwQGAgnACwFAIMYH6EEwAB51wBhvAC+BfhBABaoPoC+BfhBAA8QPUF8B7CCQB4AeEF8BzhBAD8gPACtI9wAgAGRL8XBDLCCQCYFNUXWBXhBAAsiuoLzIpwAgABjOoLjIhwAgBoF9UX+APhBADQKVRf0NUIJwAAr6L6AncRTgAAfkf1Bc0RTgAAhkd4CSyEEwCApcMLwcV8Ajac/OyCvRLeq5t8WH6Jvy8FAOBFVF3MJ2DDid0NQ8pc2o8QAwCBF15qdof79FrwbwEfTro6xCgEGQCwhp5Da/x9CQGJcOIFVGMAAPAc4cTgIYYAAwAINIQTg6MKAwAINIQTi6AKAwCwCsJJACHAAADMgHACtwIM4QUA4G2EE7iF8AIA8DbCCXwWXgguAABXEE7gM1RdAACuIJzAMKi6AACUYCP8MeTl5UlsbKxERERIUlKSlJSUtLvvypUrZdKkSdK3b1+9JCcnd7g/rBNc2lsAANbi98pJQUGBZGRkSH5+vg4mubm5kpKSImVlZdK/f/9W+2/cuFFuv/12mThxog4zy5YtkxtuuEF27twpMTExfrkH+BcVFwCwliCbzWbz5wWoQDJu3Dh58cUX9XpTU5MMHjxY5syZIwsWLDjn8Y2NjbqCoo5PS0s75/41NTXSu3dvmb3pZgnv1a1L7gHmRHABcC6Np+tl123L5eTJkxIZGemVn2H/XLryusUSGhrh8Xl+/PGM/L1osVevNSAqJw0NDVJaWioLFy50bAsODtZNNcXFxS6d4/Tp03L27Fnp169fm6/X19frpflfAqCjiguhBQACOJwcP35cVz4GDBjgtF2t79mzx6VzzJ8/XwYNGqQDTVuys7NlyZIlXXK9CAw0EwFAgPc56YxnnnlG1q1bp/uhqP4nbVFVGdWnpXnlRDUbAZ6g2gIAFg8nUVFREhISIlVVVU7b1Xp0dHSHxz733HM6nHz00UcyevTodvcLDw/XC+BNhBYAsEg4CQsLk/j4eCkqKpJp06Y5OsSq9dmzZ7d73PLly+U3v/mNbNiwQRISEnx4xYB7CC0AYMJmHdXkkp6erkNGYmKiHkpcV1cnM2fO1K+rEThqiLDqO6KoocOZmZny5z//Wc+NUllZqbf36tVLL4AZEFoAwMCTsKWmpuomGhU4xo4dK9u3b5fCwkJHJ9ny8nI5evSoY/8VK1boUT633HKLDBw40LGocwBmx0RzAIw8EaryxhtvyMiRI/X+l19+uXzwwQdOr9fW1urWjwsuuEC6d+8ucXFxei4zU1VOFHUT7TXjqM6uzR04cMBHVwUYB5UWAEaYCHXz5s16IlTVmvGLX/xCt2Kobhnbtm2TUaNG6X3U+T7++GNZu3atDj0ffvih3HfffXpk7dSpU80xCZuvMQkbAgGhBegaVp+ELcnNiVBVa4fqevHee+85to0fP163fNirIyqkqP0WLVrk2Ef1L/35z38uTz31lDmadQB0PZqHgMBVU1PjtDSfiLStiVCbzxN2rolQ1faW84qpSkvz/dXjZd599105fPiwqPrHJ598It98841+1IypmnUA+AbNQ4Bxdd91VEKDwzw+/semBv215VxeWVlZsnjx4i6ZCFUNQmlrf/vgFOWFF16Qu+66S/c5CQ0N1YFHPbT36quvdvleCCcACC2AhVRUVDg16/h6ri8VTrZs2aKrJ0OHDpXPPvtM7r///g5nc2+JcALArdBCYAGMLTIy0qU+J55MhKq2d7T/Dz/8II899pisX79epkyZorepiVLVSFw1qtbVcEKfEwBuoT8LYA1hzSZCtbNPhDphwoQ2j1Hbm++v/O1vf3Psrx7EqxbVlNOcCkHq3K6icgKgS9A0BJhPhpsToc6dO1cmT54sv/3tb3VlRD3f7osvvpCXX35Zv64qNur1efPm6TlOVLPOp59+Kq+//rrk5OS4fF2EEwBeRdMQYFypqalSXV2tJ0JVnVrVkOCWE6E2r4KokThqbpMnnnhCN99cdNFF8vbbbzvmOFFUYFEP3Z0+fbp89913OqCoR87cc889Ll8X85wAMAQCCwJ9npPkgXd3erTOR0df8uq1+gqVEwCGQIUFgB3hBIBhEViAwEQ4AWDqwEJYAayHcALA1KiuANZDOAFgOVRXAHML2HByW5+t0uu8/wyPWvP9RL9eDwDvIawA5hKw4aSlGX03t7md0AJYD2EFMDbCiQehhcACWAthBTAWwkkXBBbCCmAthBXAvwgnXYDqCmBthBXAtwgnXkJgAayLsAJ4F+HEh2gOAqwfVggqQOcRTvyIsAJYD1UVoPMIJwZCUxBgPYQVwH2EE4OjugJYC01AwLkRTkyGsAJYB1UVoG2EE5MjrADWQVgB/o1wYjGEFcA6aAJCoCKcWBxhBbAGqioIJISTAA4rBBXAvKiqwMoIJwGMqgpgDVRVYDWEEzgQVgBroKoCsyOcoF00AQHmR1UFZkQ4gUuoqgDWQFiBGRBO4BGqKoA10AQEIyKcoNOoqgDWQFUFRkE4QZejqgJYA1UV+AvhBF5FUAGsgaoKfIlwAp+h+QewDqoq8CbCCfyGqgpgDVRV0NUIJzAEggpgHVRV0FmEExgOQQWwDoIKPEE4gaERVADroPkHriKcwDQIKoC1UFVBewgnMCVG/gDWQlBBc4QTWAJVFcA6CCognMByCCqAddBPJTARTmBpBBXAWqiqBAbCCQIGQQWwFoKKdQX7+wIAfwWVlp1qAZg7qLRsAoJ5UTlBQKOaAlgL1RRroHIC/D+qKYC1UE0xL7+Hk7y8PImNjZWIiAhJSkqSkpKSdvfduXOn/PKXv9T7BwUFSW5urk+vFYEVUggqgLVCCkGl85/DyhtvvCEjR47U+19++eXywQcfOL1us9kkMzNTBg4cKN27d5fk5GTZu3evmCacFBQUSEZGhmRlZcm2bdtkzJgxkpKSIseOHWtz/9OnT8vw4cPlmWeekejoaJ9fLwIPIQWwFoJK5z6HN2/eLLfffrvceeed8uWXX8q0adP08vXXXzv2Wb58ufz+97+X/Px82bp1q/Ts2VOf88yZM+KqIJuKOH6iEtq4cePkxRdf1OtNTU0yePBgmTNnjixYsKDDY1XKe/DBB/XijpqaGundu7ds+nqQ9DrP74UjmBB9UwBr6ahvSuPpetl123I5efKkREZGeuXn2z+XkgfeLaHBYR6f58emBvno6EtuXau7n8OpqalSV1cn7733nmPb+PHjZezYsTqMqEgxaNAgefjhh+WRRx7Rr6vrGTBggKxevVpuu+02Y3eIbWhokNLSUlm4cKFjW3BwsC7/FBcXd9nPqa+v14ud+kNS6mqbuuxnILDcHLpJf113IsnflwKgC0zu9+9/9X986KI2w4nii3/H/2hrEGnq5PH/H3aaCw8P10tXfA6r7arS0pyqirz99tv6+2+//VYqKyv1OexU8FIhSB1r+HBy/PhxaWxs1GmqObW+Z8+eLvs52dnZsmTJklbbU8ZXdtnPQKBa7+8LAOAj//rXv/SHrDeEhYXprgobK1/t9Ll69eqlKx/NqSabxYsXd8nnsAoebe2vtttft29rbx9XWH4osUqEzVPeiRMnZOjQoVJeXu61v2i+oJKx+gtYUVHhtVKjt1nhHhTuwziscA9WuQ8r3IO92j5kyBDp16+f136G6liqKg6qktFZNptNDxhprq2qidH5LZxERUVJSEiIVFVVOW1X613Z2bW9cpYKJmb+hbFT92D2+7DCPSjch3FY4R6sch9WuAd7c4c3qYCiFqN/DkdHR3e4v/2r2qZG6zTfR/VLcZXfeoSqMlZ8fLwUFRU5tqmOOGp9woQJ/rosAAACQpgHn8Nqe/P9lb/97W+O/YcNG6YDSvN9VBVNjdpx57Pdr806qrklPT1dEhISJDExUc9bonoBz5w5U7+elpYmMTExut+Iokpeu3btcnx/+PBh2b59u25jGzFihD9vBQAA08lw83N47ty5MnnyZPntb38rU6ZMkXXr1skXX3whL7/8sn5dNSmpUbRPPfWUXHTRRTqsLFq0SI/gUUOOXWbzsxdeeME2ZMgQW1hYmC0xMdG2ZcsWx2uTJ0+2paenO9a//fZb1V261aL2c9WZM2dsWVlZ+quZWeE+rHAPCvdhHFa4B6vchxXuwUr30VWfw8pf//pX28UXX6z3v+yyy2zvv/++0+tNTU22RYsW2QYMGGALDw+3XXfddbaysjKbO/w6zwkAAEBLzEIGAAAMhXACAAAMhXACAAAMhXACAAAMxdTh5LPPPpMbb7xRD1FSw5fsc/t39rHN7j4+2oj3oaYqVudqvqhHXPvzPt566y254YYb5Pzzz9evq2HgrjjX47mNfg/qYVct3wtvT7bU0X2cPXtW5s+fr/8s1dNC1T5quOCRI0cM9bvhjXsw4u+FuiZ1Deo++vbtq3+/1ZwQZvv/lCf34ev341z30Nw999yj91FDa432XgQCU4cTNRZbPd5Z/cVoiyePbXb38dFGvQ/lsssuk6NHjzqWTZv+/cA6f92Hev2qq66SZcuWuXxOVx7PbfR7UNQMmc3fi4MHD4o3dXQfp0+f1n+31dwD6qsKXGVlZTJ16tQOz+nr3w1v3IMRfy8uvvhi/UTYHTt26GtRH3IqAFdXV5vq/1Oe3Iev349z3YPd+vXrZcuWLTrEnIs/3ouAYLMIdSvr1693GmcdHR1te/bZZx3bTpw4ocdc/+Uvf2n3PGqM9/333+9Yb2xstA0aNMiWnZ1tM9N9qHH5Y8aMsflLy/tozj5fzZdffnnO89x66622KVOmOG1LSkqy3X333Taz3MOrr75q6927t81fOroPu5KSEr3fwYMHDfm70VX3YOTfC7uTJ0/q/T766CPT/H/K0/vw5/vR3j0cOnTIFhMTY/v6669tQ4cOtf3ud7/r8Dz+fi+sytSVk46c67HNbbE/Prr5Med6fLQR78NONf2o5D98+HCZPn26ftih2ah7bH7vivpXib/eD0/V1tbqB06qB6HddNNNsnPnTjHaw81UCbtPnz6m+d1w9x7M8Huh/pzVTJvqd1z9C9ys74Ur92HE90NN3T5jxgyZN2+eruicixneC7OybDjx5LHNHT0+2p1HPXclTx8/rcKL6utQWFgoK1as0CFn0qRJcurUKTGTcz2e2wwuueQSWbVqlbzzzjuydu1a/T/AiRMnyqFDh8QIVPOg6r+hms/ae0CbEX833L0HI/9evPfee/oxHKrPwu9+9zv9rBL1UDazvRfu3IcR3w/VXBsaGioPPPCAS/sb+b0wO78+Wwfe8/Of/9zx/ejRo/X/BNS/3P/617/q/hvwHfWwq+YPvFLB5NJLL5WXXnpJnnzySb9em+pYeuutt+pO1+rDwYzcuQej/l5ce+21unO1+rBbuXKlvh/Vv6x///5iJu7eh5HeD1UBef7553W/EVWBg39ZtnLS/LHNrj4K2pPHRxvxPtqiSt2qw9q+ffvETM71eG4z6tatm1xxxRV+fy/sH+qqc676F25HFQcj/m64ew9G/r1QndzVw0vHjx8vr7zyiv7Xu/pqpvfC3fsw2vvx+eef606sQ4YM0detFvX36uGHH9ade832XpidZcOJJ49t9uTx0d7WVY+fVn0e9u/fr4cjm8m5Hs9tRqoMrEY0+PO9sH+oq/b+jz76SA+N7ogRfzfcvQcz/V6oP9v6+nrTvBee3IfR3g/V1+Srr77SlR/7ovrCqP4nGzZsMP17YTo2Ezt16pQeLaEWdSs5OTn6e3tv/WeeecbWp08f2zvvvGP76quvbDfddJNt2LBhth9++MFxjp/97Gf6iYx269at0yNhVq9ebdu1a5ftrrvu0ueorKw01X08/PDDto0bN+pRJX//+99tycnJtqioKNuxY8f8dh//+te/9Lp6gqV6Xf1Zq/WjR486zjFjxgzbggULHOvq2kNDQ23PPfecbffu3bp3f7du3Ww7duwwzT0sWbLEtmHDBtv+/fttpaWltttuu80WERFh27lzp1fu4Vz30dDQYJs6dartggsusG3fvl1fu32pr683zO+GN+7BaL8XtbW1toULF9qKi4ttBw4csH3xxRe2mTNn6j9nNVrELP+f8vQ+fP1+nOv3u6W2RusY4b0IBKYOJ5988on+C9ZysT/e2ZXHNqu/fOoDz9XHR5vlPlJTU20DBw7U96CGxan1ffv2+fU+1JDatl5vft2ePJ7b6Pfw4IMPOv4+qffwv/7rv2zbtm3z2j2c6z7sw6DbWtRxRvnd8MY9GO33Qv0D4+abb9ZDT9U1qWtToUsNi27O3++Ft+7D1+/HuX6/XQknRngvAkGQ+o+/qzcAAACW73MCAADMiXACAAAMhXACAAAMhXACAAAMhXACAAAMhXACAAAMhXACAAAMhXACAAAMhXACAAAMhXACAAAMJdTfFwDAf6655hoZNWqU/n7NmjXSrVs3uffee2Xp0qUSFBTk78sDEKConAAB7rXXXpPQ0FApKSmR559/XnJycuSPf/yjvy8LQADjwX9AgFdOjh07Jjt37nRUShYsWCDvvvuu7Nq1y9+XByBAUTkBAtz48eOdmnAmTJgge/fulcbGRr9eF4DARTgBAACGQjgBAtzWrVud1rds2SIXXXSRhISE+O2aAAQ2wgkQ4MrLyyUjI0PKysrkL3/5i7zwwgsyd+5cf18WgADGUGIgwKWlpckPP/wgiYmJulqigsldd93l78sCEMAIJ0CAU3Ob5ObmyooVK/x9KQCg0awDAAAMhXACAAAMhUnYAACAoVA5AQAAhkI4AQAAhkI4AQAAhkI4AQAAhkI4AQAAhkI4AQAAhkI4AQAAhkI4AQAAhkI4AQAAYiT/BymkprrzK026AAAAAElFTkSuQmCC",
                        "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)]].reshape(len(p), len(e))),\n",
                ")\n",
                "plt.colorbar(cb)\n",
                "plt.xlabel(\"p\")\n",
                "plt.ylabel(\"e\")\n",
                "plt.show()"
            ]
        }
    ],
    "metadata": {
        "kernelspec": {
            "display_name": "few2.0rc1",
            "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.9"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 2
}