diff --git a/code_new/ringdown_overtones_example.ipynb b/code_new/ringdown_overtones_example.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..b643c00101ad0349c9ce804ada78cae2067a912e --- /dev/null +++ b/code_new/ringdown_overtones_example.ipynb @@ -0,0 +1,407 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy\n", + "from pycbc.waveform import generator\n", + "from matplotlib import pyplot as plt\n", + "from pycbc import psd, filter\n", + "import numpy as np\n", + "from random import seed\n", + "from random import gauss\n", + "import pycbc\n", + "import pylab" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "genclass = generator.TDomainFreqTauRingdownGenerator\n", + "detectors = ['H1', 'L1']\n", + "\n", + "delta_f = 1./8\n", + "delta_t = 1./2048\n", + "\n", + "f_lower = 20.\n", + "f_final = 1024.\n", + "\n", + "params = {}\n", + "params['amp220'] = 5.e-21\n", + "params['amp221'] = 2\n", + "params['f_220'] = 256.607043264397\n", + "params['f_221'] = 250.9439959543337\n", + "params['tau_220'] = 0.004047225085452862\n", + "params['tau_221'] = 0.001338559280517331\n", + "params['dec'] = -1.25\n", + "params['distance'] = 410.0\n", + "#params['f_lower'] = 18.0\n", + "#params['f_ref'] = 20.0\n", + "#params['final_mass'] = 63.0\n", + "#params['final_spin'] = 0.64\n", + "params['inclination'] = 2.5\n", + "params['injtype'] = 'ringdown'\n", + "params['lmns'] = '222'\n", + "params['phi220'] = 2.0\n", + "params['phi221'] = 1.1\n", + "params['polarization'] = 1.75\n", + "params['ra'] = 2.2\n", + "params['t_final'] = 1.\n", + "params['tc'] = 0.\n", + "params['approximant'] = 'TdQNMfromFreqTau'\n", + "\n", + "tc = params['tc']\n", + "\n", + "gen = generator.FDomainDetFrameGenerator(genclass, \n", + " epoch=tc,\n", + " detectors=detectors,\n", + " variable_args=params.keys(),\n", + " f_lower = f_lower,\n", + " f_final = f_final,\n", + " delta_f=delta_f,\n", + " delta_t=delta_t)\n", + "\n", + "signal = gen.generate(**params)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAIbCAYAAABhdZgEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAABcSAAAXEgFnn9JSAABlS0lEQVR4nO3dd3hj1bX38d9WcRuX6ePpjemNYegdQicklCSkEEKAdHLTbtrNTS73Te8BElJuICSBdAJpQOi9zMA0ZobpvY+nuY2LpP3+IVk+R5ZsyZIs2ef7eR49to6OpG0d20dLe621jbVWAAAAAJAPvkIPAAAAAMDARcABAAAAIG8IOAAAAADkDQEHAAAAgLwh4AAAAACQNwQcAAAAAPKGgAMAAABA3hBwAAAAAMgbAg4AAAAAeUPAAQAAACBvCDgAAAAA5A0BBwAAAIC8IeAAAAAAkDcEHAAAAADyxtMBhzFmkTHmC8aYvxpjdhpjrDHGFnpciYwxE40xHzfGPGKM2WuMaTfG1MWuv6Wb+51jjPkfY8y/jDEHYj/f1j4cOgAAADzOWFt076/7jDHmQUlvTdxurTV9P5rUjDHPSzpDUquklyXtlTRF0kmxXX5orf10kvstl7QgYfM2a+2kvA0WAAAAcPB6wPF5SYMkLYldtkoqLcKA4w+Snpf0a2ttg2P75ZIelBSQdLG19tGE+31H0lFFf7adklaLgAMAAAB9yNMBRyJjTIuKMODojjHm55I+KOkea+37u9mvVtIeEXAAAACgD3m6hqM3jDEVxpgvGmOWGWMaY5eXjTHvK9CQVsS+jinQ8wMAAAApBQo9gP7EGDNS0mOS5itaR/GMJCPpdEn3GGNOtNZ+vI+HNSX2dW8fPy8AAADQI2Y4MvMrRYON2yRNstZebq29TNIMSa9KusUYc0lfDcYYM1jS9bGrf+ur5wUAAADSRcCRJmPM8ZIuU7QA+9PW2taO26y1+xSto5Ckj/ThsH4maYSinase6MPnBQAAANJCSlX6Lop9fdBaG0m80Vq7zBjTKOlk53ZjzI8lXZDhc33RWtttABHrsHWtpEOS3mOp/gcAAEARIuBI36TY168bY77ezX5lCdfHKZpylYma7m40xlwn6ZuSmiRdbq3dnOHjAwAAAH2CgCN9Helnz0valO6drLVX5nIQxpg3K1pL0i7pamvty7l8fAAAACCXCDjStzP29UFr7fcLMQBjzDmS/qxoZ6x3Jy70BwAAABQbisbT91js61WFeHJjzAmS/i6pVNLN1tr7CzEOAAAAIBMEHGmy1r6iaNBxhjHmJ8aY6sR9jDEL8tEW1xgzQ9IjkqolfcJae0+unwMAAADIB+Pl5kbGmMslfdmx6WRF05VecWz7qrX2X7H9Ryr6xn+hpCOSlkvarWiR93xJ4yXdZq39ZI7HuUzS8ZIOSHooxW5rrbXfSrjfzZJujl0NSjpBUpukZY7dPmqtXZrL8QIAAAAdvF7DMULSKUm2n5KwjyTJWrvfGHO6pA9IeqeigcfpkvZJ2izpdkl/yMM4hzjG8r4U+zwj6VsJ28ap689XkrCty0wNAAAAkCuenuEAAAAAkF/UcAAAAADIGwIOAAAAAHlDwAEAAAAgbzxVNG6M2SupQtKOQo8FAAAAyKHxkpqttbWFHkgiTxWNG2PqS0tLq6ZOnVrooQAAAAA5s2nTJrW2tjZYa4uuA6mnZjgk7Zg6ders1atXF3ocAAAAQM7MmTNHa9asKcosHmo4AAAAAOQNAQcAAACAvCHgAAAAAJA3BBwAAAAA8oaAAwAAAEDeEHAAAAAAyBsCDgAAAAB5Q8ABAAAAIG8IOAAAAADkDQEHAAAAgLwh4AAAAACQNwQcAAAAAPKGgAMAAABA3hBwAAAAAMgbAg4AAAAAeUPAAQAAACBvCDgAAAAA5A0BBwAgK/vqW/TTpzdp2fbDhR4KAKAIBQo9AABA/3bL75ZqydbDKgv69NinztH4oRWFHhIAoIgwwwEA6LVwxGrp9iOSpJb2iO5+YUthBwQAKDoEHACAXjvY2KpwxMav/2nJDh091l7AEQEAig0BBwCg1/bWt7iuN7WF9ccl2ws0GgBAMSLgAAD02r761i7b7nlhq9rDkQKMBgBQjAg4AAC9ljjDIUm7j7boodf3FGA0AIBiRMABAOi1fUe7BhySdNfzW2StTXobAMBbCDgAAL3mnOE4d8aI+Pcrdx7Vkq2sywEAIOAAAGRhnyPgOG/GSJ0+dVj8+i+f21yIIQEAigwBBwCg15wBx6jqMt181uT49cfe2KetdU2FGBYAoIgQcAAAem3vUWfAUapzp4/UlBGDJEnWSr9iIUAA8DwCDgBArxxrC6u+JRS/XltTJp/P6KYzO2c5/vTqTh1tZiFAAPAyAg4AQK8406mMkUZUlkqSrl44TkMqgpKkY+1h/W4xCwECgJcRcAAAesXZoWp4ZakC/ugppbzEr+tOnRi/7cFlu/p8bACA4lE0AYcxpsIYc6Ux5i5jzDpjTIsxpskYs8IY8xVjTGWhxwgA6OSc4aitLnPddt7MkfHv9zckX6sDAOANRRNwSHq3pAck3SgpLOnvkp6TNFnS/0paYowZmfruAIC+lNihyqm6LBj/vqElxCKAAOBhxRRwtEv6haTZ1trZ1tp3WGsvkTRD0jJJMyX9qIDjAwA47D3aGv++tqbUdVt1WSD+fShi1dIe6bNxAQCKS9EEHNbaX1trP2StfSNh+x5JH4tdvdoYU9L3owMAJHLNcFS5ZziqHDMcktTQQqcqAPCqogk4erAi9rVU0rDudgQA9A1XwFHjDjjKgj4F/SZ+vZ6AAwA8q78EHFNiX9slHSrkQAAAUXu7KRo3xrhmOZzrdQAAvCXQ8y5F4ROxr49Ya1u73VOSMWZ1ipum5m5IAOBd1lrtr+/8d5xYNC5JVWUBHWpqkxQtHAcAeFPRz3AYYy6TdJOisxtfLvBwAACSDje3qy3cWQieOMMhRQOODtRwAIB3FfUMhzFmpqR7JRlJn7XWrujhLpIka+2cFI+3WtLs3I0QALxp79HOdKqyoE/V5V1PJ1Wl7ta4AABvKtoZDmPMWEmPSBoi6QfW2tsKPCQAQEzion/GmC77OIMQZjgAwLuKMuAwxgyV9KikiZJ+Jek/CzsiAICTM+AYmSSdSnK3xq0/xgwHAHhV0QUcxphKSQ8rmvr0V0kfsCxRCwBFpbsOVR2o4QAASEUWcBhjSiX9TdLJkv4t6V3W2nBhRwUASORKqarpeYaDGg4A8K6iCTiMMX5Jv5d0vqTnJF1trW0r7KgAAMk4i8ZHVpUm3afaMcPBOhwA4F3F1KXqFklXxb6vk3RnsiJESf9pra3rs1EBALrY51iDI9UMR7VrhoOUKgDwqmIKOIY4vr8q5V7SrYoGJACAAknsUpVMFTMcAAAVUUqVtfZWa61J47K10GMFAC9rDYV1sKkz4zXZKuNSYg0HMxwA4FVFE3AAAPqHAw2trusjq5PXcLi7VDHDAQBeRcABAMiIM51q6KASlQb8SfdzBhyNrSHR4RwAvImAAwCQkb1HO2c4UqVTSVJ1eWdKVThi1dxGl3MA8CICDgBARpwzHKNSpFNJ7hkOSaqnjgMAPImAAwCQkXQ6VElSacCvkkDnaYY6DgDwJgIOAEBG9rpmOFIHHJJ78T86VQGANxFwAAAy4lxlvKeAw9kal7U4AMCbCDgAABnZ3+BcZTx1DYeUOMNBwAEAXkTAAQBIm7W21zMcpFQBgDcRcAAA0lbfEtKx9s72tt0VjUvuTlX1x5jhAAAvIuAAAKRtv6NgPOg3GlJR0u3+VRSNA4DnEXAAANLm7FA1sqpMPp/pdn93ShUzHADgRQQcAIC0Oes3amu6T6eSpGpqOADA8wg4AABpc3ao6m6V8Q5VdKkCAM8j4AAApC2TDlVSQtE4MxwA4EkEHACAtDlrOHrqUCVRwwEAIOAAAGRgX31mMxws/AcAIOAAAKQt44CjvHOGg5QqAPAmAg4AQFpC4YgOOIrG0+lS5azhaGwNKRKxeRkbAKB4EXAAANJS19gmZ7yQXpeqzhkOa6WmNtKqAMBrCDgAAGlxplNVlQVUURLoZu/O/ZzqqeMAAM8h4AAApCXTDlWSFPT7VBbsPNWw+B8AeA8BBwAgLc4ZjnTqNzpU0xoXADyNgAMAkBZnwDGyKv2Aw73aODMcAOA1BBwAgLQcbGyLfz+iqueC8Q4s/gcA3kbAAQBIS1soEv++POhP+37OGY76Y8xwAIDXEHAAANLSGu4MOIIBk/b9nDUcdKkCAO8h4AAApKXdMcNR4k//9FFd7qzhIOAAAK8h4AAApKXdMcNREkj/9OGu4SClCgC8hoADAJCWNmdKVQYzHFWlzHAAgJcRcAAA0tIesvHvMwo4nEXjzHAAgOcQcAAA0tKWk5QqZjgAwGsIOAAAaWlzFY1n0KWqnBoOAPAyAg4AQFrae1vDUUYNBwB4GQEHACAtve9SRcABAF5GwAEASIszpSqTGQ7nwn+NrSGFI7abvQEAAw0BBwAgLW3h7LtUSVIjsxwA4CkEHACAtDhTqkozSKmqLHUHHLTGBQBvCfS8S98xxiySdKGkk2OXsZJkrU2/HQoAIC96m1IV8Ps0qMSvprawJOo4AMBriirgkPRlSW8t9CAAAF31tmhciq7F0RlwMMMBAF5SbAHHS5JWSloSu2yVVFrIAQEApEjEKhRx1nBkNvFcVRbQ3vro98xwAIC3FFXAYa39tvO6MWRSAUAxcK4yLkklGaRUSe7CcWo4AMBbKBoHAPSoPTHgyDClyr3aODMcAOAlBBwAgB45C8alzIrGpWgNRwdqOADAWwg4AAA9ag+7F+vLPOBgtXEA8KqiquHIFWPM6hQ3Te3TgQDAAJGYUtWbovEO9QQcAOApzHAAAHrU6kipKvH7Mm7qUe1IqaJoHAC8ZUDOcFhr5yTbHpv5mN3HwwGAfs85w5Hp7IYkVZNSBQCexQwHAKBH2Sz6J1E0DgBeRsABAOiRs0tVpgXjEkXjAOBlBBwAgB61hbMNOJjhAACvIuAAAPTI2Ra3tBcpVdXlzHAAgFcVVdG4MeZySV92bCqJbX/Zse2r1tp/9enAAMDjsk+p6pzhaG4Lqz0c6dXjAAD6n6IKOCSNkHRKku2nJOwDAOhDri5Vgcy7VDlrOCSpsSWkIYNKsh4XAKD4FdXHS9bae6y1pofLPYUeJwB4jatLVS9mJipLAnIu3UFaFQB4R1EFHACA4tSaZUqVz2dUWeJcbZzCcQDwCgIOAECPsl2HQ5Kqy52dqpjhAACvIOAAAPSoPZRdSpXkruNghgMAvIOAAwDQo2zX4ZBY/A8AvIqAAwDQI+c6HMFeplSx+B8AeBMBBwCgR205TqlihgMAvIOAAwDQozZX0Xjm63BIUjUzHADgSQQcAIAe5bponBkOAPAOAg4AQI9yUzTeOcNBlyoA8A4CDgDwgJ2Hm3X93Yv10fte69Wbfec6HL0vGmeGAwC8KNDzLgCA/u7/nt2sZ9cfkCSdNnW43nvqxIzu3xbq7FKVm3U4CDgAwCuY4QAAD9h4oDH+/YH6lozv35bzlcYLn1IViVhX9y0AQH4wwwEAHrDnaGeQEYrYbvZMzlk0HvT3tktV4VOqrLVasfOofvfKNv1jxR6VBHz6ybtP0JnThhdkPADgBQQcADDAWWu1N9uAI5yLLlWOovFjfTvD0dga0t+W79J9L2/Xmj318e3H2sO65fdL9a//OEtjB5f36ZgAwCsIOABggGtoDam5LRy/7gwe0tWW46Lx1lBEbaFIr9OzMvGTpzbqzqc2qsnxGjgdaW7Xx+5bqj996LQ+GQ8AeA3/WQFggNt31F2zEe7FDEdbKLdtcaW+qePYUtek7/57XZdg4+RJQ3XTmZPj15fvOKJvPvxG3scDAF7EDAcADHB7EgKO9nB2KVWlvZwFGFTil89IHfFOQ0tIwypLe/VY6XrDkT5VFvTpnSdN0HtOmaBpo6okSU2tIf1hyQ5J0q9e2KqTJw3VpfNG53VMAOA1zHAAwAC3tz5xhiPLlKpeznAYY1yzHH1ROL6lrin+/VnTRujWt8yJBxuSdOtb5mjW6Or49c/9ZaW2Ou4DAMgeAQcADHB7E2Y4Qr2Z4XCsw9HbgENKXIujb1KqOkwePqjL7WVBv+58zwmqLI2Oq6E1pI/ct1Qt7cnrPQAAmSPgAIABLnGGoz3bLlVZFFa7ZzgKH3B0bP/O2+bHr7+xp163/n113scGAF5BwAEAA1ziDEdvUqpac7AOh9T3q41vTSPgkKTL5o3WDadPil//w5IdWrO7PuX+AID0EXAAwACXGHBkWzTe23U4JKm6D2s4jja362BTW/x6dwGHJP3XZbM0bWRl/PrLmw/mbWwA4CUEHAAwwHUtGi9cSpV7tfH8plRtOdg5u1FR4tfIqu47YpUEfDp7+oj49ZU7j+RraADgKQQcADCAtbSHdcjxKb/Uy4X/crAOhyRVl3fOcBxpzm/A4UynmjRskIzpORVs/ria+Pcrdx7Ny7gAwGsIOABgANtf39plW+9mOHLTpWqEY5bhQEPXseXSZmf9xoju06k6LBg32HX/o8fyX9gOAAMdAQcADGCJ6VRS5m1xrbWudTh6u/Cf5A449jd0HVsuOWc4pvRQv9Fh4rAK1ThmYVbtYpYDALJFwAEAA9ieo8e6bGvPsEtVYpF5NjMcI10BR35nOLYkpFSlwxjjSqtaQR0HAGSNgAMABrB9SWY4Mk2pSqz5yKYt7siqsvj3++tbZW3m6V3psNa61+BIM6VKSqjj2MEMBwBki4ADAAawPUe7BhyZtsV1FoxL2XWpGlndOcNxrD2sxtb8tMata2xzPXa6KVWSNN9Rx0GnKgDIHgEHAAxgyWc4Mk2pSpzh6P2pY2hFiQK+zhmSfKVVOWc3BlcENbiiJO37OgvHdx9tyXtxOwAMdAQcADCAJZvhyLRovC0h4Mhm4T+fz2h4paOOI0kXrVzYUtcY/76nBf8S1daUuWpNXt91JFfDAgBPIuAAgAFsX7KUqgxnOJwpVQGfkc/X+xoOyZ1Wla9OVVvqmuPfZxpwSO46jhUFqOMIR6yeXX9Ai7ccUqgX66YAQDEJ9LwLAKAYbT7QqL8t363zZo7U8eMHd7k9HLHalyQdKJzhDEeu1uDoMLIP1uJwzXCk2aHKaf64wXr8jf2SClPH8Y2H3tBdz2+RJNWUB3X+zJG6YNYonT19uKrKgj3cGwCKCwEHAPRDkYjVzb9+VZvrmvSrF7bo+S+cr+qEN6IHG1uTdqRqz6JLVTYF4x1GODtV5Sng2Oqc4cigQ1WHxBXHrbVprVSeC3WNrfrtS9vi148ea9cDy3bpgWW7FPQbnTplmC6dO1qXzK3V0EHp16YAQKGQUgUA/dCaPfXxlbTrW0JJ27cmq9+QMm+L2+pIqcr1DMf+JEXt2YpErLYcdLTE7VVK1eD49web2rTrSNf1TPLlj0t2dKmb6dAetnpuQ53+64HXddLXH9f1dy/Wn17dwYroAIoaAQcA9EPPbjjgur5qd9eAI9kq41LXrlM9cc1wZLEGRwd3DUfuZzh2Hz3mqjtJd9E/p6GDSjR+aHn8+sqdfVPHEQpHdN/LnbMbHzpnin76nhN09cKxGlzhnsHqqPP43F9W6sSvPaYP/fbVpAs9AkChkVIFAP3Qs+vdAcfq3fVd9tnrmOEYVV2qfbGOUJl2qcp1SpVz8b9kbXuz5UynGlVdqkGlvTvVzR83WDsORd/Ar9h5RJfNG52T8XXnibX7tTt23AI+o5vOmKyR1WW6dN5ohcIRvbrtsB5ZtVf/en2Pq/6lPWz179X71BqK6J73n5z3cQJAJpjhAIAiEYlYffffa3Xzr5do/b6GlPs1tYb02rbDrm2re5jhGDekIv59pilVbflMqcrDDEc2LXGdFhRgxXFn7cbFc2s1srozOAv4fTp1yjDd+pY5evmLb9LvP3Cq3nPKBFcdx9PrDmjpdvfvBgAUGgEHABSJB5fv0k+e2qTH39ivz/1lZcr9Xtp0sMtq4VvqmtSUsGq3c4Zj3JDO9KBM2+I6ZzhyEnA4UqoaWkJqaQ9n/ZhOm+uyq9/o4KzjWLXrqCIZBmqZ2nSgUc9vrItfv/7UiSn39fuMTps6TF+/ap4W/9ebtMDRpey2xzfkc5jd2t/QolW7jtLKF4BL0QUcxphyY8z/M8asN8a0GGN2G2PuNsaMLfTYACBT++tb9Jk/rdB3Hlnb7Zswa63+77kt8evLdxzRpgONSfdNrN+I3l96Y487rSpVwGGtMnrz3OYIbnKRUjW8slTOhk+5Xvxva44Cjrlja+LjbGgNuQrR88E5uzFjVJVOnjw0rfsF/D598oJp8evPrD/QZQasL7y4qU7nfvdpvfmO53XKN57Qlx9cpVc2H8x7oAag+BVVDYcxpkzSk5JOlbRH0t8kTZL0fklvNsacaq3dXLgRAoC0YscR/XXpTl2xYIxOnJT6TWEkYvWhe1/Tsu1HJEnlQb8+/qZpSfd9adPBLgHDg8t26TMXzeiy73Mb6rpsk6J1HM7xpEqpkqKzHKU+f8qxOzlTqrJZZbxD0O/TsEElqmtskxT9VHzCsIoe7pW+LY6AozcF4x0qSwOaOqJSG/dHA7+VO49o6ojKrMeXTFNrSPe/tjN+/b2nTcyoDe+500fo+PGDtXzHEUnSbU9s0G9u7LtajlW7juqDv3lNzW3R2aqDTW367cvb9NuXt6m2ukyXzx+tGaOqVN/SroaWUPxreziiK+aP0QWzR/XZWAH0vWKb4fhvRYONlyRNt9Zea609RdJnJI2QdHchBweg8NJNv2kNhfXIqj1ak6SYOlFdY6ved/diXfiDZ/TSpoPd7rtmd72u/cVL+vVL2/TuX76idXtT11r8YcmOeLAhSb94drMON7Ul3feXz2/psu2BZbtkrfvT4R2Hml1vqC90vFFz1nFYa1POcEiZFY67UqoCuVmLIl9rcbSHI9pxuLNT05RerMHh1Fcrjj+4fJcaYilxVaUBXbUws0l9Y4xrluPZPpzl2FLXpPfdvViNCSl9HfbWt+iu57foc/ev1Nf+9YZue2KDfvXCVv3ltZ362/Lduvk3r+peR2euvnawsVX3v7ZTS7cfZjYGyJOiCTiMMSWSbold/Zi1Np5LYK39gaSVks4xxiwqxPhQfI61hbu8GRtorLVpv8G21mrn4ea099+4v1Grdx9N6zXcuL9Bf3p1R1prJmw+0KhvPbxWf3ltZ4/FyTsPN+vTf1quG+9ZEv9kNpX99S26/u7FmvM//9bHfrdUDS2p1x041NSma3/+sj5871Jddvtz8RWbkzna3K733rVYz6w/oA37G3Xzr5ekDFKONrfrw/e+ppb26BvwtlBEn/zjcrWGur7mdY2t+vYja13bGlpD+tkzm7rsu3F/o55cu7/L9p2Hj3V50/iMozvVcSMrEwKOznHXHwvpmON3IXGGI5TBGyt3W9zcnDbytRbHjkPN8d87n5HGD81u5mSBo44jXyuOW2td6VTXLBrXq85a50wfoYUTBsev/+jx9bkYXrf21bfovXe9ooOxQNrvM7rjXQv1nbfN11nThsvvSy9A/e8HV+n/nu3bBIb2cER3Pb9F5373aX3mzyt09Z0v6uRvPKEv3L9ST7yxL/6/1FqrHYea9ciqvfrBY+v1kXtf0//8bZU2p0h57Cu7jxzTlrqmAX8exMBQTClVZ0iqkbTJWrssye1/kTRf0hWSXuvLgaE4femB17VsxxFdOrdWl80brTljqpOmIBxsbNW/V+/TM+v3K+j36aI5tbpg1khVlHT99W8NhfXU2gP6x8rdOlDfqlOmDNVVC8dqSpI0Cmutlu04oj+/ulPLth/WuCEVeseJ43TezJFJC2t3Hm7WHxbv0F+X7pQxRm9bNE7vPmWCRjm60HRobA3pL6/u0K9f2qYtdU06e/oIffCsKTrjuGFdfsZIxOrRNXv106c3acXOoxpSEdRNZ07W9adP6rLytBR90/TDx9brqXXRN66nTRmmz10yQwsnDOmy7776Fn3/0XX682s7Za1UUeLXLecfp5vOnKzSgDsdp6U9rDuf2qifPbM5vmjZfa9s03euma9po6q6vHZ/eW2n/vcfa+Kfij67/oD+8+IZ+uBZU+RLeJPyyuaDuuX3y+JtQP+1co827mvUL993Ypc3lDsONev6uxe7ZgG++s81amkP62PnHdfldX7frxa7Upma2sK68Z4leuBjp2t0TeesQCRi9ak/Ldf2Q82ux3hjT71++NgGfeHSma7t33xobdLF2O55cavef8Zk1dZ0Hve7X+gMiGaMqtKwyhK9GJtp+euyXa40KWc73LOmDdecMdXx6+v3NagtFFFJwOdKpyrx+zTC8QZfUkZFvbnuUiXlr1OV87iPHVLe5fc0U84ZjtW769UejuTsNeiwZOthrXXMlF3XTbF4d4wx+tQF03X93YslRVPvXtt2SIsmplcLkqmjze26/q7F2umYUfre2+frigVjJEnvOHG8Dja26uFVe/XEG/vU3BZWVVlQ1WUBVZUFVFUW1P1Ld8YXqPz6Q2+oqS2kT7xpWt5XdX9+Q51u/cfqeLpch7rGVv1hyQ79YckOlQf9ml5bpc0HGtXQ0nX25t5Xtuvti8bpExdMc/2vyEY4YuUzSvnzH25q0z9W7tb9r+3UitjaMCOqSnXmccN1+tRhOnPacI2uKVckYrXxQKOWbz+iZTuOaPmOIzrU1KqTJg3VFQvG6JzpI1QWTP63EY5Y7a1v0dCKEpWXdP/3E4lYrd3boOa2kKbXViU950jR88PiLYf0/MY6ba1r0pQRlTpx4hAtmjhEQxxd1jqef8P+Bq3ccVTr9zVoyKASnTZ1mOaNrUn5t9caCmvd3ga1hiKaNbpald0E7JGI1ea6Rm2ta9a4oeWaOqIy5eM2t4W0cudRrd5dr4oSv+aNrdH0UVVJa9kiEautB5u0Zk+9mtvCmlVbrRm1yfft+Dm3HmzSzsPHNHZwuSYPH9RtgF7X2Kr1extUEvBp2qgq1ZQnf62LVTEFHAtiX5emuL1j+/w+GAuKXFsoosfe2KeGlpDufHqT7nx6k8YNKddl80brkrm1mjRskB5dHe1V/+Kmg65P2v+5co/Kg369adZIvWXBGJ09fYSWbT+ivy3fpYde36N6x4ll8dZDuuPJjVo4YbCuPmGcrpg/Wq2hiP66dJf+8toObTrQ+eZm7d4GPf7GPo2oKtXbFo3TtSeO1/ihFXpq7X7d98o2Pb3+gJwfRN32xAb95KmNunhura4/daJOnjxUOw4d0z0vbtWfX90RT6+Qom8yn11/QLNGV+sDZ03Wm+dHT+gPLtulnz27SZsd4zjc3K7vPbpeP392s95/+iTdeOZkDa4o0erdR/XDxzbo8Tf2uV7LlzYf1FV3vqiL54zSZy+eoeNGVqmpNaRfPLtZv3h2s+tT8ua2sL7zyDr9cckOffny2XrTrJEyxuiZ9Qf0lb+t0raD7jfjy7Yf0eW3P69bzj9OHz5nqkoCPtU1tuqLf31dj61xjyMUsfrWw2v1wsY6ff8dCzSyqkzWWv3yuS361iNru8yWrNvXoLf8+Hn97LpFOmXKMEnRPPIbfrVEdY1d38B+99/r1BqK6FMXRN/MHGsL66YUMyt761v0/l8t0Z8/fJqqYifQO57c6JqFmDJ8ULwb0s+f3aTzZ46MF/m+svmg7l/amY//2Ytn6J4Xt+pAQ6taQxHd9sQGffPqeZKiszHO3P2bzposnzHxgONfK/fof66YrdKAX+3hSHy7JJ09fYSmjaxS0G/UHrZqD1ut39eguWNrXAvA1daUKZiwYF8mrXHbXClVOQo48rT43xZXwXj29RazRlcr4DMKRaxaQxGt39egOWNqer5jBn7z0tb492ccN0zHjez9uM+aNlwnTBispbFUvh89vkG/vemULEfY1bG2sG789RKtc7Rv/vKbZ+uqheNc+w2rLNV1p05MGURde9J4veeXr8QD+R89vkFNrSH912Wzsgo66lvatX5vg/w+o5KAT6UBn0oDfjW3hfXDx9brkdV7XftXlPjVGoq4/i6OtYe1opuZ13DE6g9Lduivy3bp+lMn6qPnHachFUHVNbZpw74Grd/XoPX7G9XUGtIpk4fpojmjNLyytMvjtIbCemzNPv1xyQ69uOmgKoJ+HTeqUtNGVuq4kZWaNrJK7eGIHli2S4+/sa9Lh7oDDa16YNkuPbBslyRpwtAKHWpqS5ri9s+Ve/TPlXtUVRrQhXNG6YoFYzRp2CC9vuuoXt95RCt3HtWqXUfV1BZW0G+0cPwQnX7cMJ0+dbiOHz9YJQGfDja26rkNdXomdl466EgTnTSsQnPG1mjOmGrNGFWljfsb9dyGOi3eesj1wYXU+f9/yohBOnHiEFWVBbVy5xGt2lXvOu90GFTi10mTh+q0KcM0d2yNNtc1adXOo3p9VzQwCTlmNqePqtLCCUO0cMJgLRg3WAebWrV022G9tu2wlm4/4vowqCTg04xRVZozplqzx1SrLOjX8h1HtGz7Ea3bW6/Ef5Ulfp9mjq7S3LE1mjqiUtsONmn17nq9EQs0XPsGfJo1uloLxtVo3tgatYQiWhPbd93eBtfPWR70a+boKs0eXa05Y2o0uCKoN/bUa/Xueq3efTS+jlKHMTVlmlFbpRm11Zo0rELbDzVr95HiXfjTFMtUnDHmB5I+JemH1tpPJ7l9gaTlkpZaa7tNqzLGrE5x09TZs2eXrl6d6mb0F0u2HtLbf/ZSTh7L7zNpv/kK+qP7pvtebUhFUIebU6f+OI0bUq5dR44pnT/J2tisSKqVpJ0Glfh1/ITBemFj97UJUvQf9aXzRmvJlkNd3gQaoy5jO2f6CFWWBvSv1/e4tleXBVyBmyTNrK3SO08ar9uf3KhDjhNUVVlAg0oCrp9leGWJvvrWufrHyt166PXONwZBv9EFs0bp4VXubV+7cq7GDC7Xh3/7mpoc//A/fM5U/Xv1Xteb0A+dPUWfvmi6PvTb1/T0us7ZgveeOlHTa6v05QdXxbedPX2E7nrfiXp+Y51uvGdJ/Oe/YNZI3fbOhbrijufjQcfYweV65JNnqTTg1+W3P6cNsU9O54yp1t8+doZ+v2RH/LH9PqPHPnW2poyo1B1PbND3H1sf+7lL9cIXzlNbKKKTvv54PHXr5+9dpIvn1Lp+70sCPq34ykUqL4k+X0c61Xeuma93nDRef1yyXZ+//3VJ0smThur3HzxVU//rofjP9sIXztfYwel9Knvb4xv0w1h6ztULx+oH1x6f1v268+sXt+p//h79X3z29BE5K3D+0gOv675XtkuSbjh9km59y5ysH/PNdzynVbuir+83r56nd508IevH7LC/vkWnf+vJ+Juln123SJfMrc3qMZ/bcEDvvWtx/PpfPnxat80NMnH0WLv+tnyXfvvStvjvuCR97Lyp+uzFM7u5Z2r76lt03S9fcT3eO08ar4vn1CocsQpFrCLWKhyxmjJiULcBX1soortf2KLbn9jQ5c1fMsZI7zxpgj578Qz5jPTUuv16bM0+PbPugOt/iSSNH1qu2aOrNXVEpf71+p4uH7AMKvGrNOh3/X9z8hnp5MlDdenc0bp4Tq2OHGvTH5fs0IPLdqV9niik8qBf44eWa8P+xrTOU+h7u3/5UbUf3L7GWpv9P74cK6YZjo6PdJpT3N7xjqEqxe3wkJMmDdULXzhfj6zaq0dW7dGr2w53+w+wtrpMl86rVVNrSA+v2uuaHk8MNsqCPl00u1bTR1Xqnyv3uFIdEj9ZkqRTpwzVZfNGa/GWQ3p09T7Xp8GJJ5FJwyr07lMmqD1sdd/L2+IrCktypSVI0Tfi1544XidNHqrfL97uenOcGGgYI102b7SuP3WiXthYp1+9uDX+Mza1hbsEG/PG1ujTF05XwG/0nUfW6fVd0an5iI1+ou40bFCJPnnhdF02t1Z3PLlRv315W/w1eyZhtWu/z+j9p0/SJy+crsVbDupLD6yKp0us3dugW/+xxrX/WdOG6ztvm6/yoF+fv3+l/r06+qlXXWObPnKfe7JzTE2Z7rxukY4fP1gPLtulz92/Um2hiNrDVp+//3X5jOKBoN9n9K2r5+ntJ47XjWdO0nW/fEXr90XfzPz82c16eNVeV2rUNSeM0/++ZY58PqMdh5r1i1gu+bPrD+hTf1yu5zbUxX+/Jg2r0PffcbwGlQb0w2uP1zU/fVGhiNWuI8d069/XaOrIQfE3TsZIX79qngJ+n649cbz+79nN2h6rMfjBY+v1/Xcs0K8dufvXnzZRpQG/SgN+XTi7Vv9YsVuS9MDSXbp4Tq0rnerkSUPj6Q5zx9TEA45o4fj4+OsuRWc4/D7jChozSqkKd77xyktKVQ5rONwdqnLT+Wr+uMHxgGPlziM5DTjueHJjPNgYU1OmC2aNzPoxzzxuuBZNHBKv//nh4+t1702n9HrGwFqr17Yd1u8X79C/Xt8dD4Q7vOvkCfrPJN3U0jWqukx//NBpeu9dr8R/jzvSmpJZOGGwbj5zii6eM0oBx+/jixvr9JW/d02TSmXRxCH637fM0dyxnQHMVQvH6aqF49QaCuvlzYe0+8gxTRk+SDNHV7tSWD514XT96dUduu3xDfEPZ5rawl2CFKeIlV7efEgvbz4UD7Z7Y3hlqa5aOEbXLBqnUVVlemnzQT2/sU4vbKxzBUEdn/QfP36wjh8/WNXlQT2yaq8eW7Mv6QxCT461h+P/R518RqooCaRsGtBhUIlfp00dpjljarR2b71e23Y43qkukd9nNH1U9NP+HYebtXz7Ede5NZnSgE8lfp8rOyAVY6IfEu092tJjPVtlaUBzx1arqTWattXdOIyJtuKuLA1o7Z7u9+3Yf1RVmfY3tKT1QebwyhK1hiJJU/yKXTEFHDmTKrKLzXzM7uPhIE/GDi7XTWdO1k1nTta++hb9e/VePfz6Xr2y5aAiNvqG5rJ5o/Xm+aN1woQh8bqAr145V8+sO6B/rNyjx2P/eH1GOmvaCF25cIwuml0bL9i85fxpWrO7Xvcv3am/Ld8V/+c4dnC5rlk0Tm87YVy8nef1p03SoaY2/XXpTv1xyY74m86Az+jiObV69ykTdNqUYfFxfOjsKXpi7X799qVtrsW+pgwfpBvOmKRrTugsHL14Tq3W7W3QL5/brL8t3x3/J1bi9+maReP0wbOnxNcbOGXKMN101hT95sWtuuuFLTriCHpm1lbp0xdO14WzR8XffJwxdbgeXrVX33t0nevNWmnAp5vOnKyPnDs1nlZ061vm6F0nT9D//mO1K7VHkk6YMFhfu3KeZsdqCs6fOUqPfmqovv3IWt378nbXvmVBn7502Sxdd2pn68+fXbdI972yXV/95xq1htz/pM+aNly3vXNhfEXlKxeO1YRhFfrgb16Lp091/LMuD/p153Un6LwZ0TduI6vK9IcPnqbrfvmK1sRqNZzBxmXzavXta+bFj8sXLpmpHYea47Mo/3QEYOVBv3723kXxNx4Lxg/Wf7xpmn4Qm6G4f+lOV+rSu0+eoONjC7KVBHz6zEXT9Yk/LI8/7vDK0vj4SwM+veeUzjeyVy8cGw84nly7X0eb27vUb3SYM7ZaejX6/arYG7Z99e6AQ5KCPl/8dyezovHOfXPVpcqZUnUghylVrjU4ctTCdsG4Gv3ulej3uexUtXzHEd37SmfAecMZk1xvoHuro5bjuruig35h40G97Wcv6T/eNE1nTxueUeCxatdR/eefV7g+eHF6x4nj9LUr52ZdczF0UIl+94FTdeM9S3rsrrVs+xF97HdLNXZwuW44fZLOmzlStz2xIf730qE86FdFiV9toYhaQ5H47/74oeX69IXTdeXxY1OOuzTg1znTR6QcQ9Dv03tOmairF47Tr1/aqp8+vcmVpjNsUImmjarUjFj92qNr9rk+BEhUFvTpsnmj9fZF41USMNqwr1Eb9jdqY+zS0NKuM6cN19sWjdPZ00a4fk8umzdal80bLSlaw7Zy51ENHVSi+eNqujQfuGzeaDW3hfTk2v3654o9enLdfrWFIpoyfJDmxdJ+5o8brBmx2pUXNx3Ui5vq9OrWw67/y6OqS3XO9BE6Z/pInXnccFWVBbTjcLNW767Xql1HtWp3vTbtb9SIqlKdPW24zpwWbWjg/MDCWqvth5r16tbDWrr9sFraI5ozploLxtdo9ugaV/3Isbawlm4/rJc2HdRLmw9qS12TJgyt0LyxNfFxHzeyUn5jtOlAo5ZtP6JlOw5H06L2NaiyJKCFE4do0YRo3ciC8TWqKguqNRTWhn3RBiqrd9drze56tYTCmjumRsePH6yFE4ZEHzd2fmiLpVW+viuayrX9YLPGDSmPpWPVaGZtVfw1bwtFtG5vg1bsPKKVO49ozZ56lQb8mjW6SrNGV2vW6GrNrK1SRUlAx9rCWrs3mj61JpZG1djSrumxVK85Y6JpaiOro6nGe462aN3eBq3d26B1e+u14/Axja4p0/2VJdrbczJDQQzIlKpunmP17NmzZ5NSNbAdamrTkeY2TRzWfQGWFC0I27CvUWOHlCfNrXUKhSNasfOI/D6f5o+t6VLY7GSt1YqdR7WvvkULJwzWyKquheFO0S5F+zR9VJXOnjai28feX9+iP7+2U5GI1bUnjdfIJEXnHRpbQ/r9K9v1+q6jumRurS6ZU5vysdvDEf3p1R16YOkuTRtVqVvOn5Yy5cZaq3+vjnZsaWoN65bzj9O1J45P+divbD6oL/71dW2ua9KiiUP0vbcvSLkg27q9Dfr475fGP0n7+PnH6ZMXTE96LHcfOaYP/ObV+KeiwwaV6O4bTnKtutzhaHO7rv/VYldO9nkzRujn7z2xS1FfS3tY7/q/l10tbSXptncer7ce725XGgpH9Pafv9Rl3+GVJXri0+eqpqLzU9FIxOqy259L+ubtXSdPiNd1SNHjceo3nojnR3/24hn63qPr4jMUD3/iLM0aHQ3uXtt2SNf8NJpqVVHi1+u3Xqybfr0kPiv2lTfP1o1nTtasLz8S/2Tz3588WzNq05swvvXvq3XPi1slSTeeMVlfuSL7z212HGrWWd95Kn59w9cvzXr25FhbWLO+8kj8+nOfOy/rLlWStHZvvS750XPx64988izNrK3u5h49C4UjeutPXoj/7h43slIP/cdZOVlYUYr+jV7785e1eOsh1/YF4wfrP84/TufPHNljkLB+X4Pe8fOXXB9aSNFU0bctGqdrT5qQVb1JMk2tIX3/0fV6Zn20Xirg88nnMwr4jFpDyT9hT2SM9J5TorMugys6C5KttWoLR1Ti9+W8KP3osXa9tu2QyoMBTR9VqWEJ55NIxGrFziN6eNVePbxqj3Ycis5oLxhXo3ecNF5XLBiTsuA6n6KzxJEeu6K1tEff8O86fExzx0bfWOe7sD9X2sMR+Y3p9rw6UMyZM0dr1qwhpaoHHR+Bjktxe8f2wjXrRr8wdFBJ/JPwnlSUBJK+OU0m4Pel3e3FGBP/ZDsdx8WKA9MxsrqsS7elVCpLA/rA2VPS2rfjE7v3nNJzhxxjjC6ZO1qXzB2d1mOfMmWYHv/0OdrX0KLa6rJuT1Qzaqv091vO1NPr9mvisEHxN9XJjBlcrj9/+DT9/JnN2t/Qqo+cMzXlAnI1FUHde9PJ+si9S/X8xjqdN2OEfnrdoqRv8MqCfv3y+hN11Z0vxmdDbjh9UpdgQ4r+XvzwHcfrstufc+WM/9dls1zBhiT5fEafu2SGbrzn1S6Pc9OZk1zXg36frlgwJv5G//YnNsSDjZFVpZrpCBZm1lbH06Wa28LaerDJtQZHxwxHwG+k2HvHUCSTlKp8rMPhfkNW19iadaefrY6VwEv8Po1Js0alJzNGRT+R7Ohm9tOnN+m2dy7M6jF/+/I2Vxvjr105N2fBhhT9G/3xexbqU39c7kqpXLHjiG769auaM6ZaX7x0ls50zJQ57TzcrOvvWuwKNk6fOkzvOnmCLpozKuvuX6kMKg3EAtrkQe3rO4/qruc3658r9ySdpVswfrC++tY5mu9oZ9zBGJO3cdeUB3X+zNSLF/p8JlbEPERfvHSmttQ1Kej35SQgzkZJwJfW711Z0K/Tpyb/XSl2ue4qh94ppqOwIvb1hBS3d2xf2QdjAZBjPp/R6JrytD4VKwv6dcnc0d0GGx0qSgL61IXT9c2r5/W4WnVVWVC/velkLf7Sm3T3DSelbAspRbvr/O4Dp+jti8bpkxdM05cun5Vy30nDB+nLb+58g3TalGEpF247b8ZInTTJ3YL43BkjdNzIrrMNzsdwpjOcNW2E63UcVBrQFMeM0erd9a46n3jA4fiEL6OF/3K80rgUPcbVZZ2fee2vzz6typlONWFYRdprQPTEGKOPnDs1fv0fK3Zr+8FU5YY923u0Rd9/tHONjGtOGKdTY53WcmlkVZnuu/lU3f+R07qkB63eXa/r7npFt/59dZe1ew42tur6uxbHf4eMke58zwn63QdO1RULxuTtTXs65o2r0Y/euVDPff48feTcqfHfocEVQX3z6nl64COnJw02iokxRlNGVBY82AD6UjHNcLwg6aikqcaY4621yxNuf1vs6z/6dFQABhRjTI8pbh3GDanQd9++oOcdFU2JqiwNaPOBJt145qSUgZUxRp+7ZKary9rNZyafhZo/rsbVfrfD2dO7ftI4Z0xNvE3za1sPuT6Z7uhq5sz77u0MR64CDik6W1ffEk2RyUVr3M2ulrjZrTCe6LK5tfr+sAptO9isiJV+8dwmfe3KeT3fMYmv/rNz/Zma8qD+67LedXdK16KJQ/XrG0/W8h1HdMcTG/SEo73zPS9u1fMb6/Sja4/X3LE1amwN6YZfLXG9ll+/cl68RqBYjK4p1+cvmamPn3+cXt95VLPHVMdrzQAUn6KZ4bDWtkn6cezqT4wx8bOFMebTiq6/8Yy1lkX/ABSlKxaM0ScumNbjG5+TJg3VR8+dqqDf6NoTx+uM45J/um2M6TJTYky0C1Ei5wKAT67rfEPpM53pS8HeznDkYR0OKXHxv+w7VW3NY8AR8Pv0obM7Zzn+9OrOXo356XX7XW2kv3jpzC75/vly/PjBuuuGk/T3W85wpeRt3N+oK3/ygn7y1EZ98DevxrvWSdJ/XjRd7z4ld125cq2iJKBTpgwj2ACKXNEEHDFfk/SKpNMlbTDG/NEY87Kk70s6IOnGQg4OAHLlc5fM1IavX6Zvv21+t2lmVyYEHHPH1CR9g+pcm6CjIFWKttDsyGH2OzpoZdKlqi3k6FKVyxkOV2vc7Gc4tuQx4JCkq08YGw/e2kIR3f381ozu39Ie1lf+1tm05IQJg/WOE8fncohpmT9usB782Bn6wFmT1fGrF4pYffff61zd595/xqS068UAoDtFFXBYa1sknSfpq4qux3GlpImS7pF0grV2c8EGBwAFMH5ohavmI1k6leSe4XAaXdOZPhb0OVOqerfSeC4Lm50d1rJNqbLWauOBzg5G+Qg4yoJ+3Xzm5Pj1e1/e5mqF2pOfPLUx3oTA7zP6+lXzCtY5pyzo15cun637bj5FY2q6phhetXCsvnz57H7TiQhAcSuqgEOSrLXHrLVfsdYeZ60ttdaOtta+31q7s9BjA4BC+MKlMzWkIqgpwwfpfadPSrrPkEElSdsYj3K8qfe7UqrSr+FwF43n7g2oc4bjQJYpVfsbWl11Kx3rH+Tae06dGC9UbmwN6d6X02ucuK++RT9/pvMzs5vOnJxWU4R8O33qcD38ybN15fFj4tvOmzFC33nbfE+0EQXQN4qpaBwAkMSiiUP16n9fKJ9Rt584zx5TrV1H3CvWO2c4nEXj7b2s4cjlDIezNe6+LFOqOlrWStFAZkiarbEzVVka0PtOn6Q7ntwoSbr7+S268YzJrkXKkrnnxa3xmaLa6jJ94k3T8jK+3qgpD+pH71yo9542SQcaWnTh7NqcdfgCAKkIZzgAAF35fabH9JZkaVWjnClVjtmJcC9TqnJbw+FMqcpuhmOdYzHFmXmeObjh9EkqC0Zfh4NNbfrTqzu63b+pNaT7HDMhN505uceF1gph0cQhumTuaIINADlHwAEAA4SzcLyDc4bDlVKVSVvcUJ4CjurOGY66xraMgqBEztXbZ6a5gnpvDass1TtP6uzc9ItnN7tmgRLdv3Sn6luibXArSwO69uS+LxQHgEIi4ACAAWLu2K6f7NdWd9Z1uIrGiyClylnDEY5YHWpq6/Vj9WXAIUkfOHtKfCHFXUeO6YFlu5LuF45Y3fX8lvj1a08ar2pauALwGAIOABggaqvLNDShdqE2FzMceVr4r7I0oApH7UNv06rawxFt2t/ZoWpGHwQcYweX663Hd7Ys/sZDb2h/fdfxP/7GPm2LrUruM9F0LADwGgIOABggjDFd6jhqq51F450BR0ZF43lahyO66rtz8b/eFY5vqWuKB0V+n9FxIytzMr6efPKCafGA6Uhzuz5//0pZ635d73quc3bj0rmjNX5oRZ+MDQCKCQEHAAwgsx0BR0150NU9yRksZFIvka+UKsldOH6gl52qnOlUU4YPUmmg+45RuTJ+aIX++/LZ8etPrTugPyzpLCBfseOIFm89FL9+81mTBQBeRMABAAPI/LGD49+PH+pel8OZUtVdkXMid9F4bjsYjah2znD0LqVqraMlbl+kUzm96+TxOnfGiPj1r/5zjbbHUqh+6ajdWDRxiBZOGNLl/gDgBQQcADCAXDh7lE6bMkxVZQF99NzjXLfloi1uLms4JOUkpcrZErevF9Mzxug718zX4IpoIXhzW1if+fNy7TjUrIde3xPfz7lCOQB4TfE1AgcA9FpJwKfff/BUhSO2y3oKfmeXqiJMqdqfg5SqfK0w3p2R1WX62pVzdcvvlkmSlmw9rOvueiUe1I0fWq6L5tT2+bgAoFgwwwEAA1CyxduCvUipCoUjcsYmuSwalxJnODJPqapvaXetrj5zdN8HHJL05vlj9JYFY+LXOzpTSdKNZ0xmMT0AnkbAAQAe4XzTm25KVWI3q5wHHNXZpVStd8xuVJUGNHZweTd759f/e+scjXL8PJJUVRbQ209koT8A3kbAAQAeEXAEC+m2xW1LmAnJa0pVQ2uXtrI9ecOZTlVbJWMKN5MwuKJE33nbAte2d588QZWlZC8D8DYCDgDwCHfReHopVc4OVVJ+i8bbQhHVHwtldP91ewvXoSqZc6aP0MfOmyopujjgTbTCBQCKxgHAK9xtcdNNqXIHHLluizu4IqgSvy8+k7K/oUU1sY5P6XB2qJpZBAGHJH324pl650kTNGRQCbMbACBmOADAM5z1F6FezHD4jDstKxeMMRrRy9a41lpXh6qZfdwStzvjh1YQbABADAEHAHhE74rGnYv+5eeUMaKXnap2H21RQ0tnCtb0ArTEBQD0jIADADwi2IuUqnwu+tfB1Ro3g7U4nCuMjx1crpry9FOxAAB9h4ADADzCmQ6V7gyHM6Uq1x2qOvS2Ne7ahA5VAIDiRMABAB7h78XCf86ZkHylVCW2xk3X2iIsGAcAdEXAAQAe4ewwFepFl6pgID9rXDhTqvbVp1/DUWwtcQEAyRFwAIBH+H3OLlW9SKnK1wyHI6XqQJozHK2hsDYfaIpfn1lbPB2qAABuBBwA4BGuGY502+L2QZcqV0pVmjMcm/Y3xYOmoN9oyohBeRkbACB7BBwA4BEBX+ZF486UqrwVjTtSqprawmpq7Xm18XX7OtOppo6ozFswBADIHv+hAcAjAr0oGu+LlKphlaVyDC2twvG1ezoLxmcV0YJ/AICuCDgAwCMC/uJc+M/vc682vuvwsR7vQ0tcAOg/CDgAwCP8vVr4r3O/fKVUSdG0qA4b9jd0s2fUOlriAkC/QcABAB7hnKFIu2g8lP8ZDkmaPqozaFi/r7HbfY80t2mvo7icDlUAUNwIOADAI5w1HL1Zh6MkT+twSNK0UY4Zjn3dz3A406lqyoMa5WirCwAoPgQcAOARAVdb3DQDjj4oGpcSZzgaZG3q8a3d09mhamZtlYzJXyAEAMgeAQcAeERv2uL2xTockjR9ZGfAUd8S6rZT1Rt7qN8AgP6EgAMAPKJXbXGdAUcei8ZrKoKu9TjWd5NWtWr30fj3c8bW5G1MAIDcIOAAAI8IOIvG063hCDm6VOV5cb10CsdbQ2FXMDKPgAMAih4BBwB4RG9qONrC4fj3+WyLK0nHjey5cHz93sZ4S9+SgM91HwBAcSLgAACPcHWpSrMtrnOGI+jPb3F2YuF4Ms50qlm1VXmtKwEA5Ab/qQHAI1xF471pi+v353xMTtNdrXEbk3aqen1XZ8Axl3QqAOgXCDgAwCOcKVXtac5wtLqKxvM7wzHNMcPR0BpyLe7XYTUBBwD0OwQcAOARvVr4r4/W4ZC6LuKXWDjeHo7oDceif3PHEHAAQH9QFAGHMWaQMea9xpg7jDGvGGNajTHWGHNroccGAAOFs94hFLHdLq7Xoc210nj+TxnOOo7EwvEN+xrVFguAgn6j6bUUjANAfxAo9ABipkn6TaEHAQADmd/nTokKR6wrzSqZ9j5a+K/DtJFVem5DnaSuhePOgvHpo6pUGshvTQkAIDeKYoZDUoOkuyR9WNIiSV8p7HAAYOBJDC7SaY3r7lLVFzMcnbMWiSlVrvoN0qkAoN8oihkOa+0mSTd3XDfGXFTA4QDAgOTsUiWlF3C09nFKlbNwfOP+aKcqY6KBkqtD1TgCDgDoL4plhgMAkGddZjjCPXeqcheN57dLlSRNc8xwNLaGtPtotFNVOGK1Zk99/La5Y6rzPhYAQG4QcACARwR7McPR1zUc1WVBja4pi1/vqOPYfKBRLe3Rsfh9RrNGE3AAQH9RFClVuWaMWZ3ipql9OhAAKCKJRePptMbt6y5VUjStak9sZmPDvgadN2OkK51q2shKlQUpGAeA/oIZDgDwiEBiwJHG4n/OlKq+mOGQpOkjuxaOr9rVmU41h4JxAOhXcjLDYYx5QNKsDO92vbV2cS6eP5G1dk6y7bGZj9n5eE4AKHY+n5HPSB2ZVOnNcPRtlyop+Vocq1wrjJNOBQD9Sa5SqiZLmpHhfSpy9NwAgDQF/L744nnpzHC0hcLx70v7LKWqc4Zjw/5GhSNWqx1rcMwbywwHAPQnOQk4rLXH5+JxAAD5FfAZtcW+T69ovO9nOJytcZvbwnphY52a2qKBjzGiYBwA+hlqOADAQ5x1HOmkVLUXoGi8sjSgsYPL49cfWLYr/v2U4YM0qHRA9jsBgAGLgAMAPMQ5S9HTDEckYl37BPtgHY4OzrSqR1btjX9POhUA9D8EHADgIX7XDEf3NRxtCbeX9FFKleQuHD/W3llHMpeAAwD6naKZl451uhoduzom9vVmY8wlse/3WGuv6vuRAcDA4ZzhaO8hpao9MeDoo5QqKbrWRjK0xAWA/qdoAg5JCyVNTNg2NnaRpG19OxwAGHicMxzhHlKq2kLugKOvisYl9wyH0xxa4gJAv1M0AYe1dlKhxwAAA13AUYfR3kNb3MQZkL4MOI5LMsMxaViFqsuCfTYGAEBuUMMBAB4S9HX+2w9nmFLVl0Xjg0oDGjek3LVtDvUbANAvEXAAgIe4isZ7mOFodaRUlfh9MqbvAg6pax0HHaoAoH8i4AAAD3HOUmRSNN6XsxsdEus45lIwDgD9EgEHAHhIJkXjhVj0z2laQsAxZwwF4wDQHxFwAICHBFxtcXtYhyPknOHo+9PFwgmD1ZHFNX1UpYYMKunzMQAAslc0XaoAAPnnTI3qsS1uuLABx9QRlfralXP13Po6ffjcqX3+/ACA3CDgAAAP8Tu6VLVnsA5HaQFSqiTpPadM1HtOSVyiCQDQn5BSBQAeEnR2qeohpcpZVF6IGQ4AwMDAGQQAPKS3RePBQN93qQIADAwEHADgIUFX0Xj6KVUlzHAAAHqJMwgAeEjAVTTeQ5eqAheNAwAGBs4gAOAhzpSqTBb+K8Q6HACAgYEzCAB4SNDRpSrU0wwHKVUAgBzgDAIAHuJ3pFSFMikaJ+AAAPQSZxAA8BB3W9yeAg5HW1xSqgAAvcQZBAA8JOCYqeipLW4rKVUAgBzgDAIAHhJwFY33tPCfs2icdTgAAL1DwAEAHuJsi9tjSlWIGg4AQPY4gwCAh/hdXap6WPgvTEoVACB7nEEAwENcReM9tMV1damiaBwA0EucQQDAQ5xF4z3OcIQ6b2eGAwDQW5xBAMBDAq62uD0s/MdK4wCAHOAMAgAe4iwa76ktrrtonC5VAIDeIeAAAA9xt8VNf6VxUqoAAL3FGQQAPMRdw5F+ShVF4wCA3uIMAgAe4q7h6KlonHU4AADZ4wwCAB7iWvgvg3U4SpnhAAD0EmcQAPCQQAYL/7nW4WCGAwDQS5xBAMBDMmmL2+5Yh4OAAwDQW5xBAMBDXEXjPdVwsA4HACAHOIMAgIe4azh66FLFOhwAgBwg4AAAD3GlVGVQw8E6HACA3uIMAgAe4ioaJ6UKANAHOIMAgIdkklLVzjocAIAc4AwCAB6SycJ/7WG6VAEAsscZBAA8JN11OKy1LPwHAMgJziAA4CGulKpu1uFoT5j9YIYDANBbnEEAwEOC/vS6VLUnBCO0xQUA9FZRBBzGmJnGmM8bY54yxtQZY9qNMXuNMX81xpxV6PEBwEDhTzOlyrkGh0SXKgBA7wUKPYCYxyWNldQo6WVJhyTNlnSVpCuNMZ+21v6ocMMDgIHBWTQejlhZa2VM19mLrjMcBBwAgN4pljPIWknXSxphrb3QWnuttXaepA9LMpK+Z4yZXdARAsAAEEhIjUo1y9GWEHCw8B8AoLeK4gxirb3AWvtba21LwvafS3pUkl/S2wsyOAAYQJxdqqTUrXGdKVUBn5HPRw0HAKB3iiLg6MGK2NcxBR0FAAwAicXfqRb/Yw0OAECu9IezyJTY170FHQUADAD+hJmKVDMczhoOCsYBANkolqLxpIwxUyW9OXb17xncb3WKm6ZmPSgA6McSZyvaU8xwtDpSqpjhAABko2jPIsaYgKR7JJVK+qO19rXCjggA+r/EGY5wiqJx1wwHa3AAALKQkxkOY8wDkmZleLfrrbWLu7n9dklnStos6aOZPLC1dk6y7bGZD7pdAfCsAClVAIA+lquUqsmSZmR4n4pUNxhjviTpI5L2SbrYWnsoi7EBAGKMMQr4TLwdbsq2uKRUAQByJCcBh7X2+Fw8jiQZYz4s6WuSjkq6xFq7MVePDQCIplXFA45wqi5VBBwAgNwoqrOIMeadkn4iqVnS5dba5YUdEQAMPM4Aoj1FSpWzaJyUKgBANormLGKMuUzSbySFJF1lrX2hwEMCgAHJWTieumi8czurjAMAslEUbXGNMWdI+oskI+kd1tpHCzwkABiwnIv/pWqL60qpCtClCgDQe0URcEj6p6RySVskXWmMuTLJPs9ba3/Zp6MCgAEo4OucsUg1w+EsGmeGAwCQjWIJOAbHvk6OXVIh4ACALDlTqtopGgcA5FlRBBzWWubrAaCPOFOqUq3D0eZKqSLgAAD0HmcRAPCYdIrGnSlVpcxwAACywFkEADzG3RaXlCoAQH5xFgEAjwn4M2uLS5cqAEA2CDgAwGP8ji5V7Wl1qfLnfUwAgIGLgAMAPCbocxaNJ0+pamMdDgBAjhBwAIDHOIvGQ6lSqliHAwCQI5xFAMBjnEXg6bTFJeAAAGSDswgAeIy7aDyNLlWswwEAyAJnEQDwmIBrpfFUReOOLlXMcAAAssBZBAA8JuDoUhVKMcPhSqlihgMAkAXOIgDgMX5/pkXjdKkCAPQeAQcAeIy7LW6qhf9YaRwAkBucRQDAYwLOLlWpFv4jpQoAkCOcRQDAYwLpLPwXYoYDAJAbnEUAwGMCadRwMMMBAMgVziIA4DGuLlVp1HCw8B8AIBucRQDAY1wpVakW/mMdDgBAjnAWAQCPSactLilVAIBc4SwCAB4TdKVUpZrhcBaNsw4HAKD3CDgAwGMyLhonpQoAkAXOIgDgMYEeFv6z1pJSBQDIGc4iAOAx7oX/uqZUhSNW1hGHUDQOAMgGZxEA8JieZjjaE7YRcAAAssFZBAA8xt0Wt2vA4VxlXCKlCgCQHc4iAOAx7pSqJAFHQucqisYBANngLAIAHuNOqepaw9GesI22uACAbBBwAIDHuGY4ktRwOFOqfMa9PwAAmeIsAgAe467h6H6Gg4JxAEC2OJMAgMf0tPAfi/4BAHKJMwkAeEzAl35KFR2qAADZ4kwCAB7Tc0pVZxBCShUAIFucSQDAY1wpVUkX/nPUcAToUAUAyA4BBwB4jCulKkkNR2soHP+eGg4AQLY4kwCAx7hnOLqmVDW3dQYcFSWBPhkTAGDgIuAAAI8J9tClyh1w+PtkTACAgYuAAwA8xt9DSlVzayj+PQEHACBbBBwA4DHOLlXtyVKq2kmpAgDkDgEHAHiMs4YjnGSG4xgpVQCAHCqKgMMYM98Y82NjzMvGmN3GmFZjzFFjzEvGmI8bY4KFHiMADBQ9LfzX1ErAAQDInWKZKz9b0sckbZO0RtIBSSMknSHpVEnXGGMusta2FW6IADAwuIvGu6ZUHWvvrOEoJ6UKAJClYjmTPCTpIWvtZudGY8woSY9LOkfSByX9uABjA4ABxe+o4YhYKRKx8jm20aUKAJBLRZFSZa3dnBhsxLbvk/Tt2NXz+3ZUADAwBRMW82tPmOUgpQoAkEtFEXD0oD32lXQqAMgB5wyH1LVw3JlSRZcqAEC2ijrgMMYMkfSZ2NV/FXIsADBQBH0JMxwJheOkVAEAcqmoProyxkyT9CVFA6FRkk6XVCnpZ5Luy+BxVqe4aWq2YwSA/s7ZFlfqOsPRTEoVACCHiirgUDTIeF/Cttslfdla27WVCgAgY4kpVaGExf+aSakCAORQTs4kxpgHJM3K8G7XW2sXOzdYa5+PPpzxS5og6SpJ/yPp0lhb3K3pPLC1dk6Kca6WNDvDcQLAgNK1aDyhhsORUlXODAcAIEu5+uhqsqQZGd6nItUN1tqwpC2SfmCM2Srpfkl3SLqitwMEAEQlTHAo3E0Nx6BSAg4AQHZyEnBYa4/PxeOk8ICkRkmXGGNKWPwPALJjjFHQb+LF4s62uJGIdReNB0mpAgBkp6i7VEmStdZKOqRocDSkwMMBgAEh4OhU5SwabwmFXfuRUgUAyFbRBxzGmCmSxkuql1RX4OEAwIAQcORVtTuKxp2zGxIpVQCA7BVFwGGM+bgxpjbJ9hmSfifJSPpNrLYDAJAlZ2vckKOGw9kSV5LKAgQcAIDsFEty7mck/cgYs0LSRkUDjImSFikaFD0r6YuFGx4ADCx+R0pVyJFS5WyJWx70y5dYYQ4AQIaKJeD4kqTLJJ0o6WJJ5YrWbTwm6feSfss6HACQO0HXDEfylCoW/QMA5EJRBBzW2vuUwUriAIDsOFOqnEXjrlXGqd8AAORAUdRwAAD6lrNLlXPhv+Y2xyrjtMQFAOQAAQcAeJCzS5UzpepYO6uMAwByi4ADADzI7ww4XDMcrDIOAMgtAg4A8KCg39GlytEWt6nV2aWKlCoAQPYIOADAg9wzHI6UKrpUAQByjIADADwomGrhv3ZSqgAAuUXAAQAeFHAt/OdYh4OUKgBAjhFwAIAHOdfhSFU0TkoVACAXCDgAwIPcbXGTp1Sx8B8AIBcIOADAg/yulCrnSuPOhf8IOAAA2SPgAAAPcheNO2o4XClV1HAAALJHwAEAHhTwJ5/hOEZKFQAgxwg4AMCDUtVwOBf+o2gcAJALBBwA4EGBNBb+oy0uACAXCDgAwINStsVtpy0uACC3CDgAwINcC/+lKBpnpXEAQC4QcACAByWb4QiFI2oLdQYf5XSpAgDkAAEHAHhQsqJxZzqVxDocAIDcIOAAAA9yt8WNzmo4C8Yl2uICAHKDgAMAPCjZDIezJa7fZ1Ti5xQBAMgeZxMA8CBX0XishsO1ynjQL2NMl/sBAJApAg4A8KBkReOsMg4AyAcCDgDwIHdKVbSGw73KOB2qAAC5QcABAB7kLBpvj9VwuFcZZ4YDAJAbBBwA4EHOGY5wrEsVi/4BAPKBgAMAPChZDUdzW2dKFYv+AQByhYADADwo6OxSFU7epQoAgFwg4AAAD/I7i8aTpFRVlBBwAAByg4ADADzImVIVLxqnLS4AIA8IOADAg5wL/4UjXVcapy0uACBXCDgAwIPcMxzRlCra4gIA8oGAAwA8KOh3tsXtWjROW1wAQK4QcACAB/mdXao6UqpoiwsAyAMCDgDwoKCv+5Qq2uICAHKFgAMAPMjvI6UKANA3CDgAwIMC/s5//+1hVhoHAOQPAQcAeJC7aJyF/wAA+UPAAQAe5FppvGPhPwIOAEAeFG3AYYz5sjHGxi7XFXo8ADCQBJ0pVZGIrLWuLlUs/AcAyJWiDDiMMTMkfUmSLfRYAGAgSiwabw1FFHH8x2WGAwCQK0UXcBhjjKRfSDoi6e+FHQ0ADExBn7to3JlOJRFwAAByp+gCDkk3Szpb0mcUDToAADnmdxSNS1Jja8h1nZQqAECuFFXAYYyplfQdSU9Ya+8r9HgAYKByLvwnSUePtce/Lwn4XClXAABko6gCDkm3SyqX9JFCDwQABjLnOhySVN/SGXCQTgUAyKWiCTiMMW+W9HZJ37DWbij0eABgIEucwag/1plSNYh0KgBADhXFWcUYUynpTknrJX07B4+3OsVNU7N9bAAYCIIJNRzOGY5yZjgAADmUk4DDGPOApFkZ3u16a+3i2PffkDRe0pusta25GBMAILWuMxykVAEA8iNXMxyTJc3I8D4VkmSMOVnSxyT91lr7ZC4GY62dk2x7bOZjdi6eAwD6M2dbXEmqb3Eu+kfAAQDInZwEHNba47O4+2WK1pLMM8Y8nXDbzNjXLxljbpb0iLX2W1k8FwBAks9n5DOKL/bnnuEoimxbAMAAUUxnleO7uW1m7LK1T0YCAB4Q8PnUFo5IooYDAJA/Be9SZa291Vprkl0k/Tq223tj224o4FABYEAJOArHG1qcXaoIOAAAuVPwgAMAUBjOwnFSqgAA+ULAAQAeFXQs/ucsGielCgCQSwQcAOBRgRQzHKRUAQByqajnzWM1GzcUeBgAMCC5Ag5X0XhRnxoAAP0MMxwA4FEBR0pVYyvrcAAA8oOAAwA8yjnDYW3ndgIOAEAuEXAAgEc52+I60aUKAJBLBBwA4FEBX/JTADMcAIBcIuAAAI9KNcNBW1wAQC4RcACARzlrOJwGkVIFAMghAg4A8ChSqgAAfYGAAwA8ipQqAEBfIOAAAI9yrsPhREoVACCXCDgAwKOS1XAYI5UFOTUAAHKHswoAeFSygKM86JcxyVOtAADoDQIOAPCoZDUcLPoHAMg1Ag4A8KhkXaroUAUAyDUCDgDwqOQzHAQcAIDcIuAAAI9KVsNBwAEAyDUCDgDwqGRtcanhAADkGgEHAHhU0i5VzHAAAHKMgAMAPIqicQBAXyDgAACPCtIWFwDQBwg4AMCj/BSNAwD6AAEHAHhU8qJxAg4AQG4RcACARyVvi0tKFQAgtwg4AMCjWPgPANAXCDgAwKOCSbpU0RYXAJBrBBwA4FHJisYHkVIFAMgxAg4A8KjkbXGZ4QAA5BYBBwB4lJ+UKgBAHyDgAACPSlY0TkoVACDXCDgAwKOStcVlhgMAkGsEHADgUSz8BwDoCwQcAOBRwaQL/xFwAAByi4ADADwqWVtcVhoHAOQaAQcAeFQwIaUq4DMqCXBaAADkFmcWAPCoxBkOCsYBAPlAwAEAHpXYFpeWuACAfCDgAACPSkypomAcAJAPBBwA4FGkVAEA+gIBBwB4VNDnPgWQUgUAyAcCDgDwKGY4AAB9oSgCDmPMucYY283l5UKPEQAGmmBC0Tg1HACAfCi2+fNNkp5PsR0AkEOBLkXjxXZKAAAMBMV2dnneWntDoQcBAF4Q8DHDAQDIv6JIqQIA9L3EdTgIOAAA+UDAAQAelVg0TkoVACAfiu3sMs0Y801JwyTVKVrP8Yi1NlLYYQHAwJPYFpcZDgBAPhRbwHF67OL0ujHmGmvthnQfxBizOsVNU3s9MgAYYBJTqmiLCwDIh2JJqToq6buSTlV0dmOYpDdJelnSPEmPGmNqCjc8ABh4AsxwAAD6QE5mOIwxD0ialeHdrrfWLpYka+0yScsSbn/SGHOmpKcknSXpo5K+mc4DW2vnpBjnakmzMxwnAAxIXYvGi23SGwAwEOTq7DJZ0owM71PR0w7W2rAx5tuKBhwXK82AAwDQM9riAgD6Qk4CDmvt8bl4nBQ6ajdG5/E5AMBzjDHy+4zCESuJgAMAkB/FUsPRnSGxr00FHQUADEBBR1oVKVUAgHzoDwHHNbGvSws6CgAYgM6ZPkKSNHZwuaaPqizwaAAAA1FRfJxljPmkpPuttTsc24ykD0r6lCQr6aeFGR0ADFy3v2uhXtp0UAvHD1HA3x8+gwIA9DdFEXBI+qSk7xljlkraIqlM0Xa4kyVFJP2Htfa1wg0PAAam0oBf584YWehhAAAGsGIJOL4v6SJJcxRtWxuUtEfSvZJut9YuKeDYAAAAAPRSUQQc1to7JN1R6HEAAAAAyC0SdgEAAADkDQEHAAAAgLwh4AAAAACQNwQcAAAAAPKGgAMAAABA3hBwAAAAAMgbAg4AAAAAeUPAAQAAACBvCDgAAAAA5A0BBwAAAIC8IeAAAAAAkDcEHAAAAADyhoADAAAAQN4QcAAAAADIG2OtLfQY+owxpr60tLRq6tSphR4KAAAAkDObNm1Sa2trg7W2utBjSeS1gKNd0VmdtYUeC3KmI3rcVNBRIFc4ngMLx3Ng4XgOPBzTgWWmpIi1NljogSQKFHoAfWy9JFlr5xR6IMgNY8xqiWM6UHA8BxaO58DC8Rx4OKYDS8fxLEbUcAAAAADIGwIOAAAAAHlDwAEAAAAgbwg4AAAAAOQNAQcAAACAvPFUW1wAAAAAfYsZDgAAAAB5Q8ABAAAAIG8IOAAAAADkDQEHAAAAgLwh4AAAAACQNwQcAAAAAPKGgAMAAABA3hBwAAAAAMibARFwGGPOMMY8ZIw5ZIxpNMYsNsZcn8XjXWGMecYYUx+7PG2MuTzFvn5jzDuMMd8zxjxrjGkyxlhjzD29/oEGOGNMuTHm/xlj1htjWowxu40xdxtjxvbisYYYY24zxmwzxrTGvv7IGDO4m/v4jTGfMsa8bow5Zow5YIz5kzFmVlY/mEcV8ngaY2bEjuXvjTGbYn971hgzKdufy6sKdTyNMUFjzEXGmB8bY1YZY5pjf59vxP6/jsjJD+hBBf4bvcEY84fYcTxkjGmLPf9fjDFnZP3DeVChz6EJ9y8xxqyJ/d8NZfzDoNB/n/c4zpvJLh/O+gfsYK3t1xdJ10gKSYpIelrSXyQdlmQlfa8Xj/fJ2H3bJT0s6UFJzbFttyTZf3DstsTLPYV+bYrxIqlM0kux12i3pD9KeiV2fb+kKRk81nBJG2L33RR7rFWx6+skDU1yH5+kv8b2ORz7fXk69vvTJOnkQr9G/elSBMfzRyn+/iYV+rXpj5dCHk9JFziO3xZJ90v6u6QDsW17JM0o9GvU3y5F8Df6qqLn06Wx4/knSSti94lI+nChX6P+dCn08UzyGLfGjqOVFCr069PfLoU+npLuid3+SOz7xMt5OftZC/1iZ3mghko6GnuxrnZsH+V40c/N4PFmKBq8tEg6zbF9uqS62D/N4xLuM0jSbyT9h6TTJN0gAo7uXuOvxV6fFyVVOrZ/Orb96Qwe697Yfe6XFHBsvz3VMZB0c+y29ZJGObZfE9u+wflYXIr+eN4k6Vux4zdR0loRcPTL4ynp/NgJ8uSE7TWxk6GV9GKhX6P+dimCv9FTJFUl2f6W2Pn2mKThhX6d+sul0Mcz4f6zJLVK+rkIOPrl8VRnwHFu3n/WQr/YWR6oz8VeqAeT3HZV7LZ/ZPB4d8bu86Mkt30qdtsdPTzGO9P5Q/XiRVKJpCOx12dhkts7PvValMZjjZYUjv2zG5VwW6minwyEJI1MuG1N7DmuTPKYf4vddk2hX6v+cCmG45nkcQg4BtDxdNxnjDpnPyYW+rXqL5diPqax+z0ee/63FPq16g+XYjqekoyk5yTtkzREBBz98niqDwOO/l7D0VFX8Zckt/1L0ZmKC4wxZTl4vI5tV6Q/PCQ4Q9FPKzdZa5cluT2T1/gSRdOjnrPW7nPeYK1tlfQPSX5Jl3VsN8ZMVvQTmWOK/n5k8/wo8PFEzhXt8bTW7lY0tUqKBh9IT9Ee05j22Ne2DO7jZcV0PD8k6UxJn7HWHk7j+dBVMR3PvOvvAceC2NeliTdYa9sUzV0rUzQlqluxgpoJsatdDry1doeiaVUTjTHVvRyv16U8Xgnb5+fpsTrus8pa266uMnl+FP54IreK9njG/j8PiV3dm859IKm4j+mbFE2jOyzp5XTug+I4nsaY0Yqmsj5hrb03jedCckVxPGOuNsbcYYy50xjzWWPMzDSeMyOBXD9gX4m96a+JXd2ZYredkk5UNLd7ZQ8P2RFsHLbWNnXzeMNjj/d6+qNFTMdr3N3xkqKvbz4eK5fPj8IfT+RWMR/Pjyl6vnrdWrslzfugiI6pMeb9ks5R9EPAqYqem49Kepe19kgaz4/iOZ4/VvQ4fjSN50FqxXI8JenjCde/bYz5qaRPWGtz0n2s3wYckiod3zen2KcjcKjK4PFSPVamj4euenqNc3m8kj1WLp8fhT+eyK2iPJ7GmIWS/jt29fNpPDc6FdMxPUPS+xzXD0n6gLX232k8N6IKfjyNMW+VdLWk/7XWrk/jeZBawY+nohk9L0l6UtGgpFbSpYoWs39U0XTHT6Xx/D0qaMBhjHlA0Zz6TFxvrV2cj/EAAIqHMWaUom2syxRt5vFwgYeEXrLW3izpZmNMpaIdIT8n6X5jzP9Zaz9Y2NEhHcaYKkVnN9ZL+maBh4McsNbelrBpi6Q7jTHPKJqGdYsx5gexsoKsFHqGY7Ki/3gyURH72piwrT7JvoNiXxvSeNyOx6voZp9MHg9d9fQa5/J4JXusXD4/Cn88kVtFdTxjb24ekjRJ0p8lfSaN54VbUR1TSbLWNkp6TdK1sYYuHzDG/Ntae38aY/C6Qh/Pb0gaJ+mCWCEyslPo45mStXa1Mebvkt4m6U2KdrPKSkGLxq21x1trTYaXp2P3rVc0/1OK/gEk07F9WxrD2R77OsQYMyjFPpk8HrrqeI1zebwyeaxcPj8KfzyRW0VzPGNvRP8u6QRJj0q6zlobSeN54VY0xzSFjoLjt2ZwHy8r9PG8QtHun182xjztvMRu9zu2HZ/GGLyu0MezJxtiX0dncJ+UCj3Dka0Vks5W9KS0xnmDMSYoaa6ifxw95hlaa48YY7YrWnizUNLzCY83XtGC8W2xYAeZWxH7ekKK2zu291Tg39vH6rjPXGNMMEmnqkyeH4U/nsitojiexpiAogsAnqvoYlhXx7oOInNFcUy7URf7OiKD+3hZMRzPMkWL/1PpuG1wGmPwumI4nt3p6AyYqpFSRvp7W9yOtRTeluS2Nyv6h/G4tbYlB4/Xse0f6Q8PCV5QdFZqaopPPzJ5jR+RFJF0ljFmpPMGY0ypop/EhBVNyZAkxbrbvCGpXJ1rrvT2+VHg44mcK/jxNMYYSb9SdBXq5ZIu76ZrIHpW8GPag443p5syuI+XFfocOilV9klsl3BiNgq6VbR/n7H7dLxPStVqNzP5XlkwnxdJQxU9WFbRT8E6to9UdCoo6eqJiq5GvFbS2ITtMxRdibFF0qmO7dMU/SSmXdJxPYyJlca7f32+Fnt9XpA0yLH907HtTyfsf0vsWH0zyWPdG7vPXyQFHNtvS3UMJN0cu229HCtuKtp1w8Z+bwLZ/IxeuhT6eCZ5DFYa78fH03HbG5JGFPr1GAiXQh5TRZvCvENSScJ2EztXNiv6JumkQr9O/eVS6L/RbsZlxUrj/ep4Spop6b2SShO2j5D0QOw+yyWZnPyshX6xc3CwrlE0aoso2tbrz4ouJGQlfT/FfaxSvClRtP2XVTS4eEjSg7F/ilbSx1M83p2KLlz0sjoDnQOObS8X+nUqlouis04vx16j3YqmTnRc3y9pSsL+t6b6x6doitvG2O0bJf1B0fVROgKKoUnu41O0641VtC3jnyU9Ffv9aZZ0SqFfo/50KYLjeYLz70zRVeStoq3+OrbdXOjXqb9cCnk8Fc3j7/jf/KiiRYrJLjML/Tr1p0uBj+m5sduOSHpC0n2KZhJsiW0PS/pkoV+j/nQp9P/cbsZFwNHPjqfj7/NQ7H/ufYq+H6qPbd8haXrOftZCv9g5OmBnSHpY0UCjSdISSe/rZv+UAUfs9iskPatoNX9D7Ps3d/N4TzseM+ml0K9RMV0UTWn6f7E/iFZJexRNoxiXZN+Uf1yx24dKul3RgqnW2NfbJA3u5vn9in56sErRN6h1igYeswv92vTHSyGPp+MfZneXWwv9GvWnS6GOp6Qb0jiWVklmrbkU7TEdIenLigYbOxTNHmhW9M3PXZJOKPRr0x8vhT6HpngcAo5+djwljZH0Q0XX4dij6JobDYp2kbtV0pBc/pwm9qQAAAAAkHP9vWgcAAAAQBEj4AAAAACQNwQcAAAAAPKGgAMAAABA3hBwAAAAAMgbAg4AAAAAeUPAAQAAACBvCDgAAAAA5A0BBwAAAIC8IeAAAAAAkDcEHAAAAADyhoADAAAAQN4QcAAAAADIGwIOAAAAAHlDwAEAAAAgbwg4AAAAAOQNAQcAAACAvPn/YOtMgqVNdjsAAAAASUVORK5CYII=\n", + "text/plain": [ + "<Figure size 900x600 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig = plt.figure(0)\n", + "\n", + "plt.plot(signal['H1'].to_timeseries().sample_times, signal['H1'].to_timeseries())\n", + "\n", + "plt.xlim(-0.01,0.05)\n", + "\n", + "fig.set_dpi(150)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "metadata": {}, + "outputs": [], + "source": [ + "adv_psd = {}\n", + "for det in detectors:\n", + " adv_psd[det] = psd.analytical.aLIGOZeroDetHighPower(length=1024*8+1, delta_f=1./8, low_freq_cutoff=20)" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Help on function aLIGOZeroDetHighPower in module pycbc.psd.analytical:\n", + "\n", + "aLIGOZeroDetHighPower(length, delta_f, low_freq_cutoff)\n", + " Return a FrequencySeries containing the aLIGOZeroDetHighPower PSD from LALSimulation.\n", + "\n" + ] + } + ], + "source": [ + "help(psd.analytical.aLIGOZeroDetHighPower)" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "metadata": {}, + "outputs": [ + { + "ename": "AttributeError", + "evalue": "module 'pycbc.psd' has no attribute 'delta_f'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m<ipython-input-58-e64cb88395a2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mdelta_t\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1.0\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;36m4096\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mtsamples\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m32\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mdelta_t\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mts\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpycbc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnoise\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnoise_from_psd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtsamples\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdelta_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpsd\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mseed\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m127\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;32m/work/francisco.jimenez/sio/pyenv/versions/3.6.6/lib/python3.6/site-packages/pycbc/noise/gaussian.py\u001b[0m in \u001b[0;36mnoise_from_psd\u001b[0;34m(length, delta_t, psd, seed)\u001b[0m\n\u001b[1;32m 102\u001b[0m \u001b[0mrandomness\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlal\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgsl_rng\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"ranlux\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mseed\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 103\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 104\u001b[0;31m \u001b[0mN\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mint\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1.0\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mdelta_t\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mpsd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdelta_f\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 105\u001b[0m \u001b[0mn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m//\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 106\u001b[0m \u001b[0mstride\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m//\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mAttributeError\u001b[0m: module 'pycbc.psd' has no attribute 'delta_f'" + ] + } + ], + "source": [ + "delta_t = 1.0 / 4096\n", + "tsamples = int(32 / delta_t)\n", + "ts = pycbc.noise.noise_from_psd(tsamples, delta_t, psd, seed=127)" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "72.31060179760523\n" + ] + } + ], + "source": [ + "optsnr = {}\n", + "\n", + "snrsq = {}\n", + "addsnr = 0\n", + "for det in signal.keys():\n", + " snrsq[det] = filter.sigmasq(signal[det], psd=adv_psd[det],\n", + " low_frequency_cutoff=f_lower,\n", + " high_frequency_cutoff=f_final)\n", + " addsnr += snrsq[det]\n", + "optsnr = numpy.sqrt(addsnr)\n", + "\n", + "print(optsnr)" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "metadata": {}, + "outputs": [ + { + "ename": "NameError", + "evalue": "name 'i' is not defined", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m<ipython-input-56-8af46ede6ee5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0m_\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgauss\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0madv_psd\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdet\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 3\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mNameError\u001b[0m: name 'i' is not defined" + ] + } + ], + "source": [ + "for _ in range(10):\n", + " value = gauss(0, np.sqrt(adv_psd[det][i]))\n", + " print(value)" + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "metadata": {}, + "outputs": [], + "source": [ + "# The color of the noise matches a PSD which you provide\n", + "flen = int(2048 / delta_f) + 1\n", + "psd = pycbc.psd.aLIGOZeroDetHighPower(len(signal['H1'].to_timeseries()), delta_f, f_lower)\n", + "\n", + "# Generate 4 seconds of noise at 4096 Hz\n", + "tsamples = int(4 / delta_t)\n", + "ts = pycbc.noise.noise_from_psd(tsamples, delta_t, psd, seed=127)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "lengths do not match", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m<ipython-input-78-d70d9b4b69c6>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mpylab\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mts\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msample_times\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mts\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0msignal\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'H1'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_timeseries\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mpylab\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mylabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Strain'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mpylab\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mxlabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Time (s)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mpylab\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m<decorator-gen-154>\u001b[0m in \u001b[0;36m__add__\u001b[0;34m(self, other)\u001b[0m\n", + "\u001b[0;32m/work/francisco.jimenez/sio/pyenv/versions/3.6.6/lib/python3.6/site-packages/pycbc/types/array.py\u001b[0m in \u001b[0;36m_returntype\u001b[0;34m(fn, self, *args)\u001b[0m\n\u001b[1;32m 230\u001b[0m \u001b[0;34m@\u001b[0m\u001b[0mdecorator\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 231\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_returntype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 232\u001b[0;31m \u001b[0mary\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# pylint:disable=not-callable\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 233\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mary\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mNotImplemented\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 234\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mNotImplemented\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m<decorator-gen-153>\u001b[0m in \u001b[0;36m__add__\u001b[0;34m(self, other)\u001b[0m\n", + "\u001b[0;32m/work/francisco.jimenez/sio/pyenv/versions/3.6.6/lib/python3.6/site-packages/pycbc/types/array.py\u001b[0m in \u001b[0;36m_convert\u001b[0;34m(fn, self, *args)\u001b[0m\n\u001b[1;32m 64\u001b[0m \u001b[0;31m# Convert this array to the current processing scheme\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 65\u001b[0m \u001b[0m_convert_to_scheme\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 66\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 67\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 68\u001b[0m \u001b[0;34m@\u001b[0m\u001b[0mdecorator\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m<decorator-gen-152>\u001b[0m in \u001b[0;36m__add__\u001b[0;34m(self, other)\u001b[0m\n", + "\u001b[0;32m/work/francisco.jimenez/sio/pyenv/versions/3.6.6/lib/python3.6/site-packages/pycbc/types/array.py\u001b[0m in \u001b[0;36m_checkother\u001b[0;34m(fn, self, *args)\u001b[0m\n\u001b[1;32m 251\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mtype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mArray\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 252\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 253\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'lengths do not match'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 254\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprecision\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprecision\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 255\u001b[0m \u001b[0m_convert_to_scheme\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mValueError\u001b[0m: lengths do not match" + ] + } + ], + "source": [ + "pylab.plot(ts.sample_times, ts+signal['H1'].to_timeseries())\n", + "pylab.ylabel('Strain')\n", + "pylab.xlabel('Time (s)')\n", + "pylab.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "ename": "NameError", + "evalue": "name 'get_lm_f0tau' is not defined", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m<ipython-input-6-32b2ef174fcb>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mget_lm_f0tau\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmass\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mspin\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ml\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnmodes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;31mNameError\u001b[0m: name 'get_lm_f0tau' is not defined" + ] + } + ], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "from pycbc.conversions import get_lm_f0tau\n", + "import qnm" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "def wRD_to_f_Phys(f,M):\n", + " c=2.99792458*10**8;G=6.67259*10**(-11);MS=1.9885*10**30;\n", + " return (c**3/(M*MS*G*2*np.pi))*f\n", + "\n", + "def tauRD_to_t_Phys(tau,M):\n", + " c=2.99792458*10**8;G=6.67259*10**(-11);MS=1.9885*10**30;\n", + " return ((M*MS*G)/c**3)*tau" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "omegas = [qnm.modes_cache(s=-2,l=2,m=2,n=i)(a=0.7)[0] for i in range (0,2)]\n", + "w = (np.real(omegas))/1\n", + "tau=-1/(np.imag(omegas))*1" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.004205658110224193" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "tauRD_to_t_Phys(tau[0],69)" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [], + "source": [ + "frompycbc=get_lm_f0tau(69, 0.7, 2, 2, 2)" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.00420654, 0.00139151])" + ] + }, + "execution_count": 35, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "frompycbc[1]" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.0034787750000000004" + ] + }, + "execution_count": 38, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "0.00139151*2.5" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.6.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}