{ "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": 1, "metadata": {}, "outputs": [], "source": [ "from few.amplitude.ampinterp2d import AmpInterpKerrEqEcc\n", "\n", "kerr_amp_spline = AmpInterpKerrEqEcc()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1, 6993)" ] }, "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", "\n", "amplitudes = kerr_amp_spline(a, p, e, xI)\n", "amplitudes.shape # has shape (N_points, N_modes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The data grids are not rectilinear in $(a, p, e)$ space, but may be easily mapped to from uniform distributions with provided utility functions:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total trajectory points: 10201\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAF4CAYAAAC/wIoGAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVkVJREFUeJzt3X9cVHW+P/AXP2SwFH+E/JQkLX+tCi5cWTRXLYxrXcrbt43VVojUbgZ7Xee2X6UUNE0sjeh2MTaVbCvTzXWtu7pYUWy5Ut5Qvte6ieEvyATltgriyihzvn+4jAzMwJyZ8/u8no/HPB7O8Zw5nzMz53Ne8+ZzzvETBEEAEREREZFB+KvdACIiIiIiKTHgEhEREZGhMOASERERkaEw4BIRERGRoTDgEhEREZGhMOASERERkaEw4BIRERGRoTDgEhEREZGhMOASERERkaEw4BIRERGRoTDgEhEREZEsPv30U6SlpSEqKgp+fn7YvXt3r8u0tbXhmWeewbBhw2CxWBAbG4vS0lJR6w30sr1ERERERD1qbW1FXFwcHnvsMTz44IMeLfPwww+jsbERW7Zswe23346zZ8/CbreLWi8DLhERERHJYtasWZg1a5bH85eVleHPf/4zTpw4gcGDBwMAYmNjRa/XdAHXbrfj+++/R//+/eHn56d2c4jIgARBQEtLC6KiouDvb7yRYOxHidTnSz9z5coV2Gw2r9fbdb+3WCywWCxevV5X77//PhITE/HCCy/gzTffxM0334z7778fq1evRt++fT1+HdMF3O+//x4xMTFqN4OITKC+vh5Dhw5VuxmSYz9KpB1i+5krV67g1ltvxvnz4v7k36Ffv364dOmS07T8/HysXLnSq9fr6sSJE9i/fz+Cg4Pxhz/8AU1NTXjyySfxv//7v3j99dc9fh3TBdz+/fsDACb/ZCkCA51/bQQd+96j17CNjPJ4fZ6+plTrk4IUbSYys2t2Gyqafuvob4ymY7vq6+sREhKicmuIjOf/3PFvvc7jbT9js9lw/rwdFV+EoV8/cX+BuXRJwPSkc932famqt8D1vxD5+fnh7bffxoABAwAAhYWFeOihh7Bx40aPq7imC7gdZfXAQAsCA4Od/i/QP8ij17B3Wc6doKPfAR6+pju20UMV/ZCkaDMRXWfUP993bFdISAgDLpEEZkVmOz33NI8A3vcz/fr5oV9/sUOorld95dz3IyMjER0d7Qi3ADBmzBgIgoDvvvsOd9xxh0evY7qA6yvbaOX+3KjkuoiIiEgZXQMt3TBlyhS8++67uHTpEvr16wcAOHbsGPz9/UUNxWDAFUFM4Aw6+p2MLZGHHttMRESkdWYOtJcuXUJtba3j+cmTJ1FdXY3Bgwfj1ltvRW5uLs6cOYPf/va3AIC5c+di9erVyMrKwqpVq9DU1IRf//rXeOyxx3iSmRGwektERKRfZg61nX355ZeYMWOG47nVagUAZGZmYuvWrTh79izq6uoc/9+vXz98+OGH+OUvf4nExETccsstePjhh7FmzRpR62XA9ZCS1Vs1wi2rt0RERN5joHVt+vTpEATB7f9v3bq127TRo0fjww8/9Gm9DLh/x4BHREREnmKg1TYGXI1h9ZaIiEh7GGj1hQHXA0Y/uYyIiIi6Y6jVLwZcDWH1loiISD16CrS///ZFDBiwWe1maBYDLhEREZmWnkIteY4BtxdKDU9g9ZaIiEh+DLTmwICrAQy3RERE8mGoNR8GXIkwMBIREWkDAy0x4ELdcMrqrfo8/Qz4vhERaRdDLXXGgNsD3i7XOKT4LN29BoMvEZE6GGrJHQZcCXgbcFi9lY+S723ndZnl/SUiUgtDLXmCAZcMQwsVd4ZdIiJpMdCSN0wfcNUKIazeSkMLodadjrYZ8X0nIpITQy35yvQB1x2eeKRtWg62XTHoEhH1jqGWpGTqgMvqrf7oKdh2ZRs91DCfAxGRFBhqvfOns8Vobm5WuxmaZtqAG3Tse8A/SO1mkAf0HGq7YjWXiMyOoZaUYNqAqxZWb8UxUrjtjNVcIjIThlpSGgOuCxx/qz6jBtvOGHKJyMgYaklNDLgKYvW2d2YItp0x5BKRkTDUklYw4BqYnoKT2YJtZwy5RKRnDLWkRQy4XhIbSMwc4HrD94Yhl4j0h8GWtIwB16D0EJYYbJ0x5BKR1jHUkl4w4HYhR+hikOuO74lrDLlEpDUMtaRH/mo3QI+0HkC03D7b6KEMt0QyKC4uRmxsLIKDg5GUlISDBw/2OH9RURFGjRqFvn37IiYmBkuWLMGVK1cUai1p3azIbMeDSI9YwZUZw9wNfC88wyouibVjxw5YrVaUlJQgKSkJRUVFSE1NRU1NDcLCwrrNv23bNixbtgylpaWYPHkyjh07hkcffRR+fn4oLCxUYQtIKxhoySgYcA1Gq8GI4VYchlwSo7CwEAsXLkRWVhYAoKSkBHv27EFpaSmWLVvWbf4DBw5gypQpmDt3LgAgNjYWc+bMwRdffKFou0kbGGrJiDhEoROGMOlxSAKRvGw2G6qqqpCSkuKY5u/vj5SUFFRWVrpcZvLkyaiqqnIMYzhx4gT27t2Le++91+X8bW1taG5udnqQ/nEIgj796Wyx2k3QBVZwRRJTVVM62Gmt4sdgSyS/pqYmtLe3Izw83Gl6eHg4jh496nKZuXPnoqmpCXfeeScEQcC1a9fwxBNP4Omnn3Y5f0FBAVatWiV520l5DLRkFqzgkiwYbn3H95DkUlFRgbVr12Ljxo04dOgQdu3ahT179mD16tUu58/NzcXFixcdj/r6eoVbTL7gCWOkpk8//RRpaWmIioqCn58fdu/e7fGyf/nLXxAYGIj4+HjR62UF9++kDhNmrt4ymBEpJzQ0FAEBAWhsbHSa3tjYiIiICJfLrFixAvPmzcOCBQsAAOPHj0draysef/xxPPPMM/D3d659WCwWWCwWeTaAZMNAS1rQ2tqKuLg4PPbYY3jwwQc9Xu7ChQvIyMjA3Xff3a1/8wQDrghaCpFaxGBLpLygoCAkJCSgvLwcs2fPBgDY7XaUl5cjJyfH5TKXL1/uFmIDAgIAAIIgyNpekhdDLSmh6zj8nn4Ez5o1C7NmzRK9jieeeAJz585FQECAqKpvBwZcsHorBYZbefBqCuQJq9WKzMxMJCYmYtKkSSgqKkJra6vjqgoZGRmIjo5GQUEBACAtLQ2FhYWYOHEikpKSUFtbixUrViAtLc0RdElfGGxJrN0tcQgW+oha5sqlqwA+QExMjNP0/Px8rFy5UrK2vf766zhx4gTeeustrFmzxqvXYMAlnzHcEqkrPT0d58+fR15eHhoaGhAfH4+ysjLHiWd1dXVOFdvly5fDz88Py5cvx5kzZzBkyBCkpaXhueeeU2sTyEsMtqSG+vp6hISEOJ5LOYTp22+/xbJly/DZZ58hMND7mGr6gOtpONNqFU3tdjHcEmlDTk6O2yEJFRUVTs8DAwORn5+P/Px8BVpGUmOoJbWFhIQ4BVyptLe3Y+7cuVi1ahVGjhzp02uZPuBKzUyBz0zbSkSkNgZbMrqWlhZ8+eWXOHz4sOMHu91uhyAICAwMxAcffIC77rrLo9cydcDVe0BTs3qr9/eOiEgvGGzJLEJCQnDkyBGnaRs3bsTHH3+MnTt34rbbbvP4tUwdcD3laZA0S+gzy3YSEamFoZaM4tKlS6itrXU8P3nyJKqrqzF48GDceuutyM3NxZkzZ/Db3/4W/v7+GDdunNPyYWFhCA4O7ja9Nwy4OqVW9ZbhlohIPgy21BM93qb3yy+/xIwZMxzPrVYrACAzMxNbt27F2bNnUVdXJ/l6GXB7ofZJXK4w3BIRGQuDLRnV9OnTe7y+9tatW3tcfuXKlV5dgsy0Adc2MkrSjTd6+DP69hERqYHBlkgepg24eqVG9ZbhlohIWgy2RPJiwO0BTy4z9rYRESmJoZZIOQy4OqJ09ZbhlojIdwy2RMpjwCWXGG61QYsnORKRZxhsidTDgOuG1oYnKBl0GG6JiLzHYEukPgZccsJwS0TkHQZbIu1gwNUBpaq3DLdEROIx2BJpDwOuD4wUCI20LUbB8bdE2sZgS6RdDLguaClYaKktRETEYEvq0ONtetXEgEus3moQf9gQaQ+DLZF+MOB6SYlQqETIYbglIuoZgy2R/vir3YDi4mLExsYiODgYSUlJOHjwYI/zFxUVYdSoUejbty9iYmKwZMkSXLlyRaHWGgvDLRFRzxhuifRJ1Qrujh07YLVaUVJSgqSkJBQVFSE1NRU1NTUICwvrNv+2bduwbNkylJaWYvLkyTh27BgeffRR+Pn5obCwULF2GyEYGmEbjIrDE4jUx2BLpG+qBtzCwkIsXLgQWVlZAICSkhLs2bMHpaWlWLZsWbf5Dxw4gClTpmDu3LkAgNjYWMyZMwdffPGF23W0tbWhra3N8by5ubnHNmklXMjZDoZb7dLK94/IrBhsiYxBtSEKNpsNVVVVSElJudEYf3+kpKSgsrLS5TKTJ09GVVWVYxjDiRMnsHfvXtx7771u11NQUIABAwY4HjExMdJuCBER6d6syGyGWyIDUa2C29TUhPb2doSHhztNDw8Px9GjR10uM3fuXDQ1NeHOO++EIAi4du0annjiCTz99NNu15Obmwur1ep43tzc7FPI1fvJZazeahert0TKY6glMibVTzITo6KiAmvXrsXGjRtx6NAh7Nq1C3v27MHq1avdLmOxWBASEuL0MCuGWyKiGxhuiYxLtYAbGhqKgIAANDY2Ok1vbGxERESEy2VWrFiBefPmYcGCBRg/fjz++Z//GWvXrkVBQQHsdrsSzZadXFU8hlttY/WWfCXmijTTp0+Hn59ft8d9992nYIvVw+EIRManWsANCgpCQkICysvLHdPsdjvKy8uRnJzscpnLly/D39+5yQEBAQAAQRB8b1MvIUOvIVGv7TYLhlvyVccVafLz83Ho0CHExcUhNTUV586dczn/rl27cPbsWcfjq6++QkBAAH72s58p3HJlMdgSmYeqV1GwWq3IzMxEYmIiJk2ahKKiIrS2tjquqpCRkYHo6GgUFBQAANLS0lBYWIiJEyciKSkJtbW1WLFiBdLS0hxBV88YdIjIG2KvSDN48GCn59u3b8dNN91k2IDLUEtkPqoG3PT0dJw/fx55eXloaGhAfHw8ysrKHCee1dXVOVVsly9fDj8/PyxfvhxnzpzBkCFDkJaWhueee06tTdA8Vm+1jT9qyFcdV6TJzc11TOvtijRdbdmyBT//+c9x8803u/x/sZdb1BKGWzKCP50tVrsJuqP6rXpzcnKQk5Pj8v8qKiqcngcGBiI/Px/5+fkKtMyZ3EFRjqDDcKttDLckBW+uSNPZwYMH8dVXX2HLli1u5ykoKMCqVat8bquSGGyJzE1XV1EgIiJpbdmyBePHj8ekSZPczpObm4uLFy86HvX19Qq2UByOsyUiQAMVXK1Qs5rG6q35sHpLUvHmijQdWltbsX37djz77LM9zmexWGCxWHxuq9zkDrau+lXuy0TaxIDrAb2FRb2112x4QCQpdb4izezZswHcuCKNu+FfHd599120tbXhF7/4hQItlY9cwdaTvrTrPNy/ibSBAVdlUneGDLfaxoMfyUHsFWk6bNmyBbNnz8Ytt9yiRrMlIXW49bUP7Vie+zqRujgGl0ghPOCRXNLT07Fhwwbk5eUhPj4e1dXV3a5Ic/bsWadlampqsH//fsyfP1+NJvtMjrG2UhYIWGwguu7TTz9FWloaoqKi4Ofnh927d/c4/65duzBz5kwMGTIEISEhSE5Oxr59+0SvlxXcXuipk9JTW82G4ZbkJuaKNAAwatQoSW6QowYtB9uur8t9n8yutbUVcXFxeOyxx/Dggw/2Ov+nn36KmTNnYu3atRg4cCBef/11pKWl4YsvvsDEiRM9Xi8DLtQLH1Kul+FWu3iAI5KGXoJt13WwDyAzmzVrFmbNmuXx/EVFRU7P165di/feew//+Z//yYBL8ro4wrOzqQccb+t9JiIiD+gx3BIZVdebvch5pRW73Y6WlpZud2DsDQNuD+TsAPVYvfU02Hae3+whl5UbIt/cPeP6iXFBEr6m0uGWVVzSok8aRyLwkrjj+rXWNgAfICYmxml6fn4+Vq5cKV3jOtmwYQMuXbqEhx9+WNRyDLg6p0RHLTbYdl3WrCGXBzQi33SEWymxckt6o8Xb9NbX1yMkJMTxXK7q7bZt27Bq1Sq89957CAsLE7Ws6QOuGiFEL8HHl2Db9XXMFnL18hkTaVHXYCvV/sRwSySNkJAQp4Arh+3bt2PBggV49913kZKSInp50wdcPZOrs5Yq2JoVwy2R9+So2gIMt0R68s477+Cxxx7D9u3bcd9993n1GqYOuD0FEbk6Q61XIhhufcNwS+Q9I4dbjsMls7p06RJqa2sdz0+ePInq6moMHjwYt956K3Jzc3HmzBn89re/BXB9WEJmZiZefvllJCUloaGhAQDQt29fDBgwwOP18kYP5CBnuDVDcObBi8g7syKzDR1uiczsyy+/xMSJEx2X+LJarZg4cSLy8vIAAGfPnkVdXZ1j/tdeew3Xrl1DdnY2IiMjHY/FixeLWq9pK7hBx74H/KU8L1c5UnfYZgifcmO4JfLO3TMKerxCgi/7FsMtkfqmT5/e401ltm7d6vTc1Y1pvMEKrgtaHp7AcKs9DLdE3pGragsw3BKZHQOuiTHc+o7hlkh7GG6JiAFXIVqr3jLc+o7hlog8xf6CSFkMuDqh93BrtOvg8mBFJB13/Zs3+xmrt0QEMOB2I0fnqKUwxMqt77T0eRLpFe9SRkRyYsDVAak6bYZb3zHcEklL69cGJyJ9Mu1lwsyG4dY3DLZERKS0P50tVrsJusUKrsx8DUZSVCUYbn3DcEukPDH7ndart+xDiJTHgNuJ1jpJo4RbPZ9gxgMTkbZprd8mIm1gwJWR2uFIC+FWz9T+/IjMgAGViOTAgKtRvnb6DLe+Ybgl0j49hGP2JUTqYMA1IC2FWz0OT+ABiYiISN8YcP9O6kqALyHJl7ZoKdzqEcMtkfo82Q/1UL0lIvUw4BoIw61vGG6JlGX0fc7o20ekZQy4MlCjeqvFcKun4Qk8EBHpB6u3RNQbBlwNYaetDoZbMoLi4mLExsYiODgYSUlJOHjwYI/zX7hwAdnZ2YiMjITFYsHIkSOxd+9ehVp7nZH7PPYrROrincyg/06W1Vvv8ABERrFjxw5YrVaUlJQgKSkJRUVFSE1NRU1NDcLCwrrNb7PZMHPmTISFhWHnzp2Ijo7G6dOnMXDgQEXae/eMArf/19t+qff+moiUwYArMW9Dk5GGJugBwy0ZSWFhIRYuXIisrCwAQElJCfbs2YPS0lIsW7as2/ylpaX44YcfcODAAfTp0wcAEBsb6/b129ra0NZ240drc3OztBtgMOxfSAq8Ta9vOERBxxhuvcODDxmJzWZDVVUVUlJSHNP8/f2RkpKCyspKl8u8//77SE5ORnZ2NsLDwzFu3DisXbsW7e3tLucvKCjAgAEDHI+YmBhZtqU3rN4SkadMH3Cl7DCVrN5qOdxqeXgCwy0ZTVNTE9rb2xEeHu40PTw8HA0NDS6XOXHiBHbu3In29nbs3bsXK1aswIsvvog1a9a4nD83NxcXL150POrr6yXfDqNgH0OkDRyiQKbBAw/RdXa7HWFhYXjttdcQEBCAhIQEnDlzBuvXr0d+fn63+S0WCywWdX9Us3pLRGIw4KqM1VtlMNySUYWGhiIgIACNjY1O0xsbGxEREeFymcjISPTp0wcBAQGOaWPGjEFDQwNsNhuCgoJkbbM7et9P9d5+IiMx9RAFLQxPEEvL4VareNAhIwsKCkJCQgLKy8sd0+x2O8rLy5GcnOxymSlTpqC2thZ2u90x7dixY4iMjFQt3PaE1VsiEsvUAVdtRuu0tVi9ZbglM7Bardi0aRPeeOMNfPPNN1i0aBFaW1sdV1XIyMhAbm6uY/5Fixbhhx9+wOLFi3Hs2DHs2bMHa9euRXZ2tlqboHvsa4i0xbRDFGwjo3S38azeisMDDplFeno6zp8/j7y8PDQ0NCA+Ph5lZWWOE8/q6urg73+jnhETE4N9+/ZhyZIlmDBhAqKjo7F48WIsXbpUrU1wy2iFACJSht4yniZ5E6TEdtpaD7daqt4y2JIZ5eTkICcnx+X/VVRUdJuWnJyMzz//XOZWeU7P+62e205kVByioAMMt57jgYaIlMQ+h0ibGHBVwD+5yYMHGiLtmxXp+Thf9pVE+vfpp58iLS0NUVFR8PPzw+7du3tdpqKiAj/+8Y9hsVhw++23Y+vWraLXy4DrI7lDFau3nmG4JdIHI4VW9jtEvWttbUVcXByKiz279fDJkydx3333YcaMGaiursavfvUrLFiwAPv27RO1Xo7BVZiYzl3r4VYreJAhMh4jBWEisf501rMwqAezZs3CrFmzPJ6/pKQEt912G1588UUA16/RvX//frz00ktITU31+HVYwSWvaaF6y3BLpH963I/12GYiqTQ3Nzs92tqkywOVlZVISUlxmpaamorKykpRr8MKrg/EdnCs3kqLBxgiUgP7HjKCuu9D4d83WNQy9r9dAXD9UoOd5efnY+XKlZK0q6GhwXGJww7h4eFobm7G3/72N/Tt29ej12HAJa+oXb3lAYbIuDg8gUjb6uvrERIS4nhusWivKMeA6yUzV28ZbonIrNj/EAEhISFOAVdKERERaGxsdJrW2NiIkJAQj6u3AMfgao7Ww63aeHAhIrWw/yGSX3JyMsrLy52mffjhh0hOThb1OqoH3OLiYsTGxiI4OBhJSUk4ePBgj/NfuHAB2dnZiIyMhMViwciRI7F3716FWusdI/25Tc3qLQ8uRMbTdb82Un9JRMClS5dQXV2N6upqANcvA1ZdXY26ujoAQG5uLjIyMhzzP/HEEzhx4gT+7//9vzh69Cg2btyI3/3ud1iyZImo9ao6RGHHjh2wWq0oKSlBUlISioqKkJqaipqaGoSFhXWb32azYebMmQgLC8POnTsRHR2N06dPY+DAgYq2W66gxeqtewy3RKQm9kFE3vnyyy8xY8YMx3Or1QoAyMzMxNatW3H27FlH2AWA2267DXv27MGSJUvw8ssvY+jQodi8ebOoS4QBKgfcwsJCLFy4EFlZWQCuX/tsz549KC0txbJly7rNX1paih9++AEHDhxAnz59AACxsbFKNlk0T6sRegi3alRveVAhIrWxHyLy3vTp0yEIgtv/d3WXsunTp+Pw4cM+rVe1IQo2mw1VVVVO1zrz9/dHSkqK22udvf/++0hOTkZ2djbCw8Mxbtw4rF27Fu3t7W7X09bW1u16bSQewy0RyY3DE4hIKqoF3KamJrS3t7u81llDQ4PLZU6cOIGdO3eivb0de/fuxYoVK/Diiy9izZo1btdTUFCAAQMGOB5dr90mlpjQZaTqrdIYbomM5+4ZBWo3QTT2RUT6pPpJZmLY7XaEhYXhtddeQ0JCAtLT0/HMM8+gpKTE7TK5ubm4ePGi41FfX69gi41B6eotDyhEpAXsi4j0S7UxuKGhoQgICHB5rbOIiAiXy0RGRqJPnz4ICAhwTBszZgwaGhpgs9kQFBTUbRmLxaLKBYhZvfUODyhE5tF5f+fwBKLr/nS2WO0mGIJqFdygoCAkJCQ4XevMbrejvLzc7bXOpkyZgtraWtjtdse0Y8eOITIy0mW4lbzNEocvPYRbJau3DLdEpBXsj4j0TdUhClarFZs2bcIbb7yBb775BosWLUJra6vjqgoZGRnIzc11zL9o0SL88MMPWLx4MY4dO4Y9e/Zg7dq1yM7OVmsTDI3hlojMiP0Rkf6pepmw9PR0nD9/Hnl5eWhoaEB8fDzKysocJ57V1dXB3/9GBo+JicG+ffuwZMkSTJgwAdHR0Vi8eDGWLl2q1ia45Mmf2vRQvVUKDyZEREQkJVUDLgDk5OQgJyfH5f9VVFR0m5acnIzPP/9c5lZ1Z7YQplT11mzvKxF1p6Xxt+yTiIxBV1dR0AMjVG8ZbonIjNgnERkHAy6pggcSInNjH0BEclJ9iILaLo6w9Fqx9LQjZvXWMzywEVGHnvrNrv2lnP0T+yUiY2EFlxTFgwiRPIqLixEbG4vg4GAkJSXh4MGDbufdunUr/Pz8nB7BwcEKtrZ3ShYD2C8RGQ8DroLMXr3lQYRIHjt27IDVakV+fj4OHTqEuLg4pKam4ty5c26XCQkJwdmzZx2P06dPK9hi7WC/RGRMDLi9kHJ4gpYx3BLpV2FhIRYuXIisrCyMHTsWJSUluOmmm1BaWup2GT8/P0RERDgeHZdnlMvdMwo8nlfrxQAi0j7TB1ylrhhg5g6b4ZZIPjabDVVVVUhJSXFM8/f3R0pKCiorK90ud+nSJQwbNgwxMTF44IEH8PXXX7udt62tDc3NzU4PSdruojDAoQlEJAXTB1wpsHrrHg8gRPJqampCe3t7twpseHg4GhoaXC4zatQolJaW4r333sNbb70Fu92OyZMn47vvXO+vBQUFGDBggOMRExPjU5u10C9ooQ1EXf3pbLHaTTAMBtweSNUBarl6y3BLZD7JycnIyMhAfHw8pk2bhl27dmHIkCH4zW9+43L+3NxcXLx40fGor6+XpV1a7iuJSF9Mf5kwkgfDLZEyQkNDERAQgMbGRqfpjY2NiIiI8Og1+vTpg4kTJ6K2ttbl/1ssFlgsxgmf7J+IjI8VXB/1NjxByxUJuaq3PHgQKScoKAgJCQkoLy93TLPb7SgvL0dycrJHr9He3o4jR44gMjJSrmZqBvsnInNgBdcNo3eCDLdExmG1WpGZmYnExERMmjQJRUVFaG1tRVZWFgAgIyMD0dHRKCi4fiWDZ599Fj/5yU9w++2348KFC1i/fj1Onz6NBQsWKNbmrsUBLRcDiEh/GHB9oOfqrRwYbonUkZ6ejvPnzyMvLw8NDQ2Ij49HWVmZ48Szuro6+Pvf+IPdX//6VyxcuBANDQ0YNGgQEhIScODAAYwdO1atTeiVFD/K2UcRmQcDrgnJUb3lgYNIXTk5OcjJyXH5fxUVFU7PX3rpJbz00ksKtMozShQD2EcRmQvH4MpEq9VbhlsiUhP7CyJSAgOuC550wHq/9q1UeLAiIrGU7j/ZTxGZDwOuDMxUvSUi8oVW+0si0jcGXC+wensdqyJEpARffpyznyIyJ55kJjGtViOkrt7yoKEcuX5Q8TMkpcyKzAZUKAzwO05kXgy4XRixQ2S41Sa1/xIgZv38zMkX7r5rWi0IEKnhT2eL1W6CoXg1RCEjIwOvv/46jh8/LnV7NK+nUGCGzppBRxzb6KFuH3rS03bobVtIXUp9X9hXEWlHcXExYmNjERwcjKSkJBw8eLDH+YuKijBq1Cj07dsXMTExWLJkCa5cuSJqnV5VcIOCglBQUID58+cjOjoa06ZNw/Tp0zFt2jTccccd3rwkyUTK6i0PGD0zc9Bzt+38zlBnQUe/E72feNOH8XtHpB07duyA1WpFSUkJkpKSUFRUhNTUVNTU1CAsLKzb/Nu2bcOyZctQWlqKyZMn49ixY3j00Ufh5+eHwsJCj9frVQV38+bNOHbsGOrr6/HCCy+gX79+ePHFFzF69GgMHWrOg7wWq7cMt/JhFdMzfJ969tlnn+EXv/gFkpOTcebMGQDAm2++if3796vcMvlpsc8kIukVFhZi4cKFyMrKwtixY1FSUoKbbroJpaWlLuc/cOAApkyZgrlz5yI2Nhb33HMP5syZ02vVtyufrqIwaNAg3HLLLRg0aBAGDhyIwMBADBkyxJeXVFVvIc6sB2eG2+5BjbzH9/K63//+90hNTUXfvn1x+PBhtLVd/0F68eJFrF27VuXW6Rf7KyL5NTc3Oz06+q+ubDYbqqqqkJKS4pjm7++PlJQUVFZWulxm8uTJqKqqcgTaEydOYO/evbj33ntFtdGrIQpPP/00KioqcPjwYYwZMwbTpk3DsmXL8NOf/hSDBg3y5iVJYlJVb816sDBz8FJa1/faLN+5NWvWoKSkBBkZGdi+fbtj+pQpU7BmzRoVW0ZEZhBUH4SA4CBRy7RfsQMAYmJinKbn5+dj5cqV3eZvampCe3s7wsPDnaaHh4fj6NGjLtcxd+5cNDU14c4774QgCLh27RqeeOIJPP3006La6lXAXbduHYYMGYL8/Hw8+OCDGDlypDcvYxha+1Mbw613GGq1wSyBt6amBj/96U+7TR8wYAAuXLigfIMUIPf4W6N+V4i0pr6+HiEhIY7nFot0OaiiogJr167Fxo0bkZSUhNraWixevBirV6/GihUrPH4drwLu4cOH8ec//xkVFRV48cUXERQU5DjRbPr06boMvByeYE78XLXPqIE3IiICtbW1iI2NdZq+f/9+DB8+XJ1GKUSOooBRvhdEehASEuIUcN0JDQ1FQEAAGhsbnaY3NjYiIiLC5TIrVqzAvHnzsGDBAgDA+PHj0draiscffxzPPPMM/P09G13r1RjcuLg4/Ou//it27dqF8+fPY+/evQgKCkJ2djbGjBnjzUvqFqu3+sPxn/pmlM9v4cKFWLx4Mb744gv4+fnh+++/x9tvv42nnnoKixYtUrt5REQ+CwoKQkJCAsrLyx3T7HY7ysvLkZyc7HKZy5cvdwuxAQEBAABBEDxet1cVXEEQcPjwYVRUVKCiogL79+9Hc3MzJkyYgGnTpnnzkiQBhtue6T0QUXedP1O9fW+XLVsGu92Ou+++G5cvX8ZPf/pTWCwWPPXUU/jlL3+pdvNUJ6Y/09tnT2QmVqsVmZmZSExMxKRJk1BUVITW1lZkZWUBuH5vhejoaBQUFAAA0tLSUFhYiIkTJzqGKKxYsQJpaWmOoOsJrwLu4MGDcenSJcTFxWHatGlYuHAhpk6dioEDB3rzcppnpmBktAOFmT47s9PbUAY/Pz8888wz+PWvf43a2lpcunQJY8eORb9+/dRumiy8uQYuEelfeno6zp8/j7y8PDQ0NCA+Ph5lZWWOE8/q6uqcKrbLly+Hn58fli9fjjNnzmDIkCFIS0vDc889J2q9XgXct956C1OnTvVo/IUeeHsg1NLwBCmqt1oPBJ7iQZQA/VR3g4KCMHbsWLWbIZu7ZxQ4PZe639TyZ0tE1+Xk5CAnJ8fl/1VUVDg9DwwMRH5+PvLz831ap1cB97777vNppXpilrBkhIOEWT4rEq/ju2GE77lecf8kIiV5FXDJWNVbvR/0eeAkT+mlqmt2nvZp/AzJKP50tljtJhgOA67OSXk7Xr1hsCVfMOzqGz8zIuqJ6QNuT52kuwClpeqtr/R4kGCwJalxCIMyjNR3EpG2mTrg6v1gZrahCQy2JDdWddXnSb/Gz4aIemPqgKtnZgq3DLakBlZ1yejU7Fu5X5HcTBtwg459D/gHiV7OCH9i00vHwmBLWqCXqm5xcTHWr1+PhoYGxMXF4ZVXXsGkSZN6XW779u2YM2cOHnjgAezevVv+hvpIy5+BWvTYV0rVZn4fyB3TBtzeaLnD8KV6q4fOQMvvPZmbVqu6O3bsgNVqRUlJCZKSklBUVITU1FTU1NQgLCzM7XKnTp3CU089halTpyrYWhKLfaJ7Yt8bre27JB8GXNIMduKkF1oLuoWFhVi4cKHj1pclJSXYs2cPSktLsWzZMpfLtLe345FHHsGqVavw2Wef4cKFC7K20ZO/fvX2410r77dc2AfKz9P32OjfNTNgwBVBC8MTjFi9ZadOeqWF4Qs2mw1VVVXIzc11TPP390dKSgoqKyvdLvfss88iLCwM8+fPx2effdbjOtra2tDWdqPvaW5u9r3hJsd+T9s8+Xy0ekyl6xhwXdBqx8NwS6RdalV1m5qa0N7e7rive4fw8HAcPXrU5TL79+/Hli1bUF1d7dE6CgoKsGrVKl+b6hOt9mGeYl9nPL19pnr/zuodA66HtFC99ZYWdzJ29mRUttFDce3aFeCc2i1xraWlBfPmzcOmTZsQGhrq0TK5ubmwWq2O583NzYiJiel1ubtnFHjcLqPdtIZ9HDEAq4sBVye87fy1uAOx4yeSTmhoKAICAtDY2Og0vbGxEREREd3mP378OE6dOoW0tDTHNLvdDgAIDAxETU0NRowY4bSMxWKBxaLej3wt9mOusG8jMfh9kRcDbhf8wsmH7y2R9IKCgpCQkIDy8nLMnj0bwPXAWl5ejpycnG7zjx49GkeOHHGatnz5crS0tODll1/2qDJLN7BfI1+Vf5Lb+0wkGgOuB9QenqD36i0PAETyslqtyMzMRGJiIiZNmoSioiK0trY6rqqQkZGB6OhoFBQUIDg4GOPGjXNafuDAgQDQbbpUfOlDtdKPdcY+jUj7GHA1juGWiHqTnp6O8+fPIy8vDw0NDYiPj0dZWZnjxLO6ujr4+/ur3Er39DD+ln0Zkb4w4BqQFsItDwZEysrJyXE5JAEAKioqelx269at0jdIAuzLiMhbDLiduOrI1ByeoIeqhitGOiCo8fnr9XMnMhIj9WNEZsSAazBqVzz0eFBQe4x1Vz21h+GXjMbdd1qNvkyP/RcRucaAq1HeBBk1w61eDgxaC7NiuWs/gy+Rb/TShxGRZzRx1kFxcTFiY2MRHByMpKQkHDx40KPltm/fDj8/P8elcXyhpeEJDLfSuDjC0u1hVGbZTtIfrX8fbaOHarYPIyLvqV7B3bFjB6xWK0pKSpCUlISioiKkpqaipqYGYWFhbpc7deoUnnrqKUydOlXB1lJnWjwoaP1gqpSu7wMrvKQnSvxg12L/RUTSUb2CW1hYiIULFyIrKwtjx45FSUkJbrrpJpSWlrpdpr29HY888ghWrVqF4cOHK9ha+emlequVgwMrl57h+0RapMYPL1ZsicxB1YBrs9lQVVWFlJQUxzR/f3+kpKSgsrLS7XLPPvsswsLCMH/+/F7X0dbWhubmZqdHt3ZoaHiCWGYMtwxqvuN7SFolZ5+mdt9FRMpRdYhCU1MT2tvbHRcj7xAeHo6jR4+6XGb//v3YsmULqqurPVpHQUEBVq1a5WtTFSG2mmGmcMsgJp/O7y2HMpBU7p5RoHYTHBhsicxH9SEKYrS0tGDevHnYtGkTQkNDPVomNzcXFy9edDzq6+tlbqV3tB4s1PqzHquMymJll9Qkx492hlsic1K1ghsaGoqAgAA0NjY6TW9sbERERES3+Y8fP45Tp04hLS3NMc1utwMAAgMDUVNTgxEjRjgtY7FYYLG4P1jrtfNTsnqr9HvEcKUNrOySnJT4Tum1fyfzKP8kV+0mGJaqATcoKAgJCQkoLy93XOrLbrejvLzc5S0nR48ejSNHjjhNW758OVpaWvDyyy8jJiZGknYpHbC0PDRByQMEg612MeySN9TapxlsiUj1IQpWqxWbNm3CG2+8gW+++QaLFi1Ca2srsrKyAAAZGRnIzb3+Cyc4OBjjxo1zegwcOBD9+/fHuHHjEBQUpOamKMKI4ZZ/EtcXfl4kB6n6NoZbIu0Re7+DCxcuIDs7G5GRkbBYLBg5ciT27t0rap2qXwc3PT0d58+fR15eHhoaGhAfH4+ysjLHiWd1dXXw95cnh2uhI9RiNUyJ94UBSf86PkMtfodJ2+T6zmihTyciZ2Lvd2Cz2TBz5kyEhYVh586diI6OxunTpzFw4EBR61U94AJATk6OyyEJAFBRUdHjslu3bpW0LVoOXka4+LmW31/yDocvkBYw3BIpp+slV3s636nz/Q4AoKSkBHv27EFpaSmWLVvWbf7S0lL88MMPOHDgAPr06QMAiI2NFd1G1YcoqMU2MkrtJogKA3oPt/yztjnwcyZv+Nq/MdwSide/XkD/0yIf9QIAICYmBgMGDHA8CgpcXxbQm/sdvP/++0hOTkZ2djbCw8Mxbtw4rF27Fu3t7aK2TxMVXDPSWqVLrgMEw445cfgCKYXhlkh59fX1CAkJcTx3V7315n4HJ06cwMcff4xHHnkEe/fuRW1tLZ588klcvXoV+fn5HreRAbcTrYYxuau3DLckFwZd6qrrd8Hb/o3Blkg9ISEhTgFXSna7HWFhYXjttdcQEBCAhIQEnDlzBuvXr2fA1TotDU2Q4yDBYEtdMeiaD/sBIhJ7vwMAiIyMRJ8+fRAQEOCYNmbMGDQ0NMBms3l8xSzTjsEl6cMtx19Sb/gdISmwekukD53vd9Ch434HycnJLpeZMmUKamtrHTfyAoBjx44hMjJS1OVgGXD/TqmDrhaqt1LfdpehhcTid4YA7/o4hlsifRFzvwMAWLRoEX744QcsXrwYx44dw549e7B27VpkZ2eLWi+HKGiUnOFWSgwp5AsOXTCmu2cUADL0DQy3RPoj9n4HMTEx2LdvH5YsWYIJEyYgOjoaixcvxtKlS0WtlwFXQZ4exPUQbhlsSUoMuubgy+erpXArZf/H7zyZgdj7HSQnJ+Pzzz/3aZ0MuFAmrKndiTHckh5cHGFRfV8hZYj5Ia90uFWyj/NkXdwnjKn8k9zeZyKvcQyuxshRvZXq4MBxk6QEfs+8I+Ze77t27UJiYiIGDhyIm2++GfHx8XjzzTcVbK3nlLp1eOeH1nRtn5bbSqQVrOAqQM2hCVKGWyIlcdiC58Te633w4MF45plnMHr0aAQFBeGPf/wjsrKyEBYWhtTUVBW2QFlG6s9cbQv3GSIGXEN1dF1JEW719v60DPOTfR39Twuyr4Nu4LCF3om91/v06dOdni9evBhvvPEG9u/fL1vA7fwZevpjnkOrvMPQS8SAKzu1qrdGD7dKBFlv1s3wKw9Wc93ruNd758vs9Hav984EQcDHH3+MmpoaPP/88y7naWtrQ1vbjfe+ubm5x9eUou/gX5+k1fV94L5ERmfqgKuVjk9r4VYr70sHNcOsWF3bysArLVZzu/PmXu8AcPHiRURHR6OtrQ0BAQHYuHEjZs6c6XLegoICrFq1StJ298ToP9C1gIGXjM7UAVduanQYRgi3egq0vWHglR6rudLo378/qqurcenSJZSXl8NqtWL48OHdhi8AQG5uLqxWq+N5c3MzYmJiFGyt57TQh+kRAy8ZjWkDbvNtFgT0PpvspKze6jXcGinQ9oaBVzqs5l7nzb3egevDGG6//XYAQHx8PL755hsUFBS4DLgWiwUWizT9Q299nrf9GIOttBh4Se94mTCZeNIZaCXcqnG5mZZhfo6HmfF98A0vleTdvd5dsdvtTuNspeRpOPKmH+N3QBm8NBnpjWkruEbia7hVCkNczzq/P6zsimP2aq7VakVmZiYSExMxadIkFBUVdbvXe3R0NAoKCgBcH1ObmJiIESNGoK2tDXv37sWbb76JV199Vc3NEI1hSx2d33cz73ekbQy4MlCyequHcMtgKx7DrnhmHpsr9l7vra2tePLJJ/Hdd9+hb9++GD16NN566y2kp6fL2s6e+j0xfRmDrXYw7JJWMeBKTMkdXMtj1RhqpdPxXjLoesas1Vwx93pfs2YN1qxZI0s7ZkVmA/eNkOW1AYZbLWPYJS1hwFWBHHcs85TcBwcGW/kw6HrOrCFXC7z54e3pMgy3+sGwS2pjwJWQ1ocmyHlwYLBVDocveMbMQxa0pLf3n+HW+Bh2uyv/JLf3mcgnDLg6pKVwy2CrLlZ1e8dqrrLc9TW+/LhnuDUOhl1SCgOuRJSq3mol3DLYaguDbs8YcrXJk/6M4da4GHZJTgy4ClEj3DLYmg+DrnscsqA/DLfmwf2TpMaAKwEldki1wy2Drb4w6LrHaq429NanMdyaE6u6JBXeyUwBvlZv1Qy3vMuWvvGzc43hSRkdAUVsH8jPhwDePY18w4Dro95+Yeo13DLYGgc/S9d40NQmfi7kCoMuicWAq2FqhlsyHgbd7njAVIe7vo2fB/WGQZc8xYDrA7mrt2JIscMzAJkDP2NnPFhqAz8HEoPDF6g3PMlMo5S8LzsDj/nwJDRnPINbPl1/6Ht7i3Eid7j/kisMuF6Ss3rLcEtKaRnmx5DbCa+woA4t3GWR+4H+MehSZwy4MtBDuGWwpQ6s5jpjyPXe3TMKgB5Chqv+TSvnDvS0PPcNfeGlxghgwPWKXDsMwy2pidXcGxhyvaP0eEil+jJX6+G+og+s6poXA67EvK3eKhFutR5s2261Sf6alrogyV/TyFjNvYEh13e99Yfe9GVa6cc6t4P7i/ZpKeiWf5KrdhNMgVdREEmOncNM4bbtVpvbh5Lro55p5fuiNp6hLZ2u/ZzY91bLV3npaJtW20c38MoL6iguLkZsbCyCg4ORlJSEgwcPerTc9u3b4efnh9mzZ4teJwOuhOS+LJi31Q41O12tBkuG3t7xYH0dD4bq09N3kWFXHxh0lbNjxw5YrVbk5+fj0KFDiIuLQ2pqKs6dO9fjcqdOncJTTz2FqVOnerVeBlwR1Kze6uVPeXoOjXptt5x4kL6OB0LfeFu91XtQ1Hv7zYBBV36FhYVYuHAhsrKyMHbsWJSUlOCmm25CaWmp22Xa29vxyCOPYNWqVRg+fLhX62XAlYg31VujhFsjBkMjbpO3eJC+jgdBZRnpO8eqrvYx6IrT3Nzs9Ghrc10AtNlsqKqqQkpKimOav78/UlJSUFlZ6fb1n332WYSFhWH+/Plet5EnmXmop+qtlsKtUh2omYJf520180lrvMoCTzzz1IDjbW77RU/6NCMHQZ7IqW1aOhlNbiEn2xAYKG5fu3bt+vsSExPjND0/Px8rV67sNn9TUxPa29sRHh7uND08PBxHjx51uY79+/djy5YtqK6uFtW2rhhwVaDXcGumUOuO2cMuQy5DrlhiTqI1crDtikFX28wUdL1RX1+PkJAQx3OLRaLrWbe0YN68edi0aRNCQ0N9ei0GXA9IXb31hFbCLUOtex3vjdmCLkMuQ663eurXzBRuO2PQ1TYGXddCQkKcAq47oaGhCAgIQGNjo9P0xsZGREREdJv/+PHjOHXqFNLS0hzT7HY7ACAwMBA1NTUYMWKER23kGFyFeVLN0EK45fhTz5nxvTJrGOmM4/Wkw+8Tx7prHcfoeicoKAgJCQkoLy93TLPb7SgvL0dycnK3+UePHo0jR46gurra8bj//vsxY8YMVFdXdxsa0RNWcHshZfVW6nArV7Al75itostKLiu5Yrjr2xjqnLGiq22s6IpntVqRmZmJxMRETJo0CUVFRWhtbUVWVhYAICMjA9HR0SgoKEBwcDDGjRvntPzAgQMBoNv03rCCqxCth1szViHlYqb3klUn7VRyxVxIfdOmTZg6dSoGDRqEQYMGISUlxeMLr4shZvwtOeO+pW2s6HouPT0dGzZsQF5eHuLj41FdXY2ysjLHiWd1dXU4e/as5OtlBbcHSo69VSvcmiWIqcFMFV1Wc9XVcSH1kpISJCUloaioCKmpqaipqUFYWFi3+SsqKjBnzhxMnjwZwcHBeP7553HPPffg66+/RnR0tM/t8bR/ZIDrHfctbWNF1zM5OTnIyclx+X8VFRU9Lrt161av1skKrhekHpqgRrg1U5VRbWZ5n80cVtSu5Ii9kPrbb7+NJ598EvHx8Rg9ejQ2b97sGBfnSltbW7frXnaVuKAQgOsDvav3x8zfF7FYzdU+tfsA6o4B1w2pfo1pLdwy2KrDLO+7mQ/Cah3gvL2QemeXL1/G1atXMXjwYJf/X1BQgAEDBjgenpzo0VPfZ+bviS8YdLWNwxa0hQFXJDHVWy2FW7MELK0zw2dg5gOwGge3ni6k3tDQ4NFrLF26FFFRUU4hubPc3FxcvHjR8aivr/e4fTzgS8/M+5ge9BR0yz/JVbg15sUxuC4oMZZGzL3YfWWGUKUnZhiba+ZxgxdHWHBzzRW1m+GxdevWYfv27aioqEBwcLDLeSwWi3QXcmc4kwSvtqB9HJ+rLk1UcLV49q8rUlVvlQq3rNpqGz8b42q+TbmqpdgLqXe2YcMGrFu3Dh988AEmTJggZzMBMNzKge+p9vGvGOpQPeB2nP2bn5+PQ4cOIS4uDqmpqTh37pzL+TvO/v3kk09QWVmJmJgY3HPPPThz5owk7ZHil5YUl8bxpdNisNUPI39WPPAqQ+yF1Du88MILWL16NcrKypCYmChde45+5+gDeWBXBsfmah/H5ypP9SEKnc/+BYCSkhLs2bMHpaWlWLZsWbf53377bafnmzdvxu9//3uUl5cjIyNDtnZ6Wr2VYtytr+GW9KftVpshhyyYeaiCksRcSB0Ann/+eeTl5WHbtm2IjY11jNXt168f+vXrJ0sbvenXvO3PjLgveYL7G9ENqgbcjrN/c3NvDLqW+uzftrY2tLXdqMq6urxNB3fVW6mueStnuGWw1T+GXPJWeno6zp8/j7y8PDQ0NCA+Pr7bhdT9/W/8we7VV1+FzWbDQw895PQ6+fn5WLlypZJN70aKvqzzaxhxn+oJx+YSXadqwO3p7N+jR4969Bq9nf1bUFCAVatW+dxWT/g67pbhlhhyyVtiLqR+6tQpWdrQuUjQuc/zpG+Tqx8za9jlPkdmp/oYXF90nP37hz/8we3Zv55e3sbX6q0a4dbI4zfNzKifK8cImoM35yAo9X036r7lDvc5MjNVK7hSnP370Ucf9Xj2r5SXt/GGnOFWK2KHnpfstU59N0Sy19I7I1ZzWVUyn576N7X6MTNcqq8D9zkyK1UDbuezf2fPng3gxtm/7v7UBlw/+/e5557Dvn37JD37t1v7fKzeyhFu1Qy2UgZZMeswc+g1Ysgl4/Ok79PCj3SzBF2OyyUzUv0qClo4+9eXS4P5ckkwrYdbJQKtJ7q2w2yB12ghlxUl4+p8iTDAfR+nhXDbmdH2MXe475GZqB5wtXr2ryfVW1/G3YoJt0odDLQSaHtjxsBrtAMwD7TmpbVw28FM1Vzue2QGqgdcQL2zf+W6fZ6ewq1eQm1POrbB6EGXIZe0KnFBIYDe/6Kl1XDbmdH2M1e475EZaCLgqiHkZBsQ6PrKC75Ub/UQbo0Qal3pvF1GDbtGO/jyQGs8HX2gns/gN0M1l/seGZ1pA64v9BhujRpq3TFyVddoIZfMQQ/V266Mvq8x5Crry81WtZtgKrq+Dq4cpLprWVeehlupr9MYO/S86cJtZ0bdfj2GBXf0XOkjz+j5+6rntnuC+x8ZFQOuSN5Ub8WEW6kYNdh5y4jvh5EOvDzI6l/ncxqM9nkaaV9zpWWYn+E+MyIG3E56q97qIdwaMchJyWjvjdEPvKQvrvpCo3xHjbIdPWHIJSPhGFwfaSHcGi20yc3I43P1jOMBjcdoodDoY3IB7odkHKzg/p031Vtfw60U420Zbr1nlPfOSCGCFST9M/pnKPV5EkQkDwZceD80wR1Pw60vOBRBGkZ5H3nAJa0x+nfSyNtn9B8pZA4MuF5yV72VO9waJZBpjRHeU6MccHlwJVIf90PSO9OPwZVyaIKc4dYIAUzrYoee1/24XKOMEeQ4QP25OMIiOhT11q9pfX80yv7mDvdD0jNTB1w9hFsGW2XxBDQi3/XU14np0/Rwd0KGXCJtMu0QhaBj3yu6PoZbfdHze8+hCqRVvuxXWh6eZZR9zh3ui6RHpg24vZGyeiu289NyR24mev4MjH7AJe1y992Tan/S6n5p9H2OIZd8UVxcjNjYWAQHByMpKQkHDx50O++mTZswdepUDBo0CIMGDUJKSkqP87vDgOuCVOHWm8vJaLXzNit+HuriQVU/evqspN6PtFoEYMgl6m7Hjh2wWq3Iz8/HoUOHEBcXh9TUVJw7d87l/BUVFZgzZw4++eQTVFZWIiYmBvfccw/OnDkjar0MuB7wNtyKodUOm/Qbco1+sCVtSFxQ2OP/y7n/6HXfJDKTwsJCLFy4EFlZWRg7dixKSkpw0003obS01OX8b7/9Np588knEx8dj9OjR2Lx5M+x2O8rLy0WtlwG3C0+veStVuGWw1Qe9fkZGCLmsGulH1++bEvuN1vZNI+xzPeH+SADQ3Nzs9Ghra3M5n81mQ1VVFVJSUhzT/P39kZKSgsrKSo/WdfnyZVy9ehWDBw8W1UZTX0WhK0+HJkgZbkk/jHAZMSKpDTje5rKfVLJ/09q+ySsrUFdfbraq3YRugo59j0B/cd9Tf/v1jBMTE+M0PT8/HytXruw2f1NTE9rb2xEeHu40PTw8HEePHvVonUuXLkVUVJRTSPYEA24PjBZuZ0Z49mUS48OG0ZK/ppZp7UDqCSMcbHlA7V1xcTHWr1+PhoYGxMXF4ZVXXsGkSZNczvv1118jLy8PVVVVOH36NF566SX86le/UrbBEtPavmmE/a4n3CfNrb6+HiEhIY7nFovroZy+WrduHbZv346KigoEBweLWpYB9++6Vm/lCrdKBVs5wqyn6zF66NXagZSo4ySOkpISJCUloaioCKmpqaipqUFYWFi3+S9fvozhw4fjZz/7GZYsWeLz+luG+Tn1gWr9dYr7JpEyQkJCnAKuO6GhoQgICEBjY6PT9MbGRkRERPS47IYNG7Bu3Tp89NFHmDBhgug2cgwuPBt3q/VwOzPiqNNDTVppB91ghHGBHPvnntiTOP7hH/4B69evx89//nPZKi9q0dLQLyPsdz3hPkm9CQoKQkJCgtMJYh0njCUnJ7td7oUXXsDq1atRVlaGxMREr9bNCq4L7q6a4IonHZhcHa4eAmTnNhqpsqvHSpHR/2RqVh0nceTm5jqmiT2JozdtbW1OJ5E0Nze7nVdLAVMLjL7fcagC9cZqtSIzMxOJiYmYNGkSioqK0NraiqysLABARkYGoqOjUVBQAAB4/vnnkZeXh23btiE2NhYNDQ0AgH79+qFfv34er9f0AdeXoQlqhFs9hFp3OtpulKCrx5CrdzyYdifFSRy9KSgowKpVq1z+n5iCgFK4bxJpR3p6Os6fP4+8vDw0NDQgPj4eZWVljj6rrq4O/v43BhS8+uqrsNlseOihh5xex92JbO6YOuDqJdzqOdS6YqSgq7cDqdGrSSSP3NxcWK03zgJvbm52Oou6oz/UUvVWS/um0fc7/vCk3uTk5CAnJ8fl/1VUVDg9P3XqlCTrNG3AtY2M6nXjvQ23UnTyRgu1rsyMOGqIkEukJl9O4vCUxWIx3FhdkhZDLmkNTzL7u67VW7XCrdlOzjLC9mqpauUJvZ/4whNbnHl7EodUOj4PLe4HWmqT3vc7Ir1hwIU2wq0Rgp4v9L7tWjqQkvlYrVZs2rQJb7zxBr755hssWrSo20kcnU9Cs9lsqK6uRnV1NWw2G86cOYPq6mrU1taqtQmy4b6pHP74JC0x7RCFDp6eINFTuPU12NJ1RhqbS/Lin0OdiT2J4/vvv8fEiRMdzzds2IANGzZg2rRp3cbD+cqM18p2x+hjcYm0xNQB19OTyuQItwy27ul1bK6WTmrpDQ+0xiPmJI7Y2FgIgjQ/ENputbnsB3vq45T+MaunfVPv+OOTtMK0QxSab1Mn3Jp9KIKn+B4R6Zen+68Z93OOxSVShmkDbldyh1sGW/H0+H7pabyf3g+0HO+nvv9zx791myZ2v1VqP9fTvql33DdJCxhwIS7cxg49L6qjZLD1Dd87Iu2yjYxyeu7t/mq2/VzvPy6J9MD0AVdsuPUUg6109PY+slJEZtF8m0VXN7ThvqkcVnFv+HKztfeZSHKmDrhyhFsGW3nwPZWH3itJPIiqqyXmxvsvxT7K/ZyIpGLagNu5Y+4gRbgl+ejp/WWliIh6ovcfl0RaZ9qA25Uv4ZZVW+XwfZYeD7TkLVvM9e+OlPul3Ps4f3wqh39hITUx4ML3cEvkCg+kRERE6jB9wPU23LJqqx6+79QZq0TqkmN/NMs+boa/nnD/JLWYNuDaYmxehVsGW23gZ0CkvlujmtRuglf41xUi4zP1rXpd6S3cKu2hkENu/29n848VbAl5Qy+3COWte4mIyEgYcP/OkyEJcugpwPqyrBnC78yIo4rdy56IlGWW/Zs/LonkwYALZau2vgRab9djhrBL5tYyzA/9TwtqN8N0ZoQfA9BH7WZ4RS9/XTEC7p+kBtMHXLnDrVKB1tM2GC3s6qHKwwMpERGRskwbcG+NakLgzRa3/+9LuNVCqHWno21GC7rkO/6plLRGDz9giUibTBtw3fE22Go51LryUMghw4RcHgSJ1NW1/zNK30JE+mXay4S54k24fSjkkO7CbQe9tluPeFkiMipX/Qj7FnF4PVwi6bGCC/HB1kidN4csEJEcjPRXIiLSH9NWcGeEHxN90wY9V2t7o/ft4o0fiBUi5c3u//9kX4ec+zb/skJkXKYNuGIYOdgSEcmF/SaZ3ZebrWo3wbQYcHtgtmBrpm0l18wwFpCIiIyPY3C7MHvI47g5+fB6uERERMpgwAVDrVHwcmFEREQEmHiIwuz+/890QxA8xfeEiIiIpFJcXIzY2FgEBwcjKSkJBw8e7HH+d999F6NHj0ZwcDDGjx+PvXv3il6naQMuEREREclrx44dsFqtyM/Px6FDhxAXF4fU1FScO3fO5fwHDhzAnDlzMH/+fBw+fBizZ8/G7Nmz8dVXX4laLwMuERERyY6X8jOnwsJCLFy4EFlZWRg7dixKSkpw0003obS01OX8L7/8Mv7xH/8Rv/71rzFmzBisXr0aP/7xj/Ef//EfotZrujG4giAAAC5dsqvcEm37R78vsbslTu1miHattU3tJvTI/rcrajehV+1X9LtvtNsEtZsAAGi3Xf+cO/oboxHTj165dNXn9cm5X2tln9TzfieGVvZRpTQ3N8v+2t72M9cEGyDya3dNsDmtu4PFYoHFYuk2v81mQ1VVFXJzcx3T/P39kZKSgsrKSpfrqKyshNXqfHm11NRU7N69W1RbTRdwW1paAADTk3iB7959oHYDvKDHNpNRtbS0YMCAAWo3Q3Li+lEp9knu16RPA958RvZ1iO1ngoKCEBERgYqG33q1vn79+iEmJsZpWn5+PlauXNlt3qamJrS3tyM8PNxpenh4OI4edX0Tl4aGBpfzNzQ0iGqn6QJuVFQU6uvr0b9/f/j5GfPPJc3NzYiJiUF9fT1CQkLUbo7kuH36Z/RtFAQBLS0tiIqKUrspsjBaP2r07yO3T9/cbZ+3/UxwcDBOnjwJm827654LgtBtv3dVvVWb6QKuv78/hg4dqnYzFBESEmLInb0Dt0//jLyNRqzcdjBqP2rk7yPA7dM7V9vnbT8THByM4OBgKZrVo9DQUAQEBKCxsdFpemNjIyIiIlwuExERIWp+d3iSGRERERFJLigoCAkJCSgvL3dMs9vtKC8vR3JysstlkpOTneYHgA8//NDt/O6YroJLRERERMqwWq3IzMxEYmIiJk2ahKKiIrS2tiIrKwsAkJGRgejoaBQUFAAAFi9ejGnTpuHFF1/Efffdh+3bt+PLL7/Ea6+9Jmq9DLgGZLFYkJ+fr8kxMVLg9umfGbaR9MPo30dun77pffvS09Nx/vx55OXloaGhAfHx8SgrK3OcSFZXVwd//xsDCiZPnoxt27Zh+fLlePrpp3HHHXdg9+7dGDdunKj1+glGvY4NEREREZkSx+ASERERkaEw4BIRERGRoTDgEhEREZGhMOASERERkaEw4OpUcXExYmNjERwcjKSkJBw8eNDtvJs2bcLUqVMxaNAgDBo0CCkpKT3OrwVitq+z7du3w8/PD7Nnz5a3gT4Su30XLlxAdnY2IiMjYbFYMHLkSOzdu1eh1npH7DYWFRVh1KhR6Nu3L2JiYrBkyRJcuXJFodaSUa1cuRJ+fn5Oj9GjR/e4zLvvvovRo0cjODgY48eP1/S+Fhsb2237/Pz8kJ2d7XL+rVu3dptXiQv+e+rTTz9FWloaoqKi4Ofnh927dzv9vyAIyMvLQ2RkJPr27YuUlBR8++23vb6ut8cUqfW0fVevXsXSpUsxfvx43HzzzYiKikJGRga+//77Hl/Tm++4GTDg6tCOHTtgtVqRn5+PQ4cOIS4uDqmpqTh37pzL+SsqKjBnzhx88sknqKysRExMDO655x6cOXNG4ZZ7Ruz2dTh16hSeeuopTJ06VaGWekfs9tlsNsycOROnTp3Czp07UVNTg02bNiE6OlrhlntO7DZu27YNy5YtQ35+Pr755hts2bIFO3bswNNPP61wy8mIfvSjH+Hs2bOOx/79+93Oe+DAAcyZMwfz58/H4cOHMXv2bMyePRtfffWVgi323H/91385bduHH34IAPjZz37mdpmQkBCnZU6fPq1Uc3vV2tqKuLg4FBcXu/z/F154Af/+7/+OkpISfPHFF7j55puRmpra449hb48pcuhp+y5fvoxDhw5hxYoVOHToEHbt2oWamhrcf//9vb6umO+4aQikO5MmTRKys7Mdz9vb24WoqCihoKDAo+WvXbsm9O/fX3jjjTfkaqJPvNm+a9euCZMnTxY2b94sZGZmCg888IACLfWO2O179dVXheHDhws2m02pJvpM7DZmZ2cLd911l9M0q9UqTJkyRdZ2kvHl5+cLcXFxHs//8MMPC/fdd5/TtKSkJOFf/uVfJG6ZPBYvXiyMGDFCsNvtLv//9ddfFwYMGKBso7wEQPjDH/7geG6324WIiAhh/fr1jmkXLlwQLBaL8M4777h9HV+PmXLpun2uHDx4UAAgnD592u08Yr/jZsEKrs7YbDZUVVUhJSXFMc3f3x8pKSmorKz06DUuX76Mq1evYvDgwXI102vebt+zzz6LsLAwzJ8/X4lmes2b7Xv//feRnJyM7OxshIeHY9y4cVi7di3a29uVarYo3mzj5MmTUVVV5fiz4YkTJ7B3717ce++9irSZjO3bb79FVFQUhg8fjkceeQR1dXVu562srHT67gJAamqqx/2rmmw2G9566y089thj8PPzczvfpUuXMGzYMMTExOCBBx7A119/rWArvXfy5Ek0NDQ4fT4DBgxAUlKS289HimOmmi5evAg/Pz8MHDiwx/nEfMfNggFXZ5qamtDe3u64A0iH8PBwNDQ0ePQaS5cuRVRUVLdOXAu82b79+/djy5Yt2LRpkxJN9Ik323fixAns3LkT7e3t2Lt3L1asWIEXX3wRa9asUaLJonmzjXPnzsWzzz6LO++8E3369MGIESMwffp0DlEgnyUlJWHr1q0oKyvDq6++ipMnT2Lq1KloaWlxOX9DQ4NP/auadu/ejQsXLuDRRx91O8+oUaNQWlqK9957D2+99RbsdjsmT56M7777TrmGeqnjMxDz+UhxzFTLlStXsHTpUsyZMwchISFu5xP7HTcL3qrXZNatW4ft27ejoqJCUycWeKulpQXz5s3Dpk2bEBoaqnZzZGG32xEWFobXXnsNAQEBSEhIwJkzZ7B+/Xrk5+er3TxJVFRUYO3atdi4cSOSkpJQW1uLxYsXY/Xq1VixYoXazSMdmzVrluPfEyZMQFJSEoYNG4bf/e53mv+Lj1hbtmzBrFmzEBUV5Xae5ORkJCcnO55PnjwZY8aMwW9+8xusXr1aiWaSB65evYqHH34YgiDg1Vdf7XFeM33HxWDA1ZnQ0FAEBASgsbHRaXpjYyMiIiJ6XHbDhg1Yt24dPvroI0yYMEHOZnpN7PYdP34cp06dQlpammOa3W4HAAQGBqKmpgYjRoyQt9EiePP5RUZGok+fPggICHBMGzNmDBoaGmCz2RAUFCRrm8XyZhtXrFiBefPmYcGCBQCA8ePHo7W1FY8//jieeeYZp/uUE/li4MCBGDlyJGpra13+f0REhFf9q9pOnz6Njz76CLt27RK1XJ8+fTBx4kS374eWdHwGjY2NiIyMdExvbGxEfHy8y2V8OWaqpSPcnj59Gh9//HGP1VtXevuOmwWPGjoTFBSEhIQElJeXO6bZ7XaUl5c7/Srv6oUXXsDq1atRVlaGxMREJZrqFbHbN3r0aBw5cgTV1dWOx/33348ZM2aguroaMTExSja/V958flOmTEFtba0juAPAsWPHEBkZqblwC3i3jZcvX+4WYjsCvSAI8jWWTOfSpUs4fvy4U0DqLDk52em7CwAffvhhj/2rFrz++usICwvDfffdJ2q59vZ2HDlyxO37oSW33XYbIiIinD6f5uZmfPHFF24/H2+PmWrpCLfffvstPvroI9xyyy2iX6O377hpqH2WG4m3fft2wWKxCFu3bhX+53/+R3j88ceFgQMHCg0NDYIgCMK8efOEZcuWOeZft26dEBQUJOzcuVM4e/as49HS0qLWJvRI7PZ1pfWrKIjdvrq6OqF///5CTk6OUFNTI/zxj38UwsLChDVr1qi1Cb0Su435+flC//79hXfeeUc4ceKE8MEHHwgjRowQHn74YbU2gQzi3/7t34SKigrh5MmTwl/+8hchJSVFCA0NFc6dOycIQvfv4l/+8hchMDBQ2LBhg/DNN98I+fn5Qp8+fYQjR46otQm9am9vF2699VZh6dKl3f6v6/atWrVK2Ldvn3D8+HGhqqpK+PnPfy4EBwcLX3/9tZJNdqulpUU4fPiwcPjwYQGAUFhYKBw+fNhxFYF169YJAwcOFN577z3hv//7v4UHHnhAuO2224S//e1vjte46667hFdeecXxvLf+SCvbZ7PZhPvvv18YOnSoUF1d7XS8bmtrc7t9vX3HzYoBV6deeeUV4dZbbxWCgoKESZMmCZ9//rnj/6ZNmyZkZmY6ng8bNkwA0O2Rn5+vfMM9JGb7utJ6wBUE8dt34MABISkpSbBYLMLw4cOF5557Trh27ZrCrRZHzDZevXpVWLlypTBixAghODhYiImJEZ588knhr3/9q/INJ0NJT08XIiMjhaCgICE6OlpIT08XamtrHf/van/73e9+J4wcOVIICgoSfvSjHwl79uxRuNXi7Nu3TwAg1NTUdPu/rtv3q1/9yrFfhoeHC/fee69w6NAhBVvbs08++cTl8apjG+x2u7BixQohPDxcsFgswt13391tu4cNG9bt+NZTf6Sknrbv5MmTLv8PgPDJJ584XqPr9vX2HTcrP0Hg3/+IiIiIyDg4BpeIiIiIDIUBl4iIiIgMhQGXiIiIiAyFAZeIiIiIDIUBl4iIiIgMhQGXiIiIiAyFAZeIiIiIDIUBl4iIiIgMhQGXiIiIiAyFAZeIiIiIDIUBl4iIiIgMhQGXyI3Y2FgUFRU5TYuPj8fKlStVaQ8Rka+mT5+OnJwc5OTkYMCAAQgNDcWKFSsgCILaTSOSFAMuERGRibzxxhsIDAzEwYMH8fLLL6OwsBCbN29Wu1lEkgpUuwFERESknJiYGLz00kvw8/PDqFGjcOTIEbz00ktYuHCh2k0jkgwruERERCbyk5/8BH5+fo7nycnJ+Pbbb9He3q5iq4ikxYBL5Ia/v3+3cWlXr15VqTVERETkKQZcIjeGDBmCs2fPOp43Nzfj5MmTKraIiMh3X3zxhdPzzz//HHfccQcCAgJUahGR9Bhwidy466678Oabb+Kzzz7DkSNHkJmZyQMAEeleXV0drFYrampq8M477+CVV17B4sWL1W4WkaR4khmRG7m5uTh58iT+6Z/+CQMGDMDq1atZwSUi3cvIyMDf/vY3TJo0CQEBAVi8eDEef/xxtZtFJCk/gRe/IyIiMoXp06cjPj6+2zW+iYyGQxSIiIiIyFAYcImIiIjIUDhEgYiIiIgMhRVcIiIiIjIUBlwiIiIiMhQGXCIiIiIyFAZcIiIiIjIUBlwiIiIiMhQGXCIiIiIyFAZcIiIiIjIUBlwiIiIiMpT/D7A8ARuQPFU+AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from few.utils.mappings import apex_of_uwyz, z_of_a, y_of_x\n", "\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", "u_grid, w_grid = np.asarray(np.meshgrid(uv, wv, indexing=\"ij\")).reshape(2, -1)\n", "a_map, p_map, e_map, xI_map = 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])\n", "\n", "# a must be a float\n", "teuk_modes = kerr_amp_spline(a, p_map, e_map, xI_map).reshape(\n", " uv.size, wv.size, -1\n", ") # a, p, e, xI\n", "\n", "# look at the contours of the (2,2,0) mode\n", "\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_modes[:, :, 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_map.reshape(uv.size, wv.size),\n", " e_map.reshape(uv.size, wv.size),\n", " np.abs(teuk_modes[:, :, 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": "markdown", "metadata": {}, "source": [ "Specific modes can be selected by providing a list of tuples of $(\\ell, m, n)$ values as the `specific_modes` kwarg:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Indices of interest: [1165 3497]\n", "True\n", "True\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAGCCAYAAAALybB5AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVpFJREFUeJzt3Xt4FPW9P/D3JrAbEBLAmBusRJCrkARDiRuKxBpNkVI5PYemeAnmaDzHJueHbm0hCoSbhFYI8dhIFIjYWgqVKvYpNzE1WkoUTcgpciDINRGyAU4lNyWLu/P7g7Jkk93Nzu7szszO+/U8+zzsZGbnOyH7nfd89jvf1QmCIICIiIiIKESEyd0AIiIiIiIpMeASERERUUhhwCUiIiKikMKAS0REREQhhQGXiIiIiEIKAy4RERERhRQGXCIiIiIKKQy4RERERBRSGHCJiIiIKKQw4BIRERFRSGHAJSIiIlKojz76CLNmzUJCQgJ0Oh127NgR0P0tXboUOp3O6TF27NiA7jMQGHCJiIiIFKqjowPJyckoKysL2j7vuOMONDU1OR779+8P2r6l0kfuBhARERGRazNmzMCMGTPc/ryzsxPPP/88fv/73+Py5cuYMGECfvnLXyIjI8Pnffbp0wdxcXE+b68ErOASERERqVRBQQGqq6uxdetW/P3vf8ecOXPw/e9/H1988YXPr/nFF18gISEBI0aMwMMPP4yGhgYJWxwcOkEQBLkbQURERESe6XQ6vPPOO5g9ezYAoKGhASNGjEBDQwMSEhIc62VmZmLKlClYtWqV6H3s3r0b7e3tGDNmDJqamrBs2TKcO3cOn3/+OQYOHCjVoQQchygQERERqdDhw4dhs9kwevRop+WdnZ24+eabAQDHjh3DuHHjPL7OggULsHr1agBwGg6RlJSEtLQ0DB8+HH/4wx/w+OOPS3wEgcOAS0RERKRC7e3tCA8PR01NDcLDw51+NmDAAADAiBEjcPToUY+vcz0MuzJo0CCMHj0aJ06c8L/BQcSAS0RERKRCkyZNgs1mw4ULFzBt2jSX6+j1er+m+Wpvb8fJkyfx6KOP+vwacmDAJSIiIlKo9vZ2p+rp6dOnUVdXhyFDhmD06NF4+OGHkZOTg7Vr12LSpEm4ePEiKisrkZSUhJkzZ4re37PPPotZs2Zh+PDhOH/+PIqKihAeHo65c+dKeVgBx5vMiIiIiBSqqqoK99xzT4/l8+bNw+bNm3H16lWsXLkSv/nNb3Du3DlER0fjrrvuwrJlyzBx4kTR+/vJT36Cjz76CP/3f/+HW265Bd/97nfxwgsvYOTIkVIcTtAw4BIRERFRwK1evRqFhYWYP38+SktL3a731ltvYfHixThz5gxGjRqFX/7yl3jggQdE7Yvz4BIRERFRQH366ad49dVXkZSU5HG9AwcOYO7cuXj88cdx6NAhzJ49G7Nnz8bnn38uan+s4BIRERFRwLS3t+POO+/EK6+8gpUrVyIlJcVtBTc7OxsdHR3485//7Fh21113ISUlBeXl5V7vU3M3mdntdpw/fx4DBw6ETqeTuzlEFIIEQUBbWxsSEhIQFhZ6H5SxHyWSnz/9zJUrV2C1Wn3eb/f3vcFggMFgcLtNfn4+Zs6ciczMTKxcudLj61dXV8NsNjsty8rKwo4dO0S1U3MB9/z58zAajXI3g4g0oLGxEcOGDZO7GZJjP0qkHGL7mStXruDWW2/CxYt2n/Y3YMAAtLe3Oy0rKirC0qVLXa6/detW1NbW4tNPP/Xq9S0WC2JjY52WxcbGwmKxiGqn5gLu9a+Zy4jOQZ8wfcD2Yx2d0GOZ/vh5yV9TCv62K5A8HbOS203a9q3diqpLv1HV11qKcf24GhsbERkZKXNriELPv476Wa/r+NrPWK1WXLxoR9UnMRgwQNwnMO3tAjLSLvR477ur3jY2NmL+/PnYt28fIiIiRO3LX5oLuNfL6n3C9AELuNaxw1z+Yv3Zn7vX9Jf+2JdAAIO+WNaxzlehno7ZPn6E49/6Y18GqEVEvgvVj++vH1dkZCQDLpEEZsTnOz0Xkxd87WcGDNBhwECxQ6iuVX29fe/X1NTgwoULuPPOOx3LbDYbPvroI/z6179GZ2dnj29gi4uLQ3Nzs9Oy5uZmxMXFiWqp5gIu3aCUUNg91PrzGko5JiIiIne6B9pQde+99+Lw4cNOy3JzczF27FgsWLCgR7gFAJPJhMrKSjz99NOOZfv27YPJZBK1bwZcibkLa/4ELykCoBIF4rgYdImISGm0Emi7GzhwICZMmOC07KabbsLNN9/sWJ6Tk4OhQ4eiuLgYADB//nxMnz4da9euxcyZM7F161Z89tlneO2110TtmwFXo+QMgMEI7NaxwxhyiYhINloNtWI1NDQ4zQKRnp6OLVu2YNGiRXjuuecwatQo7Nixo0dQ7g0DrsIFIgzKFfxCtRJNRETEQOudqqoqj88BYM6cOZgzZ45f+2HAlVAghieECjnCLau4REQUKAy0ysaAq2ChUL2Vu2rLkEtERFJgoFUXBlwKGLnDLRERkT8YatWLATfAfK0eqr16y3BLRERqo6ZA+8cv1iIqaqPczVAsBlyNCFa4ZbAlIiI1UVOoJe8x4CqQWkOiWttNRETawUCrDQy4GhCM6i3DLRERKRVDrfYw4AaQL8FSjUFRjW0mIqLQxUBLDLghLtDVW4ZbIiJSAoZa6ooBN4Qx3BIRUShjqCV3GHADJNSHJ6iprUREFDoYaskbDLgBoIRvzgpkGxhuiYgoWBhoyRcMuAqhltColnYSEZF6MdSSvxhwJaSEyi0QuHaoMdwq5f+EiIg8Y6glKTHgklfUGG6JiEjZGGp9s7upDK2trXI3Q9EYcBVAyvDIiiURESkZQy0FAwMu9Uqt1VuGfSIiZWCopWBjwA0hgQh0ag23REQkL4ZakhMDrsykCpAMt85YvSUiCj6GWlIKBlxySc3hloiIgoehlpSIATcEsFrpjL8PIqLAY7AlJWPAlZFSq6RKbRcREcmLoZbUggFX5aSuVqo93LJ6S0QkLYZaUqMwuRtAyqH2cEukZWVlZUhMTERERATS0tJw8OBBj+uXlpZizJgx6NevH4xGI5555hlcuXIlSK0lpZsRn+94EKkRK7gykSJMslrpjL8P0qpt27bBbDajvLwcaWlpKC0tRVZWFurr6xETE9Nj/S1btmDhwoWoqKhAeno6jh8/jsceeww6nQ4lJSUyHAEpBQMthQoGXAKg/uotwy1pWUlJCfLy8pCbmwsAKC8vx86dO1FRUYGFCxf2WP/AgQOYOnUqHnroIQBAYmIi5s6di08++SSo7SZlYKilUMQhCiolZaBTe7gl0jKr1YqamhpkZmY6loWFhSEzMxPV1dUut0lPT0dNTY1jGMOpU6ewa9cuPPDAAy7X7+zsRGtrq9OD1I9DENRpd1OZ3E1QBQZcUj1Wb0nLLl26BJvNhtjYWKflsbGxsFgsLrd56KGHsHz5cnz3u99F3759MXLkSGRkZOC5555zuX5xcTGioqIcD6PRKPlxUHBwbC0F0/r165GUlITIyEhERkbCZDJh9+7dbtffvHkzdDqd0yMiIsKnfTPgysDfiimrt0Tkj6qqKqxatQqvvPIKamtr8fbbb2Pnzp1YsWKFy/ULCwvR0tLieDQ2Nga5xeQPhlqSy7Bhw7B69WrU1NTgs88+w/e+9z08+OCDOHLkiNttIiMj0dTU5HicPXvWp31zDK6GhUK4ZfWWtC46Ohrh4eFobm52Wt7c3Iy4uDiX2yxevBiPPvoonnjiCQDAxIkT0dHRgSeffBLPP/88wsKcax8GgwEGgyEwB0ABw0BLcps1a5bT8xdeeAHr16/Hxx9/jDvuuMPlNjqdzm3fJQYruCrDQHcDfxdEgF6vR2pqKiorKx3L7HY7KisrYTKZXG7z9ddf9wix4eHhAABBEALXWAo4VmspGLqPye/s7Ox1G5vNhq1bt6Kjo8Nt3wQA7e3tGD58OIxGY6/VXk80XcG9XsEMZlBSStVUKe3wFcMt0Q1msxnz5s3D5MmTMWXKFJSWlqKjo8Mxq0JOTg6GDh2K4uJiANeqKiUlJZg0aRLS0tJw4sQJLF68GLNmzXIEXVIXBloSa0dbMiKEvqK2udJ+FcB7PcbhFxUVYenSpS63OXz4MEwmE65cuYIBAwbgnXfewfjx412uO2bMGFRUVCApKQktLS1Ys2YN0tPTceTIEQwbJi63aDrgqo1UoU7t4ZaInGVnZ+PixYtYsmQJLBYLUlJSsGfPHseNZw0NDU4V20WLFkGn02HRokU4d+4cbrnlFsyaNQsvvPCCXIdAPmKwJTk0NjYiMjLS8dzTEKYxY8agrq4OLS0t2L59O+bNm4cPP/zQZcg1mUxO1d309HSMGzcOr776qtt7BNxhwCXVYfWWqKeCggIUFBS4/FlVVZXT8z59+qCoqAhFRUVBaBlJjaGW5HZ9VgRv6PV63H777QCA1NRUfPrpp3jppZfw6quv9rpt3759MWnSJJw4cUJ0GxlwVYLV22sYbolIqxhsKRTY7XavxuwC18btHj582O0c3Z4w4AaR2sMlEREFH4MtqVVhYSFmzJiBW2+9FW1tbdiyZQuqqqqwd+9eAD3vD1i+fDnuuusu3H777bh8+TJefPFFnD171jHjixgMuLgWPJVcGWT19hol/x8REUmJoZZCwYULF5CTk4OmpiZERUUhKSkJe/fuxX333Qeg5/0BX331FfLy8mCxWDB48GCkpqbiwIEDbm9K84QBl1SB4ZaItIDBljxR29f0btq0yePPu98fsG7dOqxbt06SfWs24FpHJwT14OWunsq9f38w3BJRqGOwJZKWZgOuWkgR7hhuiYiUicGWKDAYcImIiIKMwZYosBhwg8DXCiqrt6zeElHoYKglCh4GXFIkhlsiChUMtkTBx4CrUFqu3jLcElEoYLAlkg8DLgIbqNQaMomIyDcMtkTyY8ANUWoN1qzeEpFaMdgSKQcDrgJpNeRp9biJSN0YbImUhwG3F10roWIDmFxVVDVWbxluiUhtGGyJlIsBl2THcEtEasJgS3JQ29f0yo0BV2H8DXtqrN4SEakBgy2ReoTJ3YCysjIkJiYiIiICaWlpOHjwoMf1S0tLMWbMGPTr1w9GoxHPPPMMrly5EqTWeo9B0zus3hKR0s2Iz2e4JVIZWQPutm3bYDabUVRUhNraWiQnJyMrKwsXLlxwuf6WLVuwcOFCFBUV4ejRo9i0aRO2bduG5557LsgtDwythT2tHS8RqQ+DLZE6yRpwS0pKkJeXh9zcXIwfPx7l5eXo378/KioqXK5/4MABTJ06FQ899BASExNx//33Y+7cub1WfT3xFLLUVoVVU3sZbolIyVi1JVI32QKu1WpFTU0NMjMzbzQmLAyZmZmorq52uU16ejpqamocgfbUqVPYtWsXHnjgAbf76ezsRGtrq9Mj0NQUNImI6AYGW6LQINtNZpcuXYLNZkNsbKzT8tjYWBw7dszlNg899BAuXbqE7373uxAEAd9++y3+8z//0+MQheLiYixbtszlz8RUEQNdcdTSzWWs3hKR0jDUEoUW2W8yE6OqqgqrVq3CK6+8gtraWrz99tvYuXMnVqxY4XabwsJCtLS0OB6NjY1BbDF1x3BLRErCii1RaJKtghsdHY3w8HA0Nzc7LW9ubkZcXJzLbRYvXoxHH30UTzzxBABg4sSJ6OjowJNPPonnn38eYWE987rBYIDBYOixXH/8PBCml+BInMlRSVVL9ZbhloiUhMGWKHTJVsHV6/VITU1FZWWlY5ndbkdlZSVMJpPLbb7++useITY8PBwAIAhC4BobYAx+ROQvMVMuZmRkQKfT9XjMnDkziC2WD6u2RKFP1i96MJvNmDdvHiZPnowpU6agtLQUHR0dyM3NBQDk5ORg6NChKC4uBgDMmjULJSUlmDRpEtLS0nDixAksXrwYs2bNcgRdrWH1loiuT7lYXl6OtLQ0lJaWIisrC/X19YiJiemx/ttvvw2r1ep4/n//939ITk7GnDlzgtnsoGOoJdIOWQNudnY2Ll68iCVLlsBisSAlJQV79uxx3HjW0NDgVLFdtGgRdDodFi1ahHPnzuGWW27BrFmz8MILL8h1CE7UEjaDjeGWKLC6TrkIAOXl5di5cycqKiqwcOHCHusPGTLE6fnWrVvRv3//kA24DLZE2iP7V/UWFBSgoKDA5c+qqqqcnvfp0wdFRUUoKioKQstuCGRA8+e11RCoGW6JAuv6lIuFhYWOZb1Nudjdpk2b8JOf/AQ33XSTy593dnais7PT8TwY0y1KheGWQsHupjK5m6A6qppFIdgYzohI6TxNuWixWHrd/uDBg/j8888dN++6UlxcjKioKMfDaDT63e5A4zhbIm2TvYKrdAy5vuPvjkj5Nm3ahIkTJ2LKlClu1yksLITZbHY8b21tVWzI7RpqXX3KxX6JSBtYwZVIsIcLKH14Ak8iRMHhy5SL13V0dGDr1q14/PHHPa5nMBgQGRnp9FCi3sKtp+VEFFoYcGXEEEhE/vJlysXr3nrrLXR2duKRRx4JdDMDqutwBOvYYQyxRMQhCmqk9M6bwZ0ouMROuXjdpk2bMHv2bNx8881yNFsS3lRtiUh7WMElSTHcEgVfdnY21qxZgyVLliAlJQV1dXU9plxsampy2qa+vh779+/vdXiCUvEmMiLlW79+PZKSkhxDm0wmE3bv3u1xm7feegtjx45FREQEJk6ciF27dvm0b1ZwJeBL1YBBkIikJGbKRQAYM2aMar8B0lWwZfWWSHmGDRuG1atXY9SoURAEAW+88QYefPBBHDp0CHfccUeP9Q8cOIC5c+eiuLgYP/jBD7BlyxbMnj0btbW1mDBhgqh9s4KrMkruxBnaiSiQ3FVtldwvEmnZrFmz8MADD2DUqFEYPXo0XnjhBQwYMAAff/yxy/VfeuklfP/738fPf/5zjBs3DitWrMCdd96JX//616L3zYBLkmC4JaJAcjccgeGWKPhaW1udHl2/CMYdm82GrVu3oqOjw+0NsNXV1cjMzHRalpWV5fWX1nTFIQr/dL2TFBvUgtm5siMnIq3hOFuiwPigeTT6tBtEbfNtRyeA93rMg11UVISlS5e63Obw4cMwmUy4cuUKBgwYgHfeeQfjx493ua7FYvH5S2u6Y8CFc3C0jh0W8GpkqFU7Q+14iEgZ7r2nGPDQJ/Oin7RAiV/T29jY6DQftsHgPiiPGTMGdXV1aGlpwfbt2zFv3jx8+OGHbkOuVDQfcNlB+ofhloikdu89N6YzYx9DpDxivvBFr9fj9ttvBwCkpqbi008/xUsvvYRXX321x7pxcXE+fWmNK5oeg6umcKumthIR+apruPWEfSKROtntdrdjdk0mk9OX1gDAvn37ev3SGlc0X8H1h9Y7WFZWiEhK3oZbf7DfIgqewsJCzJgxA7feeiva2tqwZcsWVFVVYe/evQB6fgnN/PnzMX36dKxduxYzZ87E1q1b8dlnn+G1114TvW8G3CDzpXPVepCWi7e/d54wifzjuJHMy/cc+0Qidbhw4QJycnLQ1NSEqKgoJCUlYe/evbjvvvsAXPsSmrCwG4MJ0tPTsWXLFixatAjPPfccRo0ahR07doieAxfQcMC1jk7Q7sFLIJRCnb8nS1fbh9LvhyiQPN1IxvcRkbpt2rTJ489dfQnNnDlzMGfOHL/3zYznIy1XEELhpBPo/7+urx8Kvy+iQJgRn+8It96+J/157/K9SKQdDLgKp+UgLTW5fpcMu0RERMHFgBtEoRBu1HgMSrpI8PULRYhIWe9lIlI2BlwfBKuTZWfuHyX//hh0iTyT+r3B9xqRtmh6HlwSRy0nCOvYYYoOt12ppZ1EUpsRn8+/fyIKGAZcChlqCrZdqbHNRFLy5j3A9wkRicEhCiL52smKrX4qrTNXcvVWab8rX3DIAlHg8H1FarW7qUzuJqgWK7jdsCNUl1AIt12F2vEQ+aJ7P8z3BRGJxYBLvVJi6FfrcARvhOpxEbmixP6FiNSPAVcEzp6gDFr4/WjhGImCgQGaSJsYcINAzR2s0tqupeCnpWMlcofvAyLyBQOul1i9lVcoD0nwRIvHTNrS/W9cyotqpV2gE1HwMOCSW0o5OWg95Gn9+ImIiMRiwCVFY7gj0i6+/4nIVwy4XvCnkxVTBVVSZ66E6q2Sfh9y4++CelNWVobExEREREQgLS0NBw8e9Lj+5cuXkZ+fj/j4eBgMBowePRq7du0KUmsDTwl9GBHJhwGXFImBrif+Tsidbdu2wWw2o6ioCLW1tUhOTkZWVhYuXLjgcn2r1Yr77rsPZ86cwfbt21FfX48NGzZg6NChQW03QygRBQq/yawLOTtbJYUXuU86SvpdEKlBSUkJ8vLykJubCwAoLy/Hzp07UVFRgYULF/ZYv6KiAv/4xz9w4MAB9O3bFwCQmJgYtPbOiM93ubxr3xOsT86IKDSxgtsLdrLBxXDrGX8/1J3VakVNTQ0yMzMdy8LCwpCZmYnq6mqX2/zpT3+CyWRCfn4+YmNjMWHCBKxatQo2m83l+p2dnWhtbXV6EFFg8Wt6/cOAS05YxSZSl0uXLsFmsyE2NtZpeWxsLCwWi8ttTp06he3bt8Nms2HXrl1YvHgx1q5di5UrV7pcv7i4GFFRUY6H0WiUpO2BeM+zsEBEAAOuIjDY8XcgBn9X5C+73Y6YmBi89tprSE1NRXZ2Np5//nmUl5e7XL+wsBAtLS2OR2NjY5BbTEQkDsfgeqC1ICFX5UNrv2ciKUVHRyM8PBzNzc1Oy5ubmxEXF+dym/j4ePTt2xfh4eGOZePGjYPFYoHVaoVer3da32AwwGAwSN94N9gnEJG/WMENEG/DotY7cq0fv6/4e6Pr9Ho9UlNTUVlZ6Vhmt9tRWVkJk8nkcpupU6fixIkTsNvtjmXHjx9HfHx8j3AbLFJcYHN4AhFdx4D7T1rvGOU4foY0ImmYzWZs2LABb7zxBo4ePYqnnnoKHR0djlkVcnJyUFhY6Fj/qaeewj/+8Q/Mnz8fx48fx86dO7Fq1Srk57ue3YCISG04RIFkwXDrP+vYYZq/MKNrsrOzcfHiRSxZsgQWiwUpKSnYs2eP48azhoYGhIXdqGcYjUbs3bsXzzzzDJKSkjB06FDMnz8fCxYskOsQHHztG/heIKKuGHADgMMTPNPqcRMFUkFBAQoKClz+rKqqqscyk8mEjz/+OMCtIiKSB4couKGlEMbKBxEREYUSBlyZaClAd6XV4w4U/j5Jrbr+7fp7kc2LdCLqjgFX44J5YmAYIyJ3X9MLsI8gCjXFxcX4zne+g4EDByImJgazZ89GfX29x202b94MnU7n9IiIiBC9bwZc9Ax5/Hpe6fHERUSBwD6XSLk+/PBD5Ofn4+OPP8a+fftw9epV3H///ejo6PC4XWRkJJqamhyPs2fPit43bzKTgVLCXrBODEo5XiKSH2f/IOrd7qYyuZsgiT179jg937x5M2JiYlBTU4O7777b7XY6nc7tF9V4ixVcIiJSJQZlInm0trY6PTo7O73arqWlBQAwZMgQj+u1t7dj+PDhMBqNePDBB3HkyBHRbWQFV0LsbHti9TbwWBEjIiKxGs5HI6yfuLGt9m+uALg2l3ZXRUVFWLp0qedt7XY8/fTTmDp1KiZMmOB2vTFjxqCiogJJSUloaWnBmjVrkJ6ejiNHjmDYMO8zheYDrpTjb72hlMAXjECklGMlIuW63hexvyBSj8bGRkRGRjqeGwyGXrfJz8/H559/jv3793tcz2QyOX3NeHp6OsaNG4dXX30VK1as8LqNmg+4atAy8tofTtRJ7z4CUAKerIgokPipBZF8IiMjnQJubwoKCvDnP/8ZH330kagqLAD07dsXkyZNwokTJ0Rtp+kxuMHuIH0JfdfDbfd/+4MnBiKSEy+AibRBEAQUFBTgnXfewV/+8hfcdtttol/DZrPh8OHDiI+PF7Wd7AG3rKwMiYmJiIiIQFpaGg4ePOhx/cuXLyM/Px/x8fEwGAwYPXo0du3aFaTWuheI0ChVoA02nryIiIgoPz8fb775JrZs2YKBAwfCYrHAYrHgm2++cayTk5ODwsJCx/Ply5fjvffew6lTp1BbW4tHHnkEZ8+exRNPPCFq37IOUdi2bRvMZjPKy8uRlpaG0tJSZGVlob6+HjExMT3Wt1qtuO+++xATE4Pt27dj6NChOHv2LAYNGiR63/rj54EwvQRHERgMt0SkFWL7DX4KRaQO69evBwBkZGQ4LX/99dfx2GOPAQAaGhoQFnaj3vrVV18hLy8PFosFgwcPRmpqKg4cOIDx48eL2resAbekpAR5eXnIzc0FAJSXl2Pnzp2oqKjAwoULe6xfUVGBf/zjHzhw4AD69u0LAEhMTJSsPYEMZ2JeO5DhlicGIiIiCgZBEHpdp6qqyun5unXrsG7dOr/3LdsQBavVipqaGmRmZt5oTFgYMjMzUV1d7XKbP/3pTzCZTMjPz0dsbCwmTJiAVatWwWazud1PZ2dnj/nalKy3cKvkyi6rt0Qkhi8X3LxIJyJvyBZwL126BJvNhtjYWKflsbGxsFgsLrc5deoUtm/fDpvNhl27dmHx4sVYu3YtVq5c6XY/xcXFiIqKcjy6z90mhd46XG+DX6DDayBPDAy3RNSbe+8plrsJRKQRqpomzG63IyYmBq+99hrCw8ORmpqKc+fO4cUXX0RRUZHLbQoLC2E2mx3PW1tbXYZctQQ0NU0VRkSkRt6eD1hNJlIu2QJudHQ0wsPD0dzc7LS8ubnZ7fcPx8fHo2/fvggPD3csGzduHCwWC6xWK/T6njeNGQwGryYglpuShx70Ri0XB0SkPGL6D6V9CtV1G4ZdksrupjK5mxASZBuioNfrkZqaisrKSscyu92OyspKp2+w6Grq1Kk4ceIE7Ha7Y9nx48cRHx/vMtwGgxTDE4IRbgPV+TLcEpGaWccOk6Qfk+p1iEgass6DazabsWHDBrzxxhs4evQonnrqKXR0dDhmVeg+N9pTTz2Ff/zjH5g/fz6OHz+OnTt3YtWqVcjPz/erHeyUiIi0JVCBlOcTImWQdQxudnY2Ll68iCVLlsBisSAlJQV79uxx3HjWfW40o9GIvXv34plnnkFSUhKGDh2K+fPnY8GCBXIdgkdKqd4GCjtyIvKV2E+VpPwUKtB9l3XsMA5ZIJKZ7DeZFRQUoKCgwOXPus+NBgAmkwkff/xxgFsVWgLR0TLcEpHaBLPfYsglkpfsX9UrN386PH87LzVXb4mI1ESOi3IWAojko/mAGyi9dWzBCres3hKRUgVrOi72WUTaw4BLovBEQURqInefJff+ibRK0wE3UB1PoKq3/JIHInKnrKwMiYmJiIiIQFpaGg4ePOh23c2bN0On0zk9IiIigtja4FBKuFRKO4i0RPabzNRKDTcPSN1GuTppby8ItHoBoIa/RQqsbdu2wWw2o7y8HGlpaSgtLUVWVhbq6+sRExPjcpvIyEjU19c7nut0umA1VxRf/74ZKom0TbMVXOvoBFn2q9Yby4J9smgZaXA8xG5DpDUlJSXIy8tDbm4uxo8fj/LycvTv3x8VFRVut9HpdIiLi3M8rk/PGCgz4v2br1wMJYZbJbaJKJRpNuAGiqdOjOGrd1KEVP6eSUusVitqamqQmZnpWBYWFobMzExUV1e73a69vR3Dhw+H0WjEgw8+iCNHjrhdt7OzE62trU4Pf+iPfcnAR0QBxYDrAzV8JKymSdEB6auvDLmkFZcuXYLNZutRgY2NjYXFYnG5zZgxY1BRUYF3330Xb775Jux2O9LT0/Hll677jeLiYkRFRTkeRqNR8uNwxZd+jMGZ1Gx3U5ncTQgZDLhB4m/gCuXxpYEKowy5RK6ZTCbk5OQgJSUF06dPx9tvv41bbrkFr776qsv1CwsL0dLS4ng0NjYGucXeUXq4VXr7iEIJbzKTkFI6L7VUb4MRQFtGGkL64oAoOjoa4eHhaG5udlre3NyMuLg4r16jb9++mDRpEk6cOOHy5waDAQaDsi8YldL/EpEysIL7T4EMW2qtJKo93GqBGobLUGDp9XqkpqaisrLSscxut6OyshImk8mr17DZbDh8+DDi4+MD1UzR/Qn/tonIHwy4XXgTutx1uqweeC/Y4ZZhmkKd2WzGhg0b8MYbb+Do0aN46qmn0NHRgdzcXABATk4OCgsLHesvX74c7733Hk6dOoXa2lo88sgjOHv2LJ544gm5DsEv7H+JqDsOUYDyq7diPmKXquoRiBMGgyZRYGRnZ+PixYtYsmQJLBYLUlJSsGfPHseNZw0NDQgLu1HP+Oqrr5CXlweLxYLBgwcjNTUVBw4cwPjx44PSXin7F4ZbInJF8wG3e+jimM3AYLiVHj/Cpa4KCgpQUFDg8mdVVVVOz9etW4d169YFoVVERPLQ9BAFsaFL7PCEYIc6pVZvlRBuldAGIvKun/K2L2P1lojc0XTApZ5CMdwSESkFQzlRcDDgBgiDHX8HRBQ4DIpE5IlmA27rbYEdniAVb8cDSzE8QcpjYbgNLI6/JSIick+zAdcdKW4wU2O400K4VWq7iELdvfcUO/7dW1/jzcUbq7dE1BsG3BCgpGoeQyQREZF4u5vK5G5CSPEp4Obk5OD111/HyZMnpW6PIil9eIK/pDoOhtvgUNIFDVGwqb16y/cvaUlxcTG+853vYODAgYiJicHs2bNRX1/f63ZvvfUWxo4di4iICEycOBG7du0SvW+fAq5er0dxcTFGjRoFo9GIRx55BBs3bsQXX3zhy8uFFK2GPK0eN5E//vrXv+KRRx6ByWTCuXPnAAC//e1vsX//fplbRkTkvw8//BD5+fn4+OOPsW/fPly9ehX3338/Ojo63G5z4MABzJ07F48//jgOHTqE2bNnY/bs2fj8889F7dungLtx40YcP34cjY2N+NWvfoUBAwZg7dq1GDt2LIYNU/fVtdr4Ww2QohrCcEsk3h//+EdkZWWhX79+OHToEDo7r31i09LSglWrVsncOnn01p+pvXpLpDV79uzBY489hjvuuAPJycnYvHkzGhoaUFNT43abl156Cd///vfx85//HOPGjcOKFStw55134te//rWoffs1Bnfw4MG4+eabMXjwYAwaNAh9+vTBLbfc4s9LKo6Y4QlqC3oMt+rDjzdDx8qVK1FeXo4NGzagb9++juVTp05FbW2tjC0LHP79EoWG1tZWp8f1C/TetLS0AACGDBnidp3q6mpkZmY6LcvKykJ1dbWoNvr0Vb3PPfccqqqqcOjQIYwbNw7Tp0/HwoULcffdd2Pw4MG+vCR14834W7lPFgy3RL6rr6/H3Xff3WN5VFQULl++HPwGKRyrt0TS0jfqER6hF7WN7YodAGA0Gp2WFxUVYenSpR63tdvtePrppzF16lRMmDDB7XoWiwWxsbFOy2JjY2GxWES11aeAu3r1atxyyy0oKirCj370I4wePdqXl1EtVm/Vd7yhQO4LGpJWXFwcTpw4gcTERKfl+/fvx4gRI+RpVIAxpBKFhsbGRkRGRjqeGwy9Z4L8/Hx8/vnnQbvHwKeAe+jQIXz44YeoqqrC2rVrodfrMX36dGRkZCAjI0O1gbd71ZSBwjWGWyL/5eXlYf78+aioqIBOp8P58+dRXV2NZ599FosXL5a7eUHnqb9lMCZSlsjISKeA25uCggL8+c9/xkcffdTrvVpxcXFobm52Wtbc3Iy4uDhRbfQp4CYnJyM5ORn/7//9PwDA//zP/2DdunXIz8+H3W6HzWbz5WVVS+rAF+jhCf6cLBhuiaSxcOFC2O123Hvvvfj6669x9913w2Aw4Nlnn8V//dd/yd08IiK/CYKA//qv/8I777yDqqoq3Hbbbb1uYzKZUFlZiaefftqxbN++fTCZTKL27VPAFQQBhw4dQlVVFaqqqrB//360trYiKSkJ06dP9+UlVUPtlQSGW3XipwmhR6fT4fnnn8fPf/5znDhxAu3t7Rg/fjwGDBggd9MoQPg+Jq3Jz8/Hli1b8O6772LgwIGOcbRRUVHo168fgGvfrTB06FAUF1/7xsP58+dj+vTpWLt2LWbOnImtW7fis88+w2uvvSZq3z4F3CFDhqC9vR3JycmYPn068vLyMG3aNAwaNMiXlyOR5OgkGW6JAkOv12P8+PFyN0Ox1F5UINKy9evXAwAyMjKclr/++ut47LHHAAANDQ0IC7sxqVd6ejq2bNmCRYsW4bnnnsOoUaOwY8cOjzemueJTwH3zzTcxbdo0UeMv1MbbECnH8ARf8UShTqz6kJrde0+x3E0gIpkIgtDrOlVVVT2WzZkzB3PmzPFr3z7Ngztz5syQDrfuaDUgsnpLRP7ydKEm19ehE1Ho8uuLHkKJL5VTOYKfr9U8X08UDLdEJAWth1V+EkOe7G4qk7sJIYcB10vB6JwDOTzBFwy38uNJkYiISDwGXBeUGiqCWb1luCUiOWm94ktE/mHA9ZFaAiDDrXop9UKLiIhI6RhwvaCESkKwwg7DLREFkxYu5LRwjERKw4AL8WNfAxECAzH+VgnBnIiCo6ysDImJiYiIiEBaWhoOHjzo1XZbt26FTqfD7NmzA9tAEdh3EZG/GHC70fKVNqu3yqHlv0MSb9u2bTCbzSgqKkJtbS2Sk5ORlZWFCxcueNzuzJkzePbZZzFt2rSgtJN9DBEFCwNuL7pXEuSo3voSdsRWQHjiIVKvkpIS5OXlITc3F+PHj0d5eTn69++PiooKt9vYbDY8/PDDWLZsGUaMGBG0tmqtr+HFKpE8GHBDEMOtuvGESGJYrVbU1NQgMzPTsSwsLAyZmZmorq52u93y5csRExODxx9/vNd9dHZ2orW11ekRKByeQERSYMDtQonBItBtYrglUrdLly7BZrMhNjbWaXlsbCwsFovLbfbv349NmzZhw4YNXu2juLgYUVFRjofRaBTdzut9TddPrJTY5xJRaNB8wPU0PEAJwxPEElP90Fq4VdoXabjCEz4FWltbGx599FFs2LAB0dHRXm1TWFiIlpYWx6OxsdGr7WbE5/vTVI9aRhp6PJSG72ci+fSRuwFyUnrgYedIRL2Jjo5GeHg4mpubnZY3NzcjLi6ux/onT57EmTNnMGvWLMcyu90OAOjTpw/q6+sxcuRIp20MBgMMBvEBUuxwA2/XdxdmW0YaFN+vE1FwaL6C6y0lVge6Y/VW3XhBQ77Q6/VITU1FZWWlY5ndbkdlZSVMJlOP9ceOHYvDhw+jrq7O8fjhD3+Ie+65B3V1dT4NPwgmtfRdfD+Tt3Y3lcndhJCk2Qpu5OlOoE+E43n3zigYNzp4qjSI7RwZbom0y2w2Y968eZg8eTKmTJmC0tJSdHR0IDc3FwCQk5ODoUOHori4GBEREZgwYYLT9oMGDQKAHsuVxpu+i1VcIgI0HHC1iuGWKPRkZ2fj4sWLWLJkCSwWC1JSUrBnzx7HjWcNDQ0IC5PvAztX/Y7YogL7LiISgwHXC8HuWANVveUJQrn4cSb5q6CgAAUFBS5/VlVV5XHbzZs3S98gDeP7mUh+HIPrgtzDE8TgnJFEFOp4cU5EYjHgwvPVdqh0rKFyHKGI1R4i99TWd/H9TKQMigi4ZWVlSExMREREBNLS0nDw4EGvttu6dSt0Oh1mz54d2AZ60DZch7bhOlHbSHVzGYcmEFEo4CdRRCQ12QPutm3bYDabUVRUhNraWiQnJyMrKwsXLlzwuN2ZM2fw7LPPYtq0aZK2R0xHKzbYyoHh9hql3lXNag9pidj3oa/9F/s9IpI94JaUlCAvLw+5ubkYP348ysvL0b9/f1RUVLjdxmaz4eGHH8ayZcswYsQIv/bv6/CEQITbQFRviYjk5M0MCqEiVI+LSI1kDbhWqxU1NTXIzMx0LAsLC0NmZiaqq6vdbrd8+XLExMTg8ccfD0Yze/An3EpRSeTQhNDAkyGR+/6M/RcR+UPWacIuXboEm83mmKvxutjYWBw7dszlNvv378emTZtQV1fn1T46OzvR2XkjVLa2tjr+raSAIXVbAn1yuB7yB54VArofIiI1UNL5hIhUNg9uW1sbHn30UWzYsAHR0dFebVNcXIxly5Z5tW7XSoK7gKiG6m0gwq2743a1nKGXiO69p9jnbVm9JSJ/yRpwo6OjER4ejubmZqflzc3NiIuL67H+yZMncebMGcyaNcuxzG63AwD69OmD+vp6jBw50mmbwsJCmM1mx/PW1lYYjUboj58HwvRSHo7PlHrl72uYbxuuY8jthVL/z4mIiEKBrAFXr9cjNTUVlZWVjqm+7HY7KisrXX4jz9ixY3H48GGnZYsWLUJbWxteeuklGI3GHtsYDAYYDNJUA+SeNSFY1VspjpNDGIi0rWWkweWnVl0v7lz1aVL0YcGeNYUXrOSr3U1lcjchZMk+i4LZbMaGDRvwxhtv4OjRo3jqqafQ0dGB3NxcAEBOTg4KCwsBABEREZgwYYLTY9CgQRg4cCAmTJgAvd73iqw3wxP84a7D9bZjVFO4DeTr+UJpU4TxZEhawaEGRPTRRx9h1qxZSEhIgE6nw44dOzyuX1VVBZ1O1+NhsVhE7Vf2MbjZ2dm4ePEilixZAovFgpSUFOzZs8dx41lDQwPCwmTP4YoIap74eyJR+vERkbp5c6GpxkDMC1Yizzo6OpCcnIx///d/x49+9COvt6uvr0dkZKTjeUxMjKj9yh5wAaCgoMDlkATgWpL3ZPPmzdI3SGKBrt4qPdxyTO4NPBkSXcO5vIm0YcaMGZgxY4bo7WJiYjBo0CCf9yt/aVQBehueEKrVTV++ZtiffRERucPqLZF6tLa2Oj26TscqlZSUFMTHx+O+++7D3/72N9HbK6KCq3a+VCflrt4ycBIREWnXwEYB4Xpx+cVmvbZ+95v6i4qKsHTpUknaFR8fj/LyckyePBmdnZ3YuHEjMjIy8Mknn+DOO+/0+nUYcHvhbxD05wYnhtvQwmoPaV0ovQdC6ViIxGpsbHQaHyvVbFUAMGbMGIwZM8bxPD09HSdPnsS6devw29/+1uvX0XzADfTsCa7I1THKHWw5FpeIgJ4X72ocnkCkZZGRkU4BN9CmTJmC/fv3i9qGY3D95CmwKal6K3e4lYtSpghjtYeIiMg3dXV1iI+PF7WNpiu4vYXIQIRCKYIOwy0RhQo1Vm95wUrkvfb2dpw4ccLx/PTp06irq8OQIUNw6623orCwEOfOncNvfvMbAEBpaSluu+023HHHHbhy5Qo2btyIv/zlL3jvvfdE7VfTAbcrJXWynoK32sMthykQaY9SPkkhouD77LPPcM899ziem81mAMC8efOwefNmNDU1oaGhwfFzq9WKn/3sZzh37hz69++PpKQkvP/++06v4Q3NBlzr6ISAHryvHbqUc0MqLdxqFas9pDWTnyhx+7NAz38b6DDN9zOROBkZGRAE94Wt7t9n8Itf/AK/+MUv/N4vx+C64U04FFuJ9LdjFFO9ZbglIiVx1f8p6ZMzIgotDLiQvpMNRPWW4VadWO0hCh18P5OUdjeVyd2EkMaA64LSbi4LtXAbrDZy3B9pSVlZGRITExEREYG0tDQcPHjQ7bpvv/02Jk+ejEGDBuGmm25CSkqKqPkliYiUTvMB19fqrbvhCXKOvVVDuCUi6W3btg1msxlFRUWora1FcnIysrKycOHCBZfrDxkyBM8//zyqq6vx97//Hbm5ucjNzcXevXsD2k455h0nIm3SfMANht6qt1IMTWC4VR5+nEnBUlJSgry8POTm5mL8+PEoLy9H//79UVFR4XL9jIwM/Mu//AvGjRuHkSNHYv78+UhKShI9kbqW8P1MpC6aDriuwqOSgiLDLRH1xmq1oqamBpmZmY5lYWFhyMzMRHV1da/bC4KAyspK1NfX4+6773a5TmdnJ1pbW50e3uJQISKSg6YDrq/EDE/wp3rrDYZbZWK1h4Ll0qVLsNlsiI2NdVoeGxsLi8XidruWlhYMGDAAer0eM2fOxMsvv4z77rvP5brFxcWIiopyPIxGo6g2dn8/BGp4QqDCNN/PROqj2YDbepuyx395cwJguHWPVSMizwYOHIi6ujp8+umneOGFF2A2m1FVVeVy3cLCQrS0tDgejY2NovcX6PlviYi60uwXPbjiT2CUsnrLcEtE3oqOjkZ4eDiam5udljc3NyMuLs7tdmFhYbj99tsBACkpKTh69CiKi4uRkZHRY12DwQCDQdlFgUBh9ZZInTRbwfWVVF8zy3AbunhCpGDS6/VITU1FZWWlY5ndbkdlZSVMJpPXr2O329HZGfhPPjh7AhEFAyu4/xTs6q2vGG6JqDuz2Yx58+Zh8uTJmDJlCkpLS9HR0YHc3FwAQE5ODoYOHYri4mIA18bUTp48GSNHjkRnZyd27dqF3/72t1i/fr2ch0FEJBkGXBn4Wr1luCUiV7Kzs3Hx4kUsWbIEFosFKSkp2LNnj+PGs4aGBoSF3fjArqOjAz/96U/x5Zdfol+/fhg7dizefPNNZGdny3UIisRPY4jUiwEX3gdHb4cn+NIpMtyGBp4QSS4FBQUoKChw+bPuN4+tXLkSK1euDEg7Jj9R4vRcf+xLx0V9IIcn8MZSIupK82Nw/Q2OYjtVX+4kZrgVhyc6Ivl0ff+p+b3Ii1UiddN8wPWWFNVbX4YmMNyqB0+IRETkjd1NZXI3IeRpeohCsKu3rgTrjuLOW61uf2Zo0AelDUQU+tz1aZw9gYiCSdMB11uBrN564m8A9xRqXa3HoEtEUmgbrnP0m576vu59nFTTMPqLn8YQqZ9mhyi0GYNXvQ3m0ITOW62Ohy/bqp1cY/54QiQSx1UfxyFZRCQVzQZcqYkNOFKHW19DravXISLyVdf+K1gXflJe2PJilSg0MOD2wt+PzFxVbwMRbqXEkEtEUure57FSS0SBxoDrg+7VAqmu+MV2+lJVbd29NnmHFR+intxVVXvr5xh+iUgKDLgeBLt6KwYDKBEpnS8318qJF6tEoYMBVyRvq7eBHJoQrHCrthAtxw1mPCES9eSuOMDqLBEFC6cJUwBvO321BU4iIs5/S0RyYAXXDVcViEBUbxluiYh8J9UnN/w0hii0MOAGgZrDbSD2rZTJ3P3FEyKRsxnx+Y5/d39/cHgCEQUTA64LUldv/cHKrXfk+oIHIrrh4syRjn+r6QYzXqwShR4G3ADzp3rLcEtEana9/2P1luiG3U1lcjdBExhwveBr9ZbhNnSx4kPU0/W+jZ+oENF1H330EWbNmoWEhATodDrs2LGj122qqqpw5513wmAw4Pbbb8fmzZtF75cBtxtfx4d6+3Fcb+E2kF/eQEREznixShRYHR0dSE5ORlmZd5Xr06dPY+bMmbjnnntQV1eHp59+Gk888QT27t0rar+cJqwXvn5rmavqrTfhlsQLdrWIJ0Qi7wVyeAIrxUTKN2PGDMyYMcPr9cvLy3Hbbbdh7dq1AIBx48Zh//79WLduHbKysrx+HVZwJeDt0ARPGG6JKJRw/lui0NXa2ur06OyU7mKzuroamZmZTsuysrJQXV0t6nVYwe2i+/AEX6u3rniqYmgp3IbKFGFEpH78NIa0LPJ0J/r0EfcJy7ffXstFRqPRaXlRURGWLl0qSbssFgtiY2OdlsXGxqK1tRXffPMN+vXr59XrMOD6yZvqLcNt4HB4ApGydL2I5ewJRKGpsbERkZGRjucGg/I+sWHA/afeKouugg3DLREREWlNZGSkU8CVUlxcHJqbm52WNTc3IzIy0uvqLcAxuG4FujLIcEtE5Dt/+2h+GkOkTCaTCZWVlU7L9u3bB5PJJOp1GHC9EIjqrRZJPf6WwxOIbigrK0NiYiIiIiKQlpaGgwcPul13w4YNmDZtGgYPHozBgwcjMzPT4/piRJ3s5A1mROTQ3t6Ouro61NXVAbg2DVhdXR0aGhoAAIWFhcjJyXGs/5//+Z84deoUfvGLX+DYsWN45ZVX8Ic//AHPPPOMqP0y4KL3m8t6w6EJRCSnbdu2wWw2o6ioCLW1tUhOTkZWVhYuXLjgcv2qqirMnTsXH3zwAaqrq2E0GnH//ffj3LlzkrWJF/hEBACfffYZJk2ahEmTJgEAzGYzJk2ahCVLlgAAmpqaHGEXAG677Tbs3LkT+/btQ3JyMtauXYuNGzeKmiIM4BjcXnlTve2O4Tb0sHpLSlZSUoK8vDzk5uYCuDaP5M6dO1FRUYGFCxf2WP93v/ud0/ONGzfij3/8IyorK50qKaGK72ei4MnIyIAguP8E19W3lGVkZODQoUN+7VfzAVds9dafOW+lCLeJwy46/n3my1v8fj214gTvRNdYrVbU1NSgsLDQsSwsLAyZmZlezxv59ddf4+rVqxgyZIjLn3d2djrNc9na2tpjnTEr1rl9/a59n6FB71WbiIj8oemA68vMCb1xV731Ndx2DbSefqb0sMv5b4kC49KlS7DZbC7njTx27JhXr7FgwQIkJCT0mFz9uuLiYixbtsyn9ont+7zpK3iBS0S90ewY3IGNPTtRf6u3UoXbxGEXHQ8x21Bg8ONMCmWrV6/G1q1b8c477yAiIsLlOoWFhWhpaXE8Ghsb3b5ey0iDoodp8f1MctrdVCZ3EzRD0xVcT3rrBAMZbn2VOOyi4iu5UmD1huiG6OhohIeHu5w3Mi4uzuO2a9aswerVq/H+++8jKSnJ7XoGg8GnidzlDrNEpF2KqOAqYXobsdVbb4jp3MVWbINFivFyah6ewGoPKZ1er0dqaqrTvJF2ux2VlZUe54381a9+hRUrVmDPnj2YPHmyJG1R83udiEKL7AFXidPbSFG99TbcSh1slRiSiSiwzGYzNmzYgDfeeANHjx7FU089hY6ODsesCjk5OU43of3yl7/E4sWLUVFRgcTERFgsFlgsFrS3t0vWpkBVb/35BIcXrETaIXvA7Tq9zfjx41FeXo7+/fujoqLC5fq/+93v8NOf/hQpKSkYO3YsNm7c6KhWBErX6q23QxO8wTAqHocnEPWUnZ2NNWvWYMmSJUhJSUFdXR327NnjuPGsoaEBTU1NjvXXr18Pq9WKf/u3f0N8fLzjsWbNGr/bwvlviUgJZB2Dq5TpbbqGpu5X+IEamqCVcKvmjyxZ7SE1KSgoQEFBgcufVVVVOT0/c+ZMQNoQdbITbcNd36gmN76fibRF1gqup+ltLBaLV6/hzfQ2UVFRjofRaPS5vVIMTVDqWFtXlDZfJau3RMp2vY/kzWVEJDfZhyj4w9/pbaJOdnpdvZUq3JI6sNpDFHxq/sSHiJRF1iEKck5vE3m6E+jj20dpDLfe4cmKiKTk66c4vGAl0h5ZK7hKmt5GTPW2u1AMt1oensCTIZFveIMZESmF7F/0YDabMW/ePEyePBlTpkxBaWlpj+lthg4diuLiYgDXprdZsmQJtmzZ4pjeBgAGDBiAAQMG+NQGf4YmhGK4lQKrt0TaxPG3RKQEsgfc7OxsXLx4EUuWLIHFYkFKSkqP6W3Cwm4UmrtOb9NVUVERli5dGtC2iq1OaDXcSo3VWyLyFd/TRNoke8AF5J3eJlBDE9QcbpU2PIGIlG3MinUB3wdnUSEiMRQRcOUi5sre26EJag62UpFyeAKrt0TawGFNRCQlVU8T5g/98fM9lrmr3jLcEhF5NvCsoLjxt7xoJSXZ3VQmdxM0RbMB1xM1h9szX97i1/b+Dk9g9ZZImziDAhEpCQPuP/nylbzXKSXcEhEpnS8X0Rx/S0RiaXoM7nX+DE0IpXDL6i0R+atrn+jvJ0r+4PuaSNtYwfWSGsKtnCcTtYZbIpKGP+NveYMZEUlN8wHXm+qt0m6cCAStTg3GKg+R/wIZUHnBS0S+0HTA9XZogitKq97KSa3VW4ZbImkppV/ke5uINBtwraMTHP8OhXG3/gxP0Gr1loiIiEKTZgOuN9QSbv3BG8uIyF/+TBEm9fAGvreJCOAsCqLG3YoJt/fFHXP7s32WsV6/jjfkurlMreGWiKTnTf8o9oKa/QIR+UrTFVwx4269Dbf3xR3zGG6vr6MEWh2awAoPkbS0cCMuEamLZgNu623ej7v1Jtx6E2wDwdfqLYcmEFGo4fubSJnKysqQmJiIiIgIpKWl4eDBg27X3bx5M3Q6ndMjIiJC9D41P0ShO7HhVinV2GBSa7glImXx1JewbyAKDdu2bYPZbEZ5eTnS0tJQWlqKrKws1NfXIyYmxuU2kZGRqK+vdzzX6cSP89dsBfc6T+NulR5u5ajeqjncsrpDJL177ymWuwkA+P4mUqqSkhLk5eUhNzcX48ePR3l5Ofr374+Kigq32+h0OsTFxTkesbGxover2YDbZtT5HG7lGo7QldrDbbDx5EcUOIG4wYyIlKu1tdXp0dnpumBltVpRU1ODzMxMx7KwsDBkZmaiurra7eu3t7dj+PDhMBqNePDBB3HkyBHRbdRswO1KbLhVKyWdYDjulkhaYsa4HTlyBP/6r/+KxMRE6HQ6lJaW+rzfrjfrSkVs/8D3OCnd7qYyuZvQg/74eeiPfSnucfw8AMBoNCIqKsrxKC52/UnOpUuXYLPZelRgY2NjYbFYXG4zZswYVFRU4N1338Wbb74Ju92O9PR0fPmluPe55sfgqjHc+lK91epNZURaIHaM29dff40RI0Zgzpw5eOaZZ/zad9twHQb4sJ2aPxEi0rrGxkZERkY6nhsM0l3omkwmmEwmx/P09HSMGzcOr776KlasWOH162i6guttuJV6SII/8+DKMeetmsMtKzukBWLHuH3nO9/Biy++iJ/85CeSnpi6kmt+biIKvMjISKeHu34kOjoa4eHhaG5udlre3NyMuLg4r/bVt29fTJo0CSdOnBDVRs0GXKvR+3CrFGofd8twSyQ9X8e4ScWbOXDF9EEcnkAUOvR6PVJTU1FZWelYZrfbUVlZ6VSl9cRms+Hw4cOIj48XtW/ND1EAGG7dYbglUj5PY9yOHZOmD+vs7HS6iaS1tdWv1+PwBCLtMJvNmDdvHiZPnowpU6agtLQUHR0dyM3NBQDk5ORg6NChjnG8y5cvx1133YXbb78dly9fxosvvoizZ8/iiSeeELVfzQfcYIdbqb+m1xOGWyKSQnFxMZYtW+byZ2K+wrw3rN4ShZ7s7GxcvHgRS5YsgcViQUpKCvbs2eO4KG9oaEBY2I0BBV999RXy8vJgsVgwePBgpKam4sCBAxg/fryo/Wo24N6acAl9bnI9ZkRp4TbYN5Ux3BKphxRj3HpTWFgIs9nseN7a2gqj0eh4fr3PDOYFPBGpR0FBAQoKClz+rKqqyun5unXrsG7dOr/3qdkxuK4oYX7b7hhuvcdwS1okxRi33hgMhh43lXjStd/q3h9J1cfw/U5Enmi2gttdoIOtL5UNhlvv8WRHWiZ2jJvVasX//u//Ov597tw51NXVYcCAAbj99ttF7VuqvpPTBxKRlDQfcINRsQ1GuNVqsAUYbonEjnE7f/48Jk2a5Hi+Zs0arFmzBtOnT+/xcaES8T1PRL3RbMC9J/Y4Igb0lbsZLjHceo8nOqJrxIxxS0xMhCAEfiYDb4cnsHpLRFLjGNwAE1u9Zbj1HsMtkXIE6wYzvu+JyBuareAGg1LDrdRzUDLcEtF17voxKaq3fN8TkbcYcANETLhl1dZ7PMERKceELevwH9PkbgURUU8coiAzhlvvMdwSKcutCZfc/sybvorVWyIKFE1XcP8tstbjz7e33unT63pbvVXjkAQGWyK67p7Y4wB6v1nX3z6IfQCp2e6mMrmboEmaDbizB/4PAlHADkS41WqwBXhiI1ITMf0aZ04gokDSbMANBG/CLau23mGwJVKvrv0Wq7dEJAcGXIlIHW61GmwBntCI1EbsjDHe9i3sC4jIVwy4bogZfytluGWwJSI16616y6EJRBQMDLh+UELVlsGWiJRA7PCr3rBfICJ/MOC64E31Vu6qLYMtEclt9sD/wR5hssuf+VO9Zd9ARP5iwO1GinCrhmAr58eEPHkRhSZPfRrDLREFEwNuF8EKt1oMtjxpEWmHr/0V+wkikgoDLpQbbNUeagGesIi0oLf+z5s+iH0FEUlJswF3R1syIoTev4GHwVY8nqiItON6H3m9n+vefzHcEpEcNBtwveEp3EodbBlqiSjUMNwSkVwYcF1QW7BlqCUiObmr3vaGfQgRBQoD7j/5OxQhmKGWgZaIlOLMl7fA4GJ5b/0U+xMiCiTNBtwPmkejT7urbtmZVMFWjaGWJyAi8mRHW7Lj3137OIZbIpKbZgOuJ0oItQy0RKQG3ftDT30X+xgiChYGXEg3rlZNgZYnGiLy1wfNowHc6PsYbomc7W4qk7sJmqXZgNtwPhph/SLc/jyQVdpghlmeVIgoUBrOR6PfP//trl9jH0REctBswO0qFMIsTyJEJAdP/SL7JSKSi2YDrr5Rj/AI98FWiWGWJwsiUhJ947U+tHvfx76KiOSm2YA7sFFAuF45IZYnBCJSm4GNAqIab/SJ7MeIyJWysjK8+OKLsFgsSE5Oxssvv4wpU6a4Xf+tt97C4sWLcebMGYwaNQq//OUv8cADD4jap2YDrjuBCLDs9IkoFEWe7gT6RLCPIyK3tm3bBrPZjPLycqSlpaG0tBRZWVmor69HTExMj/UPHDiAuXPnori4GD/4wQ+wZcsWzJ49G7W1tZgwYYLX+9VswI083Yk+fXSitmEnTkTkjP0iEXlSUlKCvLw85ObmAgDKy8uxc+dOVFRUYOHChT3Wf+mll/D9738fP//5zwEAK1aswL59+/DrX/8a5eXlXu9XcwFXEK4NSwg7dgZhYd5/+xgAfBuIBhFRyPnWbgVwo78JNV370W9F9qNEWtLa2hrw1/a1n/lWsAJ2H7ZBz+MyGAwwGHp+eZbVakVNTQ0KCwsdy8LCwpCZmYnq6mqX+6iurobZbHZalpWVhR07dohqq+YCbltbGwCg6tJvZG4JEYW6trY2REVFyd0MybEfJfJOVNTGgO9DbD+j1+sRFxeHKotv798BAwbAaDQ6LSsqKsLSpUt7rHvp0iXYbDbExsY6LY+NjcWxY8dcvr7FYnG5vsViEdVOzQXchIQENDY2YuDAgdDpxA1RUIvW1lYYjUY0NjYiMjJS7uZIjsenfqF+jIIgoK2tDQkJCXI3JSBCrR8N9b9HHp+6uTs+X/uZiIgInD59Glar1af2CILQ433vqnorN80F3LCwMAwbNkzuZgRFZGRkSL7Zr+PxqV8oH2MoVm6vC9V+NJT/HgEen9q5Oj5f+5mIiAhERLj/siupREdHIzw8HM3NzU7Lm5ubERcX53KbuLg4Ueu7EyauqUREREREvdPr9UhNTUVlZaVjmd1uR2VlJUwmk8ttTCaT0/oAsG/fPrfru6O5Ci4RERERBYfZbMa8efMwefJkTJkyBaWlpejo6HDMqpCTk4OhQ4eiuLgYADB//nxMnz4da9euxcyZM7F161Z89tlneO2110TtlwE3BBkMBhQVFSlyTIwUeHzqp4VjJPUI9b9HHp+6qf34srOzcfHiRSxZsgQWiwUpKSnYs2eP40ayhoYGhIXdGFCQnp6OLVu2YNGiRXjuuecwatQo7NixQ9QcuACgE0J1HhsiIiIi0iSOwSUiIiKikMKAS0REREQhhQGXiIiIiEIKAy4RERERhRQGXJUqKytDYmIiIiIikJaWhoMHD7pdd8OGDZg2bRoGDx6MwYMHIzMz0+P6SiDm+LraunUrdDodZs+eHdgG+kns8V2+fBn5+fmIj4+HwWDA6NGjsWvXriC11jdij7G0tBRjxoxBv379YDQa8cwzz+DKlStBai2FqqVLl0Kn0zk9xo4d63Gbt956C2PHjkVERAQmTpyo6PdaYmJij+PT6XTIz893uf7mzZt7rBuMCf+99dFHH2HWrFlISEiATqfDjh07nH4uCAKWLFmC+Ph49OvXD5mZmfjiiy96fV1fzylS83R8V69exYIFCzBx4kTcdNNNSEhIQE5ODs6fP+/xNX35G9cCBlwV2rZtG8xmM4qKilBbW4vk5GRkZWXhwoULLtevqqrC3Llz8cEHH6C6uhpGoxH3338/zp07F+SWe0fs8V135swZPPvss5g2bVqQWuobscdntVpx33334cyZM9i+fTvq6+uxYcMGDB06NMgt957YY9yyZQsWLlyIoqIiHD16FJs2bcK2bdvw3HPPBbnlFIruuOMONDU1OR779+93u+6BAwcwd+5cPP744zh06BBmz56N2bNn4/PPPw9ii7336aefOh3bvn37AABz5sxxu01kZKTTNmfPng1Wc3vV0dGB5ORklJWVufz5r371K/z3f/83ysvL8cknn+Cmm25CVlaWx4thX88pgeDp+L7++mvU1tZi8eLFqK2txdtvv436+nr88Ic/7PV1xfyNa4ZAqjNlyhQhPz/f8dxmswkJCQlCcXGxV9t/++23wsCBA4U33ngjUE30iy/H9+233wrp6enCxo0bhXnz5gkPPvhgEFrqG7HHt379emHEiBGC1WoNVhP9JvYY8/Pzhe9973tOy8xmszB16tSAtpNCX1FRkZCcnOz1+j/+8Y+FmTNnOi1LS0sT/uM//kPilgXG/PnzhZEjRwp2u93lz19//XUhKioquI3yEQDhnXfecTy32+1CXFyc8OKLLzqWXb58WTAYDMLvf/97t6/j7zkzULofnysHDx4UAAhnz551u47Yv3GtYAVXZaxWK2pqapCZmelYFhYWhszMTFRXV3v1Gl9//TWuXr2KIUOGBKqZPvP1+JYvX46YmBg8/vjjwWimz3w5vj/96U8wmUzIz89HbGwsJkyYgFWrVsFmswWr2aL4cozp6emoqalxfGx46tQp7Nq1Cw888EBQ2kyh7YsvvkBCQgJGjBiBhx9+GA0NDW7Xra6udvrbBYCsrCyv+1c5Wa1WvPnmm/j3f/936HQ6t+u1t7dj+PDhMBqNePDBB3HkyJEgttJ3p0+fhsVicfr/iYqKQlpamtv/HynOmXJqaWmBTqfDoEGDPK4n5m9cKxhwVebSpUuw2WyObwC5LjY2FhaLxavXWLBgARISEnp04krgy/Ht378fmzZtwoYNG4LRRL/4cnynTp3C9u3bYbPZsGvXLixevBhr167FypUrg9Fk0Xw5xoceegjLly/Hd7/7XfTt2xcjR45ERkYGhyiQ39LS0rB582bs2bMH69evx+nTpzFt2jS0tbW5XN9isfjVv8ppx44duHz5Mh577DG364wZMwYVFRV499138eabb8JutyM9PR1ffvll8Brqo+v/B2L+f6Q4Z8rlypUrWLBgAebOnYvIyEi364n9G9cKflWvxqxevRpbt25FVVWVom4s8FVbWxseffRRbNiwAdHR0XI3JyDsdjtiYmLw2muvITw8HKmpqTh37hxefPFFFBUVyd08SVRVVWHVqlV45ZVXkJaWhhMnTmD+/PlYsWIFFi9eLHfzSMVmzJjh+HdSUhLS0tIwfPhw/OEPf1D8Jz5ibdq0CTNmzEBCQoLbdUwmE0wmk+N5eno6xo0bh1dffRUrVqwIRjPJC1evXsWPf/xjCIKA9evXe1xXS3/jYjDgqkx0dDTCw8PR3NzstLy5uRlxcXEet12zZg1Wr16N999/H0lJSYFsps/EHt/Jkydx5swZzJo1y7HMbrcDAPr06YP6+nqMHDkysI0WwZf/v/j4ePTt2xfh4eGOZePGjYPFYoHVaoVerw9om8Xy5RgXL16MRx99FE888QQAYOLEiejo6MCTTz6J559/3ul7yon8MWjQIIwePRonTpxw+fO4uDif+le5nT17Fu+//z7efvttUdv17dsXkyZNcvv7UJLr/wfNzc2Ij493LG9ubkZKSorLbfw5Z8rlerg9e/Ys/vKXv3is3rrS29+4VvCsoTJ6vR6pqamorKx0LLPb7aisrHS6Ku/uV7/6FVasWIE9e/Zg8uTJwWiqT8Qe39ixY3H48GHU1dU5Hj/84Q9xzz33oK6uDkajMZjN75Uv/39Tp07FiRMnHMEdAI4fP474+HjFhVvAt2P8+uuve4TY64FeEITANZY0p729HSdPnnQKSF2ZTCanv10A2Ldvn8f+VQlef/11xMTEYObMmaK2s9lsOHz4sNvfh5LcdtttiIuLc/r/aW1txSeffOL2/8fXc6ZcrofbL774Au+//z5uvvlm0a/R29+4Zsh9lxuJt3XrVsFgMAibN28W/vd//1d48sknhUGDBgkWi0UQBEF49NFHhYULFzrWX716taDX64Xt27cLTU1NjkdbW5tch+CR2OPrTumzKIg9voaGBmHgwIFCQUGBUF9fL/z5z38WYmJihJUrV8p1CL0Se4xFRUXCwIEDhd///vfCqVOnhPfee08YOXKk8OMf/1iuQ6AQ8bOf/UyoqqoSTp8+Lfztb38TMjMzhejoaOHChQuCIPT8W/zb3/4m9OnTR1izZo1w9OhRoaioSOjbt69w+PBhuQ6hVzabTbj11luFBQsW9PhZ9+NbtmyZsHfvXuHkyZNCTU2N8JOf/ESIiIgQjhw5Eswmu9XW1iYcOnRIOHTokABAKCkpEQ4dOuSYRWD16tXCoEGDhHfffVf4+9//Ljz44IPCbbfdJnzzzTeO1/je974nvPzyy47nvfVHSjk+q9Uq/PCHPxSGDRsm1NXVOZ2vOzs73R5fb3/jWsWAq1Ivv/yycOuttwp6vV6YMmWK8PHHHzt+Nn36dGHevHmO58OHDxcA9HgUFRUFv+FeEnN83Sk94AqC+OM7cOCAkJaWJhgMBmHEiBHCCy+8IHz77bdBbrU4Yo7x6tWrwtKlS4WRI0cKERERgtFoFH76058KX331VfAbTiElOztbiI+PF/R6vTB06FAhOztbOHHihOPnrt5vf/jDH4TRo0cLer1euOOOO4SdO3cGudXi7N27VwAg1NfX9/hZ9+N7+umnHe/L2NhY4YEHHhBqa2uD2FrPPvjgA5fnq+vHYLfbhcWLFwuxsbGCwWAQ7r333h7HPXz48B7nN0/9UTB5Or7Tp0+7/BkA4YMPPnC8Rvfj6+1vXKt0gsDP/4iIiIgodHAMLhERERGFFAZcIiIiIgopDLhEREREFFIYcImIiIgopDDgEhEREVFIYcAlIiIiopDCgEtEREREIYUBl4iIiIhCCgMuEREREYUUBlwiIiIiCikMuERuJCYmorS01GlZSkoKli5dKkt7iIj8lZGRgYKCAhQUFCAqKgrR0dFYvHgxBEGQu2lEkmLAJSIi0pA33ngDffr0wcGDB/HSSy+hpKQEGzdulLtZRJLqI3cDiIiIKHiMRiPWrVsHnU6HMWPG4PDhw1i3bh3y8vLkbhqRZFjBJSIi0pC77roLOp3O8dxkMuGLL76AzWaTsVVE0mLAJXIjLCysx7i0q1evytQaIiIi8hYDLpEbt9xyC5qamhzPW1tbcfr0aRlbRETkv08++cTp+ccff4xRo0YhPDxcphYRSY8Bl8iN733ve/jtb3+Lv/71rzh8+DDmzZvHEwARqV5DQwPMZjPq6+vx+9//Hi+//DLmz58vd7OIJMWbzIjcKCwsxOnTp/GDH/wAUVFRWLFiBSu4RKR6OTk5+OabbzBlyhSEh4dj/vz5ePLJJ+VuFpGkdAInvyMiItKEjIwMpKSk9JjjmyjUcIgCEREREYUUBlwiIiIiCikcokBEREREIYUVXCIiIiIKKQy4RERERBRSGHCJiIiIKKQw4BIRERFRSGHAJSIiIqKQwoBLRERERCGFAZeIiIiIQgoDLhERERGFFAZcIiIiIgop/x8f+Nw1VS149AAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# (2, 2, 0) and (7, -3, 1) modes\n", "specific_modes = [(2, 2, 0), (7, -3, 1)]\n", "\n", "# notice this returns a dictionary with keys as the mode tuple and values as the mode values at all trajectory points\n", "specific_teuk_modes = kerr_amp_spline(\n", " a, p_map, e_map, xI_map, specific_modes=specific_modes\n", ")\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", "# make sure they are the same\n", "print(\n", " np.allclose(\n", " specific_teuk_modes[(2, 2, 0)].reshape(uv.size, wv.size),\n", " teuk_modes[:, :, inds[0]],\n", " )\n", ")\n", "\n", "# to check -m modes we need to take the conjugate\n", "print(\n", " np.allclose(\n", " specific_teuk_modes[(7, -3, 1)].reshape(uv.size, wv.size),\n", " np.conj(teuk_modes[:, :, inds[1]]),\n", " )\n", ")\n", "\n", "# look at the contours of the (7,-3, 1) mode\n", "\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(specific_teuk_modes[(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_map.reshape(uv.size, wv.size),\n", " e_map.reshape(uv.size, wv.size),\n", " np.abs(specific_teuk_modes[(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": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total trajectory points: 2500\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAG2CAYAAACkgiamAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAO9ZJREFUeJzt3Xt0FPX9//FXLiTLxQQkJTdDAqJiRIhNSAxewG9TcvqlIp62Rg8lnByLVUTRVQqxkgBaA0ppLFKiFKTF2tB+qy1f9RsvkbSlBFITcgSUCFYIIglQJYEgCWbn9we/bFlzYbPJ7s7uPh/nzNH97Gdm3h+XZV5+ZnYmyDAMQwAAACYR7O0CAAAALkQ4AQAApkI4AQAApkI4AQAApkI4AQAApkI4AQAApkI4AQAApkI4AQAApkI4AQAApkI4AQAApmKKcLJmzRolJSXJYrEoIyNDVVVV3fadMmWKgoKCOi3Tpk3zYMUAAPiH3hyDN27c2On4a7FYOvX78MMPNX36dEVGRmrw4MGaOHGi6uvrna7J6+Fk8+bNslqtKiwsVE1NjSZMmKDs7GwdO3asy/6vvPKKjh49al/27NmjkJAQ/eAHP/Bw5QAA+LbeHoMlKSIiwuE4fOjQIYf3P/74Y914440aO3asKioq9P7772vx4sVdhpjuBHn7wX8ZGRmaOHGinnvuOUmSzWZTQkKCHnjgAS1atOii6xcXF6ugoEBHjx7V4MGD3V0uAAB+o7fH4I0bN+qhhx7SyZMnu93mnXfeqQEDBmjTpk0u1xXq8pr9oK2tTdXV1crPz7e3BQcHKysrS5WVlU5tY/369brzzju7DSatra1qbW21v7bZbPr88881fPhwBQUF9W0AAAC/ZhiGTp06pbi4OAUHu+9kw9mzZ9XW1tbn7RiG0enYFh4ervDw8E59XT0Gnz59WomJibLZbPrmN7+pp556Stdcc42k88fY119/XT/5yU+UnZ2tXbt2adSoUcrPz9eMGTN6NRCvOXLkiCHJ2L59u0P7ggULjPT09Iuuv3PnTkOSsXPnzm77FBYWGpJYWFhYWFhcXg4fPtznY153vvzySyPqG8H9UueQIUM6tRUWFna5X1eOwdu3bzd+85vfGLt27TIqKiqM7373u0ZERIT9v8/Ro0cNScagQYOMVatWGbt27TKKioqMoKAgo6Kiwun/Jl6dOemr9evX69prr1V6enq3ffLz82W1Wu2vm5qaNHLkSGVMXqTQUOfPfwFAXw384Ki3S0AvfWVrU0Xji7rkkkvcto+2tjadOG7TmztiNHiI67MzLadtyr6+QYcPH1ZERIS9vatZE1dlZmYqMzPT/nrSpEm6+uqr9fzzz+uJJ56QzWaTJN122216+OGHJUkpKSnavn27SkpKNHnyZKf249VwEhUVpZCQEDU2Njq0NzY2KiYmpsd1W1paVFpaqmXLlvXYr7vprNBQC+EEgMcM3HNECg7zdhlwkScuAxg8JFhDLun7qaOIiAiHcNKdvhyDOwwYMEDXXXedDhw4YN9maGiokpOTHfpdffXV2rZtm5Mj8PKvdcLCwpSamqry8nJ7m81mU3l5uUMy68of//hHtba26oc//KG7ywSAPhm454i3SwA66csxuEN7e7t2796t2NhY+zYnTpyouro6h34fffSREhMTna7N66d1rFarZs+erbS0NKWnp6u4uFgtLS3Ky8uTJOXm5io+Pl5FRUUO661fv14zZszQ8OHDvVE2AFwUoQRm19tj8LJly3T99ddrzJgxOnnypJ555hkdOnRIP/rRj+zbXLBggXJycnTzzTfrlltuUVlZmf73f/9XFRUVTtfl9XCSk5Oj48ePq6CgQA0NDUpJSVFZWZmio6MlSfX19Z2ukK6rq9O2bdv01ltveaNkALgoggl8QW+PwV988YXmzJmjhoYGDRs2TKmpqdq+fbvDaZzbb79dJSUlKioq0oMPPqirrrpKf/rTn3TjjTc6XZfX73Piac3NzYqMjNQN31rCNScA3IJg4j++srXpnaPPq6mpyanrOFzRcVzatieuT9ecnD5l043jPnNrrZ7i9TvEAoA/IZgAfef10zoA4A8IJUD/YeYEAPqIYAL0L2ZOAMBFhBLAPZg5AQAXEEwA92HmBAB6gVACuB8zJwDgJIIJ4BnMnADARRBKAM9i5gQAekAwATyPmRMA6AKhBPAewgkA/H8EEsAcCCcAAh6hBDAXwgmAgEUoAcyJcAIg4BBKAHMjnAAIGIQSwDcQTgD4PUIJ4FsIJwD8FqEE8E2EEwB+h1AC+DbCCQC/QCAB/AfhBIBPI5QA/odwAsDnEEgA/0Y4AeAzCCVAYCCcADA9QgkQWAgnAEyJQAIELsIJAFMhlAAgnAAwBUIJgA6EEwBeQyAB0BXCCQCPI5QA6AnhBIDHEEoAOINwAsCtCCQAeotwAsAtCCUAXEU4AdCvCCUA+opwAqDPCCQA+hPhBIBLCCQA3IVwAsBpBBIAnkA4AdAjAgkATyOcAOiEQALAmwgnAOwIJQDMgHACBDgCCQCzIZwAAYhAAsDMCCdAgCCQAPAVhBPATxFGAPiqYG8XAKD/DNxzxL4AgDPWrFmjpKQkWSwWZWRkqKqqyqn1SktLFRQUpBkzZji0nz59WvPmzdNll12mgQMHKjk5WSUlJb2qiZkTwMcRRAC4avPmzbJarSopKVFGRoaKi4uVnZ2turo6jRgxotv1Dh48qEcffVQ33XRTp/esVqveffddvfTSS0pKStJbb72luXPnKi4uTtOnT3eqLmZOAB/EDAmA/rBq1SrNmTNHeXl59hmOQYMGacOGDd2u097erpkzZ2rp0qUaPXp0p/e3b9+u2bNna8qUKUpKStI999yjCRMmOD0jIxFOAJ9BIAHgjObmZoeltbW1y35tbW2qrq5WVlaWvS04OFhZWVmqrKzsdvvLli3TiBEjdPfdd3f5/qRJk7RlyxYdOXJEhmFo69at+uijjzR16lSnx8BpHcDECCJA4Cg9maHwrwa4vH7r6XOSXlVCQoJDe2FhoZYsWdKp/4kTJ9Te3q7o6GiH9ujoaO3bt6/LfWzbtk3r169XbW1tt3WsXr1a99xzjy677DKFhoYqODhY69at08033+z0WAgngMkQSAD0xeHDhxUREWF/HR4e3i/bPXXqlGbNmqV169YpKiqq236rV6/Wjh07tGXLFiUmJupvf/ub7r//fsXFxTnM0vSEcAJ4GWEEQH+KiIhwCCfdiYqKUkhIiBobGx3aGxsbFRMT06n/xx9/rIMHD+rWW2+1t9lsNklSaGio6urqFBcXp8cee0yvvvqqpk2bJkkaP368amtrtXLlSsIJYGYEEgDeFhYWptTUVJWXl9t/Dmyz2VReXq558+Z16j927Fjt3r3boe3xxx/XqVOn9OyzzyohIUFnz57VuXPnFBzseElrSEiIPcg4g3ACeABhBIAZWa1WzZ49W2lpaUpPT1dxcbFaWlqUl5cnScrNzVV8fLyKiopksVg0btw4h/WHDh0qSfb2sLAwTZ48WQsWLNDAgQOVmJiov/71r/rtb3+rVatWOV2X13+t09ubv5w8eVL333+/YmNjFR4eriuvvFJvvPGGh6oFnMevawCYXU5OjlauXKmCggKlpKSotrZWZWVl9otk6+vrdfTo0V5ts7S0VBMnTtTMmTOVnJys5cuX62c/+5nuvfdep7fh1ZmT3t78pa2tTd/+9rc1YsQI/c///I/i4+N16NAhe3IDvI0gAsDXzJs3r8vTOJJUUVHR47obN27s1BYTE6MXX3yxTzV5NZxcePMXSSopKdHrr7+uDRs2aNGiRZ36b9iwQZ9//rm2b9+uAQPO/9wqKSnJkyUDnRBIAKB/ee20jis3f9myZYsyMzN1//33Kzo6WuPGjdNTTz2l9vb2bvfT2tra6YY0QF9xygYA3MdrMyeu3PzlX//6l959913NnDlTb7zxhg4cOKC5c+fq3LlzKiws7HKdoqIiLV26tN/rR2AhhACA5/jUr3VsNptGjBihF154QSEhIUpNTdWRI0f0zDPPdBtO8vPzZbVa7a+bm5s73T0P6AqBBAC8w2vhpLc3f5Gk2NhYDRgwQCEhIfa2q6++Wg0NDWpra1NYWFindcLDw/vt7njwb4QRADAHr11zcuHNXzp03PwlMzOzy3VuuOEGHThwwOFGLh999JFiY2O7DCbAxXDtCACYj1fvc2K1WrVu3Tr95je/0Ycffqj77ruv081f8vPz7f3vu+8+ff7555o/f74++ugjvf7663rqqad0//33e2sI8DEXhhECCQCYk1evOcnJydHx48dVUFCghoYGpaSkdLr5y4W3wE1ISNCbb76phx9+WOPHj1d8fLzmz5+vhQsXemsI8AGEEADwLV6/ILa3N3/JzMzUjh073FwVfB2BBAB8l9fDCdBfCCQA4B8IJ/BpBBIA8D+EE/gUwggA+D/CCUyPQAIAgYVwAlMikABA4CKcwDQIJAAAiXACLyOQAAC+jnACjyOQAAB6QjiB2xFGAAC9QTiBWxBIAACuIpyg3xBIAAD9gXCCPiGQAAD6G+EEvUYgAQC4E+EETiGQAAA8hXCCLhFGAADeQjiBHYEEAGAGhJMARyABAJgN4SQAEUgAAGZGOAkAhBEAgC8hnPgpAgkAwFcRTvwIgQQA4A8IJz6OQAIA8DeEEx9DGAEA+DvCiQ8gkAAAAgnhxIQIIwCAQEY4MQkCCQAA5xFOvIhAAgBAZ4QTDyOQAADQM8KJBxBIAABwHuHETQgkAAC4hnDSjwgkAAD0HeGkDwgjAAD0P8JJLxFIAABwr2BvF+ALBu45Yl8AAPAna9asUVJSkiwWizIyMlRVVeXUeqWlpQoKCtKMGTMc2g3DUEFBgWJjYzVw4EBlZWVp//79vaqJcNINAgkAwN9t3rxZVqtVhYWFqqmp0YQJE5Sdna1jx471uN7Bgwf16KOP6qabbur03tNPP61f/vKXKikp0c6dOzV48GBlZ2fr7NmzTtdFOLkAgQQAEEhWrVqlOXPmKC8vT8nJySopKdGgQYO0YcOGbtdpb2/XzJkztXTpUo0ePdrhPcMwVFxcrMcff1y33Xabxo8fr9/+9rf67LPP9Oc//9npugI+nBBIAAD+pLm52WFpbW3tsl9bW5uqq6uVlZVlbwsODlZWVpYqKyu73f6yZcs0YsQI3X333Z3e++STT9TQ0OCwzcjISGVkZPS4za8L2AtiB35wVKHBYd4uAwAASdK7n16hkEHhLq/ffuZ8CElISHBoLyws1JIlSzr1P3HihNrb2xUdHe3QHh0drX379nW5j23btmn9+vWqra3t8v2Ghgb7Nr6+zY73nBGw4QQAAH90+PBhRURE2F+Hh7seeC506tQpzZo1S+vWrVNUVFS/bLM7hBMAAPxIRESEQzjpTlRUlEJCQtTY2OjQ3tjYqJiYmE79P/74Yx08eFC33nqrvc1ms0mSQkNDVVdXZ1+vsbFRsbGxDttMSUlxegwBf80JAACBKCwsTKmpqSovL7e32Ww2lZeXKzMzs1P/sWPHavfu3aqtrbUv06dP1y233KLa2lolJCRo1KhRiomJcdhmc3Ozdu7c2eU2u8PMCQAAAcpqtWr27NlKS0tTenq6iouL1dLSory8PElSbm6u4uPjVVRUJIvFonHjxjmsP3ToUElyaH/ooYf05JNP6oorrtCoUaO0ePFixcXFdbofSk8IJwAABKicnBwdP35cBQUFamhoUEpKisrKyuwXtNbX1ys4uHcnWX7yk5+opaVF99xzj06ePKkbb7xRZWVlslgsTm8jyDAMo1d79XHNzc2KjIxUVuyP+bUOAKBHX9na9M7R59XU1OTUdRyu6DguJZf+pM+/1vngzqfdWquncM0JAAAwFcIJAAAwFcIJAAAwFcIJAAAwFcIJAAAwFcIJAAAwFVOEkzVr1igpKUkWi0UZGRmqqqrqtu/GjRsVFBTksPTmt9MAAMDcvB5ONm/eLKvVqsLCQtXU1GjChAnKzs7WsWPHul0nIiJCR48etS+HDh3yYMUAAMCdvB5OVq1apTlz5igvL0/JyckqKSnRoEGDtGHDhm7XCQoKUkxMjH35+qOZAQCA7/JqOGlra1N1dbWysrLsbcHBwcrKylJlZWW3650+fVqJiYlKSEjQbbfdpr1793bbt7W1Vc3NzQ4LAAAwL6+GkxMnTqi9vb3TzEd0dLQaGhq6XOeqq67Shg0b9Je//EUvvfSSbDabJk2apE8//bTL/kVFRYqMjLQvCQkJ/T4OAADQf7x+Wqe3MjMzlZubq5SUFE2ePFmvvPKKvvGNb+j555/vsn9+fr6amprsy+HDhz1cMQAA6A2vPpU4KipKISEhamxsdGhvbGxUTEyMU9sYMGCArrvuOh04cKDL98PDwxUe7vqDlAAAgGd5deYkLCxMqampKi8vt7fZbDaVl5crMzPTqW20t7dr9+7dio2NdVeZAADAg7w6cyJJVqtVs2fPVlpamtLT01VcXKyWlhbl5eVJknJzcxUfH6+ioiJJ0rJly3T99ddrzJgxOnnypJ555hkdOnRIP/rRj7w5DAAA0E+8Hk5ycnJ0/PhxFRQUqKGhQSkpKSorK7NfJFtfX6/g4P9M8HzxxReaM2eOGhoaNGzYMKWmpmr79u1KTk721hAAAEA/CjIMw/B2EZ7U3NysyMhIZcX+WKHBYd4uBwBgYl/Z2vTO0efV1NSkiIgIt+yj47iUXPoThQxy/RrJ9jOt+uDOp91aq6f43K91AACAfyOcAAAAUyGcAAAAUyGcAAAAUyGcAAAAUyGcAAAAUyGcAAAAUyGcAAAAUyGcAAAAUyGcAAAAUyGcAAAAUyGcAAAAUyGcAAAAUyGcAAAAUwn1dgEAALjiy3Hxbt/HV1+dlY66fTf4GsIJAMAtPBEe4J8IJwAQ4AgRMBvCCQD4EIIEAgHhBADcjEAB9A7hBAC6QagAvINwAsAvESwA30U4AWBKhAsgcBFOALgNAQOAKwgnAC6KkAHAk7h9PRCAvhwX36sFgP9as2aNkpKSZLFYlJGRoaqqqm77vvLKK0pLS9PQoUM1ePBgpaSkaNOmTfb3z507p4ULF+raa6/V4MGDFRcXp9zcXH322We9qomZE8APECAAuGLz5s2yWq0qKSlRRkaGiouLlZ2drbq6Oo0YMaJT/0svvVQ//elPNXbsWIWFhem1115TXl6eRowYoezsbJ05c0Y1NTVavHixJkyYoC+++ELz58/X9OnT9d577zldV5BhGEZ/DtTsmpubFRkZqazYHys0OMzb5QBdImwA5vDVV2f1j/IlampqUkREhFv20XFcSi79iUIGhbu8nfYzrfrgzqd7VWtGRoYmTpyo5557TpJks9mUkJCgBx54QIsWLXJqG9/85jc1bdo0PfHEE12+/89//lPp6ek6dOiQRo4c6dQ2mTkBPIjQAcAs2traVF1drfz8fHtbcHCwsrKyVFlZedH1DcPQu+++q7q6Oq1YsaLbfk1NTQoKCtLQoUOdro1wAvQRgQOAmTQ3Nzu8Dg8PV3h45xmZEydOqL29XdHR0Q7t0dHR2rdvX7fbb2pqUnx8vFpbWxUSEqJf/epX+va3v91l37Nnz2rhwoW66667ejXzRDgBukDgAOBpLYciFGyxuLy+7exZSVJCQoJDe2FhoZYsWdKX0hxccsklqq2t1enTp1VeXi6r1arRo0drypQpDv3OnTunO+64Q4ZhaO3atb3aB+EEAYPAASAQHD582GGWoqtZE0mKiopSSEiIGhsbHdobGxsVExPT7faDg4M1ZswYSVJKSoo+/PBDFRUVOYSTjmBy6NAhvfvuu72+XodwAp9G4AAARxEREU6FgbCwMKWmpqq8vFwzZsyQdP6C2PLycs2bN8/p/dlsNrW2ttpfdwST/fv3a+vWrRo+fHivx0A4gSkROgDA/axWq2bPnq20tDSlp6eruLhYLS0tysvLkyTl5uYqPj5eRUVFkqSioiKlpaXp8ssvV2trq9544w1t2rTJftrm3Llz+v73v6+amhq99tpram9vV0NDg6TzP0MOC3PuV7KEE3gUoQMAzCMnJ0fHjx9XQUGBGhoalJKSorKyMvtFsvX19QoO/s/9WltaWjR37lx9+umnGjhwoMaOHauXXnpJOTk5kqQjR45oy5Ytks6f8rnQ1q1bO12X0h3uc4J+Q/AA4G88eZ+TxBVP9vmC2EMLH3drrZ7CzAmcQvAAAHgK4QSSCB8AAPMgnAQAggcAwJcQTnwcwQMA4G8IJyZH+AAABBrCiZcRPgAAcEQ4cTPCBwAAvUM46QcEEAAA+g/hxEkEEAAAPINwcgECCAAA3hew4eTL5FiFhrp+m2AAAOAewRfvAgAA4DmEEwAAYCqEEwAAYCqEEwAAYCqEEwAAYCqmCCdr1qxRUlKSLBaLMjIyVFVV5dR6paWlCgoK0owZM9xbIAAA8Bivh5PNmzfLarWqsLBQNTU1mjBhgrKzs3Xs2LEe1zt48KAeffRR3XTTTR6qFAAAeILL4eTvf/+7fvjDHyozM1NHjhyRJG3atEnbtm3r1XZWrVqlOXPmKC8vT8nJySopKdGgQYO0YcOGbtdpb2/XzJkztXTpUo0ePdrVIQAAABNyKZz86U9/UnZ2tgYOHKhdu3aptbVVktTU1KSnnnrK6e20tbWpurpaWVlZ/ykoOFhZWVmqrKzsdr1ly5ZpxIgRuvvuuy+6j9bWVjU3NzssAADAvFwKJ08++aRKSkq0bt06DRgwwN5+ww03qKamxuntnDhxQu3t7YqOjnZoj46OVkNDQ5frbNu2TevXr9e6deuc2kdRUZEiIyPtS0JCgtP1AQAAz3MpnNTV1enmm2/u1B4ZGamTJ0/2taZunTp1SrNmzdK6desUFRXl1Dr5+flqamqyL4cPH3ZbfQAAoO9cerZOTEyMDhw4oKSkJIf2bdu29eoakKioKIWEhKixsdGhvbGxUTExMZ36f/zxxzp48KBuvfVWe5vNZpMkhYaGqq6uTpdffrnDOuHh4QoPD3e6JgAA4F0uzZzMmTNH8+fP186dOxUUFKTPPvtMv/vd7/Too4/qvvvuc3o7YWFhSk1NVXl5ub3NZrOpvLxcmZmZnfqPHTtWu3fvVm1trX2ZPn26brnlFtXW1nLKBgAAP+DSzMmiRYtks9n0rW99S2fOnNHNN9+s8PBwPfroo3rggQd6tS2r1arZs2crLS1N6enpKi4uVktLi/Ly8iRJubm5io+PV1FRkSwWi8aNG+ew/tChQyWpUzsAAPBNLoWToKAg/fSnP9WCBQt04MABnT59WsnJyRoyZEivt5WTk6Pjx4+roKBADQ0NSklJUVlZmf0i2fr6egUHe/12LAAAwENcCicdwsLClJyc3Oci5s2bp3nz5nX5XkVFRY/rbty4sc/7B+CfvrgqzNslwMe1t9qk8ov3Q//qUzgBEBg4yAPwJMIJYGKEAgCBiHACOImgAACeQTiB3yA8AIB/IJzA4wgRAICeEE7gNEIFAMATCCcBhoABADA7womPImQAAPwV4cQECBoAAPwH4cQNCBsAALiOcNILhA4AANwv4MMJgQMAAHMJ2HByckyYQsIJJgAAmE2wtwsAAAC4EOEEAACYCuEEAIAAtmbNGiUlJclisSgjI0NVVVXd9l23bp1uuukmDRs2TMOGDVNWVlaP/e+9914FBQWpuLi4VzURTgAACFCbN2+W1WpVYWGhampqNGHCBGVnZ+vYsWNd9q+oqNBdd92lrVu3qrKyUgkJCZo6daqOHDnSqe+rr76qHTt2KC4urtd1EU4AAAhQq1at0pw5c5SXl6fk5GSVlJRo0KBB2rBhQ5f9f/e732nu3LlKSUnR2LFj9etf/1o2m03l5eUO/Y4cOaIHHnhAv/vd7zRgwIBe10U4AQDAjzQ3Nzssra2tXfZra2tTdXW1srKy7G3BwcHKyspSZWWlU/s6c+aMzp07p0svvdTeZrPZNGvWLC1YsEDXXHONS2MI2J8SAwBgJkP+FayQcNfnDNpbz6+bkJDg0F5YWKglS5Z06n/ixAm1t7crOjraoT06Olr79u1zap8LFy5UXFycQ8BZsWKFQkND9eCDD/ZyBP9BOAEAwI8cPnxYERER9tfh4eFu2c/y5ctVWlqqiooKWSwWSVJ1dbWeffZZ1dTUKCgoyOVtc1oHAAA/EhER4bB0F06ioqIUEhKixsZGh/bGxkbFxMT0uI+VK1dq+fLleuuttzR+/Hh7+9///ncdO3ZMI0eOVGhoqEJDQ3Xo0CE98sgjSkpKcnoMhBMAAAJQWFiYUlNTHS5m7bi4NTMzs9v1nn76aT3xxBMqKytTWlqaw3uzZs3S+++/r9raWvsSFxenBQsW6M0333S6Nk7rAAAQoKxWq2bPnq20tDSlp6eruLhYLS0tysvLkyTl5uYqPj5eRUVFks5fT1JQUKCXX35ZSUlJamhokCQNGTJEQ4YM0fDhwzV8+HCHfQwYMEAxMTG66qqrnK6LcAIAQIDKycnR8ePHVVBQoIaGBqWkpKisrMx+kWx9fb2Cg/9zkmXt2rVqa2vT97//fYftdHfRrasIJwAABLB58+Zp3rx5Xb5XUVHh8PrgwYO93r4r63DNCQAAMBXCCQAAMBXCCQAAMBXCCQAAMBXCCQAAMBXCCQAAMBXCCQAAMBXCCQAAMBXCCQAAMBXCCQAAMBXCCQAAMBWerYMenbrc5u0SAsYlH/P/CgAgEU7cjoM7nOWrf1YIVQD6W8CGk9OjbQq2+ObBADATM4QqAhLgXwI2nADwH54KSIQgwDMIJwDgpP4MQQQdoHuEEwDwgv4KOoQc+CPCCQD4sP4+pUXYgRkQTgAAdpy6ghkQTgAAbkHQgasIJwAA0+uPoEPA8R2EEwBAQHAl4NjOev8+PoGIGAkAAEyFcAIAAEzFFOFkzZo1SkpKksViUUZGhqqqqrrt+8orrygtLU1Dhw7V4MGDlZKSok2bNnmwWgAA4E5eDyebN2+W1WpVYWGhampqNGHCBGVnZ+vYsWNd9r/00kv105/+VJWVlXr//feVl5envLw8vfnmmx6uHAAAuIPXw8mqVas0Z84c5eXlKTk5WSUlJRo0aJA2bNjQZf8pU6bo9ttv19VXX63LL79c8+fP1/jx47Vt2zYPVw4AANzBq+Gkra1N1dXVysrKsrcFBwcrKytLlZWVF13fMAyVl5errq5ON998sztLBQAAHuLVnxKfOHFC7e3tio6OdmiPjo7Wvn37ul2vqalJ8fHxam1tVUhIiH71q1/p29/+dpd9W1tb1draan/d3NzcP8UDAAC38Mn7nFxyySWqra3V6dOnVV5eLqvVqtGjR2vKlCmd+hYVFWnp0qWeLxIAALjEq+EkKipKISEhamxsdGhvbGxUTExMt+sFBwdrzJgxkqSUlBR9+OGHKioq6jKc5Ofny2q12l83NzcrISGhfwYAAAD6nVevOQkLC1NqaqrKy8vtbTabTeXl5crMzHR6OzabzeHUzYXCw8MVERHhsAAAAPPy+mkdq9Wq2bNnKy0tTenp6SouLlZLS4vy8vIkSbm5uYqPj1dRUZGk86dp0tLSdPnll6u1tVVvvPGGNm3apLVr13pzGAAAoJ94PZzk5OTo+PHjKigoUENDg1JSUlRWVma/SLa+vl7Bwf+Z4GlpadHcuXP16aefauDAgRo7dqxeeukl5eTkeGsIAACgHwUZhmF4uwhPam5uVmRkpBJXPKlgi8Xb5QAATMx29qwOLXxcTU1NbrssoOO4dM2Pn1JIuOvHpfbWs9r7/GNurdVTvH4TNgAAgAsRTgAAgKkQTgAAgKkQTgAAgKkQTgAAgKkQTgAAgKkQTgAAgKkQTgAACGBr1qxRUlKSLBaLMjIyVFVV1W3fvXv36nvf+56SkpIUFBSk4uLiLvsdOXJEP/zhDzV8+HANHDhQ1157rd577z2na/L6HWK9ZXBis0IGdf08Hl91+pNIb5cAAPAhmzdvltVqVUlJiTIyMlRcXKzs7GzV1dVpxIgRnfqfOXNGo0eP1g9+8AM9/PDDXW7ziy++0A033KBbbrlF//d//6dvfOMb2r9/v4YNG+Z0XQEbTvzRkFFNHt0fYQgAfNuqVas0Z84c+/PsSkpK9Prrr2vDhg1atGhRp/4TJ07UxIkTJanL9yVpxYoVSkhI0IsvvmhvGzVqVK/qIpzAZf0Zhgg6ANA/mpubHV6Hh4crPDy8U7+2tjZVV1crPz/f3hYcHKysrCxVVla6vP8tW7YoOztbP/jBD/TXv/5V8fHxmjt3rubMmeP0NggnMIX+CDoEHAC+bOiBNoWGun4p6FdftUmSEhISHNoLCwu1ZMmSTv1PnDih9vZ2+4N2O0RHR2vfvn0u1/Gvf/1La9euldVq1WOPPaZ//vOfevDBBxUWFqbZs2c7tQ3CCfxGXwMO4QaAPzh8+LDDg/+6mjVxJ5vNprS0ND311FOSpOuuu0579uxRSUkJ4QToLWZvAPiDiIgIp55KHBUVpZCQEDU2Njq0NzY2KiYmxuX9x8bGKjk52aHt6quv1p/+9Cent0E4AfoRszcAfEVYWJhSU1NVXl6uGTNmSDo/61FeXq558+a5vN0bbrhBdXV1Dm0fffSREhMTnd4G4QQwkb6EG4INgN6yWq2aPXu20tLSlJ6eruLiYrW0tNh/vZObm6v4+HgVFRVJOn8R7QcffGD/9yNHjqi2tlZDhgzRmDFjJEkPP/ywJk2apKeeekp33HGHqqqq9MILL+iFF15wui7CCeAnCDYAeisnJ0fHjx9XQUGBGhoalJKSorKyMvtFsvX19QoO/s9Fup999pmuu+46++uVK1dq5cqVmjx5sioqKiSd/7nxq6++qvz8fC1btkyjRo1ScXGxZs6c6XRdQYZhGP0zRN/Q3NysyMhIJZf+RCGDPHuREOBPCDQIBLazZ3Vo4eNqampy6joOV3Qcl2741hKFhlpc3s5XX53VP8qXuLVWT2HmBIBLejNTQ5AB0BuEEwBu50yQIcAA6EA4AWAKBBgAHQgnAHwGAQYIDIQTAH6FAAP4PsIJgIBzsQBDeAG8i3ACAF/D7AvgXYQTAHABsy+A+xBOAMANCC+A6wgnAOAFhBege4QTADAhrntBICOcAICPYvYF/opwAgB+itkX+CrCCQAEMGZfYEaEEwBAt5h9gTcQTgAAfcLsC/ob4QQA4FbMvqC3CCcAAK9j9gUXIpwAAEyP8BJYCCcAAJ/XU3ghuPiegA0n/3XZfoUPGaC36q/ydikAADdi1sX3BGw46TB1ZJ1T/QgxAOCfegovzR+Ge7ASdAj4cOIsZ0OMRJABAH8xOLHZ2yUEJMKJGzAbAwCA6wgnXuRMiCHAAAACDeHE5JiFAQAEGsKJn2AWBgDgLwgnAYQAAwDwBYQTOLhYgCG8AADcjXCCXiG8AADcjXCCftVTeCG4AACcQTiBxzDrAgBwBuEEpsGsCwBAkoK9XYAkrVmzRklJSbJYLMrIyFBVVVW3fdetW6ebbrpJw4YN07Bhw5SVldVjf/iHqSPrul0AAP7F6zMnmzdvltVqVUlJiTIyMlRcXKzs7GzV1dVpxIgRnfpXVFTorrvu0qRJk2SxWLRixQpNnTpVe/fuVXx8vBdGAG9jxgUA/EuQYRiGNwvIyMjQxIkT9dxzz0mSbDabEhIS9MADD2jRokUXXb+9vV3Dhg3Tc889p9zc3Iv2b25uVmRkpOZtu13hQwb0uX74LoILgItpP9OqD+58Wk1NTYqIiHDLPjqOSzd8a4lCQy0ub+err87qH+VL3Fqrp3h15qStrU3V1dXKz8+3twUHBysrK0uVlZVObePMmTM6d+6cLr300i7fb21tVWtrq/11czNPmMR53c24EFoAwLu8Gk5OnDih9vZ2RUdHO7RHR0dr3759Tm1j4cKFiouLU1ZWVpfvFxUVaenSpX2uFYGD00QA4F1ev+akL5YvX67S0lJVVFTIYul6Kiw/P19Wq9X+urm5WQkJCZ4qEX6G2RYAcD+vhpOoqCiFhISosbHRob2xsVExMTE9rrty5UotX75c77zzjsaPH99tv/DwcIWHh/dLvUB3CC0A0H+8Gk7CwsKUmpqq8vJyzZgxQ9L5C2LLy8s1b968btd7+umn9bOf/Uxvvvmm0tLSPFQt0HuEFgDoPa+f1rFarZo9e7bS0tKUnp6u4uJitbS0KC8vT5KUm5ur+Ph4FRUVSZJWrFihgoICvfzyy0pKSlJDQ4MkaciQIRoyZIjXxgH0BqEFALrn9Zuw5eTkaOXKlSooKFBKSopqa2tVVlZmv0i2vr5eR48etfdfu3at2tra9P3vf1+xsbH2ZeXKld4aAtBvuNEcAE/rzY1QJemPf/yjxo4dK4vFomuvvVZvvPGGw/unT5/WvHnzdNlll2ngwIFKTk5WSUlJr2ry+syJJM2bN6/b0zgVFRUOrw8ePOj+ggCTYaYFgDv09kao27dv11133aWioiJ997vf1csvv6wZM2aopqZG48aNk3T+jMi7776rl156SUlJSXrrrbc0d+5cxcXFafr06U7V5fWbsHkaN2FDICC0AP3D32/C1tsboebk5KilpUWvvfaave36669XSkqKfXZk3LhxysnJ0eLFi+19UlNT9Z3vfEdPPvmkU3V5/bQOgP7H6SEgcDU3NzssF96I9EIdN0K98D5hF7sRamVlZaf7imVnZzv0nzRpkrZs2aIjR47IMAxt3bpVH330kaZOner0GExxWgeAZ3B6CDCvgR8cVWhwmMvrf2Vrk6RO9/IqLCzUkiVLOvV35UaoDQ0NXfbv+HGKJK1evVr33HOPLrvsMoWGhio4OFjr1q3TzTff7PRYCCcACC2AHzl8+LDDaR1P3+tr9erV2rFjh7Zs2aLExET97W9/0/3339/j3dy/jnACoFtdhRYCC2BuERERTl1z4sqNUGNiYnrs/+WXX+qxxx7Tq6++qmnTpkmSxo8fr9raWq1cudLpcMI1JwB6hetZAP9w4Y1QO3TcCDUzM7PLdTIzMx36S9Lbb79t73/u3DmdO3dOwcGO8SIkJEQ2m83p2pg5AdAvODUE+J7e3gh1/vz5mjx5sn7+859r2rRpKi0t1XvvvacXXnhB0vlZm8mTJ2vBggUaOHCgEhMT9de//lW//e1vtWrVKqfrIpwAcCtODQHmlZOTo+PHj6ugoEANDQ1KSUnpdCPUC2dBJk2apJdfflmPP/64HnvsMV1xxRX685//bL/HiSSVlpYqPz9fM2fO1Oeff67ExET97Gc/07333ut0XdznBIApEFhgRp68z0lW7I/7/Gudd44+79ZaPYWZEwCmwAwLgA6EEwCmRWABAhPhBIBP+XpgIawA/odwAsCnMbsC+B/CCQC/w+wK4NsCNpzcOXSnhlzyn59HbfpikherAeBOhBXAtwRsOPm6WcO2d9lOaAH8D2EFMDfCyUV0FVoILIB/IawA5kI4ccHXAwthBfAvhBXAuwgn/YDZFcC/EVYAzyKcuAmBBfBfhBXAvQgnHsTpIMA/XRhWCCpA3xFOvIiwAvgfZlWAviOcmAinggD/Q1gBeo9wYnLMrgD+hVNAwMURTnwMYQXwH8yqAF0jnPg4wgrgPwgrwHmEEz9DWAH8B6eAEKgIJ36OsAL4B2ZVEEgIJwHmwrBCUAF8F7Mq8GeEkwDGrArgH5hVgb8hnMCOsAL4B2ZV4OsIJ+gWp4AA38esCnwR4QROYVYF8A+EFfgCwglcwqwK4B84BQQzIpygz5hVAfwDsyowC8IJ+h2zKoB/YFYF3kI4gVsRVAD/wKwKPIlwAo/h9A/gP5hVgTsRTuA1zKoA/oFZFfQ3wglMgaAC+A9mVdBXhBOYDkEF8B8EFbiCcAJTI6gA/oPTP3AW4QQ+g6AC+BdmVdAdwgl8Er/8AfwLQQUXIpzALzCrAvgPggoIJ/A7BBXAf3CdSmAinMCvEVQA/8KsSmAgnCBgEFQA/0JQ8V/B3i4A8IZZw7Z3uqgWgO+aOrKu0ykg+C5mThDQmE0B/AuzKf6BmRPg/2M2BfAvzKb4Lq+HkzVr1igpKUkWi0UZGRmqqqrqtu/evXv1ve99T0lJSQoKClJxcbHnCkXA6AgpBBXAP3SEFIJK13pzHJakP/7xjxo7dqwsFouuvfZavfHGGw7vG4ahgoICxcbGauDAgcrKytL+/ft7VZNXw8nmzZtltVpVWFiompoaTZgwQdnZ2Tp27FiX/c+cOaPRo0dr+fLliomJ8XC1CESEFMC/EFQc9fY4vH37dt111126++67tWvXLs2YMUMzZszQnj177H2efvpp/fKXv1RJSYl27typwYMHKzs7W2fPnnW6riDDMIw+j85FGRkZmjhxop577jlJks1mU0JCgh544AEtWrSox3WTkpL00EMP6aGHHurVPpubmxUZGalte+I05BKvTxzBB3FtCuBfero2pf1Mqz6482k1NTUpIiLCLfvvOC5lxf5YocFhLm/nK1ub3jn6fK9q7e1xOCcnRy0tLXrttdfsbddff71SUlJUUlIiwzAUFxenRx55RI8++qgkqampSdHR0dq4caPuvPNOp+ry2gWxbW1tqq6uVn5+vr0tODhYWVlZqqys7Lf9tLa2qrW11f66qalJktRy2tZv+0BguT10mySp9GSGlysB0B8mX3r+//rf/fSKTu+1nzl//PDE/8d/ZbRJfTg0fWW0STofdi4UHh6u8PDwTv1dOQ5XVlbKarU6tGVnZ+vPf/6zJOmTTz5RQ0ODsrKy7O9HRkYqIyNDlZWV5g8nJ06cUHt7u6Kjox3ao6OjtW/fvn7bT1FRkZYuXdqpPfv6hn7bBwLVq94uAICH/Pvf/1ZkZKRbth0WFqaYmBhVNLzY520NGTJECQkJDm2FhYVasmRJp76uHIcbGhq67N/Q0GB/v6Otuz7O8PufEufn5zukvJMnTyoxMVH19fVu+4PmCc3NzUpISNDhw4fdNtXobv4wBolxmIk/jEHyj3H4wxik87PtI0eO1KWXXuq2fVgsFn3yySdqa2vr87YMw1BQUJBDW1ezJmbntXASFRWlkJAQNTY2OrQ3Njb268Wu3U1nRUZG+vQXpkNERITPj8MfxiAxDjPxhzFI/jEOfxiDdP50hztZLBZZLBa37uPrXDkOx8TE9Ni/45+NjY2KjY116JOSkuJ0bV67IjQsLEypqakqLy+3t9lsNpWXlyszM9NbZQEAEBBcOQ5nZmY69Jekt99+295/1KhRiomJcejT3NysnTt39urY7tXTOlarVbNnz1ZaWprS09NVXFyslpYW5eXlSZJyc3MVHx+voqIiSecv3vnggw/s/37kyBHV1tZqyJAhGjNmjNfGAQCAL+rtcXj+/PmaPHmyfv7zn2vatGkqLS3Ve++9pxdeeEGSFBQUpIceekhPPvmkrrjiCo0aNUqLFy9WXFycZsyY4XxhhpetXr3aGDlypBEWFmakp6cbO3bssL83efJkY/bs2fbXn3zyiSGp0zJ58mSn93f27FmjsLDQOHv2bD+OwvP8YRz+MAbDYBxm4g9jMAz/GIc/jMEw/GccPenNcdgwDOMPf/iDceWVVxphYWHGNddcY7z++usO79tsNmPx4sVGdHS0ER4ebnzrW98y6urqelWTV+9zAgAA8HXchQwAAJgK4QQAAJgK4QQAAJgK4QQAAJiKT4eTv/3tb7r11lsVFxenoKAg+739OxguPra5t4+P7it3jGPJkiUKCgpyWMaOHevGUVx8HK+88oqmTp2q4cOHKygoSLW1tU5t92KP5+5P7hjDxo0bO30W7r7ZUk/jOHfunBYuXKhrr71WgwcPVlxcnHJzc/XZZ59ddLue/G64Ywxm/F4sWbJEY8eO1eDBgzVs2DBlZWVp586dF92u2f6ecmUcnv48LjaGC917770KCgpScXHxRbfr6c8iEPh0OGlpadGECRO0Zs2aLt935bHNvX18dH9wxzgk6ZprrtHRo0fty7Zt29xRvt3FxtHS0qIbb7xRK1ascHqbzjyeuz+5YwzS+btkXvhZHDp0qD/K7VZP4zhz5oxqamq0ePFi1dTU6JVXXlFdXZ2mT5/e4zY9/d1wxxgk830vrrzySj333HPavXu3tm3bpqSkJE2dOlXHjx/vdptm/HvKlXFInv08LjaGDq+++qp27NihuLi4i27TG59FQHDpR9EmJMl49dVX7a9tNpsRExNjPPPMM/a2kydPGuHh4cbvf//7breTnp5u3H///fbX7e3tRlxcnFFUVOSWur+uv8ZRWFhoTJgwwY2V9uzr47hQx/1qdu3addHt3HHHHca0adMc2jIyMowf//jH/VBlz/prDC+++KIRGRnZr7X1Rk/j6FBVVWVIMg4dOtRtH29+N/prDGb+XnRoamoyJBnvvPNOt33M9vdUV5wZhzc/j+7G8Omnnxrx8fHGnj17jMTEROMXv/hFj9vx9mfhr3x65qQnF3tsc1c6Hh994ToXe3y0u7kyjg779+9XXFycRo8erZkzZ6q+vt7d5fa7yspKh7FL5x/P7a3Pw1WnT59WYmKiEhISdNttt2nv3r3eLslBU1OTgoKCNHTo0C7fN+N34+suNoYOZv5etLW16YUXXlBkZKQmTJjQbR+zfxbOjKODmT4Pm82mWbNmacGCBbrmmmsu2t8XPgtf5bfhxJXHNvf0+OjePOq5P7n6+OmMjAxt3LhRZWVlWrt2rT755BPddNNNOnXqlFvr7W8Xezy3L7jqqqu0YcMG/eUvf9FLL70km82mSZMm6dNPP/V2aZKks2fPauHChbrrrru6fUCbGb8bF3JmDJJ5vxevvfaahgwZIovFol/84hd6++23FRUV1WVfM38WvRmHZL7PY8WKFQoNDdWDDz7oVH8zfxa+zqvP1oH7fOc737H/+/jx45WRkaHExET94Q9/0N133+3FygJPZmamwwOvJk2apKuvvlrPP/+8nnjiCS9Wdv7C0jvuuEOGYWjt2rVercVVvRmDWb8Xt9xyi2pra3XixAmtW7dOd9xxh3bu3KkRI0Z4rSZX9HYcZvo8qqur9eyzz6qmpkZBQUEe3Tc689uZkwsf23yhnh4F7crjo93NlXF0ZejQobryyit14MCBfq3P3S72eG5fNGDAAF133XVe/yw6DuqHDh3S22+/3eOMgxm/G1LvxtAVs3wvBg8erDFjxuj666/X+vXrFRoaqvXr13fZ16yfhdS7cXTFm5/H3//+dx07dkwjR45UaGioQkNDdejQIT3yyCNKSkrqch0zfxa+zm/DiSuPbXbl8dHu1l+Pnz59+rQ+/vhjxcbGuqNMt7nY47l9UXt7u3bv3u3Vz6LjoL5//3698847Gj58eI/9zfjd6O0YumLW74XNZlNra2uX75nxs+hOT+Poijc/j1mzZun9999XbW2tfYmLi9OCBQv05ptvdrmOL30WPsfbV+T2xalTp4xdu3YZu3btMiQZq1atMnbt2mW/Wn/58uXG0KFDjb/85S/G+++/b9x2223GqFGjjC+//NK+jf/6r/8yVq9ebX9dWlpqhIeHGxs3bjQ++OAD45577jGGDh1qNDQ0+NQ4HnnkEaOiosL45JNPjH/84x9GVlaWERUVZRw7dsxr4/j3v/9t7Nq1y3j99dcNSUZpaamxa9cu4+jRo/ZtzJo1y1i0aJH99T/+8Q8jNDTUWLlypfHhhx8ahYWFxoABA4zdu3f7zBiWLl1qvPnmm8bHH39sVFdXG3feeadhsViMvXv3umUMFxtHW1ubMX36dOOyyy4zamtrjaNHj9qX1tZW+za8/d1wxxjM9r04ffq0kZ+fb1RWVhoHDx403nvvPSMvL88IDw839uzZ0+04zPb3lKvj8PTncbHv99d19WsdM3wWgcCnw8nWrVsNSZ2Wjsc7O/PY5sTERKOwsNChrafHR/vKOHJycozY2FgjLCzMiI+PN3JycowDBw54dRwvvvhil+9fWLcrj+c2+xgeeugh+5+n6Oho47//+7+Nmpoat43hYuPo+Bl0V8vWrVvt2/D2d8MdYzDb9+LLL780br/9diMuLs4ICwszYmNjjenTpxtVVVUO2/D2Z+GucXj687jY9/vrugonZvgsAkGQYRiGKzMuAAAA7uC315wAAADfRDgBAACmQjgBAACmQjgBAACmQjgBAACmQjgBAACmQjgBAACmQjgBAACmQjgBAACmQjgBAACmEurtAgB4z5QpUzRu3DhJ0qZNmzRgwADdd999WrZsmYKCgrxcHYBAxcwJEOB+85vfKDQ0VFVVVXr22We1atUq/frXv/Z2WQACGA/+AwLYlClTdOzYMe3du9c+U7Jo0SJt2bJFH3zwgZerAxComDkBAtz111/vcAonMzNT+/fvV3t7uxerAhDICCcAAMBUCCdAgNu5c6fD6x07duiKK65QSEiIlyoCEOgIJ0CAq6+vl9VqVV1dnX7/+99r9erVmj9/vrfLAhDA+CkxEOByc3P15ZdfKj09XSEhIZo/f77uueceb5cFIIARToAAN2DAABUXF2vt2rXeLgUAJHFaBwAAmAzhBAAAmAo3YQMAAKbCzAkAADAVwgkAADAVwgkAADAVwgkAADAVwgkAADAVwgkAADAVwgkAADAVwgkAADAVwgkAADCV/wcppKa6ZP69MAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "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": "few-venv", "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.3" } }, "nbformat": 4, "nbformat_minor": 2 }