PW.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3c90ac17",
   "metadata": {},
   "source": [
    "# All-Optical Radiation Reaction at 10^21 W/cm^2\n",
    "\n",
    "https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.113.134801\n",
    "\n",
    "Setup: Classical Radiation Reaction in a Plane Wave. Check validity of equation 3 of paper with OSIRIS-QED 1D simulations\n",
    "\n",
    "Initial electron energy $\\gamma_0 = 2000$\n",
    "\n",
    "- PW different $a_0$, frontal collision\n",
    "- PW $a_0$=20, varying angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8d933037",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import h5py\n",
    "import glob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "54e996ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "g0 = 2000 #[] initial electron energy\n",
    "\n",
    "def gf(g0,a0,th):\n",
    "    \"\"\" equation 3 \"\"\"\n",
    "    k = 4.4199e-5 * (a0/10)**2 * (1+np.cos(th*np.pi/180) ) / 2\n",
    "    return g0/(1+k*g0)\n",
    "\n",
    "# parameter scan: a0 (angle=0)\n",
    "files = glob.glob('PW/a*103.h5')\n",
    "files.sort()\n",
    "gf_a0_sim = np.zeros(len(files))\n",
    "a0lst_sim = np.zeros_like(gf_a0_sim)\n",
    "for i in range(len(files)):\n",
    "    a0 = int(int(files[i][4:6]))\n",
    "    a0lst_sim[i] = a0\n",
    "    filef = h5py.File(files[i], 'r')\n",
    "    enelst1 = np.array(filef['ene'])\n",
    "    gf_a0_sim[i] = np.mean(enelst1)\n",
    "a0lst_th = np.linspace(np.min(a0lst_sim),np.max(a0lst_sim),200)\n",
    "gf_a0_th = gf(g0,a0lst_th,0)\n",
    "    \n",
    "# parameter scan: angle (a0=20)\n",
    "files = glob.glob('PW/th*103.h5')\n",
    "files.sort()\n",
    "thlst_sim = np.array([0,15,30,45,60,75,90])\n",
    "gf_th_sim = np.zeros(len(thlst_sim))\n",
    "for i in range(len(files)):\n",
    "    # sim\n",
    "    filef = h5py.File(files[i], 'r')\n",
    "    enelst = np.array(filef['ene'])\n",
    "    qlst = np.array(filef['q'])\n",
    "    gf_th_sim[i] = np.sum(qlst*enelst)/np.sum(qlst)\n",
    "thlst_th = np.linspace(0,90,200)\n",
    "gf_th_th = gf(g0,20,thlst_th)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "97cd11b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA94AAAF1CAYAAAAeKTIvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB2rUlEQVR4nO3dd1wT9/8H8FfCFNkCAoooDhAERHDgtiKI1lG3oLitFrWKrdS9vharFWtbq7VuxVGtW6viAgfKEAQXFaVOhlUBUdn3+6M/06aiohKOJK/n43EPyN0nyevucvnknbvcSQRBEEBERERERERECiEVOwARERERERGRKmPhTURERERERKRALLyJiIiIiIiIFIiFNxEREREREZECsfAmIiIiIiIiUiAW3kREREREREQKxMKbiIiIiIiISIFYeBMREREREREpEAtvIiIiIiIiIgVi4U1ERERERESkQKIX3iEhIWjatCkMDAxgYWGBnj17Ijk5Wa5NXl4eAgMDUa1aNejr66N3797IyMiQa3Pnzh107doVenp6sLCwwJdffomioiK5NqdOnUKTJk2go6ODevXqYf369YqePSIiIpXA/pqIiOj9iV54R0REIDAwEOfPn0d4eDgKCwvh7e2NZ8+eydpMmjQJ+/fvx44dOxAREYEHDx6gV69esunFxcXo2rUrCgoKcO7cOWzYsAHr16/HrFmzZG1SU1PRtWtXdOjQAQkJCZg4cSJGjhyJI0eOVOj8EhERKSP210RERO9PIgiCIHaIf3v48CEsLCwQERGBtm3bIjs7G+bm5tiyZQv69OkDALh+/ToaNmyIqKgotGjRAr///js+/vhjPHjwANWrVwcArFy5EsHBwXj48CG0tbURHByMgwcP4vLly7LnGjBgALKysnD48GFR5pWIiEhZsb8mIiIqO02xA/xXdnY2AMDU1BQAEBcXh8LCQnh5ecnaODg4oFatWrKOPCoqCs7OzrJOHAB8fHwwduxYXLlyBW5uboiKipJ7jJdtJk6cWGqO/Px85Ofny26XlJTg8ePHqFatGiQSSXnNLhERqRlBEPD06VNYW1tDKhX9wLP3Vln6a4B9NhERlb/y7q8rVeFdUlKCiRMnolWrVmjUqBEAID09Hdra2jA2NpZrW716daSnp8va/LsTfzn95bQ3tcnJycGLFy9QpUoVuWkhISGYO3duuc0bERHRv929exc1a9YUO8Z7qUz9NcA+m4iIFKe8+utKVXgHBgbi8uXLOHPmjNhRMHXqVAQFBcluZ2dno1atWrh79y4MDQ1FTEZERMosJycHNjY2MDAwEDvKe6tM/TXAPpuIiMpfeffXlabwHjduHA4cOIDIyEi5bxQsLS1RUFCArKwsuW/RMzIyYGlpKWsTHR0t93gvz6L67zb/PbNqRkYGDA0NS/32XEdHBzo6Oq+MNzQ0ZCdOREQfTFkPga5s/TXAPpuIiBSnvPpr0X9cJggCxo0bh927d+PEiROoU6eO3HR3d3doaWnh+PHjsnHJycm4c+cOPD09AQCenp5ISkpCZmamrE14eDgMDQ3h6Ogoa/Pvx3jZ5uVjEBER0euxvyYiInp/op/V/LPPPsOWLVuwd+9e2Nvby8YbGRnJvtkeO3YsDh06hPXr18PQ0BDjx48HAJw7dw7A35cnady4MaytrbFo0SKkp6dj8ODBGDlyJL7++msAf1+epFGjRggMDMTw4cNx4sQJTJgwAQcPHoSPj89bc+bk5MDIyAjZ2dn89pyIiN6bsvYnytJfA8q7jImIqPIo975EEBmAUod169bJ2rx48UL47LPPBBMTE0FPT0/45JNPhLS0NLnH+fPPPwVfX1+hSpUqgpmZmTB58mShsLBQrs3JkyeFxo0bC9ra2oKdnZ3cc7xNdna2AEDIzs7+kNklIiI1p6z9ibL014KgvMuYiIgqj/LuS0Tf460s+O05ESmTkpISFBQUiB1DLWlpaUFDQ+O109mfKF5Zl3FxcTEKCwsrMBkBb99GiIgqg/LuryvNydWIiKh8FBQUIDU1FSUlJWJHUVvGxsawtLRU2hOoqTpBEJCeno6srCyxo6gtbiNEpG5YeL+jv54UgTsoiKiyEgQBaWlp0NDQgI2NDaRS0c+hqVYEQcDz589lJw+zsrISORGV5mXRbWFhAT09PRZ/FYjbCBGpKxbe72jY/Af4YogGurUxEjsKEdErioqK8Pz5c1hbW0NPT0/sOGrp5YnGMjMzYWFhwUNqK5ni4mJZ0V2tWjWx46glbiNEpI64K+QdCZAgNOwxPp88B5cvXxY7DhGRnOLiYgCAtra2yEnU28svPfj74crn5TrhF1Pi4jZCROqGhfd7kEg1sHHbETg7O8PDwwM//PADHj9+LHYsIiIZHjorLi7/yo/rSFxc/kSkblh4vwcJBLRtYQ9NTU3ExcVhwoQJqFGjBoYOHYqoqCjwRPFERERERET0EgvvdySVApP9q2Hvb+uRlpaG77//Hq6ursjLy8OGDRvQsmVLNG7cGCtWrEBOTo7YcYmIVMKpU6cgkUh4Fmqi1+A2QkRUubHwfkdrp1uhSyt9AICZmRnGjx+P+Ph4REVFYejQodDV1UViYiI+++wzWFtbIzAwEMnJySKnJiKqvCQSyRuHOXPmiB2RSFTcRoiIlB8L73dkZvLqieAlEglatGiBdevW4cGDB/juu+/g4OCAZ8+e4aeffoKDgwO6du2K8PBwHoZORPQfaWlpsuG7776DoaGh3LgvvvhCtGwFBQWiPTfRS9xGiIiUHwvvcmZiYoLPP/8cV69exbFjx9CtWzdIJBIcOnQI3t7ecHZ2xi+//IIXL16IHZWIqFKwtLSUDUZGRpBIJHLj9PX1ZW3j4uLg4eEBPT09tGzZ8pUjivbu3YsmTZpAV1cXdnZ2mDt3LoqKimTT79y5gx49ekBfXx+Ghobo168fMjIyZNPnzJmDxo0bY/Xq1ahTpw50dXWxceNGVKtWDfn5+XLP1bNnTwwePFhBS4XoH9xGiIiUHwtvBZFIJOjYsSP27duH5ORkjB8/HlWrVsWVK1cwevRo1K5dGyEhIfwtFhEplCAIePbsmSiDIo7wmT59OpYsWYLY2Fhoampi+PDhsmmnT59GQECA7MvPn3/+GevXr8eCBQsAACUlJejRowceP36MiIgIhIeH49atW+jfv7/cc6SkpOC3337Drl27kJCQgL59+6K4uBj79u2TtcnMzMTBgwflnp+UE7cRbiNERBVCoDLJzs4WAAjZ2dnv/RhPnjwRlixZItSqVUsAIAAQDAwMhClTpggPHjwox7REpK5evHghXL16VXjx4oUgCIKQm5sre7+p6CE3N/ed869bt04wMjJ6ZfzJkycFAMKxY8dk4w4ePCgAkM1rx44dha+//lrufps2bRKsrKwEQRCEo0ePChoaGsKdO3dk069cuSIAEKKjowVBEITZs2cLWlpaQmZmptzjjB07VvD19ZXdXrJkiWBnZyeUlJSUOh//XQ//Vh79Cb3Zm5Yxt5HKv40QEVUG5d1fc493BTI2NkZQUBBSUlKwadMmODk54enTp1i0aBFq166NTz/9FDdv3hQ7JhFRpeXi4iL738rKCsDfe9YA4NKlS5g3bx709fVlw6hRo5CWlobnz5/j2rVrsLGxgY2NjewxHB0dYWxsjGvXrsnG2drawtzcXO55R40ahaNHj+L+/fsAgPXr12Po0KG8FjFVOtxGiIgqp1fPFEYKp6WlhUGDBsHPzw+HDh1CSEgIzp07h1WrVmHNmjUICAjAjBkzYGdnJ3ZUIlJyenp6yM3NFe25y5uWlpbs/5cf6EtKSgAAubm5mDt3Lnr16vXK/XR1dcv8HFWrVn1lnJubG1xdXbFx40Z4e3vjypUrOHjw4LvGp0qI28jfuI0QESkWC28RSaVSfPzxx/j4449x5swZLFiwAIcPH8a6deuwadMmDBkyBDNmzEDt2rXFjkpESkoikZT6IVkVNWnSBMnJyahXr16p0xs2bIi7d+/i7t27sj16V69eRVZWFhwdHd/6+CNHjsR3332H+/fvw8vLS26vICkvbiP/4DZCRKQ4PNS8kmjdujV+//13nDt3Dj4+PigqKsKaNWtQv359jB49Grdv3xY7IhFRpTZr1ixs3LgRc+fOxZUrV3Dt2jVs27YNM2bMAAB4eXnB2dkZ/v7+uHjxIqKjoxEQEIB27drBw8PjrY/v5+eHe/fu4ZdffuEJo0gpcRshIhIPC+9KxtPTE4cPH8bZs2fRqVMnFBUV4ZdffkH9+vUxZswY3L17V+yIRESVko+PDw4cOICjR4+iadOmaNGiBZYuXQpbW1sAf+/Z3Lt3L0xMTNC2bVt4eXnBzs4O27dvL9PjGxkZoXfv3tDX10fPnj0VOCdEisFthIhIPBJBUMC1LFRQTk4OjIyMkJ2dDUNDwwp73jNnzmDOnDk4fvw4AEBHRwfjx4/H1KlTYWpqWmE5iEg55OXlITU1VXZ9XSpfHTt2hJOTE77//vs3tnvTehCrP1Enb1rG3EYUqzy2ESKiyqC8+2vu8a7kWrdujWPHjiEyMhLt2rVDfn4+vv32W9StWxeLFi3CixcvxI5IRKTynjx5gt27d+PUqVMIDAwUOw5RpcNthIjozVh4K4k2bdrg5MmTOHToEJydnZGVlYXg4GDUr18fa9euRXFxsdgRiYhUlpubG4YOHYpvvvkG9vb2YschqnS4jRARvRkLbyUikUjg6+uL+Ph4bNiwAbVq1cL9+/cxYsQIuLi4YP/+/eAvB4iIyt+ff/6J7OxsfPHFF2JHIaqUuI0QEb0ZC28lpKGhgYCAACQnJ+Pbb7+FiYkJrl69iu7du6Ndu3a4ePGi2BGJiIiIiIjo/7HwVmK6urqYPHkybt26ha+++gq6uro4ffo0PDw8MHLkSGRkZIgdkYiIiIiISO2x8FYBxsbGCAkJwR9//IGBAwdCEATZNcAXL16M/Px8sSMSERERERGpLRbeKsTGxgZbtmzBmTNn4O7ujqdPn2LKlClwcnLCvn37+PtvIiIiIiIiEbDwVkGtWrVCdHQ01q1bB0tLS9y8eRM9evSAt7c3rly5InY8IiIiIiIitcLCW0VJpVIMHToUf/zxB7766itoa2vj2LFjcHV1xRdffIGnT5+KHZGIiIiIiEgtsPBWcQYGBggJCcG1a9fQs2dPFBcXY8mSJWjYsCF27NjBw8+JiIiIiIgUjIW3mrCzs8Pu3btx6NAh2NnZ4f79++jXrx86d+6MGzduiB2PiAgAcPfuXQwfPhzW1tbQ1taGra0tPv/8czx69EjWJjU1FX5+frC2toauri5q1qyJHj164Pr167I2EokEe/bskbv9cjA0NETTpk2xd+9euedev349jI2NZbeLi4uxcOFCODg4oEqVKjA1NUXz5s2xevVqhc0/0dtwGyEiUk4svNWMr68vLl++jNmzZ0NHRwdHjx5Fo0aNMGvWLLx48QIPnxQhPjkPD58UiR2ViNTMrVu34OHhgRs3bmDr1q1ISUnBypUrcfz4cXh6euLx48coLCxEp06dkJ2djV27diE5ORnbt2+Hs7MzsrKy3vj469atQ1paGmJjY9GqVSv06dMHSUlJr20/d+5cLF26FPPnz8fVq1dx8uRJjB49+q3PQ6Qo3EaIiJSXptgBqOJVqVIFc+bMwaBBgzBu3DgcOXIE8+fPx56Tj2Hu/CUESCCVAEF+pujSSl/suEQkkodPinAvswg1LTRhbqL47iIwMBDa2to4evQoqlSpAgCoVasW3NzcULduXUyfPh2ffvopbt68iePHj8PW1hYAYGtri1atWr318Y2NjWFpaQlLS0vMnz8fy5Ytw8mTJ+Hs7Fxq+3379uGzzz5D3759ZeNcXV3LYU5JVXAb4TZCRFRW3OOtxurVq4fff/8dO3fuhG3dxqjmNBkCJACAEgEI3fqYe76J1NShs7kYOOMBJi/LxMAZD3DobK5Cn+/x48c4cuQIPvvsM1lB8ZKlpSX8/f2xfft2mJubQyqVYufOnSguLn6v5yoqKsKaNWsAANra2q9tZ2lpiRMnTuDhw4fv9Tyk2riNcBshInoXLLzVnEQiQe/evbFl5wlIpBpy00pKgHuZhSIlIyKxPHxShNAtj1Hy/+derIgv4m7cuAFBENCwYcNSpzds2BBPnjyBlpYWvv/+e8yaNQsmJib46KOPMH/+fNy6deutzzFw4EDo6+tDR0cHkyZNQu3atdGvX7/Xtg8NDcXDhw9haWkJFxcXjBkzBr///vt7zyOpDm4jf+M2QkRUdiy8CQBQ39YAUon8OKGkCJMnDEZKSoo4oYhIFPcyi2QFxUslJcD9h4o/AqYsV1oIDAxEeno6wsLC4OnpiR07dsDJyQnh4eFvvN/SpUuRkJCA33//HY6Ojli9ejVMTU1f297R0RGXL1/G+fPnMXz4cGRmZqJbt24YOXLkO88XqRZuI3/jNkJEVHYsvAkAYG6iiSA/U0hlr4gS3IqajZNHf4OzszMWLlyIwkLu/SZSBzUtNF/5Ik4qBWqYK+43rPXq1YNEIsG1a9dKnX7t2jWYmJjA3NwcwN+XSuzWrRsWLFiAS5cuoU2bNvjf//73xuewtLREvXr14O3tjXXr1qF///7IzMx8432kUimaNm2KiRMnYteuXVi/fj3WrFmD1NTU95tRUgncRv7BbYSIqGxYeJNMl1b62DrfGqETLbB9QU2c3D0XXl5eyMvLw9SpU9G0aVPExsaKHZOIFOy/X8RJpUDQQFOFnjyqWrVq6NSpE3766Se8ePFCbtrLPXf9+/eHRCJ55b4SiQQODg549uxZmZ+vWbNmcHd3x4IFC94pp6OjIwC803OR6uE28nrcRoiISsfCm+SYm2iicQNdmJtows7ODkePHsWGDRtgamqKS5cuoXnz5ggODkZeXp7YUYlIgf79RdzW+dYVcoWDH3/8Efn5+fDx8UFkZCTu3r2Lw4cPo1OnTqhRowYWLFiAhIQE9OjRAzt37sTVq1eRkpKCNWvWYO3atejRo8c7Pd/EiRPx888/4/79+6VO79OnD5YuXYoLFy7g9u3bOHXqFAIDA9GgQQM4ODiUxyyTEuM2wm2EiOhdsPCmN5JIJAgICMC1a9cwcOBAlJSUYNGiRXBzc8OFCxfEjkdECvTvL+IqQv369REbGws7Ozv069cPdevWxejRo9GhQwdERUXB1NQUNWvWRO3atTF37lw0b94cTZo0wbJlyzB37lxMnz79nZ6vc+fOqFOnzmv36Pn4+GD//v3o1q0bGjRogCFDhsDBwQFHjx6FpiavxkncRriNEBGVnUQoy1k6CDk5OTAyMkJ2djYMDQ3FjiOavXv34tNPP0VGRgakUim++OILzJ07F7q6umJHIyIAeXl5SE1NRZ06dbhdiuhN64H9ieK9aRlzG6kcuB6IqLIr7/6ae7zpnfTo0QNXr17FoEGDuPebiIiIiIioDFh40zszNTXFpk2bsHfvXlhaWuL69eto2bIlf/tNRERERERUChbe9N66d++OK1euvLL3OyYmRuxoRERERERElQYLb/ogpe399vT0xNy5c3ndbyIiIiIiIrDwpnLycu/3gAEDUFxcjDlz5qB169b4448/xI5GpJZ43kxxcflXflxH4uLyJyJ1w8Kbyo2pqSm2bt2KLVu2wNjYGNHR0WjcuDF++ukndrBEFURDQwMAUFBQIHIS9fb8+XMAgJaWlshJ6L9erpOX64jEwW2EiNQNL7JI5W7gwIFo06YNhg0bhmPHjiEwMBD79u3D2rVrYW1tLXY8IpWmqakJPT09PHz4EFpaWpBK+f1qRRIEAc+fP0dmZiaMjY1lX4RQ5aGhoQFjY2NkZmYCAPT09CCRSEROpT64jRCRuuJ1vMuI1119dyUlJVi+fDmmTJmCvLw8mJqaYuXKlejbt6/Y0YhUWkFBAVJTU1FSUiJ2FLVlbGwMS0vLUgs69ieK97ZlLAgC0tPTkZWVVfHhCMCbtxEiosqgvPtrFt5lxA9K7+/69esYNGgQ4uLiAAD+/v748ccfYWxsLG4wIhVWUlLCw81FoqWl9ca9eOxPFK+sy7i4uJgnAhXB27YRIqLKoLz7ax5qTgrn4OCAqKgo/O9//8OCBQsQFhaGiIgIbNy4ER06dBA7HpFKkkql0NXVFTsGUaWmoaHBApCIiCqE6D/+i4yMRLdu3WBtbQ2JRII9e/bITZdIJKUOixcvlrWpXbv2K9MXLlwo9ziJiYlo06YNdHV1YWNjg0WLFlXE7NH/09LSwty5c3HmzBnUq1cP9+7dQ8eOHTFt2jTubSAiIiIiIpUmeuH97NkzuLq6Yvny5aVOT0tLkxvWrl0LiUSC3r17y7WbN2+eXLvx48fLpuXk5MDb2xu2traIi4vD4sWLMWfOHKxatUqh80avatGiBRISEjBq1CgIgoCQkBC0bt0aN2/eFDsaERERERGRQoh+qLmvry98fX1fO93S0lLu9t69e9GhQwfY2dnJjTcwMHil7UthYWEoKCjA2rVroa2tDScnJyQkJCA0NBSjR4/+8Jmgd1K1alWsWrUK3t7eGDVqFKKjo+Hm5oaffvoJgwYNEjseERERERFRuRJ9j/e7yMjIwMGDBzFixIhXpi1cuBDVqlWDm5sbFi9ejKKiItm0qKgotG3bFtra2rJxPj4+SE5OxpMnT0p9rvz8fOTk5MgNVL769OmDxMREtG3bFk+fPsXgwYMxaNAgLmsiIiIiIlIpSlV4b9iwAQYGBujVq5fc+AkTJmDbtm04efIkPv30U3z99deYMmWKbHp6ejqqV68ud5+Xt9PT00t9rpCQEBgZGckGGxubcp4bAgAbGxucOHEC8+bNg4aGBsLCwtC4cWOcP39e7GhERERERETlQqkK77Vr18Lf3/+VM/UGBQWhffv2cHFxwZgxY7BkyRL88MMPyM/Pf+/nmjp1KrKzs2XD3bt3PzQ+vYaGhgZmzpyJyMhI1K5dG6mpqWjdujUWLFiA4uJiseMRERERERF9EKUpvE+fPo3k5GSMHDnyrW2bN2+OoqIi/PnnnwD+/p14RkaGXJuXt1/3u3AdHR0YGhrKDaRYLVu2REJCAgYOHIji4mLMmDEDHTt2xL1798SORkRERERE9N6UpvBes2YN3N3d4erq+ta2CQkJkEqlsLCwAAB4enoiMjJS7rJV4eHhsLe3h4mJicIy07szMjJCWFgYNmzYAH19fURERKBx48Y4dOiQ2NGIiIiIiIjei+iFd25uLhISEpCQkAAASE1NRUJCAu7cuSNrk5OTgx07dpS6tzsqKgrfffcdLl26hFu3biEsLAyTJk3CoEGDZEW1n58ftLW1MWLECFy5cgXbt2/HsmXLEBQUVCHzSO9GIpEgICAA8fHxaNKkCR49eoSuXbtiypQpvOY3EREREREpHdEL79jYWLi5ucHNzQ3A37/XdnNzw6xZs2Rttm3bBkEQMHDgwFfur6Ojg23btqFdu3ZwcnLCggULMGnSJLlrdBsZGeHo0aNITU2Fu7s7Jk+ejFmzZvFSYpVcvXr1cO7cOdk12RcvXox27drJfSlDRERERERU2UkEQRDEDqEMcnJyYGRkhOzsbP7eWwS7du3C8OHDkZ2dDRMTE6xfvx7du3cXOxYR0Ttjf6J4XMZERPShyrsvEX2PN1FZ9OrVC/Hx8WjatCmePHmCHj16ICgoCAUFBWJHIyIiIiIieiMW3qQ06tSpgzNnzmDSpEkAgKVLl6JNmzZITU0VORkREREREdHrsfAmpaKtrY3Q0FDs3bsXJiYmiI6OhpubG3bv3i12NCIiIiIiolKx8Cal1L17d8THx6NFixbIzs5Gr169MGHCBOTn54sdjYiIiIiISA4Lb1Jatra2iIyMxJdffgkA+OGHH9CqVSvcvHlT5GRERERERET/YOFNSk1LSwuLFi3CgQMHUK1aNcTFxcHd3R179uwROxoREREREREAFt6kIrp27YqEhAS0bNkS2dnZ+OSTTzBlyhQUFRWJHY2IiIiIiNQcC29SGTVr1sSpU6dkZz1fvHgxOnbsiLS0NJGTERERERGROmPhTSpFS0sLoaGh2LFjBwwMDBAZGQk3NzdERESIHY2IiIiIiNQUC29SSX369EFsbCwaNWqEjIwMfPTRR/jmm28gCILY0YiIiIiISM2w8CaV1aBBA5w/fx6DBw9GSUkJvvrqK/Ts2RNZWVl4+KQI8cl5ePiEvwEnIiIiIiLF0hQ7AJEiVa1aFRs2bEDr1q0xfvx47Nu3D+26z4SZ8xcQBAmkEiDIzxRdWumLHZWIiIiIiFQU93iTypNIJBg9ejTOnTuHug5NYeoYBEGQAABKBCB062Pu+SYiIiIiIoVh4U1qw93dHWs3HYBEqiE3vqQEuP+QhTcRERERESkGC29SKw3rmkIqkR8nCMUoevFAnEBERERERKTyWHiTWjE30USQnymk///KF0qKkXxqKjq1d8fvv/8ubjgiIiIiIlJJLLxJ7XRppY+t860ROtECyyZowtb4NrKystC1a1fMmzcPJSUlYkckIiIiIiIVwsKb1JK5iSYaN9CFi6MNTp06hTFjxkAQBMyePRs9evRAVlaW2BGJiIiIiEhFsPAmtaejo4MVK1Zg7dq10NHRwYEDB9C0aVMkJSWJHY2IiIiIiFQAC2+i/zds2DCcPXsWtWrVQkpKClq0aIFt27aJHYuIiIiIiJQcC2+if3F3d0dcXBy8vLzw/PlzDBw4EEFBQSgsLBQ7GhERERERKSkW3kT/YWZmhsOHD+Orr74CACxduhSdOnVCRkaGyMmIiIiIiEgZsfAmKoWGhgZCQkLw22+/QV9fHxEREXB3d8f58+fFjkZEREREREqGhTfRG/Tq1QsxMTFwcHDA/fv30bZtW6xcuRKCIIgdjYiIiIiIlAQLb6K3cHBwQHR0NHr16oXCwkKMHTsWI0aMwIsXL8SORkRERERESoCFN1EZGBgYYOfOnVi4cCGkUinWrVuH1q1b486dO2JHIyIiIiKiSo6FN1EZSSQSBAcH48iRI6hWrRouXrwIDw8PREREiB2NiIiIiIgqMRbeRO/Iy8sLcXFxcHNzw8OHD9GxY0f88MMP/N03ERERERGVioU30XuwtbXFmTNn4Ofnh+LiYkyYMAHDhg1DXl6e2NGIiIiIiKiSYeFN9J709PSwefNmfPvtt5BKpdiwYQPatm2Le/fuiR2NiIiIiIgqERbeRB9AIpFg8uTJOHLkCExNTRETEwN3d3ecOXNG7GhEROUuMjIS3bp1g7W1NSQSCfbs2SM3XSKRlDosXrxY1qZ27dqvTF+4cKHc4yQmJqJNmzbQ1dWFjY0NFi1aVBGzR0REpDAsvInKgZeXF2JjY+Hi4oLMzEx06NABK1as4O++iUilPHv2DK6urli+fHmp09PS0uSGtWvXQiKRoHfv3nLt5s2bJ9du/Pjxsmk5OTnw9vaGra0t4uLisHjxYsyZMwerVq1S6LwREREpkqbYAYhURZ06dXDu3DmMGDEC27dvx2effYa4uDgsX74cOjo6YscjIvpgvr6+8PX1fe10S0tLudt79+5Fhw4dYGdnJzfewMDglbYvhYWFoaCgAGvXroW2tjacnJyQkJCA0NBQjB49+sNngoiISATc401UjqpWrYqtW7fim2++gVQqxZo1a9C+fXs8ePBA7GhERBUqIyMDBw8exIgRI16ZtnDhQlSrVg1ubm5YvHgxioqKZNOioqLQtm1baGtry8b5+PggOTkZT548KfW58vPzkZOTIzcQERFVJiy8icqZRCLBlClTcOjQIRgbG+P8+fNwd3dHVFSU2NGIiCrMhg0bYGBggF69esmNnzBhArZt24aTJ0/i008/xddff40pU6bIpqenp6N69epy93l5Oz09vdTnCgkJgZGRkWywsbEp57khIiL6MCy8iRTEx8cHsbGxaNSoEdLT09GuXTv88ssvYsciIqoQa9euhb+/P3R1deXGBwUFoX379nBxccGYMWOwZMkS/PDDD8jPz3/v55o6dSqys7Nlw927dz80PhERUbli4U2kQHXr1kVUVBR69+6NwsJCjB49GmPGjEFBQYHY0YiIFOb06dNITk7GyJEj39q2efPmKCoqwp9//gng79+JZ2RkyLV5eft1vwvX0dGBoaGh3EBERFSZsPAmUjB9fX3s2LEDCxYsgEQiwc8//4yPPvrotYdMEhEpuzVr1sDd3R2urq5vbZuQkACpVAoLCwsAgKenJyIjI1FYWChrEx4eDnt7e5iYmCgsMxERkSKx8CaqABKJBNOmTcOBAwdgZGSEs2fPwt3dHdHR0WJHIyIqs9zcXCQkJCAhIQEAkJqaioSEBNy5c0fWJicnBzt27Ch1b3dUVBS+++47XLp0Cbdu3UJYWBgmTZqEQYMGyYpqPz8/aGtrY8SIEbhy5Qq2b9+OZcuWISgoqELmkYiISBFYeBNVoC5duiAmJgYNGzbEgwcP0KZNG6xbt07sWEREZRIbGws3Nze4ubkB+Pv32m5ubpg1a5aszbZt2yAIAgYOHPjK/XV0dLBt2za0a9cOTk5OWLBgASZNmiR3jW4jIyMcPXoUqampcHd3x+TJkzFr1ixeSoyIiJSaRBAEQewQyiAnJwdGRkbIzs7mb8fogz19+hQBAQHYs2cPACAwMBBLly6FlpaWuMGISOHYnygelzEREX2o8u5LuMebSAQGBgb47bffMG/ePADA8uXL0bFjR2RmZoqcjIiIiIiIyhsLbyKRSKVSzJw5E/v27YOhoSFOnz4NDw8PxMXFiR2NiIiIiIjKEQtvIpF169YN0dHRsLe3x927d9G6dWuEhYWJHYuIiIiIiMoJC2+iSsDe3h4XLlzAxx9/jLy8PAwaNAiTJ09GUVGR2NGIiIiIiOgDsfAmqiSMjIywd+9ezJgxAwAQGhoKX19fPHr0SORkRERERET0IVh4E1UiUqkU8+fPx86dO1G1alUcO3YMTZs2RWJiotjRiIiIiIjoPbHwJqqEevfujaioKNjZ2SE1NRWenp7YuXOn2LGIiIiIiOg9sPAmqqScnZ0RExODTp064fnz5+jbty+mT5+O4uJisaMREREREdE7YOFNVImZmpri0KFD+OKLLwAAX3/9NXr06IHs7GyRkxERERERUVmJXnhHRkaiW7dusLa2hkQiwZ49e+SmDx06FBKJRG7o3LmzXJvHjx/D398fhoaGMDY2xogRI5CbmyvXJjExEW3atIGuri5sbGywaNEiRc8aUbnQ1NTE4sWLsXnzZujq6uLgwYNo1qwZrl+/LnY0IiIiIiIqA9EL72fPnsHV1RXLly9/bZvOnTsjLS1NNmzdulVuur+/P65cuYLw8HAcOHAAkZGRGD16tGx6Tk4OvL29YWtri7i4OCxevBhz5szBqlWrFDZfROXN398fZ8+ehY2NDf744w80a9YM+/fvFzsWEREREZHK+etJ+V7WV7NcH+09+Pr6wtfX941tdHR0YGlpWeq0a9eu4fDhw4iJiYGHhwcA4IcffkCXLl3w7bffwtraGmFhYSgoKMDatWuhra0NJycnJCQkIDQ0VK5AJ6rsmjRpgtjYWPTt2xeRkZHo3r075s2bh+nTp0MqFf17NCIiIiIipXfobC4WrU8r18dUik/qp06dgoWFBezt7TF27Fi56xpHRUXB2NhYVnQDgJeXF6RSKS5cuCBr07ZtW2hra8va+Pj4IDk5GU+ePCn1OfPz85GTkyM3EFUGFhYWOHbsGMaNGwcAmDVrFvr27YunT5+KnIyIiIiISLk9eJiHJWGPUCKU7+NW+sK7c+fO2LhxI44fP45vvvkGERER8PX1lZ3ZOT09HRYWFnL30dTUhKmpKdLT02VtqlevLtfm5e2Xbf4rJCQERkZGssHGxqa8Z43ovWlpaeGHH37AmjVroK2tjV27dsHT0xM3b94UOxoRERERkdJJTU3F9OnT0eajPhAgKffHr/SF94ABA9C9e3c4OzujZ8+eOHDgAGJiYnDq1CmFPu/UqVORnZ0tG+7evavQ5yN6H8OHD0dERASsrKxw5coVeHh44OjRo2LHIiIiIiKq9AoLC7Fr1y507twZdevWxddff437qXEQhPK/fG+lL7z/y87ODmZmZkhJSQEAWFpaIjMzU65NUVERHj9+LPtduKWlJTIyMuTavLz9ut+O6+jowNDQUG4gqoxatGiB2NhYtGjRAllZWfD19cW3334LQSjn42OIiIiIiFTArVu3MG3aNNjY2KB37944cuQIBEGAl5cXNq//AZMGmqC8T5+kdIX3vXv38OjRI1hZWQEAPD09kZWVhbi4OFmbEydOoKSkBM2bN5e1iYyMRGFhoaxNeHg47O3tYWJiUrEzQKQA1tbWOHXqFEaMGIGSkhJ8+eWXGDBoLM4nZuFhOZ+RkYiIiIhI2eTn5+PXX39Fp06dULduXYSEhCAjIwMWFhb46quvkJKSgvDwcPTp0wfd2xpj7XSrcn1+0c9qnpubK9t7Dfx9bH1CQgJMTU1hamqKuXPnonfv3rC0tMTNmzcxZcoU1KtXDz4+PgCAhg0bonPnzhg1ahRWrlyJwsJCjBs3DgMGDIC1tTUAwM/PD3PnzsWIESMQHByMy5cvY9myZVi6dKko80ykCDo6Ovjll1/QpEkTfP3TGWQaBGPayhxIJNmY7FcNXVrpix2RiIiIiKhCXblyBatXr8amTZvkTtLt7e2N0aNHo1u3bnIn4X7JzKR8S2WJIPLxqKdOnUKHDh1eGT9kyBCsWLECPXv2RHx8PLKysmBtbQ1vb2/Mnz9f7mRpjx8/xrhx47B//35IpVL07t0b33//PfT1/yk0EhMTERgYiJiYGJiZmWH8+PEIDg4uc86cnBwYGRkhOzubh51TpfbwSREGTL8vd1IIiUTAtv/VgHk5v4EQ0btjf6J4XMZEROotNzcX27dvx+rVq3H+/HnZeGtrawwfPhzDhg2DnZ3dGx+jvPsS0QtvZcFOnJRFfHIeJi/LfGV8G9uzmDNlACSS8j9LIxGVHfsTxeMyJiJSP4IgIDo6GqtXr8a2bduQm5sLANDQ0EC3bt0wcuRI+Pj4QFOzbDuiyrsv4e4vIhVT00ITUgnkrj0olBThm/lBuJscjp9++gm6urriBSQiIiIiKiePHj3C5s2bsXr1aly+fFk2vn79+hgxYgSGDBny2hNqVyQW3kQqxtxEE0F+pgjd+hglJYBUCjSySEDki0ysW7cOV69exW+//YYaNWqIHZWIiIiI6J2VlJTg5MmTWL16NXbt2oWCggIAgK6uLvr27YuRI0eiTZs2lepITxbeRCqoSyt9NHXUxf2HRahhrglzk174uPUR9OvXDxcuXICHhwd+++03tGzZUuyoRERERERlcv/+faxfvx5r1qxBamqqbHzjxo0xatQo+Pn5wdjYWLyAb8DCm0hFmZtoyp1MzcvLC7GxsejZsyeSkpLQvn17/PTTTxg5cqSIKYmIiIiIXq+wsBCHDh3C6tWrcejQIZSUlAAADA0N4e/vj5EjR6JJkyYip3w7Ft5EasTOzg7nzp3DsGHDsHPnTowaNQoXL17Ed999V+plFIiIiIiIxHDt2jWsW7cOmzZtQnp6umx8mzZtMHLkSPTp0wd6enoiJnw3LLyJ1Iy+vj5+/fVXhISEYMaMGVixYgUuX76MHTt2yF2mj4iIiIioImVlZWHbtm1Yt24doqOjZeMtLCwwdOhQDB8+HPb29iImfH8svInUkEQiwbRp0+Dq6go/Pz+cPn0aHh4e2L17Nzw8PMSOR0RERERqori4GCdOnMC6deuwe/du5OXlAfj7MmBdu3bFsGHD0KVLF6U/OpOFN5Ea69q1K6Kjo9GjRw8kJyejTZs2WLVqFQYPHix2NCIiIiJSYTdv3sT69euxYcMG3L17VzbeyckJw4YNw6BBg1TqaEwW3kRqzt7eHhcuXMDgwYOxf/9+BAQEID4+HosWLYKmJt8iiIiIiKh85ObmYufOnVi3bh0iIyNl442NjeHn54dhw4bB3d29Ul0GrLxIxQ5AROIzMjLCnj17MHPmTADA0qVL0blzZzx69EjkZERERESkzARBQGRkJIYNGwZLS0sMGzYMkZGRkEgk8PHxwbZt25CWlobly5fDw8NDJYtugHu8iej/SaVSzJs3D40bN0ZAQACOHz8ODw8P7N27Fy4uLmLHIyIiIiIlcvfuXWzYsAHr16/HzZs3ZePr1auHYcOGISAgADVr1hQxYcVi4U1Ecnr16oUGDRqgR48euHXrFjw9PbF+/Xr07dtX7GhEREREVInl5uZi165d2LhxI06cOAFBEAD8fVWdfv36YdiwYWjVqpXK7tV+Ex5qTkSvaNSoEWJiYuDt7Y3nz5+jX79+mDZtGoqLi8WORkRERESVSHFxMcLDwzF48GBUr14dQ4YMwfHjxyEIAtq1a4f169cjPT0da9asQevWrdWy6Aa4x5uIXsPU1BQHDx7EtGnTsHjxYoSEhODSpUsICwuDsbGx2PGIiIiISESXL1/Gxo0bERYWhgcPHsjG169fH4MHD8agQYNQp04dERNWLiy8iei1NDU1sWjRIjRu3BgjRozAoUOH0KxZM+zduxcNGzYUOx4RERERVaCMjAxs2bIFmzZtQnx8vGy8iYkJBgwYgICAADRv3lxt92q/CQtvInorPz8/ODg44JNPPsGNGzfQvHlzbN68Gd27dxc7GhEREREp0IsXL7Bv3z5s3LgRR44ckf30UEtLC127dkVAQAC6dOkCHR0dkZNWbiy8iahMmjRpgtjYWPTt2xcRERHo0aMH5s6dixkzZkAq5ekiiIiIiJTdwydFuJdZBGszKZKvnMfGjRuxY8cO5OTkyNo0b94cAQEB6N+/P6pVqyZiWuXCwpuIyszc3Bzh4eEICgrCjz/+iNmzZyMhIQEbNmyAgYGB2PGIiIiI6D0dOpuLJVseQRAkEIRiJJ/6CWnXtwMAbG1tZb/btre3Fzmpcnqv3VQhISHlnYOIlISWlhZ++OEHrFmzBtra2ti9ezdatGiBlJQUsaMR0X+wvyYiore5f/8+FixagcWb/4Ig/P3bbIlEA/btQjBkxEScOnUKt27dwvz581l0f4D3Krx37twJAGjdunW5hiEi5TF8+HBERETAysoKV69eRdOmTXHkyBGxYxHRv7C/JiKi0mRlZWHNmjX46KOPYGNjg2+XbYJEIl8aSqQamDglBO3atePPCsvBey1BDw8P+Pr64s6dO9i1axf++OMP2cXRiUh9tGjRAnFxcfD09ERWVha6dOmCRYsW8f2AqJJgf01ERC+9ePECO3fuRK9evVC9enWMHDkSJ0+ehCAIcGpgBkC+f5BKgRrm/GVyeZEI79kDJyYmokuXLhg4cCCSkpJw48YNmJiYoFGjRli/fn05xxRfTk4OjIyMkJ2dDUNDQ7HjEFUq+fn5GDduHFavXg0AGDBgANasWQM9PT2RkxFVPhXdn6hbfw2wzyYieqmoqAgnT57Eli1bsGvXLrmTpDk5OcHf3x8DBw5E7dq1cehsLkK3PkZJyd9Fd9BAU3RppS9ienGVd1/yTl9hTJ8+HdOnT4eenh5cXFwQHh4udy3fR48eISkp6YNDEZFy0dHRwapVq9CkSRNMmDAB27Ztw/Xr17F7927Url1b7HhEaof9NRGR+hIEATExMdiyZQu2bduGjIwM2bRatWph4MCB8Pf3h7Ozs9z9urTSR1NHXdx/WIQa5powN+He7vL0Tnu8W7Rogbt372LBggUYOnSoAmNVPvz2nKhsIiMj0adPHzx8+BBmZmb49ddf0aFDB7FjEVUaFdGfqHN/DbDPJiL1dOXKFWzfvh1bt26VO+mtqakp+vXrB39/f7Rs2ZK/1y6j8u5L3vlQ882bN2PatGkwNzfHd999hzZt2nxwCGXATpyo7O7evYtPPvkEcXFx0NDQQGhoKMaPHw+JRCJ2NCLRVVR/oq79NcA+m4jUR3JyMrZv345ff/0VV65ckY3X09NDjx494O/vj06dOkFbW1vElMqpvPuSd/66Y9CgQUhOTka3bt3g6+uL3r1749atWx8chIhUh42NDU6fPo1BgwahuLgYn3/+OQICAvD8+XOxoxGpDfbXRESq6datWwgJCUHjxo3h4OCA2bNn48qVK9DS0kK3bt2wefNmZGRkYMuWLejatSuL7krivY4zqFKlCubMmYPk5GTo6emhUaNGCA4OxuXLl1FcXFzeGYlICVWpUgUbN25EaGgoNDQ0sHnzZnh6euLmzZtiRyNSG+yviYhUw507d/Dtt9+iadOmqFu3LqZNm4ZLly5BU1MTnTt3xrp165CZmYl9+/bB398f+vrqe1K0yuqdDjXPz8/H2bNncf36dSQnJyM5ORnXr1/H3bt3Afz9Q34dHR04OjoiLi5OYaHFwMPWiN5fREQE+vXrh8zMTBgbGyMsLAxdunQROxaRKCqiP1Hn/hpgn01EquHBgwfYsWMHtm/fjqioKNl4qVSKjz76CP3798cnn3yCatWqiZhSdYl6VvMOHTogPj4erq6uaNCgAdq0aYMRI0agQYMGaNCgAfLy8pCQkIDExMQPDkZEqqNdu3a4ePEi+vTpg/Pnz+Pjjz/G7NmzMXPmTJ7gg0gB2F8TESmnjIwM7Ny5E7/++itOnz6Nl/tIJRIJ2rZti/79+6N3796wsLAQOSm9q3fa421vb4/t27ejcePGCoxUOfHbc6IPV1BQgEmTJuGnn34CAHTt2hWbNm2CiYmJyMmIKk5F9Cfq3F8D7LOJSLncv38fu3btwm+//YbTp0+jpKRENq1ly5bo378/+vTpA2traxFTqh/Rz2qurtiJE5WfjRs34tNPP0VeXh7s7Oywe/duuLi4iB2LqEKwP1E8LmMiquz+/PNP/Pbbb9i5cyfOnz8vN61p06bo378/+vXrBxsbG5ESkqiHmhMRlYeAgAA4OzujV69euHXrFlq0aIFffvkF/v7+YkcjIiIiUog//vgDv/32G3777bdXzq/RqlUr9O7dG7169YKtra1ICUmRWHgTkSjc3NwQFxcHf39/HD58GIMGDcKFCxfw7bff8rIXREREpPQEQcCVK1dke7YvX74smyaVStGuXTv07t0bn3zyCQ8jVwMsvIlINKampjhw4ADmzp2L+fPn44cffsDFixexY8cOWFlZiR2PiIiI6J0IgoD4+HhZsf3HH3/IpmlqaqJjx47o3bs3evbsCXNzcxGTUkVj4U1EotLQ0MC8efPQtGlTDB48GGfPnkWTJk2wY8cOtG7dWux4RERERG9UVFSEs2fPYu/evdizZw9SU1Nl03R0dODt7Y3evXuje/fuPKGsGmPhTUSVQrdu3RATE4NevXrh8uXL6NChA5YsWYLx48dDIpGIHY+IiIhI5vnz5zh69Cj27NmDAwcO4NGjR7JpVapUQZcuXdC7d2907dqVJ3kkACy8iagSqV+/Ps6fP49Ro0Zh69at+Pzzz3HhwgWsWrUKVatWFTseERERqbGHDx/iwIED2LNnD44ePYq8vDzZNFNTU3Tr1g09evSAt7c3P7fQK1h4E1GlUrVqVYSFhaFZs2b44osvsGXLFiQlJWHXrl2oV6+e2PGIiIhIjaSkpGDv3r3Yu3cvzp49K3eN7dq1a6Nnz57o2bMnWrVqBU1Nllb0enx1EFGlI5FIMHHiRDRp0gT9+vVDUlISPDw8sHnzZnz88cdixyMiIiIV8vBJEe5lFqGmhSbMjDUQFxeHPXv2YM+ePbhy5Ypc2yZNmqBHjx7o2bMnnJ2d+XM4KjMW3kRUabVt2xZxcXHo27cvoqKi0K1bN8yYMQNz5syBhoaG2PGIiIhIyR06m4slWx5DEAAIJUi/tBDXon6WTdfQ0ED79u3Rs2dPdO/eHbVq1RIvLCk1iSAIgtghlEFOTg6MjIyQnZ3NEyQQVbCCggIEBQVh+fLlAICOHTtiy5YtsLCwEDkZ0btjf6J4XMZE9Da3b9/Gjj3Hcej6R4BEKhsvlBTh0m5vdGzrhp49e6JLly48E7maKu++hHu8iajS09bWxo8//oiWLVti1KhROH78OJo0aYLt27ejgWNz2eFh5iZ8SyMiIqJXFRcX4/z58zhw4AAOHjyIpKQkGFt7wq2Hl1w7iVQTh49fRLNG/NKOyhc/pRKR0vDz80Pjxo3Ru3dvXL9+Hf0/XY4GbWsCkEAqAYL8TNGllb7YMYmIiKgSePLkCQ4fPoyDBw/i999/x+PHj2XTpFIpGjUwByAAkPxrPFCnhl7FhyWVx8KbiJSKo6MjYmJiMHz0F3hoGIyXnWWJAIRufYymjrrc801ERKSGBEHAtWvXcODAARw4cADnzp1DcXGxbLqJiQk6d+6Mjz/+GD4+PqhWrRoOnc1F6NbHKCn5u+gOGmjKzxGkEHxVEZHS0dfXx1ezl+KLZQ/lxpeUAPcfFrHDJCIiUhM5OTk4fvw4Dh8+jCNHjuD27dty052cnNC1a1d8/PHH8PT0fOWSX11a6aOpoy7uPyxCDXP+bI0Uh68sIlJKNhZakEr+3tP9klBSjLOndqNxg4HiBSMiIiKFKSkpQUJCAg4fPozDhw8jKioKRUVFsuk6Ojro0KEDPv74Y3Tt2hW1a9d+62Oam7DgJsXjK4yIlJK5iSaC/Exlh4dBKEZyxFSM+3k7EmJO4Pvvv0eVKlXEjklEREQf6OHDhzh69CiOHDmCI0eOIDMzU256/fr10blzZ3Tu3Bnt2rVD1apVRUpK9HosvIlIaf378DCralL8sqIRZs36FatXr0ZcXBx27NiBunXrih2TiIiI3kFRUREuXLgg26sdFxeHf18BuWrVqujYsSM6d+4MHx8f2NnZiZiWqGxYeBORUvv34WEzZsxAixYtMHDgQMTHx8Pd3R0bNmxAjx49RE5JREREryMIAq5fv45jx47h2LFjOHXqFHJycuTauLq6yvZqt2zZEtra2iKlJXo/0rc3UazIyEh069YN1tbWkEgk2LNnj2xaYWEhgoOD4ezsjKpVq8La2hoBAQF48OCB3GPUrl0bEolEbli4cKFcm8TERLRp0wa6urqwsbHBokWLKmL2iKiCeXl5IT4+Hi1btkR2djZ69uyJKVOmyP3+i4iIiMSVlpaGzZs3Y+jQobCxsYGjoyMmTJiAffv2IScnB6amphgwYADWr1+PBw8eICEhAQsXLkT79u1ZdJNSEn2P97Nnz+Dq6orhw4ejV69ectOeP3+OixcvYubMmXB1dcWTJ0/w+eefo3v37oiNjZVrO2/ePIwaNUp228DAQPZ/Tk4OvL294eXlhZUrVyIpKQnDhw+HsbExRo8erdgZJKIKV7NmTZw6dQrBwcFYunQpFi9ejPPnz2P79u2wsrISOx4REZHaefr0KSIiImR7ta9cuSI3XUdHB23atIGXlxe8vLzQuHFjaGhoiJSWqPyJXnj7+vrC19e31GlGRkYIDw+XG/fjjz+iWbNmuHPnDmrVqiUbb2BgAEtLy1IfJywsDAUFBVi7di20tbXh5OSEhIQEhIaGsvAmUlFaWloIDQ1Fq1atMGzYMJw+fRpubm7YunUrOnToIHY8IiIilVZQUIALFy7g+PHjOHbsGC5cuCB39JlEIkGTJk1khXarVq14UlRSaaIX3u8qOzsbEokExsbGcuMXLlyI+fPno1atWvDz88OkSZNk1+mLiopC27Zt5Q5L8fHxwTfffIMnT57AxMTklefJz89Hfn6+7PZ/f2dCRMqhd+/ecHZ2Rp8+fZCUlAQvLy/Mnj0b06dP5zfpRERE5aSgoAAxMTE4deoUTp48iXPnzuHFixdyberWrSsrtDt06IBq1aqJlJao4ilV4Z2Xl4fg4GAMHDgQhoaGsvETJkxAkyZNYGpqinPnzmHq1KlIS0tDaGgoACA9PR116tSRe6zq1avLppVWeIeEhGDu3LkKnBsiqigNGjTA+fPnMX78eKxduxazZ89GREQEwsLCXnukDBEREb1eQUEBYmNjZYX22bNnXym0zc3N0aFDB3Tq1AkdO3Z85fM4kTpRmsK7sLAQ/fr1gyAIWLFihdy0oKAg2f8uLi7Q1tbGp59+ipCQEOjo6LzX802dOlXucXNycmBjY/N+4YlIdHp6elizZg3at2+PMWPG4MSJE2jcuDG2bNmCjz76SOx4RERElVphYaGs0D516hTOnDmD58+fy7UxMzND+/bt0aFDB7Rv3x4NGzaERCIRKTFR5aIUhffLovv27ds4ceKE3N7u0jRv3hxFRUX4888/YW9vD0tLS2RkZMi1eXn7dXu7dHR03rtoJ6LKa/DgwfDw8EC/fv1w+fJleHl5YdasWZg5cyYPPSciIvp/L168QHR0NE6fPo3Tp0/j7NmzePbsmVybatWqyRXajo6OLLSJXqPSF94vi+4bN27g5MmTZfotSEJCAqRSKSwsLAAAnp6emD59OgoLC6GlpQUACA8Ph729famHmRORamvYsCEuXLiAzz//HKtXr8bcuXMRGRmJLVu28NBzIiJSCw+fFOFeZhFqWmjC3EQTjx49wtmzZ3HmzBmcOXMGsbGxKCwslLtPtWrV0K5dO1mx7ejoCKlU9KsTEykF0Qvv3NxcpKSkyG6npqYiISEBpqamsLKyQp8+fXDx4kUcOHAAxcXFSE9PBwCYmppCW1sbUVFRuHDhAjp06AADAwNERUVh0qRJGDRokKyo9vPzw9y5czFixAgEBwfj8uXLWLZsGZYuXSrKPBOR+PT09PDLL7+gXbt2GDNmDE6ePAlXV1eEhYXBy8tL7HhEREQKc/DsU4RueQxBkABCCbL/WIaLJ757pZ2VlRXatGmD1q1bo3379nBycmKhTfSeJIIgCGIGOHXqVKmX9hkyZAjmzJnz2pMwnDx5Eu3bt8fFixfx2Wef4fr168jPz0edOnUwePBgBAUFyR0qnpiYiMDAQMTExMDMzAzjx49HcHBwmXPm5OTAyMgI2dnZbz3UnYiUy/Xr19GvXz8kJSVBIpFg5syZmDVrFg89J4Vgf6J4XMZE8oqKipCUlISoqCicOpOIh0bBkEj+6eOEkiJEbW6FOjbGskK7TZs2qF27Ng8dJ7VV3n2J6IW3smAnTqTaXrx4gc8//xy//PILAKB9+/bYsmULrKysRE5Gqob9ieJxGZO6y8jIwPnz5xEVFYXz588jJiZGdiI0Y2tPuPXY9sp9Zg3VRvtm/LkV0Uvl3ZeIfqg5EVFlUKVKFaxatQrt27fHp59+ilOnTqFx48bYvHkzOnXqJHY8IiKiUhUWFuLSpUuyIjsqKgqpqamvtDMyMkKzZs3g4emNCw8FCPhnT7ZUCjjVN6vI2ERqh4U3EdG/+Pn5wd3dHf369UNiYiJ8fHwwffp0zJ49G5qafMskIiLxCIKA+/fvIzo6WlZox8bGIi8vT66dRCKBo6MjPD090aJFC3h6esLBwUH2++xDZ3MRuvUxSkr+LrqDBprC3IR9HJEi8VDzMuJha0Tq5cWLF5g0aRJ+/vlnAECbNm0QFhYGGxsbkZORsmN/onhcxqQqHj58iJiYGMTGxsr+vjzR8L+ZmJigRYsWsiK7WbNmMDIyevNjPynC/YdFqGGuyaKbqBQ81JyIqAJUqVIFK1euRPv27TFq1CicPn0ajRs3xtq1a9GjRw+x4xERkYrJyspCXFycXKF9586dV9ppaGjAyclJbm92/fr13/ls4+YmLLiJKhK3NiKiNxgwYAA8PDwwcOBAxMbGomfPnggMDMTixYtRpUoVseMREZESys7OxqVLlxAXFycrsm/cuPFKO4lEAnt7e3h4eKBp06bw8PBA48aNoaenJ0JqIvoQLLyJiN6iXr16OHv2LKZPn45vv/0Wy5cvR2RkJLZt2wZHR0ex4xERUSUlCAIePHiAhIQExMfHy/7eunWr1PZ16tSRFdhNmzZFkyZN+HMJIhXBwpuIqAy0tbWxePFieHl5ISAgAElJSfDw8MCyZcswcuRIXueUiEjNFRcX48aNG3IFdkJCAh4+fFhq+1q1aqFx48Zo2rSprNiuVq1aBacmoorCwpuI6B34+Pjg0qVLCAgIQHh4OEaPHo3w8HCsWrUKxsbGYscjIqIK8OjRI1y+fBlJSUm4fPkyEhISkJiYiBcvXrzSVkNDAw4ODnBzc0Pjxo3h5uYGV1dXFtlEaoaFNxHRO7K0tMThw4exZMkSTJs2DTt27EB0dDS2bNmCli1bih2PiIjKyfPnz3H16lVZgf3yb1paWqnt9fT04OrqKiuwGzdujEaNGvGcIETEwpuI6H1IpVJ8+eWXaNeuHQYOHIhbt26hbdu2mDt3Lr766itoaGiIHZGIiP7l4ZMi3MssQk2LV8/mXVBQgBs3buDKlStyRfatW7fwuivv1q5dG87OzmjUqBFcXFzg5uaGevXq8f2fiErFwpuI6AM0a9YM8fHxGDt2LLZs2YIZM2bg+PHj2Lx5M6ytrcWOR1SuIiMjsXjxYsTFxSEtLQ27d+9Gz549ZdOHDh2KDRs2yN3Hx8cHhw8flt1+/Pgxxo8fj/3790MqlaJ3795YtmwZ9PX1ZW0SExMRGBiImJgYmJubY/z48ZgyZYrC549U16GzuQjd8hglAiCRCGhR+zqK/jqB69ev4/r167h58yaKi4tLva+FhQUaNWokK7IbNWoEJycnGBgYVPBcEJEyY+FNRPSBDA0NsXnzZnh7eyMwMBAnT56Ei4sL1q9fj48//ljseETl5tmzZ3B1dcXw4cPRq1evUtt07twZ69atk93W0dGRm+7v74+0tDSEh4ejsLAQw4YNw+jRo7FlyxYAQE5ODry9veHl5YWVK1ciKSkJw4cPh7GxMUaPHq24mSOVUVJSgnv37uH69eu4du0akq7fx018Bonk7+tcC4IE527WR9Tmoch/li67n4GBARwdHWUF9su/FhYWYs0KEakQifC642dITk5ODoyMjJCdnc3LOhDRayUnJ2PgwIGIj48HAEyYMAHffPMNdHV1RU5GlYWq9CcSiaTUPd5ZWVnYs2dPqfe5du0aHB0dERMTAw8PDwDA4cOH0aVLF9y7dw/W1tZYsWIFpk+fjvT0dGhrawMAvvrqK+zZswfXr18vUzZVWcb0eiUlJXjw4AFSUlJw48YNpKSkyA3Pnz+XtTW29oRbj22vPEZt6Ta4N6wKBwcHNGzYEFZWVrxCBRHJlHdfwj3eRETlyN7eHlFRUfjqq6/w3Xff4fvvv8eJEyewZcsWODs7ix2PSOFOnToFCwsLmJiY4KOPPsL//vc/2dmbo6KiYGxsLCu6AcDLywtSqRQXLlzAJ598gqioKLRt21ZWdAN/H67+zTff4MmTJzAxMXnlOfPz85Gfny+7nZOTo8A5pIpSXFyMe/fulVpY37x5E3l5ea+9r5aWFurXrw8HBwfUaeCO+GwBAv4pqqVS4Jv5Qa/81puISFH4bkNEVM50dHSwdOlSdOrUCcOGDcPly5fRtGlTLFy4EBMmTIBUKhU7IpFCdO7cGb169UKdOnVw8+ZNTJs2Db6+voiKioKGhgbS09NfOWxXU1MTpqamSE//+5Df9PR01KlTR65N9erVZdNKK7xDQkIwd+5cBc0VKUpJSQnS0tLw559/yobbt2/L/V9QUPDa+2tqaqJOnTqoX78+6tWrJxvq1q2LOnXqQEtLS9b20NlchG59jJKSv4vuoIGmLLqJqELxHYeISEG6dOki+33qwYMHMWnSJBw6dAjr16/niddIJQ0YMED2v7OzM1xcXFC3bl2cOnUKHTt2VNjzTp06FUFBQbLbOTk5sLGxUdjzUdkUFhYiLS0Nd+7ceaWo/vPPP3Hnzp03FtYAoK2tjbp168oV1i+HWrVqQVOzbB9lu7TSR1NHXdx/WIQa5q+e1ZyISNH4rkNEpEAWFhbYv38/Vq5cicmTJyM8PBwuLi5YvXq13G9jiVSRnZ0dzMzMkJKSgo4dO8LS0hKZmZlybYqKivD48WNYWloCACwtLZGRkSHX5uXtl23+S0dH55WTuJFiPX/+HPfv38f9+/dx79493Lt3T+7/e/fuISMj47WX4npJQ0MDtWrVQu3atWFra4vatWvLDTVr1iy3y3OZm7DgJiLx8N2HiEjBJBIJxo4di/bt28Pf3x/x8fH45JNPMHLkSCxdulTuMkpEquTevXt49OgRrKysAACenp7IyspCXFwc3N3dAQAnTpxASUkJmjdvLmszffp0FBYWyg4VDg8Ph729famHmdObvena1f8lCAJycnKQkZGB9PT0V/4+ePBAVlw/fvy4TM+vpaWFGjVqoE6dOq8U1bVr14a1tXWZ91oTESkzntW8jHiGVCIqDwUFBZg5cyYWL14MQRBQv359hIWFoWnTpmJHowqizP1Jbm4uUlJSAABubm4IDQ1Fhw4dYGpqClNTU8ydOxe9e/eGpaUlbt68iSlTpuDp06dISkqS7ZH29fVFRkYGVq5cKbucmIeHh+xyYtnZ2bC3t4e3tzeCg4Nx+fJlDB8+HEuXLi3z5cSUeRmXp0Nnc7Fky2MI/3/t6o/d/4KNYSoyMzNLLawzMjLeeMKy/9LT00PNmjXlhho1asjdNjMz43ktiEgplXdfwsK7jNiJE1F5OnnyJAICAnDv3j1oampizpw5+Oqrr8rtkEqqvJS5Pzl16hQ6dOjwyvghQ4ZgxYoV6NmzJ+Lj45GVlQVra2t4e3tj/vz5spOjAcDjx48xbtw47N+/H1KpFL1798b3338vd+RHYmIiAgMDERMTAzMzM4wfPx7BwcFlzqnMy7g0giDgxYsXyM7OlhsePXqER48e4a+//nrl/6xnUtRq/ysk0n/eU4SSIkRtbiV37erSGBgYoHr16rC0tJT7a2VlJVdUGxkZ8fJbRKSyWHiLRNU6cSIS35MnTzBmzBj8+uuvAIDWrVtj06ZNqF27trjBSKHYnyjey2V8889HsLM1rfDnFwQBeXl5eP78OZ49e4Znz57J/v/vuNzc3FcK6tKGoqKid8rwumtXF6TMg4X+41IL65d/9fT0ymtREBEpLV7Hm4hIRZiYmGDbtm3o2rUrAgMDcebMGbi6uuKnn36Cv7+/2PGIlN7Q+Q/gWTcKdUzuoqSkRDYUFxfL3f73UFhYiIKCAuTn57/T37y8PFkx/fz5c5SUlJT7/EilUhgaGsLIyAhGRkaoVq2a3GBmZib7X1PXHIt2CRAE+WtX796+kicYIyISAd95iYhEJJFIEBAQgNatW2PQoEGIiorCoEGDcPDgQSxfvpwnkyL6IBJE326IZQtGvvXwakXR0dGBnp4eqlatKhv+fVtPT09WSL9t0NfXf6dDu4t1ee1qIqLKgoealxEPDSQiRSsqKsKCBQswf/58FBcXo0aNGli7di28vb3FjkbliP2J4r1cxm1GXIamtgH0H/8IneJUSKVS2aChoSF3+9/jX16eTFtb+53+/re41tPTE/2M3Q+fFPHa1URE74GHmhMRqShNTU3Mnj0bPj4+CAgIwI0bN+Dj44NRY4Mx7NOvUK+WPj84E70jqRRYs+Jrtd12eO1qIqLKgdd3ICKqZFq0aIH4+HiMGzcOVg79caPkU0z/OQcDpt/HobO5YscjUho8vJqIiCoL9kRERJVQ1apVMWveUlwtvg8Bf/+mU4AE34b9Bdf6GqhhUUXkhESV39rpVrCz1X97QyIiIgXjHm8iokrqXmaRrOj+hxQ9+oxCUlKSKJmIlIkZ93QTEVElwcKbiKiSqmmhCel/6m6hpBiXLx6Hh4cHFi1ahOLiYnHCEREREVGZsfAmIqqkzE00EeRnCun/v1NLpcDoHtrw/qgpCgoKEBwcjHbt2uHmzZviBiUiIiKiN+IxWERElViXVvpo6qgrdzmgAZ33Yt26dZg4cSLOnj0LV1dXLFmyBKNHj36na/wSERERUcXgHm8iokrO3EQTjRvoys7MLJFIMHz4cCQmJqJdu3Z49uwZxowZgy5duuDBgwcipyUiIiKi/2LhTUSkpGrXro0TJ04gNDQUOjo6OHz4MBo1aoSwsDAIgiB2PCIiIiL6fyy8iYiUmFQqxaRJk3Dx4kW4u7vjyZMnGDRoED755BOkp6eLHY+IiIiIwMKbiEglODo6IioqCvPmzYOWlhb27t0LR0dH7v0mIiIiqgRYeBMRqQgtLS3MnDkTsbGxcHNzk+397tmzJ9LS0sSOR0RERKS2WHgTEakYFxcXXLhwAfPnz4eWlhb27dsHJycnbN68mXu/iYiIiETAwpuISAVpaWlhxowZiI2NRZMmTfDkyRMMHjyYe7+JiIiIRMDCm4hIhbm4uOD8+fP43//+x73fRERERCJh4U1EpOK0tLQwffp0xMXFce83ERERkQhYeBMRqQlnZ+dX9n47Ojpiw4YN3PtNREREpEAsvImI1MjLvd8vr/udlZWFoUOHonPnzvjzzz/FjkdERESkklh4ExGpoUaNGuH8+fNYuHAhdHR0cPToUTRq1AjLli1DcXGx2PGIiIiIVAoLbyIiNaWpqYng4GAkJiaibdu2ePbsGSZOnIjWrVvj6tWrYscjIiIiUhksvImI1FyDBg1w8uRJrFixAgYGBjh//jzc3Nwwb948FBQUiB2PiIiISOmx8CYiIkilUowZMwZXr15F165dUVBQgNmzZ8Pd3R3R0dFixyMiIiJSaiy8iYhIpmbNmti/fz+2bNkCMzMzXL58GZ6enpg8eTKePXsmdjwiIiIipcTCm4iI5EgkEgwcOBDXrl2Dv78/SkpKEBoaCmdnZxw/flzseERERERKh4U3ERGVyszMDJs3b8bBgwdhY2OD1NRUeHl5Yfjw4Xj06JHY8YiIiIiUBgtvIiJ6oy5duuDKlSsIDAwEAKxbtw4ODg7YtGkTBEEQOR0RERFR5Sd64R0ZGYlu3brB2toaEokEe/bskZsuCAJmzZoFKysrVKlSBV5eXrhx44Zcm8ePH8Pf3x+GhoYwNjbGiBEjkJubK9cmMTERbdq0ga6uLmxsbLBo0SJFzxoRkcowMDDAjz/+iLNnz8LJyQl//fUXAgIC0KlTp1fek4mIiIhInuiF97Nnz+Dq6orly5eXOn3RokX4/vvvsXLlSly4cAFVq1aFj48P8vLyZG38/f1x5coVhIeH48CBA4iMjMTo0aNl03NycuDt7Q1bW1vExcVh8eLFmDNnDlatWqXw+SMiUiUtW7bExYsX8fXXX0NXVxfHjx+Hs7MzFixYwEuPEREREb2GRKhExwlKJBLs3r0bPXv2BPD33m5ra2tMnjwZX3zxBQAgOzsb1atXx/r16zFgwABcu3YNjo6OiImJgYeHBwDg8OHD6NKlC+7duwdra2usWLEC06dPR3p6OrS1tQEAX331Ffbs2YPr16+XKVtOTg6MjIyQnZ0NQ0PD8p95IiIlc/PmTYwdOxbh4eEAgIYNG2LVqlVo3bq1yMkqN/YnisdlTEREH6q8+xLR93i/SWpqKtLT0+Hl5SUbZ2RkhObNmyMqKgoAEBUVBWNjY1nRDQBeXl6QSqW4cOGCrE3btm1lRTcA+Pj4IDk5GU+ePCn1ufPz85GTkyM3EBHRP+rWrYsjR45g8+bNMDc3x7Vr19CmTRuMHj36te+tREREROqoUhfe6enpAIDq1avLja9evbpsWnp6OiwsLOSma2pqwtTUVK5NaY/x7+f4r5CQEBgZGckGGxubD58hIiIVI5FI4O/vj+vXr2PkyJEAgF9++QUODg7YunUrT75GREREhEpeeItp6tSpyM7Olg13794VOxIRUaVlamqKX375BREREXBwcEBmZib8/Pzg6+uLW7duiR2PiIiISFSVuvC2tLQEAGRkZMiNz8jIkE2ztLREZmam3PSioiI8fvxYrk1pj/Hv5/gvHR0dGBoayg1ERPRmbdu2RUJCAubNmwdtbW0cOXIETk5OmD9/vtxJMYmIiIjUSaUuvOvUqQNLS0scP35cNi4nJwcXLlyAp6cnAMDT0xNZWVmIi4uTtTlx4gRKSkrQvHlzWZvIyEgUFhbK2oSHh8Pe3h4mJiYVNDdEROpBR0cHM2fORFJSEj766CPk5eVh1qxZcHZ2xtGjR8WOR0RERFThRC+8c3NzkZCQgISEBAB/n1AtISEBd+7cgUQiwcSJE/G///0P+/btQ1JSEgICAmBtbS0783nDhg3RuXNnjBo1CtHR0Th79izGjRuHAQMGwNraGgDg5+cHbW1tjBgxAleuXMH27duxbNkyBAUFiTTXRESqr0GDBjh27Bi2bt0KKysrpKSkwMfHB/369cO9e/fEjkdERERUYUQvvGNjY+Hm5gY3NzcAQFBQENzc3DBr1iwAwJQpUzB+/HiMHj0aTZs2RW5uLg4fPgxdXV3ZY4SFhcHBwQEdO3ZEly5d0Lp1a7lrdBsZGeHo0aNITU2Fu7s7Jk+ejFmzZsld65uIiMqfRCLBgAEDcP36dUycOBFSqRQ7duyAg4MDvv32W7kjkYiIiIhUVaW6jndlxmuCEhF9uEuXLuGzzz7DuXPnAABOTk5YsWIF2rRpI3KyisP+RPG4jImI6EOp1XW8iYhItbi6uuL06dNYs2YNzMzMcOXKFbRt2xZDhgzBtRsZiE/Ow8MnRWLHJCIiIipXLLyJiKhCSaVSDB8+HMnJyfj0008hkUgQHp2Pz0KfY/KyTAyc8QCHzuaKHZOIiIio3LDwJiIiUZiammLlypU4cjwa9u0XQiLRAACUCMCSLY+455uIiIhUBgtvIiISlZl1I0gk8t2RIEgwLmg+0tLSREpFREREVH5YeBMRkahqWmhCKpEfJ5QUYe+OVbC3t8e3336LgoICccIRERERlQMW3kREJCpzE00E+ZlC+v89klQK9GnzDI0b2eLp06f48ssv4eLigiNHjogblIiIiOg98XJiZcRLkxARKdbDJ0W4/7AINcw1YW6iiZKSEmzcuBHBwcHIzMwEAHTv3h1Lly6FnZ2dyGnfH/sTxeMyJiKiD8XLiRERkUoyN9FE4wa6MDfRBPD32c+HDh2KP/74A0FBQdDU1MS+ffvg6OiIGTNm4NmzZyInJiIiIiobFt5ERFSpGRkZYcmSJUhMTISXlxfy8/OxYMECODg4YPv27eCBW0RERFTZsfAmIiKl0LBhQxw9ehS7du1C7dq1ce/ePQwYMAAdOnRAfHy82PGIiIiIXouFNxERKQ2JRIJPPvkEV69exdy5c6Grq4uIiAi4u7tjxIgRSE9PFzsiERER0StYeBMRkdKpUqUKZs2aheTkZAwcOBCCIGDt2rWoX78+QkJCkJeXJ3ZEIiIiIhkW3kREpLRq1aqFLVu24OzZs2jWrBlyc3Mxbdo0/v6biIiIKhUW3kREpPRatmyJqKgobN68GTVq1MDt27cxYMAAtGnTBjExMWLHIyIiIjXHwpuIiFSCVCqFv78/kpOTMWfOHOjp6cn2hAcEBOD+/ftiRyQiIiI1xcKbiIhUStWqVTF79mwkJydj8ODBAIBNmzahQYMGmDt3Lp4/fy5yQiIiIlI3LLyJiEgl1axZExs3bsSFCxfQsmVLPH/+HHPmzEGDBg2wYcMGFBcXix2RiIiI1AQLbyIiUmnNmjXDmTNnsG3bNtja2uL+/fsYOnQo3N3dcfToUbHjERERkRpg4U1ERCpPIpGgf//+uH79Or755hsYGRnh0qVL8PHxgY+PDy5duiR2RCIiIlJhLLyJiEht6OrqYsqUKUhJScHnn38OLS0tHD16FG5ubhg6dCju3bsndkQiIiJSQSy8iYhI7ZiZmeG7777DtWvX0K9fPwiCgA0bNqB+/fqYNm0asrOzxY5IREREKoSFNxERqa26deti+/btOH/+PFq3bo28vDyEhISgXr16+PHHH1FQUCB2RCIiIlIBLLyJiEjtNW/eHJGRkdizZw/s7e3x119/Yfz48XBycsJvv/0GQRDEjkhERERKjIU3ERER/j4BW48ePZCUlISffvoJFhYWSElJQZ8+fdCiRQucPHlS7IhERESkpFh4ExER/YuWlhbGjh2LlJQUzJw5E3p6eoiOjsZHH30EHx8fxMXFiR2RiIiIlAwLbyIiolIYGBhg3rx5uHnzJgIDA2VnQPfw8ED//v3xxx9/iB2RiIiIlAQLbyIiojewtLTEjz/+iOvXr2PQoEGQSCT49ddf4ejoiE8//RT3798XOyIRERFVciy8iYiIysDOzg6bNm1CQkICPv74YxQXF2PVqlWoV68epkyZgsePH4sdkYiIiCopFt5ERETvwMXFBfv378fp06fRqlUr5OXlYfHixbCzs8PXX3+NZ8+eiR2RiIiIKhkW3kRERO+hdevWOH36NA4cOAAXFxdkZ2dj+vTpqFu3Ln788Ufk5+eLHZGIiIgqCRbeRERE70kikaBr166Ij49HWFgY7OzskJGRgfHjx6N+/fpYvXo1CgsLxY5JREREImPhTURE9IGkUin8/Pxw7do1/PTTT7C2tsbdu3cxatQoODo6IiwsDMXFxWLHJCIiIpGw8CYiIion2trasmuAh4aGwtzcHCkpKRg0aBBcXFzw22+/IfNxgdgxiYiIqIKx8CYiIipnVapUwaRJk3Dr1i18/fXXMDY2xtWrVzF+xg4Mm58mdjwiIiKqYCy8iYiIFERfXx9Tp05FamoqgqcvhH27EEgk7HqJiIjUDXt/IiIiBTM2Nkb/wZ9DItUQOwoRERGJgIU3ERFRBahpoQmpROwUREREJAYW3kRERBXA3EQTQX6mkLLnJSIiUjuaYgcgIiJSF11a6cOhphXq/iJ2EiIiIqpI/N6diIioApmZ8DtvIiIidcPCm4iIiIiIiEiBWHgTERERERERKRALbyIiIiIiIiIFYuFNREREREREpEAsvImIiIiIiIgUiIU3ERERERERkQKx8CYiIiIiIiJSIBbeRERERERERArEwpuIiIiIiIhIgZSi8K5duzYkEskrQ2BgIACgffv2r0wbM2aM3GPcuXMHXbt2hZ6eHiwsLPDll1+iqKhIjNkhIiJSSpGRkejWrRusra0hkUiwZ88e2bTCwkIEBwfD2dkZVatWhbW1NQICAvDgwQO5xyitT1+4cKFcm8TERLRp0wa6urqwsbHBokWLKmL2iIiIFEZT7ABlERMTg+LiYtnty5cvo1OnTujbt69s3KhRozBv3jzZbT09Pdn/xcXF6Nq1KywtLXHu3DmkpaUhICAAWlpa+PrrrytmJoiIiJTcs2fP4OrqiuHDh6NXr15y054/f46LFy9i5syZcHV1xZMnT/D555+je/fuiI2NlWs7b948jBo1SnbbwMBA9n9OTg68vb3h5eWFlStXIikpCcOHD4exsTFGjx6t2BkkIiJSEKUovM3NzeVuL1y4EHXr1kW7du1k4/T09GBpaVnq/Y8ePYqrV6/i2LFjqF69Oho3boz58+cjODgYc+bMgba2tkLzExERqQJfX1/4+vqWOs3IyAjh4eFy43788Uc0a9YMd+7cQa1atWTjDQwMXttnh4WFoaCgAGvXroW2tjacnJyQkJCA0NBQFt5ERKS0lKLw/reCggJs3rwZQUFBkEgksvFhYWHYvHkzLC0t0a1bN8ycOVO21zsqKgrOzs6oXr26rL2Pjw/Gjh2LK1euwM3N7ZXnyc/PR35+vux2dnY2gL+/iSciInpfL/sRQRBETqJ42dnZkEgkMDY2lhu/cOFCzJ8/H7Vq1YKfnx8mTZoETc2/P5JERUWhbdu2cl+K+/j44JtvvsGTJ09gYmLyyvOwzyYiovJW3v210hXee/bsQVZWFoYOHSob5+fnB1tbW1hbWyMxMRHBwcFITk7Grl27AADp6elyRTcA2e309PRSnyckJARz5859ZbyNjU05zQkREamzR48ewcjISOwYCpOXl4fg4GAMHDgQhoaGsvETJkxAkyZNYGpqinPnzmHq1KlIS0tDaGgogL/75Tp16sg91r/77NIKb/bZRESkKOXVXytd4b1mzRr4+vrC2tpaNu7fh545OzvDysoKHTt2xM2bN1G3bt33ep6pU6ciKChIdjsrKwu2tra4c+eOSn9QysnJgY2NDe7evSv3QUnVqMN8qsM8Auoxn+owj4D6zGd2djZq1aoFU1NTsaMoTGFhIfr16wdBELBixQq5af/uW11cXKCtrY1PP/0UISEh0NHRea/nU9c+uyzUZbsqCy4LeVwe/+Cy+AeXxT/Ku79WqsL79u3bOHbsmGxP9us0b94cAJCSkoK6devC0tIS0dHRcm0yMjIA4LW/MdPR0Sn1A4CRkZFavAgNDQ05nypCHeYRUI/5VId5BNRnPqVSpbiwyDt7WXTfvn0bJ06ceOu6bN68OYqKivDnn3/C3t4elpaWsj76JfbZH05dtquy4LKQx+XxDy6Lf3BZ/KO8+mul6vXXrVsHCwsLdO3a9Y3tEhISAABWVlYAAE9PTyQlJSEzM1PWJjw8HIaGhnB0dFRYXiIiInXysui+ceMGjh07hmrVqr31PgkJCZBKpbCwsADwd58dGRmJwsJCWZvw8HDY29uXepg5ERGRMlCaPd4lJSVYt24dhgwZIjsBCwDcvHkTW7ZsQZcuXVCtWjUkJiZi0qRJaNu2LVxcXAAA3t7ecHR0xODBg7Fo0SKkp6djxowZCAwMfO/D2oiIiNRNbm4uUlJSZLdTU1ORkJAAU1NTWFlZoU+fPrh48SIOHDiA4uJi2XlUTE1Noa2tjaioKFy4cAEdOnSAgYEBoqKiMGnSJAwaNEhWVPv5+WHu3LkYMWIEgoODcfnyZSxbtgxLly4VZZ6JiIjKhaAkjhw5IgAQkpOT5cbfuXNHaNu2rWBqairo6OgI9erVE7788kshOztbrt2ff/4p+Pr6ClWqVBHMzMyEyZMnC4WFhWV+/ry8PGH27NlCXl5eucxPZcX5VB3qMI+CoB7zqQ7zKAicT2Vw8uRJAcArw5AhQ4TU1NRSpwEQTp48KQiCIMTFxQnNmzcXjIyMBF1dXaFhw4bC119//cqyuHTpktC6dWtBR0dHqFGjhrBw4cJ3yqnMy7i8cVn8g8tCHpfHP7gs/sFl8Y/yXhYSQVCD65kQERERERERiUSpfuNNREREREREpGxYeBMREREREREpEAtvIiIiIiIiIgVi4U1ERERERESkQCy8iYiIiIiIiBSIhfdbzJkzBxKJRG5wcHAQO9YHi4yMRLdu3WBtbQ2JRII9e/bITRcEAbNmzYKVlRWqVKkCLy8v3LhxQ5yw7+lt8zh06NBX1m3nzp3FCfsBQkJC0LRpUxgYGMDCwgI9e/ZEcnKyXJu8vDwEBgaiWrVq0NfXR+/evZGRkSFS4ndXlnls3779K+tzzJgxIiV+PytWrICLiwsMDQ1haGgIT09P/P7777Lpyr4egbfPoyqsx9IsXLgQEokEEydOlI1ThfVZGS1fvhy1a9eGrq4umjdvjujoaLEjVQh16fPKQh3eS99Hae9Dqvqe+zr379/HoEGDUK1aNVSpUgXOzs6IjY2VTVeFz79lUbt27VfWu0QiQWBgIAD1e108ffoUEydOhK2tLapUqYKWLVsiJiZGNr28XhcsvMvAyckJaWlpsuHMmTNiR/pgz549g6urK5YvX17q9EWLFuH777/HypUrceHCBVStWhU+Pj7Iy8ur4KTv723zCACdO3eWW7dbt26twITlIyIiAoGBgTh//jzCw8NRWFgIb29vPHv2TNZm0qRJ2L9/P3bs2IGIiAg8ePAAvXr1EjH1uynLPALAqFGj5NbnokWLREr8fmrWrImFCxciLi4OsbGx+Oijj9CjRw9cuXIFgPKvR+Dt8wgo/3r8r5iYGPz8889wcXGRG68K67Oy2b59O4KCgjB79mxcvHgRrq6u8PHxQWZmptjRFE5d+ryyUIf30nf1uvchQPXec1/nyZMnaNWqFbS0tPD777/j6tWrWLJkCUxMTGRtVOHzb1nExMTIrfPw8HAAQN++fWVt1OV1AQAjR45EeHg4Nm3ahKSkJHh7e8PLywv3798HUI6vi3K5GrgKmz17tuDq6ip2DIUCIOzevVt2u6SkRLC0tBQWL14sG5eVlSXo6OgIW7duFSHhh/vvPAqCIAwZMkTo0aOHKHkUKTMzUwAgRERECILw97rT0tISduzYIWtz7do1AYAQFRUlVswP8t95FARBaNeunfD555+LF0pBTExMhNWrV6vkenzp5TwKguqtx6dPnwr169cXwsPD5eZNldenmJo1ayYEBgbKbhcXFwvW1tZCSEiIiKkqnjr1eWWlDu+lr/O69yFBUL333DcJDg4WWrdu/drpqvj5t6w+//xzoW7dukJJSYkgCOr1unj+/LmgoaEhHDhwQG58kyZNhOnTp5fr64J7vMvgxo0bsLa2hp2dHfz9/XHnzh2xIylUamoq0tPT4eXlJRtnZGSE5s2bIyoqSsRk5e/UqVOwsLCAvb09xo4di0ePHokd6YNlZ2cDAExNTQEAcXFxKCwslFufDg4OqFWrltKuz//O40thYWEwMzNDo0aNMHXqVDx//lyMeOWiuLgY27Ztw7Nnz+Dp6amS6/G/8/iSKq3HwMBAdO3aVW69Aaq5XYqtoKAAcXFxcstUKpXCy8uLy/T/qWKf9zbq8F76Nq97H3pJld5z32Tfvn3w8PBA3759YWFhATc3N/zyyy+y6er0+fffCgoKsHnzZgwfPhwSiUQ2Xl1eF0VFRSguLoaurq7c+CpVquDMmTPl+rrQLJfEKqx58+ZYv3497O3tkZaWhrlz56JNmza4fPkyDAwMxI6nEOnp6QCA6tWry42vXr26bJoq6Ny5M3r16oU6derg5s2bmDZtGnx9fREVFQUNDQ2x472XkpISTJw4Ea1atUKjRo0A/L0+tbW1YWxsLNdWWddnafMIAH5+frC1tYW1tTUSExMRHByM5ORk7Nq1S8S07y4pKQmenp7Iy8uDvr4+du/eDUdHRyQkJKjMenzdPAKqsx4BYNu2bbh48aLc78ReUrXtsjL466+/UFxcXGrfdf36dZFSVR6q2Oe9iTq8l5bFm96HANV6z32bW7duYcWKFQgKCsK0adMQExODCRMmQFtbG0OGDFGbz7//tWfPHmRlZWHo0KGycer0ujAwMICnpyfmz5+Phg0bonr16ti6dSuioqJQr169cn1dsPB+C19fX9n/Li4uaN68OWxtbfHrr79ixIgRIiajDzVgwADZ/87OznBxcUHdunVx6tQpdOzYUcRk7y8wMBCXL19WifMQvM7r5nH06NGy/52dnWFlZYWOHTvi5s2bqFu3bkXHfG/29vZISEhAdnY2du7ciSFDhiAiIkLsWOXqdfPo6OioMuvx7t27+PzzzxEeHv7Kt+hEYlDFPu9N1OG99G3K8j6kKu+5ZVFSUgIPDw98/fXXAAA3NzdcvnwZK1euxJAhQ0ROJ541a9bA19cX1tbWsnHq9LoAgE2bNmH48OGoUaMGNDQ00KRJEwwcOBBxcXHl+jw81PwdGRsbo0GDBkhJSRE7isJYWloCwCtn+MzIyJBNU0V2dnYwMzNT2nU7btw4HDhwACdPnkTNmjVl4y0tLVFQUICsrCy59sq4Pl83j6Vp3rw5ACjd+tTW1ka9evXg7u6OkJAQuLq6YtmyZSq1Hl83j6VR1vUYFxeHzMxMNGnSBJqamtDU1ERERAS+//57aGpqonr16iqzPisLMzMzaGhoqF3f9b6Uvc97G3V4L32bt70PFRcXv3IfZX3PLQsrKyvZ0VUvNWzYUPYTUnX8/Hv79m0cO3YMI0eOfGM7VX5dAEDdunURERGB3Nxc3L17F9HR0SgsLISdnV25vi5YeL+j3Nxc3Lx5E1ZWVmJHUZg6derA0tISx48fl43LycnBhQsX5H6HqWru3buHR48eKd26FQQB48aNw+7du3HixAnUqVNHbrq7uzu0tLTk1mdycjLu3LmjNOvzbfNYmoSEBABQuvX5XyUlJcjPz1eJ9fg6L+exNMq6Hjt27IikpCQkJCTIBg8PD/j7+8v+V9X1KRZtbW24u7vLLdOSkhIcP36cy7QUytrnvS91eC/9r7e9D5X2EwNlfc8ti1atWr1yKdI//vgDtra2ANTz8++6detgYWGBrl27vrGdKr8u/q1q1aqwsrLCkydPcOTIEfTo0aN8XxflcDI4lTZ58mTh1KlTQmpqqnD27FnBy8tLMDMzEzIzM8WO9kGePn0qxMfHC/Hx8QIAITQ0VIiPjxdu374tCIIgLFy4UDA2Nhb27t0rJCYmCj169BDq1KkjvHjxQuTkZfemeXz69KnwxRdfCFFRUUJqaqpw7NgxoUmTJkL9+vWFvLw8saO/k7FjxwpGRkbCqVOnhLS0NNnw/PlzWZsxY8YItWrVEk6cOCHExsYKnp6egqenp4ip383b5jElJUWYN2+eEBsbK6Smpgp79+4V7OzshLZt24qc/N189dVXQkREhJCamiokJiYKX331lSCRSISjR48KgqD861EQ3jyPqrIeX+e/Z4lVhfVZ2Wzbtk3Q0dER1q9fL1y9elUYPXq0YGxsLKSnp4sdTeHUpc8rC3V4L31f/34fUvX33P+Kjo4WNDU1hQULFgg3btwQwsLCBD09PWHz5s2yNqrw+besiouLhVq1agnBwcFy49XtdSEIgnD48GHh999/F27duiUcPXpUcHV1FZo3by4UFBQIglB+rwsW3m/Rv39/wcrKStDW1hZq1Kgh9O/fX0hJSRE71gc7efKkAOCVYciQIYIg/H1JhZkzZwrVq1cXdHR0hI4dOwrJycnihn5Hb5rH58+fC97e3oK5ubmgpaUl2NraCqNGjVLKD2elzSMAYd26dbI2L168ED777DPBxMRE0NPTEz755BMhLS1NvNDv6G3zeOfOHaFt27aCqampoKOjI9SrV0/48ssvhezsbHGDv6Phw4cLtra2gra2tmBubi507NhR9kFREJR/PQrCm+dRVdbj6/y38FaF9VkZ/fDDD0KtWrUEbW1toVmzZsL58+fFjlQh1KXPKwt1eC99X/9+H1L199zS7N+/X2jUqJGgo6MjODg4CKtWrZKbrgqff8vqyJEjAoBX5k8dXxfbt28X7OzsBG1tbcHS0lIIDAwUsrKyZNPL63UhEQRBeNfd8ERERERERERUNvyNNxEREREREZECsfAmIiIiIiIiUiAW3kREREREREQKxMKbiIiIiIiISIFYeBMREREREREpEAtvIiIiIiIiIgVi4U1ERERERESkQCy8iYiIiIiIiBSIhTcRERERERGRArHwJqL3Nnv2bDg7O6Nq1aqoXr06xo4di8LCQrFjERERqR1BEBAaGoo6depAT08PPXv2RHZ2NgCgffv2kEgkkEgkSEhIkLtffn4+xo0bh/HjxyM/P7/Uxx46dKjs/nv27FHwnBCpJhbeRPReBEGAIAj4+eefcfXqVaxfvx6//fYbVq9eLXY0IiIitfPll19ixYoV2LBhA06fPo24uDjMmTNHNn3UqFFIS0tDo0aN5O63ZcsWeHl5oVOnTti0aVOpj71s2TKkpaUpMj6RymPhTUTvRSKRYN68eWjZsiVsbW3h6+sLLy8vJCcnAwAOHDgAe3t71K9fn8U4ERGRAl24cAGhoaHYvn072rZtC3d3d4waNQqHDh2StdHT04OlpSU0NTXl7ltSUoKSkhIAf3+pXhojIyNYWloqbgaI1AALbyJ6L7dv30ZgYCAaNWoEExMT6Ovr49dff0XNmjVRVFSEoKAgnDhxAvHx8Vi8eDEePXokdmQiIiKV9O2336Jjx45o0qSJbFz16tXx119/vfW+/v7+OHHiBMLDwzF48GBFxiRSayy8ieidPXz4EE2bNsWjR48QGhqKM2fO4Ny5c5BKpXB1dUV0dDScnJxQo0YN6Ovrw9fXF0ePHhU7NhERkcrJz8/HwYMH8cknn8iNz8vLg5GR0Vvvr6urix9//BE//PADdHV1FRWTSO1pvr0JEZG8/fv3o7i4GFu3boVEIgEA/PjjjygsLETjxo0RERGBGjVqyNrXqFED9+/fFysuERGRyrp48SJevHiByZMnY8qUKbLxhYWF6NChg4jJiOjfWHgT0TurVq0acnJysG/fPjg6OmL//v0ICQlBjRo1YG5uLnY8IiIitfHHH3+gatWqr5ytvGvXrmjVqpU4oYjoFSy8ieiddevWDSNGjMDgwYNRpUoVDBo0CP369cPt27cBANbW1nJ7uO/fv49mzZqJFZeIiEhl5eTkwMzMDPXq1ZONu337Nm7cuIHevXuLmIyI/o2FNxG9M6lUipUrV2LlypWlTm/WrBkuX76M+/fvw8jICL///jtmzpxZwSmJiIhUn5mZGbKzsyEIguznXwsWLECXLl3g6OgocjoieomFNxGVO01NTSxZsgQdOnRASUkJpkyZgmrVqokdi4iISOV89NFHyMvLw8KFCzFgwACEhYVh//79iI6OFjsaEf0Lz2pORArRvXt3/PHHH0hJScHo0aPFjkNERKSSqlevjvXr12PFihVwcnLC+fPncebMGdjY2IgdjYj+hXu8iYiIiIiUWP/+/dG/f3+xYxDRG3CPNxERERGRivvpp5+gr6+PpKSkd77vmDFjoK+vr4BUROpDIgiCIHYIIiIiIiJSjPv37+PFixcAgFq1akFbW/ud7p+ZmYmcnBwAgJWVFapWrVruGYlUHQtvIiIiIiIiIgXioeZERERERERECsTCm4iIiIiIiEiBWHgTERERERERKRALbyIiIiIiIiIFYuFNREREREREpEAsvImIiIiIiIgUiIU3ERERERERkQKx8CYiIiIiIiJSIBbeRERERERERAr0f9O+wNJeE3DHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1000x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10,5))\n",
    "\n",
    "axes[0].plot(a0lst_th, gf_a0_th,'k-',label='Theory')\n",
    "axes[0].plot(a0lst_sim, gf_a0_sim,'.',color='royalblue',label='OSIRIS')\n",
    "axes[0].legend()\n",
    "axes[0].set_xlabel(r'$a_0$')\n",
    "axes[0].set_ylabel(r'$\\gamma_f$')\n",
    "#axes[0].set_title(r'PW: $\\gamma_e=2000,\\lambda=1\\mu \\mathrm{m},\\tau_\\mathrm{rise}=50~\\omega_p^{-1}$')\n",
    "axes[0].set_xticks([5,10,15,20,25,30,35,40])\n",
    "axes[0].set_xlim([5,40])\n",
    "axes[0].set_yticks([750,1000,1250,1500,1750,2000])\n",
    "axes[0].set_ylim([750,2000])\n",
    "\n",
    "\n",
    "axes[1].plot(thlst_th, gf_th_th,'k-',label='Theory')\n",
    "axes[1].plot(thlst_sim, gf_th_sim,'.',color='royalblue',label='OSIRIS')\n",
    "axes[1].legend(loc='upper left')\n",
    "axes[1].set_xlabel(r'$\\theta[^\\circ]$')\n",
    "axes[1].set_ylabel(r'$\\gamma_f$')\n",
    "#axes[1].set_title(r'PW: $\\gamma_e=2000,a_0=20,\\lambda=1\\mu \\mathrm{m},\\tau_\\mathrm{rise}=50~\\omega_p^{-1}$')\n",
    "axes[1].set_xticks([0,15,30,45,60,75,90])\n",
    "axes[1].set_xlim([0,90])\n",
    "axes[1].set_yticks([1250,1500,1750,2000])\n",
    "\n",
    "# aspect\n",
    "for i in range(2):\n",
    "    xmin, xmax = axes[i].get_xlim()\n",
    "    ymin, ymax = axes[i].get_ylim()\n",
    "    axesratio = (ymax-ymin)/(xmax-xmin)\n",
    "    axes[i].set_aspect(aspect=3/4/axesratio)\n",
    "\n",
    "plt.tight_layout()\n",
    "#plt.savefig('PW.pdf',transparent=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cfde2c3-3aa0-4f67-b4c3-9e5d61cedddf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}