{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dynamics of multiple spin enembles: two driven-dissipative ensembles \n",
    "\n",
    "Notebook author: Nathan Shammah (nathan.shammah at gmail.com)\n",
    "\n",
    "We use the Permutational Invariant Quantum Solver (PIQS) library, imported in QuTiP as $\\texttt{qutip.piqs}$ to study the driven-dissipative open quantum dynamics of multiple two-level-system (TLS), or spin, ensembles.\n",
    "\n",
    "We consider a system of two TLS ensembles with populations $N_1$ and $N_2$ with identical frequency $\\omega_{0}$ with collective pumping and collective emission at identical rates, $\\gamma_\\text{CE}=(1+\\bar{n})\\gamma_0$ and $\\gamma_\\text{CP}=\\bar{n}\\gamma_0$, respectively, with $\\bar{n}=\\frac{1}{e^{\\hbar\\omega_0/k_\\mathrm{B}T}-1}$ and \n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\dot{\\rho} &=& \n",
    "-i\\lbrack \\omega_{0}\\left(J_z^{(1)}+J_z^{(2)}\\right),\\rho \\rbrack\n",
    "+\\frac{\\gamma_\\text {CE}}{2}\\mathcal{L}_{J_{-}^{(1)}+ J_{-}^{(2)}}[\\rho]\n",
    "+\\frac{\\gamma_\\text {CP}}{2}\\mathcal{L}_{J_{+}^{(1)}+J_{+}^{(2)}}[\\rho]\n",
    "\\end{eqnarray}\n",
    "\n",
    "Ref. [2] has shown that for $N_1<N_2$, if the system is initialized in the state $|{\\psi_0}\\rangle=|{\\downarrow\\cdots\\downarrow}\\rangle_1\\otimes|{\\uparrow\\cdots\\uparrow}\\rangle_2$, the system relaxes to a steady state for which the first subsystem is excited, i.e. $\\langle J_z^{(1)}(\\infty)\\rangle>0$ and for some parameters  $\\frac{\\langle J_z^{(1)}(\\infty)\\rangle}{(N_1/2)}\\rightarrow 0.5$, also in the limit of zero temperature, $T\\rightarrow 0$.  \n",
    "\n",
    "Notice that $\\mathcal{L}_{J_{-}^{(1)}+ J_{-}^{(2)}}[\\rho]\\neq \\mathcal{L}_{J_{-}^{(1)}}[\\rho]+\\mathcal{L}_{ J_{-}^{(2)}}[\\rho]$, which is a case treated in Ref. [3] two obtain syncronized ensembles of atoms. \n",
    "\n",
    "Here we explore what happens when to the master equation of Eq. (1) one adds also collective and local terms relative to single ensembles, \n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\dot{\\rho} &=& \n",
    "-i\\lbrack \\omega_{0}\\left(J_z^{(1)}+J_z^{(2)}\\right),\\rho \\rbrack\n",
    "+\\frac{\\gamma_\\text{CE}}{2}\\mathcal{L}_{J_{-}^{(1)}+ J_{-}^{(2)}}[\\rho]\n",
    "+\\frac{\\gamma_\\text{CP}}{2}\\mathcal{L}_{J_{+}^{(1)}+J_{+}^{(2)}}[\\rho]\\\\\n",
    "&& +\\frac{\\gamma_\\text{CEi}}{2}\\mathcal{L}_{J_{-}^{(1)}}[\\rho]\n",
    "+\\frac{\\gamma_\\text{CEi}}{2}\\mathcal{L}_{J_{-}^{(2)}}[\\rho]\n",
    "+\\sum_{n}^{N_1}\\frac{\\gamma_\\text{E}}{2}\\mathcal{L}_{J_{-,n}^{(1)}}[\\rho]+\\frac{\\gamma_\\text{D}}{2}\\mathcal{L}_{J_{z,n}^{(1)}}[\\rho]+\\sum_{n}^{N_2}\\frac{\\gamma_\\text{E}}{2}\\mathcal{L}_{J_{-,n}^{(2)}}[\\rho]+\\frac{\\gamma_\\text{D}}{2}\\mathcal{L}_{J_{z,n}^{(2)}}[\\rho]\n",
    "\\end{eqnarray}\n",
    "\n",
    "where $\\gamma_\\text {CEi}$ is the rate of superradiant decay for the individual ensembles of TLSs, $\\gamma_\\text{E}$ and $\\gamma_\\text{D}$ are the rates of local emission and dephasing.\n",
    "\n",
    "Firstly, we will show how the collective dynamics of Eq. (1) can be investigated in a simple way using QuTiP's [4] $\\texttt{jmat}$ function, which defines collective spins for maximally symmetric states in a Hilbert space of dimension $N_i+1$.\n",
    "\n",
    "Secondly, we will exploit the permutational invariance of the local processes in Eq. (2) to investigate the exact dynamics using the Dicke basis, $\\rho = \\sum_{j,m,m'}p_{jmm'}|j,m\\rangle\\langle j,m'|$ [1], where $p_{jmm'}$ is a probability density. We will do so numerically using the PIQS library [1]. \n",
    "\n",
    "In the following we might use in plots thefollowing equivalent notation $\\gamma_\\text {CE}=\\gamma_\\Downarrow$ (gCE),\n",
    "$\\gamma_\\text {CP}=\\gamma_\\Uparrow$ (gCP), $\\gamma_\\text {E}=\\gamma_\\downarrow$ (gE), $\\gamma_\\text {P}=\\gamma_\\uparrow$ (gP), and \n",
    "$\\gamma_\\text {D}=\\gamma_\\phi$ (gD)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qutip import *\n",
    "from qutip.piqs import *\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1) Collective processes only (QuTiP $\\texttt{jmat}$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### System properties - QuTiP jmat()\n",
    "QuTiP's jmat() functions span the symmetric (N+1)-dimensional Hilbert space. They can be used to efficiently investigate the collective dynamics only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n0 = 1309202.45774\n",
      "gCE = 1309203.45774\n",
      "gCP = 1309202.45774\n"
     ]
    }
   ],
   "source": [
    "# Number of TLSs in the two ensembles\n",
    "N1 = 1\n",
    "N2 = 4\n",
    "N = N1 + N2\n",
    "\n",
    "# TLSs bare frequency\n",
    "w0 = 1\n",
    "\n",
    "# Bose-Einstein distribution determines the occupation number\n",
    "giga = 10**(6)\n",
    "frequency_hertz = w0*10*giga \n",
    "temperature_kelvin = 10**(2)\n",
    "x = (frequency_hertz / temperature_kelvin) * (constants.hbar / constants.Boltzmann)\n",
    "n0 = 1/(np.exp(x)-1)\n",
    "print(\"n0 =\",n0)\n",
    "\n",
    "# set collective pumping and collective emission rates (coupled ensembles) \n",
    "g0 = 1\n",
    "gCE = g0 * (1 + n0)\n",
    "gCP = g0 * n0\n",
    "print(\"gCE =\", gCE)\n",
    "print(\"gCP =\", gCP)\n",
    "\n",
    "# define identity operators and norms in the tensor space \n",
    "dim1_mat = N1 + 1\n",
    "dim2_mat = N2 + 1\n",
    "id1_mat = qeye(dim1_mat)\n",
    "id2_mat = qeye(dim2_mat)\n",
    "norm2 = id2_mat.tr()\n",
    "norm1 = id1_mat.tr()\n",
    "\n",
    "# build collective spin operators for N1 and N2  \n",
    "jx1_mat = jmat(N1/2,\"x\")\n",
    "jx2_mat = jmat(N2/2,\"x\")\n",
    "jy1_mat = jmat(N1/2,\"y\")\n",
    "jy2_mat = jmat(N2/2,\"y\")\n",
    "jz1_mat = jmat(N1/2,\"z\")\n",
    "jz2_mat = jmat(N2/2,\"z\")\n",
    "jm1_mat = jmat(N1/2,\"-\")\n",
    "jm2_mat = jmat(N2/2,\"-\")\n",
    "\n",
    "# place collective spin operators in tensor space (N1 + N2)  \n",
    "jz1_tot = tensor(jz1_mat, id2_mat)\n",
    "jz2_tot = tensor(id1_mat, jz2_mat)\n",
    "jx12_mat = tensor(jx1_mat, id2_mat) + tensor(id1_mat, jx2_mat)\n",
    "jy12_mat = tensor(jy1_mat, id2_mat) + tensor(id1_mat, jy2_mat)\n",
    "jz12_mat = tensor(jz1_mat, id2_mat) + tensor(id1_mat, jz2_mat)\n",
    "jm12_mat = tensor(jm1_mat, id2_mat) + tensor(id1_mat, jm2_mat) \n",
    "jp12_mat = jm12_mat.dag()\n",
    "\n",
    "# define Hamiltonian  \n",
    "h1_mat = w0 * jz1_mat\n",
    "h2_mat = w0 * jz2_mat\n",
    "htot = tensor(h1_mat, id2_mat) + tensor(id1_mat, h2_mat) \n",
    "\n",
    "# build Liouvillian using QuTiP\n",
    "collapse_operators = [np.sqrt(gCE)*jm12_mat, np.sqrt(gCP)*jp12_mat]\n",
    "L_collective = liouvillian(htot, collapse_operators)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "#Check the algebra of the spin operators in the tensor space\n",
    "\n",
    "print(jp12_mat*jm12_mat - jm12_mat*jp12_mat == 2*jz12_mat)\n",
    "print(jx12_mat*jy12_mat - jy12_mat*jx12_mat == 1j*jz12_mat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set superradiant delay time for the excited ensemble (N2)\n",
    "td0 = np.log(N2)/(N2*gCE)\n",
    "tmax = 30 * td0\n",
    "nt = 1001\n",
    "t = np.linspace(0, tmax, nt)\n",
    "\n",
    "#set initial tensor state for spins (Use QuTiP's jmat() basis)\n",
    "excited1 = np.zeros(jz1_mat.shape)\n",
    "excited2 = np.zeros(jz2_mat.shape)\n",
    "ground1 = np.zeros(jz1_mat.shape)\n",
    "ground2 = np.zeros(jz2_mat.shape)\n",
    "excited1[0,0] = 1\n",
    "excited2[0,0] = 1\n",
    "ground1[-1,-1] = 1\n",
    "ground2[-1,-1] = 1\n",
    "\n",
    "excited1 = Qobj(excited1)\n",
    "excited2 = Qobj(excited2)\n",
    "ground1 = Qobj(ground1)\n",
    "ground2 = Qobj(ground2)\n",
    "\n",
    "sdp =  tensor(excited1, excited2)\n",
    "sdap =  tensor(ground1, excited2)\n",
    "ground12 =  tensor(ground1, ground2)\n",
    "\n",
    "rho0 = sdap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#solve using qutip (using QuTiP's jmat() basis)\n",
    "result = mesolve(L_collective, rho0, t, [], \n",
    "                 e_ops = [jz12_mat, jz1_tot, jz2_tot], \n",
    "                 options = Options(store_states=True))\n",
    "rhot = result.states\n",
    "jzt = result.expect[0]\n",
    "jz1t = result.expect[1]\n",
    "jz2t = result.expect[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvUAAAGICAYAAADBFh8PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3Xl4m1eZ///3keQ9tmUlabYmTZykW5oujrvQvY1TCrRl2sbN0E6nwJS4AwwzMJAAM3wpP2boxDAzDLvdgQLD0sRugUIpYHeBbrRx3CVdaJMoafbVlu3YsS1L5/fHI8u7LdtSZNmf13X5snT0nKNbPk586+h+zmOstYiIiIiISOpyJTsAEREREREZHyX1IiIiIiIpTkm9iIiIiEiKU1IvIiIiIpLilNSLiIiIiKQ4JfUiIiIiIilOSb2IiIiISIpTUi8iIiIikuKU1IuIiIiIpDhPsgNIRTNmzLALFy5MdhgiIiIiMslt2bLlqLV25kjHKakfg4ULF1JXV5fsMERERERkkjPGvBPLcSq/ERERERFJcUrqRURERERSnJJ6EREREZEUp6ReRERERCTFKakXEREREUlxSupFRERERFLcpNnS0hizGghYa2tjPL4IKAb8QCHgj7WviIiIiMhEMimSemNMCXA/UBrj8YXABmvtql5tVcYYv7XWn6AwRUREREQSIqWT+khyvh7YAjSMomsZUNGvrQLYQIxvDERERETioaOjg4aGBlpaWgiFQskORxLI7XaTm5uLz+cjIyMjrmOndFIfWVUvAzDGrB9F19UMTOrrgJo4hSYiIiIyoo6ODnbv3k1BQQELFy4kLS0NY0yyw5IEsNYSDAZpbm5m9+7dLFiwIK6J/ZQ7UdYY48Wpoe+zsm+tDUQeL0xGXCIiIjL1NDQ0UFBQwIwZM0hPT1dCP4kZY0hPT2fGjBkUFBTQ0DCaIpORTbmkHvBBTxI/CCX1IiIiclK0tLSQl5eX7DDkJMvLy6OlpSWuY07FpN6b7ADGwlrLvuP72HJoS7JDERERkTgJhUKkpaUlOww5ydLS0uJ+/kRK19SfTMaYtcBagAULFpz05993fB/vefg9+DJ9PHXbU/p4TkREZJLQ3/SpJxFzPhVX6oFobX3MrLWV1tpia23xzJkzExXWkOZNm0d+Rj4N7Q0caD1w0p9fRERERCauqZjUd+9D7+vd2CvJn5D71BtjWDZ9GQCvH3s9ydGIiIiIxJff76e8vDzZYcSkrKws2SEMMOWS+sgJsn4G1tb7cK5IOyGTeiCa1L929LUkRyIiIiISXxs2bKCoqGjE48rKyigtLcUYgzGG0tJS/P6Tm755vV6qq6tP6nOOZMol9RG1QHG/tqJI+4R1zoxzAHj9qFbqRUREZHKpra2lpKRkxOMqKiqoqqrC6/WyevVqqqqqKCw8uZsXlpWVUVHR/5JHyTWZknofg+xsY4zxGmN2RE507baegVeOLYu0T1jdK/VvHHuDsA0nORoRERGR+Kiurmb16tUxH+/3+wkEAqxZsyaBUQ2tsLCQhoYGAoGhdkg/+VI6qY8k7BuMMVU4Cf0GY0yFMab/b0Wf+vlICc76SN/Vxph1wIaJXHoDMCtnFjOzZtISbGF38+5khyMiIiISFxUVFaOqU6+tdYorYlnZT5SysjIqKyuT9vz9pfSWlt3JeQzHFAzSXg/UJyi0hFk2YxlP7XmK1469xsL8hckOR0RERGRcuuvhR1NCU1NTQ2FhIV5v8i4/tHbtWlasWMG6deuSFkNvKb1SPxWdM1119SIiIjJ5jHaVHmKvv0+0wsJC6usnxhqxkvoUs2yGtrUUERGRyWOs9fSrVq1KYFSxmUgnzCqpTzHdJ8u+eexNusJdSY5GREREZOzGsuI+Eerpu5WUlFBXV5fsMIAUr6mfigoyC5g3bR77ju9jR2AHZ/jOSHZIIiIiEmcLP/toskMYlV3/8b4x9auoqOBzn/vcqPqMt57e7/ePqn5/pOPXrFkz6k8bEkFJfQr62PkfI8Odwdxpc5MdioiIiMi4jDY5r62t5bbbbhvTc9XX148qqY/1+IaGhjHFE08qv0lBNy6+kesWXkduem6yQxEREREZs+Fq0gOBAIsXL+5z5dbx1NMHAgHuu+++uB+/ceNG1q5dO+JxiaakXkRERESSoqSkJFoj39+mTZtoaGjos0peUVERvZJsf+Xl5SxevDia+Hdvldmtrq6OQCBATU0NlZWV0cfr6+uprKyktra2T/tQx/dWX19PcXHxmF9/PKn8JkU9suMRXjzwIv9c/M8UZA7Yhl9ERERS2Fhr1FNRSUnJoDXpxcXF3H///RQVFQFOAl1eXk5NTc2AMQKBACUlJaxbt47a2lp8Pt+AkpmSkhKqqqpYtWpV9Ln8fj/r16/vM+aKFSvYsmXLoMf3N5btOBNFK/Up6qG3H+JXO37Fq0deTXYoIiIiImM2VAlOUVERDQ0NlJeXU15ezn333ceOHTsG3fXG6/VSVFREIBCguLg45jr96urq6JuGboWFhUN+etBfXV3dgP7JopX6FLXmjDWsOm0VSwuWJjsUERERkTHrXlEPBAIDkvHR1KoHAgHASfAHG6u/oS4aFQgEomP1P753Al9ZWcmaNWtiji/RtFKfot5b+F7+5uy/0Q44IiIikvJKS0uprKwcc3+/3x/dLz4QCLBp06Zoe29erze6U01DQwOrV68ekNw3NDREPw3of3xvFRUVE+IE2W5K6kVEREQkqdauXcvGjRvH1DcQCODz+SguLmblypWUlpZGt7xcv34969evjx5bVlbGli1bqK6uprCwkMLCQtavX095eTnV1dWUl5dTVVUVXeXvf3w3v9+Pz+cb8175iWCstcmOIeUUFxfbiXD1sM0HN/P8/ue5afFNLMxfmOxwREREZJTefPNNzjrrrGSHMSGUlZVRVlYW9xr1sVy1diTr16/nwgsvHNcFp2Kde2PMFmvtiFvsaKU+hT287WHu33o/fz7w52SHIiIiIjIuZWVlMZ+gOhqD1cePl9/vT/oVZPtTUp/Czpt5HgAvHX4pyZGIiIiIjE9RURHr1q2L65iJWKUHqKqqivuY46Xdb1JY0Szn46kth7ZgrcUYk+SIRERERCaORCT0E5VW6lPYEu8S8jPyOdR2iP2t+5MdjoiIiIgkiZL6FOYyLopO6VmtFxEREZGpSUl9ilsxawUAdQeTvxuPiIiIiCSHkvoUVzzL2eFIK/UiIiIiU5eS+hR3hu8Msj3Z7G7ZzeG2w8kOR0RERESSQEl9ivO4PFxwygUA1B+qH+FoEREREZmMlNRPAsWznRKcukOqqxcRERGZipTUTwLdJ8uqrl5ERERkalJSPwksm76MDHcGDe0NtAXbkh2OiIiIiJxkSuongXR3Or94/y948rYnyU7LTnY4IiIiIqPm9/spLy9PagxlZWVJff7xUFI/SczPnY/LaDpFREQkNW3YsIGioqIRjysrK6O0tBRjDMYYSktL8fv9cYnB6/VSXV0dl7FONk+yA5D4ag22kpOWk+wwREREREaltraWioqKEY/rPqagoICSkhKqqqriFkNZWRllZWWsXr06bmOeLFranUQ+WvtRLvv5Zexp2ZPsUERERERiVl1dPapE2u/3EwgEWLNmTVzjKCwspKGhgUAgENdxTwYl9ZNIhjsDgLcb305yJCIiIiKxq6ioGFU9e21tLQAlJSVxj6WsrIzKysq4j5toSuonkXUXruPZDzzLygUrkx2KiIiISEy66+ELCwtj7lNTU0NhYSFerzfu8axdu5aNGzfGfdxEU1I/icyZNkf19CIiIpJSRrtKD85KfSJW6bsVFhZSX1+fsPETQUn9JBS2YTpCHckOQ0RERGREY62nX7VqVcJiKisri+mk3YlESf0ks+mtTVy98Wp+9ubPkh2KiIiIyLDGsuKeyHr6biUlJdTV1SVs/ETQlpaTTG56Lo0djTy771k+dM6Hkh2OiIiIjMe9+UM/dsPXoTjyt77uAfjNPw0zTlPP7Yor4cArgx9XdBfc9A3n9v6XoPLqocdc+xTMvWDox2NQUVHB5z73uVH1iUc9/apVq6ipqRn2mDVr1oz6U4RkUlI/yVw691JcxsWWw1u0Z72IiIhMeKNNzmtra7ntttvG9Fy1tbX4/f7oav9IGhoaxvQ8yWCstcmOIeUUFxfbifyRzJ2/vZOXj7zM16/5unbCERERmcDefPNNzjrrrGSHkTS1tbXU1NSwYcOGAY8FAgFWrFjBhg0boqvlfr+fxYsXU1VVNa4VdGMMI+XAK1asYMuWLWN+jpHEOvfGmC3W2uKRjlNN/SR0xalXAPDMvmeSHImIiIjI0EpKSoZcNd+0aRMNDQ19trqsqKjA6/UOmtCXl5ezePHi6Im03VtljkV9fT3FxSPm0ROKkvpJ6PJ5lwPw9N6nR3wXKiIiIpJMJSUlVFdXD2gvLi7m/vvvp6ioCHAS7fLycqqqqgYcGwgEKCkpYceOHWPa976/sWyzmWyqqZ+EzvSdyYysGRxqO8T2wHaWFixNdkgiIiIigyorK6OsrGzA6ntRURF1dXWUl5cDsHnzZnbs2DFosu71eikqKiIQCFBcXDzui1LV1dWl3JaWSuonIZdxcdncy/jVjl/x9L6nldSLiIjIhNWdpAcCgQHJ+Nq1a2MeJxAIAE6CP9hYsaqsrGTNmjVj6ptMKr+ZpC4/tacER0RERGQiKy0tpbKycsz9/X5/dF/5QCDApk2bou2jVVFRMao3ExOFkvpJ6tK5l+I2bl46/BJNHU0jdxARERFJkrVr17Jx48Yx9Q0EAvh8PoqLi1m5ciWlpaXRLS/Xr1/P+vXro8d21+V3P9b/JF2/34/P5xt3+U4yqPxmkspLz6N4VjEvHHyBP+39EzcuvjHZIYmIiIgMqbi4mPr6+uiJsbHqnYD334KyqqqqT+JeVFREUVER69atG3SsVDxBtptW6iexaxZcA8CTe55MciQiIiIiwysrK4v5olCj0V1rHwu/358yV5DtTyv1k9jKBSs51HaIVQtWJTsUERERkWF1r6LHU21tLSUlJTEfP9h2malCSf0kNjtnNp9a8alkhyEiIiKSFKNJ6FOdym9ERERERFKckvoUEg5bXvAfoz0YirlPV7iLn735Mz755CcJhWPvJyIiIiKpQ+U3KeC1fU08XL+PR7fu51BzB5V3ruC6ZbNj6us2bn78xo/Zd3wfrxx5haJZ8a1VExEREZHkU1KfAjbV7eHHz78Tvf+bVw/EnNQbYyg7twyPy8MZvjMSFaKIiIiIJJGS+hRw43lz+yT1NW8coq2zi+z02Kbv5qU3Jyo0EREREZkAVFOfAlYsKGBOfmb0/olgiCf+cjiJEYmIiIjIRKKkPgW4XIb3LZ/Tp+3Xr+wf1Ri7mnZx73P38s2XvhnP0ERERERkAlBSnyJuPG9un/tPvnWE5vZgzP1bu1p5aNtDPPT2Q9oFR0RERGSSUVKfIs49NZ8Fvuzo/c6uMDWvH4q5/9m+s5mfO59j7ceoO1SXiBBFRERExszv91NeXp7UGMrKypL6/OOhpD5FGGO48bx+JTivxl6CY4zh+oXXA/C7Xb+La2wiIiIi47VhwwaKikbeerusrIzS0lKMMRhjKC0txe/3xyUGr9dLdXV1XMY62ZTUp5Abzu1bgvPMtqM0tHbG3P/6RU5S/4ddf6AzFHs/ERERkUSrra2lpKRkxOMqKiqoqqrC6/WyevVqqqqqKCwsjEsMZWVlVFRUxGWsk01JfQo5c3YuS06ZFr3fFbY88vK+mPufXnA6ZxScQXNnM3/a+6dEhCgiIiIyatXV1axevTrm4/1+P4FAgDVr1sQ1jsLCQhoaGggEAnEd92RQUp9CjDHcfMG8Pm0P1cee1APcuPhGAB7Z8Ujc4hIREREZj4qKilHVs9fW1gLEtLI/WmVlZVRWVsZ93ESbFEm9MabIGLPWGFPS/T2GPl5jzDpjTGHkdqExZkMsfZPp5gvmYUzP/a37mnjrYEvM/d9X+D7cxs3Te5+mob0hARGKiIiIxK67Hn40JTQ1NTUUFhbi9XrjHs/atWvZuHFj3MdNtJRP6o0xhcAGa22ltbbWWlsJlEXah+MDNgA7gEZgC7DZWlub2IjHZ643i8sWz+jT9lD93pj7z8iawaVzL6XLdvHYzsfiHZ6IiIjIqIx2lR5ir78fq8LCQurr6xM2fiKkfFIPlAH9z2iowEnYR7IKKAAWW2sLrLUpcbrz6hWn9rn/i5f20RUKx9z/psU3AfDrHb+Oa1wiIiIiozXWevpVq1YlLKZUPGF2MiT1q4H+b6XqIu0jstYGrLXx2QfpJHn3stlMy/BE7x9p6eDpbUdj7n/1/KvJTcvl9WOvs71xeyJCFBERERnRWFbcE1lP362kpIS6utS6rk9KJ/XGGC9QCPQpDrfWBiKPx2d/owkmK93N+5b33bP+wc27Y+6f6cmMbm/50LaH4hqbiIiISKzGUnoznnr6+vp6KisrKS8vH3F/+zVr1qTUnvWekQ+Z0HzQk8QPohAYbhW+0BjTvaLvAxpSpQSntPhUNtbtid6vffMwh5rbmZWXGVP/2864DY/Lwy1Lb0lUiCIiIjJOy3+0fFTHn+U7i003bhrQf+tdW6Ntt/36Nt5seHNU4w7W/8EbHmTZ9GWjGmcwo03Oa2true2220b9PIFAgLq6OtauXRsdZ9WqVezYsWPIPg0NqbOpSEqv1APjOeW5AcBaWx35qgTW9Ery+4jsqlNnjKk7cuTIOJ42PlacVsDSXnvWh8KWTZv3DNOjrzN9Z/L5iz/P0oKliQhPREREZETD1a4HAgEWL17cZ7V8PPX0fr+fDRt6TrksLi6OjjeYjRs3Rt8ApIJUX6kfs8jqfv9NSCsiXwNW6yNJfyVAcXGxTXiAIzDGcMfFC7j3129E237+4m4+es0S3C4zTE8RERFJFb1XyOPVv/dK/liMt39vJSUlrF+/fvDn2bSJhoaGPltdVlRURK8k2195eTkVFRXU1NTg8/kG9C0qKqKmpiZ6v66uDq/XO+gnBfX19RQXF4/npZ10qb5SD0Rr6+PBj1OSE/9NTxPg5qJTyUzrmcL9Te089dbhUY3x2M7HuOO3d+iEWREREUmKkpKSQWvXi4uLuf/++ykqKgKcRLu8vJyqqqoBxwYCAUpKStixY8ew+973f4Nw//33DxrTWGr9ky3Vk/ruenlf78ZeSfmQ9fTGmHWDNHcXTqXECbb5WWncdN7cPm0/fSH2E2YBNh/czKtHXuVXO34Vz9BEREREYjJUCU5RURENDQ2Ul5dTXl7Offfdx44dOwbd9cbr9VJUVEQgEKC4uHjEOv3KykrWrFkz5FaadXV10TcTqSKly2+stQFjjJ+BtfU+YMitKrsvWGWMqe53TPebg5TZ4vL2i09jU13PxaeefOsw7xxr5bTpObH1P/N2zplxDu9d9N5EhSgiIiIypO7V80AgMCAZH01Ne3dtvNfrHXSsbrW1tRQWFg65JWZ3wp9qUn2lHqAW6F/0VBRpH1QkkS8bJOkvAeqH2U1nwjnv1HzOmZcXvW8t/PC5XTH3X1KwhFuW3kKmJ7Zdc0RERETirbS0lMrK/qc6xs7v90f3lQ8EAmzatCna3lt9fT0+ny+a0A9W9lNRUZFSJ8h2mwxJ/XqgtF9bWaQdcMpxjDE7jDG9Z6ih9z72kZKdMuAjiQw23owxfOjSRX3aNm3eQ3N7cNRjtQZbCYVD8QpNREREJCZr165l48aNY+obCATw+XwUFxezcuVKSktLo1terl+/Pnoirt/vZ+XKlaxYsQJjDMaYASfp+v1+fD7fmPbAT7aUT+ojq+rrjTEbjDGrI7XyGwZZhff161cNFBlj1hljNgAbgFJrbf+r0054N5w3h5m5GdH7rZ2hUW1vCfD9rd9nZdVK/rT3T/EOT0RERGRExcXF1NePPg3r3sHG6/WyZcsWampqokl5VVVVdPvLwsJCGhsbsdZGv/rvUZ+KJ8h2S+ma+m6RRHzI34JI4l8wSHtKXGhqJBkeN3dechr/VfN2tO2BZ3fxwUsX4nHH9r7N4/LQGmzlp2/+lGsWXJOoUEVEREQGVVZWRm1tbdxPUB1qH/rB9N/LPpWk/Eq9OO64eAHpnp7p3Bc4we9fPxRz/5uX3kyWJ4sXDr7Am8dGd5U5ERERkfEqKipi3brBNiccu9ra2iFPiB3MYNtlpgol9ZPE9GkZ3Hz+vD5t33lqO9bGdp2svPQ8bl16KwDff+37cY9PRERE5GQrKSlJyfr4sVBSP4l85MpFmF4Xk319fzN/fPtIzP3vWnYXHpeHmndqeKf5nQREKCIiIiKJoKR+EllySi7vPnt2n7ZvPxn7lWJn58zmxsIbCdswD7z2QLzDExEREZEEUVI/yXzsmiV97m/e1ciLOxuGOHqgD5/zYQyGR3Y8wuG2w/EOT0REREQSQEn9JLP81HyuPH1mn7ZvPrEt5v4L8xey6rRVBMNBrdaLiIiIpAgl9ZPQx65e3Of+09uO8mf/sZj7rz3XuUbXprc2cbD1YFxjExERkb5i3dRCJo9EzLmS+kno4sLpXLyoz7W2+Nrv34r5F+gM3xlcv/B6OsOdVLxakYgQRUREBHC73QSDo78KvKS2YDCI2+2O65hK6iepz7z7jD73695p5Mm3Yq+R/+j5H8VlXPxy2y/Z0zy6q9OKiIhIbHJzc2lubk52GHKSNTc3k5ubG9cxldRPUsULfVx75il92r76+7cJh2NbrV+Uv4gbC2/krOln0drVmogQRUREpjyfz0djYyNHjx6ls7NTpTiTmLWWzs5Ojh49SmNjIz6fb+ROo2D0yzN6xcXFtq6uLtlhjOiN/c289xtP92n7xgcu4Kbz5sbUvy3YRpYnC9N783sRERGJq46ODhoaGmhpaSEUCiU7HEkgt9tNbm4uPp+PjIyMmPoYY7ZYa4tHOs4z7uhkwjp7bh43njeXX7+yP9r2td+/xbuXzSLDM3IdV3ZadiLDExERESAjI4M5c+YwZ86cZIciKUzlN5Pcp1adjtvVs9K+u6GN7z+zc1Rj7G3Zy6ee+hTP7Xsu3uGJiIiISBwoqZ/kFs3I4faLFvRp+9YT2znU3B7zGDXv1FDzTg3feOkbqvUTERERmYCU1E8Bn1p1OvlZadH7bZ0hNjz2l5j733HWHdx59p18/Zqvq75eREREZAJSUj8FFOSk88/Xnd6n7eGX9rHlncaY+qe701l34Tpm58xORHgiIiIiMk5K6qeI2y9awJmz++6H+sVHXqMrFB7VOF3hLp7f/3w8QxMRERGRcVJSP0V43C7+341n92l7bV8zP3g29pNmQ+EQH/rdh1hbs5b6Q/XxDlFERERExkhJ/RRy6eIZ3HBu3+2y/vMPb7PraGwXl3K73Fw852IA7n3+XjpDnXGPUURERERGT0n9FPPFG5f1OWm2oyvMZx9+NeZdbT5y7kdYmLeQnU07uX/r/YkKU0RERERGQUn9FDMzN4Mv3NC3DOfP/gYe3Lwnpv4Z7gzuvfReAP536/+yrXFbvEMUERERkVFSUj8F3Vo0jyuWzujT9u+PvsnuY20x9V8xawWlp5fSFe7iX5/9V4KhYCLCFBEREZEYKamfgowxfOXm5WSnu6Ntxzu6+MSDLxGMcTecT674JHNz5vLGsTf4zivfSVSoIiIiIhIDJfVT1HxfNp97z5l92l7eE+Abj8dWTpObnst9V9yHy7j4/tbvs/ng5kSEKSIiIiIxUFI/hf3NJadRctYpfdq+/eR2XvAfi6l/0awi7l5+NxbL55/5PE0dTYkIU0RERERGoKR+CjPGsOHWc5mZmxFtC1v4xwdf5ujxjpjGuOe8e1g+YzkHWw/ypee/FPMuOiIiIiISP0rqp7jp0zL4r9vO69N2sLmdj/20Pqb6+jRXGvddcR85aTnUvFPDj17/UaJCFREREZEhKKkXrlg6k7IrC/u0vbCzga/89s2Y+p+Wdxr/ftm/MztnNitmrUhEiCIiIiIyDE+yA5CJ4TPvPoOt+5p4bkdPPf0Dz+5i+bx8bik6dcT+K09byWXzLiPTk5nIMEVERERkEFqpFwA8bhff/MAFzPNm9Wn/7ENbYz5xtndC/8TuJ2gLxrbvvYiIiIiMj5J6iZo+LYPv/c0K0j09vxadoTAf+XEd2w61xDzOg395kH988h/59B8/TdjGtu+9iIiIiIydknrpY/mp+Wy4dXmftub2Lj74wGYONbfHNMYlcy7Bl+nj8nmX4zL6FRMRERFJNGVcMsDNF5zKp687vU/bvsAJ7vrBizS2do7Yf2H+Qh69+VFuP+v2RIUoIiIiIr0oqZdBfeyaJXzgogV92v5ysIU7f/ACTSeCI/aflj4tent743aq366Oe4wiIiIi4lBSL4MyxvDl9y9j5Zl9rzj72r5m/vYHL9LcPnJiD9DU0cTf/eHv+NLzX6Lq7apEhCoiIiIy5SmplyF53C6+dXsRlxT6+rS/sifA337/RQJtI5fi5Gfks/bctQB8+fkv88vtv0xIrCIiIiJTmZJ6GVZWupvv33UhFy4s6NP+8p4At1U8z8GmkU+eveOsO/jUik9hsXzh2S/w0zd/mqhwRURERKYkJfUyopwMDw986CKKFnj7tL996Di3fvc5dh5tHXGMD53zIT5d/GkA/uPF/+C7r3wXa21C4hURERGZapTUS0ymZXj40Ycv4uJFfUtx9gVOcMt3no3pAlV3LbuLL136JVzGxXde/g73vXgfoXAoUSGLiIiITBlK6iVmuZlp/OjDF7Hq7Fl92hvbgtzxvy/w4Iu7RxzjlqW38LWrvkaaK42f/+XnfPyJj3O883iiQhYRERGZEpTUy6hkprn57h1FlK44tU97V9jy2Ye3cu8jr9PZNfxVZFedtorKVZV4M7w8s+8Z7nzsTva27E1k2CIiIiKTmpJ6GTWP20X56nP51KrTBzz2w+d2Ufq959jT0DbsGMWzi/nZe39GYX4h2wPb+etH/5qXD7+cqJBFREREJrUxJ/XGmFuMMXfHMxhJHcYYPrFyKd+5o4jMtL6/Rq/sbeK933ia3249MOwY8/Pm85P3/oQr5l1BhjuD0/JOS2TIIiIiIpOWGc0OJMaY84FSTLMAAAAgAElEQVR7gBXARsAAa4DNQIW1dkostRYXF9u6urpkhzFhvLavibL/28K+wIkBj9103lzuvWkZvpz0IfuHbZhDrYeYM20OAB2hDg63HWZ+7vyExSwiIiKSCowxW6y1xSMdN+JKvTEmzxjzGWNMHfA5oMpae6G19mvW2q9GnqQSuMcYs9kY82ljTN74X4KkinPm5fPoJy6n5KxZAx575JX9XPfff+R3rx0csr/LuKIJPcC3X/42tz5yK7/1/zYh8YqIiIhMNkMm9caYW40xfwAeBxqttcXW2jXW2sf7H2utfclae4+19kKgCXjCGPN7Y8wtiQtdJhJvdjr3/+0KvnDD2aS5TZ/Hjh7v5J6fbOEjP65j97Hha+2ttRxpO0JHqIN5ufMSGbKIiIjIpDFo+Y0x5ntAI1Bprd05poGNWQSUAYustWvGFeUEo/Kb4b2xv5lPV73CGweaBzyW7nFxz5WF/P3VS8hKdw85xvbG7SwpWBK9/9jOx7hm/jVkejITErOIiIjIRBRr+c2oaurFoaR+ZMFQmO88uYNvPrGNrvDA37E5+Zl8YuVSVq84lTT38FVgz+17jrLaMuZNm8dnLvwM186/FmPMsH1EREREJoO41dSLjEWa28U/lizlkY9fzgULvAMeP9DUzuce3sqq//ojv3xpH6FBEv9u09KncXrB6ew7vo9/evKfuOt3d1F3UG+qRERERLqNaaU+ciJsIRAAGqy1A+ssJjGt1I9OOGx5+KV9/Mdjf+Ho8Y5BjymckcPdVxRyS9E8MtMGluV0hbvY9NYmvvfK92jsaATgsnmX8Q/n/wPLZixLaPwiIiIiyRLX8htjzEeAUqAYsIAfp+YenOS+INK2Eai21u4aW9ipQUn92LS0B/nmE9v54XO7hrzq7PScdP72XQu5812nDboN5vHO4/zfm//Hj17/Ea3BVgAunnMxH172Yd41910qyxEREZFJZdxJfWQ1/vPABUAVzlaWTSM86Uqc5L8AZ9/6J0YbeCpQUj8++wMn+OYT26mq2zNovT1AutvFu8+ZzV9fOJ93FU7H5eqbrDe2N/LAaw+w8a2NtHU5O+qc6TuTO866g3cvfDdZnqyEvw4RERGRRBtXUm+MuQBYC5SPY/ebjwD51tqvjaX/RKakPj7eOdbK/zy+jUde3j9kcg8w35fFbSvmc+N5c1k4I6fPY82dzWx6axM/eeMnHGs/BsDHz/84ZeeVJTR2ERERkZNhzEl9ZCvKImvtQ3EIIh9Yaa19eLxjTSRK6uPrQNMJfvjsLn72wm5aOrqGPfaceXnccO5c3rd8DvN92dH2jlAHv/X/lupt1fznVf/J7JzZAPxq+69o62rj+oXXU5BZkNDXISIiIhJv2tIygZTUJ0ZLe5AHX9zDT194h10jXKQK4MzZuVxz5ilcc8YpFC3w4um3Naa1lpt+eRO7mndRUVLBpfMuBZyTbj0uT0Jeg4iIiEg8JSSpN8acb619eVyRTQJK6hMrHLb8eecxNm7ew2OvHRzypNre8jI9XL50BpcUTueiRT5OPyUXS5jf7fodf9z7R75y+VeiifzHH/84gY4AV516FVfNv4ql3qU6wVZEREQmpEQl9d+11v59v7Z8nPr7mqmS8CupP3kCbZ38+tUD/OaV/by4q4FYf1292WlcuNDHxYt8XLDAy9lz8slKdxMMBblq41W0BFuix87KnsVFsy/iwtkXctGci5g3bV6CXo2IiIjI6MQlqTfGbMfZurIGqAVK+yf1vY5dCdjJuuNNb0rqk+NQczuPbT3Ao1sPsOWdRoY5t3YAl4HTZ+WyfF4+Z8xNJ5zxFv62zTx/4Bka2hv6HDtv2jwuOOUCls9YzvIZyznDdwbp7oHba4qIiIgkWryS+guANUAJUISzR309ToJfA9T1vvCUMeaWyXZS7GCU1CdfoK2TP759hKfeOsJTbx2msS04pnFm5aUzf1YTmXk7OeF6m73tr9HWdbzPMdfOv5b/ufZ/AGgLtrEtsI0l3iXkpOUMNqSIiIhI3CSq/OYPOMn8KpxEvzvJrwO2AKustWvGFPE4GGOKcC6M5ce5GJbfWlubqH5K6ieWUNiydV8TL/iP8cLOBjbvaqClffhddIYWxpV5gIycPUzLOwAZu1mcdRWr5t3OPG8mR7pep/yVT3LujPP46ft+AkBnqJNH/Y+yIG8BC3IXMCNrhmr0RUREJC4SldR/z1p7T6/7RcBK4MJI033W2pdGG+x4GGMKcS50tapXWxWw3lrrj3c/UFI/0YXClr8cbObFnQ1seaeR1/Y1xbSbztAs4CTp7py3yTjlMUJti8htLWV2fgZ5uY1sNf8SPTrNlcEpmfOYlzOf+Xmnclr+XObnzWVW9ixm5cxieuZ03C73+F6kiIiITAmJSupXWmsfH1dkcWaM2QBsttZW92orAcqstaXx7gdK6lNRU1uQ1/Y3sXVfE1v3NvH2oRZ2Hm0d9qJXsTLpR8iY8Tiu9GOYtGO4PCO8gbAuZrGKszLuYFqGB5PWwN7gU8zJXkjxjGvJzfSQne4mSBPTs3xMy8ggK91NpscV+e4ecIVdERERmZymzD71xpgdOGU//l5tXqDRWjtk5jPWfqCkfrLo7ArzzrFWth0+ztuHWth2+DjvHGvlnWNt4yjfAVxtuNIbcKUfxZUWwHgCmLQmXGlNGE8TLk8rHUdK6DxaAoA75y2yFzxA1/GlnNjzd5ExTpB7xpcAsOE0bCgTG86CyHdjM3HbbNJMz1eWazZezibN7cLtDhN07yXTlUOeZy5pbhfpHoPHZUj3uJ37bkOa24XH7SItctvtMs6XMbhcBrcLXKZfW+S+K9pGnzaPy/T0cRlcBowxGHq+u4yhu0LJmJ77hu7vkWN73XZFHsf06hPpH+3Xb4zouCqHEhGRFBVrUj/gCjyRK8peEI8TXo0xecBt1tr/He9YQ4zvxamF77N9ibU24CQEpnCwUpqx9pPJJd3jYumsXJbOyuW9y+f0eaypLcjuhrY+X4ea2znY1M7hlnaOHu8ceuBwNuH2bMLtpw7+uAkCPXvvh4M+Oo6UEA72XPHWuE8Q7srBuNswriDGFQRaBgwVjHwBHG05m217nddhPE1MW3of4Y48Wl//fPT4nCVfwbjasTYdrBusBxv5TrjXbevGWg+2K5eOQ+/v+ZnNqAUTpvPo1WCdHYE8uVsxnibABdYFuLC253aftsj9UGsh3f/9uNIPg6uDcOdMCGdG4m/GuFsBA9YABhspgeq+3/sxMNhQNti0yA+wA+PqwobTe9oIYVxdPX0xmD5rGmZA8h95CxK55xzsHGL6HdOn03B3Gez9Rf8x+h8zcIz+cQ4YcJDnGGEMve+RCUi/ljLRfOyaJdx9RWGywxjUgKTeWrvTGOM1xnwX2GCt3TWWgY0xHwEKrbWfG2eMw/GBk4wP8Xghzkmw8eonU0R+dhrLs/NZfmr+oI93doU5cryDg03tHGp2vhrbgjS2dtLQ1ul8b+2ksc35Hgz1yh6jSWbkbufM6Kp9tC3oo3XbFwALphPjbse4T2Bc7eBux7hOOG0upx1XO+GO3m9MLKH2udiuvjv0GFcHxt2JYZg3Jb14OnNZdvgs0k0X6QR51fdHQu4g4WPvoiuS1J9SUEtLzqGYxuv28R1LyAkb3IT5+Zx97Mxup2D3+9nd+i4AzvZtZPf0HaMa8x/3ZXH2CTcuY3nI18HvCrooPHImrxz9IADn5D/GO3OfGdWYdxwLc0vAmbvncgz/OdvF0pY86vc6b5Rmp22ndUnPmoWJfPLZnYiEcUXfjLgi9y47HmbdYeet2EEP/P38DKZ3udi2477oOPOXrudE3wskR8cN4qHLOudkuE2YTILMDoap3NcaPeYDp06j1WVo3b6OlvB0AC6e/0V2Z7UP+jq7cHOC9Eh/yzSc4366p4VpkfefX5iVxauZHnz7ruf11msAuPyUb7LDu3fQMcMYWsmM3p/GCQzw/x1q5fz2EAA/LMjgF3npLDx2Os8duxuAorxHODh76HlqITv6M82mHQ8h/raxg1ubnd/pZ7I9fHVmFota83l2n3Ouy5y07ZhFlUOO2UomIZyfaQadZBDk0rYu1h85AcABj+GeedPwdbl4zf8f0X5LlqznhGvwT7vbSacj8m8kzXSRTTundIWp3Nezu9Yd83NpdRmat6+PztMl87/InqwTg44ZxEObdX6mLsLkGqfU7yd7WpgWKSX8wqxstmZ68O59D6+3OfN0xcxv4i/YM+iYYVy02Ozo/TzTisHypUOtnBeZpx95M/hFfgYLj53Bs8ecTxNX5P2Kg7OfHfwHCjTbHMKReZpmTkTmqb3XPKU583Q8n2f2O/M0N307ZuEw82SzCEbmKct0kkEnl7UGWXfU+Xkd9Lgoi8zT1l7ztHTJetqGnKcM2iP/J6fTRbZpZ3ZXmMp9PYsot8/Po9VlaOo1T+8abp6sh+NkAeAmTJ5x/m3+ZE8zuZF5+n+zcnglMk+vRebpyhHmqcn2/H/uNcej83R+5JPlH3kzeTg/g4XHTueZyL+nFXmPcGiYf09Ndlp0nnJNW6956gC65ymbhf3myTXMPB23WQQj6WUWHWSa7nlyfl8PeFzcMy8XX9DFqzt75un0YebphM2gPfJ/VDpBckw7s5I0T+3B0JCvPdkGJPUA1tqXjDGfBT4X2dayCtjUe/vKwRhjzgfuARbhvCG4P94B9+M9Wf2MMWtxLrLFggULxvi0Mlmke1zM82Yxz5s14rHWWlo7QzS2dtLcHuR4exct7V0c7+iipaMrcj/I8e7bHV20dnTRHgxxIhimIxjiRDBEZ2cnnq7jpLe3cNQu5EQkYSoyb3OeawdZHCLLs4lsOsiiney9Lpqsiy/2iuVnuzrIcbVgXJ1gQmBCdBkIGvjf0HU8Yi8BE+Ii12t8MO1R0jnCNRn3Rvv/KJBLm8vwHRuKfm5wQ2sAd1czXcYQAkK9vu+209lilwJhsswJLnG/TsjAne4nyY78572jq4CcjnTabQe7I2POD58go7MTi5McWpx18rCBIG722RlAGAzMNscAS5HZRZHb+UP0Enl4Q7n4bEc09mzbQVY4HB0z3G8JsAs33X9O3JFnPcU0cbrLeaV+VzYwg8xeb4jcpu9/7jay3N3zZ8lGonfaQhgyTTunuY4B4HJ5CLrmOnPROxZj6XANktVDZBTn+DDQhqHTFWaW6VmjOOGaRqvb0PsToZCri+PuodY9wxh6Ev7jkT/y000T+ZGPM0KudI6705hperaPdZnOYcakz5it0QSvlZnGmRe38dLizsRten6mHhOkxT3Uawdoj75pcv5Uu0hztTPTOPOU5cqmxZ2D7RWnMaERxuyMjtkJdOLCmiAzTZPzjMZDizuPnHDfq1u3uizHhxy3C0NX5BY04yIvHGKm6fkzetyVF4mrZ9wuVxfNQ44ZxuAkRjYyJoDPNEfnqcuVSbPbxXRXr9fv6hxmTKJjArREPsXKMW3ReXK5vDS7s3CZnn9PbjNcnAAnop9xOaO78Lg6mBGZp0xXNs3uHMKu3tsRh0YYsyM6ZkfkXsjVxYzIz/SE8dDsziO73zwdH3aegrgin3f2nqfppidZbHHlD5inoKuLpmHmyYWTIFqgKRJ1gWmJzlOHK4smtwtf79fv6hxmTKJjAjRH5inbtDG9e15cbprcWRjT+/+o4LBjQlv0Z9oaeRa3qyP6+jNc2TS5pw2Yp+HH7MBFR+SWc6/L1RUds824aXLnk9lvnlpGMU9NuJgWDidnniawmGrqjTG34uxX371XfYCe0pXFQAHOynYdUHmydsCJ7L6zZbAaeGOMxamZH7BF5Vj7dVNNvYybtdDWAK1HnK+2o5BVAIVXO48374fqv4MTjdDRDO1N0NmzwtfxgYdoO/VKTgRDZP3p3yio/9agT3Miaza/v+4JOkNhgqEwt9ReRVawsW8oGLpMGo/PuZs/zridcNhS2PIiNx38Fl146DJpBE0aXXgIkkbQeLjf9xnaTBahsOX65mrmdO2lCxdd1k0IQ9C6CeFmh2sBj3uuBAtZ4eO8r6uGkHURwhU5xrkdxsWfzbkcwYcFloR3sdDuJYTBWgjhinw3nLAZPM+5WGsJW7iYVzFYLIawNYQjx4WtYY+dyX7rrNjkc5zTzKFIuu78QQzT81/A2/bU6GrtfHOI7MgfJdur9MfirEIdxBkzjSDzzaE+P0sisQAcst7omy8vzeSZNpy1dHfkeIs1IULWzQE7OzrOHHOgz5uB3uMfN5kcx1lZzaSdvO6ELNyzKh52dQCWY+GZ2MjaTb7rKGkM/oepAw8tOKuAbsIUdL9lC2dES4NsZMzmcAGdkdc0zTSRaQZfBQthaCQ3en96JA0hnIbpfv0mCKaLNjuNNuscm2HayI0k04M5Ri428sc3j1bS6QLrwURWW63pAhOkw2bSYgsirylIgevokGM2k00nTv8c2smiA6wbE1lpt4TB1UEIN43hU6L9fK6DGDP439A2MjhBhvOaCDKNE2ANJtyzEGBdJ8BYGkI985TrOka6GfyTtE48tETm3k0YL5H/E0JZfefJhGgJeemMrELmmCYyXYOfwB/CEOg1T77ueQpl9JqnTnB10RbO4YTNc16TaSPHNfQ8NTKN7tK1abSSRsiZ+97z5OokaDNpCTvz5KILr/vIkGO2kE1XdAW43XmDbd2YsPNztoTB3U7YugiEZ0X7FbgPRv899XeCdNojv8/pdDr/7vvPk7sNsAT6zZNniHkK4qE18rM3hMnvTsZD2b3mqR1MiOOhAoKR5882TWQMMU9hDM1Mi97Pj7z9IpTZb56CnAhPoz3676mVbNfQ67EBcuiZJ2elnnB6v3nqIGgzOd5rnvKHmafjZNPV/YkKHWTQ6fwb7TNPJ5xPH0I98+Qddp4y6Oi1Up9NO2AwoZ5PmU7WPH36utP5+LVLh3z9iTDmmvrBWGsfAh7qNXg+TolKAGiw1g79L/skMMZ4hymliXs/kRF1HIfGndD4DiwpgbRIwvW7z8NbjzpJe6jffzRLVvUk9e4M2P1cv0ENZOZBZj4ZbkNGTjoFAGdeCZ5OSM+BtGznKz0b0nLIyvLyV6fP6xli6ZPgToe0LPBkgDsD4/aQBlwf+XKcC9w95Mvr+z/LJcP+KP6lz70bhj12bK4f8QhrLf3XL/r/6ei/wDHYn5aBY9hhHx+pfyxjjDbOQWMYZ9wiyaBfS5mIstIm7pbUMSX1/UWS+JO6H/0QuuvefThvMIDoibC9H49XP5GBmvbBKz+Dhp3Q4He+jveqMf/YZph5unP7RCM07nJuZ+ZDzimQM8P5mlvU0yerAO76tfM90+scmz4NBivJOP3dzlcsfIvG9BJTnem1284wR52MUERERBJiTEl9f8aYhWM9oXY8IrvV+BlYI+8DAkPtYDPWfjKFNR+AAy/Dodfh8Bsw+1y4/J+cx040whP/1vd4dzoULIKC0+iz3nTVOrjinyFvjrOyPhSXCxZdGfeXISIiIpNTXJJ6oNEYk2etbY6U5pQSw4m1cVKLUw1Q36utKNKeiH4yVbz8c/jLb2BfPbTs7/vY8cM9Sb1vEVz6CfAV9nzlzYXBrho7RVfKRUREJLHildQXAIuNMZuttU3GmCqcHXNirAkYl/WR5+q9v1JZ5AuIltVswdmRpzLWfjJFdLbC7j/Drqfh4nsgN3LC4p4/O0k9QEY+zD0fZp0Ds852Vuq7pefAdV8++XGLiIiIRMQrqa+x1i41xuRHLjgFRLaHSLBIKc16Y8wGYDPOCbwbBimh8Y2xn0w21sLRbfD2Y/D272HPixCO7Aoy8yw4b41z+/w7YP4lMG8FTF8yeD27iIiIyAQQr6S+0hhzvrX25Uj5zQpgZZzGHpG1tp6+ZTT9Hw/gfJowqn4yCXV1wncvhWPbetqMC+ZeAAuvcFbhu82/yPkSERERmeCGTOpHc/KrtfarkWSeSPnNTrQblSSbtXDwVXj7D3Dlp8EY8KRDzkxnX/il18Hp18PiayFrrNcxExEREUm+4VbqN+BccComvfeqt9bu7E7yRU66tgZ45UHY8kM4+pbTtvgaODWyu3rpA5A9A9zx+qBKREREJLmGy2pKI1dX3Qg8PtqdbJJ9QSqZYqyFvXVQ9wN4/WHoilyePnsGLLvZ2e+9W+7swccQERERSVEjLVVeB9wG2Mi+7rXAHxhDki+SUMET8NPV0B65ltjia6H4w055jTstubGJiIiIJNhwSX21tfY2Y0whUBL5WoOz5WOfJN9a+4vEhyrSS6gLtlbB2Tc5W0qmZzsXdWo7Cis+6OwVLyIiIjJFDJfUNwBEtnisjHxhjFkErKJXkm+MaQRWW2ufTGy4MuWFw/DGL+HJf4dj26G9CS65x3nssk8kNzYRERGRJBkyqbfW3jNE+04GJvn3ALXGmEJr7TuJCFQE/1Pwhy84O9oAFCyC/FOTGpKIiIjIRDDu7T8iSf56Y8xGoJxR7JgjEpPAHvjDv8Abv3Lu586Bq9bBBXeqXl5ERESE+F18CmttvTGmLl7jiUTtesZJ6NOynf3mL/kopGUlOyoRERGRCWPcSb0x5hbAb619GV1wSuKlvQkyI5c6OHcNNPhhxV0qtxEREREZhCsOY5QDj0eSe+1NL+MTCsIfy+Hry+HYDqfN5YJr/0UJvYiIiMgQ4lF+swJnL3u/tfbhOIwnU1XDTqj+MOyvd+5vr4Xpi5Mbk4iIiEgKiMeJsk3A/XGIRaayN38Dv/wodDRB/nx4/7eg8OpkRyUiIiKSEuJ2oqzImIS6oPaL8Py3nPtn3gDv/zZkeZMbl4iIiEgKGbSm3hjzGWPMteMd3BhzvjHm0+MdRyaxI3+BFyvB5YF3fwXW/EQJvYiIiMgoDbVSX42z93xF5HaFtXZXLAMaY/KAMpz96jcTuUiVyKBmnwM3f8/Ze/60S5MdjYiIiEhKGjSpj1xQ6h4AY8ytQKUxxuIk94OeDBs5rgxnW8sqa21xYkKWlLevHtoDsDjyYdA5tyY3HhEREZEUN2JNvbX2IeAhY8wioMwYUw7UABVAACf5XwnUAmWRNwQig9v5NPz8r8GG4e5amLUs2RGJiIiIpLyYT5SNJOufBT5rjFkJfB5nVb7SWvvZBMUnk8lffgtVH4RQBywvhRmnJzsiERERkUlhTLvfWGsfBx6Pcywymb35a9h0F9gQXHg3vOerzkWlRERERGTclFVJ4r39e6j6kJPQX/ZP8N6vKaEXERERiSNlVpJYJwLw0N0QDsK7Pg4l94IxyY5KREREZFLRxacksbK8UPoA7HgSrvs3JfQiIiIiCaCkXhIjHO4psVlS4nyJiIiISEKo/Ebir/UY3H8NbK9NdiQiIiIiU4KSeomvYDs8eDsceBke/7KzYi8iIiIiCaWkXuLrsc/Anj9D3jz4wM+1y42IiIjISaCMS+Jnyw+h/sfgyXQS+ry5yY5IREREZEpQUi/xsXcL/PYzzu0b/hvmnJfceERERESmECX1Mn5dHVD9QQh1OleLPf/2ZEckIiIiMqUoqZfx82TAzZWw7BZ4933JjkZERERkytE+9RIfp73L+RIRERGRk04r9TJ2R7fDtppkRyEiIiIy5Smpl7EJBeHhu+Gnq+GVjcmORkRERGRKU1IvY/PMf8P+lyB/AZxxfbKjEREREZnSlNTL6B15G/70Vef2X30bMvOTG4+IiIjIFKekXkYnHIZff8LZvvKCO2HRlcmOSERERGTKU1Ivo1P/Q9j9POScAtd9OdnRiIiIiAhK6mU0wiF4/jvO7feWQ1ZBcuMREREREUD71MtouNxwdw28WgVn/1WyoxERERGRCCX1MjpZBXDx2mRHISIiIiK9qPxGRhYOw5YfQbA92ZGIiIiIyCCU1MvIXt3o7Hjz45vA2mRHIyIiIiL9KKmX4XUch9p7ndsrPgTGJDUcERERERlISb0M7/lvwfGDMG8FnLsm2dGIiIiIyCCU1MvQ2hrguW85t1d9GVz6dRERERGZiJSlydCe/Tp0tsDia2HhZcmORkRERESGoKReBnf8MLxQ6dy+9l+TG4uIiIiIDEv71MvgsmfAX30H9tY59fQiIiIiMmEpqZfBuVxwzi3Ol4iIiIhMaCq/kYE6W5MdgYiIiIiMgpJ66avlEPznWfCbTzpXkhURERGRCU9JvfT1529DR5Nzoqy2sBQRERFJCcrapMeJAGz+gXP78k8lNxYRERERiZmSeumx+X5nX/pFV8Kp2vFGREREJFUoqRdHsB3+/D3ntlbpRURERFKKknpxvFYNbUdh9nIovDrZ0YiIiIjIKCipF8fBrc73i/8ejEluLCIiIiIyKrr4lDjeswGKPwze05IdiYiIiIiMkpJ66THzjGRHICIiIiJjkNJJvTGmCCgG/EAh4LfW1o7QxwusBaqBBsAHlAE1I/WdlJr2QdNemH+Rym5EREREUlTK1tQbYwqBDdbaSmttrbW2EiiLtA/HB2wAdgCNwBZg85RM6MHZxvIH18GT/57sSERERERkjFI2qcdZXa/o11aBk7CPZBVQACy21hZYa6vjHVxK6OqEl37i3F5SktxYRERERGTMUjmpXw3U92uri7SPyFobsNb64x5VKnnrt9B6BGaeCfMvTnY0IiIiIjJGKZnUR+riC3Fq4qOstYHI4yOV4AjAlh8631d8SPX0IiIiIiksVU+U9UFPEj+IQpyTZ4dSaIzpXtH3AQ0jleAYY9binGDLggULRhftRNTgB/+T4MmE89YkOxoRERERGYdUTeq94+jbANA7iTfGVBljGC6xj5yIWwlQXFxsx/H8E0P9j53vy26GrILkxiIiIiIi45KS5TfjEamlr+zXHOsJtpOHJwuyfLDig8mORERERETGKekr9aOsf2/oXXJjjPEOU4IzGn6ckpx4jTfxXb0eLv8kuNOSHYmIiIiIjFNSk/ruveZH0WUzUE5PvbwP6JPkR24OWU9vjFlnrS3v10/UkaIAAA0wSURBVNx9wm0hA3fUmbw86cmOQERERETiIKlJfWRLydIx9AsYY/wMrK33AUNuVdn9JsIYU93vGF/k++Tf4rK9Gd74JZz9fsjMT3Y0IiIiIhIHqVxTXwsU92srirQPKpLIlw2S9JcA9VOi9OaNX8Ej/wCb7kp2JCIiIiISJ6mc1K9n4Cp/WaQdcMpxjDE7IttRdmvoXccfKdkpAz6SyGAnjFd+7nxfHtM1ukREREQkBST9RNmxipTgrDfGbMCptS8ENgyyCu/r16/aGLM6sk/9dJwSntIpcXXZxnfgnWednW/OuinZ0YiIiIhInKRsUg9gra1nmBNbI+U0AzZhH+lCU5PW6w873894D2TmJTcWEREREYmbVC6/kdF6/RfO93NuTW4cIiIiIhJXSuqnimM74MArkJ4LS0qSHY2IiIiIxFFKl9/IKLQHYP4l4FsEaZnJjkZERERE4khJ/VQxbwX83e8hHEp2JCIiIiISZyq/mWpc7mRHICIiIiJxpqR+KtjxJOx/6f9v735i7LruOoB/T+JgTJt2MvkjQlqHjhsoC1rhuEENoCDVw4YdsmsQCMQiniVsiJUFa2RLSLBBmmSHkFAVr0BUojaopaUSxB4kglBLmolLStKoiT1pFPKnaQ6LeycZ3ryZeWP7zb33zecjje7Meffd97Nk3fd95/3uuUmtXVcCAMAUCPX7wZf/OHniV5PVr3RdCQAAUyDUz7prV5KXn2lWvbn/4a6rAQBgCoT6WffNLzXbBxaTAwe7rQUAgKkQ6mfdN/+u2X7q17utAwCAqRHqZ9kbryb//Y3kltuamXoAAGaSUD/Lnv37pL6XfOJXkh//aNfVAAAwJUL9LHvnjeTQvNYbAIAZ546ys+yhR5MHfz95792uKwEAYIqE+ll364HmBwCAmaX9ZlZ975nkrde6rgIAgD0g1M+iWpMv/k5ybiF5+T+7rgYAgCkT6mfRq881d5I9eHty9892XQ0AAFMm1M+ib19otkc+n9xya7e1AAAwdUL9LHr2y832gV/rtg4AAPaEUD9r3nkjufLPSUryyc93XQ0AAHtAqJ81z38t+dHbyX1Hkw/d1XU1AADsAaF+1nzvmWb7ycVu6wAAYM+4K9GseeSPkgd/L0npuhIAAPaIUD+LPnxP1xUAALCHtN/Mkh++1XUFAAB0QKifJX/7B8mf/Xzy3D92XQkAAHtI+82sqDV5/qvJ6y8lH/7JrqsBAGAPmamfFa/8VxPoP3RPcs/PdV0NAAB7SKifFatfabYLjyTFyjcAAPuJUD8rVr/abD/xSLd1AACw54T6WfDej5IrX29+XxDqAQD2G6F+Frz8H8nbryVz9ydzh7uuBgCAPWb1m1lw188kv/s3yds/6LoSAAA6INTPgtsOabsBANjHtN8AAMDACfVD98qzyV//VrLyl11XAgBAR4T6oXv+n5JvfemDdeoBANh3hPqh+843mu39D3dbBwAAnRHqh6zWDaH+l7qtBQCAzgj1Q3btSvL6i8lP3Jnc/amuqwEAoCNC/ZCtz9If/lxSSre1AADQGaF+yF74l2b78V/stg4AADrl5lNDdv/DyZtXk5/+5a4rAQCgQ0L9kH3mN5sfAAD2Ne03AAAwcGbqh+rK15N3304+/lBy8PauqwEAoENm6ofqa3+a/NVvJN++2HUlAAB0TKgfovfeS757ufn9Yw91WwsAAJ0T6ofolW8lb7+WfOS+5KP3dV0NAAAdE+qH6IV/bbYf+2y3dQAA0AtC/RCth3o3nQIAIEL9MH336WZrph4AgAj1w/PuO8m7byW33Jbc++muqwEAoAesUz80B34s+cN/T95cSw4c7LoaAAB6wEz9UB2a67oCAAB6Qqgfmh++2XUFAAD0jFA/NH/xueTPP5OsvdB1JQAA9ISe+iH536vJteeTA4eS2+/tuhoAAHrCTP2QvLjSbO/9dHKrz2MAADRmIhmWUk4kWau1Xpxw/6NJjiVZTbKQZHXS53bqxX9rtj/1C93WAQBArww+1JdSjid5MsnJCfdfSHK21rq4YeypUspqrXV1SmXeHP+zHuqPdlsHAAC9Mtj2m1LKQillOc1M+9VdPHUpyfLI2HKSszertqlZb7+5T6gHAOADgw31tdbVWutSrfWJXT71RJKVkbFL7Xh//eCl5PWXkoMfSeaPdF0NAAA9Mvj2m90opcxlzMx+rXWtlJJSysJWLTillNNJTifJ4cOHp17rJofmkt8+n7zx/eSWwX4WAwBgCvZVqE8ynzQhfovHF9JcPLtJ+43AE0ly7NixOpXqtnPboeSBxZ33AwBg39lvU75zXRcAAAA3234L9QAAMHM6b79pl5ic1NVtWmd285pzN+M4AADQB52G+vU143fxlKeTnLuBl1zvl59P8n6oby+g3fg4AAAMRqehvl1pZqKbRt2k11srpaxmc2/9fJo70gr1AAAMzn7sqb+Y5NjI2NF2HAAABmdWQv18xqxsU0qZK6U8164xv+5MNn87sNSOAwDA4HR+oez1avvgH0+ztvxckrOllMUkF2qt5zfsOr/xeW0LzplSytk0PfoLSc5qvQEAYKgGG+rb1Wu2nV1v97ljzPhKkpUplQYAAHtqVtpvAABg3xLqAQBg4IR6AAAYOKEeAAAGrtRau65hcEop30/ynY5e/q4kr3T02gA3wvkLGKouz1/311rv3mknoX5gSimXaq2jN88C6D3nL2CohnD+0n4DAAADJ9QDAMDACfXD80TXBQBcJ+cvYKh6f/7SUw8AAANnph4AAAZOqAcAgIE70HUB7KyUcjTJsSSrSRaSrNZaL3ZbFcB4pZQTSdZGz1OllLkkp5OcT3I1yXySpSQXnNOArpVSjidZTPJqkiNJLtdanxjZp7eZTKjvuVLKQpKztdbFDWNPlVJWa62rHZYGsEn7pvhkkpNjHp5Pcrb9SZK1JI/25Q0R2L/ac1dqrWc2jF0upczVWs+1f/c6k2m/6b+lJMsjY8v54E0RoHOllIVSynKamaur2+y6mOSOJEdqrXfUWs/vSYEA21saM3ZxZLzXmUyo778TSVZGxi614wC9UGtdrbUujX5VvcW+a32Y1QIYsThmbG3D773OZNpveqztP90061VrXSulpJSy4I0RAODG1FrHtQyeSDszP4RMJtT323zS/IfZ4vGFNBdqAAzFQnshbdKc465qwQH6ppRyOsnKej99BpDJhPp+m+u6AICb6GqSbAzx7UVmEeyBPmgnHRaTTbP3vc9keuoB2BNtL/1oz31vLjIDqLWer7UuJTnTrn5ztOuaJiXUD0DbxwUwi1bTtOQ4zwG90bbZLCf5h43jfT5XCfX9tt6bNb9xcMN/KP30wGCUUh4bM7x+0dnCXtYCMIGLSebaNex7n8mE+h5rPyWuZnMf13yauzV2/h8IYBLrN21ptxutv0E6nwGdaO+zcW2bVpu5IWQyob7/Lqa5HfFGR9txgEFo3/CWxrzxHU+zwsRWK0oATNtcmsA+en5an4RYX5u+15lMqO+/M9l8u/Wldhygj+YzfqWIqxtn6tuvrZeSPLpXhQGMqrWuJPnimIfOJDm3YTKi15ms1Fq7roEdtF8HnUrydJpPjSu11l58KgRI3g/oj6c5R51IM+N1McmFkSUsT7T73Jkm+J/tw9fWAO3a9EeSvNpuL4+u2NXnTCbUAwDAwGm/AQCAgRPqAQBg4IR6AAAYOKEeAAAGTqgHYKrauzECMEVCPQBT0y5hudVdGgG4SQ50XQAAM+1U3FwKYOrM1AOwK6WUo6WU5ybcfa7WunaDxwBgB0I9ALt1KsmmoD6qbb156kaOAcBktN8AsFvHk0xyW/SlJCdv8BgATMBMPQC7dTTJhe12KKXMJcm41ptJjwHA5IR6AHZUSjleSlkupawH8ZPt31utbPOFjLTeXMcxAJhQqbV2XQMAA1FKeSzJqVrrgzvsdyHJyS0ukp3oGABMzkw9ALuxmB164dvWm7VtWm92PAYAu2OmHoCJlVJqksVa65ahvJ2JX621nr+eY5RSTidZTnIuyfqyl0eSvFprPXcj9QPMKqEegIm0ve+Xk9yxzSx8SikXaq2L13uMdqb/Wq21jIw/luTOWuuZ6/03AMwq7TcATOp4kpX1ML6+ws1GpZSFbL/+/I7HyNbLXZ5Pcnq3RQPsB0I9AJP6bJJLG/4eF7BPpGmduZFjLGb8cpdz7Q8AI4R6AHbjctIsT5nxs+mntuu3n/AYX0gzKz9qYYv9AfY9oR6ASf1JksVSyokkqbWubHywbb25NO6JuzxGaq2rY557Km5YBTCWC2UBuClKKWeTXJhgpn67Y5xOszLOyZHxo0metLY9wHhm6gG4WY7fSKBvbeqnb2fvH09ycuwzAMiBrgsAYPjamfSdWm+2e/5cmotmTyRZbWfs3zc6cw/A/6f9BoAbVkpZTrI82iMPwN7QfgPAzTAv0AN0x0w9AAAMnJl6AAAYOKEeAAAGTqgHAICBE+oBAGDghHoAABg4oR4AAAZOqAcAgIH7P6cPBaU57MfvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a1c39de80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot jz1t, jz2t, jz12t\n",
    "\n",
    "j2max = (0.5 * N + 1) * (0.5 * N)\n",
    "jmax = 0.5 * N \n",
    "j1max = 0.5 * N1\n",
    "j2max = 0.5 * N2 \n",
    "\n",
    "label_size = 20\n",
    "plt.rc('text', usetex = True)\n",
    "plt.rc('xtick', labelsize = label_size) \n",
    "plt.rc('ytick', labelsize = label_size)\n",
    "\n",
    "fig_size = (12, 6)\n",
    "lw = 2\n",
    "\n",
    "fig1 = plt.figure(figsize = fig_size)\n",
    "plt.plot(t/td0, jzt/jmax, '-', label = r\"$\\langle J_{z,\\mathrm{ tot}}\\rangle$\", linewidth = 2*lw)\n",
    "plt.plot(t/td0, jz1t/j1max, '--', label = r\"$\\langle J_{z,1}\\rangle$\", linewidth = lw)\n",
    "plt.plot(t/td0, jz2t/j2max, '-.', label = r\"$\\langle J_{z,2}\\rangle$\", linewidth = lw)\n",
    "plt.xlabel(r'$t/t_\\text{D}$', fontsize = label_size)\n",
    "plt.ylabel(r'$\\langle J_z(t)\\rangle$', fontsize = label_size)\n",
    "plt.xticks([0, (tmax/2)/td0, tmax/td0])\n",
    "plt.legend(fontsize = label_size)\n",
    "plt.show()\n",
    "plt.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# check partial traces\n",
    "print(jz12_mat.ptrace(0)/norm2 == jz1_mat)\n",
    "print(jz12_mat.ptrace(1)/norm1 == jz2_mat)\n",
    "\n",
    "rho1pt = rho0.ptrace(0)\n",
    "rho2pt = rho0.ptrace(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2) Local-collective processes in the Dicke basis (PIQS + QuTiP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### System general and collective properties - QuTiP in the Dicke basis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "occupation number, n0 =  0.000481909935526\n"
     ]
    }
   ],
   "source": [
    "# Number of TLSs in the two ensembles\n",
    "N1 = 5\n",
    "N2 = 15\n",
    "N = N1 + N2\n",
    "# local-collective simulations with this system size take approx 5 minutes on a MacBook Pro for time integration\n",
    "\n",
    "# TLSs bare frequency\n",
    "w0 = 1\n",
    "\n",
    "# Bose-Einstein distribution determines the occupation number\n",
    "# low temperature limit\n",
    "frequency_hertz = 10**(13)\n",
    "temperature_kelvin = 10**(1)\n",
    "x = (frequency_hertz / temperature_kelvin) * (constants.hbar / constants.Boltzmann)\n",
    "n0 = 1/(np.exp(x) -1)\n",
    "print(\"occupation number, n0 = \",n0)\n",
    "\n",
    "# set collective pumping and collective emission rates (coupled ensembles) \n",
    "g0 = 1\n",
    "gCE = g0 * (1 + n0)\n",
    "gCP = g0 * n0\n",
    "\n",
    "# Local rates\n",
    "gE = 1\n",
    "gD = 1\n",
    "\n",
    "# Collective rates of the single ensembles\n",
    "gCEi = 1\n",
    "\n",
    "# Algebra in the Dicke basis\n",
    "[jx1_dicke, jy1_dicke, jz1_dicke] = jspin(N1)\n",
    "jp1_dicke  = jspin(N1,\"+\")\n",
    "jm1_dicke = jp1_dicke.dag()\n",
    "[jx2_dicke, jy2_dicke, jz2_dicke] = jspin(N2)\n",
    "jp2_dicke  = jspin(N2,\"+\")\n",
    "jm2_dicke = jp2_dicke.dag()\n",
    "# Bulding the tensor space for N1 + N2\n",
    "dim1_dicke = num_dicke_states(N1)\n",
    "dim2_dicke = num_dicke_states(N2)\n",
    "id1_dicke = qeye(dim1_dicke)\n",
    "id2_dicke = qeye(dim2_dicke)\n",
    "norm2_dicke = id2_dicke.tr()\n",
    "norm1_dicke = id1_dicke.tr()\n",
    "\n",
    "# Place operators of a single ensemble (N1 or N2) in the tensor space\n",
    "jz1_dicke_tot = tensor(jz1_dicke, id2_dicke)\n",
    "jz2_dicke_tot = tensor(id1_dicke, jz2_dicke)\n",
    "\n",
    "# Place operators of two ensemble (N1 + N2) in the tensor space\n",
    "jx12_dicke = tensor(jx1_dicke, id2_dicke) + tensor(id1_dicke, jx2_dicke)\n",
    "jy12_dicke = tensor(jy1_dicke, id2_dicke) + tensor(id1_dicke, jy2_dicke)\n",
    "jz12_dicke = tensor(jz1_dicke, id2_dicke) + tensor(id1_dicke, jz2_dicke)\n",
    "jm12_dicke = tensor(jm1_dicke, id2_dicke) + tensor(id1_dicke, jm2_dicke)\n",
    "jp12_dicke = jm12_dicke.dag()\n",
    "\n",
    "h1_dicke = w0 * jz1_dicke\n",
    "h2_dicke = w0 * jz2_dicke\n",
    "\n",
    "htot = tensor(h1_dicke, id2_dicke) + tensor(id1_dicke, h2_dicke) \n",
    "\n",
    "# Build the collective Liovillian (Hamiltonian + collective Lindbladian)\n",
    "L_collective_dicke = liouvillian(htot,[np.sqrt(gCE)*jm12_dicke, np.sqrt(gCP)*jp12_dicke])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# check algebra relations in tensor space\n",
    "print(jp12_dicke * jm12_dicke - jm12_dicke * jp12_dicke == 2*jz12_dicke)\n",
    "print(jx12_dicke * jy12_dicke - jy12_dicke * jx12_dicke == 1j*jz12_dicke)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### System local properties - Building local Lindbladians with PIQS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "case 1\n",
      "case 2\n",
      "case 3\n",
      "case 4\n"
     ]
    }
   ],
   "source": [
    "## Define Piqs objects\n",
    "\n",
    "# case 1: only collective coupled processes (already defined above)\n",
    "system1 = Dicke(N = N1)\n",
    "system2 = Dicke(N = N2)\n",
    "\n",
    "# case 2: collective coupled processes + dephasing\n",
    "system1gD = Dicke(N = N1)\n",
    "system2gD = Dicke(N = N2)\n",
    "\n",
    "system1gD.dephasing = gD\n",
    "system2gD.dephasing = gD\n",
    "\n",
    "# case 3: collective coupled processes + local emission\n",
    "system1gE = Dicke(N = N1)\n",
    "system2gE = Dicke(N = N2)\n",
    "\n",
    "system1gE.emission = gE\n",
    "system2gE.emission = gE\n",
    "\n",
    "# case 4: collective coupled processes + collective emission of single ensembles\n",
    "system1gCEi = Dicke(N = N1)\n",
    "system2gCEi = Dicke(N = N2)\n",
    "\n",
    "system1gCEi.collective_emission = gCEi\n",
    "system2gCEi.collective_emission = gCEi\n",
    "\n",
    "# Define identity operators in tensor space\n",
    "id_tls1 = to_super(qeye(dim1_dicke))\n",
    "id_tls2 = to_super(qeye(dim2_dicke))\n",
    "\n",
    "###Build the Lindbladians\n",
    "\n",
    "## case 1\n",
    "L1_local_dicke = system1.liouvillian()\n",
    "L2_local_dicke = system2.liouvillian()\n",
    "print(\"case 1\")\n",
    "# Build local Lindbladians in tensor space\n",
    "L_local_dicke = super_tensor(L1_local_dicke, id_tls2) + super_tensor(id_tls1, L2_local_dicke)\n",
    "\n",
    "# Total local-collective Liouvillian in tensor space\n",
    "L_dicke_tot = L_collective_dicke + L_local_dicke \n",
    "\n",
    "\n",
    "## case 2\n",
    "L1gD_local_dicke = system1gD.liouvillian()\n",
    "L2gD_local_dicke = system2gD.liouvillian()\n",
    "print(\"case 2\")\n",
    "# Build local Lindbladians in tensor space\n",
    "LgD_local_dicke = super_tensor(L1gD_local_dicke, id_tls2) + super_tensor(id_tls1, L2gD_local_dicke)\n",
    "\n",
    "# Total local-collective Liouvillian in tensor space\n",
    "LgD_dicke_tot = L_collective_dicke + LgD_local_dicke \n",
    "\n",
    "\n",
    "## case 3\n",
    "L1gE_local_dicke = system1gE.liouvillian()\n",
    "L2gE_local_dicke = system2gE.liouvillian()\n",
    "print(\"case 3\")\n",
    "# Build local Lindbladians in tensor space\n",
    "LgE_local_dicke = super_tensor(L1gE_local_dicke, id_tls2) + super_tensor(id_tls1, L2gE_local_dicke)\n",
    "\n",
    "# Total local-collective Liouvillian in tensor space\n",
    "LgE_dicke_tot = L_collective_dicke + LgE_local_dicke \n",
    "\n",
    "\n",
    "## case 4\n",
    "L1gCEi_local_dicke = system1gCEi.liouvillian()\n",
    "L2gCEi_local_dicke = system2gCEi.liouvillian()\n",
    "\n",
    "# Build local Lindbladians in tensor space\n",
    "LgCEi_local_dicke = super_tensor(L1gCEi_local_dicke, id_tls2) + super_tensor(id_tls1, L2gCEi_local_dicke)\n",
    "\n",
    "# Total local-collective Liouvillian in tensor space\n",
    "LgCEi_dicke_tot = L_collective_dicke + LgCEi_local_dicke \n",
    "print(\"case 4\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Initial conditions\n",
    "# set superradiant delay time for the excited ensemble (N2)\n",
    "td0 = np.log(N2)/(N2*gCE)\n",
    "tmax = 30 * td0\n",
    "nt = 1001\n",
    "t = np.linspace(0, tmax, nt)\n",
    "\n",
    "# set initial tensor state for spins (Use QuTiP's jmat() basis)\n",
    "excited1_dicke = excited(N1)\n",
    "excited2_dicke = excited(N2)\n",
    "ground1_dicke = ground(N1)\n",
    "ground2_dicke = ground(N2)\n",
    "\n",
    "sdp_dicke =  tensor(excited1_dicke, excited2_dicke)\n",
    "sdap_dicke =  tensor(ground1_dicke, excited2_dicke)\n",
    "ground12_dicke =  tensor(ground1_dicke, ground2_dicke)\n",
    "\n",
    "rho0_dicke = sdap_dicke"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "case 1\n",
      "case 2\n",
      "case 3\n",
      "case 4\n"
     ]
    }
   ],
   "source": [
    "## Solve using qutip (using the Dicke basis)\n",
    "\n",
    "# case 1\n",
    "result_0 = mesolve(L_dicke_tot, rho0_dicke, t, [], \n",
    "                   e_ops = [jz12_dicke, jz1_dicke_tot, jz2_dicke_tot], \n",
    "                   options = Options(store_states=True))\n",
    "rhot_0 = result_0.states\n",
    "jzt_0 = result_0.expect[0]\n",
    "jz1t_0 = result_0.expect[1]\n",
    "jz2t_0 = result_0.expect[2]\n",
    "print(\"case 1\")\n",
    "# case 2\n",
    "result_gD = mesolve(LgD_dicke_tot, rho0_dicke, t, [], \n",
    "                    e_ops = [jz12_dicke, jz1_dicke_tot, jz2_dicke_tot], \n",
    "                    options = Options(store_states=True))\n",
    "rhot_gD = result_gD.states\n",
    "jzt_gD = result_gD.expect[0]\n",
    "jz1t_gD = result_gD.expect[1]\n",
    "jz2t_gD = result_gD.expect[2]\n",
    "print(\"case 2\")\n",
    "# case 3\n",
    "result_gE = mesolve(LgE_dicke_tot, rho0_dicke, t, [], \n",
    "                    e_ops = [jz12_dicke, jz1_dicke_tot, jz2_dicke_tot], \n",
    "                    options = Options(store_states=True))\n",
    "rhot_gE = result_gE.states\n",
    "jzt_gE = result_gE.expect[0]\n",
    "jz1t_gE = result_gE.expect[1]\n",
    "jz2t_gE = result_gE.expect[2]\n",
    "print(\"case 3\")\n",
    "# case 4\n",
    "result_gCEi = mesolve(LgCEi_dicke_tot, rho0_dicke, t, [], \n",
    "                      e_ops = [jz12_dicke, jz1_dicke_tot, jz2_dicke_tot], \n",
    "                      options = Options(store_states=True))\n",
    "rhot_gCEi = result_gCEi.states\n",
    "jzt_gCEi = result_gCEi.expect[0]\n",
    "jz1t_gCEi = result_gCEi.expect[1]\n",
    "jz2t_gCEi = result_gCEi.expect[2]\n",
    "print(\"case 4\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualization with parameter dependence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAEsCAYAAAD0PAtkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3Xl8VNXd+PHPSQIhgYRJIGBUtgACigGToIiiVQJ1a5VNcKnLU0mgtX30pxBxKT76IAatra1WE3zcW0UC1VatkMStUtSQsMqegKACAslA2EJIzu+Pc2cySSbbZLYk3/frdV935s6995xB5Dvn3HO+R2mtEUIIIYT/hQS6AkIIIURHJUFYCCGECBAJwkIIIUSASBAWQgghAkSCsBBCCBEgEoSFEEKIAJEgLIQPKaWmKKV0c7dA17etUkplNufPr7nnCeEvYYGugBDtXB4wvs6xdGCKtS/xe42EEEFDgrAQPqS1tmMCsZNSKgkThPO01hKEhejApDtaCCGECBAJwkIIIUSASBAWIsgopbKUUmV1jmVag7cSXI4lWcdSXd7nWscKlVJpLSx3inWd4/rUOuWXuSkjyc09il0Gm+XWqXNzy6h1jlIqwaXcYqXUlAa+g6N+ZdZ5qe7Oa8H3bvT7CNFaEoSFCD65gK1OgEsF7NbeYRqA1jrPCkqFmIFe44HFQJZSKqs5BVrXL8E8v55q3aduwLFZ52RhBpUlWO8d93C8zwGSXe5j86CMBXXKcGxTrfOW4N4iq34Z1vtGg2ZjdWrq+wjhFVpr2WSTzY8bMAfQQEIj52gg03pts95nAbku5xQ63gNlQFade6Ra16U2o05ljvLq3N9Rh0zrXlNcPnccs1nvp7i+b0UZqW7KmONyzFFOUt3z6tw7wfHn1tB5jdWpqe8jm2ze2KQlLERwKqKm1ZuK1UJzHFNK2YAkYInVfWrDBGknrXUepvU8lUZYLW4bMKfOnOUka3PlOtK7uIHP8pVSc1xb8i0sY7WbMlzLdYwoj23se2kz8jwPSHH3eTPq1OD3EcJbZIqSEMFpMaY1BqZ7OQ8TmB3BwxGA8qgJ1u6mO5VgWoSNcXyejAnarkpd32gz5cotrbVdKTXeqnemVdc8rfV4L5TRYLlNKKGBINxUnZr4PkJ4hbSEhQhOOQBWKzcF0+1cQk1QGQ+UuBwD98E2gaYTgjg/11qX1NlaFPy01nla62QgBvNcNlUpNcebZbRQY9+/yTo18n2E8AoJwkIEISu42qnfNerIwOV6zNF9m+56D2vQkWOgU2NlFVllZdb9zOr29qT+dq31QkzrfZQvymiKdd9U6iRLcaljs+tU9/t4v7aio5LuaCGC1zuYwFrk0lp0jBK2YaXDtLpN0zGjoWMxXdkJmOCSYz0bbsoMzPNlx+hnGzVpNdMbu9DBmhKVYV3v6AZPouZZdavLaEYdllh1cHx/O2akdUMarJNSqrCJ7yNEq0kQFiJ45QKOwAY4pyPZHK9djmcrpUowgWcJJmhkWK23Jmmtc5RSydb1uZjg9Y5r2c3wDjCQmqlFJcBCrXW2F8toymLr3glY046aeI7dVJ0a/D5CeIPSWhYUEUIIIQJBngkLIYQQASJBWAghhAgQCcJCCCFEgEgQFkIIIQJEgrAQQggRIDJFqZV69uyp+/fvH+hqCCGECCKFhYUHtdZxTZ0nQbiV+vfvz+rVq5s+UQghRIehlPq2OedJd7QQQggRIO2yJWzlzLU3M12fY1WaFGpS05U091ohhBDCU+0uCFurziyiiTVUXc5PwCzqPd7l2BKllGOFGiGEEMIn2k13tFIqQSmVhWnJljZ1vot06idkz8LNyipCCCGEN7WbIGytAZruQXL1KViLpbtYbR0XQgghfKbdBGFPWKvR1Gs5O1ZdsbqqhRBCCJ/o0EEYiIWaoOuG74Pw6dOQnw+Z0vsthBAdTbsbmNVCtkBXgMpKuPZaqKgw++HDA10jIYQQftLRW8IeUUqlKaVWK6VWHzhwoHU3i4iAyEjz+uWXW185IYQQbYYEYZzPhptNa52ttU7RWqfExTWZlaxp555r9itWtP5eQggh2oyOHoQd84BjXQ+6BGX/zBO+4gqz37wZysv9UqQQQojA69BB2BqQVUL9Z8OxmIxb/gnCyclmX10Nn3/ulyKFEEIEXocOwpY8TMpKV0nWcf8YOrTm9ccf+61YIYQQgdVeg3AsbkY+K6VsSqlipVSay+EM6qe4TLeO+8fAgRAaal7n5/utWCGEEIHVbqYoWc9x52Lm9tqATKXUeCBXa53jcmqt579aa7tSKkMplQkUWNdn+jVvdKdOMHgwbNkC69ZBaSnExjZ9nRBCiDat3QRh6/luo61X65wYN8eLqJ+60r/Wr4cJE+DTT+Gzz2DixIBWRwghhO+11+7otqdTp5pR0p98Eti6CCGE8AsJwsHkkkvM/quvAlsPIYQQfiFBOFhs2QK3325er1kDJ08Gtj5CCCF8ToJwsDjjDPj+e1DK5JNeuzbQNRJCCOFjEoSDhc1mArHW5v2XXwa2PkIIIXxOgnAwcU3aIc+FhRCi3ZMgHEwGDap5LS1hIYRo9yQIB5PBg80+LAx27TJJO4QQQrRbEoSDiaMl3K2b2a9bF7i6CCGE8DkJwsHkwgvhmWdg7FjzXkZICyFEuyZBOJicfTbcey9cc415L0FYCCHaNQnCwWjkSLOXICyEEO2aBOFg8/nnsHy5eb15M5w6Fdj6CCGE8BkJwsHmjTfg0UehVy+TOWvTpkDXSAghhI9IEA42jmlK3bubvXRJCyFEuyVBONg4gnCI9Z9GpikJIUS7JUE42DjmCh89avbSEhZCiHZLgnCwGTjQ7PfvN3t5JiyEEO2WBOFgExlp5gufPm1e//gjHDwY6FoJIYTwAQnCwWjoUEhIgP79zfvNmwNaHSGEEL4hQTgYrVgBxcUmjSVIl7QQQrRTEoSDkVJmf+65Zv/NN4GrixBCCJ/xTRBWahJK3eWTe3ckjulK0hIWQoh2Kcxrd1JqJDATSAYWAwqlVgMFQBZay1yb5iopgTFjoGtX816CsBBCtEutC8JKRQPpwDSgGMhG65kuZzyFUhcAM1HKEZyz0fpIq8pt7844w0xR6tQJIiJg714oK4OYmEDXLKhorTl58iSnTp3i9OnTtbaqqqp6x06fPk11dTVa63qb436tOS5qkz8T0VZddtllREdH+6Usz4KwUpMxwTcG08pNafBcrddgWsig1AzgY5Q6ZF23zKPy27vISDjrLPj+exgyBDZuNCOkx4wJdM184siRI+zevZvdu3ezZ88eDh06RGlpqXNfWlpKeXk5J06c4Pjx47X2QgjhbWvXrmXEiBF+KavlQVipF4EyIB2td7boWq0XAYtQagCQjlLT0Hpai+vQEQwaZIJw794mCG/a1OaD8MmTJykoKGDt2rWsX7+e9evXs3XrVg4fPuzxPTt37kyXLl0ICwsjNDSUsLCwBrfQ0FBCQ0MBUErV27xxXNQmfyaiLYqKivJbWS0PwrW7mz1jgvcDrb5PezZoEHz2GXTrZt630efCmzZtYtmyZeTn57Nq1SoqKirqnRMREUG/fv3o27cvffr0IS4ujtjY2FpbVFQUkZGRREZGEhERQWRkJF26dHEGVSGEaIu8NzBLeJcjh3R1tdlv2xa4urRQWVkZL7/8Mq+88grf1JlelZiYyKhRo0hMTCQxMZFzzz2XuLg4aTEJITok7wdhM1grAbADpTIIy0N1F3JoA0F47969zJ8/n1deeYXjx48DEBMTww033MDPfvYzLrvsMnr06BHgWgohRPBofRA2g62mAimABkowz4wBElAqxjq2GMhB612tLrMjuOgi+NOfYNgw+PRTM23p1Cno3DnQNavn5MmTLFiwgKefftoZfFNTU7n77ru55ppr6NSpU4BrKIQQwcnT0dHRwIPABcASYCpaNz66RqlxwANWUM5C6489Kruj6NMHfvMb87p/f9i502xDhgS0WnUVFBRw++23s9nKb33DDTfw+OOPM3z48ADXTAghgp8no6MvANKAhS0aHa11PpBv3WMGSiWh9dMtLr8jOuccE4C3bQuaIKy15tlnn+X++++nqqqKIUOGsGjRIsaOHRvoqgkhRJvRsrSVZmpRAlrPavH0JFc1U5UmeXyPjiA/H+bPh549zfutWwNbH8vp06eZMWMG9957L1VVVdx7772sWbNGArAQQrRQy1rCJvB6Hnxr3+sw4PVkHUqpJMzz6RLMALESrXVeE9fYMK37HKAUiMUkI8lt6lqfev11s918s3kfBIOzKisrufnmm8nJySEiIoJXX32VG2+8MdDVEkKINsm7uaMDnB9aKZUAZGqtx7scW6KUKtFalzRyaSyQaW1gRnbPCGgAhpoR0lVVZh/gIFxVVcVNN93E0qVL6d69Ox999BGjR48OaJ2EEKIt8+YqSun1jijVHaVmW4s7+EM6kFXnWBY1wbUx4zFpOAdqrWO01jnerlyLOYLwEWuWV4C7o2fPns3SpUux2Wzk5eVJABZCiFbyPAgrtQOlClDqCZS60u05Wh9G66eAHg2e411TgKI6x1Zbx5uktbY30WL2L0cQ/uEH6NIF9u2rCch+lpWVxR/+8Ac6derEu+++S0pKw+nChRBCNE9rWsJTMaOdJwB5QJoVlBeg1JXWNCbDjIy2taqmTbCe6yZgnum6FK3t1ucJvizfJxxBeMeOmtfbt/u9GuvXr+e3v/0tAIsWLeLyyy/3ex2EEKI98vyZsFkdaY3zvVIrgFxMt24GoFGqCNMSLbSO+3LVpFhTLRN03UjADNZqSIJSytFijgVKA94lHRMDsbFQWgr9+pmFHLZuheRkv1Xh+PHjTJ8+nVOnTpGWlsbtt9/ut7KFEKK982bayhKr6/kpAMwo5XHAKExreYEXy3KnNS3tUgDXoGsN6CLggXjYMPjxR4iPN+/9PDjrkUceYfPmzQwbNow//OEPfi1bCCHaO28G4SW13mldRP3ns0HJaj1n1zmcZW31grBSKg0zpYm+ffv6tnL//jcoBa+9Bi+95NcgvG7dOp599llCQkJ44403iIyM9FvZQgjREXhvdLR57htw1rNhbyjBdFHXu5/WOltrnaK1TomLi/NScQ1wrC50zjlm76cgXF1dzaxZs6iqquLuu+8m2Y9d4EII0VG0PGOWt7JcKRWNUnd55V6G43lvbO1inEG0wefBSqk5bg47BngFx4Cufv3MfutW0Nrnxb311lusWrWK+Ph4Hn/8cZ+XJ4QQHVHLgrAjY5ZSL6BUf49LNSsvzUXrlzy+Rx1Wl3IJ9Z8NxwINTj1yJPhwM3raEcwDO2VpwwaIi4Of/9wM0jp61ExV8qHKykp+97vfAfC///u/REdHN3GFEEIIT7S8O9qMin4AmIlSy1HqrlrTkRqi1EiUehGllgPFaD23xWU3LQ+TstJVknXcLSs4p7sJ0qlAUSOjrf2jd284eNBMUxo82Bzz8TSll19+mZKSEoYMGcJtt93m07KEEKIj8+yZsEnC8QBa/xSzdvBLVvKO7dZc4eXWtgOlDqFUAY5sVlr/1IfLGGZg5i+7SreOA6Z7WilVbA2ucih1bQlbXdjpwAwf1bP54uIgKgoOH67pkvbhc+GKigoee+wxAB5//HHCwrw5dk8IIYSr1v8Lq/VSYKnzvVLdMc9R7UBpk+sMe5HW2q6UylBKZQIFVj0y3bRyY+tcl6OUmmLNE+6B6dKeGhTZs5QyiTrWrIHu3c0xH7aE//a3v/HDDz9w/vnnM3nyZJ+VI4QQwrtTlAwTdNc0eZ6P6CamRlndyzFujgc+V3RDHEE4PNy891EQ1lrzzDPPAHD//fcTEuLN1OJCCCHqkn9l2wJHysrKSrP3UXd0bm4uGzduJD4+nunTp/ukDCGEEDUkCLcFjiBst8aI7dgB1dVeL8aREevuu++mc+fOXr+/EEKI2prujlZqAGaQUhLgOkHVDixGa1/mgxYAl10Gzz9vckZ/9pmZorRnT81ALS/49ttvWb58OeHh4aSn11+VUgghhPc1HITNAKu5wA7MqOadbs4Zh1IvAiskGPvQoEE1reHBg00Q3r7dq0H41VdfRWvNxIkT6dGjh9fuK4QQomHuu6NNAJ5qTUN6yW0ABpOqUuuZwBqvZdISjfNB+srq6mpeeeUVAH75y1967b5CCCEa574lbEY4Nz+blSOTlvCdDz+EggKTNQu8OkI6Pz+fb7/9lv79+3PllVd67b5CCCEa17yBWc3JiCV866WX4NFHoarKvPdiEH7ttdcAuPPOO2VakhBC+FHj/+KaBRt2AHYr89VEl88mo5Rv8yeKGj6apnTixAnee+89AG6++Wav3FMIIUTzNNXseRIzMjoGmAbc5Hz2azJlDfRp7UQNRxAuKzP7nTtrAnIrLF++nKNHj5KUlMQgRxlCCCH8oqkgvNoafHUYrfPQ+kagB0o5Hhz6fk09YTgC5K5d0LcvnD5tXrfSO++8A8C0adNafS8hhBAt09Q84forCGm9yJqadIFvqiTccgThHTvgvPNg927zXNixspIHTpw4wT/+8Q8Apk6tu+6FEKKiooLS0lLKy8upcozHEB1OaGgoUVFRxMbGEu5IH+wlTQXh1dbav3OA8Wi9CzBTk0wQ9tviDB3e2Web3NH79sF115ljrRyc9dFHH3Hs2DFGjRrFgAEDvFBJIdqPiooKdu/eTUxMDP3796dTp04opQJdLeFnWmsqKys5cuQIu3fvpm/fvl4NxI0HYa3XoFQJplt6l5vPkhu9Xqn+9a4TngkJgWHDoKLCrDEMrR6c9c9//hOAiRMnNnGmEB1PaWkpMTEx9OzZM9BVEQGklKJz587OvwelpaXEx8d77f4NJevo73xtnge7XxXJNYmH6zU1ypzTm5TqjlJ3yXSnVigqgk2bYPRo874VLeHq6mo+/PBDAK5ztKyFEE7l5eVER8s/V6JGdHQ05eXlXr2n+yCs9S6Umt1AYK3PdFnb3HwSA4xCqWgrAcgSaxOecHSFeSFrVlFREfv376dv374MHz7cC5UTon2pqqqiU6dOga6GCCKdOnXy+tiAhrujtX4KpWagVBJQjFmj17HIvQ1IACZgRkhnoPURN3fJRevBVivY8ZNSEhO3htZwxhkQGmoGZ508CV26tPg277//PgDXXnutPOcSogHy/4Zw5Yu/D41PUdJ6EVrPAvKBZGAm8ABmzjCY4DurgQAMkI1SI61WsLLuMc4rNe+IVq0yaSuvvx4GDDABubjYo1t98MEHgAnCQggRDHJyckhOTkYpRUxMDFOnTqWkpKTpCxuRkZER1D+mmpejUOs1aP2UtaDDTCAbrZdZwbWx657CkVPanLsTmVvsuV69zJrCO3bUTE3y4Lnw/v37Wb16NREREZIrWggRFNLT05k6dSoJCQlkZWUxd+5cSkpKGDhwIDk5OYGuns80vZ6we09S0xpunGug1nqntUKT8ES/fhAWBt99Bz//uTnmQRDOz88H4LLLLiMiIsKbNRRCiBbLyMggOzubJUuWMGXKFOfxOXPmOINzYWEhSUlJAaylb3iarX8qSr2NUhNbPNq5qdazaFhYmOmGBuhu/ZbxYHDWxx9/DMC4cfJkQAgRWHa7nYULF5KWllYrADtkZWWRkJDAjBkzAlA732vNkjkTgKWYaUjbUeoFj4KyaJlhw8ze8YyjFS1hCcJCiEBbsGABAJmZmQ2ek5mZSVFREUVFRf6qlt94GoRz0DoWGATMAtZguqfrB2XhXUOHmv3x42bfwpbwzp072bVrFzExMYwYMcLLlRNCtDfZ2dnExMQQExNDRkZGvc9b+7y2qKgIm82GzeZulqvh6IZevXq181hGRgYxMTEUFRUxfvx4lFIkJyc3Gqgdg7Ts9toZmRcuXOj2uD94GoRLAdC6BK2z0fpGKygPpG5QNksgXuGV2oqalvDevdC5s9kfPdrsyx2t4CuuuILQ0FBf1FAI0U4sXLiQJUuWkJmZydy5c8nJyakVdLOzs912IbfE6tWrSUhIaPQcx+eFhYW1jtvtdqZOnUp6ejpZWVmUlJQ0mgc/PT3dWW9XixcvJjU1tdEfAr7iWRA2I6TdHd/pJii/BOShVD+PaylqXH45vPAC3H9/7UUdmskRhGVUtBCiMSUlJRw6dIjc3FzS0tKYM2cOhYWFLF682HmOa8sxLy/P7WtvaKyFmpmZyZQpU0hLSyMtLY2SkpIGz09ISCApKYmsrKxa9y4qKgrYIjateSbcNBOUM4BRwEKfltVRDBgAM2dCSkrNNKVmdklrrfnkk08ACcJCiMYVFRXVe05rs9mYNm0adrudnJycWq3gJUtqkiHm5uY2u5yUlJQmn/U65gonJ9dfriA1NdX5euDAppe4d0x9cpTpaBWnpaU1u87e5Nsg7KB1EbC6yfNEy7RwrnBxcTH79+8nLi6OoY5ny0KIFlNKtYmtNRrqZk5NTSUvL4/S0tImu5Gbw/G8t7FA7GhZuyuvpV3IU6ZMwWazOVv0ixcvbnWXemv4NggrNQmlRlrvJEmHt+Tnw8MPm6UNodkt4ZUrVwJwySWXBHUGGSFE8LLZbOTm5nolAEPNc1rHKGl3FixYQFJSUq1Wb2vceOONZGdnO7uiHXUIBF+3hBcC+Sg1CVl72Hv+/neYPx9Kzfi45raEXYOwEMJzWus2sflKXl6e1wJiQkICc+bMIScnp96AKYCpU6dit9tZtGiRV8oDM0rabreTnZ2NzWbz2nfxhKcZs5orGbgRKEHrZT4uq+NwdCUftn7XeNASFkIIT3k7aGVmZmK320lPTyc3N5fx48djt9tZvHgxRUVFLFmyxKvZshwDtBYsWMCNN97otft6wtcDsw5bi0Cs9Wk5HY1jmtK330LXrnDoEJSVNXpJaWkpmzZtIjw8vF2mfhNC+IdjXq63ZWVlsWTJEkpKSkhPT2fBggXExsZSXFzsk2e26enpzsAfSP4ZmCW8yxGEt2ypmabURJf0qlWrABg1ahThjmfJQgjRQnW7ohtb5ailKyBNmTKFwsJCtNaUlZU1+uw5MzOzXpd7WloaWutag7XcnQdmapKjRRxILQ/CSs1GqdbPb1FqJErd3+r7dETx8RAVZVrA/azp1010SUtXtBDCG4qLi2sFOUfLtdQxRgUT4DIyMoI6zWRWVlbAW8Hg2TPhHCADpbKs11lovatZV5q80umYbFoFQP2n8KJpSpnW8Ndfm2AMTbaEJQgLIbzBNdiCeT6ckpJSKxf96tWryc/PD0gGqsY45gcXFBRQUlISsLnBrloehLXeCZiMWUpNBrJRSmOCsfvBV+a8dMw0pSVoneJhfYVDUhKcOGFaxdBoS7iyspKvv/4agDFjxvijdkKIdshutzNtWv1VbG02G/n5+SxcuJCEhASWLFkSdAEYTP0dqzEFSx1bNzpa66WY/NADgHSUWgjkAlmAHROsxwF5QLoVwIU3/OUvpkX8n//A00832hLeuHEjJ0+eZPDgwfTo0cOPlRRCtCc2m63BQVKOgFZQUBDQ5BeNSUpKoqyJQaz+5p0pSia4PgA8gFLjgAcxrd5stH7AK2U0k1IqCUgBSoAEoERr3WQiU0+vCxhHsg3XrFla1xx3UVBQAJhBWUII4StpaWn1uqtF47w/T1jrfCDf6/dtBqVUApCptR7vcmyJUqpEa93gMD1Prws4reHkSeje3cwZ/vFH6N273mkShIUQ/tDUkoSivvY2RSkd0xXuKgtoeLXo1l0XWOedB3371oyQbqBL2rEGZ0qKPIoXQohg0t6C8BSg7pj41dZxX1wXWGefbfaOX55uBmedOHGCDRs2EBISwgUXXODHygkhhGhKuwnCSikb5llurQcSWmu79bnbGd+eXhcUEhPNPjTU7N20hNetW0dVVRXnnnsuXbt29WPlhBBCNKXdBGEgFmqCpxsNBVNPrwu88883++PHzd5NS1ieBwshRPBqT0HY09EALb5OKZWmlFqtlFp94MABD4v1AkdLeP9+s3fTEpbnwUIIEbzaUxD2G611ttY6RWudEhcXF7iKDBtmuqJ37zbvd+yA6upap0hLWAghgle7C8LWM16/XRdQXbrAOeeYwBsTYzJoff+98+Py8nK2bNlCp06dSHS0moUQQgSN9hSEHfN5Y10PugTXhub7enpdcHj2WfjiCxgyxLx36ZJet24dWmuGDx8uKycJIUQQajdB2BpYVUL9Z7yxgL2hpBueXhc0xo+HSy6BoUPNe5fBWRs3bgSQVrAQQgSpdhOELXmY1JOukqzjvrgueJxzjtm7tIQdQXj48OGBqJEQQogmtLcgnAFMrXMs3ToOmG5mpVSxUiqtJdcFrepq+N3vYOlS894lCG/YsAGA8x1TmYQQQgQV7+eODiCttV0plaGUysSsV+zICV23SznWw+uCT0gI/PWvUGJVdetWALTW0hIWQogg55eWsFJqpFLqLqXUSF+XpbUu0lpnaK1ztNYL666EpLW2a61jtNbZLbkuqDnmACsFxcVQUcG+ffsoLS3FZrNx5plnBrZ+Qog2Kzs7m5iYGGJiYsjIqN85mJOTE4BaudeW6urgr5bwg5g8zCVKKQ1kaK2X+ans9i85Gd55B6KjzWpKW7ey0Urgcf7556PcLG8ohBBNWbhwIbm5uWRmZmK328nKymLUqFHO9YKzs7NJS0tr4i7uTZ06Fbu9oUSF9WVlZZGQ0HACQ1/W1Zf8FYQzgAVa6zUASqnufiq3Y0hONvsQq2Pjm2/YsHcvIF3RQgjPlJSUcOjQIXJzc53H0tLSmDFjhjOw1Q2iJSUljQZKV0uWLAlYXfPy8khNTa33OhD80h2ttd7pCMDW+8P+KLfDSEoy+/Jys//mG3keLIRolaKiIjIza6/marPZmDZtGna7nZycHGeAc8jKqrsirH+0tK6uPwBcA3cgtLfR0R1TTAwMGgSnT5v3LkFYRkYL4RtKqQa37OyaISfZ2dmNnusqOTm5wfNcu1ILCwsbvWdhYWGrv1/dAOuQmppKXl4epaWlzW71+lpbqmtd7Wp0dIc2aRJs2gTvv4/+5hu+sdJXnnfeeQGumBCiPbHZbOTm5jJ1at1ZnS3j7WfC7nirrr7kkyCslJoExGqtX/LF/YW0d5YgAAAgAElEQVQbmZlQUQFdu0JxMVXV1Zx55pnExsY2fa0QosW01s06Ly0trdkDgprbgk1OTm52+b6Ql5fX6q5nbz4Tbow36upLXuuOtqYhvaiUcsyzjbGW+3vBH1OTBBAeDoMGoaqrGYp0RQshfCOQA5laKtjr2qqWsFIqGpNZahpQDGRrrWe6nPKUUuoCYKZSKhlYbJ1zpDXligbs32+eDwPnAfEyKEsI4WVFRUWMHz8+0NVolrZQV49awkqpyUqpFUA+UGatrTtNa51f91yt9Rqt9Uyt9SjgMPCxUmq51WUtvOnWW+HLLwEThGVktBDC2wI9pacl6ta1pKThJIiNfeZLLQ7CSqkXMYsdpGutR7Xkua/WepHWOgWYCVyolFrc0vJFI8aMcb6UICyE8IXi4mJstrax/HrdupaUlJCenk5paanzmN1uJyMjg6KiokBUseXd0XW6mz2itd4JPNDa+4g66gThM889N3B1EUK0S64BDMwUrOLiYuc8XbvdzoIFC2plqwqUunVNTU0lJSWFcePGOY+tXr2a/Pz8gP2wkClK7cno0eiQEKiuJgGZBC6E8C673c60adNqHUtLS8NutztbmAsWLGDu3LkBby27qyuYaUv5+fksXLiQhIQElixZEtC6ej0IW4O1EgA7UCqDsPyoe3dK+/enR0kJCmDLlppsWkII0Uo2m81t69Zms5GRkcHAgQMpKysLeAB21KmhlrijfgUFBQFvrbc6CCulZmDW4k0BNFAClFkfJyilYqxji4EcrfWu1pYpGraxVy8udwww2LBBgrAQwi8SEhLIysoKigDcHGlpafW6qwPBoyBstXYfBC4AlgBTm8oHrZQaBzxgBeUsrfXHnpQtGpevNWOxuqLXrQtwbYQQHUkwrlLUEJvNFhQ/GFochK15v2nAQmuAVbNY05fyrXvMUEolaa2fbmn5onHv7N/PVky3A2vXBrg2QgghGtOiIKyUGgAkaK1ntaZQrfUipVR3pdQkWVfYe8rLy9m6axcVnTpBZaUJwlqDrCcshBBBqUUDaK0lCZd6o2Ct9WEJwN61adMmAKKHDoUePaCsDPbsCXCthBBCNMSruaO9dS/hGcfyheP69YOjR83BNWsauUIIIUQgeXMqaXrdA1aX82wJ0P7hCMK9R4+u6YLOywtgjYRo2wK5UpEIPr74++BxEFZK7VBKFSilnlBKXenuHKvL+SmgR0PnCO/ZsGEDAOeNGAEXXGAOfvpp4CokRBsWGhpKZWVloKshgkhlZSWhoaFevWdrWsJTMaOdJwB5QJoVlBcopa60pjEBzpHRgR8L3s45WsLDhw+HG24wB3fsCGCNhGi7oqKiOHJEcg2JGkeOHCEqKsqr9/Q4CFurIz1graAUggnI7wDJmKBcZgXlF5RSd2GWOxQ+cuDAAfbv30+3bt3o168f3HGH+eDkSfj220avLSoq4oEHHmDs2LH07t2bbt260atXLx588EHfV1yIIBUbG0tZWRkHDx7k1KlT0jXdQWmtOXXqFAcPHqSsrIzY2Fiv3t+baStLrK7npwCUUknAOGAUprW8wItliTpcW8FKKejVC7p1MwO0Xn0V5s2rd43WmhtuuIF//OMf9T47duwYVVVVzvdbtmxh9uzZzJo1i6uvvtqUIUQ7Fh4eTt++fSktLWXXrl21/n8QHUtoaChRUVH07duX8PBwr97bm0F4iesbrXUREJi1oTqgWl3RDueeC19/Df/6l9sgrJRiyJAhREdH84tf/ILrrruO4cOH0717d44ePVrr2cebb77J+++/z/vvv09iYiIPPvggU6dOJSRElokQ7Vd4eDjx8fHEx8cHuiqinfLav6DWc18RII4gfP7559ccdCQm79fPeaiqqoo9LnOHH3zwQXbu3Mlzzz3HVVddxdlnn01UVBTx8fH06tXLed5vf/tbFi5cSHx8POvXr2f69OlcfPHF/Oc///HtFxNCiHasRUFYKTVAKTXJGwUrpaKtZ8XCCxwjo2u1hCdMMPvVqwHT/Txz5kwuvvhivrWeE9tstmY94+jVqxezZ89m586dZGVlER8fz9dff80ll1zCY4895t0vI4QQHUSLM2YBO63BVv09LdRaeWmu1volT+8hamit3XdHn3ceREZCSQl8+y2ZmZm89NJLHDp0iO+//96jssLDw0lLS2Pbtm088sgjREdHc/3113vjawghRIejPBnxp5TqDsylZhWld5paN9hK2DETGABktpdVlFJSUvRqq6UZKLt376Zfv3706tWL/fv31/7w0kth5UpOh4fTo6KCcqXIyclh0iSvdGhQXl5ea8j+E088wS233GJGaAshRAellCrUWqc0dZ5Hz4StJBwPaK1/ilk7+CUrecd2a1rScmvboZQ6pJQqwGTUytJa/7S9BOBg4bYr2uHiizkNhFVUMBnzDNhbARioFYCXLVvGQw89RGJiIn/961+9VoYQQrRXrR4dbS3o4FzUwWolJwB2oLSpdYZF67ntina46CLeBm4F7unalWFuRkl7y6WXXsr111/Pe++9x6233soHH3zA888/T0xMjM/KFEKItszr80usVvIaa8UlCcB+4HZktKW4d29mAUeAxGPH6LRtm8/q0atXL/7+97/z0ksv0bVrV9566y1GjBjBJ5984rMyhRCiLZNJnu1AY93RAy+9lLdsNkocB17y7Vg4pRS//OUvWbt2LRdddBF79uxh3LhxLFsmq1YKIURdEoTbuNOnT7N582YAzj333PonKMV1l1+Ocxmr11+HEyd8Xq9BgwbxxRdfMG/ePIYNG8ZPf/pTn5cphBBtTZNB2Job/KRSaoXLgKvlSqnF3pozLDy3Y8cOTp06Rb9+/YiOjq71mTMpx0UXmX1cHBw5AqtW+aVuYWFhPProoxQVFdG1a1cAjh49yhtvvCF5eIUQgkaCsLUW8JOY/M9ZWusJ1shmxzYNOKyUelGCceA09Dz4008/ZcCAAdxzzz1wySXmYPfusGsXXOnfVSVdc63ec8893HbbbUyaNImDBw/6tR5CCBFs3AZha4TzVGsa0ktWko56tNb5WuuZwJpAB2KlVJJSKk0plerYN+Mam1JqjlIqwXqdoJTKbM61waKh58FPPfUUVVVVZmTyhRdCeLhZ1jAiIhDVdEpNTaV79+68++67JCYmkpubG9D6CCFEILkNwtYI52aP4LFGQgds5I1SKgGTACRba52ntc4G0q3jjYkFMoFizHznQqBAa53n2xp7j7vpSTt27OBf//oX4eHh/PrXv4YuXWq6pL/4AqqqzKIOAegSnj59OuvWrePSSy9l7969TJgwgfvuu4+Kigq/10UIIQKtvQzMSgey6hzLwgTYpowHYoCBWusYrXWOtyvnS+66o1944QW01tx000307NnTHLzsMrP/7DMYNw6uuQZWrPB3dQHo168fn376KY8//jihoaE888wzXHTRRZSXlwekPkIIESitCsJKKf8+XGzYFOovm7jaOt4krbVda13S9JnB5cSJE+zYsYPQ0FCGDBkCmHWAX375ZQDuvvvumpMdQfjf/4ZrrzWvH344IK1hMOtzPvzww6xcuZKEhARGjBhRK/uWEEJ0BK1tCY/3Si1aQSllw2ToKnU9rrW2W5831SXdZm3evJnq6mrOOecc5+CnxYsXY7fbGT16NMnJyTUnjxkDYWFQVAS33QZnnGFWV3r33QDV3rjoootYu3Ytzz33nPPY1q1bOXDgQABrJYQQ/tFoELamJVU1sFUDc/xUz8bEQk3QdaOpIJyglJpibWlKqSZbz9Z5q5VSqwMZLNx1RW/fvp2QkBDuuqvOKpFdu0JysnkevHataQUDzJ4NJ0/6q8puRUVFOVvBx48fZ+LEiSQmJrJ8+fKA1ksIIXyt0SCstZ4AzNRah7rZQoAH/FPNRtlacW0pgNY6x9qygWlNBWJrAFiK1jolLi6uFcW3jruR0QsWLOC7775j+vTp9S+44gqzz82FtDQ491woLobf/94f1W2W8vJy4uLi2LdvH1dddRX33nsvJwP8I0EIIXylOd3RjT0rbTOjiN2xngVn1znc3AFdAdfQwg3x8fHO5Bi1TJhg9itWQKdO8Oc/m/fz58MPP/iyqs3Wu3dvPv74Y+bPn09YWBh//OMfSU5O5quvvgp01YQQwuuaDMJa6/xGPlvT2LVKqf4tqYw1T7e5m63Ota1pEbsqwXRRe+t+PuPaHa215j//+Q/V1dUNXzBmDERGwoYNJuheeSXceSc89ZR5RhwkQkNDefDBB1m5ciXnnHMOmzZtYsyYMTz22GOBrpoQQniV26UMlVL9tda7WnKjBq4pU0pFa62POBKAAO9orY+4uT6BlrVAC4CF1LTUYzHLJzru5wiiDbbklVJztNYL6xx2DPBKoP6I66BRVlbGd999R0REBAMGDKCwsJBLLrmE0aNHs6qhtJTh4aZL+oMPTGv4jjvAGkkdjC688ELWrl3Lo48+ytNPP82AAQMCXSUhhPAqt0FYa71LKTUbWNKcYKyUmoEJinXFAAOVUgVa68NKqSXAEqBeNn9ritDUllTeus6ulCqh/rPhWKDBqUeOoK+UyqlzTqy1D+opS9988w1gFm0IDQ0lJ8dMb05JSWn8wp/+1ATh5ctNEHa1c6cZuDVokA9q7JmIiAgyMzO58847ndOwAFasWMHFF18s05qEEG2a2yAMoLV+Sik1QymVhMkoVURNYHJMC5oAaCDDXesWyNVaD7byUDtWF+jhveo75QEp1G65JtHIM2utdYlSKt1NkE4FihoZbR0UXAdlaa1ZunQpAJMnT278Qsdz4dxcqK6GEOuJxKefws9+Bn37wpdfQpAFt6FDhzpfb926lZ///Of07NmTP/7xj0yePBmlVABrJ4QQnmlqdPQirfUsIB9IBmZiRkRPs07J0FrPaiAAA2QrpUZqrQ8DyrrHOO9UvZYM6rei063jgDNPdLFSKs3lnFLXecRWF3Y6MMMHdfQqRxBOTExk/fr17Nixg7i4OMaOHdv4heecA/37w6FD8PXXNceTk6FPH9i0yTwnDuJVjrTWJCYm8v333zN16lSuvvpqtm/fHuhqCSFEizUrWYfWeo3W+ilrQYeZQLbWepkVXBu77ilgp/X6sPXa6/+6W63WDGvxhSlKqTmYXNJ1W7mxda7LAZKsRRwyMc+kp2qtg/ZZsMP69esBE4QdreCJEycSGhra+IVKwc9/bl7//e81x6OiTOKO6GhYuhSCeBDU0KFDWbVqFS+88AI2m43ly5czfPhw5s2bxwk/rJUshBBeo7Vu8QYs9uQ669runl4bjFtycrL2t+rqah0dHa0BvX//fj1s2DAN6OXLlzfvBp9+qjVoPWiQ1tXVtT/75z+1Dgkxn//lL96vvJft379f33HHHRrz407feeedga6SEEJoYLVuRgzxNG3lVKXU20qpiS7Pepsb9BttPYum7d69myNHjtCrVy+UUpSVlRETE8MVjmQcTbn0UujZ0yxtaA3wcrruOnjxRfP617+Gt9/2buW9rFevXrzyyit8/vnnJCcn88ADNfljKisrA1gzIYRoWmtyR08AlmKmIW1XSr3gSVAWLefaFR0XF8d3333H119/TadOnZp3g9DQmi7pZW5WoJwxwyTwCA0NugFaDRk7diwFBQWcc845gOnhSU1N5Ze//CV79+4NcO2EEMI9T4NwjtY6FhgEzALWYAZr1QvKXqqncOEahMEktxjU0mlFkyaZvbsgDDB3rskx7VhxqQ1wHSG9ceNGVq1axcsvv8zgwYN57LHHZKlEIUTQ8TQIO3Iul2iTR/lGKygPpE5QVkodUko1s59UNIcjCA8dOhS73cOZVOPGQffusG5d/S5pMAO4zjuv5v3y5XDvvXD6tGfl+dn555/PN998w8SJEzl27Bjz5s0jISGBZ555RgZvCSGChkdBWJsR0u6O73QTlF8C8pRS/VpRT+HCEYQPHTpEz549+X//7/+1/CZdusCNN5rXb7zR+LnHj5vlD//4R7j6aigtbfz8IDF48GCWLVvGZ599xpgxYzh48CD33XcfKSkpjaf3FEIIP2ntesKNsoJyBjAKk2JStNLJkyfZtm0boaGhbNiwgaqqKs4880zPbnbbbWb/5psmU1ZDIiNNt3WvXpCXByNHwuefe1ZmAFx22WV88cUXfPDBB4wcOZLp06cTYiUpqaio4NSpUwGuoRCio/JpEHbQZt7tan+U1d5t2rSJ6upqBg8ezIoVKwC41tPntpdcAgMGwPffwyefNH1uQQFceCHs2WNyUD/8MLSREchKKa655hoKCwvJyHDmcOHFF19k8ODBPPfcc9JNLYTwO58GYaXUJKXUSOtt8KZgakMcXdFnnXUWBw8eZMCAAbVSOraIUjWt4Vdeafr8vn3hiy/goYdMRq358+GWWzwrO0BCQkLo3Lmz8/2HH37I7t27+c1vfkP//v158sknOXxYZtEJIfzD1y3hhUC+UmoSIP+yecG6desAOG0NkLr22mtblzf5jjtM/uglS2DfvqbP79QJ/vd/Ta7pQYPgN7/xvOwg8K9//Ytly5aRnJzMjz/+yNy5c+nXrx9z587lhyBZY1kI0X75OggnY3JNl2itF/m4rA6hsLAQgD179gBw3XXXte6G/fubhRsqKyE7u/nXXXYZbN4Mrrmqf/UrePppOHmydXXyo5CQECZOnEhBQQErVqzg8ssv5/Dhwzz55JN89dVXga6eEKKdUzqIE/W3BSkpKXr1av887q6urqZ79+4cPXoUpRQREREcOnSILl26tO7GH39spizFx8OuXeDSXdtsmzbVTGnq2xcefxxuvhnCGlyoK2itWrWKN998kz/96U/OXNxPPPEECQkJTJ48uflJUYQQHZZSqlBr3cTasn4amCW8Y9u2bRw9epSzzz6bvXv38t5777U+AIMZZHXeebB3r+mW9sSwYfDhh3D++bB7N9x+OwwdCosWQUVF6+voRxdffDHPP/+8MwDv27ePRx99lJtuuom+ffvy0EMPsXPnzgDXUgjRHrQ4CCulZiulrmxtwUqpkUqp+1t7n47E0RWdkpJC7969SU1N9c6NlYJ77jGv589vfLpSY/e4+mpYs8YM8kpIgOJiSEuDc8+FNjwNKDo6mj/96U8MGzaMffv2OVvFEyZMYMmSJTLFSQjhMU9awjnAjVZqygVKqf7NvVApFW0F8dWYdXvzPSi/w3IE4aSkJO/f/LbboF8/85zX09YwmHzTd9wBW7fC3/4Gw4ebrm5HF3dlpem6bkMiIyOZOXMm33zzDZ9//jm/+MUv6NKlC7m5uUybNo3vv/8+0FUUQrRRrXomrJSajAmmGsjSWrtNRFznvCVa65c8LjTI+POZ8OWXX87nn39OXFwc999/P3PmzPFuAYsWmZbr0KGwcaMJqK1VXW0ybnXrZt4vXgzTp5su8JkzzUIS3uhS97OysjLefPNNtm3bxp///GfALBoxceJExo4dy0033eR5EhUhRJvX3GfCXllTFxgAPAnsAF4ARgL9rWMFwAJggDfKCrbNX+sJV1VV6aioKOe6uQ8++KD3C6mo0Lp/f7OW8KJF3r+/1lr/8Y9ad+1qygCtbTatZ87U+ssv669t3MasWrXK+d8nJCREjx8/Xr/22mv68OHDga6aEMLPaOZ6wl4fHa2UGkdNqzdba92uu5z91RLeunUrQ4cOJTQ0lKqqKlauXMmYMWO8X9Bbb5lRzT17wrZtEBPj/TIOH4bXX4dXX4WioprjkydDTo73y/OTiooKPvzwQ958803ef/9957Pizp07M2HCBP7v//6PXr16BbiWQgh/CNjoaK11vjYLOExr7wHYnwoKCgCoqqqiR48eXHTRRb4paPp0Mwf44EF45BHflNG9u0nyUVgI69fDffdB795w6aU152zeDE88YX4ItBHh4eFMnDiRpUuXsnfvXrKysrjsssuorKzkq6++IjY21nnuhx9+yI8//hjA2gohgoHME24lf7WEZ82axYsvvgjALbfcwptvvum7wtavh6Qk8zz3s89qJ+TwldOnzaCtiAjz/uGHzUhtgMRE00r+2c/M4hGtyRAWAPv27WPr1q1cfvnlANjtdnr16kVVVRWjR4/m2muv5dprryUxMbF12c+EEEFD5gm3MytXrnS+9njBhuZKTISMDPPU9he/MN3HvhYWVhOAAVJTTdnR0eZHwbx55ofB2WfD3Lm+r48XnXHGGc4ADGYJytTUVEJDQ/nPf/7DQw89xMiRI+nbty/p6ens3r07gLUVQviTtIRbyR8tYbvdToz1bDYkJIQDBw7U6tr0iVOnYMwY02V8003w178GpgVaUWGWT3zvPfjgA/jhBzMFyrHgxMGD8NprZhpUYqLJg91GlJeXk5ubywcffMCHH37IPit39w8//EB8fDxguq179uxJcnKyM3mIECL4Nbcl3PZyCnZAjhzGw4YN41e/+pXvAzCYeb1vvgkpKWaw1vnnB6YFGh4O115rNq1h3braaTXz8+F+K+dLz55m6tO4cWYbODCou66joqKYNGkSkyZNorq6mjVr1lBQUOAMwAC//vWv2bVrF927d+cnP/kJ48aNIzU1laFDh0rXtRDtQXOGUMsW2ClKjzzyiAb0fffd5/Oy6nn3Xa2VMtOJ3n7b/+U3ZdUqre+4Q+s+fWqmPTm2gQO1rqwMdA09dvLkSZ2WlqYHDhzonPrk2Hr37q3fDsb/HkIIrXXzpyhJS7gN+Pe//w3gmylJTbn+enjySfOM+JZbTHfv1Kn+r0dDRo82m9awY4dpGefnm0UpbLaaBSS0NvmtBw0yo7AvuQRGjQrqRCHh4eFkZWUB8O2335Kfn09eXh4ff/wx+/fvp0ePHs5zn3/+eZYtW8all17KpZdeykUXXUR0dHSgqi6EaCZ5JtxKvn4mfOzYMWw2G6dPn+bVV1/l9ttv91lZDdLaTFeaP99k0Xr+eUhP9389WqK6GvbvNytDAWzfDuecU/ucTp3Mc+RRo8yUqXPP9X89PaC1Zvv27fTp04cIazDbDTfcwHvvvVfrvCFDhjBq1CjGjx/PbbfdFoiqCtFhyTPhduLzzz/n9OnTAOzYsSMwlVDKLE0YGgqPPWbSTW7YYNYODtaWZEhITQAG0wLetQtWroQvvjD7DRvMwLPCQrjzzppzs7NN2s6UFDMie8gQE7CDhFKKc+r8oMjKyuK2225j5cqVfPHFF6xdu5atW7eydetWjh075gzCdrud2bNnk5KSwsiRIxk+fDhdu3YNxNcQQiAt4VbzdUv4nnvu4dlnnwVg7dq1jBgxwmdlNctrr5n80qdOmeUPX3sNkpMDWydPHTliVn0qKIC77675QTFhAuTm1pzXubPpyk5MhJ/+1HTLB7mKigo2bNhAQUEB/fr145prrgEgLy+P8ePHO89TSjFw4EBGjBhBYmIis2bNIi4uLlDVFqLdaG5LWIJwK/k6CA8YMIBdu3Zx5pln8t133wXHiNivvzZzeLdtM63k2283LeWzzw50zbzj889NS7mgwIzGLimp+ezOO+Hll83rXbvgv/7L/BgZOrRmO/PMoB2VvXv3bpYuXUphYSHr169n8+bNzp4WgB9//NEZhB955BF++OEHhg4d6twGDBhAWJh0oAnRFAnCfuLLIPzdd9/Rp08fAO6//36eeuopn5TjkePH4Xe/gz/9yWS66tTJzCf+7/+GCy4I2iDkkfJy+OYbkzRk0CC40lpO+913YeLE+ud362a6sJctg759zbG9e03ikSDr+j116hRbtmxh3bp1bN++nccee8z52ZAhQ9hWJ21o586dGTx4MHfddRf3WGtQnzx5kqNHj9KjR4/g+JEoRBCQIOwnvgzCv//977nfmgO7adMmhg0b5pNyWmXHDnjoIbMGsePv0pAhJgf1NdeYZ6rtteVUWgqrVpk811u2mDWUt2wxCUQAjh2DyEjzevx4k3Skd29ISKi9JSWZru4g89lnn7Fp0ya2bNnCli1b2Lx5M3v27AHg0UcfZd68eQAsX76cq666iujoaBISEhg4cKBzP3DgQMaOHUt4eHggv4oQfidB2E98GYSHDh3K1q1bGThwYOAGZTVXSQn8+c8mwYcjCAFERZkpQcnJJtAkJsKAAbUTbrQ3hw5BcTFceGHNsSuvNF3c1spKtaSng5UXnK1b4be/hT593G8BbkkfPXqUbdu20aNHD/r16wfAO++8w4wZMzhy5Ijba+x2O927dwdMj86+ffvo27cvffr0oU+fPs7XNptNWtKi3ZAg7Ce+CsI//PADZ511FmFhYbz++uvcdNNNXi/DJ06fNnN0ly41e3c/HpQyI5f79jVbXBz06AGxsWaLiTFduhERZouMrHkdEWFGaYeFtakUlYCZNvXDD+YHS0mJCdQlJXD11XDrreac996DG25o+B7bt5sucTBLQX73HZxxRs0WHw+9evl9NLfWmtLSUoqLiykuLqakpITi4mJ+/PFH3n//fed5jh+W7syaNYu//OUvgJkXvWjRIuLj4znjjDOIj493vo5wzTEuRJCSIOwnvgrC8+fP5+GHH2bSpEksXbrU6/f3m+++M1OC1q83U4I2bIA9e0xA8oawsJqgHBpa/7WjZaVU7dcN7VtzjjecPg0nTpjn7I6VpVxfDxlS8+Nj1y7zbN6d7t3hrLNq7nnokPlzcffn5PqdfOyzY8coqaxkj7XtdtnP6dmT31mDwv5VXs41Vtd3va8WEsLahAT6W70pfzt8mO8qK+kZGkpcWBg9Q0PNFhaGLSREWtei5ZYuNYMsW0GCsJ/4IghXVVUxYMAA9uzZw/Lly5kwYYJX7x9wp0/D99/D7t1mO3TIbKWlNdvx4yYYOfaO1ydPQlWVuYdoNzRQRU3igq3A28BeYJ+1d7yuBI4AUda5qUBDC5dfD7xrvT4A/AqIA3paWyxgs7ZEQHKMCQDWroVWTgftsMk6lFJTALvWOq+Z51m9lEYAAA26SURBVCcBKUAJkACUNPdaX3n22WfZs2cPMTExjBs3LpBV8Y2wMOjXz2ytUV1dE5Dd7auqzHmObNKO1w3tW3NOMKm7NvO+ffD+++ZZ/cGDYLdDWVnNfulS85weTHrSDz5wf98RI8xqWmDuf/HF5pm/Y4uOrnl9ww1m0Q+Ab7+FnTvNY4WuXc3meB0RASEhKGr/YzQEmOemClprSg8fJspmcx679d13GbF9O4fsdg6UlnLQbudAWRkHy8roMWGCmT4H7N26lZzJkxv8Y/v8tdcYa815z3jmGV55911sUVHYoqKIiY42r6OjGdSnD7P/67+c1/3jk0+I7NKFqK5dieralW6RkURFRhLVtatM52qrBg70W1Ht6m+IUioVWAQ0K7mxUioByNRaj3c5tkQpVaK1LmnkUp85deoUjzzyCADDhw+X5esaExJitiDKZhWUzjvPrCrljuNHhKPL9n/+B26+uSZgHzoEBw6YYD1smLkXmM9OnjTbgQP17ztpUs25H3xggrs7UVEmaYrDxImmzG7dzGfdutVsV1yBSk2lB5gfFp9/DpGR3DF6tBn45jp+ID4ewsKorq52dt/3OfNM3n77bQ4ePMjBgwc5cOAAdrudsrIy7HY7vUePdqY2PaA1B0pLOVBaWq/KycnJzP797wGorq7m+uHDG/qT58UXXyTdSvH67rvv8vTTTxMVFUVUVBTdunUjMjLSuc2bN8/Zdb5ixQpOnDhBZGQkERERtc6z2WySF7wdaRdB2AqmGUAhUP//moalA1l1jmUBmTQzkHvbtGnTOH78OCEhIbz11luBqILoSOo+L01Obl4GtB494OhR05p2t7neo18/M12tvNxc47qPiqp936+/NoPX3AkJgdRU83rNGpg2reH67doF/foREhJi5q8vX05MZCTTHIP7wsPNdumlZkQ/mB8Dt98O4eEs6tKFP8+cyQmtOVFVxfGqKraffz57unQxI703boTNm6kKCWHexRdTfvIk9hMnsJ88SdmJExSdPEl5eblJCbpvH4SEsPebb1i3ciWnMV3qVS7V7dy5M48++qjz/X333cfGjRvdfrVf/epXPP/889Yf19dcffXVREREEB4eTpcuXQgPD3dur7zyCgOsXo5Fixbx1Vdf1frcsfXv3985+LO6upqlS5c679epU6da24ABA5zrm5eXl3P06FHCwsLqnSfP4punXQRhq9WaDqCUauAnt1tTqB+EVwO5bs71uXvuuYd33zVPsGbPns1ZjoE1QgQbpWq6lpv6ezptWsMBs+4AvQ8/NEG8brA+dswETIe4OJg82f24gRMnauZng2nFO7a6XFai4sgReP11AEKBrtbmMOSf/4TrrjNv5s2Dxx6jE/Bo3XvGx4Pd7lyqjjPPhP37mQXMqnPql1dfzWeXX25a7MuXmx8BnTrx6ZEjnOzalUqgUmsqtea/4uP5/tQpk9Fs3jwoLGRAWRmLSkupAue2Gvijdf/TBw/CU09BaCh9VqygfNu2Wue+DmwDxo4dy01nnw2rVlFZVcXKBx+sdV454GgSvPHGG9zarRucOsWnH37I/732GlVAtbVtB4ox62UfKS42P1iU4u7f/pa9+/ahQkMJCQsjNCyM7VFRqPBwbr31Vv77uuugvJxtO3bw+z/8gRDrvJDQUE6Fh3O4WzdCQ0N5asECzrR6ON5+5x3Wb9jgPC8kLIxTkZHQuTMDBw7kl7fcAqdOcer0aZ7/y1+c54V26mRed+5MaGgoP/nJTxjoxy5oV+0iCHtCKWXDPAOu1XLWWtuVUiilEvzRJa215uabb+a9997jxIkTAIwePZonn3zS10ULEXh1p5k1dzBMSgrk5DTv3L//3QRn18F+FRVmc3m2jM0Gr7xS81ndzTE1DEzX/OTJtT8/dcrse/YETF5upZT5wVBdXXu0e2UlVFcz+pJLGO3oql+2zKz8Bbj8NHD6ctUqM/0M4KqrYPly4oBJdc679vLLmfrEE1RUVHB2dDS88IK5xNpc9bvpJjb27WvmfC9fDvPnE05NEHfY37kzO0aMoLKykp49e5q0tQcP8jPgZ3XOnRcSwmPV1ea7f/GFeTQBPOfmO/XCDJj7yU9+YvK3f/QR51C/ZfQu4MhN9/iMGXDFFQBMtzZX1wD/su75y507Yf58OgP31jlvD2Dls+Nvf/sbAy++2DxaWbMGRo50U1vf6LBBGDMwEq21vYHPEzCDtepRSqUBaQB9HWkJPaSU4t///rczAE+cOJGc5v7jIoRomqMLuoe70OaiWze4447m3XP6dLM1x4YN7o9XV9ce2HfNNaYrvqHpaVYXMGAGm/361zUDEB1bdTW2s86qWXu8vByee875Wd3zb50+vebHxUcfwf3317qX43Vvm42vMzNryr/+ejh8uPb9tIbqav7nttt4dPp0k5P8yy/hJz+B6mpOnTpFtTVoUldXo6uryfvznzkZEUHv3r1h4UJITOT06dOcOHbMlG/9GQ0fPpy/3XEHVVVV9HLkGNCakydOUOVIgGOd+4vp07lkwADzb/OuXRAVha6uprKiAqU1yvoz79K1K3dNm0ZVVRUJCQn1x0f4SbuboqSUKgbSmxrhbI2KLtRa1/sTV0ppYHxzRkl7Y4rSRx99xL59+7jmmmvo5filK4QQwr+8GIg77BSltuiqq+p2EgkhhPC7AAwmC6ogbI1ybq7SRrqSW1KmzRv3EUL8//bu7saJJArD8HcyMGb3eiVPBrO7GYwzMEwE4AywiGA1k4EhAgQZzGwEDM4AS3u9y0AGZy/qtNW0239ju6vbfh8JgWvahRGlOq6/UwB21ZogXJzZ3eEtnyXd7vFXFuu9fUmLIBwbtso/BwDgKFoThGMncmNnc2MX9FwpY11ZXynjFkEYAHBUHbuG5uDulVJWll1GOQAAR3WKQbjIyf4TM+uZ2dc4XlSYaHn0PY5yAACOqjXT0fuIddy3Smd7e5JuzGwo6c7dy4du++X3xZT0xMxulNaYi1zSTEUDAI7uJIJw7G5eO3qNZ57VlM8kzY700QAAWOnkknU0zcz+lfTPAar6RdJ/B6gHyIl2jK47VBv+zd1/3fQQQbglzOxhm+wqQJvRjtF1TbfhU9yYBQBAJxCEAQDIhCDcHu9yfwDgAGjH6LpG2zBrwgAAZMJIGACATE7inHCXxb3GfyhdGDGQNN/mHmMgFzMbKeVXv6+U9yS9lvRJ0qNScpyxUtIc2jRaw8yuJA0lfZN0oXS3/LvKM430zQThjIqbo9x9WCr7aGZzsnahjaLzeq/6y1b6SjehFbeh/ZD0igCMNok2LHeflMq+xLW2t/G6sb6Z6ei8xpKmlbKpdrvSETg6MxuY2VRpRPC45tGhUma6C3d/VkkbC7TBuKbsvlLeWN9MEM5rpOWUmQ9RDrSGu8/dfVydslvxLFeBou2GNWU/Sn9urG9mOjqTWD9bGlXEpRIyswEdGQAclrvXLaWMFCPfpvtmgnA+fWlxsUSdgdKGAKBLBrFxS0pt/JEpabRZXG87K9aD1XDfTBDOZ+nOY6DjHiWpHHRjM4sIxGib+LI4lJZGx432zawJAziIWAuurhmz0RCt5O6f3H0saRK7oy9zfA6CcGax/gCcqrnSFDXtHK0U085TSX+Xy5tqswThfIo1hX65sPQfz3owOsXM3tQUF5tbBk1+FmBH95J6cYa40b6ZIJxJfPuaa3n9oa+UjYggjM4okhvE72VFR0Z7RnZx3v37mqnnXtN9M0E4r3ultGhll1EOdEZ0TOOaDupKaefpqp2mQJN6SgG22k6LL4/F2eDG+maCcF4TLaf/G0c50FZ91e8gfSyPhGP6bizpVVMfDFjH3WeSPtT8aCLptvQlsrG+masMM4tpkWtJn5W+jc3ItYu2iYD6VqmNjpRGEvdKlzOUjySN4pnnSoH6hqUVtE2cDb7Q5gscjt43E4QBAMiE6WgAADIhCAMAkAlBGACATAjCAABkQhAGsJfIMgTgCQjCAJ4sjiRlSXwPnAKuMgSwj2uRjAN4MkbCAH5iZpdm9nXLx3t1KSl3rAM4WwRhAFXXkjbmeo6p6I/71AGcO6ajAVRdabtE9WMt59fdtQ7grDESBlB1Kelu3QPF3aprbkfaWAcAgjAApWNGZjY1syJwvojXq3Y+v1RlKvoJdQBnjwscACyY2RtJ1+7++4bn7iS9WLEpa6s6ADASBvCzoTas5cZU9I81U9Eb6wCQMBIGsGBmLmm47t7UGOnOy/cI71JH3OU6lXQrqTjGdCHpm7vf7vP5ga4hCAOQtLjE/IukZ2tGuTKzO3cfPrWOGEl/d3erlL+R9NzdJ0/9NwBdw3Q0gMKVpFkRPIsd0GVmNtD6878b69Dq40ufJL3e9UMDXUYQBlD4U9JD6XVdQBwpTSXvU8dQ9ceXevELOBsEYQBlX6TFzUh1o9XrdevFW9bxUmnUWzVY8TxwsgjCAAp/SRpGOkq5+6z8w5iKfqh74451yN3nNe+9Fgk+cGbYmAVgK2Z2I+lui5HwujpeK+2cflEpv5T0nrPFODeMhAFs62qfAByW1oNjdPxWq/NQAyeLCxwAbBQj1U1T0eve31PapDWSNI8R8UJ1ZAycC6ajAWxkZlNJ0+oaL4D9MB0NYBt9AjBweIyEAQDIhJEwAACZEIQBAMiEIAwAQCYEYQAAMiEIAwCQCUEYAIBMCMIAAGRCEAYAIJP/Aec5rpe3x62CAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a1e6d5470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAEsCAYAAAD0PAtkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3Xl4lNXZ+PHvSUjCGiZhBw0QUEQWIQFBwQUIYFlKVRB9bd+qBYK1r21/KnHDrS4N1r619a0m2mrV2goBFdGqJIKCyhICyi4k7KuQjSUkITm/P84zk5msM5mV5P5c11wh9zznec7MAPec85xFaa0RQgghROCFBbsCQgghRHMlSVgIIYQIEknCQgghRJBIEhZCCCGCRJKwEEIIESSShIUQQoggkSQshB8ppaYrpbS7j2DX90KllEp15/1z9zghAqVFsCsgRBOXCYyvFksGpls/8wJeIyFEyJAkLIQfaa0LMYnYQSmVgEnCmVprScJCNGPSHS2EEEIEiSRhIYQQIkgkCQsRYpRSaUqpgmqxVGvwVrxTLMGKJTn9vtyKbVBKzfHwutOtcvbySdWuX1DLNRJqOUeu02Cz5dXq7O41XI5RSsU7XTdXKTW9jtdgr1+BdVxSbcd58LrrfT1CeEuSsBChZzlgq5bgkoBC66fdTACtdaaVlDZgBnqNB94F0pRSae5c0Cq/CHP/eoZ1nuoJx2Ydk4YZVBZv/W4/h/33DCDR6Ty2RlzjuWrXsD9mWMctonavWvVLsX6vN2nWV6eGXo8QPqG1loc85BHABzAP0EB8PcdoINX6s836PQ1Y7nTMBvvvQAGQVu0cSVa5JDfqVGC/XrXz2+uQap1rutPz9pjN+n268+9eXCOplmvMc4rZr5NQ/bhq5463v291HVdfnRp6PfKQhy8e0hIWIjTlUNXqTcJqodljSikbkAAssrpPbZgk7aC1zsS0nmdQD6vFbQPmVZuznGA9nDmP9M6t47kspdQ855a8h9fIruUazte1jyiPre91aTPyPBMYVtvzbtSpztcjhK/IFCUhQtO7mNYYmO7lTExiticPewLKpCpZ1zbdKQ/TIqyP/flETNJ2lu/8izZTrmqltS5USo236p1q1TVTaz3eB9eo87oNyKOOJNxQnRp4PUL4hLSEhQhNGQBWK3cYpts5j6qkMh7Ic4pB7ck2noYXBHE8r7XOq/bwKPlprTO11olADOa+bJJSap4vr+Gh+l5/g3Wq5/UI4ROShIUIQVZyLaRm16h9BS7nmL37Ntn5HNagI/tAp/qulWNdK7X6c1a3d2PqX6i1XoBpvQ/3xzUaYp03iWqLpTjV0e06VX89vq+taK6kO1qI0LUQk1hznFqL9lHCNqzlMK1u02TMaOhYTFd2PCa5ZFj3hhsyG3N/2T762UbVsprJ9RW0s6ZEpVjl7d3gCVTdq/b6Gm7UYZFVB/vrL8SMtK5LnXVSSm1o4PUI4TVJwkKEruWAPbEBjulINvufneLpSqk8TOJZhEkaKVbrrUFa6wylVKJVfjkmeS10vrYbFgJ9qJpalAcs0Fqn+/AaDXnXOnc81rSjBu5jN1SnOl+PEL6gtJYNRYQQQohgkHvCQgghRJBIEhZCCCGCRJKwEEIIESSShIUQQoggkSQshBBCBIlMUfJSx44dda9evYJdDSGEECFkw4YNJ7TWnRo6TpKwl3r16kV2dnbDBwohhGg2lFL73DlOuqOFEEKIIGmSLWFrzdxCN5frs+9KM4yqpeny3C0rhBBCNFaTS8LWrjOv0sAeqk7Hx2M29R7vFFuklLLvUCOEEEL4RZPpjlZKxSul0jAt2fyGjneSTM0F2dOoZWcVIYQQwpeaTBK29gBNbsTi6tOxNkt3km3FhRBCCL9pMkm4MazdaGq0nO27rlhd1UIIIYRfNOskDMRCVdKthSRhIYQQftPck7At2BUA4M9/hrFjYevWYNdECCFEADW50dGBoJSag9lsnbi4OO9P+PjjUFgIf/87vPCC9+cTQtSrtLSU/Px8Tp06RUVFRbCrI0JceHg47dq1IzY2lqioKJ+eW5Iw5t5wPV3SNViDv9IBhg0bpr2uwMCBsHo1tG3r9amEEPUrLS1l//79xMTE0KtXLyIiIlBKBbtaIkRprSkvL6e4uJj9+/cTFxfn00Tc3Luj7fOAY52D1oAt5+f9a+xY87O8PCCXE6I5y8/PJyYmho4dOxIZGSkJWNRLKUVkZCQdO3YkJiaG/HxPZsA2rFknYav1m0fNe8OxmBW3ApOEL7nE/Ny1KyCXE6I5O3XqFNHR0cGuhrgARUdHc+rUKZ+es1knYUsmZslKZwlWPDDsSTgnB5YuDdhlhWiOKioqiIiICHY1xAUoIiLC52MImmoSjqWWkc9KKZtSKtcaWGWXQs0lLpOteGD07Gl+5uWZQVpCCL+SLmjRGP74e9NkBmZZ93EfwszttQGpSqnxwHKtdYbToS73f7XWhUqpFKVUKrDeKp8a0HWjO3eGFi3g/HmYOBG0BvlPQgghmrwmk4St+7v1tl6tY2JqiedQc+nKwAkLg7g40xK+805JwEII0Uw01e7oC8/FF5ufBw4Etx5CCCECRpJwqLAn4c2bYcWK4NZFCCFEQEgSDhX2lbfuuw8mTICSkuDWRwjRpKWnpxMTE0NMTAwpKTXv5GVkZNRSKrRdiK+pydwTvuB17Wp+dukCl18OP/xQlZiFEMKHFixYwPLly0lNTaWwsJC0tDSGDx/O9OlmB9f09HTmzJnTwFm8M2PGDAoL3V6okLS0NOLj695TJxReU2NIEg4VkyZBjx5w2WUmCQshhB/k5eVx8uRJli9f7ojNmTOH2bNnOxKWc3LMzMwkKSmpxp+9tWjRIp+cB0LnNTWGdEeHij594KabJAELIfwqJyeH1NRUl5jNZmPmzJkUFhaSkZHhSFzgmiydk1x98vLymDGj+vIL/hOI1+QvkoRD1ZEjwa6BEM2SUsqjR2JiYq3lnSUmJnp83trKb9iwwevX55yMnCUlJZGZmUl+fn693b7uiI+P92lLtyGBeE3+It3RoaK83GxjWFAAH38MW7bA0aPmHrEQQviZzWZj+fLlAWvB+vqecG0C/ZoaQ5JwqGjRwixZWVYGo0aZbQ337JEkLESAae3d7qS1lfe2BeuLFrA7MjMzSUtL8+ocGRkZrF+/noceegibrcbqwQ6Bain74jX5kyThUKEUzJ8PLVvCLbfARReZlbSEECJAvB2gZL/3mpaWRn5+fr1JOFCCOejKHZKEQ8mjjwa7BkKIZionJ4fx48d7dY7p06dTWFgYMvdgffGa/E2aWkIIIWpM1cnLq3sPm/qeS09PJzk52ad1ayxfvSZ/kiQcSnbsgIULzaCsW26BXr3AxxtICyFEbXJzc126j/Py8khOTiY/P98RKywsJCUlhZycHJfjnKWlpXHLLbeExOpUjX1NgSTd0aHkjTcgNRWefhp27YJ9++C778xALSGE8CPnxATmXuqwYcMYN26cI5adnU1WVpZLYktJSSE+Pt4xT3f69OlkZ2eTkJAQmIrXo7GvKZAkCYcS+0joY8fglVegfXu45JLg1kkI0eQVFhYyc+bMGnGbzUZWVhYLFixwzP2tnqwWLVpEZmam4/fqi2YEizevKZCkOzqUdOxofubnw4gRZgnL8PDg1kkI0eTZbLY6F7ywJ6j169fXOdjKk/m+geLtawoUaQmHkthY87NaF4oQQgTTnDlzanTt2gV77eXGqu81BZIk4VDinIQrK+Gpp+D77+Htt2XOsBAiaGw2W51dthdiAob6X1Mgyf/socSehE+eNEk3PR3+9S8zQEsIIUSTIy3hUFK9O/rJJ80KWjExwauTEEIIv5EkHErsybagwHRHz54d3PoIIYTwK+mODiUtWphpSVpDUVGwayOEEMLP/JOElboJpWb55dxNnXOX9Llz8NFH8Prrwa2TEEIIv/BdElZqCEq9glLrgXggBqWyUepllBris+s0dc6Ds86dgylT4Je/hIqK4NZLCCGEz3l3T1ipaCAZmAnkAuloPdfpiOdRaigwF6USgXetY4q9um5T5twSttlgxgzo2hVKSswew0IIIZqMxiVhpW7GJN8YIA2th9V5rNYbgblWudnA5yh10iq3pFHXb8rmzIGpU6FfP/P7woXBrY8QQgi/8TwJK/UKUAAko/Uej8pq/SrwKkr1BpJRaiZa11zcszmrY5k1IYQQTY/nSdi1u7lxTPJ+0OvzNBfHj5t7xP37B7smQgghfEimKIWavDx4801YudL8vm6d2V3pttuCWi0hhBC+5/skrFS0NVK6lzVwS3hi1Sr4+c/hb38zv192mZk7bLOZ+cNCCCGaDO+TsFKzUeozlMq3BlxlAQuAdCAHpU6i1HqUuh+lenl9vaauXz/46U9h9Gjze3S0WUFr5UpQKqhVE0I0Henp6cTExBATE0NKSkqN5zMyMoJQq7pdaPV1V2NHR0cDDwNDgUXADLSuf4knpcYBD6KUfUT15426dlM3cqR5OJPkK4TwoQULFrB8+XJSU1MpLCwkLS2N4cOHO/bfTU9PZ86cOY0+/4wZMzzaYzgtLa3efX39Xd9gaszo6KHAHGCBR6Ojtc7CtJLtrecEtP6Dx9dvzrSWhCyE8EpeXh4nT55k+fLljticOXOYPXu2I6l5kkBrs2jRIq/KO/O0vs77G18Iex171h1tphbFo/XdHk9PclY1VemmRp+jqTp/HnbsgJycqth770HPnnDPPcGrlxCiScjJySE1NdUlZrPZmDlzJoWFhWRkZDiSW3WZmZmBqKILT+vr/AXAOXGHKs+SsNZ70HqxT66sdZE/FutQSiUopeYopZLsP90oY1NKzVNKxVt/jldKpbpT1ueKi81UpLFjq2KtW8P+/bBzZ8CrI0Rzo5RCVetxmjp1KkopPvzwQ0csPT0dpZRLN+jhw4dRStG9e3eX8omJiSil2LBhgyP2xBNPoJTiiSeecMQ2bNiAUorExESX8t27d0cpxeHDh71+fXUl2KSkJDIzM8nPz6+3azjQLrT6esp3WxkqNQStN/nsfI2qgooHUrXW451ii5RSeVrrvHqKxgKp1gOgEJittQ78175oa0B5cbHZzjAszAzS2rkTevcOeHWEEM2DzWZj+fLlzJgxw+tz+fqecG18Wd9gUtpX016Uehmt764Wa4+5f7w8EAlaKZUKrNdaZzjFkoBkrXWdn5SVvOOBbCC2gYTtYtiwYTo7O9uLWteiXTs4fRoKC830JCGEz2zfvp3+svBNrfr06UNubq5LzJ5Ms7OzKSwsJC8vj3nz5gWjejXUVt/k5GTS0tIASElJqdGV7S13//4opTbo+pZ0tjR+ipJSu62pR8+i1NhajzFdzs8DHeo8xremAznVYtlWvEFa60JPErDf2BOv7CkshAig6oOYcnJyyMzMxGazAaZrOCkpieTk5GBUr4ZQH3TlDm/mCc/AjHaeAGQCc6yk/BxKjXVZqMOMjLZ5VdMGKKVsmNZsvnNca11oPX/h3DSwJ2Hn7pyFC82OSsuWBadOQogmLScnh/Hjx7vEUlJSaoxAttls5OUFv61SW30vRI1PwlpvROsH0XoYWodhEvJCIBGTlAuspPwySs3CbHfoT7GmWrquGxENJeF4pdR06zFHKRW8nRSsb50uLeEtWyAjA775Jjh1EkI0adWn82RmZjru0xYWFjr+nJmZGRLJr3p96/tiEApfGuriu4FZkGd1PT8PgFIJwDhgOKa1/JwPr1Ubb1ra+QDV7iUvUkq5xAKmtu7om282q2ldeWXAqyOEaPpyc3Md3c4A8fHx5OebjsXs7GySkpLIy8tj0aJFITH1p3p98/LySE1NddQZzJeHlJQUhg8fHrIjqH2ZhF1nZ2udQ837syHJaj2nVwunWY8aSVgpNQcz4Iy4uDjfV8j+F8u5O/qKK8xDCCH8wDl5gUnC48ePJz09ndzcXAoLC8nPzw+JBAw165uUlMSwYcMYN26cI5adnU1WVpZLsg41vkvC5r5v0CmlbPV0SXsiD9NFXeN8Wut0rKQ9bNgw3++qIAOzhBABVFhYyMyZNe8Y2udA17eARzDUVV+bzUZWVhYLFiwgPj6eRYsWhXQChsasmOWrVa7MbkuzfHIuw97pH+t6GWWr9nwtVVG1jbe3f80KfB9GbfeEATIzYcECOHo04FUSQjRdNputziRbWFgYcomsvvra67p+/fqQ7YJ25llLWOs9KGVDqZeBVLTe26irKjUbs/zlQ40qX2vVdKFSKo+a94ZjgTqnHtkX+FBKZVQ7xp7MA39Hv7bR0QC//z1kZcHAgTBpUsCrJYRofvLy8i64qUBz5syp0V0dqjzvjtZ6I0o9CDxkbeawCFiI1sX1llNqCDAX6I1J4K96Xt0GZQLDcL0XnWDFa6W1zlNKJdeSpJOAHB91bXumru7om26CAQOgR4+AV0kI0TwlJCQEuwoes9lsIdd6r0vj7gmbbQsfBECpm4HXrNHQGrPko/0rSB8gBtOazMZsYbjRuyrXKwXzpcB5kFWy9bCqq2zABszylvbj8pVS8fZEbB2TDMz2Y13rNmgQzJoF11zjGv/lL4NSHSGEEP7h/cAss6FD1aYOZqnKeOzJuKF9hn3I6pJOsS9fadUjtZZWbmy1chn2OcJAB0yX9oygrZ41erR5CCGEaNJ8OUXJMEnXn63dBi5f/9Qoq3s5ppZ44OcDe0prMyhr9+6arWQhhBAXHN8nYeG9sjLYuhXKy10X59Aa+vaFs2chPx9ianyXEEK4QWtdY7tCIRrisw2PnEgSDkXHj0NCAnTrBs77h4aFwciRcO4cFBRIEhaiEcLDwykvLycyMjLYVREXmPLycsLDw316zoaTsFK9MYOU7AOv7AqBd9F6iU9rJMw84cGDoWvXms9lhcSaKEJcsNq1a0dxcTEdO3YMdlXEBaa4uJh27dr59Jx1J2EzwOohYDdmVPOeWo4Zh1KvAJ9JMvahtm3h22+DXQshmqTY2Fj2798PQHR0NBEREdI1Leqktaa8vJzi4mIKCgp8vlRx7UnYJOAZaP1gA7XLArIcK2lJIg6csjKQ7jQhPBYVFUVcXBz5+fns3buXioqKYFdJhLjw8HDatWtHXFwcUVFRPj137UnYjHB+ze2zmFZyzZayaDytoaQEoqLA+R5Ebi5cdx1ER8O2bcGrXz3Kyso4ePAgvXr1IizMrIy6c+dOjh8/zqWXXkqXLl2CXEPR3EVFRdGtWze6desW7KqIZs69taOVivZzPUR1w4dDmzbw3Xeu8a5d4dAhOHAAQuQb/Pnz5x1/rqyspEOHDvTp08dl2bgXXniBa6+9loULFzpimzZtYsqUKbz44osBra8QQoSK+pOw6WbeDRSi1EmUutHpuZtRapef69d8tWljflZfurJNG9izx6wr7eNRep7atWsXP/7xj7nuuuscsbCwMC655BLi4uI4c+aMI963b1+uvvpqLrnkEkds06ZNfPTRR6xdu9YRKy8vZ+TIkcydO5fKysrAvBAhhAiShkZH/x4zMjobGA7Mwex0vwStF6PUovqLi0aLtjofatvOsFevgFalLp07d2b16tWcOnWKgwcPctFFFwFm95Lqw/jnzZvHvHmum1VNnDiRxYsX06lTJ0dsy5YtrF27lvz8fEdXNsCTTz5Jp06duP3222lvX1tbCCEucA0l4WynfYIzgUyUmo1SY9H6c1ynLAlfCsE9hUtLS/n73//O3LlzUUrRvn17PvzwQy655BI6d+7sOM7deXTdunXjpptcd8a87LLL+OKLLzh16pQjdu7cOZ577jnKysq49dZbHfGcnBw6d+7sSP5CCHGhaSgJ19xBSOtXralJQ/1TJQHUn4TXrIHUVDOX+MknA1IdrTWTJ08mKyuLqKgo7rrrLgBGjRrl0+u0atWKa6+91iVWUVHBiy++SG5uLrGxVct+33333axbt46srCzGjh3r03oIIUQgNDQwK9tq+e5CqV6OaFXrOHSaaU1NfUn47Fl4//2ALtyhlCI5OZmePXsyaNCggF0XoE2bNiQnJ7NgwQJHrKKigu7du9OhQwdGjBjhiM+bN48xY8bw5ZdfBrSOQgjRGPW3hM3ewWYbQq331vJcYr3llepVo5xwT31JOCEB/vUvGDjQ79VwXmN3xowZTJkyhVatWvn9ug0JDw/nvffeo6KiwqX7++OPP2br1q0uiy+sWbOG3bt3M3HiRJf7z0IIEWy1t4RdW71Fde4B7LyKlnOZKgWO6U1KtUepWTLdyU31Dcyy2eDWW/2ehE+cOMF1113H7t27HbFQSMDOqt9//vLLL1m4cCEjR450xNLT0/nZz37G66+/7oidO3eO8vLygNVTCCFqU3sS1novSj1QR2KtSanZmD14q4sBhqNUtLUAyCLrIRpibwkXFwetCo8++iirVq3i/vvvD1odPBUbG8uMGTOIiIhwxK699lomTpzIj370I0fsrbfeIjY2lt///vfBqKYQQgD1dUdr/bx1PzgByMXs0Wvf5N4GxAMTMCOkU9C6tmyxHK0vsVrB9hZwB5/VvilraHT0+vWQmQlXX21W0PKDP/zhD2iteeyxx/xy/kC54447uOOOO1xiW7Zs4fTp09hsVd8dd+7cyUsvvcSNN94oA72EEAHR0D3hVwGskdBJmKRrAwqA9ZjkW9/grHSUGoLWm6z1qBOBcT6od9PXUBL+9FOYPx/uv99vSbht27akpaX55dzB9uKLL3L//ffTtm1bR2zZsmW89NJLFBcXO5Lw+fPn2bp1K4MHD5ZF/oUQPufefsLmnnDVfWGlete6q1LNcs9bydfcW1ZqDzK32D0NJeHrr4f/9/8gKcmnlz116hRvvfUWycnJPt83M9RcfPHFLr9PnDiR4uJiRo8e7YitW7eOUaNGMXr0aFatWhXoKgohmjj3knBNvwdmunWkc0tZ6z2OpCzq16EDXHkl9O1b+/OjR5uHjz311FP84Q9/YPPmzbz88ss+P38oGzhwIAOrDXY7fvw43bt3Z8iQIY7Y2bNnmTBhAuPGjeOJJ56QFrIQotEam4RnoJQG3gWy6rgfXLv6u6+FXffu4LSmcqCMHDmSXr168Ytf/CLg1w5FP/nJT5g2bRolJSWO2JdffslXX33FuXPneNJpsZQlS5aQmJhIz549g1FVIcQFSGndiN5hpSoxq2nZMN3LeZhlLT/D06R8gRs2bJjOzs4OzsUPHzbbGQ4ZAh07+uy058+fp0WLxn4/a/rOnj3LypUrqaysZMqUKQDk5+fTqVMnwsPDKSgooI21AUdlZaXLGthCiOZBKbVBaz2soeMa+79DBlrHAn2BuzH3i2cCizFzg3eh1Msuuy4Jz50/D/n5dW9ZOGsWjB8PPr5XKQm4fq1bt2bSpEmOBAxQWFjItGnTmDRpkiMBAwwfPpwJEyZw+PDhYFRVCBHiGvu/rdkoVus8IN16mAFbMB4zknomkIxSBcB0tF7hbWWbnfh4s2/w3r1QWxfniBFw6hRERXl9qbvuuosBAwZw991307p1a6/P19zEx8ezZMkSl9ixY8fYuHEjbdu2paNTT0VaWhotWrTgJz/5CR06yIw9IZqzxnVHu3121RuYC9wPxKP1Pv9dLDj82h09eDDs3w/ffAP9+/vnGsDmzZsZPHgwbdq0Ye/evS4JQ3jn+PHjbNu2jeuvvx4wy4D27NmTAwcOsGnTJq644grA7M3cpk0bunfvHsTaCiF8xd3uaP/2O5ppTCko9S6wAHdHVAvj228hACNvBw4cyLJlyzhw4IAkYB/r3LmzyzaPFRUVzJs3j7Vr17pshPHwww+TkZHBO++8w2233Qa4rtsthGiaAjNiROscIEijly5g7v4HXFQEXvRoKKWYPHkyc+fObfQ5hHtatGjBr371K9566y2XAVstW7akbdu2XHnllY7YH//4R4YOHcrChQuDUVUhRAD4NwkrdRNK2SdYyiId/nD55WZDh2PHGlW8srLSxxUSjfHWW2+Rn59PfHy8I7Zy5Uo2bdrE+fPnHbENGzZw//3388UXXwSjmkIIH/N3S3gBkIVSNyF7D3vumWegTx944426j4mOhlatzAAuD5WUlNCvXz/uu+8+ysrKGl9P4RMREREu3c8LFy4kKyuLG264wRH76KOPeOGFF1wGgZ06dYpFixZxrJFfxIQQwePvuSiJwC1AHlovaehgUU1REeTlwdGjdR/z8cemJdyIuaifffYZu3fv5osvvnDZdUiEhlatWtXYSGLy5MmcP3+eJKflSlevXs0tt9zCiBEjWLNmjSN+9OhRunbtGrD6CiE85++BWUXAq369RlPW0PrRALGxjT79tGnTWL9+PSUlJTIA6AKRmJhIYmKiSywqKorx48czatQoR+zEiRN069aN/v37s2XLFlkwRIgQJasyhLIA7Ck8bFiDI+hFiBs7dmyNFvP3339P+/bt6dKli0sCvuqqq+jatSuvvfaazFEWIgR4noSVegDYgNafe3VlM2ArCa3/4NV5mjJ3WsIFBXDnneanB4N1ZDnFpu3qq6/m5MmTnDx50hE7cuQIa9asoV27di77KD/55JOcOXOGX/7yl/Tq1SsItRWi+WrM/8IZwC3W0pTPoVQvt0sqFY1SD6BUNpAMZDXi+s1HdLT5WV8Sjo42ewt/+aXbLWatNUOHDuXWW2/lxIkTPqioCEXh4eEuc5S7du3Krl27WLRokWObSq01r776Ks8//zxnzpxxHPvpp5/y2muvcfDgwYDXW4jmxPOWsFmAw0woVepmIN3aUSmtzsFX5rhkzDSlRbixiojAvZZweDgsWQIXXwxOaxbXZ/v27Xz33XccO3bMpUUkmjalFH379qVvte0x09PTWbduHf2dVmVLS0vjvffe4/XXX+eOO+4AYN++fezdu5fExETatm0byKoL0WR5d09Y68XAYmt5ymSUWgAsB9IwuyzNBcZhdlhKthK4cJc7SRjgRz/y6LSXX345+/bt4/vvv5fNGpo5pRSTJk1i0qRJLvGpU6fSokULrrnmGkfsnXfe4eGHH+aee+7hpZdeAsyOUnv27KF///5ye0OIRvDN/8AmuT4IPIhS44CHMa3edLR+0CfXcJNSKgEYhtleMR7I01oyi4KUAAAgAElEQVRn+qucX/lxYFZcXBxxcXE+P69oGu68807uvPNOl1hMTAwJCQkuo7C/+uorJkyYwJgxY/j886phIidOnJAlUIVwg++bQVpnEaR7vUqpeCBVaz3eKbZIKZWnzY5PPi3nd+62hI8dg9deM39+5JF6D5X1iEVjzZ07t8bSpsXFxcTFxbmsg11QUECnTp247LLL2Lp1q6OFXFFR4bgXLYQwmlr/UTKmK9xZGpDqp3L+ZR+YVVxc/9rQJSXw6KNgdRHW5+mnn2bs2LGsWCE7Swrv3Xzzzezbt48//KFqksPu3btp06YNMTExLl3UgwcPZtiwYRw6dCgYVRUiJDW1JDwdyKkWy7bi/ijnXxERZknKigpwGrlaQ1wczJsHzz0HDawF/f7777NixQpKS0t9XFnRnDmvuDZ8+HCKiop47733HLGioiJ27NjBtm3b6NKliyM+a9YsrrnmGpeVvoRoTppMElZK2TD3cvOd41rrQuv5eF+WC5hu3aBr1/qTcFgYpKbCHXc0uHxlZmYmb7/9do3FHYTwpfDwcJdk2759ewoLC1m1apXLYMAvvviC1atX07JlS0fs5ZdfZsyYMS7rYwvRVDWZJAzEQlXyrEVdybSx5QIjNxeOHAGn/9C8ERMTw+23305kZKRPzieEu9q1a1djyc1vvvmGTz75hAEDBjhiX375JStXrqSgoMARW7NmDePHj+evf/1rwOorRCA0pSTc2AmvHpdTSs1RSmUrpbJ/+OGHRl7Wx4qKzKIdmcEd0C2EJzp27MjEiRNdurP/9Kc/sWzZMpfdo7766isyMzP57rvvHLHi4mImTpzI448/HtA6C+FLTSkJB4zWOl1rPUxrPaxTp07Bro6xbh3ccAM89VStT5eWljJ69Gjmz58vewiLkNalSxcmT55Mjx49HLGf/exnvP/++8yaNcsR+/bbb/nss89YtmyZS/n/+q//4p577iE/3+UOkxAhqcklYeseb8DK+d1vfws9eoDTIJdaDRwI114LTnM4na1evZqvvvqKpUuXyqIK4oLTuXNnpk2b5rLhyMCBA1myZAnz5893xM6cOcO///1v0tPTaeO0gtz8+fO588472bJlS0DrLURDmtJySfb5vLGY1boAl+Ra13zfxpYLjOJiOHwYGlrjuVu3ejdwuPrqq/nPf/5DeXm5jysoRHDExMRw4403usQiIiL49NNPOXDgAFFRUY74kiVL2LZtG3fffbcj9vbbb/PBBx9w55131lgxTIhAaTJJWGtdqJTKo+Y93ligsK5FNxpbLmCeeQaeeAK87PZu1aqVyz02IZqiyMhIxo8fXyP+t7/9jQ0bNrgsKpKZmUlGRgZjxoxxxNatW8eDDz7I1KlT+e1vfxuQOovmrckkYUsmZulJ5zm/CVbcH+X8r2tXz47Pz4fTp83cYSEEACNHjmTkyJEusYceeogxY8YwevRoRyw7O5sVK1bQs2dPR+zs2bMkJiaSkJDA22+/LSvOCZ9qajcHU4AZ1WLJVhww3cxKqVyl1BxPyl0Q3nkHOnSAFNdqf/bZZ/z617/mq6++ClLFhAg9/fr14+c//zl9+vRxxGbMmMFHH33k0m29ZcsWduzYwXfffeeSgEeNGsWoUaPYv3+/IyaDHoWnmlRL2OpaTlFKpQLrMXN8U2vpUo5tZLnAW70aXngBRo6skVxruOwyaN0aqn1Tf++993jllVfo2rWry+L7QghXnTp1qnF/eMiQIWzcuJFip41UysrKWLduHRUVFS4bVcyaNYuVK1fy0ksvOc5z5swZwsPDXRYkEcIuIElYKTUE092brbXe5M9raa1zqLkEpfPzhUCMp+WC5vhxeP99944dMsQM5Kq2SP5dd91F165dmTJlih8qKETTFhkZyZAhQ2rEDh48yPfff0/r1q0d8W3btrFnzx6i7eu+A2+88Qb33nsvKSkpPPvss4CZMrhr1y4uvfRSWTinmQtUS/hhzDrMeUopDaRorWVNOne4u5MS1Llk5fDhwxk+fLgPKyWE6NKli8vSnGCmAn7//ff07t3bETt27BhgWtl23377LSNGjGDo0KHk5FR99//kk0/o2bMn/fr1k6mEzUSgknAK8JzWeiOAUqp9gK574bN/o3YnCTurrGxwHWkhhG+1aNGCyy+/3CX21FNP8fDDD3P+/HlHrLi4mD59+rgcW1paypQpU9Bac/r0aVq1agXABx98QGVlJddddx2xsS530kQTEJAkrLXeU+13DzNKM+ZJSxjgww/h17+GCRPglVd455130Frzox/9SP4BCxEk1e8HJyUlsXv3bpeBXEVFRdxwww2cPXvWkYDBJPGcnBxWr17tGNPx8ccf8/XXXzN16lRGjBgRmBch/EKaSqHOnoSdBoXUq1072LMHNm8G4JlnnuGnP/0p27dv91MFhRCN5dzl3LlzZ5YtW8bnn3/ucsykSZO44YYbXFrNS5cu5ZlnnuGbb75xxHJycpgwYQIvvPCC/ysufKZJjY5ukjxtCY8YYRLwZZehteaOO+5g9erVck9YiAvU7373uxqx6dOn06lTJ6677jpHbNOmTSxfvtzl3vO5c+e46KKL6N+/P1988YUj6R8/fpzY2FiXbSVFcCitte9PqtRNQKzW+jWfnzzEDBs2TGdnZ/v3IlFRUFYGJSUg0xyEELU4duwYa9eupWPHjlx99dUAbN68mcGDB3PJJZfw/fffO44dMmQI27dvZ926dVxxxRUAbN26laKiIgYMGED79jJsx1tKqQ1a62ENHeez7mil1BCl1CtKKfs82xhru7+XrSlKorEaOzhLCNFsdOnShR//+MeOBAxmk4sjR46wePFiR0xrTUlJCWVlZfTq1csRf/HFFxk1ahRvv/22I7Z9+3aefvrpGl3kwne8SsJKqWil1ANKqWzgIWCR1nq41voPWuvnrW8B6cBcpdR6pdT9Sqno+s8qavC0S3rtWpg5k+9/9jN27tyJP3o7hBChTylF165dXdbMVkqxc+dOTp065dLivfjii0lISGDAgAGO2OrVq5k/fz5vvvmmI3bmzBnGjRvHr3/9a5dryf8zjdOoGwJKqZsxyzrGAGn1NbmtaUlzrXKzgc+VUietcjJX2B2eDs46fRoWLqRIKa58+22OHz9OyOx7LIQICW3btnX5ff78+S7bQgIMHjyYBx54gCuvvNIR2717N59//jmHDx92OTYhIYGzZ8/y8ccfO5YC3bdvHwAXXXQR4dUWERKGx0lYKfUKUAAkV5961BCt9avAq0qp3kCyUmqm1nqmp3VodjxtCQ8fzqn//V/++dFHjD53ThKwEKJRRowYUWMKVHx8fI1tUSsrK9m5cyclJSV07tzZEX/iiSd44403eOWVV0hOTgZgx44dLF++nCuvvFKmV9GIJKy1nuvtRa3k/aC352k2PE3C0dG0+81v+NNvfuO/OgkhmqV27drV2BY1LCyMkydPsmfPHtq1a+eIt2nThq5du9K3b19HbOXKldx7773cddddjiRcXFzMzTffzKBBg/jjH//oOLa8vJyIiAg/v6Lg8svo6OYkIKOjjx83mzLExIBMKRBCXMBWrFjBv//9b66//npuu+02wMxxTkxMZODAgWy21jgA6N+/P8XFxaxatYr4+HjAjOI+d+4c/fr1q9GlHkrcHR3t8yRsDbyKBwqBfK21mzcyL0wBScKNkLdqFb3WrCGsVSv41a+CXR0hhKhTUVERq1evRmvt2GimsrKS6Ohozpw5w+nTp2nTpg0A//3f/81bb73Fa6+9xi9+8QsANm7cyLvvvss111zD5MmTg/Y6nLmbhL1uVlmDrWZgdknSQB7mnjFAvFIqxoq9C2Rorfd6e01Rvx9++IEZ117LBkBfeilKkrAQIoS1b9++RvIMCwujqKiIQ4cOORIwmEFegwYN4rLLLnPEvv76a1JTUykoKHCc5+TJkyQkJHDFFVewdOlSx7Hbt28nJiaGLl26uOwPHSyNmqJkTU36vVLqU0zinaG1jtVad7CmKE2wHn3tMWAj8KBS6l2l1Fhfvogmb9ky+NGP4C9/cevw/fv3U9CjB+9374566CGQWw5CiAtQeHg4cXFxLrFnn32W7777zmVv9JEjR/LUU08xbdo0RywvL4/9+/c7RmjbTZ06lW7durFjxw5H7KOPPuJPf/pTUJb3bczo6KHAHGCBJ6OjtdZZQJZ1jtlKqQSt9R88vX6zdOQIfPIJdO/u1uGJiYnkHTzImTNnwOkbpBBCNEWJiYkkJia6xIYOHcquXbs4ffq0I6a1pnPnzhQWFrosVPLOO+/wzjvv8MYbb9C/f/9AVRvwMAlbU4vitdZ3e3NRrfWrSqn2SqmbZK6wGyZMMK1hpxGG7mgjCVgI0Uy1aNHCZVQ2mIVKvv766xrHTp06lZiYGBISEgJVvao6yeho74TawCytNVprs1B7WRmsWQMHDsDttwe7akII0WwEZe1oX51LNN7+/fuJiYnh5ptvhoICuO46SE4Gpw3FhRBChAZf7iecXD1gdTk/IAnaS2fOwBNPwCOPNHjo5s2bKS4uNvdBunSBH/8Y7rjDLGUphBAipDS6O1optRszFWk5kIkZIV3rvWKl1DhAa62b3FYcAemOPnvWDLCKijLbGTYwrP7o0aMUFxdz6aWX+rdeQgghahWI7ugZmNHOEzBJeI61U9JzSqmxzrslWSOjbV5cq3lr1cok4NJSk4Qb0LVrV0nAQghxAWh0EtZab9RaP6i1Hqa1DsMk5IVAIiYpF1hJ+WWl1CxANmpoLKWgQwfz55MnPS9/6hQsXw6Vlb6tlxBCCK/48p5wnrWH8AQrKQ/HJOUOmNby7314reYnNtb8zM+v85AjR44wduxYnnjiCdcnBg0y05y2bfNf/YQQQnjMl7sBLHL+RWudA+T48PzNmz0J19MSzsnJYcWKFTU3177uOtixw7SIhRBChAyfJWHrvq/wF3t3dD0t4auvvpoPP/yQyMhI1ydefx3CfNnpIYQQwhcas2LWUF+scmUN3LpFa/2at+dqFtzojo6JiXHsQOJCErAQQoQkj/53ttaK3mMNturV2ItaOy89JAnYA94MzLI7ccLsTSyEECIkeNxE0lpvBB4E5iqlPlVKzXKejlQXpdQQpdQr1s5LuVrrhxpR3+argZZwfn4+KSkpvPfee7WXf/556NwZ/vxnP1VQCCGEpxp1T1hrXYRJxCilbgZeU0olYLY1LATsmaIPYN9POBtIs5K48FQDLeGcnBwWLFjAVVddxY033ljzgAEDICICiov9WEkhhBCe8HpgltZ6MbDY/rtSqj0Qj5WMrYQtvNVASzguLo7HHnuM7nVtd5iUZNaSbt3aTxUUQgjhKV9OUQIcrWRp7fpaA0n40ksv5cknn6y7fGSkeQghhAgZPk/Cwk+GDIHMTKirpeuJo0eha1fvzyOEEMIrMnflQmGzwbhx0L9/jadOnz7NokWLyMvLq/8cZ8/C5ZdDnz5urUEthBDCvxpsCVtzg5MB+8Aru0LgXV/MGRbeycnJ4ZZbbiExMZF6d3Rq3drsxhQZCdu3Q0JC4CophBCihjqTsDXA6iFgN2ZU855ajhmnlHoF+EyScQA89RTs2wcvvght2zrCLVq0YNKkSQwePLjhcyxeDN26mZHSQgghgqrW/YStBDzD3cU0fLmSVmNZU6SGYaZDxWM2lMhsoIwNmANkYKZVxWJa/csbKmsXkP2E7Xr3hr17Ydcu6Ns3MNcUQgjhMXf3E661JWyNcHZ7NSv7SlruV8+3lFLxQKrWerxTbJFSKk9rXd+N0lgg1XqA6WKf7W4CDrj58+H8+aqR0t6oqIDCwqr5x0IIIQKuqQzMSgbSqsXSqEqu9RmPWVCkj9Y6Rmud4evK+cxdd8GcOS5JuKysjD179tTcOak+WVlw0UXwq1/5oZJCCCHc5VUSVkqN9VVFvDSdmtsmZlvxBmmtCxtoMYesnJwc4uPjueaaa9wv1KePmaa0bRtUVvqvckIIIerlbUt4fMOH+Jd1XzeeqqUyAZNYrefjg1Evv9i+HV57DVavdoR++OEHOnXqRF9P7hH36gVbtsCmTbLDkhBCBFG9U5SUUp8B4+p6GjNlKdgbMcRCVdKtRTxmsFZd4pVS9hZzLGapzXq7pJVSczADuoiLi/Ostt747DP4zW/gnntg9GgApk6dyrFjxygtLfXsXAMG+KGCQgghPFFvEtZaT1BKzdZav1rb80qpB/xTLY/YvCibD+CcdK0BXdSXiLXW6UA6mNHRXlzfM126mJ/VtiNUStGyZcvGnbOoCMrKoFMnLysnhBDCU+70RdbXigzNUcRusu4Fp1cLuzugK/A6dzY/jx0DQGvt2YCs6v75T+jRA55+2geVE0II4akGk7DWOque5+rdqEEp1cuTyiil4j142KqV9aZF7CwP00Xtq/P5TrWW8MaNG+nUqROzZ89u3PkuvxzOnIEDB3xUQSGEEJ6otTtaKdVLa73XkxPVUaZAKRWttS62LwACLNRa19jU1j7X14NLrgcWUNVSj8XM87Wfz55E62zJK6Xmaa0XVAvbB3jFU3PEdXDZk7DVEv722285efIkZ86cadz5hg6F3FyIbzpj14QQ4kJS64pZ4Ljfu8idZKyUmg2s11pvqhbvBfSxnrMn4oVa64le1rv69XMxK3zlOMXigQ1a65g6ysQDuZj5wXm1xGPqGezlENAVsyorzbrPFRVQWoqOiODAgQOUlZV5NjpaCCGEX3m1YhaA1vp5pdRsaznIXEyr0J6s7NOCJmBGSKfU1rrFLP94iVKqvVIq2or5Y4mmTMySlc4t1wTquWettc5TSiXXMj84CchxJwEHXFiYGUB19CgcP4666CLfjc7OyzPd3CNH+uZ8QgghGlTvPWGt9ata67uBLCARmAs8CMy0DknRWt9dRwIGSFdKDbGWwVTWOeqa8uSNFExXt7NkKw6Y7mmlVK41vcgu33kesdWFnQw08iZrAPToYX4eOuS7c37xBVx6KfziF7J4hxBCBFCDWxmCYwCWYxCWUqp3bbsq1VLueasLGq11kVJqD67bIfqE1rpQKZWilErF3Cu2ryVdvZUbW61chlJqujVPuAOmhT8jpFfP6tEDNmxg+9q1/Pbxx7nhhhv4zW9+4905r7oKLr7YtILPnnXZoUkIIYT/uJWEa/F7qlrD9bJawfY/77EnZV+z7gfXOZDK6l6ucX84pNeKrs1FFwGwdu1aPv30U9q1a+d9Eo6MhM2bJfkKIUSANTYJz1BKaeBdIKue7uganJOyaASrO3qSzUZGRgYdO3b0zXklAQshRMB5s3DwBGAxZhrSLqXUy0qpG50GYAl/sFrCnQsLufnmm7nuuut8e/4dO2DCBDNQSwghhF81NglnaK1jgb7A3Zj7xTOpJSn7qJ7C7ppr4PXX4be/9c/5n3kGli+HRx7xz/mFEEI4NLY72r7mch5mDeV0MAO2MDsrJWGScrJSqgCYrrVe4X11Bb17swf429/+xnUFBYwf7+ONrP74R4iJgSef9O15hRBC1NColrDWem4d8T1a63St9S1WS7kP8BqQqZTq6UU9hZNVq1bxzDPP8Morr/j+5J06wZ//bBKxEEIIv/LrZrJWUk4BhmOWmBQ+MGjvXh68+mpmTJ7s3wtVVsJLL0Fh6K1bIoQQTUFAdnS3pg8FaG3Hpm/owoU89/XX3Dp0qH8vlJIC//M/MHMmeLNbkxBCiFo19p6wW5RSNwF51prS8r+4r8yaBcXF0MEfK4A6ueceWLoU7r8flPLvtYQQohnyaxLGdEHHWBs8yPxgHzh69ChrevVi+PDh9LAvYekvvXrB1q3Qwt9/TYQQonnyd3d0Imat6Tyt9at+vlazkJmZyY033sivfvWrwFzQOQFv3AiTJ5sNJIQQQnjN3wOziqxNIDY1fLRwR0xMDOOuv55rO3eGDz4I3IW1Nt3gH38ML74YuOsKIUQTVud+wsI9Ad1P2G7nTrjsMujdm4pduwgLC0MF4p7t0aPw7LPw/PMQFWViWsv9YiGEqMbd/YQ9bgkrpR5QSo1tXLVczjNEKXW/t+dplnr1QitFxd69tImMpH379tx1110cOXLEv9ft2tXMIbYn4NJSGDYMFiyAsjL/XlsIIZqgxnRHZwC3WEtTPqeU6uVuQaVUtJXEszH79mY14vrN1smTJzl48CCr16/nIBCuNRdVVnLq1Clef/113njjjcBWaOlSyMmBN9+E8PDAXlsIIZoAj4e9WvsIzwVQSt0MpFs7KqVprZfUVsY6LhkzTWmRO010UdNbb73Fb3/7W1q3bs1SrbkY2PTvf3MkIYEPPviA++67L7AVmj4dPv3UDN6yJ+EjR2DqVDO3+IEHAlsfIYS4wHg1MEtrvVhrPQGTlK9USu22Nm4YopTqpZT6vVJqPTAMSNZaT9Rav+aLijdHJSUlREVFcfbsWQq6dweg7d69XHLJJdx///2O+8InTpxgy5Yt/q+QUmbHpbFOdyeWLoUNG2DVqqrY+fNmQ4h33pFFP5qKykrzuVZWVsXOn4fTp6GkpCqmNZw4YR7OTpyAAwfg3LmqWFER7N4Nx49XxcrLYft22LXLtXxuLmzb5lr+2DHYsgV++KEqduYMfPedOa+zLVvg229Nne327TMzAE6erIoVFJi/z7m5rq89O9vEne3cCevWmddhd/QorF0Le/dWxUpLTdnNm13L79gBmza5vn+1vaaSEnPsvn0135Ndu6Cioir2ww+wZw+cOlUVO3vWvPfOn4nWcPhwzZkPhYXm8ygtrYqdOwf5+ea9dX5PiorM5++stNQ85N993bTWPn0A44CFmL2Gx/n6/KH2SExM1IGyc+dOHRYWppVS+uDjj2sNWv/0py7HHD58WPfr10/HxcXpH374IWB1czh7VusPP9T6yy+rYtu3m7r27Ol67PTpWk+apPWBA1Wxr77S+u23td69uypWWKj1pk1a79tXFaus1PrwYfNwdvKkiZWUVMWKi7Xes0dr5/ejtFTrjRu13rzZtfy6dVpnZmp96lRVbMsWrTMyzOuwO3pU69de0/qDD1zL/9//ab1ggXkf7JYs0frRR7XesKEqtmmT1v/zP1q/+mpVrKxM6//+b63vuMP1nI8/rvWUKa7lFy/W+pprtP7jH6tihw5pfcUVWk+c6Fp+wgStL77Ytf6/+53WHTq4ll+zRus2bbQeM8a1fJcuWoeFub5/P/+5+Uz//veqWEaGid14Y1Xs3DkTi4hwPef48Sb+ySdVsT//2cTuuacqtnevicXFuZYfNMjEv/22KvbQQyb2zDNVsa+/NrGRI13Ld+hg4s6v6Y47ar6mxYtrvqbS0tpf04QJNV/TX/5S8zXt2VP7axo40Pevyd3PKZCv6eGHTezpp6tia9ZoHRmp9bXXupbv3Vvr6GjX13TvvVp37qz1O+9UxT7+WOuLLtJ6zpyqWFmZ1n37at2/v+s558wx9Vq9uir27rtaJyZqnZqqfQXI1m7kEJ9PUdJaZ2mzgcNMrbXc8/WhF198kcrKSu6880563HCDCVb7Nh0bG4vNZsNms3HG+ZtqoLRqBVOmmC0X7dq2hfnzYfZs12OzssyUJ/tAL4C//x1++lPznN3KlTBkCDjPja6ogO7d4eKLXc95++0mvsJp065//AN694bHH6+KHTkCQ4eaec/O7roLkpJc91P+5z9N1/vixVWxvDwzZevZZ13LP/YYzJvn2iJYuhSeftq0cuxyc+Evf4GPPqqKaW3ur7/9tus516yBZctMq8i5/qtWubYQz583rbtt21zLHzliWj7OrZmzZ02L7+zZqlhlpWndOMfs562sdG3NhIfXHAcQEQGtW7t+nmFhEBtbc3W3Tp2gRw/XY2026NMHOnasirVoYWYC9O3rWj4+Hvr3dy3fpQsMGOBavnVrGDTInNfZgAEweLDra4iLM3/PYmNd65SQYK5npxQkJpq4s379YPhwiHbaUr1LFxOLi6uKRUaasoMG1Sw/eLD5N1Tfa2rZ0hzrfE77e9Knj3nP7Tp2NIvutGtXFWvVyuxL7nxOMAMvu3RxjbVvbz6riIiqWFSU2eCldWvX96RdO/Nv3VlEhHk418nee+I8q6KiwgzurD7As6jIrA5YPXb8uGsvyJkzcPCga+u+stL0gFTvBdmzx/QuOP8bPXbM9E7s30/AuZOp5RH8lvChQ4d0q1atNKC3bNliWndgvj2Wl7sce/ToUX3mzJmA1Msr69dr/f77plVrl56u9a23av3FF1Wx5ctNy+fee6ti5eVad+2qdY8erue87TYTz8qqir35pvmG/sgjVbGjR7UePLhmq/Guu0xLcNeuqtg//2laDQsXVsVyc7W+807Xb/Naaz1/vtb33Wc+H7v339f6ySe1zsmpiu3erfWf/qT1smVVsYoKrd94w/QEOPvqK9PiPnq0KnbwoNYrV7rWs7TUXGPbNtfyBw6YFmVpaVXs1CnTunD+e3L+vKl39b87ZWXm/Xb+nITwRmWl69+nigrTe3XunOtxRUVaFxSY551jR4+6/j09c8b0lB0/7nrOnTvNw1lenmmZO/8bPXbM/H+0d6/3r82Cmy1hmSfspUDNEx47diwrVqxg4MCBbLa3fnv3Nveatm0zrYJaaK3Jy8ujT/WWgBBCCL/x2zxhERx5VvfolClTqoKDB5ufOTm1liktLeW//uu/SEhI4Pvvv/d3FYUQQnhIkvAFYs+ePXzzzTfMmzevKjhypPn5zTe1lomMjOT8+fMUFxdz0003BecesRBCiDpJEr5AKKUYOXIkMTExVcGrrzY/v/66zjJ///vfueyyy9i6dSuzZ89Gbj8IIUTokCR8AVi/fj3nnecz2g0fDunpNUfTOmnXrh1Lliyhbdu2/Otf/+LPf/6zH2sqhBDCE5KEQ9zhw4e58soradeuHXv27HF9snVrM+3n8svrPUf//v15/fXXAbjvvvv4yHlajDG875wAAA5ZSURBVBBCiKCRJBzicqxBV+Xl5XS3VslqjOnTp/Poo49SUVHBjBkzWLt2ra+qKIQQopEkCYe4KVOmsHfvXhYuXEiU88IEdkVFZnGIadMaPNdTTz3FnXfeSUlJCVOmTGHnzp1+qLEQQgh3yTxhLwVlP2Fn5eVmRRv7ursNzAcuLy9n2rRp/Oc//6FLly5kZmYycODAAFVWCCGaB5kn3ATk5OSwv6Fl1CIizPKHn39ecxm7Wg+PYNGiRSQlJXHs2DGuv/56gvolQgghmjFJwiHs9ttvp2fPnjz66KP1H/izn8GYMa7ru9ajTZs2fPjhh0yePJmTJ09Kt7QQQgSJx/sJi8DIy8tjx44dgFl0w22nT9dcRL0WLVu2ZMmSJSxbtoybbrrJEddaO7ZEFEII4V/SEg5Rb1tzf6+++mruuuuuhgvs2wcTJ5odjNy8zx8ZGemSgNetW8fw4cP59ttvG1VnIYQQnpEkHIK01o4k/Nhjj3HRRRc1XKh9e1i/Hr74wnV7PA88+eSTbNiwgcXOW/YJIYTwG0nCIeibb75h165ddOnShXHjxrlXyGYD+73jOXNc99V007vvvsvDDz9MSkqKI5aVlcWhQ4c8PpcQQoiGNbkkrJSarpRK8uD4BKXUHKVUkv2nP+vnDnsSHDJkCC1aeHDb/t57YfRos4n7tGlw6pRH123bti3PPPMMbdq0AaCkpITbbruN3r17M336dJYuXUp5eblH5xRCCFG3JpWErQT6qgfHxwOpWut0rXWm1jodSLbiQbF7925Wr14NQGJiomeFW7SAd96Biy82mzqMHg1btza6LkVFRYwZM4aKigoWL17MtGnT6Ny5M7fccgtvvPEG+/btkw0hhBDCC01isQ4raaYAG6yfyVrrTDfKpQLrtdYZTrEkq/wMd67ty8U6tNbceOONfPDBB4wZM4aPP/6Yli1ben6i3Fy44QazeEdYGNx+O/ziFzBqlEnUHjp06BD//Oc/efPNN9laLal36dKFESNGkJCQQL9+/bjqqqvo2bOn53UWQogmxN3FOppEEnamlMrF/SScC4zXWuc5xWxAgdbarXk6vkrCWmt+97vf8fjjj9O2bVt27NhBjx49Gn/C4mJ48EGzy1JFhYm1aQNXXgkPPwxJVq/7gQNw9Cj07AmdO5tYZaX5GVazoyQ3N5f//Oc/fPLJJ3z99dcUFBS4PP+///u//OY3vwFg2bJl/PWvf+Wmm25i1qxZAJw4cYIVK1bQvn17x6NVq1ZERUURFRVFy5YtiYqKIqyWawshxIXC3STcbOcJW8k2Hsh3jmutC5VSKKXinZOzP50/f57o6GhKSkoAePXVV71LwADR0fDXv8J998Hf/gYZGbBrF6xYYe4d2731FjzyiEnYzz1nYp98ApMnQ2QkREWZR4sWEBZGn/BwfhUezq/CwtAdO1IeE8PZ0lJWXn45/2jThmHDhsHq1fDAA3SKiOA/q1YxYMAA+OEH+MlPiCgqovvWrWigHDgB2L8GVv86uD4igueio7n33nt57Je/hBkzOKEUI/btY9SoUbz55pswYwbnf/iBDRs2AGYPZW0+v1rfluKICJ7o149hw4bxpz/9CW69lfPHj/OTkhIqY2L4+OOPYf58WLuW7777jtKyslrPU/3szwwcSGVMDB988IGj/FMtWrD6/Hn+8Y9/0G35cnj7bfLy8mp8canLP/r0YYfN5lL+g+ho/q+4mEceeYTroqLgscc4evSo24PnMrt1I6t7d5fyW9u25benT3PLLbcw68Yb4bbbKC4uZteuXW6dc0f79vyjb1+X8vlKcavWDB06lNTUVMf77O70t1MRETw7eHCN8jMrK9GxsSxZssTxPm/esoWy0lK3zvvsoEE1yj8ZHs5XFRXyOTXicyqOiOA558/pttvM51RR4fo+r1vH5s2b6/z3VN1zAwfWKP9EWBhfO39O//wnuR58Tm/26WPeA6fy77drx1+dP6fHH4eXX4b4AN6R1Fo3qQeQCyS5cVy8efm1PqfrOwcwB8gGsuPi4rQvhIeHa0A/8sgjPjlfrY4c0frDD7U+frwq9sorWickaJ2eXhV7/32tzWxj9x9PPlmj/JmkJL106VL97bffan3woMfn/MB8DjolJcVRvqRDBw3opKQkc62uXT065yHrnOPGjXMp3xW0zWYzscmTPa5rbeV/0aWLBnRubq7Wjz3m8TknW3V1Lv/u5ZdrQL/77rtaL13q8Tkfs87pXH7XZZdpQM+bN0/rQ4c8PudS65zO5c/Fxtb6PsvnJJ9TyH9OmzY1/v9YJ0C21g3nrGbbHa2USgA26Fq6nZVSGtNN3WCXtq+6oxcvXsy1115Lp06dvD6XT1RWQlkZlJbCuXOmS7uy0vys7c+dOkGXLqZsfj7s2AExMdC/v4mVloL9fXL+O2f/s1NMV1Zy/vx5Stu0oaRvX1q2bEm7yEj4+mvOVVZysGdPWrZsaeZPr1pFRUkJe/LyqKysRFdWUmk9qPZ3W2tNZUQExUOH0r59ewYNGgQrV3L+7Fm+iYwkrFUrRo0aBTk5cOIE27dvp8z65l7XvxN7vGDwYMJateL66693lF9TWkpxVBSjR4+m9eHDkJdHbm4uhYWFbn0ExX36UG6zuZTfXlbGgchIBg8eTNewMNi0yaMW1tlu3Sjp3t2l/KGyMrZGRtKrVy8ujYuDL7/0qIVV1r49p/r2dSlfdO4ca1u2JDY21vSOWO+zuy2syogICgYPrlF+RWUl4a1bM3bsWMf7vHnzZsfn1JD8QYNqlJfPyYvPqUULCq64oqr8ihXmc9La9X3+4QePeizyBw+uUX5NaSnFLVtWvc+5uR71WJzq25cy588pN9d8TlFRVe/zxo1w1VWmJ9FLck/4AkvCQgghmo4L8p6wh1OD8rXW7n1Vrf+aNl+cRwghhPBUyCRh+5xdD4qsBxZ4cUn7oKtYwJGErQFbzs8LIYQQfhEySVibkchuzc310fUKlVJ58P/bu7/bNo4gjuO/6YCm8x6A6kBJOhA7oK0KLHZgwhUEUge0KwjsDqRUEJkdmECeA8vqYPKwc8T5eDySonR7R34/gBFrdVyckcUO9+9oUPnVUNKjt7QzGgBwuk79MOadpOqc/XmUAwDwoo4xCA+1PrqVmQ3M7JuZXZWKZ1offU+jHACAF9WZ6ehDxDruB6WzvwNJ12Y2lnTrpSsplQL0SkxJz4rrK+Pz10xFAwDacBRBOHY3N45e45lXNeULSYsXejUAADY6unPCbTOz/yT9+wxV/aJ0iyPQZ7Rj9N1zteFf3X3r7UsE4Y4ws/tdDnYDXUY7Rt+13YaPcWMWAAC9QBAGACATgnB3fMz9AsAzoB2j71ptw6wJAwCQCSNhAAAyOYpzwn0WKRV/V0oYMZK03CWFIpCLmU2U7le/q5QPJF1J+iLpQelynKnSpTm0aXSGmV1IGkv6LulMKa3tx8ozrfTNBOGMisxR7j4ulX02syW3dqGLovP6pPpkK0OlTGhFNrRHSe8IwOiSaMNy91mp7Guktb2Jn1vrm5mOzmsqaV4pm2u/lI7AizOzkZnNlUYEDw2PjpVupjtz91eVa2OBLpjWlN1VylvrmwnCeU20fmXmfZQDneHuS3efVqfsNjxLKlB03bim7LH099b6ZqajM4n1s7VRRSSVkJmN6MgA4Hm5e91SykQx8m27byYI5zOUVokl6oyUNgQAfTKKjVtSauMPTEmjyyK97aJYD1bLfTNBOJ+1nMdAzz1IUjnoxmYWEYjRNfFlcSytjY5b7ZtZEwbwLGItuLpmzEZDdJK7f3H3qaRZ7I4+z/EeBOHMYv0BOFZLpSlq2jk6Kaad55L+Lpe31WYJwvkUawrDcmHpfzzrwegVM3tfU1xsbhm1+S7Anu4kDeIMcat9M0E4k/j2tdT6+sNQ6TYigjB6o7jcIP5bVnRktGdkF+fdfzRMPQ/a7psJwnndKV2LVnYe5UBvRMc0remgLpR2nm7aaQq0aaAUYKvttPjyWJwNbq1vJgjnNdP69X/TKAe6aqj6HaQP5ZFwTN9NJb1r68WAJu6+kPRXza9mkm5KXyJb65tJZZhZTItcSvpH6dvYgrt20TURUD8otdGJ0kjiTik5Q/lI0iSeea0UqK9ZWkHXxNngM21P4PDifTNBGACATJiOBgAgE4IwAACZEIQBAMiEIAwAQCYEYQAHiVuGADwBQRjAk8WRpCwX3wPHgFSGAA5xKS7jAJ6MkTCAn5jZuZl92/HxQd2VlHvWAZwsgjCAqktJW+96jqnoz4fUAZw6pqMBVF1ot4vqp1q/X3ffOoCTxkgYQNW5pNumB4rcqg3ZkbbWAYAgDEDpmJGZzc2sCJxv4udNO5/fqjIV/YQ6gJNHAgcAK2b2XtKlu/+25blbSW82bMraqQ4AjIQB/GysLWu5MRX92DAVvbUOAAkjYQArZuaSxk15U2OkuyznEd6njsjlOpd0I6k4xnQm6bu73xzy/kDfEIQBSFolMf8q6VXDKFdmduvu46fWESPpH+5ulfL3kl67++yp/wagb5iOBlC4kLQogmexA7rMzEZqPv+7tQ5tPr70RdLVvi8N9BlBGEDhD0n3pZ/rAuJEaSr5kDrGqj++NIg/wMkgCAMo+yqtMiPVjVYvm9aLd6zjrdKot2q04XngaBGEART+lDSO6yjl7ovyL2Mq+r7ug3vWIXdf1nz2UlzwgRPDxiwAOzGza0m3O4yEm+q4Uto5/aZSfi7pE2eLcWoYCQPY1cUhATisrQfH6PiDNt9DDRwtEjgA2CpGqtumops+P1DapDWRtIwR8Up1ZAycCqajAWxlZnNJ8+oaL4DDMB0NYBdDAjDw/BgJAwCQCSNhAAAyIQgDAJAJQRgAgEwIwgAAZEIQBgAgE4IwAACZEIQBAMiEIAwAQCb/A5mPR6KWm5n1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a1e71f898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Plots jz1t, jz2t, jz12t in the Dicke basis for different parameter values\n",
    "\n",
    "#spin normalization constants\n",
    "j2_max = (0.5 * N + 1) * (0.5 * N)\n",
    "jmax = 0.5 * N \n",
    "j1max = 0.5 * N1\n",
    "j2max = 0.5 * N2 \n",
    "\n",
    "#plot graphics properties\n",
    "plt.rc('text', usetex = True)\n",
    "label_size = 20\n",
    "fig_size = (14, 7)\n",
    "lw = 2\n",
    "lw1 = 1*lw\n",
    "lw2 = 1*lw\n",
    "lw3 = 1*lw\n",
    "\n",
    "fig1 = plt.figure(figsize=(7,4))\n",
    "plt.rc('xtick', labelsize = label_size) \n",
    "plt.rc('ytick', labelsize = label_size)\n",
    "\n",
    "plt.plot(t/td0, jz1t_0/j1max, '-k', label = r\"$\\gamma_\\Downarrow$ Only\", linewidth = lw)\n",
    "plt.plot(t/td0, jz2t_0/j2max, '-r', linewidth = lw)\n",
    "\n",
    "plt.plot(t/td0, jz1t_gE/j1max, '--k', label = r\"$\\gamma_\\downarrow=\\gamma_\\Downarrow$\", linewidth = lw2)\n",
    "plt.plot(t/td0, jz2t_gE/j2max, '--r', linewidth = lw2)\n",
    "\n",
    "\n",
    "plt.rcParams['text.latex.preamble']=[r\"\\usepackage{xcolor}\"]\n",
    "plt.xlabel(r'$t/t_\\text{D}$', fontsize = label_size)\n",
    "\n",
    "#make double label y-axis - STARTS\n",
    "left = -5.5\n",
    "center = 0\n",
    "yshift = -0.4\n",
    "#label Jz1\n",
    "plt.text(left, center+yshift,r'$\\langle J_{z}^{(1)}(t)\\rangle$,',\n",
    "         horizontalalignment = 'right',\n",
    "         verticalalignment='center',\n",
    "         color = \"k\", rotation='vertical',fontsize = label_size)\n",
    "\n",
    "#label Jz2\n",
    "plt.text(left, center-yshift, r' $\\langle J_{z}^{(2)}(t)\\rangle$',\n",
    "        horizontalalignment='right', verticalalignment='center',\n",
    "         color = \"r\", rotation='vertical',fontsize = label_size)\n",
    "#make double label y-axis - ENDS\n",
    "\n",
    "plt.xticks([0, (tmax/2)/td0, tmax/td0])\n",
    "plt.yticks([-1, -0.5, 0, 0.5, 1])\n",
    "plt.legend(fontsize = label_size)\n",
    "plt.title(r'Two ensembles', fontsize = label_size)\n",
    "\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "## Second Figure\n",
    "plt.rc('xtick', labelsize = label_size) \n",
    "plt.rc('ytick', labelsize = label_size)\n",
    "fig2 = plt.figure(figsize=(7,4))\n",
    "\n",
    "plt.plot(t/td0, jz1t_gCEi/j1max, '-.k', label = r\"$\\gamma_{\\Downarrow,i}=\\gamma_\\Downarrow$\",\n",
    "         linewidth = lw3)\n",
    "plt.plot(t/td0, jz2t_gCEi/j2max, '-.r', linewidth = lw3)\n",
    "\n",
    "plt.plot(t/td0, jz1t_gD/j1max, ':k', label = r\"$\\gamma_\\phi=\\gamma_\\Downarrow$\", linewidth = lw1)\n",
    "plt.plot(t/td0, jz2t_gD/j2max, ':r',linewidth = lw1)\n",
    "\n",
    "plt.rcParams['text.latex.preamble']=[r\"\\usepackage{xcolor}\"]\n",
    "plt.xlabel(r'$t/t_\\text{D}$', fontsize = label_size)\n",
    "\n",
    "#make double label y-axis - STARTS\n",
    "#label Jz1\n",
    "plt.text(left, center+yshift,r'$\\langle J_{z}^{(1)}(t)\\rangle$,',\n",
    "         horizontalalignment = 'right',\n",
    "         verticalalignment='center',\n",
    "         color = \"k\", rotation='vertical',fontsize = label_size)\n",
    "\n",
    "#label Jz2\n",
    "plt.text(left, center-yshift, r' $\\langle J_{z}^{(2)}(t)\\rangle$',\n",
    "        horizontalalignment='right', verticalalignment='center',\n",
    "         color = \"r\", rotation='vertical',fontsize = label_size)\n",
    "#make double label y-axis - ENDS\n",
    "\n",
    "plt.xticks([0, (tmax/2)/td0, tmax/td0])\n",
    "plt.yticks([-1, -0.5, 0, 0.5, 1])\n",
    "plt.legend(fontsize = label_size)\n",
    "plt.title(r'Two ensembles', fontsize = label_size)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have studied the dissipative dynamics of two ensembles of TLSs, exploring the possibility of the systems to undergo local dephasing, collective emission of the single ensembles, collective emission of the two ensembles coupled to the same reservoir and local de-excitations. We have found that in the general casse spin exchange between antisymmetrically prepared ensemble is transient [1]. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### References\n",
    "\n",
    "[1] N. Shammah, S. Ahmed, N. Lambert, S. De Liberato, and F. Nori, https://arxiv.org/abs/1805.05129\n",
    "Open quantum systems with local and collective incoherent processes: Efficient numerical simulation using permutational invariance\n",
    "    \n",
    "[2] Y. Hama, W.J. Munro, and K. Nemoto, *Phys. Rev. Lett.* **120**, 060403 (2018)\n",
    "    Relaxation to Negative Temperatures in Double Domain Systems\n",
    "    \n",
    "[3] M. Xu, D.A. Tieri, E.C. Fine, J.K. Thompson, and M.J. Holland, *Phys. Rev. Lett.* **113**, 154101 (2014)\n",
    "    Synchronization of Two Ensembles of Atoms\n",
    "\n",
    "[4] B.A. Chase and J.M Geremia, *Phys. Rev. A* **78**, 052101 (2010)\n",
    "    Collective processes of an ensemble of spin-1/2 particles\n",
    "\n",
    "[5] J.R. Johansson, P.D. Nation, and F. Nori, *Comp. Phys. Comm.* **183**, 1760 (2012) \n",
    "    http://qutip.org"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "QuTiP: Quantum Toolbox in Python\n",
      "Copyright (c) 2011 and later.\n",
      "A. J. Pitchford, P. D. Nation, R. J. Johansson, A. Grimsmo, and C. Granade\n",
      "\n",
      "QuTiP Version:      4.3.0.dev0+cdc2204a\n",
      "Numpy Version:      1.13.3\n",
      "Scipy Version:      1.0.0\n",
      "Cython Version:     0.27.3\n",
      "Matplotlib Version: 2.1.1\n",
      "Python Version:     3.6.3\n",
      "Number of CPUs:     2\n",
      "BLAS Info:          INTEL MKL\n",
      "OPENMP Installed:   False\n",
      "INTEL MKL Ext:      True\n",
      "Platform Info:      Darwin (x86_64)\n",
      "Installation path:  /Users/nathanshammah/Dropbox/GitHub/qutip/qutip\n",
      "Please cite QuTiP in your publication.\n",
      "For your convenience a bibtex file can be easily generated using `qutip.about.cite()`\n"
     ]
    }
   ],
   "source": [
    "qutip.about()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}