{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fa7ba993",
   "metadata": {},
   "source": [
    "# Gravitational wave production"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f719b13",
   "metadata": {},
   "source": [
    "## Importing the modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3bee1551",
   "metadata": {},
   "outputs": [],
   "source": [
    "#start by importing the controller and manipulate modules\n",
    "from analytical.controller import *\n",
    "from numerical.manipulate import *\n",
    "#reimport numpy (though it is pulled by numerical) for smarter syntax highlighting in vscode\n",
    "import numpy as np\n",
    "#import matplotlib too\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4579b931",
   "metadata": {},
   "source": [
    "## The Standard Model, reproducing the results of [2004.11392](https://arxiv.org/abs/2004.11392)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3369f4a6",
   "metadata": {},
   "source": [
    "The relevant part of the config file is as follows\n",
    "```ini\n",
    "[Model]\n",
    "modelpath = /Users/jacopo/NextCloud/AUTOTHERM/autotherm/analytical/models/SPSM_GW.fr\n",
    "# Symbol for the Lagrangian in the model file\n",
    "lagrangian = Ltot\n",
    "# \"Name\" of the particle whose production rate must be computed\n",
    "produced = T[1]\n",
    "# List of the particles in the thermal bath (or leave empty for SM assumption)\n",
    "# Unused for now, but still required.\n",
    "assumptions = Element[ht,Reals], Element[kappa,Reals]\n",
    "replacements =\n",
    "includeSM = yes\n",
    "noneq = kappa\n",
    "flavorexpand = \n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18a6bc76",
   "metadata": {},
   "outputs": [],
   "source": [
    "GWdict=analytical_pipeline(\"../../MyModels/grav/grav.cfg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a7ab001",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "-  $$ \\left\\vert \\mathcal{M}(1;-1,-1)\\right\\vert^2=6 \\vert h_t\\vert^2 \\kappa^{2} s + \\frac{2 \\kappa^{2} \\left(t^{2} + u^{2}\\right) \\left(5 g_{1}^{2} + 9 g_{2}^{2} + 24 g_{3}^{2}\\right)}{s}$$\n",
       "-  $$ \\left\\vert \\mathcal{M}(-1;-1,1)\\right\\vert^2=- 6 \\vert h_t\\vert^2 \\kappa^{2} t - \\frac{2 \\kappa^{2} \\left(s^{2} + u^{2}\\right) \\left(5 g_{1}^{2} + 9 g_{2}^{2} + 24 g_{3}^{2}\\right)}{t}$$\n",
       "-  $$ \\left\\vert \\mathcal{M}(-1;1,-1)\\right\\vert^2=- 6 \\vert h_t\\vert^2 \\kappa^{2} u - \\frac{2 \\kappa^{2} \\left(s^{2} + t^{2}\\right) \\left(5 g_{1}^{2} + 9 g_{2}^{2} + 24 g_{3}^{2}\\right)}{u}$$\n",
       "-  $$ \\left\\vert \\mathcal{M}(1;1,1)\\right\\vert^2=\\frac{\\kappa^{2} \\left(g_{1}^{2} + 15 g_{2}^{2} + 48 g_{3}^{2}\\right) \\left(s^{2} + t^{2} + u^{2}\\right)^{2}}{4 s t u}$$\n"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, Markdown\n",
    "out=\"\"\n",
    "s = sympy.Symbol(\"s\")\n",
    "t = sympy.Symbol(\"t\")\n",
    "u = sympy.Symbol(\"u\")\n",
    "g1 = sympy.Symbol(\"g1\")\n",
    "g2 = sympy.Symbol(\"g2\")\n",
    "g3 = sympy.Symbol(\"g3\")\n",
    "ht = sympy.Symbol(\"ht\")\n",
    "kappa=sympy.Symbol(\"kappa\")\n",
    "for key, item, in GWdict[0].items():\n",
    "    sitem=sympy.simplify(sympy.sympify(item).subs(2*t*u,(s*s-t*t-u*u)))\n",
    "    sitemsubst= sitem\n",
    "    sitemsubstcollect=0\n",
    "    for exprtemp in sympy.factor(sitemsubst.expand().as_independent(ht)).subs(t+u,-s)\\\n",
    "        .subs(2*t*t+2*t*u,2*t*t+(s*s-t*t-u*u)).subs(t*t+2*t*u,t*t+(s*s-t*t-u*u)).subs(t*t+t*u,t*t+(s*s-t*t-u*u)/2):\n",
    "        sitemsubstcollect += exprtemp.simplify()\n",
    "        # beautify the output\n",
    "    out+=f\"-  $$ \\\\left\\\\vert \\\\mathcal{{M}}({key[0]};{key[1]},{key[2]})\\\\right\\\\vert^2={sympy.latex(sitemsubstcollect).replace('ht^{2}',r'\\vert h_t\\vert^2')}$$\\n\"\n",
    "    # display(sympy.pprint(sitem))\n",
    "display(Markdown(out))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9e3fde99",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{T^{3} \\kappa^{2} \\left(11 g_{1}^{2} \\log{\\left(\\frac{24 k^{2}}{11 T^{2} g_{1}^{2}} \\right)} + 33 g_{2}^{2} \\log{\\left(\\frac{24 k^{2}}{11 T^{2} g_{2}^{2}} \\right)} + 96 g_{3}^{2} \\log{\\left(\\frac{2 k^{2}}{T^{2} g_{3}^{2}} \\right)}\\right)}{192 \\pi}$"
      ],
      "text/plain": [
       "T**3*kappa**2*(11*g1**2*log(24*k**2/(11*T**2*g1**2)) + 33*g2**2*log(24*k**2/(11*T**2*g2**2)) + 96*g3**2*log(2*k**2/(T**2*g3**2)))/(192*pi)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "GW=NumRate(*GWdict,2)\n",
    "GW.get_leadlog()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcd52177",
   "metadata": {},
   "source": [
    "make the masses explicit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "79fd6851",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{T \\kappa^{2} \\left(m_{D1}^2 \\log{\\left(\\frac{4 k^{2}}{m_{D1}^2} \\right)} + 3 m_{D2}^2 \\log{\\left(\\frac{4 k^{2}}{m_{D2}^2} \\right)} + 8 m_{D3}^2 \\log{\\left(\\frac{4 k^{2}}{m_{D3}^2} \\right)}\\right)}{32 \\pi}$"
      ],
      "text/plain": [
       "T*kappa**2*(m_{D1}^2*log(4*k**2/m_{D1}^2) + 3*m_{D2}^2*log(4*k**2/m_{D2}^2) + 8*m_{D3}^2*log(4*k**2/m_{D3}^2))/(32*pi)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "md1 = sympy.Symbol(\"m_{D1}^2\")\n",
    "md2 = sympy.Symbol(\"m_{D2}^2\")\n",
    "md3 = sympy.Symbol(\"m_{D3}^2\")\n",
    "T = sympy.S(\"T\")\n",
    "GW.get_leadlog().subs(g1*g1,6*md1/(11*T**2)).subs(g2*g2,6*md2/(11*T**2)).subs(g3*g3,md3/(2*T**2)).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01107a85",
   "metadata": {},
   "source": [
    "Illustration plot at the three temperatures in [2004.11392](https://arxiv.org/abs/2004.11392)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fb4e6e69",
   "metadata": {},
   "outputs": [],
   "source": [
    "k=numpy.geomspace(0.01,15,200)\n",
    "tempslog=[15,9,3]\n",
    "couplings=[(numpy.sqrt(0.162787),numpy.sqrt(0.275424),numpy.sqrt(0.27615),numpy.sqrt(0.184357)),\n",
    "           (numpy.sqrt(0.145664),numpy.sqrt(0.325047),numpy.sqrt(0.41907),numpy.sqrt(0.295078)),\n",
    "           (numpy.sqrt(0.131767),numpy.sqrt(0.396405),numpy.sqrt(0.872336),numpy.sqrt(0.626187))]\n",
    "GWtestrate = [GW.rate(k,numpy.sqrt(32*numpy.pi),couptuple,0) for couptuple in couplings]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "450cd79e",
   "metadata": {},
   "outputs": [],
   "source": [
    "GWtestratefallback = [GW.rate(k,numpy.sqrt(32*numpy.pi),couptuple,-1) for couptuple in couplings]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ab9eb251",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwgAAAFcCAYAAABshZnfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACRt0lEQVR4nOzdd1yV1R/A8c9lgwIqIqAioLhw4MCBe8/cOy1NyjQ1R6kttczUzJUjR7my3KZmlovce+GejFyAisje9/z+uHp/4gQZF+T7fr1ucZ/n3PN87zjX53ufMzRKKYUQQgghhBBCAEaGDkAIIYQQQgiRc0iCIIQQQgghhNCTBEEIIYQQQgihJwmCEEIIIYQQQk8SBCGEEEIIIYSeJAhCCCGEEEIIPUkQhBBCCCGEEHqSIAghhBBCCCH0TAwdgBB5jVKK5ORkUlJSDB2KEELkSaamphgbGxs6DCFyLEkQhMhGiYmJBAcHExsba+hQhBAiz9JoNBQvXpz8+fMbOhQhciSNUkoZOggh8gKtVsu1a9cwNjbG3t4eMzMzNBqNocMSQog8RSnFvXv3iI2NpXTp0nIlQYjnkCsIQmSTxMREtFotzs7OWFlZGTocIYTIs+zt7QkKCiIpKUkSBCGeQwYpC5HNjIyk2QkhhCHJ1VshXk7OVIQQQgghhBB6kiAIIYQQQggh9CRBEEJki2XLllGgQAFDhyGymKurK7NmzTJ0GFni66+/pkqVKoYOQ6TBnj170Gg0PHz40NChCJErSYIghHht/fr1o2PHjmkq26NHD65evZrmutN6ovmqcjdv3qR///4ULVoUMzMzXFxcGDZsGGFhYWmORWQujUbDpk2bsuVYclKffRo1asTw4cMNHYYQIhNIgiCEyHJJSUlYWlpSpEiRbD1uQEAAXl5eXLt2jVWrVnH9+nUWLFiAr68v3t7ePHjwIFvjEWmXmJho6BCEECLPkgRBCEOLiXnxLT4+7WXj4tJWNp3Wr19PpUqVsLS0xM7OjmbNmhETE8PXX3/N8uXL2bx5MxqNBo1Gw549ewgKCkKj0bBmzRoaNmyIhYUFv//++3O7GG3ZsoUaNWpgYWFB4cKF6dSpE6D7JfK///5jxIgR+rpfx+DBgzEzM2PHjh00bNiQEiVK0Lp1a3bt2sXt27f58ssvX6vejHr8Vjy5Ck1iom5bQsLzy2q1/9+WlPTyj8fTZdPrRe/5834h7tixI/369Uu1LSoqil69epEvXz6KFSvGvHnz9PtcXV0B6NSpExqNRn//8S/9v/zyC25ublhYWACwbds26tWrR4ECBbCzs+Ott97C398/1fFu3bpFr169KFSoEPny5cPLy4ujR4+ybNkyvvnmG86cOaP/HC1btgyAhw8f8v7772Nvb4+NjQ1NmjThzJkzqeqdMmUKDg4OWFtb4+PjQ/zTL3g2UUoRExOT7bf0LJPUr18/9u7dy48//pjqtX66zW/atClVe378vq9YsQJXV1dsbW3p2bMnUVFR+jJarZbJkyfj5uaGpaUlnp6erF+/PlW9f//9N2XKlMHS0pLGjRsTFBT0Wq+1EOIRJYTIFnFxcerixYsqLi4u9Q7deeLzb23apC5rZfXisg0bpi5buPDzy6XDnTt3lImJiZoxY4YKDAxUZ8+eVfPmzVNRUVEqKipKde/eXbVq1UoFBwer4OBglZCQoAIDAxWgXF1d1YYNG1RAQIC6c+eOWrp0qbK1tdXX/ddffyljY2M1btw4dfHiReXn56cmTZqklFIqLCxMFS9eXE2YMEFf94u4uLiomTNnPrM9LCxMaTQafZ1P++CDD1TBggWVVqtN12uSGR6/FXfv/n/bxIm6be+/n7rs47c8MPD/22bO1G17++3UZR+/5efP/3/bokXpi+1l73nDhg3VsGHDUpXv0KGD6tu3r/6+i4uLsra2VpMnT1ZXrlxRs2fPVsbGxmrHjh1KKaXu3r2rALV06VIVHBys7j56EcaPH6/y5cunWrVqpU6dOqXOnDmjlFJq/fr1asOGDeratWvq9OnTql27dqpSpUoqJSVFKaVUVFSUKlmypKpfv77av3+/unbtmlqzZo06dOiQio2NVZ988omqUKGC/nMUGxurlFKqWbNmql27dur48ePq6tWr6pNPPlF2dnYqLCxMKaXUmjVrlLm5ufrll1/U5cuX1Zdffqmsra2Vp6dn+l7QTBAdHa2AbL9FR0enOcaHDx8qb29v9cEHH+hf619++SVVm1dKqY0bN6onTz3Gjx+v8ufPrzp37qzOnTun9u3bpxwdHdUXX3yhLzNx4kRVrlw5tW3bNuXv76+WLl2qzM3N1Z49e5RSSt24cUOZm5urkSNHqsuXL6vffvtNOTg4KECFh4c/N94Xfh8LIZRSSslCaUKIFwoODiY5OZnOnTvj4uICQKVKlfT7LS0tSUhIwNHR8ZnHDh8+nM6dO7+w7u+++46ePXvyzTff6Ld5enoCUKhQIYyNjbG2tn5u3Wlx7do1lFKUL1/+ufvLly9PeHg49+7dy/auTznZq97ztKhbty6fffYZAGXKlOHgwYPMnDmT5s2bY29vD0CBAgWeeW8TExP59ddf9WUAunTpkqrMkiVLsLe35+LFi1SsWJGVK1dy7949jh8/TqFChQBwd3fXl8+fPz8mJiapjnXgwAGOHTvG3bt3MTc3B2DatGls2rSJ9evXM2DAAGbNmoWPjw8+Pj4ATJw4kV27dhnsKkJOZ2tri5mZGVZWVvrXOq0LkGm1WpYtW4a1tTUA77zzDr6+vnz33XckJCQwadIkdu3ahbe3NwAlS5bkwIEDLFy4kIYNGzJ//nxKlSrF9OnTAShbtiznzp3j+++/z4JnKkTeIAmCEIYWHf3ifU//A3v37ovLPr0AWyZcYvf09KRp06ZUqlSJli1b0qJFC7p27UrBggVf+VgvL6+X7vfz8+ODDz7IcIyvotLRTSK7PH7Ln1xQe9QoGD4cTJ76Vn78llta/n/b4MHwwQfPfjwev+VPln2q988rZeQ9f+zxidyT99My4NzFxSVVcgC6RG/cuHEcPXqU+/fvo33Uf+rGjRtUrFgRPz8/qlatqk8O0uLMmTNER0djZ2eXantcXJy++9KlS5cYOHDgM89j9+7daT5OZrGysiL6Zd8TWXjc7ODq6qpPDgCcnJy4++iDf/36dWJjY2nevHmqxyQmJlK1alVA917VqlUr1f6nP4NCiPSRBEEIQ8uXz/BlX8DY2JidO3dy6NAhduzYwZw5c/jyyy85evQobm5urzj8y49v+eRZbBZwd3dHo9Fw6dIl/diGJ126dImCBQs+c0KaHZ730piZ6W5pKWtqqrultWx6vOw9NzIyeibhSnqdQQ4v8LzPTLt27XBxceHnn3+maNGiaLVaKlasqB/E/Dqfo+joaJycnNizZ88z+3LiVLwajeaV7SknSuvnxfSpD6lGo9Engo8To61bt1KsWLFU5R5f/RFCZD4ZpCyEeCmNRkPdunX55ptvOH36NGZmZmzcuBEAMzMzUlJSXqveypUr4+vr+8L9GakbwM7OjubNm/PTTz8R99QA7pCQEH7//Xd69Ojx2gOg32Qves/t7e0JDg7Wl0tJSeH8+fPPPP7IkSPP3H+yq5epqWma3tuwsDCuXLnCV199RdOmTfXdwp5UuXJl/Pz8Xjgj1fM+R9WqVSMkJAQTExPc3d1T3QoXLgzouqAdPXr0pc9LpPb0a21vb09UVBQxT0yO4Ofnl646PTw8MDc358aNG8+8V87OzoDuvTp27Fiqx8l7JUTGSIIghHiho0ePMmnSJE6cOMGNGzf4448/uHfvnv5kz9XVlbNnz3LlyhXu37+frl+Tx48fz6pVqxg/fjyXLl16ps+wq6sr+/bt4/bt29y/f/+ldd2+fRs/P79Ut/DwcObOnUtCQgItW7Zk37593Lx5k23bttG8eXOKFSvGd99993ovzBvsZe95kyZN2Lp1K1u3buXy5csMGjTouQtRHTx4kKlTp3L16lXmzZvHunXrGDZsmH6/q6srvr6+hISEPHPC/6SCBQtiZ2fHokWLuH79Ov/++y8jR45MVaZXr144OjrSsWNHDh48SEBAABs2bODw4cP6YwUGBuLn58f9+/dJSEigWbNmeHt707FjR3bs2EFQUBCHDh3iyy+/5MSJEwAMGzaMJUuWsHTpUq5evcr48eO5cOFCJrzCby5XV1eOHj1KUFAQ9+/fp1atWlhZWfHFF1/g7+/PypUr9bNIpZW1tTWffvopI0aMYPny5fj7+3Pq1CnmzJnD8uXLARg4cCDXrl1j1KhRXLly5bWOI4R4imHHSAuRd+TGWTMuXryoWrZsqezt7ZW5ubkqU6aMmjNnjn7/3bt3VfPmzVX+/PkVoHbv3q2fxej06dOp6np6FiOllNqwYYOqUqWKMjMzU4ULF1adO3fW7zt8+LCqXLmyMjc3Vy/7qnJxcXnuDCwrVqxQSikVFBSk+vbtqxwcHJSpqalydnZWQ4cOVffv38/4C/QGetl7npiYqAYNGqQKFSqkihQpoiZPnvzcWYy++eYb1a1bN2VlZaUcHR3Vjz/+mOoYf/75p3J3d1cmJibKxcVFKaWbzeZ5MwTt3LlTlS9fXpmbm6vKlSurPXv2KEBt3LhRXyYoKEh16dJF2djYKCsrK+Xl5aWOHj2qlFIqPj5edenSRRUoUEA/e5JSSkVGRqqhQ4eqokWL6j8XvXv3Vjdu3NDX+91336nChQur/Pnzq759+6rRo0cbZBaj3OLKlSuqdu3aytLSUgEqMDBQbdy4Ubm7uytLS0v11ltvqUWLFj0zi9HTr+nMmTP1nwullNJqtWrWrFmqbNmyytTUVNnb26uWLVuqvXv36sts2bJFubu7K3Nzc1W/fn21ZMkSmcVIiAzQKJUDR/AJ8QaKj48nMDAw1RzvQgghsp98HwvxctLFSAghhBBCCKEnCYIQQgghhBBCTxIEIYQQQgghhJ4kCEIIIYQQQgg9SRCEyGYyL4AQQhiWfA8L8XKSIAiRTR6vFhobG2vgSIQQIm97vBK3sbGxgSMRImcyMXQAQuQVxsbGFChQgLt37wJgZWUlq/gKIUQ202q13Lt3DysrK0xM5DRIiOeRliFENnJ0dATQJwlCCCGyn5GRESVKlJAfaYR4AVkoTQgDSElJISkpydBhCCFEnmRmZoaRkfSyFuJFJEEQQgghhBBC6En6LIQQQgghhNCTBEEIIYQQQgihJwmCEEIIIYQQQk8SBCGEEEIIIYSeJAhCCCGEEEIIPUkQhBBCCCGEEHqSIIhca/78+VSrVg1TU1O+/vrrZ/Y3atQICwsL8ufPT/78+WndunX2BymE0Ltw4QINGjTAxsYGDw8P9uzZY+iQhBCPDBgwACcnJ2xsbKhUqRJbtmwxdEjCgCRBELmWk5MTX3/9NV26dHlhmV9++YXo6Giio6P5559/sjE6IcSTkpKS6NChA127diU8PJzZs2fTtWtXwsLCDB2aEAIYOXIkQUFBREZGsmTJEvr06SPtMw+TBEHkWh07dqR9+/YUKFDA0KEIIV7hypUrhIeH8/HHH2NsbEyzZs2oWrUqGzduNHRoQgigXLlymJubA6DRaEhMTOT27dsGjkoYiiQIecjNmzfRaDRpugUEBGTqsaOjoxk/fjytWrWiUKFCaDQali1b9sLyCQkJjBkzhqJFi2JpaUmtWrXYuXNnuo87YsQI7O3tad68OWfPnk3z4wIDAxkyZAhlypTBysoKKysrPDw8GDx4cLrqeax9+/ZYWVkRFRX1wjK9e/fGzMxMfrERem9am1VKPXP/woULaY5J2qXISd609gnw0UcfYWlpSY0aNWjSpAmVKlVKc0zSPt8sJoYOQGQfc3NzVqxYob8fFxfHgAEDaNy4Mf3799dv12g0lCxZMlOPff/+fSZMmECJEiXw9PR8Zd/jfv36sX79eoYPH07p0qVZtmwZbdq0Yffu3dSrVy9Nx5w6dSoeHh4YGxszZ84cWrduzeXLl7G2tn7p4/766y969OiBiYkJvXv3xtPTEyMjIy5fvswff/zB/PnzCQwMxMXFJa1Pn969e7NlyxY2btzIu++++8z+2NhYNm/eTKtWrbCzs0tzveLN9ia12bJly1KgQAFmzJjB0KFD8fX1Ze/evWmOW9qlyGnepPb52E8//cScOXPYs2cP58+fR6PRpCkeaZ9vICXyrOPHjytATZs2LcuPFR8fr4KDg1Mdd+nSpc8te/ToUQWoH374Qb8tLi5OlSpVSnl7ez9T/sMPP1Tjx49/ZQxly5ZVO3bseGmZ69evq3z58qny5curO3fuPLM/KSlJ/fjjj+rGjRuvPN6TYmNjlbW1tWrZsuVz969cuVIBavXq1emqV+Qtub3NnjlzRjVo0EAVKlRItWjRQvXu3Vt98803r4xF2qXIDXJ7+3zaW2+9pbZu3frKWKR9vpmki1Ee9viSX3ouIb4uc3NzHB0d01R2/fr1GBsbM2DAAP02CwsLfHx8OHz4MDdv3nytGIyMjJ7p4vC0qVOnEhMTw9KlS3Fycnpmv4mJCR9//DHOzs6ptt++fZv+/fvj4OCAubk5FSpUYMmSJfr9lpaWdO7cGV9fX+7evftMvStXrsTa2pr27du/1nMTeUNub7OVK1dm7969hIWFsX37dgICAqhZs+Yr63+ddvmqNgnSLkXmyu3t82nJyclcv379lfVL+3wzSYKQhz3+MqtcufJLyyUlJXH//v003bRabYbjOn36NGXKlMHGxibV9scnEn5+foDuyys+Pp6UlJRUfwM8fPiQnTt3kpCQQGJiIjNnzuTBgwfUqlXrpcf+66+/cHd3f2W5J4WGhlK7dm127drFkCFD+PHHH3F3d8fHx4dZs2bpy/Xu3Zvk5GTWrl2b6vEPHjxg+/btdOrUCUtLyzQfV+Q9ub3Nnj17lvj4eGJjY/nhhx/QarW0atXqlfWnt12mtU2CtEuReXJz+4yIiGDlypVER0eTnJzMunXr2L17Nw0aNHhl/dI+31CGvoQhDKdx48bK3t7+leV2796tgDTdAgMDX1nfqy6HVqhQQTVp0uSZ7RcuXFCAWrBggVJKqfHjxz9z/Md13r17V1WvXl3lz59fFSxYUDVu3FidPHnypXFFREQoQHXs2PGZfeHh4erevXv6W2xsrH6fj4+PcnJyUvfv30/1mJ49eypbW1t92eTkZOXk5PTMJd0FCxYoQG3fvv2l8QmR29vs8OHDla2trbK2tladO3fWd5F4mddpl2ltk0pJuxSZJze3z4iICNWoUSNla2urbGxsVLVq1dSGDRteeWxpn28uGaSch507dw5PT89XlvP09EzzDEJpveT5MnFxcfqp1p5kYWGh3w/w9ddfP3eBNAB7e3tOnDiRruNGRkYCkD9//mf2NWrUiDNnzujv//DDD3z66acopdiwYQPdu3dHKcX9+/f1ZVq2bMnq1as5deoUdevWxdjYmJ49ezJz5kyCgoJwdXUFdJdJHRwcaNq0abriFXlPbm+zM2fOZObMmemqO73t8pNPPklzmwSkXYpMk5vbp42NDbt370533dI+31ySIORRwcHB3L9/P019JQsWLEizZs2yISodS0tLEhISntkeHx+v358VHs9uFB0d/cy+hQsXEhUVRWhoKH369NFvv3fvHg8fPmTRokUsWrToufU+2Xeyd+/ezJw5k5UrV/LFF19w69Yt9u/fr58bXogXyattNr3tMr1tEqRdioyT9int800jCUIelda+kgCJiYk8ePAgTfXa29tnuME6OTk9d3GW4OBgAIoWLZqh+l/E1tYWJycnzp8//8y+x30rg4KCUm1/3D+0T58+9O3b97n1PvkaV69enXLlyrFq1Sq++OILVq1ahVKK3r17Z9KzEG+qvNpm09su09smQdqlyDhpn9I+3zSSIORR586dA9L2ZXbo0CEaN26cpnoDAwP1lwBfV5UqVdi9ezeRkZGpBlUdPXpUvz+rtG3bll9++YVjx46laXYVe3t7rK2tSUlJSfMvQr1792bs2LGcPXuWlStXUrp0aWrUqJHR0MUbLi+32fS0y9dpkyDtUmSMtE9pn28aSRDyqLNnz2JsbIyHh8cry2Z3f8muXbsybdo0Fi1axKeffgroVoFcunQptWrVemaK0cw0evRoVq5cSf/+/fH19cXBwSHVfvXUNKnGxsZ06dKFlStXcv78eSpWrJhq/71797C3t0+17fEX3bhx4/Dz83vhOAohnpSX22x62uXrtEmQdikyRtqntM83jUY9fcYj8oQqVaoQHx/P5cuXs+2Yc+fO5eHDh9y5c4f58+fTuXNnqlatCsDQoUOxtbXVl+3evTsbN25kxIgRuLu7s3z5co4dO4avr2+apl3LiM2bN9OrVy/MzMz0K0IqpQgMDGTlypXcvn2b33//nZ49ewK6Kdtq1arFvXv3+OCDD/Dw8ODBgwecOnWKXbt2PfdSct26dTl06BAA165dw93dPUufk8j98nqbTU+7fJ02CdIuxeuT9int841jgJmThIElJSUpMzMz1a1bt2w9rouLS5qncouLi1OffvqpcnR0VObm5qpGjRpq27Zt2Rbr9evX1aBBg5S7u7uysLBQlpaWqly5cmrgwIHKz8/vmfKhoaFq8ODBytnZWZmamipHR0fVtGlTtWjRoufWP2/ePAWomjVrZvVTEW8AabM66WmX6W2TSkm7FK9H2qeOtM83i1xBEEIIIYQQQujJSspCCCGEEEIIPUkQhBBCCCGEEHqSIAghhBBCCCH0JEEQQgghhBBC6OX4BOH27dv06dMHOzs7LC0tqVSpEidOnNDvV0oxbtw4nJycsLS0pFmzZly7ds2AEQshhBBCCJF75egEITw8nLp162Jqaso///zDxYsXmT59OgULFtSXmTp1KrNnz2bBggUcPXqUfPny0bJlS+Lj4w0YuRBCCCGEELlTjp7m9LPPPuPgwYPs37//ufuVUhQtWpRPPvlEvzpgREQEDg4OLFu2TL+QlRBCCCGEECJtcnSC4OHhQcuWLbl16xZ79+6lWLFifPTRR3zwwQcABAQEUKpUKU6fPk2VKlX0j2vYsCFVqlThxx9/fKZOrVZLUFAQpqamaDQa/XZzc3PMzc2z/DkJIf4vISGBhIQE/f2UlBQCAgIoWbIkxsbG+u3SPoUwDGmjQuRcT7dPpRRJSUm4urpiZJSxTkImGQ0uKwUEBDB//nxGjhzJF198wfHjx/n4448xMzOjb9++hISEAODg4JDqcQ4ODvp9TwsKCqJUqVJZHrsQQgghhBDZzd/fn5IlS2aojhydIGi1Wry8vJg0aRIAVatW5fz58yxYsIC+ffu+Vp2mpqYAHDt2DCcnJ/12+fXjxSIjI3F2dubmzZvY2NgYOpxcQ163V3v614/Tp0/Tvn17aZ/pJJ+19JPXLG2kjWacfNZej7xur/Z0+wwODqZmzZr6c92MyNEJgpOTEx4eHqm2lS9fng0bNgDg6OgIQGhoaKovqtDQ0FRdjp70uFuRk5MTxYsXz4Ko31w2NjbSSF+DvG5pFxwcDEj7fF3yWUs/ec3SR9ro65PP2uuR1y39nuxC/7py9CxGdevW5cqVK6m2Xb16FRcXFwDc3NxwdHTE19dXvz8yMpKjR4/i7e2drbEKIYQQQgjxJsjRVxBGjBhBnTp1mDRpEt27d+fYsWMsWrSIRYsWAboMafjw4UycOJHSpUvj5ubG2LFjKVq0KB07djRs8EIIIYQQQuRCOTpBqFGjBhs3buTzzz9nwoQJuLm5MWvWLHr37q0vM3r0aGJiYhgwYAAPHz6kXr16bNu2DQsLi+fW+biPpPSVTDtzc3PGjx8vr1k6yeuWfvnz5wekfaaXfNbST16z1yNtNP3ks/Z65HVLv8w8x83R05xmhcjISGxtbYmIiJA+bULkMLdu3cLZ2VnapxA5lLRRIXKuzDzHzdFjEIQQQgghhBDZSxIEIYQQQgghhJ4kCEIIIYQQQgg9SRCEEEIIIYQQepIgCCGEEEIIIfQkQRBCCCGEEELoSYIghBBCCCGE0JMEQQghhBBCCKEnCYIQQgghhBBCTxIEIYQQQgghhJ4kCEIIIYQQQgg9SRCEEEIIIYQQepIgCCGEEEIIIfQkQRBCCCGEEELomRg6ACGEEEIIIcTrSUlJ4dy5c+zcuTPT6pQEQQghhBBCiFwiPj6eY8eOceDAAfbv38+hQ4eIjIzM1GNIgiCEEEIIIUQOFR4ezqFDh9i/fz8HDhzg+PHjJCYmpipjbW1NzZo18fX1zZRjSoIghBBCCCFEDnHr1i19MrB//37Onz+PUipVGUdHR+rXr0+9evWoX78+lSpVIjY2Fltb20yJQRIEIYQQQgghDEApRWBgIHv27GHPnj3s37+foKCgZ8qVKVNGnwzUq1ePUqVKodFosiwuSRCEEEIIIYTIJjdu3GD37t36240bN1LtNzIyomrVqvpkoF69ejg4OGRrjJIgCCGEEEIIkUXu3LmTKiEICAhItd/ExISaNWvSuHFjGjZsSO3atbG2tjZQtI9iMujRhRBCCCGEeIOEhoayZ88efUJw9erVVPuNjY3x8vKicePGNGrUiLp165I/f34DRft8kiAIIYQQQgjxmqKioti7dy87duzA19eXixcvptqv0WioVq0ajRs3pnHjxtSrVw8bGxsDRZs2kiAIIYQQQgiRRikpKZw8eZIdO3awc+dODh06RHJycqoynp6e+oSgQYMGFChQwDDBviZJEIQQQgghhHiJoKAgdu7cqb9KEB4enmp/yZIlad68Oc2bN6dRo0bY2dkZKNLMIQmCEEIIIYQQT4iMjGT37t36qwTXrl1Ltd/W1pamTZvqk4JSpUoZKNKskaMThK+//ppvvvkm1bayZcty+fJlQLfU9CeffMLq1atJSEigZcuW/PTTT9k+FZQQQgghhMi9lFKcO3eOv//+m7///ptDhw6RkpKi329sbIy3tzfNmzenRYsWeHl5YWKSo0+jMyTHP7MKFSqwa9cu/f0n34wRI0awdetW1q1bh62tLUOGDKFz584cPHjQEKEKIYQQQohcIjo6Gl9fX7Zu3crff//N7du3U+0vU6aMPiFo1KhRjh9YnJlyfIJgYmKCo6PjM9sjIiJYvHgxK1eupEmTJgAsXbqU8uXLc+TIEWrXrp3doQohhBBCiBxKKcW1a9f4+++/2bp1K/v27SMxMVG/39LSkqZNm9KmTRtat26Nq6ur4YI1sByfIFy7do2iRYtiYWGBt7c3kydPpkSJEpw8eZKkpCSaNWumL1uuXDlKlCjB4cOHX5kgREZGprpvbm6Oubl5ljwHIcTzJSQkkJCQoL8fFRUFSPsUIqeQNipyu/j4ePbu3atPCvz9/VPtL1myJG3btqVt27Y0bNgQCwsLA0Wafk+3z6fbZUYYZVpNWaBWrVosW7aMbdu2MX/+fAIDA6lfvz5RUVGEhIRgZmb2zLRRDg4OhISEvLJuZ2dnbG1t9bfJkydn0bMQQrzI5MmTU7VDDw8PQNqnEDmFtFGRG4WFhfHrr7/SpUsXChcuTKtWrZg9ezb+/v6YmprSrFkzZsyYweXLl7l+/TqzZ8+mZcuWuSo5gGfbp7Ozc6bVrVFKqUyrLYs9fPgQFxcXZsyYgaWlJe+9916qzAnQL1X9/fffP7eOyMhIbG1tuXnzZqq+ZPLrhxDZ7+lfP27fvo2Hh4e0TyFyCGmjIrcIDAxk8+bNbN68mf3796caYFy0aFHatm1LmzZtaNq0KdbW1gaMNPM87wqCs7MzERERGR4vkeO7GD2pQIEClClThuvXr9O8eXMSExN5+PBhqqsIoaGhzx2z8DQbG5s8NdhEiJzo6ZOKx5dHpX0KkTNIGxU5lVKK06dPs2nTJjZv3szZs2dT7a9cuTIdOnSgQ4cOVKtWDY1GY6BIs05WJua5KkGIjo7G39+fd955h+rVq2Nqaoqvry9dunQB4MqVK9y4cQNvb28DRyqEEEIIITJTUlISe/fuZdOmTfz555/cvHlTv8/IyIgGDRrokwI3NzcDRprN4uPZs+AyK1dGZFqVOTpB+PTTT2nXrh0uLi7cuXOH8ePHY2xsTK9evbC1tcXHx4eRI0dSqFAhbGxsGDp0KN7e3jKDkRBCCCHEGyApKYl///2XdevWsWnTJsLCwvT7rKysaNWqFR06dKBt27a5fvXitFAKrl6F/fuhd2+wXDYfRozgeMJQfmZsph0nRycIt27dolevXoSFhWFvb0+9evU4cuQI9vb2AMycORMjIyO6dOmSaqE0IYQQQgiROyUmJuLr66tPCsLDw/X7ChcuTIcOHejYsSNNmzbF0tLSgJFmPaUg7L9oCl89BHv3wp69NLmyizthFri7QyMXF0hIoJXdCW4W/Ic51zPnuLlqkHJmeDxIOTMGcAghMtetW7cybYCVECLzSRsVWSUxMZGdO3eybt06Nm/ezMOHD/X7ihQpQufOnenWrRsNGjR4o1cwBiA2Fvbt49DvgXRZ0w3HpFucpqp+d/9Kxwko5MX48dC4VizcugWlSxMZFZVp57hv+CsshBBCCCFyoqSkJHbu3MmaNWvYvHkzERH/70Pv4OBAly5d6NatG/Xr18fY2NiAkWax2Fj+/iOeNTsL0akTdCwdCK1bUxIHQhhEJJZEFy9H/sY1oFEjljS1B5fHD7aCMmUyPSRJEIQQQgghRLbQarUcPnyYlStXsnbtWu7fv6/f5+TkpE8K6tat+8YmBXduJLPv10C6q7UY/bsLDh1in/s6fr3YHiMj6LjEA7y8cKxYkcPFt1KlT0Usyl7K1hglQRBCCCGEEFnq/PnzrFy5kpUrV/Lff//ptxcpUoTu3bvTvXt36tati5FRjl7D97VotWCkUbBgAcnbdlHmz+XEUBoP1lCZcwB0Uhsx+rw9bdsCGg0cPw6AoabdkQRBCCGEEEJkuhs3brBq1SpWrlyZap2C/Pnz07lzZ3r37k2TJk3ezDEFISEcXRPEJ+tqY2MDf/+tgV9+weTUKRrgw11jJyLqtIOeA6FZM2qVLk2tHLRUwxv4jgghhBBCCEOIjIxk3bp1LF++nP379+u3m5qa0qZNG95++23atWv3xs0+9N+VeLbPu45X2HaqnV0G589jY1KJg8lnMTeH+HiwGDQI7t1jS2N7jGtUBuOqr6zXUCRBEEIIIYQQr02r1bJ7926WLVvGhg0biIuLA0Cj0dCwYUN69+5Nly5dKFiwoIEjzTwxMZAvH/DnnzB/Pt/u7M7ilPcYyXaqcR40GspVMmVZjxAa9nDEwgJ4/30AcsPICkkQhBBCCCFEuvn7+7N8+XKWL1/OjRs39NvLlStHv3796N27N8WLFzdghJksNpZk3700/74pB46aERQExQIDYds22mDJFbMKlKtlD0PXQuPGaAoXpq+hY35NkiAIIYQQQog0iYqKYt26dSxbtixVFyJbW1t69epFv379qFmzJhpNDupQ/5piYxS+K+5wf99F3nswHfbuxSQ+nhjXUJKTi7B/P/Ts0AESE+ncujWdK1QATU1Dh50pJEEQQgghhBAvpJTixIkTLFq0iFWrVhETEwPouhC1aNGCfv360aFDhzdiXEFKChiH3IbJkzn+Rxjtg1dRGDPeZSfGaMHZmbk9D2L3fidKlQJwhVGjDBx15pMEQQghhBBCPCMyMpKVK1eyaNEiTp8+rd9epkwZ3nvvPfr06fNmdCG6f59lCxOY+nsx+vaFMT7mMH8+dbRGVGEMdYr9R/SHM7Ht0gzKl6fmG3B15FUkQRBCCCGEEHonTpxg4cKFqa4WmJub07VrVz788EPq1auXq7sQPXwIWxeH0DlpDZZ/b4CDB4l1+4FL/iP55x8YM6YwTJmCaZkynG5SCqyrGDrkbCcJghBCCCFEHhcdHc3KlStZuHAhp06d0m8vW7YsH374Ie+++y52dnYGjDATHD0KGzfiNWMQ/kku/Mku2qEbR9HJajt2v31MizaPTo3fwG5D6SEJghBCCCFEHuXv78+8efNYsmQJERERAJiZmemvFtSvXz/XXi0ID01k4VIzzp2D338HJk2CP/+kDUXxpSlJlb3g/RbQvj1OLi70MHTAOYgkCEIIIYQQeYhWq2Xnzp3MmTOHv//+G6UUAO7u7gwcOJC+fftSuHBhA0eZflotRIbEUuDwP7B+PWw5wNiEGyQna/j6ayjdsyfky8e0tk6YvVUcbMcbOuQcSxIEIYQQQog8IDIykuXLlzN37lyuXr2q3966dWuGDh1Ky5YtMTIyMmCErykqii3f+jForge1E/exPqUrAAWBUc1PUqKzF4ULA716Qa9emBk02NxBEgQhhBBCiDdYYGAgs2bNYunSpURFRQFgbW3Ne++9x+DBgylTpoyBI0wfrRYOHoQSJcDl5gFo3hyn+Arc5gQHqU2KaymMu3WGrl2ZVKM65M4eUgYlCYIQQgghxBvo+PHjTJs2jfXr16PVagHdKsdDhgzh3Xffxdra2sARplNEBGzaRP+lDVi+142vvoJvx1QBoHqpCLZ5LaXR8CoY17oGuXTcRE4hCYIQQgghxBtCq9WydetWpk2bxr59+/TbW7RowSeffELz5s1zzaBjpeDY/gTWTvuP8UzAZsd6SEigRYnP2GQ7WVcof364fBlNiRK0zCXPKzeQBEEIIYQQIpeLj49nxYoVTJ8+nStXrgBgYmLC22+/zSeffELlypUNHGE67d4Ny3+l34rPuKwtS1WgDwlQvjxdu1jR5fMUzK2MdWVdXAwa6ptIEgQhhBBCiFwqOjqahQsXMm3aNEJCQgCwsbHhww8/5OOPP841Kx1HRijmL9Bw8CBs3gyan39Gs2oV72HPGStvSnXy1q1NULkyZnKlIMtJgiCEEEIIkcs8fPiQuXPnMmvWLMLCwgAoXrw4I0aM4P3338fGxsbAEb6aUqC5egV++w2j3zYxIfQMsXFGHD0Ktfv3B1tbRr/dDurWhdw4u1IuJgmCEEIIIUQuce/ePWbNmsXcuXOJjIwEoFSpUnz++ee88847mJnl/Ek8Lx2LYvxHd0kJvMGGB00AyA+MbbyLwm+3oHx5wLYZNGtm0DjzMkkQhBBCCCFyuLt37/L999+zYMECYmNjAahQoQJffPEF3bt3x8Qk557SKQVJSWD28C4MH47xhnOsSzyHMS7cM3LAvrUX9OnDZ+3qQj5DRytAEgQhhBBCiBzrwYMH/PDDD8yePVufGFSvXp0vv/ySDh065PiFzf5eH8sXE61o0gRmTLaF7dspk/iA6Q5Tafy2E4VH+YGTo6HDFE+RBEEIIYQQIoeJiIhg5syZzJw5U9+VyMvLiwkTJtCqVascO1VpSgqkRERjtnENLF2K9no5zoT+QlgYTJtmjtGCBeDmxsjq1WWtghxMEgQhhBBCiBwiOjqauXPnMnXqVMLDwwHw9PRkwoQJtGvXLscmBgCLvrrBdz/m4/OkCQxMmA1AK80xFnz1Nd1GFNeNM+7WzbBBijTJ2delnjJlyhQ0Gg3Dhw/Xb4uPj2fw4MHY2dmRP39+unTpQmhoqOGCFEIIIYRIp+TkZBYuXIi7uzuff/454eHhlCtXjrVr13Lq1Cnat2+f45KDpCTd+AL27IFatYj5biY3ou1Yl9AeypSBKVMwuRXEh98Wp1AhQ0cr0iPXJAjHjx9n4cKFzyz0MWLECLZs2cK6devYu3cvd+7coXPnzgaKUgghhBAi7ZRSbN68mYoVKzJw4EBCQ0MpWbIkK1as4Pz583Tr1i1HjjP4fnIKLi663ACNBo4do4/JGlZ5/8jWbSZw+TKMGQNFixo6VPEact4n7jmio6Pp3bs3P//8MwULFtRvj4iIYPHixcyYMYMmTZpQvXp1li5dyqFDhzhy5IgBIxZCCCGEeLmjR4/SsGFDOnbsyJUrVyhcuDBz5szh0qVL9OnTB2NjY0OHqKcUEB8Py5dDnToErTxMcDD8/jvQoAH89BP2t/3oeWgYFi0byviCXC5XjEEYPHgwbdu2pVmzZkycOFG//eTJkyQlJdHsiXlyy5UrR4kSJTh8+DC1a9d+YZ2PB/w8Zm5ujrm5eeYHL4R4oYSEBBISEvT3o6KiAGmfQuQU0kazRlBQEKNHj2bdunUAWFhYMHLkSEaPHo2tra2Bo0tNKZgxPoJF8xLZqZpTIvwMAMMLjqLJyv107GqiSwYGDTJwpHnP0+3z6XaZETn+CsLq1as5deoUkydPfmZfSEgIZmZmFChQINV2BwcH/XLjL+Ls7Iytra3+9rz6hRBZa/LkyanaoYeHByDtU4icQtpo5oqNjeXrr7+mfPnyrFu3Do1Gw3vvvce1a9f47rvvclxywLFjaN7uxd8TT3L1gT2LwzuBszN89x1lL26kWy8TTE0NHWTe9XT7dHZ2zrS6c/QVhJs3bzJs2DB27tyJhYVFptf95DLk8suHENnv888/Z+TIkfr7t2/fxsPDQ9qnEDmEtNHMoZRiw4YNfPLJJ9y4cQOAxo0bM2vWrGfGVhpSSgps2QLLlsHq1WCxdi2sXs2XhPJ26eP0HO8JPQIgBy/Klpc83T4jIyMzLUnI0e/wyZMnuXv3LtWqVdNvS0lJYd++fcydO5ft27eTmJjIw4cPU11FCA0NxdHx5Ytu2NjYpPpyE0Jkv6e7JTy+PCrtU4icQdpoxp0/f56PP/6Y3bt3A1CiRAmmT59Oly5dctasRHfvouYu4OOFo7l514JVq+C9IUMgLIwmH38MVasaOkLxlKzs2pejE4SmTZty7ty5VNvee+89ypUrx5gxY3B2dsbU1BRfX1+6dOkCwJUrV7hx4wbe3t6GCFkIIYQQgtjYWL755humT59OSkoKFhYWjBkzhtGjR2NlZWXo8AAICIANC+/zacRYNMuXYRIfzxeV7fmv/yCaNweKu8LSpYYOUxhAjk4QrK2tqVixYqpt+fLlw87OTr/dx8eHkSNHUqhQIWxsbBg6dCje3t4vHaAshBBCCJFVtm/fzqBBgwgMDASgU6dOzJgxA1dXV8MG9oTI3Sep2LwicSmFqcNZ6hIPNWow8FM76G7o6ISh5egEIS1mzpyJkZERXbp0ISEhgZYtW/LTTz8ZOiwhhBBC5DGhoaGMGDGCVatWAbrB3PPmzaNdu3YGjkw3G9GlS+DhAfTvj83SpbzLfIJwxaKuF0z8DhrK9KRCJ9clCHv27El138LCgnnz5jFv3jzDBCSEEEKIPE0pxbJlyxg5ciQPHz7EyMiIjz/+mG+//Zb8+fMbOjwiwpJp0tyYc+c1BAZCsZo1YcUK5vU6ivHoelDxR0OHKHKYdE9zGh4ezoMHDwC4d+8ef/zxBxcuXMj0wIQQQgghcrrg4GDat29P//79efjwIVWrVuXYsWPMnDnToMmBUkBSEixZgm2tcuSLCcXUFE6cAPr2hcBAjH9dCk915RYC0pkg/PLLL1SvXh0vLy/mz59Pp06d8PX1pWfPnvzyyy9ZFaMQQgghRI6ilGLlypVUqFCBv/76CzMzM6ZMmcKxY8eoXr26weJKSIAfpydT0+0uCaUrgo8P+PuzyHI4//0HHToAlpZQvLjBYhQ5X7q6GM2ePZsLFy4QFxdHiRIlCAwMxN7enoiICBo2bMj777+fVXEKIYQQQuQI9+7dY+DAgfzxxx8AVKtWjeXLlz8zsUq2i49HLVjG96M7Eax1YDW16VvkIYwaRbmBA8HwvZ1ELpGuBMHExARLS0ssLS1xd3fH3t4eAFtb25w1l68QQgghRBbYtWsX77zzDiEhIZiYmDB27Fg+//xzTA20pHBsrG5xsx49gL59sVi7lqkcINbWiV5jvWDQfMgh06qK3CNdCYKxsTHx8fFYWFiwd+9e/fbo6OhMD0wIIYQQIqdISkpi3LhxfP/99yil8PDw4LfffqOqARcQi49Opmw5I27dNqJYMag3aBAcOkSfz+tA//5gYWGw2ETulq4xCLt27dKv2Hb8+HGUUoBuMZBFixZlfnRCCCGEEAYWGBhIgwYNmDJlCkopPvzwQ44fP26Q5ECpR/9Zvx6LGpVobX8CNzeIjkY3Tam/P3z0kSQHIkPSdQXB1tZW/3fLli0JDg6mSJEi+psQQgghxJtk06ZN9O3bl8jISAoUKMDPP/9M165dsz0OpWD1KsWUr6LZbt0Vx7M7APih2HtYXTmNaT4zQANmZtkem3jzpHua08ceXz0QQgghhHjTpKSk8OWXX9KpUyciIyOpU6cOfn5+BkkOADhxgh8HXuJsoDXTzzaD/Plh3DhsLxx6lBwIkXly3UJpQgghhBBZ6cGDB7z99tts374dgBEjRvD9999n+0DkkyehShUw/mUhmoEDmUp99hs3ZtigRBgXAI8mixEis732FQSA+fPn4+vrS3h4eGbFI4QQQghhMH5+fnh5ebF9+3YsLS35/fffmTFjRrYnBx99BF5esHQp0Lo1WFnRoI8LXwb4kH/OZEkORJbK0BWEuXPn8s0336DRaHB2dqZatWqpbo6OjpkVpxBCCCFEltqyZQu9evUiJiaGkiVLsnHjRipXrpy9QaSkwNKllDlXEI2mC9euAe+X0A0+lvMqkU0ylCBcuHCB5ORkTp8+zalTpzh16hQ///wzN2/eRKPR4OjoyO3btzMrViGEEEKILDF79mxGjBiBVqulWbNmrFmzhkKFCmXLsVNSYMkSqGFymio/vgdnzvARpjT+9QSe7zxKUCQ5ENnotROExwujFS1alKJFi9K2bVv9vrCwME6ePImfn1+GAxRCCCGEyCopKSmMGDGCOXPmAPD+++/z008/ZWuXojFDYpi+IB+NCceXM2gKFMBs3Dg8e5TLthiEeNJrJwgvm8XIzs6OFi1a0KJFi9etXgghhBA5gFKKBw8eEBgYyMmTJw0dTqaKi4ujZ8+e/PnnnwB8//33jBo1Sv8jaJZLToaffuLj3xayhm10YDPK5wM0UyZB4cLZE4MQz/HaCcK2bdtSrYsghBBCiNwpKiqKwMBAAgICCAwMJCgoKNX/o6OjDR1ipouKiqJ9+/bs2bMHc3NzVqxYQbdu3bL8uMnJMG8exMXBZ4PjYMoUSkQHE1CtK6Y/zYZatbI8BiFe5bUThAMHDmBnZ0f16tUzMx4hhBBCZDKtVsudO3cICAggICAAf3//VH/fu3fvlXU4OTnh5OTEqVOnsiHirBUWFkbr1q05fvw41tbW/PXXXzRo0CBbjr3rz1iGD7fCzAy6d7em5Ny5cPcuph98AMbG2RKDEK/y2gnCrVu3aN26NWZmZrRr14727dvTtGlTzGQFPyGEECLbxcXFPZMAPP5/YGAgCQkJL328nZ0dbm5u+purq6v+/y4uLlhaWnLr1i2cnZ2z6RlljTt37tC8eXMuXryInZ0d27dvz74fOzdvpuWQj+hd5wAN+rrh6gqU7Jw9xxYiHV47QViyZAlarZaDBw+yZcsWhg8fTnBwMM2bN6dDhw689dZb2Tb6XwghhMgL4uPj8ff35/r161y7dk1/u379Ojdv3nzpY01MTHBxcaFkyZKULFmSUqVK6f8uWbJknug2HBwcTKNGjbh27RrFihVjx44deHh4ZOkxN26EuTMS+Nu+H+YbV6MBfnPtDR8chOwa6yBEOmVomlMjIyPq169P/fr1mTp1KpcuXWLLli0sXLiQAQMGULNmTdq3b0+vXr0oVqxYZsUshBBCvLESEhIICAhIdfL/+O+bN2++dJIQW1tb/Yn/0/93dnbGxCRD/+znaqGhoTRp0oRr167h4uLC7t27cXNzy9JjRkcpBvWLIzTSink4MdLYGEaPhrFjJTkQOVqmflOUL1+e8uXLM3r0aO7evcuWLVv0MwN8+umnmXkoIYQQIle7f/8+ly9ffuYWGBiIVqt94eNsbGwoXbq0/ubu7q7/287OLvtm4MlF7t+/T7Nmzbh8+TLFixfPluSAkBDyf/ABP0WacpLqfFTlMCw9AVWqZO1xhcgEWfZTQpEiRfDx8cHHxyerDiGEEELkaCkpKQQFBelP/i9duqT/Oyws7IWPs7a2fm4CULp0aQoXLpztSYBSEBsaBQ8eZOtxM8PDhw9p0aIF58+fx8nJKUuTg+BgGDgQxoyBOgmX4K+/6GxmRudva8PI/ZCHr+CI3CXTP6k1atR47heXUgqNRsOxY8cy+5BCCCGEQUVHR3PlypVnrgZcvXqVxMTEFz7OxcWFcuXK6W/ly5enbNmyODg4ZP+VgKQkAgI1XA0wwdMTnM7vhF9+wfdMYTpdnUI5dYk/es/O3pgyKCEhgU6dOnH69GmKFCnCv//+i7u7e5Yd79sJij//1ODvD2fPNsZo1ixo0gQqVcqyYwqRFdKdIMTFxfHgwYNnxhRcuHCBChUqsH79+kwLTgghhMhJ4uLiuHTpEufOneP8+fP6261bt174GAsLC8qWLZsqEShXrhxlypTBysoqG6NHN/n+tWvcPHqH5Rvyox48YGyBueDvD//9R/8K99h7tiC//w5vx9+EtWspTGWisCYQN3jJVY+cRqvV0rdvX/bs2YO1tTXbt2+nXLksXJl4/34m7RnN7Wa7mDQrH0ZGwLBhWXc8IbJQuhKE9evXM3z4cAoXLoxWq+Xnn3+m1qMFPd555x1OnTqFg4MDCxYs4Pr161SqVAkfH588PShKCCFE7pOcnMz169dTJQLnzp3D39//heMDihQp8kwSUL58eUqUKIGRkVG2xK20inD/BxR6cF130u/vz7jggWw5bM/48dDx/u/wwQeE4clY/CjMPcbSQf/4inZ3CKtYUDcdf716MH065Z1LcdHkMq51ihKW9DPkkmlOR40axZo1azA1NWXjxo1UyYK+/3v2wL69WsaZTIZx4yig1bK52gCo8HumH0uI7JSuM/eJEydy8uRJHBwcOHnyJH379uWLL77g7bff1s+q8O6772JmZkb9+vX5559/uHTpErNmzcqK2IUQQogMUUpx8+bNZxKBS5cuvbBrkJ2dHZUqVaJSpUpUrFiRChUqUL58+Wyd2vtuqOLSZQ0lSoBb2AmYMYOLZ5OpeWEJ+UkihNr6sjfq9MDPz54LF6BjXXcoWJCSrgV47+EeSpZIRvvOEozKuEOpUsx1cgJ9z6YyMHIkZkD5x5tuRWbbc8yIOXPmMGPGDACWLl1K06ZNM/0Y165B06YKrdaI+uyiMVp4912YOzfTjyVEdktXgpCUlISDgwMA1atXZ9++fXTq1Inr16/r+0o+vvQK4OPjQ82aNTM5ZCGEECL94uPjOX/+PKdPn+b06dP4+flx/vx5oqKinls+X758VKhQgYoVK+qTgUqVKlGkSJHsGR8QF0f0uUA2rYzl9pVoxhRdAVevwtWrfOp+iBWHSjFpEnxeKxJWraIYNsSQn3gsiHFyJ1/polCqFB95RdL1c6haFSjaEB48wAZYkvXPwCD+/fdfRowYAcCUKVPo3bt3lhyn9M1/GWIRQFSsETUtz8NPS6Ffvyw5lhDZLV0JQpEiRTh79iyVK1cGoFChQuzcuZO+ffty9uxZAExNTf9feQa7Fs2fP5/58+cTFBQEQIUKFRg3bhytW7cGdF/2n3zyCatXryYhIYGWLVvy008/6ZMYIYQQeVNERAR+fn76ZOD06dNcvHiRlJSUZ8qamJhQrly5VIlAxYoVcXV1zfquQUpBaChcusTf18uw+UQxmjaF7g57oXFjElRB3kHX738wb5GfGADKlfenVKlSmJsDlSvD999jW6YMV82u4VK3OGa21/SHSP0z3Zs9BWpgYCDdunUjJSWFd955h9GjR2dq/Tt2gLc3WO/9C9q3Z6YCowoesHYvZPGCa0JkJ4162YorT7l16xYmJiY4Ojo+s+/gwYPUrVsXExMT/WVWpRQPHz6kYMGC+lmM7t69m+bgtmzZgrGxMaVLl0YpxfLly/nhhx84ffo0FSpUYNCgQWzdupVly5Zha2vLkCFDMDIy4uDBgy+sMzIyEltbWyIiIrCxsUlzLEKIrHfr1i2cnZ2lfYp0CQ4OTpUInD59moCAgOeWtbOzo2rVqlStWpUqVapQuXJlypQpg5mZWZbGqBRoYqJh925iz1yjz88NuXavAMfM6mEZEQLA+Aa7mbCvEe+/Dz9/GQRubmBjQ3ujLTgVTua7TicoXNUZypTR3aytszTm58nJbTQ6Opo6depw7tw5atSowd69e7G0tMy0+qdPh08/hd69YcXP8WgaNdTNTjR7NmT3YHMhniMzz3HT9RN/8eLFX7ivbt26gG5gV2Zp165dqvvfffcd8+fP58iRIxQvXpzFixezcuVKmjRpAuj6GZYvX54jR45Qu3bt51UphBAiFwsJCeHYsWMcO3aMkydPcvr0aUJDQ59btkSJEvpk4PGtePHiWdo9KOpBEla3r2F86TycP8/y8PZM+NuLNm1gzrAQaN8eS2A3D3hIQa7HFaaS0V1wc6N5+Vuk1IcGDYASJXRXFuzt+VMfb5Msizu3U0oxYMAAzp07h6OjIxs3bszU5ACgpksoJiZFKFhQg9bMAuPduyUxEG+sDE8vFBkZydKlSwkJCcHNzY0qVapQqVKlTG+YKSkprFu3jpiYGLy9vTl58iRJSUk0a9ZMX6ZcuXKUKFGCw4cPvzJBiIxMPdDK3Nwcc3PzTI1ZCPFyCQkJJCQk6O8/7gsu7VOA7hfhkydP6hOCY8eOcePGjWfKGRkZUbZs2VSJQJUqVbCzs8uy2GJjITwcitlEwY8/os6dp+LmiVxMcOcynSjLVV1s3k4EBHhx/jzg6gpeXmjc3VmQ4ksBd3tc26+GKqXAwoJ6QL3/PysoUiTL4k+r3NJGlyxZwqpVqzA2Nmb9+vXPTMX+upKTH61ttmMH9T/owYX3x1Fmjm58gyQHwtCebp9Pt8sMURnUtGlTVbhwYdW6dWtVoUIFZWJiokxMTFTZsmVV9+7dM1q9Onv2rMqXL58yNjZWtra2auvWrUoppX7//XdlZmb2TPkaNWqo0aNHv7C+iIgIBTxzGz9+fIZjFUKkz/jx45/bHqV95j2JiYnq1KlTasGCBap///6qYsWKysjI6JnPgkajURUqVFDvvfeemjdvnjp8+LCKiYnJsriSErXq0u5gFbV+m1JTpijVp4/6vesfSqNRqmVLpVRCglImJkqB8uKYAqX+tOimVO3aSvn4qJCfNqjdu5UKDc2yELNUbmij58+fV5aWlgpQU6ZMybR6165VqnRppYIn/qKUkZFSoJS3t+49FyIHeFH7jIiIyHDd6RqD8Dz58uVjz5491KhRA9BlM+fOncPPz48zZ84wZ86cjFRPYmIiN27cICIigvXr1/PLL7+wd+9e/Pz8eO+991JlTgA1a9akcePGfP/998+t73H/rJs3b6bqn2XoXz+EyIue/vXj9u3beHh4SPvMA+7cucPBgwc5dOgQx44d49SpU8THxz9TztnZmZo1a+pv1atXxzqL+t7fvw937kDlilpdZ3M/P6rvm8GplCpspQ1t+AeAw+X7U+fSYjw9wc8PGDMG7Oy4VrAmdjVLUahyccjuVZCzSE5vo7GxsdSoUYOLFy/SsmVL/v7770wZWJ6YCJ6eisuXNYzme77nM+jfH376CeS7SOQQz7uCkFljhDLcxahy5cqpZisyNzfHy8sLLy+vjFYNgJmZmX5Z9OrVq3P8+HF+/PFHevToQWJiIg8fPqRAgQL68qGhoc8dRP00GxubHDfASoi85umTiseXR6V9vlm0Wi2XLl3i4MGDHDhwgAMHDhAYGPhMOVtb21TJQI0aNXBycsr0eJSC/07cI7//GQr/dxLOnGHvQ08a/TMGd3e4ds0INm2CwEDKcIlLlOFu0arQwBYqVqR6peqE1n6iB9CjH6RKZ3qkhpfT2+hnn33GxYsXcXR05Ndff820WafMUuL4y3UUSy47MYFxMGUKjB79xiR+4s2QlYl5hhOEqVOnMm7cONavX58tvx5otVoSEhKoXr06pqam+Pr60qVLFwCuXLnCjRs38Pb2zvI4hBBCPF98fDzHjx/XJwSHDh0iPDw8VRmNRkPlypWpW7cu3t7e1KxZE3d390yfVjQ5GYKCwN0dmDkTdu+m104f1sR34Ef+5GN0V7kr2h5FoxmNkZGGxEQw++orMDJigVt5fvMywjjfd/o6zQDDjw4Qu3fv1vdSWLZsGUUyOGYjPBwuXoS6dRS0bUup3bv5ztQUlq2At9/OjJCFyDUynCC4uroSGRmJh4cHPXr0oHbt2lStWhXnTFiK/fPPP6d169aUKFGCqKgoVq5cyZ49e9i+fTu2trb4+PgwcuRIChUqhI2NDUOHDsXb21tmMBJCiGwUFRXFgQMH2L17NwcOHODkyZPPrEJsZWVFrVq1qFevHnXr1qV27drY2tpmahzxcYoY/xDsAo7DyZMEX46g5J+zSEmB6Ggw++cf2LmTMnhhSmvu25WD5j2hcmXsPD2JqKuwtn30C3H//gBkboQis0RFRdH/0Xs0YMAAWrZsmaH6QkKgSRO4eRMOHNDg+f77uv5jf/wBjRplPGAhcpkMJwhdunQhNDSUhg0bcujQIebPn09kZCSFChWiatWq7Nix47Xrvnv3Lu+++y7BwcHY2tpSuXJltm/fTvPmzQGYOXMmRkZGdOnSJdVCaUIIIbJObGwshw4dYvfu3fz7778cP378mQXIHB0dqVu3rj4hqFKlSqqFNDNKq9X19tDs3wc7dzJ1nRtfXnmHQaxjNsN0MQCWBWaQlGxEUBCU+eADaNuWT8vV5EuvRMztPgI+0teZ/asKiNc1atQogoKCcHV1Zdq0aRmur0ABcHKCyMhHvYjefhtatYJH6zoJkddkOEE4f/48hw8fxtPTU78tKCiI06dP61dXfl2LFy9+6X4LCwvmzZvHvHnzMnQcIYQQL5aQkMDRo0f1CcGRI0eeuULg5uZG48aNadiwIXXr1qVkyZKZu95ATAycOIE6fIQ2u0dx8LARZ86A2+rVMH8+RelNMv0JoJRu8apq1dBUr865plE4lbPFyAgo0w0Aw/ecFxmxf/9+Fi5cCOimN82MQesW/11hU+IIHm5cgnPlR+MYJTkQeViGE4QaNWoQExOTapurqyuurq506tQpo9ULIYTIZkopzpw5w7Zt2/D19eXgwYPExcWlKlOsWDGaNGlC48aNady4Ma6urpkZAFy7xolV1/hueTEKRQSxOKIrpKSgAcI8BhMVlZ+jR8GtZUuIj6edhzdBpU5RokVjyPf/H6cyZzZ8kVMkJiYyaNAgQNe1qHHjxq9d1759cO0a+NQ8B82aYX33LtZTBsOGDZkVrhC5VoYThGHDhvH111+zdu3aVLMJCSGEyD0ePHjAzp072bZtG9u2bSMkJCTVfnt7e31C0KRJE9zd3TPvCkFEBJu3mfPPbgt8fKDGgVkwciTJ1GITRyiCE4oUNMWLQ+3azGp9C+sa5ShfHjDpAB06YIuMF8gLZs6cyYULF7C3t2fy5MmvXc+VK7oeRPHxCmfrCbSIvAtVqsCCBZkXrBC5WIYThK5duwJQunRpOnXqRK1atahatSoVK1bEzMwswwEKIYTIfCkpKZw4cUKfEBw7dgytVqvfb2VlRZMmTWjevDlNmzbFw8MjUxKClBQ4u+suF/70p49aofsZ9+JFVtUKZM0RF4oXhxoNvcDcnCrVrJhmtoVaTfNDv5vgXByAOhmOQuRG//33HxMmTABg2rRpFMpAF6AyZaB/l4f4rztF/ci/wMsLduyAggUzK1whcrUMJwiBgYGcOXNGvzDapEmTCAoKwsTEhLJly2Z4HIIQQojMERERwd9//82WLVvYsWMHYWFhqfZXrFiRVq1a0apVK+rVq5cpU1dHR+tujnGB8O233Pv3ItX+O4KGwrSjFbbo5tXvVvwwxT9xoXFjoJY3REZiYWbGJxmOQLwpRo8eTWxsLA0bNuSdd97JUF2amzeYs7c+SQnBmFUqB9u3S3IgxBMynCC4uLjg4uJC+/bt9duioqLw8/OT5EAIIQzs1q1b/Pnnn2zatIk9e/aQlJSk32dra0vz5s1p1aoVLVu2pHjx4hk/oFJw+TLs28fM4/UYtawCffvC4glmsHQpjkB1TmBvk0BYl1HYtqsAdevSpUgRuugryfA/TeINc+jQIdauXYtGo2H27NmvdTXr0CHYsgUmTQLNgAFobt7ArEwZ2LlTBiQL8ZQs+Ra2tramfv361K9fPyuqF0K8IRISEggNDSUkJISQkBAuXbpk6JByPaUUFy5cYPPmzWzatIkTJ06k2l++fHk6dOhA27ZtqV27NiYmGfxnQCm4fp3Jn4ax9aAtP2vfp3z4IQDKVh9LSsoEAgOBYsVg8mTw9OR4nTJobG2Auhk7tsgTtFotI0aMAMDHx4fKlSunu46wMHjrLd1iaC4uMHDxYvjgA1i4EBwcMjtkIXI9+ZlGCJGpkpOTuXfvXqoT/8e3p7c9fPjQ0OG+EZRS+Pn5sXr1ajZs2IC/v79+n0ajwdvbm44dO9KhQwfKlCmToWMlJcGJExAYCG/31IKHB1y5wr/s4CC1+ZcqlLc4Bd7eNGlekP/+gBIlHj34s890MWUoApHXrF69mmPHjpE/f36+/fbb16rDzg5++AF+/RXefRewKgZ//525gQrxBpEEQQiRJikpKdy9e5c7d+689Hbv3j2UUmmu19TUFEdHRxwdHbGxscHX1zcLn8Wb5fLly6xevZrVq1dz5coV/XZzc3OaNWtGx44dadeuHQ4Z+IVUKUiKjMPs+EHYuZPLZxV1tk3Fygq6djXCzNkZAgIYWmY3PUsn0Krfe9B6JpiZYQGUeOURhHixxMREvvrqKwA+++wzHB0dX68ipfDx+5j3BtbDyKpHJkYoxJtJEgQh8jitVktYWNgrT/xDQkJSzXLzMkZGRtjb2+tP/B/fHBwcnrlfsGBBfX/iW7du4ezsnJVPN9cLCgpizZo1rFq1ijNnzui3W1hY0LZtW3r06EHr1q3Jnz//6x/k0TiCX78J5Js/q9Aj8TcmpYwBoAIaypWeRIXKJoSHg8PPP4O9Pe3z5cvoUxPiGUuWLCEwMBAHBwd9N6O0io2Fb7+FcePActYUmDsXo/nzoVYNKFkyiyIW4s2QqQnC8ePHqVatGsbGxplZrRDiNSUnJ3Pnzh1u3rypv926dYvbt2+nOvl/cuDqyxgZGeHo6EjRokWfe3NycsLJyYnChQun73sgMRFCQ+GJE17xf1FRUaxZs4alS5dy6NAh/XYTExNatmxJz549ad++PTY2r79G8KVL8M8/uu4Xhb8dDrNnY0RvAmjDThozqWhRaN4co+bNudg+Ho314wTENUPPTYgXiYuL03cp+vLLL7GyskrX4wcMgN9/hyu77/DH0S90G2fNkuRAiDTI1ARh5MiRXLp0CQ8PD+rUqYO3tzd16tTB3t4+Mw8jhEDX5SckJCTVif+TicDNmzfT9at/kSJFXnji//hWpEiRtJ/4KwWRkXDnDkd2ReN/MYGGBc9SPP469OzJgcSafPghOJuFss3vNbsNvMGUUhw8eJAlS5awdu1a/Yr1Go2Gxo0b07NnTzp37oydnd1r1g+akGBdP+ytW+l54XfOXrXE0RHerl0bFiygTd14NpXcQpMPS4PXLXh0pUfGEIjssHDhQu7cuYOzszMDBgxI9+N9fMB3ezLDz/bXbfj4YxgyJJOjFOLNlKkJwv79+0lOTsbPz48jR46wbt06Ro0ahVKKOnXqsHz58sw8nBBvtIcPHxIUFERQUBA3btx45irAnTt3SE5OfmU9pqamFC9enOLFi+Ps7IyzszPFixenWLFi+hN/BweHtC9sqBRERHBhXxhB56KoZuuPU6w/tGrF8YTKDBsGjoTwh19JiIsDYBT7OEB91vIj3VgP7u6YVKvJxYsQ4/ho/VsTE91IwtDQ133J3gghISH8+uuvLFmyJNW4grJly+Lj40OfPn1wcnJ6vcqVInj3ZUYOS+a8vyVn40rrT/Y7NTmBo2t93WyPjTpBu3YUyp+fDhl+RkKkX3x8PN9//z0AY8eOfa01ORp7ReFv1xCr+6ehQQOYNi2zwxTijZXpYxBMTEzw8vKicOHCFCpUiEKFCrFz506uXbuW2YcSItdSSqVKAIKCgvjvv/9S3Y+IiHhlPcbGxhQtWlR/4v/45P/J+0WKFMHIyChNcUUHR3Hr1F0ck29RIOom1KjBmfiyjB8PBeODWXqtHgQHQ1wcg9jLfhqwhkl0Zx3ky4eqUZnDh8G5SAF9coCtLTWM/TFV+bGuWgOqu4GnJxUqwK5dULSwERS/r1uk6M4dyINjEJRS7N+/nzlz5rBx40ZSUlIAyJcvH927d8fHx4c6deq81tzvd+7opnasYHIF2renwNUb/Ml9YsnHaapSrYYJvPUWX3crDOUfP8oi856cEK/h119/JSQkBGdnZ/r27Zvmx12+DIULQ2E7Be+9h9WV01C0KKxdC6amWRixEG+WTE0Qpk6dypEjR7h79y5ubm7UqlWLPn36MH36dEylYYo8Jjw8PNUJ/9O3yMjIV9Zhb2+Pq6srJUqUSHXS/zgRcHJySlOXH60WtBFRmATfhMKFCYotwuLFYHbvFmMDfeDWLbh1i9aRf3GA+qzhc91J/5w5JNQsy+bNULxIIbgboK+zgtl1YozsMC9TDjzfAXd3ypWDDRugqB1Q/Do4OYGVFTP0j6qq/8saaNoUwAx4vW4yuV1sbCwrV65kzpw5qRaW9Pb2xsfHh+7du2NtbZ3+ipWC48dZvsaCfjMq06IFbN/sAsHBWJorFngsonSb0lT56G8oKt27RM6SkpLCDz/8AOi6Lqf16ubDh7q1DhIT4e8tWiq6uYGZme5LSdY6ECJdMjVBWLFiBebm5rRp0wZvb29q165NQVm6XLyhlFIEBwfj7++Pv78/169fT/V3eHj4K+soUqQIrq6uz9wer1CeL60zw2i1PIw0Yt8+SLgZSrf7C+DmTbh1ix6Hh7MpsjG/8j49WAtz5nCv1hAmToRi9vaMvbdDX01xbmFDBDFFS4NHM3BwoEwZmD8fShROAcf9ul/jnJyYb2n56FET9I+3ATp3Bt0v0KXS+lLmOffu3WPu3LnMnTuXBw8eAGBpack777zDkCFDqFSpUrrrDA1RrJ/+H/VC1uO5bw7cuEEdl2bATuLiIMXUAuNt26BSJd55naRDiGyyceNGrl+/TsGCBXn//ffT/Lh793TDZIyMoKizsW7hg48+Aje3LIxWiDdTpiYI586dIyoqiqNHj3LkyBHmz59PWFgY7u7ueHt7M3DgwMw8nBBZLjk5mRs3bjxz8u/v709AQACxsbEvfbyDgwMuLi4vTALSMitHcjI8uBlDkbvnISgI/vuP6ZtK8e91Z4ZZ/kyLh2th4kSu1R5Khw5QtHBBut3/Wv94I/qSiDm3KQYFCkBiIm5uMGgQuBRJAdelULw4ODvza5FimBbMD3ynf3wBQNd0rYB66X4Nxf8FBAQwY8YMlixZQtyjLlhubm4MHjyY/v37v94PKkFB8NtvjJpakRVRHRmCOXO4AfnyUbq2HSH/xuBQ8lGiWadO5j0ZIbKAUko/9mDIkCHpmq63dGk4cSCeO3dNKFTo0emNJAdCvJZMH4NgbW1Ns2bNqF69OtWqVePw4cOsXLmStWvXSoIgciSlFCEhIVy5ciXV7erVqwQFBb10ILCRkREuLi6UKlWKUqVK4e7urv+7ZMmSafvHTauF4GCO/HWfM8cSqJPvDJXij0O7dvg5t8PLCxwKmnD7fm39Q07yO39TkyasoQWRcOsWLi7g5QWuTlqUwwdoSjhD8eJMMS3OFKfrOFWdAIV0nX0KAz/9BLqT/n76eqUjYNbw9/fn22+/5bffftOPL/Dy8mLMmDF06tQpXVPCKgUbN+pWhF2wABy//hqWL6cXrbisKU6VGubw+UZo2RIsLZGOFSI3+ffffzlx4gSWlpYMHTo0TY9RSj/BFrYTPsH2zBn47Tdwdc26QIV4w2VqgrBo0SKOHDnCoUOHiIuLo3bt2tSuXZvff/+datWqZeahhEi32NhYrl279txE4GXjAczNzSlZsmSqk//Hf7u4uLyyf6zSKpJDwzC94Q92dty2dOerryD6TgTr/KvrugIlJvIjK1lNL6axhkr8DAUK4PRJO1JSICzKjKRirpi6FQdXV95LSaCJ6X7qNugCdT8AZ2eK5IPjx0HXvWeR/vgumfPyidcQEBDAxIkT+fXXX/WJQcuWLRkzZgyNGjVK36DjlBTYvh3N8uVMPr+EExfz0bQpDO3bVzd+5J2etO5cFqy9sujZCJH1pk+fDkD//v3TNEV6TAy0agWffgodTLY+/uUDrl2TBEGIDMjUBOHy5cu0bt2aCRMmULx48cysWog0Cw8P5+LFi1y4cIELFy5w6dIlrly5wo0bN174GCMjI1xdXSlbtmyqW+nSpSlatOgrZwFKToYrV+D2lWhahK0Cf3/w9+eTve2Zf68LE5nESGbC6NGYjPyeZctAo7EhUd3AjCQwNqau7WWiNYcoUbEkNBoPjRpRpAjcvg0ODhqMjQP1x2ueWS+WyBL37t1jwoQJLFiwQH8FqnXr1nz99dfUrFkzzfXEx8OvP4bz15K7/BHTEpPb/wEwvNVbnGv3Di1bAmUaQ+PGWfE0hMhW/v7+bNu2DYDhw4en6TGzZsGBAxAYoKW59mOsAEaMgObyLSlERmRqgjBjxoxXFxIik0REROiTgCdvwcHBL3xMwYIFn0kCypYti7u7+6vn2U5MhIAAjm69z9694GV8miaJ26B5c+73HE7FimBklI847WDdST9gTlXisMKfUlCsGJibU6QITJoEriUUysEXyrhA0aIMMTFBt4TP//uJa9CNCRa5Q1xcHD/++COTJ0/WX5Vq0aIF33zzDbVr137Fo5+g1cI//2A0fzFfbl3IfcqynQq0LRQFffrQu39l8MyiJyGEgSxcuBClFC1btsTd3T1Nj/n0U4iIgBaX5mD1V4BuIMJ33736gUKIl8rUBOHEiRMEBATQvXt3AD788EP9TC4jRozA29s7Mw8n8oiEhAQuXLjAmTNnOHfunD4RuH379gsf4+zsjIeHBxUqVMDDw0OfCBQuXPjl3TqUgpAQSEwk1t6FTz8F/0sJbL1VBZPAa5CSwh9MYSpjGMJpmvA35MuHw7DhFC0Kjo4QbtcDhzK2UKoUgwtWpL/DNUrU6Q+2gwHdSf/nnwMYAfUz86USBqKU4s8//2TYsGH895/uV/6qVasybdo0mjRpkqY6kpJ0szEeOwYzJiVCv36Y3b/P57iidXOnxqj+8N4GsJA1CsSbJy4ujsWLFwPw0Ucfpflx5uYw9a198MNw3Yaffwb9DGtCiNeVqQnCF198wbQnVio8cOAAP/30E7GxsUycOJGtW7dm5uHEG+ju3bucOXMGPz8/zpw5w5kzZ7h8+fILBwoXK1aMChUqpLp5eHhgY2Pz0uM8fKDFJugsRtevwpUr/LzViR/8mtNRbWJq4nDo2ROL31exdCnEx5sTZKbBPSUF8uenrn0wvZL3U6uyPbRfCFWqoNHougLpTv9X/D++THtlRE4VEBDAxx9/rP9+K168OJMmTaJ3795pXqAOf3+C5/5J79nD0Wo1+PhYUGHUKAgJYeSAAVCuXBY+AyEMb926dTx48IASJUrQtm3bV5Y/eRKqVQNNQjx88IFu4wcfQMOGWRypEHlDpiYI4eHhVK5cWX/fw8ODho8a68SJEzPzUCKX02q1+Pv7c/LkSX0y4OfnR0hIyHPLFypUCE9PTypXrpwqEShQoMALj6EUhAbEEHfuOm4PT4OlJdpuPShZEv77z4gg4y64pOgW/tIygGu8zwXcdZNox8RgZASTJ4OtLRQq9itU1M3/316joX1WvCgiV0lJSWHmzJmMHTuW+Ph4TE1N+eSTT/jqq6/StH7FpUtw8o//6HNuDKxbRwmtlg/bdMahpgtFigCjR2f9kxAih/jp0eDiDz/88JWzeu3dC40aQadOsHbybUzMzXWLMk6dmg2RCpE3ZHqC8KR169bp/37RiZ948ymlCAoK4sSJE/rbyZMniYiIeKasRqPB3d0dT09PqlSpgqenJ56enhQvXvylXYNCQ+HsWagdtBrra6fgwgUWHqnCoAff0YFANvEe1KiBUY8ePD53u1b2LVxsj0PZsrRzqkQZi/2Ua1gaasfqrlsD/x8nJzPDiP+7cuUK7733HocPHwagSZMmzJs3j3Jp/KX/5C+nqfGBJ1bY8RbbKIAWWrXip/GhUFPmnRJ5i5+fH0ePHsXU1BQfH59Xlr94EUxNwd4eTMqWglOndOuBvOQHIyFE+mRqglClShWWL19O3759U21fsWIFnp4yoi6vuHPnDseOHUuVEISFhT1TzsLCAk9PT6pWrapPBCpVqvTStQPiY1K4sP0WUWcCaJS8Szf145Qp1K2rmzjI1+kvmgT/DkBp4tGgJcasINRrAjVqALB5MxQpAjY2P+rrLfroJsTLKKWYN28eo0aNIj4+Hmtra2bMmIGPj88rpyyNjASbpDDo3p1q//5LFU7iyn9Et+9NgQkDQL4jRR61dOlSADp27IiDw6tX7hg0COrXB2fnRxtMTCCNg5qFEGmTqQnCjz/+SMeOHfn111+pUqUKoPtlICIigs2bN2fmoUQOkZSUxJkzZzh8+DCHDh3i8OHD+kGaTzI1NaVy5cp4eXlRo0YNvLy88PDwwNT0xUtz3bgBp09D1cA/KHF5B5w5w7+nnWib8AcViOI8k6BgQZg8mSpVNBgZQWy1llC4AFSoQIOyFYkpE45l8YaAr75e+XdEvI4HDx7g4+PDpk2bAGjevDm//PILJUqUeOnjbt7UrUQdEADnzxbE+MEDNKamHHx3KZafD4dSnbI+eCFyqMTERFauXAlAv3790vy4ihEH4c+9MHKkDNwXIgu8doIwYMAAZsyYkerX3mLFinH8+HF27drFpUuXAN3c302bNk3fgkAix7p37x6HDx/WJwTHjx8nLi4uVRkjIyMqVqyoTwS8vLyoVKnSC6cR1Woh4HgYQbuu00zt1GUGixbh4wO7dsFCjxsMuLgQgEoUpzD3KG4bheo+AI1nZUhJYe1aE3TjQd95dNOtCiwrA4vMcOzYMbp168aNGzcwNTXlhx9+4OOPP37191pgILYTZ3PkyAwiIjQcPW5EnSVLwM4Oy1ckFkLkBf/88w/379/H0dGRFi1avLTs7NnQuTMUL6qFjz/WdS2KjIQpU7IpWiHyEPWajIyMVGho6Os+PE0mTZqkvLy8VP78+ZW9vb3q0KGDunz5cqoycXFx6qOPPlKFChVS+fLlU507d1YhISEvrDMiIkIBKiIiIktjf1PcvXtXrVu3Tg0ePFhVqFBBAc/cChYsqNq0aaO+/fZb5evrqyIjI19YX3KyUhcvKhW06bRSkycr1bmzuujURIFSVkSrZIyUAqXu3VNjxyrl6anUcp+9Sn32mVKrVyvtxUtKm5Scbc9fZK+bN2/muPa5fPlyZW5urgDl7u6uTpw48dLyDx4o9duiGKVGj1bKzEwpUFt6r1JXr2ZTwEJkocxuo506dVKA+vTTT19aztdX90+DjY1SDxb/8f87WXweIkRukpnnuK99BUEpleHk5FX27t3L4MGDqVGjBsnJyXzxxRe0aNGCixcv6mcJGTFiBFu3bmXdunXY2toyZMgQOnfuzMGDB7M8vjfRvXv32Lt3L3v27GHPnj1cuHDhmTIeHh54e3tTp04dvL29KVu27Aunc7wbGEORmyfhxAkYNIiPP7Xkp59gdJW7fO/3OQBlMMKO+5S0uMP9Nh/iUKcUGBkxYQJMmADQ4NFNN4moENkhJSWFUaNGMXPmTADat2/PihUrXjqFbvgDRRnXBO5HWVGGf6lBIjRtylujykPp7IpciNwhLCyMv/76C+CZsYtPK1hQN+7As5KWgt9/ptv4ySe6AWVCiEyXqWMQMtvjJdcfW7ZsGUWKFOHkyZM0aNCAiIgIFi9ezMqVK/WLES1dupTy5ctz5MiR9K1cmkclJiZy8OBBtm/fzrZt2zhz5swzZSpXrkyjRo1o1KgR9evXp3Dhws+tS2kVmps34NAhkvYfodziUQQkFuc2PSlKMNSpQ/XqtbGygviiJaFsD/DywtjLi1BPM4wLVgZ+yuJnLMSrxcXF0adPH/744w8Axo0bx/jx41++rsHlyxQcNIjWUe9ximoklygFP30NbdqAdLEU4hmrVq0iKSmJatWqUbFixZeWrVpVN71p0i8r4KerULgwjBiRTZEKkfe8doKg0WiyfVzB42kxCxUqBMDJkydJSkqiWbNm+jLlypWjRIkSHD58+KUJQmRkZKr75ubmL+wj/6bx9/dn27ZtbN++nX///ZeYmJhU+9OaEJCYCEqxc585X34JpbVX+P1keUDX9z8//dFQlPN2jSjaIB5MTOjdG/r2BWNjd2C1vqqXz3ot3lQJCQkkJCTo70dFRQGGbZ/h4eG0b9+eAwcOYGZmxooVK/Srwz8tKkq3VsaYMWD75ZewZw9zLC6Rb+xITD79FczMsiVmIbJKVrbRVatWAfDuu++mqbwmMQGzieN0dz77DKyt03U8Id40T7fPp9tlRmSoi1G/fv1e+YXw+Be4jNJqtQwfPpy6devqf2kICQnBzMzsmcWyHBwcXrnugrN+fjSd8ePH8/XXX2dKrDmNVqvlyJEjbNq0iU2bNnHt2rVU+x0cHGjRogWtWrWiefPm2NvbP7+ihAQ2T7/Otk1xDGQhnud+gyVLMHXqxfHjEFLEFYyNdT/11KnDWtfbOLW0w8Zjpb6KvJGCibSaPHky33zzzTPbDdU+7927R7NmzTh79iy2trZs3rxZv9jj83ToALt36xKFOTNmgEaD7Q8/gJtblscqRHbIqjZ68+ZNDh06hEajoVu3bi8st2QJxMbCgAFgtnixbhKLokXho4/SfCwh3lQvap+Z4bUThFf1F8xsgwcP5vz58xw4cCBT6rt582aqvsRv2tWDxMRE/v33XzZt2sTmzZtTJUwmJibUrVuXVq1a0bJlSzw9PZ/bdeLOHTizJ5zWAfNgzx44dIhlcb+ziU6UZA2exMPJk9T6the//QZ1a5lA0UiwsgKgbHY9WZFrff7554wcOVJ///bt23h4eBikfd69e5emTZty/vx5HB0d2bFjB5UqVXp+Ya0WZs1inK0xQW7DaN8ecHGB9euzPE4hslNWtdH1j9pK/fr1KVr0+avQxMTA55/D3bu6Fe3fqVsXunfXLaNsaZn+JyPEG+bp9hkZGflM8v66XjtBeLywSXYYMmQIf/31F/v27aN48eL67Y6OjiQmJvLw4cNUVxFCQ0NxdHR8aZ02NjYvHWyYGyUnJ/Pvv//y22+/sXnz5lSXmmxsbGjbti0dO3akVatWz33uSqvQXL4ECQn8V6gqrq5galqAB0lTyI+uG1IP638o6WRC/Q7e0O8ClC+PpQZ69wbdxylHD2sROczT3RIef2azu32GhYXRpEkTLly4gJOTE7t376Zs2WdT3MuXITwoAu+5vWHrVhoBV/bXxrRerWyLVYjslFVtdM2aNQAv7L4HuvXPxo2DtWuhZ0/A1BMePU4IkbXdb3P02ZxSiqFDh7Jx40b27NmD21OX7atXr46pqSm+vr506dIFgCtXrnDjxg28vb0NEXK2U0px6tQpfvvtN1atWkVoaKh+n6OjIx06dKBTp040btwYs+f1h46I4M8pF5m61J5aMf8yPfpDaNOGEn9tpWRJsLPTcKfMaMrUKQyNGtGzfHl6yoBL8QaJjo6mTZs2XLhwgaJFi7Jnzx5Kl352yqGDB6Fl8xRsEuLx0x6jiLk5zJqFad2aBohaiNwrKCiIo0ePotFo9P92P4+5OQwerLsJIbJXjk4QBg8ezMqVK9m8eTPW1tb6bjK2trZYWlpia2uLj48PI0eOpFChQtjY2DB06FC8vb3f+BmM7t69y9KlS1m2bBmXL1/Wb7ezs6NHjx68/fbbeHt7P9N1KCZGt/hYnTPzsd+1Cg4dIjmlPQf5gzCSmG5uDmZmaFBcvKhBl5iOy94nJ0Q2SUxMpEuXLhw7doxChQqxa9eu5yYHAFXOrcAlrjqOBKNcS8LGBfBoxXghRNqtW7cOgIYNG77yaj8A27fDH3/AqFHg7p7F0QkhgNdfKC078JxFuQC1dOlSfZnHC6UVLFhQWVlZqU6dOqng4OAX1pmbF0rTarVq//796u2331ZmZmb618PCwkL16NFDbdmyRSUkJDzvgUpduqSUUqpePd36MovK/KD7A1Rk6WpqUdNVKmDJbqViYrL1OQnxpOxcKE2r1SofHx8FqHz58qkjR448UyYu7tEfY8cqBeoOjiq5czelcuH3hxCZITPaqLe3twLUvHnznrs/IUGp3r11i6NptUqpRo10/1598slrH1OIvCBHLJSWHVQaFmOzsLBg3rx5zJs3LxsiMoy4uDiWL1/OvHnzOH/+vH57jRo1GDhwIF27dn2mL2hcrGLp2AB2bYpmXWIHjG/fgOBgWrZ04PZt0DRvBsPmQevWWLu58UF2PykhDGzmzJksXrwYIyMj1q5dS61aqccRXL2qW8Jg8mToVqUKaDQ4jR0A48fDy9ZDEEK8UGhoKEeOHAF0iw8+z6pV8PvvuhnCAv+6gNmePboZ8oYNy8ZIhcjbcnSCkNdFRETw008/MXPmTO7duweApaUlb7/9NoMGDaJ69eqpymu1YHT1Mvz+O8a/reWLoKNEUIqDlKCB1T04e5YxY5rz5Zeg0VQBqmT7cxIiJ9i5cyeffvopANOnT6dNmzbPlFm8GPz9dat5d/brjPEF3aB8IcTr27p1K0opqlevnmrSkSfVr68bd1C6NJgtmqvb2LEjZNLsLEKIV5MEIQe6e/cus2bNYt68efoZI1xcXBg+fDj9+vV7Zt2HK1d0U8HFBQbzj5/uBMYMGG32I+aVSlNm6JfQvQFYWmKazc9FiJwmODiYPn36oJTCx8eHYc/7VfLePSZdGYjJkKUMG2uDsTGSHAiRCf7880/gxVcPAEqWhLlzgYcPodivuo1DhmR9cEIIPUkQcpDo6GimTZvGtGnT9Ksbe3h48Nlnn9GzZ09MTf9/ep+cpDA5tA8AC9eGbNoESjlx26IUxZqWgz59+KJ9e/2aBEIISElJoXfv3ty9e5fKlSszZ86cVCvCX7sG7gXD0DRrhvHZs3zX4D7Y7wFk5i4hMiouLo4dO3YAL08Q9JYv162SVqECvGTBQiFE5pMEIQdITk5myZIljB8/Xj9TU40aNfjyyy9p165dqpmIDu+M5svB4ZS9f5D54b2gVi1cjhxh7lzdZdlibn6QP7+BnokQOdvEiRPZvXs3+fLlY+3atVg+sdjSP/9Ax46Kzwpt4OuQs2gcHWHhQpBpfYXIFL6+vsTFxVGiRAk8PT2f2X/uHCxdqhtq4FJCwfz5uh2DB0s7FCKbSYJgYHv37uWjjz7i4sWLAJQqVYopU6bQpUuXVL9scvMmTJ9OzMLr7I7/izM0Z7aVLaaenpCYyEcfPV7jQJIDIZ7n2LFjTJgwAYAFCxY8sxDauZOJJCaacS6kMNrCDhj7+kK5coYIVYg30pPdizTPOeH/8Ufd2J87d2D1sgTduIMNGx6vxCmEyEaSIBhIeHg4o0eP5pdffgF06xeMGzeOgQMH6hc0u3cPfvgBygX+Q//NHSApiabA1CLT6DmkMKYf/6dbf14I8VIJCQn0798frVZL79696dOnT+oCKSmMPtWTSsTRrMBJjH13gYeHYYIV4g2k1WrZsmUL8OLuRZ0765KDIUMACwuYMkU3jZhcPRAi20mCYAAbNmxgyJAh+u5EH374IVOmTHlm8PGKFboEobRTffolJWPUqBGazz9nVPPm8oUpRDpMmjSJCxcuUKRIEX788Uf99rAwKFgQjCZOhI0baW1mBn/ugsqVDRitEG8ePz8/QkJCyJ8/Pw1fMJ6gTRvdLRX5t04Ig5AEIRvFxMQwdOhQli5dCkDZsmX5+eefqV+/PqCbsCH8wh3cFn8FlSszcOBwfH1h0IdW4HgManoZMHohcqezZ88yadIkAObOnYudnR0AUVHQuDGUKgUrvv+A/Js3w5gxusE8QohMtX37dgCaNGmiv0r+Qvv2QUwMtGihW/9ACJHtJEHIJufPn6d79+5cunQJIyMjPvvsM8aOHYuFhQUA2zfH0+ftFGrEn+Vv7VIoVAirAQPYutUKMAIkORAivZRSDB48mOTkZDp16kTXrl31+44f100RfO8ePLQqSv5jx8BEvhKFyAqPE4SWLVs+s8/PD06cgJ49H82xMWEC+PrC1KkwalT2BiqEAHRnniKLrVixgho1anDp0iWcnJzw9fXlu+++0yUHSsEff+A2qBURsaYEaUsQVaclbN0qU5QKkUFr167lwIEDWFlZMXv27FQDI5tUvs+e7w6yaRMUL44kB0JkkaioKA4ePAg8P0GYORM++ABGjgRCQnRLKAM8kdALIbKXJAhZSCnFuHHjePfdd4mPj6dVq1b4+fnRqFEjIiJg94YHulFZXbpQJngv/xbpxZk1V7A+8A/Urm3o8IXI1WJjYxn16NfHzz77LPWqrUrBBx/gPaoetY7NMVCEQuQNu3fvJjk5mZIlS1KqVKln9nt5gbs79O8PrF8PWi3UqgVubtkfrBACkC5GWSY+Ph4fHx9WrlwJwOeff87EiRMxMjLiv/903Zwf3LflbMJ5SpqYwJgx1PviC7lqIEQmmTZtGjdv3qREiRJ8+umngG5Q8kcfwbTqq3HetAlMTWXMgRBZ7GXdiwCGDn1iqYNRa3Qbe/bMpuiEEM8jCUIWiImJoV27duzevRsTExMWLlxI//79dTuVwtlZQ8mSYGpqTHiPKdCztMyaIkQmun//Pj/88AMAP/zwg35BtEGDYN06uLnejYOA5rvvoEoVwwUqRB7wePXkFyUIAEZG6Nb7OXBAlyl065ZN0QkhnkcShEwWHR1N27Zt2bdvH9bW1vzxxx80a9aM4GBwvH8ezbvvYPTbb6xaVQFra8ifv4uhQxbijTN16lSio6OpVq0a3Z440fh2giLc9xQ/PBiIpmFD+OQTA0YpxJsvMDCQ69evY2JiQuPGjVPtCw6GgACoU+fR1YN163Q76teHYsWyP1ghhJ6MQchE0dHRtGnThn379mFjY8OOHTto1qwZ//wDHqUTme61Ujddw6ef4uT0aLYGIUSmCgkJYe7cuQB8++23qQYmlz36KzsfeFHF4gr88sujny2FEFnl33//BaBmzZrY2Nik2rd4MdSrB++++2jDkSO6/8vVAyEMTv51zCRJSUl07tyZ/fv3Y2try86dO6lduzYoReCC7TyMMWNTYhtSGjeD334zdLhCvLGmTJlCXFwctWvXpnXr1gCEhgLh4Y+mSQG+/lo3KlIIkaV2P5qR6OmrBwCJibofypo1e7RhzRo4exZ69crGCIUQzyMJQiZQSvHhhx+yc+dOrKys2LFjBzVr1tTNxDBsGIP+bMVy3sX3oz8w3vEPPFqoSQiRue7fv8+iRYsAmDBhAhqNhs2boWRJmLPCFn78UbdU6+NEQQiRZZRSL00QJkzQzWrao8ejDRoNVKok/0YKkQNIgpAJvv32W5YuXYqRkRFr166lSpWaTP9BS9I7/WHOHDTAu7OqYz5vhsy1LkQWmjdvHnFxcVSvXp1mj36WXLUKYmPhTrAR9OmjW2PE1NTAkQrx5rt+/Tp37tzBzMyMOnXqPLdMvnxgYYHuBzUhRI4hCUIGbd68mfHjxwMwf/582rRpS8+e8OloIz7c1VWXEKxcCcOGGThSId5ssbGxzJmjW9Ng9OjR+rEHK1fCql+T+OorQ0YnRN7z+OpB7dq19TOJgW4Zkvv3nygYG6tbrfDttyEqKpujFEI8jyQIGeDv70/fvn0BGD58OAMGDECjgX79oGBB6L6gKWzbJv0phcgGy5YtIywsDDc3Nzp37qzfbnRwPz3HuJBv/XIDRidE3vM4QWjUqFGq7adOgYMDtGunSxb491/dlEYHD8rsHULkEJIgvKb4+Hi6du1KREQEderUYerUqbodfn60bw+BgdCqkyU0bWrYQIXIA5RS+pmLRowYgYmJCfv2QUoKMHas7uTj0CHDBilEHqKUYs+ePcCz4w8OHND1KLKweDS96bZtuh1t2jzaIIQwNEkQXtNXX32Fn58f9vb2rF27ltWrTYlYuBqqVoXPP8fWRhk6RCHyjH379nHp0iXy5ctH3759uX4dmjSBqmVjiN57QjfmYOxYQ4cpRJ5x5coVQkJCMDc3183o94Rhw8DfXzdIGYBHKy3zkoXUhBDZS0bMvoYDBw4wY8YMAJYsWcKePcV4910or/HkOFbkMzaWX0GEyEY//fQTAH369MHGxoZ9+8DGBopHXiQ/MdD3fV0fZyFEtnjcvahOnTpYWFg8s79kyUd/BATA9eu68XpNmmRjhEKIl5EEIZ1iYmLo168fSin69+/PW2+9xcnt9yhmlEQH7SbydWzxxM8iQoisFhISwh9//AHAoEGDAHjrLQjc6MfDRl11i6GNGWPIEIXIcx53L3p6/MEzHl89qFNHl9ULIXIESRDSaeLEifj7++Ps7Ky7iqDVUn1GH05rT2FXsSisOCirswqRjRYvXkxycjLe3t54enrqt9vOm4QtN6BHL1kUTYhspJRi3759ADRs2DDVvqZNoVw5XY8/R0eke5EQOZScyabD5cuXmT59OgBz585Dq7WF2bNhxw7sLaIxWrtaZmAQIhsppViyZAkAAwcOJDgYLl5Et/rSo6sKfPaZ4QIUIg8KDAwkJCQEU1NT3aKhj1y/rpuwaNGiR2sfANSoAd7ekiAIkcPk+ARh3759tGvXjqJFi6LRaNi0aVOq/Uopxo0bh5OTE5aWljRr1oxr165lehxKKYYOHUpSUhJvvfUWt261o3zZFHaM3qUrMH06lC+f6ccVQrzYoUOHCAgIIH/+/HTp0oVp06BCBRg3zwH274dvv4XKlQ0dphB5ysGDBwGoVq1aqvUPihfXrVM4dSoUKPBo45df6mYYq149+wMVQrxQjk8QYmJi8PT0ZN68ec/dP3XqVGbPns2CBQs4evQo+fLlo2XLlsTHx2dqHP/88w+7du3C3NycGTN+ZNEiCL1nzMVe3+pWZ33U91kIkX1+/fVXALp06UK+fPkIC9PND1Cnrkb3q6SsjiZEtnucINStWzfVdgsL3UymI0YYIiohRHrk+DEIrVu3pnXr1s/dp5Ri1qxZfPXVV3To0AHQnTA4ODiwadMmevbsmSkxaLVavvjiCwA+/vhjSpcuyZEjsGIFvP9+VdCsyJTjCCHSLj4+njVr1gDw7rvvArBsGYwbq3ArKbOICWEohx6tOfJ0gvCMU6d044NkcLIQOU6OTxBe5nE/x2bNmum32draUqtWLQ4fPvzSBCEyMjLVfXNzc8zNzZ9bds2aNZw5cwYbGxvGjBkDDx5godHwwQcFM+eJCJFHJSQkkJCQoL8fFRUFpK19btmyhYiICJydnVPNlFJydFfdMq1ffCFTmwqRQeltow8fPuT8+fNA6gRh/344c0a3erKLC7qV0lq0gPBwXaLwxAQDQoi0ebp9Pt0uMyLHdzF6mZCQEAAcHBxSbXdwcNDvexFnZ2dsbW31t8mTJz+3XFJSEmMfLbA0dOhojh61Q332OZQpA3/+mQnPQoi8a/LkyanaoYeHB5C29vn46sHbb7/NrVtGREUBt2/Dxo0wfz4kJWXnUxHijZTeNnrkyBGUUpQqVSrVv82LF8PQoTBnzqMNFy9CWJiu39GjOoUQ6fN0+3R2ds60unP1FYSMuHnzJjZPXNZ80dWDVatW4e/vT5EiRQgPH0bbtjBEU4k5ahEUlCsIQmTE559/zsiRI/X3b9++jYeHxyvbZ0xMDH///TcA3bt3Z8QI2LEDfu5wmp5KQb164OaWPU9CiDdYetvoi8Yf1KkDQUHQvv2jDXv3/n+HqWmWxS/Em+zp9hkZGZlpSUKuThAcHR0BCA0NxcnJSb89NDSUKlWqvPSxNjY2qb7cnker1TJ16lQARowYQVRUfkyNkmmv3Qxt20L9+hl7AkLkcU93HXp8efRV7XPr1q3ExcVRsmRJKlSoytWrEB0NFU8s0xXo3TsrwxYiz0hvG31RgjBggO6m92idBJ5aJ0EIkXYv6x6fUbm6i5GbmxuOjo74+vrqt0VGRnL06FG8vb0zXP/WrVu5cOECNjY2DBo0iO+6nyFIW4Lm7IKJEzNcvxDi9axbtw6Abt26YW6u4exZOLE2gIpXNoCJCXTrZuAIhch7kpKSOHr0KPCKAcpKSYIgRA6X468gREdHc/36df39wMBA/Pz8KFSoECVKlGD48OFMnDiR0qVL4+bmxtixYylatCgdO3bM8LEfXz0YOHAgtra2MGkSRQmG7t3hFVcohBBZIyYmhq1btwK6BAF0U5tWP7tUV6B1a7CzM1R4QuRZZ86cITY2lgIFClD+iXWBzpzRLRNkZvZow3//6RYzNDXVLZQmhMhxcnyCcOLECRo3bqy//7ivVd++fVm2bBmjR48mJiaGAQMG8PDhQ+rVq8e2bduw0C/T+HpOnz7NgQMHMDU1xc5uODcO3KDE+vW6nV9+maG6hRCv759//iEuLg43NzfKlauGUroEgY0bdQV69DBofELkVY+nN61Tpw5GRroOCjExuhzAzAyuXQMnJ+DYMd0DPD2fWFJZCJGT5PgEoVGjRiilXrhfo9EwYcIEJkyYkKnHnT9/PgBNm3ZlzBgnxpkkc0NrR5EWVWVlViEM6PHVg44dOzJ6tIZdu2Dq5BQ6tG6tmzqxTRsDRyhE3vS88QfXr+tWTbaygkfDBqF2bViwAPLly/4ghRBpkuMTBEOIiIjg999/B6BHj0EkJUHhwiYUGb9Xpk4UwoC0Wi3btm0DoFWrNrz/Pty8CZb5jeGHH3Q3IUS2U0o9N0Hw9NT1JgoOfnSlD6BECfjwQwNEKYRIK0kQnuPXX38lNjaWChUq0LdvPfr1g/h4wKL8qx4qhMhCZ86cISQkhHz58tGwYX0uXNAtR9KkiaEjEyJvu3HjBrdv38bExIQaT40rMDKCYsUMFJgQ4rXk6lmMsoJSigULFgDw0UcfoQG4f1+6SQqRAzxe+6Bp06aYm5tjbQ292z7EZNc2eGI1SSFE9np89aBq1apYWVm9uOB//8GiRXDuXDZFJoR4HZIgPOXkyZNcvHgRCwsLkpN7E3fgpG5U1Tvv6KZmE0IYzD///ANAmyfHGezYoZu5SKZLFMJgnte9aPVq3QDlR0P6dHbt0nUvGjYsmyMUQqSHJAhP+fXXXwGoWLETw4bZ0qxbQUhO5v9TpQghDCE8PJzDhw8DYG/fmo4dYe1aYOdOXYE6dQwWmxB53fMShO3b4cQJ3QrKeo9nMKpZM/uCE0Kkm4xBeEJSUhKrVq0CoEGDdwgLU3QLXazb2bevASMTQuzYsQOtVouHhweHDpVg82awsVF03/coQWje3LABCpFHRUZGcu5Rl6EnE4TvvoPGjZ9aNkgSBCFyBUkQnrBt2zbu37+Pg4MD33/fnB+abie57TQoUkT3LSeEMJgnuxf16wfW1lDP7Q6s+E+34FL9+oYNUIg86siRI2i1Wtzc3HByctJvL1oU3n33iYIJCXD+vO5vL6/sDVIIkS6SIDxhxYoVALz99tuYmJjAujWYkQRdu4KJvFRCGIpWq9UnCK1bt6ZCBahQAViwRVfA2xvy5zdcgELkYc/rXvRcly/ruuwWKADOzlkfmBDitckYhEdiY2P1CzBVrNgblZgEmzbpdvbsabjAhBCcPn2au3fvkj9/furVq/f/Hbt36/7frJlhAhNCPDdBmDtXN0g5IuKJgmfO6P5fubKM6RMih5ME4ZHt27cTGxuLvb0rPj7VaNsgkv+1d+9BVdf5H8dfgBzwxnEJ5aIoXkolES8lYavhygSSZqVGbrNjbuNOdtNlGrN+4213Z53KcSjXdGrTajebLq7ZZWTaRUXLxF3NmlrXlCgvKAqmCI4o8P398T3nq9yEcyC+53iej5kzyPdceM93vi+HN5/P9/PR2bNSz55SS38VAfCzco8epKen669/dWj3bnPTZO3aZb6A6UWALWpqarR7925JVxqEy5elZ5+VZs40d1K2uBuE5OQOrhKAp5g347Jx40ZJUmLifdqzJ0jDbw2Tfv2i+T9dMH0UYCf3/ge33z5Jjz9uLip29IihPps2mU1Cg42ZAHSMr7/+WlVVVXI6nbr55pslSRcuSL/7nfTvf0sjR1714mefNZckjo62p1gArUaDIKm6uloffWTOZf7zn6cpOVm6eLGbdMOTNlcG4MyZMyosLJQkpaRM0v33S6dPS33ig6T4W7jZEbCRe3pRamqqgl1/THM6pRUrmnjxDTcwHRDwEzQIkrZu3aqKigrFxsbqtttuU3Cw1LWr3VUBkK4sbzps2DCNGxevcePYsxDwFa2+QRmAX2HujK5ML5oy5V7zLyB5edKrr0onTthcGQD39KJJkyZZx4KCJP3hD9L69ea9QgBs0bBBuHBB+u9/m2ji9+6VFi6UXHkG4NsCvkGoqanRB67Viv7+92nKypLOrdlgTqBcv97e4oAAV1dXp7y8PEnSrbdmqbzc9URVldkg/Pa30vnz9hUIBLAjR47o2LFjCgkJ0RjXxmf5+eYSxOPHN3jx1q3Sc89Jb7zR8YUC8FjANwg7d+5UeXm5IiJu0IUL4/Xdd4Yidpm/kCgtzdbagEC3b98+nT59Wt27d9eWLbcrOlrKzZX0n/9ItbVSnz6spw7YxD16MHLkSHV1zcs9ckQKD3ftU3I1107LSkrqwAoBeCvgGwT39KLp06fqwIFOeuX/jiio7LTUuTM3PwI2c08vSk9P17FjoaqtlYYNkzldQWL1IsBGTd1/8Nhj0k8/SX/8Y4MXHzhgfm3UOQDwRQHdINTV1WnTpk2SpGnTpmnIEOlX1eZ660pNlRwOG6sD4N7/ICsrS59+KhUXu6Yu7N9vvqDeGooAOlJzNyiHh5tbCFkMw9xFWZKGDOmg6gC0RUCvYlRYWKiSkhJFRERo4sSJ5sEdO8yvjSZQAuhI5eXl1vKmmZmZkqSEBNeTX35pfh0xosPrAmAuP/yVa+OzerubN+X4camyUurUSRo0qAOqA9BWAT2C4J5eFBMzWQsXhun77yXt3Gk+SYMA2CovL0+GYSgpKUm9e/e58sTFi1emKzCCANhi165dMgxDQ4cOVWxsrCRp9mxzHzT3BucW9+jBwIFSaGjHFgrAKwHbIBiGYTUIxcXTlJsrnSs+Ix07Zq6hyNxmwFYff/yxJOmXv5yi+Hhp3jzX0on/+595g3JUlNS7t71FAgFqh2u0fcKECZLMSH74oblKeKMlTpleBPidgJ1i9NVXX+mHH35Qly5d9Prrmfr8c2nEryKlsjLp22+lbt3sLhEIWJcvX7aWNw0Nnazjx6VvvnHtfzBihLmV8pEjrgMAOlrDBiE4WCookP75T8m14ukVjz4qTZkiXbrUwVUC8FbANgibN2+WZG6+NGNGF82Y4XrihhuYXgTYrLCwUGfPnlVUVJSWLx+jjAxzYTFLVJT5AGCLA65pfmmu5cCDgswVxoYNa+LFwcFSv34dVxyANgv4BmH69Ok2VwKgIffoQVZWlrp0CVFWls0FAWgkKSlJUTTqwHUpYO9BKCoqUmhomE6cuEtlZTInTU6bJj3zjLmIMwDbuPc/mDx5cv0nDEPKzpaeeoqcAjZzr/534oT0+9+bU4waOXdOuv9+afFiqa6uYwsE4LWAbRAk6Re/yFROTne9/LKkkhLpH/+QXnjBXMQZgG2KiooUHh6unTsztGqVdOaM64myMundd6WVK8kpYLO77rpLknlzcm6u9PTTTbzou++k996TXnnFnGoEwC9cN2ldvXq1EhISFB4erpSUFO3Zs6fF99x330MaOlSaOVPmjcmSuUZzvcnOAOyQlpalv/wlQk8+af4RUtKV5U0TEsgpYKOuXbtqvOt+vaQkadYs6Te/aeKFRUXmV/Y/APzKdXEPwjvvvKOcnBytXbtWKSkpys3NVUZGhg4ePKhevXo1+Z7IyEi9+GKWQkNdC6F87GoQ2AYe8AkPPPCAMjPNFRL793cddDcILJcI2GrChAlyOBySpLFjzUeT3A3CwIEdUxiAdnFdjCCsXLlSc+bM0ezZs5WYmKi1a9eqS5cuWrduXbPvyc7OlsPhuLJKonsEocklGAB0pO7du2vGjLs0b560Zs1VT7jXUx861Ja6AJjc04taRIMA+CW/H0G4dOmS9u7dq2eeecY6FhwcrPT0dH3xxRfNvi8hYbYqKiqs77t9843ZLTGCAHSY6upqVVdXW9+fP39ekpSZOUM1NTVWRsPCwhQWFsYIAtDBmsvoHXfcoXPnKvS3v4UqM7NG8fEOM6MNMcUI+Nk0zOfVv9e2ld+PIJSVlam2tlbR0dH1jkdHR+vkyZPNvu/pp+PldN4op9Mpp9Opy19+aT5BgwB0mOXLl1sZdDqdSkxMlCS9994kOZ3jrOPLly833+BuEBhBADpEcxlNTExSjx4z9MQTnXXjjbVatmxF0x/ACALws2mYz/j4+Hb77CDDaLQpul8pKSlR7969tWvXLqWmplrHFyxYoIKCAhUWFtZ7fUVFhZxOp7ZsOa6xY83dkoNOn1a3UaMUdOGCVFUlueZVAvh5Nfzrx/Hjx5WYmKgePX7SwoXhmjvX3Hk1LCxMYUFB5kaGlZXSqVNSz552lQ0EjOYyumrVCY0a1U2PPNJZ48bVaOVKNR5BqK6WIiLMHZRPn2ZzQ6CdNTWCEB8fr3PnzikiIqJNn+33U4yioqIUEhKi0tLSesdLS0sVExPT7PvGju125eRFREhnz0qlpTQHQAeypg65uIdHf/wxWBER4ZIaLGVaUWHmlF80gA7RXEanT++imJhu2rtXungxRE3NLlJYmPlHt6NHzeYeQLtqmM/25PdTjBwOh0aPHq38/HzrWF1dnfLz8+uNKLi5O62rOy5J5lJG12goAll1dbWWLl3a+JzhmjhvnnPPb272nLlzaq0uAIlrzRucM++4MxoSYp630FCpe/drvKFTJ3MZsgDOLNeadzhvnmv2d1wv+H2DIEk5OTl69dVX9cYbb+jAgQOaO3euqqqqNHv27Eavbc+TFyiqq6u1bNkyzpmHOG+eq6yslEQ+PcW15jnOmXfIqOe41rzDefMcDUID2dnZWrFihRYvXqwRI0Zo//79ysvLa3TjcrPmzZPuvVe6xqpHAGy2YoWZ048/trsSAK2Rmys9+KD06ad2VwLAQ9dFgyBJjz/+uH788UdVV1ersLBQKSkprX/zv/4lffCB5Bo6BeCDduwwc3rsmN2VAGiN7dulDRuurGQEwG/4/U3KnnIv2nT+/HnzZqu6uiv/efXqZd4EiXrcN6W15/q6gYDz5jn3/GYrn1c7dMj8Gh1NThvgWvMc58w718xoQ8XF5teoqIDOLNeadzhvnnPnsz0WKPX7ZU499f3332sg6zEDAADgOlRUVKQBAwa06TMCrkGoq6tTSUmJunfvrqAAXlUB8EW1tbU6fPiwBg0apJCQELvLAdAAGQV8l2EYOn/+vOLi4hQc3La7CAKuQQAAAADQvOvmJmUAAAAAbUeDAAAAAMBCgwAAAADAElANwurVq5WQkKDw8HClpKRoz549dpfk05YuXaqgoKB6jyFDhthdlk/ZsWOHpkyZori4OAUFBemDDz6o97xhGFq8eLFiY2PVuXNnpaen65B7uU40QkZbj3y2DhltP+TTM2S0ZeTTOy2dt4ceeqjRtZeZmenRzwiYBuGdd95RTk6OlixZon379ik5OVkZGRk6deqU3aX5tJtvvlknTpywHp999pndJfmUqqoqJScna/Xq1U0+//zzz+ull17S2rVrVVhYqK5duyojI0MXL17s4Ep9Hxn1HPlsGRltH+TTO2T02sind1o6b5KUmZlZ79p7++23PfshRoAYM2aM8dhjj1nf19bWGnFxccby5cttrMq3LVmyxEhOTra7DL8hydi0aZP1fV1dnRETE2O88MIL1rGzZ88aYWFhxttvv21Dhb6NjHqGfHqOjHqPfHqOjHqGfHqn4XkzDMOYNWuWMXXq1DZ9bkCMIFy6dEl79+5Venq6dSw4OFjp6en64osvbKzM9x06dEhxcXEaMGCAHnzwQR05csTukvxGcXGxTp48We+6czqdSklJ4bprgIx6h3y2DRltHfLpPTLqPfLZNtu3b1evXr00ePBgzZ07V+Xl5R69PyAahLKyMtXW1io6Orre8ejoaJ08edKmqnxfSkqKXn/9deXl5WnNmjUqLi7WuHHjrK28cW3ua4vrrmVk1HPks+3IaOuQT++Q0bYhn97LzMzUm2++qfz8fD333HMqKCjQpEmTVFtb2+rP6PQz1gc/N2nSJOvfw4cPV0pKivr166d3331XDz/8sI2VASCfgG8jo7DLAw88YP07KSlJw4cP18CBA7V9+3ZNnDixVZ8RECMIUVFRCgkJUWlpab3jpaWliomJsakq/9OjRw/ddNNNOnz4sN2l+AX3tcV11zIy2nbk03NktHXIZ/sgo54hn+1nwIABioqK8ujaC4gGweFwaPTo0crPz7eO1dXVKT8/X6mpqTZW5l8qKytVVFSk2NhYu0vxC/3791dMTEy9666iokKFhYVcdw2Q0bYjn54jo61DPtsHGfUM+Ww/x44dU3l5uUfXXsBMMcrJydGsWbN0yy23aMyYMcrNzVVVVZVmz55td2k+66mnntKUKVPUr18/lZSUaMmSJQoJCdHMmTPtLs1nVFZW1uvIi4uLtX//fkVGRqpv376aP3++/vSnP+nGG29U//79tWjRIsXFxemee+6xr2gfRUY9Qz5bh4y2D/LpOTLaMvLpnWudt8jISC1btkzTpk1TTEyMioqKtGDBAg0aNEgZGRmt/yFtWgPJz6xatcro27ev4XA4jDFjxhi7d++2uySflp2dbcTGxhoOh8Po3bu3kZ2dbRw+fNjusnzKtm3bDEmNHrNmzTIMw1ymbdGiRUZ0dLQRFhZmTJw40Th48KC9RfswMtp65LN1yGj7IZ+eIaMtI5/eudZ5u3DhgnHnnXcaPXv2NEJDQ41+/foZc+bMMU6ePOnRzwgyDMNoYyMDAAAA4DoREPcgAAAAAGgdGgQAAAAAFhoEAAAAABYaBAAAAAAWGgQAAAAAFhoEAAAAABYaBAAAAAAWGgQAAAAAFhoE+LS0tDTNnz/f7jIANIF8Ar6NjMJbNAjwe8uWLVOfPn0UFBR0zcf27dvtLhUIOOQT8G1kFE3pZHcBQFtt3rxZK1eu1Pjx461j8+bNU0VFhdavX28di4yMtKM8IKCRT8C3kVE0hREE+JVPPvlETqdTb731liTp6NGj+vbbb5WZmamYmBjr0blzZ4WFhdU75nA4bK4euL6RT8C3kVG0FiMI8BsbNmzQI488og0bNmjy5MmSpA8//FBpaWmKiIiwuTogsJFPwLeRUXiCEQT4hdWrV+vRRx/VRx99ZP3HJplDo3fffbeNlQEgn4BvI6PwFCMI8Hnvv/++Tp06pc8//1y33nqrdbyiokIFBQV67bXXbKwOCGzkE/BtZBTeYAQBPm/kyJHq2bOn1q1bJ8MwrONbtmxRYmKi4uPjbawOCGzkE/BtZBTeoEGAzxs4cKC2bdumzZs364knnrCOb968WVOnTrWxMgDkE/BtZBTeoEGAX7jpppu0bds2bdy4UfPnz1dNTY22bNnC3EnAB5BPwLeRUXiKexDgNwYPHqytW7cqLS1NBQUF6tatm0aNGmV3WQBEPgFfR0bhiSDj6glpgJ948sknVVNTo5dfftnuUgA0QD4B30ZG0RJGEOCXhg0bptTUVLvLANAE8gn4NjKKljCCAAAAAMDCTcoAAAAALDQIAAAAACw0CAAAAAAsNAgAAAAALDQIAAAAACw0CAAAAAAsNAgAAAAALDQIAAAAACw0CAAAAAAs/w/V6mCw8NDtHQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 900x300 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams['ytick.right'] = True \n",
    "plt.rcParams['ytick.labelright'] = False \n",
    "plt.rcParams['ytick.left'] = plt.rcParams['ytick.labelleft'] = True\n",
    "plt.rcParams[\"xtick.top\"] =  True\n",
    "plt.rcParams['xtick.direction']='in'\n",
    "plt.rcParams['ytick.direction']='in'\n",
    "\n",
    "fig=plt.figure()\n",
    "cc=1.5\n",
    "fig.set_size_inches(6*cc,2*cc)\n",
    "gs = fig.add_gridspec(1, 3, hspace=0, wspace=0)\n",
    "axs = gs.subplots(sharey=True)\n",
    "styles=[\"r--\",\"b:\",\"k\"]\n",
    "labels = [\"strict LO\",\"subtracted\",\"tuned\"]\n",
    "for i, ax in enumerate(axs):\n",
    "    if i>0:\n",
    "        labels=3*[None]\n",
    "    for j in range(0,3):\n",
    "        ax.plot(k,GWtestrate[i][j+1],styles[j],label=labels[j])\n",
    "    ax.set_xlim(0,15)\n",
    "    ax.set_ylim(-2,60)\n",
    "    ax.set_xlabel(\"k/T\")\n",
    "    if i==0:\n",
    "        ax.set_ylabel(r\"$\\Gamma_\\mathrm{GW}\\; m_\\mathrm{Pl}^2/T^3$\")\n",
    "    ax.set_title(\"$T=10^{\"+str(tempslog[i])+\"}$ GeV\")\n",
    "    ax.set_xticks(ax.get_xticks().tolist())\n",
    " \n",
    "    tlabels = ax.get_xticklabels()\n",
    "    # remove the last  label\n",
    "    if i<2:\n",
    "        tlabels[-1] =  \"\"\n",
    "    # set these new labels\n",
    "    ax.set_xticklabels(tlabels)\n",
    "\n",
    "plt.figlegend(loc='upper center',ncols=6,bbox_to_anchor=(0.5,1.1),bbox_transform=plt.gcf().transFigure)\n",
    "plt.savefig(\"../../../gw_SM_rate.pdf\",bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4c08ebc",
   "metadata": {},
   "source": [
    "## SMASH, reproducing the results of [2011.04731](https://arxiv.org/abs/2011.04731)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06905656",
   "metadata": {},
   "source": [
    "The configuration file is altogether similar to the one above, but points to `analytical/models/SMASH_GW.fr`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecb7bb43",
   "metadata": {},
   "outputs": [],
   "source": [
    "smashGW=analytical_pipeline(\"../../MyModels/smashgrav/smashgrav.cfg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "64c31bc1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "-  $$ \\left\\vert \\mathcal{M}(1;-1,-1)\\right\\vert^2=3 \\kappa^{2} s \\left(2 \\vert h_t\\vert^2 + \\vert y_Q\\vert^{2}\\right) + \\frac{2 \\kappa^{2} \\left(t^{2} + u^{2}\\right) \\left(16 g_{1}^{2} + 27 g_{2}^{2} + 84 g_{3}^{2}\\right)}{3 s}$$\n",
       "-  $$ \\left\\vert \\mathcal{M}(-1;-1,1)\\right\\vert^2=- 3 \\kappa^{2} t \\left(2 \\vert h_t\\vert^2 + \\vert y_Q\\vert^{2}\\right) - \\frac{2 \\kappa^{2} \\left(s^{2} + u^{2}\\right) \\left(16 g_{1}^{2} + 27 g_{2}^{2} + 84 g_{3}^{2}\\right)}{3 t}$$\n",
       "-  $$ \\left\\vert \\mathcal{M}(-1;1,-1)\\right\\vert^2=- 3 \\kappa^{2} u \\left(2 \\vert h_t\\vert^2 + \\vert y_Q\\vert^{2}\\right) - \\frac{2 \\kappa^{2} \\left(s^{2} + t^{2}\\right) \\left(16 g_{1}^{2} + 27 g_{2}^{2} + 84 g_{3}^{2}\\right)}{3 u}$$\n",
       "-  $$ \\left\\vert \\mathcal{M}(1;1,1)\\right\\vert^2=\\frac{\\kappa^{2} \\left(g_{1}^{2} + 15 g_{2}^{2} + 48 g_{3}^{2}\\right) \\left(s^{2} + t^{2} + u^{2}\\right)^{2}}{4 s t u}$$\n"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, Markdown\n",
    "out=\"\"\n",
    "s = sympy.Symbol(\"s\")\n",
    "t = sympy.Symbol(\"t\")\n",
    "u = sympy.Symbol(\"u\")\n",
    "yq=sympy.S('yq')\n",
    "for key, item, in smashGW[0].items():\n",
    "    sitem=sympy.simplify(sympy.sympify(item).subs(2*t*u,(s*s-t*t-u*u)))\n",
    "    sitemsubst= sitem\n",
    "    sitemsubstcollect=0\n",
    "    for exprtemp in sympy.factor(sitemsubst.expand().as_independent(ht,yq)).subs(t+u,-s)\\\n",
    "        .subs(2*t*t+2*t*u,2*t*t+(s*s-t*t-u*u)).subs(t*t+2*t*u,t*t+(s*s-t*t-u*u)).subs(t*t+t*u,t*t+(s*s-t*t-u*u)/2):\n",
    "        sitemsubstcollect += exprtemp.simplify()\n",
    "        # again, beautify output\n",
    "    out+=f\"-  $$ \\\\left\\\\vert \\\\mathcal{{M}}({key[0]};{key[1]},{key[2]})\\\\right\\\\vert^2={sympy.latex(sitemsubstcollect).replace('ht^{2}',r'\\vert h_t\\vert^2').replace('yq',r'\\vert y_Q\\vert')}$$\\n\"\n",
    "    # display(sympy.pprint(sitem))\n",
    "display(Markdown(out))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e308484",
   "metadata": {},
   "source": [
    "If we subtract off the SM contribution we have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "5320cf87",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "-  $$\\Delta \\left\\vert \\mathcal{M}(1;-1,-1)\\right\\vert^2=3 \\kappa^{2} s \\,\\vert y_Q\\vert^{2} + \\frac{2 \\kappa^{2} \\left(g_{1}^{2} + 12 g_{3}^{2}\\right) \\left(t^{2} + u^{2}\\right)}{3 s}$$\n",
       "-  $$\\Delta \\left\\vert \\mathcal{M}(-1;-1,1)\\right\\vert^2=- 3 \\kappa^{2} t \\,\\vert y_Q\\vert^{2} - \\frac{2 \\kappa^{2} \\left(g_{1}^{2} + 12 g_{3}^{2}\\right) \\left(s^{2} + u^{2}\\right)}{3 t}$$\n",
       "-  $$\\Delta \\left\\vert \\mathcal{M}(-1;1,-1)\\right\\vert^2=- 3 \\kappa^{2} u \\,\\vert y_Q\\vert^{2} - \\frac{2 \\kappa^{2} \\left(g_{1}^{2} + 12 g_{3}^{2}\\right) \\left(s^{2} + t^{2}\\right)}{3 u}$$\n",
       "-  $$\\Delta \\left\\vert \\mathcal{M}(1;1,1)\\right\\vert^2=0$$\n"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "out=\"\"\n",
    "s = sympy.Symbol(\"s\")\n",
    "t = sympy.Symbol(\"t\")\n",
    "u = sympy.Symbol(\"u\")\n",
    "for key, item, in smashGW[0].items():\n",
    "    sitem=sympy.sympify(item).subs(2*t*u,(s*s-t*t-u*u)).simplify()\n",
    "    sitemSM=sympy.sympify(GWdict[0][key]).subs(2*t*u,(s*s-t*t-u*u)).simplify()\n",
    "    diff=sympy.simplify((sitem-sitemSM))\n",
    "    diffcollect=0\n",
    "    for exprtemp in sympy.factor(diff.expand().as_independent(ht,yq)).subs(t+u,-s)\\\n",
    "        .subs(2*t*t+2*t*u,2*t*t+(s*s-t*t-u*u)).subs(t*t+2*t*u,t*t+(s*s-t*t-u*u)).subs(t*t+t*u,t*t+(s*s-t*t-u*u)/2):\n",
    "        diffcollect += exprtemp.simplify()\n",
    "    out+=f\"-  $$\\\\Delta \\\\left\\\\vert \\\\mathcal{{M}}({key[0]};{key[1]},{key[2]})\\\\right\\\\vert^2={sympy.latex(diffcollect).replace('yq',r'\\,\\vert y_Q\\vert')}$$\\n\"\n",
    "    # display(sympy.pprint(sitem))\n",
    "display(Markdown(out))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71e10eb8",
   "metadata": {},
   "source": [
    "The contribution from the BSM Yukawa $y_Q$ agrees with that of [2011.04731](https://arxiv.org/abs/2011.04731). Its contribution through U(1) and SU(3) interactions also agrees"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "142114ce",
   "metadata": {},
   "source": [
    "The thermal masses agree with [2011.04731](https://arxiv.org/abs/2011.04731)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6ecb1876",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('(35*g1^2)/18', '(11*g2^2)/6', '(13*g3^2)/6')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smashGW[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1e3ed04",
   "metadata": {},
   "source": [
    "and the leading-log term has the usual form in terms of the masses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1600ca31",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{T \\kappa^{2} \\left(m_{D1}^2 \\log{\\left(\\frac{4 k^{2}}{m_{D1}^2} \\right)} + 3 m_{D2}^2 \\log{\\left(\\frac{4 k^{2}}{m_{D2}^2} \\right)} + 8 m_{D3}^2 \\log{\\left(\\frac{4 k^{2}}{m_{D3}^2} \\right)}\\right)}{32 \\pi}$"
      ],
      "text/plain": [
       "T*kappa**2*(m_{D1}^2*log(4*k**2/m_{D1}^2) + 3*m_{D2}^2*log(4*k**2/m_{D2}^2) + 8*m_{D3}^2*log(4*k**2/m_{D3}^2))/(32*pi)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smashGWrate=NumRate(*smashGW,2)\n",
    "smashGWrate.get_leadlog().subs(g1*g1,18*md1/(35*T**2)).subs(g2*g2,6*md2/(11*T**2)).subs(g3*g3,6*md3/(13*T**2)).simplify()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv (3.12.11)",
   "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.12.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}