gas-drag.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "# Binary black hole mergers in AGN disks -- Playing With Gas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### Chiara Mingarelli, cmingarelli@flatironinstitute.org"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# RUN IN PYTHON 3.XX\n",
    "from IPython.display import display, Math, Latex\n",
    "import math\n",
    "from math import sqrt, cos, sin, pi, log10\n",
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import FormatStrFormatter, LinearLocator, NullFormatter, NullLocator, MultipleLocator\n",
    "import matplotlib.ticker\n",
    "import matplotlib.colors\n",
    "from matplotlib.font_manager import FontProperties\n",
    "from matplotlib import rc, text\n",
    "\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = \"retina\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# physical constants for natural units c = G = 1\n",
    "c=2.99792458*(10**8)\n",
    "G=6.67428*(10**(-11))\n",
    "s_mass=G*(1.98892*10**(30))/(c**3)\n",
    "solar_mass_gr = 1.989e33\n",
    "secs_in_pc = 102927133.01 #about 1e8 s per pc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Mass functions\n",
    "\n",
    "def mu(m1,m2): return s_mass*(m1*m2)/(m1+m2) # reduced mass\n",
    "\n",
    "def M(m1,m2):  return s_mass*(m1+m2) # total mass\n",
    "\n",
    "def mchirp(m1,m2): return ((mu(m1,m2))**(3./5))*((M(m1,m2))**(2./5))  # chirp mass\n",
    "\n",
    "def mchirp_q(q,Mtot): \n",
    "    \"\"\"\n",
    "    chirp mass in terms of q and M_tot. Answer in seconds.\n",
    "    \"\"\"\n",
    "    ans = (q/(1+q)**2)**(3/5)*Mtot*s_mass\n",
    "    return ans\n",
    "\n",
    "def parsec2sec(d): return d*3.08568025e16/299792458\n",
    "\n",
    "def sec2pc(sec): return 299792458*sec/3.08568025e16\n",
    "\n",
    "# For GWs: strain, GW frequency and time to coalescence \n",
    "def generic_strain(q_mass_ratio, Mtot, dist, freq):\n",
    "    strain = sqrt(32./5)*mchirp_q(q_mass_ratio,Mtot)**(5/3)*(pi*freq)**(2/3)/parsec2sec(dist*1e6)\n",
    "    return strain\n",
    "\n",
    "def generic_strain_wMc(chirp_mass, dist, freq):\n",
    "    strain = sqrt(32./5)*(chirp_mass*s_mass)**(5/3)*(pi*freq)**(2/3)/parsec2sec(dist*1e6)\n",
    "    return strain\n",
    "\n",
    "def freq_gw(q, Mtot, tc):\n",
    "    \"\"\"\n",
    "    GW frquency as a function of time to coalescence in years, total mass and mass ratio\n",
    "    Result from integration of standard df/dt for GWs\n",
    "    \"\"\"\n",
    "    ans = mchirp_q(q,Mtot)**(-5/8)/pi*(256/5*tc*31556926)**(-3/8)\n",
    "    return ans\n",
    "\n",
    "def freq_gw_wMc(chirp_mass, tc):\n",
    "    \"\"\"\n",
    "    GW frquency as a function of time to coalescence in years and chirp mass (directly)\n",
    "    Result from integration of standard df/dt for GWs\n",
    "    \"\"\"\n",
    "    ans = (chirp_mass*s_mass)**(-5/8)/pi*(256/5*tc*31556926)**(-3/8)\n",
    "    return ans\n",
    "\n",
    "def time_to_c(q, Mtot, freq):\n",
    "    \"\"\"\n",
    "    time to coalescence of a binary in seconds\n",
    "    \"\"\"\n",
    "    ans = (pi*freq)**(-8/3)*mchirp_q(q,Mtot)**(-5/3)*5/256\n",
    "    return ans\n",
    "\n",
    "def time_to_c_wMc(chirp_mass, freq):\n",
    "    \"\"\"\n",
    "    freq. in Hz, input chirp mass in solar masses, answer in years\n",
    "    \"\"\"\n",
    "    ans = (pi*freq)**(-8/3)*(chirp_mass*s_mass)**(-5/3)*5/256\n",
    "    return (ans/31556926)\n",
    "\n",
    "def a_to_f(Mtot, fgw):\n",
    "    # Kepler's 3rd law\n",
    "    # answer in seconds\n",
    "    frac = (Mtot*s_mass/(pi*fgw)**2)**(1./3)\n",
    "    return frac\n",
    "\n",
    "def fgw_to_a(Mtot, a_bh):\n",
    "    \"\"\"\n",
    "    Kepler's 3rd law\n",
    "    input BH mass in solar masses\n",
    "    a in pc\n",
    "    answer in seconds\n",
    "    \"\"\"\n",
    "    ans = (Mtot*s_mass/parsec2sec(a_bh)**3)**(0.5)/pi\n",
    "    return ans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# lisa starts at 0.5 AU\n",
    "# ends at around 1e-2 AU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Hardening and Dynamical Friction Timescales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Black hole merger timescales from galaxy merger timescale; Binney and Tremaine 1987 \n",
    "# \"Galactic Dynamics\"; also Sesana and Khan 2015\n",
    "\n",
    "def t_hard(Mstar,q,gamma,Mtot):\n",
    "    \"\"\"\n",
    "    Hardening timescale with stars, Eq. 7 Sesana & Khan 2015\n",
    "    Answer in yrs\n",
    "    \"\"\"\n",
    "    a_val = parsec2sec(r0_sol(Mstar, gamma))\n",
    "    H = 15\n",
    "    aStarGW = a_StarGW(Mstar,q,Mtot,gamma,H) #check units\n",
    "    sigma_inf = sigmaVel(Mstar)*1000/c\n",
    "    rinf_val = r_inf(Mstar,gamma,Mtot)\n",
    "    rho_inf = rho_r(Mstar, gamma, rinf_val)\n",
    "    ans = sigma_inf/(H*rho_inf*aStarGW)\n",
    "    return ans/31536000\n",
    "\n",
    "def t_gw(a, q, Mtot):\n",
    "    # ans in yr\n",
    "    # a in pc\n",
    "    const = 5/256\n",
    "    num = parsec2sec(a)**4\n",
    "    deno = (Mtot*s_mass)**3 * q/(1+q)**2\n",
    "    ans = const*num/deno\n",
    "    return ans/31556926\n",
    "\n",
    "def soundSpeed(Temp):\n",
    "    \"\"\"\n",
    "    This is the sound speed, c_s, where c_s^2 = k_B*T/m_p\n",
    "    Returns dimensionless quantity. To get physical units, multiply by c.\n",
    "    \"\"\"\n",
    "    kbTemp_GR = 3.80526e-76*Temp # seconds, for 1K\n",
    "    mProton = 4.14324e-63 #seconds\n",
    "    ans = kbTemp_GR/mProton\n",
    "    return sqrt(ans)\n",
    "    \n",
    "def Mdot_to_Sigma(a_bin, alpha, Temp, SigmaSurf, bin_Mass):\n",
    "    \"\"\"\n",
    "    Parameters from Baruteau et al. 2011.\n",
    "    SigmaSurf input in g/cm^2 (change to kg/m**2)\n",
    "    input binary mass in solar masses (converts to seconds)\n",
    "    input semimajor axis in pc (converts to seconds)\n",
    "    answer is kg/s\n",
    "    for solar masses, divide by 1.98892*1e30.\n",
    "    \"\"\"\n",
    "    SigmaSurf = SigmaSurf*(1e-3)/(1e-2)**2\n",
    "    num = 3*pi*alpha*(soundSpeed(Temp)*c)**2*SigmaSurf\n",
    "    Omega = (bin_Mass*s_mass/(parsec2sec(a_bin)**3))**0.5 ## parsec to seconds, make sure all in seconds\n",
    "    deno = Omega\n",
    "    ans = num/deno #kg/s\n",
    "    return ans\n",
    "\n",
    "def a_bin(a_bhs,qmass):\n",
    "    \"\"\"\n",
    "    Set a_bin to be the Hill radius, R_H ~ a(q/3)**1/3 for circular systems\n",
    "    a_bhs is the semimajor axis of the stellar mass BHs, in pc\n",
    "    answer in seconds\n",
    "    \"\"\"\n",
    "    return parsec2sec(a_bhs)*(qmass/3)**(1/3)\n",
    "\n",
    "def t_ACC(a_bhs,M1,M2,SigmaSurf,Temp, alpha):\n",
    "    \"\"\"\n",
    "    temperature in Kelvin\n",
    "    SigmaSurf input in g/cm^2\n",
    "    BH masses in solar masses\n",
    "    a_bhs in pc; start with 1e-3pc\n",
    "    Answer in YEARS\n",
    "    \"\"\"\n",
    "    SigmaSurf = SigmaSurf*(1/solar_mass_gr*s_mass)/(1e-2/c)**2\n",
    "    a0 = sqrt(M(M1,M2)*s_mass/SigmaSurf)\n",
    "    num = mu(M1,M2)*M(M1,M2)**0.25\n",
    "    q=M2/M1\n",
    "    deno = a_bin(a_bhs,q)*3*pi*alpha*soundSpeed(Temp)**2*SigmaSurf**0.75\n",
    "    ans = num/deno\n",
    "    return ans/31556926 #yrs\n",
    "    \n",
    "def t_baruteau(Temp, SigmaSurf, bin_Mass, SMBH_mass, a_SMBH, a_bin):\n",
    "    \"\"\"\n",
    "    temperature in Kelvin\n",
    "    SigmaSurf input in g/cm^2\n",
    "    SMBH_mass in solar masses\n",
    "    bin_mass in solar masses\n",
    "    a_SMBH in pc\n",
    "    a_bin in AU\n",
    "    Answer in YEARS\n",
    "    \"\"\"\n",
    "    const = 1/5.443310539518173e-10\n",
    "    ans = const*(Temp/2000)**3*(SigmaSurf/100)**(-1)*(SMBH_mass/3e6)**(-0.5)*(bin_Mass/2)**(-2)*(a_SMBH/0.1)**(3/2)*(a_bin/50)\n",
    "    return 900*ans\n",
    "  \n",
    "def t_sh(Mstar,q,gamma, Mtot, a_bin):\n",
    "    # answer in yr\n",
    "    sig = sigmaVel(Mstar) # M-sigma\n",
    "    H = 15\n",
    "    inf_rad = r_inf(Mstar,gamma,Mtot) # Influence radius\n",
    "    rho_influence = rho_r(Mstar, gamma, inf_rad) #density profile at influence radius\n",
    "    num = sig*1000/c\n",
    "    deno = H*parsec2sec(a_bin)*rho_influence\n",
    "    ans = num/deno\n",
    "    return ans/31556926\n",
    "\n",
    "def dy_fric_a(Mstar,M2, a):\n",
    "    \"\"\"\n",
    "    Final eq from https://webhome.weizmann.ac.il/home/iair/astrocourse/tutorial8.pdf\n",
    "    returns timescale in yr\n",
    "    Mbh should be mass of secondary\n",
    "    a in kpc\n",
    "    \"\"\"\n",
    "    # assume log(Lambda) = 10\n",
    "    vc = sqrt(2)*sigmaVel(Mstar)\n",
    "    ans = 2.64e10*(a/2)**2*(vc/250)*(1e6/M2)\n",
    "    return ans\n",
    "\n",
    "def fmax(Mtot):\n",
    "    const = 1/(np.pi*6**(3/2))\n",
    "    ans = const/(Mtot*s_mass)\n",
    "    return ans\n",
    "\n",
    "def find_nearest(array,value):\n",
    "    #returns index of matching value in array\n",
    "    idx = (np.abs(array-value)).argmin()\n",
    "    return idx \n",
    "\n",
    "def pc_to_AU(parParam):\n",
    "    \"\"\"\n",
    "    change parsecs to AU\n",
    "    input parsecs\n",
    "    output in AU\n",
    "    \"\"\"\n",
    "    return parParam*206265\n",
    "\n",
    "def AU_to_pc(auParam): \n",
    "    \"\"\"\n",
    "    change AU measurement to parsecs\n",
    "    \"\"\"\n",
    "    return auParam*4.8481e-6\n",
    "\n",
    "def a_to_Om(Mtot, semi_a):\n",
    "    # Kepler's 3rd law\n",
    "    # answer in seconds\n",
    "    semi_a = parsec2sec(semi_a) #parsecs to seconds\n",
    "    frac = Mtot*s_mass/(semi_a)**3\n",
    "    ans = frac**(1/2)\n",
    "    return ans\n",
    "\n",
    "def tempCheckwUnits(temp):\n",
    "    mp = 938.23e6/c**2 #Mev/c^2\n",
    "    kB = 8.617e-5 #eV/K\n",
    "    ans = sqrt(kB*temp/mp)\n",
    "    return ans\n",
    "\n",
    "def fdot_gw(Mchirp,fgw):\n",
    "    const = 96./5*pi**(8./3)\n",
    "    num = (Mchirp*s_mass)**(5./3)*fgw**(11./3)\n",
    "    ans = num*const\n",
    "    return ans\n",
    "\n",
    "def fdot_gw_mod(m1,m2,fgw):\n",
    "    const = 96./5*pi**(8./3)\n",
    "    num = (mchirp(m1,m2))**(5./3)*fgw**(11./3)\n",
    "    ans = num*const\n",
    "    return ans\n",
    "\n",
    "def fdot_drag(m1,m2,fgw, temperature=1e5,rho_0=1e-6, Const=1):\n",
    "    rho_0= rho_0*(1e-3)/(0.01)**3# cm -> m\n",
    "    num = 2*rho_0*Const*3*(6.67428*(10**(-11)))**(5/3.)*((m1+m2)*2e30)**(2/3.)*fgw**(2./3)\n",
    "    deno = (soundSpeed(temperature)*c)**2*2*pi**(1/3.)\n",
    "    ans = num/deno\n",
    "    return ans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gw_vec = []\n",
    "drag_vec = []\n",
    "xax = []\n",
    "m1=0.5\n",
    "m2=0.3\n",
    "temperature = 9e5\n",
    "Const=1\n",
    "rho_0 = 1e-6\n",
    "\n",
    "\n",
    "for qq in range (1,10000):\n",
    "    qq=qq/10000\n",
    "    xax.append(qq)\n",
    "    gw_vec.append(fdot_gw_mod(m1,m2,qq))\n",
    "    drag_vec.append(fdot_drag(m1,m2,qq,temperature,rho_0, Const))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAAIqCAYAAAAaWboDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3iU15n+8e8Z9U4VRQLTu+SCO4m7MQZjYyD5ORs7jtM2a2/axgFssA0G20CKEztxNk4267RNdkMxrnFvMa64SPQOEogiQL3PnN8fMxrNyBoJmJFezej+XNdcIz3zzpyHFNCt855zjLUWERERERERaeFyugEREREREZHuRkFJRERERESkFQUlERERERGRVhSUREREREREWlFQEhERERERaUVBSUREREREpBUFJRERERERkVYUlERERERERFpRUBIREREREWlFQUlERERERKQVBSUREREREZFWFJRERERERERaiXe6AQmfMWYPkAnsdbgVEREREZHuZBhQYa0dfqpvVFCKDZkpKSl9xo8f38fpRkREREREuostW7ZQW1t7Wu9VUIoNe8ePH99nw4YNTvchIiIiItJtTJ48mY8++mjv6bxXa5RERERERERaUVASERERERFpRUFJRERERESkFQUlERERERGRVhSUREREREREWlFQEhERERERaUVBSUREREREpBUFJRERERERkVYUlERERERERFpRUBIREREREWlFQakLGWP2GmNsG49nne5NRERERERaxDvdQA9zHhAX8P0gYAPwf860IyIiIiIibVFQ6kLW2qOB3xtjvg5UAH93piMREREREWlLTN96Z4yZa4x51BjzljGmwneb2587eE+uMeb3xpiDxph63+1yPzfG9I5wbwb4OvBna21NJD9bRERERETCE+szSouAM4EqoBgY197FxpiRwHogG1gHbAXOB74HTDPGTLHWHotQb1cDw4HfRejzREREREQkQmJ6Rgn4ATAGyAT+7SSufwxvSPqutXaWtXaBtfYK4GFgLPBA4MXGmGUhNmcIfFwWYqxvAh9Yaz853T+ciIiIiIh0jpieUbLWvtb8tfdOt9CMMSOAqcBe4FetXr4P+BZwizHmh9baal/950C7t/IB+9sYKxu4Abijg/eKiIiIiEStwuJysjOTGJCZ7HQrpyymg9IpusL3/KK11hP4grW20hjzNt4gdSHwiq9eCpSexlhfBeqBv53Km4wxG0K81O4thSIiIiIiXamu0c3DL2/nt2/u5srxA3j8lskdTlx0N7F+692pGOt73h7i9R2+5zHhDOLbxOEbwN+stZXhfJaIiIiISHfz7u5jTPv5m/zmjd14LLy0+TDPFpY43dYp04xSiyzfc3mI15vrvcIc5zJgNHDzqb7RWju5rbpvpumc8NoSERERETl9lXWNLH9+K395L3jlyUUj+pKXkxXiXd2XgtLJa54rtOF8iG/dVHTNO4qIiIiItOPVrYdZuHYjJeV1/lpGUjwLZ4zn/503JOpuuwMFpUDNM0ah4m5mq+tERERERHq0Y1X13P/MZtZ9cjCoftX4ASybNYmBWdG3iUMzBaUW23zPodYgjfY9h1rDJCIiIiLSI1hreerTgyx5ejPHqxv89b5piSy5YSIz8gZF5SxSIAWlFs1biU81xrgCd74zxmQAU4Ba4F0nmhMRERER6Q5KymtZtHYjr2w9ElSffXYO91w3gd5piQ51FlkKSj7W2l3GmBfxbgF+B/BowMtLgDTgNwFnKImIiIiI9Bgej+V/3t/P8ue3UlXf5K8Pzkrmgdl5XD4228HuIi+mg5IxZhYwy/ftQN/zRcaYJ3xfl1pr7wx4y+3AeuARY8yVwBbgAuByvLfcLez0pkVEREREupk9pdXMX13A+3uOB9W/ctEZzJs2jvSk2IsVsfcnCnYWcGur2gjfA2Af4A9Kvlmlc4H7gWnAdKAEeARYYq09joiIiIhID9Hk9vC7f+7h4Ze2U9/kX5nCiH5pLJ+Tz/nD+zjYXeeK6aBkrV0MLD7F9xQBt3VGPyIiIiIi0WLzwQrmrf6UjQcq/LU4l+FfLxnBd68cTXJCnIPddb6YDkoiIiIiInJq6pvc/PLVnfz69V00eVqOEJ04OJMVc/KZFIWHx54OBSUREREREQFgw77jzFtVwK6jLfuXJca7+MFVY/jG54eTEOdysLuupaAkIiIiItLDVdc38eMXtvGHd/ZiWyaROH9YHx6ak8fI/umO9eYUBSURERERkR7sje1HuXtNIQfKav21tMQ4Fkwfz5fPH4rLFd0Hx54uBSURERERkR6orKaBpc9sYfVHxUH1y8b254Eb88jpleJQZ92DgpKIiIiISA9ireX5jYe4d91GSqsa/PXeqQncN3MiN5w1GGN65ixSIAUlEREREZEe4khFHfes28gLmw4H1WeeOZj7Zk6gX3qSQ511PwpKIiIiIiIxzlrL3z8sZumzm6msa/LXB2QmsWxWHldPGOBgd92TgpKIiIiISAzbf6yGu9YW8PbOY0H1L50/lLumjyMzOcGhzro3BSURERERkRjk9lieWL+Xn7ywjdpGt79+Rt9UHpqdx8Uj+znYXfenoCQiIiIiEmO2H65k3qoCPikq89dcBr7x+RH84KoxpCTGOdhddFBQEhERERGJEQ1NHn79+i5++doOGt0tJ8eOG5jBijn5nDmkl4PdRRcFJRERERGRGPBJURnzVxWw7XClv5YQZ/jOFaP59qUjSYx3Odhd9FFQEhERERGJYrUNbn764jZ+//YePC2TSJw9tBcr5+QzekCGc81FMQUlEREREZEotX5nKQvWFLL/eI2/lpIQx4+uGcutFw8jzqWDY0+XgpKIiIiISJQpr23koee28LcPioLqnxvVj4dm5zGkT6pDncUOBSURERERkSjy4qZDLHpyI0cq6/21zOR4Fl03gS9MzsUYzSJFgoKSiIiIiEgUOFpZz+KnN/FsQUlQfdrEgdx/w0SyM5Md6iw2KSiJiIiIiHRj1lrWfnyA+5/ZTFlNo7/eLz2JpTdM5Nq8QQ52F7sUlEREREREuqkDZbXcvaaQN7YfDap/YXIuC2eMp1dqokOdxT4FJRERERGRbsbjsfz5vX2seH4r1Q1ufz23dwoP3pjHJWP6O9hdz6CgJCIiIiLSjew6WsWC1QV8sPeEv2YMfPXiYdw5dSxpSfoRvivoP2URERERkW6g0e3h8Td384tXdtDQ5PHXR2Wns2JOPpPP6O1gdz2PgpKIiIiIiMM2Hihn3qoCNpdU+GvxLsPtl43kjitGkRQf52B3PZOCkoiIiIiIQ+oa3fzilR08/uZu3B7rr+flZLFybj7jB2U62F3PpqAkIiIiIuKA9/ccZ8HqAnaXVvtrSfEufjh1DF+bMpz4OJeD3YmCkoiIiIhIF6qsa2TlP7bxp3f3BdUvGN6HFXPyGdYvzaHOJJCCkoiIiIhIF3lt6xEWri3kYHmdv5aRFM9d08dz03lDcLmMg91JIAUlEREREZFOdry6gaXPbGbtxweC6leOy2bZjZMYlJXiUGcSioKSiIiIiEgnsdbyTEEJi5/axLHqBn+9T1oii6+fyMz8QRijWaTuSEFJRERERKQTHCqvY9GTG3l5y+Gg+qyzBnPvzIn0SUt0qDM5GQpKIiIiIiIR5PFY/vZBEQ89t4XK+iZ/fVBWMg/emMfl47Id7E5OloKSiIiIiEiE7C2tZsGaAt7dfTyofvOFQ5k/bRwZyQkOdSanSkFJRERERCRMTW4Pv397Dz99cTv1TR5/fXi/NJbPzuOCEX0d7E5Oh4KSiIiIiEgYtpRUMH91AQXF5f5anMvwzc+P4PtXjSY5Ic7B7uR0KSiJiIiIiJyG+iY3v3p1J4+9vosmj/XXJwzKZOXcfCblZDnYnYRLQUlERERE5BRt2HeC+asL2Hmkyl9LjHfxvStH861LRpAQ53KwO4kEBSURERERkZNUXd/ET17cxhPr92JbJpE494zeLJ+Tz6jsdOeak4hSUBIREREROQlv7TjKXWsKKT5R66+lJcYx/9px3HzBGbhcOjg2ligodSFjTAawFLgRyAY+Br5nrf3A0cZEREREJKTymkaWPbuZv28oDqpfOqY/D9w4idzeqQ51Jp1JQalr/Q7IB24FioGbgZeNMROstQcc7UxEREREPuMfG0u4Z90mjlbW+2u9UhO497oJ3Hh2DsZoFilWKSh1EWNMCjAHmGOtfd1XXmyMmQn8G7DIqd5EREREJNiRyjruW7eJ5zceCqrPyB/E4pkT6Z+R5FBn0lViOigZY+YClwJnAWcCGcBfrLU3t/OeXOB+YBrQFygBngSWWGtPhNFOPBAH1LWq1wKfC+NzRURERCRCrLWs2lDM0mc2U1HX5K9nZySxdNYkrpk40MHupCvFdFDCO0tzJlCF91a3ce1dbIwZCazHu35oHbAVOB/4HjDNGDPFWnvsdBqx1lYaY94BFhljNgKHgC8BFwE7T+czRURERCRyio7XcPfaQt7aURpUv+m8Idw1fTxZKQkOdSZOiPUN3n8AjAEy8d7e1pHH8Iak71prZ1lrF1hrrwAeBsYCDwRebIxZZoyxHTwuC3jLLYAHb2irB74L/BVwh/nnFBEREZHT5PZYfv/PPUx9+M2gkDS0Typ/+cYFLJ+Tr5DUA8X0jJK19rXmrztaaGeMGQFMBfYCv2r18n3At4BbjDE/tNZW++o/B/7cQRv7A/rZBVxqjEkDMq21JcaY/wX2dPynEREREZFI23G4kvmrC/hof5m/5jLwtSnD+Y+pY0hNjOkfl6Ud+m++xRW+5xettZ7AF3y3zb2NN0hdCLziq5cCwXOzJ8EXtKqNMb2Ba4B5J/M+Y8yGEC+1e0uhiIiIiARraPLwmzd28eirO2lwt/zoN2ZAOivm5HP20N4OdifdgYJSi7G+5+0hXt+BNyiNwReUTpUx5hq8tztuBUYBPwa2Af99Op8nIiIiIqeuoLiMeasK2Hqo0l9LiDPccfkobr9sFInxsb46RU6GglKLLN9zeYjXm+u9whzjISAXOA6sBhZaaxtP5s3W2slt1X0zTeeE0ZeIiIhIzKttcPPwy9v53Vu78diW+llDerFybj5jBmQ415x0OwpKJ695kZNt96p2WGv/D/i/yLQjIiIiIifrnV3HWLCmgH3Havy1lIQ47rxmLF+9eBhxLh0cK8EUlFo0zxhlhXg9s9V1IiIiItLNVdQ18tBzW/nr+/uD6lNG9eWhG/MZ2jfVoc6ku1NQarHN9zwmxOujfc+h1jCJiIiISDfy8ubDLHyykMMV9f5aRnI8i2aM54vnDulwV2Tp2RSUWjRvJT7VGOMK3PnOGJMBTAFqgXedaE5ERERETk5pVT1Lnt7M058eDKpPnTCApbMmMSAz2aHOJJooKPlYa3cZY17Eu7PdHcCjAS8vAdKA3wScoSQiIiIi3Yi1lnWfHGTJ05s4UdOyV1a/9ETuv2ES104aqFkkOWkxHZSMMbOAWb5vB/qeLzLGPOH7utRae2fAW24H1gOPGGOuBLYAFwCX473lbmGnNy0iIiIip+xgWS0L1xby2rajQfU55+SyaMZ4eqclOtSZRKuYDkrAWcCtrWojfA+AfYA/KPlmlc4F7gemAdOBEuARYIm19nindywiIiIiJ83jsfzl/f2seH4rVfVN/npOrxQenJ3HpWP6O9idRLOYDkrW2sXA4lN8TxFwW2f0IyIiIiKRs/toFQtWF/L+3pbfZRsDt140jDuvGUt6Ukz/qNt9uZugsgTKi6G8CGqOwYX/5nRXp0z/6xERERGRqNLk9vDbt/bw8MvbaWjy77/FyP5prJiTz7nD+jjYXQ9QV+4LQb4g1Px1me/ryoPQsi8aYODcr0N8dN3+qKAkIiIiIlFj08Fy5q8uYOOBCn8t3mX49qUj+fcrRpGcEOdgdzHA3QRVhwKCT0AQan7Un+qxotYbnnoP64yOO42CkoiIiIh0e3WNbh59dQf/+cZu3B7rr0/KyWTFnHwmDs5ysLsoUlfRxmxQQBiqOAjWHf44adnQawhk5ULWEIiLrtkkUFASERERkW7ug73Hmb+6gN1HW05pSYp38YOrx/CNzw0nPs7lYHfdiMcNlYcCwk/r2aAi721z4YpP9gWg5sfQ4O8zcyAh+s+qUlASERERkW6pqr6Jlf/Yyh/f2RdUP394H5bPzmNE/3SHOnNIfeVnZ4PKAmeDDkRuNsgffIa0fN1riPf71L7eXTNinIKSiIiIiHQ7r207wsI1hRwsr/PX0pPiWXDtOP7l/KG4XDH2g7rHDVWHQ68LKt8fmdmguKTgENRrSPD3mYMhISX8cWKAgpKIiIiIdBsnqhtY+sxm1nx8IKh+xbhsls2axOBeUfpDfH1VO2uDirxrgzxNHX9OR1L7Ba8Nan2LXFq/HjEbFAkKSiIiIiLiOGstzxUe4r6nNlJa1eCv905NYPH1E7n+zMGY7voDvsfjnQ1qnvlpvS6ovBhqT4Q/Tlxiq1viWs0GZeVoNiiCFJRERERExFGHK+pY9ORGXtp8OKh+w1mDufe6CfRNT3KoM5+G6tBnBvlngxrDHye1Xwdrg/qBSxtXdBUFJRERERFxhLWW//2giAee20JlXcttZwMzk3ngxklcOX5A5zfh8UD1kRBrg3y1SM0GZea0szYoBxJTwx9HIkZBSURERES63L5j1dy1ppD1u44F1b98wVDmXzuOzOSEyAzUUA3lB1qFoMA1QgciNBvUN/S6oKxcSOuv2aAoo6AkIiIiIl3G7bH899t7+MmL26hr9Pjrw/qmsnxOPheO6HvyH+bxQPXR0GcGlRVB7fHwm3YleNf/fGZdUMDaoMS08MeRbkVBSURERES6xLZDlcxbXcCnRWX+msvANz8/gu9fNYaUxLjgNzTUeM8G+syZQUUt5wa5GwhbSp+2Z4N6Nc8GZWs2qAdSUBIRERGRTlXf5Oax13bx2Os7aXRbAAweLsj2sPSyLEYnbYMPX/nsbXE1xzr45JPgivetDRoSsDFCq7VBST3s4Fo5KQpKIiIiIhJ5jbVQfoCdO7aw7s33iKs8wIOmlMEJx8hxHSPXdZz4igZ4KsxxUnq3vzYoPRtccR1/jkgrCkoiIiIicmqsherS0GcGlRVBTSkAo4Afwmd/6rQnMY4rHjIHhzgzKNe7NigpI7J/NhEfBSURERERCdZY1/7aoPJicNeHP05yrzZmggLWBqUP0GyQOEZBSURERKQnsda79qcsxGxQeZF3J7kwNVkXh+jDAduPxvQc8idMJHPg8JZglJkDyZkR+AOJdA4FJREREZFY4p8NanVoamAoaqoLf5zkrKDZoO11Wfxhs4ctNVkctH05Qm8yUpK457oJzDknB2NM+GOKdCEFJREREZFo0TwbFOrMoPJiqD4S/jgmzrc2qPUmCQFf+2aDjlbWs/ipTTxbWBL0EdPzBrL4+olkZySH34+IAxSURERERLqLpvrg2aCytmaDasMfJykr9LqgrFxIHwhx7f+YaK1l9UcHWPrMZsprG/31/hlJLL1hItMmDQq/TxEHKSiJiIiIdAVrofZEq7VBrTZJqDoc/jjGBRmD2z4zyD8blBXWEEXHa7h7bSFv7SgNqn/x3FwWTp9AVmpCWJ8v0h0oKImIiIhEQlNDx2uDGmvCHycxo1UACjgzKCsXMgZ1OBt0ujweyx/f2cvKF7ZR0+D213N7p7B8dj6fG92vU8YVcYKCkoiIiEhHmmeDgoJPq62zqw5zcocDtcO4vEEn1NqgXkPCng06XTuPVDJ/dSEb9p1oadfAbRcP585rxpCaqB8rJbbof9EiIiIiTQ1QebD9tUGN1eGPk5j+2XODerWeDepet601uj385o1dPPLKThrcHn99dHY6K+bmc87Q3g52J9J5FJREREQktvlng0KcGVReDJWHCHs2COMNOqHWBWX5ZoOiaJvswuJy5q0uYEtJhb8W7zLcfvko7rh8JEnxOgxWYpeCkoiIiEQ3dyNUHAy9Lqi8GBqqwh8nIa39tUGZg7vdbNDpqmt08/DL2/ndW3twe1oC5Jm5WayYm8+4gTooVmKfgpKIiIh0X9ZCXVnbs0HNa4MqS4jMbNDAVrNAAaGo1xBI7hVVs0Gn693dx7hrTSF7SltuNUxOcHHn1LHcNmU4ca7Y/89ABBSURERExEnuRm/QCbUuqLwYGirDHych9bO3wQXODmUMhvjE8MeJYpV1jSx/fit/eW9/UP2iEX1ZPiePM/qmOdSZiDMUlERERKTz1Ja1sS4ocG1QCVhPx5/TkfSB7a8NSundI2aDTterWw+zcO1GSsrr/LWMpHgWzhjP/ztvCEb/2UkPpKAkIiIip8fd1DIbFGptUH1Fx5/TkfiUNkLQkOC1QfFJ4Y/TAx2rquf+Zzaz7pODQfWrxg9g2axJDMxKdqgzEecpKImIiEjb6srbng3yrw06GKHZoAFtnxnUvHW2ZoMizlrLU58eZMnTmzle3eCv901LZMkNE5mRN0izSNLjKSiJiIj0RO4mqDrU/tqg+vLwx4lPCX1mUFYuZOZoNqiLlZTXsmjtRl7ZeiSoPvvsHO65bgK903r2Wi2RZgpKIiIisaiuIvSZQeXF3u20rTv8cdKyA26La2Pr7NQ+mg3qJjwey18/2M9Dz22lqr7JXx+clcwDs/O4fGy2g92JdD8KSiIiItHG4/YekBpybVCR97a5cMUnhz4zqHk2KEFrWKLBntJqFqwu4L09x4PqX7noDOZNG0d6kn4kFGlN/68QERHpbuorQ58ZVF4MFQciNBvUv+0zg5pvkUvtq9mgKNfk9vBf/9zDz17aTn1Ty3qyEf3SWD4nn/OH93GwO5HuTUFJRESkK3ncUHXYF372f3ZdUPn+yMwGxSWFPjMoa4h3p7iElPDHkW5r88EK5q8uoPBAy/+e4lyGf71kBN+9cjTJCXEOdifS/SkoiYiIRFJ9VTtrg4q8a4M8TR1/TkdS+7W/Niitn2aDeqj6Jje/fHUnv359F00e669PHJzJijn5TMrJcrA7keihoCQiInKyPJ6W2aDy1rNBvjBUeyL8ceISQ58ZlDUEsnI0GyRt2rDvOPNWFbDraLW/lhjv4gdXjeEbnx9OQpzLwe5EoouCkoiISLOG6tBnBvlngxrDHye1X4gzg3yhKLUfuPQDrZy86vomfvzCNv7wzl5syyQS5w3rzfI5+Yzsn+5YbyLRSkEpgowxlwB3ApOBwcBt1tonTvUaERHpBB4PVB8JcWaQrxap2aDMnHbWBuVAYmr444j4vLH9KHevKeRAWa2/lpYYx4Jrx/HlC87A5dItmCKnQ0EpstKBjcAffY/TvUZERE5VQzWUH2gVggLXCB2I0GxQ39DrgrJyvTvJaTZIukBZTQNLn9nC6o+Kg+qXje3PAzfmkdNLt2eKhENBKYKstc8BzwEYY5443WtERKQVjweqj4Y+M6isCGqPd/w5HXEleNf/tLlltm9tUGJa+OOIhOm5whLuXbeR0qoGf613agL3zpzArLNyMNrIQyRsUR2UjDFzgUuBs4AzgQzgL9bam9t5Ty5wPzAN6AuUAE8CS6y1EbjnQkRETllDjfdsoM+sDSpqOTfI3dDx53QkpU/oM4OyciEtW7NB0q0dqajjnnUbeWHT4aD6zDMHc9/MCfRLT3KoM5HYE9VBCViENyBVAcXAuPYuNsaMBNYD2cA6YCtwPvA9YJoxZoq19lindiwi0tN4PFBT2sbaoIDb4moi8FevK963NmhIwMYIrdYGJWlBu0Qnay1//7CYZc9upqKuZXv5AZlJLJuVx9UTBjjYnUhsivag9AO8AWkn3pml1zq4/jG8Iem71tpHm4vGmJ/5PusB4NsB9WXAwg4+83Jr7eun3LmISKxorA2xNqj5+wPgrg9/nJTe7a8NSs8Glw7QlNiz/1gNd68t5J87S4PqXzp/KHdNH0dmcoJDnYnEtqgOStZafzDq6F5cY8wIYCqwF/hVq5fvA74F3GKM+aG1tvnwgZ8Df+6gjf2n0LKISHSxFqpLQ58ZVFbknS0KlyseMgeHODMo17s2KCkj/HFEoojbY3li/V5+8sI2ahvd/voZfVN5aHYeF4/s52B3IrEvqoPSKbrC9/yitdYT+IK1ttIY8zbeIHUh8IqvXgpE4CeAyDDGbAjxUru3HIqIhNRYF7w2KPDMoOZQFInZoORebcwEBawNSh+g2SCRANsPVzJvVQGfFJX5ay4D3/j8CH5w1RhSEvX/F5HO1pOC0ljf8/YQr+/AG5TG4AtKp8oYkw6M8n3rAoYaY84Cjltr95/sNSIiEWGtd+1PWYjZoPIi705y4TJxLecGtV4XlJXrfS05M/xxRHqAhiYPv359F798bQeN7paTY8cNzGDFnHzOHNLLwe5EepaeFJSyfM/lIV5vrofzN9C5BK+TWuJ7/AH46ilc0yZr7eS26r6ZpnNOp2ERiWL+2aC21gX5Hk114Y+TnNXGbFDALXIZAzUbJBIBnxSVMX9VAdsOV/prCXGG71wxmm9fOpLEeO3IKNKVelJQ6kjzIifb7lXt8G3q0O5iqZO5RkTEPxsU6syg8mKoPhL+OCbOtzaorTODfF9rNkikU9U2uPnpi9v4/dt78AT8FHL20F6smJPPmAFanyfihJ4UlJpnjLJCvJ7Z6joRkc7TVB88G/SZrbOLoak2/HGSsj67LihrSMstcukDIa4n/VMg0r2s31XKgtWF7D9e46+lJMTxo2vGcuvFw4hz6XerIk7pSf86bvM9jwnx+mjfc6g1TCIiJ8daqD3Ram1Qq62zqw53/DkdMS7IGBx6bVBWrve2ORHpdsprG1n+/Bb++n5RUP1zo/rx0Ow8hvRJdagzEWnWk4JS87qgqcYYV+DOd8aYDGAKUAu860RzIhJFmho6XhvUWNPx53QkMaNVAAo4MygrFzIGaTZIJAq9tPkwi54s5HBFy46SmcnxLLpuAl+YnNvhkSci0jV6zL+w1tpdxpgX8e5sdwfwaMDLS4A04DcBZyiJSE/UPBv0mcNTA7bOrjpMGMsZvYzLG3RCrQ3qNUSzQSIxprSqnsVPbeKZgpKg+rSJA7n/holkZyY71JmItCWqg5IxZhYwy/ftQN/zRcaYJ3xfl7LSlDMAACAASURBVFpr7wx4y+3AeuARY8yVwBbgAuByvLfcLez0pkXEWU0NUHmw/bVBjRH4fUli+md3iuvVejYoIfxxRKTbs9ay9uMD3P/MZspqGv31fulJLL1hItfmDXKwOxEJJaqDEnAWcGur2gjfA2Af4A9Kvlmlc4H7gWnAdKAEeARYYq093ukdi0jn8c8GhTgzqLwYKg8R9mwQxht0Qq0LyvLNBun2GZEe70BZLXevKeSN7cFnln1hci4LZ4ynV2qiQ52JSEeiOihZaxcDi0/xPUXAbZ3Rj4h0MncjVBwMvS6ovBgaqsIfJyGt/bVBmYM1GyQi7fJ4LH9+bx8rnt9KdYPbX8/plcJDs/O4ZEx/B7sTkZMR1UFJRGKItVBX1vZsUPPaoMoSIjMbNLDVLFCrUJTSW7NBInLadh2tYsHqAj7Ye8JfMwZuvWgYP7pmLGlJ+vFLJBro/6ki0jXcjd6g418btL+N2aDKjj+nIwmpn70NLnB2KGMwxOtWFxGJvEa3h8ff3M0vXtlBQ5N/c11G9k9j5dx8Jp/Rx8HuRORUKSiJSGTUlrWxLihwbVAJtOzKf/rSB7a/NkizQSLigI0Hypm3qoDNJRX+WrzLcPtlI7njilEkxcc52J2InA4FJRHpmLspYDYoxNqg+oqOP6cj8SlthKAhwWuD4pPCH0dEJELqGt384pUdPP7mbtyelluD83KyWDk3n/GDMh3sTkTCoaAkIlBX3vZskH9t0MEIzQYNaPvMoOatszUbJCJR5P09x1mwuoDdpS1HCiTFu/jh1DF8bcpw4uNcDnYnIuFSUBKJde4mqDoU+syg8mKoLw9/nPiU0GcGZeVCZo5mg0QkJlTWNbLyH9v407v7guoXDO/Dijn5DOuX5lBnIhJJCkoi0a6uIvSZQeXF3u20rbvjz+lIWnbAbXFtbJ2d2kezQSIS817beoSFaws5WF7nr2UkxXPX9PHcdN4QXC79PSgSKxSURLozj9t7QGrItUFF3tvmwhWf3Cr4tF4blAMJyeGPIyISpY5XN7D0mc2s/fhAUP3Kcdksu3ESg7JSHOpMRDqLgpKIk+or218bVHEgQrNB/UOfGdRrKKT21WyQiEgbrLU8U1DC4qc2cay6wV/vk5bI4usnMjN/EEZ/f4rEJAUlkc7icUPV4dBnBpUXeQ9YDVdcUugzg7KGeHeKS9BvOkVETtWh8joWPbmRl7ccDqrPOmsw986cSJ80nckmEssUlEROV31VO2uDirxrgzxN4Y+T2q/9tUFp/TQbJCISQR6P5W8fFPHQc1uorG/5e3xQVjIP3DiJK8YNcLA7EekqCkoibfF4WmaDylvPBvnCUO2J8MeJSwy9LihrCGTlaDZIRKQL7S2tZsGaAt7dfTyofvOFQ5k/bRwZyQkOdSYiXU1BSXqmhurQ64L8s0GN4Y+T2i/EmUG+UJTaD1w6Z0NExGlNbg///fZefvrSNuoaW86NG94vjeWz87hgRF8HuxMRJygoSezxeKD6SIgzg3y1SM0GZea0szYoBxJTwx9HREQ61dZDFcxfVcCnxS27iMa5DN/8/Ai+f9VokhPiHOxORJyioCTRp6Eayg+0CkGBa4QORGY2KKVPy8xP63VBWbneneQ0GyQiErXqm9z86tWdPPb6Lpo81l+fMCiTlXPzmZST5WB3IuI0BSXpXjweqD4aem1QWRHUHu/4czriSvCu/2lzy2zf2qBEnawuIhKrNuw7wfzVBew8UuWvJca7+N6Vo/nWJSNIiNMvwkR6OgUl6VoNNd6zgT6zNijg3CB3Q8ef05GUPqHPDMrKhbRszQaJiPRA1fVN/OTFbTyxfi+2ZRKJc8/ozfI5+YzKTneuORHpVhSUJHI8HqgpbWNtUMBtcTXHwh/HFe9bGzQkYGOEVmuDkvQPnYiIBHtrx1HuWlNI8Ylafy01MY7508Zxy4Vn4HLpqAURaaGgJOH56I9Q+HdfEDoA7vrwPzOld+gzg7JyIT0bXFpYKyIiJ6e8ppFlz27m7xuKg+qXjOnPgzdOIre3Nt4Rkc9SUJLwlBfDnjdP/npXPGQODnFmUK53bVBSRuf1KyIiPco/NpZwz7pNHK1s+UVer9QE7pkxgdnn5GB0YLeIhKCgJOHJyg3+PrlXGzNBAWuD0gdoNkhERDrdkco67lu3iec3Hgqqz8gfxOKZE+mfkeRQZyISLRSUJDwjr4Avr/KGoMwcSM50uiMREenBrLWs2lDM0mc2U1HX5K9nZySxdNYkrpk40MHuRCSaKChJeJpnjERERBxWdLyGu9cW8taO0qD6TecN4a7p48lKSXCoMxGJRgpKIiIiEtXcHssf39nLj1/YRk2D218f2ieVh2bnMWVUP+eaE5GopaAkIiIiUWvH4Urmry7go/1l/prLwNemDOc/po4hNVE/6ojI6dHfHiIiIhJ1Gpo8/OaNXTz66k4a3B5/fcyAdFbMyefsob0d7E5EYoGCkoiIiESVguIy5q0qYOuhSn8tIc5wx+WjuP2yUSTGuxzsTkRihYKSiIiIRIXaBjcPv7yd3721G49tqZ85pBcr5+QzdqDO4RORyFFQEhERkW7vnV3HWLCmgH3Havy15AQXd04dy21ThhPn0sGxIhJZCkoiIiLSbVXUNfLQc1v56/v7g+oXj+zL8tn5DO2b6lBnIhLrFJRERESkW3p582EWPlnI4Yp6fy0jOZ5FM8bzxXOHYIxmkUSk8ygoiYiISLdyrKqeJU9v5qlPDwbVp04YwNJZkxiQmexQZyLSkygoiYiISLdgrWXdJwdZ8vQmTtQ0+uv90hO5/4ZJXDtpoGaRRKTLKCiJiIiI4w6W1bJwbSGvbTsaVJ9zTi6LZoynd1qiQ52JSE+loCQiIiKO8Xgsf3l/Pyue30pVfZO/ntMrhQdn53HpmP4OdiciPZmCkoiIiDhi99EqFqwu5P29x/01Y+ArF57Bj6aNIz1JP6aIiHP0N5CIiIh0qSa3h9++tYeHX95OQ5PHXx/RP42Vc/I5d1gfB7sTEfFSUBIREZEus+lgOfNXF7DxQIW/Fucy/NulI/n3K0aRnBDnYHciIi0UlERERKTT1TW6efTVHfznG7txe6y/PiknkxVz8pk4OMvB7kREPsvldAOxxBhziTHmKWPMAWOMNcZ8tY1r7jDGFBhjKnyPd4wxMxxoV0REpEt8uPc40x95i1+9tssfkpLiXSy4dhxP3j5FIUlEuiXNKEVWOrAR+KPv0ZZiYD6wA29QvRV40hgz2Vpb0CVdioiIdIGq+iZ+/I+t/PHdfdiWSSTOH96H5bPzGNE/3bnmREQ6oKAUQdba54DnAIwxT4S4Zl2r0kJjzL8BFwEKSiIiEhNe33aEhWs3cqCs1l9LT4pnwbXj+Jfzh+Jy6eBYEeneojooGWPmApcCZwFnAhnAX6y1N7fznlzgfmAa0BcoAZ4EllhrT3R608G9xAFfwDsTtb4rxxYREekMJ6obWPrsZtZ8dCCofsW4bJbNmsTgXikOdSYicmqiOigBi/AGpCq8t7SNa+9iY8xIvIEkG1gHbAXOB74HTDPGTLHWHuvUjr195AHvAMl4e7/RWlvY2eOKiIh0FmstzxUe4r6nNlJa1eCv905NYPH1E7n+zMEYo1kkEYke0R6UfoA3IO3EO7P0WgfXP4Y3JH3XWvtoc9EY8zPfZz0AfDugvgxY2MFnXm6tff0U+96GdxasFzAH+IMx5jJr7cZT/BwRERHHHa6oY9GTG3lp8+Gg+vVnDua+mRPom57kUGciIqcvqoOStdYfjDr6LZUxZgQwFdgL/KrVy/cB3wJuMcb80Fpb7av/HPhzB23sP4WWAbDWNuANdwAfGmPOwxvUvn6qnyUiIuIUay3/+0ERDzy3hcq6Jn99YGYyy2ZN4qoJAxzsTkQkPFEdlE7RFb7nF621nsAXrLWVxpi38QapC4FXfPVSoLQLenMBHf66zRizIcRL7d5yKCIiEmn7jlVz15pC1u8KvmP9Xy4YyoJrx5GZnOBQZyIikdGTgtJY3/P2EK/vwBuUxuALSqfKGJMOjPJ96wKGGmPOAo5ba/f7rlkOPAsU4d184l+AywCdpSQiIt2e22P577f38JMXt1HX2PJ7x2F9U3lodj4XjezrYHciIpHTk4JS82l25SFeb673CmOMcwleJ7XE9/gD8FVfbSDe2/kG+sYsAK611r7Q0Ydbaye3VffNNJ1z2l2LiIichG2HKpm3uoBPi8r8NZeBb35+BN+/agwpiXEOdiciElk9KSh1pHmRk233qnb4NnVod7GUtfarp/v5IiIiTmho8vCr13by2Os7aXS3/DM5bmAGK+fmk58bzu8YRUS6p54UlJpnjLJCvJ7Z6joREZEe7+P9J5i/uoDth6v8tcQ4F9+9chT/eulIEuJcDnYnItJ5elJQ2uZ7HhPi9dG+51BrmERERHqMmoYmfvridn7/9h5swL0W5wztxcq5+YzKznCuORGRLtCTglLz2qGpxhhX4M53xpgMYApQC7zrRHMiIiLdxds7S1mwpoCi47X+WmpiHPOuGcstFw0jzqWDY0Uk9vWYoGSt3WWMeRHvznZ3AI8GvLwESAN+E3CGkoiISI9SXtvIg89u4X8/LAqqf350Px68MY8hfVId6kxEpOtFdVAyxswCZvm+Heh7vsgY84Tv61Jr7Z0Bb7kdWA88Yoy5EtgCXABcjveWu4Wd3rSIiEg39MKmQ9zz5EaOVNb7a1kpCdxz3QTmnJPT4cHuIiKxJqqDEnAWcGur2gjfA2Af4A9Kvlmlc4H7gWnAdKAEeARYYq093ukdi4iIdCNHK+tZ/NQmni0sCapPzxvI4usnkp2R7FBnIiLOiuqgZK1dDCw+xfcUAbd1Rj8iIiLRwlrL6o8OsPSZzZTXNvrr/TOSWHrDRKZNGuRgdyIizovqoCQiIiKnrvhEDXev3cib248G1b94bi4Lp08gKzXBoc5ERLoPBSUREZEewuOx/Ondfaz4x1ZqGtz+em7vFJbPzudzo/s52J2ISPeioCQiItID7DxSxYLVBXy474S/ZgzcdvFw7rxmDKmJ+pFARCSQ/lYUERGJYY1uD4+/uZtfvLyDBrf/CEFGZ6ezYm4+5wzt7WB3IiLdl4KSiIhIjCosLmfe6gK2lFT4a/Euw+2Xj+KOy0eSFB/nYHciIt2bgpKIiEiMqWt08/DL2/ndW3twe6y/np+bxcq5+YwbmOlgdyIi0UFBSUREJIa8u/sYd60pZE9ptb+WnODih1eP5bYpw4iPcznYnYhI9IhoUDLG/AvwhrX2QCQ/V0RERNpXWdfI8ue38pf39gfVLxzRh+Wz8xnWL82hzkREolOkZ5T+jPcA2Psj/LkiIiISwqtbD7Nw7UZKyuv8tYykeO6eMZ6bzhuCMcbB7kREolOX33pnjJkLnGmtvaerxxYREYklx6rquf+Zzaz75GBQ/arxA1g2axIDs5Id6kxEJPqFHZSMMTcBHwK7TvItE4G7AQUlERGR02Ct5alPD7Lk6c0cr27w1/umJbLkhonMyBukWSQRkTBFYkbpfwALVPmerzbGHAQ+AgqttY2trk8D6iMwroiISI9TUl7LorUbeWXrkaD67LNzuOe6CfROS3SoMxGR2BKJoHQncDYwGRgHTAEu9r3WaIzZjDc0fQxUAF8C9kVgXBERkR7D47H89YP9PPTcVqrqm/z1wVnJPDA7j8vHZjvYnYhI7Ak7KFlrf9b8tTHGAzwBvAOc43vkAWe1etu/hzuuiIhIT7GntJoFqwt4b8/xoPpXLjqDedPGkZ6k0z5ERCIt0n+zrgQ+tNauai4YY+KACXjDUl9gg7X2rQiPKyIiEnOa3B7+6597+NlL26lv8vjrI/qlsXxOPucP7+NgdyIisS2iQclau6CNmhso9D1ERETkJGw+WMH81QUUHij31+Jchn+9ZATfvXI0yQlxDnYnIhL7wgpKxphVwAa8a5A+stYejUhXIiIiPVR9k5tfvrqTX7++iyaP9dcnDMpk5dx8JuVkOdidiEjPEe6M0mzfwwIYYw7gC034ApS1tiTMMURERHqEDfuOM29VAbuOVvtrifEuvn/VaL75+REkxLkc7E5EpGcJNyidQcumDZN9z9cDM5svMMYc5rPhqSjMcUVERGJGdX0TP35hG394Zy+2ZRKJ84b1ZvmcfEb2T3esNxGRniqsoOQLPEXAuuaaMeZx4BvAFuAYMAKYDlwbcE2ptXZAOGOLiIjEgje3H+WuNYUcKKv119IS41hw7Ti+fMEZuFw6OFZExAkR3czBGPMj4MvAJdbafwbULwaWAFfiDVa27U8QERHpGcpqGlj27BZWbSgOql82tj8P3JhHTq8UhzoTERGI/PbgdwB/DQxJANba9cDVxpj5wL14D6UVERHpkZ4vLOGedZsorar313qnJnDvzAnMOisHYzSLJCLitEgHpWy8t9u1yVq7whgzG7gb+GKExxYREenWjlTUcc+6jbyw6XBQfeaZg7lv5gT6pSc51JmIiLQW6aC0A7iqg2teB74S4XFFRES6LWstf/+wmGXPbqairslfH5CZxLJZeVw9Qct2RUS6m0gHpd8DDxtj7rbWPhjimoGAjhIXEZEeYf+xGu5eW8g/d5YG1b90/hDumj6ezOQEhzoTEZH2RDoo/RKYASw1xlwF3NtqU4cZwE14N3QQERGJWW6P5Yn1e/nJC9uobXT760P7pLJ8dh4Xj+rnYHciItKRiAYla63bF4YeA74OvGGMOQbsx7t+KQcwwC8iOa6IiEh3sv1wJfNWFfBJUZm/5jLw9c8N5z+uHktKYpyD3YmIyMmI9IwS1tpG4JvGmN8C3wGuxnsQrRvYBjxkrf1TpMcVERFxWkOTh1+/votfvraDRnfLSRjjBmawYk4+Zw7p5WB3IiJyKiIelJpZa98HbgEwxiQBbmttU/vvEhERiU6fFpUxb1UB2w5X+msJcYbvXDGab186ksR4l4PdiYjIqQo7KBljfg6sAd6y1rZ5kKy1tr6tuoiISLSrbXDzs5e28V//3IMn4F/Bs4f2YsWcfMYMyHCuOREROW2RmFH6d7y32B0zxjyFNzS9bK1tiMBni4iIdFvrd5WyYHUh+4/X+GspCXH86Jqx3HrxMOJcOjhWRCRaRSIoDQZmATcCNwO3AdXGmOeAtcBz1trKdt4vIiISVcprG1n+/Bb++n7wJq6fG9WPh2bnMaRPqkOdiYhIpIQdlKy1R4DHgceNMZnAdXhD03Tgi0C9MeYVvKHpKWvt0XDHFBERccpLmw+z6MlCDle03FWemRzPousm8IXJuRijWSQRkVgQVlAyxvwIb/jZBmCtrQD+B/gf3wYO1+ANTdfhDU7/aYxZj/f2vCettfvCGV9ERKSrlFbVs/ipTTxTUBJUv2biAJbeMInszGSHOhMRkc4Q7ozSCiAFuB/AGDMG2GG96oGngKeMMS7gUmA2cAPwMPAzY8wn1trJYfYgIiLSaay1rP34APc/s5mymkZ/vV96EktvmMi1eYMc7E5ERDpLuEGpsdVnbAHuA5YFXmSt9QCv+R7fMcacizc0zQpzfBERkU5zoKyWu9cU8sb24LvG507OZdGM8fRKTXSoMxER6WzhBqWDwNkB3xugw4MirLUfAh8Cd4c5voiISMR5PJY/v7ePFc9vpbrB7a/n9Erhodl5XDKmv4PdiYhIVwg3KD0N/Lsx5lm8644A2jxLSUREJBrsOlrFgtUFfLD3hL9mDNx60TB+dM1Y0pI67ax2ERHpRsL9234hMByYAUzz1RYZY+YAH/keHwOfWGurwxxLRESk0zS6Pfz2rd38/OUdNDR5/PWR/dNYOTefyWf0cbA7ERHpamEFJd/5SDONMWOBK4FfAmXAaCAf+CreGSZrjNmBNzR9BHxkrX01nLG7I2PMJcCdwGS850vdZq19otU1i/Gu4wp02Fo7sCt6FBGRz9p4oJx5qwrYXFLhr8W7DLdfNpI7rhhFUnycg92JiIgTInL/gG978G3GmF8Cv8K7mcME4JyAx5nAWOAmvOEpFv/VSQc2An/0PULZBlwW8L07xHUiItKJ6hrd/OKVHTz+5m7cnpY7x/Nyslg5N5/xgzId7E5ERJwU6Rut5wIlvl3uNtISGgD/9uGTCd4AImZYa58DngMwxjzRzqVN1tpDXdKUiIi06f09x1mwuoDdpS13hifFu/jh1DF8bcpw4uM63JtIRERiWLgHzq4CNtByO92a9q631m4HtgN/DWfcgPHn4j2f6Sy8M1YZwF+stTe3855cvOc+TQP6AiXAk8ASa+2JUO+LsBHGmANAA/AecLe1dncXjS0i0qNV1jWy8h/b+NO7wWeeXzC8Dyvm5DOsX5pDnYmISHcS7ozSbN/DAvh++G/exGED3vBUEvrtYVuENyBVAcXAuPYuNsaMBNYD2cA6YCtwPvA9YJoxZoq19lgn9gveYPRV39jZeP8M640xE7tgbBGRHu21rUdYuLaQg+V1/lp6Ujx3TR/Hl84bistlHOxORES6k3CD0hm0rEGa7Hu+HpjZfIEx5jCfDU9FYY7b7Ad4A9JOvDNLr3Vw/WN4w8l3rbWPBvT4M99nPQB8O6C+DO/Ofu253Fr7+sk2bK19PvB7Y8y7wG7gVuBnJ/s5IiJy8o5XN7D0mc2s/fhAUP3Kcdksu3ESg7JSHOpMRES6q3B3vSsCivDOzgBgjHkc+AawBTgGjACmA9cGXFNqrR0Qzti+8f3ByJj2fwtojBkBTAX24t1wItB9wLeAW4wxPwzYyvznwJ87aGP/KbT8GdbaKmPMJrw7BYqISARZa3mmoITFT23iWHWDv94nLZHF109kZv6gDv/9EBGRnimimzkYY34EfBm4xFr7z4D6xcASvFuIF+HMobRX+J5f9G024WetrTTGvI03SF0IvOKrlwKlndmUMSYZ7y2DHc2GYYzZEOKldm85FBHpiQ6V17HoyY28vOVwUH3WWYO5d+ZE+qQlOtSZiIhEg0jvencH8NfAkARgrV0PXG2MmQ/cC0yJ8LgnY6zveXuI13fgDUpj8AWlU2WMSQdG+b51AUONMWcBx621+33X/AR4Gu9MVDZwD5AG/OF0xhQRkWDWWv72QREPPruFyvomf31QVjIP3DiJK8aFfUODiIj0AJEOStl4b7drk7V2hTFmNnA38MUIj92RLN9zeYjXm+u9whjjXIJnhpb4Hn/Au4EDQC7eXf/6AUeBd4ELrbXB2y+1wVo7ua26b6bpnNPuWkQkRuwtreauNYW8szv4n6KbLxzK/GnjyEhOcKgzERGJNpEOSjuAqzq45nXgKxEeNxKab1I/7dsCfZs6tHuzu7X2ptP9fBERaVuT28N/v72Xn760jbrGlrurh/dLY/nsPC4Y0dfB7kREJBpFOij9HnjYGHO3tfbBENcMBPpEeNyT0TxjlBXi9cxW14mISBTYeqiC+asK+LS45a/vOJfhm58fwfevGk1yQpyD3YmISLSKdFD6JTADWGqMuQq4t9WmDjOAm/Bu6NDVtvmex4R4vXnXuVBrmEREpBupb3Lzq1d38tjru2jytNwMMH5QJivn5JOXG+r3YiIiIh2LaFCy1rp9Yegx4OvAG8aYY7RsXJCD99a0X0Ry3JPUvHZoqjHGFbjznTEmA+8GE7V41wyJiEg3tmHfCeavLmDnkSp/LTHOxfeuGs23LhlBQpzLwe5ERCQWRHpGCWttI/BNY8xvge8AV+PdaMCNd1bnIWvtnyI97kn0tcsY8yLene3uAB4NeHkJ3p3nfhNwhpKIiHQz1fVN/OTFbTyxfi82YEXp5DN6s2JOPqOy051rTkREYkrEg1Iza+37wC0AxpgkwG2tbWr/XafGGDMLmOX7dqDv+SJjzBO+r0uttXcGvOV2YD3wiDHmSryH4l4AXI73lruFkexPREQi560dR7lrTSHFJ2r9tdTEOOZPG8ctF56By6WDY0VEJHI6LSgFstbWd9JHnwXc2qo2wvcA2Af4g5JvVulc4H5gGjAdKAEeAZZYa493Up8iInKaymsaWfbsZv6+oTiofsmY/jx44yRye6c61JmIiMSyLglKncVauxhYfIrvKQJu64x+REQksv6xsYR71m3iaGXL79t6pSZwz4wJzD4nB2M0iyQiIp0jqoOSiIjEpiOVddy3bhPPbzwUVJ+RP4jFMyfSPyPJoc5ERKSnUFASEZFuw1rLqg3FLHt2C+W1jf56dkYSS2dN4pqJA9t5t4iISOQoKImISLdQdLyGu9cW8taO0qD6TecN4a7p48lKSXCoMxER6YkUlERExFFuj+WP7+zlxy9so6bB7a8P6ZPC8tn5TBnVz7nmRESkx1JQEhERx+w4XMn81QV8tL/MX3MZ+NqU4fzH1DGkJuqfKRERcYb+BRIRkS7X0OThN2/s4tFXd9Lg9vjrYwaks2JOPmcP7e1gdyIiIgpKIiLSxQqKy5i3qoCthyr9tYQ4wx2Xj+L2y0aRGO9ysDsREREvBSUREekStQ1ufv7ydn771m48tqV+5pBerJyTz9iBGc41JyIi0oqCkoiIdLp3dh3jrjUF7D1W468lJ7i4c+pYbpsynDiXDo4VEZHuRUFJREQ6TUVdI8uf38r/vLc/qH7xyL4sn53P0L6pDnUmIiLSPgUlERHpFC9vPsyiJzdyqKLOX8tIjmfRjPF88dwhGKNZJBER6b4UlEREJKKOVdWz5OnNPPXpwaD61AkDWDprEgMykx3qTERE5OQpKImISERYa1n3yUGWPL2JEzWN/nq/9ESWXD+J6XkDNYskIiJRQ0FJRETCdrCsloVrC3lt29H/z959h0dVpm8c/77pAUKJdAi9l8AGBQWkC9gQgV0si8vuqruWBUQXUFFDUcBdFTuu64K72GlWRKX8FFFQEELvkV4DIYQkJJn390cmQyaVQJKZydyf68o15jnvOeeZ+Cbkzmlu9SEx9XjixjZUqxjioc5EREQujYKSnWIk/gAAIABJREFUiIhcMofD8s6afcxYvI2zaRmuer2q4Tx9azt6tazpwe5EREQunYKSiIhckj3HzzJh/kbWxCe4asbAXVc35O8DW1EpVP/EiIiI79K/YiIiUiwZmQ7e/G4vL3yzg/MZDle9SY2KPDs0misbRXqwOxERkZKhoCQiIhdt86FExs+PY9PBM65aYIDhvp5NebBPM8KCAz3YnYiISMlRUBIRkSKlpmfy8rKdzPq/PWQ6rKverl5lZgyNpm3dKh7sTkREpOQpKImISKF+jk9g3Pw49hxPdtVCgwJ46LoW3N29MUGBAR7sTkREpHQoKImISL7OpmXwjy+38d8ff8VeOIhE58aRTB/SniY1KnmuORERkVKmoCQiInms2H6Mxxdu4uDpFFetUmgQE65vxR2dGxAQoAfHiohI+aagJCIiLqeSzzPl8y0sWHfQrd67ZQ2evrU9dauGe6gzERGRsqWgJCIiWGv5YuMRnvpkEyfOnnfVq1UIJnZQWwZ1qIsxOookIiL+Q0FJRMTPHT2TysRFm/h6y1G3+qAOdXnq5jZcUSnUQ52JiIh4joKSiIifstby4c/7mfr5VpJSM1z12pXDmDq4Hf3a1PJgdyIiIp6loCQi4of2nTzHhAVxrNp90q1+R5cGTLi+FZXDgj3UmYiIiHdQUBIR8SOZDsvs7/fyz6+2k5rucNUbXVGBaUOiuabpFR7sTkRExHsoKImI+IntR5IYNz+ODftPu2oBBu65tglj+rUgPCTQg92JiIh4FwUlEZFy7nyGg1eX7+K1FbtIz7zw5NhWtSN4dlg00fWrerA7ERER76SgJCJSjv2y7xTj58ex4+hZVy0kMIC/9WnGX3s1JTgwwIPdiYiIeC8FJRGRcujc+Qye+2oH//l+L/bCQSRiGlRlxtBomteK8FxzIiIiPkBBSUSknPl+1wkmLIhjf0KKqxYeHMi4gS2565pGBAbowbEiIiJFUVASESknElPSeebzrXzw8363+rXNq/PMre2Jiqzgoc5ERER8j4KSiEg5sGTzEZ5YtIljSWmuWpXwYJ64qQ1DY+phjI4iiYiIFIeCkoiIDzuelEbsJ5v5fONht/oN7WsTO6gtNSPCPNSZiIiIb1NQEhHxQdZaFqw7yOTPtpCYku6q14gIZcotbRnYro4HuxMREfF9CkoiIj7mwKlzPLZwE9/uOO5W/92V9Xn8hjZUqRDsoc5ERETKDwUlEREf4XBY/vfjr8z4chvnzme66vWrhTN9SDTdm1f3YHciIiLli4KSiIgP2HXsLBPmx/Hzr6dcNWPgj10b88iAFlQI0Y9zERGRkqR/WUuQMaYH8AjQCagL/NFaOyfXmHigYT6rf2GtvbG0exQR35Ke6eBf3+7hxW92cj7T4ao3r1mJ6UOj6dSwmge7ExERKb8UlEpWJWAT8F/nR36uAgJzfF4HWAt8WLqtiYiv2XggkXHz49h6+IyrFhRguL93Mx7o3ZTQoMBC1hYREZHLoaBUgqy1XwBfABhj5hQwxu3qa2PMn4EzwEel3Z+I+IbU9Exe+GYH//5uL5kO66pH16/CjKHRtK5T2YPdiYiI+AefDkrGmGFAT6Aj0AGIAN6x1v6+kHXqA5OBgcAVwGFgETDJWnuqoPVKg8l6AuSfgbnW2nNluW8R8U6r95xkwoKN7D2R7KqFBQfw8HUt+WO3RgQFBniwOxEREf/h00EJmEhWQDoLHABaFTbYGNMUWAXUBD4GtgGdgdHAQGNMN2vtyVLt2N11QGPg32W4TxHxQkmp6cz4chtzf9znVr+6SSTTh0TTqHpFD3UmIiLin3w9KD1EVkDaRdaRpeVFjH+NrJA0ylr7cnbRGPO8c1tPA3/NUZ8KPF7ENntba1cUu/Ms9wA/WWvXX+L6IlIOLNt2lMcXbuJwYqqrFhEaxGM3tua2q6LIOvgsIiIiZcmng5K11hWMivpFwhjTBOgPxAOv5lr8FHAvMMIY87C1Nvucl5nA3CLa2FfE8oL6qQncAjxwKeuLiO87eTaNyZ9t4eP1h9zq/VrXYurgdtSuEuahzkRERMSng1Ix9XG+fmWtdeRcYK1NMsZ8T1aQuhpY6qyfAE6UUj8jgTTg/YtdwRiztoBFhZ5yKCLexVrLJxsOMenTLSQkn3fVr6gYwqRb2nJj+zo6iiQiIuJh/hSUWjpfdxSwfCdZQakFzqBUXMaYSkAz56cBQANjTEcgwVq7L8c4A9wNvG+tTbqUfYmIbzqcmMLEhZtYuu2YW33Ib+rxxE1tqFYxxEOdiYiISE7+FJSqOF8TC1ieXa96Gfu4EvfrpCY5P94m6whStl5Ac6DAu/Plx1rbKb+680hTTHG2JSJly+GwvPfTPqZ9sY2zaRmuet0qYTw9pD29W9b0YHciIiKSmz8FpaJkn+diCx1VCOdNHYo8X8Z5bZXOqxHxE3tPJDNhfhyr9ya41Udc3ZBxA1sSERbsoc5ERESkIP4UlLKPGFUpYHnlXONERC5LRqaDt1bu5fmvd5CWceHSyCbVKzJ9aDSdG0d6sDsREREpjD8Fpe3O1xYFLG/ufC3oGiYRkYu25dAZxs+PY+PBC397CQww/KVHE0b1bU5YcKAHuxMREZGi+FNQyr52qL8xJiDnne+MMRFANyAF+NETzYlI+ZCWkckry3bx+ordZDgunMnbpk5lnh0WTbt6BR3UFhEREW/iN0HJWrvbGPMVWXe2ewB4OcfiSUBF4I0cz1ASESmWtb8mMH7+RnYdO+uqhQQFMKZfc+65tgnBgQEe7E5ERESKw6eDkjFmMDDY+Wlt5+s1xpg5zv8+Ya19JMcq9wOrgJeMMX2BrUAXoDdZp9w9XupNi0i5k5yWwT+WbOftH+KxOW4Hc1WjakwfGk3TGpU81puIiIhcGp8OSkBH4A+5ak2cHwC/Aq6g5DyqdCUwGRgI3AAcBl4CJllrExARKYZvdxzn0QUbOXg6xVWrGBLIhOtbcWeXhgQE6AaXIiIivsing5K1NhaILeY6+4E/lkY/IuI/Tp87z9TPtzJv7QG3eq+WNXj61vbUqxruoc5ERESkJPh0UBIR8YTFGw/zxMebOXE2zVWrWiGYp25uw+CO9TBGR5FERER8nYKSiMhFOnYmlSc+3sSSzUfd6jdF1yF2UFuqVwr1UGciIiJS0hSURESKYK3lo58PMPXzLZxJzXDVa1UOZcot7ejftnYha4uIiIgvUlASESnE/oRzPLpgIyt3nXCr3945ignXt6ZKeLCHOhMREZHSpKAkIpKPTIfl7VXx/GPJdlLSM131BpEVmD6kPV2bVfdgdyIiIlLaFJRERHLZeTSJcfPj+GXfaVctwMCfuzdm7HUtCQ8J9GB3IiIiUhYUlEREnM5nOHh9xW5eWb6T9MwLT45tVTuCGUOj6RBV1YPdiYiISFlSUBIRATbsP824eXFsP5rkqgUHGv7Wpzl/7dmUkKAAD3YnIiIiZU1BSUT8Wsr5TJ7/ejtvrdyL48JBJH7ToCozhkbTolaE55oTERERj1FQEhG/tWr3CSbM38i+hHOuWnhwIH8f0JI/dG1EYIAeHCsiIuKvFJRExO8kpqQzffFW3luz363evVl1pg1pT1RkBQ91JiIiIt5CQUlE/MrXW44ycdFGjp5Jc9UqhwUx8aY2/LZTfYzRUSQRERFRUBIRP3HibBqxn2zms7jDbvUBbWsx5ZZ21Kwc5qHORERExBspKIlIuWatZdH6g0z6dAunz6W76tUrhTLllrZc376OB7sTERERb6WgJCLl1sHTKTy+cCMrth93qw/rVJ+JN7amaoUQD3UmIiIi3k5BSUTKHYfD8s7qX5m+eBvJ5zNd9XpVw5k2pD09WtTwYHciIiLiCxSURKRc2X38LBPmx/FT/ClXzRj4wzWN+PuAllQM1Y89ERERKZp+YxCRciE908Gb3+1h5jc7OZ/hcNWb1qjIs8Oi6dQw0oPdiYiIiK9RUBIRn7fpYCLj5sWx5fAZVy0owHBfr6Y82KcZoUGBHuxOREREfJGCkoj4rNT0TF5cupN/fbuHTId11dvXq8KModG0qVvZg92JiIiIL1NQEhGftGZvAhPmx7HnRLKrFhoUwNjrWvDn7o0JCgzwYHciIiLi6xSURMSnnE3LYMbibfzvx1/d6l0aRzJ9aDSNq1f0UGciIiJSnigoiYjPWL79GI8v2MihxFRXrVJoEI/e0Irbr2pAQIDxYHciIiJSnigoiYjXS0g+z5TPtrDwl4Nu9b6tajL11nbUqRLuoc5ERESkvFJQEhGvZa3ls7jDxH6ymZPJ5131yIohxA5qy83RdTBGR5FERESk5CkoiYhXOpKYysRFm/hm61G3+uCOdXny5rZEVgzxUGciIiLiDxSURMSrWGt5/6f9PPP5VpLSMlz1OlXCePrWdvRpVcuD3YmIiIi/UFASEa8RfyKZRxds5Ic9J93qv7+6AeMHtiIiLNhDnYmIiIi/UVASEY/LyHQw+/t4nvt6O6npDle9cfWKTBvSnqubXOHB7kRERMQfKSiJiEdtO3KG8fPi2HAg0VULDDDcc20TxvRrTlhwoAe7ExEREX+loCQiHpGWkcmry3bx2ordZDisq966TmWeHRpN+/pVPNidiIiI+DsFJREpc+v2nWL8vDh2HjvrqoUEBjC6X3Pu7dGE4MAAD3YnIiIioqAkImXo3PkM/rlkB7NX7cVeOIhEp4bVmDE0mmY1K3muOREREZEcFJREpEys3HmCCQviOHAqxVWrEBLI+IGtGHF1QwIC9OBYERER8R4KSiJSqhLPpfP0F1v48OcDbvUeLWrwzK3tqF+tgoc6ExERESmYgpKIlJovNx3miY83czwpzVWrWiGYJ25sw5CYehijo0giIiLinRSURKTEHUtK5amPN7N40xG3+o3t6xA7qC01IkI91JmIiIjIxVFQEpESY61l3toDTP18K4kp6a56jYhQptzSjoHtanuwOxEREZGLp6AkIiVif8I5Hlu4ke92nnCrD78yisduaE2VCsEe6kxERESk+PSwkhJijOlhjPnEGHPQGGONMSPzGRNhjJlpjPnVGJNijFlljLnKA+2KlJhMh2X293sZMPNbt5AUFRnOO3d3YcawaIUkERER8Tk6olRyKgGbgP86P/LzbyAa+ANwAPg98I0xpo219mCZdClSgnYdS2LcvDjW7TvtqgUY+FO3xozt34IKIfoRIyIiIr5Jv8WUEGvtF8AXAMaYObmXG2PCgaHAUGvtCmc51hhzM3AfMLFsOhW5fOmZDmat2M3Ly3ZxPtPhqreoVYkZQ6P5TYNqHuxORERE5PL5bFAyxgwDegIdgQ5ABPCOtfb3haxTH5gMDASuAA4Di4BJ1tpTpdxyEBAIpOaqpwDdS3nfIiUm7sBpxs2LY9uRJFctONDwQO9m3N+rGSFBOqNXREREfJ/PBiWyjsB0AM6SdRpbq8IGG2OaAquAmsDHwDagMzAaGGiM6WatPVlazVprk4wxPwATjTGbgCPA7cA1wK7S2q9ISUk5n8nMb3bw5nd7cNgL9Q5RVXl2aDQta0d4rjkRERGREubLf/p9CGgBVCbr1LWivEZWSBplrR1srZ1gre0DvAC0BJ7OOdgYM9V5U4bCPnoVs+cRgIOsYJcGjALeAzKLuR2RMvXD7pNc/+K3vPHthZAUFhzAxBtbs+C+rgpJIiIiUu747BEla+3y7P82xhQ61hjTBOgPxAOv5lr8FHAvMMIY87C1NtlZnwnMLaKNfcVoGWvtbqCnMaYiUNlae9gY8wGwtzjbESkrZ1LTmb54G++udp/qXZtewfQh0TS4ooKHOhMREREpXT4blIqpj/P1K2utI+cC5ylx35MVpK4GljrrJwD3B8KUEGcYSzbGVAMGAOMuZj1jzNoCFhV62qHIpfhmy1EmLtrEkTMXLquLCAti4o2t+d2VUUX+gUJERETEl/lLUGrpfN1RwPKdZAWlFjiDUnEZYyoBzZyfBgANjDEdgQRr7T7nmAHOZducY/8BbAdmX8o+RUrDybNpTPp0C59sOORWv65NLaYObketymEe6kxERESk7PhLUKrifE0sYHl2vepl7ONKYHmOzyc5P94GRuboYxpQH0gA5gOPW2vTL2YH1tpO+dWdR5piLqlrESdrLR+vP8SkTzdz6tyFKVm9UgiTBrXjhva1dRRJRERE/Ia/BKWiZP/2ZwsdVQjns5EK/S3SWvsh8OGl7kOktBw6ncLERZtYtu2YW31ITD2euLEN1SqGeKgzEREREc/wl6CUfcSoSgHLK+caJ+IXHA7Lu2v2MX3xNs6mZbjq9aqG8/St7ejVsqYHuxMRERHxHH8JStudry0KWN7c+VrQNUwi5c6e42eZsGAja/YmuGrGwF1XN+TvA1tRKdRffjyIiIiI5OUvvwllXzvU3xgTkPPOd8aYCKAbkAL86InmRMpSRqaDf6/cywtf7yAt48JNIJvUqMizQ6O5slGkB7sTERER8Q5+EZSstbuNMV+RdWe7B4CXcyyeBFQE3sjxDCWRcmnzoUTGz49j08EzrlpggOG+nk15sE8zwoIDPdidiIiIiPfw2aBkjBkMDHZ+Wtv5eo0xZo7zv09Yax/Jscr9wCrgJWNMX2Ar0AXoTdYpd4+XetMiHpKansnLy3Yy6//2kOm4cM+StnUr8+ywaNrWLejyPRERERH/5LNBCegI/CFXrYnzA+BXwBWUnEeVrgQmAwOBG4DDwEvAJGttAiLl0M/xCYybH8ee4xcOmIYEBfBQvxbcc21jggIDPNidiIiIiHfy2aBkrY0FYou5zn7gj6XRj4i3OZuWwT++3MZ/f/wVm+PG950bRTJ9aHua1KjkueZEREREvJzPBiURKdiK7cd4fOEmDp5OcdUqhQYx/vpW3Nm5AQEBenCsiIiISGEUlETKkVPJ55ny+RYWrDvoVu/dsgZP39qeulXDPdSZiIiIiG9RUBIpB6y1fLHxCE99sokTZ8+76tUqBBM7qC2DOtTFGB1FEhEREblYCkoiPu7omVSeWLSJr7YcdasP6lCXp25uwxWVQj3UmYiIiIjvUlAS8VHWWj78eT9TP99KUmqGq167chhTB7ejX5taHuxORERExLcpKIn4oH0nzzFhQRyrdp90q9/RpQETrm9F5bBgD3UmIiIiUj4oKIn4kEyHZfb3e/nnV9tJTXe46g2vqMD0IdFc0/QKD3YnIiIiUn4oKIn4iO1Hkhg3P44N+0+7agEG7rm2CWP6tSA8JNCD3YmIiIiULwpKIl7ufIaDV5fv4rUVu0jPvPDk2Fa1I3h2WDTR9at6sDsRERGR8klBScSL/bLvFOPnx7Hj6FlXLSQwgL/1acZfejYlJCjAg92JiIiIlF8KSiJe6Nz5DJ77agf/+X4v9sJBJGIaVGXG0Gia14rwXHMiIiIifkBBScTLrNp1ggkLNrIv4ZyrFh4cyLiBLbnrmkYEBujBsSIiIiKlTUFJxEskpqQz7YutvP/Tfrf6tc2r88yt7YmKrOChzkRERET8j4KSiBdYsvkITyzaxLGkNFetSngwT9zUhqEx9TBGR5FEREREypKCkogHHU9KI/aTzXy+8bBb/Yb2tYkd1JaaEWEe6kxERDzF4XCQkJBAUlISaWlp2JwXq4r4MWMMoaGhREREEBkZSUBA6d7USkFJxAOstSxYd5DJn20hMSXdVa9eKZSpg9sysF0dD3YnIiKe4nA42L9/P+fOnSt6sIifsdaSmppKamoqycnJREVFlWpYUlASKWMHTp3jsYWb+HbHcbf6bzvVZ+KNbahSIdhDnYmIiKclJCRw7tw5goKCqF27NhUrViz1v5qL+AqHw0FycjJHjhzh3LlzJCQkUL169VLbn4KSSBlxOCz/+/FXZny5jXPnM131+tXCmTakPdc2r+HB7kRExBskJSUBULt2bSIi9CgIkZwCAgJc3xcHDhwgKSlJQUnE1+06dpYJ8+P4+ddTrpox8MeujXm4fwsqhupbUUREIC0t66Y+FStW9HAnIt4r+/sj+/ultOi3M5FSlJ7p4F/f7uHFb3ZyPtPhqjevWYnpQ6Pp1LCaB7sTERFvk33jBp1uJ1Kw7LsBl/aNThSURErJpoOJ/H1eHFsPn3HVggIM9/duxgO9mxIaFOjB7kRERER8U1k9NkVBSaSEpaZnMvObnbz53R4yHRf+0hFdvwozhkbTuk5lD3YnIiIiIhdDQUmkBK3ec5IJCzay90SyqxYWHMDD17Xkj90aERSoUylEREREfIGCkkgJSEpNZ8aX25j74z63+tVNIpk+JJpG1XVRroiIiIgv0Z+3RS7Tsm1H6f/Ct24hKSI0iGlD2vPePVcrJImIiEipWLFiBcYYYmNjS31fxhh69epV6vvxJgpKIpfo5Nk0Rr//C3+a8zOHE1Nd9X6ta/L12J7c3rlBmV1sKCIiUl7t2LGDsWPHEhMTQ2RkJMHBwURGRtKlSxceeeQR1q5d6zb+5ptvxhjDl19+me/2WrRogTGGkSNH5rv8ySefxBjDlClTiuxtzpw5GGMK/WjUqFFx37J4CZ16J1JM1lo+2XCISZ9uISH5vKt+RcUQYge15aboOgpIIiIil8lay+TJk5k8eTIOh4OYmBiGDx9OZGQkSUlJxMXF8fLLL/Pcc8/xyiuv8MADDwDQt29fPvvsM5YuXcrAgQPdtnnw4EF27tyJMYalS5fmu99ly5a5tnOxOnTowODBg/NdVrVq1YveTnF17tyZrVu3lupDV7Nt3bqVChUqlPp+vImCkkgxHE5MYeLCTSzddsytfutv6vHETW2IrBjioc5ERETKl8mTJxMbG0tUVBTvvfce3bp1yzPm2LFjzJw5k8TERFetT58+wIXAk1N2OBo2bBgfffQRO3fupHnz5q7lycnJrFmzhkqVKtG5c+eL7rVjx45lcvpbbhUqVKBVq1Zlsq+y2o830al3IhfB4bC8s/pXrnv+W7eQVLdKGLNHXsULwzsqJImIiJSQPXv2MHXqVEJCQli8eHG+IQmgZs2aPPPMM4wbN85Va9++PTVq1GD9+vUkJCS4jV+2bBkVKlRg/Pjxrs9z+u6770hPT6dnz54EBZXO8YRGjRrRqFEjzp49y0MPPURUVBTh4eF07NiRRYsWAZCRkcEzzzxD8+bNCQsLo2nTprzyyit5tlXQNUp79uzh3nvvpVmzZoSHhxMZGUn79u3561//ysmTJ13jzp8/z0svvURMTAzVqlWjQoUKNGrUiFtuuYVvvvnGbZsFXaOUmJjIo48+SsuWLQkLC6NatWoMGDAgz/q5+12/fj033ngjVatWpUKFCvTs2ZNVq1Zdwle09CgoiRQh/kQyd/z7Rx5fuImzaRmu+oirG7LkoR70blXTg92JiIiUP7NnzyYjI4Nhw4bRtm3bIsfnDDXGGHr37o3D4WDFihVu45YtW0b37t2JiYmhRo0aeU6/y/68OKfdXYr09HSuu+46vvjiC2655RZGjBjB7t27GTp0KEuXLmX48OG89tpr9OrVi7vvvpuzZ8/yt7/9jQ8++KDIbR8+fJirrrqK2bNn07ZtW0aNGsWIESNo3Lgx//vf/zh8+LBr7MiRIxk9ejTp6encddddjBo1ih49erBx48YCr/HK6fTp03Tt2pXp06dTpUoVxowZw9ChQ/nhhx/o378/b7zxRr7r/fzzz3Tt2pXU1FTuvvtubrrpJlauXEnfvn3Zvn37xX8hS5lOvRMpQEamg7dW7uX5r3eQluFw1ZtUr8j0odF0bhzpwe5ERMQfNZrwuadbuGjx02+85HW///574MJpdMXVt29fPvzwQ5YtW8aQIUMA2LlzJ/v37+f+++/HGEPPnj1ZsWIF1lrXtcWXcn0SwPr16ws89e7qq6/Oc63UoUOHiImJYcWKFYSGhgIwYsQIevTowW9/+1uaNm3Kpk2bXNc3jR07llatWjF9+nSGDx9eaC/z5s0jISGBmTNnMnr0aLdlycnJBARkHSdJTEzk/fffp1OnTqxevZrAwEC3sTmPPBVk/PjxbNmyhXvvvZdZs2a5vo7jx4/nyiuvZNSoUQwYMCDPDS0+//xzZs+e7XZDjTfeeIO//vWvvPjii7z22mtF7rssKCiJ5GPLoTOMnx/HxoMXznkODDD8pUcTRvVtTlhwYCFri4iIyOU4cuQIAPXq1cuzLD4+njlz5rjVqlatypgxY1yf53edUvZ/Zy/r3bs38+bNY+PGjURHR3Pq1CnWr19PzZo1ad++fbH63bBhAxs2bMh32ejRo/MEJYCZM2e6QhLAtddeS+PGjdm7dy8zZsxwuwlEkyZN6NatG9999x2ZmZl5Qk1+wsPD89QqVrzwyBJjDNZaQkNDXeEppyuuuKLQ7aenpzN37lwqVarEtGnT3G5k1bx5c0aNGsXUqVP573//y5NPPum2brdu3fLcdfBPf/oTDz74IGvWrCnyvZUVBSWRHNIyMnll2S5eX7GbDId11dvUqcyzw6JpV6+KB7sTERHxD9Zm/Ruc311k4+PjmTRpklutYcOGbkGpWbNmNGjQgK1bt3L48GHq1KnDsmXLqFy5Mp06dQKyghJkBajo6GhWrFiBw+GgT58+xb577R/+8Ic84a0wVatWpWnTpnnqdevWZe/eva4ec6pXrx6ZmZkcOXIk3wCZbdCgQTz22GM88MADLFmyhAEDBtCtWzfatGnj9r4qV67MzTffzKeffkrHjh0ZOnQo1157LV26dLmou9tt27aNc+fO0a1bNyIj855l06dPH6ZOncovv/ySZ9mVV16ZpxYcHEytWrU4depUkfsuKwpKIk5rf01g/PyN7Dp21lULCQpgTL/m3HNtE4IDdUmfiIh41uWczuZL6tSpw7Zt2zh48GCeZb169XIFqYyMDIKDg/PdRt++fZk9ezbLli3jjjvuYPny5fTo0cN1NKZ169bUqlWLpUuXMmbMmDK7PgmgSpX8//Cafa1Vfsuzl6Wnpxe67YYNbCo2AAAgAElEQVQNG7JmzRpiY2P58ssvWbBgAQBRUVE88sgjjBo1yjX2gw8+YMaMGbz77rs89dRTAISFhTFs2DD++c9/UqtWrQL3k32nwTp16uS7PLt++vTpPMsKumV6UFAQmZmZhb6/sqTf/MTvJadlEPvJZobN+sEtJF3VqBqLR1/L/b2aKSSJiIiUoey73BX0rKOLkfP0u40bN3L8+HHXUaRsvXr14ttvvyUzM/OSr0/yRq1bt+aDDz7g5MmT/Pzzz0yfPh2Hw8Ho0aN56623XOPCw8OJjY1lx44d7Nu3j7lz59K9e3fmzp3LsGHDCt1HdpjLPk0yt+ybRhQUCn2BfvsTv/btjuP0f+Fb5qyKx/nHKSqGBDLllrZ8cO81NK1RybMNioiI+KGRI0cSFBTEvHnz2Lp16yVtIzvwLFu2LM/1Sdl69+7NmTNn+PTTT9m6dSuNGzemcePGl9e8FwkKCqJTp06MHz+e9957D8B1C/LcoqKiuPPOO1myZAnNmzdn5cqVhd7QoWXLllSoUIH169fne7rc8uXLAYiJiSmBd+IZCkril06fO88jH23grv+s4eDpFFe9Z4safDW2JyOuaURAQPHOTxYREZGS0bRpUyZOnMj58+e5/vrrC3y+Tn6ndWWrU6cOrVq1Ij4+nv/85z9ERkYSHR3tNib7CFP2zQbKw9GkNWvWcPTo0Tz17Fr29UfHjx9n9erVecYlJyeTlJREUFAQISEFPyMyJCSEO++8k7Nnz+a5WcPu3bt56aWXCA4OZsSIEZfzdjxK1yiJ31m88TBPfLyZE2fTXLWqFYJ56uY2DO5Yr9gXcIqIiEjJe/LJJ7HWMmXKFLp160anTp3o3LkzkZGRnD59mvj4eNdDTXv06JHvNvr27cu2bdvYuHEjQ4YMyXN3txYtWlC3bl02btzoGn8pCrs9OFDospL27rvv8uqrr9KzZ0+aNWtGtWrV2L17N59++imhoaGum14cPHiQq6++mtatWxMTE0NUVBRnzpzhs88+48iRI4waNYqIiIhC9zV9+nS+++47XnnlFX766Sd69+7NiRMn+PDDD0lKSuKVV17x6SN0CkriN46dSeXJjzfz5Wb3c2lviq5D7KC2VK8UWsCaIiIiUtaMMcTGxnL77bcza9Ysli9fzrvvvktycjIRERE0bdqU++67jxEjRhR4elffvn159dVXAfJcn5Std+/evPPOOxhjLvm5TYXdHhzKNijdfvvtpKWlsWrVKtatW0dKSgr16tXjtttu4+GHH6Zdu3YANGrUiEmTJrFixQqWL1/OiRMniIyMpGXLlkyfPp3bbrutyH1FRkbyww8/MG3aNBYsWMDzzz9PeHg4nTt35u9//zv9+/cv7bdbqkz2XUPEdxlj1sbExMSsXbvW0614JWstH609wNTPtnAmNcNVr1U5lCm3tKN/29oe7E5EROSC7OtxWrdu7eFORLzbxX6vdOrUiXXr1q2z1ua953oRdI1SCTLGPGqM+ckYc8YYc9wY86kxpl2uMT2MMZ8YYw4aY6wxZqSH2vUL+xPOMeKtNYybF+cWkm7vHMVXD/VUSBIRERGRfOnUu5LVC3gN+AkwwGTgG2NMG2ttgnNMJWAT8F/nh5SCTIfl7VXx/GPJdlLSL9yPv0FkBaYPaU/XZtU92J2IiIiIeDsFpRJkrR2Q83NjzAggEegGfOoc8wXwhXP5nDJu0S/sPJrEuPlx/LLvwp1wAgz8uXtjxl7XkvCQQA92JyIiIiK+wKeDkjFmGNAT6Ah0ACKAd6y1vy9knfpkHekZCFwBHAYWAZOstXlvAn95Isg6vbGktyv5OJ/h4PUVu3ll+U7SMy9ce9eyVgQzhkXTMSr/p0CLiIiIiOTm00EJmEhWQDoLHABaFTbYGNMUWAXUBD4GtgGdgdHAQGNMN2ttwU/WKr4XgfXADyW4TcnHhv2nGTcvju1Hk1y14EDDg72bc1+vpoQE6XI8EREREbl4vh6UHiIrIO0i68jS8iLGv0ZWSBplrX05u2iMed65raeBv+aoTwUeL2Kbva21K3IXndvsDnS31mbmWUtKRMr5TJ7/ejtvrdyLI8cNHDtGVeXZYdG0qFX4/f9FRERERPLj00HJWusKRkU9JNQY0wToD8QDr+Za/BRwLzDCGPOwtTbZWZ8JzC2ijX357OsF4DayQtSeItaXS7Rq9wkmzN/IvoRzrlp4cCCPDGjJyK6NCAzQg2NFRERE5NL4dFAqpuwniH1lrXXkXGCtTTLGfE9WkLoaWOqsnwBOFGcnxpgXyQpJvay12y67a/dtF/SgpEJPOSxvzqSmM+2Lbby3xj2jdm9WnWlD2hMVWcFDnYmIiIhIeeFPQaml83VHAct3khWUWuAMSsVljHkVGAEMBk4ZY7If0nPWWnvWOaYS0MxZDwAaGGM6AgnW2jxHp8Td11uOMnHRRo6eSXPVKocFMfGmNvy2U/0ijyyKiIiIiFwMfwpKVZyviQUsz65fzq3R7ne+5g5ak4BY539fifu1VJOcH28DIwvbeEFPFHYeaYopXqu+5cTZNGI/2cxncYfd6gPa1mLKLe2oWTnMQ52JiIiISHnkT0GpKNmHImyhowphrS3ycIbzxg867HGRrLUsWn+QSZ9u4fS5dFe9eqVQptzSluvb1/FgdyIiIiJSXvlTUMo+YlSlgOWVc40TDzt4OoXHF25kxfbjbvVhneoz8cbWVK0Q4qHORERERKS886egtN352qKA5c2drwVdwyRlxOGwvLP6V6Yv3kby+Qt3Vq9XNZxpQ9rTo0UND3YnIiIiIv7An4JS9nVB/Y0xATnvfGeMiQC6ASnAj55oTrLsPn6WCfPj+Cn+lKtmDPzhmkb8fUBLKob605QVEREREU8J8HQDZcVauxv4CmgEPJBr8SSgIvDfHM9QkjKUnungtRW7uP7F79xCUtMaFZn312uIHdRWIUlERMQP7dixg7FjxxITE0NkZCTBwcFERkbSpUsXHnnkEdauLejpKWVrzpw5GGPcPipVqkT9+vXp168fTz75JDt37vR0m1IMPv2bpzFmMFm34gbIvhX3NcaYOc7/PmGtfSTHKvcDq4CXjDF9ga1AF6A3WafcPV7qTUsemw4mMn5+HJsPnXHVggIM9/VqygO9mxEWHOjB7kRERMQTrLVMnjyZyZMn43A4iImJYfjw4URGRpKUlERcXBwvv/wyzz33HK+88goPPJD77+Ce0aFDBwYPzvr1NCUlhWPHjrF69WqmTJnC008/zd/+9jf++c9/EhTk07+G+wVf/z/UEfhDrloT5wfAr4ArKFlrdxtjrgQmAwOBG4DDwEvAJGttQql3LC6p6Zm8tHQnb3y7h0zHhZsNtq9XhRlDo2lTt3Iha4uIiEh5NnnyZGJjY4mKiuK9996jW7duecYcO3aMmTNnkpjoPffi6tixI7GxsXnqy5cvZ+TIkbz44oukpqYya9assm9OisWnT72z1sZaa00hH43yWWe/tfaP1to61toQa21Da+1ohaSy9VN8Aje8+B2vrdjtCkmhQQE8en0rFt7fVSFJRETEj+3Zs4epU6cSEhLC4sWL8w1JADVr1uSZZ55h3LhxbvUdO3YwYcIErrzySmrUqEFoaCgNGzbk3nvv5cCBA3m2Y63l7bffpmvXrtSoUYOwsDCioqIYMGAAH3zwQYm8p969e7NkyRJCQkL417/+xbp161zL4uPjMcYwcuRIduzYwfDhw6lZsyYBAQGsWLECgLVr1zJ69Gg6dOhAZGQkYWFhNG/enIcffphTp07lu8/ExETGjBlD/fr1CQsLo1WrVjz//PPs2bPHtT8pmK8fURIfczYtg2e/3MZ/f/jVrd6lcSTTh0bTuHpFD3UmIiIi3mL27NlkZGRwxx130LZt2yLH5z6NbcGCBcyaNYvevXvTtWtXQkJC2Lx5M//+97/59NNP+fnnn6lXr55r/OOPP860adNo3Lgxv/vd76hSpQqHDx/mp59+4qOPPmL48OEl8r5atWrF7373O+bOnct7771HTEyM2/Ldu3fTpUsXWrRowZ133klKSgqVK2f98fjNN99k4cKF9OzZk379+pGZmcm6det4/vnnWbx4MatXryYiIsK1rdTUVPr06cO6dev4zW9+w5133kliYiJPP/003333XYm8n/JOQUnKzPLtx3h8wUYOJaa6apVCg3j0hlbcflUDAgL0HF4REZFCxRb0OEgvFHvpp8N9//33APTp0+eS1h8xYgQPPfQQoaGhbvWvvvqK66+/nqlTp/L666+76m+88Qb16tVj06ZNVKhQwW2dEydOXFIPBenVqxdz585lzZo1eZatXLmSRx99lGeeeSbPskcffZRXX32VwED3a7ffeust7r77bl577TXGjx/vqv/jH/9g3bp13Hbbbbz77rsYk/V71uOPP54noEn+FJSk1CUkn2fKZ1tY+MtBt3rfVjWZems76lQJ91BnIiIi4o2OHDkC4HbUJ1t8fDxz5sxxq1WtWpUxY8a4Ps9vPYD+/fvTtm1blixZkmdZcHBwnhACUL169eK0XqTs3o4fP55nWa1atXjqqafyXa9hw4b51v/0pz8xduxYlixZ4haU3n77bQICApg2bZorJAFERUUxZswYJk6ceDlvwy8oKEmpsdbyWdxhYj/ZzMnk8656ZMUQnrq5DYM61HX7xhURERGBrN8hgHx/T4iPj2fSpElutYYNG7oFJWst77zzDnPmzGHDhg2cOnWKzMwLD7EPCQlxW//OO+/k5Zdfpm3btvz2t7+lZ8+eXHPNNVSpUvJH8Ap7bx06dMhzFCxbeno6b7zxBu+//z5btmwhMTERh8P1WFAOHrzwB+kzZ86we/duoqKiaNSoUZ5tde/e/TLfhX9QUJJScSQxlYmLNvHN1qNu9Vs61uXJm9pwRaX8fwiIiIhIIS7jdDZfUqdOHbZt2+b2y3+2Xr16ucJGRkYGwcHBecaMHTuWmTNnUqdOHQYMGEC9evUID886g2XOnDn8+qv7tdIvvPACTZs25T//+Q/Tp09n+vTpBAUFccMNN/Dcc8/RrFmzEntvhw4dAqBGjRp5ltWuXTtPLdvw4cNZuHAhTZo04ZZbbqF27dquUDVz5kzS0tJcY8+cyXrkSq1atfLdVkF1caegJCXKWsv7P+3nmc+3kpSW4arXrhzG07e2o29rfWOKiIhI4bp168by5ctZunQpf/rTn4q17rFjx3jppZdo164dq1atcrvBAcB7772XZ53AwEBGjx7N6NGjOXbsGCtXruT999/no48+YvPmzWzevLnAIz3FtXz5cgC6dOmSZ1lBZ9r8/PPPLFy4kH79+vHFF1+4hUOHw8Gzzz7rNj77BhBHj7r/wTpbQXVx59O3BxfvEn8imTveXM2jCza6haQ7uzTg67E9FJJERETkoowcOZKgoCDmzZvH1q1bi7Xunj17cDgc9O/fP09IOnDgAHv27Cl0/Zo1azJkyBA+/PBD+vTpw+7du9m0aVOx30N+tm3bxkcffYQxhjvuuOOi19u1axcAgwYNynMEbc2aNaSkpLjVKleuTJMmTTh48CDx8fF5trdy5criN++HFJTksmU6LG9+u4eBL37LD3tOuuqNq1fk/Xuv5ulb2xMRlvewuIiIiEh+mjZtysSJEzl//jzXX389q1atynfc6dOn89Syr8lZuXKl23VJZ8+e5Z577iEjI8NtfFpaGkuXLnWdzpctPT2dhISsx2zmvhPepfi///s/Bg4cyPnz57nvvvvo0KHDRa+b/Z6yn6mU7dixYzzwwAP5rnPXXXfhcDh49NFH3d7b/v37mTlzZrH790c69U4uy7YjZxg/L44NBy6cMx0YYLjn2iaM6decsOC8d48RERERKcqTTz6JtZYpU6bQrVs3OnXqROfOnYmMjOT06dPEx8fzzTffANCjRw/XerVr1+a2227j/fffp2PHjvTv35/ExES+/vprwsLC6NixI+vXr3eNT0lJoV+/fjRq1IguXbrQsGFDUlNT+frrr9m6dSuDBg2idevWF933+vXriY2NBbJC2NGjR1m9ejVbtmwhICCAsWPH5jlVrihXXXUV3bp1Y8GCBXTt2pXu3btz9OhRFi9eTMuWLalbt26edcaNG8eiRYt4//332b59u+vr8OGHH9KjRw8WLVpEQICOmRRGQUkumbWWhz/cwOZDZ1y11nUq8+zQaNrX96HnPIiIiIjXMcYQGxvL7bffzqxZs1i+fDnvvvsuycnJRERE0LRpU+677z5GjBiR57lAb731Fk2aNOGDDz7g1VdfpUaNGgwaNIjJkyczdOhQt7EVK1ZkxowZLF++nFWrVrFo0SLX9l9//fViXyO1YcMGNmzYAGQdiapWrRqtWrVi2LBhjBgx4pJuDBEYGMgnn3zCxIkT+eKLL3jppZeoV68ed999NxMnTqRNmzZ51gkPD2f58uU8+eSTzJs3jxdeeIHGjRvz2GOPce2117Jo0SLXtUySP5P7MKP4HmPM2piYmJi1a9eW+b5/2XeKIa+vIjgggNH9mnNvjyYEB+qvEyIiIpci+3qc4hzBECmuN998k3vvvZdZs2bxl7/8xdPtXJKL/V7p1KkT69atW2et7VTcfeiIklyW3zSoxuRb2nFNkytoVrOSp9sREREREadDhw7lOS1v//79TJkyhaCgIG666SYPdeYbFJTkso24Ov8nRYuIiIiI5wwdOpT09HQ6depE1apViY+P57PPPuPcuXNMmzaNevXqebpFr6agJCIiIiJSDo0YMYL//e9/zJ8/n8TERCpVqkSXLl148MEHGTJkiKfb83oKSiIiIiIi5dD999/P/fff7+k2fJauuhcREREREclFQUlERERERHxGWd21W0FJRERExEsYYwBwOBwe7kTEe2UHpezvl9KioCQiIiLiJUJDQwFITk72cCci3iv7+yP7+6W0KCiJiIiIeImIiAgAjhw5QlJSEg6Ho8xOMxLxZtZaHA4HSUlJHDlyBLjw/VJadNc7ERERES8RGRlJcnIy586d48CBA55uR8RrVahQgcjIyFLdh4KSiIiIiJcICAggKiqKhIQEkpKSSEtL0xElESdjDKGhoURERBAZGUlAQOmeHKegJCIiIuJFAgICqF69OtWrV/d0KyJ+TdcoiYiIiIiI5KKgJCIiIiIikouCkoiIiIiISC4KSiIiIiIiIrkoKImIiIiIiOSioCQiIiIiIpKLgpKIiIiIiEguRg8x833GmJPh4eGRrVu39nQrIiIiIiJeY+vWraSkpCRYa68o7roKSuWAMWYvUBmIz2dxAFALOAo4CtjE5Y5p5XzddtFNe5+L+Rp4+z4vd3vFXb844y92bFHjNA+9f59lPQ+Lu05J/Ewsarnmonfs73K26QvzsKgxmofesU9v/plYFv82g+fnYiPgjLW2cXFXVFAq54wxdYGDQD1r7aHSGGOMWQtgre1Ukr2XpYv5Gnj7Pi93e8VdvzjjL3ZsUeM0D71/n2U9D4u7Tkn8TLyI5ZqLXrC/y9mmL8zDosZoHnrHPr35Z2JZ/NvsXO6zc1HXKImIiIiIiOSioCQiIiIiIpKLglL5lwRMcr6W9hhf5on3V9L7vNztFXf94oy/2LFFjdM89P59lvU8LO46JfHzrrzPQyj791ga+7ucbfrCPCzuPn2RfiaW7lzUv81F0DVKctl8+dxTKT80D8VbaC6KN9A8FG/hy3NRR5RERERERERy0RElERERERGRXHRESUREREREJBcFJRERERERkVwUlERERERERHJRUBIREREREclFQUlERERERCQXBSUREREREZFcFJRERERERERyUVASr2KMecwYY40xr3i6F/E/xpgHjDFxxpgzzo8fjDE3erov8S/GmEeNMT855+BxY8ynxph2nu5L/I8xpocx5hNjzEHnv80jPd2TlG/GmPuNMXuNManGmLXGmGs92Y+CkngNY8zVwD1AnKd7Eb91ABgPxABXAsuARcaYaI92Jf6mF/Aa0BXoA2QA3xhjIj3ZlPilSsAmYDSQ4uFepJwzxgwHXgSeAX4DrAIWG2MaeKwna62n9i3iYoypAqwjKyg9CWyy1j7o2a5EwBiTADxqrX3D072IfzLGVAISgcHW2k893Y/4J2PMWeBBa+0cT/ci5ZMxZjUQZ629J0dtJzDPWvuoJ3rSESU/ZIwZZox52RjznfPUDmuMmVvEOvWNMf8xxhwyxqQZY+KNMTONMdVKqK1/kfWNsKyEtic+wEvnYvZ+Ao0xt5H1F9VVJblt8S7ePA+dIsj69/pUKWxbvIgPzEXxY6U5P40xIUAn4Ktcm/iKrKPrHhHkqR2LR00EOgBnyTrVqFVhg40xTcn6RbEm8DGwDehM1qH4gcaYbtbak5fajDHmHqAZMOJStyE+y6vmonMf7YEfgDBnX7daazdezjbF63ndPMzlRWA9WfNSyjdvn4vi30pzflYHAoGjuTZzFOhXUm+guHREyT89BLQAKgP3XcT418ia5KOstYOttROstX2AF4CWwNM5Bxtjpjr/ylDYRy/n2JZknYt6p7X2fMm9RfERXjMXc9gOdASuBl4H3taF9OWeN87D7HWfB7oDQ621mZf+FsVHeO1cFKGU56dT7muCTD61MqNrlPyc8wficuAda+3v81neBNgNxANNrbWOHMsigMNkTeKa1tpkZ706WX8ZKMw+a+05k3UHndlAzl8AAsn6pnAAFa21aZf05sSneHouFtLXN8Cv1to/F+sNiU/ypnlojHkBuA3oba3ddqnvSXyTN81F57q6RklcSnp+Ok+9Owfcbq39KMfYV4F21tqepfh2CqRT76QofZyvX+Wc5ADW2iRjzPdAf7L++r7UWT8BnLjI7S8Cfs5Vmw3sJOtIk44ySbbSnosFCQBCL3MbUn6UyTw0xrxIVkjqpZAkBfDUz0SRi1Gs+WmtPW+MWQtcB3yUY/h1wPyyaDg/OvVOitLS+bqjgOU7na8tLmXj1trT1tpNOT+AZCDB+bkOeUq2Up2LAMaY6caYa40xjYwx7Y0x08i6VfM7l7pNKXfKYh6+CvwRuB04ZYyp7fyodKnblHKpLOZiJWNMR2NMR7J+Z2zg/Nxjt2sWn3Ep8/N5YKQx5m5jTGvnH4zqArNKqcci6YiSFKWK8zWxgOXZ9apl0Iv4t7KYi7WBuc7XRLKe6XW9tXbJZWxTypeymIf3O1+X5qpPAmIvY7tSvpTFXLySrNOrsk1yfrwNjLyM7Ur5V+z5aa39wBhzBVk3jahD1jO8brDW/lpqXRZBQUkul3G+ltiRH2ttr5LalviVy56L1tqRJdOK+LGSmIem6FEiRSqJubgix3ZESlK+89Na+xpZN4HwCjr1ToqSnfirFLC8cq5xIqVFc1G8geaheAvNRfFm5WJ+KihJUbY7Xws6x7m587Wgc1BFSormongDzUPxFpqL4s3KxfxUUJKiZJ+b3N8Y4zZfnLd37AakAD+WdWPidzQXxRtoHoq30FwUb1Yu5qeCkhTKWrsb+ApoBDyQa/EkoCLw3+xnNIiUFs1F8Qaah+ItNBfFm5WX+akHzvohY8xgYLDz09rAAGAP8J2zdsJa+0iO8U2BVWQ9XfljYCvQBehN1iHTrtbak2XTvZQnmoviDTQPxVtoLoo388f5qaDkh4wxscBThQz51VrbKNc6UcBkYCBwBVlPVF4ETLLWJpROp1LeaS6KN9A8FG+huSjezB/np4KSiIiIiIhILrpGSUREREREJBcFJRERERERkVwUlERERERERHJRUBIREREREclFQUlERERERCQXBSUREREREZFcFJRERERERERyUVASERERERHJRUFJREREREQkFwUlERERERGRXBSUREREREREclFQEhERERERyUVBSUREvIIxxhbxMdLTPYqIiP8I8nQDIiIiuUwqoL6+TLsQERG/Zqy1nu5BREQEY4wFsNYaT/ciIiKiU+9ERMTnGGPmOk/Ha2CMGW2M2WiMSTHGfJNr3PXGmMXGmJPGmDRjzG5jzLPGmMoFbHeAMeZ7Y0yyMSbBGLPQGNMix/7q5xjbz1mbWMC2DhhjdhWw7E5jzApjzGljTKoxZosx5jFjTEiucUHOfXxjjKlhjPm3MeaI871sMsbcVcjXaKAx5jNjzDHn+P3GmEXGmD7O5Tc5t/2vAtYPc34NjubuS0TEH+jUOxER8WWvAt2BL4DPgfPZC4wxk4EngJPAp8BxoAPwd+B6Y0xXa21SjvHDgXeBNOAD4AjQA/gB2FJSDRtj3gbuAvYB84BEoCvwNNDHGDPAWpuZa7VIZx/ngA+BcOB3wNvGmExr7Tu59vE08BiQBCwCDgB1yfpa3QEsI+trFg/cbox5OOfXwul3QDVghrX2PCIifkZBSUREvIoxJjafcry1dk4+9Y5AR2vtr7m2cR1ZIWklcJO1NjHHsruBN4EnyQpNOI8wzQIygW7W2l9yjH8ZePAy3lLOvu4mKyR9BNxlrU3NsWwKMBH4K1kBMKffAP8C7s8OUcaYl8i6bms88E6O7dxAVkjaDfSw1h7KscyQFZiw1jqMMW8A04A7ne8/p3sBS9bXSkTE7+gaJRGR/2/v3kKsrKIAjv8XKYgXxiZFLYSIKJUKLKisxIjKHrrQUxeSCAMjKqinhorKCqKHKIIiSQky6CEQkm76ECaCGEQYWeJDSBNpUpoY2gVXD/sbO23PcS49zHjO/wfDx9ln72/vMy8z66xvr60JYWiPUgdbMvPalr7rKf/cP5SZdVBBRGwEbgYWZObuNu9/DZyVmWc3r+8F3gbWZebKqu+ZwF5gBjA/Mweb9uuBzcBTmfl8mzkGgWOZeX417wXA7Mw8XPWfRMl6fZuZV7W0/QUcAeZl5pFqzDZgCTAtM482bR8DNwG3ZubGel3V+NmUbNOuzFzc0r4I+AbYnJk3nuoektStzChJkiaUURZz2NGhfQnlEbq7ShLlJJOAeRHR12SbLqlLVDYAAANISURBVG3at7RZz8GI2AlcPYp1nSQiZgAXAfuBxzqs6xiwsE377jpIavxAeWxvJnC0absSOA58OtyaMvNARLwP3B0Rl2fm0O9zVXOts0yS1DMMlCRJp7N9Hdr7gQCeHmb8dMoeob7m9f5RzjMa/c11zjDr+rtN26Fh+p7R0tYHHBjFvqLXKfuWVgE7ImIKsILymT8Y4T0kqetY9U6SdDrr9LjeYUqwEMP8/Nj0H9rDNKfD/ea2aTveXDt96dhXvR6a44th1jS5w/1G6jdg1kgr1WXmNmAncGdE9PFvEYe1mdkuaJOknmCgJEnqRtuB2RFx4Qj7f9lcl9VvNHuULmkz5mBznd9mzAJKtuqEzDwE7AYujoiZI1zXWGyn/H1fPooxbwBTgXsoRRyOYxEHST3OQEmS1I1ebq5vRcS8+s2ImB4RV7Q0baBkYlZExOKq+2pKIYfaLkqRhdsjYlbLvacCr55iXVOAtU32pl5Xf5v5R+u1obk6fPZz2oxZTyklPkDZi/VJXUlQknqNe5QkSV0nMzc1B8E+B+xpKsF9T8nynEvJHH1GqYxHZh6OiAco5yhti4jWc5QWUsqMX1PN8UdTOnwA+CoiNgCTKZmcvbTZ75SZayLiMkrWZllEbKKcp9QPnAcspWRyxlyOPDM/iogXgceB7yJi6Byluc1n2ArcX405EhHvAA82TW+OdX5J6hYGSpKkrpSZL0TEVuARSpbkNkrWaJBSze3dqv97EXGQcr7SHZQKdJ8D9zVt/wmUGk8CvwMrKcUQ9jX3XQ3s6bCuVRHxYdP/BkrFul8owdVLlOzO/5KZA03p8IeBW4BpwM+UKoGd7r+OEigNUg7vlaSe5jlKkiQNo+XcphPnKHWbloN4n8nMZ8d7PZI03tyjJElSj2sOtn0U+BNYM87LkaQJwUfvJEnqURGxlLIP6zpgEfBKZv40vquSpInBQEmSpN61HHgC+JWyb2tgfJcjSROHe5QkSZIkqeIeJUmSJEmqGChJkiRJUsVASZIkSZIqBkqSJEmSVDFQkiRJkqSKgZIkSZIkVQyUJEmSJKlioCRJkiRJFQMlSZIkSaoYKEmSJElSxUBJkiRJkioGSpIkSZJUMVCSJEmSpMo/oAXhgH3AqpYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 277,
       "width": 421
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(xax, gw_vec, label = \"GW Emission\")\n",
    "plt.loglog(xax, drag_vec, label=\"Gas Drag\")\n",
    "plt.xlabel(\"Frequency\")\n",
    "plt.ylabel(r\"$df/dt$\")\n",
    "plt.legend(loc=4)\n",
    "plt.tight_layout()\n",
    "# plt.savefig(\"test_run_m1isOne.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fname = f\"fdot_and_f_m1_{m1:.1f}_m2_{m2:.1f}_T{temperature:.0f}.txt\"\n",
    "to_print = np.array([xax, gw_vec, drag_vec]).T\n",
    "np.savetxt(fname, to_print)\n",
    "# add header\n",
    "header = \"\"\"# These are calculated in BinaryResidenceTime_AGN_Disk-ChangeGas.ipynb\n",
    "# using \"\"\"+f\"m1={m1:.1f}\"+\" and \"+f\"m2={m2:.1f}\"+\"\"\"Msun\n",
    "# GW_frequency roughly in the LISA range, from 10^-4 to 1 Hz\n",
    "# the gw_fdot assumes pure GW-driven inspiral (Peters 1964)\n",
    "# the drag_fdot assumes no-feedback and the drag is calculated using a\n",
    "# temperature of \"\"\"+f\"{temperature:.0f} K\"+\" and density \"+f\"{rho_0} g cm^-3\"+\"\"\"\n",
    "# for the RSG envelope as representative\n",
    "#\n",
    "# GW_frequency\tgw_fdot\tdrag_fdot\n",
    "#\n",
    "\"\"\"\n",
    "# print(header)\n",
    "with open(fname, \"r+\") as f:\n",
    "    s = f.read()\n",
    "    f.seek(0)\n",
    "    f.write(header+s)\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "# don't run this now\n",
    "def t_baruteau(Temp, SigmaSurf, bin_Mass, SMBH_mass, a_SMBH, a_bin):\n",
    "    \"\"\"\n",
    "    temperature in Kelvin\n",
    "    SigmaSurf input in g/cm^2\n",
    "    SMBH_mass in solar masses\n",
    "    bin_mass in solar masses\n",
    "    a_SMBH in pc\n",
    "    a_bin in AU\n",
    "    Answer in YEARS\n",
    "    \"\"\"\n",
    "    const = 1/5.443310539518173e-10\n",
    "    ans = const*(Temp/2000)**3*(SigmaSurf/100)**(-1)*(SMBH_mass/3e6)**(-0.5)*(bin_Mass/2)**(-2)*(a_SMBH/0.1)**(3/2)*(a_bin/50)\n",
    "    return 900*ans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "M1 = 99     # solar masses  \n",
    "M2 = 1     # solar masses\n",
    "Mtot = M1+M2\n",
    "qratio = M2/M1\n",
    "Temp = 10**5  # K this was 1e5.5... \n",
    "SigmaSurf = 10**(5.5) #g/cm^2, Figure 2 of Sirko & Goodman 2003\n",
    "a_bhs = 1e-2 #pcs\n",
    "alpha = 0.03\n",
    "SMBH_mass = 1e7\n",
    "a_SMBH = 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# All in AU\n",
    "t_bteau100 = [] # distance in AU, surface density /1000\n",
    "t_bteau60 = []\n",
    "t_gws100 = [] # distance in pc\n",
    "t_gws60 = []\n",
    "t_acc60 = [] \n",
    "t_acc100 = [] \n",
    "xaxis_t = []\n",
    "t_all100 = []\n",
    "t_all200 = []\n",
    "# t_baruteau(Temp, SigmaSurf, bin_Mass, SMBH_mass, a_SMBH, a_bin) # functional form\n",
    "\n",
    "for ii in range(1000*1000):\n",
    "    ii = (ii+1)/1e5\n",
    "    xaxis_t.append(ii)\n",
    "    \n",
    "    # Equal Mass 60 Msun binary\n",
    "    t_bteau60.append(t_baruteau(Temp,SigmaSurf,60, SMBH_mass,1e-3, ii))\n",
    "    t_gws60.append(t_gw(AU_to_pc(ii), 1, 60))   \n",
    "    t_acc60.append(t_ACC(AU_to_pc(ii),30,30,SigmaSurf,Temp, alpha)) \n",
    "    \n",
    "    # Unequal Mass 100:1 Msun binary\n",
    "    t_bteau100.append(t_baruteau(Temp,SigmaSurf,100, SMBH_mass,1e-3, ii))\n",
    "    t_gws100.append(t_gw(AU_to_pc(ii), 1/99, 100))   \n",
    "    t_acc100.append(t_ACC(AU_to_pc(ii),99,1,SigmaSurf,Temp, alpha))  \n",
    "\n",
    "\n",
    "    t_all200.append(min(t_baruteau(1e5,10**5.5,200,1e7,1e-3,ii), t_gw(AU_to_pc(ii), (1), 200), t_ACC(AU_to_pc(ii),100,100,SigmaSurf,Temp, alpha) ))\n",
    "    t_all100.append(min(t_baruteau(1e5,10**5.5,100,1e7,1e-3,ii), t_gw(AU_to_pc(ii), (1/99), 100), t_ACC(AU_to_pc(ii),99,1,SigmaSurf,Temp, alpha) ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# What are LISA Band Limits? 1e-4 to 1e-2 Hz\n",
    "lisa_min = a_to_f(Mtot, 1e-4)/499.00478 #seconds converted to AU\n",
    "lisa_max = a_to_f(Mtot, 1e-2)/499.00478 #seconds\n",
    "\n",
    "ligo_min = a_to_f(Mtot, 10)/499.00478 #seconds converted to AU\n",
    "ligo_max = a_to_f(Mtot, 80)/499.00478 #seconds to AU, and also up to isco freq = 76 Hz\n",
    "\n",
    "et_min = a_to_f(Mtot, 1)/499.00478 #seconds converted to AU\n",
    "et_max = a_to_f(Mtot, 1e3)/499.00478 #seconds to AU, and also up to isco freq = 76 Hz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.ticker as tck"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda2/envs/py3/lib/python3.6/site-packages/matplotlib/font_manager.py:1328: UserWarning: findfont: Font family ['serif'] not found. Falling back to DejaVu Sans\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n",
      "/anaconda2/envs/py3/lib/python3.6/site-packages/matplotlib/ticker.py:2516: UserWarning: AutoMinorLocator does not work with logarithmic scale\n",
      "  warnings.warn('AutoMinorLocator does not work with logarithmic '\n",
      "/anaconda2/envs/py3/lib/python3.6/site-packages/matplotlib/ticker.py:2516: UserWarning: AutoMinorLocator does not work with logarithmic scale\n",
      "  warnings.warn('AutoMinorLocator does not work with logarithmic '\n",
      "/anaconda2/envs/py3/lib/python3.6/site-packages/matplotlib/font_manager.py:1328: UserWarning: findfont: Font family ['serif'] not found. Falling back to DejaVu Sans\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n",
      "/anaconda2/envs/py3/lib/python3.6/site-packages/matplotlib/ticker.py:2516: UserWarning: AutoMinorLocator does not work with logarithmic scale\n",
      "  warnings.warn('AutoMinorLocator does not work with logarithmic '\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAEwCAYAAABrH3KHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXl8m1ed7/8+Wi15kxd5zWZnc+LYTpyktBTolpRS9pJSOsDvzpSSlJnhsgw0hQ60c7nQSSl0ZhiYSaCX4XIHaFP2KTAkKVC6UJI4sRPHzmI7i3d5kWVrX87vj0eSJctbEi9xfN6vl162znnOec4jPdJH55zvIqSUKBQKhUKhmB908z0AhUKhUCgWM0qIFQqFQqGYR5QQKxQKhUIxjyghVigUCoViHlFCrFAoFArFPKKEWKFQKBSKeUQJsUKhUCgU84gSYoVCoVAo5hElxAqFQqFQzCNKiBUKhUKhmEeUECsUCoVCMY8oIVYoFAqFYh5RQqxQKBQKxTyihFihUCgUinlECbFCoVAoFPOIEmLFokYIsUMIcUAIIYUQg0KI/dHnLUKIPUII25hjW6LHSiHEnoS6ndH2MnrMzsvpe7EihHhYCHF0vsdxOQghtgkh9gohBud7LIrrAyGlnO8xKBTzihCiHGgB9kkpd0XLbMBRACnlyoRjdwD7geellPeO6WcvsBPYJaXcd7l9L0air+f22GuzUIj+CHtYSinmeyyKhY+aESsUMDC2QErpBPYC5UKIbQlVrWP+JtIyTt3l9H1dI4QoT1xFAJBSPr/QRDjK4fkegOL6QQmxQqGYdRJWARb9crxCMRYlxArFOESXlHcBB6WUB+e6byFEbXQf8kD0/wPRfeaHo/V7o3vNLdH+Yu1s0f3n/bG96DH97Y3uXQ8KIXYmtHs4tl+bsM+992r6HENs5r8t2ld5bIacuNc609cdrdsRrXs4Wlc75rV+OKHt/tjefcJYHo6OY8cU7+sVj0GxyJFSqod6LOoH2ixNAoPAAbQlZgnsH+fY2mjdnnHqHo7WbbuSvsfp70C03Y7o853RtnsBW7SsBdib0GZ/Ql1srDuB8mhfLdHnD8fGGX1+IKEPebV9TnA9Y8dqi16LnK3rjj4fTDhuJ1A7pt22Mde+J9YuoY8dieMc+/xqxqAe6mFAoVDEeE6OGlRtA/ZHZ2ubpZRj94R3jrO/O9my6+X0HaMVKJdSPh99Hps9H5DaPjNAHbAlod9aYI8QcRuiOmCllHKfEOIIkCujhmQJ7AKOJDyfiT6nRErpnMBiesauO/q/TQixU0q5Lzrm+IwXTRATVyW2M7rHvy/h3K3RNuXjvV9XOgaFAlBCrFCMh5TyoBDiXrTZ2R7g3jGH7JNS7k4siC6f7mEKptH3RKQYfkXLYsuctUCrnNz4aTzRH0Cb3cY4AuReZZ8zydVe925grxBiF3BvgpBuY8zYE0VZSrk7unz+MJA3xRivdAwKhdojVigmITZLnA3L5tnoO49kQZ0ue9H2bmNttwBPXGWfc8mkY5RSPok20y0HWhJWMiZtFxXg3dH2z87SGBQKJcQKxSTEZoWzMXuZjb5bGMclaqzL0FiiS8BPArujx35USll3NX3OMROOMWpAVSulPCilzEFbLo6tZByeoN2O6LL1HqaxwnGVY1AolBArFIyKYu6Y8r3Rv7vHOXa8Pb7Y8mXizOhy+p4uY/sCILpP60Tbf94ptAhQ+9GWwGPtUsYdnfnZoscdBnJjs+Mr7XMCEvd1Y6/R5eyVXsl15wKfSzj82eixsR8gTuBA1Jp5R7RdHaPv4Y7ofu590efl0ee50euwXc0YFApAWU2rx+J+oC0N70ezcJVoX5wx6+YDJFvUbouWxaygH06o28moRXQLmlXttPueYFyx/nagCdae6POjaEKxLToOyaiFcW20PjaORMvo2Dh2jjnXjuixLQnHyIS2l93nBNe0I+F1sEXHH+t3T0LZTF53OaMW1w+TYHmdUH8gof/E9ztWvj86jpboMe8aO+6rGYN6qIcKcalQLHKi/sK75ahFcsyieJdcmFGvFIoFhVqaVigWMdEgFeWJIgwgtT3iBZWMQaFYqCghVigWN7loFtM7x0SU2gM8N79DUygWB2ppWqFY5ERFdyfaPmgd8KzU3G0UCsUcoIRYoVAoFIp5RC1NKxQKhUIxjyghVigUCoViHlFCrFAoFArFPKKEWKFQKBSKeUQJsUKhUCgU84hKg3iVCCHswJ3Rp+cB3/yNRqFQKBSzSBqwIvr/b6WUjpnoVAnx1XMn8P/mexAKhUKhmFM+BPznTHSklqYVCoVCoZhH1Iz46jkf++f73/8+69atu6JOHA4Hg4OD5OTkjFsfCAQwmUxX1PeVIDu8yP4AIm9mzhkKhTAYJr/dYucTpZYZOeflEnv97Xb7uPU+n4+0tLQ5HtXMMZ3xD3SN4OrzkZV/bV5nIBDEZDLO6jli159bnDFjfY74QnQ6vSyzGRf0PQSL43MwEU1NTXz4wx+OPT0/U2NSQnz1xPeE161bx+bNm6+ok87OThwOx4Qi4Pf7MZvNVzbCK0BmuJE9PkThzHzggsEgRuPkX6Cx84ny9Bk55+USe/1LSkrGrfd6vVgs8/MjYSaYzvh7L7gY7HKTUzw/78FUzMXnIHb9BcuzZqS/DqeXv/run+ka8vODBzZQtSx/RvqdLxbD52CazJg9kFqaVigUillCSslD3z/KmZ4Rhn0hHvpBPb5geL6HpbjGUEKsUCgUs4QQgifuqSLdpMeoF3zy9nLSjPr5HpbiGmPRLE0LIR4GWoGtUsrd8z0ehUKxONhQms23PrQZk17HxhLrfA9HcQ2yKGbE0eTnSCmfjz4vn98RKRSK65FIRDLgDqSU37LGzk0r8+ZhRIqFwGKZER8EjgohWgGklK3zPB6FQnGd4Q+F+cz+Bho7h/jJx96IzTp3Xg6Khc2CnRELIXYIIWxjysqFEHujdXtj9VJKJ7AH+DZweB6Gq1AormOklDz4vSP8sr6TVoebnf/3KP6QMspSTI8FJ8RCCFt0v3c/kDum+gCwJ7oEvRc4Gm1TDtiAMmCPWppWKBQziRCC+7YujT9fU5SBXoh5HJFiIbHglqajs9snhRB7EssT9oFjy891QohcIcQ2oBxwSimd0Xa1aIZbCoVCMSO8o7qETqcXKWHnW8oRSogV02TBCfEkbCVVXAfQRHcf2kx4J7BSSrlvqs6EEI8Dj13OAHw+H16v93KaxPH7/QSDQfx+/4T1c4kIBhGhMDIYnJH+gqGp+4mdT87xtcaIvf4TvYc+38LO5zGd8fv9fgKT3IfzTSCQagg14+eY4j6IIaVMEdsPb9WCwUz0Wi/0ewgW/jVczfhn69qvJyEuB5xjypxAXnQWvetyOpNSPg48PtVxQojNwBGAtLS0K47YYjabMRqNk0YNmtPIWsYQ0hBGTBEN63KYMrJW9HxiDq8zkdjrP9l7uJAjCsHU4zebg5iMoTm91y6X2R5b7Pone62+9+p5/njWwb9/aDMG/eXt8C30ewgW/jVc6fhnK7Tn9STErWiz30RsQMs8jEUxhv1DP+NM8BwVaWtYl7aG1eZVWHUL+8OsWHxIKfnyC0185+U2AL74i0a+/J4NahlacVVcT0J8GNgxpiyX6GxVMb8c9R2n0d/Enzya0boOHctNS6kwr6EiTXsskXmomEOKaxkhBGbj6Ay4sdOFNxjGarqevkoVc811c/dIKZ8XQuwRQtiiRlm1wBEpZd18j22xE5ZhzgWSt+8jRGgLXKAtcIFfDx8AwIqFis41VEQqqciuoCK7gmxT9nwMWaGYkM/cuZb2QS/eQJh//sAmLCb181FxdSw4IY76Bu+MPt0hhNgX3QMG2I5mlHUU2Mxl7gsrZgcdOr5e9BXOhVpp9p2hyX+GC4GLSGTScR681HnrqWurj5eVWEpYZ1unCXNWBeWZ5Rh0C+62VVxHCCH46o4a9DqBXqeWpBVXz4L7Rou5L0UfY+taUeJ7zSGEYJlxCSutZbw16w4APBEvZ/3naPKdodl3hmb/GZzhoZS2nd5OOr2dHOo6BIBJZ+LrW7/O6qzVc3oNisXJ0QuDvNbSx9/enny/mQwLLgSD4hpmwQmx4vrAqrNQY6mixlIFaEYw3d2XaM64wOm0NpqHmjnnOkdIhpLahSIhlqYvTSobDg7z9KmnWZe9jvW29VTaKufsOhTXL78+0cUnnz2OPxTBZjXxoRuXz/eQFNcpSogV1wRCCIp0BRRnLuP28rsACIQDtAy30OxqptnZTNNQE+mGdNL0yS4Ep4dO80rvK7zS+wprstbwjTd8I6k+EA5g1BmVZati2kgp+emxDvyhCAD/fOgs791USrpZfWUqZh51VymuWUx6E+ts61hnWwfLtDJfONWhvnmoOf5/RXZFSv0P237IL9t/SUWWZgC2zraOtVlryTBmzNrYFQsbIQT/9IGNfGDfnxj2hfiPv9qqRPh6IBJB9J8DWxFk2Od7NHHUnaVYUIydDQPcVnwbueZcmoea2ZK3JaW+aaiJ4eAwh/sPc7h/NOfH0vSlmjBna8Zg6TJ9VseuWFhYTQae+R9b0esEuekqk9KCxDsIHUeh/Qi0H4b2I6T5nPC2J+EN1445kRJixYKn1FpKqbWUu5fcnVInpaTb2z1uu0vuS1xyX+JAp+Y+ZdaZWZu9li2lW6i2V1NjryHfkj+rY1dcGwx4AjQ5Rrh5eXIeGXvmtRthTDGGcAh6T8UFl/bD0H92/GPbDyshVijmCiEE3735u3R6O2ke0vaZm53NtI60EpbJaer8ET8Ngw00DDbEy0rSS6i2V1OVX8XarLVsKtmESa9mR9cTFwY9fOKFRhxuP9945wbKxMyFdVXMIsPdUdGNCm/nMQh6pmwmLbkI07W1+qWEWHHdI4SIz5rvKNbcp/xhP+eGz9HkbIoLdJ+/L6Vtp7uTTncnvzn/GwD2bd/HTSU3JR0zXvB/xcJASsmX/3CWzmHN9mD3fzfx3TsqyZnncSnGEPRBV/2o8HYchaFLU7fTGaCoCkq3wJKtsGQLPksxFqt19sd8GSghVixKzHozlbbKJFen0x2ncRgctIfbqXfUc6r/VJJxmECwIX9DUj/uoJt3/PQdbMjbwMaCjTyw4QElygsIIQT/cMdaHvjJcdyBMF+8bQ3pRhUpa16REgbbosvL0SXm7hMQmUYmuKxSWBIT3a1QXAPGMTHtrzBD3myihFihiJJrymWtfS0lJVoqu2AkyNnBs5xwnKChr4F+Tz+ZpsykNif7TtLn7eP37b+n093JR6o+klTf7e5mODBMeXY5ep36gr8WKc5M4+m7NxCRkvUFmQx2ued7SIsLnws665L3dj39U7czWKBkU4LwboGsktkf7yyghFihmACjzsj6vPWsz1vPfdw3bn7ak30n4/9X26tT6n967qd86/i3SDemsyF/A9X5mhFYlb2K3LTclOMVs0soIukZ8VGalTxLqrArV7Y5IRIGx+nkvV1HM4wJdzsueatGBbd0CxRWgv762M9XQqxQXAV/teGvuG3ZbTQ4GlielRp5qcGhGX65g25e73qd17tej9ctzVwaNwSrsdewNmctxuvki+VaxBMM8+hvm2hyjPB/7tlISdbs5JZVJDDigI4jCXu7xyAwPHW7tOzovm50tlu6GazX7w9XJcQKxVWgEzrKs8spzy4ftz7fkk9eWh79vtSltkvDl7g0fIkXWl8AtDja6/PWU22vjrtPFVoL1Z7zDPHwb07x53YtP8ynfnWSZ+7ZSIZKXzhzhALaXm774VHxHTw/dTuhg4LK5L3dvFWgWzzxvNVdqFDMIl+6+UtIKelyd9HgaKDeUU9DXwNN/U0ExxifBCIBjjuOc9xxPF5mt9jjwrx9+XaWZi4dewrFNHlg8zKOdQ4RjEhuKcvHqoyyrhwpNavl9sPQflT721UPYf/UbTMKR5eYl2yF4o1gXtxbA0qIFYpZRghBSUYJJRkl3FU2Gkf79MBpGvqi4uxooGOkI6Wtw+vg0MVDHLp4iFW2VSlC3DnSSVF6ETqxeGYPV0ptSTaP3b4GTzDMe9YXz/dwFhYBt+an234Y04XXoasORnqmbqc3aUIbF94tkL0U1CpPEkqIFYp5wKQ3UWWvospexQfXfRCAfm8/J/pO0OBooMHRwIm+E3hCowEKqvOTjcG8IS9v/8nbsRqtVNur+cbt31C5mhMIhCOY9Mk/UO5cXTBPo1lARCLQfy7ZoKq3EaSWAGPSdYScFaPLy6VboGgDGFR0sqlQn1qF4hohz5LHrUtv5daltwIQjoRpGWrhhOMEF4cvYkuzJR3f1N9ESIZwBVx0jnSmiHCLs4Uj3UeotlezOmfx5G+WUvJ/6i7xh7Z+/v3d1WoJeio8A9F4zFHR7TgCvtTc4CmYMqG0NsGgass1lUhhIaGEWKG4RtHr9KzJWcOanDXj1vf7+skx5zDoHxzXdep3l37HP9f9MwAWg4WKnAo2Fm6kJl9znyqwXp+zwz0vneMnp7T44o8eaOard63HoFNLoQCEg9DTGDWmisVjPjeNhgLsFbBkC4HCjZjK3gj2taB842eERSnEQgiblNI53+NQKK6G7cu3s23ZNtqH24kQSamPuU6Btox9zHGMY45j8bKi9CKq80cttCtyU1NILkTWJvgEB8MRAuEIhsUqGK7OpMxDdB6D0DQiS1nzkg2qSmohLQuAsNcLFssUHSguh0UjxEKIh4FYuo3N8zkWhWKmEEKwNGt8S+pbl96KSW+iwdFAl7srpb7b3U23u5vfXvgtAAZhYGXGKsr1a9moq2F9ViXFaaULzn3qveuL6XL5cXj8fP6W1Rj1i8SQLehNjsfcfgRcqQaAKegMUFSdbFCVU6YMquaQRSPEQJ6UcuV8D0KhmCvuWX0P96y+B4BeTy9HOo7Q7GqmwdFAY19jUhxtgJAMcXq4mdM082vnzwHINtqoyFzPvUs/wKachfP79WNv0IKrLLQfEdNGShhoTZjtHoaekxAJTd02e+lodKolW6G4OjUes2JOWRRCLIQoB2qFEBK4V0r5/HyPSaGYSwqsBdy25Dbutmg5m0OREGcHz2oW2n2alfZ51/mUdkNBJ68PvMrbS96VUvey4w+UWJawPH0FejE/S7+vXxrkxdY+dr9lFboE0b3uBNg3lJDgPiq+3oGp2xmtyfGYS7dAlnLdutZYsEIshNgBHEzc640K7m7gALAd2C2ldEopW4HtQoha4JAQ4qDaI1YsZgw6A+vy1rEubx33cR8AQ/4h/tj4OnWdxzgfOUvT8ClGQlo4wnWZ65PaByIB/vepxwjKIFa9lWdv+hlWw9zmeP1FUzdPvHSOcESSZTbwNzeWzen5Z41IGHqbEkJDHtHiM08rHvPq5L3dgvWgX7Bf84uGBfcOCSFswE5gD7ASSBTUA8B2KWWrEKIVOBo9BgApZZ0Q4jmgHKibu1ErFNc+2eZsbsi7kdWBKnKK04nICB3eS5wbOUuuOS/p2JaRswSlFhksx5SbIsJnhpt5/tKPWJdVybqsSlZmrMaom9k42s2OEcIRTZx+faaXD29aQpZ5AcbqHulNXmLuPAaBkanbpdmSMw+VbgaLyqS8EFlwQhydyT4phNiTWB6dIROd/cZEN1cIsU1KeXBMH0qEFYop0AkdS63LWWpNTWahQ8eb8t/CKVcj67IqU+qPO49xqPcAh3oPAGAUJlZnrmFdViXro+JcYL66ONqfftNKukb89I74efruyoUhwiE/tDcmGFQdBufFqdsJvZZtKHG2m7tyUcVjvp5ZcEI8CVuB1jFlA2h7w7XR+meBvXM9MIXiemNt1jr+YcMTSCkJRAIp9U2uxqTnQRnglOskp1wn+XG0LNeUF58xr8+qZE1mBRb99I2GDDrBl7dXIKUk/VpM3iClJrIJgTLSuuohnPp6pZBRlJwEoWQjmOZ26V8xd1yDd+8VU07yMjXR53lSyt2X25kQ4nHgsctp4/P5xs1ZOx38fj/BYBC/f/yg6ROVzxYiGESEwshgcOqDp0EwNHU/sfPJOb7WGLHXf6L30OfzjVu+UJjO+P1+P4FJ7sMJ24WSj39/8V9Qmb6B5pEmmodP0eXvTGkzEOjnlb6XeKXvJUCbZa+wlrE2cz0VGeuoyFxPadoSdEJH94ifo50utq/ITuojZiLm94cva7yTEZjiPpi4oRtd1zF0nXXouurQdR5FuB1Jh4w3/5d6M5GiaiIltURKNiNLNiMzS5Ldh8LAFX63zDSL4XMwG20n43oS4lagdkyZDWi5ks6klI8Dj091nBBiM3AEIC0tDcsVOrqbzWaMRiNm88RxWSerm2mkMYQ0hBHGmVvuM07RV+x8Yg6vM5HY6z/Ze3il7++1wlTjN5uDmIyhq77X1pnXsy531MDLGRikyXWKJlcjTcONNLtO4Ql7ktpEiNDqaaHV08Kve34JQLo+g3cWPMBP/rScfk8Qs2Elb19XclVjm4rY9U/6WkUi0Hcm2aCq91Q8HvOk5JSNznSXbEEUbkBvME0ew/ka43r/HExEWtrs5LC+noT4MLBjTFkuUZFUKBTzh82Uw035N3NT/s0AhGWYS54LnHI1auLsOsV5dytyjGWwOzzCa+f99Hu0FZWvvtLGm8rs/FfPD8g0ZFGRtZ6V6SvRz3ayC89AskFVRx34pxmPecnmuPB689ZjyVOpLBXJXDdCLKV8XgixJxa+MrovfEQZZikU1x56oWdFejkr0su5u/idALhDbs4MN3PKdTIqzo04g07+/qY7+NxgF33uIF++YxXpJsF/Xvge/oi2HP7sTT8j3zyabEBKeXV+xJEghr4G6Gkane0OTGdhTWjuQrHoVEu2Qv6a5HjM18jysuLaYsEJcYL7EsAOIcS+BJ/g7cAeIcRRtDCWu8brQ6FQXHukG9LZlLM5HsFLSkmPr5siSzFP352HLxSmNN1Aq7slLsIF5sIkEQY4PdzE442PJhmCrc5Yi1k//nK7zt2F0VGHwXEMo6OOfEcDusg09sit+WPiMW+Kx2NWKC6HBSfEMfel6GNsXStKfBWKBUsgHOFcv5v1BZkIISiyaFGgSrK0vTm/30+eKY+/Xvk/OeVqJMeUm9LHKVcjDn8vDkcvLzl+B2gz8JXpq1iXWUEV6VR73awYOIfJcQy9p3vqgemMWijIhL1dbMtVPGbFjLDghFihUFyfDPmCPPybUzQ5Rtj7nmrW2TPHPS7PnM/7lt7H+ybo58xwc0pZWIY5M3KaMyOn+Xm0zCbDVGeGqTJlUe0LUOX3kym1PepwxhL0y28YFd2iajDOjqGOQqGEWKFQXBM8dug0x7pcAHz6V4388L7N2NKmb7Uv/EMY+47zRV86fxlcQpP7LCf0YRrMZtpMqf049Xpeslp4yapZ0ApguTGflaZKakpv5I1rt7LKtgr9Yk2hqJgzlBArFIprgk+8sZwT3fUMB0L8RfUSss2TfD1FQhicZzD01mF01GF0HMcwdDZebQOqIRpFG4Z0gpNmMw1mE/XpuZwwCVwk+x5L4Hywj/PBP3DozB/4+hn40LoPsfuGyw5DcFlIKRnq9dJz3kVWXhrFq2xJ9d1tQ3SedVJUlo19eSZGk/phcL2hhHgGcTgcdHZ2kpmZSWZm8rKay+ViZESLH5uVlUVGRkZS/fDwME6nZnOWmZmZ4q82PDyMy+WasH5oaIhAQIvYk52djclkSqp3Op0Eo8E5xqsfHBwkFNJSqOXk5KT4NA6MOIlENB9JW3o2Bn3yEf3Dg8josl5uhg3dmNB7/SODGKLB53MzbehEcn2fawDp9SNcZuwkRxCKRCL09/cDWlad/Pz8Cet1Oh15eclxkUOhEIODgwDo9Xpyc3PHrXc6nRgMBkpKkv1U/X4//f39+P1+srKyUvr3+XwMDGiZcNLS0lL693q98fNbLBZycpLjAXs8nvh7b7VasdmSv4hHRkbi731GRgZZWckGQcPDwwwPa8kZJrv3/H4/drs95d5zOp14PJpPb9CX+iXvGh7C7/dF+88mzTzm3nM5CQQ046asLBtmU7JRlHNokGAwdm/mYDKOufec/YRCITKAf7xzDU5/hO2rRg2w+gf7kG4Hemcz9uFGMpwNGPrq0YW0MfeQTxgdUEAhfegZ9eWNmG1cyn4joZwK8nMquGfFjbzPkoOUknbvJZpcjTR3NtLmbqXd086AuT9pbJW5lXR2asFIhBAUFxfz6MuPEgwHqbZXc9fyuwgMadem1+spLCxMah8Khejt7QXAYDCkvDfBYJBXXjhB/cFLgJ6qG1YmCbHf76f+1WaaX+tGSD0GmU7+kgwKy7IoKsvCVpJGSO9BCHFN33sw/vde4r1ns9mwWq1J9YODg/HAKjk5OSn+vwMDA/EgG7m5uSnfi7HPLUBeXl6Kf7zD4Yh/L9rt9pRYB729vfHvxYKCAmYLJcQKhWJecPmDKfGhNxZnYRSRhJluHeGeFqSnD4As+jEwfhQtKfQEczbgL9hIyF5L0F5LOGsFIwMOkBIZChJou4S89Ed0uXZK165nadHb2GTQ4gD5IwEGTP3UdRzjvDzLaXcTVfYqcI+eIxwJc+jiIdxBN78+/2tuK70NkRAvKxQJ0ePpoSS9BCEE4XCEwW4P/R0jDHZ5wWvm3Z9IjjuUUzz6w7PnvCvluvo6RgcgIxLHxWEcF4c5+YcOIiKEPiNIXkkGhctslFdIClZkkZa+AOJuK+IoIVYoFHNKREq+8VobL7b28cx7azC4uzD2NWIYPE2O6yiWwROIyGhIVD15hMaJOxW2FhK0b8adXkMwp4JQ9ip09lL0Oj1SSiKd7YSO/jfe042EzrcQbr+A6O9HJ0eDhuiWLMezoQr9ilXoy1ZSVVXL0tzV5BSnU7A8i3A4TI+7J35861Ar7qAmjHaLncL0Qnrd2oxXSsmJtmYe//kelnsrWOpdg2Uwm7CI+Q7rMITS8bmDSUKZU5SOOd1AbnE2ZasLUvygy2ryycgxM9jpx92V+noGvSG6W5z0nBuh8UB/tE8rhWVZ5Cwxk5Ynyc5f2JGwrneElNPIcamYkMQQl0eOHGHz5s1X1E9nZycOhwO73T5uvd/vn9sQl61uZI8PUTgzlqLBYHDqEJfR84ny+QluH3v9xy4zaVrGAAAgAElEQVRNx/B6vQs6tN90xt97wcVglztpljaTiOAIP/ztC/g6/swm3Tm2GFrIkVNHqJJ6M8H8anw51cjirQTttUTSkxPcB48fIXjyOKHTjYSaG5GuaUS+Gos5DVasIeuBXZS++/aU6mAkyJnBMzQ4Ggj6wtxheTs95110tw7Rc96FZ2jqhA6597i54Yb1lGWXxbdophuExO8J0nt+mO62IXraXHS3DeF3h6ZsZzDpKFiexQ3vKKN07dWlSlwMn4OJOHr0KFu2bIk93SKlPDoTY1IzYoVCMTvICHrnWYzRQBlGxzH0g6f5n0iI/SabYB4QylxBqKCWoH0TQXstodx1oDPi9/sx6fVIv4+xCQA939tL6FTDpEPSlS7FUL6GcE8X4ZbTEB6zzO33wenUPmREcuH7P8NpLGYgmIm4uJKRzhF+Jo9N+TIMmwboyTxPb8YFejMu4Oi4RPjnITKMGVTlV1Ftr6baXk1VfhU5aZOLpNlqZOn6XJau1/aCpZQMObz0tLnoaR2iu81Ff/sIkUjyCxsKROg86xw368SpVzrJLU4nf2kGBqMyBJsPlBArFIoZQfgGkkTX4DiOLjg8ZbuIMZOQfSPB6L5u0L4RmZYaqCPw+st4fvQ9RlpOk/au95P+4N8m1RsqKpOEWGRkYaioHH2srUSXNZq9Sfp9hM6dIdR8klDTSULNJ4k4ekAIDGvWjR4nJd9/9BWGB7MBT/QxPkazHvvyTMwlEZzZXZwzn6B+5CjnnOeIjEkIMRIc4bWu13it67V42bLMZXFh3lq4lVU5qyZ97YQQ2Aqs2AqsrH1DEQDBQBjHxWF6Wl30tGkz9ZFBP0JAwfJkYyufO8jvvq/5XRvMeh782pvRG1SO47lGCbFCobh8wgEMg00Ye49Fw0PWYRi+MGUzKXSEbGs1Y6oCbbYbzl4JsSVaj5tQcxNhx2ukbX97cttgkEhUaEOnG1P6Nm59IzIYwFixAUPFBnSlSydd7hXmNIyV1Rgrq3H1eelrc9HX1k8O/eRZR5fmhRBkWkIMD47tAHKL0zXrZdmP6dVfYq8oI72mmrTqavQZW4BoHO2gm8a+Ro52HeWU8xQNjgYGfAMpY7o4fJGLwxf5r9b/4v1r3s8XbvpCUr035MVimHxZ1WjSU7LKRkmC9fXIoI+BLjdGc/KMtzfBOCy3yJoiwpeaB2h4sT1upV2wQoXwnA2UECsUismRMh6POR6Tuf8EIjx1PGafMZcTYi1rKt9CpLCWUH410qi5sMhwmPCl84Re+yWh5kZCpxsJX2jVUgwaTZhvvTMpDaehonJ0SM7BlH1VU+0NmGpvmOJSJAFvGLM1+auv4/QgJ1/S3JTk2uKUdgWlFvq7Bsh2X8ReaGTVg++lYHkmpjStH8e/fIO+3/+C/t9DP4AQmFetJK2mBkv0sXXlZqpsVVgsFqSUdIx0cKLvBA2OBhocDZwaOEUoMrrfW22vThnHp37/Kc4OnqXGXsPfbvxbym3lk15vjIycNDJyUu09zOlG1t5YRE+bi8Ly7JT6jtODnG/o43yDZrWO0AzBilfaKCzLorAsi9yidIROhfq8GpQQKxSKJETIg7H7hLa83FuHse8Yek/PlO2kzkQob0N0X3cTP+4r4UtHvUgEd/XZeWxpPuHDRwk1NxJsPknoTBN4J1jmDQYIt57FsHY0p7E+v4C0L+7Bsq4KXW7e+O3GEPCFGOh0M9Dppj/6Ny3DyFsfrEw6Lrdk1L/V5UjNkHTj/9jKGx8QCCGQkQhijJ+8t75+zIsh8Z89h//sOYae/zEAOqsV04ZK0jfVYqmpoaimmiVlb+NtZW8DwB/20zzQTIOjgROOE2wq2JTUZURGOOE4gSvg4sCFA3yq9lNJ9eFImF+1/YpqezXLMpdNy/ircEUWhX+pvcZj95UBetrGuFNJGOzyMNjl4dTL2g8XU5qeghVZFJVna/2VZWHJNKX0pZgYJcQKxWImEtFS/EXT/eW0vo59sAkhx/fVTSScsTRuTBUsqCWUux6iGY6k30ek/o+859zrrB28SOWBSzjdqUuxKeh06JeXY1hbCeMkYTdsuQndBN4DkYgWoWqgcyQuusP9vpTjAt4QQX84aZk2t9jKsspc8kozMJpT90gTl2zHijBA0Rf+Hs+x43gb6vHW1+M/fSbFECzi8eD782F8fz4cLzMuWRKfMVs21lBdsY4ae82419fr6SUYdevKMeewJHNJUn3bUBuff/nzAGSbs+OGYDX5NWywbyDLNPmysm6cWe0t96+lq2WInjbNEGygY4SxjjYBX5j25kHam0fX7rPsForKsigsy6aoPIu80gy19zwJSogVisWEdxA6jiYkuT8CPme8eiIHs4jBSih/Y3xfN2jfhLSkutq5//1pgiePE249y83hMDdPMRyRm4+hohJjRSWGig0YVlcgrNNznfK4AtGZ7ggDnW4GuzyEQ5Ep2xnT9LidfmyF1oQyA294l7bMO9jlnqjphJhWrMC0YgW2974HiIpuYyPeek2YvcfrCTkcKe2C7e0E29txvfACAMJkovDvHyXn/e9PObYovYhX73+VFmcLvZ7elBlvQ9+oodqQf4iXO17m5Y6X42Xl2eVx6+waew0rbSsx6CaXAFuhFVuhlXVv1JbrA74Q7Wf7GOzwRd2nXHhdqS5bLocXl8PLmT9rKyl6g4633L+G9TeP7xq42FFCrFBcr4RD0HsqLrjmS3+GgXPTahqyrRmd7do3EbaNJriPuIYInWxEX+hGv2xFUrtgYz3hs6nZjwAwmzGsqohaMUcNquwF01pClRFJX/sIPReHcPX66O9w4xsJTtlO6AS2Agu5JenklqSTV5JBRq55Wue8GnRWK9atW7Fu3aqNX0pC3d0MHT5M+FQT3vp6fI2NyECyiMlAAGNxqlj1PPEE+vx8LDU1rNm0ibW5a1OOyUvL482lb6ahr4Ehf6oPdetQK61Drfzs3M8AsBgsbMjfQHV+NVV2TZzzLfkp7RIxpRkoXpVNeVVR/LqGB3xR9ykXPeeH6L04TCSUPG0OhyJk5qaucPzxuTOkZ5s1Y7CV2ej1i3PWrIRYobheGO6Oim50ptt5DIKje7ATfsVZcuPp/gaN6+kLryZ7WarBEoDnP5/B+3/3ac3u/yusf/kQABcGPRxocfCBtZWEzzQBoF+6Iu42ZKioRF+2CmGY+itHSomMSHRjvpRf3n+WUGDyGa81y6QJbqkmvDmF6eiN8//lLoTAWFxM+vbtWN71LkATXd/pM6Oz5vp6gpcuYamuSmobcbsZ+P7/07YRhGDN4T+jTzBik+EwCMEtS2/hlqW3IKXk4vDFuBFYQ18DZwbOEJLJgT+8IS+Huw9zuHt0qbwkvYRqezW3LL2Fd5S/Y1rXlZVnISvPwuotWpztcDBCX/tIPOhIT9sQrn5fisV1wBei4XftILW0zg8+/ZYkIY5EJAIWhSGYEmKFYiES9EFXvSa6HUc04R26NHU7nQGKqqB0y2iu3dzyeIL74AUXkUmWZvWlS+P/h5o1F6JjnUN89jencPlD5JffxN1feTOGNevRZV6+q8vIgI9Xf9pCxY1FLKscNcgSOkFucTq9F0b9kvVGHbnFo6KbW5KOJWPhGAkJkwlL1QYsVRvgQx8EIOx0os9Otl72NjZqIgyYV61CPyZxguf112n/xCexVFfH95pLq6pYvvKdvHOl5j7lDXlp6m/iRN8J6h31NDga6BnHAK/T3UmnuxOr0ZoixC3OFkx6E3n6yQ3l9EZd3KI6hm8kiNmSLDe9513xgC65pRlxC/QY7U0D/PaZRq2vFVkURo3Brsc42kqIFYprHSlhsA3aj47OeLtPQGTqpVmySjWxXbIVv70a84obwDi98H4R5yA6W3KkJ0PFBjAYMKxcG7do/mVzDy6/Ntt6ukPw5ltrybNeviD2t4/w8vPnCHhD9He4k4QYoHRNDmmZBuxLs8grSScr33LdzZb0Y7IfAZjLyyn+31/CW1+Poagopd7b0EBkeBj3K6/gfuWVeLlp+XIsG2viLlSb1lRRWziacKLb3c2JvhOccGjifKr/FL6wZtxWlV+Vcp5/OvpP/L799+SYc9jzlj3cVHLTtK8rLSNVPHOK07ntQxV0tw2RlZe6bN1z3oXfE+Ji4wAXG0cN/WyFVs0QrDybwjLtXhi7erLQUEKsUFxr+FzQWTe6xNx+GDz9U7czWKBkU1x4WbIFskb3GyNe77RF2P/HQ4x89X+R+cj/wvTGW+LlusJicn/6IiIh1eEjt6yiw+XlgtPL1++uvCIRBhjq8xLwaoLe1z6S4ie8aksBS/3Zcxpz/VrAkJ+PbccObDt2jFvvPzv+vn/gwgUCFy4w9PNfACDS0kjbUKnNmqtryNtYw/bl29m+fDugxdE+O3iWBkdDishKKePGYIP+QQrTC1Pqn/jzE1TkVlCdX025rTwl1elY0rPNrH9TCevfNL4BV3/H+Cszzh4Pzh4PzX/qBkbjaGtBR7IpLM8iPXth3SOLRoiFEPuBbcA+KeXsZvpWKKZLJAyO08l7u45mJgzCnEjeKk1wSzdrfwsrQX91y3ZSSuSvn2Pkub0gJcP/+AWynvw3jNFgGkIIGJNv2KTX8eRd6xkJhCnNuvIkIeUb7YwM+mmr76P2rdPzg1VAyVNfxf7JT+A9Xo+3oUEzBGtqgmDyion0+fAeOYr3yGieAkNRUZL7VMX69azPWz/2FLiDbipyKzjhOIFEsiJrRVL9BdcFftj8w/jzDGMGG/I3xC20q+xV5I4TtnQy3vrRSlx95XS3uuJ7zX2XJo6j3Xl21Po/I9esiXLUhcq+7NqOoz2jQiyE+Cya2G0BBoA64Akp5fGZPM8VjGublPLe6P+DQoi9UsrW+RyTYpEy4oju6UaFt+MYBKaOx0xa9qjgxsTXenlfbFMhw2E83/oa8hf742U6e2FSfGZPMMzxziHeuDz53NlpRrLTrn7vrurWUlZvKVABIS4DIQSmpUsxLV1K9ju1fd2I34/v1Km4EZivvoFgZ2dK21B3N8Pd3Qz/938DkLdrFwWf+mTKcRmmDPZu30tERrg4cDFltpvoOgVaHO0/df2JP3X9KV62NHOpFkc7v5oaew1rctZgnOSHoxCCbLuVbPtoHO1QIIzj0ojm1xyNpT0ymBrhbWTAz7mBXs4d1VJUbrl7Rdw97VpkRoRYCLEJ+BzQCuwDnoxWlQOfF0K0SCk/NxPnuhKklAcTnh5B+5GgUMwuoYC2l5sovIPnp24ndFBQmbDEvFWb/Y4TSGKmiHg8dHzms/hefDFeZqisIfPxr8aFuM8d4NO/auRs/whfe1tlihhfDqFAmFMvd7L+TSUYTKMzFSGEEuEZQGc2Y920Ceum0ehcwd5efNEZs/d4Pd6TJ5He5Chilo2pwUQufnQn6ASWmhpy7r+fQmthyjEb8jfwidpPxC21+32pWymXhi9xafgSL7RqPtNmvZl1ueviSS5q7DUUWgsnXQkxmPQUr8ymeOXoj8ORQT8952MW2i56z7sIBZOt6xMNx65FrlqIhRBlQLmUMtUDHQ4B3xZClAkh7pFS/uRqz5dw3h3AQSmlM6GsHNgNHAC2A7vH1NuAusQyhWJGkBKG2pP3dbvqYRrxmMkoHN3TXbIVijeCOWPqdjNEqK+PSw99DN/Jk/Ey0y3byPjMF5P2gp/84zma+0YA+PyBZp6/fwv56Zcvmj53kJf3n2Wwy8NQn4+bd6waN6qTYmYxFhRg3LaNzG3bAJChEP5z5zRRjs6cLTXJQhzx+XC/9hqEQrj/8BK5f/EXyfV+P4Hz5ylbtYoHqx7U+pWSLncXDY4GzUK7r4Gm/qZ4VLAY/rCf447jHHeMLpgWWAq4c8Wd7L5h+ruHGTlmMnIKWLmpAIBwOMJAh1vLPBUNOlJUlhpH+1piJmbENinljyc7QErZJoQYEEJkSSldkx07FVEx3QnsAVYCiaJ6ANgupWwVQrQCR6PHxHi/2h9WzAgBt+anmyi8I1PHY0Zv0oQ2LrxbIHtp3H1orvG3tHBp5y6CHR2jhXd/gIyPfyIllOPDb1nF6b4Rekf8fOKNZVckwq4+L3987iyeIS2QRXfLEJ1nnSy5ymT1istHGAykVVSQVlFBzgfuG/cYX1MThDQDOtOKFZpVd8Is2tfQwIUP/3/orFbSqqqi+83VFNTUcFfZXdxVdhcAgXCA5oHmJPepjpGOlPP1envHDUbyauerODwOqu3VLM9aPqkhmF6vw74sE/uyTDbcMuFh1xRXLcRSjmbGFkLcHi17UQjxIHAv2qz0uJQy9dW9svM5gSeFEHsSy6MzZGJ7v1LKOiFEbnR/+KAQYhvwXPTYcrVHrJg2kQj0n0Pf9gr01GvC29sIcupwiuSsiO7pRme7RRvAcG1YdLr//Gfa//bjRFzR38Y6Hel/81m8m986bjzlfKuJp++upGfEz03LLn9Z2nFxmFd+fI6gLxqDWcCm7cuUCF/DWDZupPzXv9KSWoxjPxhLdhHxePC8/jqe11+P1xlLS+NGYJbqajasX0+1vZoPrtN8pvu9/UnZp070ncAT8oybdWr/6f0cvKjtMP7DG/+Be1bfMwtXO3/MtNX0Q8CD0T3jvcCdwH3AXBhrbUXbo05kAKiNzqK/DQxE9x92j3NsEkKIx4HHLmcAPp8Przc1c8t08Pv9BINB/P7xlzInKp8tRDCICIWRwWn4qk6DYGga4Qij55NzfK0xYq+/d6ATXWcdus6j6LqOoeusQ/iHmGr+J03pRIo3ESnZrD2KN0H6mHjMwQgEr+weuVp8vtEECCO//jV9jz0et6wVFgv2Pf+IZ8UmAj1efD4fPe4ARRnJPxpK0w2Uphsu+35sb3Jy/LftRMLat7neINj8jmUUlWddVl+BQGpc45kmELsPrvCzPBmJ78GCoagIc9R/2ev1Jl1DMBRCb7cTHi+OdkcHwY4OXL/6lVZgNGKqWIu5qgpzVTXm6ipuKL6BN+S/AdZp2aPaXG3kpuWmvPb1jtHsVqszVifVSynZ+budLMlYwobcDVTlVVGeXT5hHO2reQ9m6/2baSH+kZTSJYT4NvBVKeWh6B7yXFBO8jI10ed5UsongecvpzMp5ePA41MdJ4TYjGYARlpaGhbL9Pw0x2I2mzEajZP6SM6l/6Q0hpCGcFI+2KvFOEVfsfOJufQTjQQxDpzB2NtA9oXXyXKdxjB0fhoNBdgrkgyqhH0tep2ea9dJQrtH+/d9m76nn46X6e35LP23f8eyoZLwBRc6fZCnXrvI79v6eea9NSzPsU7S4+RIKWl+rZuTfxhdhkxLN/Kme1eRUzy95A5jme3PgckYwmw2X/FneSpmq9+5JHYNloceouihhwh2dye7T508mfqDOhgkcOIkgRMnGUZzddLn5cXdp9LfdDNVlamBREKREPetvY8TfSc45zzH+sL16HWjn7KOkQ7q++qp76vnhfOaIZjFYKEyr1KLoZ1fQ7W9Grt19Efxlb4HaeNkBJsJZlqIbxBCPAQgpXwkKsK7ge/M8HnGoxWoHVNmA1rm4NyKBYLO3Yup9zimnnqMvfUYHY3owtP4lWvNI1xci375GzThLamFtGvbEnMsMhSi+4uP4dw/6p5kWrWSZXv3YiwtjZc9fewCBy9pjgWfeKGR776vhhzL5e8HRyKSuv++QNvxvnhZVn4ab3r/6gUXcEExOcaiIox3FZF111sBkMFgNI728dE42hcuprQL9/cz8uKLjLz4IpFhF5bK5DzRoYEB9DYbu2p2TXjuBkdDSpk35OVIzxGO9ByJlxWnF1Ntr+ZDqz/ERsvGK73UWWFGhTgqvpsS9o3L0Zar54LDwNjQM7lEZ6uKRUjIh6nvFMbe+qjwNmBwd03dTmeAoupkg6qcMgI+34KdzYRH3PR+4pN4X301XmZ9wxtY8o1/QZ+V/IPi3eV2Xu5y4gtF2FSSRYZp/K+JUDjIkKefcDhEhOT9chmRuJ1+0leE2bBCm/kaTHoycswMyx6Gr9BvQUYkwju7hm1hfQS/s59+z8yvbUQiEXSz6IY2F0z7GvQ6qK3VHmg/BCNeL9LrJeLxaJHeEnI29y5bRn9TU1IX/tZWpN+PzmLBUFw8bi7qZXIZ/7HpP/CGvHhCHrwhb4qFdiIhR4imoaYJ6/V6PTabjYKCgqmvcYaY6YAen0HzIwZASnloJvufDCnl80KIPUIIm5TSKYSoBY5IKevmagyKeURK9K6LmGKi62jA2NeMGJNxZjxCGcUEC2oYsK7EVH4z9qrbpx0KciEQ7Onl0kMP4U/4kst+97so/tKXEKbUme6anHS+vL2CZscID24ZP8JVKBxkcKQXKSLoDAJdQm4nCURCETIS0t7phECnXxguSnq9btaM2K+HaGFXeg3CYECfmQmZmfEyGYloxpCRCIzdupISY/FoFjDd2HtVSmQ4jNDpsBosWI1W8siLVkkiMkKEiPY3wbAyTT/58nI4HKa/vx+r1UpGxty4Ec6Gsda+xIKZcFka01/MfQlghxBiX4Jf8HZgjxDiKLAZmHg9Q7GgEYFhTL0NGHsbMPXWY+ypR++fepoV0acRLNhAoKCGYEENgYIaIunaL98+hwO73X5diTBA379+I0mE8//6Y+R//OPxL9RQOIJhTND8N6/I480rJs6yM+TpR4oI6enpFBUWo9frEUIQDkUYGfAlhSFMSzdiyTQCVy9CUkYQU8QwvlrCwTB6ox6jWc2Ix2OuriHi9xNoa0OGQgijibS1a5Lqwx4PgdZWzYcf0KWlobNa0VksCKsVXYKwR2QEf9hPIBwg05g54fillPT19dHX10dvb++CFeI9wE4hRGIkq13Ax2bqBDH3JUajdyXWtaLE9/ojEsYweE4T3N56TL31GAZbEdOIxxzMLiNYqAluoKCGUO5qbel5EVH4yCP4Gk/hO3OG4n94HNv73heve+7IJb77ynl+tPNGsi3TN8wLh0PoDIKiwmIM0RzDwUAY96AfKUffF2u2OSX9nUIxHXRmM+a1azVPimDqypb0eBKeSG25O6FMGIzorBZ0FivCaiHNYsFisBCJTOx2KIQgNzeXvr6+ObHQjzHTn5BdpBpMSWZQiBXXPzpPX3yma+qtx+g4gS4hwf1ERExZBKKiGyyoJmCvQqalppVbbOjS01ny7//G8KlT2G69NV7+jUNn+dqBMwDs+v4RvvfADdPuM0IEHTr0em3WGPCGcA+NWskKIUi3mWdlVqlYPAghtO2TcbZQhNmMPjubiMczrpulDAUJu4KEY37yCHRpZvRFRegmmenGZsuJPyhnm5kW4ifGRtkSQrxvooMVCiIBjMNnMA2dwNRzHNPZJgy+1OD0Y5FCTzB3TdJsN5y9XIvTvIiRUuJ+9VUybr45qdxYUIAlYW8OYFneqFvSkDfEsG/q/fSxCCHwjQTxjozOHnQ6QUZuGnrD4n4vFLOLPjNT23NGs9KOeL1EPF4i3qghWMrMVxLx+dBPsaw+H/v4My3EKa5CU4W/VCwipETv69JEd+gkJtdJjK5mhJw62EfYah/d1y2sIZhfiTReuX/r9YgMBOj6whcZ+vnPKXr8MXI+8IFJj3/3xlLaB7283jbAtz5YS4bZQG//5QdT0RtHv7j0Bh0ZOeYFn6hdMbc4nU5aW1uprR27oDo9hNGI3miMewBIKZF+f9w6O+LxaH7NOt3cximYJjMtxN8RQhwAnp3v1IfzgcPhoHOcVGPTbet0asZGJpOJ7OzkIOUejwdXdInFarWSnp4cDGFoaGjG9jSsVitjJW7IMxyPjpVlzcBkSF4q6h8eHHcpR+/twOSsw+A6i/S0EAkO4QNGPXdTwxvm6Dzo7au1mW50xuvwGYlE+8/Ly0sytohEIvT3p2Z7uVxir39hYWF8yRU0K8qenh78fj9Wq5XCwuTsMz6fj4GBmUnopdfrJ+0/LS2N3Nzk8JLDw8MMDw8z+Nx+nC+9BDk5DHzjXykuKCCtoiJ+nN/vTwmG8Z41Ft61uhRXfy8uwO0KAcl7xW7PCB6PlqTdak0n3Zq8rKczgDlDRygQwZJpICLDREJhrgSdTodOl7ycHQ6H4veW0OkYa3wdCgenlb45kWPHj/HUU18lJ0d7LeuO1XH8+DE+8pEH+dd/+iZPfe2rfP7Rz1FbW8sPfvADnE4n999/P21tbezfv58dOzRPyQMHDnDnnXfyzW9+kwcffHDK80opJ51xGY3GpHopJcFY9DMhUoLiRCIRQqHxVzJ+8pOfcPvtt2OzTb0989RTT/Hoo4+yd+9edu7cGS+PRCLs3LmTZ555hieeeIJHHnkkqV04HCYcvrL3GqC1tZWnnnqK8vLypPO2trayZ88e7rjjDn7729/yla98hby8vKTPJUAoFOLJJ5/k0UcfTX0PdDr++uGHeeaZZ/jyl7/MZz79aWQwSCgUGnfMer0+pf+BgYF4NK2xn7uZZKaF+A4p5ZAQ4g4hxBNAP7BvJq2mFQsL48hZMi49O+kxYXMhwYxVBMUKQoVrCdRuxmQZs4fj749bRyrGx/aud+I9cgR/SwsZt92KedWqeF2LY4RfHm/nf25biy7hi14IkSJsUzHeDy5TmgHT7AQdmnF++tOf8Dcf/xte+K9fsWnjaJrArz/9NerqNG/Hz/zdZ6k7dpTW1lbKy8uRUvLNb36Tu+++m/Ly0by2W7Zs4Z577pmWCM8VTqeT73znOzz66KM0NTVNS4hj19TSkryoefDgQX784x9TVlbGZz7zmRkd57Fjx7jrrrv4zW9+w6ZNm5J+XG/fvp0DBw6wfPlyli9fzo033sjZs2cnHXtra3LU4kOHDqWMXZjNc7r3O11mdP0oIbHDYTRfhc8DR4UQ/yaEWDGT51IsDIIZK5OeS52FQFYl7uJ3M7T67+jb9C0Gar7G8MqP4bXdTii9XMtQpLhsRFoahZ97hLwHHiD/Yx9DRK2Zj5wf5JEfn+BQcx/7Xmq9qi+ioTutfHoAACAASURBVF4Pv/t+MzJ87X2ZTQen08kHP/xBPvN3n0kSYYBPf+rvuOOObfHnu3btoq6ujmPHtPhEd9xxBzabjb1798aPee65564pEQaw2WyXLZoHDx7knnvuSRIzp9PJ888/T1lZGbfffvtMD5P777+fj3zkI2zalPw+PP+8Fo04JrCbNm1icHCQgwcPpvSROPa2trakscdEeDbGPtPMRkCPlWh+vgeBe2NBPYQQnxVCRKSUX5vJc15L2O12SkpKrrqP8bBarZPG2B27lH21SNzJ/VszJzhSIy9zggw6mdlklb4ZX/oawrkbCaWvwCj0jOcoI70+RFbauEEm8vIm9mnV6XQTvm6Xi91uT1me0uv1lJSU4PV6x42slZaWdtXv+2RM1L+/tRVzeTmZmZlkxgyxSkpgTJjApgYXHQFtuvqz0x4eelsuJbbxfaV7gy4G3cnvfbo1g3RrBl0tQ/zppy2EgxFGBnXkFGnbI3r97LonJfYvx8l4ZdBP3+3qDy+9BMCO970fgyG13Ucf3Ek4qC1bbtu2DZvNxjPPPMPevXupq6sjNzeXffv2xcX44MGD7E8IGbp79262bt3Ks88+y65du9i2bVtS/zEf3Lq6Og4ePEh5eTnPPvss5eXl7NkzmlDO6XTyxBNPTHoteXl5PPzww5jG+bzEMJlMk9bH0Ov1rFq1Kr4iELuWRx55hJUrV/Lkk0/G3dQSGR4envY4Ezl48CBtbW0MDw/z8Y9/nIMHD7J792527tzJ4cOH4yIcWy7Ozc2lvr6et771rUn9GAyGpLHHrvULX/gCn//85+NjT3wNLscPejaXoxOZ6U/Qk8B+YEtiesQoh9BmytetECvGQegZWvc5gsHglEkfFNNn8EfP0v2lL1H4+c+R+8EPTnrsZ+9cy6UBD3UXBvneAzdMKMJTISOScEgTwnBYLsidgiNHDwMkLS8DHDtWx7e/820GBgf4yF9+hO133gnAzp07efLJJ9m7dy979+7lwIEDrFy5kn379rFt27akfurq6qirq2PPnj1s27aNI0fGj6578OBB9uzZw4EDBwD46Ec/miTmoM1qE4V5NmltbWXz5s3A6Gz0+eefZ/PmzXFhHvuD4mrHWVdXR3l5efwHTV1dHZs3b2bbtm20tramLKfbbLZx7UCuZuzXEjNt2nivlPK+cUQYNIvquYo7rVBcl8hIhN6vfZ3uxx+HcJieL3+FkegsbyJ0OsFT99bw7INbWF04+crGZJSstlF75zKs2Say8tLGDQP5L787y+rHfsPqx37Dv/wudU/vid80x+ufeaUtpf7vf3EyXv+jI5dS6j/9fEO8/hcNl28YubJM2yoZu5+4aZNmrfvii4e4/bY74uW7dmnxgfbt28fAwADl5eXs2LGDPXv2sHfv3ng9QG1tLQP/f3tnHt5Gde7/75H3xHZkOfuecUJIIaRIDltvWRKJnS4gxyy9t4XGFtuva2ph6C30XpogA11u2ZzA7XIvbbGV0BUKkoFelgK2FMKWkOBJSMjq2B6SOLYsW+f3h2bGWq3VlhS/n+fRY+ucM2feMyPpO+/Z3p4eGAwGiKIYUQAkSUJNTY0qvKIoQpKkjIqF3W5XHyoUexwOB+rr6+FwOMbMtkCx1ev1EARB7SVQJk4qSJKEqqqq0CoStt3pdKqvQO8/06Q76EPUpUry+PGmdJ6PICYSvsFBHLyrcSS+K4DiZctQvGyZ+l46OYgXPjiE2pXzg44tLsjD1NLUx96r9NOx4MxKHD25P3bhLGTVKr/Ibnl2M9Z9/wdBeVO0U6CrCO6KFAQBer0eFotFFU+LxQKTyaR6toG4XC5YrVYYDIawGciAf0y5urpaFSGr1RpVsJPp8k2G9vb2oHrq6uqwadMm1d7GxsaoxyZrp16vj3icTqeDTqdTvVuFnp4eVFdXp2S73W6H2WxWu6YlSYLT6cwKjzllIWaMnQ2gM9bMaMbYtZzzLamejyAmIsOShH133on+DpeaVnrxxZjz8EPQyEvZ9vWcxNd/9TbErj74OHDDOfOjVRcTz8khvOPYixWr56G4NHhIIb8wD4i90VlWIggCHnvkcdz9w0asXrVa9YRHw2KxwGq1qkuWFA+strY2qJzizdlsNphMJthstjAhVrpklfKKhxzKWHVNK+dXHgQkSVLHQZW/tbW10Gq1cXnrydppNBqh0+ngdrvVtcM9PT3quLzVaoUkSdBqtXC73aiurg5bY5yI7W63W71/gbaP1xhwLFIWYs75VnkiVi+AllBBZoythj884fgMeBDEKcbgvn3YV2/BYMCs0Iobb8CMu+9WZ0YDwM+duyB2+Sda/ehP7+OLS6ZibkXim56cPD6IN/8s4kSPB8d6BnDJTUv94hsH37pkCb51yZKo+Y2Xn47Gy0+Pmn//l87E/V86M2r+T81n4Wc1qcWSXbu2DqtWrcbdP2yErkKHKVr/RMfPpM/w+9/9Iaz8mjVrwpb1WCyWMJEF/MturFYrOjs7IwpUfX09amtrUVNTA5vNho6OjqBZ2OlAkiRs3OiPvWO321FfX68Kr+KFtra2Bj0EKMu0zGYzzGYz3G43rFYrAKCjoyPpjTZGw+FwwGq1wmQyQZIktLW1qXYqeQaDAS6XK+waJWp74Dr/Sy+9FJdeeikaGhqg1+uDHgYyBUvXmip5K8tNAKYAUDr4tfCvIz5l95pmjBkgxzzu6OhQJw4kyoEDB9ClRP+JQKTNGMYSLvaBHx4Am5GexaHxTNZSzseEyaOWGyuU6x9tBnS0WdNjSf+772LfbbdjOGCiyvQf/AC6W24O2xiizzOENc3/xK4jJ/DTNStw9VnB7YjH/u1vHMCrLbvgHRjZ8OCCa6swZ2nwrPhD0l7k52uwZPHSZJuWFKdS9CVRFGEymcJEfqxRumiTJRcjSIV63oonHE2Et8vRypYFDPsA/qGHgC7yas65C2kgbWPE8vjwZsbYIgACgJ4ok7YIgoiD421t2P/9deDyzj6ssBCzbQ+g/IorIpafXJSP//7GSnzaexKGBYl3uYnvdOEfv9upzozW5DGcc82iMBEm0oMy0YgYe/R6vfoAEjhJLF074qVK2hcAcs53AwifDkkQRNz0/M//4vD69epuYnlaLeY+9igmyU/vPh/HhweP4cw5wevHZ5QXY0Z54r0Y29r24TX7LnWryMKSfHzBXIWpc5OfZU1Ex+l0YsOGDTAajepY6Hidd6KKv9lshtPphM/nw9SpUwFkz9ImChRKEFkE9/lwxNaEnt/8Rk0rmDcP8zY2o2jRIgDAgHcY32t5B84Pj+DpunOxcmHyE058Po7X7bvw7kufqmklZQW46MalKNPlyJ6VOYjRaERvb29GzjuRMRqNWdm1nl3WEMQExjcwgP3f/k6QCJesWIGFz/xBFWEAuOfZ9/Hce4cwOOxD3W878GlvclOYvYPDeGHj+0EiXDlnMlZetZBEmCDGkbR7xIyxawGAc76FMbaWc/5kus+RrSjRl4K2HJQ5duwYTpw4AQAoLy9HaUhg6uPHj6uL2MvKylBcXByWr0RfipQfGH1pypQpYdvaSZKkRnGJlN/b26tGcamoqEDoNJWeExJ8cnxP7eQpyA/ZBjIw+pKuVBv2xNl9ohf58laFujItNCETbo4e6wHv94AdK8I0BE/WCoyuxBhTu5Ui5Ws0mrDtMIeGhlTvQ9kuL1K+JEnIz88Pm6zl8XjQ3d0Nj8eD8vLysPpjRUfq7+9Xz19SUoKKiuAx15MnT6J7714c2fAA+Hvvqa0vM5kw+8EmnBwawjE5qldpaSm+Y1yCf+w8gqMnBvGVz89BqWZIjfo12mfP4/Fg2rRpKC0txcljg3ju8XdxePcxDOcNgDMv5p5eAf3F8zF4AsHHH/8MHs+AXP8UFBcFf/YCoyPl5eWFTaZKLD8/bBLa0PCQ2kWv0WjCNhIJjL6Ul58PhpDjh0bCbPq3aRwtP3hCYWD0IwBh35tY0ZFC80O3Tw2MnhQrupJGownbZjLV/MDoSbHyI0UnSjV/aGhI/V2JlZ+fnx/2u5JqvtfrVT97kfKPHDmiXr/p06djrEirR8wYewCAHYCywK5NjsKUFTDGhNilCGJ8Gdy7Fwca78bAzp1qmu7rX8ecn/8MmuJwz3SebhKe+vpK3HvN53Dfl85Aniax8EnS4ZPY3NSBw7tHVhqedu5MnPfVKmgKqJOMIMabdHvEAvwBZtcA/olbjDEzgOhbs4wTjDE9/Ptd0xRQIquQ/vgnDB065H/DGGbcfTd0//avav7JwfB4syvmabFiXuITfA7t/gyv/LoTnj7FCwPOuXoR5q9Ib9AQgiDiJ91C7JDjEXNA3XUrK7Yu4Zy7GWNjOld9tDWo5eXlKC8vj3psWVkZBgYGoq4jLisrSyn6UqxZmaHdpRyeoPe60tGPjxp9SckvrRh1HfHUcp0afSmUWNGVYuXn5+fHnR96HQCgqKgopehLJSUlo67fXbDu+9CInTj5djvmPPQgyuQJNZxzPPXabjz56m5svv0CzIkSrCFSd3Qgymfvgzf24dXf7VKXJ+UXaHDp2jOwaMXItTnyyTEMhETeKi+bApRF/3zFir6Uan5+itGXIkVZijefMTZq9KJE85VuUgWNRjPq8WOdH6k7eDzzI0V0Gs/8WHsbjGV3dCDpFuIexlgLAM4YU8IhbkzzOQjilILl52POT3+GwU/2oCQghOFDL36ER1/2b/Zw86/eRuutF2BKSeIRrDjn2PriXvzz2ZGNI0rKCnDVHSswY2H0h0OCIMaHtA4IyZt6bIB/HXEFgHrOeVLd0owxM2NMG5ImMMaa5bzm0HyCyHb40BCkLc+Ch3hGeaWTg0QYAL6weCoK8vzjv2XFBfD5ktsF760/iUEiXDFzEszWahJhgsgSxmJDj60A1B21GGMLOed74j1eFtd6+PemrsLIdpkA4ABg4pyLjDERgEsuA8ZY+MavgMg5dybcCIIYA3wnT2L/976PE6+8gsE9ezD9e98dtfwFVVPRZD4Lzu1H8HDNChQXJLfl4qIV07CtbR+GvD7MXqLFFbcuR/Fkig1NjB2SJEEUxYzv4ZwrpHvW9OcZY88wxtqVF/xiGTecc4lz3hShbrOcL8p/3QB0jDGj/H5jhBeJMJE19Dz9NE688goAoHvjRpx4/fWYx3z17Ll45IazkxZhAJixqBymb56BKsNUfOlbn5/QItzW5sS555+DwuJ8NN5zFxrvuQtXXHkZGu+5K9OmpZ2mpiYwxtQAEIFYLBYwxtDUFPZTmxKiKMJisWDjxo1qlKnAdLvdDovFEhZvOBtszyTp9ohfgn9MOFAATWmqeyUAMSStB4A+5HwRkWdNC4wxvSziscrfB+DeRAwcGBhAf39/IoeoeDweeL1eeDyeqPnjCfN6wYaGwQPWUKaCdyh2Pcr5+Di3VUG5/tHu4YC853OyTLr+epS8/Tb6X30NU26+GZrPf149V8cnEja+tge/WLMcJSmILgD4hn3Q5AU/Y89aWoqKBfMwOOQBwidhq3g8HgyO8jlU4HI3eaTJU2OJj3NokPw5V61aBfN1Zmzd6sb6+9cD8HtvS5edBl2FDuu+v85/Du6Dz5fYsrBIiKIInU6nTpYMnaw1lixcuBAA8PHHHwed1+l0oqWlBYIgYN26dQnbFK282+2GyWSCw+FQPWGlrMlkwgsvvABBELBw4UIYDAbs2rVr3G0fzf5AOOfgnIf9FqT6GxCNdAtxC+c86NFSnryVDgQEd1NDfl8ZoWwYsvjG/c3inN8H4L5Y5QKjLxUXFycdnaeoqAgFBQWjzowe1+hLBUPg+cNgMWYVJkLM6Evy+dg4tjMQ5fqPdg9Tjb40/+c/x4n/+7+gwA1/2XYA32/ZhsFhH6x/3IEnvmZIeG2wwqc7evDy/+7A1XeuQMXM8ChWsewvKvKisGAo5meN9fvtG+tISKFokHr0JWXDEKWeigodDHoDXG4XGNNAw7j/b4rbIIqiiJUrV8LlcgVt8jJe2yu2tbXBbDZj9+7d6jklScLmzZshCAKqq6uTtiXScbW1taivrw+MTgTAH9wCABYvXgwAqK6uRk9PD1566aWoW26Ope3R7A+EMQbGWNj3JXQjpXSRbiFuZYwthN9TVagD8FAa6hbh934D0QIY3xhiBBEnfW+/jUnV1WABX3rN5Mlh0ZOkk4MYHPY/pW/d24sDUj/m6RKPI7yr4zCcv/oQvmGOvz6yDdc1VGNSefSlK4QfURThcruw/v6RvYca77kLl19+GRwOBwDAZrPBbrdjw4YNsFgssNlssFgs0Gq1sFgsaihDq9Xqr8/lUoPTK1GWlAhAkiTB4XBAEAQ1ZrESlzfwfBs3boxad7wIggC3e6QD0Gq1wmq1oqqqKmK8ZEmS1JjF0dDpdGq8XwWn06m212KxwOl0wmq1or6+Hu3t7UHd1Eodbrd71L2vE7U9l0m3EOvhn1ClTO9k8v/pEOJ2AKFBNHWQvVGCyBY45+h+4gl0/eK/ULn2m5i+bt2o5f/1/IX4tLcfzu2H8eubz0lKhAFgyrQSaPIYfMMcw14f+o8PkhCPwu133Ibdu0W0vdSGxx55HGvX1gHwi9GzW7bgoYcehNFoBGMMNpsNZrMZNTU1kCRJFUNFiAG/cJhMJjWIvSIyZrNZFZWWlhbY7XasWbMGFRUVWLlyJYxGI+x2O2w2W9D56uvro9YdC1EU1djoikdqt9thMBhUcYskglqtNqbIReradbvdEARBtc/tdsNgMMBoNEIUxbB9DLRarbotbbpsz2XSLcQWAAaMjOUy+D3ilOGc2xljNsaYlnMuyWO+HfGM9xLEeMG9Xhz88Y/xmX0zAKD7yadQtHQpplxzzajHWS8/HXesWozy4uSHAqYvKMela8/E238RceVtZ4174IaCDZnbu8fbmPhePY89+jgA/w//lVddjs1b7Hj+uReg1WqxY/tOiKKo/vAHhio0Go1JhS185plnAECdgGQ2+/0KrVaLzs7OqOdLBiX2riiKqqfqcDjQ3NwMi8UyJkIWaK9er4cgCHA6nWGeLeBvX1VVVVpsdzpHpgjpdLqcnKmdbiF2ysuXVBhj8T3CjZRXli8BgJkxtpFzrowNmwDYGGMu+AXfkqrBBJEuhk+cwP7vfBd9r72mpk067zyUXnSR+v7IsQH8/u19+NbqxUHBDTQalrAIc87DAiQsOmsqFpxZCU2SY8wTEUEQsHZtHe6+pxGiKGJK2RTc8e3b8eCDTapYpgNJkiAIAurr/T9vyl9JklBXV6d63emgvb0dDQ0N6vu6ujps2rQJANDS0oLGxsjbOyTbNa3X6yMep9PpoNPpVM9WoaenJ2wsORnbFdEOtD8XYy6nW4g7GWPrEDy72QTgtngrkEW3SX6F5okg8SWyEO+hQ9hnuRWejz5S06Z8+cuY9Z//ASZvMbjz8HHc/Kt27Jf6MeTz4fuXLk36fCd6B+D47w9x0Y1LoZsVPCmLRDhxOlwd0Gq1EAQBG5ub8c7WrRAEIeYyGwVRFCEIQsTxW6UOk8mEuro6XH/99dDr9RBFUfXmlK7dSOeLVrcy9hzqOUuSpE4OU/7W1tZCq9WqHmY0oUq2a9poNKrjvopH2tPTo9pntVpVL9/tdqO6ujqi55qI7W63O+zBRavVhkU/ywXSLcS18E+gCmQREhBigsg1Bj76CPvqLRg6fFhNm3rHHZh65x1BHusf3t6H/ZJ/OcRjr3Tiq2fPgTCtNKy+WHTtO46/PbINfZ8N4q+PbIPZmh2TspLpHk4UzlObNd3W5kSrvRWAf4wYAHbv9vsNLzz/IgBg1SWr8NOfPwyTyQSLxQK9Xo+6ujrU1vqDyjU3N8Nms6kC2NDQgJqaGhiNRqxcuRIdHR2qINXX16Ourg6NjY0wm814++23sXr1auh0OpjNZthsNoiiCJvNFna+1tbWqHULgoCamhrU19cHjRtLkoSamhoAIwJuNpthNpvhdrtVT7ajoyPtXbgOh0OddCZJEtra2tRrpOQZDAa4XK6IY92J2q6EHgX8DzkmkwkNDQ3Q6/VBDwS5AFNiMaalMsbOjtA1HZZ2KhG4fKmjo0OdZJAoBw4cQFdXV9TgBB6PZ3yXL4l94IcHwGakZ5zR6/XGXr4kn48J4ctuxgPl+kcL4BAp6MOJ117H/m9/G74+OVBCfj5m/cd/QHvtV8OOHxr24Zu/6UD7nh48epMelyxNfEP5vR904+8b34fXI8eQzWO43LIci86aGuPIyPaHcuSTY+g92IeKWaPfg0PSXuTna7BkcfJefTKkKsTxMOwdRl5BHgqKUlvPHQmfz5e25UuSJKGlpUXt4h4v0tmGVAgU28Ax9VgiHI/927dvBwAsW7YsKN3lcgV2qVdzzhPasCoaafWIowhubzrPQRDZgrR5Mw7eex+gBF4vLcXc//oFJl9wQcTy+XkaPHqTHvt6TmLZrMT3ef7wtQN45XcfqZtpFJbk4wrLmZh7eu51xRGp09HRMe4inE0oy8HMZnNQ93ygp5wrpCzEjLFrOedb5P8jrdMwAbgs1fMQRLbAOcfRX/4SRx97XE3LnzkT85qbUbz0NACAd9iHf3Z248LTgns4SovyExZhzjne+pMI198/GalHV4Sr71yBytmJd20Tpwa5NiFpLDCbzeo4uzI2nIvXJR0e8d2MMZFz/g6ASyPkJ9dXSxBZCB8cxMF//3d89qc/q2lFp5+Oec1PoGDGDADA8QEvbn/ajVd3HcVjN+lx5fJZSZ9v2OtD22+3Y1f7yPjztPlluOqOszB5SmZ2ICOIbCIXhTeUdAhxDUY28LBGGiNOwzkIIuMMHz+OvT9owMm33lLTJn/xi5jzs58hr3RkTHX9czvw6q6jAIDvPPMOls4sQ1USk7IG+rx4/on3cGDXyEzaBWdW4tK1Z6CwOO2B0wiCyBDpGHH/GIASASlsjPhUnqhFTBy8+/fj0DduDhJhbU0N5j32aJAIA0DDZUshTPWn3XZRlfp/Ihw72o8tD7qCRPiML87GlbctJxEmiFOMdHyjN3POn4yWyRgr55wfS8N5CCIj9H/wAfbdeiuGu46qadO++11U1teFbagBABWTC/Hrm8+Ba28Pvnr23ITPd3jPMfzt0W3oPz4Sser8r1bh7EvnRzwfQRC5TTqEOFbQhXqkZ69pgsgIfW+8oYowKyjArPXrMeWaq9X8I8cGML08eJnX/MpJmF+Z+J7Ru7d14cWnPsDQoH/TBE0+g/Ebn8OS6hkptIAgiGwmHV3TVsZYd5TXMIBTK0wGMeGoXLsW2hozNGVlmPfUk6oIc87xU8dOrHr4H/jgwGcpn+fdlz/F80+8p4pw0aR8fPnbZ5MIE8QpTjo8YjeAZ6LkTQVwXRrOQRAZgzGGmT/6ESb/279h8pIlarrt7x/hiX/4O4Ru+XU7/njHFzBrSuLxirmP440tH+Md5z41rXxqcdSYwgRBnFqkQ4jbOecPRstkjL2YhnMQxLjABwdx9IknoLvlFmgKR7aNZAUFKJgbPN57rX4Onn7rExwfGMLSmeUoLUru6xQqwtMXluOq28/Kim0rCYIYe9LRNR05hIYM57wtDecgiLHn+HH0rfsBun7+Cxy854eItf3raTPK0PyvBtx07nw89fVqlCUZwvDMi+agpMx/7KIVU/GV751NIkycckiSFBYOkfCTDiHuZYw9zhhLfM8+gsgiNG++ieH33wcAHPvLX3Di5VeC8n0RhPmCqqn4yVeXoyAv+a/SlGmTcOXtZ+Hzpvm43LIcBYXp3+OYiMyTT25CW5szdsExxmAwhIUKHAuamprAGFNjIgdisVjAGENTU1jgu5QQRREWiwUbN26EIAhh6Xa7HRaLZdRIV5mwezxJuWuacx5pNy2CyDl8RiMKDx3C4F//hulWK0ovuVjNe/mjI7A9tx3/W3ceppamtqOVp38IRSXBX72Zi6Zg5qIpKdVLJM5DDz+IKVot3lo9vrsziaIInU6n7pG8adOmcYkWpAhhZ2fwYhen04mWlhYIghAUCzhV3G43Vq9ejba2trD2mUwmOBwOCIIAQRBgMBjC7MqU3eNN5kNoEES2wBiKv/1tzP/tb1B58zfUNbvPtO/F2t90YMfhE/jmbzrQPzic9Ck+dh3B/9zzBg58HF+cW2LsaGtz4ic/2YCtW93YunX8ukxFUYTBYAgKTjBeIfscDgfMZjNEcSRkvCRJaG1thSAIad8uUgnVGNo+xftXBFYJa6jsG51pu8cb2qKHmLAMvvwyCr7wBbDASVl5eZh8zjlB5XSTi9Tx4qPHPeg67klqjfD2Nw7ipf/ZDnDgucffhbmhGtoZiddDpAfnS05s+MkDuPba62B70IY//C548YfT6VR/+FtbW9Ha2gqtVhsx3el0YsOGDbBYLLDZbLBYLGhoaIDdbkdPT4/q+W3YsEENcG+321UBUeL4Kl6d3W6HKIrQarVwOBywWCwwGo1wOp2wWq2wWCxwuVxoaWnBpk2bYDab4263IAhBY7VWqxVWqxVVVVWw2cJXm0qShA0bNqjvOedhG8tUVlaGeaTKdZIkCRaLRbW9vr4e7e3tQd3UgD9og9vtjiqqidqdS5AQExMOPjyMgccfx+CWZzF0ySUouedusFHik5o+NwP3XnMGWtr34le3nIPpZcnFaJ5zmhYlpQXoP+5FSWkhmIZ2ycoUoiii2rASALD+/g04/XOnqcHoAWDrVjeam5vR2toKwB+HtqOjAzqdLmK62WxGTU0NJEmCy+UPUet2u/HMM8+gtbUVa9asQUVFBQwGA9asWQPAHzkoVIyU4zZs2KDWYzQaUVVVhc7OThiNRvUhoLm5GVVVVdiwYUNcQqx44sCIR2q322EwGFSBiySCWq02SOjijUfsdrshCAKam5vV9waDQW1DYOhC5Tzd3d1pszuXmFBCzBgTANg45zWZtoXIDHxgACd/sh5Dr78OAPC+/DLyzlqOoi9/edTjvn7Bk3b1YgAAIABJREFUQnzlrGmYUpqcCANA+dQSXHXHCrz9ZxGmW85AcWlys6yzmSfeexwbP2hOS11fnH0hfnHhfwWlNb//BDZ+4J+wU3+GBbcuvy2pujc95a+jw9UOABAWCXjo4Qfx2KP+0JYt9hasXLly5LyymFit1ojpCkajURWYZ57xe9jKBKN4vdZnnnkmSKCVMVSn06nGHw7s0g21IRpK7F5RFFVP1eFwoLm5WfW4002g2Or1erUdod4t4Pe8q6qqUrY7sHtbp9ONW7d/KkwoIYY/OIU2ZinilMTX24uT9/wQwzt2qGn5F34RhZdfrr4/8JkHv3zzPfz4S2egMD/4qb8wwZnRvmEfNCHHzFhYjmu+9fkkrCfShSRJqNRVYt33f6CmVRtW4oYba7H+JxtQNrkMn0mf4fjx8C3yJUkadXZvaFlBEFTxrK+vh8/ni+vYwLFQAGHeYzK0t7cHdR/X1dVh06ZNAICWlhY0NjZGPC7Zrmm9Xh90nIJOp4NOpwubJd7T04Pq6vDVsInYrYh2oO1OpzPrPeYJI8SMMTOAFvjDNhITjOG9+3CysRG+gwfVtMKaGhRb6tVu6V3dg/jJqzvR2z8Ez9AwHq5ZkXSQhZ6DfXj+ifdwydeWYvaSirS0gUgPd9/TGCTCAHDdtf4NAJ98ahO++63vYfVqI264sRY1NTUwGo1wu93o6emByWRCTU1NWHqkH3qTyYS6ujrU1tZCr9dDFEW8+OKLuPXWWwEgqqBbLBY0NTUFdZVLkqR2aY+GMu4cKtySJEGn0wGA+re2thZarVb1MqOJVbJd00ajUR33VbxS5VpptVpYrVZIkgStVgu3243q6uow7zURu91ud1ivg1arVY/LZrJWiGXhdHLOpYA0AYAVgAOACf74xzEfTxljRgCZXyxIZIShd9/FyR/9CPzYcX+CRoPiO+9A0Ve+ElRu2+EB9PYPAQD++u5B1F8o4PSZiS+P37+zF88/8R48J4fw3OPv4boGw4TZqvLW5bcl3V0cD5Yzb8Wty29P+vjNWzbjyac2QX+2Pqj7d+tWN7RaLe6+pxErlq/Adddeh4aGBtTU+J/bGxsbVa8sUrri3TU3N8Nms0Gr1cJsNqO9vR2rV6+GTqeD2WxWPcT6+nrU1dWhsbERgiCgo6MDwMi4scPhUCdwKbODlYliStdsdXU1mpubIYqiOh6rzFIO7K6WJEm1VxF3s9kMs9kMt9sNq9UKAOjo6Eh7N25oO9ra2tSHBCXPYDDA5XKFdbEnanfgLHSTyaROftPr9UEPA9kIi7V70HjDGNPCH7HJBqCKcy4G5HUCMHHORcaYHkAr57xKzquPUJ3IOXcyxlox0iVdDWAD5zwtq78ZYwYAHYD/A6FMKkiUAwcOoKurC9OmTYuY7/F4UFSU2vrVROBiH/jhAbAZyY+JBuL1elFQMPqYqHI+JqRPtAZffhn9D9gArxxSsLgYk354DwouuCCs7JEjR/DUuwN4dc8xbPq3aqxcGPwk3d/fj5KS0feS/uitQ3jpt9vhG/Z/r/KL8nClZTnmfS7zT+Xx2H/kk2PoPdiHilmj34ND0l7k52uwZPHSdJoYE859YGxsV10Oe4eRV5CHgqL0b6wSrzeZLJIkoaWlRe0OHwvGug3JEii2iqcdmg7EZ//27dsBAMuWLQtKd7lcgd3n1ZxzVzpszzqPWPZwmxhjQfPRZQ8ZijBzzt2MMR1jzMg5d3LOw7dcGamzRq5DC7945+4WLERccM4x+Ic/YGDTSKhsVqHFpPXrkb80sngwxvCDS+ahYXIF5ukSW1bEOYfr+U/w1p9HxvYmlRfi6jtXYNr8suQaQRAJ0tHRMaYinM3o9Xp1jDiwaz7QU85Wsk6IR2ElADEkrQeAHtTtTATAh4cx8F+/xOBf/qKmaebPx+QN66GZNQsA4BnywbmzF1cu0wWNA+fnMcxOUISHh334v999hA9fHxl/rpg1GVffeRbKKxOPxkQQyZLtk5LGGrPZrM6aVsaGc+Ga5JIQCwBCx4MlAJXxViB726Z4yjLG7gNwb7x1A8DAwAD6+/sTOUTF4/HA6/XC4/FEzR9PmNcLNjQMrnTppoh3KHY9yvl4Cm3l/f3wbngAvvb2kXqXL0fBv/8Q3rIywOPBZ/1DuOeFffjgcD96+zyoOcv/EVKuf7R7ODAwEJY2ODCEl3+zE/t3jMQjnrW4HKtuWYqCEiT9eRgLItkfisfjweAon0MF7vN3vXMe3yzgdOHjHBqM7Tl9nEPDffD50r/OO95Z09lMtrdh1apVQe9D7Y3Hfs45OOdh3994vkPJkEtCLMLv/QaiBRB5c9IU4ZzfB+C+WOUCx4iLi4tjjsFFo6ioCAUFBaOOA4/rGHHBEHj+MFiMcd1EiDlGLJ+PJdlO39Gj6Lv7Hvg+/njknKtXo+QH64J2z/rD21344LD/C/b4Pw/j3IVaLJk2Sb3+o93DwLwTvQN4/pHt6N5/Qk1beu5MXPKvpyMvP/vG0ADE/HwWFXlRWDAU87PG+v0iNdbjtaFoMPZjxBrGwZhmzMZBs3F8NVFyvQ2x7GeMgTEW9n0pLk7PnJlQckmI2wGErojXQRZBYmIzvHsP+hobwY8cUdOKbroJRbfcHLYEyXL+bHx46CTeO9iHb31xDpZMS3ybyaOfnsBfH9mGPmnEc6y+aiHOuXpR0kueCIKYmOSMEHPO7YwxG2NMyzmX5FnTHZxzCnA5wRlyb0XfvfcCfX3+BI0GJd/9Dgqvuipi+aJ8DWxXC3j/UB++kETEo70fduPvG9+Hd2BYPh3DxV9bimUXzE66DQRBTFyyTogDli8BgJkxtjFgrbAJgI0x5gJgAGDJhI1EduHr6R4R4ZISTLrvXhTI2xByzrG7ZwBCyKSpKSX5SYnwh68fwD+e/gg+eYy0oDgPV1iWY96yzC9PIggiN8k6IVaWL8mv0DwRJL5ECIVGI3wHD2Hwz3/G5A3rkbd4MQBg2Mfxi1f349l3u/DQl6pw7oLEN+dQ4JzD9dxebHPsV9NKK4pw9Z0rUDmnNOU2EAQxccntEXeCkCn62k0ofXKTKsIA8MQ/D8C+rQvDHLjnud3Y3Z3cDObhIR+cv/4wSIQr55biuoZqEmGCIFKGhJjIKXhfH/p//gvwEyeC0hlj0EwJ7mpes2I6ZsgRjs5dUI7ZU5Kbjf3PP3Zi51uH1ffzz9Dh2nV6lFaM3yx2giBOXbKua5ogouHr6kJf493wiSKG9+7FZNsDoy6vmlZagIe+XAXnzl7UnTcLmiRnMxsuX4A9247is65+fO5fZuPCG05DXoKRmAiCIKJBvyZEzjD8wYfwyeHhht95B0NvvRWU7xkKX6hfVVkCy/mzkxZhACgpLcTV/28FzvnyAlx801ISYWJURFGExWKB3W6HxWIZNWxiPGXtdnvcoRebmprAGFNjIAdisVjAGENTU/p3+I1k42hti9XuTLUjU5BHTOQMBRdfhKID34TnV79GyQ/WoeBf/kXNc+7sxS9f3Y9HrluMedrUFt0fO9qP8qnBs6y10yfhzItn0xphIiYmkwkOhwOCIEAQBBgMBnR2Rt53aLSykiRh48aNsFqt6OzsjCsmsRJRKvR8TqcTLS0tEAQhLG5wKoxm42hti3WNxrsdmYYe7YmcouiGG1D65CYUXnqpmrZ5Wxd+9Pc96Orz4nt/6kTvyeS35RQ7evH0j97Ex64jsQsTRAhKOERFSJTwfMr+x4mU1Wq1CYuNw+GA2WyGKI5sy6+EURQEIe37LkezcbS2xXONxrsdmYY8YiIr4ZxjcPMWFF52KVjZSPQixhjyFiwIKnv6jEkoymfwDHHkaxgGhpIL7XnowwF8/A//emTnrz5EaUURZgqJrzUmsp+tW91oe6kNixYJaLW3YNGiRdjwkwfUfEmS8NBPR+/6rKysDBOh9vb2oDjHgD/4gNvtDhOPRMomgiAIcLtH9jmyWq2wWq2oqqqCzWYLKy9JkhonORqR2joao7Wtu7s7rnYn2o5choSYyDr44CD6H3wI3rY2eF9/3T8pK2Cv6FDOmDkZ9122EK3burD+ykUoL07uY10pFOLQ+16c6B5E5ZzJYd3TxKlBW5sTDz38IJ5/7gUAwG23W/D7p/8QVEar1Sb1Yy+KYlgXslarRXd3d0pl4z23Eg9d8TrtdjsMBoMqaJEEPtm2xrIlWttitTvZduQyJMREVsH7TuCk1YbhbdsAAMPbtmHwj39E0Zo1ox53UZUWFwpTUhrDLSjW4Pwb5mKf6yQuvH7pmASGP9XpfvxR9D7xeFrqmnThRZj9y0eD0noefwy9zU8AACpuvQ2Vt92RUJ2SJOGGm67Hzh3+wCCiKEKSJKxenZ4f9lAvTjlnVVVVSmXjQYnFK4qi2i6Hw4Hm5mZYLJZxFa9YbRstL9F2BHZp63Q66PWhsYGyHxJiImvw9R5B/yM2+A7uU9MKr7kahdddp77f3d2Px14/gHsvW4jSEKFMVIS9nuEwsS3VFWL11xcmbjyRE9jtrTDoDapHdvcPG7F61eqwcsl2Ta9cuVL14hR6enpQXV0ddnwiZeOhvb09yJ66ujps2rQJANDS0oLGxsaIx41F1/RobdPpdKO2O5F2KKId2Ban05lzHjMJMZEVDO/ZhZNP/Af4iZG4vsV1a1F4/fWqwLr2HUfj33bjxOAw7nluNx7+UhXy85LzgHs+7cebz3yKsy6fgblnJL/1JZFbuLe6sWiRf3yyrc2JzyQJ110bGtQt+e5as9kMq9UKSZKg1WrhdrtRXV2temlutxs6nQ6CIMQsO2o73G4IgqA+UEiSBJ3Ov9+58re2thZarVb1KqOJ01h0TY/WNr1eHzUvkXa43e4gEVbaohyXS5AQExnH+87b6N/4IDAohxQsKECJtQGFIQG+j5zw4sSgP+LR+4f60Nndj6XTEw9huH/7cXQ8ewC+IQ7XHw+iuCwfU+cnXg8RTuVtdyTcXZwIuttuR+XtdyZ9fN3aOtxw4/W4/sZarL9/A1xuFx59JD1d6QoOhwNWqxUGgwEulwvNzc1q3gMPPIBzzjlH9fhGK6ssDQL8nl99fb0qvIoH29raCkmSUFNTA8Df1a6IvNlshtvthtVqBQB0dHSkvdt2NBtHa1ukvETb0dPTo9ZnMplgMpnQ0NAAvV4Pt9udU13UjPPkZpgSfhhjBsgxkTs6OtRJBoly4MABdHV1Ydq0aRHzPR5PzGDt6YSLfeCHB8BmpCcQttfrRUGEXbAGX/obBn6/CeDyZhyTSjF5/X8i/6yzItbzq7cP4dn3uvDgNVVJifDHb/bgvRdHliYVTsrDebVzUTmvRL3+s2dHDmfY398fFig8l4jH/iOfHEPvwT5UzJo8arlD0l7k52uwZPHSdJoYE859YCw9qy5FUcSVV12OHdt3BqUPe4eRV5A3JnMEfD5fzKD08RLaLTtepLMNqRAotop3HZoeiXjs3759OwBg2bJlQekulytw6KCac+5KzvpgMn81iQkJ9/kw0PLfGPhdsyrCrGI6Jjc2RRVhAPjGyhn47Y3LEhZh7uN49++Hg0R4sq4AF92yAJXzcldcieTZ8uxmrIowPkzkBnq9Xh1rDpyFHegp5wrUNU2MO3zQg/6nfoYh1xtqmmbRaSi5fh3yZs0EAPR5hmF/twtfM8xAnmZkHJgxBm1JYh/bIa8PHVsO4OBHI4EidHNLcN71c1A0ib4CE5G2NieaHrRh1arVQd5UrpCLE5LGArPZrM6aVsaGc/G60K8QMa74jh9D/yP3Y7hzh5qW//lzUVK3DpD8wyRdJwbx/T934uOjA+jtH8J3Lpyb9Pk8fUP45+8/Re+BATVt9rIyVH9lFvIKqENoorJ6tRFHDh3NtBlJk4tiM1acCteChJgYN3yHD+DkL34M35GDalrh6mtQVHsLmCYPHH6x/OuHPfj4qP//lne6YFxSgTNjjFlG4vhRD9743ac4KY1sebn4fB3ONE6jPaMJgsgaSIiJcWFY/AgDj28AP3Hcn8AYitbcgiLTl8PKfn3lDOzq6seruyVYV81PSoSPfnISbz7zKbwD8iQwBqy4fAaElRWpNIMgCCLtkBATY47X/U8MbnoY8A76EwoKUbL2eygwXBCxvIYx3HvZAnx0ZBrOml2a8Pk+ff8YXH86CN+wv6s7r4Bh5bWzMWtpWYwjCYIgxp8JMUjGGKtnjHXKr17G2PjP+Z/AsJJJgE+eGV1ajknr7ldFmHMO99GTYccU5WsSFmHOOXa+3o32LQdUES6anIcvfn0+iTBBEFnLhBBiAB2c8yrOeRWAjQDCY5IRY0b+shUo+Npt0MyY41+eVHU6AMDr47h/22H8vzc/xd/2fRajltHx+Tjeee4wPmjrUtPKphbiom8uQMVsWp5EEET2MiG6pjnngTuMaznnUsaMmaDkn3cxis+7CKxgJIrSr3Z24++fHgMA2N49jEWnz8DnkthAxOsZRvvmAzj8cZ+aNnVBCc5dMxeFJRS4gSCI7CZrPWLGmJkxpg1JExhjzXJec2h+HHUKADrTaigRN4EiDAA3VVVgcZk/7cq55VgyKXqow9HY+Vp3kAjPW16OC26aRyJMEEROkHUesSyu9QBsAKoABHqvDgAmzrnIGBMBuOQyYIzVR6hO5JwHdkObAdgjlCMywOSCPDx0zhy8fPAEahZpgSOepOo5/cKp6NpzEr37B7D0i5VYdvFUWp5EEETOkHUeMedc4pyHxR9TJlhxzkX5rxuAjjFmlN9vjPAKHQteqRxPjD+SHLAhkGklBVgjVKQknHkFGpx//VxUf3UWPncJrREmsgO73Q5Jim8UrKmpCYwxNYBCIBaLBYwxNDWNHpYxEZxOJ6qqqsAYQ01NTZCdoijCYrHAbrfDYrGE5d16660R8zLRjlOFrPOIR2ElgFAR7QGgRxyTr+Ru6bhFmDF2H4B7E7APAwMD6O/vT+QQFY/HA6/XC48nslcYLX2sYF4v2NAwuNcbu3Ac/PGTXjTv+gwPVc/EGdrI48DK+XiMtkoHBjBlVlGw4OYD008rTuk6Kdc/2j0cGBiImJ4rxGO/x+PB4CifQwXu889K50qwjnHCxzk0GNtz+jiHhvvg8yX3QCdJEjZt2oS77roLu3btQnn5SJhNny+y7QsXLgQAfPzxx0FlnE4nWlpaIAgC1q1bF/X4RBBFES+++CJ27doFURRx2WWXYe3atWhpaQHgj2T0wgsvQBAELFy4EAaDAbt27VLznn/+eSxevDgsb7zbkSzxnJtzDs552G/BWP0G5JIQCwjupob8vjKeg2VP2BrvyTjn9wG4L1a5wOhLxcXFSUfnKSoqQkFBwagRlsY1+lLBEHj+MFiEiEmJ0iL24hc7/Lfunq1H0PyFeZg7OXw8WDkfi9JOzjl2/F83dvzjKJZfOh2Lz0tv3FHl+o92D3M5+hIQ2/6iIi8KC4ZiftZYv1+k0hUJKV40SF/0pajnYByMaZKOMKTT6WC1WnHXXXdBowmvJ1K9bW1tMJvN2L17t5ovSRI2b94MQRBQXV2dtohHe/bsUb3SxYsXw2q1orm5GRqNRg2isHjxYgBAdXU1enp68NJLL6ne7+LFi6HRaILylG0mx7MdqRDLBsYYGGNh35fi4vREowsll4RYhN/7DUQLmnyV9VwwYzJ+vasbn3l9mFGSj5K85L6Ie9wSdvzDvz/wey8eQWllIWYuSXzDD4LYutWNtpfasGiRgFZ7CxYtWoQNP3lAzZckCQ/9dPQu1MrKSjWucDoQBAFu98gCD6vVCqvViqqqKthstojHSJKkxiaO187QvZk7OzvVtPb2dgiCEJSv0+ngdrvR3d0dNS+wzmTaMdHJJSFuh3+yVSA6yN4okb3MnVyI9foZ2Lz3OBpXzMSk/OSEeP6KKdi77Rh6Pu3HdGESKufntndKZIa2NiceevhBPP/cCwCA22634PdP/yGojFarHTfREEVRjWOueKR2ux0Gg0EVtGiBDdJhp9vtRmtrq2pLaCQqrVaL7u7uUfNSbcdEJ2eEmHNuZ4zZGGNazrnEGNPDv1GHO+bBRMY5Q1uMz09LbXervHwNzrt+Djrf7MXpF02FJo8mZWUbHc/vgfvvn6Slrvln6HB5/fKgNNfzn8D9wl4AgP7yBai+YmFCdUqShBtuuh47d3wMwC8ekiRh9erMCYTdbofZbIYoiqo9DocDzc3NsFgsYypeTU1NsNlsqsCGerOA/5pVVVUBwKh5ibZDCV8I+D1rvT60w3PikHVCHLB8CQDMjLGNARtwmADYGGMuAAYAlkzYSIwPJz/zoqQ8P2hSVtGkfHxu1bQMWkXkMnZ7Kwx6gyo8d/+wEatXrQ4rN55d0+3t7UH11NXVYdOmTQCAlpYWNDY2Rj02ma5pBafTCb1eHySAK1euVL1ZhZ6eHlRXV0On00XNS7QdimgHtmMix1jOOiGWRbdJfoXmiSDxnRAc2nUC7ZsP4PSLKrHk/Ljm4xFETNxb3Vi0yD/O2dbmxGeShOuuDd96fqy6pt1uNxYvXqw+CEiSpAa0V/7W1tZCq9WqXuVo4pSsnYo3qtQtSRI6OjpgNpthtVohSRK0Wi3cbjeqq6tVwVbylLFhJS+Rdrjd7iARVtqhHDcRyTohJojdLgnbnjsEzoH3HV2YpC3AnGXlsQ8kMk71FQsT7i5OBMMVC1B95aKkj69bW4cbbrwe199Yi/X3b4DL7cKjjzyeRgv9oqaso7Xb7aivr1eF94EHHgBjDK2trZAkCTU1NQD8XeSCIMBsNsNsNsPtdsNq9S/y6OjoSGu3rdPphMlkCkvv7e0FADgcDlitVhgMBrhcLjQ3N6tlHA4H7rrrLlRXV6t5ibajp6dHrc9kMsFkMqGhoQF6vR5ut3tCdlEzznmmbchpApcvdXR0qJMVEuXAgQPo6urCtGmRu109Hs/4Ll8S+8APD4AlsfdzJLxeLwpiLIXyHerH9oMe7PrguJo2aUo+zr9xHsqnjX3bles/e/bsiPn9/f05vXwpHvuPfHIMvQf7UBEjBvQhaS/y8zVYsnhpOk2MCefpW74kiiKuvOpy7Ni+Myh92DuMvII8FBSlf4tUn8+HLVu2hHmEuYTP50tpCVKg2Cqed2j6WBKP/du3bwcALFu2LCjd5XKpXfEAqjnnrnTYlPkFXQQBYHiYw/Xe8SAR1s4qxkXfXDguIkxMPLY8uxmrIowPE2OLXq9Xx5oDZ2EHesoTDeqaJjLOoGcYb/2jG90Be03PXDIZK6+bg/xCelYk0k9bmxNND9qwatXqIK9srJnIE5ICMZvN6ji1MjY8ka8LCTGRUfqOD+GfL3fhxLEhNW1RtRZnXT4DGg0tTyLGhtWrjThy6Oi4n9doNGbFzlLZwEQW3lBIiImM0XvUgzdfOQrPwMjer2ecPQVLrphBgRsIgpgwkBATGeHgvn50vNaN4WH/ZEGNBtCfUYa5Z5STCBMEMaEgISbGnc4dx/Fex0j8joJCDc67eCp0NIOfIIgJCAkxMW5wzvG+S0LnjhNq2qTSPFywahpKywvAD+d2mMGJCOecejCIU4pMLOklISbGheEhDvcb3Ti4byS+Z8XUQpx38VQUFad/vSYxtmjklY/Dw8PIz6efEeLUQYlXPJ4PmDR9jxhzPAPDePOVniARnjWvBP9inEYinKPk5eXD5+M4dPgghoaGMuJFEES64Zyr65kLC8Njpo8V9ChLjDl7dvVB6vaq76tOL8WZei0YLU/KWaZMqkTviSPo6+vD7j2nVkhwzgHGMCafz1OhKz/X2xCv/dOnTx8Ha/yQEBNjzmlnlKHn6AAO7/dgebUWVaenFg6RyDz5eQWoKJ2Oz052Y3hoCD74Yh+UBriPj/kD3PCgD3mFGhQUpr+3JtdFDMj9NsSyPy8vD1qtFqWlpeNmEwkxMeYwDcPZ52lxXPJh2qz07F1NZJ78vAJUls0c13OOx57rvQf7UKGdjOkL0h9oJNf3Kwdyvw3ZaD+NERPjQl4+IxEmCIKIAAkxQRAEQWQQEmKCIAiCyCA0Rpw6an+rEsMyGbq6utDb24uKioqI+YODg+M6nZ7v7wfvHgST0nPOoaGhmOtNlfOxE5kZv1Gu/8GDByPmDwwMoLg4d7vX47G/5+AJHDs6gPKj2dnOwUEvCgtHj2udKkr7dUfTP1kn1z9DQO63IRX7Q37j03YRGK3/Sw3G2E0A/jfTdhAEQRDjytc450+noyLqmiYIgiCIDEIecYowxqYBuFR+uweAsmFyB4DqkOKhaYHvY/0fqb54iOe4aGVitSGe9qRq/2j2xcpPxz0Yrb3xkol7kM7P0Gj2xcqnexA9LVHoHkRPG8978DX5/xc5511J1hMECfEYwRjjnHM2Wlrg+1j/R6ovWTviLROrDfG0J1X742lDsvaP1p5IdufSPUjnZyieNtA9oHuQqv2ptCFb7kGyUNc0QRAEQWQQEmKCIAiCyCAkxARBEASRQUiIx44fx5H24wT/T5cd8ZaJ1YZ42pOq/fHUkaz9oe8zZf9oZZK9B+n8DMVTB92D8Pd0D2LbEW+ZXLkHSUGTtXKEsZokMF7kuv1A7rch1+0Hcr8NuW4/kPttyEb7ySMmCIIgiAxCQkwQBEEQGYSEmCAIgiAyCAkxQRAEQWQQEuLcYUxm640juW4/kPttyHX7gdxvQ67bD+R+G7LOfpo1TRAEQRAZhDxigiAIgsggJMQEQRAEkUFIiAmCIAgig5AQn0IwxvSMMTNjTJtpW1LlVGgDQSQLY0zItA0TjUxecxLiUwTGmBmAkXNu55xLmbYnGRhjDYyxTsZYZ6ZtSRbGmMAYa820HcnCGLMxxhyMsYZM25IMjLFWxlgvY8yWaVuShTGmB+DKtB3xIn9mzLn6mQEyf81JiE8B5Ce5Rs55U6ZtSZFKznmV/MrJhwkARgA56c0zxgTOuZVzbgJgybQ9icIYM3LOazjnFQDqc9Wr5Jy7AfRk2o54kB2Abs65HUCl/D5gp7VOAAAH+ElEQVTnyPQ1JyE+NTAC6JE9SkcuduvKP5p6xhjP1S+zbHdLpu1IFs65CKjeQc55lJxzZ8DbDuSImOU4tQBE+f92AKYM2pKzkBBniEhjuXK3ZrOc15yAoBoANMsesRtAY7rtjUQ628A5F2VPzABg03g8TKTTfsaYEYAzZsE0k+bPkdq7AsCablujnC+t9svHawG4x6tXZSzakEkSbI8WI0IsAciKXohcuyf5mTZgoiHf/Hr4PY4q+D+8Cg4AJs65yBgT4R+zqJKPq49QnSh7ARIAXUAdNWNkPmRbxqINAPxdRIyxFvi/0O4cst+CEfGqZow1jOVQwVjdA9krrpF/qPRyl13O2C+zhnM+5g8SY9yGcSfJ9iji60awKGeEZO9JxuGc0ysDLwAcgBDw3gygM6RML/wTsGLVZYTfIwb8H8L6XGtDyDHNuWo//D9Gjlz8HIUcYwOgzTX75e+CVv5fSIeNGfocdabDtrFuj5zXIKc1ADBnyu5U70kmrzl1TWcPKxH+NNkDQB/rQC57xfIYpZZzvnEM7IuHpNsgj2+3ym1oHgvj4iBp+7OIVO6BOeAeOHhmJsylZD+AVgAueeZ9pu5bSp8jeYxekP9mA1Hbw4MnaVXK77ORUe9Jpq85dU1nDwKCu1Egv6+M52A+Dl1xcZB0G3h2zPhO6R4AgCxemZywkso9sAPI9A9prtsPpP5ddgNg6TYqBUZtT8BvTzZc+2jEakNGrzl5xNmDiPBlL1oAubSmNtfbkOv2A7nfhly3Hzg12hDIqdCerG4DCXH20I7wGYc6+Jdh5Aq53oZctx/I/Tbkuv3AqdGGQE6F9mR1G0iIswRlbEWZUi+PVXTwMZq1Ohbkehty3X4g99uQ6/YDp0YbAjkV2pPtbaAx4nEmYHo9AJgZYxsDJsWYANgYYy7419Nm5e5Gud6GXLcfyP025Lr9wKnRhkBOhfbkahuYPG2bIAiCIIgMQF3TBEEQBJFBSIgJgiAIIoOQEBMEQRBEBiEhJgiCIIgMQkJMEARBEBmEhJggCIIgMggJMUEQBEFkEBJigiAIgsggJMQEQRAEkUFIiAmCIAgig5AQEwRBEEQGISEmCIIgiAxCQkwQo8AYMzLGXIwxzhizyS8HY8wWUs7FGDNnys5sJhuvjRwGL1J6PWPMGJJmZIx1yp+B+oB0Qf4sdMpl9Iyx0Ji3BBETir5EEDFgjDUAsHHOmfxeC2A3gA2c8yY5TZ8tsU0zjSxGPUr4uWy7NvL9DAyPF5jXCUDinBtC0psBGDnnVSHpegC1nHOr/N4MQMym9hLZD3nEBJEg8g94B4CVAWn0wwtVhF0AdEpaNl0bWSilKCJsBGAFoI/mMUdAAtCtvJED0Demw1Zi4kBCTBAJIotNNQCH/F4vd1E2yO+V7ux6xlgzY6w3sGtW7t42Kl3dcpo54JhOxliDnMYD6tXKZcJEQq7LzBhrVbpW5ff1clroeRrk9F6l/kRsi1YWgABAC39Qdn3otQmoryHANsXeUa9bHPclkj2hbALQEiXPJAtpqmIqBraXIGLCOacXveg1ygtAAwAOoBl+8eUA6kPKOAA0BLzvVcrIx7vk/7UAOgPK8cD/5bJaAFo5rRVAc0AZWwT79AAcAfUb5bTWgDQOwBxwHnPAsVz+G7dtcZQVIl0b+TyugDwhsHy06xbHPYpqT+i5ohwvBFwTIUIbmgPrDzmuISTNrNwPetErnhd5xAQRJ5xzC+fcBKAKgJUx5ohxiCj/dcMvFOCcS5zzKnmijxlQx5wVnHIZpet0A4B62Rs2Angmgl1uADrGmAt+8XACqJXrrgewBn4vL8w2+Vg3/OOfcdsWR9lo1AZcF3DORfl94ASpsOsWizjtEeDvSo6EBcBK2ZO2yDZYA/KjHRcJCf4eE4KICxJigkgQWTyaARgTnSUrC2qrXE+oOEY6lxuAE/6uUhOPMt7K/ZOLnABcsvhq4Z80tFF+1Yxyvh74u1Pjti3RdoQQes0SEblU7OmJdiyAbs65VXnBL8L1AWLeHcFuwO/hixHSI56LICJBQkwQybES/kk/kX6ER2MNAD3nXIzTgwQAG4B6BEwKCkQeFxVkATEBqIG/O7heGU+WPcX6SMfDLzDOBG2LVTba8c3wT4YSQspGG7dVkceWo9Ub03b5ISZSng0hPQYBYq5cM7tiQ8ixJvivXSAC/N48QcRHpvvG6UWvbH7B32XqwsgYsTJO7ID/hx/we0W9cpogH8Ph/4HXwj/Oq4zDCgA65bJmue5W+X/lHNoIdnRGSg+wsRN+0bAF2GWT7epEwNhywHkaQsrHbVu0snJes/zeHHptAuxtle1tCDj/aNdNGedujnINotoTUk5ZhqS8V9oWOuav2M2V8gGfhVbZxmYEjCMHHNuqtIle9IrnReuICSIHYIw1c84taaqLAzDwLFpWFA+yp7uGc74xxXrSdi0j1C0AsHB5XTFBxAMJMUFkKbLwKB62nSfeDR6t3lwVYiP3T0RLtR4Bfo810bHtWPVq4fee01ovcepDY8QEkb2sgb8rFGkUYWWM05LAGHVWkA4RlusRx0osSYSJZCCPmCAIgiAyCHnEBEEQBJFBSIgJgiAIIoOQEBMEQRBEBiEhJgiCIIgMQkJMEARBEBmEhJggCIIgMggJMUEQBEFkEBJigiAIgsggJMQEQRAEkUFIiAmCIAgig5AQEwRBEEQGISEmCIIgiAxCQkwQBEEQGeT/A/loIPwHi5E0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 244.354x151.019 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 152,
       "width": 241
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax1 = plt.subplots()\n",
    "plt.title(\"BHB merger timescales\")\n",
    "#ax1 = fig.add_subplot(111)\n",
    "#ax2 = ax1.twiny()\n",
    "\n",
    "# Data\n",
    "ax1.loglog(xaxis_t, t_gws60, label = \"GWs $q=1, M=60~M_\\odot$\", linestyle = \":\")\n",
    "ax1.loglog(xaxis_t, t_bteau60, label = \"Baruteau, $M=60~M_\\odot$\")\n",
    "ax1.loglog(xaxis_t, t_acc60, label = \"Accretion, $M=60~M_\\odot$\", linestyle=\"-.\")\n",
    "\n",
    "# Data\n",
    "#ax1.loglog(xaxis_t, t_gws100, label = \"GWs $q=1, M=60~M_\\odot$\", linestyle = \":\")\n",
    "#ax1.loglog(xaxis_t, t_bteau100, label = \"Baruteau, $M=60~M_\\odot$\")\n",
    "#ax1.loglog(xaxis_t, t_acc100, label = \"Accretion, $M=60~M_\\odot$\", linestyle=\"-.\")\n",
    "ax1.loglog(xaxis_t, t_all100, label = \"$q=0.01, M=100M_\\odot$\", linestyle=\"-.\")\n",
    "ax1.loglog(xaxis_t, t_all200, label = \"$q=1, M=200M_\\odot$\", linestyle=\"-.\")\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# GW Detectors\n",
    "ax1.axvspan(lisa_min, lisa_max, alpha=0.2,  color = 'C4' ) #LISA band\n",
    "ax1.axvspan(et_min, et_max, alpha=0.2,  color = 'C7') #ET band\n",
    "ax1.axvspan(ligo_min, ligo_max, alpha=0.3,  color = 'C6') #LIGO band\n",
    "\n",
    "# Fine tuning\n",
    "ax1.set_axisbelow(True)\n",
    "\n",
    "# Turn on the minor TICKS, which are required for the minor GRID\n",
    "ax1.minorticks_on()\n",
    "# Customize the major grid\n",
    "#ax1.grid(which='major', linestyle='-', linewidth='0.5', color='red')\n",
    "# Customize the minor grid\n",
    "#ax1.grid(which='minor', linestyle=':', linewidth='0.5', color='black')\n",
    "\n",
    "ax1.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator(10))\n",
    "#ax1.set_yscale('log')\n",
    "ax1.grid(True, alpha=0.1, linewidth=0.5)\n",
    "ax1.axhline(y=.5, alpha=0.2, linestyle = \":\", color = 'grey')\n",
    "ax1.axhline(y=10, alpha=0.2, linestyle = \"-.\", color = 'grey')\n",
    "ax1.axhline(y=1e4, alpha=0.2, linestyle = \":\", color = 'grey')\n",
    "\n",
    "\n",
    "ax1.set_ylim((0.013157894736842105 / 31556926),1e9) # ymin is 1/f_isco\n",
    "ax1.set_ylabel(\"Time (yrs)\")\n",
    "ax1.set_xlabel(\"Binary separation, a (AU)\")\n",
    "plt.legend(loc=4)\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"BHB_in_AGN_disk-25oct19-vary_qx2.pdf\")"
   ]
  }
 ],
 "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.5"
  },
  "name": "gas-drag.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 1
}