class_sz_electron_galaxy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a3461285-3f18-4ba7-9619-4ead5fa207c5",
   "metadata": {},
   "source": [
    "# Intialize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "dd77bdd4-5984-4527-9a21-82574543d1e8",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:absl:At this time, the v2.11+ optimizer `tf.keras.optimizers.Adam` runs slowly on M1/M2 Macs, please use the legacy Keras optimizer instead, located at `tf.keras.optimizers.legacy.Adam`.\n",
      "2024-02-20 18:22:59.720070: I metal_plugin/src/device/metal_device.cc:1154] Metal device set to: Apple M1 Pro\n",
      "2024-02-20 18:22:59.720093: I metal_plugin/src/device/metal_device.cc:296] systemMemory: 16.00 GB\n",
      "2024-02-20 18:22:59.720102: I metal_plugin/src/device/metal_device.cc:313] maxCacheSize: 5.33 GB\n",
      "2024-02-20 18:22:59.720137: I tensorflow/core/common_runtime/pluggable_device/pluggable_device_factory.cc:303] Could not identify NUMA node of platform GPU ID 0, defaulting to 0. Your kernel may not have been built with NUMA support.\n",
      "2024-02-20 18:22:59.720159: I tensorflow/core/common_runtime/pluggable_device/pluggable_device_factory.cc:269] Created TensorFlow device (/job:localhost/replica:0/task:0/device:GPU:0 with 0 MB memory) -> physical PluggableDevice (device: 0, name: METAL, pci bus id: <undefined>)\n",
      "WARNING:absl:At this time, the v2.11+ optimizer `tf.keras.optimizers.Adam` runs slowly on M1/M2 Macs, please use the legacy Keras optimizer instead, located at `tf.keras.optimizers.legacy.Adam`.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:Please fix your imports. Module tensorflow.python.training.tracking.data_structures has been moved to tensorflow.python.trackable.data_structures. The old module will be deleted in version 2.11.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:Please fix your imports. Module tensorflow.python.training.tracking.data_structures has been moved to tensorflow.python.trackable.data_structures. The old module will be deleted in version 2.11.\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from classy_sz import Class\n",
    "import os\n",
    "import time\n",
    "\n",
    "\n",
    "# the parameters needed for cosmology:\n",
    "# use the last column of Planck 2018 (https://arxiv.org/pdf/1807.06209.pdf)\n",
    "# TT,TE,EE+lowE+lensing+BAO\n",
    "cosmo_params = {\n",
    "'omega_b': 0.02242,\n",
    "'omega_cdm':  0.11933,\n",
    "'H0': 67.66,\n",
    "'tau_reio': 0.0561,\n",
    "'ln10^{10}A_s': 3.047,\n",
    "'n_s': 0.9665,\n",
    " \n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "common_params = {\n",
    "\n",
    "'z_min' : 0.005,\n",
    "'z_max' : 3.0,\n",
    "'M_min' : 1.0e10, \n",
    "'M_max' : 3.5e15,\n",
    "    \n",
    "\n",
    "# 'delta for galaxies' : \"200c\",\n",
    "# 'delta for matter density' : \"200c\",\n",
    "# 'delta for electron density':\"200c\",    \n",
    "\n",
    "'mass function' : 'T08M200c',\n",
    "'concentration parameter' : 'B13',\n",
    "\n",
    "'redshift_epsabs': 1.0e-40,\n",
    "'redshift_epsrel': 0.0005,\n",
    "'mass_epsabs': 1.0e-40,\n",
    "'mass_epsrel': 0.0005,\n",
    "\n",
    "\n",
    "'ell_max': 1600.0,\n",
    "'ell_min': 2.0,\n",
    "'dell': 10,\n",
    "\n",
    "'non_linear' : 'hmcode',\n",
    "\n",
    "'hm_consistency' : 1,\n",
    "# 'x_outSZ': 4.,\n",
    "# 'truncate_wrt_rvir':0,\n",
    "}\n",
    "\n",
    "\n",
    "import os \n",
    "path_to_class_sz = os.getcwd() + '/../../class_sz/'\n",
    "\n",
    "# a simple conversion from cl's to dl's\n",
    "def l_to_dl(lp):\n",
    "    return lp*(lp+1.)/2./np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "67fe81d1-95dd-4d5e-8b61-0f317c03ae42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/Users/boris/.local/lib/python3.9/site-packages/classy_sz-2.9.4-py3.9-macosx-11.0-arm64.egg/classy_sz.cpython-39-darwin.so'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import classy_sz\n",
    "classy_sz.__file__\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ebb65920-5a2f-492a-a5e5-e04f03d4407a",
   "metadata": {},
   "outputs": [],
   "source": [
    "ksz_params = {\n",
    "'output': 'tau_gal_1h,tau_gal_2h',\n",
    "\n",
    "# 'pressure profile':'B16',\n",
    "# 'gas profile mode' :'agn',\n",
    "    \n",
    "\"ell_min\" : 2,\n",
    "\"ell_max\" : 20000,\n",
    "'dell': 0,\n",
    "'dlogell': 0.2,\n",
    "    \n",
    "'M_min' : 1.0e10, \n",
    "'M_max' : 5e15,\n",
    "\n",
    "'gas profile':'B16',\n",
    "'gas profile mode' : 'custom', # important to read values of parameters\n",
    "'use_xout_in_density_profile_from_enclosed_mass' : 1,\n",
    "'n_z_m_to_xout' : 30,\n",
    "'n_mass_m_to_xout' : 30,\n",
    "\n",
    "    \n",
    "\n",
    "'n_m_density_profile' :30, # default: 100, decrease for faster\n",
    "'n_z_density_profile' :30, # default: 100, decrease for faster\n",
    "\n",
    "\n",
    "    \n",
    "'k_min_samp_fftw' : 1e-3,\n",
    "'k_max_samp_fftw' : 1e3,\n",
    "'N_samp_fftw' : 1024,\n",
    "    \n",
    "    \n",
    "'hm_consistency' : 1,\n",
    "    \n",
    "    \n",
    "'use_fft_for_profiles_transform' : 1,\n",
    "    \n",
    "    \n",
    "'x_min_gas_density_fftw' : 1e-6,\n",
    "'x_max_gas_density_fftw' : 1e5,    \n",
    "    \n",
    "}\n",
    "\n",
    "HOD_blue = {\n",
    "'sigma_log10M_HOD': 0.68660116,\n",
    "'alpha_s_HOD':    1.3039425,\n",
    "'M1_prime_HOD': 10**12.701308, # Msun/h\n",
    "'M_min_HOD': 10**11.795964, # Msun/h\n",
    "'M0_HOD' :0,\n",
    "'x_out_truncated_nfw_profile_satellite_galaxies':  1.0868995,\n",
    "'f_cen_HOD' : 1., \n",
    "    \n",
    "'galaxy_sample': 'unwise',\n",
    "'galaxy_sample_id': \"blue\",\n",
    "'UNWISE_dndz_file': path_to_class_sz +\"/class_sz_auxiliary_files/normalised_dndz_cosmos.txt\"\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2b679b49-795a-49b5-9440-af6959351f95",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.77 s, sys: 85.6 ms, total: 1.86 s\n",
      "Wall time: 451 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "M = Class()\n",
    "\n",
    "M.set(common_params)\n",
    "M.set(cosmo_params)\n",
    "M.set(ksz_params)\n",
    "M.set(HOD_blue)\n",
    "M.set({\n",
    "'use_fft_for_profiles_transform' : 1,\n",
    "'ndim_redshifts':30, ## precision parameter -- how many z in redshift grid.\n",
    "})\n",
    "M.compute_class_szfast() # fast mode works\n",
    "# M.compute() # slow mode works"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9ed95e42-d23c-4715-a45a-bf61393cea5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "l = np.asarray(M.cl_sz()['ell'])\n",
    "cl_eg_1h = np.asarray(M.cl_eg()['1h'])/l_to_dl(l)\n",
    "cl_eg_2h = np.asarray(M.cl_eg()['2h'])/l_to_dl(l)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "59dcf226-9372-43eb-b8f1-e7fe08516356",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.00000000e+00, 2.44280552e+00, 2.98364940e+00, 3.64423760e+00,\n",
       "       4.45108186e+00, 5.43656366e+00, 6.64023385e+00, 8.11039993e+00,\n",
       "       9.90606485e+00, 1.20992949e+01, 1.47781122e+01, 1.80500270e+01,\n",
       "       2.20463528e+01, 2.69274761e+01, 3.28892935e+01, 4.01710738e+01,\n",
       "       4.90650604e+01, 5.99282001e+01, 7.31964689e+01, 8.94023690e+01,\n",
       "       1.09196300e+02, 1.33372662e+02, 1.62901737e+02, 1.98968631e+02,\n",
       "       2.43020835e+02, 2.96826318e+02, 3.62544484e+02, 4.42812832e+02,\n",
       "       5.40852815e+02, 6.60599120e+02, 8.06857587e+02, 9.85498082e+02,\n",
       "       1.20369008e+03, 1.47019038e+03, 1.79569458e+03, 2.19326632e+03,\n",
       "       2.67886153e+03, 3.27196886e+03, 3.99639179e+03, 4.88120396e+03,\n",
       "       5.96191597e+03, 7.28190061e+03, 8.89413350e+03, 1.08633192e+04,\n",
       "       1.32684880e+04, 1.62061679e+04, 1.97942581e+04])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "7e211b08-8a4e-4dc8-8eef-aafc5a66ae42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 1e-10)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGsCAYAAAAGzwdbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABPiUlEQVR4nO3de3xT9f0/8NfJtU3vF3qD3rjLvaKUCih+7URleJmbiEwQFSfCplYd4gW+uH0Fb4zNoWwooj8voJsypwynFeRWQIqVO3JpaYFeaEvTe5Mmn98fadKmLYWUJuckeT0fjzyanpwk70OEvPxcJSGEABEREZFCqOQugIiIiKgthhMiIiJSFIYTIiIiUhSGEyIiIlIUhhMiIiJSFIYTIiIiUhSGEyIiIlIUhhMiIiJSFIYTIiIiUhSGEyIiIlIUrwonW7ZswZQpU5CQkABJkrB+/XrZ308IgYULFyI+Ph6BgYHIzMzEsWPH3FoXERGRL/OqcFJXV4eRI0dixYoVinm/l19+GX/5y1+wcuVK7Nq1C0FBQZg0aRIaGxs9UiMREZGvkbx14z9JkvDZZ5/h9ttvdxxramrCs88+i48++ghVVVUYNmwYXnrpJUycONEt7yeEQEJCAp544gk8+eSTAACj0YjY2FisWbMGd99992W/LxERkb/xqpaTi5k3bx5ycnKwdu1a7Nu3D7/61a9w0003ua2bJT8/HyUlJcjMzHQcCwsLQ3p6OnJyctzynkRERL7OZ8JJYWEh3nnnHXzyySeYMGEC+vXrhyeffBLjx4/HO++845b3LCkpAQDExsY6HY+NjXU8RkRERK7xmXCyf/9+WCwWDBw4EMHBwY7bd999hxMnTgAAjhw5AkmSurw9/fTTMl8JERGRf9PIXUBPqa2thVqtRm5uLtRqtdNjwcHBAIC+ffvi8OHDXb5OVFTUJb9nXFwcAKC0tBTx8fGO46WlpRg1atQlvw4RERG18plwkpaWBovFgrKyMkyYMKHTc3Q6HQYPHtxj75mamoq4uDhkZ2c7wkh1dTV27dqFOXPm9Nj7EBER+ROvCie1tbU4fvy44/f8/Hzk5eUhMjISAwcOxPTp0zFjxgy89tprSEtLw7lz55CdnY0RI0Zg8uTJPfp+SUlJkCQJjz32GP74xz9iwIABSE1NxfPPP4+EhASnWT1ERER06bxqKvHmzZtx/fXXdzg+c+ZMrFmzBmazGX/84x/x3nvv4cyZM4iOjsbYsWOxePFiDB8+vMffD7BNJ160aBH+/ve/o6qqCuPHj8cbb7yBgQMHuvx+RERE1I1wsmXLFrzyyivIzc1FcXFxh7U/OrN582ZkZWXh4MGDSExMxHPPPYf77rvvMsomIiIiX+XybB1XV2nNz8/H5MmTcf311yMvLw+PPfYYHnzwQXz11VcuF0tERES+77K6dTpbNbW9+fPn48svv8SBAwccx+6++25UVVVh48aN3X1rIiIi8lFuHxCbk5PjtIIqAEyaNAmPPfbYBZ/T1NSEpqYmx+9WqxWVlZWIioqCJEnuKpWIiIh6kBACNTU1SEhIgEp16Z01bg8nJSUlna6gWl1djYaGBgQGBnZ4zpIlS7B48WJ3l0ZEREQeUFRUhD59+lzy+YqcSrxgwQJkZWU5fjcajUhKSkJRURFCQ0N77H0efPd77DxZiZfuHI7JIxJ67HUBAC/3AyxNwMM7gIiknn3tFjknKjD7vT3o2ysIn88b75b3ICIi6q7q6mokJiYiJCTEpee5PZzExcWhtLTU6VhpaSlCQ0M7bTUBAL1eD71e3+F4aGhoj4YTbWAQVPpGGIJDevR1AQDBBqDRBBi0QE+/dot+vSWo9AZUNKl7vn4iIqIe4uqQDLfvrZORkYHs7GynY19//TUyMjLc/dby0rYEL3OD294iNjQAAFDT1Iy6pma3vQ8REZEnuRxOamtrkZeXh7y8PACtq6YWFhYCsHXJzJgxw3H+ww8/jJMnT+L3v/89jhw5gjfeeAMff/wxHn/88Z65AqWyh5PmRre9RUiAFkE62z5CJdXuex8iIiJPcjmc7NmzB2lpaUhLSwMAZGVlIS0tDQsXLgQAFBcXO4IKYNt/5ssvv8TXX3+NkSNH4rXXXsNbb72FSZMm9dAlKJTG3nJS79a3iQ2ztZ6UGhlOiIjIN7g85mTixInoamkU+7Lu7Z/zww8/uPpW3k1rCw0wuzc0xIUG4OS5OpTWMJwQEXmSEALNzc2wWCxylyIbtVoNjUbT48t8KHK2jk+wt5w0u2/MCWALJwBQYmy6yJlERNRTTCYTiouLUV/v3tZxb2AwGBAfHw+dTtdjr8lw4i6OAbHubdFwdOtwzAkRkUdYrVbk5+dDrVYjISEBOp3OLxcIFULAZDLh3LlzyM/Px4ABA1xaaK0rDCfu4ujWcW+qbm05YTghIvIEk8kEq9WKxMREGAwGucuRVWBgILRaLU6dOgWTyYSAgIAeeV23TyX2Wxr3z9YBWqcTc7YOEZFn9VQrgbdzx58D/2TdxVMDYtmtQ0REPobhxF20LU19HhoQW1bTBIu12xtMExERKYZfh5MuZkRfPo1nWk6ig3VQSYDFKlBRxxk7RETUfSkpKVi+fLncZfh3OHErrWcWYdOoVegVYtuHiINiiYioK1u2bMGUKVOQkJAASZKwfv16uUvqFMOJu9hbTtw8IBYA+kTYupCKKt3bhURERN6trq4OI0eOxIoVK+QupUucSuwuHtj4zy45yoDcU+dRUFHn9vciIqKOhBBoMHt+pdhArdqlNVZuvvlm3HzzzV2eU19fj/vvvx+ffPIJIiIi8Nxzz+Ghhx663FJdwnDiLh7Y+M8uOTIIAHCK4YSISBYNZguGLPzK4+976IVJMOh69qv8tddewx/+8Ac888wz+Mc//oE5c+bguuuuw6BBg3r0fbrCbh138dDGfwCQEm3r1imo4DLKRER0eW655RY88sgj6N+/P+bPn4/o6Ghs2rTJozWw5cRdPLTOCQAkR9laTgoZToiIZBGoVePQC5Nked+eNmLECMd9SZIQFxeHsrKyHn+frjCcuIuHNv4DgJQoW8tJSXUjGkwWBOp6/j9WIiK6MEmSerx7RS5ardbpd0mSYLVaPVoDu3XcxUMb/wFAuEGHsEDbf0yFlWw9ISIi78Zw4i5az00lBlpbTzhjh4iILqS2thZ5eXnIy8sDAOTn5yMvLw+FhYXyFtYOw4m7eHBALNA67oQzdoiI6EL27NmDtLQ0pKWlAQCysrKQlpaGhQsXylyZM9/oIFMiDw6IBdq2nLBbh4iIOjdx4kSILvZuKSgo6HDM3sriSWw5cZe2G/+5dRMfmyS2nBARkY/w63DikY3/hBWwmN34RjaOlpNytpwQEZF38+twYufK0r+XzD5bB/DIuBP7mJNiYwOamj2/hDIREVFPYThxF7UOQEvo8cCMnehgHYJ0algFcPo8NwAkIiLvxXDiLpLk0c3/JEnijB0iIvIJDCfu5MHN/4A2e+xw3AkREXkxhhN30niu5QTgWidEROQbGE7cybHWiYfCSSTXOiEiIu/HcOJOHtz8D2DLCRER+QaGE3fy4OZ/QOuYk9PnG2C2eHYHSSIi8g2SJGH9+vWy1sBw4k4e3vwvNiQAeo0KzVaBs1WcTkxERM6WLFmCq6++GiEhIYiJicHtt9+Oo0ePyl1WBwwn7uThzf9UKgnJLSvFnuK4EyIiaue7777D3LlzsXPnTnz99dcwm8248cYbUVenrOEA3PjPnTy8+R9gG3fyU2lty7iTXh57XyIivyaEx/5H1InWYFtX6xJt3LjR6fc1a9YgJiYGubm5uPbaax3Hy8vLcccdd+Crr75C79698dprr+HWW2/tsbIvhuHEndpu/uch3J2YiEgG5nrgxQTPv+8zZwFdULefbjQaAQCRkZFOxxcvXoyXX34Zr7zyCl5//XVMnz4dp06d6nCeu7Bbx5008rScAJyxQ0REXbNarXjssccwbtw4DBs2zOmx++67D9OmTUP//v3x4osvora2Frt37/ZYbWw5cSetZ6cSA3CMOWHLCRGRB2kNtlYMOd63m+bOnYsDBw5g27ZtHR4bMWKE435QUBBCQ0NRVlbW7fdyFcMJHNvz9TyNZxdhA4CUlpaTwop6WKwCapXbro6IiOwk6bK6Vzxt3rx5+OKLL7Blyxb06dOnw+Nardbpd0mSYLV6bokKduu4kwc3/rOLDwuAVi3BZLGipNpz3UlERKR8QgjMmzcPn332Gb799lukpqbKXVKnGE7cycMb/wGARq1CYkTLdOJyjjshIqJWc+fOxfvvv48PP/wQISEhKCkpQUlJCRoalLU2ll+HEwHh3jeQoVsHaB13cqqS406IiKjVm2++CaPRiIkTJyI+Pt5xW7dundylOeGYE3eSoVsHsM/YOYcCztghIqI2hLj4/5R3dk5VVZUbqrkwv245cTuNZ5evt7OvdXKqnC0nRETkfRhO3Mk+xUuWlhOw5YSIiLwSw4k7eXjjP7u2++tcShMeERGRkjCcuJNGnjEnfSIMUElAg9mCczVNHn1vIiKiy8Vw4k5aeWbr6DQq9I6wBSOuFEtE5B5smbZxx58Dw4k7ybDxn10Kx50QEbmFffXU+nr+zx/Q+ufQflXZy8GpxO4kw8Z/dslRBmw9xg0AiYh6mlqtRnh4uGOvGYPBAEnyv61ChBCor69HWVkZwsPDoVare+y1GU7cqe3Gf0LY9l7wkBTH7sRM9kREPS0uLg4APLoZnlKFh4c7/jx6CsMJ3JgZ7C0ngG3Gjj2seEAywwkRkdtIkoT4+HjExMTAbDbLXY5stFptj7aY2DGcuFPbMGJu8HA4sY13KaiogxDCL5sciYjcTa1Wu+XL2d9xQKw7qbWAqiX/eXitk6RIWzipaWzG+Xr/TfVEROR9GE7cTaa1TgK0asSH2bqVOGOHiIi8CcOJu8m0SizQdqVYhhMiIvIeDCfuJlPLCdBmrRNuAEhERF6E4cTdtPKFk9YZO2w5ISIi78Fw4m4yduukOGbssOWEiIi8h1+HE49siyBjt05SSzgprGQ4ISIi7+HX4cQjZNr8D2jt1qmsM8HYwOnERETkHRhO3E3Gzf+C9RpEB+sBAIXs2iEiIi/BcOJuMm7+B7Qdd8JBsURE5B0YTtyt7eZ/MuCMHSIi8jYMJwAkuHHfGcW0nLBbh4iIvAPDibs51jmRJxwkR7PlhIiIvAvDibs5unXYckJERHQpuhVOVqxYgZSUFAQEBCA9PR27d+/u8vzly5dj0KBBCAwMRGJiIh5//HE0NsrzZe1xGvmmEgNAcqSt5eRcTRPqTc2y1EBEROQKl8PJunXrkJWVhUWLFmHv3r0YOXIkJk2ahLKysk7P//DDD/H0009j0aJFOHz4MN5++22sW7cOzzzzzGUX7xVkbjkJM2gRbtACAE6x9YSIiLyAy+Fk2bJlmD17NmbNmoUhQ4Zg5cqVMBgMWL16dafn79ixA+PGjcM999yDlJQU3HjjjZg2bdpFW1t8hswtJwBn7BARkXdxKZyYTCbk5uYiMzOz9QVUKmRmZiInJ6fT51xzzTXIzc11hJGTJ09iw4YNuOWWWy74Pk1NTaiurna6eS37ImwyhhOOOyEiIm+iceXk8vJyWCwWxMbGOh2PjY3FkSNHOn3OPffcg/LycowfPx5CCDQ3N+Phhx/usltnyZIlWLx4sSulKZeMG//ZseWEiIi8idtn62zevBkvvvgi3njjDezduxeffvopvvzyS/zhD3+44HMWLFgAo9HouBUVFbm7TPeRceM/O0fLSTlbToiISPlcajmJjo6GWq1GaWmp0/HS0lLExcV1+pznn38e9957Lx588EEAwPDhw1FXV4eHHnoIzz77LFSqjvlIr9dDr9e7UppyseWEiIjIJS61nOh0OowePRrZ2dmOY1arFdnZ2cjIyOj0OfX19R0CiFqtBgAIIVyt1/s4xpzI12qR3NJyctbYiEazRbY6iIiILoVLLScAkJWVhZkzZ+Kqq67CmDFjsHz5ctTV1WHWrFkAgBkzZqB3795YsmQJAGDKlClYtmwZ0tLSkJ6ejuPHj+P555/HlClTHCHFp8m8fD0ARAXpEKzXoLapGUWV9RgQGyJbLURERBfjcjiZOnUqzp07h4ULF6KkpASjRo3Cxo0bHYNkCwsLnVpKnnvuOUiShOeeew5nzpxBr169MGXKFPzf//1fz11FN3mk3Ubmjf8AQJIkJEcZcPBsNU5VMJwQEZGyuRxOAGDevHmYN29ep49t3rzZ+Q00GixatAiLFi3qzlt5hOTGff+U0HICAClRQTh4thoFHHdCREQKx7113M3ecmJpAqzyjfewjzvhKrFERKR0DCfuZg8ngKwzdlJaZuyw5YSIiJSO4cTdNG3CiYxdO2w5ISIib8Fw4m4qFaDW2e7LOCg2JdrWcnL6fD1MzVbZ6iAiIroYhhNPcKwSK1/LSUyIHgFaFawCOFMlX0giIiK6GIYTT7CPO5FxITZJkpAcyXEnRESkfAwnnqCAJeyB1nEnhRx3QkRECsZw4gkK2PwPaB13wpYTIiJSMoYTT1BYywln7BARkZIxnHiCAjb/A7jWCREReQf/Diee2hRZIUvY21tOiirrYbH6wY7QRETklfw7nHiKAjb/A4D4sEBo1RLMFoGznE5MREQKxXDiCQppOVGrJCRG2lpP8svZtUNERMrEcOIJjgGx8rdWDIoNAQAcLq6WuRIiIqLOMZwAkNz9Bo4BsfKHk2G9wwAAB84ynBARkTIxnHiCo1tHOeHk4BmjzJUQERF1juHEExwDYuUdcwIAQxNCAQAny+tQ02iWuRoiIqKOGE48QSEDYgEgOliP+DBbPYfYtUNERArEcOIJClmEzW5oAsedEBGRcjGceIJClq+3G85xJ0REpGAMJ56gkI3/7Ib1to072c9wQkRECuTX4UR4av16hbWc2GfsnDhXi3pTs8zVEBEROfPrcOIxClrnBABiQwPQK0QPqwAOF9fIXQ4REZEThhNPUNA6J3bDWqYUHzzLrh0iIlIWhhNPUMjGf23Zu3b2n2Y4ISIiZWE48QQFrXNix+nERESkVAwnnqCgFWLthvexhZNjpTVoNFtkroaIiKgVw4kn2MOJQhZhA4CEsABEGLRotgr8VMpBsUREpBwMJ55g79axNgMWZUzdlSSpddwJ1zshIiIFYTgBIElufgN7ywmgyEGxB85w3AkRESkHw4kn2FtOAEUNih3WMiiW04mJiEhJGE48QZJal7BXVMuJba2TI8U1MDVbZa6GiIjIxq/DifDQ6vUAWpewV9BCbEmRBoQEaGCyWHGsjINiiYhIGfw6nHiUwjb/A1oGxdq7djjuhIiIFILhxFMUtvmfnb1r5wDHnRARkUIwnHiKwjb/s+N0YiIiUhqGE09R4OZ/QGs4OVxcjWYLB8USEZH8GE48RYGb/wFAalQQgnRqNJqtOFleJ3c5REREDCceo8DN/wBApZIwJKFl3Am7doiISAEYTjzFMSBWWS0nQOsOxRx3QkRESsBw4imOAbHKajkBgOG9OZ2YiIiUg+HEUxQ6IBZoHRR78KwRVqsnV6YjIiLqiOHEUxQ6IBYA+vUKQoBWhTqTBfkVHBRLRETy8utw4tE2AoUOiAUAjVqFK+I5KJaIiJTBr8NJK8n9b2Efc6LAlhOg7Q7FHHdCRETyYjjxFAVu/NeWYxl7tpwQEZHMGE48RYEb/7Vln0584IwRwqPbNRMRETljOPEUXZDtp6lW3jouYGBsCHRqFaobm1FUqcwARURE/oHhxFMCw20/G6rkrOKCdBoVBsWFAOAOxUREJC+GE08JjLD9bKyStYyucNwJEREpAcOJpwSE234qtOUEaF2MjcvYExGRnBhOPMXRrXMeUOiA07bTiTkoloiI5MJw4in2lhNhUeyg2EFxIVCrJFTWmVBsVN5icURE5B8YTjxFGwio9bb7Cu3aCdCqMSAmGAC7doiISD4MJ54iSc5dOwrVukMxwwkREcmD4cST7F07ip6x07IYG5exJyIimTCceJJ9OrFCu3UATicmIiL5MZzA1uPiEV7QrXNFfChUElBW04Syag6KJSIiz2M48SQv6NYx6DTo18s2KJYrxRIRkRwYTjxJ4UvY2znGnZzhuBMiIvI8hhNP8oIl7AFgaIJt3AmnExMRkRwYTjzJsYS9csecAJxOTERE8vLrcOLxJdq9pFtnSEvLyVljIypqm2SuhoiI/E23wsmKFSuQkpKCgIAApKenY/fu3V2eX1VVhblz5yI+Ph56vR4DBw7Ehg0bulWwV/OSbp2QAC1So4MA2PbZISIi8iSXw8m6deuQlZWFRYsWYe/evRg5ciQmTZqEsrKyTs83mUz42c9+hoKCAvzjH//A0aNHsWrVKvTu3fuyi/c6XtKtA3DcCRERycflcLJs2TLMnj0bs2bNwpAhQ7By5UoYDAasXr260/NXr16NyspKrF+/HuPGjUNKSgquu+46jBw58rKL9zpe0q0DtBl3wunERETkYS6FE5PJhNzcXGRmZra+gEqFzMxM5OTkdPqczz//HBkZGZg7dy5iY2MxbNgwvPjii7BYLBd8n6amJlRXVzvdfIJjnRMjYLXKWsrFcDoxERHJxaVwUl5eDovFgtjYWKfjsbGxKCkp6fQ5J0+exD/+8Q9YLBZs2LABzz//PF577TX88Y9/vOD7LFmyBGFhYY5bYmKiK2Uql73lBAJoUvaXvr1bp7CyHsZ6s8zVEBGRP3H7bB2r1YqYmBj8/e9/x+jRozF16lQ8++yzWLly5QWfs2DBAhiNRsetqKjI3WV6hkYPaA22+wofdxJu0CExMhAAu3aIiMizNK6cHB0dDbVajdLSUqfjpaWliIuL6/Q58fHx0Gq1UKvVjmNXXHEFSkpKYDKZoNPpOjxHr9dDr9e7Upr3CAgHzPWKn7EDAMMSwlBU2YADZ424pn+03OUQEZGfcKnlRKfTYfTo0cjOznYcs1qtyM7ORkZGRqfPGTduHI4fPw5rmzEWP/30E+Lj4zsNJnLw1L5/ALxiZ2I7jjshIiI5uNytk5WVhVWrVuHdd9/F4cOHMWfOHNTV1WHWrFkAgBkzZmDBggWO8+fMmYPKyko8+uij+Omnn/Dll1/ixRdfxNy5c3vuKryJF+xMbGcfd3KA04mJiMiDXOrWAYCpU6fi3LlzWLhwIUpKSjBq1Chs3LjRMUi2sLAQKlVr5klMTMRXX32Fxx9/HCNGjEDv3r3x6KOPYv78+T13Fd7EC3YmtrO3nJwsr0NNoxkhAVqZKyIiIn/gcjgBgHnz5mHevHmdPrZ58+YOxzIyMrBz587uvJXv8aK1TqKD9YgPC0CxsRGHi2swJjVS7pKIiMgP+PXeOrLwkiXs7YYm2MedsGuHiIg8g+HE07xoCXugdaVYhhMiIvIUhhNP86JuHQAY1rtlUCzXOiEiIg9hOPE0L+vWsQ+KPV5WiwbThbccICIi6ikMJ57mZd06MSF6RAfrYRXA4RKud0JERO7HcOJpjm4d7+gmkSQJw3tzvRMiIvIchhNP86J1TuyGcVAsERF5ULfWOfEVCeGBqGowI0jvwT8G+5iTpmrA0gyolf8R2KcT/1BYJW8hRETkF5T/zehGf73nSs+/aUBY6/1GIxAU5fkaXDS2byRUEnCsrBZnqhrQOzxQ7pKIiMiHsVvH09QaQBdiu+8lXTvhBh2uTLK1+Gw6UiZzNURE5OsYTuTgRTsT210/OAYAsPkowwkREbkXw4kcAlu6drxkOjEAXD/IFk62H69Ao5nrnRARkfswnMjBC2fsXBEfgthQPRrMFuzKr5S7HCIi8mEMJ3JwrHXiPS0nkiQ5Wk847oSIiNyJ4UQOXjjmBAAmDuK4EyIicj+GEzl4YbcOAIwfEA2tWkJBRT3yy+vkLoeIiHwUw4kcvGxnYrtgvQZjUiMBAN+ya4eIiNyE4UQOXrYzcVvXs2uHiIjcjOFEDl62M3Fb9nEnu05Woq6pWeZqiIjIFzGcyMFLu3UAoF+vICRFGmCyWLHjRIXc5RARkQ9iOJGDlw6IBexTinsB4LgTIiJyD4YTOTimEntftw4ATGyzlL0QQuZqiIjI1zCcyMHerWOuB5pNspbSHRl9oxCgVaHY2IijpTVyl0NERD6G4UQO+jAAku2+F3btBGjVuKZfNAB27RARUc9jOJGDSgUE2Df/q5K1lO6yjzvZfOSczJUQEZGvYTiRixfur9OWfUpxbuF5GOvNMldDRES+hOFELl48YwcAEiMNGBATDItVYOtxtp4QEVHPYTiRixevdWJ3fcusHY47ISKinsRwIhcvn04MABNbxp18d/QcrFZOKSYiop7BcCIXL+/WAYCrkiMRrNegos6E/WeMcpdDREQ+guFELj7QraPTqDBhgG1K8SZuBEhERD2E4UQuPtByArTuUryJ406IiKiHMJzIxQfGnADAdS3jTn48bcS5miaZqyEiIl/AcCIXH+jWAYDY0AAMTQgFAGz5iVOKiYjo8jGcyMVHunUA4H9aphRz3AkREfUEhhO5+Ei3DtC6WuyWn86h2WKVuRoiIvJ2DCdy8ZFuHQAYlRiOcIMW1Y3N2FtYJXc5RETk5RhO5GLv1rE0AeYGWUu5XGqVhOsG2gbGsmuHiIguF8OJXPQhgKS23feBrh3HuBNOKSYiosvEcCIXSfKprp1rB/SCJAFHSmpwtsq7W4KIiEheDCdy8qEZOxFBOqQlhgMANh/llGIiIuo+hhM5+VDLCdBmtViOOyEiosvAcCInH5pODADXt4w72X68HE3NFpmrISIib8VwIicf6tYBgKEJoYgJ0aPeZMHu/Eq5yyEiIi/FcCInH+vWkSQJE1v22tl0hONOiIioexhO5GRvOfGRbh2gddzJZo47ISKibmI4kZN9zImPdOsAwPgB0dCoJJwsr0NBeZ3c5RARkRdiOJGTj3XrAEBIgBZXp0QC4KwdIiLqHoYTOfnYgFi76wfbl7LnuBMiInIdw4mcfGwqsZ193MnOkxWoNzXLXA0REXkbhhM5+WC3DgD0jwlGn4hAmJqtyDlRIXc5RETkZRhO5NS2W0cIOSvpUZIkOVpPvuVGgERE5CKGEznZW06szYCpVtZSepp93Mnmo+cgfCh4ERGR+zGcyElrANQ6230f69rJ6BsNvUaFM1UNOFbmW8GLiIjci+FETpLkszN2AnVqZPSLAgBsYtcOERG5gOFEbj46KBYAx50QEVG3MJzIzUenEwOt4WTPqfOobjTLXA0REXkLhhO5+Wi3DgAkRRnQr1cQLFaBbcfK5S6HiIi8BMOJ3Hy4Wwdg1w4REbmO4URuPrgzcVvXD7bvUnwOViunFBMR0cUxnMjNB3cmbuvqlEgE6dQor23CwbPVcpdDRERegOFEbj7eraPTqDB+QDQA7lJMRESXplvhZMWKFUhJSUFAQADS09Oxe/fuS3re2rVrIUkSbr/99u68rW/y8W4dgONOiIjINS6Hk3Xr1iErKwuLFi3C3r17MXLkSEyaNAllZV1/8RQUFODJJ5/EhAkTul2sTwq2fXGjpkTeOtxoYks4+fF0FSpqm2SuhoiIlM7lcLJs2TLMnj0bs2bNwpAhQ7By5UoYDAasXr36gs+xWCyYPn06Fi9ejL59+15WwT4nPNn2s+qUT23+11ZcWACGxIdCCODfP56VuxwiIlI4l8KJyWRCbm4uMjMzW19ApUJmZiZycnIu+LwXXngBMTExeOCBBy7pfZqamlBdXe1081nhiQAkwFwP1FfIXY3bTBuTCABYvb0AFs7aISKiLrgUTsrLy2GxWBAbG+t0PDY2FiUlnXdLbNu2DW+//TZWrVp1ye+zZMkShIWFOW6JiYmulOldNHogJN52//wpeWtxo1+OTkS4QYvCynr896DvdmEREdHlc+tsnZqaGtx7771YtWoVoqOjL/l5CxYsgNFodNyKiorcWKUCRNi7dgpkLcOdAnVq3DvWdp2rtp6UuRoiIlIyjSsnR0dHQ61Wo7S01Ol4aWkp4uLiOpx/4sQJFBQUYMqUKY5jVqvV9sYaDY4ePYp+/fp1eJ5er4der3elNO8WngQU5vh0ywkA3JuRjL99dxJ7C6uQe6oSo5Mj5S6JiIgUyKWWE51Oh9GjRyM7O9txzGq1Ijs7GxkZGR3OHzx4MPbv34+8vDzH7dZbb8X111+PvLw83+6ucUXbQbE+LCYkAHek9QYArNqSL3M1RESkVC61nABAVlYWZs6ciauuugpjxozB8uXLUVdXh1mzZgEAZsyYgd69e2PJkiUICAjAsGHDnJ4fHh4OAB2O+zV7t46Pt5wAwIMTUrFuTxG+OlSCgvI6pEQHyV0SEREpjMvhZOrUqTh37hwWLlyIkpISjBo1Chs3bnQMki0sLIRKxYVnXeJoOSmUtw4PGBAbgusH9cKmo+ewens+XriNIZWIiJxJQih/cY3q6mqEhYXBaDQiNDRU7nJ6XlUhsHw4oNYBz5YCPh7udhwvxz1v7UKAVoWcp29ARJBO7pKIiMgNuvv97dvfgt4iJAFQaQCLCagplrsat8voF4Uh8aFoNFvxwS7f78oiIiLXMJwogVoDhNoGivr6oFgAkCQJD11rWyl4zY5TaDRbZK6IiIiUhOFEKfxoUCwATB4Rj/iwAJTXNuHzPC5pT0RErRhOlMJPphPbadUqzBqXAsC2KJsXDH0iIiIPYThRCj9rOQGAu8ckIVivwbGyWmz+6Zzc5RARkUIwnChFeIrtpx9MJ7YLDdDi7qttC/Gt2sIl7YmIyIbhRCki/Ktbx27W+FSoVRJ2nKjAwbNGucshIiIFYDhRCvuYk+ozgMUsby0e1Ds8EJOH23Zlfmsrl7QnIiKGE+UIjgE0AYCwAkYf34W5ndkTbNOK//3jWRQbG2SuhoiI5MZwohSSZNudGPCrQbEAMLxPGMb2jUSzVWDN9gK5yyEiIpkxnCiJn00nbsveevLhrkLUNPpPtxYREXXEcKIkEf6zAWB71w+KQd9eQahpasa67/2rW4uIiJwxnChJuP+tdWKnUkmO1pN3theg2WKVuSIiIpILw4mS+Ol0Yrs70nojKkiHM1UN2HCgRO5yiIhIJgwnSuKnA2LtArRqzMhIAWBblI1L2hMR+SeGEyWxd+vUlQGmenlrkcmvxyZBr1Fh/xkjduVXyl0OERHJgOFESQIjAH2o7b4fDooFgKhgPX45ug8ALmlPROSvGE6URJLaTCf2z3ACAA+MT4UkAdlHynC8rFbucoiIyMMYTpTGzwfFAkDfXsHIvCIWAPD2NraeEBH5G4YTpXFMJy6QtQy52acV/3PvGZTXNslcDREReRLDidLYZ+z4ccsJAFydEoGRieEwNVvxXo5//1kQEfkbhhOlifDfhdjakiQJsyekAgDe33kKjWaLzBUREZGnMJwojR/vr9PeTUPj0CciEJV1Jvxz72m5yyEiIg9hOFEae7dOoxFoqJK1FLlp1CrcP87WevL21nxYrVyUjYjIHzCcKI0+GDBE2+778XRiu7uuTkRIgAYny+uQfaRM7nKIiMgDGE6UiNOJHYL1GkxPt/15cFE2IiL/wHCiRH6+x057912TAo1Kwu6CSuQVVcldDhERuRnDiRJxUKyTuLAA3DoqAQCwaitbT4iIfB3DiRJxOnEHD463Lcr2n/3FKKr0z00RiYj8BcOJEnF/nQ6GJIRiwoBoWAWwenu+3OUQEZEbMZwoUUSK7WfVKUBw+qzdgy1L2q/7vgjGerPM1RARkbswnChRWB8AEmCuB+rK5a5GMa4dEI1BsSGoN1nw4W62KhER+SqGEyXS6IGQeNt9Dop1kCQJD7Ysab9mRz5MzVaZKyIiIndgOFGqCO5O3JlbRyUgJkSP0uom/PvHs3KXQ0REbsBwolScTtwpvUaNmdekALBNKxYck0NE5HMYTpQqgjN2LmR6ehICtWocKanB5qPn5C6HiIh6GMOJUoVzrZMLCTfocE+6bRXdBZ/uR2WdSeaKiIioJzGcKBX31+lS1s8Gom+vIJRUN+KpT35k9w4RkQ9hOFEqx5iTIsBqkbcWBQrSa/D6tDToNCpkHynD6u0FcpdEREQ9hOFEqUITAJUGsJqBmmK5q1GkoQlheG7yFQCApf85jH2nq+QtiIiIegTDiVKp1C2LsYHjTrpw79hkTBoaC7NF4Lcf/YCaRq4cS0Tk7RhOlIzTiS9KkiS8fOdI9A4PxKmKejzz2QGOPyEi8nIMJ0rG6cSXJMygxV+mpUGtkvDvH8/i4z1FcpdERESXgeFEyTid+JKNTo7AEzcOBAAs+vwgfiqtkbkiIiLqLoYTJWu7OzFd1MPX9sOEAdFoNFsx78O9aDBxlhMRkTdiOFGycNtCY2w5uTQqlYRld41CdLAeP5XW4oUvDsldEhERdQPDiZLZu3WqzwDNXAX1UvQK0WP51FGQJOCj3YX4Yh83ByQi8jYMJ0oWHANoAgEIDop1wfgB0XhkYj8AwIJ/7kdhRb3MFRERkSsYTpRMkoC4Ybb7p3fLW4uXeTxzIK5KjkBNUzPmfbQXpmar3CUREdElYjhRupTxtp8F2+Wtw8to1Cr8eVoawgK12HfaiFe+OiJ3SUREdIkYTpQu2R5OtspbhxfqHR6Il385AgCwams+vj1SKnNFRER0KRhOlC4pHZDUtunExtNyV+N1Jg2Nw33XpAAAnvj4R5QYG+UtiIiILorhROn0IUDCKNt9du10y4JbBmNoQijO15vx6NofYLFyeXsiIiVjOPEGyeNsP9m10y16jRqvT0uDQafGrvxKvP7tMblLIiKiLjCceAP7oNhTbDnprr69gvF/d9hmPv0l+xh2nqyQuSIiIroQhhNvkDQWkFRA5UmgmouKddcdaX3wy9F9YBXAo2t/QEVtk9wlERFRJxhOvEFAGBBnm3XCcSeXZ/GtQ9G3VxBKq5vw5Cc/QgiOPyEiUhqGE2/h6NrZJm8dXi5Ir8GKe66ETqPCpqPn8Pa2fLlLIiKidhhOvIVjMTaGk8t1RXwonv/5EADASxuP4MeiKnkLIiIiJwwn3iIpA4AEVBwHakrkrsbr/To9CTcPi4PZIvDbj35AdaNZ7pKIiKgFw4m3CAwH4obb7nPWzmWTJAlL7xyB3uGBKKysxzOf7uf4EyIihWA48Sbs2ulRYYFavH5PGjQqCV/sK8ba74vkLomIiMBw4l0ci7Gx5aSnXJkUgScnDQIA/O/nB/FTaY3MFRERUbfCyYoVK5CSkoKAgACkp6dj9+7dFzx31apVmDBhAiIiIhAREYHMzMwuz6cuJF8DQALKjwK1ZXJX4zMemtAX1w7shaZmK+Z+sBcNJovcJRER+TWXw8m6deuQlZWFRYsWYe/evRg5ciQmTZqEsrLOvyw3b96MadOmYdOmTcjJyUFiYiJuvPFGnDlz5rKL9zuGSCB2qO0+x530GJVKwrK7RqJXiB7Hymqx+N8H5S6JiMivuRxOli1bhtmzZ2PWrFkYMmQIVq5cCYPBgNWrV3d6/gcffIBHHnkEo0aNwuDBg/HWW2/BarUiOzv7sov3S+zacYvoYD2WTx0FSQLWfl+Ez3/kSrxERHJxKZyYTCbk5uYiMzOz9QVUKmRmZiInJ+eSXqO+vh5msxmRkZEXPKepqQnV1dVON2rBQbFuM65/NOZO7A8AeObT/ThVUSdzRURE/smlcFJeXg6LxYLY2Fin47GxsSgpubS1N+bPn4+EhASngNPekiVLEBYW5rglJia6UqZvs7ecnDsM1HHzup72WOYAXJ0SgdqmZvz2ox9garbKXRIRkd/x6GydpUuXYu3atfjss88QEBBwwfMWLFgAo9HouBUVcYqnQ1AUEGNb3ZTjTnqeRq3Cn+9OQ1igFvtOG7HkP4e5/gkRkYe5FE6io6OhVqtRWlrqdLy0tBRxcXFdPvfVV1/F0qVL8d///hcjRozo8ly9Xo/Q0FCnG7XhGHfCrh13SAgPxKu/GgkAeGd7AZ7+5340mjmDh4jIU1wKJzqdDqNHj3YazGof3JqRkXHB57388sv4wx/+gI0bN+Kqq67qfrVkk9ISTthy4jY/GxKL5yZfAZUErNtThKl/y8HZqga5yyIi8gsud+tkZWVh1apVePfdd3H48GHMmTMHdXV1mDVrFgBgxowZWLBggeP8l156Cc8//zxWr16NlJQUlJSUoKSkBLW1tT13Ff4muWVQbOkBoL5S3lp82IMT+uLd+8cg3KDFj6eNmPL6Nuw8yXE+RETu5nI4mTp1Kl599VUsXLgQo0aNQl5eHjZu3OgYJFtYWIji4mLH+W+++SZMJhN++ctfIj4+3nF79dVXe+4q/E1wLyDatqopTu2QtxYfN2FAL/x73ngMiQ9FRZ0J09/ahdXb8jkOhYjIjSThBf/KVldXIywsDEajkeNP7L54HNizGhj7CHDTErmr8XkNJgsWfLoP6/Ns65/cPioBS34xAoE6tcyVEREpV3e/v7m3jrdyrHeyVd46/ESgTo0/TR2FhT8fArVKwvq8s7jzzR0oqqyXuzQiIp/DcOKt7ONOSg4ADeflrcVPSJKE+8en4oMH0xEVpMOh4mpM+es2bPnpnNylERH5FIYTbxUSC0QNACCAwp1yV+NXxvaNwhe/G4+RieGoqjfjvnd2483NJzgOhYiohzCceLMUrncil/iwQKx7aCymXpUIqwBe2ngEcz/ci9qmZrlLIyLyegwn3iyZ++zIKUCrxtI7h+P/7hgGrVrChv0luGPFduSXc08eIqLLwXDizewtJyX7gEajvLX4KUmSMD09GWsfykBMiB7Hympx6+vbkH249OJPJiKiTjGceLPQBCCyLyCsHHcis9HJEfjit+NxVXIEapqa8cC7e7D8m59gtXIcChGRqxhOvB332VGMmNAAfDh7LGZkJAMAln9zDA/9vz2objTLXBkRkXdhOPF2KRNsPxlOFEGnUeGF24bh1V+NhE6jwjeHy3DbX7fjWGmN3KUREXkNhhNvZx93Uvwj0Fgtby3k8MvRffDPh69B7/BA5JfX4fYV2/Gf/cUXfyIRETGceL2wPkBEKiAswMHP5K6G2hjeJwyfzxuHa/pFoc5kwZwP9uKljUdg4TgUIqIuMZz4gqsfsP3c9ifAwnU2lCQqWI/37h+D2RNSAQBvbj6B+97Zjap6k8yVEREpF8OJL7jqfiAwEjifDxz8VO5qqB2NWoVnJw/BX6alIVCrxtZj5Zjy1204dJbdcEREnWE48QW6ICBjru3+llcBq1XeeqhTt45MwKePXIOkSAOKKhvwize34195Z+Qui4hIcRhOfMWY2YA+DCg/Chz5t9zV0AVcER+Kz+eNw3UDe6HRbMWja/Pwwr8PwWxhoCQismM48RUBYUD6b2z3t7wCcBM6xQo36LD6vqsx7/r+AIDV2/Px67d2oby2SebKiIiUgeHEl4ydA2iDgJL9wLH/yl0NdUGtkvDkpEH4272jEazXYFd+Jaa8vg0/FlXJXRoRkewYTnyJIbJ15s53L7P1xAtMGhqH9XPHoW+vIBQbG/Grv+Xg4++L5C6LiEhWDCe+JmMeoAkAzuwB8r+Tuxq6BP1jgvGvueNw45BYmJqt+P0/92H+P/ah2Nggd2lERLJgOPE1IbHAlTNt97e8Km8tdMlCArRY+evReOJnAyFJwLo9RZjw0iY8vi4PB85wx2ki8i+SEMpv+6+urkZYWBiMRiNCQ0PlLkf5jKeBP48CrGbg/q+ApLFyV0Qu2HGiHH/+5hh25Vc6jmX0jcKDE1Jx/aAYqFSSjNUREV267n5/M5z4qs9/B+x9F+j/M+DX/5C7GuqGfaer8Pa2fHyxr9ix5H3fXkF4YHwq7ryyDwK0apkrJCLqGsMJOavMB14fbdtz56HNQEKa3BVRN52tasC7Owrw4e5C1DTatieIDNLh1+lJuDcjBb1C9DJXSETUOYYT6ujT3wD71gKDfw7c/YHc1dBlqm1qxsffF2H19nycPm8bLKtTq3B7WgIenNAXA2NDZK6QiMgZwwl1dO4osCIdgADm5ACxQ+SuiHpAs8WK/x4qxaqtJ/FDYZXj+LUDe2H2hFSM7x8NSeK4FCKSH8MJde7jmcCh9cCwXwK/fFvuaqiH5Z46j7e2nsRXB0vQMiwFg+NC8MD4VNw6KgF6DcelEJF8GE6ocyX7gZXjAUkFzP0eiO4vd0XkBoUV9Vi9PR8f7ylCvckCAOgVosfMjGRMT09GRJBO5gqJyB8xnNCFfXg38NN/gFG/Bm5fIXc15EbGBjM+2l2INdsLUFLdCAAI0Krwy9F9cP+4VPTtFSxzhUTkTxhO6MJO7wHeugFQaYDf7gUikuWuiNzM1GzFhv3FWLX1JA6erQYASBJww+BYPDghFempkRyXQkRux3BCXXvvduDkJuCqB4CfL5O7GvIQIQR2nqzE29tO4pvDZY7jw3uH4cEJqbhleDy0ai4UTUTuwXBCXSvYDqy5BVDrgEf3AaHxcldEHnbiXC3e3paPf+aeRlOzFQAQHxaA+65Jwd1jkhAWqJW5QiLyNQwndHGrbwYKdwBj5wI3vSh3NSSTyjoT3t95Cu/lFKC81gQACNKpcdfVibh/XCoSIw0yV0hEvoLhhC7ueDbw/i8ATSDw2z1AWB+5KyIZNZot+DzvLN7adhI/ldYCAFQScNOwODw4oS+uTIqQuUIi8nYMJ3RxQgBvZQJn9gBRA4BZG4DgGLmrIpkJIbDlWDne2noSW4+VO45fmRSO2RP64sahcVBzs0Ei6gaGE7o0VYXAO7cAxiIgZihw3xeAIVLuqkghjpRU4+2t+fhX3lmYLLZxKYmRgbh/XCp+cWUfjkshIpcwnNClqzhhCyi1JUD8SGDG50BguNxVkYKU1TTi/+Wcwv/beQpV9WYAti6fIQmhGJMShfS+kRiTEsnF3YioSwwn5JpzR20Bpb4c6HM1cO9ngJ4bx5GzBpMF/9x7Gu/lFDjGpbQ1OC4E6amRGJMahTGpkdwhmYicMJyQ60oOAGsmA41VQPJ4YPongI4zNahzpdWN2JVfiV0nK7ArvxLHyzqGlX69gpDeNwrpqZFIT41CXFiADJUSkVIwnFD3nNkLvHcb0FQN9Psf4O6PAC2/UOjiymub8H1+JXblV2LnyQocKanpcE5ylMERVNL7RqJPBMMvkT9hOKHuK9wF/L87AHMdMPBm4K73AA3HEpBrqupN2J1fid0tgeXgWaNjp2S73uGBtrDS1xZYkqMMXEafyIcxnNDlyd8CfPAroLkRGHIbcOdqQK2RuyryYtWNZuQWnLd1BeVXYP9pI5rbpZXYUD3GpNq6gcb2jUS/XsEMK0Q+hOGELt+xb4C10wCLCRgxFbh9JaDivivUM+qamrG38Dx2nbSFlR+LjI7pynbRwTqMSbXNBErvG4VBsSFQcY0VIq/FcEI948iXwMczAGszcOVMYMqfbdvZEvWwRrMFPxRWYVd+BXadrMTewvOOPX/swg1aXJ0S6Ri3MiQhlAvCEXkRhhPqOQc+Bf75ACCswJjfADe/xIBCbtfUbMH+00bHANvcU+dRb7I4nROi1+CqlAjHjKBhvcO4qzKRgjGcUM/K+whY/7Dt/rhHgczFDCjkUWaLFQfPVjumLn+fX4mapmancww6NUYnR7QMso3CiD5h0GvUMlVMRO0xnFDP27Ma+OJx2/2JC4CJT8tbD/k1i1XgcHG1Y62V3QWVjtVr7fQaFdKSwh1Tl6+IC0W4QctBtkQyYTgh99j5JrCxJZSMfQS47vdAIHerJflZrQI/ldVg10n79OUKlNeaOpwXotcgMdKAxMhAJEUakBRpQGLLz94RgWxpIXIjhhNyn63LgOzFtvv6MCBjLjB2DhDAz4KUQwiBE+fqsCu/ArvzK7Gn4DzOVDV0+RxJAuJDA9CnJay0Dy/RwTq2uhBdBoYTcq8jG4Bv/wCUHbL9HhgBXPM7IP03gC5I3tqILqDRbMHp8/UorKxHYUU9is43oLCyHkWVtmPtB9y2F6hVt4SVQEdgsd/6RBgQqGOrC1FXGE7I/axW4NBnwKYlQMUx2zFDNDAhC7jqfkAbKG99RC4QQqCiztQaVirqUdQSZIoqG3DW2ICL/esYE6J3hJb24SUmRM81WsjvMZyQ51iagf2fAN8tBc4X2I6FxAMTngCunAFouDMteb+mZgvOVjXaWl0q63G65ae9Fab9zKH2dBoVEiOcW1zaBplgPVdgJt/HcEKeZzEDeR8CW14BjEW2Y2GJwLVPAaPuAdRaeesjchMhBIwN5taw0tLaYu8uOlPVAEv7jYXaiQrStQzUNSCpZbCuPbzEhwVysTnyCQwnJJ/mJmDve8CWV4HaEtuxiBTguqeBEXcBKvbLk39ptlhRbGxsF15ax7qcbzcFuj2tWkLv8MA24cW59SUskMGfvAPDCcnP3GBbG2Xbn4C6c7Zj0QNt66MMuYP79BC1qG40O4WVosrWgbpF5+thtnT9z3JIgAbRwXqEBWoRbtAiPFCLcIOu9XeDFuGBOoS1eSw0QAMNV9MlD2M4IeUw1QG7/w5s/zPQcN52LGYocN1TQP+fAfpgeesjUjCLVaC0urHTFpfCygaU1zZ1+7VDAjSO4BJu0LYJN21/1zkCT1jLYzoNQw11D8MJKU9jNbBrJbDjr0CT0XZMUgMJaUDKOCBlApCYzvVSiFxQb2rGmfMNqGowo6rejKp6E4z2+w2mlmOt94315osO3r0Yg07taIGxt8yEBeratNq0+b1N2AnQskvX3zGckHI1nAdyVgD71gFVhc6PSWogfiSQMt52SxoLBITJUyeRjzJbrKhuMDsCjdEpxJhhrDe1hp02vxsbzBedTt0VvUbVSRdTuy6o9q04Bh2CdGoufucjGE7IO1QVAgXbgVPbgIJtrVOR7SQVEDeiJaxMsIWVwHA5KiXye1arQE1jc2uLTEO7lpqWFhpjJ481X2S2Ulc0Ksmpmymis5Yagw7hgVqEBmph0KkRoFEjQKdCgFaNQK2au1UrBMMJeSfjaVtYKdgKnNoOVJ5sd4IExI8AkltaVpIzuLcPkcIJIVBnsqCqvqVrqV23k7ElyLS21NgeO19vhqnZ2iM1qFUSArVqBGhtgcUeWtr+bjvWGmj0bc4JbHNO298D2wWhAI0aWrXElp4LYDgh32A8YwspBS0tK5Un2p0gAXHDWlpVMoDIvkBYH7auEPmIBpPFaeyM0anVpmOXVHWDGY1mCxrMFjSaLbiMBptuU6skBGhUCNSpode0BJgOAefCQaiz4BTYPlDpvDMIMZyQb6oudg4r9mXz29OH2kKK45bYcmv5PSQeUHNFTiJfJoSAyWJFo9mKxpawYgstVjSYLGhstqCx5WeDyep4vKnteS3Pa2z5vW3waf+7EoOQ41jLOZ22/lwgCCVF9vx+UQwn5B9qSlrDyplcoKoIaKi8+PMkNRCa0HWA4awhIrpEPRuE2r+GPEHok4czcHVKZI++Zne/v/m/kuRdQuKAYXfabnamOlt3kLHINobF8dN+/wxgNbfcL7rwa+vDgPBE5wATHGtrldGHtNza3NcGAl7UvEpEPUeSJOg1tpYKd6/Y2xNBqLHd/c6CkEFBu2x3K5ysWLECr7zyCkpKSjBy5Ei8/vrrGDNmzAXP/+STT/D888+joKAAAwYMwEsvvYRbbrml20UTOdEFAb0G2m6dsVqB2tJOgstpwFho+9lw3rYWS6kRKD1wae8rqTsGlg63zh5rd0wXzNVzieiCPBmElMLlcLJu3TpkZWVh5cqVSE9Px/LlyzFp0iQcPXoUMTExHc7fsWMHpk2bhiVLluDnP/85PvzwQ9x+++3Yu3cvhg0b1iMXQdQllQoIjbfdEq/u/JymWqC6TetLVUsrS30F0FTT5lZt+ymsgLAAjVW22+XStQswOgOg1gNqHaDR2X6qtZ0ca3NzHNO3nKuz7RDd9nlqbRfHdNwHiYgUweUxJ+np6bj66qvx17/+FQBgtVqRmJiI3/72t3j66ac7nD916lTU1dXhiy++cBwbO3YsRo0ahZUrV17Se3LMCSmKEIC53hZSGqudQ0tnQeaCx6sB6+Wt3NnjJHUnYaflptLYgp6ktq1Ho1Lb7qva/e70mKrNOW0fa/98+7F259vP7ex8ldrWrdbhmP19W26QWu5LLfelNsfbHrvA4x3ObftauMjjPfle7Z6Hli5FR9ei5Hy//WOO37t6rN3v7Laky+SRMScmkwm5ublYsGCB45hKpUJmZiZycnI6fU5OTg6ysrKcjk2aNAnr16+/4Ps0NTWhqal1/wij0bb0eXV1tSvlErmZAdAbAH1c954uhG1H56ZawFRt+9lUC5hqbeNoLCbAYgasptb7zaaOx5vNLcdajluaWn62vW9q89P+Wk2AaB+Omltu9Zf5Z0O+Ser4s32A6W4g6vAYXDi3fYjqJFR1OHShul05p7P3UVAtHY5d5DVuXwHEDu3kdbrP/r3t6twbl8JJeXk5LBYLYmNjnY7HxsbiyJEjnT6npKSk0/NLSkou+D5LlizB4sWLOxxPTEx0pVwiIiK6VL+/xm0vXVNTg7CwS9+aRJGzdRYsWODU2mK1WlFZWYmoqKgeXXymuroaiYmJKCoq8tnuIl6jb+A1+gZeo2/gNV46IQRqamqQkJDg0vNcCifR0dFQq9UoLS11Ol5aWoq4uM6btuPi4lw6HwD0ej30er3TsfDwcFdKdUloaKjP/gdmx2v0DbxG38Br9A28xkvjSouJnUvzF3U6HUaPHo3s7GzHMavViuzsbGRkZHT6nIyMDKfzAeDrr7++4PlERETk31zu1snKysLMmTNx1VVXYcyYMVi+fDnq6uowa9YsAMCMGTPQu3dvLFmyBADw6KOP4rrrrsNrr72GyZMnY+3atdizZw/+/ve/9+yVEBERkU9wOZxMnToV586dw8KFC1FSUoJRo0Zh48aNjkGvhYWFULVZUOqaa67Bhx9+iOeeew7PPPMMBgwYgPXr1ytijRO9Xo9FixZ16ELyJbxG38Br9A28Rt/Aa3Q/r9hbh4iIiPwH18wmIiIiRWE4ISIiIkVhOCEiIiJFYTghIiIiRfHrcLJixQqkpKQgICAA6enp2L17t9wldWrJkiW4+uqrERISgpiYGNx+++04evSo0zkTJ06EJElOt4cfftjpnMLCQkyePBkGgwExMTF46qmn0NzsvLfK5s2bceWVV0Kv16N///5Ys2aNuy8PAPC///u/HeofPHiw4/HGxkbMnTsXUVFRCA4Oxp133tlhcT8lXx8ApKSkdLhGSZIwd+5cAN75GW7ZsgVTpkxBQkICJEnqsGeWEAILFy5EfHw8AgMDkZmZiWPHjjmdU1lZienTpyM0NBTh4eF44IEHUFtb63TOvn37MGHCBAQEBCAxMREvv/xyh1o++eQTDB48GAEBARg+fDg2bNjg9ms0m82YP38+hg8fjqCgICQkJGDGjBk4e/as02t09tkvXbrUK64RAO67774O9d90001O53jz5wig07+bkiThlVdecZyj5M/xUr4nPPnv6GV/vwo/tXbtWqHT6cTq1avFwYMHxezZs0V4eLgoLS2Vu7QOJk2aJN555x1x4MABkZeXJ2655RaRlJQkamtrHedcd911Yvbs2aK4uNhxMxqNjsebm5vFsGHDRGZmpvjhhx/Ehg0bRHR0tFiwYIHjnJMnTwqDwSCysrLEoUOHxOuvvy7UarXYuHGj269x0aJFYujQoU71nzt3zvH4ww8/LBITE0V2drbYs2ePGDt2rLjmmmu85vqEEKKsrMzp+r7++msBQGzatEkI4Z2f4YYNG8Szzz4rPv30UwFAfPbZZ06PL126VISFhYn169eLH3/8Udx6660iNTVVNDQ0OM656aabxMiRI8XOnTvF1q1bRf/+/cW0adMcjxuNRhEbGyumT58uDhw4ID766CMRGBgo/va3vznO2b59u1Cr1eLll18Whw4dEs8995zQarVi//79br3GqqoqkZmZKdatWyeOHDkicnJyxJgxY8To0aOdXiM5OVm88MILTp9t27+/Sr5GIYSYOXOmuOmmm5zqr6ysdDrHmz9HIYTTtRUXF4vVq1cLSZLEiRMnHOco+XO8lO8JT/072hPfr34bTsaMGSPmzp3r+N1isYiEhASxZMkSGau6NGVlZQKA+O677xzHrrvuOvHoo49e8DkbNmwQKpVKlJSUOI69+eabIjQ0VDQ1NQkhhPj9738vhg4d6vS8qVOnikmTJvXsBXRi0aJFYuTIkZ0+VlVVJbRarfjkk08cxw4fPiwAiJycHCGE8q+vM48++qjo16+fsFqtQgjv/wzb/4NvtVpFXFyceOWVVxzHqqqqhF6vFx999JEQQohDhw4JAOL77793nPOf//xHSJIkzpw5I4QQ4o033hARERGOaxRCiPnz54tBgwY5fr/rrrvE5MmTnepJT08Xv/nNb9x6jZ3ZvXu3ACBOnTrlOJacnCz+9Kc/XfA5Sr/GmTNnittuu+2Cz/HFz/G2224T//M//+N0zJs+x/bfE578d7Qnvl/9slvHZDIhNzcXmZmZjmMqlQqZmZnIycmRsbJLYzQaAQCRkZFOxz/44ANER0dj2LBhWLBgAerrW7e9z8nJwfDhw512iJ40aRKqq6tx8OBBxzlt/0zs53jqz+TYsWNISEhA3759MX36dBQWFgIAcnNzYTabnWobPHgwkpKSHLV5w/W1ZTKZ8P777+P+++932szS2z/DtvLz81FSUuJUT1hYGNLT050+t/DwcFx11VWOczIzM6FSqbBr1y7HOddeey10Op3jnEmTJuHo0aM4f/684xylXLfRaIQkSR32A1u6dCmioqKQlpaGV155xamp3BuucfPmzYiJicGgQYMwZ84cVFRUONXvS59jaWkpvvzySzzwwAMdHvOWz7H994Sn/h3tqe9XRe5K7G7l5eWwWCxOHwAAxMbG4siRIzJVdWmsVisee+wxjBs3zmmV3XvuuQfJyclISEjAvn37MH/+fBw9ehSffvopAKCkpKTT67U/1tU51dXVaGhoQGBgoNuuKz09HWvWrMGgQYNQXFyMxYsXY8KECThw4ABKSkqg0+k6/GMfGxt70drtj3V1jieur73169ejqqoK9913n+OYt3+G7dlr6qyetvXGxMQ4Pa7RaBAZGel0TmpqaofXsD8WERFxweu2v4anNDY2Yv78+Zg2bZrTZmm/+93vcOWVVyIyMhI7duzAggULUFxcjGXLljmuQ8nXeNNNN+EXv/gFUlNTceLECTzzzDO4+eabkZOTA7Va7XOf47vvvouQkBD84he/cDruLZ9jZ98Tnvp39Pz58z3y/eqX4cSbzZ07FwcOHMC2bducjj/00EOO+8OHD0d8fDxuuOEGnDhxAv369fN0mS67+eabHfdHjBiB9PR0JCcn4+OPP/boF6qnvP3227j55pudthH39s/Q35nNZtx1110QQuDNN990eiwrK8txf8SIEdDpdPjNb36DJUuWeMUS6Hfffbfj/vDhwzFixAj069cPmzdvxg033CBjZe6xevVqTJ8+HQEBAU7HveVzvND3hDfxy26d6OhoqNXqDqOUS0tLERcXJ1NVFzdv3jx88cUX2LRpE/r06dPluenp6QCA48ePAwDi4uI6vV77Y12dExoa6vGAEB4ejoEDB+L48eOIi4uDyWRCVVVVh9ouVrv9sa7O8fT1nTp1Ct988w0efPDBLs/z9s/QXlNXf8/i4uJQVlbm9HhzczMqKyt75LP11N9nezA5deoUvv7664tuMZ+eno7m5mYUFBQA8I5rbKtv376Ijo52+m/TFz5HANi6dSuOHj160b+fgDI/xwt9T3jq39Ge+n71y3Ci0+kwevRoZGdnO45ZrVZkZ2cjIyNDxso6J4TAvHnz8Nlnn+Hbb7/t0GzYmby8PABAfHw8ACAjIwP79+93+gfE/o/okCFDHOe0/TOxnyPHn0ltbS1OnDiB+Ph4jB49Glqt1qm2o0ePorCw0FGbN13fO++8g5iYGEyePLnL87z9M0xNTUVcXJxTPdXV1di1a5fT51ZVVYXc3FzHOd9++y2sVqsjnGVkZGDLli0wm82Oc77++msMGjQIERERjnPkum57MDl27Bi++eYbREVFXfQ5eXl5UKlUjq4QpV9je6dPn0ZFRYXTf5ve/jnavf322xg9ejRGjhx50XOV9Dle7HvCU/+O9tj36yUPnfUxa9euFXq9XqxZs0YcOnRIPPTQQyI8PNxplLJSzJkzR4SFhYnNmzc7TWGrr68XQghx/Phx8cILL4g9e/aI/Px88a9//Uv07dtXXHvttY7XsE8Ru/HGG0VeXp7YuHGj6NWrV6dTxJ566ilx+PBhsWLFCo9NtX3iiSfE5s2bRX5+vti+fbvIzMwU0dHRoqysTAhhmwKXlJQkvv32W7Fnzx6RkZEhMjIyvOb67CwWi0hKShLz5893Ou6tn2FNTY344YcfxA8//CAAiGXLlokffvjBMVNl6dKlIjw8XPzrX/8S+/btE7fddlunU4nT0tLErl27xLZt28SAAQOcpqBWVVWJ2NhYce+994oDBw6ItWvXCoPB0GF6pkajEa+++qo4fPiwWLRoUY9NQe3qGk0mk7j11ltFnz59RF5entPfT/vshh07dog//elPIi8vT5w4cUK8//77olevXmLGjBlecY01NTXiySefFDk5OSI/P19888034sorrxQDBgwQjY2Njtfw5s/Rzmg0CoPBIN58880Oz1f653ix7wkhPPfvaE98v/ptOBFCiNdff10kJSUJnU4nxowZI3bu3Cl3SZ0C0OntnXfeEUIIUVhYKK699loRGRkp9Hq96N+/v3jqqaec1sgQQoiCggJx8803i8DAQBEdHS2eeOIJYTabnc7ZtGmTGDVqlNDpdKJv376O93C3qVOnivj4eKHT6UTv3r3F1KlTxfHjxx2PNzQ0iEceeUREREQIg8Eg7rjjDlFcXOz0Gkq+PruvvvpKABBHjx51Ou6tn+GmTZs6/W9z5syZQgjbdOLnn39exMbGCr1eL2644YYO115RUSGmTZsmgoODRWhoqJg1a5aoqalxOufHH38U48ePF3q9XvTu3VssXbq0Qy0ff/yxGDhwoNDpdGLo0KHiyy+/dPs15ufnX/Dvp339mtzcXJGeni7CwsJEQECAuOKKK8SLL77o9MWu5Gusr68XN954o+jVq5fQarUiOTlZzJ49u8MXjTd/jnZ/+9vfRGBgoKiqqurwfKV/jhf7nhDCs/+OXu73q9RyUURERESK4JdjToiIiEi5GE6IiIhIURhOiIiISFEYToiIiEhRGE6IiIhIURhOiIiISFEYToiIiEhRGE6IiIhIURhOiIiISFEYToiIiEhRGE6IiIhIURhOiIiISFH+P2lT5yrUlV59AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(l,cl_eg_1h,label='1h')\n",
    "plt.plot(l,cl_eg_2h,label='2h')\n",
    "plt.legend()\n",
    "# plt.yscale('log')\n",
    "# plt.yscale('log')\n",
    "plt.ylim(0,1e-10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "6d702e10-a66e-4994-b724-1f202daa30a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\"cleg_l_1h_2h_1p2h.txt\",np.c_[l,cl_eg_1h,cl_eg_2h,cl_eg_1h+cl_eg_2h])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1151ca8-772a-4287-be8e-9bea990cae4a",
   "metadata": {},
   "source": [
    "# shock"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1e65adca-7692-446a-9837-48ab40e5edbe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.16 s, sys: 102 ms, total: 2.26 s\n",
      "Wall time: 450 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "M = Class()\n",
    "\n",
    "M.set(common_params)\n",
    "M.set(cosmo_params)\n",
    "M.set(ksz_params)\n",
    "M.set(HOD_blue)\n",
    "\n",
    "M.set({\n",
    "'use_fft_for_profiles_transform' : 1,\n",
    "'ndim_redshifts':30, ## precision parameter -- how many z in redshift grid.\n",
    "'gas profile mode' : 'shock'\n",
    "})\n",
    "\n",
    "M.compute_class_szfast() # fast mode works\n",
    "# M.compute() # slow mode works"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5f99e8da-b1a3-45e0-a34c-293d9e8c1ae3",
   "metadata": {},
   "outputs": [],
   "source": [
    "l = np.asarray(M.cl_sz()['ell'])\n",
    "cl_eg_1h = np.asarray(M.cl_eg()['1h'])/l_to_dl(l)\n",
    "cl_eg_2h = np.asarray(M.cl_eg()['2h'])/l_to_dl(l)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "910c7a54-f1fa-4c1f-b6e1-75cc349fcd84",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.00000000e+00, 2.44280552e+00, 2.98364940e+00, 3.64423760e+00,\n",
       "       4.45108186e+00, 5.43656366e+00, 6.64023385e+00, 8.11039993e+00,\n",
       "       9.90606485e+00, 1.20992949e+01, 1.47781122e+01, 1.80500270e+01,\n",
       "       2.20463528e+01, 2.69274761e+01, 3.28892935e+01, 4.01710738e+01,\n",
       "       4.90650604e+01, 5.99282001e+01, 7.31964689e+01, 8.94023690e+01,\n",
       "       1.09196300e+02, 1.33372662e+02, 1.62901737e+02, 1.98968631e+02,\n",
       "       2.43020835e+02, 2.96826318e+02, 3.62544484e+02, 4.42812832e+02,\n",
       "       5.40852815e+02, 6.60599120e+02, 8.06857587e+02, 9.85498082e+02,\n",
       "       1.20369008e+03, 1.47019038e+03, 1.79569458e+03, 2.19326632e+03,\n",
       "       2.67886153e+03, 3.27196886e+03, 3.99639179e+03, 4.88120396e+03,\n",
       "       5.96191597e+03, 7.28190061e+03, 8.89413350e+03, 1.08633192e+04,\n",
       "       1.32684880e+04, 1.62061679e+04, 1.97942581e+04])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2bab2b53-7c03-46f8-b703-820811df069d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 1e-10)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGsCAYAAAAGzwdbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABRBklEQVR4nO3deXzT9f0H8Nc3d++D3hflPuQuUCqgODsQEe+JiIJMcTLY1M4NUYR5/ARvN0VxTsTNKSgqboIoVhCVAlKo3CjQ0gI9KKXp3bTJ5/dHmrShBRpI8v0meT0fjz6afvNN8v420rz8nJIQQoCIiIhIIVRyF0BERETUFsMJERERKQrDCRERESkKwwkREREpCsMJERERKQrDCRERESkKwwkREREpCsMJERERKQrDCRERESkKwwkREREpileFk82bN2Py5MlISEiAJElYs2aN7K8nhMDChQsRHx+PgIAAZGZm4pdffnFrXURERL7Mq8JJbW0tBg8ejKVLlyrm9Z577jn8/e9/x7Jly7Bt2zYEBQVhwoQJaGho8EiNREREvkby1o3/JEnCp59+ihtvvNF+rLGxEY899hg++OADVFZWYsCAAXj22Wcxbtw4t7yeEAIJCQn405/+hIcffhgAYDQaERsbixUrVuD222+/5NclIiLyN17VcnIhc+fORU5ODlauXIndu3fjN7/5Da655hq3dbPk5+ejpKQEmZmZ9mNhYWFIT09HTk6OW16TiIjI1/lMOCksLMQ777yDjz76CGPHjkWPHj3w8MMPY8yYMXjnnXfc8polJSUAgNjYWIfjsbGx9vuIiIjIOT4TTvbs2QOz2YzevXsjODjY/vXtt9/iyJEjAICDBw9CkqTzfj3yyCMyXwkREZF/08hdgKvU1NRArVYjNzcXarXa4b7g4GAAQPfu3XHgwIHzPk+XLl06/ZpxcXEAgNLSUsTHx9uPl5aWYsiQIZ1+HiIiImrlM+Fk6NChMJvNKCsrw9ixYzs8R6fToW/fvi57zW7duiEuLg7Z2dn2MFJVVYVt27Zh9uzZLnsdIiIif+JV4aSmpgaHDx+2/5yfn4+8vDxERkaid+/emDZtGqZPn44XX3wRQ4cOxalTp5CdnY1BgwZh0qRJLn29lJQUSJKEBx98EE8//TR69eqFbt264fHHH0dCQoLDrB4iIiLqPK+aSrxp0yZcddVV7Y7PmDEDK1asQFNTE55++mn861//wokTJxAVFYVRo0bhiSeewMCBA13+eoB1OvGiRYvwj3/8A5WVlRgzZgxef/119O7d2+nXIyIioosIJ5s3b8bzzz+P3NxcFBcXt1v7oyObNm1CVlYW9u3bh+TkZCxYsAB33333JZRNREREvsrp2TrOrtKan5+PSZMm4aqrrkJeXh4efPBB3Hvvvfjyyy+dLpaIiIh83yV163S0aurZ5s2bh7Vr12Lv3r32Y7fffjsqKyuxfv36i31pIiIi8lFuHxCbk5PjsIIqAEyYMAEPPvjgOR/T2NiIxsZG+88WiwUVFRXo0qULJElyV6lERETkQkIIVFdXIyEhASpV5ztr3B5OSkpKOlxBtaqqCvX19QgICGj3mMWLF+OJJ55wd2lERETkAUVFRUhKSur0+YqcSjx//nxkZWXZfzYajUhJSUFRURFCQ0Nd9jqz3t2BnKOn8ewtAzFpUILLnhcA8OZVQMUvwNRVQOpo1z63E/6dU4Bn1x/Cr/vH4OUpQ2Wrg4iI/E9VVRWSk5MREhLi1OPcHk7i4uJQWlrqcKy0tBShoaEdtpoAgF6vh16vb3c8NDTUpeFEGxAElb4egcEhLn1eAEBoEFArAQY14OrndkLflFio9EUorlO5/hqJiIg6wdkhGW7fWycjIwPZ2dkOxzZs2ICMjAx3v7S8NC3Bq7le1jK6dgkEABRW1MGLlrQhIiI/5nQ4qampQV5eHvLy8gC0rppaWFgIwNolM336dPv5999/P44ePYq//OUvOHjwIF5//XV8+OGHeOihh1xzBUqlNVi/NzXIWkZSRCAkCahpbEZFrUnWWoiIiDrD6XCyY8cODB06FEOHWscvZGVlYejQoVi4cCEAoLi42B5UAOv+M2vXrsWGDRswePBgvPjii/jnP/+JCRMmuOgSFEprbbGQu+XEoFUjPtQalApO18laCxERUWc4PeZk3Lhx5+0esC3rfvZjdu3a5exLeTeNreVE3nACACldAnHS2IDCilqkdY2QuxwiIp8ghEBzczPMZrPcpchGrVZDo9G4fJkPRc7W8QnaljEnCggnXSODsPVoBY6x5YSIyCVMJhOKi4tRV8e/q4GBgYiPj4dOp3PZczKcuIut5aRZ3jEngLXlBAAKGU6IiC6ZxWJBfn4+1Go1EhISoNPp/HKBUCEETCYTTp06hfz8fPTq1cuphdbOx6/DiYAbZ68oqOUktUsQAKDgdK3MlRAReT+TyQSLxYLk5GQEBgbKXY6sAgICoNVqcezYMZhMJhgMBpc8r9unEvstWzhRQMtJ2+nERETkGq5qJfB27vg98DfrLvYBsfIHAlu3TnmNCTWNzTJXQ0REdH4MJ+5i79aRv+Uk1KBFRKAWAMedEBGR8jGcuIuCBsQCQNeWcSfHOO6EiIjOITU1Fa+88orcZTCcuI2CBsQCreNOjnHcCRGR39q8eTMmT56MhIQESJKENWvWyF1ShxhO3EVBA2IBoGtkSzhhtw4Rkd+qra3F4MGDsXTpUrlLOS+/nkrsVhpltZyktHTrFFawW4eIyNWEEKhv8vxKsQFatVNrrEycOBETJ0487zl1dXX47W9/i48++ggRERFYsGAB7rvvvkst1SkMJ+6iVc7y9UBrt05BOVtOiIhcrb7JjP4Lv/T46+5/cgICda79KH/xxRfx1FNP4dFHH8Xq1asxe/ZsXHnllejTp49LX+d8/Lpb5zxbBF06W8uJzBv/2djCSbGxHqZmi8zVEBGRUl177bX4/e9/j549e2LevHmIiorCxo0bPVoDW07cxd5yoowxJ9HBegTq1KgzmXH8TB26RwfLXRIRkc8I0Kqx/8kJsryuqw0aNMh+W5IkxMXFoayszOWvcz4MJ+6ibVnSWCEtJ5IkISUyEAdLqnGsguGEiMiVJElyefeKXLRarcPPkiTBYvFsi7tfd+u4lUZZLScAkGKbsVPOQbFERKRcvhHzlMg2ldjcCFjMgMr1TW/O4lonRET+raamBocPH7b/nJ+fj7y8PERGRiIlJUXGyhwxnLiLps3OjM0NgC5Ivlpa2FaJ5RL2RET+aceOHbjqqqvsP2dlZQEAZsyYgRUrVshUVXsMJ+5iazkBrF07iggnbDkhIvJn48aNgzjPVNWCgoJ2x/Ly8txX0DlwzIm7qNSAqmVQkUIGxXaNtC3EVgeLxZ3zqImIiC4ew4k72WbsKGRQbEK4ARqVBFOzBSVVyqiJiIjobAwn7mRb60QhLScatQpJEdbupgLuTkxERArl1+HErSvEAm2mEysjnABAzxjr+iY/l1TLXAkREVHH/Dqc2DizaZJTtMra/A8A+sWHAgAOFDOcEBGRMjGcuJOt5aRZOeM7+reEk/3FVTJXQkRE1DGGE3eyD4hVTstJ/wRrODlUWo1mMzcAJCIi5WE4cSet8lpOkiMCEaRTw9RswVEuY09ERArEcOJOGtuYE+UseqZSSfZxJ/tPsmuHiIiUh+HEnbTK2/wPaB0Uy3EnRER0NkmSsGbNGllrYDhxJ1vLiULWObGxjTs5wHBCRORXFi9ejBEjRiAkJAQxMTG48cYbcejQIbnLaofhxJ3sU4mV1XLSv023zvn2WCAiIt/y7bffYs6cOdi6dSs2bNiApqYmjB8/HrW1yhqDyI3/3ElhK8Ta9IkLgUoCTteaUFbdiNhQw4UfRERE5yaEPOMLtYGAE2t1rV+/3uHnFStWICYmBrm5ubjiiivsx8vLy3HTTTfhyy+/RGJiIl588UVcf/31Liv7QhhO3EmjzJYTg1aN7tHBOFxWg/3FVQwnRESXqqkOeCbB86/76MlL2vXeaDQCACIjIx2OP/HEE3juuefw/PPP49VXX8W0adNw7Nixdue5C7t13EmrvOXrbfpzxg4RkV+zWCx48MEHMXr0aAwYMMDhvrvvvhtTp05Fz5498cwzz6Cmpgbbt2/3WG1sOXEnhQ6IBawzdv7700nO2CEicgVtoLUVQ47XvUhz5szB3r178f3337e7b9CgQfbbQUFBCA0NRVlZ2UW/lrMYTtxJoQNiAc7YISJyKUm6pO4VT5s7dy4+//xzbN68GUlJSe3u12q1Dj9LkgSLxXOrirNbB4Cbtv1rDScKbDmxdevkl9eiztQsczVEROQJQgjMnTsXn376Kb755ht069ZN7pI6xHDiThplLsIGANEhekQF6yEEcKiEOxQTEfmDOXPm4L333sP777+PkJAQlJSUoKSkBPX1yvqfaIYTd9Iqb/n6tmxdOxx3QkTkH9544w0YjUaMGzcO8fHx9q9Vq1bJXZoDvx5zIuDmBcg0ytv4r63+8aHY/PMpztghIvITnVl4s6NzKisr3VDNubHlxJ3sLSfKai6z6RcfAoCDYomISFkYTtzJPiBWmS0nl7V06xwsqYbZwmXsiYhIGRhO3EmhK8TadIsKhkGrQp3JjGOnlbWvAhER+S+GE3dS6N46NmqVhD6xtq4dztghIiJlYDhxJ1vLidkEWMzy1nIOrTN2jDJXQkTkXbiru5U7fg8MJ+6kbbOhnkIHxXKPHSIi59hWT62rU+YyEZ5m+z2cvarspfDrqcRuZ2s5AayDYvXB8tVyDv3ibcvYs1uHiKgz1Go1wsPD7XvNBAYGQpLctta4YgkhUFdXh7KyMoSHh0OtVrvsuRlO3EmlAtR6wNyo2JaTvi3hpKSqAadrGtElWC9zRUREyhcXFwcAHt0MT6nCw8Ptvw9XYThxN63BGk4UOp04WK9B1y6BOHa6DgeKqzGmF8MJEdGFSJKE+Ph4xMTEoKmpSe5yZKPVal3aYmLDcALrZpJuowkAYFTsEvaAddyJNZxUYUyvKLnLISLyGmq12i0fzv6OA2LdTavczf9s7INiuVIsEREpAMOJu2kDrd8VutYJ0DooljN2iIhICRhO3E3jBS0nLWudHDlVg4YmZa7HQkRE/oPhxN3s++sot+UkPsyA8EAtmi0Ch8tq5C6HiIj8HMOJu9lbTpQbTiRJQr84du0QEZEyMJy4m63lRMHhBGi7jD3DCRERyYvhxN3s3TrKHXMCcMYOEREph1+HE4/s2eQF3TpA22Xsq7iZFRERycqvw4lHeEnLSc+YYGjVEqobmnH8jLKDFBER+TaGE3fzkpYTnUaFXjEhANi1Q0RE8mI4cTcvGRALOHbtEBERyYXhxN28pFsHaDNjh9OJiYhIRgwnACS4cec/jfe0nHDGDhERKQHDibvZNv7zhpaTlnBy/Ew9jPX+uwU4ERHJ66LCydKlS5GamgqDwYD09HRs3779vOe/8sor6NOnDwICApCcnIyHHnoIDQ3K/7B2CS9qOQkL1CIx3FrvQbaeEBGRTJwOJ6tWrUJWVhYWLVqEnTt3YvDgwZgwYQLKyso6PP/999/HI488gkWLFuHAgQN4++23sWrVKjz66KOXXLxX0HrHbB2bfuzaISIimTkdTl566SXMmjULM2fORP/+/bFs2TIEBgZi+fLlHZ6/ZcsWjB49GnfccQdSU1Mxfvx4TJ069YKtLT5Do/yN/9rqH2+dTswZO0REJBenwonJZEJubi4yMzNbn0ClQmZmJnJycjp8zOWXX47c3Fx7GDl69CjWrVuHa6+99pyv09jYiKqqKocvr2WfSuwd3VjcY4eIiOSmcebk8vJymM1mxMbGOhyPjY3FwYMHO3zMHXfcgfLycowZMwZCCDQ3N+P+++8/b7fO4sWL8cQTTzhTmnJpva3lJAwA8HNJDZrMFmjVHDNNRESe5fZPnk2bNuGZZ57B66+/jp07d+KTTz7B2rVr8dRTT53zMfPnz4fRaLR/FRUVubtM97GvEOsdLSdJEQEI1mtgMltw9FSt3OUQEZEfcqrlJCoqCmq1GqWlpQ7HS0tLERcX1+FjHn/8cdx111249957AQADBw5EbW0t7rvvPjz22GNQqdrnI71eD71e70xpyuVFi7ABgEoloV98CH4sOIP9xUb0iQuRuyQiIvIzTrWc6HQ6pKWlITs7237MYrEgOzsbGRkZHT6mrq6uXQBRq9UA4B+739pbTurkrcMJ9hk7XCmWiIhk4FTLCQBkZWVhxowZGD58OEaOHIlXXnkFtbW1mDlzJgBg+vTpSExMxOLFiwEAkydPxksvvYShQ4ciPT0dhw8fxuOPP47JkyfbQ4pPs7WcWJoBczOgdvpX7nH97XvsVMtcCRER+SOnPymnTJmCU6dOYeHChSgpKcGQIUOwfv16+yDZwsJCh5aSBQsWQJIkLFiwACdOnEB0dDQmT56M//u//3PdVSiZLZwA1kGxauV3k7SdsSOEgCS5cXl/IiKis1zU/8bPnTsXc+fO7fC+TZs2Ob6ARoNFixZh0aJFF/NS3s/WrQNYB8XqlR9OeseGQCUBFbUmlFY1Ii7McOEHERERuYhfzxP1yIgXSWoNKF4yndigVaNHdDAALsZGRESe59fhxMbtvRYa71rCHuBibEREJB+GE0/QBlq/e1E44YwdIiKSC8OJJ9g2//OStU6AtjN2GE6IiMizGE48wbb5nxe2nOSfrkWdqVnmaoiIyJ8wnHiCF7acRIfoER2ihxDAwRKud0JERJ7DcOIJXthyArR27XDcCREReRLDiSdovTSccMYOERHJgOHEE7Tetc6JTT8OiiUiIhkwnHiCvVvHe8acAK3dOgeLq2G2+MEmjUREpAgMJ57gpS0n3aKCYNCqUN9kRsHpWrnLISIiP+Hf4cRTjQFe2nKiVknoE8euHSIi8iz/DieeovW+5ettOGOHiIg8jeHEE2zL13tZtw7AGTtEROR5DCeeYN/4z7u6dQCgf3wIAHbrEBGR5zCcAHD3psT2dU68sOWkT1woJAkorWpEeU2j3OUQEZEfYDjxBC9uOQnWa5DaJQgAW0+IiMgzGE48wYtbTgCgH7t2iIjIgxhOPMFLl6+34YwdIiLyJIYTT/DSjf9sbDN2fjpulLkSIiLyBwwnnmBfIdb7xpwAQFrXSKgkIL+8FsVG7wxYRETkPRhOPMFLV4i1CQvQYmBSOADgh8On5S2GiIh8nl+HE+Gp9eu9dG+dtkb36AIA+OFwucyVEBGRr/PrcOIxXj4gFgDG9IwCYA0nQnCHYiIich+GE0/w8gGxADCsawT0GhXKqhtx5FSN3OUQEZEPYzjxBFu3jjAD5iZ5a7lIBq0aw1MjAADf/8KuHSIich+GE0+wtZwAXt16MtrWtXOEg2KJiMh9GE48QaOHfQcfL51ODACje1jDydajp9FstshcDRER+SqGE0+QJJ8YFDsgMQyhBg2qG5qx5wQXZCMiIvdgOPEU++Z/3htO1CoJGS1Tirewa4eIiNyE4QTWhg238/LN/2xGt5lSTERE5A4MJ55ibznx3jEnAHB5y7iTHcfOoKHJLHM1RETkixhOPMVHWk56RAchLtQAU7MFOwrOyF0OERH5IL8OJx5d6FTr3fvr2EiShMt7tixlf4RdO0RE5Hp+HU48yt6tUydvHS5gm1LMcSdEROQODCeeYu/W8e6WE6B1UOyeE0YY67xzxVsiIlIuhhNP8YGpxDZxYQb0iA6CEEDOUU4pJiIi12I48RQfajkBOKWYiIjch+HEU3xkKrFN6z47DCdERORaDCeeog20fveBAbEAMKp7F6gk4OipWhQbvb+rioiIlIPhxFO0LS0nPtKtExagxcDEMADAD4c57oSIiFyH4cRTNN6/8d/ZbF07WzjuhIiIXIjhxFN8rOUEcBx3Ijy6oh0REfkyhhNP8cGWk7SuEdBpVCitasSRUzVyl0NERD7Cr8OJR/9fXx9s/d5Y7clXdSuDVo0RqREAOO6EiIhcx6/DSSvJ/S9hCLd+b6h0/2t5kG2X4u857oSIiFyE4cRTAsKt3+sr5azC5WzjTrYePY1ms0XmaoiIyBcwnHiKj7acDEwMQ4hBg+qGZuw9WSV3OURE5AMYTjylbcuJD81sUaskZHTvAoBL2RMRkWswnHhKgHXgKIQZMPnWzBbus0NERK7EcOIp2gBArbferj8jby0uZgsnO46dQUOTWeZqiIjI2zGceJKPDortER2E2FA9TM0W7CjwreBFRESex3DiST46KFaSJIzuwV2KiYjINRhOPMlHW04A7rNDRESuw3DiSbZBsT7WcgK0hpPdJ4ww1jXJXA0REXkzhhNPsnXr+NiAWACICzOge3QQhAByjnIpeyIiungMJ57kw906ADDG1rXDcSdERHQJGE48yUcHxNpwnx0iInIFhhMAkgf2/QPg8y0nGd27QCUBR0/VosTYIHc5RETkpRhOPMmHB8QCQFigFgMTwwBwtVgiIrp4DCee5MMDYm0u51L2RER0iRhOPMnHu3WA1kGxPxwph/ChDQ6JiMhzGE48yccHxAJAWtcI6DQqlFY14sipWrnLISIiL3RR4WTp0qVITU2FwWBAeno6tm/fft7zKysrMWfOHMTHx0Ov16N3795Yt27dRRXsSh7/P3tby0mDEbBYPPvaHmLQqjG8q3VsDbt2iIjoYjgdTlatWoWsrCwsWrQIO3fuxODBgzFhwgSUlZV1eL7JZMKvf/1rFBQUYPXq1Th06BDeeustJCYmXnLxXsfWciIsgKla1lLcaTTHnRAR0SVwOpy89NJLmDVrFmbOnIn+/ftj2bJlCAwMxPLlyzs8f/ny5aioqMCaNWswevRopKam4sorr8TgwYMvuXivozUAmgDrbR8ed2ILJzlHT8Ns4bgTIiJyjlPhxGQyITc3F5mZma1PoFIhMzMTOTk5HT7mv//9LzIyMjBnzhzExsZiwIABeOaZZ2A2m8/5Oo2NjaiqqnL48hn2QbG+O2NnYGIYQgwaVDc0Y88Jo9zlEBGRl3EqnJSXl8NsNiM2NtbheGxsLEpKSjp8zNGjR7F69WqYzWasW7cOjz/+OF588UU8/fTT53ydxYsXIywszP6VnJzsTJnK5geDYtUqCRnduwBg1w4RETnP7bN1LBYLYmJi8I9//ANpaWmYMmUKHnvsMSxbtuycj5k/fz6MRqP9q6ioyN1leo4fTCcGWrt2uM8OERE5S+PMyVFRUVCr1SgtLXU4Xlpairi4uA4fEx8fD61WC7VabT/Wr18/lJSUwGQyQafTtXuMXq+HXq93pjTv4QctJwAwuqe15eTHgjNoaDLDoFVf4BFERERWTrWc6HQ6pKWlITs7237MYrEgOzsbGRkZHT5m9OjROHz4MCxtps7+/PPPiI+P7zCY+DzbEvY+3nLSIzoYsaF6mJotyD3mu+NriIjI9Zzu1snKysJbb72Fd999FwcOHMDs2bNRW1uLmTNnAgCmT5+O+fPn28+fPXs2Kioq8MADD+Dnn3/G2rVr8cwzz2DOnDmuu4pL5Kl9/wD4xYBYAJAkCaN7cEoxERE5z6luHQCYMmUKTp06hYULF6KkpARDhgzB+vXr7YNkCwsLoVK1Zp7k5GR8+eWXeOihhzBo0CAkJibigQcewLx581x3Fd7ET7p1AOs+O5/sOsFwQkRETnE6nADA3LlzMXfu3A7v27RpU7tjGRkZ2Lp168W8lO/xkwGxQOu4kz0njDDWNyEsQCtzRURE5A24t46n+VHLSXxYALpHB8EigK1HT8tdDhEReQmGE0/zkwGxNhx3QkREzmI48TQ/GRBrw312iIjIWQwnnuZH3ToAkNG9C1QScORULUqMDXKXQ0REXoDhxNNsLScNVUCbtV98VVigFgMSwwCw9YSIiDqH4cTTbC0nEECjf2yKZ+/a4VL2RETUCQwnnqbRAdog620/GxS75fBpCCFkroaIiJSO4UQOfjYodnhqBHQaFUqqGnDkVK3c5RARkcL5dTgZkBiGUd0jERnk4T1+/GxQrEGrRnq3SADA2t3FMldDRERKd1ErxPqKJ28YIM8L+9EqsTY3D0vEd7+U48MdRfjDr3pCpfLojkZERORF/LrlRDZ+1nICABMHxCPEoMGJynoOjCUiovNiOJGDH7acGLRq3DgkEQCw6scimashIiIlYziRg30Je/8YEGszZUQyAOCrfaU4U2uSuRoiIlIqhhM5+GG3DmAdgHxZQihMZgs+3XVC7nKIiEihGE7k4IfdOja21pNVPxZxzRMiIuoQw4kc/LTlBABuGJwInUaFQ6XV+Om4f6yQS0REzmE4kYMft5yEBWpx7YA4ABwYS0REHWM4kYN9QGylrGXI5baWrp3//XQSdaZmmashIiKlYTiRgx936wDAqG5d0LVLIGoam7liLBERtcNwIgdbt05jFWAxy1qKHFQqCbcNt7aefLiDXTtEROSI4UQOhrDW2w3+OSj0lmFJUEnAjwVncORUjdzlEBGRgjCcyEGtBXTB1tt+thCbTVyYAVf1iQEAfMiBsURE1AbDiVz8fFAs0Dow9uOdx9FktshcDRERKQXDiVzsg2L9s+UEAH7VNwZRwXqU15iQfaBM7nKIiEghGE7k4sdrndho1SrckmbdDJADY4mIyIbhRC62QbF+Op3YxjZrZ9OhMpQYG2SuhoiIlIDhRC5sOQEA9IgOxojUCFgEsDqXrSdERMRwIh/bmBM/na3T1pQRKQCAD3cch8XCzQCJiPwdw4lcbLN1/LxbBwCuHRiHYL0GhRV12Jp/Wu5yiIhIZgwncmG3jl2gToPJgxMAcDNAIiJiOJGPfSqxf64Qe7bbW9Y8+WJvCYx1TTJXQ0REcmI4kQtbThwMSgpD37gQmJot+OynE3KXQ0REMmI4kYvBtkIsB8QCgCS1bga4cju7doiI/BnDiVxsLSccEGt309BE6NQq7C+uwt4T7O4iIvJXDCdysc3WMdUAZo6xAICIIB3GXxYLgANjiYj8GcOJXGwrxAIcFNvGlJaBsWvyTqChySxzNUREJAeGE7mo1IA+1Hqbg2LtRveIQmJ4AKobmvHF3mK5yyEiIhkwnMiJq8S2o1K1Doxl1w4RkX9iOJETB8V26NbhSZAkYOvRChSU18pdDhEReRjDiZy41kmHEsMDcEWvaADAhzvYekJE5G8YTuRkXyW2Us4qFMk2MHZ17nE0my0yV0NERJ7EcCIntpycU2a/WEQG6VBW3YhNh07JXQ4REXkQw4mcOCD2nHQaFW4emggAWMWuHSIiv8JwIifbQmzs1umQrWvnm4NlKKtqkLkaIiLyFIYTObFb57x6xYZgaEo4zBaBj3dyM0AiIn/BcCInDoi9oNtbWk8+3FEEIYTM1RARkScwnMiJLScXNGlQAgJ1auSX12J7foXc5RARkQcwnMiJLScXFKzX4LpB8QA4MJaIyF8wnMjJ3nLC2TrnM2VECgBg3Z5iVDVwB2ciIl/HcCIn22ydpjqg2SRvLQo2LCUcPWOC0dBkwX/zTspdDhERuRnDiZz0YQAk62127ZyTJEmYMrx1YCwREfk2hhM5qVSAIdR6m4Niz+umYYnQqiXsPm7E/pNVcpdDRERuxHAiNw6K7ZSoYD0y+8UCYOsJEZGvYziRGwfFdpptxdhPd51AQ5NZ5mqIiMhdGE7kZhsUy26dCxrbKxoJYQYY65vw1f5SucshIiI3YTiRG7t1Ok2tknBrWhIAYNWPhTJXQ0RE7sJwIjeuEuuU3wxPhiQBPxw+jaKKOrnLISIiN2A4kRtbTpySHBmI0T2iAHBgLBGRr2I4kRsHxDrttpaBsatzj8Ns4WaARES+huFEbhwQ67Tx/WMRHqhFsbEBmw6VyV0OERG5GMOJ3GzhpO60vHV4EYNWjVuHWQfGPrf+EJrMFpkrIiIiV2I4kVuY9UMWRo6fcMacq3oiPFCLQ6XV+FfOMbnLISIiF7qocLJ06VKkpqbCYDAgPT0d27dv79TjVq5cCUmScOONN17My/qm8FTr9+pioKlB1lK8SUSQDn+Z0BcA8MqGn1FWzd8dEZGvcDqcrFq1CllZWVi0aBF27tyJwYMHY8KECSgrO3/ff0FBAR5++GGMHTv2oov1SYGRgDbIeputJ06ZMiIZg5LCUN3YjCXrDspdDhERuYjT4eSll17CrFmzMHPmTPTv3x/Lli1DYGAgli9ffs7HmM1mTJs2DU888QS6d+9+SQX7HEkCIrpab1eye8IZapWEJ28YAEkCPtl1Aj8WVMhdEhERuYBT4cRkMiE3NxeZmZmtT6BSITMzEzk5Oed83JNPPomYmBjcc889nXqdxsZGVFVVOXz5tPCWcHKG4cRZQ5LDMWW4dWrxws/2oZmDY4mIvJ5T4aS8vBxmsxmxsbEOx2NjY1FSUtLhY77//nu8/fbbeOuttzr9OosXL0ZYWJj9Kzk52ZkyvU94ivV7JZdkvxh/uaYvwgK0OFBchf9s4++QiMjbuXW2TnV1Ne666y689dZbiIqK6vTj5s+fD6PRaP8qKvLxsRjs1rkkkUE6PDy+NwDgha8OobymUeaKiIjoUmicOTkqKgpqtRqlpY47wpaWliIuLq7d+UeOHEFBQQEmT55sP2axWJvdNRoNDh06hB49erR7nF6vh16vd6Y072ZrOWG3zkW7I70rVv5YhH0nq/DsFwfx/G8Gy10SERFdJKdaTnQ6HdLS0pCdnW0/ZrFYkJ2djYyMjHbn9+3bF3v27EFeXp796/rrr8dVV12FvLw83++u6SzbmBN261w02+BYAPgo9zhyj3E7ACIib+VUywkAZGVlYcaMGRg+fDhGjhyJV155BbW1tZg5cyYAYPr06UhMTMTixYthMBgwYMAAh8eHh4cDQLvjfs3WclJXDjTWAPpgeevxUmldI3BrWhJW5x7Hov/uxWdzxkCtkuQui4iInOR0OJkyZQpOnTqFhQsXoqSkBEOGDMH69evtg2QLCwuhUnHhWacEhAOGMKDBaF3rJKaf3BV5rUcm9sWX+0qw90QVPtheiDtHdZW7JCIicpIkhFD8tq5VVVUICwuD0WhEaGio3OW4x7KxQMluYOoqoM81clfj1Vb8kI+//m8/wgK02PjwOEQG6eQuiYjIL13s5zebOJSC04ld5s5RXdE3LgTG+iY8/yVXjiUi8jYMJ0oRkWr9zunEl0yjVtkHx678sQg/FVXKWxARETmF4UQp7NOJC2Qtw1eM7BaJm4YmQghg4Wd7YbEovveSiIhaMJwoBacTu9z8iX0RrNfgp+NGfLjDxxfyIyLyIQwnSsFVYl0uJtSABzN7AQCeXX8QlXUmmSsiIqLOYDhRirCWBekajEB9payl+JIZl6eid2wwztQ14YWvDsldDhERdQLDiVLog4HAlv2H2LXjMlq1Ck9cbx0c+59thdh7wihzRUREdCEMJ0rC6cRukdGjCyYPToAQwOMcHEtEpHgMJ0rCcSdu89i1/RCkU2NXYSU+3nlc7nKIiOg8GE6UhLsTu01cmAF/vNo6OHbJFwdhrG+SuSIiIjoXhhMl4XRit5o5uht6RAfhdK0JL2/4We5yiIjoHBhOlITdOm6l07QOjv1XTgH2n6ySuSIiIuoIw4mStG05Uf5+jF5pTK8oXDswDhYBLPrvXnjBvpdERH6H4URJbGudmGqAugp5a/FhCyb1R4BWjR8LzuDTXSfkLoeIiM7CcKIkWgMQHGe9XVkgaym+LCE8AHN/1RMA8My6g6hu4OBYIiIlYThRmggOivWEe8d2Q7eoIJTXNOKVr3+RuxwiImqD4URpbONOOJ3YrfQaNf56/WUAgBVbCnCopFrmioiIyIbhRGm4SqzHXNk7GhMui4XZIrDwMw6OJSJSCoYTpeF0Yo9aMKk/9BoVtuVX4L8/nZS7HCIiAsOJ8rDlxKOSIwMx5yrb4NgDqGlslrkiIiJiOFEarnXicfdd0R1duwSitKoRr2ZzcCwRkdwYTpQmLAmQVEBzA1BTKnc1fsGgVWPR5P4AgLe/z8fhMg6OJSKSE8OJ0qi1QGii9Ta7djzmV31jkdkvBs0WgUX/3cfBsUREMmI4USJOJ5bFwusug06jwg+HT2PdnhK5yyEi8lsMJ0pkHxTLcOJJKV0Ccf+VPQAAT6/dj1oOjiUikgXDiRJxOrFsfj+uB5IiAlBsbMDSjYflLoeIyC8xnCgRpxPLxqBVY+F11sGxb313FEdP1chcERGR/2E4USKOOZHVr/vHYlyfaDSZrYNjLRYOjiUi8iSGEyWytZwYjwMWs7y1+CFJkvDXyZdBp1bhu1/K8cT/OHuHiMiTGE6UKDQBUGkBSxNQXSx3NX4pNSoIz946EJIEvJtzDM99eUjukoiI/AbDiRKp1NbF2AB27cjopqFJePrGAQCANzYd4QBZIiIPYThRKg6KVYRp6V2xYFI/AMDzXx7C8u/zZa6IiMj3MZwoFacTK8a9Y7vjwcxeAIAnP9+PldsZGImI3InhRKlsLSfs1lGEB67uhfuu6A4AmP/pHnyWd0LmioiIfBfDiVKFp1q/s1tHESRJwvyJfTEtPQVCAFkf/oSv9nGJeyIid2A4USp26yiOJEl46oYBuHloIswWgbnv78Lmn0/JXRYRkc9hOFEqW7dO1QnA3CRvLWSnUkl47tZBmDggDiazBff9ewe251fIXRYRkU9hOFGq4FhAYwCExboYGymGRq3C324finF9otHQZMFvV/yIn4oq5S6LiMhnMJwolSQBYcnW2xx3ojg6jQrL7kzDqO6RqGlsxox3tuNgSZXcZRER+QSGEyXjuBNFM2jV+OeMERiSHI7Kuibc+c/t3CiQiMgFGE6UjNOJFS9Yr8G7M0eiX3woymsaMe2f21BUUSd3WUREXo3hRMlsuxOzW0fRwgK1+Pc9I9EjOgjFxgbc+fY2lFY1yF0WEZHXYjhRMnbreI2oYD3+c+8oJEcG4NjpOtz5z204XdMod1lERF6J4UTJuL+OV4kLM+D9e0chLtSAX8pqMH35dhjrOQ2ciMhZDCdKZlsltroYaGI3gTdIjgzEe/emo0uQDvtOVmHmO9tR29gsd1lERF6F4UTJAiMBXbD19pkCWUuhzusZE4x/35OOUIMGOwsrMetfO9DQZJa7LCIir8FwomSSBMQPtt4u2iZvLeSU/gmhePe3IxGkU2PLkdP4/X92wtRskbssIiKvwHCidF1HW78f+0HeOshpQ1Mi8PbdI6DXqPDNwTI8tCoPZouQuywiIsVjOFG6rpdbvxf8AAh+sHmbUd274M270qBVS1i7pxjzPt4NCwMKEdF5MZwoXfJIQKUBqo5z1o6XGtcnBq9OHQa1SsLq3OP46//2QTBoEhGdE8OJ0umCgIRh1tvs2vFa1wyIwwu/GQRJAv6VcwxL1h9kQCEiOgeGE2/QtmuHvNZNQ5Pw9I0DAABvfnsUr31zWOaKiIiUieHEG6SOsX5ny4nXm5beFQsm9QMAvLjhZ7z9fb7MFRERKQ/DiTdITgckFXAmH6g6KXc1dInuHdsdD2X2BgA89fl+fLCdY4mIiNpiOPEGhlAgbpD1Nrt2fMIfr+6J313RHQDw6Kd78FneCZkrIiJSDoYTb2Hv2vle3jrIJSRJwiMT++LOUSkQAsj68Cd8ua9E7rKIiBSB4cRb2Bdj2yJvHeQykiThyesH4OZhiTBbBP7w/i58+/MpucsiIpIdw4m36JoBQALKfwZqyuSuhlxEpZLw3C2DMHFAHExmC3737x3YdvS03GUREcmK4cRbBEQAsZdZb3PWjk/RqFX42+1DMa5PNBqaLLjn3R34qahS7rKIiGTDcOJNbF07HBTrc3QaFZbdmYZR3SNR09iMO9/ehlU/FnKhNiLySwwn3iSV4058mUGrxj9njMCI1AhUNzRj3sd7cPs/tuLIqRq5SyMi8iiGE2+S0rJSbNk+oK5C3lrILYL1GnwwaxQeu7YfArRqbMuvwMRXvsPfvv4Fjc1mucsjIvIIhhNvEhwNRPWx3mbric/SqFWYdUV3fPXQFbiydzRMZgte/vpnXPu377A9n6GUiHzfRYWTpUuXIjU1FQaDAenp6di+ffs5z33rrbcwduxYREREICIiApmZmec9ny6AXTt+IzkyECtmjsDfpw5FVLAOR07V4rY3czD/kz0w1jXJXR4Rkds4HU5WrVqFrKwsLFq0CDt37sTgwYMxYcIElJV1PL1106ZNmDp1KjZu3IicnBwkJydj/PjxOHGCK2JeFPt6J1yMzR9IkoTrByfg66wrMWV4MgDgg+2FuPqlb/H57pMcMEtEPkkSTv51S09Px4gRI/Daa68BACwWC5KTk/GHP/wBjzzyyAUfbzabERERgddeew3Tp0/v1GtWVVUhLCwMRqMRoaGhzpTre6qKgZf6WvfamVcAGMLkrog8aOvR03j00z04eqoWAPCrvjF48obLkBQRKHNlRETtXeznt1MtJyaTCbm5ucjMzGx9ApUKmZmZyMnJ6dRz1NXVoampCZGRkec8p7GxEVVVVQ5f1CI0HojsDggLULhV7mrIw0Z174IvHhiLB67uBa1awjcHy/Drlzbjn98dRbPZInd5REQu4VQ4KS8vh9lsRmxsrMPx2NhYlJR0bl+QefPmISEhwSHgnG3x4sUICwuzfyUnJztTpu/r2jJrh4ux+SW9Ro2Hft0bXzwwFiNTI1HfZMbTaw/gxtd/wN4TRrnLIyK6ZB6drbNkyRKsXLkSn376KQwGwznPmz9/PoxGo/2rqKjIg1V6ga4tmwByMTa/1jMmBCvvG4XFNw9EiEGDvSeqcP1r3+Ppz/ejtrFZ7vKIiC6aU+EkKioKarUapaWlDsdLS0sRFxd33se+8MILWLJkCb766isMGjTovOfq9XqEhoY6fFEbthk7J3cBjVygy5+pVBKmjkxB9p+uxHWD4mERwD+/z8f4lzdj40HuwURE3smpcKLT6ZCWlobs7Gz7MYvFguzsbGRkZJzzcc899xyeeuoprF+/HsOHD7/4askqPAUISwaEGTjOadkExIQY8Nodw/DO3SOQGB6AE5X1mLniR8x9fyfKqhvkLo+IyClOd+tkZWXhrbfewrvvvosDBw5g9uzZqK2txcyZMwEA06dPx/z58+3nP/vss3j88cexfPlypKamoqSkBCUlJaip4f/xXxLus0MduKpvDL566ArcO6YbVBLw+e5iZL74LVZuL4TFwmnHROQdnA4nU6ZMwQsvvICFCxdiyJAhyMvLw/r16+2DZAsLC1FcXGw//4033oDJZMKtt96K+Ph4+9cLL7zguqvwR/bF2BhOyFGQXoMF1/XHZ3PGYEBiKKoamvHIJ9Z9eg6XVctdHhHRBTm9zokcuM5JB04fAV4dBqh1wCOFgDZA7opIgZrNFqzYUoAXv/oZ9U1maNUSfj+uJ35/VQ/oNWq5yyMiH+eRdU5IQSK7A8FxgNkEHN8hdzWkUBq1CveOte7Tc1WfaDSZBf6W/Qsm/u07bDt6Wu7yiIg6xHDirSSJXTvUacmRgVh+9wi8dsdQRAXrcfRULab8Yyse+Xg39+khIsVhOPFm9kGx3GeHLkySJFw3KAHZWVdi6kjrwoYrfyzC1S9twn9/4j49RKQcDCfezBZOjv8INDfKWwt5jbBALRbfPAgf/i4DPWOCUV5jwh8/2IWZK35EUUWd3OURETGceLXoPkBgFNDcYF2QjcgJI7tFYu0fx+ChzN7QqVXYdOgUxr+8GW9t5j49RCQvhhNvJkmt++ywa4cugl6jxgOZvbDugbEY2c26T8//rTuAG5b+gD3HuU8PEcmD4cTbdeWgWLp0PWOCsXLWKDx7y0CEGjTYd7IKNyz9Hk/+j/v0EJHnMZx4O9uMncJtgJkfInTxVCoJU0akIPtP43D94ARYBLD8B+s+Pd8cLL3wExARuQjDibeLuQwwhANNtcAxdu3QpYsO0ePvU4dixcwRSIqw7tPz2xU7cNfb2/D57pNoaDLLXSIR+TiGE2+nUgEDbrHe/uFv8tZCPmVcH+s+Pfdd0R1qlYTvfinH3Pd3YeT/fY0Fa/Ygr6iS04+JyC24fL0vOFMA/H2YdZfiWRuBxGFyV0Q+pqC8Fqtzj+PjncdRbGzd5bhHdBBuTUvGTUMTERdmkLFCIlKii/38ZjjxFZ/8Dti9Eug3GZjyntzVkI8yWwRyjpzG6twirN9XgoYm65RjlQSM7RWNW9OS8Ov+sTBouW8PETGcUNlB4PV06+3fbwNi+spbD/m86oYmrNtTjNW5x/FjwRn78RCDBpMHJ+DWtCQMTQ6HJEkyVklEcmI4IWDVncCB/wGDbgduflPuasiPFJTX4uOdx/HJzhM4UVlvP949Ogi3piXh5qFJ7PYh8kMMJwSc2Am8dRUgqYE/7gQiUuWuiPyMxSKw9ehprM49jnV7ix26fcb0isYtwxIx4bI4dvsQ+QmGE7L6903AkW+A4fcA170kdzXkx6obmvDFnhKszj2O7QUV9uMhBg2uG2Tt9hmWwm4fIl/GcEJWBd8DKyYBaj3w4G4gJE7uiohw7HQtPt55Ah/nHnfs9okKwi1pSbh5WCLiwwJkrJCI3IHhhKyEAJZPAIq2AZf/ERj/lNwVEdlZLAJb863dPl/sKUF9y4JukgSM6RmFW9OSML5/HAJ07PYh8gUMJ9Tq5y+B928DdMHAg3uAwEi5KyJqp6axGev2FOPj3OPYlt+m20evwXWD43HLsCSkdY1gtw+RF2M4oVZCAMvGAqV7gHGPAuPmyV0R0XkVnq7Dxzuti7wdP9Pa7dMtKgi3DEvETcOSkBjObh8ib8NwQo72fgKsngkERAAP7gX0wXJXRHRBFovAtvwKfLzzONbtKUadqbXbZ3QPa7fPhMvY7UPkLRhOyJHFDLw2Aqg4Aox/Grj8D3JXROSU2sZmfLG3BKtzi7D1aGu3T7Beg0kD43Hr8CQMZ7cPkaIxnFB7O/8N/HcuEBwHPPAToOUiWOSdiirq8MnOE1i9swhFFa3dPqldAnHLsCTcnMZuHyIlYjih9ppNwN+HAFUngOteBob/Vu6KiC6JxSLwY0EFVucex9qzun0u79EFtwxLwjUD4hCo08hcKREBDCd0LluXAevnAeFdgT/sBNT8o02+obaxGev3luDjncex5chp+/EgnRqTBsXj1rRkDO8aAZWK3T5EcmE4oY6Z6oBXBgJ15cDNbwGDbpO7IiKXK6qow6e7TmB17nEUVtTZj4cYNBiSHI6hKREYmhKOocnhCA/UyVgpkX9hOKFz++5FIPtJILovMDsHUKnkrojILYQQ+LHgDD5u6fapaWxud073qCAMSWkJLMnh6BsXAo2a/yaI3IHhhM6twQi8PABorAKm/Afod53cFRG5XZPZgkMl1dhVVIldhWeQV1iJo+W17c4L0KoxMCmspWUlAsNSwhETysHjRK7AcELnl/2ktQUlYRgw6xvrCEIiP3Om1oS845XYVdgSWIoqUd3QvnUlMTzA2rrS0iV0WUIod1ImuggMJ3R+NaesY0+a64G71gA9rpK7IiLZWSwCR8trsLOwNbD8XFoNy1l/FbVqCf0TwlrCirWFJTkygGusEF0Awwld2BfzgG3LgNSxwN2fy10NkSLVNDZjd0vrSl5Ll1B5jandeV2CdNag0jJ2ZVByOIL1nA1H1BbDCV2Y8Tjwt8GApRm4dTkw4Ba5KyJSPCEEjp+pt49d2VVYiX0njWgyO/7pVElA79gQe8vK0JRw9IgO5lRm8msMJ9Q5Xy0AtrwKqLTA7e8DvcfLXRGR12loMmN/cZW9K2hXYSVOVNa3Oy9Er3EYuzIkORwRQZzKTP6D4YQ6x2IGPpkF7P0Y0BiAOz8GUsfIXRWR1yuramhpXbEGlt3HjahvMrc7r1tUUOvYlZQI9IkLgZZTmclHMZxQ55mbgFV3AT9/AehCgBmfAYlpcldF5FOazRYcKq1uCSuV2FV0BkdPtZ/KbNCqMCjRFlasgSWWU5nJRzCckHOaGoD3fwPkbwYCIoC71wGx/eWuisinGeuaWqYyn7G3sFR1MJU5IczQuqptSjguSwjjVGbySgwn5LzGGuBfNwAndgDBscDML4AuPeSuishvWCwC+adrHcauHCypajeVWaOS0D8h1GHsSnJkINQcbEsKx3BCF6euAnh3MlC6FwhLAX67HghLlLsqIr9V29iMPSeM9sCys7AS5TWN7c7TqCTEhxuQHBGI5IhAJEUEIDmy9Xt0sJ4zhUh2DCd08WrKgOXXABVHgC69rC0owdFyV0VEsE5lPlFZ7zB2Zd/JKpiaLed9nE6jQlJ4AJJsgaVNgEmOCEBkkI6LyJHbMZzQpakssgaUquNA3EBgxudAQLjcVRFRBywWgdLqBhw/U4+iijrH72fqUGxsgPnsvqGzBOrUSIoIQFKENawkRQQiOdL2cyDCArUeuhryZQwndOlOH7EGlNoyIDkduOtTQBckd1VE5KRmswXFxgYUnbEGluNtgktRRT1Kqxtwob/8IQaNPbgkt219ibR+D+JquNQJDCfkGiV7gRXXWncy7j4OmLoK0HJaI5EvaWw242RlA463hJXjZ+pQ1Kb1paMxLmeLCNS2dBFZg8vZ3UecXUQAwwm50vEdwLvXA021QJ9JwG3vAmo28RL5i3qTGScqOw4uRWfqUFnXdMHniA7Rd9hdlBwZgPiwAOg0XHjOHzCckGvlbwbeuxUwNwKDpgA3LgNU/GNCREB1Q5O1u6gltNi6j2wBpqax/dotbakkIC7UgKSWbqKzu4/iwwI4TdpHMJyQ6x36Alh1p3WjwOH3AJNeBDi6n4jOQwgBY31Tm1aX9gN2G5rOP9Oo7TRpW1cRp0l7J4YTco89q4GP7wUggMFTgasXAaHxcldFRF5KCIHyGpO9u6jtuJfjZ+px4kw9TObOTZNObOk2igzSIiJQh/BAHSICtQgP1CEyyHo71KBlkJERwwm5T+67wP/+aL2tCQDSfweMedC67D0RkQtZLAJl1Y0tLS5txr1UdH6adFsqCQgLsIUXxxATEdT2mPW77TYH9LoGwwm517EtwNd/BYq2WX/WhwGj/wiMms3pxkTkMbZp0rYuopOV9aisa8KZOhPO1DWhss6EM3UmVNY2ofoCY1/OJ0CrtrfCRARpWwNN23BzVqgJMWjYSnMWhhNyPyGAn78Esp8EyvZZjwXFAFf8GUi7G9DoZC2PiKgtU7MFlfUma3ipbRteWkPM2ccq65rQ7ETLTFsqCQh3CCxah64mh2NBreFGr/HdVhqGE/IciwXYuxrY+H/AmQLrsfAUYNyjwKDbAJXv/kMjIt8mhEB1YzMqa22tMaYOWmbahJta6+1ak/miXzNQpz6r20nbLsjYwk1koA7hQVqE6DVesf0Awwl5XrMJ2PUv4NvngJpS67HofsCvFgB9J3FmDxH5jcZmM4x1TTjTEmQq60yoqG293TbY2AJPZZ2p3Q7UnaVRSQhv1zLTfmBw27E14QE6j68vw3BC8jHVAdvfBL5/2bqyLAAkDgcyFwHdrpC3NiIihbJYBKobms/TQtPmdm1ruKlvuvhWmmC95pwtNLemJSE5MtCFV8hwQkpQXwls+Tuw9Q2gqc56rPtVwNULgcRhspZGROQrGprMbYJMa6g5e2xNRZv7jPVNF9xP6ePZlyOtq2tnYTKckHJUlwKbnwdyVwCWlmWu+11v7e6J7iNraURE/shsEaiqP8/Ymbom/PFXvRAX5tq91BhOSHnOFAAbFwO7VwFo+c+sSy9rV0/3K4HUsUBgpJwVEhGRGzGckHKV7rfO7Dm4FvaQAgCQgLgBQLcrrV9dMwB9iFxVEhGRizGckPLVnwEKfrBuKpi/GTh1wPF+lQZITLO2rHS7AkgaCWhd28RIRESew3BC3qe6FCj4Dsj/Fjj6LVB5zPF+jQFITm/pBhoHxA8B1Bo5KiUioovAcELe70wBkN8SVvI3t66dYqMLAVJHt3QDXQHE9AdUnp2zT0REncdwQr5FCKD8Z2uLSv631hYW2xoqNgERQJee1tVpHb66AmFJgDZAntqJiAgAwwn5OosZKNndOl7l2JbWtVTOJSim4+ASngKEJzO8EBG5GcMJ+ZdmE1C2H6gs7ODrGGCqufBzMLwQEbnVxX5+c3QheSeNDkgYYv06mxDWmUEdBpc24aW2zPp1YkfHrxEUA4TEAYYwQB8KGEJbvoe1ud32WJvzNAbuLUREdJEuKpwsXboUzz//PEpKSjB48GC8+uqrGDly5DnP/+ijj/D444+joKAAvXr1wrPPPotrr732oosmOi9Jsi7uFhjpmvByMVRaa0jpMNicfayD+3RBgEbPHZ6JyC85HU5WrVqFrKwsLFu2DOnp6XjllVcwYcIEHDp0CDExMe3O37JlC6ZOnYrFixfjuuuuw/vvv48bb7wRO3fuxIABA1xyEURO6XR4OQbUllsH4jYYgcYqoKHqrO/G9scgrMv21522fl0KlQZQ660tRQ7f9YBa5/hdoz/3fed6jvPdd/ZzqLTWsCSp2CpERG7l9JiT9PR0jBgxAq+99hoAwGKxIDk5GX/4wx/wyCOPtDt/ypQpqK2txeeff24/NmrUKAwZMgTLli3r1GtyzAl5DYvF2urSUZhpqOxcwDFVy30VFyapW4KK2hqgVKrzH1NpWm6r2tyvbnO/+qzHd3TsQs/ZwTFJ1RqmJBUA6ayfcYH72/wMqeV4J861/9zJ8+234dxz29+Pludw620PvhbDr8/wyJgTk8mE3NxczJ8/335MpVIhMzMTOTk5HT4mJycHWVlZDscmTJiANWvWnPN1Ghsb0djYaP/ZaLROIa2qqnKmXCL5qMKAgDDgYsbUWszWmUjNTYC5EWhutLbENDcCZpN1MHBzQ/tj5pbbbX9uNgEW288m6/nNJsDS8rzmNs/R0WPNJgCWDopsbvkicrc2waXDUIY296ODYOPu+85xvN3jXHXfxdTYyftufguIHwRXsn1uOzv3xqlwUl5eDrPZjNjYWIfjsbGxOHjwYIePKSkp6fD8kpKSc77O4sWL8cQTT7Q7npyc7Ey5RERE1FmPjnXbU1dXVyMsLKzT5ytyts78+fMdWlssFgsqKirQpUsXSC5s7quqqkJycjKKiop8truI1+gbeI2+gdfoG3iNnSeEQHV1NRISEpx6nFPhJCoqCmq1GqWljsuKl5aWIi4ursPHxMXFOXU+AOj1euj1eodj4eHhzpTqlNDQUJ/9D8yG1+gbeI2+gdfoG3iNneNMi4mNUxuT6HQ6pKWlITs7237MYrEgOzsbGRkZHT4mIyPD4XwA2LBhwznPJyIiIv/mdLdOVlYWZsyYgeHDh2PkyJF45ZVXUFtbi5kzZwIApk+fjsTERCxevBgA8MADD+DKK6/Eiy++iEmTJmHlypXYsWMH/vGPf7j2SoiIiMgnOB1OpkyZglOnTmHhwoUoKSnBkCFDsH79evug18LCQqja7BR7+eWX4/3338eCBQvw6KOPolevXlizZo0i1jjR6/VYtGhRuy4kX8Jr9A28Rt/Aa/QNvEb384q9dYiIiMh/ODXmhIiIiMjdGE6IiIhIURhOiIiISFEYToiIiEhR/DqcLF26FKmpqTAYDEhPT8f27dvlLqlDixcvxogRIxASEoKYmBjceOONOHTokMM548aNgyRJDl/333+/wzmFhYWYNGkSAgMDERMTgz//+c9obnbcH2XTpk0YNmwY9Ho9evbsiRUrVrj78gAAf/3rX9vV37dvX/v9DQ0NmDNnDrp06YLg4GDccsst7Rb3U/L1AUBqamq7a5QkCXPmzAHgne/h5s2bMXnyZCQkJECSpHZ7ZgkhsHDhQsTHxyMgIACZmZn45ZdfHM6pqKjAtGnTEBoaivDwcNxzzz2oqalxOGf37t0YO3YsDAYDkpOT8dxzz7Wr5aOPPkLfvn1hMBgwcOBArFu3zu3X2NTUhHnz5mHgwIEICgpCQkICpk+fjpMnTzo8R0fv/ZIlS7ziGgHg7rvvblf/Nddc43CON7+PADr8tylJEp5//nn7OUp+HzvzOeHJv6OX/Pkq/NTKlSuFTqcTy5cvF/v27ROzZs0S4eHhorS0VO7S2pkwYYJ45513xN69e0VeXp649tprRUpKiqipqbGfc+WVV4pZs2aJ4uJi+5fRaLTf39zcLAYMGCAyMzPFrl27xLp160RUVJSYP3++/ZyjR4+KwMBAkZWVJfbv3y9effVVoVarxfr1691+jYsWLRKXXXaZQ/2nTp2y33///feL5ORkkZ2dLXbs2CFGjRolLr/8cq+5PiGEKCsrc7i+DRs2CABi48aNQgjvfA/XrVsnHnvsMfHJJ58IAOLTTz91uH/JkiUiLCxMrFmzRvz000/i+uuvF926dRP19fX2c6655hoxePBgsXXrVvHdd9+Jnj17iqlTp9rvNxqNIjY2VkybNk3s3btXfPDBByIgIEC8+eab9nN++OEHoVarxXPPPSf2798vFixYILRardizZ49br7GyslJkZmaKVatWiYMHD4qcnBwxcuRIkZaW5vAcXbt2FU8++aTDe9v236+Sr1EIIWbMmCGuueYah/orKioczvHm91EI4XBtxcXFYvny5UKSJHHkyBH7OUp+HzvzOeGpv6Ou+Hz123AycuRIMWfOHPvPZrNZJCQkiMWLF8tYVeeUlZUJAOLbb7+1H7vyyivFAw88cM7HrFu3TqhUKlFSUmI/9sYbb4jQ0FDR2NgohBDiL3/5i7jsssscHjdlyhQxYcIE115ABxYtWiQGDx7c4X2VlZVCq9WKjz76yH7swIEDAoDIyckRQij/+jrywAMPiB49egiLxSKE8P738Ow/+BaLRcTFxYnnn3/efqyyslLo9XrxwQcfCCGE2L9/vwAgfvzxR/s5X3zxhZAkSZw4cUIIIcTrr78uIiIi7NcohBDz5s0Tffr0sf982223iUmTJjnUk56eLn73u9+59Ro7sn37dgFAHDt2zH6sa9eu4uWXXz7nY5R+jTNmzBA33HDDOR/ji+/jDTfcIH71q185HPOm9/HszwlP/h11xeerX3brmEwm5ObmIjMz035MpVIhMzMTOTk5MlbWOUajEQAQGRnpcPw///kPoqKiMGDAAMyfPx91dXX2+3JycjBw4ECHHaInTJiAqqoq7Nu3z35O29+J7RxP/U5++eUXJCQkoHv37pg2bRoKCwsBALm5uWhqanKorW/fvkhJSbHX5g3X15bJZMJ7772H3/72tw6bWXr7e9hWfn4+SkpKHOoJCwtDenq6w/sWHh6O4cOH28/JzMyESqXCtm3b7OdcccUV0Ol09nMmTJiAQ4cO4cyZM/ZzlHLdRqMRkiS12w9syZIl6NKlC4YOHYrnn3/eoancG65x06ZNiImJQZ8+fTB79mycPn3aoX5feh9LS0uxdu1a3HPPPe3u85b38ezPCU/9HXXV56sidyV2t/LycpjNZoc3AABiY2Nx8OBBmarqHIvFggcffBCjR492WGX3jjvuQNeuXZGQkIDdu3dj3rx5OHToED755BMAQElJSYfXa7vvfOdUVVWhvr4eAQEBbruu9PR0rFixAn369EFxcTGeeOIJjB07Fnv37kVJSQl0Ol27P/axsbEXrN123/nO8cT1nW3NmjWorKzE3XffbT/m7e/h2Ww1dVRP23pjYmIc7tdoNIiMjHQ4p1u3bu2ew3ZfRETEOa/b9hye0tDQgHnz5mHq1KkOm6X98Y9/xLBhwxAZGYktW7Zg/vz5KC4uxksvvWS/DiVf4zXXXIObb74Z3bp1w5EjR/Doo49i4sSJyMnJgVqt9rn38d1330VISAhuvvlmh+Pe8j529Dnhqb+jZ86cccnnq1+GE282Z84c7N27F99//73D8fvuu89+e+DAgYiPj8fVV1+NI0eOoEePHp4u02kTJ0603x40aBDS09PRtWtXfPjhhx79QPWUt99+GxMnTnTYRtzb30N/19TUhNtuuw1CCLzxxhsO92VlZdlvDxo0CDqdDr/73e+wePFir1gC/fbbb7ffHjhwIAYNGoQePXpg06ZNuPrqq2WszD2WL1+OadOmwWAwOBz3lvfxXJ8T3sQvu3WioqKgVqvbjVIuLS1FXFycTFVd2Ny5c/H5559j48aNSEpKOu+56enpAIDDhw8DAOLi4jq8Xtt95zsnNDTU4wEhPDwcvXv3xuHDhxEXFweTyYTKysp2tV2odtt95zvH09d37NgxfP3117j33nvPe563v4e2ms737ywuLg5lZWUO9zc3N6OiosIl762n/j3bgsmxY8ewYcOGC24xn56ejubmZhQUFADwjmtsq3v37oiKinL4b9MX3kcA+O6773Do0KEL/vsElPk+nutzwlN/R131+eqX4USn0yEtLQ3Z2dn2YxaLBdnZ2cjIyJCxso4JITB37lx8+umn+Oabb9o1G3YkLy8PABAfHw8AyMjIwJ49exz+gNj+iPbv399+Ttvfie0cOX4nNTU1OHLkCOLj45GWlgatVutQ26FDh1BYWGivzZuu75133kFMTAwmTZp03vO8/T3s1q0b4uLiHOqpqqrCtm3bHN63yspK5Obm2s/55ptvYLFY7OEsIyMDmzdvRlNTk/2cDRs2oE+fPoiIiLCfI9d124LJL7/8gq+//hpdunS54GPy8vKgUqnsXSFKv8azHT9+HKdPn3b4b9Pb30ebt99+G2lpaRg8ePAFz1XS+3ihzwlP/R112edrp4fO+piVK1cKvV4vVqxYIfbv3y/uu+8+ER4e7jBKWSlmz54twsLCxKZNmxymsNXV1QkhhDh8+LB48sknxY4dO0R+fr747LPPRPfu3cUVV1xhfw7bFLHx48eLvLw8sX79ehEdHd3hFLE///nP4sCBA2Lp0qUem2r7pz/9SWzatEnk5+eLH374QWRmZoqoqChRVlYmhLBOgUtJSRHffPON2LFjh8jIyBAZGRlec302ZrNZpKSkiHnz5jkc99b3sLq6WuzatUvs2rVLABAvvfSS2LVrl32mypIlS0R4eLj47LPPxO7du8UNN9zQ4VTioUOHim3btonvv/9e9OrVy2EKamVlpYiNjRV33XWX2Lt3r1i5cqUIDAxsNz1To9GIF154QRw4cEAsWrTIZVNQz3eNJpNJXH/99SIpKUnk5eU5/Pu0zW7YsmWLePnll0VeXp44cuSIeO+990R0dLSYPn26V1xjdXW1ePjhh0VOTo7Iz88XX3/9tRg2bJjo1auXaGhosD+HN7+PNkajUQQGBoo33nij3eOV/j5e6HNCCM/9HXXF56vfhhMhhHj11VdFSkqK0Ol0YuTIkWLr1q1yl9QhAB1+vfPOO0IIIQoLC8UVV1whIiMjhV6vFz179hR//vOfHdbIEEKIgoICMXHiRBEQECCioqLEn/70J9HU1ORwzsaNG8WQIUOETqcT3bt3t7+Gu02ZMkXEx8cLnU4nEhMTxZQpU8Thw4ft99fX14vf//73IiIiQgQGBoqbbrpJFBcXOzyHkq/P5ssvvxQAxKFDhxyOe+t7uHHjxg7/25wxY4YQwjqd+PHHHxexsbFCr9eLq6++ut21nz59WkydOlUEBweL0NBQMXPmTFFdXe1wzk8//STGjBkj9Hq9SExMFEuWLGlXy4cffih69+4tdDqduOyyy8TatWvdfo35+fnn/PdpW78mNzdXpKeni7CwMGEwGES/fv3EM8884/DBruRrrKurE+PHjxfR0dFCq9WKrl27ilmzZrX7oPHm99HmzTffFAEBAaKysrLd45X+Pl7oc0IIz/4dvdTPV6nlooiIiIgUwS/HnBAREZFyMZwQERGRojCcEBERkaIwnBAREZGiMJwQERGRojCcEBERkaIwnBAREZGiMJwQERGRojCcEBERkaIwnBAREZGiMJwQERGRojCcEBERkaL8P5906cC81NXZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(l,cl_eg_1h,label='1h')\n",
    "plt.plot(l,cl_eg_2h,label='2h')\n",
    "plt.legend()\n",
    "# plt.yscale('log')\n",
    "# plt.yscale('log')\n",
    "plt.ylim(0,1e-10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "dcb13f70-3257-4d3e-8ed5-b10c9d12dbe5",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\"cleg_l_1h_2h_1p2h_shock.txt\",np.c_[l,cl_eg_1h,cl_eg_2h,cl_eg_1h+cl_eg_2h])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0b226fe-38e4-46a1-91ea-a676f66629be",
   "metadata": {},
   "source": [
    "# scaled NFW"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "59eac1bb-1027-4819-9984-f3c1c58dbdaa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3 s, sys: 131 ms, total: 3.13 s\n",
      "Wall time: 585 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "M = Class()\n",
    "\n",
    "M.set(common_params)\n",
    "M.set(cosmo_params)\n",
    "M.set(ksz_params)\n",
    "M.set(HOD_blue)\n",
    "\n",
    "M.set({\n",
    "'use_fft_for_profiles_transform' : 1,\n",
    "'ndim_redshifts':30, ## precision parameter -- how many z in redshift grid.\n",
    "'gas profile':'nfw'\n",
    "})\n",
    "\n",
    "M.compute_class_szfast() # fast mode works\n",
    "# M.compute() # slow mode works"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8992e420-59f2-4721-95ba-00eeae8c2892",
   "metadata": {},
   "outputs": [],
   "source": [
    "l = np.asarray(M.cl_sz()['ell'])\n",
    "cl_eg_1h = np.asarray(M.cl_eg()['1h'])/l_to_dl(l)\n",
    "cl_eg_2h = np.asarray(M.cl_eg()['2h'])/l_to_dl(l)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a19055f4-5187-4c2d-a6e5-9e666a2cbd46",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.00000000e+00, 2.44280552e+00, 2.98364940e+00, 3.64423760e+00,\n",
       "       4.45108186e+00, 5.43656366e+00, 6.64023385e+00, 8.11039993e+00,\n",
       "       9.90606485e+00, 1.20992949e+01, 1.47781122e+01, 1.80500270e+01,\n",
       "       2.20463528e+01, 2.69274761e+01, 3.28892935e+01, 4.01710738e+01,\n",
       "       4.90650604e+01, 5.99282001e+01, 7.31964689e+01, 8.94023690e+01,\n",
       "       1.09196300e+02, 1.33372662e+02, 1.62901737e+02, 1.98968631e+02,\n",
       "       2.43020835e+02, 2.96826318e+02, 3.62544484e+02, 4.42812832e+02,\n",
       "       5.40852815e+02, 6.60599120e+02, 8.06857587e+02, 9.85498082e+02,\n",
       "       1.20369008e+03, 1.47019038e+03, 1.79569458e+03, 2.19326632e+03,\n",
       "       2.67886153e+03, 3.27196886e+03, 3.99639179e+03, 4.88120396e+03,\n",
       "       5.96191597e+03, 7.28190061e+03, 8.89413350e+03, 1.08633192e+04,\n",
       "       1.32684880e+04, 1.62061679e+04, 1.97942581e+04])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "6493cd44-c7bc-4753-a629-1379b79ad8e7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 1e-10)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGsCAYAAAAGzwdbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABR7klEQVR4nO3deXzT9f0H8Ffu9D7pSUvLDXJfpSCos6MiQ5lzIjJBvCbCPDodogI/jwmeY1MU50R0TkGc4CYMxQqCUKgUKiIFOXoBPWhL07tpk8/vjzShoaU0pcn3m+T1fDzySPrNN/m+v400Lz/fz6EQQggQERERyYRS6gKIiIiIWmM4ISIiIllhOCEiIiJZYTghIiIiWWE4ISIiIllhOCEiIiJZYTghIiIiWWE4ISIiIllhOCEiIiJZYTghIiIiWXGrcLJz505Mnz4dMTExUCgU2LRpk+THE0Jg6dKliI6Oho+PD1JSUnD8+HGn1kVEROTJ3Cqc1NbWYvjw4Vi1apVsjvfSSy/hb3/7G1avXo19+/bBz88PqampaGhocEmNREREnkbhrgv/KRQKbNy4ETNmzLBta2xsxFNPPYWPP/4YlZWVGDJkCF588UVce+21TjmeEAIxMTH44x//iMceewwAYDAYEBkZibVr1+L222+/4uMSERF5G7dqObmchQsXIiMjA+vWrcOhQ4fw29/+FjfccIPTLrPk5uaiuLgYKSkptm1BQUFISkpCRkaGU45JRETk6TwmnBQUFOC9997Dhg0bMGnSJPTp0wePPfYYrr76arz33ntOOWZxcTEAIDIy0m57ZGSk7TkiIiJyjMeEkx9//BEmkwn9+/eHv7+/7fbtt9/i5MmTAICjR49CoVB0eHviiSckPhMiIiLvppa6gO5SU1MDlUqFrKwsqFQqu+f8/f0BAL1790ZOTk6H7xMWFtbpY0ZFRQEASkpKEB0dbdteUlKCESNGdPp9iIiI6AKPCScjR46EyWRCaWkpJk2a1O4+Wq0WAwcO7LZjJiYmIioqCunp6bYwUlVVhX379mH+/PnddhwiIiJv4lbhpKamBidOnLD9nJubi+zsbISGhqJ///6YPXs25syZg1dffRUjR47EuXPnkJ6ejmHDhmHatGnderz4+HgoFAo88sgjeP7559GvXz8kJiZiyZIliImJsRvVQ0RERJ3nVkOJd+zYgeuuu67N9rlz52Lt2rVoamrC888/jw8++ABnzpxBeHg4xo8fj2eeeQZDhw7t9uMBluHEy5Ytw9///ndUVlbi6quvxptvvon+/fs7fDwiIiLqQjjZuXMnXn75ZWRlZaGoqKjN3B/t2bFjB9LS0vDTTz8hLi4OTz/9NO66664rKJuIiIg8lcOjdRydpTU3NxfTpk3Dddddh+zsbDzyyCO499578eWXXzpcLBEREXm+K7qs096sqRdbtGgRNm/ejMOHD9u23X777aisrMTWrVu7emgiIiLyUE7vEJuRkWE3gyoApKam4pFHHrnkaxobG9HY2Gj72Ww2o6KiAmFhYVAoFM4qlYiIiLqREALV1dWIiYmBUtn5izVODyfFxcXtzqBaVVWF+vp6+Pj4tHnN8uXL8cwzzzi7NCIiInKBwsJC9OzZs9P7y3Io8eLFi5GWlmb72WAwID4+HoWFhQgMDOy249z3wX5knCzHiluG4lfDY7rtfXFyO/DJnUDEYOCer7rvfa/AkbMG3Pb2XoT5afDtn34hdTlEROQFqqqqEBcXh4CAAIde5/RwEhUVhZKSErttJSUlCAwMbLfVBAB0Oh10Ol2b7YGBgd0aTrQ+flDq6uHrH9Ct74uQcECnANRNQHe+7xUYrPWBUueL882AQuuDAL1G6pKIiMhLONolw+lr6yQnJyM9Pd1u27Zt25CcnOzsQ0tH0xK6mhqkraOVQL0GYX5aAEB+eZ3E1RAREV2aw+GkpqYG2dnZyM7OBnBh1tSCggIAlksyc+bMse3/wAMP4NSpU/jTn/6Eo0eP4s0338Qnn3yCRx99tHvOQI5s4UReIaBXmC8AhhMiIpI3h8PJ/v37MXLkSIwcORIAkJaWhpEjR2Lp0qUAgKKiIltQASzrz2zevBnbtm3D8OHD8eqrr+If//gHUlNTu+kUZMgaTprl03ICAAlhfgCAvPJaiSshIiK6NIf7nFx77bXoaGoU67TuF7/m4MGDjh7KfalbtZwIAchk+HOvlnCSz3BCRHTFhBBobm6GyWSSuhTJqFQqqNXqbp/mQ5ajddyeplVH3+ZGQKOXrpZWEsItl3XyeFmHiOiKGI1GFBUVoa6Of099fX0RHR0NrVbbbe/JcOIMrcNJU51swglbToiIrpzZbEZubi5UKhViYmKg1Wq9coJQIQSMRiPOnTuH3Nxc9OvXz6GJ1jrCcOIMKg2gVAPmZln1O0lo6RBbUtWIOmMzfLX8+ImIHGU0GmE2mxEXFwdfX1+py5GUj48PNBoN8vPzYTQaodd3z/+MO30osZx1fVWhTrD1O6l34kEcE+yrRZCPZX6Tggo2RRIRXYnuaiVwd874PfA36ywa+YUTAEgIt1zaOVnKSztERCRPDCfOYu1nIrNwMijKMoVwTlGVxJUQERG1j+HEWTQt1yGbZRZOoi3T6R9hOCEiooskJCRg5cqVUpfBcAI4aRoStTxbTgbHtISTswwnRETeZufOnZg+fTpiYmKgUCiwadMmqUtqF8OJs1hbTmQWTga2XNYprmpARa1R4mqIiMiVamtrMXz4cKxatUrqUjrEsaTOItM+JwF6DXqF+SK/vA45RVWY2Ddc6pKIiNyeEAL1Ta6fKdZHo3JojpWpU6di6tSpHe5TV1eHu+++Gxs2bEBISAiefvpp3H///VdaqkMYTpxFpn1OAGBwdCDyy+tw5CzDCRFRd6hvMmHw0i9dftwjz6Z2+5xVr776Kp577jk8+eST+PTTTzF//nxcc801GDBgQLcepyO8rOMsMu1zAljCCcBOsURE1NaNN96IBx98EH379sWiRYsQHh6O7du3u7QGtpw4i0znOQHYKZaIqLv5aFQ48myqJMftbsOGDbM9VigUiIqKQmlpabcfpyNeHU4EnDhFrBuEkxPnatDQZILeCf9xExF5E4VC4TFLgmg0GrufFQoFzGazS2vgZR1nsYYTGa2tYxUVqEeIrwYms8DxkhqpyyEiIrLDcOIstrV15LeGjUKhaDUZm0HiaoiIyFVqamqQnZ2N7OxsAEBubi6ys7NRUFAgbWEXYThxFttlHfm1nACtOsWy3wkRkdfYv38/Ro4ciZEjRwIA0tLSMHLkSCxdulTiyux5xgUyOdLIt+UEaNUpliN2iIi8xrXXXgshLt3fMi8vr802ayuLK7HlxFlk3OcEuBBOcoqqYTY7sWMwERGRgxhOnEXGfU4AoE8Pf2hVStQ0NuP0efmNKCIiIu/FcOIsMu9zolEp0T/KHwA7xRIRkbwwnDiLTNfWaY2dYomISI4YTpxFxmvrWHEaeyIikiOvDicddFi+cjJeW8dqcEwQALacEBGRvHh1OHEqa8uJjMPJwOgAAMBZQwPO1xolroaIiMiC4cRZ3KDPSaBeg7hQS8fdHF7aISIimWA4cZbWfU6cev3oyrDfCRERyQ3DibNY+5wAsp2IDQAGR7f0O2E4ISIiWNZf27Rpk6Q1MJw4i3WeE0DWl3Zs09izUywRkcdbvnw5xo4di4CAAERERGDGjBk4duyY1GW1wXDiLCoNoGxZusgNwsmJ0ho0NpskroaIiJzp22+/xYIFC7B3715s27YNTU1NmDJlCmpra6UuzQ4X/nMmjS/QWCXryzoxQXoE+WhgqG/C8ZIaDIkNkrokIiL3I4Q0y5VofAGFotO7b9261e7ntWvXIiIiAllZWZg8ebJte1lZGX7961/jyy+/RGxsLF599VXcdNNN3Vb25TCcOJNabwknMl1fB7BcWxwcHYiMU+U4UlTFcEJE1BVNdcALMa4/7pNnAa1fl19uMFiWLwkNDbXb/swzz+Cll17Cyy+/jNdffx2zZ89Gfn5+m/2chZd1nEnm6+tYsd8JEZH3MZvNeOSRRzBx4kQMGTLE7rm77roLs2bNQt++ffHCCy+gpqYGmZmZLqvNq1tOnD7CVyPvlYmtBnE4MRHRldH4WloxpDhuFy1YsACHDx/Gd9991+a5YcOG2R77+fkhMDAQpaWlXT6Wo7w6nFgpHLhe5xBrOJFxnxPgwlwnOUVVEEI47/dBROSpFIorurziagsXLsQXX3yBnTt3omfPnm2e12g0dj8rFAqYzWZXlcfLOk6ldo+Wk74R/tCoFKhuaMbp8/IdWURERFdGCIGFCxdi48aN+Oabb5CYmCh1Se1iOHEmN+lzolUr0S/Css4OL+0QEXmuBQsW4MMPP8RHH32EgIAAFBcXo7i4GPX18vofU4YTZ3KTPicAO8USEXmDt956CwaDAddeey2io6Ntt/Xr10tdmh32OXEmN+lzAnCNHSIibyA6MRKkvX0qKyudUM2lseXEmdhyQkRE5DCGE2dSu0efE+DCcOIzlfUw1DVJXA0REXkzhhNnsrWcyKujUXuCfDToGWKpl5d2iIhISgwnzmTrcyL/cAJwMjYiIpIHhhNncqOWE8B+MjYiIupYZzqXegNn/B4YTpxJ7WbhhJ1iiYguyzp7al2d/Ac7uIL193DxrLJXgkOJnclNW06Ol1bD2GyGVs3sSkR0MZVKheDgYNtaM76+vl657IcQAnV1dSgtLUVwcDBUKlW3vTfDiTO5WZ+TniE+CNCrUd3QjBOlNbaWFCIishcVFQUALl0MT66Cg4Ntv4/uwnACwGl5181aThQKBQZHB2JfbgWOFFUxnBARXYJCoUB0dDQiIiLQ1OS90y9oNJpubTGxYjhxJrXecu8G85xYDY5pCSdnq4DRUldDRCRvKpXKKV/O3o6dCpxJ42u5d4MZYq0uTGNvkLgSIiLyVl4dTgScPAxM09Jy4gZr61i1HrHDYXJERCQFrw4nTueGLSd9I/yhVipQ1dCMswb3CVVEROQ5GE6cyQ37nOjUKvSN8AfA+U6IiEgaDCfOZG05aa4H3OgSCSdjIyIiKTGcOJO1zwngXv1O2CmWiIgkxHDiTNbp6wG3mesEaNVywjV2iIhIAgwnzqRSA8qWtQbcKZy0tJwUVtTDUO+9kwsREZE0GE6czTZix33CSbCvFrHBllafo2w9ISIiF2M4cTbbXCfuE04AYFA0L+0QEZE0GE6czc3W17HiiB0iIpIKwwkAp650rXbTcBIdAADIKWY4ISIi12I4cTZ3bTmJDgIA/FxcgyaTWeJqiIjImzCcOJs1nLhZn5OeIT4I0KlhNJlx8lyN1OUQEZEX6VI4WbVqFRISEqDX65GUlITMzMwO91+5ciUGDBgAHx8fxMXF4dFHH0VDg/tMSnZF3LTlRKlUXOgUy34nRETkQg6Hk/Xr1yMtLQ3Lli3DgQMHMHz4cKSmpqK0tLTd/T/66CM88cQTWLZsGXJycvDuu+9i/fr1ePLJJ6+4eLdgW1/HvcIJwE6xREQkDYfDyWuvvYb77rsP8+bNw+DBg7F69Wr4+vpizZo17e6/Z88eTJw4EXfccQcSEhIwZcoUzJo167KtLR7DDec5sRrM4cRERCQBh8KJ0WhEVlYWUlJSLryBUomUlBRkZGS0+5oJEyYgKyvLFkZOnTqFLVu24MYbb7zkcRobG1FVVWV3c1tuOs8JYD+NvXCjhQuJiMi9qR3ZuaysDCaTCZGRkXbbIyMjcfTo0XZfc8cdd6CsrAxXX301hBBobm7GAw880OFlneXLl+OZZ55xpLQuccn3rRu3nPSN8IdaqUBlXROKDA2ICfa5/IuIiIiukNNH6+zYsQMvvPAC3nzzTRw4cACfffYZNm/ejOeee+6Sr1m8eDEMBoPtVlhY6OwyncfW58T9OgDrNSr0jfAHAOTw0g4REbmIQy0n4eHhUKlUKCkpsdteUlKCqKiodl+zZMkS3Hnnnbj33nsBAEOHDkVtbS3uv/9+PPXUU1Aq2+YjnU4HnU7nSGnyZWs5qZO2ji4aFB2Io8XVOHK2CtcPirz8C4iIiK6QQy0nWq0Wo0ePRnp6um2b2WxGeno6kpOT231NXV1dmwCiUqkAwDv6Mdj6nLhfywnATrFEROR6DrWcAEBaWhrmzp2LMWPGYNy4cVi5ciVqa2sxb948AMCcOXMQGxuL5cuXAwCmT5+O1157DSNHjkRSUhJOnDiBJUuWYPr06baQ4tHcvOWkdadYIiIiV3A4nMycORPnzp3D0qVLUVxcjBEjRmDr1q22TrIFBQV2LSVPP/00FAoFnn76aZw5cwY9evTA9OnT8ec//7n7zuIKKeDExXXcuM8JcGF14vzyOlQ3NCFAr5G4IiIi8nQOhxMAWLhwIRYuXNjuczt27LA/gFqNZcuWYdmyZV05lPuzzRDrni0noX5aRAfpUWRowNHiaoxNCJW6JCIi8nBcW8fZbGvruGfLCdCq3wlniiUiIhdgOHE2N11bpzVOY09ERK7EcOJsag8IJy0tJz+eMUhcCREReQOGE2fzgJaT0b1CAAA5xVWorDNKXA0REXk6hhNns/U5cd9wEhGoR78IfwgB7D1VLnU5RETk4RhOnM0DWk4AYEKfMADA7hMMJ0RE5FwMJ86mbjVax2yWtpYrMKFvOABgz8kyiSshIiJPx3DibJpWK/m68XDi8YlhUCqAk+dqUWxw3/MgIiL5YzhxNg8JJ0G+GgyJDQIAZJxi6wkRETkPw4mzKVWASmt57KazxFols98JERG5gFeHE5etiWzrFOu+LScAMLGPpd9Jxsly71hRmoiIJOHV4cRK4cR1/wC0mojNvVtOxiSEQKNS4ExlPfLL3ftciIhIvhhOXMED1tcBAF+tGiPjLROy7TnJSztEROQcDCeu4OYrE7dmm++EQ4qJiMhJGE5cwUP6nADAxJb5TvaeLIfZzH4nRETU/RhOXMFD+pwAwPCewfDRqFBea8SxkmqpyyEiIg/EcOIKHtLnBAC0aiXGJYYCYL8TIiJyDoYTV9DoLfce0HICXOh3sucE+50QEVH3YzhxBY2v5d4D+pwAF/qd7MutQLPJfdcLIiIieWI4cQW1teXEvVcmthoUHYggHw1qGptx6IxB6nKIiMjDMJy4grXlpNkzwolKqUByb8ulnQz2OyEiom7m3eHEVSNhNZ7VcgIAE/pa19lhvxMiIupe3h1OXMXW58SDwknLOjv788+jockkcTVERORJGE5cwcP6nABAnx5+iAzUwdhsxoH881KXQ0REHoThxBVs85x4TjhRKBS21hPOd0JERN2J4cQVbNPXe044AbjODhEROQfDCQCFsw/gqeGkZb6TQ6cNqG5okrgaIiLyFAwnrqD2zHASG+yDhDBfmMwCmbkVUpdDREQeguHEFTxobZ2LJbPfCRERdTOGE1fQeM6qxBebyPlOiIiomzGcuIItnHhgy0nLTLFHi6tRXtMocTVEROQJGE5cQe25LSdh/joMjAoAAGSc4qUdIiK6cl4dToSr5q/34D4nADjfCRERdSuvDicu0zqcmM3S1uIE1n4ne9jvhIiIugHDiStYwwngka0n4xJDoVIqkFdehzOVnjVcmoiIXI/hxBXUrcKJh811AgABeg2G9QwCwNYTIiK6cgwnrqBUAiqd5bEHra/TmnUq+wz2OyEioivEcOIqGs9bmbi1iS2dYnefLIMQLupoTEREHonhxFU0vpZ7Dw0no3qFQKtWoqSqEafKaqUuh4iI3BjDiauoPbvlRK9RYUyvEADsd0JERFeG4QSAwunLEuNCy4mH9jkBLvQ74XwnRER0JRhOXMXD+5wAwIS+ln4nGafKYTaz3wkREXWNV4cTl/bbtPU58bwp7K2GxQbBX6dGZV0TjhRVSV0OERG5Ka8OJy5l63PieZOwWalVSiQlhgIA9pxkvxMiIuoahhNX0Xju4n+tJbPfCRERXSGGE1fx8MX/rCa29DvJzK2Asdnz1hEiIiLnYzhxFVvLied2iAWAAZEBCPXTos5owqHTlVKXQ0REbojhxFXU3hFOlEqF7dLO7hO8tENERI5jOHEVL2k5AeynsiciInIUw4mr2PqceH44sU7GdrDgPOqNJomrISIid8Nw4ipe1HLSK8wXscE+aDIJfJ9XIXU5RETkZhhOXMVLhhIDgEKh4JBiIiLqMoYTV9EFWu4bq6Wtw0Um9rWGE/Y7ISIix3h1OHHp6i8+lhV7UX/elUeVzISWTrGHzxhgqGuSuBoiInInXh1OLnDBssReFk4iA/Xo08MPZgHszeWlHSIi6jyGE1exhZNKSctwJWvrSQb7nRARkQMYTlxFH2y5b6wCTM2SluIq1n4nu0+w3wkREXUew4mr6IMuPG4wSFeHC43vHQaFAjheWoPSas9eU4iIiLoPw4mrqNSAriWgeEm/k2BfLa6KsYxS4qUdIiLqLIYTV/IJttx7STgBLvQ72cN1doiIqJMYTlzJK8NJS78TzndCRESdxHDiStYROw2VkpbhSmMTQqFWKnD6fD0KKzx/dlwiIrpyDCeu5GVznQCAn06NkfHBADhqh4iIOqdL4WTVqlVISEiAXq9HUlISMjMzO9y/srISCxYsQHR0NHQ6Hfr3748tW7Z0qWC35oXhBACSrf1O2CmWiIg6weFwsn79eqSlpWHZsmU4cOAAhg8fjtTUVJSWlra7v9FoxC9/+Uvk5eXh008/xbFjx/DOO+8gNjb2iot3O14aTia2WgRQCJcuGkBERG5I7egLXnvtNdx3332YN28eAGD16tXYvHkz1qxZgyeeeKLN/mvWrEFFRQX27NkDjUYDAEhISLiyqt2VdSI2LwsnI+KDodcoUVbTiJ9LajAgKkDqkoiISMYcajkxGo3IyspCSkrKhTdQKpGSkoKMjIx2X/Of//wHycnJWLBgASIjIzFkyBC88MILMJlMlzxOY2Mjqqqq7G4ewQunsAcAnVqF5N6W1pPNPxZJXA0REcmdQ+GkrKwMJpMJkZGRdtsjIyNRXFzc7mtOnTqFTz/9FCaTCVu2bMGSJUvw6quv4vnnn7/kcZYvX46goCDbLS4uzpEyHaZwwbp/ALz2sg4AzBhpuYz376zTMJt5aYeIiC7N6aN1zGYzIiIi8Pe//x2jR4/GzJkz8dRTT2H16tWXfM3ixYthMBhst8LCQmeX6RpeHE5Sr4pCgF6NM5X12HuKHWOJiOjSHAon4eHhUKlUKCkpsdteUlKCqKiodl8THR2N/v37Q6VS2bYNGjQIxcXFMBqN7b5Gp9MhMDDQ7uYRvHASNiu9RoXpw2MAAJ9mnZa4GiIikjOHwolWq8Xo0aORnp5u22Y2m5Geno7k5OR2XzNx4kScOHECZrPZtu3nn39GdHQ0tFptF8t2U61bTrxw1MpvR/cEAGw5XITqhiaJqyEiIrly+LJOWloa3nnnHbz//vvIycnB/PnzUVtbaxu9M2fOHCxevNi2//z581FRUYGHH34YP//8MzZv3owXXngBCxYs6L6zcBfWcCJMgLFG2lokMCIuGH16+KGhyYwt7BhLRESX4PBQ4pkzZ+LcuXNYunQpiouLMWLECGzdutXWSbagoABK5YXMExcXhy+//BKPPvoohg0bhtjYWDz88MNYtGhR951FF7l8zg2ND6DWA80NltYTnXcNqVUoFLh1dBxe3HoUn2adxsyx8VKXREREMuRwOAGAhQsXYuHChe0+t2PHjjbbkpOTsXfv3q4cyvP4hADVRZZwEux9X863jIrFy18exfd555FbVovEcD+pSyIiIpnh2jqu5qUTsVlFBuoxuX8PAJZhxURERBdjOHE1Lx5ObPXb0ZZ5a/594DRMnPOEiIguwnDial46S2xr1w+KQJCPBkWGBuw5yZWKiYjIHsOJq7HlBHqNCjdxzhMiIroEhhNX8+KJ2Fr77RjLnCdbDxfDUM85T4iI6AKGEwCuWloHAMNJi6GxQegf6Y/GZjM2H+KcJ0REdAHDiavxsg4Ay5wn1o6xG7I8ZO0kIiLqFgwnrmYNJw0GaeuQgZtHxkClVOBgQSVOlHrfjLlERNQ+hhNXY8uJTUSAHtcNaJnz5AA7xhIRkQXDiat5+SRsF7u1ZTHAzzjnCRERtWA4cTW2nNj5xcBIhPhqUFLViF3Hz0ldDhERyQDDiatZw0lTHdDcKG0tMqBVK3HziFgAwAbOeUJERGA4cT1dIKBo+bV78SyxrVkv7Wz7qQSGOs55QkTk7RhOXE2pZL+TiwyJDcKg6EAYTWb859BZqcshIiKJMZxIgROxtWFtPfl0P+c8ISLydgwnUmCn2DZmjIiBWqnAD6cN+LmkWupyiIhIQl4dTkbFh+DaAT0Q5q9z7YFtE7FVuva4Mhbmr8MvBkYA4GKARETezqvDydO/Goy188ZhdK8Q1x6YLSftujDnyRk0m8wSV0NERFLx6nAiGXaIbdd1AyMQ5qdFWU0jdnLOEyIir8VwIgW2nLRLo1JixsiWOU/289IOEZG3YjiRAsPJJVkv7XydU4LztUaJqyEiIikwnEjBFk4qJS1DjgZFB2JIbCCaTAKfZ5+RuhwiIpIAw4kU2HLSoVtHtcx5wpWKiYi8EsOJFDgJW4duHhELjUqBw2eqkFNUJXU5RETkYgwnUmDLSYdC/LRIGRQJgHOeEBF5I4YTKdgmYTMAZs7n0Z7fjrFc2tl08AyaOOcJEZFXYTiRgnWeEwig0SBlJbI1uV8P9AjQobzWiO1HS6Uuh4iIXIjhRApqLaD1tzzmpZ12qVVK3NIy5wkv7RAReReGE6lwltjL+k3LnCffHC1FWU2jxNUQEZGrMJxIhZ1iL6t/ZACG9wxCs1ng8+yzUpdDREQuwnAiFdtw4kopq5C9W8fEAeClHSIib8JwIhW2nHTKTcNioFUrkVNUhcNn2HmYiMgbMJxIhS0nnRLkq8GUwZzzhIjImzCcSIUtJ51mXQzw8+wzMDZzzhMiIk/HcCIVhpNOm9SvByIDdThf14RvjpZIXQ4RETkZw4lUbLPEVkpahjtQKRW4xboYIC/tEBF5PIYTqbDlxCHWSzvbj51DaXWDxNUQEZEzMZxIhZOwOaRPD3+Mig+GySzw+UHOeUJE5MkYTqTClhOH3TraMufJhqxCCCEkroaIiJyF4UQqrcMJv2g75VfDo6FTK/FzSQ1+5JwnREQei+FEKtZwYjICTfXS1uImAvUa3DAkCgA7xhIReTKGE6lo/QClxvKYl3Y67cKcJ2fR0GSSuBoiInIGhhOpKBStZollOOmsCX3CEROkh6G+Cek5pVKXQ0RETsBwIiV2inVY6zlPNmQVSlwNERE5A8OJlBhOusR6aWfnz+dQUsU5T4iIPA3DiZQ4S2yXJIT7YWxCCMwC+CAjT+pyiIiomzGcSIkTsXXZPVf3BgD8Y1cuzlRytBMRkSdhOJESL+t0WepVkUhKDEVjsxkvbT0qdTlERNSNGE6kxHDSZQqFAkt+NRgKhWVY8YEC/g6JiDwFw4mUGE6uyJDYINzaMnLn2f8e4ZT2REQeguFESrZwUilpGe7s8dQB8NWqkF1Yif/8wAUBiYg8AcOJlDgJ2xWLCNTjwWv7AABe/N9R1Bs5aywRkbtjOJESW066xb2TeiM22AdnDQ34x65TUpdDRERXiOFESuxz0i30GhUWTR0IAHhzx0lOzEZE5OYYTqRkDSfGasDUJG0tbm76sGiMig9GfZMJL395TOpyiIjoCjCcSEkfdOFxg0G6OjyAdWgxAHyadRo/nubvk4jIXTGcSEmpAnQtAYWXdq7YyPgQzBgRAwB47gsOLSYiclcMJ1LjiJ1u9acbBkKvUSIzrwJbDxdLXQ4REXUBw4nU2Cm2W8UE++D+yZahxS/8LwcNTRxaTETkbhhOpMZw0u0euKY3IgN1KKyox9o9eVKXQ0REDmI4kZrtsk6llFV4FF+tGn9KtQwtfuObEzhX3ShxRURE5AiGE6mx5cQpfj0yFsN6BqGmsRmvbftZ6nKIiMgBXQonq1atQkJCAvR6PZKSkpCZmdmp161btw4KhQIzZszoymE9E8OJUyiVF4YWr/++ADlFVRJXREREneVwOFm/fj3S0tKwbNkyHDhwAMOHD0dqaipKS0s7fF1eXh4ee+wxTJo0qcvFeiSGE6cZmxCKaUOjYRbA85s5tJiIyF04HE5ee+013HfffZg3bx4GDx6M1atXw9fXF2vWrLnka0wmE2bPno1nnnkGvXv3vqKCPY41nDRUSlqGp3pi6kBoVUrsPlGO9JyOAzQREcmDQ+HEaDQiKysLKSkpF95AqURKSgoyMjIu+bpnn30WERERuOeeezp1nMbGRlRVVdndPJZvuOW+mnNyOENcqC/umZQIAPjzlhwYm80SV0RERJfjUDgpKyuDyWRCZGSk3fbIyEgUF7f/5frdd9/h3XffxTvvvNPp4yxfvhxBQUG2W1xcnCNlupeQXpb7ynxp6/BgD17bB+H+WuSW1eKfe/l7JiKSO6eO1qmursadd96Jd955B+Hh4Z1+3eLFi2EwGGy3wsJCJ1YpseB4y32Dgf1OnCRAr8EfpwwAAPz1659xvtYocUVERNQRtSM7h4eHQ6VSoaSkxG57SUkJoqKi2ux/8uRJ5OXlYfr06bZtZrOlWV2tVuPYsWPo06dPm9fpdDrodDpHSnNfWj/ALwKoLQXO51/og0Ld6rYxcXh/Tx6OFldj5dc/45mbh0hdEhERXYJDLSdarRajR49Genq6bZvZbEZ6ejqSk5Pb7D9w4ED8+OOPyM7Ott1uuukmXHfddcjOzvbsyzWOCEmw3J/Pk7IKj6ZSKrC0ZWjxh/sKcKK0WuKKiIjoUhxqOQGAtLQ0zJ07F2PGjMG4ceOwcuVK1NbWYt68eQCAOXPmIDY2FsuXL4der8eQIfb/hxocHAwAbbZ7tZBewOlMhhMnm9A3HL8cHIltR0rw/OYcrJ03TuqSiIioHQ6Hk5kzZ+LcuXNYunQpiouLMWLECGzdutXWSbagoABKJSeedQhbTlzmyRsHYcexUuw4dg7f/nwO1/TvIXVJRER0EYVwg5mpqqqqEBQUBIPBgMDAQKnL6X4HPwQ+XwD0+QVw50apq/F4z31xBO9+l4t+Ef7438OToFYxTBMROUNXv7/5V1kO2HLiUg/9oh9CfDU4XlqDjzMLpC6HiIguwnAiB8HWuU4KAbNJ2lq8QJCvBo/+sj8A4LVtP8NQ3yRxRURE1BrDiRwExgBKDWBuAqrOSl2NV7hjXDz6RfjjfF0TXk8/LnU5RETUCsOJHChVFyZj46Udl1CrlHhq2iAAwPsZecgtq5W4IiIismI4kQvrNPYMJy5z7YAIXNO/B5pMAsu35EhdDhERtWA4kQt2ipXE09MGQaVU4KsjJdhzskzqcoiICAwn8mENJ1wA0KX6RQZgdpLlktpzX+TAZJb9yHoiIo/HcCIXbDmRzCMp/RGoVyOnqAob9nvwIpNERG6C4UQugtnnRCqhflo8dH0/AMArX/2MmsZmiSsiIvJuDCdyYW05qT0HGDlyxNXmJCcgMdwPZTWNeHP7CanLISLyagwncuETDOiDLY/Ps9+Jq2nVSjx5o2Vo8T++y0VhRZ3EFREReS+GEznhcGJJpQyKwIQ+YTA2m7Fi61GpyyEi8loMJ3LCTrGSUigUWPKrwVAqgM2HivB9XoXUJREReSWGEznhcGLJDYoOxMyxcQAsqxebObSYiMjlGE7khC0nspD2ywHw16lx6LQBm7LPSF0OEZHXYTiREw4nloUeATosuK4vAGDF/47iXHWjxBUREXkXhhM5sbWc5AOClxOkNG9iAnqH+6G0uhHz1mZy7hMiIhdiOJGToDhAoQSa64GaUqmr8Wp6jQrv3jUWoX5aHD5ThQf/dQBNJrPUZREReQWGEzlRa4HAWMtjXtqRXGK4H9bcNRY+GhV2/nwOi/59CIItWkRETsdwIjfsFCsrI+KCsWr2SKiUCnx24Axe+eqY1CUREXk8hhO5sU7ExuHEsvGLgZF44ddDAACrtp/EPzPypC2IiMjDMZzIDVtOZGnm2Hg8mtIfALD0Pz9h6+FiiSsiIvJcDCdyE5xguWc4kZ2Hru+LWePiIQTw0LqDnEGWiMhJGE7kpvVwYpIVhUKB526+CimDImBsNuPe9/fjRGm11GUREXkchhO5sYaTqjNAMyf/khu1SonXZ43CyPhgGOqbMHfN9yipapC6LCIij8JwIjd+4YDGF4AAKgulroba4aNV4d25Y9E73A9nKusxd00mqhqapC6LiMhjMJzIjULBTrFuINRPi/fvHodwfx2OFlfj9x9kobHZJHVZREQegeFEjmyrE+dJWQVdRlyoL9bOGws/rQoZp8rx+IZDXMWYiKgbMJzIERcAdBtDYoOw+s7RUCsV+M8PZ7H8fzlSl0RE5PYYTuSIl3XcyqR+PfDSrcMAAO/sysU/dp2SuCIiIvfGcCJHHE7sdm4Z1ROLbhgIAHh+cw7++8NZiSsiInJfDCdy1LrlhAvNuY0HrumNucmWS3J//OQHZJwsl7giIiL3xHAiR8HxlvvGKqD+vLS1UKcpFAosnX4Vpg6JgtFkxv3/3I+jxVVSl0VE5HYYTuRI6wv4R1oes9+JW1EpFfjLzBEYlxCK6oZmzF2TiTOV9VKXRUTkVhhO5Mo2nJj9TtyNXqPCO3PGoF+EP0qqGnHXmkwY6jhJGxFRZzGcyBWHE7u1IF8N1t49DlGBehwvrcF9H+xHQxMnaSMi6gyGE7nicGK3Fxvsg7V3j0WATo3MvAo8uj4bJk7SRkR0WQwncsXhxB5hYFQg3p4zGlqVEv87XIxn//sTBEdgERF1iOFErthy4jEm9AnHq7cNBwC8n5GPt3dykjYioo4wnMhVSEufE0MhYGqWtha6YtOHx+DpaYMAACv+dxSfHTgtcUVERPLFcCJXAdGASguYm4GqM1JXQ93g3km9cd+kRADAnz49hF3Hz0lcERGRPDGcyJVSdWEyNg4n9hiLpw7C9OExaDYLPPDPLBw+Y5C6JCIi2WE4kTMOJ/Y4SqUCr/x2GJJ7h6HWaMK8td+jsKJO6rKIiGSF4UTO2CnWI+nUKrw9ZzQGRgXgXHUj5q7JREWtUeqyiIhkg+FEzjic2GMF6jV4/+5xiA32wamyWtzz/veoN3KSNiIigOFE3thy4tEiA/V4/+6xCPLR4GBBJf7w8UE0m8xSl0VEJDmGEzkLYZ8TT9c3IgD/mDsGOrUSX+eUYMnnnKSNiIjhRM6sLSd1ZUBjjaSlkPOMTQjFX28fCYUC+DizAK9/c0LqkoiIJMVwImf6IMAn1PK47Ji0tZBT3TAkCs/cdBUA4LVtP+OT7wslroiISDoMJ3LXc6zlPm+3tHWQ081JTsCD1/YBACz67BBe2JLDlYyJyCsxnMhd4iTLfd530tZBLvF46gDcNSEBQgB/33kKN73xHSdqIyKvw3Aidwkt4SR/D9fY8QIKhQL/d9NVeGfOGIT7a/FzSQ1+/eZurNp+giN5iMhrMJzIXdRQS98TYzVQ9IPU1ZCL/HJwJL58ZDJSr4pEk0ng5S+P4bdvZyC3rFbq0oiInI7hRO6UKqDX1ZbHeTulrYVcKsxfh9W/G43XbhuOAJ0aBwsqceNfd+GfGXkcbkxEHo3hxB1Y+53k7pK2DnI5hUKBW0b1xNZHJ2NCnzDUN5mw5POfMGdNJooNDVKXR0TkFAwn7iChpeWkYC9gapK2FpJEbLAPPrwnCcumD4ZOrcSu42WY8pdv8Xn2GalLIyLqdgwn7iDiKst8J021wNmDUldDElEqFZg3MRGbH5qEYT2DUNXQjIfXZWPhRwdwngsHEpEHYThxB0olkDDR8jiX/U68Xd8If/x7/gQ8ktIPKqUCXxwqQurKndh+rFTq0oiIugXDibtImGy5z2O/EwI0KiUeSemPjQ9OQJ8efiitbsS8977Hkxt/RG0jh5wTkXtjOHEX1k6xBfuA5kZpayHZGNYzGJsfmoS7JyYCAD7aV4Cpf92F/XkVEldGRNR1DCfuosdAwDccaK4HzmRJXQ3JiF6jwtLpg/HRvUmICdKjoKIOt72dgRe3HkVjM6e/JyL3w3DiLhSKC6N2OJU9tWNC33BsfXQyfjOqJ8wCeGvHSdz8xm7kFFVJXRoRkUMYTtyJbb4Tdoql9gXqNXj1tuFY/bvRCPXT4mhxNW5+YzdWf3sSJjMnbiMi99ClcLJq1SokJCRAr9cjKSkJmZmZl9z3nXfewaRJkxASEoKQkBCkpKR0uD91wNoptjATaOIEXHRpNwyJwpePTEbKoAgYTWas+N9RzHw7A/nlnP6eiOTP4XCyfv16pKWlYdmyZThw4ACGDx+O1NRUlJa2P4xxx44dmDVrFrZv346MjAzExcVhypQpOHOGk0c5LLwf4B8JmBqB099LXQ3JXI8AHd6ZMwYv3ToM/jo19uefx9S/7sJH+wo4/T0RyZpCOPhXKikpCWPHjsUbb7wBADCbzYiLi8Mf/vAHPPHEE5d9vclkQkhICN544w3MmTOnU8esqqpCUFAQDAYDAgMDHSnX83x6D3D4U+CaRcB1T0pdDbmJwoo6PLbhB+zLtYziuW5AD7z4m2GICNRLXBkRebKufn871HJiNBqRlZWFlJSUC2+gVCIlJQUZGRmdeo+6ujo0NTUhNDT0kvs0NjaiqqrK7kYt2CmWuiAu1Bcf3zceT08bBK1aie3HzmHKyp3YfKhI6tKIiNpwKJyUlZXBZDIhMjLSbntkZCSKi4s79R6LFi1CTEyMXcC52PLlyxEUFGS7xcXFOVKmZ0ts6Xdy+nugqV7aWsitKJUK3DupN774w9W4KiYQlXVNWPDRATy87iAMdVyziYjkw6WjdVasWIF169Zh48aN0Osv3Zy8ePFiGAwG262wsNCFVcpcaG8gIAYwGYHCfVJXQ26of2QANj44EX/4RV8oFcDn2WeRunIndh0/J3VpREQAHAwn4eHhUKlUKCkpsdteUlKCqKioDl/7yiuvYMWKFfjqq68wbNiwDvfV6XQIDAy0u1ELhaLVkGJOZU9do1Ur8ccpA/Dp/AlIDPdDcVUD7nw3E0s/P4w6I6e/JyJpORROtFotRo8ejfT0dNs2s9mM9PR0JCcnX/J1L730Ep577jls3boVY8aM6Xq1ZJHQEk64zg5doVHxIdj80NWYk9wLAPBBRj6m/e07HCg4L3FlROTNHL6sk5aWhnfeeQfvv/8+cnJyMH/+fNTW1mLevHkAgDlz5mDx4sW2/V988UUsWbIEa9asQUJCAoqLi1FcXIyampruOwtvY205OZMFNPL3SFfGV6vGszcPwQd3j0NUoB65ZbW49a09ePWrYzA2m6Uuj4i8kMPhZObMmXjllVewdOlSjBgxAtnZ2di6dautk2xBQQGKii6MAHjrrbdgNBpx6623Ijo62nZ75ZVXuu8svE1wLyAoDjA3s98JdZvJ/Xvgy0cmY8aIGJgF8Po3J/DrN3fj55JqqUsjIi/j8DwnUuA8J+3YOB/44SPg6keBlP+TuhryMJsPFeGpTT+isq4JWrUSj08ZgLuvToRKqZC6NCJyIy6Z54RkhJ1iyYmmDYvGV49MxnUDesDYbMaft+TgV69/h/XfF6DeyJWOici52HLirioLgZVDAIUKeCIf0AVIXRF5ICEE1n1fiOe+OIK6llASoFfj1tE98bvxvdCnh7/EFRKRnHX1+5vhxJ39dThwPg+4YwPQf4rU1ZAHq6g1YsP+Qny4Lx+FFRcm/5vYNwx3ju+FlEGRUKvYEEtE9rr6/a12Yk3kbAmTLOEkbxfDCTlVqJ8Wv7+mD+6b1BvfHj+Hf+3NR/rRUuw+UY7dJ8oRGajDrHHxmDUuHpFcr4eIrhBbTtzZD+uBjfcDMSOB+3dIXQ15mcKKOnycWYD13xeivNYIAFApFUi9KhK/G98Lyb3DoFCwAy2RN+NlHW9UdRZ4bRCgUAKL8gB9kNQVkRdqbDZh6+Fi/DMjH/vzL0ze1qeHH343vhduGdUTQT4aCSskIqkwnHirv40CKk4Cs9YBA6ZKXQ15uZyiKny4Nx8bD56xdaD10agwY2QMfje+F66KYYAm8iYMJ97qvw8DWWuB8QuAG16QuhoiAEB1QxM2HjyDf2bk43jphVmMR8UH487kXpg6JBp6jUrCConIFRhOvNWPnwL/vgeIGgY8wDlPSF6EEMjMrcA/9+Zj6+FiNJstf25CfDW4bWwcfpfUC3GhvhJXSUTOwnDiraqLgVcHAFAAfzoF+IZKXRFRu0qrG7A+sxAfZRagyNAAwLLI9rX9e+DO5F64pn8EZ6Al8jAMJ97szWSg9Agw5c/AhIVSV0PUoWaTGd8cLcU/9+Zj1/Ey2/bYYB/MHh+PmWPiEOavk7BCIuouDCfe7MAHwH/+APhFAI8cAjQ+UldE1Cm5ZbX41958bMg6DUN9EwBAq1LixqFRuDO5F0bFh3A4MpEbYzjxZs1G4PXRgKEAuOFFYPwDUldE5JCGJhP++8NZfLg3Hz+cNti2D4oOxO/Gx2PGiFj46ThnJJG7YTjxdvvXAF88CgREAw9lAxrO0knu6YfCSny4Nx//+eEsGpvNAAB/nRq/GRWL343vhX6RXEeKyF0wnHi75kbgbyOBqjPAja8A4+6TuiKiK1JZZ8SnWafx4d585JXX2baP7x2KO8cnYMpVkdBwPR8iWWM4ISDzHWDLY0BgLPDQQUDNToXk/sxmgd0ny/DPjHx8nVOCltHI6BGgw6yxcZiVFI/oIPazIpIjhhMCmhqAv40AqouAX60ExsyTuiKibnW2sh4fZxbg48xClNU0ArCs55MyKAJ3jk/AhD5hUHI4MpFsMJyQxd7VwNZFQFA88NABQMU1TcjzGJvN+OqIZT2ffbkVtu2J4X6YnRSP346OQ5Av/9snkhrDCVk01QMrhwG1pcBNbwCj7pS6IiKn+rmkGh/uzcdnB86gprEZAKDXKHHT8BjcODQaYxJC4c+RPkSSYDihC/a8AXz1FBCSACzcz9YT8go1jc34PNuyns/R4mrbdpVSgatiApGUGIpxiWEYlxDKVhUiF2E4oQuMtZbWk7oyYMZbwIg7pK6IyGWEEMjKP49P9hdiz8lynD5fb/e8QgEMiAzA+N5hGJcYinGJoQjnjLRETsFwQva+Wwl8vQwI7QMsyARUbNYm73Smsh7f51ZgX2459uVW4NS52jb79Onhh6TeYUhKDEVSYhiigjhPEFF3YDghe401wMqhQH0FcMs7wLDbpK6ISBZKqxuQmVuBzNwK7DtVgWMl1W32iQ/1bbkMFIrxvcPQM8SH0+gTdQHDCbW18xXgm+eA8P7Ag3sBpUrqiohk53ytEZl5FbbA8tNZg20uFavoIL2tz0pS71D0DvdjWCHqBIYTaquhClg5BGgwALeuAYb8RuqKiGSvqqEJWfnnse9UBTJzy3HotAHNF6WVcH+drWUlqXco+kcEcH4VonYwnFD7drwI7HgB6DEImL8HUHK6byJH1BmbcbCgEvtOWfqsHCyshLFlzR+rYF8NxiaE2vqsDIoOgJpT6xMxnNAl1Fda+p40VgG3fQAMvlnqiojcWkOTCYdOG7DvVDky8yqwP+886ptMdvv469QY3SsESb0tgWVobDC0aoYV8j4MJ3Rp3/wZ2PkSEDkE+P0utp4QdaMmkxmHzxiwr6XPyve5FahumQzOSq9RYlR8CJISLcOXR8YHQ69hHzDyfAwndGl1FZZ5T4zVwO0fAQOnSV0RkccymQVyiqpawko5MnMrcL6uyW4frUqJ4XFBtrAyulcI/DiLLXkghhPqWPqzwK5XgejhwP3fWmaiIiKnM5sFTpyrsfVZ2ZdbgXPVjXb7qJQKDIkNwviWTrZjEkIR5MNZbMn9MZxQx2rLLX1PmmqBOz4B+qdKXRGRVxJCIK+8ztJnpSWsnKlsO4vtoKhAW5+VcYlhCPXTSlQxUdcxnNDlbVsK7P4rEDsauDedrSdEMnH6fF3L0GXLTLZ55XVt9ukX4Y+k3i1zrSSGIiJAx7lWSPYYTujyas5ZWk+a64FfPgdMfEjqioioHSVVDbY+K/tOVeB4aU2bfQL1aiSE+yEhzA8J4X5IDPdFrzA/JIb5IYStLCQTDCfUObv/BmxbYnk89WUg6X5p6yGiyyqvacT3eeexr6WD7ZGiKnT0lzvIR2MJLGEtgSW8JcCE+XFFZnIphhPqHCEsU9rvetXy869WAmPmSVoSETmmocmE/PI65JbVIr+8Fnnltcgtq0VeWR2Kqxo6fG2Ir+ZCYAnzQ0K4r631hZ1wqbsxnFDnCQF89TSQ8QYABTDjTWDEHVJXRUTdoN5oQn5FLfLKapFbVme5L7eEmJKqxg5fG+qnRUKYr/3lopYAE6BncCHHMZyQY4QA/rcIyHwbUCgtKxcPvVXqqojIieqMzcgrq0NeS2tLXktrS255bZvhzRcL99eiV5gltCSG2wcYf87RQpfAcEKOEwL44hEgay2gUAG/fY/T2xN5qZrGZsslopbwYrlMVIu88jqU1VwuuOgsgaUlrLS+XMTJ5bwbwwl1jdkM/GchkP0vQKkGZn4IDJgqdVVEJCPVDU12fVxyra0vZbUorzV2+NoeATrbpSHrZaJeLT/7ahlcPB3DCXWd2QRs/D3w4wZApQVu/xjolyJ1VUTkBqoampDfcmkor6W1xdLHpQ4VlwkukYG6lstEfi2ddC0BpleoH3y0XHvIEzCc0JUxNQP/vhs48jmg1gN3rAd6Xyt1VUTkxgx1Ta36t7S6XFRei8qL1hu6WFSgHgnhvq1GFVnue4X5ctFEN8JwQlfO1AR8Mgc4tgXQ+AK/+zfQa4LUVRGRB6qsMyKvvGU0UdmFDrq5ZbWoami+5OsUCiA6UG8JK+F+6BXqix4BOoT6aRHur0OYvxahflro1AwwcsBwQt2juRFYNxs4sQ3Q+gN3bgTixkldFRF5kfO1xguXiVoCjLXVpbqD4NJagE5tCyph/jqE+WkR5q9FmJ/Odm8JNFqE+GmhUSmdfFbeieGEuk9TPfDx7cCpHYAuEJjzORA7SuqqiMjLCSFQUWu0CyyFFXUorzWirMaI8ppGVNQa0Wx2/GstyEfTElpaBxhLsAltCTbhLY9DfLVQKbmuUWcwnFD3MtYB/7oVyN8N6IOBuf8FoodJXRURUYeEEKiqb0Z5bSPKay2BxXJvREWtEWU1jbbH5bWWMONollEogFBfrS202LXG+GsR7nehxSbcX4tAvQZKLw0zDCfU/RqrgX/eApzOBHzDgDs3MaAQkUcxmQUM9U12Iaa8ttF2X3FRq8z5y3TkbY9KqbCElVZhxnpJKbQl2LR+HKBTe8yK0wwn5BwNBuCDGcDZA5aZZIfPAq75ExCSIHVlREQu12wyo6KupeWlxnihdabVY0urjKWVprN9ZFrTqpS2VhlrR1/bpSU/XZsWG1+tSrZhhuGEnKf+PLBpAXBss+VnpRoYeScw+XEgKFba2oiIZKyx2YTztU0os4WWtkGmzHqZqaYRtUaTw8fQa5R2/WRC/XQtLTEtnYFb9Z8J89O6dCg2wwk53+n9wDfPA6e2W35W6SwrGl+dBgRESlsbEZEHaGgyXbI1pqzVpabyGkvLTGOz2eFj+GlVto6+rUPMrLHxiA/z7dbzYTgh18nbDWz/s6WzLACofYCk+4EJDwN+YdLWRkTkJYQQqDOa2vST6agDcJPp0l/5/54/AaN7hXRrjQwn5FpCWIYaf/M8cGa/ZZvWHxg/H0heCPgES1kdERFdRAiB6sZmS4hp3QG45fGD1/ZBRKC+W4/JcELSEAI4/pUlpBQfsmzTBwHJfwDGPwDoAqStj4iIJNPV729OiUdXRqEA+qcCv98J3PZPoMcgywif7c8DK4cBu/9mmTOFiIiokxhOqHsoFMDgm4D5u4HfvAuE9QXqK4BtS4C/jQD2rgaaGqSukoiI3AAv65BzmJqBQ+uBb1cAlQWWbf6RQJ/rgcRJQMLVQHC8tDUSEZFTsc8JyVOzEcj+EPj2ZaD6rP1zwb1agspkS1jhnClERB6F4YTkrbkRyNsF5H0H5O4Czh4ExEWTDYX2BhImWW6Jk4CAKGlqJSKibsFwQu6loQoo3Afk7rSElqIfAHHRZEJh/S5cAkqYBPhHSFMrERF1CcMJubcGA5Cf0dK6sgsoOgTgov80ewxsaVlpCSuc8I2ISNYYTsiz1J8H8vdcuAxU8mPbfXoMtIwKCooDguMsHWyDWu59QiwjiIiISDIMJ+TZ6ios0+XntrSslB7peH+tfzuhJc7SCTcoznKJiOGFiMipGE7Iu9SWAWcOAJX5lqHKhkLLfWUhUFt6+derdJawYhdg4i88DogGlK5buZOIyBN19ftb7cSaiJzHLxzoP6X955rqAcPptqHF+rjqLGBqBMpPWG7tUaqBwBjAPwrQBwK6wIvugy6xveWm4j8tIqKu6tJf0FWrVuHll19GcXExhg8fjtdffx3jxo275P4bNmzAkiVLkJeXh379+uHFF1/EjTfe2OWiiTqk8QHC+1lu7Wk2AlVnWsJK4UUhpsDynLn5ws9dqsG3/eBiuw9q9XNA2+e0fpbWHSUncSYi7+NwOFm/fj3S0tKwevVqJCUlYeXKlUhNTcWxY8cQEdF2qOeePXswa9YsLF++HL/61a/w0UcfYcaMGThw4ACGDBnSLSdB5BC1FghNtNzaYzYB1cWWYFJ7Dmissgx9tt0bWu6r2z7XXG95j6Y6y62m+MpqVaotIUWtBdR6QKUF1LoL29q917Xar6P9HdxPpQEUKsvlLvbXISIncrjPSVJSEsaOHYs33ngDAGA2mxEXF4c//OEPeOKJJ9rsP3PmTNTW1uKLL76wbRs/fjxGjBiB1atXd+qY7HNCbqPZ2BJaDBcFmup2ws3FoadlP2ON1GdxeQrlhaCiVLc87mibuiXUtGyzPa9q9RpVx9uUastxbdus79neNuuxW25QtDxWtPOz4jLPW3+2PkYn9u3qeysuvO5y742WgNg6KLbZ7ozHlzoWLtTfnY8ZhN2aS/qcGI1GZGVlYfHixbZtSqUSKSkpyMjIaPc1GRkZSEtLs9uWmpqKTZs2XfI4jY2NaGxstP1sMBgAWE6SSP40gDoc8A8H/LvwcrMJMNYCpibLzLqmRstjU6Ml/JiMlu1m6/PGlu0tj1v/3GwEzEb711nfy9QINF/0vibjRccxAjC1U6TpEtuJnKlVILILM613uXjbxYHqcts7+ZoO369NUVf4fh3V15njdLKG3/wDiB6G7mT93nZ07I1D4aSsrAwmkwmRkZF22yMjI3H06NF2X1NcXNzu/sXFl27uXr58OZ555pk22+Pi4hwpl4iIiDrrqUlOe+vq6moEBQV1en9ZDilYvHixXWuL2WxGRUUFwsLCoOjGJr6qqirExcWhsLDQYy8X8Rw9A8/RM/AcPQPPsfOEEKiurkZMTIxDr3MonISHh0OlUqGkpMRue0lJCaKi2l+kLSoqyqH9AUCn00Gn09ltCw4OdqRUhwQGBnrsf2BWPEfPwHP0DDxHz8Bz7BxHWkysHBqnqNVqMXr0aKSnp9u2mc1mpKenIzk5ud3XJCcn2+0PANu2bbvk/kREROTdHL6sk5aWhrlz52LMmDEYN24cVq5cidraWsybNw8AMGfOHMTGxmL58uUAgIcffhjXXHMNXn31VUybNg3r1q3D/v378fe//717z4SIiIg8gsPhZObMmTh37hyWLl2K4uJijBgxAlu3brV1ei0oKICy1cRREyZMwEcffYSnn34aTz75JPr164dNmzbJYo4TnU6HZcuWtbmE5El4jp6B5+gZeI6egefofG6xtg4RERF5D86NTURERLLCcEJERESywnBCREREssJwQkRERLLi1eFk1apVSEhIgF6vR1JSEjIzM6UuqV3Lly/H2LFjERAQgIiICMyYMQPHjh2z2+faa6+FQqGwuz3wwAN2+xQUFGDatGnw9fVFREQEHn/8cTQ3N9vts2PHDowaNQo6nQ59+/bF2rVrnX16AID/+7//a1P/wIEDbc83NDRgwYIFCAsLg7+/P37zm9+0mdxPzucHAAkJCW3OUaFQYMGCBQDc8zPcuXMnpk+fjpiYGCgUijZrZgkhsHTpUkRHR8PHxwcpKSk4fvy43T4VFRWYPXs2AgMDERwcjHvuuQc1NfaLHx46dAiTJk2CXq9HXFwcXnrppTa1bNiwAQMHDoRer8fQoUOxZcsWp59jU1MTFi1ahKFDh8LPzw8xMTGYM2cOzp49a/ce7X32K1ascItzBIC77rqrTf033HCD3T7u/DkCaPffpkKhwMsvv2zbR86fY2e+J1z5d/SKv1+Fl1q3bp3QarVizZo14qeffhL33XefCA4OFiUlJVKX1kZqaqp47733xOHDh0V2dra48cYbRXx8vKipqbHtc80114j77rtPFBUV2W4Gg8H2fHNzsxgyZIhISUkRBw8eFFu2bBHh4eFi8eLFtn1OnTolfH19RVpamjhy5Ih4/fXXhUqlElu3bnX6OS5btkxcddVVdvWfO3fO9vwDDzwg4uLiRHp6uti/f78YP368mDBhgtucnxBClJaW2p3ftm3bBACxfft2IYR7foZbtmwRTz31lPjss88EALFx40a751esWCGCgoLEpk2bxA8//CBuuukmkZiYKOrr62373HDDDWL48OFi7969YteuXaJv375i1qxZtucNBoOIjIwUs2fPFocPHxYff/yx8PHxEW+//bZtn927dwuVSiVeeuklceTIEfH0008LjUYjfvzxR6eeY2VlpUhJSRHr168XR48eFRkZGWLcuHFi9OjRdu/Rq1cv8eyzz9p9tq3//cr5HIUQYu7cueKGG26wq7+iosJuH3f+HIUQdudWVFQk1qxZIxQKhTh58qRtHzl/jp35nnDV39Hu+H712nAybtw4sWDBAtvPJpNJxMTEiOXLl0tYVeeUlpYKAOLbb7+1bbvmmmvEww8/fMnXbNmyRSiVSlFcXGzb9tZbb4nAwEDR2NgohBDiT3/6k7jqqqvsXjdz5kyRmpravSfQjmXLlonhw4e3+1xlZaXQaDRiw4YNtm05OTkCgMjIyBBCyP/82vPwww+LPn36CLPZLIRw/8/w4j/4ZrNZREVFiZdfftm2rbKyUuh0OvHxxx8LIYQ4cuSIACC+//572z7/+9//hEKhEGfOnBFCCPHmm2+KkJAQ2zkKIcSiRYvEgAEDbD/fdtttYtq0aXb1JCUlid///vdOPcf2ZGZmCgAiPz/ftq1Xr17iL3/5yyVfI/dznDt3rrj55psv+RpP/Bxvvvlm8Ytf/MJumzt9jhd/T7jy72h3fL965WUdo9GIrKwspKSk2LYplUqkpKQgIyNDwso6x2AwAABCQ0Pttv/rX/9CeHg4hgwZgsWLF6Ours72XEZGBoYOHWq3QnRqaiqqqqrw008/2fZp/Tux7uOq38nx48cRExOD3r17Y/bs2SgoKAAAZGVloampya62gQMHIj4+3labO5xfa0ajER9++CHuvvtuu8Us3f0zbC03NxfFxcV29QQFBSEpKcnucwsODsaYMWNs+6SkpECpVGLfvn22fSZPngytVmvbJzU1FceOHcP58+dt+8jlvA0GAxQKRZv1wFasWIGwsDCMHDkSL7/8sl1TuTuc444dOxAREYEBAwZg/vz5KC8vt6vfkz7HkpISbN68Gffcc0+b59zlc7z4e8JVf0e76/tVlqsSO1tZWRlMJpPdBwAAkZGROHr0qERVdY7ZbMYjjzyCiRMn2s2ye8cdd6BXr16IiYnBoUOHsGjRIhw7dgyfffYZAKC4uLjd87U+19E+VVVVqK+vh4+Pj9POKykpCWvXrsWAAQNQVFSEZ555BpMmTcLhw4dRXFwMrVbb5o99ZGTkZWu3PtfRPq44v4tt2rQJlZWVuOuuu2zb3P0zvJi1pvbqaV1vRESE3fNqtRqhoaF2+yQmJrZ5D+tzISEhlzxv63u4SkNDAxYtWoRZs2bZLZb20EMPYdSoUQgNDcWePXuwePFiFBUV4bXXXrOdh5zP8YYbbsAtt9yCxMREnDx5Ek8++SSmTp2KjIwMqFQqj/sc33//fQQEBOCWW26x2+4un2N73xOu+jt6/vz5bvl+9cpw4s4WLFiAw4cP47vvvrPbfv/999seDx06FNHR0bj++utx8uRJ9OnTx9VlOmzq1Km2x8OGDUNSUhJ69eqFTz75xKVfqK7y7rvvYurUqXbLiLv7Z+jtmpqacNttt0EIgbfeesvuubS0NNvjYcOGQavV4ve//z2WL1/uFlOg33777bbHQ4cOxbBhw9CnTx/s2LED119/vYSVOceaNWswe/Zs6PV6u+3u8jle6nvCnXjlZZ3w8HCoVKo2vZRLSkoQFRUlUVWXt3DhQnzxxRfYvn07evbs2eG+SUlJAIATJ04AAKKioto9X+tzHe0TGBjo8oAQHByM/v3748SJE4iKioLRaERlZWWb2i5Xu/W5jvZx9fnl5+fj66+/xr333tvhfu7+GVpr6ujfWVRUFEpLS+2eb25uRkVFRbd8tq7692wNJvn5+di2bdtll5hPSkpCc3Mz8vLyALjHObbWu3dvhIeH2/236QmfIwDs2rULx44du+y/T0Cen+Olvidc9Xe0u75fvTKcaLVajB49Gunp6bZtZrMZ6enpSE5OlrCy9gkhsHDhQmzcuBHffPNNm2bD9mRnZwMAoqOjAQDJycn48ccf7f6AWP+IDh482LZP69+JdR8pfic1NTU4efIkoqOjMXr0aGg0Grvajh07hoKCAltt7nR+7733HiIiIjBt2rQO93P3zzAxMRFRUVF29VRVVWHfvn12n1tlZSWysrJs+3zzzTcwm822cJacnIydO3eiqanJts+2bdswYMAAhISE2PaR6rytweT48eP4+uuvERYWdtnXZGdnQ6lU2i6FyP0cL3b69GmUl5fb/bfp7p+j1bvvvovRo0dj+PDhl91XTp/j5b4nXPV3tNu+XzvdddbDrFu3Tuh0OrF27Vpx5MgRcf/994vg4GC7XspyMX/+fBEUFCR27NhhN4Strq5OCCHEiRMnxLPPPiv2798vcnNzxeeffy569+4tJk+ebHsP6xCxKVOmiOzsbLF161bRo0ePdoeIPf744yInJ0esWrXKZUNt//jHP4odO3aI3NxcsXv3bpGSkiLCw8NFaWmpEMIyBC4+Pl588803Yv/+/SI5OVkkJye7zflZmUwmER8fLxYtWmS33V0/w+rqanHw4EFx8OBBAUC89tpr4uDBg7aRKitWrBDBwcHi888/F4cOHRI333xzu0OJR44cKfbt2ye+++470a9fP7shqJWVlSIyMlLceeed4vDhw2LdunXC19e3zfBMtVotXnnlFZGTkyOWLVvWbUNQOzpHo9EobrrpJtGzZ0+RnZ1t9+/TOrphz5494i9/+YvIzs4WJ0+eFB9++KHo0aOHmDNnjlucY3V1tXjsscdERkaGyM3NFV9//bUYNWqU6Nevn2hoaLC9hzt/jlYGg0H4+vqKt956q83r5f45Xu57QgjX/R3tju9Xrw0nQgjx+uuvi/j4eKHVasW4cePE3r17pS6pXQDavb333ntCCCEKCgrE5MmTRWhoqNDpdKJv377i8ccft5sjQwgh8vLyxNSpU4WPj48IDw8Xf/zjH0VTU5PdPtu3bxcjRowQWq1W9O7d23YMZ5s5c6aIjo4WWq1WxMbGipkzZ4oTJ07Ynq+vrxcPPvigCAkJEb6+vuLXv/61KCoqsnsPOZ+f1ZdffikAiGPHjtltd9fPcPv27e3+tzl37lwhhGU48ZIlS0RkZKTQ6XTi+uuvb3Pu5eXlYtasWcLf318EBgaKefPmierqart9fvjhB3H11VcLnU4nYmNjxYoVK9rU8sknn4j+/fsLrVYrrrrqKrF582ann2Nubu4l/31a56/JysoSSUlJIigoSOj1ejFo0CDxwgsv2H2xy/kc6+rqxJQpU0SPHj2ERqMRvXr1Evfdd1+bLxp3/hyt3n77beHj4yMqKyvbvF7un+PlvieEcO3f0Sv9flW0nBQRERGRLHhlnxMiIiKSL4YTIiIikhWGEyIiIpIVhhMiIiKSFYYTIiIikhWGEyIiIpIVhhMiIiKSFYYTIiIikhWGEyIiIpIVhhMiIiKSFYYTIiIikhWGEyIiIpKV/wcntzcEI2NFFwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(l,cl_eg_1h,label='1h')\n",
    "plt.plot(l,cl_eg_2h,label='2h')\n",
    "plt.legend()\n",
    "# plt.yscale('log')\n",
    "# plt.yscale('log')\n",
    "plt.ylim(0,1e-10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "c1e730d8-fa1a-40e7-867b-29fc00bc5ef7",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\"cleg_l_1h_2h_1p2h_fb_times_nfw.txt\",np.c_[l,cl_eg_1h,cl_eg_2h,cl_eg_1h+cl_eg_2h])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e819006-c075-420a-aad6-2bed6d555cc3",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}