fig6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PhysRevA.96.022128 - Figure 6\n",
    "\n",
    "Positron number scaling law\n",
    "\n",
    "Dots - data from paper (extracted with WebPlotDigitizer)\n",
    "\n",
    "Black line is equation 24\n",
    "\n",
    "Blue line is the integral $P_\\pm ~dN_\\gamma/dw$ (without saddle-point approximation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import API\n",
    "import pwpp\n",
    "from pwpp import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'def fun(E0): return integrate.quad(lambda w: ppm(0.855*sqrt(I21*1e21/1e18)*lbd,2*pi*c/(lbd*1e-6)*hbar/e*1e-9,tau*c/(lbd*1e-6)*1e-15,w) * dNgdw(E0/m,0.855*sqrt(I21*1e21/1e18)*lbd,2*pi*c/(lbd*1e-6)*hbar/e*1e-9,tau*c/(lbd*1e-6)*1e-15,w),0.01,E0)[0];\\nfun = np.vectorize(fun)\\nfull_E0 = fun(E0lst)'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xdim = 10\n",
    "\n",
    "# tau\n",
    "E0 = 1; #[GeV]\n",
    "I21 = 5; #[]\n",
    "lbd = 1; #[\\mu m]\n",
    "taulst = np.logspace(log10(5),log10(100),xdim);\n",
    "def fun(tau): return Np(E0/m, 0.855*sqrt(I21*1e21/1e18)*lbd, 2*pi*c/(lbd*1e-6)*hbar/e*1e-9, tau*c/(lbd*1e-6)*1e-15);\n",
    "fun = np.vectorize(fun)\n",
    "sadle_tau = fun(taulst)\n",
    "\"\"\"def fun(n): return integrate.quad(lambda w: ppm(0.855*sqrt(I21*1e21/1e18)*lbd,2*pi*c/(lbd*1e-6)*hbar/e*1e-9,n,w) * dNgdw(E0/m,0.855*sqrt(I21*1e21/1e18)*lbd,2*pi*c/(lbd*1e-6)*hbar/e*1e-9,n,w), 0.01,E0)[0]\n",
    "fun = np.vectorize(fun)\n",
    "full_tau = fun(taulst*c/(lbd*1e-6)*1e-15)\"\"\"\n",
    "\n",
    "# lbd\n",
    "E0 = 1; #[GeV]\n",
    "I21 = 5; #[]\n",
    "tau = 30; #[fs]\n",
    "lbdlst = np.logspace(log10(0.1),log10(20),xdim);\n",
    "def fun(lbd): return Np(E0/m, 0.855*sqrt(I21*1e21/1e18)*lbd, 2*pi*c/(lbd*1e-6)*hbar/e*1e-9, tau*c/(lbd*1e-6)*1e-15);\n",
    "fun = np.vectorize(fun)\n",
    "saddle_lbd = fun(lbdlst)\n",
    "\"\"\"def fun(lbd): return integrate.quad(lambda w: ppm(0.855*sqrt(I21*1e21/1e18)*lbd,2*pi*c/(lbd*1e-6)*hbar/e*1e-9,tau*c/(lbd*1e-6)*1e-15,w) * dNgdw(E0/m,0.855*sqrt(I21*1e21/1e18)*lbd,2*pi*c/(lbd*1e-6)*hbar/e*1e-9,tau*c/(lbd*1e-6)*1e-15,w),0.01,E0)[0];\n",
    "fun = np.vectorize(fun)\n",
    "full_lbd = fun(lbdlst)\"\"\"\n",
    "\n",
    "# I21\n",
    "E0 = 1; #[GeV]\n",
    "tau = 30; #[fs]\n",
    "lbd = 1; #[\\mu m]\n",
    "I21lst = np.logspace(log10(1),log10(40),xdim);\n",
    "def fun(I21): return Np(E0/m, 0.855*sqrt(I21*1e21/1e18)*lbd, 2*pi*c/(lbd*1e-6)*hbar/e*1e-9, tau*c/(lbd*1e-6)*1e-15);\n",
    "fun = np.vectorize(fun)\n",
    "saddle_I21 = fun(I21lst)\n",
    "\"\"\"def fun(I21): return integrate.quad(lambda w: ppm(0.855*sqrt(I21*1e21/1e18)*lbd,2*pi*c/(lbd*1e-6)*hbar/e*1e-9,tau*c/(lbd*1e-6)*1e-15,w) * dNgdw(E0/m,0.855*sqrt(I21*1e21/1e18)*lbd,2*pi*c/(lbd*1e-6)*hbar/e*1e-9,tau*c/(lbd*1e-6)*1e-15,w),0.01,E0)[0];\n",
    "fun = np.vectorize(fun)\n",
    "full_I21 = fun(I21lst)\"\"\"\n",
    "\n",
    "# E0\n",
    "I21 = 5;\n",
    "tau = 30; #[fs]\n",
    "lbd = 1; #[\\mu m]\n",
    "E0lst = np.logspace(log10(0.5),log10(4),xdim);\n",
    "def fun(E0): return Np(E0/m, 0.855*sqrt(I21*1e21/1e18)*lbd, 2*pi*c/(lbd*1e-6)*hbar/e*1e-9, tau*c/(lbd*1e-6)*1e-15);\n",
    "fun = np.vectorize(fun)\n",
    "saddle_E0 = fun(E0lst)\n",
    "\"\"\"def fun(E0): return integrate.quad(lambda w: ppm(0.855*sqrt(I21*1e21/1e18)*lbd,2*pi*c/(lbd*1e-6)*hbar/e*1e-9,tau*c/(lbd*1e-6)*1e-15,w) * dNgdw(E0/m,0.855*sqrt(I21*1e21/1e18)*lbd,2*pi*c/(lbd*1e-6)*hbar/e*1e-9,tau*c/(lbd*1e-6)*1e-15,w),0.01,E0)[0];\n",
    "fun = np.vectorize(fun)\n",
    "full_E0 = fun(E0lst)\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAI4CAYAAABndZP2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAByZklEQVR4nO3dd3hUdfbH8fdJQgABRRAb6mJb146ICqyFJiLKqtgbYAF7RV3soi7YxY4oCiIsFhQbyo9iEGTQpVuwIkoRUYqAkEJyfn/MJIaYQMrM3Cmf1/PMM5mZO/eeTCZnznzvt5i7IyIiIpJKMoIOQERERCTaVOCIiIhIylGBIyIiIilHBY6IiIikHBU4IiIiknJU4IiIiEjKyQo6gJows4XAWqAQ2OjuLYONSESSiZk1BJ4HDgAcuNDdQ4EGJSJRkdQFTkQ7d/8t6CBEJCk9Bnzg7qeZWTawVdABiUh0pEKBIyJSZWa2DXA00BPA3fOB/CBjEpHoSfYCx4H/MzMHnnX3waUfNLPeQG+AevXqHfqPf/wjgBBFUs/MmTN/c/cmQcdRQ7sDvwIvmtnBwEzgGnf/o/RGyiMi0RePHGLJvFSDmTV19yVmtj0wHrjK3T8qb9uWLVv6jBkz4hugSIoys5nJ3ufNzFoC04F/uvsnZvYYsMbdb6/oOcojItERjxyS1KOo3H1J5Ho58CZweLARiUgSWQwsdvdPIrdfB1oEGI+IRFHSFjhmVs/MGhT/DHQCPg82KhFJFu6+DFhkZvtE7uoAfBlgSCISRcncB2cH4E0zg/DvMdLdPwg2JBFJMlcBIyIjqBYAFwQcj4hESdIWOO6+ADg46DhEJHm5+xwgqfsSiUj5kvYUlYiIiEhFVOCIiIhIylGBIyIiIilHBY6IiIikHBU4IiIiknJU4IiIiEjKUYEjIiIiKUcFjoiIiKQcFTgiIiKSclTgiIiISMpRgSMiIiIpRwWOiIiIpBwVOCIiIpJyVOCIiIhIylGBIyIiIilHBY6IiIikHBU4IiIiknJU4IiIiEjKUYEjIiIiKUcFjoiIiKQcFTgiIiKSclTgiIiISMpRgSMiIiIpRwWOiIiIpBwVOCIiIpJyVOCIiIhIylGBIyIiIilHBY6IiIikHBU4IiIiknKygg5ARCQoZrYQWAsUAhvdvWWwEYlItKjAEZF0187dfws6CBGJLp2iEhERkZSjAkdE0pkD/2dmM82sd3kbmFlvM5thZjN+/fXXOIcnItWlAkdE0tmR7t4COB64wsyOLruBuw9295bu3rJJkybxj1BEqkUFjoikLXdfErleDrwJHB5sRCISLSpwRCQtmVk9M2tQ/DPQCfg82KhEJFqSfhSVmWUCM4Al7n5i0PGISNLYAXjTzCCcC0e6+wfBhiQi0ZL0BQ5wDTAf2DroQEQkebj7AuDgoOMQkdhI6lNUZrYLcALwfNCxiIiISOJI6gIHGAjcBBQFHIeIiIgkkKQtcMzsRGC5u8/czDaav0JERCQNJW2BA/wT+FdkLZlRQHsze7n0Bpq/QkREJD0lbYHj7je7+y7u3gw4C5jk7ucFHJaIiIgkgKQtcEREREQqkgrDxHH3HCAn4DBEREQkQagFR0RERFKOChwRERFJOSpwREREJOWowBEREZGUowJHREREUo4KHBEREUk5KnBEREQk5ajAERERkZSjAkdERERSjgocERERSTkqcERERCTlqMARERGRlKMCR0RERFKOChwRERFJOSpwREREJOWowBEREZGUowJHREREUo4KHBEREUk5KnBEREQk5ajAERERkZSjAkdERERSjgocERERSTkqcEQkrZlZppnNNrN3g45FRKJHBY6IpLtrgPlBByEi0aUCR0TSlpntApwAPB90LCISXYEWOGaWYWbdzOzUIOOQ1FOwchYbvn2agpWzgg5FYigKOWQgcBNQtJlj9DazGWY249dff63mYUQk3rICPv6zwA7AfsBoMzsA+Ie7vx5sWJJsClbOYuOK6WQ1bgXA2mnnQlE+ZGTToM0IajVqEXCEEiPVziFmdiKw3N1nmlnbirZz98HAYICWLVt6NIIWkdgLusBp4+77m9nsyO2vgRcAFThSaQUrZ21S0GTvemr4Z4qgqICNK6arwEldNckh/wT+ZWZdgDrA1mb2srufF6NYRSSOgu6Ds6r0DXcvIPiiSxJIZU41bVwxfZOCBhwysoFMyKhV0qojKanaOcTdb3b3Xdy9GXAWMEnFjUjqCLqYmGlmXQEHMLP6QINgQ5JEUbZlpqJTTVmNW4ULmqICyKhF7V1Ppfaup5acstpc603pU1tq5UlKyiEiUq6gC5xbgf8Ce5jZw0BH4L1gQ5IglS44yrbMVHSqqVajFjRoM+IvhcqWCpbyCihABU9yiUoOcfccICeqkYlIoAItcNx9HdDVzI4GWgL/AV4LMiYJTtmCY6sD7tikZWZzp5pqNWpR5YKkbAGVt2g0+YveUOfkJKIcIiIViVmBY2bXAbOBOe6+enPbuvtHwEexikWSQ9mCwwtWldsyEy1lT22BqXNyAlEOEZGaiGULzsP8eV78JyKJKnI9290Xx/DYkoTKFhzFRU2sioyyp7YA8heNrlSLkcSFcoiIVFssC5ydgEMil+aR65MAA9zMVvJnsprj7iNjGIskgYr60sT6mKWPE+/jy2Yph4hItZl7/OatioxwaM6fyao5cACQ6e4x7Q/UsmVLnzFjRiwPIWlAo67CzGymu7cM4LiB5RBQHhGJlnjkkLh2Mo50CJxqZquAbYCDgVpAbjzjEKmOyg5bl9hRDhGRyorbRH9m1sbMHjCzb4B5wPXAl8DpQJNq7K+OmX1qZnPN7Asz6xflkEU2Ud6wdYmfaOcQEUltsRxFlU14ToqTga6E14v5CXgLGAN85O6FNThEHtDe3deZWS3C3+red3d96iSgVDi1U14naImdOOQQEUlhsTxF9StQH5gLDALecvc50dq5hzsPrYvcrBW5aCG8BJQqp3aC6ASd5mKaQ0QktcWywGlAuJVlDbAtcLCZFQFfuvvGaBzAzDKBmcBewFPu/kmZx3sDvQF22223aBxSqqGyMxIng80NW0+FVqoEE/McIiKpK5YFTi/+HOVwIeFvYg4UmNl8/pzTYg7hIZ5rqnqASPN0czNrCLxpZge4++elHh8MDIbw6Ica/C5SA+lwaidVWqkSTMxziIikrpgVOO4+pPRtM9ubTYd3dgZ6Fm8OZNbgWKvN7MPIPj/f0vYSX+lwaieVWqkSRTxziIiknlh2Mm7j7tOKb7v7t8C3lFonxsy2B1oQHupZ1f03AQoixU1d4Fjg/hoHLjERyxmJE0E6tFLFW6xziIiktlieoppiZsuBd4A3gYnunl96A3dfDnwQuVTVTsCwSD+cDOBVd3+3hjGLVEs6tFIFINY5RERSWCwLnKaEh3eeRDg55ZnZuMjP79X0fLm7zyPcTC2SEFK9lSoAMc0hIpLaYjbRn7svc/dB7n484Um4LgEKgWeAX83s/8zsMjPbOVYxiCSCgpWz2PDt0xSsnBV0KElFOUREaiIuMxm7+1p3H+XuZxNOVCcB3wO3AYsiMxLfHI9YROKpeHTVhvkPs3bauSpyqkk5RESqKm5LNRRz9wJ3/8DdL3P3psA/gUnA+fGORWpGLRNbpuUdok85REQqI66LbZYnsrTCdKBv0LFI5Wnel8rR6KrYUw4RkfIEXuBIctK8L5Wj0VUiIsFQgSPVopaJytPoKhGR+Au8wDGzvSMTeJV7WxKTWiYkUSiHSDLTGnaxE3iBA7xCeCbSim5LglLLhCSIpM8hlf2QS/YPw/Lir+rvVHb7qjw/0Y5flb6Myf63D0IiFDi2hdsiIpsTtxziG9ez4duno/4BWZkPuYq2i+cHdE2eX178QJUGK5Tdx1YH3MH6z++udIGQaMevbF/GaP3t001gBY6Z7UY4EdUys10Bc/efCC+aJyKyWUHkkMI/fmDD/Iej9gFZlQ+5iqYciNcHdE2fX+GUCVUYrFB2H/k/v1/p5yfi8SvblzEaf/t0FGQLzrDI9W7AS4STUnvUgiMilRP/HOJFRPMDsiofcuVtV9XRjDX5gK7p8yv8PaswWKHsPrJ3Op6NK/5Xqecn4vEr25cxGn/7dBRYgePu7QDMbHbxz8UPBRSSSGDU1Fx1geQQywAyo/YBWZUPuQq3i9MHdE2fX1H8VRmsUN4+Mrfep1LPT9TjV6YvYzT+9unI3IOtJ8xslru3qOh2tLRs2dJnzJgR7d2K1FgyTppoZjPdvWXQcUD8cgjAoc3386mvXRnVPjg1lSx9cCT6kvm1j0cOSYROxt9s4bZISlNTc41VK4eYWR3gI6A24Vz4urvfudnnZG1F3b0v3+S+qo4mjPbow5oeP97Pl+jRa795gRc47n7W5m5L8JL5W0Iy0KSJNVODHJIHtHf3dWZWC5hqZu9Hln4QkSQXeIEjia0yp09CoRA5OTm0bduW1q1bBxRp8tKkicHw8Pn5dZGbtSIX9QEUSREqcGSztnT6JBQK0aFDB/Lz88nOzmbixIkqcqpBTc3BMLNMYCawF/CUu39Szja9gd4Au+22W3wDFKkmtbxDRqx2bGbXmVlbM2sYq2NI9BWsnMWGb5+mYOUsoNTpk7IjRyJycnLIz8+nsLCQ/Px8cnJyNrv/UCjEgAEDCIVCMfoNJFXEI4e4e6G7Nwd2AQ43swPK2Wawu7d095ZNmjSJVSgiUVPc8r5h/sOsnXZuST5PN7FswXmYSHOvmf0EzAbmRK5nu/viGB5bqqH4n8KL8tlYZCxtdDuHHN1js6dP2rZtS3Z2dkkLTtu2bSvcv1p7pIrilkPcfbWZfQh0Bj6P1n5FgqCBC2GxLHB2Ag6JXJpHrk8iPAmXm9lK/kxWc9x9ZAxjkUrYuGI6XpSPUQRFzshn+pBb6++0bt26wn+O1q1bM3HixEr1wSmvtUcFjmxGTHOImTUBCiLFTV3gWOD+qEUvEhANXAiLWYHj7r8AH0QuAJhZfcKJqjl/Jq1rgExABU7Ashq3YmORQZFTUOhMmbueRpUoQlq3bl2pQqUqrT0iccghOwHDIv1wMoBX3f3dmsYtEjQNXAiLaydjd19HeCjmKmAb4GDCIxdy4xmHlK9WoxYsbXQ7I5/pw5S565m3EB6NYhFSldYekfJEM4e4+zzCRZJIytHAhTgWOGbWBjg5ctkTWA28Bwyg1Dc0CdYhR/cgt9bfaZSTw6OVLELWr1/PsmXL/nL55ZdfNvn5iCOO4JFHHuHmm2+Ow2+SmtJ5ZIRyiIhURcwKHDPLBjoSTkZdgR2An4C3gDHAR+5eGKvjS/W1bt2aww47jOXLlzNr1qxyi5fSl7Vr1/5lH2ZGkyZN2HHHHdlxxx1p1qwZb7/9Nu+//z79+/fnsssuIzMzc4uxaI6dPyXjkg41oRwiIjURyxacX4H6wFxgEPCWu8+J4fEkCtavX8/VV1/Niy++SFFR0V8e32abbUqKlhYtWpT8XPqyww470KRJE7KyNn17ff/991x++eVcddVVDB8+nGeffZbmzZtXGItGXW0qDUdGKIeISLXFssBpQHgq9DXAtsDBZlYEfOnuG2N4XClHZVpCvvrqK04//XS++OILLr30Ug466KC/FC5169atdgx77rknH3zwAf/973+57rrraNmyJddddx133XUX9erV+8v2GnW1qTQcGaEcIiLVFssCpxd/jnK4kPA3MQcKzGw+f85pMYfwEM81MYwlrZVuCWm1fz2GPHIZexxyxibf/keMGMEll1zCVlttxbhx4zj22GNjEouZcc4559C5c2f69u3LQw89xGuvvcbTTz9Nly5dNtlWo642lYYjI5RDRKTaLLwcSxwOZLY3fx3euWPkYXf3LXfIqIGWLVv6jBkzYnmIhDVgwABuv/12DtkrizH37EDt7EwyMmuHPyzr7ss111zDc889x1FHHcWoUaPYeeed4xbb1KlTueSSS/jyyy85/fTTeeyxx9hpp51KHlcfnMRkZjPdvWWcjxloDoH0ziMi0RSPHBK3UVTu/i3wLfBa8X1mtj3QgvBQT4mR4paQow6qTa0sI8McigqYP+Mtut90IXPnzuXmm2/m7rvv/ku/mVg78sgjmT17Ng8++CD33HMP48aN47777uOSSy4hIyOj0nPsSOpTDhGJnnQYkRm3Fpygpfs3r1AoxPxPXuFfe7xPBoWM+XgD1z79O9m16zB8+HCOP/74oEPk22+/5bLLLmPixIm0atWKwYMHc+CBBwYdlpQjiBacRJDueURSQyKMyIxHDonZYpvFzKyZmXWLNC9LQFq3bs2F1w6kzmHDuPmVplz0wFIOPOhgZs+enRDFDcDee+/N+PHjGT58ON999x0tWrSgb9++rF+/PujQJEDKISLRVd6IzFQU9QLHzMZEmo0xsxOAT4BLgclmdla0jyeVt2DBAo458QoGj/yIG264gZycHHbdddegw9qEmXHeeefx1Vdf0b17d+6//34OOOAAxo0bF3RoEifKISKxVTIik8yUHpEZixac3d19eeTnvsDR7t4JODRyWwLw5ptv0qJFC3744QfefvttHnzwQWrVqhV0WBVq3LgxQ4YMIScnh+zsbDp37sw555zDsmXLNvu8UCjEgAEDCIVCcYpUYkA5RCSGikdk1t33+pSeMDQWBU7piVJqu/vXAO7+c4yOJ5uRn5/PtddeS7du3dhnn32YNWsWXbt2DTqsSjvmmGOYO3cud911F6NHj2bfffdl8ODB5U5CWDwc/vbbb6dDhw4qcpKXcohIjNVq1IK6e1+essUNxCZZvGdmT0ZW/R1vZj3MLMvMugLLt/RkqZzKtFQsXLiQo446iscee4xrrrmGKVOm0KxZs/gFGSW1a9fmzjvvZN68eTRv3pxLLrmEo48+mi+++GKT7cqbGDAdFaycxYZvn6Zg5aygQ6ku5RARqbFYjAm+ifDid0uAVcBuwBBgHHBRDI6XdiqzhMHbb79Njx49cHdGjx5Nt27dAoo2evbZZx8mTZrESy+9RJ8+fWjevDnDhw/nrLPC3TI0MWBijI6IAuUQEamxqLfguHuBu99AeAKuroTPm+/g7ie4+4/ROo6Z7WpmH5rZl2b2hZldE619J7rNtVSs/+UTrru4IyeddBJ77LEHs2bNSonippiZ0aNHD7766ivatGnD+eefz3vvvQeER4pNnDiRe+65J23XrUqF0RHxyiEiktpiNqubu28APovV/oGNQB93n2VmDYCZZjbe3b+M4TETQkUtFQs+G8vZZ57Gp/M3cGGXhgx87kka7LxHsMHGyHbbbcc777xD+/btOe200xg3bhxHH3102k8MmErrVcUhh4hICotphz0zq2tm10ZaWn4xs/zI5ZfIfdea2VbV2be7/+zusyI/rwXmA02jGX+iKq+l4v333+fwo8/gy4W5PH9DEx68pDFZf8wOOtSY2nrrrfnggw/YfffdOfHEE5k5c2bQIQUu1UZHxDKHiEhqi9lMxma2KzAJaAZ8DHwBrIw83AjYD/gn8CPQwd1/qsGxmgEfAQeUXnDPzHoDvQF22223Q3/8MfVat4uKirj99tvp378/B+6/N89fVcheOxlk1EqJD7jKWLJkCUceeSTr1q3jo48+Yt999w06pJQXj1lI45lDKkszGYtER7KvRTUQ2ADs7e4Ly9sgUpiMAR4FTq3OQSIjLUYD15ZdTdjdBwODIZyYqrP/RObuXHHFFQwaNIiLL76Yxx9/nKwN81N+fZGymjZtyvjx4znyyCM59thjmTp1alKOFpO/GEgccoiIpKZYnqLqCNxaUWICiDx2R2TbKjOzWoSLmxHu/kZ19pHMbrvtNgYNGsS///1vnnvuOerWrZsWcxuUZ6+99mL8+PGsX7+ejh07bnFCQEkKMc8hIrKpFJhmokQsC5yqtJhUuXXFzIzw0NH57v5IVZ+f7B566CH69+9P7969GTBgQNDhJIQDDzyQsWPHsmzZMjp16sSqVauCDklqJqY5REQ2VTzNxIb5D7N22rlJX+TEssCZAPzHzHavaINI8/I9wPhq7P+fwPlAezObE7l0qVakSWbIkCHceOONnHHGGTz99NOEaz0BaNWqFWPGjOHrr7+mS5curFu3LuiQpPpinUNEpJRUmGaitFj2wbkW+BD4xsymA58TnrQLYFtgf6AVsBC4rqo7d/epQNp9sr/++uv07t2bzp07M3z4cDIzM4MOKeF07NiRUaNGcdppp3HKKafw7rvvUrt27b9sFwqFyMnJoW3btmk9tDyBXUsMc4iIbCqVppmA2M6Ds9jMDiI8iqkrcDLhkQ8QTlJfADcCz7n7+ljFkUrGjx/POeecQ+vWrRk9ejTZ2dlBh5SwTjnlFF544QV69uzJ2WefzauvvkpW1p9v98rMBp3qClbOSugO6cohIvFVPM1EIueFqohlC07xRF2PRS6yGVtqTQiFQpx88snst99+vPvuu2y1lab+2JIePXrw+++/c80113DxxRfzwgsvkJERPitb3mzQ6VTgJMuSDsohIvFVq1GLhMwF1RHTAkcqZ0utCfPmzaNLly7svPPOjBs3joYNGwYXbJK5+uqrWb16NXfeeSfbbLMNAwcOxMzSft2q8s61p0pSExGBGM9kXBlmdrSZTQo6jiBtbm2p77//nuOOO4569eoxfvx4dthhh+ACTVK333471157LY8//jj9+vUDtG5Vybl2MkvOtSfr8FDlEBEpTyK04DQBjgk6iCBV1JqwdOlSjj32WAryN/Dey1fRdOuVhCd1laowMx5++GF+//13+vXrR8OGDbn22mvTet2qsufagaQ4ZVWBtM8hIvJXMStwzGy3Sm7aJFYxJIvi1oTSfXBWrFjBsccey6/Lf+HNu5vQrGAka6e9nmwfPAkjIyODwYMHs2bNGq677jq22WYbLrjggqDDClTpc+0bvn064U5ZKYeISE3EsgVnIZWbfMsquV1KK92asG7dOrp06cL333/PmOcupcU275BIHzzJKisrixEjRrB27VouvvhittlmG7p16xZ0WAkhQYeHLiSGOSSy1tVLwA6R5w92d3VmFkkRsSxwNhBeAPP1LWzXksiCmAK5ubmcfPLJzJw5k9GjR9PhqF1ZO21con3wJK3atWvzxhtv0KlTJ8466yzeffddOnXqFHRYgUvQ4aGxziEbgT7uPsvMGgAzzWy8u39ZjX2JSIKJZYEzFyh09yGb28jMVqMCB4CNGzdy9tlnM3HiRIYNG8ZJJ50EkIgfPEmtXr16vPvuu7Rt25ZTTjmF8ePH06ZNm6DDClwCDg+NaQ5x95+BnyM/rzWz+UBTQAWOSAqI5SiqmcChldw27WYkLquoqIhevXoxZswYHnvsMbp3717yWLouoBlL2267Lf/3f//HzjvvzAknnMDcuXODDkn+Km45JLLkwyHAJ+U81tvMZpjZjF9//bUmhxGROIplgXMfcNaWNnL30e4e+HD1ILk7N9xwA0OHDuWuu+7i6quvDjqktLDDDjswYcIE6tevT6dOnfj222+DDkk2FZccYmb1gdHAte6+ppz9D3b3lu7eskkT9WeW9JOsU0jErLBw9yXuPjlW+08l//nPf3j00Ue5+uqrueOOO4IOJ6387W9/Y/z48RQVFdGxY0cWLVpU4bahUIgBAwYQCoXiGGH6ikcOMbNahIubEe7+RiyPJZKMknmFcXNPjwFMLVu29BkzZgQdxl889dRTXHnllXTv3p0XX3yxZCkBia9Zs2bRrl07dt55Z6ZOnUrjxo03eVxrV23KzGa6e8ug46gJMzNgGLDS3a+tzHNq167t+++/P1tvvXVUY1mzZg2rV6+mYcOGUd93rPefrPuO9f5TJfb6tdZRtH4R4YGGRsZWu5JRd+ca73/27NlL3H2XqARcgUSY6C9lbWl9qREjRnDllVfyr3/9iyFDhqi4CVCLFi1455136NSpE127dmXixInUrVu35PF0X7sqRf0TOB/4zMzmRO67xd3HVvSE/Px85s6dy8EHHxy1D5U1a9Ywd+5cioqKyMjIiOq+Y73/ZN13rPefSrEftP/e1Cvp4mZYrZodp3j/hDv0x5a7p8Xl0EMP9XiaNm2a161b1zMzM71u3bo+bdq0TR5/++23PTMz09u1a+cbNmyIa2xSsddee83NzE8++WTfuHFjyf1b+numG2CGJ8D/dbwvgGdmZnr//v2j8jq6u/fv398zMzM9FvuO9f6Tdd+x3n+qxZ6/Yqav/+Ypz18xM6r79xj/v6rJIEY2t77U5MmTOeOMM2h+4D6MevhEMtdrVGqiOO200xg4cCBjxozh6quvLv5QS/u1q+RP0V6ctXiplszMzJgs/BrL/SfrvmO9/1SLPZojeYv3Txwm+FUfnBipqM/GzJkzadeuHU132o6378ikcYOiZFz7J+XdeOONPPTQQwwYMIC+ffsGHU7CSYU+ONWxyy67+GuvvRb1AndLp7MTef/Juu9Y71+xb37/bdq0iXkfHBU4MVT2TfLFF19wzDHHUL9+fcYP70XjVS8ARUAmdfe9nrp7Xx7X+KRiRUVFnHvuuYwaNYqXXnqJ888/P+iQEkq6FjiJOlhBJNnEI4eok3EMlV5f6ptvvqFDhw5kZ2czYcIE/tZoDWunvawlGBJURkYGQ4cO5ZdffuHCCy9kp512omPHjkGHJSIilaQ+OHGwcOFCOnToQFFRERMnTmSvvfYqWfun7r7X6/RUgipet2rfffelW7duzJkzJ+iQRESkklTgxNiSJUto374969atY/z48ey7774lj2kJhsTXsGFDxo4dyzbbbMPxxx/Pjz/+GHRIIiJSCSpwYuiXX36hQ4cO/Pbbb/zf//0fBx98cNAhSTXssssufPDBB2zYsIHOnTuzcuXKoEMSEZEtUIETIytWrODYY49l0aJFjB07lsMOOyzokKQG9t9/f9566y0WLFjASSedRG5ubtAhiYgEJhnWp1KBEwOrV6+mU6dOfPPNN7z99tsceeSRQYckUXDMMccwfPhwpk6dynnnnUdhYWGF22rdKhFJVcmyPpVGUVVBZeYGWLduHV26dOGzzz5jzJgxdOjQIc5RSiydccYZLFmyhOuvv57rrruOxx57jPCSRn/SulUikso2rpgORflAERQVsHHF9ITsS6oCp5Iq86G1YcMGunbtyqeffsqrr75Kly5dAopWYum6665j0aJFPProo+y6667ceOONmzyudatEJJVlNW4FGdkJP82JCpxK2tKHVl5eHqeccgqTJ09m6KB7OP7AZRSsnJWQVa3U3EMPPcSSJUu46aabaNq0Keecc07JY8VTkRcXw9Gepl1EJEjF05xsXDGdrMatEvZzTgVOJW3uQ6ugoIAzzzyTcePGMfjxOzhx55fZMD9fSzCksIyMDIYNG8Yvv/xCz5492XHHHWnfvj3w57pVsZzqXEQkSLUatUj4zzYVOJVU0YdWYWEh559/Pm+99RZPPPEE5x2XES5uEvzcpNRcnTp1GDNmDEceeSSnnHIKU6ZM4aCDDgI2ncVaRETiTwVOFZT90CoqKuLiiy/mlVde4YEHHuDKK68M9yZPgnOTEh0NGzbk/fffp3Xr1hx//PGEQiF22223oMMSEUl7GiZeTe7OlVdeydChQ+nXr19JR1MtwZB+dt11V95//33WrVvH8ccfz6pVq4IOSUQk7anAqQZ3p0+fPjzzzDP06XUsfa84YZPHtQRD+jnwwAMZM2YM3377LSeffLImAhQRCZgKnGq44447ePTRR+l14rbcfMJ3rAudl7ATHUn8tGvXjpdeeomPPvqIHj16UFRUFHRIIiJpS31wqqh///7ce++99Dy9Df3P+RkzV2diKXHWWWexePFibrzxRpo2bcojjzwSdEgiImlJLThVMHDgQG699VbOO+88nnlqIBlZtYFMdSaWTfTp04err76aRx99tNIFjpZ2EBGJLrXgVNKgQYO47rrrOPXUU3nxxRfJysoiMwkmOpL4MzMeeeQRli5dSp8+fdhxxx03mQiwLC3tICISfUnbgmNmL5jZcjP7PNbHGvJUPy677DK6dDqSkSNHkpUVrgvVmVgqkpmZyfDhwznmmGM4//zzee655yrctrxZskVEpGaStsABhgKdY32QkS/cR++r7+KYg+vy3CXLsXUxr6ckRdSpU4exY8dy3HHH0bt3b+6///5ytyueJTszM1NLO4hI0ipYOYsN3z6dMINukvYUlbt/ZGbNYnmMt956ix6X3MoR/6jD8FuaUCerUJ2JpUq22morxowZQ8+ePenbty8rVqzg/vvv32QFci3tICLJrmDlLNZOOze8yniCLFOUtAVOZZhZb6A3UK3ZZb/55hsOOWhfRv47n3p1itSZWKolOzubl19+mW233ZYHH3yQFStW8Oyzz5ac6gQt7SAiyW3jiunh4iaBlilK6QLH3QcDgwFatmzpVX3+jTfeyDXXXIOt+1ydiaVGMjIyePLJJ9luu+24++67WbVqFSNHjqROnTpBhyYiUmNZjVsl3DJFKV3gREN2djYkwaqpkvjMjH79+tG4cWOuueYaTjjhBMaMGUODBg2CDk1EpEaKlylKpMYAFTgicXb11VfTqFEjevbsSfv27Xn//ffZbrvtgg5LRKRGaiVYY0DSjqIys/8CIWAfM1tsZhcFHZNIZZ133nm8+eabfP755xx11FEsWrSo0s/VpIAiIluWtC047n520DGI1ETXrl0ZN24cXbt25Z///Cfjx49nn3322exzNClg9JjZC8CJwHJ3PyDoeEQkupK2BUckFRx99NHk5OSQl5fHkUceycyZMze7vSYFjKqhxGEuLREJhgockYAdcsghTJ06lXr16tGuXbvNFi2aFDB63P0jYGXQcYhIbKjAEUkAe++9Nx9//DG77rornTt35q233ip3u+JJAe+55x6dnooTM+ttZjPMbMavv/4adDgiUkkqcEQSRNOmTfnoo49o3rw5p556KsOGDSt3u9atW3PzzTeruIkTdx/s7i3dvWWTJk2CDkdEKkkFjkgCady4MRMmTKBdu3b07NmTRx55JOiQRESSkgockQRTv3593n33XU477TT69OnDrbfeinuVJ+IWEQlMIiy8qQJHJAHVrl2bUaNG0atXL/r3789ll11GYWFhtfaleXPKp7m0RGKjeOHNDfMfZu20cwMrcpJ2HhyRVJeZmcmzzz5L48aNue+++1i1ahXDhw8PLx9SSZo3p2KaS0skNhJl4U214IgkMDNjwIABPPjgg7z66qt07dqVP/74o9LP17w5IhJvJQtvkhnowptqwRFJAjfccAONGjWiV69edOzYkTFjxrDDDjts8XnF8+YUt+Bo3hwRibVEWXhTBY5IkrjwwgvZdtttOfvss9lvv/14+OGH6dGjB2ZW4XOK583Jycmhbdu2fzk9FQqFKnxMRKS6EmHhTRU4IknklFNOYc6cOfTq1YsLLriAl19+mWeffZY999yzwue0bt263OJF/XNEJJWpD45IkvnHP/7B5MmTeeaZZ/j000858MADeeihh9i4cWOV9qP+OSKSylTgiCShjIwMLr30UubPn0+nTp248cYbOeKII5g9e3al96F1rUQklanAEUliTZs25c033+T1119n6dKlHHbYYfz73/9m/fr1W3yu1rUSkVSmAkckyZkZp556Kl9++SUXXHABDzzwAAcddBCTJk3a4nO1rpWIpCoVOCIpYtttt+W5555j0qRJmBkdOnTgoosuYuXKlUGHJiISdypwRFJMu3btmDdvHn379mXYsGHsu+++vPrqq1rPSkTSigockRRUt25dBgwYwIwZM9h1110588wzOemkk1i8eHHQoYlIGgpi8U0VOCIprHnz5kyfPp2HH36YCRMmsN9++/H0009TVFQUdGgikiaCWnxTBY5IisvKyuL666/n888/p1WrVlxxxRUcddRRfPnll0GHJiJpoLzFN+NBBY5Imthjjz0YN24cw4YN46uvvqJ58+b069ePvLy8oEMTkRQW1OKbKnBE0oiZ0b17d+bPn8/pp5/OXXfdRYsWLZg2bVrQoYlIiipefLPuvtfToM2IuK1RpQJHJA1tv/32jBgxgvfee4+1a9dy1FFHsWDBgqDDEpEUVatRC+rufXlcF+DUYpsiaaxLly588cUXvP/+++yxxx5BhyMiEjVqwRFJcw0aNOCMM84IOgwRkahSgSMiIiIpRwWOiIiIpBwVOCIiIpJyVOCIiIhIylGBIyIiIilHBY6IiIikHBU4IiIiknJU4IiIiEjKUYEjIiIiKUcFjoiIiKScpC5wzKyzmX1tZt+ZWd+g4xGR5KIcIpK6krbAMbNM4CngeGA/4Gwz2y/YqEQkWSiHiKS2pC1wgMOB79x9gbvnA6OAkwKOSUSSh3KISApL5gKnKbCo1O3FkftERCpDOUQkhWUFHUAsmVlvoHfkZq6ZfVGN3WwH/FaF7bcBfq/hNtV9vKqxBqkyr1OiHKe6+0jl987fqrBtUiuTR/LM7PM4Hbqm792qPr8y22/ufVLR88u7vzL3xTuf6fWO7+u9d8yP4O5JeQFaA+NK3b4ZuHkz2w+u5nFmVHH7LR5nS9tU9/Gqxhrw369af48gjqP3TmpeqppD4v061fS9W9XnV/L9V+Hvv5n31l/ur8x98X5P6vVOrte7MpdkPkX1P2BvM9vdzLKBs4C3N7P9O/EJq1LH2dI2NX08GSTS3yMe+4jWcfTeiZ6q5pB4q+nfqqrPj9Xxyru/svfFk17v+Ir58S1SSSUlM+sCDAQygRfc/T8xOMYMd28Z7f3GQjLFmg6S6e+RTLFGU1VzSLq+TsXi+fun+2sNer1rKqn74Lj7WGBsjA8zOMb7j6ZkijUdJNPfI5lijZpq5JC0fJ1Kiefvn+6vNej1rpGkbsERERERKU8y98ERERERKZcKHBEREUk5KnBEREQk5ajASWFmtoeZDTGz14OORZKLmZ1sZs+Z2Stm1inoeBKd/tdiz8zqmdmwyPvy3KDjSXWp8J5WgZNkzOwFM1tedjbV8lZF9vAaOxcFE6kkmiq+d8a4ey/gUuDMIOKNl6q8LhXR/1r1VPG17wa8Hnlf/ivuwaaAdPv8UIGTfIYCnUvfoVWRpZKGUvX3zm2Rx1PZUCr5upjZgWb2bpnL9vEPOWUMpfLvyV34c+2wwjjGmEqGkkafHypwkoy7fwSsLHO3VkWWLarKe8fC7gfed/dZ8Y41nqryurj7Z+5+YpnL8rgHnSKqmM8WEy5yQJ9d1ZJunx96k6SGcldFNrPGZjYIOMTMbg4mNElwFa2ofRXQETjNzC4NIrCAVWmlcf2vRVVFr/0bwKlm9gzBLzOQSlL28yOpZzKWzXP3FYT7UIhUibs/DjwedBzJQv9rsefufwAXBB1HukiF97RacFLDEmDXUrd3idwnsiV675RPr0tw9NrHV8q+3ipwUkOir4osiUvvnfLpdQmOXvv4StnXWwVOkjGz/wIhYB8zW2xmF7n7RuBKYBwwH3jV3b8IMk5JPHrvlE+vS3D02sdXur3eWmxTREREUo5acERERCTlqMARERGRlKMCR0RERFKOChwRERFJOSpwREREJOWowBEREZGUowJHREREUo4KHBEREUk5KnAk7sxsqJl55JJT5rELzexbM8s3s9WV3N9tpfa3OBYxi0j0RZYHcDNbbmYXVbBNhfkiHpRfkpcKHAnKMqA1cHnxHWa2MzAYmAa0BzpWcl8vRvY1Nsoxikhs/QYcC6wC7t7Mdn/JF3Gk/JKksoIOQNJWnrtPL3Pf3kAmMMzdp1Z2R+6+BFhiZr9GM0ARiS13XwtMMLMhwP1mtr27Ly9n0/LyRVwovyQvteBIVJnZLaWac8u7PF3B84YCOZGbEyPbDo089nczezPSjJ1rZj+Z2WtmpgJdJDV8Hbk+uLJPMLMMM1trZneWuX/bSP7oEbl9V+T2P8xsnJn9EckhF0QeP9/MvjKzdWb2oZntGbXfSgKlDwiJtlHApMjPZwDXAe2A3Mh9C4H7ynnePcBM4HHgCmAWUPyN6T3CTdiXEW7Sbgp0QQW6SNIzMwOuj9w8CBhfyaf+HahPOFeUdkjkenaZ+18DngMeInyq6wUz2xtoC/QFagGPASOBIyr/G0iiUoEjUeXuC4AFAGbWE1jo7jmltwnns78873szmx+5+WVxc7SZbQfsBZzk7m+XesrIqAcvIkG4nHAflz+oQgsO0CJyXV6Bkwd8Web+B939JQAzmwF0BS4Bdnf3NZH7dwIeM7O/ufuPVfotJOHoG7DE0kHAvBruYwXhguk+M+sV+cYlIinAzJoRbtH9D/AJ4ZxRWS2A5ZE+MmXv/8LdN5a5//3iH9x9FbAcmF5c3ER8FbnetQpxSIJSgSMxEWl2PgCYW5P9uLsTHmUxAxgAfGNmC8zssppHKSIBex74DugPzAH2NbNalXzuofy19QbCLThlT09B+DR3afkV3AdQp5IxSAJTgSOx8jegATVvwcHdF7h7d6AJ4eQ1CXjazI6v6b5FJBhm1hs4Gujp7gWEC5xs4B+V3EVzyhQyZrY9sE/Z+yU9qcCRWNk5cr0wWjv0sDn82SHxgGjtW0Tix8x2AR4E7nX34lbeOZHrLZ6miox0aggUlnnoKsKfa3OQtKdOxhIrf0SuTzOzrOrOYWFmBxEe2fAK4absTKAnsJE/R2uJSHJ5Dvie8KmpYvMJdw4+GBixhecfGrm+2MwWEe5P0wnoEbm/pZnNcvcN0QtZko1acCRW5gGDgN7AyzXYzzLgJ8KtNm8D/yXcOnSiu8+saZAiEl+R0ZUdCJ+aKukIHPn5cyrX0bgFsJLw8O5+wHDCp8RPB9YAZ6q4EQv34RSJn8gEfm0JD/92dy/bzFzV/Rnhlp0hQAd336WmMYpIYigvX5jZ/xH+/Do2DsdXfklSasGRoPwNKAAmRmFft0b21T0K+xKRxFM2X7QgPDFoPCi/JCm14EjcRea+2C5yc627f72ZzSuzv50Iz24MkO/uNR65JSKJoWy+IDwr+kLgDHd/LQ7HV35JUipwREREJOXoFJWIiIikHBU4IiIiknJU4IiIiEjKUYEjIiIiKUcFjoiIiKQcFTgiIiKSclTgiIiISMpRgSMiIiIpRwWOiIiIpBwVOCIiIpJyVOCIiIhIyskKOoB42W677bxZs2ZBhyGSEmbOnPmbuzcJOo54Ux4RiY545JC0KXCaNWvGjBkzgg5DJCWY2Y9BxxAE5RGR6IhHDtEpKhEREUk5KnBEREQk5ajAERERkZSjAkdERERSjgocERERSTkqcERERCTlqMARERGRlKMCR0RERFKOChwRERFJOSpwREREJOWowBEREZGUowJHREREUo4KHBEREUk5KnBEREQk5ajAERERkZSjAkdERERSjgocERERSTkqcERERCTlqMARERGRlKMCR0RERFKOChwRERFJOSpwREREJN7qxfoAWbE+gIiIiKS3UChETk4Obdu2Ze3atQB/j/UxVeCIiIhIzIRCIdq3b09+fj5mhplBHM4g6RSViIiIxMSqVavo378/ubm5FBUVUVhYyN577w1QFOtjqwVHREREaqz4NNQxxxyDuzN48GBeffVVcnNzi1ttqF27NkOGDKFNmzbfxDoeFTgiIiJSI8WnofLy8gBwdxo0aEDPnj3p1asXeXl5JX1wWrduDfBHrGNSgSMiIiLV4u58/PHHXH311eTm5pbc361bN4YNG0b9+vVL7osUNnGjAkdEREQqpfg0VIsWLZg/fz6DBw9m/vz5bLXVVmRmZuLu1K5dmxtuuGGT4iYIKnBERERki6ZNm0a7du0oKCjA3QE44ogjGDJkCGeeeSbz5s0rexoqUCpwREREpEIrVqzgpZdeYsCAAeTn5wNgZlx55ZU8/vjjJdu1bt06IQqbYipwREREBNh0Qr66devy0EMP8dprr5Gfn8/+++/P6tWrKSwspHbt2px99tlBh7tZKnBERESEUChEhw4dSkZCFRUVUb9+fXr16kXv3r056KCDNimAEqm1pjxJWeCY2YNAVyAf+B64wN1XBxqUiCQN5RCRTbk7gwYNYsOGDSX3HXvssbzyyitsu+22Jfcl2mmozUnWmYzHAwe4+0HAN8DNAccjIslFOUTSVigUYsCAAYRCITZu3MjIkSM5+OCDeemll0qWUqhTpw79+vXbpLhJNknZguPu/1fq5nTgtKBiEZHkoxwi6ar0aajMzEyaNGnC0qVL2W+//Rg2bBi77747U6dOTYpTUFuSlAVOGRcCr5T3gJn1BnoD7LbbbvGMSUSSR4U5BJRHJLV88MEH5Obm4u4UFRWRmZnJmDFj6Nq1KxkZ4ZM6Rx11VMBRRkfCFjhmNgHYsZyHbnX3tyLb3ApsBEaUtw93HwwMBmjZsqXHKFQRSUDRyCGgPCKpYfny5QwcOJAnnniiZA6b7OxsRo0aRZs2bQKOLjYStsBx946be9zMegInAh28+K8lIhKhHCLprHi00z777MOkSZMYMmQIeXl5nHrqqXTp0oVly5alxGmozUnYAmdzzKwzcBNwjLuvDzoeEUkuyiGSykKhEO3atSM/Px93Jysrix49enDjjTeyzz77BB1e3CTrKKongQbAeDObY2aDgg5IJBEVrJzFhm+fpmDlrBrva/Xq1Tz99NMceuihUYgscMohkpKmT5/ORRddRF5eHu6OmdGnTx+ef/75tCpuIElbcNx9r6BjEElEBStnsXHFdLIatwJg7bRzoSgfMrJp0GYEtRq1qNL+ioqKmDx5MkOGDGH06NHk5uZy8MEHxyL0uFIOkVRQetI9gFtuuYWcnBy23nprsrKyKCoqonbt2px00knBBhqQpCxwROSvClbO2qSgyd711PDPFEFRARtXTK90gbNkyRKGDh3KCy+8wIIFC9hmm2244IILuOiii2jRokXJaAsRCUbp4d5mRmFhITvttBOPPPIIvXr14rPPPkuaGYdjRQWOSIrYuGL6JgUNOGRkh3/OqFXSqlORgoIC3n33XYYMGcL7779PUVERbdu2pV+/fnTr1o2tttoqLr+HiGxZ6eHeAO3bt+ftt9+mXr16QHLNOBwrKnBEUkRW41abFDS1dz2V2rueWnLKqqLWm6+++ooXXniBYcOGsXz5cnbaaSf+/e9/c+GFF7LXXjqTI5JICgoKGDx4MI8//nhJH5vatWtz7733lhQ3EqYCRyRF1GrUggZtRvyloCmvsFm3bh2vvfYaQ4YM4eOPPyYrK4sTTzyRiy66iM6dO5OVpdQgkgiK+9kcc8wxrFq1ihtuuIGvvvqKdu3a0b17d37++ee0Pg21OcpiIkmidAfiilpjajVqUeFj7s4nn3xC//79GTduHPn5+eyzzz488MADdO/enR122CGW4YtIFZW3uvff//533nrrLbp27YqZBRxhYlOBI5IEynYgrsqIqPfff58nn3yS+fPn88MPPwBgZmRnZ/PCCy+k7CymIsnu7bff3qSfzYknnsgbb7xBrVq1Ao4sOWgohEiCKj2HTdkOxBtXTN/scwsLCxk3bhzt27enS5cujB07loULF3L00UeTkZGBu1NYWMjkyZPj88uISKX98ccf3H333QwcOLCkn02dOnW45ZZbVNxUgVpwRBJQ2RabrQ64o1IjolasWMGgQYN49tlnWbRoEVtttRVmhruTkZHBfvvtx//+9z/y8/PJzs4umT9DRIJRei6bI444guHDh3PrrbeyZMkSTjvtNM444wy+++479bOpBhU4IgmobIuNF6wqtwNxsddee42HH36YOXPmkJeXx7HHHsvDDz/M9ttvz/HHH19S0HTv3p3u3bun/fwYIomguI9Nfn4+mZmZ7L777nz99dccdthhjBo1iiOPPDLoEJOaChyRBFR2yHdxUVO6sHF3pk2bxq233lpyqikzM5Phw4dz3nnnlWw3ceLEvxQ0KmxEgpeTk0NeXh5FRUUUFhaybNkyRowYwVlnnaXJNKNABY5IAqpoyDeE+9e8+eabPPTQQ3zyySfUrVu35DQUwKJFizbZlyb8Ekk8K1asYPbs2RQVFQGQlZXFmDFjdNo4ilTgiASgOkO+161bx+23386wYcNYtWoVe+65J08++ST77rsvJ554ovrViCQBd2fw4MH07duXNWvWcPLJJ7PvvvvStWtXfRGJMhU4InFW1SHfS5cu5YknnuCpp55i7dq1AGRnZzN06NCSc/TlnYYSkcRQ3JF4n3324ZlnnmHChAm0a9eOxx9/nAMOOCDo8FKWChyROCjdYlPekO/yCpx58+bxyCOPMHLkSAoLC9lvv/348ssvS87XT5kypaTA0WkokcQUCoVo3749eXl5uDt169blmWee4ZJLLtFEfTGmAkckxqoy5NvdGThwIE8//TTfffcd9erV49JLL+Xaa6/ll19+KRlxoVNRIslhzJgx5Obmlty+4ooruPTSSwOMKH2owBGJscoM+c7Ly+O///0v9957L99//z0Q7nQ4evRojjvuOAD22GMPnYoSSRLuztChQ3nqqacAShbF7NatW8CRpQ8VOCIxtrkh3ytXruTZAQN44okn+Pnnn9lxxx1LRkS5O7NmzSopcECnokSSwZIlS+jduzdjx47l6KOP5sorr9RkfQFQgSMSY+UN+X799ddLJubLzc2lU6dODBs2jHr16tGxY0edhhJJMqFQiA8//JDc3Fwef/xx8vPzeeyxx7jyyis1p01AVOCIxEFxi82iRYu4oue/eOedd4DKT8wnIomruCNxcV+bgw46iNGjR7PXXnsFHFl6U4EjUkOVmdPm119/pX///jz99NMUFhZqYj6RFFE8MKC4uDEzzjjjDBU3CUAFjkgVlS5ogM3OaTN+/HgGDBjAJ598Qm5uLj179uTEE0/k3HPP1WkokST3yy+/cOmllzJmzJiS01C1a9emffv2AUcmUMUCx8xaAZ2BVsDOQF3gN+BrYDIwxt1XRTtIkURRdsh39q6nljunzYYNG7jpppt48skngfCpqJdffpmzzz4bSN/TUMohkgrcnVdeeYUrr7ySdevW8eCDD9KqVSumTJmSdv/TiaxSBY6Z9QBuAPYH1gJzgW+BDUAj4AjgfOApM3sV6OfuP8QkYpEAlR3yDb7JCCnf+lCeffZZ7r77bpYuXbrJcxcuXFjyc7qdhlIOkVQQCoV49913mTZtGjk5ORx++OEMHTqUfffdF0CrfyeYLRY4ZjYPaAK8BHQH5nhx54FNt9sGOBE4F/jSzHq6+ytRjlckUGWHfNfe9VRq73oq+b9OY/Tk37n7n+fy/fff06ZNG2699VZuuOGGtD8VpRwiqSAUCtG2bVvy8/MBuOyyy3j88cfJylJPj0RVmb/MEOBZd8/d3Ebu/jswAhhhZgcDO0YhPpGEUnbId9a2h/Dwww/z8MMPs2zZMg466CDeffddunTpgplxyCGHpOWpqDKUQySp/fbbb1xyySUlxU1GRga77rqripsEt8W/jrs/VtWduvtcwk3QIkmjMqOh4M8h35MnT+aqqw7ms88+C99fqxZPP/00//znP0u2TbdTUeVRDpFk9u6773LRRRexcuVKsrKyKCoqonbt2mnbIptMNPuQCH92Ht4w/2HWTjuXgpWzKtx25syZdO7cmbZt2/Ljjz+WLJhXVFTERx99FK+QRSSGioqK6NevH127dmXnnXdm1qxZfPTRR9x7771MnDgx7b+4JIMat6+ZWQZwMmDuPrrGEYnESVVW+A6FQrz66qvMnTuXDz/8kMaNG/PQQw9x6KGH0qVLl7TvZ1MTyiGSSEKhEOPGjePDDz/ko48+okePHgwaNIg6deoAqLBJItE4gfgssAOwHzDazA4A/uHur0dh35tlZn2Ah4Am7v5brI8nqaMqK3y/8cYbnHHGGRQWFgJw4YUX8uijj7L11lsD6TvkO4qUQyQhhEIh2rVrR15eHgDXXHMNjz76aEkrrSSXaBQ4bdx9fzObHbn9NfACENPkZGa7Ap2An2J5HElNlVnhe/ny5fznP//hqaeeKiluMjIy2GuvvUqKG1A/myhQDpGEMHjw4JLiJiMjgx122EHFTRKLRh+cTSblcvcC4jND8qPATcBfhpuKbEnJcG8yN1nhu+7el1NUb3/+85//sMcee/Dkk0/SpUsXateuTWZmpjoXxoZyiATK3bn//vsZNmwYZkZGRob+11NANJLITDPrSiRJmFl9oEEU9lshMzsJWOLuc1VdS3WUt8I3wFNPPUW/fv349ddf6datG/3792efffYhFArpNFTsKIdIYP744w8uuugiXnnlFc4880x69+7NJ598ov/1FGDlzLdVtR2Ek9F/gaMIz3fREfjQ3a+t4X4nUP48GLcCtwCd3P13M1sItCzv/LmZ9QZ6A+y2226H/vjjjzUJSVLY8uXL6dmzJ++//z4A2dnZ5OTkKMFVwMxmunvLKO0rYXNIZD/KIylq4cKFnHzyycybN48BAwZw00036ZRUnEQzh1R4jJoWOCU7MjsaaAksBl4rb6bSKB3nQGAisD5y1y7AUuBwd19W0fNatmzpM2bMiEVIkqAqM69NUVERQ4YM4d///jdr1qyhqKgIdyczM5N77rmHm2++Oc5RJ4dYJKdEzyGgPJIqQqEQQ4YM4fXXw928Ro0aRefOnQOOKr3Eo8CpzFIN1wGzCU+vvrqi7dz9IyDmk4C4+2fA9qXiW8hmvn1Jeio7SqrsKt8A8+bN49JLLyUUCnH00UdzySWXcPHFF2vId5Qph0gimTZtGm3btqWgoAAzU3GTwirTB+dh/jw3/hORRBW5nu3ui2MWnUg1bW5em0mTJnHXXXcxbdo0tt12W4YOHUr37t0xM3bffXf1tYk+5RBJCLm5uVx11VUUFBQAYGZ8//33AUclsVKZAmcn4JDIpXnk+iTAADezlfyZrOa4+8iYRFoBd28Wz+NJ4ip9SqrsopjF89o88MAD9O3bt+Q01PDhwzf59qYh3zGhHCKBW7x4Md26dWPWrFlaciFNVGYtql+ADyIXoKRTYHP+TFbNgWuATCCuyUkEyj8lVXqU1M/rtuPqC0/mrbfe2uR5s2fPVvN0jCmHSNA+/vhjTj31VP744w/GjBnD9ttvr5baNFCtYeLuvg6YamargG2Ag4FawGZXCxaJlfJOSdXd+3JocCCPPfYYd955J+7O5ZdfzgsvvEBBQYH62QRIOUTi5dlnn+Wqq67ib3/7G5MmTWK//fYDtORCOqhygWNmbQivG3MysCewGngPGECpb2gi8VTeKanBgwdz5513smzZMk488USeeOIJmjVrxnnnnadvbwFSDpFYC4VCTJgwgVmzZjFmzBg6d+7Mf//7Xxo2bBh0aBJHlRlFlU14XoqTga6E14z5CXgLGAN85O6FsQtRZMtKT9y3JmNfTr/w7pLTUdnZ2dx88800a9YMUD+beFMOkXgKhUK0b9+e3NxwY+D555/Piy++SGZmZsCRSbxVpgXnV6A+MBcYBLzl7nNiGZRIdWRtewij3vuCPn3OZMWKFZgZ7k5hYSGTJ0+mTZs2QYeYrpRDJG5efvnlkuLGzNh3331V3KSpyqxF1QDIB9YA2wIHm9lBZhaPtWJEShSsnMWGb5+mYOWsvzz21Vdf0b59e7p3786ee+7Jiy++SJ06dcjMzFRfm+Aph0hcDB8+nOeff75kPak6derofz+NVSbB9OLPUQ4XEv4m5kCBmc3nzzkt5hAe4rkmFoFKeqto4r6cnBzuvvtupkyZQv369Xn22We5+OKLycjIYO+991Zfm8SgHCIx9+STT3LVVVfRrl07brzxRubMmaP//TRXmWHiQ0rfNrO92XR4Z2egZ/HmhId5ikRVeaOknnrpI66//vqSOW1efvllTjjhhJLnqK9NYlAOkVh7+OGHueGGGzjppJN45ZVXqF27Nscff3zQYUnAKtPJuI27Tyu+7e7fAt8Cr5XaZnugBeGhniJRUdHEfbkbM7njwY946rlXNtl+3rx5mxQ4khiUQyQWQqEQOTk5/PTTTwwaNIjTTz+dESNGUKtWraBDkwRRmVNUU8xsOfAO8CYw0d3zS2/g7sspM5GXSE1UNHHfF5+OocetrzDv81c444wzePvttzWnTeJTDpGoKh4plZeXh7vTqVMnRo4cSVaWunXJnyrzbmhKeHjnSYSTU56ZjYv8/J7Ol0sslD0lVfBbiFFT63LVVQ+z1VZb8c4773DiiSeWfIvTufaEphwiUfXhhx+WFDdmxtFHH63iRv6iMn1wlhEe2jnIzBoAJxBOVM8Adc1sMuFE9Za7L41lsJI+Sp+SWv0HnH3W43w6+xvatWvHyy+/zM477wyon00yUA6RaHJ3Pvvss5Lipnbt2rRv3z7osCQBVWaYeAl3X+vuo9z9bKAJ4ST1PXAbsMjMPjWzm2MQp6SZ4on7JvzUhv16LuDT2d+QlZXF3XffXVLcSPJRDpGacHeuvvpqRo0axWmnnca9997LpEmT9CVHylWlAqc0dy9w9w/c/TJ3bwr8E5gEnB+16CRtFRUV8chz4zn7mhHk5TsQTm5TpkwJODKJFuUQqYqioiIuvfRSnnzySfr06cOrr77KLbfcouJGKlTtAqcsd5/u7n3dfb9o7VPS07JlyzjuuOPo27cvxxxzjCbsSxPKIVKRwsJCLrroIgYPHswtt9zCgw8+iJkFHZYkOPXKksCUHgZeq1ELAD744AO6d+/OunXrGDx4MBdffDHTp09XR2KRNDVlyhSuuuoq5s6dy1133cUdd9yh4kYqRQWOBKLsMPDaLYfSu8/jjBw5kj322IOcnBz22y/8RV4diUXS05QpU2jXrh2FhYVkZWXRqVMnFTdSaVE5RRWZmbTC2yKw6VpSpYeBL1iynlbHnMzIkSMxM5YuXcrvv/8edLgSR8ohUlZeXh6XXHIJhYXhhebdnZycnGCDkqQSrT44r2zhtqS54habDfMfZu20c7Fa20JGNq/lrKftdYtZsHgVGRkZuDsFBQVKZOlHOURK5Obm0q1bN+bPn0+tWrXUB0+qJVoFTtk2Q7UhyibKTty3dvXPXPPynlz66HKaH3wQL7/8X2rXrq1Elr6UQwSA9evX869//Yv333+fZ599lsmTJ3PPPfcwceJEnaqWKqlRHxwz241wIqplZrsC5u4/EV4wT6RE6Yn75i4opPf1j7Hgh8Xceeed3HbbbWRlZTFx4kR1Jk4zyiFS2rp16+jatSuTJ0/mhRdeoGfPngDKB1ItNe1kPCxyvRvwEuGk1B59+5JSivvc1N3/du55cAgPDM6hUaPGTJo0iWOOOaZkO3UmTkvKIUIoFOKDDz5gzJgxfP755wwfPpxzzz036LAkydWowHH3dgBmNrv45+KHahSVpIzivje/rdrAFY/9xoSZfwCwdu1asrOzA45OgqYcIsULZ+bm5gJw7733qriRqIhWHxwlIynXxhXT+WjO7xx9zWKmzFtPu0PqAqgjsZSlHJKmxo4dW1LcZGRkkJERtflnJc1F6530zRZuSxooPQwcwkXMXU/PoNsdP7N1vQw+eHA38gs1IkLKpRyShn799VdGjRoFhIub2rVrKy9I1ERloj93P2tztyX1lZ24r3D/QZzdKzzy4cLzT+b+Pkezza7H8ODheepILH+hHJJ+li1bRocOHVi8eDEDBw5k/fr1ygsSVZrJWKKi9DDwn5b9wUm9/8XSX9fz4osvloyEAGjdWiMiRNLdb7/9Rtu2bVm0aBFjx46lXbt2W36SSBWpwJGoKB4GPuvrtZx77zJWrCkkMyubffbZJ+jQRCSB5OXlccopp7Bw4ULGjx/PUUcdFXRIkqK22AfHzK4zs7Zm1jAO8UiSKe53AzB5zSWceMsyVv9RRGFReAVgdSQW5RCB8Gip/v37c8oppzB16lSGDRum4kZiqjItOA8TGeFgZj8Bs4E5kevZ7r44ZtFthpldBVwBFALvuftNQcSRzkr3uxn83jpuff43/v73v7Nw4UIyiwrUkViKKYekuVAoRIcOHcjNzcXd6dWrF2eeeWbQYUmKq0yBsxNwSOTSPHJ9EuGJuNzMVvJnsprj7iNjEmkpZtYuEsPB7p5nZtvH+pjyVxtXTKewII87XlzBoHfWcGKHgxj11jTmzZunjsRSmnJImsvJySEvLw93x8xo1qxZ0CFJGthigePuvwAfRC4AmFl9womqOX8mrWuATCDmyQm4DLjP3fMiMS6PwzGljNk/ZHHv/ct5/5M/uPRf2/Lo889Tp149zUgsm1AOke23356ioiIAateurU7FEhfVmgfH3de5+1TgQ2AB4ebnWkB+FGPbnL8DR5nZJ2Y22cwOK28jM+ttZjPMbMavv/4ap9DSw3vvvceRnS/l/U/+oH2Lrejd51HqNCn3zyDyF8mSQ0B5pKZ++uknbrvtNnbaaSduueUWJk2apC9AEhdVHkVlZm2AkyOXPYHVwHvAAEp9Q6spM5sA7FjOQ7cSjrsR0Ao4DHjVzPZw901mQ3X3wcBggJYtW2qm1GooXkcqq3ErajVqAcD8+fM5//zzKX65J8/N44OPl3LI0UFGKskimXIIKI/UxNq1a+natSvr168nFAqx3377BR2SpJEtFjhmlg10JJyMugI7AD8BbwFjgI/cvTDagbl7x83EdBnwRiQZfWpmRcB2gL5eRVHZyfsatBnB1Lm/061bt5JZRzdu3KjOxLJZyiHpqbCwkHPOOYcvvviC9957T8WNxF1lWnB+BeoDc4FBwFvuPieWQVXCGKAd8KGZ/R3IBn4LNKIUVHryPooKGPb8Y1x+23/Ze++9ee+99/j555/VmVgqQzkkjYRCIXJycpg7dy7vvvsuTz31FMcdd1zQYUkaqkyB0wDIA9YA2wIHR77tfOnuG2MZ3Ga8ALxgZp8TPmffo7ymZame4tNSVmtbyMjGC/N58NXfuX/kS7Rv357Ro0fTsGFDmjVrpsJGKkM5JE2UHQ5++umnc/nllwcdlqSpyhQ4vfhzlMOFhL+JOVBgZvP5c06LOYSHeK6JRaCluXs+cF6sj5OOyp6W+qnWmdzQfygTQ7/Ro0cPBg8eTHZ2dtBhSnJRDkkTpYeDAxx00EEBRyTprDLDxIeUvm1me7Pp8M7OQM/izQkP85QkVfq01Ko1G/jXdf9h0fKNZGVl0bt3bxU3UmXKIemjWbNmJcPB69SpQ4cOHQKOSNJZlUdRufu3wLfAa8X3RSbJagEcHL3QJAjFa0r99PN6zrj7ZxYvD59BcHcmT55MmzZtAo5Qkp1ySGr67bffuO2229h22225+OKLOeWUU3QKWwIVlcU2I5NkbTKRlySnWo1a8FXtmznl5ivZkFeLWtlGYWGhRkpJTCmHJLe8vDy6devGkiVLyMnJoVWrVkGHJFKpYeJvA3e6++zK7NDM6gCXA+vdfVAN45M4GzNmDOeccyk77LADOVPGsnr1ao2UkhpRDklt7k7v3r2ZMmUKo0aNUnEjCaMyLTgLgelmNgcYAUwF5pUe/WBmOwOHE57johuwFLgg2sFKbD322GNcd911HHbYYbz99tvssMMOACpspKYWohySkkKhEHfddRf/93//R79+/bSApiSUynQyvtrMHgOuBe4CtiG8QN4awkM/GxKeQ8KATyPbvRyLibskNgoLCzn77LN57bXXOOaYYxg7dixbbbVV0GFJilAOSU2hUIi2bduSn59PZmYmHTtWOK+iSCAq1QfH3b8HrjKzPkBr4AhgZ6AOsAL4ivBspD/GKlCJruK5bvLqHEzXc/oydepUzIxPPvmEuXPnqtVGoko5JPW8/PLL5OeHlw7TIARJRFXqZByZO2Jy5CJJqnium9/XbODUu35hznd5mBnuTkFBATk5OSpwJCaUQ1LDokWLeOWVVzAzzIzatWtrEIIknGqtJi7JbeOK6axdl8vp/X7m8x9yefKuU6hTpw6ZmZkaLSUim7Vu3Tq6du1KQUEBw4cP595772XixIn6UiQJp8bDxM2sBbBb5OZP7j6rpvuU6Cu9KnhenYM5+95fmPNdPi/2bcpZV97KIcfeoNFSEgjlkORRvIDm559/znvvvac1piShVbvAMbNDCY+IACg+b97MzADOcfeZNYxNoqT08gsbCrLo/nhDPpmfx9CHL+Ss7ldQq1ELWrfWaCmJL+WQ5HPTTTfxzjvvaAFNSQo1acEZDFzp7hNK32lmxwLPEZ6VVBJA8fIL+QWF9BzwMx/O+oahQ4dxfvfuQYcm6U05JEmEQiEeeOABxowZw1VXXaUFNCUp1KTAqV82MQG4+3gze7IG+5Uoy2rciryNmVz84M9MmLmBPZvWYe+99w46LBHlkCRQejh4RkYGp59+etAhiVRKTToZLzezC82sZGE8M8s0s4uB32oemkRLxjYHc/K9xthPNrDr9lksXBYeKSUSMOWQJPD666+XDAcHmDp1aoDRiFReTQqcHsBZwEoz+8rMvgJWAmdGHpMEUFRUxEUXXcSns78hKyuLpStcI6UkUSiHJLiNGzfywQfh5cEyMjI0HFySSrVPUbn7AqCTmW3HpiMg9M0rQbg7l19+OcOGDaNfv34ce+yxGiklCUM5JPENGDCAL7/8knvuuYfMzEzlDkkqNR4mHklGmyQkM9vZ3ZfWdN9Sfe7Oddddx7PPPkvfvn25/fbbMTMlJ0k4yiGJ6X//+x/9+vXj3HPP5bbbbgs6HJEqq3GBU4Hp/PmNTOIoFArx4YcfMn/+fF5++WWuvfZa+vfvT2TorUiyUA4J0Pr16zn//PPZaaedePJJ9feW5FSTeXD+tZmH61R3v1J9oVCIGy/tTOHGfKZ/mcspp5zCI488ouJGEpJySOL697//zddff82ECRNo2LBh0OGIVEtNWnDeJLyeTHmfng1qsF+pgtIzFM//5BWOPiCDASNzOaNtPY49clcVN5LIlEMSTCgUYvDgwQwdOpRrr72WDh06BB2SSLXVpMD5DrjQ3ReWfcDMFtVgv1JJpWcoJiObH77ehgEjV3Pq0fUYeEVj1mxbN+gQRTZHOSSBhEIh2rdvT25uLmbGSSedFHRIIjVSk2Hiw4HtK3js+RrsVyqpeIZiKGLo2N+4d9AnnNC6Pk9c3YRa2XXY45Azgg5RZHOUQxLIhx9+SF5eHgBmRigUCjgikZqpyTDxezfzWL/q7lcqL6txK8jI5r8TVtDnmd/o0ulIXhn2EBlrZ5LVuBW1Gmmme0lcyiGJpaCgAHfHzDTfjaSEKhU4ZnYtMAuY5e7rYhKRVFqtRi14/5ceXP3ELXQ45ghGvzWeOnXqwI5HBB2aSLmUQxLTokWLePTRRznwwAM588wzad++vaaUkKRX1RacAUA24Gb2PTCTcLKaSThh/R7l+KQCoVCIp556ilGjRnHkkUfx1ntjw8WNSGJTDkkwRUVF9OzZk40bN/Lmm2+y5557Bh2SSFRUtcDZHxhPeDr17wjPU9EV2IpwwppLeIXg59y9KJqByp9KL35nZtx+++3Uq1cv6LBEKkM5JEGEQiFycnL45ZdfmDRpEs8995yKG0kpVS1wBgFvuvsNxXeYWSPgDuBcYBrwH+A0Mzve3TdGLVIp8fzzz5csfmdm/O9//6Njx44BRyVSKcohCSAUCtGhQwfy8vIoKiriyCOP5KKLLgo6LJGoquooqiOBsaXvcPeV7n4t8Abhb2UtCX9LuzEaAcqmpkyZwogRIzAzLX4nyUg5JAHk5OSUFDcARx11lObMkpRT1QJnOXBIBY+9AnSPzGnxKOFvYxIFBStnseHbp5k6bhhdunRh991355133uHee+9l4sSJ6gwoyUQ5JAG0bduWjIxw+s/OzqZr164BRyQSfVU9RfUScJuZTXL32WUeawpsF/l5JnB3TYOriJk1J9zUXQfYCFzu7p/G6nhBKp7Mb+63azn59p/ZYfudmTBhAk2bNuWEE04IOjyRqlIOSQBFRUUUFRVx6KGH8sQTT+hLkqSkqhY49wAHA9PN7BXCTc3LgH2AO4EZke1qAwXRCrIcDwD93P19M+sSud02hscLzMYV0/lywTpOu/Nntt4qg/de6E3Tpk2DDkukupRDArZ27Vq6d+/O3/72Nz788EMaNNCqGJKaqlTguHsBcJKZ9SZ8fvw8wAmvJfMFcElk09bA91GM8y+hAFtHft4GWBrDYwVq3Kdr6XXHz2TXMsb8Zzf2OEitNpK8lEOCd/311/PDDz/w0UcfqbiRlFatmYzdfTAw2Mz2AHYCfgW+dXePbDIWyIlKhOW7FhhnZg8R7kfUpryNIkm0N8Buu+0Ww3Cir2DlLOZOGc45Fz1OwcYizj+uIRn73KHZiSUlJEsOgeTOI6WFQiGeeeYZhg8fTt++fTnyyCODDkkkpqq9VIOZZbv7AmBB2cfcfXqNogrvfwKwYzkP3Qp0AK5z99FmdgYwBPjLOOniJArQsmVLL/t4oipYOYvfp57Drff+RMHGIvbcuRYvjVvL3kcu5ZCjg45OJDqSIYdEYknKPFJa2YU0jz/++KBDEom5Ko2iMrMsM/uPmS0HNpjZr2Y2yswq/PZTXe7e0d0PKOfyFtCD8JBSgNeAw6N9/CBtXDGd+0f8wv/N2MC9FzaiS6t6ZGdnazi4JD3lkGB8+OGH5ObmAuG5sz7++OOAIxKJvaoOE78VuJ5wYvg3MJLwfBVTzKxvlGPbnKXAMZGf2wPfxvHYMffOtPU89OpqzuvYgAtPaMwBbXpoOLikCuWQAPzxxx8AWkhT0kpVT1F1B25x90dL32lmFwJPm9n8yLejWOsFPGZmWUAukfPjqeDzzz/nwivu4vBDD+CJgRdTb6ejuPAk9buRlKEcEmcLFizg8ccfp0WLFnTr1k0LaUraqGqBsxPhhfE24e4vmNnfgRuAmCcnd58KHBrr48TbqlWrOPnkk2nQoAFvvPUBDTUcXFKPckgcFRYW0r17dzIyMnjzzTeTupO0SFVV9RTVt0DzCh77gPD8FlINU6dO5YgjjuDHH39k9OjRmutGUpVySBw9+OCDfPzxxzz11FMqbiTtVLUFZwhwp5lNdvc5ZR7bBVgXlajSRMHKWWxcMZ2vfq5Huw4Xs3HjRmrVqqU1YSSVKYfEQSgUYuTIkQwaNIjTTz+dc8/VqheSfqpa4DxBuGPedDN7FXif8CykexGehfSD6IaXuoqXYKAon08/+oONG8OLJhcVFZGTk6Nz5JKqlENirPSQcIALL7xQX5okLVXpFFVkEq7TCY+EaAuMACYCzwKfA32iHF/K2rhiOhTl8/kPuVz3xHKabpdFRkaGhoNLSlMOib3ilcIBMjIymD277JJfIulhiy04ZtYRmOnuqwDcvQh4GHjYzPYhvDjej+6+OKaRppisxq1YsTaD8/svp2GDTF4c/BAzvlpP27Zt1XojKUU5JL7q16+Pu2tIuKS9ypyi+j/Azewnwiv8ziQ8CmKmu38NfB3D+FKWbX0Qlzy7Lb+sXsiH771A647nc2zQQYnEhnJInKxevZoHH3yQ3XbbjQsvvJBOnTrpC5OkrcoUOPsRHk55KNAC6As0IJywFhNJVJHLLHf/JUaxpoTijsU3PzqND6f8jxdffJHWHc8POiyRWFIOiZO+ffuydOlSpk2bxuGHp/TkzCJbtMUCx92/Ar4ifK4cgMh8FS34M2n1IbwirwOZMYk0BYTXmDqbVyet5LFnfuPK3mfRs2fPoMMSiSnlkPiYN28ezz33HFdeeaWKGxGqv5r4N8A3wKji+8xsL8IJSyqwYPar/Pj1Gq5/+jfaHFCbS7rtEnRIIoFQDomuadOm0b17d+rVq8edd94ZdDgiCaGqE/1VyN2/c/dXo7W/VPT+tJX0vP9Xttsmk2eu3Y5P5xcEHZJIwlAOqZ5QKES7du34/vvvyc3N5euv1aVJBKJY4MjmFRQUMOyN//HLykLaNq/L5Y+vZd8jzgw6LBFJchMmTCA/Px8IL82Qk5MTbEAiCaJap6ik6q6//nrmzJnDXXfdRXZ2NldoOLiIRMGKFSuA8Jw3GhYu8icVODFSPFoqq3Erho+Zw5NPPkmfPn10flxEombFihUMGzaMww8/nJNOOol27drpi5NIhAqcGCi9DMP/vi7gkpuXcthhh3HfffcFHZqIpJC77rqLNWvW8MILL7D//vsHHY5IQlEfnBgoXoZh2cp8et63lK3qOJ999hn/+9//gg5NRFLEl19+yTPPPMMll1yi4kakHCpwYiCrcSvyCrPoed+v/P5HEbtsl0lBQYE6/4lIjYVCIQYMGMDFF19M/fr16devX9AhiSQknaKKgVqNWnDHWwfwv6+/YZ/d6vL14nwtoikiNRYKhejQoQN5eXkUFRVx5ZVX0qRJk6DDEklIKnCiqLhj8YvvLOH5YW9w880307VrV3JycrSIpojUWPFK4UVFRQBsv/32AUckkrhU4ERJccfi6Z+v4drbl9K54z+55557yMzMVGEjIlHRtm1bMjMzKSoqIjs7m44dOwYdkkjCUh+cKNm4YjpLlq/ngvuXsdv2tXjhvpPJzNSSOiISPf/4xz+oW7cue+yxBx9++KG+PIlshlpwomRjvUPocd8v/JHrvHHPzmy3e/ugQxKRFHPPPfewdu1aJk+eTPPmzYMORyShqQUnCtyds3rdy+xv82jfoi7/fm4VM77OCzosEUkh33zzDU888QQXXXSRihuRSlCBU0MFK2fx6B1nMnbsWMyMd6atZ/oXf2hIuIhE1Y033kjdunW59957gw5FJCnoFFUNFKycxfQRp9L3voUcd1g91uZm8MmX6zUkXESiJhQK8fzzz/P2228zYMAAdthhh6BDEkkKKnBqYMOyj7ly4FIa1svgyaubULjTmbzx6TYaEi4iUREKhWjfvj25ubmYGa1atQo6JJGkoQKnikovovnIy18yb0E+w/ruQKOGdWhwyBnc3KFF0CGKSIoonvcGwMwIhUJqHRapJBU4VVB6Ec3PFhZy7wOL6NyhFadfeD5ZjVtRq5GKGxGJnpYtW+LuADr1LVJF6mRcBcWLaBZsLOSqgT9TO8vJ+XgOc347RMWNiETdhAkTALj88suZNGmSTn2LVEHCFjhmdrqZfWFmRWbWssxjN5vZd2b2tZkdF6+Yshq3goxsBr7+O5/9kM9OjbWIpkgiS8Q8UlkLFixg4MCB9OjRg6eeekrFjUgVJWyBA3wOdAM+Kn2nme0HnAXsD3QGnjazmE8ZXNz35vvMnjz82hq2a5jFwl9czcYiiS2h8khlhUIh/vWvf5GRkUH//v2DDkckKSVsHxx3nw/hjnVlnASMcvc84Acz+w44HAjFKpbivjcF+XlcfONStt22IcNeGsHs2bM1YkokgSVSHqms4o7E+fn5ZGVl8eOPP7LzzjsHHZZI0knYAmczmgLTS91eHLkvZor73jw2ehXzFuTx3ye707lzZzp37hzLw4pI7MQ9j1TWpEmTyM/PB6CoqIicnBx9iRKphkALHDObAOxYzkO3uvtbUdh/b6A3wG677Vbt/Xz1cz1WL8jnoVdXc+rRW3Pq2ZfWNDQRiZJkySOVtWHDhuLjUrt2bZ0CF6mmQAscd+9YjactAXYtdXuXyH3l7X8wMBigZcuWXo1jhSfaOu5SMjyPWpnGFTf014gpkQSSDHmkstatW8eQIUPYf//9Oeecc2jXrp1ab0SqKRlPUb0NjDSzR4Cdgb2BT6N9kOJOxfM/+Ya8vDzcnYyMDKZ/voZ2XaN9NBGJs7jkkaq6//77WbZsGW+++aZmLRapoYQtcMzsFOAJoAnwnpnNcffj3P0LM3sV+BLYCFzh7oXRPHbpCf32zCgEXM3FIkkoyDxSFaFQiDFjxjBw4EDOOeccFTciUZCwBY67vwm8WcFj/wH+E6tjl57Q7+rHfqbRNltx2ZXX06VLFzUXiySRIPNIZYVCITp06EBubi7uzmmnnRZ0SCIpIWELnCAVT+j32OhfmLcgn1eH/YfTu98QdFgikoKK15tyD7cUf/XVV0GHJJISEnmiv8DUatSCH7a5i4deW8OZpx6n4kZEYuboo48uWW9Kp8FFokcFThkFK2fx+bgb6Xb2VdSvvzVPPTsi6JBEJIUtXLiw5NSU1psSiR6doiqlYOUsfp96Ni//9xd+WvI7//hbXb755hslHBGJifXr19O3b18OPfRQXnnlFTIy9J1TJFpU4JSyccV0vvjhDx55/XdOOXIrmu2UrVlERSTqQqEQOTk5LFiwgMWLFzNy5EgVNyJRpgKnFN/6UK56/De2qZfB3Rdsy6WPreXB69oGHZaIpJDiUVN5eXkUFRXRvn17jjrqqKDDEkk5KnD4c1K/B178jM8W5PHEXafw6ZrdeHDQmWq9EZGoysnJIT8/n6KiIgBatNDM6CKxkPYFTvGkfl8sWEv/B5dw5qnHceWdbwQdloikqLZt25KdnU1ubi6ZmZl069Yt6JBEUlLaFzgbV0ynID+PKx/7lYb1MnnopvZBhyQiKax169ZMnDiRnJwc2rZtq1ZikRhJ+wInq3ErHntjLfMW5DPslqbssFd11u0TEam81q1bq7ARibG077Y/6r0veGDUSk5oty9n93lbK4WLiIikgLQucKZMmUKPHj0oLCxiYugHZnydF3RIIiIiEgVpXeDcfffdJVOkFxQUkJOTE2xAIpKyQqEQAwYMIBQKBR2KSFpI2z448+bNIycnh8zMTACys7O1BoyIxETx3Df5+flkZ2czceJE9cERibG0LHDW//IJ3c8+nUbbbs2wl0Ywe/ZsjWYQkZgpnvumsLCQ/Px8zZAuEgdpV+AUrJzFvX1OZO6XvzHslqZ0OHx7One+OeiwRCSFFc99U9yCo9ZikdhLuwJn9sejeeiV3zj16HqceEQdNq6YrpFTIhJTmvtGJP7SqsApKCjgkttG07B+Jvf1agIZtchq3CrosEQkDWjuG5H4SqsC57777mPOvK95ddiDND1iK7Iat1LrjYiISApKmwJn9W8L6dfvLo499lhO735D0OGIiIhIDKXNPDhLf17JtvXh9yUhzUMhIjGneW9EgpU2LTgb8p2nrt2Oz34o0BBNEYkpzXsjEry0acFpul0mnVrW5dOvizREU0Riqrx5b0QkvtKmwGnYqAnjlpzCg4M+0DcpEYmp4nlvMjMzNe+NSEDS5hRVnQZNufDagUGHISJpQPPeiAQvbQocEZF40rw3IsFKm1NUIiIikj5U4IiIRImGhoskDp2iEhGJAg0NF0ksCduCY2anm9kXZlZkZi1L3X+smc00s88i1+2DjFNEElc884iGhosklkRuwfkc6AY8W+b+34Cu7r7UzA4AxgFN4x2ciCSFuOWR4qHhxS04GhouEqyELXDcfT6AmZW9f3apm18Adc2strvnxTE8EUkC8cwjGhouklgStsCppFOBWRUlJTPrDfSO3Mwzs8/jFlnVbQP8nsD7rs4+tiP8TTka+97cNtV5rLKxBSXR3w97RyOQBFGVPJJrZl9E8djReB9W9+9Z1edVZfvKbpvo/4exlIi/eyzzTlmxzyHuHtgFmEC4Cbns5aRS2+QALct57v7A98CelTzWjCB/10rENziR912dfVT2Na/Mvje3TXUe0/sh+H1E8feJZx6J6u8djfdhdWOq6vOqsn1lt030/8NYXhLxd4/n/3U8jhVoC467d6zO88xsF+BNoLu7fx/dqALzToLvO+j4NrdNdR9LZEG/3vHYR1TEOY8kzO9dSnVjqurzqrJ9Ir5OsmXx/LvF/FgWqaQSlpnlADe4+4zI7YbAZKCfu79Rhf3McPeWW95SoiWRX/NEjk2iL1p5JAZxpfX7MJ1//3T+3eMlkYeJn2Jmi4HWwHtmNi7y0JXAXsAdZjYnctm+ErscHKtYpUKJ/JoncmwSJTHII9GW7u/DdP790/l3j4uEb8ERERERqaqEbcERERERqS4VOCIiIpJyVOCIiMSJme1hZkPM7PWgYymWiDFJcjGzemY2zMyeM7NzE+U4KnBEJKWZWaaZzTazd2uwjxfMbHl5k4WaWWcz+9rMvjOzvpvbj7svAPoArczsKzObb2bVmvI4mjG5+0XViUESh5nVMbNPzWxuZP21fjXYV2XfW0si/1fdgNfdvRfwrxgcp/g9XKXjpG2BE6+KU5KDvsWmtGuA+eU9YGbbm1mDMvftVc6mQ4HO5Tw/E3gKOB7YDzjbzPYzswPN7N0yl+JRWo8By9z9H8DBZWMLKKbAKSdXzmZepzygvbsfDDQHOptZqzLPjeZ761lgK6A+sAuwKLJZYazew2WPU87+NpFSBU5F1WA0KkFJPlV5P+hbbGqy8GR+JwDPV7DJMcAYM6sd2b4X8ETZjdz9I2BlOc8/HPgu8v7JB0YRnkH5M3c/scxluZltAxwNLIzsN9/dVwcZUwWvS1Skak42s13N7EMz+zLSWnJNDfYVlVYMD1sXuVkrcik7TDoq7y0gn3Bh8jqwA7CYcPEB4boiJu/hco6zWSlV4FBONRitSlCS0lAq/36Q1DQQuAkoKu9Bd3+N8Erir0S+DV8InF6F/TflzzwC4QS8uVXJmwP1gE5mttTMnjezekHGZGaNzWwQcIiZ3VyF41TGUFIzJ28E+rj7fkAr4IqyeSSIVgwLn46dAywHxrv7J6Ufj+J7ayDh/6vlQB3gDeBUM3sGeCeG7+FNjrOlnaRUgVNBNRiVSlCSTxXfD5JizOxEYLm7z9zcdu7+AJALPAP8q9S34Fj4A2gEHOvuO0du/6WPTDxjcvcV7n6pu+/p7gOivO+UzMnu/rO7z4r8vJbwacayRWTcWzHcvdDdm0ceP9zMDijneDV9b+1Kmf8rd//D3S9w98vcfUSUjvMX5R1ncxL6TRQlUakEJWWU+36I8bdYCcY/gX+Z2ULCHw7tzezlshuZ2VHAAYTXpbqzisdYQjjhF9slcl9FFgOLS32zfh1oEXBM8ZZSOdnMmgGHALFqLSlW6dcpctrzQ8pvGarpe2tv/vy/uppw/gz6/6pc6VDglKuqlaCktlh+i5VguPvN7r6LuzcDzgImuft5pbcxs0MIT5l/EnAB0NjM7q3CYf4H7G1mu5tZduQ4b28mpmXAIjPbJ3JXB+DLIGNKFMmYk82sPjAauNbd15R9PJ6tGGbWxMJrrGFmdYFjga/KxBuN9xbAUcDfgV+AaUH/X1UkHQqcRP82I/Gl94OUthVwhrt/7+5FQHfgx7Ibmdl/gRCwj5ktNrOLANx9I+F1rcYRPk3xqrt/sYVjXgWMMLN5hPvk9E+AmOIpJf4HzawW4eJmhFewYGucWzF2Aj6MvK/+R7gPTtmpEaL93soByivaEuI9nHJrUUWaC9919wMit7OAbwh/U1pC+A9/ToL9w0uM6P0gEqxU/B80MwOGASvd/doKtjkEGAmcCPwAjAC+d/fbytm2GaVeo8h9Sf86BS2lWnDKqwaT4NuMxIjeDyLBSuH/wX8C5xPu11W8Gn2XMtskRCtGOku5FhwRERGRlGrBEREREQEVOCIiIpKCVOCIiIhIylGBIyIiIilHBY6IiIikHBU4IiIiknJU4IiIiEjKUYEjIiIxYWYXm5lXcPnDzDK38PyhpbbPKefx1mY2KjJBXr6ZrTGz/5nZPWa2UzXiHWNmq4pXAC/n8QaRuIdGbt9WKr7FVT2exJYKHIm7Uknv0KBjEZGYOgT4A2hdzqWNuxdWYh/LIttfXvpOM+sDfAw0AW4DOhJelHEc0Bt4oRrxDgMaEl5eoTynEZ6heFjk9ouR2MZW41gSY1lBByBp6RBgI/B50IGISEw1B75w9+k12Ede2eebWTvgQeAxd7+uzPZjzWwAcHo1jvUesILwsgqjy3m8O/AT4UUmcfclwBIz+7Uax5IYUwuOBKE5MN/d84IORERiI7Ig5UHAvBjs/t/Ab5Hrv3D3P9x9aDkxHWxmb0dOQ20ws48jK34XPy8f+C9wvJk1LvPc3YBjgOGuNY6SggociatSSW920LGISEztDdQHvjCzrDKXzfa92ZzIKtvHAOMjBUlln9cCmAY0AnoBpxJurZlQ5nT5MKAW4dNdpZ0HGPBSdWOX+FKBI/FWnPRU4IiktuaR60eBgjKXGTXYb2OgDuFTRZsoW0iVefjByHPau/vr7j4WOAVYANxevJG7zwC+JHw6qrTzgenu/k0NYpc4Uh8cibfmkes5AcYgIrHXHHDgSKBsS8uq4h/MbE/CrSbbE+6Q3CtSZFSJme0I/FzmvlruvtHM6hJu9ekPFJUpfiYA55bZ3TDgfjP7u7t/Y2aHA/8ALqtqXBIcteBIvDWPXM8xs23N7F0z+8bM5prZ/5nZXsUbmtntkceKzOzkQKIVkeo6BFjg7tPcfUaZy/elthsEDHP3vwM3ASMip7IrsgLIBXYrc/9vwGGRy3NlHmsEZBJuqSnbmnQlsK2Zlf48fBko4s9WnO5AHvBKJX5vSRAqcCTeDgF+cPfVhL/dDXT3v7v7wcC7wPOlth0PdAY+inuUIlJTzYHPNreBmTUBWgFDAdx9POF+LhVOIeHuGwnnhGPNLLv0/cUFFLC0zNNWEy5YnuDPImiTi7sXldrXUsL557zIMc4E3nH3VUjSUIEj8dacyOkpd1/t7hNKPTYNaFZ8w92nu/uCeAYnIjVnZjsAO7LlEVS7AT+7e0Gp+xby19aZsh4AtgPur0w87v4HMAU4GJhVTotSeafEhgF/AwZEjjWsnG0kgakPjsRNqaRXUQfja4G34haQiMTKIZHrjWbWqpzH57n7+uru3N0nmllf4D4zO4jwyKYfCHc+/jvhEVB/EG4lLnY94ZafcWY2hHB/ne2AFkCmu/ctc5gxwBrgOmA58EF145VgqMCReCpOen8pcMzsTmAPwjOQikhyax65vjtyKc0J94lZT3hU006RzsDFrTjNKGeEVFnu/oCZfQxcQ7jzcBPCfXO+JtxXZlDpmZLdfZaZHQbcCTwObAP8Cswi3A+o7P43mNmrwMXAyMipMUkiKnAknppHrueUvtPMbgO6AJ1q8q1ORBKDu98H3FeJ7X41s0+BnsBzZnYs4T44M0tvFxn15GWXdnD3jwkv11DZuObz1/ltNrd9L8Jz5pQr0hk6MxKzJBj1wZG4cff73N3cvWRRukjLTVfCxc3vwUUnIgG5FLjAzL4hPFfNuWVmCv4b4dFOE4MIbgtuJRxb2TlzJAGYZpyWoJjZ/oTXo/oeWBe5e6O7t4w8fhfh5uEmwFrCzc+tShdIIpK6zKwZ4X4yAGvd/esAw/mLyIrlTSM38909FstSSDWpwBEREZGUo1NUIiIiknJU4IiIiEjKUYEjIiIiKUcFjoiIiKQcFTgiIiKSclTgiIiISMpRgSMiIiIpRwWOiIiIpJz/BwOioDDIL5N6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot all\n",
    "fntsz = 16\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2,figsize=(8,8))\n",
    "\n",
    "axes[0,0].set_xscale('log')\n",
    "axes[0,1].set_xscale('log')\n",
    "axes[1,0].set_xscale('log')\n",
    "axes[1,1].set_xscale('log')\n",
    "\n",
    "f = open('data/fig6a/wpd.json')\n",
    "data = json.load(f)\n",
    "datatmp = data['datasetColl'][0]['data']\n",
    "datatmp2 = np.array([datatmp[i]['value'] for i in range(len(datatmp))])\n",
    "indxs = np.argsort(datatmp2[:,0])\n",
    "axes[0,0].plot(datatmp2[:,0][indxs],datatmp2[:,1][indxs],'k.')\n",
    "datatmp = data['datasetColl'][1]['data']\n",
    "datatmp2 = np.array([datatmp[i]['value'] for i in range(len(datatmp))])\n",
    "indxs = np.argsort(datatmp2[:,0])\n",
    "axes[0,0].plot(datatmp2[:,0][indxs],datatmp2[:,1][indxs],'.',color='#EDB120')\n",
    "\n",
    "#\n",
    "f = open('data/fig6b/wpd.json')\n",
    "data = json.load(f)\n",
    "datatmp = data['datasetColl'][0]['data']\n",
    "datatmp2 = np.array([datatmp[i]['value'] for i in range(len(datatmp))])\n",
    "indxs = np.argsort(datatmp2[:,0])\n",
    "axes[0,1].plot(datatmp2[:,0][indxs],datatmp2[:,1][indxs],'.',color='#EDB120')\n",
    "\n",
    "#\n",
    "f = open('data/fig6c/wpd.json')\n",
    "data = json.load(f)\n",
    "datatmp = data['datasetColl'][0]['data']\n",
    "datatmp2 = np.array([datatmp[i]['value'] for i in range(len(datatmp))])\n",
    "indxs = np.argsort(datatmp2[:,0])\n",
    "axes[1,0].plot(datatmp2[:,0][indxs],datatmp2[:,1][indxs],'k.')\n",
    "datatmp = data['datasetColl'][1]['data']\n",
    "datatmp2 = np.array([datatmp[i]['value'] for i in range(len(datatmp))])\n",
    "indxs = np.argsort(datatmp2[:,0])\n",
    "axes[1,0].plot(datatmp2[:,0][indxs],datatmp2[:,1][indxs],'.',color='#EDB120')\n",
    "\n",
    "#\n",
    "f = open('data/fig6d/wpd.json')\n",
    "data = json.load(f)\n",
    "datatmp = data['datasetColl'][0]['data']\n",
    "datatmp2 = np.array([datatmp[i]['value'] for i in range(len(datatmp))])\n",
    "indxs = np.argsort(datatmp2[:,0])\n",
    "axes[1,1].plot(datatmp2[:,0][indxs],datatmp2[:,1][indxs],'k.')\n",
    "\n",
    "\n",
    "axes[0,0].plot(taulst,1e6*sadle_tau,'k-')\n",
    "#axes[0,0].plot(taulst,1e6*full_tau,'b-')\n",
    "axes[0,0].set_xlim([4,120])\n",
    "axes[0,0].set_ylim([0,5])\n",
    "axes[0,0].set_xlabel(r'$\\tau[\\mathrm{fs}]$',size=fntsz)\n",
    "axes[0,0].set_ylabel(r'$10^6 N_\\pm /N_e$',size=fntsz)\n",
    "axes[0,0].set_xticks([5,10,50,100])\n",
    "axes[0,0].set_yticks([0,1,2,3,4,5])\n",
    "\n",
    "axes[0,1].plot(lbdlst, lbdlst*0+1.8735083532219514,'k.')\n",
    "axes[0,1].plot(lbdlst,1e6*saddle_lbd,'k-')\n",
    "#axes[0,1].plot(lbdlst,1e6*full_lbd,'b-')\n",
    "axes[0,1].set_xlim([0.1,20])\n",
    "axes[0,1].set_ylim([0,5])\n",
    "axes[0,1].set_xlabel(r'$\\lambda[\\mu \\mathrm{m}]$',size=fntsz)\n",
    "axes[0,1].set_ylabel(r'$10^6 N_\\pm /N_e$',size=fntsz)\n",
    "axes[0,1].set_xticks([0.1,0.5,1,5,10])\n",
    "axes[0,1].set_yticks([0,1,2,3,4,5,6])\n",
    "\n",
    "axes[1,0].plot(I21lst,log10(saddle_I21),'k-')\n",
    "#axes[1,0].plot(I21lst,log10(full_I21),'b-')\n",
    "axes[1,0].set_xlim([1,40])\n",
    "axes[1,0].set_ylim([-12,0])\n",
    "axes[1,0].set_xlabel(r'$I_{21}$',size=fntsz)\n",
    "axes[1,0].set_ylabel(r'$log_{10}(N_\\pm /N_e)$',size=fntsz)\n",
    "axes[1,0].set_xticks([1,5,10])\n",
    "axes[1,0].set_yticks([-12,-10,-8,-6,-4,-2])\n",
    "\n",
    "axes[1,1].plot(E0lst,log10(saddle_E0),'k-')\n",
    "#axes[1,1].plot(E0lst,log10(full_E0),'b-')\n",
    "axes[1,1].set_xlim([0.4,4])\n",
    "axes[1,1].set_ylim([-12,0])\n",
    "axes[1,1].set_xlabel(r'$E_{0}[\\mathrm{GeV}]$',size=fntsz)\n",
    "axes[1,1].set_ylabel(r'$log_{10}(N_\\pm /N_e)$',size=fntsz)\n",
    "axes[1,1].set_xticks([0.5,1,2])\n",
    "axes[1,1].set_yticks([-12,-10,-8,-6,-4,-2])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "QPPQ",
   "language": "python",
   "name": "qppq"
  },
  "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}