{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d2652dea-87f8-416b-b430-85e3c0ebc1c0",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Part 3. Compare simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "950ceee7-079c-4707-bf19-6453009311cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4a2bf064-24f9-4fd6-8b47-8cfe94d93af3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rc(\"figure\", dpi=80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "caf39f01-cb2a-4dae-b3d7-7d3798aaf376",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_values(sim):\n",
    "    N = sim.params.N\n",
    "    Uc = sim.params.forcing.milestone.movement.periodic_uniform.speed\n",
    "    Dc = sim.params.forcing.milestone.objects.diameter\n",
    "    Lf = sim.params.forcing.milestone.movement.periodic_uniform.length\n",
    "\n",
    "    Fhc = Uc / (N * Dc)\n",
    "    Rec = Uc * Dc / sim.params.nu_2\n",
    "\n",
    "    period = sim.forcing.get_info()[\"period\"]\n",
    "    t_statio = period\n",
    "    averages = sim.output.spatial_means.get_dimless_numbers_averaged(\n",
    "        tmin=t_statio\n",
    "    )\n",
    "\n",
    "    U2 = averages[\"dimensional\"][\"Uh2\"]\n",
    "    epsK = averages[\"dimensional\"][\"epsK\"]\n",
    "    Gamma = averages[\"Gamma\"]\n",
    "    Fh = averages[\"Fh\"]\n",
    "    R2 = averages[\"R2\"]\n",
    "    R4 = averages[\"R4\"]\n",
    "\n",
    "    return N, Uc, Dc, Lf, Fhc, Rec, U2, epsK, Gamma, Fh, R2, R4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "54d1ff31-08ab-4159-b1f9-8b0999437756",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['ns3d.strat_144x144x48_V4.5x4.5x1.5_N0.2_Lf3.5_U0.01_D0.5_2021-10-05_15-35-47',\n",
       " 'ns3d.strat_144x144x48_V4.5x4.5x1.5_N0.2_Lf3.5_U0.02_D0.5_2021-10-05_16-12-36',\n",
       " 'ns3d.strat_144x144x48_V4.5x4.5x1.5_N0.2_Lf3.5_U0.04_D0.5_2021-10-05_16-49-31',\n",
       " 'ns3d.strat_144x144x48_V4.5x4.5x1.5_N0.2_Lf3.5_U0.06_D0.5_2021-10-05_17-25-48',\n",
       " 'ns3d.strat_144x144x48_V4.5x4.5x1.5_N0.2_Lf3.5_U0.08_D0.5_2021-10-05_18-01-23']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import fluidsim as fls\n",
    "\n",
    "path_dir_data = Path(fls.FLUIDSIM_PATH) / \"tutorial_parametric_study\"\n",
    "path_runs = sorted(path_dir_data.glob(\"*\"), key=lambda p: p.name)\n",
    "[p.name for p in path_runs]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52889d94-4fa8-4591-9943-ca180f6a8477",
   "metadata": {},
   "source": [
    "## Compare simulations with Pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8885d07e-0620-4fbd-a50d-42189bb1b88f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>N</th>\n",
       "      <th>Uc</th>\n",
       "      <th>Dc</th>\n",
       "      <th>Lf</th>\n",
       "      <th>Fhc</th>\n",
       "      <th>Rec</th>\n",
       "      <th>U2</th>\n",
       "      <th>epsK</th>\n",
       "      <th>Gamma</th>\n",
       "      <th>Fh</th>\n",
       "      <th>R2</th>\n",
       "      <th>R4</th>\n",
       "      <th>min_R</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.2</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.5</td>\n",
       "      <td>3.5</td>\n",
       "      <td>0.1</td>\n",
       "      <td>5000.0</td>\n",
       "      <td>0.000030</td>\n",
       "      <td>5.564985e-08</td>\n",
       "      <td>0.300576</td>\n",
       "      <td>0.009451</td>\n",
       "      <td>1.391246</td>\n",
       "      <td>0.170957</td>\n",
       "      <td>0.170957</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.2</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.5</td>\n",
       "      <td>3.5</td>\n",
       "      <td>0.2</td>\n",
       "      <td>10000.0</td>\n",
       "      <td>0.000109</td>\n",
       "      <td>3.917927e-07</td>\n",
       "      <td>0.468731</td>\n",
       "      <td>0.018019</td>\n",
       "      <td>9.794817</td>\n",
       "      <td>2.204561</td>\n",
       "      <td>2.204561</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.2</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.5</td>\n",
       "      <td>3.5</td>\n",
       "      <td>0.4</td>\n",
       "      <td>20000.0</td>\n",
       "      <td>0.000419</td>\n",
       "      <td>2.834370e-06</td>\n",
       "      <td>0.618684</td>\n",
       "      <td>0.034420</td>\n",
       "      <td>70.859251</td>\n",
       "      <td>30.298269</td>\n",
       "      <td>30.298269</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.2</td>\n",
       "      <td>0.06</td>\n",
       "      <td>0.5</td>\n",
       "      <td>3.5</td>\n",
       "      <td>0.6</td>\n",
       "      <td>30000.0</td>\n",
       "      <td>0.000914</td>\n",
       "      <td>8.980284e-06</td>\n",
       "      <td>0.629455</td>\n",
       "      <td>0.049961</td>\n",
       "      <td>224.507095</td>\n",
       "      <td>139.360780</td>\n",
       "      <td>139.360780</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.2</td>\n",
       "      <td>0.08</td>\n",
       "      <td>0.5</td>\n",
       "      <td>3.5</td>\n",
       "      <td>0.8</td>\n",
       "      <td>40000.0</td>\n",
       "      <td>0.001652</td>\n",
       "      <td>2.199885e-05</td>\n",
       "      <td>0.561550</td>\n",
       "      <td>0.067149</td>\n",
       "      <td>549.971217</td>\n",
       "      <td>463.932956</td>\n",
       "      <td>463.932956</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     N    Uc   Dc   Lf  Fhc      Rec        U2          epsK     Gamma  \\\n",
       "0  0.2  0.01  0.5  3.5  0.1   5000.0  0.000030  5.564985e-08  0.300576   \n",
       "1  0.2  0.02  0.5  3.5  0.2  10000.0  0.000109  3.917927e-07  0.468731   \n",
       "2  0.2  0.04  0.5  3.5  0.4  20000.0  0.000419  2.834370e-06  0.618684   \n",
       "3  0.2  0.06  0.5  3.5  0.6  30000.0  0.000914  8.980284e-06  0.629455   \n",
       "4  0.2  0.08  0.5  3.5  0.8  40000.0  0.001652  2.199885e-05  0.561550   \n",
       "\n",
       "         Fh          R2          R4       min_R  \n",
       "0  0.009451    1.391246    0.170957    0.170957  \n",
       "1  0.018019    9.794817    2.204561    2.204561  \n",
       "2  0.034420   70.859251   30.298269   30.298269  \n",
       "3  0.049961  224.507095  139.360780  139.360780  \n",
       "4  0.067149  549.971217  463.932956  463.932956  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from pandas import DataFrame\n",
    "\n",
    "# fmt: off\n",
    "columns = [\"N\", \"Uc\", \"Dc\", \"Lf\", \"Fhc\", \"Rec\", \"U2\", \"epsK\", \"Gamma\", \"Fh\", \"R2\", \"R4\"]\n",
    "# fmt: on\n",
    "\n",
    "values = []\n",
    "for path in path_runs:\n",
    "    sim = fls.load(path, hide_stdout=True)\n",
    "    values.append(get_values(sim))\n",
    "\n",
    "df = DataFrame(values, columns=columns)\n",
    "df[\"min_R\"] = np.array([df.R2, df.R4]).min(axis=0)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9bb6e45d-1db1-4fbe-bc95-39207fdba39b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAEhCAYAAAAedPXMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAxOAAAMTgF/d4wjAAAizUlEQVR4nO3dfZwdVZ3n8c+3050Q7ACCCWTshAAJEZERWNDBgZXxYWdVJiLBGVFcZGEIMvpyBXVdZwYdEF18aZxVZEgCY1R2RDSIGR9HBceB+BA2CSJgQpJXCB3CQ3BI6CEJ6e7f/lHVUFz6oW73faq63zev8+JW1alzf/emk1+fU6dOKSIwMzMrmo5mB2BmZjYeTmBmZlZITmBmZlZITmBmZlZITmBmZlZITmBmZlZInc0OoJ6mTJkS06dPb3YYZmYj2rZt2zMRMWWi7Zx40pR49JGB3PUf3jZwV0ScPNH3baZSJ7Dp06fT29vb7DDMzEYk6fFatPPo9gHu23RE7voHTd04sxbv20ylTmBmZu1DaFDNDqKhnMDMzMoinMDMzKyA3AMzM7NC0mCzI2gsJzAzszIIwAnMzMyKSG32cBEnMLM2Frvugwc+D30PwIHHwbzL0IvmNDssGycPIZpZW4jdD8Pq/waDzwABj/8c/n0NceoPUNcBzQ7PxmOwvbpgXkrKrF09fCvEIMnFE4DBJJk9+qMmBmUTochfysA9MLN2tW8XROXSQwH9fU0Jx2qgzYYQ3QMza1czXg+a9Px90Q/TX9uceGxiAtQfuUsZOIGZtSkdfDIc9V5QZ1I6psDLLkcvOrLZodk4CA8hmlkb0RF/SfS8A/Y8DPvPQpP2b3ZINhFtNoToBGalFXv3sPsXP2Tf5vvpnHk4U099Mx3TDmp2WC1HXdOga36zw7Aa8DR6sxKIgX6evPZvGNixHfr3sW/jPexZfRsvvmwxHS/yFHEroYiktBFfA7NSeuZ3axjY8Qj070t2DPQTe3ez59c/bW5gZnWkwfylDNwDs1IafOIxnru/KdW/L+mRmZVVSRJTXu6BWSl1HnkMDFb8be7souuoVzQnILMGqPUsREn7SbpV0gZJd0v6saS5w9SbI2lA0rpMOarWn6+SE5iVUlfPUez3mv8KHR0weQpMmkTXUa9gyiv/uNmhmdXPYBUlv6XA/Ih4JfAd4PoR6j0VEcdnyqZxfYYqeAjRSqv7z85jv1e9nv6tDzBpxkvpnD0Pqb0e+GdtpA6PU4mIPcD3M7t+CXyotu8yfk5gVmqdh/bQeWhPs8MwawChqPsvaB8g6YUN50WSVgOTgFuBqyJesFZZTXkI0cysLKobQuyW1Jspl47WtKSPAXOB/zXM4e3ASyPiZOANwGnAZbX4SKNxD8zMrCyqG0Lsi4hcwxOSPgScBbwhIp6uPB4Re4HH0te/l/SPwDuBz1QVUZXcAzMzK4uoouSU9szOAd4YEU+OUGeGpK709RSSZLd2nJ8iN/fAzMxKQoO1vQYmqQf4HLAZuD2dBLU3Il4t6Qrg4Yi4DjgVuELSAEleuQ24qqbBDMMJzMysLGq8klRE9JIsdD/cscszr28Bbqntu4+t7kOIkuZJWpXeCLda0rEj1DtO0s8k3Z+WszLHLpD0gKRNkpYNdVXNzCwVQL/ylxJoxDWwJcDSiDgauBpYXllB0v4kUzP/JiKOAV4B/Ft67AjgSpJZLXOBQ4GLGhC3mVmx1OEaWCurawKTNAM4Cbgx3bUCmDXMUiTvBH4ZEXcARMRARDyeHjsbWBkRj0REANeRXFA0M7OsQeUvJVDvHtgsYHtE9AOkCWgrMLui3suBvZK+m66h9VVJ09Njs4EHM3W3DHM+kMyWyd7X0NfXV8vPYmbW2kL5Swm0yjT6TpKb3xYBJwDbgH+otpGIWBwRPUOlu7u7xmGambWo8ONUau0hYKakzojoVzIHczZJLyxrK3B7RGwDkHQj8KPMseyqxnOGOd/MzEoyNJhXXXtgEfEYsAY4N921EOiNiI0VVW8GTpY09KjcNwN3p69XAAskHZYmwIuBm+oZt5lZIbXZJI5G3Ae2CFierqO1CzgfQNL1JJMzVkbEVkmfAlZJGiQZQrwIICI2S/o4cGfa3s9IZjaamVlWm/XA6p7AImI9cMow+y+s2P4a8LUR2lgGLKtLgGZmZVGSyRl5eSUOM7OyKMnkjLycwMzMyiBwD8zMzIpIhK+BmZlZIbkHZmZmheRrYGZmVkjugZmZWSH5GpiZmRWOZyGamVlhDTiBmZlZAXkavZmZFZOHEM3MrJDcAzMzs8LxJA4zMysmuQdmZmbFFCV5UGVeTmBmZmXhIUQzM8trMLaxb/CfCfrp6ngzk3RkE4NxAjMzsxz6B3/Ff/SfBwgI9g58gf07v0hXx582JZ5osx5YR7MDMDMrqt39lwPPAHvT//ezu/9viWjSsvCDyl9KwD0wM7NxGmTTC/YFOwh2IQ5qbDBtOI3ePTAzs3Hq4IgX7BMHIw5oQjTJEGLeUgZOYGZm4zS18++ALmBy+v9OpnZegdSkf1oHqygl4CFEM7Nx6ux4DdO6fsIzgysh9tE16QwmaW7zAipJzyovJzAzswno0Cz2m/RXzQ4DaL/V6Ovez5U0T9IqSRskrZZ07DB1Tpe0W9K6TJk61jEzMxuipAeWt5RAI3pgS4ClEbFc0tnAcuDkYeqtj4jjR2hjtGNmZm0vAmKgvaY11PXTSpoBnATcmO5aAcySmjlIbGZWUm3WA6t3up4FbI+IfoCICGArMHuYukdJWpMOM15SxbFnSbpUUu9Q6evrq9kHMTNrdRH5Sxm0yiSONUBPROyU1AN8X9KOiLh5jGPPExGLgcVD2z09PSX5YzIzy8GTOGrqIWCmpE4ASSLpfW3NVoqIXRGxM33dC3wdOG2sY2Zm9pxa38gsaT9Jt6aT8O6W9OORLgFJOkPS7yQ9IOkWSXW/m7uuCSwiHiPpQZ2b7loI9EbExmw9STOV3vknaRpwBrB2rGNmZpZRn2tgS4H5EfFK4DvA9ZUVJHUDNwBnRsQ84GHgb2vwiUbViCkri4BFkjYAHwXOB5B0vaQFaZ2FwD2S7gZ+CfwY+HKOY2ZmBhDJfWB5S64mI/ZExPfT+QuQ/Bs8Z5iqbwLWRsTv0u1rgXMm+pHGUvdrYBGxHjhlmP0XZl5fA1wzwvkjHjMzs4zqelbdknoz24vTeQSj+QBJL6zSbODBzPYW0stHQ5P46qFVJnGYmdmEVL1Ib19E9ORuXfoYMBd4fbWR1Ut73fVmZlZmdXoemKQPAWcBb4qIp4epshU4PLM9h8wtVPXiBGZmVhL1uA9M0qUk17PeGBFPjlDth8CJkl6Wbl8C3DSBj5KLhxDNzEqi1s/5Su+9/RywGbg9uROKvRHxaklXAA9HxHUR8ZSkC4Fb09umfgucV9NghuEEZmZWFjW+kTm993bYRiPi8ortlcDKmgYwBicwM7MyiNr3wFqdE5iZWVk4gZmZWRG5B2ZmZsXUZov5OoGZmZVAAINOYGZmVjzleVBlXk5gZmZlkC7m206cwMzMSsKTOMzMrJicwMzMrIjcAzMzs0KKwWZH0FhOYGZmZeEemJmZFZGHEM3MrJCcwMzMrHgCDyGamVkRyTcym5lZ8QQeQjQzs6KKZgfQWE5gZmYl0W49sI56v4GkeZJWSdogabWkY4epc7qk3ZLWZcrUzPELJD0gaZOkZZK66h23mVnRRCh3KYO6JzBgCbA0Io4GrgaWj1BvfUQcnym7ASQdAVwJnAbMBQ4FLqp/2GZmxRIDHblLGdT1U0iaAZwE3JjuWgHMkjS3imbOBlZGxCMREcB1wDm1jdTMrODCPbBamwVsj4h+gDQBbQVmD1P3KElr0mHGSzL7ZwMPZra3jHC+mVlba7cE1iqTONYAPRGxU1IP8H1JOyLi5moakXQpcOnQ9oEHHljjMM3MWle02SzEevfAHgJmSuoEkCSS3tPWbKWI2BURO9PXvcDXSa55kdY9PFN9TuX5mXYWR0TPUOnu7q7lZzEza2lF6oFJOkLSf85O2Ev3vzFvG3VNYBHxGEnv6tx010KgNyI2ZutJmimpI309DTgDWJseXgEskHRYmgAvBm6qZ9xmZsUjGKyiNDNS6V3Ar4EvARsknZI5fHXedhoxFWURsEjSBuCjwPkAkq6XtCCtsxC4R9LdwC+BHwNfBoiIzcDHgTuBjcDjJDMbzcwso0A9sA8DJ0TEcSQ54RuSXpceyx1c1dfAJE0BpgxtR8Su0epHxHrglGH2X5h5fQ1wzShtLAOWVRurmVm7KNhSUkovFxERP5H0FmClpIuoYj2R3AlM0h+R9IqOrjg0KW8bZmZWPwVKYB2SpkXEUwARcU+axL4HHJC7kSre8P8A7wF+AxwEXA58pIrzzcysXop1H9hSknuEnxUR9wFvBv5f3kaqGULsiohfSepMs+ZVklYDn6uiDTMzq5fmJ6ZcIuKLI+y/H/gveduppge2L/3/E5JOlDQdmF7F+WZmVkcxmL+0IklvlbR27JqJahLYTZIOAT4F/CvJPV4jTrwwM7PGKsoQoqRXS7pX0lOSlqb3hP0C+CLw+bzt5B5CjIihRv9F0sHAfkMX4MzMrPmanZiq8PfAtcBPgbcDvwK+DbxuaCH3PKqaRi/pGGDe0HmSiIhbqmnDzMzqo0AJbP+I+FL6+sp0/dtLImKgmkaqmUb/OeBdwH3A0JsE4ARmZtYCCpTAKhPV9mqTF1TXA3srcGREPF3tm5iZWb2pMLMQgZdJWjPSdkScmKeRahLYQ8CeKuqbmVmDRMDgYG1XB5T0BWAByYLqJ0TEumHqnA78AFif2X3KGNey3lSL+KpJYB8BvinpR2QSWUR8tRaBmJnZxNRhevy3gM8Ad4xRb31EHJ+30Yj41zz1JF0dEf9zpOPVJLCLgT8kWWgxew3MCczMrAXU+hpYRPwckgl7TfJGoCYJ7HRg/tDTla2xIoLe2x9m+6pH6J7VzdyFRzK5u6vZYZlZC2niJI6j0mtYA8CXI+LaGrU76geqJoFtHqsxq5+ff/BOHvjWZmIg6OjqYO3n7+asnyxg6kv2a3ZoZtYiqkxg3ZJ6M9uLI2LxON52DdATETsl9QDfl7QjIm4eR1uVRl2ZvtoE9jNJ3+H518C+MM7ALKff3//vbPjGJmIg+bMcGBhg94493LP0Xl71sf/U5OjMrFVUmcD6IqJn4u/53CO1IqJX0teB04BaJLBRVZPAJgMbgGMy+3I/t8XG74l7f8+kKZPof/q50dvBZwZ57K7HmxiVmbWaZgwhSpoJPBoRg5KmAWcAN9Sq+dEOVrOU1PkTj8XG45BXHMzA3uff49cxuYNDT5rRpIjMrOVE7ROYpCXAW4DDgB9Jeioi5kq6HlgZESuBhcB7JfWT5JRvkjw7Mk/7RwKz082tEbG5osqoK9NXu5TUm8gsJQUwzjFTq8LBL3sx88+Zy4abNz17DWy/Q/bjuEUvb3ZoZtYyBIM1n4W4aIT9F2ZeX0OVC7unyxJ+BZgFbE13z5b0EHB+RNybtj3qMFM1S0n9E8nw4VqeP43eGuC0z72GIxccwcOrHmHarG7mvu0IujwL0cwyCrSU1HLg6ohYkd0p6WyS3tur8jRSTQ/sRODY8axXZRMniZ7T/4Ce0/+g2aGYWYsqUAI7qDJ5AUTEtyRdlbeRatYd2QJMqaK+mZk1SJAsJ5W3NNkOSe+W9GwOktQh6TzgibyNVNMDuwz4iaSf8fxp9FdU0YaZmdVJgXpg5wFLgC9K2p7um0lyT9l78jZSTQL7NPAMsB/giy9mZi2mKAksIjYCr5c0nWQiB8BDY03aqFRNApsfEfOradzMzBqnKAlsSJqwxn1DazXXwNZLOmC8b2RmZnWU3geWt7QqSRvy1q2mB7YbWCPpX3j+NbBLxwhmHsl8/5cAO4H3DM3xH6augJ8CJ0bEQem+OcAm4J5M1YURsamK2M3MSi9qfB9YvUj6w1EOT8vbTjUJ7L60VGsJsDQilqdz/JcDJ49Q94MkyaryaZxPVfOsGTOz9qOaP9CyjtaRzGwfLuMekreRapaS+ru8dYdImgGcxHPLgawArpE0N72Il617LHAmcD7w9mrfy8ys3bXy0GCFB4FTI+LhygPpahy5VLuU1KuA40lmIgJjrkY/C9g+9AyxiAhJW0nWvno2gUnqApYBF/DcKh9ZL5K0GpgE3ApcNdwN1ZIuBZ4d0jzwwAPzfjQzs0JL7gMrTAJbCRwJvCCBAd/L20ju/qakjwHXAZ8CXgt8EviTvOeP4ePALRFx/zDHtgMvjYiTgTeQLNN/2XCNRMTiiOgZKt3d3TUKz8ys9RVlEkdEfCAi7hjh2MV526lmwPSdwGuA3ohYSHIda3CMcx4CZkrqhGcnaczmucUbh7wWeL+kLcAdwAGStkiaHhF7I+IxgIj4PfCPJEnMzMwyipLAaqWaBLYnIvYAHZIUEeuBo0Y7IU08a4Bz010LSRLgxop6p0XE4RExBzgV2BURcyLicUkz0iFGJE0BziJZUNjMzIYUcBq9pEFJAxXl95K+l85AH1VV0+jTRLIO+Gz6KOpJOc5bBCxPhyB3kUzSoOJ5MqM5FbhC0kAa721A7sUezczaRVGm0WdcTvLv+rJ0+wKSNXcfJZnB/qejnVxNAnsvyVOZLyO5DvYa4N1jnZT21E4ZZv+Fw1QnIrYAB2W2bwFuqSJOM7O21Co9qyqcGREnZbavkHRXRJwk6aKxTh4zgaU3nM2PiG+m218ADiaZ9OLV6c3MWkQLrDJfrWnpXIfHAdK1EYduZN431sl5roF9gucvb//HJFPZfwJ8pJpIzcysfgZDuUuLWAzcLekGSTeQzG/4vKRu4M6xTs4zhHh4RNyW2d4dEV8BkPSO8URsZma1FbTO5Iy8ImKJpDt57pasz0fEb9PX7xvr/DwJbHLF9sLM6xfnON/MzBqgaAkstRk4IPM6tzwJrFPSARGxCyAiNgNIOhA/F8zMrGUULYFJeg3JEoOPpLsOlbQwIn6R5/w818BuIpkG/+y6TOnrG9JjZmbWAop2HxjJNbCzI+KEiDgBOBv4fN6T8ySwq4CngV5JayWtBXpJHqnyyXEEbGZmtRbJfWB5S4uYGhHPTtaIiFVk1tody5hDiOlCvOdKmguckO5eW7mahpmZNVcL9azy6pP0hoj4CYCk1wP/kffkah6nspHMCvJmZtZaCpjAPgCsSFdagmRU8Ky8J1f1OBUzM2tNAQwU54GWAETEXeno3vx01/qIGPMG5iFOYGZmpdBSkzNGJemAil1DTyiZKmnq0Kz3sTiBmZmVRFESGPAkSacxG/DQdpBvoXgnMDOzUgiIsZ7Q2CIioiZjnU5gZmYlUaAeWE04gZmZlUQLLdLbEE5gZmYl4R6YmZkVkhOYmZkVTuAEZmZmBeVrYGZmVkgRzY6gsZzAzMzKIF2Nvp04gZmZlUJxlpKqlWKt/GhmZiMaDOUueUj6gqQtkkLS8aPUu0DSA5I2SVomqatWn2k0TmBmZiURkb/k9C3gVODBkSpIOgK4EjgNmAscClw0sU+ST90TmKR5klZJ2iBptaRjR6krSbdJerJi/xmSfpdm+FuGWcnYzKztRSh3ydde/DwieseodjawMiIeiYgArgPOmeBHyaURPbAlwNKIOBq4Glg+St0PApuyOyR1AzcAZ0bEPOBh4G/rE6qZWTEFVQ8hdkvqzZRLx/nWs3l+D21Luq/u6prAJM0ATgJuTHetAGalDzCrrHsscCbwvysOvQlYGxG/S7evpUHZ3cysSAYHlbsAfRHRkymLmx1/terdA5sFbI+IfoC0e7mViuycXvBbBiwCBiraGC67z5T0ghmUki7N/kbR19dXsw9iZtbqaj2JI6etwOGZ7Tk894DKumqVSRwfB26JiPsn0khELM7+RtHd3V2j8MzMWlz6PLC8pYZWAAskHSZJwMXATTV9hxHUO4E9RKa3lH642bwwO78WeL+kLcAdwAHp1M3pDJ/dn+3VmZlZotaTOCQtkdQL9AA/krQx3X+9pAXJe8Zmkk7IncBG4HGSuQ91V9cbmSPiMUlrgHNJJm8sBHojYmNFvdOGXkuaA6yLiDnp9g+BL0l6WXod7BIalN3NzIpiaBJHTduMWDTC/gsrtpeRXAZqqEasxLEIWC7pY8Au4HxIMjjJ1MuVo50cEU9JuhC4Ne3J/RY4r84xm5kVjtdCrLGIWA+cMsz+C4epTkRsAQ6q2LcSGDXRmZm1O69Gb2ZmheQemJmZFVK7LebrBGZmVhKD7oGZmVnhVLdIbyk4gZmZlUBQ8xU2Wp4TmJlZSbgHZmZmheQEZmZmheQhRDMzK6Q264A5gZmZlUEAA+6BmZlZEfk+MDMzK6Q2y19OYGZmZeEemJmZFVKb5S8nMDOzsnAPzMzMCidwD8zMzApqsNkBNJgTmJlZSbgHZmZmheQemJmZFZIX8zUzs0JyD8zMzArHsxDNzKyw3AMzM7NCcgIzM7NCarchxI56v4GkeZJWSdogabWkY4epc4qkdWm5V9ISSVPSY6dL2p05vk7S1HrHbWZWNINVlDKoewIDlgBLI+Jo4Gpg+TB17gZOjojjgeOAGcAlmePrI+L4TNld55jNzAolgAEidymDuiYwSTOAk4Ab010rgFmS5mbrRcTTEbEv3ZwMTKX9esNmZhMSVZQyqHcPbBawPSL6ASIigK3A7MqKkuZIuhvYAewErs0cPkrSmnQI8pLKczNtXCqpd6j09fXV9MOYmbUyDyE2SURsiYhXAocBU4Cz0kNrgJ6IOBF4G3CxpD8foY3FEdEzVLq7uxsSu5lZKwhF7lIG9U5gDwEzJXUCSBJJ72vrSCdERB9wE/CudHtXROxMX/cCXwdOq3PcZmaF4x5YDUXEYyQ9qHPTXQuB3ojYmK0naa6krvT1ZJKe1m/S7ZmSOtLX04AzgLX1jNvMrIicwGpvEbBI0gbgo8D5AJKul7QgrfM6YG16DWwt8ChwZXpsIXBPeuyXwI+BLzcgbjOzwogq/yuDut/IHBHrgVOG2X9h5vVSYOkI518DXFO3AM3MSqIsPau8vBKHmVlJlKVnlVfLzEI0M7OJqcc1sJyrKTVlxST3wMzMSiJUTeXcNYdWU1ou6WyS1ZROHqbe+nQ1pYZxD8zMrCQGidwlj7yrKTWLE5iZWQkEdRlCzL2aEjlXTKolDyGamZVElZM4uiX1ZrYXR8Ticb710IpJOyX1AN+XtCMibh5ne7k4gZmZlUSV0+j7IqJnjDrPrqYUEf0jraYUEbsyr3slDa2YVNcE5iFEM7OSqPU1sCpWU2rKiklOYGZmJRHKX6qQZzWlpqyY5CFEM7MSSB5oWfu1OHKuptSUFZOcwMzMSsJLSZmZWSHlvbZVFk5gZmYl0V7pywnMzKw0BkvypOW8nMDMzEogWYnDCczMzAqovdKXE5iZWUnkv0G5LJzAzMxKwgnMzMwKyfeBmZlZ4QRVr0ZfeE5gZmYl4SFEMzMrJN8HZmZmhdRu18D8OJWMJ57YzS9X9fLEE7ubHYqZWdVq/TywVlf3BCZpnqRVkjZIWi3p2GHqnCJpXVrulbRE0pTM8QskPSBpk6RlkrpqHednP7OKeXO+yIK33MS8OV/ks59ZVeu3MDOrq6jivzJoRA9sCbA0Io4GrgaWD1PnbuDkiDgeOA6YAVwCIOkI4EqSx1PPBQ4FLqplgL/+1TY+/ck76O8fZM+efvr7B/n0J+/g17/aVsu3MTOrm6GlpNwDqxFJM4CTgBvTXSuAWZLmZutFxNMRsS/dnAxM5blVUc4GVkbEIxERwHXAObWM84c/2ERU/HlGwI9+uKmWb2NmVkdBP4O5SxnUuwc2C9geEf0AaQLaCsyurChpTvo46h3ATuDa9NBs4MFM1S3DnZ+2camk3qHS19eXK8hDDplKZ9fzv4rOrg4OOWRqrvPNzFqBe2BNEhFbIuKVwGHAFOCscbSxOCJ6hkp3d3eu895xzrFMmTyJSZMEwKRJYsrkSfzFO15wuc7MrCUFMKj8pQzqncAeAmZK6gSQJJLe09aRToiIPuAm4F3prq3A4Zkqc0Y7fzwOecn+3H7HeSx463wOn3MgC946n9vvOI9DXrJ/Ld/GzKyu2q0HVtf7wCLiMUlrgHNJJm8sBHojYmO2XnpN7MGI2CdpMvA24Dfp4RXAHZI+ATwKXEyS4Gpq7tyD+cr/PbPWzZqZNUxZElNejRhCXAQskrQB+ChwPoCk6yUtSOu8DlibXgNbS5KorgSIiM3Ax4E7gY3A4yQzG83MLGOAyF3KQFE5/a5Eenp6ore3t9lhmJmNSNK2iOiZaDuTdFDMnHx57vrbnrmsJu/bTF5KysysJNptCNEJzMysJAZUjvu78nICMzMrhfJc28rLCczMrAQCnMDMzKyYBtrseWClnoUoaS/JtPvx6gbyrUfVehx7czj25ihq7N3A1IiYMmbNMUhaDcys4pTtEXHyRN+3mUqdwCZKUm9Rp5k69uZw7M1R1NiLGneraJm1EM3MzKrhBGZmZoXkBDa6xc0OYAIce3M49uYoauxFjbsl+BqYmZkVkntgZmZWSE5gZmZWSG2ZwCTNk7RK0gZJqyUN++hlSRdIekDSJknLJHWl++dI+pmknZLWFSju10n6taT7JN0r6TOSGvIzUIPYT5G0Li33SloiacL3zjQi9sxxSbpN0pONiDt9z4l+76dL2p357tdJmlqE2NNjx6V/V+9PS9VPem9G7JLOr/jOd0i6pRGxF0pEtF0BbgPek74+G1g9TJ0jgIeBwwABK4G/So8dDJwKvAVYV6C4TwCOTF/vB9wx1F4BYt8f6EpfdwDfBj5YhNgzdS4FlgFPFuhn5vRG/ozX4WdmM3Bquj0JmF6E2Iep+1tgYTP+HFq5ND2Ahn9gmAHsAjrTbQGPAHMr6n0YuC6z/Wbgjoo6DfvLXcu4M8euAT5RtNhJku8Pgf9RlNiBY4GfA0c1KoHVIvZmJbAaxX4h8E9FjL2i3quBx0h/gXN5rrTjEOIskiVU+gEi+QnZCsyuqDcbeDCzvWWYOo1U07glHUbym+F36xFshZrEng7d3g3sAHYC19Yx5iETjj0dFlpG8nTygTrHm1Wrn5mjJK1Jh8IuqWO8WbWI/eXAXknfTYfhvippen3DBmr/b8wFwNciYl/tQy22dkxgbU/SAcA/A5+JiLuaHU9eEbElIl5JMuQyBWjI9Ywa+DhwS0Tc3+xAxmEN0BMRJwJvAy6W9OdNjimvTuANJL84nABsA/6hqRFVSdKLgHcANzQ7llbUjgnsIWCmpE5ILqyT/NaztaLeVuDwzPacYeo0Uk3iljSNZPjtOxHRqJsoa/qdR0QfcBPwrnoEW6EWsb8WeL+kLSTXHQ+QtKUBvYEJxx4RuyJiZ/q6F/g6cFp9wwZq871vBW6PiG1pL+hG4I/qGXSqlj/vbwfujYj76hNqwTV7DLMZBfgZz7/AetcwdY7khRdY31dR53QaO4ljQnGTrHx9J3B50b5zYC7PTeKYDHwDuKoIsVfUm0NjJ3FM9HufCXSkr6elPz//vSCxzwbuBw5Itz8MfK8IsWfq/BtwYaN+XopWmh5AUz40zAd+AWwA7gKOS/dfDyzI1PtLYFNabsj8A7o/0EvyqJZn0tefLkDcfw3sA9Zlyl8X5Du/iGQm1t3AvcAXgP2KEHtFW3NobAKb6Pf+vvT7HvreP0G6gk+rx54ee3f6c/Mb4AfArALFPh94CpjWqJ+XohUvJWVmZoXUjtfAzMysBJzAzMyskJzAzMyskJzAzMyskJzAzMyskJzAzMyskDqbHYBZK0pXzdgL7E533UXy92VdRPx9k8IyswwnMLOR/UVErBvakLS8eaGYWSUPIZpV5xhJP00fVHiLpMnNDsisXTmBmY3sG5kn4r4t3Xc88GfAMcChwMJmBWfW7jyEaDayyiHEtwLfjoin0+1fkzyg0syawD0ws+rsybwewL8EmjWNE5iZmRWSE5iZmRWSH6diZmaF5B6YmZkVkhOYmZkVkhOYmZkVkhOYmZkVkhOYmZkVkhOYmZkVkhOYmZkVkhOYmZkV0v8HuSJt1oxOPO0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 480x320 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df[\"log10_R2\"] = np.log10(df.R2)\n",
    "fig, ax = plt.subplots()\n",
    "df.plot.scatter(ax=ax, x=\"Fh\", y=\"Gamma\", c=\"log10_R2\", cmap=\"plasma\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d99f4a9-c95c-4b2a-8bc0-48fb954380bd",
   "metadata": {
    "tags": []
   },
   "source": [
    "## A figure with data from different simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b3f59d22-da54-48b6-9268-6939033f4e85",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "compute mean of spectra\n",
      "tmin =  800.502 ; tmax =  1600.67\n",
      "imin =      100 ; imax =      200\n",
      "compute mean of spectra\n",
      "tmin =  400.051 ; tmax =  800.132\n",
      "imin =      100 ; imax =      200\n",
      "compute mean of spectra\n",
      "tmin =  200.093 ; tmax =  400.025\n",
      "imin =      100 ; imax =      200\n",
      "compute mean of spectra\n",
      "tmin =  133.342 ; tmax =  266.703\n",
      "imin =      101 ; imax =      201\n",
      "compute mean of spectra\n",
      "tmin =  100.004 ; tmax =  200.024\n",
      "imin =      100 ; imax =      200\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEsCAYAAABQVrO3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAxOAAAMTgF/d4wjAABhLElEQVR4nO3dd3hUVf7H8feZyaQ30iGFBELoECQUqYKCiqhg733tu7rorrqWdf3Z27p214LYUFaKBTuIiPTepBNCQnrvmXJ+f0wSE5JAyiSTSb6v55lnyr1z5ySEfHLuPed8ldYaIYQQorsyOLsBQgghhDNJEAohhOjWJAiFEEJ0axKEQgghujUJQiGEEN2aBKEQQohuTYJQCCFEt+bm7AZ0Zh4eHjo0NNTZzRBCCNFGaWlpVVprj8a2SRCeQGhoKKmpqc5uhhBCiDZSSmU3tU1OjQohhOjWJAiFEEJ0axKEQgghujUJQiGEEN2aBKEQQohuTYJQCCFEtyZBKIQQoluTeYSiRbTWWKpsaJvG5GlEKeXsJgkhRJtIEHZhWmusZhvmSitVFVbMlTU3C+Y6z6sqLHW2Wetta2w72n58ZVB4+rjh6WPC08eEh48JT19T9fPq131NeHrXfd2E0SQnIoQQnYcEYSehtcZm0Q2Dp6JOeDUINCvmOvs22FZpRdt0q9pjNBlw9zRi8rDfvHxNmDzdap8rg6Ky1ExFqZmKEjMFWWVUllqwNePz3DyM9QK0bkjWhKjHca97eLmhDNL7FEI4ngShg1mtNoqyy8nPKKMop5yq8uoAO76n1UgvrDkh0hijm6E2oEyeRnx7eNQ+d68TXqY6wWbyNOLu4Xbc6272x+4GDMaW99q01pgrrPZwrA7I+o8t9Z5Xlpkpyi6nqsJ60mMrhT0c64SlydPNHtbV9+6e9va7H/d63dda83UJIbo2CcJWqig1k59RRn5GKQWZZeRnlFGQWUZhdvkJe2EGo6oXPN4BHvV6XvZgcqv3vN7247aZPI0YO8kvd6UU7l5uuHu54R/i1ez3Wa02Kkst9YOz0SA1U1lmoSinnKwUMzZLy/9wcDMZaoOxsdB093TD09eEf4gX/iGe+AV74uljkmuhQnRhEoQnUZhdVh14ZRRklJKfaQ+88mJzvf0MBkVAmBexQ4PpEeFDYLg3AWFeeHi71euVGd06R2h1JkajAW9/d7z93Vv0PqvZRlX19c6qCgtVFdY/euB17mt641UVDfctya+kqsJywlA1eRrxD/bEL9gejv7BNSFpv3f3lP9GQrgypXXrTsd1Bz18Q/XjV31W7zUPbzd6RHgTGOFDj3Bv++Nwb/xDvTpNz0y0XE2oVpVbKS+uoii3nKKcCopzyinKraAop5ySvMpGT197+pjqBaN/iBd+wZ54+Zrw8HbDw8uEu5eclhXCmZRSaVrrqMa2yZ+yJ2A0GRh+RnR14PnQI8IbT9/Of5pM22xosxldVWW/1Tw+/r76sa2qCqrv/9hubnJ/XVUF2oby9MLg6Yny8sTg5Y3ByxPl6YnBywuDl5d9e4PX7I+Vh0en+j4aTQa8TO54+UJAqBcRfQIa7GOz2igpqKQ4p8Iejrnl1Y/toZmVUlw7orYxbh5GPKpPHdfee9d/XHN//EAiN3dDp/p+CdGVSI/wBKKionRr6hFqrbFkZ2NOTUNXVjQaPH+8Zkabq04aPA1eq36PzVwnuKq3YbG0w3fDwZRCedmD1B6m1Y+9vGof1wZszeOaYPXywugfgEe/eNx790a5dY6/56xmG8V59t5jUW4FlWX2a5pV5RYqyy1UlVXf13luMduadWyjyVAnHN3qj7L1NeET6GG/BXjgE+COm7uxnb9aIVyL9Ajbia2qCnNKCpWHDlF16DBVhw9ReTiZqkOHsJWUOOxzlLs7ymSy39d97OmJm79/nddMKJM7yt2Ewd0dTPb72v2rtylT3dcaOW7NsaofGxrZD4MBXVGBrbwcW3kFuqIcW0UFtrJy++PyCmwV5ejyOo/L7PvUbi+v3l5h327Nz8ecnm7fp7wcmvFHmnJ3xyM+Ho/+/fHon4Bn//549O+PW1CQw77/zWU0GQgMt58qby6rxWYPxjILVRX2+8oyyx8DhBoZQJRztILKshP/sePh42YPxdqAdD/uuQfe/iY5XSsE0iM8oZoeoSU/n6pDh+yBVx10lYcPYT6aCrb6f9EbQ0PwiOuDe5843KNjMPh4NwydRoPIVCeoqsPH3R3c3LrlKTGtNbqy0h6WFRXVwVlW+9ial0vl/v1U7N1L5d59WDIz673fGBqCZ4I9FD37J9iDsk8f+/e0C7DZNJVl9mAsLzZTWlhJaUElpYVV9vuCytrXLFWN9zqVAi9/d3wDPfAO8MA30AOfQHd8e3gSEOZNYJiXS1wKEKI5TtQjlCA8gZ7e3vqXYcOxFhTU32Ay4d47Bo+4ONyrQ8+jTx/c4+Iw+vk5pa3dnSU/n8q9+6jct7c2HCv370dXVv6xk5sbHn361A/HhP64hYV22V/2WmuqKqyU5leHY1H1fUFVvbAsLaxqdNqPh7cbAaFetcFoHw3tjX+ITCsRrkWCsJUi3D302osusgddnD3oPPrEYYqK6jTXpUTTtNVK1ZGU+uG4dy/mtLR6+xkDA+udWjVFR2MMCMDo64vB3x+Djw/K0LVPIdpsmvLiKsoKqyjOraAgq4yCrDIKs8opyCqjrLCqwXtMHkb8gv8YJetf/TgwzJvACG8MshKQ6EQkCFuptYNlROdmLS6mct++euFYuW8ftrKyxt+gFAY/P4x+fhj8/WsD0v7cD6OfP0Z/Pwy+9Z8bg0NwCw7qEn80VVVYakOxMKuc4tzqaSW5FZTkVjSYVuLmbiAkypfQaD9Ce/sRGuNHj54+MsVIOI0EYStJEHYf2mbDnJZGxZ49WNIzsJYUYysqxlpcVH1fjK24+r6oCGtxMVhPvjQcBgPG4CBMoWG4hYXhFhpqvw8Lwy0sFLfQ6vvgYJTRNUd62mya0oJKezjmVJCXXkp2SjHZKcX1BvUY3BQhkb6Ex/oTFudPeKw/gWHesoas6BAShK0kQSiaorVGl5VhLS7GWlRkD8miImwlJViLirAWFmLNycWSnY0lK8s+nSY7G8zmxg9oMOAWHFwnJKsDMiwMU50ANQYFucxpWq01xXkV5KSUkJVSRHZKMVlHiqko+eN74O5pJLS3P+Fx/vSKD6Rn3wDcvVy/By06HwnCVurKQWixWaiyVmG2mamyVlFlq7LfN/Fa3edmm/0XWbBXMCGeIYR6hxLiFYKPycfJX1XnprXGWlCAJas6HKsDsu69OTsLS3ZO04Hp5oZbSEhtUJpqQjMiAlNET0w9I3CLiMDg4dGxX1wzaa0pzq0gM7mIrOQiMpPtAVkzslUpCIn2o1dCIL3iA+nVLxBPH5OTWy26AglCQCl1NzBTa31Gc9/jqCDUWtuDx1Y/VMxW88lfqxtINduPC6Um39/IcWq223TzJnK3hJebF6Fe9lCsuYV6hxLsGVwbliFeIfTw6IHR4JqnATuCttmwFhb+EZbVN3NWVv0Qzclp8vSsMTgYU0QEbj3rB6SpZ0/762Fhnebapc1qIzetlGP7C2pvFaV//CEQ1MuHnvGB9IoPoGd8IH5Bnk5srXBV3T4IlVIm4G0gqj2CsNJayZGiIxwuPMzhwsMkFyVzuPAwqcWpVFgqMNvM6BOtveUABmXA3eCOu7H6Vv3YZDT98brhuOfH72swNfmah9Gj3nObtpFbnkt2eTY55Tm1t+zybHLLc8mryGuyrUZlJMgzqMmwDPUKJdgrmFCvUDzd5JdeU7TVijUvD3NGJuaMdCzpGZgzMrBkpGOueZyV1XhYGgz2nmVNUEbUBGUEpshI+1QgX9+O/6IAbdPkZZRybF8B6QcKOHagkNKCP6bB+AV5EjM4iNihIUQN6CGr6IhmkSBU6gYgDfhba4NQa01uRS7JhckcLjpcG3qHCw9zrORYg6AL8w4j2i8abzfv+iFU/djD6NEglGqD6LjXjg+hxoLKzdA5/rqvYbaZySvPqw3HuiGZXZZNTkUOOWX28KyyNRyaX8PX5FsbliGeIYR4h9QLy5rHAR4BGJRrXDvrSNpiwZKTg/lYuj0gG4Rmhr1n2cjvAbewMNz79LHPke3TB48+cbj37WvvTXbg/MGa06k1oZi2L5/CrHJ7G00GogYGETs0mNhhIfgEdM5TwsL5XC4IlVIvA+cBvYERWuutdbb1A+YBIUAhcJ3WetcJjmUA5mutL1VK/dSSIAwKD9K3L7y9NviKq4rrbXc3uNM7oDex/rHEBcTZb/5xxAbEyvWyZtJaU1RV1KBHeXxYZpdnU1RV1ORx3JQbQV5B9U7N1vQq656qDfEKkV7mcXRVFeasbCyZGfaeZOrRP5YNPHSowbQSg4/PH8vaJSTYFydISMAY0HCh8vZSkFlG8o4ckrfncOxAYe1iAGG9/YgbHkLssFCCI31kwr+o5YpBOAk4BKwCZh0XhMuBD7TW7yulLgLu01qPUkoNAl4/7lDfAQcAD631xy0NQlOQSQ/49wCCPYOJDYitDbq4AHvY9fLpJde6OlCVteqPHmZZDrkVufV6nHVP1VpsTa/FWdPLbOoW5h1GT5+e+Lo759RgZ6K1xpKVZV9W8OCh6vuDVO7bhzU/v96+bhER1YsSDMBz0CA8Bw+yLz7RzmFUUWomZXcuydtzObIzl6py+7+9X7An/ZLCSBgdQXCk/Ft2dy4XhDWUUsnUCUKlVBj2YAvSWluU/X9YOjBBa32giWM8CJwGWIEk4H6t9TtN7DsHmFPz3DfANzI1M5UAj477S1e0Xd1eZr2QrOll1gnNwsrCJo/j5+5HT5+ef9x8e9LLpxcRPhH08u1FiFdItz0dq7XGmpNDxd59VO6zL0pQsX8fVfsP2CugVDP4+eE5cGBtMHoNH44pOrrdwtFqtZG+v4DD23M4tCWbknz7tcXgKF8SRoeTMCoC3x5y+rQ76kpBOBL4RGvdv84+67GH2/JmHK9FPcKuPH1C2FVZq8irqL6WWZZNdnk2WWVZpJemk1GawbGSY2SUZTTaw3QzuBHuHU4v314NArOnT0/CvMO63SlybTZTefgwFbt3194qd/9e7/SqW0QE3qNH4TN6NN6jR7dbMGqbJv1gAXvXZXJwcxaVZRaUgpjBwQwc35PYoSEY3brnHzLdkQRhK0kQCgCbtpFTnkN6aTrpJemkl6bbA7I0w/649FiD68c1vNy8Gj39evxgnx6ePTrdgCdH0TYbVUeOULF7N+WbNlO2YT2V+/84geMWEYF3UhJeiYl4JSbi2T/BXurLgaxmG8k7c9izJoMjO3PRNo2Xn4mEMREMmRjZotJZwjV1pSBs8anRtpAgFM1VUlViD8o6PcnagT/Vp2fzK/KbnEajUPTw7EGoV2jtiOMY/xii/aKJ9osm0jcSd2PXKCEFYMnNpWzDRsrWr28QjMrDA8/Bg/EaPhzvpJH4nHoqBm/HBVVpYSV712bw++p0CjLtPdXoQUEMnRxJ76Ehslh4F9VlgrD6tRXA+3UGy9yvtU5qj8+XIBSOZLFZak/D1pt7WZb9x8CfsmwyyzJrV++poVBE+EQQ4xdDlF9UbUjWPHf1U7CW/HwqduygfOs2yrdto3z7dmzF9l628vDA59RT8T19Kn5TpuAWEuKQz9Rac2xfATt+SeXQ1hy0TeMX5MmwqVEMGt9LlnrrYlwuCJVSbwHnABFALlCstY6v3tYfeB8IBoqA67XWO9qjHRKEwhmsNitZZVkcLT5KSnEKR4uP1ruVmksbvMfH5FPbmwz1tt+HedV57B1GqFeoy/Qqtc1G1aFDlK5eTfHynynbsMG+MIBSeA0fju/kSfhMmIjn4EEOWXu1JL+S3avS2PnrMcqLqnD3NDJoYiTDpkTJSjZdhMsFYWchQSg6G601eRV59YIxrSSNrLIsssvsA32KzY1frwQI9Qqlt39vevv3Ji4grvZxlG8UJmPnXdPTWlhIycqVFC9bTumvv2Irtf8xYAwKwmf8eHwnTsB3ypQ2F8a2mK3sW5/J1h9TyM8ow2BQ9B0ZRuIZ0YT19nfElyKcRIKwlSQIhSsqM5eRU55DVlmWPSCrR8JmlWWRWpxKclEyJeaSeu8xKiORvpHEBsQS6x9LQo8E+gf1p09An07Xi9RVVZRt3Urpr79S8usqKvfsAUB5euJ/5nQCLrgQ71FJbeopapvmyK5cti07Suoe+3zJnvEBjJjem9ihwTJR3wVJELaSBKHoimqWCzxSdIQjRUdILkrmSKH9cUpxSr3rk0ZlJC4gjoQeCfVuYd4du8zaiZgzsyj5ZQWFS76gfPNmAEzR0QReMJuA887DFBnZpuPnpBazbdlR9m3IxGbRhET7MuqcOOKGh3Sa74E4OQnCVpIgFN2N1WblaPFR9uXvq3dLK0mrt5+fux99A/rSN7D6Vv3Y2QFZeegwhYsXUbBkCdbsHAC8k5LwP/88/M88E6N/609vlhZUsuWHFHb+mobVbCMk2pdxF8YTPSDIUc0X7UiCsJUkCIWwK64q5kDBAfbm7WVf/j4OFhzkYOHBBivz+Jp8SeiRwKDgQQwJGcLQkKFE+7XfSjJN0RYLJatWUfTllxQvW46urES5u+M7dSoB587Ed+JElHvrTvmWFlay5ccUdq5Iw2qx0ScxlHEXxhMQ6uXgr0I4kgRhK3WpILTZQNtAW+33Nmud5/q453W317nV3cdggoBI8JTl57qrmlOshwoOcaDgAIcK7ff78vbVG7Dj7+7P4ODBDAkZwrDQYYwIG9GhyxZai4sp/uEHCr/4krL16wEwBgTgd9ZZBJx/Pl4jElsV1EU55axedICDm7MxuhkYfno0I8/qLdMuOikJwlZqcRCW50PO/urbPihIAZv5BCHUSMDUPtdNhJbtBMF1gve0F48ACIyGgGgIiPrjcWCM/blPGDhgeLtwHTZt40jREXbm7GRX7i525uxkT94eKq32dT8Vin49+nFK2CmMjBjJyLCRhHqHdkjbzMeOUbh0KUVffkXl/v0AeA4dSvCNN+A3bRrK2PJF9NP25vPrgv3kppXg5WdizHl9GDiuJwaj/Nx3JhKErdRoENqs9oCrCbvcOsFXmt3wIMpQfTPa7w3GOq8ZjnutZp+TvKfe85a+p87rTb7HCEo10q4677VUQGGq/VaQAkVp0FjFB6OHvecYEF0nMOs89o8Et841KlE4ntlm5kD+AbZmb2VT5iY2ZW4ipzyndnusfyzjeo1jfOR4ksKT8Da1/5JnFXv3kj9/PoWLl6ArKzFFRRF03XUEXjC7xSvZ2GyaPWvSWfvFIcqLqgiO9GHSZQn06tejnVovWkqCsJWieobr1O9fsYdczj574OUeBGtl/R09AyAkofrWD4L72R/3iO0+v+RtVijJhIKjUHjUHo6FqdWPq1+rKmnkjQr8IuqEY1SdHmX1c0+Zv9XVaK05Wny0NhTXZawjozQDsNf5PCX8FCZETmBC5AT6BPRp12uMltxc8j/+hPxPPsFaUIAhIIAel15KjyuvxBQe1qJjVVVY2PzdEbb8lILNokkYHc64C+LxCZSKF84mQdhKUf4GnTqnZoKugh69/wi84Pg/HvuE2HtQomlaQ0XBH6FY05OsfXwUSrMaf69nAATE1D8FW/NvEdS3+/yx0YVprTlceJhVaav47dhvbMzYSJWtCoBI30hOiz6NSVGTSApPard5jbbycgqXLCHv/XlUHTkCJhMBM2YQdP11eA4Y0KJjFWSVsWrBfo7szMXkaeTUWX0ZPClS1jF1IgnCVooKDdSpK96v/oXbB0yy1FK7MpdDYVp1OB5tGJqNnX5VRgiKg5D+EJpQ5z4BPNq2yohwnnJLORsyNrAqbRW/HP2FY6XHAPB282Zcr3FMjp7MpKhJBHk6fuqCttkoWbGCvPfmUrZxIwA+404l6Lrr8JkwodkT9bXWJG/PYeWn+yjJryQ8zp8pVw2QIsFOIkHYSl1q1GhXUPf0a34y5OyF7L3209Z5hxqGpH+kPRBD+9e57y89eBejteZAwQF+Sf2FX47+wrbsbWg0CsXw0OGcFn0ap0Wf1i6nUMt37CTv/fcp+u47sFrxGDCA0Lvvwnfy5GZ/VlWFhfVfHmb7z0dRSpE4PYZRM2Jxc2/5wBzRehKErSRB6EIsVfYwzNkL2fvqhOR+sJTX39erx3E9yAH2x/5RMsLVBeRV5PFr6q/8kvoLv6X9RpnFXkqpb0BfZvSZwdlxZxPtF+3QzzQfO0bevA/I//RTdGUlXiNHEvbXu/FOan7hm6wjRfz80R5yjpYQEOrFlKsGENlfBtN0FAnCVpIg7AJsNvvp1Zx91cFYJyjL8+vva/L+Y7BTUB/7KdegPvabT6j0IjuhKmsVGzI2sDxlOT8e+ZH8Svu/6bCQYbWh6MjTp+bMTHJee52ChQvBasX39NMJu/cePOLimvV+m9XG1mVH2fDVYSwWGyPOiGHMeX0wmuQPsPYmQdhKEoRdmNZQmlMdjHvq9CL3QfGxhvubfKpDMdZ+36NOSPr3sk8tEU5ltplZe2wt3xz+hmUpyyi3lONmcGNy1GRmx89mfOR43AyOmexelZxM1r9fovj778HNjR6XX07oHbdjDAxs1vsLMsv4ce5uspKLCIn2ZfqNg+kR4do1JTs7CcJWkiDspqpK7dcg8w5B3uHq+0OQf9g+cOf4BQqM7vapMrXhWH0fEG0PSZn+0eHKLeUsT1nOFwe+YG36WjSaEK8Qzu17Lhf3u5hof8ecOi3bvJnMZ56hYtt2jIGBhN59N4EXX9SsiflWq42NS5PZ9G0yRpOBSZf1Z8CpEbKQdzuRIGwlCULRgKXKPoK1bjjWBGZ+sn0loeN5+NsD0T/Sfh8QVed5pH3BARnh2m7SS9L58uCXLDmwhNQS+//nsT3HcnHCxUyJmYLJ0LY6jFprir5eStazz2LJzsZz0CAi/vkIXsOHN+v9x/YX8MO7uygtqCRhTDiTL++Pu6cs0+ZoEoStJEEoWsRmtfcYa8KxMA2KjkFRqv2+MK3hwJ0aHv5/BGVgTPUo1+rBPP6RMojHAWzaxvqM9SzYu4CfU37Goi2EeIVwSf9LuCThEoK9gtt0fGtJKTlvvE7evA/AaqXHVVcRdvddGHxOfsqzosTMsnm7Sd6RS1AvH86+ZSiB4e2/uk53IkHYShKEwqG0tg/QKaoOyMLqgCxKq/NaI2FZM4gnpH+dgJTFBNoipzyHxfsXs2DfAjJKM3A3uDOjzwyuGngV/YP6t+nYlQcOkP7wI5Rv2YJbr570/Ne/8J048aTv01qz5YcU1i45iMnDyOnXDaJPYseswdodSBC2kgSh6HBa2+dK1izrVzOIJ2e/PSzrql1MoHqlo5pbSD8Z5dpMFpuFZSnL+Gj3R2zN3grAxMiJ3DzsZhLDElt9XG2zkf/JfLJffBFbWRkB559H2P3349bj5NMljv6exw/v7KKi1EzSjFhGzYyTFWkcQIKwlSQIRadSWfzHAu+100GaWEzAwx+C+9qnggTH2x/XTA1xl1NujdmZs5O5O+fy45Ef0WjGRIzh5mE3MypiVKsHsJiPHSP9X/+i9JeVGIOCiHjoQfxnzDjp+4rzKvj2zR1kpxQTMziY6TcOwsO7bdcyuzsJwlZqTRBarDZS88s5mF1Can45Vlv9729j3+3m/hs0tptu5IiN7+fY4/m4G0mM6cHgXv6YpNyMc1ktUHAEcg9ULwx/4I9bcfpxOyt7KIYPsd8ihkD4YPsIV+lBAnCo4BDv7nyXpYeWYtVWhocO5+ZhNzMxcmKrArFmME3mk09izc/H7+yziHjkkZP2Di1mK7/M38ee1en06OnDzDuH4R8sxX9bS4KwlU4UhEUVZg5ll3Iwq4SD2SX2x9klJOeWYrZ2n++pp8nAiOgeJMX2ICk2iFNiAvHzlL9cO43KYnvFlJpgzNoNmbvsr9X9o8cjwB6I4YOh5zCIGmW/JtmNB+kcLT7Kezvf44sDX2C2menfoz+3Dr+V02NOb1UgWvLyyPjnoxT/+CPGkBB6PvYYflOnnPA9Wmu2/nSU1YsO4OXnzsw7hhHWW6bjtIYEYStFRUXpNdv3cTC7hIPZpRzKLql9nF1cvxSTQUFMkDd9Q33pE+pD31BfYoK9cW+kt9T4/6GGLza2X2Nvbew/ZeP7NXa85n3u8fJKq9h0JJ+NR/LYfKSAcrMVsH8f+kf4M6o6GJN696BXoPwV2+lUlULW75C5EzJ22sMxcxdUFv6xj0cARI2EqNH2YIxKAq9ApzXZWTJLM3l/1/t8vu9zKqwVDA4ezF9O+Qun9jy1xYFY0zvM+L//w1ZUROCllxJ+399PWv/w4OYsfpy7G2VQzLh1KNEDHb/YeFcnQdhKJv8QHXn7vHqv+Xm40SfMl74hPvQN86VvndDzcOueq4uYrTZ+Ty9iY7I9GDck59f7QyEy0Ku2x5jUuwcJ4X4Y5eJ/56O1fTm6Y1vg6HpI3Wh/XLf+Zkh/iB4FMeMgdjwE9u42p1RzynN4Z8c7LNi7ALPNzOiI0cxJmsPg4MEtPpY5I4NjDzxA2Zq1uMfG0uuF5/EafOLjpB8oYOnr2zFXWZl2/WDiR7asVmJ3J0HYSj5BYfrv836u7eXFh/oS6uchKz+chNaao3nlbEjOY+ORPDYm57M/64+ivH6ebozs3YOk3vZwTIwOxNPUPf+I6PQsVZCxA1LXV4fjBntY1vCPgt7j7LfYCfaBOV38/8exkmO8se0Nvjz4JVprzu17Lnedchdh3i0LJm2zkTfvA7JffBGUIvyhBwm8+OIT/n7JSS3hq5e3UlZcxWlX9GfwxMi2fjndhgRhK8moUcfJrz6VuuFIHpuS89meWkiV1b5UmcmoGNwroN7p1GBfqejdaRUdgyOr4chv9vvsPX9s8wmDuInQ70zoNw28u+4pvL15e3l+4/OsTV+Ll5sXNw+7mWsHXYvJ2LJr5OU7dpJ2112Yjx0jYNYsIh79JwbPpmufFmaX8+XLWynKLmfM+X0YeVZv+eO8GSQIW0mCsP1UmK3sSCu0n05NzmPjkXwKy/9YnqxPiE+906lxIT7yn72zKs35IxiTf4PMHfbXlcF+fTFhOiScBWGDulxvUWvNytSVPLfxOY4UHaFvQF8eGvsQSRHNL88EYC0oIO3vf6d05a94Dh1K1KuvYgpvuodZWljJV69sIze1hBHTYjj1gr7y/+MkJAhbSYKw49hsmoPZJWyoE4wpeWW124N93O3B2DuIpNgeDO4VgLtb9x3R2KkVZ8L+H2D/93DwZ6iqPi3uHwUDZsCgWRAztktV7KiyVvHezvd4e/vbVNmqOL/v+cxJmtOiElDaZiP75ZfJffMt3MLCiHrtNbyGDmly/8oyM1+/uo2MQ0UMnxrN+IvjJQxPQIKwlSQInSuzqKJ2AM7G5Hx2HSukZlpmzbSN0XFBjI4LYkRMIN7uslBxp2OptPcW930P+76zr8MK4BsOA8+DwbO7VCimFKXw+NrHWZO+Bn93f+4eeTcX9rsQg2r+H22FX31N+oMPgsFAr6efxv+sM5vct6rCwtevbiP9QCFDT4ti4qX9JAybIEHYShKEnUtppYWtRwvYkJzHhuT60zbcDIohkQGMiQtiVKz9FiArcXQuWtunaOxaDLuX2Oc1wh+hOOQCiB7r8nMXtdZ8f+R7nl3/LNnl2SSGJvLkhCdbVPqpfPt2jt5+B9acHELvvovgW25pMuDMlVaWvr6NtL0FDJkcyaTLEiQMGyFB2EoShJ2b2WpjZ1ohG5LzWH/YPm2j5jqjUtA/3I/R1cE4Oi6IcP+mByCIDtZUKPr1sgfi4Asg8hSXvqZYUlXCq1tf5ePfP8bLzYsHRj/ArPhZzQ4p87FjHL39Dir37CFg1ix6PvYvlHvji6ybq6wsfW07aXvzGTIpkkmXSxgeT4KwlSQIXYvNptmXVcyGw3msO2wPx6w68xl7B3szOjaIUXFBjIkLIibIW35ZdAZa26do7FoEOxfZl4sD+xzFoRfBiKvshY5d1Lr0dfxj1T/IKstiavRU/jnun82+dmgrLSVtzj2U/PIL3mPHEvXyfzD6N76yjLnKyjevbyd1Tz5Dp0Qx8RI5TVqXBGErSRC6Nq01KXllrD9c02PMIzn3jwE4YX4etdcYR8cFkRDmJ6v8O5vWkLYZdi609xaLj9lfj5sMp1wDA88FN9ebWlNYWcjjax/nu+TvCPYM5rHxjzEpalKz3qstFjKeeIKC+Z/i0a8f0e++gyms8RGl5iorS1/dRtq+AhLPiGbchTKApoYEYStJEHY9mUUVtaG4/nAeezOLaxcVD/AyMSq2B2P7BDNtUDi9g09eUFW0I5sNDq+ATfNgz1KwmcErCIZfBkk3Qki8s1vYIlprvjn8DU+sfYJiczGX9r+Uv436Gx7Gkwe71prcd94h+4UXMcXEEPPee7hHNT6Zvu4AmtHnxjHqnDhHfykuSYKwlSQIu76Csio2JuezIdl+OnVnWiGW6qGpCeG+TBsUzhkDwxkeFSi9RWcqyYbtn9pDMXe//bX4aTDmVug71aUG2GSUZvDgqgdZn7Ge/j3688JpL9Dbv3ez3pv/6Wdk/OtfuIWFETP3PTz6NH7KuLLcwhf/3kJ2SjETL+3HsCnNH6jTVUkQtpIEYfdTVmXhtwO5/Lg7g2W/Z5FbWgVAqJ8HZwwMY9qgcMb1DZEl4ZxFa0heBevfsvcStc2+rNvoWyDxcvDwc3YLm8Vqs/LfHf/lja1v4OXmxWPjH+PM2KanSdRV+PVSjt13H8aAAGLmvodn//6N7ldeXMWi5zdTkFXGmTcN6fZrk0oQtpIEYfdmtWm2Hs3nx91Z/Lg7g4PZpQB4mYxMSghh2qAIpg4II8in8ZF8op0VpMD6t2HzPKgotBcjTrwSRv/JXnPRBaxLX8d9K+8jtyKXqwddzV9H/hWT4eTTfoqXLSPt7r+ivL2JeedtvIYObXS/otxyFj6zicoyC+ffnUjP+EAHfwWuQ4KwlSQIRV2Hskv46fdMftydyaYj+di0vexUUu8gzhgUxrRBEcSFyHXFDldVCtsXwPr/2ustAvSbDmNugb6nd/opGFllWdyz4h62Zm/llLBTeOG0FwjxCjnp+0p++43UO+5EmUzEvPdek6vQZKcUs+iFzRjdFBffn0RA6IlLPnVVEoStJEEompJbUsnyPVn89HsmK/fl1E7s7xvqw7RBEUwbFMaI6B5yXbEjaQ3Jv8K6t2DvN/bTphHDYPLfof85nfo6otlq5oVNL/Dx7x8T4RPBy1NeZmDwwJO+r3TNGo7eehvKw8MehkMaL+WUvCOHb17fTmC4Nxfel4SHV/dbhUmCsJUkCEVzVJit/HYgp7q3mEVOiX3uYpifB2cPieDsoT0ZFRskNRg7Uv4RWPsGbJoLlgoIGwyT/wYDz+/UgbjkwBIeW/MYRmXk8QmPN+u6YW0YenrS+4N5TV4z3PpTCr99foDoQUHMvGMYhkaKhndlEoStJEEoWspm02xLLeCH3Zl8uyO9dt5iiK8HZw0JZ8bQnoyODcKtm/0ScpriTFjzCmx4F8xlEDoATnvAvqRbJw3ErVlb+euKv5JTnsMtw27h9sTbT7pWacmvq0i9/XYMfn70/vADPPo2vEaqtWbFJ3vZ/esxhp4WxaTLEtrrS+iUJAhbSYJQtIXWmt/Ti/l2ZzpLd6RzqHqwTbCPO2cOiWDGkJ6M7SOh2CFKc2D1K/bBNeZSCB8KUx6A/jM65TXEjNIM7vr5Lnbn7ub0mNN5csKTeJtOfG2vePlyUv9yF25BQfT+5JNG5xlarTa+enkraXsLmHRZAkNPazQXuiQJwlaSIBSOorVmX2YJS3ek8+2OdPZn2UsT9fA2ceZg++nTcX2DMUkotq+SbFj9H1j/DljKodcIOPNJ6D3O2S1roMJSwSOrH+Hbw98yKHgQr53+2kkH0RR98w1p99yLKSaa2I8/xi2k4f4VpWY+f2YjRTkVnPuX4UQP6LrFk+uSIGwlCULRXvZnFvPNjgy+2ZHO3sxiwL6yzfRB9tOn4+KD8XCTuYrtpjgTVv0bNr4L1ip7jcRpj0GP5k1s7yhaa97Y9gZvbHuDSN9IXj/9dfoEnnjd1fxPPyXj0X/hMXAgvT/8AKOvb8N9Mkr5/JlNKAUX3ZdEYHjXH0kqQdhKEoSiIxzIKuG7neks3ZHB7+lFAPh5uHH6wDDOGhLB5IQwvNwlFNtF3iH44WHY8zUYPWDcnTDhr51uYv7i/Yv515p/4evuy5tnvMmQkKYL9gLkvPEG2f95Ge9TxxLz1luNVq04siuXpa9uIzDCh4vuG4m7Z9ceSSpBCCilbgEuAYzANK21+WTvkSAUHe1wTinf7kzn+50ZbEstBOxFiE9LCOPsoRFMGRCGv6fUWXS4Q7/Adw9A1i7wCYUp/4AR14Cx84TDr6m/MmfFHAzKwCtTX2F0z9FN7qu1JuOxxyiY/yn+M2fS69lnUI0MDtr8/RHWLD5IfFIY028c3KUX6O72QaiUigH+rrW+syXvkyAUzpRWUM73OzP4bmcGG47koTW4Gw2Mjw/m7CE9OWNQuKxq40hWC2z9CJY/AaVZENIfpj8OCdOd3bJamzM3c+eyO6m0VvLCaS9wWvRpTe6rrVZS77qLkp+WEXzbrYTddVfDfbTmu7d2cmhrNhMu7sfw07vumqQShErdBIwHegMrtdaPNud9EoSis8guruSH3fZQXHMwF4tNYzQoxsQFcfaQCM4d3otAbwlFh6gsht9eto8ytZTbF/c+80kI7RzTDfbm7eXmH2+mqLKIpyY+xVlxZzW5r628nCNXX0PFzp30fOopAmfParBPVbmFBU9toDi3ggvuHUl4XOP1Dl2dSwahUupl4Dzs4TVCa721zrZ+wDwgBCgErtNa7zrBsf4BxGutb1BKzQP+Xfd4TZEgFJ1RQVkVP/2exXc7M1i5P5sqiw13o4Fpg8K5KCmKSf1CZfK+IxSmwbLH7FUvDG72hb0n/x28Ap3dMg4VHuJPP/yJ7LJs/m/8/3F+/PlN7mvOyiL50suw5OQQ8847+IxpeEo1+2gxC5/ZhLe/O5c8OApPn653+t1Vg3AScAhYBcw6LgiXAx9ord9XSl0E3Ke1HqWUGgS8ftyhvgOKgRKt9Tyl1I3Vjz87WRskCEVnV1Jp4YddGfxvYyprDuUCEOHvyQWnRHLRyCj6hDYcMSha6OgG+O4+SNsEvuEw4zkY1HTwdJTU4lRu/P5G0kvTeWLCE5zb99wm963Yu48jV1wBRiOxn85vtHzTrl/TWPHxXvokhnLWLUO63PVClwzCGkqpZOoEoVIqDDgABGmtLcr+r5UOTNBaH2jiGKcAV2it71VKvQTM11qva2S/OcCcmucBAQGRBQUFjv2ChGgnR/PK+N+mVBZuSiWtoByApN49uCQpmhnDeuLr0XkGfrgcmw22fgw/PAQVBTBgJpzzAvhFOLVZqcWp3PD9DWSUZvDUxKc4p885Te5b8usqjt56K6bISOIWfIYxMLDedq01P7y7iwMbs7rkZPuuFoQjgU+01v3r7LMeuF9rvfwEx3kJSAT2aa1vbs5nS49QuCKbTbPmUC7/23iUb3dmUGmx4WUyMmNoTy4aGcWYuCBZDLy1ijPh27/D7iXgEQDTHoVTrnPqcm1Hi49y/XfXk1Oew79P+zdTYqY0uW/exx+T+X+P4zNxItFvvoEy1p+WU1luYcET6ykpqOSSB0YRHNl1zihIELaSBKFwdYXlZr7efoz/bUxl69ECAKJ6eHHBiEhmnxIlZaNaa89SWHovFB+D6LFw7ksQdvJqEe3lcOFhrvvuOkqqSnjtjNcY23Nso/tprUn/x4MULl7c5EjSzMNFLHxuE0E9vbno/iTcukgR6q4WhC0+NdpaEoSiK9mfWczCzWks2ZJGRlEFACN79+CCUyKZOawXAV5db4BEu6ooguWP2+sgGtxg4hyYeC+4OWf07p68Pdzw3Q1YtZV3z3y3yUn3tspKjlxxJRW7dhH937fwnTSpwT4blh5m/VeHGX56NBMu7tfeTe8QXSoIq19bAbxfZ7DM/VrrJEd/tgSh6IqsNs3qgzks2pzGtzvTqTDbcHczMG1gOBeOjGRSv1BZCLwlUjfBV3+BzJ0QPgTOfw16JTqlKVuytnDzDzfj6ebJvLPmNbkcW1VqKocvuBClFHGLF2Hq1avedpvVxuIXtpBxqJDz/pJI9CDXX4/UJYNQKfUWcA4QAeQCxVrr+Opt/YH3gWCgCLhea73D0W2QIBRdXUmlhW92pLNocyprD+UBEOrnweWjorl8TAw9A7yc3EIXYamCVS/CyufsBYIn3mOfamHs+F72ytSV3LX8LoK9gvnknE8I8w5rdL/i5T+TevvteA4fRuyHHzZYhq0wu5zPnliPycPIZQ+PxsvXteepumQQdgYShKI7OZpXxuItaXy24ShpBeUYDYppA8O5+tTejOsb3OWG07eL9O2w5HbI3AFRo+CCtyEorsOb8dXBr/jHqn8wMGgg75/1fpMlnLJeeIHct9+hx1VXEfHQgw2271mTzrJ5vxM3PISzbx3q0j8DbQ5CpZQvUKa1timlBgMVWuuDDm5npyNBKLojq02zYm8WH6w5wi/7sgHoE+rD1WN7c8EpUXIt8WQslfaJ+GteBXc/mPlvGHZxhzfjjW1v8PrW15kcNZn/TPkPRkPDQS/aYiHluusp27iRyBdfwH/GjPrbtebHd3exf2MWky9PYMhk151S0aYgVEo9BowD0oACIBL7BPWjWutHHNvUzkWCUHR3R3JL+XhdCgs2HqWgzIyXycisEZFcNTaGwb0CnN28zu3gclh8K5Rkwoir4exnwb3jyh1prXnot4f48uCXXDnwSu4ffX+j+5mzsjh8wYXo8nLiFi/CPSam3vbKcgufPb6esqIqLnlgFEG9XHOkcVuDcJXWeoJSygjsrpm2oJRaobU+zeGt7UQkCIWwqzBb+WrbMT5ce4Tt1VUxRsQEctWY3pwzrCeeXWSIvcOV5sCim+HgMggdCBe/D2EDOuzjzVYzN/94MxszN3L/6Pu5cuCVjTdz9WpSbrwJz0GD6D3/EwzHXS9MP1jI4uc3ERTpy8X3JWE0ud5gqhMFYbO/Gq21FXi07kttbJcQwkV4moxcnBTNl3dO4Is7xnPxyCh+Ty/inv9tY+xTy3j8690czil1djM7H58QuPJzOP0RyNkHb0+BHZ932MebjCZemvISsf6xPLvhWVamrmy8mePGEXzrLVTs2kX2Cy822N6zbwBJM2LJTS1h7ZeH2rvZHa45PcKHsa/XuVFX76yUcgfO01p33L+oE0iPUIimFZaZWbg5lY/WHeFQtj0EJ8SHcM2pvTljYLisXnO8I6vhf9dDSQaMvhmmP9Fhcw5TilK48psrqbRW8sHZHzAgqGGvVFssHLnmWso3byb6rTfxnTy53nab1cai5zeTmVzErLtHENm/R4e03VHaemp0IHA6MBJ7L3ATsExrvcfRDe1sJAiFODmt7Uu6fbw2he93ZWCxaRLCfbljSjwzh/WSShh1FWfCwhsh+VeIGg2XfgR+4R3y0ZszN3PTDzcR5BnEZzM/I9gruME+5rQ0Ds2+AGU00ufLL3ALDa23vSCzjM+eWI+nr4nLHh6Dh5frrF/rsOkTSikD9kA8AxgALNZaL3FEIzsjCUIhWiarqIJ3Vx3mw7VHKKuyEhfiw22T+zJrRCTubq53XaldWC2w7FF7vUO/XnDZxxB5Sod89JcHv+TBVQ8yKmIU/532X9wMDYOs6NtvSfvrHHynTCHq9dcaTJnYuTKNXz7ZS8KYcKZdP7hD2u0IDrlGWC0Ae6/wPa31tV05BIUQLRfm78kDMwby231T+cvp/cgtqeTvC7cz5fkVfLAmmQqz1dlNdD6jG0x/3D7HsDwP5p7dYdcNz+t7Hpf1v4wNGRv4z+b/NLqP/9ln4z9zJiU//0zhwoUNtg+e2IveQ4PZty6T/Rsz27vJHaI5p0bf11pfp5S6ArgF+BVIAr7TWr/U/k10HukRCtE2xRVmPlx7hHd/PUxuaRVhfh7cMrkvV4yOwctdRpqSthk+vdK+ePfk+2Dy/e1eycJsNXPD9zewNXsrL0x+gemx0xvsYy0s5NB552MrLibuiyW4R0fX215WVMWn/7cOm1Vz6UOj8QvybNc2O0JbrxEu11pPVUr9Apypta6oXuh6ldZ6fDu0t9OQIBTCMcqrrHyyPoW3fjlIVnElIb7u3DypD1eN7Y23u+tcZ2oXRekw/zJI3wqDZ8OsN8DUvkvbZZZmcsnXl1BhqWD+OfMbXZO0dPVqUm64Ea9TTqH3hx80KNmUvD2Hpa9vJzIhkPPuHtHpB0e19dToAKXUf4HadYKqR4+69sJzQogO4+Vu5MYJcaz8+xQeO38wJqOBJ7/Zw4RnfuaNFQcpqbQ4u4nO498Trv8WBp4HuxbDhxdAeUG7fmS4TzjPT36eSmslc1bMocxc1mAfn3Hj6HHN1ZRv3kze3LkNtscOC2HI5EjS9hWw8Zvkdm1ve2tOj7B3nafHtNZmpZQPME5r/WO7ts7JpEcoRPuotFj5fFMqr/98kLSCcgK9Tdw0IY5rxsXi79lNl3Cz2eD7f8C6NyBsMFy10B6S7ejdHe/y0uaXmB0/m8fGP9awSRUVHJ41G/OxY8QtWYxHn/o9R0uVlYXPbSIntYSZdwyn95CGI1E7i3ZddFsp1VNrnd6mg3RSEoRCtK8qi41Fm1N5bcUBjuaV4+/pxo0T+nDd+Njuuaap1vDbS/DToxAQA9csgeC+7fZxNm3j9mW381vabzw18Slm9pnZYJ+yzVs4cuWVeA0fTu+PP2pwirQwu5z/PbUBq1Vz7p+H0ys+sN3a2xbtHYR5wByt9fvVzy/XWs9v00E7CQlCITqG2Wpj8ZY0Xvv5AEdyy/DzcOP68bHcMCGOQO9ueBVmy0fw5Z/BOwSuXgwRjRfZdYS8ijwu+vIiSs2lfH7u50T7RzfYJ/Opp8ib9wFh991H8PXXNdiecaiQL/+zFZtN03tIMD37BhA1IIiQKN92a3dLtXcQ/gx8AARrrZ+vGVzTpoN2EhKEQnQsi9XGl9uO8eryAxzKKcXPw43bp8Rz/fjY7ree6e9fwec32AfOXLUIohxee7zW2vS13PzDzQwOHswHZ3+A6bg6irbycvsp0owM4hYvanCKFOxhuHrhAdIPFYIGd08jf3ppcoP9nMWR8wgbpbWeC+xVSj0FdO6hQ0KITsvNaOCCU6L4cc5kXro0kUAfE898t4epz69g4aZUbLZutMTxwHPhigVgNcOHs+Hohnb7qLE9x3L9kOvZmbuT17e93mC7wcuLnk89ia6q4tgDD6CtDeeDRvQJ4IK/jeTG5yYS2T8Qc6XrzBltVRAqpSbWeXoVgNb6K2ApkOCAdgkhujGjQTFrRCQ/zZnMwzMHUWa2cs//tjHzlVWs2p/j7OZ1nL5T4Mr/gc1SHYbr2+2j7ky8k0HBg3h3x7tszNjYYLv3KacQdO21VGzb3ugo0hqeviY8vF3r+m5re4T/VkoFA2it0wCUUhdprVdprSMd1johRLfm4WafdvHL36Zwy+Q+HMgu4ap313H93PUcyCpxdvM6RuwEexhqq31qRTuFoclo4umJT+Pp5sk/Vv2D4qriBvuE3n0X7nFxZP/nZSr372+XdjhDa4PwPuDDmidKqb8A/3BIi4QQ4jgBXiYeOHsgy++ZzPmJvfh5bzZnvbSSR7/cRUFZlbOb1/5iJ9jLObVzGMYFxHFv0r2kl6bz7IZnG2w3eHrS6+mn0FYraffdh67qGt/7VgWh1noZsEUp9ZBS6mngfOA0RzZMCCGOF9XDm/9cNoJFt49jSGQA769OZvJzK3j/t8OYrTZnN699xY6vH4Zpm9rlYy5OuJjxvcaz5MASVhxd0WC71/DhBN9yM5W7fyf7jTfapQ0drUVBqJTKUkp9rZR6BFgLXAn0wr70WlF7NFAIIY53SkwPFt02jpcuTcTLZOTRr3Zzzsu/svZQrrOb1r5qwtBmsYdhxk6Hf4RSin+N+xd+7n78c/U/yS1v+D0Nve02PAYNJPet/1K+davD29DRWtojHAC8DFiAGwBfYDLwmVLqQQe3TQghmmSoHlCz/N7J/OX0fiTnlnHZf9dy16dbyCyqcHbz2k/seLjsIzCXwYezIOeAwz8i3Cech8c+TF5FHo+ueZTjp9kpd3cin3kG5eZG2n33YSstdXgbOlKLglBrnae1/kFr/aTWerbWOhoYA8xF1h4VQjiBt7sbc6Yl8ONfJzF1QBhfbD3G1OdX8PbKQ133dGn8GXDRe1CWZx9NWnTM4R9xdtzZnNPnHFYcXcHiA4sbbPfo14+we+ZgPpJCxpNPOvzzO1Jrp0/4KaVmKKVuBqYBKVrrfzq2aUII0Xy9g31477pRvHNNEkG+7jzxze+c8/KvrD+c5+ymtY+B58J5r0Bhiv00aZnjv85/jPkHYd5hPLfhOTJKMxps73H11fhMmEDhwkUULl3q8M/vKC29RthbKTUPOAjcg/206CzgU6XUHqXUDY5vohBCNN8Zg8L58a/Vp0tzyrjkrTXcs2AbOSWVzm6a4424Eqb9H2T/bi/lZC536OH93f3556n/pMRc0vgpUoOBXk8/hTEkhPSHH6Fi7z6Hfn5HaWmP8ENgIdBLa3261vpKrfWFWutBwFlAP6XUnQ5vpRBCtICnycicaQl8/9dJTOwXwsLNqUx9fgXz16c0+GXu8sb/Bcb9GY6ug4U3gc2xK7pMiprE+X3P57e03/ji4BcNtruFhBD1n5fQZjOpt99O+Y4dDv38jtDmtUa7MllrVAjXp7Xmmx0ZPPb1LjKLKpnYL4RnLhxGr8D2LX7boWw2WHQT7FwIo/4EM54D5bjVLgsrC5n9xWwqrBUsOX8JYd5hDfYpWLiI9EceAZuNXafOIcs9jtvfPN1hbWirdl1rVCnVvgWzhBCiDZRSnDOsJz/OmcwlSVH8uj+HM/+9kgUbj3ad3qHBYK9sHzsRNrwNq19x6OEDPAJ45NRHKK4q5rE1jzX6fQu88AL6LFmM3xmnYyspdanvrSMW3d6llLqu5olS6nIHHFMIIRzK39PEsxcN573rkvByN/L3z7dz47yNXWeqhZsHXPoRhA6AHx+GnYscevjTok9jZp+Z/JL6C18f+rrRfTz69SPqlVdwC+68BXob44gg3AYopdS91c//5IBjCiFEu5g6wD6YZvaISJbvyWLai7+wcFOqS/VgmuQVaJ9w7xsBi291eMWK+0ffT7BnME+tf4rssmyHHtuZpAyTEKLbCfA28e9LE3n7miQ8TEbu+d82buoqvcPAaLjiM1AG+PQKKDjqsEPXPUX6xLonusYfD0gZJiFENzZtUDg//nUSs0dEsqy6d/jVNsdPTu9wvRLhgregNMs+raLScZU6psZM5azYs1iWsozvj3zvsOM6k5RhEkJ0a4He7rW9Q3c3A3+ev4U5C7ZSXGF2dtPaZtD5MPUhyNwJi2+xjyx1kPtH30+gRyBPrXuK/Ip8hx3XWaQMkxBCYO8dfnvXJKb0D2XR5jRmvPwrm464+C/5iffCkItgz9fw8xMOO2ywVzD3j76fvIo8nlr/lMOO6yxShkkIIaqF+nnw3nWjeOz8wWQVVXLJW2t47ecD2Gwuei1MKTj/Veh1Cvz6POxuOCG+tWbEzeC06NP49vC3/JD8g8OO6wxShkkIIepQSnHNqbF89ecJxIf68tz3e7nmvfVkFbvoQBqTl31ahU8oLLkdsvY45LBKKf556j/p4dGDx9Y+RmZppkOO6wxShkkIIRqREO7HF3eO5/LRMaw6kMOM//zKr/tddMpAQCRc/L59LdLProQKx/RbQrxCeHTcoxRWFvK3lX+jwuKafyxIGSYhhGiCp8nIUxcM5dUrRlBptnHNe+t5Zdl+1zxVGjsBpv8f5B6AL+8EB019mBozlWsHXcuWrC3c+8u97MrZRYXVtQKxzfMItdYZWuuvpQyTEKKrmjmsF1/9eQL9w/144cd93DhvAwVlVc5uVsuNvd0+mnT3F7D2dYcddk7SHGbEzeCX1F+4bOllFFW61pWyll4jHHOS7V5KqUFta5IQQnQ+sSE+LL59PBeeEsXPe7OZ+coqdqYVOrtZLaMUnPcqBMfDj49AyjqHHNagDDw18Snenv42fxr6J4wGo0OO21Fa2iO8Ryn1o1LqeqXUIKVUsFIqUik1VSn1LLAGCG+HdgohhNN5uRt5/uJhPHXBULKKKrnwjdUs3ORiFWo8/eGSD8Bggs+vh9JchxzWoAyM7TmWv5zyF9yUm0OO2VFaeo3wEuzzBccDXwJpwHbgYSAFGK+1/tnRjRRCiM5CKcXlo2NYcOupBPm4c8//tvHIFzupsjhuwnq7Cx8M57wARWmw+GaHTrZ3RS2+Rqi13qC1vklrHa+19gQu1VpP0Vq/qrUubYc2CiFEp5MYHchXf57A2D5BfLDmCFe9s47ckkpnN6v5RlwJiVfBgZ/gt5ec3RqncsSi2zOVUvcAKKVGKKXmO+CYQgjR6YX4evDRjWO4YXwc65PzOO/V3/g93YUGisx4DkIHwvLH4cgaZ7fGaRwxavRuoK9Sag9wK/BQW48phBCuws1o4JFzB/HshcPIKq7gwjdW893ODGc3q3ncve3zC43usOhmqHCxwT8O4ogK9T8AOcC9QE+gvK3HdLTq0azfKKV+UUrNdXZ7hBBdzyWjopn/p7F4uxu57eNN/HflQdcoUxQ2AM58AgpT4Ju/O7s1TuGIU6MXaK0f0Vp/DdwGvOGAYzradGCV1noyYFZKDXF2g4QQXU9SbBCLbx9PfKgvT36zhweX7MRidYGBKEk3QMLZsP1T2LnQ2a3pcC2dR7iszuMQAK11baGr6pJMFzmsdY5zGPCpfuwLdM/+vxCi3UUHefP5beOYEB/CJ+tSuGHeRkorLc5u1okpBee9Yl+P9Os5UNQFajK2QEt7hEF1Htdbbrx6IW601g4p4qWUelkplayU0kqpxOO29VNKrVZK7VNKbVBKDT7J4Q4AE5VSv9ubqB1XslkIIY4T4GVi7vWjuCQpipX7srninXXklXbylWh8Q+Hcl6GiAJbc1q2mVLTl1Kg67vmsNhyrMZ8DE4AjjWx7C/iv1joBeAZ4H6B6kv+K4273A9cCn2mtBwI5SqlxDm6rEELUYzIaeObCYdw5JZ5tRwu46M3VpOaXObtZJzZgBpxyDRxaARvecXZrOkxLg1A38RgaBmObaK1Xaq0bLNmglAoDkoCPql9aCEQrpeK11ru11qcdd3sa+9eZV71/HhDY2GcqpeYopVJrbiUlJY3tJoQQzaKU4t4z+/PouYM4nFPKRW+s4WB2J/+9cuaTENjbvgRbzn5nt6ZDtDQIE5RSHyul7gZ8lFKedbZ11PCoaCBda20B+3lO7KvaxJzgPZ8A1yulVgCJwI+N7aS1flFrHVVz8/X1dWjDhRDd03Xj43jp0kSySyq59K01nXuuoYcfzH4TLBWw+FawdvLrmw7Q0iA8A1gNDMU+TaJAKXVAKbUEiHJw2xxGa52vtZ5e3UOc7ajrmEII0VznJ0byxpWnUFRu4bL/rmXb0QJnN6lpvcfBuDshbSOsedXZrWl3LV1rdK3W+jWt9Y1a60QgALgc+BZY4vjmNeoo0FMp+6quSimFvTeY0kGfL4QQrTJ9cARvX5tEpcXKVe+sY9ORvJO/yVmmPAjB/eDnJyF7r7Nb067aNI9Qa11ZvfboW1rrmx3VqJN8ZhawGbiq+qULgVSt9YGO+HwhhGiLyQmhvH/9aKxac/W761l7yDHVHxzO5AWz3gBrFXxxJ9iszm5Ru3HEhPp2oZR6SymViv2U6/dKqbpBdwtwi1JqH3A/cL0z2iiEEK0xtk8wH944GqNSXDd3Pav25zi7SY2LHgVjb4PU9V16FGmnDUKt9S3Vg1bctNbhWuv4Otv2aq1P1VonaK2TtNY7nNlWIYRoqZG9g/jopjG4Gw3cOG8DK/dlO7tJjZv6EATGwE//goKueQWq0wahEEJ0dcOjA/nkT2PxNBm56YONnTMM3X1g5ktgLoWl94IrrJ/aQhKEQgjhREMiA/j4pjF4VYfhbwc64WnS+NNh2KWw/3vYtcjZrXE4N2c3wFXZbDbXWFle1FJKYTDI336i86kJwyveXsufPtjIhzeOYWTvHs5uVn1nPmUv4vvtfdBnCngHnfw9LkKCsIWqqqpISUnBbJapiK7IZDIRExODu7u7s5siRD1DIgOYe/1orn53HdfNXc/8P41lSGSAs5v1B59gOOtpWPQn+OEhmPW6s1vkMBKELZSSkoKfnx/BwcHYpzAKV6G1Jjc3l5SUFOLj40/+BiE62MjePXj7miSuf38D1763ns9vG0dciM/J39hRhl4M2z+DrR/D0Iug71Rnt8ghJAhbwGazYTabCQ4Oxs1NvnWuKDg4mLy8PGw2m5wmFZ3S+PgQXr18BLd+tIlr3lvHwtvGEebnefI3dgSlYOa/4bUx8PVf4bY19ir3Lk5+E7RAzTVB6Qm6rpp/O7m+Kzqz6YMjeHL2UI7mlXPtexsoruhEl2ICY+xTKvKTYcVTzm6NQ0gQdgFKKQoKCuq9Fhsby9atWx32GY8//jh9+/alb9++PPjgg03ul5WVxVlnnUW/fv0YMmQIK1eurN325JNP0r9/fwwGA0uWLHFY2wCGDRtGenq6Q48phDNdNjqGe6Yl8Ht6Ebd/vBlzZ6p0P+ZW6DUC1rwG6duc3Zo2kyAUJ7Vy5Urmz5/P9u3b2b17N99//z1Lly5tdN/777+fsWPHsn//fubOncsVV1xRO7DojDPO4Ntvv2XSpEkObd+BAwfw9fWlZ8+eDj2uEM5259R4LhsVza/7c3jki52d50yGwWgv4gvw5V9cvkKFBKE4qc8++4yrr74aHx8fPDw8uOGGG5g/f36j+y5YsIBbb70VgFGjRtGrVy9++eUXAEaPHk2fPn1O+nnvv/8+Z5xxBpdffjmDBg1i3Lhx7N69m9mzZzNw4ECmT59O3VqRixcvZtasWYC95zpw4EASExNJTEzkyJHG6joL4RqUUvzfrCFM7BfC/PVHefOXQ85u0h96DoNT74D0rbDhbWe3pk1kxEcb3TRvA0dy26fqdO9gb965dpRDj7l3714uvfTSRreNGDGCuXPnNng9JSWFCRMm1D6PjY3l008/bbBfbm4uZrOZiIiIevumpLR8WaYNGzawY8cOYmJiuPrqqzn33HNZvXo14eHhzJw5k3nz5nHHHXcAsGTJEubOnUt+fj7PP/886enpeHl5UVZWJgNihMszGQ28fuUpXPzmGp79fg99Q32YPjji5G/sCKfdD7uXwPLHYeC5ENBpq/GdkARhF9bYoJ7+/fs79Nphezn11FOJibHXWk5KSsJsNhMeHg7Ye5r799srZ2dkZFBUVERCQgJWq5V+/fpx1VVXMX36dM455xyiolzzP6YQdfl5mnj7miTOf+037v5sKwtvG8fAnv7ObpZ9+bVzXoSPL4KFN8EVn4GHP9g60eCeZpAgbCNH99haIzQ0lNzcXAIDA2tfy8nJISwsrMG+rekRxsTE1DvFmJycXBtSddVMK8nIyKjtFTa178l4ev4xXNxoNDZ4brHYr0l88cUXnHfeebWvr127ltWrV7NixQrGjh3L/PnzmThxYos/X4jOJjrImzevGsmV76zlpnkb+fLO8QT7eji7WdBvGoy9A9a+Bu+eCbETwOZa1+vlvFEXcOaZZ/LWW2/VPv/ggw/o06dPo4NHanqEjd0aC0GAiy++mA8//JDS0lIqKyt57733uOyyy5rc98033wTspzfT0tKYPHmyA77Kxi1ZsoTZs2cDUFxcTGZmJhMnTuThhx9mwoQJbNmypd0+W4iONjouiMdnDSGtoJw7P9nSeUaSnvkEnP5PyN3vktcLJQi7gJdeeon09HSGDRtGYmIin3zyCf/73/8cdvzTTjuNSy+9lKFDhzJw4ECmTZvGzJkzAdi4cSMzZsyo3feZZ55h9erV9OvXj+uuu46PPvoIk8kE2AeyREVFsWbNGm666SaioqLIzm79avtFRUXs37+fpKQkAAoLC7ngggsYOnQow4YNw2w2c+2117bhKxei87l0VAzXnNqbNYdyefKb353dHDulYOIcuGevvVKFi821Vp1mOG4nFBUVpVNTU2ufW61W9u3bR0JCAkaj0YktEwCffvopq1at4tVXX232e+TfUHQFZquNK99Zx/rDebx4yXAuOKVzXQufe+2/KfMcyh1vneHsptRSSqVprRv9RkmPULisyy67rEUhKERXUTOSNMLfkwcX72RvRrGzm+TSJAiFEMIFhfh68OoVIzBbbdz28SZKKl17UrszSRAKIYSLSooN4v6zB3Aou5QHFu3oPCvPuBgJQiGEcGE3Tohj+qBwvtp2jPnrjzq7OS5JglAIIVyYUornLhpOZKAX//pqF7+nFzm7SS5HglAIIVxcgLeJV64YgdWmufOTzZRXWZ3dJJciQSiEEF3AKTE9uGd6fw5ml/LEN7ud3RyXIkHYBbhKPcLrr7+ehIQEhg8fzvjx49mwYYPD2jdz5kzWrVvnsOMJ4YpumdSHU/sE89HaFH7anens5rgMCUJxUo6qRzh79mx2797Ntm3beOCBB7j44osd0r6SkhJ+//13Ro8e7ZDjCeGqDAbFC5cMJ8DLxH0Lt5NTUunsJrkECUJxUo6qR3jeeefh5mZf533s2LGkpaXVLp5d14oVKxgyZAi33XYbw4YNY+jQoWzfvp3rrruOoUOHMmbMGNLS0mr3//bbbznzzDNRSvHOO+8waNAgEhMTGTp0qPQSRbfTK9CLx2cNIbe0igcXy5SK5pDqE231yWWQf7h9jt0jDq5oWPevLTpLPcL//Oc/zJgxozYYj7dnzx7mzZvHG2+8wcMPP8zUqVNZtWoVAwYM4I477uCll17iueeeA+yFea+77joA7rnnHvbs2UPPnj0xm81UVspfxKL7OXd4L77blcHS7eks3pLW6ZZg62wkCLuwzlqP8KOPPmLBggX1rh8eLz4+npEjRwL2eoTx8fEMGDAAsFe6X7x4MQBms5nVq1czb948AE4//fTaQr5nn302CQkJ7fzVCNE5PX7+ENYdyuOfX+5ifHwI4f6eJ39TNyVB2FYO7rG1hivVI/zss8/417/+xbJly2oL7TamufUIly9fzvjx42srXCxcuJBNmzaxYsUKZsyYweOPP95kySghurIePu48dcFQ/vTBRh5cvIO3r0lq9I9jIdcIuwRXqUe4YMECHnroIX766adWFettTN16hBaLhYMHD5KUlMS9997LRRddxPr16x3yOUK4ommDwjlveC9++j2LL7cdc3ZzOi0Jwi7AVeoRXnnllVRUVHD++eeTmJhIYmIiubm5rW6X1prvv/+es846C7CXWLrhhhsYMmQIiYmJbNq0iTlz5rThKxfC9T163mCCfdx59Mtd5Moo0kZJPcITkHqEndvatWt5/PHH+frrr5v9Hvk3FN3RV9uO8ef5W5iV2IuXLhvR7p8n9QiF6CBjx45tUQgK0V3NHNaTqQPCWLL1GL/sy3Z2czodCUIhhOjilFL836wh+LgbeXDxDlmL9DgShEII0Q1EBnpxz/T+pOaX88ry/c5uTqciQSiEEN3ENaf2ZnAvf/678hD7Moud3ZxOQ4JQCCG6CTejgSdmD8WqNQ8t2SnLr1WTIBRCiG4kMTqQK8fEsP5wHku2pp38Dd2ABGEX4CplmGosX74co9HISy+95LD2SRkmIZrv3un9CfJx54mleyiqMDu7OU4nQShOylFlmAAKCwu5//77603CbyspwyREywR6u3P/WQPIKank3z/uc3ZznE6CUJyUo8owAdx555089NBDBAcHN/l5UoZJiPZ30cgohkcH8uGaIxzIKnF2c5xKFt1uoz8v+zNHi4+2y7Gj/aJ55fRXHHpMZ5Zh+vzzzzEYDJx33nksWrTohO2UMkxCtC+DQfHPcwdxweureWLpbuZe333PqEgQdmGdqQxTRkYGjz/+OCtWrGjW/lKGSYj2d0pMD85P7MUXW4+xYm8Wp/VvWLGmO5AgbCNH99hawxXKMG3atIn09HQSExNr2/fll1+SnZ3NE0880eBYUoZJiI5x31kD+H5XBk9+8zsT4kNwM3a/K2YShF1ATRmmZ599FmheGaaWuPjii7njjjv485//jJubG++99x6PPvpok/u++eabPProo/XKMJlMJjIzM2v3u+6660hMTOTuu+9uUVuOd3wZpuTkZJKSkkhKSiInJ4f169dLEApxAr0CvfjTxD68svwACzamcsUYx5RIcyVdLvqVUkFKqU1KqZLjXn9FKfWrUuofzmpbe3GVMkyOJmWYhHCMWyb3JcTXgxd/3EtJpcXZzelwXa4Mk1LKBPgBC7TWZ1S/lgRcqbX+q1Lqc+AOrXXmiY4DUoaps5MyTEI4zvz1KTywaAd/mRrPnOn923QsKcPkZFprs9Y677iXxwDLqx//Aozs2FaJ9iBlmIRwnItHRhEf5svbvx4mq7jC2c3pUJ0iCJVSLyulkpVSWimVeNy2fkqp1UqpfUqpDUqpwa34iECgqPpxcfVzIYQQ1dyMBv5+Zn/KzVZeXta9qlN0iiAEPgcmAEca2fYW8F+tdQLwDPA+gFJqkFJqxXG3+5s4fgHgX/3Yr/q5EEKIOqYNCmdk7x7MX3+U5JxSZzenw3SKINRar9Rapx7/ulIqDEgCPqp+aSEQrZSK11rv1lqfdtzt6SY+Yj0wpfrxJGBTYzsppeYopVJrbiUl3Xu1BSFE96KU4u9n9sdq07zYjZZe6xRBeALRQLrW2gKg7SN7UoATju9VSv0EjFBK/aSUGqK13gB4KKV+BbY1NVBGa/2i1jqq5ubr6+vYr0YIITq5MX2CmZwQypfbjrH7WNHJ39AFdPYgbBWt9Rla6+Dq+53Vr92htZ6otX7c2e0TQojO7G9n2keNvvDDXie3pGN09iA8CvRUSrkBKPuaYTHYe4VCCCHawZDIAM4Z2pNle7LYkpLv7Oa0u04dhFrrLGAzcFX1SxcCqVrrA85rVefjKvUItdY8+uijJCQkMHToUKZMmdLkcVpq2LBhpKenO+x4QnR3d5/RD6Xg3z91/RGknWKJNaXUW8A5QATwvVKqWGsdX735FuD96hVhioDrndTMbqtuPUI3NzfGjx/PuHHjOOeccxrsW1OP8LvvvmPDhg3Mnj2bw4cPYzKZePnll9m+fTs7d+7E3d2djIwMh7TvwIED+Pr6NrqknBCidfqF+3H+8F4s2XqMjcl5JMUGObtJ7aZT9Ai11rdUD1Bx01qH1wlBtNZ7tdanaq0TtNZJWusdzmxrd+SoeoTPPfccTz/9NO7u7gD1yjXV9f7773PGGWdw+eWXM2jQIMaNG8fu3buZPXs2AwcOZPr06dQd0bt48WJmzZoF2HuuAwcOJDExkcTExHqLhQshWuauMxIwKHipi/cKO0WP0JUdve12qo62zyVL9+gYot943aHHdFY9wqKiIjIzM/niiy/4/PPPAZgzZ06TbdmwYQM7duwgJiamtqzS6tWrCQ8PZ+bMmcybN4877rgDsC+8PXfuXPLz83n++edJT0/Hy8uLsrIyDIZO8beeEC4pLsSHWSMiWbQ5jQ3JeYzqor1CCcIurDPVI7RYLFgsFsrLy1m3bh3JycmMGzeOAQMGMHz48Ab7n3rqqbWlnpKSkjCbzYSHhwP2nub+/fa/UDMyMigqKiIhIQGr1Uq/fv246qqrmD59Oueccw5RUY0uLSiEaKY/T+3Hki1p/PvHfXzyp7HObk67kCBsI0f32FrDFeoRBgUF4evry1VX2cc9xcbGMn78eDZs2NBoEDa3HuEXX3zBeeedV/v62rVrWb16NStWrGDs2LHMnz+fiRMnNvr1CiFOLi7Eh1mJkSza0nV7hXLeqAuoqUdYozn1CBu7NRaCYK8x+OGHH1JaWkplZSXvvfdekzX+auoRAvXqEQJcfvnlfPfddwDk5eWxfv16hg0b1qavvW49wuLiYjIzM5k4cSIPP/wwEyZMYMuWLW06vhAC7pgaj0HRZdcglSDsAlylHuFTTz3Fd999x5AhQ5g0aRL33Xcfo0ePbnW7ioqK2L9/P0lJSQAUFhZywQUXMHToUIYNG4bZbObaa69tw1cuhADoG+rLzGG9+HV/TpecV9jl6hE6ktQj7Nw+/fRTVq1axauvvtrs98i/oRCtsz+zmOkvrWRq/zDevW7UCfeVeoRCdJDLLrusRSEohGi9fuF+nDkogmV7stiT0bXWIJUgFEII0Sy3ndYXgLd+OeTkljiWBKEQQohmGR4dyPj4YL7cdoyjeWXObo7DSBAKIYRottsmx2O1ad75tev0CiUIhRBCNNv4+GCGRPrz2caj5JZUOrs5DiFBKIQQotmUUtw6uS8VZhvz1nSNtXwlCIUQQrTI2UN60jvYm3mrkymptDi7OW0mQdgFuEo9wvXr1zN27FhGjBjBwIEDefbZZx3WvpkzZ7Ju3TqHHU8I0TSjQXHLpL4UlpuZv87166RLEIqTqluPcPfu3Xz//fcsXbq00X1r6hHu37+fuXPncsUVV2A2mwG4+eab+cc//sGWLVv47bffeP7559m9e3eb21dSUsLvv//eplVqhBAtc+HISML8PHj710NUWqzObk6byKLbbbT09e0UZpe3y7EDQr045/a2rcXpCHXrEQK19QgbK8y7YMECDhw4ANSvR3jGGWfU67mWlpbi7u5OUFDDBXxXrFjBnXfeycSJE/ntt9/QWvPxxx/z4osvsmnTJry9vVm0aBGRkZEAfPvtt5x55pkopXjnnXd48cUXcXd3x2q18s477zBmzJh2+s4I0X15uBm5aWIcT36zhwUbjnL1qbHOblKrSY+wm9m7d29t0drjb9dff32j70lJSaF37961z2tqDB7vRPUIAebOncvDDz9MTEwMCQkJPPnkk00W592zZw833XQT27dvZ9asWUydOpX777+fHTt2kJSUxEsvvVS7b93CvPfccw/Lli1j69atbN68mcGDB7f0WySEaKarxvYm2MedV38+QIXZdXuF0iNso87QY2tKZ6pHCPD000/z1FNPccUVV3Do0CEmT55MUlISgwYNarBvfHw8I0eOBOz1COPj4xkwYAAAo0ePZvHixQCYzWZWr17NvHnzADj99NNrC/meffbZJCQkdNBXJ0T34+3uxm2n9eXxpb/zyboUbpgQB4DNxZawlh5hF1BTj7CuE9UjbGmPsDX1CI/fNycnh8WLF3PFFVcA0KdPH8aOHctvv/3W6Gc2tx7h8uXLGT9+fG2Fi4ULF/L0009jNpuZMWMGn376aaPHF0I4xlVjexPm58FrPx+guMJMfmkV4FpJKEHYBbhCPcIePXrg4+PD8uXLAXtQr1u3jiFDhrTpa69bj9BisXDw4EGSkpK49957ueiii1i/fn2bji+EODFPk5F7pieQW1rFGysO8uW2Y85uUovJqdEu4KWXXuLuu+9m2LBhGAwGIiIi2q0eIcCll15arx7hI488wjfffAPY6xFeffXV9OvXD3d393r1CBcsWMDf/vY3LBYLZrOZu+++m1NPPbXV7dJa8/333/Pcc88B9hJLN9xwA3l5ebi5uREaGtpkuAshHOeikdHM/S2Zd349jL+XG7c5u0EtJPUIT0DqEXZua9eu5fHHH+frr79u9nvk31CI9rEjtZC/fb6NPRnFPHxsLxVew6QeoRDtbezYsS0KQSFE+xkaFcB3d09i5d+mAA0H6nVmEoRCCCEcJibY29lNaDEJQiGEEN2aBGEL1MzLk+uqrqvm366xOZZCiO5JRo22gMFgwGQykZubS3BwsPwydTFaa3JzczGZTBgM8jegEMJOgrCFYmJiSElJIS8vz9lNEa1gMpkaXQxACNF9SRC2kLu7O/Hx8dhsNjlF6mKUUtITFEI0IEHYSvILVQghugb5bS6EEKJbkyAUQgjheMp14kWWWDsBpVQlkO3sdrSQF9A+lYKd//mOOHZbjtHS9zZ3f0ft5wuUNOM4rkh+rtvvGO31c93cfTvq5zpUa+3R6Battdy60A1Y0FU/3xHHbssxWvre5u7vqP2AVGf+27fnTX6u2+8Y7fVz3dx9O8PPtev0XUVzOa7sROf7fEccuy3HaOl7m7u/o/fripz9tcvPdev2b86+zv63lVOjQnQVSqlU3cTq+kK4qo74uZYeoRBdx4vOboAQ7aDdf66lRyiEEKJbkx6hEEKIbk2CUAghRLcmQSiEEKJbkyAUQgjRrUkQCtGFKaWClFKblFJddcUZ0Q0ppSYopdYqpVYrpe5p6/EkCIXo2oqBacBaZzdECAc6BEzSWo8DZiqlvNtyMCnDJEQXprU2A3lKKWc3RQiH0Vofq/PUCtjacjzpEQrhIpRSLyulkpVSWimVeNy2ftWnifYppTYopQY7qZlCtEhbfq6VUtOAg1rrira0QYJQCNfxOTABONLItreA/2qtE4BngPc7sF1CtEWrfq6VUlHAA4BcIxSiu9Bar9Rapx7/ulIqDEgCPqp+aSEQrZSK78j2CdEarfm5Vkp5YA/F27TWbR4IJkEohOuLBtK11hYAbV83MQWIAVBK/QSMUEr9pJQa4rxmCtEiJ/q5vgIYBLyllFqhlIpsywfJYBkhujit9RnOboMQjqS1ngvMddTxpEcohOs7CvRUSrkBKPsQ0Rjsfz0L4ao67OdaglAIF6e1zgI2A1dVv3Qh9qreB5zXKiHapiN/rqUMkxAuQin1FnAOEAHkAsVa6/jqbf2xDx4IBoqA67XWO5zUVCGarTP8XEsQCiGE6Nbk1KgQQohuTYJQCCFEtyZBKIQQoluTIBRCCNGtSRAKIYTo1iQIhRBCdGsShEIIIbo1CUIhhBDdmgShEEKIbk2CUIhuRCkVrJTKdnY7hOhMJAiF6F6SgA3OboQQnYkEoRDdSxKwHkApdbpSapNSarKT2ySEU8mi20J0I0qpJcB7wFjgFOCa6nI3QnRb0iMUontJAp4FKoAZWusspZSfUupHAKXUf5RS05zaQiE6mAShEN2EUqonYAM2Ar5aaxuA1roYWKeUegN74dMfndhMITqcBKEQ3UcSsAm4AUhSSt1VZ9s64HSt9XNOaZkQTuTm7AYIITrMKGCz1rpKKTUb+EUplQ78Btxc/fwUrfVmp7ZSiA4mg2WE6MaUUu7AQuB2wAN4Qmt9qXNbJUTHkiAUQgjRrck1QiGEEN2aBKEQQohuTYJQCCFEtyZBKIQQoluTIBRCCNGtSRAKIYTo1iQIhRBCdGsShEIIIbo1CUIhhBDdmgShEEKIbu3/AQ3jF0Xmn0qUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 480x320 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "for path in path_runs:\n",
    "    sim = fls.load(path, hide_stdout=True)\n",
    "    period = sim.forcing.get_info()[\"period\"]\n",
    "    data = sim.output.spectra.load1d_mean(tmin=period)\n",
    "    spectrum = data[\"spectra_vx_kx\"]\n",
    "    kx = data[\"kx\"]\n",
    "    ax.loglog(\n",
    "        kx,\n",
    "        spectrum * kx ** (5 / 3),\n",
    "        label=f\"U = {sim.params.forcing.milestone.movement.periodic_uniform.speed} m/s\",\n",
    "    )\n",
    "\n",
    "ax.set_xlabel(\"$k_x$\")\n",
    "ax.set_ylabel(\"$E_x(k_z) {k_z}^{5/3}$\")\n",
    "plt.legend();"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}