notebooks/new_features_tutorial.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a target=\"_blank\" href=\"https://colab.research.google.com/github/CosmoStatGW/gwfast/blob/master/notebooks/new_features_tutorial.ipynb\"> <img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/>\n",
    "</a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TUTORIAL OF THE NEW FEATURES IMPLEMENTED AFTER <a href=https://arxiv.org/abs/2207.06910>arXiv:2207.06910</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "\n",
    "import copy\n",
    "import numpy as onp\n",
    "from astropy.cosmology import Planck18\n",
    "\n",
    "PACKAGE_PARENT = '..'\n",
    "SCRIPT_DIR = os.path.dirname(os.path.realpath(os.path.join(os.getcwd())))\n",
    "sys.path.append(SCRIPT_DIR)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rc\n",
    "rc('text', usetex=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gwfast.gwfastGlobals as glob\n",
    "import gwfast.gwfastUtils as utils\n",
    "from gwfast.waveforms import TaylorF2_RestrictedPN, IMRPhenomD_NRTidalv2, LAL_WF\n",
    "from gwfast.signal import GWSignal\n",
    "from gwfast.network import DetNet"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FEATURE EXAMPLE: computation of $\\texttt{TaylorF2\\_RestrictedPN}$ down to the Kerr ISCO of the remnant BH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select a GW170817-like event for example\n",
    "z = onp.array([0.00980])\n",
    "\n",
    "event_ex = {'Mc':onp.array([1.1859])*(1.+z), \n",
    "            'dL':Planck18.luminosity_distance(z).value/1000., \n",
    "            'theta':onp.array([onp.pi/2. + 0.4080839999999999]), \n",
    "            'phi':onp.array([3.4461599999999994]),\n",
    "            'iota':onp.array([2.545065595974997]), \n",
    "            'psi':onp.array([0.]),\n",
    "            'tcoal':0.,\n",
    "            'eta':onp.array([0.24786618323504223]), \n",
    "            'Phicoal':onp.array([0.]), \n",
    "            'chi1z':onp.array([0.005136138323169717]), \n",
    "            'chi2z':onp.array([0.003235146993487445]), \n",
    "            'Lambda1':onp.array([368.17802383555687]), \n",
    "            'Lambda2':onp.array([586.5487031450857])\n",
    "           }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
     ]
    }
   ],
   "source": [
    "TF2Schw = TaylorF2_RestrictedPN(is_tidal=True)\n",
    "TF2Kerr = TaylorF2_RestrictedPN(is_tidal=True, which_ISCO='Kerr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The cut frequency for a remnant Schwarzschild BH is 1590.1 Hz\n",
      "The cut frequency for a remnant Kerr BH is 3433.5 Hz\n"
     ]
    }
   ],
   "source": [
    "fcutSchw = TF2Schw.fcut(**event_ex)\n",
    "print('The cut frequency for a remnant Schwarzschild BH is %.1f Hz'%(fcutSchw))\n",
    "fcutKerr = TF2Kerr.fcut(**event_ex)\n",
    "print('The cut frequency for a remnant Kerr BH is %.1f Hz'%(fcutKerr))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FEATURE EXAMPLE: waveform overlap "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### As an example, we check that the gwfast and LAL implementations of $\\texttt{IMRPhenomD\\_NRTidalv2}$ are equal on 100 random events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Random sample 100 BNS events\n",
    "nevents=100\n",
    "\n",
    "zs = onp.random.uniform(1e-2, .5, nevents)\n",
    "\n",
    "dLs = Planck18.luminosity_distance(zs).value/1000\n",
    "\n",
    "Mcs = onp.random.normal(loc=1.156, scale=0.056, size=nevents) \n",
    "   \n",
    "events_rand = {'Mc': Mcs*(1.+zs), \n",
    "               'eta': onp.random.uniform(0.23, 0.25, nevents), \n",
    "               'dL': dLs, \n",
    "               'theta':onp.arccos(onp.random.uniform(-1., 1., nevents)), \n",
    "               'phi':onp.random.uniform(0., 2.*onp.pi, nevents), \n",
    "               'iota':onp.arccos(onp.random.uniform(-1., 1., nevents)), \n",
    "               'psi':onp.random.uniform(0., 2.*onp.pi, nevents), \n",
    "               'tcoal':onp.random.uniform(0., 1., nevents), \n",
    "               'Phicoal': onp.random.uniform(0., 2.*onp.pi, nevents),\n",
    "               'chi1z':onp.random.uniform(-.05, .05, nevents), \n",
    "               'chi2z':onp.random.uniform(-.05, .05, nevents), \n",
    "               'Lambda1':onp.random.uniform(0., 2000., nevents), \n",
    "               'Lambda2':onp.random.uniform(0., 2000., nevents),\n",
    "              }\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All available detectors are: ['L1', 'H1', 'Virgo', 'KAGRA', 'LIGOI', 'ETS', 'ETMR', 'CE1Id', 'CE2NM', 'CE2NSW']\n",
      "Using detectors ['L1', 'H1', 'Virgo']\n"
     ]
    }
   ],
   "source": [
    "# Initialize the LIGO-Virgo network\n",
    "\n",
    "alldetectors = copy.deepcopy(glob.detectors)\n",
    "print('All available detectors are: '+str(list(alldetectors.keys())))\n",
    "\n",
    "# select only LIGO and Virgo\n",
    "LVdetectors = {det:alldetectors[det] for det in ['L1', 'H1', 'Virgo']}\n",
    "print('Using detectors '+str(list(LVdetectors.keys())))\n",
    "\n",
    "# We use the O2 psds\n",
    "LVdetectors['L1']['psd_path'] = os.path.join(glob.detPath, 'LVC_O1O2O3', '2017-08-06_DCH_C02_L1_O2_Sensitivity_strain_asd.txt')\n",
    "LVdetectors['H1']['psd_path'] = os.path.join(glob.detPath, 'LVC_O1O2O3', '2017-06-10_DCH_C02_H1_O2_Sensitivity_strain_asd.txt')\n",
    "LVdetectors['Virgo']['psd_path'] = os.path.join(glob.detPath, 'LVC_O1O2O3', 'Hrec_hoft_V1O2Repro2A_16384Hz.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using ASD from file /Users/francesco.iacovelli/Desktop/PhD/Research/2021-09_paper_BNS_massfun_MCMC1/GWfast/psds/LVC_O1O2O3/2017-08-06_DCH_C02_L1_O2_Sensitivity_strain_asd.txt \n",
      "Initializing jax...\n",
      "Jax local device count: 1\n",
      "Jax  device count: 1\n",
      "Using ASD from file /Users/francesco.iacovelli/Desktop/PhD/Research/2021-09_paper_BNS_massfun_MCMC1/GWfast/psds/LVC_O1O2O3/2017-06-10_DCH_C02_H1_O2_Sensitivity_strain_asd.txt \n",
      "Initializing jax...\n",
      "Jax local device count: 1\n",
      "Jax  device count: 1\n",
      "Using ASD from file /Users/francesco.iacovelli/Desktop/PhD/Research/2021-09_paper_BNS_massfun_MCMC1/GWfast/psds/LVC_O1O2O3/Hrec_hoft_V1O2Repro2A_16384Hz.txt \n",
      "Initializing jax...\n",
      "Jax local device count: 1\n",
      "Jax  device count: 1\n"
     ]
    }
   ],
   "source": [
    "# Initialise the signal and network classes\n",
    "# NOTE: the initilisation waveform does not need to be one of the two waveform of which one wants to compute the overlap\n",
    "myLVSignals = {}\n",
    "\n",
    "for d in LVdetectors.keys():\n",
    "\n",
    "    myLVSignals[d] = GWSignal(TaylorF2_RestrictedPN(), \n",
    "                psd_path=LVdetectors[d]['psd_path'],\n",
    "                detector_shape = LVdetectors[d]['shape'],\n",
    "                det_lat= LVdetectors[d]['lat'],\n",
    "                det_long=LVdetectors[d]['long'],\n",
    "                det_xax=LVdetectors[d]['xax'], \n",
    "                verbose=True,\n",
    "                useEarthMotion = False,\n",
    "                fmin=10.,\n",
    "                IntTablePath=None) \n",
    "        \n",
    "myLVNet = DetNet(myLVSignals)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now compue the overlaps on the same events for the network and check they are 1\n",
    "overlap = myLVNet.WFOverlap(IMRPhenomD_NRTidalv2(), LAL_WF('IMRPhenomD_NRTidalv2', is_tidal=True), events_rand, events_rand)\n",
    "\n",
    "onp.allclose(overlap, 1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# One can do it for the single detector too\n",
    "overlap = myLVSignals['Virgo'].WFOverlap(IMRPhenomD_NRTidalv2(), LAL_WF('IMRPhenomD_NRTidalv2', is_tidal=True), events_rand, events_rand)\n",
    "\n",
    "onp.allclose(overlap, 1.)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FEATURE EXAMPLE: usage of $\\texttt{TEOBResumSPA}$ waveform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gwfast.waveforms import TEOBResumSPA_WF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use TEOBResumSPA with all modes up to l=4 for a BNS\n",
    "# For references see arXiv:2104.07533, arXiv:2012.00027, arXiv:2001.09082, arXiv:1904.09550, arXiv:1806.01772, arXiv:1506.08457, arXiv:1406.6913\n",
    "mywfTEOBHM = TEOBResumSPA_WF(is_tidal=True, modes=[[2,1], [2,2], [3,1], [3,2], [3,3], [4,1], [4,2], [4,3], [4,4]])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot only the amplitude as an example\n",
    "fgrid = onp.geomspace(5., mywfTEOBHM.fcut(**event_ex), 10000)\n",
    "AmplHM = mywfTEOBHM.Ampl(fgrid, **event_ex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAHTCAYAAAA6fiz2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAA9hAAAPYQGoP6dpAABfV0lEQVR4nO3de3xT9f0/8FeSQrmUNg0Iys32FNAhCKQt3t2UU3GyzSkJ6Kbz2lZXp9PNpp1ufpmXNvW2uVVJdaLOGyTe9hMnJOD9RtqAgqi0Pa0gXlDSNFykQJLfHyyxadI2TdKcJH09H48+NOdzLu+0H5p3P5/PeR+Fz+fzgYiIiIj6pJQ7ACIiIqJUwKSJiIiIKAJMmoiIiIgiwKSJiIiIKAJMmoiIiIgiwKSJiIiIKAJMmoiIiIgiwKSJiIiIKAJMmoiIiIgiwKSJiIiIKAJMmoiIiIgiwKSJKAUVFBRAoVAM6Ku8vDxwfG5u7oCOzc3NhcvlConDZrOhpKQk7PkKCgpgMBhCjnO5XP1ev6CgAOXl5ZAkKehYSZJ6PTY3NxclJSUwGAyD8S1POw0NDUE/u8LCwsDPy2azITc3N+j7H8nPTK/Xw+FwDCgG/7Hh+hdRslHwgb1EqcXhcKCwsBAAoNVqQ9p6267VatHU1AQA0Ov1kCSp1/27kyQJLpcLVqsVoigGtvsTJj9BEKBWq4OO8W9vbW0NOqf/w9X/odzb+wAQdF2XyxWIXZIkqNVqCIIQck21Wo1169b1+b6GssLCwsD3uPv3z08QBEiSFPS976vPuFyuoOONRiMqKyv7jaOgoCBwnMlkQllZWYzvjGiQ+YgopTQ1NfkA+ERRDNre2trqA+DTarURbe/o6PAB8KnV6j6vp9PpfAB8Vqs1aLtWqw3E0dHRETZOtVrtA+Azm80Rvw8/k8kUiK/n+f1tZWVlIe+prKys1+Poh59nuJ+b0Wj0AQh89fyZ99dnuh/f89hwul9Lp9NF/Z6IEoXTc0Qpxj+aEymNRjMocfhHHMxmc9iYtFptYOTAarUO+PxlZWUoKyuDy+VCQ0NDRMeo1erAiIXL5UJNTc2Ar5vOJEmCxWKBIAiwWq0hP7fKysqIRoh6U1lZCZPJBABB08Hh2Gw2AIBOpwt6TZTMmDQRpRhBEGA0GiNeu6NWq2E0GmE0GqO6XnV1NSorK4Om5vxTKmq1us8kzj9913NtUqT0ej0AwG63D+g4/wc2P4iD+RNdf6ISTnV1dUzXKCsrg1qtDpouDcefSJeUlEAQhJApPqJkxKSJKAX1TGLivX93Wq02JOHqvl6pL+HWywyE/7hoR8u4uHjg1Gp1oK9E+30vKioC0PfP3Z/AFRUVBa43kEXkRHJg0kREISwWS58fYE6nE0D/SZP/Q9e//0D5p3q6LzgfyHGRJIoNDQ0oLCwMugOvrq4uouvU1dWhsLAw6M6ywsLCwDm7T1H57xTzf/lH0fwMBkNQe88Y9Ho9CgsLUVJSAofDgfLycuTm5qKgoCAwomYwGFBQUIDc3Nyw78G/eLuurq7PhMZkMsFkMkW9kL6xsRFA3/3DH7NWqw38fKOZxiVKKLkXVRFRfPS24Ls3fS3qBeATBKHXY81mc8SLd3u7Rl8LwTs6OoIWLPcUyUJwAL7W1tZe42ptbQ0sZve/X0EQgl73dnxHR0fIvlqtNmhbz9itVqtPFMXANXt+T0wmU9A5en5vu8fa/brdr+U/r38BflNTU0js3b8/RqOxz+9RuPfd28+z+/vor//4f/b+vurvu30dQ5QMmDQRpYlokyb/B67/y//hnKikSa1WB12/eyJQWVkZ9pzd76zzx9z9OLVa3e/dW92Tk+53kXVP2Hr7XvoTj8rKypA70Do6OoLuLOypv8TDarX2+r31v290uyOx+x1r/iTS/73tmVT6+d9f9++XTqfzmUymPpOoeN0959+v+8/Xn+jxjkdKZhmxjlQRUeqTc8G0v5hiT/4F5i6Xq9fF5r0d232dTDj+wpk6nQ5msznkumazGSUlJbDZbLBYLCELp/3TT+EW1/sX3g90SnEgRFEMxNR9Cs0/Lenf1tsUnNlshsPhgMlkgs1mC9xVZ7FYAufv7a5I4Mj3vfv789fN8jMajX1+//0L+7ufo6ioCDabDY2NjVGvvyMabFzTRDTEqdVq+I6MOge+EkkUxZDrt7a2orq6OrBmqLcF3WVlZUHH+W+jt9lsfa5L8pcwePjhh3vdx78eKdw6m+4FH8Ot/RJFESaTKeo7FvvTs6go0PcdceFotVqYTCa0traio6MDZrM5cOebzWYLFFDtjc1mC3xJkgRBEKDT6dDU1NRv2QJ/outfMN79PXFdEyUzjjQRUYiysrJ+PzQjEe3da4IgoLKyEoIgQK/XQ6/XR/RhKooi1q1bh8LCQtTU1IT98O4+IlJaWtpv7P5Rpe6MRmPI6IxWq0VRUREKCwshiuKgVrcOt8A6lnpcarUaOp0OOp0ORqMRCxYsgMPhgMFg6HU0raOjI6pr+UsRdK8gD/wwOsYyEZTMmDQRUQj/NE9vuk+dDSadTgdBEIJGM/qj1WqhVqvhcrnCTu11T5r8CU9ferumf4pr5cqVgfi6jzr1N8UVi3DnjOQ6LpcLq1atwpIlS3rd3z89WVBQAIvFEvfRMn9SJEkSFApFSDvLDlAyY9JERAMWaSmBWOssAUeSIP+amUiSJv/1ekvo/OeIZbTE5XLB6XRCq9WGPIPNZrNh5cqVsFgsyM/PR1NTU8RxD7ZVq1ahvLwcra2tfSZD/lGgwSg26X/+Yc+RJuCHUSj/sxKJkg3XNBFRvwoKCoKm67o/mLcv3SuHRyseiVf30Qv/h7X/w7m/48K9xwULFgQ9bNbPP81lNpthMpngcrn6HbXraTBH7/znjmQKrK8F+LHwX9tqtaKpqSnoa8mSJRHHRyQHJk1E1K+eU0/+kZP+Hn3hv0uq+4LfaK8bzciD0+mExWJBYWFh0MJw/wJwvV7fa5LicDhQWFgYUoQS+CH56OtRNv73PNAkKBELoR0OR5+JSfe76OKpe38JN/rmXww+0MfmECUKkyYiioo/iSkvLw+bGHS/gy2a2++738E10IfIdr/lPtwHsH/RsyRJKCwsDEkgHA4HFixYAABYunRpr9exWCxhEyeHwxFItnomXf7n9YUb6aqrq4v44cSx6q3yucViCcQc7/VM/kX1vSVj/u0caaKkJVN9KCKKE6PRGFSUEt2KVYarCO3zHSluGG7/3r78+3XnL8Lo//JXxdZqtYFChf7tfV2/Z3HLnsd3LxDpLzzpL0zpP7bn+/QX3+x+nnBFG7tXx8b/ill2L5LZW3HInpW//cf2jL2347sXgfR/3/zHdT9erVYHvn9lZWWB62q1Wl9lZaWvqakp8PMRBCFQqNNfvNJftNL/3v3X7Vld3B9D923+4pl99Zneio+G0/Pn1rN4p79ievfzm0ymiM9PlAhMmohSXM/qzt2/wlVl7ujoCPpgjvQrXHVr/wdduPMJgtBrxez+ru9PhoxGY9Cxra2tYY/t+QHv8x1JEPxJR7h2v6amJp9Opwuc15989FXRWqfT+URRDDy2pWc18v6O9/l+eGxKz/fb0dERqIwuCELg+x7uES3+5LD7tZuamkK+R/6k0mq1+tRqta+1tdXX2trqq6ysDErYtFqtr6ysLCTB7O1nNpDq8z2PV6vVQX2jeyLZX9JJJBeFz5fgSnZEREREKYhrmoiIiIgiwKSJiIiIKAJMmoiIiIgiwKSJiIiIKAJMmoiIiIgiwKSJiIiIKAJ8YG+ceL1efPnllxgzZkzYJ3cTERFR8vH5fNizZw8mTpwIpbLvsSQmTTGqr69HfX09Dh48iNbWVrnDISIioijs2LEDkydP7nMfFreMk87OTqjVauzYsQPZ2dlyhxOWx+NBa2srCgoKoFKp5A6Hkgz7B/WF/YP6k6p9xO12Y8qUKXC5XMjJyelzX440xYl/Si47Ozupk6asrCxkZ2enVIemxGD/oL6wf1B/Ur2PRLK0hgvBiYiIiCLApClG9fX1mDlzJoqLi+UOhYiIiAYRk6YYVVRUYOvWrbDb7XKHQkRERIOISRMRERFRBJg0EREREUWASVOMuKaJiIhoaGDSFCOuaSIiIhoamDQRERERRYBJExEREVEEmDQRERERRYBJExEREVEEmDTFiHfPERERDQ1MmmLEu+eIiIiGBiZNcTbrtjXIq1qNzv2H5A6FiIiI4ihD7gDS1Zy/rgUAXL9gOm4qmSFzNERERBQrjjQNsgfWNSOvajXyqlZjb9dhucMhIiKiKDFpirMtyxbi09vPxTCVIqTNP3X30qadMkRGREREseD03CAYMUyF5jvPAwCsatyBSstHQe03PLsJNzy7CQDw4V/OQc6oYYkOkYiIiAaISdMgW1I0BUuKpmBv12Hct3YbHn2nLajdv/bp7xfNxflzJ8kRIhEREUWA03MJkpWZgb/8fCbaaxfh5oXHhbTf8Owm5FWthu6hd3HgkEeGCImIiKgvTJpiFE1xy4qzpqG9dhG23fFTXHLy1KC2xs87cPyfX0Ve1Wq0fbcv3uESERFRlJg0xSiW4pbDM5S445ez0VZzHn5zyrEh7Wfd8zryqlbj9pe3wufzxSNcIiIiihKTpiSgUCjw1/Nnob12Ed6qPCuk/V9vtyG/+hXkVa1Gy649MkRIREREXAieZKZoRqG9dhE8Xh+ue9qB/275OqhdvO9NAMAJE7Ox+voz5AiRiIhoSGLSlKRUSgUeuqQQAPC+tBsXNbwf1P7xl27kVa0GAKy98UzMmDAm4TESERENJUyaUsDJwli01y7CYY8XZ9/7BrY79we1n3P/kdGni+dPwV0XzIZCEVpYk4iIiGLDpCmFZKiUePN/a57+8+GXuP6ZjUHtz2zYgWc27AAAbLhlAcaPGZHwGImIiNLVkE6aLBYLAMBut6OkpASiKPa5PZn8Ys5E/GLORBw87MWMW/8b0j7/znUAgL8tnYtfzmPRTCIiolgN2aTJZrNBkiRUVlZCq9VCr9ejqamp1+3JaniGEu21iwAAL3/0Ja57Onj06fcrN+H3KzchM0OJxlvOliNEIiKitJASSZPD4UBpaWlI8iJJEiwWCwRBgCRJKCsrg1qtjuicoigGRpAkSUJRUVGf21PBz06ciJ+dOBG793ah8A5bUFvXYS9mLzuyrXKhCr89a7ocIRIREaWspE+a/EmRw+EIaes+CiRJEkpLS2E2mwd8DZPJBKPRGPH2ZDc2KzMw+vTIWxLuWP1JUHvdmm2oW7MNAPDZHeciM0OV8BiJiIhSTdInTTqdLux2SZKCXguCAJvth9EVi8USso//fIIgBF7X1dWhuro6aFtf21PN1WcIuPoMIezoEwAcd+urAIB79HOgK5yc6PCIiIhSRtInTb2x2WzQaDRB2zQaDRwOB7Raba/JVs9ziKIIrVYLi8USOKa37alsbFYmWu88F83Nzah9143Xt30b1P5H84f4o/lDAEDLnT9FhorF4omIiLpL2aTJ5XKF3e50OiM6XpIk6PV6CIIAl8sFURSh0+l63d5TV1cXurq6Aq/dbjcAwOPxwOPxDPwNJYDH44HX60XDJXOhUqnQ9t0+iPe/FbLftFuO3I1nKT8Z86aqExwlycXfP5K1/5K82D+oP6naRwYSb8omTb3pLZnqSRAEdHR0RLy9p5qaGixbtixke2trK7KysiKKIdG8Xi+cTidaWlqgVB4ZSXr18gL4fD48tGE3/vNJZ9D+OtORKuTnTBuDG087ikUz01y4/kHkx/5B/UnVPrJ3796I903ZpEmtVoeMKjmdzojvnotVdXU1brrppsBrt9uNKVOmoKCgANnZ2QmJYaA8Hg9aWlowbdo0qFTBi7/vnwHcD8D2yS6UPxm86H5tyx6sbTnyoOC3bv4xJqpHJipkSqC++gcR+wf1J1X7iH+mKBIpmzSJogiTyRSyPVElAjIzM5GZmYn6+nrU19cHhvdUKlVSdxalUtlnjAtnHYP22kXwen0Q/vRKSPsZd78BADhuwhisufHMQY2VEq+//kFDG/sH9ScV+8hAYk2d8TMET731vKvNX1MpUSNNfhUVFdi6dSvsdntCrzvYlEoF2msXob12Ef7ys5kh7Z99swd5VauRV7UaX3TsD3MGIiKi9JL0I002mw1WqxXAkXVExcXFgYXZZrMZBoMBxcXFsNvtUdVoov5deXo+rjw9Hx6vDwVhRp9ON74GAFh745mYMWFMosMjIiJKCIXP5/PJHUQq6z49t23bNnR2dib1mqbm5mZMnz495qHTDW1OLDG912v7lmULkZWZ9Dk5dRPP/kHph/2D+pOqfcTtdiMnJyeiz++Ump5LRuk6Pdef+fkatNcuQtOt4R9mPOu2NcirWo1nNmxPcGRERESDg0lTjOrr6zFz5kwUFxfLHYos/I9saa9dhBMmhmbo1c9vDqx9OuTxyhAhERFRfDBpitFQHWkKZ/X1Z6C9dhHevPkszJgQWqtq+i3/xfw7bXijRzVyIiKiVMBFJxR3U8eOwtobfwyv14dL/vUB3m3dHWjbtacLlz26IfC6+c6fYhgf2UJERCmASRMNGqVSgadLTwYAtH23D797xoEtO4OLiE2/5b8ozstF/a+1GD9mhBxhEhERRYRJU4x6Frek8PLHjcbLvzsDBw97MePW/wa12ds7MP/OdYHX0l3nQankI1uIiCi5cF4kRlzTNDDDM5SBheP3LZkTdh/hT68gr2o1djhZNJOIiJIHR5pINhdqJ+NC7WTs6zqME25bE9J+Rt2RopmXnDwVt58/iw8MJiIiWXGkKUZDveRAPIzOzAiMPi2afUxI+5Pvb0d+9ZHRJ9f+gzJESERExIrgcTOQiqJySaVqrfu6DuPx99pR9+pnYdsvOXkq7vjl7ARHld5SqX9Q4rF/UH9StY+wIjilvNGZGfjtT6ahreY8TNGMDGl/8v3tyKtajfrXWrC367AMERIR0VDDNU2U1BQKBd6qPBsAsGVnJ372j7eD2u9e8xnuXnNkNOqK0/Jw289PSHiMREQ0NDBpopQxa1IO2msX4bDHi0ffacOGtg7YPvkm0L7inXaseKcdAPDJX8/FyOGpMzxMRETJj9NzMeJC8MTLUClRdmYBHrmsCK/+/oyw+/zoL68ir2o1Xt3yVYKjIyKidMWkKUas0ySv44/ORnvtIrTVnBe2/ZonHcirWo1Pv3bD4+U9D0REFD1Oz1FaUCgUaK9dBAC4d+1n+Mf6lqD2c//2VuD/LdecgqI8TULjIyKi1MekidLOH845Dn845zgcOORB2b+b8Oa2b4PadcvfAwDMmZyDFytOY9FMIiKKCJMmSlsjhqnwxJXzAQDbvtmDc+5/M6j9wy86kV/9CgBgze/PxHFHj0l4jERElDqYNNGQMGPCGLTXLsLXnQdwcs26kPaFfzuSUP101tF46JLCRIdHREQpgAvBY8S751LL0TkjAgvHdYWTQ9r/u+Vr5FWtRl7Vauze2yVDhERElKyYNMWId8+lJoVCgXv0c9BeuwivXB++bEHhHTZc+q8P8NnXe8CnDREREafnaMibOfFI2YKuwx6cc/+b+Hz3/kDbW83fBabuskdkwPHnEmSo+LcGEdFQxKSJ6H8yM1R44+azAADSt3vx+5Wb8NEXnYF294HDmHbLfwEAt59/Ai49JU+OMImISCb8k5koDOGoLPznutMh3XUebl54XEj7n1/6OLD2qeuwR4YIiYgo0Zg0EfVBqVSg4qxpaK9dhAcunhd2n+NufTVQdZyIiNIXp+eIIvSLORPxizkT0bn/EOb8dW1Ie/eq4613nQeVkkUziYjSCUeaiAYoZ9QwtNcuQnvtIlw8f2rYfQr+9AryqlZj4/aOBEdHRESDhSNNMaqvr0d9fT08Hq5rGYpqLpyNmgtno/P7Q5izLHT06YIH3wUAFOflwnzNqYkOj4iI4ogjTTFinSYCgJyRwwJFM8PNytnbOwILx7/dw6KZRESpiCNNRHGkUCgg1SwCALz+2S5cviI0mS6+0wYA+H/XnY7Zk3MSGh8REUWPSRPRIPnJcePRXrsIBw97ce7f34T07b6g9p//820AwIXzJuHeJXOgUHDhOBFRMmPSRDTIhmcosf4PPwEAvPbpLlzxWPDo0/Mbd+L5jTsBAM9deyoKj81NdIhERBQBJk1ECXTW8UdGnzq/P4RH3pLwj/UtQe2LHzqycPwUYSyevPokli0gIkoiXAhOJIOckcPwh3OOQ3vtIrxYcVpI+3vS7kDZgpZde2SIkIiIeuJIE5HM5k5Ro712EVz7D2Le7Vb4fMHt4n1HHhg8YpgSn/z1XK59IiKSyZBOmiwWCwDAbrejpKQEoigGtqvValitVpSXl0MQBDnDpCFCPWo42moWwefz4Y/mj/Cc44ug9gOHvMivfgUAsPySQpw762g5wiQiGrKG7PSczWaDJEnQ6XQoLy+HwWAAALhcLtjtdoiiiOLiYhiNRpkjpaFGoVDg3iVz0F67CLabfhx2n2uebEJe1WqcdJctwdEREQ1dKTHS5HA4UFpaiqampqDtkiTBYrFAEARIkoSysjKo1eqIzimKYmBkSZIkFBUVAQDUanUgUfKPNBHJZdr4LLTXLoLX64Pwp1dC2r9xdyGvajUAYN0ffoyCo7ISHSIR0ZCR9EmTPylyOBwhbXq9PpBISZKE0tJSmM3mAV/DZDKFjCjZbDao1eqIkzCiwaRUKtBee6RopvHVT/HQ660h+yy49w0AwK9Pmoo7L5id0PiIiIaCpE+adDpd2O2SJAW9FgQBNtsPUxUWiyVkH//5uq9RqqurQ3V1dci6JVEUodFoUF5eDqvVGstbIIorw7nHw3Du8XAfOIQT/y/0eXdPfbAdT32wHQCwZdlCZGUm/T9zIqKUkLK/TW02GzQaTdA2jUYDh8MBrVbba7LV8xyiKEKr1cJisUCn06GhoQEulwuVlZVQq9VhEy8A6OrqQlfXD88Qc7vdAACPx5O0D+/1eDzwer1JGx8NzOhhSrTeeS4AoPK5zXjOsTNkn1m3rQEA3HH+Cbh4/pQ+z8f+QX1h/6D+pGofGUi8KZs0uVyusNudTmdEx0uSBL1eD0EQ4HK5IIoidDodlixZApvNBpvNBqvV2ut0X01NDZYtWxayvbW1FVlZybmuxOv1wul0oqWlBUrlkL0HIC2VnjgCpScW4Lt9h3GJ+fOQ9ltf+hi3vvQxckeq8NSSY6EMU7aA/YP6wv5B/UnVPrJ3796I91X4fD2rwiQnhUKB7qHW1dXBarUGTZ0VFBTAaDRGNMoUq3AjTVOmTIHT6UR2dvagXz8aHo8HLS0tmDZtGlQqldzh0CDy+XyYduuaPvf5T8WpOGHiD32V/YP6wv5B/UnVPuJ2u6HRaNDZ2dnv53fKjjSp1eqQUSWn05mwhduZmZnIzMwM2a5SqZK6syiVyqSPkeLDv3D8A2k3lja8H9L+i/ojj2wZPVyFLcsWAmD/oL6xf1B/UrGPDCTW1Bk/68FfLqAnf+mARKmvr8fMmTNRXFyc0OsSReokYSzaaxeh9a7zwrbvO+hBfvUrKLjlVTTv7gq7DxERpVjS1H0dU8+73fy1lhJdIqCiogJbt26F3W7vf2ciGan+V7agvXYR7tXPCbvP7/7fFyi45VWcUrMOKTJzT0SUMEk/PedfkA0cWXxdXFwcWLNkNpthMBhQXFwMu90eVY2mWNXX16O+vj7l7hagoW1x4WQsLpwM576D0N4eWlLjq84DgUe2fHjbOcgZOSzRIRIRJZ2UWQie7NxuN3JyciJaSCYXj8eD5uZmTJ8+PaXmm2nw+Xw+XLHCjte3fdvrPrcu+hGuPoPPYRyq+PuD+pOqfWQgn99JP9JERINPoVDgX5cVorm5GVnjp+D0utdD9rlj9Se4Y/UnAIBtd/wUwzNSanafiChm/K0XIy4Ep3RzTM4ItNcuQltN+IXjADDj1v8ir2o13mrufWSKiCjdcKQpRhUVFaioqAgM7xGlC4Xih+fdvbRpJ254dlPIPpf+awMAYH6+BqvKT0lkeERECcekiYj6df7cSTh/7iTsOXAIs8M8725DmxN5VasBAGtvPBMzJoxJdIhERIOO03Mx4vQcDSVjRgwLTN0tPGFC2H3Ouf9N5FWtRv1rLSxbQERphXfPxQnvnqNUF23/2OHcjzPqXutzn+Y7f4phKv6Nlsr4+4P6k6p9ZCCf3/wtRkQxmaIZFag4fvzR4aflpt9yZOH4m32UNCAiSnZc00REcaFSKvDq788EANjbndAvfy9kn988emTh+PFHj8Hq68+ASqlIaIxERLHgSFOMuKaJKFRxngbttYvw6e3nhm3/9Os9KPjTK8irWo3Pvt6T4OiIiKLDpClGfPYcUe9GDFMFnnd335Lwz7tb+LcjC8cfWNfMheNElNQ4PUdECXGhdjIu1E7GLvcBzL9rXUj7fdZtuM+6DQCwZdlCZGXy1xMRJReONBFRQo3PPlJx/LM7zsUv5kwMu8+s29Ygr2o1/vPhlwmOjoiod0yaYsQ1TUTRycxQ4YGL56G9dhFWlp0cdp/rn9mIvKrVuOzRDeg67ElwhEREwZg0xYhrmohid5IwNjD6NC4rM6T9jW3f4rhbX0Ve1Wq07OLCcSKSBxcNEFHSyMxQofFWEQDw+LvtuO0/H4fsI973JgBg2S9OwK9OmsqimUSUMPxtQ0RJ6bJT89BeuwhN/0uierrtPx9j+i3/hXjfG+jYdzDB0RHRUMSkiYiS2tisTLTXLkLLnT/FdWdNC2lv2bUX8263Iq9qNd6XdssQIRENFUyaiCglZKiU+OPC49B613l4pjT8wvGLGt5HXtVq1L36KReOE1HcMWmKEe+eI0oslVKBUwqOLBzf/H/n4MTJOSH7PPh6a2Dh+E7X9zJESUTpiElTjHj3HJF8xowYhv9cdzraaxfh/qXhK46fVrseeVWr8dKmnQmOjojSDZMmIkoLF8ybjPbaRXjz5rPCtt/w7CbkVa3GdU87OHVHRFFhyQEiSitTx45Ce+0iHPZ48eeXPsYzG7YHtb/80Vd4+aOvAABvG87CJPVIKBQKOUIlohTDkSYiSksZKiVqLpyN9tpFWH5JYdh9Tje+hvzqV7DSvj1sOxFRd0yaiCjtnTvr6MDC8YKjRoe0G57bjLyq1bjyMTvcBw7JECERpQJOzxHRkDFmxDCs+8NP4PP58Og77bj95a1B7es/3YUT/2/tkf//w4+RP240p+6IKIAjTUQ05CgUClx1ej7aas6D9cYzw+5z9r1vIL/6Fby0aSd8Pl+CIySiZMSkiYiGLIVCgekTxqC9dhG2LFuIX8yZGLLPDc9uQn71K7j1xc04cIh33RENZUyaYsTilkTpISszAw9cPA/ttYtwozgjpP3J97fj+D+/ih/f/Rqkb/fKECERyY1JU4xY3JIo/dwgTkd77SL8v+tOx1TNqKC2z3fvx9n3voG8qtX45/pmHPJ4ZYqSiBKNSRMRUS9mT87Bm5VnYeOfS3DhvEkh7fes3Ybpt/wXt7ywGS279sgQIRElEpMmIqJ+5I4ejvuWzkVbzXm464LZIe1PfbAd4n1v4i8vbcEnX7m5cJwoTbHkABFRhBQKBX510lT86qSp6Pz+EK572oG3mr8LtD/x3ud44r3PAQAP/VqLn84+Rq5QiWgQMGkiIopCzshh+PdVJ6Fz/yGYm3bgjtWfBLVf+5QDAGA493j87MRjMKXH2igiSj2cniMiikHOqGG4+gwBbTXnYcXlxZikHhnUbnz1U5xR9xr++v+2smQBUYqL20jT+vXr4XA4YLfbIUkSXC4X1Go1BEGAIAgoKCjAkiVLkJ2dHa9LxsxisQAA7HY7SkpKIIpiULter4fZbJYjNCJKMQqFAmcdPx7vVJ2N9u/2oezfjdj2zQ+lCR59pw2PvtOGM6aPw1/Pn4VJ6pEYnsG/W4lSScxJU1VVFRoaGjB27FgsWLAARUVFEEURgiDA5XLB6XTC5XJh7dq1qK2tRWFhIcrLy3H22WfHI/6o2Ww2SJKEyspKaLVa6PV6NDU1BbU7HA4ZIySiVJU3bjTW3vhj7NpzAO9LTlz/zMZA21vN3+Gse17HcRPGoE53IqaNz8LoTK6UIEoFUf9LXbduHQwGAy666CK0tbUhJycnouM2btyImpoaGI1GmM3miEaeHA4HSktLg5IaAJAkCRaLBYIgQJIklJWVQa1WRxSHKIqBkSVJklBUVBRoc7lcAABBECI6FxFROOPHjMAv5kzEotnH4LmmL3C/bRu+6jwAAPjsmz04v/4dAMD1C6bjvNlH4/ijk2cknohCRZU0Pfzww5AkCY2NjQM+dt68eVi1ahUkSYJOp0NDQwPy8vJ63d+fFIUb9ek+OiRJEkpLS6OaTjOZTDAajYHXNpsNOp0uaBsRUbRUSgWWFE/BkuIp2PxFJ1a824bnHTsD7Q+sa8YD65px+al5uOK0PBw7drSM0RJRbwacNLW1tUEQBJSWlsZ0YUEQsHbtWtx99924+eabe91Pp9OF3S5JUsj5bDZb4LXFYgnZx3++7iNIdXV1qK6uDmxzOBwha5uIiOJl9uQc3LdkLs6bdQzad+/D/dZt2HfwyALxx95tx2PvtuOcmROgL5qCkpkTZI6WiLobcNKUn5+P/Pz8uAXQV8LUF5vNBo1GE7RNo9HA4XBAq9X2mmz1PIcoitBqtbBYLIFj/MmXJEloaGhAWVlZVDESEfVG/F9CdPH8qWj9di9++5QDX3R8DwBYu/UbrN36DfLGjsKiE4/BtT+ZhiyueyKSXcr+K/SvO+rJ6XRGdLwkSdDr9YEF66IoQqfTQavVQqvVBo1ahdPV1YWurq7Aa7fbDQDweDzweJLztmKPxwOv15u08ZG82D/kMSJDgROOGYPX/3AmNu9043fPbgokT+2796P+tVZ8tMOFm0pmYFzWcEzsUdIgUdg/qD+p2kcGEm/KJk296S2Z6kkQBHR0dPTaLooiWltbe22vqanBsmXLQra3trYiKysrohgSzev1wul0oqWlBUolb3WmYOwf8hsJ4JHzJ8Lj9eHO17/Bu9v3AQDeatmNt1rew4gMBe47bxKm5AzHMJUiobGxf1B/UrWP7N27t/+d/ieuSdMjjzyCq6++Op6n7JVarQ4ZVXI6nRHfPRer6upq3HTTTYHXbrcbU6ZMQUFBQVLVourO4/GgpaUF06ZNg0qlkjscSjLsH8nl38fNAACYG79A1QtbAAAHDvvw2/98AQCYMzkHT15VjFHDE/O3L/sH9SdV+4h/pigSEf9r27RpU79TX1arNWFJkyiKMJlMIdu7lw4YTJmZmcjMzER9fT3q6+sDw3sqlSqpO4tSqUz6GEk+7B/J56KTjkVRvgZeH1D+7ya0fXdk9OnDLzoxe5kNpxaMxW9OycPCEyZAoRjc0Sf2D+pPKvaRgcQacdJ01113weFw9DqS43K50NbWFvGFo+GvMg6E1lDy11pK1EiTX0VFBSoqKuB2uyOuVUVENBDTxo8BALzw21Px4sad2O78Ho++c+T37butu/Fu627cfv4JuFA7GSOGqaBSJnbqjmioiDhpMhqNkCQJCxYs6HWfu+++Oy5BdWez2WC1WgEcWUdUXFwcuMvNbDbDYDCguLgYdrtdlkee9BxpIiIaLOpRw3H5aUfuXhZnjsedqz/Bx18emVr480sf488vfYzjjx6DOy+YjcJjc+UMlSgtKXw+ny/SndetW9dn0tRfezrzjzR1dnYm9Zqm5uZmTJ8+PaWGTikx2D9SU9dhD5aY3seHO1xB22//5SwAQOHUXMycGPvvJPYP6k+q9pGBfH4PaAVhfwnRUE2YiIjkkpmhwksVp2HbN3vwXutuPPFeO1q/3Yc/v7jlf+1K3H7+LCwpniJzpESpL+1KDiQap+eIKBnMmDAGMyaMgThzAn77ZBM+/KITANB12IvK5z6C9N0+/HzOMQCAEyZy/SVRNGIupFBdXR2POFJWRUUFtm7dCrvdLncoRESYpB6Jl647HZ/efi4+XrYQi7WTAQDL32jFogfexqIH3kb9ay0yR0mUmmJOmvoqEElERPIYMUyF0ZkZuPOCWTh92jiMHT080Hb3ms9w64ubMYAlrUSEOEzPDXZdkGTH6TkiSmYjhqnw5NUnAQAOe7x4YF0zHljfgiff346Pv3TjzOlHYVLuSPxkxlEYnz1C5miJkhvXNMWIdZqIKFVkqJS46ZzjkDlMhbvXfIaN213YuN0FABDGjcaL152G7BHD5A2SKIkxaSIiGmKu+XEBsjIz0PrtXnz8pRtNn3dA+m4fSu57A8svKcS8qazxRBQOkyYioiFGpVTgslPzAq8/3OHCFY/Z8Y27C1c/3oiKs6Zh5sRsnJSvGfJLMIi6Y9IUI65pIqJUN2eKGq/f/BNc+OC7aNm1F399eSsAYHiGEvPzNPj7RXOhHsmPC6KY754b6ndfsOQAEaWD7BHD0HBpIc6YPg5zJudgeIYSBw978XbLd7j6iUYcOMQ/DIli/tOhsLAwHnEQEZHMhKOy8O+rjtxp59x3EO9Lu1H13EfYuN2FXz74Hv5wSi6myxwjkZxiHmkqLS2NRxxERJRENKOH47zZx+CRy4oxJjMDzbv24i+2r+Daf1Du0IhkE3PS1B+32z3YlyAiokEyP1+DdX/8MaZqRmHXvsMovHM9Tq1Zh5c27ZQ7NKKEi2vSNHbsWLz22mtB2+x2O6qqqlBUVIQZM2bE83JJob6+HjNnzkRxcbHcoRARDYrxY0bgnxfPxahhR+6k+7LzAG54dhOK77Sh7IlG7HIfkDlCosRQ+OK4krugoADt7e0wGAy46667QtqnTZuGlpb0fOaRv7hlZ2cnsrOz5Q4nLI/Hg+bmZkyfPh0qlUrucCjJsH9QXzweD+ybP8WosRNh/WQX/vlaC/yfHrMmZcNyzakYMYz9ZihL1d8hA/n8jutIU1VVFdasWYPly5dj/vz5+Pzzz4PadTpdPC9HREQJpB6hwgkTs/GHc47DmzefhYd/UwTN6OHYstONW17YAq93aN9NTekvrkmTQqGAKIqQJAk5OTkQBAEvvPBCoH3s2LHxvBwREclkimYUSmZOwD8ungelAnjO8QUueOhdtH23T+7QiAZNXJMml8sFAFCr1bBaraipqcHixYtx0UUXAeDDfYmI0s1p08bhrgtmY/RwFT7c4cIlj3wA5z7eYUfpKa5Jk9VqDXpdWVmJxsZGNDY2Yvr06SwASUSUhi6aPxXr/vAT5I0dhZ2u7/FH84ecqqO0FPekaeHChXjkkUcCpQa0Wi1aWlowd+5cWCyWeF4uKfDuOSIi4OicEXjw14UYnqHE+k934ZYXN2PTDheTJ0orcX2YkNfrBYCwK9DNZjPWrVsXz8slhYqKClRUVARW3xMRDVUzJ2bjtp/PxC0vbMEzG3bgmQ07cGrBWDx0SSFyRg6TOzyimEU90rRp06Ze23pLHrgQnIgovf1q/lTcq5+Ds48fj+EqJd5t3Y1LHvkAHVznRGkg6qTJZDIl5BgiIkodCoUCiwsn49HLi/FixWkYO3o4Nu/sRMn9b+D6ZzbixY07h/yD3il1RT09Z7Va8a9//Qu5ubkR7e90OmGz2aK9HBERpZiZE7OxsvxkXPLIBnztPoD/fPgl/vPhl7C3O3HHL2fxjmpKOVEnTZIkoaysLOQvBoVC0etfEfwHQkQ0tEwbPwbr//hjfNDmxPvSbjS8KeGpD7Yjf9xoXH2GIHd4RAMSddIkCAKMRiPUanVgm8/nQ11dHQwGQ8j+HR0dqK6ujvZyRESUokYNz8BZx43HWceNx4QxI/DXl7firlc+wfFHZ+P06ePkDo8oYlEnTaIoYvHixSHbn3vuOSxYsCDsMZyeIyIa2q44LQ9bv3LD0vQFKp524MWK05A/brTcYRFFJOqF4OFGkwbjmGTHOk1ERJFTKBS445ezMGeKGp3fH8LP//E2nt2wnYvDKSVEnTTl5+cn5JhkV1FRga1bt7LaORFRhEYMU6Hh0kLMm6rG3q7DqHp+MyotH7EQJiW9qJOm9evXh93e118LvR1DRERDy4TsEbBccyr+dN7xUCkVMDd9gX+sb5E7LKI+RZ00mc3msNv7ukOut2OIiGjoUSkVKDuzADUXzgYA/GN9M7bs7JQ5KqLeRb0QfOXKlcjNzYVGowna3tjYiHvuuSdk/927d2PVqlV46KGHor0kERGlIX3hZLz26S78d8vX+NMLm/HCb0+DSskSNZR8ok6aXC4Xamtrw7Y1NTWF3c46TURE1JNCocCyX5yAt5u/w0dfdOKpDz7Hb07JkzssohBRJ01arRYPP/xwUJ2mvnR0dKCsrCzayxERURobnz0CN597HP7y0se4+9XPsPCEozEhe4TcYREFiTppKioqwrx58yLePz8/H0VFRdFeblBYLBYAgN1uR0lJCURRBHCk2rkgCHC5XAAQcWJIRETR+/VJx+K5pi/w4Red+OvLW1H/K63cIREFiXohuNFoTMgxg8Vms0GSJOh0OpSXlwfVkCovL0dhYSFqampkjJCIaGhRKRW484LZUCqA1R99BevWb+QOiShI1CNNOTk5CTkGABwOB0pLS0PWSkmSBIvFAkEQAs/Ci3RUSBTFoJGl7qNg5eXlEAQBgiBwlImIKIFmTcrBFafl419vt+G6px1Yfkkhzjp+vNxhEQGIYqSpra0NmzZtilsAzz//fJ/t/ik0h8MR0qbX61FZWQmdTgedTofS0tKoYjCZTEEjTf7puYaGBj76hYgowSrPPQ7ij8aj67AXVz/RiGc2bJc7JCIAUYw05efno6qqCo2Njbj66qtjuviSJUvwpz/9qc99dDpd2O2SJAW9FgQhKMGxWCwh+/jPJwg/PFm7rq4O1dXVQdsqKysD++r1+l7vBiQiovjLzFDhwV8Xour5j/C8Yyeqn9+ML13f46aSGbwLm2QV1fRcbW0tHn74YSxcuBDl5eW48MILIz7W7XbDZDLBZrNh+fLlUT9axWazhdSI0mg0cDgc0Gq1vSZbPc8hiiK0Wi0sFgt0Ol1grVNZWVnI+YmIKDGGZyhxr34OJueOwgPrmvGP9S3Yve8g7vzlLCZOJJuo1zSVlpZCFEUYjUYYDAZotVqUlJRAEARoNBqo1Wq4XC44nU64XC7Y7XZYrVZ0dnaivLwca9asiSlw/51tPTmdzoiOlyQJer0+cJecKIrQ6XQQRRE2mw02mw1WqxUPP/xw2OO7urrQ1dUVeO12uwEAHo8HHo9nYG8mQTweD7xeb9LGR/Ji/6C+yNU/bji7AMdkZ+KWF7fg6Q+2Qxg7CleclpfQGCgyqfo7ZCDxRp00AUem6pYvXw4AeO6552C1WrFq1SpIkhRIXjQaDQRBgFarhdFoxIIFC2K5ZL96S6Z6EgQBHR0dYdv8C8T9/w2npqYGy5YtC9ne2tqKrKysiGJINK/XC6fTiZaWFiiVUd84SWmK/YP6Imf/mJcDlBaNhcm+G7WvfopJw/YiPzczoTFQ/1L1d8jevXsj3lfh6+sJu0lEoVAEPQy4oaEBJpMpaL1Rbm4uzGZzn8lOvIQbaZoyZQqcTieys7MH/frR8Hg8aGlpwbRp06BSqeQOh5IM+wf1Re7+4fP5cM2TG2H7dBdmTczGc9ecjAxV6nwwDwVy95Foud1uaDQadHZ29vv5HdNIk5xEUYTJZArZnqgCmpmZmcjMzER9fT3q6+sDw3sqlSqpO4tSqUz6GEk+7B/UF7n7x10XzsaG+97Ali/dePTd7bj2JwWyxEG9k7uPRGMgsaZUmt596q373W7AD7WWEl1XqaKiAlu3boXdbk/odYmIhprx2SPwl5+fAAC437YN0reRT6sQxUPSJ002my1QQ6mmpiZQtwkAzGYzDAYDLBYLTCYTzGazXGESEVECLNZOwhnTx+HgYS+qntsMrzclVphQmkiZNU3Jqvv03LZt2yKaE5WLx+NBc3Mzpk+fnlJDp5QY7B/Ul2TqHzuc+7Hwb29i/0EPbv/lLFx68rGyxkNHJFMfGQi3242cnJyIPr+TfqQp2XF6jogosaZoRqFy4XEAgNpXPsEO536ZI6KhgkkTERGlnN+ckofCY3Ox76AHpU80Yl/XYblDoiGASVOM6uvrMXPmTBQXF8sdChHRkKFUKvDAxfMwLisTn369B9c+5cCBQ6lVVJFSD5OmGHF6johIHpPUI/HwbwoxcpgKb277Fpc9ugGd+w/JHRalMSZNRESUsuZNzcUTV81HVmYGPmhz4uf/fBtbv3TLHRalKSZNMeL0HBGRvIrzNFhZfjIm547Edud+XPjQO3hp0065w6I0FFVF8PXr18ctAI1Gg7lz58btfIlWUVGBioqKwC2LRESUeCdMzMHLvzsdNzy7CW9s+xY3PLsJH3/pRuXC4/i4FYqbASdNbW1tWLt2bdwCGDduXEonTURElBzUo4bj0cuLce/az/Dg661oeFPCJ1+58c9faZEzcpjc4VEaGHDSlJ+fj9ra2sGIhYiIKCYqpQKV5x6PmROzcbP5I7zV/B0ubngfT1w1H+OyMuUOj1IcxyxjxDVNRETJ52cnToTl2lMwLms4tn7lxhLTe/i684DcYVGKY9IUI5YcICJKTidMzMGq8lMwMWcEpG/34ZJ/fYCOfQflDotSGJMmIiJKW8JRWVh1zSk4OnsEWnbtxZWP27H/IKuHU3SiunsunPXr18PhcMBut0OSJLhcLqjVagiCAEEQUFBQgCVLliTtw2yJiCg9Tc4dhSeumg/98vewcbsLNzy7CaZLCqFUKuQOjVJMzCNNVVVV0Gg0KC8vR0tLC4qKilBWVobly5ejqqoKoihCo9Fg7dq10Gq1WLp0aVxLFsiNa5qIiJLfjAlj8OjlRRiuUsK69Rv8bV2z3CFRCop6pGndunUwGAy46KKL0NbWFnGNoo0bN6KmpgZGoxFmsznlR55Yp4mIKDUUHqvBXRfOxh/NH+KBdc2YecwYnDvrGLnDohQSVdL08MMPQ5IkNDY2DvjYefPmYdWqVZAkCTqdDg0NDcjLy4smDCIiogHRFU7G1i/dePSdNty06kPkjRuN449O7T/eKXEGPD3X1tYGQRBQU1MT04UFQcDatWthNptjOg8REdFA/Om843HatLHYf9CDa/7dhM7v+ZBfisyAk6b8/HwsWLAgbgHcfPPNcTsXERFRfzJUSvzjYi0mqUeiffd+/GHVJni9PrnDohQQ00LwTZs2xSkMIiKixNGMHo7llxRieIYStk924Z+vtcgdEqWAqJOmu+++G2effTauvfbaeMaTcnj3HBFRapo9OQd3/HIWAOB+2za89tkumSOiZBfTSNO6detQWVkZeH3PPffEHFCqYUVwIqLUtaRoCn510lT4fMANz2zE9t375Q6JkljUSdPu3buxatUquN3uwLbW1ta4BEVERJQot/18JuZOUcN94DDKn2zC9wc9codESSrqpKm2tha5ubmYN28eVCoVFi5ciMbGRrzwwgtBiRQREVEyy8xQ4aFLtBiXNRyffOXGn17YDJ+PC8MpVEzTc5WVlfB6vVi5ciXy8/PR1NSExYsXIzc3F2PHjsXChQtRXV3NRIqIiJLaMTkj8Y+LtVApFXhh40488d7ncodESSguD+zV6XRYvnw5ysrK4PV6YbfbUVtbi/z8fFit1qBE6tprr+Vdd0RElHROKRiL6p8eDwC4/eWtsLc7ZY6Ikk1ckia/kpISAIBWq0VpaSmWL1+OxsbGoETK5/NBp9Nh4cKFHH0iIqKkctXp+fjZicfgsNeH3z7lwC73AblDoiQS16Rp8eLFvbZ1T6RaWlqwfPlyGAyGeF6eiIgoJgqFAsbFJ2LGhCx8u6cLv33KgYOHvXKHRUkirklTpB555BE0NDSgsLBQjsvHFes0ERGll9GZGTBdWoQxmRlo/LwDd73yidwhUZKQJWnKzc1Nm2fOsU4TEVH6yR83GvctnQsAeOzddryw8Qt5A6KkIEvStHjxYrS0tODqq6+W4/JERET9Kpk5Ab87exoAoPr5zfj4y06ZIyK5yZI0ERERpYLfizPw4xlH4cAhL655sgmu/QflDolkNOCkqa2tLa4lA55//vm4nYuIiCieVEoF/n7RXEzRjMQO5/e44dlN8HhZ+HKoGnDSlJ+fj2effRaPPPJIzBdfsmQJBEGI+TxERESDRT1qOJZfUojMDCXe2PYt/m7bJndIJJOopuf89ZYWLlw44JEit9uNu+++GwsXLoTRaMTcuXOjCYGIiChhTpiYg5oLZwMAHljfAtvWb2SOiOSQEe2BpaWlEEURRqMRBoMBWq0WJSUlEAQBGo0GarUaLpcLTqcTLpcLdrsdVqsVnZ2dKC8vx5o1a+L5PoiIiAbVhdrJ2LTDhSfe+xw3rtyEF687DQVHZckdFiVQ1EkTcGSqbvny5QCA5557DlarFatWrYIkSXA6j5Sf12g0EAQBWq0WRqMRCxYsiD3qOLFYLAAAu92OkpISiKIIAHC5XKipqUFJSQk0Gg20Wq2cYRIRUZK4ddFMfPylG02fd+DKx+x4/tpTMTYrU+6wKEFiSpq6W7x4cZ8VwZONzWaDJEmorKyEVquFXq9HU1MTAECv18NqtQIA6urqmDQREREAYHiGEssvKcQFD76Dz3fvR9m/m/DU1SdhxDCV3KFRAqREyQGHwxG2ergkSairq4PFYkFdXR1cLlfE5xRFEZWVlYHzFBUVATiSTAmCAEmS4HK5AvsQEREBwFFjMrHi8mKMGZGBps878Efzh/DyjrohIeFJ00DLFfin0BwOR0ibXq9HZWUldDoddDodSktLo4rJZDIFnoMnSRIkSQJwJIFqaGiI6pxERJS+pk8YA9MlhRimUuDlj77CPWs/kzskSoC4Tc9FymQy4aGHHop4f51OF3a7P7HxEwQBNpst8NpisYTs4z9f9zIHdXV1qK6uDtqm1WohCAIEQUBubi7KysoijpeIiIaGU6eNQ82FJ+KP5g/x4OutmKoZhYvmT5U7LBpECU+afL74DGHabDZoNJqgbRqNBg6HA1qtttdkq+c5RFGEVquFxWKBTqeDKIqBtU0ulyvkGkRERH66wsnY7tyPB9Y145YXt+AY9Uj8eMZRcodFgyThSVO89LZ+yX/XXn8kSYJer4cgCHC5XBBFMTAKVVhYGBip6u3Bwl1dXejq6gq8drvdAACPxwOPxzOwN5MgHo8HXq83aeMjebF/UF/YP3p3/VkCtu/ehxc3fYnfPtmEZ0tPwsyJ2XKHlXCp2kcGEm/CkyaFQjGo5490MbggCOjo6AjbFsl0XE1NDZYtWxayvbW1FVlZyVm3w+v1wul0oqWlBUplStwDQAnE/kF9Yf/o21WzR6Dt6xH48OsDuHzFB/jbosk4anTKjktEJVX7yN69eyPed1B/op2dnTAYDIFEyefzYd26dbj22msDrxUKBYxGI7KzB5aVq9XqkFElp9MJtVodl9j7U11djZtuuinw2u12Y8qUKSgoKBjwe0kUj8eDlpYWTJs2DSoVb4+lYOwf1Bf2j/49dmw+ljR8gOZde3H7m7uxquwkjBkxTO6wEiZV+4h/pigSg5o05eTkBIpf+l1zzTUDWgjeG1EUYTKZQrb7SwcMtszMTGRmZqK+vh719fWB4T2VSpXUnUWpVCZ9jCQf9g/qC/tH33KzVHjsyvm4oP4dbPtmLyqe2YQVl8/H8IzUGXWJVSr2kYHEmvCfZCzTc92n3no+6NdfaylRI01+FRUV2Lp1K+x2e0KvS0REyWeSeiQevbwYo4er8E7LblQ9/1HcboAi+SV9+muz2QI1lGpqagJ1mwDAbDbDYDDAYrHAZDL1umh7MNXX12PmzJkoLi5O+LWJiCj5zJqUg3/+WguVUoHnHTtxv61Z7pAoThS+BKfA11xzTciUXTpwu93IyclBZ2dnUq9pam5uxvTp01Nq6JQSg/2D+sL+MXDPbNiO6uc3AwDqFp+IJcVTZI5ocKVqHxnI53fCR5rSMWEiIiLq6eL5U3HdWdMAANUvbMYb276VOSKKVdJPzyU7Ts8REVFv/nDODFwwbxI8Xh9++2QTPv6yU+6QKAaDnjR1dnbi7rvvxiOPPDLYl5IFF4ITEVFvFAoFjItPxMmCBvsOenDlY3Z86fpe7rAoSnFPmu655x5UV1fjhRdeCMwT3nzzzdDr9WmbOBEREfVmeIYSpkuKMH18Fr5xd+HKx+xwHzgkd1gUhYiSpvb29ohOtmTJEqxduxZNTU1YvHgxcnNzsXDhQrzwwgvIycmJJc6kxek5IiLqT86oYVhxRTGOGpOJT7/eg98+6cDBw165w6IBiihpslqtWL9+fb/7CYKAtWvXYu3atfB6vVizZg1ycnKwePFiqFSqwINw0wmn54iIKBKTc0dhxeXFGDVchbdbvkP185tZwynFRFQRvLS0FMXFxf0mBj1HW0RRhCiKAI6sbUrX0SYiIqJIzJqUg/pfa3H14414zvEFJueOxI0lM+QOiyIU0UjT+vXrkZ+fH7StuLgY1157LV577bXANoVC0eszXJgwERERAWcdNx63nz8LAPD3dc1Y1bhD5ogoUhElTQaDAXV1dUHbWltbYbfbsWDBAqhUqsBIlF6vx+effz4owSYjrmkiIqKB+tVJU1FxVgEA4E/Pb8ZbzazhlAoiSpp8Ph/y8vKCtomiiMbGxsDapQULFsBsNsNqtUIQBCxcuBD33HMPNm3aNAhhJw+uaSIiomj88ZzjcP7ciTjs9eHaJx345KvwMzWUPCJKmoxGY0i5gFWrVgX+XxRF1NbWoqWlBR0dHVi5ciXy8/OxfPlyaLVajB07FkuXLo1oMTkREdFQoFAoUKc7ESfla7C36zCuWGHHV52s4ZTMIkqaFixYAL1eH9EJc3JyoNPpsHz58kASVVNTA5/PB6PRGFOwRERE6SQzQ4WGS4swbXwWvnYfwBUrWMMpmUVc3DLahdw5OTkoKyvDqlWrsGbNmqjOQURElK5yRg3DYz1qOB3ysIZTMuKz52LEheBERBSrybmj8OhlrOGU7Jg0xYgLwYmIKB5mT85B/a+0UCoAS9MX+Pu6ZrlDoh4iKm7ZUzwXdGs0GsydOzdu5yMiIkpVZx0/Hnf8cjb+9MJm/M3WjEnqkdAXTZE7LPqfASdNbW1tWLt2bdwCGDduHJMmIiKi//nVSVOxo2M/Hnq9FdXPb8bROSNwxvSj5A6LEEXSlJ+fj9ra2sGIhYiIiADcfM5x2NnxPf7z4Ze49kkHVpWfgpkTs+UOa8jjmqYYcSE4ERHFm1KpwN36E3Gy8L8aTo9twJcu1nCSG5OmGHEhOBERDYbMDBVMlxZh+vgsfOPuwuUrNqDze9ZwkhOTJiIioiSVM3IYHrtyPsaPycS2b/bimn83oeuwR+6whiwmTURERElsknokVlxRjNHDVXhP2g2D5SPWcJIJkyYiIqIkd8LEHDx0SSEylAq8uOlL3L3mM7lDGpKYNBEREaWAM2cchZoLZwMAHny9FU++/7nMEQ09TJqIiIhShL5oCm4UZwAA/vLSFti2fiNzREMLkyYiIqIUcv2CaVhaNAVeH/C7Zzbiwx0uuUMaMpg0ERERpRCFQoE7LpiFM2cche8PeXDV43Zs371f7rCGBCZNMWJxSyIiSrRhKiUe/LUWJ0zMxnd7D+KyFRvg3HdQ7rDSHpOmGLG4JRERySErMwMrLi/GJPVItH23D1c/bseBQ6zhNJiYNBEREaWo8dkj8NgVxcgekQHHdhd+/+wmeLys4TRYmDQRERGlsOkTxuDh3xRhuEqJVz/+Gnes3ip3SGmLSRMREVGKO0kYi3uWzAEArHinHY+8JckcUXpi0kRERJQGfjFnIqp/ejwA4M5XPsErm7+SOaL0w6SJiIgoTZSdKeA3pxwLnw/4/cpNsLc75Q4prQzppMliscBiscBgMMBmswW26/V6uFwu+QIjIiKKgkKhwG0/PwElMyfg4GEvrn68ES279sodVtoYskmTzWaDJEnQ6XQoLy+HwWAAAEiSBJvNhsLCQhQUFCA3N5cJFBERpQyVUoEHLpqHeVPV6Pz+EC5fsQG79hyQO6y0kCF3AJFwOBwoLS1FU1NT0HZJkmCxWCAIAiRJQllZGdRqdUTnFEURoigGzlNUVAQAcLlc6OjoCFwXQMTnJCIiSgYjh6vwyG+KsPihd9G+ez+ueqwRz5adjNGZKfGxn7SSfqTJYrEA+CGB6U6v16OyshI6nQ46nQ6lpaVRXcNkMgVGmrRabWB7Y2Nj0GsiIqJUMTYrE49dMR+a0cOxeWcnrnvagcMer9xhpbSkTzl1Ol3Y7ZIUfDulIAhB65IsFkvIPv7zCYIQeF1XV4fq6uqgbQBgMBhQXV0dS+hERESyyhs3Go9cVoRfPfw+XvvsW/z5pS2464LZUCgUcoeWkpI+aeqNzWaDRqMJ2qbRaOBwOKDVantNtnqeQxRFaLVaWCyWoGNsNhuMRmPc4yYiIkok7dRcPHDRPFzzZBOe2bADk9Qjcd3Z0+UOKyWlbNLU2+JspzOy2yslSYJer4cgCHC5XBBFMaJEy6+rqwtdXV2B1263GwDg8Xjg8STns388Hg+8Xm/SxkfyYv+gvrB/pLYFxx+Fv/zsR/i///cJ7lm7DUdnZ+KCeZPieo1U7SMDiTdlk6beRHqnmyAIgQXf4fRcdN5TTU0Nli1bFrK9tbUVWVlZEcWQaF6vF06nEy0tLVAqk345GyUY+wf1hf0j9Z08FtDNUsOyxQXDc5txyP0d5k0cFbfzp2of2bs38pIMKZs0qdXqkFElp9OZsDvdqqurcdNNNwVeu91uTJkyBQUFBcjOzk5IDAPl8XjQ0tKCadOmQaVSyR0OJRn2D+oL+0d6qCmYhgPKD/HyR1/jzje+xaqyk3Dc0WPicu5U7SP+maJIpGzSJIoiTCZTyHZ/6YDBlpmZiczMTNTX16O+vj4wvKdSqZK6syiVyqSPkeTD/kF9Yf9IfSoVcO+Sufh2zwZ80ObElY834YWKU3FMzsi4nD8V+8hAYk2d8TMET731vNvNX2sp0TWVKioqsHXrVtjt9oRel4iIKBqZGSo0XFqE6eOz8LX7AC5/1A73gUNyh5USkj5pstlsgRpKNTU1gbpNAGA2m2EwGGCxWGAymWA2m+UKk4iIKGXkjBqGFVcU46gxmfjsmz249skmHDzMGk79Ufh8Pp/cQaSy7tNz27ZtQ2dnZ1KvaWpubsb06dNTauiUEoP9g/rC/pGetuzsxFLTe9h30IML503CvUvmRF3DKVX7iNvtRk5OTkSf30k/0pTsOD1HRESpatakHDx4SSFUSgWe37gT967dJndISY1JU4zq6+sxc+ZMFBcXyx0KERHRgP14xlGouWA2AOCfr7Xg6Q+2yxxR8mLSFCOONBERUapbUjwFNyw4UiX8zy9twfpPv5E5ouTEpImIiIjwe3E69IWT4fH6UPHURnz0hUvukJIOkyYiIiKCQqHAXRfOxhnTx+H7Qx5c+ZgdO5z75Q4rqTBpihHXNBERUboYplLiwV9r8aNjsvHd3oO4bMUGdOw7KHdYSYNJU4y4pomIiNLJmBHD8NgVxZiYMwLSt/tQ+kQjDhxKrYfwDhYmTURERBRkQvYIPHblfIwZkYHGzztw06pN8HpZ1pFJU4w4PUdEROloxoQxaLi0CMNVSryy+Wvc+concockOyZNMeL0HBERpatTCsbibv2JAIB/vd2GR99ukzkieTFpIiIiol6dP3cSqn56PADg9tVb8d/NX8kckXyYNBEREVGfys8UcOnJx8LnA36/chMa251yhyQLJk1ERETUJ4VCgf/7xQkQfzQBXYe9uPqJRrR+u1fusBKOSVOMuBCciIiGApVSgX9cPA9zpqjh2n8Il6/YgG/3dMkdVkIxaYoRF4ITEdFQMXK4Cv+6rAjHjh2FHc7vcdXjduw/eFjusBKGSRMRERFFbFxWJh67Yj40o4fjoy868bunN+Kwxyt3WAnBpImIiIgGJH/caDxyWREyM5RY9+ku/OU/H8PnS//il0yaiIiIaMC0U3Px94vmQaEAnv5gO5a/Kckd0qBj0kRERERROXfW0fi/n58AALhnbTPWte6ROaLBxaQpRrx7joiIhrLLTs1D2ZkCAOD+d3bh3dbdMkc0eJg0xYh3zxER0VBXde7xWDT7aBz2Atc+tRGffu2WO6RBwaSJiIiIYqJUKnD34tmYPWEE9nYdxhUr7Piq83u5w4o7Jk1EREQUs8xhKvzl7KMx7ajR+KrzAK5YYYf7wCG5w4orJk1EREQUF2MyVXj0siIcNSYTn369B9c+2YSDh9OnhhOTJiIiIoqbSbkjseLyYowarsI7LbtR9fxHaVPDiUkTERERxdWsSTl48NdaqJQKPO/Yifus2+QOKS6YNBEREVHc/eS48bjrglkAgH+sb8EzG7bLHFHsmDTFiHWaiIiIwltaPBXXL5gOALj1xS147dNdMkcUGyZNMWKdJiIiot7dKE6HrnAyPF4fKp52YPMXnXKHFDUmTURERDRoFAoFai6cjTOmj8P+gx5c8ZgdO5z75Q4rKkyaiIiIaFANUynx4K+1+NEx2fhubxcuW7EBHfsOyh3WgDFpIiIiokE3ZsQwrLi8GBNzRkD6dh9Kn2jEgUMeucMaECZNRERElBBH54zAY1fOx5gRGWj8vAN/WPUhvN7UqeHEpImIiIgSZsaEMTBdWohhKgVWb/4Kd73yidwhRYxJExERESXUqQXjcI9+DgDgkbfb8OjbbTJHFJkMuQOQk8ViAQDY7XaUlJRAFMXAdrVaDUmSUFRUBK1WK2eYREREaef8uZPwpesAjK9+ittXb8UxOSPw09nHyB1Wn4bsSJPNZoMkSdDpdCgvL4fBYAAAuFwuSJIEURRRVlaGlStXyhwpERFRerrmxwIuOXkqfD7g9ys3oelzp9wh9SklRpocDgdKS0vR1NQUtF2SJFgsFgiCAEmSUFZWBrVaHdE5RVEMjCz5R5QAQK1Ww2QyAQAEQcDSpUvj90aIiIgoQKFQ4P9+fgK+7jwA2ye7cNXjjXju2lNRcFSW3KGFlfQjTf4pNIfDEdKm1+tRWVkJnU4HnU6H0tLSqK5hMpkCI00AYDQaYbVaYTAYIk7CiIiIaOAyVEo8cPE8zJmcA9f+Q7h8xQZ8u6dL7rDCSvqRJp1OF3a7JElBrwVBgM1mC7y2WCwh+/jPJwhC4HVdXR2qq6sD22w2G5xOJ6xWK2w2G8rLy2G1WuPxVoiIiCiMUcMz8K/Li3Hhg+9iu3M/rnrcjmfLTsao4cmVpiRXNANgs9mg0WiCtmk0GjgcDmi12l6TrZ7nEEURWq0WFosFOp0uaKpOFEUmTERERAkwLisTj185Hxc++A4++qITv3t6I0yXFiJDlTyTYimbNLlcrrDbnc7IFpFJkgS9Xg9BEOByuSCKInQ6HcrKylBXVxcYpeptTVNXVxe6un4YPnS73QAAj8cDjyc5K5x6PB54vd6kjY/kxf5BfWH/oP7Eo49MzR2Bhku1uORfdqz7dBf+/OIW3H7+TCgUijhGGmwg8aZs0tSb3pKpngRBQEdHR9i2ysrKfo+vqanBsmXLQra3trYiKys5F7B5vV44nU60tLRAqUyezJ2SA/sH9YX9g/oTrz4yBkDlGUfhjte+wTP2Hcj07MNFJ+bGL9Ae9u7dG/G+KZs0qdXqkFElp9OZsIXb1dXVuOmmmwKv3W43pkyZgoKCAmRnZyckhoHyeDxoaWnBtGnToFKp5A6Hkgz7B/WF/YP6E88+Mn06oBj9Of768id4zOHEiQWTcf7ciXGKNJh/pigSKZs0iaIYKA3QnX890mDLzMxEZmYm6uvrUV9fHxjeU6lUSf0LRalUJn2MJB/2D+oL+wf1J5595MrTBXzVeQAPv9UGw/ObcXTOSJw6bVwcogw2kFhTaoy1+9Rb9zvggB9qLSW6REBFRQW2bt0Ku92e0OsSERGlu+qf/giLTjwGhzw+lP+7CZ9+Hfmo0GBI+qTJZrMFaijV1NQE6jYBgNlshsFggMVigclkgtlsTnh89fX1mDlzJoqLixN+bSIionSmVCpwr34O5udpsKfrMK5YYcdXnd/LFo/C5/P5ZLt6GnG73cjJyUFnZ2dSr2lqbm7G9OnTObxOIdg/qC/sH9Sfwewjrv0Hsfihd9H67T4cf/QYmK85BWNGDIvLuQfy+Z30I01EREQ0tKlHDcdjV8zHUWMy8enXe3Dtkw4cPOxNeBxMmmLE6TkiIqLBN0UzCisuL8ao4Sq83fIdqp7/CImeLGPSFCMuBCciIkqMWZNyUP9rLVRKBZ537MR91m0JvT6TJiIiIkoZZx03Hnf+chYA4B/rW/DMhu0JuzaTphhxeo6IiCixLpo/FdefPQ0AcOuLW/DGtm8Tcl0mTTHi9BwREVHi3VgyAxdqJ8Hj9aHiKUdCajgxaSIiIqKUo1AoUHPhbJyUr8HersO46rFG7HIfGNRrMmkiIiKilJSZoYLp0kII40Zjp+t7XP1EI/YfPDxo12PSFCOuaSIiIpKPetRwPHp5MXJHDcNHX3TixpWb4PUOTikCJk0x4pomIiIieeWNG42G3xRhuEqJNR9/g9pXPx2U6zBpIiIiopRXnKfB3foTAQANb0p46oPP434NJk1ERESUFs6fOwk3ijMAAH956eO4lyJg0hQjrmkiIiJKHtcvmIYL5/1QiuCzr/fE7dxMmmLENU1ERETJQ6FQoGbxbMz/XymCKx+zY9ee+JQiYNJEREREaSUzQwXTJYXI/18pgtLHG/H9QU/M52XSRERERGknd/RwrPhfKYIP41SKgEkTERERpaXupQhe/fhrGGMsRcCkiYiIiNJW91IEpjclPP3B9qjPxaQpRrx7joiIKLl1L0Xw55e24K3m6EoRMGmKEe+eIyIiSn7dSxH89ikHWnbtHfA5mDQRERFR2vOXIig6Nhd7DhzGVY/b0bHv4IDOwaSJiIiIhoTMDBVMlxZicu5IfL57P655sgkHD3sjPp5JExEREQ0ZY7My8ejlxcjKzMAHbU7c8fLWiI9l0kRERERDyowJY/CPX82DUgE8v3FnxMcxaSIiIqIh56zjxuOWRTMHdAyTJiIiIhqSrjwtD7rCyRHvz6SJiIiIhiSFQoFbFv0o4v2ZNMWIxS2JiIhS1zBV5KkQk6YYsbglERHR0MCkiYiIiCgCTJqIiIiIIsCkiYiIiCgCTJqIiIiIIsCkiYiIiCgCGXIHICeLxQIAsNvtKCkpgSiKAIC6ujoIggBJkqDT6SAIgpxhEhERURIYskmTzWaDJEmorKyEVquFXq9HU1NTIJHS6XRwuVwoLS2F2WyWOVoiIiKSW0okTQ6HA6WlpWhqagraLkkSLBZLYFSorKwMarU6onOKohgYWZIkCUVFRQAAp9OJ3bt3AwDUajUcDkf83ggRERGlrKRPmvxJUbjkxT86BBxJfKIdFTKZTDAajQCAJUuWoLS0FC6XCzabDU6nM7Y3QERERGkh6ZMmnU4XdrskSUGvBUGAzWYLvLZYLCH7+M/XfY1SXV0dqqurA9vUajXMZjMcDge0Wi3XMxERERGAFEiaemOz2aDRaIK2aTSaQLLTW7LV8xyiKEKr1cJisQTWMTU0NKCyshINDQ2orq4erLdAREREKSRlkyaXyxV2e6TTaZIkQa/XQxAEuFwuiKIInU4HtVoNtVodGLXqLfnq6upCV1dX4LXb7QYAeDweeDyeAbyTxPF4PPB6vUkbH8mL/YP6wv5B/UnVPjKQeFM2aepNb8lUT4IgoKOjI2xbWVlZv8fX1NRg2bJlIds//PBDZGVlRRRDonm9XnR0dGDv3r1QKlmii4Kxf1Bf2D+oP6naR/bu3QsA8Pl8/e6bskmTWq0OGVVyOp0R3z0Xq+rqatx0002B1zt37sTMmTNx1llnJeT6REREFD979uxBTk5On/sofJGkVklAoVAEZYH+6bXuZQhyc3PR1taWsMSpO6/XixkzZqCpqQkKhaLf/YuLi2G322Pap6/2cG1utxtTpkzBjh07kJ2d3W+MconkeyP3+aM5x0COYf/oHfsH+0df0rV/DOQ4OfoHkLp9xOfzYc+ePZg4cWK/I2QpNdLkcrkCCVHPu9r8tZbkSJgAQKlUYvjw4f1mqX4qlarfTtXfPn2199WWnZ2d1B06ku+N3OeP5hwDOYb9o3fsH+wffUnX/jGQ4+TsH0Bq9pFIP7uTPmmy2WywWq0AjqwjKi4uDizONpvNMBgMgaxR7srdFRUVcd23v336ah9ILMlmsGOPx/mjOQf7R3ywf7B/9CVd+8dAjmP/6Fss8afM9BzFzu12IycnB52dnUn9VwDJg/2D+sL+Qf0ZCn0kdZa3U8wyMzNx2223ITMzU+5QKAmxf1Bf2D+oP0Ohj3CkiYiIiCgCHGkiuFwuGAwGPpyYwrJYLLBYLDAYDEGPKiICjvQPm80Gg8EQ9tFVRMCRZ8WmAyZNhMbGxoiLgtLQYrPZIEkSdDodysvLYTAY5A6JkojL5YLdbocoiiguLg48+JyoO5vNljZ/lDNpSjMOhwOFhYUh2yVJQl1dHSwWC+rq6oKSJFEUZSvVQIk10P4hiiIqKysD+xQVFSUyXEqwgfYPtVodSJSsVivKy8sTGS4lWDSfL/7/71kmKGX5KG2YzWZfU1OTL9yPVavVBv6/tbXVp9PpgtorKyt9TU1Ngx4jySeW/uHz+Xw6nc7X2to6qDGSfGLpH1ar1VdZWcn+kcai7R9ms9nn8/l8oigOfpAJkPR1mihyvT1cuOc6A0EQuDZlCIqlf9TV1aG6ujp9/lqkELH0D1EUodFoUF5eHqirR+klmv7hcDggiuKgx5ZInJ4bAmw2GzQaTdA2jUaTNnPMFJv++ofNZoMoitBqtbBYLHKESDLqq380NDSgrq4OwJGpOi4EH3oi+f1hsVggSRIaGhrkCDGuONI0BPS2yNv/wOOei/S0Wm0iwqIk0Vf/8D/jURAEuFwuiKLY61+clJ766h9LliyBzWYLPLlB7qcyUOL11T/8f2yl08wGk6YhrPti33QbQqXYuVwuCIKAjo4OuUOhJOR/Fqg/iebvEOqu581Gra2t8gUTR5yeGwLUanVgVMnP6XTyjjkCwP5BfWP/oL4Mtf7BpGkI6O0vQN4+TgD7B/WN/YP6MtT6B5OmNNV9aLTnHU/+ejvp+pcA9Y/9g/rC/kF9Gcr9g2ua0oh/MSYA1NTUoLi4OLDewGw2w2AwoLi4GHa7nQs2hyD2D+oL+wf1hf3jCD6wl4iIiCgCnJ4jIiIiigCTJiIiIqIIMGkiIiIiigCTJiIiIqIIMGkiIiIiigCTJiIiIqIIMGkiIiIiigCTJiKiHiRJSstrEVFsmDQREXVjMBiCHhMx2CRJQkNDQ8KuR0TRY9JERGnFYDCgrq4OdXV1Az62oaEBBQUF0Gq1AI48Y6u8vByFhYVQKBTIzc1FeXk5bDYbgCMJT3l5OQoKCqBQKFBQUIDy8nI4HI6IrymKIlpbWwPnJKLkxceoEFHaqKurw+7du7F06VIsWLAA69atCyRA/fEnQP7na3XncDhQWFiIyspKGI3GkHaLxQK9Xg+TyYSysrKoYi8sLERTU1NUxxJRYnCkiYjShsFgwNKlS6FWq7FkyZKIEyYAMBqNKC8vD9vmf2L72LFj+2zXaDQDire7pUuXwmAwRH08EQ0+Jk1ElBYsFgsAQKvVQhAEmEymAR1vs9kCT22XQ1lZGdc2ESU5Jk1ElBasVitEUYzqWIvFMqBRqcGgVquh0Wi4tokoiWXIHQARUSwsFgusVisaGhqg1WoDC7MrKysjPofVakVJSUncY1MoFBBFESUlJYEpPKPRCKfTiY6OjpD9RVGEw+GIOvkjosHFpImIUppOp4NOp0NDQwOqq6ujmmKTJAl6vb7f/VauXInW1tawx4cjimLQwnKbzQZJksIuNgeAgoIC2O32CKMmokRj0kREKc9/i3+0IzSSJEW0iHvp0qVhR7BsNlvItJrL5QpKxPyvKysre41TrVaz2CVREuOaJiJKeY2NjVCr1YEpsIFyOp1RH9sb/x18fnq9HoIghC1Z4CcIQkILaxLRwHCkiYhSXlNTE4qKimI6x2AkK/5ErK6uDjabLezUXs8Y4p28EVH8cKSJiFJeY2NjTHe/aTQaOJ3OOEb0A4fDAYPBALPZDEEQAByZzguXpDmdzsA+RJR8mDQRUcpzOBwoLi6O+vjBnBbT6/WBxep+Vqs17IiSy+Vi0kSUxJg0EVFKi3UROHCkIGZfd635E6rdu3f32d5ztMq/ENxsNgdt760Wk91uR0FBQSQhE5EMuKaJiFJarIvAAaCkpCTsI0xcLhcMBgMaGxsBHHmgr/8uOFEUIUkSjEZjIAkyGo1oampCeXk5nE4nLBYLdDpd4OHBu3fvDlQuD8fhcODhhx+O+n0Q0eDiA3uJKKX5nxc30Mem9FRQUICmpibZFmI7HA7U1NSEjEoRUfLg9BwRpRyLxRJ4TpvNZuv1QbsDUV5eLuuz32pqauLyPoho8DBpIqKUYzAY0NraCkmSIAhCXJ4bV1lZ2Wul7sEmSRJcLhcfn0KU5Dg9R0Qpxz8i1NTUBKPRGLcpNUmSAuUBEqmkpARms5k1moiSHJMmIqJuHA4HJEmK6hl20airq4NOp2OpAaIUwKSJiIiIKAJc00REREQUASZNRERERBFg0kREREQUASZNRERERBFg0kREREQUASZNRERERBFg0kREREQUASZNRERERBFg0kREREQUgf8PoyM0p/P9oZEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(fgrid, AmplHM)\n",
    "\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlim(min(fgrid), max(fgrid))\n",
    "plt.grid(alpha=.5)\n",
    "plt.xlabel(r'$f\\ (\\rm Hz)$', fontsize=15)\n",
    "plt.ylabel(r'$|\\tilde{h}_+(f)|\\ ({\\rm Hz}^{-1})$', fontsize=15)\n",
    "plt.title(r'$\\texttt{TEOBResumSPA}$', fontsize=22)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FEATURE EXAMPLE: detector relative orientation and distance"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### As an example we verify that the LIGO Hanford and Livingston detectors are rotated by $\\sim90^\\circ$ and compute their great circle distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The angle between L1 and H1 is 89.9°\n",
      "The great circle distance between L1 and H1 is 3027.1 km\n"
     ]
    }
   ],
   "source": [
    "ang = utils.ang_btw_dets_GC(alldetectors['L1'], alldetectors['H1'])\n",
    "dist = utils.dist_btw_dets_GC(alldetectors['L1'], alldetectors['H1'])\n",
    "print('The angle between L1 and H1 is %.1f°'%ang)\n",
    "print('The great circle distance between L1 and H1 is %.1f km'%dist)"
   ]
  }
 ],
 "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.9.4 (default, Apr 16 2021, 21:48:01) \n[Clang 12.0.0 (clang-1200.0.32.29)]"
  },
  "vscode": {
   "interpreter": {
    "hash": "b9e741e012400601e1c0b02fc8f7f31880c5dc826bdb78111d362dc9671d5597"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}