{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example - 03 - Venus Aerocapture: Part 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, you will learn to create a vehicle object in AMAT and simulate its aerocapture trajectory in the Venusian atmosphere. \n",
    "\n",
    "For reference, we will re-create a few results from the paper \"Craig and Lyne, Parametric Study of Aerocapture for Missions to Venus, Journal of Spacecraft and Rockets, Vol. 42, No. 6., 2005. DOI:10.2514/1.2589"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us re-use the code from example-01, 02 to create a planet object for Venus and add an atmosphere model. In addition to Planet, import the Vehicle class from AMAT for this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from AMAT.planet import Planet\n",
    "from AMAT.vehicle import Vehicle\n",
    "\n",
    "# Create a planet object\n",
    "planet=Planet(\"VENUS\")    \n",
    "\n",
    "# Load an nominal atmospheric profile with height, temp, pressure, density data\n",
    "planet.loadAtmosphereModel('../atmdata/Venus/venus-gram-avg.dat', 0 , 1 ,2, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a vehicle object flying in the target planet atmosphere.\n",
    "# with params m = 300 kg, beta = 78.0, L/D = 0.35, A = 3.1416, AoA = 0, RN = 1.54\n",
    "\n",
    "# These values are taken from the reference article mentioned above.\n",
    "vehicle=Vehicle('Apollo', 300.0, 78.0, 0.35, 3.1416, 0.0, 1.54, planet)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set initial conditions for the vehicle at atmospheric entry interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# h0 = 180 km, LON = 0 deg, LAT = 0 deg\n",
    "# v0 = 12 km/s, HDG = 0 deg, FPA = 0 deg\n",
    "# DOWNRANGE0 = 0 deg, HEATLOAD0 = 0.\n",
    "\n",
    "# See help(vehicle) for more details.\n",
    "vehicle.setInitialState(180.0,0.0,0.0,12.0,0.0,-4.5,0.0,0.0)\n",
    "\n",
    "# Set solver tolerance = 1E-6 (recommended value)\n",
    "# Setting this too low can result in long execution times.\n",
    "vehicle.setSolverParams(1E-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the overshoot and undershoot limit EFPA\n",
    "\n",
    "# Set max. propogation time = 2400.0 secs.\n",
    "# Set max. time step = 0.1 sec.\n",
    "# Set low value for guess = -80.0 deg\n",
    "# Set high value for guess = -4.0 deg\n",
    "# Set EFPA tolerance = 1E-10 (recommended)\n",
    "# Set target apoapsis = 407 km\n",
    "\n",
    "# This calculation migt take a couple of minutes. Hang on!\n",
    "overShootLimit, exitflag_os  = vehicle.findOverShootLimit (2400.0,0.1,-80.0,-4.0,1E-10,407.0)\n",
    "underShootLimit,exitflag_us  = vehicle.findUnderShootLimit(2400.0,0.1,-80.0,-4.0,1E-10,407.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "1.0\n"
     ]
    }
   ],
   "source": [
    "# exitflag_os and exitflag_us will be set to 1 if a solution was found. Otherwise, it will be 0.\n",
    "print(exitflag_os)\n",
    "print(exitflag_us)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overshoot  limit : -7.0519 deg\n",
      "Undershoot limit : -9.4396 deg\n"
     ]
    }
   ],
   "source": [
    "# print the overshoot and undershoot limits we just computed.\n",
    "print(\"Overshoot  limit : \"+str('{:.4f}'.format(overShootLimit))+ \" deg\")\n",
    "print(\"Undershoot limit : \"+str('{:.4f}'.format(underShootLimit))+ \" deg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the limiting flight path angles for our vehicle at Venus. Let us now calculate these trajectories, and their associated deceleration and heating profiles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOQAAADaCAYAAAC2J+5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2deZyT1fX/3yeZJZPZGGaGHdkE60ZFxrpAEXdrK1rcFRWt4lqraKti1dpqa7Uu1bbitGpbRHHBr0t/7qKitC64gYqKoiAIzL4lk0yW8/vjSTIDzJLMJJNl7vv1ymvyLHmek0w+Ofeee+65oqoYDIbUwJZsAwwGQztGkAZDCmEEaTCkEEaQBkMKYQRpMKQQRpAGQwqRlWwDeqKsrEzHjh3b5fHq6mrKy8v7z6AeMPZ0j7EH3nvvvRpV7fymqprSj6lTp2p39HS8vzH2dI+xRxVYqV18302T1WBIIYwgDYYUIu0FOW/evGSbsA3Gnu4x9nSPaIJyWUVkBnC9qh4iIkOBC4BVwNTQfr+IXA5UAcWq+pfOrlNRUaErV65MiI0GQzIQkfdUtaKzYwnzkKq6HMgLbZ4KfKCqTwA5wJ4iMh0oVdVFQImI7BvL9RtbW/n1K0tZ8NLSuNrdkf/33XP8t+Z/+IK+hN3DkFn8eWkj32zp/fcl0U3WttDf14EFIrITUA98BBwFrAkd/zS0HTX1rW4qWy7h7/UL4mXrNgQ1yLyVF3DsiuPxqz8h9zBkFm996uHye+rZ98LNeNt61/Lslz6kqr4PLAUeBt5T1SBQhiVOAA8wrLPXVldXU1FREXlUVlYCMNjptE7I8pKIZndtWx0+9VGSPYg8e17PLzAMeG5Z0gjAOT8uJDdHtjlWWVkZ+Q5jffc7pV8SA0TkB0AD8BPgNRFZB1QDIVVRCNR29try8nI660MW5jpQFcTuw93mJz83O642b27dDMCwvOFxva4hM/l0fRtP/7cVR45wyeyiHY7PmzcvEkASkZqurtNfUdYDgM2qWg88AOwOPAtMDh3fDXg+lguKCBLIAaDW5YmfpSE2eyxBDnd06rgNhm340yOWdzzryAKGlth7fZ2ECVJE9gQmiMgewIPANBH5EVag51lVXQF4ROQsoCEUBIqNQC4AdS533OwOs7l1CwAjjIc09MCGrX4eesWF3QbzT9jRO8ZCwpqsqroaGN1h11Whv891OOfGvtzDFnQQpIk6d2tfLtMpWzyWIIcZD2nogdsfa8QfgFMOzmfc8L51ndI6McCuloesb42/IDeHBGmarIbu2Fof4B/PtgBw5SnFfb6eEWQXfBcK6gzPM4I0dM0djzfiaVNmHZDHHuNy+ny9tBZkFg4AmjymyWrof+qaAix8uhmAq08dFJdrprUgc7A8ZKM3AU1WE9Qx9MDd/9dES6ty2FQH+3wvNy7XTG9BiuUhm+MsyBa/iyZ/E7m2XEqyS+J6bUNm0OwO8pcnQ97xtPh4R0h3QdpCgmyL7zjklnBSgGMYItLD2YaByMJnmqlvDjJtj1xmTHbE7bppLUhHSJAtvvh6yEiE1QR0DJ3Q6g1yx2NWIsCCOHpHSHdB2i1Bunzx9ZBhQY5wmP6jYUfue7aFqoYgUyflcHhF/LwjpLkgnVnWh+H2x9dDbgkFdIYZD2nYjjaf8qdHLe949anFce/SZIQgW/3euF73O5PHauiCRS+1sLE6wO5js5l1gLPnF8RIWgsyP9uaFtUajHMfstVk6Rh2xB/QyBSrq04pxmaLf8AvrQVZkG15SG8gvh4ykhRgmqyGDjz6mouvvvMzYUQWJ8zMT8g90luQuZaH9AYTFdQZEdfrGtKXQEC5cVEDYOWsZtkTMxyW1oIsyrE8ZJvGT5D+oJ8qTxWCMNQxJG7XNaQ3S1518cVGP+OGZXH6YQUJu09aC7LYYXWqfcRPkFXeKoIEKcstI9sW3yoEhvTEH1BuetDyjgtOKyY7K3HJImkuSMtD+iV+gjQBHcP2LFlmecfxw7OYk0DvCGkuyEEhDxkgfoWuTJaOoSP96R0hzQVZEOpDqt1LW5xKp271bAXMtCuDxcOvuFi7yYqsJto7QpoLMlyeUbK8NLcG43LNLaZSgCGEP6DctLjdOyYqstqRtBZkOJeVLC8tcRLkZjMx2RDioVdcfBnyjqcdmnjvCGkuyLyQICXLS0trfPqQWyJN1qFxuZ4hPenYd7xmzqB+8Y6Q5oJ0bCPIODVZTWK5AVj8cgtffedn55FZnHpIYrJyOiO9BWkLlfi3Gw9piB8+v/L7xVbO6jWn9Z93hHQXZJw9pMvvNqU7DPzrBcs7ThyZxSn96B0hzQUZ6UNme2l2912QW0MBnaGOoaZ0xwCl1Rvkd6Gc1RvmlvSrd4Q0F6RNbNjUSm9raO17to5prhoWPtPMppoA35+QzfEHxn++Y0+ktSAB7Gp5ycY41GY1ghzYNLmC3PyQ1Xf87VklCZnv2BNpL8jsONZmbR+DNIIciPz5iSZqm4Lsv3suR+2bnDVB01+QodqsTd54NFnDgjTFrQYatY0Bbg9VkrvpZyVJiyEkcjm6GSLySoftHBG5QEQOFZFSEbGJyPUiMkdEzuztfXLDgvT0fUm6rabJOmD545JGmt3K4RWOuNZZjZWECTK03mMegFg/N3cDT6rqy6paC5yKtYjrg8D+IjK666t1Ta7NarI2xaFYcnjqlRHkwGJTjZ+/PWVVIf/d2ckd7kp0k7Ut9PdkYDBwpohcGRLoUcCa0PG1wKGdXaC6ujqyNntFRQWVlZXbHA8nmDe3xS+oY6ZeDSx+/6C1gtXsHzqZOik+a3RsT2VlZeQ7DJR1dV7CFmzdjmOBu1V1uYi8gCXAMqA+dNwDdKqC8vJyVq5c2eWF87Ic4IOWPnpIVY00WYfmGg85UPjqOx/3PdeMzQY3zI1vFfKOzJs3j3nz5gEgIjVdnddfQR0H0BR6/h9gD6AaCA/0FAK1vblwuDarq4/LCTT4GvAEPRRkFVCQ3T+Z/Ybkc90DDfgDMOfQfHYd0/f1HftKfwnyTWBK6Hk28C7wLDA5tG8S8HJvLhwuBen2981DRpqrZtrVgOHdz7w88qqL3Gz4zZmpkSqZyCjrnsAEEdkDuBOoEJHZQJuqPgcsAcaLyNnAClVd15v7FOSEiiUH4iNIE9AZGKgqv6qsA+CS2UXsNLS/em/dkzArVHU10DFyetF2xwPAgr7epzAkSE/Ag6r2evwoPO1qqBHkgOA/b7XyxiovpUU2rjylONnmREj7xABnqMmK3Yvb0/spWO2lO0xSQKbjDyhX/92KJ14zZxCDCuxJtqidtBdkXocyHk19mPFhmqwDh/ufa+GzDT7GD8/i/KMLk23ONqS9IMOLtkqWl0ZX7wUZyWPNM4LMZFpag9zwL8s73vSzEnKyU2uaXfoLMjwn0u6lyd37Jqsp/zgwuO3RRrbWB9nnezlJmV7VExkjSLK8NPXBQ24x1eYyns21fm5/zBoOv+W8wSk5CT3tBdmx8lxvqwYENECVtxrALLCTwfzmnw24PMqsA/L44Z7JSyDvjrQXZLuHbOt1UKfGW0NAA5TllJJjS362hiH+fLDWy/3Pt5Blhz+cOzjZ5nRJ+gvS1l69vLd9yHCE1YxBZiaqymV/q0MVLjq2iF1Gp+6qZukvSHvfo6zt065M/zETWbrczZurvZQV27j29NRJAuiMjBEk9t4HdbaYFa8yllZvkF/da6XI/e6skpRKAuiMtBdkPII6psmaudz2WBMbqqwqcmf/KPVn8aS9IDsmBvS+D2kqBWQiG6v93LLEqpNz2wWDsfdzjdXekPaCzMsKVQfL8vShyWqmXmUiC/5Rj9ujHDfDycy9klNFLlbSX5DhNSKze5/LapICMo//ferhoVesuY5/nJcacx2jIXMEmeXpQx/STL3KJAIB5ZK7rEDO/BOKGTssdYc5tifq+ZAiUg4cAoRHzqeo6mUJsSoG2tf38NDgCsT8em/AS11bPXaxU57bZe0hQxqx8JlmPviyjZ2G2LkqheY6RkMsHnIpMAiQ0CMl4sd2sZMjVqWw5l4US97qrQJgaO4QbJL2DYYBz5Y6P9feb83muPOiUvLz0ut/GkvFgFdVdWF4Q0ReSoA9vcKZlUebz0tzmzvmqgEmKSCz+NW99TS5lR/vl8fRB6RHIKcjsQjySxF5HKt6nADjgQMTYlWMOO15NPgaCNo8tHoVpyN6QW41/ceM4dUPW3noFReOHOHOi1JzNkdPxCLIc4HzAVdoe3r8zekd4cAO2Vb6nNMRfTPFFEfODNp87YGcq08tZtzw9AnkdCQWQS4CvlBVP4CIpEaZLtrHIiXLQ6MryPDS6F9rVrzKDO5c2sSaDT4mjsziihPTK5DTkVhEdQZwroi00N5kHZMQq2LEabdmfkt2K/UtsQ19mFo66c/6rX5ufNBa9fiuS0rJzUm/pmqYWAR5kaquCm+IyOTuTu5P2scivdQ3xybIrWYJurRGVbnwzlrcHuWEA50cNjX9AjkdiSUmvNt22+fF05C+EOlDZnlojNFDmhWv0puHl7l44d1WSgpt3HFR6k48jpZYPOThIrIaaAT+CfgSYlEvcHbI1jFN1oFDdUOAy/5qBXJuPa+EYYNTJqzRa6J+B6p6toicBPwKuBT4NmFWxYgzK9yH9MTUZG3xteAKuMizOyjOTt9AwEDl8nvqqG0KctAUB2cekfpTq6KhW0GKyGtYzVp/eBfgBn4LTARGJdK4aOnYZG2IYcbH5sgY5LC0HLMayDz3jpuHXnGRlyssvKw0Y/5/PXnISzoGcjoiIpMSYE+vaJ/x4aEhBg9p5kGmJ83uIBfeYa1e+JszBzFhRHqOOXZGt0GdrsQYOvZF/M3pHR2rBsTShzT9x/Tk2vvr+bY6wNRJOfziuKJkmxNXErkc3QwReWW7fTNF5L7Qc5uIXC8ic0TkzL7cKzIOmeWhoReCNGlz6cOrH7bylyebybLDvfNLyUqDKgCx0GtBishoEely8E5VlwN5Hc4vBQ6mfZbIqcBmVX0Q2F9ERu94lehoT52LrclaZZYwTyuaXEHOudVaDXzBaYPYa+fcJFsUf3olSBEZAXwEfE9EpnZzaluH53OwhkvCHAWsCT1fCxza2QWqq6upqKiIPCorK3c4J2+bYY/o50RGpl6ZauVpwS8X1rF+a4C9J+Zw9anpFRWvrKyMfIeBLife9mrgRlW/E5E9VXWTiPQYbxaRE4FHgY4/aWVAfei5B+g0u7u8vJyVK1d2e/3IsEeMTdatHiPIdOHZt93c91wLOdnwwJVlZGelV1N13rx5zJs3DwARqenqvFgqBlwATFTV+SLyE2ALsElVW6J4+UXAhYADGCsiPwOqgfDyQ4VAbbS2bI9zm9keSiCgUVUY22qarGlBXVOA8263vh6/nVvC7mMzd7mHWJqsuwDLAVT1P8DC7k9vR1UPVNWZwMnA86p6H/AsEM6HnQS8HIMt2xBusmbntgJEXcHceMj04Bd/qWNzbYADds/lsuMzK6q6PbEI8gPAHYqOXkaHgE1niMiewAQR2aOLU5YA40XkbGCFqq6LwZZtCAvSnuMFiKrZGs7ScdhMlk4qs2RZCw8vc+F0CPf/qiwtaqv2hVj6kCuAa4CbgY3ACd2drKqrgdHb7fsGmBt6HgAWxHD/LgkPe9iyLUFGMxa51Ws1V4c4hmRMlkemse47HxfeaTVV/3T+YHYemTkJAF0RSy7rl8BZAGJ9g2ckyqhYiUxQzraKXEWTz2qaq6mNz6/M+X01TW7lp9OdnPvjzMhV7Ymeclk/AFoBL1YfcgNWLquEXvvDRBsYDZFxSHv0fciIIHONIFORG/7dwDuftTGq3M698zMnV7UnevKQx4X7diIyT1Ujg4AicnVCLYuBsCCD9hg8pNdk6aQqr37Yyh8fbsRmg39fXc7gopSoONov9JTL2jHQUiYiE0WkSESOBGYn1rToCQ97BGweQKPrQ5oma0pS0xjgzD/UoGoVq5oxOTWXHk8UsQR1/g78DpgGbALmJcSiXmATGw6bA0/QA/a2KPuQZgwy1QgElDP+UM13tQH23z2Xa08flGyT+p1YgjrVWGUgUy6oA9aMD0/Qg2R7qGvqOX3OeMjU46bFjby40kNpkY3FC8ozLnE8GmLJ1FlLe5WAYuA74PVEGNUb8ux51PsakKxWapp69pBVoTzWIUaQKcEL77byu0UNiMCiBeXsNDT9y3H0hlje9SGqugGsqVPAdYkxqXeE81nJ9lLTGIuHNE3WZLN+q5/Tf1+NKvxm7iAOr0jvynF9IRZBjhWRsaHnTqzZGr+Jsz29puOMj9rG7j2kJ+ChwddAlmRRmpP+lcrSGW+bctJvq6hrDnLkD/JYkGazOOJNLIKcj5U+B9a0qnPib07v6SjImh76kGHvOCS33Kx4lURUlZ/fXcvKz9sYM9TOv68qw2YbeP3GjsQiyPtU9RmIBHV+AnRZ4qO/6VhXp3ZrkGBQu/znmjHI1ODu/2vm/udacOQIj14/ZECNN3ZFj4IMzeS/FNhHRH4a2m3Dytx5JoG2xUQ4n9VZ4KUxaCWYd/UPNhHW5PPCu61csdCqqXr/L8uYOinzZv/3hh4FqarfishfgAOAN8O7gapEGhYr4XzWwsI2GoGapq4FWeUxEdZk8tmGNk69sYpgEH59ejEnHpSfbJNShm47UCJyDICqfo2Vxzom9BiLNek4ZQhn6+QXWDM+uou0mqSA5FHXFODYa6todCmzf+jkugE4+N8dPXnIjoWQzwO+pr1o8vZrfSSVcB/SmW8JsrtIq6mlkxw8bUFmX1/Fl5v87DUhhwd+ZYI429OtIFX1rx02b1bVjyES1PlBIg2LlbCHdDij8ZBGkP1NIKCceXMNb672MrLMzv/9bgj5eSbCvT3RBnUOAqaLSLgPKcBpwOEJtC0mwh4yJy8kyG6ydUyTtX9RVS6/p46ly90UOYX//H4oo4cMzEycnujxJ0pVvwWagACWEAUrqBOX2f7xIiLI3Cg8ZKjJaiqW9w+3PdrEX55sJicblv52CHuOz9wiVX0lqp8pVX1SRJ5SVQ3vE5F9E2dW7ETq6uR234f0BX3UemsRhLLcLstjGuLEP59v5qq/W9U+//mrcg7aa+CmxUVDNBUDwt6x46xtwQr4TEykcbEQzmW1ZVtVA7rK1qn2VqMo5bnlZNlMsymRLFnWwrm3WTVxbrugxAxvREFP38hjVXV9ZwdEJKUiImEPqaGqATVdeEgT0OkfnnzTxZk3WxONb5g7iF8cN7BzVKOlp4oBXYlRgEsSYlEvyc+yfn2Dobo61Q2de8j2Wjqm/5gonnvHzSk3VhMIwlWnFnPNHDPWGC0xxZ1FZKiILAC+JMWSy/PtliADNkuQW+q6EGQkj9V4yETw/95yc/z1Vfj8cMnsIn53lhFjLETViRKRQ4ELgO8DH2JVm0upIpn5oT6kR61VdV0epdkdpNC57W+OabImjsdfdzHn99X4A3DBrEJuu6BkwFSLixc9ekgR+TOwCHgQq+T/M6r6XVfN2WQRbrK6/G6GlVg5rJ15yfAY5BDTZI0ri15s4dSbLDHOP6GIu34+2IixF0QzDvkLYD+sKuTnAyUAIpJSg0nh2R6ugIthpSFB1u8oyCqTNhd37nmqibNuqSEYhGtPL+aP84xn7C3RjkOuB+4KifBYEbkcmIC1olVKEPaQbr+b3cMesrYzD2kEGS+CQeWa++q59ZEmAG4+t4QrTjLR1L4Q00CcqrZhrfOIiJybEIt6STiX1R1wM6zU+nXuzEOGm6zDHJ0uR2mIEm+bcvatNTzyqgu7DRZeVspZPypMtllpT69HxlX17/E0pK/YxIbT7sQdcDO4xFq4efs+ZEADVHmrAat8h6F31DUFOO43VbyxyktBnvDIdUM4Yh+TgRMPEpZuLyIzROSV0PNRIvIfEdkgItd0OOdyETldRC6Oxz3D2TolJVb63PaCrPXWEdAAJdmDyLWbGeq94aOv2tj3ws28scrLiFI7r90xzIgxjiRMkKq6nPY1JKcAxwJTgfkiMkREpgOlqroIKIlHbmx4LLKouHNBmlo6feORV1uYfslmvt7iZ+qkHFbcPZy9djY/bPEk0RPS2gBU9RlV9Yeqn68BGrDKSK4JnfdpaHsHqqurqaioiDwqKys7Ow1oH4ssLLIEuXU7QZrSHb3D51d+eW8dp91UQ6tXOePwfF67Y5iZQhUDlZWVke8w0OWshn79REVkDPCCqraJSBlQHzrkATqNspSXl7Ny5cqorh+OtDoLWoFBbK7zb3PcLEEXO19u8jHn99Ws/LyNLDvcfuFgLphVaIY1YmTevHnMm2cthyMiNV2d12+CDOW/zgb+ENpVjVVwGaAQqO3rPcJN1tw8DzYbVDUE8fmV7CzryxOZmGyarD2iqjz4kouf311LS6uy0xA7ixaUM22PgbUaVX/TnzUUTgP+oap+ERkKPAtMDh3bDXi+rzcIB3Va1c2IUjuqsKmm3UuaWjrRUVUf4NQbqznrlhpaWpUTDnTyfuUII8Z+IJFR1j2BCSKyh4j8EWstkFdFZA2wp6quADwichbQEAoC9Yn8ULaO2++O9G/Wb23vR5qkgO5RVRa/3MIeZ2/isdfdOB3C3y8v5aFflzOowBQx7g8S1mRV1dVY6XYAV4Ye259zYzzvGclnDbjYaUgW//vEy7dVHTykSQrokvVb/Vz851qee8eaLXPI3g4WXlbKuOEpNYcg48moMFnHBPOdhli/6Ou3tgsysgSdCepEcHuC3PJII396pAlPm1KcL/zp/MHMPbLABG6SQEYJMpJg7rc8JBDxkKrKFhPUiaCqPPaamysr6/i22mrWn3RQPreeV8KIsoz6WqQVGfXJh8ch3QE3e4YW/NwQEmSDr4G2YBsFWQWR8wYiqsrL73u47oF63v3MSjGcsnMOd1w0mOl7mqBNsskwQYabrO0eckOV9etvxiBhxccerr2/nuWrrMSJIYNs/PasEs46sgD7AFw+PBXJLEHaw0EdN2M6eEhVHbBjkKrKsg883PpIIy+/ZxUAKym0ccWJRVx8bJGpHp5iZJQgw+OQbr+LonwbJYU26puDbK0PDLgxSH9AeWK5mz892sj7a62maUGecOnxRVx2XDHFBUaIqUhGCTK/Q1AHYOLILN75rI3Pv/Wz1Tkw8ljrmwP8+8UW/vpkM+s2W/3nIYNs/Hx2EecfXUhJoRlPTGUyS5DhqgEBNwC7jM7mnc/a+OJbH1VjMns9j/e+8HLvM808vMxFq9cqML/zyCzmn1DMGYfn48gxHjEdyEhBuvyWICeNtga1P9/oo2lY5jVZXa1BHl/uYuEzzZGIKcChUx2cf3QhR+/vNMGaNCNDBWk1WXcZZQnyi299yB6ZscBOMKi8sdrLohdbeHy5i5ZWyxsOKrBx5hEFnHd0IZNGmeyadCWjBNmx8hzAxLAgN/oYlOZR1nXf+Vj0UguLXnLxzZb27KN9d83lnKMKOOmgfJwO0yxNdzJKkJHEAL+VjzlxVBY2G3y92c9wT/pVLG9yBVm63MW/XmzhzdXeyP7R5XZOO6yA0w8rYJfRxhtmEhklSIfNgQ0bnqAHf9CPIyeL743O5tPvGnAH3DhsDoqyipJtZrf4/MqLK1tZsszFkyvckQCN0yHMnu7k9MMLOGgvh1kKPEPJKEGKCPlZ+TT7m3EH3BTZipg6KYfPGq25z0McQ1IyYToYVN782MuSZS0sXe6mtsPqzzMm53LG4QUcNyN/h2URDJlHRgkSrH5ks78Zl99FUXYRe0/K5eFPLEGmUnNVVfngyzaWLHPx6GsuNla3z9vcfWw2Jx+cz8kH5ZvpTwOMjBNkfpYTvO2R1r0n5mDLDwkyBfJYv9joY8kyF0uWtfDFxvbgzJihdk46qIBTDs43S34PYDJQkOHkACuws/fEHLKLrJpCRbbkCPLrzT6eeMPNo6+5eO+L9vHC8kE2Tjwwn5MOzmf/3XJTsjlt6F8yVpBhD5mXa2P0To1UAc01g/vNji83+Vi63MUTb7i3EWGhU/jpdCcnH1zAwVMcZJmBe0MHMk+QoRkfLf6WyL6yYfVUAd+uT+zioZ9/a4lw6XIXH33la7fJIfxkfyfH/dDJj/bNIy/XBGcMnZNxgizoUFcnTO6gWmiB9z8qxNMWjFtep6ry6XofS5e7eWK5i4+/aRdhoVM4en8nx8/I57AKhxGhISoyTpCF2dYKTE2+5sg+F9YCO821g3lqhZuTDiro9fWDQeW9L9p46r9unnzTzWcb2kU4qMDGrAPymD0jn8P2ziM3xzRHDbGRcYIsyLLE1tyhyRqeCxl0lfKHxY0cPyM/pqTrNp/y2kcenl7h5pn/udlU0z5EMbjQxrHTncz+oZODp+SRk21EaOg9GSfIomwrE6c55CFbA600+hrJlmxGFQ3m4298/OnRRq48pfv+5OZaP6+87+H5d1p59m03TW6NHBtVbmfWAU5mHeDkwO87IpXRDdvi8/nYuHEjHo8n2aYkDYfDwahRo8jOjm48OeMEWRjxkJYgqzyh9SAd5dx4STnH/LqKa+5roLYpyMXHFjF6iB1Pm7KhKsAHa72881kbr37Qyuqvfdtcd89x2Rx9gJNjpjnZe2KOGaKIgo0bN1JYWMjYsWMH5OelqtTW1rJx40bGjRsX1WsyUJBWHzLsISNL0OUO5cf7ObnjosHM/1sdtz/WxO2PNWGzQTC443WcDmHGZAeHTXVw9P5Oxo8wGTOx4vF4BqwYwUrlLC0tpbq6OurXZJ4gQ0GdcB8yXNwqXLrj5z8tYp9dcrj1kSaWr/JQ3xwkOwuGD7YzeUIOe0/M5YeTczlgN4cJysSBgSrGMLG+/8wTZFY4ytoEdL6ex367OVh6g1WD1B9Q7DbzxUl37r77bkaNGkV1dTW77LILBx54YEyvDwaDnHfeeVxzzTWMHTs2MUZGQQYK0upDtmznIbuqFGAyZRJP1qHfxOU6/pfHdrr/3//+N/n5+fz0pz8FYNasWUyePJmSkpKor22z2Rg5cmRM9jz99NN8//vfZ8yYMTG9rjsSJkgRmQFcr6qHiIgNuBb4CrCr6r862xeP+0airKGgTruHTM9KAYaeeeaZZ7juuusi26WlpRx66KE893fmOeYAAAnySURBVNxzvPHGG0ycOJGvv/6aLVu24Ha7KS4uZtWqVRQUFDBmzBgCgQA77bQTAM8//zxPPvkkDz/8MB9++CFNTU2sWLGCSy65hLVr10a2L774Yh5//HHq6uqYO3du3N5LIle/Wi4ieaHNU4HNqvqgiCwUkWXAgdvvU9Vv+3rfSJTVF/KQXlOxPNl05dnihaqi2j4sVVZWxs9+9jPeeustHA4H48aN44EHHuCYY45h2bJlTJkyhWAwyDnnnMMZZ5zBkUceybRp03jnnXc48sgjqaur4+OPP+Yf//gHixcvpqysjL/+9a9s2LAhsn3PPfcwfvx4Zs6cGdf3kuh8rnBW9VHAmtDztcChXezrMwWRoE64D5netXQMPXPUUUfx9ttvR7bXr1/PMcccw6JFi3A6nQSDQTweDzNnzuS8884DICvL8kU33ngjq1ev5qGHHoq83mazEQgE2Lx5M21tbZSVlVFUVLTDdiLiDv2VYFkG1Ieee4BhXezbgerqaioqKiKPysrKbm+Ub3diw4Y70Io/6O8gSOMhM5W5c+fS1NTEU089xf3338/FF1/MyJEjmTRpEtOnT6e4uJiCggKOOOII3nrrLb788ktWrVoFwA033MCUKVPYaaedWLduHV9//TXr16/nq6++4uqrr+amm25i2bJlnH/++TtsT5gwgcWLF2/jnbuisrIy8h3G+u53ikRzsd4iIq+p6kwRWQz8WVXfEZGrgDqsJus2+1R1B7VVVFToypUrY7rvxGd3o9HXyMdHfMgeL+yFDRubZn2DXUzV7v5kzZo17Lrrrsk2I+ls/zmIyHuqWtHZuf3lIZ8FJoeeTwJe7mJfXAj3I9e5vgagPLfciNGQFiRMkCKyJzBBRPYAlgDjReRsYIWqrutiX1wIj0WubV4LmOaqIX1IZJR1NTC6w64F2x0PbL8vXoSzdT5v/hyAEXkjEnEbgyHuZOSs2bCHXNNkCXJUXmwDvgZDsshIQRbnFAPwSdMngPGQhvQhIwVZnlMKQF2bNaoyymkEmenceuutnHzyybjdbtxuN+effz7r1u0YlggGg5x77rl88803UV33hhtu4PXXX4+ztV2TcbmsYEVVO2I8ZHIZ8lR8ugxVx2zq8tg+++yD1+vF6bTWd5k+fTrjx4/f4bxYc1bHjBkT1ThjmPfeew+Px8O0adOifk1HBoQgRxpBDjhUlVmzZjFnzhwqKyt57LHHePPNNwkGgyxfvpy5c+fy4osvUltby0cffcTBBx/MkiVLKCgo4Ac/+AEiQltbGyLCihUruPPOO7nllltoa2tj9erVrFmzhuOPPx6bzbbN9hNPPEFBQQEVFRXk5ubGbHdmCtLRngiRZ3cw3DE8idYYuvNs8cJut2/jyUSEwYMHc+KJJ/Lxxx+zZs0aHn74YR566CHef/99wJqtcfzxx/PBBx8wadIkRowYwY033sgvf/lLRo8ezZw5c3j66aeZNm0a48aN4+233+bVV19l4cKFVFdXc8UVV5CXl7fN9hFHHMHYsWN7JUbI0D7kkA6J5BMKJmCTjHybhg6UlZWxZcuWyLbf375Mg81mIxgM8tVXXwHtc1/r6uqYOXMmv/jFL4D2/Nb58+djs9m45ZZbtrlGIBCgpaWFmpqaSD7r9tt9zW/NSA85oWBC5PmkgolJtMTQX+y6664UFRVx1113EQwGmTJlCg8++CBffPEF69evZ926dcyePZtLL72UxsZG1q5dy7777suBBx7IKaecws4778xnn32Gz+fj5ptv5pBDDqG4uJjPP/8cVWXr1q24XC6uueYabrvtNnbbbTeuuuoqmpqattnesGEDS5cuZb/99sPhcMT8PhKayxoPepPLCnDS/07j1arXWLTvAxwx7PAEWGboCZPLahFLLmtGekiAB/b5B+vd69m16HvJNsVgiJqM7Vw5s/KMGA1pR8YK0pAatLa2xjSOl0moKq2trTG9JmObrIbkM3z4cDZt2oTP5+v55AwlOzub4cNjGHYL1yNJ1cfUqVO1O+69995uj/c3xp7uMfaoAiu1i+972jdZeyrp0d8Ye7rH2NM9aS9IgyGTSPlxSBGpBtZ3c0oZUNNP5kSDsad7jD0wRlXLOzuQ8oI0GAYSpslqMKQQRpAGQwphBGkwpBBpLUgRuVxETheRi5NtC4CITBORLSKyWUR2SZINM0TkldBzm4hcLyJzROTMZNsT2j4l9BmtF5Hu15WPvy2FIvKYiKwTkb+F9p0jImeJyC9DC0AllaQb0FtEZDpQqqqLgBIR2TfZNgEzgeGqOlxVP0+GAaq6HNhhkSNgfxEZ3fUrE2+PWJMFx6vqMFUdo6oN/WzOfsBcYA/gEBHZB5ihqg8AW4ET+tmeHUhbQbLtYj2fhraThogMAY4F1onIYcm0he4XOUoGYXv2AE4UkU9EZEp/G6GqL6mqS1XdwMdYn8/a0OFPSPJ3CNJbkFEt1tNfqGqVqu4D/AS4u7+bY12Qap/RalX9PnA+8KAkadlqESkENgA+UujzgfQWZDXgDD0vBGqTaEsEVf0EuB/YseRZ/5Oqn9EbwOtA9Escx5fTgetIwc8nnQXZcbGe3YDnk2gL2/3at2E1o5NNwhY06g3bfUbfqmpdEmw4FnhSVZuBF4HdQ4eS/h2CNBakqq4APCJyFtAQCh4kk+NF5L8icjnwuqp6kmFEFIscJdOe+SLygohcCjyaBFsuBO4AnhaRVcAhwLsi8jNgOLC4v23aHpM6ZzCkEGnrIQ2GTMQI0mBIIYwgDYYUwgjSYEghjCANhhTCCNJgSCGMIDMQEXmpv2YuiEVSEw4yCTMOmaaIyF1YeZi7AEOAV4CDVPVQEbGraqAfbenX+2UyRpBpiojsqqprRGQu8D1VvSq8L9m2GXqPabKmKZ0JLyTQySLypohkichtIvIHEbkz1Iw9VkQeFpFfA4jID0TkVBF5IJTjuQMicoiIPCMi80XkXRGZLSJ/FJFXRCQvdPwJESkTkedF5MLQ+b9O9GeQiRhBZhiquir01481F3Kzql4KNADrVPUU4IjQ6VcAdcAy2pOst+d1rDzPO4DfANNV9UpgJTAldHywqtZglVP8L3AScEbc39wAwKztkZn4O/xtCj13dXge7u/tpqrPg1Xuo7MLqapfRFpUVUWk4zVagZzQ8Y73bVBVt4jkxO/tDByMhxzYqIjMDj0/KhQxjWFlGEO8MYJMY0JVCQ4A9hGRUaF947GmO03AalLuFTo2EdhPRHbFmpI1HrgYuENE3sCav7kf8IGI5HW4RwUwMnSNacBuIjISqxxHRej4qJCQdwH2FZGpWHWO2teWN0SFibIatkFETgSeTtZ8zoGO8ZCGCKFCXZ8YMSYP4yENhhTCeEiDIYUwgjQYUggjSIMhhTCCNBhSCCNIgyGF+P8Euw9eIuIz4wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 234x234 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAN8AAADaCAYAAAAmNJefAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAeGElEQVR4nO2de5hU1ZW3319VdTfNVbBBERDBeMEoyghqPhyjotFxMmocY2IGDUnUGIckaOLExIwzzugTRz+NMYlmWqP5Mho1xhtOosZoHi8Y8BZEEBVtBS8gzf3SNN1dZ31/nFNFNfSlqqnqqtOs93n6qbPPdVV1/Wrvvdbea8vMcByn90mU2wDH2VVx8TlOmXDxOU6ZcPE5Tplw8TlOmXDxOU6ZSJXbgFzq6upsn3326fR4Y2Mjw4cP7z2DusHt6Rq3B15++eVVZtbxQ82sYv4OP/xw64rujvc2bk/XuD1mwEvWyffdm52OUyZcfI5TJiqqz9cdF1xwQblNaIfbs41169axatUqWltbs/tuvvlmFi9eXDabtqc37Bk2bBh77LFHXufKKmhs5+TJk+2ll14qtxlOD1i8eDH77LMP/fr1Q1K5zSkL6XSat956iwkTJmT3SXrZzCZ3dH6saj6nsqmtrS23CWUlmUwWdH6f6PP9/qNHOe/Fr5fbDKeX+OlPf8qDDz5IfX09Tz/9dMHXB0HA+eefz3vvvVd84wqgT9R8v33/dzy64rFym+H0Ar/+9a8ZMGAAn/vc5wA49dRTmThxIkOHDs37HolEglGjRhX03NmzZ3PooYcyduzYgq7ripKJT9JQ4MfAZOA/zexeSVOB+wEDjjWzN4vxrLS1FeM2Tgx45JFHuOKKK7Ll3XffnRNOOIFHH32UZ599lv322493332XFStW0NTUxJAhQ1iwYAEDBw5k7NixpNNp9t57bwAee+wxHnroIe6++27mz5/Phg0bmDNnDt/61rdYsmRJtjxz5kx+97vfsWbNGmbMmFG091LKmm8E8FVgT+BnwL3AscBIK7KXJ6k+UYH3GVInvLfT92j70z4d7s8EqDPU1dXxta99jblz59KvXz/GjRvHHXfcwWmnncZTTz3FpEmTCIKA8847j3PPPZeTTz6ZqVOn8sILL3DyySezZs0aFi5cyG233cZdd91FXV0dP//5z1m2bFm2fMsttzB+/HiOPfbYnX5fuZTsW5up1SSNAW6SNAI4HThP0gVm9kSxnlWVcPFVEp0JpxiccsopzJs3j4kTJwKwdOlSZs2axaxZs7jooosIgoDm5maOPfZY9ttvP5YsWUIqFX4/rrrqKm655RbWr1+fvV8ikSCdTrN8+XJaWlqoq6tj8ODBO5RbWlqK/l5K6nCRNB64BrjIzFaa2RTgs8BPJe22/fmNjY1Mnjw5+1dfX5/Xc5IqzMvkxJcZM2awYcMGHn74YW6//XZmzpzJqFGj2H///Tn66KMZMmQIAwcO5KSTTmLu3Lm8/fbbLFiwAIArr7ySSZMmsffee9PQ0MC7777L0qVLeeedd/j+97/P1VdfzVNPPcWFF164Q3nfffflrrvuIp9GW319ffY7DNR1dl7J43ySEsAC4Dgza4z2/QvwJzN7Jffcnsb5Lnzpn3ngw4dYedqHxTDZ6QGLFy9uF9/aVdn+c+gqzlfyUIOZBcBcYE3O7hbg9aI9ZBcN6jrxppTezlnAQcAc4L+BMyRdTOjtfMrMmov2LFx8TvwopcPlxu12vQjcV4pnufScONInRrh4zefEkb4hPu/zOTGkT4jP2bW47rrr+OIXv0hTUxNNTU1ceOGFNDQ07HBeoWM4r7zyyh6NFe0pfSI67c3OXYspU6awdetW+vfvD8DRRx/N+PHjdziv0DGcY8eOzSuOl+Hll1+mubmZqVOn5n1NLi4+J/aYGaeeeirTp0+nvr6e++67j+eee44gCHjmmWeYMWMGf/zjH1m9ejWvvvoqxx9/PPfccw8DBw7kiCOOQBItLS1IYs6cOdx4441ce+21tLS08Nprr7F48WLOPPNMEolEu/IDDzzAwIEDmTx5MjU1NQXb3TfE532+imLEw4XNGOiIrgZMJJPJdjWUJIYNG8ZZZ53FwoULWbx4MXfffTe/+c1veOWVcBzH7NmzOfPMM/nrX//K/vvvz1577cVVV13FpZdeypgxY5g+fTqzZ89m6tSpjBs3jnnz5vHnP/+ZX/ziFzQ2NvLd736X2traduWTTjqJffbZp0fCg74iPq/5KopSjzSqq6tjxYoV2XJb27ZZLYlEgiAIeOedd4BtP8xr1qzJjvdsbW3Njve85JJLuP/++7n22ms58MADs/doaWlh06ZNrFq1Kju+c+3ate3KO/uj7+JzYseECRMYPHgwN910E0EQMGnSJO68807eeustli5dSkNDA2eccQazZs1i/fr1LFmyhCOPPJJPf/rTnH322XziE5/gjTfeoLW1lWuuuYZp06YxZMgQ3nzzTcyMjz/+mM2bN3P55Zdz/fXXc9BBB3HZZZexYcOGduVly5Zx//33c9RRR9GvX7+C30efyOFyyfxLuXPpb3xsZxnxsZ0hFTW2szfwms+JI31CfI4TR/qE+NzbWRlUUhemHBT6/vuE+DIEFpTbhF2Wfv36sXr16l1agM3NzVRVVeV9fp/wdmb+4YEFJNSnfk9iw+jRo/nggw9obGwstyllZeTIkXmf2yfEl7Z09jXVN95S7KiqqmLcuHHlNiNWlKyakDRU0q8kLZT0hWjfdySdI2lmMZ+VaW6mvdnpxIhSttEyqQM/A3xB0tHA7mb2P8BQSUcW60EBoeiCqAZ0nDhQMvGZ2ZtR/pYxwE3AKUBmiZjXo3I7epq9LFPzZUToOOUk3+xlJe0g5aQObATWAWujQ82EyXTbMXz4cHoywiUgdLikveZzKoALLrggu1ybpFWdnVdS16CZNQDTCBMptQL9o0ODgNXFe473+Zz40ZupA38DTIx2HwQUbWWTbLPTxefEiF5LHWhmL0o6TtJXgHVm9kyxnuUOFyeO9GbqQMzsqlI8K8gE2d3h4sSIPjEcxON8ThzpE+IzMuLzZqcTH/qE+Nzh4sSRviE+d7g4MaRviM8yQXav+Zz40EfE58PLnPjRN8TnDhcnhuQlPklfkVSXUz5V0jck5T9zsIRkhpftyrOonfiRb813I/CYpMxMhMuAW4FzSmJVgXjN58SRfMU3M8o9mBkR08/M2oD9SmNWYWxzuLj4nPiQ7/CyQyTdAhwgaTMwXNIIYO/SmZY/7nBx4ki+4vtX4GTgOmA54TShk4Bfl8iugtjW7HTxOfEhL/GZ2Vbg4Zxdb0V/FUFgASmlPMjuxIo+EWowM6oSKR9e5sSKHokv6u91dXyQpPskNUi6Odq3j6SPJK2QdEJPntsZAQFJpdzh4sSKvJqd0QTYzwM1hILdCzigi0uOAmYABvxV0hTgeGCsmbXujMEdEVhAlVLZXC6OEwfydbh8Cfg2sDkqH97VyWb2RGZb0kJgJWEKwW9LutTM7uroukz2sgy5iWi6IiAgmfA+n1MZ1NfX52be2+nsZbcCy8xsE3Tf7MwgaVB03VJgmqTRwO8lvWhmOzhsepy9LKr5vNnpVAL5Zi/rVHyS5hPWdFsBAd+QZNH2aPILsJ8DXJEpmNkHkq4GDqaI3tLQ4VKVDbY7ThzoquY708zeBpA0ysw+jLYFfKq7G0s6HXjIzDZK2gNYaeHgy1rgia6vLozQ4ZL0ms+JFZ2Kz8zelrQX8AngZEmZVH8Cvgk839m1ki4CLgVWS6oG7gX+UdJ9wBwz+6hYbwAyzc4qAlx8Tnzors+3nNDLeQBh8xNCD+adXV1kZjcDN2+3++qeGJgPgQXUJGs8zufEii7FFzUTfyJpIvAfVqFzdtKEI1x8eJkTJ/INsr8O7JYpFDtIvrMEFvgIFyd25BtqOAc4VVLm2z0K2L80JhWOWUBKVe5wcWJFvuI73czeyxQkVcRUogwBXvM58SNf8R0m6XIgGZUPAKaWxqTCCaKaz+fzOXEiX/F9knCUy1jgZcKhYhVDYEYq4XE+J17k63AZDqQJZ64fClxYMot6QID3+Zz4ka/4fkUowFsJZyxcXyqDeoJlvZ0VGQlxnA7pstkZxfdWmtn8nH0P0cVI7XIQ4DPZnfjRXc33TWClpBMl3SBpjJn9BTixF2zLm8AsFJ87XJwY0Z345ppZEM3PazCz96P9r5bYroIILCCV8ClFTrzoTny5nagNOdvJ7U8sJ4EPL3NiSHfiuyHKu7IcuCln+4ZesC1vwuFlVdm08Y4TB7qL851rZrO33ynp70tkT48wAlI+n8+JGV3WfB0JL9r/+66u6yR72XnRgiuXSipqysIwyF5F2h0uTowoVd7OTPaygwlzt0wBjjGzO4CPCecIFo1s9jJvdjoxoiTiM7MnzGyzmTUBC4FTgCXR4UVRuWhk8nZ6nM+JEwWLT9Jukh6WND13zb5Ozh0ELANagbXR7mZgz47Oz6QOzPzlpF/rEp/P51QS9fX12e8wRUgdmMtWwjwsd0r6JNBpajS2ZS/7IjA02jcIWN3RyT1JHWhmGGGQvcm2FHSt45SCfFMH5l3zSdozmsc3HHgbwMwWdXF+NnsZ8EfCmREABwGPdXZdoRiGEAklSXsCJSdG5Jsu/hFgBNAU7doTeKCL87fPXnYj8KKkr0XXXrMzRucSWIAQSSW82enEinybnR+a2T9kCt3NZO8ke1lJyOTsTCjhDhcnVuQrvnWSvgWsi8qTgItLY1JhBBaQUIIkSR9e5sSKfPt8xrZFUgRUl8acwgnMSCjs83mz04kT+YrvCsJabyyhd3NmySwqkIAAkSAhucPFiRX5iu+/COfwrQImAN8pmUUFYplmp9d8TszIt8/3upndlilIOr9E9hRMQECCRORwcfE58SHfmm+UpL+VdJik84BPl9KoQmjvcPFmpxMf8q35fgx8j9DLuYhwldqKYJvDxWs+J150tTjmcDNrBDCzDcDlOccOp5MhYr1NbrPTpxQ5caKrmu9GSdPNzCS9CqwH2ghDDRWzVkNgAco6XLzZ6cSHrhbH/Kec4rlm9ipkV6YdVWrD8sUsrPl8eJkTN7p0uEgaHA0lmyZp72h7b+Dfe8O4fAgIsn0+H+HixInuHC5GOIxsMnBIzr75nV7Ry4QOlwQJ93Y6MaO7lWk3AhdL2sPMPs7sl1Q5KxTlOFzMHS5OjMg31DBT0jSghbCpmgaOK5lVBbDN4ZLwms+JFfkG2d8mHF52nZkdAzxeOpMKI2g3vMwXSnHiQ77iOwY4Axgj6SYgr+Flko6R9GRO+WxJKyQtlbRbV9fmixGQQMhrPidm5Nvs/Amw1szel/R3wC35XGRmz0iqhWyIYryZdZg8qadkHC5JfDKtEy/yrfnqgQ8AzOzRzHaetESvBwNnSVokaVJHJ/Yke1mQjfMlfZUipyIodvay1whHvKwmHOEyBfhsIQaZ2WvAoZL+FrhT0sFm7TtpPcleFsb5Eh7ncyqGfLOX5Su+t4EXITtbtbmnhpnZs5KeJkwluKan98mQ8XYmfK0GJ2bk2+y8HxgMPBNdU/Cy0FGfL8P7ZrbTwoNMnM+zlznxI1/xXUPoLDFgHnBHPhdJOgTYV9LBwCWSHpc0C/htj6ztAMuOcHHxOfEi32bnQ8DKaHsP8lwWOurnjYmKC+lBjdkd7RwuLj4nRuRb820EZki6jbDm+1HpTCqMdg4XT6DkxIi8aj4ze1jSc8C+wA8IxVgRhCNcfFaDEz+6msl+GZ3n5zyEIq+x11MyqQOTSvqy0E6s6KrmSwPPQYdtuc0d7CsLQY7DxUMNTpzoaib7dZltScMIVxeaAxxuZkV3nPSUbXk7XXxOvMjX4fJL4Pgo1LBM0o0ltKkgts3nSxLgsxqc+JCv+B4C5kbbaSqkvwftHS4+sNqJE/nG+dqAUyVNAb4K3FU6kwojXJ/P43xO/Mg31HCXpAWE6zR83sxeKa1Z+RMQOVwk7/M5sSKvZqekLwOfNLPfAntGw8YqAvP1+ZyYkm+fbwrhuuqY2R8I08dXBO0dLi4+Jz7kK77XiFallfRZYGTJLCoQd7g4cSVfh8ujwK2SDiBco6+CvJ2WE+fzms+JD/mKbxPwK8IRL4eb2esls6hAcoeXubfTiROFBNmPKzTInpu9TFJC0r9Jmh45cIpCJnWgDy9z4kZJg+xm9gxQGxW/BCw3szuBT0ka0/mV+ZNJHZhQ0jNWO7EiX/G1AadJ+iHwAoUF2TPZy04BFkfbS4ATtj+xx9nLsnE+F59TfoqavSwnyH4gPQ+y1wFro+1mYIf8nT3KXpZ1uHgCJacyyDd7Wbc1n6TTJD0O3A2cAwzpoU2NQP9oexBFWtk2E+dzh4sTN7pbn+8rhMmTHgS+C/wv8BNJZ/bgWX8AJkbb+wN/6sE9diCwACSEMAzz9RqcmNBds/PvCEMLTZkdkn5NmC7+d93dfLvsZfcA/ynpq8AcM2voudnbSFualJJIyno8U8o3guI45aO7b+ncXOEBmFmzpLz6fNtlL4Mw/0tRSVuapJIA2X5fKu/wpeOUj+6+pQMldZTHpV8pjOkJuTVdVaKKVmujhpoyW+U43dOd+L4HzCBcnwHCJaFF6Lm8rpNrepU2ayOhsOtanaimNWgBBpTXKMfJg+7E9xkzm7P9TklHlMieggksaFfztQStZbbIcfKjS29nR8KL9r9QGnMKp83aSEY1X5VStLr4nJiQ7wiXiiVtAclszVdNS9DSzRWOUxnEX3xBW7bZWZ2oos3aymyR4+RH7MXXrtmZqPKaz4kNsRdf2gKSiUzNV+19Pic29AHxtWWD7Cml3NvpxIY+IL5weBnkxvkcp/KJvfjaLE0iEl9mhIvjxIHYi699zVflNZ8TG/qE+JLZmq/a+3xObOgT4ssOL/MRLk6M6HXxSZoqaYWk5VEe0J0id2C1x/mcOFGOiW/HAiOtSFPO00E6Z4RLtY9wcWJDr9Z8kkYApwMNkk4sxj23n8/nNZ8TF3q15jOzlcAUSZ8E7pd0lJmtyxzPpA7MkJsFqjN2nM/nfT6nvNTX1+emvdy51IHFxswWSbodGA9kU1L0LHVgkPV2VieqaQ62FtNUxymYoqUOLCaSlFNsAXZ6zYeWoIWaZJjpojZZy5b0lp29peP0Cr3t7TxT0vOSvgM8bWbNO3vDlqCFqijNTP9kfza3bd7ZWzpOr9Dbfb77gPuKec+tQQvVUc3XP1VLU1NTN1c4TmUQ+yB7S9BCTSIU34DkAJraXHxOPOgT4qtKVAHQP9WfJu/zOTGhT4ivJhHm6RyQ7O81nxMb+oT4qhOZPl9/mtLtxWdmvPjGVh5/cQvrN/lCKk7lEPu86qH4ombndt7O1jbjn65uZEFDC6PqUixoaOEbpw7i+18aQm1N7H93nJgT+29gKL6w2dk/Wduu5rvqznVsbAp49dZRPHn9nrz8i7148/1W/uaCj3hmwU5HORxnp4h9zbc1p9k5uGoIG1o3ALBuU5qfPbiB1345iprqMLa/9x4p7r1iBA/P2cz0qxv53NH9ufq8oQysjf1vkBNDYv+t25remg01DKseyrrW9aQtzb1/3syJh9eyV92Ovy+nTR3A/Fv3YuOWgMPO/4gnX3EPqdP7xFp8Zsbmts0MSIULo6QSKYZUDWZtyzp+9fgmvnzywE6vHTY4ye3/MpybvjmM865bxfn/dxUfrfLpSE7vEWvxNQfNVCWqsnE+gN2rd2feuyv4aFWazxxe2+09TjmyP/NvHcXQQQkOO/8jfnDbWlat97XdndITa/FtbmtiYKr9cmC71+zOPc9/wPQTB5JMqpMr2zNkYIJrvz6MV+r3Yu3GNAd++UMuuH4VCxp8bqBTOmIuvm1NzgzDq0fwxOvLmHFS503Ozhg9PMUtF9ex+FejGLtHin/4wcdMufAjbrhvPR96k9QpMrH2dm5s28iAVHuRBevGMGz0h+w3uqqTq7pn+G5JLp++G5edPYSnFzRzz1ObOez8jxg/MsVJU2o5+YhajjiwhlSeNavjdESsxbdyayMjaoa327fo1ZGMmfCXotw/mRTHT6rl+Em1/Pzbxl9e38qj85qY+ZPVNCxv42/2q+bICTUcOaGGg8dVM27PVN5NXceJtfg+3vIxe/QbkS2/sayFZW/sS9uh/42Z0X7u7s5RlRLHTOzHMRP78aPzwzjiC2+0MG/xVm5/dBOL3mth5bqAA0ZXMWFsFfuOSjF2RIq990gxdo8UY4ansvFGx4EyiC+aSLsSGGJmP9uZey1tWsro2tFAGHb4/q1r+fa0idxraRZuWMQhQw4ugsUds9vAJJ+ZXMtnJm/zqG7aErB4aSuvL22hYXkbzyxoZunHbSxbmeaj1W3sNiDBiN2SDB+aZPiQaHu3JHVDEgwZkGBQ//B18IDotb8Y3D/htWkfpVfFJ+loYHczu17Sv0o60szm5XPt/A8/ZO6Hb7OlbSst6Vaag608uPEP/GPt97h3xSYefLaJ5WvSXHzmEEZ88M9865WL+fb+M6mrriOhRJRkSYjOv8h1NcPYd+C+PX5/A2sTTDmwhikH1uxwLJ02Pl6bpnF9QOO6dPi3PmDlujTz325hQ5OxfnPA+s0BG6PXDU0BG7cYtdVi8IAE/WtEbbXoVy361YSvtdWiNrudoF+NqE5BKimqUiKVDLdTybD2TiUyx7btz5SFkACBBCJ6zW6rXZnc412c05cZPzLV4UCOfOjtmu8UYHG0/XpUzoqvq+xl9yx6nruX/z+SVk3CUiSppv/Wo5m/YhKLEk1M2q+aX15aR7/qBF8dN4OaZA0PffgI61rWYRiBBQR0Path2ojjuOSAWUV+yyHJpNirLsVeneay6pggMDY3h8LcstXY0mI0txhbtoavzS1BdjtzvLXNaE1DW9po2gpt6YDWtvAHoC0dDjhvC8LjmXJrm2GAWc5ftmzblUPb2p+/4/W7At/83GA+f2x7GeWbvUxFyl2bF5Lqgdlm9r+S/h441cy+njk+efJkKzR7meNUMpJeNrPJHR3r7ThfI9A/2h4ErO7l5ztOxdDb4vsDMDHaPgh4rJef7zgVQ6+Kz8zmAM2SvgKsM7NnevP5jlNJ9PrwMjO7yszuMLMbCr02pxNbEbg9XeP2dE2sxnZW2ofn9nSN29M1sRKf4/QlejXU0B2SGoGlXZxSB3S68EQZcHu6xu2BsWY2vKMDFSU+x9mV8Gan45QJF5/jlAkXn+OUidiIT9J3JJ0jaWa5bQGQNFXSCknLJR1QJhuOkfRktJ2Q9G+Spkv6cjns2d6mqHx29DktlbRbL9sySNJ9khok3RztO0/SVyRdKqms3/9YiC9nKtL/AEMlHVlum4BjgZFmNtLM3iyHAdEIocyEwi8By83sTuBTksaU26ZoJeLxZranmY01s3W9bM5RwAzgYGCapCnAMWZ2B/Ax8PletqcdsRAfHU9FKhuSRgCnAw2STiynLYTLa0P7z2gJcEJ5zAG22XQwcJakRZIm9bYRZvaEmW02syZgIeFntCQ6vIgyf4/iIr46YG203QzsWUZbMLOVZjYF+Czw095uTnVCRX1GAGb2mpkdClwI3Kli5vUoAEmDgGVAKxX0GcVFfBU5FcnMFgG3A+PLbQsV+hkBmNmzwNPA0DKZcA5wBRX2GcVFfBU1FWm7X/AWwqZwucn9jPYH/lRGW4AdPqf3zWxNGWw4HXjIzDYCfwQ+GR0q+/coFuKrwKlIZ0p6PkoG9bSZlWW9MUmHAPtKOhi4Bxgv6avAHDNrqACbLpH0uKRZwG/LYMtFwI+B2ZIWANOAFyV9DRgJ3NXbNrWzz4eXOU55iEXN5zh9ERef45QJF5/jlAkXn+OUCRef45QJF5/jlAkXX8yR9ERvjc5XSNmD930Fj/PFAEk3EY5JPAAYATwJHGdmJ0hKmlmvLSLf28/ry7j4YoCkCWa2WNIM4EAzuyyzr9y2OT3Hm50xoCORRWKcKOk5SSlJ10v6kaQbo6bo6ZLulvRDAElHSPqSpDui8Y47IGmapEckXSLpRUlnSPovSU9Kqo2OPyCpTtJjki6Kzv9hqT+DvoiLL8aY2YLotY1wLt9yM5sFrAMazOxs4KTo9O8Ca4Cn2Da4eHueJhzz+GPg34Gjzex7wEvApOj4MDNbRZiC73ngC8C5RX9zuwCxXhbaAaAt53VDtL05ZzvTPzvIzB6DMOVERzcyszZJm8zMJOXeYwtQHR3Pfe46M2uSVF28t7Pr4DXfroNJOiPaPiXyXI4sq0W7OC6+mBDNlv8/wBRJo6N94wmn7+xL2Cw8LDq2H3CUpAmE04zGAzOBH0t6lnD+4VHAXyXV5jxjMjAqusdU4CBJowjTQUyOjo+ORHsAcKSkwwnz6vR8Pe1dFPd27sJIOotwpeCyzEfc1fGabxclSgK1yIVXPrzmc5wy4TWf45QJF5/jlAkXn+OUCRef45QJF5/jlIn/D827DL4kJz5PAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 234x234 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOsAAADaCAYAAABHLLXdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2dd3xUVfr/38+UdEpCQu8IIlJEI7I0WUBhXRXsDRWRjXyVda27srju6lrX+gUFN7YV609RFGxrQQHZLyioS1eKdIRAAiGkTLnP7487GRJIIDOZSTLMeb9e9zX3nnvvuc/c5DOnPec5oqoYDIaGj6O+DTAYDDXDiNVgiBGMWA2GGMGI1WCIEYxYDYYYwYjVYIgRXPVtQG3IzMzUjh07Vns+Ly+PrKysujPoGBh7jo6xB5YtW7ZHVat+qKrG7Hbaaafp0TjW+brG2HN0jD2qwFKt5v/dVIMNhhjBiNVgiBFEY9jdMDs7W5cuXVrt+dzcXHJycurQoqNTn/bs2rWL/Pz8SmkHDhygUaNG9WJPVcSbPW63m8zMTJo2bRpME5Flqppd1fXHtVgNh1izZg3dunXD6XTWtykG7L6i0tJSNm3axEknnRRMP5pYTTU4jjBCbTiICMnJySHdE5diXf7r59i/aFN9m3FcMm3aNGbPnk1ubi7z588PKw/Lsvjd737Hpk2bImtchHn22WfZu3dvnT0vLsW6/6uNFHz8U32bcdwxc+ZMUlNTueCCC8jJyeHxxx+noKAg5HwcDgdt2rQJ6Z45c+awefPmkJ8VKpZl8cwzzwAwceJEmjVrFvVnlhPTThHh4nX5EXdc/k5Flblz53LPPfcEjzt37syCBQu47777+Pjjj1m4cCFdu3bl559/5pdffqG4uJgmTZqwfPly0tLS6NChA36/n/bt2wPwySef8N577/HGG2/www8/UFhYyKJFi7j55ptZt25d8HjSpEnMmjWL/Px8xo0bB8APP/zAhg0bmD17Ns8++yyvv/46rVq14ptvvmHSpElMmDCBUaNGsXDhQp588kl++9vf8vXXX/PKK69w7rnnsmjRItatW0e3bt3o2rVrMK+HH36Y2bNnM2jQIKZOncqUKVN48MEHGTVqFAMHDuS9996jZ8+e7N27l48++ojc3NyIvd+4FOvls15j3veD6tuMesM1YlOt8/B93vGItPLB+3Lcbjdut5ubbrqJxYsXk5SURKdOnXjppZcYPXo08+bNo2/fvliWxYQJE7jmmmuC//TffPMNo0aNIj8/n5UrV/L888/z2muvkZmZyTPPPMOWLVuCxzNmzKBz584MHTo0+OwXXniBadOm8atf/Yq9e/eyadMmcnJymDNnDvv27SMjI4ObbrqJd999l1atWnHeeeexfv162rZty759+9iwYQN9+/blp59+YsGCBcG8WrRoQdu2benTpw/t2rXD4XAwefJkZsyYQdOmTRk/fjxTpkzh3HPPJTMzk6KiItLS0mr9viFOxQpQapXUtwn1RlVCiwTnnHMOS5YsoXfv3gBs3ryZAQMG4Ha7GTduHDfeeCOWZVFaWsrQoUPp2rUr69atw+Wy/w3vv/9+ZsyYwf79+4N5OhwO/H4/O3fuxOPxkJmZSePGjY849ng8lWzJy8tj9+7dZGRk8M0337Bnzx4AMjMzadSoESICEPy87rrruOuuu5g8eTJ+vx8R4ZxzzqFPnz7cfvvtwbyqaqN26dKFrVu3UlxcTGJiIrt27WLw4MH06tUr+N0iQdzWBf1Y9W3Ccce4ceMoLCzk/fff58UXX2TSpEk0bdqU1NRUunXrxqBBg2jSpAlpaWmMHDmSxYsXs379epYvXw7AvffeS9++fWnfvj0bN27k559/ZvPmzWzYsIHJkyfzwAMPMG/ePCZOnHjEcZcuXXjttdeCJfvvf/97Ro8ezSOPPMLgwYNJT0/n9ddfp0ePHrhcLjZu3MiWLVvYvXs3mzdvpmPHjrRr144+ffrQu3dvVqxYESxtK+aVlZWF1+vl22+/ZePGjWzcuBGA0aNH065du+D+oEGDmD59OklJSRF7v3E3zmqpRcs57Ziz7FH633dllCxreKxZs6bSeJ6hYXD438WMs1bAUrtEVTUlqyG2iDux+tUPgEXs1igM8UkcitUuUS1TshpijLgTq1VesoopWQ2xRdyJ1Y8tVjW9wYYYI/7EWl6ymmpwVIgV3+DPP/+chx56iNLSUh5//PFqr9u6dSszZ86Mmh2hEHditQJDVWqqwREnlnyDTzjhBMrKykhKSuL222+vNr927dpxzTXXhGRLtIg7DyZTskaPhuQb/Le//Q2v10uzZs3YtWsXbdq0oUOHDgwYMIC33nqLwsJCAFauXMmzzz7LPffcw7Rp0ygqKuKWW24J5te6dWuWLVvGgAEDeOyxx5g7dy4PPvgg119/PQsWLOCzzz7jjjvuoFu3blF/vyGVrCLScELPhYnf8gHxLdYdIrXeqqKmvsFfffUVJ554IgUFBXTu3JmePXty//33s3DhQho3bszAgQMBGDVqFEOGDAn6Bo8ePZrRo0fzzDPPVDquyje4Y8eO/PrXv+a2227jwgsvpFmzZvzwww9Mnz6dESNGcPnllwPQs2dPioqKSElJYdSoUbjdbjZt2hTMr1u3bpSVlXHmmWeSlJREUVERZ5xxBrNmzSItLY3+/fsHvZiizTHFKiIzReQhEekK3FQHNkUVv98Wq9+KX7G2Vq31VhXlvsHllPsGX3bZZbzyyiukpKRU8g2+4YYbACr5Bq9YsYLXX389mEdNfYOlih8Ql8uFx+Ph7bffZsCAAagq+fn57Nu374jrN2/ezIoVK+jRoweqWmV+5513HpMnT2bw4MH4/X7S09MZP348vXr1Cv2PEAY1KVm3qOpkIAFoeqyLGzo+X6Bktfz1bMnxR0PyDa6Y7+LFi/n0009Zu3Yt48eP5+9//ztz585l27Zt/Pjjj2zfvp2ioiLmzp3LunXrWLp0aTC/devWsWnTJlSVSy65hIyMDBISErjiiiv44x//yNixY6kzl93qYpRWqNIMAZyB/auOdX1dbuHEdV2f95Nmvdda37/9yZDvjWVWr15d3yYYquDwvwtHiRt8zA4mVV1QYf+18n0RaYld2gJkq+q7kfwRiRa+QDXY8puS1RBbhNUbLCLvAmVAaSCpCxATYrXKq8FqxGqILcIdutmqqn8oPxCR0AbF6hFfoK3qN21WQ4wRrlg/FpFJQGHguC9wa2RMii7B3uA4HroxxCbhejA9CFiABLaYca4ItlXjeOjGEJuEK9Z3VXW6qr6sqi8DD0TSqGhi+W2RxvM4azR59NFHufzyyykuLqa4uJiJEydW6TQQqv/vvffeG7avcV0xe/ZsVq1aFbX8wy0RTxSRWdjVYAE6A2dGzKooYllmPms0Of300ykrKyMlJQWAQYMG0blz5yOuC9X/t0OHDiGNZy5btozS0tKgN1Q0mTp1KjfffDMXXHBBVJ8Trlg/ARYDPiARqBsXjgiggWqwcYqoG1SV888/n7Fjx5Kbm8vbb7/N119/jWVZLFiwgHHjxvHpp5+yd+9e/vvf/zJs2DDefPNN0tLS6NevHyKCx+NBRFi0aBFPPfUU//jHP/B4PKxYsYI1a9Zw8cUX43A4Kh2/++67pKWlkZ2dTWJiIlu2bGHhwoV8+eWX3HPPPcyfP5/mzZvz73//m7/+9a+MHTuWiy66iDfeeIO33nqLYcOGVevPfNFFF1XKa86cOfTv358PP/yQYcOGMX/+fFJSUrjxxhuZNm0aZ511Fjt27OCtt97i+eefJyEh4dgvrgrCFesQIENVp4nIcCBmFlEpL1k1jqvBC2VyrfMYrA9Vme50OiuVgCJCRkYGl156KStXrmTNmjW88cYbvP7663z33XeAPcPl4osv5vvvv6dbt260bt2a+++/nzvvvJN27doxduxY5syZw8CBA+nUqRNLlizhyy+/5NlnnyUvL4877riD5OTkSscjR46kY8eOJCYmAvD0009z7733MnLkSNxuN1999RUvvPACS5cuZfHixWRkZDBu3DjmzZvHwYMHjxrruGJeKSkptG3bln79+rF69WpUlVtvvZWbbrqJtWvXcs011/DMM88wYsQIunfvzo4dO+jYsWNY7zxcsa5R1WkAqvqFiCwF3gozrzqlfMgmnqvB1QktEmRmZvLLL78Ej8vdO8Gu+lqWxYYNG4BDMXvz8/ODcYS9Xm/QV/i2227jnXfe4R//+Afdu3cP5uHxeCgqKmLPnj1B3+CCgoJKx4f79vr9ftauXcspp5zCli1b2LlzZ9DeiteX23jZZZdVG+v4iSeeCOa1bdu2I95Bo0aNSElJYcWKFZx66qns3r2b7OxsBgwYUCkmcqiEK1aviPQDioFxQP7RL284aKCDyYrjkjWanHTSSTRu3JipU6diWRZ9+/bl1Vdf5aeffmLz5s1s3LiRCy+8kFtuuYX9+/ezbt06zjjjDM4880yuuOIKTjjhBNauXYvX6+Xhhx9m+PDhNGnShB9//BFVZdeuXRw8eJApU6bw+OOP06NHD+666y4KCwsrHW/ZsoV33nmH/v37k5SUxA033MC1117LgAEDuP/++xk6dCjPP/88IkKPHj3YsGEDmzdvZseOHWzcuJEhQ4YE/ZndbnfQnzknJ+eIvNLS0vj8889Zv349lmUxdOhQrrnmGnbt2gXANddcw7Bhwxg4cCCPPfZY+C+3Oj/Eo22AG/gfYDr2+GpaiPenA/8CVgKXBdJuB64GJlW47oi0ils4vsFfffOxZr3XWv957R0h3xvLGN/ghklEfYPLEZHrsUvQL1S1EJgR/k8EzYHxQEvgaRHZDjRT1cdF5C8ickbgB6FSmqouOVqmNcGyApEiTAeTIcaosVhV9QURaQyMEJF0wAv8R1XXh/pQVf0RQETaAVOBc4A1gdOrA8fuKtIqiTUvL4/s7EPBy3NycsjJyTn6s02b1dDAyM3NrbjaXGZ114XaZk3QwOwaEXEAA0TkTEBV9cVQMhKRzsDDQB6wDygP1lOKXeJKFWmVyMrKIuTlM/xGrIaGRcVCRkT2VHddqGJ9VUSeBT5W1TLg68AWMqq6MTDssxyYD6QETjUC9gJaRVqtKa8Gmw4mQ6wRqrvhH7CroleJyCQROV9EEsN9uNoLziwGXgd6B5J7YDtdfFRFWq0prwabtW4MsUZIYlXVH1V1Z6DKOxPbc2mbiIQ0cCcit4hIrohcC/xTVRcBpSJyHbBPVRdUlRbKM6rjkLuhCUUaDWLVN3j16tVMmjQJsGNBVcex4gxHk5CqwSIyFrt6eglwCjAXuEBVQ6oKq+pTVaQd8YaqSqstlhlnjSqx6hvco0cPioqKALj77ruPml9VcYbrglDbrDOAN4GtwBfAh6paN3EYI0TQ3dBUg+sEbSC+wS+++GIwPnFCQgKdOnWisLCQ66+/nn/+8580btwYgB07dnDnnXfy0ksvMWXKFFJTUzn33HOZPXs2aWlppKWl8dxzzzFp0iTGjx/P/PnzmTp1KhMmTOCjjz5i2bJljBkzhkGDBkX8XYYq1hHlY50i4gL+KCLnAP9PA+6HDZ3g0E0cr3XT/P3aB/bYPXp7lekN1Te4c+fOWJbFhAkT+P7779m3bx8zZ84kLS2N7t27c9ZZZzF//nxat26N2+3G6XQyZswYFi9ezKpVq+jatSsdO3akT58+FBUV0b17d04++WT2799P9+7d+frrryktLaVfv35s2rSpQYg1X0R6ApcGtl3Aq8DbkTYsWpgpctULLRI0VN9gOBSf+JVXXuGRRx7h5ZdfJj8/n6SkJByOyt03e/bsYd68eQwePJhNmzZVmd/VV1/NbbfdxqOPPsqSJUtITk7m0ksvZfv26LzfUMW6GvgBuyo8XFWj91ePEmqGbqJKQ/UNXr9+PatXrwZg1apVvPzyy+zZs4cLL7yQiRMnUlZWxu7du9m6dSvr1q3D4/Hw+eefk5SUxNatW7nooot45513yMzMZPv27ZSVlTFkyBA++OADWrRowciRIxkzZgyzZs3ir3/9a8hr9dSI6vwQq9qAi0O5PtpbOL7Bb7/3gma911rvH3N1yPfGMsY3uGESim9wqEM3syL/c1G3qOlgMsQo8bfkY3mbNY47mAyxSdyJ9VDJGn9OEfH4nRsyof49ai1WEckSkf+KyFgR6VDb/KJNvK7P6na7KS0tPfaFhjpBVdm7dy9JSUk1vicS8X4LgAdU9S0R6RiB/KJLeW9wnFWDMzMza+zaZ6gbkpKSaNu2bY2vD1usItJMVfcCyar6FoCqbgo3v7riUMC0+KoSNm3alKZNY37FzrgmrGqwiDwP3Bc4PEVExkXMoihjnCIMsUq4bdaVwPuB/UXAXyJjTvQJlqximQ4XQ0wRrliLgU4iMhCYBfwUOZOiS3lvsOUg2H41GGKBcNusq7GnyF2JHSni24hZFGWC1V+HHZZUnHE3emWIUUKdz9oOuAroh+0jvBtoDDwERH6aQRQIVoOdgvpNyWqIHUISq6puFZEvsJfLKJ9wrsDzkTYsWpS7GaoT8JtOJkPsEHI1WFW/5bBqb6DEjQmCQzZOTMlqiCnCarOKyBjgt9glrADdgOivrRcByte6UaegPhPo2xA7hNvBdDLwHNABWAacHTGLok1guEYdpmQ1xBbhdoVmAX6gPdAHmBgxi6KMP9hmVdNmNcQU4Yp1C3aY/+eA/kD9xGYMg+C6rA7TG2yILcKtBo9U1ScC+38SkaxIGRRtrAq9wWpKVkMMEa5Y94nI+9hLWgh2sO/so9/SMAh6MJneYEOMEa5YP8D2CS7vTu0fGXOiT9Af2IFpsxpiirDEqqqvHJa0OQK21AmHqsGC+oxYDbFD3DnGanDoRk2b1RBThDuf9dTDjgdHxpzoYwWdIkyb1RBbhO3ILyI/BJIdwAhixJHftFkNsUrcOfIH26wOM3RjiC3iz5HfuBsaYpS4c+QvnyJnGacIQ4wRd478liqiYkeKMEM3hhgi7hz5VS2cOFEHYKrBhhgiXLH+C+iOXboOJoYc+S21cOIw46yGmCPcavAkoExVXxOR2UCXmt4oIo2AF4HTgE9U9UYRmYBdUmcCj6uqVVVamLZWwrIUBw7UIUashpiiPuIG9wfGAT2B4SJyOjBEVV/CXkn9ksAyHJXSwrTzCJRDJaupBhtiiTqPG6yqn6nqQVUtxhb9OcC6wOlVgeOzq0g7gry8PLKzs4Nbbm5uTZ4fKFlNb7ChYZCbmxv8H8auSVZJuNXgmcAE7LjBC7HbriERqA5vAbxAYSC5FGiJbXDBYWlHkJWVxdKlS0N6rqUa7GAy46yGhkBOTg45OTkAiMie6q4Lt2Rtih0zeDGQD/w9jDyuBu4B8oCUQFoj7DmyVaVFBFXLlKyGmCRcsX6I3UHUGegItAjl5oBTxXuqegD4FHvcFqAH8Anw7yrSIoKqhQun3RtsxlkNMUS41eB3VfWB8gMRaV/TG0XkRuBOYK+IJABPAd+KyPXY1d2HVdUvIpXSwrTzCCy1cIgjMM5qxGqIHWosVhGZCqQFDjNEZBZ2W1OATsDQmuSjqtOB6TW47uma2hYKSqDNKqbNaogtQilZv8cepimr4lzMhHWxVHGKcYowxB41FmtgzLM6YiasS8UOJjPOaogl4jKsi1OcqJiS1RBbxJ1YLTTgwWSGbgyxRdzFYFK1cIrL7mAyQzeGGCLuYjBZ2B1MlvENNsQYcReDCVUcps1qiEFCrgYHYjA9hj2EU140XRpJo6KJFZh1g/ENNsQY4XowfQFsx559A7bL4RPVXt2AUFWcDieWKVkNMUa4Yl2sqneUH8RUdEPKh278xt3QEFOEK9ZUEXkC2IftbtgDuCxiVkUR24PJhTp8phpsiCnCFeuPwHLAhy3W0ohZFGUUC6c48ZuhG0OMEYoj/6fYiygrdpC00RVOtwEeibBtUcEKVIN9ps1qiDFC8Q2uGBv4GWC+qvoBROS0SBsWLSp2MJlxVkMsEW41+D/AH0TkZOA74NjBjxoIVqAajBh3Q0NsEa5v8JNAAnawtF3AHUe/vAERcIowQzeGWKM2Qzcvlx+IyP9EyJ6ooyguRyCsi6kGG2KIcMXaXkSuBg4A2dgdTDMiZlUUKe9gUjDjrIaYItxq8ENAEnAWsAOIsZLVZQKmGWKOcMWaAbQFEgELe5mLmMCioruhqQYbYodwxfoJ9sybd4D9wN0RsyjqlFeDTQeTIbYIt836jaoGBSoiFwc+U1X1YEQsixKWWjid9uRz02Y1xBLhihUR+V/sSPlOYKCI9MBeBvLKCNkWFTRQDVZTDTbEGOGK9SdsZ4jy//Z5gc8G7yOsqjidZpzVEHuEJVZVrW7u6vxa2FInWOW9wZiS1RBb1Cq6YSDES8ygloXlUFxONyp6RJv1YInFix8f4Jd8Xz1ZaDBUT21DkTaOiBV1hPos1Alucdlt1sPGWW9+Op+p7xRy9p278HhNqWtoWNRWrO9ExIo6Qn0W6gCnuPBjVWqz7i7w8/6iYuY/1YoOLVzkfnCgHi01GI4k3LjBiSLSH/iPiAwRkdsibFdUUJ+F5VISHG4ssSq1Wb/8oYQhvRNpkubgz2Ob8PR7hViWKV0NDYe4CphWXrImOBLxYVVqs371Qyln9kkCoP9JiTRKdvD5d6WcnZ1cX+YaDJWIq4BpwTarw4Vf/FjeQ2JdvtHLVSPsFS1FhKvPTuONL4qMWA0NhkgETAN7lfIGHzBNfXZvsNvhxo8fLbN7fVWVtVs8dG/vDl57yZkp3DdzHyVlFsmJcbckkKEBUpuAaSsAL3bAtKrWbG14+O2SNdGRiB8Ly2PPP9hV4MflFDKbOIOXtmrm4rRuCXy4uISLz0ytL4sNhiDhFhn/BE4BrgV6EgPtVahcsvrkUMm6douXE9u5j7j+yuGpvDGvQbs6G+KIcMX6BOAmxsK6lLdZEwLVYCsg1k2/+Ojc6shKxpiBqXz5fQn7imJmBqDhOCZcsS5W1YdV9WNVncWhtmuNCQz5fBHYd4jIX0VkrIhcW11abbF7g5UERyIWFlaZLcKteX7aNj9SrE3SHAw/NZnZC4uPOGcw1DXhirW9iFwtImNE5H6gX6gZqOoCoLyr9Upgp6q+Cvwq0LtcVVqtKK8GJzjc+PRQNXh7no92Wc4q77limKkKGxoGtQ3rcjawk/DDungCn+cAawL767DXe60qrRJ5eXlkZ2cHt9zco0dEPTTOmoBffcFq8NY8P22yqu5r+80ZyXy3zsPOvcZf2BAdcnNzg//DQGZ114U768YHPFd+LCItqN30uEygILBfCrSsJq0SWVlZLF26tMYPUZ8/UA1OwIc/WA3eluejXTViTU50cP6AZN6eX8zNF8aUK7QhRsjJySEnJwcAEdlT3XU1LllF5FMRkcD+OyIyL7B9CSyspb15QEpgvxH2pPaq0mpFsDdY3PjVj+XxoaoBsVZdDQa4fFgab8wrqu3jDYZaUWOxqurZgXVuAJYAI1R1GDAMmFBLOz4Cegf2uwGfV5NWK9SvWA5witOOw+SGwv1evH5Ib1T9qxjWN4nNv/hYv91bWxMMhrAJqc0qIu1E5C6gP3C3iNwD3AM8GOqDRaQX0EVEegJvAp1FZDywSFU3VpNWK9RnoaI4xIFLXFgpDrZtL6NtppNApaFKXE7hkjNTefNL09FkqD9CarOq6tbAcIsT+Lo8GXgh1Aer6gqgYg/vnw877z88rbaoz8JyKiJil6zJDrbv9FTbuVSRq85K4+oH8/jzlU1wOKoXtsEQLULuYFLVb4FvK6bFjiO/HxU7FKlLXFhJDvbsLqN1syO9lw7n9BMTaJTs4LNlpYw83Tj3G+qesHqDRWQM8FvsElaw25QDI2hXVLA7mMCBw17vJlnYs9dL8/SkY94rItxwXiP+ObfQiNVQL4Q7znoy9tDNx8DfgVciZlEUUZ+FFWizOsWFJjnYu8dDy/Tqe4IrcsXwVBauKGNbnhlzNdQ94Yo1C3vJjPZAH2BixCyKIuXuhk5x4BInVrKDgr1eWmbUTKxpyQ6uHJbKjDkm5Iuh7glXrP/Cdlp4Drtn+LFIGRRN1GO3WQW7N1iTYP8+Dy1qWLIC/OGixjz3wQH2F5mYw4a6JVyxrsFe4+YU7PHQ5hGzKIqox4/lCFSDHS401UVRflmNS1aAzq3d/OaMZGbMKYyipQbDkYQr1i+AW4HrgHHAeZEyKJpYHj9WsDfYCakODhaEVrIC/PHyJkybXUhxqSldDXVHfMVg8vhQh+IQCThFONHiMpo1Du036+SOCQzulcRT7xTy56uaRslag6EycRWDySqzS1YHDpzipDTBQfMEf1hODg9MSGfApJ1cf06jkEtmgyEcwq0G/wh8iL22zXzsRaoaPBqoBjsCThHFTiErMbyqbJfWbq4+O417Xy449sUGQwQIV6z7gf+q6nxVnQ+sFZH/EZFWEbQt4lge3yHfYIeTIqfQzBX+mOmUq5ow9/9K+L/VDX7xPMNxQLhifQr4RETOCRz/CXsY5+qIWBUl1OPHLxbOgFPEQRGaOsOPr5TeyMmTN2aQ89heyjwmer8huoQr1kmqms2hNm9SYEJ618iYFR2sgFjd4sYlLgqBxlI7b6SLhqTQrZ2be2eGHIbKYAiJcDuYeonIDOBEETkIZIlIc2yPpgaLXbL6cYoLpzjZr0IjajdHVUSYcUszTv+fHQzpnciofinHvslgCINwS9a/AJ8A1wOLgKuAkcDMCNkVFawyHz7x43I4cTlc7LOUFKv2fr7N053MnJzF9Y/uYetu4zdsiA7hlqy/Ac4FRhOYdaOqDX7WjeUJiFVcuMTJPh8k+yIT/eHMPkncenETxvxlF18+0YrGqWbJDUNkqc2sm1xibNaNz+PDoY7grJs9FrhLIteTe/uljel/UhKX3rcbr890OBkiS1zNuvH6PDgDXznRkcBuv+IoLIlY/iLC//4+g+QEYeyDeUawhogSqVk3j0fKoGji83lxBWr+Tk2kJAn8+ZGNtu9yCm/+pTllXuWy+3abIR1DxAhXrA5V/beqFgJ3AZsiZ1L08HnKbAd+wO91k9AUfAWlqBVZh/zEBOGte5rjdAi//fMu8gvNWjmG2hNSB5Bh778AAA0MSURBVFPAYf8qoJ+I/BBIdmBHyx8UYdsiTllpGS6xv7LXk0BqYx/OVDe+/WW40yMbqiXBLbz5lyzuer6AAb/fyXt/b0739gkRfYYhvgipZFXVrdjT45Zh+wR/FTi+POKWRQFvWQWxliWQmubF3SwFX34x3lWryDvjDIpffjliz3M6hUdvyOBPlzdh6K2/8OLHBzgUetlgCI1Q4wZPxS5BXwa2As9gj62eEXnTIo+3rAyXwxZrSbGbtDQvrmYp+PYWUzh5Mgn9+7P/ttvwbdoU0ede95tGfP5YS56eXcjFf8tj+x4zFmsInVDbrP2A6aq6DTtW8GxV7QRkR9yyKFBRrEVFbtJSfbibp1H6cx6eefNo9OCDpE6cSNEjj0T82T07JfB/T7emRwc3p+bs4NH/tx+P15SyhpoTqljfUtUyEbkUaA3cF0gvi6xZ0cHj8eBy2DGC9xe6SE7xkNi+KSVf/xd33744UlNJnTSJkjffxDoQ+aBoiQnC38ens2haKxYuL6XHddt58eMDZojHUCNCFWtzEbkDeBKYoKo+ETkJ2+2wwePzeHA67d7gffvdJCZ7SWzXhOLv1pMweDAAzlatSBw+nJLXXouaHSe0cTPngRb860+ZvDHvICdft50Z7xdSVGLCxBiqJ1Sx3o098Xyoqi4MzF/tCkyKuGURRi2LMvWQ5LQDeufvc+FO9JDYrgllG3aTMOhQZ3bqxIkcnD496p1Bg3ol8dmjLXnhj5nM+76Uzldu47bpe1mz2XPsmw1xR6hr3fiAuRWOdwJzIm1UNPAf9OJtBEnOJAoPWnjL3OD0kNC6EZ48Dwm/+lXw2oThw8HnwzN/PolDh0bdtsG9khjcK4ktu3w8O/cAI/+4i5YZTq4YnsplQ1NpnRmuC7fheCJuvM19BSX4mzlJdiaxdquXtumpeKwyXJ48fI50HOnpwWtFhNQ//IGDTz5Zpza2b+HiwQnp/Px6Wx76XTqrfvbSe8IOfjVpB/e/uo/v1pWZoZ84Jm5+sn35xfgyHCQ5k1iz2UP7zBRK/KXI1lVYkoi3oKSSY0Ty1Vdz4O678f30E65u3erUVqdTGH5qMsNPTWb6LcrXK0v5aEkJYx/Io7BYGdI7kSG9kxjcO4keHdxHXa7ScPwQR2ItwZfuJMmZxOrNXjplpbDFX4pv6bckt2lJ8cpfaDK4U/B6R0oKqbfeSuGf/0zGrFn1ZneCWxjWN5lhfZN5bGIGP+/0smB5KQuWl/H4W4UcKLEY3CuJ/j0SOb17Iqd1TSA1OW4qTHFF3IjVm1+Mr4mQ5Eji+3UeLjw/k/96CvAs2UjqKTdS9P2OSmIFSLv1VnZ3707Z55+TOGJEPVlemU6t3HRq5ebakY0A2Lrbx8IVpXyzpoxZCw6y8mcvXdu46Nc9kewTE+ndJYEeHdykGQHHPHEjVs+2/XgzhVRnGt/+WMYzJ7XmySW78W/aTtPbT2Pvu2toc3Pl+fOSnEyT556jYNw4spYtw9miRT1ZXz3tmru4cngaVw5PA6DMo/x3o4dv15axaGUpuR8cYM0WLy3TnZzcyc3JHRPo2dHNSR3cdG3jNqVwDBE3Yi3dVMDBPhYlhY3p2MJF52aZHPAeQIYMIn1Udzbc/CGW14/DXTlgd9LZZ5M6cSJ7R4yg2Wef4WzZsp6+Qc1ITBD6dU+kX/fEYJrfr2zY6WPVzx5WbvLy/qJiHn7Dy4YdPjIaOTihjYuubdyc0NZN1zYuTmjjpktrF8mJRsgNibgRa8naPIp/7WPdxmQuGJyKU5xklSSwZ9TptGyeRmqvluydvYqsS3sfcW/alCkgQl6vXiSNGUPSeeeRMHQojsaN6+GbhI7TKXRr66ZbWzcXDD6UblnKtjw/67Z7Wb/dy0/bfHy9opT1271s+sVH0zQH7Vu46NjCRfsWLjoctjVKMWKuSxq8WEXkdmA30ERVnw4nD6vUy4ElW9l8eyHffdmY5/6chn/PHk5eVcTKy9rRE2g3+Uw23DyXpiNOwJ1ROUKhiNBoyhSSr7iC0vffp+ipp/BeeSXO9u1xZ2fj6tYNZ5cuuLp0wdmqFY6sLCTp2Kup1zcOh9A+IMThp1aeImhZys58P5t3+diyy8fmXT5W/uzhw8Ul9vFuHy4HtEh30iLDSct056H9DCfN051kNHLQNM1Bepr9mZQgpue6FjRosYrIIKCZqj4uIn8RkTNUdUlN7l343r9Z9cO3lBWX4tmaT8mvy1i8Zwn3NTmNZl+8TcH06Yw+ZyhP7foXyRlZNDu9GbsnJLL8wj+RMbIbSR2b4kxJRJJciAiC2KPSvTshff6AWhatyqDjrgJ869dTOns2/o0b8f/yC1ZeHpKQgCMrC0d6OpKWhqSmHtrS0pDkZMTtBre78qfLdWS6ywUilTeHAyTwz19FesVNqkqvAZmB7TQntid460PnVJWDpRYFB5S9B5SCQj8F+y32brPYXuhnxQGLohKlqMT+PFBioaqkJTlolCKkJDlIdIu9JQhJbiHBBYkJDhIShESX4HaB0yE4nYFPh/0D43SAy8ER6cGvJXYUP/v4UHr5D8Xhr+DQtVV8Vtiv6Xs7Gh0G9KD5CeE1pRq0WIFzsNeCBVgdOA6KNS8vj+zsQxN+cnJyyMnJAeCjL2byn9bLcbkEd3dwu7xc931jLl2/gNI1LpIuvphLbroJa/s7vLvtffZ596GDFF+fUrx7V2Id9KCFekQUiYpOCeenjuS28fdxOKqKHjiAlZeHVVCAHjyIFhXZnwcPYhUVoSUl4PWiXi9aUoIWFgaP8XrB5wvuq88HqkdsallVphNI1yrSiLBTRVNVarSOXgJYCj6v4ssHn1+xLDvNstT+9NuffksPmcthXw8NfoVDaVRKK0fQwDk7sfzUEW9Aq0mvkJ9UeTZ0Ntz9N5qf8NtKabm5ueTm5pYfZlZ3rzRkjxgRyQXmqOoHIvJb4HxVvaH8fHZ2ti5durT+DDQYIoyILAusdnEEDb2HIA8ob0A2AvbWoy0GQ73S0MX6EVDePdsDexUAgyEuadBiVdVFQKmIXAfsU9UF9W2TwVBfNGixAqjq/ar6kqo+Eeq9FRrtDQJjz9Ex9hydBi/W2tDQXrax5+gYe47OcS1Wg+F4okEP3RwLEckDNh/lkkxgTx2ZUxOMPUfH2AMdVDWrqhMxLVaDIZ4w1WCDIUYwYjUYYgQjVoMhRjhuxSoit4vI1SLSIGIai8hAEflFRHaKyIn1ZMMQEfkisO8Qkb+KyFgRuba+7QkcXxF4R5tFpEbzAyJoSyMReVtENorI9EDaBBG5TkTuFJF610q9GxANKkytewVIF5GGsHDWUKCVqrZS1R/rw4CAB1j5xNUrgZ2q+irwq8BynvVmj9jz1zqraktV7aCq++rYnP7AOKAnMFxETgeGqOpLwC7gkjq25wiOS7FS9dS6ekNEmgNjgI0iclZ92gKUh/uv+I7WYa+xWx+U29MTuFREVolI37o2QlU/U9WDqloMrMR+P+sCp1dRz/9DcPyKNRMoCOyXAvUaOElVd6vq6cC5wLS6ruJVQ0N7RytUtQ8wEXhV6imkhIg0ArYAXhrQ+4HjV6wNcmqdqq4CXgQ617ctNNx3tBB7oe70Y10bJa4G7qEBvp/jVawNamrdYaWEB7tqXt9UfEfdgM/r0ZbD39FWVc2vBxvGAO+p6gHgU+DkwKl6/x+C41SsDXBq3cUi8p9A8Lf5qlpaH0aISC+gi4j0BN4EOovIeGCRqm6sZ3tuE5F/i8gtwFv1YMuN2EuZzhGR5cBw4FsRuR5oBURvDdAaYtwNDYYY4bgsWQ2G4xEjVoMhRjBiNRhiBCNWgyFGMGI1GGIEI1aDIUYwYo0zROSzuppBIjb16mxxPGHGWY9DRGQqtl/riUBz4Avg16o6QkScquqvQ1vq9HnHM0asxyEicpKqrhGRcUB3Vb2rPK2+bTOEj6kGH4dUJcqAeHuLyNci4hKRx0XkIRF5KlA1HiMib4jI3QAi0k9ErhSRlwI+s0cgIsNFZK6I3CYi34rIhSLyiIh8ISLJgfPvikimiHwiIjcGrr872u/geMSINY5Q1eWBTx/2XNadqnoLsA/YqKpXACMDl98B5APzOOTQfjjzsf1mnwT+BgxS1T8BS4G+gfMZqroHO6Tnf4DLgGsi/uXigIa+Pqsh8vgqfBYG9g9W2C9vX/ZQ1U/ADgFTVUaq6hORIlVVEamYRwmQEDhf8bn7VLVYRBIi93XiB1OyGqpDReTCwP45gZ7dVvVqUZxjxHqcEohGMQA4XUTaBtI6Y09J64JdTT0lcK4r0F9ETsKeNtcZmAQ8KSILseff9ge+F5HkCs/IBtoE8hgI9BCRNtghWrID59sGRH4icIaInIYdF6tLXbyH4wnTG2yoMSJyKfZK9PUyHzfeMSWroUYEgr6tMkKtP0zJajDECKZkNRhiBCNWgyFGMGI1GGIEI1aDIUYwYjUYYoT/DwVwXcZF+F2SAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 234x234 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rcParams\n",
    "\n",
    "# Reset initial conditions and propogate overshoot trajectory\n",
    "vehicle.setInitialState(180.0,0.0,0.0,12.0,0.0,overShootLimit,0.0,0.0)\n",
    "vehicle.propogateEntry (2400.0,0.1,180.0)\n",
    "\n",
    "# Extract and save variables to plot\n",
    "t_min_os         = vehicle.t_minc\n",
    "h_km_os          = vehicle.h_kmc\n",
    "acc_net_g_os     = vehicle.acc_net_g\n",
    "q_stag_con_os    = vehicle.q_stag_con\n",
    "q_stag_rad_os    = vehicle.q_stag_rad\n",
    "\n",
    "# Reset initial conditions and propogate undershoot trajectory\n",
    "vehicle.setInitialState(180.0,0.0,0.0,12.0,0.0,underShootLimit,0.0,0.0)\n",
    "vehicle.propogateEntry (2400.0,0.1,0.0)\n",
    "\n",
    "# Extract and save variable to plot\n",
    "t_min_us         = vehicle.t_minc\n",
    "h_km_us          = vehicle.h_kmc\n",
    "acc_net_g_us     = vehicle.acc_net_g\n",
    "q_stag_con_us    = vehicle.q_stag_con\n",
    "q_stag_rad_us    = vehicle.q_stag_rad\n",
    "\n",
    "'''\n",
    "Create fig #1 - altitude history of aerocapture maneuver\n",
    "'''\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_size_inches([3.25,3.25])\n",
    "plt.rc('font',family='Times New Roman')\n",
    "params = {'mathtext.default': 'regular' }          \n",
    "plt.rcParams.update(params)\n",
    "plt.plot(t_min_os , h_km_os, linestyle='solid' , color='xkcd:blue',linewidth=2.0,  label='Overshoot')\n",
    "plt.plot(t_min_us , h_km_us, linestyle='solid' , color='xkcd:green',linewidth=2.0,  label='Undershoot')\n",
    "\n",
    "plt.xlabel('Time, min',fontsize=10)\n",
    "plt.ylabel(\"Altitude, km\",fontsize=10)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.tick_params(direction='in')\n",
    "ax.yaxis.set_ticks_position('both')\n",
    "ax.xaxis.set_ticks_position('both')\n",
    "plt.tick_params(direction='in')\n",
    "plt.tick_params(axis='x',labelsize=10)\n",
    "plt.tick_params(axis='y',labelsize=10)\n",
    "\n",
    "plt.legend(loc='lower right', fontsize=8)\n",
    "\n",
    "\n",
    "plt.savefig('../plots/craig-lyne-altitude.png',bbox_inches='tight')\n",
    "plt.savefig('../plots/craig-lyne-altitude.pdf', dpi=300,bbox_inches='tight')\n",
    "plt.savefig('../plots/craig-lyne-altitude.eps', dpi=300,bbox_inches='tight')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_size_inches([3.25,3.25])\n",
    "plt.rc('font',family='Times New Roman')\n",
    "plt.plot(t_min_os , acc_net_g_os, linestyle='solid' , color='xkcd:blue',linewidth=1.0,  label='Overshoot')\n",
    "plt.plot(t_min_us , acc_net_g_us, linestyle='solid' , color='xkcd:green',linewidth=1.0,  label='Undershoot')\n",
    "\n",
    "plt.xlabel('Time, min',fontsize=10)\n",
    "plt.ylabel(\"Deceleration, Earth g\",fontsize=10)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.tick_params(direction='in')\n",
    "ax.yaxis.set_ticks_position('both')\n",
    "ax.xaxis.set_ticks_position('both')\n",
    "plt.tick_params(direction='in')\n",
    "plt.tick_params(axis='x',labelsize=10)\n",
    "plt.tick_params(axis='y',labelsize=10)\n",
    "\n",
    "plt.legend(loc='upper right', fontsize=8)\n",
    "\n",
    "\n",
    "plt.savefig('../plots/craig-lyne-deceleration.png',bbox_inches='tight')\n",
    "plt.savefig('../plots/craig-lyne-deceleration.pdf', dpi=300,bbox_inches='tight')\n",
    "plt.savefig('../plots/craig-lyne-deceleration.eps', dpi=300,bbox_inches='tight')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_size_inches([3.25,3.25])\n",
    "plt.rc('font',family='Times New Roman')\n",
    "plt.plot(t_min_os , q_stag_con_os, linestyle='solid' , color='xkcd:blue',linewidth=1.0,  label='Overshoot convective')\n",
    "plt.plot(t_min_os , q_stag_rad_os, linestyle='solid' , color='xkcd:red',linewidth=1.0,  label='Overshoot radiative')\n",
    "plt.plot(t_min_us , q_stag_con_us, linestyle='solid' , color='xkcd:magenta',linewidth=1.0,  label='Undershoot convective')\n",
    "plt.plot(t_min_us , q_stag_rad_us, linestyle='solid' , color='xkcd:green',linewidth=1.0,  label='Undershoot radiative')\n",
    "\n",
    "plt.xlabel('Time, min',fontsize=10)\n",
    "plt.ylabel(\"Stagnation-point heat rate, \"+r'$W/cm^2$',fontsize=10)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.tick_params(direction='in')\n",
    "ax.yaxis.set_ticks_position('both')\n",
    "ax.xaxis.set_ticks_position('both')\n",
    "plt.tick_params(direction='in')\n",
    "plt.tick_params(axis='x',labelsize=10)\n",
    "plt.tick_params(axis='y',labelsize=10)\n",
    "\n",
    "plt.legend(loc='upper right', fontsize=8)\n",
    "\n",
    "\n",
    "plt.savefig('../plots/craig-lyne-heating.png',bbox_inches='tight')\n",
    "plt.savefig('../plots/craig-lyne-heating.pdf', dpi=300,bbox_inches='tight')\n",
    "plt.savefig('../plots/craig-lyne-heating.eps', dpi=300,bbox_inches='tight')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plots are now saved in plots/craig-lyne-* and should match with the results from the paper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Congratulations!** You have simulated an aerocapture vehicle flying in the Venusian atmosphere. In the next tutorial, we will compute the corridor widths for drag modulation aerocapture at Venus."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}