significance_of_abundances.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "from petitRADTRANS import Radtrans\n",
    "import petitRADTRANS.retrieval_examples.emission.master_retrieval_model as rm\n",
    "from petitRADTRANS import nat_cst as nc\n",
    "import petitRADTRANS.rebin_give_width as rgw\n",
    "from scipy.interpolate import interp1d\n",
    "import sklearn\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import pandas as pd\n",
    "import os\n",
    "import sys\n",
    "\n",
    "from torch.distributions import Independent, Distribution\n",
    "from fab_priors import BoxUniform_New\n",
    "\n",
    "import torch\n",
    "from sbi.inference import SNRE_A, SNRE, prepare_for_sbi, simulate_for_sbi\n",
    "from sbi.utils.get_nn_models import posterior_nn\n",
    "from sbi import utils as utils\n",
    "from sbi.types import Array, OneOrMore, ScalarFloat\n",
    "\n",
    "from Coverage_eval_SBI import Coverage_class\n",
    "\n",
    "np.set_printoptions(threshold=sys.maxsize)\n",
    "torch.set_printoptions(threshold=10_000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "op= '/home/mvasist/petitRADTRANS/petitRADTRANS/retrieval_examples/emission/'\n",
    "observation_files = {}\n",
    "observation_files['NIRISS SOSS'] = op +'NIRISS_SOSS_flux.dat'\n",
    "observation_files['NIRSpec G395M'] = op +'NIRSpec_G395M_flux.dat'\n",
    "observation_files['MIRI LRS'] = op +'MIRI_LRS_flux.dat'\n",
    "\n",
    "# Wavelength range of observations, fixed parameters that will not be retrieved\n",
    "WLENGTH = [0.3, 15.0] \n",
    "# LOG_G =  2.58\n",
    "R_pl =   1.84*nc.r_jup_mean\n",
    "R_star = 1.81*nc.r_sun\n",
    "gamma = 1\n",
    "t_equ= 0\n",
    "\n",
    "# Get host star spectrum to calculate F_pl / F_star later.\n",
    "T_star = 6295.\n",
    "x = nc.get_PHOENIX_spec(T_star)\n",
    "fstar = interp1d(x[:,0], x[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NIRISS SOSS\n",
      "NIRSpec G395M\n",
      "MIRI LRS\n"
     ]
    }
   ],
   "source": [
    "####################################################################################\n",
    "####################################################################################\n",
    "### READ IN OBSERVATION\n",
    "####################################################################################\n",
    "####################################################################################\n",
    "\n",
    "# Read in data, convert all to cgs! \n",
    "\n",
    "'''\n",
    "Im using only data_flux_nu_error['MIRI LRS'] from here to calculate the likelihood. \n",
    "'''\n",
    "\n",
    "data_wlen = {}\n",
    "data_flux_nu = {}\n",
    "data_flux_nu_error = {}\n",
    "data_wlen_bins = {}\n",
    "\n",
    "for name in observation_files.keys():\n",
    "    print(name)\n",
    "    dat_obs = np.genfromtxt(observation_files[name])\n",
    "    data_wlen[name] = dat_obs[:,0]*1e-4\n",
    "    data_flux_nu[name] = dat_obs[:,1]\n",
    "    data_flux_nu_error[name] = dat_obs[:,2]\n",
    "    \n",
    "    data_wlen_bins[name] = np.zeros_like(data_wlen[name])\n",
    "    data_wlen_bins[name][:-1] = np.diff(data_wlen[name])\n",
    "    data_wlen_bins[name][-1] = data_wlen_bins[name][-2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Simulator(params): \n",
    "\n",
    "    NaN_spectra = 0\n",
    "\n",
    "    atmosphere = Radtrans(line_species = ['H2O', 'CO_all_iso', \\\n",
    "                                         'CO2', 'CH4', \\\n",
    "                                          'Na', 'K'], \\\n",
    "          rayleigh_species = ['H2', 'He'], \\\n",
    "          continuum_opacities = ['H2-H2', 'H2-He'], \\\n",
    "          wlen_bords_micron = WLENGTH)#, mode='c-k')\n",
    "\n",
    "\n",
    "    pressures = np.logspace(-6, 2, 100)\n",
    "    atmosphere.setup_opa_structure(pressures)\n",
    "    temperature = 1200. * np.ones_like(pressures)\n",
    "\n",
    "    \n",
    "    t_int = params[0].numpy()                             #200.\n",
    "    log_kappa_IR = params[1].numpy()                      #-2\n",
    "    log_gravity = params[2].numpy()                       #params[5].numpy() 1e1**2.45 \n",
    "\n",
    "    gravity = np.exp(log_gravity)\n",
    "    kappa_IR = np.exp(log_kappa_IR)\n",
    "    \n",
    "    temperature = nc.guillot_global(pressures, kappa_IR, gamma, gravity, t_int, t_equ)\n",
    "    \n",
    "    abundances = {}\n",
    "    abundances['H2'] = np.exp(params[3].numpy())*np.ones_like(temperature)         #0.74 * np.ones_like(temperature) (params[3].numpy()) 0.75 * np.ones_like(temperature)\n",
    "    abundances['He'] = np.exp(params[4].numpy())*np.ones_like(temperature)         #0.25 * np.ones_like(temperature)  #0.24 * np.ones_like(temperature) (params[4].numpy())\n",
    "    abundances['H2O'] = np.exp(params[5].numpy())*np.ones_like(temperature)        #0.001 * np.ones_like(temperature)\n",
    "    abundances['CO_all_iso'] = np.exp(params[6].numpy())*np.ones_like(temperature) #0.01 * np.ones_like(temperature)\n",
    "    abundances['CO2'] = np.exp(params[7].numpy())*np.ones_like(temperature)        #0.00001 * np.ones_like(temperature)\n",
    "    abundances['CH4'] = np.exp(params[8].numpy())*np.ones_like(temperature)        #0.000001 * np.ones_like(temperature)\n",
    "    abundances['Na'] = np.exp(params[9].numpy())*np.ones_like(temperature)         #0.00001 * np.ones_like(temperature)\n",
    "    abundances['K'] = np.exp(params[10].numpy())*np.ones_like(temperature)          #0.000001 * np.ones_like(temperature)\n",
    "\n",
    "    MMW = rm.calc_MMW(abundances) * np.ones_like(temperature)\n",
    "    \n",
    "    atmosphere.calc_flux(temperature, abundances, gravity, MMW)\n",
    "\n",
    "    wlen, flux_nu = nc.c/atmosphere.freq, atmosphere.flux/1e-6\n",
    "\n",
    "    # Just to make sure that a long chain does not die\n",
    "    # unexpectedly:\n",
    "    # Return -inf if forward model returns NaN values\n",
    "    if np.sum(np.isnan(flux_nu)) > 0:\n",
    "        print(\"NaN spectrum encountered\")\n",
    "        NaN_spectra += 1\n",
    "        return torch.ones([1,371])* -np.inf\n",
    "\n",
    "    # Convert to observation for emission case\n",
    "    flux_star = fstar(wlen)\n",
    "    flux_sq   = flux_nu/flux_star*(R_pl/R_star)**2 \n",
    "\n",
    "    flux_rebinned = rgw.rebin_give_width(wlen, flux_sq, \\\n",
    "                data_wlen['MIRI LRS'], data_wlen_bins['MIRI LRS'])\n",
    "\n",
    "    #flux_rebinned = np.reshape(flux_rebinned, (371,1))    \n",
    "\n",
    "    FR= torch.Tensor(flux_rebinned)\n",
    "    \n",
    "    p1, t1 = pressures, temperature\n",
    "    \n",
    "    return FR ,wlen, flux_nu, p1, t1   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "abund = ['He', 'H2', 'H20', 'CO_all_iso', 'CO2', 'CH4', 'Na', 'K']\n",
    "arr = np.zeros((8,8))\n",
    "# abundances_values = [np.log(0.75), np.log(0.25), -3, -2, -5, -6, -5, -6] \n",
    "# np.fill_diagonal(arr, abundances_values)\n",
    "\n",
    "abundances_values_0 = [0, 0, 0, 0, 0, 0, 0, 0] \n",
    "np.fill_diagonal(arr, abundances_values_0)\n",
    "\n",
    "abundances = arr\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "\n",
      "  Read CIA opacities for H2-H2...\n",
      "  Read CIA opacities for H2-He...\n",
      " Done.\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAGtCAYAAAAs+MGvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABXKElEQVR4nO3dd3yV9fn/8deVDYQkQAKyhyAIgmxUwIELt63Wbetoba32q/5qq3ZYO6zaYd1aV12toy607oEigrKRPZQNMgKBQAhZn98f55yQhJPk3CfnnDvj/fTBg/ucc5/7viJwcuXzuT7Xx5xziIiIiEh8JfkdgIiIiEhLoKRLREREJAGUdImIiIgkgJIuERERkQRQ0iUiIiKSACl+B1Cf3Nxc16tXL7/DEBEREanX7Nmztznn8sK91uiTrl69ejFr1iy/wxARERGpl5mtqe01TS+KiIiIJICSLhEREZEEUNIlIiIikgCNvqZLREREmo/S0lLWr19PcXGx36E0SEZGBt26dSM1NTXi9yjpEhERkYRZv349bdu2pVevXpiZ3+FExTlHfn4+69evp3fv3hG/T9OLIiIikjDFxcV06NChySZcAGZGhw4dPI/WKekSERGRhGrKCVdINF+Dki4RERGRBFDSJSIiIi3Ku+++S//+/enbty933nnnAa/v27eP888/n759+zJmzBhWr14dk/sq6RIREZEWo7y8nGuuuYZ33nmHxYsX8/zzz7N48eJq5zzxxBO0a9eOlStXcsMNN3DTTTfF5N5KukRERKTFmDFjBn379qVPnz6kpaVxwQUXMGnSpGrnTJo0iR/84AcAnHvuuXz00Uc45xp8b7WMEBEREV/8/s1FLN64K6bXHNgli9+dMajW1zds2ED37t0rH3fr1o0vv/yy1nNSUlLIzs4mPz+f3NzcBsWmkS4RERGRBEj4SJeZ9QderPJUH+BW59w9iY5FRERE/FPXiFS8dO3alXXr1lU+Xr9+PV27dg17Trdu3SgrK2Pnzp106NChwfdO+EiXc26Zc26oc24oMAIoAl5LdBwiIiLS8owaNYoVK1awatUqSkpKeOGFFzjzzDOrnXPmmWfy9NNPA/Dyyy8zYcKEmPQW83t68Xjga+fcGp/jaJbeXbiJXje/xe59ZX6HIiIi0iikpKTwwAMPcPLJJ3PooYdy3nnnMWjQIG699VbeeOMNAK688kry8/Pp27cvd999d9i2ElHdOyZXid4FwPM1nzSzq4CrAHr06JHomJqN+z5aCcDqbXs4rGu2z9GIiIg0Dqeeeiqnnnpqtef+8Ic/VB5nZGTw3//+N+b39W2ky8zSgDOBA74q59yjzrmRzrmReXl5iQ+umUhOCgyFllc0fJmriIiINIyf04unAHOcc5t9jKFZSwolXTHoLSIiIiIN42fSdSFhphYldpKDNX8a6RIRkcYkFo1G/RbN1+BL0mVmbYATgVf9uL+IiIj4IyMjg/z8/CadeDnnyM/PJyMjw9P7fCmkd87tARre8ELqFBrg0kiXiIg0Ft26dWP9+vVs3brV71AaJCMjg27dunl6j9+rFyUBlHSJiEhjkZqaSu/evf0Owxd+9+mSONLqRRERkcZDSVczlhQqpG/C8+YiIiLNhZKuZiwpuGVBebmSLhEREb8p6WrGKpMujXSJiIj4TklXM6aaLhERkcZDSVczZmqOKiIi0mgo6WrGQtOLFZpeFBER8Z2SrhZAOZeIiIj/lHS1AA5lXSIiIn5T0tUCaKRLRETEf0q6WgAlXSIiIv5T0tUCKOcSERHxn5KuFsBpqEtERMR3SrpaAKVcIiIi/lPS1RIo6xIREfGdkq4WQC0jRERE/KekqwVQSZeIiIj/lHSJiIiIJICSrhZAA10iIiL+U9LVAmh6UURExH9KuloAFdKLiIj4T0lXC6CRLhEREf8p6WoBlHOJiIj4T0lXS6ChLhEREd8p6WoBlHKJiIj4T0lXC6CBLhEREf8p6WoBnLIuERER3ynpagF27yvzOwQREZEWT0lXC/C395f7HYKIiEiLp6RLREREJAGUdImIiIgkgJIuERERkQRQ0iUiIiKSAEq6mjE1ihAREWk8lHSJiIiIJICSLhEREZEEUNLVzCzauJMVmwv9DkNERERqSPE7AImt0+6bCsDqO0/DfI5FRERE9tNIl4iIiEgC+JJ0mVmOmb1sZkvNbImZHelHHCIiIiKJ4tf04r3Au865c80sDWjtUxzNmlpGiIiINB4JT7rMLBs4GrgMwDlXApQkOo6WwDmlXSIiIo2FH9OLvYGtwL/MbK6ZPW5mbXyIo9kzUym9iIhIY+FH0pUCDAceds4NA/YAN1c9wcyuMrNZZjZr69atPoTYPKQkBZKu8f1yfY5ERERE/Ei61gPrnXNfBh+/TCAJq+Sce9Q5N9I5NzIvLy/hATYXoenF9JRknyMRERGRhCddzrlvgXVm1j/41PHA4kTH0RK4MEciIiLiD79WL/4M+Hdw5eI3wOU+xdGsheroVU8vIiLiP1+SLufcPGCkH/duzmquVnQ1fhcRERH/qCN9M1JeUSPpCiZhah0hIiLiPyVdzUhpefjkqkI5l4iIiO+UdDUjFTWnF0M1XT7EIiIiItUp6WpGaiZXDk0vioiINBZKupqRAwrplWuJiIg0Gkq6mpGatVuhpKvmtKOIiIgknpKu5qRm0lU5vehDLCIiIlKNkq5mxNVSMq+kS0RExH9KupqRmsnV/tWLyrpERET8pqSrGTlw9WLwd+VcIiIivlPS1YxULZh3zlVmXUq6RERE/KekqxmpmlxVuCqF9JpeFBER8Z2SrmakanJVXuH213Qp5xIREfGdkq7mpNpI1/4UTDmXiIiI/5R0NSNVkyvn9neoV3NUERER/ynpakaqJlflVUe6lHOJiIj4TklXM+JqTi9W9ukSERERvynpakaqTS9WVHmsoS4RERHfKelqRlyN6cVQsqWUS0RExH9KupqRA6YXqxyLiIiIv5R0NSO11nQp5xIREfGdkq5mpGpzVFe1I72SLhEREd9FlXSZWRszS451MNIwNUe6Kp/3IRYRERGpLqKky8ySzOwiM3vLzLYAS4FNZrbYzP5qZn3jG6ZEomZz1P3HSrtERET8FulI12TgYOAW4CDnXHfnXEdgHPAFcJeZXRKnGCVCtY1uKecSERHxX0qE553gnCut+aRzbjvwCvCKmaXGNDLxrNr0YtUNrzXBKCIi4ruIRrrCJVzRnCPxFj650kiXiIiI/zwX0pvZTfEIRBqueh3X/mP16RIREfFfvdOLZvZS1YfAUOCueAUk0atWx4X2XhQREWlMIqnp2uWc+2HogZk9HMd4pAGqFdJXz8BERETEZ5FML95e4/Gv4xGINFy16UXCH4uIiIg/6k26nHOrAMwsN/h4e7yDkujU1ptLfbpERET856WQ/sm4RSExUW0bIKiy4bUv4YiIiEgVXpIui1sUEhO1rV5Uny4RERH/eUm69J27kas+i1hLUb2IiIj4QiNdzYirkWiFarmUdImIiPjPS9J1S9yikJiodfWisi4RERHfRZx0OecWmtlHZnZq1efN7NHYhyXRqG2Ta6VcIiIi/vO6DVBv4CYz+12V50bGMB5pgKrNUW9/ewlLvy0ENL0oIiLSGHhNugqA44FOZvammWXHPiSJVtXkasryrfuf11iXiIiI77wmXeacK3PO/RR4BZgKdIx9WBKd8MmVRrpERET8F8nei1U9Ejpwzj1lZguAa7ze1MxWA4VAOVDmnNMUZQzUllypOaqIiIj/PCVdzrl/1ng8G7giynsf55zbFuV7JYzacytlXSIiIn6LKOkys/up4zu3c+7/YhaRRK2iliEtTS+KiIj4L9KRrllVjn8P/K62EyPkgPfNzAH/dM5VazthZlcBVwH06NGjgbdqOWrLrZRziYiI+C+ipMs593To2Myur/o4SuOccxvMrCPwgZktdc5NqXK/R4FHAUaOHKmcIUK113Tpf6GIiIjfvK5ehBgMnDjnNgR/3wK8Boxu6DWl9tYQyrlERET8F03S1SBm1sbM2oaOgZOAhYmOo1mqJbnSNkAiIiL+i7SQvpD939Jbm9mu0EuAc85lebhnJ+A1Mwvd/z/OuXc9vF9qUVtrCKVcIiIi/ou0pqttrG7onPsGODxW15P9au08r6xLRETEdwmfXpT4USG9iIhI46WkqxlRywgREZHGS0lXM1LbiJYGukRERPznOekys5viEYjEQK0lXcq6RERE/FZvIb2ZvVT1ITAUuCteAUn0akuutOG1iIiI/yJZvbjLOffD0AMzeziO8UgD1DqNqKRLRETEd5FML95e4/Gv4xGINFxtSZemF0VERPxXb9LlnFsFYGa5wcfb4x2UxJYK6UVERPznpZD+ybhFITGh2UUREZHGy0vSZXGLQuJKzVFFRET85yXp0nfuJko5l4iIiP800tWMOGVXIiIijZaXpOuWuEUhcaeETERExF8RJ13OuYXxDETiSw1SRURE/OVpGyAz+56ZtQ0e/8bMXjWz4fEJTbyqK6/SSJeIiIi/vO69+FvnXKGZjQNOAJ4A1KG+CVDKJSIi4i+vSVd58PfTgEedc28BabENSaJV12CWBrpERET85TXp2mBm/wTOB942s/QoriE+UK8uERERf3lNmM4D3gNOds4VAO2BX8Q6KBEREZHmJsXLyc65IuDVKo83AZtiHZREq/bRLA10iYiI+EtTgy2EUym9iIiIr5R0NSMqpBcREWm8lHS1ECqkFxER8Ve9SZeZnWhmj5nZ0ODjq+IelcScUi4RERF/RVJIfwVwNfAbM2sPDI1rRBK1ujvSJywMERERCSOS6cVC51yBc+5G4CRgVJxjknhQ0iUiIuKrSJKutwDMLNc5dzPwTHxDknhQTZeIiIi/6k26nHOTgodPBh/fH9eIJGp1rl5MXBgiIiIShpfVixa3KCTunEa6REREfOUl6dJ37Uaurgao+sMTERHxl0a6WggNdImIiPjLS9J1S9yikLjT9KKIiIi/Ik66nHML4xmINJwK6UVERBqvSJqjVjKzkcCvgZ7B9xrgnHND4hCbxJAGukRERPzlKekC/g38AlgAVMQ+HImXuorsRUREJP68Jl1bnXNvxCUSabC60qoK5VwiIiK+8pp0/c7MHgc+AvaFnnTOvRrTqCTmVEgvIiLiL69J1+XAACCV/dOLDlDS1QjUlVgp5xIREfGX16RrlHOuf1wiEREREWnGvPTpAphmZgPjEonElUa6RERE/OV1pOsIYJ6ZrSJQ0xV1ywgzSwZmARucc6d7fb94U6GsS0RExFdek66JMbz3dcASICuG15RaKOUSERHxl6ekyzm3JhY3NbNuwGnA7cD/i8U1pW5avSgiIuIvTzVdZva0meVUedzOzJ6M4r73AL+klgarZnaVmc0ys1lbt26N4vItk7YBEhERaby8FtIPcc4VhB4453YAw7xcwMxOB7Y452bXdo5z7lHn3Ejn3Mi8vDyPIUo4GukSERHxl9ekK8nM2oUemFl7vNeFjQXONLPVwAvABDN7zuM1JIy6tvpRziUiIuIvrwnT34HpZvbf4OPvEajLiphz7hbgFgAzOxa40Tl3icc4xCPlXCIiIv7yWkj/jJnNAiYEn/quc25x7MOSWNNIl4iIiL+8jnQRTLJikmg55z4BPonFtaS+QnplXSIiIn7yWtMlTVRF2HWiIiIikihKuloIjXSJiIj4y9P0opmFa2S6E5jtnJsXk4gkanVOLyrnEhER8ZXXka6RwE+ArsFfPyawNdBjZvbLGMcmIiIi0mx4LaTvBgx3zu0GMLPfAW8BRwOzgb/ENjzxoq7BLG14LSIi4i+vI10dgX1VHpcCnZxze2s8L42EWeB35VwiIiL+8jrS9W/gSzObFHx8BvAfM2tDjNpISGwZgREw5VwiIiL+ijjpMjMDngLeIbCVD8BPnHOzgscXxzY08Src/opmBs5p70URERGfRZx0Oeecmb3tnBsMzKr3DdIoJBmUAxXKuURERHzltaZrjpmNikskEhdGsKhLE4wiIiK+8lrTNQa42MzWAHsIlgw554bEPDLxLGxapUJ6ERGRRsFr0nVyXKKQuFPOJSIi4i9PSZdzbk28ApEYCJNZJWmkS0REpFHwVNNlAZeY2a3Bxz3MbHR8QpNYSAo26lJzVBEREX95LaR/CDgSuDD4uBB4MKYRSVwo5xIREfGX50J659xwM5sL4JzbYWZpcYhLouDCzC+GRrrCvSYiIiKJ43Wkq9TMkglWD5lZHlAR86jEk9fnbqDXzW+xZdeBOzFVNoxQziUiIuIrr0nXfcBrQEczux2YCvw55lGJJy/MXAvAyq27D3xRhfQiIiKNgtfVi/82s9nA8QS+nZ/tnFsSl8gkJva3RlXWJSIi4ievNV0455YCS+MQi0Qp1HU+3GiWWe2viYiISOJ4nV6URixcXmVW+2siIiKSOEq6moHKxCrMcNb+QnqlXSIiIn5S0tWMhB/pslpfExERkcSpt6bLzC4F1gM/BcqAKc65h+MdmETO9lfLH/gatb8mIiIiiRNJIf0o4HDn3PcAzOze+IYksbS/pktZl4iIiJ8iSbp2Ad3M7EfADqBNfEOSaIVPrLR6UUREpDGIpKbrt8DrQCcgDfhZPAMS7+pKqEzNUUVERBqFepMuF/A6MNo59x/n3N74hyVe1Jl0hc5JSCQiIiJSGy+rF63+U8QPFcGsK3xzVIKvKe0SERHxk5ekS9+1G6lQPlURtk+XWkaIiIg0BhrpagYiG+lKYEAiIiJyAC9J1y1xi0IapDyUdIV5bX+mrKxLRETETxEnXc65hfEMRKJXEcynwm4DpA2vRUREGoVI+nQBYGYDgLOArsGnNgBvOOeWxCMw8SCYUVXUkVgp5xIREfFXRCNdZnYT8AKB2aoZwV8GPG9mN8cvPIlERV2F9KrpEhERaRQiHem6EhjknCut+qSZ3Q0sAu6MdWASuboK6UO0DZCIiIi/Iq3pqgC6hHm+c/A18VFlTVeY15JU0yUiItIoRDrSdT3wkZmtANYFn+sB9EXbAvnOVY501TG9mMiARERE5AARJV3OuXfN7BBgNNUL6Wc658rjFZxExrnqv1dVuQ2QhrpERER8FfHqRedcBfBFzefN7HLn3L9iGpV4EhrNCl9Ir562IiIijYGX5qi1+b2Xk80sw8xmmNl8M1tkZp7eL7ULvw1QgAa6RERE/BXRSJeZfVXbS0Anj/fcB0xwzu02s1Rgqpm945w7YBRNvAmbWFXWdCnrEhER8VOk04udgJOBHTWeN2Calxu6QHHR7uDD1OAvZQQxUGfLCP0fFhER8VWkSdf/gEzn3LyaL5jZJ15vambJwGwCqx8fdM596fUasl/lVj9hcle1jBAREWkcIqrpcs5d6ZybWstrF3m9qXOu3Dk3FOgGjDazw6q+bmZXmdksM5u1detWr5dvcUJ1W+G2Aaqs6UpUMCIiIhJWpNsA1bsELpJzanLOFQCTgYk1nn/UOTfSOTcyLy/P62VbnLpWL4aoZYSIiIi/Il29ONnMfmZmPao+aWZpZjbBzJ4GfhDJhcwsz8xygsetgBOBpR5illqEy6sqpxcTHIuIiIhUF2lN10TgCgIbXPcGCoAMIBl4H7jHOTc3wmt1Bp4O1nUlAS855/7nKWqppq6u86b5RRERkUYh0o70xcBDwEPBNg+5wN7g9KAnzrmvgGFe3yf1q2sKUS0jRERE/BVxR/oQ51wpsCkOsUgDqWWESNNXVl7B8s27WZO/h007i/l2VzH5u0vYubeUfWXlZKQm06FNGu2DvzpkptG+TTod2qTRo0NrsjJS/f4SRKQWnpMuaXwsuEaxrm2AlHOJND7OOTbuLGbe2gLmrdvBvHUFLNiwk+LSispz0lKSyG2TRlarVDJSk9lauI956wrYsaeEsjBLljtlpdOvY1v6H9SW4T3aMbJXOzplZSTyyxKRWijpagb2r14M81rwd410iTQO2/eUMGX5ViYv28K0r/PZWrgPCCRXg7pkceHoHgztnkPfjpkclJVB+zZpYfdQdc6xa28Z+Xv2sX1PCdt2l7A6fw8rNu9m5ZZCnvtiDU9MXQVAt3atGNc3l2P7d2Rcv1wy0/XRL+KHBv3LM7O/AW2Ah4O1WuKD/YlVuJGu4Gsa6xLxRVl5Ba/N3cDtby8BYNfeUiocdGiTxrh+uQzv0Y6h3XM4tHMWaSmRb4drZmS3TiW7dSp9wnTWKSmrYMmmXcxas4OZq7bz1oJNvDBzHanJxuje7TljSBdOOawz2a01HSmSKA39cec2oAz4vZl1cs590PCQJFrhRrMqky7lXCIJVRpMth74eCVrtxfRrnUqg7vlMKJHO47pn8eQrtkkJXlubxixtJQkDu+ew+Hdc7hyXG9KyyuYs2YHHy/bwvuLNnPzqwv47aSFHNu/I98Z1pUTB3YiNTnypE9EvGto0nU6ga7yrYDHgZ4Njki8q2MboFC9l3IukcQor3C8Omc99328gnXb9zK4azaPf38kxx/aMew0YaKkJicxpk8HxvTpwM0TB7Bwwy4mzdvAG/M38sHizeS1TeeCUd25YHQPuua08i1OkeasoUnXPmAKsAG4oeHhSDTq3AZIRV0iCTNj1XZ+/+YiFm3cxcDOWTzxg0FMGOBvshWOmTG4WzaDu2Vzy6mH8unyLTz3xVoemLySByevZMKATlwxrhdH9unQ6GIXacoiSrrM7Bnn3PdrPu+cey32IYlXdW0DpN6oIvG3fkcRf357CW8v+JbO2Rncc/5Qzji8C8lxnD6MleQkY8KATkwY0In1O4p4fsZaXpixjose28ygLln8cHxvThvcxVO9mYiEF+lI1+DQgZm975w7KU7xSEPU0ZJeA10iseec44WZ67j9rSWUVziuP6EfPz76YFqlJfsdWlS6tWvNL04ewM8m9OP1uRt4fOoqbnhxPne9s4zLx/biojE9aKs+YCJRizTpqvotWztQNzL7pxfrGOlS1iUSUxsL9nLzqwuYsnwrR/bpwF/OHUL39q39DismMlKTuWB0D84b2Z1Pl2/l0SnfcMc7S3lg8kouPaInV4zrTW5mut9hijQ5kSZdB5nZZcB89n8fl0ZmdX7RAc/VtS+jiHjnnOO/s9fzxzcXU1bh+ONZg7h4TM+4rkT0S1KScdyAjhw3oCNfrS/gn59+w8Offs2Tn6/i4jE9+fExfejYVo1XRSIVadJ1GzACuBzoZmYLgEXBX4udc6/EJzxpKNXRi8TO5l3F3PzKV0xetpXRvdvz13OH0LNDG7/DSogh3XJ48OLhfL11Nw9N/pqnpq3m31+u4dIjenL1sX1p3ybN7xBFGr1IN7x+tOpjM+tGoM5rCHA2oKTLRzVXF5ntT7K0DZBIbHyweDM3/nc++8rK+d0ZA/nBkb2a5ehWfQ7Oy+Tv5x3Ozyb05b6PV/DE1FU8P2MdPxzfmyvH9VbNl0gdomoZ4ZxbD6wH3oltOBKNuj72VdMl0jCl5RXc9c5SHp+6isO6ZnHfBcPok5fpd1i+65XbhrvPG8rVxxzM399fzj0fruCeD1dw8qBO3HbmIDpnq9eXSE3agKsZqNlGx9g/sqUWOyLRW7+jiJ89P5e5awv4/pE9+fVph5Ke0jRXJsZLv05teeTSEcxfV8CFj33Be4s289mKbVx6ZE9unjhAfb5EqlDS1cxVdqTXQJeIJ58u38p1L8ylrNzx4EXDOW1IZ79DatQO757D4j9MZNbq7fzm9YX889Nv+GTpVt674Wi/QxNpNDx1uzOzgWGeOzZWwUh0rK4JRm14LeKJc44HJ6/ksn/N4KCsDN782TglXB6M7NWel68+CoBlmwu58b/zfY5IpPHw2mL4JTO7yQJamdn9wB3xCEyiF244XyNdIvUrLi3np/+ew1/fW8aZh3fhtZ+OpXduy1idGEuZ6Sks/eNEDuuaxcuz1/Pnt5dQVl7hd1givvOadI0BugPTgJnARmBsrIMSj8LUdIUkqU+XSES2Fu7j8n/N5N1F3/LrUw/lnvOHNtnO8o1BRmoyr1x9FGcN7cKjU77h3Eem8/XW3X6HJeIrr0lXKbAXaAVkAKucc/rxxWd1r15UTZdIXYpLy3n8s2847m+fMHvNDv5x3lB+dHQfFYDHQHpKMvdeMIz7LxzG6vw9nHrvZzz1+SqtppYWy2sh/UxgEjAKyAUeMbNznHPfi3lkElOq6RI50Esz13HXu0vJ31PCMYfk8bszBqodRByccXgXxvRuz02vfMVtby5m8rKt/PXcIXTMUjd7aVm8jnRd6Zy71TlX6pzb5Jw7C3gjHoFJZLbsKj7guao/oCcF/4T1g6XIfs45/vnp1/zyla/YXlTCI5cM56nLRynhiqOOWRk8edko/nj2YXy5Kp+J937G+4u+9TsskYTyNNLlnJsV5rlnYxeOeDFvXQFnP/j5Ac9blU5dda5sFGmBnHNc9exsPli8mdOGdOZv5x6u2q0EMTMuPaInR/Zpz/UvzuOqZ2dzwaju/Pb0gbRJVwcjaf48/S03s1vDPe+c+0NswhEvlmzaFfG5qqEQCfw7+L8X5vHB4s0cfUge918wrEVu5eO3vh3b8urVY/nHh8t55NOv+eKbfO65YBhDu+f4HZpIXHmdXtxT5Vc5cArQK8YxSYTKK+pPpEJTjcq5ROCpaat5c/5GTjnsIJ6+fJQSLh+lpSRx08QBPP+jIygpq+Cch6fxwMcrIvpcE2mqvE4v/r3qYzP7G/BeTCOSiNX64RTm+4g+xqSlW7hhJ3e8vZThPXJ44KLhWp3YSBzRpwPvXH80v35tAX97fzlTVmzjnvOH0iVHezdK8+N1pKum1kC3WAQi3kU20qWWESLrthdx0WNf0L5NGo//YBTJGuFqVLJbpXL/hcP42/cOZ9GGnUy8ZwpvfbXJ77BEYs5rTdcC9g+aJAN5wB9jHZREpiKCTCr0rUUtI6SlWrG5kBP/MQWAp64YTfs2aT5HJOGYGeeO6MbInu247oW5XPOfOXy6vBu/O2OQiuyl2fD6N/l09iddZcBm51x5bEOShqr6M7xquqQl+3DxZn74TGDR9e/PHMTwHu18jkjq0yu3DS9ffRT3fLichz75mpmrd3DvBUMZ0i3H79BEGiyipMvMCqmlLMjMcM5lxTQqiUhtUyRVS1X2j3SJtBxbCosZfftHlY9vOWUAPziql38BiSepyUn84uQBjO+Xxw0vzuO7D03j5yf158dH99HiB2nSIqrpcs61dc5l1fYr3kFKeEm1FAKH7c2loS5pIUrKKqolXA9fPJwfH3OwjxFJtI7o04F3rhvPSYM6cde7S7nwsS/YWLDX77BEohZR0mVmzwZ/vy6+4YgXkf7AZ6aRLmk5DvnNO5XHS/84kVMGd/YxGmmonNZpPHjRcP5yzhAWbNjJqfepk700XZGuXhxhZl2AK8ysnZm1r/orngFK7SJd8m5ooEtahlteXVB5PO/WE8lIVaf55sDMOG9Ud964dhzd2rXiqmdnc/tbiyktr/A7NBFPIi2kfwT4COgDzKZ6rbYLPi8JVuv0otV8bFq9KM3eoo07eX7GWgCm/OI4clprlWJz07djJq9cfRS3vbGYxz5bxdJvC/nH+UPJzUz3OzSRiERa03Wfc+5Q4EnnXB/nXO8qv5Rw+SSSgS6HRrqkeSspq+Bnz8/ltPumAnDjSYfQo0Nrn6OSeElPSeaO7w7mju8O5stV2znrgc/5an2B32GJRMRTc1Tn3NXxCkS8i3QNj2q6pDk75Dfv8Ob8jQBkZaRw7YR+PkckiXDh6B688pOjAPjeI9OZNG+DzxGJ1K+hHemlEbIax4ZppEuapcnLtlQeXzCqO3N+e6KP0UiiDe6WzRvXjuXwbjlc98I8/vreUiq0d6M0Ykq6mrDapherFtg7AFNHeml+Kiocl/9rJgB/Ovsw7jxnCCnJ+khraTpkpvPcD8dwwajuPDj5a656dja795X5HZZIWJ4+oczsrkiek8QI248r7HloflGanT6/ervy+JIjevoYifgtLSWJO747mNvOGMjkZVuYeM8UVmwu9DsskQN4/bEw3Nj9KbEIpKX5cPFmXp8b/xoEQzVd0vy8u3D/ZsiTbzzWv0Ck0TAzLhvbm6cuH0VhcRlnPfh5Za2fSGMRaXPUq4ObXfc3s6+q/FoFLKjv/VLdizPX8sNnZnH9i/NwDSm2qm16scpxYPWiNew+Io1IYXEpP3luDgAPXjSc3rltfI5IGpPx/fJ4/4ajGdg5i/97YS7//PRrff5JoxHpSNd/gDOAN4K/h36NcM5d7OWGZtbdzCab2WIzW9QSu9w//tmqyuP1O6Lf0qLWycUD+nSpZYQ0D4XFpQy+7X0AxvXN5bQh6jYvB+qUlcFzPxzDqYd15o53lnLDi/MoLi33OyyRiPt07XTOrQYuBsYDP3DOrQEyzWy0x3uWAT93zg0EjgCuMbOBHq/RZG0s2MuKLbs5e2gXAOatK4j7PQ1NL0rzEEq4AB79/ggfI5HGLiM1mQcuGsaNJx3CpPkbOfeRaazfUeR3WNLCea3pehA4Ergw+Lgw+FzEnHObnHNzgseFwBKgq8c4mqy5awsA+MFRvUhPSWJ+HJKumiNgZmoZIU3f0X+ZXHk87eYJtE6LdEMNaanMjGsn9OOxS0eyZlsRZ9w/lc9WbPU7LGnBvCZdY5xz1wDFAM65HUDUe22YWS9gGPBljeevMrNZZjZr69bm9Q9kxZZCzGDAQVn06tCGNduj/8nL096LGuuSJuzv7y9jbfDfyuI/nEyXnFY+RyRNyQkDOzHp2rHktU3nB0/O4KFPVqrOS3zhNekqNbNkQu2fzPKAqHYcNbNM4BXgeufcrqqvOecedc6NdM6NzMvLi+byjdaKzbvp0b41rdKS6dauFRviUdMV5kR9vkhTVV7huP/jlQD89ydHaoRLotInL5PXfjqWUwd35i/vLuPq5+aon5cknNek6z7gNaCjmd0OTAX+7PWmZpZKIOH6t3PuVa/vb8pWbCmkX8dMALq2a9WgGoNI86iIkzORRujej1YAcNLATozq1d7naKQpa5Oewv0XDuM3px3KB0s2c85D01jXgNkGEa+87r34b+CXwB3AJuBs59x/vVzDAnNiTwBLnHN3e3lvU+ecY01+UeUS927tWrGruIyde0vjet9ATZeGuqTp+dt7y7jvoxW0TU/hkUtUOC8NZ2b8cHwfnr58NJt27uWsBz/ny2/y/Q5LWgjPe2Y455Y65x50zj3gnFsSxT3HApcCE8xsXvDXqVFcp8nZvqeEfWUVlfUo3dq1Boh6irG2RKpmrZeao0pT9OvXFvDA5MC04jvXjycpSWO2Ejvj+uXy+jVjyWmdysWPf8nT01brh1OJO0/FEWaWDpwD9Kr6XufcHyK9hnNuKi10xmvTzmKAyqSrU1YGAJsLixlIlufreZle1GeJNCXTvt7Gv79cC8DUm46r/AFFJJb65GXy+jVj+X8vzuN3byxi/voC/vydwWSkJvsdmjRTXke6JgFnEei1tafKL4nAhoLAiFbXYNLVsW06ANsK98X0PjV/WjMzrV6UJmPTzr1c9FhgQfMtpwxQwiVxlZWRyqOXjuSGEw7h1Tkb1M9L4srrMqBuzrmJcYmkBdgYTLo6ZwdGuHIzA0nX1t1RJl0R5lEa6ZKm5I//WwzAXecM5vxRPXyORlqCpCTjuhP6cVjXLK5/YR5n3D+VBy8azlF9c/0OTZoZryNd08xscFwiaQE2FuwlPSWJ9m0Crc1apSWTmZ7C1ihHumobvVJNlzRFzjkm/P0T3l7wLWcP7aKESxLu+EMD/bw6ZKZz6ZMzePyzb1TnJTHlNekaB8w2s2XBDa8XmNlX8QisOdq4s5guOa2qJUW5mWls210S5zurI700fg998jXfbA1UK9x17hCfo5GWKtDP6yhOPLQTf3prCdc+P5c96uclMeJ1evGUuETRQmzdtY+8YB1XSF7bdLYWFkd1vUgTqUCOp6xLGq8Vmwv563vLAJh/60mkp6iQWfzTNiOVhy8ZziOffsNf31vK8m8LeeTSERycl+l3aNLEee3TtQbYBXQCelb5JRHYtmcfeZnVk67czPQGTC+GV3V20TnVdEnj988p3wAw6ZqxZLdO9TkakUCZxtXHHsyzV44hf08JZz3wOe8u/NbvsKSJ85R0mdkPgSnAe8Dvg7/fFvuwmqdthfvIzay+VWX7NmnsKIptc9SaCZZpGyBpxFZu2c3Ls9dz4ejuHN49x+9wRKoZ2zeX//1sHAd3zOQnz83mrneXUl6hD1SJjtearuuAUcAa59xxBDarLoh1UM3RvrJydhWXVa5YDGnXOo2CohIqovhHHEkiZQaGWkZI47SvrJzrX5wLwGVH9fY5GpHwuuS04qUfH8GFo3vw8Cdf84MnZ5Af7apzadG8Jl3FzrliCDRKdc4tBfrHPqzmJz9YLN+hRtKV0zqVCgeFURRqRpJIOaeRLmm8Hvh4JQs37OKRS4bT/6C2focjUqv0lGTu+O5g/nLOEGas3s4Z909l7todfoclTYzXpGu9meUArwMfmtkkYE2sg2qOtgV/Kqo5vZjdKlC/sjOGU4w1OkYEarpidnWRhttZVMrFj3/B/R+v5IzDuzDxsM5+hyQSkfNGdefVq48iOdk475/TeWa6tg+SyHktpP+Oc67AOXcb8BvgcQId6qUeoZGu3LYHTi8C7Cjy3jYi4ulFU8sIaVx+/t95fL4yn6MO7sBtZwz0OxwRTw7rms3/rh3P+H553DppEdc+P5fC4tjW5krzFFHLCDMrZP9giYU59r5xYAsT6jpfc/ViTnClVsFe7/9gveRRqumSxmL61/l8uGQLv5zYn58e29fvcESikt06lce/P5KHP/2auz9YzqINO3no4hEM7KJvh1K7iEa6nHNtnXNZwV8HHMc7yOYgNL3Yocb0Yk5wpKsgipGuSJnmF6WRqKhw3PHOEjpnZ3DFWBXOS9OWlGRcc1xfXrjqCPaWlvOdhz7nxZlrNd0otYoo6TKzDDO73sweMLOrzMxrU9UWb1thCa3TkmmdVv1/XeVIl8earpdmruO3ry+M6FxtAySNxf8WbOKr9Tv5+Un9yUhVA1RpHkb1as9b/zeeUb3ac9MrC/j5f+dTVKIu9nKgSJOnp4FS4DPgVGAQgfYREqFtu/cd0C4CICdYSO+1pus3kyJLuCDYMkI/eYlPnp2+msnLtnLyoE788X9LGHBQW74zrKvfYYnEVG5mOk9fMZr7P17BvR+tYOGGnTx08XD6dtSqXNkv0qRroHNuMICZPQHMiF9IzVP+ngMbowKkJCfRNj3F80iX1X8KUKVlhKeri8TGsD+8X9n89+OlW0hOMu69YBjJSZH+DRZpOpKTjOtPOISRPdtz3QtzOfOBz/nzdwZztn7IkKBIk67KjMA5V2Y1exJIvbYVltCjQ+uwr+W0SfVc05Xk4c9A2wCJHz5bsbUy4frH+Yeza28ZJwzsRNecVj5HJhJf4/rl8vZ14/nZf+Zy/Yvz+HLVdn53xkBNqUvESdfhZrYreGxAq+DjwPdzFdPXa3tRCcN65IR9LaeV962A6hooqPpSZcsIT1cXabhLnwgMiP/r8lEc17+jz9GIJFanrAz+86Mx/P2D5Tz8ydfMWbOD+y4cpibALVykqxeTa6xYTNHqxcg55ygoKqlcqVhTTutUzy0j6hptrJpg7d/wWmmXJM5rc9dXHivhkpYqJTmJmyYO4OkrRpO/p4QzHpjK09PUTLUl89qRXqKwp6Sc0nJHu+BKxZqyWqVS6Dnp8nKyarokcUrKKrjhxfkAnD+yu8/RiPjvmEPyePf68Yw9uAO/e2MRP3x6lvZubKGUdCXAjj2Beq12tYx0ZWWksqvY2/Liumq6DpheBGVdkjDvLfq28viPZx/mYyQijUduZjpPXjaK284YyGcrtzHx3s+Ysnyr32FJginpSoDQysScWke6UtjlcQuJSBd/BVYvmjrSSzVl5RUs31wYl2v/7Pm5ALxx7VjSUvQRIxJiZlw2tjeTrhlLTqtUvv/kDP70v8XsKyv3OzRJEH0iJkDB3sBIV201XVkZqZSUVVBcGvk/PC8rSLV6UWr66/vLOOkfU1i1bU9Mr3vzK18BkJ6SxJBuOTG9tkhzcWjnLN782TguPaInj09dxXcenMbKLbv9DksSQElXAoRWJtZa05URWETqZbTLS5sjMyVdUt3MVduB/dtTxcKu4lJemLkOgC9/dXzMrivSHGWkJvPHsw/jse+PZNPOvZx+/2c8P0NbCDV3SroSINSDq9aRrmBX+kIPdV3eRro0vSjV7dkX2+kM5xxDbnsfgLSUpFr/rotIdScO7MS71x/NiJ7tuOXVBVz93Jy47sUr/lLSlQA79tRT05UReH6XhxWMnhYvaqSrSVq4YScbCvbG5drLgvVcsfp7MXnZlsrj2b85ITYXFWkhOmVl8OwVY7jllAF8tHQzE+/5jE9VZN8sKelKgB1FJbRNTyE1Ofz/7qxWoelFLyNd3mJQztW0VFQ4Tr9/KmPv/Djm1z789+9XHpeWVzT4es45rnhqFgBXHd2Hthnhf7gQkdolJRk/PuZgXr16LJkZKfzgyRnc8uoCdu/TxtnNiZKuBCgoKiGnTe3fiKIZ6fLCzDTS1YSs215En1+9Xfm4qCR2H7rlFY6dVf6eXfz4l1RUNOwvxxvzNwKQmmz84uT+DbqWSEs3uFs2//vZOH58dB9emLmWifdMYdrX2/wOS2JESVcCFOwtJadV7TUuoZquSArpF6zfyXUvzKWu75M1670Cj5R1NRU1Wzls2RW7Yvf3q/TQCsnfE339SEFRCde9MA+AKb88rtbRXBGJXEZqMreceigv/+RIUpKMix77ktveWMTeErWWaOr0CZkAO4pKa63nAmgbXL0YSSH9RY99waR5G9laGPk3YtV0NS1Lv62edMVyeiFcUrSlsDiqa1VUOIb+4QMAfjmxP52ztZG1SCyN6Nmed647msuO6sVT01Zzyr1TmL1mu99hSQMo6UqAgqKSWrvRA7QK7jz/zsIDRyFqKo8iezJtA9Sk/PW9ZdUe74lh0pWSHBj3HNRl/5apc9YWRHWtD5Zsrjz+6bF9GxSXiITXKi2Z284cxH9+NIayCse5j0zn9rcWa9SriVLSlQA79pTU2qML9k8Hzl9XUO+1PNbPB99j6v3ShBXF6MPVOcdl/5oJwF/OHcLUm44D4LevL/R8rU079/LjZ2cDMOUXx8UkPhGp3VEH5/Lu9UdzwagePPbZKk6+Zwqfr1StV1OjpCvOysor2FVcFnHfovpWk5VEsdpMI11NR7hl4mUNLHQPeXTKN5XHPdq3Jq9telTX2bGnhCPvCKyq/N0ZA+nRoXVM4hORumWmp3DHdwfznx+NIckCC2F+8d/56uvVhCjpirPQSrG6Rrqq2l1PXVdpeRTTi6imq6l46vNVBzxXXtHwtg5QPaFrm5FKekpy5eOSssjvMeyPgTquC0d35/KxvWMSm4hELjTqdfWxB/Pq3A2ccPen/O+rjZrRaAKUdMVZ5RZAbeoe6frb9w4HvHWlj5iZRrqaiMnLDhzpiibRrundhZuY9nU+AMf2z6t8/prjDgbgmemrI7rOll37i+5vO3NQg+MSkehkpCZz08QBvHHtWDpnt+La/8zlR8/MYtPO+DRUlthQ0hVnO4ObXWe3qnukKzPd+/6LtalZ9xUY6VLa1ZR8d3hXbjjhECDQW6uhfvLcnMrjUb3aVx6fPbQrAH96a0m91ygqKWPsXR+TkZrE+zccXW2kTET8MahLNq/99Ch+c9qhTF25jRPvnsIz01fH5HNDYk9JV5yFtgCqa/Ui7E+6YrFSzdU49tq9Xvx38ykD+O7wQEIU6+1AfjS+T+Vx9/b767Hq+5D+yXNzKC13XDG2N4d0ahvTmEQkeinJSfxwfB/ev/4YhvXI4dZJizj7wc8jWpwliaWkK852BAsc60u62qQHRg32xLD7eIhqupqerIxUkpIC2fJrczc06Fo1O86npez/Z5+Run+06o35td/n1TnrmRJM/v7v+H4NikdE4qNHh9Y8c8Vo7rtwGJt3FXP2Q5/zq9cWqNC+EVHSFWcFwZquurYBgqojXbHtvWIEtwFSVVeTkpGa7Km4vS5Lvt1Vebz0jxNrPa+0LPzfkclLt/D/XppPq9RkZv3mhGqJmog0LmbGmYd34aOfH8PlR/XmxZnrmPD3T3lp1roGb/klDaekK852FJWQkmS0DSZVtWkTw+nFqrOJDo10NRWhurtLj+gJQOfsDAA6RtnaIeSjJVsqj8MlTO9dfzQQvp5wbX4Rlz8V6O31xA9GkpvZsFhEJDHaZqRy6xkDefPacfTJbcMvX/6K7z48jdlrdvgdWouW8KTLzJ40sy1m5r0jYxMU2gKo5n6INbVJCyRddW35Eu3IR2gboNvfWsySTbvqf4P44pPg9F0o+clITaZLdgZHH5JX19vqNaWemrBDOmWSlpzE1t3Vt5b6Zutujv7rZADuvWAoR/XNbVAcIpJ4A7tk8dKPj+Rv3zucjQV7OefhaVz7nzms217kd2gtkh8jXU8Btc9xNDMFRSURNUatrOmqY3ox0m0fqg5qGYGO9DuKSnjss1Vc9NgXEV1DEu/yYLf4UJsRCBTIlkXREDekosIxq56fbM2MvLbp1fbz3FdWzoS/fwrA3ecdzlnBVY4i0vQkJRnnjujG5BuP5brj+/Hhks0c//dPueOdJTFZMS+RS3jS5ZybArSYHTt3FNW9BVBISnIS6SlJdRbSF5V6n3p0AAYVml9sMu49f2jlcWqyNahPV8He/R+ofXLb1Hpebtt0Nhbs5eR/TOHhT76m/2/eBeDXpx7Kd4d3i/r+ItJ4tElP4YYTD+GTG4/jjMO78M9Pv+HYv37Cs1+sadAPdxK5RlnTZWZXmdksM5u1dWtsl8snWkFRKdmtItsCKDM9pc6arkj34AvXpytUP1nfNKf4r2oj3dTkpHq3hqpL1dGrj288ttbz8jLT+OKb7SzbXMhd7y4F4LCuWfzo6D61vkdEmqaDsjP4+3mH8+a14+jXMZPfvr6Qifd+xuRlW9TTMc4aZdLlnHvUOTfSOTcyL69h9Sx+KygqjXgLoDb1JF3R7CofWL2o5qiN3b6y8H+2KckW9d6LBUUlnHzPFKD+baj2hakXnHTNuKjuKyJNw+Bu2bxw1RE8eukIyiscl/9rJt9/cgZLv1Xtb7w0yqSrOdlRVFLvFkAhrdOS2V1XTVep96QrsHrRKlcvapyrcaqtlq8hI10vzFxXefz6NWPrPPfQzlmVx+9cN57Vd55GcpL+tog0d2bGSYMO4r3rj+bW0wfy1fqdnHrvZ9z43/lsKNCWQrFWdx8DaZC9JeXsK6sgJ8KRrphNL9qBj0M1XZpdbJxmrQ5f5piaFH3Sdec7SyuPe3aovZ4L4PoT+jGkWzanDe6sKWiRFigtJYkrxvXmu8O78sDHK3lm+hremLeRi4/owU+P7UteA1vXSIAfLSOeB6YD/c1svZldmegYEiXSbvQhbdJTKKqjkH5vhN3qa84kmqHWqI1YUUkZVz07O+xrKclGWQM3vL73gqH1ntM6LYXTh3RRwiXSwuW0TuM3pw9k8i+O5TvDuvLM9DUc/ZfJ3PXuUnW2j4GEj3Q55y5M9D39sj/pinyka/2O2nunRDrSVZNhlfvqbdutfzSNzYdVmpfWtH1PCUu/LfR8zarJ+5mHd4kqLhFpubrmtOKuc4fw42P6cM+HK3jk0695bvoarhjXmyvG9Sa7VWTf16Q61XTFUeUWQBGPdCXX3acripouCE4vavuHRiu9yl6IH9xwdLXXNuyIrqbi7veXVx5r9EpEotUnL5P7LhzGO9eN56i+Hbj3oxWMu/Nj7n5/mUa+oqCkK468Ti+2TovN6sVw32PLtXqx0VofTKxOHNiJfp3aVnvtnBGBHlleV58+PnUVAId3z2l4gCLS4g04KIt/XjqSt/5vHOP65XLfxysZd9dk/vreUrbvUfIVKSVdcRTqLO5lenFPSVmt32CjmV4MbXitvnexMW3lNh7/7Jtqz63YXMjKLbu5+/1lldO4Xvzxf4sBOCPMNGDHrEDxariWDpF49srRUb1PRCScQV2yefiSEbx7/XiO6Z/HQ598zbi7PuaOd5aQX2MrMTmQVi/G0c7gSFe2hz5dFQ6KSytolXbgxsR1dauviwHlFcq6YuGix78EIK9tOmcN7YpzjhP/MaXy9VG92zO+X3S95cIl262DG1QXlZSH3aw6nKrLvLMyVHchIrE34KAsHrxoOMs3F/LAxyt5dMo3PDNtDReM7s5VR/ehc3Yrv0NslDTSFUc7ikppnZZMekpk3ywzg/sv1rbpdVTNUS3wK5oRGKnd+4s2A3DvRyuqPR9tI1Mg7CrF1sGN0Ota1VrTjS/NB+DH6iYvInF2SKe23HfhMD644RhOOewgnp2+hvF3TeaXL89n5Rbvi4CaOyVdcRTYdzGyei4IjHQBtdZ1FZWUe+6zZVTfBkiit63K0Hlo9HLu2oIGXzfUx+2Y/geOkLUOJuKPTfmGVdv21Hut8grH9G/yAbSFj4gkTN+Omdx9/lAm33gsF4/pwaR5Gznh7in88OmZzFi1XbuiBCnpiqOCotKIG6PC/qSrtpGuopIyMtO8zwgHarr0F76h5qzZUXm8MTiFd1jXrOonefzfvGdfWeUq19zMA5sPtg5OMz89fQ0T75lywOs1PVal3izc9URE4ql7+9b8/qzDmHbzBK47vh+z1+zgvH9O5zsPTePtBZta/PciJV1x5Hmkq3IqKfw0YlFJeWViFikzC9R06aeMBntj/kYAOrZNpyjY2qOwuHqCfMNL8zxdM7Qh9VlDw/fSal0lyY6kmD7Uhf7RS0d4ikNEJJY6ZKZzw4mHMO3m4/njWYPYUVTCT/89hwl//4Rnp6+OqlymOVDSFUfeR7oCoxp1TS+GzqmbVTtSn67Y+N9XmwAY3DW7clFDzaQrNGoVqWP/9gkAR/TpEPb11jUWVNT1U+LnK7dVHp806CBPcYiIxEOrtGQuPbIXH//8WB66eDg5rdP47aRFjL3rY+7+YDlbdhX7HWJCafViHHkd6cqMZHrR40hXgFUb6SopqyAtRfm2F1X/TFqnp1SORu7c6y3JqmpX8f73ThjQMew5XXKqrwDaU1JW64rE37y+EID//Wxc1DGJiMRDcpJx6uDOnHLYQcxYtZ1Hp3zDfR+t4KHJKzl1cGe+f2RPRvRs1+ybOes7b5yUVzh27i2NuEcXRFZI3zqimq79CVZo9WLV2cXvPPR5xDFJwE2vfAVAarLRJi258s9o087ofkr7dmcxQ257v/Jxp6yMsOflZqbzh7MGVT7+eMkWzntkOqurFNUv+7aQ9TuKWLe9iEuO6MFhXbOjiklEJN7MjDF9OvDEZaP45MZj+f6RvZi8bAvnPjKdU++byn++XFvrwENzoKQrTnbtLcU5yI5i9WJdLSO81nQFdl6sbtHGXR6vIW8FpxZn/OqEyp0DdhWX8u3O6tv0HJzXJqLrfbkqP+J7f//IXvz13CEAXP/iPGas3s5Pntu/QfbJ90xh3F2TKatwjOsbXY8wEZFE65XbhlvPGMiXvzqeP39nMM45fvXaAkbf/iE3v/IV89cVNLtVj5pejJOCvd660QO0SQvVdNVeSJ/psaYLwm8LJJH7an1B5XG7NmmBPTJLyquNVIV0bBt+xKqmBet3Vh7/bELfes/fUli903Oolqxmgl7bNKWISGPVOi2Fi8b04MLR3Zm7roDnv1zLpHkbeWHmOg7tnMWFo7tz1tCuzWKTbSVdceJ130WAlOQkWqcls3tf+DqhopIyWntevRhurEu8+KpKggSE3S0A4Mg+HSJeDj1n7f72Ez8/qX+95w+rsYfihoK99Lr5LU4c2Kna86rVE5GmyswY3qMdw3u047dnDOSNeRt5fsZabp20iD+9tYQTD+3E2cO6cswheU32s05JV5yEdl/3snoRAtu27Np74PRiaXkFpeXOcyF9aPWiRC9UoH76kM7A/tYeNaUkG8Vl9S+Dds5VbnIdqaP65oZ9/oPFgc74t50xkPGHaGpRRJqHrIxULjmiJ5cc0ZMF63fyypz1vDl/I28t2ERO61ROH9KZ7wzrxvAeOU2q+F5JV5xs3xOaXox8pAsgq1VK2BVxodVytX3Dr02okF4a7tbTBwIHtnH48dF9GNo9h5dmrYtopOvcR6ZXThd+cuOxEd//85snMPbOj3nq8lGs217EbyctAuBPZx/GJUf0jPg6IiJNyeBu2Qzuls2vTzuUqSu28drcDbw8ez3PfbGWHu1bc/awrpw9tAt98jL9DrVeSrriZPuewDfVDpkek66M1GqtBEJCe+9F1qdrPwv+J9Hrkp3BkQfn0jG4wrDmYoZj+3fkyIM78MqcDZSG2T+xqpVbdjM72Nn+b987nF65kRXeA3TNacXqO0+rfHz+qB4s31yo1Yoi0iKkJidx3ICOHDegI7v3lfHewm95fd4GHvh4Bfd9tILDu+dw9tAunDq4c60rwv3WNCdFm4D8PSWkJSd5ng7MblVb0hUc6YqipitcznXLq195uk5LtXrbHjbuLKawyp9JzT+DrFaBx6nJRnnF/q7xHyzefMCf5QWPTq88PuWwhjUwTUtJUsIlIi1SZnoK54zoxrNXjmH6Lcfz61MPpbSsgt+/uZgxf/6Icx6exuOffcP6HUV+h1qNkq44yd9dQofMNM9zzVmtUsNOL+6NMumCsDkXz89Y5/k6LVGoY3zV4vkObaqPXoaalSYnGcs372bLrmI27dzLj56ZxbkPT2Pd9sA/+vzd+9i2O1Drl5JkUf1ZiohIdZ2yMvjR0X14+7rxfPj/juHGkw5hb0k5f3prCePumsxZD0zl4U++rtbf0C9KuuJk+54S2rfxNrUIkJWREraQPtSMs00tK+dqE6jp0vRiNELJEsDvz9zfoLTmlHEo6coPJlRX/3tOZduP5Zt3M/4vk3HOMbXKNj1PXT46bnGLiLRUfTtmcu2Efrx93Xg+/cWx3HzKAADuencpx/7tE77/5Axf49OP2nGSv3tfdElXq1QKi0upqHAkJe1PlkL9mLI89ilRRVf0xv9lcuVxTpUFER3apFc7L1RnV1oemFrM372P+z9eUe2cWyctIq9t4H2nD+nMuH7hVyOKiEhs9OzQhp8cczA/OeZg1u8o4r1Fm0n2+Ruikq44yd9TEtVKiuxWqVQ42F1jj71QM8za9t2rS5KyLs+qjnLVbPtRsz9MSnLgcahdxOr8IlbnV68jePaLNQD0zm3DAxcNj3m8IiJSu27tWnPluN5+h6HpxXiJfnox8A1+V426rlAhd6hou25VVtBpejEqVUe5vrjl+FrPq7pQYuGG+rdX6tmhdcMCExGRJktJVxzsLSmnqKQ8yunFwDfxmnVdu4IjXZFsg+Cq51yaXvToL+8urTw+bXBnMlIPrKM76uAOAMz89QmVzx3Rp3291+6c3TiXMYuISPwp6YqD/GCPrlyPPbqgykhXcc2RrjJSk42MlPoL6Q/oFKWsy5OHPvm68vhXpx0a9pznrhzDN38+tdqqxueuHFPtnGuOO5iVt59Cr+Do1vh+ufzujEGIiEjLpJquOAitYmtfo+A6EqFC+ZptIwqLS2mbkVqtuL42VXdlNzNWbN7tOY7GZs++MjJSk0mOc4HajFXbK4+//vOptd4v3J9DqLYr5Opj+5KSnMSTl43i46Vb+OH4PrENVkREmhQlXXGwfU8o6fI+0hWaPjywpquMthnR9ehasGFnvec1Zg9OXslf31sGwDd/PjWixDMaFRWO8/4ZaF56wwmHRJXgDemWzVfrd3LuiG6V9V598jKbxPYUIiISX0q64iA/mHRFM72YHVwptyO4YXZIYKQrsj+u+nf/azqKSsoqEy6AT5ZvYcKATnG51zPTV1ceXzuhb1TXeOPacTGKRkREmhvVdMVB/u5ATVc0I11t01NIS06qTNxCCovLaJseWbuIaoX0dQzWVESwObPfngu2WnjhqiPolJXOM9PXxO1eHy7ZUnkc72lMERFpeZR0xUH+nhLSUrzvuwiBGqzczDS2FYZJuiId6apa01XHeWWNPOlak7+HP7+9lH4dMzmiTwfOPLwL01bmU1xaHvN7rdteVNkx/i/nDIn59UVERJR0xcHmXcV0ykqPuj9Wh8z0yhWQIaFC+khEmkqVN/Kk68mpqwD41amBFYRH9OlASXkFc9cWxPQ+z05fzYWPfVH5+Hsju8X0+iIiIqCkKy427yqmU9vo+zF1yEyrXAEZsqOolHatvXejryvx+3JVvufrJUp5heODxZs55pA8jhvQEYCRvdpjFtu4p329jd9OWsT6HXsrn1MzWRERiQclXXGwZdc+OmVFn3TlZqazbff+ka69JeXsLS2nfaSF+RHWdF32r5mUBfcLbGzeXfgtG3cWc8Go7pXPZbdKZWDnLL74JjZJl3OOX7+2EIC0ZP1TEBGR+NJ3mjjYvKuYjlnee3SFhEa6QrVZoanGDh4L808dfBC/PjV8c8+QmnsENgYlZRX87o1FHJzXhpMGHVTttTG9OzB3bQH7yhpe1/XSrHWs2raHP519GMtvP4VHLhnBracPbPB1RUREwlHSVUVJWQWz12yv/8Q67N5Xxp6S8oaNdLVJp6S8gsJ9ga1/9vf9iiyRCw103XnOEDpkpjOwc1at5y79tv79AhPtV68tYNvuffzi5P4HrCIc06c9+8oqmL+uYb3Hyiscd7wT2O7n3BGBGq6Jhx3EFY1gQ1QREWmelHRVce9Hyznn4enMXrMj6mts3lUMQKcGjHTltg2MaG0rDIxw5XtsthoaIQulK0N75NR67pJNjSvpmrN2B6/OWc/grtlh+3GN7hXY33BGA+u6/v3lGgqKSrn19IFh91YUERGJNTVHrWLqysA38mkrtzGiZzumrthGl5wMT93EK5OuBhTSh9777c5i+uRlVhbVRzK9ePnYXrw4cx2wvyC8rn5cGwuKo44z1jbvKuaKp2aSmZ7CPy8dQVrKgT8TtGuTRnpKEn97fzkl5Y5Zq7czd20BB2VncP0J/WiTlsKYPu3rXOk5Z+0Obp20iEM7Z3H52F5x/IpERET2U9JFYHPpV2evpzC4yfSqbXuYv66AS574kqyMFKbePKFyI+r6bAiuguuc0yrqeLq3D2yQvG5HoN5qS2EgMcprW//o2a9OPZQXZgSTruBzdbWGeHP+Rs4d0Y2xfXOjjjcWNhTs5ernZlNcWs4b146jSx3///aVBYr/7/toBW3SktlbWs6qbXu47oV51c4b1zeXP519GL1y21Q+V1ZewS9f/gqAX07sr5WKIiKSMJpeBJ76fDW3vbmYb7buAeDrbXv4ZNlWAHYVl/Hq7PURX2vd9iKSDLo2IOk6KDuDJIN12wMJ3IYde8lpnUqbCJqtplSpgQrlE3UlXWUVjosf/5I7g/VNfiguLef/np/Lyi27+cd5QzmkU9s6z3/+R0cAMGFARz795XF8ddtJla+F9q4EmLpyG8f+7RN63fwWVz0zC4A35m9k5Zbd/P7MQRzXv2McvhoREZHwNNIFrN1efQXfmvw9LN9cSM8OrWmbkcJLs9Zz2djICqzXbi+ic3arsFNjkUpNTqJzdqvKka4NBXvrTOLSkpMoCbZ+MDNcjfaopRE0QX3k06+5eEyPylG2RNm+p4TL/zWD+et3cs/5QzllcOd633PkwR1Y+seJ1WqxXr9mLId0yqR1WgoVFY6563Zw6RMzKCoJrHJ8f/FmZq7ezn++XEvv3DZ8/8iecfuaREREwvFlpMvMJprZMjNbaWY3+xFDVbv2llYe9+zQmoKiUmat2U6/jm05d3g3Fm/axfLNhbW+/8WZa3klOBq2ZnsRPWKQuBzcMZPlm3cHrplf+zUf//5Ilt9+SrXnQrsAWXCCsV/H2mvSslulMr5fYGpx/F8mNzTsiDnneGfBJo7962SWbCrkn5eO4OxhXSN+f83i96Hdc2idFvgZIinJGNGzPYv/MJFFvz+Z7wav+4MnZzBrzQ4uHN1d04oiIpJwCU+6zCwZeBA4BRgIXGhmvjZHqjrSNax7DgCbd+2jX6dMThvSheQkY9K8DWHf+8KMtdz0ygJ+/t/5LFi/k6WbCul/UN3TY5EY2i2bZd/uYkthMau27eHQMG0fXrn6SE4YGFjhV3Wfx+evOoILR/cgIzXwx3vNcX0rX5t+ywQO6ZRJ15xW3Hr6QObdeiIPXDS88vVLHv+SDQV7mbl6O9O/zmfJpl2s2FzImvw95O/e1+Ctg/aVlfPuwm+55Ikvufrfc2jfJo3XrxnLyTX6ccVKm/QU7j5/KD8a35uiknIy01M4d0T3+t8oIiISY35ML44GVjrnvgEwsxeAs4DFPsSCc6560tWjHa/P2wjAgIPaktc2nbF9c3nui7Ukm9G+TRonH3YQ6SnJfLx0C79+fSHDe+SweNMurn1+DntLyxnSLbvBcQ3r0Y4KB395dxlAtWt++P+OoUObNNpVWc340c+PYWuwxcTwHu0Y3qNd5WvJScY7141nz74yOme34v0bjql2r+xWqRzXP4/Jy7YydeU2xt75cb3x5Wam0SY9hTZpKbRKS6ZtRgrZrVJpnZZMZnoKbTNSadc6lby26aSnJrNow05mrdnB7NU7KNxXRse26fzi5P78aHyfBk3FRupnx/fjsK7ZjOjZLuLWGyIiIrHkR9LVFVhX5fF6YIwPcQCwt7Scow7uwIdLtgAwvl8urVIDK+KGdQ8kLucM78qU5Vu57+OVANz25v78cGTPdjx9xWj+/v5ynvx8FekpSRx9SF6D4xrbN5euOa14efZ6OmWlc+TBHSpf6xtmurBTVkadDVnDjZRV9eRlo1i2uZCKisB+hO8v2syJAzvRMSudJDNKyioo2FvKuws30Tu3DSnJSewuLmPPvjKKy8rJ313CN1v3sHNvKfvKyikuPXB7oUM6ZXL64V04aWAnxvfLJSWBW+9kZaRy1tDIpy9FRERizUKNNBN2Q7NzgYnOuR8GH18KjHHOXVvlnKuAqwB69OgxYs2aNXGPa0PBXhZt2MlJgw5i4YadbC3cV7nRckWF460Fm+jarhVpyUnMXbuDsgpH15xWHDegI6nJSewtKeeN+Rs4tHMWQ7rlxCymj5Zs5uh+edXaHjQF+8rK2bm3lK2F+9izr5xDOmWS01ojTCIi0ryZ2Wzn3Miwr/mQdB0J3OacOzn4+BYA59wd4c4fOXKkmzVrVgIjFBEREYlOXUmXH6sXZwL9zKy3maUBFwBv+BCHiIiISMIkvKbLOVdmZtcC7wHJwJPOuUWJjkNEREQkkXxpjuqcext42497i4iIiPhB2wCJiIiIJICSLhEREZEEUNIlIiIikgBKukREREQSQEmXiIiISAIo6RIRERFJACVdIiIiIgmgpEtEREQkAZR0iYiIiCSAki4RERGRBFDSJSIiIpIASrpEREREEsCcc37HUCcz2wqsScCtsoGdCbiPF37EFO97xvr6sbheQ64R7XtzgW1R3lMi1xj/XTdEY/16mttnVWP8nGrodfRZlTg9nXN54V5o9ElXopjZo865q/yOoyo/Yor3PWN9/VhcryHXiPa9ZjbLOTcymntK5Brjv+uGaKxfT3P7rGqMn1MNvY4+qxoHTS/u96bfAYThR0zxvmesrx+L6zXkGo3x743s19z+fBrr19PcPqsa4+dUQ6/TWP/utCga6RLxgX56FJGmQJ9VsaWRLhF/POp3ACIiEdBnVQxppEtEREQkATTSJSIiIpIASrpEREREEkBJl4iIiEgCKOkSERERSQAlXSI+M7M2Zva0mT1mZhf7HY+ISDhm1sfMnjCzl/2OpalS0iUSB2b2pJltMbOFNZ6faGbLzGylmd0cfPq7wMvOuR8BZyY8WBFpsbx8VjnnvnHOXelPpM2Dki6R+HgKmFj1CTNLBh4ETgEGAhea2UCgG7AueFp5AmMUEXmKyD+rpIGUdInEgXNuCrC9xtOjgZXBnxZLgBeAs4D1BBIv0L9JEUkgj59V0kD6gBdJnK7sH9GCQLLVFXgVOMfMHkb7o4mI/8J+VplZBzN7BBhmZrf4E1rTluJ3ACItnXNuD3C533GIiNTFOZcP/MTvOJoyjXSJJM4GoHuVx92Cz4mINCb6rIoTJV0iiTMT6Gdmvc0sDbgAeMPnmEREatJnVZwo6RKJAzN7HpgO9Dez9WZ2pXOuDLgWeA9YArzknFvkZ5wi0rLpsyqxzDnndwwiIiIizZ5GukREREQSQEmXiIiISAIo6RIRERFJACVdIiIiIgmgpEtEREQkAZR0iYiIiCSAki4RqZWZ/cPMrq/y+D0ze7zK47+b2f+L4f2eMrNzY3W9Ktf9VZXjXma2MML3XW9m3/d4r2le42soM8szs3cTfV8R8UZJl4jU5XPgKAAzSwJygUFVXj8KSHiSEYVf1X9KdWaWAlwB/MfL+5xzR3m8R4M557YCm8xsbCyuJyLxoaRLROoyDTgyeDwIWAgUmlk7M0sHDgXmmNmtZjbTzBaa2aMWMMDMZoQuFBxhWhA8HmFmn5rZ7ODoWeeaN67tHDP7xMzuMrMZZrbczMYHn29tZi+Z2WIze83MvjSzkWZ2J9DKzOaZ2b+Dl082s8fMbJGZvW9mrcJ87ROAOcHu3KH7/sPMZpnZEjMbZWavmtkKM/tTlbh3Vzm+ycwWmNn8YByh69xjZrOA68zseDObGzzvyeD/V8xstZn93szmBF8bEHz+mODXMi/4vrbB270OXBz5H62IJJqSLhGplXNuI1BmZj0IjGpNB74kkIiNBBY450qAB5xzo5xzhwGtgNOdc0uBNDPrHbzc+cCLZpYK3A+c65wbATwJ3F71vhGck+KcGw1cD/wu+NxPgR3OuYHAb4ERwa/hZmCvc26ocy6UlPQDHnTODQIKgHPCfPljgdk1nitxzo0EHgEmAdcAhwGXmVmHGl/DKcBZwBjn3OHAX6q8nBa8zoPAU8D5zrnBQApwdZXztjnnhgMPAzcGn7sRuMY5NxQYD+wNPj8r+FhEGiklXSJSn2kEEq5Q0jW9yuPPg+ccFxxZWkBghCg0BfkSgWSL4O8vAv0JJCofmNk84DdAtxr3rO+cV4O/zwZ6BY/HAS8AOOcWAl/V8TWtcs7NC3ONqjoDW2s8F9r0dwGwyDm3yTm3D/gG6F7j3BOAfznnioIxba/y2ovB3/sHY1kefPw0cHSV88J9nZ8Dd5vZ/wE5oZE4YAvQJczXISKNREzqCUSkWQvVdQ0mML24Dvg5sAv4l5llAA8BI51z68zsNiAj+N4Xgf+a2auAc86tMLPBBBKWI6md1XPOvuDv5UT3ObavynE5gdG5mvay/+uo+b6KGteo8BjHngjPO+DrdM7daWZvAacCn5vZycFRxQz2j3qJSCOkkS4Rqc804HRgu3OuPDhik0NginEa+xOTbWaWCVSuPnTOfU0gYfgt+0d3lgF5ZnYkBKYSzaxqcX6k59T0OXBe8PyBBJLEkNLglKUXS4C+Ht9T1QfA5WbWOhhT+zDnLAN6mVnoPpcCn9Z1UTM72Dm3wDl3FzATGBB86RACSbGINFJKukSkPgsIrFr8osZzO51z25xzBcBjBL7hv0cgEajqReASAlONBGvAzgXuMrP5wDyCKyRDIjknjIcIJGqLgT8Bi4CdwdceBb6qUkgfiXeoPtXniXPuXQLTkbOCU6Q3hjmnGLicwGjgAgIjZo/Uc+nrgwsWvgJKg3ECHAe8FW28IhJ/5pzzOwYRkQYzs2Qg1TlXbGYHAx8C/YMJXLTXfA34pXNuRazijBczmwKc5Zzb4XcsIhKearpEpLloDUwOTiMa8NOGJFxBNxMoqG/USZeZ5QF3K+ESadw00iUiIiKSAKrpEhEREUkAJV0iIiIiCaCkS0RERCQBlHSJiIiIJICSLhEREZEE+P/FmQBkoEKFYAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for ind, ab in enumerate(abundances):\n",
    "    print(ab)\n",
    "    observation1, wlen, flux_nu, p1, t1 = Simulator([torch.ones(1)*2000.,torch.ones(1)*-2.0, torch.ones(1)*2.58,\\\n",
    "                                                    torch.ones(1)*ab[0], torch.ones(1)*ab[1], torch.ones(1)*ab[2],\\\n",
    "                                                    torch.ones(1)*ab[3], torch.ones(1)*ab[4], torch.ones(1)*ab[5],\\\n",
    "                                                    torch.ones(1)*ab[6], torch.ones(1)*ab[7]])\n",
    "    \n",
    "    # Spectrum\n",
    "    \n",
    "    plt.figure(figsize=(10,7))\n",
    "    #plt.title(abund[ind])\n",
    "    plt.plot(wlen/1e-4, flux_nu, label= str(np.round(ab[ind], 1)))\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel('Wavelength (microns)')\n",
    "    plt.ylabel(r'Planet flux $F_\\nu$ (10$^{-6}$ erg cm$^{-2}$ s$^{-1}$ Hz$^{-1}$)')\n",
    "    \n",
    "    # PT profile\n",
    "    \n",
    "#     plt.plot(t1, p1, label= str(np.round(ab, 1)))\n",
    "#     plt.yscale('log')\n",
    "#     plt.ylim([1e2, 1e-6])\n",
    "#     plt.ylabel('')\n",
    "#     plt.xlabel('T (K)')\n",
    "#     plt.ylabel('P (bar)')\n",
    "    \n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    plt.clf()\n",
    "    break\n",
    "    \n",
    "# plt.legend()\n",
    "# plt.show()\n",
    "# plt.clf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:petitRT]",
   "language": "python",
   "name": "conda-env-petitRT-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}