{ "cells": [ { "cell_type": "markdown", "id": "fdccdd1c", "metadata": {}, "source": [ "# Working with Data and Models From the Literature" ] }, { "cell_type": "markdown", "id": "dc7ef8da", "metadata": {}, "source": [ "A very incomplete set of data from from the literature exist in ``$ARES/input/litdata``. Each file, named using the convention ``.py``, is composed of dictionaries containing the information most useful to *ARES* (at least when first transcribed). \n", "\n", "**NOTE:** May be worth interfacing with [CoReCon](https://github.com/EGaraldi/corecon) at some point! \n", "\n", "To see a complete listing of options, consult the following list:" ] }, { "cell_type": "code", "execution_count": 1, "id": "ef0cfb56", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "%pylab inline\n", "import ares\n", "import numpy as np\n", "import matplotlib.pyplot as pl" ] }, { "cell_type": "code", "execution_count": 2, "id": "e9c0ee03", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['ueda2003',\n", " 'mirocha2016',\n", " 'test_schaerer2002',\n", " 'song2016',\n", " 'whitaker2012',\n", " 'blue_tides',\n", " 'atek2015',\n", " 'mortlock2011',\n", " 'bouwens2017',\n", " 'sazonov2004',\n", " 'dunne2009',\n", " 'parsec',\n", " 'eldridge2009',\n", " 'furlanetto2017',\n", " 'mirocha2017',\n", " 'finkelstein2012',\n", " 'noeske2007',\n", " 'leitherer1999',\n", " 'alavi2016',\n", " 'haardt2012',\n", " 'tomczak2014',\n", " 'feulner2005',\n", " 'bpass_v1',\n", " 'aird2015',\n", " 'mirocha2018',\n", " 'marchesini2009_10',\n", " 'emma',\n", " 'vanderburg2010',\n", " 'mirocha2019',\n", " 'duncan2014',\n", " 'daddi2007',\n", " 'starburst99',\n", " 'park2019',\n", " 'oesch2014',\n", " 'eldridge2017',\n", " 'bowler2020',\n", " 'calzetti1994',\n", " 'oesch2013',\n", " 'sanders2015',\n", " 'kajisawa2010',\n", " 'lee2011',\n", " 'stefanon2017',\n", " 'stark2011',\n", " 'mcbride2009',\n", " 'madau2014',\n", " 'parsa2016',\n", " 'stark2010',\n", " 'oesch2016',\n", " 'morishita2018',\n", " 'perez2008',\n", " 'ferland1980',\n", " 'bpass_v2',\n", " 'bowman2018',\n", " 'sun2020',\n", " 'gruppioni2020',\n", " 'rojasruiz2020',\n", " 'stefanon2019',\n", " 'kusakabe2020',\n", " 'finkelstein2015',\n", " 'moustakas2013',\n", " 'bouwens2015',\n", " 'mirocha2020',\n", " 'reddy2009',\n", " 'mclure2013',\n", " 'robertson2015',\n", " 'mesinger2016',\n", " 'kroupa2001',\n", " 'karim2011',\n", " 'schneider',\n", " 'behroozi2013',\n", " 'schaerer2002',\n", " 'weisz2014',\n", " 'oesch2018',\n", " 'gonzalez2012',\n", " 'inoue2011',\n", " 'bouwens2014',\n", " 'ueda2014']" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ares.util.lit_options" ] }, { "cell_type": "markdown", "id": "03c4ff49", "metadata": {}, "source": [ "If any of these papers ring a bell, you can check out the contents in the following way:" ] }, { "cell_type": "code", "execution_count": 3, "id": "86bde8a8", "metadata": {}, "outputs": [], "source": [ "litdata = ares.util.read_lit('mirocha2017') # a self-serving example" ] }, { "cell_type": "markdown", "id": "29d4f067", "metadata": {}, "source": [ "or, look directly at the source code, which lives in ``$ARES/input/litdata``. Hopefully the contents of these files are fairly self-explanatory! \n", "\n", "We'll cover a few options below that I've used often enough to warrant the development of special routines to interface with the data and/or to plot the results nicely." ] }, { "cell_type": "markdown", "id": "a6184073", "metadata": {}, "source": [ "## The high-z galaxy luminosity function" ] }, { "cell_type": "markdown", "id": "b0381cdf", "metadata": {}, "source": [ "Measured luminosity functions from the following works are included in *ARES*:\n", "\n", "* Bouwens et al. (2015)\n", "* Finkelstein et al. (2015)\n", "* Parsa et al. (2016)\n", "* van der Burg et al. (2010)\n", "* ... many more now" ] }, { "cell_type": "markdown", "id": "3e6468ef", "metadata": {}, "source": [ "## Stellar population models" ] }, { "cell_type": "markdown", "id": "eadae667", "metadata": {}, "source": [ "Currently, *ARES* can easily handle both the *starburst99* original dataset and the *BPASS* version 1.0 models (both of which are downloaded automatically). You can access the data via" ] }, { "cell_type": "code", "execution_count": 4, "id": "5c0027b4", "metadata": {}, "outputs": [], "source": [ "s99 = ares.util.read_lit('leitherer1999')\n", "bpass = ares.util.read_lit('eldridge2009')" ] }, { "cell_type": "markdown", "id": "3e0e0c7b", "metadata": {}, "source": [ "or, to create more useful objects for handling these data" ] }, { "cell_type": "code", "execution_count": 5, "id": "c26233c9", "metadata": {}, "outputs": [], "source": [ "s99 = ares.sources.SynthesisModel(source_sed='leitherer1999')\n", "bpass = ares.sources.SynthesisModel(source_sed='eldridge2009')" ] }, { "cell_type": "markdown", "id": "42e3c77e", "metadata": {}, "source": [ "The spectra for these models are stored in the exact same way to facilitate comparison and uniform use throughout *ARES*. The most important attributes are ``wavelengths`` (or ``energies`` or ``frequencies``), ``times``, and ``data`` (a 2-D array with shape (``wavelengths``, ``times``)). So, to compare the spectra for continuous star formation in the steady-state limit (*ARES* assumes continuous star formation by default), you could do:" ] }, { "cell_type": "code", "execution_count": 6, "id": "63f4c2df", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "# Loaded fig8b.dat\n", "# Loaded $ARES/input/bpass_v1/SEDS/sed.bpass.constant.nocont.sin.z020\n" ] }, { "data": { "text/plain": [ "(1e+33, 1e+41)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVsAAAD8CAYAAADZoQcPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAArNklEQVR4nO3dd3xUVfrH8c+TSSUQIPQQINRQlCYCoiCwgGDBtmvBsiqKsotlXX+Kbe3ddRVlRSxrRZZl1aUJAoKAIkhRWightFBDDSGEtPP7YwYySSbJhGTm3pl53q8Xr8w9987Mc5nkm5Nz7z1XjDEopZTyrTCrC1BKqVCgYauUUn6gYauUUn6gYauUUn6gYauUUn6gYauUUn4QbnUBntSvX98kJSVZXYZSSlXKypUrDxpjGnhaZ8uwTUpKYsWKFVaXoZRSlSIiO8pap8MISinlBxq2SinlBxq2SinlBxq2SinlBxq2qkLZufnMT9nP6p1HuPKdJRzKOmV1SUoFHFuejaDs5e7PVrJ4y0HiosPJzMnn+ok/M+eBfjjCxOrSlAoY2rNVFVq85SAAmTn5AKQeyKL1Y7PIOpVvZVlKBRQNW3XW7v9ytdUlKBUwNGxVufKzDvFpxEs04GipdfM3HmDUpytIGjuT1+Zs9H9xSgUQDVtVrtzl/6KfYy0vRnxAkuwFYNWTg8+s/27DfgDGL9jKq7M1cJUqi4atKtfJvEIABjtWsTDqr2yPHkF8RD7N42uU2vafC7eSk1fg7xKVCggatqpcHsPzxSZEhRd960y4+Ty6Na8DQPsnZ7PzULafqlMqcGjYqnKd7tmWNDdzOOdIGgDnNI3jlWs7n1l33XtLyS/w/DylQpWGrSpX3snMMtfNiHqCEY75xJ9Io3ZMxJn2fZk5fL16tz/KUypgaNiqcnXYPKHc9S9GfEiNDy6kwboPOE82cXNHZ+j+39Q1pB/R4QSlThNjjNU1lNKjRw+j89naxNO1K7W5iahBy+MfFGtb+uhAmtSOqc6qlLIlEVlpjOnhaZ32bFW1krxs+rSuV6ztgpe+5/CJXIsqUsoeNGxVtXuq0Y9sjx7BTe2Lpt7o/txcCgvt91eUUv6iYRsCcvIK2JqRVX0veOPkclcnr3oWgBf6RjBtzIVn2oe9tbj6alAqwGjYhoDxC1L53d9/YMv+47wxdzP3fLaSBZsO8PK3G5m3YT8z1uwhJ6+AvIJCKKzgooSajSF5GETWKmq7YpznbU0hnRPrsPG5oQBs2n+cpLEzOaJDCioE6RSLISAjM4cYcnh+Zgo/bt5LTU4yf306eSU+/u7h2/gq/HG45WtoPdDzi93rOnBp3EK56wg4lQnfPVF828+vhRH/IbrdEF665lwe/WotAN2em8us+/rSMSGuunZRKdvzSc9WRGJFZKWIXO62/ImIvC8iN/niPZUHx/dD7gkuOjadlOg7MKnzSI2+lV+j7+Y/kU+X2rybSQHg18nPMmXFLjyeqRLl6tGec01RmyMCwqM91zDpD3BgIzf2bM4zwzudab503GIyjusk5Cp0eHXql4h8BFwOHDDGnOPWPhR4C3AAHxhjXna1PwucANYbY2aIyC3AUWPMdBH5tzHm+vLeT0/9Onv5BYWEO1y/Qys4beuDLv8m61QhG9YsZ3DYSmrHRDAkbz4AvXPeZh/12B49oviTnj7m/FqQB8fSIcwBdZrDtsXwyeVlv9llb8D5I1my5SA3f7jsTHPai5cSppOQqyBR3qlf3g4jfAy8A3zq9qIOYDwwGEgHfhGRaUACsAFw7+okAmtdj3WmEh/5du1eRn+xilt6t2Dz/uP8u4Lt7/zN9Tsv0tWQV7QuUTLYZ+qVes4ZjgiIb1m03LIvXPEWdLwSXkkqvf3MB6H7H7mobX0WPzyAvq8uAKDHC/OKzSKmVLDyahjBGLMIOFyiuSeQaoxJM8bkApOBK4EBQG9gBHCXiIThDOPEyrynqrzpa/YA8NnPO1i2reTHVTlTo55lZfSfijdGVXCBw3m3QUzdstc/Vw+MoVl8DVKedR40O3wil12H9UozFfyqEnxNgV1uy+lAU2PM48aYB4BJwPvGmELgK+BaEXkXmO7pxURklIisEJEVGRkZVSgrNBUUGlZsSefGRGfItqhXegrEyqpXcsLwjsMr9wJjd5Vue6YOZB8mJtLBe7ecB0DfVxdwQm+xo4JcVcLW00DbmQFgY8zHxpgZrscnjDG3G2NGG2O+8PRixpiJxpgexpgeDRo0qEJZoSf9SDZr0o/ycsEbvHRwDL9ed4rZ911kdVkQFg5/nA7XfVa8/dWWsOdX6teMOtPU6ak55ObrTGEqeFUlbNOBZm7LicCeqpWjKquw0DDszcVc/c+f6BG2GYA6024n5qX61ftG7S+HAY9X7jmOCGjZz3OPeOLFdI8/xXU9Es80tXviW+a67vygVLCpStj+ArQVkZYiEgncAEyrnrKUt/Zm5nDcH3+C3/AFxDXxbttaCc6vjqJpF8+cxeBG/p7Mq1d35IFBbc+03fXpClL2lj2to1KByquwFZEvgaVAsoiki8hIY0w+MAaYA6QAU4wx631XqvJkx6ETZx7HiQ8ONEWcxdjvnfPg5v+Wbn/6GIzdWbztufrcN7At3/2lH5d1dob5sLcWMz+l/B7uyVw9qUUFFm/PRrjRGNPEGBNhjEk0xnzoap9ljGlnjGltjHnBt6UqT3a4bkHj3jusVveuhDu/r9xzajeFNoM8r4uu7QzdJ4oOgoYd2kK7RrUYP6L7mbaRn6zgp60HSz09LSOLpLEz6fC32XywOK1ydSllIT0NK8BtP3SCCIdw70AfhW1cAiSeV/2vGx4JN7l6v+PPh72/AXBjz6LDAHd+soJlaYeKPe2ycUvOPH5+Zgoz1uhhAhUYNGwD3KZ9x/kt4g4cn19V9ReLLnEebY+RVX/N8tRvU/T4vX6wcRZP963JvAcvZvljvyOhTgw3vP8zs9ftO7PZSdcNKBtyBIAxk1aTNHYme46e9G2tSlWRhm2AOpqdy4HMHNamH6MGJyFt4dm9UBe3y3EjYouvEx9/e9RNgjFul2VPvpGo8d1o0yCWhnHRfDayJw4RRn+xstgY7tVhi1ke/WfmRT5EomTwdsQ4hrw8g8ycvNLvoZRN6KxfAeLEqXxGfvILDWtFU1BoWLjpACdOHyQqYw4Yr/QaBb9N8rwuzFGFF/ZS/bbw1FHnxQ6nPVMH+j9Gk4sfZumjv+PWj5YzZtJq/nF9V2qQwz8i3wWgTdgelkTdD8AVjp+5+aV9fPLUGBw614KyIe3Z2lxOXgEPTF5Np6fm8HPaYab9toeZa/eyPux65jZ4k1Z1qvgRhoVD/0chLhG3a1KcxA9hCyDiPGh2z49FbQtfhOn30yA8m0/uOJ+k+rHc8/lKNkTfUebLfC5PMOFvt3qerUwpi2nY2lBaRhZz1u9jf2YO/1mZzje/Fh0E6pkUz3ODGwPQ9vhyvi8oI3we2+MMsOYXlP9mYeHQfyw86OGsPfFzD7HxOcWXV38GryTRcOo1TB5Z9kG6/MveJFOcUz/+OXwa8kwdCubryTHKXnQYwYaemraexVuKTntyhAnPDO9EYt0Y+ic3hOzDcPoOM3llnFt7+vzYO2ZDYSGkzoVJ15XeLsztW6Bkj9AfwwglDXoa5j0ND2+D3Svhi9/Djh+p/W5nNiafAzvctr19Nnw6nPDWF3M4LJ64guNnVjkWvwqLX3XOQnbth8UvsFDKAhq2NnQ0O49a0eE8OLgd2bkFdGtehz6t3S6/rejWNeffVbxXGhZWPFRbDYC0Ba515QSqI6rsdb5y4QPQ+08QHgVtB8OQF2Duk3Aig+gTC4pv2+ICeNJ5vm5qsz+QtP11suI7UfOwWy99w//gQAr8aZnz/0Epi2jY2tCp/AIubF2f2y9s6XkDU0HYtuhTus09fNsNLQrb8sZlz7m2/PfxBRFn0J7WZwx0vxXW/RdmPFDm0wbc8jg7D91H84bO09eueHsJTffOZULkm3BwMzxbF3qOgktf8239SpVBw9aGTuYVEBNZTggu+Uf5L+B+yxpP3E/pKtazdRtG8DCXgWWi46DH7dCsJyAQUweO7S62icMRdiZoAabfexEvzIznoiUtz5yxwPKJzn+3ToPDadDtZh1eUH6jYWtDJ3PLCNu8k/BC4/KfHNuw4jfw94Gv6tKo6B5mxCVUuPnjl3UkKtxB0oJJxW/v86lrFrIa8c7JzpP6Bu7/iQoYGrY2lJ1bQI0ID2F7cEvFTy70YgawYr1Zt5AJwlOmHrokmSWpB7l19yNkmlgmNZ9Gjf2uCymm3Or82vhcaNrDORVkWXcVVqqK9IiBzRhjSg8jbP7OefPGozvKfuJpPW4v64WLHrsPI/j6KjEb+ObPF7KosAu/mjZ03/FnNv9hISR0K9pg31pY+S/47Grn/7M3v9SUqqTg/0kLMKfyCzGG4mG78mPn190rK36BpL4Vb+MesBExbiuCr2dbUg5RDPlsD6+1mIB56ii0uLD0Ru/0gMk3QY7Oq6uqjw4j2MzxHOcwQGyk20dz+pSlAm+u/fcmMF1DB407Ow82hYDJo3pTMyqcr1fv5sMl2xi/YCurdx5l0l2zIGU61GwEG2fAj285n7BxBvz7ONw01TlDmVJVpGFrM/uO5QDQKM5twoPT58h6E7amrPt4eQjhxp1LbBK8PdverZy3ZT+naW0eHNyOWz9azk9bD5E0diarn7yEurGRzrMdBj/rfMKvk+Cb0fCfPzpv0V7TiwOPSpVDhxFsZu8x51SBCXXcwvb0ubA7l1b8ArHe3CyzglC9c74XrxG4YqPC+XxkrzPL3Z6bS9LYmcXvf9Z1BAx7DbZ8B6+3hSVv+r9QFVQ0bG1mr6tn26S221jq6Z7tvjUVv0CTLp7b67pdIHHuH6DrzTD4Gc/b1mnuRaWBLSbSwfaXLyvWdtenK3j867VFDb1GwagfnI/nPQULXwnq3r/yLQ1bm0k/kk1keBj1Yt3GCatjjoJ6rYseR8TAVeMhtprvwBuAJt3Zq9jyF8t20ulvs4saGp8DY1Y6D6QtfBFeTIBcH9zrTQU9DVubScs4Qav6sYS5z8la3ulZna8venz3It8VFqT6tKnPvAf78dil7c+0ncgt4PkZG4o2qt8GbpsJzXo5J/7511A4utPDqylVNg1bH8o6lU9Bofd/dhYUGtbtOUbrhjWLGvNzy39STF3n2OL9a8oeQvBWiF5F1aZhLUb1a82PYwdygetA2gdLtjFzzd6ijURg5HcwYorzfmlvnguHt1U8KZBSLno2QjU4lHWKzfuzSM3IInX/cbYcyCL1QBYHjp+iVnQ4w85pzFVdm9KrVb1y7yKwaEsG+zNPMewct0tyn/figFevUdWwF6ppnRi+HNWbnYeyufvzlfx50irCHedxSSe3z6PdJXDxWPjhZRjXtaj97sXQpHOp11TqNA3bKvpgcRovzEo5c9ykZlQ4rRvWpF+7BrSsH8vWjCxmrtnLlBXpNIqL4vLOCQzvkkDnxNpIiZ7kpGU7qV8zkiEdK5j/AKDPffDTuModsGnYEQ5sKHv9TVNhxUdQI7THcpvXq8E/b+rOgNcXcvdnK3n9D134/XmJRRv0HwsFuc7ZxDbOcLa91xeGvwPdb7GmaGV7GrZVsPfYSV6dvYmL2tRnVL9WtGlYk8Zx0aVC9ORVBcxL2c+03/bw2dIdfLhkG+0a1eT685vTIr4GDoeweudR5qXs556LWxMZ7sXoTlxT14NKhO2d88o/uJPQFYaP8/71gljL+rF8/ac+XP3Pn3joP7+x89AJHhyS7FwpAoOecj4+sh3ecg3fTL8fEntAww6W1KzsTcO2CiYuSqPAGF68+lyaxdcoc7uYSAdXdEngii4JHDuZx6y1e/li2Q6em1G8l3l55ybcN7Ctd29+OtAr07ONjHX+U17p1rwuP/xffy5+bSHjvk+lZ8t6XNS2RK+/bpJzOsoTh+DtbjDnMbj5q5Ad/1Zl07A9SxnHT/Hl8p1c3a1puUFbUu2YCG7s2Zwbzm/G7qMnOXwil/xCQ62ocNo0rFmqV1y209vpeZ++1KJeLFPuvoDr3lvKzR8u49M7etKvnYdx9Nh6zrHcOY9CyjTn7XiUcqNnI5ylD5ds41R+IX/q37rijT0QERLr1qBzYh26N69L20a1KhG0nF3PVp2Vni3jWfqoc+rFWz9azrrdZUysfv6dzjNCZj4EJ4/4sUIVCDRsz8LR7Fw+W7qdyzsn0KpBzYqf4AtdR0CXETDgcWveP8Q0qR1DX9cQwhXvLOHFWSmlNwqPhOFvQ/Yh+Hq0/iJUxWjYnoWtGVnUjA7nzwPOrldbZTdNdY69Xv2u889X5RefjezF0kcH0qxuDSYuSiNp7ExMyUBt0gUueRE2fwvLJlhTqLIlDduzcF6LeH58ZCDtG8f59o2yD3tubzvYt++rytSkdgzTx1x0Zrnbc3PZffRk8Y163Q3thsHcvzkvgFAKDduzFu7ww3/dqk99/x6q0mrXiGDbS5dyV9+WHM3O48KXv2fsf90mCRKBK8dDjXrw71sg64B1xSrb0LC1s/xTVlegyiAiPH5ZxzPLk3/ZxTvfu91OJ7YeXP+FM2j/OxIKy5pnWIUKDVs7y8+xugJVgdQXhrHgof4AvP7dZo5lu03wnngeDHsFti2CZe9aU6CyDQ1bOyuoYBIaZblwRxgt68fy5V29AXjsm7XFN+h+q3P89vvndTghxGnY2pn2bAPGBa3rMapfK2au2cvCTW6hKgIXPeCcmnH155bVp6ynYWtH+aec9xsrdN78kQi3S2z7P2ZNTapCDw1JJqleDe75fCV5BW5jtFG1nF/nl3FnDBUSNGzt6PmGMLF/0bIjAh7bA08egv6PWFaWKl9keBh/GdyOnLxCrhr/Y9GKhkUH0vSgZ+jSsLWr/evAEeV8fNsM50UMDp3Kwu6Gd0kAYP2eTDbvP+5sdL8M+7snLKhK2YGGrZ0V5jnnlm18rtWVKC+JCKueHEztmAju+PgXsnNdQ0G3u+5rtnyidcUpS2nY2llBnnMIQQWU+NhI3rqhK+lHTvLcDNccCi0ugCjXFYeZe8t+sgpaGrZ2VpgPYRq2gah/ckOu6prAl8t3Fp2dcMvXzq97VllXmLKMhq2dFeTpOG0Ae+mazrRqEMtD//mNnLyCortrTB4BR3ZYW5zyOw1bOyvM055tAIuJdPDk5R05mJXL63M2Qc2GRSs/vxZyypgXVwUlDVs7K8jXMdsANyC5Ib1bxfPJ0u3sOnoKHkyB6z6Fw2nwzZ90ztsQomFrZycOQJgOIwS6R4a2J6/A8MlP2yEuwXnLnMHPOO/M+9PbVpen/KTaw1ZEOojIBBGZKiKjXW3NRWSaiHwkImOr+z2DVvov2rMNArFRzl+YHyzZVtR4wRjoMBzmPQ3bf/T8RBVUvApbV0geEJF1JdqHisgmEUk9HaLGmBRjzD3AdUAP16btgJnGmDuAjijv6d1wA167RrXOPD6e45oV7PSct3WTYOodzrvzqqDmbc/2Y2Coe4OIOIDxwDCcAXqjiHR0rRsOLAHmuzZfDdwgIt8DC6pedggJj7G6AlWNzn36u6KF6Dj4w8fOA2d6g8ig51XYGmMWASXv0dITSDXGpBljcoHJwJWu7acZY/oAN7m2vR14yhgzELjM03uIyCgRWSEiKzIyMs5iV5Syrwk3n3fm8ZQVu4pWNOkMdy+C+m0sqEr5U1XGbJsCbt81pANNRaS/iIwTkfeAWa51s4H7RGQCsN3TixljJhpjehhjejRo0KAKZQUZozP8B4Oh5zSmUZxzrouHp65h7zG3+5ZV5hb2KmBV5VC3p+8QY4xZCCws0bgO+H0V3it0adgGjWWPDSJp7EwALnjpe7a8MIwIf9zLTtlCVT7pdKCZ23IisKdq5ahSNGyDyraXLj3z+Kb3l1lYifK3qoTtL0BbEWkpIpHADcC06ilLnaFhG1REhK0vOgN3+fbDTP9N+yehwttTv74ElgLJIpIuIiONMfnAGGAOkAJMMcas912pSgUHR5iw6fmhnNu0Nn/73zoOZumE4qHA27MRbjTGNDHGRBhjEo0xH7raZxlj2hljWhtjXvBtqUoFj6hwB29c14UTpwp4/Ou1GL1sN+jp6LxSFmnbqBYPDmnHnPX7mbthv9XlKB/TsFXKQiMvaklyo1o8O2ODcxpGFbQ0bJWyUIQjjKeHdyL9yEne+yHN6nKUD2nY2p6O5QW7C1rX47LOTfjnwlR2Hz1Z8RNUQNKwVcoGHr+0A4BzknEVlDRs7S66ttUVKD9IqBPDHRe15OvVu1m3W+/gEIw0bO3u8jetrkD5yej+rYmPjeSFmSl6KlgQ0rC1uxrxVleg/CQuOoL7f9eWpWmHWHD6jrwqaGjYKmUjI3o1p2X9WF6atZH8Ar1UO5ho2CplIxGOMB4Z2p4tB7L48pddFT9BBQwNW6Vs5pJOjeiZFM/Ls1I4kJljdTmqmmjYKmUzIsLL155LbkEhr+mpYEFDw1YpG2rVoCa3X9iSqavS9VSwIKFhq5RNjRnYhvgakTw7Y4OeChYENGyVsqm46Aj+Mrgdy7cdZva6fVaXo6pIw9ZOjqXD842srkLZyA3nN6Ndo5q8+G0KeXoqWEDTsLWTZRMgX48+qyLhjjBG9WvNrsMn2X7whNXlqCrQsLWLfWvhp7etrkLZUMv6sQBMWaHn3QYyDVu7OLrT6gqUTXVOdE5GtOuwTr8YyDRs7SIs3OoKlE1FOMK4pFMjft11VM9KCGAatnYR5rC6AmVj/ZMbsi8zh2l66/OApWFrF2ERpds6Xe3/OpQtDUhuCMD9k3/V3m2A0rC1i5LDCK0GwDXvW1OLsp3GtaNpFBcFwCuz9RLeQKRhaxfh0cWXh74EDg+9XRWy5j14MQATftjKzkPZFlejKkvD1i4cJXq2DTtYU4eyrVrREUwe1RuAV+dstLgaVVkatnah43DKC71b1ePegW2YsWYvU3S+24CiYWsXRi/FVN4ZM7ANrRrE8vB/17BhT6bV5SgvadjahlvP9srx1pWhbC8q3MEr13YGYOxXaygs1L+KAoGGrV24/7wkXWRZGSownJ8UzyvXnsua9GN0e26u1eUoL2jY2oZb2ope4KAqdl2PZgAcO5nHyh1HLK5GVUTD1i7cD5Dp1WTKCyLCyicGERvp4NGv1pCbr+P+dqZhaxvas1WVV69mFG/d0I3N+7N4Z0Gq1eWocmjY2kWxnq1OSqO8N6hjI67qmsC7C1PZmpFldTmqDBq2tuEetvqxqMp5/LKOREc4eHraep07wab0p9ou3H9AImtaV4cKSA1qRfHg4HYs3nKQOev3W12O8kDD1jbcwlbnRFBn4ZbeLWjfuBbPzdjAydwCq8tRJWjY2oX+6aeqKNwRxjPDO7H76EneXagHy+xGw9Y2NGxV1fVqVY8ruyYwYVEaOw7pDSLtRMPWLrZ8Z3UFKkg8dmkHIsKEZ6dvsLoU5UbD1g5OZcGPb1ldhQoSjeKiuX9QW+ZvPMDCTQesLke5aNjawTf3WF2BCjK39WlJUr0avDgrhfwCvbLMDjRs7WD3aqsrUEEmMjyMscM6sHl/Fv9eofPe2oGGrVJB6pJOjeiZFM8b323meE6e1eWEPA1bO9DzapUPiAhPXN6BQydymfDDVqvLCXnVHrYi0kFEJojIVBEZ7WoLE5EXRORtEfljdb9nwNOwVT7SObEOV3ZN4MMl29h3LMfqckKaV2ErIh+JyAERWVeifaiIbBKRVBEZC2CMSTHG3ANcB/RwbXol0BTIA9Krr/wgEaZhq3znoSHJFBQa3py32epSQpq3PduPgaHuDSLiAMYDw4COwI0i0tG1bjiwBJjv2jwZWGqMeRAYXfWyg0zJO+sqVY2axdfg5t4tmLJiF1v2H7e6nJDlVdgaYxYBh0s09wRSjTFpxphcYDLOHizGmGnGmD7ATa5t04HTU8l7vGhbREaJyAoRWZGRkVHJ3Qhga6fC3t+srkIFuXsHtqVGZDivzN5kdSkhqypjtk0B93NK0oGmItJfRMaJyHvALNe6r4BLRORtYJGnFzPGTDTG9DDG9GjQoEEVygow/x1pdQUqBMTHRnLPxa2Yl7KfX7aX7Dcpf6jK36/ioc0YYxYCC0s0ZgOaKkpZaORFrfjs5x28OCuFr0b3QcTTj7Dylar0bNOBZm7LicCeqpWjlPKVmEgHfxnUjtU7jzJn/T6rywk5VQnbX4C2ItJSRCKBG4Bp1VOWUsoXfn9eIm0a1uTV2ZvI08t4/crbU7++BJYCySKSLiIjjTH5wBhgDpACTDHGrPddqUqpqgp3hPHI0PakHTzB5F/0Ml5/8mrM1hhzYxntsyg6CKaUCgCDOjTk/KS6vDVvC9d0a0pslJ566A96ua5SIUZEePTSDhzMOsX7i9OsLidkaNhaydOtcPo97P86VMjp3rwuw85pzMRFaWQcP2V1OSFBw9ZKhfml2zpe6f86VEj6v0uSOZVfyLj5W6wuJSRo2FqpwMO0d3ruo/KTVg1qcmPPZkxavpO0jCyrywl6GrZW2jzbQ6OGrfKf+3/XjqjwMF6bo5fx+pqGrb/lZMKOpc7H25eUXi/6kSj/aVArilH9WvHtun2s2nmk4ieos6Y/2f72nz/Cv4bC4TRY8WHp9TqMoPzsrr6tqF8zipdnbcR4OmirqoWGrb+dnuFrxb/K2EDDVvlXbFQ49w9qy/Lth5mXonfj9RUNW3873XNIme55vfZslQVuOL8ZrerH8srsjXo3Xh/RsLXKkW1lrNCwVf4X4Qjj4aHJpB7IYupKvZmKL2jY2o32bJVFLunUmO7N6/CPeZvJzvVwDriqEg1bpRRQdBnv/sxTfLSkrL+81NnSsFVKnXF+UjyDOzZiwg9pHMrSy3irk4at3+mpNcreHhmaTHZuPm9/n2p1KUFFw9Zu9DxHZbE2DWtx/fnN+GLZDnYcOmF1OUFDw9ZujJ52o6z3wKB2OMJEL+OtRhq2/lZhz1V7tsp6jeKiuatvK2as2ctvu45aXU5Q0LD1t5yjxZcHPF58WXu2yiZG9WtFfGwkL3+rl/FWBw1bfzp1vHTbxSUmC9dvamUTtaIjuG9gG5amHWLh5gyrywl4Grb+dGhr8eXYBqW30Z6tspERvVrQol4NXp61kYJC7QhUhYatPx3dWXw5obuHjfQbWtlHZHgYDw1JZtP+43y1Si/jrQoNW386uqP4cpjD+fXpY9DoXOdj7dkqm7ns3CZ0SazNG3M3k5NXYHU5AUvD1p9K9mzdJwq/9n3ofD007OTfmpSqQFiYMHZYB/Yey+Hjn7ZbXU7A0rD1pyMleraH3a4/b9gBrpkIjnD/1qSUFy5oXY8ByQ0YvyCVIydyrS4nIGnY+lPJnu2BDdbUodRZeGRYe7JO5TN+gV7GezY0bP3FmNJjtn3/ak0tSp2F9o3j+H33RD5duoNdh7OtLifgaNj6y4kMyCvxDRrfyppalDpLDw5phwi8MXez1aUEHA1bfyk5XnvLN9B1hCWlKHW2mtSO4fYLW/LNr7tZt/uY1eUEFA1bfzmyvfhy6wF6VwYVkEb3b03tmAhemb3R6lICioatvxzd7vx69yK4d5WlpShVFbVjIhgzoA2Ltxxk8Ra9jNdbGrb+cmQ71GwETbpAvdZWV6NUldxyQQsS68bw8rcbKdTLeL2iYesvR3ZAnRZWV6FUtYgKd/DQkGTW78lk2m97rC4nIGjY+suRHVA3yeoqlKo2w7sk0CkhjtfmbNLLeL2gYesPBXmQma5hq4KK8zLe9uw+epLPf95R8RNCnIatPxzb5Zxgpq4OI6jg0rdtAy5qU5/xC1LJzMmzuhxb07D1h9OnfWnPVgWhR4a250h2Hu8vSrO6FFvTsPWH0xc06AEyFYTOTazN5Z2b8MHibRzIzLG6HNvSsPWHI9shLALiEqyuRCmfeGhIMnkFhYz7fovVpdiWhq0/bJ4DtZsWTRauVJBJqh/LDT2bMXn5LrYdPGF1ObakYetrBzZCRkrpy3WVCjL3/a4tEY4wXv9uk9Wl2JKGra8d0rk/VWhoWCuaO/u2ZOaavaxN10lqStKw9bXTc9jeOd/aOpTyg1H9WlG3hk5S44mGra8d2grRdSCxh9WVKOVztaIjGDOwLUtSdZKakjRsfe1QKtRrY3UVSvnNzb2b07RODK/M1klq3FV72IpIBxGZICJTRWS0W3usiKwUkcur+z1t7XCazvKlQkpUuIMHB7dj3e5MZq7da3U5tuFV2IrIRyJyQETWlWgfKiKbRCRVRMYCGGNSjDH3ANcB7n87PwJMqa7CA0JejvNS3XgNWxVarurWlPaNa/H6d5vIKyi0uhxb8LZn+zEw1L1BRBzAeGAY0BG4UUQ6utYNB5YA813Lg4ANwP5qqTpQHHdNPVc70do6lPIzR5jw8NBkdhzKZvLynRU/IQR4FbbGmEXA4RLNPYFUY0yaMSYXmAxc6dp+mjGmD3CTa9sBQG9gBHCXiITGWHGm60+ouCbW1qGUBQYkN6RnUjxvzU/lxKl8q8uxXFVCrymwy205HWgqIv1FZJyIvAfMAjDGPG6MeQCYBLxvjCn1d4WIjBKRFSKyIiMjSI5iHneFbS29TFeFHhHhkWHtOZh1io+WbLO6HMuFV+G5nu5WaIwxC4GFnp5gjPm4rBczxkwEJgL06NEjOA5hZrqGEbRnq0LUeS3qMqRjI95blMaIXs2pVzPK6pIsU5WebTrQzG05EdD7Y7g7vhciYiEqzupKlLLMw0OTyc7NZ/yCrVaXYqmqhO0vQFsRaSkikcANwLTqKStIZO5x9mr1luUqhLVpWIvfn5fI5z/vIP1IttXlWMbbU7++BJYCySKSLiIjjTH5wBhgDpACTDHGrPddqQHo+D6opUMISj0wqB0i8MbczVaXYhmvxmyNMTeW0T4L10Ew5cHxPdCst9VVKGW5hDox3NYniYmL07irbys6NAm9obXQOAXLCjmZcGy3ThiulMvo/q2pFRXOa3NCcwpGDVtf2bcGTAE0v8DqSpSyhTo1Ihndvw3fbzzAsrRDVpfjdxq2vnL6tK/4VtbWoZSN3NYniUZxUbw8eyPGBMcZnt7SsPWVHNfkyTF1ra1DKRuJiXTwwKB2rN55lO82hNbV+xq2vpJ30vk1IsbaOpSymT+cl0irBrG8NmcT+SE0SY2Gra9o2CrlUbgjjIcvSSb1QBb/XZVudTl+o2HrK3nZ4IjUO+oq5cElnRrTtVkd/jF3Czl5BVaX4xcatr6SdxLCtVerlCciwiND27MvM4dPftpudTl+oWHrKwW5EB5pdRVK2dYFrevRP7kB4xekciw7z+pyfE7D1lcK8yGsKpOqKRX8Hr6kPcdP5fPuD8E/SY2Gra+YQg1bpSrQMSGOK7sk8K8ft7HvWI7V5fiUhq2vFOZDiNyQQqmq+OuQZAqN4a35wT1JjaaBrxQWaM9WKS80i6/BTb1a8O9fdpF6IMvqcnxGw9ZXCvP1tC+lvDRmYBtiIhy8HsST1GjY+orRnq1S3qpfM4pR/Voze/0+Vu08YnU5PqFh6yuFBdqzVaoS7uzbkvo1I3nl2+CcpEbD1lcK80E0bJXyVmxUOPcObMuybYdZuDlI7rDtRsPWV/QAmVKVdmPP5jSPr8Er326ksDC4ercatr6iB8iUqrTI8DD+OqQdG/cd53+/7ba6nGqlYesrelGDUmflis4JdEqI4+/fbeZUfvBMUqNh6yt6UYNSZyUszDlJTfqRk0xattPqcqqN2PGon4hkADtci7WBY148raLtPK33ps19uazH9YGDXtRYHm/3s6Jty1pX3n5VtHz6cXXsZ3k1VmY7b/fTU5vdPlP93i1/XSB977YwxjTwuMYYY+t/wMTq2M7Tem/a3JfLebzCX/tZ0bZlrStvv7zd7+rYz+r6TL3dz0D4TPV7N/i+dz39C4S/c6dX03ae1nvTNt2Lx9WhMq9X3rZlrStvvypatmpfq2M/PbXZ7TPV793y1wXi924pthxGCDQissIY08PqOnwtVPYTQmdfdT/9JxB6toFgotUF+Emo7CeEzr7qfvqJ9myVUsoPtGerlFJ+oGGrlFJ+oGGrlFJ+oGHrAyJylYi8LyL/E5EhVtfjKyLSQUQmiMhUERltdT2+JCKxIrJSRC63uhZfEpH+IrLY9bn2t7oeXxGRMBF5QUTeFpE/+uM9NWy9JCIficgBEVlXon2oiGwSkVQRGQtgjPnGGHMXcBtwvQXlnrVK7meKMeYe4DogoE4fqsx+ujwCTPFvldWjkvtqgCwgGkj3d61VUcn9vBJoCuThr/2sjqsqQuEf0A/oDqxza3MAW4FWQCTwG9DRbf3fge5W1+7L/QSGAz8BI6yu3Vf7CQwCbsD5y/Nyq2v38b6GudY3Ar6wunYf7udY4G7XNlP9UZ/2bL1kjFkEHC7R3BNINcakGWNygcnAleL0CvCtMWaVv2utisrsp2v7acaYPsBN/q20aiq5nwOA3sAI4C6RwJphqDL7aowpdK0/AkT5scwqq+Rnmo5zHwH8MrWYzgFYNU2BXW7L6UAv4F6cvaHaItLGGDPBiuKqkcf9dI3pXYPzh3KW/8uqdh730xgzBkBEbgMOugVSICvrM70GuASoA7xjQV3Vrayf0beAt0WkL7DIH4Vo2FaNeGgzxphxwDh/F+NDZe3nQmChf0vxKY/7eeaBMR/7rxSfK+sz/Qr4yt/F+FBZ+5kNjPRnIQH155ANpQPN3JYTgT0W1eJLup/BJ1T21Tb7qWFbNb8AbUWkpYhE4jyIMs3imnxB9zP4hMq+2mY/NWy9JCJfAkuBZBFJF5GRxph8YAwwB0gBphhj1ltZZ1XpfgbXfkLo7Kvd91MnolFKKT/Qnq1SSvmBhq1SSvmBhq1SSvmBhq1SSvmBhq1SSvmBhq1SSvmBhq1SSvmBhq1SSvmBhq1SSvnB/wOQzmLlC0AeOAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pl.loglog(s99.wavelengths, s99.data[:,-1])\n", "pl.loglog(bpass.wavelengths, bpass.data[:,-1])\n", "pl.ylim(1e33, 1e41)" ] }, { "cell_type": "markdown", "id": "092cb425", "metadata": {}, "source": [ "The most common options for these models are: ``pop_Z``, ``pop_ssp``, ``pop_binaries``, ``pop_imf``, and ``pop_nebular``. See :doc:`params_populations` for a description of each of these parameters.\n", "\n" ] }, { "cell_type": "markdown", "id": "2c6bb59e", "metadata": {}, "source": [ "## Parametric SEDs for galaxies and quasars" ] }, { "cell_type": "markdown", "id": "eea60e34", "metadata": {}, "source": [ "So far, there is only one litdata module in this category: the multi-wavelength AGN template described in Sazonov et al. 2004." ] }, { "cell_type": "markdown", "id": "11871951", "metadata": {}, "source": [ "## Reproducing Models from *ARES* Papers" ] }, { "cell_type": "markdown", "id": "dee1ba07", "metadata": {}, "source": [ "If you're interested in reproducing a model from a paper exactly, you can either (1) contact me directly for the model of interest, or preferably (someday) download it from my website, or (2) re-compute it yourself. In the latter case, you just need to make sure you supply the required parameters. To facilitate this, I store \"parameter files\" (just dictionaries) in the litdata framework as well. You can access them like any other dataset from the literature, e.g., " ] }, { "cell_type": "code", "execution_count": 7, "id": "a52e2268", "metadata": {}, "outputs": [], "source": [ "m17 = ares.util.read_lit('mirocha2017')" ] }, { "cell_type": "markdown", "id": "cfeb474c", "metadata": {}, "source": [ "A few of the models we focused on most get their own dictionary, for example our reference double power law model for the star-formation efficiency is stored in the ``dpl`` variable:" ] }, { "cell_type": "code", "execution_count": 8, "id": "4a8cc5be", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "# Loaded $ARES/input/inits/inits_planck_TTTEEE_lowl_lowE_best.txt.\n", "\n", "############################################################################\n", "## ARES Simulation: Overview ##\n", "############################################################################\n", "## ---------------------------------------------------------------------- ##\n", "## Source Populations ##\n", "## ---------------------------------------------------------------------- ##\n", "## sfrd sed radio O/IR Lya LW LyC Xray RTE ##\n", "## pop #0 : sfe-func yes - - x x x - x ##\n", "## pop #1 : sfrd->0 yes - - - - - x x ##\n", "## ---------------------------------------------------------------------- ##\n", "## Notes ##\n", "## ---------------------------------------------------------------------- ##\n", "## + pop_calib_lum != None, which means changes to pop_Z ##\n", "## will *not* affect UVLF. Set pop_calib_lum=None to restore ##\n", "## \"normal\" behavior (see S3.4 in Mirocha et al. 2017). ##\n", "## ---------------------------------------------------------------------- ##\n", "## Physics ##\n", "## ---------------------------------------------------------------------- ##\n", "## cgm_initial_temperature : 20000 ##\n", "## clumping_factor : 3 ##\n", "## secondary_ionization : 3 ##\n", "## approx_Salpha : 3 ##\n", "## include_He : True ##\n", "## feedback_LW : False ##\n", "############################################################################\n", "# Loaded $ARES/input/bpass_v1/SEDS/sed.bpass.constant.nocont.sin.z020\n", "# WARNING: revisit scalability wrt fesc.\n", "# Loaded $ARES/input/hmf/hmf_ST_planck_TTTEEE_lowl_lowE_best_logM_1400_4-18_z_1201_0-60.hdf5.\n", "# Loaded /Users/jordanmirocha/Dropbox/work/mods/ares/input/optical_depth/optical_depth_planck_TTTEEE_lowl_lowE_best_He_1000x2158_z_5-60_logE_2.3-4.5.hdf5.\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "gs-21cm: 100% |#############################################| Time: 0:00:03 \n" ] }, { "data": { "text/plain": [ "(,\n", " )" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEsCAYAAADTvkjJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAy60lEQVR4nO3dd3xUdb7/8dcnPSSBAEkooYUaqjQRFFREBXt3dV0s6CKu7uru77rquveu2667etdddW3YsCsWrg2kiqhIVYTQQiehhBJ6Qkgyn98fc7ibxSSkzMyZOfN5Ph7zYHLOmTnvHE4+OfnO9/s9oqoYY4yJfDFuBzDGGBMYVtCNMcYjrKAbY4xHWEE3xhiPsIJujDEeYQXdGGM8wgq6McZ4hBV0Y4zxCCvoxhjjEVbQjQkQEblLRLSax1S3s5noIDb035jAEJFUILXKoquBvwKXq+pMd1KZaGIF3ZggEJGfAE8BV6rqbLfzmOhgTS7GBJiI3AY8CVxkxdyEkhV0YwJIRO7C38wyWlW/cjuPiS5xbgcwxitE5D+A+4BzVfU7t/OY6GNt6MYEgIg8ANwPXAmsrLLqsKoedieViTZW0I1pJBERYD/QtJrV96jq46FNZKKVFXRjjPEI+1DUGGM8wgq6McZ4hBV0Y4zxCCvoASYivxSRlSKSJyJviUiSiLQQkZkiss75t7lL2ZJEZJGIfO9k/L2zPCT5RKS9iHwuIqud/d8dqv2LyEsisktE8qose0hEtonIMudxYQj2+aiIrBGR5SIyRUTSQ7DPsDj/alLdz4zbmUJNRNJF5D3n3FgtIsMa8j5W0ANIRLKBXwCDVbUPEAtch78722xV7QbMdr52QxlwjqqeAvQHxojI0BDmqwD+n6r2BIYCd4pIrxDtfxIwpprlf1fV/s4j0JNoVbfPmUAfVe0H5AMPhGCf4XL+/UAtPzPR5nHgM1XNBU4BVjfkTaygB14ckCwicUATYDtwGfCKs/4V4HI3gqnf8T7R8c5DCVE+Vd2hqt86zw/hP2mzQ7F/VZ0HFAf6feu7T1WdoaoVzpcLgHbB3idhcv7VorqfmaghIk2BM4EXAVT1mKrub8h7WUEPIFXdBvwPsBXYARxQ1RlAK1Xd4WyzA8hyK6OIxIrIMmAXMFNVF7qRT0Q6AQMAV/ZfxV1O88dLLjRFjAOmhWA/YXP+naiWn5lo0hnYDbwsIt+JyAsiktKQN7KCHkBOQbgMyAHaAinOrHthQ1UrVbU//ivDISLSJ9QZnGlm38c/6OZgqPdfxTNAF/zNTzuAv4VqxyLyIP4mqDdCtc9wFAk/MyEQBwwEnlHVAcARGtgsZgU9sM4FNqnqblUtBz4ATgeKRKQNgPPvLhczAuD8STcXf3tryPKJSDz+Yv6Gqn7gLHbl+KhqkfMLzgc8DwwJxX5F5CbgYuAGDc3IvrA7/6qo6WcmmhQChc5fywDv4S/w9WYFPbC2AkNFpIkzHHwU/nbij4CbnG1uAj50I5yIZB7vVSEiyfh/mNaEKp9zTF4EVqvqY1VWuXJ8jhc5xxVAXk3bBnCfY/BP4HWpqpYEe3+OsDj/alDTz0zUUNWdQIGI9HAWjQJWNeS9bOh/gDldAX+E/8/p74Db8N/FZjLQAf8JfI2qhvQDOidbP/wfisXi/2U+WVX/ICItQ5FPRIYDXwIrAJ+z+Df429GDun8ReQs4G8gAioDfOV/3x//B8Gbg9uNtzUHc5wNAIrDX2WyBqk4I8j7/lzA4/2pS3c+Mqpa5myq0RKQ/8AKQAGwEblHVffV+HyvoxhjjDdbkYowxHmEF3RhjPMIKujHGeIQVdGOM8YiILOgisllEVjgTKi1xloXlBEQiMt7tDCfjZsZo23e07LO+IiFjsAXiGERkQXeMdCZUGux8Ha4TEDX4P0lELgnkdrVw84cp2vYdLfusr0jIGGxRXdBPFO4TEDVEXQt1Ywu6McYDIrIfuohsAvbhHxDynKpOFJH9qppeZZt9qvqDZhfnz5rjvwkHNWnSJKhZKyoqiIuLa9BrKysriY2NDdh2NWlMxsaKtn1Hyz7rKxIyBltdj0FJSYmqarUX45Fa0Nuq6nYRycI/v/TPgY/qUtCrSklJ0SNHjgQ3rDHGBJCIlKhqtbMxRmSTi6pud/7dBUzBP6lSOE9AZIwxQRdxBV1EUkQk7fhz4Hz8kyqF8wRExhgTdJHYaNUKmOKfmI044E1V/UxEFgOTReRWnAmIXMxojDEhF5Ft6IFibejGmEjjuTZ0Y4wxP2QF3RhjPMIKujHGeIQVdGOM8Qgr6MYY4xFW0I0xxiOsoBtjjEdYQTfGGI+wgm6MMR5hBd0YYzzCCroxxniEFXRjjPEIK+jGGOMRVtCNMcYjrKAbY4xHWEE3xhiPsIJujDEeYQXdGGM8wgq6McZ4hBV0Y4zxCCvoxhjjEVbQjTHGI6ygG2OMR1hBN8YYj7CCbowxHmEF3RhjPMIKujHGeISnCrqIjBGRtSKyXkTudzuPMcaEkmcKuojEAk8BFwC9gOtFpJe7qYwxJnTi3A4QQEOA9aq6EUBE3gYuA1bV9AKfKjsOlNIsOZ4mCV46FMaY6uw9XMaRssqTbqfoybc5+SbOe9Vxu7q+YS28VMWygYIqXxcCp9X2grJyH8MenkOMwB1nd+E/zu+BiAQ1pDEmtFSVuWt3M2n+Zr7I3+12nKDyUkGvrhL/4FeeiIwHxgPExSfw8JV9+Xr9Hp76fAPpyQn89MzOwc4ZEVSVfSXlHCwt59DRCkrLK0mMiyEtKY7s5skkxsW6HdGYWqkqs1bv4h+z8lm5/SBZaYn84pyudGyZUqfX1/Xari7bSbXlqWHvdflfa17npYJeCLSv8nU7YPuJG6nqRGAiQEpKil4/pAM/Gtye8kofj0xfw4juGeS2bhqaxGFAVdlaXELetoOs2HaAdUWH2FpcQuG+UkrLq//TVATaNE2if4d0TstpyZg+rWnVNCnEyY2p2cKNe/njp6vI23aQji2b8OjV/bh8QDbxsZ752LBaEoh2m3AgInFAPjAK2AYsBn6sqitrek1KSooeOXIE8Letnff3eXTNSuWd8UM93fSy40Ap8/J3My9/D19v2MP+knIA4mOFLpmpdGjRhA4tmpDdPJn0JvGkJsaTHB9LWUUlB0rLKSguZeOewyzZvI9t+0uJETizeyZ3juzKqZ1auPzdmWi2fX8pD09bw8ffbyc7PZl7zu3GFQOyifNQIReRElWt9s8MzxR0ABG5EPgHEAu8pKp/rm37qgUd4NVvNvNfH67k5ZtPZWRuVlCzhlrRwaN8/P12Pv5+O98XHgCgVdNERnTLZFDH5vTNbka3Vqn1bkrZsPswH363jTcXFbDncBnn9mzFny7vQ+tmdsVuQsfnU15fuIWHp67Bp8qEs7ow4awuJCd4r2kwagp6fZ1Y0MsrfZz32BckJ8Qx9RfDI/4q3edTvsjfzSvf+D8MUoU+2U25uF9bRvbIonur1IB9j6XHKpk0fzOPz84nPjaGR68+hTF9WgfkvY2pTUFxCfe9v5z5G/ZyZvdM/nx5H9q3aOJ2rKCxgl6DEws6wLtLCrj3veVMuuVUzu4RmVfpxyp8TF5SwMR5G9laXEJWWiLXDenA5f3b0jkzNaj73rznCHe/s4zvC/Zz35hc7ji7S1D3Z6Lb1BU7+PV7y1FVfntxL647tX3EX4idjBX0GlRX0I9V+Djr0c/p2LIJb48f5lKyhqn0Kf/73Tb+MTufguJSBnRI59bhOYzu3TqkHwYdLa/k3veW8/H32/l/53Xn56O6hWzfJjqUV/r4y7Q1vPjVJvq3T+fJ6wd4+qq8qtoKupd6uQREQlwMtw7P4U+frua7rfsY0KG525Hq5Iv83fzxk1Ws33WYPtlN+cMtfTi7e6YrVytJ8bH840f9iY8R/jYzn7SkOG4+IyfkOYw37TlcxoTXlrJkyz5uGtaRBy/qRUKcdz70bAy7Qj/hCh3gSFkFp/9lDsM6t+TZsYNcSFZ32/aX8sePV/HZyp3kZKTw69E9GNOndVj82VnpU+54fSmz1+zitXFDOL1rhtuRTIRbv+swt0xaxO5DZfz1qn5c1j/b7UghV9sVuv1aq0ZKYhxjh3Zk+qqdbNrzw4IfLj76fjuj/z6PL/J3c+/oHnx2zwgu6NsmLIo5QGyM8NiP+tM5I4U73/yWHQdK3Y5kItjCjXu56pn5lJRV8tZPh0ZlMT8ZK+g1uOn0TsTHxvDClxvdjvIDR8sreeCDFfzire/o0TqNGb88kztHdg3L0ZupiXE8N3YQR8t93Pf+ioDMV2Giz8xVRYx9cREZqQlM+dkZEdMUGmpW0GuQmZbIVQOzeXdpIXsOl7kd5/+s33WIy/75NW8t2srPzu7C2+OHhv2HQZ0zU7n/glzm5e9m8pKCk7/AmCqmrtjBHa8vpWebNN6/43Q6tAzv891NVtBrcduIzpRX+nh1/ma3owDw4bJtXPLk1+w5XMYr44bw6zG5ETOUeezQjgzt3II/fbKa3YfC5xekCW8fLtvGz9/6jv7t03n9ttNIb5LgdqSwFhnVwCVdMlM5t2crXl2whZJjFa7l8PmURz5bw91vL6Nvu2ZMvXsEZ3XPdC1PQ8TECH+6vC+l5ZU8NjPf7TgmAnyyfDv3vLOMUzs155VxQ0hLinc7Utizgn4St5/Zmf0l5by7pNCV/Zceq+T215fy9NwNXD+kPa/felrEToTVNSuVscM68s7irazaftDtOCaMfZG/m1++s4zBHZvz8s1DSEm0HtZ1YQX9JAZ3asHADum88NVGKip9Id33waPl3PTSImatLuK/Lu7Ff1/RN+L7294zqjtpSfE8NnOt21FMmFq6ZR8TXltK16w0XrjpVE/OxxIskV0dQmT8mV0oKC7ls5U7Q7bP4iPHuOH5hXy7dR9PXDeAccNzwqY7YmM0axLPbcNzmLV6F3nbDrgdx4SZ9bsOM27SYlo1TeTVcUNolmzNLPVhBb0OzuvVipyMFJ77YmNIut0VHTzKj577hvyiQ0y8cRCXnNI26PsMpZvO6ETTpDiemL3O7SgmjOw7coxbX1lMfKzw2q2nkZmW6HakiGMFvQ5iY4Q7zurCim0HmL6yKKj7Kigu4Zpnv2H7/lIm3TKEc3JbBXV/bmiaFM+44TnMWFVkbekG8M+hNOH1pew4cJTnxg4O+6644coKeh1dOTCbLpkpPDp9TdDa0tfvOsTVz87nQGk5b/x0KMO6tAzKfsLBLWfkkJIQywtfhd/ALRNaqsrvPspj4aZiHrmqH4M62qChhrKCXkdxsTHcOzqXDbuP8N7SwPd4ydt2gGufW0ClD965fSj926cHfB/hpFlyPFcPascn3++wfulR7vUFW3hrUQF3jezK5QNsOH9jWEGvh9G9WzGgQzqPzczn4NHygL3vks3FXD9xAcnxsbw7YVjU3NP0xtM7cazSx5sLt7odxbhkeeF+/vjJakb2yORX53V3O07Es4JeDyLC7y/tzZ7DZTw8dU1A3nP2av8cFZlpibw7YRg5GXW7I7kXdMlM5ewemby+cAvHKkLbJdS470BJOT9741sy0xJ57Nr+xMREfi8ut1lBr6d+7dK5bURn3lq0lfnr9zTqvd5YuIWfvrrEf2Pq24fRNj05QCkjx82nd2L3oTKm5e1wO4oJIZ9P+X/vLqPo4FGeumEgzVNsSH8gWEFvgF+e253OmSn84m3/CVlfqsr/TF/Lg1PyOKt7Jm+PHxq1XbTO7JZJ+xbJNmlXlHnp603MWr2LBy/s6fnPi0LJCnoDJCfE8uxPBlFyrII7Xl9K6bHKOr/24NFy7nrzO/75+XquO7U9z984OKqHNcfECFcPbM/X6/dSUFzidhwTAmt2HuSRz9Zybs9W3HR6J7fjeIoV9Abq3iqN/7nmFL4r2M+4SYs5UnbyybuWFeznoie+5LOVO7n/glwevrIvcREyW2IwXTUoGxGC0nvIhJeyikrueXsZTZPj+MtVfT0x+jmcWDVphAv7tuHv1/Zn4aa9XPD4l0xfuROf799Hkqoqq7Yf5Hcf5nHVM/Px+WDy7cOYcFYXO5kd7Zo3YXjXDN5bWviD42e85bGZ+azZeYi/XtWPjNTobGYMpuj9Wz9ALh+QTZtmSdz/wQpuf20pGakJ9G/fnKbJcew9fIz8okPsOHCUGIEfn9aBe8/PpVkTm5/iRNcMbs8v3vqObzbu5Qy796gnLdy4l4nzNnL9kA6M6um9EdDhwG4SXc1NohuiotLHZyt3MmtVEat3HOJwWQXpTeLJyUhhRLcMRuZmkZUWmdPehsLR8kpO/fMszuvViseu7e92HBNgR8srGf2PeajCtLtHRPXnRo1V202i7agGSFxsDBf3a8vF/bw1kVaoJMXHMqZ3az7L28nR8kqS4m3KVC95fPY6tuwt4c3bTrNiHkTWhm7CxiWntOVQWQVz1+52O4oJoJXbDzBx3kauHdyO0605LaisoJuwcXqXlrRMSeDj5dvdjmICpKLSx/3vr6B5kwR+c2FPt+N4XkQVdBF5SES2icgy53FhlXUPiMh6EVkrIqPdzGkaJi42hgv7tmH26qI6dQM14e/lrzezYtsBHrq0l93gOQQiqqA7/q6q/Z3HVAAR6QVcB/QGxgBPi4g1wkagS05py9FyH7NWB3feeRN82/eX8tjMfEblZnFR3zZux4kKkVjQq3MZ8LaqlqnqJmA9MMTlTKYBBndsTuumSXyy3OZ2iXR/nroanyoPXdrbxlyESCQW9LtEZLmIvCQix2fCzwaqTgZS6Cz7AREZLyJLRGRJRYX9WR9uYmKE0b1bMS9/NyXH7P8nUs3fsIdPl+/gjrO72N2HQijsCrqIzBKRvGoelwHPAF2A/sAO4G/HX1bNW1XbwV5VJ6rqYFUdHBdn3afC0ejerSmr8DEv33q7RKKKSh+//2gV2enJTDiri9txokrYVTRVPbcu24nI88AnzpeFQPsqq9sB1lUiQg3JaUGz5HimryxiTB9re400ry/YwtqiQzz7k4E2niDEwu4KvTYiUvWn+wogz3n+EXCdiCSKSA7QDVgU6nwmMOJiYxjVM4vZq4soD9L9W01w7D1cxmMz8xneNYPRvVu7HSfqRFRBBx4RkRUishwYCfwSQFVXApOBVcBnwJ2qWvc5bU3YGd27NQePVrBwY7HbUUw9PD57HUeOVfLQpb3sg1AXhF2TS21UdWwt6/4M/DmEcUwQndktk6T4GKav3Mnwbja6MBJs3H2YNxdu5foh7emaleZ2nKhU7yt0EekhIheIyJUiMkJEUoMRzES35IRYzuqeyYxVP5yS2ISnR6evJSEuhrtH2c2e3VKngi4inUTkERHZjr9Z41PgPeALoFhE5ojItWJ/Y5kAGt27NUUHy1i+7YDbUcxJLN2yj2l5O7n9zC5RezvFcHDSgi4ij+L/8LEH8BugD9AMSATaABcC84G/AMtEZGDQ0pqock5uFrExwoyVO92OYmqhqjw8dTWZaYncNiLH7ThRrS5X6GlAd1W9TFUnqepqVT2kquWqWqSqs1T1t6raGfgTYDPwmIBIb5LAoI7NmbNml9tRTC1mrCpiyZZ9/PLc7jY1rstOWtBVdYKq1qlPt6q+q6pvND6WMX6jcrNYs/MQ2/eXuh3FVKOi0sdfP1tDl8wUrh3czu04Ua+ubegnbRQTkf6NTmPMCc7JzQLg87V2lR6Opny3jY27j3Dv6Fy74XkYqOv/wJu1rXRmO5ze+DjG/LuuWam0a57M59bsEnbKK308MWcdfbKbMrq33SM0HNS1oJ8uIk9Wt0JEugKzsJGZJghEhHNys/h6/V6OlttYsXDy/tJCCopL+dV53W0QUZioa0G/CLhJRB6oulBEOgCzgdXAVQHOZgwAI3OzKC2vZMHGvW5HMY5jFT6enLOeU9qnM7JHlttxjKNOBV1VvwWuBn4nImPh/+ZVmQ1sAy5R1WNBS2mi2rDOLUmKj7FmlzAyeUkB2/bb1Xm4qfOnGKo6A/gp8LyI3IC/meUgcIGqlgQpnzEkxccyvGsGc9buQtVGjbrtaHkl/5yznsEdm3OmTcsQVur1sbSqvgb8J/Aq/vnGz1NVG8Zngm5kbhYFxaWs33XY7ShR7+1FW9l58KhdnYehOo0CEJEZJywqdx5vV/0PVdXzAxfNmH853k47Z80uurWyiZ/ccrS8kqfnbuC0nBYM69LS7TjmBHUd1rXthK/fCnQQY2rTNj2Z3NZpzFmzi9vtLjiueXdpIbsOlfGP6/rb1XkYqlNBV9Vbgh3EmJM5JzeL5+Zt5EBpOc2S492OE3XKK30898UGBnZIZ1hnuzoPRza0y0SMc3KzqPQpX66ze4264ePvt1O4r5Q7R3a1q/MwVe+ZdERkNDAKyOKEXwiqemOAchnzAwM6NCe9STyfr9nNxf3auh0nqvh8ytNzN5DbOu3/pmMw4adeBV1E/oR/Ct3lwE78PV2MCYnYGOHMbpl8kb8Ln0+JibGrxFCZsaqI9bsO88T1A+zqPIzV9wp9PHCzqr4ajDDGnMw5uVl89P128rYfoF+7dLfjRAVV5em56+nUsgkX9W1z8hcY19S3Dd2H/2YWxrjizO6ZiGBzpIfQV+v3sLzwABPO6kKs/VUU1upb0J8GbgtGEGPqokVKAv3bp/P5WvtgNFSe+nw9rZsmccXAbLejmJOob5PLH4FPROR7/O3o5VVXquq4QAUzpiYje2Tx91n57DlcRkaq3b8ymJZuKWbBxmL+8+JeJMbFuh3HnER9r9D/AFwAxOK/n2j7Ex7GBN05uVmowrx8u0oPtme/2EjzJvFcP8R+vCNBfa/Q7wLGqeqkIGQxpk56tWlKZloic9bs4sqBdtuzYNm85wizVhdx18iuNEmwe4VGgvpeoR8DvgpGEGPqKiZGOLt7JvPyd1NR6XM7jmdNmr+ZuBhh7NCObkcxdVTfgj4RuDUYQYypj5G5WRw8WsF3BfvdjuJJB0rLmbykgEtOaUtW0yS345g6qu/fUW2Aq5zRot/zww9FxwcqmDG1Gd4tg7gYYc6aXZzaqYXbcTzn7UVbKTlWya3Dc9yOYuqhvlfoXYBlwAGgE9CtyqNroEKJyDUislJEfCIy+IR1D4jIehFZ6/xiOb58kIiscNY9ITaczdOaJsUzuFNzu4tREJRX+pg0fzPDOrekd9tmbscx9VCvK3RVHRmsICfIA64Enqu6UER6AdcBvYG2wCwR6a6qlcAz+EeyLgCmAmOAaSHKa1wwskcWD09bw44DpbRplux2HM+YlreTHQeO8sfL+rgdxdTTSa/QRWRQXd9MRJJEpGfjIoGqrlbVtdWsugx4W1XLVHUTsB4Y4tzftKmqfqP+e5S9Clze2BwmvI10Jomaa4OMAkZVefHLjeRkpNgkXBGoLk0uH4rIFBEZLSLVbi8i2SLyALAOOCOgCf9dNlBQ5etCZ1m28/zE5T8gIuNFZImILKmoqAhaUBN83bJSyU5PtmkAAmjpln18X3iAcWd0ssnPIlBdmlx6APcDrwNJIvId/jsYHQVa4G/+yAHmAterap26NYrILKB1NaseVNUPa3pZNcu0luU/XKg6EX9vHVJSUmy2yAgmIozMzeSDb7dRVlFpIxkD4MWvNtEsOZ6rBln//kh00it0VT2iqv8JtAPGAkuAJPw9Xg4CTwG9VXVUXYu5877nqmqfah41FXPwX3lXHbLWDtjuLG9XzXLjcSN7ZFFyrJLFm/a5HSXiFRSXMH3lTn58WgcbSBSh6vy/pqplwP86D7d8BLwpIo/h/1C0G7BIVStF5JCIDAUWAjcCT7qY04TIsC4tSYiLYc6aXQzvluF2nIj28tebiRHhpmGd3I5iGigsb0EnIleISCEwDPhURKYDqOpKYDKwCvgMuNPp4QJwB/AC/g9KN2A9XKJCk4Q4hnVuydy11o7eGAeP+gcSXdyvDa2b2UCiSBWWf1ep6hRgSg3r/gz8uZrlSwDrZxWFRvbI5KGPV7F5zxE6ZaS4HSciTV5cwOGyCm4d3tntKKYRwvIK3Zj6+Ff3RbtKb4iKSh8vf72ZITkt6NvOBhJFMivoJuJ1bJlC58wU5lh/9AaZsaqIbftLbZi/B1hBN54wskcWCzbupeSYjS2orxe+3EiHFk04t2crt6OYRrKCbjxhZI8sjlX4+GbDXrejRJRvt+7j2637GXdGJ7tfqAc0qqCLyPki8q2IbBGRD0VkYKCCGVMfp+Y0p0lCrI0aracXv9pEWlIc1wy2OxJ5QWOv0J8B7sbfu+TvwGMiMrbRqYypp8S4WIZ3zWDu2t34p/MxJ7Ntfymf5e3k+iEdSEkMyw5vpp4aW9B3qeqXqnpIVecCFwL3NT6WMfU3MjeLbftLWbfrsNtRIsIr8zcDcNPpnVzNYQKnQQVdRN4RkV8DX4vIH0Uk3llVCZQFLJ0x9XB2j0wAa3apg8NlFby1cCsX9GlNdrpNPewVDb1CfxIoAdKBC4ANIjIXyMc/F7kxIdemWTK5rdOYs9oK+sm8u6SAQ2UV3DbCBhJ5SZ0azkTkeeBuVS0BcCbh+qrK+ligJzAA6BeEnMbUyfm9W/PknHXsOVxGRmqi23HCUqVPeenrTQzq2Jz+7dPdjmMCqK5X6OOA1JpWqmqlquap6muqem9gohlTf6N7t0IVZq8ucjtK2Jq5qoiCYhtI5EV1LejWQdVEhF5tmpKdnsz0lVbQa/LSV5to1zyZ83vZQCKvqU8buvUFM2FPRBjduzVfrdvD4TIbNXqi5YX7WbS5mJtP70RcrI0r9Jr6/I8+JyL3ishIEWkatETGNNLo3q04Vumzybqq8eJXm0hNjONHp9pAIi+qT0FvC/wWmA3sE5F8EXlTRH4lImeJSFpwIhpTP4M7taBlSoI1u5xgx4FSPl2+gx+d2p60pPiTv8BEnPoMD7sU2I3/HqODgMHO4xIgBfDV8/2MCYrYGOHcnq34dMUOu9doFa/M34JPlZttIJFn1fUKXQHUb42qvqGqv1TVEUBToC9wS7BCGlNf5/duxeGyCpusy1FyrIK3Fm1lTJ/WtG/RxO04Jkga3cvFKfIrVfW1AGUyptHO6JpBSkKsNbs43l1SyIHScuuq6HF1LeijgP1BzGFMQCXFx3J2jyxmriqi0hfdHbSODyQa0CGdQR1buB3HBFFdC/pdwGUAIpIrItNFJE9E3haRC4MXz5iGu6Bva/YcLmPRpmK3o7hq5qoituwt4Ta7X6jn1bWgnwUsd56/CRQDXwKdgU9E5CURscFHJqyck5tFk4RYPl6+3e0ornrxq420a57M6N42kMjr6lrQU4GDItIPeEFVr1fVO1R1CP5ifwEwPlghjWmIJglxjOrZimkrdlBe6XM7jiuWFexn8eZ93HJGjg0kigJ1/R8uAtoAI4HJVVeo6pfAPcDtAU1mTABc0q8N+0rK+Xr9HrejuOKFLzeSZgOJokZdC/pU4HngXqB3NeuXAl0DFcqYQDmrRyZpSXF8snyH21FCrnBfCdPydnL9aR1ItTsSRYW6FvT7gG+Bz4FcEZkgIlVnX7wBiL6fGBP2EuNiOb9Xa6bn7aSsotLtOCE16evNADaQKIrUqaCr6kFV/amqjlXV54COwB4RWSsihcB/AU8FM6gxDXXJKW04VFbBF2t3ux0lZA4eLeftxQVc1LcNbe2ORFGjQZ+SqOoDQH/gFfxt6teq6hMBzGVMwJzRNYPmTeKjqtnlnUUFHC6r4Kd2R6Ko0uCPvZ0pAP5bVX+lqu8HMpSIXCMiK0XEJyKDqyzvJCKlIrLMeTxbZd0gEVkhIutF5AnrRmmOi4+N4YK+bZi5qigqptStqPQxaf5mTstpQd92zdyOY0IoXPsx5QFXAvOqWbdBVfs7jwlVlj+Dv+tkN+cxJvgxTaS4amA2peWVTF3h/av0T5bvYNv+UrtfaBQKy4KuqqtVdW1dtxeRNkBTVf1GVRV4Fbg8WPlM5BnYoTk5GSm8t7TQ7ShB5fMpz8zdQLesVEblZrkdx4RYWBb0k8gRke9E5AsRGeEsywaq/qQWOst+QETGi8gSEVlSUeH9P7+Nn4hw9aB2LNpUzNa9JW7HCZo5a3axtugQd5zdhZgYa3WMNq4VdBGZ5cwHc+LjslpetgPooKoDgF8Bbzp3T6ruzK12RiZVnaiqg1V1cFyc9c2NJlcMyEYE3v/Wm1fpqsrTc9eTnZ7MJae0dTuOcYFrFU1Vz23Aa8qAMuf5UhHZAHTHf0Xersqm7YDonsDD/EDb9GTO6JLB+98Wcveobp67gl24qZhvt+7nD5f1Jt6G+UeliPpfF5FMEYl1nnfG/+HnRlXdARwSkaFO75YbgQ9djGrC1NWD2lG4r5SFHpyB8em5G8hITeDawTbMP1qFZUEXkSucAUvDgE9FZLqz6kxguYh8D7wHTFDV4z+ZdwAvAOuBDcC0EMc2EWB079akJcbx9uKtbkcJqLxtB5iXv5tbzsghKd5uuRetwrIRWVWnAFOqWf4+UG2fd1VdAvQJcjQT4ZITYrlqUDveXLiV/7q4jJapiW5HCohn5m4gLTGOscM6uh3FuCgsr9CNCaYbTuvAsUofk5d448PR/KJDTM3bwdhhHWmaFO92HOMiK+gm6nRrlcbQzi14c9EWT9ye7vFZ62gSH2vD/I0VdBOdfjK0IwXFpczLj+wJu9bsPMinK3Zwyxk5NE9JcDuOcZkVdBOVzu/VmozURF5fsMXtKI3y+Kx1pCXGcduIHLejmDBgBd1EpYS4GK4f0p45a3exec8Rt+M0yMrtB5iWt5NbhueQ3sSuzo0VdBPFxg7tSHxMDC98tdHtKA3y+Kx1pCXFcetwuzo3flbQTdTKaprElQOzeXdJIXsOl7kdp16+L9jPjFVF3Do8h2bJ1rPF+FlBN1HtthGdKavw8eo3kdOWrqr899TVtExJsKtz82+soJuo1jUrlXN7tuLVbzZTciwyZt+cs2YXCzcVc8+53UizfuemCivoJupNOKsz+0vKeWdxgdtRTqqi0sdfpq0hJyOF64Z0cDuOCTNW0E3UG9ypBUM6teCZuRs4Wl7pdpxavbe0kHW7DnPfmB42o6L5ATsjjAF+dX53dh0qC+t+6YeOlvO3mfkM7JDO6N6t3Y5jwpAVdGOAoZ1bMrxrBs/M3cCRML2R9D9mrWPP4TIeurQ3dg90Ux0r6MY4fnV+d/YeOcak+ZvdjvIDa3ceYtL8zVx3agf6tUt3O44JU1bQjXEM7NCcUblZPDt3A7sPhU+/dFXldx/lkZYUx69H93A7jgljVtCNqeI3F/WktLySv81Y63aU//PBt9tYsLGYe0f3sAm4TK2soBtTRZfMVG4+vRPvLCkgb9sBt+Ow6+BRfv/xSgZ3bM51p1o3RVM7K+jGnODno7rRokkCD320Ep+L86WrKr+ZkkdZhY9Hru5HrMduam0Czwq6MSdolhzPfRfksmTLPl5f6F43xo++386s1UX8x/k96JyZ6loOEzmsoBtTjWsGtWNEtwz+Mm0NBcUlId//1r0l/PZ/8xjQIZ1xNl+LqSMr6MZUQ0T4y1X9EOD+D5aHtOmlrKKSu976FgGeuG6ANbWYOrOCbkwNstOTefCiXny9fi/PztsQsv0+PHUNywsP8Og1p9C+RZOQ7ddEPivoxtTi+iHtubhfG/42I59Fm4qDvr93Fm9l0vzNjDsjx4b3m3qzgm5MLUSEh6/sS4cWTfjZG98GtT19/vo9PDgljxHdMvjNhblB24/xLivoxpxEWlI8E8cO4lhFJbdMWsyBkvKA72NZwX5uf20pnTNTeOqGgcTZTIqmAeysMaYOurVKY+KNg9m6t4SbJy3iQGngivrywv2MfXEhzVMSmHTLEJraTStMA1lBN6aOhnZuyZM/HkDetgPc8MICio8ca/R7zsvfzQ3PL6RZcjxvjR9K2/TkACQ10SosC7qIPCoia0RkuYhMEZH0KuseEJH1IrJWREZXWT5IRFY4654Qm1/UBMHo3q2ZOHYw+UWHufSfXzV4egCfT3nxq03cMmkx2c2TmXz7MLKtmJtGCsuCDswE+qhqPyAfeABARHoB1wG9gTHA0yIS67zmGWA80M15jAl1aBMdRuZmMfn2YVT6lKuemc8zczdQXumr8+u37D3CTS8v4o+frGJkjyzeu+N0uzI3ASGq7s1VURcicgVwtareICIPAKjqw8666cBDwGbgc1XNdZZfD5ytqrfX9t4pKSl65MiRIKY3XrbncBm/+WAFM1YV0SUzhQlndeGSU9qSFB9b7fbrdx3itW+28OaircTFxPDgRT254bQOdrMKUy8iUqKqKdWtiwt1mAYYB7zjPM8GFlRZV+gsK3een7jcmKDJSE1k4o2Dmb26iEc+W8u97y3ndx+t5PQuGeS2TiO9STwVPqVwXwlLNu9jzc5DxMUIVw9qxy/P606rpklufwvGY1wr6CIyC6hu5MSDqvqhs82DQAXwxvGXVbO91rK8uv2Ox980Q0KCzS1tGm9Uz1ack5vF1+v3MjVvBws27mXOmiKOzxaQlhhHv/bNePDCnlwxMJuM1ER3AxvPcq2gq+q5ta0XkZuAi4FR+q92oUKgfZXN2gHbneXtqlle3X4nAhPB3+TSoPDGnEBEGN4tg+HdMgA4VuGjtLyS2BghNTES/hA2XhCWH4qKyBjgPuBSVa06NO8j4DoRSRSRHPwffi5S1R3AIREZ6vRuuRH4MOTBjXEkxMXQLDneirkJqXA92/4JJAIznQ+MFqjqBFVdKSKTgVX4m2LuVNVK5zV3AJOAZGCa8zDGmKgR9r1cgsl6uRhjIk1tvVzCssnFGGNM/VlBN8YYj7CCbowxHmEF3RhjPMIKujHGeIQVdGOM8Qgr6MYY4xFW0I0xxiOsoBtjjEdYQTfGGI+wgm6MMR5hBd0YYzzCCroxxniEFXRjjPEIK+jGGOMRVtCNMcYjrKAbY4xHWEE3xhiPsIJujDEeYQXdGGM8wgq6McZ4hBV0Y4zxCCvoxhjjEVbQjTHGI6ygG2OMR1hBN8YYj7CCbowxHmEF3RhjPCIsC7qIPCoia0RkuYhMEZF0Z3knESkVkWXO49kqrxkkIitEZL2IPCEi4to3YIwxLgjLgg7MBPqoaj8gH3igyroNqtrfeUyosvwZYDzQzXmMCVlaY4wJA2FZ0FV1hqpWOF8uANrVtr2ItAGaquo3qqrAq8DlwU1pjDHhJc7tAHUwDninytc5IvIdcBD4rap+CWQDhVW2KXSW/YCIjMd/JX/865IG5ooFKiPgdY15bRxQcdKtAre/xrzWjk/47dOOT+0aenySa1yjqq48gFlAXjWPy6ps8yAwBRDn60SgpfN8EFAANAVOBWZVed0I4OM6ZFjSiPwTI+F1jdxng46PS1nt+ITfPu34BOH41PZw7QpdVc+tbb2I3ARcDIxS57tX1TKgzHm+VEQ2AN3xX5FXbZZpB2wPRu4qPo6Q1zX2taHenx2f4LzWjk/47TPgjl/5hhURGQM8BpylqrurLM8EilW1UkQ6A18CfVW1WEQWAz8HFgJTgSdVdepJ9rNEVQcH7RuJcHZ8amfHp3Z2fGoXjOMTrm3o/8TfvDLT6X24QP09Ws4E/iAiFfjbrCaoarHzmjuASfjbl6Y5j5OZGODcXmPHp3Z2fGpnx6d2AT8+YXmFbowxpv7CstuiMcaY+rOCbowxHhG1BV1ExojIWmeqgPvdzhMORGSzM33CMhFZ4ixrISIzRWSd829zt3OGioi8JCK7RCSvyrIaj4eIPOCcT2tFZLQ7qUOnhuPzkIhsqzI9x4VV1kXN8RGR9iLyuYisFpGVInK3szy450+g+0FGwgP/QIANQGcgAfge6OV2LrcfwGYg44RljwD3O8/vB/7qds4QHo8zgYFA3smOB9DLOY8SgRzn/Ip1+3tw4fg8BPxHNdtG1fEB2gADnedp+Kcw6RXs8ydar9CHAOtVdaOqHgPeBi5zOVO4ugx4xXn+ClE0pYKqzgOKT1hc0/G4DHhbVctUdROwHv955lk1HJ+aRNXxUdUdqvqt8/wQsBr/6PWgnj/RWtCz8Y8yPa7GqQKijAIzRGSpM0UCQCtV3QH+kxTIci1deKjpeNg59S93OTOlvlSlSSFqj4+IdAIG4B8jE9TzJ1oLenVT61r/TThDVQcCFwB3isiZbgeKIHZO+T0DdAH6AzuAvznLo/L4iEgq8D5wj6oerG3TapbV+/hEa0EvBNpX+ToUUwWEPVXd7vy7C/8cOkOAImc2y+OzWu5yL2FYqOl42DkFqGqRqlaqqg94nn81G0Td8RGRePzF/A1V/cBZHNTzJ1oL+mKgm4jkiEgCcB3wkcuZXCUiKSKSdvw5cD7+ydI+Am5yNrsJ+NCdhGGjpuPxEXCdiCSKSA7+OfkXuZDPVceLleMK/OcQRNnxcW6w8yKwWlUfq7IqqOdPuA79DypVrRCRu4Dp+Hu8vKSqK12O5bZWwBRnqoU44E1V/cyZI2eyiNwKbAWucTFjSInIW8DZQIaIFAK/A/5CNcdDVVeKyGRgFf4pUe9U1YZOxxoRajg+Z4tIf/zNBZuB2yEqj88ZwFhghYgsc5b9hiCfPzb03xhjPCJam1yMMcZzrKAbY4xHWEE3xhiPsIJujDEeYQXdGGM8wgq6McZ4hBV0Y4zxCCvoxtSBiKQ683yf6nYWABF5TkT+x+0cJrxYQTembu4Dlqjq4uMLRGSSiKiIvH/ixiJyubOu4oTtZ1X35s62P6lHnj8Ad4hI53q8xnicFXRjTkJEkoA7gOeqWb0VuEREWp2wfDywJViZVHUbMBv4WbD2YSKPFXQTdUTkXRGZXuXrTBEpEZFBNbxkDJAMzKhm3TpgAXBzlffrAJwHvNzAfGc7V+wnPjafsOkUoD5X9cbjrKCbaPRvU5Wq6m78Nx+4vIbtzwK+U9WKGtZPBG5zZtgDuA3/1XNDr9Dn47+F2fFHb/xTqX5+wnYLgVYi0rOB+zEeYwXdRKMT554GOEzNd2PKAbbV8n7vAS3wzzQYC4zDX+Src7aIHD7xUXUDVT2mqjtVdSewF3gK2AhMqOb7AP+9cY2JzulzTdQrBFJFpJmqHnDmfx8O3FnD9snAgZreTFWPishrwE/x3xA4DvgYuKGazRfyr/mwq1pXw9s/g/+Xz1BVLTth3dEq+Yyxgm6i0vEr2/b4C/XDwE7ggxq2343/Crw2zwHfAR2Al1W1/F8tMP+mVFXXn7iwum1F5NfAlcAwVd1TzXsdz7T7JNlMlLCCbqLR8YLezulXPhYYrqpHa9j+W+Cu2t5QVVc7NwM5g+qvwOtFRC7H3zVxjKqurWGzvkAl/l8kxlhBN1FpG+AD7sZ/z8vzT3LHqmnA30SkvaoW1LLdaCBJVYsbE05EegOvAw8Ba0SktbOq0vkA97izga9OcvNhE0XsQ1ETdZzeKkXAKcCoqoOFath+NTAX/5V8bduVNLaYO04FUvA3Be2o8qg6qEmAH1N933gTpewWdMbUgYiMAN4GuqlqSRjkuRb4T6C/x+/NaerBrtCNqQNV/RL4Pf4ujOEgEbjFirmpyq7QjTHGI+wK3RhjPMIKujHGeIQVdGOM8Qgr6MYY4xFW0I0xxiOsoBtjjEf8fz2+Ouoc9G+mAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sim = ares.simulations.Global21cm(**m17.dpl)\n", "sim.run()\n", "sim.GlobalSignature() # voila!" ] }, { "cell_type": "markdown", "id": "837bd2db", "metadata": {}, "source": [ "Hopefully this results *exactly* in the solid black curve from Figure 2 of `Mirocha, Furlanetto, & Sun (2017) `_, provided you're using *ARES* version 0.2. If it doesn't, please contact me.\n", "\n", "**NOTE:** If using a newer version of the code, this may shift slightly due to changes in the default cosmological parameters and halo mass function. If you pull a recent version and the result is significantly different, please get in touch with me.\n", "\n", "Alternatively, you can use the ``ParameterBundle`` framework, which also taps into our collection of data from the literature. To access the set of parameters for the \"dpl\" model, you simply do:" ] }, { "cell_type": "code", "execution_count": 9, "id": "16d2b83c", "metadata": {}, "outputs": [], "source": [ "pars = ares.util.ParameterBundle('mirocha2017:dpl')" ] }, { "cell_type": "markdown", "id": "50b64998", "metadata": {}, "source": [ "This tells *ARES* to retrieve the ``dpl`` variable within the ``mirocha2017`` module. See :doc:`param_bundles` for more on these objects." ] }, { "cell_type": "markdown", "id": "4bdba976", "metadata": {}, "source": [ "### [Mirocha, Furlanetto, & Sun (2017)](http://adsabs.harvard.edu/abs/2017MNRAS.464.1365M) (``mirocha2017``)" ] }, { "cell_type": "markdown", "id": "877ad4fc", "metadata": {}, "source": [ "This model has a few options: ``dpl``, and the extensions ``floor`` and ``steep``, as explored in the paper. \n", "\n", "Non-standard pre-requisites:\n", " * High resolution optical depth table for X-ray background. To generate one for yourself, navigate to ``$ARES/input/optical_depth`` and open the ``generate_optical_depth_tables.py`` file. Between lines 35 and 45 there are a block of parameters that set the resolution of the table. Make sure that ``helium=1``, ``zi=50``, ``zf=5``, and ``Nz=[1e3]``. It should only take a few minutes to generate this table.\n", " \n", "The following parameters are uncertain and typically treated as free parameters within the ranges denoted by brackets (not all of which are hard limits):\n", "\n", " * ``pop_Z{0}``, :math:`[0.001, 0.04]`\n", " * ``pop_Tmin{0}``, :math:`[300, \\sim \\mathrm{few} \\times 10^5]` (``pop_Tmin{1}`` is tied to this value by default).\n", " * ``pop_fesc{0}``, in general can lie in range :math:`[0, 1]`, but for consistency with observational constraints on :math:`\\tau_e` (from, e.g., *Planck*), it's probably best to limit values to :math:`\\lesssim 0.2`.\n", " * ``pop_fesc_LW{0}``, :math:`[0, 1]`\n", " * ``pop_rad_yield{1}``, :math:`[10^{38}, 10^{42}]` :math:`2.6 \\times 10^{39}` by default\n", " * ``pop_logN{1}``, :math:`[18, 23]`, :math:`-\\infty` by default.\n", "\n", "**NOTE:** Changes in the metallicity (``pop_Z{0}``) in general affect the luminosity function (LF). However, by default, the normalization of the star formation efficiency will automatically be adjusted to guarantee that the LF does *not* change upon changes to ``pop_Z{0}``. Set the ``pop_calib_lum{0}`` parameter to ``None`` to remove this behavior.\n", "\n", "To re-make the right-hand panel of Figure 1 from the paper, you could do something like:" ] }, { "cell_type": "code", "execution_count": 10, "id": "adf032c8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "# WARNING: revisit scalability wrt fesc.\n", "# WARNING: revisit scalability wrt fesc.\n", "# WARNING: revisit scalability wrt fesc.\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEsCAYAAADTvkjJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABSL0lEQVR4nO3dd3zURfrA8c+TSgoEQggttJDQQZrSu1RBQcWzAfaC+rPceerZ7/Qsdydnb9gFAVEEFUQp0nvvIYEAodcE0rM7vz++yxkg2bRv2u7zfr32JPudnZl8b/NkMjvzjBhjUEopVfn5lHcHlFJK2UMDulJKeQgN6Eop5SE0oCullIfQgK6UUh5CA7pSSnkIDehKKeUhNKArpZSH0ICulFIeQgO6UjYRkfUiYvJ4ZJZ335R38CvvDijlQa4DGuX6+g2gA/Bl+XRHeRvRXC5K2U9E5gEDgA+NMfeVd3+Ud9ApF6VsJiLLsIL5fzSYq7KkAV0pG4nIeqA78HdjzF/Kuz/Ku+gculI2EZFtQCvgKWPMq+XdH+V9dA5dKRuISBwQC3wITMl1aZ8xZm/59Ep5Gw3oSpWQiAjgzOfyJmNM+zLsjvJiGtCVUspD6IeiSinlITSgK6WUh9CArpRSHkIDus1EZIuIOF2PVBGJFJEOIpLlyuuRJSLtyqlvkSLiyNW/RNfzZdI/ERkuItm52t9YVu2LSLKrfmeu5xLPP+d6fFUGbR7O1V66iLQtgzYrxPsvP3n9zJR3n8qaiLQVkbRc9+GV4tSjAd1GIjIEaAPUN8b4AAJ8C/wM7DXGCLAXmFNOXTwOxLj6FgY0EJGXy7B/6cDzrvabAZeJyPgyav9T4Ok8nt9ujPFxPcaUQZuzgVDXPTgF/FIGbVaU998l3PzMeJtlwAbXPQgHJhenEg3opaOOiAQDvkA8UBcY77o2HqhXHp0ylvNrosOwfniclFH/jDHzjTH/dP07HsjESl5V6u0bYx4F4uyut6htGmPuNMakub5cifX/Q6m2SQV5/xXg4p8ZryEizYGqQE8AY8wZY8yWYtWlyxbt5ZpGuMz1ZaoxJlREjGt0dL7MBV+Xcf8CsUbKAhwyxtQvj/6JyPVYI7EWwM6yaN/V5jTXKAjXlFMjwABngV7GmM2l2eZF1zKApcaYK0uzzYr0/stLXj8z5didMicifwVeAVKBUCAFuMwYs6+odekI3UauucnWWLk8qgF+rqx7FYYxJtP1g94OiBSR/yvrPohIE2Aq8IMxZldZt5/LWKAKEID1Q7SkrBoWkd1Yv0gGllWbFVFl+JkpA1WwYvGbrp/NLGBxcSrSgG6vl7FGGCuMMWeB5VhTCojIgNz/LW+uP+kOAPdC2fVPRKoCu4AdxphRuZ4v8/tjjFns+gWXAzyK9WdvqRORhUA0EGvK6E/kivb+yyXfnxkvshjAGPOs6+uvgTrFqUgDur02AtVEJMq1HbwzsB84ArznKvMecLg8Oici3c6vqhCRukBDYFNZ9c91Tw4BZ4wxbXJdKpf7c1Fwex5rKqq02/wc6AP0NMYklXZ7LhXi/ZePjeT9M+M1jDG/Aw4Rucv11LVYH5gXqzJ92PgAErE+aHRizcvWxHqTZmP9iZ0NdCinvj2Wq29OIMH1fJn0D/i3q43cffiqLNrHmp80uR5LXf//nO9HBjCgDNo0F92DU2XQZoV4/7np8yU/M+Xdp3K4B38FHLnei+2KU49+KKqUUh5Cp1yUUspDaEBXSikPoQFdKaU8hAZ0pZTyEJUyoLuSKm0RkY0istb1XLiI/CYiu13/rVHe/QQQkZXl3YeClGcfva1tb2mzqCpDH0ubHfegUq5ycW3Z7myMOZHrudexloC9KiJPAjWMMU+UVx9z9cthjPEt5mu/McbcZFc5N68vdh9Lytva9pY2i6oy9LG02XEPKuUIPR/XAF+4/v0FMLL8umKbITaXU0p5sMo6Qt8LnMbaKPGhMeYjETljjKmeq8xpY8wl0y4icg9wj+vLTmXRX6WUspPJJ7laZQ3o9Ywxh1yJ8H8DHgJmFSagX1SPqYzfv1LKe4lIvgG9Uk65GGMOuf57DJgBXAEcdeUnOZ+n5Fj59VAppcpepQvoIhLiytiHiIQAg4CtwCxgnKvYOGBm+fRQKaXKR6WbchGRaKxROYAfMNkY87KI1ASmYWUQ3A+MNsa4zVimUy5KqcrG3ZRLpQvodtKArpSqbDxuDl0ppdSlNKArpZSH8CvvDpQHERkBjCjvfiillJ10Dt2Lv3+lVOWjc+hKKeUFNKArpZSH0ICulFIeQgO6Ukp5CA3oSinlIXTZolJKeQhdtujF379SqvLRZYtKKeUFNKArpZSH0ICulFIeQgO6Ukp5CA3oSinlIXTZolJKeQhdtujF379SqvLRZYtKKeUFNKArpZSH0ICulFIeQgO6Ukp5CA3oSinlIXTZolJKeQhdtujF379SqvLRZYtKKeUFNKArpZSH8KiALiJDRGSXiMSLyJPl3R+llCpLHjOHLiK+QBwwEEgC1gA3GWO2u3mNzqErpSoVd3PonrTK5Qog3hizB0BEpgDXAPkGdICt8auoG9GImtXrlEEXlVJlxRjD0ZRMsnKc/3vOJ/0EkpWaV+l8v8525nAi8+SlV3INBt2PC43rf/MvZIxB3JYoHE8K6PWBA7m+TgK6FPSim5bdhY8xDHY25tVxM/Hx9S21DnoKYwwZjgzOZp393yMlK4XU7FSyndlkO7Kt/zqzcTgd+Pv6E+AbQKBvIMF+wUQGR1I7uDYRQRH4+uj9VvY6k5bFN6sPMH3dARKOpwKGvj4buc33V/r6bip0PWdFmFqtKpOrhXLcr3KEysrRy8LJ60+QS37hicg9wD3nvx4b0J2dZzczJ3AfoZPG8NzYyaXZx0rBGMPR1CPsObaZvad2si8lkaPpJziWlcyxzGROZaeQ48wpcTv+Pv7EVI+hZc2WtK/Vnh71exAZHGnDd6C8UcLxc3y6dC/frU8iI9vJFY1q8GSXvXTe+yE1UnaQHliLHQ3v41xIw3zrEIEc42BB1jZ+zFjPOZNBa78orvJviq/8MfjIK9hIns+6e8WlbZsCagH4P7bmX4enzCGLSDfgBWPMYNfXTwEYY15x8xpjjCEnJ5u7PunBpsA03mj9Iv0uv66Mel3+HI4cdiYuYuWuX9l1cgsJ2Uc5IFmk5/q4PNTppE5ODpE5DiIdDmo6nFTzD6Vq9UZUrdWaao17U7V6Y0ICQgjwCcDPxw9/H3/8fPzw8/Ejx5lDpiOTTEcmqdmpHEs7xpHUIySdTWLnqZ3sPLWT05mnAWhTsw2jYkdxVfRVhPiHlNNdUZXJ/pNp/HdeHDM2HsTf14dR7eszPvoIjda8BIc3QY0m0PtxaHcD+Pq7rWvh/oW8vuZ1ks4l0aVuFx7t9Cita7Yuo++kcNzNoXtSQPfD+lB0AHAQ60PRm40x29y85n8fiu45sI2xv91AvRx/pty5zqOnXnbsWcWctZ+z9cwGtvufI9XHem+EOJ00zhTCM4Pxy6qBIzMSZ2Ytshxh1KgWRu/oavRr5E/V9ENwcjfsXwnJB0B8oOkA6PVnaNStyP0xxrD7zG4WJy1mzt45xJ2OI9Q/lHGtxzG21ViC/YPtvgXKA5w8l8mEeXFMWX0AXx9hXPfG3Nc+kPDlL8HW7yCsAfR9Ctr9CXzdT0YcTzvOK6tf4bd9vxFTPYa/dP4L3et1R6TgkXVZ84qADiAiw4D/Ar7Ap8aYlwsof8Eql1cm3c7knLX8JeJGxl31dKn2tazt2LOO6csmsC5jMwkB1vccmeMkNjucRsGtiKnbm5iY/tSuUZ2QQOvNfy4zh4On09mcdIaFu46xLP4koYF+PNg/hrt6NsHP1wdO7IbN02Dd55B6DJoPg6v+A9XqFaufxhg2n9jMp1s+ZcGBBdQKqsWzXZ+lX8N+dt0KVck5nYYpaw7w2i87Sc3M4aYrGvJgv2hq75oEvz0Hxgk9HoEeD0NAwYOBn/f8zMsrXybTkcn97e9nXOtx+Pu4H8mXJ68J6EV1cUBPy0hl1NddCDDCzDs2VvpRek5ONl/PfY3fkmawJTATI0LTTKFdQEt6t7yRvp2uxs+v8N9j/LFzvDpnB/N2HKNrdDhv39SRWlUDrYtZabD6Q/j9NfANgJHvQsuSpcvZdHwTf1/xd+JOx3Fd7HX8rcvfCPANKFGdqnLbeSSFJ77bwqYDZ+jSJJyXRrYhNuAUzHoQ9i62/lIcPgFqNCqwrrTsNF5Z/Qo/xP9Ah8gO/KPHP2hUreDXlTcN6PnIax36G1Mf5LOMRfy11i2MGVY59yalpp3l7R8eZUHaSg77CzVznHSVGEZ0vJ8e7YeUuP7p65J49oet1AmrwqS7ulCvetAfF08mwPd3w8F1cOUL0PPRErWV7cjmnY3v8OnWT+kQ2YEJfSdQM6hmyb4BVek4nIaJS/bwn1/jqBbkxzNXteKa9vWQ7TNh5oOAgcEvQ8dx1qeLBdhzZg+P/P4IicmJ3N3ubu6/7H78fCrHGhEN6PnIK6Cnpp1lxDddqeUIYOo9G8qpZ8WTlZXJez88zs8pCzjiL8Rm+jAgYjC3DX2ekCB7P2Bct+8Ut326hhohAcwY352aoYF/XMzOgJnjrXnMfs9An8dL3N4ve3/hmWXP0KBqAyYOmqhB3YsknU7jsambWJ14iiGt6/DyqDbUDPKB356Hle9C/c5w/aeFGpUDLDu4jL8s+gsBvgG83vt1utQtcHVzhaIB/SK50ufendf3/4+vbmWacxP/bPIoI3rfUeb9K44v57zK1KSv2R8gRGfBdfVv4tbBT5TqtNH6/ae56aOVtIsK4+u7uhCYe/rG6YAfxsPmKTD0dehyb4nbW314NQ/Mf4CoqlF8OvhTalSpUeI6VcW2KO44D0/ZQI7D8OLVrbm2Y30k9QRMvRUOrIQr7oFBL4NfwVNxxhgm75zM62teJ7Z6LG/3f5u6oXXL4Luwlwb0fOS39f/46UNcM2MgzbND+ezeVeXQs8LbHLecNxY+wroq6dTPNoyKuJo7h/8dvzLaCPHjpkM89M0G7u0TzVNDW1540emAaWNh1xwYMwOi+5S4vVWHVzF+3njaRLRh4qCJ+BewDE1VTk6n4Z2F8UyYF0fz2lV5/9ZONIkIgeNxMOl6OHcMrnkH2l5fqPqMMUxYP4HPtn5G3wZ9ea3Xa5V29ZSmzy2iWjXq0Yto1gWmsnLLr+XdnXy9891fuHvZ3WwNSGMULZh+0zLuHfnPMgvmACMuq8dNVzTgo8V7WL331IUXfXxh1AcQEQvf3gbJB0vcXpe6XfhHj3+w/th6/rHyHyWuT1U86VkOxk9azxu/xXHNZfX4fnx3K5gnLoNPBkJWKtz2U6GDudM4eWnlS3y29TNuaHYD/+3730obzAuiAT0fdw54GX8DX6/Md19SuTmdfILxH/fhw3NzaZDtx0fd3ufv474lNCSsXPrzzFWtaFAjmCe/23xB3gwAAqvCnyZBTgbMeqigpBeFMix6GHe3vZsZ8TP4MeHHEtenKo6T5zK56eOVzN1+hGeuasmEP7UnOMAPds6Gr0ZCaCTcNQ+iOheqvhxnDn9b+jemxU3jjjZ38EzXZzw63YQG9Hw0a9SOrjk1WeF3nPj9+W+1LWtLN/zEmGn9WBJwiqE5Dflq7HI6tuhVrn0KCfTjhatbsedEKl+uSLy0QEQMXPkiJMyHDV/Z0ub49uPpGNmRl1e9zIGzBwp+garw9p5I5dr3l7PjcArv39KRu3pFWxt7tv0A08ZA7TZwx1wIb1Ko+pzGyfPLn+fnPT/zfx3+j0c7PVohNwrZSQO6G7dc8STZAp/M+1t5dwWAt6c/xqMbnyTZ18kTkbfy+p0/E1SlYmyP79c8kt7NavHm/N2cPJd5aYHL74LGvWDu09b8Zwn5+fjxSq9X8MGH55Y9hzd/FuQJdh05y+gPVnA2I4dv7unKkDauDys3fwvT77BWsoydCcHhharPGMMrq15hVsIsHmj/AHe3u7sUe19xaEB3o/tlQ+mQGcRiEjh55ki59SMnJ5u/fnIVH6X+RpNsPz7q+yW3Dn2i3PqTFxHh2ataci4zh4lL915awMcHrnoDstNgodsNvIVWL7Qej3Z+lLVH1/LjHp16qay2Hkzmxo9W4CMw7d6udGzoWr209XtrT0Oj7nDrd1ClWqHrfHP9m0zZNYXbWt/Gve1KvsKqsvDKgC4iI0Tko8KUHdXsLlJ8ffjo5/LZZHQ6+Th3fdKdOX776ZMZzudjltAyulO59KUgsbWrMqxtXb5cnsiZtKxLC9RqBpffDeu/hCNbbGnzutjraF+rPf9e82/OZJyxpU5VdjYeOMPNH68kOMCPafd2IyayqnUhfh58fw807Ao3T4PA0ELXOXnHZD7Z+gmjm43msU6Pefw0S25eGdCNMT8aY+4puCSM7HcvzTJ9WJC+lozMtNLu2gUOH9/PPVMHsj4wnT/5dOCtOxcQHFS1TPtQVA/1jyE1y8GnyxLzLtD3CQisBgvsGaX7iA/PdnuW5KxkPtpSqN/RqoLYcTiFsZ+sIizYn6n3dqVxhGv68MBqmDoGarWAm6YUKh/LeUuSlvDamtfoG9WXp7s87VXBHLw0oBfV4NpXc8Rf+OznF8uszb0Hd3L/jOHEB+RwX7VhPDPmy0qRW6ZFnWoMalWbL5YnkpHtuLRAUA3o9iDEzbFSm9qgWY1mjIoZxTc7vyHpbJItdarStfdEKmM+WU1IoB+T7+pKVA1X0D4eB5NGQ9U6MOZ7CKpe6DrjTsfx+OLHaVajGa/1fs2jV7PkRwN6Idw29BnqZRvmnpiD05FHkLLZjj3reGD29ST5O3m41p8Yf+3rpd6mnW7v0YTk9GxmbTqUd4Eu90CVMFhk3/c1vv14/MSPt9a/ZVudqnQcOpPOrRNX4TSGr+7sQoNwVzBPOwWTXTnLx8ywligW0qmMUzw4/0GC/YJ5u//bHrvOvCAa0AshICCQwVV7kxBomPjTC6Xa1vrti3h4wThO+hn+GnUPt131bKm2Vxq6RofTrHYoX65IzHv1SZUw6Doedv5k21x6ZHAkY1qNYU7iHBLOJNhSp7Jfcno2Yz9dTUp6Nl/ecQUxka658Zwsa5ol5RDcOBlqNC50nQ6ng6eWPMXJ9JO83f9t6oR47/nAGtALafzIfxOVbZhxfEapzaUv3fATf14xnnM+hudiH+eGgf9XKu2UNhFhTNdGbD2YwoYDZ/Iu1OU+CAiFFe/a1u6YVmMI8gviky2f2Fansk+2w8n4SevYdzKVj8Z2pk1910Y4Y2D2X2DfUms7f4MrilTvh5s/ZPmh5TzV5SlaR1Ss04XKmgb0QqoSGMyoWqNI8hfemfGY7fXPXTGZJzc8gQN4ue3fuarnbba3UZZGdYwiJMCXb1btz7tAUHVof7OVkdGGdekANarUYHSz0czeO1s3G1Uwxhie/WEry+JP8sq17ejWNFe2zDUTYf0X0Osv1jFxRbDs4DI+2PQBVze9mutivefoyPx4ZUAvyrLF3O4a/gLNMn34KXUJh47vs60/Pyz8kOd3vkyggdevmOARZ5qGBvoxrG1d5mw9QnpWPp87XHEPOLJg7We2tTuu9Th8xIfPttpXpyq5DxfvYcqaAzzYL4brO0X9ceHgepj7N4gdBP2KdkrYsbRjPLnkSWJqxPBM12e8bkVLXrwyoBdl2WJuPr6+3NX6L5z2FV6eMc6Wvnz20z94OfFtqjuECb0/oWvbQbbUWxFc2zGKc5k5/Lo9n01ZEbEQMxDWfmLNodogMjiSkTEjmRk/k1MZpwp+gSp1v+86xmu/7GR4u7o8NrDZHxfST8O34yC0Noz60Np8VkhO4+TZZc+S6cjkP33+Q5BfUMEv8gJeGdBLYmiPMVyZE8XiwJNM+XVCiep6ffLd/PfEVOpn+/DWwMm0i+1qUy8rhi5NwqlfPYjv17vJstjlPjh3FHbMsq3dW1reQpYzi+93f29bnap4DpxK45GpG2leuyr/uv4yfFwHkuN0woz7IeUwjP6i0Fv6z5uycwrLDy3nL53/QpOwwuV28QYa0IvhmdFfEZVteO/ARHbsWVfk1zsdDp74dARfZa+kbWYVPh79C80atSuFnpYvHx9hZId6LNl9nGMpGXkXatofqjeydo/apGn1pnSp24Wpu6aS48yxrV5VNBnZDu6ftA6H0/DhmE4EBeRaF77yPWsvwuCXIapoO5/3nNnDG+veoFf9XoxuNtrmXlduGtCLoUZYLZ7s8DIZPvDMvDs4nXy80K89dHwfd03swWzfRHplhjPxtsXUqlGvFHtbvkZ1iMJp4Octh/Mu4OMD7W+BvYvgtH2fS9zU4iaOpB7h9wO/21anKprnZ25j68EU/vun9jSqmSuJ3NFtMP9FaD7M+hylCLKd2Ty55EmC/YL5e4+/67z5RTSgF1OfTtdwZ/Wr2R3g4P6pgzh+Op9NNLn8vPRzxs28ivWB57iO1rxz5wKqBHr2BoiYyFCa1Q5lzlY3yc3a3wQIbJxsW7t9o/pSL6Qek3faV6cqvOnrkpi69gAP9Y9hQMvaf1zIybRytFQJgxFvFepA59y+2PYFO07t4LluzxERFGFzrys/DeglcO+oV7g7ZCA7ArK55btBTJz1HDk52ReUcToczFv1LQ9/PICn4/+NwfBik0d4YdyUSrGV3w5D29RlTeIpjp3NZ9qlekOI7msFdKcz7zJF5Ovjy+jmo1lzZA37U/JZOqlKReKJVJ6buZWu0eE8cmWzCy8ufBmOboWr34HQWkWqd1/KPt7f+D4DGw3kykZX2thjz6Fnitrw/U+f/y4f732fQ/5CDYeTxtnBBEkg50w6h/0yOO7ng48x9M6uxV9HfEKDOtE29L7y2HkkhSH/XcJLI9twa9d8TmbfMh2+u9PKeR3d15Z2j6YeZdB3g7ir7V081OEhW+pU7mU7nFz//nIST6bxyyO9qBuWa/VJ4jL4/CroNA5GvFmkep3GyZ1z72TXqV3MHDmTWsFF+2XgSfRM0YsUdx16fq4f8AAzx6xmfNVhNM8J45RPBvG+pznnk0UTRzVuC+zJt/2+5u27F3pdMAdoXrsq0REhzNmazzw6QIvhEBgGm6bY1m7tkNp0q9eNWQmzcDhLPwePggm/xbEpKZlXr217YTDPToeZD0CNRjCo6Jk2Z+yewdqja/lz5z97dTAvSNmdJlyBGGN+BH4UEduOMakSGMz9175mV3UeRUQY2rYOHyzaw6nULMJDAi4t5F8FWo6wli9mZ1hf22BkzEgeX/Q4q46sonu97rbUqfK2IuEk7y9K4E+dGzC0bd0LLy56DU7vhbGzipTbHOBE+gn+s+4/dK7dmWtjr7Wxx57HK0foquwNbl0Hh9Pw+y432/zbXAuZKRD/m23t9mvQj2oB1fgh/gfb6lSXOpuRzV++3UTjmiE8N6LVhRcPb4Zlb0GHWyG6T5Hrfmv9W6TnpPNct+d0VUsBNKCrMtGmXhgRoYEs2OkmoDfpA8ERVn4XmwT6BjKsyTAW7F9ASlaKbfWqC70yZyeHktP59+jLCAnM9Ye/IwdmPQTBNWHgP4pc75bjW5gRP4MxLcfoBqJCqFQBXUReEJGDIrLR9RiW69pTIhIvIrtEZHB59lNdysdH6N+iFovijpPtyGcli68ftB4Ju36BzHO2tX1106vJdGQyf9982+pUf1gef4LJq/ZzZ48mdGpU48KLq96Hwxth6GtF3g3qNE5eXf0qEUER3NOuyJk6vFKlCuguE4wx7V2P2QAi0gq4EWgNDAHeExHvWBNYifRvUZuzGTms23c6/0JtroOcdIj7xbZ220S0ISo0irmJc22rU1lSM3P463ebaVwzmD8Pan7hxeQkWPhPaDYEWo8qct0/JvzI5hObebTTo4QGFG3e3VtVxoCel2uAKcaYTGPMXiAeKFpSZVXqesZGEODr437apUFXqFrPOvHdJiLCkCZDWHl4pSbsstnrv+zk4Jl0Xr/+sgu39gP8+gwYJwx9vcgbiM5lnWPCugm0i2jH8OjhNvbYs1XGgP6giGwWkU9F5Pzfd/WB3Amwk1zPXUJE7hGRtSKytrQ7qi4UGuhHl+hw5u84mn8hHx9oORwS5kNWqm1tD2k8BIdxMG/fPNvq9HZrE0/xxYp9jOvWmCuaXDSdsncxbJsBPR+1lioW0SdbP+Fkxkme6vIUPlIZw1T5qHB3SkTmicjWPB7XAO8DTYH2wGHgP+dflkdVee4YMsZ8ZIzpbIzpXBr9V+71bxFJwvFU9p10E6xbDIecDIi3b867WY1mRIdFM2fvHNvq9GbZDifP/LCVemFVeHzwRVMtjhyY8wSENYQeDxe57qOpR/l6+9dcFX0VbSLa2NRj71DhArox5kpjTJs8HjONMUeNMQ5jjBP4mD+mVZKABrmqiQIKTq6iylz/FtbBv26nXRr1gCrVrTNHbXJ+2mXd0XUcTXXzF4IqlC+WJ7LzyFmeG9H6wlUtYOW3P7bdyqToX/Q85e9vep8ck8OD7R+0qbfeo8IFdHdEJPduhVHAVte/ZwE3ikigiDQBYoHVZd0/VbBGNUNoEhHCkt0n8i/k6wfNh1ofjDqy8y9XREMaD8Fg+G2ffevcvdHh5HQm/BZH/xaRDG5d+8KLqSesfC3Rfa2NYkW0J3kPM+JncGPzG4mqGlXwC9QFKlVAB14XkS0ishnoBzwKYIzZBkwDtgO/AA8YY3SvdwXVKzaClXtOkpXjJhFXi+GQkQyJS21rt0lYE2Kqx7DgwALb6vRG//hpOzlOwwsjWl+60ef3V60lp8X4IBSsTURBfkHc3c62TdxepVIFdGPMGGNMW2NMO2PM1caYw7muvWyMaWqMaW6M0YnSCqxnTARpWQ7W73ezfLFpf/ALsnXaBaydo+uOruN0hpu2Vb5+33WM2VuO8FD/GBrWvCj184l4WPcZdLoNajXP8/XubDy2kfn753NHmzsIr1K0NevKUuSALiLNRWSoiFwrIr1ERBeIqiLp1rQmvj7Ckt1uDgYJCIaYAbDzZ9tS6gIMaDQAp3GyKGmRbXV6i8wcB8/P2kZ0rRDu7p1Hkrn5L4JvIPR9slj1v7n+TSKCIri15a0l7Kn3KlRAF5HGIvK6iBzCmtb4GZgOLAJOicgCEblBKkmiBbuzLaqiqVrFn44Nq7PU3Tw6WHOwZw/DoQ22td0qvBV1Quowf7/uGi2qz5clsu9kGi9e3ZpAv4vWnB9YbSVW6/EwhEYWue41R9aw9uha7m57N8H+nn3oS2kqMKCLyL+wPnxsDvwNaAOEAYFAXWAYsBx4FdgoIh1Lrbc2Mcb8aIzRvcTlqGdMLTYfTOZ0alb+hWIHgfjavtqlf4P+rDi0grTsNNvq9XTHz2by9oJ4BrSIpFfsRelrjYFfn4XQ2tDtgWLV/97G94gMiuS6ZtfZ0FvvVZgRelWgmTHmGmPM58aYHcaYs8aYbNcywnnGmGeMMdHAS0DL0u2y8gS9mkVgDCxPOJl/oeBwaNgVdv9qa9sDGg4g05HJ8kPLba3Xk73xWxwZ2Q7+dlUeP947f4YDK6HvU0VOjQt/jM7vaHsHgb6BNvTWexUY0I0x9xljCrWm2xjzrTFmUsm7pTxdu/phVK3i534eHaDZYOvIsuQk29ruWLsjYYFhOu1SSDsOpzB1zX7GdmtM01oXBWxHDsx7ASKaQYcxxar/vY3vUSuoFtc3u77knfVyhZ1DL/DXpoi0L3FvlNfw8/WhR9MIluw+gdtjAGNdiTNtHKX7+fjRJ6oPi5IWke20b527JzLG8NLP26kW5M/DA2IvLbB5KpzcDQOes/YPFNH50fmdbe/U0bkNCrvKxe3R6a5sh5rKThVJr2YRHDyTzt4TbtIA1GpuHSIdZ/+0y9mss6w9oil93Jm34xjL4k/y6JXNCAv2v/CiI9s6iajuZda+gWLQ0bm9ChvQu4vI23ldEJEYYB66M1MVUa8Y68M1t7tGRaxR+t5F1tF0NulatysBPgEsTlpsW52eJsfh5JXZO2haK4SbuzS8tMDGyXBmH/R7ulibiHR0br/CBvSrgHEi8lTuJ0WkITAf2AFUmo+nddlixdCwZjCNagYXbh49O83WXaPB/sFcXvdylhxcYludnubbdUnsOZHKE0Na4O97UajIyYLF/4L6nazVSMXwydZPCK8SznWxlSZ0VHiFCujGmPXA9cDzIjIG/pdXZT5wEBhhjHGz/qxi0WWLFUfPmAhW7jmV/ylGAI17WrtGd9s7q9e7fm/2pexjX8o+W+v1BOlZDv47L45OjWowsFXtSwts+AqSD0C/vxVrdL7r1C6WHVzGLS1voYqfPQeCqyLsFDXG/ArcDXwsIrdgTbOkAEONMbqgVxVLz5gIzmXmsDnpTP6F/IOsZE9xc601zzbpHdUbQKdd8vD58kSOpmTyxJAWl+Zryc6Axf+2DiNpOqBY9X+69VOC/YL5U/M/2dBbdV6Rtv4bY74CngW+xMo3PtAYk1waHVPeoVvTmojA0t1u1qMDNBtkzdce32Vb21FVo4gOi9Y0ABdJTsvm/d/j6de81qUHVwCs/wLOHir26PzguYPMTZzL9c2uJywwzIYeq/MKu2zx1/MPYCCQ7XpMueiaUkVSPTiAdvXDWBpfwDz6+Xlau6ddonqz7ug6UrPtOx2psnt/UQJnM3P465AWl17MzoAlb0CjntCkd7Hq/2LbF4gIY1oVb926yl9hR+gHL3p8A2zM43mliqxHTAQb9p/hXGZO/oXCoqB2G9uXL/aO6k2OM4cVh1bYWm9ldSQ5g8+W7WVk+/q0rFvt0gIbv4ZzR6DvE8UanZ/KOMWM3TMYHj2cOiF1bOixyq1QOwGMMbeXdkeU9+oZE8F7vyeweu9J+rfI4wO482IHwbI3If0MBFW3pe32ke2p6l+VxUmLubLRlbbUWZm9OT8OpzE8NrDZpRcd2db9j7oCGvcqVv3f7PyGDEcGt7fWkFIaKlU+dOWZOjaqQRV/H/fr0cFavmgckGDfARX+Pv50q9eNJQeX4DT2pemtjBJPpDJtbRK3dGlEg/A8Mh5u/Q7O7Idefy7W6DwtO41vdn5Dvwb9iK6eR/pdVWJF3qsrIoOBAUAkF/1CMMaMtalfpUpERgBFPx9LlYoq/r5c3jicZfEFBPSoyyGoBuz+Ddpca1v7vaN68+u+X9lxageta7a2rd7K5q0Fu/H3Fcb3a3rpRafTmjuv3cb6xVoMPyb8SHJmMre30dF5aSnSCF1EXgLmAIOAOkCtix6Vgq5Dr3h6xkQQd/Qcx1Lc7Ab18bWWycX/ZuuhFz3r90QQr16+uOf4OX7YcJBbuzQismoe68J3/QwndkHPR4s1OncaJ1/v+Jo2NdvQvlb7kndY5amoUy73ALcZY9obY4YYY4bmfpRGB5V36BETAcDSgkbpzQZD6nE4vNG2tmsG1aRtRFuWJHnvrtG3F8QT4OfDvX3yGJ0bA0v+A+HR0HpUsepfdnAZiSmJ3NLqlkvXtSvbFDWgO7EOs1DKVq3qViM8JKDggN50ACC250jvGdWTrSe2cirjlK31VgYJx88xc+NBxnZrTK2qeeRU2bPQOjWqxyPWX0nFMGnHJGoF1WJwo+JN16jCKWpAfw+4qzQ6orybj4/QvWlNlsUXkE43pCZEdbY9oPeu3xuDYdnBZbbWWxm8PX83gX6+3JPXOaFgzZ1XrQeX3Vis+hPOJLDs0DJubHEj/r7+Bb9AFVtRA/o/gLYisklEvhKRT3M/SqODynv0jIngaEom8cfOuS8YOwgOrodzBWxGKoKWNVsSXiXc65J1xR87x6xNhxjbrRERoXmMzvevgsQl0P0h8CteRsRJOyYR4BOgKXLLQFED+t+BoYAv1nmiDS56KFVsPWMLOY8eOwgwkGDfiUM+4kPP+j1Zfmg5DqfDtnorurfm76aKv5vR+bI3ISgcOo0rVv3Jmcn8mPAjw5sOJ7xKHmkElK2KGtAfBO4wxrQxxlxpjBmY+1EaHSwNmj63YoqqEUzjmsEFL1+s0846kDjO3jQAver3IjkzmS0ntthab0W1++hZftx8iLHdGlMzr9H5yQTYNRsuvxMCQorVxvS46WQ4Mril5S0l7K0qjKIG9CzAvqTU5USXLVZcPQqTTtfHB2IGWiN0h5t0AUXUrV43fMTHa6Zd3loQT7C70fmqD8HHDy4v3sdm2c5svtn5DV3qdKFZjTx2nirbFTWgfwTcWRodUQqgV6yVTnfTgTPuC8YOhIxkSFpjW9thgWFcVusyr1i+mHD8HD9tPsTY7o0JDwm4tED6GdjwNbS9HqoWL+fK/P3zOZp2lFtb3VqyzqpCK2pArwvcJyLrReQzEfko96M0Oqi8S7foCCudboHLF/tZo0ebsy/2qt+LHad2cCK9gPYruQ9+TyDQz4c7ezbJu8D6LyA7FbqOL3YbU3dOpX5ofXrVL17eF1V0RQ3oTbGyLCYDjYHYXI8YuzolIqNFZJuIOEWk80XXnhKReBHZ5UpDcP75TiKyxXXtLdHdC5VSWLC/lU63oLwuVcKgYTcrDYCNekVZwWfpwUo/s5ivpNNpzNhwkBsvb5j3yhZHtjXd0rgX1G1XrDYSziSw9uhabmh+A77FXLuuiq6oB1z0c/Pob2O/tgLXAhfsxRaRVsCNQGtgCPCeiJx/t7yPtZP1/C+YITb2R5WhHjERbDhwhrMZ2e4Lxg6Eo1sh2b7Mzc1rNKdWUC2Pnnb5ePEegPznzrfPhJSD0O2BYrfxbdy3+Pn4MTJmZLHrUEVXYEAXkU6FrUxEqohIy5J1CYwxO4wxeR1Ncw0wxRiTaYzZC8QDV7jON61mjFlhrF0pXwIjS9oPVT56xkbgcBpW7y1g1+b5Qy/i7Ruliwi9onqx4tAKcpz2feBaURw/m8mUNQe4tmN96lUPurSAMbDiXQhvCrHF29WZnpPOrPhZDGw0UJcqlrHCjNBnisgMERksInmWF5H6IvIUsBvoYWsPL1QfOJDr6yTXc/Vd/774+UuIyD0islZE1pZaL1WJdGxYyHS6tVpAWEPbD73oWb8nZ7PPsun4JlvrrQg+XbaXLIeT+/LK2QJwYBUcWg9d77dWExXDL3t/4Wz2WW5odkMJeqqKozDpc5sDTwJfA1VEZAPW6UQZQDjW9EcT4HfgJmNMoSYfRWQeVsbGiz1tjJmZ38vyeM64ef7SJ435CGu1DiJi34nDyjaFTqcrYk27bJoCOZnF3sl4sa51u+InfixJWkKn2oX+A7XCS07P5qsV+xjWti7RtULzLrTiXahSHdrfXOx2pu2aRtOwph517yqLAn8FG2NSjTHPAlHAGGAtUAVrxUsK8C7Q2hgzoLDB3FXvla4NShc/8gvmYI28c+9IjQIOuZ6PyuN5VUn1io1g97FzHEl2k04XrGmX7FTYZ1/OuKoBVelQu4PHrUf/cnki5zJzeKBvPusXTifCzp+g8+3F3ki07eQ2tp7cyujmozWrYjko9N9UrnnrH4wxjxljRrnS544xxkwwxuwszU7mMgu4UUQCRaQJ1oefq40xh4GzItLVtbplLODuF4Oq4M6n0y1wlN6kF/gG2p6sq1f9XsSdjuNI6hFb6y0vaVk5fLpsL/1bRNKqXh5nhYK1skV84Iri77n7dte3VPGtwoimen5MeaiQR9CJyCgRSQK6AT+LyFwAY8w2YBqwHfgFeMAYcz7xxv3ARKwPShOwDuJQlVTLOtWoGRJQcEAPCLGCut3pdOv3BPCY7IvfrD7A6bRsHsjrNCKwNmmt/wpaXwvV6hWrjbNZZ5m9dzZDmwylWkA+vzRUqaqQAd0YM8MYE2WMCTTG1DbGDM517WVjTFNjTHNjzJxcz691Tdk0NcY8aNzmYFUVnY+P0D0mgqUFpdMFa9rlZLyVe8QmMdVjqBNSxyOmXTJzHHy8eA9dmoTTqVE+q07WfwVZZ6Fb8TcS/bTnJ9Jz0vlT8z8Vuw5VMhUyoCsF0DOmJsfOZrK7wHS6rrxw8fNsa1tE6FXfWr6Y7ShgPXwF9/36gxxJyeCBfvnMnTtyrOmWRj2gXoditWGMYdquabSq2YrWEd57Lmt588qArtkWK4f/HUtX0PLF8GioGVsq2RfTctJYf2y9rfWWpRyHkw8WJdAuKoxervTEl9j5EyTvL9E2/43HNxJ/Jl6XKpYzrwzomm2xcoiqEUyTiJCC87qANe2SuBSyUm1rv0vdLvj7+FfqNAA/bznMvpNpjO8bk/+qkxXvQo3G0Lz4xwJP3TWVUP9QhjbRo4XLk1cGdFV59Iipyco9J92n0wVr2sWRCXvtm/MO9g+mU+1OlTYNgNNpeG9hAjGRoQxqVTvvQgfWQNJqa3RezJwrpzNO82vir4xoOoJg/+AS9FiVVIkCuogMcmVe3CciM0Wko10dUwqsY+nSshys33fafcFG3cE/pFSyLyYkJ3DwnH35YsrKgp3H2HX0LOP7NsXHJ5/R+cp3ITAM2hf/AIqZ8TPJdmYzutnoYteh7FHSEfr7wMNAG2AC8IaIjClxr5Ry6R4TgZ+PsHBXAeeH+gVaKXV3/2blI7HJ/7IvJlWuaRdjDO8sjCeqRhAjLstnGeKZA7B9FnQaC4H57BwtgNM4+TbuWzpGdiS2RmwJeqzsUNKAfswYs8QYc9YY8zswDHii5N1SylKtij9XNAlnwc6jBReOHQjJB+C4ffvcGldrTP3Q+pVuHn1Fwkk2HjjDfX2a4u+bz4/56g+t/15xb7HbWXl4JfvP7md0cx2dVwTFCugiMlVE/gosE5F/iIi/65IDyLStd0oB/VtEEnf0HAdOpbkvGONavmjjapfzyxdXHVlFpqPyvLXf/T2eWlUDub5TVN4FMs/Cui+g1TVQvfjnu3+761tqBNZgUKNBxa5D2ae4I/S3gTSgOjAUSBCR34E4YLYtPVPKpX+LSAAW7jrmvmBYfajdxv7li1G9SM9JZ92RdbbWW1o27D/NsviT3N2rCVX88/mgc8MkyEyBbg8Wu51jacdYeGAhI2NGEuCbxzF2qswVKqCLyMci8r+Pr40xS40x7xhj7jLGdMbKtvgg8AxW4q4KTdehVy7RtUJpEhHC/B0FBHSAFlfB/hVwroA59yK4vM7lBPoGVppdo+8uTCAsyJ+buzTKu4DTASvfgwZdIKr4GRG/3/09DuPg+mbXF7sOZa/CjtDvAPL91MQY4zDGbDXGfGWMedyerpUeXYde+fRvEcmKPSdJyyrg0IkWwwEDcfal8gnyC6Jznc6VYh5955EU5u04ym3dGxMamE927F2z4cy+Em0kynHmMD1uOt3qdqNhtYbFrkfZq7ABXfNgqnLVv0UkWTlOlsWfdF+wTlvr0IsdP9nafq/6vUhMSWR/yn5b67Xb+78nEBzgy+09GudfaMV7UL2h65df8SxJWsLRtKOat6WCKcocuia7UuXm8sbhhAb6FbzaRQRaDoc9C60P/mxy/uT6ijztsu9kKj9uOsStXRtRPTifOe2D62H/cuhyH/gW5nybvE2Lm0ZkUCS9G/Qudh3KfkUJ6B+KyOMi0k9ENDemKlMBfj70bhbBgp3HCs6+2GI4OLKsNek2aVitIY2rNa7Q0y4fLNqDn48Pd/Vskn+hle9BQFXoUPztIklnk1h2cBnXNrsWfx//gl+gykxRAno9rA895wOnRSRORCaLyGMi0kdEqpZOF5Wy9GseydGUTLYdSnFfsGFXCI6wkk7ZqGf9nqw5sob0nHRb67XDkeQMvluXxOjOUURWy2ddQvJB2DYDOo6FKsUfk02Pm46IcF3sdcWuQ5WOogT0q7GWKbbCOhHoZ6zj4F4EFgIF7M1WqmT6tYhEBH7dXsC0i4+vlWgq7lfrrFGb9Krfi0xHJmuOrLGtTrt8vGQPDmPyP/wZYPVHYJzQpfgbibId2cyIn0GfqD7UCcnrSGBVngob0A2Asew0xkwyxjxqjOkFVAPaAreXViftpssWK6eI0EAubxzOL1sPF1y4xXDrwAYbk3V1qtOJIL+gCpes61RqFpNX7eeay+rRIDyf5FhZqbDuc2g5Amrks5yxEObvn8+pjFPc0FzT5FZEJV7l4gry24wxX9nUp1KnyxYrr2Ft6hB39BwJxws49CK6LwSEws4fbWs70DeQLnW6sPTg0oLn8cvQ58v2kp7t4P6+bkbnGyZBxhno+kCJ2poWN436ofXpXq97iepRpaOwAX0AcKYU+6FUoQxpUxeAX7YWcHizfxWIuRJ2zrY20tikV1Qvks4lkZiSaFudJXE2I5vPlycyuHVtYmvn8zHW+Y1EUZdDwy7FbmvPmT2sObKG65tdj49o5u2KqLD/rzwIXAMgIi1EZK6IbBWRKSIyrPS6p9SF6oRVoUPD6szeUohpl1bXQOox2LfctvbPHx5dUaZdvl65n5SMnPyPlwNrI9HpvdCt5KNzPx8/RsWMKlE9qvQUNqD3ATa7/j0ZOAUsAaKBn0TkU8n3OBSl7DW0TR22HUph/8kCknU1G2zlSN/6nW1t1wutR9Owpiw+uNi2OosrI9vBJ0v30Cs2gnZR1fMvuOJd10aiEcVuKy07jVnxsxjYaCA1g2oWux5Vugob0EOBFBFpB0w0xtxkjLnfGHMFVrAfCuictCoTQ89Pu2wrYJQeEALNh8D2mWDjQc99G/Rl3ZF1JGcm21ZncUxbe4AT57Lcj86T1lm5bbrcX6KNRL8k/sLZ7LO6M7SCK2xAPwrUBfoB03JfMMYsAR4Bir8WSqkiaBAeTJv61ZhT0Dw6QJvrIP0U7FlkW/tXNrqSHJPDoiT76iyqrBwnHy7aQ6dGNejSJDz/givegcBq0LFk585M3TWVmOoxdIzUQ8kqssIG9NnAx8DjQOs8rq8D3AwTKhZdtlj5DW1Tlw37z3DoTAGbfGKutI5Y2/a9bW23rtma2sG1mbdvnm11FtX365M4eCadB/u7Ofz5zH7rr5NO4yCw+Pv+tp7YyvaT27mh+Q35t6UqhMIG9CeA9VgbiFqIyH0ikjv74i1AIT6lqhh02WLlN7ydNe0ya9Mh9wX9Aq2Uujt+tG2TkYhwZaMrWX5oOWnZBczjl4KsHCfvLIznsgbV6dusVv4FV7lOJOpyX4nam7prKkF+QYyILv4cvCobhQroxpgUY8zdxpgxxpgPgUbACRHZJSJJwHPAu6XZUaVya1QzhI4Nq/PDhkIc3tzmOuswh3j7RtQDGg4g05FZLrldvl+fRNLpdB65Mjb/EXNGsnUiUetREJbPqUWFkJyZzC97f+Gq6KsIDSjeuaOq7BRrMakx5imgPfAF1pz6DcaYt2zsl1IFGtmhPjuPnGXH4QJyu0T3gaBw2GrftEvHyI6EVwln3v6ynXYp9Oh8/ZfWTtnuxT+RCGBWwiwyHBn6YWglUezdAa4UAP80xjxmjLFvXRggIqNFZJuIOEWkc67nG4tIuohsdD0+yHWtk4hsEZF4EXlLl1F6vqva1sXPR/hhYwGjdF9/a036rtm2pdT19fGlX4N+LE5aTJYjy5Y6C6NQo3NHjjXd0qgn1OtQ7LaMMUzbNY12tdrRIrxFsetRZaeibvfaClwL5LXYN8EY0971yD05+D7W0slY12NI6XdTlaeaoYH0blaLWRsP4XQWsBX/spsgO836kNAmAxoOIDU7lZWHV9pWpzuFHp1v+x6SD5R4I9HqI6tJTEnU0XklUiEDujFmhzFmV2HLi0hdoJoxZoWxkmx8CYwsrf6pimNkh/ocTs5g5d4CTjJqcAWEN4WNk21ru0vdLoT6h5bZapdCjc6dTlg6AWq1gGYlG9NM3TWVsMAwBjceXKJ6VNmpkAG9AE1EZIOILBKRXq7n6gNJucokuZ67hIjcIyJrRWRtaXdUlb6BLWsTGujH9HVJ7guKQPubYd8yOLXXlrYDfAPo06AP8/fPJ9vGjUt5KfTofPdcOLYdej4KPsX/8T6WdoyF+xcysulIAn0Di12PKlvlFtBFZJ4rH8zFj2vcvOww0NAY0wF4DJjsOj0pr+FKnn+DG2M+MsZ0NsZ0zuu6qlyCAny5un09Zm85THJ6AUH1shsBgU1TbGt/aOOhpGSlsPyQffli8jJ9XSFG58bAkjesM1XblOzwiam7puIwDp1uqWTKLaAbY640xrTJ45HvJKcxJtMYc9L173VAAtAMa0See21WFFDAAmXlKW66vCEZ2U5mFvThaFiUteJl02RrasIG3et1JywwjNl7Z9tSX14ysh28OT+OTo1quB+d71sGSauhx/9ZHwQXU6Yjk+lx0+kT1YcG1RoUux5V9irVlIuI1BIRX9e/o7E+/NxjjDkMnBWRrq7VLWMB+z79UhVa26gwWterxjerDxScp7z9LdYOyn3LbGnb39efgY0GsvDAwlI7mu6L5YkcTcnkiSEt3O/UXPIGhNSCDreWqL05e+dwKuMUt7S6pUT1qLJXIQO6iIxybVjqBvwsInNdl3oDm0VkEzAduM8Yc8p17X5gIhCPNXKfU8bdVuXoxisasuNwCpuTCkiY1WK4ldtk/Re2tT2syTDSc9JZdMD+3C7J6dm893sC/ZrX4gp3OVsObYSE+dD1fvAPKnZ7xhgm7ZhETPUYutQpfu50VT4qZEA3xswwxkQZYwKNMbWNMYNdz39njGltjLnMGNPRGPNjrtesdU3ZNDXGPGgq0pEyqtRd074eQf6+TFmz333BgGBrCeP2mZB6wpa2O0Z2JDIoslSmXT5clEByejaPDy5gHfjSCdYvqsvvKlF7646uY+epndzc8mbN21IJVciArlRRVaviz4jL6vLDhkOcSStgo0/nO8CRBRvsOTXR18eXwU0Gs+TgEltT6h5LyeDTZXu5pn09WtWr5qbgDusX1OV3QZWwErU5acckwgLDGB49vET1qPKhAV15jNu6NyE928GUNQfcF4xsAY17wdrPbDuebliTYeQ4c2xdk/7Wgt3kOAyPDWzmvuDvr1q537s/VKL2Dp47yIIDC7gu9jqC/Io/baPKj1cGdE2f65la1atG96Y1+WJ5ItmOAlaxdL4DzuyD+Pm2tN26ZmuahDXhh/gfbKkv/tg5pqw+wE1XNKRRzZD8Cx7dBtt/sDIqBruZYy+EKTunIAg3Nr+xRPWo8uOVAV3T53quO3o04XByRsGHSLcYDiGRsPYTW9oVEUbFjGLj8Y3sSd5T4vpe+nk7QQG+PHJlrPuCv79qzZ2XcJv/2ayzTI+bzoCGA6gbWrdEdany45UBXXmu/i0iaVwzmE+XFbAb1C/AOvghbi6cTLCl7RFNR+ArviUepS/cdYzfdx3n4QGx1Ax1s0vz8GbYMcta2VLC0fn0uOmcyz7HHW3vKFE9qnxpQFcexcdHuL1HEzbsP8OqPQXkd7n8LmsDzgp7UvlHBEXQO6o3s+Jnke0sXiqAbIeTl37aTpOIEMZ2a+y+8KLXrNOYuo4vVlvnZTmy+Hr713Sp24XWNfM6kExVFhrQlce5oXMDIkIDeXtBvPuCVetY6QA2ToJzx21pe1TMKE5mnGRpUvEOvvh65T4Sjqfy9LCWBPi5+fE8uA52/gTdxkNQ9eJ11uXnPT9zLP0Yd7TW0XllpwFdeZygAF/u7R3N0vgTrNt3yn3hbg9BTgas+diWtntG9aRmlZrMiJ9R5NeeOJfJf+ftpldsBANaRuZf0Bj49TkIjijx6NxpnHy27TNahLegW71uJapLlT8N6Moj3dK1IeEhAbw5v4BReq1m0HwYrP4IslJL3K6/jz9XN72axUmLOZZ2rEiv/efPO0jLyuH5Ea3cb+qJmwv7lkLfJ6GKm/XphfD7gd/Zm7yX21vfrhuJPIBXBnRdtuj5ggP8uLtXNIvjjrN+/2n3hXs8DOmnYb09G41GNxuN0ziZtmtaoV+zLP4E3284yP19mhITWTX/go4cmPe8ldu9020l6qcxhk+2fkL90PoMajyoRHWpisErA7ouW/QOY7o1omZIAK/O3uk+aVfDrtCwu7V9PrvkCbYaVGtA76jefBv3baGOp8vIdvD0jC00rhnM+H4x7gtvnATHd8KVL5QooyLA8kPL2Xx8M3e0uQM/H78S1aUqBq8M6Mo7hAb68cjAZqxOPMXcbUfdF+7/NJw7AmvsWZd+c8ubOZVxirmJcwss+97CeBJPpvHSyLZU8ffNv2BGCix8GaKugJYjStQ/YwzvbXqPOiF1GBUzqkR1qYpDA7ryaDdd3oCYyFBenbODrBw3u0cb94TovtYoPfNcidvtVrcbTcKaMGnHJLd/HWxJSua93xO4tkN9esZGuK/091fh3DEY+pp1AlMJLDu0jM3HN3N327vxL+FIX1UcGtCVR/Pz9eHpYS1JPJnGVyv3uS/c7xlIOwGrPyxxuyLCzS1uZtvJbWw6vinPMhnZDh6dtpGI0ECeH1HA+u+j22HVB9ZmqPodS9Q3Ywzvb3yfuiF1dXTuYTSgK4/Xt3ktesVGMOG3OI4kZ+RfsMHl1sHKS/9rjYRL6OqmV1MtoBqfbv00z+v/mruL+GPn+NfodoQFuxklGwNz/mqtaBnwfIn7teTgEjaf2Mzd7XR07mk0oCuPJyK8NLIN2Q4nz8/a6r7woJcgOw0W/KPE7Qb7B3Nry1tZeGAhcafjLri2LP4Enyzdy9hujegV6+ZYObDOQE1cAgOeK/EWf4fTwYR1E2hQtQEjm44sUV2q4vHKgK7LFr1Po5ohPHxlLHO3HeWXrYfzLxgRa2UuXP8VHM57qqQobm55M8F+wUzcMvF/zx1JzuDhKRuIiQzlyaEFHFxx9gj88gQ06Aodx5W4P7MSZhF/Jp6HOz6so3MP5JUBXZcteqe7e0XTsm41/jZjK8dS3Ey99H4cgmvC7L+W+DDpsMAw/tTiT8xNnMu+lH1kO5w8OHk9aVkOPri1I8EBbpYLGgM/PQo5mXDNu+DjZgVMIaTnpPPOhndoV6sdgxrpunNP5JUBXXknf18f3rqxPWlZOfz52004nfmsPgmqDgNfhAMrbUmvO7bVWPx9/Plo80e8/PMO1u47zavXtXO/gQhgy3TYNRv6PwMRBaxPL4Svtn/FsfRj/LnTn3VXqIfSgK68Smztqjw7vBVLdp/gg8Vu0ua2vwWa9offnofTiSVqMyIoghub38iPCT/x5bqV3NGjCVdfVs/9i07thZ//DFGXlzhfC8Dhc4eZuGUiVza8ko61S7ZKRlVcGtCV17n5ioYMb1eXf83dxbzt+Ww4EoERb1n/nfV/JZ56iQm8BqcjgPrR83n6qpbuC+dkwvTbQYDrJpZ4qgXgtTWvAfD45Y+XuC5VcWlAV15HRPjX9ZfRpl4YD0/ZwPZDKXkXrN7AWvWydxEs+2+x21sWf4InpsVTM3sYZ9jCmqOr3L/gt+fg0Aa45j2o0bjY7Z63OGkx8/fP595291IvtIC/DFSlpgFdeaWgAF8+HtuZqlX8GfPJKuKPnc27YKfboPW1sOAl2Le8yO0siz/BHZ+voUlECN/c+GfqhdTjtdWv5X8AxvovrQ1EXcdDy+FFbu9iqdmp/HPVP4kOi2Zsq7Elrk9VbF4Z0HXZogKoE1aFyXd3wcdHuOnjVcQdzSOoi8CIN62R8rRxcLqA3aa5zNx4kNs/s4L5pLu6UC+sGk9c8QTxZ+L5YtsXl75gzyJrVUvT/jCw5OvgAf615l8cTj3Mi91f1GWKXsArA7ouW1TnRdcKZfJdXQC47r3lLNmdx8lFVarBjZPBkQmTRlupdt1wOA0Tfovj4Skbad+wOt/c3fV/Z4P2b9ifAQ0H8MGmDziQcuCPFyWtg6m3Qs1YGP05+JY8++HipMV8t/s7bmt9G+0j25e4PlXxidu0oh5ORIw3f//qDwfPpHPn52vYfewcD/SL4aH+Mfj7XjTeSVwKX42CupfBLdPzPPrt4Jl0/jxtIyv3nOK6jlH889o2BPpd+KHm0dSjXDPzGlqGt2TioIn4Ht4EX46E4Bpw288QFlXi7+d42nFG/zia8KBwplw1hQDfgBLXqSoGEcEYk+e6Uw3oXvz9qwudzcjmuZnbmLHhIK3rVeNvw1rSI+aiDIg7foJvb4PareDWGRBSE4DUzBwmLtnLB4sS8BF44erWXN8pKt/13j/E/8Czy57lwYbDuHflZOuXw22zrQ9iSyjbkc2dv97JzlM7mTRsErE1Yktcp6o4Kl1AF5F/ASOALCABuN0Yc8Z17SngTsAB/J8xZq7r+U7A50AQMBt4uKBorQFd5WXOlsP846ftHErOoHOjGtxweQP6NY+kVlVr2oS4X2HqrZiqtdnW8z2+OxzO9LVJnM3M4aq2dXlyaAsahAe7bcM4HDw5czS/pMQxMTOEy2/8zpaRuTGGl1e9zNRdU/lX738xpMmQEtepKpbKGNAHAQuMMTki8hqAMeYJEWkFfANcAdQD5gHNjDEOEVkNPAysxArobxlj5hTQjgZ0laeMbAffrN7PVyv2seeEddZo/epB1A2rgo+PUOfsNp46+zI1OMtbjus51OpOxvWMoUPDGgVXfmoP/Pxnzu1ZyM2NYzgZEMCkYZNpHNa4xP2euGUib65/k9ta38afO/+5xPWpiqfSBfTcRGQUcL0x5hbX6BxjzCuua3OBF4BEYKExpoXr+ZuAvsaYewuoWwO6cssYw5aDyaxIOMmOwykcTcnEaQyhgX60Csvk1uMTqH1oHkQ0gx6PQJvrwL9K3pUd3wVrJsLaz6zj4wa9xIHmg7h1zq2E+Ifw2eDPqB1Su9h9nbpzKi+teonh0cN5uefL+IhXrnnweJU9oP8ITDXGfC0i7wArjTFfu659AszBCuivGmOudD3fC3jCGON2Ia8GdGWLXb/A/Bfh2HbwD4HoPhDZykp168iGM/vhwCo4uhV8/KD9zdD3b1CtLgCbj2/mnt/uoXpgdSYOmkhU1aJNvRhj+HjLx7y94W36RPVhQr8J+PvoEkVP5S6gl9vJsCIyD6iTx6WnjTEzXWWeBnKASedflkd54+b5vNq9B9Ali8o+zYdAs8Gw53fYPtNaDRP3CxhXuoDAalCvg7XrtN2NEHph/vN2tdoxcdBE7v3tXm6ZfQuv9nqVbvW6Farpc1nneHHFi/yS+AvDo4fz9x5/12DuxSrsCF1ExgH3AQOMMWmu53TKRVUOOVnWQRk+vhBYQFZFlz1n9vDY74+xJ3kPN7a4kfGXjad6lep5ljXGMGfvHCasn8DxtOM80P4B7mx7p06zeIFKN+UiIkOAN4A+xpjjuZ5vDUzmjw9F5wOxrg9F1wAPAauwPhR92xgzu4B2NKCrCiUtO4031r3Bt3HfEugbyODGg+lRvwdNqjXB38efo2lH2XhsIz/t+Yn9Z/fTIrwFT3d5WjcOeZHKGNDjgUDgpOuplcaY+1zXngbuwJqKeeT8ShYR6cwfyxbnAA/pskVVWcWfjufrHV8zZ+8c0nLSLrgmCB0iO3BD8xsY0ngIvjZkY1SVR6UL6GVFA7qq6LKd2cSdiuNw6mGyHFmEB4XTMrwlYYFh5d01VU40oOdDA7pSqrJxF9D1ExSllPIQ5bZssTyJyAis1AJKKeUxdMrFi79/pVTlo1MuSinlBTSgK6WUh9CArpRSHkIDulJKeQgN6Eop5SF02aJSSnkIXbboxd+/Uqry0WWLSinlBTSgK6WUh9CArpRSHkIDulJKeQgN6Eop5SF02aJSSnkIXbboxd+/Uqry0WWLSinlBTSgK6WUh9CArpRSHkIDulJKeQgN6Eop5SE0oCullIfQdehKKeUhdB26F3//SqnKR9ehK6WUF9CArpRSHqJCBnQR+ZeI7BSRzSIyQ0Squ55vLCLpIrLR9fgg12s6icgWEYkXkbdEJM8/SZRSylNVyIAO/Aa0Mca0A+KAp3JdSzDGtHc97sv1/PvAPUCs6zGkzHqrlFIVQIUM6MaYX40xOa4vVwJR7sqLSF2gmjFmhetTzi+BkaXbS6WUqlgqw7LFO4Cpub5uIiIbgBTgGWPMEqA+kJSrTJLruUuIyD1YI/nzX9veYaWUKhfGmHJ5APOArXk8rslV5mlgBn8srwwEarr+3Qk4AFQDLgfm5XpdL+DHQvRhbQn6/1FleF0J2yzW/Smnvur9qXht6v0phfvj7lFuI3RjzJXurovIOGA4MMC4vntjTCaQ6fr3OhFJAJphjchzT8tEAYdKo9+5/FhJXlfS15Z1e3p/Sue1en8qXpu2q5Abi0RkCPAG0McYczzX87WAU8YYh4hEA0uAtsaYUyKyBngIWAXMBt42xswuoJ21xpjOpfaNVHJ6f9zT++Oe3h/3SuP+VNQ59Hewpld+c81xrzTWipbewN9FJAdwAPcZY065XnM/8DkQBMxxPQrykc399jR6f9zT++Oe3h/3bL8/FXKErpRSqugq5LJFpZRSRacBXSmlPITXBnQRGSIiu1ypAp4s7/5UBCKS6EqfsFFE1rqeCxeR30Rkt+u/Ncq7n2VFRD4VkWMisjXXc/neDxF5yvV+2iUig8un12Unn/vzgogczJWeY1iua15zf0SkgYgsFJEdIrJNRB52PV+67x+710FWhgfgCyQA0UAAsAloVd79Ku8HkAhEXPTc68CTrn8/CbxW3v0sw/vRG+gIbC3ofgCtXO+jQKCJ6/3lW97fQzncnxeAv+RR1qvuD1AX6Oj6d1WsFCatSvv9460j9CuAeGPMHmNMFjAFuKac+1RRXQN84fr3F3hRSgVjzGLg1EVP53c/rgGmGGMyjTF7gXis95nHyuf+5Mer7o8x5rAxZr3r32eBHVi710v1/eOtAb0+1i7T8/JNFeBlDPCriKxzpUgAqG2MOQzWmxSILLfeVQz53Q99T/3hQVem1E9zTSl47f0RkcZAB6w9MqX6/vHWgJ5XAhddvwk9jDEdgaHAAyLSu7w7VInoe8ryPtAUaA8cBv7jet4r74+IhALfAY8YY1LcFc3juSLfH28N6ElAg1xfl0WqgArPGHPI9d9jWDl0rgCOurJZns9qeaz8elgh5Hc/9D0FGGOOGmMcxhgn8DF/TBt43f0REX+sYD7JGPO96+lSff94a0BfA8SKSBMRCQBuBGaVc5/KlYiEiEjV8/8GBmElS5sFjHMVGwfMLJ8eVhj53Y9ZwI0iEigiTbBy8q8uh/6Vq/PBymUU1nsIvOz+uA7Y+QTYYYx5I9elUn3/VNSt/6XKGJMjIg8Cc7FWvHxqjNlWzt0qb7WBGa5UC37AZGPML64cOdNE5E5gPzC6HPtYpkTkG6AvECEiScDzwKvkcT+MMdtEZBqwHcgBHjDGOMql42Ukn/vTV0TaY00XJAL3glfenx7AGGCLiGx0Pfc3Svn9o1v/lVLKQ3jrlItSSnkcDehKKeUhNKArpZSH0ICulFIeQgO6Ukp5CA3oSinlITSgK6WUh9CArlQhiEioK8/35eXdFwAR+VBE/l3e/VAViwZ0pQrnCWCtMWbN+SdE5HMRMSLy3cWFRWSk61rOReXn5VW5q+ytRejP34H7RSS6CK9RHk4DulIFEJEqwP3Ah3lc3g+MEJHaFz1/D7CvtPpkjDkIzAfGl1YbqvLRgK68joh8KyJzc31dS0TSRKRTPi8ZAgQBv+ZxbTewErgtV30NgYHAZ8XsX1/XiP3iR+JFRWcARRnVKw+nAV15owtSlRpjjmMdPjAyn/J9gA3GmJx8rn8E3OXKsAdwF9boubgj9OVYR5idf7TGSqW68KJyq4DaItKymO0oD6MBXXmji3NPA5wj/9OYmgAH3dQ3HQjHyjToC9yBFeTz0ldEzl38yF3AGJNljDlijDkCnATeBfYA9+XxfYB1Nq5S3pk+V3m9JCBURMKMMcmu/O89gQfyKR8EJOdXmTEmQ0S+Au7GOhDYD/gRuCWP4qv4Ix92brvzqf59rF8+XY0xmRddy8jVP6U0oCuvdH5k2wArUL8CHAG+z6f8cawRuDsfAhuAhsBnxpjsP2ZgLpBujIm/+Mm8yorIX4FrgW7GmBN51HW+T8cL6JvyEhrQlTc6H9CjXOvKxwA9jTEZ+ZRfDzzorkJjzA7XYSA9yHsEXiQiMhJraeIQY8yufIq1BRxYv0iU0oCuvNJBwAk8jHXm5aACTqyaA/xHRBoYYw64KTcYqGKMOVWSzolIa+Br4AVgp4jUcV1yuD7APa8vsLSAw4eVF9EPRZXXca1WOQpcBgzIvVkon/I7gN+xRvLuyqWVNJi7XA6EYE0FHc71yL2pSYCbyXttvPJSegSdUoUgIr2AKUCsMSatAvTnBuBZoL2Hn82pikBH6EoVgjFmCfAi1hLGiiAQuF2DucpNR+hKKeUhdISulFIeQgO6Ukp5CA3oSinlITSgK6WUh9CArpRSHkIDulJKeYj/B1we/aei5tNqAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "dpl = ares.util.ParameterBundle('mirocha2017:dpl')\n", " \n", "ax = None\n", "for model in ['floor', 'dpl', 'steep']:\n", " pars = dpl + ares.util.ParameterBundle('mirocha2017:%s' % model)\n", " sim = ares.simulations.Global21cm(verbose=False, progress_bar=False, **pars)\n", " sim.run()\n", " ax, zax = sim.GlobalSignature(ax=ax)" ] }, { "cell_type": "markdown", "id": "88cc3ee6", "metadata": {}, "source": [ "For more thorough parameter space explorations, you might want to consider using the ``ModelGrid`` (see [the example](example_grid) or ``ModelSample`` (:doc:`example_mc_sampling`) machinery. If you'd like to do some forecasting or fitting with these models, check out :doc:`example_mcmc_gs` and :doc:`example_mcmc_lf`.\n", "\n", "**NOTE:** Notice that the ``floor`` and ``steep`` options are defined *relative* to the ``dpl`` model, i.e., they only contain the parameters that are different from the ``dpl`` model, which is why we updated the parameter dictionary rather than creating a new one just with the ``steep`` or ``floor`` parameters.\n" ] }, { "cell_type": "markdown", "id": "d29e0d30", "metadata": {}, "source": [ "### [Furlanetto et al. (2017)](https://arxiv.org/abs/1611.01169) (``furlanetto2017``)" ] }, { "cell_type": "markdown", "id": "e48d30c7", "metadata": {}, "source": [ "The main options in this model are whether to use momentum-driven or energy-driven feedback, what are accessible separately via, e.g., " ] }, { "cell_type": "code", "execution_count": 11, "id": "8808211f", "metadata": {}, "outputs": [], "source": [ "E = ares.util.ParameterBundle('furlanetto2017:energy')\n", "p = ares.util.ParameterBundle('furlanetto2017:momentum')\n", "fshock = ares.util.ParameterBundle('furlanetto2017:fshock')" ] }, { "cell_type": "markdown", "id": "59576548", "metadata": {}, "source": [ "The only difference is the assumed slope of the star formation efficiency in low-mass halos, which is defined in the parameter ``'pq_func_par2[0]'``, i.e., the third parameter (``par2``) of the first parameterized quantity (``[0]``), in addition to a power-law index that describes the rate of redshift evolution, ``pq_func_par2[1]``.\n", "\n", "To do a quick comparison, you could simply do: " ] }, { "cell_type": "code", "execution_count": 12, "id": "87c3ab7e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$f_{\\\\ast}$')" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEWCAYAAACQdqdGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtZElEQVR4nO3dd3hUVf7H8fdJo/cikFBD7yV0EFAsq6CIrIiuBVFEsaG7imXXdVd/lnV1V7FhxUIHBSs2EKVIDb0mlCSEEkqAFNLO748b1xgTSJnkTvm8nicPM2cmd74zJB8O5557jrHWIiIivinI7QJERKTkFOIiIj5MIS4i4sMU4iIiPkwhLiLiwxTiIiI+TCEuIuLDFOIiIj7M50PcGNPCGPO2MWau27WIiJQ3V0PcGPOOMeawMWZzvvZLjTE7jDG7jTGTz3YMa22stXZc2VYqIuKdQlx+/feAKcD7vzQYY4KBV4CLgHhgtTFmIRAMPJ3v+2+x1h4un1JFRLyPqyFurV1qjGmWr7kXsNtaGwtgjJkJXGmtfRoYVs4lioh4Nbd74gUJB+Ly3I8Hehf2ZGNMHeApoJsx5uHcsM//nPHAeIAqVar0aNu2rWcrFhEpY2vXrk2y1tbL3+6NIW4KaCt0qUVr7VFgwtkOaK2dCkwFiIqKsmvWrClVgSIi5c0Ys6+gdm+cnRIPNM5zPwI44FItIiJezRtDfDXQyhjT3BgTBlwLLHS5JhERr+T2FMMZwAqgjTEm3hgzzlqbBdwFLAK2AbOttVs88FrDjTFTk5OTS3soERGvYQJtZx+NiYuILzLGrLXWRuVv98bhFBERKSKFuIiID1OIi4j4MG+cJ14mjDHDgeEtW7Z0uxQRr5CTY4mOP8GizQf5eushkk6foV2D6rRvVJ32DZ0/W9avSsXQYLdLlbPQiU2RAJKVncPK2GN8tSWRr7cc4vCpM4QEwbjwOLqG7iPxVA77k7OYm9GH01QmMugQfRpa2rXtQN8u7YmsX93ttxCwCjuxGTA9cZFAtv3gSeavS+CT9QkcPnWGiqFBDG5dn2sbHWLg9icJPpJnFm8QjBs/no2nq1N11YsMin8DlkHGT8EcCKpHRtVGJF3+Dt1aNyM4aQecOQW1m0PlOmAKuuBaypJCXMRPHTl1hgXRCcxfl8DWxJOEBBkGt6nP6I5VGNi4AhXrt4BjsbAzCK58BdoOA5sDWWdoXKUejYNDoNkkODiU44kxxO/dyalDewhNTmT0e5upUy2G16pPI+po7rV4YdWcMK/XFkZOdQL9ZCJUqAYVqrr7YfgxDaeI+BFrLWv3HWfain18tTmRzGxLl4gajOwewRVNs6i18U1Y9z40Px+um/XLNxWrB336TBaLtx/mi02J7Ni+mWY5++lY6Sj9a5+ibdgRqodmY27+3Hnyh1fD7m+hRmMn3Ou3hfAe0OGqMnj3/q2w4RSFuIgfSMvIZkF0AtNW7GNb4kmqVQxhVI8IruvVhFbEwbL/wKa5Tlh3ugb63QXndSj166acyeK77Yf5fOMBFu84QkZWDs3qVGZk9wiu6hZO42Mr4MA6OLIDDm+HpB0Q0QvG5ob8zOshOBQadILzOjl/VmugYZkCBHyI55mdctuuXbvcLkfEI5JOn+Gdn/bw4cp9nEzPom2DatzQtykjuoZTJSzYCcMlz8Ky/0KPm6HvnVAjokxqOZWeyVebDzJ/XQIrYo8C0KtZbUZ2D+fyzg2pVjEUsrMg/QRUqev8D2DuWEhYCyf2/3qg7jfCFS87t3d982uwB7iAD/FfqCcu/uDAiTSmLo1l5ur9nMnK4ZL2DRjbvxm9mtXCxHwPPz4PvSdAhxGQfhJysqBy7XKrL/54KguiDzBvXTyxR1KoHBbMsM4NGd2zCd2b1MTk72mnJ8OhLZC4EepEQquLnPH0F3LX/q/WECKiIKIntLkc6gbeVGGFeC6FuPiyvUkpvLYkhvnr47EWRnQLZ8KgSFrWrQzbP4Mf/w2J0VA9Ai550vWxZ2st6+NOMHt1HAs3HCA1I5vW51VldM8mjOwWTq0qYYV/c9YZSFjnvJ+EtRC/Bo7vgaumQpfRkLQLVk2FJn2gST+o3rDc3pcbFOK5FOLiiw4mp/Pf73Yya3UcIcFBXNuzMePPb0FErcrOE6ZfCzu/hNqRMGASdB4NIWcJSBecPpPFZxsOMGN1HBviThAWEsTwzo24oW9TukTU+H3vvCApSRBSwZnxsv0LmHcrZKY4j9VsCk36woV/LbMhIzcpxHMpxMWXJKdl8sYPMbyzbA/ZOZY/9WnKHYMjqV85BDbPg3bDIawybP8cMtOcnneQ919huS3xJNN/3s/8dfGkZGTTKbwGN/RpyvAujagUVoz6szPh4CbYvxL2r4C4VTBxJVSqBavfhn3LoflAaDYQarfw6ROmCvFcCnHxBemZ2Xy4ch9TFu/mRGomI7o24oGL29C4RihsnAVLn3eGFq58Fbpd73a5JXb6TBYfr0/ggxV72XnoNNUrhjC6Z2Nu7NuMxrUrF/+AeadL/vQirHwNTh9y7lcPh1YXw7AXfTLMFeK5FOLi7RZvP8zjC7ew/1gqA1vV5aFL29KxYTXYMB1+eA5O7IOGXWDQZGjzB58MpPystazac4z3V+7jq80HsdZyUfvzGNu/Ob2b1y7aUEvBB4aju2HvjxD7g3Mx0+gPnMc+vsM52Rs5BJoOgNCKnntDZSDgQ1xTDMXbxR9P5YlPt/LN1kNE1qvC36/owMBWuZubWwvvXOKc7Bv8MLS+xC/CuyAHTqTxwcp9zFi1nxOpmbRvWJ2x/ZtxRddGVAjx0FBRTg58NAr2/gTZZyCkEjQbAD3HOf8weqGAD/FfqCcu3uZMVjZv/biHl7/fhcFwz4WtGNevCWE7FsDyl+D6uVC1PqQcdXqOfhre+aVlZPNJdALvLtvDzkOnqVetAjf1bcr1vZuefVZLcWSkwr5lzlWlu76BXrdBnzucz3rpv5x/LJv294qTxArxXApx8SY/xx7l4fmbiE1K4Q8dG/DY5e0IP/gdfP8UHNkG9dvDVW9Aw85ul+oaay0/7U7izR/3sHTnESqFBvPHqAjGDWhO0zpVPPtiOdnOieHYJTB9NGSlO2vCtLwAWl/qrC9T0Z2VHBXiuRTi4g3SMrJ5btF23lu+l8a1KvOPKzswOLIGvHsZJKyBOq1gyMPQ/ioI0t4tv9h+8CRv/biHBdEJZOVYLmnfgAmDI+nauKbnXywjBfYshR1fws5FcPog3LPemeVyZIcz1bFWM8+/biEU4rkU4uK21XuP8Zc5G9h7NJWb+jZlcu9QKjVo4zy46FGo3w46XwvBWmS0MIdOpjNt+d7/LTfQp0VtJgyKZFDreiU/CXo21sLhrb+uNzP7Rti6AM7rCG0vh3ZXOI+V4VCXQjyXQlzckp6ZzfOLdvD2sj2E16zElAsr0HXnS04v745lHlmQKtCcPpPFjJ/38/ZPezh4Mp22DaoxYVAkwzo3JCS4DP8Hc2wP7PgCtn3mzE/HQsuL4E9zy+wlFeK5FOLihq0HTnL3jHXEHElhYrcw7gueQ+jm2VChOgy411nnJMzD47sBJCMrh0+iE3jjhxhijqTQuHYlJgyK5OruEWW/vdzpI86SByEVoOt1kJUBb10ILQZBh5HQqJtHeugBH+KaYihusNby4cp9/PPzbdSsFMqLI1vTf8FA5+rK3rc7l8iX48JU/i4nx/LttkO8siSGDXEnqF+tArcNbMF1vZtQpUI5DU+dTIRP74GYxZCT6YybdxgJUWOhZpMSHzbgQ/wX6olLeUlOzeSheRtZsmUf9zfaysixf6ZutYqw7VNo1B1qhLtdot+y1rIi5iivLNnNst1HqVEplJv7NeOW/s2pUTm0fIpIO+4Mt2yZ71xoNPZLaNK7xIdTiOdSiEt5WLvvOPdNX0P/lK/5a+WPqZJxBG79HiJ6uF1awFm//zivLI7h222HqFohhBv6NuXWAc2pU7VC+RWRkgSVapdqppFCPJdCXMqStZa3f4xlxaKZPBI2k0i7H8Kj4OJ/QtN+bpcX0LYlnuSVxbv5fFMiFUKCuL53U24/vwX1q3v35fa/UIjnUohLWUnLyObh+Rv5MnovP1d5gGrVqhN80d+h/ZUBc5WlL9h9+DSvLtnNgugDBAcZru3ZmAmDImlUs5LbpZ2VQjyXQlzKwoG4GJZ/9BSPJF/JPRe1584OWQTVifSKy7WlYPuOOhtszF0bT5Ax/DEqgjuHtCTcS8NcIZ5LIS4edeY0cZ8/Q92NUwm22Wy88H2izr/c7aqkGOKPp/Lakhhmr4kDYFSPCO4c3LJkS+GWIYV4LoW4eERODnbDdNK+/DuVM46wOKQ/LcY8T9PI9m5XJiV04EQary2JYdbqOHKsZVSPCCYO8Z4wV4jnUoiLJ2RlZnL4hX4cSsnhq/C7mXjT9VSvWE5T16RMHUxO5/UfYpi+aj85OdYZZvGCnrlCPJdCXErs+F744TlSBj/B3Z/sZcP2XYwa2JWH/tCOoCCduPQ3B5PTeW3Jbmascnrmbod5wIe4rtiUEjtzGn56AZZPwZogHq80mQ+TWvHEFR24oW8zt6uTMpaYnMbrS2KYsSoOi+WPUY25a0jLcp/NEvAh/gv1xKXIrHX2s/zmcTh9kJOtr+aGfX9gV1o1plzXjQvanud2hVKOEpPTeHVxDDNX78dguLZXY+4c3JIGNcpnnrlCPJdCXIplxnVwKpGNnR/h+q9yqBgazLs396RjeA23KxOXxB9P5ZXFMcxZE0dQkOH63k24Y3Ak9auVbZgrxHMpxOWsTh2C7/8B/SdB3ZaQnsxXu05zz8yNNKtbmXfH9vLaecRSvuKOpfLy97uYty6B0GDDjX2bcfv5Lcrscv7CQlxbhogAZGfC8pfh5R6wYRbErwZg9uaT3Dk9mo7h1Zlzez8FuPxP49qVeW5UF767fxCXdWzIWz/GMvC5xfxr0XaSUzPLrQ71xEViFsOXD0LSTmdh/0ufgboteevHWJ78fBsDW9XljRt6UDlMO+1I4XYfPsV/vt3FZxsTqVYhhFsHtuCWAc2o5qGppxpOyaUQl9/56hHY8bkT3q0vxQIvfLOTl7/fzWWdGvDi6K5UCCnjjQXEb2w/eJIXv9nJoi2HqFk5lAmDIrmxb9NSdwIU4rkU4kJmujN00qQPNB/obIhrgiG0Ijk5lic+3cK0FfsYHdWY/xvZiWDNAZcS2BSfzL+/2cGSHUeoW7UCE4dEMqZXkxLvNKQxcRGAXd/Aq31g8ZOw+1unLawKhFYkO8fy0LyNTFuxj/Hnt+CZqxXgUnKdImrw3thezJ3Ql5b1q/DEp1sZ8vwS5q2N9+jraJBPAsPxffDVw86wSZ1WcMMnEDnkfw9n51genLuReeviuffCVtw3tFXZ7JouASeqWW1mju/L8t1JPP/1DmKTTnv0+ApxCQw7v4LYxTD079Bn4m+WiM3Osfxlzgbmr09g0tDW3Du0lXt1it/q17Iu8yLrkJnt2SFshbj4r5jvISMV2g2DqHHQdtjv9rXMys7hz3M28En0AR64qDV3X6gAl7JjjCEsxLP/wwuYMXFjzHBjzNTk5GS3S5GydvIAzBkLH1wFy19yLp8PDikwwO+f7QT4Xy5powAXnxQwIW6t/dRaO75GDV0u7beys2DFKzClJ2z/HIY8CjcuLHBrtKzsHCbN3sDCDQd46NK2TBzS0oWCRUpPwyniP/YuhUWPOBfsXPYc1G5R4NNyciwPztvIpxsOMPkPbZkwKLKcCxXxHIW4+La04xC3ClpfApEXwNivnPnfhcwssdbyt4Wbmb/OOYmpABdfFzDDKeJnrHXWOJnS0xn/TjvutDfte9YAf/rL7Xy4cj+3D2rBPRdqCEV8n3ri4nuSdsPnk2DPUgiPgmEvQqVa5/y2/3y7i6lLY7mpb1MmX9pW88DFLyjExbekJMEbAyEoFC5/AXqMhaBz/4fyjR9i+O93u/hjjwgeH95BAS5+QyEuviFpF9RtBVXqwrD/QIvBUK1oO+t8sHIfT3+5nWGdG/LM1Z21H6b4FY2Ji3dLPQafTHTGvvctd9q6jC5ygH+28QB/W7CZC9vW58XRXbUWivgd9cTFO1kLG2fDoochPRkGTIKGXYt1iJ92JTFpVjRRTWsx5bruhAarzyL+RyEu3mn2jbBtIUT0hOH/hfM6FOvbN8af4PYP1hBZrypv3diTSmFaD1z8k0JcvEd2FgQFO1MEW14Izc+HqFuctmKIPXKase+uplaVMKbd0osalT2zs4qIN9L/L8U7JKyFqYNh4yznfo+boddtxQ7wQyfTueHtVQC8f0svzqtetjuQi7hNPXFx15nTsPgp+Pl1qHpekeZ7FyY5LZMb317FidQMZo7vS4t6VT1YqIh3UoiLe2J/gAV3QfJ+Z6nYoY9DxZItUHYmK5vbP1hDbNJp3hvbi04RWuhMAoNCXNyTkQKhleCWRc56JyVkreWhuRtZGXuM/4zuSv+WdT1YpIh3U4hL+fll2mDaMehzB7S9DFpdBMGlO/H4/Nc7/rcm+Ihu4ef+BhE/ohObUj5O7IePRsHH4521vnNynPZSBvj0n/fzyuIYxvRqzJ2DtSKhBJ6A6YkbY4YDw1u21Mp15SonG1a/Bd8+4dy/9NncWSel7z8s3n6Yvy7YzKDW9fjnlR21HooEpIDpiWtnH5cc2QFfTXbGvCeuhD4Tij1tsCCb4pOZOH0dbRtU45XruxOiqzElQAVMT1zKUXams0lx60vgvPZw22Jo2KXQdb6L68CJNG6ZtppalcN45+aeVK2gH2MJXOq+iGcdWO9ctDP9Gji0xWlr1NVjAZ5yJotx09aQlpHNOzf31MU8EvDUhRHPyEyDJc/A8pehSj0Y/VGx1zs5l+wcy70zo9lx8CTv3NyTNg2qefT4Ir5IIS6ll5MD71wCiRug2w1w8ZNQqabHX+aZL7fx7bZDPHFFBwa3qe/x44v4IoW4lFxGqnOxTlAQ9L3L6YFHDimTl5qxaj9v/riHG/s25aZ+zcrkNUR8kcbEpWRiFsOrvWHDDOd+52vKLMCX7U7ir59s5vzW9fjbsPZl8hoivkohLsWTngwL74YPRkBwGNRuUaYvF3PkNHd8uJbmdasw5bpumkooko+GU6ToYhbDJ3fC6YPQ/14Y/LAznFJGktMyuW3aGkKCg3jn5p5Ur6h1wUXyU4hL0WWmOasMXvshhPco05fKzrHcPWM9+4+lMv22PjSuXblMX0/EVynE5ey2fw7JCdB7fO6CVRdDcNn/2Dz9xTaW7jzC0yM70at57TJ/PRFfpQFGKVjqMZh3K8y8zjl5mZ3ltJdDgM9ZE8dbP+3hpr5NGdOrSZm/nogvU09cfm/bp/DZJEg77ox7D7i/XMIbYO2+4zz68Wb6RdbhMc1EETknhbj81rFYmH2Tc7XlDZ9Ag47l9tKJyWnc/sFaGtSoyCvXdSdUM1FEzkkhLo7EDc4iVbVbwA0fQ9N+pV7ruzjSM7MZ//5a0jKymH5bb2pVCSu31xbxZerqBLqUozBnLLxxPuxb4bS1GFSuAW6t5ZH5m9iUkMx/ru1G6/O0JopIUaknHsj+N/Z9AoY8BhFRrpTx7rK9zF+fwKShrbmo/Xmu1CDiqxTigerTe2Hte9Cgc7mPfee1IuYoT32xjYvan8fdF2jXJZHiUogHqoZdYcijMGBSuQ6d5JVwIo2J09fRrE5lXrimC0FB2l5NpLgU4oEi9ZizTVrzQdDteoga62o56ZnZTPhgLZlZOUy9MYpquqRepEQU4oFg5yJYeA+kJnl8o4aSsNbyyMfOicy3bowisl5Vt0sS8VkKcX+WdgIWPQrRH0L9DnDdLGerNJdNW76X+eucE5lDdSJTpFQU4v4sbpVzyfzAP8OgByGkgtsVsXrvMZ78fBtD2+lEpognKMT9zZlTsH8ltLoIWl8M96yDWs3crgqAwyfTufOjdUTUqsQLo3UiU8QTfP5iH2PMCGPMm8aYBcaYi92ux1V7foTX+sGsP8Hpw06blwR4ZnYOE6ev43R6Fm/cEKW1wUU8xNUQN8a8Y4w5bIzZnK/9UmPMDmPMbmPM5LMdw1r7ibX2NuBmYHQZluu9MlLhiwdh2jAICoEbF0BV79pI+P++2Mbqvcd55upO2qVexIPcHk55D5gCvP9LgzEmGHgFuAiIB1YbYxYCwcDT+b7/FmttbpeTx3K/L7BknYGpgyFpB/S6HYY+DmFV3K7qNxZEJ/Dusr2M7d+MK7uGu12OiF9xNcSttUuNMc3yNfcCdltrYwGMMTOBK621TwPD8h/DGGOAZ4AvrbXryrhk75GTDUHBzsnKqLHO1MHm57td1e9sP3iSyfM20bNZLR65rJ3b5Yj4HW8cEw8H4vLcj89tK8zdwFBglDFmQkFPMMaMN8asMcasOXLkiOcqdcuBaHh9AMR879zvc4dXBvjJ9EwmfLCWqhVDtLSsSBlxezilIAVNWbCFPdla+xLw0tkOaK2dCkwFiIqKKvRYXi87E378Nyz9F1SpR8EflXew1vLn2RuIO57GjNv6UL96RbdLEvFL3hji8UDjPPcjgAMu1eI9Dm+Hj2+HxGjoPBr+8CxUquV2VYV688dYvt56iMcub6c9MkXKkDeG+GqglTGmOZAAXAtc525JXmDfT5AcB9d8AO2vcLuas1q15xjPfrWDP3RswLgBzd0uR8SvuT3FcAawAmhjjIk3xoyz1mYBdwGLgG3AbGvtFjfrdM2xPbDrW+d21Di4a43XB/jhU+ncNX0dTWpX5rlRnXHOO4tIWSl1T9wYE2qtzcx/uyistWMKaf8C+KK0teVljBkODG/Z0gcu9bYW1r4Lix6DyrXhnvXOcrGVvXtYIis7h3tmrOdkeibvj+ullQlFykGpeuK50/uWGWNqGWOuBJ7zTFmeZ6391Fo7vkaNGm6XcnYnD8BHo5wddxr3hFu+cm297+L69zc7WRl7jKdGdKJtg+pulyMSEErVE7fWWmPMq8CzOPO7L/NIVYHq1CF4ta9zAc9lzztDKEG+MS3v262HeG1JDGN6NebqHhFulyMSMEoc4saYbsA9OPPcrseZ2/2kMQZr7S0eqi8wZGc6ve1q58GA+6DtcKjrA8M+ueKOpXL/7Gg6hlfn8eHur1cuEkhK082Lxblkfg6QBJzOvR94l76Xxs5F8FI3OLjJuT9gkk8F+JmsbO6avg4LvHpdDyqGBrtdkkhAKXGIW2uTrbVrgZE4PfKNQERum9cxxgw3xkxNTk52uxTHmVPObjvTr4EK1cH4Zvg9/cV2NsQn869RXWhSp7Lb5YgEnHOGuDFmmTGmwLNUuSc2Y6y1c4AHgToers9jvOrE5r7l8Fp/WP8B9L8Pxi+G89q7XVWxfbEpkfeW72XcgOZc2rGB2+WIBKSijIn3ACoBJ40xW4B+1tpkcE5sAv+Xe/sA8G5ZFepXdn8HxsDYL6FJH7erKZG9SSk8OHcjXRvX5KFL27pdjkjAKspwSiwwPPd2O8D9Pb58UeJG2LfCuT3oIZiwzGcDPD0zmzs/WkdIsOGV67sTFuIbM2hE/FFRfvseB6YYYxJwFqKabIwZZYyJLNvS/ER2Fix9Ht68AL5+1LmQJyQMKvjuDu//+GwrWxNP8sI1XQivWcntckQC2jmHU6y1c4wxS3HW8n4TGADcClQxxpwCooH1wDpr7QdlWGupuHLF5tEY+HgCxK+CDlfB5S84wyg+bEF0AtN/3s+EQZFc0FY71Yu4zTjD2kV8sjEbgAustUeNMa2Bbnm+ulprvf63Oioqyq5Zs6bsX+jwdnhziDP/+7J/Q6dRPh/gsUdOM/zln2jfqDozbutDiNYHFyk3xpi11tqo/O3FutjHWtslz+2dwE5gVunL8yPZWRAcAvXaQN+J0GMs1PD9LcnSM7OZOH09YSFBvDSmmwJcxEvoN9GTtnwMU6IgOd7pdV/wmF8EODjj4NsST/LC6K40rKFxcBFvoRD3hLQTMO82mHOzs1FDdobbFXnUpxsOMP3n/dw+qAVD2tR3uxwRycMbN4XwLbFL4JM74dRBGPwwDHzAZ1YdLIo9SSk8PH8TPZrW4s8Xt3G7HBHJRyFeWhtnQ2gluPUbCO/hdjUelZ6ZzcTc+eAvj+mmjY5FvFDAhLhHpxgmboDgMKjfztnr0gRDmP+tG/LU59vYmniSt2+KopHmg4t4pYDpWnlk7ZTfXLjzmNNWoZpfBviXmxL5YOU+bhvYnAvbef3MUZGAFTA98VI7FutcuBP3868X7vipuGOpPDhvI10a1+Qvl2hdFBFvphAvivi1MG04BIXAyLf84sKdwmRk5XDXjPUATBnTTeuiiHg5hXhRNOgEXcc4GzbU8O+tx57/egcb4k7w6vXdaVzb/4aJRPyNQrwoQsLg8n+7XUWZW7zjMFOXxnJ97yZc1qmh2+WISBHo/8oCwMHkdB6YvYG2Darx12G+t0GFSKBSiAvZOZb7Zq0nLSObKdd11z6ZIj4kYELc6/bY9CIvf7+LlbHH+MeVHWhZ33fXORcJRAET4l61x6YX+Tn2KC99t4sRXRsxqod/n7QV8UcBE+Lye8dTMrh3ZjRNalfmyas6Yfx02qSIP9PslABlreUvczdyNOUM8+/oT9UK+lEQ8UXqiQeo91fs49tth3jo0rZ0itAQk4ivUogHoC0Hknnq821c0LY+4wY0d7scESkFhXiASc3I4u4Z66lZOZR/jeqscXARH6eB0ADz+IIt7ElK4aNbe1OnagW3yxGRUlJPPIAsiE5gztp47hrSkn6Rdd0uR0Q8IGBCPNAv9tl/NJVHP95Mj6a1uPfCVm6XIyIeEjAhHsgX+2Rm53D3zPUYA/+9tish2mZNxG9oTDwAvPDNTjbEneCV67oTUUvLy4r4E3XJ/NxPu5J4/YcYxvRqzOWdtbysiL9RiPuxo6fPMGl2NJH1qvK3YR3cLkdEyoCGU/yUtZY/z9lAclom08b2olKYlpcV8Ufqifupd5ftZfGOIzx6WTvaN6rudjkiUkYU4n5oc0Iyz3y5naHt6nNj36ZulyMiZUgh7mdSM7K4Z+Z6alUJ5blRXXRZvYif05i4n/nHp1vZk5TCh+N6U7tKmNvliEgZU0/cj3yxKZGZq+O4/fxI+rfUZfUigUAh7icSTqQxed5GukTU4IGLW7tdjoiUk4AJcX9eOyU7xzJpZjTZOZaXxnQjVJfViwSMgPlt9+e1U6Z8v5tVe4/xzxEdaVqnitvliEg5CpgQ91dr9h7jv9/tZETXRozsrt3qRQKNQtyHJadlcu/MaMJrVeKfIzq6XY6IuEBTDH2UtZbHPtnMwZPpzJnQl2oVQ90uSURcoJ64j5q/LoFPNxxg0tBWdG9Sy+1yRMQlCnEftDcphb8t2Eyv5rW5Y3BLt8sRERcpxH1MZnYO985cT3CQ4cXRXQkO0mX1IoFMY+I+5j/f7mRDfDKvXt+d8JqV3C5HRFymnrgPWRFzlFeXxDA6qjGXddIuPSKiEPcZJ1IzuH92NM3rVOFvw9u7XY6IeAkNp/gAay0Pz99E0ukzzL+jP1Uq6K9NRBzqifuA2Wvi+HLzQR64uA2dIvxv2QARKTmFuJeLOXKavy/cSr/IOowf2MLtckTEyyjEvVhGVg73zYymQmgQL1zTlSBNJxSRfDS46sX+/c0ONiUk8/qfetCgRkW3yxERL6SeuJdavjuJqUtjGdOrMZd2bOB2OSLipRTiXuh4SgaTZkfTvG4V/jpM0wlFpHABE+K+srOPtZbJ8zdyLCWDl67tRuUwjXiJSOECJsR9ZWefWavjWLTlEH++uA0dw727VhFxX8CEuC+IPXKaJz51phPepumEIlIECnEvkZGVw72aTigixaQBVy/x4rc7NZ1QRIpNPXEvsCLmKK//EMO1PTWdUESKRyHusuTUTO6fHU2zOppOKCLFp+EUF1lreeTjTRw5dYZ5d/TT6oQiUmzqibto3roEPt+UyP0Xt6ZL45pulyMiPkgh7pJ9R1N4fMFmejevze3nR7pdjoj4KIW4C5zNjqO12bGIlJoGYV3w8ve7iY47wctjutFImx2LSCmoJ17O1uw9xpTvdzGyezjDuzRyuxwR8XEK8XJ0Mj2T+2ZFE16rEk9c0cHtckTED2g4pRw9vmALicnpzL69L9Uqhrpdjoj4AfXEy8mC6AQ+Xp/A3Re0pEfTWm6XIyJ+QiFeDuKPp/LYJ5vp3qQmdw1p6XY5IuJHFOJlLDvHcv+sDeTkWP4zuhshwfrIRcRzNCZexl7/IYZVe4/x/B+70KROZbfLERE/o25hGdoYf4IXv9nJ5Z0acnX3cLfLERE/pBAvI6kZWdw3M5q6VSvw1FUdMUZXZYqI52k4pYz887Nt7Dmawke39qZm5TC3yxERP6WeeBn4estBZqzaz/iBLegXWdftckTEjynEPezwyXQmz99E+4bVuf/i1m6XIyJ+TiHuQdZa/jJ3IylnsnhpTFcqhAS7XZKI+DmfD3FjTDtjzOvGmLnGmDvcrGXa8r38sPMIj13ejpb1q7lZiogECFdD3BjzjjHmsDFmc772S40xO4wxu40xk892DGvtNmvtBOAaIKos6z2bnYdO8X9fbueCtvX5U5+mbpUhIgHG7Z74e8CleRuMMcHAK8AfgPbAGGNMe2NMJ2PMZ/m+6ud+zxXAT8B35Vu+40xWNvfMWE+1CiE8e3VnTScUkXLj6hRDa+1SY0yzfM29gN3W2lgAY8xM4Epr7dPAsEKOsxBYaIz5HJhehiUX6PlFO9h+8BRv3xRFvWoVyvvlRSSAeeM88XAgLs/9eKB3YU82xgwGRgIVgC8Kec54YDxAkyZNPFSmY9nuJN78cQ9/6tOEC9ud59Fji4icizeGeEFjEbawJ1trlwBLznZAa+1UYCpAVFRUoccqrhOpGTwwewOR9arw6GXtPXVYEZEic3tMvCDxQOM89yOAAy7VUihrLY98vImjKWf477XdqBSm6YQiUv68McRXA62MMc2NMWHAtcBCl2v6nXnrEvhi00Huv6gNHcNruF2OiAQot6cYzgBWAG2MMfHGmHHW2izgLmARsA2Yba3d4mad+e0/msrjCzbTu3ltxp/fwu1yRCSAuT07ZUwh7V9QyEnKkjLGDAeGt2xZup11srJzuG/WeoKCDC+M7kpwkKYTioh7vHE4pUxYaz+11o6vUaN0Qx+vLolh3f4TPDmiI+E1K3moOhGRkgmYEPeE9fuP89/vdjGiayOu7KpNHkTEfQrxIko5k8WkWdE0qF6RJ67s6HY5IiKAd84T90r//Gwr+46lMvO2PtSoFOp2OSIiQAD1xI0xw40xU5OTk4v9vV9tPsjM1XFMGBRJ7xZ1yqA6EZGSCZgQL+mJzZwcy/Nf76BjeHUmDdUmDyLiXTSccg5BQYbpt/UmLSObsJCA+TdPRHyEQrwI6ler6HYJIiIFUtdSRMSHKcRFRHxYwIR4aWaniIh4q4AJcU9ddi8i4k0CJsRFRPyRQlxExIcpxEVEfJix1mNbTvoEY8wRYF+ephpA8llu//JnXSCpBC+Z95jFeTx/+9nu5681b1t51n2utsLeg6fq9tRnnb/Nn35G8t5W3cV73O3fyabW2nq/a7XWBvQXMPVst/P8uaa0xy/O4/nbz3Y/f61u1X2utsLeg6fq9tRnXZSfC1/9GVHdvvs7WdiXhlPg03PczttW2uMX5/H87We7X1CtbtR9rrbC3oOn6vbUZ52/zZ9+RvLeVt3Fe9zt38kCBdxwSkkZY9ZYa6PcrqO4VHf58cWaQXWXN0/XrZ540U11u4ASUt3lxxdrBtVd3jxat3riIiI+TD1xEREfphAXEfFhCnERER+mEC8hY8xAY8zrxpi3jDHL3a6nKIwxTYwxC40x7xhjJrtdT1EZY9obY2YbY14zxoxyu55zMca0MMa8bYyZm6etijFmmjHmTWPM9W7WV5hC6v5dm7cppO4RuZ/1AmPMxW7WV5BCam6XmylzjTF3FPlgnpx07utfwDvAYWBzvvZLgR3AbmByvsdGALf7Qs3A0F9qBd73lc8aeAAYmHt7obfXm+exuXlu3wAMz709y1fqPlubj9RdC3jbx2oOKk7N5faX4gtfwPlA97x/CUAwEAO0AMKADUD7PI/PBqr7Qs1AHWAx8D0w1lc+a6A+8ArwL2CZt9eb5/G8If4w0DX39nRfqftsbT5S97+B7r5SM3AFsBy4rqivq+GUPKy1S4Fj+Zp7AbuttbHW2gxgJnAlOMMTQLK19mT5VvqrYtY8FnjcWnsBcHn5VvpbxanbWnvYWjsRmEzJ1pwoteL+bBQgHojIvV1uv3ceqNsVpa3bOJ4FvrTWrivbah2e+KyttQuttf2AIg+5KcTPLRyIy3M/PrcNYBzwbrlXdG6F1fwVcI8x5nVgrwt1nUuBdRtjmhljpgLv4/TGvUVh9dbJ/Yy7GWMezn1sPnC1MeY1yujy62Ioct2FvBe3FOfzvhtn+HCUMWZCOdeZV3E+68HGmJeMMW8AXxT1BbTb/bmZAtosgLX28XKupagKrNlauxnw5hODhdW9FxhfzrUURWH1HgUm5GtMwfmfkDcoTt2/a3NRcep+CXipXKo6u+LUvARYUtwXUE/83OKBxnnuRwAHXKqlqHyxZvC9un2t3l+o7vJT5jUrxM9tNdDKGNPcGBMGXAssdLmmc/HFmsH36va1en+hustP2ddcnmecvf0LmAEkApk4/4KOy22/DNiJc5b5Ubfr9PWafbFuX6tXdQdOzVoAS0TEh2k4RUTEhynERUR8mEJcRMSHKcRFRHyYQlxExIcpxEVEfJhCXETEhynERUR8mEJcApYxZo4xxhpjphfw2IO5jx0txvGeMMbML6vjixREV2xKwDLG7Mm9mWqt7ZCnPRzYBiQDW621lxTxeBuBf1lrPyiL44sURD1xCUjGmDpAM+BtoLUxpkKeh18EPsZZA2N1EY8XCbQDPiuL44sURiEugapn7p/TgBycbeAwxgwFLgaeAZoDa4p4vJHAEmvt8TI6vkiBFOISqKKAeGttHLAZ6JK7VOgU4K/8ugb0avjfriszznK8q3B61yU6vkhJaWcfCVQ9+bUXvB7oDDQC0oFXgYeAg9bahNzndM193u8YYxri7KX4x1IcH2PM+Tg7GNUGvgZettZml+ZNiv9TT1wCVRS/9oLX4QxxPAJMzA3OKH471NEFqGuM+cEYc8gYMzzPYyOANXkDubjHN8aMBO4H/o6zSW42zvrUImelEJeAk9tzbsRve8odgLnW2mW5bXlDGJye+Alr7SCcvTJvzPPYSPIMpZTw+A8AY6y1u621x621LwP7jTEXlOa9iv9TiEsg+uWk4y8h+zNQj9zNmI0x9XHGrNfk3g8B6gPP5j4/BDia+1gtYDC/HQ8v7vFrAInW2jRjzDBjzJLc71sCdCrlexU/pzFxCURRwB5r7TEAa20OkJTvcfg1hNsDm/KMT3cGNuXeHgbstNbuLMXxTwINjDEGWAxszG1vC8SW6B1KwNDFPiLnYIy5AehorX0o9/5snJOOPxpjPga2WGsfK+VrPAJUAf5qrc0xxvQBngOGWmszSvkWxI+pJy5ybl347cyUzvzaW14BzP/ddxTf0zgnNpc5HXL2AH9UgMu5qCcuIuLDdGJTRMSHKcRFRHyYQlxExIcpxEVEfJhCXETEhynERUR8mEJcRMSHKcRFRHzY/wP9r/z0TWLTpgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import ares\n", " \n", "ls = ['-', '--']\n", "for i, model in enumerate([E, p]):\n", " pars = model + fshock\n", " pop = ares.populations.GalaxyPopulation(**pars)\n", " M = np.logspace(7, 13)\n", " pl.loglog(M, pop.fstar(z=6, Mh=M), ls=ls[i])\n", " \n", "pl.xlabel(r'$M_h / M_{\\odot}$')\n", "pl.ylabel(r'$f_{\\ast}$')" ] }, { "cell_type": "markdown", "id": "356a4b00", "metadata": {}, "source": [ "### Creating your own" ] }, { "cell_type": "markdown", "id": "b46ea127", "metadata": {}, "source": [ "As with parameter bundles, you can write your own litdata modules without modifying the *ARES* source code. Just create a new ``.py`` file and stick it in one of the following places (searched in this order!):\n", "\n", "* Your current working directory.\n", "* ``$HOME/.ares``\n", "* ``$ARES/input/litdata``\n", "\n", "For example, if I created the following file (``junk_lf.py``; which you'll notice resembles the other LF litdata modules) in my current directory:" ] }, { "cell_type": "code", "execution_count": 13, "id": "bd31aba5", "metadata": {}, "outputs": [], "source": [ "redshifts = [4, 5]\n", "wavelength = 1600.\n", "units = {'phi': 1} # i.e., not abundances not recorded as log10 values\n", "\n", "data = {}\n", "data['lf'] = \\\n", "{\n", " 4: {\n", " 'M': [-23, -22, -21, -20],\n", " 'phi': list(np.random.rand(4) * 1e-4),\n", " 'err': [tuple(np.random.rand(2) * 1e-7) for i in range(4)]\n", " },\n", " 5: {\n", " 'M': [-23, -22, -21, -20],\n", " 'phi': list(np.random.rand(4) * 1e-4),\n", " 'err': [tuple(np.random.rand(2) * 1e-7) for i in range(4)],\n", " }\n", "}" ] }, { "cell_type": "markdown", "id": "f3894008", "metadata": {}, "source": [ "then the built-in plotting routines will automatically find it. For example, you could compare this completely made-up LF with the rest via the commands `obslf = ares.analysis.GalaxyPopulation(); ax = obslf.Plot(z=4, sources='junk_lf')`." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 5 }