diff --git a/include/background.h b/include/background.h index 809c12601..602114a6c 100644 --- a/include/background.h +++ b/include/background.h @@ -89,6 +89,10 @@ struct background /* the following parameters help to define the analytical ncdm phase space distributions (p-s-d) */ double * ncdm_psd_parameters; /**< list of parameters for specifying/modifying ncdm p.s.d.'s, to be customized for given model (could be e.g. mixing angles) */ + int output_ncdm_binning; + int collective_ncdm; /**< regards all neutrinos as having the same distribution function, + thus solving only one hierarcy */ + int collective_ncdm_N; /**< amount of collective ncdm species */ double * M_ncdm; /**< vector of masses of non-cold relic: dimensionless ratios m_ncdm/T_ncdm */ double * m_ncdm_in_eV; /**< list of ncdm masses in eV (inferred from M_ncdm and other parameters above) */ double * Omega0_ncdm, Omega0_ncdm_tot; /**< Omega0_ncdm for each species and for the total Omega0_ncdm */ @@ -189,6 +193,8 @@ struct background int index_bg_rho_ncdm1; /**< density of first ncdm species (others contiguous) */ int index_bg_p_ncdm1; /**< pressure of first ncdm species (others contiguous) */ int index_bg_pseudo_p_ncdm1;/**< another statistical momentum useful in ncdma approximation */ + int index_bg_q_ncdm1; + int index_bg_w_ncdm1; int index_bg_rho_tot; /**< Total density */ int index_bg_p_tot; /**< Total pressure */ diff --git a/notebooks/grand_neutrinos.ipynb b/notebooks/grand_neutrinos.ipynb new file mode 100644 index 000000000..2d935b49f --- /dev/null +++ b/notebooks/grand_neutrinos.ipynb @@ -0,0 +1,530 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "db2fd51a", + "metadata": {}, + "source": [ + "# Test of the \"grand neutrino\" collectivisation implementation" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "7cc69add", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from classy import Class\n", + "from time import perf_counter\n", + "from copy import deepcopy\n", + "from scipy.interpolate import interp1d\n", + "def interp_onto(xdata, ydata, xonto):\n", + " interp = interp1d(xdata, ydata)\n", + " return interp(xonto)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "7664519b", + "metadata": {}, + "outputs": [ + { + "ename": "CosmoComputationError", + "evalue": "\n\nError in Class: perturbations_init(L:1016) :error in perturbations_solve(ppr, pba, pth, ppt, index_md, index_ic, index_k, pppw[thread]);\n=>perturbations_solve(L:3335) :error in generic_evolver(perturbations_derivs, interval_limit[index_interval], interval_limit[index_interval+1], ppw->pv->y, ppw->pv->used_in_sources, ppw->pv->pt_size, &ppaw, ppr->tol_perturbations_integration, ppr->smallest_allowed_variation, perturbations_timescale, ppr->perturbations_integration_stepsize, ppt->tau_sampling, tau_actual_size, perturbations_sources, perhaps_print_variables, ppt->error_message);\n=>evolver_ndf15(L:187) :error in initialize_jacobian(&jac,neq,error_message);\n=>initialize_jacobian(L:1574) :could not allocate jac->Ci with size -8", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mCosmoComputationError\u001b[0m Traceback (most recent call last)", + "Input \u001b[0;32mIn [20]\u001b[0m, in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 26\u001b[0m reference\u001b[38;5;241m.\u001b[39mset(standard_settings \u001b[38;5;241m|\u001b[39m {\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mNumber of momentum bins\u001b[39m\u001b[38;5;124m'\u001b[39m: \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m10000, 10000, 10000\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mQuadrature strategy\u001b[39m\u001b[38;5;124m'\u001b[39m: \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m2, 2, 2\u001b[39m\u001b[38;5;124m'\u001b[39m})\n\u001b[1;32m 27\u001b[0m tic \u001b[38;5;241m=\u001b[39m perf_counter()\n\u001b[0;32m---> 28\u001b[0m \u001b[43mreference\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcompute\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 29\u001b[0m toc \u001b[38;5;241m=\u001b[39m perf_counter()\n\u001b[1;32m 30\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mComputed reference in \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mtoc\u001b[38;5;250m \u001b[39m\u001b[38;5;241m-\u001b[39m\u001b[38;5;250m \u001b[39mtic\u001b[38;5;132;01m:\u001b[39;00m\u001b[38;5;124m.5\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m s\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n", + "File \u001b[0;32mclassy.pyx:398\u001b[0m, in \u001b[0;36mclassy.Class.compute\u001b[0;34m()\u001b[0m\n", + "\u001b[0;31mCosmoComputationError\u001b[0m: \n\nError in Class: perturbations_init(L:1016) :error in perturbations_solve(ppr, pba, pth, ppt, index_md, index_ic, index_k, pppw[thread]);\n=>perturbations_solve(L:3335) :error in generic_evolver(perturbations_derivs, interval_limit[index_interval], interval_limit[index_interval+1], ppw->pv->y, ppw->pv->used_in_sources, ppw->pv->pt_size, &ppaw, ppr->tol_perturbations_integration, ppr->smallest_allowed_variation, perturbations_timescale, ppr->perturbations_integration_stepsize, ppt->tau_sampling, tau_actual_size, perturbations_sources, perhaps_print_variables, ppt->error_message);\n=>evolver_ndf15(L:187) :error in initialize_jacobian(&jac,neq,error_message);\n=>initialize_jacobian(L:1574) :could not allocate jac->Ci with size -8" + ] + } + ], + "source": [ + "# Near degenerate case\n", + "masses = [0.0584, 0.0312, 0.03]\n", + "#masses = [50, 5, 0.5]\n", + "degs = [1, 1, 1]\n", + "Nbins = 5\n", + "\n", + "standard_settings = {\n", + " 'H0': 70.0,\n", + " 'Omega_cdm': 0.26,\n", + " # etc...\n", + " \n", + " 'N_ncdm': len(masses),\n", + " 'm_ncdm': ', '.join([str(i) for i in masses]),\n", + " 'deg_ncdm': ', '.join([str(i) for i in degs]),\n", + "\n", + " 'Quadrature strategy': '3, 3, 3',\n", + " 'Maximum_q': '12, 12, 12',\n", + " 'Number of momentum bins': f'{Nbins}, {Nbins}, {Nbins}',\n", + " \n", + "# 'output_ncdm_binning': 0,\n", + " 'output': 'tCl',\n", + "# 'ncdm_fluid_approximation': 0,\n", + "}\n", + "\n", + "reference = Class()\n", + "reference.set(standard_settings | {'Number of momentum bins': '10000, 10000, 10000', 'Quadrature strategy': '2, 2, 2'})\n", + "tic = perf_counter()\n", + "reference.compute()\n", + "toc = perf_counter()\n", + "print(f\"Computed reference in {toc - tic:.5} s\")\n", + "\n", + "bg_ref = reference.get_background()\n", + "a_ref = 1/(1 + bg_ref['z'])" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "87ca073a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Computed 3 separate species in 0.068621 s\n", + "Computed 3 collective species in 0.04097 s\n" + ] + } + ], + "source": [ + "separate = Class()\n", + "separate.set(standard_settings)\n", + "tic = perf_counter()\n", + "separate.compute()\n", + "toc = perf_counter()\n", + "print(f\"Computed {len(masses)} separate species in {toc - tic:.5} s\")\n", + "\n", + "collective = Class()\n", + "collective.set(standard_settings | {'collective_ncdm': 1})\n", + "tic = perf_counter()\n", + "collective.compute()\n", + "toc = perf_counter()\n", + "print(f\"Computed {len(masses)} collective species in {toc - tic:.5} s\")" + ] + }, + { + "cell_type": "markdown", + "id": "cb2ee85f", + "metadata": {}, + "source": [ + "Pretty OK speed-up, about factor 2 in background." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "7359e129", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAArEElEQVR4nO3deXxU1f3/8ddnJgkhCTsJOyQIsskWwqIWRanFVgQrimht5Vv9WvypX61f9evS1q2t/bb9aX9Vq0VF1FatYom4VUGl4sqigAjIJkhQtgCBAAlZzu+PGWLIOhNmMneS99PHPJh75txzPzeZycdz75lzzDmHiIiI1/hiHYCIiEhNlKBERMSTlKBERMSTlKBERMSTlKBERMSTlKBERMSTlKBE4oSZbTaz79by2lgz+6KxYxKJJiUokSgys7+Z2Tdmtt/M1pnZFVVeH29ma83skJm9Y2a9GnIc59wi51y/yEQt4g1KUCLRdS+Q6ZxrDUwCfm1mIwDMrCPwT+CXQHtgKfCPWAUq4jVKUCJR5Jz73DlXfHQz+DghuH0+8Llz7gXnXBFwJzDUzPrX0eRIM1ttZnvN7AkzSwYws3Fmlne0UvBy4I1mttLMCszsH5XqdjSzV8xsn5ntMbNFZqa/BeI5elOKRJmZ/cXMDgFrgW+A14IvDQJWHK3nnDsIbAyW1+ZHwAQCSe5E4Bd11J0KnA1kAUOA6cHy/wbygHSgE3AbgcQp4ilKUCJR5pz7P0ArYCyBS3pHe1RpQEGV6gXBurV50Dm31Tm3B/gNcHEddf/snPs6WPdlYFiwvAToAvRyzpUE718pQYnnKEGJNALnXJlz7j2gO3BVsLgQaF2lamvgQB1Nba30fAvQtY662ys9P0QgIQL8AdgAvGlmm8zslnrCF4kJJSiRxpXAt/egPgeGHn3BzFKDr31ex/49Kj3vCXwdbgDOuQPOuf92zvUmMHDjBjMbH247ItGmBCUSJWaWYWbTzCzNzPxmNoHAJbm3glXmAieZ2ZTgAIZfASudc2vraPZqM+tuZu2B22nAqD8zm2hmfczMCFxSLAPKw21HJNqUoESixxG4nJcH7AX+CFzvnJsH4JzbBUwhcC9pLzAamFZPm88AbwKbCAyo+HUD4uoLLCBwifFD4C/OuXca0I5IVJnujYqIiBepByUiIp6kBCUiIp6kBCUiIp6kBCUiIp6UEOsAGkPHjh1dZmZmrMMQEZEaLFu2bLdzLr1qebNIUJmZmSxdujTWYYiISA3MbEtN5brEJyIinqQEJSIinqQEJSIinqQEJSIinqQEJSIinqQEJSIinqQEJSIintQsvgd1vPbl5uKKj9Rf0UJs0OqvaCHUCaWdYMXItRVSXKE21chxhR5YBJtq+ufY6O/niP4cInM88/vxtWyJtWyJLyUFX0oqvpSW+NLSQjt3qZESVAjW/v0RfIWH6q4UwVVLLJQlUByh5Z1Q4nIhNhbCSYZ0vFA5F9qfmVj87OttKLR6IR0vVCG1ZfXWC/nPaaixh/JziPQxI9RWKO9nM8PKyrGiI9iRkmNfS0wkIT2dhIwMEjIySMrKokWfPrQ4sS8tevfGEhMbGn2zoAQVgrnXDqWwpLDOOqGsq+VC+KSGUie0KiEcKx5jbsTzCrWtkNppqr+PRozZa78LgJLyEorLijlcepiiUj/FJZB0pJwWJdCyGNKPJNKjuJRuh/fQregwnZdupPVLc/Bvz8dSUkgZPpyUUaNIGzeOFif2VW+rimaxYGFOTo7TVEciEm3OOUrLSyksKST/cD75RfnkH85nx6EdbNm/hU0Fm9i4byNHDhQwdl8nTtnZlr6bDpO0ZjNJvXrRasIE2v7wPJKa2dyhZrbMOZdTtVw9KBGRCDEzEv2JtPO3o11yO/rQp1od5xyb92/mkx2fsHj7Yu7LW0Tq/hQu3tmKUR8uJH/mTFLGjKbdRdNoNf7MZn0ZUD0oEZEYKikrYfH2xczbOI/5W+aTfTiDn6zrTKd3V+NLbknbaRfRbupUEjp2jHWoUVNbD0oJSkTEI/YU7SF3Qy5/W/03Eo84rts9nBMXrKN0y1e0/sH3aXfppbQcPDjWYUacEpQSlIjEieKyYnLX5/L4qscpLi3ias5g1Hu7OPz2v0kefBLtL/0xrSd8D0tKinWoEaEEpQQlInGmpLyEVza+wqOfPcr+I/u5PH0S3/2kjMMvzoPEBNpNvYi2F00lMSMj1qEeFyUoJSgRiVOl5aW8/uXrzFw5k/zD+Vx6wlTO+7IjRc//k+J162l1xjjSTj+d1LFjSezUKdbhhk0JSglKROJcWXkZ87fM568r/8pX+79ifI8zmVI8iB4fb+HQovcp2bqVxF49aTloEMmDTiIpK4vE7t1I6tYNX2pqrMOvlRKUEpSINBHlrpyl25eSuyGX+Vvmk+hLZEyX0XzH9aHvtnLafbmH0rXrOLJlC2X5+QD4UlPxt2mDv21b/G3b4EtNxRKTsKQkrEWLwL9JiZg/AXyG+XyAgc8HRmDbfIGpnyq97m/fnrY/PO+4zkcJSglKRJqgQyWHWLpjKYvyFrFs5zI27duEw9EltQudUzvTIyGd7oUtaFvoaFVkpBY5kg+WkFBUgq+0HF9JGf6ScnylZdiRUqysHHBQ7gJTPbnywMwjzlU8rDyYN8rLSejWld5/uP+4zkFf1BURaYJSElM4rftpnNb9NCAwAnDD3g1sLdzK9sLtfH3wazYk7uVA6gEOHDnA/iP7OVR6iNLyUkrKSwL/lpVQ6kobdPy+7Ur4ZyRPqBIlKBGRJqSFvwWDOg5iUMdBYe13dJqmcspxzuFwFXMSHn1e8Z/7dttn0Vu1SQlKREQqpmnyEi1YKCIinqQEJSIinqQEJSIinqQEJSIinhSXgyTM7DzgHKA18Lhz7s3YRiQiIpHW6D0oM5tlZjvNbFWV8rPN7Asz22Bmt9TVhnMu1zn3n8AM4KJoxisiIrERix7UbOBB4KmjBWbmBx4CzgLygCVmNg/wA/dW2f+nzrmdwee/CO4nIiJNTKMnKOfcu2aWWaV4FLDBObcJwMyeAyY75+4FJlZtw8wM+B3wunPuk5qOY2ZXAlcC9OzZM3InICIijcIrgyS6AVsrbecFy2pzLfBd4AIzm1FTBefcTOdcjnMuJz09PXKRiohIo4jLQRLOuT8Df451HCIiEj1e6UFtA3pU2u4eLBMRkWbKKwlqCdDXzLLMLAmYBsyLcUwiIhJDsRhm/izwIdDPzPLM7HLnXClwDfAGsAZ43jn3eWPHJiIi3hGLUXwX11L+GvBaI4cjIiIe5ZVLfCIiIsdQghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9SghIREU9KaMhOZvbfQE5w/8XOuT9ENCoREWn2GtqD8jnnLnbOXQh0imRAIiIi0MAeFDDIzC4J7t/FzH7gnHstgnGJiEgz19AE9Q6QGHz+JtAxMuGIiIgENChBOeeejHQgIiIilTXoHpSZXWdmjwaf/zKyIYmIiDR8kMQJwNbg81YRikVERKRCQxOUA1qa2UlA1wjGEzIzSzWzpWY2MRbHFxGR6Gpogvq/gAF3AbeGs6OZzTKznWa2qkr52Wb2hZltMLNbQmjqf4Dnwzm2iIjEj4YmqGuBB4EjwL1h7jsbOLtygZn5gYeA7wMDgYvNbKCZDTazV6o8MszsLGA1sLOB8YuIiMc1dJh5a2ASgeR0aTg7OufeNbPMKsWjgA3OuU0AZvYcMNk5dy9Q7RKemY0DUgkks8Nm9ppzrrxKnSuBKwF69uwZTogiIuIBDU1QC4HezrmVZrY+AnF049tBFwB5wOjaKjvnbgcws+nA7qrJKVhnJjATICcnx0UgRhERaUQN/R7Us5WePxq5cMKOY3asji0iItHlldnMtwE9Km13D5aJiEgzddwJysx+GrwndDyWAH3NLMvMkoBpwLzjjU1EROJXJHpQrwAtQq1sZs8CHwL9zCzPzC53zpUC1wBvAGuA551zn0cgNhERiVMNXQ+q8uzlY5xzIfd2nHMX11L+GqAZ0UVEBGhAgjKzcwh8TwkCX9a9BF2OExGRCGtID6ojcDj4bznw24hGJCIiQgMSlHPuSTPbS+ALug44AKyMdGAiItK8NXSQxETn3FTn3EVUmbZIREQkEho6k0RLMzs6f1BqpIIRERE5qqEJ6k4CE8YC3B2ZUERERL7V0ATVzzl3E4CZTQI2RC4kERERDTMXERGPOp5h5ulomLmIiERJyAnKzNoDPwcygFXAPOfc3mgFJiIizVs4w8yfI/Cdp5eBFOA9MxsVlahERKTZC+cSX7pz7vfB56+Y2T+AZ4AxkQ9LRESau3B6UHvMbPDRjeDy7CmRD0lERCS8HtTVwBwzWwR8BgwCNkYlKhERafZC7kE559YC2cA7BAZKLAdqXDpDRETkeIUziu9M4EfAPgKj+FYSmCxWREQk4sK5xDcLuB5IBIYA5xG4zNcn4lGJiEizF06C2uKcyw0+fyEKsYiIiFSo9x6UmT1lZtcDH5nZDdEPSUREJLRBErMJzLnXCfixmW0xs3lmdo+ZXRjV6EREpNmq9xKfc+5tM1sI3OKcm25mCcAAYCgwEl3uE5Embv++fNa8OQv7Zjllye3pevp/0Kt/dqzDavJCGmbunCsHJgaflzrnPnPO/c05d3NUoxMRibGlrz5K2Z+G0WPNTMCRtmcVXZ49iyW5D8U6tCYvnEESK83sDuCeYMISEWmyysvK+Pix/2L41/9geb/rGDn1VromBP5kLsl9iKGf/pLVnU9g4JizYxxp0xVOgmoPnA5cZWYfE/ge1ErnnC7xiUiTUlx0iM8eupT+BxazeeJzjBn53WNeH3ne1Xz0zUq6vHkjJdlnkJjUIkaRNm3hzCQx1Tk3AOgF3EVgFV3NZi4iTUrB3t1suG8CXQpXUfijV+lfJTkdddKPfkfr8n2s+NesRo6w+QhlmPnJFlw+F8A5V+yc+8Q59+TRZd9FRJqCndu+ZM8D40kqP0yLn71Fj75Da62b1roda7tfSNsVjzZihM1LKD2onwDLzOw5M5tuZp2jHZSISGNbu3g+7tHxHEhKp+t1C+jYuUe9+2RNuJo+ZRvZsvaTRoiw+ak3QTnnrnLOZQN3Au2A2Wb2oZn91sxOMzN/tIMUEYmWgwf28eGj13HCqxfxZafvMeCGV0lt1TakfTv36MPaxIF8/f4z0Q2ymQp5kERwNvO1wP1m1hI4A7gQuA/IiU54IiKRV1pyhE2ffUD+khfo/808ulsaa854lDHjpoTd1r7M75O+KTfyQUpYo/gqO8M59xrwmplNjmRAXvTho9djJQdrr+DqntTdQpn0vZ42IjJxfES+HVBPHF75WdTbRt0aJ876j2GN8r6IwO/sON8XoR2nvmPUewgSyg7T+sh2Msp20psSSloMYmP2bQyZ8B/0aJFcfwM16DRsAr3W3ce+3dtp27HmOyD79+Wz/rGf0uHQJvaMvpnsCT9u0LFCVVZaSv6Orez95ksO79tB6aF9lB3eT3lRARQfwEqLsfISrOxI4N/yUnzlJfhcCebKAQfOERh84AAXfC+6it/T0XoHU3ow6vpno3IeYScoMzsHuDg4bsIIrAn1UoTj8hRf0V78dSUowH07jqQOEagTwnFcSMepL4wInU897dQYa9V9Ivazra+Jb694H296q/vlCP0O620nEm2E2E69TYSzeHfNQvuM1X2c4oQWFLfrQWGn3mQNGcugEC/l1SVzwEj2WGu2LHmd7O//R4111j3+n3Q4tJntXc8i+4Pr2dA5iz5Dv3PcxwbI35HH5qX/4kjep6Ts/YJORZvp6PLJsHLaOx/7rDWHLYXD/jSK/amUJKRR7m+B8yVS7kvEJabh/Ik4XxL4E4K/Kwu+N4L/mmH4cEf/7Jsv8K7w+fCnpUfkPGrSkB5UR+Bw8N9y4LcRjciDRl/7ZKxDEBGPMp+PzWnDKd24CKieoDat+pjs/W/z1bQFjBmQw5L7t9Hytdth6L8bfMztWzfw5dtPkPHVq5xQ9iWOtmxNGcihjoPZ1mUa+7v0oV3nTNpndKdjQkMvlMVe2JE75540s73ApQT+R/MAgS/tiog0SyWdh9Nu8+s1vrbrnb9wIGUkQwcEbtV3nXwn3Z4aw+Y1S8kcEN7t+63rV7Dj5XsYXrCAgwmZ7MqajH/EufTql01H3/H3Ur2moal1onNuKoCZPUwTv8QnIlKXVr1H0Wv9g5QcKT5mVonysjJOyF/I5pG/qCjr1nsAq1oMo3DhYyEnqCPFRSz7+y8ZseVx9qSMYP25L9I/ZzwnRPxMvKWhCaqlmfUMPk+NVDAiIvGo56AxJL1RypdffMoJg8dUlG/87AN6uwISTz7vmPoH+0yi2+rHQmp7765v+GbmBZxQ8jVrTv8LQ8+cFsnQPa2hfcI7gWuDj7sjFo2ISBxq1aY9ef6u7Fn/8THlu1e8zvqkAbRpf+xAgl5jJtPDfU3ehlV1trtn5zb2P/xdjHL8Vy1qVskJGt6D6nd0miMzm0RgXr5GY2Y+4B6gNbDUOadRDCISU7tb9qZs5xfHlLXYuZJ97atPl9S5Rx82+3qy45PX6d7npBrbO3zwAHv+ei7FCe044fp/kZySFpW4vSzsHlRwmPklZvaDo8/D3H+Wme00s1VVys82sy/MbIOZ3VJPM5OB7kAJkBfO8UVEoqG4bR9aFhz7/+oZh9bj7zqkxvo72w7Ft21pre2tnHU1LcoP0/Pqec0yOUHDLvEdHWaeTmAJjnCHmc8GjllAJThd0kPA94GBBL5nNdDMBpvZK1UeGUA/4APn3A3AVQ04BxGRiErs1J+ORZsrtg8U7KG7207HPiNq3qFbDp0OfFbjS5+//yrZu1+h+LxHadWmfRSijQ8hX+Izs/bAz4EMYBUwzzm3N9wDOufeNbPMKsWjgA3OuU3BYz0HTHbO3UtwJd8qseQBR4KbZbXEeyVwJUDPnj1rqiIiEjFtew2my7KdFB0qJDkljW1rl9LbJdC977Aa62cM/A49P7uDgr27adOuY0V5eVkZyW//gmWdpjAmQl/mjVfh9KCeI/Cdp5eBFOA9M4vUelDdgK2VtvOCZbX5JzDBzB4A3q2pgnNupnMuxzmXk54evW86i4gAdOsTuJS3bUPga6EFmz/hq4RetS5m2L3PEI44P9+s//SY8pULX6BL6Tb6T70nugHHgXAGSaQ7534ffP6Kmf0DeAYYU8c+UeGcOwRc3tjHFRGpTcvUVnzty2Dvls9gyCnY9s/Y2+rEWusnJCax2d+N/VtXwaizKsp9S2byWfo5jK5lXr/mJJwe1B4zG3x0I3g5LiVCcWwDKi++0j1YJiISN3YlZ1KyfQ0AbQ+soyyj5hF6R+1JyaS80si/3V9v4aTDn5A+bkZU44wX4fSgrgbmmNki4DNgELAxQnEsAfqaWRaBxDSNMEcHiojE2uE2fUjet57SkiP0LPmSkl7D66xf3O5EUnatqNje+P4cinydyBo4MtqhxoWQe1DB9aCygXcIDJRYTmAm87CY2bPAh0A/M8szs8udc6XANcAbwBrgeefc5+G2LSISS/7OA+lw+Eu2bVxFspXQfUDdt+n97TNpe2R7xXbSlwvIyxiHNcF59RoinFF8ZwI/AvYRGMW3kgasSuCcqzGpHV1fKtz2RES8ok3PwXT99Bs+WfM+LUinc/u6B2ildOxJx7JduPLAWm2Zh1axediPGiPUuBDOJb5ZwPVAIjAEOI/AZb4+EY9KRCQOdes7lAQrJ/mLuXyTciL1DXNo0zmTVCuioGAPhXt30Y39lJ00tlFijQfhJKgtzrnc4PMXohCLiEhcS23Vlm9IZ0jRMj7MrH8OgY5dewOw5+uN7NmyigQ60Klrr2iHGTfqvdBpZk+Z2fXAR2Z2Q/RDEhGJX1+3Cgx2bjfgzHrrtkxtxV5asX/HFo7sXM+uFj3q3ac5CeVO3GwCaz53An5sZlvMbJ6Z3WNmF0Y1OhGRONPl/Hv5qO8N9Bv53ZDq7/e14cj+XSTs28zBNPWeKqs3QTnn3gb+H7DOOTccOAG4HfgC0FhIEZFKumb1Z8yP7gh5JN4hf2vKDuaTdvArXLveUY4uvoR0D8o5V25mE4HfBoeEfxZ8iIjIcShKbIs7tId2JTs51F49qMrCGWy/0szuCK7FJCIiEVCS1BZf0R7auAJattP0RpWFk2zaE5jh4Wsze0n3oEREjl95cjtaHvyalnaE1HadYh2Op4Q8zNw5NxXAzFoQ+P7TYALLZGjIuYhIA7mWbemw4ysA2nTsGuNovKXeBGVmJwMfOeccgHOuGPgk+BARkeNgSal0Kd9JKT5at9PSQJWFconvJ8AyM3vOzKabmS6SiohEiC8pFZ85DlgqPr8/1uF4Sr09KOfcVQBm1p/AkuyzzawNgUlj/wW875yrcVVbERGpmy85FYBikmMcifeENZu5c+5+59zZwHigCLgI+DhawYmINHUJScEE5VOCqiqcufgws+EElti4CNgO9HfOtYlGYCIizUFCsAd1RAmqmlAGSZxIICldAhwgMGpvnHPuSzP7MsrxiYg0aQkt0wAlqJqE0oNaS2DF2wucc1Vnjwh7PSgREflWi2CCcqYBElWFkqDOJ/AF3TfNbAHwPPAv51xJVCOLspKSEvLy8igqKop1KHIckpOT6d69O4mJibEORaRBfAktACg3vYerCmUUXy6Qa2apwGTgSuAxM3sNaB3d8KInLy+PVq1akZmZiZnFOhxpAOcc+fn55OXlkZWVFetwRBokIfg/V2U+JaiqwhnFd9A594xz7lygP/AhgWXf41JRUREdOnRQcopjZkaHDh3UC5a45k9IinUIntWgiV+dc3udczOdc/WvyOVhSk7xT79DiXcJwQRluqVfjWYmFxGJIX9isAfllKCqUoKKI7m5uaxevTrWYYhIBB1NUOpBVacEFUeUoESansSKHlR5bAPxICWoGDvvvPMYMWIEgwYNYubMmQCkpaVVvD5nzhymT5/OBx98wLx587jpppsYNmwYGzduZPny5YwZM4YhQ4bwwx/+kL1798bqNESkgRIqelBSVVhTHTVVZeWO3YXFEW+3Y1oL/L6633azZs2iffv2HD58mJEjRzJlypQa651yyilMmjSJiRMncsEFFwAwZMgQHnjgAU4//XR+9atfcdddd/GnP/0p0qchIlGUkHB0eLku8VWlBAXsLixm9G/fini7H982nk6t656+5M9//jNz584FYOvWraxfvz6ktgsKCti3bx+nn346AJdddhkXXqgFjkXijfkCF7JMl/iqUYIi0NP5+LbxUWm3LgsXLmTBggV8+OGHpKSkMG7cOIqKio4ZOq3v+Ig0D850x6UqJSjA77N6ezrRUFBQQLt27UhJSWHt2rV89NFHAHTq1Ik1a9bQr18/5s6dS6tWrQBo1aoVBw4cAKBNmza0a9eORYsWMXbsWJ5++umK3pSIxB+nu1DVKGXH0Nlnn01paSkDBgzglltuYcyYMQD87ne/Y+LEiZxyyil06dKlov60adP4wx/+wPDhw9m4cSNPPvkkN910E0OGDGH58uX86le/itWpiMjxUg+qGvWgYqhFixa8/vrrNb52dCBEZaeeemq1YeZHe10iEt/KNZt5NUrZIiKeoEt8VSlBiYh4geaVrEYJSkTEAzRIojolKBERD9BcfNUpQYmIeIFmM69GCUpExBOUoKpSgoqxHTt2cMkll9C7d29GjBjBySefXDH1UWOZPn06c+bMadRjiojUJy4TlJn1NLNcM5tlZrfEOp6Gcs5x3nnncdppp7Fp0yaWLVvGc889R15eXqxDa3SlpaV1bos0dboHVV2jJ6hgUtlpZquqlJ9tZl+Y2YYQks5gYI5z7qfA8KgFG2Vvv/02SUlJzJgxo6KsV69eXHvttWzevJmxY8eSnZ1NdnY2H3zwARCYv2/ixIkV9a+55hpmz54NwC233MLAgQMZMmQIN954IwAvvPACJ510EkOHDuW0004LObbCwkLGjx9PdnY2gwcP5qWXXqp47Z577qFfv3585zvf4eKLL+aPf/wjAOPGjWPp0qUA7N69m8zMTADKysq46aabGDlyJEOGDOGvf/1rxbmMHTuWSZMmMXDgwGrbte0n0hRpstjqYjGTxGzgQeCpowVm5gceAs4C8oAlZjYP8AP3Vtn/p8BHwBwz+ynw9HFHVF4GhTuPu5lq0jLAV/u3wz///HOys7NrfC0jI4P58+eTnJzM+vXrufjiiyv++NckPz+fuXPnsnbtWsyMffv2AXD33Xfzxhtv0K1bt4qyUCQnJzN37lxat27N7t27GTNmDJMmTWLp0qW8+OKLrFixgpKSErKzsxkxYkSdbT3++OO0adOGJUuWUFxczKmnnsr3vvc9AD755BNWrVpFVlYWCxcuPGZ75syZNe6XlZUV8nmIxA/1oKpq9ATlnHvXzDKrFI8CNjjnNgGY2XPAZOfcvcDEKnUxsxuBO4JtzQGeqKHOlcCVAD179qw7qMKdcF//8E+mPjeshdZd6q8XdPXVV/Pee++RlJTEggULuOaaa1i+fDl+v59169bVuW+bNm1ITk7m8ssvZ+LEiRW9rFNPPZXp06czdepUzj///JBjcc5x22238e677+Lz+di2bRs7duzg/fffZ/LkySQnJ5OcnMy5555bb1tvvvkmK1eurLjPVVBQwPr160lKSmLUqFHHJJzK27XtpwQlTZJG8VXjlbn4ugFbK23nAaPrqP8v4E4zuwTYXFMF59xMYCZATk5O3b/5tIxAMom0tIw6Xx40aBAvvvhixfZDDz3E7t27ycnJ4f7776dTp06sWLGC8vJykpMDs60nJCRQXv7tpYCjy3EkJCSwePFi3nrrLebMmcODDz7I22+/zSOPPMLHH3/Mq6++yogRI1i2bBk33ngjn376KV27duW1116rMba///3v7Nq1i2XLlpGYmEhmZma9S39Ujq1yXeccDzzwABMmTDim/sKFC0lNTT2mrPJ2bfuJNEX6mm51cTlIwjm3yjl3gXNuhnPuxuNu0OcP9HQi/ajj8h7AmWeeSVFREQ8//HBF2aFDh4BAb6FLly74fD6efvppysrKgMA9qtWrV1NcXMy+fft4663AQouFhYUUFBTwgx/8gPvvv58VK1YAsHHjRkaPHs3dd99Neno6W7du5YknnmD58uW1Jqejx8/IyCAxMZF33nmHLVu2AIEe2csvv0xRURGFhYW88sorFftkZmaybNkygGNGBU6YMIGHH36YkpISANatW8fBgwfr/bU0dD+ReKT+U3Ve6UFtA3pU2u4eLGvSzIzc3Fx+/vOf8/vf/5709HRSU1P53//9X7Kzs5kyZQpPPfUUZ599dkXPokePHkydOpWTTjqJrKwshg8PjBE5cOAAkydPpqioCOcc9913HwA33XQT69evxznH+PHjGTp0aI2x/OxnP+P666+vOMbLL7/Mueeey+DBg8nJyaF//8Al0JEjRzJp0iSGDBlCp06dGDx4MG3atAHgxhtvZOrUqcycOZNzzjmnou0rrriCzZs3k52djXOO9PR0cnNz6/35NHQ/EWkazMXgumfwHtQrzrmTgtsJwDpgPIHEtAS4xDn3eSSOl5OT46oOMFizZg0DBgyIRPPNTmFhIWlpaRw6dIjTTjuNmTNn1jrYozHodylx7842rEgeydBbFsQ6kpgws2XOuZyq5Y3egzKzZ4FxQEczyyMw2OFxM7sGeIPAyL1ZkUpOEnlXXnklq1evpqioiMsuuyymyUmkydBs5tXEYhTfxbWUvwbUflNEPOOZZ56JdQgi0gzE5SAJERFp+pSgRETEk5SgRETEk5SgREQ8QYMkqlKCipH8/HyGDRvGsGHD6Ny5M926davYPnLkSNSOe8UVV7B69eo66+zatYvRo0czfPhwFi1aFLVYRETq4pUv6jY7HTp0YPny5QDceeedpKWlVcxADoHlJhISIv/reeyxx+qt89ZbbzF48OCQ6oqIRIt6UB4yffp0ZsyYwejRo7n55ptZvHgxJ598MsOHD+eUU07hiy++AGD27NlMnjyZcePG0bdvX+666y4AHnnkkYpeWFZWFmeccUa1Y1ReEiMtLY3bb7+doUOHMmbMGHbs2MHy5cu5+eabeemllxg2bBiHDx9uvB+AiEgl6kEBZeVl5BflR7zdDskd8NczH19VeXl5fPDBB/j9fvbv38+iRYtISEhgwYIF3HbbbRWTyy5evJhVq1aRkpLCyJEjOeecc5gxYwYzZsygpKSEM888kxtuuKHOYx08eJAxY8bwm9/8hptvvplHH32UX/ziF9x9990sXbqUBx98sMHnLiJyvJSggPyifMa/MD7i7b514VtkpNQ9o3lVF154IX5/IKkVFBRw2WWXsX79esysYtJUgLPOOosOHToAcP755/Pee++RkxOYKeS6667jzDPPrHcpjKSkpIplOUaMGMH8+fPDilVEIsdpkEQ1SlAEejpvXfhWVNoNV+XlJn75y19yxhlnMHfuXDZv3sy4ceMqXrMq06Ic3Z49ezZbtmwJqfeTmJhYsZ/f79cy6yLiKUpQgN/nD7un0xgKCgro1q0bQMWy7kfNnz+fPXv20LJlS3Jzc5k1axbLli3jj3/8I4sWLcLn0+1FEYlv+ivmYTfffDO33norw4cPr9a7GTVqFFOmTGHIkCFMmTKFnJwcHnzwQfbs2cMZZ5zBsGHDuOKKK2IUuYjI8YvJchuNrakttzF79mwNYqgknn+XIgDc2YblLccw7H/eiHUkMVHbchvqQYmIeIGW26hG96Di0PTp05k+fXqswxARiSr1oEREPEDDzKtTghIR8QCj6Y8HCJcSlIiIeJISlIiIB+gSX3UaJBEj+fn5jB8fmF5p+/bt+P1+0tPTgcA8e0lJSbEMT0Qk5pSgYiRWy22IiMQL/QX0kOnTp5OcnMynn37KqaeeyrRp07juuusoKiqiZcuWPPHEE/Tr14/Zs2czd+5cCgoK2LZtG5deeil33HEHjzzyCI888ggQmCYpMzOTd955J8ZnJSLSMEpQgCsro3R35JfbSOjYAfN7d7kNEREvU4ICSnfns+H00yPebp9//5vETt5dbkNExMuUoAj0dPr8+99RaTdcjbnchoiIlylBAeb3h93TaQxabkOkOdEw86r0V8zDtNyGiDRn6kF5wJ133llj+cknn8y6desqtn/9619XPO/evTu5ubnH1H/iiSeiEZ6IRNmK5JG0OO2/Yh2G5yhBiYjE2NBbFsQ6BE9SgopDWm5DRJqDZn0PqjmsJtzU6Xco0nQ12wSVnJxMfn6+/sDFMecc+fn5JCcnxzoUEYmCZnuJr3v37uTl5bFr165YhyLHITk5me7du8c6DBGJgmaboBITE8nKyop1GCIiUotme4lPRES8TQlKREQ8SQlKREQ8yZrDKDYz2wVsiXUcdegI7I51EA2k2BtfvMYNij1WvB57L+dcetXCZpGgvM7MljrncmIdR0Mo9sYXr3GDYo+VeI1dl/hERMSTlKBERMSTlKC8YWasAzgOir3xxWvcoNhjJS5j1z0oERHxJPWgRETEk5SgRETEk5SgRETEk5SgRETEk5SgPMrMepvZ42Y2p1LZADN7xMzmmNlVsYyvNrXEXa3My8xsoJk9b2YPm9kFsY4nHGY2NvgeeczMPoh1POEws3FmtigY/7hYxxOOePhs1sbLn08lqCgws1lmttPMVlUpP9vMvjCzDWZ2S11tOOc2Oecur1K2xjk3A5gKnBpHcVcri5ZInAPwfeAB59xVwE+iFmwVEfr5Lwq+R14BnoxmvJVF6OfugEIgGciLVqxVRejnHtXPZm2i9Zn1DOecHhF+AKcB2cCqSmV+YCPQG0gCVgADgcEE/phUfmRU2m9OlbYnAa8Dl8RT3LWVefEcgo+HgD8A78fp++Z5oFU8xQ74gvt1Av4eT7EH94naZ7OR3jNR/3yGfX6xDqCpPoDMKm+ak4E3Km3fCtwaQjs1vmmAV+Mt7sb6AETwHPzAS/H2vgF6Ao82ZtwR/rknNfYfy0jFHqwblc9mI/zcPZegdImv8XQDtlbazguW1cjMOpjZI8BwM7s1WDbOzP5sZn8FXotqtN+KRNzVyhpZuOeQaWYzgacI9KJiKazYgy4HnohaRKEL9+d+fvC9/TTwYJRjq0+4scfis1mb4/7MekWzXfLd65xz+cCMKmULgYWxiCdUtcRdrczLnHObgStjHUdDOefuiHUMDeGc+yfwz1jH0RDx8NmsjZc/n+pBNZ5tQI9K292DZV4Xr3FXFs/noNhjQ7F7gBJU41kC9DWzLDNLAqYB82IcUyjiNe7K4vkcFHtsKHYviPVNsKb4AJ4FvgFKCFz/vTxY/gNgHYERNrfHOs6mEndTOQfFrtibU+yhPDSbuYiIeJIu8YmIiCcpQYmIiCcpQYmIiCcpQYmIiCcpQYmIiCcpQYmIiCcpQYmIiCcpQYmIiCcpQYnEETPLNbNlZva5mcXthLYiodBMEiJxxMzaO+f2mFlLAnOune4Cs1GLNDlabkMkvvyXmf0w+LwH0BdQgpImSQlKJE6Y2Tjgu8DJzrlDZrYQSI5lTCLRpHtQIvGjDbA3mJz6A2NiHZBINClBicSPfwEJZrYG+B3wUYzjEYkqDZIQERFPUg9KREQ8SQlKREQ8SQlKREQ8SQlKREQ8SQlKREQ8SQlKREQ8SQlKREQ86f8DgOTPpsliiwAAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "def get_density_error(Nbins, strategy=1):\n", + " collective = Class()\n", + " collective.set(standard_settings | {'collective_ncdm': 1, 'Number of momentum bins': f'{Nbins}, {Nbins}, {Nbins}', 'Quadrature strategy': f'{strategy}, {strategy}, {strategy}'})\n", + " collective.compute()\n", + " bg_col = collective.get_background()\n", + " a_col = 1/(1 + bg_col['z'])\n", + " rho_tot_ref = interp_onto(a_ref, np.sum(np.vstack([bg_ref[f'(.)rho_ncdm[{idx}]'] for idx in range(len(masses))]), axis=0), a_col)\n", + " error = (bg_col['(.)rho_ncdm[0]'] - rho_tot_ref)/rho_tot_ref\n", + " return a_col, error\n", + "\n", + "bins = 10\n", + "\n", + "fig, ax = plt.subplots(1, 1, figsize=(6, 4))\n", + "ax.set_title(f'{Nbins} bins')\n", + "ax.set(xlabel='a', ylabel=r'$\\Delta /rho_\\mathrm{tot}/rho_\\mathrm{tot,sep}$', xscale='log', yscale='log')\n", + "\n", + "a_auto, error_auto = get_density_error(bins, strategy=0)\n", + "ax.plot(a_auto, np.abs(error_auto), '-', lw=1.2, label='auto')\n", + "\n", + "a_lag, error_lag = get_density_error(bins, strategy=1)\n", + "ax.plot(a_lag, np.abs(error_lag), '-', lw=1.2, label='Gauss-Laguerre')\n", + "\n", + "a_trapz_inf, error_trapz_inf = get_density_error(bins, strategy=2)\n", + "ax.plot(a_trapz_inf, np.abs(error_trapz_inf), '-', lw=1.2, label='Trapz inf')\n", + "\n", + "a_trapz, error_trapz = get_density_error(bins, strategy=3)\n", + "ax.plot(a_trapz, np.abs(error_trapz), '-', lw=1.2, label='Trapz')\n", + "\n", + "ax.legend(loc='lower left')\n", + "fig.tight_layout()\n", + "fig.savefig('density_error.pdf')" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "221a04ef", + "metadata": {}, + "outputs": [], + "source": [ + "def get_ini_rho_error(Nbins, strategy=1, auto_tol=None):\n", + " collective = Class()\n", + " if strategy == 0:\n", + " # +10 bins is approximately required to reduce tol_ncdm_bg by 1 order of mag\n", + " if auto_tol is None:\n", + " collective.set(standard_settings | {'collective_ncdm': 1, 'output_ncdm_binning': 1, 'Quadrature strategy': f'{strategy}, {strategy}, {strategy}'})\n", + " elif auto_tol == 'bin':\n", + " auto_tol = 1e-5*0.1*(Nbins/10)\n", + " collective.set(standard_settings | {'collective_ncdm': 1, 'output_ncdm_binning': 1, 'Quadrature strategy': f'{strategy}, {strategy}, {strategy}', 'tol_ncdm_bg': auto_tol})\n", + " else:\n", + " collective.set(standard_settings | {'collective_ncdm': 1, 'output_ncdm_binning': 1, 'Quadrature strategy': f'{strategy}, {strategy}, {strategy}', 'tol_ncdm_bg': auto_tol})\n", + " else:\n", + " collective.set(standard_settings | {'collective_ncdm': 1, 'Number of momentum bins': f'{Nbins}, {Nbins}, {Nbins}', 'Quadrature strategy': f'{strategy}, {strategy}, {strategy}'})\n", + " bins = Nbins\n", + " collective.compute()\n", + "\n", + " bg_col = collective.get_background()\n", + " if strategy == 0:\n", + " # Get amount of bins\n", + " max_bin_idx = np.argmax([bg_col[f'q_ncdm[0][{idx}]'][0] for idx in range(1000) if f'q_ncdm[0][{idx}]' in bg_col])\n", + " print(f\"Automatic quadrature has {max_bin_idx} bins at a tol_ncdm_bg={auto_tol}.\")\n", + " bins = max_bin_idx\n", + " a_col = 1/(1 + bg_col['z'])\n", + " rho_tot_ref = interp_onto(a_ref, np.sum(np.vstack([bg_ref[f'(.)rho_ncdm[{idx}]'] for idx in range(len(masses))]), axis=0), a_col)\n", + "\n", + " error = (bg_col['(.)rho_ncdm[0]'] - rho_tot_ref)/rho_tot_ref\n", + " return bins, error[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "1ebe6456", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Automatic quadrature has 9 bins at a tol_ncdm_bg=None.\n", + "Automatic quadrature has 5 bins at a tol_ncdm_bg=0.01.\n" + ] + } + ], + "source": [ + "binlist = [5, 10, 50, 90, 250, 500, 1000]\n", + "laguerre_cut = 110\n", + "\n", + "results = {}\n", + "strategies = ['auto', 'auto_tol', 'laguerre', 'trapz_inf', 'trapz']\n", + "auto_tol = 1e-2\n", + "for kw in strategies:\n", + " results[kw] = {\n", + " 'bins': [],\n", + " 'error': []\n", + " }\n", + "\n", + "auto_bins, auto_error = get_ini_rho_error(bin, strategy=0)\n", + "results['auto']['bins'].append(auto_bins)\n", + "results['auto']['error'].append(auto_error)\n", + "\n", + "auto_tol_bins, auto_tol_error = get_ini_rho_error(bin, strategy=0, auto_tol=1e-2)\n", + "results['auto_tol']['bins'].append(auto_tol_bins)\n", + "results['auto_tol']['error'].append(auto_tol_error)\n", + "\n", + "for bin in binlist:\n", + " if bin < laguerre_cut:\n", + " laguerre_bins, laguerre_error = get_ini_rho_error(bin, strategy=1)\n", + " results['laguerre']['bins'].append(laguerre_bins)\n", + " results['laguerre']['error'].append(laguerre_error)\n", + "\n", + " trapz_inf_bins, trapz_inf_error = get_ini_rho_error(bin, strategy=2)\n", + " results['trapz_inf']['bins'].append(trapz_inf_bins)\n", + " results['trapz_inf']['error'].append(trapz_inf_error)\n", + "\n", + " trapz_bins, trapz_error = get_ini_rho_error(bin, strategy=3)\n", + " results['trapz']['bins'].append(trapz_bins)\n", + " results['trapz']['error'].append(trapz_error)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "6102f969", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAADkCAYAAABdTMmsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAABN5UlEQVR4nO3dd3xUZdbA8d+ZyaQXIDSpCUUQREGpglIUBAnqrhUVwaUvKoqK2FFRsazrgq4IiBTbi66iRJFiQVCkKNhABSEhoRMghfTkvH/cCQ4hZZJMMpPk+fKZT+b2M8kwZ557n3seUVUMwzAMoyJs3g7AMAzDqP5MMjEMwzAqzCQTwzAMo8JMMjEMwzAqzCQTwzAMo8JMMjEMwzAqzCQTHyUiK0RkZAnL54jII27u6ysRGePmuv1EJNHdOMtDRB4UkfmVeYyaSERURE6KyFPejsXwDhG5TETSRCRfRC7zdjyuTDKpQiIS5+4bQFWHqOoi53ajRGR9oeUTVPXJyoizsqnq06o6BkBEopwfkn7ejquaOF9VHyqYEJHOIvK9iKQ7f3YubkMRqSciHzoTUryI3OSyrJ/zAyrN5THSZXmUiHwqIsdF5KCIvFzU30xEbnX+Pd368uKy7y+dr+G3kv6PiEiAiCwQkRRnHFNclnUQkS3OGI+LyBoR6eCyvL/zOMkiEldovw1F5B0R2e9c/o2I9CjDayg2rmLWv9u5XopzuwCXZV+KyBHnsh9F5KqCZaq6RlVDgb3uxlZVTDIxao1iPvzKlMR8KemJiD/wEfAmUBdYBHzknF+UV4BsoBFwM/CqiHR0Wb5fVUNdHotclv0XOAycBXQG+gL/LBRPXeBB4NcyvpR3gK1AJPAQ8L6INChm3elAW6Al0B+YKiKDC+IHrgXqAfWBj4F3XbY9CSwA7itiv6HAZuBC5/aLgE9EJNTN11BSXKcRkcuBacClzvVbAY+7rDIZOEtVw4FxwJsicpabcXiPqppHFT2AOOAy5/NRwHrgBeA4sAcY4rLuV8AY4BwgE8gD0oATzuULgRnO53WBWOCIc1+xQLPC+yompiDnvo4D27H+oyW6LG8C/M+57z3AnS7LpgNLgcVAKtaHSFeX5fcD+5zLfgcuddnuTefzvYA6X1sa1ofUMaCTy34aAulAg2Jewz+AHc7XsBJo6bJMgUnATmf8/YBEZ2wHgSVAAPAS1ofRfufzAOf2Z6zvxt95IdaH7wrna/oGaOzc73HgN6BLOd4/CrRxmR7k/P2Ky7y9wOAitg3BSiRnu8xbAsx0fZ0lHHsHcIXL9PPAa4XWmYOVYIp9vxWx37OBLCDMZd46YEIx6+8HBrlMPwm8W8R6fs6/e3oRyy4D4tyILQW40M3X4VZczmVvA0+7TF8KHCxm3e5Y//+7F5ofh/OzxFcepmXiXT2wPmTrA88Br4uIuK6gqjuACcAGtb4t1iliPzbgDaxvOS2ADOBlN2N4DGjtfFwOuJ7asAHLgR+Bplhv+ruc36wKXIn17a8O1jfBl53btgNuB7qpaphz33FFHP8S5886zte31rm/W1zWGQ58rqpHCm/sPAXwIPB3oAHWB9E7hVa7Gut3XXDKozHWt8+WWN/8HgJ6Yn3jPh/rP/DDLtsXXt8d1zv3UR/rw3ID8INz+n3gRZfXECsiJ4p5xJZwjI7AT+r8dHH6yTm/sLOBXFX9w2Xej4XWbSgih0Rkj4j8W0RCXJa9BNwoIsEi0hQYAnzm8hq6A12xEkpZdAR2q2pqCXEVHKMuVsvox5LWFZETWB/As4GnyxhPwT46A/7ALjfWdSsuFx2LWLeRiES67DNWRDKBjVjJeUsZwvcKk0y8K15V56lqHlaz+iysUxBloqpJqvo/VU13/qd8CusbvjuuB55S1WOqmgDMclnWDas18ISqZqvqbmAecKPLOutV9VPna1iC9WEMVksqAOggIg5VjVPVP92MaREw3CWxjnDuuygTgGdUdYeq5mJ9eHQWkZYu6zzjfH0Zzul84DFVzXLOuxl4QlUPOxPW485jUsz67vhQVb9X1UzgQyBTVRc7f0//B3QpWFFVY1S1TjGPmBKOEQokF5qXDIQVs25KCev+hpVMzwIGYJ3uedFl3a+xPgRTsFpqW4BlACJix2qJ3a6q+SXE64nXULC82HWdX7gisL7MbC1jPIhIONb77XFVLRxbUdyKq9D6hdfFdX3n3z0MuAJYVY7fa5UzycS7DhY8UdV051N3z9Ge4vy2+JrzomoK1n/8Os7/5KVpAiS4TMe7PG8JNHH9pozVCnBNeAddnqcDgSLip6q7gLuwTmkdFpF3RaSJO69HVTc699VPRNoDbbBaPUVpCfzHJb5jgGC1pAokFNrmiPNDvkATTn/d8c55xa3vjkMuzzOKmC7z37kIaUB4oXnhWKcVy7Suqh5U1e2qmq+qe4CpwDVwqoX6GfAB1umy+linVp917uefWC2k76rgNRQsL3FdVT2J1UpaLCIN3Q1GRIKwWuPfqeozbm7mdlwu6xdel8Lrq2qOqq4ABonIlW7G4jUmmVQPpZV2vgdoB/RQ66JdwakjKX6TUw4AzV2mW7g8TwD2FPqmHKaqV7gVtOrbqtoH6wNf+evD57TVitl8EdaprhHA+yV8mCcA4wvFGKSq35ZwjMLT+50xFmjhnFdajB4hVjfwtGIeK0rY9FfgvEKnRs+j6AvgfwB+ItLWZd75xawL1msu+Hyoh/U7ednZOkvCOq1a8D64FPibs3fSQeAi4F8i4s6p1l+BViLi+i2+yLhU9TjW+/X80tZ1sgHBnP7FoljOHlXLsFpe493Zppxx/VrEuoecv9ei+GGdhvZpJplUD4eAZiX00gnD+rZ7QkTqYV0HcddS4AERqSsizYA7XJZtAlJF5H4RCRIRu4icKyLdStupiLQTkQHO/6CZzviKaqofcc5vVWj+m8DfsBLK4hIONccZf0fncSNE5LrS4ivkHeBhEWkgIvWBR53HL5ZY3V/7lfE4RVKrG3hoMY8hJWz6FdbpxDudXVNvd87/oohjnMRqWTwhIiEi0hu4CufpQ7G6zbYUS3NgJlZPMVT1KFbnhYki4icidbCurf3k3P0orI4inZ2PLVinCh9y7nuUFOqK6xLXH8A24DERCRSRv2ElxP8V85oXY/2t6jpbrWOxOjwgIgNFpIvzfRqOdZruOFbnAUTEJiKBgMOalMCC/1Mi4sC6lpUBjCx8Wkn+6sIeVda4ill3tFhdmetgXVsreA3tRWSI8/+bQ0RuwfpyuLaYffkMk0yqhy+wvs0cFJGjRSx/CatX1lHgO1wujLrhcazTOnuAVbhcm3Ce34/B+oDY49z/fKzz0aUJwPpAOop1Kqwh8EDhlZyn954CvnGequrpnJ+AdcFasS6qF0lVP8Rq8bzrPMX3C9bF4bKYgfUB+BPws/O4M4pb2flhm+pc12tUNRurc8GtwAmsXm1XO+cX3Bzq2rL5J9b75DBWAp2oqgXfnrsA32J1n/0W67Xd6bLt34HBWMl/F5AD3O2M44TzNNlBVT2I1WssxeV6Q3OsHm3FuRHr4v1xrPfMtQWdLUTkZhFx/Yb/GPAn1nt2LfC8qha83+s4X1eyc53WWD3bClq1l2Ali0/5q6PKKueyi7De64OwvpQVtAwvdnkN8Vi954pSbFwi0sK5rxbO39dnWB1uvsTqfRfPX18ABeepYazf9WTgBlX9oYTfn0+Q0zuCGIbvEJEFWPc+PFzqylXI+W2xo6qekRwr+biZWD3DZqmqW9UPfIGIrAImO3smVksi8jDWtbPXvBzHpVittgCsrtpfejMeVyaZGD7JeTphG9b9GHu8G41hGKUxp7kMnyMiT2KdrnreJBLDqB5My8QwDMOoMNMyMQzDMCrMJBPDMAyjwnymAmpVq1+/vkZFRXk7DMMwjGrl+++/P6qqZ1R1rrXJJCoqii1bfL52mmEYhk8Rkfii5pvTXIZhGEaFmWRiGIZhVJhJJoZhGEaFmWRSRulbt3L0tbmkby3zMAmGYRg1Vq29AF8e6Vu3snfkKDQ3F/Hzo8mzMwnu1g1baCgSEMDplcANwzBqD5NMyiB902ZOBDbleJ221D2xE717yl8L7XZsISHYQkKwh4ZgCw45NW0LDf3reUgIttAQ7MUtCwnFFhyE2MrXaIxbsZm9m+Jo0T2KqCGlVoo3Kujg7mT2/XGcpmfXpXErd4op1+w4qoKvv1Zfjq8yYzPJpAzSWnThh85RqNgQVdq0dRDerB752dlodjaanUN+dg6a43zu/Kk5OeQfyUH35aA5GWhOKvk5OWhOLpqTg+bmgoKeatgI4ucAhwNx+Dt/OsDPgTj8EIcD8fOzpv0c4HyenpLJgeMBQD1+XnaCZlu/JrTZGd3BK6TS216VeABP7zo9JZv4X5PQfBAbtDy3PsHhxQ05U3nSU7KJ/+XoX3F08k4cVSE9JZv4n333tfpyfK6x+TlsXHV3F48mFJNMyiBJGqA2a2RNFTicGkDWUUXEAeKwPqz8gQDrY6vgrJeIgDg/zAQEwU84/dMtLw/ycq3EkpsLeXnO5zmQa81X5/P83BzITLemc3JO/Txpj4TQltaBVTi57yh1WjaAcrZyPDW+oEerv3molpwn9pKemk3BEEqaD+nJWd5JJslZp8dxIovgMN/4APO09BO+/Vp9OT7X2PLy8tn3x3GTTLyl6dl18XPYyMvLx263cdmoDj7VjI1bsZkVHySRL3ZsmkfrbYuo+8MRwocMITwmhqAunc11HQ86uDuZj/699dT74eIbzvbK+8FX4qgKvv5afTm+wrE1PbuuR/dfa6sGd+3aVctzB7wvnw+F06+ZtOh3Lmlffkny8ljS1q3D0bgx4TFDiYiJIaBNG2+HWiP4yvvBV+KoCr7+Wn05Pk/EJiLfq2rXM+abZFI75B4/TurKVSTHLidjy/cEnHMOETExhMcMxdGokbfDMwyjmjDJpJDalkxc5ezbR/Inn5KyfDlZu3YR3K0b4cNiCL/8cuzh4d4OzzAMH2aSSSG1OZm4yvz9D1Jil5Mc+wl5R48S2q8v4THDCO3XF1tAgLfDMwzDx9ToZCIiVwNDgXDgdVVdVdo2JpmcTvPzyfjhB5KXx5L62WdoXh5hgwYRMSyG4O7dEbvd2yEahuEDfDaZiMgCIAY4rKrnuswfDPwHsAPzVXWmG/uqC7ygqqNLW9ckk+JpdjZp69eTEhtL6udfYA8PJ/yKKwgfNozAjh1MjzDDqMV8OZlcAqQBiwuSiYjYgT+AgUAisBkYjpVYnim0i3+o6mHndv8C3lLVH0o7rkkm7slLO0nqmtWkLI/l5IYN+LdsafUIGzYM/xYtvB2eYRhVzGeTCYCIRAGxLsmkFzBdVS93Tj8AoKqFE0nB9gLMBFar6hp3jlneZLLt8Da2HNpC10Zd6dywc5m3r85yjx4l5dMVJMfGkvnTTwSefx4RQ2MIv2IIfvXrezs8wzCqQHHJxFdvWmwKJLhMJwI9Slj/DuAyIEJE2qjqnKJWEpFxwDiAFuX4Vr3t8DZGfTaKfM0nwB7AvEHzalVC8atfn3q3jqDerSPIjo8nOTaW42+/zaFnnyWkVy8ihsUQeull2ENDvB2qYRhVrEaUoFfVWap6oapOKC6RONebq6pdVbVrgwZlr1m15dAWVPNRlKy8TDYf3FyhuKsz/5YtaTBpEq1WfErUu+8S0Lo1h154gZ19+rBvyhRSv/gSzc72dpiGYVQRX22Z7AOau0w3c87zqq6NuuIvdrLzc8gX4YfDP5Cbn4ufzVd/jZVPRAjqdC5Bnc6l4dT7SN+4keTYT9g/dSpitxM2eDARw2IIuuCCcldCNgzD95V6zURE3D0fdEJVU8oVxJnXTPywLsBfipVENgM3qeqv5dl/Ucp9zWThZWw5/jtN6rbhXyE2OtXvxHOXPIe/3TeKufmK/MxM0r5aS3Lsck6u/Rp7g/rW9ZWYGALbne3t8AzDKKdyX4AXkS+xiqyW1B9UgYWqurgcgb0D9APqA4eAx1T1dRG5AngJqwfXAlV9qqz7Lkm5kklmMjzfFvKywB5A4sSvGbv2blqEt+Df/f5NsCPYkyHWGHnJyaSsWkXK8ljSN28moG1bwofFEDF0KI4mTbwdnmEYZeDTvbm8oVzJ5Md34ZN7IDsN/ENh6IscbtufcavGER4QziuXvkKYf1jlBFxD5Bw8SMonn5Ac+wlZO3YQ1PVCImKGET74cux16ng7PMMwSmGSSSHlSiYLh0Lc+r+moy6GUbEczzzOhDUTUFXmDJxDvcB6ng22hsratYvk5bGkxMaSc/gwoX36WD3C+vfHFhTk7fAMwyiCR5OJiNwDdMW6gL9JVZ+veIhVq8RksnQkbF925ny7P+RlFzmdJsKkxg04YbMzN7wLjW54x/NB11CqSsbWbaTELidlxWdoVhZhAy8jPGYYIb16WqNKGobhE4pLJuXtXmNT1eGqeh1Q8+qXD3wCWvWHwtdA8rKLnQ5VZU5SGmf5RzDSdoSE1AQM94gIwRd0ofGjj9L267U0/feLaF4+iXfeyc5+/Tn41NNk/PQTtbUVbRjVQXlbJguBVVgtk8uxSph86tnQKpdbp7l2roaPboesFMhJL349RxAERMBVL5PTqh/3r7ufHw//yNxBc2ldp7VnA69F8k+eJPWLL60eYeu/wdGsqdUjbFgMAdHR3g7PMGolT5/mGsnpPby0PD25vMntayY5mfDFk7B5PuRmnrncLxC6jYEBj4AjEIDc/FymfzudtYlrmTNwDh0jO3o4+ton99gxUlasICX2EzK2biXw3HMJjxlK+BVX4GjY0NvhGUat4enTXHWAi1V1EdCiuiWSMnEEQlhjKK5Srtgg7KxTiQTAz+bHE72f4IroKxizcgzfH/q+ioKtufzq1aPezTcT9c7btF6zmrDLLuXEe++zq19/9v7jH5z43wfkpaZ6O0zDqLXKm0xa81ftrJrdF1YVvvsv5GRY0zY/qzVScNd7Trq1vFALzyY2pnWfxvD2w5mwegLf7PumigOvufybNaP+hAm0il1O9P/eJ6D9ORyZNYudvfuQOPkuUtesId+UcjGMKlXeZKJAkIicC9Tsu87iv7VuVgTrgnz0JXDrR9bPggv0mSdg74YzNhUR7rzgTiZ2nsgdX9zBmni3ChobbhIRAs85h0ZT76PNF5/TfO5c7OFh7H/wIXb2uZgDjzzCyY2b0Px8b4dqGDVeea+ZtAD+iXXN5BVV3evpwCqb29dM3r0Zfv8UQhrCVS9D24F/LSu4QH/yMLQbCje+Wexulv6+lGc2PcPjFz3Ola2v9MArMIqTn51N2tq1pMR+QtqXX2KvW5fwoUOJGBZDQPv2ZnAvw6gAT5egP6Cq00SkNVCuelzVyiVToc/dp10XAazEMvlHWP8iHCq5bNj17a4n2BHMI988wsmckwxvP7wSA67dbP7+hA8cSPjAgeSlppK6ajXJscs59sYb+LduRUSMVSPMv1kzb4dqGDVGeVsmz2PVzXoByFPVWzwcV6Xz1kiLX+z9gvvW3sfEzhMZ02lMlR+/Nss5dJiUFZ+SsjyWzF9/JahLF8KHxRA+ZAh+det6OzzDqBY83TX4NeBHYD1wi6pOrXiIVcubw/Zu2L+ByV9OZnj74dx1wV3mtIsXZO3eQ0psLMmxseTs309I74uIiBlG2KUDsAWbgp2GURxPdw3+Cqirqj8BOysSmKeISIiIbBGRGG/HUppeTXoxd+Bc3vv9PZ7a+BT5ai4QV7WAVtE0uPMOWq/8jKi33sS/eQsOzZzJH737sO/e+0hbuxbNyfF2mIZRbXi90KOILABigMMF45k45w8G/oNVgn6+qs4sZT9PAGnAdlWNLe243myZFNiRtIPxq8fTp2kfnuj9RK0eZMsXaG4uJzd8R0rsclJXr0ECAggfMoTwmBiCunQ2LUjDwIerBovIJVhJYLHL4Fh2rMGxBmKN/74ZGI6VWJ4ptIt/AOcDkUAgcLS6JBOA3cm7GbtqrBlky8fkZ2SQ9uWXJC+PJW3dOhyNGxMeM5SImBgC2rTxdniG4TWVlkxE5B/AblX9qgL7iOL0kRZ7AdNV9XLn9AMAqlo4kRRs/xQQAnQAMoC/qZ557khExgHjAFq0aHFhfHx8eUP2qH1p+xi7aizNQpvxUv+XzCBbPib3+HFSV64iOXY5GVu+J+Ccc5w9woaSs38/6Zs2E9y9G8Fdung7VMOodJ6+AH9FQWFHEbkea8jeVRUILorTk8m1wGBVHeOcHgH0UNXbS9nPKKpZy6TAkfQjjFs9jjD/MF6+9GXC/cO9HZJRhJx9+0j+5FNSli8na+dOsNlAFQkIoMUbC0xCMWo8j12AF5GhwE0icoXz+d8rkkg8SVUXupNIfFGD4Aa8cfkbZOdlM3rlaI5lHvN2SEYRHE2bUn/cWFot/5i6I0ZYZXRU0cxM0tau9XZ4huE15enNVR/rVFIDoB7wtEcjsuwDmrtMN3POq9HqBNZh/qD5hDhCGPXZKA6ePOjtkIwShF8xBAkIsFonIhx/913S1psabEbt5PUL8FDkaS4/rAvwl2Ilkc3ATapa8m3mZeBrp7lcZeRmcPdXdxOXHMe8gfNoHt689I0Mr0jfutW6ZnLhBZzcuJGjr86h3s030WDKFGwBAd4OzzA8zqP3mYjIMBFZKiLvishVFQzsHWAD0E5EEkVktKrmArcDK4EdwFJPJhJfF+QXxOz+s+kQ2YGRn41k1/Fd3g7JKEZwly7UHz+O4K5daTBpEi2XLCb18y+Iu+56Mv/4w9vhGUaVKe9Ni8NU9XpVvREYXJEAnMP/nqWqDlVtpqqvO+d/qqpnq2prVX2qIseojhx2B89d8hwXNbmI21bexq9JtSaXVmvBXboQvexDAjt0IO7a6zi25E0z3LBRK5Q3mQSJSAtn9eAQTwZk/KVgkK2hrYYyeuVothz0zdNyxunsoaE0mfkMTWY+w5HZs0kYN57cI0e8HZZhVKryJpPvsE5D3QH85LlwjMJsYuP+bvdzyzm3MHHNRNbvW+/tkAw3hV9xBa2WfYhmZrL7qqtJ/fJLb4dkGJWmvF2De2HV5/oKOONCjOFZIsLtXW7nn53/yZ1f3MmquFWkZOYwbvEWUjJN/Shf5mjShBYL36DeqFEk3jmZA48/Tn5GhrfDMgyPq0jX4PpUXtdgowi3nXsb07pP4/6v7+f5b95k1fZDfL7jkLfDMkohdjv1x40l6u23Sf92A3uuuZbM7du9HZZheFSZk4mqLlLVsaq6WFWXOCsHG1Xk+nbXM6PPDD5K+A+Out+wdHOit0My3BTU6VyiP/yA4AsvJO6GG0l6fYEZUtioMdwuUysi9YC7gYbAdqzCjMcrKzCjeBc3GUT2/h0EnPUW247lkJJxAeFBpkBkdWALDuasJ58g5JKLOfjIo6StX0eTmTNxNGrk7dAMo0LK0jJ5F0gFlgPBwHoR6V4pURklWrP9EPbMc8lIGIm93ufc/+UzpvtpNRM+cCDRH32EiI09V15FyiqfqEhkGOVWlmTSQFWfU9VYZ/XeYcCsSorLKMHSLQmczM4jL70t6XvH8O3hT5jx3QwzyFY142jUkObz5xE5cQL777mX/Q8/TP7Jk94OyzDKxe1yKiLyOXCXqv7sMu8nVT2vsoKrTL5cTqXApLd+4JOfD5wx32EXcvL++rv5Bx/A0XQ+eSfbkn28B37BceSmt2JIm568cvMFVRmyUU6Zv/3GvnvvRXNyaPrCCwR16uTtkAyjSJ4opzIJeEdEXhWRf4rIK8CfHovQOMO0Ie25uE19ghz20+a7JhKA7PSzSI+fgD1kJ8Et5+LfYBWhLedzdU/Tbbi6CGzfnuj33ye0z8XEDb+Jo3NeQ/PyvB2WYbjN7WSiqr8BFwBfYl2E34Y1+qFRSZrXC2bJmB7895YLaBgWcEZScaXZDdDUbgiKiGKz5ROf/nOx6xu+xxYYSONHHqbZy7M5tmQJ8SNHkrOvxhfLNmoIt5OJiAwAXgV6APHAD4BPXPUVEZuIPCUis0VkpLfj8bT+7Rry9dT+3NKzBYGOov9kgQ4bl0ddSoCfVak2X/NpX7d9VYZpeEhYv360+mgZtpAQdl/9N5JjP/F2SIZRqrKc5lqA1ZPrO6AV8ChQ4eqDIrJARA6LyC+F5g8Wkd9FZJeITCtlN1dhjXmSgzVmfI0T6LDTMCwQKWa5AB0jz2P+oPmMP288beq04YUtL3DopLmpsTryq1+f5nPm0OCuyRx46CH2TZ1KXmqqt8MyjGKV5QL8WlXt6/EARC4B0rDuWykYz8SONZ7JQKzksBnrlJodKDwO/D+cj+Oq+pqIvK+q15Z23OpwAd6VqtJr5hccTM4EwG4T/GxCbr6Sl2/9Dc+KCOTbaQMQEdJz0pmydgp7Tuxh3qB5tAhv4c3wjQrI2rmTfffeR35aGk2ef47gC0ynCsN7yn0BXkQWi8hdwHciMsXTganq10DhMWq7A7tUdbeqZmPd43KVqv6sqjGFHoexEk7BDZQ18qrlpj3HSMmwLqgHOez0bh3JW2N6cFHryFPXUpIzctgcZ/0agh3BzO4/m/ManMetK27l92O/ey12o2IC2rYl6r2lhA0cSPytIzkyazaam+vtsAzjNO6c5lqIdRalETBCROJF5GMReVJErqukuJoCCS7Tic55xfkAuFxEZgNfF7eSiIwTkS0isuVINSsJ/vo3e8jIyaNhWAD/veUCFo/uQdeoeiwZ/dcF+oycPBas331qG4fdwcyLZzKgxQBu++w2th7e6sVXYFSEzd+fRtPup/lrczjx3nvE33wL2QkJpW9oGFWk1GSiql8A/wH+UNUuQGvgIeB3oFvlhuceVU1X1dGqeoeqvlLCenNVtauqdm3QoEFVhlhxCpMHtOXrqf3p367haYsKLtDfOaDtGT0i7DY7j/R8hBva38C4VeNMCftqLrR3b6I//gh7g/rsuepqTny4zFQ/MHxCWa6ZfKuqF1VKEGeOAd8LmK6qlzunHwBw3nnvEdXtmoknvPHLG8zaOotn+jzD4OgKDZBpeJmqcuK99zj0zExC+/WlzjXXkPnrdoK7dyO4Sxdvh2fUYMVdM3G70CPwk4g8BjypWul1OzYDbUUkGtgH3AjcVMnHrPFuO/c2wv3DeWDdA6Rkp3B9u+u9HZJRTiJC3euvJ7hrNxInTSJhzFiw2RB/f1q8scAkFKPKlaVrcD2sD/X9IvKRp66ZiMg7wAagnYgkishoVc3FGslxJbADWKqqZhB0D7jm7Gt49pJnmblpJvN/nm9OkVRzAa2iiRg2zJrIz0ezs0nftNm7QRm1UqktE+cpp+9U9XrndADQEeiEdQPjexUJQFWLvIteVT8FPq3Ivo2iDYoaRKh/KHd9eRfJWclMuXAKIsXdwWL4uuBePZG5c9GsLMjPxz+qpbdDMmqhUq+ZiEjBXe9/AJ8Bn6nqwSqIrVLVxmsmhf145Ef+ueafXNbyMh7t+Sh2W/HlWgzflr51K+kbN3Ly2w3kHDpI1Ntv4xcZ6e2wjBqouGsmZbkA3x4YAlwORGDV6PoM+EZVq929HSaZWHYe38n41ePp3LAzMy+eib/dDLJVneVnZZEwegz5WVm0XPgGtpAQb4dk1DAVrhqsqr+p6r9VdTAwAFgPXAds9FyYRlVrW7cti4YsYkfSDiZ9Pon0nHRvh2RUgC0ggGb/fQXNyiJx8l1ojqkcbVSNMo8B79RfVT9V1TuAJz0ZkFH1moc1Z/GQxRzNOMrYVWNJzkr2dkhGBdjDw2k+by5Zu//kwMMPm3HmjSpR5mQiIkOBm0TkCudzU4a+BmgQ3ICFgxeCwKjPRnE4/bC3QzIqwNGoES3mzyftq7UcefFFb4dj1ALlaZnUBzKABljdhZ/2aESG10QERDBv4DwaBjfk1hW3kpBiynVUZwGtWtH8tTkce/Mtji1a5O1wjBquzMlEVRep6ljnzyWq+lNlBGZ4R7AjmNkDZtMxsiO3fmYKRFZ3QZ070/TfL3LohX+ZcVGMSlWuayYiMkxElorIuyJylaeDMrzL3+7Pc5c8R99mfblt5W1sO7zN2yEZFRDWvz9nTZ/OgQce4OS333o7HKOGKu8F+GGqer2q3giYIk81kN1m57Fej3Ht2dcybvU4vtn3jbdDMiqgzjV/p/6kSSTefgcZv5piEobnlTeZBIlICxFpAZiO7DWUiDDlwimMP288d3xxByvjVno7JKMCIsePI+Lqq0kYN57svXu9HY5Rw5Sl0KOr6cAdzudPeCYUw1eN7jSaiIAIpn09jdTsVK49u9SBLA0fJCI0euhBcpOS2Dt2rLlL3vCo8rZM2qnqfap6H9DBkwGVh7OVtMw5nnxp48Ub5XDt2dfyzCXP8PTGp3n959e9HY5RTmK30+S5Z3E0akzCuPHknzzp7ZCMGsLtciqnNvjr3pK3sUZgHOG8dlK+AEQWADHA4YLxTJzzB2MNymUH5qvqzFJiqquqb4rI/6nqDaUdt6hyKjk5OSQmJpKZmVnOV1PzZeVmcSzrGCGOEML9w70WR2BgIM2aNcPhcHgthuosLzWV+FtG4Fe/Ps1f/S/ib8roGO6pcG0ulx2NBPpglVPJB36sSPdgEbkESAMWuwyOZccqLDkQa8jezVgJzA4UHiDrH1jjvr8PKLBEVd8o7bhFJZM9e/YQFhZGZGSk21V0UzJzuHfpj7xw/fmEB9aOD7b0nHTiU+KJCIjgrJCzqrzisKqSlJREamoq0dHRVXrsmiTn0GHihw8n6MILafLsTMRW3hMVRm1S4dpcIlJPRJ4ELgJ+AT72xH0mqvo1cKzQ7O7ALlXdrarZwLvAVar6s6rGFHocBm4DHlPVAcDQ8saSmZlZpkQCsGb7IVZtP8TnOw6V97DVTrAjmOiIaFKzU0lMSyS/0sdKO52IEBkZaVqQFeRo1JDm8+dzct06Dr/wL2+HY1RzZfkq8i6QCiwHgoH1ItK9UqKCpoDr7deJznnF+Qy4U0TmAHHFrSQi40Rki4hsOXLkSHHrlCnQpVusMJduTizTdtVdoF8g0RHRZORmkJCaQF5+1RaONuOveEZAq2iavzaH4++8Q9IbC70djlGNlSWZNFDV51Q11jkW+zBgViXFVSaq+ouqXquqE1T13hLWm6uqXVW1a4MGDSp83JTMHH6IPwHA9/HHSc2s+gqty5YtY/v27eXe/quvvuJbN25kmz59Oi+88MJp8/zt/kSHR5OTn0N8ajy5+bkAJCUl0b9/f0JDQ7n99tvLFZfdbqdz58507tyZK6+8slz7MNwTdP75NP33ixx+8UWSl8d6OxyjmipLMjkmIp0KJlR1N1YLpTLsA5q7TDdzzvMpa7YfwmG3viE77MIaL5zqqqpkUhyH3UFUeBQoxKXEkZOXQ2BgIE8++eQZyacsgoKC2LZtG9u2bePjjz8u934M94T168dZjz/O/gcfJO0bc4OqUXZlSSaTgHdE5FUR+aeIvAL8WUlxbQbaiki0iPhjjT3vc58oS7ckcDLbOr1zMjvPI6e6rr76ai688EI6duzI3LlzT80PDQ099fz9999n1KhRfPvtt3z88cfcd999dO7cmT///JNt27bRs2dPzjvvPP72t79x/PjxYo8VFxfHnDlz+Pe//03nzp1Zt24dcXFxDBgwgPPOO49LL72UvW7c3OZn86NleEv8xI89KXtwBDro06cPgYGBZ6y7atUqevXqxQUXXMB1111HWlpaGX9DRmWp8/e/0eCOO9h3x51k/GLukjfKxu2bFlX1NxG5ALga696SbcA9FQ1ARN4B+gH1RSQR60L66yJyO7ASqwfXAlX12rt70ls/8MnPB86YX9AqKbAl/hhR084spje001m8cvMFbh1rwYIF1KtXj4yMDLp168Y111xDZDE3ll100UVceeWVxMTEcO211o2E5513HrNnz6Zv3748+uijPP7447z00ktFbh8VFcWECRMIDQ3l3nuts4PDhg1j5MiRjBw5kgULFnDnnXeybNmy07Z7/vnneeutt87Y38UXX8zUp6eyJ3kPLcPPHIf86NGjzJgxgzVr1hASEsKzzz7Liy++yKOPPnrGupmZmXTt2hU/Pz+mTZvG1VdfXcJvzfCUyLFjyD18mIRx44h65238W5rx5A33uJ1MRGQAcDNwAqs3109YXXErRFWLHA9FVT8FPq3o/ssrL185mpYFwJiLozmSlsVPiSfIzPmr51JO3ukvv/B0oMPG+c3qMPaSaA6lWD2P6ocGYLcVf/F41qxZfPjhhwAkJCSwc+fOYpNJYcnJyZw4cYK+ffsCMHLkSK677jq3ti2wYcMGPvjgAwBGjBjB1KlTz1jnvvvu47777itye1Vl/8n9xCXHkZWXddqy7777ju3bt9O7d28AsrOz6dWrV5H7iY+Pp2nTpuzevZsBAwbQqVMnWrduXabXYpSdiNDowQfITTrK3jFjiXrnbfzq1/d2WEY1UJZyKguAuwAHcB5WC6Uj0MbjUfmAo2lZ9Hj68wrtIzMnn417jnH1K39dk9j44KU0Cj/z9A9Y1y/WrFnDhg0bCA4Opl+/fqe6v7r2XvJ2l9jiWiaXXHIJs2bNoklIEw7JIZIyksjJ+6tTgqoycOBA3nnnndO227hxI+PHjwfgiSee4Morr6RpU6vzXqtWrejXrx9bt241yaSKiM1Gk2efJWHsOBLGjafF4sXYQ00JPqNkZUkm8aq6zPn8vUqIxafUDw1g44OXnjE/KzePOWv/5H/f7yMr98z7KwL8bFxzYVMm9G1NgJ+9yP0WJzk5mbp16xIcHMxvv/3Gd999d2pZo0aN2LFjB+3atePDDz8kLCwMgLCwMFJTUwGIiIigbt26rFu3josvvpglS5acaqW8/PLLAGf0rgoLCyMlJeXU9EUXXcS7777LiBEjeOutt7j44ovPiLOklgk4v90GNyLMEUZqTiopWSmEB4TTs2dPJk2axK5du2jTpg0nT55k37599OjRg23btp3a/vjx4wQHBxMQEMDRo0f55ptvimwhGZXH5u9Ps5dnEz/iVvbdeQfN58wxd8kbJSo1mYjIYuAH4DsRmaKqtWIMULtNim1BREeGUtyZKptAq/qhtKhX9m9ygwcPZs6cOZxzzjm0a9eOnj17nlo2c+ZMYmJiaNCgAV27dj114frGG29k7NixzJo1i/fff59FixYxYcIE0tPTadWqFW+8YRUD+O23306dXnI1bNgwrr32Wj766CNmz57N7Nmzue2223j++edp0KDBqe3LKjo6mpSUFLKzs/ni0y/48JMP6dWlFwsXLmT48OFkZVmnwGbMmMHZZ5992rY7duxg/Pjx2Gw28vPzmTZtGh06eL0EXK1jDwuj+dzXiB9+E/sffIgmzz1r7pI3ilVqORXntZLzXR71gB+dj59UtVq2Uooqp7Jjxw7OOeecErdTVXrN/IKDydapJrtN8LMJuflKXr71uzwrIpBvpw3wqRvrYmJi+OCDD/D3wrfL5Kxk9qXto2FwQ+oHee78uzt/L6PisvbsIf6mm4m46ioaTbvf2+EYXlbuciqq+gVWwcU/VLUL0Bp4CPgd6ObpQH3dpj3HSMmwrgMEOez0bh3JW2N6cFHrSIIc1mmt5IwcNscV3yXXG2JjY72SSMAaW755WHMOpx/m0MlDlLUenOFdAdHOu+T/7/9Ien2Bt8MxfJRbbVZVzceq7Iuq5jprZL2pqrXuRPbr3+whIyePhmEB/PeWC1g8ugddo+qxZHQP/nvLBTQMCyAjJ48F63d7O1SfEuYfRlR4FMcyj3Hg5AGTUKqZoPPOo9l/XuLwSy+RbG4iNYpQlhOgP4nIYyJSu0+aKkwe0Javp/anf7uGpy3q364hX0/tz50D2la8z3QN5O0CkUbFhF5yCWc9+QT7H3qYtHXrvR2O4WPK0purHtAXmCgiG7HuM6m210zKa+6thU4VZibDsolw9asQGEGgw87dA88uemODQL9AoiKiiE+JJyE1geZhzbHV8u8n1Umdq68m98gREidPpuWiRQR1Orf0jYxaodT/xSLSS0REVa9X1XOAlsDjwC6gR2UH6PN+XwG/fQK/f+btSKqNAHuAVSAyL4e4lLhTBSKN6iFyzBjqXHsNCePHkx0X5+1wDB/hzlfCW4HvReRdERmFNaLhD6q6qKQKvbXG1iWn/zTc4rA7iIqI+qtAZH7VV1w2ykdEaDRtGiE9e7B3zFhyixnOwahd3OnNNVFVLwCmA3WBhSKyQUSeFpFLnKMi1k6ZyZCw2XqesAkyU0pevxJ4swR9cdwtQe9aIDIuOY7svGwAtm3bRq9evejYsSPnnXce//d//+feizGqjNhsnDVzJo7mzdg7fjx5pmBnref2yWpV/U1V/62qg4EBWMP2XgdsrKzgiiMirUTkdRF532VeiIgsEpF5InJzlQTy+wqwO4fqtTus6Srm7RL0RSlLCXq7zU6L8BYE2APYk7yHzNxMgoODWbx4Mb/++iufffYZd911FydOnPBojEbF2fz9aTZ7Nigk3nEH+dnZ3g7J8KJyXflU1QysURdtRd28UhIRWSAih0Xkl0LzB4vI7yKyS0SmlXL83ao6utDsvwPvq+pYoGpGU9q6BLKd38iy0zxyqqs6lqAvLCQkpEwl6G1io3lYc0IdocSlxNEsuhlt27YFoEmTJjRs2JDiRsY0vMseGkqLua+Rk5DIgWkPoPmmh15tVZbeXIhIF+Am4HrgINAea5yTslgIvAwsdtmvHXgFGIg1RO9mEfkYq/z8M4W2/4dz3PfCmgE/O597dgzZpSNh+7Iz59sL3QSYsBGmR5y5Xoer4fpFbh2qOpegLyj0WJySStCLCE1Cm3Ao/RDxKfFWcvEPZdOmTWRnZ5sijz7Mr0EDWsyfR9xNN3P42WdpOG2aT1V/MKqGO7W5zgaGYyWRVKwij/1UdY+I7CnrAVX1axGJKjS7O7DLOXojIvIucJVzeOAYN3ediJVQtlHOFtdp8vMgzZmzet0OaUdg//eQ61KxN69Qs77wtF8gNO0KF90BKc7xUEIbgq34y0zVvQR9SUorQV9QINIudvam7sUv1Y8RI0awaNEibKYmlE/zj4qi+WtziB85Cr+GDYkcXfjEgVHTudMy+Q1r5MNrVfXnQss8dW9eUyDBZTqRErodi0gk8BTQRUQecCadD4CXRWQosLyY7cYB4wBatGhRckRph+HF9mV4CUXIzYT49TDfpfrwlN8g/KwiV68pJeiLU5YS9CdTT3LFlVfwwGMPnFbw0vBdQZ060ew//yHhn//Er359Iq66ytshGVXInWTyd6xhc1eJyBpgKfCZqnqtL6eqJgETCs07CdxWynZzgblgFXos8SChDa0P/sJys+Cbl+DHd05vpRTwC4Tzh0Pvu8CviHLzoQ3PnOdUU0rQF8fdEvTZ2dmMvmk0I24dQffLu3M046hHC0QalSf04j40mfEk+x9+BHu9eoQW8f4xaqZSk4lzDJNlIhICXIX1zX6+iHwKhHsojn1Ac5fpZs553mOzF9uCILI1FHdOWGwQ2QbqRZX5kDWpBH1UVNSpEvTLli1j1apVdOjQwa0S9EuXLuXrr78mKSmJd5a8Q3Z+Nq/MfYUBPX2rErNRtIirriL36FES75xMy0ULCTrvPG+HZFSBUkvQF7mRSF2sbsE3quqAcmwfBcSq6rnOaT/gD+BSrCSyGbipMsd9L28JelTh3x0gZb81bfOzHvm51gMgvCnc/WvxCccLvFmCvqJO5pxkb8peIgIiOCvkrFMJxZSg912qyuGZz5L88ce0fPstAqKjvR2S4SHlLkFfFFU9rqpzy5lI3gE2AO1EJFFERqtqLnA7sBLYASytzERSIfHfWjcrAjiCIfoSuPUj66cj2JqfeQL2bvBaiEXxZgn6igpxhBAVEUVKdgr70vaZApHVgIjQ8P6phFx0EQnmLvlawZ3eXKVcqT7lhKqWegu4qg4vZv6nwKduHst7vvsv5GRAaGO46mVoO9CaP+JD2LkaProdTh6GDf+Flhd5N9YaJMgviOiIaOKT/yoQafg2sdlo8szTJEyYwN5x42m5ZDF2l/uljJrFnQvwi7B6bZV0zkax7h9ZXMI6NcclU6HP3eAodFNe24Ew+UdY/yIc8s2GVXUWYA+wEkpKPPEp8aaFUg2Ivz9NZ81m7623knj7HTSf+xq2atpCNkpWrmsmNUG5r5kYXpebn0t8Sjz7/txHhw4dTE+vaiD36FHiht9E4Lkdafqvf5mx5Ksxj14zMQxv8rP5ERUehU1sjFwxkv1p+70dklEKv/r1aTF/HukbN3HomZlmpM0ayCQTo1qy2+zUC6xH6zqtGbFiBH+e+NPbIRml8G/ZkuavvcaJ//2PpPnzvR2O4WEmmfiYUHOB0m0iwov9XqTnWT0Z9dkofjn6S+kbGV4V1Olcms2axZFZsznx4TJvh2N4kEkmxhlyc3NLnPYlfjY/nuz9JDGtYhi9cjSbDmzydkhGKUL79KbJ009x4NFHSfv6a2+HY3iIO12DUym6BpcAqqqeugveKMby5cuZMWMG2dnZREZG8tZbb9GoUSOOHDnCTTfdxP79++nVqxerV6/m+++/Jy0tjZiYGH75xfqm/sILL5CWlsb06dP5888/mTRpEkeOHCE4OJh58+bRvn17Ro0aRWBgIFu3bqV3794cO3bstOlJkyYVuZ0vsImNqd2mUiegDhPXTOT5vs8zoEWZb4EyqlDEsGHkHjlK4uS7aLnwDYLOP9/bIRkVpaq18nHhhRdqYdu3bz9jXmm2Htqq836ap1sPbS3ztkUJCQk5Y96xY8c0Pz9fVVXnzZunU6ZMUVXVSZMm6dNPP62qqitWrFBAjxw5onv27NGOHTue2v7555/Xxx57TFVVBwwYoH/88Yeqqn733Xfav39/VVUdOXKkDh06VHNzc4ucLm47byrq7/X2jre18+LOumznMi9EZJTVwZnP6u89emrm7t3eDsVwE7BFi/hMLet4JnWBtsCpGyxUtUa2U/Py80jKTCpxnV+TfuW+tfeRk5eDw+7g+b7P0zGyY4nbRAZGYi+hBH1REhMTueGGGzhw4ADZ2dlEO0tTrF+//lS5+sGDB1O3bt0S95OWlsa33357Wln6ghpZANdddx12u/2M6dK28yXD2w8nzD+MR755hJTsFEZ0GOHtkIwSNLzvXnKPHiVhzFhavvM2jobFF0I1fJvbyURExgCT+WvMkJ5YZVFq5PmEpMwkLn3v0tJXdMrKy+LOL+4sdb3Pr/uchsFl+w9zxx13MGXKFK688kq++uorpk+fXuL6fn5+5LuMeFdQsj4/P586deqcVqHXVUhISJHTpW3na2JaxRDmCOOetfeQnJXMpM6TTIFIHyU2G02emkHChIkkFNwl76yIbVQvZWmZTAa6Ad+pan8RaQ88XTlheV9kYCSfX/d5ieuUt2VSVsnJyTRt2hSARYv+GrGxd+/eLF26lPvvv59Vq1adGqK3UaNGHD58mKSkJEJDQ4mNjWXw4MGEh4cTHR3Ne++9x3XXXYeq8tNPP3F+Keery7udN/Vt3pc5l83hji/uIDkrmQd6PIBNTH8TX2TdJT+LvSNHkjjpdprPn2fukq+GypJMMlU1U0QQkQBV/U1E2lVaZF5mt9lLbUE0DG7I/EHz2XJoC10bdaVzw84VPm56ejrNmjU7NT1lyhSmT5/OddddR926dRkwYAB79lgDXD722GMMHz6cJUuW0KtXLxo3bkxYWBgOh4NHH32U7t2707Rp09MulL/11ltMnDiRGTNmkJOTw4033uhWUijvdt7UtXFXXr/8dSaumUhKdgoz+szAYXN4OyyjCPbQEJrPfY244Texf+r9NP3XC4i9bKeDDe9yu5yKiHyINfjUXVinto4DDlW9otKiKz6WVsBDQISqXuucdzUwFGuMlddVdVVJ+6gJ5VSysrKw2+34+fmxYcMGJk6cWG1ORXmCu3+vuOQ4xq4ey9l1z+Zfff9FoF9gqdsY3pG9dy9xw28ifPBgGj38kDk96YMqXE5FVf+mqidUdTrwCPA6cHU5AlkgIodF5JdC8weLyO8isktEppUSy25VHV1o3jJVHYs1AuMNZY2rOtq7dy/dunXj/PPP584772TevHneDsknRUVEsWTIEhJSExi/ejyp2aneDskohn+LFjSf+xrJH35I0lzzfq5OyjueyVpV/VhVs8ux+UJgsOsMEbEDrwBDgA7AcBHpICKdRCS20KO0q9cPO/dV47Vt25atW7fy448/snnzZrp16+btkHxW45DGLBy8kMy8TEavHE1SRsk99QzvCerYkaazZ3Hk5Zc58b8PvB2O4aZSk4mIrHf+TBWRFJdHqoiUOn5JYc6uxMcKze4O7HK2OLKBd4GrVPVnVY0p9DhcTJwiIs8CK1T1h2LWGSciW0RkyxEzWE+tUy+wHq8Pep0QRwijPhvFgbQD3g7JKEZo7940eeYZDkyfTupXX3k7HMMNpSYTVe3jfPqqqoa7PMKAOR6KoymQ4DKd6JxXJBGJFJE5QBcRecA5+w7gMuBaEZlQ1HZqjQ7ZVVW7NmjQwEOhG9VJqH8or172KlHhUYxYMYLdybu9HZJRjIiYoTS8Zwr77rqbjFp0LbC6KstprsuKmDe4iHmVTlWTVHWCqrZW1Wec82ap6oXO+Z5KckYNFOgXyIv9X6R74+6MWjGKX5PMQGa+KnLUKOrefBMJ4yeQtdskfl/mzmmuiSLyM9BeRH5yeewBfvZQHPsA13FYmznnGUalcNgczOgzgyHRQxi9cjSbD272dkhGMRrecw8hfS8hYcxYcg4VeZbb8AHutEzeBoYBHzl/FjwuVNWbPRTHZqCtiESLiD9wI/Cxh/ZdbSQlJdG5c2c6d+5M48aNadq06anp7Ozy9HVwz5gxY9i+fXuJ6xw5coQePXrQpUsX1q1bV2mxVCWb2JjWfRojO4xkwuoJfLn3S2+HZBTBukv+KfxbtSJh7FjS1n/D0dfmkr51q7dDM1yUetOiqiYDycBwTxxQRN4B+gH1RSQReExVXxeR24GVgB1YoKq17txDZGTkqftEpk+fTmhoKPfee++p5bm5ufj5lamcmlvmuzFQ0eeff06nTp3cWrc6EREmdp5IeEA4U9ZO4YmLnmBY62HeDqvMth3e5tGbZ32NOBw0+89L7LnuehLGjgURxN+fFm8sILhLF2+HZ+BeCfr1qtqniFL05SpBr6pFJiVV/RT4tCz7qg0Kl4a/8cYbmTx5MpmZmQQFBfHGG2/Qrl07Fi5cyIcffkhycjL79u3jlltu4bHHHmPOnDnMmWNdQkpOTiYqKoovvzz9G3i/fv144YUX6Nq1K6GhoUyePJnY2FiCgoL46KOPOHDgAFOnTiUjI4MtW7awYcMGgoKCvPHrqDQ3n3Mz4f7hPPrto6Rkp3DzOZ5qdFe+b/Z9w+1f3E5efh4B9gDmDZpXIxOKLSSEsIEDSXrtNVBFs7JIXbPGJBMf4U7LpI/zp6m+VoT0rVtJ37SZ4O7dKu1NnZiYyLfffovdbiclJYV169bh5+fHmjVrePDBB/nf//4HwKZNm/jll18IDg6mW7duDB06lAkTJjBhwgRycnIYMGAAU6ZMKfFYJ0+epGfPnjz11FNMnTqVefPm8fDDD/PEE0+wZcsWXn755Up5jb5gWOthhPmHce/ae0nOSmbi+RN99g7sE5kn+CLhC1bFrWLD/g3kYxX2zMnPYcuhLTUymQCE9uvLsUWL0OxsUOXYgjfISdxHvZEjCerS2Wf/XrWB58+Z1BCal0fu0ZJvbMv45Rf2T5mC5uQgDgdNXnyRoHPPLXEbv/qRZa455FoaPjk5mZEjR7Jz505EhJycnFPrDRw4kMhIq5Dk3//+d9avX0/XrlbVg8mTJzNgwACGDSv5FI6/vz8xMTEAXHjhhaxevbpMsVZ3/Zr349XLXj1VIPL+7vf7TIFI1wSy8cBG6gfXZ1DLQVzW8jJmbppJbn4uDpuDro3OqHRRYwR36UKLNxac+gIndjvHFi0m/tZbCezQgXojbyV80CDEYWqwVbWylKAPAK4Boly3U9UnPB+W9+UeTWJX375ur69ZWeybNKnU9dqsXYujUdlK0LuWhn/kkUfo378/H374IXFxcfTr1+/UssLfygqmFy5cSHx8vFutCofDcWo7u93u00P2VpZujbtZBSJXWwUin+j9hNcKRJ7IPMHnez9nVfwqNh3YRIPgBgxsOZB/dv4nnep3OvW3alOnTY2+ZuIquEuX084CNP3XCzS8716Ov/UWB594ksPPv0Ddm2+i7vXXY4+I8GKktUtZWiYfYV2I/x7wzZGRPMivfiRt1q4tcZ3ytkwqwrUc/cKFC09btnr1ao4dO0ZQUBDLli1jwYIFfP/997zwwgusW7cOm803vmFXBx0jO7JwyELGrx7PlC+n8Hzf56usQKRrAtl4YCMNgxsyqOUgJnWedFoCcdW5Yecan0RK4mjcmIb33EP9iRNJ/ugjji1azNH/vkqdv11N3REjCHAOKGdUnrIkk2aq6pWbFL1B7PZSWxCORgPwW/hGpV8zcTV16lRGjhzJjBkzGDp06GnLunfvzjXXXENiYiK33HILXbt25bbbbuPYsWP0798fgK5du9a4HlmVpVVEKxYPXsy41eOYuGYiswfMJtQ/tFKOdTzzuJVA4lax6eAmGgU3YmDLgdze+XbOrX+uuRbgJltwMHWHD6fODTeQ9vXXHFu0iN1XDCW0b1/qjRpJcI8e5ndZScpSgn4uMFtVPXWjolfVhBL0rhYuXFjjL5AXVlV/r6SMJCaumQjAnIFzqBdYzyP7LSqBDIoaxKCWg0wC8aDM3//g2OJFpCyPxT86mnq33kr4sBgzAFc5VbgEPdAH+MFZJv4nEflZRH7yXIiG4ZsigyJ5/fLXCXYEM3LFSA6ePFjufR3PPM77f7zP2FVj6b+0P/N+mke7eu1484o3+eyaz7in6z10alD0qSyjfALbnU2Tp56izZdfEDZwIIdffJFd/Qdw5OVXyE0y1aM9pSwtkxY47y1xna+qeyshrkpX01omtVFV/70yczO5d+29/H78d+YOnEt0hHvn4Y9lHjvVAtl8cDONQxozqOUgBkUNomNkR5M4qlh+VhYpn3zKsUWLyN6zh/BhMdS7dSSB7c72dmjVQnEtk1KTSRE3K55aRDluWvQVJplUf974e+Xk5/DIN4+wYf8GXr3sVTpEdihyvYIEsjJuJVsObjmVQC6PupwOkR1MAvEBqkr6xo0cW7iItLVrCe7Zg3ojRxJ6ySWI6axSrOKSiTs3LZqbFQ3DyWFz8HSfp5m5aSb/WPkP7rrgLtJy0ujaqCstwluwJn4Nq+JX/ZVAogZx9wV3mwTig0SEkJ49CenZk6w9ezi+ZAn77p6Co3Fj6t06goirrsIWHOztMKsNt09z1TSmZVL9efPvpao8+u2jLNu1DHH+U5QmoU24POpyBkUNokM9k0Cqm7wTJzjx/vsce/Mt8jMyqHv99dS95WYcjRp5OzSf4YkL8D5DRFqJyOsi8n6h+SHOkRRjvBWbUTuICC3DW55KIopyQ/sbWPH3Fdx94d3mWkg1Za9Th8gxY2izehWNH32Ekxs3suvSy9h3z71k/FwjOrJWmipPJiKyQEQOi8gvheYPdvYU2yUi00rah3N439FFLLofWOrJeKuSt0rQG+XTtVFXAuwB2MVOgD2AodFDTQKpIcThIGLoUKL+711aLl6M5uURd8ONxN10MykrV6F5ed4O0ed4ozbXQuBlYHHBDBGxA68AA7GG7N0sIh9jlaN/ptD2/yhqHHgRGQhsB6rmNuVK4K0S9Eb5dG7YmXmD5tWaMia1kYgQfEEXgi/oQs6+fRx78y0OPPQQh597jrq33EKda6/BHmYuK4MXkomqfi0iUYVmdwd2qepuABF5F7jKOSSvu6es+gEhQAcgQ0Q+VdV8z0TtPVVRgt4ov9pexqQ2cTRtSqP7p1J/0iSSP/iAY0uWcPTll4m45u/UGzEC/+bNS99JDeYrX3ObAgku04lAj+JWFpFI4Cmgi4g8oKrPqOpDzmWjgKNFJRIRGQeMA2jRooVHAj+4O5l9fxyn6dl1adyqcorKVWUJesMwSmYPDaHerSOoe/NNpH31FccWLuLPywcTdukA6o0cCTY76ZurrsSSr/CVZFImqpoETChm2cIStpsLzAWrN1dJx8jPVzJSSr5OcTg+hZXzfyUvNx+7n43Lx3SkYcuSb7sJCvfHZivbefWqLEFvGIZ7xG4n7NJLCbv0UjJ+/ZXjixcTP3IU5OeDKvj50fihBwkfMgR7nTreDrfS+Uoy2Qe4thGbOed5TUZKNgunfeP2+nk5+Xz6aum9PUbN7E1InYAyxVKVJegNwyi7oI4dCXr2WeyNGnFsnrOQam4uB2c8xcHHn8DeoD4BbdoQ0Lat9bNNWwLatqlR11t8JZlsBtqKSDRWErkRuMmbAQWF+zNqZu8S1ylomeTn5mMrQ8ukIkwJesPwXWH9+3N88ZJTw1K0eH0+9nr1yNq1i+xdu8jauYvjG74jKy4OcnLwa9ToryTTtg0Bbdrg37oN9tCQUo/la6o8mYjIO1gXy+uLSCLwmKq+LiK3AyuxenAtUNVfqzo2VzablNqCiK7TgKvv7lLp10xcmRL0huG7Co8EWXDNJCA6GgYOPLWe5uSQvXcvWTt3kbXLeqStW0d2fDzk5uLX5CyXFoyzNdO6lU/fkW/ugHdRne+ANyXoDaP60+xssuLirFbMrl1k7dxJ1s5dZO/dC6o4mjY9syXTqhW2wKq7I6LctbkMwzCMqiH+/gSefTaBZ59ewTg/K4vsPXv+asns3EnKypXkJCSACI7mzaxWjEui8Y+OrtIxW0wyqSFGjRrFqFGjvB2GYRiVwBYQQGD79gS2b3/a/PyMDLJ273Zej9lJ1m+/kRIbS86+fWC349+ihTPBWElGc3PJSUgkuFdPj3dbNsnEMAyjmrIFBVk9yTp2PG1+/smTZO3eTdYfO8natYuMn3/h+P8tJc85GJjMnUuLNxZ4NKGYZFKIqpr6StVAbb3WZxjusIWEENSpE0GdOp2ad/S1uRz5z38gPx/NybE6CXgwmZi+oi4CAwNJSkoyH1Q+TlVJSkoisAovOhpGdRfcvRvi7w92O+JwENy9m0f3b1omLpo1a0ZiYiJHjhzxdihGKQIDA2nWrJm3wzCMaqO4bsueYpKJC4fDQXS0e+N6G4ZhVDfBXbpUWr0wc5rLMAzDqDCTTAzDMIwKM8nEMAzDqLBaW05FRI4A8UAEkFzGzd3dxp31SlunpOXFLasPHC01Ou8pz++8qvZdme8Hd9ctz9+8pGXm/VD1+/CFz4jKej+0VNUGZ8xV1Vr9AOZW1jburFfaOiUtL24ZsMXbv1dP/86rat+V+X7wxHvCvB+qft/V9TOiqt8P5jQXLK/EbdxZr7R1Slpenth9QWXGXdF9V+b7wd11y/s3N++Hytl3df2MqNL3Q609zVWTicgWLaKqp1E7mfeD4aqy3g+mZVIzzfV2AIZPMe8Hw1WlvB9My8QwDMOoMNMyMQzDMCrMJBPDMAyjwkwyMQzDMCrMJJMaTkRaicjrIvK+t2MxfIOIXC0i80Tk/0RkkLfjMbxLRM4RkTki8r6ITCzvfkwyqYZEZIGIHBaRXwrNHywiv4vILhGZBqCqu1V1tHciNapKGd8Ty1R1LDABuMEb8RqVq4zvhx2qOgG4Huhd3mOaZFI9LQQGu84QETvwCjAE6AAMF5EOVR+a4SULKft74mHncqPmWUgZ3g8iciXwCfBpeQ9okkk1pKpfA8cKze4O7HK2RLKBd4Grqjw4wyvK8p4Qy7PAClX9oapjNSpfWT8jVPVjVR0C3FzeY5pkUnM0BRJcphOBpiISKSJzgC4i8oB3QjO8pMj3BHAHcBlwrYhM8EZghlcU9xnRT0RmichrVKBlYkZarOFUNQnr3LhhAKCqs4BZ3o7D8A2q+hXwVUX3Y1omNcc+oLnLdDPnPKP2Mu8Jw1Wlvh9MMqk5NgNtRSRaRPyBG4GPvRyT4V3mPWG4qtT3g0km1ZCIvANsANqJSKKIjFbVXOB2YCWwA1iqqr96M06j6pj3hOHKG+8HU+jRMAzDqDDTMjEMwzAqzCQTwzAMo8JMMjEMwzAqzCQTwzAMo8JMMjEMwzAqzCQTwzAMo8JMMjGMSiQiz4hIf+cYIkXWRhORhSJybRHzu4qIKXtiVAsmmRhG5eoBfAf0Bb4uy4aqukVV76yUqAzDw0wyMYxKICLPi8hPQDesO5HHAK+KyKPFbHKZiGwRkT9EJMa5j34iEut8Pt054NFXIrJbRO50zg8RkU9E5EcR+UVEzGBXhleYqsGGUQlU9T4RWQrcCkwBvlLVkkaxi8Iab6I18KWItClinfZAfyAM+F1EXsUaAGm/qg4FEJEIz70Kw3CfaZkYRuW5APgRKwnsKGXdpaqar6o7gd3ObQr7RFWzVPUocBhoBPwMDBSRZ0XkYlVN9mD8huE20zIxDA8Tkc5Yw6Y2A44CwdZs2Qb0UtWMIjYrXCSvqKJ5WS7P8wA/Vf1DRC4ArgBmiMjnqvpExV6BYZSdaZkYhoep6jZV7Qz8gTXW9hfA5arauZhEAnCdiNhEpDXQCvjdnWOJSBMgXVXfBJ7Hag0ZRpUzLRPDqAQi0gA4rqr5ItJeVbeXssleYBMQDkxQ1UwRcedQnYDnRSQfyAEmViRuwygvU4LeMAzDqDBzmsswDMOoMJNMDMMwjAozycQwDMOoMJNMDMMwjAozycQwDMOoMJNMDMMwjAozycQwDMOoMJNMDMMwjAr7fzJSORvbhYGfAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(1, 1, figsize=(6, 3))\n", + "ax.set_title(f'Initial density error, m={masses}')\n", + "ax.set(xlabel='# bins', ylabel=r'Initial $|\\Delta /rho_\\mathrm{tot}/rho_\\mathrm{tot,sep}|$', xscale='log', yscale='log')\n", + "\n", + "ax.plot(results['auto']['bins'], np.abs(results['auto']['error']), '*-', lw=1.2, ms=16, label='auto, tol=1e-5')\n", + "ax.plot(results['auto_tol']['bins'], np.abs(results['auto_tol']['error']), '*-', ms=16, lw=1.2, label='auto, tol=1e-2')\n", + "ax.plot(results['laguerre']['bins'], np.abs(results['laguerre']['error']), '.-', lw=1.2, label='Laguerre')\n", + "ax.plot(results['trapz_inf']['bins'], np.abs(results['trapz_inf']['error']), '.-', lw=1.2, label='Trapz inf')\n", + "ax.plot(results['trapz']['bins'], np.abs(results['trapz']['error']), '.-', lw=1.2, label='Trapz')\n", + "\n", + "ax.legend()\n", + "fig.savefig('precision.pdf')" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "b6e41f34", + "metadata": {}, + "outputs": [], + "source": [ + "def get_binning(bins, strategy=3, auto_tol=None):\n", + " lag_cut = 90\n", + " collective = Class()\n", + " if strategy == 0:\n", + " # +10 bins is approximately required to reduce tol_ncdm_bg by 1 order of mag\n", + " if auto_tol == 'bin':\n", + " auto_tol = 1e-5*0.1*(bins/10)\n", + " collective.set(standard_settings | {'collective_ncdm': 1, 'output_ncdm_binning': 1, 'Quadrature strategy': f'{strategy}, {strategy}, {strategy}', 'tol_ncdm_bg': auto_tol})\n", + " elif auto_tol is None:\n", + " new_settings = deepcopy(standard_settings) | {'collective_ncdm': 1, 'output_ncdm_binning': 1, 'Quadrature strategy': f'{strategy}, {strategy}, {strategy}'}\n", + " del new_settings['Number of momentum bins']\n", + " del new_settings['Maximum_q']\n", + " collective.set(new_settings)\n", + " else:\n", + " collective.set(standard_settings | {'collective_ncdm': 1, 'output_ncdm_binning': 1, 'Quadrature strategy': f'{strategy}, {strategy}, {strategy}', 'tol_ncdm_bg': auto_tol})\n", + " else:\n", + " if strategy == 1:\n", + " if bins > lag_cut:\n", + " bins = lag_cut\n", + " qmax = 1e+2\n", + " collective.set(standard_settings | {'collective_ncdm': 1, 'output_ncdm_binning': 1, 'Number of momentum bins': f'{bins}, {bins}, {bins}', 'Quadrature strategy': f'{strategy}, {strategy}, {strategy}', 'Maximum_q': f'{qmax}, {qmax}, {qmax}'})\n", + " collective.compute()\n", + " bg = collective.get_background()\n", + " print(bg.keys())\n", + " if strategy == 0:\n", + " # Get amount of bins\n", + " max_bin_idx = np.argmax([bg[f'q_ncdm[0][{idx}]'][0] for idx in range(1000) if f'q_ncdm[0][{idx}]' in bg])\n", + " print(f\"Automatic quadrature has {max_bin_idx} bins at target {bins} bins.\")\n", + " qlist = [bg[f'q_ncdm[0][{idx}]'][0] for idx in range(max_bin_idx)]\n", + " wlist = [bg[f'w_ncdm[0][{idx}]'][0] for idx in range(max_bin_idx)]\n", + " return np.array(qlist), np.array(wlist), max_bin_idx\n", + " else:\n", + " qlist = [bg[f'q_ncdm[0][{idx}]'][0] for idx in range(bins)]\n", + " wlist = [bg[f'w_ncdm[0][{idx}]'][0] for idx in range(bins)]\n", + " return np.array(qlist), np.array(wlist)\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "ce2b24de", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "dict_keys(['z', 'proper time [Gyr]', 'conf. time [Mpc]', 'H [1/Mpc]', 'comov. dist.', 'ang.diam.dist.', 'lum. dist.', 'comov.snd.hrz.', '(.)rho_g', '(.)rho_b', '(.)rho_cdm', '(.)rho_ncdm[0]', '(.)p_ncdm[0]', 'q_ncdm[0][0]', 'w_ncdm[0][0]', 'q_ncdm[0][1]', 'w_ncdm[0][1]', 'q_ncdm[0][2]', 'w_ncdm[0][2]', 'q_ncdm[0][3]', 'w_ncdm[0][3]', 'q_ncdm[0][4]', 'w_ncdm[0][4]', 'q_ncdm[0][5]', 'w_ncdm[0][5]', 'q_ncdm[0][6]', 'w_ncdm[0][6]', 'q_ncdm[0][7]', 'w_ncdm[0][7]', 'q_ncdm[0][8]', 'w_ncdm[0][8]', 'q_ncdm[0][9]', 'w_ncdm[0][9]', '(.)rho_lambda', '(.)rho_ur', '(.)rho_crit', '(.)rho_tot', '(.)p_tot', '(.)p_tot_prime', 'gr.fac. D', 'gr.fac. f'])\n", + "Automatic quadrature has 9 bins at target 253 bins.\n", + "dict_keys(['z', 'proper time [Gyr]', 'conf. time [Mpc]', 'H [1/Mpc]', 'comov. dist.', 'ang.diam.dist.', 'lum. dist.', 'comov.snd.hrz.', '(.)rho_g', '(.)rho_b', '(.)rho_cdm', '(.)rho_ncdm[0]', '(.)p_ncdm[0]', 'q_ncdm[0][0]', 'w_ncdm[0][0]', 'q_ncdm[0][1]', 'w_ncdm[0][1]', 'q_ncdm[0][2]', 'w_ncdm[0][2]', 'q_ncdm[0][3]', 'w_ncdm[0][3]', 'q_ncdm[0][4]', 'w_ncdm[0][4]', 'q_ncdm[0][5]', 'w_ncdm[0][5]', '(.)rho_lambda', '(.)rho_ur', '(.)rho_crit', '(.)rho_tot', '(.)p_tot', '(.)p_tot_prime', 'gr.fac. D', 'gr.fac. f'])\n", + "Automatic quadrature has 5 bins at target 9 bins.\n", + "dict_keys(['z', 'proper time [Gyr]', 'conf. time [Mpc]', 'H [1/Mpc]', 'comov. dist.', 'ang.diam.dist.', 'lum. dist.', 'comov.snd.hrz.', '(.)rho_g', '(.)rho_b', '(.)rho_cdm', '(.)rho_ncdm[0]', '(.)p_ncdm[0]', 'q_ncdm[0][0]', 'w_ncdm[0][0]', 'q_ncdm[0][1]', 'w_ncdm[0][1]', 'q_ncdm[0][2]', 'w_ncdm[0][2]', 'q_ncdm[0][3]', 'w_ncdm[0][3]', 'q_ncdm[0][4]', 'w_ncdm[0][4]', 'q_ncdm[0][5]', 'w_ncdm[0][5]', 'q_ncdm[0][6]', 'w_ncdm[0][6]', 'q_ncdm[0][7]', 'w_ncdm[0][7]', 'q_ncdm[0][8]', 'w_ncdm[0][8]', '(.)rho_lambda', '(.)rho_ur', '(.)rho_crit', '(.)rho_tot', '(.)p_tot', '(.)p_tot_prime', 'gr.fac. D', 'gr.fac. f'])\n", + "dict_keys(['z', 'proper time [Gyr]', 'conf. time [Mpc]', 'H [1/Mpc]', 'comov. dist.', 'ang.diam.dist.', 'lum. dist.', 'comov.snd.hrz.', '(.)rho_g', '(.)rho_b', '(.)rho_cdm', '(.)rho_ncdm[0]', '(.)p_ncdm[0]', 'q_ncdm[0][0]', 'w_ncdm[0][0]', 'q_ncdm[0][1]', 'w_ncdm[0][1]', 'q_ncdm[0][2]', 'w_ncdm[0][2]', 'q_ncdm[0][3]', 'w_ncdm[0][3]', 'q_ncdm[0][4]', 'w_ncdm[0][4]', 'q_ncdm[0][5]', 'w_ncdm[0][5]', 'q_ncdm[0][6]', 'w_ncdm[0][6]', 'q_ncdm[0][7]', 'w_ncdm[0][7]', 'q_ncdm[0][8]', 'w_ncdm[0][8]', '(.)rho_lambda', '(.)rho_ur', '(.)rho_crit', '(.)rho_tot', '(.)p_tot', '(.)p_tot_prime', 'gr.fac. D', 'gr.fac. f'])\n", + "dict_keys(['z', 'proper time [Gyr]', 'conf. time [Mpc]', 'H [1/Mpc]', 'comov. dist.', 'ang.diam.dist.', 'lum. dist.', 'comov.snd.hrz.', '(.)rho_g', '(.)rho_b', '(.)rho_cdm', '(.)rho_ncdm[0]', '(.)p_ncdm[0]', 'q_ncdm[0][0]', 'w_ncdm[0][0]', 'q_ncdm[0][1]', 'w_ncdm[0][1]', 'q_ncdm[0][2]', 'w_ncdm[0][2]', 'q_ncdm[0][3]', 'w_ncdm[0][3]', 'q_ncdm[0][4]', 'w_ncdm[0][4]', 'q_ncdm[0][5]', 'w_ncdm[0][5]', 'q_ncdm[0][6]', 'w_ncdm[0][6]', 'q_ncdm[0][7]', 'w_ncdm[0][7]', 'q_ncdm[0][8]', 'w_ncdm[0][8]', '(.)rho_lambda', '(.)rho_ur', '(.)rho_crit', '(.)rho_tot', '(.)p_tot', '(.)p_tot_prime', 'gr.fac. D', 'gr.fac. f'])\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbYAAADQCAYAAABr00SDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAABJE0lEQVR4nO3deXxU1fn48c+TSSZ7AiQhgbAkyL4mJGwKyKKALOJWKy4VFSmt1q1Wq60/rdVvaW2rtVotKlK74M6OLAoIiAgBwr4jS9gTyL7PnN8fdxJDzIpJJpk879drXsnc9blzZ+5zz73nniPGGJRSSilP4eXuAJRSSqm6pIlNKaWUR9HEppRSyqNoYlNKKeVRNLEppZTyKJrYlFJKeRSPTmwiEiMiRkS8Xe/XiMi0eljPZyJyd10vt4L1lMYvIneIyIo6XPZuERnh+v85EflPHS77aRF5u66WV27ZL4hIqoicqY/lV7HeN0XkmYZcZ10p87vIFpHp7o5HuYeI3Of6DhgR6VwPy58jIi9UMT5bRDrV9XrBwxNbfajooG+Muc4Y86+GjMMY819jzJjqpqvuy1Vmeb2MMWt+aFwiMkJEUsot+/+MMfVxQtEB+CXQ0xgTVdfLL7OeqSKyvuwwY8wMY8zv62udDaSFMWZWyRsRGS0i+0QkV0RWi0jHymZ0JcfVrmn3icg1ZcZNFRGH68BV8hpRZnyciKwTkQwRSansBEFE/p/roHtNReMrmSdORLa44toiInFVTNtKROaJSI6IHBOR28uMGykiO0UkXUTSXNNFlxl/q4hscK1nTbnldhWRBSJyXkQuiMhyEelWi22oNK4KphUR+aMrxjTX/+IaFy4iX7mGp4vI1yJyVcm8xph3jDFBNY3LtcxJIrLLtU83iEjP2sxfljEmyBhz5HLnr4omtmaupDTbRHUA0owx59wdSFMnIuHAp8AzQCsgCfigilnmAtuAMOA3wMciElFm/NeuA1fJa02Zcf8D1rrWczXwcxG5vlw8VwA/Ak7XYhvswALgP0BL4F/AAtfwirwOFAKRwB3AGyLSyzVuDzDWGNMCaAscBN4oM+8F4BVgZgXLbQEsBLq5lr3JFVdNVRVXedOBG4B+QF9gEvBT17hs4F4gAuvz+COw6HJ/8yLSBfgvMANrGxcBCxvlMcQY0yReQHusH955IA14zTXcC/gtcAw4B7wHhLrGxQAG8Ha9XwNMK7PMe4G9wEVgOdCxzLhewEqsL/BZ4GlgHNYXrgjrS7O97HIBXyAd6F1mORFAHtDa9X4ikOyabgPQt4ptvhbYB2QArwFflsQPTAXWu/4X4GXX9mcCO4HeWF/6IlfM2cAi1/RHgSeBHUAB4O0ado1r/HPAx1gHtixgK9CvTFwG6Fzm/RzgBSDQta1O1/qysQ4KzwH/KTP99cBu12ewBuhRZtxR4HFXbBmuGPwq+GyuKbeuOcAIIKXcdOW360Os70iWK4bEqr5jQA8gH3C41pNedpvLzHs/cAjr+7IQaFvu85qBdXBMxzpwSQ2/92tcn+2Gkn2IlUz+69rXm4GYWv6WYijzu3ANmw5sKPO+ZF92r2D+rq7vTXCZYeuAGeW/m5WsPxerlF3y/iPgqXLTLAPGl91/NdiuMcDJsp8tcBwYV8G0gVi/i65lhv0bmFnBtL7AH4A9FYybBqypJq5Wrs87rAbbUOO4XOM2ANPLvL8P2FjBdF5YSc/gOhZV9nuuIrYHgSXllpkHjK5k+jnAm1jH0Sys41fHitbrmvZ1YIlr2m+AK1zjKjy+VRVrkyixiYgNWIyVvGKAaOB91+iprtdIoBMQhHVAqm6Zk7GS1U1YyWcd1lkoIhIMfI7142oLdAa+MMYsA/4P+MBYZ6H9yi7TGFOAdWCcUmbwrcCXxphzIhIPzMY6owoD/ol1xuNbQXwlZ9C/BcKBw8BV5adzGQMMxzrghLrWmWasy0z/Bf7kindSmXmmABOwLkcVV7DMyVgHnFZYZ9jzRcSnkvWXbH8OcB1wynx3pn6q3HZ1xfqcH8H63JdinUWWPau+FeskIhbrLHRqBev6vNy6vjdNJa7H+u60wEpAr7niqvA7ZozZi5WUSkogLcovUERGYR34bgXauJbxfrnJJgIDXNtzKzC2hvEC3Abc5YrpCuBr4F2sfbMXeLZMLDtcl50qev2jinX0AraXvHHty8Ou4RVNe8QYk1Vm2PZy08aLde/zgIg8U+6s/hXgJyLi47pENwTr91ayDT8CCowxS6uIt7Jt2GFcR0OXHZVsQ1eg2BhzoLJtEJEOIpKOdfB+HPhTLeMpMRw4Y4xJq8G01cZVziX7raJpRWQH1snZQuBt88OucEi5/wXrJLoydwC/xzqGJWMdjypzG/A7rNLlIeBF1/AKj29VBdkkEhswECvB/MoYk2OMyTfGlNzzuAP4qzHmiDEmG3gKuK0GxeMZwB+MMXtdB/b/A+Jc9xUmYn0R/+JaV5Yx5psaxvo/rB1U4nbXMLDOiv9pjPnGGOMw1n25AmBwBcsZD+w2xnxsjCnCOhhUVkGiCAgGumOdre41xlR3CedVY8wJY0xeJeO3lFn3XwG/SuKsrR9jnfWtdC37z4A/cGW52E4ZYy5glVDi6mC9JdYbY5YaYxxYZ8IlJydVfceqcwcw2xiz1XVy8xQwRERiykwz0xiTbow5Dqymdtv0rjHmsDEmA/gMOGyM+dz1vf0IiC+Z0BjT1xjTopLXz6tYRxBWCbmsDKzvVW2nXYt1sGsN3Ix1EvWrMtMuBm7BShj7gHeMMZuh9KTy/4CHq4i1rrYhs6ppjTHHXScy4VgnmPtqG5CItMMqiTxWw1mqjauC6TPKTRtUcp8NrO8EEIJ1LKrpd7oinwNXi3Uf3Y5VMLADAVXMs8QYs9b1u/gN1u+ifSXTzjPGbHJ9r//Ld7+RWh/fmkpiaw8cq6Rk0RbrDLnEMaxLa5HVLLMj8LeSs1msS0iCdVbcHuts9XKsBgJEZJDrwBYHzCuzzl+WPYt2rattBctpC5woeeM6Cz1RwXQYY1ZhlTxeB86JyCwRCakmzgqXVdF4Y4wTSKkkztq6ZH+5ln0C63MvUTaB52L9eOtK+WX7uU6CqvqOVaf8NmVjnVHW1TadLfN/XgXv6+LzycY6+JUVgnVZqFbTuk4yvzXGOI0xO4HnsRIZItIK60rI81gnS+2BsSJSknSfA/5tjDnqzm0oy3WCVXK/rsb3k1z3HFcA/zDGzK3hbLXZhoqmDwGyy5VacZ2ozQV+LSKXXGmqKWPMPuBurGPNaayEvwfr2FCZsseRbKzjbGXHkQp/I5dzfGsqie0E0KGSL9UprIRRogNQzKU//sqW+dNyZ7T+xpgNrnGVVUOtsjsEV0ngQ6yz1CnA4jKXbE4AL5ZbZ0AlX/rTWD96wKr9VPZ9Bet91RiTAPTEKrKXnCFXFm913TqUXbcX0A7rswbrS1f2LK1sjcTqlnvJ/iqzXSerma8mcsrG5bq8GFH55Jeo6jtW220KxLrUXBfbVCtiPbaRXcnrzSpm3c13pdeSbbjCNbyiaTu5Slcl+lUyLVifX0kJohPgMMa8Z4wpNsakYF22He8aPxp4SETOiPUIR3vgQxF5suotL42rb9nSCtal34riOgB4uypE1GQbvLFKoNWdMAIgIi2xktpCY8yL1U3/A+K6ZL9VMy2AD5Uf26rluorT2xgThnUJPAbrPm9lyh5HgrAun5+qfPJK11vZ8a1CTSWxbcI60M8UkUAR8StTbXUu8KiIxLo+uJJ7YNWdeb8JPFVS20hEQl3X9sG6VNJGRB4REV8RCRaRQa5xZ4EY18G+Mv/DuuR2B99dhgR4C5jhKs2Ja1smlDtAlFgC9BKRm1wH24e4NIGUEpEBrmX6YB3c87EqVZTEezlf5IQy634E65LpRte4ZOB2EbGJyDismm0lzgJhIhJayXI/BCaIVbXcB6u6fgHWTfAf6gBWCWyCa9m/xbrxXxNVfcfOAu2k8tp1c4F7xKpq7ov1HfymJqUO+e6ZspgaxlklYz22EVTJa0YVs84DeovIzSLiB/w/rPtV37v85rr/kww86/qcbsRKIJ+4tuk6EYl0/d8dq6ZlSa3AA9ZguV1EvEQkCuu3ssM1fjTWZcw41+sU1j3p113Le07KVa8vYw1WJZ+HXL/bB13DV1WwDTlY97Cfd+3vq7DuK//btZ6bRKSbK8YIrMvx21ylN1zffT+shOfl+hx8XONCsCqjfWWM+XX5dbsu5VV4slRdXBV4D3hMRKJFpC3W72mOaz2DRWSoiNhFxN91chCJVTGjQtV8vohIgmvbI4BZWIm7qku040tiwLrXttEYU93VovLrrOr4VqEmkdhcpaBJWJU4jmMVfX/sGj0ba6evBb7F2uhf1GCZ87Cqv74vIpnALqzKCLhKWNe61nkGqzbbSNesH7n+ponI1kqW/Q3WDmiLdU+kZHgSVu2517BqYh6igooRrmlTsao7z8S6rNUF+KqSzQnBSpoXsS6JpQEvuca9A/QU69Ln/Ermr8gCrM/4IlbFhZtc98TAuv8xCauG3x1A6XJdX/K5wBHXOi+57GCM2Q/cCfwdSHUtZ5IxprAWsVXIdQ/q58DbWKWlHKq+TFJ23qq+Y6uwzoLPiEhqBfN+jnXw/gQrOV7BpfdZq9Iea581eOmuLGPMeaz7YS9i7fNBlNkGsR5IL1viuw1IdE07E7jFtQywktMOEcnBqhz0KVayxxiTiVVh61HXvMlYv70XXOPTjDFnSl5Yieqi6zIWWJ9Xhb8D13foBuAnWN/Ne4EbSr5bYjUU8FmZWX6OdX/3HNZ39mfGmJLSTjTWJdMsrFp4TuDGMvPehXUZ+A1gmOv/t1zjbsSqKHRPuRJzhzLbUNWJXKVxicgwEckuM+0/se5D78T6HJe4hoF1Uvc61vHgJFapeIIpV6GrnEo/X5e/YX22+7H23/1VTAvWif2zWJcgE7B++7VV1fGtQlLuUqxSqgGJyG+B88aYf1Y7cd2tsyPWgSkfq7LMW9XM0miISDJW9fKa1DBslMRqhecjY8xyN8dxD1Y1ej+sxy+OeMLnC5rYlFJKeZgmcSlSKaWUqilNbEoppTyKJjallFIepfE1XtmAwsPDTUxMjLvDUEopj7Rly5ZUY0xNnyWtM806scXExJCUlOTuMJRSyiOJyLHqp6p7eilSKaWUR9HEppRSyqNoYlNKKeVRNLEppZTyKJrYlFJKeRRNbEoppTyKJjallFIeRRObUkopj6KJTSmllEfRxKaUUsqjaGJTSinlUTwqsYnIDSLyloh8ICJj3B2PUm53YhOs+4v1t4Zyt20j9Z+zyN22rR4DU6r+NJpGkEVkNjAROGeM6V1m+Djgb4ANeNsYM7OyZRhj5gPzRaQl8GdgRb0GrVRjU5QHqQch9QAcXgXb54Jxgtig+wSI6AZ+LcAvFPxbXPJ/cvYJ9n+9kj6//xQpdiDe3rR96U/4x8Ujdh+87HbEbgdvb0SkQTcrd9s2cjdtJmDgAALi4xt03arpaTSJDZgDvAa8VzJARGzA68C1QAqwWUQWYiW5P5Sb/15jzDnX/791zaeUZ8rPtJLX+X1wfv93/188Bhho0QFsdjDGmt44rXHOYshLJycrlcyMDIou5pN3wcm5LDvnsnzpehqk2DVLYSEnH37k++sWQXy8rZevHS+7L+Lri/j6W3/tdrxMPlKUiYRGIi2iELvdSo6u8eJjL532kuGlL9/SZFrw7VHOPPssprgYsdvp8O5sTW6qSo0msRlj1opITLnBA4FDxpgjACLyPjDZGPMHrNLdJcQ6jZwJfGaM2VrRekRkOjAdoEOHDnW3AUrVoeRzySSdTSIxtAtxTpuVvM7vh9T9cP4ABdmnyLLZyWrZgeyWHcgKbUtG2xs45x3CGeycz88lN/0woV4Z5DnAnmHDebCQ1ucPEJ1aQIdUB+GZNoq9AjkZBscjhBOdhC2dDT9eZ7A5wGGD+beG4xMVSLRDiCp00LqgiIjcXELzcpC8dEx+LsZhMA7BOAQHdoocPngX52GcgrmwH9NyLE7faExhIcXZOZjCQterAGdhIaawCFNQUDrcWVhgDSssxBQUgNNZ+rmYoiKr5KaJTVWh0SS2SkQDJ8q8TwEGVTH9L4BrgFAR6WyMebP8BMaYWcAsgMTERFOHsSp1WbYcvcCuvXsYEppGtPmWT099yV9yDuHA4AX0KCzGePuS4eVNlkBuuC/F4SUnZYVQeBjOnoQiO5EXfOhwDmIuOOlzoZiY1GDC0/MRAxnhkB3dnmPR4bzTJoQjQW05E9iG4d06ERmRxpLzz+KgmP3tIP5Qd7a3bsP2AEHy0vHySccWmI6EXgRxADb8pSMh3hGE28NoYw8lqNiPnfsymFC8iynOJIKwElKy7zySwtqRGJlIXPebIGYo+AbX+PPJTUri+H3TrBKbjw8BAwfU+T5QnkWMaTzHdleJbXHJPTYRuQUYZ4yZ5np/FzDIGPNgXawvMTHRaEejqsE4HXDxaJmS136OpexkZ+4JdvsLSX7+HLD7IEZwiAEBDBTndKY4uxeB3kGE+obQ70I6/Y8ex6tlOKF2PyLOnyHo5DG8U44jRYXYWrfGr0sXfLt0wbdrV+vvFZ3wCggAYMuxi9zx9kaKip34eHvx32mDSejY8rtSYmQica3jAHA4DRdzC0nLLiQtu4Bz2fmkZJzjROZJzuSeJjX/LBlF58h1plJAGuJ9EbHlA9DS4SDUYThu98Fg8EF463wW/XMzod0A6DQSOo2A6ASwVX2OrffYmiYR2WKMSWzw9TbyxDYEeM4YM9b1/ikA16XIH0wTm6oXxQWQdrg0eZW8TNohTkkxm4LC+Mo/lK02OO9TgK3YF3I7kJPXlW6hfUmICWbh2d/hNEV4e/nw52H/YIh/DPkbNpCxcCG5X39duip7t24E9I/Ht0uX0mRma9Gi2hC3HLvIxiNpDO4URkLHlnWy2VuOXeSOtzZSaHLpZd/OgJaHWeGdQXrAMUrqmniLnaEt+zBafBl6/ijhJ7aAjz/EDIMrRlrJLuwKuMzKKacOpXP6UDrRXVsS1Sm0TrZLXT5NbFSY2LyBA8Bo4CSwGbjdGLO7LtaniU1dthOb4NAXENoObD6X3gO78C0YB86QdhwM68iX3gGsL85nrzONfK9sKGpJmE934sL7071lX175LJ2iYnNp6enUFvavX0TPgwUEbTlA/p49eLdpg3d4OPm7dlmVQmw2Ih56iPCfTnf3p1GqfMJMPpfMfcunUeQsQrDhkzmOLMc5fIMPYLwvEmmPZURwe8Y5iohL2Y73+X0Q0g6uGPFdiS4wvEbrProzlSWv7wABb28vJj8ar8nNzZp9YhORucAIIBw4CzxrjHlHRMYDr2DVhJxtjHmxrtapiU3V2vn98PU/YOu/ANdvJzga2vajKKwzO/1DWJ6VwVcZx0kpOIBDcjCFkUR49yC+dTwTu17F8E5d8LZ99whpSTIYEuKk07c7yF63npwNGzD5+QQMSCRw2HCChg3FfsUV5CUnc/yeezFFRYiPT5OoIVj+Euep9Dw2fZvGqsO7SDq/gYvOXXgHHsEmdjr79WCsfyuGXjxLh5MbCSxMg6g+VpK7YiR0GGKV8E5sgqPrrJJe+4EAJH12lG8WHLFWKjB4cicSxsW4b8OVJjZ30MSmqmUMnN0NexbA3oVWlfqQaMg8RZ7Adj8/vowZwWpj41T+PpwUYwqiae3Tg4TI/kzqdiVXxnbEx/b9thBMYSG5W7eRs34d2WvXUXDgAD7t2xM0bBiBw4YSOGhQ6X2xsjztftPFnEK+OnKKzw59RXLq16SzEy/7BRz5UbTIbctTMWGMyzmOz7GvSPb2Iim8HYnnjxFXUAA2X7h7IbQfyJkjGSx4eRvFRVallaE/6kK/0e3dvHXNmyY2N9DEpipkDJzaZiWyPQusS4sdhpDRbQxJLdux+sh6jpxaxl67DzYDtvy2BNgTSIjsz8RugxkS2xZfb1vp4somIp/Wrclet57s9evI3fA1xuEgYNBAgoYOI2j4MHw6dmzwh58bm5dX7ue1dRvxCtyPd9B+bAFHsNt86N6yG3sv7sTpdGA3Tt46c464QgeM+g0M+yUAZ45kcPLARbIu5LNvwxnG/bQ3MX1qdilT1T13JbbGXt1fqYbhdELKZlcyWwiZKZyNGUxSj2tY7YSktL2kHZoNDn8ceTH0so3iEaeT3t0m03PAWPzttgoXm5uUxPF778UUFlkVIozBHhtL4LChRP/tbwQMSMTLz6+BN7ZxG961Nf9c25qi9HBM1jBuv6otK498TXLuYmz+DhAoQFgfEECco8C6HOkS1Sm09L5ay8hAPvvnTsbe15tO8RHu2hzlBlpi0xJb8+V0wPGvYc8CzN5FHCu4QFK7vmwIaMHm/AukO85DcSjFuTFE2nswqE0iY7r2Y1BsOEG+lZ8TGoeD3M2byVz6GRkLF2LyrerviNDqnqlEPvFEA21g01W+Eooxhvd3rGNm8qM4TBGCwS427u5wHVOveooQe0iFy9n1ZQrrPjzI6Kk96DogqoG3QumlSDfQxNb8JJ9JImnfJ/S/cBq/Y1+xRQrZ3CKaJK9iskweFEVQmNORSJ+eDG6byDVdejCoUxih/j5VLtc4HORu2ULWsmVkLl+BMyeHoKuvxrd7d9LefLP04eKmUNmjMUs+l8zqY19z5kxbVhzYT3HICuy+2fy42538PH4qQfag782zd8Mp1vxnPyPu7E6PK9u4IermSxObG2hiayaKC+HIGr7a/i4P5uzEagrRuo9lK2pDXnYsrX16MCQ6kVFdrmBQbBgtA+3VLtY4neRt20bmZ8vIWr4cR0YGgcOHEXLddQSPGIFXYCDgeZU9GouCYgdLdpzkH5s/4bQsxO6bx4SOt/HrK6d9L8Ed2HyGL+bsZdiPu9J7eLSbIm5+NLG5gSY2D1aUB4dXcWLT+ySdX8/KADvr/e1WBX0BjNDF90am9vwpgzu1onVwze5zGWPIS062SmbLllN84QJBQ4cSct04gkaNwhb0/RKDqn/bjqfyp/UfsCPnI2ze+fQNnky/4PGM6t6h9AH0w9vOseLt3Vx5U2etLdlANLG5gSY2D1OYAwdXkL97PmtPrGGpvy9f+vshTjv5mXF0atGVk7YPcBoHdpsPb415q7TpqKoYY8jftYvMpZ+RuXwZxefOE3jVlYSMu47g0aOwhVR8f0c1vLOZOTz+2Ry2Zn4IXgU4Ll7N7Bsf4sorrFLa0R2pLJu1iwETY/QZtwagtSKVuhz5mXBgOUW75/H1yfUsCfDnC38/iluGY8uPJ+tEH4pzO2ETL66L68ZVvcZ8rz3EiuRu3UbGksU4s7LI27qNotOnCRw8mIgHHiB49OgaNVulGl5kSCCDW4/lq+T2eIVsxTd8FT9bu46fZU3j7t5TiOkbzoSf92XpGztwFDkZMDG22T9e4Yk0sammJ/cC7P8Mx54FJJ38iiVBwSz39yM/LAy/ojhGRY7lpwOuIz3XyR1vb8SI1djv4E5hxLXuXGVCKzx+nNRZs8j45FPreTYvL1rdey9h992Ld8u6aVNR1a/BncKwe/tQlDmA4pwEQtvs5B9b5vDvPf/i+s7jCbGH0PPORJL/d4LiIidDbrxCk5uH0cSmmoacVNi3GLN7PttPbWRJSEuW+vmTFRGBd2Evrmx1LdMTJ9I3+tKHcf87bXC1jf0WX7hA5mefkbloMXnJyXhHlakWLoItOFiTWhOS0LHlJfu9X7sJvLrqBmbtfJ3/FP0HALuXnZfvepPd/zmJo8jJ0Fu7aHLzIJrYVOOVedpKZnvms//UZha1iGCxnz8Xolrjld+dhOBR3JcwkStjoys9KCV0bFlhQnPm5ZG1ahWZCxeR/dVX+LRtS+ikSbSd+QeKL168pD1G7f+r6Sm/3x+7tic5ft358PA6EEOhs5DP8j7l0UefYuHfkikudjJiSjfES5ObJ9DEphqX9BOwdxHsWcCRM1tY3LItC/z8ONc2EvI709v/ap6Om8Q13Tphq+VByDgc5GzcSObCRWStXIn4+hIyfjwx//k3fv36lSZHe0wMHd6drVX0PcykbsNYcOxfFDiKMAYWH1lKbEgnbnz0Nha9moyzyMnIn/TAS5Nbk6e1IrVWpPtdOGI1Y7V3IafOJrMorCPz/Pw46ZWNM68jXYOGcWfvSUzs3e2SNhhrwhhD/p49ZC5cRObSpTiysggePZrQ6ycReOWViE/VD14rz1LS04BPUWf+vnoXRWHvMTz6an7T82k+e3U3bbu04Jp7emKroNFqVXta3d8NNLG50fkDsHcB7FnA+fO7WdL6Cj61+/OtVzrO/LZ09L2KW3tO4Ef9+hFYRfNV5ZU8DG2PjaHwyBEyFi6i8OhRAocMIfT6SQSNvgZbUGA9bphqKtJzC3ls/ko25bxEdEgEbwx+hXVvHiWiQzBjp/XG5qPJ7YfSxOYGmtgakDFwbo/VWv6ehaSn7WdZVDc+9vVnv0nDWRhOlG0wN3SdwF0JA2gRUH3LH+VlrV1Lys8fgGKrbRF7p060/PGthIwfj3eENoKrKvb+ln28uPnX+Phd5M/9X+Hkh3m0iAzkup/2xruSxq1VzWhicwNNbPXMGDidbF1m3LOA7IvfsrJdLz70CWC38yyOolBamYFMuGI8UxOHEBXqX+tVOPPzyV6zhoyFi8hes8ZqpR/Ay4uIhx9uVL1Lq8brxMUs7pr/FGl8zdigXxK7qR0hYX7c9mgCPr6a3C6XJjY30MRWD5xOOLkF9syHvQvJz0xhdbs4PvAJYJvjFMXF/gQVJ3Bth7HcN2AEsRG1b4KqpPX8jEWLyFq+ArHZCB5/HX7dunN25swm1bu0ajwcDif3fPoyW3Pew/v0BG5KGUVM22BueywBu7/Ws7sc2vKIarqcDji+sbQvs6Kcc6zv2J/3W1/BNy18KXbmYs/qyqh29zItcRQ927S4rGeG8vfvJ2PhQjIXL8GRnk7QqJG0/dOfCBp6FWK3Ll36duuqtRnVZbHZvBgUPpmv93lB9P/41Pc8087dxcJXk5n0i374BmhFo6ZCS2xaYrs8jmI4tt66Z7Z3MY78dL6JGcj7PoGsLzhGodOBLa8Pg1uP5t6EsQyMCb+sZFZ0+jSZS5aQsXARBQcPEjBoEKGTJhE85lpswcH1sGGqOdty7CJ3vL2RIq8z+LX7F1E+bZh26iG8jI3rH47DP6j2936bM70U6Qaa2GqpuBC+/dJKZvuW4CzKZXunIcy1B7Mq5wgFzlzI7Ulcq5Hc0/86hndug/dlVJt2ZGaStWIFGQsXkbt5M77duhE6aRIhEyfgExlZDxum1HdKOjkVrxxe3/McYQG5TD/1DLYiHyY/Ek9AiCa3mtLE5gaa2GqgKB8Or7KS2f7PMM5i9nQeygf2YJZnHSbXmY4ztys9Q67mzr7XMa5nx1o9a1baV1n/eByZmVYlkNWrsYWFETpxIiGTJuLXtWs9bqBSlVuyM4XHV/2O4JDd/OzM7zE5NhwTvmXAFfE16hmiudPEVgdEJBD4EnjOGLO4uuk1sZVxYhMcXQcxwyCyFxxcad0zO7AcxMaRzlfzgV8IizMOkFl8FkfuFXTyH8qU3uOZ3LcLQbV41qxE7tatHJ96D6awEAAJCCB0wnhCJk0iIDER8dLniJT7zduWwlMr38Q/fBljD9xLaF44y/u8xSvXv6TJrRrNuvKIiMwGJgLnjDG9ywwfB/wNsAFvG2NmVrOoJ4EP6y1QT3ViE/xrEhQXgAh4eYM9kJQuo/i4/23MT9tHWm4yjtSORPuM4Cc9JvCjuJ60qkEv0+UZYyjYu5eMJUtI//Cj0qSGlxdh0+4j4uc/r+ONU+qHuTG+HdkF9/N/m07yWbdZXHvgHsbt+CmbuyVrYmukGkViA+YArwHvlQwQERvwOnAtkAJsFpGFWEnuD+XmvxfoB+wBatYVsvrO0XVWUsNwzsuLxVcMYa6XN2eyNuM4F02EDOKerv+P2/v3Iyr08j7ewqNHyViyhMwlS62WQAYPouXtU7gw51+Y4mLEx4fAIUPqdruUqiN3De7I4Yxb+PjUNlZ2nc21h6YStjCOi11yaBmlLdk0No0isRlj1opITLnBA4FDxpgjACLyPjDZGPMHrNLdJURkBBAI9ATyRGSpMcZZn3F7iott41gZGspn/r5s8bPjk30Rn+Kr+FHsw9w1IJHY8Mv74RadOWP1Or1kCfm7d+Pfrx8tb7uNkOvGlbYEEjRihFbPV03Cc2MnkLvIwWcn32ZDj08ZkXc18/6ylcmPxBMWXfvnMVX9aTT32FyJbXHJpUgRuQUYZ4yZ5np/FzDIGPNgNcuZCqRWdo9NRKYD0wE6dOiQcOzYsTrbhqYkqzCLZUdWMnfPIg5mbcWnKJA+eS3oG3EdY4fdRa+2oZdVPb/44kWylq8gc8kScpOS8O3cmZCJEwmZMB57u3b1sCVKNRxjDM8u3Mn8M8/SIyqI6ZlP8W1yKpMfjieigz5+Ul6zvsdWl4wxc6oZPwuYBVblkYaIqbHIK87ji2Nr+N+uhexO/wZHsT+23Hiubvs77u4/jISOrS6ryw5nTg5Zq1aRsXgxOV9twCcqipAJE4h85rdao1F5FBHhuUl9SP/4Udak/obVnecxxvsW5r+8jStvuoL8nCKiu7YkqlOou0Nt1hpzYjsJtC/zvp1rmKqFQkch61O+4j+75rM19SscDhvk9CMx/Al+MngkwzpHXNazZs7CQnLWriVjyRKyV6/BKyiIkOuuI+LnP8evb1/tjVh5LC8v4eVbhjH1fw+x4NsX6ZHYm465nVnz3/0g4O3txeRH4zW5uVGtEpuIHAJ2AduBHcB2Y8yh+ggM2Ax0EZFYrIR2G3B7Pa3LoxQ7i9l0ehP/3b2Ar898SZHDgSO7F31bPMBd/a9hdPe2+PnUvmFX43CQ+803ZCxZQtaKlSBC8Jhraf/GPwgYOBCxaWOxqnmweQnv3HYzN/77W2YmPc/vg/5hjTBWm5MnD1zUxOZGtS2xzQfswBms2or/FpFUrMSzwxjz08sJQkTmAiOAcBFJAZ41xrwjIg8Cy7FqQs42xuy+nOV7uuRzyWw+s5lQeyibT+9hTcrn5DtyKc7qQbegqUzpM4bxvdsT7Fe7tu6sh6c3YWvRgoKDh8hctgxndjbBo0bSduYfCBw2DC+7tsKgmie7txcf3/EYY/59mDfOvM5k73swxQYQoru2dHd4zVptE9sYY0zfkjci8i5wI1ZV/X6XG4QxZkolw5cCSy93uc1B8rlkpi6bisM4wEBxbkc6+PyIW3uN5Ya+sYQF+dZ6mcYY0j/5lDPPPgsOBwD+8fFEPvkkwaNG4hWo1ZuVAvC3ezPvtr9w7fu38lGHefQ7MZH2RTaOFxcR5e7gmrHaJrY0EelnjNkOYIz5RkRmGWN+jfWsmapnJe3YdY8K5sSFXN7d8wnFPg5EQMSL+xMn8NjAn9V6ucYY8nftInPZMrKWr6Do5EmrPzUALy+CRowgdNL3nrJQqtkLDwxicpun+JjHWYYvdxwdy+blRxnYNdzdoTVbtU1s04F/ichuIBnoAeTVdVCqYluOXWTKrI0UOqzH89qG+nF134GsuPAZDlOMj5cPo2Jq/pCzMYb8HTvIXLacrOXLKTp7lsBBgwibfj/eUVGcfOjh0r7NAgYOqK/NUqrJm9SrNx9sn4J3m3f5JrMnEw54k5NRQGBo7a+YqB+uRolNRMYDycaYgyIyFLgBiAcOAc/WX3iqrI1H0ih29RDtJXDH4A48MLILPz7XnqSzSSRGJlbbxI8xhvzt28lctpzMFcspPneewMGDCf/5zwgaPRrvlt/dG+jw7mx9eFqpGkjo2JL/3H43r2/LY2/3dwg9/ieSlh7l6ind3B1as1TTEtuNwPMiEgnsw6oVmez6m14vkanvGdwpDLu3F0XFTny8vRjcybrUEdc6rsqEZpxO8pK3k7V8GZkrVlJ8/jyBQ4YQ8cADBI0adUkyKysgPl4TmlI1lNCxJW93eJKHV51ijeND4tZNpt+o9rSIDHB3aM1OrVoeEZGngWhgN5AATAEuGGOaZJMSTbF1/5J7bIM7hZHQsfKaV1YyS7buma1YSXFaGoFDBhMydhzBo0dha9Gi4YJWqhnJLMxkyuIpjN49lZ6RPRg3vXf1M3moptLyyI+NMaW1H0XkH8Cv6jYkVZWEji0rTWjG6SRv61Yyl68ga8UKii9cIOjKK4l4+GGCR43EFuq5z9UUFRWRkpJCfn6+u0NRP4Cfnx/t2rXDx6d2j6Y0JiH2EF4e+TIPpv4S322PcPbbDkTGhrg7rGaltoktU0QSjDFbAIwxW0RE20xyI+NwWMls2XKyVqzAkZ5O4FVXEfHoIwSPGoUtpHn8oFJSUggODiYmJkZbPWmijDGkpaWRkpJCbGysu8P5Qbq27Moj1/yM5SlJrPzQzh1PDNXvZQOqbWK7D/hURDYDW4A+QFGdR6WqZBwOcrdsIWvZcjJXrsCZnkHg0KG0fvyXBI0ahS24+TXGmp+fr0mtiRMRwsLCOH/+vLtDqRPjO41n54h9XPg4j93bjtG7f4y7Q2o2apXYjDEHRKQ/Vq3IPsBe4Ol6iEuVYxwOcjcnkbl8GVkrP8eZmUngsGFEPvEEQSNGNMtkVp4mtabP0/bhYyN+we+2vsGy97Po2jcau3fTvcTalNS6EWRjTCFWL9XaU3U9M8XF5CYlWRVAVn6OMyuLwOHDiHzySYJGjsAWpH1AKdWY+Xj5MOPuH/HR77byx/++QZv+ATV6LEf9MLVv1l3Vi9xt20j95yxyk5LI2bCB0//vWQ4Ov5oTP52BIy2NyKeeosuGDbR/7TVCJ03UpNZIzZ8/HxFh3759NZr+lVdeITc397LXt337doYMGUKfPn2YNGkSmZmZ35tmzZo1TJxYcasx48ePJz09/bLXr6rXLqINLQc78dnSlteT/sH9K+4n+Vyyu8PyaJrYGoHcbds4PvUezr/8MsfuvIvjM36GIz2dqN/+hq4bvqLd3/9O6MQJ2IK0jcbGbu7cuQwdOpS5c+fWaPofmtimTZvGzJkz2blzJzfeeCMvvfRSreZfunQpLfTRj3pX1OscgtDzzFUUOYtIOtu0HjNqai4rsYlIoIhoHyV1JHfTZkyRqw6OlxB+//20e/VvhIwfrw0O16Mtxy7y+upDbDl2sU6Wl52dzfr163nnnXd4//33S4eXLzE9+OCDzJkzh1dffZVTp04xcuRIRo4cCViJsU+fPvTu3Zsnn3yy2nUeOHCA4cOHA3DttdfyySefVDhdZmYmEyZMoFu3bsyYMQOnqwWbmJgYUlNTOXr0KD169OD++++nV69ejBkzhrw8q7W8V199lZ49e9K3b19uu+22y/twmrkB0Qkkd/yc/ifHEOgMJjGywR/talZqlNhExEtEbheRJSJyDqv1kdMiskdEXhKRzvUbpmcLGDgAsdvBZkPsvgQOvcrdITVpDqfhbGZ+la/P95zl9rc28pcV+7n9rY18vudstfM4nFU3ZrBgwQLGjRtH165dCQsLY8uWLVVO/9BDD9G2bVtWr17N6tWrOXXqFE8++SSrVq0iOTmZzZs3M3/+/CqX0atXLxYsWADARx99xIkTJyqcbtOmTfz9739nz549HD58mE8//fR70xw8eJAHHniA3bt306JFi9IkOXPmTLZt28aOHTt48803q4xHVSyudRxP3DEd/xAffuXzR73HVs9qWnlkNfA58BSwyxjjBBCRVsBI4I8iMs8Y85/6CdOzBcTHa7uMdSg1u4BB//dFjacvKHYy7b3qLw198/RoIkP8Kh0/d+5cHn74YQBuu+025s6dS0JCQo3j2Lx5MyNGjCAiIgKAO+64g7Vr13LDDTdUOs/s2bN56KGH+P3vf8/111+PvZL+8QYOHEinTp0AmDJlCuvXr+eWW265ZJrY2Fji4uIASEhI4OjRowD07duXO+64gxtuuKHKWFTV4qPiaXFbKsvf2kX2+AKCWmoDyfWlpontGmPM955XM8ZcAD4BPhERrcf6A2i7jHUnPMiXb54eXeU0O1MyeOB/Wyl2OPG2efH67f3p067qllnCq+jb7sKFC6xatYqdO3ciIjgcDkSEl156CW9v79JLf0Cdto7SvXt3VqxYAViXJZcsWVLhdOWr0VdUrd7X97vts9lspZcilyxZwtq1a1m0aBEvvvgiO3fuxNu71hWqFRDTJ4zWHYPZvPgII+/q4e5wPFaNLkWWJDUROSQi80XkdyJyc9lLkBUlPqXcweYlRIb4Vfm6pmck/7t/MI+N6cb/7h/MNT0jq53H5lX5M1Yff/wxd911F8eOHePo0aOcOHGC2NhY1q1bR8eOHdmzZw8FBQWkp6fzxRfflSaDg4PJysoCrFLVl19+SWpqKg6Hg7lz53L11VcD8JOf/IRNmzZ9b73nzp0DwOl08sILLzBjxowK49u0aRPffvstTqeTDz74gKFDh9bos3Q6nZw4cYKRI0fyxz/+kYyMDLKzs2s0r/o+EeHKmzqz9+szXDiV4+5wPFZtK4/MB44DZ4BrgR0iclxEvhaRf9Z1cErVp4SOLXlgZOcqG5Ouqblz53LjjTdeMuzmm29m7ty5tG/fnltvvZXevXtz6623El+mZD59+nTGjRvHyJEjadOmDTNnzmTkyJH069ePhIQEJk+eDMCOHTto27Zthevt2rUr3bt3p23bttxzzz0VxjdgwAAefPBBevToQWxs7PdirYzD4eDOO++kT58+xMfH89BDD2ktyh8oqlMosX3D2bjgsLtD8Vi1bd1/hzGmb5n3g7C6tHkN6GeMqfg6SCPVFFv3VxXbu3cvPXp45qWdzMxM7rvvPj766CN3h9IgPHlflrh4Joe5z2/ixsfiadO5hbvDqTfuat2/tiW2NBEpbd3fGPMNcJ0xJqWpJTWlmoqQkJBmk9Sai5ZRgfS4qg0bPj1MbQoXqmZqm9juB94QkbdE5AEReQ3Iq4e4lFLKow2cEEvqiSy+3Z7q7lA8Tq0SmzHmEDAU+AyIAg4BE+ohLqWU8miBLXzpd017Ns4/jNPhrH4GVWM1fUC7tDqYMcZpjPnUGPOMMeYVY0xa+WncwfUQ+Ysi8ncRududsSilVE3Ej+lIXlYR+74+4+5QPEpNS2yrReQXItKh7EARsYvIKBH5F3DZyUREZovIORHZVW74OBHZ73rM4NfVLGYy0A6rf7iUy41FKaUaiq+/N4njY9i06AhFhQ53h+MxaprYxgEOYK6InHI1pfUtcBCYArxijJnzA+KY41pHKVdblK8D1wE9gSki0lNE+ojI4nKv1kA3YIMx5jHgZz8gFqWUajC9h0fj5e3FjlUVN4emaq+mie2Pxph/YD271hEYDcQbYzoaY+43xmz7IUEYY9YCF8oNHggcMsYccfUB9z4w2Riz0xgzsdzrHFYpraQ120pPfURkuogkiUiSp/TUqxqPhu625rnnniM6Opq4uDji4uJYunTp96bRbmsaN5uPF4Ou78TW5cfJz9Z2LupCTRPbcNffdcaYImPMaWNMej3FVCIaKHsKk+IaVplPgbEi8ndgbWUTGWNmGWMSjTGJJW3yKVVXGrrbGoBHH32U5ORkkpOTGT9+fK3m1W5rGoeuAyIJDvMjadlRd4fiEWqa2L4Qka+BKBG5V0QSRKRRteBpjMk1xtxnjPmFMeZ1d8ejmoATm2DdX6y/dcAd3dbUlHZb07iJl3DljVewc00KmWn6BNUPVdO2Ih8H7sS6xBcLPAPsEpHdIvJBPcV2Emhf5n071zClquZ0QObpql/7lsK/JsGqF6y/+5ZWP4+z6pv77ui2BuC1116jb9++3HvvvVy8WHHfctptTePXvmcr2nZuwaZF37o7lCavxk10G2MOi8g1xpgDJcNEJAjoXS+RwWagi4jEYiW024Db62ldypNkn4O/dq/59MX58P6U6qd7bB+EtKl0tDu6rfnZz37GM888g4jwzDPP8Mtf/pLZs2d/bzrttqbxExGG3HgFH/0hCbu/ja4DoojqVHWPE6pitep7omxSc73PBjb+0CBEZC4wAggXkRTgWWPMOyLyILAcsAGzjTG7f+i6VDMQ1NpKQlU5tQ0+vgccRWDzgVvehbbVdBsU1LrSUe7qtiYyMrL0//vvv7/SSiLabU3T4HQYRGDn6pPsXX+ayY/Ga3K7DNVeihSRa11NaMW53k+v6yCMMVOMMW2MMT7GmHbGmHdcw5caY7oaY64wxrxY1+tVHsrLZpWsqnp1Hw93L4JRv7H+dh9f/TxetkpX6a5ua06fPl36/7x58+jdu+ILKNptTdNw8sBFSlqOdBQ7OXmg4kvLqmo1Oe26F+u5sN+6esyOq9eIlGoo7Qdarzowd+7c71X2KOm25o033ijttiY2NrbCbmtK7rWVdFtjjGHChAnVdlvzxBNPkJycjIgQExPDP/9Zce9RJd3WHDp0iJEjR9a625qMjAyMMdptTT2L7toSb28vioucIEJ01x/epVJzVG23NSIyyxgz3fX/TGC0MWZAQwRX37TbGs/hyV2daLc1zcuZIxlsW3mcUwcucs9Lw/CqooPbxq4xd1tT2h2NMebXwHv1F45SqjzttqZ5ieoUyui7e1Bc6OTkPr0UeTmqvRRpjFlQ8r+IHMKq5h8O7AC2u1r8V0opVUfsft7ExkWw75vTtO/Zyt3hNDm17Y9tPnAcOIPVvNYOETkuIl+LSMUX95VSStVa98FRHNl2nsL8YneH0uTUts7uGGNM35I3IvIucCPwGtCv0rmUUkrVSrvuLbH7e3Nk23m6D6n8+Un1fbUtsaWJSGkCM8Z8A1xnjEkxxiypYj6llFK14GXzouvAKPZt1L7aaqu2Jbb7gfdEZDeQDPQAtGEzpZSqB90HR5H8+XGyLuQT3MrP3eE0GbUqsbkqigwFPgOigEPAhHqIS6km6ezZs9x+++106tSJhIQEhgwZwrx58xo0hqlTp/Lxxx836DpV/QiLDiK8XRAHNmmprTZqeykSY4zTGPOpMeYZY8wrxpi0+ghMqabGGMMNN9zA8OHDOXLkCFu2bOH9998nJaX5deheXFxc5XtVc90GRbF/4xmqe+ZYfafWiU0pT5F8Lpm3d75N8rnkOlneqlWrsNvtzJgxo3RYx44d+cUvfsHRo0cZNmwY/fv3p3///mzYsAGovEsbgF//+tel3cU8/vjjAHz00Uf07t2bfv36MXz4cGoqOzub0aNH079/f/r06cOCBaVP8fD73/+ebt26MXToUKZMmcKf//xnAEaMGEFJAwapqanExMQAVmskv/rVrxgwYAB9+/Ytbe1kzZo1DBs2jOuvv56ePXt+731l86mqdRkQSfq5PM4dy3J3KE2GtmSqPI7D6SAtv+oLCbvTdvOrL39FkaMIH5sPL139Er3CelU5T5hfGLYq2ovcvXs3/fv3r3Bc69atWblyJX5+fhw8eJApU6ZQVas3aWlpzJs3j3379iEipb1cP//88yxfvpzo6Oha9Xzt5+fHvHnzCAkJITU1lcGDB3P99deTlJTEJ598wvbt2ykqKqJ///7V9kjwzjvvEBoayubNmykoKOCqq65izJgxAGzdupVdu3YRGxvLmjVrLnk/a9asCueLjY2t8XY0R4GhvnTo2Yr9G88QGRPi7nCaBE1syuOk5acx+qPRNZ6+wFHAQ6seqna6L370Ba0DKm/hv7wHHniA9evXY7fb+fzzz3nwwQdJTk7GZrNx4MCBKucNDQ3Fz8+P++67j4kTJ5aW6q666iqmTp3Krbfeyk033VTjWIwxPP3006xduxYvLy9OnjzJ2bNn+eqrr5g8eTJ+fn74+fkxadKkape1YsUKduzYUXofLyMjg4MHD2K32xk4cOAliars+8rm08RWvW6Do1j7/gGuuqUzNm+90FYdTWzK44T5hfHFj76ocprLLbFVpVevXqWdcwK8/vrrpKamkpiYyMsvv0xkZCTbt2/H6XTi52fVcKusSxtvb282bdrEF198wccff8xrr73GqlWrePPNN/nmm29YsmQJCQkJbNmyhccff5xt27bRtm1bli5dWmFs//3vfzl//jxbtmzBx8eHmJiYarvPKRtb2WmNMfz9739n7Nixl0y/Zs0aAgMDLxlW9n1l86nqxfYNZ81/9nF8dxqx/SLcHU6jp6lfeRybl43WAa2rfI1sP5K3x7zNL/r/grfHvM3I9iOrnaeqy5AAo0aNIj8/nzfeeKN0WG5uLmCVTtq0aYOXlxf//ve/cTis3rgr69ImOzubjIwMxo8fz8svv8z27dsBOHz4MIMGDeL5558nIiKCEydO8O6775KcnFxpUitZf+vWrfHx8WH16tUcO3YMsEqAixYtIj8/n+zsbBYvXlw6T0xMTGkv4GVrWY4dO5Y33niDoqIiAA4cOEBOTk61++Vy51PgbbfROaE1+/WZthrREptqtuJaxxHXOq7OlicizJ8/n0cffZQ//elPREREEBgYyB//+Ef69+/PzTffzHvvvce4ceNKSzLt27evsEubrKwsJk+eTH5+PsYY/vrXvwLwq1/9ioMHD2KMYfTo0fTrV3GDPz/96U955JFHStexaNEiJk2aRJ8+fUhMTKR7d6uH8QEDBnD99dfTt29fIiMj6dOnD6GhVseWjz/+OLfeeiuzZs1iwoTvnuqZNm0aR48epX///hhjiIiIYP78+dV+Ppc7n7J0G9yGBX/bRn5OEX6BPu4Op1GrttsaT6bd1niO5t7VyQ+RnZ1NUFAQubm5DB8+nFmzZlVaCaYh6L6smHEa/vP/viZ+TEd6D492dzg10pi7rVFKebDp06cTFxdXWqp0Z1JTlRMvoeugKPZvPF39xM2cXopUqpn73//+5+4QVA11GxRF0pKjpJ/NpUVkgLvDabS0xKaUUk1Ei9YBRHUKYf83WomkKprYlFKqCek2uA37vzmDcTbf+hHV8ZjEJiIdRGS+iMwWkV+7Ox6llKoPnRNak5NRwOnD6e4OpdFqFInNlYzOiciucsPHich+ETlUg2TVB/jYGHMvEF9vwSqllBv5BfoQ2ydc+2mrQqNIbMAcYFzZASJiA14HrgN6AlNEpKeI9BGRxeVerYGNwH0isgpY1sDxq2YuLS2NuLg44uLiiIqKIjo6uvR9YWFhva132rRp7Nmzp8ppzp8/z6BBg4iPj2fdunX1FotqON0GR3FoyzmKCx3uDqVRahS1Io0xa0UkptzggcAhY8wRABF5H5hsjPkDMLHctIjI48CzrmV9DLxbz2ErVSosLIzk5GQAnnvuOYKCgkpb5Aer2xZv77r/ub399tvVTvPFF1/Qp0+fGk2rmoYOvcKweXvx7fZUugyIdHc4jU5jKbFVJBo4UeZ9imtYZZYBD4nIm8DRyiYSkekikiQiSefPn6+TQFXTlLttG6n/nEXutm31svypU6cyY8YMBg0axBNPPMGmTZsYMmQI8fHxXHnllezfvx+AOXPmMHnyZEaMGEGXLl343e9+B8Cbb75ZWuqLjY1l5MiR31tH2a5lgoKC+M1vfkO/fv0YPHgwZ8+eJTk5mSeeeIIFCxYQFxdHXp52eO8JbN5edBkQqZcjK9EoSmx1wRizC7ilBtPNAmaB1fJIfcelGp5xOChOrbrbmrxduzj12GOYoiLEx4e2f/0r/r17VzmPd3gYYqu6vcjyUlJS2LBhAzabjczMTNatW4e3tzeff/45Tz/9dGmjyZs2bWLXrl0EBAQwYMAAJkyYwIwZM5gxYwZFRUWMGjWKxx57rMp15eTkMHjwYF588UWeeOIJ3nrrLX7729/y/PPPk5SUxGuvvVar2FXj1n1wFB/PTCIno4DAUF93h9OoNObEdhJoX+Z9O9cwpapUnJrGoauvrvH0pqCAkw88UO10nb/8Ep/ImndbA/CjH/0ImysZZmRkcPfdd3Pw4EFEpLQxYIBrr72WsDCr94CbbrqJ9evXk5hotUT08MMPM2rUqGq7lLHb7aXd2yQkJLBy5cpaxaqalogOwbSIDODg5rPEXdPB3eE0Ko05sW0GuohILFZCuw243b0hqabAOzyMzl9+WeU0l1tiq62y3bY888wzjBw5knnz5nH06FFGjBhROk5ELpmv5P2cOXM4duxYjUpbPj4+pfPZbDaKi4trHa9qOkSEboOj2LfxjCa2chpFYhORucAIIFxEUrAqgbwjIg8CywEbMNsYs9uNYaomQmy2aktWPpGj8J7zLrmbNhMwcAAB8fX/hEhGRgbR0dZt4jlz5lwybuXKlVy4cAF/f3/mz5/P7Nmz2bJlC3/+859Zt24dXl6N+Xa4cpeuA6PYuOAIqSnZhLcLcnc4jUajSGzGmCmVDF8KVN7JlFI/QEB8fIMktBJPPPEEd999Ny+88MIl3cCA1dP0zTffTEpKCnfeeSeJiYncc889XLhwobTSSGJiotZsVJcIbuVHu24t2f/NGcLbdXZ3OI2Gdluj3dZ4hKbc1cmcOXO0ckcZTXlfusO+jaf5+tPD3P2HK/GyNa6SvXZbo5RSqtY6xUVQmF9Myr6L7g6l0WgUlyKVas6mTp3K1KlT3R2GaqLsft5cEd+abZ8f5/yJLKK7tiSqU6i7w3IrTWxKKdXERXQMZv+HZzi57yI2by8mPxrfrJObXopUSqkmrijfajPSGHA4nJw80LwvS2piU0qpJq5d95aUPApps3kR3bWlewNyM01sSinVxEV1CiV+bEd8A7y5/pG4Zn0ZEvQem1J1Ii0tjdGjRwNw5swZbDYbERERgNUOpN1ud2d4qhnoc3U0W5cdwz9Iv2ua2JSqA+7qtkapEkEt/WjVNpDjey7QIjLA3eG4lV6KVM3WmSMZbFl2lDNHMupl+Q3RbY1SZXXo2YoTe6ru2aI50FNI5XGcTkNeZtW9Vp87lsnyt3fjKHZi8/Zi7LRetO4YUuU8/iF2vLykymnKa8hua5Tq0DOMXWtP4ihyYvNpvuUWTWzK4+RlFjLn11/VeHpHkZOlb+ysdrqpM68isEXt+r1qyG5rlGrTJRQMnD6cTrvurdwdjttoYlMexz/EztSZV1U5TUmJzVnsxKsWJbbaashua5Ty9rHRtmsLju+5oIlNKU/i5SXVlqxiW0Rww6PxnDxwscGaINJua1RD6NAzjL1fn+bKm9wdifvor0U1W1GdQkkYF9Ngz/w88cQTPPXUU8THx3+vE9CSbmv69u3LzTffTGJiIq+99lpptzVxcXFMmzatQeJUTVuHXq1IS8kmJ6PA3aG4jZbYlKpjzz33XIXDhwwZwoEDB0rfv/DCC6X/t2vXjvnz518y/bvvvlsf4SkP1yIygKBWvpzYc4HuQ9q4Oxy30BKbUkp5EBGhQ88wju+54O5Q3EYTm1JuNnXqVK0coupUh16tOLHnAk5n8+xIWhOb8hjNuTd4T6H7sG60696Kgrxizh/PcncobqGJTXkEPz8/0tLS9MDYhBljSEtLw8/Pz92hNHm+/t5ExYY021ZItPKI8gjt2rUjJSWF8+fPuzsU9QP4+fnRrl07d4fhEdr3bMXxPRdIHB/r7lAanCY25RF8fHyIjW1+P2ClKtOhVxiblxylIK8YX//mdahvspciRaSTiLwjIh+XGRYoIv8SkbdE5A53xqeUUu4U0SEYX39vUvY1v9qRbklsIjJbRM6JyK5yw8eJyH4ROSQiv65qGcaYI8aY+8oNvgn42BhzP3B9HYetlFJNhpeX0L5Hy2ZZ7d9dJbY5wLiyA0TEBrwOXAf0BKaISE8R6SMii8u9Wley3HbACdf/jnqKXSmlmoQOvcI4sftCs6tU5ZYLr8aYtSISU27wQOCQMeYIgIi8D0w2xvwBmFjDRadgJbdkKknaIjIdmO56my8iu6tYXihQWWddlY2rzfBwILWK9denqratIZZV03mqm66u9lFj2z9Qd/vInfunumma8j5qMr+hu//QIL+hioZ3q0F8dc8Y45YXEAPsKvP+FuDtMu/vAl6rYv4w4E3gMPCUa1gg8C7wBnBHDWKYdbnjKxtXm+FAkhs//yq3vb6XVdN5GmofNbb9U5f7yJ37x5P3kf6Gqh/urv3TZKvKGGPSgBnlhuUA99RiMYt+wPjKxtV2uLvUZTyXs6yaztNQ+6ix7R+ou5jcuX+qm6Yp7yP9DdV+XQ1CXFm14VdsXYpcbIzp7Xo/BHjOGDPW9f4pAGNdivRIIpJkjEl0dxyqYrp/Gj/dR42bu/ZPY6ruvxnoIiKxImIHbgMWujmm+jbL3QGoKun+afx0HzVubtk/bimxichcYATWjd+zwLPGmHdEZDzwCmADZhtjXmzw4JRSSjVpbrsUqZRSStWHxnQpUimllPrBNLEppZTyKJrYlFJKeRRNbI1QRQ08K/fTRrYbN/3dNH4icoPr9/OBiIypr/VoYqtj9djAs6oHtdxf2sh2A6vN/tHfjXvUch/Nd/1+ZgA/rq+YNLHVvTnUTwPPqn7MoYb7C21k2x3mUPP9o9xjDrXfR791ja8XTbZJrcbK1F8Dz6oe1GZ/UYNGtlXdquX+2dPA4Slqt49EZC8wE/jMGLO1vmLSH2fDiOa7M32wDpDRlU0sImEi8iYQX9K0mGpQle2vT4GbReQNGkmbeM1UhftHfzeNSmW/oV8A1wC3iMiMimasC1pia4QqauBZud9lNLKtGpD+bho/Y8yrwKv1vR4tsTWMk0D7Mu/buYapxkn3V+Om+6fxc+s+0sTWMJpjA89Nme6vxk33T+Pn1n2kia2OuRp4/hroJiIpInKfMaYYeBBYDuwFPjTGVNVzt2ogur8aN90/jV9j3EfaCLJSSimPoiU2pZRSHkUTm1JKKY+iiU0ppZRH0cSmlFLKo2hiU0op5VE0sSmllPIomtiUUkp5FE1sSimlPIomNqWaABHpIiJrRGSXiPxFRA67OyalGitNbEo1cq5OG98DHjPG9AbsgDYhpVQltNsapRq/G4A9ZTpm3Aukuy0apRo5LbEp1fjFY/XaXaIfsMM9oSjV+GliU6rxSwO6A4jIIOAnwHa3RqRUI6at+yvVyIlIOLAECACWAncAHYwxTrcGplQjpYlNqSZERNoDHxtjBrk7FqUaK70UqVTTovfXlKqGltiUUkp5FC2xKaWU8iia2JRSSnkUTWxKKaU8iiY2pZRSHkUTm1JKKY+iiU0ppZRH0cSmlFLKo/x/DCR2pptjG/AAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "bins = 253\n", + "\n", + "fig, ax = plt.subplots(1, 1, figsize=(6, 3))\n", + "\n", + "ax.set(xlabel=r'$q$', ylabel=r'$q^2 f(q) \\ dq$', xscale='linear', yscale='linear')\n", + "\n", + "q_auto, w_auto, auto_bins = get_binning(bins, strategy=0)\n", + "ax.plot(q_auto, q_auto**2*w_auto, '.-', lw=1.2, label=f'Auto, {auto_bins} bins')\n", + "\n", + "bins = auto_bins\n", + "ax.set_title(f'collective distribution function, m={masses}, {bins} bins')\n", + "\n", + "q_auto_tol, w_auto_tol, auto_bins_tol = get_binning(bins, strategy=0, auto_tol=1e-2)\n", + "ax.plot(q_auto_tol, q_auto_tol**2*w_auto_tol, '.-', lw=1.2, ms=6, label=f'Auto, {auto_bins_tol} bins')\n", + "\n", + "q_lag, w_lag = get_binning(bins, strategy=1)\n", + "ax.plot(q_lag, q_lag**2*w_lag, '.-', lw=1.2, label='Gauss-Laguerre')\n", + "\n", + "q_trapz_inf, w_trapz_inf = get_binning(bins, strategy=2)\n", + "ax.plot(q_trapz_inf, q_trapz_inf**2*w_trapz_inf, '.-', lw=1.2, label='Trapz inf')\n", + "\n", + "q_trapz, w_trapz = get_binning(bins, strategy=3)\n", + "ax.plot(q_trapz, q_trapz**2*w_trapz, '.-', lw=1.2, label='Trapz')\n", + "\n", + "#ax.set(xlim=[0.05, 5000], xscale='log')\n", + "#ax.set(ylim=[4e-9, 1e-2], yscale='log')\n", + "\n", + "#ax.set(xlim=[0.05, 1e+6], xscale='log')\n", + "#ax.set(ylim=[1e-12, 1e-2], yscale='log')\n", + "\n", + "ax.set(xlim=[0.05, 1e+2], xscale='log')\n", + "ax.set(ylim=[1e-10, 1e-1], yscale='log')\n", + "\n", + "ax.legend()\n", + "fig.tight_layout()\n", + "fig.savefig('binning.pdf')" + ] + }, + { + "cell_type": "markdown", + "id": "04263ac6", + "metadata": {}, + "source": [ + "# Perturbations" + ] + }, + { + "cell_type": "markdown", + "id": "6ce9d9e6", + "metadata": {}, + "source": [ + "Massive speed-up in perturbations! About a factor 4 when including perturbations!" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "fa4a5a30", + "metadata": {}, + "outputs": [], + "source": [ + "def get_delta_nu(bins, k_output_value=0.2, strategy=1):\n", + " collective = Class()\n", + " extra_settings = {\n", + " 'collective_ncdm': 1, \n", + " 'Number of momentum bins': f'{bins}, {bins}, {bins}', \n", + " 'Quadrature strategy': f'{strategy}, {strategy}, {strategy}',\n", + " 'output': 'tCl, mPk',\n", + " 'k_output_values': k_output_value\n", + " }\n", + " collective.set(standard_settings | extra_settings)\n", + " collective.compute()\n", + " pt = collective.get_perturbations()['scalar'][0]\n", + " return pt['tau [Mpc]'], pt['delta_ncdm[0]']\n", + "\n", + "klist = np.logspace(-3, 1, 10)\n", + "delta_nu_list = []\n", + "tau_list = []\n", + "for k in klist:\n", + " tau, delta_nu = get_delta_nu(10, k_output_value=k)\n", + " delta_nu_list.append(delta_nu)\n", + " tau_list.append(tau)" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "id": "f2b6d36f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAADkCAYAAABUp2nWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAABLyElEQVR4nO29e5xkV3Xf+131fp5TPaOZkTQzPZIsITNIwMwoggTHwQ9iQQjYDrElfGMTy+YqMSaBhAsYfwy+Nw5w8fUDG1uWDRF2eJiPLgEZy+CbgCM7JjHStAQSQiAE0zMjaaZnpuucer/Oun+cc3pqevpRVV3d9drfz6c+3bXrnFO7dleftddea/+WqCoGg8FgMAyTyKg7YDAYDIbpwxgXg8FgMAwdY1wMBoPBMHSMcTEYDAbD0DHGxWAwGAxDxxgXg8FgMAwdY1wMhjUQkV8SkT8adT92ChFREbl+SNeaF5GyiESHcT3DZGKMi2GmEJGXi8ipzY5T1f+oqj+3E33aCiLyBhH5mxH34bsi8sPhc1VdVNWcqnZG2S/DaDHGxWBYhYjERt2HXhhGPyflsxomD2NcDGNJMBv+9yLyVRFxRORPRSTV9fqrReQRESmKyN+KyAu7XrtkiUdE7hWR/yAiWeAvgKuDZZuyiFwtIu8RkftE5D+LiAu8IWj7z8H51wTX/BkRWRSRcyLyrq7rJ0Xkt0TkmeDxWyKSXOdzvUFE/oeI/G7wub4hIj/U9botIh8WkWdF5HTQ7+iqc39TRM4DfwrcDfz94LMUg+P+SkR+btV7/k3XcxWRXxCRbwHf6ureq0Tk6eDzfUBEIsHx3yMiXxSR88FrHxORQvDanwDzwJ8Fffg/usYrFhxztYjcLyIXROQpEfn5rr68R0Q+JSJ/LCIlEXlcRG7pev3twTiUROTJ7rEyjDfGuBjGmZ8AbgOuBV4IvAFARI4AHwH+d2A38AfA/evd0ENUtQK8EngmWLbJqeozwcuvBe4DCsDH1rnE9wE3Aj8E/IqIPD9ofxfwUuDFwIuAW4Ff3qArLwG+DVwBvBv4tIjsCl67F2gD1wNHgH8M/Nyqc58G9gH/G3AX8OXgsxQ2+vyr+NHgWoe72n4MuAU4ij8ePxu0C/Be4Grg+cBB4D0AqvovgEXgnwZ9+L/XeK9PAqeC818H/EcR+cGu118THFMA7gd+F0BEbgTeBPw9Vc0DPwJ8t4/PaBghxrgYxpkPquozqnoB+DP8mzfAG4E/UNX/paodVf0o0MC/wQ/Kl1X1M6rqqWptnWN+VVVrqvoo8Ci+IQH4KeD/VNWzqroE/CrwLzZ4r7PAb6lqS1X/FHgS+Ccisg94FfBvVbWiqmeB3wRu7zr3GVX9HVVtb9DPXnivql5YdY33B22LwG8BdwCo6lOq+v+paiP4fL8B/KNe3kREDgIvA96uqnVVfQT4I+Cnuw77G1V9IIjR/AkXx7UDJIHDIhJX1e+q6rcH/sSGHcWstxrGmee6fq/iz3wBDgE/IyK/2PV6ouv1QTg5QH9ywe9XAye6XjuxSV9O66WKseHxh4A48KyIhK9FVvWtl372wlrX6W5b+QyB0ftt4B8C+aBPyz2+z9XABVUtrbr2LV3PV49rSkRiqvqUiPxbfC/pBSLyBeCtXd6mYYwxnothEjkJ/JqqFroeGVX9RPB6Fch0HX9l1+/ryYBvRR78GXzDEDIftK3HfumyHl3Hn8T3wK7o+lyWqr5gg36u1e8K63/+jc47uEafAP5jcPzNqmrhL8d193+jsXsG2CUi+VXXPr3BORcvrPpxVf0+/PFV4P29nGcYPca4GCaRPwTuEpGXiE9WRP5J1w3sEeD1IhIVkdu4dAnnDLBbROwh9ucTwC+LyB4RuQL4FeA/b3D8XuDNIhIXkX+OH8d4QFWfBf4S+H9ExBKRSBBM32gJ6gxwQEQSXW2PAD8uIpkgseHOHj/H20RkLljK+jf4CQPgeytlwBGR/cDb1ujDdWtdUFVPAn8LvFdEUuInXtzJxuMD+DEXEfnBIJZWB2qA1+NnMYwYY1wME4eqPgT8PH7gdxl4iiDYH/BvgH8KFPHjIZ/pOvcb+MbgafEzzbaylBbyH4CHgK8CXwOOB23r8b+AG4BzwK8Br1PV88FrP42/xPd1/M92H3DVBtf6IvA48JyInAvafhNo4t/0P8r6CQqr+SzwML5x+nPgw0H7r+IH+Z2g/dOrznsvvnEtisi/X+O6dwDX4Hsx/wV4t6r+1x76kwTehz9Oz+Eb5Xf2+FkMI0ZMsTCDYecQkTcAPxcs9RgMU4vxXAwGg8EwdIxxMRgMBsPQMctiBoPBYBg6xnMxGAwGw9AxxsVgMBgMQ8fs0A+44oor9Jprrhl1NwwGg2FiePjhh8+p6p61XjPGJeCaa67hoYceGnU3DAaDYWIQkRPrvWaWxQwGg8EwdIxxMRgMBsPQMcbFYDAYDEPHGBeDwWAwDB0T0DcYDLNFqwaqEEtCJDrq3kwtxrgYpo4HTz3Ig6ce5JdfulGl4dnBU1+lPiKzt1DhqcdyfZnnqs9xpnKG5559mPP/83fpoHgIKhE8ieBJFC8S9X8Gz1UidIKfK8dIZOV5RyIo4WuCR/ia/7v/U/xzCH8XNGhHIsRjcRKxOPFYjGQ8QSIWIyoRRISIRIhIBMH/XURWfu9uX/2aiBDhYntEIkSIgLBmu520+bEbfmzoY2+Mi2Hq+NKJL/Jfvv0Z3vWSd3FpTa7Z5Mf+6BjfTXS4Ir2HKzJXcEX6Cvak91z8mfF/hm3xaHzUXe6JywxH5TnOVC/+PFM5w5nqGVpe65LzIrZfTFMC5SvpegBEFIRO8GgSCV6LqCIKghIBfJMBUdWL5wZtAsGxwblBPTVZub5PA6iJXwVNAS9414uGCd8wifjHiATP6XruF7nxgvfwgoeil7Svx3x6rzEuBkMvuM8+TEc7VFoVconc5idMMc8VqyzGW1zbSnPz9f+ApdoSZypneOzcYyzXl9E1bjx20l7X+OxO7/Z/z+whG89uW7899bhQv3DRWASGotuAnK2evcxwxCIx9mX2sS+zj5v33Mwrsq9gX2YfV2av5MrMlZz/7/+F7/vaB/jCj/x35vYdJCKCSGAEBKD7uay0S9i+1u8oEa+JeE2k0yLiNYh6LSKdJuI1iLSb/uudRnBMk0inideqU6/XqNVq1Bs1GvUGzUaNVrNBu1mn3Wrgtep47RZ0GsS1TZw2CVrEpU0C/5GUi4+EdIjTIqZt4tpcGZeLBigYX/ENmQLMbc/SoDEuU4Cq8l8X/ysvP/hy4pHJmHVuJ0696P+sXZh54/LQt79DW4S/V4/xrpf9X5e81vbaXKhfYKm2xLnqOZZqSyu/n6v5j4fchzhXO3fZTRwgHUtfNEKZVcanyygVkoVLPMi1DEfoffRqOF6050Xsy+67xHDsy+5jV2rXhst/x2ufIALc+vzr2F2wtja4O4iq4tbbLFeaXKg2/Z+VJsvVJhcqrUvbg5/FepOodkhw0SglaJOQNrmYxxUp2JWCq/N5/t029NkYlyngGxe+wVv/6q389g/8Nj84/4Oj7s7IOd+qQQSc0mn22/Oj7s5IeezEtwAotJqXvRaLxNib2cvezF7Yvf41VBWn4fiGJzA6S7UllqpLK78/eeFJ/qb2N1RalTXf54r0Fcwl53Cbbs+GIzQYV2avZF9mc8PRE7VlqprEzk/WpENEsNNx7HSca+jNY+x4ilNrdRmh5irj1GK52qST2h4zYIzLFHCh+F3/59ITYIwLrvo3rpNLJzh84O+PuDej5elnvwM52NWuDXwNEaGQKlBIFbhh7oYNj622qheNT22J87XzLFX93y/UL3Ctfe32GY5ePku9iEuWK6PTn9wQjQi7sgl2ZRMjeX9jXKYA9/yTwc9vjrgn40FZ/JXlbz5zgh85MuLOjJBWx+Ns8TTkYE9ncOPSD5l4hvn4PPPWeHqMsaZDOZIfdTdmguk33zOAU13yfzaKo+3IGND22lSDb/XihWdH25kR841nS8TEAeCqTgX1vE3OmH4SLYda1BiXnWDijYuIXCciHxaR+7raDovIp0Tk90XkdaPs307g1Jf9n013xD0ZPaVmaeX3pdLSCHsyeo4vLhOP+uOxy2tTrZjvR6rtUo9NTiB/khmpcRGRj4jIWRF5bFX7bSLypIg8JSLv2Ogaqvq0qt65qvmVwO+o6r8CfnrI3R47Qo/FbZVH25ExwK2cXfm91CxSb3VG2JvRsrC4TCbpL4fZnkdpebaNLUC6U6aVsEfdjZlg1J7LvcBt3Q0iEgU+hG8gDgN3BJ7IzSLyuVWPvetc90+A20XkA2yYBzMduMFs3W1XR9yT0eOUTl98Eqnx+DPO6DozYhZOFsmmmsRUSatSKZ4bdZdGTl5LtI1x2RFGGtBX1QdF5JpVzbcCT6nq0wAi8kngtar6XuDVPV73LPALgaH69BC7PJY4gcfidBoj7snoueA+s/K7F63z8Illjh3aNcIejYZz5QYnzlc5uruB5XkIUCvNtnHRVp00DTQ1N+quzASj9lzWYj9wsuv5qaBtTURkt4jcDRwRkXcGbdeIyD3AHwMf2ODcN4rIQyLy0NLS5C4ZuJ06AI5evtFt1lhynwNgV6eDxlocP1EcbYdGxCOLRQDakSZ2xw/kN0vnR9ij0VMvXQBAMoXRdmRGmPhUZFU9D9y1qu27wBt7OPce4B6AW265ZWMBnjHG9XyPxcVkA50v+zGXA60OTrTDw4vLqOrMaYwtnFwmFhFq0iAv/r95u7I84l6NlnLxLGkgmpk9T3YUjKPncho42PX8QNBmWAdH2wBUhDVlOmaJYtWfnV5FnFqkw1KpwanlndnjMU4cP1Hk8NUWrtfCjqUB8CoXRtyr0VJ1fM8tkTPLYjvBOBqXrwA3iMi1IpIAbgfuH3GfxhoHJR3sYXDrsxvABj9zLut52HGbcsQfk+OLszVj73jKo6eKHDlYwNUOdjxHU6MQaK7NKjU3MC75qc/xGQtGnYr8CeDLwI0ickpE7lTVNvAm4AvAE8CnVPXxUfZznKm36zQFDrZ878WtPDfiHo2WUquE1fGwkwWqApkEHD8xW8blm2dKVJsdjszP4Ypix3O4kkdm3Li0yr7nlraNcdkJRp0tdsc67Q8AD+xwdyYSp+7/wxwkzjcJUnH33DTaTo2QcqdK1oNCyoYyvGB/jIdnzHNZCIL5L9qfpfRVwU7kqURyxJvFkfZr1LSDZcGsvWfEPZkNxnFZzNAH4b6Og0HuvluebcmTstcg04mwKwja3nRVhyeeLVFttkfcs53j+OIyu7MJCjH/ZmolbWrRPInWbO/Q96oX8FSw564YdVdmAmNcJhy35BuTg9krgYs6Y7NKlRZJjWGn/aWP6+Yqfgzi5OzEohYWlzkyX8ANJh5WaheNmEWqXdrkzOnGV0TOkE2amkc7gTEuE44TxFgO2tcC4FZne6NcWdqkvDh2xl/62JP2N5jOSlDfqbb49lLFj7eU/e+GndlDM2GT9Wbbc4nUi5QkP3Np6aPCGJcJx634nsqB3YeBizGYWURVKUeUJCns7D4AWu3zXLcnOzNB/YWT/uc8Ml/AqfrGxcrupZMokPNmW3su3nSpRCarSNgkY4zLhOPUfE9l197D5DseTmN2ln9WU+/UaYmQlCxW7moA3No5js3PcTzYTDntLCwWiQi88EABJ/Bi7exVaLpAXmq016hIOSsk2w41o4i8YxjjMuG4jSIxVTKF67DUw23N7rq6G6TaJiM58vn9QdsyRw/NsVxt8Z1zl5fgnTYWThZ53r48uWQMt+bv67Dy+5G0v3GwVJxdCZhUu0TTGJcdwxiXCcdpur4wYWYXlkZwWrOrjBxmzqXjNrHsbnKe78kdO+TfWI8HKbrTiucpC4u+MQVWvFjLPkAs67eVi2fXPX/ayXklI7e/gxjjMuE4zRKWB0Rj2BLDCUQsZxGn5CsiZ5K7IJHD9jzcZonr9+TIp2I8POVxl6fPlSnV2xw5WADAbbpkPSUezxDP+dlzVWdGEz5UyVHGSxVG3ZOZwRiXCcdt17CJAmBHkrje7K6pLzm+ccml9oCI78m1K0QiwpH5ORamPGMs9MyOzAeeS6uMFYSZkoHkSagMPGt06i4xPNQYlx3DGJcJx/HqWBFfaMGKpXCZ3cqL59wzAOSDTDFLoiue3NH5Ak+eKeHWp1fYc2FxGTsd57orskAw8RB/4pEJJE9aMyq7Xyn6WZWRjBGt3CmMcZlwHK+FHU0CYMdzOHgzkRW1FstBdlTBOgBc6skdOzSHKjx6sjiq7m07C4tFXnywQCTi7+NwvQaW+BsGswV/30+nOt3e23qUgyqcsayR298ppsK4iEhERH5NRH5HRH4maLtORD4sIveNun/biUsHO5oBwI7n6YhQndFyx8X6MlFVdhX8TDErmlopR/DigwVEmNriYeVGmyfPlDgyX1hp8yceKQCsgi954s2ocam5vnFJGkXkHWPkxkVEPiIiZ0XksVXtt4nIkyLylIi8Y5PLvBa/7ksLv3Ilqvq0qt65Pb0eDzpeh5KAFfeXQayknwnjVGYzaFtqutieRy7QjrLjGVwUVSWfivO8vfmpFbF89GQRVTg6f3HZx8XDivkTj3giSVnTSH06P/9mhLGmZN54LjvFyI0LcC9wW3eDiESBDwGvBA4Dd4jIYRG5WUQ+t+qxF7gR+FtVfSvwr3a4/yOj1PDlPOzAqNgp/x/HKZ8aWZ9GSaldIddR7FweACuepy1Qa/vFwo4e8oP6njd9y4ZhssKLgkwxVcWJgBW/uCO9JDmiM7rJth3I7WeMIvKOMXLjoqoPAqtTWG4Fngq8jybwSeC1qvo1VX31qsdZfG8lnJLNTETbrfoBbCvIgLEy/ow9FLOcNcqdGllPKGT8OENodN1gtn50vkCp3uappemTQVlYLHL93hx22v/s9XqRpgh28uKmwWo0T7w5m8alE8jt5wpGEXmnGLlxWYf9wMmu56eCtvX4NPAjIvI7wIMAIrJbRO4GjojIO9c6SUTeKCIPichDS0uTpyYcbhoMPRY7u9dvr5wZWZ9GSUUbpLwoqbifIWWFnlwwTiubKadsv4uqsnCyyNGueIvr+v8+VvLiMlktapFsz6Z4pdaWaWgcO2926O8UIy0WNixUtQrcuartPHDXJufdA9wDcMstt0zcWokTeChWxjcqVu4qv706m7uwK9Jmr3fxKx3K7oebK6+9IstcJs7DJ5a5/db5kfRxOzhxvsqFSnNlfwtc/MzhGAA04xa7a9/Z8f6NA5F6EYcce4OJh2H7GVfP5TRwsOv5gaDN0EVoROygloud91NwndpsbpQri0dSkyvPQ2XksPSziHA0ELGcJrqVkENCuX0rczHG0E7a5LzZ1J6LNh3KRhF5RxlX4/IV4AYRuVZEEsDtwP0j7tPYEdZusXK+cUnl9pHwFLdRHGGvRoOnHhVRkpJeaQvHxemqcXP00BzfXqpQrE6PksHCYpFcMsYNe/MrbRcnHvtW2jrJAnkto563430cNYmmQzWa3/xAw9AYuXERkU8AXwZuFJFTInKnqraBNwFfAJ4APqWqj4+yn+NIWLvFCjYNSqqA5Xm4zdlbVy81S2ggtx9ir8juX9yVHqbqLkyRiOXxxWVedNAmGrlYBMsNSjHY+a5QZXqOhLSpVWfPe0m2S9SNcdlRRm5cVPUOVb1KVeOqekBVPxy0P6Cqz1PV71HVXxt1P8cRp1Ek63nEs8HSRySCrb6m1KwRFk1LdUmqp3NXEVPF6drbEd6Ep0XEstbs8MSzJY4cvFTWJMyQCycecFH6pFScvX1QmY5LK1EYdTdmipEbF8PguK0ytqcQz6y02RLFnUFl5DAjLBMvrLRJZu4yTy6TiPH8q/JTE3f52mmHjqccPVS4pN1pOERVyXYti8UD2f3KDBqXnJZpJ0ym2E5ijMsE47YqWCrQVRPckgROpzHCXo2GC25gXFJd8h7ROJZ3uSd3dH6OR04WaXcmP/YQGskXr/ZcmiUsBYlc/BdPhMrI7owZl06LLDU0ZUQrdxJjXCYYp1PDlkuzya1oEkenV/l3Pc46flp2Pr3vknabCM4qrbVjh+aoNjs8eWbyYw8Li8tcszvDrmziknanXcVWuaQtZfkbCGdNdr8eKkGnCyPtx6xhjMsE43pNrMilNxU7lsFl8mfk/XK+7Mdc7PxVl7TbkTjuKk8uDOpP+mZKVeX4YvGS/S0hTqeGtWrikQ1k99vl2ZLdDxWRo1njuewkxrhMMI62saLpS9rsRJ6qQKszW95LserfMHfZBy5ptyJJ3FWe3IG5NHvyyYkve3y6WGOp1LhkZ36IP/FIXtKWC2T3vdpkG9V+CatvxrNG+mUnMcZlQlFVHBS7K5gPYAVBS6deHEGvRofbKJL2PKy5VZ5LLH1ZATV/M2Vh4jPGFlZVnuzG0Q52LHVJWy5foK0RdMZk92uuP/FIGEXkHcUYlwml1q7REt9T6cYOgpZuebYEDdxWCcvzsOxLb7RWPEtJ/PIE3Rw7NMfihSpLpclNflhYLJKKR7jxysv3b7iiWLHsJW0SieBKjsiMKSM3S77nkrGN57KTGOMyobiBxIuVLFzSHmpJuYG21KxQ7lTJdaCQuXS2bid8ZeSwPEHIStxlglOSjy8u88IDBeLRS/+NO+0WpTUmHgBlyRObMQWHdsX/G4cxJ8POYIzLhLJaETkkFLF0Aj2tWaHi1Ul7EfKpVdlzQTkCp3qpUvRN+23iUZlY49Jod/j6M+4lemIh5fIzqAh28vLXatE8idZseS5h9c38nPFcdhJjXCYUt3y56i1cFLF0KrOljFyhRdqLrtSPDwmN72pPLhWP8oKr7YnNGHv8GZdmx7tsZz6AW/KLxVmpy2MM9bhFqj35Kdh9US9S0jT5dGrzYw1DY+KMi4hcJyIfFpH7utqyIvJREflDEfkpEXm+iNwtIveJyFRWpnSDmi1WdtW+jvzlelqzQEU6JDVxWXuoCuyULy+gduzQHF895dBsT17qdmgU18oUC0sx2JnLl4FacYvsjCkjR+tFXMkhIpsfbBgaO2pcROQjInJWRB5b1X6biDwpIk+JyDs2ukZQnfLOVc0/Dtynqj8PvEZVn1DVu4CfAF421A8xJoSeiZ27NDsqZ/nGJRS1nBXK4pHk8plpaHzXKqB2dH6ORtvj689OntDnwski+wtp9lqXf+aViUdm32WvdZIFcjpb2nOxZpFKxEi/7DQ77bncC9zW3SAiUeBDwCuBw8AdInJYRG4Wkc+teuxd57oHuFi5shNc9zXAnwMPbMPnGDmh8Qg9lZBoejf5joczQ8rIjU6DRkRISuay10LjG6oEdxPqcU3i0tgji8U14y0ATi3YUBqUHOjGS81hUaHTbm9n98aKRMulZhSRd5wdNS6q+iCwekp9K/BU4JE0gU8Cr1XVr6nqq1c91gsknMI3MBB8JlW9X1VfCfzUNnyUkePWl4mpkl7luRBPYaviNmdn6cOtFQFIrlEMKlQFdtbYOHiVneZqO8XDExbUP+PWOV2srWS8rWYlkzB/eWVwCSRQZkkZOd1xacSN57LTjEPMZT8XvQ7wDcXl/xUBIrJbRO4GjojIO4PmTwP/TER+H/gzEXm5iHxQRP6ADTwXEXmjiDwkIg8tLS1t/ZPsIE7TxfK8lZtFN5Zerqc1zThBckM6Zl/2WjxzBRnPW7eA2tFDcyxMmOeysHh55clunMblcvsh0Ywf5C/PkHHJemVaicu/G4btJbb5IeOFqp4H7lrVVgH+5apD/6qHa90D3ANwyy236JC6uCM4rTK2ApHLa4LbkdhMye4Xg0ywTHKNHdjJPJbn4bTW9uSOzs/xua8+y7NOjavs9JrHjBsLi0US0QiHr157Nu40XNKekkhevhQUD3apV2dFGVmVvJboJI1x2WnGwXM5DRzsen4gaDNsgNuuYXO5YQGwI0lcb3a0xZYc/+uSS+25/EURbI3gtiprnnvsUChiWdyu7g2d44vL3LTfIhlb++/vl2JY+9zUiuz+bGQTeo0KcToQ7Hcy7BzjYFy+AtwgIteKSAK4Hbh/xH0ae5xOHUvia75mR9M4zE7A9pzrZ0fls5dnR8HGBdSef5VFMhaZGJ2xVsfjq6ecNfXEQpzO+hOPdCCB0pwRZeRyIFopaaOIvNPsdCryJ4AvAzeKyCkRuVNV28CbgC8ATwCfUtXHd7Jfk4irLexocs3XrHgGB8XTydu/MQgXArn9XWvEGACsSALHW1tDLBGL8MID9sTs1P/GsyUabW/deAuA22lgRdaeeGQD49KpzEaqemXZ/27Eska0cqfZ0ZiLqt6xTvsDTGnK8Hbh4GHHLk+9BbATFl5DqLQq5NfQl5o2nHoRUWVubm3jYkdTuN76qdlHD83xkb/5DvVWh1R87Rn/uLBwMtw8uYHnoi0OxS7PnAOw5gLZ/RlRRr6oiGx0xXaacVgWM/RJ22tTFrDi69xAgvVld0YkYNym4ysiz60RcwGsmO/JrcfR+TlaHeWx0+OvuXX8xDL7rCRX2etLmbjaWXfikUimqGgKmZGSDPVAETlpGeOy0xjjMoGECr9Wcu1soVBTKhS3nHbK7TJ5T7Fz6xjbRJ6mQL29dtxlkhSSF04WOXJwbkMpE3eDiQdASXJEZ0R2v1n2l/+ylhGt3GmMcZlAwn0ddnLtpRE77f8jOTMiu1/xamQ7gpVeO85gBXsc1pPE2ZNPMr8rM/ZB/fPlBifOVzeMtzTqDvWIYCfW3zRYieSJNWfDuHhV/28eVuE07BzGuEwgocKvlV7b1beDrCm3OhvLYmVtkvai68ZL7CBTaKMaN8cOzXF8sYjq+G53CitPHj20frzFdQNF5DXk9kPqsTypGZHd12qRlkbJW2afy05jjMsE4pT9Wi12Zp0YQyAJ41RnY6NclRZJb/3clF48uaPzBZZKDU4t14bev2GxcHKZWES46er1b5ThUuh6Ew+ARtwm3ZkNeSCpF3HIkUpM3H7xiccYlwkkLHwV1m5ZjR2k5Lr1GdnLIB4pLpfbD1mR3d+ggFroDYxz3GVhscjzr7JIJ9bPaFuZeGTXXwZqJ+yZkd2PNoqUZf34k2H7MMZlAglrtVirFJFDUrl9JDfQ05omPPUoR5Qk60u32NlAGbmyvn7cjfvyZBLRsY27dDzl0ZPFNeu3dBMuha6u89ONl7TJz4jsfrzlUDWKyCPBGJcJxFlRvV17XwdJC9tTnMb0z04rrQqeyJqKyCGhEd7Ik4tFI7z4YGFsPZdvnilRaXY23JkP4ASlBVbX+elG03OkpEW9Ov0GJtVyqcWMIvIoMMZlAnEbDjnPI5ZdZ11dBAvBaa+tpzVNhEXT0hvMTnP5q4mq4myyt+Po/BxPPFui2hw/6ZwwmL9RphiAW99k4gFEMr6BKi1PlhL4IKS9Ei0jtz8SjHGZQNxWCdtTiK+/FGQRxW2Pb3B6WIQZYOlEYd1jJD2H5Xm4m6TfHjs0Fyw/jV8m1fHFZXZnE8zvWntzZIjTcIioktvAcwmlUELdrWkm55WM3P6IMMZlAnFaFaxN/nR2JL6untY0cS5QRM4mN9iBHUtgeeA2N14GCr2CcVwaW1hc5sh8YdM68E6zhKVKJLp+dlQohVKbduPidchTRY0i8kiYOOMiIteJyIdF5L5V7dmg8Nerg2Jhfy0id4vIy0fT0+3DV0TeOLXSiqZwdPyWd4bNOfdZAHLp9Spg+9hsXkCtkEnwPXuyY1f22Km2+PZSZdN4C4Ry+xv/W6cDKZRGabqzCRvB7nyjiDwaBjIuIvLrIvL7IvLCPs/7iIicFZHHVrXfJiJPishTIvKOja4RlEO+c42X3g58KjwMKAMp/MqWU4XjNbFl/dRbADuWxWX6VZHPlfyYS2GdzLmQfCSO24Mnd3R+juOLy2O1mTIUq9ws3gLgdurYm0w8MrafptyacmXksNqmZIxxGQWDei7vAd4C/JSIvKKP8+4FbutuEJEo8CHglcBh4A4ROSwiN4vI51Y91pyeBn34OhBuSf9rVX0lvsH51T76NxG4dLBiG1dNtOM5agLNTnOHejUaikHm3C573crYANiRBI63+VgcOzTHcrXFd86NTzLEwmKRiMALDxQ2Pdb1mliRjSceuYK/qVSn3LhUin7CQjxn5PZHwaDbVl+NXzEyDfwRcKiXk1T1QRG5ZlXzrcBTqvo0gIh8Enitqr43eJ9eeDmQxTdONS7K9y8Daxc9mVBU1Zfbj2c3PM5K2lAFt77MFRvseZh03PoySc+jMLexcbFi6U2XxeDiZsqHTyxz3Z7x2Hy3cLLI8/blySU3/3d1tM2B6MbZUbn8HG2N4NWKQ+rheBJW20zmjCLyKBjUc2kADwLvB67bYh/2Aye7np8K2tZERHaLyN3AERF5J4CqvktV/y3wceAPgR8VkT8A/gT43Q2u9cYgTvPQ0tJkpGXW2jXaItib1GmxU/5NctqVkd1WyZfbL2w8O7XjOUo9FFC7fk+OfCrG8SD1d9R4nrKwuLyhnlg3juimE49INEJJskTq4xVbGjZhzCVtG+MyCjadConIzwIfU9WGiLwWuBL4mqr+3bb3bg1U9Txw1zqv3dv19NM9XOse4B6AW265ZXwW2TfArfpG0EpunF4Z6mm5pWfgyqPb3q9RUelUyXlgZzd2UO1EHq0LpWYJe4Oxi0SEI/NzYxPUf/pcmVK9zZGDhU2P9TptSuKXGNiMsuSJNdcvoDYNtINSzll742QPw/bQi+fybwLD8h7grcC1wLtF5G9FZG1xq/44DRzsen4gaDOsgeP6Q2OnNp6pr+hpldfX05oGKl6dTCdCbhNhwlAlODTOG3Fsfo5vni3h1lvD6OKWOL6yeXJzz6VceQ5PBLuHfR3VaJ74lCsje9UiAPmCqeUyCnoxLmEU9FXAD6rqO1T1R4BfA35vCH34CnCDiFwrIgngduD+IVx3Klmp5ZLZ+B/Gzvl235ly2f0KTVJelEhk4/0fdto3xqEk/UYcPVRAFR4Zg6WxhcUidjrOdVdsvNQFXXL7qc0NUSNmkW5Pt+dC7QJlTZHPbpz8YtgeejEuJ0XkXmAvXFQHVNU/x/diekZEPgF8GbhRRE6JyJ2q2gbeBHwBeAL4lKo+3s91Z4kVYcLMxkF6K++Hrdx1CmRNC2XpkNKNs6MArHToyT276bEvPlhAZDw2Uy4sLvPig4VNjSd0Tzw2L4zVnAHZ/UjD8WNLPYydYfj0ki32BuCfAb8J/L8i8nngceAIF72anlDVO9Zpf4CLGV6GDQg9kdAzWY98fj+iijPlQduKeCR6SAhcKaBW2dyTy6fi3LgvP3KF5HKjzZNnStx2U2+rz07ZL8Vg9WBc2kmbfGm6jUus6VCOGEXkUbGp56Kqrqr+J1V9FPjn+AbpDcA88JPb2z3DajZVRA6IZHaR9zycKa6V3uq0qEWElGystwXdBdR6ywo8emiORxaLeN7o8jy+erKIam/xFrgYT9pIEXmFVIG8VvE6na10caxJtBxqRm5/ZPSVihwYmg+o6utV9RdU9bvb1C/DOrj1IjFV0pt4LkTj2Oqn6k4rocpxagO5/ZB+lwmPzs9RarT51tnRydKHy3Iv7iFTDC6WFLA3USsAX3Y/IkrJmd5l01S7RCNmRCtHxcRpi806TtPF7nhIavN/GpsITmt6lZHdIMaQjm8+FsncPtJ9eHLHxqAy5cJikev35rDT8Z6OdwNja1kbe7UAsaz/+UKJlGkk67lGEXmEGOMyYTitMjZAZPM/nSVxXK++7X0aFcUgLTuT6EHeI5kn7yluj3s7rtmdYVc2MbK4i6qycLLY0/6WEKfpkvKUZA8qwLFg13rVmYzNw4OQ0wqdTfaDGbYPY1wmDLddw2b9Gurd9KqnNamcdXzPJZvuYR+DCLYKTqs3zTAR4cgIK1OeOF/lQqXZ88588CceVo8holQgu193p9Nz0WaVFE20h7Rsw/ZgjMuE4XqNTYUJQ6xYBlenN2B7zvU3iNqZ3rKpLInidnpfJjx6aI6nlyosV3beQPejhBzitqtY0tvEIxPI7jfL0xlzqQS1aiRdGG1HZhhjXCYMR1vYkd60OK14Flc219OaVJYr/g2kYG0ewIawgFrvhuJokKUV3uh3koXFIrlkjBv29p7t5HQam8rth2QLfrpye0qVkcNYUjRrFJFHhTEuE4aDhxXbPPUWwE5YeCKUm9OZMRbK7V8xd3CTI336LaD2ooM20YiMJO5yfHF55f17xdUWViTV07ErsvuBRMq0UQ08l4SR2x8ZxrhMEC2vRUXASvQmBW+Helo9bBycRNymS77jYc31JkxoxzKU+iiglknEeP5VeY6fKA7Yw8GoNTs88WyJIwf7ixc42sHepM5PSCqdoapJqE2n5xJW2UzmjSLyqDDGZYIo1f002l6ECeGinpZTmrpinACU22XynoeV37h+SYgdz/ddQO3Y/ByPnirS7uzc0uLXTjt0PO0r3gLgir8U2islyRGZ0k22zVKoiGxEK0eFMS4ThFMKtKN6zIBZUUYuba6nNYlUvBpZT3reBxKWKXD7kMQ5emiOarPDN57buaXFMEOt1535AM1GmVpEsJO9GVqAaiRPvDmdxqVT9ccwjC0Zdp5BK1GODBG5DngXYKvq64K2eeCDwAXgm8D/AH4K//MdVtV/MKLuDhW37O/rsNK9uforelpTqoxc0QbpToRErLc5UncBtV6rc64E9ReXuWn/zuyZWFhcXtln0ytu4J2GS6G9UI3lSU6pMrJXXaatEexNisgZto8d9VxE5CMiclZEHlvVfpuIPCkiT4nIOza6hqo+rap3rmq+GbhPVX8WOKKqf62qdwGfAz461A8xQsLaLHamxxhDoDHVSw2TSaRCi5T25rXARaPs9qCMHHJgLs2efHLHgvqqyvHFYl9eC4AbVBy1Ur3HGJoxi3R7OpM9pF6kRJbUJnV+DNvHTi+L3Qvc1t0gIlHgQ8ArgcPAHSJyWERuFpHPrXqsd1f9n8CdIvJF4PNd7a/HL308FTiV3hSRQ0JxS2dKg7YV8Uj2ILcfEhrlXmT3Q0SEY/NzO1b2+HSxxlKpwdE+4y3h0udmdX66aSVsst50ei7RRpGS9Jb4YtgedtS4qOqD+EtX3dwKPBV4JE3gk8BrVfVrqvrqVY/11nf+JfBuVf1B4J/AylKZo6rrTs1E5I0i8pCIPLS0NP6z+1CY0OpF9RZfTyvlebhTuK6uqpQiSoreC0FZuXCZsL+/9dFDBRYvVFkqNfo6bxAW+qg82U1YisHqcbkPoJMqkNfRCXNuJ/GmS9UoIo+UcQjo7wdOdj0/FbStiYjsFpG7gSMi8s6g+fPAm4P27wZtdwL/aaM3VtV7VPUWVb1lz57xD/yFir75TeT2V0hksTzFmcJ9LtV2lY4IyUhve37gojKyUzvf13vtpIjlwmKRVDzCjVf2d2N0q/6+jl4UkUMkNUdamtRrvUniTBLJtkMt2ntyg2H4TNyCpKqeB+5a1fYY8LpVbe/eyX7tBE7DId/xiGV7XFcXwUJwWtM3Ow1vpuk+biD53NUDFVB7wdU28ahw/MQyP/KC3pYkB+X44jIvPFAgHu1v3ufUe6vzcwmZQBl5+RypdO8pzJNAplPifLqPsTAMnXHwXE4D3VusDwRthlU4zTKWKsR6k38BsInidKZPGdkJ6sWn44WezwkLqLl97u1IxaPctN/eds+l0e7w9Wfcvve3ALgNB1El34fnEg9k90Mdrmki65WM3P6IGQfj8hXgBhG5VkQSwO3A/SPu01jitqtYff7J7EgC19v+WMFOc97x5x/ZPrKjiCWxlYE8uaPzczx6yqHZ3r7NlI8/49LseH3vzAdfbj+vSiTae/ZcPJTdnzZlZM8jrxW0h5pHhu1jp1ORPwF8GbhRRE6JyJ2q2gbeBHwBeAL4lKo+vpP9mhScTr1nYcIQO5ruS09rUjhT9LOj8une0rJDLCK47Wrf73fs0BzNtsfXn92+7KowmN9vphiA065gaX//zulg93rD7S8GNe60ag4RUUgbuf1RsqMxF1W9Y532B4AHdrIvk4jjNbmyR0XkECuWwW0Wt6dDI+RCOUjLzvcXA7ElhjvAMmG4mfLhE8s9lx3ul+OLy+wvpNlr9SY+2Y3brmP3KLcfkrF849KaMmXk0vJZdgGSMRsoR8k4LIsZesTVNlaPwoQhdiJPXaDRma6lseUgoL+7MN/XeVYkieu1+n6/K+0U+wvpbY27PLJYHCjeAuB6Tewe6/yE5AJplM6UGZdK0ffEwpiSYTQY4zIhqCquKHaPcvshdqinVZ2upQ+nUSSmSqHQ256fEDuaxmGwZcIj8wWOb9NO/TNundPFWt/7W0JcbWFF+/N48vYcngrUigO957hSC2JICaOIPFImLhV5Vqm2q7RFsBP95e5bXXpae/rIJBp3yq0SdsfDKvSnemvFs7idc6gqIr3XSgE/7vK5rz7LM8UaVxf68yA3YyHwiAaJtwA4otix/tKJI9EoRckifaZmjzvNcpCmbm383Wi1Wpw6dYp6ffqyKYdNKpXiwIEDxOO9J4wY4zIhhDVZrD6ECQGsoL58P3pak0C5UyHnKYVsfzEoO5Gn0xAqrQq5HuvihIRxl+OLy9tgXIokohEOX93/xj/1PF9uv8/PA1CWPNEpk91vlQNF5E3k9k+dOkU+n+eaa67pe6IxS6gq58+f59SpU1x77bU9n2eWxSaEcF9HWKOlV+xsqKf13ND7NEoqXp2MFyGX7G9+ZKUKADiVM32/5+GrLVLxyLYUD1tYLPKC/RbJWH9BeYBK5Qwdkb4UkUOq0TyJ1nTpi3lVP4aU28Srrdfr7N692xiWTRARdu/e3beHZ4zLhBCKLVp9CBMCWIHIpTNlsvsVmqS9aN83BivlG2c3qI3TD/FohBfuL/DwkIP6rY7HV08XVzyjfgmLwfXr1QLUYxap9nR5LlotUtMEVn5zCR1jWHpjkHEyxmVCcCq+59GPMCGAPaCe1rhToUOyD7n9EDtYJnTK/RsX8IuHff0Zh3qrM9D5a/GNZ0vUW97gmWKBoex34gHQjNukO9MlDxRpFHElRzQy3oajWCzye7/3e9v6Hu9973u5/vrrufHGG/nCF76w5jHf+c53eMlLXsL111/PT/7kT9Js9l6pdSOMcZkQVoQJe1REDsnlriaiitsobkOvRkc54pHS/veDhMY5jGH1y9H5Aq2O8rXTw5vtL5zsv/JkN/3W+emmnbTIrS8cPpHEmg7lCZDb327j8vWvf51PfvKTPP7443z+85/nX//rf02nc/mk6O1vfztvectbeOqpp5ibm+PDH/7wUN7fGJcJIRQmtHPrCkavSSQ9R97zcKYoaNv22lQjQlL6D6qHxtmpDVZi4WiokDzElOTjJ5bZZyW52u7fWMLFEgJWj3V+utHkHJaW8da46UwqiaZDdQIUkd/xjnfw7W9/mxe/+MW85S1v4Yd+6Ic4evQoN998M5/97GcB+O53v8tNN920cs6v//qv8573vKen63/2s5/l9ttvJ5lMcu2113L99dfzd3/3d5cco6p88Ytf5HWv83V/f+ZnfobPfOYzQ/l8JltsQnAaReKqpPq9gURj2AruFCkjl+q+oUxF+p+d2la4TDjYxsErckkO7c4MtTLlwskiRw7ODbz+Hy55hkug/SCZAlFR3FIRqzAd+0JSbZdyvD8PfxS8733v47HHHuORRx6h3W5TrVaxLItz587x0pe+lNe85jUbnv+BD3yAj33sY5e1f//3fz8f/OAHOX36NC996UtX2g8cOMDp05dqAp8/f55CoUAsFlv3mEGZOOMiIj+KXxDMAj6sqn+5ug2/Jsx7gPPAf1PV+0bS2SHiNl3sjoek+p+R2UQH0tMaV5yKH2PIxPsXJkzlriThbW2Z8Nj8HA9+a7C9Mqs5X25w4nyV19/an9JAN07DN3R9ye0HRAKJlPLy0tQYl4xXopX43r7O+dU/e5yvPzPcrLnDV1u8+5++oKdjVZVf+qVf4sEHHyQSiXD69GnOnNk4o/Ftb3sbb3vb24bR1W1hR42LiHwEeDVwVlVv6mq/DfhtIAr8kaq+b71rqOpngM+IyBzw68BfrtH2deB3VPWvReR+YPKNS6uCjcAANzNLYjhTJP9SDBSR04kBYhSJPLbn4TYHv5EcOTTHpxdOc/JCjfnd/SkmrGZFrPLQ4FIlbtMloUoqSLPuh3jONy7TpIyc1zLt5GQpIn/sYx9jaWmJhx9+mHg8zjXXXEO9XicWi+F5F5W4u9OBN/Nc9u/fz8mTF+swnjp1iv37L/Vud+/eTbFYpN1uE4vF1jxmUHbac7kX+F3gj8MGEYkCHwJege9xfCUwCFHgvavO/9muUse/HJzXTdh2Cni3iLwGmIrpmNOuYdP/HggAK5pkUadnF/KZom9c8qkBqodGIljqG+tBCXfRH19c3rpxOblMLCLcdPXgN0O3Wcb2QCL9h1CTgXGpOdORTaitOmkaaKo/Y92rhzFM8vk8pZKfTOE4Dnv37iUej/OlL32JEydOALBv3z7Onj3L+fPnyeVyfO5zn+O2224DNvdcXvOa1/D617+et771rTzzzDN861vf4tZbb73kGBHhB37gB7jvvvu4/fbb+ehHP8prX/vaoXy+HQ3oq+qDwOrF7luBp1T1aVVtAp8EXquqX1PVV696nBWf9wN/oarHAVa3qepZVf0F4B3AVEzJXK+BFek/9RbAjmVwdXoCtuddf7kgn+svLTvElhhOuzbw+9+4L082ER1K3GVhscjzr7JIJwabOMBgdX5C0rZvoFvl6TAu1aB8gKQLo+1ID+zevZuXvexl3HTTTTzyyCM89NBD3HzzzfzxH/8x3/u9/rJePB7nV37lV7j11lt5xStesdLeCy94wQv4iZ/4CQ4fPsxtt93Ghz70IaJR/3v2qle9imee8ZeX3//+9/Mbv/EbXH/99Zw/f54777xzKJ9vHGIu+4GTXc9PAS/Z4PhfBH4YsEXkelW9e3Ub8Hngl4As8IH1LiQibwTeCDA/P/ia907gaJvnRQdLr7TjOdzWWTz1iMjkJwheqPjZUbsG1EqzInGe20IBtVg0wosOFraskNzxlEdPFnndsa2V43U6dewBJx7hLvZpkd0vO+fIAtHsZMjtf/zjH9/0mDe/+c28+c1vHuj673rXu3jXu951WfsDD1yscHLdddddlkU2DMbBuPSFqn4Q+OBmbQRGY5Nr3QPcA3DLLbfosPq4HTh4fQsThliJPFoTSo0S9hRU5ysGadm75gabEFjRFE96W8ueO3Zojg996SkqjTbZPiVoQr55pkSl2Rl4f0uIoy2uig723cgHxiWUTJl0KkV/4hHPTsVq+EQzDtPY08DBrucHgjZDQMtrURWw4gN6Lkn/5uVWpkO8stRwyXoe9tyAy2KxDC5bK1d8dH4OT+HRU8WBrxEG8wfdmR/iaqfvOj8hqUyOmiaQKZHdb5R8I5nMT4bnMs2Mg3H5CnCDiFwrIgngduD+EfdprHCDf3x7wAwYO+3P4lx3Omx2qV0m3/GwrcE2ylnxHBXxjfaghAYhNBCDsLC4zO5sgvldW0sKcATs+GCeC0BJckSmRMGhWfJjLplNFJEN28+OGhcR+QTwZeBGETklIneqaht4E/AF4AngU6r6+E72a9wJhQntPjNgQqyMH7R1pkR2v9qpkvUEOz1ggkMiKKC2hdl6IZPge/ZktxTUP764zJH5wpb2yrRaVaoRweqzzk83lUie2BZSs8eJsKpmts86P4bhs6MxF1W9Y532B4AH1nrN0CVMmB7sH8YO9LQGkZkfR8o0yHgR4tHB5kah7L5bfobd2QHSmQOOHZrjL79+ZqDNlE61xbeXKvz40a0F810nmHgMoIgcUovmSbamQx7Iqy3jqWBPyYbQSWYclsUMm+BWQmHCwW6EdpBV5U6JMnKVNilv8HlRuEzolLa2THh0fo5itcXT5/rfM/NIEKs5crCwpT44Zf8zhKUEBqERt0i3p8NzkVqREhnSycSouzLzGOMyATiBgq89gDAhgBXK7tenIyOoLB1S2l8Fym7CAmruFj25Y1sQsTx+YpmIwAu3aFzCCqODTjwAWokCmS1mz40LkYZDSXITUafFSO4bRk4oTGgNIEwIkMjtI+15OFMQtFVVyhElKYMpCANYmXCZcGsF1L5nTw4rFRtov8vCySLP25fvu5LmakIDafdZ56ebTtLGmhLZ/XizSCWyeZGwccBI7htGjlv3b155a0DNn3gay1Pc5uTfQOqdOi0RUjJ4dlSoHuzWt7ZMGIkIR+bn+i577HnKI4vLW97fAuBUg4lHn3V+LiFVICMNmo3JlwhKtkrUYpNhXIzkvmHkOE2HfMcjmh58Xd1GcNqD62mNC05QuyQVHfwGEhppdwjLhEfn5/it//ZN3HoLK9Vb9trT58q49faKRtlWCA2kbQ2eGBDJBPuglpe44sqDmxw93qQ7Lsupwb24ncRI7htGjtMqY6tCdLDUWwiUkduTPzN13FARuTDwNaKZ3eQ7wymgduzQHKrwyGKR739eb3GP4yubJ4fguQSfIT/gkilclEqpFCffuGS9Eq1B0rL/4h3w3NeG25krb4ZXrivwfglGct8wEnxhwsGFDQHsSIIT3nACdaPkXNFPvc0lt5BqGktiqeI2tx7EftFBGxF4+MRyz8ZlYbGIlYpx3RWDL+2FuM0SeU+JxgbPjkrk/bGsuROeTahKXst4A5QeGDVGct8wEpxOHVu29qeyo2ncdnE4HRohS46fHZUboF58NxYRnCEUUMun4ty4L99XUH8hiLdEIlvPaHJaZSzd2nWSgXGplybbuLRrDjHxYJDNxj16GMPESO4bRo7rNbGjW8vbt2IZHNmantY4cL7sZ3jNDZiWHWJLDLcznGXCo4fmeGSxiOdtrn1abrR58kxpy3piIW6njiVb82qzgVRKqzzZqerlYiC3n9n6cuNOYCT3DSPHFybcmv6UnbBoNIV6u04qNnga76i5UPFvIHOFre1styJJzgxpmfDY/Bwf/1+LfOtsmRuv3DjR4Ksni6gOJ94C4HQa2JGtTTxCZeTOhMvul4tLFID4hMjtw3RL7hvPZcxRVRzRLQkTAliB6KVbnezaaaXGMlFVdm/RuNjR1NAKqIUlinvRGQuXz168xc2TIY62saKDbygFyNm78VTQ2taLn42SsFBYwigijwUTZ1xE5EdF5A9F5E9F5B8Hbf9QRO4WkT8Skb8VketE5MMict+o+7tVKs0yHRHsLQgTwkU9Lcc9NYRejY5Sq4TleVhzWxMmtOJZHPFQ3XoZn2t2Z9iVTfQUd1lYLHL93tzAopurcWXwOj8h0ViMkmSI1ItD6dOoaJb8iVPSMqKV48BOqyJ/RETOishjq9pvE5EnReQpEXnHRtdQ1c+o6s8DdwE/GbT9tareBXwO+GhQMnk4C4cjxgl0xawtCBMC2OlAGXnCa7pU2hVyHShkt7a0ZyfytEWobaHccYiIcHS+sKkMjKqycLK4ZT2xlet5Hu4W6vx0U5Yc0QlXcAhjRlnbiFaOAzvtudwL3NbdICJR4EPAK4HDwB0iclhEbhaRz616dKcI/XJwXjevBzZfxJwg3EBc0drCBkro0tMqT7YyclnrZD0hs4Wa89C1TLhFCZiQo4fmePpchQuV9eM4J85XuVBpriyjbZVa9RxtkYHr/HRTjVjEW5MtXulVfeOeKwyus2YYHjtqXFT1QWB11PBW4KnA22gCnwReq6pfU9VXr3qcFZ/3A3+hqsfDi4jIPOCoTolIUoAzBGFCADuQB3Gqw7mZjoqqNkl50S0LE9qp4SgjhxwNAvQLGyyNLZz0XxtWppjj+nsYBq3z0009lic14cZFq8s0NIad39oSsmE4jEPMZT9wsuv5qaBtPX4R+GHgdSJyV1f7ncB/AhCR3SJyN3BERN653oVE5I0i8pCIPLS0tDTwB9hOwhosVnaL+zqmRHa/Im1SuvV4hZ3x1+XdIRVQe9GBAtGIbBh3WVgskkvGuGHvcLSvnHJQ5ye19WWgRsIm05ls4yKNIi45YrGtebWG4TAOxqUvVPWDqnpMVe9S1bu72t+tqn8b/H4+eP17VPW9G1zrHlW9RVVv2bNnPF3pMLvLzl29pevk8vuJqk68MnI54m1Jbj/ECjZhhjGtrZJORDl8lbVhxtjxxWVedNAmOoTNk9CtiLy1iQdAJ2GT1cmW3Y81HMqRrcefdgojub/9nAa6BY0OBG0GLtZg2YowIYCkClieh9uY3Nlpx+tQEUhGtrbnB8DO+8uEw0zNPnZojkdPOrQ7l29WrTU7PPFsiSMHh7fBLywZYG1Bbj/ESxawtIx6k7vRNtFyqUYnZ0nMSO5vP18BbhCRa0UkAdwO3D/iPo0NbqNIwlNSW9yRTiSCpYLTmtzZablZQkVIRbauyWUFnqBTG97GwSPzBWqtDt947vKw39dOO3Q8HVq8BS4ucdpbEK1cITNHTDwq5cktd5xqO9Rjk2NcjOT+EBGRTwAvB64QkVPAu1X1wyLyJuALQBT4iKo+vpP9GmfcZgnb8yA+hNk6EdzO1vW0RoUbxBgysa1nR6VzVxFTxW0Wt3ytkDCof3xxmZv2X9rHhSFvngRwgjo/lr11JeNYIJlSWj5LzpoM+ZTVpDtlWunrR92NnjGS+0NEVe9Yp/0B4IG1Xpt1nFYFG4EhlG21JM5ypzGEXo2GC46/ATST2PoObEnZWJ6HM8RlwgNzafbmkzx8Ypmf/vvXXPLa8cVlrtmdYXdu6/GiEKfpEFMlndr6eESzflJApXgODt245euNgpyWaA+Ylv3+v3s/37jwjaH253t3fS9vv/XtPR1rJPcNO47TqWEN6c9kRVN8Vye3YNjZoj+jyg0hO4pIBFsZ6jKhv5ly7rKMMVXl+GKR77t+uDvH3VYZ2wOJbH11O2lNtjKytpvkqKFb3Gw8KozkvmHHcb0GV21RmDDEjqUnutTxkhOoFWS3GH8KsIniDmGHfjfHDs3x+cef42ypzt68ryLwjFNnqdQYarwFwGlVsYYUNk0HelzNCTUutdIFMgDpwZb0evUwhomR3DeMFMdrY29RmDDEjucoidLxhiPYuNNcqPh7kebsraVlh1iRBI433GXCo4cKABw/UVxpC2Vhjg5JCTnE9erYkeHMD7Nzfip+uzKZ4pXlov/diGYnJ15kJPcNI2UYisghdtJGa0K54WIPOMMbJW6Q2bW7MJxSvHYkybeHvEz4gqttEtEIxxeXue0m38NaWCySikc2lePvF9drsTeaHsq1rEAypVOdTONScXyPK56bLF0xI7lvGAmtTouagBUfzk1pRRl5SJInO43bdEl7HnO7rhrK9axYBpfh7utIxaO8YL91iYjlwsllXri/QDw63H83RztYQzIuqXSWhsZhiKnZO0ndDRSRc0Zuf1wwxmWMccJ9DEMQJgSw0/6sLkzpnTTK7RKWp9jWkMYjkaMkDH2Z8Nj8HF897dBsezTaHR4/7XIkWC4bJu4QvVqJRHAlR7QxmftcmoEicto2cvvjgjEuY4wb1F4ZhjAhdMnuD0lPa6cpd2pkO2ClhpQ9F9TIKQ35hnr00BzNtsfjzzg8/oxLs+MNdWc+QLtVpxwRrMTwltoqkRyxCTUunbI/EcsWjHEZF0zMZYxZESZMD2cdOcyymlRl5Ko2SHsRYkNaXrICo+2UTlPYYkmDbo4dCjdTFlfajg45U6wUlmIYYuptLWqRaE+mPJAXVNG05vrTCFTVLStszwKDFNUznssYsyJMmNm6dhR0KSNPaKnjCi3SOrz5UCi775aGu0y4z0qxv5Dm+Illji8us7+QZq+1teJmqwkritpDmngANOIW6Qk1LtSKlDRNNtV7ZmUqleL8+fNDqUY6zagq58+fJ5Xq7ztsPJcxJvQwwlosWyUUvwxlQyaNSqRDUocTwAawMmF1zuEoI3dz9NAcD333AhGRoe9vga4KpenhLQM14zbZ6lNDu95OEmkUKUmOfB9eyIEDBzh16hTjWm5jnEilUhw40J94rjEuY4xT9deRrfxwjEs8u4e0502s7H5ZlBTD8wDsQAzU3YZlwqPzBf7sUd8j+tnvu3bo178otz8crxZ8ZeTchMrux5sulUh/8ad4PM611w7/b2PwMctiY4zbWEZUyVnDkWMglsT2FHcClZEbnQaNyHAUkUOsleqcw18mPNZVynjY8RYAp+rPtq2tqmV3oakCOanRak6e/lyy5VCLDncfkWFrGOMyxjgNl7znER1isNkmgtOaPGXkME6UHuINxAqWCd1tWCZ8/lUWqXiERDTC4auHLwPvDqnOTzcSKiMXJy8ml+qUaMQnR25/FjDGZYxxWiUsBSIbl22tNtv8u089ytlSfcPjAGyJ4XaGq6e1ExSDAHYqPpw9LgDxzBVkPQ9nG9Jv49EIt167m2OH5khuQ9ldp14EIJ8fjhQOQCzYgBhKqUwS2U6J1hC/G4atIyZTwkdEloATwVMb6L7jbPb8CmC7pnur32uY52x03HqvrdVuxmtr47W6zYxX/23dz8147dx4HVLVtfO/VdU8Vj2Ae/p8/tBO9WWY52x03HqvrdVuxmtr47W6zYzX1r5zZrzGY7zMstja/Fmfz7eTQd6r13M2Om6919ZqN+O1tfFa3WbGq/+2nRozM149YpbFhoCIPKSqt4y6H5OCGa/+MOPVH2a8+mO7xst4LsPhnlF3YMIw49UfZrz6w4xXf2zLeBnPxWAwGAxDx3guBoPBYBg6xrgYDAaDYegY42IwGAyGoWOMyzYiIs8XkbtF5D4R+Vej7s8kICI/KiJ/KCJ/KiL/eNT9GXdE5DoR+bCI3DfqvowrIpIVkY8G36ufGnV/xp1hfaeMcVkHEfmIiJwVkcdWtd8mIk+KyFMi8o6NrqGqT6jqXcBPAC/bzv6OA0Mas8+o6s8DdwE/uZ39HTVDGq+nVfXO7e3p+NHn2P04cF/wvXrNjnd2DOhnvIb1nTLGZX3uBW7rbhCRKPAh4JXAYeAOETksIjeLyOdWPfYG57wG+HPggZ3t/ki4lyGMWcAvB+dNM/cyvPGaNe6lx7EDDgAng8M6O9jHceJeeh+voWDquayDqj4oItesar4VeEpVnwYQkU8Cr1XV9wKvXuc69wP3i8ifAx/fxi6PnGGMmfg1Z98H/IWqHt/mLo+UYX3HZpF+xg44hW9gHmFGJ9R9jtfXh/GeMznQW2A/F2dA4H9p1y22IiIvF5EPisgfMBuey1r0NWbALwI/DLxORO7azo6NKf1+x3aLyN3AERF553Z3bsxZb+w+DfwzEfl9dlZaZ9xZc7yG9Z0ynss2oqp/BfzViLsxUajqB4EPjrofk4KqnsePTxnWQVUrwL8cdT8mhWF9p4zn0h+ngYNdzw8EbYb1MWPWH2a8BseMXX9s63gZ49IfXwFuEJFrRSQB3A7cP+I+jTtmzPrDjNfgmLHrj20dL2Nc1kFEPgF8GbhRRE6JyJ2q2gbeBHwBeAL4lKo+Psp+jhNmzPrDjNfgmLHrj1GMlxGuNBgMBsPQMZ6LwWAwGIaOMS4Gg8FgGDrGuBgMBoNh6BjjYjAYDIahY4yLwWAwGIaOMS4Gg8FgGDrGuBgMBoNh6BjjYjCMABG5ZnVtjTWOebmIOCLyQPD8+0TkYRF5XEQ+KyLJoP1LIlIWkVt2ou8GQy8Y42IwjDd/raqvCn6vA69U1RcAVeCfA6jqDwAPjah/BsOaGONiMIwY8cvKLojI39voOFV9SFXPBk+T+MbGYBhLjOS+wTBCRORG4JPAG1T10R7PuRO4EvjsdvbNYNgKxrgYDKNjD76B+HFV7an6n4jsAd4NHFXV1nZ2zmDYCmZZzGAYHQ6wCHxfH+fcCHxNVc9tT5cMhuFgPBeDYXQ0gR8DviAiZVX9eA/nfBN43/Z2y2DYOsZzMRhGSFCC99XAW0TkNT2cMk+QJWYwjDPGczEYRoCqfhe4Kfi9CGyYKdZ13kOYtGPDBGA8F4NhfGkCN4WbKNdDRL4EXAeYAL9hbDCVKA0Gg8EwdIznYjAYDIahY4yLwWAwGIaOMS4Gg8FgGDrGuBgMBoNh6BjjYjAYDIah8/8DMMEJB8n2On4AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "tau_indices = [0, 250, -1]\n", + "\n", + "fig, ax = plt.subplots(1, 1, figsize=(6, 3))\n", + "ax.set_title(f'neutrino perturbations')\n", + "ax.set(xlabel='k [?]', ylabel=r'$\\delta_\\nu$', xscale='log', yscale='log')\n", + "for tau_idx in tau_indices:\n", + " delta_nu_k = [delta_nu_list[k_idx][tau_idx] for k_idx in range(len(klist))]\n", + " ax.plot(klist, np.abs(delta_nu_k), label=f'tau={tau_list[0][tau_idx]}')\n", + "ax.legend()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c9a3be8", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.10.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/source/background.c b/source/background.c index 22246e23e..9bfb7cbf9 100644 --- a/source/background.c +++ b/source/background.c @@ -493,6 +493,7 @@ int background_functions( if (pba->has_ncdm == _TRUE_) { /* Loop over species: */ + int cumulative_q_size = 0; for (n_ncdm=0; n_ncdmN_ncdm; n_ncdm++) { /* function returning background ncdm[n_ncdm] quantities (only @@ -511,6 +512,14 @@ int background_functions( &pseudo_p_ncdm), pba->error_message, pba->error_message); + + if (pba->output_ncdm_binning == _TRUE_) { + for (int index_q = 0; index_q < pba->q_size_ncdm_bg[n_ncdm]; index_q++) { + pvecback[pba->index_bg_q_ncdm1 + cumulative_q_size + index_q] = pba->q_ncdm_bg[n_ncdm][index_q]; + pvecback[pba->index_bg_w_ncdm1 + cumulative_q_size + index_q] = pba->w_ncdm_bg[n_ncdm][index_q]; + } + cumulative_q_size += pba->q_size_ncdm_bg[n_ncdm]; + } pvecback[pba->index_bg_rho_ncdm1+n_ncdm] = rho_ncdm; rho_tot += rho_ncdm; @@ -1052,6 +1061,12 @@ int background_indices( class_define_index(pba->index_bg_rho_ncdm1,pba->has_ncdm,index_bg,pba->N_ncdm); class_define_index(pba->index_bg_p_ncdm1,pba->has_ncdm,index_bg,pba->N_ncdm); class_define_index(pba->index_bg_pseudo_p_ncdm1,pba->has_ncdm,index_bg,pba->N_ncdm); + int total_q_size = 0; + for (int n = 0; n < pba->N_ncdm; n++) { + total_q_size += pba->q_size_ncdm_bg[n]; + } + class_define_index(pba->index_bg_q_ncdm1,pba->output_ncdm_binning,index_bg,total_q_size); + class_define_index(pba->index_bg_w_ncdm1,pba->output_ncdm_binning,index_bg,total_q_size); /* - index for dcdm */ class_define_index(pba->index_bg_rho_dcdm,pba->has_dcdm,index_bg,1); @@ -1258,7 +1273,21 @@ int background_ncdm_distribution( } /** - b) deal now with case of reading analytical function */ - else{ + else if (pba->collective_ncdm == _TRUE_) { + /* The collective homogeneous distribution of all species */ + *f0 = 0.; + for (int n = 0; n < pba->collective_ncdm_N; n++) { + // Add the contribution of the n'th species + // The first species is the reference species + // Remember: Here, q is actually q/T + double mass_ratio = pba->m_ncdm_in_eV[n]/pba->m_ncdm_in_eV[0]; + double temp_ratio = pba->T_ncdm[n]/pba->T_ncdm[0]; + double deg_ratio = pba->deg_ncdm[n]/pba->deg_ncdm[0]; + *f0 += deg_ratio*pow(mass_ratio/temp_ratio, 4.)*(1./(exp(mass_ratio/temp_ratio*q - ksi) + 1.) + 1./(exp(mass_ratio/temp_ratio*q + ksi) + 1.)); + } + *f0 *= 1.0/pow(2*_PI_,3); + } + else { /** Next enter your analytic expression(s) for the p.s.d.'s. If you need different p.s.d.'s for different species, put each @@ -2444,6 +2473,14 @@ int background_output_titles( class_store_columntitle(titles,tmp,_TRUE_); sprintf(tmp,"(.)p_ncdm[%d]",n); class_store_columntitle(titles,tmp,_TRUE_); + if (pba->output_ncdm_binning == _TRUE_) { + for (int index_q = 0; index_q < pba->q_size_ncdm_bg[n]; index_q++) { + sprintf(tmp,"q_ncdm[%d][%d]",n,index_q); + class_store_columntitle(titles,tmp,_TRUE_); + sprintf(tmp,"w_ncdm[%d][%d]",n,index_q); + class_store_columntitle(titles,tmp,_TRUE_); + } + } } } class_store_columntitle(titles,"(.)rho_lambda",pba->has_lambda); @@ -2514,9 +2551,18 @@ int background_output_data( class_store_double(dataptr,pvecback[pba->index_bg_rho_cdm],pba->has_cdm,storeidx); class_store_double(dataptr,pvecback[pba->index_bg_rho_idm],pba->has_idm,storeidx); if (pba->has_ncdm == _TRUE_) { + int cumulative_size = 0; for (n=0; nN_ncdm; n++) { class_store_double(dataptr,pvecback[pba->index_bg_rho_ncdm1+n],_TRUE_,storeidx); class_store_double(dataptr,pvecback[pba->index_bg_p_ncdm1+n],_TRUE_,storeidx); + if (pba->output_ncdm_binning == _TRUE_) { + int q_size = pba->q_size_ncdm_bg[n]; + for (int index_q = 0; index_q < pba->q_size_ncdm_bg[n]; index_q++) { + class_store_double(dataptr,pvecback[pba->index_bg_q_ncdm1 + cumulative_size + index_q],_TRUE_,storeidx); + class_store_double(dataptr,pvecback[pba->index_bg_w_ncdm1 + cumulative_size + index_q],_TRUE_,storeidx); + } + cumulative_size += q_size; + } } } class_store_double(dataptr,pvecback[pba->index_bg_rho_lambda],pba->has_lambda,storeidx); diff --git a/source/input.c b/source/input.c index c376d1d8a..e9c7bcee2 100644 --- a/source/input.c +++ b/source/input.c @@ -2536,6 +2536,25 @@ int input_read_parameters_species(struct file_content * pfc, /** 5.a) Number of non-cold relics */ /* Read */ class_read_int("N_ncdm",N_ncdm); + int collective_ncdm = 0; + pba->collective_ncdm_N = 0; + class_read_int("collective_ncdm", collective_ncdm); + if (collective_ncdm != 0) { + pba->collective_ncdm = _TRUE_; + } + else { + pba->collective_ncdm = _FALSE_; + } + + int output_ncdm_binning = 0; + class_read_int("output_ncdm_binning", output_ncdm_binning); + if (output_ncdm_binning != 0) { + pba->output_ncdm_binning = _TRUE_; + } + else { + pba->output_ncdm_binning = _FALSE_; + } + /* Complete set of parameters */ if (N_ncdm > 0){ pba->N_ncdm = N_ncdm; @@ -2556,6 +2575,9 @@ int input_read_parameters_species(struct file_content * pfc, } } if (fileentries > 0) { + if (pba->collective_ncdm == _TRUE_) { + class_test(_TRUE_, errmsg, "You cannot use psd files with collective_ncdm."); + } /** 5.b.1) Check if filenames for interpolation tables are given */ /* Read */ @@ -2582,6 +2604,7 @@ int input_read_parameters_species(struct file_content * pfc, class_read_list_of_doubles_or_default("Omega_ncdm",pba->Omega0_ncdm,0.0,N_ncdm); class_read_list_of_doubles_or_default("omega_ncdm",pba->M_ncdm,0.0,N_ncdm); for (n=0; ncollective_ncdm == _TRUE_) && (pba->m_ncdm_in_eV[n] == 0)), errmsg, "You must input m_ncdm when using the collective_ncdm feature.") if (pba->M_ncdm[n]!=0.0){ /* Test */ class_test(pba->Omega0_ncdm[n]!=0,errmsg, @@ -2589,6 +2612,7 @@ int input_read_parameters_species(struct file_content * pfc, /* Complete set of parameters */ pba->Omega0_ncdm[n] = pba->M_ncdm[n]/pba->h/pba->h; } + /* Set default value this is the right place for passing the default value of the mass (all parameters must have a default value; most of them are defined @@ -2643,6 +2667,14 @@ int input_read_parameters_species(struct file_content * pfc, else { class_read_list_of_integers_or_default("ncdm_N_momentum_bins", pba->ncdm_input_q_size, 150, N_ncdm); } + + // Collapse all species into a single species + if (pba->collective_ncdm == _TRUE_) { + // Set N_ncdm back to 1 to only use one hierarchy + pba->collective_ncdm_N = N_ncdm; + N_ncdm = 1; + pba->N_ncdm = 1; + } /** Last step of 5) (i.e. NCDM) -- Calculate the masses and momenta */ class_call(background_ncdm_init(ppr,pba), @@ -2681,7 +2713,7 @@ int input_read_parameters_species(struct file_content * pfc, pba->deg_ncdm[n] *=fnu_factor; } } - else{ + else { /* Case of only Omega/omega: */ class_call(background_ncdm_M_from_Omega(ppr,pba,n), pba->error_message,