{ "cells": [ { "cell_type": "markdown", "id": "3ba84ba1", "metadata": {}, "source": [ "# The `ParameterizedQuantity` Framework" ] }, { "cell_type": "markdown", "id": "172ff4ed", "metadata": {}, "source": [ "The goal of this framework is to make it possible to take any parameter in *ARES* (the ones that take on numerical values, anyways), and convert it from a constant to a function of an arbitrary set of variables. \n", "\n", "One approach is to allow the user to supply a function of their own to each parameter in place of a numerical value. *ARES* permits this functionality for some parameters (e.g., ``pop_sed``), however, it is often advantageous to retain access to the parameters of the user-supplied function, for example to allow these parameters to vary in some fit. This is the primary motivation for the `ParameterizedQuantity` framework.\n", "\n", "We have already seen PQs in use in the [More Realistic Galaxy Populations](example_pop_galaxy) example, which showed how to make the efficiency of star formation a function of halo mass:\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "5e8baf59", "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": "01076048", "metadata": {}, "outputs": [], "source": [ "pars = \\\n", "{\n", " 'pop_sfr_model': 'sfe-func',\n", " 'pop_sed': 'eldridge2009',\n", " \n", " 'pop_fstar': 'pq',\n", " 'pq_func': 'dpl',\n", " 'pq_func_var': 'Mh',\n", " 'pq_func_par0': 0.05,\n", " 'pq_func_par1': 2.8e11,\n", " 'pq_func_par2': 0.51,\n", " 'pq_func_par3': -0.61,\n", " 'pq_func_par4': 1e10, # Halo mass at which fstar is normalized\n", "}\n", " \n", "pop = ares.populations.GalaxyPopulation(**pars)" ] }, { "cell_type": "markdown", "id": "6f7fefac", "metadata": {}, "source": [ "There are three important steps shown above:\n", "\n", " + Setting ``pop_fstar='pq'`` tells *ARES* that this quantity will be represented by a PQ object. With no ID number supplied, it is assumed that parameters with the prefix ``pq`` are used to construct this object.\n", " + The function adopted is a double power-law, ``'dpl'``, which is a four-parameter model. The parameters ``pq_func_par0``, ``pq_func_par1``, etc. store the values of these parameters. The meaning of the parameters is of course different depending on what function we choose -- see [this listing](params_pq) of the available functions and their corresponding parameters.\n", " + The independent variable is halo mass, ``Mh``. Note that this name is important, i.e., just ``M`` will cause an error. This is a convention of the :class:`GalaxyCohort ` class.\n", "\n", "Let's plot it just for a sanity check:" ] }, { "cell_type": "code", "execution_count": 3, "id": "9ea153dc", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAD8CAYAAAA7fRx2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAh8UlEQVR4nO3dfVyO9+IH8M+3J4QhMs8HR23zmjBhnBk7k+UxFRrJbIbM0/hhzLONbWfTHEMoeTqj1WSeRtrhUFaWbIydDduwPKwOwpbSw/f3R1wTZdV939f3uu/78369er3u73VX1+e76rPL9SiklCAiIstzUB2AiMhesHCJiHTCwiUi0gkLl4hIJyxcIiKdsHCJiHTipDrAw9SpU0c2bdpUdQwiojJLS0v7n5TSvaT3DF24TZs2xZEjR1THICIqMyHEudLe4y4FIiKdsHCJiHTCwiUi0gkLl4hIJyxcIiKdsHCJrIiUEtnZ2SgsLFQdhSrA0KeFEREwY8YMfPXVV7hw4QLS09Px+++/w8HBAbVq1YKbmxvc3NxQu3ZtNGjQAAMHDkSPHj1UR6ZSsHCJDCAnJweJiYlo3rw5/vrXvxZ77/Dhw9i/f3+xZYWFhbhy5QquXLlSbHnjxo0fKFwpJYQQlglO5cJdCkQKpaSkYPDgwXBzc0OPHj2wfv36Bz6nUaNGxcbOzs6lfr/nnnvugWUhISEYNmwYUlNTTQ9MJhFGfuKDt7e35JVmZGvy8vKwZcsWLFmyBIcPHy72XocOHR5YlpqaiitXrqBhw4Zo1KgRatasifz8fFy7dg1Xr17VtnRPnjyJyZMno1KlStrX/vLLL2jWrBkKCgoAAE8//TQmTJiAwMBAuLi4WH6ydkgIkSal9C7pPe5SINLJzZs3sXz5cixbtgwXLlx44H1PT0906tTpgV0A7du3f+BznZ2dUbduXdStW1db1q9fvwc+b/v27VrZAkVb1CkpKahfvz5mzpyJ0aNHw8mJNaAbKaVhP9q1ayeJrF1BQYFct26drFevngRQ7MPFxUW+/PLL8ptvvrHY+lNTU+WwYcOki4vLA+tv3bq1TEpKsti67RGAI7KUTuMuBSILu3LlClq0aIGsrCxt2aOPPorXXnsNo0ePxqOPPqpLjl9//RWrV69GeHg4Ll26VOy9kJAQvPfee6hfv74uWWzZw3Yp8KAZkYXVrl0b8+bNAwA0aNAAkZGROHfuHObMmaNb2QJFJT979mz8/PPPWLhwIapUqaK9t3HjRvz73//WLYu94hYukRlJKXHq1Ck89thjxZbn5eXhn//8J0JDQ1GtWjVF6Yo7f/48Jk+ejC1btqBLly44cOAATx8zg4dt4bJwiczk2rVrCA0NxbZt23DkyBE8+eSTqiOVSUJCAurXr/9A3vz8fB5QqwDuUiCysH379sHLywsxMTHIzc1FcHAwcnNzVccqEx8fnxLLtnv37pg3b16xsxzINCxcIhPk5uZi6tSp6N69O9LT07XlHTp0QH5+vsJkppkzZw4OHDiA+fPnw8fHB5cvX1YdySawcIkq6MKFC+jSpQs++OAD3N01V7t2bWzduhURERGoWrWq4oQVk5+fj6+++kob79+/H0899RSOHTumMJVtYOESVUBycjK8vb2LXS77wgsv4Ntvv0X//v3VBTMDJycnxMfHY+7cudpBtEuXLqFLly744osvFKezbixconKKiopCt27dtH9mOzo6IiwsDJ9//rnNnMfq6OiIefPmIT4+HjVq1ABQdKVcz549sXHjRsXprBcLl6gcjh8/jhEjRuD27dsAinYhJCQkYNKkSXBwsL0/Jx8fHyQmJmo30MnPz8ewYcOwaNEiGPkMJ6Oyvd8QIgvy8vLC3Llztdepqakl3qHLlrRq1QrJyclo1aqVtmzmzJkYM2aMVR8YVIGFS1ROc+bMweLFi3Ho0CE0a9ZMdRxdNGrUCImJifj73/+uLdu9ezcyMjIUprI+uhWuEKK5EGKNEOJTvdZJZKrvv/8eN2/eLLbMwcEBkydPNswVY3qpUaMGdu/ejeDgYDRp0gT/+c9/0KBBA9WxrEqZClcIESWEyBBCnLhvua8Q4gchxBkhxPSHfQ8p5U9SyhGmhCXSU1JSEjp16oTAwEBtn629c3FxwYYNG3DkyBG72bo3p7Ju4a4D4HvvAiGEI4DlAHoCaAlgsBCipRCilRBi530fdR/8lkTGtW3bNvj4+CArKwsJCQkYP3686kiG4eDgAHd39weWHzt2DHl5eQoSWY8yFa6U8iCAq/ct7gDgzJ0t19sAogH4SSm/lVL2ue+jzDt6hBCjhBBHhBBHMjMzyzwRInOJiIhAQEAAcnJyABTdZSs0NFRxKmPbv38/OnXqhODgYB5IewhT9uE2BPDLPeP0O8tKJISoLYRYCaCtEGJGaZ8npVwtpfSWUnqX9H9RIktasmQJRo0apT2GvEWLFvjyyy/Rtm1bxcmM69SpU+jduzdu3bqF2NhYDBkyhKVbClMKt6T7uJV6Yp6U8oqUMlRK+Vcp5TsmrJfIIt5//31MmjRJG7dr1w6HDh1C8+bNFaYyPg8PD4wcOVIbx8bGIjQ0lOfplsCUwk0H0PiecSMAF02LQ6TG22+/jWnTpmnjZ555Bvv27Sv2zDAqmRACS5YswYQJE7Rla9aswaxZsxSmMiZTCjcVgIcQopkQwgXAiwC2mycWkT6klJgzZw5mz56tLevWrRt2796NRx55RGEy63K3dF966SVt2aJFi7B06VKFqYynrKeFbQaQDOAxIUS6EGKElDIfwDgA8QD+CyBGSnnSclGJzO/mzZv45JNPtHH37t2xa9cuuzvH1hyEEIiIiEDv3r21Za+//jqio6MVpjIWPvGB7N6FCxfQtWtXeHp6Ii4uDpUrV1YdyaplZ2eje/fuSE5OBlD0SPddu3bBx8dHcTJ9WN0jdoQQfQH0bdGixcjTp0+rjkN24OLFi6hduzYqVaqkOopNuHr1Krp06YLvvvsOAFC1alV8++23dnGxhNU9YkdKuUNKOerubeGIzOncuXMPLGvQoAHL1ozc3NwQHx+Pxo2LjqtPnToVTZs2VRvKAAxZuESWEhYWhscffxy7d+9WHcXmNWrUCPHx8YiIiCh2M3N7ZshdCndxHy6Z04oVKzB27FgARfsVExIS0LVrV8WpyNZY3S4FInOLiorSyhYAOnbsCG/vEv8myMLy8vIQGxurOoYSLFyyeTExMXj11Ve1cceOHbFr1y6rfcijNcvKykLv3r0xaNAgrFixQnUc3bFwyabt3bsXQ4cO1S4zbdu2Lfbs2cOLGhR5++23kZCQAACYMGEC4uPjFSfSlyELVwjRVwix+vr166qjkBVLSUmBv7+/dsvAJ554Anv37kXNmjXVBrNjCxYsQPv27QEABQUFGDRoEP773/8qTqUfQxYuTwsjU508eRK9evVCdnY2AKBJkybYu3cv6tSpoziZfXN1dcW2bdu0h1LeuHEDfn5+yMrKUhtMJ4YsXCJTXLt2DT169MC1a9cAAHXq1MHevXu1P3JSq379+ti5cydcXV0BAKdPn8bQoUO1W2LaMhYu2ZxatWph4sSJAIBq1aphz549eOyxxxSnonu1bt0aa9eu1ca7du3SnoZsy5xUByCyhGnTpqF27dpo3rw52rVrpzoOlWDQoEFIS0vDP/7xDwBFB9Tatm2LgIAAxckshxc+EJEyBQUF6NWrF/bu3Qug6F8khw8fRsuWLRUnqzhe+EA2TUqJlStXIjc3V3UUKidHR0ds3rxZe6pGkyZN4OzsrDiV5RiycHlaGJXHm2++iTFjxsDX19dujnbbEjc3N2zduhXBwcFISUmBh4eH6kgWw10KZNWWLVtW7BHmCxYsKPb0BiK9cZcC2aS4uLhiz9Hq06cPZswo9YHQZIVu376tOoJZsXDJKiUlJWHIkCHaJbsdO3ZEdHQ0nJx44o2tWL9+PTw9PXHhwgXVUcyGhUtW59SpU/Dz89MOknl4eGDnzp28GY0NmTFjBoYPH45z584hKChIuzzb2rFwyapkZmaiZ8+euHr1KgCgbt262LNnDy/ZtTG+vr5wcCiqp0OHDmH69OmKE5kHC5esxq1bt9CvXz/89NNPAIAqVapg586d2ilFZDu6du2KRYsWaeOwsDDExcUpTGQeLFyyGlOmTEFKSgqAokdyb968WbvzFNmeqVOnol+/ftr45ZdfhrU/VJaFS1Zj5syZaNOmDQBgyZIl8PPzUxuILMrBwQHr1q3TnvR748YNDBw4EDk5OYqTVRwLl6xGgwYNcPDgQURGRhY7HYxsV61atfDpp59qT1Q+duwYpkyZojhVxRmycHmlGZWmevXqGDFihOoYpKOnnnoKYWFh2nj58uXYsmWLwkQVZ8jC5Q3ICQCOHz+u3dSE7NuYMWMQGBiojUeOHAlr3CAzZOESXb58GX369EGvXr2wcuVK1XFIMSEEIiMj0bRpU9StWxfR0dGwxg0y3kuBDCc7OxvdunVDamoqgKLdCD/88APq16+vOBmpdvz4cbi7uxv6d+Fh91LgdZBkKIWFhRg2bJhWtg4ODoiJiTH0Hxjpx8vLS3UEk3CXAhnKrFmzih0QWbp0KXx9fRUmIqPLyMjAxYsXVccoExYuGcbatWvxzjvvaOMJEyZg7NixChOR0e3btw+tW7fGkCFDUFBQoDrOn2LhkiEcPHgQo0eP1sa9e/cudioQ0f3S09Ph6+uLy5cv48CBA3j33XdVR/pTLFxS7scff0RAQIB2R6hWrVph8+bNcHR0VJyMjKxRo0Z48803tfHcuXORnJysMNGfY+GSUrm5uejXrx+uXLkCoOjuXzt27ED16tUVJyNrMGvWLPztb38DUPRAyiFDhhj6/FxDFi6vNLMflSpVwtSpU+Hs7IxKlSph27Zt+Mtf/qI6FlkJJycnfPzxx9o5uWfPnkVoaCiMerorz8MlQ0hKSsKlS5cwcOBA1VHICsXGxmLQoEHaeO3atRg+fLiSLA87D5eFS0Q2YeTIkYiMjAQAVK1aFUePHoWnp6fuOfgQSTKUtLQ03Lx5U3UMsjFLlizB448/DgD4/fffMXjwYMM9hJKFS7o6ffo0fHx80LlzZ/z888+q45ANqVq1KjZv3gwXFxcAwNGjRxEdHa04VXEsXNLN9evX0a9fP1y7dg0nTpxA//79UVhYqDoW2ZA2bdrgvffeg6urKyIjIxESEqI6UjG8lwLpoqCgAC+++CK+//57AEDlypURERGhPSiQyFwmTJgAf39/Q57twt920sW0adOwZ88ebRwVFYUOHTooTES2ysHBwZBlC7BwSQdRUVHFLtOdOXMmBg8erDAR2Ztbt27h8OHDqmOwcMmykpKSEBoaqo379++PBQsWKExE9ub48eNo3749unfvjh9//FFpFhYuWcy5c+eK3SPBy8sLGzdu5H5b0k1hYSFCQkJw8uRJ/PbbbwgJCUF+fr6yPPzNJ4v4/fff4efnh8zMTACAu7s7tm/fjmrVqilORvbEwcEBa9asgZNT0fkBycnJWLhwobo8ytZMNi0vLw/u7u4AAGdnZ8TFxRn2QAbZNm9v72K7sd566y1l+3MNeWmvEKIvgL4tWrQYefr0adVxqILy8/MxZcoUPPnkk3j11VdVxyE7VlBQgOeeew6JiYkAAA8PD3z99deoWrWq2dfFeykQkd07e/YsvLy8tMvKR48ebZEnQvNeCqQL3h+BjKxp06b46KOPtPGqVauwc+dOXTOwcMksMjMz4eXlhdmzZ/NyXTKsYcOGITAwUBuPGDECGRkZuq2fhUsmy8vLw4ABA3D27Fm8/fbbyu5DSvRnhBBYuXIl6tWrB6Doib96PjuPhUsmmzhxIg4ePAig6Bf63htBExlNnTp1sHbtWjg5OWHevHl46623dFs3b15DJlm1ahXCw8O18cKFC9GnTx+FiYj+nK+vL86cOaP7qYrcwqUKS0xMxLhx47RxUFAQpk+frjARUdmpOC+chUsVcv78eQQGBmqXSbZt2xZRUVEQQihORlQxhYWFsPRpqCxcKrfs7Gz079+/2GW7n332GVxdXRUnI6qY8+fPa08i+eabbyy2HhYulYuUEiNGjMDXX38NoOiy3S1btqBJkyaKkxFV3MiRI7Fv3z7k5eUhJCQEOTk5FlkPC5fK5fTp09i1a5c2XrZsGbp06aIwEZHpPvroI1SpUgUAcOLECcyePdsi62HhUrl4enoiJSUFHh4eCA0NxahRo1RHIjKZp6cnPvjgA228ePFiHDhwwOzr4b0UqEKysrLg6uqqPSGVyNpJKdGzZ0/Ex8cDKDqL4fjx43jkkUfK9X14LwUyu5o1a7JsyaYIIRAVFYVatWoBANLT07F//36zroOFSw9VUFCA0NBQix65JTKKBg0aIDw8HJ6envjyyy/h5+dn1u/PwqWHmjNnDlatWoXOnTvjk08+UR2HyOKCgoJw/PhxizxV2pCFK4ToK4RYff36ddVR7FpMTAwWLVoEoOipp0ePHlWciEgflSpVssj35UEzKtGxY8fQuXNnZGdnAwB69uyJHTt2wNHRUXEyImPjQTMql//973/o37+/VrYeHh7YtGkTy5bIRCxcKiY/Px9BQUE4e/YsAKB69erYtm0batasqTQXkS1g4VIxU6dOxb59+7Txv/71LzzxxBMKExHZDhYuaTZs2IAlS5Zo4/nz56Nfv37qAhHZGBYuAQBOnTpV7DJdf39/zJo1S2EiItvDwiUARQfG5s2bByEEWrZsifXr18PBgb8eRObER+wQgKLLGqdPn47WrVvDw8MD1atXVx2JyOawcKmYnj17qo5AZLP4b0Y7lpSUhNu3b6uOQWQ3WLh2Kjk5Gc8//zyef/55XL58WXUcIrvAwrVDFy9eRGBgIG7fvo2kpCSMGDFCdSQiu8DCtTO5ubkIDAzEpUuXAABubm5YtmyZ4lRE9oGFa0eklBg3bhxSUlIAAA4ODoiJiUGzZs0UJyOyDyxcO7Jq1SpERkZq4w8++ADPP/+8wkRE9oWFayeSkpIwfvx4bRwcHIzXX39dXSAiO8TCtQPp6ekYMGAA8vPzAQBPPfUUIiIiIIRQnIzIvrBwbVxOTg4CAwPx66+/AgDq1KmDrVu3okqVKoqTEdkfFq6Ny8rK0rZsHR0dERsbiyZNmihORWSfeGmvjatXrx6SkpIwatQotG/fHt26dVMdichusXDtQJUqVbBhwwbVMYjsHncp2KCCgoIHlgkheJCMSDEWro25desWunTpgmXLlsHIT2QmskeGLFwhRF8hxOrr16+rjmJVpJQYM2YMkpOTMX78eIwePVp1JCK6hyELV0q5Q0o5qkaNGqqjWJVly5Zh/fr12rhVq1YK0xDR/QxZuFR+Bw4cwKRJk7TxSy+9hHHjxilMRET3Y+HagPPnz2PgwIHawTJvb2+sXLmSB8mIDIaFa+VycnIQEBCAzMxMAEDdunURFxeHypUrK05GRPdj4VoxKSVCQ0ORlpYGAHByckJsbCwaN26sOBkRlYSFa8WWL19e7CDZhx9+iGeffVZhIiJ6GBaulcrMzMQbb7yhjYcPH46xY8cqTEREf4aFa6Xc3d3xxRdfoH79+vD29kZ4eDgPkhEZHO+lYMU6deqEtLQ0FBQU8CAZkRVg4Vq5+vXrq45ARGXEXQpWJDo6Gj/99JPqGERUQSxcK3Hw4EGEhITA29sbCQkJquMQUQWwcK1Aeno6Bg4ciPz8fFy7dg2zZ89GYWGh6lhEVE4sXIO7eyVZRkYGgKKzE2JjY+HgwB8dkbXhX62BSSnx2muvITU1FcAfzyTjlWRE1omFa2Dh4eFYu3atNg4LC0PXrl0VJiIiU7BwDSopKQkTJ07UxsOGDcP48eMVJiIiU7FwDejChQsYMGCA9njzdu3a8XaLRDaAhWswubm5CAgIwK+//gqg6CBZXFwcqlSpojgZEZmKhWswjo6O6Ny5s/Y6JiYGTZo0UZyKiMyBl/YajJOTEz788EO0a9cON27cQLdu3VRHIiIzYeEa1NChQ1VHICIz4y4FA8jKyuKVY0R2gIWrWG5uLnr16gV/f3/cuHFDdRwisiDuUlBs4sSJSE5OBgB07twZR48ehYuLi+JURGQJ3MJVKCIiAqtWrdLGr7zyCsuWyIaxcBVJSUnBuHHjtPHgwYMxadIkhYmIyNJYuApcunQJAQEBuH37NgCgdevWiIyM5JVkRDaOhauz27dvY+DAgbh06RIAwM3NDVu3boWrq6viZERkaSxcnU2aNAmHDh0CADg4OCA6OhrNmjVTnIqI9MDC1VFUVBRWrFihjd999134+PgoTEREemLh6iQ/Px9hYWHaOCgoCFOmTFGYiIj0xsLViZOTExITE9GzZ0+0atUKa9as4UEyIjvDCx90VKtWLezYsQNXr15F1apVVcchIp1xC1dnjo6OcHd3Vx2DiBRg4VrQpk2bEBMTozoGERkEdylYSFpaGkaMGIGcnBykpaVh4cKFcHLif24ie6bbFq4Qor8QIkIIsU0I0UOv9aqQmZmJgIAA5OTkAAB27typvSYi+1WmwhVCRAkhMoQQJ+5b7iuE+EEIcUYIMf1h30NK+ZmUciSA4QCCKpzY4PLz8xEUFITz588DAGrUqIHPPvsM1apVU5yMiFQr679x1wFYBmDD3QVCCEcAywH4AEgHkCqE2A7AEcA79339K1LKjDuvZ935Opv0xhtvYP/+/QAAIQQ+/vhjeHh4KE5FREZQpsKVUh4UQjS9b3EHAGeklD8BgBAiGoCflPIdAH3u/x6i6KTTdwHsllIeLW1dQohRAEYBsLqHJ27atKnYxQ3z589H7969FSYiIiMxZR9uQwC/3DNOv7OsNOMBdAcwQAgRWtonSSlXSym9pZTe1nT61DfffINXX31VG/v5+WHmzJkKExGR0Zhy2Lyky6RkaZ8spVwKYKkJ6zOsK1euwN/fH7du3QIAPPbYY9iwYQMcHHjWHRH9wZRGSAfQ+J5xIwAXTYtjfaSUCA4OxtmzZwEA1atXx2effYZHHnlEbTAiMhxTCjcVgIcQopkQwgXAiwC2myeW9RBCYPz48ahRowYAYOPGjXj88ccVpyIiIyrTLgUhxGYA3QDUEUKkA5grpVwjhBgHIB5FZyZESSlPWiypgfXu3Rupqan44osv4OfnpzoOERmUkLLU3a7KCCH6AujbokWLkadPn1Ydh4iozIQQaVJK75LeM+RRHSnlDinlqLv/TDeaa9euISMj488/kYjoHoYsXCMrKCjAkCFD0K5dO6SmpqqOQ0RWhIVbTnPnzsWePXuQnp6OZ599Fr/88suffxEREVi45RIXF4eFCxdq48mTJ6Nx48YP+Qoioj+wcMvou+++w0svvaSNfX19sWDBAoWJiMjasHDL4Pr16+jfvz9+++03AEDz5s2xadMmODo6Kk5GRNbEkIUrhOgrhFh9/fp11VFQWFiIoUOH4u7paa6urti6dStq1aqlOBkRWRtDFq6RTgtbsGABdu7cqY2joqLg5eWlMBERWStDFq5RbN++HfPnz9fGU6dORVCQzd47nYgsjIX7EJ9//rn2unv37li0aJHCNERk7fhUw4cIDw+Hp6cnVqxYgejoaD4EkohMYsh7Kdzl7e0tjxw5ojoGcnJyULlyZdUxiMgKWN29FIyGZUtE5sDCvcfu3bsxatQoPtKciCzCkDsl77k9o27rPHPmDIYMGYKsrCwcO3YMcXFxaNjwYY9oIyIqH0Nu4ep9Hu5vv/0Gf39/ZGVlAQAuXrwIZ2dnXdZNRPbDkIWrJyklXnnlFZw4cQIAUKlSJcTFxaFu3bqKkxGRrbH7wn3//fcRGxurjVeuXIn27dsrTEREtsquCzchIQEzZszQxmPHjsXw4cPVBSIim2a3hfvzzz/jxRdfRGFhIQDgmWeeQVhYmOJURGTL7LJws7Oz4e/vj6tXrwIAGjRogNjYWLi4uChORkS2zC4Ld9q0aTh27BgAwNnZGVu2bEG9evUUpyIiW2eXhTtjxgx07NgRALB8+XI8/fTTihMRkT2wywsfGjZsiAMHDuDTTz9FcHCwRdZBRHQ/3ryGiMiM7P7mNbdu3cKhQ4dUxyAiO2fzhSulxJgxY/Dss89i8eLFMPIWPRHZNpsv3BUrVmD9+vUoLCzElClTEB8frzoSEdkpmy7cxMREvP7669p4+PDheOGFF9QFIiK7ZrOFm56ejgEDBiA/Px8A4O3tjfDwcAghFCcjIntlk4Wbm5uLwMBAZGRkAADc3d0RFxfHJzcQkVI2Wbjjxo3DV199BQBwdHRETEwMGjdurDgVEdk7myvc1atXIzIyUhsvXrwY3bp1UxeIiOgOmyrc5ORkjBs3ThsPHToUEyZMUJiIiOgPhixcIURfIcTq69evl+vrzp8/DweHoim1adMGq1at4kEyIjIMQxZuRZ9pFhQUhKSkJLRp0wZbt26Fq6urhRISEZWfIW9eYwpvb28cPXqUW7ZEZDiG3MI1FcuWiIzIJguXiMiIWLhERDph4RIR6YSFS0SkExYuEZFOWLhERDph4RIR6cTQD5EUQmQCyAJw7zW+Ne4Zl/S6DoD/mbjqe79vRT+vpPfuX/awMedWfpxb6Z9X1uWcm+lz+4uU0r3Ed6SUhv4AsLq0cUmvARwx9zor8nklvfewuXBunJsl51bW5Zybeed2/4c17FLY8ZBxaa/Nvc6KfF5J7z1sLvePObfy49xK/7yyLufczDu3Ygy9S6EihBBHZCnPhLd2nJt14tyskyXmZg1buOW1WnUAC+LcrBPnZp3MPjeb28IlIjIqW9zCJSIyJBYuEZFOWLhERDqx+cIVQjQRQmwXQkQJIaarzmNOQoiWQogYIUS4EGKA6jzmIIRoLoRYI4T49J5lVYUQ64UQEUKIYJX5TFHK3B5YZo1KmVv/Oz+zbUKIHirzmaKUuT0hhFgphPhUCDGmzN/M3Cf26vEBIApABoAT9y33BfADgDMApt9Z1h3A6DuvN6jObua5/R+ALndeb1ed3Rxzuue9T+95HQKg753Xn6iejznn9rBlqj/MOLdaANaono+F5uZQnrlZ6xbuOhT9h9EIIRwBLAfQE0BLAIOFEC0BfA3gRSHEPgD7dc5ZEetQ9rltRNHc3gdQW+ec5bEOZZ9TSRoB+OXO6wILZayodTBtbka2DuaZ26w7X2Mk62Di3IQQ/QAkAfh3WVdqlYUrpTwI4Op9izsAOCOl/ElKeRtANAA/AC8DmCul/DuA3vomLb/yzE1KmSGlHAtgOky/nt1iyvnzKkk6ikoXMNjvrBnmZlimzk0UeQ/AbinlUcumLR9z/NyklNullJ0BlHk3l6F+eU3UEH9sBQFFf6QNAewBMEEIsRLAWQW5zKHEuQkhmgohVgPYAOB9JckqrrQ51b7zs2orhJhx5704AIFCiHBY8LJLMyrz3EqZr5GV5+c2HkW79AYIIUJ1zlkR5fm5dRNCLBVCrALweVlXYEuPSS/pUb1SSnkCgLUfUCptbmcBjNI5i7mUNqcrAELvW/g7iv6lYi3KM7cHlhlceea2FMBSXVKZR3nm9h8A/ynvCmxpCzcdQON7xo0AXFSUxdxscW62OKe7ODfrZPG52VLhpgLwEEI0E0K4AHgRwHbFmczFFudmi3O6i3OzTpafm+rTMyp4SsdmAJcA5KHo/0oj7izvBeAUgB8BzFSdk3Oz3TlxbpxbRT548xoiIp3Y0i4FIiJDY+ESEemEhUtEpBMWLhGRTli4REQ6YeESEemEhUtEpBMWLhGRTli4REQ6+X+dPAlFYPVmLAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "Mh = np.logspace(8, 13)\n", "sfe_dpl = pop.fstar(z=6, Mh=Mh)\n", "\n", "pl.loglog(Mh, sfe_dpl, color='k', ls='--', lw=3)" ] }, { "cell_type": "markdown", "id": "1b482bb6", "metadata": {}, "source": [ "## Multi-Variable `ParameterizedQuantity`" ] }, { "cell_type": "markdown", "id": "56251542", "metadata": {}, "source": [ "More complicated models are also available. For example, say we wanted to allow the normalization of the SFE to evolve with redshift, i.e.,\n", "\n", "$$ f_{\\ast}(M_h) = \\frac{2 f_{\\ast,0} \\left(\\frac{1+z}{7}\\right)^{\\gamma_z}} {\\left(\\frac{M_h}{M_{\\text{p}}} \\right)^{\\gamma_{\\text{lo}}} + \\left(\\frac{M_h}{M_{\\text{p}}} \\right)^{\\gamma_{\\text{hi}}}} $$\n" ] }, { "cell_type": "markdown", "id": "e37ae8b6", "metadata": {}, "source": [ "Starting from the pure ``dpl`` model above, we can make a few modifications:" ] }, { "cell_type": "code", "execution_count": 4, "id": "e1726b7a", "metadata": {}, "outputs": [], "source": [ "# Extra multiplicative boost with redshift, par0 * (var / par1)**par2\n", "pars = \\\n", "{\n", " 'pop_sfr_model': 'sfe-func',\n", " 'pop_sed': 'eldridge2009',\n", "\n", " 'pop_fstar': 'pq[0]', # Give it an ID this time, since we'll add another\n", " 'pq_func[0]': 'dpl_evolN', # dpl w/ evolution in the Normalization \n", " 'pq_func_var[0]': 'Mh',\n", " 'pq_func_var2[0]': '1+z', # indicate 1+z as the second indep. variable\n", "\t \n", " # Old parameters that we still need\n", " 'pq_func_par0[0]': 0.05,\n", " 'pq_func_par1[0]': 2.8e11, \n", " 'pq_func_par2[0]': 0.51,\n", " 'pq_func_par3[0]': -0.61,\n", " 'pq_func_par4[0]': 1e10,\n", " 'pq_func_par5[0]': 7., # New param: \"pivot\" redshift\n", " 'pq_func_par6[0]': 1., # New param: PL evolution index\n", "}" ] }, { "cell_type": "markdown", "id": "8eec37f7", "metadata": {}, "source": [ "To verify that this has worked, let's again plot the SFE, now as a function of redshift, and compare to the previous $z$-independent model:\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "a22e18cf", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAD8CAYAAAA7fRx2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABNuUlEQVR4nO3dd1hURxfA4d+l9y4dxV6jWLDEEmM0ltjFXhB7iRpN1JjEmB6TmC8x9o69d2yxa6wo9q7YAAVBel3Y+/1xjQiLBhWWBebN4xPYOzs7I3i4zJ45I8myjCAIgpD/9Ap6AIIgCMWFCLiCIAhaIgKuIAiCloiAKwiCoCUi4AqCIGiJCLiCIAhaYlDQA3gVBwcH2dPTs6CHIQiCkGtnz56NlGW5RE7XdDrgenp6cubMmYIehiAIQq5JknT/ZdfEkoIgCIKWiIArCIKgJSLgCoIgaIlOr+HmRKVSERISQkpKSkEPRStMTExwd3fH0NCwoIciCMJbKnQBNyQkBEtLSzw9PZEkqaCHk69kWSYqKoqQkBBKly5d0MMRBOEtFbolhZSUFOzt7Yt8sAWQJAl7e/ticzcv/DdZlklMTCQjI6OghyK8gUJ3hwsUi2D7r+I0VyHT1atX2bt3L7dv3+bRo0eEPQ4jPCGcp6lPSTdKR06TMdUzxdLQktbvt2bm/2aK75VCoFAGXEEoqmRZJjg2mAmrJnAx5iLGzsYYVDDAwMoAs2f/ZXeEI9RcXpMSZiWoUaIGNR1rcu/YPfav2U+jdxvx4Ycf0rBhQxGQdYAIuIJQQMLDw/H39yfDOIPybcpz+tFpAh8HEpUSBeXBNMqUlPspJN1OQhWjIj0mnfTYdKRkiQwpA31zffTN9WnSognN2jQjJD6Ec0/OsefeHjCEjK4ZLA9ezuzvZmMfYc+g3oPo168fTk5OBT31YksE3AKWkZFBnTp1cHNzIyAgoKCHI+QztVrNvn37mDd/Hvtu7MO6sTVWta2QTkg4mjpS37U+dZ3rUpKSTPt6Gt51vHF1dcXFxeX5H0tLSzIyMoiNjSU6OhozMzNcXV2fv8bjxMc08GlAqn0qZhXNcOntglql5q+zf/FD8x/4oMIHDB40mJYtW6KnV+jexinURMAtYNOnT6dy5crExcUV9FCEfKRSqViyZAlTZ0wlvmQ8tu/ZUrJdSdIT0onaH8UX7b/g464fZ/m1f83qNS/tz8DAAHt7e+zt7TWuOZs7c2T+EY4fP86+ffvYtGoTxrWNsXnXBpv6NlyJuILvbF+cv3bml+9+oVWrVmK5QUvEj7c31KxZM7y8vPDy8sLExIT169e/dh8hISHs2LGDQYMG5cMIBV2gVqtZs2YNVWpXYcrBKZiOMcW5mzPp0ek8nPsQ+/X2/Nr6VwZ1GJSnQc/Dw4Pu3buzYMECHpx5wE/Nf8J2rS0P5z4kLTINJx8nVP1U9P2jLxevXMyz1xVerVDf4X67/QpXw/L2zrCKqxVT2lX9z3YHDhwAYM6cORw8eJDOnTs/v9a4cWPi4+M1njNt2jSaN2/+/PNPPvmEX3/9Nce2QuG3d+9ePvviMx65PcJhhAP2xvbEHIsh7Z80erfuzeBFg6lSpUq+j8Pc3Jz+/fvTv39/rl27xoIFC1gydQk27W1w6e3CxBsTGW02mlalW6EniXuw/FSoA25BW7ZsGbt27WLjxo3o6+s/f/zo0aP/+dyAgAAcHR2pXbs2hw4dysdRCgUhNSOV5VeXk9YvDScrJ2LPxJK8N5lxvuMYPX005ubmBTKuypUr87///Y8JEybw7Xff8kHVD1gTtoaJRyfif8Wf8d7jsU+yp2TJkhgYiPCQ1yRdPia9Tp06cvbyjNeuXaNy5coFNKJM69evZ+HChWzduhUTE5Ms13Jzhztp0iSWL1+OgYEBKSkpxMXF0blzZ1asWKHxPF2Zs5A7gY8D+fbEt9yPu0/arTQiNkYwrOMwJkyYgK2t7dt1LsugSgIDE9DT/+/2uaCW1ewI3sHMczMJSwwj7WQa1hetWbZoGZUqVcqT1yhOJEk6K8tynRyviYD7+gICAvj9998JCAjIkzuVQ4cOMW3atJdmKejCnIVXO3nyJPau9qwKXcWGmxtwt3Bncv3JqG6rqFChAi4uLrnvTJYh4hrc+weibkH8Y0gIz/x/+rOdh8ZWYGIDptbK/61cwd0bPOqBU9XXDsjJ6cm0/7U9j90ek/o4lSdLnzBlyBRGjx4tshlew6sCrvid4Q34+vpiZ2dHw4YNARg1ahQDBw4s4FEJBSEhIYEvv/wS/2P+lB5cGtlUpn/V/ozwGoGpgSm45aITWYao23D3MNw9qgTapEjlmrE1WDqBhRN41FX+b+4AqhRIiYHkGEiJVT4OPgwX1yrPM7IE9zpQsj5UbAMu1f9zGCb6JjSTmjF12lSc+zvjNt6NqTunsnnbZvwX+Yt6HnlA3OEWAsVxzoXB2bNn6e7XnbQmaVh7W5P8IJmxVcYyssvI3HWQEgsX18HZpRB+SXnMyg08G0PpxuDZCGw9cz8gWYaY+/DgFDw8CQ9PQ/gVQAYXL6jVD97xARPrV3Zz8eJF+g7qS2S1SOzesyPlYQrRK6JZ/sdyWrdunfvxFFNiSaGQK45z1mWyLDNjxgwmz5+MyyAX9C30idgSQa2MWsybM49XnsMnyxByBs76w5VNynqsc3Wo2RfKfQB2ZSAvc2KTnsKl9UpQj7gCBqZQtRPU9lXufl8iLS2N77//nhk7ZuDc3xk9Iz1CF4UyqcskJk6cKPJ2X0EE3EKuOM5ZVz19+hS/AX6c4ASOHRxJC08jyj+KP774g759+748EMky3NgFh36Cx5fAyEK526zdH1xr5v/AZRnCgiBoGVzaCGnxUOZ9aPEtuNR46dNOnTpFtwHdMOhigFlZMyK2RdBEvwn+S/wLLNNC170q4IqVcEHIpWPHjlGzcU0uVrqIUycnYk7EYLHFglMBp+jXr9/Lg+3do7DoQ1jTE9KSoO2f8Ol1aDddO8EWlLtmt9rKa352A1r+BI8uwLwmsHEQRN/L8Wn16tUj8GAgbv+4EX0kGsf2jpx2O83shbO1M+4iRgRcQfgPsiwzffp02oxsg9kwM0xLmxIyPwQfMx9OHD5BuXLlcn5i2DlY3gmWtoXYECXYjTwFdfzA2FK7k3iRkTk0GAljzkOjcXAtAGbUgV2fQ2KURnNHR0f279lPR7OOhC0Pw6q6FUdKHiE4Nlj7Yy/ktBZwJUkqI0nSIkmSNmjrNQUhL6Snp7Pq9io8xnigilER+Uck/hP8+fPPPzE2NtZ8QmIUbBoK85tC2Hn48AcYHaQsH+jr0FFJJtbQfIoyNq+ecHoezKwDV7ZoNDU0NGTmjJnMGjSLGY1nEJcWR+8dvQl8HKj9cRdiuQq4kiQtliQpQpKky9kebyVJ0g1Jkm5LkvT5q/qQZTlYlmWROyUUKinpKXx14iuS6iSRfjUd+z32nPn7DB06dNBsLMtwcT3M8obLG6HxpzDmArw7CgxNtT/43LJyhfYzYNgxsCkJ632VZYakpxpNe/ToQdNyTVnz0RoczRwZtncYe+/tZf/+/QUw8MInt3m4/sBMYNm/D0iSpA/MAloAIUCgJEnbAH3g52zPHyDLcsRbj1YQtCg8MZwxB8dwJeoKo2uOplmjZrh976axsxCAmAcQMA5u7wW3OkoAc8q7OgmqDDVP4lOJiE8lOjENE0N9rEwNsDY1xNrUEAtjg7fPHHCqAoP2wdH/wZFflXzg9jOgfAuNpi4WLixttZSRB0Yy7tA4Qv1DGX10NFOmTBEZDK+Qq4Ary/IRSZI8sz1cF7gty3IwgCRJa4AOsiz/DLTN01EKghatW7eOM6FnCHQLJEmVxPT3p9OsZLOcG6vVcHo+7P9O+bzVL1B38Btvu5VlmdsRCZwIjuJkcBTBTxKJiE/laWLaK5+nryfhZGlMzZK21C5lSx1PWyq7WGGo/5qrhvqG0HQiVGgJm4fBSh+o5au8yWZskaWpjYkNzSKbceziMdz83Ji9cTbhI8OZOWNmltoiQqa32WnmBjx84fMQoN7LGkuSZA/8CNSUJGnSs8CcU7shwBCAkiVLvsXwdJ+npyeWlpbo6+tjYGBA9hQ4QbtkWea3337jp00/4TbIDfsUe1a0W0F52/I5PyEhAjYPhTsHoFxzaPuH8iv5a4pMSGXPlcecuBPFyeCnRCakAuBqbUIVV2tqlbLF0dIYR0sTHC2NsbMwIkWVQVyyithkFXHJ6cQmq7j/NImz956y49IjAEwN9fHysKFFFSc61XTD1two94Ny9YKhh+Hgj3DsLwgJhB4rlTzhF3Tv3J11q9dxJeEKTl2c2LpvK096PGHVylUYGb3G6xUTuc7DfXaHGyDLcrVnn3cFWsqyPOjZ532BurIsj8qrwRX1PFxPT0/OnDmDg4PDK9sVpTnrqoyMDMZ8Moa1d9bi3N2ZxOuJ2B61JfBIYM51BIIPwaYhym6xVj9Dbb/X2rCgVsucCI5i1ekH/H3lMaoMGScrYxqUsadBWXsalHHAw870jX49D4tJ5uz9aM7ej+ZkcBTXH8djpK9Hy2rO9PD2oEEZe/T0XqPfOwdgvZ/ysc8i5YfLC9LS0vDt78sh9SEcWjsQczyGmo9rsmH9hmIZdPOrlkII4PHC5+5A2Fv0V6g0a9aMp0+VNxWuX7/O8uXL6dq1awGPSngTycnJ9Ordi1Nmp3Du7kzsqVg8r3mydftWzWCbkQ6Hfoajv4NDBei7WSkUk0tP4lPZcDaENYEPuB+VhI2ZIX3re9LN252KTpZ5sv7pamOKq40p7Woox+5cDYtj3ZmHbAoKYfuFMDzsTOnhXZJ+DUphaZKLrImyzWDIIVjTG1Z2hQ+mQMMxz3/AGBkZsXLFSsaNG8eqDatw8nHi3Ilz+HT1KbZB92Xe5g7XALgJfACEAoFAL1mWr+TV4P7zDnfX58qunbzk/A60nprr5v8WIF+9evXzdavcFiAvXbo0tra2SJLE0KFDGTJkSI6vIe5w809UVBTtO7fnYfWHWNexJnJ3JI3VjVm2dJnmm2OxobBxIDw4ATX7QOtflZzW3LxOQiqzD91h+cn7pKWrqVfajl71StKyqjMmhtpZ70xRZbDnymPWBj7k+J0o7MyN+Pj9cvSuXxJjg1yMIS0RtoyAq1ugamfoMDPL/GVZZtKkSSy6vAjnrs7EnIih5qPid6f71ne4kiStBpoCDpIkhQBTZFleJEnSx8AelMyExXkZbAuDtylADsrOJVdXVyIiImjRogWVKlWiSZMm+TVcIZt79+7RqmMr0tqkYVXeikerHtGvaj+mTZumeWd77xis6wvpqdB5AVTvlqvXSEhNZ+HRYBYcCSZZlUHnWu4Me68s5Rwt/vvJeczEUJ8OXm508HLjYkgMU3dd57uAqyw5fpfPPqxIu+qur15qMDKHrv5w7E/Y9y1E3oI+G8DSGQBJkvj5559hEixarwRdcaeblU7WUpAkqR3Qrly5coNv3bqV5Zqu3O29bQHy7L755hssLCz47LPPNK7pypyLkuvXr9OicwtM+phg5GhE6PxQpvScwtixYzUbn1kCOz8D29LQczU4vORNtBekqDJYcfI+sw/d4WliGq2qOvPphxUo71SAO8yykWWZI7cimbrrOtcexVHNzYqvPqpC/TKaB1NquLUX1vmCuT303QL2ZbP0++KdbmJgIluGbqH6O/9dIrIoeNUdLrIs6+yf2rVry9ldvXpV4zFt2759u9y0aVM5ISHhjftISEiQ4+Linn/coEEDedeuXTm21YU5FyWJiYmyezV3ucK0CnLl2ZVl62rW8tq1azUbpqfJcsCnsjzFSpaXd5blpOhc9X/iTqTc9LeDcqmJAXLvBSfl8w9y97yCkpGhljcFPZTf/Xm/XGpigPz1lktyYqrqv58YclaWfykty7+UkeXQoCyX1Gq1PHHiRNmts5tczb+aPOHwBDk9Iz2fZqBbgDPyS2KaTt7h/ktXsxTs7e2xs7N7Xi3pTQqQBwcH06lTJ0DZOtqrVy++/PLLHNvqwpyLkpvRN/Hd5ktMfAzhs8LZMGsDzZply7NNjFJ2XN07quwUa/7tf+bWxiarmLrrGqtPP6SknRk/dKxGkwol8nEmeSs5LYNf91xnybF7eNqb8Xu3GtQuZffqJ0XeVupFJD+F7iug7PvPL8myzIMHD9gbt5fpQdPpUr4LUxoU/Y0RojxjIVcc55xfLj25xLB9wzDRN6FtclsaVGxA/frZ6sJGXIdV3ZQjbdpNV+oM/Ifdlx/z9dbLRCakMqhxGcY2r4CpUeFM/j9+J5Lx6y/yKDaZIU3KMrZF+Ve/qRb3CFZ0gcib0Hk+VOus0WTGuRnMvzgfv2p++JXxw8bGpsgGXnHEjlDsqVQqzkWeY9SBUdiZ2LHgwwW4W7prNrx3TCmjqG8MfjuVY2peISohlS83X2b3lcdUdrFika8377i/+kQFXfduWQd2f9KYH3dcY+7hOxy8HsGMXjWp8LL1ZysX5e9qdU/YMACSo8E76298H3t9TGxqLEsuL2H2/2bjW8mXyZMna2E2ukWUZxSKvFWrVlHLpxbD9g3D1cKVpa2X5hxsr2yG5R2Vc8MG7fvPYHviThRt/jrKgRsRTGhVkW0fNyz0wfZfliaGTO1SnSX9vYlKTKPjrGPsvvzo5U8wtYG+m5QtwTvGwekFWS5LkoSPrQ/JQcmYfmjKH/v/YMaMGfk7CR0kAq5QpC1btowR/xuBur2a9EfpTH93Oo5mjpoNT8xWdlO51oIBe8C21Ev7zFDL/LH3Jr0XnsTcyIDNI95lRNNyr1+3oBB4v5IjAaMaUcHJkmErgvj97xuo1S9ZhjQ0hW7LoUJrJasjW9At7VmaMjfKEHc+DldfV75a8RUrVqzQwix0h05+h0iS1E6SpPmxsbEFPRShEFuyZAmjZozCfbg7KfdTkNfJmMjZNjOo1bDnS9gzCSq3hX5bwOzlbxQ9jk2h14KTTN9/i45ebmwf1YiqrkXjrvZlnK1NWDu0Pt3quDPjwG0GLztDXIoq58YGRtBtWY5B19jYmM0bN+N61pWkW0l4DPFg5LSRBAQEaGkmBU8nA64sy9tlWR5ibV20v5GF/LNgwQI+mfMJHsM9SLqbhPkucw7sOoCj4wt3t+lpsGkQnJgJdYdC16WvrFt76EYEbf46yqXQWH7vWoP/dffC3Lh4vA1ibKDPL12q812Hqhy++YSOM49xOyIh58avCLrm5ubs3LoT87/NSQlNwX2EO33H9yUwsHgUMtfJgCsIb2Pu3LmMXzQej2EeJN1OwupvKw7uPkiJEi+kaKUlKW+OXd6opHy1/uWlaV+yLDPn0B38/ANxtDRm28eN6FI7hzXgIk6SJPo18GTloHrEJqvoNOsYJ+5oHskDvDLo2tra8vf2v1GvV5Mel47zcGfa92tPcHDRP7JHBFyhSJk1axafL/sc96HuJN1KwvagLQd2H8De/oXdUymxShrT7f1Kge1Gn7y00leKKoNP1p7nl93X+egdFzaPaFgg23J1Sb0y9mwb1QgnaxN8l5x++Ztp2YNu0PLnl1xcXNi5fifRC6NBHyx9LWnduTVRUS8J4EWECLhCkTF37ly+XPkl7oPdSbyeiMMRB/bv2o+d3QtrsolRsLQdhJwGn8VQq99L+3scm0K3eSfYdiGM8S0rMqNnzUKbW5vX3GxMWT+0AVVdrRixMohVpx7k3PDfoFv2A9g+Gq5ufX6pUqVKbFq4iUezH2HoYIiqrYp2ndqRnJyspVlonwi4BSgmJgYfHx8qVapE5cqVOXHiREEPqdDasmULn/t//jzYOh13Yu/OvdjY2GQ2iguDJa3hyQ3osTrHBP1/nXsQTbuZ/3AnIoH5fesw8v1yRTZR/03ZmhuxclA93qtQgi82X+Kv/bfIcSOVgRF0X64cPbRxkFJf95lGjRqx4JsFhM4PxaysGbQBtazW4iy0SwTcAjRmzBhatWrF9evXuXDhgthN9jYqgsdQD5JuJlHinxLsCdhDljddn96Fxa2UoNtnI1T48KVdbT0fSvf5JzE11GfzyIa0qOKkhQkUTmZGBszvV4fOtdz4396bTNl2hYyc0saMzKH3OrAvD2v6wMPMN8m6du3K1z2+ptLjSsS7xDPz8kwtzkC7dDLgFoa0sGbNmuHl5YWXlxcmJiasX7/+tZ4fFxfHkSNHntdgMDIyyno3JuTa/gf7+TbwW6qXqE7LxJb8vePvrME28jYsaQOpceC7DTwb5diPLMvMPXyHMWvOU9PDhq0jG758d5XwnKG+HtN8ajCkSRmWnbjPuHXnSc/I4S7V1FbZHGHhqJyVFn71+aVx48axYdIG+lTuw4prK1hxtWjm5+pkTossy9uB7XXq1Bn8qna/nP6F60+v5+lrV7KrxMS6E/+z3YEDyq9F/xYg79w589fT3JRnDA4OpkSJEvj5+XHhwgVq167N9OnTnxfEEXLn8MPDfHb4M6o4VGFe83lYtM32htaTm8qarTod+u946ekMGWqZ7wOu4n/8Hm2ru/B7txq5K8otAKCnJ/FFm8pYmxry254bAPzetQYG2TeDWDoruc6LWipFbwbsBrvSzy+P9x7Po8RH/HbmN0palaQMZXB3LzoZIToZcAuLtylAnp6eTlBQEDNmzKBevXqMGTOGqVOn8v333+fnkIuMzZs3cybqDLuMd1HRtiJzm8/FwihbsI24rgRbZOgfAI45L9mkqDL4ZM15dl95zODGpZnUuvLrnfklPDfy/XIArw66tp5K0F3SWgm6A/eChZKypyfp8VOjn/Dd5cvov0cT8msIRzYdoUqVvDtyviAV6oCbmzvR/LJ+/XpWrlzJ1q1bMTTMei5Ubu5w3d3dcXd3p1495aBjHx8fpk7N/dE+xdnu3bvxm+KH22g3HFWOzGsxD0ujbL/6h19Vgq2ePvgGQImKOfYVnZjG4GVnOPsgmsltqzCwUekc2wm5l6ug61gZeq1XvkaruytfIyMzAMwMzTDdbUpKhRTsB9vTtntbTh04lTWPupAq1AG3oAQEBDB79mwCAgI0z70id3e4zs7OeHh4cOPGDSpWrMj+/fuLzE/x/HTo0CF6ftoT19GupIWnEbY+DKOe2Y5ueXwJlnUAfSPw3f7SExrCYpLpu+gUD58mM7NnLT6q7qKFGRQPuQq6Ht7QZSGs7aOcFdd9xfPNJ+OGjKNZz2a4jHUhvWM6nbp2Yv+e/RgbG2t1HnlNJ98003W+vr6EhITQsGFDvLy8WLRo0Rv1M2PGDHr37k316tU5f/48X3zxRR6PtGg5efIknYd3xuVjF9Jj0klfnc6+7fswNX1hO+7jS8pdk4GJsmb7kmB7LzKRrnNPEBGXyrKBdUWwzQcj3y/H+JYV2Xo+jE/XX8g5e6FyW2jzG9zYCTvHw7O0Mi8vL/x/8SdkQQhm5cy4X+k+g4cMzjntrBARd7hvIK92w3h5eZG9wLqQs3PnztG2f1ucPnYiIymD5KXJHNp+CA8Pj8xG4VdgaXswNFPWbO3K5NjX9cdx9F10mgy1zOoh9anmJmp25JcX73RNDPSZ2uUdzXzmuoMh5gEc/wtsPKCRcq5c+/btmXxzMr+s/wXnrs7s3rSbqVOnMmnSJG1PI8+IO1xB5129epXWvVrjMNwBOV0mbmEcf2/8m7JlMw8uVN4gaw8GxsoywkuC7fmHMXSfdxI9CdYNFcFWG0a+X45Rzcqx9sxDft51Pee71ObfQrUusO8buJiZYvnpp5/S3rk90f9E49TZiambprJ9+3btDT6P6WTALQx5uIJ2BAcH82GXD7EebI2kJxE1N4rda3Zn3SQSeeuFN8i2ZzlB9kXH70TSe8FJrE0N2TDsXco5ihxbbRnXogJ965di/pFg5hy+o9lATw86zoFSjWDLcLirvA8iSRJzZs+h7J2yJN1Kwn2wO/0n9Ofq1auafRQCOhlwRXlGASAsLIzmHZpj1t8MfVN9ImZHELA0gBo1amQ2irqjBFtZDf22vXTN9sD1cPovCcTN1pQNwxrgYWempVkIoATOb9tXpYOXK7/uvsHKU/c1GxkYQ48Vym8na/soX1uUTUEb1m5AvUVNRmIGDoMc6NCzA9HR0VqexdvTyYD7Xwr7wvnrKE5zze7wqcPod9PH0NaQRzMfsWnupudpdABE31OWEdJTlTtbx0o59rPnymOGLj9LJWdL1g5pgKOVZmaJkP/09CSmda1Bs0qOfLXlMtsvhGk2MrWFXmtB0lMO8kxWgmqJEiXYsnIL4fPDMbAyQGovcf1W3m560oZCF3BNTEyIiooqFoFIlmWioqJyTD0r6pJUSeww3oF5SXMiFkSw4tcVvPfee5kNYkPAvx2oEpXtuk45p9TtuvSIkSuDqOZmzYpB9bA1N8qxnaAdhvp6zOpVC+9Sdoxde56DNyI0G9mVhh6rlDfS1vWDDOV0iRo1arDwh4Xo7dXDuJwx+9X7tTz6t1fojklXqVSEhISQkpJSQKPSLhMTE9zd3TU2VxRlaRlpjDowipOPTjLtvWnUNK+Jg4NDZoP4cGWXUuITJdi61syxn4CLYYxZcx4vDxv8/byxNCk+f4e6Li5FRc/5Jwl+ksjaofWp7m6j2ej8atgyDGr3h7Z/Pq9ZnJ6ezl/n/2LJlSVMrj+ZbhW7aXPo/6lIHZNuaGhI6dJiN1BRpFKpiEuI48cLP3I87DjfvfsdLUq1yNoo6alysm78I+i75aXBduv5UMauPU+dUnYs9vPGopgchVNYWJkYssTPm86zjzPAP5BNwxtS0j7burpXT4i6BUd/B4cK0GAkAAYGBoypNYZbMbf4+dTPlLEuQ22n2oWifGahW1IQiia1Wo3fAD+a/tSUvff3Mr7OeDqV75S1UUqssvc+6g70XA0l6+XY16agEMauPU/d0nb4DxDBVlc5Wprg71cXVYZM/yWniU5M02z0/ldQub1y0OeN3c8f1tfT55cmv+Bu6c7ofaPx/sCb8PBwLY7+zYiAKxQ4WZYZO24s+1X7oRqkH03nI+ePsjZKS4SV3SD8slLMukzTHPvaFBTCp+sv0KCsPUv618XMSARbXVbO0YKFvnUIiUlm0LIzpKgysjbQ04NO88ClhrL994WSjlZGVjSNbUpMQgzRjaPx6eGDSvWS04R1hAi4QoH76aefWHVrFQ6tHYjaF0UTgyZZj8VRpcDqnsqxOF0WQoWWOfaz9Xwon62/QIMy9izs5y2OwykkvD3t+LO7F0EPovlkzXnNLcBGZtBzDRhZwOoeyrLSM9VcqhG6IBSzMmbcLXOXcZ+O0/LoX48IuEKBmjt3Lr/v/R3nbs7EnIihYWpDZs+anbkel6GCDX5w9zB0mA1VO+XYz46Ljxi37gLennYs8hXBtrBp844LX7apzO4rj/lhRw6bGqxcoMdKZe1+vS9kpAPQsmVLJnWdRMS2CGyb2LLy8kqWLVum5dHnnk4GXLHTrHhYt24dExdNxNXXlfiL8VS6V4kVy1dk1hZWq2HLCKWwSZtpypsoOdhz5TFj1pyjVkkbFvcXwbawGtS4DAMalmbJsXss/ueuZgP3OtBuOtw9An9/+fzhiRMn0livMfEX4nHp7cKYX8dw9uxZLY4893Qy4IqdZkXf33//zZAfhijHmd9JwvGkI5s3bM4svyfLsGs8XFoHzSYrBU5ysP9aOB+vCuIdd2uW+NXFXLxBVqh99VFlWlZ14ocdVzlwPYc3wbx6Qf2RcGru82PXJUnCf4k/FoctUEWqcBrsRBffLjx58kTLo/9vOhlwhaLt1KlTdB/THdcRrqQ9TsNkhwk7t+7EwuKFExsO/ACBC+Hd0dD40xz7OXQjguErgqjsYsXSAXVFNkIRoKcn8Ud3Lyq7WDFq1TmuP47TbNTiOyjzPgSMhQenALCwsGDzms1E+0ejZ6KHfkd9uvfqTnp6upZn8Goi4ApadfPmTdr5tsN5hDMZ8Rmo1qj4e/vfWTc2HJsOR6dBLV/lH1cO+ZXHb0cyZPlZyjtZsHxAPazEpoYiw8zIgEW+3liYGDDQ/wxP4lOzNtA3gK5LlFKOa/souw6B8uXL4z/Nn7BFYZiVN+O6y3W++uqrApjBy4mAK2iVnpUeLqOUYt+xi2LZs3FP1kMCz/rD3q+VN8fa/pFjsD17/ymDlp3B096M5QPrYW0mgm1R42xtwsJ+3kQlpjJkeQ7pYqa20GM1qJKVoKtSdp5+9NFHjGs7jic7n2D/gT3zjs7j9OnTBTCDnImAK2hNbGos40+Nx8TGhKrXq7J9+XYqVnzhrLErW2D7J1CuBXSa//y4lRddDo2l/+JAHC2NWTGoHnaiNkKR9Y67NX929+LcgxgmbLioWT/FsRJ0ng9h52DHuOenRXz11VfUSqxF4vVESg0uhWVZ3SnDKQKuoBXJ6cmMOjCK+3H3md5sOhtmbcDb2zuzwZ0DsHEQeNSDbsvAQDOQ3gyPp++iU1iZGrJycH0cLYtfUZ/iplU1Fya0qsi2C2H8tf+2ZoNKbeC9iXB+pbLmD+jp6bFy+UrmfjQXezN7xh0aR2yqbmQ8iYAr5CuVSsWOXTsYf3g85yPO83Pjn6nvUj9ro5CzsKaPcrJur7XPT2990d3IRHovPIWhvh6rBtfDzcZUo41QNA1/ryxdarnzx76b7Lr0SLPBe59DhVaw+3O4fwIAW1tbPmz0Ib83/Z3HSY/58p8vUctqLY9ckwi4Qr5Rq9UMHDSQoRuGcjjkMJPqTqKlZ7ZdYk9uwMouYFEC+mwEUxuNfkJjkum94CQZapmVg+pRyt5cOxMQdIIkSfzUuRq1Stowbt0FroZly1z4d/uvTSmlnGNcZp1dL0cvJnhP4HDIYWYHzWbv3r1aHn1WIuAK+WbSpEn8nfQ3tk1sidgSgfXdbHnVMQ9gWUflOPO+W8DSWaOPJ/Gp9Fl4ivjUdJYNqEt5J91ZjxO0x9hAn7l9a2NtasjgZWeISsiWuWBqo9TQVSXB2r5KUfpnelTsQVPHpsy7NA+fz3w4evSodgf/AhFwhXzx559/sjBoISXaluDpwae0c2hH69atMxskRiqVv9ISoc8mpeh0NrFJKvotPs3j2BT8/bzFgY/FnKOlCfP71SYyIZURK4NQZWRbInCspJyLFnpGOXL9GUmSuDX7FskPk3Ed4kqPIT14/PixlkevEAFXyHOrV69myuopuPRyIe5sHLWjazN3ztzM+gip8bDSB2JDlTVb52oafSSlpePnf5rbEfHM71eb2qXsNNoIxU91dxt+9anOqbtP+Xb7Fc0GVdorG2WCliophs8smL2AhFUJSHoSJt1M6NG7R4FsitDJgCtqKRRee/fuZfjPw3Ef4q6csnrBnbVr1mJg8GwXWHqqkjf56CJ09YdSDTT6SE3PYOjys5x/GMNfPWrSuHwJ7U5C0GkdvNwY+l4ZVpx8kPNhlO9/CWWbKXe5oUpNBQ8PD5b/tZzQxUplsZtuN5k8ebKWR66jAVfUUiiczp49S4/RPZQtu4/SMN1tyvbN2zE1fZZRoFbD5mEQfAg6zISKrTT6SM9Q88ma8xy9FcnULtVp/Y6LdichFAoTWlbi/YolmLL1CqeCo7Je1NOHLovAwgnW+UKicr1FixZ82u5TIndFYv+BPXMOzWHbtm1aHbdOBlyh8Llz5w5t+7SlxPASZCRmkLoqld1bd2fWtZVl2D0RrmxStut69dLoQ62WmbTpErsuP2Zy2yp0q+Oh5VkIhYW+nsT0njUpaW/GiJVBhMUkZ21gZqfkcyeEw6ZBoFZ2qn311VfUiKtB4o1E3PzcGDhxIMHBwVobtwi4wltTqVR81PUjLPtbIulLRC+IZveG3Xh4vBAwj0yD0/OhwcfQcIxGH7Is89POa6w/G8LoD8ozsJE4t054NSsTQ+b3rUNquprhK85qbv91q6WU9bxzAA79DGRuilBvU5ORkoGtry0+vXy0diitCLjCW1OhovSnpTG0M+TxrMdsXrSZKlVeOLb8zGI4+ANU7wEtvs+xjzmH77Dwn7v4NijF2ObltTRyobAr52jB791qcCEklslbLmtu/63tCzX7wpHf4MYuAOzt7Vm3eB2P5j/CyMmIyFqRjBo9SivjFQFXeCsqtYpxh8cRlhHGmDJjWP7rcho2bJjZ4Np22PEplP9QWbfV0/yWW336Ab/uvkH7Gq5MaVe1UJy+KuiOllWdGd2sHOvPhrDi1APNBm2mKWeibRqqHEAK1K1bl59H/Ez4xnCs61kTJAeRlJSU72MVAVd4Y7Is883xbzgWeoyv63/NkA+G0LZt28wG947BhoHgVlvJSNDXrOq189Ijvtx8iaYVSzCtaw309ESwFV7fJ80r8H7FEny77QqB955mvWhoAt2WKz/s1/aFNCWwDh8+nBbWLXCIdUDdRM3txBxqNeQxSeMWXIfUqVNHPnPmTEEPQ8jBzJkzia0ey5q7axjhNYLhNYZnbRB+BRa3BksnGLBHeRMjm2O3I/FbEsg77tasGFhPHI0jvJXYZBUdZv5DQmoGAaMa4WydrbjRrX1K/neNntBxNkgSarWaeFU83QO6o5bVrG+3Hmvjt8uOkiTprCzLdXK6Ju5whdc2Y8YMpmyZwpq7a2jh2IJh1YdlbRDzAFZ0UYrQ9NmUY7C98DCGIcvOUNrBnMXi0EchD1ibGjKvbx2S0tIZvvIsqenZ3kQr3xzemwAXVkGQctCknp4e1sbWTHtvGpHJkXzxzxf5WuRGBFzhtWzYsIGvln2FS28X4oLiiFwTmXXNNekpLO+s/NrWZ6NSlT+bO08S6L/kNHYWRiwbWFcUEBfyTEVnS6Z1rcG5BzF8H5DD6b/vTczcFBF2/vnD1RyqMd57PEdCjtB/Vn/mzZuXL+MTAVfItSNHjjDom0HPD350DXJl7py5mQ3SEmFVN+UOt+dqcKqq0cfj2BT6LTqNvp7E8gH1cLISNW2FvNXmHReGNlF2om04G5L1op4+dF4I5g5KZbHk6OeXWju3xvyhOUHmQUycOTFfTv4VAVfIlStXrtBlSBdlF9mTNIx2GBGwJQAzs2e1azPSYb2fspXSZxF4NtToIyYpjX6LTxGbrMLfry6eDqLMopA/xresSP0ydny5+RJXwrKVCDC3h65LlTKOm4crOyABExMTUralkBaehvNAZ7r270p0dHQOvb85EXCF/xQSEkLrbq2xG2KHOk1N0rIkdm/ejb29vdJAliFgDNzao6TgVG6n0UdyWgYDl57hXmQS8/vWFpW/hHxloK/HjJ61sDEzZPiKIGKTVFkbeHhDyx/h5i449iegBNwNqzbwdOlT9M310Wunx8VLF/N0XCLgCq8UExND646tMe5pjL6ZPk/mPWHH6h14enpmNjr4I5xboayPeQ/U6EOVoebjVUEEPYjmzx5evFvOQaONIOS1EpbGzO5dm0exyYxddx61OltGVt0hUK0LHPge7h4BoFy5csz/cT76B/UpUbUEHjXydnu5CLjCS6WmptKhcwcSmidg5GJE6KxQ1s5Yi5eXV2ajwIXKLp5a/aDpJI0+ZFmpj7D/egTfdahGG1GMRtCi2qVsmdy2CgeuRzDzYLY8W0mCdn+BfTklXzxeqZHbuXNngpYFsavLLspYl8nT8ehkwBXlGXXDtN+nEVwpGIvKFoQuCmXWxFk0b948s8G17bDjM+U8qY9yPtL8l9032HA2hDEflKdv/VJaHL0gKPrWL0Wnmm78se8mh25EZL1obKEUuUmNV4JuhlIj18DAABsTmzwfi04GXFGeUTfoN9PHpp4Nj9c+5ouOX9CnT5/Mi/dPKN+g7nXAZwnoG2g8f/E/d5l7+A696pXkE1EfQSggkiTxU6d3qOhkyZg15wmJzraF17EytP0D7v8Dh37K17HoZMAVCt6yK8tYeWMlPSv2ZP7g+Ywfn3lkCRHXYHV3sCkJPXM+ZXfbhTC+C7hKy6pOfN+hmqiPIBQoUyN95vSpjVotM2JlkOamCK+eyrLY0d/h5t/5Ng4RcAUNu+/t5rczv9G8ZHMm1p1Ip06dMgNmbKiyi8zARNnYYG6v8fxjtyP5dN156nraMb1HTfRFfQRBB5R2MOe3rtW5GBLLDwHXNBu0/hWc3oHNQyDmYb6MQQRc4bnDhw8z9n9j+eLoF9RyrMXPjX9GX++FLbfJMcpe9JQ46L0BbDXXZC+HxjJ0+VnKOFiwwLcOJoZiy66gO1pVc2Fw49IsP3mfredDs140NIVuS5/llPeH9LQ8f30RcAUALl++TJehXdhjsgeDRAP+eO8PTAxe2AX271lkkbegxwpwqa7Rx4OoJPovCcTa1JClA+pibSq27Aq6Z0KrSnh72vL5xkvcDI/PetG+rFJGNPQM7P06z19bBFyB0NBQ2nRrg8MQB9RpakKnh5IS80IF/H/PIrt3VKmyVKapRh+RCan0W3yKDLWapQPqalZqEgQdYaivx8xetTA3NmDYirMkpGY7vbdqR6g3DE7NVd6vyEMi4BZzsbGxtOrQCuNexuiZ6fFkrrKxwdXVNbPR3smZZ5FV76bRR2JqOgP8A3kcl8JCX2/KOVpocQaC8PqcrEyY0bMm9yITmbjxouZJES2+h/4BSgZDHhIBtxhLS0ujk08nEj5IwNjFmNA5oayZsSbrxobjM+HETOUn/rujNfpQZagZsTKIK2FxzOpVi9qlbLU3AUF4Cw3K2vNZy4rsuPiIpcfvZb1oYASejfL8NUXALabUajV+A/y4VeYWFlUsCF0cyoxPZ9CiRYvMRpc3wt9fQpUO0PInjY0NsiwzceNFDt98wo8dq/FBZSctz0IQ3s6wJmVpXtmRH3de49yDvC1UkxMRcIupL7/8kv1p+7F514bwDeGM/2g8/fr1y2xw96iyblvyXeg0Xylrl81ve26wKSiUsc0r0KNuSS2OXhDyhp6exO9dvXCyMuHjVeeITsz7zIQsr5evvQs6afbs2cwPnE+Jj0oQdSCKTm6dmDTphToI4VdgTW+wLQ09VipnQmWz9Pg9Zh9SdpGN/qCcFkcvCHnL2syQWb1q8SQ+lXE5FbnJQyLgFjNRUVF8t+o75cSGc3HUjqnNrJmzsm1s8Hl2PM7GHI/H2XXpEd9sv0KLKmIXmVA01PCw4au2lTl44wlzj9zJt9cRAbeYeZjxENchrqQ/TMcp0Ik1q9ZgYPCsDsK/GxtS46H3+hyPxzkVHMWYteepVdKWGT3FLjKh6OhbvxRtq7swbc8NTgZH5ctriIBbjNyLvceoA6NwsXBhdbfVBGwOwNz82akL2Tc2OL+j8fyb4fEMXnYGD1tTFvYTu8iEokWSJKZ2qY6nvTmjVp/jSXxqnr+GCLjFRFRyFMP3DUdCYk7zOdSsVBNHR0floloNW4a/cmPDo9hkfBefxsRQn6UD6mJrbqTdCQiCFlgYGzC7Ty3iU1SMWXOOjDxezxUBt4hLTEykc/fODNwxkMjkSGZ+MJOSVtkyCvZ9raSANf8mx40Nsckq+i8OJD4lnSV+3rjbalYHE4SiopKzFd93qMa1R3HcjUzM0741i5gKRUZ6ejrde3bnUtlLWMZbMsx5GNVLZKuBcHIuHJ8B3oOh4ScafaSmZzBk2RmCIxPw96tLVVdRo1go+rrW8aBFFSdszPL2Nzlxh1tEybLMiJEjOGd3DisvK8KWh2FwP9vP16tbYffnUKkttP5FY2ODWi0zbt0FTt19yrSuNWgoziITipG8DrYgAm6R9fPPP7MlbAt279vxZPsThtUfxrBhwzIbPDgJGweDuzd0WZjjxoYfd15jx8VHTGpdiQ5eblocvSAUTSLgFkHLly/ntx2/4eTjRMzxGD4w+YAff/wxs8GTm7Cqu5L21XONUgc0m4VHg1n0z136v+vJkCZ5e5CeIBRXYg23iNm3bx8f//oxHp94kHA1gXJ3yrFk55LMzQnx4cqJDfqGShHxHE5sCLgYxg87rtG6mjOT21YRGxsEIY/o5B2uOLX3zVy8eJHuH3fHbYQbqY9SsdhvwZaNWzAyerYWlRoPq7pCUhT0Wgd2pTX6OBkcxbi1F/D2tOWP7l5iY4Mg5CGdDLji1N7X9/DhQz7q8RElhpZAnaImZXUKu7bs4vnfYYZKOTbk8WXo6g9utTT6uBkez5BlZyhpb8YCsbFBEPKcWFIoImYumIlJLxP0TPSI+CuCg+sO4u7urlyUZQj4BG7vg3Z/QYUPNZ7/ODbl+cYGfz/vfHmHVhCKO528wxVejypDxZN3n2DqbkrYnDDWzlzLO++8sDX38C9wbgU0mQC1fTWeH5eiov+S02JjgyDkM3GHW8jJsszXx78mMDyQnxr/RMV3K1KhQoXMBudWwKGfoUYveP8LjeenpasZtvwstyMSWOLnLTY2CEI+EgG3kJtxbgYBwQGMrjmadmXbZb14ex9sGw1l3of2f+W4sWHChgscvxPF/7rVoHH5ElocuSAUPyLgFlLz5s1jV+gu7pS9g08FHwa9Myhrg0cXYJ0vOFWBbsuUNLBsft1zgy3nwxjfsiKda7lraeSCUHyJNdxCKCAggAnzJnC79G1MH5nySbVPsubKxjyAlV3BxAZ6rQcTK40+lp24x9zDd+hdryQjmpbV3uAFoRgTAbeQCQwMpN+EfrgPcyflfgoZWzPQl15I30qOVk5sSE+BPhvAykWjjz1XHjNl2xWaV3bk2/ZVxcYGQdASsaRQiAQHB9O+X3uchjuRHpeOeqOaHXt3YGFhoTRQpShnkUXfhb6bwbGyRh9n70czevU5arjbMKNnLQz0xc9cQdAWEXALiaioKFp3bo2lryXoQ8yiGI5sOYKzs7PSQK2GLcPg/jHosgg8G2n0EfwkgUFLA3GxNmGRbx1MjcTGBkHQJnF7UwgkJyfTvnN7VG1VGDoY8mj2IzYt3ETFihUzG+37Gq5shhbfwTs+Gn08iU/Fd8lp9CSJpQPqYm9hrMUZCIIAIuDqPLVaTd9+fXn4zkPMypoROj+Uhd8upGHDhpmNTs1TiojXHQLvjtboIyktnYFLA3kSn8qi/t6UsjfX4gwEQfiXCLg67rPPPuOYwTGsva15vPYxU3pOwcfnhTvYa9th10SliHirqRq5tukZaj5edY7LobHM6FkLLw8b7U5AEITnRMDVYREREWy6twmHVg5E/h1J70q9GTt2bGaDh6dh4yBwrwOdF2gUEZdlmclbL3PgegTfdahGiypOWp6BIAgvEgFXh11IvoBVByvUN9Q0SGnA/37/X+bFyNtKEXErV6WIuJFm/YNZB2+z+vRDRjQtS5/6pbQ4ckEQciKyFHTU+YjzTDo6ieolqvNbm9+wsbBBX//ZHWzCE1jZRVk+6L0BzDXPGtt4NoRpf9+kU003xresqHFdEATtEwFXx8iyzP24+4w6MAonMydmNJuBrYltZoO0RFjdXTm5oX8A2GvuEvvnViQTN16kYTl7fulSXWxsEAQdIQKuDomIiKCNTxssh1oi6UvMaT4na7BVZyhrtmHnoPsKZe02m6thcQxbcZZyjhbM6VMbIwOxaiQIukL8a9QRiYmJtO3Uluim0UQkRdCJTpS0KpnZQJZh1wS4sRNa/wqVPtLoIywmGT//01gYG7DEzxsrE82CNYIgFBwRcHVARkYGPXv1JLxOOKalTQmZG0JZs2xLBcemQ+BCJc+27mCNPmKTVfgtCSQpNQP/Ad64WGuexCsIQsESAbeAybLMqNGjOGN9BqtaVjxa+YifBvxE+/btMxtd2gD7pkC1LtD8W40+/i0iHhyZwNy+tankrFkdTBCEgicCbgGbNm0a6++ux/4De57sfMLAWgMZOXJkZoN7/8CW4VCqIXScA3pZv2SyrBQRPxEcxS9dqtOwnGbGgiAIukEE3AK0Zs0aftjwA87dnYk5GcN70ntMnTo1s0HEdVjTC2xLQ4+VYKBZ/+A3UURcEAoNkaVQQI4cOcKwH4bh/ok7iTcS8bzmydLdS9H79w42/jGs9AEDE+i9HkxtNfpYcfI+sw/doWddUURcEAoDEXALwLVr1+gypAuuo11Je5KG8S5jtu7fiomJidIgNV45sSHpKfjtBFvNXWL7robz9dbLNKvkyPcdRBFxQSgMRMAtAOt2rsNusB3qdDWJyxI5vus4dnZ2ysUMFazvD+FXlC27rl4az7/wMIZRq89R1dWaGT1riiLiglBIiH+pWpakSuJKhSuY2ZsROT+S7Su24+npqVyUZQgYq5y22/YPqPChxvMfRCUxcGkg9hZGLOpfB3Nj8TNTEAoL8a9Vi9LV6Xx6+FNuRN9gRvMZVG1XFXt7+8wGR6bBueXQZDzU9tV4/tPENPovOY0qQ2bNkLo4WppocfSCILwtEXC1QJZlMjIy+OH0D/wT+g9fN/iaJu5NsjY6vxoO/gDVe8D7X2r0kaLKYPCyM4TEJLNyUD3KOVpoafSCIOQVEXC1YOrUqex6uovoatEMemcQXSt0zdrgzkHY9jGUfg/az9AoIp6hlhm79jxBD6KZ2bMW3p52Why9IAh5Razh5rOVK1fyy/ZfiK4Wjd4NPXqX7J21QfgVWNsXHCpC9+VgYKTRx487rrHr8mO+bFOZj6prHnsuCELhIAJuPjp48CAjfxmJ20A3Eq4lYHfaDiurF7bdxobCCh8wtlRybU2sNfpY9M9dFh+7S/93PRnYqLQWRy8IQl7TWsCVJKmjJEkLJEnaKkmS5tvvRcyVK1foNqIbrsNdSQtPw3yvOZs3bMbI6NkdbEqskmublqAEW2s3jT52XXrEDzuu0rKqE5PbVhG5toJQyOUq4EqStFiSpAhJki5ne7yVJEk3JEm6LUnS56/qQ5blLbIsDwb6A93feMSFQFhYGG26tVFybVPVpKxIYefmndjY2CgN0tOUZYTIG9BtGThX0+jjzL2njFl7npoeNkzvURN9PRFsBaGwy+2bZv7ATGDZvw9IkqQPzAJaACFAoCRJ2wB94Odszx8gy3LEs4+/eva8Iik+Pp42Hdtg2N0QfTN9wqeHs3/1fkqWfFbbVpZh+2i4e1gpRlP2fY0+gp8kMHjZGVytTVjo642Job5GG0EQCp9cBVxZlo9IkuSZ7eG6wG1ZloMBJElaA3SQZflnoG32PiTl9+GpwC5ZloPeatQ6SqVS4dPNh5gmMVi4W/Dwr4es+3MdXl5emY0O/gQXViupX169NPqITEil/5JAJEnC368uduaab6IJglA4vc0arhvw8IXPQ5499jKjgOaAjyRJw17WSJKkIZIknZEk6cyTJ0/eYnjaJcsyw4YP44rrFSzfsSRsaRh/jP6Dli1bZjY6uxSO/Ao1+yqbG7JJSktnoH8gEfEpLPKtg6eDuRZnIAhCfnubPNycFhXllzWWZfkv4K//6lSW5fnAfIA6deq8tD9dEx0dzQn5BLZNbInYGsGopqMYMGBAZoNb+5Rtu2U/ULbt5pBrO3r1eS6FxjK3T21qltSsDiYIQuH2NgE3BPB44XN3IOzthlN4HX16FP0m+hjeNKSlZUu++eabzIuPLsB6X3CqAt2Wgn7Ws8ZkWeabbVfYdy2c7zpU5cOqztodvCAIWvE2ATcQKC9JUmkgFOgBaC5KFgMnwk7wzfFvqOdSjz98/sDY0DgzhSvmgZL+ZWIDvdYrObfZzD8SzPKT9xnapAz9GnhqdeyCIGhPbtPCVgMngIqSJIVIkjRQluV04GNgD3ANWCfL8pX8G6ruSUtL48bTG4w9NBZPa0/+aPoHluaWmbm2yTFKsFWlQJ8NYKW5S2z7hTB+3nWddjVcmdiqknYnIAiCVuU2S6HnSx7fCezM0xEBkiS1A9qVK1cur7vOMyEhITRt1xTHTxwxNzNnTvM5WBq9cPeangpr+0DUHei7CRwra/RxKjiKT9ddoK6nHdO6VkdP5NoKQpGmk1t7ZVneLsvyEGtrza2uuiAuLo42ndogdZGISY6hUXgjnM1fWHeVZdg6Eu4dhY6zoXQTjT5uR8QzeNkZPOxMmd+vNsYGItdWEIo6nQy4ukylUtGlaxfim8Zj7GJM6JxQGlVslLXR/u/g0npoNhmqd9PoIyIuBd/FgRgZ6OPvVxcbM5FrKwjFgQi4r0GWZQYPGcw1j2tYVLMg1D+U6WOn06JFi8xGZxbDP/+DWr7Q+FONPhJT0xmwNJDopDSW9PfGw85MizMQBKEgiYD7Gr777jt2Ru/EtrEt4ZvDGdtiLP37989scPNv2PEplGsBH/1PI9c2PUPNyFVBXA2LY1avWrzjrptLJoIg5A8RcHPJ39+fP/f9iVMnJ6KPRPOR3UdMnjw5s0HYOeXwR+d3oKs/6Gd9P1KWZSZvvcyhG0/4oeM7vF/JUavjFwSh4IkTH3Jh7969jPljDB6jPUi4nEDlkMrMD5ifmWsbfR9WdgMze+i1Dow1j7+ZfegOq08/ZOT7ZelVr6SWZyAIgi7QyTtcSZLaSZI0PzY2tqCHwoULF+gxugduw91ICUvB6rAVG9dvxNDw2W6x5GhY6QMZqUpdW0vNXWKbz4Xw254bdKrpxmcfVtTyDARB0BU6GXB1KS3scNBhSgwpgTpJTdrqNHZu3pl5akN6KqzpDdH3oMcqcNTcuHD8diQTNlykQRl7fulSXRQRF4RiTCcDrq6IT4vnoN1BLOwsSFiWwI61O3Bze1YQTa2GLcPh/jGlrq1nI43n33gcz9DlZyntYM7cvrUxMhB/3YJQnIk13JdQZagYe2gs92LvMefDOdToVgNTU9PMBvu/hcsbofk38I6PxvMfx6bQf8lpTI30WeJXF2tTQ402giAULyLgZiPLMsnJyfwY9COnHp3ix0Y/Ut+lftZGgQvh2J9QZyA0/ESjj/gUFX7+gcQlq1g3rAFuNqYabQRBKH5EwM3mm2++YVvUNtK90xnpNZL2ZdtnbXBjN+wcDxVaQetfNXJtVRlqRqwM4mZ4PIv7e1PVteDXoQVB0A1iUfEFixcvZsahGaR7p6MKUvGR3UdZG4QGwQY/cK4OPotzzLX9cvMljt6K5KdO1XivQgktjl4QBF2nkwG3INLC9uzZw9i/xuLq60r8pXhK3SyFs/MLKV7R92BVNzB3UHJtjTSPv5lx4DbrzoQwulk5unuLXFtBELLSyYCr7bSw8+fP0+uTXrgNcyMlJAWbIzZsWLchM9c26alS1zZDBb03gqWTRh8bzobwv7036VzLjbEtKmhl3IIgFC7Ffg334cOHtOvVjhJDS5CRlIFqrYqde3diafmstq0qBdb0Uu5w+22FEprB9J9bkXy+8SINy9kztbPItRUEIWfFOuDGxsbSplMbTHqboGesx5O/nnBowyFcXV2VBmo1bBkGD04oa7al3tXo49qjOIatOEs5Rwvm9BG5toIgvFyxjQ5paWl07tqZhA8SMHI2ImxOGGtnraVq1aqZjfZNgSubocX3UK2LRh+PYpPxWxKIhbEBS/y8sTIRubaCILxcsQy4siwzaPAgbpa+iUUVC0IXhTJzwkyaNWuW2ej0Ajj+F3gPhndHafQRl6LCb0kgCanpLO7vjYu1yLUVBOHViuWSQmJiIldtr2JT04bwjeGM/2g8ffr0yWxwfQfsmgAVWkPrX3LOtV0RxO2IBJb4eVPF1UrLMxAEoTAqlgF3d9huUmumYv3AmnpO9fjiiy8yL4achQ0DwcULfBaBXtazxmRZZtKmS/xzO5LffKrTuLzItRUEIXd0MuDm56m9R0OO8sPJH2jo1pDpvaejj35mVsHTu0qurYXjS3Nt/9x3iw1nQ/ikeXm61vHI8/EJglB06eQabn7k4cbHx3M16iqfHv6U8rbl+f293zE2MMbA4NnPnKSnSl1bOQP6bAQLzTvXdWceMn3/LbrWdmfMB+XzbGyCIBQPOnmHm9cePHhAozaNcB7njI2VDbM+mIW54Qt3r6pkWN0DYh4qubYOmsH0yM0nfLHpEo3LO/BT53dErq0gCK9NJ+9w81JMTAytO7bGuJcxCakJVLhUAUezF84TU6th81B4eAo6z4NSDTT6uBoWx4iVQZR3smR271oY6hf5vzZBEPJBkY4cqampdOzSkaQPkzByMiJ0dig+72erXbt3MlzdCh/+AFU7afQRFpOMn/9pLE0MWNLfG0uRaysIwhsqsgFXlmUGDBzAnfJ3sKis5NrO/nw2TZs2zWx0ah6cmAl1h0KDjzX6+DfXNik1gyV+3jhbm2hvAoIgFDlFNuBOnjyZfSn7sGlgQ/iGcCa2m0jv3r0zG1zfAbsmQsWPoNXPGrm2aelqhq84y50nCcztW5tKziLXVhCEt1MkA+6CBQuY/c9sHNs58vTQUzq6duTzzz/PbBByRsm1dasFXRbmmGv7+caLHLsdxS9dqtOwnIOWZyAIQlFU5LIUdu3axWezP8NjlAfxF+Kp8aQGs7fMfiHXNljJtbV0gp5rwchMo48/9t5k07lQPm1RgS613bU8A0EQiqoiFXCDgoLoPa43bp+4kfIwBbvjdqzbty4z1zYxClb4gKxW6trmkGu7NvABfx24Tfc6HnzcLO83XgiCUHzp5JLCm574cDXkKk7DnchIyCBjfQY7t+zEwsJCuahKhjU9ITYEeq4BB81geuhGBF9svkyTCiX4oVM1kWsrCEKe0smA+6Y7zcrULoODvQNsgh3rdmQekaNWw6Yh8PA0dJ4PJetrPPdyaCwjVwZRUeTaCoKQT4rUksK7bu+yr8c+DHoYoK//whtheyfDtW3Q8ieo2lHjeaExyQzwD8Ta1JAlft5YGBepvxZBEHREkYssxvrGWR84OTcz17b+CI32sckq/JacJlmVwYZh7+JkJXJtBUHIH0X79+Zr22H351Cp7UtzbYctP8vdyETm9alNRWfLAhqoIAjFQZG7w33uYSBsHARutaHzghxzbSduvMiJ4Cj+6F6Dd0WurSAI+axo3uFG3YHV3cHSBXrlnGv7+9832fws17ZTTZFrKwhC/it6ATcx6lldW1mpa2uueee6+vQDZh68TQ9vkWsrCIL2FK0lhX/r2saFge92sC+r0eTgjQi+2nKZ9yqU4PuOItdWEATtKVoB9+o2CAmEbkvBo67G5RdzbWeJXFtBELSsaAXcGt3BqSo4V9O4FBKdhJ9/ILZmRiLXVhCEAlH0ok4OwVbJtQ0kRZXBykH1RK6tIAgFQid/p37TWgo5SU3PYOjyM9yLSmRe39pUcBK5toIgFAydDLh5dWqvLMtM2HCRk8FP+dWnOu+WFbm2giAUHJ0MuHll2t832Ho+jPEtK4pcW0EQClyRDbirTj1g1sE79KzrwYimmulhgiAI2lYkA+7B6xFM3nqZphVL8H0HkWsrCIJuKHIB93JoLCNXBVHZxZJZvWphIHJtBUHQEUUqGr2Ya7vY1xtzkWsrCIIOKVIBN/hJInoS+Pt54yhybQVB0DFF6hawSYUSHB7/PiaG+v/dWBAEQcuK1B0uIIKtIAg6q8gFXEEQBF0lAq4gCIKWiIArCIKgJSLgCoIgaIkIuIIgCFqikwE3L8szCoIg6AqdDLh5VZ5REARBl0iyLBf0GF5KkqQnQAzw4q2u9Quf5/SxAxD5li/9Yr9v2i6na9kfe9XnYm6vT8zt5e1y+7iY29vPrZQsyyVyvCLLsk7/Aea/7POcPgbO5PVrvkm7nK69ai5ibmJu+Tm33D4u5pa3c8v+RyeXFLLZ/orPX/ZxXr/mm7TL6dqr5pL9czG31yfm9vJ2uX1czC1v55aFTi8pvAlJks7IslynoMeRH8TcCicxt8IpP+ZWGO5wX9f8gh5APhJzK5zE3AqnPJ9bkbvDFQRB0FVF8Q5XEARBJ4mAKwiCoCUi4AqCIGhJkQ+4kiSVlCRpmyRJiyVJ+rygx5OXJEmqIknSOkmS5kiS5FPQ48kLkiSVkSRpkSRJG154zFySpKWSJC2QJKl3QY7vbbxkbhqPFUYvmVvHZ1+zrZIkfViQ43sbL5lbZUmS5kqStEGSpOG57iyvE3u18QdYDEQAl7M93gq4AdwGPn/2WHNg6LOPlxX02PN4bp8CjZ99vK2gx54Xc3rh2oYXPu4LtHv28dqCnk9ezu1VjxX0nzycmy2wqKDnk09z03uduRXWO1x/lL+Y5yRJ0gdmAa2BKkBPSZKqAOeAHpIkHQAOanmcb8Kf3M9tOcrcfgPstTzO1+FP7ueUE3fg4bOPM/JpjG/Kn7ebmy7zJ2/m9tWz5+gSf95ybpIktQf+Afbn9kULZcCVZfkI8DTbw3WB27IsB8uynAasAToAfsAUWZabAR9pd6Sv73XmJstyhCzLI4HPefv97PnmNb9eOQlBCbqgY9+zeTA3nfW2c5MUvwC7ZFkOyt/Rvp68+LrJsrxNluV3gVwvc+nUN+9bciPzLgiUf6RuwG5gtCRJc4F7BTCuvJDj3CRJ8pQkaT6wDPitQEb25l42J/tnX6uakiRNenZtE9BFkqQ55OO2yzyU67m9ZL667HW+bqNQlvR8JEkapuVxvonX+bo1lSTpL0mS5gE7c/sCRemYdCmHx2RZli8Dhf0NpZfN7R4wRMtjySsvm1MUMCzbg4kov6kUFq8zN43HdNzrzO0v4C+tjCpvvM7cDgGHXvcFitIdbgjg8cLn7kBYAY0lrxXFuRXFOf1LzK1wyve5FaWAGwiUlySptCRJRkAPYFsBjymvFMW5FcU5/UvMrXDK/7kVdHrGG6Z0rAYeASqUn0oDnz3eBrgJ3AG+LOhxirkV3TmJuYm5vckfUbxGEARBS4rSkoIgCIJOEwFXEARBS0TAFQRB0BIRcAVBELREBFxBEAQtEQFXEARBS0TAFQRB0BIRcAVBELREBFxBEAQt+T/ksPN/wt/JIgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pop = ares.populations.GalaxyPopulation(**pars)\n", "\n", "redshifts = [4,5,6]\n", "Mh = np.logspace(8, 13)\n", "\n", "pl.loglog(Mh, sfe_dpl, color='k', ls='--', lw=3)\n", "for z in redshifts:\n", " fstar = pop.SFE(z=z, Mh=Mh)\n", " pl.loglog(Mh, fstar, label=r'$z={}$'.format(z))\n", "\n", "pl.legend()" ] }, { "cell_type": "markdown", "id": "869a3ad7", "metadata": {}, "source": [ "**NOTE:** The only method of ParameterizedQuantity objects ever called is the ``__call__`` method, which accepts ``**kwargs``. As a result, we must always supply arguments accordingly (i.e., supplying positional arguments only will not suffice), hence the ``z=z, Mh=Mh`` usage above." ] }, { "cell_type": "markdown", "id": "7675f654", "metadata": {}, "source": [ "## Multiple `ParameterizedQuantity` objects at once" ] }, { "cell_type": "markdown", "id": "f8ea1b23", "metadata": {}, "source": [ "In general, we can use the same approach outlined above to parameterize other quantities as a function of halo mass and/or redshift. For example, we can use a double power-law SFE model and set the escape fraction to be a step function in halo mass, " ] }, { "cell_type": "code", "execution_count": 6, "id": "06f8c6ae", "metadata": {}, "outputs": [], "source": [ "pars = \\\n", "{\n", " 'pop_sfr_model': 'sfe-func',\n", " 'pop_sed': 'eldridge2009',\n", " \n", " 'pop_fstar': 'pq[0]',\n", " 'pq_func[0]': 'dpl',\n", " 'pq_func_par0[0]': 0.05,\n", " 'pq_func_par1[0]': 2.8e11,\n", " 'pq_func_par2[0]': 0.5,\n", " 'pq_func_par3[0]': -0.5,\n", " 'pq_func_par4[0]': 1e10,\n", "\n", " 'pop_fesc': 'pq[1]',\n", " 'pq_func[1]': 'step_abs',\n", " 'pq_func_par0[1]': 0.02,\n", " 'pq_func_par1[1]': 0.2,\n", " 'pq_func_par2[1]': 1e10,\n", "}" ] }, { "cell_type": "markdown", "id": "24613643", "metadata": {}, "source": [ "Note that here we gave ID numbers for each PQ in square brackets, both when identifying the parameters to be treated as PQs (``pop_fstar`` and ``pop_fesc``) and when setting the values of their sub-parameters (e.g., ``pq_func[0]``, ``pq_func_par0[0]``, etc. \n", "\n", "To check the result:" ] }, { "cell_type": "code", "execution_count": 7, "id": "e924c668", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEQCAYAAADlK+DYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA8NUlEQVR4nO3de5xN9f748dfbuMu4hAjnkERO5ZLSqaRSp6gTSkhUlBkaP12cSqXbSemgSE3GiEp1dIS+OZKoXFJxjEg0uaYME6JyGczt/ftjr9E27TF7Zl/W3rPfz8djP2bvtT6ftd972PPen8967/URVcUYY4yJNOXcDsAYY4zxxRKUMcaYiGQJyhhjTESyBGWMMSYiWYIyxhgTkSxBGWOMiUjl3Q4gXMqVK6dVqlRxOwxjjPFbVlaWqmrMDiRiJkFVqVKFw4cPux2GMcb4TUSOuB2Dm1zLzCJyrYhsFJEtIjLCx/6WIvKliBwTkX/42B8nImtEZF54IjbGGBNOriQoEYkDkoEuQCvgFhFpVajZfmAYMK6Iw9wDpIcsSGOMMa5yawR1IbBFVbepajbwDtDNu4Gq7lHVVUBO4c4i0gi4Dng1HMEaY4wJP7cSVENgh9fjDGebvyYADwL5J2skIgkikiYiabm5uSUO0hhjjHvcSlDiY5tfV60VkeuBPaq6uri2qpqqqu1VtX358jFTD2KMiSEiXCvCRhG2iODjfD63irDOuX0hQuvi+opQW4RFImx2ftYK1+vx5laCygAaez1uBOzys+8lwA0ish3P1OCVIvJWcMMzxpjIJ4KP8/kUPp//PdBJlfOAp4FUP/qOAD5RpTnwifM47NwaVqwCmotIU2An0Afo609HVX0YeBhARC4H/qGq/UITpjH+UfXcYpWqIuJrYsQUEPHcgsw5n8825zkKzud/W9BAlS+82q/AMyAorm834HKn3RvAEuChoEdfDFcSlKrmishQ4CMgDpimqhtEZLCzP0VE6gNpQDyQLyL3Aq1U9YAbMRtTlJ9+ghYt4EBM/8+05FScadNgwICgH9bX+fwOJ2l/J/ChH31PUyUTQJVMEeoFJ9ySce3EjKrOB+YX2pbidf8nfs/0RR1jCZ7Mboxrdu70JKdbboGWLd2OJvSys7P59NNP+fLLLwHltNPq06xZM0SE7Oxsdu3axc6dGcTFxdGu3flcc83fqFChotthu65t29L0qlNehDSvDamqnik6h9/n80W4Ak+CurSkfd1ilQPGBKigQLRfP+ja1d1YQm3p0qUMHDiQbdu2MXjwYB566CGaNGnyh3bffPMNr7zyCpMn38ChQy35z3/+w7nnnhv+gKPez7mqtD9JA7/O54twHp6v5XRRZZ8ffXeL0MAZPTUA9pT2FQQiZq/xZEywFCSosl4oOmvWLK666irKlSvHkiVLmDRpks/kBHDuuecyadIkFi5cyC+//MIFF1zAO++8E96AY4NzPp+mIlTEcz5/rncDEf4EzAH6q7LJz75zgdud+7cD74fwNRSpjL+ljAm9vDzPz7g4d+MIpRkzZtC/f386dOjAhx9+SHx8vF/9rrrqKr7++mtuvvlm+vXrR+XKlenevXtog40hquSKUOh8PhtEcM7nkwI8DpwKvOIUaeSq0r6ovs6hnwNminAn8CNwczhfVwHRGCk9qlatmtrFYk0ofPopdO4MS5ZAp05uRxN8M2bMoF+/fnTs2JF58+ZxyimnlPgYBw8e5Oqrr2bNmjXMmzePq6++OgSRlj0ikqWq1dyOwy02xWdMgMryCOqrr75iwIABdOzYkfnz55cqOQFUr16dDz/8kLPPPptu3bqxenWx37M3xhKUMYEqq+eg9u3bx0033US9evV49913qVq1akDHq1WrFgsXLqROnTr06tWL3377LUiRmrLKEpQxASqLCSovL49bb72VXbt2MWvWLOrWrRuU49arV48ZM2bwww8/kJCQQKycYjClYwnKmACVxSm+cePG8dFHHzFx4kQuvPDCoB77kksuYdSoUcycOZMpU6YE9dimbLEiCWMCNGsW3HwzrFsHZeGrPps3b+bcc8+la9euzJ49OySXMMrPz6dLly4sW7aMb775hjPPPDPoz1EWWJGEMSYgZWkElZ+fz6BBg6hcuTIvv/xyyK6vV65cOV577TUqVKhAUlKSTfUZnyxBGROgsnQOaurUqSxdupRx48Zx+umnh/S5Tj/9dJ555hkWLlzIzJkzQ/pcJjrZFJ8xAXrjDbjjDti6Fc44w+1oSu+nn36iZcuWtG3blk8//TQsVyfPy8ujQ4cO7Ny5k++++44aNWqE/DmjiU3xGWMCUlam+B5//HGysrKYPHly2JbOiIuLIyUlhd27d/Poo4+G5TlN9LAEZUyAysIU3/r165k6dSpJSUmcddZZYX3u9u3bM2TIEFJSUvjuu+/C+twmslmCMiZAZWEE9cADDxAfH8/IkSNdef4nnniCKlWq2CjKnMASlDEBivYR1MKFC1mwYAEjR47k1FNPdSWGevXq8cADDzBnzhxWrFjhSgwm8liRhDEBGj8e7r8ffvkFatZ0O5qSyc/Pp23bthw8eJD09HQqVarkWiyHDh3izDPP5KyzzmLp0qW2hDxWJGEjKGMCFM1TfLNnz2bdunWMGjXK1eQEcMopp/DEE0/w2Wef8cEHH7gai4kMNoIyJkDPPQcPPwxZWVClitvR+C8/P5/WrVuTm5vL+vXriYuADJuTk0OrVq2oXr06q1evjvlRlI2gjDEBidYR1Jw5c1i/fj2PP/54RCQngAoVKvDII4+wZs0a5s+f73Y4xmWuJSgRuVZENorIFhEZ4WN/SxH5UkSOicg/vLY3FpHFIpIuIhtE5J7wRm7MiaKxSCI/P59//vOftGzZkl69erkdzgn69etHkyZNePrpp+0SSDHOlQQlInFAMtAFaAXcIiKtCjXbDwwDxhXangsMV9WzgYuAJB99jQmb3FwQgXJRNB/x3nvv8c033/DYY49FzOipQIUKFRgxYgQrV67kk08+cTsc4yK33lIXAltUdZuqZgPvAN28G6jqHlVdBeQU2p6pql859w8C6UDD8IRtzB/l5UXX9J6qMmrUKFq0aEHv3r3dDsenO+64g0aNGvH000+7HYpxkVsJqiGww+txBqVIMiLSBGgLrAxOWMaUXG5udE3vffzxx6xdu5aHHnoo4kZPBSpVqsSDDz7IsmXLWLZsmdvhRDQRrhVhowhbRPBxuoSWInwpwjERvE6X0EKEtV63AyLc6+x7UoSdXvu6hvElHedWgvJVmlOiyWYROQWYDdyrqgeKaJMgImkikpZbcKLAmCCLthHUmDFjaNCgAX379nU7lJO66667qFu3LmPHjnU7lIglgo/TJfh1ukSVjaq0UaUNcD6QBbzn1WR8wX5VXKlYcStBZQCNvR43Anb521lEKuBJTm+r6pyi2qlqqqq2V9X25aPpI66JKtE0gvrqq6/4+OOPuffee13/3lNxqlSpwtChQ5k3bx7p6eluhxOpnNMlbFOliNMl7FHlD6dLCukMbFXlh9CFWnJuJahVQHMRaSoiFYE+wFx/OornixFTgXRVfSGEMRrjl2hKUOPGjaN69eokJia6HYpfhgwZQuXKlRk/frzbobikTnkR0rxuCYUaBOV0CZ6/wTMKbRsqwjoRpolQqxTHDJgrCUpVc4GhwEd4ihxmquoGERksIoMBRKS+iGQA9wMjRSRDROKBS4D+wJUista5uTI/agxEzxTf999/z8yZM0lMTIyadZfq1q3LHXfcwfTp09m9e7fb4bjg51xV2nvdUgs1CMLpEioCNwDvem2eBDQD2gCZwPMlOWawuFYYq6rzVfUsVW2mqs8421JUNcW5/5OqNlLVeFWt6dw/oKrLVVVU9TxVbePc7Bt9xjXRMoJ68cUXERHuuSe6vjp43333kZ2dTXJystuhRKKATpc4ugBfqXL8E4Aqu1XJUyUfmIJnKjHsouibG8ZEpmgYQR04cIBp06bRq1cvGjVq5HY4JXLWWWfRrVs3kpOTycrKcjucSOOcLqGpMxLy+3SJl1soNL0nQgOvhz2A9QFFWUqWoIwJUDSMoF5//XUOHjwYdaOnAvfffz/79+/n7bffdjuUiKKKj9MlbBBhsAjO6RLqi+B1uoQMEeKdfVWBq4HCxWZjRPhGhHXAFcB9YXpJJ7CLxRoToL59IS0NNm1yOxLf8vPzadGiBXXq1OHLL790O5xSUVXatWtHbm4u69ati5mLyNrFYo0xAcnNjewpvg8//JAtW7ZE7egJQEQYNmwY69evZ+nSpW6HY8LEEpQxAYr0Kb4XX3yRhg0bctNNN7kdSkD69OnDqaeeysSJE90OxYSJJShjAhTJRRLp6eksWrSIu+++mwoVKrgdTkCqVKlCQkIC77//Pj/8EFHfJzUhYgnKmABF8ggqOTmZihUrMmjQILdDCYohQ4YgIrzyyituh2LCwBKUMQHKy4vMBHXw4EGmT59O7969qVu3rtvhBEXjxo3p3r07U6ZM4ciRI26HY0LMEpQxAYrUIom3336bgwcPkpSU5HYoQZWUlMQvv/zCzJkz3Q7FhJiVmRsToMsvB1WIpOIyVeW8886jYsWKpKWllamybFWlVatWxMfHs3Jl2V5px8rMjTEBicQiieXLl7N+/XqSkpLKVHICT8n5kCFD+N///sfq1avdDseEkCUoYwIUiUUSr7zyCjVr1qRPnz5uhxISt912G1WrVmXSpEluh2JCyBKUMQGKtCKJ3bt3M3v2bAYMGEDVqlXdDickatasSd++ffn3v//NL7/84nY4JkQsQRkToEgrkpg2bRo5OTkMHjzY7VBCasiQIRw5coTp06e7HYoJESuSMCZA550HzZrBe+8V3zbU8vPzadasGU2bNuXTTz91O5yQu+iii/j1119JT08vc+fawIokbARlTIAiqUhi4cKFbN++vcyPngoMHjyYjRs3smzZMrdDMSFgCcqYAEVSkURKSgr16tWje/fubocSFr169aJGjRpMnjzZ7VBMCFiCMiZAkVIkkZGRwbx58xg4cCAVK1Z0O5ywqFq1KrfddhuzZ8/m559/djscE2SWoIwJUKQUSUydOpW8vLwyc909fyUmJpKdnc3rr7/udigmyCxBGROgSJjiy83N5dVXX+Vvf/sbZ5xxhrvBhNlf/vIXLr30UlJTU4mVoq9Y4VqCEpFrRWSjiGwRkRE+9rcUkS9F5JiI/KMkfY0Jp0gokliwYAEZGRkkJia6G4hLEhMT2bx5M4sXL3Y7FBNEriQoEYkDkoEuQCvgFhFpVajZfmAYMK4UfY0Jm0gYQaWmplK/fn3+/ve/uxuIS3r27Ent2rVJTU11O5SwE+FaETaKsEUEHx/2aSnClyIcE6HQh322i/CNCGtFSPPaXluERSJsdn7WCsdrKcytEdSFwBZV3aaq2cA7QDfvBqq6R1VXATkl7WtMOLldJJGRkcEHH3zAwIEDo35RwtKqXLkyt912G3PmzGHv3r1uhxM2Ivj4wI5fH/a9XKFKG1Xae20bAXyiSnPgE+dx2LmVoBoCO7weZzjbQt3XmKBzu0hi2rRp5Ofnc9ddd7kXRAQYNGgQOTk5vPHGG26HEk7OB3a2qVLEh332qOLrw/7JdAMKfpFvAN2DEGuJuZWgfH3l29+zm373FZEEEUkTkbTc3Fy/gzOmJNyc4svLyzteHNG0aVN3gogQrVq1isViiUA/sCuwUITVIiR4bT9NlUwA52e9gCMtBbcSVAbQ2OtxI2BXsPuqaqqqtlfV9uXdPklgyiw3iyQ++ugjduzYQUJCQvGNY0BCQgKbN29maSQtzhWQOuVFSPO6Ff6HDuTDPsAlqrTDM0WYJMJlpQ41BNxKUKuA5iLSVEQqAn2AuWHoa0zQuTmCSk1N5bTTTuOGG25wJ4AI07NnT2rWrFmGiiV+zlWlvdet8AsL5MM+qp62quwB3sMzZQiwW4QGAM7PPaV9BYFwJUGpai4wFPgISAdmquoGERksIoMBRKS+iGQA9wMjRSRDROKL6uvG6zBGFfLz3UlQO3fuPH7liFgtjiisSpUqsXZlCecDO01FKNEHdhGqiVC94D7wN2C9s3sucLtz/3bg/aBG7SfXvgelqvNV9SxVbaaqzzjbUlQ1xbn/k6o2UtV4Va3p3D9QVF9j3JCX5/npxhTfa6+9Rl5eXswXRxQ2aNAgsrOzY2IZDlV8fGBngwiDRXA+7FNfBK8P+2SIEA+cBiwX4Wvgf8AHqixwDv0ccLUIm4GrncdhZ8ttGBOAo0ehShV49ll4+OHwPW9eXh7NmjXjrLPOYuHCheF74ihx8cUXs3///qhfhsOW2zDGlJpbI6hFixbxww8/xNx19/yVkJDAxo0b+eyzz9wOxQTAEpQxASj49kK4z0GlpqZSt25dunWz76j7UrAMR9kplohNlqCMCUDBCCqcCSozM5O5c+cyYMCAmFlWo6SqVq1Kv379mDVrFvv27XM7HFNKlqCMCUDBCCqcU3xWHOGfQYMGcezYMd588023QzGlZAnKmACEe4ovPz+fKVOmcOWVV9K8efPwPGmUat26NR06dIi1K0uUKZagjAlAuIskFi5cyPbt22N2WY2SSkhIID09neXLl7sdiikFS1DGBCDcI6iC4oju3buH5wmjXO/evYmPj2fy5Mluh2JKwRKUMQEIZ5HErl27rDiihKpVq0b//v2tWCJKWYIyJgDhLJKYNm0aeXl59t2nEkpMTOTYsWMxcWWJssYSlDEBCNcUX8GyGp07d+bMM88M7ZOVMeeeey5//etfmTx5shVLRBlLUMYEIFxFEgsXLuSHH36wZTVKKTExkY0bN7Js2TK3QzElYAnKmACEawSVkpLCaaedZsURpdSrVy9q1qxpxRJRxhKUMQEIR5HEjh07mDdvHnfeeacVR5SS9zIce/fudTsc4ydLUMYEIBxFEq+++iqqasURARo8eDDZ2dm89tprbodi/GQJypgAhHoElZOTw5QpU+jSpQtNmjQJzZPEiLPPPptOnToxefJk8vPz3Q7H+MESlDEBCPUIat68eWRmZjJ48ODQPEGMGTJkCNu2bbM1tKKEJShjAhDqIomUlBQaN25M165dQ/MEMaZHjx7Uq1ePSZMmuR2K8YMlKGMCEMopvs2bN7Nw4ULuuusu4txYU74MqlixInfeeSfz5s1jx44dbodjimEJypgAhHKKb9KkSZQvX96KI4IsISEBVWXKlCluhxIUIlwrwkYRtogwwsf+liJ8KcIxEf7htb2xCItFSBdhgwj3eO17UoSdIqx1bq4M4V1LUCJyrYhsFJEtIuLjlyoiIhOd/etEpJ3XvvtEZIOIrBeRGSJSObzRG+MRqhFUVlYWr732GjfddBMNGjQI7sFjXJMmTejatSupqalkZ2e7HU5ARIgDkoEuQCvgFhFaFWq2HxgGjCu0PRcYrsrZwEVAUqG+41Vp49zmh+YVnJwrCUpEfPxSpfAvtQvQ3LklAJOcvg3x/LLbq+o5QBzQJ0yhG3OCUI2gZsyYwa+//srdd98d3AMbAJKSkti9ezdz5sxxO5RAXQhsUWWbKtnAO0A37waq7FFlFZBTaHumKl859w8C6UDD8ITtH7dGUM4vVbepqs9fqvN4unqsAGqKSMFHyfJAFREpD1QFdoUrcGO8haJIQlVJTk7mnHPOoWPHjsE7sDnummuuoVmzZiQnJ7sdSqAaAt4n0zIoRZIRoQnQFljptXmoCOtEmCZCrYCiLCW3EpQ/v1SfbVR1J56h6o9AJvCbqlrNqHFFKKb4Vq5cyZo1a7j77rsRkeAd2BxXrlw5hgwZwvLly1m3bp3b4ZxEnfIipHndCl+M0dd/kBJdEVeEU4DZwL2qHHA2TwKaAW3w/J19vmRxB4dbCcqfX6rPNiJSC8/oqilwOlBNRPr5fBKRBBFJE5G03IKPusYEUSim+JKTk6levTr9+vn8b22CZMCAAVSuXDnCR1E/56rS3uuWWqhBBtDY63EjSjCjJEIFPMnpbVWOz3eqsluVPFXygSl4Zr3Czq0E5c8vtag2VwHfq+peVc0B5gAX+3oSVU1V1faq2r58uJY8NTEl2COo3bt3M3PmTG6//XaqV68enIMan2rXrk3fvn156623+PXXX90Op7RWAc1FaCpCRTzn4+f601EEAaYC6aq8UGifd2VOD2B9kOItEbcSlPNLlaYiUtQvdS5wm1PNdxGeqbxMPFN7F4lIVfHMf3TGc3LPmLAL9giqoLJs6NChwTmgOamkpCSysrJ4/fXX3Q6lVFTJBYYCH+H5OzhTlQ0iDBZhMIAI9UXIAO4HRoqQIUI8cAnQH7jSRzn5GBG+EWEdcAVwX7hfG4C4tYCXiHQFJuCpwpumqs+IyGAAVU1xks/LwLVAFjBAVdOcvk8BvfGUSa4B7lLVYyd7vmrVqunhw4dD9XJMjJo0Ce6+GzIzoX79wI6Vk5PDn//8Z8477zwWLFgQnABNsS699FIyMzPZtGlTxH0hWkSyVLWa23G4xbV5L1WdDyfW1qtqitd9BZKK6PsE8ERIAzTGD8Gc4ps9ezaZmZll5guk0WLYsGH07t2bDz/8kOuvv97tcIwXu5KEMQEI5hTfSy+9RLNmzejSpUvgBzN+69GjBw0bNuTFF190OxRTiCUoYwIQrBHU6tWr+eKLLxg6dCjlytnbMpwqVKhAUlISH3/8Md9++63b4Rgv9k4wJgDBGkFNnDiRatWqMWDAgMCDMiU2aNAgKlWqxEsvveR2KGWKCDNEqOn1uJYIb/vb3xKUMQEIxpUkMjMzmTFjBgMGDKBGjRrBCcyUSJ06dbj11luZPn06+/fvdzucsqSlKr8WPFDlF/jDtQKLZAnKmAAEY4ovOTmZ3Nxc7rnnnuIbm5AZNmwYWVlZVqQSXOW8L5Mkwql4Krf96xzos4tIBV/3jYkFBSOo0p42OnLkCCkpKdxwww2ceeaZwQvMlFjr1q3p3LkzEydOjPqrnEeQccAXIjwtwtPAcuA5fzsHlKCc7yp9LiK1RKQbMCaQ4xkTbfLyAhs9vfnmm+zbt4/77nPle5CmkPvvv59du3Yxc+ZMt0MpE1R5E7gJ2A38BPRQ5d/+9g/4i7oicgeeSw1dCHRV1Yi8srh9UdeEwkMPwYsvwtGjJe+rqvzlL3+hcuXKrF692i4MGwHy8/M555xzIubfJNq/qCvCIOA/qhwQYQLQDhipyjJ/+pd6BCUibUXkNeByYAAQD4wSkWmlPaYx0SY3t/QjqAULFpCens59993n+h9C41GuXDnuu+8+1qxZw9KlS90Opyz4f05yuhxoCTxMCa6MHsgU3zY8lyJ6F/gZOOQ8juRLAxsTVIFM8Y0dO5aGDRvSu3fv4AZlAtKvXz/q1KnD88+7ssJEWZPv/LwaeF2VzynBFYxKnaBU9TdVXQ3ciGeF23VAI2ebMTEhN7d034FKS0tj8eLF3HvvvVSsWDH4gZlSq1KlCklJScybN4/0dLsOdYB2iDAZz7VTPxShMiXIO8U2FJHPRSS+iH0CbFXVd4EHgVP9fWJjyoLSjqDGjh1LfHw8CQmF158zkSApKYkqVaowduxYt0OJdr3xXGn9GlV+A2oD//C3sz+Z7HygCoCIbBCR498kdJZjf9a5v0tVXytJ5MZEu9KMoLZt28asWbMYMmQI8fE+P/sZl9WtW5eBAwfy1ltvkZGR4XY40SwbOA2OrwRcCTjpyhPe/ElQ24C/O/fPdp7AGEPpiiReeOEF4uLiGDZsWGiCMkExfPhw8vPz7SKygUnFs/p5QQ7ZB4z3t7M/CeoJ4GUR2YlnWfYRItJTRJqVNFJjypqSTvHt3buXadOm0b9/f04//fTQBWYC1rRpU3r16sXkyZOjecVdt7VR5UGcUZMqBwhmkYRzfunPwOOAAJcC04BNIvKriCwRkfEi0r800RsTzUo6xTdhwgSOHj3KAw88ELqgTNA8+OCDHDx4kEmTJrkdSrTKFiEOz+AGEU4D8vzt7Fc1haruVtWpwDdAF1WNxzPdlwiswHPxv3ElDNyYqFeSEdSvv/7Kyy+/zE033UTLli1DG5gJijZt2nDNNdcwYcIEsrKy3A4nIonwubOEvC8TgPeBeiI8C3wGPOvvsUtUZq6qrVV1n3N/k6r+R1VHqOo1qnpaSY5lTFlQkhFUcnIyBw4c4JFHHgltUCaoRo4cyZ49e+wiskXzKqRjgwhehXS8AwwHRgO7gG6qzPL3wHY1c2MC4G+RxOHDh5kwYQJdu3albdu2oQ/MBM2ll15Kp06dGDNmDMeO+V2AFjYiXCvCRhG2iDDCx/6WInwpwjGRE0u8i+orQm0RFomw2flZq/BxvZy0kE6Vjaokq/KyKiX6YpklKGMC4O8U35QpU/j5559t9BSlRo4cya5du3j99dfdDuUEzvmdZKALnlMtt4j8Yb2l/XgupjCuBH1HAJ+o0hz4xHlcFKeQDq9COnqKEHAhnSUoYwLgzxTf0aNHGTt2LJ06deKSSy4JT2AmqDp37kyHDh147rnnyMnJcTscbxcCW1TZpko28A7QzbuBKntUWQUUDvxkfbsBbzj33wC6FxWAKicppONXEZaIMF6EEhfSuZagRORaEdkoIltExMewVEREJjr714lIO699NUVkloh8JyLpIvLX8EZvjIc/I6gpU6awa9cunnjiifAEZYJORBg5ciTbt2/n7bf9XrE8HBoCO7weZzjbAu17miqZAM7Peic7kCq7VfEqpCMohXSuJCgR8TG0lMLD0i5Ac+eWAHjXeb4ILFDVlkBrKNm8pjHBUtwI6ujRo4wePZrLLruMyy+/PGxxmeC77rrraNeuHaNGjQrjKKpOeRHSvG6Fr43l6zL4/q6hFEhfn1RprYpTSMcmVf6jyghVrlGlxIV0bo2gnKGlblNVn8NS5/F053JKK4CaItLAuS7gZcBUAFXNVtVfwxi7MccVVyQxZcoUMjMzefLJJ21JjSgnIjz11FNs3bqVN998M0zP+nOuKu29bqmFGmQAjb0eN8JTLeePk/XdLUIDAOfnnpLHHji3EpQ/w9Ki2pwB7AVeE5E1IvKqiPhc0EtEEkQkTUTScgvW5jYmiE42xXfkyBFGjx5Np06duOKKK8IbmAmJ6667jgsuuICnn346UpaFXwU0F6GpCBWBPsDcIPSdC9zu3L8dz3eZws6tBOXP0LKoNuXxrMo4SVXbAocposJEVVNVtb2qti8fyLrcxhThZFN8qampx0dPpmwoGEVt3749Iir6VMkFhuK5Yng6MFOVDSIMFmEwgAj1RcgA7gdGipAhQnxRfZ1DPwdcLcJmPGs5PRfeV+YR8JLvpXpST1HDk6p6jfP4YQBVHe3VZjKwRFVnOI834lm9V4EVqtrE2d4RGKGq153sOW3JdxMK558Pp58O//3vidsPHTpEs2bNaNWqFYsXL3YnOBMSqsrFF1/Mzp072bx5M5Uqhe762dG+5Hug3BpBOUNLaSoiRQ1L5wK3OdV8FwG/qWqmqv4E7BCRFk67zsC3YYvcGC9FjaAmTpzInj17ePZZv6/qYqKEiPDPf/6THTt2MHnyZLfDKdNcGUEBiEhXPNdpigOmqeozIjIYQFVTnMUQXwauBbKAAaqa5vRtA7wKVMTzLeYBqvrLyZ7PRlAmFM45B1q2hFleF2/Zv38/Z5xxBpdddhlz5/p7OsBEE1Wlc+fOrF+/nq1bt1K9evWQPE+sj6BcOzGjqvOB+YW2pXjdVyCpiL5rgfahjM8Yf/iq4hszZgwHDhzgmWeecScoE3IiwujRo7nooosYP348jz/+uNshlUl2JQljAlB4ii8zM5OJEyfSt29fzj33XPcCMyHXoUMHbrzxRsaNG8fevXvdDqdMsgRlTAAKl5k/+eST5OTk8NRTT7kXlAmbUaNGcfjwYTvXGCKWoIwJgPcIasOGDbz66qvcfffdNGtmC07HgrPPPpsBAwbwyiuv8P3337sdTpljCcqYAHiPoB588EGqV6/OY4895m5QJqyeeuopypcvz4gRJ7vgtykNS1DGBKCgSOLjjz9m/vz5jBw5kjp16rgdlgmjhg0b8sADDzBz5ky++OILt8MpU1wrMw83KzM3oVC7NvTtq3z+eTt+/fVX0tPTqVy5stthmTA7fPgwzZs3509/+hNffvll0K67GOtl5jaCMiYAeXmwceMG1q5dy7PPPmvJKUZVq1aNZ555hpUrV/Kf//zH7XDKDBtBGROAatUU1Um0a/dvPvvsM7tieQzLy8ujffv27Nu3j++++46qVasGfEwbQRljSu3YsVyOHDnIxIkTLTnFuLi4OCZOnMiOHTsYPXp08R1MsSxBGVNK3377LXl5Qtu259KuXbviO5gyr2PHjvTt25exY8eydetWt8OJepagjCkFVeWee+4FynPllZ3cDsdEkLFjx1KhQgXuu+8+t0OJepagjCmFd999l48//gSAGjVi9hSB8eH000/n8ccf57///S8ffPCB2+FENSuSMKaEfvvtN1q2bEmDBk1Ys+ZLnnkGHnnE7ahMJMnOzqZ169YcO3aM9evXl7pgwookjDElMnLkSPbs2cPEiclA0Uu+m9hVsWJFUlJS+P777+26jAGwBGVMCaxatYrk5GSSkpI47zxPYYQlKONLp06dGDhwIM8//zxff/212+FEJUtQxvgpJyeHhIQE6tevz9NPP01urme7rxV1jQFPwUTt2rVJTEwkLy8vJM8hwrUibBRhiwh/uCCgCCLCRGf/OhHaOdtbiLDW63ZAhHudfU+KsNNrX9eQBF8MS1DG+GnMmDGsXbuW5ORkatSoQcHfGxtBmaLUrl2b8ePHs3LlSqZMmRL044sQByQDXYBWwC0itCrUrAvQ3LklAJMAVNmoShtV2gDn41m5/D2vfuML9queuLhsuNhbyxg/bNiwgX/+85/06tWLHj16ANgIyvilb9++/PTTT/Ts2TMUh78Q2KLKNgAR3gG6Ad96tekGTFdFgRUi1BShgSqZXm06A1tV+SEUQZaWjaCMKUZeXh4DBw6kevXqvPTSS17bPT9tBGVORkQYPnx4qK5y3xDY4fU4w9lW0jZ9gBmFtg11pgSniVArGMGWlCUoY4rxwgsv8L///Y+XXnqJevXqHd9eMIKyBGVCp055EdK8bgmFGvi6vlbh7w6dtI0IFYEbgHe99k8CmgFtgEzg+ZJGHgyuJSgRuVZENorIFhHxcWJPREQmOvvXiUi7QvvjRGSNiMwLX9Qm1qxbt46RI0fSvXt3+vTpc8I+m+IzofdzrirtvW6phRpkAI29HjcCdpWwTRfgK1V2F2xQZbcqearkA1PwTCWGnSsJSkR8nNgTv07sebkHSA9xqCaGHTt2jH79+lGrVi1SU1P/cDFYm+IzEWAV0FyEps5IqA8wt1CbucBtTjXfRcBvhc4/3UKh6T0RGng97AGsD37oxXPrreWc2FPnxJ6c5MSeOif2pKaINFDVTBFpBFwHPAPcH+bYTYx47LHH+Oabb5g3bx5169b9w34bQRm3qZIrwlDgIyAOmKbKBhEGO/tTgPlAV2ALnkq9AQX9RagKXA0kFjr0GBHa4JkK3O5jf1i4laB8nbTr4EebhnjmQycADwLVT/YkIpKAZ/RFxYoVAwrYxJbFixczbtw4EhMTue6663y2sRGUiQROCfj8QttSvO4rkFRE3yzgVB/b+wc5zFJx6xxUqU/sicj1wB5VXV3ck6hqqqq2V9X25e2viPHTnj17uPXWWznrrLMYN25cke2sSMKY0HLrrRXIib2ewA0i0hWoDMSLyFuq2i+E8ZoYkZ+fz+23387+/fv58MMPOeWUU4psa1N8xoSWWyMo58SeNBWRYk7siYiIc2JPM1X1YVVtpKpNnH6fWnIywfL888+zYMECxo8fT+vWrU/a1qb4jAktV95aqporIoVO7OkGEXFO7OlJT+wZEwrLly/nkUce4aabbmLw4MHFtrcRlDGhZetBGQNkZmbSrl07qlevzqpVq6hRo0axfZYsgSuugMWL4fLLQx6iiUGxvh6UTU6YmJednc3NN9/MgQMHWLRokV/JCaxIwphQs7eWiXn3338/n3/+Oe+88w7nnHOO3/1sis+Y0LJr8ZmYlpKSQnJyMsOHD6d3794l6mtFEsaEliUoE7M+/fRThg4dSteuXfnXv/5V4v42gjImtCxBmZi0adMmevbsydlnn82MGTOIK0WWsRGUMaFlCcrEnN27d9OlSxfi4uL473//S3x8fKmOY0USxoSWvbVMTDl06BDXX389mZmZLF68mCZNmpT6WDbFZ0xoWYIyMSMnJ4ebb76ZNWvW8H//93906FD4+sQlY1N8xoSWvbVMTMjLy+OOO+5gwYIFTJkyheuvvz7gY9oIypjQsnNQpsxTVQYPHsy///1vRo8ezV133RWU49oIypjQsgRlyjRV5f777+fVV1/l0UcfZcSIEUE7thVJGBNalqBMmaWqPPDAA0yYMIFhw4bx9NNPB/X4BSMom+IzJjTss58pkwpGThMmTGDo0KFMmDABEV9rYJaejaCMCS0bQZkyJz8/n2HDhjFhwgTuueceJk6cGPTkBFYkYUyoWYIyZUpOTg633347L7/8MsOHD2f8+PEhSU5gRRImMohwrQgbRdgiwh9OsoogIkx09q8ToZ3Xvu0ifCPCWhHSvLbXFmGRCJudn7XC9Xq8WYIyZUZWVhY9evTgrbfeYtSoUYwdOzZkyQlsis+4T4Q4IBnoArQCbhGhVaFmXYDmzi0BmFRo/xWqtFGlvde2EcAnqjQHPnEeh50lKFMm7N27l6uvvpr58+eTkpLCo48+GtLkBFYkYSLChcAWVbapkg28A3Qr1KYbMF0VVWUFUFOEBsUctxvwhnP/DaB7EGP2myUoE/U2bdrEX//6V7766iveffddEhMTw/K8BSOocvYuMiFTp7wIaV63hEINGgI7vB5nONv8baPAQhFWFzr2aapkAjg/6wX6SkrDJidMVFu6dCk33ngjcXFxLF68mIsuuihsz52b6xk9hXigZmLaz7mFpt4K8/W/T0vQ5hJVdolQD1gkwneqLCtNpKFgn/1MVFJVkpOTueqqq6hXrx4rVqwIa3ICzxSfnX8yLssAGns9bgTs8reN6vGfe4D38EwZAuwumAZ0fu4JeuR+cC1Bici1IrJRRLaIiI/KExERmejsXyci7ZztjUVksYiki8gGEbkn/NEbNx09epRBgwYxdOhQunTpwooVKzjjjDPCHkduriUo47pVQHMRmopQEegDzC3UZi5wm1PNdxHwmyqZIlQToTqACNWAvwHrvfrc7ty/HXg/1C/EF1feXiJSUHlyNZ7svkpE5qrqt17NvCtPOuCpPOkA5ALDVfUrEakOrBaRRYX6mjJq69atx69I/thjj/Hkk09SzqWTQHl5ViBh3KVKrghDgY+AOGCaKhtEGOzsTwHmA12BLUAWMMDpfhrwnjNFXR74tyoLnH3PATNFuBP4Ebg5PK/oRG59/nMqT3QbgIgUVJ54Jxmn8kQVWCEiNUWkgapmQsHJOz0oIul4TvhZgirjZs+ezcCBA48vNBiMK5IHwkZQJhKoMh9PEvLeluJ1X4EkH/22Aa2LOOY+oHNwIy05t6b4Aq08AUBEmgBtgZXBD9FEikOHDjFo0KDjS7SvWbPG9eQEvxdJGGNCw60EFWjlCSJyCjAbuFdVD/h8EpEEEUkTkbTcgppgE1VWrlxJ27ZtmTp1Kg8//DDLli3jz3/+s9thAVYkYUyouZWgAqo8EZEKeJLT26o6p6gnUdVUVW2vqu3L21+SqHLkyBEefPBBLr74YnJycli6dCnPPvssFStWdDu042yKz5jQcitBOZUn0lREiqk8ERERp/JEM8VzeYCpQLqqvhDesE04LF++nDZt2jB27FjuvPNOvv76azp27Oh2WH9gRRLGhJYrCUpVc+F45Uk6MFNVN4jIYBEZ7DSbD2zDU3kyBbjb2X4J0B+4UkTWOreu4X0FJhR+/vln7rzzTjp27MixY8dYtGgRqamp1KhRw+3QfLIRlDGhJZ4iubKvWrVqevjwYbfDMD7k5eUxZcoUHn30UQ4cOMDw4cN57LHHqFatmtuhnVTv3vD11/Ddd25HYsoqEclS1ch+I4SQff4zrlq8eDH33nsv69at47LLLiM5OZlzzjnH7bD8YkUSxoSWXerIuGLdunVcd911XHnllfz222+8++67LFmyJGqSE9gUnzGhZgnKhNXmzZvp378/bdq04YsvvuC5554jPT2dnj17hnx5jGCzIgljQss+/5mw2Lx5M6NGjeKtt96iUqVK/OMf/+Dhhx+mVi1XFuoMChtBGRNa9vYyIbV69Wr+9a9/MWvWLCpXrsy9997LAw88QP369d0OLWB2JQljQssSlAm6vLw8PvjgAyZMmMDixYupUaMGI0aM4J577uG0005zO7ygsSIJY0LL3l4maPbt28drr71GSkoKW7dupXHjxowZM4bExETi4+PdDi/obIrPmNCyt5cJSH5+PkuWLGHatGnMmjWLY8eOcemll/Lcc8/RvXt3yvIlpvLyoFIlt6Mwpuwqu389TEht2rSJt99+m+nTp7N9+3Zq1KjBoEGDSExMjKpS8UDYCMqY0LK3l/HbDz/8wOzZs5kxYwZpaWmICJ07d+aZZ56hR48eVKlSxe0Qw8qKJIwJLUtQ5qS+++473n//febMmcP//vc/ANq2bcu4cePo06cPDRsWXsYrdliRhDGhZW8vc4Jjx47x2WefMX/+fD744AM2bdoEwPnnn8/o0aPp2bMnZ555pstRRgab4jMmtOztFePy8/NZv349H3/8MYsWLWLZsmVkZWVRqVIlOnXqxLBhw7jhhhto3Lhx8QeLMXYlCRMJRLgWeBGIA15V5blC+8XZ3xXIAu5Q5SsRGgPTgfpAPpCqyotOnyeBQcBe5zCPOEvLh5UlqBiTk5PD119/zfLly1m2bBlLly5l//79ALRo0YKBAwdyzTXXcMUVV0T81cTdZiMo4zYR4oBk4Go8i7yuEmGuKt96NesCNHduHYBJzs9cYLiTrKoDq0VY5NV3vCrjwvVafLG3VxmmquzcuZNVq1axcuVKVqxYwapVq8jKygKgadOmdOvWjcsvv5wrrrjCRkklZEUSJgJcCGxRZRuACO8A3eCEBNUNmK6KAitEqClCA1UygUwAVQ6KkA40LNTXVZagyoi8vDy2bdvG119/zdq1a1mzZg2rV69m9+7dAFSoUIE2bdpw5513cumll3LJJZfEdIFDMFiRhIkADYEdXo8z8IyOimvTECc5AYjQBGgLrPRqN1SE24A0PCOtX4IXtn/s7RVl8vLy2L59O9999x3ffvst3377LRs2bGD9+vUcOXIEgLi4OFq1asU111zDBRdcwPnnn0+bNm1irgw81GyKz4RenfIipHltSFUl1euxryUACq9Ce9I2IpwCzAbuVeWAs3kS8LTT7mngeWBgCYMPmL29IlBOTg4//vgj27ZtY+vWrWzdupXNmzezefNmtmzZQnZ29vG2DRo04OyzzyYxMZFzzz2X8847j3POOYfKlSu7+ApigxVJmND7OVeV9idpkAF4z803Anb520aECniS09uqzClooMrugvsiTAHmlSr8AFmCCrP8/Hz27dtHRkYGO3fuJCMjgx9//PH4bfv27ezcuZP8/PzjfSpVqkSzZs1o3rw5119/PS1atKBFixa0atUqqperiHY2gjIRYBXQXISmwE6gD9C3UJu5eKbr3sEz/febKplOdd9UIF2VF7w7eJ2jAugBrA/liyiKvb2CIC8vj/3797N3797jt927d7Nnzx52797NTz/9xE8//URmZiaZmZnk5OSc0L98+fI0atSIP/3pT1xxxRU0adKEJk2acMYZZ3DGGWfQsGFDypWztSUjjRVJGLepkivCUOAjPGXm01TZIMJgZ38KMB9PifkWPGXmA5zulwD9gW9EWOtsKygnHyNCGzxTfNuBxLC8oEJEtfB0ZZieWKRQ7b4Wqt0XH7X7+pU/fX2pVq2aHj582O/4jh7N4dFHV5OVdZisrCwOHTrM4cOHOXz4EIcO/X47ePAQnuP6+j0K8fHVqVmzFjVr1qRmzZrUrl2LWrVqU6tWLerUOZXatWtTs2ZNypWzv3TRpn9/GDoUxrlaiGvKMhHJUtWY/b6HKyMoEfFRuy9zVbXY2n0/+wbs6NE8XnjhooCPc+CA5/bjj0EIykScunXdjsCYssutKT6ndl+d2n05Se2+OrX7UlNEGgBN/OgbsPj4ysyevZH4+HiqV69O1apVbZrNnKBcOWjRwu0ojCm73EpQgdTu+9M3YOXKwY032l8fY4xxi1sJKpDafX/6eg4gkgAkAFSsWLEk8RljjHGZWwkqkNr9in70BUBVU8HzpbZq1aq5Uw1ijDGmVNw6qeLU7ktTEamIp3Z/bqE2c4HbxOMi4DdVzfSzrzHGmCjnyghKVXNFpFDtvm4QEad2X4us3S+qrwsvwxhjTAi59j2ocCvp96CMMcZtsf49KKubNsYYE5EsQRljjIlIMTPFJyL5wBGvTXFAXqH7vraVx7PyZEl5H8vf/cVtK+6+97bSxF2amIvaXpq4w/m79rX9ZI8t7uLjKm5/aeKOpvek9+NgxV1FVWN3IKGqMXkDUgvfL2JbWqDH93d/cduKu19oW4njLk3MwYw7nL9rX9tP9tjidifuaHpPhiruWL7FbmaG//q472tbMI7v7/7ithV3342Yi9oejXGf7LHFXfTz+bu/NHFH03vS+3Ew445ZMTPFV1oikqaqJ1swLCJFY9zRGDNY3OFmcceOWB5B+Su1+CYRKRrjjsaYweION4s7RtgIyhhjTESyEZQxxpiIZAnKGGNMRLIEZYwxJiJZgioBEekoIiki8qqIfOF2PP4SkT+JyFwRmSYiI9yOx18i0kpEZorIJBHp6XY8xRGRM0RkqojM8tpWTUTeEJEpInKrm/H5UkTMf9gWaYqIu7vze35fRP7mZnxFKSLus52/K7NEZIib8UUct7+I5fYNmAbsAdYX2n4tsBHP1dRHFNrXHUiMlriBqwriBaZHUdzDgY7O/bmRHq/Xvlle9/sDf3fu/ycaYj7ZtiiJuxYwNQrjLhfOuKPhZiMoeB3Pf6TjRCQOSAa6AK2AW0SklVeTvsCMcAVYhNfxP+41QB8R+RRYHOY4C3sd/+N+E0/cY4FTwxxngdcp+f8Pb42AHc79k11mJ5heJ7CY3fI6wYl7pNMnXF4nwLhF5AZgOfBJ6MKMPjGfoFR1GbC/0OYLgS2quk1Vs4F3gG7gmS7Ds3jigfBGeqISxj0AeEJVrwSuC2+kJypJ3Kq6R1WTgBHAz2EOFSj5/w8fMvAkKQjT+y0IMbsi0LidxU3/BXyoql+FNtrfBeP3rapzVfViIOKmgd0U8wmqCA35/VMveP7INHTu3wm8FvaI/FNU3AuAYSKSAmx3Ia7i+IxbRJqISCowHRjrSmS+FRXvqc7vuK2IPOzsmwPcJCKTcPeSN37HXMTrcEtJftf/D890ds+CxU9dVJLf9+UiMlFEJuNZqNU4XFlRNwqIj20KoKpPhDmWkvAZt6quByK5yKCouLcDCWGOxR9FxbsPGFxo42Gc1aBdVpKY/7DNRSWJeyIwMSxRFa8kcS8BloQhpqhjIyjfMoDGXo8bAbtciqUkLO7wiLZ4ITpjBos7plmC8m0V0FxEmopIRaAPMNflmPxhcYdHtMUL0RkzWNyxze0yQrdveKrxMoEcPJ967nS2dwU2AVuBR92O0+K2eMtyzBa33Xzd7GKxxhhjIpJN8RljjIlIlqCMMcZEJEtQxhhjIpIlKGOMMRHJEpQxxpiIZAnKGGNMRLIEZYwxJiJZgjLGGBORLEGZmCMi74qIisi/fex70Nm3rwTHe0pE5oTq+MbEKruShIk5IvK9czdLVf/itb0hkA78Bnyrqtf4ebx1wFhVfTMUxzcmVtkIysQUETkVaAJMBc4SkUpeu8cD7+G5ptoqP4/XDDgbmBeK4xsTyyxBmVhzgfPzDSAfz3LciMhVwN+A54CmQJqfx7sRWKKqv4To+MbELEtQJta0BzJUdQewHmjtLIfwMvAYv6/hswqOr3Y64yTH64FnVFSq4xtjimYr6ppYcwG/j17WAOcBpwNHgVeAh4CfVHWn06aN0+4PRKQBcCFwcwDHR0Quw7NycG1gIfCSquYF8iKNKQtsBGViTXt+H718hWfa7REgyUkK7Tlx+q01UEdElorIbhH5u9e+7kCad7Ip6fFF5EbgfuBJ4FYgD8/6QsbEPEtQJmY4I57TOXGE8xdglqp+7mzzTjDgGUH9qqqdgAHAbV77bsRreq+Uxx8O3KKqW1T1F1V9CfhRRK4M5LUaUxZYgjKxpKCAoSCBrATq4pleQ0Tq4TlHlOY8Lg/UA/7ltC8P7HP21QIu58TzTyU9fg0gU1WPiMj1IrLE6bcEODfA12pM1LNzUCaWtAe+V9X9AKqaD/xcaD/8nmBaAd94nQ86D/jGuX89sElVNwVw/ANAfRERYDGwztneEthWqldoTBliX9Q1pggi0h84R1Ufch7PxFPA8JmIvAdsUNWRAT7HI0A14DFVzReRi4AxwFWqmh3gSzAmqtkIypiitebECr7z+H2U8yUw5w89Sm40niKJzz0DKb4HbrbkZIyNoIwxxkQoK5IwxhgTkSxBGWOMiUiWoIwxxkQkS1DGGGMikiUoY4wxEckSlDHGmIhkCcoYY0xEsgRljDEmIv1/c/jprKTw25oAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pop = ares.populations.GalaxyPopulation(**pars)\n", "\n", "Mh = np.logspace(7, 13, 100)\n", "\n", "fig, ax1 = pl.subplots(num=2)\n", "\n", "ax1.semilogx(Mh, pop.fstar(z=6, Mh=Mh), color='k')\n", "ax1.set_ylabel(r'$f_{\\ast}$')\n", "ax1.set_xlabel(r'$M_h / M_{\\odot}$')\n", " \n", "ax2 = ax1.twinx()\n", "ax2.tick_params('y', colors='b')\n", "ax2.set_ylabel(r'$f_{\\mathrm{esc}}$', color='b')\n", "ax2.semilogx(Mh, pop.fesc(z=6, Mh=Mh), color='b')" ] }, { "cell_type": "markdown", "id": "65f292e9", "metadata": {}, "source": [ "## Allowed Parameters\n", "\n", "An incomplete list so far:\n", "\n", "* ``pop_fstar``\n", "* ``pop_fesc``\n", "* ``pop_focc``" ] } ], "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 }