{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import solve_for_masses as em\n", "import mass_loss as ms\n", "from numpy.random import normal\n", "%matplotlib inline\n", "\n", "from importlib import reload" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "reload(ms)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As an example we work with the famous Kepler-36 system, which contains a rocky and gaseous planet. Planet parameters are taken from Carter et al. (2012, Science 337 556). \n", "\n", "To start with we need to make some basic choices. First we need to assume the composition of the solid cores. This uses the Fortney et al. (2007, ApJ 659 1661) mass-radius relation. They can either be iron-rock or water-rock (both the iron fraction and ice fraction cannot be non-zero). Xiron=1/3 will consider an \"Earth-like\" composition of 1/3 iron and 2/3 silicate rock." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "Xiron = 1./3.\n", "Xice = 0." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As discussed in the paper (Owen & Campos Estrada, 2020), one needs to choose a Kelvin-Helmholtz timescale for the H/He atmospheres at which to do the comparision (i.e. at what age is mass-loss important), the answer is very weakly dependent on this value, here we pick the standard value of 100 Myr." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "Tkh_Myr=100." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we input the planetary and stellar parameters and their errors. Any errors in the Carter et al. (2012) paper that are not symmetric we just make symmetric crudely using addition in quadrature. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "K36s_M = 1.071 # stellar mass, solar masses\n", "K36s_Mer = 0.043 # stellar mass error, solar masses\n", "K36s_Teff = 5911 # Stellar effective temperature, K\n", "K36s_Teffer = 66 # stellar effective temperature error, K\n", "K36s_R = 1.626 # stellar radius, solar radii - note we should really use the MS radii but this is and example.\n", "K36s_Rer=0.019 # stellar radius error, solar radii \n", "K36s_age=6800 # stellar age, Myr\n", "K36s_age_er=1000. #stellar age error, Myr\n", "\n", "#radius in earth unit and period in days + errors for planet b\n", "K36b_R = 1.486\n", "K36b_Rer=0.035\n", "K36b_P = 13.83989\n", "K36b_Per = np.sqrt(0.00082**2.+0.00060**2.) # make period errors symmetric\n", "\n", "#radius in earth unit and period in days + errors for planet c\n", "K36c_R = 3.679\n", "K36c_Rer=0.054\n", "K36c_P = 16.23855\n", "K36c_Per = np.sqrt(0.00038**2.+0.00054**2.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A note on the efficiency parameter. The code contains three options for the efficiency parameter in the mass-loss models.\n", "\n", "1. Option 1, is to use a constant efficiency parameter.\n", "2. Option 2, is to use the scaling $\\eta \\propto v_{\\rm esc}^{-2}$, where $v_{\\rm esc}$ is the escape velocity from the planet's surface. This fit was provided in Owen & Wu (2017) as an approximation of the hydrodynamic models. It works well for typical sub-neptunes; however, fails when the escape temperature from the planet is smaller than the outflow temperature. In this case it overestimates the efficiency. This occurs for larger planets, or planet's where the minimum required core-mass is low. \n", "3. Option 3 (Default) is a full fit to the Owen & Jackson (2012) mass-loss rates. It provides a normalised value of the efficiency (scaled to the highest value in the table). For very puffy planets, photoevaporation does not occur, but rather Roche Lobe overflow (grey region in Figure 5 of Owen & Jackson 2012). If this happens the efficiency value is simply extrapolated at a constant value for the last hydrodyanmic one. \n", "4. Additional options - you are free to modify the \"efficiency\" function in mass_loss.py to include your own efficiency function" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# set the efficiency option \n", "eff_option = 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now do the calculation for the minimum mass of Kepler-36c to be consistent with photoevaporation assuming the mean values for all parameters." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The location of the radius valley can be chosen. The default is set to 1.8 Rearth, but this is generally only true for sun-like stars. One should change the location of the valley accordingly (for example for M-dwarfs the valley is generally at smaller radii - see van Eylen et al. 2021). " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "valley_loc=1.8" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "8.412384425520381 Mearth\n" ] } ], "source": [ "# evaluate minimum mass for the mean value and print it\n", "\n", "system = em.psystem('Kepler36')\n", "system.add_planet('36b',K36b_R,K36b_Rer,K36b_P,K36b_Per)\n", "system.add_planet('36c',K36c_R,K36c_Rer,K36c_P,K36c_Per)\n", "system.star.mass=K36s_M\n", "system.star.radius=K36s_R\n", "system.star.Teff=K36s_Teff\n", "system.star.age = K36s_age\n", "system.update_planet_info()\n", "system.above_or_below_valley()\n", "system.mass_rocky(Xiron,Xice)\n", "\n", "ms.find_hardest_rocky(system,Tkh_Myr,Xiron,Xice,eff_option)\n", "\n", "Mout, flag = ms.min_mass_gaseous(system.planets[0],system.planets[1],Tkh_Myr,Xiron,Xice,system.star.age,eff_option)\n", "Mmedian = np.copy(Mout)\n", "\n", "print(Mmedian,'Mearth')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, given the stellar and planetary parameters contain errors these must be included, therefore we assume the errors are gaussian and randonly sample them. Here we use 1000 samples, but a value should be chosen such that the upper-limit is converged. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "Mout_error= np.zeros(1000)\n", "flag_out = np.zeros(1000,dtype=np.int8)\n", "\n", "for i in range(np.size(Mout_error)):\n", " \n", " K36b_R_use = normal(K36b_R,K36b_Rer,1)\n", " K36b_P_use = normal(K36b_P,K36b_Per,1)\n", " \n", " K36c_R_use = normal(K36c_R,K36c_Rer,1)\n", " K36c_P_use = normal(K36c_P,K36c_Per,1)\n", " \n", " system = em.psystem('Kepler36_%d' %i)\n", " system.add_planet('36b',K36b_R_use,K36b_Rer,K36b_P_use,K36b_Per)\n", " system.add_planet('36c',K36c_R_use,K36c_Rer,K36c_P_use,K36c_Per)\n", " system.star.mass=normal(K36s_M,K36s_Mer,1)\n", " system.star.radius=normal(K36s_R,K36s_Rer,1)\n", " system.star.Teff=normal(K36s_Teff,K36s_Teffer,1)\n", " system.star.age = normal(K36s_age,K36s_age_er,1)\n", " system.update_planet_info()\n", " system.above_or_below_valley()\n", " system.mass_rocky(Xiron,Xice)\n", " \n", " ms.find_hardest_rocky(system,Tkh_Myr,Xiron,Xice,eff_option)\n", "\n", " Mout, flag = ms.min_mass_gaseous(system.planets[0],system.planets[1],Tkh_Myr,Xiron,Xice,system.star.age,eff_option)\n", " \n", " Mout_error[i] = Mout\n", " flag_out[i]=flag" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now plot the resulting distribution and calculate the 95% upper-limit, the 95% upperlimit is plotted as a dashed line, the actual measured mass is plotted as the dotted line. This indicates (as known from previous work - Lopez & Fortney 2013, Owen & Morton 2016) that Kepler-36 is consistent with the photoevaporation model. The value to quote, as we are concerned with the minimum mass to be consistent is something like the 95% upper-limit, not the mean or median value. " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The 95% upper limit to be consistent with photoevaporation is\n", "6.930341420699428 Mearth\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAG6CAYAAAAF5Ty4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJXklEQVR4nO3deVhUZfsH8O8szAwDuCR7iqBZKpQFZC6Zmqm5UO6WpYJmGVopaW7llonZm2mumdtrmi+/cinTUjTTTCvcyq1FUzEFES22gRmYOb8/RoYQzJlhDmc48/1c17l6OHPOzD0PNnPznOc8t0IQBAFEREREMqKUOgAiIiIiV2OCQ0RERLLDBIeIiIhkhwkOERERyQ4THCIiIpIdJjhEREQkO0xwiIiISHaY4BAREZHsMMEhIiIi2WGCQ0RERLIjeYKzb98+xMXFITQ0FAqFAlu2bLntOXv37kVMTAx0Oh0aNWqEZcuWiR8oERER1RiSJzgFBQVo0aIFFi1aZNfx586dQ/fu3dGuXTscPXoUkydPxssvv4yNGzeKHCkRERHVFAp3KrapUCiwefNm9OrV65bHTJgwAZ9//jlOnz5t2zdy5Ej89NNPOHjwYDVESURERO5OLXUAjjp48CC6dOlSbl/Xrl2xcuVKFBcXw8vLq8I5RqMRRqPR9rPFYsH169dRr149KBQK0WMmIiKiqhMEAXl5eQgNDYVS+e8XoWpcgpOZmYmgoKBy+4KCglBSUoLs7GyEhIRUOCc5ORkzZsyorhCJiIhIRBcvXkT9+vX/9Zgal+AAqDDqUnqV7VajMZMmTUJSUpLt55ycHISFheHixYuoVauWeIGSQwoKChAaGgoAuHz5Mnx8fCSOiIiI3Elubi4aNGgAPz+/2x5b4xKc4OBgZGZmltuXlZUFtVqNevXqVXqOVquFVqutsL9WrVpMcNyISqWytWvVqsUEh4iIKmXP9BLJ76JyVOvWrZGamlpu386dOxEbG1vp/BsiKs9kAt55x7qZTFJHQ0QkDskTnPz8fBw7dgzHjh0DYL0N/NixY0hPTwdgvbw0ZMgQ2/EjR47EhQsXkJSUhNOnT2PVqlVYuXIlxo0bJ0X4RDVOcTHw2mvWrbhY6miIiMQh+SWqQ4cOoWPHjrafS+fKDB06FGvWrEFGRoYt2QGAiIgIbN++HWPHjsXixYsRGhqK999/H3379q322Mm11Go1hg4damuTONRq4EY3g91MRHLlVuvgVJfc3FzUrl0bOTk5nINDRERUQzjy/S35JSoiIiIiV+MANbkNQRBgMBgAAHq9noswEhGR0ziCQ27DYDDA19cXvr6+tkSHXK+gAKhTx7oVFEgdDRGRODiCQ+SBcnKkjoCISFxMcIg8jLc38NtvZW0iIjligkPkYZRKoEkTqaMgIhIX5+AQERGR7HAEh8jDFBcDy5db288/D7DCCRHJERMcIg9jMgGjR1vb8fFMcIhInpjgkNtQqVTo16+frU3iUKmAG90MdjMRyRVLNbBUAxERUY3AUg1ERETk0ZjgEBERkewwwSG3UVBQAIVCAYVCgQLWEBCNwQDcead1Y0UMIpIrTjIm8jCCAFy+XNYmIpIjJjhEHkanA44eLWsTEckRExwiD6NSAfffL3UURETi4hwcIiIikh2O4BB5mOJiYP16a/uZZ7iSMRHJExMcIg9jMgEJCdZ2//5McIhInpjgkNtQqVTo3r27rU3iUKmAG93MUg1EJFss1cBSDURERDUCSzUQERGRR2OCQ0RERLLDBIfcRkFBAXx8fODj48NSDSIyGIAmTawbSzUQkVxxkjG5FQO/cUUnCMCZM2VtIiI5YoJD5GF0OmD//rI2EZEcMcEh8jAqFdC2rdRREBGJi3NwiIiISHY4gkPkYUpKgM2bre3evQE1PwWISIb40UbkYYxGYMAAazs/nwkOEckTP9rIbSiVSrRv397WJnEolcCNbga7mYjkiqUaWKqBiIioRmCpBiIiIvJoTHCIiIhIdpjgkNsoKChAQEAAAgICWKpBRIWFwP33W7fCQqmjISISBycZk1vJzs6WOgTZs1iAn34qaxMRyRETHCIPo9MBO3eWtYmI5IgJDpGHUamAzp2ljoKISFycg0NERESywxEcIg9TUgLs2GFtd+3KlYyJSJ740UbkYYxGoGdPa5ulGohIrvjRRm5DqVQiNjbW1iZxKJXAjW5mqQYiki0mOOQ2vL29kZaWJnUYsuftDbCbiUju+PcbERERyQ4THCIiIpIdJjjkNgwGA8LDwxEeHg6DwSB1OLJVWAi0bWvdWKqBiOSKc3DIbQiCgAsXLtjaJA6LBThwoKxNRCRHTHCIPIxWC2zeXNYmIpIjJjhEHkatBnr1kjoKIiJxcQ4OERERyQ5HcIg8jNkMfPuttd2unbX4JhGR3DDBIfIwRUVAx47Wdn4+4OMjbTxERGJggkNuQ6FQoHnz5rY2iUOhAG50M9jNRCRXTHDIbej1epw8eVLqMGRPrwfYzUQkd5xkTERERLLDBIeIiIhkhwkOuQ2DwYDIyEhERkayVIOICguBzp2tG0s1EJFccQ4OuQ1BEHDq1Clbm8RhsQC7dpW1iYjkiAkOkYfRaoF168raRERyxASHyMOo1cAzz0gdBRGRuDgHh4iIiGSHIzhEHsZsBo4csbajo1mqgYjkyS1GcJYsWYKIiAjodDrExMTg29JCObewfv16tGjRAnq9HiEhIUhISMC1a9eqKVqimq2oCGjZ0roVFUkdDRGROCRPcFJSUjBmzBhMmTIFR48eRbt27dCtWzekp6dXevz+/fsxZMgQDB8+HCdPnsQnn3yCtLQ0PPfcc9UcObmaQqFAw4YN0bBhQ5ZqEJFCATRsaN3YzUQkVwpB4vtxH3roIURHR2Pp0qW2fc2aNUOvXr2QnJxc4fj//Oc/WLp0Kc6ePWvbt3DhQsydOxcXL1606zVzc3NRu3Zt5OTkoFatWlV/E0RERCQ6R76/JR3BMZlMOHz4MLp06VJuf5cuXXDgwIFKz2nTpg3+/PNPbN++HYIg4MqVK/j000/Ro0ePW76O0WhEbm5uuY2IiIjkS9IEJzs7G2azGUFBQeX2BwUFITMzs9Jz2rRpg/Xr12PgwIHQaDQIDg5GnTp1sHDhwlu+TnJyMmrXrm3bGjRo4NL3QURERO5F8jk4ACrMtxAE4ZZzME6dOoWXX34ZU6dOxeHDh/HVV1/h3LlzGDly5C2ff9KkScjJybFt9l7KoupVWFiIBx98EA8++CAKWUNANEVFQK9e1o2TjIlIriS9Tdzf3x8qlarCaE1WVlaFUZ1SycnJaNu2LcaPHw8AuO++++Dj44N27dph1qxZCAkJqXCOVquFlku2uj2LxYJDhw7Z2iQOsxn47LOyNhGRHEk6gqPRaBATE4PU1NRy+1NTU9GmTZtKzzEYDFAqy4eturGQB+sXEd2eRgMsX27dNBqpoyEiEofkC/0lJSVh8ODBiI2NRevWrbF8+XKkp6fbLjlNmjQJly5dwtq1awEAcXFxGDFiBJYuXYquXbsiIyMDY8aMQcuWLREaGirlWyGqEby8gBEjpI6CiEhckic4AwcOxLVr1zBz5kxkZGQgKioK27dvR8OGDQEAGRkZ5dbEiY+PR15eHhYtWoRXX30VderUwaOPPoq3335bqrdAREREbkbydXCkwHVw3FNBQQF8fX0BAPn5+fDx8ZE4InmyWIDTp63tZs0ApVvcakBEdHuOfH9LPoJDRNWrsBCIirK28/MB5pFEJEdMcMit+Pv7Sx2CR2A3E5HcMcEht+Hj44OrV69KHYbs+fgA7GYikjtefSciIiLZYYJDREREssMEh9xGYWEhOnTogA4dOrBUg4iKioBnnrFuLNVARHLF28R5m7jb4G3i1aOgALjRzbyLiohqFN4mTkS3pNEA771X1iYikiMmOEQexssLGDNG6iiIiMTFOThEREQkOxzBIfIwFgtQWt4tLIylGohInpjgEHmYwkIgIsLa5iRjIpIrJjjkVvR6vdQheAR2MxHJHRMcchs+Pj4oKCiQOgzZ8/Gx3ipORCRnvPpOREREssMEh4iIiGSHCQ65jaKiIvTo0QM9evRAEWsIiMZoBEaMsG5Go9TREBGJg6UaWKrBbbBUQ/VgqQYiqqlYqoGIbsnLC5g1q6xNRCRHTHCIPIxGA0yZInUURETi4hwcIiIikh2O4BB5GEEAsrOtbX9/QKGQNh4iIjEwwSHyMAYDEBhobXOSMRHJFS9RERERkexwBIfcho+PDzxw1YJq5+NjvUxFRCRnHMEhIiIi2WGCQ0RERLLDBIfcRlFREfr374/+/fuzVIOIjEZgzBjrxlINRCRXDpdqeOmllzB69Gjcc889YsUkOpZqcE8s1VA9WKqBiGoqR76/HR7BWbt2LZo3b47OnTvjs88+46RQohrGywuYPNm6sVQDEcmVwwnO5cuXsXDhQmRkZKB3794IDw/HnDlzkF26chgRuTWNBnjrLeum0UgdDRGROBxOcHx8fJCYmIgTJ05g165diImJwRtvvIEGDRogPj4ehw4dEiNOIiIiIrtVaZLxo48+ik2bNuHcuXNo06YNPvroIzz00EN46KGHsHXrVlfFSEQuJAjWeTgFBVwPh4jkq0oJTmFhIVasWIG4uDjs2bMHzZo1w7Rp02A2m9GrVy+8+eabroqTiFzEYLBOMvb1tbaJiOTIqQTn7NmzSEpKwp133omRI0eifv362LlzJ06cOIGpU6fi0KFDmDBhAhYuXOjqeImIiIhuy+FSDd26dUNqaip8fHyQkJCAl156CY0bN65wXFxcHObMmeOSIMkz6PV65Ofn29okDr3eent4aZuISI4cTnDOnj2L9957DwkJCbY1SyoTFRWFPXv2VCk48iwKhYJr31QDhYJr3xCR/Dmc4Pz22292Hefn54f27ds7HBARERFRVTk8B0elUuHHH3+s9LHDhw9DpVJVOSjyTEajEfHx8YiPj4eRNQREYzIBU6ZYN5NJ6miIiMThcILzbysXWywWKBSKKgVEnqukpAT//e9/8d///hclJSVShyNbxcXA7NnWrbhY6miIiMTh8CUqALdMYg4fPozatWtXKSAiEpdaDbzySlmbiEiO7Pp4W7BgARYsWADAmtz06tULWq223DGFhYXIyspCv379XB8lEbmMVgvMny91FERE4rIrwQkMDERkZCQA4Pz582jUqBHq1KlT7hitVot7770Xr5T+aUhEREQkEYXgYDnwjh07YunSpWjatKlYMYnOkXLrVH0KCgpsSw/k5+fzlnEiIirHke9vh6/Ac20bopqtoMBapgGwLvjHPJKI5MiuBCc9PR0hISHw8vJCenr6bY8PCwurcmBEREREzrIrwYmIiMDBgwfRsmVLhIeH3/ZWcLPZ7JLgyLPo9XpkZWXZ2iQOvR640c0s1UBEsmVXgrNq1SpbvalVq1ZxrRsShUKhQEBAgNRhyJ5CAbCbiUjuHJ5kLAecZExERFTzOPL97fBKxpUpKirCL7/8wktTVCVGoxGjRo3CqFGjWKpBRCYT8NZb1o2lGohIrhwewVm4cCH+/vtvvPHGGwCsqxc//vjjuH79OsLDw/HNN9+gQYMGogTrKhzBcU+8Tbx68C4qIqqpRB3BWbFiRblF/iZMmIA77rgD7733HgRBwKxZsxwOmIiqj1oNPPecdWOpBiKSK4c/3tLT022L/OXl5WHfvn343//+hz59+qBu3bqYOnWqy4MkItfRaoEPP5Q6CiIicTk8gmM0GuHl5QUAOHjwICwWCx577DEAQHh4ODIzM10bIREREZGDHE5wwsLC8O233wIAPvvsM9x///2262BXr17lnBYiIiKSnMMJzrPPPouZM2ciJiYGH3zwAZ599lnbY4cOHcLdd9/t0gCJyLUKCqwTi318rG0iIjlyeA7OlClToFarceDAAfTu3Rsvv/yy7bETJ06gb9++Lg2QiFzPYJA6AiIicXGhP15ScxsWi8VW6ywsLAxKpUuWaaKbWCxAaUm5sDCA3UxENYWo1cSJxKJUKhEeHi51GLKnVALsZiKSO6cSnHXr1uHjjz/GhQsXUFhYWO4xhUKBs2fPuiQ4IiIiImc4nOC8/fbbmDRpEpo3b44WLVpAq9WKERd5IJPJhClTpgAA3nrrLWg0GokjkqfiYmDxYmt71CjgxqoPRESy4vAcnMaNG6N79+5YuHChWDGJjnNw3BNLNVQPlmogoppK1Dk4mZmZ6N27t9PBEZG0VCpg0KCyNhGRHDl8/0RMTIzL59gsWbIEERER0Ol0iImJsS0keCtGoxFTpkxBw4YNodVq0bhxY6xatcqlMRHJlU4HrF9v3XQ6qaMhIhKHwyM48+bNw7PPPovo6GjExMRUOYCUlBSMGTMGS5YsQdu2bfHBBx+gW7duOHXqFMLCwio9Z8CAAbhy5QpWrlyJu+66C1lZWSgpKalyLERERCQPDs/Buffee5GZmYnr168jODgY9erVK/+ECgV++uknu5/voYceQnR0NJYuXWrb16xZM/Tq1QvJyckVjv/qq6/w1FNP4Y8//sAdd9zhSOg2nIPjnjgHh4iI/o2oc3Dq1asHf39/p4P7J5PJhMOHD2PixInl9nfp0gUHDhyo9JzPP/8csbGxmDt3Lj766CP4+PjgiSeewJtvvglvb+9KzzEajTAajbafc3NzXRI/UU1UUFC2Ds7585xkTETy5HCC880337jsxbOzs2E2mxEUFFRuf1BQ0C2rkv/xxx/Yv38/dDodNm/ejOzsbCQmJuL69eu3nIeTnJyMGTNmuCxuopouO1vqCIiIxOUWKxkrFIpyPwuCUGFfKYvFAoVCgfXr16N27doArPOC+vXrh8WLF1c6ijNp0iQkJSXZfs7NzUWDBg1c+A7IFby9vXHixAlbm8Th7Q3c6Gawm4lIrpyqQnP16lVMmjQJrVu3RpMmTXDy5EkAwAcffICjR4/a/Tz+/v5QqVQVRmuysrIqjOqUCgkJwZ133mlLbgDrnB1BEPDnn39Weo5Wq0WtWrXKbeR+lEolIiMjERkZyTpUIlIqgchI68ZuJiK5cvjj7dy5c2jRogXef/99KBQK/PHHH7b5LT///DPef/99u59Lo9EgJiYGqamp5fanpqaiTZs2lZ7Ttm1bXL58Gfn5+bZ9v/32G5RKJerXr+/o2yEiIiIZcjjBee2111CnTh38/vvv2LdvH/55E9bDDz+M7777zqHnS0pKwooVK7Bq1SqcPn0aY8eORXp6OkaOHAnAenlpyJAhtuMHDRqEevXqISEhAadOncK+ffswfvx4DBs2jJc1ajiTyYTp06dj+vTpMJlMUocjW8XFwIcfWrfiYqmjISISh8NzcHbv3o2lS5ciNDQUZrO53GMhISG4fPmyQ883cOBAXLt2DTNnzkRGRgaioqKwfft2NGzYEACQkZGB9PR02/G+vr5ITU3FSy+9hNjYWNSrVw8DBgzArFmzHH0r5GaKi4ttk8HHjx/PWlQiMZmA55+3tgcNYi0qIpInhxOcoqKiW64/U1BQ4NTcicTERCQmJlb62Jo1ayrsa9q0aYXLWkRkH5UKePLJsjYRkRw5nI3cc8892LVrV6WP7du3D1FRUVUOiojEo9MBW7ZYN5ZqICK5cngEZ8SIEUhKSkJoaCieeeYZANa5E59++imWLFmCRYsWuTxIIiIiIkc4XKoBAJ5//nmsWLECSqUSFosFSqUSgiBgxIgRWLZsmRhxuhRLNbgnlmogIqJ/I2qpBgBYvnw5hg0bhm3btuHKlSvw9/dHz549b3lrNxG5D4MBaN7c2j51CtDrpY2HiEgMTq9k3KpVK7Rq1cqVsRBRNRAE4MKFsjYRkRy5RakGIgDQ6XT48ccfbW0Sh04H3OhmTjImItmyK8FRKpW3rA1VmZvXxyGyh0qlwoMPPih1GLKnUgHsZiKSO7sSnKlTp5ZLcFavXo38/HzExcUhODgYGRkZ+OKLL+Dj44Nhw4aJFiwRERGRPexKcKZPn25rv/vuuwgODsauXbtsd7wAQF5eHh577DHoOWORnGQymbBgwQIAwCuvvMKVjEVSUgKkpFjbAwcCal6oJiIZcvg28caNG+Odd95Bnz59Kjy2ceNGjBs3DufOnXNZgGLgbeLuibeJV4+CAqD0b5P8fIDdTEQ1hai3iV+6dAnqW/zJp1arkZmZ6ehTElE1UiqBxx4raxMRyZHDH2/NmjXDvHnzUHxTGWKTyYR3330XTZs2dVlwROR63t5Aaqp18/aWOhoiInE4PIIza9Ys9OrVC40aNUKfPn0QHByMzMxMbNq0CZmZmdiyZYsIYRIRERHZz+EEp0ePHvjqq68wZcoULF68GBaLBQqFAi1btsTq1avxWOnYNxEREZFEnLp/olOnTujUqRMMBgP++usv1K1bl3dPEdUQBkPZOjhpaSzVQETyVKUbRPV6PRMbohpGEKw1qErb/yZu4X5czTNW6fUC/LTY+tLDVXoOIiJHcQUMchs6nQ579uyxtUkcOh1wo5tvW6rhap4RmblF4gdFRORiTHDIbahUKnTo0EHqMGRPpQIc7WalAgj0cyzpzMorgoXFPIlIIkxwiOi2Av10+H5yJ4fOaTV7N0d/iEgyTHDIbRQXF2P58uUAgOeffx5eXl4SRyRPJSXAF19Y2z17slQDEcmTwx9tmZmZCA4OFiMW8nAmkwmjR48GAMTHxzPBEYnRCPTubW3n5zPBISJ5cngl47CwMDz99NP47rvvxIiHiESmVAJt2lg3lmogIrly+OPt9ddfx7fffotHHnkE999/P1auXInCwkIxYiMiEXh7A999Z91YqoGI5MrhBGfq1Km4cOECNmzYgFq1amHEiBGoX78+xo0bh7Nnz4oRIxEREZFDnBqgVqlUGDBgAPbt24djx46hb9++WLZsGe655x707NkTO3bscHWcRERERHar8hX4e++9F926dUNUVBQsFgt2796N7t27IzY2Fr/99psrYiQiFyostJZqePBBa5uISI6cTnCys7ORnJyMiIgI9OvXD2q1GikpKcjNzcWWLVuQl5eH+Ph4F4ZKRK5gsQCHDlk3i0XqaIiIxOHwDaI//PADFi9ejE8++QSCIGDgwIF45ZVXEB0dbTsmLi4OarUavXr1cmWsJHNarRZf3FigRavVShyNfGm1ZevgsJuJSK4cTnBat26N4OBgTJw4ES+++CICAwMrPS48PBxt2rSpcoDkOdRqNXr06CF1GLKnVgPsZiKSO4cTnLVr12LgwIG3XYStWbNmtsKJRERERNXJ4Tk4f/zxB65evVrpYxkZGZg5c2aVgyLPVFxcjDVr1mDNmjUoLi6WOhzZMpuB1FTrZjZLHQ0RkTgcTnBmzJiBP//8s9LHLl++jBkzZlQ5KPJMJpMJCQkJSEhIgMlkkjoc2SoqArp0sW5FrIVJRDLl8CUqQRBu+Vh+fj7rBxG5OaUSaNGirE1EJEd2JTg///wzjh07Zvt5+/bt+OWXX8odU1hYiPXr16Nx48YuDZCIXMvbG/jH/85ERLJkV4KzefNm26UnhUJxy3k23t7eWL16teuiIyIiInKCXQnO888/j549e0IQBLRs2RKrV69GVFRUuWO0Wi0aN24Mb1bvIyIiIonZleCEhIQgJCQEALBnzx5ER0fDz89P1MCISByFhUC3btb2l1+yojgRyZPDk4zbt28vRhxEVE0sFmDv3rI2EZEc2ZXgDBs2DG+88QYiIiIwbNiwfz1WoVBg5cqVLgmOPItWq8X//d//2dokDq0WuNHNLNVARLJlV4KzZ88evPLKKwCAr7/+GgqF4pbH/ttjRP9GrVajf//+Uoche2o1wG4mIrmzK8E5d+6crX3+/HmxYiEiIiJyCYfn4BCJpaSkBJs3bwYA9O7dG2o1/3mKwWwGvv/e2m7VClCppI2HiEgM/AYht2E0GjFgwAAA1lWxmeCIo6gIePhhazs/H/DxkTYeIiIx2PUNEhERYffcGoVCgbNnz1YpKCISj0IB3HVXWZuISI7sSnDat2/PycNEMqHXA7//LnUURETisivBWbNmjchhEBEREbkOawkTERGR7Ng1gpOeno6QkBB4eXkhPT39tseHhYVVOTAiEkdREdC3r7W9cSOg00kbDxGRGOyeZHzw4EG0bNkS4eHht52PYzabXRIcEbme2Qxs317WJiKSI7sSnFWrVqFx48a2Nicckxg0Gg1Wr15ta5NrxC3cj6t5RtvPFrMCjfoGAwA6zsuEUiXc8tysvCLR4yMiEoNdCc7QoUNt7fj4eLFiIQ/n5eXFf18iuJpnRGbuTYnKXdbVybMKJAiIiKgaVGklNUEQkJ+fD19fX47qELk5pQII9HNuwk2AH6tyElHN4lSC88MPP2Dq1KnYt28fTCYTNBoNHnnkEcyYMQOtWrVydYzkIUpKSrBjxw4AQNeuXbmSsYsF+unw/eROMJuB48et++69l6UaiEieHP4G+frrr9GtWzf4+fnhqaeeQnBwMDIzM7F161a0b98e27dvR6dOncSIlWTOaDSiZ8+eAFiqQUxFRcADD1jbLNVARHLl8DfIhAkT8MADD2DXrl3w9fW17c/Ly0OnTp0wceJEpKWluTRIInIdhQIIDS1rExHJkcMJzokTJ7B+/fpyyQ0A+Pn5YcKECXj22WddFhwRuZ5eD1y6JHUURETicngl48DAQCiVlZ+mUqkQEBBQ5aCIiIiIqsLhBOeFF17Ae++9h+Li4nL7TSYT5s2bh+eff95lwRERERE5w65LVPPmzbO1NRoNzp8/j0aNGqFPnz62ScabNm2CSqWCt7e3aMESUdUVFQGDB1vbH33EUg1EJE92JTjjxo2rdP/ChQsr7Hvttdfw6quvVi0qIhKN2Qx8+qm1vWaNpKEQEYnGrgTn3LlzYsdBBI1Gg0WLFtnaJA6NBrjRzWA3E5Fc2ZXgNGzYUOw4iODl5YVRo0ZJHYbseXkB7GYikjuupEZEbufmAqHOCPDTYutLD7soIiKqaZxKcPbt24f3338fp0+fRmFhYbnHFAoFzp4969DzLVmyBO+88w4yMjIQGRmJ+fPno127drc977vvvkP79u0RFRWFY8eOOfSa5H7MZjO+/fZbAEC7du2gYg0BUVgsQOn/oo0bA7dY9UFSlRYIJSJygMMJzv79+9GpUyd06NABp0+fxuOPP468vDwcPHgQjRo1Qtu2bR16vpSUFIwZMwZLlixB27Zt8cEHH6Bbt244deoUwsLCbnleTk4OhgwZgk6dOuHKlSuOvg1yQ0VFRejYsSMAa6kGH9YQEEVhIXD33da2u5dqcKZAaFZeESyC9b+tZu926nU5+kNU8zmc4EybNg0JCQlYunQpvLy8MGvWLERHR+Pnn3/G448/jj59+jj0fPPmzcPw4cPx3HPPAQDmz5+PHTt2YOnSpUhOTr7leS+88AIGDRoElUqFLVu2OPo2iDxa7dpSR2Cf0gKhjmg1ezcyc61JDkeBiDyXU6Uaxo0bB8WNIjZmsxkAcN999+GNN97AzJkzERcXZ9dzmUwmHD58GBMnTiy3v0uXLjhw4MAtz1u9ejXOnj2LdevWYdasWbd9HaPRCKOx7Hp+bm6uXfERyZGPD/D331JHIZ4AP63T55aO/hBRzedwgmMwGODr6wulUgmtVovs7GzbY02bNsWpU6fsfq7s7GyYzWYEBQWV2x8UFITMzMxKz/n9998xceJEfPvtt3ZXm05OTsaMGTPsjouIaq6qXFoqHf0hoprP4emFYWFhtjkvzZs3x7Zt22yP7d27F/Xq1XM4CMVNJY0FQaiwD7COFg0aNAgzZszA3aWTCOwwadIk5OTk2LaLFy86HCMRERHVHA6P4HTo0AHffPMN+vXrhxEjRiAxMRGnT5+GVqvFzp07HVrF2N/fHyqVqsJoTVZWVoVRHQDIy8vDoUOHcPToUYwePRoAYLFYIAgC1Go1du7ciUcffbTCeVqtFlqt88PWRHJiNAIvvGBtf/ABwP81iEiOHE5wZsyYgevXrwMARo4cCYPBgPXr10OhUOD111/HlClT7H4ujUaDmJgYpKamonfv3rb9qampePLJJyscX6tWLRw/frzcviVLluDrr7/Gp59+ioiICEffDpHHKSkB/vtfa3vxYiY4RCRPDic4/v7+8Pf3t/2clJSEpKQkpwNISkrC4MGDERsbi9atW2P58uVIT0/HyJEjAVgvL126dAlr166FUqlEVFRUufMDAwOh0+kq7Keax8vLC3PnzrW1SRxeXsCNbga7mYjkqkorGV++fBnXrl1DvXr1EBoa6tRzDBw4ENeuXcPMmTORkZGBqKgobN++3VYeIiMjA+np6VUJk2oIjUaD8ePHSx2G7Gk0QHV2szPr0WTlcaIvEVWNUwnOpk2bMGnSJJw5c8a2r3Hjxpg9ezb69evn8PMlJiYiMTGx0sfW3Kbc8fTp0zF9+nSHX5OIqgfXoyEiKTic4KSkpODpp59G06ZNMXXqVAQHByMjIwMpKSkYOHAgPv74YwwcOFCMWEnmzGYzjhw5AgCIjo5mqQaRWCxARoa1HRIiXqmGqqxH48rnICLPpBAEwaFlrSIjIxEeHo6tW7dC+Y9PRovFgh49eiA9PR0nT550eaCulJubi9q1ayMnJwe1atWSOhy6oaCgAL6+vgBYqsGVStd2Ca5lXRW4oAC40c1uX6qhut3cV0TkXhz5/nb4b7ezZ88iMTGxXHIDAEqlEomJiQ4X2iSi6qdWWzciIrly+COuYcOGMBgMlT5mMBjQoEGDKgdFROLx8QGKi6WOgohIXA6P4Lz66quYOXNmuRINgHVxvlmzZmHcuHEuC46IiIjIGXaN4Lz88svlfs7NzUV4eDg6deqE4OBgZGZmYvfu3fD393eoFhURERGRGOxKcBYtWlTp/q1bt5b7OT09HYsWLcKCBQuqHhkRicJoBErX5pw3jysZE5E82ZXgWCwWseMgompSUgIsWWJtz53LBIeI5In3UZDb8PLywrRp02xtEoeXF3Cjm1mqgYhky+kEZ/fu3di9ezeuXbsGf39/dOrUqdJK3kT20mg0XJW6Gmg0ALuZiOTO4QTHZDKhb9++2L59OwRBgFqtRklJCebMmYMePXpg48aN/OubiIiIJOXwbeIzZ87Ejh07MGfOHFy5cgUmkwlXrlzB22+/jR07dmDmzJlixEkewGKx4OTJkzh58iTnfYlIEIC//7Zujq1jTkRUczg8grNhwwZMnjy5XNXngIAAjBs3Dvn5+Vi7di3efPNNlwZJnqGwsBBRUVEAWKpBTAYDULeutc1SDUQkVw6P4Pz5559o165dpY+1a9cOly5dqnJQRERERFXhcIITEBCA48ePV/rY8ePHERAQUOWgiEg8ej1gMlk3vV7qaIiIxOHwJaonnngCU6dORVhYGPr06WPb/9lnn2H69Ol45plnXBogkVzELdyPq3nGKj1HgJ8WW196uErPoVDw9nAikj+HE5y33noL3333Hfr37w8fHx8EBwfjypUryM/Px7333ou33npLjDiJaryreUZk5hZJHQYRkUdwOMGpW7cufvzxR6xZswZ79uzBtWvXEB0djU6dOmHIkCHQcllUon+lVACBfjqHzsnKK4LFRXc8mUzAlCnW9ltvWdfFISKSG4cSnMLCQjz22GOYMWMGXnjhBbzwwgtixUUkW4F+Onw/uZND57Savdtloz/FxcB//mNtT5/OBIeI5MmhBMfb2xvHjx+HWs0KD+R6Xl5eGDdunK1N4vDyAm50M+fiEJFsOZyptG7dGj/++CM6dOggQjjkyTQaDd555x2pw5A9jQZgNxOR3Dmc4Lz77rt48sknERwcjD59+sDX11eMuIjcUlXuhMrK4wRjIqLq4tQIjslkQkJCAhISEqDX66FQKGyPKxQK5OTkuDRI8gwWiwXp6ekAgLCwMCiVDi/TJDo53AklCEBJibWtVltvGycikhuHE5y+ffuWS2iIXKWwsBAREREA3L9UgzN3QpUK8JP2TkODASgdeGWpBiKSK4cTnDVr1ogQBlHN4sydUEREVH3sTnAKCwuxZcsWXLhwAYGBgYiLi2NZBqIaSK8H/vqrrE1EJEd2JTiXL1/GI488gnPnzkEQrKuN1a5dG19++SVatWolaoBE5FoKBVCnjtRREBGJy64E5/XXX8elS5fw+uuvo1WrVvj999/x1ltv4cUXX8TRo0fFjpGIbsjKK0Kr2bsdPoeIyNPYleCkpqZi8uTJeOONNwAA3bp1Q+PGjfHEE0/gypUrCAoKEjVIIrKyCKjyXVwmEzB7trU9eTJXMiYiebIrwcnMzMQjjzxSbl+HDh0gCAITHKJq4Io7r0qfo7gYmDHDum/8eCY4RCRPdiU4ZrMZ3t7e5fbpdNZbZEtKF9QgqiK1Wo3ExERbm8psfelhlz2XWg3c6Gawm4lIruz+ePv111/LfemYzWYAwC+//FLh2OjoaBeERp5Gq9Vi8eLFUoche1otwG4mIrmzO8GJj4+vdP/gwYNtbUEQoFAobMkPEVFN5Mxk7lIBflqXjrgRkXPsSnBWr14tdhxEEAQB2dnZAAB/f3+umE2SccVkbiKSll0JztChQ8WOgwgGgwGBgYEA3L9UQ01WUFC2Ds7ff7NUwz9VZTJ3Vl4RLIILgyGiKuEUQyIPxHsDKleVS0utZu/mqA+RG2GCQ+RhvL2BP/8saxMRyRETHCIPo1QCd94pdRREROJSSh0AERERkatxBIfIw5hMwIIF1vYrr3AlYyKSJyY4RB6muBh47TVrOzGRCQ4RyRMTHHIbarXatiQBSzWIR60GSld+YDcTkVzx443chlarxZo1a6QOQ/a0WoDdTERyx0nGREREJDscwSG3IQgCDAYDAECv17NUAxEROY0jOOQ2DAYDfH194evra0t0yPVKSzXUqWNtExHJEUdwiDxQTo7UEcgXK5ETuQcmOEQextsb+O23sja5FiuRE7kHJjhEHkapBJo0kToK+WElciL3wgSHiMgFWImcyL0wwSHyMMXFwPLl1vbzzwNeXtLGQ0QkBiY4RB7GZAJGj7a24+OZ4BCRPDHBIbehUqnQr18/W5vEoVIBN7oZ7GYikismOOQ2dDodPvnkE6nDkD2dDmA3E5HccaE/IiIikh0mOERERCQ7THDIbRQUFEChUEChUKCANQREYzAAd95p3VgRg4jkinNwiDyMIACXL5e1iYjkiAkOkYfR6YCjR8vaRERyxASHyMOoVMD990sdBRGRuDgHh4iIiGSHIzhEHqa4GFi/3tp+5hmuZExE8sQEh8jDmExAQoK13b8/Exwikie3uES1ZMkSREREQKfTISYmBt9+++0tj920aRM6d+6MgIAA1KpVC61bt8aOHTuqMVoSi0qlQvfu3dG9e3eWahCRSgV0727d2M1EJFeSJzgpKSkYM2YMpkyZgqNHj6Jdu3bo1q0b0tPTKz1+37596Ny5M7Zv347Dhw+jY8eOiIuLw9HS20KoxtLpdNi2bRu2bdsGHW/vEY1OB2zbZt3YzUQkV5InOPPmzcPw4cPx3HPPoVmzZpg/fz4aNGiApUuXVnr8/Pnz8dprr+HBBx9EkyZNMHv2bDRp0gRbt26t5siJiIjIXUma4JhMJhw+fBhdunQpt79Lly44cOCAXc9hsViQl5eHO+6445bHGI1G5ObmltuIiIhIviRNcLKzs2E2mxEUFFRuf1BQEDIzM+16jnfffRcFBQUYMGDALY9JTk5G7dq1bVuDBg2qFDeJo6CgAD4+PvDx8WGpBhEZDECTJtaNpRqISK4kv0QFAAqFotzPgiBU2FeZDRs2YPr06UhJSUFgYOAtj5s0aRJycnJs28WLF6scM4nDYDDAwG9dUQkCcOaMdWOpBiKSK0lvE/f394dKpaowWpOVlVVhVOdmKSkpGD58OD755BM89thj/3qsVquFVqutcrxEcqDTAfv3l7XJfWTlFaHV7N1Onx/gp8XWlx52YURENZekCY5Go0FMTAxSU1PRu3dv2/7U1FQ8+eSTtzxvw4YNGDZsGDZs2IAePXpUR6hEsqFSAW3bSh0FVcYiAJm5RVKHQSQLki/0l5SUhMGDByM2NhatW7fG8uXLkZ6ejpEjRwKwXl66dOkS1q5dC8Ca3AwZMgQLFixAq1atbKM/3t7eqF27tmTvg4jIWQF+VRthzsorgoWXG4nKkTzBGThwIK5du4aZM2ciIyMDUVFR2L59Oxo2bAgAyMjIKLcmzgcffICSkhKMGjUKo0aNsu0fOnQo1qxZU93hE9U4JSXA5s3Wdu/egFryTwGq6mWlVrN3c+SH6CZu8dGWmJiIxMTESh+7OWn55ptvxA+ISMaMRqD0psP8fCY4RCRP/Ggjt6FUKtG+fXtbm8ShVAI3uhnsZiKSKyY45Da8vb05QlcNvL0BdjMRyR3/fiMiIiLZYYJDREREssNLVOQ2CgoKEB4eDgA4f/48fHx8RHmduIX7cTXP6NS5WXk1/06VwkKgdWtr++BB6yUrIiK5YYJDbiU7O1v017iaZ/ToW2otFuCnn8raRERyxASHPJZSAQT6OVeroKoLs0lJpwN27ixrExHJERMc8liBfjp8P7mT1GFUO5UK6NxZ6iiIiMTFBIckc/NcGLOp0Nbu8M4eqDS3nxzC4oJERFQZJjgkmZvnwlhMZcnOlVwjlBqFFGHJXkkJsGOHtd21K1cyJiJ54kcbSa50LozZJODijX1BtbRQaW49QYTFBZ1nNAI9e1rbLNVARHLFjzaSXOlcmMLCQjyyLRYAsG9CJ3j/y/3LLC7oPKUSiI0taxMRyRETHHIb3t7eSEtLkzoM2fP2BtjNRCR3/PuNiIiIZIcJDhEREckOExxyGwaDAeHh4QgPD4fBYJA6HNkqLATatrVuhYW3P56IqCbiHBxyG4Ig4MKFC7a2PbLyitBq9m6HXkcO9aSqwmIBDhwoaxMRyRETHKrRLAJ4N5WDtFpg8+ayNhGRHDHBoRrJFbWganI9qapQq4FevaSOgohIXExwqEZieQYiIvo3THCIPIzZDHz7rbXdrp21+CbJgzNz0kqxrhvJDRMcIg9TVAR07Ght5+cDPj7SxkOuwzlpRGWY4JDbUCgUaN68ua1N4lAogBvdDHazPFRlPllpXTeO/pDcMMEht6HX63Hy5Empw5A9vR5gN8tLVZKL0rpuHP0huWGCQ0TkwVwx+kPkjpjgEBF5MFeM/hC5I5ZqILdhMBgQGRmJyMhIlmoQUWEh0LmzdWOpBiKSK47gUJXELdyPq3lGp869uWSCIAg4deqUrU3isFiAXbvK2kRSqMpnRylObqZ/wwSHquRqnpFD1DWMVgusW1fWJpICPztIbExwyCWUCiDQT+fUuZ5aMkEqajXwzDNSR0Fk5cxnB29tJ3swwSGXCPTT4fvJnaQOg4hqGGc+O3hrO9mDCQ6RhzGbgSNHrO3oaJZqoJqHt7aTPZjgEHmYoiKgZUtrm6UaqCbire1kDyY45DYUCgUaNmxoa5M4FArgRjezVAO5hDNzYW6+i5LI1ZjgkNvQ6/U4f/681GHInl4PsJvJlTgXhtwRExxy6Vo2ROQ5XHEHJO+iJLEwwSGuR0FETuFt1uTOmOCQjdRr2RQWFuKRRx4BAOzbtw/e3t5Vfk6qqKgIeOopa/t//wN0zv3KiYjcGhMcspF6LRuLxYJDhw7Z2iQOsxn47LOyNhGRHDHBIfIwGg2wfHlZm4hIjpjgEHkYLy9gxAipoyAiEpdS6gCIiIiIXI0jOEQexmIBTp+2tps1A5T8M4eIZIgJDpGHKSwEoqKsbZZqICK5YoJDbsXf31/qEDwCu5mI5I4JDrkNHx8fXL16VeowZM/HB2A3E5Hc8eo7ERERyQ4THCIiIpIdJjjkNgoLC9GhQwd06NABhYWFUocjW0VFwDPPWLciliAjIpniHBxyGxaLBXv37rW1SRxmM/Dxx9Z26YrGRERywwSHyMNoNMB775W1iTxRVl4RWs3e7dS5AX5aVlKvAZjgEHkYLy9gzBipoyCSlkUAMnN5jVbOmOAQEZHHCPDTOn1uVl4RLIILgyFRMcEh8jAWC5Cebm2HhbFUA3mWqlxaajV7N0d9ahAmOEQeprAQiIiwtlmqgYjkigkOuRW9Xi91CB6B3UxEcscEh9yGj48PCgoKpA5D9nx8AHYzEckdr74TERGR7HAERybiFu7H1TyjU+dm5XHSHBERyQsTHJm4mmes8bP7i4qK0LdvXwDAxo0bodPpJI5InoxGYPRoa3vRIkDr/F2zRERuiwmOzCgVQKCfc4lBVdaHcAWz2Yzt27fb2iSOkhJgxQpre/58JjhEjuIqyDUDExyZCfTT4fvJnaQOg9yYlxcwa1ZZm4gcw1WQawYmOG6E82ioOmg0wJQpUkdBVPNItQpyVb4bXKGmjjq5RYKzZMkSvPPOO8jIyEBkZCTmz5+Pdu3a3fL4vXv3IikpCSdPnkRoaChee+01jBw5shojFocc5tEQEcmVVKsg87vBOZInOCkpKRgzZgyWLFmCtm3b4oMPPkC3bt1w6tQphIWFVTj+3Llz6N69O0aMGIF169bhu+++Q2JiIgICAmwTVKXkilGYmjyPhtyfIADZ2da2vz+gUEgbD5EncWb+jiu+G5xROupUU+ccSZ7gzJs3D8OHD8dzzz0HAJg/fz527NiBpUuXIjk5ucLxy5YtQ1hYGObPnw8AaNasGQ4dOoT//Oc/bpHguCLT5jwaEpPBAAQGWtss1UBUvaoyf6e6vxtKR51q6pwjSRMck8mEw4cPY+LEieX2d+nSBQcOHKj0nIMHD6JLly7l9nXt2hUrV65EcXExvCqZNWk0GmE0lo2q5OTkAAByc3Or+hYqKCkqgMVohFIB+Ps6N5pSR10iSmzu7p+rGOfm5vJOKpH8cxXj3FyA3UwkvjrqEpRoqvY/W3V/N1Ql5ux8IywCUFJkdmnMpc8lCLef0CRpgpOdnQ2z2YygoKBy+4OCgpCZmVnpOZmZmZUeX1JSguzsbISEhFQ4Jzk5GTNmzKiwv0GDBlWI/vYuVOHc2hNcFkaNFBoaKnUIHoHdTFSz1LTvhosAar/p+ufNy8tD7dq1//UYyS9RAYDipkkAgiBU2He74yvbX2rSpElISkqy/WyxWHD9+nXUq1fvX19HDLm5uWjQoAEuXryIWrVqVetr12TsN+ew35zDfnMc+8w57DfHCIKAvLw8u/4IljTB8ff3h0qlqjBak5WVVWGUplRwcHClx6vVatSrV6/Sc7RaLbQ3rWZWp04d5wN3gVq1avEfsxPYb85hvzmH/eY49plz2G/2u93ITSlJi21qNBrExMQgNTW13P7U1FS0adOm0nNat25d4fidO3ciNja20vk3RERE5HkkryaelJSEFStWYNWqVTh9+jTGjh2L9PR027o2kyZNwpAhQ2zHjxw5EhcuXEBSUhJOnz6NVatWYeXKlRg3bpxUb4GIiIjcjORzcAYOHIhr165h5syZyMjIQFRUFLZv346GDRsCADIyMpCenm47PiIiAtu3b8fYsWOxePFihIaG4v3333eLW8TtodVqMW3atAqXzOjfsd+cw35zDvvNcewz57DfxKMQ7LnXioiIiKgGkfwSFREREZGrMcEhIiIi2WGCQ0RERLLDBIeIiIhkhwlONbl06RKeffZZ1KtXD3q9Hvfffz8OHz4sdVhuLTw8HAqFosI2atQoqUNzayUlJXj99dcREREBb29vNGrUCDNnzoTFYpE6NLeXl5eHMWPGoGHDhvD29kabNm2QlpYmdVhuZd++fYiLi0NoaCgUCgW2bNlS7nFBEDB9+nSEhobC29sbHTp0wMmTJ6UJ1o3crt82bdqErl27wt/fHwqFAseOHZMkTjlhglMN/vrrL7Rt2xZeXl748ssvcerUKbz77ruSr6bs7tLS0pCRkWHbShd47N+/v8SRube3334by5Ytw6JFi3D69GnMnTsX77zzDhYuXCh1aG7vueeeQ2pqKj766CMcP34cXbp0wWOPPYZLly5JHZrbKCgoQIsWLbBo0aJKH587dy7mzZuHRYsWIS0tDcHBwejcuTPy8vKqOVL3crt+KygoQNu2bTFnzpxqjkzGBBLdhAkThIcffljqMGq8V155RWjcuLFgsVikDsWt9ejRQxg2bFi5fX369BGeffZZiSKqGQwGg6BSqYQvvvii3P4WLVoIU6ZMkSgq9wZA2Lx5s+1ni8UiBAcHC3PmzLHtKyoqEmrXri0sW7ZMggjd08399k/nzp0TAAhHjx6t1pjkiCM41eDzzz9HbGws+vfvj8DAQDzwwAP48MMPpQ6rRjGZTFi3bh2GDRtW7QVSa5qHH34Yu3fvxm+//QYA+Omnn7B//350795d4sjcW0lJCcxmM3Q6Xbn93t7e2L9/v0RR1Sznzp1DZmYmunTpYtun1WrRvn17HDhwQMLIyBMxwakGf/zxB5YuXYomTZpgx44dGDlyJF5++WWsXbtW6tBqjC1btuDvv/9GfHy81KG4vQkTJuDpp59G06ZN4eXlhQceeABjxozB008/LXVobs3Pzw+tW7fGm2++icuXL8NsNmPdunX44YcfkJGRIXV4NUJpIeSbiyUHBQVVKJJMJDbJSzV4AovFgtjYWMyePRsA8MADD+DkyZNYunRpuTpbdGsrV65Et27dEBoaKnUobi8lJQXr1q3Dxx9/jMjISBw7dgxjxoxBaGgohg4dKnV4bu2jjz7CsGHDcOedd0KlUiE6OhqDBg3CkSNHpA6tRrl5lFUQBI68UrXjCE41CAkJQfPmzcvta9asWbkaW3RrFy5cwK5du/Dcc89JHUqNMH78eEycOBFPPfUU7r33XgwePBhjx45FcnKy1KG5vcaNG2Pv3r3Iz8/HxYsX8eOPP6K4uBgRERFSh1YjBAcHA0CF0ZqsrKwKozpEYmOCUw3atm2LX3/9tdy+3377zVZQlP7d6tWrERgYiB49ekgdSo1gMBigVJb/X1ulUvE2cQf4+PggJCQEf/31F3bs2IEnn3xS6pBqhIiICAQHB9vueASs8+f27t2LNm3aSBgZeSJeoqoGY8eORZs2bTB79mwMGDAAP/74I5YvX47ly5dLHZrbs1gsWL16NYYOHQq1mv9c7REXF4e33noLYWFhiIyMxNGjRzFv3jwMGzZM6tDc3o4dOyAIAu655x6cOXMG48ePxz333IOEhASpQ3Mb+fn5OHPmjO3nc+fO4dixY7jjjjsQFhaGMWPGYPbs2WjSpAmaNGmC2bNnQ6/XY9CgQRJGLb3b9dv169eRnp6Oy5cvA4Dtj+Lg4GDbyBg5SOrbuDzF1q1bhaioKEGr1QpNmzYVli9fLnVINcKOHTsEAMKvv/4qdSg1Rm5urvDKK68IYWFhgk6nExo1aiRMmTJFMBqNUofm9lJSUoRGjRoJGo1GCA4OFkaNGiX8/fffUoflVvbs2SMAqLANHTpUEATrreLTpk0TgoODBa1WKzzyyCPC8ePHpQ3aDdyu31avXl3p49OmTZM07ppMIQiCUP1pFREREZF4OAeHiIiIZIcJDhEREckOExwiIiKSHSY4REREJDtMcIiIiEh2mOAQERGR7DDBISIiItlhgkNERESywwSHiIiIZIcJDpED1qxZA4VCAYVCgW+++abC44Ig4K677oJCoUCHDh0qnHf+/HmnX9OZc93dzz//jISEBERERECn08HX1xfR0dGYO3curl+/LnV4Ns7+3t1FfHy8Lf6oqCjbfine15YtW2yvqVAocOjQIZc8L9HNmOAQOcHPzw8rV66ssH/v3r04e/Ys/Pz8yu3v0aMHDh48iJCQEIdfqyrnurMPP/wQMTExSEtLw/jx4/HVV19h8+bN6N+/P5YtW4bhw4dLHWIFjv7e3UlwcDAOHjyIjz/+uMJj1fm+2rdvj4MHD+L111932XMSVYblmYmcMHDgQKxfvx6LFy9GrVq1bPtXrlyJ1q1bIzc3t9zxAQEBCAgIcOq1qnKuuzp48CBefPFFdO7cGVu2bIFWq7U91rlzZ7z66qv46quvXPJaBoMBer3eJc/l6O/dnWi1WrRq1arSx6rzfdWtWxetWrXCL7/84rLnJKoMR3CInPD0008DADZs2GDbl5OTg40bN2LYsGEVjq/sMtP06dOhUChw8uRJPP3006hduzaCgoIwbNgw5OTk2HXuzz//jP79+6N27dq44447kJSUhJKSEvz66694/PHH4efnh/DwcMydO7dcPPHx8QgPD68QZ+nzVrbP2deqzOzZs6FQKLB8+fJyyU0pjUaDJ554oty+/fv3o1OnTvDz84Ner0ebNm2wbdu2SmM9cuQI+vXrh7p166Jx48a2x3///XcMGjQIgYGB0Gq1aNasGRYvXnzbeEs5+nsHgDNnziAhIQFNmjSBXq/HnXfeibi4OBw/frzccVevXsXzzz+PBg0aQKvVIiAgAG3btsWuXbscOsYZzrwvInfHBIfICbVq1UK/fv2watUq274NGzZAqVRi4MCBDj1X3759cffdd2Pjxo2YOHEiPv74Y4wdO9aucwcMGIAWLVpg48aNGDFiBN577z2MHTsWvXr1Qo8ePbB582Y8+uijmDBhAjZt2uRQXGK9ltlsxtdff42YmBg0aNDArtfeu3cvHn30UeTk5GDlypXYsGED/Pz8EBcXh5SUlArH9+nTB3fddRc++eQTLFu2DABw6tQpPPjggzhx4gTeffddfPHFF+jRowdefvllzJgxw644nPm9X758GfXq1cOcOXPw1VdfYfHixVCr1XjooYfw66+/2o4bPHgwtmzZgqlTp2Lnzp1YsWIFHnvsMVy7ds2hY5zhyn/PRG5DICK7rV69WgAgpKWlCXv27BEACCdOnBAEQRAefPBBIT4+XhAEQYiMjBTat29f4bxz587Z9k2bNk0AIMydO7fcayQmJgo6nU6wWCy3Pffdd98td+79998vABA2bdpk21dcXCwEBAQIffr0se0bOnSo0LBhwwrvr/R5K9vn7GvdLDMzUwAgPPXUU7c85matWrUSAgMDhby8PNu+kpISISoqSqhfv76tr0pjnTp1aoXn6Nq1q1C/fn0hJyen3P7Ro0cLOp1OuH79+i1f39nfe2VKSkoEk8kkNGnSRBg7dqxtv6+vrzBmzJh/PdeeYypzq9+3K9+Xo/752kRi4AgOkZPat2+Pxo0bY9WqVTh+/DjS0tKcGs6/+VLMfffdh6KiImRlZd323J49e5b7uVmzZlAoFOjWrZttn1qtxl133YULFy44HJtUr/VPBQUF+OGHH9CvXz/4+vra9qtUKgwePBh//vlnuZEQwDoq9k9FRUXYvXs3evfuDb1ej5KSEtvWvXt3FBUV4fvvv7crHkd/7yUlJZg9ezaaN28OjUYDtVoNjUaD33//HadPn7Yd17JlS6xZswazZs3C999/j+Li4grPZc8xzqrKv+fMzEy8/PLLiI6ORvv27bFkyRJYLBaXxUbkDCY4RE5SKBRISEjAunXrsGzZMtx9991o166dw89Tr169cj+XzkkpLCy87bl33HFHuZ81Gg30ej10Ol2F/UVFRQ7HJsZr+fv7Q6/X49y5c3a97l9//QVBECq9iyw0NBQAKlyiufnYa9euoaSkBAsXLoSXl1e5rXv37gCA7Oxsu+Jx9PeelJSEN954A7169cLWrVvxww8/IC0tDS1atCj3O05JScHQoUOxYsUKtG7dGnfccQeGDBmCzMxMh45xlrP/ni9duoRHH30U999/PyZMmIBJkybh8OHDGDJkSJVjIqoKJjhEVRAfH4/s7GwsW7YMCQkJUodjN51OB6PRWGG/vV/yVaFSqdCpUyccPnwYf/75522Pr1u3LpRKJTIyMio8dvnyZQDWpOmfbp4oXbduXahUKsTHxyMtLa3SrTTRsYcjv/d169ZhyJAhmD17Nrp27YqWLVsiNja2Ql/7+/tj/vz5OH/+PC5cuIDk5GRs2rQJ8fHxDh1TFc78e54wYQImT56MYcOG4csvv8Tff/+NlStX4tq1a9ixY4dL4iJyBhMcoiq48847MX78eMTFxWHo0KFSh2O38PBwZGVl4cqVK7Z9JpOp2r6QJk2aBEEQMGLECJhMpgqPFxcXY+vWrQAAHx8fPPTQQ9i0aVO5EQ+LxYJ169ahfv36uPvuu//19fR6PTp27IijR4/ivvvuQ2xsbIXt5pG0f+PI712hUFS4U2zbtm24dOnSLc8JCwvD6NGj0blzZxw5csTpYxzlzL/nb775xnYX1j8NHToUX375pUviInIG18EhqqI5c+ZIHYLDBg4ciKlTp+Kpp57C+PHjUVRUhPfffx9ms7laXr9169ZYunQpEhMTERMTgxdffBGRkZEoLi7G0aNHsXz5ckRFRSEuLg4AkJycjM6dO6Njx44YN24cNBoNlixZghMnTmDDhg0VRmwqs2DBAjz88MNo164dXnzxRYSHhyMvLw9nzpzB1q1b8fXXXzv0Huz9vffs2RNr1qxB06ZNcd999+Hw4cN45513UL9+fdsxOTk56NixIwYNGoSmTZvCz88PaWlp+Oqrr9CnTx+7j3EFZ/49Dx8+HP/73/9QXFyMjz/+GPHx8ejRoweCg4NdFheRo5jgEHmgiIgIfPbZZ5g8eTL69euHkJAQJCUl4erVq3bfMl1VI0aMQMuWLfHee+/h7bffRmZmJry8vHD33Xdj0KBBGD16tO3Y9u3b4+uvv8a0adMQHx8Pi8WCFi1a4PPPP68w+flWmjdvjiNHjuDNN9/E66+/jqysLNSpUwdNmjRx6PKUoxYsWAAvLy8kJycjPz8f0dHR2LRpU7mVfHU6HR566CF89NFHOH/+PIqLixEWFoYJEybgtddes/sYKTzwwAMYNGgQ1qxZg/j4eDz++ON46qmnMHToULRv316yuIgUgiAIUgdBRETiiY+PxzfffIMzZ85AoVBApVK57LlPnjyJvn374sMPP8R3332HBx54AEeOHMHOnTuxa9euCq8lCALMZjPWrl2L4cOHIy0tDbGxsS6Lh6gUR3CIiDzAhQsX4OXlhcjISJw4ccJlzxsZGYn/+7//w6uvvorTp09DqVSid+/e2LZtW6WJ1GeffYbevXu77PWJboUjOEREMnf+/HnbXVve3t6IjIyULJa///4bZ86csf3cvHlzl9UKI/onJjhEREQkO7xNnIiIiGSHCQ4RERHJDhMcIiIikh0mOERERCQ7THCIiIhIdpjgEBERkewwwSEiIiLZYYJDREREssMEh4iIiGSHCQ4RERHJzv8D5UHsXvoT6xkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.hist(Mout_error,bins=30,density=True,histtype='step',lw=2) \n", "#note if density=True this simply means are under the histogram is 1, not the height of the bars sums to 1!\n", "plt.plot([8.08,8.08],[0.,1.],':',color='b', lw=1.5) #Must set this to the meassured mass of the planet\n", "plt.plot([(np.percentile(Mout_error,5.)),(np.percentile(Mout_error,5.))],[0.,1.],'--',color='k', lw=1.5)\n", "plt.ylim((0.,1.))\n", "plt.ylabel('Probability density',fontsize=12)\n", "plt.xlabel(r'Minimum Core Mass [M$_\\oplus$]',fontsize=12)\n", "print('The 95% upper limit to be consistent with photoevaporation is')\n", "print((np.percentile(Mout_error,5.)),'Mearth')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we will run for the two other efficiency options to show the effect on the result" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "eff_option = 1\n", "\n", "Mout_error_1= np.zeros(1000)\n", "flag_out_1= np.zeros(1000,dtype=np.int8)\n", "\n", "for i in range(np.size(Mout_error)):\n", " \n", " K36b_R_use = normal(K36b_R,K36b_Rer,1)\n", " K36b_P_use = normal(K36b_P,K36b_Per,1)\n", " \n", " K36c_R_use = normal(K36c_R,K36c_Rer,1)\n", " K36c_P_use = normal(K36c_P,K36c_Per,1)\n", " \n", " system = em.psystem('Kepler36_%d' %i)\n", " system.add_planet('36b',K36b_R_use,K36b_Rer,K36b_P_use,K36b_Per)\n", " system.add_planet('36c',K36c_R_use,K36c_Rer,K36c_P_use,K36c_Per)\n", " system.star.mass=normal(K36s_M,K36s_Mer,1)\n", " system.star.radius=normal(K36s_R,K36s_Rer,1)\n", " system.star.Teff=normal(K36s_Teff,K36s_Teffer,1)\n", " system.star.age = normal(K36s_age,K36s_age_er,1)\n", " system.update_planet_info()\n", " system.above_or_below_valley()\n", " system.mass_rocky(Xiron,Xice)\n", " \n", " ms.find_hardest_rocky(system,Tkh_Myr,Xiron,Xice,eff_option)\n", "\n", " Mout, flag = ms.min_mass_gaseous(system.planets[0],system.planets[1],Tkh_Myr,Xiron,Xice,system.star.age,eff_option)\n", " \n", " Mout_error_1[i] = Mout\n", " flag_out_1[i]=flag" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "eff_option = 2\n", "\n", "Mout_error_2= np.zeros(1000)\n", "flag_out_2= np.zeros(1000,dtype=np.int8)\n", "\n", "for i in range(np.size(Mout_error)):\n", " \n", " K36b_R_use = normal(K36b_R,K36b_Rer,1)\n", " K36b_P_use = normal(K36b_P,K36b_Per,1)\n", " \n", " K36c_R_use = normal(K36c_R,K36c_Rer,1)\n", " K36c_P_use = normal(K36c_P,K36c_Per,1)\n", " \n", " system = em.psystem('Kepler36_%d' %i)\n", " system.add_planet('36b',K36b_R_use,K36b_Rer,K36b_P_use,K36b_Per)\n", " system.add_planet('36c',K36c_R_use,K36c_Rer,K36c_P_use,K36c_Per)\n", " system.star.mass=normal(K36s_M,K36s_Mer,1)\n", " system.star.radius=normal(K36s_R,K36s_Rer,1)\n", " system.star.Teff=normal(K36s_Teff,K36s_Teffer,1)\n", " system.star.age = normal(K36s_age,K36s_age_er,1)\n", " system.update_planet_info()\n", " system.above_or_below_valley()\n", " system.mass_rocky(Xiron,Xice)\n", " \n", " ms.find_hardest_rocky(system,Tkh_Myr,Xiron,Xice,eff_option)\n", "\n", " Mout, flag = ms.min_mass_gaseous(system.planets[0],system.planets[1],Tkh_Myr,Xiron,Xice,system.star.age,eff_option)\n", " \n", " Mout_error_2[i] = Mout\n", " flag_out_2[i]=flag" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAG6CAYAAAAF5Ty4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABqy0lEQVR4nO3dd1gUV9sG8HsBWUCKShFQQGyxYC8RGxpU7DE2jEbBFg1Go9gr9hqV2EsUXxM1xFhib1ij0RhLNIomKooFglgoShH2fH8Q9su6lJ1lYWG5f9e1V9wzZ2aeGcjuw5lTZEIIASIiIiIDYqTvAIiIiIh0jQkOERERGRwmOERERGRwmOAQERGRwWGCQ0RERAaHCQ4REREZHCY4REREZHCY4BAREZHBYYJDREREBocJDhERERkcvSc4Z8+eRZcuXeDs7AyZTIa9e/fmus+ZM2fQoEEDmJmZoWLFili3bl3+B0pERERFht4TnDdv3qBOnTpYtWqVRvUjIiLQsWNHtGjRAteuXcOUKVMwatQo7Nq1K58jJSIioqJCVpgW25TJZNizZw+6deuWbZ2JEydi3759CA8PV5YNHz4cf/zxB3799dcCiJKIiIgKOxN9ByDVr7/+inbt2qmU+fj4YNOmTXj37h1KlCihtk9KSgpSUlKU7xUKBV6+fAlbW1vIZLJ8j5mIiIjyTgiBhIQEODs7w8go54dQRS7BiY6ORtmyZVXKypYti7S0NMTGxsLJyUltnwULFmDWrFkFFSIRERHlo8ePH6N8+fI51ilyCQ4AtVaXzKds2bXGTJ48GYGBgcr3cXFxcHV1xePHj2FtbZ1/gRIREZHOxMfHw8XFBVZWVrnWLXIJjqOjI6Kjo1XKYmJiYGJiAltb2yz3kcvlkMvlauXW1tZMcIiIiIoYTbqX6H0UlVSenp44fvy4StmxY8fQsGHDLPvfEBERUfGj9wQnMTER169fx/Xr1wFkDAO/fv06IiMjAWQ8XhowYICy/vDhw/Ho0SMEBgYiPDwcmzdvxqZNmzBu3Dh9hE9ERESFkN4fUf3+++9o3bq18n1mXxk/Pz9s2bIFUVFRymQHANzd3XHo0CGMGTMGq1evhrOzM1asWIEePXoUeOxERERUOBWqeXAKSnx8PGxsbBAXF8c+OEQSKRQKpKam6jsMIjJQpqam2Q4Bl/L9rfcWHCIqOlJTUxEREQGFQqHvUIjIQBkZGcHd3R2mpqZ5Og4THCLSiBACUVFRMDY2houLS66TbBERSaVQKPDs2TNERUXB1dU1T5PxMsEhIo2kpaXh7du3cHZ2hoWFhb7DISIDZW9vj2fPniEtLS1Po6P5JxgRaSQ9PR0A8txsTESUk8zPmMzPHG0xwSEiSbh+GxHlJ119xjDBISIiIoPDPjhEpLUuK3/B84SUAj+vvZUc+0c2L/DzFjczZ87E3r17lROxFhapqamoUaMG/ve//6FZs2b6DkfNuHHjkJqaihUrVuT5WHv37sW4ceMQERGBkSNHIjg4WK2sbt26GD16NF6/fq3RMStUqIDRo0dj9OjReY6vMGMLDhFp7XlCCqLjkwv8JTWpevz4MQYPHgxnZ2eYmprCzc0NX331FV68eJFPdyZvjhw5gnr16sHc3BzlypVDQECARvv5+/ujW7du+RtcIbBhwwa4ubkpk5uHDx9i8ODBcHd3h7m5OSpVqoSgoCC1+ZoiIyPRpUsXlCxZEnZ2dhg1apRKneTkZPj7+6NWrVowMTHJ8l7u3r0bbdu2hb29PaytreHp6YmjR4+q1JkwYQJCQkIQERGR52sdNmwYevbsicePH2POnDlZlvn6+uKvv/7S+JiXL1/G559/nufYCju24BBRnhnJAAcrs3w/T0xCMhQSpyZ98OABPD09UbVqVezYsQPu7u64desWxo8fj8OHD+PixYsoU6ZM/gSsheTkZHTv3h29evXCzp07kZiYiN9++03fYRUqK1euxMyZM5Xv79y5A4VCgfXr16Ny5cr4888/MXToULx58wZff/01gIwOq506dYK9vT1++eUXvHjxAn5+fhBCYOXKlco65ubmGDVqFHbt2pXluc+ePYu2bdti/vz5KFWqFEJCQtClSxdcunQJ9erVAwA4ODigXbt2WLduHRYtWqT1dSYmJiImJgY+Pj5wdnbOtgwAzM3NNT6uvb291jEVKaIYiouLEwBEXFycvkMhKjKSkpLE7du3RVJSkrLsw3knhNvEA+LDeScKJAZtzte+fXtRvnx58fbtW5XyqKgoYWFhIYYPHy6EEGLFihXCw8NDuX3Pnj0CgFi1apWyrF27dmLSpEnK9/v27RP169cXcrlcuLu7i5kzZ4p3794ptwMQGzduFN26dRPm5uaicuXK4ueff84x3qSkJGFpaSmOHj2q8TVm8vPzEx9//LHy/eHDh0WzZs2EjY2NKFOmjOjUqZO4d++eyj6PHz8Wvr6+onTp0sLCwkI0aNBAXLx4UQghRFBQkKhTp46y7oMHD0SlSpXE8OHDRXp6unj48KHo3LmzKFWqlLCwsBA1atQQBw8eVNY/ffq0aNSokTA1NRWOjo5i4sSJKvfHy8tLjBw5UowfP16ULl1alC1bVgQFBeV4jVeuXBFGRka5fn4vXrxYuLu7K98fOnRIGBkZiadPnyrLduzYIeRyeZbHev9e5qRGjRpi1qxZKmVbtmwRLi4uOe6XkpIixo8fL5ydnYWFhYVo3LixOHXqlBBCiFOnTgkAKq/sykJCQoSNjY3KsX/++WfRoEEDIZfLha2trfjkk0+U29zc3MTy5cuV71+/fi2GDh0q7O3thZWVlWjdurW4fv26cnvm78HWrVuFm5ubsLa2Fr6+viI+Pl5ZJz09XSxcuFBUqlRJmJqaChcXFzF37lwhhBCtW7cWI0aMUIkvNjZWmJqairCwMLX7ktVnTSYp3998REVEBuvly5c4evQoAgIC1P7CdXR0RL9+/RAaGgohBFq1aoVbt24hNjYWAHDmzBnY2dnhzJkzADLmAbpw4QK8vLwAAEePHsVnn32GUaNG4fbt21i/fj22bNmCefPmqZxn1qxZ6N27N27cuIGOHTuiX79+ePnyZbYxm5mZwcfHBxMmTMixnibevHmDwMBAXL58GWFhYTAyMsInn3yinIk6MTERXl5eePbsGfbt24c//vgDEyZMyHKm6j///BPNmjVDr169sHbtWhgZGWHEiBFISUnB2bNncfPmTSxatAiWlpYAgKdPn6Jjx45o1KgR/vjjD6xduxabNm3C3LlzVY77v//9DyVLlsSlS5ewePFizJ49G8ePH8/2ms6ePYuqVavmOk1/XFycSsvcr7/+Cg8PD5VWDx8fH6SkpODKlSu538xsKBQKJCQkqLUCNm7cGI8fP8ajR4+y3XfgwIE4f/48fvjhB9y4cQO9evVC+/bt8ffff6Np06a4e/cuAGDXrl2IiorKtux9Bw8eRPfu3dGpUydcu3YNYWFhaNiwYZYxCCHQqVMnREdH49ChQ7hy5Qrq168Pb29vld+/+/fvY+/evThw4AAOHDiAM2fOYOHChcrtkydPxqJFizB9+nTcvn0b27dvR9myZQEAQ4YMwfbt25GS8v+Plrdt2wZnZ2eVtSh1LtcUyACxBYdIuqLYgnPx4kUBQOzZsyfL7cuWLRMAxD///CMUCoWws7MTP/30kxBCiLp164oFCxYIBwcHIYQQFy5cECYmJiIhIUEIIUSLFi3E/PnzVY733XffCScnJ+V7AGLatGnK94mJiUImk4nDhw9nG/PMmTNFxYoVxdSpU4WHh4dKi8OIESNE586ds903t1aHmJgYAUDcvHlTCCHE+vXrhZWVlXjx4kWW9TP/cr9w4YIoU6aMWLJkicr2WrVqiZkzZ2a575QpU8QHH3wgFAqFsmz16tXC0tJSpKenCyEyWnCaN2+usl+jRo3ExIkTs72Gr776Snz00UfZbhdCiHv37glra2uxceNGZdnQoUNF27Zt1eqampqK7du3q5Vr2oKzePFiUaZMGfHPP/+olGd+z5w+fTrbGGUymcrPVwghvL29xeTJk4UQQrx69UrZSpMpq7L3W3A8PT1Fv379so35vy04YWFhwtraWiQnJ6vUqVSpkli/fr0QIuP3wMLCQqXFZvz48eLDDz8UQggRHx8v5HK5yv3+r+TkZFGmTBkRGhqqLKtbt262vztswSEiyiPx71rDMpkMMpkMLVu2xOnTp/H69WvcunULw4cPR3p6OsLDw3H69GnUr19f2UJx5coVzJ49G5aWlsrX0KFDERUVhbdv3yrPUbt2beW/S5YsCSsrK8TExGQZz6tXr7BgwQKsXLkSc+fOxSeffIJmzZrh77//BpDRitK8ueajx+7fv4++ffuiYsWKsLa2hru7O4CMzrYAcP36ddSrVy/HPkiRkZFo06YNpk2bhnHjxqlsGzVqFObOnYtmzZohKCgIN27cUG4LDw+Hp6enypwmzZo1Q2JiIp48eaIs++/9AQAnJ6ds7w8AJCUlwcws+/5ez549Q/v27dGrVy8MGTJEZVtW86sIIbSed2XHjh2YOXMmQkND4eDgoLIts8Xwv78L/3X16lUIIVC1alWV36EzZ87g/v37WsWT6fr16/D29tao7pUrV5CYmAhbW1uVOCIiIlTiqFChAqysrJTv//tzCg8PR0pKSrbnlMvl+Oyzz7B582ZlfH/88Qf8/f21vELNsJMxERmsypUrQyaT4fbt21mOiLlz5w5Kly4NOzs7AECrVq2wYcMGnDt3DnXq1EGpUqXQsmVLnDlzBqdPn0arVq2U+yoUCsyaNQvdu3dXO+5/v4Dfn2peJpNlu1jp3bt3kZKSouysOnv2bMTHx6N58+YIDg7GxYsXsW3bNo2vv0uXLnBxccHGjRvh7OwMhUIBDw8P5cghTTqm2tvbw9nZGT/88AMGDx6s8mhoyJAh8PHxwcGDB3Hs2DEsWLAAS5cuxciRI7NMHP6bUGaScn8AwM7ODjdv3sxy27Nnz9C6dWt4enpiw4YNKtscHR1x6dIllbJXr17h3bt3ykcpUoSGhmLw4MHYuXMn2rRpo7Y98/FOdh16FQoFjI2NceXKFRgbG6tsy0yitSWlw7FCoYCTkxNOnz6ttq1UqVLKf+f0c9LkfEOGDEHdunXx5MkTbN68Gd7e3nBzc9M4Tm2wBYeIDJatrS3atm2LNWvWICkpSWVbdHQ0tm3bBl9fX+UXbmY/nJ9++kmZzHh5eeHEiRMq/W8AoH79+rh79y4qV66s9tJ2IdJy5coByOhnkmn58uXo0qUL+vbti2HDhinr5ObFixcIDw/HtGnT4O3tjerVq+PVq1cqdWrXro3r16/n2NfH3NwcBw4cUPYNSkhIUNnu4uKC4cOHY/fu3Rg7diw2btwIAKhRowYuXLigTGoA4MKFC7CystL4GrJSr1493LlzR+W4QEafn1atWqF+/foICQlR+xl4enrizz//RFRUlLLs2LFjkMvlaNCggaQYduzYAX9/f2zfvh2dOnXKss6ff/6JEiVKoGbNmtleR3p6OmJiYtR+fxwdHSXF877atWsjLCxMo7r169dHdHQ0TExM1OLITPxzU6VKFZibm+d4zlq1aqFhw4bYuHEjtm/fjkGDBml07LxggkNEBm3VqlVISUmBj48Pzp49i8ePH+PIkSNo27YtypUrp9Ip2MPDA7a2tti2bZsywWnVqhX27t2LpKQklcdDM2bMwNatWzFz5kzcunUL4eHhCA0NxbRp07SO1cXFBX369MGIESOwdetW3L9/HydOnMCNGzdQsmRJ7Nu3L8fHN/9VunRp2NraYsOGDbh37x5OnjyJwMBAlTqffvopHB0d0a1bN5w/fx4PHjzArl278Ouvv6rUK1myJA4ePAgTExN06NABiYmJAIDRo0fj6NGjiIiIwNWrV3Hy5ElUr14dABAQEIDHjx9j5MiRuHPnDn7++WcEBQUhMDAwTyvRt27dGm/evMGtW7eUZc+ePUOrVq3g4uKCr7/+Gs+fP0d0dDSio6OVddq1a4caNWqgf//+yo6348aNw9ChQ1VapW7fvq1M+uLi4nD9+nWViQ537NiBAQMGYOnSpWjSpInyPHFxcSpxnjt3Di1atMi2daNq1aro168fBgwYgN27dyMiIgKXL1/GokWLcOjQIa3vDwAEBQVhx44dCAoKQnh4OG7evInFixdnWbdNmzbw9PREt27dcPToUTx8+BAXLlzAtGnT8Pvvv2t0PjMzM0ycOBETJkxQ/t5evHgRmzZtUqk3ZMgQLFy4EOnp6fjkk0/ydI2a4CMqIsqzmIRkNJmv2V+MeT2PVFWqVMHvv/+OmTNnwtfXFy9evFB+qQcFBan0P5HJZPDy8sLevXvRokULABl/DdvY2Cj7sWTy8fHBgQMHMHv2bCxevBglSpRAtWrV1Pp9SPW///0PS5Yswbx58/Do0SOUK1cOn332GQ4fPgxvb2907doVp06dyvKLU6FQwMQk42PdyMgIP/zwA0aNGgUPDw988MEHWLFihcpjNlNTUxw7dgxjx45Fx44dkZaWhho1amD16tVqx7a0tMThw4fh4+ODjh074vDhw0hPT8eIESPw5MkTWFtbo3379li+fDmAjNaoQ4cOYfz48ahTpw7KlCmDwYMH5ykBBDJa5bp3745t27ZhwYIFADJaYu7du4d79+6hfPnyKvUzW3qMjY1x8OBBBAQEoFmzZjA3N0ffvn2V8+Rk6tixo8rIp8zHhZnHWb9+PdLS0jBixAiMGDFCWc/Pzw9btmxRvt+xYwdmzZqV47WEhIRg7ty5GDt2LJ4+fQpbW1t4enqiY8eOEu+KqlatWmHnzp2YM2cOFi5cCGtra7Rs2TLLujKZDIcOHcLUqVMxaNAgPH/+HI6OjmjZsqWkR3fTp0+HiYkJZsyYgWfPnsHJyQnDhw9XqfPpp59i9OjR6Nu3b479qHRFJt5v5ysG4uPjYWNjg7i4uFyHGhJRhuTkZERERMDd3V354dRkfhii46UnHXnlaG2Gi1M060RZnLRv3x6VK1fGqlWr9B1Kvrp58ybatGmDe/fuqXR8LSwOHjyI8ePH48aNG8qEkzJmFK9QoQIuX76M+vXrZ1svq8+aTFK+v3nniUhr9lbyYnXewurVq1e4cOECTp8+rfZXsyGqVasWFi9ejIcPH6JWrVr6DkfNmzdvEBISwuTmX+/evUNUVBQmTZqEJk2a5Jjc6BLvPhFpjQteFg6DBg3C5cuXMXbsWHz88cf6DqdA+Pn56TuEbPXu3VvfIRQq58+fR+vWrVG1alX89NNPBXZeJjhEREXcnj179B0CUbZatWqlNuqtIHAUFRERERkcJjhERERkcJjgEBERkcFhgkNEREQGhwkOERERGRwmOERERGRwOEyciLS33gtI1GxtJJ2ydACGnSn48xYDd+7cgb+/P65fv45q1arh+vXramV79+6Fu7s7rl27hrp16+Z6TH9/f7x+/Rp79+7N9/iJMjHBISLtJcYACc/0HUWuoqOjMW/ePBw8eBBPnz6Fg4MD6tati9GjR8Pbu2CWfMjPL/lWrVqhbt26CA4OzrXemTPqieGwYcOwbt06ABkLNZYsWRJ3796FpaVllmWlSpVCVFSUxqtNf/PNN3qZB4WKNyY4RJR3MiPA0jH/z5MYDQiFpF0ePnyIZs2aoVSpUli8eDFq166Nd+/e4ejRoxgxYgTu3LmTT8EWTkOHDsXs2bNVyiwsLJT/vn//Pjp16gQ3N7ccyxwdNf9529jY5CFiIi2JYiguLk4AEHFxcfoOhajISEpKErdv3xZJSUn/X/h1NSGCrDP+WxC0OF+HDh1EuXLlRGJiotq2V69eKf/96NEj0bVrV1GyZElhZWUlevXqJaKjo5Xbg4KCRJ06dcTWrVuFm5ubsLa2Fr6+viI+Pl5ZZ+fOncLDw0OYmZmJMmXKCG9vb5GYmCiCgoIEAJXXqVOnhBBCTJgwQVSpUkWYm5sLd3d3MW3aNJGamqrxef38/NSOHRERkeW98PLyEl999VW29+r942QVd1BQkIiIiBAAxLVr15T7/vnnn6Jjx47CyspKWFpaiubNm4t79+4pY/z444+VdRUKhVi0aJFwd3cXZmZmonbt2mLnzp3K7adOnRIAxIkTJ0SDBg2Eubm58PT0FHfu3FGJ9+effxYNGjQQcrlc2Nraik8++UQIIcSsWbOEh4eH2vXVr19fTJ8+Pdvrp8Ihy8+af0n5/mYnYyIyWC9fvsSRI0cwYsQIlCxZUm17qVKlAABCCHTr1g0vX77EmTNncPz4cdy/fx++vr4q9e/fv4+9e/fiwIEDOHDgAM6cOYOFCxcCAKKiovDpp59i0KBBCA8Px+nTp9G9e3cIITBu3Dj07t0b7du3R1RUFKKiotC0aVMAgJWVFbZs2YLbt2/jm2++wcaNG7F8+XKNz/vNN9/A09MTQ4cOVR7bxcVFq/sVFRWFmjVrYuzYsYiKisK4ceOyLHvf06dP0bJlS5iZmeHkyZO4cuUKBg0ahLS0tCzPM23aNISEhGDt2rW4desWxowZg88++0zt8dnUqVOxdOlS/P777zAxMcGgQYOU2w4ePIju3bujU6dOuHbtGsLCwtCwYUMAGWtz3b59G5cvX1bWv3HjBq5duwZ/f3+t7g0VPXxERUQG6969exBCoFq1ajnWO3HiBG7cuIGIiAhlcvDdd9+hZs2auHz5Mho1agQAUCgU2LJlC6ysrAAA/fv3R1hYGObNm4eoqCikpaWhe/fuykc5/13p2tzcHCkpKWqPdqZNm6b8d4UKFTB27FiEhoZiwoQJyvKczmtjYwNTU1NYWFho9NhozZo1+Pbbb1XKVq9eDT8/Pzg6OsLExASWlpbKY1laWqqVxcbGqu1vY2ODH374ASVKlAAAVK1aNcvzv3nzBsuWLcPJkyfh6ekJAKhYsSJ++eUXrF+/Hl5eXsq68+bNU76fNGkSOnXqhOTkZJiZmWHevHno06cPZs2apaxfp04dAED58uXh4+ODkJAQ5c8uJCQEXl5eqFixYq73iAwDExwiMlji346tMpksx3rh4eFwcXFRafmoUaMGSpUqhfDwcOWXZIUKFZRJBgA4OTkhJiZjFFmdOnXg7e2NWrVqwcfHB+3atUPPnj1RunTpHM/9008/ITg4GPfu3UNiYiLS0tJgbW2tUien80rVr18/TJ06VaXMwcFBq2Nlun79Olq0aKFMbnJy+/ZtJCcno23btirlqampqFevnkpZ7dq1lf92cnICAMTExMDV1RXXr1/H0KFDsz3P0KFDMWjQICxbtgzGxsbYtm0bli5dKuWyqIhjgkNEBqtKlSqQyWQIDw9Ht27dsq0nhMgyCXq//P0vcJlMBoUio9OzsbExjh8/jgsXLuDYsWNYuXIlpk6dikuXLsHd3T3L8168eFHZCuHj46NsBXn/izin80plY2ODypUra7VvdszNzTWumxn3wYMHUa5cOZVtcrlc5f1/rzvz55C5f27n7NKlC+RyOfbs2QO5XI6UlBT06NFD4zip6GMfHCIyWGXKlIGPjw9Wr16NN2/eqG1//fo1gIzWmsjISDx+/Fi57fbt24iLi0P16tU1Pp9MJkOzZs0wa9YsXLt2DaamptizZw8AwNTUFOnp6Sr1z58/Dzc3N0ydOhUNGzZElSpV8OjRI8nXmdWxC1Lt2rVx7tw5vHv3Lte6NWrUgFwuR2RkJCpXrqzyktJ3qHbt2ggLC8t2u4mJCfz8/BASEoKQkBD06dNHZbQYGT624BCRQVuzZg2aNm2Kxo0bY/bs2ahduzbS0tJw/PhxrF27FuHh4WjTpg1q166Nfv36ITg4GGlpaQgICICXl5ey42puLl26hLCwMLRr1w4ODg64dOkSnj9/rkyQKlSogKNHj+Lu3buwtbVVtqRERkbihx9+QKNGjXDw4EFlQiRFhQoVcOnSJTx8+BCWlpYoU6YMjIyy/vv17du3iI6OVimTy+W5PkrLyZdffomVK1eiT58+mDx5MmxsbHDx4kU0btwYH3zwgUpdKysrjBs3DmPGjIFCoUDz5s0RHx+PCxcuwNLSEn5+fhqdMygoCN7e3qhUqRL69OmDtLQ0HD58WKXv0pAhQ5T3//z581pfHxVNTHCIKO8So4Glmrd05Ok8Erm7u+Pq1auYN2+eciSQvb09GjRogLVr1wLIaHnZu3cvRo4ciZYtW8LIyAjt27fHypUrNT6PtbU1zp49i+DgYMTHx8PNzQ1Lly5Fhw4dAGT0CTl9+jQaNmyIxMREnDp1Ch9//DHGjBmDL7/8EikpKejUqROmT5+OmTNnSrrGcePGwc/PDzVq1EBSUhIiIiJQoUKFLOtu3LgRGzduVCnz8fHBkSNHJJ3zv2xtbXHy5EmMHz8eXl5eMDY2Rt26ddGsWbMs68+ZMwcODg5YsGABHjx4gFKlSqF+/fqYMmWKxuds1aoVdu7ciTlz5mDhwoWwtrZGy5YtVepUqVIFTZs2xYsXL/Dhhx9qfX1UNMmEKH7TS8bHx8PGxgZxcXFqnfmIKGvJycmIiIiAu7s7zMzMMgqXVtfPTMZWzsDY8II/LxUpmSPohg0bhsDAQH2HQxrK8rPmX1K+v9mCQ0Tas8zb6Jsid14qMmJiYvDdd9/h6dOnGDhwoL7DIT1ggkNE2uOCl1RIlS1bFnZ2dtiwYUOe+hdR0cUEh4iIDE4x7H1B7+EwcSIiIjI4THCIiIjI4DDBISIiIoPDBIeIiIgMDhMcIiIiMjhMcIiIiMjgcJg4EWnN94AvYpNiC/y8duZ2CO0cWuDnNSQPHz6Eu7s7rl27hrp16xa64xWUVq1aoW7duggODgaQsa7X6NGjMXr0aL3GRXnHBIeItBabFIuYtzH6DiNH/v7++N///qdWntf1l0iVi4sLoqKiYGdnBwA4ffo0WrdujVevXqFUqVL6DU6Cy5cvo2TJkvoOg3SACQ4R5ZmRzAh25nb5fp7YpFgohELyfu3bt0dISIhKmVwu11VYBMDY2BiOjo76DiPP7O3t9R0C6YoohuLi4gQAERcXp+9QiIqMpKQkcfv2bZGUlKQs++jHj4THFg/x0Y8fFUgM2pzPz89PfPzxxznWCQoKEi4uLsLU1FQ4OTmJkSNHKrclJyeL8ePHi/LlywtTU1NRuXJl8e233wohhEhLSxODBg0SFSpUEGZmZqJq1aoiODg4y/PPnDlT2NvbCysrK/H555+LlJQUZR2FQiEWLVok3N3dhZmZmahdu7bYuXNntvFOmjRJfPjhh2rltWrVEjNmzFC+37x5s6hWrZqQy+Xigw8+EKtXr1Zui4iIEADEtWvXlGWnT58WjRo1EqampsLR0VFMnDhRvHv3Trk9PT1dLFy4UFSqVEmYmpoKFxcXMXfuXLXjZf77vy8/Pz/xv//9T5QpU0YkJyerxN29e3fRv3//LK81JSVFjBgxQjg6Ogq5XC7c3NzE/PnzldtfvXolhg4dKhwcHIRcLhc1a9YU+/fvF0IIERsbK/r06SPKlSsnzM3NhYeHh9i+fbvK8b28vMRXX32lfO/m5iaWL1+ufA9AbNy4UXTr1k2Ym5uLypUri59//lnlGD///LOoXLmyMDMzE61atRJbtmwRAMSrV6+yvCbKWVafNZmkfH+zBYeIirWffvoJy5cvxw8//ICaNWsiOjoaf/zxh3L7gAED8Ouvv2LFihWoU6cOIiIiEBub0e9IoVCgfPny+PHHH2FnZ4cLFy7g888/h5OTE3r37q08RlhYGMzMzHDq1Ck8fPgQAwcOhJ2dHebNmwcAmDZtGnbv3o21a9eiSpUqOHv2LD777DPY29vDy8tLLeZ+/fph4cKFuH//PipVqgQAuHXrFm7evImffvoJALBx40YEBQVh1apVqFevHq5du4ahQ4eiZMmS8PPzUzvm06dP0bFjR/j7+2Pr1q24c+cOhg4dCjMzM8ycORMAMHnyZGzcuBHLly9H8+bNERUVhTt37qgdy8XFBbt27UKPHj1w9+5dWFtbw9zcHKamphg1ahT27duHXr16AQBiY2Nx4MCBbB8XrlixAvv27cOPP/4IV1dXPH78GI8fP1be/w4dOiAhIQHff/89KlWqhNu3b8PY2BhAxqrUDRo0wMSJE2FtbY2DBw+if//+qFixIj788MMcfitUzZo1C4sXL8aSJUuwcuVK9OvXD48ePUKZMmXw8OFD9OzZE1999RWGDBmCa9euYdy4cRofm/JRfmRfhR1bcIikK8otOMbGxqJkyZIqr9mzZwshhFi6dKmoWrWqSE1NVdv37t27AoA4fvy4xucLCAgQPXr0UDl/mTJlxJs3b5Rla9euFZaWliI9PV0kJiYKMzMzceHCBZXjDB48WHz66afZnqd27drKaxBCiMmTJ4tGjRop37u4uKi1VsyZM0d4enoKIdRbcKZMmSI++OADoVAolPVXr16tjDM+Pl7I5XKxcePGLON5/3inTp3KshXjiy++EB06dFC+Dw4OFhUrVlQ573+NHDlSfPTRR1luP3r0qDAyMhJ3797Nct+sdOzYUYwdO1b5XpMWnGnTpinfJyYmCplMJg4fPiyEEGLixInCw8ND5RxTp05lC04esAWHiEhDrVu3xtq1a1XKypQpAwDo1asXgoODUbFiRbRv3x4dO3ZEly5dYGJiguvXr8PY2DjLVpRM69atw7fffotHjx4hKSkJqampaqOI6tSpAwsLC+V7T09PJCYm4vHjx4iJiUFycjLatm2rsk9qairq1auX7Xn79euHzZs3Y/r06RBCYMeOHcqRP8+fP8fjx48xePBgDB06VLlPWloabGxssjxeeHg4PD09IZPJlGXNmjVDYmIinjx5gujoaKSkpMDb2zvbmDQxdOhQNGrUCE+fPkW5cuUQEhICf39/lfP+l7+/P9q2bYsPPvgA7du3R+fOndGuXTsAwPXr11G+fHlUrVo1y33T09OxcOFChIaG4unTp0hJSUFKSorkTsS1a9dW/rtkyZKwsrJCTExG5/q7d++iUaNGKvUbN24s6fiUP5jgEJHBK1myJCpXrpzlNhcXF9y9exfHjx/HiRMnEBAQgCVLluDMmTMwNzfP8bg//vgjxowZg6VLl8LT0xNWVlZYsmQJLl26pFFcMpkMCkVGp+mDBw+iXLlyKttz6gjdt29fTJo0CVevXkVSUhIeP36MPn36AIDymBs3blR7FJP5+OZ9Qgi1JEP8uyK3TCbL9V5oql69eqhTpw62bt0KHx8f3Lx5E/v378+2fv369REREYHDhw/jxIkT6N27N9q0aYOffvop15iWLl2K5cuXIzg4GLVq1ULJkiUxevRopKamSoq5RIkSKu//+3PL6b6RfjHBIaJiz9zcHF27dkXXrl0xYsQIVKtWDTdv3kStWrWgUChw5swZtGnTRm2/c+fOoWnTpggICFCW3b9/X63eH3/8gaSkJOUX8sWLF2FpaYny5cujdOnSkMvliIyMzLGl6H3ly5dHy5YtsW3bNiQlJaFNmzYoW7YsAKBs2bIoV64cHjx4gH79+ml0vBo1amDXrl0qX9gXLlyAlZUVypUrB3t7e5ibmyMsLAxDhgzJ9XimpqYAMlpR3jdkyBAsX74cT58+RZs2beDi4pLjsaytreHr6wtfX1/07NkT7du3x8uXL1G7dm08efIEf/31V5atOOfOncPHH3+Mzz77DEBG4vf333+jevXqucavqWrVquHQoUMqZb///rvOjk/aY4JDRAYvJSUF0dHRKmUmJiaws7PDli1bkJ6ejg8//BAWFhb47rvvYG5uDjc3N9ja2sLPzw+DBg1SdjJ+9OgRYmJi0Lt3b1SuXBlbt27F0aNH4e7uju+++w6XL1+Gu7u7yrlSU1MxePBgTJs2DY8ePUJQUBC+/PJLGBkZwcrKCuPGjcOYMWOgUCjQvHlzxMfH48KFC7C0tMyyQ3Cmfv36YebMmUhNTcXy5ctVts2cOROjRo2CtbU1OnTogJSUFPz+++949eoVAgMD1Y4VEBCA4OBgjBw5El9++SXu3r2LoKAgBAYGwsjICGZmZpg4cSImTJgAU1NTNGvWDM+fP8etW7cwePBgteO5ublBJpPhwIED6NixI8zNzWFpaamMe9y4cdi4cSO2bt2a489u+fLlcHJyQt26dWFkZISdO3fC0dERpUqVgpeXF1q2bIkePXpg2bJlqFy5Mu7cuQOZTIb27dujcuXK2LVrFy5cuIDSpUtj2bJliI6O1mmCM2zYMCxbtgwTJ07E4MGDcf36dWzZsgUAsn3sRgWDCQ4R5VlsUiy8d+atb4am59HGkSNH4OTkpFL2wQcf4M6dOyhVqhQWLlyIwMBApKeno1atWti/fz9sbW0BAGvXrsWUKVMQEBCAFy9ewNXVFVOmTAEADB8+HNevX4evry9kMhk+/fRTBAQE4PDhwyrn8vb2RpUqVdCyZUukpKSgT58+ypFJADBnzhw4ODhgwYIFePDgAUqVKoX69esrz5OdXr16YeTIkTA2Nka3bt1Utg0ZMgQWFhZYsmQJJkyYgJIlS6JWrVrZztBbrlw5HDp0COPHj0edOnVQpkwZZVKWafr06TAxMcGMGTPw7NkzODk5Yfjw4dkeb9asWZg0aRIGDhyIAQMGKL/4ra2t0aNHDxw8eFAt7vdZWlpi0aJF+Pvvv2FsbIxGjRrh0KFDMDLKWGlo165dGDduHD799FO8efMGlStXxsKFC5XxRkREwMfHBxYWFvj888/RrVs3xMXF5XhOKdzd3fHTTz9h7Nix+Oabb+Dp6YmpU6fiiy++4FxLeiYTxfBhYXx8PGxsbBAXFwdra2t9h0NUJCQnJyMiIgLu7u4wMzMDAHjv9NbLTMYOFg4I6xVW4OfVhr+/P16/fo29e/fqO5RCpW3btqhevTpWrFih71B0bt68eVi3bp1yODtJk9VnTSYp399swSEirRXE7MWF6byUdy9fvsSxY8dw8uRJrFq1St/h6MSaNWvQqFEj2Nra4vz581iyZAm+/PJLfYdV7DHBISKtccFLkqp+/fp49eoVFi1ahA8++EDf4ejE33//jblz5+Lly5dwdXXF2LFjMXnyZH2HVezxERUfURFpJKdmYyIiXdHVIyqj/AxSU2vWrFFeSIMGDXDu3Lkc62/btk05cZaTkxMGDhyIFy9eFFC0REREVNjpPcEJDQ3F6NGjMXXqVFy7dg0tWrRAhw4dEBkZmWX9X375BQMGDMDgwYNx69Yt7Ny5E5cvX9ZoXgYiyrti2OhLRAVIV58xek9wli1bhsGDB2PIkCGoXr06goOD4eLiojateqaLFy+iQoUKGDVqFNzd3dG8eXMMGzaMEysR5bPMGXClzgJLRCRF5mdMdrNua0qvnYxTU1Nx5coVTJo0SaW8Xbt2uHDhQpb7NG3aFFOnTsWhQ4fQoUMHxMTE4KeffkKnTp2yPU/m+iOZ4uPjdXMBRMWIiYkJLCws8Pz5c5QoUUI5DwkRka4oFAo8f/4cFhYWMDHJW4qi1wQnNjYW6enpyunFM5UtW1Zt1tFMTZs2xbZt2+Dr64vk5GSkpaWha9euWLlyZbbnWbBgAWbNmqXT2ImKG5lMBicnJ0RERODRo0f6DoeIDJSRkRFcXV3zPBN0oRgmntVCZdld2O3btzFq1CjMmDEDPj4+iIqKwvjx4zF8+HBs2rQpy30mT56sMjV5fHx8rmufEJE6U1NTVKlShY+piCjfmJqa6qSFWK8Jjp2dHYyNjdVaa2JiYtRadTItWLAAzZo1w/jx4wFkLGNfsmRJtGjRAnPnzlWbjh3IWJGXU2YT6UbmukRERIWZXh+im5qaokGDBjh+/LhK+fHjx9G0adMs93n79q1aZpfZEYmjO4iIiAgoBKOoAgMD8e2332Lz5s0IDw/HmDFjEBkZqVzAbfLkyRgwYICyfpcuXbB7926sXbsWDx48wPnz5zFq1Cg0btwYzs7O+roMIiIiKkT03gfH19cXL168wOzZsxEVFQUPDw8cOnQIbm5uAICoqCiVOXH8/f2RkJCAVatWYezYsShVqhQ++ugjLFq0SF+XQERERIUMl2rgUg1ERERFQpFbqoGIiIhIl5jgEBERkcFhgkNEREQGhwkOERERGRwmOERERGRwmOAQERGRwWGCQ0RERAaHCQ4REREZHCY4REREZHCY4BAREZHBYYJDREREBocJDhERERkcJjhERERkcJjgEBERkcFhgkNEREQGhwkOERERGRwmOERERGRwmOAQERGRwWGCQ0RERAaHCQ4REREZHCY4REREZHCY4BAREZHBYYJDREREBocJDhERERkcJjhERERkcJjgEBERkcGRnOCkpqbmRxxEREREOiM5wSlXrhwmT56MyMjI/IiHiIiIKM8kJzhdunTBihUrUKlSJXzyyScICwvLj7iIiIiItCY5wdm8eTOePHmCefPm4Y8//kC7du1QvXp1rFq1CgkJCfkRIxEREZEkWnUyLl26NCZMmID79+9jz549cHFxwVdffYVy5crhyy+/xJ07d3QdJxEREZHG8jSKSiaToWvXrli0aBG8vLyQmJiINWvWoGbNmujRowdiYmJ0FScRERGRxrROcNLS0rBjxw40b94cDRs2xIMHD7Bo0SI8fPgQwcHBOHfuHAYMGKDLWImIiIg0YiJ1h6dPn2L9+vXYuHEj/vnnH7Ro0QI//vgjPvnkExgZZeRLI0eORLly5fDZZ5/pPGAiIiKi3EhOcCpUqAATExP06dMHX331FerWrZtlvYoVK6Js2bJ5jY+IiIhIMskJTlBQEIYNGwZ7e/sc69WtWxcRERFaB0ZERESkLcl9cFxdXZWPot738uVLbN26Nc9BEREREeWF5ARn4MCBuH//fpbbIiIiMHDgwDwHRURERJQXkhMcIUS225KTk2FsbJyngIiIiIjySqM+OJGRkXj48KHy/bVr15CcnKxSJykpCRs2bICrq6tOAyQiIiKSSqMEJyQkBLNmzYJMJoNMJkNAQIBancyWnW+++Ua3ERIRERFJpFGC07t3b3h4eEAIgd69e2P+/PmoUqWKSh25XA4PDw9UqFAhP+IkIiIi0phGCU716tVRvXp1ABmtOZ07d4atrW2+BkZERESkLcnz4Pj5+eVHHEREREQ6o1GCM3v2bAwZMgTOzs6YPXt2jnVlMhmmT5+uk+CIiIiItCETOY37/peRkREuXryIxo0bZzvJn/KAMhnS09N1FmB+iI+Ph42NDeLi4mBtba3vcIiIiEgDUr6/NWrBUSgUWf6biIiIqDCSPNEfERERUWEnOcFJTk5GfHy8StmPP/6ISZMm4cSJEzoLjIiIiEhbkhOc/v37Y9SoUcr3K1asQJ8+fbB48WL4+Pjg0KFDOg2QiIiISCrJCc5vv/2G9u3bK9+vWLECn332GV6/fo3u3bvj66+/1mmARERERFJJTnCeP3+OcuXKAchYPfzBgwcYOXIkrK2tMXjwYPz55586D5KIiIhICskJjoWFBeLi4gAA586dg6WlJRo2bAgAMDMzQ2Jiom4jJCIiIpJI8kzGtWrVwurVq+Hm5oY1a9agdevWkMlkADJWHXd0dNR5kERERERSSE5wpk+fjs6dO6Nu3bowNTVVGTl18OBB1K9fX6cBEhEREUklOcH56KOPEB4ejitXrqBu3bqoWLGiyra6devqMj4iIiIiyTRaqsHQcKkGIiKiokfnSzVkJSYmBo8ePUJSUpLatpYtW2p7WCIiIqI8k5zgREVFoX///jh16pTaNiFEkVhsk4iIiAyb5ATnyy+/xLVr17Bo0SLUrl0bcrk8P+IiIgBY7wUkxmhe39IBGHYm/+IhIioiJCc4Z86cwddff42BAwfmRzxE9F+JMUDCM31HQURU5EhOcGQyGVxcXPIjFiLKjswIsMxhjqnEaEAoCi4eIqJCTnKC06tXLxw4cABt2rTJj3iIKCuWjsDY8Oy3L63Olh4iov+QvFRD7969cfDgQYwaNQrHjx/H1atX1V5SrVmzBu7u7jAzM0ODBg1w7ty5HOunpKRg6tSpcHNzg1wuR6VKlbB582bJ5yUiIiLDpNVEfwCwatUqrF69WmWbNqOoQkNDMXr0aKxZswbNmjXD+vXr0aFDB9y+fRuurq5Z7tO7d2/8888/2LRpEypXroyYmBikpaVJvRQiIiIyUJITnJCQEJ0GsGzZMgwePBhDhgwBAAQHB+Po0aNYu3YtFixYoFb/yJEjOHPmDB48eIAyZcoAACpUqKDTmIiIiKhok5zg+Pn56ezkqampuHLlCiZNmqRS3q5dO1y4cCHLffbt24eGDRti8eLF+O6771CyZEl07doVc+bMgbm5eZb7pKSkICUlRfk+Pj5eZ9dAREREhY/WMxkDwN27dxEbG4u6deuiZMmSkvePjY1Feno6ypYtq1JetmxZREdHZ7nPgwcP8Msvv8DMzAx79uxBbGwsAgIC8PLly2z74SxYsACzZs2SHB8REREVTZI7GQPA1q1bUb58edSoUQMtW7bE3bt3AWT0jdm4caPk48lkMpX3mX15sqJQKCCTybBt2zY0btwYHTt2xLJly7Bly5Ysl40AgMmTJyMuLk75evz4seQYiYiIqOiQ3IKzc+dO+Pv7o3PnzujQoQNGjBih3Fa/fn38+OOPGDp0qEbHsrOzg7GxsVprTUxMjFqrTiYnJyeUK1cONjY2yrLq1atDCIEnT56gSpUqavvI5XLOuEza42zCRERFjuQWnAULFmDgwIHYt28fPv/8c5Vt1atXx+3btzU+lqmpKRo0aIDjx4+rlB8/fhxNmzbNcp9mzZrh2bNnSExMVJb99ddfMDIyQvny5SVcCZGGMmcT1vQlJRkiIqJ8IbkFJzw8HIsWLcpyW5kyZfDixQtJxwsMDET//v3RsGFDeHp6YsOGDYiMjMTw4cMBZDxeevr0KbZu3QoA6Nu3L+bMmYOBAwdi1qxZiI2Nxfjx4zFo0KBsOxkT6YSuZhOW0iKUmHVfNCIiypnkBMfCwgJxcXFZbnv69ClKly4t6Xi+vr548eIFZs+ejaioKHh4eODQoUNwc3MDkLF6eWRkpLK+paUljh8/jpEjR6Jhw4awtbVF7969MXfuXKmXQiSNrmYT5vpSRET5TnKC06xZM6xatQo9evRQ27Zlyxa0atVKchABAQEICAjIctuWLVvUyqpVq6b2WIuoyMmtRei/LB3yNxYiIgMjOcGZMWMGmjdvjsaNG6Nv376QyWTYvXs3goKCcPbsWfz222/5ESeR4cmtRYiIiLQmuZNxw4YNcfjwYSQmJmLs2LEQQmD+/Pn466+/cOjQIXh4eORHnEREREQa02qiv9atWyM8PBz379/HP//8Azs7O1StWlXXsRERERFpJU8zGVeqVAmVKlXSVSxEREREOiEpwXn+/DnWr1+Ps2fP4tmzjFEgzs7OaN26NT7//HPY2trmS5BERNrwPeCL2KRYyfvZmdshtHNoPkRERAVF4wQnLCwMPXr0QHx8PIyNjWFnZwchBO7evYsTJ07g66+/xp49e9CyZcv8jJdIdzSdj4Zz0RRZsUmxiHnLiReJiiONEpznz5/D19cXNjY2+Pbbb9GxY0dYWFgAAN6+fYsDBw5g3Lhx6NmzJ8LDw9mSQ0UD56MpNoxkRrAzt8u1XmxSLBSaTNZIRIWeRgnOpk2bkJ6ejvPnz6sth2BhYYHevXujSZMmqFOnDjZt2oQJEybkS7BE+ULT+Wg4F02RZWduh7BeYbnW897pzRYfIgOhUYJz7NgxDBo0KMe1nlxdXTFw4EAcOXKECQ4VLZyPhojI4Gg0D054eDiaN2+ea70WLVogPJxfFERERKRfGiU4r1+/hoND7s3zDg4OeP36dV5jIiIiIsoTjRKclJQUlChRItd6JiYmSE1NzXNQRERERHmh8TDxu3fvwsQk5+p37tzJc0BEREREeaVxguPv759rHSEEZDJZXuIhIiIiyjONEpyQkJD8joOIiIhIZzRKcPz8/PI7DiIiIiKd0aiTMREREVFRwgSHiIiIDI6k1cSJigQuoklEVOwxwSHDw0U0iYiKPSY4ZLi4iCYRUbElOcGJjo6Go6MGXxpE+sZFNImIii3JnYxdXV3x6aef4vz58/kRDxEREVGeSU5wpk2bhnPnzqFly5aoW7cuNm3ahKSkpPyIjYiIiEgrkhOcGTNm4NGjR9ixYwesra0xdOhQlC9fHuPGjcP9+/fzI0YiIiIiSbTqZGxsbIzevXujd+/euHHjBlatWoV169YhODgY7du3x8iRI+Hj46PrWImICkRsUiy8d3pL2sfO3A6hnUPzKSIikirPo6hq1aqFDh064MaNG/jtt98QFhaGw4cPo169eti+fTuqVq2qiziJiAqMQigQ81aDuZSIqNDSeibj2NhYLFiwAO7u7ujZsydMTEwQGhqK+Ph47N27FwkJCRqtQE5EVFjYmdvBwcJB0stIxgnhiQojyS04ly5dwurVq7Fz504IIeDr64uvvvoK9evXV9bp0qULTExM0K1bN13GSkS5SYwGllbPvZ6lAzDsTP7HU8Ro84jJe6c3W3uICiHJCY6npyccHR0xadIkfPHFF3BwyHqStAoVKqBp06Z5DpCIJBAKzuJMRAQtEpytW7fC19cXJUqUyLFe9erVcerUKa0DIyIJNJ2NOTE6IwkiIjJwkhOcBw8e4Pnz53B2dlbbFhUVhY0bN2LGjBk6CY6INKTp46al1dnCQ0TFguTecbNmzcKTJ0+y3Pbs2TPMmjUrz0ERERER5YXkFhwhRLbbEhMTc310RWTwcuvomxhdcLFQoeN7wBexSbGS9+M8O0TSaJTg3LhxA9evX1e+P3ToEO7cuaNSJykpCdu2bUOlSpV0GiBRkcOOvpSD2KRYjroiKgAaJTh79uxRPnqSyWSYPXt2lvXMzc0REhKiu+iIihJNO/pqW58MipHMCHbmdrnWi02KhYIdw4kk0yjB+fzzz9G5c2cIIdC4cWOEhITAw8NDpY5cLkelSpVgbm6eL4ESFXqcV4YksDO3Q1ivsFzrcZ4dIu1olOA4OTnByckJAHDq1CnUr18fVlZW+RoYERERkbYkdzL28vLKjziIiIiIdEajBGfQoEGYPn063N3dMWjQoBzrymQybNq0SSfBERFl0mb0kTajlYjIMGiU4Jw6dQpfffUVAODkyZOQyWTZ1s1pGxGRtjj6iIik0CjBiYiIUP774cOH+RULEVGuNB199F9S6xNR0Se5Dw4RkT5pOvqIiIo3yUs1EBERERV2GrXguLu7a9y3RiaT4f79+3kKioiIiCgvNEpwvLy82HmYiIiIigyNEpwtW7bkcxhEhUeXlb/geUKK1vvbW8mxf2RzHUZERERSsZMx0XueJ6QgOj5Z32FQERObFAvvnd4a1SOi/KdRghMZGQknJyeUKFECkZGRudZ3dXXNc2BE+mYkAxyszDSuH5OQDIXIx4CoUFMIBefpISpENO5k/Ouvv6Jx48aoUKFCrv1x0tPTdRIckT45WJnh4pTc/yLP1GR+WIG2/GjzKG1fajIckJGMDV75Cx+l6YC2c+xwbh6i/KVRgrN582ZUqlRJ+W92OCbSP20epaXLAciAdIE89TOi/xfaOVTfIRBRFjRKcPz8/JT/9vf3z69YiEgLUh6lGafmczBERIVEnjoZCyGQmJgIS0tLtuoQ6YmkR2lLzYCE/I2HiKgw0Gom40uXLsHHxwcWFhYoVaoULCws4OPjg4sXL+o6PiIiIiLJJLfgnDx5Eh06dICVlRX69OkDR0dHREdHY//+/fDy8sKhQ4fg7a15x0wiIiIiXZOc4EycOBH16tXDiRMnYGlpqSxPSEiAt7c3Jk2ahMuXL+s0SCIiIiIpJCc4f/75J7Zt26aS3ACAlZUVJk6ciM8++0xnwRFR/nDAK+xLHZLRJycnlg7AsDMFExQRkQ5JTnAcHBxgZJR11x1jY2PY29vnOSiioiwmIRlN5odJ3q8gl3gwlgk44CU7HBORwZKc4AwbNgzLly9Hp06dUKJECWV5amoqli1bhs8//1ynARIVNQqBwrvUg6UDYhKSkS4A45yGlydGA0JRsLEREemQRgnOsmXLlP82NTXFw4cPUbFiRXTv3l3ZyXj37t0wNjaGubl5vgVLVJjZW8m12q9Al3gYdgZd/51x2dHaDBfHZjMgYGl1IOFZAQVFRKR7GiU448aNy7J85cqVamUTJkzA2LFj8xYVURGk7eOlzCUepD7aikkopK1ERESFgEYJTkRERH7HQVTsFepHW0RERYxGCY6bm1t+x0FUbGn7aEtX+xMRGaI8LdWgK2vWrMGSJUsQFRWFmjVrIjg4GC1atMh1v/Pnz8PLywseHh64fv16/gdKlA+4ojcRke5pleCcPXsWK1asQHh4OJKSklS2yWQy3L9/X+NjhYaGYvTo0VizZg2aNWuG9evXo0OHDrh9+zZcXV2z3S8uLg4DBgyAt7c3/vnnH20ug4gk6LLylzyvQF6QQ+GJqHiTvBbVL7/8Am9vb8TFxSE8PBzVqlVDuXLlEBkZCRMTE7Rs2VLS8ZYtW4bBgwdjyJAhqF69OoKDg+Hi4oK1a9fmuN+wYcPQt29feHp6Sr0EItLC84QURMcn5+mV1wSJiEhTkhOcoKAgDBw4EEeOHAEAzJ07F+fOncPVq1eRmJiI7t27a3ys1NRUXLlyBe3atVMpb9euHS5cuJDtfiEhIbh//z6CgoI0Ok9KSgri4+NVXkSkHSMZ4GhtJullJNN31ERU3Gi1VMO4ceMgk2V8YqWnpwMAateujenTp2P27Nno0qWLRseKjY1Feno6ypYtq1JetmxZREdHZ7nP33//jUmTJuHcuXMwMdEs/AULFmDWrFka1SWinDlYmeHiFGkL6mYOhSciKiiSW3Devn0LS0tLGBkZQS6XIzY2VrmtWrVquH37tuQgMpOlTEIItTIgI5nq27cvZs2ahapVq2p8/MmTJyMuLk75evz4seQYiYiIqOiQ3ILj6uqq7NRbo0YNHDx4EB06dAAAnDlzBra2thofy87ODsbGxmqtNTExMWqtOkDGiuW///47rl27hi+//BIAoFAoIISAiYkJjh07ho8++khtP7lcDrmcQ2mJiIiKC8kJTqtWrXD69Gn07NkTQ4cORUBAAMLDwyGXy3Hs2DFJsxibmpqiQYMGOH78OD755BNl+fHjx/Hxxx+r1be2tsbNmzdVytasWYOTJ0/ip59+gru7u9TLISIqEmKTYuG9U9qjQTtzO4R2Ds2niIgKN8kJzqxZs/Dy5UsAwPDhw/H27Vts27YNMpkM06ZNw9SpUyUdLzAwEP3790fDhg3h6emJDRs2IDIyEsOHDweQ8Xjp6dOn2Lp1K4yMjODh4aGyv4ODA8zMzNTKiYgMiUIoEPM2Rt9hEBUZkhMcOzs72NnZKd8HBgYiMDBQ6wB8fX3x4sULzJ49G1FRUfDw8MChQ4eUsydHRUUhMjJS6+MTERVlduZ2uVd6T2xSLBRcDZ6KuTzNZPzs2TO8ePECtra2cHZ21vo4AQEBCAgIyHLbli1bctx35syZmDlzptbnJiIqzLR5xOS905utPVTsSR5FBQC7d+/GBx98ABcXF9StWxcuLi6oWrUqfvrpJ13HR0RERCSZ5AQnNDQUPXv2hLGxMWbMmIE1a9Zg+vTpMDY2hq+vL0JD2aGNiIiI9EvyI6rZs2ejQ4cO2L9/P4yM/j8/mjFjBjp16oTZs2fD19dXp0ESERERSSE5wbl//z4WL16sktwAgJGREQICAtCrVy+dBUdE+ScmIRlN5odluW1fajIc/q3T9d86MQmciZiIig7JCY6bmxvevn2b5ba3b9/CxcUlz0ERUf5TCGS7fEK6HIAMSM+hDhFRYSY5wRk7dixmz56N1q1bqwwXj4mJwdy5czFu3DidBkhEumVvlfus3sapGf91kL3Gb2ZfqmwzSpMBS987hqUDMOyMpDh8D/giNik294r/klKXiEijBGfUqFEq7+Pj41GhQgV4e3vD0dER0dHRCAsLg52dnVZrURFRwdk/snnulZaaAQmAMRRwwEvVbQoACXmPIzYplkOZiSjfaJTgrFq1Ksvy/fv3q7yPjIzEqlWr8M033+Q9MiLSH0sHzeolRgN5nFDOSGYkaTI7bSa+I6LiR6MER6HgjJhExYqmj5uWVgcSnuXpVHbmdgjrlXVnZyIibWk10R8RERFRYab1Ug1hYWEICwvDixcvYGdnB29vb3z00Ue6jI2IiIhIK5ITnNTUVPTo0QOHDh2CEAImJiZIS0vDwoUL0alTJ+zatQslSpTIj1iJiIiINCL5EdXs2bNx9OhRLFy4EP/88w9SU1Pxzz//YNGiRTh69Chmz56dH3ESERERaUxyC86OHTswZcoUjB8/Xllmb2+PcePGITExEVu3bsWcOXN0GiQRAGC9F5CowbDixOj8j4XU5DSvzVvHFJR0EHhrJIP3zvkAOK8NEeUvyQnOkydP0KJFiyy3tWjRAgsWLMhzUERZSozJ84gdyj85zmtjDBgZAwJATNYToRMR6ZTkBMfe3h43b96Et7e32rabN2/C3t5eJ4ERZUtmBFg65l5P07lcSKeymtfmeUIK0hUCxkYytZmUOa8NEeUHyQlO165dMWPGDLi6uqJ79+7K8p9//hkzZ85Ev379dBogkRpLR2BsuL6joGxkNa9Nk/lhiI5PhqO1GcIGqf9xRESka5ITnHnz5uH8+fPo1asXSpYsCUdHR/zzzz9ITExErVq1MG/evPyIk4iIiEhjkhOc0qVL47fffsOWLVtw6tQpvHjxAvXr14e3tzcGDBgAuTz3hfyIiIiI8pOkBCcpKQlt2rTBrFmzMGzYMAwbNiy/4iIiojyKTYqF907NHwnamdshtHNoPkZEVHAkJTjm5ua4efMmTEy0ngCZiIgKiEIouGI7FVuSMxVPT0/89ttvaNWqVT6EQ0REeSV1ZFpsUiwUeVwVnqiwkZzgLF26FB9//DEcHR3RvXt3WFpa5kdcRESkJamPmbx3erOlhwyO5KUaPD098eTJEwwcOBA2NjawsrKCtbW18mVjY5MfcRIRERFpTHILTo8ePSCTyfIjFiIiIiKdkJzgbNmyJR/CICIiItIdjROcpKQk7N27F48ePYKDgwO6dOnCZRmIiIioUNIowXn27BlatmyJiIgICCEAADY2Njh8+DCaNGmSrwESERERSaVRgjNt2jQ8ffoU06ZNQ5MmTfD3339j3rx5+OKLL3Dt2rX8jpGIDERMQjKazA/LveJ77K3k2D+yeT5ERESGSqME5/jx45gyZQqmT58OAOjQoQMqVaqErl274p9//kHZsmXzNUgiMgwKAUTHJ+s7DCIqBjRKcKKjo9GyZUuVslatWkEIwQSHCq0uK3/B84QUyfvFJPALWNfsrbRboy4mIRkKoeNgiKhY0CjBSU9Ph7m5uUqZmZkZACAtLU33URHpwPOEFLYWFBLaPl5qMj+MP0Mi0orGo6ju3r2rsgZVeno6AODOnTtqdevXr6+D0Ih0w0gGOFiZSd5P21YHIiLSP40THH9//yzL+/fvr/y3EAIymUyZ/BAVBg5WZrg4RfMVlYmIqOjTKMEJCQnJ7ziIiIiIdEajBMfPzy+/4yAiIiLSGcmLbRIREREVdkxwiIiIyOBIXmyTSOfWewGJMbnXS4zO/1iIiMggMMEh/UuMARKe6TsKIiIyIExwqPCQGQGWjrnXs3TI/1iIiKhIY4JDhYelIzA2XN9RUCHERTqJSComOERU6HGRTiKSigkOERVaXKSTiLTFBIeICi0u0klE2uI8OERERGRw2IJDRPlP07mOMlk6AMPO5F88RGTwmOAQUf7jXEdEVMCY4BBRwcltrqPEaEAoCi4eUhGbFAvvnd6S9rEzt0No59B8iohIe0xwiKjg5DbX0dLqbOnRI4VQIOathEeJRIUYExwiomLOztxO8j6xSbFQsLWNCjEmOERExZw2j5i8d3qztYcKNSY4RKTG94AvYpNic69Y2hgo5QzIjDWrT0RUQJjgEJGa2KRYzf46N5ZB+THCxxVEVIgwwaFCr8vKX/A8IUXyfjEJnMk2r4xkRjn3z0j4BxDpgMwYsCoLQLv+HEREusYEhwq95wkpnHZfT+zM7RDWK4dVvDNHPVk5A4Okr/ZNRJRfmOBQkWEkAxyszCTvp+2CjSRBYnRGspPTdiKiAsQEh4oMByszXJwibRIyKiBCwflriKhQYYJDRNqzdMjf+nkUk5CMJvOlPzqzt5JrvZI5ERUOTHCISHuFfEFMhQD7bxEVU0xwiMjgaNvvKiYhGQqh42CISC+Y4FCByW64977UZDgg48ulaxaPEzjcm6TS9vFSk/lhbPGRSOoCnVyckwoKExwqMNkN906XA5AB6XycQFTkcIFOKqwKRYKzZs0aLFmyBFFRUahZsyaCg4PRokWLLOvu3r0ba9euxfXr15GSkoKaNWti5syZ8PHxKeCoKVfrvYDE///g25eanJHM4N8JcP9li9fKMkfr7IeBc7g3UeEhdUJHLs5JBU3vCU5oaChGjx6NNWvWoFmzZli/fj06dOiA27dvw9XVVa3+2bNn0bZtW8yfPx+lSpVCSEgIunTpgkuXLqFevXp6uALKVmKMytBhBwCQZVs7Yxj4WA4DJyoKpD5m4uKcVND0nuAsW7YMgwcPxpAhQwAAwcHBOHr0KNauXYsFCxao1Q8ODlZ5P3/+fPz888/Yv38/E5zCSmYEWDoiJiEZ6SKjpSbLCfsKeAgxEREZLr0mOKmpqbhy5QomTZqkUt6uXTtcuHBBo2MoFAokJCSgTJky2dZJSUlBSsr/d26Nj4/XLmDSjqUjMDYcXf/twOlozZYaykVuMyNnsnQo9EPViUg/9JrgxMbGIj09HWXLllUpL1u2LKKjNZvafenSpXjz5g169+6dbZ0FCxZg1qxZeYqViAoQZ0YmojzS+yMqAJDJVDtmCCHUyrKyY8cOzJw5Ez///DMcHLJ/vDF58mQEBgYq38fHx8PFxUX7gIkof2j6mDIxOiMJIiLKhl4THDs7OxgbG6u11sTExKi16rwvNDQUgwcPxs6dO9GmTZsc68rlcsjlHIFDVOhp+rgpcxVzIqJsGOnz5KampmjQoAGOHz+uUn78+HE0bdo02/127NgBf39/bN++HZ06dcrvMImIiKiI0fsjqsDAQPTv3x8NGzaEp6cnNmzYgMjISAwfPhxAxuOlp0+fYuvWrQAykpsBAwbgm2++QZMmTZStP+bm5rCxsdHbdRQb781tk6NEzfpRUf7zPeCL2KRYjetLqUtEVBjpPcHx9fXFixcvMHv2bERFRcHDwwOHDh2Cm5sbACAqKgqRkZHK+uvXr0daWhpGjBiBESNGKMv9/PywZcuWgg6/+HlvbhsqGmKTYjkHCREVK3pPcAAgICAAAQEBWW57P2k5ffp0/gdEuft3bhuNcH6bQsNIZiRpBlqps9USERUWhSLBoSLo37ltqGixM7dDWC/1BU1JVUxCMppksfBrbuyt5Fov9ElEusUEh4joPQou/JpvuPo4FRQmOERE/9J2QdeYhGQohI6DMVBcfZwKChMcIqJ/aft4qcm/y5BQ9rj6OBU0JjhERJTvuPo4FTS9TvRHRERElB/YgkNEpCMcfUVUeDDBISLSEY6+Iio8mOAQEeURR18RFT5McIiI8oijr4gKH3YyJiIiIoPDBIeIiIgMDhMcIiIiMjhMcIiIiMjgMMEhIiIig8MEh4iIiAwOh4kTFTG+B3wRmxQraR+p9YmIijomOERFTGxSLBchJCLKBRMcoiLKSGYEO3M7SftIrU9EVFQxwSEqouzM7RDWS/rCjkRExQE7GRMREZHBYYJDREREBoePqEiS2DcpsEPGKshd50t7PBKTwEUFiYioYDDBIUkUCgEASBfgKshERFRoMcEhrTlam2m1n72VXMeREBERqWKCQ1oxlgEXp3jrOwwiIqIsMcEh0jOpMxNzVmIiotwxwSHSM85MTESke0xwiAoJqTMTc1ZiAInRwNLqmtW1dACGncnfeIio0GCCQ1RIcGZiLQgFkPBM31EQUSHEBIcyrPcCEnN/TGKL1/kfC1FuLB00r5sYnZEIEVGxwgSHMiTGaPSXsHEBhEKUKymPmpZWZysPUTHEBIdUyYwAS8dsN8ckJCNdAK+NSkPC39BEREQFigkOqbJ0BMaGZ7u56/wwRMcnw9HaDBcLMCwiIiIpuNgmERERGRwmOERERGRwmOAQERGRwWEfnGKqy8pf8DwhRfl+X2oyHJDRibjr/OznYolJ4AriRLoWk5CMJjn8f5cdeys59o9sng8RFR6xSbHw3ilt3Ts7czuEdg7Np4ioqGCCU0w9T0hBdPz/JyvpcgAyIF1ApZyI8p9Cy//vouMNPzFSCAWXMiGtMMEp5oxkgIOVGYxTM94bywBHa7Nc97O3kudzZESGT9v/j/6bDBnqHyTaLEUSmxQLBSd1pH8xwSnmHKzMcHGKN7DUDEj49/1Yac3BRKQdbVtR3n/ErKmYhGQohFanLHDaPGLy3unN1h5SYoJDRFTEaJsYNfl3Hiui4oAJThGXl7/kiIiIDBUTnCLu/c7CRJSNxOiMdal0xdJB2ppYVGCkjrziqCvDxATHQGR2FpaKnYWp2BAKLrpZTHDkFQFMcAyGsrMw6Z3vAV/EJsVqXF9KXdKCpY6XhU2MzkiWqNCROvKKo64MGxMcIh2LTYrlX4+Fia4fIy2tzpagQkrqYyaOujJsTHCI8omRzEjSX5TazPtBRHmXn7MlS23RlXp8yh4THKJ8Ymduh7Be0meZpSJC007L7Ixc6OVnnx226OoPExwiIm2w03KRl5fZkjVt9clsvdG0RZf9gnSHCU5Rtd4LSIzBvtRkpMuRsdTCUumjqJQSo3UWGpFB07TTMjsjF3p5mS1ZaquPpi267BekO0xwiqrEGCDhGRwAQPZvWYIe4yEqLjR93MTOyAZJ275y7GNX8JjgFHHpMEKMKAVjLefBUaPrIbVERAaEHX+LDiY4Oqbt0gmZ7K3kktaZeYFS8ExZBUdrLpJJRESUiQmOjnHpBCIiIv1jgpNPpC6dEJOQDIXIx4CIiIiKESY4+UTq0glN5ochOj4ZIe/GA0uTct+Bo56IiIiyxQSnkCmleAUkvNS4fjpbfSTRdlZRKbi2FBV2MQnJaDJfu0kopfYTJNIXJjiFVLqQIQalc633XNgUQDSGg7OKEgEKAfYVJIPHBKeQsLeSA/h3wj4AL2Sl8Ynptxrt62j2//uTZqSuE6UNzntBhU1ePifYT5CKGiY4hYSyyXepGZDwbx8eDvvON1wnioqjvDxayuwnWFDyOuWGtvgIznAwwSko/y6tkCt2HiYyLFyUUyuccoPyiglOQfl3aQUiKmYMbFHOvHRQlnoeQPqUG3k5Hx/BGZZCkeCsWbMGS5YsQVRUFGrWrIng4GC0aNEi2/pnzpxBYGAgbt26BWdnZ0yYMAHDhw8vwIjzQGYEWDrmXo9LJhAVbYV9UU5NW5X/tS81GdGmNuiaOq9AW1akTrmhrYJ+BEf5T+8JTmhoKEaPHo01a9agWbNmWL9+PTp06IDbt2/D1dVVrX5ERAQ6duyIoUOH4vvvv8f58+cREBAAe3t79OjRQw9XIJGlIzA2XN9RFEocwk0GpbAvyimxVdkBAIwAR+v8b035Lw6gUKWL5YAs3FcW2Gehnbmd3tbv0nuCs2zZMgwePBhDhgwBAAQHB+Po0aNYu3YtFixYoFZ/3bp1cHV1RXBwMACgevXq+P333/H1118XjQSHssUh3ER6oEmr8r+tTBz8oH+66JtUsph81uo1wUlNTcWVK1cwadIklfJ27drhwoULWe7z66+/ol27diplPj4+2LRpE969e4cSJUqo7ZOSkoKUlP/PeOPi4gAA8fHxeb0ENWnJb6BISUFacrrq8ZPTgRQBlEgH8uG8hiDtbRrSk9JhJDOCrZltvp7LRtjky8+fSLLMz4aUKGBe1YI7b2IMAAFY2gNDL+Vcd2XDjCTHgD+/Mj+7o2PfouGMfXqLI6nsGwjjdES9jUbN1U3UtqdbCZhZZfzb2Eim8XHT/+1glGAkQ+LLRCiEIl8/a18kv4BCKJCGNJ1+1mYeS4jcO0zpNcGJjY1Feno6ypYtq1JetmxZREdnPZooOjo6y/ppaWmIjY2Fk5OT2j4LFizArFmz1MpdXFzyEH3OHgOwmZPVlgRgGifnKwz2YI++QyB6T5wezinlM6l4fH490ncABiQc4bDx1/3vTEJCAmxscj6u3h9RAYBMppqFCiHUynKrn1V5psmTJyMwMFD5XqFQ4OXLl7C1tc3xPLmJj4+Hi4sLHj9+DGtra62PY2h4X9TxnmSN9yVrvC9Z431RV9zuiRACCQkJcHZ2zrWuXhMcOzs7GBsbq7XWxMTEqLXSZHJ0dMyyvomJCWxts25qk8vlkMtVO6qVKlVK+8DfY21tXSx+saTifVHHe5I13pes8b5kjfdFXXG6J7m13GQyyuc4cmRqaooGDRrg+PHjKuXHjx9H06ZNs9zH09NTrf6xY8fQsGHDLPvfEBERUfGj1wQHAAIDA/Htt99i8+bNCA8Px5gxYxAZGamc12by5MkYMGCAsv7w4cPx6NEjBAYGIjw8HJs3b8amTZswbtw4fV0CERERFTJ674Pj6+uLFy9eYPbs2YiKioKHhwcOHToENzc3AEBUVBQiIyOV9d3d3XHo0CGMGTMGq1evhrOzM1asWKGXIeJyuRxBQUFqj7+KO94XdbwnWeN9yRrvS9Z4X9TxnmRPJjQZa0VERERUhOj9ERURERGRrjHBISIiIoPDBIeIiIgMDhMcIiIiMjhMcLT09OlTfPbZZ7C1tYWFhQXq1q2LK1eu6DssvalQoQJkMpnaa8SIEfoOTa/S0tIwbdo0uLu7w9zcHBUrVsTs2bOhUCj0HZreJSQkYPTo0XBzc4O5uTmaNm2Ky5cv6zusAnX27Fl06dIFzs7OkMlk2Lt3r8p2IQRmzpwJZ2dnmJubo1WrVrh165Z+gi0gud2T3bt3w8fHB3Z2dpDJZLh+/bpe4ixoOd2Xd+/eYeLEiahVqxZKliwJZ2dnDBgwAM+e6WGV+kKECY4WXr16hWbNmqFEiRI4fPgwbt++jaVLl+p0duSi5vLly4iKilK+Midj7NWrl54j069FixZh3bp1WLVqFcLDw7F48WIsWbIEK1eu1HdoejdkyBAcP34c3333HW7evIl27dqhTZs2ePr0qb5DKzBv3rxBnTp1sGrVqiy3L168GMuWLcOqVatw+fJlODo6om3btkhISCjgSAtObvfkzZs3aNasGRYuXFjAkelXTvfl7du3uHr1KqZPn46rV69i9+7d+Ouvv9C1a1c9RFqICJJs4sSJonnz5voOo1D76quvRKVKlYRCodB3KHrVqVMnMWjQIJWy7t27i88++0xPERUOb9++FcbGxuLAgQMq5XXq1BFTp07VU1T6BUDs2bNH+V6hUAhHR0excOFCZVlycrKwsbER69at00OEBe/9e/JfERERAoC4du1agcZUGOR0XzL99ttvAoB49OhRwQRVCLEFRwv79u1Dw4YN0atXLzg4OKBevXrYuHGjvsMqNFJTU/H9999j0KBBeVrM1BA0b94cYWFh+OuvvwAAf/zxB3755Rd07NhRz5HpV1paGtLT02FmZqZSbm5ujl9++UVPURUuERERiI6ORrt27ZRlcrkcXl5euHDhgh4jo6IgLi4OMpmsWD9ZYIKjhQcPHmDt2rWoUqUKjh49iuHDh2PUqFHYunWrvkMrFPbu3YvXr1/D399f36Ho3cSJE/Hpp5+iWrVqKFGiBOrVq4fRo0fj008/1XdoemVlZQVPT0/MmTMHz549Q3p6Or7//ntcunQJUVFR+g6vUMhcVPj9hYfLli2rtuAw0X8lJydj0qRJ6Nu3b7FZgDMrel+qoShSKBRo2LAh5s+fDwCoV68ebt26hbVr16qsm1Vcbdq0CR06dNBoOXtDFxoaiu+//x7bt29HzZo1cf36dYwePRrOzs7w8/PTd3h69d1332HQoEEoV64cjI2NUb9+ffTt2xdXr17Vd2iFyvutoEKIYt8yStl79+4d+vTpA4VCgTVr1ug7HL1iC44WnJycUKNGDZWy6tWrq6yZVVw9evQIJ06cwJAhQ/QdSqEwfvx4TJo0CX369EGtWrXQv39/jBkzBgsWLNB3aHpXqVIlnDlzBomJiXj8+DF+++03vHv3Du7u7voOrVBwdHQEALXWmpiYGLVWHSIgI7np3bs3IiIicPz48WLdegMwwdFKs2bNcPfuXZWyv/76S7lAaHEWEhICBwcHdOrUSd+hFApv376FkZHq/2bGxsYcJv4fJUuWhJOTE169eoWjR4/i448/1ndIhYK7uzscHR2VIxKBjP5tZ86cQdOmTfUYGRVGmcnN33//jRMnTsDW1lbfIekdH1FpYcyYMWjatCnmz5+P3r1747fffsOGDRuwYcMGfYemVwqFAiEhIfDz84OJCX+1AKBLly6YN28eXF1dUbNmTVy7dg3Lli3DoEGD9B2a3h09ehRCCHzwwQe4d+8exo8fjw8++AADBw7Ud2gFJjExEffu3VO+j4iIwPXr11GmTBm4urpi9OjRmD9/PqpUqYIqVapg/vz5sLCwQN++ffUYdf7K7Z68fPkSkZGRyjleMv/YdHR0VLZ6GaKc7ouzszN69uyJq1ev4sCBA0hPT1e2/JUpUwampqb6Clu/9D2Mq6jav3+/8PDwEHK5XFSrVk1s2LBB3yHp3dGjRwUAcffuXX2HUmjEx8eLr776Sri6ugozMzNRsWJFMXXqVJGSkqLv0PQuNDRUVKxYUZiamgpHR0cxYsQI8fr1a32HVaBOnTolAKi9/Pz8hBAZQ8WDgoKEo6OjkMvlomXLluLmzZv6DTqf5XZPQkJCstweFBSk17jzW073JXPIfFavU6dO6Tt0vZEJIUQB5VJEREREBYJ9cIiIiMjgMMEhIiIig8MEh4iIiAwOExwiIiIyOExwiIiIyOAwwSEiIiKDwwSHiIiIDA4THCIiIjI4THCIiIjI4DDBIZJgy5YtkMlkkMlkOH36tNp2IQQqV64MmUyGVq1aqe338OFDrc+pzb6F3Y0bNzBw4EC4u7vDzMwMlpaWqF+/PhYvXoyXL1/qOzwlbX/uhYW/v78yfg8PD2W5Pq5r7969ynPKZDL8/vvvOjku0fuY4BBpwcrKCps2bVIrP3PmDO7fvw8rKyuV8k6dOuHXX3+Fk5OT5HPlZd/CbOPGjWjQoAEuX76M8ePH48iRI9izZw969eqFdevWYfDgwfoOUY3Un3th4ujoiF9//RXbt29X21aQ1+Xl5YVff/0V06ZN09kxibLCJZ+JtODr64tt27Zh9erVsLa2VpZv2rQJnp6eiI+PV6lvb28Pe3t7rc6Vl30Lq19//RVffPEF2rZti71790Iulyu3tW3bFmPHjsWRI0d0cq63b9/CwsJCJ8eS+nMvTORyOZo0aZLltoK8rtKlS6NJkya4c+eOzo5JlBW24BBp4dNPPwUA7NixQ1kWFxeHXbt2YdCgQWr1s3rMNHPmTMhkMty6dQuffvopbGxsULZsWQwaNAhxcXEa7Xvjxg306tULNjY2KFOmDAIDA5GWloa7d++iffv2sLKyQoUKFbB48WKVePz9/VGhQgW1ODOPm1WZtufKyvz58yGTybBhwwaV5CaTqakpunbtqlL2yy+/wNvbG1ZWVrCwsEDTpk1x8ODBLGO9evUqevbsidKlS6NSpUrK7X///Tf69u0LBwcHyOVyVK9eHatXr8413kxSf+4AcO/ePQwcOBBVqlSBhYUFypUrhy5duuDmzZsq9Z4/f47PP/8cLi4ukMvlsLe3R7NmzXDixAlJdbShzXURFXZMcIi0YG1tjZ49e2Lz5s3Ksh07dsDIyAi+vr6SjtWjRw9UrVoVu3btwqRJk7B9+3aMGTNGo3179+6NOnXqYNeuXRg6dCiWL1+OMWPGoFu3bujUqRP27NmDjz76CBMnTsTu3bslxZVf50pPT8fJkyfRoEEDuLi4aHTuM2fO4KOPPkJcXBw2bdqEHTt2wMrKCl26dEFoaKha/e7du6Ny5crYuXMn1q1bBwC4ffs2GjVqhD///BNLly7FgQMH0KlTJ4waNQqzZs3SKA5tfu7Pnj2Dra0tFi5ciCNHjmD16tUwMTHBhx9+iLt37yrr9e/fH3v37sWMGTNw7NgxfPvtt2jTpg1evHghqY42dPn7TFRoCCLSWEhIiAAgLl++LE6dOiUAiD///FMIIUSjRo2Ev7+/EEKImjVrCi8vL7X9IiIilGVBQUECgFi8eLHKOQICAoSZmZlQKBS57rt06VKVfevWrSsAiN27dyvL3r17J+zt7UX37t2VZX5+fsLNzU3t+jKPm1WZtud6X3R0tAAg+vTpk22d9zVp0kQ4ODiIhIQEZVlaWprw8PAQ5cuXV96rzFhnzJihdgwfHx9Rvnx5ERcXp1L+5ZdfCjMzM/Hy5ctsz6/tzz0raWlpIjU1VVSpUkWMGTNGWW5paSlGjx6d476a1MlKdj9vXV6XVP89N1F+YAsOkZa8vLxQqVIlbN68GTdv3sTly5e1as5//1FM7dq1kZycjJiYmFz37dy5s8r76tWrQyaToUOHDsoyExMTVK5cGY8ePZIcm77O9V9v3rzBpUuX0LNnT1haWirLjY2N0b9/fzx58kSlJQTIaBX7r+TkZISFheGTTz6BhYUF0tLSlK+OHTsiOTkZFy9e1CgeqT/3tLQ0zJ8/HzVq1ICpqSlMTExgamqKv//+G+Hh4cp6jRs3xpYtWzB37lxcvHgR7969UzuWJnW0lZff5+joaIwaNQr169eHl5cX1qxZA4VCobPYiLTBBIdISzKZDAMHDsT333+PdevWoWrVqmjRooXk49ja2qq8z+yTkpSUlOu+ZcqUUXlvamoKCwsLmJmZqZUnJydLji0/zmVnZwcLCwtERERodN5Xr15BCJHlKDJnZ2cAUHtE837dFy9eIC0tDStXrkSJEiVUXh07dgQAxMbGahSP1J97YGAgpk+fjm7dumH//v24dOkSLl++jDp16qj8jENDQ+Hn54dvv/0Wnp6eKFOmDAYMGIDo6GhJdbSl7e/z06dP8dFHH6Fu3bqYOHEiJk+ejCtXrmDAgAF5jokoL5jgEOWBv78/YmNjsW7dOgwcOFDf4WjMzMwMKSkpauWafsnnhbGxMby9vXHlyhU8efIk1/qlS5eGkZERoqKi1LY9e/YMQEbS9F/vd5QuXbo0jI2N4e/vj8uXL2f5ykx0NCHl5/79999jwIABmD9/Pnx8fNC4cWM0bNhQ7V7b2dkhODgYDx8+xKNHj7BgwQLs3r0b/v7+kurkhTa/zxMnTsSUKVMwaNAgHD58GK9fv8amTZvw4sULHD16VCdxEWmDCQ5RHpQrVw7jx49Hly5d4Ofnp+9wNFahQgXExMTgn3/+UZalpqYW2BfS5MmTIYTA0KFDkZqaqrb93bt32L9/PwCgZMmS+PDDD7F7926VFg+FQoHvv/8e5cuXR9WqVXM8n4WFBVq3bo1r166hdu3aaNiwodrr/Za0nEj5uctkMrWRYgcPHsTTp0+z3cfV1RVffvkl2rZti6tXr2pdRyptfp9Pnz6tHIX1X35+fjh8+LBO4iLSBufBIcqjhQsX6jsEyXx9fTFjxgz06dMH48ePR3JyMlasWIH09PQCOb+npyfWrl2LgIAANGjQAF988QVq1qyJd+/e4dq1a9iwYQM8PDzQpUsXAMCCBQvQtm1btG7dGuPGjYOpqSnWrFmDP//8Ezt27FBrscnKN998g+bNm6NFixb44osvUKFCBSQkJODevXvYv38/Tp48KekaNP25d+7cGVu2bEG1atVQu3ZtXLlyBUuWLEH58uWVdeLi4tC6dWv07dsX1apVg5WVFS5fvowjR46ge/fuGtfRBW1+nwcPHowffvgB7969w/bt2+Hv749OnTrB0dFRZ3ERScUEh6gYcnd3x88//4wpU6agZ8+ecHJyQmBgIJ4/f67xkOm8Gjp0KBo3bozly5dj0aJFiI6ORokSJVC1alX07dsXX375pbKul5cXTp48iaCgIPj7+0OhUKBOnTrYt2+fWufn7NSoUQNXr17FnDlzMG3aNMTExKBUqVKoUqWKpMdTUn3zzTcoUaIEFixYgMTERNSvXx+7d+9WmcnXzMwMH374Ib777js8fPgQ7969g6urKyZOnIgJEyZoXEcf6tWrh759+2LLli3w9/dH+/bt0adPH/j5+cHLy0tvcRHJhBBC30EQEVH+8ff3x+nTp3Hv3j3IZDIYGxvr7Ni3bt1Cjx49sHHjRpw/fx716tXD1atXcezYMZw4cULtXEIIpKenY+vWrRg8eDAuX76Mhg0b6iweokxswSEiKgYePXqEEiVKoGbNmvjzzz91dtyaNWvixx9/xNixYxEeHg4jIyN88sknOHjwYJaJ1M8//4xPPvlEZ+cnyg5bcIiIDNzDhw+Vo7bMzc1Rs2ZNvcXy+vVr3Lt3T/m+Ro0aOlsrjOi/mOAQERGRweEwcSIiIjI4THCIiIjI4DDBISIiIoPDBIeIiIgMDhMcIiIiMjhMcIiIiMjgMMEhIiIig8MEh4iIiAwOExwiIiIyOExwiIiIyOD8HyJJwa0UXnVEAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "## Now plot\n", "plt.hist(Mout_error,bins=30,density=True,histtype='step',lw=2,label='Owen & Jackson (2012) efficiency')\n", "plt.hist(Mout_error_1,bins=30,density=True,histtype='step',lw=2,label='Constant Efficiency')\n", "plt.hist(Mout_error_2,bins=30,density=True,histtype='step',lw=2,label='Escape velocity scaling')\n", "plt.ylim((0.,1.))\n", "plt.ylabel('Probability Density',fontsize=12)\n", "plt.xlabel(r'Minimum Core Mass [M$_\\oplus$]',fontsize=12)\n", "plt.legend(loc=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see for planets like Kepler-36c, which contain a large amount of H/He today (relative to other sub-neptunes ~10% compared to 1-2%) the choice of efficiency parameter matters. At young ages these inflated planets with massive H/He envelopes are large. At this point advective energy losses become important the efficiency of the outflow drops with falling escape velocity. " ] }, { "cell_type": "code", "execution_count": null, "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": 4 }