{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boson cloud timescales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import gwaxion\n",
    "from matplotlib import ticker\n",
    "\n",
    "DAYSID_SI = 86164.09053133354\n",
    "YRSID_SI = 31558149.7635456"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will look at the two key timescales governing the evolution of a boson cloud around a black hole (BH): the superradiance instability timescale, and the gravitational-wave dissipation timescale.\n",
    "\n",
    "For concreteness, let's look at a BH consistent with the GW150914 remnant, with $M = 60\\, M_\\odot$ and $\\chi = 0.7$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mbh = 60\n",
    "chi = 0.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first focus on the $e$-folding time of the cloud superradiant growth. We can obtain a quick estimate using closed-form approximations from [arXiv:1706.06311](https://arxiv.org/abs/1706.06311), e.g. for two example values of the boson mass (i.e. two values of $\\alpha$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Instability timescale\n",
      "---------------------\n",
      "alpha = 0.100 :\t231.4 days\n",
      "alpha = 0.176 :\t1.4 days\n"
     ]
    }
   ],
   "source": [
    "print('Instability timescale\\n---------------------')\n",
    "for a in [0.1, 0.176]:\n",
    "    print('alpha = %.3f :\\t%.1f days' % (a, gwaxion.tinst_approx(mbh, a, chi) / DAYSID_SI))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare this analytic approximation to numerical results, in which we evolve the cloud using differential equations as in [arXiv:1411.0686](http://stacks.iop.org/0264-9381/32/i=13/a=134001?key=crossref.8a4e2e36ca36fb28e1d0a6868573e646).\n",
    "\n",
    "Note that `gwaxion` provides both `amplitude_growth_time` and `number_growth_time`, which respectively refer to the field amplitude and occupation number, and differ by a factor of two."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/richardbrito/opt/anaconda3/envs/gwaxion/lib/python3.7/site-packages/gwaxion/physics.py:424: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  self.reduced_compton_wavelength = HBAR_SI / (mass*C_SI)\n",
      "/Users/richardbrito/opt/anaconda3/envs/gwaxion/lib/python3.7/site-packages/gwaxion/physics.py:1126: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  self.nr)\n"
     ]
    }
   ],
   "source": [
    "# create an array of alphas\n",
    "alphas = np.linspace(0, 0.2, 100)\n",
    "\n",
    "# for each of those values, compute the instability timescale numerically\n",
    "tinsts_amp = []\n",
    "tinsts_num = []\n",
    "for a in alphas:\n",
    "    bhb = gwaxion.BlackHoleBoson.from_parameters(m_bh=mbh, chi_bh=chi, alpha=a)\n",
    "    c = bhb.cloud(1, 1, 0)\n",
    "    tinsts_amp.append(c.amplitude_growth_time)\n",
    "    tinsts_num.append(c.number_growth_time)\n",
    "#     tgws.append(c.get_life_time())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/richardbrito/opt/anaconda3/envs/gwaxion/lib/python3.7/site-packages/gwaxion/physics.py:193: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  t =  27. * DAYSID_SI * (m/(10*MSUN_SI)) * (0.1/alpha)**9 / chi\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEGCAYAAABcolNbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABKD0lEQVR4nO3dd3xUVdrA8d+ZyZT0HkgIKfQWSEITG0VBLNhdddUVyyIqrvqqu+i7q6y76+padm1rRdG1YUNdy2tDxU4RpBMINZBCei8zc94/bmYyCQkkksmkPN/PZz4zc+tzL5c8c8499xyltUYIIYTwFZO/AxBCCNG7SaIRQgjhU5JohBBC+JQkGiGEED4liUYIIYRPBfg7AF+JiYnRKSkp/g5DCCF6lDVr1hRqrWM7c5u9NtGkpKSwevVqf4chhBA9ilJqT2dvU6rOhBBC+JQkGiGEED4liUYIIYRP9dp7NEL4S0NDAzk5OdTW1vo7FCHaZLfbSUxMxGKx+HxfkmiE6GQ5OTmEhoaSkpKCUsrf4QhxCK01RUVF5OTkkJqa6vP9SdWZEJ2straW6OhoSTKi21JKER0d3WWlbkk0QviAJBnR3XXlNSqJRgghhE/1/kRTUwovnQfPn+7vSIToMnl5eVx00UUMHjyYUaNGcdppp5GVlcXu3bsZM2aMz/a7ZMkSFixY0Oq8N954g5EjRzJ9+nSf7b87ee+997j33ns7ZVv/+te/qK6u9nw/7bTTKC0t7ZRtd4Xen2isIZC9HPZ8A446f0cjhM9prTnnnHOYNm0a2dnZbN68mXvuuYf8/Hy/xrV48WL+/e9/88UXX7RreYfD4eOIOsbpdHZo+TPPPJOFCxd2yr5bJpoPP/yQiIiITtl2V+j1iWbnnr3U2WIAKN+/zc/RCOF7X3zxBRaLhfnz53umpaenc8IJJzRbrra2liuuuIK0tDQyMjI8CWDTpk1MmjSJ9PR0xo4dy/bt2wF46aWXPNOvueYazx/e559/nmHDhjF16lS+/fbbVmO6++67+eabb5g/fz633XZbm/tesmQJF1xwAXPmzGHWrFksWbKEs88+mzlz5pCamspjjz3GQw89REZGBscccwzFxcWH7Ou///0vkydPJiMjg5NPPtmTYBctWsRll13GjBkzGDp0KM888wwAX375JSeeeCLnnHMOo0aNYv78+bhcLgBCQkK48847mTx5Mt9//z0PPfQQY8aMYcyYMfzrX/8C4KGHHuLKK68EYMOGDYwZM4bq6upmpbu5c+dy7bXXMn36dAYNGsRXX33FlVdeyciRI5k7d64n9muvvZYJEyYwevRo7rrrLgAeeeQRDhw4wPTp0z2lwZSUFAoLCz37bxnT7t27GTlyJL/97W8ZPXo0s2bNoqampu2Lxsd6ffPmkvIqTK5IUijAUbQbksf6OyTRh3z55Zc+2/a0adNanb5x40bGjx9/xPUff/xxwPjjuHXrVmbNmkVWVhZPPvkkN954I5dccgn19fU4nU62bNnC0qVL+fbbb7FYLFx33XW8/PLLzJw5k7vuuos1a9YQHh7O9OnTycjIOGRfd955J8uXL+eBBx5gwoQJPPjgg63uG+D7779n/fr1REVFsWTJEjZu3MjatWupra1lyJAh3Hfffaxdu5abb76ZF198kZtuuqnZvo4//nh++OEHlFI8++yz/OMf//Dsb/369fzwww9UVVWRkZHB6acbVeorV65k8+bNJCcnM3v2bN5++23OP/98qqqqGDNmDHfffTdr1qzh+eef58cff0RrzeTJk5k6dSo33XQT06ZNY9myZfztb3/jqaeeIigo6JBzUFJSwvLly3nvvfeYM2cO3377Lc8++ywTJ05k3bp1pKen87e//Y2oqCicTicnnXQS69ev53e/+x0PPfQQX3zxBTExMc222VZMkZGRbN++nVdffZVnnnmGX/3qV7z11ltceumlR7wufKHXl2jsdhvra6IAcJbs9XM0QnQf33zzDZdddhkAI0aMIDk5maysLKZMmcI999zDfffdx549ewgMDOTzzz9nzZo1TJw4kfT0dD7//HN27tzJjz/+yLRp04iNjcVqtXLhhRce1b4BZs6cSVRUlGfZ6dOnExoaSmxsLOHh4cyZMweAtLQ0du/efci2c3JyOOWUU0hLS+P+++9n06ZNnnlnnXUWgYGBxMTEMH36dFauXAnApEmTGDRoEGazmYsvvphvvvkGALPZzHnnneeJ+ZxzziE4OJiQkBDOPfdcvv76a0wmE0uWLOGyyy5j6tSpHHfcca0e85w5c1BKkZaWRr9+/UhLS8NkMjF69GjPcbz++utkZmaSkZHBpk2b2Lx58xHPY2sxAaSmppKeng7A+PHjWz1XXaXXJ5rwkCAKTcavgJpCSTSi9xs9ejRr1qw54nJa61an//rXv+a9994jMDCQU045heXLl6O15vLLL2fdunWsW7eObdu2sWjRIqD1ZrJOp5P09HTS09O58847271vgODg4GbfbTab57PJZPJ8N5lMrd7HueGGG1iwYAEbNmzgqaeeavasSMtY3d/bmm632zGbzUeMefv27YSEhHDgwIE2l/GOu+UxORwOdu3axQMPPMDnn3/O+vXrOf3004/4nMvhYvLeh9ls9us9r15fdWaz2agIiAYnNJTm+Dsc0ce0Vb3lSzNmzOCOO+7gmWee4be//S0Aq1atorq6muTkZM9yJ554Ii+//DIzZswgKyuLvXv3Mnz4cHbu3MmgQYP43e9+x86dO1m/fj2zZs3irLPO4uabbyYuLo7i4mIqKiqYPHkyN954I0VFRYSFhfHGG28wbtw4zGYz69atazPGtvb9008/HfXxl5WVMWDAAABeeOGFZvPeffddbr/9dqqqqvjyyy+59957ycrKYuXKlezatYvk5GSWLl3KvHnzWo157ty5LFy4EK01y5Yt4z//+Q9lZWXceOONrFixggULFvDmm29y/vnndzju8vJygoODCQ8PJz8/n48++shz/YSGhlJRUXFI1VlbMXU3vb5EY7PZyLcP4g3HiWy3jPB3OEL4nFKKZcuW8emnnzJ48GBGjx7NokWLSEhIaLbcddddh9PpJC0tjQsvvJAlS5Zgs9lYunQpY8aMIT09na1bt/Kb3/yGUaNG8de//pVZs2YxduxYZs6cSW5uLvHx8SxatIgpU6Zw8sknk5mZ2a4Y29p3Z1i0aBEXXHABJ5xwwiF/mCdNmsTpp5/OMcccw5/+9CfPOZkyZQoLFy5kzJgxpKamcs455xyy3czMTObOncukSZOYPHkyV199NRkZGdx8881cd911DBs2jMWLF7Nw4UIKCgo6HPe4cePIyMhg9OjRXHnllc2q4ObNm8epp556SNPwtmLqbtThil492YQJE/Tq1aupr6/n+sVf8OkeB78eaeeey0/yd2iil9uyZQsjR470dxiihUWLFhESEsKtt97abPqXX37JAw88wPvvv++nyPyntWtVKbVGaz2hM/fT60s0FouF6EDjMA9WO7pd23whhOjtev09GqUUCWE2Mu17SG8opq5yPAERMUdeUQjRq7gbL7Q0bdo0v9xL60t6faIBOHFQGFft+CfBxTmUFRwPkmiEEKLL9PqqMzAaBLh7B3AU7/FzNEII0bf0iURjt9s9icYlD20KIUSX6hOJxmaz8X8HIwCokoc2hRCiS/WJRGO32ylu7B3AWZbr52iE8D2z2ex5Mj89PZ3du3dz7LHHHnG9kJCQVqfPnTuXN998s7PDbJP3/q6++mpPVyz33HNPh7fl66ERxJH1icYANpuNGlsMNIClOs/f4Qjhc4GBgYc8mf/dd9/5J5ij9Oyzz3o+33PPPdxxxx1+jEb8Et2yRKOUOlsp9YxS6l2l1KzGaYOUUouVUh3+WWWz2XAFxQIQXHfwsP0DCdFbeZdW7r//fiZOnMjYsWM93dF701qzYMECRo0axemnn97mk+7PPPMMEydOZNy4cZx33nmeMVPa2y1+SEgIt9xyC5mZmZx00kkcPHjwkH1MmzaN1atXs3DhQmpqakhPT+eSSy45pKTywAMPeJowr1mzhnHjxjFlyhRPL9Vg9MF22223eY79qaee6tA5FL9Ml5VolFLPAWcABVrrMV7TZwMPA2bgWa31vVrrd4B3lFKRwAPAJ1rrncBVvyTRBAQEQGg8s/bdR2p8Px5taMBqtXbGYQlxRCkLP2hz3j3npPHryUkAvPLjXu5YtqHNZXff2/5RYt1/kMHoxXfZsmWeeZ988gnbt29n5cqVaK0588wzWbFiBSeeeKJnmWXLlrFt2zY2bNhAfn4+o0aN8oy54u3cc8/19Kf2xz/+kcWLF3PDDTcA7esWv6qqiszMTB588EHuvvtu/vznP/PYY4+1ekz33nsvjz32mKekdrjeiK+44goeffRRpk6dym233eaZvnjxYsLDw1m1ahV1dXUcd9xxzJo1i9TU1HadV/HLdGXV2RLgMeBF9wSllBl4HJgJ5ACrlFLvaa3dfWP/sXH+UYsLDyZLD8RaZ6aurk4SjejVWqs6c/vkk0/45JNPPH1iVVZWsn379maJZsWKFVx88cWYzWYSEhKYMWNGq9vauHEjf/zjHyktLaWyspJTTjnFM6+1bvEBT7f46enpmEwmz9ACl156Keeee+5RH3tZWRmlpaVMnToVgMsuu4yPPvrIc+zr16/33P8pKytj+/btkmh8rMsSjdZ6hVIqpcXkScCOxtIKSqnXgLOUUluAe4GPtNbt7s5VKTUPmAeQlJTUbF5ChB0oo6jWRW1tLaGhob/4WIToiPaWRH49OclTuvElrTW3334711xzzWGXa637/5bmzp3LO++8w7hx41iyZEmzgd6O1C3+L92nW0BAgGckTMDTpb7Wus3taK159NFHmyVE4Xv+vkczANjn9T2ncdoNwMnA+Uqp+QBKqWil1JNAhlLq9tY2prV+Wms9QWs9ITY2ttm85OgQHklYzmv2eyH7c18cixA9wimnnMJzzz1HZWUlAPv37z/kHsyJJ57Ia6+9htPpJDc31zPUcksVFRXEx8fT0NDAyy+/3OFYXC6Xp3TxyiuvcPzxxx92eYvFQkNDAwD9+vWjoKCAoqIi6urqPJ1iRkREEB4e7hm8zDuuU045hSeeeMKzjaysLKqqqjoct+gYf7c6a+1nh9ZaPwI80mJiETC/leXbJSIkkOjAXBKLf6bw4NZfuhkherxZs2axZcsWpkyZAhg35F966SXi4uI8y5xzzjksX76ctLQ0hg0b5qmGaukvf/kLkydPJjk5mbS0NCoqKjoUS3BwMJs2bWL8+PGEh4ezdOnSwy4/b948xo4dS2ZmJi+//DJ33nknkydPJjU1lREjmoYBef7557nyyisJCgpqVnq5+uqr2b17N5mZmWitiY2N5Z133ulQzKLjunSYgMaqs/fdjQGUUlOARVrrUxq/3w6gtf770e7LPUyAW35+PpUf38PgnUsoGvoroi955mh3IUSrZJiA9gsJCfGUrETX6yvDBKwChiqlUpVSVuAi4D1f7Mhms5HdEA1AbfF+X+xCCCFEK7os0SilXgW+B4YrpXKUUldprR3AAuBjYAvwutZ6ky/2b7PZ+LE80oilUh7aFKI7kNJM39CVrc4ubmP6h8CHvt6/8dBmDJRBSEMhLpcLk8nfBTohhOj9+sxfWpPJhD08lgZtJsxVRn1Vmb9DEkKIPsHfrc66VP/wIN53HYPdauHYqnLsoZH+DkkIIXq9PpVoBkQGcXXD9USZFR8pO+H+DkgIIfqAPlN1BpAQGQRAaZ2msrrGz9EI4Ts5OTmcddZZDB06lMGDB3PjjTdSX1/v15jeeecdT3f/AHfeeSefffaZHyM6sraGTRAd06cSTVhwELF2TabtAPU7v/V3OEL4hNaac889l7PPPpvt27eTlZVFZWUl//u//+vXuFommrvvvpuTTz7ZjxH5Vlvd7PRFfSrR2Gw2Hp9QyJvcyuA1f/F3OEL4xPLly7Hb7VxxxRWAMQjaP//5T5577jmqq6txOp3ceuutpKWlMXbsWB599FEAVq1axbHHHsu4ceOYNGkSFRUVLFmyhAULFni2fcYZZ3j6M2uri//Whg747rvveO+997jttttIT08nOzu72eBmn3/+ORkZGaSlpXHllVdSV1cHQEpKCnfddReZmZmkpaWxdeuhvXosWbKEc889l9mzZzN06FB+//vfe+Z5l0jefPNNzxAF7R3GAGj1GLOzs5k9ezbjx4/nhBNO8MQ1d+5c/ud//ofp06fzhz/84Rf9+/VGfSrRBAYGUhM0EKfJirU6F6qL/R2S6AsWhbf9Wv1803Krnz/8su3k7tLFW1hYGElJSezYsYOnn36aXbt2sXbtWtavX88ll1xCfX09F154IQ8//DA///wzn332GYGBgYfdj7uL/59++ompU6fy5z//GTCGDli1ahU///wzI0eOZPHixRx77LGceeaZ3H///axbt47Bgwd7tlNbW8vcuXNZunQpGzZswOFw8MQTT3jmx8TE8NNPP3HttdfywAMPtBrLunXrPOsvXbqUffv2tbqcN/cwBv/85z+ZM2cON998M5s2bWLDhg2enq/bOsZ58+bx6KOPsmbNGh544AGuu+46z3azsrL47LPPePDBB48YQ1/RpxJNUFAQKsBCZYjRJbgjp90dQwvRY7TVe7F7+meffcb8+fONcZqAqKgotm3bRnx8PBMnTgSMxOSe35aWXfy7O7HcuHEjJ5xwAmlpabz88sts2nT4Z7C3bdtGamoqw4YNA+Dyyy9nxYoVnvnuoQPGjx/f5hg0J510EuHh4djtdkaNGsWePXsOu09ofRgDk8nkGcagrWOsrKzku+++44ILLiA9PZ1rrrmG3NymIeIvuOACzGbzEfffl/SpVmcmk4kSp42PypK5SG2jfvdKAob13jpi0U0sauczWxOuMF5HafTo0bz11lvNppWXl7Nv3z4GDx7caiJqKzm11RV/a9zrH27ogNYcqb9F9xADZrO5zfse3sMQeC/nfUwtY/+lwxi4XC4iIiLaHO8nODj4sMfTF/WpEg1Aar9w1jhSAHDulxKN6H1OOukkqqurefFFY4xBp9PJLbfcwty5cwkKCmLWrFk8+eSTnj+mxcXFjBgxggMHDrBq1SrA6P7f4XCQkpLCunXrcLlc7Nu3j5UrV3r201YX/20NHRAaGtpq784jRoxg9+7d7NixA4D//Oc/bfYW3VH9+vVjy5YtuFyuZqOMtldrxxgWFkZqaipvvPEGYCTKn3/+uVPi7a36XKKJiQgj3z4IAJW/0c/RCNH5lFIsW7aMN954g6FDhzJs2DDsdjv33HMPYHSVn5SUxNixYxk3bhyvvPIKVquVpUuXcsMNNzBu3DhmzpxJbW0txx13HKmpqaSlpXHrrbeSmZnp2Y93F//Lly/nzjvvBJqGDpg5c2azrvsvuugi7r//fjIyMsjOzvZMt9vtPP/881xwwQWe6qv583/xiCDN3HvvvZxxxhnMmDGD+Pj4Dq/f1jG+/PLLLF68mHHjxjF69GjefffdTom3t+rSYQK6UsthAtwqKiq47rmveab4NwSYIGDhbrBJW3nRefrKMAHSxX/P11XDBPSpezRg/EJJirBwWt7fSYpP4BmTDYu/gxJCiF6sz1WdmUwmRvYLZqdOILvCJL/IhPiF5P+OaK8+l2gA0hIjUcCBSheFJdKLs+h8vbVKWvQeXXmN9slEExsVzoIRVXwd9TcG/PdCf4cjehm73U5RUZEkG9Ftaa0pKirCbrd3yf763D0aMG5iTkiKov+eLExVDVBbBnbpy1l0jsTERHJycjzdlQjRHdntdhITE7tkX3020WC2UBmSQljFdpw5azEPmebvsEQvYbFYSE1N9XcYQnQbfbLqzGQyoS2BbHClAFC/Z+XhVxBCCPGL9clEAxAeGsJ/S1MAcO790b/BCCFEL9ZnE01cVDjbbaMBsO7/AZwydoQQQvhCn000ISEhWMPj2enqj9VRCfvX+DskIYTolfpkYwAwEk1ymJknCs5keJiJuREpffdkCCGED3XLEo1SapBSarFS6k2vaUlKqfeUUs8ppRYe7T7MZjPpA4J5wzmNJyqnUlov40cIIYQvdFmiaUwQBUqpjS2mz1ZKbVNK7XAnEK31Tq31VS02MQz4QGt9JTCqM2IaPyiOUAsU1Wp+3pXfGZsUQgjRQleWaJYAs70nKKXMwOPAqRjJ42KlVFtJZC1wkVJqOfBFZwQUHRXFMQkBXBqbTfJPf4esTzpjs0IIIbx0WaLRWq8AiltMngTsaCzB1AOvAWe1sYkrgLu01jOA01tbQCk1Tym1Wim1uj1PZUdERHDZ6EDmxu1kSO5/cax/vd3HI4QQon38fY9mALDP63sOMEApFa2UehLIUErd3jjv/4DfNU7f3drGtNZPa60naK0nxMbGHnHnJpOJiIgIiqOMwZxU9hfgNWytEEKIo+fvhlaHDlIOWmtdBMxvMXEjcH5nBxAVFcWG/EQqAqIIrSmE/I0QP7azdyOEEH2Wv0s0OcBAr++JwIGuDCAqKooPdjn4oNZILq7tn3bl7oUQotfzd6JZBQxVSqUqpazARcB7XRlAYGAgmfF2vnKNA8C1TRoECCFEZ+rK5s2vAt8Dw5VSOUqpq7TWDmAB8DGwBXhda72pq2JqjIvJQ+JYrcbg0CZMB1ZBdcs2C0IIIX6pLrtHo7W+uI3pHwIfdlUcrYmLiSYpJoxlRcczMCaGY1xOf4YjhBC9ir8bA3QLkZGRpMUEcFv+fMY7zLxqDcfq76CEEKKX8Pc9mm4hICCAKanGCJubCp3kHyzyc0RCCNF7SKJpNGJgHEmhJoYGVVP1zROw9iV/hySEEL2CVJ01io2N5fbJdhLKtzJi/cPovFRU+iWgWnvURwghRHtJiaZRUFAQcZFhlESOpc4aiSrZBTmr/B2WEEL0eJJovPTr1w+UmezIE40JP7/m34CEEKIXkETjJS4ujrUFTm7ZewwAeuPb4Kj3c1RCCNGzSaLxYrVaOXZwNNmmZLa4BqJqSyDr//wdlhBC9GiSaFpISYxnfJyZ153TANA/PunfgIQQooeTRNNCTEwMxw+08YZzKj8yhpq0S0Frf4clhBA9liSaFsxmM1OH98NsC+bC2jv4rD5NmjgLIcRRkETTioT4/kyJNx4xWromF5cMhiaEEL+YJJpWREZGctrQIOxmSCaXmqVXwZf3+jssIYTokaRngFYopUgfmsTDZBNbWUvw2rfReyJQx94A1mB/hyeEED2KlGjaEB8fT6DFTHn4cMrChqNqS2HdK/4OSwghehxJNG2wWCz0798fh0vzefAcY+K3j4Cjzr+BCSFEDyOJ5jASExNZV+Dk1l3pZJMIZXthzQv+DksIIXoUSTSHERQUxMkj44i0m7mv/gJj4or7ob7Kv4EJIUQPIonmCJKTBjIz2cInrglsVUPQVQdh51f+DksIIXoMSTRHEBERwanDwwixKG6pvZIPJ70AI07zd1hCCNFjSKI5AqUUIwancPogK5t0Cvets1JbKw0ChBCivSTRtENcXBxzRoYTaVPsLXfx2tebYfunULbf36EJIUS31y0TjVJqkFJqsVLqTa9p05RSXyulnlRKTevieBg+OJULR1j5zSgrU3OehJfPh49v78owhBCiR+qyRKOUek4pVaCU2thi+myl1Dal1A6l1EIArfVOrfVVLTahgUrADuR0TdRNYmJimDkskhlJFvLiT8YVEAib34Xtn3V1KEII0aN0ZYlmCTDbe4JSygw8DpwKjAIuVkqNamP9r7XWpwJ/AP7swzhbpZQiNTUVgDp7LFsH/MqY8eGt0FDb1eEIIUSP0WWJRmu9AihuMXkSsKOxBFMPvAac1cb67i6USwBba8sopeYppVYrpVYfPHiwkyJvEhkZSXh4OF/lNPCrHSeRb02Ckl3w7b86fV9CCNFb+PsezQBgn9f3HGCAUipaKfUkkKGUuh1AKXWuUuop4D/AY61tTGv9tNZ6gtZ6QmxsbKcH6y7VJASbqHQGcFPl5caMrx+CouxO358QQvQG/k40rY0oprXWRVrr+VrrwVrrvzdOfFtrfY3W+kKt9ZddG2aTiIgIjhsezwkDAvjeNZLPLVPBWWckGyGEEIfwd6LJAQZ6fU8EDvgplnYbPHgwF4+yE2yBWysuZv3Ay+D0B/wdlhBCdEv+TjSrgKFKqVSllBW4CHjPzzEdkc1mY8zQVC4YZqWEMH6z9zSKqmUUTiGEaE2HE41SKrixtVhH13sV+B4YrpTKUUpdpbV2AAuAj4EtwOta600d3bY/JCYmMntYGIPCTZTWaRa9/RPUlsPyv8lQAkII4eWII2wqpUwYJY1LgIlAHWBTSh0EPgSe1lpvP9J2tNYXtzH9w8bt9Cgmk4nhw4ZxRcFalmys45jYBupfvgjrvm+NezYz7/Z3iEII0S20p0TzBTAYuB3or7UeqLWOA04AfgDuVUpd6sMYu63IyEgmDInnj8fYGRBiYmv8eWhlgm8fhqyP/R2eEEJ0C+1JNCdrrf+itV7v9SwLWutirfVbWuvzgKW+C7F7GzJkCDab8VhPUeAg1idfacx4+7dQstt/gQkhRDdxxESjtW4AUEpdoJQKbfz8J6XU20qpTO9l+iKLxcLw4cMBeG5jPWdvnca+mBOgtgyWXgYNNX6OUAgh/KsjjQH+pLWuUEodD8wCXgCe8E1YPUt0dDTx8fEMiTShMXFu3m+oD0uGvPVGFzVCCNGHdSTROBvfTwee0Fq/C1g7P6SeafDgwcwcFMyEfmYOOoK5puZ6dIAdLEHgkqbPQoi+qyOJZn9jFzC/Aj5UStk6uH6vFhAQwKhRo7gqzUb/YMUXFYn8Mfph9Kn/AJOcJiFE39WRv4C/wnjeZbbWuhSIAm7zRVA9VUREBCOHpHJDuh2rGV7eE8qTn282ZlYVQcEW/wYohBB+cMREo5RSAFrr6sb+xrY3fs/VWn/ivYyA5ORk0pJjuHKM0RJt2Zp91OTtgMUnw3/OkVE5hRB9Trueo1FK3aCUSvKeqJSyKqVmKKVeAC73TXg9j1KKESNGMDUlmN+mWbllgo0tew+iQ/pDRS68dK5RuhFCiD6iPYlmNkZDgFeVUgeUUpuVUjuB7cDFwD+11kt8GGOPY7VaGTVqFMcnWrGZFeXVdWxJ/xPOmBFwcKuRbGrL/R2mEEJ0CaW1bv/CSlmAGKCm8T5NtzVhwgS9evVqv8Zw4MABsrKyqHdqnl5fR7CuYonpblTpbkg6Fi59C6xBfo1RCCG8KaXWaK0ndOY2O9QcSmvd0HhvprQzg+itEhISGDBgAFUNmuxSF18V2Lk96C50aALs/Q6WXgpOh7/DFEIIn5J2tz42ZMgQBifEcNN4G4EB8NpOG3+P/is6KAYGjAfzEfs1FUKIHk0SjY8ppRg1ahQj+4dyU6Ydqwme3mrnvuSn0dNu93d4Qgjhc+1ONEqp73wZSG8WEBDA2LFjGRsfxIIMG2YFT66t5R//1/hcTVkOvHmV0T+aEEL0Mh0p0dhbTlBKndCJsfRqdrudsWPHkhlvZ/44Gwr4eksOlVXV8O71sPFNeP40KM/1d6hCCNGpOnKDYLhSahmwCdgI5APPYoxVI9ohODiYsWPH4nKtw26GYVFmNm5YT/opD2B//WLI3wiLZxmt0WKH+TtcIYToFB0p0ewC7gGygfHA1cCffRFUbxYWFsaYMWMYG2fBZlbU1tayOruQ18c8iU6cCGV74blZsPdHf4cqhBCdoiOJpl5rvUpr/bzW+jat9SVa6xd9FlkvFhUVxZgxY3D33PPoyjJ+/3E+C4PuxjX0FKgpgRfmwPo3/BypEEIcvY4kmqk+i6IPio6O9iSb4wYEYDXD0vUlXFn9OxoyrwBnHVTm+ztMIYQ4au1ONFrrCl8G0he5k824OAt/mGgnxAJfZpdx7u5zKbvgTZhyfdPCHejBQQghuhN5jsbP3MlmaJSF/50cSEygYsOBCs74r4kdB6uMhYqyYfFMOLjNv8EKIcQvIImmG4iOjmbcuHEMjLDyx8l2ksNM7Cut5dFPG8ey+eJvkLMKnp4OG9/2b7BCCNFB3TbRKKUGKaUWK6Xe9Jp2tlLqGaXUu0qpWf6Mr7OFh4eTnp5OXJidOybZmTPIwun9q9i3bx96ziMw5nxoqII3r4APb4OGWn+HLIQQ7dKliUYp9ZxSqkAptbHF9NlKqW1KqR1KqYUAWuudWuurvJfTWr+jtf4tMBe4sMsC7yIhISFkZGQQERrEecOsWM2K7Oxs1u/I4R/Bt1I/814wWWDl0/DsyXAwy98hCyHEEXV1iWYJxvg2HkopM/A4cCowCrhYKTXqCNv5Y+M6vU5gYCCZmZmEh4d7pt2/fA///monZ64cxYHz3oOoQZC/wXjepk7aaAghurcuTTRa6xVAcYvJk4AdjSWYeuA14KzW1leG+4CPtNY/tTJ/nlJqtVJq9cGDBzs7/C5jsVgYN24ccXFxAMxOsdAvSLE1v5JTXi/n0xPegLEXwrQ7wBbq52iFEOLwusM9mgHAPq/vOcAApVS0UupJIEMp5e7m+AbgZOB8pdT8lhvSWj+ttZ6gtZ4QGxvr88B9yWQyMXLkSFJSUkgMNXHXlEDG9zNTUefkt0u38UfTDdRmXNm0wsa3YNMy/wUshBBt6A6DoahWpmmtdREwv8XER4BHuiSqbkApRUpKCiEhIWzZsoUF6fDpHgevb6vnpR/2snJnMW9ffxwh9YXw35ugrtxINqfeD6H9/B2+EEIA3aNEkwMM9PqeCBzwUyzdUkxMDJmZmQQFBTErxcKfptjpH6QYGOwkQDsgpB+cfBdYgmHzu/D4JFj7sjzkKYToFpTu4j9GSqkU4H2t9ZjG7wFAFnASsB9YBfxaa73paPYzYcIEvXr16qOMtntxOBxs3bqVwsJCah0apSDYZmH48OEcbLARUZ9Lwjd3wI7PjBVST4TTHpSeoIUQ7aaUWqO1ntCZ2+zq5s2vAt9jDDmQo5S6SmvtABYAHwNbgNePNsn0VgEBAYwePZrBgwcTaDFhMyscDgc/rd/IvCU/MPO5Xfxn0AO4zn4KAqNg1wp49zop2Qgh/KrLSzRdpTeWaLyVlZWxefNm6urqqGrQPL+xjtX5TgAmpURx32kDSF13P2ReDomNP07qq8ESCKq122JCCNELSjSi84SHhzNx4kTi4uIItiiuT7dxfbqNcKti5e5iTnl6E48G/466/hlNK71zrTH8QL4UGIUQXUcSTQ8WEBDAyJEjGTFiBAEBAUzsH8Dfjg/khAEB1DtcPPhpFpctXonWGioPwq6vYPfX8OTx8N8boUKGIRBC+J4kmh5OKUX//v2ZMGEC4eHhhFgVV6XZ+MNEO/HBipNT7EaiCYmFG36CSdcACtYsgUcy4Mt7oa7S34chhOjF5B5NL6K1Zv/+/ezcuROXy4XDpTErCA4OZtiwYbz000FcGq4Z5cD25V9g2wfGimGJcMMasNj9ewBCCL/zxT2a7vDApugkSikSExOJiopi27ZtlJWVAVBdXc2KH9fy8IoaGpyapasC+cOpDzBnynWoT++EhIymJONygXaBWS4NIUTnkKqzXigoKIj09HSGDRuG2WwGIMym+J9MG0mhJvaX1vC7V9dy9vvw44zXYdZfm1bevAwenwjrXgWnw09HIIToTaTqrJerq6tj+/btFBYWAuDSmhU5Dt7JdlBa6wLg5JFxPPbrTOwWM7xyEWR9ZKwcNRhOuAXG/grMFn8dghCiC0nzZtFhNpuNMWPGkJaWht1ux6QU0wZauPd4O2cPsWAPUNQ1OI0kA3DhS3D2ExCZCsXZxgOfj2TCymegoca/ByOE6JGkIr6PiI6OJiIigr1797J3717sAXD2ECvTB1pwUMeePXtITExkw4FKXt0xhgWXrGBgzgfwzUNQmAUf3gomM0y48sg7E0IIL1J11gfV1NSQnZ3tqU5zs9lsPLbeybe7yggwKc4fn8h1UweRVPA5rH4OLnoVrEHGwlkfQ+wIiEz2wxEIIXzFF1Vnkmj6sOLiYrKzs6mqqvJMy6ty8dEezdf7anFpMJsUZ45L4Lppgxnar3GQtboKeGg01FfAiDNgyvUwcLJ0bSNELyCJpgMk0bSP1prc3Fx2795NfX29Z3pelYuP98GKPTU4Gy+RBy8Yx3njE6EiDz69yxhszdVgzIxPh8nXwOhz5XkcIXowSTQdIImmYxwOB/v27SMnJwen0+mZfrDaxfJcE1/vreXzW6bRL8xIIvtLa4hXJZhWLzaq1WoaR+gOjIJrv4OweH8chhDiKEmi6QBJNL9MfX09e/bs4cCBA3hfG7UOTUJcNCkpKQSHhDLjwS8xKcUVx6VwXlo0wdvfg5VPgTLDvC+aNrj7Gxh4jDwAKkQPIYmmAyTRHJ2amhr27NlDXl7eIfNqA0L584oScsvrAAizB3DRpCQunZREUnA9BEYaCxZsgX8fAyH9IeMSyLgMolK78jCEEB0kiaYDJNF0jurqanbv3k1BQUGz6U6XZluVnY93N/Dz/grAaAtw0og47j1vLDEhNmPgtfdvhqIdTSumngjpl8LIOU0t2IQQ3YYkmg6QRNO5qqqq2Lt3L/n5hw4tkFdv46tcxedZJcSEWFnx++kEmI1ngWvqHATm/gg/vQib3wFHrbFSaALctEGq1IToZqRTTeE3wcHBjBw5kuTkZE/Ccf9I6W+t48JkODslHFdQNApjemFlHdPu/5LpI+K4aOLfmDL7Xkybl8HalyFuRFOSaaiBbx+GMedDzBB/HaIQwkekRCN+kdraWvbt20dubi4ul6vZvICAABISEvi5JID/eXMD7kssKSqIX01I5NzMRBJCA5r6T9vwJrx1lfE5Ph3SLoAx50nLNSH8QKrOOkASTdeor69n//79HDhwgIaGhmbzlFLowEh+PKh4d30BB8pqG6fDiUNjWXz5BKOKLfdn+PEp2Pye8RCosRQkHwujz4EJV4FJuuUToitIoukASTRdy+l0kpeXR05ODjU1h3a+GRQSwn5HKJ9lV/Lp5gKOHRLNkismAeByaVbtLmbigEBM2Z/Chjcg6xNw1kH/sTD/66YNVRVCcExXHZYQfY4kmg6QROMfWmuKiorIycmhtLT0kPkWi4WgiFgCI6IZMSAagO+zi7j4mR+ID7czZ1wCZ6UnMCoK1Lb/A0sgjDrTWLlgq9FcOmkKjDzD6P5G+loTolP16USjlBoFLAKKgM+11m8ebnlJNP5XUVHB/v37KSgoOOQ+DkBUVBQJCQmszmvgL+9vYX9pU0locGwwZ4xNYM64eIbENfaxtuFNeOdacDZ1lUP8OBh+Oow4DfqNkf7WhDhKvS7RKKWeA84ACrTWY7ymzwYeBszAs1rre5VStwArtdZfK6Xe01qfebhtS6LpPhoaGsjNzWX//v3U1dUdMt9qtdKvXz/yncF8tLmQDzbkUlxlJJP4cDvf/mEGJlNjAqkth+2fwJb/wvZPoaGxQ1BbGNyWDQFW47vLaQxrIITokN6YaE4EKoEX3YlGKWUGsoCZQA6wCrgYKATuAqqBY7XWxx1u25Jouh93tdqBAwcoLi5udZmIiAhi4/qRVa74cEM+iZFB3HjyUAAOVtRx2eIfmTW6P7NH92dkrAW1awVs/QAC7HDaP4yN1FfBv8ZC8hQYegoMnQmh/bvqMIXo0XpdogFQSqUA73slminAIq31KY3fbwfQWv+98bsZeFtrfVYr25oHzANISkoav2fPni45BtFxNTU15OXlkZub26zXaLeAgABiY2Pp378/YWFhKKV45ce93LFsg2eZpKggThndj5mj+jM+ORKzu9ST/QX85+zmG+w/1kg4Q06GxEnyoKgQbegrieZ8YLbW+urG75cBk4EHgDuAYOAJrfU3h9uulGh6BncpJzc3l6KiolaXCQwMpF+/fkTGxLJ2fxUfb8rn0815FFY2Jah+YTa+/v0MrAGNzaBL9xpVa9s/gZ1fgcOrJdxNGyFioPG5tgzs4b46PCF6nL7SM0Brd3O11no3jaUV0XsopYiJiSEmJoa6ujry8vLIy8tr1kS6pqaG3bt3s3v3bsLDw7nhmH7cdfpw1h+o5JNNeXy6JZ+BkUGeJON0aX73QSGTB53E9NkXMzDUBHu+gR2fQ/GupiQD8MwM437O4OkwaBqknABBUV18FoTo3bpjoskBvP4SkAgc8FMsogvZbDaSk5NJSkqivLycvLw8CgoKmo2PU1ZWRllZGUptJyoqinmT+rFw9jBqHE0l87V7S/hgQy4fbMgFNjEkLoTpwxOYNvwWJqREYnMvWF0MVQeNUs3qXca4OiijJVvqiUZv07HDuvIUCNErdceqswCMxgAnAfsxGgP8Wmu9qSPblaqz3sHpdFJUVEReXh4lJSW0dr2aTCZiYmKIi4sjKiqKilonn2zO44ttBXydVUhFncOzbKDFzCc3n8jAqMaeo50OyF0HO78wqtj2/tA0auilbxn3dAD2rTQaGQycBNZgHx+1EP7T66rOlFKvAtOAGKVUDnCX1nqxUmoB8DFG8+bnOppkRO9hNpuJi4sjLi6O+vp6CgoKKCgooLy83LOMy+XyTDebzcTExDA9JZbzMjNwalizp4QvthXw1baDFFfVkxgZ6Fn3+qXrCbNbOW7IJRx7wY1EWRyw7wdjiIOkKU2BfPswbH0fTAGQkGF0j5N8HAycDIERXXhGhOh5/F6i8RUp0fRuNTU15OfnU1BQQHV1davLmM1moqOjiY2NJSoqCrPZTGWdgxCb8fuqrLqBjL98gqvxv4BSMDohjGMHx3Ds4GgmpkQR3LgsX/3DaEadtx6098OnCibNa2parbU8NCp6tF7Z6sxXJNH0DVprqqqqOHjwIAUFBa32swZG9VpUVBQxMTFER0djsVhwuTQbD5TxzY5Cvt1RyKrdJdQ7mpLIIxdncOa4BMB4hifEFkCgqwr2/Qh7voU938P+NXDSnXDc74yVdn4J71wPSZON0s7ASUaPBe6eqoXo5iTRdIAkmr5Ha01lZSUHDx7k4MGDbSYdpRTh4eGe1m52ux2A2gYnq3YX8312Ed9mF/HsbyYQG2o0HbjtjZ95Z91+xiZGMDk1ismDohmfHEmIqQFcDrA1dpOz4gFY/pfmOwwIhAGZkDgBZtwpz/CIbk0STQdIounbvEs6hYWFVFVVtblsSEgI0dHRREdHExoaimql6mvei6v5dEs+3v9dTApGJ4TzqwmJXDYlxZjocsHBrZCz0mhAsPcHKM425oUnwc1ND5zywS0QPtBIQvHpYA87+gMX4ij1usYAQviKUoqQkBBCQkJITU2lurqawsJCDh48SEVFRbNlKysrqaysZM+ePVitVqKjo4mKiiIyMpKAAOO/yNO/mUBZTQNr9hTz485ifthVzKb9ZWzYX8bUYbGebW0tqOSZFQ2MT57B+GPOY+iZIZhqimH/aqivbNppVRGsetY7YogZZjQ0GJAJw09r/ryPED2YlGhEn1NXV0dRURGFhYVtNpkGI1lFRER4Ek9gYGCz0k51vYO1e0tJiAgkNcZo8vzcN7u4+/3NnmVC7QGkD4wgIymSzKQIThgaa3SVU1dhDPR24CfjPk/exqZm1QCXvAVDG5tWb/8Uinca3ej0TwNbSOefFCEaSdVZB0iiEe3hcDgoKSmhsLCQ4uLiQ0YJ9Wa324mKivKUdszmQ3uH3lVYxZfbCvhpbyk/7SlpNvRBeKCFdXfO9CSr934+QHJUECPiQ7HhgPyNcGAt7F8LM++GYGO8Ht6YC5uWNW5FQfTgpqSTdIzR1FqITiKJpgMk0YiO0lpTXl5OUVERRUVFh72v425Q4E46ISEhrd7byS2rYe3eUtbuLUEpxR2njQSMhgej7/oYp0tjNZsYGR/K2MQIxiaGMzYxgiFxIU2dhG5402jNlvszFGxpXvIZfjpc/IrxuaYUvn4A+qVBv9FGVZx72AQh2kkSTQdIohFHq7a2luLiYoqLiykpKWnWFU5LFouFyMhIz8vdkq0tByvq+PuHW/g5p5Tsg4cmtMd+ncEZY42m1TsPVtLg1AyODSZAO4zGBnnrIXe90V1OxiXGSru+hhfOaNqIKcBINnEjIW4UjL+iqZQkRBsk0XSAJBrRmVwuF2VlZZ7Ec7jSDkBQUBARERFERkYSERGBxdL2czTltQ1syCljfU4ZG/aXsj6njFd/e4ynm5yFb63ntVX7sAWYGBEfxugE9yuc4f1CCbQ2VuEVZRuln/yNkL/JuK+D1//vW7IgtJ/x+fO/QEUuxI4wElHscAhLBJPpaE6T6AUk0XSAJBrhS3V1dZSUlHhKO4e7twNGE2p30gkPD/e0ZmuPez/ayocbctlbfGgPCMcPieGlqycbMTmcfJddxMj+YfQLs6EaqqFgKxRsMpLQyYuaei14fLJRMvJmCYaYoZB+CUxu7Cjd6QC0PHDah0ii6QBJNKKruB8ULSkpoaSkhLKyMlwuV5vLK6UIDQ0lIiLCk3haa1jQUll1A5sOlLE5t5xNB8rZfKCcaSNiuf1U477Phpwy5jxmDNMUEWRhRP9QRvQPY3j/UIb3D2V0Qhi2gMb97P3BKPkUbDUSzsFtUFVgzDvx9zDjf43PO7+Cl86FyBSjGi56iPGKGQrRQyE4Rrrc6WUk0XSAJBrhL+5qtpKSEkpLS5t1ANoad+IJDw/vcIlHa+1phPDT3hLu+2grW/MqKKs5tIT11W3TSI42mmF/uCGX2gYnw/qFMjg2xKh+qy42Ek5InNGyDWDdq/DOtTSrgvN2286m+z7r3zD6gYsaZKwfGClJqAeSRNMBkmhEd+FwOJolnsrKyiOuExISQnh4uCf5WK3tbz2mtSavvJatuRVszatgW145uwqrWHbdcZgaW7Kd9dg3/JxTBhi5YGBkEEPjQhgSF8IJQ2M5fmhM0wYbaoyqt6LtULij8X071BTDjT83LffYRCjMavpuC4eoVOM16mwYfbYx3ekwdmo6cilOdD3pGUCIHiggIMDTxQ1AQ0MDpaWlnldrDQvcvRXs378fMIazdieesLAwgoKCWm1ODUYJKT48kPjwQKaPiGt1mVPT4hkQGUhWfiW7C6vYW1zN3uJqPt9agEtrT6LZuL+Mv36wmUGxIQyOTWdQ/+MYNCaYxMigpubXbmPOM5pfF2dD8W6oKzPG+sldZzQ6cNv1JbxyEUQkQWQyRCQ3vicZn+PHyT2hXkYSjRBdzGKxEBsbS2ys0XVNQ0MDZWVllJaWUlZWdkgXOWAMi1BTU0NeXh5gJC/vxBMaGtqu+zxu86cO9nyud7jYU1TFjoJKthdUMiEl0jNva14FP+ws5oedxc3Wt5pNJEcH8fo1U4gMNkpb20deT2iGxWiIAFBVaLR8K9llPNfjVpFnPAtUnN3UD5y3hfuaEs1nf4bqIqM7nvCBEDYAwhMhLAECbIeuK7olqToToptxOByUl5d7kk9FRcVhGxdAU99uYWFhhIWFER4ejs1ma7PU017FVfVs2F9GdkEl2Qcr2Xmwip2FleSX12G3mNj859me6rgzHv2ajfvLCbSYSY4OIjk6iJToYJKig8hMimRkvFenofXVULoHSvY0f68phSs+aFru0fFQtKP14I65HmbfY3wuy4ENb0BoAoTFN73LaKjt56iDAJtUnQnRFwQEBHi6ugGjcUFFRQVlZWWUlZVRXl5+SHNqrTUVFRVUVFR4qtusVqsn8fySUg9AVLCVqcNim3UcClBZ5yCvrMaTZIxlbUQFWymuqmdrnnF/yG3+1MGeRPPzvlLu/3gbA6OCGBg1hIGRYxk4OojEyECig600S42n3gfFu4xE4v2qyG0+smn+Jvhs0aEHYAuD0P5w2TsQPsCYlvUJ1JZBSCyE9DNe9oi++QxRZQFkfwE/vWA03rjoZZ/sRhKNEN2cyWTyVJOBkVRqamo8SaesrKzVUUbr6+spLCyksLDQMy04ONiTdMLCwggODv5FpZ4QWwBD4kKbTXvxykmA0Qx7T3EVu4uq2VNYxZ7iaiZ6Vcdl5VfwzY5CWhNoMfP97TOICDKq4z6sGU2DdSSJwwJJiAgkLtRu3BtyOpp3xRM2AKYsgPIDRhIqP2BU0dWVGy/vIRi+f9QYqtubKQCCYowGC6feZ0yrKTF62A6KbnoFRhl/kIOiembVXfkBoweJ3HXGOcjf2DTPHm509uoDkmiE6GGUUgQFBREUFER8fDxg3OcpLy9v9mqty5yqqiqqqqrIzc0FjCTmrnILDQ0lNDT0kF6qOyo8yMLYoAjGJka0On/a8DienzuRvcXV7CuuZl9JNfuKa9hXUk2D00V4YFNDgEc+396sZGQ2KfqH2YkPt3NWeoJnHKCKiOFkj7qN/mF2YkNtRjLS2kgWFXlNA9MBDJ5hJI3Kg1CZb/yqryuDyjxo8ErYpftg+V/bPtArPzY6NQX48SnIXm6UoOzhxv7cr/BEGH5q03q568ESCAF2491sNZKW2dr+5uAthwwv3gW1pVBbDtWFxjAU1YVQth+Gzmxq8bf7W1g2r2m9gEBInmIMSzH2wubnqRNJohGiF7BYLM1atmmtqa6ubpZ4Wmvd5nK5PPPdAgICPEknNDSUkJAQ7Hb7Ud/vcYsNtbXZGq6yztFsP7NG9yc1JpgDpTXsL62lsLKO/aU17C+tYfKgKM9yP+8r49LFPwLGgHSxoTb6h9mJC7PTL8zGzSfXER1ilEAOjJlPwDhFVJCVAHNjdZmjzmi8oLyqzwIj4bibjGbc1cVGo4SaEuNzTbGRUNwOrIWs/2v9gAeMb0o0zgZ46oS2T85Zj0PGpcbn1c/Dp3eBwniMSbtAO8HlNBpL3LG/ab2Xzmu9YQUYw0q4E01COoycYwwvnnwsJE4it1rz2sp9mL7J58aTfTP4niQaIXohpRTBwcEEBwd7Sj1Op5OKigrKy8s973V1dYes6x46oaSkxDPNO/mEhIQQGhraqcnHLcTW/E/S/8wc1ux7ncNJXlktB0pr6RfWVHWlFIwZEEZeWS2FlfXkl9eRX14HGM8K3TJzuGfZ29/ewFdZB1EKooKsxIbaiA6xEhNiY3JqNL82evShNjiBTcNuJDLISlSwlTC7pemeVMtGVMf+Dkaeadz7cVfX1VVAXaXRQs7N2WD0rt1QDY5a4xklZ72R6FwNRhWem6POKGm1xuVo/j1upNHwwRZmVOsFxxhVgWHxkJDpWUxHDyF/9rNs3F/GtxsL+fadH8nKr/Sc+6tPSG19f0dJWp0J0YfV19d7Eo/7daR+29wCAgI8o5i6k8/hnu/pKvUOFwUVteSX11JQXkd+eS2XH5viiev6l3/ih51FFFfXH5IvzstM5MFfjQNge34FM//ZdC/HpCAiyEpEoIWIIAv3nJvGiP5GCWBF1kGy8isIs1sItQcQYg8gxGa8woMsxIUevjdvwEheWjc1SnDUNa/KU2ajxGUyG9VsrTzw6nRpSqvrKa6qp6iqntyyGgbHhniqMV9fvY/fv7m+2TqBFjPThsfymykpHDMoCpPJ1HdbnSmlTMBfgDBgtdb6BT+HJESPZ7VaiYmJISbGeEBTa01dXZ2nxFNZWUlFRQUOh+OQdR0Oh+ehUzeTyURwcHCzBBQSEtLh1m5HwxpgIjEyiMTIoFbnP36J8Qvf4XRRXFVPQUUdRVX1FFXWkRAR6FnOpWFcYjgl1Q2UVNVTUeeguMr4Iw7g3eL8/fUHeH11Tqv7SxsQzn9vOB4wEkHG3Z9gt5gbXyasASasZuN9/tTBTBtuVCsu31HKayv3eW7FOF3g0hqnS2M2KZ6bO9Gzj3P//S1Z+ZVU1TsOSZ7XnDjIk2iGxIUQHmj0gzc5NYpjh8SQkRTR1Aeej/g10SilngPOAAq01mO8ps8GHgbMwLNa63uBs4ABQDHQ+r+oEOKoKKWw2+3Y7XbPA6Vaa2pra6moqPAknsrKylZLPu6m2C0fOg0MDGyWeIKDgzvlOZ+jEWA2Edd4H6c1w/uH8u6C4z3f6x0uymoaKKupp6S6wTN8N8AJQ2MJtgVQXuOgoraByjqH55UY2ZS8ahuclNc6KK89NHEDXDB+oOfznqJqPtmc3+pyFnPz81bT4KKyzthmRJCFqGAr0cFW+oXZGZXQdN8lY2BEs1Feu4pfq86UUicClcCL7kSjlDIDWcBMjISyCrgYOBMo0Vo/pZR6U2t9/uG2LVVnQviOu+TjnXgqKytbvefTFnfVm7sE5L6n1JWln66mtaaspoHaBhc1DU7qHE7qHS7Pa0hciCfx7SmqYktuuaeEYjIpzEphNiksZlOz/ugKK+uwBpgItgYc2jVQB/W6Bza11iuUUiktJk8CdmitdwIopV7DKM3sA+obl2l7qEMhhM95l3zc1W5g3POprKykqqrKk4Bqampo7Qdta1VvYJR+3EnHnYCOtsl1d6GU8jwjdCTJ0cGe3raPJCakez/T0x3v0QzASCpuOcBkjKq0R5VSJwArWltRKTUPmAeQlJTk4zCFEC1ZrdZmvRqA0dqturraU+pxJ6LW7vtAU79u3g+amkwmgoKCPAnI/fJ39Zton+6YaFq7arTWuhq46nAraq2fBp4Go+rMB7EJITrIbDZ7mka7eVe9VVVVed5b6+EAjHs/7iTVctstE1BQUJAkoG6mOyaaHGCg1/dE4ICfYhFC+EBbVW8ul8vTe4H7VVlZSX19favbcT8b1LLxgdls9iQd93tQUJBPnv0RR9YdE80qYKhSKhXYD1wE/Nq/IQkhuoLJZDqk9ANGFzstE1B1dXWbz/w4nc5Dejxwb9+ddLyTUGBgIKa+2KlmF/F38+ZXgWlAjFIqB7hLa71YKbUA+BijefNzWutNfgxTCOFnFouFiIgIIiIiPNO01ockoOrq6sPe/2mrCs5dwvJOQu6XxSKDsB0t6RlACNGreCcgd+Jxv7e31wNvFouFwMDAZsknMDCw15aCel3zZiGE6GxKKaxWK1arlcjIyGbzGhoamiUf96u2trbN7TU0NHh6x265H7vd3iwJuT9brVa5F+RFEo0Qos+wWCzNxvZxczqd1NTUNEs+7ldbo5u6xwWqqamhuLj5UNcmk8mTdNylH/d3i8XS55KQJBohRJ9nNps93eN4czfDrq6uPiQRHa4XBO/Wc63tyzv5eL96a0lIEo0QQrTBuxl2S+5SkDsBeSeithojuNdrrUECGEnIXR3n/XLH0FOTkCQaIYT4BdoqBYFxX6dlEnJ/bm3kUzen09lmScj7npA7+Xi/d+c+4iTRCCFEJ7NYLFgsFsLCmo9Y6W4R5048LV+HKwl53xNqjdVqPST5uEtC/u4pQRKNEEJ0Ee8WcS0bJGitcTgchySf2tpaampq2uwdwa2+vt4zkF1r+3UnHe8E5P4eEBDg00QkiUYIIboBpVSbJSFofk/InXzcn2tra1vtIdvNuzTkPUS3m/veUGv3ojqDJBohhOgBDndPyN06zjsJuRNQTU3NER9UPdy9oc4giUYIIXq4w7WOAyOReCcg70RUW1t72AYKnUESjRBC9HLu3qyDgw8dSM19b8iddHxBEo0QQvRh3veGWvaa3Vl6X49wQgghuhVJNEIIIXxKEo0QQgifkkQjhBDCpyTRCCGE8ClJNEIIIXxKEo0QQgifUofrH6cnU0pVANv8HUc7xACF/g6iHSTOziVxdq6eEGdPiBFguNa6Ux+o6c0PbG7TWk/wdxBHopRaLXF2Homzc0mcnacnxAhGnJ29Tak6E0II4VOSaIQQQvhUb040T/s7gHaSODuXxNm5JM7O0xNiBB/E2WsbAwghhOgeenOJRgghRDcgiUYIIYRP9ZhEo5SarZTappTaoZRa2Mp8pZR6pHH+eqVU5pHWVUpFKaU+VUptb3yP9FecSqmBSqkvlFJblFKblFI3eq2zSCm1Xym1rvF1mj9ibJy3Wym1oTGO1V7Tu9O5HO51rtYppcqVUjc1zuvUc9nOOEcopb5XStUppW5tz7p+Op+txtmV1+bRxNk4rztdn22dz+52fV7S+P9nvVLqO6XUuCOt2+HzqbXu9i/ADGQDgwAr8DMwqsUypwEfAQo4BvjxSOsC/wAWNn5eCNznxzjjgczGz6FAlleci4Bb/X0uG+ftBmJa2W63OZetbCcPSO7sc9mBOOOAicDfvPfdDa/NtuLskmvzaOPshtdnm3F2s+vzWCCy8fOp+OBvZ08p0UwCdmitd2qt64HXgLNaLHMW8KI2/ABEKKXij7DuWcALjZ9fAM72V5xa61yt9U8AWusKYAsw4Cjj6dQYj7DdbnMuWyxzEpCttd5zlPH84ji11gVa61VAQwfW7fLz2VacXXhtHlWcR9BtzmcL3eH6/E5rXdL49QcgsR3rduh89pREMwDY5/U9h0Mv9LaWOdy6/bTWuWD8Z8L4BeKvOD2UUilABvCj1+QFjUXb546y2H+0MWrgE6XUGqXUPK9luuW5BC4CXm0xrbPOZXtj+CXr+uN8HpGPr004+ji70/XZHt3t+rwKo5bgSOt26Hz2lESjWpnWsl12W8u0Z93OcjRxGjOVCgHeAm7SWpc3Tn4CGAykA7nAg36M8TitdSZGEft6pdSJRxHL4XTGubQCZwJveM3vzHN5xBh8uG5HHfW+uuDahKOPsztdn4ffQDe7PpVS0zESzR86uu6R9JREkwMM9PqeCBxo5zKHWzffXdXS+F7gxzhRSlkw/iO/rLV+272A1jpfa+3UWruAZzCKtH6JUWvtfi8AlnnF0q3OZaNTgZ+01vnuCZ18Ltsb5y9Z1x/ns01ddG0edZzd7Po8km5zfSqlxgLPAmdprYvasW6HzmdPSTSrgKFKqdTGXwEXAe+1WOY94DfKcAxQ1likO9y67wGXN36+HHjXX3EqpRSwGNiitX7Ie4UW9x3OATb6KcZgpVRoY0zBwCyvWLrNufSafzEtqiU6+Vy2N85fsq4/zmeruvDaPNo4u9v1eSTd4vpUSiUBbwOXaa2z2rlux85ne1oudIcXRgujLIxWEP/bOG0+ML/xswIeb5y/AZhwuHUbp0cDnwPbG9+j/BUncDxGsXQ9sK7xdVrjvP80Lru+8R843k8xDsJoefIzsKm7nsvGeUFAERDeYpudei7bGWd/jF+H5UBp4+ewbnhtthpnV16bRxlnd7s+D/fv3p2uz2eBEq9/29WHW/eXnE/pgkYIIYRP9ZSqMyGEED2UJBohhBA+JYlGCCGET0miEUII4VOSaIQQQviUJBohhBA+JYlGCCGETwX4OwAh+hKl1GjgYSAJ4+G8OIweqFf5NTAhfEge2BSiiyil7MBPwAXATmArsEZrfa5fAxPCx6REI0TXORlYq7XeBJ7ee4+2d14huj25RyNE18nAKNGglEoAKrXW3/o3JCF8TxKNEF2njqbRC/+OMTyuEL2eJBohus4rwIlKqW0YPQx/r5T6l39DEsL3pDGAEEIIn5ISjRBCCJ+SRCOEEMKnJNEIIYTwKUk0QgghfEoSjRBCCJ+SRCOEEMKnJNEIIYTwqf8HxW3REFk7qMEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tinsts_ana = gwaxion.tinst_approx(mbh, alphas, chi)\n",
    "\n",
    "plot(alphas, tinsts_ana, label='Closed-form approximation', lw=3, c='gray', alpha=0.5)\n",
    "plot(alphas, tinsts_amp, label='Field amplitude', lw=2, ls='--')\n",
    "plot(alphas, tinsts_num, label='Occupation number', lw=2, ls='--')\n",
    "\n",
    "xlim(0, 0.2);\n",
    "yscale('log');\n",
    "ylabel(r'$\\tau$ (s)');\n",
    "xlabel(r'$\\alpha$');\n",
    "\n",
    "legend(loc='best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The closed-form expression underestimates the instability time, especially for higher values of $\\alpha$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dissipation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now look at the $e$-folding time of the cloud dissipation through gravitational-wave emission. As above, we can obtain a quick estimate using closed-form approximations from [arXiv:1706.06311](https://arxiv.org/abs/1706.06311), e.g. for two example values of the boson mass (i.e. two values of $\\alpha$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dissipation timescale\n",
      "---------------------\n",
      "alpha = 0.100 :\t557142.9 years\n",
      "alpha = 0.176 :\t115.7 years\n"
     ]
    }
   ],
   "source": [
    "print('Dissipation timescale\\n---------------------')\n",
    "for a in [0.1, 0.176]:\n",
    "    print('alpha = %.3f :\\t%.1f years' % (a, gwaxion.tgw_approx(mbh, a, chi) / YRSID_SI))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare this analytic approximation to numerical results, based on the estimates for the GW power from [arXiv:1706.06311](https://arxiv.org/abs/1706.06311). This can be obtained from a `BosonCloud` object through the `get_life_time()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/richardbrito/opt/anaconda3/envs/gwaxion/lib/python3.7/site-packages/gwaxion/physics.py:424: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  self.reduced_compton_wavelength = HBAR_SI / (mass*C_SI)\n",
      "/Users/richardbrito/opt/anaconda3/envs/gwaxion/lib/python3.7/site-packages/gwaxion/physics.py:1042: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  epsilon = 1./bhb_0.boson.omega\n",
      "/Users/richardbrito/opt/anaconda3/envs/gwaxion/lib/python3.7/site-packages/gwaxion/physics.py:1126: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  self.nr)\n",
      "/Users/richardbrito/opt/anaconda3/envs/gwaxion/lib/python3.7/site-packages/gwaxion/physics.py:1052: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  wR_0 = bhb_0.level_omega_re(n) * epsilon\n",
      "/Users/richardbrito/opt/anaconda3/envs/gwaxion/lib/python3.7/site-packages/gwaxion/physics.py:1053: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  wI_0 = bhb_0.level_omega_im(l, m, nr) * epsilon\n",
      "/Users/richardbrito/opt/anaconda3/envs/gwaxion/lib/python3.7/site-packages/gwaxion/physics.py:1055: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  dimless_m_accretion_rate = m_accretion_rate * C_SI**2 / gamma\n",
      "/Users/richardbrito/opt/anaconda3/envs/gwaxion/lib/python3.7/site-packages/gwaxion/physics.py:1056: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  dimless_j_accretion_rate = j_accretion_rate * epsilon / beta\n"
     ]
    }
   ],
   "source": [
    "# create an array of alphas\n",
    "alphas = np.linspace(0, 0.2, 100)\n",
    "\n",
    "# for each of those values, compute the instability timescale numerically\n",
    "tgws = []\n",
    "for a in alphas:\n",
    "    bhb = gwaxion.BlackHoleBoson.from_parameters(m_bh=mbh, chi_bh=chi, alpha=a)\n",
    "    c = bhb.cloud(1, 1, 0)\n",
    "    tgws.append(c.get_life_time())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/richardbrito/opt/anaconda3/envs/gwaxion/lib/python3.7/site-packages/gwaxion/physics.py:199: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  t = (6.5E4) * YRSID_SI * (m/(10*MSUN_SI)) * (0.1/alpha)**15 / chi\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEGCAYAAABcolNbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/20lEQVR4nO3deXiU5bn48e89k2Qm62QPIYEk7EuIAcLiDi6IWsVd0dZaqxQtbr/29GDPqXLssdXWWmu1tW6l9Vi0bpW2tLUuqCgooMi+CxKy73sySZ7fH5MZJiGBBDJLkvtzXbky87zbPS8DN8/zPosYY1BKKaV8xRLoAJRSSg1ummiUUkr5lCYapZRSPqWJRimllE9polFKKeVTIYEOwFcSExNNZmZmoMNQSqkBZePGjWXGmKT+POegTTSZmZls2LAh0GEopdSAIiIH+/uc2nSmlFLKpzTRKKWU8ilNNEoppXwqKJ/RiMhE4C4gEXjHGPNbEZkELAPKO8peDWCISvXI6XSSn59PU1NToENRqkd2u5309HRCQ0N9fi2/JRoReR74GlBijMn2Kp8P/AqwAs8aYx4yxuwAFouIBXimY9cLgV8bYz4UkZWAJhoVlPLz84mOjiYzMxMRCXQ4Sh3FGEN5eTn5+flkZWX5/Hr+bDpbDsz3LhARK/AkriQyCVjYUXNBRC4F1gDvdOz+AnCdiPwcSPBTzEr1WVNTEwkJCZpkVNASERISEvxW6/ZbojHGfABUdCmeCew1xuw3xrQALwELOvZfaYw5Dbih432JMea7wFKgrLtriMgiEdkgIhtKS0sBVzNGbW0tZWVlOJ1On3w2pbrSJKOCnT+/o4F+RpMGHPJ6nw/MEpE5wBWADVgFICKZwA+BSODn3Z3MGPM08DRAXl6eAdi+fTuVlZUATJkyhYQErQwppZQ/BbrXWXcp1RhjVhtj7jTGfMcY82RH4QFjzCJjzA3GmDW9vYCTEB7Z0MRPP2mkubm53wJXKpgVFRVx3XXXMXr0aCZNmsRFF13E7t27OXDgANnZ2cc/wQlavnw5S5Ys6XbbK6+8wsSJE5k7d67Prh9MVq5cyUMPPdQv53rsscdoaGjwvL/ooouoqqrql3P7Q6ATTT4wwut9OlDQnxeIjQpnW1kbuyvbqWvQXkBq8DPGcPnllzNnzhz27dvH9u3b+clPfkJxcXFA43ruuef4zW9+w3vvvder/VtbW30cUd+0tbX1af9LL72UpUuX9su1uyaaVatWERsb2y/n9odAJ5r1wFgRyRKRMOA6YGV/XiAi3E6cXTDA4Yq6/jy1UkHpvffeIzQ0lMWLF3vKcnNzOfPMMzvt19TUxLe+9S2mTJnC1KlTPQlg27ZtzJw5k9zcXHJyctizZw8A//d//+cp/853vuP5h/f3v/8948aN4+yzz+ajjz7qNqYHHniANWvWsHjxYv7jP/6jx2svX76cq6++mksuuYR58+axfPlyLrvsMi655BKysrJ44oknePTRR5k6dSqzZ8+moqLrY1/461//yqxZs5g6dSrnnXeeJ8EuW7aMb3zjG5xzzjmMHTuWZ55xdWhdvXo1Z511FpdffjmTJk1i8eLFtLe3AxAVFcV9993HrFmzWLt2LY8++ijZ2dlkZ2fz2GOPAfDoo49y8803A7Blyxays7NpaGjoVLu76aabuO2225g7dy6jRo3i/fff5+abb2bixIncdNNNnthvu+028vLymDx5Mvfffz8Ajz/+OAUFBcydO9dTG8zMzKSsrMxz/a4xHThwgIkTJ3LrrbcyefJk5s2bR2NjY89fGh/zZ/fmFcAcIFFE8oH7jTHPicgS4F+4ujc/b4zZ1p/XtdvtxNuFiiZDQVXgbrQamlavXu2zc8+ZM6fb8q1btzJ9+vTjHv/kk08Crn8cd+7cybx589i9ezdPPfUUd911FzfccAMtLS20tbWxY8cOXn75ZT766CNCQ0O5/fbbefHFFzn//PO5//772bhxIw6Hg7lz5zJ16tSjrnXffffx7rvv8sgjj5CXl8cvfvGLbq8NsHbtWjZv3kx8fDzLly9n69atfP755zQ1NTFmzBgefvhhPv/8c+655x7++Mc/cvfdd3e61hlnnMG6desQEZ599ll+9rOfea63efNm1q1bR319PVOnTuXiiy8G4NNPP2X79u1kZGQwf/58Xn/9da666irq6+vJzs7mgQceYOPGjfz+97/nk08+wRjDrFmzOPvss7n77ruZM2cOb7zxBg8++CC/+93viIiIOOoeVFZW8u6777Jy5UouueQSPvroI5599llmzJjBpk2byM3N5cEHHyQ+Pp62tjbOPfdcNm/ezJ133smjjz7Ke++9R2JiYqdz9hRTXFwce/bsYcWKFTzzzDNcc801vPbaa3z9618/7vfCF/yWaIwxC3soX0XHA39fsNlsxNtdj4IKa7TpTCm3NWvWcMcddwAwYcIEMjIy2L17N6eeeioPPvgg+fn5XHHFFYwdO5Z33nmHjRs3MmPGDAAaGxtJTk7mk08+Yc6cOSQluSb7vfbaaz0J40SuDXD++ecTHx/v2Xfu3LlER0cTHR2Nw+HgkksuAVydezZv3nzUufPz87n22mspLCykpaWl0ziRBQsWEB4eTnh4OHPnzuXTTz8lNjaWmTNnMmrUKAAWLlzImjVruOqqq7BarVx55ZWemC+//HIiIyMBuOKKK/jwww+ZOnUqy5cvJycnh+985zucfvrp3X7mSy65BBFhypQppKSkMGXKFAAmT57MgQMHyM3N5c9//jNPP/00ra2tFBYWsn37dnJyco55H7uL6dJLLyUrK4vc3FwApk+fzoEDB479h+JDgW468znvRFNc24IxJsARKeVbkydPZuPGjcfdr6e/C9dffz0rV64kPDycCy64gHfffRdjDN/85jfZtGkTmzZtYteuXSxbtgzovptsW1sbubm55Obmct999/X62oDnH003m83meW2xWDzvLRZLt89x7rjjDpYsWcKWLVv43e9+12msSNdY3e97Krfb7Vit1uPGvGfPHqKioigo6PkRs3fcXT9Ta2srX375JY888gjvvPMOmzdv5uKLLz7uOJdjxeR9DavVGtBnXoHu3uxzISEhJEZYgVbKG9txOp2EhYUFOiw1RPTUvOVL55xzDj/84Q955plnuPXWWwFYv349DQ0NZGRkePY766yzePHFFznnnHPYvXs3X331FePHj2f//v2MGjWKO++8k/3797N582bmzZvHggULuOeee0hOTqaiooLa2lpmzZrFXXfdRXl5OTExMbzyyiuccsopWK1WNm3a1GOMPV37s88+O+nPX11dTVpaGgB/+MMfOm178803uffee6mvr2f16tU89NBD7N69m08//ZQvv/ySjIwMXn75ZRYtWtRtzDfddBNLly7FGMMbb7zBCy+8QHV1NXfddRcffPABS5Ys4dVXX+Wqq67qc9w1NTVERkbicDgoLi7mH//4h+f7Ex0dTW1t7VFNZz3FFGwGfY0GYHxSOGekhTA21qpdnNWgJyK88cYb/Pvf/2b06NFMnjyZZcuWMXz48E773X777bS1tTFlyhSuvfZali9fjs1m4+WXXyY7O5vc3Fx27tzJjTfeyKRJk/jf//1f5s2bR05ODueffz6FhYWkpqaybNkyTj31VM477zymTZvWqxh7unZ/WLZsGVdffTVnnnnmUf8wz5w5k4svvpjZs2fzox/9yHNPTj31VJYuXUp2djZZWVlcfvnlR5132rRp3HTTTcycOZNZs2Zxyy23MHXqVO655x5uv/12xo0bx3PPPcfSpUspKSnpc9ynnHIKU6dOZfLkydx8882dmuAWLVrEhRdeeFTX8J5iCjYyWJuS8vLyjHvhs82bN3t6p2RnZx/15VOqP+3YsYOJEycGOgzVxbJly4iKiuL73/9+p/LVq1fzyCOP8Le//S1AkQVOd99VEdlojMnrz+sMiRqN9/+UtEajlFL+Neif0YAr0RTWtVPW2E7K8AbSAh2QUsrv3J0XupozZ05AnqUNJUMi0djtdh7d2ERpo2F8ej3aqKGUUv4zZJrO3F2cD+ugTaWU8quhk2jCXYmmqEaf0SillD8NmUSTYHd91OJapw7aVEopPxoSicZqtZIU6XocVdHUTktLS4AjUsq3RITvfe97nvePPPJIjw/DfWXDhg3ceeedJ3TsnDlzcA9P8DXvpRM2bdrEqlU+mxFryBoSiQZgWIyri3N5Y7t2cVaDns1m4/XXX/fM8Otvra2t5OXl8fjjj/vsGn2dtr83NNH4xpBJNKkOOwAVTUYTjRr0QkJCWLRoEb/85S+P2nbTTTfx6quvet5HRUUBroGLZ599Ntdccw3jxo1j6dKlvPjii8ycOZMpU6awb98+AEpLS7nyyiuZMWMGM2bM8CwNsGzZMhYtWsS8efO48cYbWb16NV/72tcAqKur8ywLkJOTw2uvvQZ0Py3+sWRmZvLAAw9wxhln8Morr/DWW29x6qmnMm3aNK6++mrq6lxLgSxdupRJkyaRk5PjGaDZ0+d2a2lp4b777uPll18mNzeXl19+uXc3Wx3XkOjeDDAqKZIfnx5OvF2OO1GdUv0pc+nfe9z2k8uncP2skQD86ZOv+OEbW3rc98BDF/fput/97nfJycnhBz/4Qa+P+eKLL9ixYwfx8fGMGjWKW265hU8//ZRf/epX/PrXv+axxx7jrrvu4p577uGMM87gq6++4oILLmDHjh2Aa9r6NWvWEB4e3mmJhB//+Mc4HA62bHF9Pvfy6t1Ni3+s2YrBNVxhzZo1lJWVccUVV/D2228TGRnJww8/zKOPPsqSJUt444032LlzJyLS65Uow8LCeOCBB9iwYQNPPPFEr++ZOr4hk2hiIiMYEe2qwGmNRg0FMTEx3HjjjTz++OOEh4f36pgZM2aQmpoKwOjRo5k3bx7gmpLfvTjZ22+/zfbt2z3H1NTUUFtbC7hWlezuWm+//TYvvfSS531cXBxAn6fFB9dSBADr1q1j+/btnjnBWlpaOPXUU4mJicFut3PLLbdw8cUXe2pVKnCGTKLRaWhUoPS2JnL9rJGe2k1/ufvuu5k2bRrf+ta3PGUhISGeFSSNMZ06x/RmSv729nbWrl3bbULpOsW/mzHmqKn43dPir1+/nri4OG666aZetTa4r2GM4fzzz2fFihVH7fPpp5/yzjvv8NJLL/HEE0/w7rvvHvNzK98Kymc0IjJRRJ4SkVdF5LaOsjki8mFH+Zy+ntNms/HWASePrG/is/ya/g5ZqaAUHx/PNddcw3PPPecpy8zM9KxX8+abb+J0Ovt0znnz5nVqWjrWcgA9HVNZWdnttPh9MXv2bD766CP27t0LQENDA7t376auro7q6mouuugiHnvsMU98vfnc7un4Vf/yW6IRkedFpEREtnYpny8iu0Rkr4gsBTDG7DDGLAauAdyziBqgDrAD+X29vs1mI7+una3lbXxZrs9o1NDxve99r1Pvs1tvvZX333+fmTNn8sknn/RYC+nJ448/zoYNG8jJyWHSpEk89dRTxz3mv//7v6msrCQ7O5tTTjmF995775jT4vdGUlISy5cvZ+HCheTk5DB79mx27txJbW0tX/va18jJyeHss8/2dIjozeeeO3cu27dv184A/cxvywSIyFm4EsUfjTHZHWVWYDdwPq7ksR5YaIzZLiKXAkuBJ4wxfxIRizGmXURSgEeNMTcc63reywSAq7r/veff5o29Tr42KpTHbzkPiyUoK3RqgNNlAtRAMeiWCTDGfABUdCmeCew1xuw3xrQALwELOvZfaYw5Dbih4317xzGVQLcrJInIIhHZICIbSktLO22zWCwkRYUCUK6DNpVSym8C3RkgDTjk9T4fmNXxDOYKXAllFYCIXAFcAMQC3fY9NMY8DTwNrhpN1+2uQZsNVDS6xtLY7fb++hxKKaV6EOhEI92UGWPMamB1l8LXgddP5mLDHeFAJRVNhqamJhwOx8mcTqkeddfLSqlg4s85HwP9kCIfGOH1Ph0o8NXF0hNcD/8qmwyNOmhT+Yjdbqe8vFwnb1VByxhDeXm531p1Al2jWQ+MFZEs4DBwHXC9ry7miAxneoqViBChtkETjfKN9PR08vPz6fqcUKlgYrfbSU9P98u1/JZoRGQFMAdIFJF84H5jzHMisgT4F2AFnjfGbPNVDDabjTumdmTwVu0MoHwjNDSUrKysQIehVNDwW6IxxizsoXwVHQ/8fc17JHNDQ4M/LqmUUkNeoJ/R+FVkZCQGOFzXzraiOs+UGkoppXxnSCUai8XCvrpQ/mtNIyt2tHimFFdKKeU7QyrRAExJjwXgYG071TU6p5FSSvnakEs0aYkOEuxCSxvsPNx1ogKllFL9bcglmujoaDJiXB97W4HO4qyUUr425BJNZGSkJ9HsLm30ybrjSimljhhyiSYkJIRxSa5uzgdq2qivrw9wREopNbgNuUQDMCXNNcdZQV27LnKklFI+FugpaAJiZHIsPz69jNRI0S7OSinlY0My0URHRzMi2lWZ00SjlFK+NSSbzqKiojyv6+vraW9vP8beSimlTsaQrNGEhoZS3Gzl+U11RIcJ06Y1dEo+Siml+s+QrNEAJDmi2VXZzq6KNu0QoJRSPjRkE01WioPIUKh1wv6iykCHo5RSg9aQTTQxMTFkdgzc3HK4KrDBKKXUIDZkE01UVBQZMVYAdhbV67K7SinlI0M20YSFhTEuIQyAXRWtOkOAUkr5SFAmGhGZKCJPicirInJbR9koEXlORF7tp2swe1Q8AuypbKegpLw/TquUUqoLvyUaEXleREpEZGuX8vkisktE9orIUgBjzA5jzGLgGiCvo2y/Mebb/RlTRmoSl4wO5VvZYVRXaYcApZTyBX/WaJYD870LRMQKPAlcCEwCForIpI5tlwJrgHd8FVBcXBxXjA3j9LRQGutqdOCmUkr5gN8SjTHmA6DrSmMzgb0dtZUW4CVgQcf+K40xpwE39PYaIrJIRDaIyIbS0tLj7m+324mIiACgvb2d6urq3l5KKaVULwX6GU0acMjrfT6QJiJzRORxEfkdsApARBJE5Clgqojc293JjDFPG2PyjDF5SUlJvQogLi6OzaWt/GFbMzsOFp/cp1FKKXWUQE9BI92UGWPMamB1l8JyYHF/BxAfH8/7+fvYWNzG2N0lnJY7ob8voZRSQ1qgazT5wAiv9+lAgT8DiI2NZUqiK99uPNxAS0uLPy+vlFKDXqATzXpgrIhkiUgYcB2w0p8BWK1WTs2KBWB7eRul5V0fIymllDoZ/uzevAJYC4wXkXwR+bYxphVYAvwL2AH82RizzV8xuU0YkURKhNDYCmt3F/r78kopNaj57RmNMWZhD+Wr6HjgHyjx8fFkJ1op/qqVj/dVcOVZBpHuHh8ppZTqq0A3nQWFqKgoclNsAGwqbtHpaJRSqh9posE1Hc1pYxIZHiWckmSlvFyno1FKqf4S6O7NQSMtJZGfnFEGQGlpKRkZGQGOSCmlBget0XRITEzEYnHdjrq6Om0+U0qpfqKJpkNISAgJCQmUN7azan8Ln+05dPyDlFJKHZcmGi/Jycn8Za+TP+928vpnBboYmlJK9QNNNF4SEhI4Y4Sr99maQ03U1NQEOCKllBr4NNF4sVgsnDl+GHE2obTR8N6Wg4EOSSmlBjxNNF2kDkthVqqrM97fthZr85lSSp0kTTRdxMbGcnaGHYB1h1soKdMxNUopdTI00XQhIuSNSSUtSqhzwt8/OxDokJRSakDTRNONYcOGcWZaKKckWbE01+J0OgMdklJKDVg6M0A3oqKiuDI71jNos6ioiBEjRhznKKWUUt3RGk03RIT09HTP+8OHD2unAKWUOkGaaHqQnJxMaGgouyraeGRtFXu+0nVqlFLqRGii6YHVaiU1NZVVXzrZUNzG7z/cG+iQlFJqQArKRCMil4nIMyLypojM6yg7U0SeEpFnReRjf8SRlpbGBZmhAKzaU0dFlc4UoJRSfeXPpZyfF5ESEdnapXy+iOwSkb0ishTAGPMXY8ytwE3AtR1lHxpjFgN/A/7gj5htNhtnjU8hPUqobjas+GiXPy6rlFKDij9rNMuB+d4FImIFngQuBCYBC0Vkktcu/92x3dv1wArfhdnZiBEjmNdRq1nxeRlNTc3+urRSSg0Kfks0xpgPgIouxTOBvcaY/caYFuAlYIG4PAz8wxjzmXtnERkJVBtjum3DEpFFIrJBRDaUlpb2S9wxMTHMGxdHvF3Ir2vnhfe398t5lVJqqAj0M5o0wHvhl/yOsjuA84CrRGSx1/ZvA7/v6WTGmKeNMXnGmLykpKR+CVBEGJ2VwYLRrlrNCxuLaG7WWo1SSvVWoAdsSjdlxhjzOPB4Nxvu931IR0tMTOSC8bFUNlcyd0Qohw4dYsyYMYEIRSmlBpxA12jyAe8h9+lAQYBi6ZGIMGZUFpeNCcNhEwoKCrRWo5RSvRToRLMeGCsiWSISBlwHrAxwTN1KSEggOjoagJbWNj78Yk+AI1JKqYHBn92bVwBrgfEiki8i3zbGtAJLgH8BO4A/G2O2+SumvhARMjMzqWxq57/WNPL//vYV5dV1gQ5LKaWCnt+e0RhjFvZQvgpY5a84TkZ8fDwjkxxEhTZT3NDOQ3/9gp9//fRAh6WUUkEt0E1nA4qIMGrUKG6YFIYAr2+r4vN9OgeaUkodiyaaPoqLi2PGqGTOHhFCm4Ef/WUr7e3tgQ5LKaWCliaaEzBmzBiuGW8nKhS2lrbwfx/sCHRISikVtDTRnIDw8HAmjh7J1ePDAHj0vYNU1TcFOCqllApOmmhO0MiRIzlvVCRTk61cPS6U0oJDxz9IKaWGoEDPDDBgWa1Wxo4Zw10trrnPDh8+TEpKCjExMQGOTCmlgkufazQiEtkx6/KQl5SURFxcnOf9vz/dSlW9zhiglFLejptoRMQiIteLyN9FpATYCRSKyDYR+bmIjPV9mMFJRBg3bhxWq5V1Ba385zuV/OCl9YEOSymlgkpvajTvAaOBe4FhxpgRxphk4ExgHfCQiHzdhzEGtfDwcEaNGkVGjAUReGtPNW+s3x/osJRSKmj0JtGcZ4z5sTFmszHGM2DEGFNhjHnNGHMl8LLvQgx+w4cPZ2J6PNeMc/VCu++vO8mvqA9wVEopFRyOm2iMMU4AEblaRKI7Xv9IRF4XkWne+wxVIsKECROYN8rGlEQrtS2GW36/jpZWHciplFJ96QzwI2NMrYicAcwD/gD81jdhDTx2u52xY8awKMdGvF3YWdrEj1777PgHKqXUINeXRNPW8fti4LfGmDeBsP4PaeBKTU1lVFoy3821YRVYuaWYgyVVgQ5LKaUCqi/jaA6LyO9wLbH8sIjY0AGfnYgI48ePp7a2lltzDJkxFsoO7WNE4lQsFr1VSqmhqS//+l2Da92Y+caYKiAe+A9fBDWQhYSEMGnSJE4dHsqwSAu1tbXs3buXtnYT6NCUUiogejOORgCMMQ3GmNeNMXs63hcaY97y3ke5xMTEMHr0aM/7P649yOW/fp8mZ9sxjlJKqcGpV+NoROQOERnpXSgiYSJyjoj8AfhmfwYlIpeJyDMi8qaIzOsomygiT4nIqyJyW39ezxfS0tJISkqisdXwjy+dbC6s564/baBdazZKqSGmN4lmPq6OACtEpEBEtovIfmAPsBD4pTFm+fFOIiLPi0iJiGztUj5fRHaJyF4RWQpgjPmLMeZW4Cbg2o6yHcaYxbia8PJ6/xEDw93lOSk2mrun27Fb4V87yvjp34NypWqllPKZ3oyjaTLG/MYYczqQAZwLTDPGZBhjbjXGbOrltZbjSloeHXOmPQlcCEwCForIJK9d/rtju3v/S4E1wDu9vGZAWa1WsrOzGZ1g5/ZcGxaBZz46yFOr9wQ6NKWU8ps+dYUyxjg7ns1U9fVCxpgPgIouxTOBvcaY/caYFuAlYIG4PAz8wxjzmdc5VhpjTgNu6Ov1A8VutzN58mROSQ7llik2AB76525WfHIwwJEppZR/BLrPbRrgvZBLfkfZHbi6UV8lIosBRGSOiDze0cV6VXcnE5FFIrJBRDaUlpb6OPTeczgcTJgwgdOGh/D1ia6hR4/8awf1zUN6QgWl1BAR6PVouuutZowxjwOPdylcDaw+1smMMU8DTwPk5eUF1VP3lJQUmpubOQ/XhJtTEq0cPvglY8eORTvtKaUGs14nGhH5uKPZqj/lAyO83qcDBf18jaAxYsQIWlpaOI98AAoKCggJCaE9KokxydEBjk4ppXyjL01n9q4FInLmSV5/PTBWRLJEJAy4Dlh5kucMWiLC6NGjSUlJ8ZQ988E+5v3yA/7y+eEARqaUUr7Tl6az8SLyBrAN2AoUA8/iWqvmuERkBTAHSBSRfOB+Y8xzIrIE14wDVuB5Y8yg7v/rnqamtbWV8vJyGlsN7QbueXkTjc42Fs4cefyTKKXUANKXRPMl8BMgG5gODAf+p7cHG2MW9lC+ih4e7g9WFouFyZMns3XrVi4bU4FV4LU9Tu59fQtltc0sOWeMPrdRSg0afUk0LcaY9biau9RJslgsZGdns2XLFi4ZXUlkqPDC9hZ+8e/dlNQ2s+zSyVgtmmyUUgNfX57RnO2zKIYod7KJi4vjnJGh3J5rI0TghXUHefifOwMdnlJK9YteJxpjTK0vAxmqrFYrU6ZMISEhgRnDQvhenp2UCOHURCft7bpCp1Jq4Av0gE3FkWc2ycnJTEyw8pMzwpGGCrZs2YLT6eRQRUOgQ1RKqROmiSZIWCwWJk6cyPDhwz3PZiorK/mfP3/MvF++z5ubtPuzUmpg0kQTRESEsWPHkpmZCYAxhoMVTTQ627nrpU38z1+34WzT5jSl1MCiiSbIiAiZmZlMmDABi8XCzdlhfGNSGFaB3390gOufWUdJTVOgw1RKqV7TRBOkhg0bxpQpUwgJCeHckaHcO9NOrE1Yf6CSC3/1Iat3lQQ6RKWU6hVNNEEsPj6eadOmYbfbGRNn5X9OC2dSgoXy+hZ+/c4ejAmqeUOVUqpbmmiCXGRkJNOnTyc2NhaHTfh+np1rxoVy8yQLDQ2u3miacJRSwUwTzQAQGhpKTk4OaWlpWES4aFQYkTSzceNGioqK+M4LG/nj2gOacJRSQUkTzQBhsVgYO3YsEydOxGJx/bG1t7fz6pptvLW9mPve3MaNz39KQVVjgCNVSqnONNEMMCkpKUyfPp2IiAgAshOtfDfXRlSY8OGeMub98gP+b91B2tu1dqOUCg6aaAagyMhIpk2bRnJyMgAzhoXwv6fZmZ5ipa65lf/+y1YWPrOOL8vqAxypUkppohmwQkJCmDhxIhMmTMBqtRJrt7Ak18btuTYcNguffFnBqi2FgQ5TKaX6tEyACjIiwrBhw4iJiWH79u3U1dUxc1gIk+Kt/PNgK5eMDccYg4hQVtdMYpQt0CErpYYgrdEMAhEREUybNo2MjAwAosKEq8aGsn/vHrZs2UJheQ3nP/o+t7+4kfxKnaBTKeVfQZloRGSUiDwnIq96lU0SkT+LyG9F5KpAxheMLBYLWVlZTJs2jfDwcE95RUUFr763gcaWNlZtKeLcX7zPY2/vprGlLYDRKqWGEr8lGhF5XkRKRGRrl/L5IrJLRPaKyFIAY8x+Y8y3u5ziQuDXxpjbgBv9FPaAExMTQ15eHunp6Z6yKYkWHjzdxhkj7DS3tvPY23s45xereW1jvvZOU0r5nD9rNMuB+d4FImIFnsSVRCYBC0VkUg/HvwBcJyI/BxJ8GOeAZ7VaGTNmDFOnTvV0g04It3DLZCv3zgpnTIKNwuomvvfKF/zna5sDHK1SarDzW6IxxnwAVHQpngns7ajBtAAvAQt6OL7EGPNdYClQ1t0+IrJIRDaIyIbS0tJ+jH5gcjgc5OXlkZGRgYhrjZvxcRZ+mGfltqmRJEeFctX0IzWfVl2CQCnlA4F+RpMGHPJ6nw+kiUiCiDwFTBWRewFEJFNEngb+CPy8u5MZY542xuQZY/KSkpJ8HfuA4H52M2PGDGJjY11lIsxKgZ+cFoqt5hD19a7xNne9vInbX9zInmJdtVsp1X8C3b1ZuikzxphyYHGXwgPAIn8ENRhFRERwyimnUFxczL59+3A6nYRYhKqqKjZs2IDNkcQ7O4ppcrbzj61FXJIznDvPHcuY5KhAh66UGuACXaPJB0Z4vU8HCgIUy6DnHncza9Ys0tPTPc1pxhiaqkp4+KwILsuOJ8QirPyigHm/fJ+7XvpcazhKqZMS6ESzHhgrIlkiEgZcB6wMcEyDXkhICGPGjCEvL4+4uDhPuSPUcFl6M78818FlUxKxiPDmpgIu/vUaqhpaAhixUmog81vTmYisAOYAiSKSD9xvjHlORJYA/wKswPPGmG3+immoi4yMJCcnh4qKCvbt2+dZ3ybK4uSyNCcXpCfwbpGVcLud2IgwANraDRsOVDAzK95TI1JKqWPxW6IxxizsoXwVsMpfcajORISEhATi4uIoKCjg4MGDOJ1OAMJNIxengMNhpaqqitjYWN7aVsRtL35GTrqDW88cxYXZwwixBrpirJQKZjJYF8vKy8szGzZsCHQYA05rayuHDh3i0KFDtLd37u4cFxfHltpwfvneQcrrXU1pabHhfOv0TK6dMYJoe2ggQlZK9SMR2WiMyevXc2qiUd1pbm7m4MGDFBYWHrVyZ0S0gy9qwlmxsYj9HUsRRIZZuX3uGL47d0wgwlVK9RNfJJpAd29WQcpmszFu3DhGjhzJgQMHKCoq8mxrqK1mrFTzi/NiOegcwUufl7BufwW2kCNNaA0trYRZLdqsppTSGo3qnYaGBg4ePEhxcfFR22JiYmiJSGJiRgqOcFengV+8tYvXNuazcOZIrp0xguQYu79DVkqdAG066wNNNL7hTjglJSVHNalFRkYyYsQIkpOTufbpdaw/UAlAiEU4b2IKC2eN5MwxiVgs2ltNqWCliaYPNNH4VmNjI1999RVFRUVHJRybzcbwtDT214exYn0+7+wsoa1jlui02HB+9LWJzM9ODUTYSqnj8EWi0QZ0dULCw8MZP348s2fPJj09HavV6tnW3NzMl/v3E1K2hx+cFsu7d5/Gf1wwnvS4cA5XNRLj1TvtcFUjDS2tgfgISik/0RqN6hdOp5OCggLy8/M943C8JSQkMDwtjd2V7czKSvA0n317+XrW7S/noimpXDEtnVlZ8dq0plQAaa8zFbRCQ0PJyMhgxIgRFBcXc+jQIc9MAwDl5eWUl5cTGRlJob2ZlJQUEAu1Ta3Ut7TxysZ8XtmYT6rDzqW5w7ksN42JqTEB/ERKqf6iNRrlE8YYKisryc/Pp6Ki6zJErsXZhg0bxvDhwylqMLzx2WH+sukw+ZWNnn1+dlUO1+SNOOpYpZTvaI1GDRgiQnx8PPHx8TQ0NHD48GGKiopoa2sDoK2tjcOHD3P48GEcDgffyB3O3eeO4fP8at74/DD/2lrEnHFH1hR6Yd1B6ptbuXhKKiPiIwL1sZRSJ0BrNMpvWltbKSoqoqCgoFOzmltISAgpKSmkpqYSHhGJ1XJkGYMzf/aep7ZzSrqD+dmpXJg9jMzESL9+BqUGO+3e3AeaaIKXMYaqqioKCgooKys7qns0QHR0NKmpqSQnJyMWK//eXszftxTyzo5iGlraPPtNTI3hBxeMZ+6EZH9+BKUGLW06U4OCiBAXF0dcXBwtLS0UFhZSWFhIU1OTZ5/a2lpqa2vZu3cviYmJzBw+jAsm59LkbOf93aX8Y2sh7+woYUdhDaFe09x8/lUlLa3tTM+I0+lvlAoSWqNRQcFdyyksLKS0tLTbWk5YWBgpKSkMGzaMyMhImpxtfLyvjDPHJnmSzS1/WM/bO0qIiwhl7oRkzpuYwpljE3VmaaV6SZvO+kATzcDldDopLi6mqKiIurq6bveJiooiJSWF5ORkbDabp/wXb+3ir18UcKD8yDOgUKswKyuBb5yawQWTh/k8fqUGsiGTaERkFPBfgMMYc1VH2Rzgx8A24CVjzOpjnUMTzeBQW1tLcXExxcXF3Q4EBYiNjSUlJYXExERCQ0MxxrCvtI63d5Twzo5iNh6spN3A/ZdM4lunZwGwv7SOg+UNzB6VQHiYtdvzKjUUDehnNCLyPPA1oMQYk+1VPh/4Fa6lnJ81xjxkjNkPfFtEXvU6hQHqADuQ76+4VWBFR0cTHR3NqFGjqKyspKioiPLy8k6LslVVVVFVVcXu3buJj48nOTmZzIQEFp89msVnj6aivoX3d5cwKyvBc8wrG/P57ep9hIVYmJUVz5ljEzlrXBLjU6J1iWql+pnfajQichauRPFHd6IRESuwGzgfV/JYDyw0xmzv2P6qV43GYoxpF5EU4FFjzA3Hup7WaAav1tZWSktLKS4upqqqqtt9LBYL8fHxJCUlkZCQQEhI5/9T/emTr3h5/VdsPlyN91+B5GgbV05P5z/nT/DhJ1AqeA3oGo0x5gMRyexSPBPY21GDQUReAhYA27s53v1f2ErA1nW7GjpCQkJITU0lNTWV5uZmT9Kpra317NPe3k5ZWRllZWWepJOYmEhCQgKhoaFcP2sk188aSXldM2v2lvHB7jI+3FNKSW0zVQ1HmuhKapv4zXv7OH1MIrNGxXeaEFQp1TuB7t6cBhzyep8PzBKRBOBBYKqI3GuM+amIXAFcAMQCT3R3MhFZBCwCGDlypC/jVkHCZrORnp5Oeno6jY2NlJSUUFpa2qkTgXfSERFiY2M9NZ2EKBsLctNYkJuGMYZdxbWdukt/vLec5R8fYPnHB7AITElzMHt0AqeOSmBGZjyRtkD/FVIq+Pm1M0BHjeZvXk1nVwMXGGNu6Xj/DWCmMeaOk72WNp0NbQ0NDZSWlh6VdLqKiYkhMTGRxMREIiKOntpmb0ktK78oZO2+Mj7/qorW9iN/X+yhFjbdNw97qKszQXNrG7YQ7VigBrYB3XTWg3zAe9bEdKAgQLGoQSQiIoKMjAwyMjJoaGigrKyM0tLSTs1rADU1NdTU1LB//34iIiJctZyEBBwOByLCmORo/t/50XD+OOqbW9lwsJK1+8pZu6+MsBCLJ8m0txtOf+hdUh3hzMqKZ0ZWPHkZcSREaSuvUoGu0YTg6gxwLnAYV2eA640x2072WlqjUd1pamryNKNVV1d3OzAUXM+B3EknLi6O0NCjn820tRvPfGwHy+s59xfvd6rxAIxOiiQvI55bzsxibEp0/38gpfrZgB5HIyIrgDlAIlAM3G+MeU5ELgIew9W9+XljzIP9cT1NNOp4nE4n5eXllJWVUVFR0anLtDcRISYmhoSEBOLj44mMjOy2C3R9cyuff1XFJ1+Ws+FAJZ8fqqTJ6TrnyiWnk5MeC8BfPj9MUU0T0zPimJLm8NSKlAoGAzrR+JsmGtUXbW1tVFVVUVZWRnl5OS0tLT3ua7PZPEsgxMXFHdV12s3Z1s7Ww9VsPFjJN0/L9HQyuO7ptazb71qjJ9QqTEqNYerIOKaOjCUvM5602PD+/4BK9ZImmj7QRKNOlDGGuro6z6qgXZ/reHPXdtxJJzr6+AM+//pFAWv3l/PZwUp2Fdd2GsdzTV46P7vqFAAq61vYWlBNTlosjgjtVq38QxNNH2iiUf2lpaWFiooKysvLqayspLW1tcd9Q0JCPDNTx8fHY7fbj3nu2iYnm/Or+exgJZ8fqmJB7nAW5KYB8PfNhXz3T58BkJUYSU66g5z0WHJHOJiU6tCpc5RPaKLpA000yheMMdTU1HiSzrFqOwDh4eHExcURGxvbY6eCnry1rYin3t/H1oIaWlo7Pz8KC7GwZdk8T3fqvSW1pMVGaPJRJ20wdm9WakARERwOBw6HA3DVdiorK6moqKCysvKoZzuNjY00NjZSUODqtR8VFeVJPA6Ho8fnOwDzJg9j3uRhONva2VVUyxf5VWw+VM3mw9WEWMSTZIwxXPXUWmoanYxNjmZyWgzZwx1kpzmYmBqtSySogNMajVL9xBhDfX09lZWVVFZWUlVV1WNPNnAlrejoaGJjYz2Jx2rtXY3Eu2t1dYOTa59ey56SOtraj/77/MjVp3DV9HTANaVOezukxNh08lDVLa3RKBXERISoqCiioqIYMWIE7e3t1NTUeBJPbW1tp3E77ma4mpoavvrqq06Jx11r6qnG404yAI6IUP5591k0trSxs6iGrQU1bM2vZmtBNXuK68hMODLjwR8/PsgT7+0lPjKMCcOimZgaw8TUGCYMi2ZsSpTObKB8QhONUj5isVg8tZWsrCxaW1uprq721Ha6To3jnXjcoqOjPUnH4XAQFhbW4/XCw6wd3aTjPGUtre145STajcERHkpFfQsf7yvn433lnm0TU2P4x11nemJ5b1cJY5OjSY8L19qPOinadKZUgDidTs9aOlVVVdTX1x/3mIiIiE6Jx2639zkJGGMoqG5iR0ENO4tq2FFYy46iGnLTY3n02lwACqsbOfWn7wIQZQthbEoUE4ZFMy7F9XPKiFiidELRQUl7nfWBJho10DidTqqrqz2J51iTgbqFhYXhcDiIiYnB4XAQFRWFxWI57nHdMcZ4ktbekjruX7mVXUW1lNUdPXj1tdtOZXpGPAD/3FpEYXUjY5OjGZMcpc9/Bjh9RqPUIBYaGuqZSRrwNLW5k0/XZzzg6vXmnqUaXM110dHRnsQTExNzzOY2b97JYUxyFC/eMhuAsrpmdhfVsrOolj0ltewurmNM8pF5217ZcIh3dpZ43kfbQhidHMWY5ChOH5PA5VPTT+yGqEFDE41SQcp7Yk9wTZNTW1vrST7V1dW0tbV1Oqa9vd2z7dAh11JPdrudmJgYz09faz2JUTYSx9g4bUxit9svmpJKcoyNvSV17C6uo7rRyaZDVWw6VEVrW7sn0RyuauSm5z9ldFIUo5IiGdXxe3RilM58MMhpolFqgLBarZ7OBXCkO3V1dTU1NTVUV1fT1NR01HFNTU00NTVRUuKqdbh7t7lrPjExMSf0rMftyunpXNnRfdoYQ1ldC3tL6thXWsfI+CM93vYU17KnpI49JUc3CSZGhfHSotmemtKOwhpEIDMhUicdHQT0GY1Sg0hLS4sn8dTU1FBbW3vMsTxuoaGhnZJPdHR0r5vceqvJ2cae4jr2l9Wxr7SefaV17C+t58uyOpqc7WxeNs+zVPZNv/+U1btczYGpDjuZCZFkJkaSmRDB1JFxzMyK79fY1BH6jEYpdUxhYWEkJSWRlJQEuJrS6urqOiWexsbGo45zOp1UVFRQUVHhKbPZbJ0ST1RUVJ+m0OnKHmplSrqDKemOTuXt7Ybi2iZPkgFIdYSTlRjJoYoGCqubKKxuYu1+V1fsq6anexLNwfJ6vv/KF4yMjyQjIYKMhAhGxEcwMj6ChMgw7ZQQJDTRKDWIWSwWT/OYW0tLC7W1tZ7EU1NT0+1Eoc3NzTQ3N1NWVuYps9vtnppPfyQfV4xCqqPz0gg/vWIKAK1t7eRXNnKgvJ6D5Q18WVbP9Iwj44T2l9az/kAl6w9UHnXeiDArf73jDEYnRQHw0d4y6ptbGRHvSkbaPdt/9E4rNcSEhYV16mRgjKGpqalTraeurq7bJjf38x53Lzc4knzcsyL0Z7NbiNXiajJLjOx2+7SMOP50yywOVjRwsLyBQxUNfNXxU93oJCXmyOzZT72/jw/3HEmajvBQ0mLDSY8L56xxSXx9dgbgSm7VjU7itUbUbzTRKDXEiQjh4eGEh4eTkpICHOloUFtb6/mpq6vrdunr7pJPWFiYJ+m4E9DJdDjoiSM8lNPGJHJaN9uqG52dai15GfFYLcKhigbyKxupbnRS3ehke2ENcRFHEuP+snrm/fIDbCEWhseGMzzWTqojnOEOO6mx4VwweRjxkf37/GqwC8pEIyKjgP8CHMaYq3oqU0r5hve8bampqYDreY938qmrq+sx+bjX8PF+5mO1Wj3ndP9ERkae8ADT43GEd27Su+u8sZ7X7t5xh6saya9sINVxpOZT1eDEER5KdaOTL8vq+bKs84wNeRlxnkTzo79sZc3eMobF2EmJsZESYye543VmQiTZaZ2fRw1Vfks0IvI88DWgxBiT7VU+H/gVYAWeNcY8ZIzZD3xbRF5179ddmVLKf9yDQaOjjwzWdCefurq6Tsmnu2a3trY2zxgfN3dtyjvxREVFERbm22YrESEp2kZStI3cEbGdts3MiueL++dR2+SksLqJgqrGTr9TvZbadieirskI4OxxSfzh5pmAa7XUy3/zEUnRNte4pCgbCVFhJETZSIgMY2ZWPIlRNqDzDA3+1t3s3/3BnzWa5cATwB/dBSJiBZ4EzgfygfUistIYs92PcSmlTpB38nHXfIwxNDQ0eJKO+8fpdB51vHvfhoYGzzgfcA1W9U48kZGRREZG9noZhf4QbQ8l2h7KuJToHvf57denUVTdRFFNE8U1zRTXNHl+ctJjPfuV1jVzoLyBA+UN3Z7nxVtmkTjGlWge/ucuXvzkILERoTjCj/xE20IZmRDBd+eO8Ry3akshIRbBHmrFFmLBFmolxCKEWIVhMXZiO5oEK+tbyK9spKGllQZnG/XNrZTVNlNS24wB/nP+BIBOE7D2J78lGmPMByKS2aV4JrC3o7aCiLwELABOKNGIyCJgEcDIkSNPPFil1AkTEU9i8H7m09LSclTy6a6rNbim33HP+ebNbrd3SjyRkZGEh4f7rPnteNzJaOwxkhG4Bp6+/f/Opqyu2fVT20xZXQvl9S1U1DeT5lVLqmpoobapldqmVg7R+f5MGBbtSTTGGO5Y8XmPtZAfL5jMN07NBOCf24q49/Ut3e5nEbjznLGEh1l9VpMK9DOaNOCQ1/t8YJaIJAAPAlNF5F5jzE+7K+t6MmPM08DT4Bqw6fvwlVK9ISLYbDZsNpuntxu4mtPcTW91dXWe112n1nFzdzzw7nItIkRERHgSj/t1eHjwLG8QFmJhTMf8b8fzk8unsPTCCVQ3OqlqcHo6LdQ2tRJpO1Kja2s3zM8eRrOzjSZnO03ONppb22lrN7S1G+K8OizERYSSnRZDRGgI9jArkWFWEqLCSI62kxRto83HA/cDnWi6+xYYY0w5sLhL4VFlSqmBzWq1HjXOxxhDc3OzJ/G4k09jY2O3HQ/cPeS6LrNgsViIiIg4KgkFUwLqjsUixEaEERsRRkZCz/uFWC08ef20Xp1zfnYq87NT+ynCvgt0oskHRni9TwcKAhSLUioIiAh2ux273e6ZyRpcHQ8aGho8ScWdgJqbm7s9j3tWhK7LLbhrQO4E5H4dERERsCa4wS7QiWY9MFZEsoDDwHXA9YENSSkVjCwWi6d3mrfW1tajElB9fT0tLUevowOda0DeY38AwsPDOyUe98/Jzn4w1Pmze/MKYA6QKCL5wP3GmOdEZAnwL1zdm583xmzzV0xKqYEvJCTkqOY3cM3f5k5A3r97qgEBNDY20tjYSHl5eafysLAwT7ObdwLyxSDUwUhnb1ZKDSmtra2epOOdgLpbYuF43OOA3InH+/VArQXp7M1KKXWSQkJCcDgcOBydR+23tbXR2NjoSUDuJNTY2NjjUgve44C6u4538vH+7c/xQMFAE41SStF5ihxv7klH3QnFOxn19BwIXDUn90SlXYWFhXVKPN4/gzEJaaJRSqlj8J501HsMEBzpiOCdfNyvj7XgXEtLi2eRuq5sNttRySc8PBy73U5IyMD8J3tgRq2UUkGgp44I7pkQ3InHOxE1NTV1Ox7Izb0OUNdZEcBVE7Lb7d0modDQ0KDtmKCJRiml+pn3TAhxcXGdtrmb4tzJx52IGhsbj5uE3DWh7prjrFarJ+l4/w4PD8dmswV0jJAmGqWU8iPvprj4+PhO29rb22lubu6UfNwJ6XhJqK2trdsBqm7uQbBdE5E/akOaaJRSKkhYLBZPEurKPTVP1yTkrh31ND+cm3ueuJ6u6048vqCJRimlBgDvqXm6a45zOp2exONOPu73xxqkCkem9+mum3Z/0ESjlFIDnIgQFhZGWFjYUeODwJVIuktAva0NnSxNNEopNch5z2TdlTGG1tbWYzatnSxNNEopNYSJCKGhoYSGhnZaprs/6ZzYSimlfEoTjVJKKZ/SRKOUUsqnNNEopZTyKU00SimlfEoTjVJKKZ/SRKOUUsqnBu1SziJSC+wKdBy9kAiUBTqIXtA4+5fG2b8GQpwDIUaA8caYfh1QM5gHbO7q73WvfUFENmic/Ufj7F8aZ/8ZCDGCK87+Pqc2nSmllPIpTTRKKaV8ajAnmqcDHUAvaZz9S+PsXxpn/xkIMYIP4hy0nQGUUkoFh8Fco1FKKRUENNEopZTyqQGTaERkvojsEpG9IrK0m+0iIo93bN8sItOOd6yIxIvIv0VkT8fvuK7n9VecIjJCRN4TkR0isk1E7vI6ZpmIHBaRTR0/FwUixo5tB0RkS0ccG7zKg+lejve6V5tEpEZE7u7Y1q/3spdxThCRtSLSLCLf782xAbqf3cbpz+/mycTZsS2Yvp893c9g+37e0PH3Z7OIfCwipxzv2D7fT2NM0P8AVmAfMAoIA74AJnXZ5yLgH4AAs4FPjncs8DNgacfrpcDDAYwzFZjW8Toa2O0V5zLg+4G+lx3bDgCJ3Zw3aO5lN+cpAjL6+172Ic5kYAbwoPe1g/C72VOcfvlunmycQfj97DHOIPt+ngbEdby+EB/82zlQajQzgb3GmP3GmBbgJWBBl30WAH80LuuAWBFJPc6xC4A/dLz+A3BZoOI0xhQaYz4DMMbUAjuAtJOMp19jPM55g+ZedtnnXGCfMebgScZzwnEaY0qMMesBZx+O9fv97ClOP343TyrO4wia+9lFMHw/PzbGVHa8XQek9+LYPt3PgZJo0oBDXu/zOfqL3tM+xzo2xRhTCK6/TLj+BxKoOD1EJBOYCnziVbyko2r7/ElW+082RgO8JSIbRWSR1z5BeS+B64AVXcr66172NoYTOTYQ9/O4fPzdhJOPM5i+n70RbN/Pb+NqJTjesX26nwMl0Ug3ZV37Zfe0T2+O7S8nE6dro0gU8BpwtzGmpqP4t8BoIBcoBH4RwBhPN8ZMw1XF/q6InHUSsRxLf9zLMOBS4BWv7f15L48bgw+P7auTvpYfvptw8nEG0/fz2CcIsu+niMzFlWj+s6/HHs9ASTT5wAiv9+lAQS/3Odaxxe6mlo7fJQGMExEJxfUX+UVjzOvuHYwxxcaYNmNMO/AMriptQGI0xrh/lwBveMUSVPeyw4XAZ8aYYndBP9/L3sZ5IscG4n72yE/fzZOOM8i+n8cTNN9PEckBngUWGGPKe3Fsn+7nQEk064GxIpLV8b+A64CVXfZZCdwoLrOB6o4q3bGOXQl8s+P1N4E3AxWniAjwHLDDGPOo9wFdnjtcDmwNUIyRIhLdEVMkMM8rlqC5l17bF9KlWaKf72Vv4zyRYwNxP7vlx+/mycYZbN/P4wmK76eIjAReB75hjNndy2P7dj9703MhGH5w9TDajasXxH91lC0GFne8FuDJju1bgLxjHdtRngC8A+zp+B0fqDiBM3BVSzcDmzp+LurY9kLHvps7/oBTAxTjKFw9T74AtgXrvezYFgGUA44u5+zXe9nLOIfh+t9hDVDV8TomCL+b3cbpz+/mScYZbN/PY/25B9P381mg0uvPdsOxjj2R+6lT0CillPKpgdJ0ppRSaoDSRKOUUsqnNNEopZTyKU00SimlfEoTjVJKKZ/SRKOUUsqnNNEopZTyqZBAB6DUUCIik4FfASNxDc5LxjUD9fqABqaUD+mATaX8RETswGfA1cB+YCew0RhzRUADU8rHtEajlP+cB3xujNkGntl7T3Z2XqWCnj6jUcp/puKq0SAiw4E6Y8xHgQ1JKd/TRKOU/zRzZPXCn+JaHlepQU8TjVL+8yfgLBHZhWuG4bUi8lhgQ1LK97QzgFJKKZ/SGo1SSimf0kSjlFLKpzTRKKWU8ilNNEoppXxKE41SSimf0kSjlFLKpzTRKKWU8qn/D6cP18vWCeZmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tgws_ana = gwaxion.tgw_approx(mbh, alphas, chi)\n",
    "\n",
    "plot(alphas, tgws_ana, label='Closed-form approximation', lw=3, c='gray', alpha=0.5)\n",
    "plot(alphas, tgws, label='Numerical result', lw=2, ls='--')\n",
    "\n",
    "xlim(0, 0.2);\n",
    "yscale('log');\n",
    "ylabel(r'$\\tau$ (s)');\n",
    "xlabel(r'$\\alpha$');\n",
    "\n",
    "legend(loc='best');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gwaxion",
   "language": "python",
   "name": "gwaxion"
  },
  "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.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}