age_PDFS.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "age_PDFS.ipynb",
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Rogerod23/Rogelio/blob/main/age_PDFS.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 311
        },
        "id": "Db88UTqmZ5DL",
        "outputId": "d1feb14f-329d-440b-ad69-d323b21090a1"
      },
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "from matplotlib import rcParams\n",
        "### THIS NOTEBOOK ONLY PLOTS CLUSTERS 2732 and 816.\n",
        "\n",
        "#############################\n",
        "#### LOAD DATA   ############\n",
        "#############################\n",
        "PDF = pd.read_csv(\"https://raw.githubusercontent.com/Rogerod23/Rogelio/main/PDF.csv\")\n",
        "percentiles = pd.read_csv(\"https://raw.githubusercontent.com/Rogerod23/Rogelio/main/percentiles_all.csv\")\n",
        "\n",
        "del PDF[\"Unnamed: 0\"]\n",
        "del percentiles[\"Unnamed: 0\"]\n",
        "\n",
        "\n",
        "cluster_list = [\"11-E\", \"90-E\", \"93-E\", \"383-W\", \"417-W\", \"451-W\", \"531-E\", \"534-E\",\n",
        "                \"589-E\", \"816-E\", \"894-E\", \"1252-W\", \"1381-W\", \"1564-W\", \"1576-W\",\n",
        "                \"1949-W\", \"2449-W\", \"2732-W\", \"2740-W\"]\n",
        "\n",
        "PDF_columns = PDF.columns\n",
        "percentiles_columns = percentiles.columns\n",
        "\n",
        "# Getting name of the age_pdf of cb2018\n",
        "a = np.where(PDF_columns.str.contains(\"tpdf_cb2018\") == True)\n",
        "age_CB2018 = PDF_columns[a]\n",
        "\n",
        "\n",
        "a = np.where(percentiles_columns.str.contains(\"cb2018\") == True)\n",
        "percentiles_CB2018 = percentiles_columns[a]\n",
        "\n",
        "a = np.where(percentiles_CB2018.str.contains(\"t_\") == True)\n",
        "percentiles_CB2018 = percentiles_CB2018[a]\n",
        "\n",
        "\n",
        "\n",
        "# Getting name of the age_pdf of Krumholz2015\n",
        "a = np.where(PDF_columns.str.contains(\"tpdf_modp020\") == True)\n",
        "age_Krumholz2015 = PDF_columns[a]\n",
        "\n",
        "a = np.where(percentiles_columns.str.contains(\"modp020\") == True)\n",
        "percentiles_Krumholz2015 = percentiles_columns[a]\n",
        "\n",
        "a = np.where(percentiles_Krumholz2015.str.contains(\"t_\") == True)\n",
        "percentiles_Krumholz2015 = percentiles_Krumholz2015[a]\n",
        "\n",
        "\n",
        "\n",
        "#PDF COLUMN INFO TO SEE WHAT PLOT YOU WANT TO SEE\n",
        "#print(percentiles_columns)\n",
        "#print(PDF_columns)\n",
        "\n",
        "#########################################################\n",
        "# Colors \n",
        "#########################################################\n",
        "\n",
        "alpha=1\n",
        "blue=(7/255,80/255,155/255,alpha)\n",
        "lblue=(106/255,175/255,210/255,alpha)\n",
        "pistache=(190/255,250/255,215/255,alpha)\n",
        "yellow=(248/255,250/255,3/255,alpha)\n",
        "red=(251/255,105/255,73/255,alpha)\n",
        "wine=(160/255,55/255,98/255,alpha)\n",
        "dgreen=(5/255,149/255,74/255,alpha)\n",
        "gray=(203/255,202/255,201/255,alpha)\n",
        "black = (0,0,0,1)\n",
        "\n",
        "mycolor=[blue,red,dgreen]\n",
        "\n",
        "#########################################################\n",
        "# Font parameters \n",
        "#########################################################\n",
        "fontsize=12\n",
        "font = {'family' : 'serif',\n",
        "        'weight' : 'regular',\n",
        "        'size'   : fontsize}\n",
        "plt.rc('font', **font)\n",
        "\n",
        "mycolor=[blue,red,dgreen]\n",
        "\n",
        "mark = \"None\"\n",
        "lw_ = 1.5\n",
        "linestyle_ = \"-\"\n",
        "\n",
        "#Figure size and subplots\n",
        "i = 0\n",
        "for CB2018, K2015 in zip(age_CB2018, age_Krumholz2015):\n",
        "#get cluster ID\n",
        "  if len(CB2018) == 25:\n",
        "    cid = CB2018[22:25]\n",
        "  else:\n",
        "    cid = CB2018[22:26]\n",
        "#Get percentiles column name\n",
        "  percentile_16th_CB2018 = \"t_16th_cb2018_z014_u3_{}\".format(cid)\n",
        "  percentile_50th_CB2018 = \"t_50th_cb2018_z014_u3_{}\".format(cid)\n",
        "  precentile_84th_CB2018 = \"t_84th_cb2018_z014_u3_{}\".format(cid)\n",
        "\n",
        "\n",
        "  percentile_16th_KRUM2015 = \"t_16th_modp020_chabrier_MW_phi0.50_{}\".format(cid)\n",
        "  percentile_50th_KRUM2015 = \"t_50th_modp020_chabrier_MW_phi0.50_{}\".format(cid)\n",
        "  precentile_84th_KRUM2015 = \"t_84th_modp020_chabrier_MW_phi0.50_{}\".format(cid)\n",
        "  fig, axes = plt.subplots(1,2,figsize=(10,4))\n",
        "\n",
        "  plt.subplots_adjust(hspace=0.5)\n",
        "  plt.subplots_adjust(wspace=0.25)\n",
        "  code=['GALAXEV','SLUG']\n",
        "#figure title\n",
        "  for c in cluster_list:\n",
        "    if cid == c[0:2]:\n",
        "      fig.suptitle('Cluster {}'.format(c), size = 12)\n",
        "    elif cid == c[0:3]:\n",
        "      fig.suptitle('Cluster {}'.format(c), size = 12)\n",
        "    elif cid == c[0:4]:\n",
        "      fig.suptitle('Cluster {}'.format(c), size = 12)\n",
        "\n",
        "\n",
        "  for icode in range(2):\n",
        "      axes[icode].set_ylabel('')\n",
        "      axes[icode].minorticks_on() #...add minor ticks to x and y axes\n",
        "      axes[icode].tick_params(which='minor', length = 2, width = 1, #...format minor ticks\n",
        "                                      direction = 'in', top = True, right = True)\n",
        "      axes[icode].tick_params(which='major', length = 4, width = 1, #...format minor ticks  \n",
        "                                        direction = 'in', top = True, right = True)\n",
        "  #    if icode==0 and icluster == 0:\n",
        "  #        axes[icluster][icode].annotate('Cluster 816-E', \n",
        "  #                                        xy = (0.4,0.95), xytext = (0.4,0.95), \n",
        "  #                                        xycoords='figure fraction', size = 12)\n",
        "  #    if icode==0 and icluster == 1:\n",
        "  #        axes[icluster][icode].annotate('Cluster 2732-W', \n",
        "  #                                        xy = (0.4,0.47), xytext = (0.4,0.47), \n",
        "  #                                        xycoords='figure fraction', size = 12)            \n",
        "      axes[icode].set_xlabel('Log(T) (yr)')\n",
        "      \n",
        "  axes[0].set_title(code[0],fontsize=12)\n",
        "  axes[1].set_title(code[1],fontsize=12)\n",
        "  axes[1].set_xlim(0, 8)\n",
        "\n",
        "\n",
        "  #Extinction plots\n",
        "\n",
        "\n",
        "\n",
        "  axes[0].plot(PDF.loc[:, PDF_columns[2]], PDF.loc[:, CB2018]/PDF.loc[:, CB2018].max(), color = black, marker = mark, \n",
        "          markevery= 10, linestyle = linestyle_, lw=lw_)\n",
        "\n",
        "  axes[1].plot(PDF.loc[:, PDF_columns[6]], PDF.loc[:, K2015]/PDF.loc[:, K2015].max(), color = black, marker = mark, \n",
        "          markevery= 10, linestyle = linestyle_, lw=lw_)\n",
        "  \n",
        "\n",
        "  axes[0].axvline(percentiles.loc[0, percentile_16th_CB2018], \n",
        "                   color =mycolor[0], linestyle = '-', linewidth = lw_,label='16th percentile')\n",
        "  axes[0].axvline(percentiles.loc[0, percentile_50th_CB2018], color =mycolor[1],\n",
        "                   linestyle =  '-', linewidth = lw_,label='50th percentile')\n",
        "  axes[0].axvline(percentiles.loc[0, precentile_84th_CB2018], color =mycolor[2],\n",
        "                   linestyle =  '-', linewidth = lw_,label='84th percentile')\n",
        "  axes[0].legend(loc = \"upper right\")\n",
        "\n",
        "  axes[1].axvline(percentiles.loc[0, percentile_16th_KRUM2015], \n",
        "                   color =mycolor[0], linestyle = '-', linewidth = lw_,label='16th percentile')\n",
        "  axes[1].axvline(percentiles.loc[0, percentile_50th_KRUM2015], color =mycolor[1],\n",
        "                   linestyle =  '-', linewidth = lw_,label='50th percentile')\n",
        "  axes[1].axvline(percentiles.loc[0, precentile_84th_KRUM2015], color =mycolor[2],\n",
        "                   linestyle =  '-', linewidth = lw_,label='84th percentile')\n",
        "\n",
        "\n",
        "  if len(CB2018) == 25:\n",
        "\n",
        "    plt.savefig(\"age_pdf_{}.png\".format(CB2018[20:25]), bbox_inches = 'tight',dpi=100)\n",
        "  else:\n",
        "    plt.savefig(\"age_pdf_{}.png\".format(CB2018[20:26]), bbox_inches = 'tight',dpi=100)\n",
        "  if i == 0:\n",
        "    plt.show()\n",
        "  plt.close()\n",
        "  i += 1\n",
        "\n"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAEmCAYAAACgb95DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXxU1fn48c+ThCRAQhIgrCaCVhahQlEEQoAEhaqAohaFFje02urXKrYqtn4Fl35F/ZrqV6sWtUZRxKWg4sJPRbFCRRGhuIBUJQSyQEIMZCHbzPn9cWeGSZgkM8lsyTzv12teMHfOvfcM4PHJc859jhhjUEoppZRS7RcV6g4opZRSSnUWGlgppZRSSvmJBlZKKaWUUn6igZVSSimllJ9oYKWUUkop5ScaWCmllFJK+YkGVkopRGSKiLwnIutF5GPH7y8TkVgRuVpE8kQkN9T99ERExorIahF5V0S2ishLItLT7fPhIvKhiGxwfH6r22f9RORvIrJORP7l+O6jPdwjU0T+IyIVIrLS7fg9IvKK2/ubRKRYRNYE8jsrpcKXBlZKRTgROQt4DvitMSbLGDMJyAGeAk42xiwDcv10rzwRyfLHtdycB7xljJkOnAakAve4ff4s8C9jTCZwJnCr4zsDjAS6GWPOMMZkAB8AK5rewBizAbgWMMB8t48mO17OdjnAh8aYWf76ckqpjkUDK6UimIhEAY8BfzbGfOc8box5B3g1ZB3zzbPACwDGGBuwBRjk9vkI4BPH5weB/wA/c3y2FVjk1vYz4Phm7rMRiAdOBRCRrkAZ0EtETnYc6w5UtvcLKaU6Lg2slIpsY4DBwLsePlsAfNn0oIjcLyLlInK54/3fRKTGmYkSkb4i8o7b9NutjuPPAP2AhxxTjs4A5VIR2SQiH4nIChHp4Tj+mOM+d4vIKhHJF5ElTftjjPmPMeaI45z+wAzgabcmbwGzHJ+fgBVofeo496AxpsDxWVfgUuBJT39QxphqYDMwxXFoPPD/gC/cjmVgBWBKqQilgZVSke1Ex68FTT8wxhxxZICaHr8F2Ob2/hqg2K3J74H1xphs4Oc4ghpjzBWOdjc6phy3iMhErGnHWcaYKY5+5DjaX+u4zzhgDjAN2NfcF3GsddoFvGiM+YfbR1cCQ0Tke6xs1k3GmA+anHsLsBfoAvyxuXsAHwFZjt9PAdY7jk1xO/ZRC+crpTo5DayUUv5WBpwtIiOMMVXA9BbaXg6sMcaUON6vAH4lIuLW5i1jjM0Y860x5qnmLmSMmQP0ByaIyP+6fbQK+MQYcyJwCnCTiIxpcu79WGuzvgQ+EJHoZm6zHpjo+HyYMWYHjQOrQcaY3S18X6VUJ6eBlVKRzbmuaqAfr/kA8A/gJRHZhjU115zjgKmOqcH1wCPAfqCXW5tD3t7YGFMJ3AvcKCLxIjIca8H6Q47P9wLvATd7ONcAd2Gtv8oWkdHOfonIQ45mG4FuwASg1nFsA9DH8TRhtbd9VUp1ThpYKRXZtmIt5j4mqyQiTzfN7LipA+Lc3ie7/b6PMeYRY8xI4A/AchE5Ec/2Am86pgazHE/unWaMKfX2C4jI1U0OVQHRWAvNYx3H6t0+rwec67hmisgA5weOqc9aoLsxZptbv250fF4FfA7cytF1WuXAvx3HNnnbb6VU56SBlVIRzBhjB64GbhORnziPi8gVwE9xW0vVxG6sUgWIyBSsLI7TvW61oD7FCsKcU3sVQDcRyRaRG7DKOMwQkRTHtYYCvtaAmudYq4VjCvEqYIsj4NmJtW5rnuPzROBcrLIKYJVnuNzte1/s6OunLdzvI2AmjddSfQRchK6vUiriiZX9VkpFMhGZBNwJxGD9wPUtcKsxptSREfojVgboEWPMn0VkGPAK1nqqNcB/AeVYC8X7YWVvGoAk4FljzP857vNfWPWgDgNXGmO+FpH5jvOrsYKw3xljdonI/VhBXzGwzFEjylPfZwI3Od4mYgVSNxljfnB8fhrWgngBEoD3gduMMQ0iMhJYDPR1fO8o4Pami9ub3O8s4DljTB+3Y7MdfzZprfxRK6U6OQ2slFJKKaX8RKcClVJKKaX8RAMrpZRSSik/0cAqQjk2lV3reJT8n44NaG9xLiJ2a/eeiDze5Fiq4zwjIv8Wkakt3GeiiNgcFbGdxwY6Km0bEdkh1gbAsY5rlovIcx7usd7tVS7W5rlbHJ9vEpEkt+vfLyI/ishz/vwzU0p1HiJyvIi87hj/1jnGluscn/3UbYz6VBxbFjk+6+dp/BORnm7j0wtu7XuJyF9FZKPj889EZFkLT9yqDk7XWEUgEfk5sAz4uTFmp+PYT7GeaLrUGPOm41gq8A1gBwY0rcItIgaYZox5v4V7PQxcAiw2xjzS5LOPgVpjzJmO9/OxNsRd1tI9RGS9MSZLRLoBJViLoFe6fd4VeMkYc66vfzZKqcggIh8CrxhjHnO8z8J6AOGnjveDsJ5+Pcl9H0238z2Of87xyfH7BKz9J59we4CjG/AiUG+M+UVAvpwKKc1YRRixNt19HGvT3Z3O48aYL7GKM7rX+7kQa3uSJI5u4+HrvU7ECuIu8tDkeqxCjBeKtT/chTSzT1sTtzj6XI21D9zFTT4/B3jb1/4qpSLK6ViV9AEwxqzHsZm3H90MlDuDKsd9qrHGsCo/30uFCQ2sIk+zm+4aYxYbY/6f26EzgZXAOxwbvHgjE2vgehFrG5BG1b2NMduwAqkHgaXAXaaVFKrjp8HP3A69jLV9Sg+3YxdiVf5WSqnm7AFuEZHuzgPGmKV+vseFeB5rvzXGXObne6kwoYFV5Gl20113ItIXOGCMqcMKjM4XkRgf73UhsNIY82+sQo1zPLT5E1YV7H7GmK3NXOdBObrlSVNvYdVLOs/R7+5Y04klHtoqpZTT9VgbhBeIyDOOQrf+diKtjLWq89HAKnK5NrkVkVMdgctXcnTz2guxskFgFYCMB87w+uLWNGAfY8w+x6EVeM56VQOlwFkicnwzl/u9c2uRph8YY44Ab7pde6bjvVJKNcsYsw5Ix9p2aRDwoYj8zR+XbvJeGr0RecWxMD7PD/dSYUgDq8jzg+NX17ScMWaLI2jZBvR2HD4XuNORJXoHKwDytE6qOZOBU90yTbOAcSKS3qTdbY5XPtaUYIs8BVdYAeB0xxON5wOrfOinUipCGWOqjDFPGWOygWzgKhE5wcvT6zi6FyUAIhLH0c25wRpvmy6BmAPcDjT3g6Tq4DSwijxfAHnAWc01EJF+wHfGmClumaLfYU0Hdmnp4iLyS8dvfwFMdtvEdhzWNikXubX9CXCCMeYfwHXAhSLSalbMQ8r+bazB7FIg1hhT1to1lFKRrWkZGWPMR8BBrId1WjrPOcbtxVqv6m6w47jTP2hhrFWdkwZWEcZRMuG3wJ/E2icNsGpLAUOwSiv8Amvtkrt3gO7AtFZucbWIRAPpxpjiJp+9RuOs1585+oTfOqy1XP/nxVquO5t8pxqs6cq78X0DX6VUZDpTRE53vnH8wGbHWg/akqsdv67GynD1dJwfhfUU9etubR8AEkXk9273iQEmceyUoeokfF2MrDoBY8xaEZkH/EVE4rHWT0VjTakVAfcBP4jIu261q/6MtUj8MRG5A2uqEOAOEfmN2+VTsTa5HSoiVxpjngZrHZfjnJNFZDXWprc/AW7EejLnp8AI4GTgMxFZiLW41NM9Rnj4Wi9hBYSr2/anopSKMPcDD4iIHWv8swOzjTFHxNpk3Lnp9yMi4qk0wmKsB28+E5EDjmOrnXUAAYwxFWJtcH63iHyKVc4mGWvmYFRAvpUKOS0QqpRSSinlJzoVqJRSSinlJxpYKaWUUkr5iQZWSimllFJ+ooGVUkoppZSfaGCllFJKKeUnISu30Lt3bzNo0CCv2hYWFjJgwACvr+1r+0i+R7j2K5D32FVwCIAEqerQ3yMQ9zD7rR2IpO9xrmO7yguoq61lZF9vC1LDli1bSo0xqT51qhPo1q2bOfnkk9t0blv+Hv11vt5b7x3o870Zdz2NP67zy60tFxOqJWTfe8uWLdXGmO6ttQtZuYXTTjvNfP755161FRF86aev7SP5HuHar0DeY9rtVpmZ9/88q0N/j0Dco+EvtwIQs/A+17EzX1vE+vXraXhoky/X32KMOc2nTnUCImLaOqa25e/RX+frvfXegT7fm3HX0/jjdOZriwBYd/59of7e0lq7TjkVuHjx4oCf05Z7+CoY36Ot5wT6Hvo9AsvXe3ibXVZt196/9/acH4x/c4G6t37v4N+7PSLhe3fKjFW40u8RHtqTsQpH/vz7aC5jBfD+7KW+9EkzVhGko48JbaXf23vOcfe9e2Y228abjJUv45C/daqMVSija3/S7xFe9Hsof+vfv3+ouxASkfpvUL93xCnyppHXgZWIjBWR70Tk8lba/UpEtojI5yLyoIi0Gt21ZsmSJe29RFjQ7xFe9Hsof49r7VkY25FF6r9B/d4Rp9CbRl49FSgi5wNzgEOttBsJPAiMBMqAdcC1wF+9uY8KL/X19ezbt4+amppQd8Wv/jTD+p/fjh07QtyT8GMyzgNA3P5sbk+z9tt2//OKjo4mOTmZ3r17ExXVIRLfx9BxTSkVCN6WW9hsjFktIutbaXcl8LYxphRARP4O3IIOQB3Svn37SExMZNCgQfgh8Rg2ogvKARgyMDnEPQk/nh53jvrROjY0xTpmjKG+vp79+/ezb98+0tPTg99R/9BxTSnld179qGmM2efl9cYCO93efwOMEJGuvnZMhV5NTQ29evXqVEGVaj8RITY2loEDB1JVVRXq7rSZjmtKqUDwdw6/L43T6uWAAL39fB8VJBpUqfr6OioqKmhoaGh0vKNOAbaBjmtKBcne/HwmTZoU6m60SyBGRk/PYB7zf+fCwkJExPWK4MVwSoU151PVLQVSS5YsafTfM9DZVnHruKZUEBQWFrJhwwby8/ND3ZU2j2v+DqwOAO4LV5KxBqSSpg0HDBiAMcb10gFI+Wrz5s385Cc/ITc395jPDh06xFVXXUVmZiajRo3irLPOOibj4rRs2TIGDRrE5ZdfHtgOh4G8vDycW6589tlnjB492ovinwaR1gMr9/+e8fLpmQ5CxzWlgqCyosL1sNTGjRtD3Ju2j2v+Dqw2A0Pd3p8MfG2MOeLn+6gIt3r1av7yl7+QlJTk8fNf/epXnH/++WzYsIGtW7cSHx+P3W4H4JEHl7Jo4bWutldffXWnDaoGDRrE+vXrG73/17/+BcDpp5/OQw891Oo1rPEkoqeEdVxTKghKSkoQEbp16xYWgVVbtSuwEpHeIvKxiPR0HHoKOEdEeolIFHA58EQ7+6jUMcaOHcuKFStITEw85rMtW7ZQUlLCjBkzACvT8tprrxEbGxvsboal5GQfn4Y0VsYqUui4plRoHCwrIyk5mYyMDDZs2BDq7rSZV4GViJzqeCR5NLBIRFY5PuoKDAO6ARhjvgL+ALwLfApsAx7zc5+V4rjjjt393GndunUMGjSI66+/noyMDGbNmsVXX30FwEsvvcTqV1bw8fp1ZGVl8ec//9l1Xk1NDddccw1jxozhnHPO8Vi/q66ujqysLESEe++9l+nTpzNixAiWLj26zUJ9fT0333wzGRkZZGZmcuedd2KMaXTuX//6V2bMmEFiYiLr16+nsrKSX//612RmZpKZmcmVV15JSYk10/T9998zffp0pkyZwqRJk1wZpzfeeINhw4YxZcoUbr31VsaPH8/EiRM5cOAAAFdccQXFxcXceOONZGVlsWXLFs444wxEhLy8PI9/dpVVVSy48fdkZmaSkZHBE088gcHQGTNWOq4pFV5qamro3q0bmZmZbN++nUOHWiwxF7a8qmNljNkCZHk4vhdIbXLsBeAFf3ROhZc/PP0J/959MKD3GDW4F/975YR2XSMvL49Vq1bxyiuv8Mgjj/D888+TnZ3NDz/8wMUXX8yGT7dSsC+fVS+vaHTexx9/zFdffUVSUhKjRo1i9erVzJs3r1Gb2NhY1q9fj4hQXl7Ou+++S1lZGSNGjGDMmDFMnz6d+++/ny1btvDxxx9jt9uZMmUKJ554IvPnz3ede/DgQd566y1efPFFevTowU033YTNZmPDhg3Y7XZmz57N119/TWZmJjNnzuTmm29mwYIFbN++nalTp7J7927OPfdcysrKuO6668jNzeW+++7jnHPO4emnn+a2227jmWee4cMPP+Shhx4iKysLsILOlp70XHjHndhtdjZs2EBFRQWjRo3ijn5JjBo7pl1/J+FIxzWlwkd1dTW2hgbi4uIYM2YMxhh27tzJuHHjQt01n0XM89IqctTW1pKens7s2bMBmD9/PrGxsaxZs6bF88aNG0dKSgpRUVGMHDmS3bt3t9h+7ty5APTs2ZNzzjmHlStXApCbm8tll11GdHQ0Xbp0Yc6cOSxfvrzRueedZ1U4nzdvHqNHj+a5555zrfOKioriwQcf5OSTT+bTTz/l+++/55JLLgHglFNOYeDAgbz55puuaw0dOpTBgwe7Pm+t382x2+08/+oqrph3EQCJiYnMmjWLN199nc6YsVJKhY+iImsbvtjYWNfWUM5jHY23ldeVancmKVhSUlLo27dvo2MDBw5k376W60H26NHD9fu4uDjq6upavY9Tr169+PLLLwGrYn1OTg7PPPMMAJWVlcesa3JfdF9SUkJtbS2pqUeTJCeddBIAH374ISLCtGnTXJ/V1tY2SpG79zs+Pr7VfjfH2Y9b7/4fuj5sFRUvLy8nfcgJEbXGSikVfIWF1gN3sbGxrs3MNbBSKkyMHj2a1157rdGxkpISv2+QW1ZW5ipVUFpa6hoM0tLSuP3225kzZw5gZYLKy8ubvU5qaipxcXGUlJQwfPhwwBpkoqKiSEtLo0uXLo2e7KuqqgpIcU5nPx75n7s5/ayZgLVe7JPv/o1mrJRSgVRUVEQvIDYujpQ+fRCRDhtY6VSg6nTOO+88KioqXMHIhx9+SEVFBTNnWsFC94QEjhw5gjGG888/v833efXVVwE4ePAgb7/9tmtq8PLLL2fFihXYbDYAnn322UaL5JuKiori0ksvddXjstvtXHnllRQXFzNu3DjS09NZtcpaV93Q0MDs2bPZtWuXV31MTEykurqaDz/8kIcffrjFtlFRUVwy50KWv7rKdeyee+7hzVdf04yVUiqgnEFUXGwsMTEx9OnTRwMrXzkrFGsBPdUWW7ZsISsri23btrF06VIuuOAC12eJiYmsWbOGRYsWMXHiRJYsWcI777zjmo47Y/o5fLdrJxkZGWRkZLBixQpyc3NZu3Ytjz/+OMuWLWPt2rXk5uayYsWK5rpA3759+fnPf86kSZO44YYbmD59OgA333wzQ4cOZeLEiWRnZ/PBBx9w9913A7jazJ07lw8++MB1rZycHKKiolxPBc6cOZPRo0cTHR3NmjVrWLZsGVOmTCE7O5t58+YxatQoPvjgA5YuXcq2bdtYvHgxb7zxhut75OTkAHDNNdfwhz/8gdtuu40zzzyTM844w3X/jz76iBtvvJHi4mJXdi3nzjuoPnKEjIwMpkyZwqFDh5hz6Ty8zVg5KxXT+Sqve0XHNaXapqioCImKIqZLFwD69+8fNoGVr+OaGONpp4bAO+2008znn38eknsr7+zYscM1NdWZ7CqwpuWGDPSxnpMbEWH37t1eVC3vWMx+ax2a9D1azmLL3p106RLDKf1+ckz75v6NiMgWY8xpgetpeNJxTSnPpt1uPXDz3j0zPX5+6aWXcnXtHsaPH0/MQusJ5/3797NlyxYAznxtEQDvz17q8fxg8HZc06lApVQrDKJrrJRSAVRUVEScWxHn/v37U1xcHMIetZ0GVkr5wFnkE6zptIKCgtB2KMCsPbLQtetKqYAqKipqtDtGv3792L9/v2utakeiTwUq5QNngdBI4dxfUTNWSqlAKioqIjZuiOt9//79sdlslJaWHlM+J9xpxkop1SxnYKVxlVIqUGpqaigrKztmKhA6Zi0rDayUUs3SjJVSKtCc+5t20cBKKdXZHc1YaWCllAoMZwHlLjFHVydpYKWU6pSOZqyUUiownIFVjFtg5VxXtX///pD0qT00sFJKNUszVkqpQPvxxx+BxoFV165diYuLc33WkWjlddUhXX755WRlZTV6VVZWuj7ftm0bEyZMYOLEicyaNYuDBw82e6277rqLfv36RcS/xQ0bNriqv7/xxhsMGzbMVT7CE+ejzt6GVVp5Xcc1pXzlKWMlIvTs2ZOysrJQdcvF13EtZIHVgAEDMMboAKTabP369Y1eCQkJgFVr6rzzzuN//ud/2LhxI2PGjOE3v/mN67xFC6/lkQePVu+94447OOuss4Le/2AQEfLy8lzvJ06cyCuvvALAueeey6JFi1o839eM1ZIlS3Ds5lDYhu52eDquKeU7T4EVEFaBlS/jmtaxUp3OO++8Q3R0NNnZ2QBcddVVDBo0iJKSElJTU0Pcu9ASEZKSkrxur2uslFKBFu6Bla90jZXqsH7zm98wadIkzj77bD788EPX8c2bNzNs2DDX+7S0NLp168YXX3zBww8/zMfr17H6lRVkZWXx9NNPu9qVlpZyySWXMHLkSC677DKP98zPz2f8+PGICA8//DBnnHEGI0eO5Nlnn3W1qaysZMGCBWRmZpKRkcETTzxxzLm5ublMnz6duLg48vLyOHDgAHPmzGHy5MlMmDCBhQsXcuTIEcDacHry5MlMmTKFM844g507dwKwbNkyBg0axNy5c7nmmmsYM2YM55xzDjU1NQCcffbZgFUhPisriz179rju35wDJaVcuODXTJ48mYyMDN5809rfS9dYKaUCpby8nMTExGPGmY4aWGnGSnnN9srfMPt+COg95LgTiJ5zTavthg8fztSpUxk7diybN28mOzubDRs2MHr0aPbv339MViY5OZkDBw5www038NHGTxl4XDqP5DTezHPjxo1s2rQJYwwDBw7kk08+YcKECY3apKens3LlSgYPHkxsbCzr1q1j9+7djBgxgnHjxjFs2DAWLlyIzWZjw4YNVFRUMGrUKEaOHElmZqbr3KioKN59911ycnKIi4tj/vz5rmm6uro6xo8fz/79+0lJSeGss87ipZdeYurUqbz11lucd9557Nixg6uvvprCwkKefPJJvvrqK5KSkjjllFNYvXo18+bN45133kFEWLlypWuzaOf9mzP/ut8x/tQx3J3zEEVFRYwYMYLcNS8z9OShXv4NKqWUb8rLy0lOTj7meM+ePemIm5prxkp1SLfeeitjx44FYOzYscyYMYO//e1vrs89ZWUcc+TNys7OJi4ujvj4eIYMGcLu3btbbD937lwABg8ezIQJE3j55Zex2+0sX76cBQsWAJCYmMisWbNYvnx5o3PPO+88AG666SbsdjvvvfceV1xxBWBtm/P3v/+d1NRU3nzzTRISEpg6dSoAM2bMoLi4mE8//dR1rXHjxpGSkkJUVBQjR45std/NKSgo4P1/fsyVv7S+V//+/Rk7diz/74230MlApVSgtBRYacZKdWreZJJCJT09na+//hqAPn36sHXr1kafl5eX06dPnxav0aNHD9fv4+LiqKura7F9SkqK6/e9evWiqKiIkpISamtrueWWW+jatavr3qNHj250rntGbd++fQCN1n852+/bt4+ysrJGT+6lpqY2esrRvd/x8fGt9rs5zn5cdv2NSFw8YBXnO25o8xkupZRqr5YCqyNHjriWRXQUGlipDun+++/nlltucb3fv38/AwZYT8KOHTuWlStXuj7bu3cv1dXVnHrqqX7tQ1lZGT179gSs9VkjR44kNTWVuLg4Hn30UVdGrb6+nurq6mavk5aWBkBJSQnp6ekA/PDDD6SkpJCWlsZxxx3XaOPnw4cPEx8f79fv4t6PV576G31GjAJg165d7LdX+P1eSinlVF5e7hp/3DnH145Wy0qnAlWHlJOT49pfavfu3bz++utccsklgLVou6GhgY8++giAv//971xwwQWujFD37gkcOVJNVVUVv/rVr9rch1dffRWwgqBNmzZx0UUXERUVxaWXXtpo6u+ee+7hueeea/Y6AwYMYNq0aeTm5gJQW1vLRRddRF1dHTNnzqS0tJTNmzcDUFVVRXZ2NocOHfKqjwkJCVRXV/P888+7+ttiP6ZMYvmr/3Ad++Mf/8iWTz7z6l5KKdUWLWWsgA43HagFQlWH9Ic//IHzzz+fKVOmMGfOHB599FGmTJkCWNN4r732GosWLSIzM5PPP//c9WQewMzzf8EH760lOzubGTNmkJOTw9q1a8nNzWX16tXcddddbNu2jaVLl/LBBx802wcR4cwzz2TWrFk8/vjjricRc3JyqK6uJiMjgylTpnDo0CGuvfZaysrKXOuysrKy+Oabb1zXev755/nyyy+ZPHkyU6dOZdGiRfTt25cePXrw9ttvc/PNNzNlyhTOPvts7rzzTlJTU1mxYgW5ubmsXbuWxx9/nGXLlrm+x4oVKwC47rrrmDt3Lk888QTjxo1rdP/XX3+dpUuXsm3bNq6//noAlj/6f3yyeQuTJk1i0qRJHH/88UyeNtXrvxctEKrjmlK+CvfAytdxTVpb0Bsop512mumIq/0jyY4dOxg+fHiou+F3uwqsmilDBh77H7I38vLyGDx4cKuL4Tsis99aZyV9jwPgP//5D5Xd7HTr1p2hKccd0765fyMissUYc1pgext+dFxTyrNpt1ulW967Z2aj43a7nZiYGG6//XbuSKkFIGbhfQB88cUXnHrqqaxevZpH2QTA+7MbP80dTN6OazoVqJRqlt1uR/SJQKVUgFRUVGCMCeuMla80sFLKB/n5+Y2m05x76XVWdrtdKy0opQLGWXW9MwVW+lSgUj5IT09n06ZNoe5G0FgZKx0mlFKB0SiwavJMTmJiItHR0R0usPIqYyUi8SKSKyKbRORzEZneTLtEEVkuIptF5DPHOQn+7bJSKlg6c8ZKxzWlQq+ljJWIdMgiod5OBS7BWug+HvglsFJE+npo99/A8cB4x+t44HY/9FMpFQKdfI3VEnRcUyqkWgqsoGNWX281sBKRKOAq4GkAY8wuYCsw30PzEZhiML0AACAASURBVMBnxhibMcYOfAb8zH/dVUoFU2fNWOm4plR4aC2wSklJ6XyBFXAC0AvY6XbsG8DTI4dvA1NFpLuIdAOmAp96aKeU6gDsdrvHfRc7AR3XlAoDrQVWycnJXhdEDhfeBFbO1Lj7NysHjtl4zRjzV2A98AOQB3wB3NWuHiqlQsJutztqdXXKwErHNaXCgDOwct8/1V1ycrKrTUfhS7mFptUQjxltReRPwCisNQjpWD8VXu7pYs4Kxc6XVipWvtixYwfZ2dlkZmbys5/9jPvuu89ju8WLFzNo0KAWr3XXXXfRr1+/iPg3uGHDBqZPt9Zov/HGGwwbNqzRBs/u7HY7AN4krJyViZ0vOk7ldR3XlAqh8vJyevToQXR0tMfPQ5mxauu45k1gdcDxq3ueLtntuLvfAX8zxtQYY2qAJ4C7PV10wIABGGNcLx2AlC8uu+wyMjIy2LBhA++//z733Xcfa9eubdRm//79PPXUU8ecu2jhtTzy4NHqvXfccQdnnXVWwPscCiJCXl6e6/3EiRN55ZVXADj33HNZtGhRs+c6AytvMlZLlixp9N8zUNj2XgeFjmtKhYHmtrNxSkpKcmSsgr/TRVvHNW8Cq++BMmCo27GTgc0e2sYC9W7v64FEbzqilC++/vprJkyYAECvXr046aST2Lp1a6M2d999N7/97W9D0b2wJSLNptyb8iVj1QHpuKZUGPjxxx9bDKySk5Opr6/HZrM32ybctBpYOZ6CeRJYACAiJwGjgRdEZLiIrBMRZw7vfeBiccB6hPnDwHRdRbIZM2awZs0aAH744Qe+/vprxo0b5/r8u+++o7CwkMzMzEbnPfzww3y8fh2rX1lBVlYWTz/9tOuz0tJSLrnkEkaOHMlll13m8b75+fmMHz8eEeHhhx/mjDPOYOTIkTz77LOuNpWVlSxYsIDMzEwyMjJcG0C7n5ubm8v06dOJi4sjLy+PAwcOMGfOHCZPnsyECRNYuHAhR44cAWDLli1MnjyZKVOmcMYZZ7Bzp7XeetmyZQwaNIi5c+dyzTXXMGbMGM455xxqamoAOPvsswGYO3cuWVlZ7Nmzx3X/5hwoKeXCBb9m8uTJZGVlsX79ejrjGisd15QKD61lrJyfNTQ0BKtL7eZtSeUlwBMisslxzjxjTLGIDAKGAV0AG3Ad8AiwCWs0zgeu9m+XVaj8fsMy/l36Q0DvMar3CTyY2fo/maeffppzzz2XE088kbKyMnJycpg6darr8zvuuIM777yTgwcPNjrvhhtu4KONnzLwuHQeyWm8mefGjRvZtGkTxhgGDhzIJ5984sqKOaWnp7Ny5UoGDx5MbGws69atY/fu3YwYMYJx48YxbNgwFi5ciM1mY8OGDVRUVDBq1ChGjhxJZmam69yoqCjeffddcnJyiIuLY/78+a5purq6OsaPH8/+/ftJSUnhrLPO4qWXXmLq1Km89dZbnHfeeezYsYOrr76awsJCnnzySb766iuSkpI45ZRTWL16NfPmzeOdd95BRFi5cqVrnZnz/s2Zf93vGH/qGO7OeYjvvvuO0047jVFTx5M+eFCrfycd0BJ0XFMqpMrLyzn++OOb/dw9sIqLiwtWt9rFq8XrjrUFlxtjxhtjTjPGvOs4vskYM9Cx7gBjzAFjzMXGmHHGmNONMb8wxhQF8guoyHTBBRcwYcIEvv/+e7Zv305OTg5ffPEFAJ9//jldunThpz/9qU/XzM7OJi4ujvj4eIYMGcLu3btbbO/cM3Dw4MFMmDCBl19+GbvdzvLly1mwYAFgbckwa9Ysli9f3ujc8847D4CbbroJu93Oe++9xxVXXAFAbGwsf//730lNTeXNN98kISHBFTTOmDGD4uJiPv306NP+48aNIyUlhaioKEaOHNlqv5tTUFDA+//8mCt/aX2vvn37MmrUKN5e9UabrhfudFxTKvS8zVjZbJ0vY6WUV5mkYNixYwfvv/8+L7zwAgBpaWlMmzaNBx54gBdffJHFixfz17/+1efr9ujRw/X7uLg46urqWmyfkpLi+n2vXr0oKiqipKSE2tpabrnlFrp27QpYA8fo0aMbneu+zmnfvn0ApKamuo452+/bt4+ysrJGT+6lpqY2ysS59zs+Pr7VfjfH2Y/Lrr8RiYvHZrNRWFjIsOpT2nQ9pZRqTSRPBSoVNpyBQ5cuXVzHunTpwuHDh6moqOC7777j8ssvB6z/aIuLi8nKyuLiiy/262L2srIy1+7rpaWljBw5ktTUVOLi4nj00UcZO3YsAPX19VRXVzd7nbS0NABKSkpIT08HrHVjKSkppKWlcdxxxznWOlkOHz5MfHy8375H03688tTf6DNiFIcOHeKrr74ibmBKK2cqpZTv7HY7hw8fbvWpQOhYgZUvdayUCgvDhg1j4MCBvPjiiwBUVFTwxhtvMHXqVBITE/n2229Zv34969ev56GHHqJfv36sX7/eFVR1757AkSPVVFVV8atf/arN/Xj11VcBKwjatGkTF110EVFRUVx66aWNpv7uuecennvuuWavM2DAAKZNm0Zubi4AtbW1XHTRRdTV1TFz5kxKS0vZvNl6WK2qqors7Gyv67okJCRQXV3N888/7+pvi/2YMonlr/4DAGMM9957L59t/MSreymllC8OHz6MMUYzVkqFWlxcHK+99ho33XQTL774IpWVlZx//vnccMMNjdrdddddrFq1ypWxuv/++zn99NOZef4vuO2m/yI7exM33ngjOTk5rF27lvj4eEaNGsWXX37Jtm3bWLp0Kenp6Y0WxbsTEc4880yKiop4/PHHGTZsGAA5OTnceOONZGRk0KVLF372s59xxx13UFZW5lqXlZWVxWOPPcbJJ58MwPPPP891113H5MmTsdlsLFq0iL59reLgb7/9Nr///e9dtVTuvPNOUlNTWbFiBbm5udTU1PD4448THR3t+h5Dhgzhl7/8Jddddx1z586lR48evPjii8yZM8d1/4ULF7J06VKKi4u5/vrreeSRR1j+6P9x3aI/MWnSJBoaGhgzZgzZZ00LyN+jUiqytbadjftnHSmwalT8Kpiv/v37G8AsXrzYqPD0zTffhLoLAfHtvh/Nt/t+bPP5u3fvNtZ/Op2PvXivsRfvNcYYU1paajZv3my+Kd1jdpbt9di+6b+RxYsXG6xKfoUmRGNLKF86rinl2Zl/WmPO/NOaRse2bt1qALNq1SpjjDH1ObeY+pxbGrWx2+0mNjbWDH5wjjlj9a1B6687X8e1kE0FOisUa2VipcKTsSoN+1Qg1FmpmPCvvB4QOq4p5T1vMlYiQnJyckgzVr6Oa7rGSikf5OfnN5rOs9lsIe5R4DgDq85YIFQpFXreBFZgLWDvSFOBusZKKR+kp6ezadOmUHcjKI4GVkop5X/eBlbJycmUd6DASjNWSimPXFOBIe6HUqpz8iWw6kgZKw2sVIs0axG5nJswNxdZuT5XSqk2cAZW7kWOPdHASnUa8fHxHDx4UIOrCNXcGitjDHV1dRQUFNC9e/fgd0wp1SmUl5fTo0cPoqOjW2zX0QIrXWOlmnXcccexb98+SkpKQt0Vv9pffgQA2+GuIe5J+DGHfwRAyiooLy/n0KFDVIpV6d7ercLVLiYmhqSkJHr37h2SfiqlOr7WtrNx0sBKdRpdunRh8ODBoe6G3/3u9jcBeO+emSHuSfhp+MutAMQsvI/bbruNnJwcJr20EID3Zy8NZdeUUp2Mt4FVUlIS9hJ7h1l+oFOBSimP6urqGu3HqJRS/uRLxgo6TvX1kAVWhYWFiIgW0lMqTNXX1xMbG+vTOUuWLEGsiqIDAtKpMKfjmlLeKy8vJyWl9U3eQx1Y+TquaeV1pZRHdXV1bQqstPK6jmtKeaOjZKy08rpSyi90KlApFUgdJbDylQZWSimP2jIVqJRS3rDb7Rw+fNjrxeuggZVSqoPTjJVSKlAOHz6MMUYzVkqpyKEZK6VUoPz4o1Uzz6fAyqaBlVKqA2vL4nWllPKGt/sEAnTv3h0RwaYZK6VUR1ZfX69TgUqpgPAlsBIRYmJidCpQKdWxacZKKRUovgRWgAZWSqmOTxevK6UCpW2BlS2QXfIbrbyulPJIK6/7Tsc1pbzTkTJWWnldKeUXbclYaeV1HdeU8kZ5eTkiQo8ePbxqHxMTrZXXlVIdm5ZbUEoFSnl5OT169CAqyrswRNdYKaU6PF28rpQKFG+3s3HSwEop1eHp4nWlVKD4GlhFx8Rgs9mor68PYK/8w6vASkTiRSRXRDaJyOciMr2Fthkisk5E/iki34jI9f7rrlIqWDr7VKCOa0qFTlsyVgCHDh0KVJf8JsbLdksAMcaMF5EhwCYRGW6M2e/eSEQGAw8B5xhjSkVkBLDArz1WSgVFBGSslqDjmlIhUV5ezuDBg71uHxMd4zqvd+/egeqWX7SasRKRKOAq4GkAY8wuYCsw30Pzm4BnjDGljrZfG2N+77/uKqWCpTNnrHRcUyq0ysvLSUpK8rp9R8pYeTMVeALQC9jpduwb4DQPbc8AYkXkbRHZKCL/KyLxfuinUirIOvnidR3XlAqhgwcP+pR5cgZWzvpX4cybwKqv41f3MLEc6OOh7SDgN8ClQBZwMlYKXSnVwXTyvQJ1XFMqRI4cOUJVVVVEB1ZOpsl78dAmDlhhjCk1xtRjDT6XO9LujTgrFDtfWlBPqfBhs9mw2WxeZ6yclYmdLzpO5XUd15QKspKSEgBSU1O9PicUgVVbxzVvAqsDjl/dl+8nux139yPgvvBzH9agdExY6qxQ7HzpAKRU+HA+0uxtxspZmdj5Ivwrr+u4plSIlJaWAh0jsGrLuOZNYPU9UAYMdTt2MrDZQ9ttNE6lpwJ1wEFvOqOUCg/OwKoTr7HScU2pEHFmrHyZCoyOjkakkyxeN8bYgSdxPF4sIicBo4EXRGS4o7ZLtKP5k8A8EenmeL8AeN4Y0zG2pFZKAdbCdei8gZWOa0qFTlumAkWE6OiYDrHGypc6Vk+IyCbHOfOMMcUiMggYBnQBbMaYl0TkBKx6MBXAt8CNfu+1UiqgfJ0K7KCWoOOaUkHXlqlAsKYDyw90ksDKGFMDXO7h+CZgYJNj9wL3+qNzSqnQ6OwZK9BxTalQKSkpITo62qc6VuAIrDpAxkr3ClRKHcMZWHXyjJVSKgRKSkro3bs3UVG+hSAaWCmlOqwIWLyulAqR0tJSn6cBwQqsOsXidaVU5ImEqUClVGg4M1a+0oyVUqrDipDF60qpEGhPxkoDqxY4KxRrAT2lwk9bM1bOSsV0nMrrfqXjmlKta0/G6vDhw9hswa104uu4FrLAylmhWAcgpcJPWzNWzkrFhH/l9YDQcU2pltlsNsrKytqYsbJKyx0+fNjf3WqRr+OaTgWqiFVUVOTKzKjGdI2VUioQysrKMMa0eSoQwr/6ugZWKiLZ7XaGDh3K448/HuquhCUNrJRSgXDggLUdZ1unAiG4+wW2hQZWKiJVVFRQUVHBvn37Qt2VsKSL15VSgeAcc4877jifz9XASqkwVlFRAUBlZWWIexKeNGOllAqE/Px8ANLS0nw+VwMrpcJYxWErsHIGWKoxzVgppQJh7969iAgDBvj+4HB0tAZWSoWtwxXWUyUaWHmmGSulVCDs3buX/v37t+mHNs1YKRWm6uvrqampAXQqsDmasVJKBcLevXvbNA0I+lSgUmHLmaVKSEjQjFUzNGOllAqE9gRWIkJiYqJmrJqjFYpVqDiDqczMTA2smqGV19tGxzWlWrZ3717S09PbfH5ycnLQAyutvK5UKw4frqBbt+70799fpwKboZXX20bHNaWa11DfwJEjR9qcsYLQBVZaeV2pFtTV1dK1azyJiYmasWqGTgUqpfytptZa29qewCopKUmnApUKNzabjeioaFdg5fhJRLmpr69HRIiOjg51V5RSnURtbS3QvsAqOTlZF68rFW5sNjtR0VEkJCRgt9tdTwiqo+rq6jRbpZTyq5oa/wRWmrFSKszY7Taio62MFWgtK080sFJK+duR6moSEhLo169fm6+hgZVSYchmsxEVpYFVS+rr67WGlVLKr6qqqxk+fLjzCbs2cU4F2u12P/bMvzSwUhGloaEBYwzRjqlA0CKhnmjGSinlb9XVVQwfPrxd10hOTsZut4f1uK2BlYooVVVVADoV2ArNWCml/MnWYKOurq7dgVVSUhIQ3tXXtUCoiijV1dUAOhXYirZmrLRAqI5rSnniHHv9kbGC4O4XqAVClWrB0YyVTgW2pL6+vs2BlRYI1XFNqaaqqq2xt6MGVlogVKlm6FSgd+rq6nQqUCnlN9XV1YgIJ5xwQruuE4rAylcaWKmIcnQqMEoDqxbo4nWllD9VV1fTtWs3YmJi2nUdDayUCjPuGSudCmyeLl5XSvlTZWUVCQnd232dTrN4XUTiRSRXRDaJyOciMr2V9l1E5D8issQvvVTKT9wDq9jYWGJjYzVj5UEkZKx0XFMqOA4dOkRtbQ3duye0+1rOwCqcM1be5uSWAGKMGS8iQ4BNIjLcGLO/mfZXA3380UGl/Ml9KhDQjZibUV9fT/fu7f/pMswtQcc1pQLuyy+/BPBLxio2NpZu3bqFdWDVasZKRKKAq4CnAYwxu4CtwPxm2icAFwFv+K+bSvmHe8YKrMBKpwKP1dkXr+u4plTwbN++HcBvP6yF+7Y23kwFngD0Ana6HfsGOK2Z9r8HHgFs7euaUv7nCqyirMAqISFBM1YetLXcQgei45pSQbJ9+3ZiYmKIi4vzy/U6Q2DV1/Gr+0qxcjykxEUkFZhsjHnVD31Tyu9cU4HROhXYks6esULHNaWCZvv27X5ZX+XUGQIrJ9PkvaddFP8b+LM3F3NWKHa+tKCeCgZnxsp9jZVOBR7L18XrzsrEzhcdp/K6jmtKBZAxhi+//JIEP67ZTEpKCspTgW0d17wJrA44fk12O5bsdhwAETkBGGyM+cCbGzsrFDtfOgCpYKiurnatrwKdCmyOr+UWnJWJnS/Cv/K6jmtKBUFhYSGVlZV0697Nb9cMVsaqreOaN08Ffg+UAUM5OuicDLzdpN0UYICIrHe8HwbUiEgWcIUxZrc3HVIqkKqqqoiK6up6r1OBnkVAuQUd15QKgm+//RaArl27ttLSex1+KtAYYweeBBYAiMhJwGjgBREZLiLrRCTaGPOMMeZUY0yWMSYLWAvkOt7r4KPCQlVVFdHRR//Z61SgZ5198bqOa0oFhzOw6tbN/xkrRxYp7PhSx+oJEdnkOGeeMaZYRAZh/QTXBcfTMiISC7zL0Z/sTjHGXODnfivVJjoV6J0IWLwOOq4pFXDffvst3bp189sTgWAFVg0NDVRXV4dlvT2vAitjTA1wuYfjm4CBTY7VAVl+6JtSfldVVUWUW2CVmJhIfX09tbW1fv0Pv6OLgKlAHdeUCoJdu3YxZMgQv17TfVubcAysdK9AFVE8TQWC7hfYlO4VqJTyh2+//ZahQ4f69ZrhvhGzBlYqoniaCgR0OtCNMQa73a4ZPKVUu9TW1pKXl6eBlVKdmfVUYOOpQNDAyp3dbgf8+xSPUiryfPfdd9jtdg2slOrMdCqwdXabtWuLBlZKqfbYvdt6cPaEE07w63U1sGqGs0KxFtBTwVRdXe3aJxB0KtCT9mSsnJWK6TiV1/1KxzWljsrPzwfg+OOP9+t1nYFVMKqvg+/jWsgCK2eFYh2AVDA1fSrQ+USJc6sbdTSwio+P9/lcZ6Viwr/yekDouKbUUXv27CE2Npa+ffu23tgHzqcCg5Wx8nVc06lAFTEaGhqoq6trNBWogdWxbLrGSinlB/n5+aSlpbn2ZvWX+Ph44uLidCpQqVCrrq4GaDQVqIHVsXTxulLKH/bs2UN6enpArh3O29poYKUihjN4itKMVYt08bpSyh/y8/P9vr7KSQMrpcKAK2Ola6xapBkrpVR71dXVUVhYqBkrpTozZ/DkPhUYExNDbGysBlZu2rN4XSmlAAoKCjDGBCxjlZSUFLSnAn2lgZWKGJ6mAsHKWmlgdZQuXldKtZez1IJmrJTqxDxNBYIGVk3pVKBSqr327NkDaGClVKfmaSoQoFu3bq6gS+nidaVU+zkzVmlpaQG5vgZWHmiFYhVsOhXonfYWCNXK6zquKVVYWEivXr0C9gNacnIytbW1HDlyJCDXd6eV15Vqhk4Fekcrr7edjmtKWQoLCxkwIHA/X/Xu3RuAgwcPBuweTlp5XalmNDcVqIFVY3a7nbi4OL9XS1ZKRY6CgoKgBFalpaUBu0db6cipIoZOBXrHZrfr+iqlVLsEK2NVUlISsHu0lQZWKmJUV1cjIsdkYjSwasxut2kNK6VUm9lsNoqLixk4cGDA7qEZK6XCQFVVFd26dTvmuAZWjdltmrFSSrXdgQMHsNvtOhWoVGdXWVlJYmLiMcc1sGrMrlOBSql2KCgoAAhoYNWzZ09ERAMrpUKpoqKChISEY453796d6upq19NwkU4DK6VUexQWWg/PBTKwio6OpmfPnhpYKRVKFRUVzWasgKDUQ+kIbHa7rrFSSrWZM7AK5BorsKYDNbBSKoRamgoEdDrQQTNWSqn2KCwsJCoqij59+gT0Pr1799anAt1phWIVbK1lrDSwsthttjYHVlp5Xcc1pQoKCujbty8xMTEBvU+wMlZaeV2pZrS0xgo0sHJqT8ZKK6/ruKZUoGtYOQUzsNLK60p5oBkr79h1jZVSqh2CHVg5gp6woYGVihjNBVbO2lYaWFm08rpSqj0KCwsDvnAdIDU1lfr6eioqKgJ+L19oYKUigt1up6qqSjNWXtDF60qptqqtraW0tDRoGSsIvyKhXgVWIhIvIrkisklEPheR6c20myUi74rIByLyhYj8zr/dVaptnEGTrrFqXXsWr3ckOq4p5X9FRUVAYGtYOYXrfoHeLtlfAogxZryIDAE2ichwY8z+Ju3+AswxxmwVkX7A1yKyxxjzuh/7rJTPnKnixMREONT4Mw2s3BiDMSYiAit0XFPK74JRHNQpXAOrVjNWIhIFXAU8DWCM2QVsBeZ7aP5XY8xWR7ti4EPA40+BSgVTo8CqCQ2sjnJWn+/si9d1XFMqMIJVHBSgX79+ABQXFwf8Xr7wZirwBKAXsNPt2DfAaU0bGmP+0uRQPBBeoaSKSJWVlYAGVq2xOQKrCMhY6bimVAAEM2PVkQOrvo5f3SdQyoEWS6qKSA9gLPBM27qmlP+0lLGKj49HRDSwwlpfBRERWOm4plQAFBQU0KVLF3r16hXwe8XFxZGSkuJa1xUufHkqsGmhCGml/X3AXcaYPZ4+dFYodr60oJ4KJGdg5WnxuojQvXt3Daw4OhXoa2DlrEzsfNFxKq/ruKaUHzlrWDnGgYDr379/wDJWbR3XvAmsDjh+TXY7lux2/BgicjVQb4z5a3NtnBWKnS8dgFQgtZSxAjSwcmjrGitnZWLni/CvvK7jmlIBEKzioE79+vULWMaqreOaN4HV90AZMNTt2MnAZk+NReR8IBu4wfH+JG86olQgtbTGCjSwcmprxqoD0nFNqQAIVnFQp0BmrNqq1cDKGGMHngQWgGtAGQ28ICLDRWSdiEQ7PpsM/A74L6C7iCQAtweq80p5y5uMVXV1dTC7FJYiZfG6jmtKBUaoMlbhtK2Nt2uslgAiIpuAF4F5jseOk4BhQBdHuxeBLKAUqHC8jvdjf5VqE2dg5XwCsCnNWFkiaPE66LimlF9VVFRw+PDhoAZW/fv3p6amhsOHDwftnq3xqkCoMaYGuNzD8U3AQLf3wcv/KeWDiooKunXrRnR0tMfPNbCyREodK9BxTSl/27PHeqbj+OOD93OHs+RCUVERSUlJQbtvS3SvQBURKisrm50GBGsjZg2sImqNlVLKz/Lz84HgBlb9+/cHwquWlQZWKiJUVFS0GFhpxsqigZVSqq2cGav09PSg3dM9YxUuNLBSEaGiosJjDSsnDawskbJ4XSnlf/n5+cTExLiCnWDQjJVSIaIZK+84F69HwhorpZR/5efnk5aW1uxa1kBITk4mLi7OtZVOOAhZYOWsUKwF9FQweBtYhdMju6HQ3qlAZ6ViOk7ldb/ScU1Fsj179gR1GhCsnTPS0tJc67sCwddxLWSBlbNCsQ5AKhhaW7zevXt3bDYbdXV1QexV+LHb7UhUVJt/4nRWKib8K68HhI5rKpLl5+cHdeG606BBg1zruwLB13FNpwJVRGgtY5WSkgJAWVlZsLoUlux2O1FROiwopXxTX19PQUFB0DNWYD2FmJeXF/T7NkdHUBURWlu83qdPHwBKSkqC1aWwZLPbidbASinlo8LCQux2e8gCq/3791NTUxP0e3uiI6jq9Ox2e6tTgc7A6sCBZvfgjQh2u00zVkopn4Wi1ILToEGDAAK6zsoXOoKqTs+5B6AGVq2z23QqUCnlu//85z8AnHRS8Pcnd67rCpfpQB1BVafX2gbMAKmpqYAGVna7naggPiqtlOocdu7cSVxcXMgWrwMBXcDuCw2sVKfnDKxaWmOVkpJCdHR0xK+x0sXrSqm22LlzJ0OGDAlqDSunAQMGEB0drRkrpYLFm4xVVFQUqampEZ+xsmlgpZRqgx07djBs2LCQ3DsmJoa0tDTNWCkVLJWVlUDLgRVY66wiPbCqr6sjtkuXUHdDKdWB1NTUsHv37pAFVhBeJRe08rrq9LzJWIEGVsYYauvqiI2La/M1tPK6jmsq8nz33XfY7XaGDx8esj785Cc/4dtvvw3ItbXyulJNOIt+Jicnt9guNTU1otdYlZeXY7fZiIuNbfM1tPK6jmsq8uzcuRMgpBmrkSNHUlpaGpAfjrXyulJNOOfd09LSWmwX6Rkr5yamce3IWCmlIs+OHTsAGDJkSMj6MGLECAC+/vrrkPXBSQMr1ent2bOHfv36ER8f32K7Pn36UFFRwZEjR4LUs/BSUFAA0K6pQKVU5Nm8eTNDhgyhqVdOugAAFUhJREFUe/fuIeuDBlZKAV988YXrf+aBlJeX51VtlUjf1sb5d9GeqUClVGQxxvCvf/2LiRMnhrQf/fv3Jzk5WQMrFblKS0uZMGECQ4cO5YEHHsButwfsXnv27HEVkGuJs0hopAZWzqlAzVgppby1a9cuDh48GPLASkQYMWKEBlYqcq1cuZK6ujrGjh3LLbfcwsqVKwNyH7vdTn5+vk8Zq0hdZ1VQUEBMly5ax0op5bWNGzcCkJGREeKe4AqsHAvNQ0ZHUBUSzz33HKNHj2bdunX079+fVatWBeQ+xcXF1NXVeZWxivTAqrCwUKcBlVI++de//kXPnj0ZOnRoqLvCT3/6U8rKyti7d29I+6GBlQq6HTt2sHnzZi699FKioqI499xzWbt2LTU1NX6/l7NgnGasWldQUKBPBCqlvGaMYd26dUycODEsMt2ZmZkA/POf/wxpP0L/J6EizvLly4mOjmbevHkAzJ49m6qqKtatW+f3ezlLLXiTsUpISCA+Pj6i11jFxmnGSinlnU8//ZS8vDwuvPDCUHcFgFNOOYWUlBTWr18f0n5o5XUVVMYYVq5cybRp0+jXrx8A2dnZJCYm8vrrr/v9fr5krEQkYvcLbGhooLi4mLjY9mWstPK6jmsqcqxcuZK4uDhmz54d6q4A1p6vkydP9ntgpZXXVVj797//ze7du/nFL37hOhYXF8fZZ5/N66+/js1m8+v99uzZQ+/evb2urxKpRUL379+P3W5v9xOBWnldxzUVGRoaGnj55Zc555xzSEpKCnV3XKZMmcL333/v13VWWnldhbVVq1a51lW5u+CCCzhw4AAfffSRX++Xl5fn1TSgU58+fVxlByKJq+q6Ll5XSnnhmWeeoaioiAULFoS6K41kZ2cD8M4774SsDxpYqaBatWoVkyZNctWMcjr33HPp0aMHzz77rF/vt2fPHq+mAZ3GjRvH9u3b2b9/v1/7Ee5cxUF18bpSqhXV1dUsXryYjIwMZsyYEeruNDJq1ChGjBjBU089FbI+aGClgmbXrl18/fXXXHDBBcd81rVrVy666CL+8Y9/UFlZ6Zf7GWO8Lg7qNHv2bIwxrFmzxi996Ci0OKhSyhvGGK699lqKioq47777nGuPwoaIcPXVV7N582a2bt0akj54FViJSLyI5IrIJhH5XESmt9D29yKyxfG62X9dVR3dc889B9DsQsfLLruMqqoqv9W0ysvL48iRIz4FVqeccgqDBw/mtdde80sfOordu3cTHR1NbJcuoe5K0Oi4ppRvbDYbN998M88++yx33nmnq7xBuLnkkkvo2rUr9913X0ju723GagkgxpjxwC+BlSLSt2kjETkL+DUw0fG6UkTalSe84YYbWLp0Kdu3bw95NVXVdjt37uSBBx7g4osvJj093WObiRMncuKJJ/L444/T0NDQ7ns+9NBDxMTEMGvWLK/PERFmz57Ne++9R0VFRbv70BFUVlaSm5vLtGnTIMx++gywJYRoXFOqo9m8eTPZ2dk8+OCD/Pa3v+W///u/Q92lZqWkpLBo0SJeeuklXn311aDfv9XASkSigKuApwGMMbuArcB8D82vAVYYY2qMMTXAC8Bv2tq5hoYGPvnkE2677TZGjRpFdnY2Gzdu1ACrg7HZbPz617+me/fuPPzww822ExH+9Kc/sWnTJq655pp2/T0fOHCAJ598kvnz5/u0xgqsjFpdXR1vv/12m+/fkTz22GOUlpayePHiUHclaEI5rikV7mw2GwUFBbz55pv88Y9/5NRTT+X0009nx44dPPPMMzz22GNhNwXY1G233cbYsWO57LLLWLFiRVDjhhgv2pwA9AJ2uh37BjjNQ9uxwItN2l3b5s7FxPDZZ59RVFTEyy+/zL333ktmZiZDhw7l7LPPZtCgQSQnJxMVFdXoJSJh/5fe2RljqKqqYu/evTz77LP88MMPPPPMM/Tte0xCoJErrriC3bt3c/fdd5Ofn8/s2bNdFdF9sWbNGmpqali0aJHP506cOJGBAweyYMECvvzyS0aNGuXzNToKYwwPPPAAP//5zxk/fjwNn6wOdZeCJWTjmlLt8eOPP/LOO+9gt9ux2WzYbDavft/0fV1dHVVVVVRXV7t+PXjwIPv27aO4uNhV+iYmJoZT/3979x8sVXnfcfz94QJXuVwEiVFA5A5iECStUTQaS8IMEDGYpLaYFNEZrtIwSbD+SkYwmXpvLEmDIDCZpDRqg1dpy5gyhIHRREHEWqnij8TGYlMtChaDqIAoVoVv/zjPknXdvfecZc+evXu/r5md3bP7nPN8n3v2PufZ5zznOWefzdKlS2ltbWXAgAEZ/wXi6dOnD2vXrmX69OnMnDmTW265hUmTJnHKKacwYMAA+vfvf2S2+Fx7Ib/dUPhekjZFnIZV7ki4L++9vcDYEmkL0xU9KuYm0su5+eabS879MmTIEK655hpmz57NPffcw6pVq1i+fHkqt0BxlTdhwgQWLlxYdNB6Me3t7TQ2NnLnnXcyd+7csvOdMWNGWfevamhoYPPmzcyfP58FCxaUnX93UYkJLdva2mhvb89/q9YnCM28XnOuHDt27GDmzJlHtY2Ghgb69u1LU1MT/fr1O/I8aNAgJk+ezLBhwxg2bBjjxo1j/Pjx9OvXr0LRV9dJJ53Exo0bWblyJR0dHXR0dBztEI9Y9VqchlVOYT9aqeZbrP62oUOHJp4vqKmpiTlz5hw5TbRnzx4OHDjA4cOHP/Jw2WtqamLw4ME0NzcnWi93SvCmm27i5ZdfLvsqwVGjRpW1HsDIkSNZtWoVixYtYv/+/WVvpztobm4uOe4trra2tg81ICR1l8nAMq/XnEti9OjRbNu2jYaGBnr16kVDQ0Oi17VwT79q6tu3L62trbS2tmJmvPPOO+zfv58DBw5gZkdOEeafKix8L/c8bty4WP/ccRpWuWmoB5Z4XZh2YN7yQCCVG6/lbj9SOB+Sqx+SEo+PqrThw4dnmr9LTU3Wa851pbGxsayeeBcdU5qammLfiaNccZquLwBvAPl7cizwRJG0T8RM55xzWfJ6zTmXii4bVmZ2GLgduBJA0mnAmcBKSWMkbZDUEJIvB2aE+WGOIbqEeXk6oTvnXHm8XnPOpSXuGKs2YLmkLWGdGWb2qqQW4HSgD3DIzO6XdAbwaFjvH8xsfWVDds65imjD6zXnXIXFGsUW5m+ZZWbnmdl4M/tVeH+LmQ0Lc7vk0i42s7PD49ZKBFkvV9V4OWqLl6NnS6Ne66kD13vqd9DL3ePEuiqwW1weUHAZd7fl5agtXg5Xabt27co6hEz01O+gl7vHGRInUbdoWCVVTms66TrVaLFXoxzlrpN2Hl6OdCXNY/v27anE4f6gEnOJZZX30fByV399L3e6lNXtYcaPH29bt26NlVZSounok6bvyXnUalxp5jHlu+sAeHDBF7t1OdLI44MlNwLQ+7o/3Lx08pp5bNq0iQ+Wbkmy/SfNrNgs5nVNkpVbp5azHyu1vufteae9fpx6t1j9kzN5TXQXjQ2X/DDrcnc5BXtmDStJrwEvxUw+FEgyeCFp+p6cRznreB6eR1dGmFmPm2RO0tvAf5a5ejn7sVLre96et+fdtTFm1uUkWJk1rJxzzjnn6k1djrFyzjnnnMuCN6ycc8455yqkZhpWkmZLelTSv0r6jaTPlUg3U9KTkrZKWqz8W8nXgDjlkPQnktaF2Z2flLRAUs3sC4i/P/LSb5K0okrhxZbgezVG0npJj0h6TtL3qx1rV2J+t3pLWibpKUn/JmmNpJOyiLcnCbOyr5C0JdRNn886prRJ6iPp2vC//7CkxyRNyjquapF0mqT3JU3MOpZqSnpsqAfh+PBQKPPTkm7sLH3cmddTJelSYBLwWTM7JKkV+MjBQNI4YDEwjug+XxuAbwA/rmK4JcUtB/A3wG1mtlZSf2Ar0Y1el1Uv2tISlCOXfhrwKWB7dSKMJ8H36jhgJTDdzF6UdAKwpLrRdi7BPvkaMA34pJkdlNQBLAX+onrR9khtRGNWz5P0CWCLpDFm9vuM40rTMOAa4Ewz2ydpCvALSaPN7JWMY6uG7wHvZR1ENSU9NtSRu4AHzOw7kgYDv5P0azO7v2hqM8v8ATxLNNq+q3RLiG4nkVu+Ang26/jLKMe1hAsHwvIiYH3W8SctR0jbC3go7JsVWcde5v64Afhh1vFWqCw/Bn6et/x14Pms46/nR/gf2EN0sMm9twG4IevYUi73YOCKgvf2AJdmHVsVyn4O0b0mtwMTs46niuWOfWyopwfwNnBx3vK/A/NLpc/89JOkjwNjgDNDV9sjkuaUSH4OsC1v+TngDEnHph1nV5KUw8yWWtg7wTHAa9WIsysJ9wdEjdv1wL6qBBhTwnJMAg5KWh26uG+XNKh60XYuYVnuAz4t6WOSegMXEVUCLj0jiRoZhXVTXc/jZWavm9ndueUwLKMvNVKXpaw9PHqMMo4N9WQ98EUASSOBM+ikXq2FU4EtgIBLgMnAx4HHJe0zs38uSHsiHz6A7w3rfgzYkX6onWohfjmOkNQAXAhcXo0gY2ghZjkkNQJ/GdLNq3KcXWkh/v5oAc4Gzgf+B/h74J+AqdUKtgstxCyLma0Lp6K2AQeBp4l6rVx6TgzPhXXT2AxiydLniOYm3Jx1IGmSdBHwWzPbqdoa4pu2Fso4xtWJq4C1kl4AjgeuN7ONpRJn3mMFNBLF8SMzO2Rmu4C7gStLpC828VYtfLuTliPnW8AaM6uVXoUk5ZhLdGr23SKfZS1JORqBdWb2YuhJXARcKCnWDTerIHZZJF0OtAKnASOIDvC11uitV4V1Uy3US1Uh6Rjg+8AsMzucdTxpCRcZ3Qj8IOtYMlDuMa4erAYeM7NTgT8Crpd0VqnEtdCwejM85w/y3AmcXCTtbmBg3vJAosqsFrqek5QDAElfIDpdUEsHvljlCAO+/5RoUF8tSrI/3iySjhJps5CkLFcDd5vZm+EAtwyYJ6nL2YJd2XaH58K6aXeRtHUnnAL8KbDEzJ7MOp6UXQb80szeyDqQDCQ+xtUDSWOIeuiWApjZDuAB4Nul1qmFhtXvgHeIuhVzTqD4tPNPAKPzlscSdckeTC+82JKUA0nnEV1Rc4VFV1ecln6IscQtx7nAAGCDpE3ALGBquPT6/CrE2ZUk++OZIukokTYLScrSF3g/b/l9olP+/VKLzr1AdJVyYd30RDbhVN1i4HEzu1dSo6RTsg4oRROAi0M9t4noirilkn6RbVhVkegYV0f6hufCenVAqRUyb1iZ2f8BHUTnMAm/rL8K3BUG4D4i6fiQ/A7gC5IGhy7ZWcDyDML+iCTlCC3gZcBsoHeYcuGvs4n8w+KWw8weMLM/NrOJZjYRWAHcH5Yfyyr+nITfqzuJKstchXEVsNHMdhZuNwsJy/Ig8OeS+oTly4H/MLNa6NWtS6Fn8HbCKZHwI+lMoik86lqYz6c3sCLUY6dSx6eGzGyOmV2QV++9ClxrZl/OOLTUdVYPZRlXFWwDXgFmAEhqBr4ElBxjVRP3Cgw7aDnRr7wPgH8BbiXqYnwK+FTuICdpJnA9cJhokOS3rBYKQfxySHoU+EzB6i+ZWUsVwy0pyf4I6VcDZxFd3bgN+LyZZT6/S8LvVSvRNBj7iU4tf91qaA6iBN+tJuA2osH47wEHgKvN7PlMAu8hwhij5cDpRA2Nm8zsV9lGla5wkUSx71W7mbVVOZyqknQusBA4j6jOW21m38s2qvSVqodq5RicFknjiepVAf2JfsDON7MPiqav87+Hc84551zVZH4q0DnnnHOuXnjDyjnnnHOuQrxh5ZxzzjlXId6wcs4555yrEG9YOeecc85ViDesnHPOOecqxBtWdU7SFEnPSDJJD+dNJFnJPIZKapc0UdL23KzEIc9fh9fPSGqT1BpnZnZJvSTdKOnUMmO6PkzE6pyrM16vuVrm81j1AJImAg8BfUpNaHYU2z4GeJjoHlrDgYm5yQElGTDFzB4MMUwE/pZocrUrzey/Otnud4FDZlbWzU4lDQZ+CUw1sz3lbMM5V7u8XvN6rVb1zjoA1+1dDTxlZi9I2kvp+0Y9C/yvmb0r6UfAEmBasYSShgJ/BZR9zzEze13SvUA78M1yt+Oc65G8XnNl81OBDknnhO70zeH5nLzPTg33ontM0kpJq0O3+NdCkq8Q3VoIM3u91K+1gs8eBqZ00n1/CfB0qKyOlfSspIOS7ggxfVvSHkmLJf1E0l5Jt4TYXpbUlpfPdEkNR/UHcs51O16vuax4j1UPJ+k44D5gupltkjQBuE/SKDPbC/wjsMbMfiBpGPAcsMTMfipJRDebfTFJnmb2qqT3gE8SVRKFzspt08wOSpoE7ADawufLgPPN7IZQhrHAp4GLgFFEd6AH+G+iO7GfDLyUJEbnXPfl9ZrLkvdYuYuB/Wa2CcDMHgHeBL4kaQRwLnBP+OwVPlxhDCZqnL9VRr5vASeW+OzE/G2a2W7gAeDy8NY0okoz33ozO2Rmz5vZHXl55LbnnOs5vF5zmfGGlTsZeK3gvdfC+0PCcv4gyTfyXis8l3MFhFH6+6ci2+wArgivvwLcW/D5vhJ50Ek+zrn65PWay4zvmB5K0nGSphF1RZ9Q8PEJwE5gV95yzuC8168Dh4DmMkJoBn5f4rPdRba5Fhgq6ULgcOjOj5MHneTjnKsjXq+5WuANq55rEHApsA5olvRZAEkXhM/WmtlLwOOEX1RhLMIFuQ2Y2WHgt8CIJBmHS4aPDesW85vCbZrZu0S/5n4GrIqZVQtR9//OJPE557otr9dc5nzwep0Lgzbbw+KqMAcLQD9gt5ntlzQVWCypF1E380V5v5wuA+6SdDGwjWiulvzu7J8TVUpHurDDVTGrw+JiSdeZ2ca8dSYAm8MYg2JWA/Mk9TGz9/Pe7wD+jLxxCJIWEg00nSfpeDO7rSCfNQXbcM51c16veb1Wy3yCUNep8E/9Rt7yemCdmf1dWG4mGvg5zcx2ldhM/vZ6Ew3YvM7Mnukk3UJgu5n9JO+9McBcM+ty/hZJ/YFNwJfD4FTnnAO8XnPp8lOBritLJZ0OIGk48BlgQ+5DM3sL+CpwXcztzQIWdlb5BN8BRkgaKWlGqLhmAXfFzOebwDe88nHOFeH1mkuN91i5Tkm6jKhyeRvoD9xqZnHHAlQqhluBqcDjZnZVNfN2ztUfr9dcmrxh5ZxzzjlXIX4q0DnnnHOuQrxh5ZxzzjlXId6wcs4555yrEG9YOeecc85ViDesnHPOOecqxBtWzjnnnHMV8v81oCFp00ewUwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 720x288 with 2 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 315
        },
        "id": "ZzeJmHqOZTCo",
        "outputId": "b9bbfe87-3294-4ec0-bfd3-c01a890df468"
      },
      "source": [
        "#Routine to make extinction + age plots\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "from matplotlib import rcParams\n",
        "### THIS NOTEBOOK ONLY PLOTS CLUSTERS 2732 and 816.\n",
        "\n",
        "#############################\n",
        "#### LOAD DATA   ############\n",
        "#############################\n",
        "PDF = pd.read_csv(\"https://raw.githubusercontent.com/Rogerod23/Rogelio/main/PDF.csv\")\n",
        "percentiles = pd.read_csv(\"https://raw.githubusercontent.com/Rogerod23/Rogelio/main/percentiles_all.csv\")\n",
        "\n",
        "del PDF[\"Unnamed: 0\"]\n",
        "del percentiles[\"Unnamed: 0\"]\n",
        "\n",
        "\n",
        "cluster_list = [\"11-E\", \"90-E\", \"93-E\", \"383-W\", \"417-W\", \"451-W\", \"531-E\", \"534-E\",\n",
        "                \"589-E\", \"816-E\", \"894-E\", \"1252-W\", \"1381-W\", \"1564-W\", \"1576-W\",\n",
        "                \"1949-W\", \"2449-W\", \"2732-W\", \"2740-W\"]\n",
        "\n",
        "PDF_columns = PDF.columns\n",
        "percentiles_columns = percentiles.columns\n",
        "\n",
        "# Getting name of the age_pdf of cb2018\n",
        "a = np.where(PDF_columns.str.contains(\"AVpdf_cb2018\") == True)\n",
        "ext_CB2018 = PDF_columns[a]\n",
        "\n",
        "\n",
        "a = np.where(percentiles_columns.str.contains(\"cb2018\") == True)\n",
        "percentiles_CB2018 = percentiles_columns[a]\n",
        "\n",
        "a = np.where(percentiles_CB2018.str.contains(\"AV_\") == True)\n",
        "percentiles_CB2018 = percentiles_CB2018[a]\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "# Getting name of the age_pdf of Krumholz2015\n",
        "a = np.where(PDF_columns.str.contains(\"AVpdf_modp020\") == True)\n",
        "ext_Krumholz2015 = PDF_columns[a]\n",
        "\n",
        "a = np.where(percentiles_columns.str.contains(\"modp020\") == True)\n",
        "percentiles_Krumholz2015 = percentiles_columns[a]\n",
        "\n",
        "a = np.where(percentiles_Krumholz2015.str.contains(\"AV_\") == True)\n",
        "percentiles_Krumholz2015 = percentiles_Krumholz2015[a]\n",
        "\n",
        "\n",
        "# Getting name of the age_pdf of Krumholz2015\n",
        "a = np.where(PDF_columns.str.contains(\"tpdf_modp020\") == True)\n",
        "age_Krumholz2015 = PDF_columns[a]\n",
        "\n",
        "a = np.where(percentiles_columns.str.contains(\"modp020\") == True)\n",
        "percentiles_Krumholz2015_age = percentiles_columns[a]\n",
        "\n",
        "a = np.where(percentiles_Krumholz2015_age.str.contains(\"t_\") == True)\n",
        "percentiles_Krumholz2015_age = percentiles_Krumholz2015_age[a]\n",
        "\n",
        "\n",
        "#PDF COLUMN INFO TO SEE WHAT PLOT YOU WANT TO SEE\n",
        "#print(percentiles_columns)\n",
        "#print(PDF_columns)\n",
        "\n",
        "#########################################################\n",
        "# Colors \n",
        "#########################################################\n",
        "\n",
        "alpha=1\n",
        "blue=(7/255,80/255,155/255,alpha)\n",
        "lblue=(106/255,175/255,210/255,alpha)\n",
        "pistache=(190/255,250/255,215/255,alpha)\n",
        "yellow=(248/255,250/255,3/255,alpha)\n",
        "red=(251/255,105/255,73/255,alpha)\n",
        "wine=(160/255,55/255,98/255,alpha)\n",
        "dgreen=(5/255,149/255,74/255,alpha)\n",
        "gray=(203/255,202/255,201/255,alpha)\n",
        "black = (0,0,0,1)\n",
        "\n",
        "mycolor=[blue,red,dgreen]\n",
        "\n",
        "#########################################################\n",
        "# Font parameters \n",
        "#########################################################\n",
        "fontsize=12\n",
        "font = {'family' : 'serif',\n",
        "        'weight' : 'regular',\n",
        "        'size'   : fontsize}\n",
        "plt.rc('font', **font)\n",
        "\n",
        "mycolor=[blue,red,dgreen]\n",
        "\n",
        "mark = \"None\"\n",
        "lw_ = 1.5\n",
        "linestyle_ = \"-\"\n",
        "\n",
        "#Figure size and subplots\n",
        "i = 0\n",
        "for CB2018, K2015, K2015_age in zip(ext_CB2018, ext_Krumholz2015, age_Krumholz2015):\n",
        "#get cluster ID\n",
        "  if len(CB2018) == 26:\n",
        "    cid = CB2018[23:26]\n",
        "  else:\n",
        "    cid = CB2018[23:27]\n",
        "#Get percentiles column name\n",
        "  percentile_16th_CB2018 = \"AV_16th_cb2018_z014_u3_{}\".format(cid)\n",
        "  percentile_50th_CB2018 = \"AV_50th_cb2018_z014_u3_{}\".format(cid)\n",
        "  precentile_84th_CB2018 = \"AV_84th_cb2018_z014_u3_{}\".format(cid)\n",
        "\n",
        "\n",
        "  percentile_16th_KRUM2015 = \"AV_16th_modp020_chabrier_MW_phi0.50_{}\".format(cid)\n",
        "  percentile_50th_KRUM2015 = \"AV_50th_modp020_chabrier_MW_phi0.50_{}\".format(cid)\n",
        "  precentile_84th_KRUM2015 = \"AV_84th_modp020_chabrier_MW_phi0.50_{}\".format(cid)\n",
        "\n",
        "  percentile_16th_KRUM2015_age = \"t_16th_modp020_chabrier_MW_phi0.50_{}\".format(cid)\n",
        "  percentile_50th_KRUM2015_age = \"t_50th_modp020_chabrier_MW_phi0.50_{}\".format(cid)\n",
        "  precentile_84th_KRUM2015_age = \"t_84th_modp020_chabrier_MW_phi0.50_{}\".format(cid)\n",
        "  fig, axes = plt.subplots(1,3,figsize=(10,4))\n",
        "\n",
        "  plt.subplots_adjust(hspace=0.5)\n",
        "  plt.subplots_adjust(wspace=0.25)\n",
        "  code=['GALAXEV','SLUG']\n",
        "#figure title\n",
        "  for c in cluster_list:\n",
        "    if cid == c[0:2]:\n",
        "      fig.suptitle('Cluster {}'.format(c), size = 12)\n",
        "    elif cid == c[0:3]:\n",
        "      fig.suptitle('Cluster {}'.format(c),  x = 0.52, size = 12)\n",
        "    elif cid == c[0:4]:\n",
        "      fig.suptitle('Cluster {}'.format(c), size = 12)\n",
        "\n",
        "\n",
        "  for icode in range(3):\n",
        "      axes[icode].set_ylabel('')\n",
        "      axes[icode].minorticks_on() #...add minor ticks to x and y axes\n",
        "      axes[icode].tick_params(which='minor', length = 2, width = 1, #...format minor ticks\n",
        "                                      direction = 'in', top = True, right = True)\n",
        "      axes[icode].tick_params(which='major', length = 4, width = 1, #...format minor ticks  \n",
        "                                        direction = 'in', top = True, right = True)\n",
        "  #    if icode==0 and icluster == 0:\n",
        "  #        axes[icluster][icode].annotate('Cluster 816-E', \n",
        "  #                                        xy = (0.4,0.95), xytext = (0.4,0.95), \n",
        "  #                                        xycoords='figure fraction', size = 12)\n",
        "  #    if icode==0 and icluster == 1:\n",
        "  #        axes[icluster][icode].annotate('Cluster 2732-W', \n",
        "  #                                        xy = (0.4,0.47), xytext = (0.4,0.47), \n",
        "  #                                        xycoords='figure fraction', size = 12)\n",
        "\n",
        "      if icode == 2:\n",
        "        axes[icode].set_xlabel('Log(Age/yr)')\n",
        "      else:          \n",
        "        axes[icode].set_xlabel(r'A$_{V}$ (mag)')\n",
        "      \n",
        "  axes[0].set_title(code[0],fontsize=12)\n",
        "  axes[1].set_title(code[1],fontsize=12)\n",
        "  axes[2].set_title(code[1],fontsize=12)\n",
        "  \n",
        "\n",
        "\n",
        "  #Extinction plots\n",
        "\n",
        "\n",
        "\n",
        "  axes[0].plot(PDF.loc[:, PDF_columns[0]], PDF.loc[:, CB2018]/PDF.loc[:, CB2018].max(), color = black, marker = mark, \n",
        "          markevery= 10, linestyle = linestyle_, lw=lw_)\n",
        "\n",
        "  axes[1].plot(PDF.loc[:, PDF_columns[4]], PDF.loc[:, K2015]/PDF.loc[:, K2015].max(), color = black, marker = mark, \n",
        "          markevery= 10, linestyle = linestyle_, lw=lw_)\n",
        "  \n",
        "  axes[2].plot(PDF.loc[:, PDF_columns[6]], PDF.loc[:, K2015_age]/PDF.loc[:, K2015_age].max(), color = black, marker = mark, \n",
        "          markevery= 10, linestyle = linestyle_, lw=lw_)\n",
        "  \n",
        "\n",
        "  axes[0].axvline(percentiles.loc[0, percentile_16th_CB2018], \n",
        "                   color =mycolor[0], linestyle = '-', linewidth = lw_,label='16th percentile')\n",
        "  axes[0].axvline(percentiles.loc[0, percentile_50th_CB2018], color =mycolor[1],\n",
        "                   linestyle =  '-', linewidth = lw_,label='50th percentile')\n",
        "  axes[0].axvline(percentiles.loc[0, precentile_84th_CB2018], color =mycolor[2],\n",
        "                   linestyle =  '-', linewidth = lw_,label='84th percentile')\n",
        "  axes[0].legend(loc = \"upper right\")\n",
        "\n",
        "  axes[1].axvline(percentiles.loc[0, percentile_16th_KRUM2015], \n",
        "                   color =mycolor[0], linestyle = '-', linewidth = lw_,label='16th percentile')\n",
        "  axes[1].axvline(percentiles.loc[0, percentile_50th_KRUM2015], color =mycolor[1],\n",
        "                   linestyle =  '-', linewidth = lw_,label='50th percentile')\n",
        "  axes[1].axvline(percentiles.loc[0, precentile_84th_KRUM2015], color =mycolor[2],\n",
        "                   linestyle =  '-', linewidth = lw_,label='84th percentile')\n",
        "  \n",
        "\n",
        "  axes[2].axvline(percentiles.loc[0, percentile_16th_KRUM2015_age], \n",
        "                   color =mycolor[0], linestyle = '-', linewidth = lw_,label='16th percentile')\n",
        "  axes[2].axvline(percentiles.loc[0, percentile_50th_KRUM2015_age], color =mycolor[1],\n",
        "                   linestyle =  '-', linewidth = lw_,label='50th percentile')\n",
        "  axes[2].axvline(percentiles.loc[0, precentile_84th_KRUM2015_age], color =mycolor[2],\n",
        "                   linestyle =  '-', linewidth = lw_,label='84th percentile')\n",
        "\n",
        "\n",
        "  if len(CB2018) == 25:\n",
        "    plt.savefig(\"ext_pdf_{}.png\".format(CB2018[20:25]), bbox_inches = 'tight',dpi=100)\n",
        "  else:\n",
        "   plt.savefig(\"ext_pdf_{}.png\".format(CB2018[20:26]), bbox_inches = 'tight',dpi=100)\n",
        "  if i == 0:\n",
        "    plt.show()\n",
        "  plt.close()\n",
        "  i += 1\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAEqCAYAAADeRr5CAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXxU5b348c93shJIIAlhC6RBkTUKKhSIEQIoVwVE7Q+L1SqiV73687q0br/rT6O1LdqWam1rS2uLS3G9grvXDfyJikWEVhGtbELAkISQBULWeX5/nDnDZJ9JZjkz+b5fr3nBnHPmnGdCvsx3nuc530eMMSillFJKqcC5It0ApZRSSqlopYmUUkoppVQ3aSKllFJKKdVNmkgppZRSSnWTJlJKKaWUUt2kiZRSSimlVDdpIqWU6jYRmSkib4nIOhF53/P3y0QkUUSuEpHdIrIy0u1sj4hMEZHVIvKmiGwWkWdEJMNn/zgRWSsi6z37b/PZN0RE/igi74jIh573PqmdaxSIyNciUiMiT/tsv09EnvN5frOIlIjIy6F8z0qp4NNESinVLSJyFvA48B/GmEJjzOnAcuDPwHhjzApgZZCutVtECoNxLh8LgVeNMXOByUAWcJ/P/seAD40xBcAZwG2e9wyQB6QYY+YYY/KBd4FVrS9gjFkPXAsY4BKfXTM8D/u45cBaY8yCYL05pVR4aCKllAqYiLiA3wM/NcZst7cbY14Hno9YwwLzGPA3AGNMM7AJyPXZPwH4yLP/IPA1cLJn32bgdp9j/w58p4PrfAAkA6cCiEgfoALIFJHxnm19gcM9fUNKqfDTREop1R2nACOBN9vZtxT4rPVGEXlARCpFZInn+R9FpM7uaRKRwSLyus9w2m2e7X8FhgAPeoYQ7YTkUhHZICLvicgqEUnzbP+95zo/EZEXRGSPiBS1bo8x5mtjzFHPa4YC84BHfQ55FVjg2X8cVmL1see1B40x+zz7+gCXAn9q7wdljKkFNgIzPZumAf8DfOqzLR8r4VJKRRlNpJRS3XG85899rXcYY456enhab78V2OLz/GqgxOeQHwHrjDGzgH/Dk8QYYy73HHejZwhxk4ichjWMuMAYM9PTjuWe46/1XGcqsAg4Eyju6I145ir9C3jKGPPfPruuAEaLyA6s3qqbjTHvtnrtrcBeIAH4Px1dA3gPKPT8fSawzrNtps+29zp5vVLKoTSRUko5RQVwtohMMMYcAeZ2cuwS4GVjTJnn+SrgYhERn2NeNcY0G2O+Msb8uaMTGWMWAUOB6SLyS59dLwAfGWOOB04CbhaRU1q99gGsuVWfAe+KSFwHl1kHnObZP9YYs42WiVSuMWZXJ+9XKeVQmkgppbrDnheVHcRz/gL4b+AZEdmCNdTWkeHAbM9Q3zrgYeAAkOlzTJW/FzbGHAZ+DtwoIskiMg5rgvmDnv17gbeAW9p5rQHuxZo/NUtEJtntEpEHPYd9AKQA04F6z7b1wCDP3X61/rZVKeUsmkgppbpjM9bk6za9RiLyaOueGx8NQJLP8wE+fx9kjHnYGJMH/Bh4QkSOp317gVc8Q32FnjvrJhtjyv19AyJyVatNR4A4rInhiZ5tjT77GwF7HtZ8ERlm7/AMZdYDfY0xW3zadaNn/xHgE+A2js2zqgT+4dm2wd92K6WcRRMppVTAjDFu4CrgDhEZZW8XkcuBE/GZC9XKLqzSAYjITKxeGtvPfWoxfYyVdNlDdTVAiojMEpEbsMoqzBORdM+5xgCB1mC6yDPXCs+Q4JXAJk+C8yXWvKuLPPtTgXOxyhyAVS5hic/7/r6nrR93cr33gPm0nAv1HnAhOj9KqaglVq+0UkoFTkROB+4B4rG+mH0F3GaMKff0+PwfrB6eh40xPxWRscBzWPOhXgb+N1CJNbF7CFbvTBPQH3jMGPMbz3X+N1Y9pmrgCmPMVhG5xPP6Wqyk6z+NMf8SkQewkrwSYIWnRlN7bZ8P3Ox5moqVON1sjNnp2T8ZawK7AP2At4E7jDFNIpIH3A0M9rxvF3Bn68nora53FvC4MWaQz7bzPD+bEV38qJVSDqWJlFJKKaVUN+nQnlJKKaVUN2kipZRSSinVTZpIBYlncdI3PLc8/z/PQqa32pNhfY57S0QeabUty/M6IyL/EJHZnVznNBFp9lRitrdleyo8GxHZJtZCsomec1aKyOPtXGOdz6NSrEVYN3n2bxCR/j7nf0BEDonI48H8mSnVmoh8R0Re9MTQO57fz+s8+070+T3/WDzLq3j2DWkvhkQkw+d3/G8+x2eKyO9E5APP/r+LyIpO7jZUKqI0NhzMGKOPHj6wqjB/g1Voz952ItaE2vk+27KAMqx6N3HtnMcAZ3RxrYc8572+nX3vA2/7PL8EuKqra2BVkwbrDqojwOJW+/sAL0X656yP2H8Aa4FrfZ4XAp/5PM/1/A6P6uD17caQ/Tvu+Xs/4Ausyen2thTgReD5SP8M9KGP9h4aG859aI9UD4m1eOsjWIu3fmlvN8Z8hlUk0LcOzfewlsHoz7HlIgK91vHACqxbplu7Hqsg4PfEWnfse3Sw/lcrt3raXIu1vtj3W+0/B3gt0PYq1Q3fxaoCDoAxZh2ehYWD6Bag0njuCPRcpxYrDo4E+VpKBYvGhkNpItVzHS7eaoy52xjzPz6bzgCeBl6nbbLijwKsQHoKa7mJFlWljTFbsBKnXwHLgHuN5ytFR0RknTHm7z6bnsVapiPNZ9v3sCpOKxVq3wC3ikhfe4MxZlmQr/E92o/Xr4wxlwX5WkoFi8aGQ2ki1XMdLt7qS0QGA6XGmAasROh8EYkP8FrfA542xvwDq2DgonaO+S+s6stDjDGbOzjPr+TY0hqtvYpVx2ehp919gRRzbE0zpULpeqzFiveJyF/FKtoZbMfTRbwq5UAaGw6liVTweBdLFZFTPYnK53JsEdTvYfX2gFWIMBmY4/fJrWG9QcYYexX7VbTfq1ULlANnich3Ojjdj4xnCYvWO4wxR4FXfM493/NcqZAzxrwD5GAtEZMLrBWRPwbj1K2eS4snIs95JuvuDsK1lAo6jQ3n0kSq53Z6/vQOsxljNnmSlC3AQM/mc4F7PL1Ar2MlPO3Nc+rIDOBUn56kBcBUEclpddwdnscerCG+TrWXTGElfHM9dxyeD7wQQDuV6hFjzBFjzJ+NMbOAWcCVInKcny9v4Ng6eQCISBLHFgoGK2ZbD4svAu4EOvryoVTEaWw4kyZSPfcpsBs4q6MDRGQIsN0YM9OnJ+g/sYb3Ejo7uYj8wPPX/wXMMMcWQ52KtRzHhT7HjgKOM8b8N3Ad8D0R6bLXq50u4tewgutSINEYU9HVOZQKBmlVGsQY8x5wEOsGjc5eZ8fJXqw5i75Gerbb/ptO4lUpJ9LYcC5NpHrIWKu+/wfwX2KtvwVYtZ2A0YAbKwl6tdVLXwf6Amd2cYmrRCQOyDHGlLTat4aWvVo/5dgdeO9gzcX6jR9zse5p9Z7qsIYff0LgC8Eq1RNniMh37SeeJN+NNSewM1d5/lyN9S09w/N6F9adsi/6HPsLIFVEfuRznXjgdNoOcyjlFBobDhXoZGfVDmPMGyJyEfBrEUnGmv8UhzVE9i1wP7BTRN70JF5gJT1NwO9F5C6soT+Au0TkGp/TZ2EtljpGRK4wxjwK1jwsz2vGi8hqrMVTRwE3Yt3ZcSIwARgP/F1EbsKarNjeNSa087aewUoAV3fvp6JUtzwA/EJE3Fgx5AbOM8YcFWvBY3sB4odFpL3bse/Gutni7yJS6tm22hjjnednjKkRa7Hln4jIx1glSgZg9S5PDMm7UqrnNDYcShctVkoppZTqJh3aU0oppZTqJk2klFJKKaW6SRMppZRSSqlu0kRKKaWUUqqbNJFSSimllOqmiJU/GDhwoMnNzW133/79+zlsrHUZR2e3rDVmDlgrpMjg4ezfv59hw4Z59/2r0loiaPSAY4VZWx/TkWAep9eM/HGdHbNp06ZyY0xWlxeKkK5iw35f/9pXBbQfIw0NDSSNaFvw2DdGnPbvoteM/DVjJTYCEa7X2fHaT460+zrfz7Y2r63cR0N9PXmD/S1i3v126us61lF8RKz8weTJk80nn3zS7j4R4Yz/supAvnXf/Bb7mn59GwDxN92PiODb/jPW3A7A2+cdWxC79TEdCeZxes3IH9fZMSKyyRgzucsLRUhXsWG/rzPvtMq/tBcj761bx5wXP27zet8Ycdq/i14z8teMldgIRLheZ8fr2z9d0O7rfD/bWjtjze2sW7eOpgc3hLyd+rpOX9NufET10N7dd98dlGNCcVy4rxnt7Q/Fcb1ZR9/afTn530VjI/DjentcdPf9h/vn1t3r+RPTwbxed4X738EJv/dR3SPVWns9UpHQ3QzZSaL9PfSGb92d9UiBM2Mk2n+vIPrfQ6zGRrh1t0eqdbzanBa30f573pWY75FyQobZE9Hefoj+9xDt7e9ItL+vaG8/RP97iPb2O0Ws/xz1/fnP70RKRKaIyHYRWdLFcReLyCYR+UREfiUiEmijioqKAn2Jo0R7+yH630M426+x4b9obz9E/3uI1dgIt2j/PeiKvj//+XXXnoicDywCqro4Lg/4FZAHVADvANcCv+tZM9tyu92Ul5dTWVlJc7O1DvCdI6x1f7dt2xbsy6kok5yczPDhw0lISAjpdZwYG0o5gcaG6i38LX+w0RizWkTWdXHcFcBrxphyABH5C3ArIQiI4uJiRITc3FwSEhIQEVyHrNtHx6S3vX1U9R7GGA4ePEhxcTEjR44M9eUcFxtKOYTGhuoV/BraM8YU+3m+KcCXPs+/ACaISJ9AG9aVI0eOkJ2dTWJiIlHQC6zCSETIzMykrq4u5NdyYmwo5QQaG6q3CPZk88G07MatBAQY2JOT/vOf/2TIkCG0vlvD5XLkXHnlAA5MrkMSG6098MADbPrkE9xudzBPq1QohSU2okFjQwM5OTmsX78+0k1RAQhFJtLe/YRtPtX279+PiHgfnU38uu+++zhw4AC/+c1vgthM1VsUFRW1+F0DAi+DGxxBjw1fH374IXfccQeHDx+m9MCBHjZV9Qa9JTaiRXl5OXv37uWVV16JdFMU/sdHsBOpUmCAz/MBWAFS1vrAYcOGYYzxPjoKiK+//prnn3+e/v378+yzz9LU2BjkJofOxo0bGTVqFCtXrmyzr6qqiiuvvJKCggImTpzIWWedRVNTU7vnWbFiBbm5uSxZsiS0DXaA3bt3M378eAD+/ve/M2nSpG4XorMVFRW1+F0D9ve4oYELemy09sADDzB48GD69u1LcXFxTNeAUcHRW2IjWhysOAhY//epyPM3PoKdSG0Exvg8Hw9sNcYc7e4J33jjDYwxPP7449TX11NeXt7jRobD6tWr+fWvf03//v3b3X/xxRdz/vnns379ejZv3kxycrJ3OKaoqKhF0nTVVVfFbBKVm5vLunXrWjz/8MMPAfjud7/Lgw8+GKGWBV3QY6O1TZs2MXv2bLKzszly5IjevaqiRchjIxq43W4OHapERPhEh+ejSo8SKREZKCLvi0iGZ9OfgXNEJFNEXMAS4A89ucb27dvp168f8+fPp2/fvhw5cqQnpwubKVOmsGrVKlJTU9vs27RpE2VlZcybNw+w5nqtWbOGxMTEcDfTkQYMGND1QQ4XjtjwVVZWRnFxMSeffDJpaWkAfPbZZ8E6vVJBE+7YiBbV1dW4m5tZtGgRNTU1fPXVV5FukvKTX4mUiJzquYV1EnC7iLzg2dUHGAukABhjPgd+DLwJfAxsAX7fkwZ+/fXXjBo1CpfLxejRo6mtre3J6cJm+PCOSzC888475Obmcv3115Ofn8+CBQv4/PPPAXjmmWdYuXIlb7zxBoWFhfz0pz/1vq6uro6rr76aU045hXPOOafdu9IaGhooLCxERPj5z3/O3LlzmTBhAsuWHVtaoLGxkVtuuYX8/HwKCgq45557MMa0eO3vfvc75s2bR2pqKuvWrePw4cP8+7//OwUFBRQUFHDFFVdQVmb1vO/YsYO5c+cyc+ZMTj/9dG+P0ksvvcTYsWOZOXMmt912G9OmTeO0006jtLQUgMsvv5ySkhJuvPFGCgsL2bRpE3PmzEFE2L17d7s/u8OHD7N06VIKCgrIz8/nD3+I7P+3kYwNX5s3bwbglFNOISUlBRHx/k4pFQlOiY1oUe/5//ySSy4BdHgvmvhVR8oYswkobGf7XiCr1ba/AX8LRuPASqROPvlkAMaMGUPt0W/bPe7Hj37Ehq+t4cs+8VuCdfkWJo7M5JdXTO/xeXbv3s0LL7zAc889x8MPP8yTTz7JrFmz2LlzJ9///vfZtm0bu3fvbjO36v333+fzzz+nf//+TJw4kdWrV3PRRRe1OCYxMZF169YhIlRWVvLmm29SUVHBhAkTOOWUU5g7d651Z9emTbz//vu43W5mzpzJ8ccfzyWXXOJ97cGDB3n11Vd56qmnSEtL4+abb6a5uZn169fjdrs577zz2Lp1KwUFBcyfP59bbrmFpUuX8s9//pPZs2eza9cuzj33XCoqKrjuuutYuXIl999/P+eccw6PPvood9xxB3/9619Zu3YtDz74IIWFhYCVZHZ2x91NN93kbUdNTQ0TJ04kLy+PgoKCHv+7dEckY8OXnUhNmjQJ2fIGfVJStEdKRZRTYiNaNHjm/86YMYOEhAQdmo8ijq4fYIxh9+7djBo1CrASqbq6OkyUjx3X19eTk5PDeeedB1jfQBITE3n55Zc7fd3UqVNJT0/H5XKRl5fHrl27Oj1+8eLFAGRkZHDOOefw9NNPA7By5Uouu+wy4uLiSEhIYNGiRTzxxBMtXrtw4UIALrroIiZNmsTjjz/unaflcrn41a9+xfjx4/n444/ZsWMHP/zhDwE46aSTyM7ObnHXyZgxY7yFMU866aQu290Rt9vNE088wdKlSwFITU1lwYIFbdreG23evJnc3FzS09MB6Ns3RXuklIoijQ0NuFwu0tLSGDRokLfnXjmfv5XNI6Kuro6mpiZOOOEEwPpAZvPrHD3adg7iL6+YzldRUtk8PT2dwYMHt9iWnZ1NcXHn9evsuS8ASUlJNDQ0dHkdW2ZmpreHori4mOXLl/PXv/4VsIbLWs9L8p0kX1ZWRn19PVlZx75E2v8ma9euRUQ488wzvfvq6+upqjpWFsa33cnJyV22uyN2O2699Vb69LFq9VVWVjJp0qRunS+WbN26lZNOOsn7vG/ffuzcuZMjR47Qt2/fCLZMKeWPhoYGb4HpQYMGcUBLmEQNRydSdsLkm0gdgqiZJ9WRSZMmsWbNmhbbysrKGDYsuCVcKioqvKUDysvLGTp0KAAjRozgzjvvZNGiRYDV01NZWdnhebKyskhKSqKsrIxx48YBVj0Xl8vFiBEjSEhIaHHn3ZEjR0JSLNVux29/+1umTJkCWPO9ov33IRj27t3LrFmzvM/79u2LMYatW7fy3e9+N4ItU0r5o6GxkQTPDUeDBw/WHqko4uihPTuRsof2Ro8eDUR/IrVw4UJqamq8ycfatWupqalh/vz5gDVkVVtbizGG888/v9vXef755wE4ePAgr732mneob8mSJaxatcq72PNjjz3WYlJ7ay6Xi0svvdQ7Z8vtdnPFFVdQUlLC1KlTycnJ4YUXrHmkTU1NnHfeefzrX//yq432e127di0PPfRQp8fa7fAdyrvvvvt4/PHH/bpWrKqurqaqqooRI0Z4t6V4euy6O4yqlAqvxoYGEhOtRda1Ryq6RCyRsivUdlZQrb6unqSkJO8wWL9+/UhMTAzLGmo9tWnTJgoLC9myZQvLli3jggsu8O5LTU3l5Zdf5vbbb+e0006jqKiI119/3Tu8tnDhQrZu3Up+fj75+fmsWrXKeyffI488wooVK3jjjTdYuXIlq1at6rANgwcP5t/+7d84/fTTueGGG5g7dy4At9xyC2PGjOG0005j1qxZvPvuu/zkJz8B8B6zePFi3n33Xe+5li9fjsvl8t61N3/+fCZNmkRcXBwvv/wyK1asYObMmcyaNYuLLrqIiRMn8u6777Js2TK2bNnC3XffzUsvveR9H8uXLwfg6quv5sc//jF33HEHZ5xxBnPmzPFe/7333uPGG2+kpKTE23u2fPlyamtryc/PZ+bMmVRVVXHttdf69W9iV6klctWb/eJPbPjau3cvQItEKikpqcU+pToTq7ERTRoaGkhMaNkjpUV1naGr+JBI/UNNnjzZtF47z9eZd77Cl19+iWvjw3zzzTfe7R9fcgaJiYkk3/Ir7zCTLVrmSIWDiLBr164eVwWPZtu2bWvzOwIgIpuMMZMj0CS/dBUbtjPvtCb0/6ggnrPPPpsPPviA/Px8mn59GwAZdz/C5Zdf3qKn74w1twPw9nnL2p5Q9XqxEhtOZMfrW/fNb7PP7Xaz7oJ8RuTkcMJvnuWXv/wlt9xyC1VVVaSlpWncOkRH8eHoob3GxkYGDRrUYltSYiIN9fURapFSzrNnzx6gZY+U/Vx7pJRyvoqKCowxJCZYQ3v2KIwO70UHRydSDQ0NbRKpxKQk6rt511dvYBfVBGt4bN++fZFtkAq5vXv34nK5vDcT2DSRUio62AmTvbqF/bmnE86jg6Pv2mtsbGTQiLY9Uo0NDTp23AG7IKfqPfbu3cuwYcOIj28ZziNGjGDLltAUp1VKBU/rREp7pKJL1PVI2ZNo7TvOlOrt9u7d22ZYD6xE6sCBA9TrULhSjmYnTAkJx+7aA+2RihaOTaSam5oxxrQd2vNk7JpIKWXZs2dPh4kUoMO7SjmcvW6p/flmFz/WHqno4NhEqqHRmgelPVJKdW7//v3tLpJtJ1I6T0opZ7NXgojzDM8nJCSQkZGhPVJRwrGJVGODtYBjRz1STU1NYW+TUk7jbnZTW1vbJk5AEymlokV1dTWuuLgWC7YPHDjQ21OlnM2xiVRHPVIJCQmIiPZIKcWxOPFdB9FmLzn07bffhrVNSqnAVFVVER8X12Jbenp6izVLlXM5trJ5Rz1SiJCYlOT4RGrJkiUUFha2eBw+fNi7f8uWLUyfPp3TTjuNBQsWcPDgwQ7Pde+99zJkyJCYrObb2vr1673V1V966SXGjh3rLefQU7FYvbmxsYM4wVoJoE+fPjrPQnUpFmMjmlRXV3uH9WwDBgzodA1UFT5dxUfEEqlhw4ZhjOkwIDr7pp2YmOj4RApg3bp1LR79+vUDrLsRFy5cyM9+9jM++OADTjnlFK655hrv65YsWdLi53LXXXdx1llnhbv5YSEi7N692/v8tNNO47nnngPg3HPP5fbbbw/atYqKiuyyGfuDdtIQ6Co2fNlfONqLE3sVeZ1noboSi7ERTaqqqtqUL0lPT+fQoUMRapHy1VV8OLaOVGNDI/Hx8d45Ub4SExI47HZHoFXB8frrrxMXF8esWbMAuPLKK8nNzaWsrKzdD8TeRETo379/pJsRNTr7wgFWPRrtkVLK2aqrq4kfnNpim/ZIRQ8Hz5Fq9NbUaC0hMSEqeqSuueYaTj/9dM4++2zWrl3r3b5x40bGjh3rfT5ixAhSUlL49NNPeeihh7wLEhcWFvLoo496jysvL+eHP/wheXl5XHbZZe1ec8+ePUybNg0R4aGHHmLOnDnk5eXx2GOPeY85fPgwS5cupaCggPz8fP7whz+0ee3KlSuZO3cuSUlJ7N69m9LSUhYtWsSMGTOYPn06N910E0ePHgWsBZpnzJjBzJkzmTNnDl9++SUAK1asIDc3l8WLF3P11VdzyimncM4553gXnT777LMBqwJ7YWEh33zzjff6HSktLeWCCy5gxowZ5Ofns2bNmoD+TWJNZ0N79nbtkVLK2aqqqjoc2tPi087n2B6ppqZG4uM7SKQSrKE9Y4z3Q7f5uT8ybPc267XxSSFpkww/jrhFV/t17Lhx45g9ezZTpkxh48aNzJo1i/Xr1zNp0iQOHDjQptdlwIABlJaWcsMNN7B582Zyc3PbdF9/8MEHbNiwAWMM2dnZfPTRR0yfPr3FMTk5OTz99NOMHDmSxMRE3nnnHXbt2sWECROYOnUqY8eO5aabbqK5uZn169dTU1PDxIkTycvLo6CgwPtal8vFm2++yfLly0lKSuKSSy7xDrs1NDQwbdo0Dhw4QHp6OmeddRbPPPMMs2fP5tVXX2XhwoVs27aNq666iv379/OnP/2Jzz//nP79+3PSSSexevVqLrroIl5//XVEhKefftq7uLJ9/Y5cfPHFTJ8+nXvvvZdvv/2WvLw8Jk2a1GsXZ25saCQ5OZm+ffu2u3/w4MFs2rQpzK1SSgWiurq6zWTzAQMG0NjY6P3CqpzLsT1STU1NbcaMbfbCjk7ulbrtttuYMmUKAFOmTGHevHn88Y9/9O5vr9elq28es2bNIikpieTkZEaPHs2uXbs6PX7x4sUAjBw5kunTp/Pss8/idrt54oknWLp0KQCpqaksWLCAJ554osVrFy5cCMDNN9+M2+3mrbfe4vLLLwesOWp/+ctfyMrK4pVXXqFfv37Mnj0bgHnz5lFSUsLHH3/sPdfUqVNJT0/H5XKRl5fXZbs7sm/fPt5++22uuOIKAIYOHUpBQQFPPfVUt84XCxoaG8nKyuqwF8/ukXJH8VC4UrGuozlSgM6TigIO7pFqIjm5T7v77CE/32QrbtHV7D9UDMCY9LbFCSMtJyeHrVu3AtaH2+bNm1vsr6ys7HB4xpaWlub9e1JSEg1dLN5sByJAZmYm3377LWVlZdTX13PrrbfSp08f77UnTZrU4rW+PWbFxdbP1Xcejn18cXExFRUVLe6sy8rKanEXom+7k5OTu2x3R+x2XHrppd7Eoby8nBNPPLFb54sFjY0NDO/k92bw4ME0Nzdz6NAhMjMzw9gypZQ/3G43NTU1xLXTIwXoPKko4OhEKqGDHqmEKCjK+cADD3Drrbd6nx84cMBb12fKlCk8/fTT3n179+6ltraWU089NahtqKioICMjA7ASjry8PLKyskhKSuK3v/2tt8essbGR2traDs9jF3YsKysjJycHgJ07d5Kens6IESMYPnx4i4WSq6urSU5ODup78W3H888/703q6urqHP17EGqNnh6pjvgufqsTIygAACAASURBVKqJlFLOc/jwYYwxbXqkNJGKHo4c2jPG+DW05+QP0OXLl3sn+e7atYsXX3yRH/7wh4A1ybqpqYn33nsPgL/85S9ccMEF3g/E1NRUamtrOXLkCBdffHG32/D8888DVtKzYcMGLrzwQlwuF5deemmLobz77ruPxx9/vMPzDBs2jDPPPJOVK1cCUF9fz4UXXkhDQwPz58+nvLycjRs3AnDkyBFmzZrldyG5fv36UVtby5NPPultb1ft8G37Nddc02Iif2/T0NB5IqWLnyrlbNXV1QDtTjYHTaSigSMLctbW1loZekIHPVKeRMq+Y8mJfvzjH3P++eczc+ZMFi1axG9/+1tmzpwJWMNya9as4fbbb6egoIBPPvnEe+ccwA9+8ANeeuklZs2axbx581i+fLn3Tr7Vq1dz7733smXLFpYtW8a7777bYRtEhDPOOIMFCxbwyCOPeO8UXL58ObW1teTn5zNz5kyqqqq49tprqaio8M6rKiws5IsvvvCe68knn+Szzz5jxowZzJ49m9tvv53BgweTlpbGa6+9xi233MLMmTM5++yzueeee8jKymLVqlWsXLmSN954g0ceeYQVK1Z438eqVasAuO6661i8eDF/+MMfmDp1aovrv/jiiyxbtowtW7Zw/fXXe9vx4Ycfcvrpp3P66aczevRoFixY4Ne/SSwWHWxsbPC7R0qpjsRibEQL+0tne5XNQRMpJ+gqPiRSt1ZOnjzZfPLJJ+3u27dvHyMXP8jo0aP5/NF/b7Gv6de34W5u5t28M5g4cSJDhw717vvKwXOkwmn37t2MHDmy1982u23bNsaNG9dmu4hsMsZMjkCT/NJZbPiac8dLvL/+fX4yfyC33Xabd3vTr62/x990P2VlZQwaNIjf/OY3XH/99Zyxxipw+vZ5y0LTeBXVYiU2nOjMO18B4K375rfY/tFHH5Gfn0/pnVeSnpFB/E33A3hj9+GHH2bNcOuzTeM2sjqKD0cO7dl3KXQ0tGcv7ujkoT2lQq2xyeqR7WzuU2ZmJi6XS3uklHIou0eq9dCefcOP9kg5nyMTKfsXp6NECiAuLs7RQ3uRsmfPnhbDY04uEaF6pqnR+iJh31DQHpfLRWZmJuXl5eFqllIqAPYcqdZDe4mJiaSkpGgiFQUceddeVz1SYH1AaI9UWzk5OWzYsCHSzVBhYPdI+Za5aE9mZiYVFRXhaJJSKkDeOVLtfN6lp6d7EqmBYW6VCoRfPVIikiwiK0Vkg4h8IiJzOzguVUSeEJGNIvJ3z2v6Bdoof3ukNJFSkRbu2PDlT4+Uvd+3rpdS4RLJ+IgW3rv2WvVIgXXnnhbkdD5/h/aKsCamTwN+ADwtIoPbOe7/At8Bpnke3wHuDLRR/iRS2iOlHKKIMMaGr8Ym/xIp7ZFSEVREhOIjWlRVVSEiHSZSOrTnfF0mUiLiAq4EHgUwxvwL2Axc0s7hE4C/G2OajTFu4O/AyYE2Sof2VDSIRGz4avJzaE97pFQkRDo+okV1dTWpqanQzjJPmkhFB396pI4DMoEvfbZ9AbR3i+xrwGwR6SsiKcBs4ON2jutUZWUlcZ478zoSFxeH2+3WNcRUJIU9Nnw1NTYhIh0uWGzLyMjQHikVCRGNj2hRVVXVZhF7W1pamnfoTzmXP4mU3Q3rW6q6EmizwJcx5nfAOmAnsBv4FLg30EYdOnSo094osHqkwNnVzVXMC3ts+GpsaiQ+Pr7TLxxgDe0dOXKE+vr6nlxOqUBFND6iRXV1dYv1SH1pIhUdAil/0Lq6Y5v/vUXkv4CJWOPbOVjfSJa0dzK7Qq398K1UW1lZ6Xci5dTb+7dt28asWbMoKCjg5JNP5v7772/3uLvvvpvc3NxOz3XvvfcyZMiQmKrm25H169czd641H/Wll15i7NixLRZE7g67Kq39IPjVm8MWG76ampq8Vf47Y8+h0l4p1VoYYgOCGB/+xkY06axHqn///ppIRZC/8eFPImUv0jXAZ9sAn+2+/hP4ozGmzhhTB/wB+El7Jx02bBjGGO8j0ETKnpjn1B6pyy67jPz8fNavX8/bb7/N/fffzxtvvNHimAMHDvDnP/+5zWuXLFnS4udx1113cdZZZ4W6yREhIuzevdv7/LTTTuO5554D4Nxzz+X222/v8TWKiopa/K4B+3t8UkvYY8NXY2PH61H6sgt2aiKlWgthbEAI4sPf2IgmVVVVnfZI1dXV6RSWCPE3PvxJpHYAFcAYn23jgY3tHJsI+FbJbARS/Wqxj0B6pJyaSG3dupXp06cD1gfZCSecwObNm1sc85Of/IT/+I//iETzHEtEOvx25kBhjw1fTU2NxMf73yOlE85VmEU0PqJFdXV1p3OkwLkjL8rSZSLluYPiT8BSABE5AZgE/E1ExonIOyJi37f5NvB98cC63XVtoI2qrq4mLi66E6l58+bx8ssvA7Bz5062bt3K1KlTvfu3b9/O/v37KSgoaPG6hx56yLuwb2FhIY8++qh3X3l5OT/84Q/Jy8vjsssua/e6e/bsYdq0aYgIDz30EHPmzCEvL4/HHnvMe8zhw4dZunQpBQUF5OfnexdM9n3typUrmTt3LklJSezevZvS0lIWLVrEjBkzmD59OjfddBNHjx4FYNOmTcyYMYOZM2cyZ84cvvzSmlu6YsUKcnNzWbx4MVdffTWnnHIK55xzDnV1dQCcffbZACxevJjCwkK++eYb7/U7UlpaygUXXMCMGTPIz89nzZo1/v2DhEAkYsNXU1MTCdojpRwq0vERLbqabA7O/ZxTFn8rmxcBfxCRDZ7XXGSMKRGRXGAskAA0A9cBDwMbsMbB9wBXBdqo6upq+sa3ranhq/UcqR+tX8FHJdsASIlPCvSSfpk48Dh+VeDf23n00Uc599xzOf7446moqGD58uXMnj3bu/+uu+7innvuadNLcMMNN7B582Zyc3PbdFt/8MEHbNiwAWMM2dnZfPTRR95eL1tOTg5PP/00I0eOJDExkXfeeYddu3YxYcIEpk6dytixY7nppptobm5m/fr11NTUMHHiRPLy8igoKPC+1uVy8eabb7J8+XKSkpK45JJLvMNuDQ0NTJs2jQMHDpCens5ZZ53FM888w+zZs3n11VdZuHAh27Zt46qrrmL//v386U9/4vPPP6d///6cdNJJrF69mosuuojXX38dEeHpp5/2zhOzr9+Riy++mOnTp3Pvvffy7bffkpeXx6RJk7qcZxZCRYQxNnw1NjYSn9B1CLfokep4WT6lQqGICMVHtOhssrmdYGmPlLP5NdncM269xBgzzRgz2Rjzpmf7BmNMtmdMG2NMqTHm+8aYqcaY7xpj/pcx5ttAG1VTU+NXj5STFy6+4IILmD59Ojt27OCf//wny5cv59NPPwXgk08+ISEhgRNPPDGgc86aNYukpCSSk5MZPXo0u3bt6vR4e829kSNHMn36dJ599lncbjdPPPEES5cuBSA1NZUFCxbwxBNPtHjtwoULAbj55ptxu9289dZbXH755YC1BtRf/vIXsrKyeOWVV+jXr583SZw3bx4lJSV8/PGxO5enTp1Keno6LpeLvLy8LtvdkX379vH2229zxRVXADB06FAKCgp46qmnunW+YAh3bNiamppobm4OaGhPe6RUuEUqPqJFY2MjR48e1R6pKOe4tfbq6+upr68nvoseKbAKdtq/YL8quIqvDhUDMCZ9eEjb2JVt27bx9ttv87e//Q2AESNGcOaZZ/KLX/yCp556irvvvpvf/e53AZ/X91tLUlISDQ0NnR7vW6gxMzOTb7/9lrKyMurr67n11lvp06cPYM1JmzRpUovX+gZ2cbH1c83KyvJus48vLi6moqKixZ11WVlZLXrafNudnJzcZbs7Yrfj0ksv9Q7/lZeXB5yQxgK7SF+CHz1S/fr1IyEhQRMppRzGviMvLS0Nmmva7D82R0oTKSdzXCJVU2P9MnXVI2Ud48z19uxEwffW9ISEBKqrq6mpqWH79u0sWbIEsD4QS0pKKCws5Pvf/35QJ59XVFR4eyPKy8vJy8sjKyuLpKQkfvvb3zJlyhTA+lZUW1vb4XlGjBgBQFlZGTk5OYA17ys9PZ0RI0YwfPhw1q1b5z2+urqa5OTkoL2P1u14/vnnvUldXV2dI38HQs2f6v82EfGpbt75cjJKqfCxFyzu378/VOxrs9/+QtvUpEN7ThZIHamwsDN0f3uknDh2PHbsWLKzs71DTjU1Nbz00kvMnj2b1NRUvvrqK9atW8e6det48MEHGTJkCOvWrfMmUampqdTW1nLkyBEuvvjibrfj+eefB6ykZ8OGDVx44YW4XC4uvfTSFkN59913H48//niH5xk2bBhnnnkmK1euBKxewwsvvJCGhgbmz59PeXk5GzdaN+IcOXKEWbNmef+D6Eq/fv2ora3lySef9La3q3b4tv2aa65h7dpeMSe1Bbt3yZ9ECnS9PaWcqEWPVDu0Ryo6OK5HqrOVsFuLj493ZLXmpKQk1qxZw80338xTTz3F4cOHOf/887nhhhtaHHfvvffywgsveHukHnjgAb773e/ygx/8gMsvv5x169Zx4403snz5ct544w2Sk5OZOHEin332GVu2bGHZsmXk5OS0mMTuS0Q444wz+Pbbb3nkkUcYO3YsAMuXL+fGG28kPz+fhIQETj75ZO666y4qKiq886oKCwv5/e9/z/jx4wF48sknue6665gxYwbNzc3cfvvtDB5sFS5+7bXX+NGPfuSttXHPPfeQlZXFqlWrWLlyJXV1dTzyyCPExcV538fo0aP5wQ9+wHXXXcfixYtJS0vjqaeeYtGiRd7r33TTTSxbtoySkhKuv/56Hn74YZ588kmuvfZaVq9eDVh3/i1YsCD4/4gO58/C3r7sHqk4RoWyWUqpALTokWrHsTlSzuswUD58i02F8zF06FADmLvvvtv4eu+99wxgJl/7hDnjv142rTUuv9U0Lr/VfPHFF2bXrl1my5Yt3n1fVuw1X1bsbfOa3mbXrl3G+qft3b744osWz++++26DVWV5v4nQ770/j45iw9ezzz5rEgp+bKbf+HSbfXaM+Fq4cKE56aSTzJzVt5k5q2/r/Aenep1Yig2nOuO/Xm7zmfbSSy8ZwGzcuLHduHW73SY+Pt6MXL5I4zaCuoqPiA3t2RVqW9/i7+2R8mNoz6lzpJQz2VVqCW715qDrKDZ82d9kA+mR0qE91ZFYio1o0lWPlF2guFl7pCKqq/hw7hwpPyabx8fHY4xx5DypSNmzZ0+L4Tn92cQmbyLlR5wAPpPNlVJO0dUcKXufdhg4m3PnSPk52RysYmX+zKnqDXJyctiwYUOkm6FCzE6k/IkTsCabHz16FLe7GZdLY0UpJ+iqRwqsRKpEJ5s7WtT3SIEWK1O9T1VVVUBfHuwyGI2NGitKOUV1dTUJCQkkJXW8Gkf//v11srnDOTKRcrlcuOK6bprT19tTkeMZz45ZVVVVfs+PgmPr7TU1NXZxpFIqXOx19jpbXzQtLU2naDic4xKpmpqaTseLbX379qW8vBywCkoqZTPGcPDgwZAUBXWKQBMp7ZFSynk6W2fPpnOknM+Rc6RSU1O7PG748OHs2LGDnTt34na7KS0t5UCtVe3ZndK21L7qXZKTkxk+PLJLBYVSVVUVcX26k0jplw6lnMLukepMWloaTc17wtQi1R2OTKT86ZFyuVwMHz6c0aNH87Of/Yw77riD69fcDsDb5y0LdTOViqiqqiriU/2fI3VsaE+/2SrlFP583vXv35/mBh3aczLHDe35m0gB9OnThz59+mh9HNXrdH9oT3uklHIKf3uk3G43brc7TK1SgYpYIrV//35EpN2CnP4mUqD1cZT/ioqK7EmdwyLdls50FBu+qqqq/K4hBZCSkkJSUpJONlftiqXYiCb+zpECXW8vkrqKD0dWNg8kkdLFWJW/YqV6szHGmiMVQI+UiJCRkaGTzVW7YiU2oo0/PVL2fi2BEDlRWdlce6SU6lhdXR2NjY3E+1mM05aZmUmj9kgp5QjGGO2RihGOS6Ti4uJIT0/3+3jtkVK9TaDr7NkyMjJo0h4ppRzh6NGjNDU1+TVHCrRHyskcd9fenj3WbZ5n3vmKX8drj5TqbQJdZ89mDe1Vh6JJSqkA+bPOnu9+vePWuRzXIxUou0cq1itZK2ULdJ09W0ZGhv5nrJRD+LPOnu9+rW7uXFGfSFnfshs5fPhwpJuiVFh0d2gvPT1d50gp5RDaIxU7oj6RsgsN6vCe6i16MkfK3az1aJRyAn97pHSyufNFfSJlFxrUCeeqt+juHCn7Jg79ZqtU5PnbI5WcnIzL5dLJ5g4W9YmU3SOliZTqLbo7R8pOpLS6uVKR52+PFFh3s2uPlHM5rrJ5oLRHSvkrVqo39+SuPdAeKdVWrMRGNPG3RwqsYXztkYqcqKtsHihNpJS/YqV6c1VVFampqSCBnVeH9lRHYiU2oon9hcifRCouLo4m7ZGKmKirbB4oO5HSyeaqt/BnWYn26NCeUs5RXV1NSkqKXzeNxMfH06w9Uo4V9YlUUlISffv21R4p1Wt0N5HSoT2lnCOQOI6Lj9O4dbCoT6TA+oDQREr1Ft1NpI4tfqr/ISsVaYGsKxsfF68FOR3Mr0RKRJJFZKWIbBCRT0RkbifH5ovIOyLy/0TkCxG5PnjNbZ8mUipSIhEb1dXV3Uqk4uLiiI+P16KcKmyc/tkRSYF8IYrXHilH8/e2nyJAjDHTRGQ0sEFExhljDvgeJCIjgQeBc4wx5SIyAVga1Ba3w15vLznUF1KqrSLCHBtVVVWMGjWK8m68Nj4hXhcuVuFUhIM/OyIpkB6pOO2RcrQue6RExAVcCTwKYIz5F7AZuKSdw28G/mqMKfccu9UY86PgNbd99np7SoVTpGKju0N7AAnxCfrNVoVFNHx2RFJVVZX/Q3vxcbjdburq6kLcKtUd/gztHQdkAl/6bPsCmNzOsXOARBF5TUQ+EJFfikjIO4p0aE9FSERioyeJlA7tqTBy/GdHJB06dMh7A0hX4jw14+zaU8pZ/EmkBnv+rPLZVgkMaufYXOAa4FKgEBiP1V0bUscSKRPqSynlK+yxUV9fT319ffcTKR3aU+Hj+M+OSDHGcPDgQb8TKbtEgiZSzhTIXXuts5T2ygEmAauMMeXGmEasQFji6eJtwa5Qaz96UmAtIyODxsZGHUNW7bKr0toPgl+9OWyxEUgRv/YkxMfr0J7yCkNsQBDjI5ifG5FUW1tLQ0ODd4mzrsTFWctB2fGvwsPf+PAnkSr1/DnAZ9sAn+2+DgG+kwiLsQJkYOsD7Qq19qMnAWH/MjbqN23VDrsqrf0geNWbwx4bgazP1Z74+ASamhrtn4Pq5UIYGxCC+Ajm50Yk2QWktUfK2fyND38SqR1ABTDGZ9t4YGM7x26hZbdtFtAAhLTs+LFCgzr3Q4VV2GOjx4lUQjxut6G2trZbr1cqAI7/7IgUO5EKtEdKEyln6jKRMsa4gT/huRVVRE4AJgF/E5Fxnrof9jL0fwIuEpEUz/OlwJPGmJCOudmJlPZIqXCKRGz0NJFK8HyzPXToULder5S/ouGzI1Lsm6P8TaS0R8rZ/J0jVQSIiGwAngIuMsaUAP2BsUACgDHmGeBJrFohHwDNwI3BbnRrxxIp7ZFSYVdEGGMjGEN7oIt8q7ApwsGfHZES+NCezpFyMr8Kchpj6oAl7WzfAGS32vZz4OfBaJy/7KxeJ9GqcAt3bLRMpIoDfn289kipMHL6Z0ekBD60pz1SThYTa+15V7XXOVIqxvV4aC9BEymlIs3uEfa3R8rlcuFyuTSRcqiYSKSSk5NJSUnR+jgq5vW0/IEO7SkVeQcPHqRv374kJSX5/Zq4+Dgd2nOomEik4FgtKaViWVVVFX379vUO0QVKh/aUiryDBw/6Paxni4+L1x4ph4qZRCozM1PLH6iY15PlYcCatCoimkgpFUEVFRV+D+vZ4uLjNJFyqIglUnaF2mAVVLN6pHRoT3XMrlJLaKo3B01nsdHTRAqE+Ph4HdpTLcRCbEQT7ZGKLl3FR8QSKbtCbTATKb1rT3XGrlJLcKs3B11nsdHzRMoa3tMeKeUrFmIjmnQrkdI5UhHTVXzEzNCezpFSvUEwEqmEhARNpJSKoG4N7WmPlGPFVCJlzZHSNcRU7ApWj5QO7SkVGc3NzVRUVHSjR0oTKaeKmUQqMzMTt9vQ3OyOdFOUChkd2lMqupWVleF2uxk6dGhAr7PLH+iC484TM4mULhOjeoPgDO1pIqVUpOzfb02zyc7O7uLIluLj4mlqaqKuri4UzVI9EHOJlJZAULGqoaGBurq6IPRIWXOk3G7tvVUq3OxEatiwwG6QjPOst6fDe84Tc4mUlkBQsaqny8PY4uPjcbvd1NTUBKNZSqkA7Nu3Dwg8kYrX9fYcK2YSqWMLF2uPlIpNwUqkdL09pSLHroU1ePDggF4X7+mR0hIIzhMziZT2SKlYZ38TDcbQHuh6e0pFwv79+xk8eDAJCQkBvS5Oe6QcK2Yqm9s9UjrZXHUk2qs3B3NoD7RHSh0T7bERTfbv3x/wsB4ci1tNpMKv11Q2T0pKIi4+ThMp1aFor96sQ3sqVKI9NqJJdxOpuDidbB4pvaayOVgVmzWRUrEqeD1SOrSnVKTs27cv4NIHcKxHSudIOY8mUkpFCR3aUyq6GWMoKyvTHqkYo4mUUlHCTqTS0tJ6dJ64OBeJiYmaSCkVZg31DUDgpQ8AXC4XycnJmkg5kCZSSkWJqqoqUlJSAr7bpy0hPT1dh/aUCrOjR48CcNxxx3Xr9WlpaTq050CaSCkVJYKxPIwtPT1de6SUCrNaTyJ1wgkndOv1aWlp2iPlQDGVSCUmJNDc3OzN+pWKJSLCkCFDgnKujIwMTaSUCrOjR2tJTk7u1mRzsOZHaiLlPPGRbkAw2UMeBw8eZPjw4RFujVLBtWLFiqCdKz093bvml1IqPI4ePcqoUaNwubrXh6E9Us4UMwU54VgiVV5eHrRzqtihRQeP0aE95UtjIzyOHj3a7WE90DlSkdJrCnKCJlKqc1p08Bgd2lO+NDZCzxgTlERKe6TCr9cV5ARNpJTqSnp6OlVVVTQ3N0e6KUr1CvX19RhjGD16dLfPoXOknEkTKaV6ofT0dAAqKysj3BKleofa2lqg+3fswbGhPU/viHKImEqk4uMTEIGysrJIN0UpR8vIyAC0urlS4WInUuPGjev2OdLS0vTOdAfyK5ESkWQRWSkiG0TkExGZ28XxCSLytYgUBaWVfhIR4uMTtEdKhU20xEZrdo+UFuVUoRSt8REKtbW1xMcnkJWV1e1z2HXkdHjPWfwtf1AEiDFmmoiMBjaIyDhjzIEOjr8KGBSMBgYqIUETKRVWRURJbPiye6Q0kVIhVkQUxkcoHDlSS9++KT06h708VHV1ddBqyqme67JHSkRcwJXAowDGmH8Bm4FLOji+H3Ah8FLwmuk/TaRUuERbbPjKzMwErJprSoVCNMdHsBljqK09QkpKcBIpLYHgLP4M7R0HZAJf+mz7ApjcwfE/Ah4GInI7UEJCvCZSKlyiKjZ8DRw4ENBESoVU1MZHsJWWltLU1ETflL49Oo9vj5RyDn8SqcGeP31T4Era6X4VkSxghjHm+SC0rVu0R0qFUVTFhq8BAwYgIhorKpSiNj6Cbdu2bQCk9HBoT+dIOVMgd+21vt9S2jnm/wI/9edkdoVa+xGsAmsJCYmUl5fr7aHKy65Kaz8IfvXmqIgNX3FxcaSnp2uPVC8XhtiAIMZHOGIjFLyJlA7tRRV/48OfRKrU8+cAn20DfLYDICLHASONMe/600C7Qq39CF4ilUBDQwOHDx8OyvlU9LOr0toPgle9Oapio7WBAwdqItXLhTA2IATxEa7YCLbt27fjcrlISkrq0Xl0aC+8/I0Pf+7a2wFUAGM4FgDjgddaHTcTGCYi6zzPxwJ1IlIIXG6M2RXQO+gm36Kcqamp4bik6r2iKjZay8zM1KE9FUpRHR/BtGPHDpKTJ/b4PJpIOVOXPVLGGDfwJ2ApgIicAEwC/iYi40TkHRGJM8b81RhzqjGm0BhTCLwBrPQ8D1sgaHVzFS7RFhutZWZmao+UCploj49g2rFjB3369OnxeRISEujTp48mUg7j7xypIkBEZAPwFHCRMaYE6I/17SHBPlBEEj3fLM4ClojIC0FtcRc0kVJhVkSUxEZrAwcO1DhRoVZElMZHsBhj2LlzZ1ASKTi2TIxyDr8Kchpj6oAl7WzfAGS32tYAFAahbd1iJ1K6TIwKh2iKjda0R0qFWjTHR7CUlJRQW1tLn+TgJVLaI+UsMbXWHmiPlFL+yszM5OjRo941wJRSwbd9+3YAkvskB+V8/fv310TKYWIukYqPjyMuLk4TKaW6oEU5lQq9HTt2AOjQXgyLuUQKROd+KOUHXSZGqdDbsWMHLpeL5OTg9Ejp0J7zxGAiZX3T1jlSSnVOe6SUCr09e/aQnZ1tF3TsMR3ac56IJVJ2hdpQFFQbNGgQpaWlXR+oehW7Si2hqd4cNKGMDV92j5T23iqNjdDZt28f2dnZXR/oJ+2RCr+u4iNiiZRdoTYUATF48GAOHDgQ9POq6GZXqSW41ZuDLpSx4cvukdJESmlshE5xcTHDhw8P2vnsREqXQQufruIjJof2NJFSqmt2j5QOgysVOqHokWpubta7bR0kZhOpw4cP6y+aUp2Ij48nMzNTh8GVCpHq6moOHz4c1B6p/v37e8+tnCFmEylAe6WU6oLOJ1QqdIqLiwGC3iMFaAkEB9FESqleLCsrS4f2lAqRffv2AQR9jhRoj5STaCKlVC+mPVJKhU4oeqR0aM95NJFSqhfLysrSREqpELF7pIYNC15VCe2Rcp6YTKQGDRoEaCKlVFcGDRpERUUFTU1NkW6KUjGnuLiYgQMHBq2qOegcKSeKyUQqwB6i0wAAHtBJREFUKSmJAQMGaCKlVBfsLx1aS0qp4Nu/f39Qh/VAe6ScKCYrm4PWklJtafXmtrKysgCtJdXbaWyERklJCUOHDg3qOe1EqrKyMqjnVR3rlZXNQRMp1ZZWb27L7pHSeVK9m8ZGaJSUlHjn7AZLfHw8aWlpHDp0KKjnVR3rlZXNAYYMGaKJlFJd0ERKqdAwxlBSUsKQIUOCfu6MjAwqKiqCfl7VPTGbSA0dOpRvv/020s1QytF0aE+p0Dh06BCNjY2aSPUCMZtIDRs2jJqaGg4fPhzppijlWBkZGcTFxWnvrVJBZsdUKBKpzMxMDh48GPTzqu6J2UTKnuCnvVJKdczlcpGVlaWJlFJBVlJSAhD0OVKgPVJOE7OJlF0Abf9+R8+dVCridBhcqeCzEykd2ot9MZtIaY+UUv7RREqp4Avl0J6dSLnd7qCfWwUu5hMp7ZFSqnOaSCkVfCUlJSQmJjJgwICgnzsjIwO3261FOR0iZhOpAQMGkJycrB8QSnVh6NChlJaW0tzcHOmmKBUz7BpSnkKOQZWZmQmgw3sOEbOVzUWEoUOHao+U8tLqze0bOnQobrdba0n1YhobwReqGlJg9UiBJlLh0msrm9vX0B4pZdPqze3T+YRKYyP4Dhw4EJI79kATqXDrtZXNAe2RUsoPmkgpFXzh6JHSWlLOENOJ1LBhwzSRUqoLdiJl366tlOqZ5uZmSktLdWivl4jpRCo7O5uamhq9s0GpTtj/2WuPlFLBcfDgQdxutyZSvYRfiZSIJIvIShHZICKfiMjcDo5bICJvisi7IvKpiPxncJsbmJycHAD27t0byWaoGBatseErKSmJjIwMTaRU0MVCfHRHKKuaAyQkJJCamqqJlEPE+3lcESDGmGkiMhrYICLjjDGt15X4NbDIGLNZRIYAW0XkG2PMi0Fss99GjBgBwJ49e5gwYUIkmqBiXxFRGBut6TC4CpEiYiA+AhXKqua2jIwMysvLQ3Z+5b8ue6RExAVcCTwKYIz5F7AZuKSdw39njNnsOa4EWAu0+w0kHLRHSoVSNMdGayNGjNA4UUEVS/ERqHAkUoMGDaKsrCxk51f+82do7zggE/jSZ9sXwOTWBxpjft1qUzIQsX/poUOH4nK59ANChUrUxkZrmkipEIiZ+AiUvTxMqIb2wEqktPabM/iTSNm/CVU+2yqBQZ29SETSgCnAX7vXtJ6Lj49n2LBh7NmzJ1JNULEtamOjtZycHEpLSzl69Gikm6JiR8zER6BKSkpISUmhX79+IbvG4MGDNZFyiEDu2jOtnndV9/5+4F5jzDft7bQr1NqPUBVYy8nJ0W/avZxdldZ+EPzqzVEZG77s+YTFxcUhv5ZyjjDEBgQxPiIRG91h15AKxfIwNrtHylMoUoWAv/HhTyJlp7y+Ky8O8NnehohcBTQaY37X0TF2hVr7EaqA0CELZVeltR8Er3pzVMeGL51P2DuFMDYgBPERidjojgMHDoR0fhRYiVRjYyOVlZUhvU5v5m98+JNI7QAqgDE+28YDG9s7WETOB2YBN3ienxBAu4POTqQ0a1chENWx4cv3DlelgiRm4iNQ9oLFoTRokDVCqsN7kddlImWMcQN/ApaC95d7EvA3ERknIu+ISJxn3wzgP4H/DfQVkX7AnaFqvD9ycnKor6/XXzYVdNEeG76GDx8OaI+UCp5Yio9AhXJ5GJsmUs7h7xypIkBEZAPwFHCR5xbV/sBYIMFz3FNAIVAO1Hge3wliewOWm5sLwO7duyPZDBW7iojS2PCVlJTE4MGDtUdKBVsRMRAfgWhoaODgwYOaSPUifhXkNMbUAUva2b4ByPZ5nt36mEg77rjjANi5cydTp06NcGtUrInm2GgtJydHEykVVLEUH/6yC9vavbyhYg8daiIVeTG91h7AyJEjASuRUkp1LCcnh2++afdGQqWUn+w7X7OzQ5sbDhw4ENBEygliPpFKSUlh6NCh7NixI9JNUcrRjj/+eHbt2kVzc3Okm6JU1Nq3bx8Q+h6p+Ph4MjMzNZFygJhPpMAa3tMeKaU6N2rUKBoaGnTCuVI9EK4eKbDmSdlV1FXkaCKllAKsRApg+/btEW6JUtFr3759pKSk0L9//5BfSxMpZ4hYImVXqA1HQbXjjz+e4uJi6uvrQ34t5Vx2lVpCU705aMIZG740keq9NDaCZ9++fQwfPjykVc1t2dnZ3qFEFTpdxUfEEim7Qm04AuK4447DGKMTaXs5u0otwa3eHHThjA1f2dnZJCUlaSLVC2lsBE9xcXFYhvXAmoe1b98+3G53WK7XW3UVH71iaO/4448H4Ouvv45wS5RyLpfLxfHHH6+JlFI9sG/fvrAlUiNGjKChoYHy8vKwXE+1r1ckUmPHjgVg27ZtEW6JUs42atQoTaSU6ia32+0d2gsH+zq62Hhk9YpEKiMjg0GDBmkipVQXTjjhBLZv364lEJTqhrKyMpqamsLaIwW6tFOk9YpECmDcuHGaSCnVhQkTJlBfX69115TqBvvucLsQdKhpj5Qz9LpEyjNhTCnVjhNPPBGAzz77LMItUSr62MPi9h2woZaVlUVCQoImUhHWqxKpyspKrbmhVCfGjx+PiGgipVQ3fP3117hcrrD1SLlcLoYPH65DexHWqxIp0AnnSnUmJSWFUaNGaSKlVDds376d73znOyQmJobtmsOHD9ceqQjrNYnU+PHjAdi6dWuEW6KUs5100kmaSCnVDdu3bw/bsJ5Ne6Qir1dUNgerkFtWVhabNm0Ky/WU82j1Zv+ceOKJbN++ndra2ohcX4WfxkbPGWP4+uuvOeGEE8J63ZEjR7Jnzx4aGxvDet3eRCube4gIkydP1kSqF9Pqzf455ZRTMMbw6aefRuT6Kvw0NnquoqKCysrKsPdIjRkzhqamJl1PNoS0srmPU089la1bt+o3baU6MW3aNAA+/PDDCLdEqehhr5wRiUQK4KuvvgrrddUxvSqRmjx5Mm63m3/84x+RbopSjpWVlcUJJ5ygiZRSAdiyZQsAeXl5Yb2uJlKR16sSqVNPPRWAjRs3RrglSjlbfn4+H374odZdU8pPGzduJDMzk9zc3LBed8CAAQwaNEgTqQjqVYlUdnY2w4YN02/aSnUhPz+fsrIyXXdPKT9t3LiRKVOm2JOSw2r06NGaSEVQr0qkRITZs2ezdu1a/aatVCdmzJgBwNtvvx3hlijlfLW1tWzdupUpU6ZE5PpjxozRRCqCelUiBTBr1ixKS0v54osvIt0UpRxrzJgxHHfccbz66quRbopSjrd582bcbnfEEqlx48ZRVlZGaWlpRK7f2/W6RGr27NkAvPvuuxFuiVLOJSLMnz+fd955R+9yVaoL77//PkDEEqmpU6cCsGHDhohcv7frdYlUbm4uI0eO5H/+538i3RSlHG3evHnU1dXxzjvvRLopSjnaiy++yOT/3969h0lRnXkc//6GAQYYEDIgBLJGUSGIJKBA1NHARiBkFSRkkZVklHiBLBs0YjRqIg67asz6YNzHxBCSrAqCaAQk8RYwiCIXgcTbhnhN1ESIIpdhRgWH4d0/qnpSjD3Q3fRMdU+/n+fpp7uqq855e6bf7lOnTp0ePJgePXrEUv/JJ59McXEx69ati6X+QlcwM5tHjR8/nuXLl7Njx45mr9vFx2dvTs/w4cPp2rUrd999d6xxuKbnuZG5rVu3sn79esaNGxdbDO3ateOkk07yC6maiM9snsSkSZOora1l8eLFzV63i4/P3pyeNm3acMEFF7Bs2TLeeeedWGNxTctzI3PLli0DiLUhBcGVths3bvSfimkCPrN5EoMGDaJPnz7Mnz8/7lCcy2mXXHIJ+/btY+7cuXGH4lzO2b9/P7fffjsDBgzghBNOiDWW8vJyPvzwQx8nFYOCbEhJYurUqaxevdq7Qp07iL59+zJmzBhmz57N9u3b4w7HuZyydOlSNm/ezLXXXhvL/FFRo0ePpn379ixYsCDWOApRQTakAKZOnUrXrl2ZNWuWzynl3EH84Ac/oLq6mpkzZ8YdinM5Y8eOHVxxxRX07duXCRMmxB0OpaWljBs3jvvvv5+9e/fGHU5BSakhJalE0l2S1kvaJGnUQba9QtLvw9uV2Qs1uzp06MA111zD8uXLWbhwYdzhuDzVEnOjof79+3PppZdyxx13sGTJkrjDcXmkpeZHTU0NEydOZMuWLcybN49WrVrFHRIAFRUV7Ny5k0WLFsUdSkFJtUeqEpCZnQJMAhZJ6t5wI0mjgUuA8vB2kaSzshRr1l122WWUl5czbdo0NmzYEHc4Lj9V0gJzo6Gbb76ZIUOGMGnSJJYuXRp3OC5/VNLC8mP9+vWUl5ezcuVK5s6dy9ChQ+MOqd6oUaMYOnQoV111FTt37ow7nIJxyIaUpCLgYuCXAGb2CvAs8PUkm08FFprZHjPbAywAvpm9cLOrVatWLFy4kLKyMr74xS/y4x//mI8++ijusFyeaMm50VDbtm155JFH+NznPsf48eOpqKjg+eef99PirlEtJT92797N2rVruemmmzj99NM59dRTee+993jooYeYPHly3OEdoKioiDlz5rB9+3bGjBnDtm3b4g6pIBSnsE1voAx4KbJuMzA4ybZDgHsbbDct4+iawVFHHcXTTz/N+eefz/Tp06msrKS8vJy+ffvSq1cvSkpKKCkpoW3bthkPJsxkv7gHLua7fv36ceKJJzZ1NS06Nxrq2rUrq1at4oYbbmD27Nncc889HHnkkZxyyin07t2bbt260alTJ1q3bl1/Ky4upqgo+0MxmyI/CiXnunTpwogRI5qjqmbLj7/85S+sXbuW/fv3U1dXl/attraWmpoa3n//fWpqaqiurmbr1q289dZbVFVV1dczYMAAbrnlFqZOnUrHjh3T+mM0l0GDBnHvvfdSUVHBMcccwznnnEP//v3p0qULHTt2rM/HxPv9UPeFZuzYsbRt2zatfVJpSCW6Yasi63YBya717J5kuyOTFZqYWC3h+uuvj21ukJ49e7JixQpWrFjBvHnzeO6553jssce8dyqPXXfddfUNqcrKSmbNmhV9OluTDrb43GioXbt23HjjjcyYMYPFixezdu1aNmzYwMqVK6mpqYk7PJeCwYMH1zekmjA3oAnyo7HcWLNmDRUVFRkFWVRUROvWrSktLaW0tJQOHTpQWlrKMcccw7BhwzjqqKM47rjjOOOMM+jatWtGdTS3CRMm0K9fP2bPnu3jgNO0bdu2+oZUqvmRSkMqoWEffmPN1ZT6+nv27MmWLbkz95skRo0axahRwVjIuro6du3axd69e9m7dy979uzJqNxMTn346ZLD161bt/rHlZWVBzREJGX7jdeicyOZsrIypkyZwpQpU+rX7dmzh+rqampraw+4ZVtT5Ech5VxJSUn942bIDchifjSWG2PHjuXll1+mqKiIVq1apXwrKipqsT0vJ554InfeeScQ5GZVVRW7d+/GzOrf74e6L0SdO3euf5xqfqTSkEr8nHTnRh433LZzZLkzkJcnaVu1akVZWVncYbjcVpC50ZjEaXDnQs2WH506daJTp06ZxFgQErnZvfvHxvm7LEhl8MLrwA6gb2TdCcDGJNtuTHE751oCzw3nGuf54QrCIRtSZrYf+DlwIYCk44GBwAJJ/ST9TlJiEo05wHnh3CElBJe7zmma0J2Ll+eGc43z/HCFItUxUpXAHEnrw33OM7O/Szoa+AzQGqgzs8ck9QfWhPv9r5k9nN2QncsplXhuONeYSjw/XAuX0nXJ4dwek83sFDMbbGbLw/XrzaxXOO9HYtvZZnZyeLslk6By5QqlTOV7/JD/r6G54vfcSE++xw/5/xqaM/7mzo/mlO/vg0Px15e6nPytvQaXG+adfI8f8v815Hv8jcn315Xv8UP+v4Z8jz9XtPS/o7++1OVkQypVqbQoU211Znu75q4z3+Nviu0K2ZtvvHHIbXL5/+K5kf52hZ4Xmb7+5v67ZVrfGynkdDbry1Rz/x9y4X2vuOaLGDx4sG3atCnpc5IY8b3fALDihrMPeG7fj74LQPHlP0TSAfNdjHjwagAeH3fzAWWl8hqzuZ3XGf92B9tG0u/NLNnsyjnhULmReF0jv/8QkDxHnly1ijOXPfOx/aM5kmv/F68z/jpbSm6ko7n2S+Tr4zeOSbpf9LutoREPXs2qVavYd9v6Jo/T9zvoPknzI7aGlKRtwJuNPN0TSGViuFS2y2ZZXmf8daa63cG2+bSZdWvkudh5bnidMdZZCLnh+/l+me6XND9ia0g555xzzuW7vB4j5ZxzzjkXJ29IOeecc85lKKcaUuGstndJWi9pk6RRcceULklDJL0maXLcsaRLUmtJ35a0StKTktZJOjPuuNIhaXo4Y/JKSS9Kmh53TNnguREvzw0XJeliSWskPS3pBUnD4o4pW8JZ558IX9uzkr4bd0yHq7HPHkkDw1xeI+k3kjL6gd1UZzZvLpUE47ZOkdQHWC+pn5m9E3NcKZH0FWACUBV3LBnqBVwGDDSzKkkjgWWS+prZ2zHHlqqLgZFm9m74kxR/kvSMmW2IO7DDVInnRpw8NxwAkiYAZwJfMLM6Sd8AesQcVjbdDawws++FDYtXJT1vZo/FHVgmGvvskdQGWAZMNrMnJM0i+FmiCenWkTM9UpKKCBL9lwBm9grwLPD1OONK00YzmwRUxx1IhqqBmWZWBWBmK4A9wGmxRpWeCjN7F8DMXgV2AkfHGtFh8tzICZ4bLmEm8J9mVgdgZnea2X0xx5RN/YF1AGa2HXgVGBRrRIensc+eLxP8PNET4fIvgPGS0r5qNWcaUkBvoAx4KbJuM5Czc5o0ZGZ/izuGw2Fm281sfmJZkoA2wLb4okqPmb2QeCzpq0ANsDy+iLLCcyNmnhsOQNKRQD9gYHj6a7WkqXHHlWUPA2MAJPUmaFh9fFK6PHGQz54hRD5TzeyvwAfASenWkUun9rqH99Hut13ACTHE4gLDCOZseSruQNIhaQBwH9ABONfMdsUc0uHy3Mg9nhuF6WhAwFeAEcCRwAZJVWa2KM7Asugi4NeSXgc+Acwws5Uxx9QUuvPxoQa7CP6nacmlHqmEhhNbKZYoCpykEuAmgvPH++OOJx1m9qKZnQCMIxjHMjTumLLEcyMHeG4UtLYE35u3m1mdmW0F5gMXxhtWVi0B1pnZscBngRmS0u6lyRPJJtJM+3M1lxpS74b3nSPrOkfWu2YSnraYC/zIzH4fdzyZMrNnCbqpL407lsPkuZEjPDcK3s7wPnqRx9+AT8UQS9ZJ6kfQ03Yb1J/uWgFcGWdcTeRdDvxMhQw/V3OpIfU6sAPoG1l3ArAxnnAK2mxgg5n9SlJbSUfFHVAqJJWFYz+i3ic4jZHPPDdyh+dGYXuVYBxN9PRPNzL7iZJc1Ca8r42sqwU6xRBLU9tI5DNV0j8B7YG0D5BypiEVdpH/nLCLNLw8dyCwIM64Ck04Z0gxcJekUuBY8qfbuiNwnaT2AOHVF+OA38Ua1WHy3MgNnhvOzPYC8wjGESGpAzCRYMqAluAl4G3gPABJHYGxQEscI/UoUByZA+xCYImZpX0BSU791l449mAO8BmCD6xrzSxvriqRdDLBEetA4O/AZjMbH29UqQvnJ3o5yVOzzKyymcNJW/j++R5B1/RHBN20SwkuVc6rsSwNeW7Ey3PDJYSNpzkEvcL7gMXALZZLX6aHQdJg4FaCsUKlwOPANWa2L9bAMnSwzx5Jg4A7gDqC07aTwykf0qujhfzvnXPOOeeaXc6c2nPOOeecyzfekHLOOeecy5A3pJxzzjnnMuQNKeecc865DHlDyjnnnHMuQ96Qcs4555zLkDeknHPOuSYgaaSk5ySZpCclfSLL5feUNCvJ+hWSfprNusJyJenNcKLOdPf9hqRTsx1TLvB5pJxzzrkmImk48ATQOpuTWoaTrD4JTDKz1yPruwGbgf1ATzOry2Kdg4BbzeyfM4z3ceBCM3slWzHlAu+RaiKSviZpl6Q2h9660TKOkDQ7m3GlUOfRkq6W5O8N1yQ8N5zLiunAH6KNqNBXgSuAI4DhWa5zNMFPq6TNzPYAtwM/ympEOcA/EJrOOUBr4EuZ7Bz+yvz9BD8/0GzM7I3w4bXNWa8rKJ4bzgGShoSn/J4K74dEnjtW0mpJ6yQtkLRE0huSpoSbnAs8laTYEcAiggbPxAb1dZP0iKQNkpZK+oWkv0uaGalzeRjLakmnNSj7y8Bjkh4PT1c+FO5XIWmLpAckzQzLvF3SQkkvS7or3P9JYGS2T3HGzsz8luUbwZHAvQRv5vkZljEBeDSm+EuAdwm6hWP/e/qt5dw8N/xWaDeCXiEDihusPwJ4DxgeLp8RLncOl58h+I07gF5AFVAZLguoBT7foMzuwB3h43OBbdF6CQ5AfhY+7gT8GbgrXC4G/kRw6g3gs2E8HSPbvxE+bh0+d1qk7N/wj+FCdwEvAu3D1/n9yHY1wLC4/y/ZvHmPVNM4h+Bo+V7gHEltMyjjgKMNSVPCo5FFkn4m6f8k3SPpeEm/kvSapGmR7WdKWhneHpLUM/LcYEmbJK2R9FNJT0t6SdJYqO+C/QNBF7Fz2eS54VzgbGC3ma0CMLPVBD+cO1bSp4GhwD3hc28T9OYklBE0fKoblPlVgsYSBA2bEuBMAEmtgK9EytwNPBTZ9/PAscD88PkXgLfDOCHo6fpd+FwtwcHQ+WHZg4AXLGwphR43sw/MrMrMboisryZo8LUYxXEH0EJ9CbiIYLBfHUF36INplnES8EBiwczmhh/4FwEnAh8CfwU+IPhiGQQ8KWmuBQMadwJnmplJmgz8EKgIx6UsBa4ys3slDQQ2AReb2a8j9b8WxuBcNnluOBf4FEGPUdS2cP0nw+X3Is/tiDxWeN/warGxwERJifWJHPgt0I3gO79hmZ0i8RiwIjh7DkBbgh4lCE/rRfadR3Ca7zKgApjbIJYqkjNa2LAib0hlmaTOQE145IqkxQTnqR+UNAIYZWZXhc+VEHTflptZTYOiuvPxow2ADWZWFe7/KvBi+IXwAlAKHAlsIfgieSIcGNsJSAzsPTXc5n4AM3tO0uYk9VQDvTP5GziXjOeGcwFJRwB/I2jcRHUL12+NLL8VPi4D3ggfbyc4EKmfhkBSD+A1M/tWZN1E4KeSvknQSNsXlvmnSJkJfwVqzWx4ZP8OBAc9EPRsXZl4zsw2SHqHoJerj5m9lNqrpyPwTorb5oUW1SrMEeOAUyWtkrQKGAKcLakd8DLQJ7LtJcA9Sb4oIDjiSDY3RfQLZF9i2f5xWW0bSccTfBlcaWZfAL5NcK4agiOdXXbgJbHRI52EFnfU4GLnueFcoAvBabiOkr4AIKk8XP9rM3sT2EDQ04OkXkB5Ymcz2w/8Efh0pMx/BR5uUM+jQAdgZPi+XhIpsxPBVXgJzwBvSRofPl9M0FvcR1J/YIuZ7WpQ/nzgVmB5Ki9aUhnQLoy9xfAeqez7EnByeA6Z8HTBNuAsgrEhPSPrJxMMMEzmXSJHG2kaRHDufWO43Dry3Fags6TiyBdMGR/X4o4aXOw8N1xBkXQGkJgw877IKbf2BO/j0cDssHfUgC9HGiuTgLslnQ28RDAHU/QA4gGCxtWvwp6n7wN/lrQ8cjBwI8FBxR2SLgW+FZa5kaB3awVBwwYzq5M0BvhJeLquCLjTzJ6X9B0OPK2XcA9wPcF4qcRrnhG+rj2S2id6mUNnAE+Z2bup/P3yhTekskhSF2Bf4osCwMw+kvQIcK6ZPSCpLkyaC4GFZvZBI8W9wIFHG+l4DegiqY8FE59FjzrWESTwRGBBOA7k+CRlHA2syrB+5w7gueEKUTiAfNghNhveyPqdZnZ6YkHSwwTvz4TbCMb+fdLM7gPuS1L/dIL5phJldAbODnu0kPSTaJkWzEk1umE54bqrk6zfDTwcbRiZ2a0EvVQHCHu4LgMuT/pq85h3T2dJeM77CaBc0lmR9WcRXH0xVtIdBEcBxwPfAA42hX/iaCNRziSCo/TRkv5dwbwfA4GrJZ0qaUm46aKwjpuA5ZKWEYwP6SFpnpl9BIwHZkh6GriAoAvZInW1DmNu1nl6XMvkueFcRm6T9BkASf8EnEZ41RyAmVUTNPrTaZh8BxgZltmZ4Iq83x5qJzMbYWabEsuS/kVSV+A8Ir1RhzAZ+G8zey6NePOC/0RMM5N0PTAAWGNmjc7wGl6qugqYZmYvZjmGT5jZjsjyH4HvmNmj4fI04OgGXbLONSnPDef+ITxAuBx4n6DBf0vY83Q4ZY4E/gvYQ3CKep6Z/U8G5UwH/gN4BRhvWfzpm3zkDalmJunfCKbI721mHx5i2x4EE7B9M8sxPAKcb2bvSTqZ4Dz5cWa2Q8H8JdMIJlCrPWhBzmWR54ZzLh95Q6oASboc+BrBkU5b4FozWxlvVM7Fz3PDOZcub0g555xzzmXIB5s755xzzmXIG1LOOeeccxnyhpRzzjnnXIa8IeWcc845lyFvSDnnnHPOZcgbUs4555xzGfKGlHPOOedchrwh5ZxzzjmXof8HYwThGy4rm5wAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 720x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "O8-QDCMD2d_w",
        "outputId": "82e6f155-e093-4d6c-ad97-6f20d6757caa",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        }
      },
      "source": [
        "c[0:2]"
      ],
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            },
            "text/plain": [
              "'27'"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 15
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 311
        },
        "id": "sk1-JeWaZ_eO",
        "outputId": "e7681021-2b1c-4bf5-bd6d-e7e26f789ebc"
      },
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "from matplotlib import rcParams\n",
        "### THIS NOTEBOOK ONLY PLOTS CLUSTERS 2732 and 816.\n",
        "\n",
        "#############################\n",
        "#### LOAD DATA   ############\n",
        "#############################\n",
        "PDF = pd.read_csv(\"https://raw.githubusercontent.com/Rogerod23/Rogelio/main/PDF.csv\")\n",
        "percentiles = pd.read_csv(\"https://raw.githubusercontent.com/Rogerod23/Rogelio/main/percentiles_all.csv\")\n",
        "\n",
        "del PDF[\"Unnamed: 0\"]\n",
        "del percentiles[\"Unnamed: 0\"]\n",
        "\n",
        "\n",
        "cluster_list = [\"11-E\", \"90-E\", \"93-E\", \"383-W\", \"417-W\", \"451-W\", \"531-E\", \"534-E\",\n",
        "                \"589-E\", \"816-E\", \"894-E\", \"1252-W\", \"1381-W\", \"1564-W\", \"1576-W\",\n",
        "                \"1949-W\", \"2449-W\", \"2732-W\", \"2740-W\"]\n",
        "\n",
        "PDF_columns = PDF.columns\n",
        "percentiles_columns = percentiles.columns\n",
        "\n",
        "# Getting name of the age_pdf of cb2018\n",
        "a = np.where(PDF_columns.str.contains(\"tpdf_cb2018\") == True)\n",
        "age_CB2018 = PDF_columns[a]\n",
        "\n",
        "\n",
        "a = np.where(percentiles_columns.str.contains(\"cb2018\") == True)\n",
        "percentiles_CB2018 = percentiles_columns[a]\n",
        "\n",
        "a = np.where(percentiles_CB2018.str.contains(\"t_\") == True)\n",
        "percentiles_CB2018 = percentiles_CB2018[a]\n",
        "\n",
        "\n",
        "\n",
        "# Getting name of the age_pdf of Krumholz2015\n",
        "a = np.where(PDF_columns.str.contains(\"tpdf_modp020\") == True)\n",
        "age_Krumholz2015 = PDF_columns[a]\n",
        "\n",
        "a = np.where(percentiles_columns.str.contains(\"modp020\") == True)\n",
        "percentiles_Krumholz2015 = percentiles_columns[a]\n",
        "\n",
        "a = np.where(percentiles_Krumholz2015.str.contains(\"t_\") == True)\n",
        "percentiles_Krumholz2015 = percentiles_Krumholz2015[a]\n",
        "\n",
        "\n",
        "\n",
        "#PDF COLUMN INFO TO SEE WHAT PLOT YOU WANT TO SEE\n",
        "#print(percentiles_columns)\n",
        "#print(PDF_columns)\n",
        "\n",
        "#########################################################\n",
        "# Colors \n",
        "#########################################################\n",
        "\n",
        "alpha=1\n",
        "blue=(7/255,80/255,155/255,alpha)\n",
        "lblue=(106/255,175/255,210/255,alpha)\n",
        "pistache=(190/255,250/255,215/255,alpha)\n",
        "yellow=(248/255,250/255,3/255,alpha)\n",
        "red=(251/255,105/255,73/255,alpha)\n",
        "wine=(160/255,55/255,98/255,alpha)\n",
        "dgreen=(5/255,149/255,74/255,alpha)\n",
        "gray=(203/255,202/255,201/255,alpha)\n",
        "black = (0,0,0,1)\n",
        "\n",
        "mycolor=[blue,red,dgreen]\n",
        "\n",
        "#########################################################\n",
        "# Font parameters \n",
        "#########################################################\n",
        "fontsize=12\n",
        "font = {'family' : 'serif',\n",
        "        'weight' : 'regular',\n",
        "        'size'   : fontsize}\n",
        "plt.rc('font', **font)\n",
        "\n",
        "mycolor=[blue,red,dgreen]\n",
        "\n",
        "mark = \"None\"\n",
        "lw_ = 1.5\n",
        "linestyle_ = \"-\"\n",
        "\n",
        "#Figure size and subplots\n",
        "i = 0\n",
        "for CB2018, K2015 in zip(age_CB2018, age_Krumholz2015):\n",
        "#get cluster ID\n",
        "  if len(CB2018) == 25:\n",
        "    cid = CB2018[22:25]\n",
        "  else:\n",
        "    cid = CB2018[22:26]\n",
        "#Get percentiles column name\n",
        "  percentile_16th_CB2018 = \"t_16th_cb2018_z014_u3_{}\".format(cid)\n",
        "  percentile_50th_CB2018 = \"t_50th_cb2018_z014_u3_{}\".format(cid)\n",
        "  precentile_84th_CB2018 = \"t_84th_cb2018_z014_u3_{}\".format(cid)\n",
        "\n",
        "\n",
        "  percentile_16th_KRUM2015 = \"t_16th_modp020_chabrier_MW_phi0.50_{}\".format(cid)\n",
        "  percentile_50th_KRUM2015 = \"t_50th_modp020_chabrier_MW_phi0.50_{}\".format(cid)\n",
        "  precentile_84th_KRUM2015 = \"t_84th_modp020_chabrier_MW_phi0.50_{}\".format(cid)\n",
        "  fig, axes = plt.subplots(1,1,figsize=(4,4))\n",
        "\n",
        "  plt.subplots_adjust(hspace=0.5)\n",
        "  plt.subplots_adjust(wspace=0.25)\n",
        "  code=['SLUG']\n",
        "#figure title\n",
        "  for c in cluster_list:\n",
        "    if cid == c[0:2]:\n",
        "      fig.suptitle('Cluster {}'.format(c), size = 12)\n",
        "    elif cid == c[0:3]:\n",
        "      fig.suptitle('Cluster {}'.format(c), size = 12)\n",
        "    elif cid == c[0:4]:\n",
        "      fig.suptitle('Cluster {}'.format(c), size = 12)\n",
        "\n",
        "\n",
        "  for icode in range(1):\n",
        "      axes.set_ylabel('')\n",
        "      axes.minorticks_on() #...add minor ticks to x and y axes\n",
        "      axes.tick_params(which='minor', length = 2, width = 1, #...format minor ticks\n",
        "                                      direction = 'in', top = True, right = True)\n",
        "      axes.tick_params(which='major', length = 4, width = 1, #...format minor ticks  \n",
        "                                        direction = 'in', top = True, right = True)\n",
        "  #    if icode==0 and icluster == 0:\n",
        "  #        axes[icluster][icode].annotate('Cluster 816-E', \n",
        "  #                                        xy = (0.4,0.95), xytext = (0.4,0.95), \n",
        "  #                                        xycoords='figure fraction', size = 12)\n",
        "  #    if icode==0 and icluster == 1:\n",
        "  #        axes[icluster][icode].annotate('Cluster 2732-W', \n",
        "  #                                        xy = (0.4,0.47), xytext = (0.4,0.47), \n",
        "  #                                        xycoords='figure fraction', size = 12)            \n",
        "      axes.set_xlabel('Log(T) (yr)')\n",
        "      \n",
        "  axes.set_title(code[0],fontsize=12)\n",
        "\n",
        "\n",
        "  #Extinction plots\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "  axes.plot(PDF.loc[:, PDF_columns[6]], PDF.loc[:, K2015], color = black, marker = mark, \n",
        "          markevery= 10, linestyle = linestyle_, lw=lw_)\n",
        "  \n",
        "\n",
        "\n",
        "\n",
        "\n",
        "  axes.axvline(percentiles.loc[0, percentile_16th_KRUM2015], \n",
        "                   color =mycolor[0], linestyle = '-', linewidth = lw_,label='16th percentile')\n",
        "  axes.axvline(percentiles.loc[0, percentile_50th_KRUM2015], color =mycolor[1],\n",
        "                   linestyle =  '-', linewidth = lw_,label='50th percentile')\n",
        "  axes.axvline(percentiles.loc[0, precentile_84th_KRUM2015], color =mycolor[2],\n",
        "                   linestyle =  '-', linewidth = lw_,label='84th percentile')\n",
        "  \n",
        "\n",
        "  axes.legend(loc = \"upper right\")\n",
        "\n",
        "\n",
        "  if len(CB2018) == 25:\n",
        "\n",
        "    plt.savefig(\"age_pdf_slug_{}.png\".format(CB2018[20:25]), bbox_inches = 'tight',dpi=100)\n",
        "  else:\n",
        "    plt.savefig(\"age_pdf_slug_{}.png\".format(CB2018[20:26]), bbox_inches = 'tight',dpi=100)\n",
        "  if i == 0:\n",
        "    plt.show()\n",
        "  plt.close()\n",
        "  i += 1\n",
        "#axes[0][0].axvline(percentiles.loc[0, percentiles_columns[6]], \n",
        "#                   color =mycolor[0], linestyle = '-', linewidth = lw_,label='16th percentile')\n",
        "#axes[0][0].axvline(percentiles.loc[0, percentiles_columns[7]], color =mycolor[1],\n",
        "#                   linestyle =  '-', linewidth = lw_,label='50th percentile')\n",
        "#axes[0][0].axvline(percentiles.loc[0, percentiles_columns[8]], color =mycolor[2],\n",
        "#                   linestyle =  '-', linewidth = lw_,label='84th percentile')\n",
        "#axes[0][0].legend()\n",
        "\n",
        "\n",
        "#plt.savefig(\"fig11_age_pdf.png\", bbox_inches = 'tight',dpi=100)\n",
        "#plt.show()\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAEmCAYAAAB1Z7W6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2deXxU1fn/308WCGuAAJGwi8oiyC5bIAGVqiCLVesuil/1q7buS3+1Aq1tXQrVr5WqrYoLLsUi7isSFRWECFURWQMIEUiCQAKEbM/vj3vvcDOZyUySSeYmOe/X674y95xzz3kmyXzmrM8jqorBYDCEIibaBhgMhvqBEQuDwRAWRiwMBkNYGLEwGAxhYcTCYDCEhRELg8EQFkYsPI6IpInIhyKSISKf2a+vEJEmInKNiGwTkQXRtjMQIjJcRF4TkQ9EZI2IvCIi7Vz5fUVkmYgst/PvcuUdJyJPiMhSEfnCfu+DArSRKiKbRCRfRF52pd8nIotc97eKyG4RebM233ODRlXN5dELOBPYDpzgSjsLKAYG2fezgQURaGsbkB5h++8DrrZfxwIfA/Nd+V8Bf7JfJwH7gDPt+9OB511l5wDfB2nnDOAgEOdK+xTY41fupWj/TevzZXoWHkVEYoD5WB+mzU66qr4LvBo1w6rGs8BCAFUtBTKBHq78k4Ev7fw8YBMw2M5bA9ztKvsV0D1IO58DCcBQABFphiU8SSLSz05rARTU9A01ZoxYeJchQE/ggwB5VwHf+ieKyIMisl9EZtj3T4hIoYik2/fJIvKuq+t/l53+DHAc8LA93HE+dJeLyAoR+UREXhSR1nb6fLudP4rIYhHZISKz/e1R1U2qesR+phMwCXjKVeRt4Bw7/3gs8VhpP5unqrvsvGbA5cA/A/2iVPUwsApIs5NGAu8DX7vSRmOJiqGaGLHwLr3sn7v8M1T1iP1N7Z9+J7DWdX8tsNtV5DYgQ1XHA7/A/qCq6pV2uZtVNV1VM0VkDDAPOEdV02w75tnlr7fbGQGcjzUM2BnsjdhzBxuxhgH/cWXNBE4SkS1YvY5bVfVjv2fvBH4E4oH/F6wN4BMg3X6dBmTYaWmutE8qed4QAiMWjYt9wFkicrKqHgImVlJ2BvCmqubY9y8Cl4iIuMq8raqlqrpBVf8VrCJVPR/oBIwSkb+6shYDX6pqL+AU4FYRGeL37INAB6ye1MciEhukmQxgjJ3fR1XXU14seqhqViXv1xACIxbexZmn6BzBOh8C/gO8IiJrsYYFwegCTLCHJRnAo8AerIlIhwPhNqyqBcBfgJtFJEFE+mJNYj5s5/8IfAjcEeBZBf6ANZ8xXkQGOXaJyMN2sc+B5sAo4KidthzoaK+iHA7XVkNgjFh4lzVYE34Vvv1F5Cn/b2AXRUBT130b1+uOqvqoqvYHbgeeF5FeBOZH4C17WJKuqqnAMFXNDfcNiMg1fkmHsFZFEoAmdlqxK78YcOZFJotIipNhD7uOAi1Uda3Lrpvt/EPAauAujs177Af+a6etCNduQ2CMWHgUVS0DrgF+KyInOOkiciUwANfchB9ZQH+7bBrWt63DX1x7FVZiCYszrMgHmovIeBG5CVgATBKRtnZdvYGq7lG4yJ77wB6+XA1k2h/iH7DmQS6y81sBU7CWVwGGYQ2FnPf9K9vWlZW09wkwmfJzE58AF2DmK2pOtNduzVX5BYzF+gB9itWtfgpob+ddg7U/YjfwOzutD9b4/hOs3sM2LGEZijXs+NSuLxP4jaudG4Hvsb6BT7bTLrXvPwbeA06y0x8EnA/8rZXYPtl+9mOs1YolwPGu/GG2PZ9h9aQewt4rgSV4i1zv+wtgQojf1ZnAXr+0acCP0f47NoRL7F+owWAwVIoZhhgMhrAwYmEwGMLCiIWhHCLSXUReF5FP7UNcGSJyg503wN7RqSKy0tlKbecdZ5dVEfmviEyw09vZ6ftFZKGrfJKIPCYin9v5X4nIk5Ws8hiijJmzMJRDRJYBi1R1vn2fDjyqqgPs+x5YKy4nquvMiut5Bc5Q1Y/80jNUNd1+3RLrrMfjqvp/dlpz4CWgWFXPq5U3Z6gRpmdh8OdUrN2QAKhqBvZhsAhyB7DfEQq7ncPAnVh7MQwexIiFwZ/twJ32KU0AVPX+CLfxSwIckFNr2/gVEW7LECGMWBj8+TXWAbNdIvKMvbEr0vQiwAE5g7cxYmEoh6ouBbphbejqASwTkSciUbXfvZS7EVlkT55ui0BbhlrAiIWhAqp6SFX/pdZR9vHA1ba/iXAo4ti5DwBEpCnHDncBbMXvgJxaJ1PvIbiDG0OUMWJhKIeI/MN9r6qfAHlAYojnLrZf/ojltMdNTzvd4T9YW7MN9QgjFgZ/TheRU50be86iDOscSGU4J0xfw+qJtLOfj8FyuvO6q+xDQCsRuc3VThzWORizlu9R4qJtgMFzPAg8JCJlWMfJy4BpqnpERPpge8sCHhWRQMucs7COmX8lInvttNdU9S2ngKrmi8hY4I8ishLraHobLDd4A2vlXRlqjNmUZTAYwsIMQwwGQ1gYsTAYDGFhxMJgMISFEQuDwRAWRiwMBkNYeHLptH379tqjRw8AsrOzSUlJCViusrz6nu/kbdxleds/qXNiwHw3useK8yPJXSqte+P+XRQdPUr/5OCbMiNhe23UXdv5XrYtVH6k6s7MzMxV1Q7++Z5cOh02bJiuXr0aABEhmI2V5dX3fCfvjHus7Qkf3jc55LMlf7OCkMfd8kCldZ++5G4yMjIoeTi4d/xI2F4bddd2vpdtC5UfqbpFJFNVh/nne7JnkZ2djYgwa9YsZs2aVe16Qj1b2/mhqOx5L9sWKt/LtoXK97Jt4eTX5Nm0tDTsgHMBux9h9yxEZDiWJ6P7VHVBkDILKB8lGyAVOE5Vc4PkT1YrWpUPd88ihE2VKmW0iYR9wXoWgXD3LCojnJ5FNPHy37Ux2FajnoWITMcKgBsqXN1hx3Wa/dxgYJ66oli582tKTb8Fahsv2+fMCXkRL//eGrNt4a6GrFLVi7GiVgVFrejabi4BXqiOYeEwe/bs2qo6InjZPi+LhZd/b43ZtrB6Fqq6s6oV26cNpwF/9Et/HDgZKAAeVNVlVa3bEJri4mJ27txJYWFhhbx7uk4BYP369XVtlsEDxMbG0qZNG9q3b09MTPi7J2pzgjMdK66le+iyHvhYVVfZcyDLRCRVVYPF7TRUk507d9KqVSt69OjhTFr5iPnZ0v7ebbtEwzRDFFFViouL2bNnDzt37qRbt25hP1ubm7IuBZ53J6jqA6q6yn69CngbuNb/QWc1xLm83PXzKoWFhSQlJVUQCkPjRkRo0qQJnTt35tAhy8PA7Nmzy33eCLIaUitiISIJwDisYLqVsYMAbtRSUlLKBWQ1YhE+WlbGgAEDyM/PN0JhCIp7+DF79uxynzcgO+AztWTLZOBdVS1xJ4rInX7lkoMZZqgeOTk5fPfdd75vDYMhUtRILESkvYh85rhQc3EJfkMQm1tFpKP9bE9gapByhmqyK9vS3qNHj1JcXBxla4KzatUqTjjhBBYsWFAh78CBA1x99dWkpqYycOBAzjzzTEpKSipWAjz55JP06NGDGTNm1K7BHmDbtm3062dFjPzqq68YNGhQna5qhSUWIjJURDKAQcDdIrLYzmoG9AGau8q2BU5Q1a8CVPVX4DUR+QRYBNxoO4Q1RIBDBQUcPHCAiy+2fOcePHgwyhYF5rXXXuNvf/sbiYmBfQBfcsklTJ8+neXLl7NmzRoSEhIoKysDrC6zWxiuueaaBisUPXr0ICMjo9z9F198AcCpp57Kww8/XKf2hLt0mom1uuGf/iPQwS/tZ2BAkHr+iiUYhlpgz549xMTE8H//939kZmayf/9+kpKSom1WBYYPH8706dNJT0+vkJeZmUlOTg6TJk0CrLH1kiVL6thC79KmTZuote3JI+rOaoiZ2Kwahw8fplnz5iQlJdG8eXMOHjzo+0b2El26BF+yXbp0KT169ODXv/41o0eP5pxzzuG7774D4JVXXmHBggW89957pKen86c//cn3XGFhIddeey1Dhgzh7LPPDri/pKioiPT0dESEv/zlL0ycOJGTTz6Z++8/Fp2xuLiYO+64g9GjR5OamsqcOXNQ1XLPPvbYY0yaNIlWrVqRkZFBQUEB//M//0NqaiqpqanMnDmTnJwcALZs2cLEiRNJS0tj7Nixvp7BG2+8QZ8+fUhLS+Ouu+5i5MiRjBkzhr17LR/HV155Jbt37+bmm28mPT2dzMxMTjvtNESEbdu2BfzdFRQUcNVVV5Gamsro0aN5/PHHq/R3cVZFCLIa4smDZCkpKWRnm3nPqlJ49CgJCQkANG3alNLSUo4ePUqzZs24/akv+W9WHgBHSqx4P83iIru9ZWDPJP46c1SN6ti2bRuLFy9m0aJFPProo7zwwguMHz+erVu38qtf/Yr169ezbdu2CnMdn332Gd999x2JiYkMHDiQ1157jYsuuqhcmSZNmpCRkYGIsH//fj744AP27dvHySefzJAhQ5g4cSIPPvggmZmZfPbZZ5SVlZGWlkavXr249NJLfc/m5eXx9ttv89JLL9G6dWtuvfVWSktLWb58OWVlZUybNo1169aRmprK5MmTueOOO7jqqqv45ptvmDBhAllZWUyZMoV9+/Zxww03sGDBAh544AHOPvtsnnrqKX7729/yzDPPsGzZMh5++GFfD2zp0qWVrnDdcsstPjvy8/MZOHAg/fv3JzU1Nazf/ezZsx3BqNPVEEMdo6oUFhb6xCI+Ph6wvk3rE0ePHqVbt25MmzYNgEsvvZQmTZrw5ptvVvrciBEjaNu2LTExMfTv35+srKxKy1944YUAtGvXjrPPPpuXX34ZgAULFnDFFVcQGxtLfHw8559/Ps8/X34OfurUqQBcdNFFDBo0iOeee843bxITE8PcuXPp168fK1euZMuWLVx22WUAnHLKKXTu3Jm33vJFRaB379707NnTlx/K7mCUlZXx/PPPc9VVVwHQqlUrzjnnnAq21wRP9iwMVefnn3+mtKTEJxaxsbGUlJRw9KjVi3B/42/w8A7Otm3bkpycXC6tc+fO7NxZ+YmD1q1b+143bdo0pEi2bdvW9zopKYlvv/0WsHa+zps3j2eeeQawuvb+8wTuidmcnByOHj1Khw7Hpu5OPPFEAJYtW4aIcMYZZ/jyjh49yoEDxzY1u+1OSEiotrg7dtx55500a9YMgP379zNo0KBq1RcIIxYNBGcc6xYLEal3PYtBgwZVmNDMycmp1MNTddi3b59v2TE3N5dOnToB0LVrV+655x7OP/98wPrG3r9/f9B6OnToQNOmTcnJyaFv376ANecWExND165diY+PL7eicejQoSqdxwgXx46///3vDB8+HLDmXw4fPhyxNjw5DDETnFXHXyycbb1Oz6K+MHXqVPLz830fsGXLlpGfn8/kyZY/j1atWnH48GFUlenTp1e7nVdffRWAvLw83nnnHd+wZMaMGbz44ouUlpYC8Oyzz5abSPUnJiaGyy+/3DeHUlZWxsyZM9m9ezcjRoygW7duLF5s7TQoKSlh2rRpbNy4MSwbnfe6bNkyHnnkkUrLOna4hx333Xcfzz33XFhtQegJznLbPL1yDR06VA0Wp//uTT39d2+GLDd37lz9aMqpeuTBW1VV9fvvv9cffvhBv//++wplf9j3o/6w78eI2xouq1ev1rS0NE1MTNTevXvr9OnTy+WvXLlSR4wYoaNHj9Zx48bpV1995cvbtGmT9uvXT0eOHKkPPvigLly4ULt3767Jyck6f/58feKJJzQ5OVm7d++uCxcuDNg+oA8//LBOnDhR+/btq3/+8599eUVFRXrXXXfpiBEjND09XS+99FI9dOiQqqqeccYZCuiIESN06dKlvmfy8/N15syZOmbMGB01apT+/e9/9+Vt3rxZf/GLX+i4ceM0NTVVn3rqKVVVXbp0qfbu3VsTExP13nvv1ddff933PubOnauqqo8++qj27dtXR4wYod99951OmDDB135GRoYOHDhQmzZtquedd145O0aNGqXjxo3Tm266SUtKSoL+HQL9b9i/n9Ua4HPpeR+cjZ1wPWX95je/4dxdmaSmphJ3ywOsX7+eZs2aceDAAQYOLB8+1MtzFnWBiJCVleVpnx51wfr1631DJzfBPGV5chhiqDrbt2/3DUEcmjRpQnFxsSf3WhjqH0YsGgjbtm0jIaFpubQmTZoA9W/5tLZwNlaBtXS6a9eu6BpUzzCrIQ0AVbXEYkyvculNm1riUVRUVKHX0RhxNmUZqocnexZmNaRq7N+/n4MHDwYchgD1bkXEEB1CrYZ4Uiwc5zdGLMJjx44dAEHFwgxDDOHgOMGhjp3fGOoQ5/BRvC0ODvV1r4XBmxixaAA4Jxyb2OdB3DgrIgZDTTFi0QDw9SwCiEV8fLwRC0NEMGLRAMjJySE2Npa4uIqLW3FxcUFd0kWLGTNmkJ6eXu4qKDgWwXLt2rWMGjWKMWPGcM4555CXlxe0rj/84Q8cd9xxjWJ+a/ny5UycOBE45g8jkAOh2sKTYmFWQ6pGTk4O7du3hwC+DuLj4ykpKfHcxqyMjIxyV8uWLQFrMnbq1Kn8+c9/5vPPP2fIkCFcd911vudmzJhR7v/i3nvv5cwzz6xr8+sEf0c3Y8aMYdGiRQBMmTKFu+++O6LtRWw1RESGi8hmEZlRSZl0EflBRDJc13mu/AQRWSAiK0RktYhMDFSPWQ2pGnv37qVjx44B85zehtd6F8F49913iY2NZfz48QBcffXVLF682Dcv05gRkaB+SyNBRFZD7MDItxA6MDLA/aqa7rpedduDFbl9JHAx8LKIJAesxRA2OTk55fwpuHHmMbwmFtdddx1jx47lrLPOYtmyYxEsV61aRZ8+fXz3Xbt2pXnz5nz99dc88sgjvPfeeyxYsID09HSeeuopX7nc3Fwuu+wy+vfvzxVXXBGwzR07djBy5EhEhEceeYTTTjuN/v378+yzz/rKBHNN5352wYIFTJw4kaZNm7Jt2zb27t3L+eefz7hx4xg1ahS33HILR44cASyfouPGjSMtLY3TTjuNH374ATjmlfzCCy8M6A7wrLPOAqydpunp6Wzfvt3XfjD27t3Lueeey7hx4xg9enTEfZeGu4Nzlaq+Znv4rhZ27NOrgXMBVHWjiKzBilw2t7r1GiyxGDp0aMA8Ryxk8b8oybG2N6fYbvVK4poGfKa6SJfjiT2/QoC5CvTt25cJEyYwfPhwVq1axfjx41m+fDmDBg1iz549Fb4927Rpw969e7nppptYs2YNPXr0qNDr/Pzzz1mxYgWqSufOnfnyyy8ZNaq8i79u3brx8ssv07NnT5o0acLSpUvJysri5JNPZsSIEfTp06dS13TOszExMXzwwQfMmzePpk2bcumll/qGCEVFRYwcOZI9e/bQtm1bzjzzTF555RUmTJjA22+/zdSpU1m/fj3XXHMN2dnZ/POf//S5AzzllFN87gDfffddRISXX37Zd+DNaT8Yl1xyCaNGjeIPf/gDP/30E/37949ouICwehZatcDI00RkmR1P5B4RcQTpeCAJ+MFV9nugwuk2Q9WorGfhDEO8NGdx1113+Ry0DB8+nEmTJvHEE0/48gN9e4Y6HT1+/HiaNm1KQkICJ510Uthu9Xr27MmoUaP497//HbZrOset3q233kpZWRkffvghV155JWAtVT/99NN06NCBt956i5YtWzJhwgQAJk2axO7du1m5cqWvrqq6AwzGrl27+Oijj5g5cyYAnTp1IjU1lZdeeqla9QUi0mdDDgBfYvUUEoC3gLbAbVjRx5wyDvuBfv6VOBOcDrNmzTLzF0EoKipi//79tlgUVMh3ehb5Ey+kxXHHAZDtsSPq3bp1Y926dQB07NiRNWvWlMvfv39/0DkZh5q61fvpp5/Cdk3n7vk47v7cYu2U37lzJ/v27Su3YtGhQ4dyqzuRcqvn2HH55Zf7Pju5ubkMGBAwKkc5Zs+ezZw5c9xJte/dW1XXAM5fukBEHgBeFZHb3cX8HqvwNWK8e4dPbm4uYH3IOFxRLGJiYoiJifHUXosHH3yQO+88Fslyz549Prd5w4cP9znPBfjxxx85fPhw0GFWddm3bx/t2lmB9HJzc+nfv3+1XNN17doVsHp3TkTyrVu30rZtW7p27UqXLl3KHV4LdIYnEjh2vPrqqz7hKiwsDGuuyvHq7RAt7947sKKVdQD22mlu76dtXOmGauCsEgQbhoiI5/ZazJs3z7eRLCsri9dff93nAfuss86ipKSETz6xAtU9/fTTnHvuub7357iaO3ToEJdcckm1bXDc6m3dupUVK1ZwwQUXVMs1XUpKCmeccYbPrd7Ro0e54IILKCoqYvLkyeTm5rJq1SrA8r85fvz4cg57K6Nly5YcPnyYF154wWdvKDvctl933XXlJo9rSkR7FiLyG+BJVXUivCQDRUAekAvsA3pzTCD6Ae9E0obGhvOh69ChA2wPXMZruzhvv/12pk+fTlxcHIcOHeLvf/87aWlpgDWEWLJkCddffz2xsbG0bdu2XIyQiy++mCuvvJKMjAxuvvlm5s2bx3vvvUdCQgIDBw7k22+/Ze3atdx///1069bNN1/gj4hw+umn89NPP/GPf/zDtwIzb948br75ZkaPHk18fDyDBw/m3nvvZd++fb55jvT0dObPn++LO/rCCy9www03MG7cOEpLS7n77rt9HsrfeecdbrvtNp9rujlz5tChQwdefPFFFixYQGFhIf/4xz+IjY31vY+TTjqJiy++mBtuuIELL7yQ1q1b89JLL/mcCKenp3PLLbdw//33s3v3bn7961/7Yqxcf/31vPbaa4AlvOecc07k/nCBfO0Fu4AMYIbrvj3wGdDOvl8AXG6/jgVeA55xlb/fuQdOxBKR4/zbMT44jxHKB+fChQsV0B9++EGL592pxfPuVNXy/hU3btyo3333ne8+2j44o0lWVpZa//aGqvrgjHRg5H8CF4rIMmAFkAPc5KpqtlWdrABeAi5S1d3hyZohEKGGIXBsF6fBUBPCXTrNVGuDVRtV7aOqzl6JH1W1g9pLq6r6uaqerarjVXW4ql6jqgdd9RSq6gxVHamqw1T1g0Dtme3e4bN3715iY2MrDZjrDEPUg86Z65IdO3aUG0o47v4NFibWaQPH2WMRExNDsJ0U7i3fgU6mNha6devGihUrom2GZzGxThs4lW3IcvDqlm9D/cKIRT2nKmLhpRURQ/3DiEU9Jzc31zqeXgn17eSpwZsYsajn5OXlkZSUVGkZIxaGSOBJsTCrIeFRVlbGvn37QopFbGwsYMTCUDkmFEAD5sCBA5SVlYUUi5iYGGJjYz0jFuvXr2f8+PGkpqYyePBgHnjggYDlZs2aFfJ4tXGrlx6x+k0ogAaMc3oxlFiAt3xxXnHFFYwePZrly5fz0Ucf8cADD/Dee++VK7Nnzx7+9a9/VXjWuNWrPbd6oTBiUY+pr2Kxbt06n2OapKQkTjzxxArH0v/4xz/yv//7v9Ewz7PUtlu9UBixqMfUV7GYNGkSb775JmCd+ly3bh0jRozw5W/evJns7GxSU1PLPWfc6tUPt3oGD1IVsXhg4xLW7d9Oi40tOWy71WseYbd6A9sfz9zUa0KWe+qpp5gyZQq9evVi3759zJs3r9zp0HvvvZc5c+ZUCAFg3OrVA7d6dY1ZDQmPqohFTEwMZWXeOBty7rnnMmrUKLZs2cI333zDvHnz+PrrrwFYvXo18fHxYXl4cmPc6tXcrZ45G9KAyc3NJSYmJqxx7OyBF5Gdnc2QIUPYdMD63UbDrd769ev56KOPWLhwIWB5eDrjjDN46KGHeOmll5g1axaPPfZYles1bvWq71bPIdTZEE+KhSE88vLyaNeuHTExoTuIXtmY5XwY3Afa4uPjOXjwIPn5+WzevJkZM2YA1gd19+7dpKen86tf/SqiE57GrV7V8eQwxBAeeXl5Ibd6OzhiEe1j2X369KFz586+7nF+fj5vvPEGEyZMoFWrVmzYsMEXpezhhx/muOOOIyMjwycUxq1e5XbUplu9KnnKqqvLeMo6RmWesiZMmKBjxozx3QfzlKWqeuDAAV21apUePHgw6p6yVq1apWPHjtXU1FQdNGiQ3n777VpcXFyuzJw5c3TgwIHatGlTTUtL05UrV6qq6hdffKG9e/fW4cOH68KFC3Xu3LmanJys3bt318WLF+ucOXM0MTFRe/furUuXLq3QtuMp68knn9TTTjtN+/XrpwsWLPDl5+fn68yZM3XUqFE6btw4vemmm7SkpETz8vJ0xIgRCmhaWpquW7fO98yePXv0vPPO07Fjx+ro0aN10aJFvrzVq1drWlqajhs3TseOHatvvmn9LRcuXKjdu3fX5ORknT9/vj7xxBO+97Fw4UJVVb3rrrt0wIABOmbMGN2xY0e59pcsWaK9e/fWxMREvfHGG312/PKXv9TU1FRNTU3VP/3pT5X+HarqKSvqwhDo6tSpkwI6a9asSt9sY6AysRg4cKBOmTLFd1+ZWBw6dEhXrVql+/bti7pYRBPjVu8Y/v8js2bNUizv+9laXbd6EHas01QReUtElopIpoj8yY5E5uQv8IuDmiEiLf3rMdu9wyOcQ2QOXpmzMHiXUNu9w5rgtGOdnk/oWKf3AfNU9Q1bBFZjefJ+xCmgqunhtGkITbXFopFOa/u71Vu6dKnvkJ0hNJGOdboEeBNAVQtE5C1gIi6xMESGI0eOcOTIkbDFwgk2ZPUsgu8CbMgYt3o1I6KxTlX1YdVyXmETsDx8+xCRx+04qO+KyPjwTTW4qcqGLAcvbfk21D9qrUMqIrHAL7CipDusBz5W1VUiMhxYJiKpqrq2tuxoqIQjFqpa7izBMbFovE57DRbVCZRdm/ssbgeWqKpvb6uqPqCqq+zXq4C3gWv9H3S2ezuXmeisSCixSEhIIC8vr5z7f9OzMKgqRUVF7Nq1ixYtWgDHtnk7F3W53VtEzgaGAReGKLoDONk/0Wz3Dk0osejSpQs7d+70BSECa5dhUVERcYetmFBlzfNr31CD54iLiyMxMdG3oS/cwL4fkJoAAB2BSURBVMgRFwsRGYkVhWyqqpaKyImqusnOu1NVH3QVTybIMo2hcpzo6cF2cMbHx1c4ofjYY4/x4osvMuRp62ToR9Pur10jDQ2KGg1DRKS9PVnZzr7vi7XycTUQZy+f3ut65FYR6WiX7QlMBZ7HUGWcHkNVJjjbt2/P/v37G31kMkP1CHefxVBgLsdinU5RK4ShO9bpPuBfwKlYwwsHd2zvvwKviUgJ0AK4UVU/qfG7aITk5ubSpk2bKkUYS0pKQlUpKSkmPr5JLVpnaIiEJRaqmgmkB0j/Eejguh8Top6/YgmGoYaEE1zIH2fIUlxcYsTCUGU8eerUOL8JTU5OTtgnTh2cIYuJTGYIhAkF0EDJzc2tcs/CiIWhMkwogAZKzYYhRiwMVceIRT1EVcOKceqPr2dRYsTCUHWMWNRD8vPzKSoqqnLPokWLFjRp0oQS07MwVANPioWZ4KwcZ49FVXsWIkL79u0pLjZbvg0VMROcDRBn92ZVexZgDUXMnIUhEGaCswHi9CyMWBjqEiMW9ZDqDkOcZ4xYGKqDEYt6SI2HIWY1xFANjFjUQ3JyckhISPD5I6gK7du3p6S4xBwmM1QZT4qFWQ2pHGerd2URtYPhHCYrLTUrIobymNWQBkh1tno7HNvybcTCUB6zGtIAqc4hMgez5dtQXYxY1EMi07MwYmGoGkYs6iE16VkYsTBUFyMW9YzDhw9z8OBBOnXqVK3nfcMQs3xqqCKeFAuzGhKcn376CbAmgatDYmIiImJ6FoYKRGw1JJzAyHa5S+ygyKtFZK641vdEpJ2IvCEiy0VkhYgMCVSHWQ0Jzq5du4Dqi4WIEB8fb06eGioQkdUQOzDyLYQIjCwi/bEc+/4Cy3HvEOB6V5H5wBpVTQV+C7wuIk3DscFg4cRTqa5YAMTFx5mlU0OVCbdnsUpVLwZCRaWZCbyjqrmqWgY8DVwHVq8CKxL7UwCqugwoAiZXx/DGSiTEIj4+3gxDDFUmooGRgeHAD67774GTRaQZVi/jqKru8MsfFmbdBiyxaNasGYmJidWuw4iFoTpEeoIzmfJDlf2AAO0D5Dn5Hf0rMbFOg5OdnU1KSkq1tno7xMcZsTAcI9xYp7WxGhLohJKEkefDmeB0LiMWx3DEoiY4PQtzmMwAxyY2nYs62u69F2jjum+DJRA5dp5/37mNnW4Ik0iJhaqSn28CIxvCJ9JisQro7brvB6xT1SPA10CCiHT1y18VYRsaLKrKrl27IiIWcMwvhsEQDhENjIwV6/RsEUkSkRhgBvA4gKrmAYuAq+xn04AmwNs1saExcfDgQQ4fPhwxscjLy4uEWYZGQrj7LIaKSAbHAiMvtrPcgZFR1e+A24EPgJXAWqy9FQ7XA0NFZDlwPzBVVQsj8D4aBZFYNgXTszBUj3CXTjNVNV1V26hqHzuCOqr6o6p2cC+tqupCVR2qqsNV9TZ1zaKp6j5VnaKqqao6SlW/DtSe2e4dmMiJhRUP24iFwY1xftOAiHTPwgxDDG6M85sGhCMW1T1x6hAbG4eImJ6FoUoYsahH7Ny5k1atWtGqVasa1SMixMXHmZ6FoUoYsahHbN68mV69ekWkrvj4eNOzMFQJIxb1iM2bN3PiiSdGpK74OCMWhqrhSbEwqyEVKS4uJisrK3JiER9vhiGGcpjVkAbCtm3bKC0tjahYmJ6FwY1ZDWkgbN68GSDiPQtzmMwQLkYs6gmbNm0C4IQTTohIffHx8RQVFVFQUBCR+gwNHyMW9YRNmzbRqlUrOnas4P6jWphdnIaqYsSinrBp0yZOPPHEGjm9cWN2cRqqiifFwqyGVCSSy6ZgDpMZKmJWQxoAqsq2bdsiKhZxcUYsDOUxqyENgMIjhRFdNgUzDDFUHSMW9YD8Asv93cCBAyNWZ1xcHDExMaZnYQgbIxb1gPz8AhISEujXr1/E6hQR2rVrZ3oWhrDxpFiYCc7y5OcfZNCgQb6hQ6Ro37696VkYfISa4IwLpxIRScDypdnHfub/qeoHAcq9ByS4kuKBAaraOkh+gapWiEiWkpLi893Q6FEoKChg2LDIx2IyYmFwM3v2bEcwAn74whILYDYgqjpSRE4CVohIX1Xd41duq6r6YpvaMVKnuPJ3q+qM8M03HD5ymNLS0loRi6SkJLKysiJer6FhEnIYYnvpvppjMUo3AmuAS/3LuoXC5hLghZqb2XhxYnuYnoUh2oQzZ3E8kETFGKaV/veKSCKWN/BlruQWIrLQDh+wREQiN73fQMnPzycmJoY+ffpEvO6kpCRzmMwQNuGIRbL90z+GaahDCr8E/mNHU3fYAsxS1bFYPZXlIlIzh5INnPz8fFq2bEVsbGzE627fvj1Hjx7l0KFDEa/b0PCoymqI/9dPqEMKlwDPl6tA9W5V3Wy/fhOrh3K5/4MmMLLF0aNHyc/PJ7F161qpv3379oDZxdnYiWRgZCcWqX8M06AxSkWkM9DGDjpUGTuA7v6JJjCyRWZmJqpK68TaEYukpCTA7OJs7EQyMPIWYB8VY5hWFqP0IuBFd4KIdBSRGX7lkoMZZoAvvvgCgMTW/vGkI4PpWRiqQkixsOcc/smxGKUnYk1cLhSRviKyVET8B9QX4icWWCEObxORFnY9w4HhwCs1ewsNly+++IKEhGbEN4nsZiwHRyxMz8IQDlXZZ/G4iKywn7lIVXeLSA+sjVrxQCmAiPQD8lT1J786dgNLgA9FpAQrKPJ0Vd1U0zfREFFVvvjiCxJTR9daG84wxPQsDOEQlljYwYtnBEhfAXT2S/se+EWQOn5vX4YQZGVlsWfPHvrV0hAEoE2bNuYwmSFszNkQj/Lll18C1NpKCEBsbCxt27Y1wxADYJzf1Fu+/vprEhISaN6iea22Y3ZxGhyM85t6ypo1axgwYEDEfG4Go0OHDuTk5NRqG4aGgRELD6KqrFmzhiFDhtR6Wx07dmTv3qBbZgwGH0YsPMi2bdvYv38/gwcPrvW2jFgYwsWTYtHYJzjXrFkDUCdikZycTG5uLiUlJbXelsHbmAnOesiaNWuIjY1lwIABtd5Wx44dUVWzImIwE5z1kTVr1tC3b1+aNWtW6205Ec7MUMQQCiMWHmTNmjV1MgQBIxaG8DFi4TEOHDhAdnY2/fv3r5P2kpMtdyV79vh7SDQYymPEwmNs3LgRgN69e4coGRlMz8IQLp4Ui8a8GrJhwwYATjrppDppr02bNsTFxRmxMJjVkPrGxo0biYmJoVevXnXSnojQsWNHMwwxmNWQ+saGDRvo2bMnTZo0qbM2k5OTTc/CEBIjFh5j48aNdTZf4WB2cRrCwYiFhygrK2Pjxo11Nl/hYMTCEA5GLDxEdnY2hw8fjkrPYs+ePSZ+iKFSPCkWjXU1pK5XQhySk5M5cuSIiR/SyInIaoiIJIjIAhFZISKrRWRikHIzRGStiGS4rlRXfjsReUNEltt1BTyD3VhXQ+p6j4WD2WthgNCrIZEOjAxws6pmBKlnPrBGVWeJyHjgdRE5QVWPhmlHg2bjxo00b96clJSAwl5ruMXi+OOPr9O2DfWHiAZGDlFPO+B8Vz3LgCJgctVMbrhs3bqV448/vta9Y/njbPnevXt3nbZrqF/URmDka0XkE3sI4o6qPgQ4qqo7wqyn0ZGVlUXPnj3rvN1Onaxwsz/95B+9wWA4RqQDI+8BPgDSgQuAG0Xk1656DviVD1hPY4x1qqpRE4uOHTsSExNjxKKREslYpw4hAyOr6ruq+oxa7AUeBdy9i0BrcxXqaYyxTvPy8igoKIiKWMTGxpKcnEx2tokk2RiJZKzTKgdGduEOfLwX8I+YE249DZ6srCyAqE0wdurUyfQsDJUS0cDIInKnX5I78PHXQIKIdA1VT2PEEYto9CzA6s2ZnoWhMiIdGPksEUmzyzUHrgWet+vJAxa56knDinf6dkTfUT1l69atQPTEolOnTkYsDJUS6cDIc4F77eXWlsBS4C+ueq4HFojIciAWmGrHQG30ZGVl0b59e1q2bBmV9lNSUsjJyaG4uJj4+NqJ2m6o34Q1wamqhao6Q1VHquowVf3ATl+hqp2dD7yqvqWqp6nqeFUdrqp3q2qRq559qjpFVVNVdZSqfh2ovca43TtaKyEOnTp1QlWNX4tGjHF+U0+Itlg4u0bNJGfjxTi/qQeUlpayffv2qPcsADNvYQiKEQsPkJ2dTXFxcVTPZZiehSEURiw8QLRXQsDaxSkipmdhCIoRCw8Q7T0WAHFxcSQnJ5uehSEonhSLxrYakpWVhYjQrVu3qNph9lo0bsxqSD0gKyuLLl261KlH70CkpKSYnkUjxqyG1AOivWzqkJKSws6dO6NthsGjGLHwAFlZWZ7wUNWtWzdycnI4cuRItE0xeBAjFlGmsLCQ7OxsT/Qsune3Dgjv2LEjRElDY8SIRZTZvn07quopsdi+fXuULTF4EU+KRWNaDfHCsqmDEYvGjVkN8TheEovOnTsTGxtrxKKRYlZDPE5WVhZNmzb1nc2IJnFxcXTu3NmIhSEgRiyiTFZWFt27dycmxht/iu7du5sJTkNAvPEf2ojZsmWLJ5ZNHbp37256FoaAeFIsGssEp6qycePGOg9XWBndu3dn586dlJSURNsUQx1T17FOzxGRD0TkYxH5WkR+45f/nl8c1LcC1dNYJjizs7M5dOiQ58SitLTUnBFphNR1rNO/Aeer6hoROQ5YJyLbVfV1O3+3qs6o8rtooEQranpluJdPo32wzeAtIh3r9DFVXWOX2w0sAwL2QgzRi5peGWavhSEYEY11qqp/80tKAHJc9y1EZKGIfCYiS0RkYFUNbkhEK2p6ZTi9iW3btkXXEIPniHSsUx8i0hoYDjzjSt4CzFLVsVg9leUiEv0NBlFiw4YNnHjiiZ5ZNgVo1qwZXbt29fV6DAaHiMY69eMB4A+q6uvP2qEBNtuv38TqoVzu/2BjCYzstZUQh969e/PDDz+ELmhoEEQyMHKVY52KyDVAsao+FqJudyxUH40hMHJRURFZWVmemtx06NOnDxs2bHBmxg0NnEgGRg471imAiEwHxgM32fcn2j87isgMv+LuWKiNiq1bt1JaWurJnkWfPn04ePAgu3fvjrYpBg8R0VinIjIO+A1wI9ZkZkvgHruq5sBtItLCLjsca07jlci+pfqBF5dNHRwBM0MRg5tw5yxmA2LHOn0JO9YpkMixWKfYeelALpBvX84wYzewBPhQRD4FHgWmq+qmmr+N+sfatWsREfr16xdtUyrQp08fwIiFoTxhbcqyY5nOCJC+Aujsuu/sX8avjt/bV6MnMzOTPn36RC0QcmV07tyZFi1a+Ho/BgOYsyFRIzMzkyFDhkTbjICIiFkRaYQY5zceZPfu3WRnZzN06NBomxKUPn36GLFoZBjnNx7k66+/BvC0WPTu3Zvt27dz+PDhaJti8AhGLKJAZmYmAIMHD46yJcE55ZRTAPjvf/8bZUsMXsGIRRTIzMzkpJNOolWrVtE2JSgjRowAYOXKlVG2xOAVPCkWDXmCU1VZvXq1p4cgYMU97dKlixGLRoSZ4PQYGzZsYNeuXaSmpkbblJCMGDGCr776KtpmGOoIM8HpMd555x0AJk2aFGVLQnPqqaeydetWcnJyQhc2NHiMWNQx7777LieffLLPyYyXceYtVq0KeAzI0MgwYlGHFBQU8Mknn3D22WdH25SwGDp0KDExMWbewgAYsahTli5dSnFxcb0Ri5YtWzJgwAA+/fTTaJti8ACeFIuGuhry8ssvk5iYyJgxY6JtSthMmjSJzz77jLy8vGibYqhlzGqIR9i5cyevvvoqM2fOJD4+PvQDHmH69OmUlpby1lsBozYYGhBmNcQjzJ8/n7KyMm688cZom1Ilhg4dSpcuXXjttdeibYohyhixqAMKCgp44oknmDZtmieipVcFEWHatGm8//77HDp0KNrmGKKIEYs64I477uDnn3/mzjvvjLYp1eLcc8+lsLCQRYsWRdsUQxQxYlHLvP/++zz++OPceuutvn0L9Y309HQGDx7MfffdZ2KgNmI8KRYNZTVk+fLlXHDBBfTr14/77rsv2uZUGxFhzpw5bNmyheeffz7a5hhqiToNjGyXvU1EMu3rDr+8HiKyzI5IliEiAQfw9X01pKysjCeeeIKJEyfSqVMn3nvvPRISEqJtVo2YPHkyw4YN45577uGnn36KtjmGWiBSqyGzsQMjAxcDL4tIsn8hETkT+B9gjH3NFBH3IYiXgOfsiGTP0sA8excWFrJw4UKGDx/Oddddx8iRI/n000/p2rVrtE2rMSLCk08+yYEDB5g8eTIFBQXRNslQx0Q6MPK1wIuqWmg76F0IXGfXMxA7hIBddiEwQES8fVY7AKrKgQMH+P777/n3v//N73//eyZNmkRSUhKXXnopBQUFPP/88yxdupSOHSuN8livGDx4MK+88gpr165l8ODBvP/++5SVlUXbLEMdEY5377ADI2PFAXnJr9z1rrwsVS0CUNUiEdls15NZRbvJzs4mIyODsrIy31VaWlruvrIrWNmSkhIKCws5cuQIhYWFvtdHjhwhLy+PvXv3kpOTQ1FRkc+WmJgYevfuzZVXXsm0adOYMGGCp+KXRpJJkybx4Ycfcs0113DmmWfSuXNnxo0bR69evUhMTKRFixa0aNGCuLhj/1r2OLharw2R4+yzz6ZFixbVfj4csQgWGDlQwIvkAOU6Bsnzz/fhTHA6zJo1q8L8xTfffMMll1wS2vowERFiYmKIjY2lWbNmJCQk+C7nPiUlhUGDBtGhQwc6duzIcccdR79+/ejbt2+9n5OoChMmTODbb7/l1VdfZcmSJXz55Ze8/PLLJtyhx9m6dWvAfT6zZ89mzpw57qSAE5xhxQ2xCTcwcmX/MYHyKtSTkpJCdnblUQ3Hjh3Lhg0biImJKXfFxsZWSAt2OWVdAWENYdKsWTMuu+wyLrvsMgBKS0s5dOiQ73KGJ24BqeprQ2Tp3DlwWJ/Zs2eX+zIWkYAfvnDEwh0YOdBr/7L+AZRzguRVVk9IWrRo4cnQf42V2NhYWrduTevWraNtiqGWiHRg5FWVlFsF9BSRJgD2z15B6jEYDB4jooGRgceBi+x9GQlYy6yP2/WsBf4LXGSXvQhYp6pVntw0GAx1T7hzFrOBx+3AyHHYgZFFpAfHAiOXqup7InIy8Ln93NOq+rarnguBp0XkaqAUuKDmb8FgMNQFYa3v2fsmZqjqSFUdpqof2OkrVLWzvafCKTtXVYfa10N+9WxT1QmqOlZV01U1K1B74W739voOTy/bt23btmibEBQv/94asm2htnuLF2efhw0bpqtXrw5ZTkQ8PXseCfvOuMdyOvPhfZNDli35210AxN3yQKXlTl9yNxkZGZQ8vKJGttUWXv67NgbbRCRTVSvso/L8zqGaqGVNeya13bOp7Hkv2xYq38u2hcr3sm3h5Ndm3Z7vWVSmlqGUtD7nO3nBehaBnnX3LCqrO5yeRSRsr426azvfy7aFyo9U3cF6Fp4UCxHJAbbbtykEOQUXIq++53vZtlD5XrYtVL6XbQuVH6m6u6tqB/9MT4qFwWDwHp6fszAYDN7AiIXBYAiLeikWtteuDL+rZbTtciMiV4vI5yKyXES+EZE0D9i02+93tk5E3oi2XQ72juBl9u9sjYjcFW2bHEQkWURetb28fSsiV0bZnuEisllEZvilDxKRL+3/vTdFJClSbVbl1KmnUNX0aNsQDBE5HzgNGKeqpfY/1nFRNgtgsao6/kUQkUeAr6Nojz/PAh+q6u/sf/JNIvJfVX0v2oZh+Wn5RFXPE5FuwHcisl5V63yziohMB87Hz+WDfd7qdWCGqi4TkTlYxy3Oj0jDqlrvLmBBtG0IYd+3QN9o2xHCxlhgK9Aq2ra4bDoETHbdrwR+6wG7umC5V+jqSvsP8M9o2WP/zLCFwUmfCmx13XfFOlbRIRLt1sthCICIPG53Cd8VkfHRtsdBRDoCfYFBLufE10bbrgCcDqxQ1fxoG+LibeAcABE5HjgZSzCijeMIIteVtgfrQGWdo6o7g2QNx+XRTlV/BA4DQyLRbn0dhqwHPlbVVSIyHFgmIqlqnWyNNj2wHPpMx/pAdgS+EpEDqvpyNA3z4xLAa379ZwJviMgWoB1wq6p+HGWbAJwPZ2dgs/36OKBtdMwJStje6KpDvexZqOoDqrrKfr0K6xvJK9/eTbF+r4+qaqmq/oT1obwqumYdQ0SaA6OBD6Jtix+LgS9VtRdwCnCriETkW7EmqOou4CPgZrHoizUnVVj5k1EhLG901aFeikUAdgDdo22Ezc/2zz2utJ1Y416vMBV4S1VLo22Ig/0BPB14GHxd6A+BOyp7rg65EGiO5X5hFvAvju0y9goR9UbnT70UCxHxDxqaTOXbWOuSTVjjRHfXrwPesQ+8OQRpYv8sdqUVA57w06eqeap6laqOVtULsT6EXhgiuSnnqU5EumIJXEQcTNVLscDqnnYEECuq2VQ88s+vqkeB57DG34hIC+BXWMuCUUdE2mPN6nvNQ9kPwC5sT2oi0gqYgkc+kCLyvP27Q0T6AalYvQsv8S4Q59rTcxXWcnlOJc+ETb08GyIit2NNIJYALYBHVNUTYgE+gXgcywdpCdYy20PqgV+2iFwPtFHVP0fbFn9EZBgwD2uM3RJrnuC3qhr1aMwi8hcs8foZKMCafP0+SrYMBeZircbsBr5X1XPtvMHAfKwl05+xllbzItKuB/5/DQZDPaC+DkMMBkMdY8TCYDCEhRELg8EQFkYsDAZDWBixMBgMYWHEwmAwhIURiwaGiJwhImtFREXkExFpVwttpIjIHBFJF5FtLmc6KiL/tV+vFZHZInKliIwKo84YEblLRHpV06Zb7S3jhlrC7LNogIhIOrAMiI/0hiaxYth+ghXHtiuQrqqz7TwFzlDVj2wb0oH7sTZXXaWqGyup9x6sEJh/qaZdScD7wJmqmhuqvKHq1Ncj6obo8Wvga1XdIiL7CX7m5VsgW1ULReRR4G/ApEAFRSQF+A3QrbpGqWqeiCwC5gA3VLceQ3DMMKQRYvtv/EREPrV/Dnfl9bId9nwpIgtFZLE91LjGLnIB8Cn4DlcF7C345X0CnFHJkGg6sMYWlma2j8sjIvIv26Y7RCRXROaKyHwR2S8if7Rt2yEis13tnCcisTX6BRkCEw23YOaq3Qur+69AXIC8RCyPT+n2/Vj7vo1973Nlh+Xs5QAw274XrJOgI4K0q8DpQfIKgLQgeU8B/3DddwSOcsx9XBOsA1FOfgaWL45YrFOWV9vp7W0bukf7b9AQL9OzaHxMBg6qagaAqn6GdeBoioh0B04FXrDzdmF9WzskYQ1dq+OKLx/LlUAgkt11qupeLF8Wl9pJk7BOVLp5Wy3nQhtU1Tn96dQRrB1DDTBi0fjoAvgfWc6x0zvZ9+4Jwn2u147HperMiivB/98kQJ3PAZfZry8AFvnl+7uPc9tl/q9rAfNLbSSISKKITAJ+xHLG46YDljevn1z3Du64E3lYR59bVcOEVpT3HuZmb4A63wBSROQXQJmq7g+zDSppx1ADjFg0HtpixY94C2glIuMARGSMnfeGqm4HvsL+RheRzsAYpwJVLQPWUUUXhvayZjP72UB841+nqhZi9SaeAV4Js6keWEOqYN6vDTXALJ02MERkLNbyIcAr9t4HsNyr7VXVgyJyJjBXRGKwuu5nub65LwaeFZHJWN6rPqL8EOFVLAHxDQvsVY7F9u1cEblFy3vlHgt8as9FBGIxcLeIxKuq263ec8C5uOYrRORBLKcvd4tIO1Wd59fOEr86DBHCbMoylMP+AO5z3b+N5dz3H/Z9K6xJz0lqeS4PVV8c1mTlLVpJqAZbBLap6nxXWl/gRlUNuW9CrPCVGcBUe2LWEGHMMMTgz8Mi0gd8Dl9HA0udTLWCEv0KuCXM+mYAD1YmFDa/A7qLyPEicpEtMjMI33fpDcD1RihqD9OzMJRDRC7GEoJDWH4wH1LVcOcMImXDQ8CZwFeqOrMu2zYEx4iFwWAICzMMMRgMYWHEwmAwhIURC4PBEBZGLAwGQ1gYsTAYDGFhxMJgMITF/wdtYQGuMQOclAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 288x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "id": "Dd2qIbjNMzHc",
        "outputId": "86968e2d-977a-4936-bbcd-602d404db574"
      },
      "source": [
        "PDF_columns"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            },
            "text/plain": [
              "'816'"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 52
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "25dx6ykZM8is",
        "outputId": "4b972f49-bd85-4f9e-ecdd-362cf459e84b"
      },
      "source": [
        "len(cluster_list[0])"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "4"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 49
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "zI88JQGkaaUT",
        "outputId": "797eb8bf-b430-4c29-bf84-7f17478fd903"
      },
      "source": [
        "a = np.where(PDF_columns.str.contains(\"tpdf_cb2018\") == True)\n",
        "age_CB2018 = PDF_columns[a]"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<bound method Index.sort of Index(['tpdf_cb2018_z014_u3_ID383', 'tpdf_cb2018_z014_u3_ID451',\n",
              "       'tpdf_cb2018_z014_u3_ID1252', 'tpdf_cb2018_z014_u3_ID1381',\n",
              "       'tpdf_cb2018_z014_u3_ID1564', 'tpdf_cb2018_z014_u3_ID1576',\n",
              "       'tpdf_cb2018_z014_u3_ID1949', 'tpdf_cb2018_z014_u3_ID2449',\n",
              "       'tpdf_cb2018_z014_u3_ID2732', 'tpdf_cb2018_z014_u3_ID2740',\n",
              "       'tpdf_cb2018_z014_u3_ID534', 'tpdf_cb2018_z014_u3_ID417',\n",
              "       'tpdf_cb2018_z014_u3_ID589', 'tpdf_cb2018_z014_u3_ID11',\n",
              "       'tpdf_cb2018_z014_u3_ID894', 'tpdf_cb2018_z014_u3_ID90',\n",
              "       'tpdf_cb2018_z014_u3_ID93', 'tpdf_cb2018_z014_u3_ID531',\n",
              "       'tpdf_cb2018_z014_u3_ID816'],\n",
              "      dtype='object')>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 14
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "GuiJ0kDaa6xI",
        "outputId": "05fc2b4d-a74f-49bc-f35a-91e7b3ccd7be"
      },
      "source": [
        "PDF_columns"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Index(['AVt_cb2018_z014_u3_ID383', 'AVpdf_cb2018_z014_u3_ID383',\n",
              "       'logt_cb2018_z014_u3_ID383', 'tpdf_cb2018_z014_u3_ID383',\n",
              "       'AVt_modp020_chabrier_MW_phi0.50_ID383',\n",
              "       'AVpdf_modp020_chabrier_MW_phi0.50_ID383',\n",
              "       'logt_modp020_chabrier_MW_phi0.50_ID383',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID383', 'AVt_cb2018_z014_u3_ID451',\n",
              "       'AVpdf_cb2018_z014_u3_ID451',\n",
              "       ...\n",
              "       'logt_modp020_chabrier_MW_phi0.50_ID531',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID531', 'AVt_cb2018_z014_u3_ID816',\n",
              "       'AVpdf_cb2018_z014_u3_ID816', 'logt_cb2018_z014_u3_ID816',\n",
              "       'tpdf_cb2018_z014_u3_ID816', 'AVt_modp020_chabrier_MW_phi0.50_ID816',\n",
              "       'AVpdf_modp020_chabrier_MW_phi0.50_ID816',\n",
              "       'logt_modp020_chabrier_MW_phi0.50_ID816',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID816'],\n",
              "      dtype='object', length=152)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 15
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "uZN1R6Albhyw",
        "outputId": "6eed90be-350a-4bc5-d8e1-c439374b7efd"
      },
      "source": [
        "axes[1]"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.axes._subplots.AxesSubplot at 0x7f96dcf208d0>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 21
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "CVQSIYxvcFGJ",
        "outputId": "2b99bfc9-a03f-4363-e4e5-d7cb0a4fb7c9"
      },
      "source": [
        "for i in range(2):\n",
        "  print(i)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "0\n",
            "1\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "tj0sch4ZcI3u",
        "outputId": "3c62c835-de56-494d-b1a0-5cabdfe75310"
      },
      "source": [
        "for CB2018, K2015 in zip(age_CB2018, age_Krumholz2018):\n",
        "  print(CB2018)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "t_16th_cb2018_z014_u3_383\n",
            "t_50th_cb2018_z014_u3_383\n",
            "t_84th_cb2018_z014_u3_383\n",
            "t_16th_cb2018_z014_u3_451\n",
            "t_50th_cb2018_z014_u3_451\n",
            "t_84th_cb2018_z014_u3_451\n",
            "t_16th_cb2018_z014_u3_1252\n",
            "t_50th_cb2018_z014_u3_1252\n",
            "t_84th_cb2018_z014_u3_1252\n",
            "t_16th_cb2018_z014_u3_1381\n",
            "t_50th_cb2018_z014_u3_1381\n",
            "t_84th_cb2018_z014_u3_1381\n",
            "t_16th_cb2018_z014_u3_1564\n",
            "t_50th_cb2018_z014_u3_1564\n",
            "t_84th_cb2018_z014_u3_1564\n",
            "t_16th_cb2018_z014_u3_1576\n",
            "t_50th_cb2018_z014_u3_1576\n",
            "t_84th_cb2018_z014_u3_1576\n",
            "t_16th_cb2018_z014_u3_1949\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7HHMbshmca2R"
      },
      "source": [
        "a = np.where(percentiles_columns.str.contains(\"cb2018\") == True)\n",
        "age_CB2018 = percentiles_columns[a]\n",
        "\n",
        "a = np.where(age_CB2018.str.contains(\"t_\") == True)\n",
        "age_CB2018 = age_CB2018[a]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "kUq6xHB6cq1d",
        "outputId": "3ad2c7df-7e63-4822-f67b-7af391254d23"
      },
      "source": [
        "age_CB2018"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Index(['t_16th_cb2018_z014_u3_383', 't_50th_cb2018_z014_u3_383',\n",
              "       't_84th_cb2018_z014_u3_383', 't_16th_cb2018_z014_u3_451',\n",
              "       't_50th_cb2018_z014_u3_451', 't_84th_cb2018_z014_u3_451',\n",
              "       't_16th_cb2018_z014_u3_1252', 't_50th_cb2018_z014_u3_1252',\n",
              "       't_84th_cb2018_z014_u3_1252', 't_16th_cb2018_z014_u3_1381',\n",
              "       't_50th_cb2018_z014_u3_1381', 't_84th_cb2018_z014_u3_1381',\n",
              "       't_16th_cb2018_z014_u3_1564', 't_50th_cb2018_z014_u3_1564',\n",
              "       't_84th_cb2018_z014_u3_1564', 't_16th_cb2018_z014_u3_1576',\n",
              "       't_50th_cb2018_z014_u3_1576', 't_84th_cb2018_z014_u3_1576',\n",
              "       't_16th_cb2018_z014_u3_1949', 't_50th_cb2018_z014_u3_1949',\n",
              "       't_84th_cb2018_z014_u3_1949', 't_16th_cb2018_z014_u3_2449',\n",
              "       't_50th_cb2018_z014_u3_2449', 't_84th_cb2018_z014_u3_2449',\n",
              "       't_16th_cb2018_z014_u3_2732', 't_50th_cb2018_z014_u3_2732',\n",
              "       't_84th_cb2018_z014_u3_2732', 't_16th_cb2018_z014_u3_2740',\n",
              "       't_50th_cb2018_z014_u3_2740', 't_84th_cb2018_z014_u3_2740',\n",
              "       't_16th_cb2018_z014_u3_534', 't_50th_cb2018_z014_u3_534',\n",
              "       't_84th_cb2018_z014_u3_534', 't_16th_cb2018_z014_u3_417',\n",
              "       't_50th_cb2018_z014_u3_417', 't_84th_cb2018_z014_u3_417',\n",
              "       't_16th_cb2018_z014_u3_589', 't_50th_cb2018_z014_u3_589',\n",
              "       't_84th_cb2018_z014_u3_589', 't_16th_cb2018_z014_u3_11',\n",
              "       't_50th_cb2018_z014_u3_11', 't_84th_cb2018_z014_u3_11',\n",
              "       't_16th_cb2018_z014_u3_894', 't_50th_cb2018_z014_u3_894',\n",
              "       't_84th_cb2018_z014_u3_894', 't_16th_cb2018_z014_u3_90',\n",
              "       't_50th_cb2018_z014_u3_90', 't_84th_cb2018_z014_u3_90',\n",
              "       't_16th_cb2018_z014_u3_93', 't_50th_cb2018_z014_u3_93',\n",
              "       't_84th_cb2018_z014_u3_93', 't_16th_cb2018_z014_u3_531',\n",
              "       't_50th_cb2018_z014_u3_531', 't_84th_cb2018_z014_u3_531',\n",
              "       't_16th_cb2018_z014_u3_816', 't_50th_cb2018_z014_u3_816',\n",
              "       't_84th_cb2018_z014_u3_816'],\n",
              "      dtype='object')"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 28
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "w044j3aMcrsm",
        "outputId": "f0ce93bb-e123-4fc7-e1e4-d4b152621ac8"
      },
      "source": [
        "a = np.where(percentiles_columns.str.contains(\"modp020\") == True)\n",
        "percentiles_Krumholz2015 = percentiles_columns[a]\n",
        "\n",
        "a = np.where(percentiles_Krumholz2015.str.contains(\"t_\") == True)\n",
        "percentiles_Krumholz2015 = percentiles_Krumholz2015[a]\n",
        "\n",
        "percentiles_Krumholz2015"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Index(['t_16th_modp020_chabrier_MW_phi0.50_383',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_383',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_383',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_451',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_451',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_451',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_1252',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_1252',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_1252',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_1381',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_1381',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_1381',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_1564',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_1564',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_1564',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_1576',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_1576',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_1576',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_1949',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_1949',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_1949',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_2449',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_2449',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_2449',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_2732',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_2732',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_2732',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_2740',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_2740',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_2740',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_534',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_534',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_534',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_417',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_417',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_417',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_589',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_589',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_589',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_11',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_11',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_11',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_894',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_894',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_894',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_90',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_90',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_90',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_93',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_93',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_93',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_531',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_531',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_531',\n",
              "       't_16th_modp020_chabrier_MW_phi0.50_816',\n",
              "       't_50th_modp020_chabrier_MW_phi0.50_816',\n",
              "       't_84th_modp020_chabrier_MW_phi0.50_816'],\n",
              "      dtype='object')"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 29
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 163
        },
        "id": "amAJET1adQjj",
        "outputId": "138707e9-3153-4f16-883b-a037dde266e8"
      },
      "source": [
        "PDF.loc[:, CB2018]"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "error",
          "ename": "NameError",
          "evalue": "ignored",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
            "\u001b[0;32m<ipython-input-40-038e147797e2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mPDF\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mCB2018\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
            "\u001b[0;31mNameError\u001b[0m: name 'CB2018' is not defined"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BriYBosMdp-F"
      },
      "source": [
        "a = np.where(PDF_columns.str.contains(\"tpdf_modp020\") == True)\n",
        "age_Krumholz2015 = PDF_columns[a]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Wp-KxAtjdxA9",
        "outputId": "d82c0cb2-d327-4073-de65-9369990e7c56"
      },
      "source": [
        "age_Krumholz2015"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Index(['tpdf_modp020_chabrier_MW_phi0.50_ID383',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID451',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID1252',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID1381',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID1564',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID1576',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID1949',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID2449',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID2732',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID2740',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID534',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID417',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID589',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID11',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID894',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID90',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID93',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID531',\n",
              "       'tpdf_modp020_chabrier_MW_phi0.50_ID816'],\n",
              "      dtype='object')"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 33
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "id": "5a-fz_UQdyH6",
        "outputId": "2757a90f-3e47-41d9-b52a-723cd778adf7"
      },
      "source": [
        "PDF_columns[2]\n",
        "PDF_columns[6]"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            },
            "text/plain": [
              "'logt_modp020_chabrier_MW_phi0.50_ID383'"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 35
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gu5HiU_hd0Ai"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}