{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "508e45d7",
   "metadata": {},
   "source": [
    "This notebook illustrates the runaway greenhouse phenomenon\n",
    "for a grey gas.  The atmosphere consists of a one-component\n",
    "saturated condensible atmosphere, with T(p) determined by the\n",
    "Clausius-Clapeyron relation.  The surface pressure (and hence\n",
    "optical thickness of the atmosphere) increases with T.  The\n",
    "script generates OLR(Tsurf) for this case. It assumes constant\n",
    "specific absorption cross-section kappa, but more general\n",
    "cases can be treated by editing the integrand function f\n",
    "and the expression for tauInf\n",
    " \n",
    "In this case, we find the OLR by carrying out the definite\n",
    "integral giving the solution to the Schwarzschild equations.\n",
    "The integral is carried out by integrating from the top of the\n",
    "atmosphere down to the ground, using delTau = tauInf-tau as the\n",
    "integration variable.\n",
    " \n",
    "**ToDo:  Turn the OLR computation into a function that\n",
    "         returns a Curve object. That will make it easier\n",
    "         for the students to modify the script to make graphs\n",
    "         exploring the parameter dependence\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "8dd8727f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Data on section of text which this script is associated with\n",
    "Chapter = '4.**'\n",
    "Figure = '**'\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "28fcda64",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ClimateUtilities import *\n",
    "import phys\n",
    "from math import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "0464a6d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Temperature profile. The argument is the ratio of\n",
    "#pressure to the reference pressure\n",
    "#psMax/g is the available ocean mass path\n",
    "def T(pp0,Ts):\n",
    "    pp0 = max(pp0,1.e-20) #Avoids math range error at top of atmosphere \n",
    "    Tad = Ts*((pp0*p0)/psMax)**Rcp\n",
    "    Tdew = T0/(1. - RTL*log(pp0)) #Dew/frost point\n",
    "    #**ToDo: Switch to adiabat when necessary\n",
    "    #   (Haven't implemented this properly yet) This is just here for debugging\n",
    "    #   (Reason OLR comes out so small is that atm is extremely optically thick)\n",
    "    return max(Tad,Tdew)\n",
    "#Saturation pressure as a function of temperature (the inverse of above)\n",
    "#  Just needed for making the graph of OLR(T), and determining tauInf\n",
    "def pp0Fun(T):#Returns saturation vapor pressure\n",
    "    logpp0 = (1-T0/T)/RTL\n",
    "    return exp(logpp0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e1255e9",
   "metadata": {},
   "source": [
    "First, let's make some plots of the T(p) profile, first for some cases with the ocean still existing, so that the atmosphere is saturated at the surface, and then for some cases after the ocean has completely evaporated into the atmosphere, in which case there is a subsaturated layer starting from the dry surface, and extending up to altitudes where the atmosphere starts to condense."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "27f56045",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define a few necessary constants\n",
    "#Any T-p point on the phase boundary \n",
    "T0 = 300.\n",
    "p0 = phys.satvpg(T0)\n",
    "RTL = phys.water.R*T0/phys.water.L_vaporization\n",
    "#\n",
    "Rcp = 2/7\n",
    "#"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4c20687",
   "metadata": {},
   "source": [
    "First do the calculation without any constraint from mass of available ocean. We plot each surface temperature case separately, since in this case T is a function of p independent of Ts, so the curves all lay on top of each other, and are just extended to higher pressure at the bottom. In this regime, as surface temperature increases, you just add more mass at the bottom of the atmosphere, leaving the upper atmosphere temperature unchanged. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "4d002c0a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc4AAAD4CAYAAABlsga0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuuUlEQVR4nO3deXhU9dnG8e8dIKKoqMgiIOACKuLCIqIIuBdURKkoal2ASmm11apV6WtrbbUuVVu11q0gXagUENyLVkXUiiKoFBCwiKCALK4ILoB53j/OxKYxgUy2k0nuz3XNlcxZ5twThjw55/wWRQRmZmZWNnlpBzAzM8slLpxmZmZZcOE0MzPLggunmZlZFlw4zczMslA/7QA12c477xzt2rVLO4aZ5YBZs2a9HxFN085hVc+FczPatWvHzJkz045hZjlA0tK0M1j18KVaMzOzLLhwmpmZZcGF08zMLAsunGZmZllw4TQzM8tCnSqcknaXNErSxLSzmJlZbsqZwilptKTVkuYWW95X0kJJiyRdsbnXiIjFETGsapOamVltljOFExgD9C26QFI94A6gH9AROF1SR0n7SXq02KNZVYZbvBhGjoRNm6ryKGZmlracGQAhIp6T1K7Y4u7AoohYDCBpHDAgIq4DTijPcSQNB4YDtGnTpsz7PfQQXH89zJ4N48bB9tuX5+hmZlbT5dIZZ0laAe8Web4ss6xEkppIugvoLGlkSdtExD0R0S0iujVtWvbRs378Y7j7bnjySTjsMHjnnTLvamZmOSTXC6dKWBalbRwRH0TEiIjYI3NWWqmGD4cpU5Ki2b07vPJKZR/BzMzSluuFcxmwa5HnrYEVKWUB4Oij4cUXYeutoU8feOCBNNOYmVlly/XC+QrQXtJukvKBwcDDKWeiY0d4+WU48EA45RS48UaIUs+Dzcwsl+RM4ZR0PzAd2EvSMknDImITcAHwBDAfGB8R89LMWahZM3jmGRg8GC6/HM47DzZsSDuVmZlVVC61qj29lOWPA49Xc5wyadgQxo6F9u3hV7+Ct9+GiRNhxx3TTmZmZuWVM2ecuSovD375S/jTn+D55+GQQ+Ctt9JOZWZm5eXCWU3OPhueegrWrIGDD4YXXkg7kZmZlYcLZzXq3RteegmaNIGjjkou45qZWW5x4axm7dvD9Olw6KHwne/AVVe5xa2ZWS5x4UzBTjvBE0/AkCHJ/c/TT4fPP087lZmZlYULZ0ry82HUKLjhBhg/Ho44AlauTDuVmZltiQtniiS47LJkdKE5c5JGQ3PmpJ3KzMw2x4WzBjj55KSryqZNyb3Px2tkr1QzMwMXzhqjSxeYMSNpPNS/P9x2mxsNmZnVRC6cNUirVsmZ54knwoUXwvnnw8aNaacyM7OiXDhrmEaNknuel10Gd94Jxx8PH3+cdiozMyvkwlkD5eUlrW1HjYKpU5P7nosXp53KzMzAhbNGGzoU/vlPWLUqmRjbw/SZmaXPhbOGO/zw/x2m7y9/STuRmVnd5sKZAwqH6evZMxks/soroaAg7VRmZnWTC2eO2GknmDIFhg2Da69NJsj+7LO0U5mZ1T0unDkkPx/uvRduuimZEPvww+G999JOZWZWt7hw5hgJLrkEJk+GefOSYfpmz047lZlZ3eHCmaMGDEha2RYUJPc+H3kk7URmZnWDC2cO69w5GaZvn32SQnrLLR6mz8ysqrlw5riWLWHatGSg+EsugREjPEyfmVlVcuGsBbbZBiZMgJEj4Z57oF8/+OijtFOZmdVOLpy1RF4e/PrXMGYMPPccHHIILFqUdiozs9rHhbOWOecceOopWLMmaXH73HNpJzIzq11cOGuh3r3h5ZehaVM4+ujkLNTMzCqHC2ctteeeyTB9vXvDkCHJ/U8P02dmVnEunLXYjjvCP/4Bw4fD9dfDoEGwfn3aqczMcpsLZy3XoAHcdVfSx3Py5OQMdNmytFOZmeUuF846QIIf/zgZXeg//0nm9pwxI+1UZma5yYWzDjn++OS+Z8OG0KcPjBuXdiIzs9zjwlnH7Ltv0uK2Wzc4/XS46io3GjIzy4YLZx3UtGnS13PIEPjlLz23p5lZNlw466ittoJRo+A3v0nm9uzdG5YvTzuVmVnNV+cKp6TdJY2SNDHtLGmT4NJL4eGHYeFCOOggmDkz7VRmZjVbmQqnpNGSVkuaW8r6XSVNlTRf0jxJFxZZt0TSHEmvS6rQr+XSckjqK2mhpEWSrtjca0TE4ogYVpEctc0JJ8CLL0J+fnLmOX582onMzGqusp5xjgH6bmb9JuCSiNgH6AGcL6ljkfVHRMSBEdGt+I6SmknartiyPcuaQ1I94A6gH9AROF1SR0n7SXq02KPZ5t9m3bXffkkXlS5d4LTT4OqrPbenmVlJylQ4I+I54MPNrH8vIl7NfP8pMB9oVcYMfYCHJDUEkHQecFsWOboDizJnkhuAccCAiJgTEScUe6wuY6Y6qVkzePppOPts+MUvkla3n3+ediozs5ql0u9xSmoHdAZeziwK4ElJsyQNL759REwApgDjJJ0JDAVOzeKQrYB3izxfxmaKtqQmku4COksaWco2/SXd88knn2QRo3bYaqtkUPgbbkgu2fbpAytWpJ3KzKzmqNTCKWlb4AHgoohYm1ncMyK6kFxKPV9S7+L7RcSNwBfAncCJEbEum8OWsKzUi4wR8UFEjIiIPSLiulK2eSQihjdu3DiLGLWHBJddBg8+CG+8kYw09OqraacyM6sZKq1wSmpAUjTHRsSkwuURsSLzdTUwmeTSavF9ewGdMuuvyvLQy4BdizxvDfgcqRKceGLSaKhePTjsMHjggbQTmZmlr1IKpyQBo4D5EXFLkeWNChv+SGoEHAsUbxHbGbgXGAAMAXaSdE0Wh38FaC9pN0n5wGDg4Yq8H/uv/fdPGg0deCCccgpcc40bDZlZ3VbW7ij3A9OBvSQtkzQss/xxSS2BnsBZwJGZbievSzoOaA68IGk2MAN4LCKmFHv5bYBBEfFWRBQA5wBLy5ojIjYBFwBPkDRKGh8R87L6KdhmNW8OzzwDZ50FP/sZnHmmGw2ZWd2l8OlDqbp16xYzPSLA1yKSRkMjR8LBByf3QFu0SDuVWc0gaVZJXe6s9qlzIwdZ+UlwxRUwaRLMmZOMNPTaa2mnMjOrXi6clrWTT4Z//SsppIcdlkyQbWZWV7hwWrkceGDSaGi//WDgQPj1r91oyMzqBhdOK7cWLeDZZ+GMM+D//i9pPPTFF2mnMjOrWi6cViENG8Jf/wrXXgtjx8IRR8CqVWmnMjOrOi6cVmES/PSnyQAJ//530mho9uy0U5mZVQ0XTqs0AwfC889DQQH07Jl0VzEzq21cOK1SdekCr7wC++6bFNLrr3ejITOrXVw4rdLtskvSaGjw4GSwhLPO8khDZlZ7uHBaldh666SxUGGjIU9PZma1hQunVZnCRkOF05MddBB4BEMzy3UunFblBgyA6dMhPx969YL77087kZlZ+blwWrXYb79kpKGDDvrvgAkFBWmnMjPLngunVZumTeGpp+C885Ih+gYOhE8/TTuVmVl2XDitWuXnw913w223waOPwqGHwttvp53KzKzsXDit2knwwx/ClCmwbBl07w7PPZd2KjOzsnHhtNQcfXRy37NJEzjqKLj33rQTmZltmQunpap9e3jppaSIDh8OP/oRbNqUdiozs9K5cFrqdtghud958cVw++3Qrx98+GHaqczMSubCaTVCvXpw881w333J/c6DD4YFC9JOZWb2TS6cVqOcey5MnQpr1ybF8x//SDuRmdn/cuG0GufQQ5MZVnbfHU44AW65xTOsmFnN4cJpNVKbNvDCC3DyyXDJJTB0KHz5ZdqpzMxcOK0Ga9QIxo+HX/wCxoyBI4+EVavSTmVmdZ0Lp9VoeXlw1VUwYQK89loy1u1rr6WdyszqMhdOywmnnAL/+lfy/WGHwcSJ6eYxs7qrzhVOSbtLGiXJv3pzTOfOSaOhAw6AQYOSS7ieYcXMqluFCqek0ZJWS5q7mW2WSJoj6XVJFZrGuLTjSeoraaGkRZKu2NxrRMTiiBhWkRyWnubNk+4q554LV18Np54K69enncrM6pKKnnGOAfqWYbsjIuLAiOhWfIWkZpK2K7Zsz7IeT1I94A6gH9AROF1Sx8y6/SQ9WuzRrAx5rQbbaisYPRpuugkmT04u3b7zTtqpzKyuqFDhjIjngIoOjtYHeEhSQwBJ5wG3ZXG87sCizJnkBmAcMCCz/ZyIOKHYY/WWAknqL+meTz75pAJvy6qSlHRTefRRWLw4aTT04otppzKzuqA67nEG8KSkWZKGf2NlxARgCjBO0pnAUODULF6/FfBukefLMstKJKmJpLuAzpJGlhg44pGIGN64ceMsYlga+vVLBonffns44oik24qZWVWqXw3H6BkRKzKXSP8paUHmzPFrEXGjpHHAncAeEbEui9dXCctKHWcmIj4ARmTx+lbD7bMPvPwynHYaDBkC//433Hgj1K+OT7eZ1TlVfsYZESsyX1cDk0kurf4PSb2ATpn1V2V5iGXArkWetwZWlCus5ayddkrGtf3Rj+C3v4XjjvMMK2ZWNaq0cEpqVNjwR1Ij4FigeIvYzsC9JPclhwA7Sbomi8O8ArSXtJukfGAw8HBl5LfcUr8+3HorjBoFzz4L3bvDvHlppzKz2qai3VHuB6YDe0laJmlYZvnjkloCzYEXJM0GZgCPRcSUYi+zDTAoIt6KiALgHGBpWY8XEZuAC4AngPnA+Ijwr8s6bOjQpHCuWwc9esDD/jPKzCqRwtNOlKpbt24xc2aFup5aipYtSwaJnzULfvUr+OlPk9a4ZlVB0qySutxZ7VPnRg6yuqN162RS7DPOgCuvTBoPebAEM6soF06r1bbeGv7yF/jNb+CBB6BnT1ha4o0AM7OyceG0Wk+CSy+Fxx6DJUugW7fkTNTMrDxcOK3O6NsXZsyAJk3gqKPgrrvSTmRmuchdxK1O6dAhGSzhjDPg+9+H11+H226D/Py0k1ltNGvWrGb169f/I0k/dZ+o5IYCYO6mTZu+27Vr1xKHaHXhtDqnceOki8qVV8L118MbbyTzezbz8P9WyerXr//HFi1a7NO0adOP8vLy3IUhBxQUFGjNmjUdV65c+UfgxJK28V9AVifVqwfXXQd/+1syx+dBB8Frr6WdymqhTk2bNl3ropk78vLyomnTpp+QXCUoeZtqzGNW45x+OrzwQjIhds+e8Pe/p53Iapk8F83ck/k3K7U+unBande1K8ycCV26wODB8H//lxRSM7OSuHCaAc2bw9NPw3e/C7/+NZx0Eqxdm3Yqs4pZuXJlvb333rvj3nvv3XHnnXc+oFmzZvsXPv/iiy82O47W9773vda77bbbvh06dOh4zDHH7PH+++/XA/jyyy81cODAdh06dOi4++677zty5MgWhfs8//zz23To0KFjmzZtOp177rm7FpTyF+jIkSNbtGnTplO7du06PfDAA9tvaf/PP/9cxx9//O5t2rTptP/++++9cOHCr5vz3X777U3atm3bqW3btp1uv/32JoXLFyxYkL///vvv3bZt207HH3/87oXvt6CggHPPPXfXNm3adOrQoUPHF154YZtsf64unGYZW20F99wDv/89PP54Ms7tf/6Tdiqz8mvRosVXCxYseGPBggVvnH322WtGjBixqvB5w4YNN3sJ+Vvf+tbaN998c96bb775xp577vnFz372sxYA9913344bNmzIe/PNN9+YPXv2/D//+c9NCwvZD37wg7Z/+MMfli5ZsmTu4sWLG06cOHH74q87a9ashpMmTdpp4cKF86ZMmfLmRRdd1GbTpk1sbv9bb71158aNG29655135l5wwQWrLr744tYAq1atqnfDDTe0nDFjxvyZM2fOv+GGG1quWbOmHsDFF1/c+oILLli1dOnSuY0bN95066237gwwYcKExosXL264ZMmSuXfeeefSH/zgB22y/bm6cJoVIcH558NTT8Hq1ckMK08+mXYqs+o3cODAtQ0aNADgkEMOWb98+fJ8AEl89tlneRs3bmT9+vVq0KBB7LDDDl8tXbq0wbp16/KOPvro9Xl5eZx55pkfPPjggzsWf92JEyfuMHDgwA+33nrr2HvvvTe0bdv2y2effbbR5vZ/9NFHdxg6dOgHAEOGDPnoxRdf3K6goIAHH3ywce/evdc2b978q6ZNm37Vu3fvtZMmTWpcUFDA9OnTtxsyZMhHAEOHDv3gkUce2QHgoYce2uHMM8/8IC8vj6OOOmr92rVr6y9durRBNj8bd0cxK8HhhyetbQcMgH79kiH7fvxjDxJv5Td0KLvOnUvWlwU3p1MnPhs9mnfLs2/Xrl33Wr9+fb3iy6+//vp3TzrppE+LLhszZszOp5xyyocA55577kePPPLIDs2aNTvgiy++yPvVr371bvPmzb967rnnttpll102Fu7Ttm3bDe+99943CtLy5cvze/Tosa7wecuWLTe8++67+fn5+VHa/qtWrcrfbbfdNgA0aNCAbbfd9qtVq1bVX758eYPWrVtvKNynVatWG5YvX95g1apV9bfbbruvCgt/u3btNqxatSof4L333mvQrl27r/fZZZddNixdurRB27Ztvz72lrhwmpVit93gxRfhnHPgkktg9my4+25o2DDtZGYVN2vWrIVl2e7yyy9vUa9evRgxYsSHANOmTdsmLy8vVq5c+e/333+/Xs+ePfc+7rjj1pY005ZK+EuzlO1ic/tnu09p22/mtb6xbHNcOM02Y9ttYcIEuPZa+PnPYcECmDwZWrZMO5nlmvKeGVaVspxx3n777U2eeOKJHZ5//vk38/KSO3t/+ctfmnzrW9/6ZKuttopWrVptOuigg9a9+OKLjY455ph1Rc8wly5dmt+iRYtvnMW1bt16w7vvvvt1454VK1bkt27demO7du02lrZ/ixYtNrz99tv5e+yxx8aNGzeybt26es2aNfuqdevWG6dNm7Zd4T7Lly/P79Onz6ctWrTY9Omnn9bbuHEjDRo0YMmSJfnNmjXbCNCyZcuNS5Ys+fr47733Xn6bNm3KfLYJvsdptkV5efCznyUF8403kkHiX3op7VRmFTNr1qyFhQ2Fij4Ki+bEiRO3/93vftfi8ccfX7Tddtt93Ty2TZs2G6ZOnbp9QUEBa9euzXv11Vcb7bfffl+0bdt2Y6NGjQqefvrpRgUFBYwdO7bJgAEDPi5+3G9/+9sfT5o0aafPP/9cCxYsyF+yZEnDww8/fP3m9j/++OM/Hj16dBNIGicdcsghn+bl5XHSSSd9Mm3atO3XrFlTb82aNfWmTZu2/UknnfRJXl4ePXr0+PS+++7bEWD06NFNTjjhhI8BTjzxxI/Hjh3bpKCggKeffrrRdttt91U2l2nBhdOszE46CaZPTy7V9ukDf/pT2onMqs7FF1/cZv369fWOPPLIDnvvvXfHM844ow3AZZddtnr9+vV5HTp02Ldz5877nHHGGe8ffPDBnwP84Q9/WDpixIh2bdu27dSuXbsvBw0a9AnA2LFjG1900UUtAbp16/bFSSed9GGHDh327du3b4dbbrllaf36ycXP0va/8MIL3//oo4/qt2nTptPtt9/e4qabbloG0Lx5869+8pOfrOjates+Xbt23eeyyy5b0bx5868Abr755mW33357izZt2nT66KOP6l944YXvA5x66qmftG3b9su2bdt2+v73v9/2jjvuyHqiQZV0vdcS3bp1i5kzZ6Ydw2qYDz6AU0+FZ56Biy5KGg7V902POk/SrIjoVnTZ7NmzlxxwwAHvp5XJym/27Nk7H3DAAe1KWuczTrMsNWkCTzwBP/oR/O53SavbDz9MO5WZVRcXTrNyqF8fbr0VRo1KJsU+6CCYMyftVGZWHepc4ZS0u6RRkiamncVy39ChMG0afP45HHJIMj2ZWREFBQUF7v2bYzL/ZqWOWF2mwilptKTVkuaWsG4vSa8XeayVdFGR9UskzcmsK/cNwy1k6CtpoaRFkq7Y3OtExOKIGFbeHGbF9egBs2bB/vvDoEHJIPFffZV2Kqsh5q5Zs6axi2fuyMzH2Rj4Rq0pVNYmDWOA3wN/Lr4iIhYCBwJIqgcsByYX2+yIiCjxBrmkZsDnEfFpkWV7RsSismTIHPMO4BhgGfCKpIeBesB1xV5jaESUOKO3WUXssgtMnQoXXJAMEj97Nvz1r7DDDmknszRt2rTpuytXrvzjypUrO1EHr/DlqAJg7qZNm75b2gZlKpwR8ZykdmXY9CjgrYjIpnlvH+D7ko6LiC8knQecDBxXxgzdgUURsRhA0jhgQERcB5yQRY6vSeoP9N9zzz3Ls7vVUYWDxHftCj/8IRx8MDz4IOyzT9rJLC1du3ZdDZyYdg6rXJX9F9Bg4P5iywJ4UtIsScOL7xARE4ApwDhJZwJDgVOzOGYr+J8ROZZllpVIUhNJdwGdJY0saZuIeCQihjdu3DiLGGbJWLYjRiRdVT7+OCmeDz+cdiozq0yVVjgl5ZP8ZTWh2KqeEdEF6AecL6l38X0j4kbgC+BO4MSIWFd8m80duoRlpXZOjYgPImJEROyROSs1q3S9eiWTY++1VzJQ/C9/6cmxzWqLyjzj7Ae8GhGrii6MiBWZr6tJ7n12L76jpF5Ap8z6q7I87jJg1yLPWwMrsnwNs0q3665JV5Wzz4arroJvfxs+/XTL+5lZzVaZhfN0il2mldRI0naF3wPHUqylkqTOwL3AAGAIsJOka7I47itAe0m7Zc56BwO+OGY1wtZbw5gxyUAJjzziybHNaoOydke5H5gO7CVpmaRhmeWPS2opaRuSVq2Tiu3aHHhB0mxgBvBYREwpts02wKCIeCsiCoBzgG80LiotQ0RsAi4AngDmA+MjYl5Z3pdZdZDgwguTCbFXrUoGS/jHP9JOZWbl5bFqN8Nj1VplW7IkGSz+3/9Ouq1cfrknx64tShqr1mon9ysyq0bt2iWTY592GowcCYMHw/r1aacys2y4cJpVs222gb/9DW68MRmi79BD4e23005lZmXlwmmWAgl+8hN4/HF4551kcuynn047lZmVhQunWYq+9S145ZVkyL5jj4VbbgE3OzCr2Vw4zVK2557w0ktJo6FLLoGzzkpmWzGzmsmF06wG2Hbb5H7nNdck9z8POyy5hGtmNY8Lp1kNISVTkj38MCxalNz3nDYt7VRmVpwLp1kNc8IJMGMG7LQTHH00/P73vu9pVpO4cJrVQHvtBS+/DP36JVOUffe78MUXaacyM3DhNKuxGjdO5vP8+c9h9Gg4/HBYvjztVGbmwmlWg+XlwdVXw6RJMG9ect/zxRfTTmVWt7lwmuWAk09Ouqxsu21y5nnvvWknMqu7XDjNcsS++yaNho46CoYPh+99D778Mu1UZnWPC6dZDtlxR3j00WSA+HvuSc4+V3jadrNq5cJplmPq1UumJJs4EebMgS5d4IUX0k5lVne4cJrlqG9/O+mysv32cMQR8Ic/uL+nWXWok4VT0u6SRkmamHYWs4oovO/Zty+cfz4MHer+nmZVrcoLp6RdJU2VNF/SPEkXVuC1RktaLWluCev6SlooaZGkKzb3OhGxOCKGlTeHWU2yww7w0ENw1VUwZgz06uVxbs2qUnWccW4CLomIfYAewPmSOhbdQFIzSdsVW7ZnCa81BuhbfKGkesAdQD+gI3C6pI6S9pP0aLFHs8p5W2Y1R14e/OIXSQF9803o2hWmTk07lVntVOWFMyLei4hXM99/CswHWhXbrA/wkKSGAJLOA24r4bWeAz4s4TDdgUWZM8kNwDhgQETMiYgTij1WbymzpP6S7vnkk0+yeatmqTvxxOTSbdOmcMwx8Nvf+r6nWWWr1nucktoBnYGXiy6PiAnAFGCcpDOBocCpWbx0K+DdIs+X8c3iXDRHE0l3AZ0ljSy+PiIeiYjhjRs3ziKCWc1QOM7tiSfCxRfDd74Dn32Wdiqz2qPaCqekbYEHgIsiYm3x9RFxI/AFcCdwYkSsy+blS1hW6t/ZEfFBRIyIiD0i4rosjmOWE7bbLumucu21cP/9cOih8Pbbaacyqx2qpXBKakBSNMdGxKRStukFdAImA1dleYhlwK5FnrcG3C3c6rS8PPjpT+Gxx2Dp0mSc2yefTDuVWe6rjla1AkYB8yPillK26QzcCwwAhgA7Sbomi8O8ArSXtJukfGAw8HDFkpvVDv36wcyZ0KpV8v0NN/i+p1lFVMcZZ0/gLOBISa9nHscV22YbYFBEvBURBcA5wNLiLyTpfmA6sJekZZKGAUTEJuAC4AmSxkfjI2Je1b0ls9yyxx4wfToMGgRXXAGnngrrsrkZYmZfU/hPz1J169YtZs6cmXYMs0oTATffDJdfDvvsA5MnQ/v2aaeqHSTNiohuaeewqlcnRw4yq6skuPRSeOIJWLkSDjoouQdqZmXnwmlWBx19dHLfc/fdoX9/+OUvoaAg7VRmucGF06yOatcO/vWvpJ/nVVclk2V7zA+zLXPhNKvDtt4a/vQnuO225JLtwQfD/PlppzKr2Vw4zeo4CX74Q3j6afjoI+jePWk0ZGYlc+E0MwD69IFZs6BjRxg4EK68Er76Ku1UZjWPC6eZfa11a5g2DYYNS4br698/OQs1s/9y4TSz/9GwIdx7L9x1Fzz1VNJlZe43ZsA1q7tcOM3sGyT43veSs8/PPksaDY0fn3Yqs5rBhdPMSnXIIcl9zwMPhNNOg5/8BDZtSjuVWbpcOM1ss3bZBaZOhR/8AG66CY49FlZvcTp4s9rLhdPMtig/H+64I+nzOX06dO2aTJZtVhe5cJpZmZ19Nrz4ItSvD717w913e4oyq3vqXOGUtLukUZImpp3FLBd17pzc9zzySBgxIum68vnnaacyqz5bLJySdpU0VdJ8SfMkXZjNNpKWSJqTmYezQnN0SRotabWkucWW95W0UNIiSVds7jUiYnFEDKtIDrO6bqed4NFH4Wc/g/vug8MOgyVL0k5lVj3Kcsa5CbgkIvYBegDnS+qY5TZHRMSBJc1VJ6mZpO2KLduzlCxjgL7Ftq0H3AH0AzoCp0vqKGk/SY8WezQrw/s1szKoVy+ZVeXhh+Gtt5L7nk8+mXYqs6q3xcIZEe9FxKuZ7z8F5gOtst1mM/oAD0lqCCDpPOC2UrI8B3xYbHF3YFHmTHIDMA4YEBFzIuKEYo8ytQWU1F/SPZ94qgizLerfH155BVq2hL594brrfN/Tares7nFKagd0BkptT1fCNgE8KWmWpOHFt4+ICcAUYJykM4GhwKlZxGoFvFvk+TI2U7QlNZF0F9BZ0siStomIRyJieOPGjbOIYVZ3tW8PL70EgwfDT3+ajHW7dm3aqcyqRv2ybihpW+AB4KKIKPG/RCnb9IyIFZnLpP+UtCBz5vi1iLhR0jjgTmCPiFiXxXtQCctK/Xs3Ij4ARmTx+mZWBo0awdixyShDl1ySDNU3eXIyaLxZbVKmM05JDUgK4tiImJTNNhGxIvN1NTCZ5NJq8X17AZ0y66/K8j0sA3Yt8rw1sCLL1zCzSiDBhRfCM88kk2J37w4TJqSdyqxylaVVrYBRwPyIuCWbbSQ1Kmz4I6kRcCxQvEVsZ+BeYAAwBNhJ0jVZvIdXgPaSdpOUDwwGHs5ifzOrZL17w6uvwv77w6mnwqWXeqg+qz3KcsbZEzgLODLTpeR1SccBSHpcUsvNbNMceEHSbGAG8FhETCn2+tsAgyLirYgoAM4BlpYURNL9wHRgL0nLJA2LiE3ABcATJI2SxkfEvOx+DGZW2Vq2hGefhfPPh5tvhmOO8VB9Vjso3PytVN26dYuZMyvU9dTMgL/8BYYPhyZNYOJE6NEj7USVT9KskrrcWe1T50YOMrPqd9ZZyRi3+fnJZdy77nKXFctdLpxmVi0OPBBmzoSjj4bvfx+GDvVQfZabXDjNrNoUDtX385/DmDEeqs9ykwunmVWrvDy4+moP1We5y4XTzFLRv39y6bZwqL5rr4WCgrRTmW2ZC6eZpWbPPf87VN+VVyZD9XmIaKvpXDjNLFWFQ/X97nfw2GPJUH1z525xN7PUuHCaWeqKDtW3dm0y3u3f/552KrOSuXCaWY3Rq1cyVN+BByaXby++GDZuTDuV2f9y4TSzGqVlS5g6FS64AH7726Tf58qVaacy+686Vzgl7S5plKSJaWcxs5Ll58PttydD9b3yCnTpAi+8kHYqs0S5C6ekXSVNlTRf0jxJF5ay3RJJczIDv1do4FdJoyWtllR8hpW+khZKWiTpis29RkQsjohhFclhZtXjO99JWt02agSHH56cgXqoPktbRc44NwGXRMQ+QA/gfEmlTVl7REQcWNIAyJKaFU49VmTZnqW8zhigb7Ft6wF3AP2AjsDphTkk7Sfp0WKPZlm8RzNL2f77J/09+/dP7nkOHgyffpp2KqvLyl04I+K9iHg18/2nJFN6tSrHS/UBHpLUEEDSecBtpRzzOeDDYou7A4syZ5IbgHEkc3sSEXMi4oRijy1ObCSpv6R7PnGHMrMaoXFjmDQJbrghmV3l4INh/vy0U1ldVSn3OCW1AzoDL5ewOoAnJc2SNPwbKyMmAFOAcZLOBIYCp2Zx+FbAu0WeL2MzBVxSE0l3AZ0ljSxpm4h4JCKGN27cOIsYZlaVJLjsMnjqKfjgA+jeHcaPTzuV1UUVLpyStgUeAC6KiLUlbNIzIrqQXEo9X1Lv4htExI3AF8CdwIkRsS6bCCUsK/UuSER8EBEjImKPiLgui+OYWQ1wxBFJl5X99oPTToMf/9hdVqx6VahwSmpAUjTHRsSkkraJiBWZr6uBySSXVou/Ti+gU2b9VVnGWAbsWuR5a2BFlq9hZjmkVSt49ln40Y+SEYeOOAJW+H+9VZOKtKoVMAqYHxG3lLJNo8KGP5IaAccCxVvEdgbuJbkvOQTYSdI1WUR5BWgvaTdJ+cBg4OFs34+Z5Zb8fLj1Vrj/fnj99aTLyrRpaaeyuqAiZ5w9gbOAIzNdTV6XdByApMcltQSaAy9Img3MAB6LiCnFXmcbYFBEvBURBcA5wNKSDijpfmA6sJekZZKGRcQm4ALgCZIGSuMjYl4F3peZ5ZDBg+Hll2GHHeCoo+Cmm9xlxaqWwp+wUnXr1i1mzqxQ11MzqyZr18KwYUmr24ED4b77YPvtq+/4kmaV1OXOap86N3KQmdVO22+ftLK9+WZ46CHo1s2zrFjVcOE0s1pDSgZJeOaZZJCEgw+Gv/0t7VRW27hwmlmt07t30mWlSxc480z44Q9hw4a0U1lt4cJpZrXSLrskZ54XXwy//z306QPLlqWdymoDF04zq7UaNEjueY4fn9zv7NwZnn467VSW61w4zazWGzQomZ6saVM49li47jooKEg7leUqF04zqxP23htmzEiK6E9/CiefDB9/nHYqy0UunGZWZ2y7bTLS0K23wuOPJ11WZs9OO5XlmjpXOCXtLmmUpIlpZzGz6iclY9w++yx8/jn06AF/+lPaqSyXbLFwStpV0lRJ8yXNk3RhsfV7FRly73VJayVdVGT9EklzMuvKPQyPpNGSVkv6RpdmSX0lLZS0SNIVm3udzLydw8qbw8xqh549ky4rPXrAuefCiBHw5Zdpp7JcUJYzzk3AJRGxD9CDZGqwjoUrI2JhRBwYEQcCXYHPSGY5KeqIzDbfGI5KUrPCgeCLLNuzhBxjgL4l7F8PuINk2rKOwOmSOkraT9KjxR7NyvB+zayOaN4c/vlPuPxyuPtu6NUL3nkn7VRW022xcEbEexHxaub7T0kGUi9touijgLciosRB2kvRB3hIUkMASecBt5WQ4zngwxL27w4sypxJbgDGAQMiYk5EnFDssTqLXGZWB9SvD9dfD5Mnw8KFyaAJTz6ZdiqrybK6xympHdAZeLmUTQYD9xdbFsCTkmZJGl58h4iYAEwBxkk6ExgKnJpFrFbAu0WeL6P0wo6kJpLuAjpLGlnKNv0l3fPJJ59kEcPMctlJJ8HMmcnACX37Jv0/zUpS5sIpaVuSSasvioi1JazPB04EJhRb1TMiupBcSj1fUu/i+0bEjcAXwJ3AiRGxruxvAZWwrNQpXyLig4gYERF7RMR1pWzzSEQMb9y4cRYxzCzXtW8PL70E3/kOdOiQdhqrqcpUOCU1ICmaYyNiUimb9QNejYhVRRdGxIrM19Uk9z67l/D6vYBOmfVXlTl9Yhmwa5HnrQHPBW9m5dKoEfz5z9C/f9pJrKYqS6taAaOA+RFxy2Y2PZ1il2klNSps+COpEXAsMLfYNp2Be4EBwBBgJ0nXZPEeXgHaS9otc9Y7GHg4i/3NzMzKrCxnnD2Bs4Aji3Q5OQ5A0uOSWkraBjgGKH422hx4QdJsYAbwWERMKbbNNsCgiHgrIgqAc4BvNC6SdD8wHdhL0jJJwwAiYhNwAfAEScOl8RExr0zv3szMLEuKKPV2YJ3XrVu3mDmz3F1PzawOkTSrpC53VvvUuZGDzMzMKsKF08zMLAsunGZmZllw4TQzM8uCC6eZmVkW3Kp2MyStoYSuMVuwM/B+FcQpj5qSpabkAGcpSU3JAbmdpW1ENK2qMFZzuHBWMkkza0qT9JqSpabkAGepyTnAWSw3+FKtmZlZFlw4zczMsuDCWfnuSTtAETUlS03JAc5SkpqSA5zFcoDvcZqZmWXBZ5xmZmZZcOE0MzPLggtnFiTtKmmqpPmS5km6sMi6H0pamFl+Y5Hl+0uanlk+R1LDNLJIaiDpT5kM8yWNrIwcm8si6e9FpqJbIun1IvuMlLQok/NbaeSQdIykWZmfySxJR1ZGjvJkKbJfG0nrJF2aZpaq+NyW498njc/sgZJeymSZKal7kX0q/TNrOSoi/CjjA9gF6JL5fjvgTaAjcATwFLBVZl2zzNf6wL+BAzLPmwD1UspyBjAu8/02wBKgXVVmKbbNzcDPM993BGYDWwG7AW9Vxs+lHDk6Ay0z33cCllf1Z6W0LEWWPQBMAC5NK0tVfW7LkaPaP7PAk0C/zPLjgGer8jPrR24+6mNlFhHvAe9lvv9U0nygFXAecH1EfJlZtzqzy7HAvyNidmb5BylmCaCRpPrA1sAGYG0VZ3kDQJKAU4HCM7oBJL8QvwTelrQI6E4yUXm15YiI14rsPg9oKGmrwp9ddWbJLDsJWAysr+jxK5ilSj635ciRxmc2gO0zmzUGVmS+r5LPrOUmX6otJ0ntSM5YXgY6AL0kvSxpmqSDMpt1AELSE5JelXRZilkmkvxCfg94B7gpIj6s4iyFegGrIuI/meetgHeLrF+WWVbdOYr6NvBaZRTN8mSR1Ai4HLi6so+fbRaq4XNbxhxpfGYvAn4j6V3gJqDw8nCVf2Ytd/iMsxwkbUtySe2iiFib+Yt4R6AHcBAwXtLuJD/fwzLLPgOeVjJL/NMpZOkOfAW0zKx/XtJTEbG4qrIUWXU6cH/RTUvYvdL6RWWRo3D7fYEbSM60KlUWWa4GfhsR65ITr8qXRZYq/dxmkaPaP7OSrgF+HBEPSDoVGAUcTRV/Zi23uHBmSVIDkv9oYyNiUmbxMmBSRAQwQ1IByQDRy4BpEfF+Zt/HgS5AZf0CyibLGcCUiNgIrJb0L6AbyaXBqspCppAPBLoW2XwZsGuR56357yWx6syBpNbAZODsiHirMjKUM8vBwClKGnPtABRI+iIifp9Clir73GaZI43P7DlAYUO7CcAfM99X2WfWco8v1WYhcw9mFDA/Im4psupBMvdlJHUA8klmVXgC2F/SNplfDH3I3M9JIcs7wJFKNCI5I11QxVkg+Wt9QUQsK7LsYWCwpK0k7Qa0B2ZUdw5JOwCPASMj4l8VPX5FskREr4hoFxHtgN8Bv67Eopntv0+VfG7LkSONz+wKkvcLyf+jwsvGVfKZtRyVduukXHqQXL4KkhaHr2cex5EUp78Cc4FXgSOL7PMdkoYnc4Eb08oCbEvyF/Q8kl+CP6nqLJl1Y4ARJezzfyQtExeSacVY3TmAK0nuob1e5NEsrZ9JkX1/QeW2qi3Pv0+lf27L8e9T7Z/ZzPJZJC1oXwa6VuVn1o/cfHjIPTMzsyz4Uq2ZmVkWXDjNzMyy4MJpZmaWBRdOMzOzLLhwmpmZZcGF08zMLAsunGZmZln4fz75XWeZdMDkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAD4CAYAAABfYrnHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhGUlEQVR4nO3deXhU1eHG8e/JJnvY15gEVEyBgggoogIiZZHFKNaquAAqUkRRqBa0P5fWjVpBRaGggkvRIosg1B0VNxRJLYoIFJEIAWLUikDBBHJ+f5yJpjEJIcnMmeX9PE8ekjtzM6+X4Jt759xzjLUWERGRSBDnO4CIiEhFqbRERCRiqLRERCRiqLRERCRiqLRERCRiJPgOUJ7GjRvb9PR03zFEJAJkZWV9ba1t4juHBFdYlpYxZggw5Nhjj2XNmjW+44hIBDDGZPvOIMEXlpcHrbXLrLWjk5OTfUcREZEwEpalJSIiUhqVloiIRAyVloiIRAyVloiIRIyQlZYxpo0x5jFjzMJQvaaIiESXCpWWMWaOMeYrY8y6EtsHGGM2GmM2G2Mmlfc9rLVbrLWXVyWsiIjEtoqeaT0ODCi+wRgTDzwMDATaARcaY9oZY35pjFle4qNptaYuxbRp8M47wX4VERHxqUKlZa19C/i2xOaTgM2BM6h84O/A2dbaT6y1g0t8fFXRQMaY0caYNcaYNXl5eRXa57//hVmz4Fe/guefr+griYhIpKnKe1qtgG3Fvt4e2FYqY0wjY8xfgc7GmMllPc9aO9ta29Va27VJk4rNyFKrljvL6tgRzjkH5syp4H+BiIhElKpM42RK2VbmMsjW2m+AMVV4vXI1bgwrVsB558Hll0NuLkyaBKa0lCIiEpGqcqa1HTi62NcpwI6qxamaOnXc5cHhw+Gmm+D666Gw0GciERGpTlU50/oQOM4Y0xrIAS4ALqqWVFWQlARPPglNm7rBGV99BY8/7raLiEhkq+iQ92eAVcDxxpjtxpjLrbUHgXHAy8BnwLPW2k+DF7Xi4uLgvvtgyhR45hkYPBj27PGdSkREqqpCZ1rW2gvL2P4C8EK1JqomxsCNN7ozriuugD594IUXoIJjO0REJAxF/TROI0bAkiXw6adw6qmwdavnQCIiUmlRX1rgLg++9hp8/TX06AEff+w7kYiIVEZMlBa4snr7bfd+V8+e8NZbvhOJiMiRipnSAmjfHt57D1q0gH793GVDERGJHDFVWgCpqW72jBNOgGHD4NFHfScSEZGKirnSAmjUyM2e0b8/XHkl3Hkn2DLn8hARkXARk6UFULs2LF0Kl1wCf/gDXHutZs8QEQl3VZkRI+IlJrrZMpo1g7/8xc2e8eSTcNRRvpOJiEhpYrq0wI0mvPdeaN4cfvc7V1xLlkBysu9kIiJSUsxeHixp4kSYNw/efRdOPx1ycnwnEhGRklRaxVx0Ebz4ops145RTYP1634lERKQ4lVYJZ57pbjwuKHDTPr3zju9EIiJSRKVVihNOgFWr3ACNvn1h8WLfiUREBFRaZUpPd+9vdeniVkN++GHfiURERKVVjkaN3ES7Q4fCuHFuNWTdhCwi4o9K6zBq1oRFi2DMGLj7brfUSX6+71QiIrEp5u/Tqoj4eJgxA1JS3OwZu3bBwoVQt67vZCIisUVnWhVkDNx8M8yd6+Yt7NXLlZeIiISOSusIjRgBy5bBpk3uXq5Nm3wnEhGJHSqtShg4EN54A/btc4tLvv++70QiIrFBpVVJ3bq5e7nq14c+fdzZl4iIBJdKqwqOOcathNy+PWRmwiOP+E4kIhLdVFpV1LSpu1TYvz+MHg233aZ7uUREgkWlVQ3q1HELSo4aBbff7lZDLijwnUpEJProPq1qkpgIjz4KrVrBn/7kljZ59lndyyUiUp10plWNjIE//tG9t/Xqq9CzJ+zY4TuViEj0UGkFwRVXwPLlsHkzdO8On37qO5GISHRQaQXJgAFuXa6DB926XG+84TuRiEjkU2kFUefO7sbjlBQ3uvBvf/OdSEQksqm0giw11a1+fNppcMklcOedGhIvIlJZKq0QqF8fXnoJLr7YzRJ/1VXusqGIiBwZDXkPkaQkePJJSEtzZ1vbtmlIvIjIkdKZVggZA3fcAbNnuyHxvXppSLyIyJFQaXlw5ZVugt1//9stb6Ih8SIiFaPS8mTgQDckPj9fQ+JFRCpKpeVR0ZD4Vq3ckPh583wnEhEJbyotz9LS4N133dnWxRfDXXdpSLyISFlUWmGgaEj88OFw880aEi8iUhYNeQ8TRx0FTz3lzrzuugu2b4f58zUkXkSkOJ1phRFj3D1cs2bBK69oSLyISEkqrTA0evRPQ+JPPhnWrvWdSEQkPKi0wtTAgW7OQmvdvIUvvug7kYiIfyqtMNapE3zwARx7LAweDDNn+k4kIuKXSivMtWoFb7/tzrzGjoWJE+HQId+pRET8UGlFgDp1YOlSuOYamDoVzjsP9u3znUpEJPRUWhEiPh4efBAeeACefx5694Zdu3ynEhEJLZVWhLn2WliyBNavdyML163znUhEJHRUWhFoyBD3PldBgZv+6ZVXfCcSEQkNlVaEOvFEN7IwPR3OOgseecR3IhGR4FNpRbCjj3b3cvXr525I/v3vobDQdyoRkeBRaUW4unXdwIzf/hb+/Gc4/3zYv993KhGR4NCEuVEgIQEefhiOO87dx7V9uxsi36yZ72QiItVLZ1pRwhi4/npYvBg+/hi6d3cjDEVEoklYlpYxZogxZvbu3bt9R4k4mZmwcqW7RNijB6xY4TuRiEj1CcvSstYus9aOTk5O9h0lInXr5kYWpqTAgAEwZ47vRCIi1SMsS0uqLi0N3n0X+vSByy+HSZM0slBEIp9KK4olJ8Py5TBmDEyZojkLRSTyqbSiXGIizJjh5ixcuhROPx1ycnynEhGpHJVWDDDGzVn4/PNuNeSTToKsLN+pRESOnEorhgwaBO+9586+evZ0w+NFRCKJSivG/PKXbmRhx44wbBjccw9Y6zuViEjFqLRiULNm8MYbcOGFMHkyjBoF+fm+U4mIHJ6mcYpRNWrAvHmQkQG33gpbtsCiRdC4se9kIiJl05lWDDMGbrkFnnnGXTLs3h02bPCdSkSkbCot4YIL4M03Yc8eV1yvveY7kYhI6VRaAriyWr0aUlPd1E9//avvRCIiP6fSkh+lpblFJfv3d+tzXX89HDrkO5WIyE9UWvI/6tVzNyFfdx3cfz+cfba7bCgiEg5UWvIz8fEwbRrMnAkvvQSnngrZ2b5TiYiotKQcY8a40vrySzf10/vv+04kIrFOpSXl6tvXlVXdutC7Nzz9tO9EIhLLVFpyWBkZP93HNXw43Hyz1uYSET9UWlIhjRrBK6/AlVfCXXfBuefC3r2+U4lIrFFpSYUlJcGsWfDgg7BsmQZoiEjoqbTkiBgD11wDL77oCqtbN3dvl4hIKKi0pFL69XPvczVoAH36wNy5vhOJSCxQaUmlHX+8G1nYu7db3mTiRM2gISLBpdKSKmnQAF54Aa69FqZOhSFDYPdu36lEJFqptKTKEhLggQfcII1XX3VD4zdv9p1KRKKRSkuqzejRrrS++srNoPH6674TiUi0UWlJterdGz78EFq2dIM1Zs70nUhEoolKS6pdmzbw3nswcCCMHQtXXw0FBb5TiUg0UGlJUNSrB0uWwI03wowZbmHJb7/1nUpEIl2C7wASveLjYcoUaN/eTf908slura5f/MJ3MokVWVlZTRMSEh4FOqBf0iNFIbDu4MGDV3Tp0uWrkg+qtCToLr0UjjsOMjPdyML5892Zl0iwJSQkPNq8efNfNGnS5D9xcXHWdx45vMLCQpOXl9du165djwJDSz6u3zwkJE45xQ3QaNMGBg1yi0xa/S9Egq9DkyZNvldhRY64uDjbpEmT3biz458/HuI8EsNSU908hZmZMGECjBwJBw74TiVRLk6FFXkCf2el9pNKS0Kqdm1YsABuvx2eeAJ69YIdO3ynEpFIEZalZYwZYoyZvVvzAUWluDi45RZ47jlYvx66dnVzGIpEk127dsVnZGS0y8jIaNe4ceNOTZs27Vj09YEDB0x5+44fP75l27Zt22VkZLQ79dRTj9u6dWti0WOTJ09unpqa2iE9Pb3DokWL6hVtf/vtt2u1bdu2XWpqaocRI0YcXVjGSq1Huv/+/fvNoEGD2qSmpnbo2LFjxsaNG5OK9pk+fXqjtLS0DmlpaR2mT5/eqGj7hg0bkjp27JiRlpbWYdCgQW2K/nsLCwsZMWLE0ampqR3atm3b7p133ql1pMc1LEvLWrvMWjs6OTnZdxQJosxMWLUKatZ0Z1yPP+47kUj1ad68+aENGzas37Bhw/pLL700b8yYMblFX9eoUaPcS5a33nrrrk2bNq3fsGHD+oEDB+6+6aabWgBkZWXVWLx4ccONGzd++tJLL2267rrrUg8ePAjA2LFj02bMmJG9devWdVu2bKmxcOHCeiW/b2X2f+CBBxonJycf/PLLL9eNGzcud8KECSkAubm58VOmTGm5evXqz9asWfPZlClTWubl5cUDTJgwIWXcuHG52dnZ65KTkw8+8MADjQEWLFiQvGXLlhpbt25dN3PmzOyxY8emHulxDcvSktjRoYMboNGzp3uP67rrIPBvSCRmNWzY8MfTpH379sUZ407MFi5cWP/cc8/9tmbNmjYjIyM/LS3thzfffLN2dnZ24t69e+P69u27Ly4ujuHDh3+zZMmSBiW/b2X2X758ef1Ro0Z9AzBy5Mj/vPfee3ULCwtZsmRJcs+ePb9v1qzZoSZNmhzq2bPn94sXL04uLCxk1apVdUeOHPkfgFGjRn2zbNmy+gBLly6tP3z48G/i4uI488wz933//fcJ2dnZiSVzlkdD3sW7hg3dopI33AD33w/r1rlh8Y0aHXZXkQobNYqj163jiC9HladDB/47Zw7bjnS/Ll26HL9v3774ktvvueeebZmZmXsArrnmmlYLFixoVLdu3UMrV67cCJCTk5PUvXv3vUXPb9myZf62bduSkpKSbIsWLX6cdyYtLS1/586dPyuDyuyfm5ub1Lp163yAxMRE6tSpcyg3NzchJycnMSUlJb9on1atWuXn5OQk5ubmJtStW/dQYqJ7+fT09Pzc3NwkgJ07dyamp6f/uE+LFi3ys7OzE9PS0io8Z45KS8JCQoIbBt+pE1x1lZtwd+lSdyYmEm2ysrI2Hu4506dPz5k+fXrO5MmTm997771Np02btsOWcp+IMcaWsf1n2yqz/5HuU9bzy/leP9tWHpWWhJURI9yMGeec425Efuop97lIVVXmjChYKnKmVWTkyJHfDho06Lhp06btSElJyd+2bduPAyF27NiRlJKSUpCenl5Q/MwqOzs7qXnz5j87e6nM/s2bN8//4osvko455piCgoIC9u7dG9+0adNDKSkpBStXrqxbtE9OTk5Sr1699jRv3vzgnj174gsKCkhMTGTr1q1JTZs2LQBo2bJlwdatW398/Z07dyalpqYe0cykek9Lws7JJ8OaNW76p3PPdcPjyxgIJRKRsrKyNhYNyij+UVRYn3zyyVFFz12wYEH9Y445Zj/AsGHDvlu8eHHD/fv3mw0bNiRt3bq1Ru/evfelpaUV1K5du3DFihW1CwsLmTdvXqOzzz77u5KvW5n9Bw0a9N2cOXMaAcydO7fBKaecsicuLo7MzMzdK1eurJeXlxefl5cXv3LlynqZmZm74+Li6N69+565c+c2AJgzZ06jwYMHfwcwdOjQ7+bNm9eosLCQFStW1K5bt+6hI7k0CDrTkjDVsiWsXAljxsBtt8Hate6+rrp1D7urSMT73e9+l7Jly5YaxhibkpKS/9hjj2UDdO3a9UBmZua3bdu2bR8fH8/UqVOzExLc/8ZnzJiRffnll7c+cOCAOeOMM77/9a9/vRtg3rx5yR9++GHt+++/f0dl9h8/fvzXw4YNa52amtohOTn50Pz58z8HaNas2aEbbrhhR5cuXX4BcOONN+5o1qzZIYD77rtv+29+85tj7rjjjlbt27f/7/jx478GOP/883f/4x//SE5LS+tQs2bNwkcffXTrkR4bU9o1xnDRtWtXu2bNGt8xxCNr3arIEydCu3bufa42bXynknBkjMmy1nYtvm3t2rVbO3Xq9LWvTFJ5a9eubdypU6f0ktt1eVDCmjFuGPzLL0NODnTrBitW+E4lIr6otCQi9O3r7udq0QL693dnX2F8kUBEgkSlJRHjmGPcDBpDhrizr1GjNOGuHFZhYWHhkY2pFu8Cf2elDr9SaUlEqVsXFi2CW2910z717q0Jd6Vc6/Ly8pJVXJEjsJ5WMrCutMc1elAiTlycG1HYsaNbYLJrV1i4EHr08J1Mws3Bgwev2LVr16O7du3SysWR48eVi0t7UKMHJaJ98om7+fjLL2H6dDebhsSm0kYPSvTRbx4S0X75SzdA48wz3T1do0fDDz/4TiUiwaLSkojXoAEsXw433QSPPOKWOcnJ8Z1KRIJBpSVRIT4e7rzTvbe1bp17n+vdd32nEpHqptKSqDJsGHzwAdSp40YWzpyp+7lEoolKS6JO+/bufa5+/WDsWLjiCt3PJRItVFoSlerXh2XL4P/+D+bMcSsjb9/uO5WIVJVKS6JWXBz88Y+weDF89hl06QJvveU7lYhUhUpLot4557j3uerXd0PjH3pI73OJRCqVlsSEdu1g9WoYMACuuUbzFopEqrAsLWPMEGPM7N27d/uOIlEkOdmtx1U0b+Hpp8O2sFmAXUQqIixLy1q7zFo7Ojk52XcUiTJF8xYuXQobN7r3ud5803cqEamosCwtkWAbOtRdLmzY0K3VpfW5RCKDSktiVkaGK65Bg9z6XJddBvv3+04lIuVRaUlMq1cPnnsObr8dnnoKTj0VvvjCdyoRKYtKS2JeXBzccou7GXnLFjdv4csv+04lIqVRaYkEDB4Ma9ZASgoMHAh33AGFpS74LSK+qLREijn2WFi1Ci66yE0BdfbZ8N13vlOJSBGVlkgJtWq597emT4eXXnKXCz/+2HcqEQGVlkipjIFx42DlSvjvf6F7d5g3z3cqEVFpiZSjRw/45z+hWze4+GK49lrIz/edSiR2qbREDqN5c3jtNZgwwV0yPOMM2LHDdyqR2KTSEqmAxES47z74+99h7Vo3/dPbb/tOJRJ7VFoiR+A3v3HLnNSr58647r9f0z+JhJJKS+QItW/vpn8aMgSuv94Nj9+713cqkdig0hKphORktyLyPffAs8+60YWbNvlOJRL9VFoilWQM/P73bsqn3Fw3wnDpUt+pRKKbSkukivr2hawsOP54yMyEm26CQ4d8pxKJTiotkWqQmgpvvQWjR8Pdd8OAAfD1175TiUSfsCwtY8wQY8zs3bt3+44iUmE1asCsWfDYY244fJcubsCGiFSfsCwta+0ya+3o5ORk31FEjtioUfDuu+49r9NOgxkzNCxepLqEZWmJRLouXdz0T/36wdVXuymgNCxepOpUWiJB0rAhPP883Hmnm0nj5JPhs898pxKJbCotkSCKi3OjCV95BfLy3LD4+fN9pxKJXCotkRA480z46CM44QS44ALNFi9SWSotkRBp1QreeOOn2eJ79oQvv/SdSiSyqLREQqhotviFC2H9ejjxRDejhohUjEpLxINhw2DNGmjZEgYOhNtu0ywaIhWh0hLxpG1beP99uOQSuP12OOsszaIhcjgqLRGPatWCxx+H2bNh5Up3ufCDD3ynEglfKi0Rz4yBK690s2jEx8Ppp7uBGppFQ+TnVFoiYaJoFo3+/d2QeC0uKfJzKi2RMNKggVuT66673OKS3bq5UYYi4qi0RMJMXBxMngyvvQbffuuK6+mnfacSCQ8qLZEwdcYZbhaNE0+E4cNh7Fg4cMB3KhG/wrK0tJ6WiNOyJbz+OkycCDNnQo8e8PnnvlOJ+BOWpaX1tER+kpgIf/mLe6/riy/cmdeiRb5TifgRlqUlIj83dKi7XJiRAeedB+PHww8/+E4lEloqLZEIkp4Ob7/tCuvBB909XV984TuVSOiotEQiTFIS3H+/u0S4aZO7XLh0qe9UIqGh0hKJUOee625GbtMGMjPdYI2CAt+pRIJLpSUSwdq0gffeg6uvhqlTtUaXRD+VlkiEO+ooeOghmD8fPv0UOneGf/zDdyqR4FBpiUSJ88+HrCw4+mgYPBgmTdLlQok+Ki2RKHLccbBqFYweDVOmQJ8+kJPjO5VI9VFpiUSZmjVh1iyYN8/d13XCCfDyy75TiVQPlZZIlLroIlizBpo3h4ED4Q9/gIMHfacSqRqVlkgUy8hwKyGPHAl33gl9+8LOnb5TiVReWJaWJswVqT61asFjj8Hjj8Pq1e5y4YoVvlOJVE5YlpYmzBWpfpddBh9+CI0awa9+BbfcosuFEnnCsrREJDjat3fFddll8Kc/wZlnanShRBaVlkiMqV0b5s6FJ59093V16qSbkSVyqLREYtQll7jSatXK3Yx8ww2Qn+87lUj5VFoiMez4493owrFj3UKTWupEwp1KSyTG1agBDz8MCxbAhg1u7kKtjCzhSqUlIoBbDflf/3JnX+ed586+DhzwnUrkf6m0RORHrVu7lZEnToSZM6F7d9i40XcqkZ+otETkfyQlufe3li+H7duhSxd46infqUQclZaIlGrQIHe5sEsXuPRSNxXUvn2+U0msU2mJSJlSUtyUT7fcAk88AV27wscf+04lsSwsS0tzD4qEj4QEuP12eO01+O47OPlkmD0brPWdTGJRWJaW5h4UCT99+sDatdCzJ1x1FVxwAej3Sgm1sCwtEQlPTZvCiy/C3Xe7e7lOPNGt2SUSKiotETkicXEwaRK89ZabJb5HD5g2TZcLJTRUWiJSKT16wEcfwVlnwYQJbv7CvDzfqSTaqbREpNIaNoTnnoOHHnKjDDt2dAM2RIJFpSUiVWIMXH21WxW5QQPo189dPiwo8J1MopFKS0SqRceOblDGlVfClClw2mnw+ee+U0m0UWmJSLWpVQtmzXIzxm/a5GaMf/pp36kkmqi0RKTaFc0Y37EjDB8OI0bAnj2+U0k0CMvS0owYIpEvLQ3efNNNAfXUU24Ow6ws36kk0oVlaWlGDJHoUDQF1Ouvw/79cMopMHUqFBb6TiaRKixLS0SiS69ebgqowYPdWl2DBkFuru9UEolUWiISEg0buqmfZs50lw07dYJXXvGdSiKNSktEQsYYGDMGPvwQGjeG/v3hxhshP993MokUKi0RCbkOHdzNyGPGwL33wqmnwubNvlNJJAhZaRljahtjnjDGPGKMGR6q1xWR8FSrlrtUuHixuwm5c2c3ylCkPFUqLWPMHGPMV8aYdSW2DzDGbDTGbDbGTApsPhdYaK29EhhaldcVkehxzjnunq7OneHSS+GSS3RPl5StqmdajwMDim8wxsQDDwMDgXbAhcaYdkAKsC3wtENVfF0RiSKpqW5Y/G23uRk0OnfW5UIpXZVKy1r7FvBtic0nAZuttVustfnA34Gzge244ir3dY0xo40xa4wxa/K0zoFIzEhIgFtvhZUroW1baNXKdyIJR8F4T6sVP51RgSurVsBiYJgxZiawrKydrbWzrbVdrbVdmzRpEoR4IhLOTjsNXngBatb0nUTCUUIQvqcpZZu11u4DRgbh9UREJEYE40xrO3B0sa9TgB1BeB0REYkxwSitD4HjjDGtjTFJwAXA80F4HRERiTFVHfL+DLAKON4Ys90Yc7m19iAwDngZ+Ax41lr7adWjiohIrKvSe1rW2gvL2P4C8EJVvreIiEhJmsZJREQihkpLREQihkpLREQihrHW+s5QJmNMHpBdYnNj4GsPccqiPGULpyygPIcTTnkqkyXNWqsZCaJcWJdWaYwxa6y1XX3nKKI8ZQunLKA8hxNOecIpi4QXXR4UEZGIodISEZGIEYmlNdt3gBKUp2zhlAWU53DCKU84ZZEwEnHvaYmISOyKxDMtERGJUSotERGJGGFVWsaYo40xbxhjPjPGfGqMGV/ssWuMMRsD2/8c2JZujNlvjPlX4OOvochjjJlf7DW3GmP+VWyfycaYzYGs/X3mCebxKSfLCcaY9wOvt8YYc1KxfXwcm1LzePzZ6WSMWWWM+cQYs8wYU6/YPj6OT6l5QnB8ahhjVhtj1gby3B7Y3tAY86ox5t+BPxsU2ydox0ciiLU2bD6AFsCJgc/rApuAdsAZwGvAUYHHmgb+TAfWhTpPiefcB9wS+LwdsBY4CmgNfA7Ee8wTtONTzt/VK8DAwPazgDd9Hpty8nj52cEt3dMrsH0U8CfPx6esPME+PgaoE/g8EfgA6A78GZgU2D4JmBKK46OPyPkIqzMta+1Oa+0/A5/vwS1t0gr4LXCPtfaHwGNfec4DgDHGAOcDzwQ2nQ383Vr7g7X2C2AzcBLVpBJ5gqacLBYoOntI5qcFQH0dm7LyBFU5eY4H3go87VVgWOBzX8enrDxBZZ29gS8TAx8WdxyeCGx/AsgMfB7U4yORI6xKqzhjTDrQGfcbWFvgdGPMB8aYlcaYbsWe2toY81Fg++khylPkdCDXWvvvwNetgG3FHt9OsVLxkAdCcHxKZLkOuNcYsw34CzA58DRfx6asPODnZ2cdMDTw0K/5aZVvX8enrDwQ5ONjjIkPXMr+CnjVWvsB0MxauxNc0QJNA08P2fGR8BaWpWWMqQMsAq6z1n6PW/erAe7ywQ3As4Gzip1AqrW2MzABeLr4ewRBzFPkQv73rMaUsnu131NwBHmCfnxKyfJb4Hpr7dHA9cBjRU8tZfdQHJuy8vj62RkFXG2MycJdpssvemopu4fi+JSVJ+jHx1p7yFp7ApACnGSM6VBe9NK+RXXmkcgQdqVljEnE/aOaZ61dHNi8HVgcuKSwGigEGgcuFXwDYK3Nwl3nbhuCPBhjEoBzgfnFnr6d//1NNYVqvhx1JHmCfXzKyHIZUPT5An66hOPr2JSax9fPjrV2g7W2n7W2C+4XjM8DT/dyfMrKE4rjU8Ra+x3wJjAAyDXGtAjkbYE7C4MQHB+JDGFVWoGzp8eAz6y1U4s9tAToE3hOWyAJ+NoY08QYEx/Y3gY4DtgSgjwAfYEN1trtxbY9D1xgjDnKGNM6kGe1rzzBPD7lZNkB9Ap83gcoulTp69iUmsfXz44xpmngzzjgD0DRqDwvx6esPCE4Pk2MMfUDn9ck8POLOw6XBZ52GbA08HlQj49EEN8jQYp/AKfhTvk/Bv4V+DgLV1J/w11//yfQJ/D8YcCnuFFF/wSGhCJP4LHHgTGl7HMz7rfSjQRGrfnKE8zjU87f1WlAVuA1PwC6+Dw2ZeXx9bMDjMeN3NsE3ENgVhqPx6fUPCE4Ph2BjwJ51vHTiNdGwArcLxcrgIahOD76iJwPTeMkIiIRI6wuD4qIiJRHpSUiIhFDpSUiIhFDpSUiIhFDpSUiIhFDpSUiIhFDpSUiIhHj/wEyuD7FHYVsogAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAD4CAYAAABfYrnHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhyUlEQVR4nO3deXhV1bnH8e+biYShUSEQICRYC2VGJeKAAlW0UGRoqShilaHFCUURrUPFWr2t1IpQVKoXUVQUAbkiOHGlglBxICqKENRCIrNRCgKCQLLuHytRLgYIJCd7n3N+n+fJwzk7e+e8G56HX/be71rLnHOIiIhEg4SgCxAREakohZaIiEQNhZaIiEQNhZaIiEQNhZaIiESNpKALOJR69eq5pk2bBl2GiESBvLy8L51zGUHXIZEVytAys15Ar5/85CcsXbo06HJEJAqYWWHQNUjkhfL2oHNujnNuWHp6etCliIhIiIQytERERMqj0BIRkaih0BIRkagRytAys15m9si2bduCLkVEREKk2kLLzH5sZo+a2czD7atGDBERKU+FQsvMJpvZF2a2/IDt3c1slZl9ZmY3H+pnOOdWO+eGVqZYERGJbxW90noc6L7/BjNLBB4EegCtgAFm1srM2prZ3AO+6ldp1eUYPx4WLYr0p4iISJAqNLjYOfeGmTU9YHNH4DPn3GoAM5sG9HHO/QU4/2gLMrNhwDCA7OzsCh3zzTcwcSJs3Aivvw4nn3y0ny4iImFWmWdajYG1+71fV7qtXGZW18z+AZxkZrccbD/n3CPOuVznXG5GRsVmZKlZE157DY49Fn7+c1i5soJnICIiUaUyoWXlbDvoMsjOua+cc1c4504ovRqrUllZPriSkqBbN1izpqo/QUREglaZ0FoHNNnvfRawoXLlVM5PfgLz5sGuXT64NgRajYiIVLXKhNa7QDMzO97MUoCLgBeqpqyj17YtvPIKfPEFnHsufPll0BWJiEhVqWjL+zPAEuCnZrbOzIY65/YBw4FXgZXAdOfcx5ErteI6doQXXoDVq/0zrq1bg65IRESqgjl30MdQgcvNzXWVWZrk5ZehTx/o0MHfNqxTpwqLE5FQMbM851xu0HVIZIVyGqeq0qMHPPssvPsu9OrlW+NFRCR6xXRoAfzyl/Dkk/DGG/71t98GXZGIiBytmA8tgAEDYNIkf4uwf3/YuzfoikRE5GjERWgBDBkCDzzgGzQuuQSKi4OuSEREjlSFpnGKFVdfDbt3w6hRkJoKjz0GCXET2yIi0S+uQgvghht8Q8bo0ZCW5ucstPLm9hARkdCJu9AC+MMffHDdc48PrrFjFVwiItEgLkPLDP78Zx9c48ZBrVpw991BVyUiIocTl6EFPrjGjfPzFP7Xf/krrttuC7oqERE5lLgNLfDBNXGiD64//MEH18iRQVclIiIHE9ehBZCY6LsId+/2TRrJyXDNNUFXJSIi5Yn70AK/BtfTT/uxW9de64PsqquCrkpERA6kUUqlkpNh2jTo3duP53r44aArEhGRAym09pOSAtOnQ8+ecMUVfuonEREJD4XWAWrUgOee8zPEDxvmn3eJiEg4KLTKUaMGzJrlVz4eOhSeeCLoikREBBRaB5WaCs8/D+ecA4MGwdSpQVckIiIKrUNIS4PZs6FrV7j0UnjmmaArEhGJbwqtw6hZE+bMgbPO8kuaTJ8edEUiIvFLoVUBtWrB3LnQqRNcfDE8+2zQFYmIxCeFVgXVrg0vvghnnOGD6+mng65IRCT+hDK0zKyXmT2ybdu2oEv5f+rUgZdfhs6d4Te/gaeeCroiEZH4EsrQcs7Ncc4NS09PD7qUH6hVy19xlTVnTJkSdEUiIvEjlKEVdmXNGd26weDBMHly0BWJiMQHhdZRqlnTt8Ofd54fgPzII0FXJCIS+xRalZCW5gcg/+IXcPnlfm0uERGJHIVWJaWm+imfevXyy5k88EDQFYmIxC6FVhWoUQNmzoQ+ffwCkuPGBV2RiEhsUmhVkZQUmDEDfvUruP56uO++oCsSEYk9Cq0qVLaQ5AUXwKhR8Ne/Bl2RiEhsSQq6gFiTnOxny0hMhN//Hvbtg1tvDboqEZHYoNCKgKQkePJJH1y33Qa7d8Odd4JZ0JWJiEQ3hVaEJCX52TJSU+Guu2DnTvjb3xRcIiKVodCKoMREP+i4Zk0YOxa++QYefBAS9CRRROSoKLQiLCEBxo/3wTVmDOzaBZMm+SsxERE5MvqvsxqYwV/+4ifbHT3aB9dTT/mmDRERqTiFVjUxg9tv91dco0b54Jo+3T/zEhGRitHTlWp2ww3+udacOdC7t3/OJSIiFRPK0ArrIpBV5aqr4LHHYP586N4dvv466IpERKJDKEMrzItAVpVBg/wg5CVL4NxzYcuWoCsSEQm/UIZWvLjwQnjuOfjgAzj7bPjii6ArEhEJN4VWwHr39s+3PvkEunSB9euDrkhEJLwUWiFw3nnwyis+sM48Ez77LOiKRETCSaEVEp07wz//Cdu3++D68MOgKxIRCR+FVojk5sKiRX62jC5d4M03g65IRCRcFFoh07Il/OtfUK+e7yqcNy/oikREwkOhFUI5ObB4MTRrBuefDzNnBl2RiEg4KLRCqkEDWLAAOnb0rfGPPhp0RSIiwVNohdgxx8Crr/ruwt/+1q/HJSISz0IZWrE+jdORqFULZs+G/v3hxhvh1lvBuaCrEhEJRihDKx6mcToSKSl+yqdhw/wSJ1dfDSUlQVclIlL9tDRJlEhMhH/8A4491i8muXUrTJmiNblEJL4otKKIGdxzjw+um2+Gbdtgxgy/RpeISDwI5e1BObTf/x4efhhefhm6dYOvvgq6IhGR6qHQilLDhvmrrLw8OOssWLs26IpERCJPoRXF+vXzM2asXw9nnAEffxx0RSIikaXQinJdusAbb8C+ff6KS/MVikgsU2jFgPbtfVjVqwfnnOPX5xIRiUUKrRhx/PF+ot02beCXv4TJk4OuSESk6oUytDQjxtHJyIDXX/dXW0OH+oHImj1DRGJJKENLM2Icvdq1/e3Biy/2Uz5dd51mzxCR2KHBxTEoJQWefBLq14dx42DzZj97Ro0aQVcmIlI5Cq0YlZAAY8dCw4Z+MPJXX8GsWVCnTtCViYgcvVDeHpSqYQY33QSPP+6fdXXpAhs3Bl2ViMjRU2jFgcsu88+5PvkETj8dVq4MuiIRkaOj0IoTPXrAwoWwezd06gSLFgVdkYjIkQtlaKnlPTI6dIAlS3yDRrduMH160BWJiByZUIaWWt4jp2wQ8imnwIUX+mYNjeUSkWgRytCSyKpbF157zU+4e8MNfixXcXHQVYmIHJ5CK06lpvrbg9dfD3//O/TvD7t2BV2ViMihKbTiWNlYrrFj4X/+xz/n+vLLoKsSETk4hZZw/fX+qisvz3cWrl4ddEUiIuULZWipe7D6/frX/jnXl1/6sVzvvht0RSIiPxTK0FL3YDDOPNN3FtasCV27wty5QVckIvL/hTK0JDgtWvixXC1bQp8+MGFC0BWJiHxPoSU/kJnpZ8/o1QuuvdZ/qSVeRMJAoSXlqlULnnsORo70V1t9+8KOHUFXJSLxTqElB5WYCPfdBw89BC+9BGedBevXB12ViMSzUIaWugfD5corfVPGZ5/BqafCBx8EXZGIxKtQhpa6B8OnRw/fWZiQ4LsM1VkoIkEIZWhJOLVrB2+/7TsM+/Tx0z+JiFQnhZYckYYNv+8sHDFCnYUiUr1CGVp6phVu6iwUkaCYq6bFlMysFvAQsAdY4JyberhjcnNz3dKlSyNemxy9iRPhmmugbVuYMweysoKuSOKVmeU553L335aXl1c/KSlpEtCGkP6SLj9QAizft2/fbzt06PDFgd9MqsxPNrPJwPnAF865Nvtt7w6MBxKBSc65e4BfATOdc3PM7FngsKEl4XfllX5hyf79oWNHeP55/6dIGCQlJU3KzMxsmZGR8Z+EhAQtdxoFSkpKrKioqNWmTZsmAb0P/H5lf/N4HOi+/wYzSwQeBHoArYABZtYKyALWlu6mpyAxpHt3P/VTaip06QLTpgVdkch32mRkZHytwIoeCQkJLiMjYxv+6viH36/MD3fOvQFsOWBzR+Az59xq59weYBrQB1iHD65Kf66ET+vWvrPwlFNgwAAYPRpKSoKuSoQEBVb0Kf03KzcnIhEejfn+igp8WDUGZgH9zGwiMOdgB5vZMDNbamZLi4qKIlCeREpGhl/eZMgQuOsuf8tw586gqxKRWBKJ0LJytjnn3E7n3GDn3JWHasJwzj3inMt1zuVmZGREoDyJpJQUmDTJT/80axZ07gzr1gVdlUj127RpU2KLFi1atWjRolW9evXa169fv13Z+927d5f3/+R3RowY0ah58+atWrRo0apTp07NCgoKkgFWrVqVkpqaenLZz7n44ouzy45ZtGhRzebNm7fKzs5uM2jQoCYlB7nVccstt2RmZ2e3adq0aZvnnnvuR4c7fteuXdazZ88fZ2dnt2nXrl2LVatWpZQdM2HChLo5OTltcnJy2kyYMKFu2fb8/PyUdu3atcjJyWnTs2fPH5edb0lJCYMGDWqSnZ3dpnnz5q0WL15c80j/XiMRWuuAJvu9zwI2ROBzJKTMfDv8nDnw6af+luE77wRdlUj1yszMLM7Pz1+Rn5+/4tJLLy264oorNpe9T01NPeQtyzvuuGPTJ598siI/P39Fjx49tt16660Ny77XpEmTb8t+ztNPP/152farrroq56GHHiosKChYvnr16tSZM2f+6MCfm5eXlzpr1qzjVq1a9fErr7zyyXXXXZe9b9++Qx4/fvz4eunp6fs+//zz5cOHD988cuTILIDNmzcnjhkzptE777yzcunSpSvHjBnTqKioKBFg5MiRWcOHD99cWFi4PD09fd/48ePrAcyYMSN99erVqQUFBcsnTpxYeNVVV2UfWOPhRCK03gWamdnxZpYCXAS8EIHPkZDr2dM3aKSlqUFD5Egcd9xx310m7dy5M8HskBdmFBYWJu/YsSOhW7duOxMSEhg4cOBXzz///LEH7jdz5sxjfvWrX21JS0tzLVq02JOTk/PtggULah3q+Llz5x4zZMiQrwAGDx78nzfffLNOSUkJzz//fHrnzp2/btCgQXFGRkZx586dv541a1Z6SUkJS5YsqTN48OD/AAwZMuSrOXPmHAMwe/bsYwYOHPhVQkIC55xzzs6vv/46qbCwMPlI/m4q2/L+DNAVqGdm64A7nHOPmtlw4FV8y/tk59zHlfkciV5lDRr9+vkGjRUr4I9/9HMYilSnIUNosnw5R3w76lDatOGbyZP/3zP8CunQocNPd+7cmXjg9nvuuWdt3759twNcc801jWfMmFG3Tp06xQsXLlxVts+6detSWrZs2ap27drFd9111/ru3bvvKCwsTG7YsOHesn1ycnL2bNy48QdhsH79+pTTTjvtu6kAGjVqtGft2rUpKSkp7mDHb968OeX444/fA5CcnEzt2rWLN2/enLR+/frkrKysPWXHNG7ceM/69euTN2/enFSnTp3i5GT/8U2bNt2zefPmFICNGzcmN23a9LtjGjZsuKewsDA5Jyfnu88+nEqFlnNuwEG2vwS8VJmfLbGjrEHjyit9g8aKFTBlip9ZQyQe5eXlrTrcPhMmTFg/YcKE9bfcckvmvffeW//+++/fkJ2dvXfNmjUfZmZmFi9atKjmBRdc8JMVK1YsL2+SiPKuzg6ynzvU8Ud6zMH2P8TP+sG2Q6lUaIlUVFmDRuvWMGoUrFkDs2drBg2pPkdzRRQpFbnSKjN48OAtPXv2bHb//fdvSEtLc2lpacUAZ5111jfZ2dnfLl++PLVp06Z797+yKiwsTMnMzPzB1UtWVtaetWvXftdIsWHDhpSsrKy9hzo+MzNzz5o1a1JOOOGEvXv37mXHjh2J9evXL87Kytq7cOHCOmXHrF+/PqVLly7bMzMz923fvj1x7969JCcnU1BQkFK/fv29AI0aNdpbUFDw3edv3LgxJTs7u8JXWaDxUlKNymvQePPNoKsSqX55eXmrypop9v8qC6yPPvqoRtm+M2bMOOaEE07YBbBhw4akssaJFStWpBQUFNT46U9/+m1OTs7eWrVqlcyfP79WSUkJU6dOrdunT5+tB35uv379ts6aNeu4Xbt2WX5+fkpBQUFq165ddx7q+J49e26dPHlyXYDHHnvs2NNPP317QkICffv23bZw4cIfFRUVJRYVFSUuXLjwR3379t2WkJDAaaedtv2xxx47FmDy5Ml1zz///K0AvXv33jp16tS6JSUlzJ8/v1adOnWKj+TWIOhKSwJQ1qDRpw907epXRv7tb4OuSiQ8Ro0albV69epUM3NZWVl7Hn300UKAefPm1b777rsbJyYmusTERDdu3LjCBg0aFAM89NBDhUOHDj1+9+7d9rOf/ezrCy64YBvA1KlT0999991a48aN25Cbm7u7b9++W5o3b946MTGRsWPHFiYl+Rg42PEjRoz4sl+/fsdnZ2e3SU9PL3722Wf/DdCgQYPiG2+8cUOHDh1aAtx0000bymq577771l144YUn3H333Y1bt279zYgRI74E6N+//7YXX3wxPScnp01aWlrJpEmTCo7076baJsw9GpowN7Zt2eKbM+bNg6uvhvvvh+Qj6iMS+V55E+YuW7asoH379l8GVZMcvWXLltVr37590wO36/agBOa44+DFF/0zrgcfhHPPBU2CIiKHotCSQCUlwb33wlNP+db43Fx4//2gqxKRsFJoSSgMHAiLF/tJdjt10kBkqTIlJSUlR9ZTLYEr/Tcrdx4qhZaERocOsHSp/3PAALj5ZijWIjZSOcuLiorSFVzRo3Q9rXRgeXnfV/eghEqDBjB/PowYAWPGwIcfwtNPwzHHBF2ZRKN9+/b9dtOmTZM2bdqklYujx3crF5f3TYWWhE5KCkycCCeeCMOH+5WQZ8+Gli2DrkyiTely7T9Y/Vail37zkNC6/HJ4/XXYtg1OPdUPShaR+KbQklA780z/nKtZMz8Y+U9/0orIIvFMoSWh16SJ7ywcOBDuuAP69oWtW4OuSkSCoNCSqJCWBk88ARMmwMsv+3kLl5fbWyQisUyhJVHDzDdmLFgAO3b451wazyUSXxRaEnU6dYL33oOTTvLjuUaOhL1HNE+0iEQrhZZEpYYN4Z//hGuu8RPtnnsubN4cdFUiEmkKLYlaKSnw97/Dk0/CO+/4mTTeeivoqkQkkhRaEvUuucSvz5WSAp07w8MPQ4hX3BGRSlBoSUxo396P5zrnHLjiChg6FHbtCroqEalqCi2JGccdB3Pnwu23w2OP+YHJhYVBVyUiVUmhJTElMdHPmvHCC/Dvf/vnXK++GnRVIlJVQhlaZtbLzB7Ztm1b0KVIlOrVC959Fxo1gh494I9/1DInIrEglKHlnJvjnBuWnp4edCkSxZo1892El14Kd94J3btDUVHQVYlIZYQytESqSs2a/vnWpEmwaJEfkPyvfwVdlYgcLYWWxDwz30341luQmgpdu8LYsWqLF4lGCi2JGyeeCHl5/nnXDTdAv35+rS4RiR4KLYkr6enw3HNw332+w7BDB/jgg6CrEpGKUmhJ3DHzk+wuWOAHIJ9+OkyeHHRVIlIRCi2JW2eeCe+/72eNHzoUBg+Gb74JuioRORSFlsS1+vX94OPbb4cpU+C00+CTT4KuSkQORqElca9sFo2XXoL16yE3F559NuiqRKQ8Ci2RUt27+9uFrVvDRRfB5Zdr0l2RsFFoiewnOxveeANuugkeeQROPRXy84OuSkTKKLREDpCcDGPG+NuFGzf624VPPhl0VSICCi2Rg+rRw4/h6tDBz184eDDs3Bl0VSLxTaElcgiNG8P8+d93F55yCixfHnRVIvFLoSVyGElJvrtw3jzYssUH16RJmrtQJAihDC2tpyVh1K2bv13YqRP87ndwySWwfXvQVYnEl1CGltbTkrDKzPSDke++G6ZN09yFItUtlKElEmaJiXDbbfD6637ap9NOg4ce0u1Ckeqg0BI5Sp07+6uss8+Gq6/2S51s2RJ0VSKxTaElUgn16sHcufC3v/k/27eHhQuDrkokdim0RCopIcEvKrlkCaSlwc9+5lvk9+0LujKR2KPQEqkiHTrAe+/BoEG+UaNzZygoCLoqkdii0BKpQrVr+wUlp02Djz/2twunTQu6KpHYodASiYALL4Rly6BNGxgwwE8BtWNH0FWJRD+FlkiENG3qmzJGj4YnnoCTT4alS4OuSiS6KbREIigpCe6804/p2r0bzjgD7r0XSkqCrkwkOim0RKpB587+dmHv3n6tru7d/bInInJkQhlamntQYtGxx8KMGX5xycWLfZPGiy8GXZVIdAllaGnuQYlVZn6y3ffeg0aN4Pzz/Wwa33wTdGUi0SGUoSUS61q0gLffhlGjYOJENWmIVJRCSyQgNWr4pozXXvMrIp9+Ovz5z1BcHHRlIuGl0BIJ2Nlnw4cfwq9/7WeP79IF1qwJuiqRcFJoiYTAscfCM8/A1Knw0Ue+SWPKFC13InIghZZIiFx8sb/qOvlkP4dh//7w1VdBVyUSHgotkZDJyYH582HMGJg9G9q2hXnzgq5KJBwUWiIhlJjoByG//ba/dfjzn8OIEbBrV9CViQQrlKGlwcUi3kkn+Vb4a6+Fv/8dcnP9aski8SqUoaXBxSLfS0uD8ePh1VfhP/+Bjh3hnnvUGi/xKZShJSI/dN55vrOwTx+45RY46yz49NOgqxKpXgotkShSty5Mn+5b41eu9K3xDzygWeMlfii0RKKMmW+NX77cD0S+5hp/Ffb550FXJhJ5Ci2RKNW4Mbz0Ejz8MLz1lm+Nf/xxDUiW2KbQEoliZjBsmB+QfOKJMHgw9O0LmzcHXZlIZIQytNTyLnJkfvxjvzryfff5LsPWrWHmzKCrEql6oQwttbyLHLmEBBg50q/V1bQpXHABDBwIW7YEXZlI1QllaInI0WvVCpYsgTvv9J2GbdvCyy8HXZVI1VBoicSg5GQYPfr7aaB+8Qu4/HLYvj3oykQqR6ElEsPKVkS+8Ub47/+GNm38opMi0UqhJRLjUlPhr3+FxYv9lFDnnus7Dr/+OujKRI6cQkskTpxxBrz/vr/qevRRf9X16qtBVyVyZEIZWmp5F4mMtDR/1fXmm1C7NnTvDkOHwtatQVcmUjGhDC21vItE1qmn+tb4m2/2s2i0aeNn1xAJu1CGlohEXmoq/OUvfgqoY46Bnj1h0CC//IlIWCm0ROLcKadAXh7cdhs89ZSfTWPOnKCrEimfQktEqFED7r7bj+uqVw9694bf/EazaUj4hDK01IghEowOHfy4rtGjYdo0P7vG888HXZXI90IZWmrEEAlOSoqfAuqddyAzE375S7jwQs0cL+EQytASkeCddJIPrj/9yV9ttWoFTzyh9bokWAotETmolBS4/Xb44ANo0QIuu8yP7SooCLoyiVcKLRE5rJYtYdEimDDBD0xu3RrGjYPi4qArk3gTytBSI4ZI+CQkwPDh8PHH0KULXH89dOoEy5cHXZnEk2oLLTP7sZk9amaHXU9VjRgi4ZWdDS++CFOnwr//7WeSv+MO+PbboCuTeFCp0DKzyWb2hZktP2B7dzNbZWafmdnNAM651c65oZX5PBEJBzO4+GJYsQL69/fNGief7BefFImkyl5pPQ5033+DmSUCDwI9gFbAADNrVcnPEZEQysjws2i89JJfYLJTJ7j2WtixI+jKJFZVKrScc28AB46Z7wh8VnpltQeYBvSp6M80s2FmttTMlhYVFVWmPBGpJj16+GddV18NDzzgGzVeeSXoqiQWReKZVmNg7X7v1wGNzayumf0DOMnMbjnYwc65R5xzuc653IyMjAiUJyKRUKeO7y5cvBhq1fJBdskl8MUXQVcmsSQSoWXlbHPOua+cc1c4505wzv0lAp8rIiFQttjk6NEwfbof3zVpEpSUBF2ZxIJIhNY6oMl+77OADRH4HBEJqRo1/FRQy5ZB27bwu99B166+cUOkMiIRWu8CzczseDNLAS4CXojA54hIyLVsCQsWwOTJ/pnXiSf6GTZ27w66MolWlW15fwZYAvzUzNaZ2VDn3D5gOPAqsBKY7pz7uPKlikg0MoPBgyE/Hy66yC+B0rYtzJ8fdGUSjSrbPTjAOdfQOZfsnMtyzj1auv0l51zz0udX/1U1pYpINMvI8BPuvvaaf9+tG1x6KahJWI5EKKdxEpHYdc458NFH/jbhtGm+UWPyZM0eLxWj0BKRapea6mfR+OADv+TJ0KG+USM/P+jKJOwUWiISmFatYOFC3xL/0UfQrp2fx1CNGnIwCi0RCVRCgr/Sys//fh7Ddu3g00+DrkzCSKElIqFQv76fx3DePGjWDJo0OfwxEn+Sgi5ARGR/557rv0TKoystERGJGgotERGJGgotERGJGgotERGJGgotERGJGgotERGJGgotERGJGgotERGJGuZCPLWymRUBhUHXUcXqAV8GXUSExcM5gs4zbHKccxlBFyGRFerQikVmttQ5lxt0HZEUD+cIOk+RIOj2oIiIRA2FloiIRA2FVvV7JOgCqkE8nCPoPEWqnZ5piYhI1NCVloiIRA2FloiIRA2FVhUysyZm9rqZrTSzj81sROn2E83sLTP7wMyWmlnH/Y65xcw+M7NVZvbz4KqvODNLNbN3zGxZ6XneWbr9ODP7XzP7tPTPY/c7JpbO814zyzezD83sf8zsmP2OiZnz3O/7o8zMmVm9/bZF3XlKjHDO6auKvoCGwMmlr+sAnwCtgHlAj9LtvwAWlL5uBSwDagDHA/8GEoM+jwqcpwG1S18nA28DpwF/BW4u3X4zMCZGz/M8IKl0+5hYPc/S902AV/GD/OtF83nqKza+dKVVhZxzG51z75W+3g6sBBoDDvhR6W7pwIbS132Aac65b51za4DPgI6EnPN2lL5NLv1y+POZUrp9CtC39HVMnadzbp5zbl/p9reArNLXMXWepe/vB27a7z1E6XlKbFBoRYiZNQVOwv/Weh1wr5mtBf4G3FK6W2Ng7X6HrSvdFnpmlmhmHwBfAP/rnHsbaOCc2wg+wIH6pbvH2nnubwjwcunrmDpPM+sNrHfOLTtg96g9T4l+Cq0IMLPawHPAdc65r4Ergeudc02A64FHy3Yt5/CoGIPgnCt2zp2Iv8roaGZtDrF7TJ6nmd0G7AOmlm0q70dEvMgqUM55tgNuA0aXs3vUnqdEP4VWFTOzZHxgTXXOzSrdfBlQ9noG399KWYd/ZlAmi+9vHUYF59xWYAHQHdhsZg0BSv/8onS3WDtPzOwy4HxgoHOu7D/sWDrPPvjnVcvMrAB/Lu+ZWSYxcJ4SvRRaVcjMDH8VtdI5N3a/b20AupS+Phv4tPT1C8BFZlbDzI4HmgHvVFe9R8vMMso65swsDegG5OPP57LS3S4DZpe+jqnzNLPuwO+B3s65b/Y7JJbO833nXH3nXFPnXFN8UJ3snNtElJ6nxIakoAuIMZ2A3wAflT4fALgV+B0w3sySgN3AMADn3MdmNh1Ygb/NdLVzrrjaqz5yDYEpZpaI/8VnunNurpktAaab2VDgc+ACiMnz/AzfOfe//vcU3nLOXRFr53mwnaP4PCUGaBonERGJGro9KCIiUUOhJSIiUUOhJSIiUUOhJSIiUUOhJSIiUUOhJSIiUUOhJSIiUeP/AJhSiW3GUfwhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAD4CAYAAABfYrnHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjgElEQVR4nO3deXhV1b3G8e8vAQQUokKYSYIokwgKkaIoongxCgItakVbFakUFUWt9YrV3tZaK6JWUaG1FIdKqYg4AM6gqFeKEC2KAkqRCAFDEBkLMmTdP1aQ3JCEkORk7ZPzfp4nT3P2mV7Dad7svddey5xziIiIxIOk0AFERETKS6UlIiJxQ6UlIiJxQ6UlIiJxQ6UlIiJxo1boAGVp3Lixy8jICB1DROJAdnb2BudcaugcEluRLq2MjAwWLVoUOoaIxAEzywmdQWJPhwdFRCRuqLRERCRuRLK0zOx8M3ts8+bNoaOIiEiERLK0nHMznXMjUlJSQkcREZEIiWRpiYiIlKTaSsvMjjGzv5rZ9HI8VocHRUTkAOUqLTObbGbrzWxJse1ZZrbczFaY2a1lvYZzbqVzbnh53k+HB0VEpCTl3dN6AsgqusHMkoFHgXOBTsBQM+tkZieY2axiX02qNHUxzsH48fDOO7F8FxERCa1cFxc7594xs4xim3sAK5xzKwHM7B/AIOfcH4ABFQ1kZiOAEQBpaWnles6OHTBxIqxfDwsWwLHHVvTdRUQkyipzTqslsLrI7TWF20pkZo3M7E/ASWY2prTHOecec85lOucyU1PLNyNL/fowaxaYwYAB8O235fwvEBGRuFKZ0rIStpW6DLJz7hvn3EjnXNvCvbEq1bYtPP88rFwJF1wAu3dX9TuIiEholSmtNUDrIrdbAWsrF6dyTj8dJk2CuXNh1Ch/rktERGqOykyYuxA4zszaALnAxcAlVZKqEi67DJYvh7vvhvbt4aabQicSEZGqUt4h71OB+UB7M1tjZsOdc3uAUcBrwFJgmnPu09hFLb/f/c4fIrz5Zpg5M3QaERGpKuYifAwtMzPTVXRpkv/8B844A5YuhffegxNPrNpsIhItZpbtnMsMnUNiq8ZO41S/Prz0Ehx1FJx/PqxbFzqRiIhUVo0tLYDmzf3hwW+/hYED/d6XiIjErxpdWuAPC06dCtnZcPnlUFAQOpGIiFRUjS8t8IcH77sPpk+HO+4InUZERCqqMkPe48qNN8KyZfuHwl92WehEIiJyqBJiTwv8FE+PPgp9+8LPfgbvvhs6kYiIHKqEKS2A2rXh2WfhmGPghz+Ef/87dCIRETkUCVVa4IfAz5rlp3gaMAA2bQqdSEREyivhSgv80iXPP+/3tC68UJPriojEi4QsLYDeveEvf4E334TrrtPkuiIi8SBhRg+W5PLL/eS6f/gDdOgAN9wQOpGIiJQloUsL4K674PPP/WzwbdrAoEGhE4mISGkS9vDgPklJ8NRTcPLJMHQozJ8fOpGIiJQmkqVlZueb2WObN2+ulverX9+PKGzZ0s+esXx5tbytiIgcokiWlnNupnNuREpKSrW9Z2oqvPoqJCdDVpZmhRcRiaJIllYobdvC7NmQnw/nnQdbtoROJCIiRam0isnM9BPrfvIJDBkCu3aFTiQiIvuotEqQlQWTJvlruIYP1zVcIiJRkfBD3ktzxRWQmwu33w6tWvlruUREJCyVVhluuw3WrIF77vEjC0eNCp1IRCSxqbTKYAaPPOJHEl5/PTRv7s9ziYhIGDqndRDJyfD3v0PPnnDppVqHS0QkJJVWOdSvDzNnQkYGDBwIn34aOpGISGJSaZVTo0b+4uO6daFfP1i1KnQiEZHEo9I6BBkZ8Npr8J//+OJavz50IhGRxBLJ0qruuQcPRZcufp7CNWv89VyaNUNEpPpEsrRCzD14KHr12j9rxsCBsHNn6EQiIokhkqUVD847D558EubNg4svhj17QicSEan5VFqVcMklMH48vPgiXHWVpnsSEYk1XVxcSdddB998A7/9rR9hOG6cvyhZRESqnkqrCvzP/8CGDXD//dC4Mdx6a+hEIiI1k0qrCpj5w4QbN8KYMdCwIVxzTehUIiI1j0qriiQl+YEZ27bBtddCgwbw05+GTiUiUrNoIEYVql0bpk2Ds86CYcPg+edDJxIRqVkiWVpRvrj4YOrW9aMJTz7ZD4V//fXQiUREao5IllbULy4+mCOOgJdfho4dYfBgeO+90IlERGqGSJZWTXDUUX4vKy0N+veH7OzQiURE4p9KK4aaNIE33vAF1q+fn/ZJREQqTqUVY61bw5w5UK8e9O0Ln30WOpGISPxSaVWDtm19cSUn++L6/PPQiURE4lMkSyueRw+Wpn17X1x79/oh8f/+d+hEIiLxJ5KlFe+jB0vTqRO8+Sbs2OGLKycndCIRkfgSydKqybp08YMztmzxxbVmTehEIiLxQ6UVQLdu8NprkJ/vi2vdutCJRETig0orkB494NVXYe1aPzhj/frQiUREoi+SpVUTB2KU5NRT/cwZq1bB2Wf7dblERKR0kSytmjoQoyS9e8PMmfDFF/Bf/+WXNxERkZJFsrQSTd++fkb4Tz/132uPS0SkZCqtiMjK8rPDL13qi2vDhtCJRESiR6UVIVlZ8NJLsHy5H1WYnx86kYhItESytBJlIEZJ+vXz57hWrPDFpVGFIiL7RbK0EmkgRknOPhtmzfJTPZ15JuTlhU4kIhINkSwt8XtZr7zih8P36aMLkEVEQKUVaWec4Ytr9WpfXGvXhk4kIhJWJEsrkc9pFde79/6ZM/r0gdzc0IlERMKJZGkl+jmt4k47zc9V+PXXfu9r9erQiUREwohkacmBTj0VXn/dD4Pv0we++ip0IhGR6ldtpWVmh5vZk2b2FzO79CCP1eHBEvTs6Zc1+eYbOP10PyxeRCSRVKq0zGyyma03syXFtmeZ2XIzW2FmtxZu/hEw3Tl3FTCwrNfV4cHS9egBc+fC9u3+fNdnn4VOJCJSfSq7p/UEkFV0g5klA48C5wKdgKFm1gloBew7G7O3ku+b0Lp1g7ffBuf8Oa6PPgqdSESkelSqtJxz7wDF5yXvAaxwzq10zu0C/gEMAtbgi6vM9zWzEWa2yMwW5Wseo1J17gzvvAP16vkLkOfPD51IRCT2YnFOqyX796jAl1VLYAYwxMwmAjNLe7Jz7jHnXKZzLjM1NTUG8WqO446Dd9+Fxo39siZvvx06kYhIbMWitKyEbc45t905N8w5d7VzbkoM3jchpaf74kpPh3PP9dd0iYjUVLEorTVA6yK3WwGayyGGmjf3e1kdO8LAgX5tLhGRmigWpbUQOM7M2phZHeBi4KUYvI8UkZrqRxVmZsKFF8Lf/hY6kYhI1avskPepwHygvZmtMbPhzrk9wCjgNWApMM0592nlo8rBHHmkvwD5jDPgsstg/PjQiUREqlatyjzZOTe0lO0vAy9X5rWlYo44AmbPhqFDYfRofyHyb34DVtKZRhGROKNpnGqgunXh2Wfhyivhzjvh+uuhoCB0KhGRyqvUnpZEV61aMGkSHH003Hef3+N68kmoXTt0MhGRilNp1WBmcO+90KgRjBkDmzf7PbD69UMnExGpGB0erOHM4NZb4c9/9gtKnnMObNoUOpWISMWotBLEiBHwzDOwYIFf2uTrr0MnEhE5dCqtBHLhhTBrFnzxhV9Y8ssvQycSETk0Kq0E068fzJkDGzdCr16wZMnBnyMiEhUqrQTUs6efr9DMLyb5zjuhE4mIlI9KK0Edfzy8/z40a+b3vqZPD51IROTgVFoJLD0d3nsPuneHiy6CRx4JnUhEpGwqrQTXqBG8+aafHf666/z1XM6FTiUiUjKVllCvHjz3HIwcCffcA5dfDrt2hU4lInIgzYghACQnw4QJ0LIl3HEH5OX581wNGoROJiKyXyT3tMzsfDN7bPPmzaGjJBQzuP12mDzZD4vXRcgiEjWRLC3n3Ezn3IiUlJTQURLSsGHw0kuwbBmceqq/GFlEJAoiWVoS3nnnwVtvwdatcMop8L//GzqRiIjOaUkZevSA+fN9gfXtC0895YfGi8SL7OzsJrVq1ZoEdEZ/pMeLAmDJnj17fta9e/f1xe9UaUmZjj3WF9fgwfDjH/v5Cm+5RSshS3yoVavWpGbNmnVMTU39NikpSRdzxIGCggLLz8/v9PXXX08CBha/X395yEE1agRvvAEXX+yXOfn5z2H37tCpRMqlc2pq6hYVVvxISkpyqampm/F7xwfQnpaUS926MGUKHHMM3H03fPUVTJsGDRuGTiZSpiQVVvwp/DcrcadKe1pSbklJ8Pvfw6RJfkj8aafB6tWhU4lIIlFpySEbPtyvgpyTAz/4AXz4YehEItHz9ddfJ3fo0KFThw4dOjVu3LhrkyZNuuy7vXPnznKdFf71r3/d1My6r1u37vujYmPGjGmWlpbWOSMjo/Nzzz33/bGOd999t367du06paWldb7iiitaFxQUlPiah/r8HTt2WP/+/Y9JS0vr3KVLlw7Lly+vs+85Dz/8cKP09PTO6enpnR9++OFG+7YvW7asTpcuXTqkp6d37t+//zH7/nsLCgq44oorWqelpXVu165dp/fee69++X+inkpLKuTss/0w+Fq1oHdvmD07dCKRaGnWrNneZcuWfbZs2bLPLrvssvyRI0fm7btdt27dgx6yXLFiRe25c+c2bN68+feTqmVnZ9edMWPG0cuXL//01Vdf/fyGG25I27NnDwDXXHNN+oQJE3JWrVq1ZOXKlXWnT59+wMH7ijz/oYceapySkrLnq6++WjJq1Ki8m266qRVAXl5e8tixY1t88MEHSxctWrR07NixLfLz85MBbrrpplajRo3Ky8nJWZKSkrLnoYceagzw7LPPpqxcubLuqlWrlkycODHnmmuuSTvUn6tKSyqsc2dYsADat/cT7mqWeJGqM2rUqNbjxo1bY0WG6k6fPv3IH/3oRxvr1avnOnTosCs9Pf27t99++/CcnJza27ZtSzr77LO3JyUlcemll37zwgsvHFX8NSvy/FmzZh155ZVXfgMwbNiwb99///0GBQUFvPDCCym9e/fe0rRp072pqal7e/fuvWXGjBkpBQUFzJ8/v8GwYcO+Bbjyyiu/mTlz5pEAL7744pGXXnrpN0lJSfTt23f7li1bauXk5NQ+lJ+LBmJIpTRv7heRvOQSP0v80qXw4INQ+5A+hiKxd+WVtF6yhEM+HFWWzp35z+TJHPKZ3e7du7ffvn17cvHt99xzz+rBgwdvnTJlSkrz5s13n3LKKTuK3p+bm1unZ8+e2/bdbtGixa7Vq1fXqVOnjmvevPn3Y3rT09N3rVu37oD/F1bk+Xl5eXXatGmzC6B27docccQRe/Py8mrl5ubWbtWq1fd7gS1bttyVm5tbOy8vr1aDBg321i78JZCRkbErLy+vDsC6detqZ2RkfP+c5s2b78rJyamdnp5e7vHIkSwtMzsfOP/YY48NHUXK4fDDYcYMuO02uPde+PxzP7LwqAP+zhMRgOzs7OWl3bd169aksWPHNn/rrbcOmEDNlbBukJm5UrYfsK0izz/U55T2+DJe64BtZYlkaTnnZgIzMzMzrwqdRconORnGjoWOHWHECOjZE2bOhHbtQicT8SqyRxQrZe1ptWjRYs+aNWsO69KlSyfwezrdunXruGDBgqWtWrXatXr16u8HQqxdu7ZOq1atdmdkZOwuumeVk5NTp1mzZgfsvVTk+c2aNdv15Zdf1mnbtu3u3bt3s23btuQmTZrsbdWq1e558+Z9vw5Ebm5unTPOOGNrs2bN9mzdujV59+7d1K5dm1WrVtVp0qTJboAWLVrsXrVq1ffvv27dujppaWmHdNWnzmlJlbriCpg7FzZu9CML58wJnUgkerKzs5fvG5RR9Gvw4MFbe/TosWPjxo2Lc3NzP8nNzf2kadOmuz788MOlaWlpe4YMGbJpxowZR+/YscOWLVtWZ9WqVXX79OmzPT09fffhhx9eMGfOnMMLCgqYMmVKo0GDBm0q/r4VeX7//v03TZ48uRHA448/ftQpp5yyNSkpicGDB2+eN29ew/z8/OT8/PzkefPmNRw8ePDmpKQkevbsufXxxx8/CmDy5MmNBgwYsAlg4MCBm6ZMmdKooKCAOXPmHN6gQYO9h3JoECK6pyXx7bTT4IMP4Pzz4Zxz4NFH/SwaIlI5mZmZOwcPHryxXbt2xycnJ/PAAw/k1Krlf41PmDAhZ/jw4W127txpZ5555pYLL7xwM8CUKVNSFi5cePiDDz64tiLPHz169IYhQ4a0SUtL65ySkrL3mWee+TdA06ZN9/7yl79c2717944At9xyy9qmTZvuBbj//vvX/PjHP2571113tTz++OP/M3r06A0AF1100ebZs2enpKend65Xr17BpEmTVh3qz8BKOsYYFZmZmW7RokWhY0gFbdkCQ4fCyy/D9dfD/ff7IfIisWBm2c65zKLbFi9evKpr164bQmWSilu8eHHjrl27ZhTfrsODEjMNG/p1uW68EcaPhwEDQOt6ikhlqLQkppKT4YEH4LHH/Pmtnj396EIRkYpQaUm1uOoqP1P8hg1+na5XXgmdSBJEQUFBgRbSiTOF/2YlzkMVydIys/PN7LHNOpZUo/TpAwsXQkYG9O/vh8hH+JSq1AxL8vPzU1Rc8aNwPa0UYElJ92sghlS77dv9pLvPPOPX6PrrX6F+lc5TIImopIEYWrk4LmnlYomWww+HqVPhxBP9LBrLlsELL0B6euhkUtMU/tI7YPVbiV/6y0OCMPOrIM+aBV9+CZmZMG9e6FQiEnUqLQnqvPP8hciNG/vlTh59VOe5RKR0Ki0Jrl07v8TJuefCqFF+pOF334VOJSJRpNKSSGjY0J/Xuv12PzCjTx9YuzZ0KhGJmkiWloa8J6akJPjd72D6dPjkE+jWza/VJSKyTyRLyzk30zk3IiUlJXQUCWDIEH+eKyUFzjrLz6ih81wiAhEtLZFOnfyFyAMHwi9+4a/n2rbt4M8TkZpNpSWR1bAhPPecnzlj+nQ//dOyZaFTiUhIKi2JNDO45Zb98xaefLIvMhFJTJEsLQ3EkOLOOguys/1hwwsu8EW2Z0/oVCJS3SJZWhqIISVp3dqPJrz6ahg3Dvr1g/UHzEwmIjVZJEtLpDSHHQYTJsATT8D8+X5Y/Pz5oVOJSHVRaUlcuvxyX1Z16kDv3hoWL5IoVFoSt048ET78EAYM8MPif/hD+Pbb0KlEJJYiWVoaiCHldeSRMGMG/PGPMHu2P1y4cGHoVCISK5EsLQ3EkENhBjfcAO++CwUF0KsXPPywDheK1ESRLC2RiujZEz76CM45B66/Hi66CLSzLlKzqLSkRjn6aHjxRbj3Xnj+eeje3Z/3EpGaIZKlpXNaUhlJSfDLX/qVkHfuhFNOgT/9SYcLRWqCSJaWzmlJVejVC/71Lz+bxtVXwyWXwNatoVOJSGVEsrREqkrjxn5U4d13w7RpfnThokWhU4lIRVVbaZnZMWb2VzObXo7H6vCgVJmkJBgzBt5+G777Dk49Fe6/3480FJH4UqnSMrPJZrbezJYU255lZsvNbIWZ3QrgnFvpnBtentfV4UGJhdNP94cLBwyAm2+G/v0hLy90KhE5FJXd03oCyCq6wcySgUeBc4FOwFAz61TJ9xGpEkcf7Zc2mTjR73l17eqXPRGR+FCp0nLOvQNsLLa5B7CicM9qF/APYFB5X9PMRpjZIjNblJ+fX5l4IiUyg5Ej/cwZjRr52eL/+79h167QyUTkYGJxTqslsLrI7TVASzNrZGZ/Ak4yszGlPdk595hzLtM5l5mamhqDeCJe586+uH7+c39d1+mnw8qVoVOJSFliUVpWwjbnnPvGOTfSOdfWOfeHGLyvyCGrX99fwzV9Onz+uZ+Ed+rU0KlEpDSxKK01QOsit1sBa2PwPiJVZsgQP0ijSxd/PdewYbBtW+hUIlJcLEprIXCcmbUxszrAxcBLMXgfkSqVnu4HZ9xxBzz5JJx0EixYEDqViBRV2SHvU4H5QHszW2Nmw51ze4BRwGvAUmCac+7TykcVib1ateDOO/0UULt3+1k17rwT9uwJnUxEAMxFeEK2zMxMt0jTF0ggmzfDqFHw9NN+Bvmnn4a2bUOnktKYWbZzLjN0DoktTeMkUoqUFPjb3/zAjGXL/DVdkydr4l2RkFRaIgdx8cXw8cfQowcMH+4HbWzYEDqVSGJSaYmUQ+vW8OabcN99fgLeLl3gtddCpxJJPCotkXJKSoJf/AI++MBPB5WVBaNHw44doZOJJA6Vlsgh6trVz6QxejSMHw+ZmfDRR6FTiSQGlZZIBdSrBw8+6A8RfvutP9/1u9/5YfIiEjsqLZFK6NcPliyBiy6CX//ar9X12WehU4nUXCotkUo6+miYMgWefRZWrfKrI99/P+zdGzqZSM2j0hKpIhdc4Pe6srL8IpNnnAErVoROJVKzqLREqlDTpvD88/DUU77AunaFCROgoCB0MpGaQaUlUsXM4Kc/9aV12mlw7bVwzjnw1Vehk4nEP5WWSIy0agWvvurX65o/H044AZ54QtNAiVSGSkskhsz8ysgff+wPFQ4bBgMHQm5u6GQi8SmSpWVm55vZY5s3bw4dRaRKHHOMX6vrgQdgzhw4/nhNvitSEZEsLefcTOfciJSUlNBRRKpMUhLceKPf6zrxRD/57jnnQE5O6GQi8SOSpSVSkx17LMyd60cVzp8PnTtrhKFIeam0RAJISoKrr/YjDE891Y8wPPNMXdclcjAqLZGA0tP9CMPJk2HxYr/kyQMPaDYNkdKotEQCM/OjCj/7DM4+2y9/0quX5jAUKYlKSyQiWrSAF1+Ev//dHyY86ST4/e81c7xIUSotkQgxg6FD/V7W4MFw++1+va4FC0InE4kGlZZIBDVpAs88Ay+8AN98A6ecAtdfD1u3hk4mEpZKSyTCBg3ye13XXguPPAKdOsHMmaFTiYSj0hKJuIYN4eGH4f334cgj/TRQF14I69aFTiZS/VRaInGiZ0/IzvaDM2bOhI4d4c9/1kXJklgiWVqae1CkZHXqwG23wSef+BWSR470i00uXRo6mUj1iGRpae5BkbIdd5yfePfxx/05r65d4Te/ge++C51MJLYiWVoicnBmcMUVfi/roovgt7/1E/G+9VboZCKxo9ISiXNNmsDTT8Mrr/g9rbPO8isn5+WFTiZS9VRaIjVEVhZ8+inccQdMmwbt2/vZ4zWPodQkKi2RGqRePbjzTr9mV2amv76rZ09YtCh0MpGqodISqYHat4c33oCpU2HNGujRA0aNgk2bQicTqRyVlkgNZQYXXwzLlvnCmjgROnSAKVPAudDpRComkqWl67REqk5KCowfDwsX+vW7fvIT6NvXl5lIvIlkaek6LZGq162bnwrqT3+Cjz7yC07+6lewfXvoZCLlF8nSEpHYSE6Gn/8cli+HSy6Bu+/200E9+6wOGUp8UGmJJKAmTeCJJ+Ddd6FRI39xct++fsi8SJSptEQS2Gmn+eHwEybA4sV+OqgbbtAoQ4kulZZIgktOhquvhs8/h6uu8oM22reHyZM1g7xEj0pLRAB/mHDiRL/ndeyxMHy4XzH5gw9CJxPZL5KlpSHvIuF06wbvvQdPPQVffQU/+IEvsPXrQycTiWhpaci7SFhmftLd5cvh5pt9gbVrBw89BHv2hE4niSySpSUi0dCwIYwb5xed/MEP/CCNrl3htddCJ5NEpdISkYPq0AFefRVeeMEvf5KVBf37a1YNqX4qLREpFzMYNMhfy3Xfff68V+fOMHo0bNwYOp0kikiWlgZiiETXYYfBL34BK1b4IfKPPOJHG44fD7t3h04nNV0kS0sDMUSiLzXVD5H/17+ge3e/x3XCCfDyy5oSSmInkqUlIvHjhBPg9ddh5kxfVv37w7nnakooiQ2VlohUmhkMGOBHGf7xj7BggR9leO21sGFD6HRSk6i0RKTK1Knjh8WvWOGnhvrzn/35rnHjYOfO0OmkJohkaWkghkh8a9QIHn4YPv4YevWCW27Zv2qy5jOUyohkaWkghkjN0KkTzJ4Nc+b4IvvJT+Dkk2Hu3NDJJF5FsrREpGY56yxYuBCeftqf4+rbF847D5YsCZ1M4o1KS0SqRVISXHqpn89w3DiYP98P1vjZzyA3N3Q6iReRLC2d0xKpuerW9ZPwrljhr+166ik47ji44w7YujV0Oom6SJaWzmmJ1HyNGsEDD/g9r0GD4K67/EjDiRM1s4aULpKlJSKJo00bmDrVLzbZsSNcc42/YPmLL0InkyiqttIys8Fm9hcze9HM+h3ksTo8KJJgTj4Z3noLXnoJ2raFtLTQiSSKzFVikjAzmwwMANY75zoX2Z4FPAQkA5Occ/cUue8o4D7n3PCDvX5mZqZbtGhRhfOJSOIws2znXGboHBJbld3TegLIKrrBzJKBR4FzgU7AUDPrVOQhtxfeLyIickgqVVrOuXeA4ivp9ABWOOdWOud2Af8ABpk3FnjFOfdhaa9pZiPMbJGZLcrPz69MPBERqWFicU6rJbC6yO01hduuA84GLjCzkaU92Tn3mHMu0zmXmZqaGoN4IiISr2rF4DWthG3OOTceGB+D9xMRkQQRiz2tNUDrIrdbAWtj8D4iIpJgYlFaC4HjzKyNmdUBLgZeisH7iIhIgqlUaZnZVGA+0N7M1pjZcOfcHmAU8BqwFJjmnNMapiIiUmmVOqflnBtayvaXgZcr89oiIiLFVeri4lgzs3xgOxClBbsbE608EL1MynNwUctUE/KkO+c05LiGi3RpAZjZoihd5R61PBC9TMpzcFHLpDwSLzRhroiIxA2VloiIxI14KK3HQgcoJmp5IHqZlOfgopZJeSQuRP6cloiIyD7xsKclIiICqLRERCSOBC0tM2ttZm+Z2VIz+9TMRhe57zozW164/d4i28eY2YrC+86prkxmdqKZ/dPM/lW4dEqP6shkZnXN7AMzW1yY57eF2482szfM7IvC/z0qcJ5xZrbMzD42s+fN7MjqyFNWpiL332xmzswaV0emsvKE+FyX8W8W5DNd5D2SzewjM5tVeDvIZ1rijHMu2BfQHOhW+H0D4HP8wpFnAm8ChxXe16TwfzsBi4HDgDbAv4Hkasr0OnBu4fbzgLerIxN+1vwjCr+vDSwAegL3ArcWbr8VGBs4Tz+gVuH2sdWVp6xMhbdb46cUywEaB/4ZBflcl5EnyGe6SK6bgL8DswpvB/lM6yu+voLuaTnn1rnCBSGdc1vxcxW2BK4G7nHOfVd43/rCpwwC/uGc+8459yWwAr/oZHVkckDDwoelsH/m+phmct62wpu1C79c4fs+Wbj9SWBwyDzOudedn3cS4J/42f1jnqesTIW3/wjcUuR2zDOVkSfI57qMPEE+0wBm1groD0wqsjnIZ1riS2TOaZlZBnAS/q/AdsDpZrbAzOaZ2cmFDyttgcnqyHQDMM7MVgP3AWOqK1PhYZR/AeuBN5xzC4Cmzrl14IsWaBI4T1FXAq9UV57SMpnZQCDXObe42MND/YyCfa5LyXMDgT7TwIP4PyYKimwL9pmW+BGJ0jKzI4DngBucc1vwE/kehT+E8UtgmpkZpSwwWU2ZrgZudM61Bm4E/rrvobHO5Jzb65w7Eb/30sPMOpfx8KB5zOxXwB5gSnXlKSVTF+BXwK9LeHion1Gwz3UpeYJ8ps1sALDeOZdd3qfEMo/El+ClZWa18eUwxTk3o3DzGmBG4WGND/B/jTWmmhaYLCXT5cC+759l/+GJalv00jm3CXgbyALyzKx5Yd7m+L+gQ+bBzC4HBgCXOuf2/VKp1kVBi2QahD//sdjMVhW+74dm1qw6MxX7GQX9XJeQJ9RnuhcwsPDf5R/AWWb2NBH4TEscCHUyrfB3mgFPAQ8W2z4SuLPw+3b4QwMGHM//PyG7ktic1C8p01KgT+H3fYHswu9jmglIBY4s/L4e8C6+GMbx/09a3xs4TxbwGZBa7PHV8W9WYqZij1nF/oEYoX5GQT7XZeQJ8pkulq0P+wdiBPlM6yu+viq1nlYV6AX8FPik8Hg7wG3AZGCymS0BdgGXO+cc8KmZTcP/ctwDXOuc21tNma4CHjKzWsBOYASAcy7WmZoDT5pZMn7PeJpzbpaZzccfXhoOfAVcGDjPCvwvlTf8ES/+6ZwbWQ15Ss1U2oMD/ozqEOZzXVqeTYT5TJfmHsJ8piWOaBonERGJG8HPaYmIiJSXSktEROKGSktEROKGSktEROKGSktEROKGSktEROKGSktEROLG/wHHo71ooyOiAwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAD4CAYAAABfYrnHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjaUlEQVR4nO3deXxW1Z3H8c8vCSmIGhDCGpLgwi7IkOJWcUP2CsoIMlJlUaqAgrhUrbV1qlNplVERUETQWpDFIiDiUrWijrSQ4IYsymAiOxFRREWWnPnjBJuJScj25N6bfN+vV17kuXnuc3+X+4Jvzj3nnmPOOURERKIgLugCRERESkuhJSIikaHQEhGRyFBoiYhIZCi0REQkMhKCLqAkDRs2dOnp6UGXISIRkJWV9blzLjnoOiS2Qh1a6enpZGZmBl2GiESAmeUEXYPEnm4PiohIZIQytMzs52Y2/auvvgq6FBERCZFQhpZz7nnn3KikpKSgSxERkRAJZWiJiIgUJZShpduDIiJSlCoLLTM70cyeMLNnj/Ze3R4UEZGilCq0zGymme0yszWFtvcysw1mttHMbivpM5xzm5xzI0t5PLW0RETkR0rb0noS6FVwg5nFA1OA3kA7YIiZtTOzU81saaGvRmUpqqwtLedg8mR4/fWyHEVERKKmVA8XO+feNLP0Qpu7Ahudc5sAzGwu0N859wegX3kLMrNRwCiA1NTUUu2zfz9MnQp79sB770GTJuU9uoiIhFlF+rSaA5sLvN6Sv61IZtbAzB4FOpvZ7cW9zzk33TmX4ZzLSE4u3YwsderAggWwdy/8x3/A4cOlPAMREYmUioSWFbGt2GWQnXO7nXPXOudOym+NVaoOHXxr6+9/h7vvruxPFxGRMKhIaG0BWhR4nQJsq1g5FTNsmP+65x545ZUgKxERkVioSGitAk4xs5ZmlghcDiypnLLKb8oUaNcOhg6FbYFGqIiIVLbSDnl/BlgBtDazLWY20jl3CBgLvAysA+Y75z6KXamlc8wxvn/r22/h8svh0KGgKxIRkcpSqtByzg1xzjV1ztVyzqU4557I377MOdcqv5/q3tiWWnpt28Kjj8Jbb8FddwVdjYiIVJZQTuNUGYYOhauvhj/8AV58MehqRESkMlTb0AJ4+GHo2BF+8QvYvPno7xcRkXCr1qF15Pmt77/3/VsHDwZdkYiIVES1Di2AVq1gxgx45x24446gqxERkYqo9qEFMHgwXHcd3H8/PP980NWIiEh51YjQApg0CTp3hquugpycoKsREZHyqDGhVbu27986fBgGDYIDB4KuSEREyqrGhBbASSfBzJmwciX86ldBVyMiImUVytCK5SKQAwfC9dfDgw/Cc89V+seLiEgMmXPFTsweuIyMDJeZmVnpn/v993DOOfDxx7B6NZx4YqUfQkSqmJllOecygq5DYiuULa1Y+8lPYN48MPP9W99/H3RFIiJSGjUytABatoRZsyArC26+OehqRESkNGpsaAEMGAA33giPPOJHFoqISLjV6NACuO8+OP10GDkSNm4MuhoRESlJjQ+txESYPx9q1YLLLoP9+4OuSEREilPjQwsgNRX+/Gd47z0YPz7oakREpDihDK1YPqdVnL594dZb4bHH4JlnquywIiJSBjXyOa3iHDwI558P778PmZnQunWVHVpEKkjPadUMoWxpBaVWLZg7189TeNll8O23QVckIiIFKbQKSUmBp5+GDz+EG24IuhoRESlIoVWEXr38gpFPPOEDTEREwkGhVYy774Zzz4Vrr4W1a4OuRkREIKShFcTowcISEmDOHKhb1/dvffNNYKWIiEi+UIaWc+5559yopKSkQOto1swH17p1MGZMoKWIiAghDa0w6d4dfvMbeOopP8GuiIgER6FVCnfdBRdc4FtbH34YdDUiIjWXQqsU4uP9bcKkJN+/9fXXQVckIlIzhTK0wjAQo7DGjf30Tp984kcUhngiERGRaiuUoRWWgRiFnXeeHwo/Zw48/njQ1YiI1DyhDK0wu+MO6NHDz5bx3ntBVyMiUrMotMooLg7+8hdo0MD3b+3dG3RFIiI1RyhDK4x9WgUlJ/uJdT/9FK65Rv1bIiJVJZShFdY+rYLOOQfuucevejxtWtDViIjUDKEMrai49Vbo0wduvBGysoKuRkSk+gtlaIX99uARcXF+poxGjWDQIAh5uSIikRfK0IrC7cEjGjaEefPgs89gxAj1b4mIxFIoQytqzjoL7rsPFi6EyZODrkZEpPoKZWhF5fZgQRMmwMUXw803w8qVQVcjIlI9maui+1lmVheYChwA3nDOzT7aPhkZGS4zMzPmtVWWPXugc2f//erVcMIJwdYjUpOYWZZzLiPoOiS2KtTSMrOZZrbLzNYU2t7LzDaY2UYzuy1/86XAs865a4CLj/K5kWtpAdSv74fAb9sGw4erf0tEpLJV9Pbgk0CvghvMLB6YAvQG2gFDzKwdkAJszn/b4ZI+NEoDMQrr2hX+9CdYsgQmTQq6GhGR6qVCoeWcexP4otDmrsBG59wm59wBYC7QH9iCD64Sj2tmo8ws08wyc3NzK1JeYG64AS69FG67DVasCLoaEZHqIxYDMZrzrxYV+LBqDiwEBprZNOD54nZ2zk13zmU45zKSk5NjUF7smcETT0CLFjB4MOzeHXRFIiLVQyxCy4rY5pxz3zjnhjvnrivNIIyoq1cPFiyAnTvhyishLy/oikREoi8WobUFaFHgdQqwLQbHCb0uXXy/1rJlvp9LREQqJhahtQo4xcxamlkicDmwJAbHiYTRo/0SJr/+Nbz9dtDViIhEW0WHvD8DrABam9kWMxvpnDsEjAVeBtYB851zH1W81GgygxkzoGVL378V0bElIiKhkFCRnZ1zQ4rZvgxYVpHPrk6OP973b51xBgwdCi++6CfbFRGRstF/nVXktNPg4YfhlVfgjjuCrkZEJJoq1NKSsrnmGj+908SJ0KyZf55LRERKT6FVhcxgyhQ/DH78eGjSxK/DJSIipaPbg1UsPh7mzIGf/Qx+8Qt4/fWgKxIRiQ6FVgDq1IHFi6FVKxgwAN59N+iKRESiQaEVkPr1/SjCevWgd2/YtCnoikREwk+hFaCUFHj5ZTh4EHr2hF27gq5IRCTcFFoBa9sWli6FrVuhb1/Yty/oikREwiuUoRXVRSDL68wz/eKR774LAwfCgQNBVyQiEk6hDK0oLwJZXv36wWOP+YePR4zQrPAiIkXRc1ohMnIk7NgBd94JTZtqZngRkcIUWiFzxx2wfTvcf78PrgkTgq5IRCQ8FFohYwYPPeRnzbjpJmjcGK64IuiqRETCQaEVQvHx8PTT8PnnMGwYJCdDjx5BVyUiErxQDsQQqF0bFi2C9u3hkku0gKSICCi0Qi0pyT98nJICffrAqlVBVyQiEqxQhlZNe06rJI0bw2uvQcOGftaMDz4IuiIRkeCEMrRq4nNaJUlJ8cF1zDFw0UWwfn3QFYmIBCOUoSU/1rKlDy6ACy/UBLsiUjMptCKkdWt49VXYv98H1+bNQVckIlK1FFoRc+qpfqqnL77wwbVjR9AViYhUHYVWBHXp4tfi2rbN93F9/nnQFYmIVI1QhpZGDx7dWWfBkiXwySc+uHbvDroiEZHYC2VoafRg6VxwASxeDOvW+eD64ougKxIRia1QhpaUXs+efuaMtWuhe3cFl4hUbwqtaqBXLwWXiNQMCq1qQsElIjVBKENLAzHKR8ElItVdKENLAzHKr2BwnX8+7NoVdEUiIpUnlKElFdOrFzz/vB8Of+65/nkuEZHqIJShpduDFXfRRfDSS7BlC3TrBjk5QVckIlJxoQwt3R6sHN26+bkKd++Gc87xLS8RkSgLZWhJ5Tn9dPj73+G773yIrV0bdEUiIuWn0KoBTjsNli/33597Lrz7bqDliIiUWyhDS31ala9dO3jzTahTx0//9M9/Bl2RiEjZmXMu6BqKlZGR4TIzM4Muo1rJyfFLmuzcCUuX+paXSHVgZlnOuYyC27KysholJCTMADoQ0l/S5UfygDWHDh26ukuXLj96aCchgIKOysx+Dvz85JNPDrqUaictzbe4uneH3r3h2WehT5+gqxKJjYSEhBlNmjRpm5ycvCcuLi68v6HLD/Ly8iw3N7fdjh07ZgAXF/55KH/z0OjB2GrWzPdxtW0L/fvD7NlBVyQSMx2Sk5P3KrCiIy4uziUnJ3+Fbx3/+OdVVYiZnWhmT5jZs1V1TClecrIfVfizn8HQofDww0FXJBITcQqs6Mm/ZkXmU4VCy8xmmtkuM1tTaHsvM9tgZhvN7DYA59wm59zIUn6uBmJUgeOP9ysgDxgA48bBXXdBiLs4RUQq3NJ6EuhVcIOZxQNTgN5AO2CImbUry4fq9mDVqV0bFiyAESPg97+HMWPg8OGgqxKJvh07dsS3adOmXZs2bdo1bNiwU6NGjToeeb1//34rad8JEyY0K/j+efPm/fCf4e23394kNTW1Q3p6eoe//vWvxx/Z/tZbbx3TqlWrdqmpqR2GDRvWIi8vr8jPLuv+3333nfXt2/fE1NTUDh07dmyzYcOGxCP7TJ48uUFaWlqHtLS0DpMnT25wZPv69esTO3bs2CYtLa1D3759Tzxyvnl5eQwbNqxFampqh1atWrV7++23jynr32uFQss59yZQeC7xrsDG/JbVAWAu0L+0n2lmo8ws08wyc3NzK1KelFJCAsyYAbfeCtOmwRVXwIEDQVclEm1NmjQ5vH79+rXr169fe+WVV+Zee+21O4+8rl279lHvaRR8/+DBg78CyMrKqr1w4cITNmzY8NFLL7308fjx41MPHToEwOjRo9OmTp2ak52dvWbTpk21n3322eMLf2Z59n/ooYcaJiUlHfrss8/WjB07dueECRNSAHbu3Bk/ceLEZitXrlyXmZm5buLEic1yc3PjASZMmJAyduzYnTk5OWuSkpIOPfTQQw0BFixYkLRp06ba2dnZa6ZNm5YzevTo1LL+vcaiT6s5sLnA6y1AczNrYGaPAp3N7PbidnbOTXfOZTjnMpKTk2NQnhTFDCZOhD/9CebNg759Ye/eoKsSkYKeffbZepdeeukXderUcW3atDmQlpb2/RtvvFE3Jyen1r59++K6d+/+TVxcHFdcccXuRYsW1a+M/ZcuXVpvxIgRuwGGDx++55133jkuLy+PRYsWJXXr1m1v48aNDycnJx/u1q3b3oULFybl5eWxYsWK44YPH74HYMSIEbuff/75egCLFy+ud8UVV+yOi4vjwgsv/Gbv3r0JOTk5tcrydxCLIe9FNXudc243cG0MjieV6Oab/SCNq6/28xW++KIfbSgSdSNG0GLNGsp8O6okHTrw7cyZ/++X9FLp0qVL62+++Sa+8Pb77rtv84ABA74GeOKJJxrNnTu3QadOnb6dOnXq5uTk5MNbt25NPOOMM/YdeX+zZs0ObN68OTExMdE1bdr04JHtaWlpB7Zv3/6jMCjP/jt37kxs2bLlAYBatWpx7LHHHt65c2fC1q1ba6WkpPxwT6Z58+YHtm7dWmvnzp0Jxx133OFatfzh09PTD+zcuTMRYPv27bXS09N/2Kdp06YHcnJyaqWlpf1w7KOJRWhtAVoUeJ0CaHGMCLnqKmjaFAYOhDPO8LPFtytTr6SIlCQrK2tDST+/8cYbd/3xj3/cZmaMHz+++ejRo1ssWLAgu6jJIMzMFbP9R9vKs39Z9ynu/SV81o+2lSQWobUKOMXMWgJbgcuB/4jBcSSGevTwDyH36QNnnw2LF/sJd0Wiqjwtolg5WkurRYsWh45sGzt2bG6/fv1OAUhJSTmwefPmHwZCbNu2LTElJeVgenr6wYItq5ycnMQmTZr8qPVSnv2bNGly4NNPP0086aSTDh48eJB9+/bFN2rU6HBKSsrB5cuXH3dkn61btyaee+65Xzdp0uTQ119/HX/w4EFq1apFdnZ2YqNGjQ4CNGvW7GB2dvYPx9++fXtiampqqVtZUPEh788AK4DWZrbFzEY65w4BY4GXgXXAfOfcRxU5jgSjc2dYsQKaNPHrc82fH3RFItVDVlbWhiODLAp+Hbk1WLCfZ+7cufVat279HcDAgQO/XLhw4QnfffedrV+/PjE7O7v2eeed901aWtrBunXr5r322mt18/LymD17doP+/ft/Wfi45dm/b9++X86cObMBwKxZs+qfeeaZX8fFxTFgwICvli9ffnxubm58bm5u/PLly48fMGDAV3FxcZxxxhlfz5o1qz7AzJkzG/Tr1+9LgIsvvvjL2bNnN8jLy+O1116re9xxxx0uy61BqGBLyzk3pJjty4BlFflsCYf0dPif//EzZwweDFu3wo03Bl2VSPU2bty4lLVr19YB3zqaNWtWDkBGRsb+AQMGfNGqVav28fHxTJo0KSchwf83PnXq1JyRI0e23L9/v51//vl7L7vssq8AZs+enbRq1aq6Dz744Lby7D9u3LjPBw4c2DI1NbVDUlLS4Xnz5v0vQOPGjQ/fcsst27p06dIW4NZbb93WuHHjwwAPPPDAlsGDB590zz33NG/fvv2348aN+xxg0KBBX73wwgtJaWlpHerUqZM3Y8aM7LL+3WjCXCmV/fv9zBl//at/EPmBByD+Rzc3RIJT1IS577//fnanTp0+D6omKb/333+/YadOndILbw/l3IMSPrVr+6Hw48bBQw/BpZfCvn1H309EpDIptKTU4uPhwQdh8mS/rMk558CWLUFXJSI1iUJLymzsWHj+edi4EU4/HVavDroikWLl5eXllW1MtQQu/5oVOQ+VQkvKpU8fP0AjPt63uBYvDroikSKtyc3NTVJwRUf+elpJwJqifh7KRSAlGjp2hJUr4eKL4ZJL/BRQEyb4KaFEwuDQoUNX79ixY8aOHTu0cnF0/LBycVE/VGhJhTRpAm+84WfRuPlm+PhjeOQRqFWm2cREYiN/ufYfrX4r0aXfPKTCjjnGjyy84w6YPh169oTdu4OuSkSqo1CGlhaBjJ64OLj3XnjqKXjnHfjpT+HDD4OuSkSqm1CGlhaBjK4rr4Tly/3DyGeeCc89F3RFIlKdhDK0JNpOPx0yM6F9e/8Q8n/+JxSziKqISJkotCQmmjXzLa5f/AJ++1sYNEgzaIhIxSm0JGZq1/Z9XPff728Tnn02ZGcHXZWIRJlCS2LKDG66CZYtg5wcyMiAv/0t6KpEJKoUWlIlevaEVav8isi9esEf/gAhXmBAREJKoSVV5pRT4B//8P1bd9zhB2ns3Rt0VSISJaEMLT2nVX3VrQtz5sB//7efdPenP4W1a4OuSkSiIpShpee0qjczGD8eXn8dvvoKunaFBQuCrkpEoiCUoSU1Q7dukJXlJ94dNMjPXXjoUNBViUiYKbQkUM2b+wl3x4yBBx6A7t1h+/agqxKRsFJoSeASE/3M8E8/7UcYnnYavPpq0FWJSBgptCQ0hg71odWwIfToAb/7HRw+HHRVIhImoQwtjR6sudq18wtLXnkl3H23D68dO4KuSkTCIpShpdGDNVvduvDkkzBzJqxYAZ07w9//HnRVIhIGoQwtEYDhw32rq149P0Dj97/X7UKRmk6hJaHWoYPv5xoyBO66y08HtW1b0FWJSFAUWhJ6xx7rRxbOmOFXRe7UCZYuDboqEQlCKENLAzGkMDMYORJWr4aUFPj5z+GGG/wKySJSc4QytDQQQ4rTpo2fdHf8eJg82U8BpbkLRWqOUIaWSEl+8hM/4e4LL/jh8F26wKOPaqkTkZoglKGl24NSGn36wAcf+DkMr7vOL3Wye3fQVYlILIUytHR7UEqrSRN48UW4/37f8urUSc90iVRnoQwtkbKIi4ObbvIPItetCxdcABMmaJCGSHWk0JJqo0sXP7pw9Gjf59WlC7z7btBViUhlCmVoqU9LyqtuXZgyBV56Cfbs8aML771X63SJVBehDC31aUlF9ewJa9bAwIFw551+sMbGjUFXJSIVFcrQUktLKsMJJ8DcuTBnDqxb5wdpPPaYhsaLRFkoQ0stLalMQ4bAhx/C2WfDtddC375aHVkkqqostMxsgJk9bmaLzaxHVR1XBPzUTy+95GfReOMNv27Xn/+sVpdI1FQotMxsppntMrM1hbb3MrMNZrbRzG4DcM4tcs5dAwwDBh/lc3V7UCpdXByMHQvvvQft28NVV0G/frB1a9CViUhpVbSl9STQq+AGM4sHpgC9gXbAEDNrV+Atd+b/vFi6PSix1KoVLF8ODz7oH0Ru3x5mzVKrSyQKKhRazrk3gS8Kbe4KbHTObXLOHQDmAv3Nmwi86JxbXdxnmtkoM8s0s8zc3NyKlCdSrPh4GDfOTwPVqROMGAG9e8PmzUFXJiIliUWfVnOg4D/9Lfnbrge6A/9uZtcWt7NzbrpzLsM5l5GcnByD8kT+5eSTfWvrkUfg7bd9q+vxx9XqEgmrWISWFbHNOeceds51cc5d65x7NAbHFSmXuDgYM8aPMMzIgFGjoEcPyM4OujIRKSwWobUFaFHgdQqgBdIl9Fq2hFdfhWnT/Jpd7dvDpEmaTUMkTGIRWquAU8yspZklApcDS2JwHJFKFxfnn+X66CM/8e5NN8Hpp/s5DUUkeBUd8v4MsAJobWZbzGykc+4QMBZ4GVgHzHfOfVTxUkWqTmoqLFkC8+f7IfFdu8Itt8A33wRdmUjNZi7EPc4ZGRkuMzMz6DKkhtuzB371Kz9Ao2VLf/uwZ8+gq5LCzCzLOZcRdB0SW6GcxkkkTOrXh+nT/bNdiYnQqxcMHQq7dgVdmUjNo9ASKaVu3fxsGnfd5W8btm2rh5JFqppCS6QMateGu+/2i0u2aeMfSu7WzQ+XF5HYU2iJlEP79vDWWzBjhl/2pHNnP9Lw66+DrkykelNoiZRTXByMHAkbNvgW16RJvvU1f75uGYrEikJLpIIaNPADNVasgMaNYfBgP6PGxx8HXZlI9aPQEqkkZ5wBq1b5NbtWrYJTT4U774Rvvw26MpHqI5ShpfW0JKri4/2aXRs2+BbXvff6/q9Fi3TLUKQyhDK0tJ6WRF3jxn5l5DfegLp14ZJL/APJa9cGXZlItIUytESqi3PP9c92Pfywv2XYsaNfx2vPnqArE4kmhZZIjCUkwPXXwyefwDXX+LW7TjkFHn0UDh8OujqRaFFoiVSRhg39vIWrV0OHDnDddfBv/+anhxKR0lFoiVSxTp38asnz58OXX8J558GgQZCTE3RlIuGn0BIJgBlcdhmsX++nhVq61D+Y/JvfaFYNkZIotEQCVKeOn4B3/XoYMADuuQdOPhkee0wrJosURaElEgKpqfDMM/CPf0CrVn715E6dYNkyPd8lUlAoQ0sPF0tNdfrp8OabsHAhHDgAffvCRRf5YfMiEtLQ0sPFUpOZ+YeRP/rIP9/13nt+lOGwYbBlS9DViQQrlKElIn6V5Ouvh40b4ZZbYO5cf+vwzjs1WENqLoWWSMjVqwcTJ/rBGpdc4uczPOkk3wr7/vugqxOpWgotkYhIT4fZs2HlSj+D/Lhx0Lq1n+NQM2tITaHQEomYn/4UXn0VXnnFz7Jx1VV+pOHixRppKNVfKENLowdFSmbmRxWuWgULFsDBg/45r7PO0rRQUr2FMrQ0elCkdMzg3//djzR8/HHYvNlPC9W7N7z7btDViVS+UIaWiJRNQgJcfbWfSf5Pf4J//tMPkx8yBD7+OOjqRCqPQkukGqlTB26+GTZtgl//GpYsgbZt/TNe//u/QVcnUnGhDC31aYlUTL16fh7DTZv8KMN58/xIw6uvhuzsoKsTKb9Qhpb6tEQqR+PGMGmSD68xY+Avf/ELUP7yl/DZZ0FXJ1J2oQwtEalcTZvCQw/52TVGjYJZs/xs8mPGaGooiRaFlkgNkpICU6b48BoxAqZP9+F1ww2wfXvQ1YkcXShDS31aIrGVmgqPPupHGw4dClOnwokn+v4vtbwkzEIZWurTEqka6ekwYwZs2OCHxx8Jr1/+0veDiYRNKENLRKrWSSfBzJm+5XX11fDkk35G+auu8hP1ioRFKENLtwdFgpGe7ltbn37q+7kWLIB27WDwYPjgg6CrEwlpaOn2oEiwmjXzQ+Wzs+G22+DFF/2kvP37+/kORYISytASkXBo1Aj+678gJwfuvhveegu6doVevfzEvJpVXqpaKENLtwdFwqV+fbjrLh9eEyf6yXjPOw/OPBOeew7y8oKuUGqKUIaWbg+KhNNxx8Gtt/rbhlOnQm4uXHqpn99wxgytpCyxV2WhZWZtzexRM3vWzK6rquOKSOWrUweuu84PlZ83D449Fq65xg/kmDgRdJNEYqVCoWVmM81sl5mtKbS9l5ltMLONZnYbgHNunXPuWmAQkHGUz9XtQZEISEiAQYMgM9OvpnzqqX7gRosWvkW2dWvQFUp1U9GW1pNAr4IbzCwemAL0BtoBQ8ysXf7PLgbeBl4r6UN1e1AkWszgwgvhlVdg9Wro2xceeABatvTTRa1bF3SFUl1UKLScc28CXxTa3BXY6Jzb5Jw7AMwF+ue/f4lz7izgiuI+08xGmVmmmWXm5uZWpDwRCUDnzvDMM/+anHfuXP+sV9++vjWmEYdSEbHo02oObC7wegvQ3MzOM7OHzewxYFlxOzvnpjvnMpxzGcnJyTEoT0SqQsuW8MgjfsTh737nbyFedJF/3mvWLNi/P+gKJYpiEVpWxDbnnHvDOXeDc+6XzrkpMTiuiIRQcjL89rc+vJ54wre0RoyAtDT/7NeuXUFXKFESi9DaArQo8DoF2BaD44hIhNSu7cPqgw/gb3+DjAzfAktNhZEjYc2ao36ESExCaxVwipm1NLNE4HJgSQyOIyIRZAbdu8MLL/gBGsOH+z6wU0+FHj38lFF6WFmKU9Eh788AK4DWZrbFzEY65w4BY4GXgXXAfOfcRxUvVUSqmzZtYNo02LzZTxe1Zg306QPt2/sZ50UKMxfioTwZGRkuMzMz6DJEpIocOADz58OcObBoESQmln5fM8tyzpX4DKhEXyincRKRmikx0a+kvGxZ2QJLag6FloiIRIZCS0REIkOhJSIikaHQEhGRyFBoiYhIZCi0REQkMkIZWlpPS0REihLK0NJ6WiIiUpRQz4hhZrlATtB1lFND4POgi6hEOp9w0/lAmnNO6xlVc6EOrSgzs8zqNKWMzifcdD5SU4Ty9qCIiEhRFFoiIhIZCq3YmR50AZVM5xNuOh+pEdSnJSIikaGWloiIRIZCS0REIkOhVQ5mVtvMVprZ+2b2kZndnb/9BDP7m5l9kv9n/QL73G5mG81sg5n1DK76HyvhfH5nZlvN7L38rz4F9gnt+RxhZvFm9q6ZLc1/Hcnrc0QR5xP165NtZh/m156Zvy3S10hiT31a5WBmBtR1zu0zs1rA28A44FLgC+fcfWZ2G1DfOfcrM2sHPAN0BZoBrwKtnHOHAzqF/6eE8+kF7HPO3V/o/aE+nyPMbAKQARzvnOtnZn8kgtfniCLO53dE+/pkAxnOuc8LbIv0NZLYU0urHJy3L/9lrfwvB/QHnsrf/hQwIP/7/sBc59z3zrlPgY34f3yhUML5FCfU5wNgZilAX2BGgc2RvD5Q7PkUJ/TnU4LIXiOpGgqtcsq/VfMesAv4m3Pun0Bj59x2gPw/G+W/vTmwucDuW/K3hUYx5wMw1sw+MLOZBW7VhP58gAeBW4G8Atsie30o+nwgutcH/C9Gr5hZlpmNyt8W5WskVUChVU7OucPOudOAFKCrmXUo4e1W1EfEpLByKuZ8pgEnAacB24EH8t8e6vMxs37ALudcVml3KWJbFM4nktengLOdc/8G9AbGmFm3Et4blXOSGFNoVZBz7kvgDXz/z04zawqQ/+eu/LdtAVoU2C0F2FZ1VZZewfNxzu3MD7M84HH+dTsm7OdzNnBxfp/JXOACM/sL0b0+RZ5PhK8PAM65bfl/7gKew9cf1WskVUShVQ5mlmxm9fK/rwN0B9YDS4Cr8t92FbA4//slwOVm9hMzawmcAqys0qJLUNz5HPnPI98lwJr870N9Ps65251zKc65dOBy4HXn3FAien2KO5+oXh8AM6trZscd+R7oga8/ktdIqk5C0AVEVFPgKTOLxwf/fOfcUjNbAcw3s5HAZ8BlAM65j8xsPrAWOASMCdmop+LO52kzOw1/GyYb+CVE4nyKcx/RvD7F+WOEr09j4Dk/cJUEYI5z7iUzW0X1ukZSyTTkXUREIkO3B0VEJDIUWiIiEhkKLRERiQyFloiIRIZCS0REIkOhJSIikaHQEhGRyPg/Popr6Ow4cgIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#**ToDo: Need to find way to plot all of these together and indicate ocean surface\n",
    "#        (Could do that prettying it up in Illustrator afterwards)\n",
    "psMax = 1.e10 # Essentially infinite ocean mass\n",
    "for Ts in [280.,300.,350.,400.,500.]:\n",
    "    pp0List = numpy.linspace(min(pp0Fun(Ts),psMax/p0),.1,1000)\n",
    "    TList = [T(pp0,Ts) for pp0 in pp0List]\n",
    "    c = Curve()\n",
    "    c.addCurve(pp0List)\n",
    "    c.addCurve(TList,'T=%f'%Ts)\n",
    "    c.switchXY = c.reverseY =c.YlogAxis = True\n",
    "    plot(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "136ce7ab",
   "metadata": {},
   "source": [
    "Now do a case with a finite ocean mass\n",
    "\n",
    "**ToDo: Should we move this to the post-runaway section?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "497d4765",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcsAAAEGCAYAAAAKdL4tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABHx0lEQVR4nO3deXxU1fn48c+ZZLKTBMhGEsIEshECCSbsIgpoQUUBtSraWiuK2Gq1xVZRu/1sa63ULt8WK6DFFrGW4lpXVFyQLQFCCIQ9gSSEJCSQEEgmmTm/PyahEbNNMpOZTJ736zUvMueee+7jgPPk3HsWpbVGCCGEEO0zuDoAIYQQwt1JshRCCCE6IclSCCGE6IQkSyGEEKITkiyFEEKITni7OoDeEBYWpk0mk6vDEEL0ATk5OZVa63BXxyHci0cnS6XUHGBOQkIC2dnZrg5HCNEHKKWKXB2DcD8efRtWa/2W1vqekJAQV4cihBCiD/PoZCmEEEI4giRLIYQQohMenSyVUnOUUs+fOXPG1aEIIYTowzw6WcozSyGEEI7Q55KlUmq4UmqVUmqdq2MRQgjRP/RqslRKvaCUKldK7bmofJZSar9S6pBS6pGO2tBaH9Fa39XF68ltWCGEED3W2z3LvwOzWhcopbyAvwCzgVTgVqVUqlJqtFLq7YteEfZcrDu3Yf95/XH++2iFPZcRQgjh4Xp1UQKt9WdKKdNFxeOBQ1rrIwBKqVeA67XWvwGu7e61lFL3APcAxMXFdfm8gM/3UHi4Bn5zc3cvLYQQwsO4wzPLGOB4q/fFzWVtUkoNVko9B4xVSj3aXj2t9fNa6yytdVZ4eNdXrqoMqMSruqHL9YUQQng+d1juTrVRpturrLU+BdzrrGBODWjAVCqjZ4UQQvyPO/Qsi4Ghrd7HAqUuioWAIF8izg6isaHJVSEIIYRwM+6QLLcDiUqpeKWUD3AL8KbLookIxsvqRcGmYpeFIIQQwr309tSRtcBmIFkpVayUuktr3QR8H3gf2Ae8qrXO7824Whs1wAhA02cnXBWCEEIIN9Pbo2Fvbaf8HeCd3oylXQkNQACHt29hLJNcHY0QQgg34A63Yd1KxNQRAFTKOgZCCCGaSbK8SOr0VMxGzRkGuzoUIYQQbkKS5UW8jd746HIy8lw2IFcIIYSbkWTZFkMFfvUDXB2FEEIINyHJsg0V4b6cM3Z9iTwhhBCeTZJlG4pNYXg3BMnCBEIIIQBJlm2KjBiFsUlRut/s6lCEEEK4AUmWbagYdhiAzZ/nuTgSIYQQ7kCSZRvCYmyr+JRsK3BxJEIIIdyBJMs2DEuzbel1vlyeWQohhJBk2aYJM1IxezXSdKat3cOEEEL0Nx6dLJVSc5RSz585Y9/adUajkcqganyrrE6KTAghRF/i0clSa/2W1vqekBD7N3PWAd7E1YY7ISohhBB9jUcny544HeGNX0OAq8MQQgjhBiRZtmNkRTlhp6CxodHVoQghhHAxSZbtOBVWisbI/o/3ujoUIYQQLibJsh3VSWEAZH9xwMWRCCGEcDVJlu1Qo6IAqD5S7eJIhBBCuJoky3ZMnpwIQGx+mYsjEUII4WqSLNuRMj0VhRnfUzLXUggh+jtJlu3w8vaiNthMVeBIV4cihBDCxTw6WXZ3BZ8WlYMNeNUPcHBUQggh+hqPTpY9WcEHoJ5IfMtlYQIhhOjvPDpZ9lTZkFwGNjbS2GhxdShCCCFcSJJlB/wHWDBajWz9aI+rQxFCCOFCkiw7oCJtm0AfyD7u4kiEEEK4kiTLDpgmjQLgtEy1FEKIfk2SZQfSp8QBcLawe6NphRBCeAaPTpY9nTqSkBpOg6+FsDJfB0cmhBCiL/HoZNnTqSNeXgZqQhoIrVAOjkwIIURf4u3qANxd5JljqEqZaymEEP2ZR/csHaF2QBnaOsjVYQghhHAhSZadODAyDCtB7N9V5OpQhBBCuIhHJ8ueDvABsEZ4AZD/zm5HhSWEEKKP8ehk2dMBPgBJcUEAWLdKz1IIIforj06WjmC6yrYwAUdkrqUQQvRXkiw7kXxFGvW+Vsr8ZV9LIYToryRZdsLH14eaAY0Ey1ZdQgjRb3l0snTEAB+ACqMRTsvCBEII0V95dLJ0xAAfgPKogwyqb3JQVEIIIfoaj06WjuI7SBHUEMjRvSdcHYoQQggXkGTZBSrCtirgzs+OujgSIYQQruDRydJRzyxjJ6UAULJfpo8IIUR/5NHJ0lHPLNMvNQFQt092gRZCiP5Idh3pglGjI/nQbz/+pwe6OhQhhAvk5OREeHt7rwTS8PBORj9lBfY0NTUtzMzMLG+rgkcnS6XUHGBOQkJCj9rx8vbC7FvLsKPnHBOYEKJP8fb2XhkVFTUyPDy82mAwaFfHIxzLarWqioqK1LKyspXAdW3V8ejfkBx1GxYguKGYsEqzA6ISQvRBaeHh4TWSKD2TwWDQ4eHhZ7DdOWi7Ti/G4xBKqUCl1Gql1Aql1G29dd3aAVVYdRiWJktvXVII4T4Mkig9W/Pfb7s50S2SpVLqBaVUuVJqz0Xls5RS+5VSh5RSjzQXzwfWaa3vpp3ucqvzHTIaFiB/zBCsOogtn+T1uC0hhBB9i1skS+DvwKzWBUopL+AvwGwgFbhVKZUKxALHm6t12M1z5G3YgOH1AFQ/I3MthRCiv3GLZKm1/gyouqh4PHBIa31Ea20GXgGuB4qxJUzoxfhHXTaUCD4k8KNg6ovqe+uyQggBwLJly8JSUlJSU1JSUg0GQ2bLzwsXLozt/Oyva2pqYuTIkalXXHHFhRGQMTExo5OSklJTUlJS09LSvrLV0rp164JNJlNaXFxc2tKlS6PaarOjOu0ds7e8u+f0lFsky3bE8L8eJNiSZAywHrhBKbUceKu9k5VS9yilspVS2RUVFT0OZvQVkxkcsBatNUcePdLj9oQQwh4/+tGPKgsKCva+++67B6OioswFBQV7CwoK9q5cubK4O+09+eSTkQkJCecvLv/0008PFBQU7N2zZ8++lrKmpiYeeuihuHfeeefAgQMH8v/zn/8MysnJ8Wt9Xkd12jtmb3l32urOZ9MWd06WbW3zobXWdVrrO7XWi7XWa9o7WWv9vNY6S2udFR4e3uNgImLCSF+9jC2X1VC+tpwzW2Q1HyFE79uxY4d/cnLy15KcPQ4fPmx8//33Q+6+++7KrtTfuHFj4LBhwxpSU1PNfn5+ev78+VXr1q0L7Wqd9o7ZW96dtnryObXmzvMsi4Ghrd7HAqUuigWAqf8xUxAymMujfDj8w8OM3TQWpWTrLiH6k+9+l6F79uDQDW7T0jj3wgtfuZPWrtzcXP+RI0d2mCwzMzOT6+rqvC4uf+qpp47PnTu39nvf+97Qp59+uvjMmTNfqzNjxoxEpRR33nlnxZIlSyoBjh8/7hMTE3Nh7lxsbKx569atQa3P66hOe8fsLe9OWx19TvZw52S5HUhUSsUDJcAtwAJXBnQs+Dx7fV/kt08+w/6F+6l4tYKImyNcGZIQop/Jz8/3nzlzZk1HdXJycva3d2zt2rUhYWFhTVOnTj339ttvD2h9bNOmTQUmk6mxpKTEe/r06UmjRo2qnz179lmtvz5rRin1lcKO6rR3zN7y7rT1tcJucotkqZRaC1wOhCmlioGfaa1XKaW+D7wPeAEvaK3zXRgmAwYeJL5pK/7zAwn8cyCHHjpE6BWh+ET4uDIsIUQv6moP0FkKCgr8lyxZcrLl/aFDh4yPP/54dHBwsOXqq68+M3fu3NqOepZffPFF0IcffhgaExMT0tDQYKirqzNcf/318W+88cZRk8nUCBATE9N0zTXXnN68eXPg7Nmzz8bFxZlLSkoufNEVFxf7REdHN7Zuu6M67R2zt7w7bXXnM26LWyRLrfWt7ZS/A7zTy+G0a5h/EB9NW8kHm7Zy9epJ5EzIoeCOAkb/dzTKILdjhRDOZbFYKCoq8svIyLgwJD8vL8/faDTqhx9+uDwxMdEMHfcs586dW/uXv/ylBODtt98esGzZssg33njjaE1NjcFisTBw4EBrTU2N4ZNPPgl+7LHHSgGmTZtWV1hY6FdQUOBjMpka169fP2jNmjVfGenYUZ32jqWnp9fbU96dthz12btFsuwrRsaPZOyGn3Eg6UG+eW0QCc8mcPC+gxx/5jhxP45zdXhCCA+Xn5/vGxkZafb3979we3HevHk1JpPJvGjRorhVq1YVxcfHd6s3VVxc7D1v3rwEAIvFom644YZTN954Yw2A0Whk2bJlx2bNmpVksVhYsGBBZVZWVj3AtGnTElavXl1kMpka26vT0fn2lnenLUdQbd3n9TRZWVk6Ozu7x+3kHCwh6+VYbhmwnLU/vBetNXu/uZfK1yvJ+DyDkIk9X/xACOFaSqkcrXVW67Lc3NzC9PT0Lo0c7W2LFy+OsVgsqr6+3rBixYrjvr6+nv+l7iS5ublh6enppraOSc/SDhkjhjC5aALeAbbHBUopklYkUZtdy95b9pK1MwvjQKOLoxRC9CfLly8vcXUM/YE7z7N0O14GAycv/wHHh8ZcKDOGGkn9VyrmEjMF3ylAW+WXOiGE8DSSLO308IufsmTVq18pCx4fzIhlIzj15imO/faYiyITQgjhLJIs7RRWd5grjn5EU6P1K+Ux98cQcWsERx8/StWGi5e5FUII0ZdJsrRTSVI0H42fyOZNO75SrpQi6fkkAkYGsO/WfdQfk8XWhRDCU0iytFNZ1kSu/9Wv+Dh359eOeQd5k7Y+DWuDlfwb8rGcl42ihRDCE0iytNOVqaPYct99RGz7ss3jAUkBjPzHSGqza9l/1/42l2YSQgjRt0iytNOl0yaScWAfQUf2tlsn7Pow4n8dT/nacop+VdSL0QkhhHAGmWdpJ6OvD/8cP5Hj4aYO68U9Ese5vecofKKQwJGBhN/Q823ChBBCuIZHJ0ul1BxgTkJCQqd17fG3+QupGOzH0o6vTdKKJM4fOs++b+/Db7gfA8YO6OAMIYQQ7sqjb8Nqrd/SWt8TEuLYZegmFU2i8ZEJWDtZgMDLz4tRr43CONjInuv20HCiwaFxCCH6j2XLloWlpKSkpqSkpBoMhsyWnxcuXBhrb1uVlZVes2bNGh4fHz9q+PDhozZs2BAIsG7dumCTyZQWFxeXtnTp0qjW53R0rCt12jtmb3l3z+kpWRu2G+569l+8UHkX2d/ZT2ZiTKf1a3fVsnPKTgJHB5KxMQMvv6/tnCOEcBPuvjbs0aNHjVOmTEkpLS3N624b8+fPN1166aVnf/jDH1bW19ers2fPGkJDQy3x8fFp77///oHhw4c3pqenj3z55ZePZGZm1jc1NdHesZY2O6rT3rH09PR6e8q701brGDvT0dqwHt2zdJaEUCPTrTPZ8PHGLtUfkDGAkf8cSe3WWlkSTwjRIzt27PBPTk4+393zq6qqDFu3bh3w4IMPVgL4+fnpsLAwy8aNGwOHDRvWkJqaavbz89Pz58+vWrduXShAR8dadOd8e8u701Z3P6eLefQzS2cZGRvJ0vgHif30iy6fEz4vnOFPD+fIj49wOOYwCcsc+xxVCNF7xq8Yn3xx2fyR86seufSRitqGWsOMl2YkXnz89jG3Vz4w4YFTJ2pPeF//yvUjWh/bdve2dvefvFhubq7/yJEjO0yWHW3+HBER0TRo0KCmm266ybR3796AMWPG1K1YseL48ePHfWJiYswtdWNjY81bt24NAujoWIvunG9veXfa6uhzsocky274xqXjCHkyjgLDt+06b+iSoTQUN1D8+2J8Y3wZ+sOhTopQCOGp8vPz/WfOnFnTUZ2ONn/+7LPPAvbt2xfwxz/+8dj06dPr7rzzzqFPPPFEVHp6+tcSsFJKA23OF2851qKjOu0ds7e8O219rbCbJFl2g7+/D+fPD6KYQ3adp5Qi4dkEzCfMHP7RYXyifYi8JdJJUQohnKWjnuAA3wHWjo4PGTCkyZ6e5MUKCgr8lyxZcrLl/aFDh4yPP/54dHBwsOXqq68+M3fu3NqOepaXXHLJ+cjISPP06dPrAG6++ebqp556KmrOnDlnSkpKfFrqFhcX+0RHRzcCxMXFmds71qKjOu0ds7e8O23Z89l2RJJlN117Zg742v/xKYMi5aUUzCfNFHy7AJ8IHwZOH+iECIUQnsZisVBUVOSXkZFxYdBKXl6ev9Fo1A8//HB5YmKiGTruWQJERUWZc3NzfdPT0xs++OCD4OTk5Ppp06bVFRYW+hUUFPiYTKbG9evXD1qzZs0RgI6OtejO+enp6fX2lHenLUd99pIsu8lrWDIbMsJpamjC286k6eXnRdrraeycupM98/aQ8WkGAzJkDqYQomP5+fm+kZGRZn9//wu3F+fNm1djMpnMixYtilu1alVRfHx8p72pP//5z8duu+224WazWcXFxTWsXbu20Gg0smzZsmOzZs1KslgsLFiwoDIrK6seoKNj06ZNS1i9enWRyWRq7M759pZ3py1HkKkj3fThnSuZsfoejn68jxGXf+1Zf5fUF9ezc/JOrPVWxn4+loDkAIfGKISwn7tPHbnY4sWLYywWi6qvrzesWLHiuK+vr+d/qTtJR1NHPLpn6awVfADOJp3HoDVbP1jPiMsf7VYbfrF+pG9IZ+fUneTOzGXsF2PxG+bn4EiFEJ5s+fLlJa6OoT/w6HmWzlrBB2DEFVN5dOFCNp871aN2ApICSP8wHctZC7kzc2kok1V+hBDC3Xh0snSmUVmjWXnN1ZT49/wjDBoTxOh3R9NwooHdV+6mscphA7iEEEI4gCTLbvLy9uLtxXdw/5qVDmkvZGIIo98czbmD58i9MpfGakmYQgjhLiRZ9kDpoHASTp9xWHsDpw8k7bU06vbUkTtTEqYQQrgLSZY9kD/xmzz84GMUH6t2WJuDZw8m7fXmhCk9TCGEcAt2J0ul1AKl1CtKqTVKqZeVUrc6I7C+wCd1CluTU9mVXerQdgfPHmzrYeZJwhRCCHfQnZ7lNK31LVrr27TWC4BLHR2Uoyil5iilnj9zxnG3Slub/Y1ZFP4olcM1Dlur94LBVw8mbb0tYe6+arckTCGEcKHuJEtfpdQ1SqkxSqmrAX9HB+Uozpw6AmAaboVFl/Dvoyuc0v7ga2wJ8+zus7Ye5ilJmEII4QrdSZb3AQOBq5v//L5DI+pDBgT4MF3fxwCLufPK3dSSMOv21LFz2k4aTsg8TCGE6G12J0ut9Tmt9T+11k9prddorc85I7C+wuwfRm3A1xb3d6jB1wxmzDtjqC+sZ+fUndQXOWy5QyFEH7Fs2bKwlJSU1JSUlFSDwZDZ8vPChQtj7WknNzfXt+XclJSU1KCgoLG//OUvIwBiYmJGJyUlpaakpKSmpaWNbH3eunXrgk0mU1pcXFza0qVLo9pqu6M67R2zt7y75/RUl9eGVUrNBn4OhAK5wO+11lscGYyzOGNt2Bbjlv6YbK8/Yv7ZOYzezk2aZ7acIW92Hl5BXqRvSJe1ZIVwAndfG/bo0aPGKVOmpJSWlub1tK2mpiaioqLSv/zyy31JSUnmmJiY0dnZ2fuGDBnSdHG9+Pj4tPfff//A8OHDG9PT00e+/PLLRzIzM+u7Uqe9Y+np6fX2lHenrdYxdqajtWHt6Vn+FfghMBF4HnimP4+EbZEamQLeZr7cW+j0a4VMDCFjYwbWBis7p+6kdlet068phHAvO3bs8E9OTv7aRs3d8eabbwbHxcU1JCUldfgsaePGjYHDhg1rSE1NNfv5+en58+dXrVu3LrSrddo7Zm95d9pyxOcE9i2kflJrvan55w1Kqc3AVmCto4JxNGcupN7iitDh5Hs/R95He5g2ZoTTrtMiKD2IsZ+PJXdmLrum7SLt9TQGXiH7YQrRq8aP73yroVmzTvPLX568UP/22yt54IFTnDjhzfXXf/XLYlvXN4POzc31HzlyZIfJsqPNn+fOnXvht+y1a9cOuvHGG7+ywPWMGTMSlVLceeedFUuWLKkEOH78uE9MTMyFhBobG2veunXrV6YBdFSnvWP2lnenrY4+J3vYkywLlVJPAr/UWpuBRsCtuzZa67eAt7Kysu521jWmTBrFC+/m0nTwhLMu8TUByQGM/XIsu2ftZves3Yz8x0givhnRa9cXQrhOfn6+/8yZM2s6qtPZ5s8A9fX1asOGDSG///3vi1vKNm3aVGAymRpLSkq8p0+fnjRq1Kj62bNnn23rcZ1S6iuFHdVp75i95d1p62uF3WRPstTAfOBupdRBIA5Yo5RK1FofdFRAfU1iUjj/HvULdiWMB+7ttev6DfVj7Odj2XP9HvbeshfzSTOx99v1nF8I0V129AS/Vn/IkCa7z2+loKDAf8mSJSdb3h86dMj4+OOPRwcHB1uuvvrqM3Pnzq3tSs9y3bp1IampqeeGDh164fmkyWRqBIiJiWm65pprTm/evDlw9uzZZ+Pi4swlJSU+LfWKi4t9oqOjvzKXraM67R2zt7w7bdnz2Xak02SplJoEbNFa39r83g9IA9KbXyuVUsO11kMdFVRfczDSilkVoLVGKdVr1zUOMjLmgzHsW7CPQw8cwlxqJv7X8b0agxCi91gsFoqKivwyMjIuDFrJy8vzNxqN+uGHHy5PTEw0Q9d6lq+88sqgb37zm1Ut72tqagwWi4WBAwdaa2pqDJ988knwY489Vgowbdq0usLCQr+CggIfk8nUuH79+kFr1qw50rq9juq0dyw9Pb3envLutOWYT75rPcs7gL8opQ4A7wHvaa2zAecML3Wg3nhmCbA1LZOfPfAjlpWUsig2xqnXupiXvxej1o3iwPcOcOypYzSUNpC8IhmDjyz7K4Snyc/P942MjDT7+/tfuL04b968GpPJZF60aFHcqlWriuLj4zvtTdXW1hq++OKL4NWrVxe1lBUXF3vPmzcvAcBisagbbrjh1I033lgDYDQaWbZs2bFZs2YlWSwWFixYUJmVlVUPMG3atITVq1cXmUymxvbqdHS+veXdacsR7Jk6kgLMBr4BhACfYEuem7TWFkcF5AzOnDoC8OLcGymIHMwVt93DrMsynXadjmitKXqyiMKfFhJ6eSij/jMK4yCjS2IRoi9z96kjF1u8eHGMxWJR9fX1hhUrVhz39fV12HO6/qajqSNdfmaptS4ACoBnlVL+wBXATcDvgayOzvV0wZdM4rc/W8LL4cPRUy/hvNVKgJdz51xeTCmF6QkTfvF+7L9rPzsm7mD0f0cTkChzMYXwZMuXLy9xdQz9QZfv1SmlZiultiql9gOrgWqt9f0X/wbWH4255hsA1OVu4f6DB5m1ezeNVqtLYom6PYr0j9Jpqm5ix4QdVG903PZhQgjRX/VkUYLfufuiBM7edaRF4iWjOOUVSuLBQUwOCeGK0FC8XDjIJvTSUC7Zegk+kT7svnI3J17ovWktQgjhiexJlie11pu01tVa6w3Ynl0+5qS4HMLZu45coBTvxtzHFvMEFkRG8ov4eAxKYe3i82Bn8B/uz9jNYwm9PJT9d+3n8E8Ooy3yKEMIIbrDnmRZqJR6UinVMo/F7Rcl6E3r517Dz0YdxGq1JaSCujoysrPZVeu6j8gYamT0O6MZsmgIx58+Tt6cPBpPyzZfQghhL3uSZcuiBMeVUl8Ah4CNSqlEp0TmAL11GxbAN243vmnPkL3nEAADjUZ8DQYaXdi7BDAYDSQtTyJxeSLVG6rZMW4Hdfl1Lo1JCCH6mi4nS631rVrrVGAY8CDwCyAQ26IEx50TXs/02m1YYF6jpuYpOLz2FQAifXzYdskljAsObonF6TG0RylFzL0xZHySQVNtEzkTcqhYX+GyeIQQoq/pzn6W9VrrbK31Kq31A1rraf159Z4Wl1w7i6XTYWPj2QtlLSvp/Km4mAX79rn0GSZAyJQQsnKyCEwLJP+GfI48fkSeYwohRBf0uWVelFLDlVKrlFLrulC3127DJqTF81RmFJ8YTn7tWIPVSr3VitlF00la843xZeynY4m6K4pjvzpG3nV5NFbJc0whhOiIWyRLpdQLSqlypdSei8pnKaX2K6UOKaUeAdBaH9Fa39WVdnvzNixA/OlEYgsLvla+ZOhQ/jNqFH69vFBBewy+BpJXJJP410SqP6wm+5JsarZ1uImBEEL0a26RLIG/A7NaFyilvIC/YFtiLxW4VSmV2vuhdd3z2zJ46995WJu+2oNUSmFQiqrGRmbv3s2np0+7JsCLYopZHMPYL8YCsPPSnRT/udilz1aFEO1btmxZWEpKSmpKSkqqwWDIbPl54cKFdm83dNNNN5kGDRqUnpiYOKp1+bp164JNJlNaXFxc2tKlS6O6cqyjc3pyviNj6WqMHdJau8ULMAF7Wr2fBLzf6v2jwKOt3q/rpL17sC32nh0XF6d7w8ZvrdAa9LGNh9s8fsps1unbtulXTp7slXi6ynzKrHdfu1t/wid6z417dOPpRleHJITLANn6ou+TXbt2FWrbWA2Xv44cOZI7ZMiQhp608c477xR8/vnnexMSEs63lDU2NmbHxsbW5+fn7z5//nxOUlLSuezs7D0dHevonNYve893ZCxdjVFrnd3899xmTnGXnmVbYoDWo2yLgRil1GCl1HPAWKXUo+2drLV+XmudpbXOCg8Pd3asABgm2jZg3vzmv9o8PshoJDszk5sjIlpi7JW4OmMcZCTtjTSGPz2citcqyMnKoXaXTKEVoj3jc3KS/1RcPBigwWpV43Nykv9aUjIIoLapyTA+Jyd5RWnpQIBTjY1e43NykleXlYUCnGho8B6fk5P88smTIQDH6uvt2VeYHTt2+CcnJ5/vSfyzZ88+Gx4e3tS6bOPGjYHDhg1rSE1NNfv5+en58+dXrVu3LrSjYx2d05W2u9OuI9uyhzsny7bWi9Na61Na63u11iO01r/p9ag6MOJq2y3N2tyN7dbxNtg+8o+rq5m0YwdVje4xuEYZFHEPx5GxMQPLOQs7Ju6gZHmJ2yR0IYRNbm6u/8iRIztMlpmZmcktt2lbv15//fUB7Z1z/Phxn5iYGHPL+9jY2AubKbd3rKNzutJ2d9p1ZFv2sOs3ml5WDLSekhILlLooli6JNg2laIAXg4o73wTdADRpzXk3GCHbWuiloWTtyqLg2wUcvO8gVe9VkbwqGZ8wu/9tCeGxtmVmXvif3Ndg0K3fD/D2trZ+P9hotLR+P8TXt6n1+zg/v6/08DqTn5/vP3PmzA5H5HVl8+eLtfWLsVJKd3Sso3O60nZ32nVkW/Zw557ldiBRKRXfvMTeLcCbLo6pU/vDwkiqLOu03uUDB7ItM5MYX1+01jS5UdL0Cfdh9H9HM+LZEVS9V0X26GyqPqzq/EQhhNMVFBT4Z2RkXOhZHjp0yHjLLbcMu+eee2Jbeo7d6VnGxcV9pcdVXFzsEx0d3djRsY7O6Urb3WnXkW3Zwy2SpVJqLbAZSFZKFSul7tJaNwHfB94H9gGvaq3zXRlnV5TFJpJc1cCpk51vjWVQCq013z94kNv37cPiRrc8lUEx9MGhZG7LxHuQN7uv2s2hJYewNrhPUheiv7FYLBQVFfllZGTUt5Tl5eX5G41G/fDDD5fPnTu3Fmw9y4KCgr0Xv1qOt2XatGl1hYWFfgUFBT719fVq/fr1g2644YbTHR3r6JyutN2ddh3Zlj3c4jas1rrNrb601u8A7/RyOD3iNXoy3p9/Qf5bX3LZwms6ra+UwuTnR6CXl3v85nKRoPQgMrdncnjJYYqXFXP649OMfHkkgSmBrg5NiH4nPz/fNzIy0uzv73/hN+t58+bVmEwm86JFi+JWrVpVFB8f32mvac6cOfFbtmwZUF1d7R0ZGTnmkUceKX3ooYcqly1bdmzWrFlJFouFBQsWVGZlZdUDGI1G2jvWXvm0adMSVq9eXWQymRq7c74jY2mv3B6qPwzgyMrK0tnZ2b1yraINBxl2ZRKffecFLnvxTrvPP9PURLCX14Wl8txJ5ZuVFHy3AOs5KyOeGUH0vdEog/vFKURPKKVy9EWb2ufm5hamp6dXuiqmjixevDjGYrGo+vp6w4oVK477+vp6/pe6k+Tm5oalp6eb2jrmFj1LTzL08hHc7PMaCd6TuczOcyvNZibt3Mm3IiP5qcnkjPB6JOy6MMbljaPgOwUc/N5BKtZXkPJCCn5xfq4OTYh+a/ny5SWujqE/cMc7f32awdvAxtvf529+C+0+d7DRyHWDB3PlwIFOiMwxfIf4Mua9MSQ9l0TNlhq2p23nxAsnZIqJEMKjSbJ0gktUNbfu34DVYt9gGKUUyxISmNS8lu2es2c7OcM1lFJEL4pmXN44gi4JYv9d+8m7No+G0gZXhyaEEE4hydIJ5lfBnz88z+7Pt3W7jU1nzpCenc1LZZ1PQ3EV/3h/Mj7OIOGPCZz+5DTb07Zzcs1J6WUKITyOJEtnuPE7RP0IXi/t/niAicHB/Hb4cG7opaX6uksZFLEPxJK1K4uAlAD23b6P/Pn50ssUQngUSZZOcNX0yZwcAJuP5na7DS+lWBIXR6CXFw1WK2tOunePLSApgLGfj2X408Opeq+KbanbKH2+FG1135iFEKKrJFk6wbCoYO5751KufvuAQ9pbdeIEt+/bR06tey9urrxs68tm5WUx4JIBHFh0gF1X7OLc/nOuDk0IIXpEkqWT3HY8kutzNjukrXujo/koPZ2s4GCHtOdsAQkBpH+UTvKqZOp217E9fTtFvyrCapbVf4QQfZMkSydpSEknrvEQ1SWne9yWQSmmN08n2X32LIsPHKDRjdaSbYtSiiHfHcK4feMIuz6Mo48fJSczh5qtHa7/LIQQbsmjk6VSao5S6vkzZ870+rXLMiwY0Hy69t8ObffT06d5+9Qpyt1ka6/O+Eb5Mupfo0h7M42m003smLSDA98/QOPpvhG/EO5i2bJlYS0LohsMhsyWnxcuXBhrb1s33XSTadCgQemJiYmjWpfHxMSMTkpKSk1JSUlNS0sb2frYunXrgk0mU1pcXFza0qVLozorv5i953fUriPb6rL2doX2pFdmZqbubbs+3qQ16BVzb3Z426cbG7XWWlutVn2+qcnh7TtL45lGfeD+A/oTwyf6i4gv9InVJ7TVanV1WEJ8BZCtL/oO2bVrV6HWOtsdXkeOHMkdMmRIQ0/aeOeddwo+//zzvQkJCedbl0dHRzeUlpbuurh+Y2NjdmxsbH1+fv7u8+fP5yQlJZ3Lzs7e0155T8/vqF1HtnXxq/nvuc084tE9S1cafdlEjoYYmP7pe+Dgnm2It22Vwt8eO8bUXbuoabJrOzyX8Q72JvFPiWRmZ+I/3J+COwrYddkuzua55+ILQrijHTt2+CcnJ3e4+XNnZs+efTY8PLzLXxwbN24MHDZsWENqaqrZz89Pz58/v2rdunWh7ZX39PyO2nVkW/aQtWGdxOBl4KGZl/Dv17Lh1lvhrbfAy8uh10gLDOTg+fMEObhdZxswdgBjN42l7MUyDv/kMNljs4m9PxbTL0x4B8s/SeHeCr5bMLRuT12AI9sMTAs8l/JCyvGu1M3NzfUfOXJkh8kyMzMzua6u7mtfDE899dTxjrbpApgxY0aiUoo777yzYsmSJZUAx48f94mJiTG31ImNjTVv3bo1qL3yi9u09/yO2nVkW/aQbyYnqki9ku+f38Hf3nkXliyBZ591aPvXhoVxbViY7VpmM4X19YzrIyNmlUEx5K4hhM0L48jSIxT/sZjyf5UzYtkIIm6JcMtdV4RwB/n5+f4zZ87scKRcTk7O/u60vWnTpgKTydRYUlLiPX369KRRo0bVz549+6xuY463Ukq3V35xmb3nd9SuI9uyhyRLJ7o+6QZ+8u5QfvrNPGJWrIAHH4Rhw5xyrR8cOsT7VVUUTpzIAO++89dqHGQk+blkhtw1hIP3HWTfgn2UPldKwh8SGDC23U3dhXCZrvYAnaWgoMB/yZIlJ1veHzp0yPj4449HBwcHW66++uozc+fOre1uz9JkMjUCxMTENF1zzTWnN2/eHDh79uyzcXFx5pKSEp+WesXFxT7R0dGN7ZVf3K6953fUriPbsoc8s3SiGydnQvZi3rniz5CT47RECfCHhAT+lZrapxJla8HjgrlkyyUkPZfEub3nyMnMoWBhAQ1lsmyeEC0sFgtFRUV+GRkZFzYvzsvL8zcajfrhhx8ub0mEOTk5+wsKCvZe/OooUdbU1Biqq6sNLT9/8sknwWPGjDkPMG3atLrCwkK/goICn/r6erV+/fpBN9xww+n2yi9u297zO2rXkW3Zo29+s/YR8fEQFHeYd/OruPvecbbC556DzEwYN86h14rw8WHmoEEAvF1ZyaaaGn4VH4+hD93OVF623UzCbw6n6MkiSv5UQsW/Koh7LI7YB2Px8utbz2aFcLT8/HzfyMhIs7+//4XbiPPmzasxmUzmRYsWxa1ataooPj6+017TnDlz4rds2TKgurraOzIycswjjzxS+o1vfKNm3rx5CQAWi0XdcMMNp2688cYaAKPRyLJly47NmjUryWKxsGDBgsqsrKx6oN3yadOmJaxevbrIZDI1duf89sod2ZY9VFv3cz1NVlaWzs7Odsm1Bz04kwbDaep+nw1nz8KYMXDFFbBqldOuueTQIT45fZovxo7Fv48N/mnt3MFzHH74MKfeOIVfvB8jfjeCsPlh8jxTOJVSKkdrndW6LDc3tzA9Pb37OyM40eLFi2MsFouqr683rFix4rivr6/nf6k7SW5ublh6erqprWPSs3SyhIBMtnv9gbp6M4FBQfD55xDVvTmxXfW7ESOos1jw9/LCbLVSbjYT6+fn1Gs6Q0BiAKNfH031R9UcevAQ+TfmE3JZCAnPJjDgEnmeKQTA8uXLS1wdQ3/g0c8sXbmCT4sJcZeAt5n/bs+3FcTE2KaQnDgBd98N5xy/yLhSiqDmZ5dLDh8mKyeHqj6y4k9bBs4YSObOzK88z9x7+17OF/ZoqpkQQnSZRydLrfVbWut7QkJCXBbDNZdkAvDurpyvHsjOhhdegBtvBLO5jTMd477oaJYOG8Ygo9Fp1+gNBm8D0YuiGX9wPHGPxFH5n0q2JW/j0EOHMFc67/MTQgjw8GTpDmZmDoeGYLYXX5Qs58yxDfZ591244w6wWJxy/ZTAQB6ItS0dubeujp8dPYqlDz+nNoYaGf6b4Yw/OJ7Ib0VS/Kdito7YStGvi7Ccc85nKIQQkiydzNvLwNh9b+P15eNfP3j33fDb38Irr8D994OTk9i6igr+VlpKZR++JdvCL9aPlJUpjNs9jtDLQzn62FG2JmyldEUp1ib33pFFCNH3SLLsBVclT2Xf1hga2poy+OMf217Ll8NPf+rUOH5qMrErK4tIHx+01hw53/ef+QWOCmT0G6PJ+CwDP5MfB+45QPbobCpeq2hzpQ8huslqtVplGLYHa/77bfc3bUmWvSD5kgoaM//A218earvCU0/BwoXw5JPwzDNOjSXK1xeAF8rKSN22jZ21HS4T2WeETg1l7KaxjHptFFpr8ufnkzMuh1PvnJKkKRxhT0VFRYgkTM9ktVpVRUVFCLCnvToydaQXpIyug1kP8a8cH264IuHrFZSyPb+sqYGHHwZvb9vSeE503eDBFA8bRnqQ3esJuy2lFOFzwxl87WBO/vMkRb8oIu+aPIInBWP6pYmBMwbKHE3RLU1NTQvLyspWlpWVpSGdDE9kBfY0NTUtbK+CLErQC6xWjfHRKExNszm87O/tV2xstO1QEhAAq1fbkmgvqGpsZPGBAzwzYgRD++B8zPZYzVbKXiyj6MkiGoobCJkWQvz/iyd0aqirQxNurK1FCYTw6N+Q3GGeJYDBoAhvHE+x3tpxRaMR1q6FF1+0Jcpeeqa4p66Oj6qrOeHEKSyuYPD533SThD8lcH7/eXZdtovcq3Kp2drhpg1CCPEVHp0s3WGeZYu0gRMwhxRQdLKTxG00/m/RgvR0WLnS6bFdFhpK4cSJjG/e3mtDVRVmq+eMKPXy8yL2/lgmHJ7A8N8N5+zOs+yYuIPd1+6WpCmE6BKPTpbuZEbKeLB68eaXe7t2wqBBkJUFaWnODaxZy4o/R86fZ9bu3TxZVNQr1+1NXgFexC2JY8KRCcT/Kp6azTXsmLiD3KtyOf3ZaVeHJ4RwY/LMspecrDQTFd3Er38RwKOPdqOBnBzbbiW94O3KSi4LDSXY25sGqxVfg2f+TtVU20Tpc6Ucf+Y4jeWNhEwNYdgTwxg4UwYC9WfyzFK0xTO/Bd1QZJgPiaYAtm3rxslvvmnrZf72tw6Pqy3XhoUR7O2NRWtm7d7Njw8f7pXr9jbvAd7EPRzHxMKJtmeaR86z+6rd7Ji4g8q3KmXKiRDiAo9Olu4ywKdF1PR1vBdyg/1fwrNn20bJPvIIPPGE01f6aWHVmqwBAxgZENAr13MVL3/bM82JhyeS9LckGssb2XPdHrLHZlO+rhxtlaQpRH/n0cnSnQb4AEQOL6c+fj3bDxy370SjEf7xD7jrLtvCBT/6Ua8kTKPBwO9GjODOIUMAeKuykt8UFfXptWU7YvA1EH1PNOMPjCdldQrW81b23rSXbanbKF1ZirXBcwY9CSHs49HJ0t1cnT4egH9t6mQKSVu8vOD55+GBB+DZZ+Hee522+Hp73q2q4tWKCpo8NFm2MBgNRH07ivF7x5P6SipeAV4cuPsAW0xbKHqqiMbTfX9tXSGEfSRZ9qIbp46BJl8+PdyNZAlgMMAf/gBLl9oS5y23QH29Q2PsyF8SE9mYkYGvwUCD1cr6Cs9ef1V5KSJujiAzJ5P0DekEjgnk6KNH2TJ0C4d+dIj647332QshXEuSZS8aEODDgLNj2Vf/ERZrN3uFSsGvfgXLlsG6dbbnmb30TFYpRUjzFJMVpaXckJ/PjrNne+XarqSUYuCMgaS/n07mzkwGXzeY4j8Ws3X4VvbdsY+zeZ7/GQjR33l0snS3AT4AU/zv4fzJGE6f7WGv5Ic/hDVrID8fSkocE5wdFsfE8N6YMWQOGADgETuYdMWAjAGkrkll4uGJRH8vmop1FWSPyWb31bup/rjao3vaQvRnMs+yl73+Osybb2XTFwYmTLTgZfDqWYNnz0JQkG3AT3k5REY6JE57FNXXk7ptGz81mfhJXFyvX9+VGk81UrK8hJI/ldBY0UjgmEBifxBLxIIIvPx6+HcrXELmWYq2eHTP0h1NmgRoAx9uqmTSqkms37e+Zw227Bryxz/CqFFw5EiPY7RXjI8PPzeZuC0iAoAGq7Xf9LCMg42YHjcxsWgiyauSQcP+u/azJW4LR396lIYTbW1iKoToazw6WbrjbdjISBgxAnK2+uNt8ObW/9zKR0c+6nnD114L3/0umEw9b8tO3gYDD8fFEdu8Y8l3Cwq4de/efpMwwTZXc8h3h5CVm0X6R+kETwym6Mkitgzbwr5v76M2xzP2DRWiv/LoZOlu8yxbTJ4M274I5K1b3yZpcBLXv3I920q6s7RPKwkJ8PTTthGzR4/appe4IFlprUkPCmLsgAH9csk4pRQDpw9k9JujGX9gPNGLo6l8rZKcrBx2Tt1JxX8qsDbJfE0h+hqPTpbuavJkOHkSzpQN4oPbPyAiMILZa2azt6KLi6x3ZsUK2wCgu+6CXt52SynFj+PiLjy73HLmDDN37aK4F6e4uIuAhAAS/5jIpOJJjPj9CBpKGsi/MZ+tCVs59ttjmCs8a0s0ITyZJEsXmDzZ9ueXX8KQAUP48FsfMjpiNEE+QY65wK9+BT/9qW1fzKuugqoqx7TbDSVmM+WNjRemnPRH3iHeDH1oKBMOTmDUa6Pwj/fnyCNH2By7mX3f2seZzWf61S1rIfoijx4Nq5SaA8xJSEi4++DBg64O5wKLBQYOhNtvh7/+9aJjVgunzp8iIjCi5xdas8b2HHPYMPjvfyExsedtdoNVawxKYdWahfv3c9eQIUxxs1vjva1ubx2ly0spW12GpdZCUEYQ0fdFE7kgEq9AGUXrSjIaVrTFo3uW7vrM0ssLJk609Swvdv+79zN51WSKa4p7fqHbboOPPoLqahg/Ht5/v+dtdoOh+dllcUMDH1dXc/DcOZfE4U4CUwNJ/HMik0onkfRcEtqqOXDPAb6M/pKDDxykbl+dq0MUQrTS55KlUmquUmqFUuoNpdRVro6nuyZPhrw8qKn5avm3079NeV0501dPp7S2tOcXuvRS2LoVhg6Fq6+2bfPlorsJcX5+7B0/njuiogB4vaKCF06cwOrBdzc64x3kTfSiaLJ2ZTH2i7EMvnYwpc+Vsj11O7um76J8XTnWRhkQJISruUWyVEq9oJQqV0rtuah8llJqv1LqkFLqEQCt9eta67uB7wA3d9Ku200daTFlClitsHnzV8snxk7kvdvf48TZE0xfPZ2ys2U9v9jw4bYL3XijbZuvt9/ueZvdFODldWGU7Mvl5fy1pIT+myr/RylFyJQQUtekMql4EvG/juf84fPsvWkvm2M3c/gnhzl3UHrkQriKWzyzVEpdBpwFXtJapzWXeQEHgCuBYmA7cKvWem/z8WXAGq31js7ad6cVfFqcPQuhofCTn9jG41zs86LPmbVmFqnhqWxduBWDcsDvNVrbEuW119rWmG1stG3/5SJaayobGwn38eGcxcIvCwtZMnQoYT4+LovJnWiL5tS7pzix4gSn/nsKLBB6eShD7h5C2PwwWSHISeSZpWiLWyRLAKWUCXi7VbKcBPxca/2N5vePNld9qvn1odZ6Qwft3QPcAxAXF5dZVFTkxOi7Z+JE8PaGL75o+/jGwo0AXG663PEX37/ftgj7Sy/ZbtW62HunTjFnzx42pKczLTTU1eG4nYbSBsr+XsaJlSeoP1qP90BvIr8VyZCFQwga7aBR1AKQZCna5ha3YdsRA7TeJbm4uex+YCZwo1Lq3vZO1lo/r7XO0lpnhYeHOzfSbrrsMti2Ddpbg/xy0+UXEuVLuS9RdNqBCd9ohPh4cJO1XGcNHszRCRMuJMqXysrY5Ia3z13FN9qXYUuHMeHQBNI3pDPoG4Mofa6U7DHZ5EzMoXRlKU1nm1wdphAey52TZVvLv2it9Z+01pla63u11s/1elQONG2a7U7oli0d16s6X8VD7z/E1BencvCUg6bADB9uGykbF2d7ePqLX0BFhWPa7qaW5fKarFZ+XVTEs8ePd3JG/6MMtu3CUtemMqlkEiOeHYGl1sKBuw+wechm9t+9nzObZN6mEI7mzsmyGBja6n0s4IDhoe5jyhTbo8PPPuu43iD/QXz07Y8433Sey/5+Gfnl+Y4NZPdu+M1vYOzY9u8J9yJvg4HszEz+kpQEQGlDA0uPHKGmSXpOrfmE+TD0waGM2zOOsV+OJfymcE6+fJKdl+5kW/I2Cp8spP5Y/1s5SQhncOdkuR1IVErFK6V8gFuAN10ck0OFhkJGBnz6aed1M6Iy+Ow7n6FQTPv7NHac6HRcU9dlZNi6t/7+cPnl8Lvf2XqbLhTk7U1k80Cfd6uq+P3x41Q0Nro0JnellCJkUggpL6Qw+eRkUv6egm+ML4VPFLLFtIVdM3ZR9o8yLHXd3HBcCOEeyVIptRbYDCQrpYqVUndprZuA7wPvA/uAV7XWDu5Sud60abZZHV1ZwnVk+Eg+v/NzgnyC+PJ4Gysa9ERGBmRnw7x58OMfw/XXu3SZvNbuGjKEoxMnMsLfH4D/V1jIB24Sm7vxDvIm6o4oMj7JYMLRCZh+bqK+sJ6CbxfwZdSXFHy3gNOfnkZb5TatEPZwm9GwzuSOU0davPYazJ9vu/s5ZUrXzqltqGWA7wAATp07xeCAwY4LSGv4v/+DH/0IoqLgH/+wZXQ3cc5iISM7m+vDwvjdiBGuDqdP0Fpz5oszlK0uo+LVCiy1FvxMfkTeEUnUt6PwH+7v6hDdioyGFW1xi55lfzZ1qu3PrtyKbdGSKPNO5jHiTyNYvn254wJSCu6/37YWn58fXHEFPPaYbSSSGwjw8mLPuHH8vHnfzh21tXxr3z5O9vLuKn2JUorQqaGkrExhctlkRv5zJP6J/hT9soitI7ayY8oOSv5aIrugCNEBSZYuFhYGaWmdD/Jpy4hBI5g6bCr3vXMfT3z8hGNHQGZlwY4dcOedsHGjLYm6CR+DgUAv24T8vLo6Np4+jZ/B9k+5P9wp6QmvAC8ib4sk/YN0JhZNJP7X8TSdaeLg9w6yOXozu6/Zzck1J2UaihAXkduwbuB737OtDVBdbVukwB5N1iYWvbWIF3a9wMKxC1l+7XK8DQ7eDuvcOQgIsE0t+eADWLDArZJng9WKr8GA1po5eXnMHDiQB4cO7fxEccHZ3Wc5+fJJyl8up+F4A4YAA2HXhxF5WyQDrxqIwdh/fq+W27CiLf3n/wA3Nm2abfm7nBz7z/U2eLPyupU8PvVxVu5cyd+y/+b4AAMCbH/++c+2DaXdbDUk3+Ze5XmrlQAvL3xa9TIbXDyqt68IGhPEiKdGMLFwIhmfZRD5rUiq3q8i79o8vhzyJQfuO2CbvykDg0Q/JT1LN1BRARERtjVily7tfjuvF7zONYnXYPQyorW+sGC5w1gstluz48bZ3u/aZRtF62Za/ttfr6jgB4cO8VF6OgktCV90mdVsper9KspfLqfyjUqs5634DvMlckEkETdHEDgm0PH/xtyA9CxFW6Rn6QbCw2HMGNuCOj0xN2UuRi8j5XXlTH1xKjml3eiqdsTL63+J8r//tS1icM89UFvr2Ov0UMsXeLiPD5OCgzE1rwx0+Px5GqWn2WUGHwNhc8JIXZtqm7/5jxQCRwZy7OljZGdksy1lG0efOMrZvLPyrFh4PEmWbmLmTNi0qf11Yu1Rdb6K4zXHuezvl/F6wes9b7AtM2bY5mOuXAmjR8MnnzjnOj0wJSSEV0aNwttgoMlqZdbu3dy8d6+rw+qTvAd4E3V7FGPeHcPkE5NJ+lsSvkN9Kfp1Edljstmeup2jPztKXb5sWi08kyRLNzFjBjQ02BJmT6WEpbB14VbSItKY/6/5PPPlM47/zd/Pz7aR9BdfgI8PTJ8ODzwAde75ZemlFMtGjOCBmBgA6i0WXjl5Eov0iOzmE+5D9D3RZGzIYPKJySQuT8RniA9FTxaxPW0721K3cfTnR6nb657/FoToDnlm6SbOnoWBA2HJEtsyrY5wvvE833njO7ya/yrPfuNZHpz4oGMavti5c/Doo/CnP0FCArz4olts+9WRl8rKuKOggM8yMpgqW4I5RENZA5XrK6n4dwWnPz0NGgJGBRDxzQjCvxlOYEqgq0PsEnlmKdoiydKNTJ1q611u2+a4Nq3ayu83/57vZHyHsIAwxzXclo0bbfMyCwth8WJ4+mkIcs+9Fq1as/H0aaYPHAjAitJSjEpxR1SURw5a6W0NJ2yJs/zVcs58fgY0BKYFEjYvjLD5YQSlB7nt5yzJUrRFkqUb+fnP4Ze/hFOnbL1MRzNbzHzrtW/xgwk/YPLQyY6/ANi6yE88Ae++Czt32hZn7wOuzM3Fz2DgrdGjAVsyNbjpl3lf01DaQMW6CirWV9gSpxX8TH6EzQ8jfH44wZOCUQb3+awlWYq2ePQzS6XUHKXU82f6yCbCM2bYlmbduNE57ZedLSOnNIfL/345z+c875yLBAXBs8/+L1HW1dmWzysrc871HOSDMWP458iRAFSYzSRs3cq7p065OCrP4BvtS+wDsYzdONY2OGhFEgGpAZT8Xwk7L93Jl9Ffsv/e/VS9X4XVLKOVhXvy6GSptX5La31PSEiIq0PpkgkTIDAQNmxwTvtxIXFsv3s70+Ons+jtRdz79r2YLU5aD7SlR7lli23E7IEDzrmOgyilCGlePulMUxNJ/v7EN085KTebOeUma+P2dT4RPkQvjGbMf8cwpWIKI9eOJPSyUE7+8yS7Z+1mU8Qm9t6+l4r1FbKlmHArchvWzVx9NRw5AgUFzruGxWrh8Y8f56lNT3H7mNv5x7x/OO9iAOXltlUXAJ5/3rZkUXKyc6/pQIsPHODV8nKOT5pEQPOatMKxLOctVG+opvK1SirfqKSpqgmDv4FB3xhE2LwwBl8zGONgY6/EIrdhRVskWbqZZctsI2KPHIH4eOde69X8V0kenEx6VLpzVvy5WHW1bbRsba3tP/Lxx/+3lJ4by6+rY1tNDXcOGQLAn4uLuTw0lNFuOnipr7M2WTnz2RkqX6uk4rUKzCVmMEDIlBAGzxlM2HVhBCQ779+NJEvRFkmWbubwYdsc/8RE+PhjGOzArSo7cu/b9xIRGMHPpv0ML4MTe08nT9oWM3jpJYiLgz/+0bbRdB8ZTHO6sZFhW7bwvZgYfj18OEDv/KLRT2mrpjanllNvnaLyzUrqcm1zN/0T/S8kzuApwRi8HfdESZKlaItHP7Psi0aMgDffhP374corbZ0xZ7NYLZgtZv7fZ/+PK/9xJWVnnTgYJzISVq+27UkWHAzz5sGcObaudB8QajRSOHEiDzfvarLlzBkysrPZ56aLMfR1yqAIHhdM/C/jGbdrHBOLJpL4f4n4Dfej5P9K2HX5Lr6M+JK9t++l/F/lNJ2RrcWEc0jP0k29956tw5WeDh9+CL0xRunFnS/yvXe+R4hfCGtvWMvlpsude8HGRttCBj//OTQ12RY2+PGPbasD9RGfVFfzxNGjvDdmDEHe3hw8d44IH58Lg4WE8zTVNlH9QTWVb1Vy6u1TNJ1qQnkrQqaFEDYnjMFzBuM/3P6pS9KzFG2RZOnG3noL5s+3rV3+wQe9M78/72QeN/37Jk7WnaTwB4WE+PVCli4pgR/+EF591fYs83e/c/41nWTqzp3UNDWxKytLbs32Im3R1GypsSXON09xbt85wLaC0Og3R9uVNCVZirZIsnRz69fbbsuuXGn/xtDdVdtQy+6Tu5kSN6V3Ltjiww9tXemWkbN90I7aWsrNZmb11sNm0aZzh85x6q1TVG+oJu21NAw+XX/iJMlStEWSpRBCtCLJUrTFowf49LUVfIQQQrgnj06WfW0FHyGEEO7Jo5OlEEII4QiSLIUQQohOSLIUQgghOiHJUgghhOiEJEshhBCiEx6dLGXqiBBCCEfoF4sSKKUqgCJXx9GGMKDS1UF0wJ3jc+fYwL3jc+fYwPXxDdNah7vw+sIN9Ytk6a6UUtnuvFKIO8fnzrGBe8fnzrGB+8cn+iePvg0rhBBCOIIkSyGEEKITkixd63lXB9AJd47PnWMD947PnWMD949P9EPyzFIIIYTohPQshRBCiE5IshRCCCE6IcnSyZRSXkqpnUqpt5vfD1JKfaiUOtj858BWdR9VSh1SSu1XSn2jF2ILVUqtU0oVKKX2KaUmuUt8SqmHlFL5Sqk9Sqm1Sik/V8amlHpBKVWulNrTqszueJRSmUqpvOZjf1JKKSfG97vmv9vdSqnXlFKhroivrdhaHVuilNJKqTBXxCZEl2mt5eXEF/BD4GXg7eb3TwOPNP/8CPDb5p9TgVzAF4gHDgNeTo5tNbCw+WcfINQd4gNigKOAf/P7V4HvuDI24DLgEmBPqzK74wG2AZMABbwLzHZifFcB3s0//9ZV8bUVW3P5UOB9bAuGhLnqs5OXvLrykp6lEymlYoFrgJWtiq/HlqRo/nNuq/JXtNYNWuujwCFgvBNjC8b2JbYKQGtt1lqfdpf4AG/AXynlDQQApa6MTWv9GVB1UbFd8SilhgDBWuvNWmsNvNTqHIfHp7X+QGvd1Px2CxDrivja+ewAngV+DLQeZdjrn50QXSHJ0rn+gO3LwNqqLFJrfQKg+c+I5vIY4HiresXNZc4yHKgAXmy+TbxSKRXoDvFprUuAZ4BjwAngjNb6A3eI7SL2xhPT/HNvxwnwXWy9MXCD+JRS1wElWuvciw65PDYh2iLJ0kmUUtcC5VrrnK6e0kaZM+f1eGO7NbZcaz0WqMN2K7E9vRZf87O/67HdhosGApVSt7tDbF3UXjwuiVMp9RjQBKxpKWonjl6JTykVADwG/LStw+3E4G5/x6KfkWTpPFOA65RShcArwHSl1D+Bk823lGj+s7y5fjG2ZzgtYrHdenSWYqBYa721+f06bMnTHeKbCRzVWldorRuB9cBkN4mtNXvjKeZ/t0J7JU6l1B3AtcBtzbcv3SG+Edh+Ecpt/v8jFtihlIpyg9iEaJMkSyfRWj+qtY7VWpuAW4CPtda3A28CdzRXuwN4o/nnN4FblFK+Sql4IBHbgAZnxVcGHFdKJTcXzQD2ukl8x4CJSqmA5hGPM4B9bhJba3bF03yrtlYpNbH5v+vbrc5xOKXULOAnwHVa63MXxe2y+LTWeVrrCK21qfn/j2LgkuZ/k27x2QnxNa4eYdQfXsDl/G807GDgI+Bg85+DWtV7DNvov/30wkg/IAPIBnYDrwMD3SU+4BdAAbAH+Ae20ZEuiw1Yi+35aSO2L/e7uhMPkNX833QY+D+aV9FyUnyHsD3/29X8es4V8bUV20XHC2keDeuKz05e8urKS5a7E0IIIToht2GFEEKITkiyFEIIITohyVIIIYTohCRLIYQQohOSLIUQQohOeLs6ACF6k1KqZboHQBRgwbbsH8B4rbXZJYEJIdyaTB0R/ZZS6ufAWa31M66ORQjh3uQ2rBBCCNEJSZZCCCFEJyRZCiGEEJ2QZCmEEEJ0QpKlEEII0QlJlkIIIUQnZOqIEEII0QnpWQohhBCdkGQphBBCdEKSpRBCCNEJSZZCCCFEJyRZCiGEEJ2QZCmEEEJ0QpKlEEII0Yn/DxWIyNdKF+CfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<ClimateGraphicsMPL.plotObj at 0x7fb279548d60>"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psMax = 10.e5 #A small 10 bar ocean\n",
    "cList = []\n",
    "for Ts in [450.,600.,750.,1000.,1500.]:\n",
    "    pp0List = numpy.linspace(min(pp0Fun(Ts),psMax/p0),.1,1000)\n",
    "    TList = [T(pp0,Ts) for pp0 in pp0List]\n",
    "    c = Curve()\n",
    "    c.addCurve(pp0List)\n",
    "    c.addCurve(TList,'$T_s$=%f'%Ts)\n",
    "    c.switchXY = c.reverseY =c.YlogAxis = True\n",
    "    c.Ylabel = 'T'\n",
    "    c.Xlabel = '$p/p_0$'\n",
    "    cList.append(c)\n",
    "plot(*cList)\n",
    "#Equivalent to: plot(cList[0],cList[1],cList[2],cList[3],cList[4])\n",
    "#In python argument lists, \"*\" is the unpacking operator\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "88f0aa88",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Integrand function, for OLR computation\n",
    "#Changed integrand to work with quadrature\n",
    "#delTau is optical depth relative to TOA (positive)\n",
    "def f(delTau,Ts):\n",
    "    delTau = min(delTau,100.)\n",
    "    pp0 = delTau/tau0\n",
    "    return phys.sigma*T(pp0,Ts)**4.*exp(-delTau)\n",
    "m = romberg(f) #Quadrature integrator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "8da8e57b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Function to compute OLR. \n",
    "def OLR(Ts):#Ts is surface temperature\n",
    "    ps = min(p0*pp0Fun(Ts),psMax)\n",
    "    tauInf = kappa*ps/g\n",
    "    tauInf = min(tauInf,10.) #Avoids underflow in optically thick case\n",
    "    #Carry out quadrature. Note integral is over delTau.\n",
    "    OLR1 = m([0.,tauInf],Ts,tolerance = 1e-3) #Set tolerance to speed things up\n",
    "    #Now we have to add in the contribution from the surface\n",
    "    return OLR1 + phys.sigma*T(ps/p0,Ts)**4.*exp(-tauInf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "1fa5d317",
   "metadata": {},
   "outputs": [],
   "source": [
    "#--------Main script starts here-------------------------------\n",
    "#Constants for Clausius-Clapeyron. These are globals\n",
    "#Note that (p0,T0) can be any point on the phase boundary\n",
    "T0 = 300. #Reference temperature\n",
    "p0 = phys.satvpg(T0) #Water vapor saturation pressure\n",
    "RTL = phys.water.R*T0/phys.water.L_vaporization\n",
    "\n",
    "#Absorption constants. Globals (used in OLR function)\n",
    "g = 10.\n",
    "kappa = .1 #Very roughly appropriate for water vapor\n",
    "\n",
    "tau0 = kappa*p0/g\n",
    "#Note: If you want to set p0 so that kappa p0/g = 1, as in the text,\n",
    "#and then compute the corresponding T(p0), you can do the following instead\n",
    "#The answer should be the same as using the fixed reference pressure above.\n",
    "#p0 = g/kappa\n",
    "#The next line estimates the temperature at the pressure p0, using\n",
    "#                         the simplified form of Clausius-Clapeyron\n",
    "#T0 = 300./(1.-(phys.water.R*300./phys.water.L_vaporization)*log(p0/phys.satvpg(300.)))\n",
    "#tau0 = 1.\n",
    "#RTL= phys.water.R*T0/phys.water.L_vaporization\n",
    "\n",
    "#Mass path of available ocean, expressed as psMax.  Ocean mass path is psMax/g\n",
    "psMax = 100e5 #Earth ocean would be about 200e5 Pa \n",
    "Rcp = 2./7. #Needed for dry adiabat after ocean is used up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "1a219ec8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWkAAAD4CAYAAAAuNhccAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAe70lEQVR4nO3de5RU1Zn38e/DnSgqhIa0gIIGM7aKKB3U+MYrS9HEYJYLhYlIRg0ZB2IcYxJMMonGYcVkvOTVN+qgosxogkTjgMYbokaNjtgoIheRVhGRBoqLiATQ7n7eP/bpUEDTXd1V1efU6d9nrVrn1O5T1U8Xq3+c3mefvc3dERGRZOoQdwEiIrJ3CmkRkQRTSIuIJJhCWkQkwRTSIiIJ1inuAgB69+7tAwcOjLsMESkB8+fPX+/uZXHX0VYSEdIDBw6kqqoq7jJEpASY2ftx19CW1N0hIpJgCmkRkQRTSIuIJJhCWkQkwRTSIiIJppAWEUkwhbSISIIppEUkNvfeCw8+CJoxee8U0iISi5UrYdIkuOuuuCtJNoW0iLQ5d5g4MWzvuAPM4q4ouRJxW7iItC9//CM8+ijceCNo2p6m6UxaRNrUpk1w+eUwbFjYStNyDmkz62hmr5vZo9HzXmY2x8yWR9ueWcdebWbVZrbMzM4sRuEiUpquugrWrw990Z30t3yzWnIm/X1gadbzycBcdx8MzI2eY2YVwBjgCGAkcJuZdSxMuSJSyp58EqZNC0E9dGjc1ZSGnELazPoDXwOyr8OOAqZH+9OBc7PaZ7j7Dnd/D6gGhhekWhEpWZs3w6WXQkUFXHNN3NWUjlzPpH8L/Aioz2rr6+41ANG2T9TeD/gg67hVUdsuzGyCmVWZWVUmk2lp3SJSYn7wA1i9Gu65B7p1i7ua0tFsSJvZ14F17j4/x/dsbDDNHkPV3X2qu1e6e2VZWbtZZEGkXXrySbj7bvjhD2G4/q5ukVy67U8EvmFmZwPdgP3M7D5grZmVu3uNmZUD66LjVwEDsl7fH1hdyKJFpHQ0dHMcfri6OVqj2TNpd7/a3fu7+0DCBcFn3P1CYDYwPjpsPDAr2p8NjDGzrmY2CBgMzCt45SJSEq66KnRz3HuvujlaI58BMNcDM83sEmAlMBrA3Reb2UxgCVALTHT3urwrFZGS8+STYajdj3+sbo7WMk/AzCaVlZWuhWhF0mXDBjjqKDjgAHjttcKdRZvZfHevLMy7JZ+GkotIwbnDd74Tblp57DF1c+RDIS0iBTdtGjz8MPzHf+imlXxp7g4RKajly+H734fTToMrr4y7mtKnkBaRgvnsM7jwQujSBaZPhw5KmLypu0NECua662DePJg5E/r3j7uadND/cyJSEC++CFOmwPjxMHp03NWkh0JaRPK2fj2MHRsm8L/llrirSRd1d4hIXurr4aKLYN06ePll2G+/uCtKF4W0iOTlN7+Bxx+H3/0Ojj027mrSR90dItJqL7wAP/sZXHABXHZZ3NWkk0JaRFolk4ExY2DQIJg6VSt+F4u6O0SkxerrYdy4MD/Hn/+sfuhiUkiLSItdd12Y4e6OO3Tbd7Gpu0NEWmT27DB5/0UXwYQJcVeTfgppEcnZW2+F276HDQtn0eqHLr5c1jjsZmbzzOwNM1tsZtdG7deY2YdmtiB6nJ31mqvNrNrMlpnZmcX8AUSkbWzeDOeeG6Ydffhh6N497orah1z6pHcAp7n7J2bWGXjRzB6Pvnazu9+QfbCZVRCW2ToCOBB42swO0+osIqWr4ULhO+/A00/DgAHNv0YKI5c1Dt3dP4medo4eTS3nMgqY4e473P09oBrQwjkiJeyXv4RHHoGbboKTT467mvYlpz5pM+toZgsIK4LPcfdXoi9NMrOFZjbNzHpGbf2AD7Jevipq2/09J5hZlZlVZTKZ1v8EIlJUDz0E114bJk6aNCnuatqfnELa3evcfSjQHxhuZkcCtwOHAkOBGuDG6PDGLiXscebt7lPdvdLdK8vKylpRuogU27x54ULhCSfoQmFcWjS6w90/Ap4DRrr72ii864E72dmlsQrI7rHqD6zOv1QRaUvvvw/f+AaUl8OsWVqnMC65jO4oM7MDov3uwAjgLTMrzzrsm8CiaH82MMbMuprZIGAwMK+gVYtIUW3eDF//OmzfHu4o1B+78clldEc5MN3MOhJCfaa7P2pm/21mQwldGSuA7wK4+2IzmwksAWqBiRrZIVI6amvDhElvvQVPPAGHHx53Re1bsyHt7guBYxppH9fEa6YAU/IrTUTamjtcfnm45fvOO+H00+OuSHTHoYj83a9+BbffDj/6EVx6adzVCCikRSRy113w05/Ct74VwlqSQSEtIvzP/8B3vwsjR8K0adBByZAY+qcQaeeefz5M3v/lL8ODD0KXLnFXJNkU0iLt2MKFYSz0oEFhqN0++8RdkexOIS3STr39Npx5JvToEUZzfP7zcVckjdHKLCLt0LvvwmmnQV0dzJ0LBx0Ud0WyNwppkXbm/fdDQG/bBs8+CxUVcVckTVFIi7Qjq1aFgN68OZxBDxkSd0XSHIW0SDtRUxPuIMxkwsT9xx4bd0WSC4W0SDuwejWMGAEffhguEg7XMhwlQyEtknLvvx/OoNeuhccegxNPjLsiaQmFtEiKVVeHPugtW2DOHDj++LgrkpZSSIuk1JIloYvjs8/gmWfgmD3mspRSoJtZRFLo9dfDgrHu8Je/KKBLWS4rs3Qzs3lm9oaZLTaza6P2XmY2x8yWR9ueWa+52syqzWyZmZ1ZzB9ARHb1zDNwyinQvXuYl0PjoEtbLmfSO4DT3P1owqKzI83seGAyMNfdBwNzo+eYWQUwBjgCGAncFq3qIiJFNmNGmMluwAD4619h8OC4K5J8NRvSHnwSPe0cPRwYBUyP2qcD50b7o4AZ7r7D3d8Dqtm5SK2IFMnNN8PYseHi4AsvhKCW0pdTn7SZdTSzBcA6YI67vwL0dfcagGjbJzq8H/BB1stXRW0iUgT19XDVVXDllXDeefDUU9CzZ/Ovk9KQU0i7e527DwX6A8PN7MgmDrfG3mKPg8wmmFmVmVVlMpmcihWRXW3bBv/4j3DjjTBxIjzwAHTrFndVUkgtGt3h7h8BzxH6mteaWTlAtF0XHbYKyP5Dqz+wupH3murule5eWab14kVabPVqOOkkmDkTfv1ruPVW6KirP6mTy+iOMjM7INrvDowA3gJmA+Ojw8YDs6L92cAYM+tqZoOAwcC8Atct0q5VVYWVVJYuDUtf/ehHYI39DSslL5ebWcqB6dEIjQ7ATHd/1MxeBmaa2SXASmA0gLsvNrOZwBKgFpjo7nXFKV+k/Zk5E779bSgrg5de0kx2aWfue3QXt7nKykqvqqqKuwyRRKurg2uugX//d/jKV+BPf4K+feOuqu2Z2Xx3r4y7jrai28JFSkAmEy4QPv00/NM/wW236QJhe6GQFkm4V16B0aNh3Tq46y645JK4K5K2pLk7RBLKPZwxf/WrYdTGSy8poNsjhbRIAm3aBBdcEMY+n3EGzJ+vlVTaK4W0SMK88AIcfTQ8/DD86lcwezb06hV3VRIXhbRIQtTWws9/Hmaw69o1dG9Mngwd9FvarunCoUgCVFfDRRfByy+HMdC33AI9esRdlSSB/o8WiVF9Pfz2t+GGlCVL4A9/gHvuUUDLTjqTFonJ8uVw8cXw4ovwta/Bf/4n9NN8kbIbnUmLtLG6ujD385AhsGgRTJ8OjzyigJbG6UxapA29+ir88z/Da6/B178ezp4PPDDuqiTJdCYt0gY++gj+5V/guOOgpibM+zx7tgJamqeQFikid7jvPviHfwhnzd/7Hrz1Fpx/vqYWldyou0OkSF5+Gf71X8PcG8OHw2OP6a5BaTmdSYsU2Pvvw5gxYTrRlSth2rRwY4oCWlpDZ9IiBbJ5c1jG6qabwl2C//ZvYcWUffeNuzIpZbksnzXAzJ41s6VmttjMvh+1X2NmH5rZguhxdtZrrjazajNbZmZnFvMHEInb1q1w/fUwaFCYa2P0aFi2DH75SwW05C+XM+la4Afu/pqZ9QDmm9mc6Gs3u/sN2QebWQUwBjgCOBB42swO0xJakjbbt8Mdd4RgXrcOzj47BPOwYXFXJmnS7Jm0u9e4+2vR/hZgKdDUsPtRwAx33+Hu7wHVwPBCFCuSBH/7G/zud/DFL4YLg0cdFfqc//xnBbQUXosuHJrZQOAY4JWoaZKZLTSzaWbWM2rrB3yQ9bJVNBLqZjbBzKrMrCqTybS8cpE2tmkTTJkCAwfCpElh+8wzYUmrE06IuzpJq5xD2sz2BR4CrnD3j4HbgUOBoUANcGPDoY28fI/Vbt19qrtXuntlWVlZS+sWaTOrV8MPfwgHHQQ/+xl8+cvw/PNh3udTT427Okm7nEZ3mFlnQkDf7+5/AnD3tVlfvxN4NHq6ChiQ9fL+wOqCVCvSRtzD+OZbb4U//jHMtzFmTBitcfTRcVcn7UkuozsMuBtY6u43ZbWXZx32TWBRtD8bGGNmXc1sEDAYmFe4kkWKZ/v2MOHR8OGhC+ORR+Cyy8KMdfffr4CWtpfLmfSJwDjgTTNbELX9BBhrZkMJXRkrgO8CuPtiM5sJLCGMDJmokR2SdEuXhnmc77kH1q+Hww8PFwfHjdPczhKvZkPa3V+k8X7mx5p4zRRgSh51iRTdxx+HiY6mTYP//V/o1CnMTDdpEpx2mubWkGTQHYfSrnz2WRiR8fvfh77mbdvCWfMNN8CFF0LfvnFXKLIrhbSkXl1dGI3xwAPw0EOhO6NHj9CVcfHFof9ZZ82SVAppSaVPPw1D5GbNCmfMa9bA5z4H55wDF1wAZ50F3brFXaVI8xTSkhobNsDjj4cRGU88Efqcu3YN6wdecEHY7rNP3FWKtIxCWkpWbS3Mnw9z54ZQ/utfw+rbffuGSY7OOQdGjFAwtwfz58/v06lTp7uAIymtKZjrgUW1tbWXDhs2bF1jByikpWS4h6Fyc+eGx3PPhelBIYxf/slPQjBXVoapQqX96NSp011f+MIXDi8rK9vUoUOHPe5wTqr6+nrLZDIVa9asuQv4RmPHKKQlsXbsCAu2vvRSOEt+6SVYG93nOmhQOFseMSLcmt2nT7y1SuyOLLWABujQoYOXlZVtXrNmzZF7O0YhLYngDu++C6+/Hm7HfuklqKoKFwABDj0UzjgDTjoJTj89hLRIlg6lFtANorr3+refQlraXG0tvP12OEt+7bUQzK+/vrProkuX0GVx+eVhCaqvfEXjl6X9UkhL0WzfHsJ46VJYsmTn9u23w00lEIbBHX00jB0b1gA85pgwP3PXrvHWLlIIX/3qVwcvWLBgn8rKyk+effbZ6ta8h0Ja8vLRR/DOO6GrYvfHihVhtAWEC3mHHAIVFeHW6yOOCKH8pS+F27FF0uiqq65as3Xr1g533nlnq+dj1q+HNModNm6EDz/c+2PlyjARfrbevUMYDx8O3/pWCOWKCjjsMN08Iul12WWX9Tv44IM/nTx5cgbgyiuvPLBHjx5111577dpHH300rym6FNLtQF1duLHj449Dv++mTeHW6EwmbBvbz2RCd8Xu+vSBfv1gwIDQV3zIIeGi3iGHhIt5++3X9j+fSLaLL2bAokV8rpDveeSR/G3atF1WnNrFhRdeuPGKK644qCGkZ82a1fOJJ55YXojvrZCOSV1d6JdteNTW7vp8x44w+c+2bWFNvYb9vT3/299CADcEccN28+awmnVT9tsPysrCWXB5OQwZEp4feGAI5IZHeXm4qCciuzrxxBO3bdiwodOKFSs619TUdNp///3rBg8e/Gkh3rvdhPTGjfCLX8B774UQdA/9pdnb1rQ19fXs0N09iL1Ag4W6doXu3cO8FPvvv/Nx0EEhfPfff8/tAQfsDOXevRW8ki5NnfEW0znnnLPpvvvu67lmzZrO55133sZCvW+zIW1mA4D/Ar5AuIVxqrv/XzPrBTwADCRM+n++u2+KXnM1cAlQB1zu7k8WquDW+t73wvSUxxwTQqlDhzDzWcO2Yb9jxz3bdt/Pta1z510fnTrt2ba3r3XpEoK3IYC7d99zv3t33VknkhTjxo3b+J3vfGfgpk2bOv3lL39ZVqj3zeVMuhb4gbu/ZmY9gPlmNgf4NjDX3a83s8nAZODHZlYBjAGOAA4Enjazw+JcnWXLFnjwwTCZ+623xlWFiKRZZWXl9q1bt3bo27fvpwcffPBnAMOGDfvSu+++223btm0d+/btO+S2225bcd55533ckvfNZWWWGsJq4Lj7FjNbCvQDRgGnRIdNB54Dfhy1z3D3HcB7ZlYNDAdebklhhfTGG+HOtZEj46pARNqDt99+e0n28/nz5+d9Rt2iP5bNbCBwDPAK0DcK8IYgb5g9oR/s0ie0Kmrb/b0mmFmVmVVlMplWlJ67N94IWy0iKiKlJueQNrN9gYeAK9y9qdP1xta42OMymbtPdfdKd68sK2v1OO+cLF4cLpb12+O/ChGRZMsppM2sMyGg73f3P0XNa82sPPp6OdAwF+oqYEDWy/sDqwtTbuusXAkDB2qJJJEUq6+vry/J3/Co7vq9fb3ZkDYzA+4Glrr7TVlfmg2Mj/bHA7Oy2seYWVczGwQMBua1ovaCWbkyDEkTkdRalMlk9i+1oI7mk94fWLS3Y3IZ3XEiMA5408wWRG0/Aa4HZprZJcBKYDSAuy82s5nAEsLIkIlxjuyAENInnxxnBSJSTLW1tZeuWbPmrmhe5lIamPr3lVn2dkAuoztepPF+ZoDT9/KaKcCUXCostk8+CXfd9e8fdyUiUizR0lONrmxS6krpf5xWaRg4ovmIRaQUpT6k168P2969461DRKQ1FNIiIgmmkBYRSTCFtIhIgrWLkO7YMUzRKSJSatpFSH/+87rbUERKU7sIaXV1iEipSn1Ib9igkBaR0pX6kNaZtIiUMoW0iEiCpTqk3RXSIlLaUh3SmzdDXZ1CWkRKV6pDWjeyiEipU0iLiCRYLiuzTDOzdWa2KKvtGjP70MwWRI+zs752tZlVm9kyMzuzWIXnQiEtIqUulzPpe4GRjbTf7O5Do8djAGZWAYwBjohec5uZdSxUsS2lkBaRUtdsSLv788DGHN9vFDDD3Xe4+3tANTA8j/ryopAWkVKXT5/0JDNbGHWH9Iza+gEfZB2zKmrbg5lNMLMqM6vKNCyfUmDr10OXLrDvvkV5exGRomttSN8OHAoMBWqAG6P2xqYx8sbewN2nunulu1eWlZW1soymaXIlESl1rQppd1/r7nXuXg/cyc4ujVXAgKxD+wOr8yux9XQji4iUulaFtJmVZz39JtAw8mM2MMbMuprZIGAwMC+/EltPIS0ipa5TcweY2R+AU4DeZrYK+AVwipkNJXRlrAC+C+Dui81sJrAEqAUmuntdUSrPwfr1MGRIXN9dRCR/zYa0u49tpPnuJo6fAkzJp6hC0Zm0iJS61N5xWFcHGzeGC4ciIqUqtSG9cWOYBa9IA0dERNpEakO64UYWhbSIlLLUhnTD/TEKaREpZQppEZEES31Ia3SHiJQyhbSISIKlNqTXr4f99oOuXeOuRESk9VIb0pmM+qNFpPSlOqTV1SEipS7VIa0zaREpdakN6fXrFdIiUvpSGdLuOpMWkXRIZUhv2QKffqo+aREpfakMad1tKCJp0WxIRwvNrjOzRVltvcxsjpktj7Y9s752tZlVm9kyMzuzWIU3RSEtImmRy5n0vcDI3domA3PdfTAwN3qOmVUAY4AjotfcZmYdC1ZtjjQDnoikRbMh7e7PAxt3ax4FTI/2pwPnZrXPcPcd7v4eUM3ORWrbjG4JF5G0aG2fdF93rwGItn2i9n7AB1nHrYra9mBmE8ysysyqMg2pWiDq7hCRtCj0hUNrpM0bO9Ddp7p7pbtXlhU4TWtq4HOfg333Lejbioi0udaG9FozKweItuui9lXAgKzj+gOrW19e69TUQHk5WGP/ZYiIlJDWhvRsYHy0Px6YldU+xsy6mtkgYDAwL78SW27NGjjwwLb+riIihdepuQPM7A/AKUBvM1sF/AK4HphpZpcAK4HRAO6+2MxmAkuAWmCiu9cVqfa9qqmBoUPb+ruKiBResyHt7mP38qXT93L8FGBKPkXlq6YGzjorzgpERAojdXccbt0abgsvL4+7EhGR/KUupGtqwlYhLSJpoJAWEUkwhbSISIIppEVEEix1If3BB+Fuw1694q5ERCR/qQvplSthwADdbSgi6ZDKkD7ooLirEBEpDIW0iEiCpSqkd+wIFw4V0iKSFqkK6Q8/DFuFtIikRapCeuXKsFVIi0haKKRFRBIsVSH97rth6N2AAc0fKyJSClIV0tXV4Sy6a9e4KxERKYxm55NuipmtALYAdUCtu1eaWS/gAWAgsAI439035Vdmbqqr4YtfbIvvJCLSNgpxJn2quw9198ro+WRgrrsPBuZGz9vE8uUKaRFJl2J0d4wCpkf704Fzi/A99rBxY3gopEUkTfINaQeeMrP5ZjYhauvr7jUA0bZPYy80swlmVmVmVZlMJs8y4J13wnbw4LzfSkQkMfLqkwZOdPfVZtYHmGNmb+X6QnefCkwFqKys9DzrYPnysNWZtIikSV5n0u6+OtquAx4GhgNrzawcINquy7fIXCxdCh07wqGHtsV3ExFpG60OaTPbx8x6NOwDZwCLgNnA+Oiw8cCsfIvMxZtvhq6Obt3a4ruJiLSNfLo7+gIPW5i4uRPwe3d/wsxeBWaa2SXASmB0/mU2b9EiOPbYtvhOIiJtp9Uh7e7vAkc30r4BOD2folpq69Zwt+FFF7XldxURKb5U3HG4ZAm4w5FHxl2JiEhhpSKkFy4M26OOircOEZFCS0VIv/IKHHCARnaISPqkJqSPOw46pOKnERHZqeRj7ZNPwsiO446LuxIRkcIr+ZCuqoL6ejj++LgrEREpvJIP6WefDd0cCmkRSaOSD+knn4Thw6Fnz7grEREpvJIO6Y0b4dVX4Ywz4q5ERKQ4Sjqkn3oq9EcrpEUkrUo6pH//e+jXD044Ie5KRESKo2RDesMGePxxGDtW46NFJL1KNt7uvhtqa2HcuLgrEREpnpIM6e3b4eabYcQIGDIk7mpERIon3+WzYnHDDbBmTeiTFhFJs6KdSZvZSDNbZmbVZja5UO+7dClcdx2cfz6cemqh3lVEJJmKEtJm1hH4HXAWUAGMNbOKQry3O5x8MtxySyHeTUQk2YrV3TEcqI5Wb8HMZgCjgCX5vnFFRRgfLSLSHhSru6Mf8EHW81VRm4iItECxQtoaafNdDjCbYGZVZlaVyWSKVIaISGkrVkivAgZkPe8PrM4+wN2nunulu1eWlZUVqQwRkdJWrJB+FRhsZoPMrAswBphdpO8lIpJaRblw6O61ZjYJeBLoCExz98XF+F4iImlWtJtZ3P0x4LFivb+ISHtQkreFi4i0FwppEZEEM3dv/qhiF2GWAd7PauoNrI+pnJZQnYWlOgsrrXUe7O7tZkhYIkJ6d2ZW5e6VcdfRHNVZWKqzsFRnOqi7Q0QkwRTSIiIJltSQnhp3ATlSnYWlOgtLdaZAIvukRUQkSOqZtIiIoJAWEUm0xIV0sZbdamUtK8zsTTNbYGZVUVsvM5tjZsujbc+s46+O6l5mZmcWsa5pZrbOzBZltbW4LjMbFv181WZ2i5k1NsVsoeu8xsw+jD7TBWZ2dgLqHGBmz5rZUjNbbGbfj9oT9Zk2UWeiPlMz62Zm88zsjajOa6P2RH2eJcPdE/MgTMb0DnAI0AV4A6iIsZ4VQO/d2n4DTI72JwO/jvYronq7AoOin6Njkeo6CTgWWJRPXcA84ATC/N+PA2e1QZ3XAFc1cmycdZYDx0b7PYC3o3oS9Zk2UWeiPtPoPfeN9jsDrwDHJ+3zLJVH0s6k/77slrt/CjQsu5Uko4Dp0f504Nys9hnuvsPd3wOqCT9Pwbn788DGfOoys3JgP3d/2cNvw39lvaaYde5NnHXWuPtr0f4WYClhJaFEfaZN1Lk3cdXp7v5J9LRz9HAS9nmWiqSFdNKW3XLgKTObb2YTora+7l4D4ZcG6BO1x117S+vqF+3v3t4WJpnZwqg7pOFP3kTUaWYDgWMIZ3+J/Ux3qxMS9pmaWUczWwCsA+a4e6I/zyRLWkg3u+xWGzvR3Y8lrHo+0cxOauLYpNXeYG91xVXv7cChwFCgBrgxao+9TjPbF3gIuMLdP27q0L3U1Ca1NlJn4j5Td69z96GEVZmGm9mRTRwe+799kiUtpJtddqstufvqaLsOeJjQfbE2+jOMaLsuOjzu2lta16pof/f2onL3tdEvcD1wJzu7hGKt08w6E4Lvfnf/U9ScuM+0sTqT+plGtX0EPAeMJIGfZylIWkgnZtktM9vHzHo07ANnAIuiesZHh40HZkX7s4ExZtbVzAYBgwkXPdpKi+qK/tzcYmbHR1fML8p6TdE0/JJGvkn4TGOtM3rfu4Gl7n5T1pcS9Znurc6kfaZmVmZmB0T73YERwFsk7PMsGXFfudz9AZxNuGr9DvDTGOs4hHDF+Q1gcUMtwOeBucDyaNsr6zU/jepeRhGvQgN/IPxZ+xnhbOOS1tQFVBJ+od8B/h/RHahFrvO/gTeBhYRfzvIE1Pl/CH9GLwQWRI+zk/aZNlFnoj5TYAjwelTPIuDnrf3dKfa/fSk8dFu4iEiCJa27Q0REsiikRUQSTCEtIpJgCmkRkQRTSIuIJJhCWkQkwRTSIiIJ9v8BkJnj9HwL1VUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#List of T and ps we want to do the computation for\n",
    "TsList = numpy.linspace(100.,3000,500)\n",
    "psList = [p0*pp0Fun(TT) for TT in TsList] #Used to determine mass of ocean used\n",
    "\n",
    "#The following statement computes the OLR for each surface pressure ps.\n",
    "#By Clausius-Clapeyron, ps also determines Ts. \n",
    "OLRList = [OLR(TT) for TT in TsList]\n",
    "\n",
    "#Plot results\n",
    "c = Curve()\n",
    "c.addCurve(TsList)\n",
    "c.addCurve(OLRList)\n",
    "w = plot(c)\n",
    "#**ToDo: Plot ps. ps/g is mass of water that has come from ocean\n",
    "#        A double-axis plot would be best, since it could show\n",
    "#        decrease of ocean volume, but matplotlib would be needed for double-axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "fadc2c3f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWkAAAD8CAYAAAC1p1UKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdo0lEQVR4nO3df5RVZb3H8feXHyIp4q9xRCChHG8CKepZaJE/ESUywEibWyYlSRmm3nQl3HKZ3WtZ95bmLXWBEpQWUZJMpCGiiFr8GBAR8BcJwggDo6Io0sjMfO8fzyaOMDBnZs45e59zPq+1ztpn9uy95zt61odnnv3s5zF3R0REkqlD3AWIiMi+KaRFRBJMIS0ikmAKaRGRBFNIi4gkmEJaRCTBOsVdAMCRRx7pffr0ibsMESkAS5cufd3dy+KuI18SEdJ9+vShuro67jJEpACY2atx15BP6u4QEUkwhbSISIIppEVEEkwhLSKSYAppEZEEU0iLiCSYQlpEYjN/PixZEncVyZaIcdIiUppuuAEOOggeeyzuSpJLLWkRiUVjI6xcCSeeGHclyaaQFpFY/OMf8N57cNJJcVeSbAppEYnFihVhq5b0/mUc0mbW0cyeMbPZ0deHm9lcM3s52h6WduxEM1tjZi+a2QW5KFxECtvy5dCxI/TrF3clydaalvQ1wPNpX08A5rl7BTAv+hoz6wdUAv2BYcCdZtYxO+WKSLFYtiwEdNeucVeSbBmFtJn1Aj4D3JO2eyQwLXo/DRiVtn+6u9e7+1pgDTAoK9WKSFFwh6VL4ZRT4q4k+TJtSd8OfAdoSttX7u6bAKLtUdH+nsCGtONqon0iIgBs2gRbtiikM9FiSJvZhcAWd1+a4TWtmX3ezHXHmVm1mVXX1dVleGkRKQbLloWtQrplmbSkBwMjzGwdMB0418zuAzabWQ+AaLslOr4G6J12fi9g454XdfdJ7p5y91RZWckssiAiwDPPgJmG32WixZB294nu3svd+xBuCD7m7pcCVcCY6LAxwKzofRVQaWZdzKwvUAEsznrlIlKwli2D44+Hbt3iriT52vNY+K3ADDMbC6wHLgZw91VmNgNYDTQA4929sd2VikjRWLYMBg+Ou4rC0KqQdvf5wPzo/RvAkH0cdwtwSztrE5Ei9PrrsH49fOtbcVdSGPTEoYjk1TPPhK1uGmZGIS0ieVVdHbYDB8ZaRsFQSItIXi1aFG4aHn543JUUBoW0iOSNOyxcCKedFnclhUMhLSJ5s349bN4Mp58edyWFQyEtInmzcGHYKqQzp5AWkbxZuDDMevfxj8ddSeFQSItI3ixcCKkUdO4cdyWFQyEtInlRXx+eNFRXR+sopEUkL5Yvh/ffV0i3lkJaRPJi101DDb9rHYW0iOTFk0/CscdCTy0B0ioKaRHJOXdYsADOOivuSgqPQlpEcu6FF6CuDs48M+5KCo9CWkRy7oknwlYt6dbLZI3DA81ssZk9a2arzOzmaP/3zew1M1sevYannTPRzNaY2YtmdkEufwERSb4FC6BHD/joR+OupPBkMul/PXCuu79rZp2Bp8zs4eh7t7n7/6YfbGb9CMts9QeOAR41s+O1OotIaXIPLemzzgrrGkrrZLLGobv7u9GXnaPXXqt/pxkJTHf3endfC6wBBrW7UhEpSK+8Ahs3qj+6rTLqkzazjma2nLAi+Fx3XxR96yozW2FmU8zssGhfT2BD2uk10T4RKUHqj26fjELa3RvdfSDQCxhkZgOAu4CPAgOBTcBPo8Ob+4Nmr5a3mY0zs2ozq66rq2tD6SJSCB59FMrL4YQT4q6kMLVqdIe7v0VYiHaYu2+OwrsJmMzuLo0aoHfaab2Ajc1ca5K7p9w9VVZW1pbaRSThmppg7lwYOlT90W2VyeiOMjM7NHrfFTgPeMHMeqQddhGwMnpfBVSaWRcz6wtUAIuzWrWIFITly8Pq4OefH3clhSuT0R09gGlm1pEQ6jPcfbaZ/cbMBhK6MtYBXwdw91VmNgNYDTQA4zWyQ6Q0zZ0btuedF28dhczc9zdQIz9SqZRX71pCWESKxpAhoSX97LPZu6aZLXX3VPaumGx64lBEcmL7dnjqKXV1tJdCWkRyYsGCMH/00KFxV1LYFNIikhMPPxzWMzzjjLgrKWwKaRHJOneYNSu0ort2jbuawqaQFpGsW7EC1q+HESPirqTwKaRFJOuqqsLDKxdeGHclhU8hLSJZV1UVFpwtL4+7ksKnkBaRrHrtNaiuVldHtiikRSSr/vznsFVIZ4dCWkSy6g9/gOOP16x32aKQFpGsqa2F+fOhslKz3mWLQlpEsuaPfwzTk37hC3FXUjwU0iKSNdOnw8c/Dv36xV1J8VBIi0hWbNgATz+tVnS2KaRFJCumTw9bhXR2ZbIyy4FmttjMnjWzVWZ2c7T/cDOba2YvR9vD0s6ZaGZrzOxFM7sgl7+AiMTPHaZOhU98Ao47Lu5qiksmLel64Fx3P4mw6OwwMzsdmADMc/cKYF70NWbWD6gE+gPDgDujVV1EpEgtXgyrV8Pll8ddSfFpMaQ9eDf6snP0cmAkMC3aPw0YFb0fCUx393p3XwusYfcitSJShKZMCbPdXXJJ3JUUn4z6pM2so5ktB7YAc919EVDu7psAou1R0eE9gQ1pp9dE+0SkCL33XuiPvvhiOOSQuKspPhmFtLs3uvtAoBcwyMwG7Ofw5oaw77WQopmNM7NqM6uuq6vLqFgRSZ6ZM2HbNnV15EqrRne4+1vAfEJf82Yz6wEQbbdEh9UAvdNO6wVsbOZak9w95e6psrKy1lcuIonwy19CRYVWYMmVTEZ3lJnZodH7rsB5wAtAFTAmOmwMMCt6XwVUmlkXM+sLVACLs1y3iCTAkiWwcCF861vQQQN6c6JTBsf0AKZFIzQ6ADPcfbaZ/R2YYWZjgfXAxQDuvsrMZgCrgQZgvLs35qZ8EYnT//0fdOsGY8a0fKy0TYsh7e4rgJOb2f8GMGQf59wC3NLu6kQksWprww3DK6/UDcNc0h8oItImd98NO3fC+PFxV1LcFNIi0mrvvBO6OkaMCHNHS+4opEWk1e6+G958E7773bgrKX4KaRFplR074Kc/hfPOg0F6ljjnMhndISLyL1OmwObNu2e9k9xSS1pEMvbee3DLLTB4MJx1VtzVlAa1pEUkY7ffDps2wYwZWsMwX9SSFpGMvP46/PjHYUTHpz4VdzWlQyEtIhn5r/+Cd9+FH/0o7kpKi0JaRFr03HNhIqUrrtAis/mmkBaR/XKHb34TDj003DSU/NKNQxHZr1//Gp56Cu65B444Iu5qSo9a0iKyT7W1cN11YYHZr3417mpKk0JaRJrlHvqgt28PD7Bovuh4qLtDRJo1dSrMng233QYf+1jc1ZQu/dsoInt58UW4+urwVOHVV8ddTWnLZPms3mb2uJk9b2arzOyaaP/3zew1M1sevYannTPRzNaY2YtmdkEufwERya7t22H0aDjwQPjNb9TNEbdMujsagOvcfZmZdQOWmtnc6Hu3ufv/ph9sZv2ASqA/cAzwqJkdryW0RJLPHb7xDVi9GubMgd69Wz5HcqvFfyPdfZO7L4vevwM8D/TczykjgenuXu/ua4E1gCY0FCkAP/kJ3HcffP/7MHRo3NUItLJP2sz6ENY7XBTtusrMVpjZFDM7LNrXE9iQdloNzYS6mY0zs2ozq66rq2t95SKSVdOnw4QJUFkJ3/te3NXILhmHtJkdDDwAXOvu24C7gI8CA4FNwE93HdrM6b7XDvdJ7p5y91RZWVlr6xaRLHr88bDi95lnhlEd6odOjoz+V5hZZ0JA3+/uMwHcfbO7N7p7EzCZ3V0aNUB6T1YvYGP2ShaRbJo/Hz7zGaiogD/9Cbp0ibsiSZfJ6A4D7gWed/efpe3vkXbYRcDK6H0VUGlmXcysL1ABLM5eySKSLU88EQK6b1947DE4/PC4K5I9ZTK6YzDwZeA5M1se7ftP4N/NbCChK2Md8HUAd19lZjOA1YSRIeM1skMkeR54AL70JfjIR0JAH3VU3BVJc1oMaXd/iub7mR/azzm3AJovSySh7rgDrr0WTj8dqqrgyCPjrkj2RbcHRErIjh1w+eVwzTUwahTMm6eATjqFtEiJWLMmzGb3q1+FIXZ/+AN07Rp3VdISTbAkUuSamsKqKhMmhEe9H3oIPv3puKuSTKklLVLEVq+Gc84JkySdeSY8+6wCutCoJS1ShN54Izzafddd0K1bmA/6K18Ba24IQBFYunTpUZ06dboHGEBhNT6bgJUNDQ1fO/XUU7c0d4BCWqSIbN0aRm7cfjts2xYmS7r55uK/OdipU6d7jj766BPKysq2dujQYa8nnJOqqanJ6urq+tXW1t4DjGjuGIW0SBGoqYFf/ALuvBPeeQdGjAiLxg4YEHdleTOg0AIaoEOHDl5WVvZ2bW3tPv9PKaRFClRjY5hO9O674S9/CdOMfuELMHEinHhi3NXlXYdCC+hdorr32UWjkBYpII2N8OSTMGMGzJwJmzdDeXkYufG1r4XHu6W4KKRFEm7TJpg7N7weeQS2bAnjmz/zmdByHjECDjgg7iqlOWeccUbF8uXLD0qlUu8+/vjja9pyDYW0SII0NMCqVbBoUXgtXBiG0UG4+XfeeXDRRSGgDzoo3lqlZddff33t9u3bO0yePLnN8zErpEVi0NgYbvatXh1Cedd21aqwxiCEGekGDYLLLoPzz4eTTtI8z0l15ZVX9jz22GPfnzBhQh3At7/97WO6devWePPNN2+ePXt2t/ZcWyEtkmX19VBXF7oltmyB2lpYvx7WrYNXXw3bDRtg587d55SXQ//+YV6NQYPgtNPguOOKd1xzLl1+Ob1XruRD2bzmgAG8N2XKB1ac+oBLL730zWuvvfbDu0J61qxZh/31r399ORs/WyEtJc09dDH885/7fu3YEYa1bdu29+vtt8P2rbd2B/O2bc3/rB49oE+fEMKXXBJu8vXrByecAEcckc/fWrJt8ODBO954441O69at67xp06ZO3bt3b6yoqHg/G9dWSBeJt9+Gp56CF14ILbn33w8tNY8GJbl/8H1z+7JxbC5/RlNT6CZoaGj/9v33d4dwU1Pr/3t/6ENwyCG7X927QyoV5mRu7tWzZ5g3Q3Jvfy3eXPrsZz+79b777justra28+jRo9/M1nVbDGkz6w38Gjia8AjjJHf/uZkdDvwe6EOY9P8Sd98anTMRGAs0Ale7+5xsFSx7mzoVvvnN0OJL17Fj+HN515/Me75vbl++j23NeR06QKdO4fdqbtuly/6/n77t0iWEZqavbt0+GMqd1LyRPXz5y19+84orruizdevWTk888cSL2bpuJh+1BuA6d19mZt2ApWY2F/gKMM/dbzWzCcAE4AYz6wdUAv2BY4BHzex4rc6SGzfeCP/932ESnRtvhJNPDq28zp3VnymST6lU6p/bt2/vUF5e/v6xxx67E+DUU0/9t1deeeXAHTt2dCwvLz/xzjvvXDd69Oh9dIg1L5OVWTYRVgPH3d8xs+eBnsBI4OzosGnAfOCGaP90d68H1prZGsIitX9vTWHSsiVLwqO/l10G996r1p1I3F566aXV6V8vXbq03S3qVg3oMbM+wMnAIqA8CvBdQb5rhbSe8IE+oZpo357XGmdm1WZWXVdX14bSS1tDA3z963D00WFCHQW0SHHKOKTN7GDgAeBad99fc725P7L3eqbe3Se5e8rdU2VlbR7nXbLuvx+eeQZ+/vNw00pEilNGIW1mnQkBfb+7z4x2bzazHtH3ewC75kKtAXqnnd4L2JidcmWXu++Gj30MPv/5uCsRSYSmpqamgrwLE9W9zzFGLYa0mRlwL/C8u/8s7VtVwJjo/RhgVtr+SjPrYmZ9gQpgcRtql31YsSI8LjxunG4OikRW1tXVdS+0oI7mk+4OrNzXMZn0ZA4Gvgw8Z2bLo33/CdwKzDCzscB64GIAd19lZjOA1YSRIeM1siO7Jk8OQ8guuyzuSkSSoaGh4Wu1tbX3RPMyF9LD8/9amWVfB5j7Xt3FeZdKpby6ujruMgpCfX14hPjCC+G+++KuRiT/zGypu6firiNfCulfHAEeeyw8XfjFL8ZdiYjkg0K6wDz4IBx8MAwZEnclIpIPCukC0tQEs2bB8OGhT1pEip9CuoAsWhSWSxo1Ku5KRCRfFNIF5MEHw5wcw4fHXYmI5ItCuoDMng1nnaUnDEVKiUK6QGzYEJZYGjYs7kpEJJ8U0gXikUfC9oIL4q1DRPJLIV0g5swJq3v07x93JSKSTwrpAtDYCI8+GlaM1lwdIqVFIV0Aqqth69YQ0iJSWhTSBWDOnNCCHjo07kpEJN8U0gVgzpywEvURR8RdiYjkm0I64d56KzxpqFEdIqVJIZ1w8+aFG4cKaZHSlMnKLFPMbIuZrUzb930ze83Mlkev4Wnfm2hma8zsRTNTtLTTI49At25w2mlxVyIiccikJT0VaO45t9vcfWD0egjAzPoBlUD/6Jw7zaxjtootNe6hP3rIkDBnh4iUnhZD2t0XAG9meL2RwHR3r3f3tcAaYFA76itpL70Er76qrg6RUtaePumrzGxF1B1yWLSvJ7Ah7ZiaaN9ezGycmVWbWXVdXV07yihec+aErUJapHS1NaTvAj4KDAQ2AT+N9jf3PFyziyi6+yR3T7l7qqysrI1lFLc5c6CiAvr2jbsSEYlLm0La3Te7e6O7NwGT2d2lUQP0Tju0F7CxfSWWpvp6mD9fTxmKlLo2hbSZ9Uj78iJg18iPKqDSzLqYWV+gAljcvhJL09NPw3vvqatDpNR1aukAM/sdcDZwpJnVADcBZ5vZQEJXxjrg6wDuvsrMZgCrgQZgvLs35qTyIjdnThjRcc45cVciInEy92a7jPMqlUp5dXV13GUkysCBcNhh8PjjcVcikixmttTdU3HXkS964jCBamvh2WfV1SEiCulEeuihsNVSWSKikE6g2bOhVy846aS4KxGRuCmkE6a+PszXceGFWoVFRBTSifPEE7B9ewhpERGFdMLMng1du8K558ZdiYgkgUI6QdxDSA8ZEoJaREQhnSDPPQdr16qrQ0R2U0gnyAMPhJuFo0bFXYmIJIVCOkFmzoQzzoDy8rgrEZGkUEgnxEsvwcqVMHp03JWISJIopBPigQfC9qKL4q1DRJJFIZ0QM2fCoEHQu3fLx4pI6VBIJ8A//gHV1fD5z8ddiYgkjUI6Ae6/P4zqqKyMuxIRSZoWQzpaaHaLma1M23e4mc01s5ej7WFp35toZmvM7EUz02SbLXAPIX3WWerqEJG9ZdKSngrsOWnmBGCeu1cA86KvMbN+QCXQPzrnTjPrmLVqi9DSpWFkx5e+FHclIpJELYa0uy8A3txj90hgWvR+GjAqbf90d69397XAGnYvUivNuP9+OOAA9UeLSPPa2idd7u6bAKLtUdH+nsCGtONqon3SjIYG+N3vwmPghx4adzUikkTZvnHY3AzIzS6iaGbjzKzazKrr6uqyXEZh+MtfYPNmGDMm7kpEJKnaGtKbzawHQLTdEu2vAdJvf/UCNjZ3AXef5O4pd0+VlZW1sYzCNnkyHHMMDB8edyUiklRtDekqYFf7bwwwK21/pZl1MbO+QAWwuH0lFqcNG+Dhh+Hyy6FTp7irEZGkajEezOx3wNnAkWZWA9wE3ArMMLOxwHrgYgB3X2VmM4DVQAMw3t0bc1R7QZsyJQy/Gzs27kpEJMnMvdku47xKpVJeXV0ddxl509gIffvCCSfAnDlxVyNSWMxsqbun4q4jX/TEYQwefDB0d3zjG3FXIiJJp5COwc9+Bh/5CIwYEXclIpJ0Cuk8W7gQ/vY3uOYa6KhnMUWkBQrpPLvtNujeHb761bgrEZFCoJDOo7Vr4Y9/hHHjoFu3uKsRkUKgkM6jH/4wjIm++uq4KxGRQqGQzpO1a2Hq1NCK7tUr7mpEpFAopPPkhz+EDh1gwoS4KxGRQqKQzoP0VnRPzQkoIq2gkM6D73439EWrFS0iraWQzrG//z3MGX399WpFi0jrKaRzqKkJ/uM/oEcPuOGGuKsRkUKkSTJzaPp0WLQozHh38MFxVyMihUgt6RzZuhWuuw5OOUUrr4hI26klnSPXXw91dfDQQ2HonYhIW7QrpM1sHfAO0Ag0uHvKzA4Hfg/0AdYBl7j71vaVWVjmzQtdHDfcACefHHc1IlLIstHGO8fdB6ZNwj0BmOfuFcC86OuSsW0bXHEFHHcc3HRT3NWISKHLxR/iI4Fp0ftpwKgc/IxEcocrr4RXX4Vf/Qq6do27IhEpdO0NaQceMbOlZjYu2lfu7psAou1R7fwZBWPaNPjtb+Hmm+FTn4q7GhEpBu29cTjY3Tea2VHAXDN7IdMTo1AfB/DhD3+4nWXEb/VqGD8ezjkHJk6MuxoRKRbtakm7+8ZouwX4EzAI2GxmPQCi7ZZ9nDvJ3VPuniorK2tPGbF7/XX47GfDHNH33acVV0Qke9oc0mZ2kJl12/UeOB9YCVQBu0YGjwFmtbfIJKuvh899Dl57DWbNgmOOibsiESkm7enuKAf+ZGa7rvNbd/+rmS0BZpjZWGA9cHH7y0ympqYwkuPJJ8P8HKedFndFIlJs2hzS7v4KcFIz+98AhrSnqELgDlddBb/5DfzgB1BZGXdFIlKM9CxcG7iHiZPuugu+8x343vfirkhEipUeC2+lxka49lr4xS/gmmvg1lsh9PiIiGSfQroVduyASy+FmTPD5En/8z8KaBHJLYV0hjZvhtGj4W9/g9tvD61oEZFcU0hnYMGCcGNw69YwR/Qll8RdkYiUCt043I/GxtDnfO65YdL+RYsU0CKSX2pJ78PKlXD55bBkCVx8MdxzDxxySNxViUipUUt6D+++CzfeGFZUWbs2PKTy+98roEUkHmpJRxoawkT9N90EtbXwxS+GG4QFPq2IiBS4kg/p+vowKdJPfgIvvQSf/GQYYveJT8RdmYhICYf0li1hYv477oCNG2HgQHjgAbjoIo19FpHkKKmQ3rkTHnkE7r0X/vzn0MVx7rkhrIcOVTiLSPIUfUjv2AFz54ZWclUVvPVW6Ge+5hoYOxZOOCHuCkVE9q3oQnrnTnjmmbBi96OPwtNPh37nQw+FESPC3M+f/jQccEDclYqItKygQ3rnTnj5ZVixAhYvDg+bLFsG//xn+P6JJ4YlrYYNg7PPhs6dYy1XRKTVchbSZjYM+DnQEbjH3W/NxnWfew5+9KPwsMkLL4SgBjjwwDC2+cor4fTTQygfVTJL4IpIscpJSJtZR+CXwFCgBlhiZlXuvrq91965M0xyNGAADB8etgMGQP/+aimLSPHJVUt6ELAmWr0FM5sOjATaHdKnnALr1rX3KiIihSFXj4X3BDakfV0T7RMRkVbIVUg3N+LYP3CA2Tgzqzaz6rq6uhyVISJS2HIV0jVA77SvewEb0w9w90nunnL3VJkmyBARaVauQnoJUGFmfc3sAKASqMrRzxIRKVo5uXHo7g1mdhUwhzAEb4q7r8rFzxIRKWY5Gyft7g8BD+Xq+iIipUCT/ouIJJhCWkQkwczdWz4q10WY1QGv7rG7O/D2Pk45Eng9p0Xlx/5+x0L6me29ZlvOb805mRzb3mP0mczfzzzW3UtnSJi7J/IFTNrP96rjri/Xv2Mh/cz2XrMt57fmnEyObe8x+kwW1s8spFeSuzv+HHcBeRDH75iLn9nea7bl/Nack8mx2Tqm0BXLZ7JoJKK7o7XMrNrdU3HXIbKLPpOSK0luSe/PpLgLENmDPpOSEwXZkhYRKRWF2pIWESkJCmkRkQRTSIuIJFhRhLSZjTKzyWY2y8zOj7seETM7wczuNrM/mtmVcdcjhSuxIW1mU8xsi5mt3GP/MDN70czWmNkEAHd/0N2vAL4CfCGGcqUEtPIz+by7fwO4BNDQPGmzxIY0MBUYlr4jbYHbTwP9gH83s35ph3wv+r5ILkylFZ9JMxsBPAXMy2+ZUkwSG9LuvgB4c4/d/1rg1t3fB6YDIy34MfCwuy/Ld61SGlrzmYyOr3L3TwJfym+lUkxyNp90jjS3wO1pwLeA84DuZnacu98dR3FSkpr9TJrZ2cDngC5oXnVph0IL6WYXuHX3O4A78l2MCPv+TM4H5ue3FClGie3u2IcWF7gVyTN9JiWnCi2ktcCtJI0+k5JTiQ1pM/sd8Hfg38ysxszGunsDsGuB2+eBGa4FbiVP9JmUOGiCJRGRBEtsS1pERBTSIiKJppAWEUkwhbSISIIppEVEEkwhLSKSYAppEZEEU0iLiCSYQlpEJMH+H2WnPwQHnxrGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<ClimateGraphicsMPL.plotObj at 0x7fb2180fb730>"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c.XlogAxis = True\n",
    "plot(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91c3eab3",
   "metadata": {},
   "source": [
    "## Post-runaway behaviour, Part I\n",
    "\n",
    "After the ocean has either completely evaporated into the atmosphere, or has gone supercritical so that the ocean surface is no longer defined, the atmosphere is on a dry adiabat until altitude where condensation sets in and the $T(p)$ profile follows the single-component saturated condensing adiabat.  The pressure at that boundary goes to zero as the surface is made hotter. In the following calculation, we assume the atmosphere to be on the dry adiabat from the surface up to the condensing layer, and compute the resulting OLR curve that continues on from the previous case once the ocean is used up.  The mass of the ocean determines when we switch to the profile with the dry adiabatic unsaturated layer.\n",
    "\n",
    "**ToDo: Discuss issue of behaviour when ocean goes supercritical before it is used up. The challenge in presenting the results is that the \"surface\" jumps from what used to be the ocean surface to what used to be the bottom of the ocean. One solution is to hold the pressure of the former \"surface\" fixed, and extrapolate both upwards and downwards along the adiabat. \n",
    "\n",
    "**ToDo: Another issue when we go to temperatures that encounter the critical point is that the constant L form of the saturated adiabat is no longer valid. Maybe we should limit the calculation to the subcritical case, and do the critical case later. Note that the Earth's ocean would just touch the critical point at the point ocean is exhausted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "560e5d13",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define the ocean mass, in terms of equivalent surface pressure. For a pure\n",
    "#condensible atmosphere, the mass per square m of the ocean is psMax/g\n",
    "#When ps exceeds psMax, we pin the surface pressure at psMax and then use the\n",
    "#dry adiabat for the lower atmosphere\n",
    "psMax = 100.e5 #about half the mass of the Earth ocean\n",
    "\n",
    "#Define R/cp\n",
    "Rcp = phys.water.Rcp #Note: Should use the temperature-dependent adiabat.\n",
    "#Function to compute the dry adiabat starting from the surface temperature.\n",
    "#Because the ocean has been depleted, the surface pressure is now fixed.\n",
    "#Caution: The use of pp0Fun in OLR(Tin) would need changing also. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "7a3dc133",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100.0"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psMax/1.e5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "bddf7219",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7377643630055348"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad4da421",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
