{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Figure 2 - Optimization pipeline\n",
    "\n",
    "Create the figure panels describing the model simplification pipeline."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import FormatStrFormatter\n",
    "from matplotlib.patches import Rectangle\n",
    "from PySONIC.utils import logger, rescale, si_format\n",
    "from PySONIC.plt import CompTimeSeries, cm2inch\n",
    "from PySONIC.constants import NPC_DENSE\n",
    "from PySONIC.neurons import getPointNeuron\n",
    "from PySONIC.core import BilayerSonophore, NeuronalBilayerSonophore, PulsedProtocol, AcousticDrive\n",
    "from utils import saveFigsAsPDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Plot parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "figindex = 2\n",
    "fs = 12\n",
    "lw = 2\n",
    "ps = 15\n",
    "figs = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Simulation parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "pneuron = getPointNeuron('RS')\n",
    "a = 32e-9  # m\n",
    "drive = AcousticDrive(\n",
    "    500e3,  # Hz\n",
    "    100e3)  # Pa\n",
    "Qm = -71.9e-5  # C/cm2\n",
    "bls = BilayerSonophore(a, pneuron.Cm0, pneuron.Qm0)\n",
    "nbls = NeuronalBilayerSonophore(a, pneuron)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Panel A: approximation of intermolecular pressure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAL8AAAC/CAYAAACv6g0GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAVnUlEQVR4nO3de3RU1fXA8e+emUxIwlMMCgoSEKEKiCSAiCJg8QFVKirWxw+1vrAv0PrGtlq6tPhTCz/9SWu11dpa1GK10FZ8oFQrVdCE8AgoiDyNEF4BEkKS2f3jDhggGSZzJzmB2Z+17kpu7s3ZB9fO8dwz554jqooxqSjgugLGuGLJb1KWJb9JWZb8JmVZ8puUlfTkFxEVERtCMk1eQ7b8Wtex7KN/KCI6/8Ef1nmPHXYk6aiTo26PeF804ia8MThKfglEw9oHbMYhN8kfDHrfWPIbh9wkv3hhNVLtIrwxgKPk11Ytuf9s2NKjs4vwxgCuHnhbtuSBobDl5M5OwhsDrro91RE6lILsKnMR3hjAUfKHNm9l/WPQZeb7LsIbAzh+4LXRHuOS03F+tQ+5jENux/kj1vIbdxx3e6zlT1Vbt25FRGjevDmZmZl06NCBKVOmNGod3Ax1ZmXx43Oh+NQTnYQ37hUUFJCdnc3OnTspKytj2rRp3Hrrraxbt67R6uCm5c/I4LEzYFOPji7CmyagoKCAfv367TsfMGAAAHv27Gm0OrhJ/ohyUgmEt+90Ed40Afn5+fTv3x+Abdu2MXHiRHJzc8nJyWm0Okiyly7Z+yJLrHK/3PAp7Y/rzrwfjWbg1BlJjW9qN+H1CRQUFzRojD7H9mHK+fH123v37s3nn39OOBymTZs2DB06lEmTJtG+fftkV0vquhBKdqR4BALR0Z7Y7xqYI1RFRQVFRUWsWrWK448/3lk9nCQ/4v0xqg11Npp4W+TGsHjxYrKysupM/L59+zJ48GBmz57Nvffey/vvv88777zDr371K0aOHJm0ejhJfpvPn9ry8/M55ZRTar1WUlLC9u3bmTRpEmPHjuXqq69m/vz5FBYWMn369KQmv5MH3r3dHonYOH8qKigooGfPnrVeKyws5IorrqBFixaUlJQwZswYsrKyKCkpoVOnTkmth5uWPy3MTd+CIQN7MNBFBYxTTzzxRJ3XCgsLOe200wBYuHAhubm5+77fOxyaLG6GOkMhfpsHm7od5yK8acIWLVpEnz59AC/h9/4hLFq0iN69eyc1lpOhzq1lWxg+oS03X3g/N174s6TGN+YAdQ51uunzKyz4LZzyt3kuwhsD2NIlJoXZyywmZTkd6sSGOo1Djlp+IQLW8hun3CQ/wpWXQNGQ2j/oMKYxOGv5X+wFG7se6yK8MYCrPr8EGPwFtFpf4iK8MYDDbs8bz8OpM+e7CG8M4PKBV0DsgTelPfTQQ+Tk5PD6668zaNCgRo/v5mUWCXivsdhQZ0rLz89nzpw5XHPNNZxzzjmNHv+Qc3tEJBsYAXQHqoFlwExVLa3j/kPO7VFVdqUHKLx4IGe8+EGCVTcmLvWf2yMi6SLyv8BCYDSQCTQDLgaWiMhkEclIqDYiqGAtv3EqVp//r3iJn6Oqo1R1gqreoaqXAl2ApdF7EjLmMvj4gj6J/ro5zDX1RavGqOofVbXiwAuqWqmqzwGXJBr4jZMCbOycneivm8Nck160SlX3LaojIkeJSEcR6SQiOSIyPHrPrkQDn/8ZZK8sTvTXzWGuKSxadcjRHhH5OXBP9LQKCON1eXr5Cfz8jAjLd+TDjX5KMYerprBoVTxDnWOBTsBjwB3AUMD3K/TRISG/xZj6GDLk4J+NGQPf+x6UlcGIEQdfv/Za7ygpgUsvPfj6LbfA5ZfD2rXQMf7lJwsKCnjllVeYOnXqvkWrZs6ciUidgzNJF0/yb1TVL0WkCDhVVZ8Xkbv9Bo5I7OFQc+Q6nBatqhSRrsBy4CwRmY035OmLCogtUd643n237muZmbGvH3107Ov1aPUPp0WrHgKeAi4CJgHXALP8BvY+4fVbijkcNZVFqw6Z/Ko6i2iyi0gfoBtQ6DfwlVc1Y/iAvrZuTwpq8otWicjRwDS8aQ1zgHtVtQzvgy/fPuoc4tTjWiejKHOYORwWrfotsBq4CzgGmJzMwBcuqaLT0g3JLNIcAZrEolUiskhVe0W/zwA+VNVDRo9nYhvA+lYB1p7+DU6fvaTelTamHhJatGrfR22qWo43ozNpVLBxfuNUrOQ/8C8mqZmqIoitz28cijXa005EbqvrXFUf8xM4ImA7sxiXYiX/m+w/f6fmue+sVcFafuNUncmvqtcBiEhQVffr74uI79lHY69txVkn59Hfb0HGJCieF9ifrXkiIjcAH/sNvLx9GpuPaeG3GGMSFk/yp4vIVBHJFpHXgAnAeX4DX1xQQff8tX6LMSZh8ST/FUBHYCXey+u5qup7wZ3b39jJoDeW+S3GmITFmt4wusbpy8DpeA+6I0UEVX3FT2AVsXV7jFOxRnt+eMD5cmBA9FDAV/JHbLTHOBZrtGdoQwb2Wn6b02zcibVuz6si0jfG9X4i8rdEA1cHrOU3bsXq9nwfeCq6YtssYAXeH0tX4AJgKzAu0cA3X9+O3M4DqPOvy5gGFqvbsx7v4bY/cBneqI/i9f3Hq+qHfgJvaBumS5uEFnwzJinieZPrI+CjZAe+6ONddFm7ylv80BgHnKzSDDD2ve2ktVrhKrwxbtbnB4gEBKm2B17jjtPkD9gqzcahuJJfRC4VkV+ISKaIXJGMwBF7mcU4dsjkj67OdgswBsgAfiYiP/EbOBKwbYmMW/G0/N/B25lll6puxpvjc6XfwLd/93h+eYet2mPciSf5K2uu0a+q24BKv4F3tkhne3Nng03GxDXUuVZERgIqIunA7Xjr+fgy8uNSOulKuNxvScYkJp7k/wHwPNAb2AX8hyR0e85bsI1OW7f6LcaYhMWT/P1U9RwRyQSCqrojGYE1IEgkqUsBGVMv8fT5HwRQ1bJkJT6ABgI22mOciqflXyQiE4H3gH37dKnqJ34Ca0AI2Di/cSie5N/79tYNNX6meNuRJiwSCNiHXMapeGZ1NsgOYY9c34NIdTVvNkThxsQhnt0Yb6vt536XK6xOD1NeVe6nCGN8iafbU3PJwjBwNvC238BD52+i3erN8F2/JRmTmHi6PdfVPBeRDsAzfgPnLt5Cbv5Gv8UYk7B6T2lW1Q1AZ7+BbajTuFbfPr8AeYDvJlsDARvqNE7Vt8+vwBq8ndj9CQYI2rssxqF69flFJAwcq6rr/AbWYMC2JTJOxfMyy8Ui8riItAA+BRaKyHi/gf9w7WkMeLib32KMSVg8D7z34O3AfgkwDzgB+B+/gYOBENVqE9uMO/Ekv6jqIuCbwD9VtTTO34up/0fruX/6V36LMSZh8SRxRETG4G1I8YaIjAB8P6p2XbmVMf9J2iRRY+otnuT/MXATcK+qFgMTgR/5jxwgaM+7xqF4Rnvex+vy7B3tuUJV1/gNrMEAAcUb8ZE6N8k2psE4G+0hEPS+2sJVxhFnoz3V6WmUpgtU24iPccPZaM+cS3Pp/EBrCIf9FmVMQpyN9gQDQRvnN07VZ7RnYjJHe3p+tJpnXtgF5fZCi3Ej7tEeEWkdPR+UjMDHbNjOhYuqoaICMmyHFtP44hnt6S4iS4ElInKciBSJSA/fkYPR0R574DWOxNPteRwYD2yM7tP1ON7ojy8ajIa2oU7jSDzJ31ZV9y2yoKpPAi39BpaAtfzGrXiSX0WkGd6LLIjIsUDQb+DKrGasbwER24jaOBJP8k8DZgPtROQhvIVqn/QbeOl5fTn+xxA5pp3fooxJSDyjPc+IyGfASCANuLFmNyhRAfH+7qoj1YQCtk6/aXzxvMD+tqqeA/wrmYE7LVrLrD9B5KpVcKL/wSNj6iuebk9rEclKduDm28sY+RlEtm5JdtHGxCWe/sYuYLWIFLL/Ks0X+Yqc5s3piVTs9lWMMYmKJ/l9r85WG03zQmvlnoYo3phDipn8ItIT2AF8GP2AK3miszkju63lN27U2ecXkevwHnLvwnuB5dxkBq7KymRZW6gO20iPcSPWA++PgJ6qOgC4ELg7mYG3npLDN34Iu/udlsxijYlbzNGe6KK0qOo8IDuZgYPifUhsc/qNK7GS/8C1FaqSGbjNV6XM/R2E3n43mcUaE7f6vI6Y1IVG0iPC4DWgxV8ms1hj4hbrabO3iJTWOM+MngugquprZmco3XuBpXqPjfYYN2Ilf9eGDByMJn+koqIhwxhTpzqTX1VXN2jgvS2/fcJrHPG9BEmiQpnNWdAedrdO+rQhY+LiLPmDrVrT72ZYN+JMV1UwKc5Z8oeD3vSGPdU2t8e4EVfyi8hIESkUkeUi8rKI+H6HNy2YxvzfQMc/vOa3KGMSEs/SJdnA74FLVLU78DnwS7+Bw8EwvTZCuNj24jVuxNPynwvMV9XPoufTgKtE/K0rHg6G2RMErbBuj3FD9BA7IorI3UBnVR0XPQ8BlUCr6KK1B95vW06YJkVVa22o42n5A9Q+tWG/GWkicpOILEigbsY4Ec9k+jXAgBrnxwFbVXVXzZtU9Sngqb0t/6H+j7Jt9zZeO70N3QZcwBm/+Uf9at0IVJUt5VsoKSvZd2wu3/z192WbKSkvobSi9KCD8t20LYPMyv2Pf3eC3WnQ50s4cw1k7YHmVULL6hAtq0M8OqodkRZZjP6knAsWbCc9IqRXQbgK0qqVhx8eRaBZBqNe+JiB/1xMqLKaUGU1wcoqNBDg6Xn/TzgYZvADz9L1tf3XG6hs1YL8JW8RCoTocuNdtP77W/tdrzq+A5uLPiEUCNHy25eT9ubb+/8HOflkWLLE+37QIPjgg6+vicDAgfDvf3vnAwfC0qUQCHjLUgYCMGQIvPSSd/3ss2H9+v2vDx8OU6Z41889F0pL979+/vlwzz3e9VGjvMXOgkEYPpy8Z59lwYI62906u+fxdHvaAYuAM1X1s+jaPcfW3Jz6gPvjSv6yyjKyHsxi8jcnc+egO2Pe2xB27tnJyi0rWbFlBV9s+4L1O9Z7R+l61peuo3TzBlqVVpJdBtm7ILsMXj8RilvA2evTuHV+iLZ7QrTaI7TcrTQvr+bJ+0eys3sOQ/9RxPlTZh4U89N5swie1IPsac/R8r5J3g9FIDPTW6x38WI45hh48kl4+mnvbbf0dO8Ih+Hll737pk+Hd975+tre4777vPLmzIFPP4VQCNLSvK8ZGTB6tBczPx82bfJ+vveejAzo08e7vm4d7N799e+GQl78Nm286xUVXpy9iel4W6m8vLyGSX6A6Jr8DwFhYCUwVlVrXXYh3uSvilSRNimNSUMncd/g+w5Zh0Tt3LOTwq8KWVi8kILiAopKivhi46ccteorOm+DjqXQcTvk7Azyl2Ht2djnRC5YIdz5i3cOKuurl35P81GXkfnWXGTCBGjVav9j4kTo2hWKiuC99yAryzsyM72jb1/v686dXgJlZkKzZs6T53DXoMlfH/Emv6rym34BhoVO5KT/fBbz3vrYsGMDc7+Yy9zVc5m3ci6hpcvI3QC9NsLcns356sw+DNt2FA/c9rev6xIOQ8eOyKOPev9LXbvWa12zs/c/OnTwWljTpCSa/M5eoBUR2pRD6x3+1+0p2lTEjKIZ/LXoFT4pzqdVObz9pwCPFytp0VdwIi2a84Pv/By57lZvQ4zjZkKXLtCpE5KdvX/r27Ej3HGH73qZps3p2+PlIUgrT2ycf3fVbqYvns5Lcx6n89ufMGYJ5HY4msKHf8k5OcM4bdUjBE44AfLyIC+PQE7O1wmekQFjxiTxX2IOR06TvywkpFVU1ut39lTv4elPnmbO73/KNW9vZuZnEFSo7N6NtItv5IIzoy32iy82QI3NkcRt8gcDhHfHn/zvr36PG2fdxLKSZbywsgPnbm5N4K5xcOWVpPXsaQ+Opl6cJv+StgE+6ns0Z0Yi3pBZHaoj1fzfc98j9/6n6DuiHY/cOosR4wchGRn2AGoS5jT5/9ItndW5PZgTI/F3Vezk+Rv6ccufl1GVmc4zw6bS7KSRjVhLc6RymvyBqgA79uyo8/qO7Zv4YERPxn2wkdVn9eaEl2bDscc2Yg3NkczZyywAucXKW+M/gTcP3uuiOlLN03cM47wPNrLo+5dxwtwCS3yTVE6TvywQpFV5BDYePKf/7rfu5rYOi5nxh3vo9cRL9jBrks5pt2dLKBq+uPjrH6ry+fixzK78I7d86xYuGfmgm8qZI57Tln9HMMjmDNDly/f9bPcTU+jy+B8Zt6otj577qMPamSOd0+QPVaaxNBuqlizyfvCvf5F26+3M6gan/fpVMtIyXFbPHOGcJn9aeRqv9oBNI4fAmjVUXzKaFW2UV39yKQM725ImpmE5Tf5wRZjHzoDFlw+FyZOpKCvlsitDTPzWwy6rZVKE85YfvGnIRXddz+nXVnHOBd8np02Oy2qZFOF0tCetIo1wMMySjUuYUTSD1R1bMHHwRJdVMinE7Se8kQD9j+vPI/MeAWDyNydzdObRLqtkUojTbg/AVb2uAiCvQx7jB4x3XBuTSpxvhTgubxyDOg7ixKNOJD1kMzRN42nId3gPOR9BRBaoal5SK2BSTqJ5lPSWP56kr+GpZMc3KSmhPEp6y2/M4cL5A68xrljym5RlyW9SlpPkb4idXoypr0ZP/oba6cWY+nLR8jfITi8mtYhIQESmisiHIrJURIpEZFB9ynCR/B2BtTXO1wEtgRYO6mIOXwOADsBAVT0ZeA64uz4FuJjeENdOL8bEoqrzROQ+4GYR6QoMAepeB6cWLlr+NXh/sXvVutOLMbGIyEjg79HT14BfE2M58tq4SP43gNNFpFv0fBxe5Y2pj+HATFWdBiwAvg0E61OAk+kN9dnpxZjaiEgP4M94XfcQXqN6CdBJVSNxlWFze0yqsk94Tcqy5Dcpy5LfpCxLfpOyLPlNyrLkbwAikiciBQccm0RkRR33B0Vklogck6T4j4rIkGSUdSSzoc5GICLfAN4DrlfVgz7QE5E7gYiqPpKkeK2A94H+qlqejDKPRJb8DUxEjgLmA0+r6kO1XM8EVgC9VHWziNwPdAbaAycA64GrVfVLEfkCeAEYBrQBHgYGAblAJXCRqm6IlvtroEhVpzboP/AwZt2eBiQiIeBlYF5tiR81DPhUVTfX+NlZwGWq2gPYhTcFZK9mqno68FO8VQumquqpeDNlr61x3xvA6KT8Q45QlvwNayqQBdwQ454eeC1/Te+qamn0+3zgqBrXZkS/rgSKVXVhjfOa960CuidS6VThfMW2I5WI3AxcBPRT1d0xblUOboTKD7hec7ZiRY3vY+3gXYlNE4/JWv4GICKD8SbujVLV4kPcvhzo2gDVyAGWNUC5Rwxr+RvGz6Jff1fL25m5qlqzRX4LeEZEWqvqtiTW4Xy85w1TBxvtaQJE5F6gSlWTsiVNdDWMD4C8Q3S5Upp1e5qGR4BhIpKsXbbvByZY4sdmLb9JWdbym5RlyW9SliW/SVmW/CZlWfKblGXJb1LWfwFxahjSTXw00gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 198.425x198.425 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Z = np.linspace(-0.4 * bls.Delta_, bls.a, 1000)\n",
    "fig, ax = plt.subplots(figsize=cm2inch(7, 7))\n",
    "for key in ['right', 'top']:\n",
    "    ax.spines[key].set_visible(False)\n",
    "for key in ['bottom', 'left']:\n",
    "    ax.spines[key].set_linewidth(2)\n",
    "ax.spines['bottom'].set_position('zero')\n",
    "ax.set_xlabel('Z (nm)', fontsize=fs)\n",
    "ax.set_ylabel('Pressure (kPa)', fontsize=fs, labelpad=-10)\n",
    "ax.set_xticks([0, bls.a * 1e9])\n",
    "ax.set_xticklabels(['0', 'a'])\n",
    "ax.tick_params(axis='x', which='major', length=25, pad=5)\n",
    "ax.set_yticks([0])\n",
    "ax.set_ylim([-10, 50])\n",
    "for item in ax.get_xticklabels() + ax.get_yticklabels():\n",
    "    item.set_fontsize(fs)\n",
    "ax.plot(Z * 1e9, bls.v_PMavg(Z, bls.v_curvrad(Z), bls.surface(Z)) * 1e-3, c='g', label='$P_m$')\n",
    "ax.plot(Z * 1e9, bls.PMavgpred(Z) * 1e-3, '--', c='r', label='$\\~P_m$')\n",
    "ax.axhline(y=0, color='k')\n",
    "ax.legend(fontsize=fs, frameon=False)\n",
    "fig.tight_layout()\n",
    "figs['a'] = fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Panel B: recasting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAACOCAYAAABgxUxiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2deXhU1d3HP2cmOyRAIBBQhApCQdpKAogL+iiKgoDUpVoVHxe2KrhvZXlZjFBqBSGiBGjr2lLkqdpYtL4W9K2tYBPckEVQIIQECAGyTzKZOe8fkztOkplk7uTeuZPM+TzPfZKcu/x+Z27mfM/vrEJKiUKhUCgUkYrNagcUCoVCoWgJJVQKhUKhiGiUUCkUCoUiolFCpVAoFIqIRgmVQqFQKCIaJVQKhUKhiGhMESohxHNCiAIhxBcNx18a0hOFEH8QQuwSQnzT8HuiGT4oFAqFomNgVkR1MXCrlPKChuOWhvR5QAzw04YjEfi1ST4oFAqFogMQY/QDhRDxwHDgCSHEAOBb4GEpZQHwf8AhKaW74drPgfON9kGhUCgUHQczIqo+wFZgPp6oaTvwjhBCSCk/kFJ+CyCE6Ac8BLzp7yFCCOl7mOCnQqFQKNoDUkpTD0AA5cCPfNIygQJgbgv3Sd/DzzUKRXvA9O9Ya0fT79LEiRPlxIkTZdP0nJycZmlr165tljZ9+nSZkZHRKK13795y4cKFQT1z4cKFsnfv3o3SMjIy5PTp04Oyf/nllzdL05Onyy67rE158me/rXny52c487Rq1ao25cmffT15uuKKKwKV80gpEdKAtf6EEEuAyQ1/fg+8JaV8reGcJlRDpJSFQohbgReB2VLKPwXxbM83TUrR5FTbHVcozKfp/234HRBCtuV7ft999zX6+8UXX9T9jGeffZaDBw/y2GOPce655+q+v6KigieffLLNfmRlZVFUVMS8efM466yzdN9fWlrKggUL2uzHggULKC0tZcmSJfTo0UP3/UVFRWRlZbXZjyeeeILKykqWL19OcnKy7vsPHDjAihUr2uzHAw88QHZ2tr9yHjCoj0pK+T/A/wAIIYYB/xJCfCKlPAj8CviqQaQmAauBcVLKPCNsKxSKlsnJybHaBYUiGGYGOmH4YAop5S4hxBwgVwhhBwqBXzac/h2eGuYGT6AFwL+llPcb7YeRVFZWsn//fo4fP05dXR1xcXGkp6fTv39/unbtaqrtkydPcujQIU6cOIHb7fbaHjhwIElJSabZra+vp6ioiOPHj3Ps2DHsdjvx8fH06tWLgQMHkpCQYJrturo6vv/+e4qLi6moqMBms5GWlkZaWhp9+vQx1bbD4eD777/n+PHjOBwOnE4nqamp9OzZkwEDBmC3202zbRYzZ85kxowZVruh6IAUFBSQl5eHy+XCbrdjt9tZvXo1CQkJJCQkkJiYSEpKivfo0qULKSkpdO7cGZut2RCJHGCdPzuGCxWAlPJ14HU/6YPNsGcGx48fZ/Hixa1ed84553DVVVeRmZmJj/i2CafTyeOPP05dXV3Aa4QQDBo0iIkTJzJgwABD7AJUV1fz2GOPtXiNzWZj2LBhXHfddfTt29cw2/v372flypUtXhMbG0tGRgYTJkwgLS3NELtut5vXX3+d7du3t3hdcnIyI0aMYNy4cXTp0sUQ2wpFe2Xbtm1MnjyZqqqqRumtfY80kpOTvQJ26tSpFq81RajaO1988QXr1jUW9tTUVIYPH05iYiIOh4PCwkIOHTpEQUEBf/jDH9i+fTt33303nTp1CtnuZ599xssvv9wsvW/fvgwePJi4uDhqa2s5fPgwBw8eZN++fezbt49Ro0Zx2223ERcXF5LdY8eOsWTJEr/nzjvvPPr3709MTAw1NTUcOnSIw4cP89VXX/HVV18xfvx4rrvuOn+1o6A4cuQIX375JVu2bGl2bsSIEaSnp+N0Ojl58iQnTpzgyJEj7Nixg7y8PH7xi18wZsyYkOz6Mnv27GZpI0eOJDU1FbvdzunTp/nuu+84ceIE27Zt49NPP+XOO+/kggsuaLNthaI98s033zBlyhSqqqro168fvXr1wuVyUV9fz5gxY6itrcXhcFBdXU15eXmzo6KiwnscPXq0VXsdcgml/v37I4QgMzMTgBkzZiCE8B5FRUXk5uY2StOE6c4772wmUiUlJTzzzDPcdNNNXHfdddxwww0kJCSwevVqPv74YxwOB7t37+Z3v/sd8fHx3mdqzS1atCWEoE+fPgAsWrSokf3c3Fy/IjV37lyys7O58cYbmTRpEkuXLuWRRx6htraWvLw8nE4nn332GatXr+btt9/2myfftEmTJgEwadIkb1ogkRo7diy33HILN9xwA5MnT+aWW27hwIEDLFu2jOLiYtxuN++99x6TJk1CStksT/n5+eTn5zdKW7RoEQB9+vRBCMGyZcv8itTvf/97rr32WlwuF1OmTGHatGnMnTuXIUOGMGrUKFwuF3/+85/JzMz0myctul23bl2zz7moqAghBGPGjGk2UKCuro6cnBzmzZvH9ddfT15eHnfeeSeLFy9m8+bNnH322TgcDtatW8egQYPo379/i/+LCkVHo6amhptvvpny8nKuueYaxo8fT2ZmJqNGjeLiiy8mOzubdevW8eqrr7J582Y++OADtm/fzu7duyksLKS8vByXy0VZWRkFBQXs2rWLm266qUWbHS6i8m1+27lzp9/mOH+jfWbOnMnMmQH78gLa2L17N7t37/Z73fr161m/fn2jtOLiYr8+TZ48uVkawEsvvdTo76Z5+u9//xvQT395evfdd5vZb2qjpXR/eQL8RlQjRoxolrZ48eJGTaqBbIP/9/TAAw80Swv0nlv7nD/55BM++eQTv7b9vaeSkhLmzZsX0N9IZe3atVa7oOhAzJ8/nz179jB48GCeffbZFr/DgbDZbN5mP4DevXsDzAp4fajOKhQKD06n02oXWmTWrIDff4VCF//+979ZuXIldrud1157jcREQ5dqDVij6nAR1ZYtW8jNzSUxMZHf/va3ukZpNW0Gmjt3LmeffXbQ92/dupXNmzczePBgHnzwwaDv82d71apVxMbGBn3/pk2b+Oijj7j44ou544472mRbzzwIKSVr1qxh9+7dTJkyhXHjxoXNttvtZunSpRQVFXHPPff4jeACIaXk/vsbDzbVY9vhcLB8+XI+/vjjkPvnFIr2hMvl4v7770dKyVNPPcXIkSM5cOBAWGx3uG/Y+PHjOX36NDU1Nezfvz/k5yQkJOgSKYDRo0djs9nYv38/FRUVIdvOysrSJVIAl19+OQD5+fnU19eHbPu5557Tdb3W1wOQl9e2qXFNJ1K2hs1m89resWOHrntdLlejvxcuXKjr/oSEBObOncu2bdva5ZB1hUIv69at48svv6Rfv35hbwLvcEIFcPDgQQDuuuuuoDv1hRCNmnBWrlwZ1OAL37RbbrmFH//4x7jdbkaNGhV0p74QwltwJiYmcvXVV7c6+KJpntLT07Hb7dTW1pKRkeFND2ZASXl5OQADBw4kMTGx1cEXTfOUkZFBbW0thYWF7Nmzp9F1rQ0o0di+fTtFRUVBvSffPE2YMAHwjNQsKCgI6j1NmjSpkVDt2rWL9PT0oN6Tb55Gjx5Nv379QvkXVSjaFaWlpcyfPx/wVGYNbvJrlQ4pVFrn3u23346UkszMTDIzMxutHaUVgEVFRd40LSoZOXJko2v79OnjHdmmHVph5ZuWm5vL0KFDAU+zoc9aa8yYMaPRtZMmTaJPnz7ev88/37OI/D333EN+fr43vaioCPAU6r73+8uTFl389re/9abl5+cDngI4UJ4GDhzo/bz85QkgNze36fpx3jy5XC6GDx8OeJrEfK/ThCJQnrTIcfv27UG/J988VVVVkZ6eTmxsLG63O+j35CtUH3/8cdDvqWmeDh06FMJ/aHiZPn261S4o2jkLFizg1KlTjB07lhtuuMEsM81HaTXQIYXqnHPOATyzpvWscaZdO3hw6POStRp2YWGhrvu0gjMmJvRuQy3fodpuSxOWlu+CgoKw29YmHevJt2a3U6dOhk3UjlQ0YVcoQuGLL74gJycHu93OqlWrzPy+ZAY60SGFKjU1lYSEBKqqqprNmm4JrfBqS+d4z549Ac/SR3pE0u12A20rsHv16uW1rQcj8q01nZWUlAR9j5TSm28jbOvJtxEVg/bCzp07rXZB0U6RUvLAAw/gdruZM2eOt+XHJDICneiQQiWEoFu3bgCcPn066PuMEIvOnTsTHx9PTU0NNTU1Qd9nRGTRvXt3QL9QGZHv1NRUwNOWHSy+AtmWWpr2rltbhsWfbTUQQqEIzMaNG/nXv/5FWlqa7gFHRtIhhQoISaiMKLyEEN5JbHpG/hkR1SQnJ2Oz2aiurtY18s+IfGuf95kzZ8JqF34QyXC/a4WiI1NZWeld93PZsmWmL8DdEh1WqLS9VcLd9AeeqAo8LzpYjIhqhBDetQbDnW/NbnV1ddD3GJHnUG0b9a7bAw2z/hUKXTzzzDMUFRUxcuRI7r777nCYLA50osN+S7Xhk3qa34wqODWhCndEBYQkVEbkOy4ujpiYGJxOZ4urvvtiVFSjvWs9QqVFnNHQR6W2+FDoZf/+/d45ldnZ2eGq0Pnd4gN0CpUQ4i0hxNi2+2M+2l5NVtSy2xLVGBVd6InmjMi3EEL3Z25UntvyrqOh6S+Y7WoUCl8eeughnE4nd999NxdeeGG4zAbsBNNbMn0CvCCE2COEmCOE0L93cZjQatkOhyPoe4yKqOLj4wGCjiyMtG1lvkMVqrZWDOLj47HZbDidzqDX3YsmoVIo9PDmm2+yZcsWUlJSWLZsmdXuADqFSkr5nJRyCJ5Vbi8EvhNCvCiEMHXMYiiE0hxkVMGp7QulR6iMth1sge12u5FSIoQIe76NEgshhG6BVkKlUDSnpKTEu/7m8uXLvVNerCakBnop5cfAx0KIVGAq8KoQokxKeaWh3rWBUKIaowovKyMqbaWHcItFpNhWEZUiGigvL2ffvn04nU5sNht2u52NGzcSGxtLUlISKSkpJCcne3fRTU5ObnVj1bq6Om699VZOnjzJlVdeGfS2R+GgrT3JtUAVUA70aLs7xqF1kusZpm2VWIBxBadVUQ1YKxZKqAKTk5NjtQsKA9mzZw/XXHMNR44caZS+devWFu+Li4trJFzJyckUFRVhs9k4cOAA+/bto6CggF69evHyyy9bsWJLQGUMSaiEEJcA04Drgf8FFjVEWRGD3oILorfpzwi7odhWQhUeZs6cqUb+dRBOnDjB+PHjOXLkCKmpqaSnp+N2u3G73QwfPpy6ujqqq6upqKhotOV7eXk5dXV1lJaW+p2U/9133wGe3dH/+te/epclCzM5BBj5p0uohBBPAPcAnYANwPlSyoBj361Ei6j0CJWVgymMKjijtflNbwQdTUKl6BhIKZk6dSqHDx9m+PDhZGZmNtoOqKX91KSU1NbWNhOw559/noqKCqZPn86AAQP42c9+FpHfCb0R1bXAfOAtKaWrtYutRHuBoazQYGVE1dZ5PaE2/RkRUVkpkqHmOxrmUSk6Bm+++SYffPABqampbNiwgQ0bNgR9rxCChIQEEhISSEtL86bn5uZSWVnJtdde610kIRLR9S31HSwhhEgHUpuc322QX20mlKY/o/uogrVt1OKsodiO9ogqGlamULR/KisreeSRRwDPckbaQszRQqh9VCuA+/EMotCQQE8jnDKCtvRRtbXgbEuh2dYOTL1RjVHiDKqPKlKxcjFRhTE8/fTTHD16lJEjR3Lvvfdy/Phxq10yg4Az00OtTt4A9JFSpvkcESNSENqoP6Nq2Vb2l2gFdtOt1luzbWTTnxKqyELb6FHRPtmzZw8rVqxACMGaNWs68v9swBE/oZZO3wLBL5NtAVYOprBSqLRnBGvbyIhKy7dekbRCqKJprb/i4ogc76QIAiklc+bMob6+nunTpzNy5EirXTKTgKsnh/otXY1nwu82wFsySCmXhPg8w7FyMIWVQtWebHcUkVQozGLz5s3885//JDU1laVLl1rtjmWEKlRP4emfsm6DklaIhIhKbzOUEbV7KwcV6I3mjIxq1PB0RUejsrKShx9+GPAMoNA2Ro1GQi0hOkkpLzXUE4PRW8P2vdYoobKidq+3wO4oYqE9I9jP3MhoLtLJyAi4w7cigsnKyuLo0aOMGDGCe++912p3wsHOQCdCrUbvE0L8NMR7w4I2gs536HdrdISmP70FdkexHapAR4NQZWZmWu2CQid79+7lueeeQwjBiy++GBX/p0B+oBOhlsjnAHlCiH1CiK+0I8RnmYaewsvIuUxWNr+FWmD7znAPt20jojkrRTLSWb9+vdUuKHTgO4Bi2rRpHX0AhS/TA50ItYT4dYj3hRW73Y7T6Qyq8PIVqfYcUUWC7fYQUUWTUCnaF5s3b+bDDz+M+gEUvrRlm4+IR0/hpQ18sLKvxgrbZkQ17cG2EipFJOK7AsXSpUvp0SOiNqWwjA69foye5iAjxcK30JRSBm1bRTWho7fpz8gmT4XCKBYtWkRhYSEjRoxg2rRpVrsTMXRoodJTcBpZu9c2MgN9ItneR/2pwRSRydq1a612QREEn3/+Oc8//zw2m421a9dGxf9mE2YFOtFmoRJCLGrrM8xCT8Fp9EoFegpOK1emiASxMFIkgx3hGU0rU8yaFfD7r4gQnE4nM2bMwOVyMWfOnGgdqRmwRmXEt3QysMiA5xhOKBGVUbWYmJgYamtrwy5U0TqPysp8K6IHKSWnTp3C4XB4W0727NlDfHw8SUlJJCcnk5SUpHtx6Xnz5pGXl0ffvn15+umnTfK+/WLEtzTs+xUHi1V9VL7P0SNURg5Pt7L5TW8/kRXNjkqoFHqpqKjgl7/8JR9++GGj9M2bNzf622azNdvyvenvO3fupL6+nuXLl7N//37eeecd7HY7f/rTnyJ6XyirMOJb+jcDnmEKVvVR6bVtxqg/KxZnjdZmR0XHp76+nuuvv55t27YRGxtL9+7dcbvduFwuevTogcPh8G4B73A4KCsro6ysrNXn7tixA/D8H27YsIFLL43oBX8so83fUillxG52o6fgNEuoghEMM/qo3G43bre71Sito4w4tHJofKSTk5NjtQvtnoULF7Jt2zZ69uzJVVddRZcuXbznmm4B73Q6vVu9N936Xfv9rbfeoqysjLFjxzJw4ECuvfZa+vfvH+ZcRRwzA53o0N9SPQWnGX1Uwdo2ssAWQhATE0N9fT0ul6tVobKyj8oM26rprzkzZ85kxoyAW/0oWmHr1q0sW7YMm83G+vXref/991u8PjY2ltTUVFJTUwNec/LkSUpLS3nqqafUXKkfyAH8bp7WoYenR0JEFe7BFHptGynQVjb9hWo7GoRKETp1dXXcd999SClZsGABl1xyidUuRSUdWqhCiWqsECqjozk9hbYZ/WN6m/6ssG3kSiSKjkt2djb79u3jvPPO49e/bhcrx3VIWhQqIYRdCHGDEGKKECLGJ/1m811rO6FEVFZGNe1dJEPtJ1J9VIpI5NixYyxevBiAVatWER8fb7FH0UtrEdWrQAZwAfCJEGJgQ/qvTPXKIELpozJqSR09gym0a+Li4gy1He5IMlqHxkc6EydOtNqFdslTTz1FRUUFEydOZPz48Va7Ew28G+hEa9/SPlLK2wGEEK8Af4zklSiaEkoTmBURldHNUHpE0qyoRkrZ6qTHSBgaHw1CpdDP9u3beeWVV4iLi2PlypVWuxP1tBZRxQsh4gGklAeBScBjwDCzHTOCUCIqo8UiGNtGR1RWTXT23SIl3J+5iqgC8+67ASuqCj+43W7mzJkDwKOPPsrAgQNbuUNhEAFD/9aE6hGgm/aHlLICuL4hPeJpL6P+NKEyutnRymguGNt1dXWAMQKtVqZQGMUf//hH8vLyOOuss5g7d67V7ihopelPSrndT5oLeN00jwwkEiKqSBcqI8VCs11XV6crkjQi33qESkoZVaunRxLBbHtjJWfOnPGO7nv22Wfp3LmzxR4pIIQJv0KILsADwHCg0VuUUo4zyC9D0BNRWdlPpImFlRGV0QM5wi2Soa6taMT6ioqOw6JFiygpKWHMmDHceuutVrsTNbTWnx1KqfwmYAfeAmpCuD9s6ImozIgswNqIyop866kcGCnQeiIqo8U50lFLKAXHN998wwsvvIDNZmP16tW6V0BXtBlDl1AaDXSXUga36qmFhFJoGjVXIhKEyspoLtxNf3o+79raWiB6hEotodQ6UkrmzJmDy+XiV7/6FRdccIHVLkUjhi6h9AkwpE3uhIlojaisjC6CrRy4XC5cLpd3bUKj7Lpcrlb7QYx+14r2z8aNG9m2bRvdu3cnKyvLancUTQilhLgL2CKE2AEc9z0hpVxihFNGYVUzFFgXWfjatmowBbSeb1+BNKKJRQiB3W73CmBL4qeESuFLeXk5jz76KADLly9vcTFZhTWEIlTPAH2BQ0CKT3rEDecJJaLqSE1/kWzb6IqBZtvlclFfX9+iUGlNf2pJHAXA3LlzKS4uZvTo0dx9991Wu6PwQyhCdSswSEpZbLQzRmNlv0UoSyiFWyyklJYNpjBjQEOwTZ7RNphi7dq1VrsQsbz//vusWbOGmJgYXnrpJTUK1FpmBToRilB9D0T8QAqwtq/Gyua3YPOtLXUUExNj2Bc02CjWrIgKWv/Mo20wxaxZs5g5M+CAqnbBhg0bePXVV6mpqfFOK9i0aRPJycl069aN1NTURj+7du1KSkoKXbp0ISUlxXuUlJRQWlrK119/zYcffshDDz0EwNNPP60GUFjPWjwDKpoRilC9BvxNCJFN8z6qrSE8zzQioa8mGNsOhwMIf7OjlWJhZUSl+qjaF7/5zW8abbHhdrsBKC0tpbS0lEOHDul+5qZNm7y///znP+fxxx9vs58K8whFqO5v+Lm0SboEzm2bO8aiJ6Iyq+kvmOY3TagSEhLCattMsbAiqtEr0EqoIp/333+fuXPnIoTg8ssv57zzzkNKidvtZtmyZVRUVHD69GlOnTrF6dOnvb+XlZVRXl5OeXl5o98LCwtxu9307t2bAQMGcPXVVzNr1iy1QkmEo1uopJQ/MsMRMwhl1J8Vo9+klMTGxoZ9+SYzBhUEa9tocYbgKyZqMEX7oKSkhDvuuAMpJfPnz6e0tLTR+Z49e9KzZ09dz8zKyqKoqIh58+Zx1llnGemuwkR0d0wIIVYLIS5uknaxEOJ549wyhvYwj6qmxrO4hxViodm2QizMECq9/WPRElFNnz7dahdCYt68eZSWljJ27Fgee+wxq91RmM/6QCdC6UH/JZDXJC0fuC2EZ5lKsBGVlNJbaCcmJhpiO9hRf1rt3ii7vraDFaqOYltvs2O0RFT5+flWu6Cb/Px8NmzYQExMDGvWrFGj8aKDzEAnQnn70s999hCfZSrB1rBra2uRUhIXFxf2jRONHkgB7SOqsdJ2dXU1YKxIRjI7d+602gVdaPtBSSl56KGHGDx4sNUuKcJDRqAToYjLv4AsIYQNoOHnoob0iCLYGraVkYWZBXYkRzVW2tY+82gRqvbGG2+8waeffkp6ejoLFiyw2h1FBBBK7/2DePa2LxZCHAbOAYrx7P4bUQRbwzajr6Y9CVVHy3ewEVVSUpJhthXGUF5ezhNPPAF4ljNKSUlp5Q5FNBDKqL9CIUQGMArPUkpHgM+klG6jnWsremv3RhZckdAEZmVUE2zlwEjbWj+Giqga07t3b6tdCJqsrCyOHTvG6NGjueOOO6x2RxFeAq52FNJ46AZR2t5wRCztYVCBlcO0zSiwtcnDVuRbG8WnjeoLRLT1UbWXLT727dvH888/jxCC7OxsNYAi+vC7xQfo6KMSQhwRQqwTQkwRQnQyxi9zCbbg6mhNYMGKhZXD082oHGgDUrRRfeG0HcksXrzYahdaRUrJgw8+iNPp5N5772XEiBFWu6QIPwsDndBTZRkF7ACmAoeEEP8rhHhYCDGord6ZhVYAB1twGdn0p4lFa8PTzajdWzmowEqBVkLVfsnNzeUf//gHXbp0YenSpoveKKKdoJv+GlZL/z3weyFEDHAZMAF4WwgRB2xpOLZJKVsuKcKEtthqfX09Lpcr4NBzMyOL1grsqqoqADp1Mi5ItbLJUxNoK6LYYITK5XJRW1uLECJq5lFFOg6Hg4cffhiAJUuWkJaWZrFHikgjpEZgKWW9lHKrlPIxKeVQ4CpgHzCn4YgIfAujlgovK/uoKisrAejcuXPYbWsiaWQkqQmPFjEFQst3cnKyYbaDede+kZzqA4kMVqxYwffff8/555/PfffdZ7U7ighETx/VJUKI5QFOzwLypZTXSSl/Z4xrxqAVXi0VnFqBbaRY2O12hBDeBTTDadtKkQxGqOrr63E4HNhstrBHVGaIc6STk+N354SI4PDhwzzzzDMArF692rD1LhXtkoB70eipUs4F/i/AuY+AeTqeFTaCKbwqKioAYwtsCE4wzGj6CybPUkrLhMrXrpFRjZbvlpodtXdtZCQX6UTyXlQPPvgg1dXV/OIXv+DKK6+02h2FtQSsUekpJS4A3g9w7kNaWKfJSvQIldGFVzBCZYZYaE2YWpOmP2pra3E6ncTGxhraV6NXqIzEynet0E9ubi7vvPMOycnJrFy50mp3FBGMHqFKAQItNx0LROQ3P5jCy4z+EvhheHwg21JKUyIqX6GSUvq9xlcshBCG2Q5mpKVm28g8g7XvWqGPU6dOefujlixZQp8+fSz2SBHJ6GkQ3guMA97xc25cw/mIw8padlJSEmVlZdTU1NCtW7dm52tqanC73SQkJBjaNm+324mLi6Ouro7a2lq//UBmRzXBRFRmVQxasq0iKn0UFxezY8cO7zQKIQQPP/wwXbp08R5du3b1/tSOLl26BPyfrq2tZerUqRQWFjJ69Ghmz54dziwp2iF6SseVQI4Qwg68LaV0NyxIOwVYAzxihoNtRSukAzWDOZ1Ob8e+0fNqtA57LWpqill9Y+CJqurq6qipqfErVGYV2NpnaEXTn/Z5t9TkGY1CtXBhwHmULbJ3714uu+wySkpKGqXv2bMnqPs7derkFa7KyrW+8sQAAA0kSURBVEpqamr46KOPOHr0KOXl5XTr1o033nhDDaBQaAScma5nHtWfhBDpwCtAvBDiJNADcAALpZR/brObJqAVhoHEwrd2b2QTGPxQaGu10aacOXMGgK5duxpqFzwC3VI0Z5ZYxMbGYrPZcDqdAeeulZeXA8aLhZYXLW/+iMamv3Xr1rFo0SJd95SWljJu3DhKSkro3bs3gwZ55vVLKZk8eTJlZWXNjjNnzniPsrIyqqqqqKqq4ujRo97nnjhxAoChQ4fy2muvce655xqWT0W7ZwaenTiaoasqI6VcIYTYAFwEdAdKgU+llOVt9dAsWiu8ysrKAHMKLq0Pxgqhai26MMu2EIKkpCQqKyuprKykS5cuYbOtfd6VlZVIKf1WPMx835FKcXHAtT79IqVk+vTpHDlyhAsvvJCf/OQn3oncAI8++mhQz6isrPSK1ksvvURxcTF33XUXI0eOJD093fCKoaLdE3D15FDGBjuAfsClwO3AC0KIV4UQr4bonKm0JlSnT58G8Bt1tBVNLAIJlVZo+ivM20pr0ZyZ+dZEQGtma4pZQhUbG0tCQgJutztg06OZlYOOwiuvvMJbb71FSkoKGzdubCRSwSKEIDk5mb59+zJs2DD69+9Peno6559/Pr1791YipdBFKEL1KvAQUAF81+SIOFoTKq3gskKozCw0NbHQmtmssG1F5cA3qmqKlNJU2x2BM2fO8PjjjwOQnZ1N//79rXVIoSC0bT6uAX4kpTxjtDNmYGVEpdkOFFmYGVFpz2xNqMyMqALZ1vJthkh27tyZ0tJSKisrm60ZV1lZSX19PUlJSVG1zl9GRsAdvpuxZMkSTp48yZgxY5g6daqJXikUzdgZ6EQoEVUB0G6+5VYKlSYWmig05dSpU6bZ1nZGDSQWWr7NjKj8CXRNTQ0Oh4O4uDhTVi9v6X2b+XlHMpmZwc3F37t3L9nZ2QghWLVqlWqeU4Sb/EAnQm36e0cI8UshxJW+R+j+mUewQmVGga09U4sgmnLy5EkAevToYbjtloTK4XBQVVVFTEyMKUPjWxIqLc/du3c3pSBsqekvWpv91q9fH9R1jzzyCPX19UybNo3hw4eb7JVC0YzpgU6E0vSnzc5rummMBCJurGmnTp2w2WxUVlZ6lwzyRRsu27NnT8Nta0LlL6Kqrq6mqqqK+Ph4r6gYiRbN+RNJLc9paWmmrCDeklD52jYD7TPXRMmfbTMqBu2dLVu28N5775GSkkJWVpbV7igUjdAtVFLKH5nhiFnY7Xa6detGaWkpp06dolevXt5zVVVVVFZWEh8fb0o/kTY3q6Kigvr6+kYTG7VJlGlpaaZEFpr4tSRUZogz/CAWWlNbOG13794d8MwDasrx48cBGv0PKDyL+Gr7QS1cuNC0d6NQhEpUbMgTqPDyLTTNEAu73R5QMDShMqt2rzVvnTlzptm28GaLhRYtac18vmj5NluoWhJJJVSNeeGFF/j2228ZNGiQWs5IEZFElVA1LTi1GraZNUjNdtNlaMxuhoqLi6Nr1664XK5mhbbZQpWamooQglOnTgUUSbOa/lRE1Zy1a9cGPHfixAkWL/asXLNy5UrveokKhQXMCnQiKoRKE4OmhVc4hCo9PR2AY8eONUrXlpUxc9VorUDWxEFD88WsfMfGxtKtWzeklI0+cymld5UEs8QiNTUVoJlI+jbzRttk31mzAn7/mT9/PuXl5YwfP54JEyaE0SuFohkBa1RRIVRaodh0KZkjR44AcPbZZ5tmWxMqTRQ1CgsLTbetCZGvUNXX11NUVGS6ba1y4Jvv0tJSampqSE5ONqVPEDwi2aNHD9xudyPb2rtWqyL8wOeff86GDRuIiYlhxYoVVrujUAQkKoRKK5C1wgo8tfvDhw8DcM4555hmWxMqX5GsqamhpKSEmJgY73kzbfsuClpcXEx9fT09e/Y0ZR6Thr/PXPu9b9++porFWWedBfxQGQC877pfv36m2W1PSCmZPXs2UkrmzJnDj3/8Y6tdUigCEhVC1aNHD+Lj4ykrK/MOmS4tLaWqqorOnTt7m4vMQCuwCwoKcLvdABw8eBDwFKhmbnGgFcqHDh3ypmm/mynOvs8vKCjwpmn5Ntu2P5HU8q2EysPrr7/Of/7zH3r16hXyNiAKRbiICqGy2Wz07dsXgO++8yxJuHevZ5/HAQMGmFq779q1K927d8fhcHgjm2+//RbAu3WCWfTt2xe73U5xcbF3FXUt3+edd56ptn1FUttleN++fYD5+dZsa+/a7Xazf/9+gKjcViInJ6fR3+Xl5d71/JYvX25aM6xCoZOZgU5EhVAB3qYNbdO33bt3AzBkyBDTbQ8YMAD4oaDWbJtdYMfGxtKvXz+klOzduxeXy+X1weymnp49e9K1a1fKy8spLCz0/oyJiTFdLAYOHIjdbufw4cNUV1dz6NAhqqurSUtLi8o5QjNnNv7+P/nkkxw/fpzRo0er9fwUkUROoBNRI1RDhw4F4Msvv6SsrIxdu3YhhGDYsGGm2/7pT38KQF5eHkVFRRQWFpKYmGi6UAH87Gc/AyA/P5+vv/6a6upqevfubfrqDL6f7c6dO/nss8+QUjJkyBDTh0AnJCQwYMAApJTs3LmTHTt2AITlXUc6W7ZsYe3atcTGxpKTk2PKyiQKhdFEzR7Q/fr1o0+fPhQVFfGb3/yG+vp6hg4damr/lMawYcNITEykoKCA7OxsAEaMGBHSPj96ycjI4J133mHnzp0cOHAAgIsuuigsI98uvPBCPvnkE7Zu3erd6feiiy4y3S7AxRdfzLfffsvbb79NbW0tAJdccklYbEcif/nLXzh8+LB3eaQlS5Z4K1AKRaQTNdUpIYR3nkhZWRlCCK677rqw2I6Li+Oaa67x2o6Pj2fcuHFhsd29e3cuvfRSwNM30a1bN+/fZjNgwACGDBmC0+nE4XDQv3//sBWOGRkZpKenU11djcvlIjMz09Q5a5HOrbfeypNPPklFRQW33XYbTzzxhNUuKRRBEzURFXgKr6lTp7Jr1y5Gjx7Nj34UvmULr7rqKmw2GwUFBVxxxRXeFRTCwY033khCQgJnzpxhwoQJJCQkhM32vffey9///necTieTJk0KW1NTTEwMs2fP5r333iMpKSmqJ7P26NGDK664gtTUVK688kpuuukm1eSniETeDXQiZKESQtwBPI5n1fRq4AEpZZ4Qwg68AFzecOkW4HEppRRCpALZwFAgEXhGSvlaqD6EwkUXXRS25idfbDYbV111VdjtgmdQxZQpUyyxnZSUxM0332yJ7dTUVG6//XZLbEcSo0ePZtOmTVa7oVCETEhCJYQYDDwLZEgpi4UQE4C/AucAU4HBwE/wNC3+B7gJeBN4GdgjpbxdCHE28LUQYpuUstCPGYVCYQDvvhuwoqpQRBITA50INaKqBaZJKbXlFvKAdCFEHGAHOuHZBdgGxAGOhmjqauBWAClloRDiQqD5MtcKhUKhUDQQUkO1lPKQlPLvAMIzfGwF8DcpZR2eqOk0cBQoBg5IKXOBgQ1/PyKE+LcQIg9PRFbtz4YQQgohZCj+KRQKRUtok9AV7QQpZcgHnsjpTWA70LUhbQnwGp5IqguwDXgUuARPf9YDDddpwpUZ4NnS92iLn+pQRzQfTb9LQG7D0TR9hp+0mX7S1gH5TdKKgEVBPnNRw/W+afkNzw3Gvr9nqjx1gDwF+h8WDf/IrSKEWAJMbvjzb8CGBoN7gLullDUN1+0C5kgptzX8fReePqo5wPdAipSyouHcm8BWKeVLQTmhUCgUiqgjaKFqdJMQycCXwCtSysVNzr0KVEkpfyWEiAX+AuyUUmYJIfKBP0opXxBC9AJ2AlOklP/VaV+/0wqFRUgp1b4iCkUbCFWofg1kAV83OTW24ecLwHDABfwTeExKWSeEOAdYA5yLp3/seSllwPWdWrCvhErRblBCpVC0jZCESqFQKBSKcKGmpysUCoUiolFCpVAoFIqIRglVFCA8vCKEeMwnzS6EeF4IsVcIcUAIMctKHxUKhSIQSqg6OEKIIXgGtNzU5NRMYBAwDBgJPCSEGBVm9xQKhaJVlFB1fO7HM+ftzSbpP8czVaBeSnka2Ajc0fRmIcSlQojPhBD5Qog8IcSN5rusUCgUP6CEqoMjpZwtpfyTn1N9gSM+fxcCZ/u5bjGwQkqZCdwDXGm8lwqFQhEYJVTRiw3P0iUaAs+8t6ZsAtYIId4AMoG5YfBNoVAovCihil4KAN8tb/vgiaoa0TAh+yfA/wLXAF8JIcK386JCoYh6lFBFL+8A9wghYoQQXfFsv/J204uEEP8BhkspX8azcGRXID2cjioUiugmqraiVzTiJWAAnjUb44AcKeXHfq57AlglhMjC01S4WEp5KGxeKhSKqEctoaRQKBSKiEY1/SkUCoUiolFCpVAoFIqIRgmVQqFQKCIaJVQKhUKhiGiUUCkUCoUiolFCpVAoFIqIRgmVQqFQKCIaJVQKhUKhiGj+H7XbiDVmzgwcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 481.89x141.732 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run effective simulation\n",
    "data, _ = nbls.simulate(drive, PulsedProtocol(5 / drive.f, 0.), method='full')\n",
    "t, Qm, Vm = [data[key].values for key in ['t', 'Qm', 'Vm']]\n",
    "t *= 1e6  # us\n",
    "Qm *= 1e5  # nC/cm2\n",
    "Qrange = (Qm.min(), Qm.max())\n",
    "dQ = Qrange[1] - Qrange[0]\n",
    "\n",
    "# Create figure and axes\n",
    "fig, axes = plt.subplots(1, 2, figsize=cm2inch(17, 5))\n",
    "for ax in axes:\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "\n",
    "# Plot Q-trace and V-trace\n",
    "ax = axes[0]\n",
    "for key in ['top', 'right']:\n",
    "    ax.spines[key].set_visible(False)\n",
    "for key in ['bottom', 'left']:\n",
    "    ax.spines[key].set_position(('axes', -0.03))\n",
    "    ax.spines[key].set_linewidth(2)\n",
    "ax.plot(t, Vm, label='Vm', c='dimgrey', linewidth=lw)\n",
    "ax.plot(t, Qm, label='Qm', c='k', linewidth=lw)\n",
    "ax.add_patch(Rectangle(\n",
    "    (t[0], Qrange[0] - 5), t[-1], dQ + 10,\n",
    "    fill=False, edgecolor='k', linestyle='--', linewidth=1\n",
    "))\n",
    "ax.yaxis.set_tick_params(width=2)\n",
    "ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))\n",
    "# ax.set_xlim((t.min(), t.max()))\n",
    "ax.set_xticks([])\n",
    "ax.set_xlabel('{}s'.format(si_format((t.max()), space=' ')), fontsize=fs)\n",
    "ax.set_ylabel('$\\\\rm nC/cm^2$ - mV', fontsize=fs, labelpad=-15)\n",
    "ax.set_yticks(ax.get_ylim())\n",
    "for item in ax.get_yticklabels():\n",
    "    item.set_fontsize(fs)\n",
    "\n",
    "# Plot inset on Q-trace\n",
    "ax = axes[1]\n",
    "for key in ['top', 'right', 'bottom', 'left']:\n",
    "    ax.spines[key].set_linewidth(1)\n",
    "    ax.spines[key].set_linestyle('--')\n",
    "ax.plot(t, Vm, label='Vm', c='dimgrey', linewidth=lw)\n",
    "ax.plot(t, Qm, label='Qm', c='k', linewidth=lw)\n",
    "ax.set_xlim((t.min(), t.max()))\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "delta = 0.05\n",
    "ax.set_ylim(Qrange[0] - delta * dQ, Qrange[1] + delta * dQ)\n",
    "\n",
    "figs['b'] = fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Panel C: mechanical simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMgAAAC/CAYAAAC7f5PnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2deZgU1dW439vr7DPM4jCswy6b7AqiCALiAgYExJioIJGQRjQ/l8QPQ5LPJMZEyacmjEZEcRdUVBTcF1AWEdABlGVE9mX2faanp7vv74+qamaGWXqrnp6x3+fpp6urbt06XVXnnnPP3YSUkggRIjSOobUFiBAhnIkoSIQIzRBRkAgRmiGiIBEiNENEQSJEaIaIgkSI0AxhoyBCiEwhhBRCbGzk2Cr1WKqfef9ZCPGfwKX05PeAEOJmP8+dK4R414f01wohHvfnWk3k1+J9FEJcI4R4QI/rtzVMrS1AA+xAPyFEdynlUQAhRCwwtnXFqo+U8o8hvNY6YF2orqcyCkhuxeuHDWFjQVRcwGrgF3X2XQe8XTeREGKaEOIrIcQ3QojNQogx6n6TEOJfQoiDQojvhRBPCyEs6mnnCyE+E0LsF0JsFEJkqOdMFUJsEULsEEIcE0L8Rd0/Xs37BfU6e4UQY9Vjq4QQ96jbF6my7BVC7BJCXK7uv7WOjEeFEL9p7o8LIToKIT5U89hVRw6PxRFCfC6EWKbKlSOEuFf9vUMIsU8IMbhOull18q73W90XK4R4XgixVb1fO4UQ/YQQFwELgTlCiL81uH4XIcQ7Qog96v+9V92fKYQ4JIT4txBiuyrbDK+eeLgjpQyLD5AJVAAjgH119n8MDAIkkAr0AfYAKerxgcBpIBa4A9gIRKMo/2rgJuDPwI9AmnrOW8BSQACfAX3U/Z0Ap3qd8er2UPXY3cBGdXsVcA9gVq99jbp/hCpbArC1joyjgXJ1ey7wbiP/fynwpLodC7wKJNZND3wOvKFuX6Tek2nq7/8DnqqTbladvD2/69zHWcDjddI8Cfxb3f4z8J+G8qr39i51OxHIBm5Qn50EpqrHZgJHW/udCsYn3FwspJQ7hRAuIcQIIA+Il1LuFUJoSSYDGcAndfa5gd7AJOAFKWW1un8OKHUQ4CMpZb66Pxs4T0ophRDTgKlCiBuB/ihKE6umOyql/Fbd3oXystRlMOCSUq7XZFf3IYSYClwjhOgDDAXiWvjr7wMbhBDdUAqF+6SUpXX+o8Za9ftQnfO03+NbuIYHKeXrQogfhRCLUe7deBSlbpQ6ru4V6vmlQohVwFXANqAW2KAm34XqorV1ws3F0ngB+CVK6f9Cg2NG4BMp5VDtg1JC70Up8T2dy4QQ6ZorhfIANaRyWMQC3wDDUR7qvWo67a2sbnhOA1nqXU+95iAhRBfgW6A78CXwh5b+sJTya6AH8BRKibxdLSQaUtPgvNpG0jSU1dIwgeryrQSqgJeBVzj3/9XF0MhxA4oVBXBIKd1NXL/NEq4K8iIwG8UCvNzg2CfAFUKI8wGEEFcDu1Hcqo+BG4UQViGEAXgC+Hkz1+mD4g79QUr5DkopakVRQm84AEghxGRVluHAp8DFQD7wV+BDYKp6vMl8hRAPAUullG8BdwLfobiW/pAPjFTzHQBc0EiaKcAqKeVK9X9M4+z/dnL2xQdASlmOYikWqfkmAjcDH/kpY5sgLBVESnkS2AfkSCmLGhz7HlgAvCqEyAb+AlwrpawA/gvsVD97UOoHzYUodwPvAvuFEPtQXpLvUVwOb+SsQQki/EkI8S2KH3+dmucJlBdvH9AN5aVtLt9HgaFCiL3ADuAwSj3EH/6KUojsBR4ANjWS5hHg10KI3cAXKBZUk+9TYIoQ4t8NzvkFMFEIsQfYjuLurfJTxjaBkJHu7hEiNElYWpAIEcKFiIJEiNAMEQWJEKEZIgoSIoTCc1oLvLovSQixRd1+QghxWAjxt9aTMkJDwq6hsD0ihOgPLEdp/d5T59BUzjau/RroJqU8EWLxIjRDxIKEhkXA08BrDfb/DHhLCPEFSsPae0KIS0MtXISmiYR5Q4jaNWOvlPIRIYQV+ErtCYAQQqL0FStoTRkj1CdiQVqPiSi9AiKEMREFaT2m06Abf4TwI6IgrYBQuuiOATa3tiwRmieiIK3DRcAOKaWrtQWJ0DzBrKT/JGv7JSUlLFu2jCFDhjBr1qyWTwiARx55BKvVyoIFC7BarbpeK4wJbTf6II6++klRXFwslyxZImNjY6XRaJQLFy7U/ZojR46UFotFJicny8cff1za7XbdrxmGhHREYcSC+MGKFStYuHBh3eGydOnShQEDBuh63W3btlFWVgaAEAKTycTmzZsZNWqUrtcNM0JqQSIt6X5w1VVXcfXVV/Phhx/idrtxOp0MGTKExYsX63rdxYsXU1ZWhtVqxWAwcP311zNokL9jqiJ4RRDN0U+O06dPS5vNJk0mk7z99tt1v97YsWNlbGysfOCBB2RZWZnu1wtTIi5WuFFZWcnWrVsxGo2MHTsWi6X+EO/c3FwMBgNpaWm6ynH06FGSk5OJj4+vt7+8vJwtW7YQHR3NxRdfjMkUXo5BSUkJ27ZtIz4+ntGjR2M0ejuiuVEilfRw4p133pFpaWkSpQCQXbt2lV9++WVri+Vh9erVskOHDh75evfuLXfs2NHaYnl45plnZHx8vEe+gQMHyr179waSZUgtSERBmuGdd96RRqNRAnLQoEGyb9++EpAxMTFy27ZtrS2efOWVV6Tah0sOGzZM9uzZUwIyMTFR7t69u7XFk0899ZRHMUaNGiW7du0qAZmamioPHjzob7YRBQkHjh8/7imZf/e730m32y0dDoe8+eabJSC7desmS0tLW02+nJwcGRsbKwH517/+Vbrdbmm32+V1110nAdmvXz9ZVVXVavJlZ2dLi8UiAfn4449LKaWsrKyUU6ZMkYAcPny4dDgc/mQdUZBwYObMmRKQV199tXS73fK3v/2tfOSRR6TD4ZAjR46UQEgq5k0xefJkCcgbbrhBut1uz/6qqirZv39/Ccj777+/VWRzu93ywgsvlIBcsGBBvWNlZWWye/fuEpD/+Mc//Mk+oiCtzZYtWyQgo6Oj5YkTJ2RZWZk0m80yISFB1tbWyt27d0shhDSZTPLQoUMhl++DDz6QgExKSpIFBQXnHN+8ebNH/lOnToVcvjVr1khApqeny/Ly8nOOv//++x75i4qKfM0+oiCtzcSJEyUglyxZIqWU8i9/+Ys0mUzSYrHI5557TkopPa7WLbfcElLZ3G63HDVqVIsl8PTp0yUgFy9eHELppHS5XPL888+XgHzyyScbTeN2u+WECRMkIP/whz/4eomIgrQmu3fvloCMjY2VxcXFsqysTMbFxXkqm506dZK1tbXy0KFD0mAwSLPZHNJSWrMOKSkpsrKy0uv/ESo2bNjgifY1V8f44osvvPofjRBSBYn05m3AY489BsC8efNISkriueeeo7q6GovFQlRUFGfOnOHdd9+lZ8+eTJ8+ndraWp544omQyffoo48CsHDhQmJiYppMN3jwYCZNmkRlZSUrV64MlXge+W6//XbMZnOT6caOHcuoUaMoLCzkpZdeCpV4vhNEbWvzFBYWyqioKAl4wpBnzpyRGzZskDNmzJDz58+XGzZskCUlJVJKKTdu3CgBmZaWJmtqanSX7/jx49JoNEqTySRPnjzZYvp33nlHArJHjx7S5XLpLt++ffs8YfDCwsIW07/44osSkIMHD64XaGiBiIvVWmhx+8mTJ59z7K677pKPPPJIvX1ut1sOHDhQAvKtt97SXb6HH35YAnLWrFlepXe5XJ6I0WeffaavcFLKP/zhDxKQt956q1fpa2pqZGpqqgTkzp07vb1MxMVqLVavXg3Az3/e3ITwZxFCMHfuXABWrVqlk1Rn8VU+g8HAzTcrSynqLZ+U0mf5LBYLN954IxCa++cXQdS2Ns3p06c9le7GKrWNWRAppTx16pQ0GAzSZDLJvLw83eTLycmRgIyPj/epAVA7LzY2ttGQa7DYuXOnBOR5550na2trfT4vJSXFWzc1YkFag9dffx23281VV11FUlKS1+dlZGRw5ZVX4nQ6eeWVV3STTyudp0+fTnR0tNfn9e7dm0suuYTKykpef/11vcTj1VeVlRpmz57tU2fJYcOGMXjwYAoLC1m/fr1e4vlNREFUtAd8ww03+HzuTTfdBJx9ifVAk2/OnDk+n6u3fLKOe+Xr/RNChOT++U1LJgZlOTAXypJi2icbuLVB2jbLsWPHPC3PTbkhTblYUkpZXl7uiX4dP3486PLt3btXArJDhw5+Rcvy8vI80S9voku+ovU86NKli1/RssOHD3vcQC/aRMLSxaqW9dcEvBpYJoRobGmvNseaNWsAmDp1KnFxLa21eS5xcXFcc801ALq4MVrJOnPmzHPGonhDWloal19+OU6nk7feeivY4nms2/XXX4/B4LtTkpmZyYUXXkhlZSXvvfdesMULCL9cLKkskZYD9A2uOK1DIO6VxuzZs4GzyhYspJRBle+11xpODxwYLpfL85/D8f4FTEsmBnX98gb7xgBFQNc6+9sk3kaHmnOxpFTcrOjoaAnIY8eOBU0+f6NDDcnPz/e4WY11cPSXTz/9VAKyZ8+evjT2ncORI0c8jYwtuFlh6WJFCyG+VT97gb8Dv5BSHtcS9OvXT9coiV74Gx1qiF5ulr/RoYakpqYyceLEoLtZda1bI2u6e0337t256KKLqKqqYsOGDS2fECK8vePVUp2FvDGEECNuu+02Xn311XMmT+vduzdOpzMQGc/h8ssv55lnnglKXv5EX1588UUefvhhhBDExMTw+OOPM3LkSGbPns3rr7/OXXfdRadOnfyKONVFBhAdaozZs2fz4YcfsmbNGubPnx9wfrW1tbzxxhtA8OT76quvWLNmje6T8HlNSyaGRlysRtI8tGnTJjlhwgRPPyUNk8nk6QkbzM+ePXt8teLn4Et0SHOx9u/fLzt27Ojpwbt+/XrZtWtXKaWUFRUV0mAwSEAOHTo0YPkCjQ41pKCgQJpMJmk0GmV+fn7A+b333nsSkP379w/IvdLQ3KzOnTs3506G1MUK1vQXl44dO5apU6eydu1a5s2b5zmQk5MTpEso/PGPf+SFF15g9erVAc8J5U90yGq18vTTT5ORkQHAyJEjOXPmDA6Hg9zcXI+bcfDgQbZt28bo0aP9li/Q6FBDUlJSmDRpEu+//z5r165lwYIFAeUXLPdKo3v37nzxxRdceOGF4TMzS6AaBgwDjk+ZMkWOGzeu0Y5+weTjjz+WgOzTp09ApZbb7ZZ9+vSRgPz4449bTN9UZ8Vf/OIXcubMmVJKKe+99145evRoTw/f66+/3m/5nE6n7NixowTk9u3b/c6nIc8++6wE5OWXXx5QPtXV1TIhIUECcv/+/UGSzitCakGCoSB/A6Zq0o8ePVqePn066HdFo7a21jMNz65du/zOR4sOpaenS6fT2WL6hgpSUVEhZ82aJS+66CJZXFws7Xa7TE1Nla+99ppngJXJZPI7ohWs6FBDioqKpNlslgaDQZ45c8bvfN58803PbCohJiyjWM1xJfCx9mPGjBm6dhkwmUyeClwg19Hcg1mzZvk8kdmxY8e4+OKLMRqNfPbZZyQlJbFmzRqKi4u5++67tYIDKSX//ve/A5Jvzpw5QXFfNDp06MCUKVNwu90BRdsC6frSpgiitoWMzz//XAIyMzPTr9LV5XLJLl26SMDrSeA0C1JWViZ79Ogh//znP9c7PmbMGPnHP/5RSinlunXrPJXrpKQkWVFR4ZN8NTU1nimH9Jjf6oUXXpCAvPTSS/06v6yszNPmc+TIkSBL1yJty8Wq8wkZTqdTZmRkSEB+9dVXPp/vj4JpCvLggw9Kg8EghwwZ4vkAMioqyhMZstvtMjEx0aMk//nPf3yST1OwwYMH+/zfvKG0tFRarVYphJAnTpzw+XxNwS655BIdpGuRNudihRyj0ejpmuCPm/Xyyy8DysAeX92X//mf/8HlcvHtt996PlJKqqurSU1NBZRI14wZMwC47bbbWLRokV/yaYOJgk1CQgJXX301Ukq/up7oLV9YEURtCyna7B6+thHUdV98aUtpqatJQ/xtI6jbZUVP9+XVV1+VgLzooot8Oq9uz+BgtKX4QcSCeMPo0aPJzMzkxIkTfPzxxy2foLJhwwaKi4u54IILdF1bY+LEiaSlpbFv3z62b9/u9Xlr166lurqaSy65hO7du+smn9Zz+auvvuL777/3+rxXX30Vl8vFlClTPBazPdNmFcRgMHi6Szz11FNen6el1cZq64XZbPaMVw9H+WJjYz0u0ooVK7w6R0oZMvnChiCao5Bz4sQJj7n3JqZ/+PBhKYSQVqvV5x6tvrpYUkp54MABTw/Vhl1wGmPPnj0SkAkJCbqOH9fYsWOHBGRycrKsrq5uMf2XX37p6VkcimmOmiDiYnlL586dmTp1Kk6n06vOiytWrEBKyaxZs0hJSdFdvr59+zJhwgSqqqp48cUXW0z/5JNPAvDLX/7Sr4FbvjJixAiGDx9OUVGRV20imny33nqrXwO32iRB1LZWQZsIuWPHjs2O5ygtLfVUzv1ZAMcfCyKlssAN6uRtzY3nyMvL81TOQ7m2hzYX2ODBg5sNdhw5ckSaTCZpMBjkjz/+GDL5GiGkFqTNK4jb7ZbDhw+XQLPtDQ899FBAjWP+KojT6fT0+Xr++eebTLdkyRIJyKlTpzaZRg/sdrvs1KmTBOTbb7/dZLpFixZJQN54440hlK5RIgriK2vXrvVMLN3Y4pYFBQUyOTlZAvL999/3KW+73S5zcnLkvHnz5H333SdzcnK86rtVF62DYK9evRr19U+cOOHpv7V161af8g4Gjz76qATkBRdc0KiVy8nJ8SyGE+DyacEgoiC+4nK5PAu23HHHHecc/9WvfiUBOWHCBJ+7pixbtswzoZzW+rx69Wqf8nA4HHLAgAESkEuXLq13zO12y9mzZ0tATp8+3ad8g0VVVZXMzMyUQKM9lq+88koJoV/qoQkiCuIPu3bt8gxWWrNmjWe/VnqbzWb53Xff+Zxvbm6up24Ayvp//ixtpk33bzAY5IYNGzz7tdI7JiZGHj582Od8g8X69eslIC0Wi9y0aZNn/wMPPCBBWewmkN6/QSSiIP6iTe5sMBikzWaTt912m2eRyyeeeMLvfBcvXuyxIsuWLfM7H21yZ7PZLO+88055yy23eBTvpZde8jvfYLF48WJPv7J77rlHzpkzxyPfunXrWls8jYiC+Ivb7ZZLly6tNzRXCCH//ve/BzSmIjc3VxqNRhkVFRXQwpgul0v+9re/rSef0Wj0uTOjXtTW1nrcUe1jsVjkqlWrWlu0uoRUQYSUMmgR42BlFCi7du1i7dq1GAwGZs+ezeDBgwPOc8aMGXTr1s2zwE4gbNu2jXXr1mG1Wrnhhhvo169fwHkGk40bN/Lee+8RFxfHjTfeSM+ePVtbpLoEb3CMNxdrjwoSoV0TUgVp0y3pESLoTURBdODRRx9l7dq1CCE8n2nTprFp0yb69evn2ZeamsqmTZuYN29evbQrVqxgxYoV9fbNmzePTZs2kZqa6tnXr18/Nm3axLRp0+qlXbt2LX//+9/r7bvnnnvYtGkTQgg6derU2reozRBxsXRACEF2dnZri9EkQ4YMIYjPPdREXKwI+jJmzJjWFqHNECazcwWX40VVvLbzBPnlNYzK7MDPhnbGaAhpwdMsp8pr+eTHSspr3AxOj2Jst2gMQZy5pCUeeuihZo/n5Jbzxq6TVNTUcmmfNK4YkB7UmVXaEu3Oxfr8QB6LXtpFpcPl2TcqswMr544iIarpdbuDyT333NPkgKItx6v415YiHK6zt2tYxyiWjEvBagqNQV+4cCFbtmxp9NjaXSf4/Ru7qa0j35SB6fz758OxhEi+Foi4WP7yQ16FRzkm9U9nydXnk55g5esjxSx6aRdud2h0+Nprr210/6EiB49sLsThklzaLZpbhiaSaDXwzRk7/9paFLJ6wdatWxvdv/1wEfe+rijHzOFduHdKPxKiTHzwXS5L3twTEtnCjXajIFJK/rRuL5UOF9OGdGLFzSNYMK4Xry+8mORYC1/kFPDS9mMhkeWyyy47Z59bSp74uhinG6b0juWesSnMHJDA3yedR4xZsPV4NRuPVIVEvsaodbm5/809uNySBeN6suz6ISya0JuXbxtNlNnA6ztP8Mm+3FaTr7VoNwqy82gxm38oJCHKxAPXDvT4zF2TY/jrdGVyhsc+PkhFTXCXYvCWb0/bOVjoIDnawLxhSR75uiSauXWYsqrui7tL67k2oeSD786Qk1dBt+QY7pp8duGwQZ0TuecKpaX/H+/vD5kVDhfajYK89JViHW4a050OsfWHg141qCNDuyZRUOHg9R3HGztdd9YfrADgmr7xxJjr3/aJPWPpmmAir9LFpqP6W5GNGzees+/5LUcBuO3SHkSZ60/FetOY7nROiuZgbgUffv/TsiLtQkFqnC7e33sGgBtGdTvnuBCC2y5V+hO99NUx3X39hmHUqlo3u07bMQiY3Cv2nPRGg2BG/3gAPvihQlfZANatW1fvd16Zne1HirCaDMwY3uWc9FaTkfmX9ADg1a9D46aGC+1CQXYcKaa61sX5HePpmhzTaJorBqaTFm8lJ6+C3SdKdZWnYRj129N2XBLOT7WQFNX4RNlju8UQbRLsL3BworRWV/mWLVtW7/fnB/IVGXqnEmdtPPI/Y1hnLEYDGw/mc6bUrqt84US7UJAvcgoAuKxvWpNpzEYDVw3qCMCH35/RVZ777ruv3u/sXOWFGtGp6TUQo80GxnRVjm87Ua2fcI2w+ZBy/8b3a/r+dYi1ML5fGlLCRz+hynq7UJA9J0sAGJmZ3Gy6KwYoCvLBd/o+4IZh1B8KFYvQL7X5qXIu7KwoyPaToVWQPScVizqsa4dm003qnw7wk4pmtXkFkVLy/akyAAZ2Smg27UU9k4m1GPkhr4LcstC4CU635EiJA4BeHZpXkKEZUZgMcKDAQYXDrZtMDz74oGe7osbJ4YJKLEYDfTs2PxfXhPPPA2DLoUJqnK5m07YX2ryCnCmzU1xVS1KMmYzEqGbTmo0GRqhWZvvholCIx6kyJ7Vu6BhnJNbS/O2OMRvonWxBAgcKanSTqe4ArQNnypES+qTHYTU1v5BQWryVvulxOJxu9p7Utx4XLrR5BTlSoIRFe6fFedVf6MJMxY34+oh+ClI3jHqmQml3yYjzrptL/zQrAPvyHcEXTGXmzJme7eNFyv3LTDk3utYYI7orBcyOI8XBFywMafMKcrJE8dc7d2i6AlwXrZ7yzbES3WSqG0bNrVQUJD3Ou2Xe+qv1lH06WpC6aArSJdnL+9ddKWB2Ho0oSJvgZLGqIEnePeD+GUo95WBuOS6dWoXrhlFzKzQF8a7jdO8URUGOFNeGpG/W8WJFQbp2aDw83pALuiQCsO9MmW4yhRNtX0FKlAfcyUsFSYw20zkpmhqnmyOFlXqKBkBepVKZ9VZBUqKNxJoF5Q43xXZ9KupTp071bB8vUgqYptqPGpKZGovZKDheVE1lK3XbCSVtXkFOq41W3loQgPM7Kq3W+0+X6yJTXUrsioJ0aKKBsCFCCLonKfWVoyX6NBjee++9nu3CSsWVS4uzenWu2WigV5oS7crJ07/Vv7Vp8wpSWKFUZlPivJ+Ov6+qIAdz9VGQemFUNVwbb/X+VndLVBTkmE4t6rfddptnu6RKuUaHWO/HyvRNV+/fGf0LmNamzStISZWiIB1ivFeQ7qo7ofnfwaZuGLW8RlWQFkK8demcoLysWgQs2Bw8eBBQ2pA0BUmK9v7+9UhVIl7Hilqve36oaPMKUuwpAb1/wJq/faJInxZrLYwqpaRctSBxPihIWqzijuXppCAaVQ4XDpcbq8lAtMU7FxCgixoxPKFTARNOtGkFsde6qK51YTYKYn14wFrERi8LolFVK3FLiDYJzEbvR4qmxyoVeq2CH2y01bVKqtXCxQfrC9BFvX8nikPbJaY1aNMK4nEPYiw+TSqQkRSF0SA4U2bXtctEuR/1D4DzNAtS6dQl1Lt27VoAiisV9zQpxrex+mctSERBwppitf6R7GMJaDYaSIuzIiXklwe/QU4Lo5bXKMrnS/0DFHcs2iSodkpd+mQ9++yzQN0CxjcFyUiMwiAgt9yOw6lfn7FwoG0riJ8lIEBqvKJUBRXB79KhhVE9FXQfLYgQguRoxYro0RayatUqAEqqfQ9wAJiMBlLUAqaoUr8uMeFA21aQKv98aIBUNe5foIMF0cKoZ10s7+tHGolRyqMptevnAhbXcVF9JSVWK2BC0yWmtWjTCqKVgH5ZEFVB8nV4wFoY1Z8Qr0ai2rBYqlNrOkBJIBZYK2AiChK+lKpRmMRAHrAOFkTD30o6QJJqQUp0sCBPPfWUkrcniuX7/dMaZgt1cFHDibatIH40cmmkxetXAmph1IAsiOqWldboZ0G0IId/LpZy/7SuKu2Vtq0gmgWJ9qME1HxoHSqZWhjV30o6QJJaSdfDgixYsEDJ21PAhFeQI5xo0wqiPWB/FERrede6qgQTLYxapipIgh8KkqBanfIQWJBkH3ohaGiBkeJIFCt80SyIP5VMze8urgx+h0AtjFqmtoP4oyCa1SnTUUG0EK0v3XQ0PAqiQwETTrRpBSkJwMXSHrAeFkSj1GNBfA/zakqlhwWZO3cucFZBUvxSELWAqdJ3Dq/Wpk0rSKFawfarBIzVSkD9HnAgLpZmQcodwa+DzJs3j1qXm3K7E6NB+LUsxNn7F7EgYYnD6Sa/ogaDgPR47wb71CXWYsRsFFTXurDXBvclfOqpp6h0uHG4JBajIMrk+5IWWuSrrMYd9P5Y1113naeLTYcYCwY/Fhc6a4EjFiQsOVNqR0o4Lz4Kk9H3vyGE8IQ39SgFT3tmMzH5tTqT1WTAahQ43VDtDK6CFBYWeoYbd0/xbqhtQ7R6X0mVo13P+B7QEmyZ960fCTwNZydDyC2zEx9lYuO9EwKXrhm0SQN6pnk3XU1jdIgxk19eQ3FlLRmJ3klwceoAABYcSURBVA/ZbYkFCxbw6FvKCk4Z8f7f4nirgZoqF+U17nNmhA+UnFxluKy30/00xGw0EG81UV7jpMxe61dbSlsgIAU58tA1O4Ch6k/5Q145s5/cyv1X9w9cshbQJlwe3q356TKbQ8+K+g51+tD+af6/OPFWAwVVLsodbtKDJRjQt29fPjuQB8DITP/vX1KsmfIapzpxX/tUkKAVSyVVDm5dtYNfXdqTKwZ2DFa2jfLp/lzWqOt8/Gyo/2t+a27Copd3eULGwSB2wHg2H6/GIODirv65MABxqtV4bGsRVbXBiWZJKfn50uV8fiAfq8nA5AH+q16sRSlfl761lypH8EY/Ol3h04U+KIt4Zt633nRxr5Ta8+KtPHrDsHrH5vx3K063REpldJ0E8GxLpETZlsr22X1STXt2W0rl5p1SZzL5zfhe/P7K8/2We8j/fniOYtw1uS8LxvU8ZxGZpii313K4oJLDBZUcyqvgswP5nsmgrx+YwC+HJPot37Uvn7vYz63Dkrimb5zXIxQrHG5OldVystzJibJavj5p54g6W8qSq89nwbhefsuXed/6c/Y9OGMws0Z08XrBz5IqBz8WVHI4v5JD+RVs/qGArskx/OfG4U2dEtJFPIO1DPRjVQ4XD8284JwDO44WB32CNpNBYJvQmzsn9gkon+4pMeesFfKvjw7yr4+U3rjdkmNIjDZjNgpMBgMOlxuH0+35rqxxUthIS3KMCeYO68CU3v7Xj0BxsRq2gzzzTQnPfKPMCtkp3kSMWWAyCAxC4HRLal2SWvW72ikbbWh0VZXy2NxxzBjWOSD5GmPJm3s8C372TIslzmrCbDRgFIIa7f45lbHw5XZno1GwkyV2XG4ZFkt3B2xBMu9b/2vgD9vvn9jlvPhzJ4/W5sBV/qtACDAIgQDPtnJcOdbwuPBsCwwCBILEGLNfjYMNOVpYyWUPfx5QHlaTgR6psfRIjaVnWiyDOiVy06ThfPrRB4HLV+Jg8YbAlhqwGgWdEkx0jjfRKd5MnxQLC669BOkM3KX89ngJ05dvDiiPWIuRHmmx9EiNo0dqLBd0TmRs79TmJpEIqdYEpCCZ960fB7wFTDry0DU7gyZVK+J2S0qqaymoqKGsuhaT0YDT5cbplpiNBiX8ajJgMRmIMhtJi7Oe044ghCA7O1sf+aSkvMZNUbWLGqfEIMApFblNRoHZoEwQYTYo7S+JUQZPIaQxZMgQ3aY1dbklxVUOzpTacbjcGIWgVr1/FpMBi+ceGom2GEmN820+AdqYi/Un9fuZqx77ot6BdxdfEhYm0lcMBkFyrMWvDnyhwCAEiVFGz4Aqf3jjjTeCKFF9jAZBapzVM96mrRNomHdinZ/tt7XIR/r27dtyolbkwIEDrS1C20GJHgXlE0Fl7ty5EqXAkIB85ZVX5CuvvFJv38KFC2V2drZMS0vz7Ovfv7/Mzs6WM2fOrJf2o48+ko899li9fUuXLpXZ2dn19o0bN05mZ2fLcePG1dufnZ0tly5d6vndsWPH1r5FgRDMd7bFT1DCvJquBSujts6mTZtISkpqbTGapKSkhHHjxrW2GP4SUr+9zfbFihAhFASrHSSsqCor5ctXnuPwtzuxREVzwaQrGXbVNAwG/yu2wcReVsoPn39I4eEfsETH0GXEaLoMv9CvTo16UJqXy6aXV3Fy315ikjow4uqfMWDc5WEjXyhpdy5WZUkxryy9h9K8+u0H54+9jKtuvyskStKci1VdUszXz/2Xmor6KzR1HjaK/ldND8lL2JyLVXjyOKv/9Huqy+vLN3LadVz2y1t1l80LIi6Wv0gp2fDvRyjNy+W8Hr246R+Pc+3dS7BER7N/80Z2rn+7deVzu9nz5qvUVJSR1KU7YxbcyaDpczCYTJz85mtOfvN1q8rncjp59/8eorq8jO4XDGPuv57gioV3YDAa2fHOWvZv2dSq8rUG7UpBfti+lWN7s4mKi+e6+/7MeZk96XPhxVy9WJkKdMvqFykryG81+U7t+YbSU8exxiUwdM7NxKWlkzFwCAOuuQ6AnE/fw1Gl/7JwTZH94XoKjh8lKT2Dn919PymduzJ4whVcPu/XAHz67H9xVLf/JQ/q0q4U5Ot3lel2Lp59I7FJZ7tx9xpxIX1HX4Kz1sHX6/RrJGsOKSVHtymNqb0nXIE56uz4k44Dh5DSsw/OmhqObQ+s64a/uF0udm5QLOxlN83HHHW229AFE68ko08/qstKyf7ovVaRr7VoNwqSf/Qwpw/uxxoTy6Dxk885PmbmDQDs/fRDalqhlC45cZTKgjwscfF0HFi/U6cQgp6XXg7A8Z3bcAWhn5SvHMneRVl+Hh0yOtFrxIX15TMYGDPrRgB2bXgbt1u/+YLDjXajIDnbtwLQb8yl9Uo/jdRumXQZMAhnrYMDW78MtXjkHfgegIyBQzAYzw0eJnXpTnx6J5x2OwU5+0Mtnuf+9b90AsJw7muROWQ4SR0zqCgu4tgeffqZhSPtRkEOqxXcniNGNZlm4GWTANi/eWNIZKpL/sF9AKT1bXq0ZacLlLE0p/d+GxKZNKTbzY+7tgPQe+ToRtMIIRgwTrFy+778PFSitTrtQkFqqio58+MPGE0mug0c0mS63iNHI4SBk/u/o6YqdJXNmopyqosLMVosJHXp3mS6884fBEDR4UO4naFbg7z4zCmqSkuI7ZBMarfMJtP1vegSAA5/uxPpDp9Rf3rSLhQk7/AhkJLUbj0ada80ouLiyOh7Pm6Xi2MhLKVLTyojAxMyujTqvmhEJSQSl5aOq9ZByYmjoRKP0zlK58WM3n2bbYdJ7tyF+NQ0qstKyTvyY6jEa1XahYKc+fEHANJ7tjx8tMcQZSjn0RD60WWnTgCQ2Llri2lTeik9gQsP/6CrTHU5c0gZQdmxd79m0wkhyFTv35Hd3+guVzjQLhQk7/AhANJ79m4xbad+AwA480PounyX550BIKFjyxNMaC5Y2alzx6PrRf7RI4B396+z5/4d1FOksKFdKEjx6ZMApDTj32t07NUbhCD/6GGcjtBMm1lVXABATHJqi2kTO3UBoPTUSaQMjZ9fcuYUAMkZLY9R79hbsXCa1WnvtAsFKck9DUBSesvTDVmiY0jt0g23yxUSP9rtdlFdXAxATHJKi+mt8QlY4xNxOWqoLCzQWzwc1VVUlhRjNJuJT2lZgZMzOmOJjqGiqJCKokLd5Wtt2ryC2CsqqKmsxGyNIibRuzEYaZk9ASg8cUxP0QCwl5Yg3S6s8QkYzd4N443vmAFAZX5gEzZ4Q/FpxXokpWc0G0DQEAYD56n3r+B46AIJrUWbVxDNeiSmd/S6J2yy6sYUqZVnPakuVmZ1ienQsvXQiE1JA6AyBP3GStX6UZKqlN4QyvvX2rR5BdG6tSee5/1sjsmd1Qd8Uv+KcE2lMgeuNT7e63M8ClKkv4tVWaK4f3Edkr0+x3P/Tp3URaZwos0rSFWp9oC9n2PWUwKe1L8EdKgKYon1Q0FCYEGqSpVJ6Lx1TwE6dFIq88UhjLS1Fu1AQXx/wInnKfPRlhXk694ifFZB4rw+R6vMV5cU6SJTXSo998/7AqaDGq4uydW/jtTatHkF8ecBm61RRMXF43Y5qSorbfmEAPBHQcwxMQijEae9GletvqForYCJ9aGAiVMVuLK4sN13OWnzClJVqrzgvjxgwBPS1DtU6Y+CCGHAGqest2JvMPQ12FSV+G6BzdYoomLjcDmd5wzNbW+0AwVR6iDRib7Noq4pSLnObQ01FWol3QcFAYiKVxSkRucX0GOBfZymSLMi5e28LaQdKIjvLgLUVRB9K8L+WBAAa4L+CiKl9Pv+xYXIArc2bV5B/KmDKOmVF6KqTM8X0E2tOnrREuvbUgiai6WngtTaq3E6ajBZrPWGAHuDFhauLNY/kNCatGkFqbXbcdbUYDJbsET79oCjNB+/QscXsKoKKd2YoqIbHUXYHOZoZWWqWnu1HqIBdQuXJJ+nG4qKU8LWdtVCtlfatIJoDzg6MdHnBxwdp7g81eXlQZdLQ3OvfK1/AJhVha+t1k9BtAq6r+4VnFWQSCU9jNEq6H494HjNguinIDV+1j/grAVx2vUb+VjlZwUdIFqzIBURCxK2VPrRSKgRpVoQPRXE3wo64KkT6GlBwv3+hQNtWkGq1TYQXyvoANFqHSQULlYgFkRXFysQC+ypg0QUJGypVB9wjI9tIABR8ZqLoJ8PHZiCqBYkFC5WIAoScbHCF39j+ADW6BgQAkd1NW6XPhOh1ajuhzXO+46KGqYozYLopyAV6kCuWB968mpEXKw2gD/dJDSEwUBUjNI2oVeoMiAFsVpACFwOh24KXKl2hoxN8kNBVKvY3sO8PgXnl82ZOh74G/AjMAgwA7++e/W7rTKhbCCVTABrXBz2ygrsFRXEJPjuprWEQ1UQix8KIoQBc1Q0tdVVOO3VfrlpLVGpWZAk3+tw5qhoDEYjzpoanLW1mMyBL8vdGKWnT3J0xzai4hOoKi6isqhgL/DrsfMXheSd88eCXAQsu3v1u8OAZ4EHgyuS95SqQ1IT0tL9Oj9KHaNRo0MpKKXEXq4EEfyxIAAmLZJltwdNLg23y+UZLBXrw1gaDSGEp31Hj/tXl/L8XDoNHsrQGXMgxO+cPytMHb179bvarGu7gLnBE8d7amvsVBQWYDAaSUhN8ysPjx+twwN2VFbgtNsxRUX5Xfqbo6KpRp+KeknuadwuJ/GpaZgt/i3ZHBUbR3VZKfaKCr+skNfXiYsnLsXzjP1+5zavXJ4MfAi8CjyjbY+dv+iRps7xx4LUjTtKQrzij0buIWViteTOXTEY/Vs1yqqjH12mTkUUl5ru96pRWiTLqUOoN/dQDgCpXbr5nUeoKuoNuukE8s4NBD5SFaLudpMEZY3CZXOmZmb0PR+zxUqt3c41d/7OM2pPL77/4lMAujVYSsAXotQOhDU6hCq1CaiTe7Q8GVtTnHWxgqsgUkr2bf4cgO7qhNn+0JoV9c0rl2cCLwNVQBxww9j5i440SNMJWAVEA68DPwN6bV65/FPgf7TtsfMXfdTUdYK1iOf/G3/TfDr26stz995e70CN2mEPCRJ5dgRanbWo620jkW71W6Keq2yjnp+zfSt7Pv0QIQxcMOlKv4XWLMgnzzxBQtp5ZA4dHvAahtUV5Rzf8jm53+9GGI10Uqfq9AeTVXF99r69BktMLMk9eiFEYIFHR2UFxzZ9TN7ebzBbo+h/yfiA5Xvrnw9w/R8fpMuAwQGtseiormLvhrcpyT1NWrdM+k1s9tn+P+Bu4GtgdxNp7gP+NHb+oq2bVy5/E8gCho6dv+ijzSuX1wKTmlMO8FFB7l797uco0at6v5fNmbohrXsPDEYjad171Dtn+fwbdBuWOWHubaQE4CLUdQ3e/Mf/AhCfmsaEW24jNqkD1phYLNExWKJjMBgNIMHlcuJ0OHA6HFSXl1JeWEB5QQHlhXmcOZTDmR8O4nI6QQj6X/kzov1o5ffIp0bpAHa98iwAcWnp9LpsMpa4OEwWKyZrFEaLFYM6p5WUbtxOJ7U1dhwV5dSUl1FTUY69tITSk8cpO3MS6XZjMBqZ8ps7/Y4AAhQcPzuv2JoHlgDQqW9/Lpw+m5iERCzRMVhjYrBERytzbkllHUSnowanw0FVWYl6//IpK8jn0NdbKSsqJKFDMtaoKA588j6Dp17HsJk/91xn7PxFnwODNq9cvgHIHjt/kXPzyuVNTbTcB/jH5pXL3UAHoOXJkRsQLAtyOP/oETr27kNhg8nErDGxuF0uhBBK6aJ+PL9R2iQEqPsNIKiXXiDUfQYQgqi4OEZfN4eew5peC8Qbhky6ij2ffFBvX3lBPuuWBRAkEYKErpn0u/xKkgJQXoDMiy+j8Mecevsq8nPJfv3FgORLzOzFtF8v9mou3ua49MZbWPfI3+rtO3VwH2/984GA5BNCIN3ulqY9OgxcsHnl8q+pU2g34EfgybHzF+3ZvHK5DTgNeD9BGcFTkIc3vfSMzWSxIqWsV2letPKVIF0i+KT37M3tz64m56stFJ06wf4tm0hKz8BgNFJTVYmjqgpHdRU11dVIlwuEwGA0KCW3xUpUbBzxqanEJacSn5JKatfudO43gO27djW5DLQvJHfvyfi7lpK7fy9Vhfmc/i6bODWk7bTbcTlqlHYIR41ipdWCxWAyY7JYsMTFY41LwBqvfMenZ5DUpRsV1faAlQOgz6gx/GbFS+R8tZmiUyfZ9+XnpPfsjdvloqayEoe9GkdVJY7qalU+gdFkwmSxYLJYiIpLID4llfjUVOJT0qg8fRIDEqPRCEIQ2/xcxg8DL6AEjQyAc/PK5QOAa8fOX/SQmuYh4OnNK5cnAjuA7339j0FZJ33ZnKnTFmStWhfXIZkXfn8HN/7tX5gs3k2z2R5pbp30cKC5ddJbE0d1FQc+eZ/KogJik1PpN/FKLGqnzToIgM0rl08DvgFOqd8XqccWjJ2/6LFgyRQsC5L39iN/xe1yccHkq3/SyhHBfyzRMQyeep23yfOAt1De4SfHzl9k37xyeTywIpgyBcWCqAQto7ZOxILoSkjb3YJlQSLUIT4+npKSkpYTthLxPswT/FMnYkEitDVCakHadHf3CBH0JqIgESI0Q0RBAuDgwYP85S9/AeCdd95h27ZtPp3/4IMPUlVVRXV1NY8++qgeInp4/vnn+eijZntVRGiESCU9SEybNs3nc5YsUbpnFBYWcvRo+1/OrC0SUZAg8fzzz5ORkcHkyZO54447mDhxIvv376empoZrrrmGXbt2cerUKRITE/nNb36D1WrFZrPxz3/+k+effx6Hw8GDDz7Ifffd5+lXBWC321mzZg0//vgjBoOBIUOGcOWVV7JkyRJ+97vfkZ6utKw/9thjjB8/nn79+p2T/tprr60n6+nTp3nttdeorKxESsn48eO5+OKLQ3q/2grBjGL95LDZbOOB/2RlZQ2y2WyrgL1ZWVmP2Gw2CdyZlZX1uM1m+z1wP9AfpS/Q18CyrKysl9V0aSjdtfdmZWWdM7LKZrP9C8gAfgkYgY+APwHTAUdWVtbvbDZbL+AzoAdKF4zG0s8F9gKPAtnATVlZWbtsNlsisBW4NSsryzcf8SdAxILoxxvq9yFgT1ZW1kkAm812GPBlloRJwF1ZWVkuwAVcpuZzCthks9nuBxYAT2dlZblsNltT6eeq+fUFegHP2Gw27RrRwDAgoiANiCiIftTU2a4NIB8nddqYbDZbV6AqKyvroM1m240yCOhGlL5ITaavk58RKM3KyhpaJ006oO9SW22USBQrPHACRpvN1lgj2MfALTabzWCz2awoI+MuU48tR3GptmdlZZ3yIj3AAaDaZrP9EjwKtBcYEew/1R6IKEh4cBrYDnxns9kajlf4X8CBUm/4BtiQlZW1Vj32Lkr95Ukv05OVleVAsTq/Ui3Qh8DSrKysVpm6KdyJVNLbMDabbQzwNDAoKysr8iB1IFIHaaPYbLbngPHAnIhy6EfEgkSI0AyROkiECM0QUZAIEZohoiARIjRDREEiRGiGiIJEiNAMEQWJEKEZ/j/RqZN8edzi2gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 198.425x198.425 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run mechanical simulation\n",
    "Qm = -71.9e-5  # C/cm2\n",
    "data, _ = bls.simulate(drive, Qm)\n",
    "t, Z, ng = [data[key].values for key in ['t', 'Z', 'ng']]\n",
    "\n",
    "# Create figure\n",
    "fig, ax = plt.subplots(figsize=cm2inch(7, 7))\n",
    "fig.suptitle('Mechanical simulation', fontsize=12)\n",
    "for skey in ['bottom', 'left', 'right', 'top']:\n",
    "    ax.spines[skey].set_visible(False)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "\n",
    "# Plot variables and labels\n",
    "t_plot = np.insert(t, 0, -1e-6) * 1e6\n",
    "Pac = drive.compute(t)  # Pa\n",
    "yvars = {'P_A': Pac * 1e-3, 'Z': Z * 1e9, 'n_g': ng * 1e22}\n",
    "colors = {'P_A': 'k', 'Z': 'C0', 'n_g': 'C5'}\n",
    "dy = 1.2\n",
    "for i, ykey in enumerate(yvars.keys()):\n",
    "    y = yvars[ykey]\n",
    "    y_plot = rescale(np.insert(y, 0, y[0])) - dy * i\n",
    "    ax.plot(t_plot, y_plot, color=colors[ykey], linewidth=lw)\n",
    "    ax.text(t_plot[0] - 0.1, y_plot[0], '$\\mathregular{{{}}}$'.format(ykey), fontsize=fs,\n",
    "            horizontalalignment='right', verticalalignment='center', color=colors[ykey])\n",
    "\n",
    "# Acoustic pressure annotations\n",
    "ax.annotate(s='', xy=(1.5, 1.1), xytext=(3.5, 1.1),\n",
    "            arrowprops=dict(arrowstyle='<|-|>', color='k'))\n",
    "ax.text(2.5, 1.12, '1/f', fontsize=fs, color='k',\n",
    "        horizontalalignment='center', verticalalignment='bottom')\n",
    "ax.annotate(s='', xy=(1.5, -0.1), xytext=(1.5, 1),\n",
    "            arrowprops=dict(arrowstyle='<|-|>', color='k'))\n",
    "ax.text(1.55, 0.4, '2A', fontsize=fs, color='k',\n",
    "        horizontalalignment='left', verticalalignment='center')\n",
    "\n",
    "# Periodic stabilization patch\n",
    "ax.add_patch(Rectangle((2, -2 * dy - 0.1), 2, 2 * dy, color='dimgrey', alpha=0.3))\n",
    "ax.text(3, -2 * dy - 0.2, 'limit cycle', fontsize=fs, color='dimgrey',\n",
    "        horizontalalignment='center', verticalalignment='top')\n",
    "# Z_last patch\n",
    "ax.add_patch(Rectangle((2, -dy - 0.1), 2, dy, edgecolor='k', facecolor='none', linestyle='--'))\n",
    "\n",
    "# ngeff annotations\n",
    "c = plt.get_cmap('tab20').colors[11]\n",
    "ax.text(t_plot[-1] + 0.1, y_plot[-1], '$\\mathregular{n_{g,eff}}$', fontsize=fs, color=c,\n",
    "        horizontalalignment='left', verticalalignment='center')\n",
    "ax.scatter([t_plot[-1]], [y_plot[-1]], color=c, s=ps)\n",
    "\n",
    "figs['c_mechsim'] = fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Panel C: cycle-averaging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJsAAAGACAYAAAC6BxBgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2deXxU1d3/32ey7ySEJBCyEMKubIrsAVTcqoLiVhWXarCWVh+tT7V9+tP7dLO0im1tfGzjLrXWUhXFXUQSBQUEWcOaBAg7IWRfZ87vjzsTh8xkMklm7iw579drXvfOPfec872TT75nvecIKSUKhRGYfG2Aou+gxKYwDCU2hWEosSkMQ4lNYRhKbArD6LNiE0KECCEeFEJsFEJ8K4TYKYRYIoSI6EFaLwkhHvKGnUYihPiVEOI2b6Uf6q2EA4D/AxKBi6SU1UKIGOAfwHPAQp9a5iOklI96M/0+6dmEENnALcBdUspqACllPfBD4EMhRJUQYrjd/Z8KIeYJIWKFEC8KIfZYPeHvhBCiQ9qjhBAfCyG+sXrMH3Rig0kI8WchxNfWtEqEENOFEAlCiBohRJrdvV8LIS4XQoQLIZ4SQmwSQmyxetR46z3lQoh/WdO5RghxpRBirdVzHxRC/NouvUeEEHut6fxJCFFuvd7uoYUQTUIIzZpGmRDiXuv1ECHEUiHEPuszPiOE+NytH15K2ec+wAJgvYvwPwF/sJ4PBQ4CIcBS4J/W83BgDTAbeAl4CL2k2AFMtMZNAHYCU5zkMRX4N2Cyfn8EeNd6/jLwkPV8FHAA3TE8CvwRENaw3wHPWM/Lgf9nPRfAamCY9fsgoA1IBi4FdgH9rPc9D5Rb73vJLl8J/Nh6fh7QBEQC91ifO9L6G3wEfO7O794nPRtgwbVXfwa4TQgRBiwCnpNSmoGLgeellGYpZYuUcpaU8nO7eMPRxfmCEOJb9D9KFDChYwZSynXAL4F7hBBPANcBsdbg54Dbred3Ai9IKS3AlcA8YLM1/fnAaLtki61pS+Aq4DwhxGPo/yQCiAGuAP4tpTxjva/Axe+wwnrcBETYxX9FStkkpWwB/uYi/ln01Trb18AoIUSclLLWdlEIkQ78Hf0PvxX9D3szMNl6Sxv6f7zt/gygwS7dEKBaSjne7p5UoFoI8Svgauvld6w2/Bl4Ev2Pugu4FUBKWSyECBVCXGDNf6pd+vdLKT+wph2L7mFs1FmvxwCbgbfQBfgCujCF9Rnsi36zi9+p0WqPtNYWuhv/LPqkZ5NSHkFvDLxgV+eJR/dolVLKRvT/+D+iF7dHrFE/BW631rcigOXALLukdwONQohbrWlmANuB86SUj0opx1s/jwJz0YvN/wM2ooshxC6t54Cnga1SykPWax8BP7bW3UxAIfC4k0ccBsQDv5RSvote1EdY038PWCCESLDeexd2/0Bu8B5wqxAiQggRCtzhbvw+KTYrP0KvT621FklfW7/fbQ1fiV6sPWsX53+BFmALuud4X0r5pi3QWqzMA+4WQmwFPkavR33pJP9ngdlCiG3oxdR+YIhVRKDX28aji87Gr9HrZputtgrgp07S3mq1f5cQogS9SN0J5EopP0MX6TohxEb0emWDkzQ64yX032ozsNb6e7gV31bRVHRACDEV/Q99jgyiH0kIcT4wTUr5F+v3B4HJUsob3Yx/CZAipVxm/f5noElK+XCXcYPod/QYQoiX0YueG6WUX/nYHI9irS48j97Klegt7UVSysNuxk9H926p6MXyFuBeae1CchlXiU1hFH25zqYwGCU2hWEosSkMQ4lNYRhKbArDUGJTGIYSm8IwlNgUhqHEpjAMJTaFYSixKQxDiU1hGEpsCsNQYlMYhhKbwjCU2BSGocSmMAwlNoVhKLEpDEOJTWEYSmwKw1BiUxiGEpvCMJTYFIahxKYwDCU2hWEosSkMQ4lNYRhKbArDUGJTGIYSm8Iw+uoCzn1qUTqLxUJDQwNVVVXU1tbSv39/IiMjSUhI6DpyzxDOLvZVsQUdLS0tfPHFFzz88MNs3LjR7XgpKSlkZ2czcuRIhg8fTl5eHpMnTyY8PNzjNvbVlScD/qHNZjMrV67k4YcfZvfu3R5NOyYmhlmzZnHxxRdz+eWXM3LkyO4m4dSz+Xy3FR99AhKLxSLXrVsn58+fL9H/Ybr83H777XLNmjWysrJSWiwWabFY2tOrr6+XpaWlcs2aNfLZZ5+VP/nJT+To0aMd0pgxY4Z85ZVXZHNzs7umOt85p7MAf/+gLx78IPoeAt+iL72+BIhwI35AsGTJEpmbmyujoqLcFtfLL78s29raepXv4cOH5SuvvCIXLlwoY2Nj29POysqShYWF7qQfdGL7O/reTwnW7zHA28CrbsT3e0pLS90W2EUXXSTr6uq8Ykdtba0sLCw8y+NNnTpVbt++3VW04BEbkA3UA/EdrqcBC9xIw+/Zvn27W0IrKyszxJ62tjb5j3/8Qw4cOFACMiwsTD755JNnFct2BJXYXO6q58anna1bt8qKigr3fnED2bBhQ5dCq66uNtyuqqoqec8997TbsHDhQtnQ0NDxtqAS2zXAxl6kIaWU8uuvv5YJCQny3HPPlVVVVd371b3MZ5995lJoa9as8al9y5cvl9HR0RKQ06ZN6yg4p797oI4gtO+qZ39RCJEuhHhPCBHlTiJDhw5l4MCBbNu2jdtvv90mZL+gtra207CysjLy8vIMtMaRBQsWsG7dOjIyMpg0aRKRkZFdR+pMhf7+4bsGQrz1ezz6VoqvABrwIvr2h78ClqG3VsfKDsVoWVmZjI+Pl4B8/vnne/qP7nGWLVvm1KNdeuml0mw2+9q8do4fP+7MnuApRqUurlD0XfK2o3d9lKBvhxhuFdtP0Md+jwFxwLXAz2UHsUkp5auvvioBGRsbK0tLS7v/i3uBZ5991kFoy5cv97VZ7hJUxShSyjYp5WNSynOkvofnKCnlz6W+DSPALqnvPrxL6hvY1qDvuenALbfcwvXXX09dXR2PPPKIUY/gko7FaHl5OQsWLPCRNZ4hYMXmBm5XwIQQLF26lPDwcP7973+zY8cOb9rlFvZiq6ioICsry4fWeIZgFlu3GDx4MHfddRdSSgoKXG2bbgw2sT3xxBOkp6f72BrPEJRik1JqUspPreezrcdPpZSaq3iLFi0C4I033qC1tdXLVrrGJrbY2Fif2uFJglJsPWXcuHGMHj2ayspKVq9e7VNb6urqAIiLi+vizsBBic0OIQRXXXUVAJ999plPbbF5NiW2IGbOnDkAPvdsgSC2mpoali5dygMPPMDSpUupqalxeb+aqQs89dRT7edms5l58+ZRVlZGS4vei+KswTBlyhSmTp1KXV0dhYWFDuEzZ87k/PPP5/Tp07z88ssO4RdddBFjx47l+PHjvPbaaw7hQujzD1taWs6yz8bVV1/N0KFD2b9/P++8845D+HXXXUdGRga7du3igw8+cAi/+eabSU1NZevWraxatcoh/PbbbycpKYmNGzdSXFzsEJ6fn89zzz3H/v37kVJSWlrKc889x4MPPuhwrw3l2ToQEhJCVFQUZrPZ4zNgu0NjYyMA0dHRPrOhKw4dOtQ+xGexWDh06JDL+9W0cCcsWLCAN998k2XLlnHLLbcYZdNZpKWlcfz4cQ4fPsygQYN8YkNXLF26lNLSUiwWCyaTiZycHJtnczotXHk2J4wdOxaArVu3+iR/KSVVVVUAJCYm+sQGd7j77rvJyckhIiKCnJwc7r77bpf3qzqbE4YPHw7A/v37fZJ/Y2MjLS0tREREEBXl1gQWnxAfH++yjtYR5dmcMGTIEECfyuMLzpw5A/i3V+sJSmxO8LXYbEVov379fJK/t1Bic0JKSgrR0dFUVVW1exkjCYT6Wk9QYnOCEMKn3u306dOAElufwSY2XzQSDh8+DBA0sz1sKLF1wtChQwEoLS01PO+KigpAn/YUTCixdUJubi4A+/btMzxvJbY+hs2z+aIYteUZDLNz7VFi6wRfiU1Kyfbt2wE455xzDM3b26ix0U5oaWlp771vbGz0ynplzjh48CBZWVkkJydz4sSJ9tkfAYYaG+0O4eHhZGZmYrFYDPVua9asAWDSpEmBKrROUWJzwYQJEwDYsGGDYXl++OGHAMydO9ewPI1Cic0FU6ZMAWDdunWG5FddXc3bb78NwJVXXmlInkaixOaCGTNmALq3MaJuu3TpUhoaGpgzZw7Dhg3zen6G09mr8kH+cQuz2SwHDRokAfn555+7G61HbNy4UUZEREhAFhUVeTUvAwiu5ReMwGQykZ+fD8AjjzzitXdJ161bx2WXXUZzczOLFi1i5syZXsnH53SmwiD/uM2ZM2dkWlqaBOT111/v0QX4Kioq5H333SdNJpME5OWXXy6bmpo8lr4Pcfq7+/qP7vdik1LKtWvXti9knJycLB9++GH51VdfdWf1bCmlvlTotm3b5N/+9jd5+eWXy7CwMAlIk8kkf/rTn8qWlpbumuavOP3dVaeum+zcuZNFixbx5Zdftl+LiIhg5MiRZGVlMXjwYOLi4oiOjiY8PJympiYaGxupra2loqKCQ4cOsXfvXurr69vjm0wmrrvuOh555JH2bpYgwWkHoRJbdyJJSXFxMa+//jqrV69m165d3U4jKyuLSZMmcdlll3H11VczYMCAnpji7yix2eGRhz5z5gz79u3j4MGDVFRUUF9fT0NDA83NzURGRhIVFUVsbCzp6elkZGQwZMgQkpOTPZG1v6PEZkeffGgDURul2RFcg44BgupnUxiGEpvCMJTYFIahxKYwDCU2hWEosSkMIyDFpmnaak3THHbH0DTtp5qmrfCFTYquCUixAc8AP3ByPR/w/SYGCqcEqtjeAmI0TWuf+KVp2iz0ztpP3Ijv9i7FgfyxWCyyqalJNjY2ysbGRllbWyvr6upkTU2NrK+vlw0NDbKpqUm2trbqUzU8l7dTAnIEQdO0Nk3TCoG70HfeA1gEPKNpWtAMRdXW1tLY2IjZbAb01wvr6+tpampi9erVTreJNJlMWCyWHuUXFhZGeHg4kZGRREZGEh0dTXx8PPHx8cTFxREfH0+/fv1ISkrCZOq+nwpIsVn5O7BT07Q4IAy4FPhRdxPZtm0bmZmZJCQkeNq+XtHU1MSTTz7Z7XgdhRYWFoYQAiEEoaGh+rwyIbBYLEgpsVgsmM1mzGYzra2ttLa2njUNyhmhoaGkpqaSlpZGZmYmQ4YMIT4+vkvbAnogXtO0fwMfAzHASE3TfuhmVAn620xPPfUUQgiGDRvGrFmz/GbloMOHDztd8t4ZKSkpDBgwgPT0dNLT00lMTCQyMpLQ0FC3PZDFYqGtrY3m5maam5tpbGykrq6O2tpaampq2o+nT5+murr6rLiJiYncf//99peCciC+AH1v0QTgNttFTdM0IAvIBVYDOcBE4CZN09pXZW5ra2P06NHs3r2bPXv2sGfPHiZMmMAVV1xBWFiYgY/hSGc7KUdGRjJr1izOO+88j76lbzKZCA8PJzw8vMuNPhoaGjh+/DhHjhyhvLyc/v37u5VHQItN07TPNU3rD5zWNG1bh+BN6HW6I8AwYC7wPaBdbP379+eGG26grq6OtWvXsn79ejZv3kxDQwM33nhjj+olnsLeewwePJgZM2aQmZnpF/siREdHM2TIEIYMGcL06dPdjheordF2NE0719oS7cguTdMs1qPLzW1jY2O55JJLuOuuu4iMjGT37t1s2rTJm2Z3iW1rngsvvJC7776bkSNH+oXQekPAi80F3a6MDhw4sP1N9KKioh636jyBzbP5W8OlNwSz2HrE6NGjSUxMpKamhvLycp/ZYRObO628QCGgW6O9wOVDf/zxx6xdu5aZM2dy0UUXGWXTWdh2ubvvvvtISkryiQ29QC2Z5S62FR8PHDjgk/zNZnN7azSYPJsSmxNsG5OdOHECX3j+uro6pJTExsYSGhrQHQZnocTmhNjYWCIjI2lqamrfPttIgrFxAEpsThFCtL/fefLkScPzD8bGAQR4p643SUpKoqKiwifbCdnE5u97V7W0tLQvKRETE8OwYcNcjmoosQE7duxwuGYbgjlz5ozT8AEDBpCSkkJrayt79uxxCE9NTSU5OZnm5maneykMGjSIxMREGhsbHTb2kFLSv39/+vXrR319vdMumMzMTOLi4qitreXgwYMO4dnZ2cTExFBdXd2+r4I9OTk5REVFUVVVxZEjRxzCc3NziYiI4NSpUxw/ftwhfPjw4ezdu7e9IVNbW8vevXsZM2aMw702VDHaCbaVwjsOOhtBW1sb4P+erePskK5mi6h+tk7Yt28fy5YtIysrizvvvNMIm9p5+umnqays5N577yU1NdXQvLvDjh07zpowEBcXZ/Nsqp+tO9i8itGezWKxBEydbdiwYcTFxWEymYiLi+tyHWBVZ+sEW7dDTU0NFovFsBkg1dXVtLW1ERsbS0SE03kDfkN4eLjLOlpHlGfrhLCwMGJiYrBYLJ3OLfMGJ06cAPQJkcGGEpsLbMWYkd0fSmx9FNtOxrZttI1Aia2PYpttYdtG2whsfV7+3ArtKUpsLjBabPX19VRWVhIaGkpaWpoheRqJEpsLjBbboUOHAH1v+JCQEEPyNBIlNhcYLTbb/LnMzExD8jMaJTYXxMTEtO9pYMRUI9sYq20X52BDic0FQoj2utPRo0e9mldlZSWVlZVERkaSkZHh1bx8hRJbFwwcOBDwvthKSkoAfbZFMNbXQImtS4wQm5SSb7/9FoBzzjnHa/n4GiW2Lhg8eDCgV9699R5pRUUFp06dap+AGKwosXVB//79SUhIoKGhgWPHjnklj7Vr1wIwfvz4oC1CQYmtS4QQ7a3DvXv3ejz9U6dOUVJSQkhICJMnT/Z4+v6EEpsbjBgxAoCtW7d6/NW+VatWATBu3Lige8GlI0psbpCbm0tsbCyVlZXtvfyeoLS0lJKSEsLCwpg9e7bH0vVXlNjcICQkpH3z2aKiIo+k2dTUxDvvvANAXl5e0Hs1UGJzmylTphAeHs6+ffvYv39/r9KSUrJy5UrOnDlDWloaU6dO9ZCV/o0Sm5vExMQwc6a+OPmKFStoaGjocVqfffYZ27dvJywsjOuuuy6ollhwhRJbN5g2bRrp6enU1NTw+uuv09ra2q34UkpWrVpFcXExQgiuv/76vrKzMqBe5es21dXVPPfcc9TW1pKens5NN93U5Rq0oK9D++6771JSUoIQgvnz5zNu3LiemuHvqG277ejVQ588eZJly5ZRXV1NZGQks2fPZuLEiU6XHmhubmbz5s0UFxdTX19PeHg41113HcOHD++NCf6OEpsdvX7ouro6VqxY0d7RGx4eTk5ODikpKe0rIB0/fpzS0tL24jYzM5P58+cH4uJ+3UWJzQ6PPLSUkt27d/PFF184XU/DRkZGBtOnT2fEiBEI0Se2p1dis8PjD33mzBnKysqoqqqipaWF8PBwkpOTyczM9Ps3272AEpsdffKhDSQod3jpKX2iLPM3VD+bwjCU2BSGocSmMAwlNoVhKLEpDEOJTWEYSmwKw1BiUxhGX+3U7ZMjCLbRIgPGZ9UIQl/ly5W7eO0PxQ7XH3/7FuKTjNudWY2NBjGtLWb+6+IXXN5z5V3ncfntEz2dtRqItyPoH9rcZuG+C593696Z80dx04MzPJm92nSjr+BKaP0GxDDr2tFnXSt+u4SvP3Tcf8vTKM8WhPzp/pXs3ey46tJf19x9VuNg0+pSnn9sVfv3Je8sJLZfpCdMUJ6tL1Cxr9Kp0AqK8h1aoRPn5PDA01e2f3/46le9apsSWxDR0tTG4z940+H6nz7pfKO33HEDueru89u/f7lyl1dsAyW2oOLheY6e6YG/XkVYhOserstum9B+/tofijG3eWcdOiW2IOHUkRpaGtscrueOdW8/hT++d1v7+Zo3HTfz9QRKbEGAlJLHbvqXw/Xfr7jV7TSi4yK45Wf68hKf/nMrzY3de9vfHZTYgoBP/7nV4VpUbDhxiVHdSmfq90YweFh/qisbvOLdlNgCnJamNt5+dr3D9d8s/3630xJCcOUPzgN0ATc1eNa7KbEFOG/86Uun1yOjHZeCcIdzpmUyZEwK9TXNfPWBZzt6ldgCGHObhXXvOwri0X9c3+M0hRBceOO5AKz5z3YsFs/1fyuxBTCf/2e70+upGb17A3/cjGwSU2M5UVHDzq89t6yrEluAUr7zBG8WfO1w/eb/ntnrtENCTcy6Rh8//Xy5c0H3BCW2AERKyR9/uMJp2LQrR3gkj2lXjiAsIoSSDYc5Wu6ZnaSV2AKQAyUnnV7vNyDGY7NwY+IjmXyZvoacp7ybElsA8rdffOz0+uInLvNoPnMWjAFg/cf7aKht7nV6SmwBRmuLmZrTjU7DBg3x7CKDadmJDJ84iJamNjZ+2rsV0kGJLeAwt5mdXp8xb5RX8ptx1UhAnw3S27mPSmwBRF11U6cNgymXeWc3v7Ezs4mJj6BibyUHd5/qVVpKbAGCxSL595/Xcqz8jNPwIWNSvZJvWHhIe0Phi3dKepWWEpuf09LURkNtMy//ZnWn9aZLF473qg0zrtaL0o2r9tNY19LjdNR7o37Mjq8O8czPPuzyPvuZtt4gNbMfwyYMZO/mo2z4ZB9514zuOpITlGfzU86cqndLaGDIG+7MvFpvgBSvKOlxQ0GJzU+p7aR7oyMz53unFdqRcXnZxPaL5Ejpacp2nOhRGkpsAU7uuIGG5BMaFsLUK/ShsC9W6A2F1uY2Sr4s55sPdlPyZTmtzY7T0s9Kw+tWBgC71h5wuJY0KI6U7CTMbRb2rnec+ZCckUByRj9am9vY/81hh/ABWYn0T4+nubGVss1HHMLTcpLolxZHY10zB7Y67j0fGhWm25ESwwVzchzCvyku5+SRWoaPS3Nqf+aYVKITIqk+Wc/RvY5dFllj04iKjeDMsVqOlZ52CB8yYRARUWFUHq7h5AF9bHToiP6E3HguwgRVJ+s4tucUdVYPXHe6kX0bKxg1PdshLRvKs/kp7lSLktPjESbjVtkPiwglOj4CaYH1H+2loebsIayO3zui3oj3Uw6UnOQP97zt8h4vLQrjkpL1Ffz1oQ/oNyCG798/hcbq7wQWmxRl82zqjfhAoq2TYSl75t5s/BaSIyelk5rZjzMn62loMhObFIUp1ERsUhS55w92GVeJzU/p6kXhc6dlEhoWYpA13yGEYLZ1NkjxihJGTc/mvMtHMGp6dpcvQwd0A2FxXuFjwM1AM/AScCPwq4Ki/Pd8aZcnMLe6FpspxHc7Il1w6TBW/H09pduOc2DXSbJGDnArXsB6tsV5hVcBNwGTgHHAWGA8sNqXdnkKs9m12EJCjfdqNiKjw5j2PX0I64OXN7sdL2DFBswBlhcU5dcUFOVL4AVgW0FRfkN3EjlRUc3B3c5nvvqStlbXdTZf1Nfsufj7YwmLCGHblwco3+leJ28gi61jizIMcN2r2IG93x7lt7cv5+Xffu61xVR6SpuLYvSBv15F5ohkA61xJKF/NLMXnAPAu89tdCtOINfZPgGeXJxXuARoBO4Bxi7OK4wAfg5kAbnoxWoOMBG4qaAov32tguzRKSSmxHKs/Ayrl2/n4pvGGv4QneGqNz42PsJASzpn7s1jWff+blIzEzC3WQgJde27AtazFRTlfwj8A9gAbAc2ApsA27tsm4BZwCLgXuCXwPfs0wgLD+H6+6cB8P6Lmzhzqt4Q292htbnzYjQhOcZASzonJj6SX7/xfW74r+ldCg0C27NRUJT/O+B3dpf+ALA4r3AGsKugKN+yOK9wV0FRfu3ivMIawMEljJmSwbiZWWwpPsCbBV/zg8cuNMb4LmjpxLMt/PksomJ7trSCNwiPdF9CAevZ3MDtUYIFP55KWEQI36zaz6p/Oa4I5AucFaMPF85nyuXDfWCNZwhoz+Yp+g+M4+b/nsnLv/mcNwu+Zv/W4wyfOIiwcOfdC/ZDfN+dOrlmJ3dbnLNHBzuPs9ruXc1Jl+Ry80Mzu+VF/BE1NmrHug/28K8nv6C1peuhIqNY+PNZgejN1KYbdnT60FXH69j42X5OHanF0rFjVdL+MzqbHXvWJesX+2uiPXIXcaxfk9JimXP9uZgMnNnhIZTY7OiTD20gaqM0OwLOVQQDwdwaVfgZSmwKw1BiUxiGEpvCMJTYFIahxKYwDCU2hWEosSkMo6926npvBMFiAXMLhEZ0GIvqU6gRBK8gJRwogu2vQ9lnUFUGllYwhULySMicCed+HzJn9GXxAWpstHcc2QTv/QgOd9j8whSmC86e9Mkw9w+QneeRrP0cNRBvR+8eWkr46s/w8UMgzRCdDOfdAyOugpRzITwaWhrg+BbY/Q5seh4arG9wTX0QLnocQv1ntq0XUGKzo3cPvep/oNg6G33y/XDhbyAitvP7m+vgyyVQ/Lguzpy5cMNyiIzvlRl+jBKbHT1/6A3Pwnv3ggiBa1/V62PucugreH0e1J+AjOlw2ycQ1r0NaAMEJTY7evbQJ3fBs+PB3AzzXoAJd3Y/jdOl8NIsqKmAkdfoHs4UdD1QahWjXiElrPyhLrTxd/RMaABJOXDrhxDZD3a9BV//2aNm+jNKbO5S9hkcWAORiXDpU71LK2UMzH9JP//kYTi+rdfmBQJKbO6y9gn9OO0hiOrd5rEAjJwH5/9Q7yJ5/yfuLTUZ4CixuUP1Idj3EYSE6wLxFBf9Tu82ObAGdi73XLp+ihKbO+x4A5AwYh5Ee3Dnu6hEmPMr/XzNr4PeuymxucM+6+YXo671fNoTfgBxg+DENtgT8GsYukSJrStaGuBAMSAg52LPpx8aoY8qAHz9F8+n70cosXVFxTq9u2PgBIjx0ppo4+/U64Oln+r1wyBFia0rjm7Sj4OneC+P6CQYOR+QsOVV7+XjY5TYusImtoFe3m/g3Fv04y7Xex8EMkpsXdEutvO8m0/OxRAaCUc2QO1R7+blI5TYXNHWDJV79UH3AT3bY9NtwqNhyEX6+d73vZuXj1Bic0VVGSChX5Yx88+GXqIfDxR5Py8foMTmitP79GNSrjH5ZU7Xj4fWGpOfwSixucJosaWOg7AYPd+648bkaSBKbLp+3UYAABUTSURBVK6wiS1xqDH5hYRC+gX6+eENxuRpIEpsrjDaswGkWXduObHd9X0BiBKbK86U68fEIcblOUDf8U6JrS8hJVQf1M8TsozLN0XfokeJrS9RfxLaGvXp20a+BZVi9Wyndulv1wcR6o34zqg+oB+N9GoAEXEQk6K/gVV3FOLTjc2/G8iWWuSOF6HuMMSmI8bciQiP6/R+JTbAsvlph2vCrL8iJPtlI52Fp12AGDgZ2VKn/+Adw9OnI1ImIpuqkCXLHMMz5iCSz0E2HEfufuPswDFXIPauQpwpRwqJ3PeWY/ycKxEJQ5DVZcjSlY7hudcg4gYjT+9GHvjYMXzEDYjoVOSp7chDjlu0ilG3IiITkSc2IQ9/6Rg+5k79uautHd/V5cgdLyIm3Odwrw1VjHZGY6V+TMg0Pu9Q67ukVWXG590d6g7z3VuRFuv3zlHvjXbG+/fB+qfhkidg2k8NMMmOTx7R36Cf82uY9Utj8+4Gls1/gepywAKYICEbk+7Z1Huj3cLW7WF0nQ2gX/bZNvgpYsydkJANIRGQkK1/d4Gqs3VG5W792H+Y8Xnbim5b14ufIsLjXNbROqI8mzPamuH0fhAm6O+DTcpsYqsJriniSmzOOLVbX22o3xDfLPySkKEfqw8G1et9SmzOOPiFfrQNihtNZAJExENrAzSe9o0NXkCJzRmln+rHLB+uEhlv827BU5QGdgNBExcCS4AwoAVYhCa/7VWadcdh73t6fW34lR4wsockZMLJHXpROnC87+zwIIHr2TQRBvwHXWDjgd8Cz/QqTSnhwwf01b6HXwUJgz1gaA+x1duCqJEQuGKDVCACTW5GEwIYDPT8L9NwGlbcBdv/CWHRMHeJp+zsGbYW6erHoLnWt7Z4iEAuRkMBC5oYBGwE6oB53U6lci9sfgE2PAPNNfpQ0XX/guQRHja3m9jqbI2V8Lh11sn4O/WXmUde7Tu77LFY9KK+rRnSz+/y9sAdrtJENrAdTcZav88BXgfGo8muXrzUH/pAMbxo1wjImQuXPgmp53rc3G5T9jm8PKfz8IzpMOV+GHYFhMd41xaLGeqOQc1hqNoPx77VP0e+0f8ZcubCbWcN9gflphuRaCIETZrR5Go0sROYDri32Nngqfr7BZkz4LxFkDnNq8Z2i7hBrsMPfal/7AmN0p8lNg1iBkD0AF2IIkRv8ISEA1L/bmnT+xLNrXodta0JWuv1Iru5GhqrdIHVHtaPspO5dfGDITHHrUcKdM9WBtyNJp9HE2nAFmA2cCOQBeQCq4EcYCJwE5rciv1AvMUMphADDXcTKfUlUDc9p//xQ8J1MTVV+caemBSIS9ffoU0b/90nIdPZzjVB6dnqgSw0sQ5IAh5FkyVoAmATcBdwBBgGzAW+B2w9KwV/FBrof8BL/qB/OmKx6AJsqITaI/oSEdUHAal3Brc1614LoKXe6pWkft1k9WqmUN3DhYTpO9KERemvEUbEQUSC3rFsE1jcII+8pB3oYgNNPgo86iRkF5q0oIldaLIWTdQAEQZb5x1MJn3VyqhE6J8bMFsUBb7YOsdV/aBv71jmIwJXbJosB1zs4aPwNwK3gaAIOAJ5BEERYCixKQxDiU1hGEpsCsNQYlMYhhKbwjACt5+td6j+Hu8SlGOjChfs/uYwf3nA9crj6blJ/GjJZfQb4OVpSvTdTt2gfejW5jaeum8lB0pOdite5ohkbv5ZHhnD+nvCDLVHvB1B99DmNgt/un8lpdtcL/xsChEk9I+m6kS90/DJlw5j/r0XEJ8U3RtzlNjsCKqH/uaz/bygfeY0LKF/NL94aQGxCZGdxt+18TAfvrqZvZv1Cc4xCREs/Pkszp3W43VOlNjsCIqHllLyq1ve4ERFjdPwJz64nagY9+ehHdh1kv/8dR37t+re8cq7zuOy2yYgHCdHdoUSmx0B/9Bbisv5+/984jTsuvumMue6c3qUrsUi+eS1LbxbuAEpYfpVI7npwemYQrrVS6bEZkdAP3TRWzv511OOq0EC/OSpKxh5Xu+XRt1SXM6L//sZrS1mLrg0l4WPzOqO4NT6bMHAiUPVnQrtPg8JDWDczGx+/OTlhEeFsv6jfby+9Et665iUZwsQak438OErm1nz5k6n4fMWTeKSWz2/TMPeb49S8NAHtLaYueru87nstgnuRFPFqB0B8dBffbCHdwo3UH2qweV9KRkJPPaPG7xmx5bicgp/+QlSwr1LLuWcqV2uM6yK0UCi8lgtrz6+pkuhATy67Hqv2jJuZjZX3q2/8b7s90XUVjX2KJ2gENvivMLli/MKx/jaDk/SUNPs1n0LfjylJ10T3eaSW8YzfMJAaqsaWfb7NT2qvwW82BbnFUYAQwuK8nf42haP4qaAZvewi6O7mEyC234xm+i4CLavO0Tx2yXdTsNvBuIX5xXORl9r7QAwEmgE7igoyu/qqS4GVi3OK3wRmAyMALYBSwuK8l/xnsW+59afz8JkMu6txMTUWL7/0Ayef2wVbz7zFedMyyQp1f0X3PzNs50PPF1QlD8WeBF41Y0484EVBUX5dwJXAI0FRfnjA11o5lazy/CYhAimXGb8SuYT5+QwcU4Orc1m/vPXr7oV19/EtqWgKL/Yev4CMGFxXmGn0xAW5xUKYArgvOMpgGnrQmwpgxMMqas549rFkwmPDOXbNWWUbKhwO56/ia3N7tz2S7r61acCGwqK8oNr+zqgrdX1Iw2f2MUqR14kMSWWy639bW/+9SssZvd+fn8T2/jFeYVjreeLgLUFRflnbIGL8wpjF+cV2pcd84C3O0tscV6htjiv8B+L8wrXLc4rfLxjeENtMx++upkPX93sKfs9RmtLW6dhyenxzLnemIZBZ8y5/hyS0mI5UlbFVx/scSuOv4ntGPDbxXmF29DrYgsX5xWOXpxX+NbivMI/Ab8H7licV2hbM34u4Hw0+jvWFxTlTwXGLc4rTLMPOHWklncLN/Lxsi001LrX1WAUbS2de4s7/t8c4vr5YH8GO8IiQpm3aBIA7724qctiH/yoNWqlpqAo/yr7C4vzCucCH6HXzTYCO4FBAAVF+RPt7y0oyi/Hcf0PWy12C5CNLmhAn506fOIg9mw6wtqVu7n4+2PxF1z98WLi/WMxpokXDqV0+3GmXD6C0LCulx7zN8/mQEFR/ifABuBHwFF0b+Z6Yv3Z2Dp7zwUc/L1tKs66D3b3eqDZk7gSW3ySb72aDZNJcMN/TSdzRLJb9/uNZysoyv8ccFoRKSjK/8Z6+h/rpzt8f3Fe4d3AyoKifIftUsZMySA2IZJj5Wc4tOcUmSMGdDN579CZ2O587EIio3u/MJ8naG1uY9/GChpqmomOjyD3/MGERXQuKb8RmxdZUlCU/2lngSGhJq5cOI4zJxs4sOUoDZX6WGTSoDhSspMwt1nYu95xxfvkjASSM/rR2tzG/m8cN3UdkJVI//R4mhtbKdt8xCE8LSeJfmlxNNY1c2DrMYdw2abX2ZJSYrhgjr5mbb8B0cRGhbJr7QHSRw4gLima2tMNHN7l+HJL5phUohMiqT5Zz9G9pxzCs8amERUbwZljtRwrddyyaMiEQUREhVF5uIaTBxyXVh16Xjr7NlZQd1ofJ6073ci+jRWMmp7tcK8Nvy9GjSDOWizVVjUhLf5RlJrbzm4ghEWEMGBwgo+scU7H8duuxnPVFCP0ufy/veM/HC2r4p7fXcLYGT7Y0LYD77+0ifde+Kb9+6/+dRP9B8b50CJHSr4sb/dsALFJUTbPpqYYdYYQgsmX6t13X3/oXp+Rt7Gvs11z72S/ExpA7vmDiU2KwhRqIjYpitzzXW+/1BfqbG4x6ZJcVvx9A9vWHqSxroWoWN9WwttadLFlDE9m1gL/nD0VFhHqso7WEeXZrPRLjiHnnBTMbRZ2O6nwG43Ns02+bBhh4X66fH43UWKzY9QkvRjYud79wWVvYRsbdaezNFBQYrNj9AX65mTdmcngLWzFaGiQeDVQYjuLjBHJRMWGc/pYHVXH63xqS6tVbGFhwfMnCp4n8QAmk2DImBQASne4XqDF27Q067M+wiPDfGqHJ1Fi68CQMakAlG0/4VM7WppsYgueDgMltg7knKOLbf92xyEkI2m1eTYXY42BhhJbB7JGDUAIOLzvdPsf3Be0e7YoJbagJSomnNSsfpjbLFTsq/SZHc1NyrP1CbJH642E8p3dWyrUkzQ3tgLKswU92aOsYivxTSNBSkl9dRMAMfGdrxgZaCixOSF7lD6BsrybiyB7isa6FixmSWR0WNAMVYESm1MG5SQRFhHCqcM11Fk9jJHUWCdwxiX6x/RvT6HE5oSQUFP7vPruLvHuCY6U6TNjUzP9a7Jkbwme2qeHGXpuGvu3HuejZZupqWokLDzkuxmB1pOz3ki3nUq9ziUlICUWi/4dKbFYr0mLRIJ+lLZwfT1baZHtc+psDZVgQc3U7YSq43X85o7lNNW3GmGPA3FJUfzihWt7ux+Br1ArT9rh1kMfP3SGde/toaayoX1gXNqi2qXQ/hNKCUIgTCAQCJNACN0D2l9DgEnoRyGs95gEAv1a/7Q4Js3NJSE5IIUGSmxn0Scf2kDURml2+Gb5nz6Oao0qDEOJTWEYSmwKw1BiUxiGEpvCMJTYFIahxKYwDCU2hWH01U5dNYLgXdQqRl6nqRpevwYe7wdvLoSWrjc560uosVFP8trVsOfd776Puw2uedkrWfk5yrN5lX0f60KLiIcblkNoJGx5BSrW+9oyv0GJzVN8/Rf9OOMRGL0AJt+nf//qKd/Z5GcosXmC6grY9wGYwmBivn7tgh+DMMHO/0C9714J9CeU2DzBzn+DtMCIqyHGuidAQgYMvQQsrbD7Hd/a5ycosXmCPe/px1HXnn3d9r3kLWPt8VOU2HpLcy0cKAIE5F56dtiIq/XrpZ/o9/VxlNh6S+mnelE5eApEd9gaNTYVBk8GcwuUr/GNfX6EEltvKf9cP+Ze5jw852L9WLbKEHP8GSW23lJh3fQvc7rz8CEX6cdSJTYltt7Q2gRHNwMCBk1yfk/GVAiNghPboM63q1n6GiW23nBss15fSxkDkfHO7wmNgMwZ+nn5auNs80OU2HqDrQgdPMX1fVl5+vHgl961x89RYusN7orNVp87+IV37fFzlNh6g7tiS78ATKFwfEuf7m9TYuspNUeg+qA+yyN5lOt7w2Ng4ER9SMsm0D6IEltPOfy1fky/AExu/IwZtqK079bblNh6irtFqA1bi7QP19uU2HpKt8U2/bt4Zt/tr+BLlNh6grkNDm/Qz9MnuxcnNhWScqG1Xm8o9EGU2HrCiW3Q1qiLxzZ/zR0y+nYXiBJbT+huEWqjj9fblNh6Qk/FljVTPx4oslsbte+gxNYTeiq2/sMhNg3qT8CpXZ63y89RYusuDZVQuUd/VS91bPfiCgHZs/Vz2zy4PoQSW3epsHbmDjofQnqwy7ESm8JtKtbpx+4WoTbsxdbH6m1KbN3lkFVsGdN6Fr8P19uU2LqDue27MdHBU3uWhn29raxvTaZUYusOJ7ZDSx0k5kBcWs/TyZ6jH/d/5Bm7AgQltu5waK1+7GkRamP49/Tj/k+gtbF3aQUQSmzdwVNii0/X57e1NULZZ723K0BQYnMXKeFgsX7e0/qaPcOv0o+733V9nx8jW2qxbP4LluKHsWz+C7LF9SxktRggYNn8tMMNImUCIn0G0tyC3Po33QsdXq+vVJQxDZF2AWLgZGRLHXLHi47x06cjUiYim6qQJcscw2MyEa/NR/YfipxxHx3XzxNZlyCSRiBrK5D7HNcKETlXIhKGIKvLkKUrHcNzr0HEDUae3o088LFj+IgbENGpyFPbkYccGypi1K2IyETkiU3Iw44TPsWYO5E7XoDqMvSf0wQJ2Zgm3IfDw1hRns1dGvXdjYns55n0EnMgeaQ+ItF42jNpGk3dYb77v7VYv3eO+55NE0uAVWjyYzQhgJeAbWjyiR4b2100EQe8AVyLJhvtri8HHkOTO9xMqfvu/PVrYddbcNXf4bz8bkd3yhd/gE8fhpHXwE1veiZNA7Fs/gtUlwMWPOfZNDEFGGUV2ihgFXCdRyzuDpqsBf4J/NrOtghgaDeE1n1am/QFZABy5nou3XELQYTo67edLvVcugYhxtwJCdkQEgEJ2fp3F7i7NL0G/NV6vhh4DjjYUyPRxGxgCXAAGAk0AnegyRI3Yr8BLEETf0STx4GLgVVo4kVgMjAC2AYsRZOv9NhGe/Z9AC21egsyMdsjSQIQN1AX3LcvQfFvYd7znkvbAER4HGLCfW7f37Vn00Q/YCag1zI1+WM0+ZqL+69AE21OPrd1uPN84Gk0ORZ4EXjVLYs12QRsAK6wXpkPrECTd1qvNaLJ8R4TGsBWawV/zI0eS7Kdmb/QvdvmF+HgWs+n70e449lygaNossWtFDX5vpvpbkGT1r4EXgAK0ER/NFnpRtwyYIS17jgFuMct23rCyV36ypGmMBh7q+fT7z8Mpv8Mvngclt8IdxZB4hDP5+MHuCMKa7vW49i/YmSrUJrdjNtqvXcqsAFNWnpkQUuDPvvC0ub4kWZoa4YNzwASJtwJ8YN6lE2XzH5M78M7+AU8OwHO/yGkjdNXPzKF6AtB+zORiZDZdUe3O2LbD6SiiUhrEeYpxqOJsWhyK7AIWIsmz5x1hyZigYFocm+HuEOAt4F51qNzNKEBw4Ac4HM0+fOzwhtOwmvf69rSpFy4+Pdd39dTQiPg5pXw1m16Y+HLJd7LyxtkTIO7un75umuxafIMmigG5gAf9N6ydo4Bv0UT2cAJYKGenxgN/Ba98RAKVKOJFWhyvTU8HL3ovAt4AL3x4or1aPIWNPE+mkhDk8faQ8JiIPdyfR2O9k/I2d8TsmDSvRCV6LEHd0pkAtz0Nhwohj0r4UyZvjyqxYzfb7XV1fITVtxtjf4K+B/sxabJO7prUwdq0ORVTq6nAx+hC2ojsBOwL79uBl6zesGJZ8XUZDkQ2yE92+IaW4BsdJHrxCTDre/31H7PIwRk5+mfIMTNfja5FtiNJjpZONaDaPIT9Nbmj4CjwFxAV4RerN5M197MnjHW47nAHuu5UB+vfpwS3GOjep1tOhADrESTv/OtQX2bvrDf6BI0+amvjVCogXiFgQR3MarwK5RnUxiGEpvCMJTYFIahxKYwDCU2hWEosSkMQ4lNYRhKbArDUGJTGMb/B8qA1D6y79s/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 113.386x425.197 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run mechanical simulation\n",
    "data, _ = bls.simulate(drive, Qm)\n",
    "t, Z, ng = [data[key].values for key in ['t', 'Z', 'ng']]\n",
    "\n",
    "# Compute variables evolution over last acoustic cycle\n",
    "t_last = t[-NPC_DENSE:] * 1e6  # us\n",
    "Z_last = Z[-NPC_DENSE:]  # m\n",
    "Cm = bls.v_capacitance(Z_last) * 1e2  # uF/m2\n",
    "Vm = Qm / Cm * 1e5  # mV\n",
    "yvars = {\n",
    "    'C_m': Cm,  # uF/cm2\n",
    "    'V_m': Vm,  # mV\n",
    "    '\\\\alpha_m': pneuron.alpham(Vm) * 1e3,  # ms-1\n",
    "    '\\\\beta_m': pneuron.betam(Vm) * 1e3,  # ms-1\n",
    "    'p_\\\\infty / \\\\tau_p': pneuron.pinf(Vm) / pneuron.taup(Vm) * 1e3,  # ms-1\n",
    "    '(1-p_\\\\infty) / \\\\tau_p': (1 - pneuron.pinf(Vm)) / pneuron.taup(Vm) * 1e3  # ms-1\n",
    "}\n",
    "\n",
    "# Determine colors\n",
    "violets = plt.get_cmap('Paired').colors[8:10][::-1]\n",
    "oranges = plt.get_cmap('Paired').colors[6:8][::-1]\n",
    "colors = {\n",
    "    'C_m': ['k', 'dimgrey'],\n",
    "    'V_m': plt.get_cmap('tab20').colors[14:16],\n",
    "    '\\\\alpha_m': violets,\n",
    "    '\\\\beta_m': oranges,\n",
    "    'p_\\\\infty / \\\\tau_p': violets,\n",
    "    '(1-p_\\\\infty) / \\\\tau_p': oranges\n",
    "}\n",
    "\n",
    "# Create figure and axes\n",
    "fig, axes = plt.subplots(6, 1, figsize=cm2inch(4, 15))\n",
    "fig.suptitle('Cycle-averaging', fontsize=fs)\n",
    "for ax in axes:\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    for skey in ['bottom', 'left', 'right', 'top']:\n",
    "        ax.spines[skey].set_visible(False)\n",
    "\n",
    "# Plot variables\n",
    "for ax, ykey in zip(axes, yvars.keys()):\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    for skey in ['bottom', 'left', 'right', 'top']:\n",
    "        ax.spines[skey].set_visible(False)\n",
    "    y = yvars[ykey]\n",
    "    ax.plot(t_last, y, color=colors[ykey][0], linewidth=lw)\n",
    "    ax.plot([t_last[0], t_last[-1]], [np.mean(y)] * 2, '--', color=colors[ykey][1])\n",
    "    ax.scatter([t_last[-1]], [np.mean(y)], s=ps, color=colors[ykey][1])\n",
    "    ax.text(t_last[0] - 0.1, y[0], '$\\mathregular{{{}}}$'.format(ykey), fontsize=fs,\n",
    "            horizontalalignment='right', verticalalignment='center', color=colors[ykey][0])\n",
    "\n",
    "figs['c_cycleavg'] = fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Panel D: hybrid integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No handles with labels found to put in legend.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyEAAAEGCAYAAABy2JpUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU5d338e8v+8oWwLAZdhfQSsV612orbrjcraLVsolrrQ9WQXtX6q3W2y5qa8VWq9jq4wOUqBVEwLuCWtC61KKISxUU1BB2SEgg+0JyPX8kM04ggUkyZ+bAfN6v13nlzJnJ5Bo8Jud7ftdizjkBAAAAQLQkxLoBAAAAAOILIQQAAABAVBFCAAAAAEQVIQQAAABAVBFCAAAAAEQVIQQAAABAVPkmhJjZcWb2mpm9b2arzOzE5uOJZvZ7M/vUzD43s+tj3VYAAAAAHZcU6wZIkpllSHpZ0jXOuRfN7EJJ+ZKOlvQjScMljZSULeltM1vtnHsnZg0GAAAA0GExq4SYWegqiedI+sI592Lz4yWSLmveHyfp/znn9jrnSiU9I2ly9FoKAAAAIJL80h1ruKTtZvZ/zWyVpFf0VZVmgKRNIa/dLKl/a29iZi5087TFAAAAADokqt2xzGyspN+EPP6gefd5SedLGuOcW9ncHetFM8tTU1AKDRQmqSFKTd6Xb4LN9u3b2/09ubm5UftZHf150fxZHf15/DvG9ucdCp+Nf8fI/Dz+HSPz8zr6swAgAqytJ6JaCXHOveScO8E5d0Lz48D+RklrnXMrm48vlpQoaXDzc31D3qavmqohrb2/Oefa/LAAAAAAYs8v3bGWShoUMiPWt9VUdSiQtFjS1WaWZGbdJI2XtChmLQUAAADQKTGbHSu0YuGc225mF0l61MwyJdVKutg5V2NmsyQNkfShpBRJf3LO/SMmjQYAAADQab6YoleSnHOvSzq5leN7JU2PfosAAAAAeMEv3bEAAAAAxAlCCAAAAICoIoQAAAAAiCpCCAAAAICoIoQAAAAAiCpCCAAAAICoIoQAAAAAiCpCCAAAAICoIoQAAAAAiCpCCAAAAICoIoQAAAAAiCpCCAAAAICoIoQAAAAAiCpCCAAAAICoIoQAAAAAiCpCCAAAAICoIoQAAAAAiCpCCADfueWWW3T66aervr4+1k1BnNu+fbuuvPJKvf3227FuCgAcVpJi3QAA2NeDDz4oSXrrrbd0+umnx7YxiGs/+tGPtGTJEr333nv697//HevmAMBhg0oIAF9pbGwM7ldXV8ewJYC0efNmSdLHH38c45YAwOGFEALAV2pra1vdB2IhKysr1k0AgMMSIQSArzQ0NAT3KysrY9gSQMrOzo51EwDgsEQIAeArhBD4CSEEALxBCAHgK4QQ+ElmZmasmwAAhyVCCABf2bt3b3C/qqoqhi0BpOTk5Fg3AQAOS4QQAL5CJQR+QggBAG8QQgD4SmgIYXYsxFpS0lfLaTnnYtgSADi8EEIA+EpoCGHFdPhJaFdBAEDnEEIA+EpoCKmrq4thS4CWi2dSmQOAyCGEAPAVQgj8JLQLFucjAEQOIQSAr9AdC35CCAEAbxBCAPhKaPcXLvoQa3THAgBvEEIA+Ap3nuEnnI8A4A1CCABf4aIPfkJlDgC8QQgB4CuhIYQxIYi10POR7lgAEDmEEAC+QiUEfkIlBAC8QQgB4CuEEPgJ5yMAeIMQAsBX6I4FP2F2LADwBiEEgK9w5xl+wvkIAN4ghADwFS764CdUQgDAG4QQAL5CCIGf0D0QALxBCAHgK1z0wU9CKyGcjwAQOYQQAL5CJQR+QigGAG8QQgD4CiEEfkIlBAC8QQgB4CuEEPgJlRAA8AYhBICvhF707d27N4YtAaiEAIBXCCEAfIU7z/ATQjEAeIMQAsBXQi/6GhsbW9yJBqKNUAwA3iCEAPCV0Is+ibvPiC26YwGANwghAHxl3xDChR9iiUoIAHiDEALAVwgh8BMqIQDgDUIIAF+hOxb8hEoIAHiDEALAV6iEwE+ohACANwghAHyFEAI/oRICAN4ghADwFbpjwU9CKyGciwAQOYQQAL5CJQR+QiUEALxBCAHgK4QQ+AljQgDAG4QQAL5Cdyz4CZUQAPAGIQSAr1AJgZ9QCQEAbxBCAPgKIQR+QiUEALxBCAHgK4QQ+AmVEADwBiEEgK8wJgR+QiUEALxBCAHgK1RC4CesEwIA3iCEAPAVQgj8hEoIAHiDEALAV+iOBT9hTAgAeIMQAsBXqITAT6iEAIA3CCEAfIUQAj8hhACANwghAHyFEAI/oTsWAHiDEALAVxgTAj+hEgIA3iCEAPAVKiHwEyohAOANQggAXyGEwE9Cz0eqcgAQOYQQAL5Cdyz4CZUQAPAGIQSAr1AJgZ8wJgQAvEEIAeArhBD4CZUQAPAGIQSAr9AdC35CJQQAvEEIAeArVELgJ6GVkMbGxhaPAQAd16EQYmaZZpYY6cYAACEEfsL5CADeCCuEmFmCmU00s7+Z2U5Jn0raZmafmNn9ZjbM22YCiBdc9MFP9q18cD4CQGSEWwl5VdIQSbdJynXODXDO9ZZ0mqR/SbrPzCZ71EYAcYQxIfATzkcA8EZSmK87yzm33+0f51yJpOckPWdmyRFtGYC4RCUEfkIlBAC8EVYlpLUA0pHXAMDBEELgJ5yPAOCNg4YQMzvbzB43sxOaH1/nfbMAxCu6v8BPqIQAgDfC6Y41VdJVku4wsx6STvC2SQDiWSCEmJmcc1z0IaY4HwHAG+F0xypyzu12zv2XpHMkneRxmwDEscBFX0pKiiTuPCO2ApUQzkcAiKxwQsjfAjvOuZ9JmutdcwDEu0AISU5umuuCiz7EUuB8TE1NlcT5CACRctAQ4pxbvM/jh71rDoB4t28lhDEhiCUqcwDgjXCn6JUkmdloSbdLymv+XpPknHPHe9A2AHGIiz74SaA7FpUQAIisdoUQSfmSfirp35IaD/JaAGg3umPBT6jMAYA32htCipxzSzxpCQCIiz74C5UQAPBGe0PIXWb2hKTlkmoDB51zCyPaKgBxi+5Y8BPORwDwRntDyFWSjpaUrK+6YzlJhBAAEcFFH/yEKXoBwBvtDSFfc84d50lLAED7jwmhOxZiiVAMAN4IZ52QUP8ys2M9aQkAiIs++AtjQgDAG+2thJwq6QozK1DTmBCm6AUQUYQQ+AnnIwB4o70h5FxPWgEAzZiiF35CJQQAvNGuEOKcK/SqIQAgMUUv/IXzEQC80a4xIWY2x8y6hTzubmZPRr5ZAOIV3V/gJ1RCAMAb7R2YfrxzbnfggXOuVNKoyDYJQDwjhMBPOB8BwBvtDSEJZtY98MDMeqj940oAoE1M0Qs/YZ0QAPBGewPEA5LeNrP5alqk8DJJv454qwDELe48w084HwHAG2GFEDP7pqR/OefmmtkqSWeoaXrei51za7xsIID4wkUf/IQxIQDgjXArIVdIesTM1klaJmmBc267d80CEK9a647lnJOZxbJZiFOEYgDwRlghxDl3vSSZ2dGSzpM028y6SnpVTaHkLedcg2etBBA3Ahd9CQkJSkhIUGNjoxoaGpSUxPAzRB9jQgDAG+0amO6c+9Q596Bz7lw1dcl6U9KlklZ60TgA8Sdw0ZeYmMiChYi5QCgOdMdiogQAiIywQoiZDTWzb4Uec85VS6qQ9Hvn3GgvGgcg/gRCSEJCAiEEMUd3LADwRriVkN9LKm/leFXzcwAQEaEhJNAFi7vPiIVAAJFEIAaACAs3hAx0zn2070Hn3CpJAyPaIgBxjUoI/CIQQsyMcxEAIizcEJJ2gOfSI9EQAJAIIfCPwLlICAGAyAs3hLxrZj/c96CZXSPpvcg2CUA8ay2E0B0LsRA6UxshBAAiK9w5L6dLet7MJumr0DFaUoqkcV40DEB8am1MCBd+iAUqIQDgnXDXCdkh6RQzGyNpZPPhvznnVnjWMgBxie5Y8AsqIQDgnfau/vVPSX0l5Uk61cxOlSTn3C8i3TAA8YkQAr9orRJC10AAiIz2hpDFknZLWi2pNvLNARDvmKIXfkElBAC8094Q0r95tXQA8ASVEPhFaCWE8UkAEFnhzo4V8E8zO86TlgCACCHwD85FAPBOeyshp0q60swK1NQdyyQ559zxEW8ZgLjEFL3wi4aGBklSYmIiIQQAIqy9IeQ8T1oBAM2Yohd+EQghSUlJhBAAiLCwQoiZmWtSeLDXRK5pAOIRXWDgF4EKHJUQAIi8cMeEvGpmN5rZkaEHzSzFzM4wszmSroh88wDEGxaIg1/QHQsAvBNud6xzJV0t6WkzG6SmaXrTJCVKelnSg865D7xpIoB4whS98IvWQgjnIgBERrgrptdIelTSo2aWLKmnpGrn3G4vGwcg/nD3GX4RCBxJSUmMTwKACGvvwHQ55+olbfOgLQAQvMhLTk4mhCCmCMQA4J32rhMCAJ5q7e4zXWAQC4QQAPAOIQSAr1AJgV8wRS8AeIcQAsBXQishXPghlpiiFwC80+4xIaHMLEnSCEmpkuSceycSjQIQv1qrhNAdC7FAIAYA73QqhEh6VtI7kuolueZ9AOiwwIVfcnIyMxIhpmpqaiRJ6enphBAAiLDOhpBPnHP3RaQlAKCvLvK4+4xYC4SQtLS0FpMkOOdkZrFsGgAc8jobQurN7BVJRZLknJvY+SYBiGd1dXWSpNTUVEIIYio0hCQkJCghIUGNjY1qaGgIhhIAQMd09rdornPu7Ii0BADU9t1nINpCz0WpqYtgbW2t6uvrCSEA0Emd/S2aYWbjJZVJknPuxc43CUA8C73woxKCWDpQCElPT49l0wDgkNeuEGJmGZKGNj/8TNKrapoZq1eE2wUgThFC4BethRCJ8xEAIiGsEGJmyZLulzRFUoGa1hfpLemPzrl7zWyUc+5975oJIF5UVlZKojsWYo8QAgDeCbcS8oCkDEl5zrlySTKzLpJ+Z2azJJ0raZA3TQQQT8rKyiRJXbt25aIPMUUIAQDvhBtCzpc0zDnnAgecc2Vm9n8kFUs6z4vGAYgvzjnt2bNHEiEEsVdeXi5JyszMlEQIAYBISgjzdY2hASTAOdcgqcg596/INgtAPKqqqlJDQ4PS09OVkpKilJQUSV9N2wtEU0lJiSQpJydHkugeCAARFG4IWWNmU/Y9aGaTJa2NbJMAxKvQKojUtFaIJNXW1sasTYhfu3btkiT16NFDEpUQAIikcLtj3SBpoZldLek9SU7SSZLSJY3zqG0A4syOHTskST179pQkKiGIqUAlhBACAJEXVghxzm2RdLKZnSFphCSTtNQ5t9zLxgGIL4WFhZKkvLw8SVRCEFtbtmyRJOXm5koihABAJLVrnRDn3ApJKzxqC4A4t3HjRknSkUceKYlKCGLHOacNGzZIkgYOHCiJEAIAkRTumBAA8NyXX34p6asQQiUEsbJt2zZVVFSoW7du6t69uyRCCABEEiEEgG+sWrVKknTCCSdIohKC2Fm9erUkadSoUTIzSYQQAIgkQggAX6iqqtJ7770nSRo9erSkryohhBBE2+uvvy7pq3NRIoQAQCQRQgD4wrJly1RTU6NvfOMb+82ORXcsRJNzTosWLZIkXXDBBcHjrBMCAJFDCAEQc845/eEPf5AkXXbZZcHjVEIQCy+//LLWr1+v3Nxcfetb3woepxICAJFDCAEQcwsWLNDrr7+u7t2764c//GHwOJUQRFt1dbVuvvlmSdK0adOC1Q+JEAIAkUQIARBTb7zxhq666ipJ0q9//Wt16dIl+ByVEERTdXW1xo8fr7Vr12r48OGaPn16i+cJIQAQOYQQADFRXFysu+++W2eeeaYqKyt1xRVX6Prrr2/xGiohiIbq6mrNnz9f//Ef/6ElS5aoW7duWrhwodLS0lq8jhACAJHTrsUKAaAzGhsb9dZbb2nevHlaunRp8GLupptu0syZM4NToQYwRS+8tHbtWuXn52vBggXas2ePJGnQoEFavHixRowYsd/rCSEAEDmEEACeKy4u1l//+lfl5+eroKBAkpSQkKDvfve7mjZtms4888xWv4/FChFplZWVWrx4sfLz84NrgUhNa9NMnTpVEyZMUFZWVqvfy/kIAJFDCAHgibaqHn379tXEiRM1fvx4nXjiiQd8j8Cg4IaGBjU0NCgxMdHzduPw9OGHHyo/P1/PP/+8KioqJEldunTRJZdcokmTJmnEiBHKzc094Hukp6dLkmpqajxvLwAc7gghACJq8+bNevbZZ/XMM89o06ZNkpqqHuecc44mT56sM844I+wwYWZKTU1VbW2t6urqgheBQDjKy8u1cOFCzZs3Tx9//HHw+EknnaTJkyfrP//zP5WRkRH2+wXGiFRXV0e8rQAQbwghAFRXV6dXXnlFBQUFysrKUrdu3fbbunbtqm7duiktLW2/sRs1NTV66aWX9NRTT+mNN96Qc06S1L9/f40fP14TJkxQ3759O9S2lJQUQkic2bVrlxYtWqSVK1cqOTlZ2dnZB9yysrKUkZER3FavXq2nnnpKixcvDgaG7t276/vf/74mTZqko446qkPtCoQQKiEA0HmEECBOOef00Ucf6dlnn9Xzzz+v0tLSsL4vJSVF3bp1U/fu3dW9e3d17dpVK1eu1O7duyU19Zs///zzNX78eJ166qlKSOjcJHypqakqLy+nH/5hrq6uTitWrNBf//pXLV++PGKDv0855RRNnjxZ55133n6zXbUX3bEAIHIIIUCcKSoq0vz58/Xss8/qs88+Cx4/5phjNG7cOO3du1elpaXas2ePdu/evd9WV1ennTt3aufOnS3ed+TIkZo4caIuuugide/ePWLtZYasw9unn36q/Px8LVy4UCUlJZKauu+dccYZuuKKK5SRkaGysjKVl5ervLxcFRUVwf3QY9XV1aqsrFRVVZUSExN14YUXauLEiRo8eHDE2kolBAAihxAC+FhxcbFSU1M7fVHf2NioN998U3/5y1+0bNky7d27V5LUo0cPXXzxxbrssss0cuRI9enT56DvVV1drd27d6u0tDS4ZWRk6Nhjj+1UG9vCjET+4JzThx9+KEkaMmSIBg4cGJyytr2qqqr0v//7v5o3b57efffd4PGjjjpKl112mS6++GLl5uYedKB4W7Zv396h7zsYxoQAQOQQQgCf2bt3r/7+979rzpw5eu211yRJOTk5Gjp0qIYOHaohQ4Zo8ODBwa1Pnz5tdnkqLi7WM888o/z8fG3YsEGSlJiYqLFjx2r8+PE644wzgpWGcKWnpys9Pb1FYPHqok/6KoRw9zk2Kioq9Nxzz2nu3Llas2ZN8HhiYqLy8vJanI+DBg0K7rcWnD/99FPNmzevxbocWVlZuuSSSzRhwgQdf/zx+4038hO6YwFA5BBCAJ/YuXOn8vPzNW/ePG3dulVS0wV4UlKSdu3apV27dmnlypX7fV9qaqoGDhzYIpikpqbq5Zdf3m9q3EmTJmnChAlhVTz8IjMzU1LT3XNEz5o1azR37lwtWLBAlZWVkprC8HHHHacvv/xSmzZt0pdffqkvv/yy1e/v2rVr8HzMy8vTG2+80aLqMWrUKE2ePFkXXXRRu2aoiiW6YwFA5BBCgBhyzuntt9/WnDlz9OKLLwa7SQ0aNEhTpkzRZZddpmOOOUbbt2/X559/HtwKCgqCF4BFRUX67LPPWozvCAhMjXv55ZdrzJgxh+Q6G4EL1MCFMLxTW1urF154QXPnzm0RGE4++WRdccUVOv/885WXlxd8bUFBgb744osW52Ngf8+ePXr//ff1/vvvB98nOzs7uC7HyJEjo/75OovuWAAQOYQQoB0qKir0yiuvqFevXho2bFjwLn171dTUaOHChXr88cf16aefSmoKDOeee66uvPJKnXbaacEuVmamPn36qE+fPjrttNP2e6/y8vLghV/g6/r16zVy5EhNmDBB/fr16/gH9oHAvzEhZH/r16/XsmXL1LNnTw0bNkzDhw/XoEGD2j1Wo7i4WHPmzNHs2bNVXFwsqSkwXHrppbr88st19NFH7/c9qampOvroo1t9zjmn4uLi4DlZUFCgzMzMdq/L4Td0xwKAyCGEAGHYsGGD/vSnP+mZZ55pcQHSt29fDRs2rMU2fPhwDRkypNXpQIuKijR79mzNmTNHu3btkiT17t1bkydP1qRJkzq0lkZ2draOP/54HX/88cFjXo7RiDZCSEvOOb355pt67LHHtGLFiv2eT0xM1KBBg4LnYujXAQMGtKiGrVmzRo8//rgWLlwYnH1sxIgRuvLKKzVu3LgOh2wzU69evdSrVy+dfPLJkg6Pc5LuWAAQOYQQ4ADee+89zZo1Sy+++GJwAb5Ro0apurpaX3zxhbZu3aqtW7fqH//4R4vvMzMNGDCgRTh555139Pzzzwcv9kaOHKkf/ehH+t73vtfuweHxhBDSpL6+XkuWLNFjjz0WXP07LS1N3/3ud9WrVy+tX79e69at08aNG4Pd9pYuXdriPVJTUzVkyBANGzZMRUVF+uc//ymp6XwdO3asrrvuOn3zm9/09eDwWApUQhifBACdRwgB9tHY2KiXX35Zs2bN0jvvvCNJSk5O1iWXXKLrr79e3/nOdyRJDQ0N2rhxY/Dib/369cGtoKBAGzdu1MaNG7V8+fLge3Ox137xPjC9rKxM+fn5euKJJ4ITFvTs2VNXX321pkyZopycnBZT2dbU1OiLL77Y77xct26dtm3bpjVr1gRnucrIyNCECRN0zTXXaNCgQTH5fIeS7OxsSU3dMgEAnUMIAZrt3btXCxcu1EMPPaQvvvhCktSlSxdNmTJF11xzzX5rFgS6vQwaNEjnnHNOi+fq6+tVUFDQIpiYmSZNmsTFXjvFayWkpKREjz/+uJ588kmVlZVJkoYNG6brr79eF198cZurf6elpWnEiBEaMWLEfs+Vl5fr888/1/r167V9+3aNHTtWXbt29fRzHE66dOkiScH/HgCAjiOE4JBWV1enZcuWaevWrerfv7/y8vKUl5enI488ss2LtNbeY8GCBXrooYdUWFgoSerXr5+uu+46TZw4UVlZWe1uV3JysoYPH67hw4cHjx0OfeJj4VAKIZs2bdLs2bO1atUq9e7dO3g+5uXlaeDAgcrLy1Pv3r0PWAErKirSY489ptmzZwerP6eccoqmTp2qMWPGtLkmTDiys7M1atQojRo1ivOxAwKVkLKyMjnnDolK5tq1a/WnP/1JWVlZOvLII4O/H3v37t2pcwkAOosQgkNSXV2dnn76af3hD3/Qtm3bWn3NvheBgT++gf2MjAzNmTNHDz/8sLZs2SJJGjx4sKZNm6Zx48Z1eDVoRNahEEK2bNmihx56SE8//XRwXZa2pKWltTgPA+GkX79+ys/PV35+fnDg85gxY3TzzTfrpJNOisbHwEGkpqYqNTVVtbW1qqmpCY4R8aMPP/xQ999/f4vuoKFSU1M1YMCAFsEkdH/AgAFh38gBgI4ghOCQsnfvXs2fP18PPvigNm3aJEkaPny4LrroIu3YsUOFhYUqLCzUpk2btHPnTu3cubPFegehkpKSgutyDBs2TNOnT9eFF154SK6lcTgLdIEJrLDtJzt27NBDDz2kefPmqa6uTmamcePG6frrr1dlZaUKCwu1YcOG4HlZWFiokpISrVu3TuvWrWvzfceOHavp06frhBNOiOKnQTi6dOmioqIilZWV+TKErF27Vr/97W+1bNkySU3jfiZNmqSMjAwVFhZq48aNKiws1K5du4ITGLTliCOO2C+cBPYHDBignj17HhLVIAD+RAjBIaGhoUGLFi3SzJkzgys0Dxs2TD/96U91wQUX7De1bUNDg7Zv397i4i/0D3BhYaEqKip07LHHavr06brgggvomuBTPXr0kNQ0RsIviouL9fDDD2vu3LnBqsWFF16oW265RcOHD99v/FCo8vLyFudkIKRs3LhRffr00Q033NDqeA74Q2gIOeKII2LdnKD169frgQce0JIlS+ScU1pamq666ipNnTq11YUhKysrtWnTpuC5F9gCjzdt2qQdO3Zox44dbd7ISU1NVd++fdW/f3/169dP/fv3b7Hfr18/9enTR0lJXGoA2B+/GeBrzjm99NJLuvfee4N3jgcNGqRbbrlF48aNa7NqkZiYqH79+qlfv3465ZRTWn3fyspKlZeXcyfP53JyciQpuK5KLFVUVGjWrFmaNWtWcNXs8847Tz/96U91zDHHhPUe2dnZGjlyZKsXhozT8D+/DU7fvHmz7r//fi1YsECNjY1KSUnR5ZdfrhtvvPGAISkzM7PNxSalr27k7BtOAvubNm1SaWlpcDHKtiQkJCg3N7dFMGktrPixqgTAW4QQ+Nb777+vu+++WytXrpQk9e/fX7fccosuvfTSTt9ZMzNlZWUx1eYhwA8hpL6+Xk899ZR+97vfBVcUP/vss/Vf//VfLRaJxOGve/fukmIfisvKyvTwww/r8ccfV21trZKSkjRp0iRNmzZN/fr16/T7h97I+eY3v9nqa6qqqrRlyxZt3rw5+HXf/R07dgTXU2qroiI1VTz79u2rPn36KDc3V3369Gl168hEIQD8iRCCiKutrdVLL72kpKSk4KrJgS07O/uglYeNGzfq3nvv1aJFiyQ1/dH/yU9+ossvv5xF/eJQZ0PIJ598olmzZmnHjh3q1auXevfu3ebXnJycFgHXOaelS5fq17/+dXDa5hNPPFF33nlncCVwxJdAV7vOVK3effddbdq0ST169Gjx+zEnJ+egv+Pq6uo0d+5czZw5U6WlpZKkiy66SD/72c+Ul5fX4TZ1REZGRnAx1rbU19dr27ZtBwwqW7duVUlJiUpKSoILcbYlKytrv2DSWmjp0aMHVW7A5wghiJjGxkY999xzuvfee9ucsSolJWW/YBK6vffee5o9e7bq6uqUlpama6+9VjfeeGOwCwTiT2gICaxaH46ioiLdd999evrpp8P+PjMLXhj27t1bFRUVWr16taSmboD//d//rQsuuICLmzjWp08fSR0LIQUFBbr77rv10ksvtfmarl27tvn7MSUlRTNnzgx2fzr55JN11113adSoUR37MFGQnJwcHNDelsbGRhUVFWnr1q3avn27tm3b1uZWUSxu/+oAAA3eSURBVFERXHvpQFJSUpSbm6vc3Fz17t1bvXv31hFHHBHcD328780HANHB/3WIiA8++EB33nmnVq1aJUk65phjdNJJJ6m4uFhFRUXBrbKyUlu2bAlOidsaM9P3v/99zZgxQ/3794/WR4BPpaenq0ePHiopKdHOnTsP+vq6ujo9+eSTmjlzpsrLy5WcnKwpU6ZowoQJKikpUVFRkXbu3Nnq1127dgW3Tz/9VFJTN5FAJY5pmxGohLR1o6U1FRUV+v3vf68///nPqq+vV2Zmpi655BJVV1e3+P1YXFysPXv2aM+ePQectWrIkCG64447NHbs2MMiECckJOiII4446EB/55x2796tbdu2HTSslJWVBcewHIyZKScnp82Qsu+xzMzMw+LfHYg1Qgg6pbi4WPfee2/wbnOvXr10++2369JLL91vxipJ+/3R3Xerrq7WlVdeqeOOOy4GnwZ+NXDgQJWUlGjDhg0H7HKyfPly3XXXXcGuU2eddZb+53/+R0OGDDngjFUBe/fuDYadnTt3qrKyUkcffXRwkTogUAnZvHnzQV/b2Nio+fPn65577gkG6B/84Ae67bbb9LWvfa3V15eWlrb5+7GkpETHHXecJkyYEJeB2MzUvXt3de/eXccee+wBX1tVVaXt27drx44d2rlzZ/Br6BY4tmvXLhUXF6u4uFiffPLJQduRnp6unj17BrecnJw2Hwf2GXgP7I8Qgg6pr6/Xk08+qQceeCB4t/naa6/VzTfffMALtvT09AOW5pkdCK0ZOHCgVq9erYKCglZDSEFBge644w6tWLFCUtOd4rvvvltnnnlmu35OUlJS8I5nAOckQh111FGSFKyUteWDDz7Qbbfdpg8++ECS9PWvf12/+tWvDth1KiEhQTk5OcrJyWlz1irOx/BkZGRo8ODBGjx48EFfu3fvXu3atWu/oNJWcKmurtamTZuCa1WF2562Asq+jwNBK5wxlMChjBCCdnv33Xd16623Bv8IjxkzRr/4xS80dOjQGLcMh6vABdlHH32k008/PXi8trZWjzzyiB566CHV1tYqOztbt9xyi66++momMYAnAufiunXrgoudhiorK9N9992n2bNnyzmn3Nxc3X777br44otZi8inkpKSwuoOJn01vXugchJaRQnd3/dxVVVV2N3DAhITE4OBpEePHsH9fR+3tp+enk6Age8RQhC23bt365577tFf/vIXSVJeXp5++ctf6qyzzuKXHTz1jW98Q5L09ttv66abbpIkvfXWW5oxY0aw69Wll16qO++8U7169YpZO3H4y8rK0rBhw7R+/XqtWrVKAwcOlNR0cfrCCy/o5z//uXbs2KHExERdd911+slPfqLMzMzYNhoRE5jePSsrK/jf/mCcc6qoqDhgSAl9XFpaqpKSElVVVQWfa6+UlJRWA0q3bt3UtWtXde3atcX+vhshBtFACMFBOee0ePFi/fznP1dRUZGSkpI0depUTZ8+nX6uiIpTTz1VycnJev311/Xxxx/rz3/+s+bPny+pqevVb37zG33rW9+KcSsRL84++2ytX79ezz77rG699VYVFhbqtttu06uvvipJGj16tH7zm98cdNwC4oOZKTs7W9nZ2WEHF6lpko3S0tJgKGltv63namtrgyved0RSUlKbAeVAWyDYZGdnM4AfB0UIiQPOOb366qt65ZVX5JxTamqqUlJSWv2677Hk5GQ98sgjwT+uJ510ku6///5gv2ggGnJycvS9731Pzz33nM4++2xJUmpqqqZNm6apU6cqNTU1xi1EPLnmmmv06KOP6pFHHtHatWv12muvqaamRl27dtXtt9+uSZMm0fUKnZaSkhJ2N7F9VVdXtxpQdu/eHZyBrbUt8HxtbW1wpsCOMjNlZmYGA1hrW1ZWVtjPZWRkEGoOM4SQw1xRUZFmzJihpUuXdup9unbtqjvuuEMTJ07kjyti4r777tPLL7+s8vJyffvb39Z9992nQYMGxbpZiENf//rXde211+qJJ57QsmXLJEmXXHKJ7rrrLroDwhfS09OVnp7e6iyV4aitrT1gWDnYVl5erqqqKlVUVKiioqJdU1q3JSEhoUUwCexnZmYecMvIyDjoa1gnJjb4Vz9MOee0aNEi3X777SotLVVWVpZuvfVW9e/fX7W1taqrq2vxtbVjga89e/bUtGnT+OOKmBo6dKhWrVqlt99+W+eccw53xBBTs2bN0tChQ/X888/rxhtvbPdMbICfpaam7jdTYHs1NDSooqJC5eXlwW3fx61tbb2mpqZGZWVlKisri+AnbZKcnNyuMBMIeenp6crIyGjx+EBbYmJixNt+KCOEHIaKior0s5/9TC+++KIk6dvf/rYeeOABjR49ukPvx5SQ8Ivhw4erS5cusW4GoKSkJM2YMUNXXHFFrJsC+FJiYmJwrEgk1NfXtwgogf3KyspWt6qqqjaf2/c19fX12r17t3bv3h2RtrYlOTn5gCElnECTlpamtLQ0paamBr8ebD8lJcWXvVgIIYcR55yWLFmi2267TaWlpcrMzNRdd92lyZMnc9cYAAAcspKTk4MzfUWSc061tbXtCjPV1dVhbVVVVS0e19fXq76+3pNqzsGEjv0NN7yEE24CY4jb2j/Q4tOEkMNEaWmpZsyYoRdeeEGSdNppp+mBBx7QgAEDYtwyAAAAfzKzYHUhJyfHs58TCDvhBpi2Qk1NTU2wG/3B9gNf6+rqglt5eblnn7Gtz90WQshh4M0339RNN92kbdu2KTMzU3feeaemTJlC9QMAAMAHQsNOpKs5B9PY2NhiDHC44eVg+4Fgc6BxxQdCCDmE1dbW6re//a1mzZol55xGjx6thx9+uF3zkAMAAODwlZCQEAxAfkIIOUStW7dON9xwgz7++GMlJibq5ptv1rRp05hmDgAAAL7HFWs7LF++XM654CYpJo83btyoX/7yl6qpqVFeXp7++Mc/dnjmKwAAACDaCCHtcNZZZ8W6CS384Ac/0K9+9StlZWXFuikAAABA2Agh7TBmzBiZWYtN0gEfe/GaxMREnX766TrvvPOi+vkBAACASCCEtMOKFSti3YQgFhAEAADAocp/yycCAAAAOKwRQgAAAABEFSEEAAAAQFQRQgAAAABEFSEEAAAAQFQRQgAAAABEFSEEAAAAQFQRQgAAAABEFSEEAAAAQFQRQgAAAABEFSEEAAAAQFQRQgAAAABEFSEEAAAAQFQRQgAAAABEFSEEAAAAQFQRQgAAAABEFSEEAAAAQFQRQgAAAABEFSEEAAAAQFQRQgAAAABEFSEEAAAAQFQRQgAAAABEFSEEAAAAQFT5JoSY2Y1m9pmZfWBmT5tZj+bj6Wb2pJl9bGafNO+nx7q9AAAAADrGFyHEzMZImiHpTOfcCZJelPTn5qdvl5Qk6fjmLV3SbbFoJwAAAIDOS4p1A5qdKOnvzrnNzY8XSnrCzFIkvS5pg3OuUZLM7H1JI2LTTAAAAACd5YtKiKSVks4ws7zmx1dJSpGU45x72Tm3TpKan58uaX5rb2JmzsxcNBoMAAAAoGPMOX9cs5vZNZJukNQo6UlJv5Q03Dm3q/n5EyU9L+kx59w9bbxHiw/jnDNPGw0AAACg3WIWQszsF5K+1/zwVUmPOOc+b36un6SPJPV0zjkzGy/pUUk/ds49FZMGAwAAAIgIX1RCzOwoScslHeucKzOzRyU1Oud+bGbflfR/JZ3vnFsV04YCAAAA6DRfhBBJMrMfq6k7VoKkN9VU9ag2s88k9ZC0JeTlbznnbohBMwEAAAB0km9CCAAAAID44JfZsQAAAADECb+sE4J2YBpiAAAAHAramq2WSggAAACAqGJMCAAAAICoohICAAAAIKoIIQAAAACiihACAIgJM5tsZh+a2Qdm9k8zGx3yXHHz8cA2qfn4MDN73czWmNk7ZnZ07D4BAKCjmB0LABB1ZnaUpPslfd05t83Mzpe0UNKRzc+VOOdOaOVb8yX93jn3lJmdJ2mBmR3nGOAIAIcUKiEAgFiolXStc25b8+NVknLNLEXSKZIazOwNM/vIzH5uZolm1k/S0ZKekSTn3FJJWZJG7fvmZlZjZvc0V0vWmNllZjbfzD41sxVmltn8urubf8YqM3vJzPpE4bMDQNyjEgIAiDrn3AZJGyTJzEzSTElLnHN1ZpYk6e+SfiYpWdLfJJVJ+pekrc65xpC32iypv6TV+/yIVEnbnXPfMLMZkp6QdIykbZLelXShmb0habqk3s65WjP7iaSTJS2K/CcGAIQihAAAYqa5IjFb0gBJ50qSc+7xfV4zU9JNkt6RtG+3K5PU0MbbP9f89QtJ/3bObWl+vwJJPSRtkfShpNVmtlTSUufc8k5+JABAGOiOBQCICTM7UtI/1RQixjjndjcfv9zMjg99qaR6SRsl9WmunAT0VVM1pDW1Ifv1+z7ZXFH5jqQrJe2S9KCZ/bZjnwYA0B6EEABA1JlZtqTXJC10zo13zlWHPD1S0i+ax4GkS/qxpL865zZL+lzSD5rfY6ykRkn/7mAbvibpY0lrnXP3SnpQ0kkd/EgAgHagOxYAIBZ+LClP0jgzGxdy/ExJd0v6o5rCRbKk+Woa0yFJEyQ9bmZ3SKqRdOk+Y0TC5pz70MyelbTKzCokVaup2xcAwGPGrIYAAAAAoonuWAAAAACiihACAAAAIKoIIQAAAACiihACAAAAIKoIIQAAAACiihACAAAAIKoIIQAAAACiihACAAAAIKr+P1EBECtQ1/5pAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 792x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run effective simulation\n",
    "pp = PulsedProtocol(150e-3, 100e-3, 100., 0.5)\n",
    "data, meta = nbls.simulate(drive, pp, method='sonic')\n",
    "\n",
    "# Plot charge density response\n",
    "fig = CompTimeSeries([(data, meta)], 'Qm').render()\n",
    "ax = fig.axes[0]\n",
    "line = ax.lines[0]\n",
    "line.set_color('k')\n",
    "line.set_label(None)\n",
    "ax.legend(fontsize=fs, frameon=False)\n",
    "t = data['t'].values\n",
    "ax.set_title('')\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "for key in ['top', 'right']:\n",
    "    ax.spines[key].set_visible(False)\n",
    "for key in ['bottom', 'left']:\n",
    "    ax.spines[key].set_position(('axes', -0.03))\n",
    "    ax.spines[key].set_linewidth(2)\n",
    "ax.spines['bottom'].set_bounds(low=t.min() * 1e3, high=t.max() * 1e3)\n",
    "ax.yaxis.set_tick_params(width=2)\n",
    "ax.yaxis.set_major_formatter(FormatStrFormatter('%+.0f'))\n",
    "ax.set_xticks([])\n",
    "ax.set_xlabel('{}s'.format(si_format((t.ptp()), space=' ')), fontsize=fs)\n",
    "ax.set_yticks(ax.get_ylim())\n",
    "for item in ax.get_yticklabels():\n",
    "    item.set_fontsize(fs)\n",
    "    \n",
    "figs['d'] = fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Save figure panels\n",
    "\n",
    "Save figure panels as **pdf** in the *figs* sub-folder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "saveFigsAsPDF(figs, figindex)"
   ]
  }
 ],
 "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"
  },
  "papermill": {
   "duration": 0.044879,
   "end_time": "2020-04-28T19:48:13.775289",
   "environment_variables": {},
   "exception": null,
   "input_path": "figure_02.ipynb",
   "output_path": "figure_02.ipynb",
   "parameters": {},
   "start_time": "2020-04-28T19:48:13.730410",
   "version": "2.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}