DiffEq_GrowthFunction.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "7a1f5668",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import solve_ivp # as suggested by Nikhil\n",
    "import warnings; warnings.filterwarnings('ignore')\n",
    "from scipy.interpolate import UnivariateSpline\n",
    "plt.rc('axes', lw =2); plt.rc('font',family = 'serif'); plt.rc('text', usetex = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "id": "c08ab9af",
   "metadata": {},
   "outputs": [],
   "source": [
    "global OmegaM0\n",
    "OmegaM0 = 0.3175\n",
    "\n",
    "global OmegaDE0\n",
    "OmegaDE0 = 1- 0.3175"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "456e362a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Omega_DE(a):\n",
    "    numerator = OmegaDE0\n",
    "    denominator = (OmegaM0*a**(-3) + OmegaDE0)\n",
    "    return numerator/denominator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "4cdcd088",
   "metadata": {},
   "outputs": [],
   "source": [
    "def diffEq(ln_a, sols):\n",
    "    a = np.exp(ln_a)\n",
    "    g, gdot = sols\n",
    "    gddot = -(2.5 + 1.5 * Omega_DE(a)) * gdot - (3 * Omega_DE(a) * g) \n",
    "    return [gdot, gddot]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "id": "526fa56c",
   "metadata": {},
   "outputs": [],
   "source": [
    "ln_a_vals = np.linspace(-5,5,5000)\n",
    "a_vals = np.exp(ln_a_vals)\n",
    "z_vals = 1/a_vals - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "id": "93dcdcae",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = solve_ivp(fun = diffEq, t_span = [-5,5], y0 = [1,0], t_eval = ln_a_vals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "id": "5221cb2c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g at a=1 is 0.78968\n"
     ]
    }
   ],
   "source": [
    "calc_g = res.y[0][np.where(np.abs(a_vals - 1) < 1e-3)]\n",
    "print('g at a=1 is %.5f'%calc_g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "id": "2c391d67",
   "metadata": {},
   "outputs": [],
   "source": [
    "normalized_D = a_vals*res.y[0]/calc_g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "id": "a828a200",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAF0CAYAAAAgtIKnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/+UlEQVR4nO3dd3xb9b3/8deR7cSxnURxljMIRCYBCrTUdhpooYPI3e3tba2ke1yIXdre29uBTdre0RlkOn7d2Omg9LYlsbroJBZQCrRAZEEppRCwwshwhmXFiUfiod8fX0uxvOIhWbLP+/l4nEdi6Vj6Skre56vv+Z7P14pGo1FERGTWc6S7ASIiMj0U+CIiNqHAFxGxCQW+iIhNKPBFRGxCgS8iYhMKfBERm1Dgi4jYhAJfRMQmMjLwfT7fWfcJBoPU1tZSW1uLx+MhFApNQ8tERGau7HQ3YDCfz0c4HKaqqoqxKj5EIhH8fj/V1dXx3ysvL6e5uXm6mioiMuNYmVhLx7KsMQPf7/dTXl4e3ycSibBo0SKam5txuVzD9u/o6Bjz+fLz86fWYBGRGSCjevjj5Xa7aWpqiv8cCAQARgz7/Px8Ojs7x3y8VatWJbeBIiIpdOTIEfr6+sjNzT1rh3awGdnDH6q8vByPx0NlZeWw+7Kysujv709m80REMoLD4aCvr2/c+8/4wK+vrwcYMewB5syZQ09PDw6Hg6KiolGfT0Rkpjh06BD9/f3k5ORw+vTpcf/ejBzSifH7/RQWFlJRUTHqPsuWLePAgQOsWLGC/fv3T2PrRERSY/Xq1Rw4cIBly5ZN6PcyclrmSILBYMLUy2AwCBAP+/r6eiKRSDqaJiIyI2RU4Pv9fmprawGoqanB7/fH79u+fXt8fn4oFKK0tJTy8nIsy8KyLGpqanA6nelotojIjJCRY/jJFPvqs2rVKg3piMisMNlcm9Fj+JJZ+vr66OnpSXczRGaFnJwcsrKykvqYCnyZsmg0SktLi86hiCSZ0+mkqKgoaTMJFfgyZbGwX7ZsGXl5eZrmKjJF0WiUzs5Ojhw5AsCKFSuS8rgKfJmSvr6+eNgvXrw43c0RmTXmzZsHmKtqly1blpThnYyapSMzT2zMPi8vL80tEZl9Yv+vknVuTIEvSaFhHJHkS/b/K9sEfktLCy6Xa8QCayIidmCbwBeR6XO2BYk0oys9bBP4RUVFhEIhrYwlMsDn81FTU5P08K2trT3rN+ldu3Yp9NNAs3Qkuf7xJegff7nWlHFkwcWfGvXuYDDIzp07qa2tpa6uLl5tNRQK4fV62bVrF16vd9QqrNPB7/dTVVVFVVVVfHW30tJStm3bNmbBwPGqqKhg+/btbNmyhZKSkmH3B4NB6urqqK+vp7q6muLiYiKRCK2trRQXF4/43vh8Ptxu91mfu7KykpqaGrxe77jbO7g9lZWVVFVVjdjukX5vPPsly0ifW6ZQ4Ety9fdBNAMC/yxLIJSUlOByuVi8eDFVVVW43e74OZ66ujpKS0vTGvZgFvqpqqpKuM3r9VJWVpa05ygsLBz1vpKSErxeL/X19Wzbti2hVlVNTQ1VVVXU1dUl/M7OnTtpaGgY13MXFxdPKIxLSkqoqamhvr4er9c7rtpZoVCIQCAwrYE/0ueWKWwzpCMykliolZeXp7sp4+J2uzOiSGDsW9DgAofBYJANGzaM+zHcbjc7d+6c0POOdYAaKhKJUFNTM6HHn+3Uwxfbq66uprGxccQeK5ggCwQCuFwugsEgFRUVuFyuhK/uTqeTuro6vF5vvPcbG8eO3R4r6d3Y2JjwPLF1HQKBAM3NzaMOcwSDQbZu3UpVVRWVlZX4/X7q6urYsmULTqcTj8dDQ0MDbrcbv99PKBSisLCQPXv2xHvowWAQv98fb1s4HJ70+1ZWVhZ/vtjrGNyTHqt9QPw9nAq/3z/i+93Q0EAgECASidDY2AgQ/xY32nvj9/sJBoM4nU6am5vZsGEDe/bsoby8fNjn3NTUNO7PLZMo8EWAhoYG1q5di8fjSRiDDoVC1NTUJIRGaWkpd955J263m4qKCvbs2UNDQwMulwu3282WLVtobGyMD214vV58Pl98PDd2ojTWU/d4PDQ1NcXHpX0+34hj9CUlJWzZsiV+sjMSicSfo6amhs2bN+N2uwkGg3i93nibwZQX37ZtGx6Ph+bm5oTbJ8vpdCZMgmhtbU3ogY/WvmQa7f0OBoO43W4aGxsTzjeM9t54vV48Hg9tbW2AGW6qqqqKfw5DP2cY/+eWSRT4Ipjw2rFjBx6Ph3379sVv9/l8w8Z/XS4Xu3btiodIbBhjcJgNHtpwOp0Js1acTifhcDge+E1NTQSDQYLBIOFweNwzyQb3rH0+XzzId+7cGe+xxixevJhdu3YNey0TGSIZKhKJJLyuobNuRmtfKgx9v0f75jLaezO07bGD2eDXN/Rznuznlk4KfEkuR9ZZT5hOC8fE645UVFTQ2NiIx+PB4/EAptd6NuO5mG+0cfdIJBIf6nC5XBMKDafTSSQSoaqqKt5jjQ0bFRYWDutNx9Z/TpZAIJAwRj70NY7Wvuk8gTrYWO8NmM8x9s0rEokM22fowW2yn1s6KfAlucaYCplpRuoF1tXVxacfAlRVVcXDPyYYDLJjx44xH2eo0eacx8aBY2HS3NxMcXHxuIcHPB4PNTU18d8PBAJUVVUNOwnt9/vZvHnzsHHmyQZVbW0tbrd72LeaoTNiRmrf4PsHf8MIBoO4XK6knpQeXNAvFAqN+t7EhoZ27doFkDDkEzP4c57q55YuCnyxpWAwSE1NDeFwGK/XmxBcDQ0N8a/8LpcLr9cbv5goNo4bGxaIneiLjd/H5veD+eofCoXiPcuSkpL4/l6vN/68O3fujC/f6fF4qKuri58gHvxYQPznioqK+DrPhYWF+Hw+6urq4icv6+rqqKmpSRiGcDqdNDQ0UFtbS0lJSbw3GxvDHvpNJTbvHcw4d+xA2NzcTGlp6bDplxUVFQnTEX0+34jti/H7/QkH0+3bt1NYWDjiifOR2rNly5aE92Tw+x17DysqKuJTS2MnbUd6bwD27NkTP2nb0NCAx+OJnxwf+jlP5HNL1zeaEUVnuVWrVkWB6KpVq9LdlFmpq6sr+vjjj0e7urrS3RTJAJWVlRPat62tLeG2hoaGJLdofBobG6Nerzf+c3Nzc7SysjJt7YkZ7f/XZHPNNvPwVTxNJPVi34bOJhgMUl5ePmz4ZirTRKciNqsnxuVyUV5enrb2pIqGdEQkaZxOZ3xoa6yhjEAgMOxK5th5hnSorq6mtrY24ZxGOBxO+9XWyWZFo9FouhuRSpNd3V3Gp7u7m3379rF27Vpyc3PT3RyRWWW0/1+TzTXbDOmIiNidAl9ExCYU+CIiNqHAFxGxCQW+iIhNKPBFRGxCgS8iSadFzDOTAl/EprSIuf0o8MWWYsXTLMtKKBscq6i4aNGipJcTnii/309xcXFCqYLS0tJ4wa6pqqioiK/+NJJgMEhVVRWWZcXXkq2trY3/fSQTWcR8vIuvxJ4vXevEBoNBamtrqa+vp76+Hp/Pl7LPJNVUWkGSLhqN0tmZnqL4eXkOLMs6635axNzI9EXMYweWdJVcABJWPAPTpsHrJCT7M0kl2wR+rHgaTL4GuIxPZ2c/BQV/S8tznzz5IvLzx7/4yeBFzFO5IlOyJHuJwMnyer0sWrQoYUnIySxiXldXN2bgRyIRiouLcTqdaTkAh0KhYXlRWVmZ0MPPlM9kPDSkI7ZXXV2Ny+UadcggGAxSX1+P3+9PKLA1eMilvr6e0tJS/H4/paWl8f1jNd9DoRA+nw+fzzfseWK11uvr6xNWkBqpHbHHjv2ex+PB5/Ph9/tZtGhRvI6/3++PDz8MHqePDU/E2pKMRcwHv46hi5iP1j44+yLmwWAwvv5A7P0c+n7H9hvt85noZzFUbEGW2LBN7LFj6xOP9JlM9TlTKmmFmzOU6uGn1kj1uvv7+6MnT/amZevv7x9329va2qKNjY3xvzudzvjPdXV10WjU1EV3u90Jv1dSUhKv415dXR2tqKiIRqPR+O96vd74bdFoNOp2uxNqrbtcroQ68E6nM9rc3ByNRqPDarB7vd6E3x388+D9qqur47Xom5qaEtrc0NAQra6ujra1tUVdLtew19LU1DTmewQMq1sfjUajFRUVCc9TXV2d8FijtW/o849l6Osf+n6f7fOZ6Gcxmurq6qjL5YoCCY8/UhuT9ZzRqOrhywxgWRb5+Vlp2cYzfj+SwYuYD549MtYi5jGTXcQ8JrYYdqzHPdlFzGPj6YMX6o71oDNhEfPRVrKaqMHv90Q+Hzj7ZzEar9dLc3MzbW1tuFwutm7dOq42TuU5U8E2Y/giZ6NFzCcmXYuYD36/x/P5DDWRNXMHL08Z+90dO3bEh5NS8ZyppB6+2NZoi5iHQqGERcyHjjMHg8GEWSPJXsQcGPc0v9EWMR/a5tjiIrEAi0nFIuZna99gk/mGMfj9Hs/nM9RE5/8PnT4aCAQmfKJ2rOeMrU08HWzUw5/V67zIBGkR88xbxHykzyj2emMnT4e+32N9PhP9LEbrhW/YsCH+7SgSidDa2hp/XwY/R0VFBZFIZMLPuX37djZs2BA/EZxK9lnxqmgx+w8dS3dzZh2teCWDjTQ3f6x9xwpa0YpXkxftgd6udLdCZFab6iLmkloZGfjjGb8MhULx+cS1tbVnH5eLRiEcGHsfEZmSwYuYjyUQCFBRUTFNrZKYjBrDj01Lq6qq4mwjTR6Ph6amJuDMTIfBlz+PqPVBWHIFODLqZYvMKuOZhZPushV2lVE9/IqKinH9Qxjae3A6nQQCgbP28k93huk69CAdHR0Jm4iIHWRU4I9XIBAYNp2rsLBw2JSvwdoiXXzpS9u5/gNXUlBQkLDJ1M3yc/8iaZHs/1czMvBH6snHLvIYzenePACWzIf1K1LUMBvKyckBoLOzM80tEZl9Yv+vYv/PpmpGDmaPdGY/Nqd4NL39ubzjfdezdlUn2+adQ++5701dA20kKysLp9PJkSNHAMjLy5t0eQMRMaLRKJ2dnRw5cgSn00lW1vgrwI5lRgZ+WVnZsLm+4XD4rJe4f/u2C/nGDY9B72HmEob8c1LZTNsoKioCiIe+iCSH0+mM//9KhhkT+MFgMF6EKHaVYEysiNPZAv/7v1zD/3zwSRY7e+DYXyB/S4pbbQ+WZbFixQqWLVtGT09PupsjMivk5OQkrWcfk1GBH7vsGMwqM+Xl5fFLyodefhy7RHzw5dRjycnup7M7m+/uOo/PVD4F7U/AqVaYuzi1L8pGsrKykv4PVESSxzalFRYtcNLWfifLCk/x7B8byZ3bD4vLYNUb091EEZEJUWmFs8jL7WPNik6OhOfy49+uNje2PQI9J9PaLhGR6WKbwAf4z3eZEqRfubWY/n6gv9dcfSsiYgO2Cvxr3/osC+f38OQz8/nNPQNnvlsfgr7u9DZMRGQa2Crw5+f3cZ3nGQBuuqXY3Nh3Clr3pK9RIiLTxD6BP1Aw7d/fESInu5/7H1nMX/+2yNx37K/Qr+mEIjK72Sjw5wKwctkp3v0Gc1b7xh+cb+7r7YTw2OVcRURmOtsEfsvho2z7wrf4+je+Ts2/PYVlRbn9Tyv4+1PzzQ5H74f+vvQ2UkQkhWwT+AB7njFLhF1wXgcV5QcB+NL31pk7e9oh8mi6miYiknK2CfyioiL8Dz7PR6u/CMCnrnkKgF27V/HUs/lmp6P3QbQ/XU0UEUkp2wQ+AJYFS68E4LIL23nDVS3091tnxvJPtcLxf6axgSIiqWOvwAdwXgJzzOycT281vfxbf3sOzx2aZ+4/eq9Z/1ZEZJaxX+BbDlj6MgCueFEbr9pwlN5ex5l5+V0tcOLpNDZQRCQ17Bf4AIsugxwzOyfWy//eL8/lcKuZusmRe9TLF5FZx56B78iGJVcAcPVLjrHx0jDdp7L46q0D9fQ798PJ5jQ2UEQk+ewZ+ACFZZA9D8uCT19revnf2bWW8PGBtSMP/0m9fBGZVewb+FlzYPHlALzxFYd54frjnOzM5us/GdzLD6WxgSIiyWXfwAdYshGycrEs+MzWvQD8v5+4zvTyj/xJvXwRmTXsHfhZufGx/Le5D3HpunbaT+bw1R8PzNjpeF69fBGZNewd+BDv5Tsc8NnrngDg6z9xcaxtjrlfvXwRmSVsE/gtLS24XC5cLlfiHYN6+W+5uoUXXxjhZGc2X/6RevkiMrvYJvDHNGgs/7PXPQnAt25by5FW9fJFZPawTeAXFRURCoUIhUborQ/q5b/xFYcpu7iNjq5sbvrRQI0d9fJFZBawTeCf1Qi9/G/vPI+WYwNX3x6+W718EZnRFPgxWbmwxMzLf92VR9h4aZiu7mxqfzjQy+/cDyf2prGBIiJTo8AfbMnl8V7+5z5kevnfbTiPg0cGevktd6pevojMWAr8wQaN5ZdfcZSXXmZq7Hxhx3pzf/cRiDyWxgaKiEyeAn+oJZdDdh6WBV/8iFkMZccvzuXp5wZWxTp8t9a+FZEZSYE/VNZcWHoVAK/c0MprX3aY3l4H//XtC8z9p9ug7eE0NlBEZHIU+CNZvAFyFgCw/T9ML/+2P64m+M+F5v4j90B/T7paJyIyKQr8kTiyYfkrAbP27btevx+AbV+/yNzfcwJa96SpcSIik6PAH43zRTB3MQCf//AT5GT3s/uvy7jrwSXm/qP3Qt+pNDZQRGRiFPijcWTB8lcBsHZ1Jx/0PAPADV+/yFx/1dsFR+9PX/tERCbINoE/avG0sSy8GOYtB0y9/IK8Xvb8YxE/968w9x/7K/S0p6C1IiLJZ5vAnxTLguWbAFi2+DSfeK9Z5/bT37yI3l7LnLg9fHc6WygiMm62Cfwxi6eNZf46yD8HgE+892mWLjrF3mcL+N4v1pj72x6BrsPJbayISArYJvAnzbKg6NUAzM/v47+rTMmF//7OhRw/kW0KqrU0prOFIiLjosAfj/xzYOELAKiqeJYL157gaNtcvvS9deb+E0/DieY0NlBE5OwU+ONV5AYri5ycKF/++D8As+B5aH+eub9ltwqriUhGU+CP19xCcwUu8PqrjuC+/Cine7K4IXYxVtdhiDyaxgaKiIxNgT8Ry14eL5/81U8+hsMRpWH3Ku4LFpr7W+6CvtPpbaOIyCgU+BORnQfLTGG1S9ed4Np/fRaAj335Yvr7MXPyj+liLBHJTAr8iVq8EeY4Afjch59kfn4PgX8s4ie/W23uP3o/nI6krXkiIqPJuMAPhULU1tbi8/mora0lEomMuW99fT319fXU1NQQDAZT30BHNhSZi7GWLz7Fp655CoBt37iIjs4s6O+FQ7tT3w4RkQnKTncDhvJ4PDQ1NQEQiUTweDw0No48z72urg6v1xv/uaqqirq6utQ3cuElkPcgdO7nP98dos53Ls8czMf7w/P53IefhOOPw8l9ULA29W0RERmnjOrhD+2hO51OAoHAqL18n8837itno9EoHR0dI24TZlmw8nUA5M7t58ufeByA2lvOp/n5gWmaB/+gaZoiklEyKvADgQCFhYUJtxUWFhIIBEbcv6qqiuLiYmpra+PDOqM5ePAgBQUFI26TkrcKFl0GwFs3HaL88iOcOp3FR72XmPu7j0DryO0WEUmHjBrSGakn73Q6R+3hV1ZW0trays6dO4lEIpSVlaW2gUMVuaH9n1h9p/jmtr9z6dtexe/uLeI3f1rOm155GI7cDc5LzOweEZE0y6gevtPpHHZbJBIZ9faamhq8Xi9NTU3U1NSwadOmUQ8OK1eu5OTJkyNuk5ZTAMteAcAF53Xw8feY8gofrb2Erm6HqZl/+K7JP76ISBJlVOCXlZURDocTbguHwyPWsPf7/ZSXl8d/rqyspLKyctThH8uyyM/PH3GbksUb4ytjfaZyL6uXd7HvQD61t5w/8AKaoPPA1J5DRCQJMirwS0pKEnrokUgkYdGSYDAYP0nrcrnYs2f4urLTPqzjyIIVrwWgIK+Pr37yMQC2f3+dqbMTjcKB3+oEroikXUaN4QM0NDRQW1sbD/SGhob4fdu3b2fDhg1UV1dTUlISn7MfG+ffsmXLiMM/Kbdgnambf+IpKsoPsWnjUe58cCn/WXsJt3/jIeg6ZBY9X7Jx+tsmIjLAikaj0XQ3IpVWr17NgQMHWLVqFfv370/dE51qhae+C/29/DNUwAs9r6S318HtX3/QnMDNmgvrPwI581PXBhGxhcnmWkYN6cxocxfD0isBuMh1kk8MnMD98PZLOdmZBX2n4NAd6WyhiNicAj+Zll5pyigD/121l7WrOni+JY/PfOtCc3/kMbNYiohIGijwk8mRDSvfAEDevD5u/oypj/+Nn7p46O9Os8/B35vFz0VEppkCP9nmF5uLrYBXv/Qo737D80SjFpWffxE9PRacCsORP6e5kSJiRwr8VFjxGnOSFvjqJ//BYucp/vbkQr72f8Xm/qP3m5k7IiLTSIGfCjnzYbkpoby08DRf/aRZA/d/b15viqtF+2H/r6G/L52tFBGbsU3gt7S0JFzElXKLyyBvJQDveeN+Nm08Sld3Ntd94YVEo0BXCxz7y/S0RUQEGwX+tLMcsOrNYDmwLLj5M4+SO7ePxgeWcetvzjH7HLkHuo+lt50iYhu2CfyioiJCodC46+cnxbyi+Nz889d08NnrngRMcbUDh3PN6lgHfm3KL4iIpJhtAj9tlr0ccpcC8PH3NLPx0jDHT+Sw9XMvMjnf8Ty0PpTeNoqILSjwU82RPTC0Y5GdHeWWzz/C3Dl9/OG+5fzwVwNDOy1+M11TRCSFFPjTIf8cU0YZuHDtSb7wkScA+NiXL+G5Q/PMhVj7f6WKmiKSUgr86bL8apjjBOBj727miheFaT+Zw7X/GxvaeQ6OataOiKSOAn+6ZM2B1W82f82CWz73cHzWzo6fn2v2OXw3dB1OYyNFZDZT4E+nApeZnw+sP6+D7f/xTwA+8ZWL2bc/D6J9sP8XZvaOiEiSKfCnW9Gr4xU1/+OdIa4qaeVkZzbv+fSL6e21TA//yJ/S20YRmZUU+NMtaw6sfgtYFg4H/OjzDzM/v4f7H1nM9u+vM/scvd+M6YuIJJECPx3y18DSlwGwdnUn3/nU3wH4bN16Hnh0kbkQ6/lfQF93OlspIrOMbQJ/2mvpnM2yV5krcYF3vWE/73jdfvr6HLxrWwknOrLgdAQO/EZX4YpI0tgm8DOOIwtW/ytYWVgWfOdTj3Luyk5C+/P5jxsvNftE/gFtwfS2U0RmDdsEflpq6ZzNvOVQdDUAzgW9/PiLQRyOKLfcvoZdd5hKmxz8I3QfSWMjRWS2sE3gZ6wlV0DBWgCuKgnzqWv2AlD1+ReeuQr3OZ+WRRSRKVPgp5vlgHPeCtl5gFn8fOOlYSIn5rClupTTPZbp4R+6I80NFZGZToGfCXLmm/F8ICcnys9uDOKcf5oHHi1k29dfYPZpDZgxfRGRSVLgZ4oF68zwDmaq5i2ffwSAr/64mF/fbWbzcODXWjBFRCZNgZ9JijbBvBUA/MurWvj4e5oBeN9/vZjQ/jzoOw3P7TR/iohMkAI/kziyYU2FuRoXuPGjj3PFi8yCKZuvL+PUaQd0H4UDt2t+vohMWPZUH+CRRx4hHA7HpzvGLm4677zzpvrQ9jR3Max8Izz/C3JyouysDXDZ5lfS9LiTT37lBXxz22MQeQzyzoElG9PdWhGZQSYc+O3t7dTV1bFz50727dvH2rVrKSwsxOl0AhCJRAiHw+zbt4+ysjI8Hg+bN29mwYIFyW777LXohdD5PLTu4Zyibn78xSBv+MjlfOs2Fy+9rI13vO6AmbUzb6VZXEVEZBwmFPg33XQTu3btoqqqioaGBtauXTvm/vv27cPv93P11VfzwQ9+kGuvvXZKjbWVFa+BroPQeYDXX3WEbdfsZfv313PN/76Ii9ae4LIL2+G5XXB+pZnlIyJyFuMawz9+/DibN2+mtLSUPXv2cO2115417AHWrl3L1q1bCQQCrF27luuuu4729vYpN9oWHNmwZnN8fv7nP/wEr33ZYbq6s3nLx17CsbY50HMCnr1N9fNFZFzGFfi7du1ix44dXH311ZN+ok2bNnHjjTeya9euST/GVGRc8bTxmLMQznkbWBZZWfDTG4MUn9PBswfz2FJdaurndx7QSVwRGZdxBf7WrVtZuHDhlJ9s4cKFGtaZqPnFsOyVACxa0MOvvvYQ+fN6ueuhpVR/beCirLZH4ZjWwxWRsdlmWmZGFk8br2UvhwXrAbhk3Qlu/cLDAHzt/4r5v9+uNvu0+KF9b7paKCIzwJQC/6abbtKY/HSwLDO0k7sUgLe6D/HprSbct37uRQT+sXBg0ZSfm3n6IiIjmFLgX3/99WzatIlHHnkk4fYdO3awbdu2qTy0DJU1F859B2TPA+Cz1z3BG65qoftUFm/6j42msmbfKXjmp9BzMs2NFZFMNOUhnc2bN3Pbbbdx1113xW/bunUrbW1tfOpTn5rqw8tgcwthjQcsR/wk7qXr2mk5lssb/30j7Sez4XQbPPszlV8QkWGmFPhbtmzBsixuvPFGdu/enRD65eXl1NXVTbmBMkSBC1a8FoAFBb389psPUrSkm78/tYC31wyaufP8zyHan+bGikgmmVLgNzY2xqc53njjjTQ1NcVDPxQKUVpaOvUWynCLN0CheW/XrOji9q8/xLzcXv5w33I+dtPFZp/2J81qWZquKSIDphT4N9xwA42NjfGfr7/++nhP37Isdu/ePeUGyggsC1a+HgrOA2DDJRH+74sPY1lRvnWbi2/+dOCiuNaH4NgD6WuniGSUKQV+dXU1FRUVCTN1Yj39wQcCSQFHFqzZAnOXAGbmjvc/HwfgP2+65EwN/ZbdWjhFRIAknLTdtGnTsMJo119/PZWVlXz5y1+e6sPLWLLnwdp3QU4BAJ98XzNb3/YM/f0Wb68p5d5g4cB0zV/AieY0N1ZE0s2KRjNrkDcUCuHz+XC5XIRCISorK+OVOEfi8/kIh8MUFhYCUFFRkXD/6tWrOXDgAKtWrWL//v2pbHr6dB6E0A+hv4feXou3faKM2/+0Auf80/z5h/dz6boT4MiBte9VdU2RWWCyuTauwG9vb09aeeOzPVZpaSlNTU2AKbXs8XhGHR6qr68nEolQXV1NKBSivLyc5ubEnmzsjVm5ciV79458JWp+fv4kX00Gad9rpmNGo3R2ZfHqD17O/Y8sZuXSLv5y632cu7ILsnLB9QGYtzzdrRWRKZhs4I9rSGfnzp0JUy4n66677hqzeFowGEz42el0EggEiEQiI+5fU1NDdXU1YBZeiR0oRnLw4EEKCgpG3GaFBevNiVwgb14ft3/jIS4ubufg0Xm8+oNXmOqafd3wzI/hVDjNjRWRdBh38bTm5ma2bdvGM888M+EneeaZZ/jgBz9IMBgcs3haIBCID83EFBYWEggEhu0bOzj4/X78fj81NTWEwzYPssUbTN0doHBhD3/87gOcU9TJ3mcLeMNHNnKyM8tchbvvx6a0sojYyrgXQNm6dSv79u2jsrISy7IoLy/H7XbjcrmGDdG0t7cTCoXw+/3s3r0by7K4+eabz1pDf6SevNPpHPH2WM8/VvK4rKyM0tLSYUM6MWMN6cwqy18FfV3QuofVy7u547sPcOX7r+Shxxbxpn/fyO++9SB5tEHoR+B6f/yEr4jMfhOapbN27Vp2797Nrl27WLhwIdXV1ZSUlOBwOMjKyopva9eupaamBqfTSUNDA3fccce4FkwZ6eRsJBIZ8XaXy4XT6Yxf+OV0OgmHw8OGhWIsyyI/P3/EbVaJzdF3XgLARa6T/OHbDzA/v4c/BZbw1o9vMIuhnzoG+34EvR1pbrCITJdJTctcuHAhW7duZffu3Tz99NP09/fT19dHOBymr6+P1tZW7rjjDq699toJ1dEvKysbNiwTDodHXLTE5XIN6/mPNtZvO5YFq/8V5q8D4CWXRvj9tx4kL7eXO/6yjC3VpfT0WKayZkihL2IXSa2HP9VFUkpKShJCe/CQDZhx+1g9e5fLRUlJSfznUCgUv00wF2aduxny1wBwZUmY27/xEHPn9PHru1fw7k+VmLo73Udg363Q25nmBotIqk068Nvb2yd1AvdsGhoaqK2txefzsX37dhoaGuL3bd++HZ/Pl7BvXV0d9fX1eL1eXd07lCMHznsnzDNX3W7aeIxffm0POdn97Nq9imv+9zL6+4Guwwp9ERuY0IVX7e3t1NTUUF9fn3C7y+XC4/Fwww03JG2+frLY4sKrs+ntNOP1XYcB+OWdRXiuL6Ovz8H73vwc3//fR8jKAnKXmYuzdCJXJKOl9MIrgOPHj1NaWorT6WTLli3xK2Gbm5sJBAIEg0Esy6K+vp5rrrlm0i8k2RT4A3o7zHh99xEAdt2xknduK6Gvz8E7XrefW7/wMNnZUZi72IT+nKmvYSwiqTHZXBv3tMwbbriBqqoqrr/++hHvP378OF/60pfYunUrTU1NfOc73xl3I2QaZOeD630QugW6j7L5NQfJzoqypaaUn/1hNT29Dn66vYkcWk2ZhrXvNQuuiMisMe4x/Obm5lHDHswJW6/XSzgc5qGHHuIrX/lKUhooSZSdD2vfl1Bh8xdf2cOcnD58jSvxXF9mpmyejpjQ1/q4IrPKuAN/rAJmQ/cLBAL87Gc/S8lJXZminALT05+7GIA3vfIwv/76mdk7b/34Brq6HeZK3NAPTWE2EZkVxh34Q0senM2uXbvwer0TbpBMg5z5poha7jIAXvuyo/zuWw8yL7eX39+7nNd+6HKOn8geONl7i0ori8wS4w58y7Im9MAul4vW1tYJNyhVWlpaEub0215OgSmtkLcSMFM27/juAywo6OHPTUt4xTUvo+XYXLMY+jM/gbZH09teEZmycQd+IBDgxAkV3JpVsvPMmH7+uQBcVRLmzz+4n+WLu/nbkwu58v1XEtqfZxZDf/4XcPR+rZErMoONO/CbmppwOp2sW7eO6667jl/+8pcJSxuOZKLDQKlUVFREKBSKX5krA7Lmwtp3w/xiAF50QTv3/+g+XKs7aH4+n5e970oe3TtwbcWhRjikhdFFZqpxB35JSQk333wzV199NY2NjbztbW9j0aJFLF68mC1btvD9739/2EnaiQ4DSZo4cuDcd8CCCwEoPqeT+265jxeuP07LsVxe/m8v464Hzcwejj0IzzVAf08aGywikzHuwHe73WzdupW6ujqefvpp2tra2LlzZ7xW/tatWykuLiYrK4vXvOY1fPnLXx6xjr1kKEe2qb1TWArAiqWnuOf793NVSSvHT+Twmg9dzg9/NbA84vHHzQwe1dQXmVGSuqat3++nsbGRhx9+GL/fj2VZ9PX1JevhJ0VX2k5QNApH/gyH7wag+5SDD/z3Zdz2x9UAfHrrXj73oSdwOICcBXDeO2DeijQ2WMR+Un6l7Xi43W7cbnf857KysmQ+vEwHy4LlrzBTNw/8hty5/fxke5Diczr54o71fHHHepqfz+OHn3uEXNqh+Qdwzlth4UXpbrmInEVSyyMPpcCfwQpL4Ny3gyMHhwO+8JEn+MFnHyY7u5/b/riaTZUv5Wh4jhnLf3YnHLlPJ3NFMlxKA//mm29O5cNLqi24wFyVm50HwAfe8jx3fOcBFs7v4S+PFLLx3VedmcHT4ofnf27m7YtIRkpp4MsskLcazt8avyr36o3H+Out97J2VQf7DuRzxXuuxNc4MIYfeQyavwenMueCOxE5Q4EvZzdnERRfE18y8SLXSfb89M+4Lz9KZ3c2nk9u4NPfvJC+Pkz55ad3QLsNFowXmWEU+DI+WXPNjJwllwOw2NnDH779AJ9839MAfOl763nzRzcSac+Gvm545qdmpo/G9UUyhgJfxs9ywMrXwqo3guUgOzvKTR9/nJ9sbyJ3bh+/v3c5L3nXy3nsqflm/8P3mODX0okiGcE2ga/iaUm0uMyUYxg4mfvO1x/g/h/dx5oVnTz1XAEvefdV3PLrgYu0TjwFT90MHc+mscEiAjYKfEmyAhecXxmvtlly0XECP/0zr77iCF3d2Xzgv1/M+//rMjo6s6Cn3ay0dfgeU4hNRNLCNoGv4mkpMMcJrn+DwhcDsLTwNH/4zgN84SP/xOGI8qPb1/CSd13F480FZiz/8N2w78cqySCSJrYJfEkRRzasevPAuH4WDgd8eutT3Fn/F4qWdPN4aAEb3vlybv2NKc3AyX1miKf9qfS2W8SGFPgydZZlxvWLP2Dq6wCv3NDKI7vuYdNGM3XzfZ8p4Z03lJhZPL0dZlGVA7/ThVoi00iBL8mTtxrWVcXn6y9ffIo7vvtXPvehJ8jK6udnf1jNCz2v5E97zHq6tO6Bp+ugU0XtRKaDAl+SKzsfznsnrHg1WA6ysuC/qvZy/y33cf6akzzfksfVW19Kzdcu4tRph7kqt/kHZny/P72VVUVmOwW+JJ9lwdKXQvG/mRO7wMYXRnh45z1c+9ZniUYtam9Zx+Xvjp3Q7TczeELfh+6j6W27yCymwJfUyVsN6z4IC18AQEFeHzv+52/88msPsdh5ikeeXEjpO17BV35UbMoydB40J3SP/Fm9fZEUUOBLamXlwhoPrH6TWUoReMvVLfzd9yde+7LDdJ/K4pNfvZgr338l/wwVQLQPWu6C5npzABCRpFHgS+pZllk6cd11ptePWULx999+kO/9zyMsKOjhgUcLefGWV3Dj98+nt9eCrsPQvMMsnK71c0WSQoEv02duoRnXL7oaLAeWBde89Tn+8Yu7ef1Vhzl1Oott33gBl7/nKv7+1HxzsdbR+80wz8l96W69yIynwJfpZTlg2cuh+FrIXQrA6uXd/PabD/KjLwRxzj9N0+NOSt/+Cv772xfQ1T0wkyf0I3j+F9BzMs0vQGTmsk3gq3hahslbaWrxLLkCLAvLgve+aT+P//Ju/uVVh+jpdfD5+gu4tOJV7P6LOTDQ9ijs/SYce0g1eUQmwTaBLxnIkQMrX2MWVxno7a9Yeopffm0Pvq/sYeXSLpqfz+c1113BO2pKOHR0LvSdgoO/N4usdB5I8wsQmVlsE/gqnpbB8lbD+VWw/JXxsf23uQ/xxK/v4qPvasbhiHLbH1dz4Vuu5tu3nWemcHYdMssp7v+NKdUgImdlm8CXDOfINoF/flW85PL8/D7+X/U/2PPTP1N2cRvtJ3P4yPYXcsV7ruKhvzvNSd1wEzz5DTj6F+jvTetLEMl0CnzJLPOWmxO6K14NWXMAU2v/gR/fy7e2PcqCgh72/GMRG9/9cj7wX5fRcmxgmOfQbnjqO3D8CS2rKDIKBb5kHsthSjOs/0j8Kt2sLPjw25/hiV/dxXve+DwAt9y+hvVvvpqbbinmdI8Fp8Lw7G2w71Yzj19EEijwJXPlLIBzN5vlFOcWAuak7q1ffJi//vheNlzcxomOHKq/djGXvPVV/O7Py8zvndwHT98Mz/8KTkfS1nyRTKPAl8w3/3xY9yEzxu/IBuDyF7bxwP/dyw8++zDLF3fz1HMFvPHfL+f1H97Ik8/km2GdtkfgyW/CwTu0kLoICnyZKWInddd9CBZcYG5ywAfe8jx7b7+L69//NDnZ/fzhvuVc8rZX8dmb15vfi/bBsb/Ck183Rdm04IrYWMYFfigUora2Fp/PR21tLZFIZFy/V1VVNe59ZQabWwjnvQNc74VcM4SzoKCX2o89zmM/v5s3vryF3l4HywqHBHvfKVOUbe83zIVbmtEjNmRFo5k1paG0tJSmpiYAIpEIHo+HxsbGMX/H7/fj8XjYt28fTqcz4b7Vq1dz4MABVq1axf79WllpVon2QzgIh+9KGLK5J7CYl10WJjt7jH/aOQtg2VWw6MXxYSKRmWKyuZZR/9KDwWDCz06nk0AgQCQSGRbkMbFefWFh4ZiPHY1G6egY+QKd/Pz8CbdVMoDlMGvpOi8xwzXHHoRoH68oaz377/a0mzV1j9yr4BfbyKghnUAgMCy4CwsLCQQCo/7Orl27cLvdZ33sgwcPUlBQMOImM1xWrpm3v/7DJvwnIhb8T2qoR2a/jAr8kcbgnU7nqGPzfr+fzZs3p7ZRMnPMLYQ1FQMLqZ8/sd/taTc1ep78Bhz9q07uyqyUUd9hRxq2GW04J3YQGG2oZ6iVK1eyd+/eyTdOZo55K8zc/ZP7oMU/sSJrPe1w6A44+mco3ABLNpqF2UVmgYwK/LKyMurq6hJuC4fDI5Y09vv9hMNh6uvrATO7p76+HrfbTUlJybD9LcvSWL3dFKw1ZRran4CWO+HUsfH/bm/XwHmBv5jx/aUvhTmLUtdWkWmQcbN0iouLaW5uBkwvftOmTfFZO8FgEKfTOeIBwLIsmpubh92nWToCmBk9x/8Bh++ZWPDHWBYsvBiWvsx8gxBJo8nmWkaN4QM0NDTE5+Fv376dhoaG+H3bt2/H5/Ml7B+JRKitrQXA6/Wq/LGMzHKA81JY/yFY8zaYu2Rivx+NQuQxeKoOQrfA8X9qERaZcTKuh59s6uHLiKL9cPxxOHIPdB+d3GPMccLil5ghn+x5SW2eyFhmxTx8kWljOcwUzoUvMMF/9N6JV9g8HTFlmQ/fDYsuM+E/sHKXSCZS4Iu9xYP/YjjxNBy9Dzqendhj9PdA6x6zzS+GwjKYvx4cWalps8gkKfBFwJyUXbDObJ374ch9cOLJiS+mcqLZbDnzzVBPYYkZ+hHJAAp8kaHyVsN5b4fuY3Dsfmh71FTdnIieE2Za59F7oeB8UwJi/jrzjUIkTRT4IqPJXQKr/wWWXw2tD5n1cydaVz8ahRNPmS1ngenxL7pMvX5JCwW+yNnkzIeiTbDs5WZqZusDk1tCsacdDv/JbAVrTfAvuCi+dq9Iqtkm8FtaWuIXZWmuvkyKIwcKX2yCuuNZaH3QXMU7mZnNJ/eZLet35oTxossgb405lyCSIrYJfJGksSwoOM9spyNmuKftYVOOYaL6TkP4YbPNWWSCf9ELVcZBUsI2gV9UVKSevSTfHKcpzbz8VWY+f7gJOp6b3GOdbjNz+g/fDfnnwMKB6aI5KuEtyWGbwBdJKUcOLHqR2bqPmOBv+xv0dU/u8TqeN9uhP0KBayD8LzK1/0UmSYEvkmy5y2Dl66DIbU7yhpvM3P7JiEbPzO0/8FsztdN5KSxYbw4yIhOgwBdJldhJ3sIXm1k9bY9A5FHoHXmpzbOK9pmTxO1PmMeev86Uhpi/DrLmJrXpMjsp8EWmw7zlMO81ptd/8mkT/u17J35BV0x/jzlncPxxsxZvgWsg/C9QITcZlQJfZDo5smDBBWbr7YTI3034dx2a/GP295qDR/tecyVv/nkm/BdcqBO+kkCBL5Iu2XlmCcUlG82QT+RROP4YnD4++ceM9sPJkNkO/s6UiZh/gRnzn7tU8/xtToEvkgnmLYd55WbIp/M5c7L3+D8mXsphsGj0zGyfFr+Z279gvankmX+uGQoSW9EnLpJJLMuEcf65sOK10LHPDPu0PwF9p6b22Kfb4NiDZsuaY4q6LVhvTvpqoXZbUOCLZCpHFsw/32z9PaYAW+Qx82d/z9Qeu+/0mZO+lgXzVpnnKSiGvFWq6jlLKfBFZgJHjjkRu/AFA+HfDO2PmxO1k724KyYaNdcJdO43hd2ycs2sn/nF5lvAnIVJeQmSfrYJfBVPk1nDkQMLLzRbf58Z9jn+uBn2mcqYf0xf95neP5gF32Phn3+uqnvOYLYJfJFZafCwT/SNpo5P+z/h+D9NOeZkOHXMbMceBCsL8teYqZ8Fa81QkJZynDFsE/gqniaznuU4U8VzxWuh+zCc2AvtT0LngeQ8R7TvTGnnw3ebbxv555rwzz8P5q3Q+H8Gs03gi9iKZcG8IrMtezn0nBxYeWsvnGw2J22Tob/HLP5+4mnzc1buwCyj88xBIHe55v5nEAW+iB3kFJyp69PfCx3PmBO+J/aamv7J0tdtvlG0P2l+zp4HeeeYg0DeGpi3UkNAaaTAF7EbR/agcf/XwalW0+s/8bQ5EEx1yudgvV1nyj7EnnveqoFvAWvMwUCF36aNAl/EzizLLNaeu8SUeOjvNVf6nmg2Rd4ms3bvWPp7zfKQHc8Oev7lZ74B5K8xawhLSijwReSMWOXNAhdQDj0nTF2eE0+bPydb2nk00Sh0tZiNB81tOQtMDaDYNm+Fav8niQJfREaXM//MSl7RKJw6emaWTsczU7/oayQ97YnXAViOgRPQq84cBOYU6mTwJCjwRWR8LMus5pW7zAz/RPvN1M/YAaDz2eTN/hks2g+dB83Wusfclj0P5sW+AayEvJWqBzQOCnwRmRzLYYZb5q2ApS81V/12HzpzAOjan5oDAJiTwSeeMltMzgIT/vNW6CAwCgW+iCSHI+vMkMuyq0zPvKvFnATueNZcBZzscwCD9bSbrf2JM7fpIJBAgS8iqWE5TMDmrYQll5tzAKfDJvg7njUHglPh1LZhzINAkZkhNK8Icpy2OCdgm8BX8TSRNLMsmLvYbIUvNrf1nDAHgM7nTLXOrpbJr/M7XiMdBLLmmvCPHQByl8PcZbOuUJxtAl9EMlDOfHBebDYw8/S7Dg2E/wHzZzKvBB5N36mBbx7PnbnNssxsoMEHgdwi8w1hhn4bsE3gq3iayAzgyIb8c8wW03Ny4AAwULO/62DqTgYPFo2aq5BPtZ6ZIgqmXlDuUvMNIHepmbU0dylkF2T8gcA2gS8iM1ROwZn6/zAwHfTowAHgoJkZ1HU49UNBMX3dZ9YKHix7ngn+2AEgNoU1g04SK/BFZGaxHAOLvi+HwlJzW38vdB8x4Z+OgwCYqaK9Q4aFALLzzhwE5i41ZSzmLoHs+dP+jUCBLyIznyP7zIygTDoIgFmF7OQzZhssaw7MWWzCf27sz4G/p6iUhAJfRGanEQ8CfXDqiJkN1H0Yugf+7O2a/vb1nTYnqLsODb9vzsIzB4A5i5P2rUCBLyL24cg6c3VwTDRqpmkOPgB0tZhrBqLR9LTz9HGznWhOvD1rjpk51De5A5QCX0TszbJMj3rOQliw/sztfafNt4HYAaD7sNn6TqWvrX2nTVsmuWaBAl9EZCRZc86UioiJRqH3xMC5gaMDB4SjpopoOg8E45RxgR8KhfD5fLhcLkKhEJWVlTidzhH3DQaD+P1+APbs2YPX641fTSsiknSWZS68yllgVgyLiQ8LHTHhH/vz1NHpuWZgnDIu8D0eD01NTQBEIhE8Hg+NjY3D9otEIvj9fqqrqwHw+XyUl5fT3Nw8bF8RkZQaPCzEujO3R6PQc3zQAaAVTh0zW2/ntDczowI/GAwm/Ox0OgkEAkQikWG9/EAgQE1NTTzw3W43oVCIUCg0Yi8/Go3S0TFypb78/My5MEJEZhHLgjlOs7E+8b7ezsQDQOzvp8Pm4rIUyKjADwQCFBYWJtxWWFhIIBDA7XYn3O52u+PfBGK/C4w6pHPw4EEKCgpGvC+arjPxImJf2XlmG1xGAszU0Z7I8ANBEr4VZFTgRyKRYbc5nc4RbwcoKSmJ/93r9VJXV5eilomITBNH1pmqolyQeF9vF5xug+zvAicm/NAZFfgjnZwdaThnqPr6ejweD5WVlaPus3LlSvbu3TvFFoqIpFH2PLNZk7sS15Hk5kxJWVkZ4XDiggjhcHjMmTd+v5/CwsIxwx7Asizy8/NH3ERE7CCjAr+kpCRh+CYSieByueKBHwwGE0ocx07yVlRUAKanP9rwj4iI3WXUkA5AQ0MDtbW1uFwu9uzZQ0NDQ/y+7du3s2HDBqqrqwmFQpSWlib8rtPpPGtPX0TErqzoLJ+isnr1ag4cOMCqVavYv39/upsjIjJlk821jBrSERGR1FHgi4jYhAJfRMQmFPgiIjahwBcRsQkFvoiITdgm8FtaWhIu4hIRsRvbBL6IiN1l3JW2qVJUVJRQlkFExG7UwxcRsQkFvoiITSjwRURsQoEvImITCnwREZtQ4IuI2IQCX0TEJhT4IiI2ocAXEbEJ2wS+aumIiN3ZJvBFROxOtXRERGxCPXwREZtQ4IuI2IQCX0TEJhT4IiI2ocAXEbEJBb6IiE0o8EVEbEKBLyJiEwp8ERGbUOCLiNiEbQJfxdNExO5sE/giInan4mkiIjahHr6IiE0o8EVEbEKBLyJiE7M68Ds6Ojhw4AAA0Wg0za2ZPh0dHViWhWVZdHR0pLs500avW6/bDqaSa7M68EVE5AwFvoiITSjwRURsQoE/golekTuZK3in4zkmSq87NftP13NMlF53avafrueYjIy78CoUCuHz+XC5XIRCISorK3E6nVPeV0TE7jIu8D0eD01NTQBEIhE8Hg+NjY1T3jcajY56Jj8/Pz8JLRcRyWwZFfjBYDDhZ6fTSSAQIBKJDOu5T2RfgEOHDlFQUDDi865atSrh55aWFgBWr149rnZPdP9UP8fgqVrr16/Hsqy0t2k6nkOvW6873W2ajucY/LqPHj067jbFfjlj1NXVRd1ud8JtLpcr2tjYOOl9HQ5HFNCmTZu2Wbc5HI4JZWxG9fAjkciw25xO54i3j3ff3Nxcuru7ycrKYunSpSM+73h7BiIimeDIkSP09fWRm5s7od/LqMAfaShmtCGa8e5rpyvwRETGklHTMsvKygiHwwm3hcPhEacrTWRfERHJsMAvKSlJGJKJRCIJ81ODwWC8pv3Z9hURkURWNJpZVcWCwSB+vx+Xy8WePXuoqqqKh7jH42HDhg1UV1efdV+78/l8VFRUpLsZItOiqqoKr9er63DOIuMCP5nseGGWz+cjHA5TVVVlqwqhcKYDALBnzx68Xq8tOgB+v59IJEI4HKaxsZFt27ZRUlKS7mZNG7/fj8fjYd++fbP+/zeYg1t9fT1gRjp27Ngx7s87o07aJttELsyaLWK9+qqqqjS3ZHpFIhH8fn/825/P56O8vJzm5uY0tyz1ysvLaWpqin/2Ho/HFq8bzszWKywsTG9DplFpaSltbW3AyJNXxpJRY/jJNNaFWTL7BAIBampq4j+73W5CoZAt1jFubm5O6OHZoZcbs2vXLtxud7qbMe2cTuekPudZG/iBQGDYUb+wsJBAIJCmFkkqud3u+Lc5IP4522FIZ/BrbG5upqGhIY2tmT5+v5/NmzenuxnTLhKJ4PP58Pv91NTUTKhTM2uHdCZyEZfMDoN7uV6vl7q6ujS2ZnqFQiHq6uoIBoO2+Dcee412+jYTU1FRET/IFxYWJgxdn82s7eFP5CIumV3q6+vxeDxUVlamuynTxuVy4fV6KS8vZ9OmTbM+9P1+P6FQiPr6eurr6+N/HzqUOxsN/kbncrkmdpCfUCGGGaSpqSlaUlKScJvT6Yw2NzenqUXTaxZ/tGNqbGyMNjQ0pLsZ06a5uTlaXV2d8DMwYv2p2Qywxf/tpqamqNPpTLgNiLa1tY3r92dtD9+uF2b5/X5qa2sBqKmpiU9TtINY7y42W6W+vn7W93RDoVDCZxwKhXA6nZSVlaWxVdMnEonE/717vd5Zf5Le5XKxbdu2+M9+vx+32z3ukYtZPQ9fF2bZRygUori4OOE2p9MZn742m8XmZAO2nIdvN36/n2AwiNPppLm5mW3btinwRUQk0awd0hERkUQKfBERm1Dgi4jYhAJfRMQmFPgiIjahwBcRsQkFvoiITSjwRURsQoEvImITCnwREZtQ4IuI2IQCX0TEJhT4IklSW1vLokWL8Hg8CbeJZApVyxRJAr/fT0NDA+Xl5fGa7Dt37uTOO+/UKmuSMRT4IkkwdPnMmpoatmzZorr0klE0pCOSBEPDvry8XGEvGUc9fJEkqq2tpaSkBLfbne6miAyjHr5IktTU1CjsJaMp8EWSoLa2lvLy8oSwH7zWrEgmyE53A0RmOp/PR2NjIwDBYBCXy0VjY6Nm50jGUeCLTFEoFIoHfnl5OYFAgMrKSrxeb5pbJpJIJ21FRGxCY/giIjahwBcRsQkFvoiITSjwRURsQoEvImITCnwREZtQ4IuI2IQCX0TEJhT4IiI2ocAXEbEJBb6IiE0o8EVEbOL/A8K84yPTPv/0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (4,4))\n",
    "\n",
    "## normalize such that D(z = 0) = D(a = 1) = 1\n",
    "plt.plot(1/np.exp(res.t)-1, normalized_D, color = 'orange', label = 'Normalized D(z), Integral',\n",
    "         lw = 5, alpha = 0.5)\n",
    "\n",
    "D_normalized = [1, 0.76872625, 0.60653086, 0.41693915] ## from other notebook\n",
    "redshifts = [0,0.5,1,2]\n",
    "\n",
    "spl = UnivariateSpline(redshifts, D_normalized)\n",
    "test_z = np.linspace(0,2,200)\n",
    "\n",
    "plt.plot(test_z, spl(test_z), color = 'mediumblue', label = 'Normalized D(z) from Sim.')\n",
    "\n",
    "\n",
    "# aesthetics\n",
    "ax = plt.gca()\n",
    "ax.minorticks_on()\n",
    "ax.tick_params(which =  'major', direction = 'in', length = 6, width = 1.5)\n",
    "ax.tick_params(which =  'minor', direction = 'in', length = 3, width = 1.25)\n",
    "plt.xlim(0,5); plt.ylim(0,1.2); plt.xlabel('$z$', fontsize = 14); plt.ylabel('$D(z)$', fontsize = 14)\n",
    "plt.legend(loc = 'upper right')\n",
    "plt.savefig('CompareGrowthFunctions.pdf', bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fdd8309",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c243e88",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}