{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fourier Formalism\n",
    "\n",
    "The Fourier formalism is described by [Marsh et al 2022](https://ui.adsabs.harvard.edu/abs/2021arXiv210708040M/abstract) and in this documentation.\n",
    "\n",
    "ALPro includes a _numerical_ implementation, as described in section 6 of [Marsh et al 2022](https://ui.adsabs.harvard.edu/abs/2021arXiv210708040M/abstract). Standard discrete cosine/sine transforms, as implemented in [scipy.fftpack](https://docs.scipy.org/doc/scipy/reference/fftpack.html), are used to obtain first-order solutions to the Schrodinger-like equation of motion for a photon-ALP beam. This notebook gives a couple of basic examples of how to use the scheme. \n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "Warning: the scheme is _perturbative_ and so only accurate to first-order! For $P_{\\gamma a} \\gtrsim 0.1$ the calculation will give inaccurate results. \n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "Note: The Fourier formalism is well tested in general, but the specific numerical implementation in ALPro should at this stage be considered a beta feature. If you discover any problems, please open a Github issue. \n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "Warning: one of the idiosyncrasies of the scheme is that the use of a discrete cosine transform requires regular samples in either $z$, $l$ or $\\varphi$ space. As a result, it can sometimes be impossible to repdroduce exactly results from models with variable cell sizes. It is therefore recommended that you work with models that can be broken up into regular spacing grids and avoid step changes at irregular intervals. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Massive ALPs\n",
    "\n",
    "This example sets up a cell model with uniform cell size and finds the conversion probability. The massless formalism is specified automatically by setting the mass to a value below the minimum $\\omega_{\\rm pl}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdcAAAFECAYAAAB1dnQsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXxcdb3/8ddntux70nTfC7RsLeUKZWsF6wIKKiquiOtFrrijonj1uuAVVBAU4SL3pyAqV3EBFdkLCFVsWUoXoLSla5JmafZkkpn5/v6YycxkT5qZJNO8n49HHnznnM/3O5851nxyzpzz/ZpzDhEREUkdz0QnICIicqRRcRUREUkxFVcREZEUU3EVERFJMd9EJyAjY2bPAhXAKxOdi4iIsBiodc6tGGinimvmqCgqKpq1fPnyWROdiIjIVPfcc8/R1NQ06H4V18zxyvLly2etW7duovMQEZny1qxZw2OPPTbolUR95yoiIpJiKq4iIiIppuIqIiKSYiquIiIiKabiKiIikmIqriIiIimm4ioiIpJiKq4iIiIppuIqIkes5s5u1u+oZ3d920SnIlOMZmgSkSPSrro2LrplPQdbgng9xn+//XjeefKciU5LpgiduYrIEenK32/iYEsQgHDE8dU/bmZ/Q8sEZyVThYqriBxxtuyr5x87G+Kvj7ed3GLfZc9dV0xgVjKV6LKwiBxxmv92NQ8G7uPByEr2+eZwNT8B4GDNXsKhEF6ffvVJeunMVUSOOBVVj7LEs5/LfPdw4ckLqKcIgGk0sO0ff5vg7GQqUHEVkSNKU2M9C0M7AQg5D0efdSGvlJ0d39+y9cGJSk2mEBVXETmi7N70BB5zALzqW0B+YSmBoxLFtaT2nxOVmkwhKq4ickRpfWV9vF1fciIAC05aG9+2qOtl2loaxz0vmVpUXEXkiJJVvyXe9sw+CYDiihns8swHwG9hdmx8eCJSkylExVVEjijlHbvi7ZJ5J8bbB8tOjrfbdv5jXHOSqUfFVUSOGMFgB7PCB+KvZy1ZHm97Y2exANl1L4xrXjL1qLiKyBFj/47N+CwCwAGbRk5+YXxfxVGnxNuz2l8a99xkalFxFZEjRsOuTfF2Xfb8XvtmLzmRDhcAos+71lXvHc/UZIrRNCUicsToqnk53m4vWtxrn9fn587CD/N8vY8X3AK+fsjPa6ePd4YyVai4isgR4968t/P14ELm2EHevujUfvv3Lv4A99buBmDzgRZeu1TVVdJDl4VF5IixqynCdjebRyInkT/3xH77j5mR+A52R23reKYmU8ykLK5mdpmZ7TKzTjPbaGZnDhO/OhbXaWY7zezS0Y5pZllmdqOZ1ZlZm5ndY2az+8TMNbN7Y/vrzOwGMwuMJhcz85rZt5Jy2WVm3zYzXUUQGaO9DR3x9pzS3H77F1Xkx9s7arWAuqTPpCuuZnYR8CPgamAF8BRwn5nNHSR+AfDXWNwK4LvAjWZ24SjHvB64EHgPcCZQCPzZzLyxMbzAX4CC2P73AO8AfjCaXIAvAf8BfAo4Bvh07PWVozhMItJHdzhCVVOiuM4qzukXs7AiDwAjQkftLlwkMm75ydQyGc+WPgf83Dl3a+z15Wb2RuATDFyALgUOOOcuj73eZmanAF8A7h7JmGZWBHwE+JBz7kEAM/sAsBt4HXA/8HrgWGCec25vLOaLwM/M7KvOueYR5nIacK9z7t7Y61fN7B4g8ZxAEjNbF2suH2i/iETV1B9ihqulmlLKCnLJ9nv7xZTlBfht9rc5zr1CjnVxsOpMps2aP/7JyhFvUp25xi6xrgQe6LPrAaJFaSCrBoi/HzjZzPwjHHMl4E+OiRXQbUkxq4BtPYU16X2yYv2HzSX2+u/Aa83smNhnXgacTfSMV0QOU+PLT/Fk9qd5MesSfuT54YAxZkahL0yOdQFwMOnRHZFUmlTFFSgHvEBNn+01wGC39U0fJN4XG28kY04HwkDdMDF9x6iL9RsqJjkXgO8BdwBbzawb2AL8wjl300Afzjm3xjm3BnhuoP0iEtV+MLrMnN/C+LP6XxLu0Zy/IN5u2/9i2vOSqWmyFdep4CLgYuC9wEmx9mVm9pEJzUokw4Ub9sTb3QWzB48rWRRvu7qXB40TGYvJVlx7zgQr+2yvBKoH6VM9SHwoNt5IxqwmenZbPkxM3zF6zoqHiknOBeBa4PvOud84515wzt0B/BDd0CQyJt62xK8IT9GsQeOyph8Tb+e27Bo0TmQsJlVxdc51ARuBtX12rSV6B+5A1g8Sv8E51z3CMTcC3ckxscdwlibFrAeW9nk8Zy0QjPUfNpfY61yixT5ZmEn2v4VIpsnqPJholwxeXEvnHRtvV3TuTmtOMnVNxruFfwjcYWZPA08SvQN3JnAzgJndDuCcuzgWfzPwSTO7HrgFOB24hOijMiMa0znXZGa3AdeY2UGgPtZnE/BQbIwHiH4/eruZfR4oI3oWemvsTuGR5nIv8GUz2xUbbwXRu5lvP7zDJSIA+V218XZe+eCXhWfMX0bYGV5zVLo6OtpayMkrGI8UZQqZdMXVOXeXmZUBVwEzgM3Auc65nj8x5/aJ32Vm5wLXEX205gDwKefc3aMYE+AzRC/f3gXkAA8DFzvnwrExwmZ2HnAT0QLdAdwJXDGaXIDLgW/FxpkGVAG3At88jMMlIjEl4fp4u7hy3qBxgewc9nmmM9tV4THHgZ1bWHR8/6kSRcZi0hVXgNids4PePTvAtseI3hx0WGPG9geJFr7Lh4jZA7x5mPcZMhfnXAvRQv6ZocYRkZELBjsoJXoBKeyM0mmDXxYGqM+ey+yOKgAO7dkCKq6SYvqeT0QyXkPS8nENVozX5x8iGjqKEncMd9fojmFJPRVXEcl4TQcTj+E0+sqGjfdWLAGg3WXR0tqStrxk6pqUl4VFREajubGBdpdFrgVpDUwbNt573NtY9a9yqinhWCvm9eOQo0wtKq4ikvG25P4b7wr+L/l08L5jp7FimPj5s2dRRfQMd2dtG845zCz9icqUocvCIpLxqps7AaOVXPJLZw4bX5oXoDg3+r1se1c41l8kdVRcRSTjHWwOxtuVhdkj6pO8tutOre0qKabLwiKS8epaE8W1oiBrRH2OLvXSsedVFloVja84WDzkU3Yio6LiKiIZb0bDBlZYN3UUUp7bfx3Xgbwl+GeuzroBgH/uqGKYR9hFRkXFVUQy3hVt11KRdQiAajsLGP5xnOzpR0N0lTpym3emMTuZivSdq4hktEg4THF8em8oLp8xon6l846LtyuCe4eIFBk9FVcRyWjNh+rwW3ShqRZyyM7JG1G/mQuOIeSivwKnU0t7a1PacpSpR8VVRDJaU/2BeLvRSkbczx/IpsozPf66aueWlOYlU5uKq4hktNb6qkTbVzyqvvU5idVzDu1VcZXUUXEVkYzW2ViTaAdKR9e3cGG8rQn8JZVUXEUko3U3J4prV/bwdwkn81YcFW8HDr2SspxEVFxFJKO51tpEO6d8VH0LZi+Lt4vad6csJxEVVxHJaJ6Ounjb8itG1Xf6wuPj7ZmhfUTC4ZTlJVObJpEQkYzm76yPt32FlaPqW1wxg5eZy8FwATvdTNY2NDKjYnSXlkUGouIqIhmtOlLM9sgsyqyJ7JKRTSCR7CuVt7Bhd3R2pwWNEWaM7uRXZEAqriKS0a71fpRdXdFVbR5adMao+y+syIsX1521bZy5RNVVxk7fuYpIRkteEacsb2Qr4iRLXnpuR21rSnISUXEVkYwVDIVp6QwB4PUYRTn+UY+xUOu6ShrosrCIZKz61q54uzQvgMdjox5jcYmH93kfYpEdYMaBDuCeFGYoU5WKq4hkrOaDe/ig934OuQI8WfOG7zCA2WWFfMP3i+jk/xFoa2kkr2B00yiK9KXiKiIZq7tqC//l/wUAm7uXAx8f9Rj+QBa7vTOYF9kHwIEdL7Bk+ZmpTFOmIH3nKiIZK9icmJ2pMzDyFXH6asiZH2837t48lpREABVXEclgodbE7EyhrMMvrsHixYlxal4cU04ioOIqIpmsLTE7k8s5/JmVfJXHxNvZjdvHlJIIqLiKSCbrbIg3Lffwi2vxvMQcw2Wdr44lIxFAxVVEMpi/81C87SsY3Yo4yWYuSprAP1xFV7BzTHmJqLiKSMbK6m6MtwOFhz9tYW5+EVUW7e+zCAd2bhlzbjK1qbiKSMbKTSquuUVjmxO4Nmt+vN2w+4UxjSWi4ioiGSs/0pxol45uubm+2osWxdvBqm1jGktEk0iISEZykQjFrhliMx4WjbG4dsw7mx/v7+CVyCxKPKezKgU5ytSl4ioiGamts5P7I6dSQisFnk7+La9gTOMVLFvL9x+PTuK/tLEwFSnKFKbiKiIZ6VAnfL77MgBmFmXz1BjHWzwteXWcVsIRh/cwFgIQAX3nKiIZqqEtsSJOSV5gzOMV5wYoz4+uBxsMRdh3qH3MY8rUpeIqIhnpUHtScc0de3EFOKoyevZqRHh5f+0w0SKDU3EVkYzUq7im4MwV4O1Z/+LuwNd5Ieuj5G/4SUrGlKlJ37mKSEbK2b2Oq3x/5pAroCTyWmDFmMecV+BY6YnOLZxVr8dx5PCpuIpIRiqu/Rdv9N0HwPpgMfDeMY9ZsvAkeD7antb+ypjHk6lLl4VFJCN5OhKT9nvyDn9e4WSzj1pByEV/Lc6MVNPa0jhMD5GBqbiKSEbyBxOT9nvzU1Ncs3Py2OedBYDHHPte3JCScWXqmZTF1cwuM7NdZtZpZhvN7Mxh4lfH4jrNbKeZXTraMc0sy8xuNLM6M2szs3vMbHafmLlmdm9sf52Z3WBmgT4xI8llhpn9wsxqY3FbzWz16I6SyNSWPGl/9hjnFU5Wn39UvN2067mUjStTy6QrrmZ2EfAj4Gqidyg8BdxnZnMHiV8A/DUWtwL4LnCjmV04yjGvBy4E3gOcCRQCfzYzb2wML/AXoCC2/z3AO4AfjDKXYuBJopO2nQcsBS4HDo7uSIlMbbmhpng7u2haysbtLluWeFGjCfzl8EzGG5o+B/zcOXdr7PXlZvZG4BPAlQPEXwoccM5dHnu9zcxOAb4A3D2SMc2sCPgI8CHn3IMAZvYBYDfwOuB+4PXAscA859zeWMwXgZ+Z2Vedc80jzOWLQJVz7uKkz7BrlMdIZMorSJq0v2CM8wony523PP7/yJKmrSkbV6aWSXXmGrvEuhJ4oM+uB4DTBum2aoD4+4GTzcw/wjFXAv7kmFgB3ZYUswrY1lNYk94nK9Z/2Fxir98K/NPM7jKzg2b2nJl90swGnGfNzNaZ2Tpg+YCfXmQKcpEIRS5RXMc6aX+yucedEW/P795JMNiRsrFl6phUxRUoB7xATZ/tNcD0QfpMHyTeFxtvJGNOB8JA3TAxfceoi/UbKiY5F4CFwGXATuANRC9X/zfwH4N8PhHpo6W5AZ9FAGhz2WTn5KZs7OLy6ey3aLEOWIjdW/6VsrFl6phsxXUq8ADPOOeudM4965z7f8ANDFJcnXNrnHNrAN1ZIRLTUl8dbzd5Ur+CTXX+sdH3cTkc2KPnXWX0Jltx7TkT7HuNpxKo7h8Ose0DxYdi441kzGqiZ7d97+fvG9N3jJ6z4qFiknMBqAL6fpGzDRjwhi0R6a/1UOICUZunKOXj71p2GecEr+WE4K3c27Vy+A4ifUyq4uqc6wI2Amv77FoLg64otX6Q+A3Oue4RjrkR6E6OiT2GszQpZj2wtM/jOWuBYKz/sLnEXj8JHN0n5iiiN0+JyAg0UML1obfzi9Bans0f8km9wzJv6Up2uFk4PDy/TxNJyOhNxruFfwjcYWZPEy1ElwIzgZsBzOx2gKS7bW8GPmlm1wO3AKcDlxB9VGZEYzrnmszsNuAaMzsI1Mf6bAIeio3xALAFuN3MPg+UAdcCt8buFB5pLtcBT5nZV4G7iD6y8yngK4d3uESmngM2jetD7wDgghkzeVeKxz92ZiFejxGOOHbUttLS2U1Btn/4jgL1O6CjEcJBCAUJd3cS6uok3NVJqLuTULCTSKiTSHeQhtnn0Fa4iIiDSMThgOmbb8HfUQcOnAuDixB9EcFFXLyNi7Br8YdoKViIwxFx4JzjxGe+hjfcAa4nLvZfEm2LtZ9aehUVsxfx2qNT9yhXj0lXXJ1zd5lZGXAVMAPYDJzrnOs5s5vbJ36XmZ1LtGh9AjgAfMo5d/coxgT4DNHLt3cBOcDDwMXOuXBsjLCZnQfcRLRAdwB3AleMMpd/mdlbiT5z+zVgT+y/Nx3mIROZctKx3Fyy3ICPJdPyebG6Bedg074mTl+cmlmgJr3WWmjaQ7B+Dy0H9xBsrCLS0YTrbMTT2YSnqwV/dzOBUAt/KrmE+/yvp7mzm47uMF2hCDd2foUVLrHogTf2M5BvP9HCPZHeD4I8HLid+Z6qEaX6pZeW8FSk95WFTVl/pdBGdof3LQ+9wNLjs6ZGcQVwzt3EIMUmdnNP322PAScd7pix/UGikzlcPkTMHuDNw7zPSHL5C9EJKUTkMCQvlF6aouXm+lo5t5jWmh2c6tlG3YYDsPhjaXmfcReJQNNeuqq2sidczObQXHbXt1Pd3EFVUyeX7v8qp4aeJovoc4ZDqTqwn/Xh+l7bWv3ewatpHwHr7rctMopvKw3Xb5tjwKcaB+3v+g+REpOyuIqIDCUda7n29dbsZ/hO1mcA2LLzBCBDi2tzFa3bn6Bh2+P4q5+hpG0H2a6TAPB46I18M3Rxr/CzfUWcOsLKUGht/bbtdDMojLTThY8u56cLH90WIOTxE7YAYU+AiMdPyLKIFB3DiqxiPGZ4DMyMh9rfxdOuBcxweMAMzIPDMPP02j4nbzlv8lfgMYuGmfH7livwEcbhwXr6Wk9fD848sVgvr80/kSVzUjd1ZjIVVxHJOGfu+Slr/C/T4AqY1f0ZYF7K32P+8jXwdLS9KLiNzo72lD5PmzbBFlq23E/Ds/eSX/M0ZV0HyAfyBwhdYvv6bdvhZrIlMo8aymgOVNCRNY1wVjEuuwjLLsKbW4w/t4RAQTELC8q5M7+AgmwfuQEvAa+XgO8cAj4PAZ+HLJ8Hn8cYZI6cQb4rXzWGDw+pWNc3FVRcRSTjLGh7nqXeLQBs9nw0Le9RMXM++2wGs10V2dbN1uefYNmpb0jLe6XCztpWHtxaQ96/buT9rf+PgiFiG1w+L7s57MpeyusWV7JoWh6zinOYXpjNjKIzqCzOZmluAI9n5JdYpTcVVxHJOHnhxKT9uSWpvxmlR1XRScxujN4e0bT1EZhMxTUSoXXPM/yppoL/27CP5/dGb+xZZMt4f9KXpZ3Oz/NuMbvzTyQ06zUUL1zJ3DnzWV5ZwKl+LxcPMryMjYqriGScXpP2lww2M+rY2aLVsDFaXEsOPJa29xmV7g6qH/9/eP/5Y0qDB7gxeAPVlMV373CzuD98Ms2FS2DxWuYedxonzJvGKYER3mUkKaHiKiIZJRIOU+ha6LkptLA0PTekACxedQHhDV/Ca46jul/kUG0VJRUz0vZ+Qwq2svdv11H4/M+Y3vP4icGF3if4Sfit+L3GmqOnsXZZJSuP+Svl+cPd6yvppOIqIhml5dBBiiz6/ESzy6UwKydt71VcPp0XA0s5pnsrHnPsWP8nTj7/0rS934DC3Rx45GZy1/+AOZFDvXY1uVzKCnP4+hnLuGD5rLQ9liSjp+IqIhmlqb6antmEmz1FpH7a/t4OzVoNr8amA9/+ANEJ3sbHoc0P0nXPZ5nZtbfX9v2ujCfLL2LJmy7jw4vnjFs+MnIqriKSUdoOJdbwaPUWp/39pp30Fnj1pwAc0/wUne2tZOcO9GBL6oTbGnj1l59kUVXvuWaqXClPzPooJ59/Ge+aXpLWHGRsJtXE/SIiw+lsPBhvd/jTX1wXHreKfRb9njXfOtjy+N3D9BibzfubuOi2DeQcWB/f1uxy+F3px2n7+D9518e/ykIV1klPZ64iklG6W2rj7a6s0rS/n3k87Jj3Lp585Xn+GDmdooPHkI5F6MIRx0/XvcL1D20nFHF8w/NB/idwHY/4ziT//Gt4xwnL0vCuki4qriKSUSKtieIazikbIjJ1Zp/3RT74g+ijOIGXGmhs76I4hQsGHNixmcsfaGbj7sQNS495TuE3y+/gbW8+jyyfHqPJNCquIpJRnsk7g992Qak1c/yMs8blPRdV5HPC7CI27WuiKxzhtxv28bGzFo55XNfdyfZffYHFO39JdveXgeMBOGluMT9813Lml+eN+T1kYug7VxHJKC9HZnN35CxuDb+ZrlljnYd25N53SmK1y1+sf5VQODKm8Vr2bWXf90/nqF134DHHD/w3U+5p5fNrj+L//n2VCmuGU3EVkYxSn7TcXNk4Ptd5wfJZlOT6yaWT85t/zYY/33p4AznHroduxfuz1zIn+Ep88y7fIn7xoZO5/Jwl+Lz61ZzpdFlYRDJK8lqu6VpubiDZfi9fOqGdc579DBXWTPWzj9C59r1k5w41RX5v4Y5mtv/vxzmm9r74tqDzcf+sy1n7wavIydKv5COF/jwSkYzSMEFnrgBvPmcN3ti0i9Op4/k7vjTivrUvraf2+6/pVVh3MZMNa3/L+R//hgrrEUbFVUQyhotE+G3nx7kv8GV+6f8OpTnj+yssv6CYl4/7fPz1vx34FZsfvWvIPpFgO1t++QWKf30e08NV8e2P5ryerP94gtPPODtt+crEUXEVkYzR1nKI2VbHUs8eVnheITd7/Cenf83bPsWmQHRBbo85Fq67nC2P/t+AsRtebeDDtz7K7O134icMQIvL4b6jv81ZV/wfMyvKxy1vGV+6DiEiGaO5roqeiQebrIg8G//FvD1eDzM+fAcHbj6bmRwk14Ic+9jHeGHjL3BLXg/503jK9xru21wdX2P1Gu+7+Y7/f3nBczTurbfwphNWjHveMr5UXEUkY7Q2JM8rXDREZHpVTJ/Dznf/lgO/eScziU7HeHzr3+HZv9PmsrggeBsu6cLg73gdpxy1hLUXfoycLP9EpS3jSJeFRSRjdDbVxNvt/omdX3fhMcvx/vvD/DN3Ta/teRZkJvUA+L3GO1fO5v7PruH8916mwjqF6MxVRDJGV3Ni0v7xmFd4OJUz5lL5xT+x5bl/ULvhDwQOvYJFurlgSQXzjjqetcuma43VKUrFVUQyRri1Lt4OZU98ce1x7PJTYfmp8dfjN2+UTFa6LCwiGcPaEsWVXN1pK5OXiquIZAxvZ0O87clXcZXJS8VVRDJGIJgorv6iygnMRGRoKq4ikjFyQo3xdm7RtAnMRGRoY76hycxOBW4ASoCtwDM9P865/WMdX0Skx1W+zxFpq6LUWvjSjGMmOh2RQaXibuGfAI8DdwEPAS3AVbF9eqhLRFJmc0c5ra4YHHyvtGKi0xEZVCqK61HAKc65kJl1O+feb2abIPYUtYhICgRDYVqDIQB8HqMwR08SyuSViu9cm0icobabmR/4MfD5wbuIiIxOfWvvdVxtAuYVFhmpVBTXJ4A3xNpbgdMAA+amYGwREQDqGlvIIlpgK/LHfzUckdFIxXWVDwM9M2hfB/wGOAj8MwVji4gAENnxCC9l/zvNLodnu1YDZ050SiKDGnNxdc51AB2x9l/N7EJgBfCrsY4tItKjqzG60HihdZDr1yVhmdzGfFnYzO43s6U9r51zTznnfuKcOzTWsUVEeoSbE8vNhXJ0p7BMbqn4znUd8KSZ/djMylIwnohIP562xIo4VqDZmWRyG3Nxdc59F1gG5AIvmdlnzUz3yItISgU6auNtX9H0CcxEZHgpmf7QOVftnPsw8HrgrcBWMzs/FWOLiADkdCUenc8umTWBmYgMLyXF1cw8ZrYYmAHcCwSAP5jZi2b2KzP7YireR0SmroJQorgWlM+cwExEhpeKuYV/SvRxnDbgFWA78PPYf5uB44GTxvo+IjKFOUdJpDH6BD1QXDFnYvMRGUYqvht9F3Cyc+6FQfbfm4L3EJEprL21kVwLAtDhAhQWFU9wRiJDG/aysJl9xMz+ZmZPmdl1Ztb3ekz5EIVVRGTMGmv2xdsNnhLMo9UyZXIb8szVzD4G3JK06VTg3WZ2unNuJ4BzzqUxPxERmhuq6PmrvsVbOqG5iIzEcH/+XQbsBVYBc4CPEX3k5gfpTMrMLjOzXWbWaWYbzWzIec7MbHUsrtPMdprZpaMd08yyzOxGM6szszYzu8fMZveJmWtm98b215nZDWYWGG0uSbFXmpkzsx+P7MiITE2v5h7PUZ2/YFXnjdw+/cqJTkdkWMMV10XAjc65fzrn9jvnbgO+AZxnZnnpSMjMLgJ+BFxNdBrFp4D7zGzAhQDMbAHw11jcCuC7wI2xaRhHM+b1wIXAe4hOWloI/NnMvLExvMBfgILY/vcA7yDpD42R5JIUeyrwcWDTyI+OyNRU2xKkCz9VlOFKF0x0OiLDGq645gMH+mz7G9HLyUenJSP4HPBz59ytzrltzrnLgSrgE4PEXwoccM5dHou/FfgF8IWRjmlmRcBHgCuccw86554BPgCcALwuNsbrgWOBDzjnnnHOPQh8EfiYmRWOIpee97uT6F3WmiZSZBi1ScvNaUUcyQSHc1dAzxxkhUNGHYbYJdaVwAN9dj1AdCm7gawaIP5+4GQz849wzJVE16SNxzjn9gLbkmJWAdti25PfJyvWf9hckrb9D/A759yjg3ymODNbZ2brgOXDxYocqWpbgvF2RYGKq0x+Iymu55nZm8ys70zZ6bhdrxzwAjV9ttcAg813Nn2QeF9svJGMOR0IA3XDxPQdoy7Wb6iY5Fx6bhJbDFw1yOcRkT6y6reyyPZTSCsV+YHhO4hMsJE85/oe4N0AZrYX2Aw4YJmZPavVb0bOzI4m+r3vGc657pH0cc6tifVdB6xOW3Iik9gHar7HN7J2ALAt+HtAMzTJ5Dbc2WcR8Fqi3xn+GmgH3kh0npQfAXVmttfM/mJm3zWzd48xn54zwb5LXlQC1f3DIbZ9oPhQbLyRjFlN9Oy2fJiYvmP0nBUPFZOcy6pYny1mFjKzENGCeVnsta53iQygKNQQbxeUzx4iUmRyGLK4OudanHOPOeeuc8693zm3jGjBPRP4LHAH0Ej0Zp8vEb1J57A557qAjYpVmDkAACAASURBVMDaPrvWEr0DdyDrB4nf4JzrHuGYG4Hu5JjYYzhLk2LWA0v7PJ6zFgjG+g+bC/BHotNBLk/62QD8JtbuQkR6iYTDFLum+OvySk19KJPfqKc/dM61AU/GfgAwsxzgRFIzh/APgTvM7OnYe1xK9BrQzbH3uj2Wx8Wx+JuBT5rZ9UQnvDgduITo5ewRjemcazKz24BrzOwgUB/rswl4KDbGA8AW4HYz+zxQBlwL3Oqcax5JLs65RqJ/jMSZWRvQ4JzbfHiHS+TI1lC7n3KLANBIPsU5uROckcjwUrLuqnOuA/hH7GesY90VW3T9KqKr7GwGznXO7Y6FzO0Tv8vMzgWuI/pozQHgU865u0cxJsBniF6+vQvIAR4GLnbOhWNjhM3sPOAmogW6g+iZ+hWjyUVERqexenf8+5oGTzmaVVgywaRc1Nw5dxPRIjbQvjUDbHuMYc6ahxoztj8IXB77GSxmD/DmYd5n2Fz6xK8ZaazIVNRWl3j6rSXQ97YIkclJs1+LyKQWbEgU187svvcLikxOKq4iMqm5psQkcZH8GROYicjIqbiKyKTmbUs8hWdFer5VMoOKq4iMi1B7E9t/+m52X3sG9Ts2Dt8hJrvzYLydVapnXCUzqLiKyLjY/IdrWFJzH/PaXqDj15fACJeCrg/nUu1KCDujoGLAxbFEJp1JebewiBx58l59MN6eHdpD7atbqFhw3LD9PhX6NE3BbryE+cd8rV8hmUFnriKSdi4SYUb3nl7bqnc8N2y/zu4wTR3RabjN46MsPzst+YmkmoqriKTdoZo95NPRa1tn1bZh+1U3dcbb0wqy8Hgs5bmJpIOKq4ik3b7uPF4XvIa/hl8T3+at3z5sv+rmRHGtLNJZq2QOFVcRSbua1givuNn8KnxOfFt2x2ALXSV07H2ed3kfZbXneZblNA4bLzJZ6IYmEUm7njPQelcY35bXPfxS0Lm7H+Ea/60APNXxbuC8tOQnkmoqriKSdgdjxbXWFfNCZD4NrpAazyzmDdPPWqoSLwo0gYRkDhVXEUm75pZmPESoo4i3dF0NgD9ivNM5zAa/SSnQXhNv+0tmpT1PkVTRd64iknYX7PoWr2R9gOeyPsbZnmcA6A47WoKhIfvlBhOzM+WUaXYmyRwqriKSdoHuJjzmKLY2uvDHtze0dg3ZrzhUF28XVQ53EVlk8lBxFZG0ywk1x9tN5MXb9W2DF9dIKESpS9z0VDZdUx9K5tB3riKSdnnhlnh7TUEVZ7RuptSaCe3NgnmvG7BPY91+Si0CQAMFlObmDRgnMhmpuIpI2hW4Fojdt7Tat5mT/Y8BsKF6JTBwcT1UvZvSnranPN4WyQS6LCwiadXd1UmeRR/FCTkPXQWJy7uh1rrButF8MDEXcUugIn0JiqSBiquIpFXzodpE2/IhL1Eorb1+0H5dDXsT7dzK9CQnkia6LCwiadXeVE9ZT9vy8BaUx/f5Ogcvrnsi5dSFX8NMa6Cj+Og0ZymSWiquIpJWnW2JO4U7PbkECqfFX2d1DT4F4hN2Mvd0R2dlunbZCelLUCQNdFlYRNIq2JaYcD/ozSOnOHGJN3eI+YWrmhJL1M0szklPciJpouIqImnV3Z44c+325VFQkiiuBeGmQfsdaEwsNzddy81JhtFlYRFJq11lZ/LBzv+hgA7OWTSDL1XMiO8rck24SATz9P47Pxxx1CSt5TqzSGeukllUXEUkrdq6HM3k00w+4fyZ5OYV0O6yyLUgAQvT1tJIXlHvp1gbag/wHc8tVHtKqfbPJiegpeYks6i4ikhaJU/On5/tw8xoskJyiT6i01xf1a+4Htq/nYt86wDY4VkAfGu80hVJCRVXEUmrtuTiGoj+ynky+yw6Wxupp5Bzuv3M6NOn9eCr8XZLQM+4SuZRcRWRtLKWKmbbQVpdDgWB6LZ7pl3K443RM9fjXVG/PskTSARzp49LniKppOIqImm1eu9P+ULW/QD8q+5bwBLK8gLx/fUDLTvXtD/ejBRqkXTJPHoUR0TSyhdqjbe9OYUAlCYV14YBlp0LtB2It/0lc9KYnUh6qLiKSFr5Q22Jdm70EvBwxTWvsybezinXIumSeXRZWETSKhBuj7f9udEz18Wh7Vzt+xml1kLk1eXA93v1KQkdjLeLpy8YlzxFUknFVUTSKiuSKK45ecUATKOBN/geAWBTi/WKD3V3UeYOgUHEGeUzdeYqmUeXhUUkrXKSimt2QfSycPL8wjl95heuq9qN1xwADVZEVnbuOGQpkloqriKSVrkuMQF/Tn70zDWvJPF4TX64sVd8Y/WueLvBNw2RTKTLwiKSNpFwhFwScwTn5UfPXAvLk+YXjjT36rM3XMZvu9/HTGugsHQuR41PqiIppeIqImnT3t5MfuwSb4cLkOPzA1BYWEKX8xGwELkWpLO9hezcAgB2dRdzWzg6l/Al8+fzzolJXWRMdFlYRNKmoyWxpFy7JVa2MY+HRiuMv26sq4q3k5eam6Gl5iRDqbiKSNp0tLfS7HIIO6PDei8b1+wtSbTrEpNGJC+SPkOLpEuG0mVhEUmbpuzZnBW8DXCcOD2HPyXta/OXQngHAO2HquPbq5qS13HVmatkJhVXEUmb1viKONbvkZpgVhk99zp1NyaK6zfrPk+L30eVK2Nm7mvGKVOR1FJxFZG0SV5uLi/L22tfOKcMYl/JhlqjMzIFO9tZzovghbAzXEn/FXNEMoGKq4ikTVtXcnHt/eumbsZq/ntvmDqKWJhzOquA+gOvMrNnv5VS6Q8gkokm5Q1NZnaZme0ys04z22hmZw4TvzoW12lmO83s0tGOaWZZZnajmdWZWZuZ3WNms/vEzDWze2P768zsBjML9IkZMhczu9LM/mVmzWZWGxvvuNEfJZHJzzXuYZVnCyfYDmZ4ek8W0TX7dG4On8/vwqt5MRwtqfUHdsb3H/JrAgnJXJOuuJrZRcCPgKuBFcBTwH1mNneQ+AXAX2NxK4DvAjea2YWjHPN64ELgPcCZQCHwZzPzxsbwAn8BCmL73wO8A/jBaHIB1gA3AacBZwMh4CEzKx3FYRLJCJX77ufXge9wT9bXOKfhN732lRdkxdv1bUEAWmt2xLd15MxAJFNNxsvCnwN+7py7Nfb6cjN7I/AJ4MoB4i8FDjjnLo+93mZmpwBfAO4eyZhmVgR8BPiQc+5BADP7ALAbeB1wP/B64FhgnnNubyzmi8DPzOyrzrnmkeTinHtDcvKx92kCTgfuHeWxEpncgom1XAnk99pVnp+46FPXEl12ztUlimtXkVbDkcw1qc5cY5dYVwIP9Nn1ANEzvYGsGiD+fuBkM/OPcMyVgD85JlZAtyXFrAK29RTWpPfJivUfNpdB8i8g+r/DoYF2mtk6M1sHLB+kv8jk1ZVUXLN6F9eKpDPX6ubobcOB5lfj2/wVi9Oamkg6TbYz13LAC9T02V5D9AxyINOBhwaI98XGsxGMOR0IA3UDxExPiuk7Rl2sX3LMULlU0d+PgOeA9QPsE8lonqTi6sku6LWvPC+LmwI3MI8qZkTqaWvZSnFn4m/X/JmaVVgy12QrrlOKmf0QOAM4wzkXHijGObcmFrsOWD1uyYmkgDfUlmj3Ka4ej7Hcu4uZLvo366v7dzAjdCD65zAwbd6ycctTJNUm1WVhEmeClX22VwLV/cMhtn2g+FBsvJGMWU307LZ8mJi+Y/ScaQ8Vk5xLnJldR/SmqLOdczsROQL5ehXXwn77G5PuCD744j/Is+jl4RZyKCqb3i9eJFNMquLqnOsCNgJr++xaS/QO3IGsHyR+g3Oue4RjbgS6k2Nij+EsTYpZDyzt83jOWiAY6z9sLklj/4hEYX1xkM8lkvECocRC6f7c/sW1PSdRQN2ux+Ptau8szDOpfj2JjMpkvCz8Q+AOM3saeJLoHbgzgZsBzOx2AOfcxbH4m4FPmtn1wC1E77q9hGjxGtGYzrkmM7sNuMbMDgL1sT6bSHyH+gCwBbjdzD4PlAHXArfG7hQeUS5m9hPgA8BbgUNm1vPbpdU5l3T3h0jmC0QSxTVrgOIaKpgdn6Vp2yHjm93fYYFVc/LCaSwZryRF0mDSFVfn3F1mVgZcBcwANgPnOud2x0Lm9onfZWbnAtcRfbTmAPAp59zdoxgT4DNEL9/eBeQADwMX93wX6pwLm9l5RJ9RfRLoAO4ErhhNLsBlsf8+3Oej/xfwjREdJJEMkRVJrHATyOtfXP2VR8O+aHuGO8gWt4AtbgEnHaPvWyWzTbriCuCcu4loERto35oBtj0GnHS4Y8b2B4HLYz+DxewB3jzM+wyZi3POhuovciTJcYkz15z84n77S+edEP9SZYnti29fNrN/IRbJJJOyuIrIkWG3q6TVZZFnHeTk9y+YM5ecQMQZHnPMsxoKaaPV8lRcJeOpuIpIWoTCES4Mfj3+emd+Sb+YrJwCXvEvZnFoO3+PHM8cqyV33mwKswebc0UkM+h2PBFJi/buxKPbeQEvHs/A34g0zl4DwGrvJv43cA0fqdw+HumJpJWKq4ikRe+1XAe/SLb0gi9w0FMBQMSXy+rXnpv23ETSTZeFRSQtkotr/hDFNa9kOtmf/xc1Lz5JxTGn48vrf/lYJNOouIpIWnQequIi76O0uWxyPXOJrrY4MG9eCZUrh7wRXySjqLiKSHrUvsj3/NFVHrd2Hg98fGLzERlH+s5VRNKiq7053u725k1gJiLjT8VVRNIi3NESb4f8Kq4ytai4ikhaRDoTxTXsU3GVqUXFVUTSIhJMFNdIIH8CMxEZfyquIpIewaRFnlRcZYpRcRWRtLCupOKapeIqU4uKq4ikhac7UVw9WQUTmInI+FNxFZG08Iba4m1PjoqrTC0qriKSFr5QYi1Xv4qrTDGaoUlE0uIV70Kawx3kWyf+/MqJTkdkXKm4ikha3Op/Py92Rx/H+fOskyY4G5HxpcvCIpIWrSNcFUfkSKTiKiJpMdL1XEWORCquIpJy4YijqaM7/ro41z+B2YiMP/05KSIp19TYwJe9d9LgCmj0V+D3njfRKYmMKxVXERnQ7le2UH33l+kMlLDiYzdRmD/yWZZaavfycd9fANhv04FvpylLkclJl4VFZEAP/+kOFrU/x+qmP/HCr/9zVH3bGw/G263eolSnJjLpqbiKSD/Nnd38tPY4Lgh+i+93v5PpVY+Mqn9HUnHt9JekOj2RSU+XhUWkn+01rQTx86/s/wCgPZyFi0Qwz+B/jz//yG/wb7iV7gXn0O3JiW/vylJxlalHxVVE+nnlYAvN5NHscii0DnItSFN9FUUVswaMr9q5meMeuxSvOdiygZcCx8X3hbNLxyttkUlDl4VFpJ89De2Asc9Ni2+r3fvK4PFP3BktrEDQ+ZkXfCmxs2B6utIUmbRUXEWkn6J96zjf8xTT7FB8W0v14MU1u3pjvH1l90f4R2RZYt+0RelJUmQSU3EVkX5W1fyaGwI/ptya49u66l8dNH5mx8vx9vNuEXMscUNT0cwlaclRZDJTcRWRfkq6a/pts8bdA8a2tTRSQfQMt8t52e/KexXXynlHpSdJkUlMxVVEeomEw1RE6vptz2nbP2B8ze4X4+29bhoFtPPv3Z9jfXgZ+2wGufnFactVZLJScRWRXhpq95Nl3f22FwcPANDS2kpjY0N8e9P+xM1LzTmzCedVUu8KWeXdyt5pa9Ker8hkpOIqIr3U70t8f7rPEnf6TovU8swjd+O7dgG+65bx7MO/BSB4cEc8pqtgLj9+7wouKdzAxtwzWPqub45f4iKTiJ5zFZFeWqt3xdsHcxZzQ/tF7AgWs89V8JvHvkKOdQFQ/OQ3cWe/o9d3sa50AactKocrbx/3vEUmE525ikgvXfWJ4tqVP5stpWvZ6I4m3zpYYNXxfQsie6h69SW6O1oJOwMgq3zBuOcrMhnpzFVEevE07U28KJnH/EgeWw40s9qzqVfcs5HFdB6o4j89l7MneDEzrJ6fHX32OGcrMjnpzFVEeslt2xdvZ1csYF5ZLgCnezbHt1/V/SHe1vVNNgZns+9QByF87HWVzJo+rd94IlORzlxFpJeSrqp4u2jGIua35+EnxDneZ+PbH4+cAMAT2+voCkei/XL95GfpV4oIqLiKHFGcc+DckKvX9LX7xWep3vI4npwi5vzbefxf9+mcZNs53rOTynlHs6K6mm1Zl8Tj90Yq2OMqAfjnrsQjOXNKc1P2OUQynYqryBGio7WJ3Teey+LOrWw5+nJOfO/wj8FsuPcWlm/4MvMsQpvLYs3jfmp5OwDHVObzt/wiFs7Po86KqSRaSHeXngax+5oW2gHmWg173TQWFmv1G5EeKq4iR4hNf72FU4KbwWDpSz+mpekzFBQNXvCa6w9y9Mav47PoZd2fhc+jlsRsSicviPb1+nzUrPk+tu4LNPrKOfqd/8UpP3+MyraXuSHwk3j8P9o+AKxKz4cTyTAqriIZZsfWZwjefSnO42f2h2+naEZ01Rn/rkfiMQELs+nJP+KWXsD2g22ctqiM+eV5vcbZes8POJUOANpcFvfmXAAt0X1ej3HRyXPjsSesuRDWXMg058CMr3o/ygmBDb3G85TOS8fHFclIKq4ik4wLBTFf1qD72//4WY4PvwRh2Pz7/6ToP+6kMxhkSftzYIm4k5/+LF978nnuCL+evICX3156GstmFgLQ2dHGkt2/jsduPem/+PO5b+HnT73Kc3saecfK2Rw/u6j/m1v0DdqKlkBn7+Kao6XlROIm5aM4ZnaZme0ys04z22hmZw4TvzoW12lmO83s0tGOaWZZZnajmdWZWZuZ3WNms/vEzDWze2P768zsBjMLpDoXOTKFYnfVJtu/cxtP/fQTPP/o74iEIzz7w7fjvlXJv26+lObObiKR6ALkkXCYFx7/A0/e9kWO73ou3v+42j/zl29fyMXfuJEC6+g3/ud8v6OIVtq6wnzp7k3xHF6472eU0QRADWUsf9OHyfZ7uXT1Im7+wEpet6xyyM/inX5sv23Fs7S0nEiPSXfmamYXAT8CLgP+HvvvfWa2zDm3Z4D4BcBfgf8F3g+cAdxkZrXOubtHMeb1wAXAe4B64IfAn81spXMubGZe4C+xfWcCZcAviJ4rXJ7iXFJue00L3o56PMEmuosX9t7Z1Yq/rWrgjoBzSW1vFt2Fc3vtt+AhfG0HGYjr8zriyyNUMKvX2N6OOnwddTjz0OvUC3B9RohkFRLKreyVl6+9Bm+wOd7DIiFwYcAwiwaFukPg8RLKmUZnTgURB3kBH93hCFkN28gKd9DeFcK6WgmFusgpnUVHUy3OfIQiUJibRV0oi13ehQR8Hjq6w8wvy6WzqZa8527DeXzR4+v1YfNOw+vPIvLML6lsf5lpro6fhi/gDncuJWUVvGPlbJbabvKevoGTWx9lFkDNr+Cxj7CC6CHIP/B33vSNO8mdtoDPrj2K/Ps/y1mt9w14jM8LPcRBr5cruz/CaZ4tvMX7j/i+EmvlLM8m7o2cxgv7m/jRHx7nTbODzHrhx/GYnYveT2Vg8DPlgUw75jR4PvG6mVxmzj9mVGOIHMnMub6//iaWmf0T2OSc+1jStu3A75xzVw4Q/z3g7c65JUnbfgYc65xbNZIxzawIqAU+5Jy7M7Z/DrAbeJNz7n4zexPR4jrPObc3FvN+4GfANOdccypyGeDzrYs1l69evbpo3bp1fUNGZMGVf+Fyz+95n+8hTgne1GvfOZ6N3Bb4wYjG2RRZwPld3+m17b3eh7naf9uI+j8UXsFHu6/ote2zvt/xad/vR9T/V6Gz+Uroo722Xe27lff6Hh1R/5+Ezufa0Lt7bbsr8E1O8bw4SI+EJpfLquCPaSc7vm0WtTyZ/ekh+0Wc8Znuy7gncnp82xrPc/w8cM2gffa5cs4I3hB/vdJe4s7A1WQPsFrND7rfwS3ht9CFH4Abz3Ks2HAFhuPgqq/xlO8Urr0/unLND/038Xbv3+N921w2oU9voai0fJhP35uLRKj61jHMdNF1XzcWvJaVn//jqMYQyWRr1qzhsccee8w5t2ag/ZPqsnDsEutK4IE+ux4AThuk26oB4u8HTjYz/wjHXAn4k2NiBXRbUswqYFtPYU16n6xY/1TlImnU7rL7bWsbYNtAiqydi7wjK+LJPOa4zn8Tl3rvoedcfl3kRJ6OHD1on9lWR/J5/0Z3NJd1f5qdnvnsttn8s/g8NhSczcbIEu6NrKILP2Zw3UUn8pZz38zs/9zGrP98kRVr38fHz1rIaYvKKKOJ8zz/6PU+m5d9btSFFcA8HhrOvpY6itnlmUflBd8e9RgiR7LJdlm4HPACNX221wCvG6TPdOChAeJ9sfFsBGNOB8JA3xWia2L7emL6jlEX65ccM9Zceun5qyh2Brt6oJiRWFyRj6+rnJrumSwuyu+1rzhUwp7O2YP07K3JN4OjSnr3z+sqZXdwzoj6dwSmc0xZQa9tvuA0dgfnYP0uIvfn8itZmlMYf20AHZXs6UrkH8FL2KJ/NzpnYGAYhsObP4cTc4oxoL0rRMDnoabtaLZ2d2I4PAbdnizyQs10mw8PjqBl43ddhM3HTG8XywoLyQl4cc5RglF9sJx2Tx5ZdNHkLQOPl0C4jcb8xcxpfJpmXxnt/lKKs2bxkUUL2F7bRn1rkPUlX8T8LxPJKqRo652UdVdxIG8pc9/xXdqDIdaVHc3//P1VNu1rJNfv47x/+xALV34NgJ77cmuaO1n5t5dYm+fnijccQ8DX/+9lv9fDzz/0Gh58qIGXnz2B3FAjjTlzyTr1I5xyxvkj+t9tIMedeQGceQFlkcioJq0QmQomW3GVNHnwc6vpqc19q390++UjGmcu/U+7o/2vGlH/ecBb+m09C/jeiPq/L/bT28jvB/tc7Odw+y8FPtZv67nx1qx++6Dn1qAThxr4bZ8EoCL2siT236vfdvyQ+VQWZvODdw05MgABn4fz3nguvPHcYWNHS4VVpL/J9v+KnjPBvrcqVhKfE6af6kHiQ7HxRjJmNdEzyr7Xx/rG9B2j50x7qJjR5iIiIhluUhVX51wXsBFY22fXWuCpQbqtHyR+g3Oue4RjbgS6k2Nij+EsTYpZDyzt83jOWiAY65+qXEREJMNNxsvCPwTuMLOngSeBS4GZwM0AZnY7gHPu4lj8zcAnzex64BbgdOASoo/UjGhM51yTmd0GXGNmB0k8irOJxFXUB4AtwO1m9nmij+JcC9zqnOt5DmTMuYiISOabdMXVOXeXmZUR/RJvBrAZONc5tzsWMrdP/C4zOxe4DvgEcAD4VM9zpSMcE+AzRC/f3gXkAA8DFzvnwrExwmZ2HnAT0aLYAdwJXJH0PqnKRUREMtike85VBmZm61avXr36cJ9zFRGR1Mmo51xFRESOBCquIiIiKabiKiIikmL6zjVDmNm+oqKiWcuXL5/oVEREprznnnuOpqam/c65Aae3U3HNEGb2LDCH6N3FPcqIPjaU3O6pvs+ROsnvM9bYofYPtK/vtqFe6xjoGCS3dQwy9xiMdPtEHoMSoNY5t2LACOecfjLkB/ifwV73tIF1wLp0vu9YYofaP9C+oT6zjoGOgY7BkXkMRrp9Mh8DfeeaWe4d4nXffel837HEDrV/oH1Dfea+r3UMdAxG8t5joWMwPsdgpNsn7THQZeEjTM/6r26QZ6+mAh0DHQPQMQAdA5i4Y6DiKiIikmK6LCwiIpJiKq4iIiIppuIqIiKSYiquIiIiKabiKiIikmIqrlOMmS0ws0fNbKuZvWBmeROd03gys6PN7Lmknw4ze+tE5zWezOyzZrYl9m/gBjOzic5pvJnZF2LHYLOZvX+i8xkvZvYHMztkZr8bYN+bzewlM9tuZh+diPzGwzDHYNB9o34fPYoztZjZY8BVzrknzKwUaHbOhSY6r4lgZvnAq8A851zbBKczLsysAvgHcCzQDTwOfME5t35CExtHZnY88AvgNMCAR4E3OucaJzSxcWBma4AC4IPOuXckbfcBW4HXAk3ARuA059xIpznMGIMdg+H2jZbOXKcQMzsW6HbOPQHgnGuYqoU15nzg4alSWJP4gGzAH/s5OLHpjLulwHrnXKdzrgN4HnjjBOc0Lpxz64CWAXa9BtjinNvvnGsF7gNeP565jZchjsGQ+0ZLxTWDmNlZZnaPme03M2dmlwwQc5mZ7TKzTjPbaGZnJu1eArSa2b1m9oyZfWXckk+RFByDZO8C7kprwik21s/vnKsFvg/sAQ4ADznndozbB0iBFPwb2AysMbNiMysB1gCzxif7w5fif/t9zQT2J73ezyQ8Jmk+Biml4ppZ8on+Yvg00NF3p5ldBPwIuBpYATwF3Gdmc2MhPuBM4DJgFbDWzNaOQ96pNNZj0BNXSPSy4F/TnXCKjenzx4rJm4H5RH95nmZmZ41L5qkzpmPgnNsK3AA8Avye6GXy8LhkPjYp+bef4TLnGKRypQD9jN8P0Apc0mfbP4Fb+2zbDnw31l4F3J+07wrgion+LON5DJK2fQD45UR/hgn4N/BO4Cd9/g18caI/y0T8G0ja9zPgvIn+LOP1uYmeqf+uz7bTgD8kvb4eeO9Ef87xPAYj2TeaH525HiHMLACsBB7os+sBov/HAfgXMM3MSszMA5wFbBu/LNNrhMegR8ZdEh7OCD//XqJnq9lm5iX6i+SlcUsyzUb6b8DMpsX+ezTR7xvvH68c02GU//YH8jRwnJnNit3o9yYy7Jik4BiklIrrkaMc8AI1fbbXANMBXPTmpa8QvUN0E7DdOffn8UwyzYY9BgD/v727CZmqiuM4/v2pYZoYmqW5kBYWpkbmE5KRQfaitCgpoRa+YKvQItJFkItWQVEKUkQtNE1xYwlaoaZYEGIpqZBGEWYvkklmGkFq5b/FuWPjcGecZ547Mz7P/D5wuc8z595zzj0M859z7rlnJF1NH/hAzVHPe+Az0lD4ftJ74DCwuYV1G0Bu3AAABPtJREFUbLa63gPAJklfAeuABdH7J/bV+97fAWwAHpR0VNJUuPDZsIQ0c/oAsCx630zhHrXBpdK6a0CjJ1rvFBFbSDMBO1ZEnAZGtrse7RIRS4Gl7a5HO0VEwx+avVlE3FcjbTN964tWrku0QdW07nLPte84QZqUURk0RgK/tL46bdHpbdDp1w+d2wadet3lLqs2cHDtIyLiHOnB78rZv/eTZsz1eZ3eBp1+/dC5bdCp113ucmsDDwv3ItlEg7HZv/2AMZImAScj4kdgObBW0h5gF/Ak6fm1N9tR32bo9Dbo9OuHzm2DTr3ucr2qDdo9ndpb/RtpZmfkbKvLjllIWtLvLOlb3N3trrfbwNfvNvB1d1obeG1hMzOzgvmeq5mZWcEcXM3MzArm4GpmZlYwB1czM7OCObiamZkVzMHVzMysYA6uZmZmBXNwNTMzK5iDq5mZWcEcXM2sKSRNlPSPpMqF1JtV3sOSzkm6sRXlmdXi4GpmF5E0VNJ5SVFju6uOrJYDuyJiewN12JCVM6nGMZJ0RNIpSYMiYhPwJfByd8szK5p/FcfMKk0GBKwHtlQ5Zm+tDCRNJf3U16wG67ASmA0sAJ6pcsw9wA3AWxHxV/baCmCNpAkRcajBss16zAv3m9lFJC0GlgEzIuKjBvNYC8wERkfE3w2c34/0yyaDszzOVSljDjAlIvZmrw0BjgOrIuLpRupuVgQPC5tZpS7Sz3jV7J1WI2kAqce6Iy+wShoo6XlJhySdyYZ135d0W+mYiDgPrAauAR7KyWMo8ChwsBRYs/P+BD4l9XrN2sbB1cwqTQZ+APpLGlG51XF+FzAE2FOZIOkKYCvwArAbeBZ4CRgP7JJ0e9nhb5OC/IKcMh4HBpGGjyvtBkZJGldHXc2awvdczeyCbFj1JtIX719zDjkGjL5ENuOz/eGctKdIP3g9MyK2lZX7BnAQeDVLJyKOSPoYmCHp+og4VpbPAuAcsC6njFK5E4CvL1FXs6ZwcDWzcpNIgXUF8EFO+u915HFttj+ZkzaHFPC+yOkFbwfmZzN/SxOUVgLTgXlks4CzHukdwLsRcSKnjN+y/XV11NWsKRxczaxcV7bfHBE7G8yjNEtSOWk3k4Zz83rFJSOAn7K/NwKnSD3V0iM2T2T7VVXOL5Xr2ZrWNg6uZlZucrbvyWMspcA5PCdNpGdRF9dxPhFxRtJ6YKGkO4HPgbnAUWBblfNL5dYK4GZN5eBqZuW6gBMRcbwHeRzM9nkrJX1LGjbemc0IrsdKYCGp9zocGAW8WOP8sRX1MGs5zxY2MwAkDQbG0bNeK8B+4A/SfdFK75CCY27PVdLIytciYh9wAHgMWEQa7q02JExW7vGI+KZ71TYrjnuuZlZyK9AfQNKcKsd8GBE1JzVFxL+SNgKzJA2MiLNlyStIKze9Imk6sJMUiMcA9wJnSCsvVVoJvEZamOKTiPgur+xstvM0agdfs6bzCk1mBoCkRcDrNQ4JYFhEnK4jrymk+6OzI+K9irQBpGHeufz/2M7PpOdi1+StCiVpWHbMlcC8iFhbpdz5pMUnbokIDwtb2zi4mllTSNoKXBUR01pY5j7g+4h4pFVlmuXxPVcza5YlwFRJD7SiMEmzgInAc60oz6wW91zNzMwK5p6rmZlZwRxczczMCubgamZmVjAHVzMzs4I5uJqZmRXMwdXMzKxgDq5mZmYF+w/0QxLWAK3cJQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 460.8x345.6 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import alpro\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt \n",
    "import numpy as np \n",
    "alpro.util.set_default_plot_params()\n",
    "\n",
    "s1 = alpro.Survival(\"1275b\")\n",
    "s1.init_model()\n",
    "s1.set_params(1e-14 * 1e-9, 1e-9)\n",
    "s1.set_coherence_r0(None)\n",
    "s1.domain.create_box_array(100.0, 0, 10.0, r0=0)\n",
    "\n",
    "s2 = alpro.Survival(\"1275b\")\n",
    "s2.init_model()\n",
    "s2.set_params(1e-14 * 1e-9, 1e-9)\n",
    "s2.set_coherence_r0(None)\n",
    "s2.domain.create_box_array(100.0, 0, 10.0, r0=0)\n",
    "\n",
    "E, P2 = s1.propagate_fourier(s1.domain, pol=\"both\", mode=\"auto\", N=100000, f_res = 2000)\n",
    "_ = s1.default_plot()\n",
    "\n",
    "P1, _ = s2.propagate(s2.domain, energies=E, pol=\"both\")\n",
    "_ = plt.plot(E, P1, ls=\"--\", c=\"C1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Massless ALPs\n",
    "\n",
    "This example sets up a cell model with uniform cell size and and finds the conversion probability. The massless formalism is specified automatically by setting the mass to a value below the minimum $\\omega_{\\rm pl}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">Warning: I've found the massless scheme implementation in ALPro to be a little less stable than the massive one. This is not surprising giving the remapping from $z\\to\\varphi$, but users should exercise caution, particularly at the high $E$ (low $\\lambda$) end.</div> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAEdCAYAAAB67qLTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd5xcVf3/8dfnTt2a3gsJBAgE0mgJRUBBEBVQVKqIggoIgsLPylcR/IKNIiAE+SK9Se8dEnpLIaSQ3nvZvlPvPb8/7uzcmdmZ3dnNhp0kn+fjsY/cuXPumTuhvPece4oYY1BKKaV2VlZ334BSSim1PWnQKaWU2qlp0CmllNqpadAppZTaqWnQKaWU2qlp0CmllNqpadAppZTaqZVk0InIhSKyTESiIjJdRI5op/yRqXJREVkqIud3pE4R6S0iN4vI5yISEZFVInKbiPTJqaOXiNwnInWpn/tEpGfXfXOllFJdreSCTkROBf4JXANMAN4DXhSR4QXKjwReSJWbAFwL3Cwip3SgzsHAEOBXwP7AWcCXgIdyPu5BYCJwfOpnInDftn1jpZRS25OU2sooIvIhMNsY8+OMc4uAx4wxv81T/q/At40xe2ac+z9gjDFmcmfqTL1/AvAc0NMYUy8i+wDzgMONMe+myhwOvA2MNsYsaOM7zQT6AYuL/XtQSinFKGCTMWbCtlTi76Kb6RIiEgQOAP6R89YrwKEFLpucej/Ty8APRCQASCfqBKgGYkBzxuc04rYGW7wLNKXqaRV0IjI1dbhvjx49guPHjx/SxucppZTKMGvWLOrq6ra5npIKOqAv4AM25JzfABxT4JqBwGt5yvtT9UlH60w9d7sauMMYk8z4nE0mowlsjDEisjH1Xlsi48ePD06dOrWdYkoppVocddRRTJs2bZt7wkruGV13E5FK4FlgDe4zu04zxhxljDkKmNUFt6aUUqoTSq1FtxmwgQE55wcA6wtcs75A+WSqPim2zlTIvZB6+Q1jTDTnc/qJiLS06kREgP5t3JtSSqluVlItOmNMHJgOHJvz1rFkPxvL9H6B8p8YYxLF1ikiVcBLuN2cJxhjGvN8TiXus7oWk4GKNu5NKaVUNyu1Fh3A9cB9IvIR7mCP83GH/08BEJF7AYwxZ6fKTwEuEpEbgduBw4BzgNM7UGcV7uCUauBkoEJEKlLXbjXGxI0x80XkJeB2EflJ6r3bgefaGnGplFKqe5Vc0BljHklN1L4CGATMwW1hrUgVGZ5TfllqKsANwAXAWuDnxpjHO1DnAcCk1PHCnFs6GpiaOj4DuBl3VCfAM8BFnf+2SimltreSCzoAY8ytwK0F3jsqz7lpuJO3O1vnVNxnee3dVw3uZHKllFI7iJJ6RqeU2kFsnA/3nAgv/gYcp7vvRqk2adAppTrunRtg2TT48DaYqavgqdKmQaeU6rjZj3jHb+UuOqRUadGgU0ptm7qV4NjdfRdKFaRBp5TquJ45m4msndk996FUETTolFIdNz5n8HHtivzllCoBGnRKqY476tcw/kzvdby5cFmlupkGnVKqc4IV3nG8qfvuQ6l2lOSEcaXUDmDwRNjvFDfw+u3V3XejVEEadEqpzhl/uvujVInToFNKdczmxfDm/0K4GvqNhkkXdPcdKdUmDTqlVMfUr4a5T7jHI47QoFMlTwejKKU6JtbgHYequu8+lCqSBp1SqmOi9d7xghfg2UvhvZu7736UaocGnVKqYzJbdADT74IFL3XPvShVBA06pVTHxOpbn4s3fvH3oVSRNOiUUh2TN+h0wrgqXRp0SqmOieYJuoQuAaZKlwadUqpjcp/RgXZdqpKmQaeU6phEpPU57bpUJUyDTinVMclo63NOEpLxL/5elCqCroyilOqYA34AI78EdhymXuudjzeCv3f33ZdSBWjQKaU6Zsy3vOMZ90L9Gvc43gTlGnSq9GjXpVKq83RPOrUD0BadUqrzJnwfIjVu4GlrTpUoDTqlVOcd9vPuvgOl2qVBp5TqmIfPdOfS+cPwrSnaklMlT4NOKdUxK951uysBjOnee1GqCDoYRSnVMZnz5fzB7rsPpYqkLTqlVMdkThif/V9Y8gYkY3DgD2H017vvvpQqQINOKVU8xwZjp14IbF4Enz/nvtz9qG66KaXapl2XSqniJWPesT8EgXDGe3mWBlOqBGjQKaWKlxFmti/Eyjrbe8/WtS5VadKgU0oVLyPMaqLw0IyN3nvaolMlSoNOKVW8jK7LKAFiBDLe0xadKk0adEqp4mUEXcwEiGeOZ9MWnSpRGnRKqeLZXtDF8We36PQZnSpROr1AKVW8HkPhlDu58skZrI8HCJL03tMWnSpR2qJTShWvrBf2mFO4N3IYLzkHE8vquowVvk6pbqRBp5TqkNrmOE5qicvswSgadKo0adelUqpDtjR5z+IWOUP5Y+IH7DWkL2dOOqwb70qpwkqyRSciF4rIMhGJish0ETminfJHpspFRWSpiJzf0TpF5Cci8qaI1IqIEZEReepYnnov8+cv2/p9ldqR1EUS6eM19OMe+zheCh0Ho77SjXelVGElF3QicirwT+AaYALwHvCiiAwvUH4k8EKq3ATgWuBmETmlg3WWA68AV7Zzi1cBgzJ+/tyxb6jUDmz+c4x57CheCv6ay/z/TZ+uzwg/pUpNyQUd8EvgbmPMHcaY+caYi4F1wAUFyp8PrDXGXJwqfwdwD3B5R+o0xtxojLkWeKed+2swxqzP+Gns5PdUascT2Up54wpGW6sYQE36dH002cZFSnWvkgo6EQkCB+C2rDK9Ahxa4LLJecq/DBwoIoFO1tmWy0Vki4jMEpHfp+rPS0SmishUYHwnPkep0pPMnkfXQlt0qpSV2mCUvoAP2JBzfgNwTIFrBgKv5SnvT9UnnaizkJuAmcAW4GDgL8BI4LwO1qPUjikr6AL0pp47gtcRTiYwdw9Fznm+G29OqfxKLehKmjHm+oyXs0WkHnhERH5tjNmSp/xR4LbsgCO/kJtUanvKWBmlZWrBAdYiAMyG+m65JaXaU1Jdl8BmwAYG5JwfAKwvcM36AuWTqfo6U2exPkz9OWob61Fqx5CxcHOc7LUujc6jUyWqpILOGBMHpgPH5rx1LO5IyXzeL1D+E2NMopN1Fqvl2du6baxHqR1DZovO+InhPaIWW4NOlaZS7Lq8HrhPRD4C3sUdVTkYmAIgIvcCGGPOTpWfAlwkIjcCtwOHAecApxdbZ6regbjP+/ZKndpXRHoCK40xW0VkMjAJeBOoAw4CbgCeMcas7Mq/AKVKVs4zugS+9GtxkuDYYPnyXalUtym5oDPGPCIifYArcOepzQFOMMasSBUZnlN+mYicgBs6FwBrgZ8bYx7vQJ3ght8fM163PFX/IXA3EANOTZUJASuAO4C/bet3VmqHkcx9RidETYCwJLz3g+Xdc29KFVByQQdgjLkVuLXAe0flOTcNmNjZOlPvX0kbk8WNMTNwW3RK7brs7BZdddhPzAQI0xJ0UQ06VXJK6hmdUqrE2d58ubjx06cyRFz3pFMlriRbdEqpEnXE5Vy+eH821jawwBnGkPIAsYbMHQx0TzpVejTolFLF67cXH7KaVU4EgP0rQsRMwF2WAbKmHyhVKrTrUinVIZG4kz7uXRHI7rrUFp0qQdqiU0p1SDRhp497V4S4OnkWZcT46tjhnNprRPfdmFIFaNAppYpmjCGSFXQB3nP2A2BY+QgIV3fTnSlVmAadUqpo5q4TmBWYSQIfP0z+lqrw/un3InG7jSuV6j4adEqpoplYA1XiDkQJ+X2UBbxVUJoTGnSqNOlgFKVU0UzGqEorEKIs6AWdtuhUqdIWnVKqaCZjwrjlD1IW8PF7//180/c+FSuS8NmNsP93uvEOlWpNg04pVbyMFp0/EKI86KOSCAOlBhwg3th996ZUAdp1qZQqXsYSX75AiHDAl7UnnU4YV6VIg04pVTzH67r0B4KUB31Ze9LphHFVijTolFJFk4wWnT/oDkbJbtHp5quq9GjQKaWKJlktujDlAXebnjTdZVyVIA06pVRxjMHneC26QDBIOGjpWpeq5GnQKaWK4yTTh0ljEQwECPosEuI9o7MTGnSq9Oj0AqVUcSw/Txz+HP96fT4BbA72W4gIji+ULmLHo/jaqEKp7qBBp5Qqjgibg4NZYuoBOMKf6hDyByG1KIod1xadKj3adamUKlos4e1FF/Kn2m6+cPqcdl2qUqQtOqVU0eK2F3TBVItubngCJzZcTZwANx3yZfbqrptTqgANOqVUcRwbidZSRpQEfkKpoEuEejHb7AFAQ7Bfd96hUnlp0CmlirNlMb+c8VV+GYbFzmDe8b8AQHnWDgZOoauV6jb6jE4pVZyMVVES+Aml9qLL2pMunmx1mVLdTVt0SqniZARdHD9Bn/t7srsnnSGATTyiuxeo0qNBp5QqTsZedG6Lzg263ZxVLA39AEsMdW+OhANnddcdKpWXdl0qpYqTsWBzwvjT0wsCwTCWGAAsXetSlSANOqVUcbJadL70qEtfyJtHZzm6H50qPRp0Sqni5D6jSwVdMFiePq9Bp0qRBp1Sqji5oy5TQRcIl6XP+zToVAnSoFNKFSej6zJOIP2MLhjygs6vQadKkAadUqo4OS26lq7LcChE0rjHFg7YOpdOlRYNOqVUcTKf0RlvMEp50EdMN19VJUzn0SmlijPxbA57sT8NTU3YWLyZmkdXFvARJ0AFqakFtnZfqtKiQaeUKo7lo97200AF4G3TU6YtOlXitOtSKVW0eDJzPzqvRRczGnSqdGmLTilVFGMMsYyga1nrsjzo4zvxP+JgMbhvD57tOaKb7lCp/LY56ERkEnAT0AuYB8xo+THGrNnW+pVSpSFet449ZA0J/DRaPbAsASAc8LGJXgCEEmGwtKNIlZauaNH9C3gLeAR4DWgArki9Fyh0kVJqB/PeLbwe+hcA15szgVOA7P3omhN2d9yZUm3qiqDbCzjEGJMUkYQx5iwRmQ1s6YK6lVIlwkl6oykdy/sdtixr41UNOlV6uqKPoQ6v5dYsIgHgFuCyzlYoIheKyDIRiYrIdBE5op3yR6bKRUVkqYic39E6ReQnIvKmiNSKiBGREXnq6CUi94lIXernPhHp2dnvqdSOxE5k7Ezg84Iu7PdRTROD2MIgew12pL4b7k6pwroi6N4GjksdzwMOBQQY3pnKRORU4J/ANcAE4D3gRRHJW5+IjAReSJWbAFwL3Cwip3SwznLgFeDKNm7vQWAicHzqZyJwX4e/pFI7oMwWnfEF08eWJVwfup33wxczNXQZiUVvdMftKVVQV3Rd/gjokTq+AXgY2Ah82Mn6fgncbYy5I/X6YhE5HrgA+G2e8ucDa40xF6dezxeRQ4DLgceLrdMYcyOAiByY76ZEZB/ccDvcGPN+6txPgbdFZG9jzIJOfl+ldggmI+iwglnvOVYQ3C3piMeihFGqdGxzi84YEzHGrE8dv4D7hPrfwHc6WpeIBIEDcFtWmV7BbSnmMzlP+ZeBA0Uk0Mk6C31OI25rsMW7QFOhekRkqohMBcZ34HOUKklZQefPHmdmixd8iVjki7olpYqyzUEnIi+nWjsAGGPeM8b8yxhT04nq+gI+YEPO+Q3AwALXDCxQ3p+qrzN1FvqcTcYY03Iidbyxg/UotUPK7LrEl9Oiy3idiOuEcVVauuIZ3VTgXRG5RUT6dEF9Ow1jzFHGmKOAWd19L0pts4w1LKVV0IXSx0lt0akS0xVdl9cC++IO5lggIr8Qkc4++9sM2MCAnPMDgPUFrllfoHwyVV9n6iz0Of1ERFpOpI77d7AepXZIJjPo/NlBlzk4xY5r0KnS0iVLGBhj1htjfgR8FTgZmCciJ3ainjgwHTg2561jyX42lun9AuU/McYkOllnoc+pxH1W12IyUNHBepTaMWUFXSjrLePzhp/YCe26VKWlS4JORCwRGQUMAp4FgsCTIvK5iDwoIr/qQHXXA+eIyHkiso+I/BMYDExJfda9InJvRvkpwBARuTFV/jzgHOAfxdaZqnegiIzHnQAPsK+IjBeR3gDGmPnAS8DtIjJZRCYDtwPP6YhLtStIWiHqTRlRE2jVoiMj+Gx9RqdKTFesdXkb7hSDJmAxsAi4O/VnPbA/7nyzohhjHkk967sCNzjnACcYY1akigzPKb9MRE7AndpwAbAW+Lkx5vEO1AnuNIU/Zrx+PvXnD1PfB+AM4GbcUZ0AzwAXFfvdlNqRvXngFH7zxGcAnFo1LOu9zBaek4yhVCnpinl03wMONMZ8VuD9ZztaoTHmVuDWAu8dlefcNNoJ07bqTL1/JW1PFic1kvSstsootbOK2xlb9ASyO4Mk4HVdGu26VCWm3aATkXOB7wLVuJPA/26MWZtRpG/mkHul1M4plmi9RU8Lyx/CMUKMABk7+ShVEtoMOhH5Me5zqBaTgNNE5DBjzFJIzyVTSu3kYklvwebcFt2CwSfzswX7A8IlI/Zkvy/43pRqS3uDUS4EVuGOLhwG/Bh3GsF12/m+lFIlpt+Wj5lszeUAWUCZld1sC4cCuEvcQkS36lElpr2uyz2Aq4wxLetW3iki1cBfRaTCGNO0fW9PKVUqvrngt5warAXgHr6U9V55IGNPunjyC70vpdrTXouuEncUY6aXcANy7+1yR0qpkmQ53jw6fzB7Hl150PuduVn3pFMlpjOjLjem/qzuyhtRSpU2n/Faar5A9v4ElVac/WUpYeIMqK9B1zFXpaSYoPu6iNTgrjSyKeN8l0w2V0rtAIzBZxLpl4Fg9oTxPrGVPBu6AoCVG3ZHZ+GoUlJM0J0OnAYgIqtwJ1sb3JVDZnZylwKl1I7EsbFSG84ljUUomL1NTyBckT72OzphXJWW9oKuB+5E7Im4e7q17K4tuDt2/1NE1gKzUz+fGmMe3n63q5TqFhnrXCbwE/L7st4OlXlBFzAadKq0tBl0xpgGYFrqBwARqcDtgG8JvgNwF3P+Gm5LT4NOqZ1Nq6DLfnIRLqtMHwc16FSJ6fBglNSUgndTPwCISBkwjg6saamU2oHY3vO5eL6gK/dadEETR6lS0hVrXWKMiQAfpH6UUjubnBZdMCfoysq9Fl0ZMTAGvK0blepWOnJSKdU+2+uOTBhfq2d05aEAMZPxe7PuYKBKSJe06JRSOzmxWGrthknGWW96MzBnrcuQ36KBICHcuXbJWDP+nLl2SnUXDTqldhWJCCx5E4ZMhKqBHbu21wi+H7yRNc0RAN7J6boUEWIEgWYAIpFGqip7d8VdK7XNtOtSqV3FY+fCw6fDQ6eB0/G9dLJ2L8jpugSIiTeJPBZp7tw9KrUdaItOqV2B48CC56HXCOi7FzhJsILtXpYplrHRXO5gFIBNVj9M0hAlSFgXdlYlRINOqV1B/Rr3z5rlEGsEf8dCDrKDLnd6AcBvq//C5+sbAHi+fFinblOp7UGDTqldwdYl3nHv3Tt8ualdxYnmTeKWn3WmDyH/Ca3KlAW97syI7mCgSogGnVK7gro13nGv3Tp8eWLNbP4RuB2AN50JiFzWqkx5MHNPOg06VTp0MIpSu4JYffrw/bU2dtPWDl2eSHjz4pJWIG+ZsoDuSadKk7bolNoFOJG69G+1k7c8gfOPZ+EPm4pevSQZj3p1Sf6g28tZSsCaQ5g4Vk1PoINTGLqAMYZ1dVG2NMYZ0quM3hUdfxapdj4adErtAupqt9Ar47VlEhCtg7KeRV1vJ7wlwAoF3Zfrn+DA4IsAfLBhAHBYZ2+3U6Yt3MRfX/yceeu81uvE4T256Muj+PLoAV/ovajSol2XSu0CmuvzdFXGG4u+PpnIaNFZ+X8/Nn5vJRQnHin+5rrAHW8t5Qf/+Sgr5ABmrKzlR3d/ws8fmqkDZHZhGnRK7QKSzbWtT0brW58rwI57z+hsK5S/kL8sfegkvrige3rWGv73hfnp1yG/xV4DKvFbXrfsM5+u5TtT3mNDfTRfFWonp0Gn1K4glifUYg1FX+7EvZVObF+BNSwDXtCZL6hFt7Y2wu+fnJN+ffCI3rz7my/zyi+O5P3ffoXvHDA0/d7ctfWc/u8P2Khht8vRoFNqF9BkwtSYyqxzJl/4FeBkdl368g/wkIygI/nFBN0/XllAY8xdhWVEn3LuPOdA+la6Lc5+VSH+8d1xXPOt/dOtu6Wbmzj9jg/Y2qR75u1KNOiU2gXcOuBPTIj9m+ftg9PnYo15ujMLMFlBl7/r0gpmBF1i+7ea5q2t58mZ3vzAv5wylqpw64EyZxwynFvOmIgvFXZLNjVx/n3Ts9buVDs3DTqldgE1zW4LpsGUp881N9QUfX120OXvuvSHvLpJbv+g+/dbSzDGPT5mn/5M2r1PwbLH7zeQf542Pj2b4qPlW/ndE3MwLRWonZoGnVK7gPqI273XiNfqinagRVdfMZy37P350BlNQzD//LhARtDJdg66TQ0xnv9sXfr1z7+yZ7vXfGPsYH5z/Oj068dnrGbKtKXb5f5UadF5dErtAppirYMu3lR80C0cdiqXfeCGxLd7DclbJhj2gs6yt2/QzXr+du62HmYq45gx+HTGDi1uPuBPvrQ7SzY18t9PVgPw15c+Z2Tfco7fb9D2vF3VzTTolNrBrNzSzAUPTMdnCbeeOZGhvcrbvsAYJkbeY7Blsa+sSJ9ORoufR5e1c0Egf0dQqKwifbzdgs5xMK9fxbGf3wA+OMw3ly3JebDlP9Bnj3YvFxH+fPL+rNjSzIfL3LmFlz4yi0d7lrP/0B7b555Vt9OuS6V2MPe+v5y5a+uZvbou3TJpk53gH/ZfuT94LV/1Tefw2D85IHobM0ZdVPRntrfpKkCwxwBmOqN4396XhWZ40XUXzbHhyZ8i796QdbpPzSyYcgRMvxuKeOYW9FtMOesAduvj/oIQTTicd+/HrK/TaQc7Kw06pXYwd767LH38ytz17ZbPnAPXYMpYbfqxhR7UJ4vv0Im3s+kqQGi3Q/hW/CpOT1zB3+3Tiq67aGJBmbeQ2WfOCGxSoZtogmcvgecuLSrselUEufMHB1Eddv8ONtTHOPeej9NdvGrnokGn1A4m8//jLaMp2xKJeF2UUbzh9x3ZYWD31U9xoe9pzvW9QB97c94yFSGvpdcYS3b9iEYROP5a3vFP5oHkVzg5fjUfHfMo9N3bKzP9bvjkzqKqG9W/ktvOOiA9x27u2noueXgmtqMjMXc2GnRK7cACvvb/E440ZQSd8ebAdaT1Mmbd4/wq8Aj/E7if3vamvGX8PouygBt2xmyfrXrW1Mf5YeMF/D75I3z+AOMOPhJ+Og3GfMsr9OJvoGZF4UoyHDaqL1efvF/69WvzN/Ln5+d19W2rbqZBp9QOzGe1v81OpNkLujh+etLAILbga1xb9OdYttdytDJXQMlREfK6Qxu3Qzfg1AUbSeAHhEm796E86HeXHjv5Nhg41l1v8+RbO7S57OkHD+f8I72BLHe9u5y7MrqH1Y5PR10qtZOLZgTdHtY6ZoV/CsDyZfsCXymqDr/jDdTwhQoEnZ3kdN/rRH21hEjQED2SAdWdvm1PvAkC5SDCm597rcmj9+7nlQmUwan3u+t3DtwvTyVt+9Vxe7NyaxMvfOY+87zquXn0LA/wrQlD27lS7Qg06JTawfhJ8lXrE6Y7e+GYwq2rFolo9mCUKnHXofR1YAqAz8ls0RVY1FmEy2K3QgAcI8yO/q3o+tv09M9g3WzsfU5iyZJRgDtn7ui9+2eX60ArLpdlCdd/bzzr6j5g5spajIHLH51NRdDPV8d88RvIqq6lXZdK7WAu9z/KDYFbiRFID0xZUxvhz8/N4+lZa1qVT2bsJFCHN9ct4BQfdP6MoPMHCwSd5SMibvBaYog2Fb9odEHxZlj4Mmxdgu/d6wkk3B0XhvUuY0TfinYuBurWwKwHi/qocMDHXeccxN4DqgCwHcNFD85k2sL8zyTVjkODTqkdzMHWfF5xDuQQaz7D7VUA/O6Jz/i/d5ZxycOzWLIpeyJ4ZtA1SFX6OODEKJY/o2zWmpY5opb3XqSxruj6C1r8KiTcFmlN+QgWGrcr8ZCRhde1TKtdBXefAE9dAB/dUdTH9SwPct95BzMiNccubjucd8/HvJCx3Jja8ZRk0InIhSKyTESiIjJdRI5op/yRqXJREVkqIud3tE4RCYnIzSKyWUSaROQZERmaU8bk+Wn1WUptL8YY+kstB1kLuD14Iwnj/iec2er4ZHn2buJ23Gu5RXzeQ7OgKT7oAiazRVe4uzTu81pZHVlirKC5T6UPp/kPA9zBNweP7N3+ta9dCTXL3eMXLodP/lPUR/avCnP/eYcwpKf7PRO24aIHZ/DIxys7cOOqlJRc0InIqcA/gWuACcB7wIsiknepBREZCbyQKjcBuBa4WURO6WCdNwKnAKcDRwDVwHMikrsMxI+BQRk/92zL91WqI2JJh37UMlDcnQfWJKtJ2E5WmYZo9mjHiJQxxxnBImcIm0LD0ueLDjrHxo9bp2OEQLDADuNA3O/teZdo3sYWXSLidlum3FM3IX08qZgW3TdugCEHeq+f+wVML+4/16G9ynn0/Mns3s8NbsfArx//jJtfX6Q7HuyASi7ogF8Cdxtj7jDGzDfGXAysAy4oUP58YK0x5uJU+Ttww+fyYusUkR7AucD/M8a8aoyZAXwfGAsck/N5tcaY9Rk/X8wOk0oBsViUkLihkzQWNckgtc2JrDK5k8iX9z6Mb8Sv4dj433lq0C/S50Mk3GW12pP0AjFGgLJg4TFsdsALuuS2Bt3KD9wVT4BYj92ZGXMXXh5YHWZY7/YH4RCuhrMeh8FeQPLsJTDz/qI+fnDPMh796WT2G+K1gq97dSE/e3CGrqCygympoBORIHAA8ErOW68Ahxa4bHKe8i8DB4pIoMg6DwACmWWMMauA+Xk+95+p7s2PReR8ESn4dygiU0VkKjC+UBmlOiKescqJXxyOsD9sNV8tlshu4UUTXpj1qQrRnDFpnEQRv6dlbLkTJUhZMP9alwBOZtBFtnEwyprp6cOllRNp6bY8ZPfeiLQ/fxCAsp7w/Sdh0LjUCQNPXwTv3ACx9he17lMZ4qEfT2Jyxl53L3y2npP+9S5z13bBM0j1hSipoAP6ApfuVrgAACAASURBVD5gQ875DUChMb4DC5T3p+orps6BgA3krm2U+7l/AE7FbeU9DFwH/K7gt1Gqi8UjTVmvrwtMobY5zo98LzIlcAPvhi4mEMn+Vz0S94KvZ1mQCMH062Qsu768LB/P+o/lcftwXrQPcidpFxLygs6JbGMQrJmRPvwoPjJ9fNCIIp7PZSrrBd9/Cgbunzph3Od3N4yB16+CSNsb0FaFA9x77sGcc+iI9LnFGxs5+V/vMmXaEpI5Xceq9Og8ug4wxlyd8XJW6vnd74E/Fyh/FLgtO+DI7X1/aueXG0zlRKlpivF13wccYC0CoLppeVaZaMbOA+GARTwj6Jqbm6hub1J3uAdXyflsSrhdmB+20aLzlXmVObGGdipugzGw5pP0y5dqvD3wJg7vle+KtpX3hrOfgXtPhPWfueeitfDBFDj04nYvD/gsrjxxDGMGV/OHp+cSSdgkbMNfXvycp2et5aqTxnQ8gNUXptSCbjNuy2pAzvkBQKFl2tcXKJ9M1SdF1Lket9XXF9iUU+btNu73Q6BaRAYYY3JbjEp1uUQ0O+h8YqipbwDjjXaURHNWmT41n3Kq72NiJsDAeIiYuF2XzSZENNJEMYuXRDLWrWyr69IX9vZ0k9g2dF3Wr4VG9z8pJ1DOhw193c8O+NhrQGVbVxZW3hvOfRVmPQDv3QI1y2Di2Vk7IgAw+79gxyFUBaHq1E8VhKr47tjeHLjb4Vz630/5dJU7qnT+unq+O+V9Thw3mEuOHsEe0Xnuai7xhtSfTW43abzlp8n788zHwCr896m6RkkFnTEmLiLTgWOBRzPeOhZ4vMBl7wPfyjl3LPCJMSYBUESd04FE6tyDqWuGAvvgjtAsZDwQBbpgHLVS7bPzPFdqrK8jhPfcLTfo9t7yOucFHgFgek0ZP624iQVbYoDwWsVIctYXacUYQ3Pcew5YHij8P2YZtB/PfXoIjaaMeXbnVyqhfi1UD4H6NdT2GIPT4D5l2X9oD/xFLGRdUKAMDjoPDvghfP48DM7z+PyNq6G28FSCkWLxVLCSZHmCExLXsijh/g0+8+laps5exOzQj4u/n0SzG6JquyqpoEu5HrhPRD4C3sUdVTkYmAIgIvcCGGPOTpWfAlwkIjcCtwOHAefgThMoqk5jTJ2I3An8TUQ2AltS18wGXkt97jdxn9e9D0SAo4GrgH8b04EJSUptAzvW3OrcmMW3caDvw/RryRlgIrb3r6cVCBMMlwHuyMxiRg/Gkg4tO9cEfVabQePb79tc9KzbquuXDHFVu7UXMOwg+OU8aFjPY6/MgNXuDUwY1rOzNWazfLDvifnfa6/L1ThIrJ4AcP/Z+/OnjyS9RmaTKbBqTCHxJg26L0DJBZ0x5hER6QNcgTtPbQ5wgjGmZd+N4Tnll4nICcANuNMF1gI/N8Y83oE6AS7F7e58BCgDXgfONsa09NkkgAtxA9ACluIOTvlXV313pbLEGt1RkZXe4sVOnqDzR7ZkvfbZuUGXuU5lKKvrsZitdOKrZvI3/+00UM4S/yjgawXL9iz39rura05gjCl+hGQ+VQN5fXMvwJ0EP76rgq4t486AyFY38GL1qT8zfjJazANCSW49cxLTV9Rw8xuLmLpgE+/Z+xInQCNhmk2YuK+Mvn36MLR/P0YO7k9FVU8IVkCwEjK6etX2U3JBB2CMuRW4tcB7R+U5Nw2Y2Nk6U+/HgItTP/nefwl4qa3PUGpbGGOIJGx3VOPSqfDAd903vv8kjDgcgPX9j+DM6L95JfQrBojbYx6OZ6+E4ktmB52V0aLzBcsozwq69lt0yU2L+J5/GgCvWYVm+bjCAR/hgEU04RC3HZrjdtbWPR1lO4bP1nijN8cP/wKC7vhr2n7fTriBZ/ndsAIO2K0Xd//wYGavrmXKtFt4dd4GEnaqGZzE/fV7LcinMG5oTybv0YdxQ3swrhkGVm/jLwOqXSUZdErtauoiCb5z23us2NrMTaeN5/jXL3UHRADOjPuxUkEXcyzqqGSxM4QBPjfoyhPZj4j9dttBN9iqYT9ZShlx7PohtB6nlS3Z7NUf87W/kHLPsiDrE+7cu9pIYpuCbuGGhnSrc2B1mEE9ipgovr35Au7AljzGDu3JrWcewJbGGE/OXMMjH69i0UbvuaoxMGtVLbNWeX+n/apCjBvagz0HVLFb73KG9ylneO9yBvUoK2q/QdU+DTqlSsCzr02lctNcRhDmkbf9HF/jbfwZmf9Kes+BWGqqQDPes6BKJyfocnYl8GUuyBwMc1LN3VwTegGAj1dbcEjb6xlkTvyOtxd0jZv4hTxA3F9HA+XUNh+eXjOyaCvehw1zYMhEZq/0uva+kG7LLtKnMsR5R+zOuYePZMGGBt74fCNvfr6R6Stq0s87W2xqiPHa/I28Nn9j1vmgz2JorzIG9QzTqzzo/lQE6V0eoFeF+7p3RZAeZYF0Szoc8BW16/yuRoNOqW2Q+PwlIs/8isYhhzH49Fs6PVT8uE8v5qxQaoX8nIkqDbYvHXTxpDs5uTljlGVvkz0xO5DTovPlbLFj/F7w2PHWz/xyOZlB529n4ES8kVNjj4MfVpu+LG9KtF0+nzmPw8fubgO9+p9Ly+awX0i3ZRcTEUYPrGb0wGouPGoUtc1x3luyhU9Trbo5a+poKvCcNG47LN3cxNLNRUzqz+CzhJDfDb1w6s9gy+uMMLQELBEsESR17LO845b3JePYElKvU+eKLO+zJGsX9y+aBp1S26D+ycvpE1tF9aIVLJ9xEiMOPL5T9WTuDpCrKZkRns2bGSqbaDYhljiDqKGK/WVpev1LaL3PXOZecoFQmTvEPsUUEXQmYz5c0t/OHLaQNyuvkgibGzsxIDlj6a+pjd7Ys3FDd7ygy9WzPMgJ+w/ihP3ddTttx7B0UyOfralj+ZZmVm5pYuXWZlZubWZzY+F/J9piO4bmuF3UQKMvigadUjuwPrFV6eN105/tdNCJKbyMlJOxUt8+S+/indC9APwlcRpT7BO5J/AXjvTNBuBZexIPcQyZA+czQzQQKkcygy5RxOarUS/o7GB7Qee9X0mEjfUdXPPcTsCGuemXL9W4zw8tgbFDd74Rij5L2HNAFXsOaN1SbowlWbmlmc2NMWqa42xtilPTFKemOcHWZvd4a1Oc+kiCaNIhlrCJJh3s3L7REtDdjxo16JTqLDu7Wy6+enanqxIKB91d9vG0jAOUjBGVkVT3ZVi8IHvAPoYPzR5Zw/r9xrvPQKgMCWZsnJpov0Vnxb3BFE6wna5Lf4ikBPGbOH5x2FrbwdVRNi2A1OCZWMVgtkbdFuJeA6q2aVDLjqgy5GffwcWsW5MtYTtEEzaxpPtnNNHy2iaWcIgmbeJJAxhsBxxjcIzBmJZjUq+9Y8e4o4IdJ/P9Isob9zM06JTaUeVMLH7OmcSwTY3s3q/jS1TlBt0RZU+wuaaWCmJsppqrHeM+P8nZSQAgjBd0ERPEGPf5Tsjvdnl+bEazyulFiAT7V/TMCrrcyeX5WHHve5oiJjcnA5X4U1MeGuq3tlM6x3rvl4W14b3SxzvSQJTuFvBZBHwWOg3do0GnVGdltIbWm148ah/Fl9c3dDLostVGkkQIE0mNrowkbCpDfnyZQWfcydlleM/BWlp5saQXdFckfpjuzlrYZzd8mUGX7FjQ+cPttzBMqApSQRepy90QpB3rPk0ffuaMSB+P06BT20CDTqnOSmSGjtu6iiU7t2WLZbIHDuTuEt4cS1IZ8mPZ3md+0/c+PaWJO+0T+CCxD2HirDDuM62v3/Q2PhH+edqEdMhZAgGf4At5UwQy6yvEn/C6LoMV7QeOqRwIDe6iQ77GQmuxF7DOa9G91TAofawtOrUtNOhUSdvmJaS2IyfenB4m0rLHW6KTe5O1HoxiyGzntYygy5z8fYxvJsf4ZvKifRBv2OP5Z/BfhIlTRwWnbv0DAP98fVG6fDjgQ0Twh7zBKL4igm5q7++yetUyKonQq7rtyeUA/l7DYJ279maoeV3x/wwdx9tCB3i70d2apyzgY8/+ndyxQCk06FQJe3rWGv707DyOGzOAa789trtvp5VE3JvNFksHXSdGvBlDkOyh5JOs+Vzpv4c4fupMBc3xIwDwOa2DqYIoIIy3lgCwyXjdi2987k1CDvndWA6EvdDwFxF0L5edwCtJd3LfbT36tls+0Gto+riPvYmtTXH6VIbauCKlZpm7tQ0QC/VmQ9TdPmebdyxQuzwNOlWynv/vHVxqzeGuj49n3qQRnRqB1pUWb2zg8kdns1ufcq777jhiVSO4LH4xp/veYKhs4kr/3VTV/Iicdcfbl4zhyxmM8nAwey/f6bEoUI3fbj0vrUxiWRPIyzOe2QkOv/Q/RswECEo58FUC5T3YZKqJmhBbaH/Ifl3EG7XZoyzQRsnUZ/bwgm6gbGXl1ubigi7j+dya8J5Q57YCu2zHArXL0qBTpalpC7cFbsQnhnHWEjY2nMi+GVuEGmPabCmsqY3w4mfrOHp0f/boxOCQfH5498es2hph1qpaTpk4lH0G9eI5ZzLH+z7mMGsu51iv8EbD0R2v2I6xnEGMYF3BItFmd3WM3MngACNkPX4yNkYlTkvXZ4gEF/ufAiBmB4GbsIYdyEGxKQAMrSjjy+3cXmbQVRcRdAyZyOs9vs0Hm8PMdvbgjK3NTChmV/B+o+FL/w/Wzeb99V5Y6kAUta006FRpWjYNn7jdgOOspbyb0XVljOGsOz/k3cVbuOzYvbj4K3u2uvzC+6fz6eo67np3OW/96uguWRx31VZvhOLr8zewR+q5UQJv5RLH7sSSV+EeHJu4nqHOWmwsjrU+4X8CD2QViUXcoPM7rVt0/aSeT8M/Sb+2xLCfLGOwbOFjZ+/0+bgECIG7O0JKpIjVM+o72KJjyAF8uPevuGPDUgAmF7uE1YB9YcC+2I7hmitfhlR460AUta006FRpyhmc0fJ86a2Fm3h57nreXezuwXbdqwvzBt2nq931H9fURthQH2VwsQsLb14ET//M3d362/92V6rPw2dZ6XUnEybjP6Nkx5dtiicdErZhGe4owyZa32si6j67srFIGgu/tD3o5bnQFQA8bx/s3VrqOWJ5R/ajm/1f/hP9M0sCA3jRPoTqsq+2/4WAURmt6AXr29nINMfijY3p9R/7V4UY1KODm5kqlUODTpUkY5ysuWUBn8XctXWc/Z+P2r3WyVkCqbGIXbTTnv4ZrErt1j18MhzitZTGyWIOtBYwVDbjbziReHIYAMmMFp3pRIsut1XlDi7Jloi6raKf9fk/ZqyspScNzAr/tN26v+7z/r4SVmreXcC730jCxnEMVoEWb2LdZ4yWFYz2rWApQ6kqcnWSfQZ53czz13VsdZRPM7awGTesZ8mOulU7Dg06VZJsO5n1L6djDHc+/gL3B25hnhnBs/ZkhslGXnUObHVtPGeIf+Yzpna1hBzA/Geygu4U39uc7X8VgJei+1Ex90GmBq9jhOVtN9CZoGvK2fy0PE/QJaPu5PSWeXr1tL8vXKs6xH2eaWE4LjAbsaOEiBNJHNdqea36aIIXP1vH4cvmMSR1rrZseMFAzLXngEp8lmA7hhVbm6mPJqgOt9Ht6Thgua32mRlBp92Wqito0KmSlOi3f6ugu6T2WnbzreRw5vIT//MA/G/iDBL2N7P24IonbSbKQsISJ0iCuubWYViMSGNtuhMxvnpmOuQAfMlmpHlTVsgBjKj7gK0vXE2vw89DqgdRjPiWFZxovUccPzWmipFW60nWyVSLrqW71MEiagKEpe1gzSxjW8H0+dt9f6GlIbox9j9ZQWc7hjPu+IClazbyRmh6ejpftHpkUd8HIPzxbfy34nGqYhv4Y/IcPlx6IMfu28YcvBd/BSvehRGHs3XZWMDd2FRHXKquoEGnSlK8eng6ZCImSMI27GavbFXu94EHqU/clBN0Dk+Erky/fjb+g07dg7XFm2wd3bScYMZ7Em/Cibfee26vuvfgo/eYPecdxv7q5aI+R1Z/zE3BW9os4ySyW3QATYQJ03bQPWAfw7n+F4GMoLMsYgQJpebuNTY00L/aWxbsncWbmbOmnqv8DzFQagCoMZUk++5b1PcBYP1nHJCYARYMk428tXBT20G37C3YvAA2zqM+/nugNz5LGKtBp7qAzsJUJSlBgJ/Ff85P45dyYeIS4kmHl4PH5i0bS2R3VTpbl2e9XrmhpujPbTDeQJCQ8boQE/WbsspJshmnjQWRxzZ/kN4NvD12c/v358SaaW5uYp/op4yXxYyS1bzljE2vd1lIJCOebcubihEVb4BHc1P2M7R3Fm1ijCzPasH+KXE2Q/oVMUWgRc9h6cMxspznZq8t/PfRuNENOcCxAsx0RgHutjyVu9iOBWr70KBTJSlhhOedSbzsHMybzgTqIgkSJv/u3S3/AzXGHYTiWz416/1FazbmXlJQQ54RjwCJpuww8iWaId72gshNseKCzle7rN0yHy1czXF/fozbnT/yVOgP3Bv8C792LuaSxEXpMptMD86J/yrruojxws3xecdxy/uekcbsoPtkRQ3H+j5Jv37DHs9TzmEcMrJ3Ud8HgBFHpA+/5vuYuuYY93/QukUOwPJ30ocry/Ylmpr8fugefYr/PKXaoL8uqZKUzFlKa3NjjJ7kD7r1dVHOv386kbjNf845iPKcAIpGipzHBZzt+xuDo4toNiHmmRHMTNqE/L5WrS6/3YxJtj3svSkSo7lmPW8+PoXVvSbx81NPyLunWvnW+QXr+EPiBzSaMj4xexPM3KXAhOhfHSJc55171xnDVGc8q01fhoq7a0AUr8VnfF7rLm6VtUxTI9KUsQ2PMSxc38DlsiB97in7cA4b1ZeDRnQk6A6Hin7QtIn+UsvB1uf8/WU/k3fv03qFm4ygeycxOn186B7tLzemVDG0RadKUmjOQ7wY/A1PBa/gbN/LbGmMYxcIuiufncucNfUs2dTEufd8gpOza/bIHvmvy6dX/yG85YzjEzOaZsLpEZtONLvV47ebs/Zy22R68ELGnDWA5sY6fPd+g+/X/IuzFv+Sl2bnadEkIvTZ3HrKxAxnFGOj/+Ze+ziecL7ESjMg63lcjCD9q0KUZWy62rKDQjXe9kEtE88fTX6JNb0y5tT5vWdysWbvu62pjRCLx5hoec8nLzrnbO7+4cFFj7gEwPLBvielX55svUM04XD+/dOpj+Y8V1z+dvrw+Ua32zLoszhgtw50lSrVBg06VZKkcQP7WCsZby1loNQQSdh8J/FM3rJz1nj/o16ysZ7gxk+z3k8Wsbloi5wpeNRH3KH/JpbdKgzYEcjYG+5/Ej/kwsSlbDTe4IlIUx0bon6m2uOYYfZk1ZI8Lbea5a226AGYZo9j8MDsUZuhzA1WCdK3MkQoI/zcvegM7zvZg0buSh7H/0uez7zhZ3jf0+91XSYi3jY8izY0soes9QK0eih77bV31mCfoo35dvrwNP9U/uS/i/Vb6/jj03O9Mg0bYPNCAGwJMMNxJ/8fsntvwoHif0FRqi3adalKkkl6XXIX+p/hmugFecvFTIBhsoE/++9itelHlTTTa/n72XW18ywtU8J2EBwu9D1DlTTTULcf9K/ExLODLuREsvZya+kibDIhbzh+cyPGhDjK5wbvbY0n0crW/M/nyiXK2KE9WLZ+C9U0s5us57HQVd7nmSCHJj+mXuq4O/lVXnUOYGFgH0gKP038klPst7gu6K5nWS3uvbesLgPgBLx5eMmI94vCwg0NGIQX7IM5qHw9/Qbt3+7fWUHDJ8PQg2G122L9gf9VxllLuXDmJTy9dz9OGj8EPn/O++zAaGIRt1V63JiBnf9cpXJo0KmS5CSzu7eSydZrPAJ8ZkZyS+BmxllLC9cVby74Xtr6z+DzFxgR7cNX/O9wif9JAF54sTcT9rwZSWQHXdhEsTLuqWUARTPec7vG+jq2mn5MTr3ul2w9Py6yZVXe4S/VNLP/kB5snDGXu4N/b/X+VqqYEH2fcf6nAViYGMaksq2cxkOUE2ON6ZtVF5DVQjIBr+uyZY4euMt1LTTDuDBxKVcfNYbvT9otz90VybLgrMfd1Wbmu63x8dYS/hS4h18+OZADBgYYOvUv6eJPNnuh2uZUBKU6SINOlaTcFUaSifxrSF4Qv4SPwz9ru65EO3uubVkCUw4H4CbI+q/ihK33AjcjCS8sFzlDeMo+lLNsb5PQludjT9iH85YzliYTpqIpTMT0S5cpj7QOuppNq9NBt5Z+DMadxnC6/02WOZ/w/TwhBxD0+wiVe4M6yokyqKon/VLPEm3jtd6qpXXQSchbi9KONcL6OTDvaRrW7gFUAbDXgCrY1uW3wtXwvXvhg1sxr/6BGezDL6IX0EiSD++7gqFN7ojYhmB/7oseA7ijLQdU6/qWquto0KmSlBt0psBiyYli/hVORGiMJbl92hJ6lAU49/CR7vqJ0ToQH/z37IKXfuTszcGAP+kF3S8SFzDH7I5/0Al8OG8xk615XOp/PN1FucAZylnx33FIfYieqdAACMazR26uq4sQ2eptzbOuYl8GN00DoM7qwaDywpu4xgcdSEWV92yvUiL06DUKNrUue4j1OX/y38XQLd8FTnVP9hrJrOW702zC1MQE7v46RGv5s+nF21xPlBB7D6xqXVlniMDknyG7HYoVHUjk/2aCY/jP1rEc1HcRwxtm8Jf494ikWsNnT96GVqRSeWjQqZKUG3S+eP4V8JMFRmJmqk0GufOtpbzx5qtsND3pVxXipJ7L4L5vgT8MscKLDh9sLSA++3Ge6H8RsxctpYIoK01/AD5LDuUDJ0hPGrkk+ET6mr2t1dhYDJl/F78NPJQ+H0q4OyqwZjpz3niIS+bvzUTpxaHWYfSjliUDT2Z8xWCY9zQVx/4BX2X2qMPZvY+nb3wVawIjOOx7l1E++570e/vISnYf8hVYmHptrcq69gf+V1lcv1/6dfNBF3PyhxMA+J/IyxB115ccIDV81zeNVyq+Sc/yIF1q8AQmAJd+pZHrXl3IXDOSL226jK+G5vNqzJ1WMLx3OV/ZR7stVdfSoFMlKTfowon8q4e016K7OH4R7/nGEp46hedDdxI1Ab708I2ctO+jYMfdn3YEn/gRS52LedmZnHV+8UZ3tGK+sD3Cmk1vyQ7QMrse4k0k7j2F/WI13OTfja/Hr+VR+ygALhu8F76v/Ai+dZvbClr1cdb1W4Z9lbHfOpfBLSeqvedwx/k+wYlMavN7ZHZX9q3yQmx0bHZWuasDdzOmPAArymG3Q9usszN+dvQopq+sYeqCTYDwSswbJfqHb+zbuRGeSrVB/41SpSkn6MoStXmLfc83tWAVC50hbKWKSMLm2sCdAIQlwe8CD8DSwtflM9xZ3ercyq1ud2a+sP114GF+mlp4ukWV00Bs1UwCMTe0x1grsnYq6FuVWrmk5blYzqLQo/bKWWuyLHsCt9V/b9pihb2uyD4V3iopN8ZOpPnAC7PKnlZ7B3w4pc36OsuyhH+dMZGj9/aeX4rA704YzTE6CEVtB9qiUyXJODlBl6zLW+6P/nt5xx7DAKllT2tN+vxxsb+wwAwHQBI2eP9f52TfezSZEBWSfyRnPo2pNTAv9T/GAGqolAh/TJzDVqqJ5/nPaI4zksG+rVnnqmlk3YoFjMg4d5n/Ua5Ofh+AfpUhsi8YgtN7T6yti4iV9WfY6JxdGCpyVg7pN5q2SEbQBf0Wg3uEWVsX5WNnb/Z9Z29eDj7D3lZGoA/chqkF7agI+fnPOQfx8fIalm9p4sDderF7xmatSnUlDTrVJV6eu56rn5vHl0f356qT9mv/gvbYOXu0FQg6vzick/g1SfzMDP2EXuJ2JzpYnGy9Q5U0U2NaD6rIDDnHCJYUHvgB3q7fJ1nvMjK1Nc83fR/QaMJUSutRnUdYn7U614sGPl+zKOvcuf4X00HXaqShCNZp98Ocxwjt883Wu50P3B+qBkPDWgj1gIFjoc+esCX7M1r4yzKW3mrazI/K32Zz4zoaKOcB+xjus4/lz9ZdgLtcmIw/s/BfSBcQEQ4e2ZuDO7KGplKdoF2Xqkv89L7prK6JcO/7K1iyqbH9C1KaYkkuuH8659z1ERsbvMCYt/sP2WR6pF9XOvmDDuBIyx3tGM5YOWSIbOLG4K1cHbibs/2vtHkP7YUcwDWBO1kePiMdcun7yhNyQNbSXADT7LE8a0/G3ry4Vdnl4TP4rm8qA3qEWr1H/9Hw5Stg0LjW7/lD8P0n4OCfwGn3QyAMpz8Me32N5km/5Kz4b7OLl3l/n9Sv5byaG/hN4GG+73N3KXjYPpoP/AcBIMf8CaoHo9TOQFt0qsvNX1fPHu11Qzk2ROu49e2NvDjHnV921bPzuOWMiVC/lsOmX0qleOHW4IT5RPbjQDOnVVV3Bq/jhsQpWeFiMn6HG4A3kCVmAsTwUy3Fr5bSnmYToryNbtBn7UlcnPg5YWKcEZvLMjvGd3xvZZU5zTeVvhV5gq49/feBEzLm2vUdBWc8DPEk70x9mcXOYEZZawEIlGW0bMu9VlQvaWBwjzDP/fwIepefCMakd/tWameg/zarbWaMIUSc3tQTJMFvH3y31WCSZZubuG3qEtbWRiAZg1snw99H0fjR/XzXN5VpwUsZNvd2Fm9swPn30VRG3P8515hKvh77Xx7ieO6Tbxa8h18EHs96fU/wr+nj3Sxvm56tVHFT8ttZZW/MeQ0U3BIonwar7flmS8xg+lLHtNAvuCJyHe/b+3JpPHvwR52vV8cWTW5H2O/ef0VGizNQnrGJacZAlgFSy5uDb6V30xJ3VIiGnNrJaItObbPE9AdYkLs6ydXnMaPqaEb97FEAfn7df9jHWsExL01m3ik16Y02/2TfTMtOMr+2Huag67/Ex2FvBZEFZhhTAjfSM9nIHcmvQ9v7jLbpNXsCx/hm8mVrZtb5S/1PtCo7w+zJIfJ5UfU2iZroAAAAEwlJREFUWdVgby74/kbTix/7n2OAuCNHfxF4jDPjv8sqs0W6dqV+yxJ6hqA33vzDYGVG12Ww3J1DmFqYOrTsNbD+mluNUjsFDTq1beY9Q/C5/EtwLatN8sEHK9iwegVPBv+AXxyOsD5j8eJxjCpQ3Zn+17JeT7K8Ff9/GXhsm271fvsY+ksth/rmAfBJ9TEcWP9aq3L3Jo/l2uTpHGPN4ObgLVnvLej7VfbenP3Mb0tgELvbhdfaDPUazD513lY8Q2Uz6032AAxfZb/cy7bZaN8aEraPEAlqTCU9K3pmF/CHsnZgoPfILr8HpUqB9lGo1gosoJzXYz8q+JZguOm1BYxddS9+cQB3pGJy/gsFr8nXuuoqp/qmMtbydgswPYblLddDmogQ5lnnUP4ncU76/ILQflDeuuXVGOrHuiHH563LMcK4fhZf8mWPwgyk/j5a7LXHHsV+jaIZx05t3QNvOOPdZc8yleV8F0u3xVE7Jw065XFs+M/XMH8dydr3HsKY1qMRjTEs29zEog0NcGUPyJnvlunbvncYZK/jlNhTWedH5yxP9UX5mi97pZHy/vlbMCf53mPs0B5M2r03swd8my0HXMrW4cex+9m3IeEerco3mDIG/fgR/jv5aa4deAMfBQ9Jv3d+4lL2KG+9e8I5Q9aQPOgnAMSCPRlz9Knb8tXymuOM4O7kcTyUPJprEnmmChxxuXc84awu/3ylSkXJdl2KyIXA/wMGAXOBS40xb7dR/kjgemAMsBb4mzFmSk6ZNusUkRDwD+B0oAx4HbjQGLM6o8xw4F/Al4EI8CBwuTGm/bWkulEiUs87yyOMGlDFsN7lkIxhXzsMnx3jW7E/8bdLz2XImucpX/keAgx+5XyueX81pw1cy/Bv/hb7gVMJbZiBAB3p4HozdNl2+kbwuj2Br/hmtl+wgOH7/P/27jxMqurM4/j3193QgAiySDfIpoIiaEDaDZXVgIREJaIxGlHJjJOIEuOWROPETJ4xJkFRnkTHDO4YMxmTTNxFJgiK4oYyCsG4BBQ3FBFEsaGXd/64t5rq6lvd1d21UfV+nqeerrrn1Dnn5Rb11r333HsPh5XRZRP77uSiGZPCV+MblpdEbNFV1wQXV/7G8ROACay5cXHDxZW763NKBn0l+LTF+fKYwykbPQeGTaO8574Zmcpf0a0zN26aDkCXjhFba6POgHdfgC0bGic95wpMXiY6SacB84HZwPLw7yOShpvZ2xH19wUeBm4DzgSOBW6S9JGZ/akVbd4AnESQ6D4mSJwPSqoyszpJpcBDYdlYoBdwJ8GtNuek/1+iHcxA4ok16xl370g6ABOByTt+RR99wpiSv3FBWbCL8n/Kr2Lw9UP5jw638JW478Mrtv07bAPm3ZZXH5Rf1nyT960nL9mQZhPd72sn8rr15ycdFkaW79l/ROTy7VbO6eNHRZZ16Nr42Nq6+gq67nd4o2U1nXcdb+vDFroPOZK1fb/O/u/dz5u2D59PmcthVeF1M/efmHT87XX6EQO5+uHgGOe3j4n4eVJSCifMz1j/zuWLfPr+incxcIeZLQhfz5E0FTgPuDyi/neB98wslmzWSjoSuBSIzTtvtk1J3YF/AmaZ2WIASTOBt4AvA4uAKQRbjIPMbENY5wfALZJ+bGbJL4PfRv/6l9UsfOYtAG7o8Bumlz7dqDx2ZY619QOYV3sqCzrOa1Q+LqG9xeU/iOxnfacz0jbmTHu6fjg1laO58sDN8Ex0nTPsat7Z4wCqti2JLP+iyz507tSdtzsNY2B1MLuyep+j6fTu03QcMp7K3tFX69hjwEjuqZ3IZ3Tha6UrWFI/mimTGv/brel3Kj9+/UA2Wg8+Zk++370z/c69neWvfUDfHl05rCJNt79pwTnHDGbz9p3U1NYz57hk03+cK3x5l+gkdQSqCHYhxnsMSHYp9TFhebxFwNmSOhBscbXUZhXB5PWGdsxsg6S1YZ1FYT9rY0kurp/y8P2PJ8SyNHwavXnQgkdXf9CQ5GaVPtIkycGuK3McVLKhSZLbHVxZM4vLy+6JvO7kgtppnFvWdOJKn74DWPC9sfDqQ2x5sZJt1XUMKNl1I7atHStY+KPzKS0Rbyz/FJpOrKTzpMsA6D3r96xZeid7HjKNgQeOho/foKz3AUnH23PoUfyi7Dw+ra7l57Xf4p5/PpIBvRqfHN+1cj/WhL95hlXu2XA1/nHD+jZpL5M6lJbww6nNX//SuWKQj5NRegOlwMaE5RuByiTvqUxSvyxsL5U2K4E6IPGEqMQ6iW1sCt+XbGxt9sDL7zU8vyrJ7rd8cV9dy7dzOaT6FnYeveuY3QU757B8r5Pg5AUNy/5SfhLr9p/J5jMXM/D063lm5pvUz1nF450nA/BS/RBOmRTu9us1hO58zoCSj/j0iIsa2rCZ91Eannw95LApjcbw8awVMOtRGH02AF0q9mPEaf/GwOFHBteS7HNQs7MPS0rEVSeMoG/3Tpw4sh9H7derSZ2pB1dy+OAedOlYyoXHDW3x38U5l1l5t0VXSMxsAjRs2Y1vtnKEE77Uj4deDu5APWvnZdzecW6z9XdYB8rVdBbk4roqNlk3BnSp4didyxuWX1NzeqMbg8bMtbM4ueMKrv98Kpd1uo9B9cEhzFfrB7DGBtGnV296jZlJRf/BrH1/G4cOH8ZxpR34ZOXd1C67jh9tOxVTCdedsC/UbKfHXy/l7lG/49lpU+lYP4V6q2EzXZkx6BzmH9CHUsHOupv5cPMWTphwLqVlwcfy+IYR9WbcJf/Fg0uWUL3XEE4eHv6m2PtANGcl7NhGt95DYNpPAWh0tlinblR/51nWLb4ZHXUewwY1fyubVMyo6s+Mqv5Jy8vLSvnv74yhtt783mrO5QFFTSHPpXDX5XbgdDO7N275jcDBZtYkYUh6AnjFzM6PW3YqwYzILgS7LpttU9IkglmWfczso7g6a4A/mtlVkn4GzDCzEXHlewMfApPMrNGuy7g6S8ePHz9+6dKlbfgXCab01xsNWynOOVcMJkyYwLJly5bFNhraKu9+bobT9FcCkxOKJgNND1IFViSp/4KZ1aTY5kqgJr6OpP7AQXF1VgAHhcvj29hB0onq7SfJk5xzzrVRvu66nAcslPQc8BTBrMp+wM0Aku4CMLOzwvo3AxdIugH4LXAMcA7BaQIptWlmWyXdCvxK0ofsOr3gZXZNZ3iM4IyouyRdQnB6wVxgQSZmXDrnnGu/vEx0ZvYHSb2AKwlO7l4NTDOzt8IqAxPqr5M0Dbie4HSB94Dvxc6hS7FNgO8DtcAf2HXC+FlmVhe2USfpq8BNBMnyC+B3BCehO+ecy0N5megAzOwmgoQSVTYhYtkyYHRb2wzLdxCc+J305O/w5PKvNdePc865/JF3x+icc865dPJE55xzrqDl3ekFhUjSO927d99n1Kg2XSDFOeeK0qpVq9i6deu7Zpb8xNUUeKLLAkkvAXsDb0QU9yKY4dmS5upFlaWyLP511PNYZl6Vwviak6kYo5a35nUvIHZTumzF2FJdX5epvc73GJOVt7QsW5/XZGNpbb1Mr8vPgI/M7NAUxpmcmfkjhw/gP9tbL6oslWXxr6OeA0uBpfkaY0sxpRJztmP0dZmedZnvMbZ1XWbr81os6zL28GN0ufdAGupFlaWy7IEUnqdDpmKMWt6a1+mMszVt+bpMbfnuHGOy8paWZevz2pr2dud1CfiuS9eM2N0XrJ2X38lnxRAjFEecHmPhSHecnuicc84VNN916ZxzrqB5onPOOVfQPNE555wraJ7onHPOFTRPdK7VJA2QtFTS3yS9HN7ktqBI2kvSC5JWSVot6dxcjymTJHWR9Jaka3M9lkyQtD78rK6SFHmD5EIgaV9Jj4f/N1+RtEeux5ROkg4M12Hs8YWk6S2+z2ddutaS1BeoMLNVkioJbjp7gJl9nuOhpY2kUqDczLaHXxargcPMLNWrn+xWJF0NDAE2mNmluR5PuklaDxxsZp/leiyZJGkZcKWZPSmpJ/CpmdXmelyZIKkrsB4Y1NJ3j2/RuVYzs/fNbFX4/ANgE9Azt6NKLzOrM7Pt4ctyQOGj4EgaCgwDHsn1WFzbSRoB1JjZkwBmtrlQk1zoROCvqfzA9kRXhCSNk3S/pHclmaRzIurMlrROUrWklZLGJmmrCig1sw2ZHndrpCPGcPfl/wHvAHPNbFOWhp+yNK3La4HLszLgNkhTjAYsk/S8pG9lZeCtlIY4hwKfSXpA0ouSrsja4FOUzu8e4BsEN8lukSe64tSVYFfchQR3SW9E0mnAfODnwKHA08AjkgYm1OsJ3AX8S6YH3AbtjtHMtpjZSGBf4AxJFdkYeCu1K05JJwGvmdlrWRtx66Xj83qsmVURbAVcIelLGR9167U3zjJgLDAbGANMljQ5C+NujXR993QDjgYeTqnXdF000x+754Pg6uDnJCx7FliQsOx14Jq41+XAE8DMXMeQqRgTym4CTsl1LOmOE7gG2EBwrGMTsBX4Sa5jyfC6nJvYRr492rguxwCL4souAy7LdSyZWJfATODuVPvyLTrXiKSOQBXwWELRYwS/oJAk4A5giZktzOoA0yDFGCsk7Rk+7w6MA/6ezXG2VypxmtnlZjbAzAYDlxJ8yfwsqwNthxTX5R5x67IrMAlYk81xtlcqcQLPA30k9ZBUQvCZXZu9UbZPijHGpLzbEnzXpWuqN1AKbExYvhGoDJ8fA5wGTI+b5ntIFsfYXqnEOAh4MjxG9yTwazN7JXtDTItU4tzdpRJjBbA8XJfPAHeZ2fPZG2JatBinBRNPriDY0/Iy8LqZPZjNQbZTSp/X8IfnEcCiVBsuS8foXHExs+UU+I8kM3uOXTd/LApmdkeux5AJZvYPYGSux5ENZvYIBT571sy2Evx4SVlBf1m5NtkE1NH0g1QBfJD94WREMcQIxRFnMcQIxRFnxmL0ROcaMbOdBCeAJ87WmkwwA2q3VwwxQnHEWQwxQnHEmckYfddlEQoPyA8JX5YAAyWNAjab2dvAPGChpOeAp4DvAv2Am3Mx3rYohhihOOIshhihOOLMWYy5nmLqj+w/gAkEJ9AmPu6IqzObYMr5DoJfWeNyPW6PsTjjLIYYiyXOXMXo17p0zjlX0PwYnXPOuYLmic4551xB80TnnHOuoHmic845V9A80TnnnCtonuicc84VNE90zjnnCponOueccwXNE51zzrmC5onOuSIh6WBJtZISL5qbqf5OkrRT0tBs9OdcMp7onMtzkrpJqpdkzTyOTaGpecBTZra4DWO4N+wn6T36FFgnaYukzmZ2H/AK8MvW9udcOvndC5zLf6MBAfeQ/Kaazd4xW9IYgtudTG/jGG4FTgFmARcmqTMRGAz81sy+CJfNB+6UNMLM1rSxb+faxS/q7Fyek3QxcB1wvJk91sY2FgJTgX5mVtOG95cQXFG+S9jGziR9nAkcYWbPh8u6AhuB28xsTlvG7lx7+a5L5/JfFcGtTJrdaktGUhnBltz/RiU5SeWSrpC0RlJ1uOvxAUmHxuqYWT1wB9ALODGijW7ADGB1LMmF7/sMeJJga9C5nPBE51z+Gw28BZRK6p34SOH9VUBX4LnEAkkdgEeBq4AVwEXAL4DhwFOSDourfjtBwp0V0cc3gc4EuzgTrQAqJQ1LYazOpZ0fo3Muj4W7/g4g+FH6UUSV9wnuwNyc4eHfNyPKLiC4GeZUM1sU1+9NwGrg2rAcM1sn6XHgeEl9zez9uHZmATuBuyP6iPU7Ani1hbE6l3ae6JzLb6MIktx84MGI8k9SaGPv8O/miLIzCZLPyoitw8XA2eEMytjkkluBScBZhLMpwy21o4A/mtmmiD4+Dv/2SWGszqWdJzrn8ltV+Pd+M1vSxjZiM84UUXYQwS7HqK3FmN7AhvD5n4EtBFtwsdMGvh3+vS3J+2P9+sw3lxOe6JzLb6PDv+2Zmh9LYj0jykRwrtvFKbwfM6uWdA8wW9LRwLPATOAdYFGS98f6bS6ZOpcxnuicy29VwCYz29iONlaHf6OuUPI6wa7NJeHMylTcCswm2KrrCVQCVzfz/iEJ43Auq3zWpXN5SlIXYBjt25oDeAn4lOA4WqK7CBJV5BadpIrEZWb2IrAKOA04n2CXZLLdloT9bjSzv7du2M6lh2/ROZe/RgKlAJLOTFLnITNrdkKKmdVJ+jMwXVK5me2IK55PcMWUuZImAUsIkuJA4DigmuCKJ4luBX5NcBL6UjP7R1Tf4azRsTSfCJ3LKL8yinN5StL5wG+aqWJADzPbmkJbRxAcTzvFzP6UUFZGsCtyJrtORXiP4Ly7O6OuxiKpR1inE3CWmS1M0u/ZBCeaH2JmvuvS5YQnOueKhKRHgT3MbGwW+3wRWG9mJ2erT+cS+TE654rHJcAYSVOy0Zmk6cDBwA+z0Z9zyfgWnXPOuYLmW3TOOecKmic655xzBc0TnXPOuYLmic4551xB80TnnHOuoHmic845V9A80TnnnCtonuicc84VtP8HzO6IFm39JpMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# set up a class for the standard calculation \n",
    "s1 = alpro.Survival(\"1275b\")\n",
    "s1.init_model()\n",
    "s1.set_params(1e-13 * 1e-9, 1e-13)\n",
    "s1.set_coherence_r0(None)\n",
    "s1.domain.create_box_array(1800.0, 0, 10.0, r0=0)\n",
    "\n",
    "# set up a class for the Fourier-like calculation \n",
    "s2 = alpro.Survival(\"1275b\")\n",
    "s2.init_model()\n",
    "s2.set_params(1e-13 * 1e-9, 1e-13)\n",
    "s2.set_coherence_r0(None)\n",
    "s2.domain.create_box_array(1800.0, 0, 10.0, r0=0)\n",
    "\n",
    "N = 200000\n",
    "pad_factor = 40\n",
    "f_res = N//pad_factor\n",
    "E, P2 = s2.propagate_fourier(s2.domain, pol=\"both\", mode=\"auto\", N=N, f_res = f_res)\n",
    "_ = s2.default_plot(mode=\"conversion\")\n",
    "\n",
    "P1, _ = s1.propagate(s1.domain, energies=E[::10], pol=\"both\")\n",
    "_ = plt.plot(s1.energies, P1, ls=\"--\", c=\"C1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Fourier scheme can be used for many more types of models and the code can also be used to inspect autocorrelation functions and so on. I will endeavour to include more detail here in time. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
