{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Blockade Interaction in a Magnetic Field"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interaction between Rydberg atoms is strongly influenced by external electric and magnetic fields. A small magnetic field for instance lifts the Zeeman degeneracy and thus strengthens the Rydberg blockade, especially if there is a non-zero angle between the interatomic and the quantization axis. This has been discussed in M. Saffman, T. G. Walker, and K. Mølmer, “Quantum information with Rydberg atoms”, [Rev. Mod. Phys. 82, 2313 (2010)](https://dx.doi.org/10.1103/RevModPhys.82.2313).  Here we show how to reproduce [Fig. 13](https://journals.aps.org/rmp/article/10.1103/RevModPhys.82.2313/figures/13/medium) using pairinteraction. This [Jupyter notebook](https://github.com/pairinteraction/pairinteraction/blob/master/doc/sphinx/examples_python/comparison_to_saffman_fig13.ipynb) and the final [Python script](https://github.com/pairinteraction/pairinteraction/blob/master/doc/sphinx/examples_python/comparison_to_saffman_fig13.py) are available on GitHub."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described in the [introduction](introduction.ipynb), we start our code with some preparations. We will make use of pairinteraction's parallel capacities which is why we load the `multiprocessing` module if supported by the operating system (in Windows, the module only works with methods defined outside an IPython notebook)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "# Arrays\n",
    "import numpy as np\n",
    "\n",
    "# Plotting\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Operating system interfaces\n",
    "import os\n",
    "import sys\n",
    "\n",
    "# Parallel computing\n",
    "if sys.platform != \"win32\":\n",
    "    from multiprocessing import Pool\n",
    "from functools import partial\n",
    "\n",
    "# pairinteraction :-)\n",
    "from pairinteraction import pireal as pi\n",
    "\n",
    "# Create cache for matrix elements\n",
    "if not os.path.exists(\"./cache\"):\n",
    "    os.makedirs(\"./cache\")\n",
    "cache = pi.MatrixElementCache(\"./cache\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by defining some constants of our calculation: the spatial separation of the Rydberg atoms and a range of magnetic field we want to iterate over.  The units of the respective quantities are given as comments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance = 10  # µm\n",
    "bfields = np.linspace(0, 20, 200)  # Gauss"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we use pairinteraction's `StateOne` class to define the single-atom state $\\left|43d_{5/2},m_j=1/2\\right\\rangle$ of a Rubudium atom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "state_one = pi.StateOne(\"Rb\", 43, 2, 2.5, 0.5)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we define how to set up the single atom system.  We do this using a function, so we can easily create systems with the magnetic field as a parameter.  Inside the function we create a new system by passing the `state_one` and the cache directory we created to `SystemOne`.\n",
    "\n",
    "To limit the size of the basis, we have to choose cutoffs on states which can couple to `state_one`.  This is done by means of the `restrict...` functions in `SystemOne`.\n",
    "\n",
    "Finally, we set the magnetic field to point in $z$-direction with the magnitude given by the argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup_system_one(bfield):\n",
    "    system_one = pi.SystemOne(state_one.getSpecies(), cache)\n",
    "    system_one.restrictEnergy(state_one.getEnergy() - 100, state_one.getEnergy() + 100)\n",
    "    system_one.restrictN(state_one.getN() - 2, state_one.getN() + 2)\n",
    "    system_one.restrictL(state_one.getL() - 2, state_one.getL() + 2)\n",
    "    system_one.setBfield([0, 0, bfield])\n",
    "    return system_one"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To investigate the $\\left|43d_{5/2},m_j=1/2;43d_{5/2},m_j=1/2\\right\\rangle$ pair state, we easily combine the same single-atom state twice into a pair state using `StateTwo`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "state_two = pi.StateTwo(state_one, state_one)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Akin to the single atom system, we now define how to create a two atom system.  We want to parametrize this in terms of the single atom system and the interaction angle.\n",
    "\n",
    "We compose a `SystemTwo` from two `system_one` because we are looking at two identical atoms.  Again we have to restrict the energy range for coupling.  Then we proceed to set the distance between the two atoms and the interaction angle.\n",
    "\n",
    "To speed up the calculation, we can tell pairinteraction that this system will have some symmetries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup_system_two(system_one, angle):\n",
    "    system_two = pi.SystemTwo(system_one, system_one, cache)\n",
    "    system_two.restrictEnergy(state_two.getEnergy() - 5, state_two.getEnergy() + 5)\n",
    "    system_two.setDistance(10)\n",
    "    system_two.setAngle(angle)\n",
    "    if angle == 0:\n",
    "        system_two.setConservedMomentaUnderRotation([int(2 * state_one.getM())])\n",
    "    system_two.setConservedParityUnderInversion(pi.ODD)\n",
    "    system_two.setConservedParityUnderPermutation(pi.ODD)\n",
    "    return system_two"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can use the definitions from above to compose our calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getEnergies(bfield, angle):\n",
    "    # Set up one atom system\n",
    "    system_one = setup_system_one(bfield)\n",
    "    system_one.diagonalize(1e-3)\n",
    "\n",
    "    # Calculate Zeeman shift\n",
    "    zeemanshift = (\n",
    "        2\n",
    "        * system_one.getHamiltonian().diagonal()[\n",
    "            system_one.getBasisvectorIndex(state_one)\n",
    "        ]\n",
    "    )  # GHz\n",
    "\n",
    "    # Set up two atom system\n",
    "    system_two = setup_system_two(system_one, angle)\n",
    "    system_two.diagonalize(1e-3)\n",
    "\n",
    "    # Calculate blockade interaction\n",
    "    eigenenergies = (system_two.getHamiltonian().diagonal() - zeemanshift) * 1e3  # MHz\n",
    "    overlaps = system_two.getOverlap(state_two)\n",
    "    blockade = 1 / np.sqrt(np.sum(overlaps / eigenenergies**2))\n",
    "\n",
    "    return blockade"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a little boiler-plate, we can then calculate and plot the result with `matplotlib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAENCAYAAACy40S2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzsnWd4VNXWgN+dhBR6C4r03gVCACkCIlUQuCBKvXQsoCheAQUBAQsoAkoUUIoKgoCIICCCH0VqCj3U0Duhk97292PPYAgpk2RmzmRmv88zz2TO2WefFUKyzupCSolGo9FoNBrr4Ga0ABqNRqPROBNasWo0Go1GY0W0YtVoNBqNxopoxarRaDQajRXRilWj0Wg0GiuiFatGo9FoNFbEpopVCNFOCHFCCBEmhBiTzrpuQggphPBPdux903UnhBBtbSmnRqPRaDTWwsNWGwsh3IEAoDVwCQgSQqyRUh5NsS4fMALYm+xYdaAHUAN4CtgshKgspUy0lbwajUaj0VgDW1qsDYAwKeUZKWUcsAzonMq6ycBUICbZsc7AMillrJTyLBBm2k+j0Wg0GofGZhYrUAK4mOzzJaBh8gVCCD+glJRynRDivRTX7klxbYmUNxBCDAWGAuTJk6de1apVrSS6RqPRuAYhISE3pZS+aZwr5uHh8T1QE52TYyYJOJKQkDC4Xr16N1JbYEvFmi5CCDfgS6B/VveQUs4D5gH4+/vL4OBg6win0Wg0LoIQ4nxa5zw8PL5/8sknq/n6+t5xc3PT/W+BpKQkER4eXv3atWvfA51SW2PLJ5DLQKlkn0uajpnJh3oK2iqEOAc8A6wxJTBldK1Go9FobE9NX1/f+1qp/oubm5v09fW9h9Jfqa+x4f2DgEpCiHJCCE9UMtIa80kp5T0pZVEpZVkpZVmU67eTlDLYtK6HEMJLCFEOqAQE2lBWjUaj0TyOm1aqj2P6N0lTf9rMFSylTBBCDAc2Au7AAillqBBiEhAspVyTzrWhQojlwFEgARimM4I1Go1GkxOwaYxVSrkeWJ/i2Pg01rZI8flj4GObCafRaDQajQ3QWV4ajUaj0VgRrVg1Go1G49AkJCQwYMCAUhUrVqxRuXLl6kePHvXMyj4rV67MX7Zs2ZqlS5eu+cEHHzxpbTnNaMWq0Wg0Gofmgw8+KF6+fPnYsLCw0FdfffXGzJkzi2V2j4SEBN55553S69evP3ny5MnQX3/9tXBISIi3LeTVilWj0Wg0Dsv9+/fd1q1bV/DDDz+8AVChQoXYM2fOeGV2n61bt+YpU6ZMbPXq1eO8vb1l165db69cubKg9SU2sEGERqPRaHIOAwdS6sgRcltzz5o1iVqw4JEOfY+xZs2a/FevXvWsWrVqdYB79+65N23a9EHyNfXq1asSGRnpnvLazz777GKXLl0eAFy8eNGzRIkSceZzJUuWjNu7d29e63wnj6IVq0aj0Wgclv379/uMGTPmyqhRo8IBXnnllTK1atWKTr4mJCTkhDHSpY5WrBqNRqPJkIwsS1tx584dj3LlysUBxMfHs3379vwfffTR1eRrLLFYS5UqFXf58uWHSU+XLl16xIK1JlqxajQajcZhqVy5csyePXvyDB8+/NakSZOeaNmy5b2qVas+ohAtsVibN28eee7cOe/jx497li1bNn7VqlWFlyxZcsYWMmvFqtFoNBqHZdCgQbdbtWpVqXTp0jX9/PwiFy9efC4r++TKlYvp06dfaNeuXeXExER69ep109/fPybjKzOPVqwajUajcVh8fX0TDx48eNwae73yyiv3XnnllXvW2Cs9dLmNRqPRaDRWRCtWjUaj0WisiFasGo1Go9FYEa1YNRqNRqOxIlqxajQajUZjRbRi1Wg0Go3GimjFqtFoNBqNFdGKVaPRaDQaK2JTxSqEaCeEOCGECBNCjEnl/GtCiMNCiANCiB1CiOqm42WFENGm4weEEHNsKadGo9FoHBdrDTpPj169epVet25d3oYNG1auUKFCjYoVK9aYPHlypue+gg0VqxDCHQgA2gPVgZ5mxZmMn6WUtaSUdYBpwJfJzp2WUtYxvV6zlZwajUajcWysMeg8I/bt25e3evXqsdOnT790+vTp0KCgoGPz588vlpVh6LZsadgACJNSngEQQiwDOgNHzQuklPeTrc8DSBvKo9FoNJochnnQeWho6DFQg843bNhQILP7hIaGejVt2rRq7ty5k/Lly5d49epVz/z58yccPHjw6Llz5zzLly8fU65cufhy5crFAxQqVCipQoUK0RcuXPCsV69epnoK21KxloBHxgxdAhqmXCSEGAaMBDyBlslOlRNC7AfuA+OklP/YUFaNRqPRpMPA3weWOnLjiHUHnRerGbWg8wK7DDqvUaNGrL+/f8S77757vV27dhENGjSoMnv27AuFCxdO+uqrrwq0adPmkR7CJ06c8Dx69Gju5s2bR2T2+zK8Cb+UMgAIEEL0AsYB/YCrQGkp5S0hRD1gtRCiRgoLFyHEUGAoQOnSpe0suUaj0WhsjTUHnZ86dcqnXr160QCnT5/2fvrpp2MANm/enD/51Jx79+65de3atcJnn312sXDhwkmZldmWivUyUCrZ55KmY2mxDPgWQEoZC8Savg4RQpwGKgPByS+QUs4D5gH4+/trN7JGo9HYiIwsS1thrUHnERERIjY2Vvj6+iaGhYXlKlSoUIK3t7d88OCB2/37993Lli0bDxAbGys6dOhQoXv37rf79et3Nysy21KxBgGVhBDlUAq1B9Ar+QIhRCUp5SnTxw7AKdNxX+C2lDJRCFEeqATYZCCtRqPRaBwXaw06379/v0+FChViAA4dOuRTsWLFaIB169blM7uWk5KS6NGjR5nKlSvHTJw48XpWZbZZVrCUMgEYDmwEjgHLpZShQohJQohOpmXDhRChQogDqDhrP9PxZsAh0/GVwGtSytu2klWj0Wg0jsmgQYNuHz58OHfp0qVrHj58OPfcuXOzZDkfPHjQp1q1atEAuXPnTjpy5Eie/fv3e69fv77ACy+8cB9g06ZNeVevXl1kx44d+apWrVq9atWq1X/55ZdMJ0oJKZ3Dg+rv7y+Dg4MzXqjRaDSahwghQqSU/qmdO3jw4LnatWvftLdM9qR69erV9u/ff9zLyytTyvDgwYNFa9euXTa1c4YnL2k0Go1GYxRHjx49Zu09dUtDjUaj0WisiFasGo1Go9FYEa1YNRqNRpMWSUlJScJoIRwN079JmvWtWrFqNBqNJi2OhIeHF9DK9V+SkpJEeHh4AeBIWmt08pJGo9FoUiUhIWHwtWvXvr927VpNtCFmJgk4kpCQMDitBVqxOhBJMol5IfPoV7sfPrl8jBZHo9G4OPXq1bsBdMpwoeYR9BOIA3Hw2kFeX/c6G8I2GC2KRqPRaLKItlgdiPuxasZAZFykwZJoNI5PfDwcPw5nzsCNGxARAe7ukDcvlCunXiVLgof+K6exM/q/nAMRGa8UanRCdAYrNRrX5PRpWLkS1q+HwECIyWBKpqcnNGoEzz+vXg0bKuWr0dgSrVgdCLOlGh2vFatGY0ZKpUinT4ctW9SxOnXg9dfB3x8qVYInn4R8+SApCe7ehXPnlCV77Bhs3QoTJsD48fDUU9C3L/TrB9WqGfldaZwZrVgdCLPFGhUfZbAkGo1jsG0bvPsuhIRA6dIwZYpSjOmNXy5cGMqXh5Yt/z126xZs3gyLF8MXX8DUqfDsszBmDLRvD0IXk2isiE5eciAi4tSgeu0K1rg6t25Bnz7QooWKny5cCGFhMHZs+ko1LYoUgVdegbVr4fJl+PxzZdV26AB168Ly5coy1misgVasDoR2BWs0sGkT1KgBv/wCH34IJ05A//6QK5d19n/iCfjf/5SiXrgQYmOV0n3mGfjnH+vcQ+PaaMXqQOjkJY0rIyVMngxt2yoLMzgYJk0CHxuVdHt6KoV95AgsWqQs2WbNoGtXZc1qNFlFK1YHQlusGlclPh4GD1YJRr17q4zf2rXtc293d5XMdPKkiuH+9ZeymD//XMml0WQWrVgdiIfJSwk6eUnjOkREQOfOsGCBcv3++CPkyWN/OXLnVjHco0ehVSsYNQrq1VOWs0aTGbRidSAeJi9pi1XjIty/r+pLN26EuXOV69foDN3SpeH33+G331QS1TPPwEcfaetVYzlasToQOsaqcSWio+HFF2HfPli1CoYONVqiR+nSRcVfe/SAiROhSRPV6UmjyQibKlYhRDshxAkhRJgQYkwq518TQhwWQhwQQuwQQlRPdu5903UnhBBtbSmno6BjrBpXIS4OXnpJZeEuXqxcwY5IoUJKvuXLVdenunXhq69UIwqNJi1spliFEO5AANAeqA70TK44TfwspawlpawDTAO+NF1bHegB1ADaAd+Y9nNqdIMIjSuQmAj//a/qpjR3rip1cXS6d1fW6/PPw4gRytK+edNoqTSOii0t1gZAmJTyjJQyDlgGPPJcKqW8n+xjHsBcot0ZWCaljJVSngXCTPs5NQ8tVu0K1liRm1E3+fnwz8zaM4sVoSu4F3PPUHnGj1c1qtOmwZAhhoqSKYoXVw0mZs9WXZzq1IHt242WSuOI2LKlYQngYrLPl4CGKRcJIYYBIwFPwNyErASwJ8W1JVK5digwFKB0VtqxOBg6eUljTaLioxj3f+OYtXcWSfJf32V+r/y83fBtxjUbRy53K3VdsJAVK+CTT1Rpzf/+Z9dbWwUhYNgwaNwYXn4ZnntOJTa9/75u7q/5F8OTl6SUAVLKCsBoYFwmr50npfSXUvr7+vraRkA7opOXNNbi6oOr+M/zZ8aeGQyuO5jAwYGEvxfOjgE7aFOhDZO2T6L9kvbcib5jN5kOHlQNGRo1Ulaf0dm/2aFuXZV09corqkSoXTvVelGjAdsq1stAqWSfS5qOpcUyoEsWr3UKdPKSxhpcj7hOyx9bcvH+Rf7q8xdzX5xL/RL1KZq7KE1KN2FF9xUs7LyQ7ee303ZxW7vE9G/dUlm2BQvCr7+Cl5fNb2lz8uWDJUvgu+9gxw5V87p3r9FSaRwBWyrWIKCSEKKcEMITlYy0JvkCIUSlZB87AKdMX68BegghvIQQ5YBKQKANZTUcKaVOXtJkm/jEeLou78qFexdY12sdrSu0TnVd/zr9WdF9BcFXgumzqs8jrmJrIyUMGABXrqja0OLFbXYruyOEcmvv2qUGqjdrphKydEN/18ZmilVKmQAMBzYCx4DlUspQIcQkIUQn07LhQohQIcQBVJy1n+naUGA5cBT4ExgmpUy0layOQGxiLEkyCW8Pb+KT4klMcupvV2Mjxm8Zz66Lu5jfaT7NyjRLd23nqp2Z3mY6vx3/jVl7ZtlMpm+/VUk/U6dCAydNQaxbV422a9kSXnsNBg5Udboa10RIJ3m08vf3l8E5uPfYzaib+H7uS6n8pbh4/yIP3n9AXs+8RoulyUFsP7+d5ouaM9RvKHNfnGvRNVJKuvzShY1hG9n/6n6q+Vp3+ndoqBpG3qIFrFsHboZnddiWxESVzDR5slK2v/4K5coZLVX6CCFCpJT+RsvhTDj5f/Ocgzm+6ptHJWHpOKsmM8QlxvHaH69RtmBZZrSbYfF1QgjmdZxHPq989P+9v1VdwjEx0LMn5M+vpsc4u1IFlRk8aZKy0M+cUXHXP/80WiqNvXGB/+o5A3N8tWjuooDODNZkjhm7Z3Ds5jG+bv81uXPlztS1T+R9gpltZxJ4OZAfDvxgNZnGjoXDh9XM0yeesNq2OYKOHVXz/lKl4IUXlAWruzW5DlqxOghmi9WsWHUCk8ZSbkTeYMo/U+hcpTMdK3fM0h69avWicanGvP/3+9yPvZ/xBRmwZw/MmAGvv64UiytSsSLs3g29eqmmGF26wN27RkulsQdasToIDy1WH5PFql3BGguZsn0K0fHRTGs9Lct7CCH4qt1X3Ii8waf/fJoteWJjVfJOyZIqYcmVyZ0bfvpJ9RfesEElbx05YrRUGlujFauDYO669DDGql3BGgs4e+csc4LnMKjuICoXqZytveo9VY8eNXvwVeBXXIu4luV9pkyBY8dg3jxV6+nqCAFvvglbtsCDB2oM3fLlRkulsSVasToIKV3B2mLVWMKU7VNwd3NnfPPxVtnvoxYfEZsQm2Wr9eBB+Owz1WS/XTuriOQ0NG2qSnJq11Ydm957DxISjJZKYwu0YnUQzK5g39zKYtUxVk1GnL97nh8P/cgQvyGUyP9YK+0sUalIJfrX6c+ckDlcvp+5ZmeJiapZQuHCKr6qeZynnlKW6xtvwBdfQJs2EB5utFQaa6MVq4PwmMWqXcGaDJi2cxoCwXuN37Pqvh88+wEJSQnM2pu5phHz5qlM2FmzlHLVpI6nJwQEqBKk3btVSU5QkNFSaayJVqwOwmPlNi7qCpZSxaE06XMj8gbz98+nX+1+lCpQKuMLMkH5QuV5ucbLzAmew90Yy9JYw8Phgw9U56GcMF/VEejXD3buVPW9TZvC/PlGS6SxFlqxOggRcREIBIV91KO+q1qsU6eqmkc95zJ95oXMIzYxlncbv2uT/Uc1HsWDuAfMCZ5j0foxYyAiIudPrbE3fn4q7tq8uXKjv/qqyqrW5Gy0YnUQIuMiyeuZF59cPoBrWqyxsSo2Fx0NnTrBoUNGS+SYxCfGMyd4Dq3Lt6Zq0ao2uUfd4nVpXb41M/fMJCYhJt21u3fDggUwciRUs25HRJegSBFVivP++8qd3qwZXLpktFSa7KAVq4MQGR9JHs88+HgoxeqKyUu//KJmWs6fD7lyqZ6rmsf5/cTvXH5wmTcbvGnT+4xuMprrkdf58eCPaa5JTFSJOCVKqLmkmqzh7q4GwP/6Kxw9quKuW7caLZUmq2SoWIUQjYQQAUKIQ0KIcCHEBSHEeiHEMCFEAXsI6QpExkeSJ1cePN09cRNuLucKllIV0VevrkaMtWqlXGSax5kdOJuyBcvyQiXbtjRqWa4l9YrX4/Ndn6c5bWnhQjhwAL78EvLqmRHZpmtXCAxUyV+tWikPjpPMSXEp0lWsQogNwGDU6Ld2QHGgOjAO8AZ+TzYCTpMNIuIiyOOZByEEPh4+LucKPnlSKdLXX1cxurp14fx5NSBb8y+Hrh9i2/ltvOH/Bu5u7ja9lxCC0U1GE3Y7jN9P/P7Y+QcPYNw4aNIEune3qSguRbVqamB6p07Kvd6zp07oy2lkZLH2lVIOklKukVJekVImSCkjpJT7pJTTpZQtgF12kNPpiYxTFiuATy4fl7NYd+xQ761Nc7n9/NT7/v3GyOOoBAQG4O3hzcC6A+1yv67VulKmQBm+2vvVY+emTYPr12H6dJ2wZG3y51du4U8/hRUrlGv44EGjpdJYSrqKVUp5E0AI8aEQ4pGcfiHE0ORrNNnDHGMFlMXqgoq1aFGobOrKV7eueteK9V/uRN9h8eHF9K7VmyK5i9jlnu5u7gyrP4xt57dx6Pq/2WQXLyqF2rMnNGxoF1FcDiFUtvWWLSrjumFDldykXcOOj6XJS28Cfwohnkt27DUbyOOypLRYXS15accOVctntnyKFIEyZWDfPmPlciQWHVhEVHwUw+oPs+t9B/kNwsfDh6/3fv3w2Nixagzap9nr16+xgGbNVBy7eXNVjtOnj3YNOzqWKtbLQHvgMyGEuc2Ldv5YkeQWa+5cuV0qxnrtGoSFKcWaHD8/rVjNJMkkAoICaFKqCXWL17XrvQv7FKbP031YfHgxt6JuERKiJra8/bZ6+NHYnmLFVEnOlCmwbBn4++tyNEfG4nIbKeUFoDlQXQixAvDJ6BohRDshxAkhRJgQYkwq50cKIY6aMo7/FkKUSXYuUQhxwPRaY6mcOZWo+Kh/LVYXcwXv3KnemzR59Lifn0pqup/98aA5nj/D/uT0ndMMbzDckPu/2eBNYhJi+H7ffN59F3x9Vd2lxn64uSlPwd9/q9+Jhg3h+++1a9gRsVSxBgNIKWOklAOArYBnehcIIdyBAJSlWx3oKYSonmLZfsBfSvk0sBJIPlAyWkpZx/Ry+szjx5KXXMhi3bEDvL3/TVgyY46zHj5sf5kcjdmBsymetzhdq3U15P61nqjFc2WfY/r2ALb9k8BHH0EBXWxnCC1aKNdwkyYwZAj07asfPh0NixSrlHJIis8BUsryGVzWAAiTUp6RUsYBy4DOKfbZIqU0BxP3ACUtE9u5kFISGR9J7ly5AWWxulKMdfduqF9fNSdPTqVK6v3MGfvL5EiE3Q5jQ9gGXq33Kp7u6T7P2pQ36r1JePwFSrRcy5AhGa/X2I4nnoCNG1UTlaVLoU4d9XukcQwyqmM9bHLTpvrKYO8SwMVkny+ZjqXFIGBDss/eQohgIcQeIUSXNOQbaloTHJ6DZy/FJsaSJJMejbG6iCs4Lk49fTdo8Pi50qXV+/nz9pXJ0QgIDMDDzYOh9YYaKsedvS/C3TIUavcVHh6GiqJBdWsaP1711U5KgmefhcmTVTcsjbFkZLF2BF4EOgHupq+Tv6yCEKIP4A98nuxwGSmlP9ALmCmEqJDyOinlPCmlv5TS39fX11ri2B2zdeqKruAjR1SP4Pr1Hz/n7a2ezF1ZsUbERbDwwEJeqv4SxfMVN0yOuDj4eLIHJa8O40jE1kdKbzTG0qSJqnF95RWlaFu0cO3fGUcgozrW86bXOSA22efzUsqMfnSXgeS1ryVNxx5BCNEKGAt0klI+nOsgpbxsej+DiunaNxXSjphnsbpiHWtgoHpPTbGCyjp15T8SSw4t4V7sPZv3Bc6I+fPVz2FGX1V6MztwtqHyaB6lQAFYskRlax88CLVrq+xhjTHYsgl/EFBJCFFOCOEJ9AAeye4VQtQF5qKU6o1kxwsJIbxMXxcFmgBHbSiroZhnsZpjrN4e3i5jsQYFqZrVcuVSP1+mDJw7Z1eRHAYpJbODZlP3ybo0KtnIMDmio1WZR5Mm0O0FU+nNIVV6o3Es+vRRoZVq1VTzjj594M4do6VyPTKKsfqZX4CPEKJuimNpIqVMAIaj+gwfA5ZLKUOFEJOS9Rf+HMgLrEhRVlMNCBZCHAS2AJ9JKZ1WsaZ0BXu5exGXGGekSHYjKEhZq2m1xCtTBi5cUDEkV2Pb+W0cuXGE4Q2GIwzsGTh3Lly5opSrEKr0Jjohmvn79WRuR6R8efjnH5g4UVmttWrBpk1GS+VaZJSCMD3Z19eAL5N9lkDL9C6WUq4H1qc4Nj7Z163SuG4XUCsD2ZyGlK5gLw8vYhNjkVIa+gfV1kRGQmgodEk1NU1RtqyKwd64AU8+aTfRHILZgbMp7FOYnjV7GiZDZKTqrtSypYrdgSq9aVG2BQFBAYxsNBIPN53J5Gh4eMCECdChgyrHadNGKdoJE4yWzDXIKMb6XDqvdJWqxnLMruDkFitAfFK8YTLZg/37lSWaVnwV/u3s42px1ov3LrL6+GoG1x2MT64Me7HYjNmz1UPN5MmPHn+rwVtcuHeBtSfWGiOYxiL8/VX3srffTv/3TGNd0n3UFEKkW40upVxlXXFcE7PFao6xmmsV4xLjDK1btDVBQerdUsXqSs3e54bMJUkm8Xr91w2T4f59NcGmfXto3PjRcy9WeZEyBcowa+8s/lPtP8YIqLEIHx8111VjPzLy4awEDphe8Gh/YAloxWoFHsZYk7mCAWITYsnr6bzTo4OCoGTJ9F28rmixxiTEMC9kHi9WeZGyBcsaJsfMmXD79uPWKoCHmwfDGwznvU3vceDaAeo8Wcf+Amo0DkpGWcFdgZPA08BZ4GMp5QDTyz4DIV2AtFzBsYmxaV7jDJgTl9Ijf34oWNC1FOuK0BWER4UzvL4xfYFBKdTp01X8u1691NcM9htMnlx5mLlnpn2F02gcnIxirKullD1QzfdPA9OFEDuEEM3tIp2LkFryEiiL1Vm5c0dNtLEk7uNqtayzg2ZTpUgVWpVPNbfPLkyfrkaTTZqU9pqC3gUZUGcAS48s5VrENfsJp9E4OJbWscYA94D7qPIYb5tJ5IKYXcGpxVidleBg9W6pYnWVWtbAy4EEXg40tMQmPBxmzYKXX1alGunxVsO3iE+M59ugb+0jnEaTA8iojrWlEGIeEAI8B8wyTZvZaBfpXITI+Eg83T0fli24givYnLjk75/x2tKl4dIl28rjKMwOnE1ez7z8t/Z/DZNh6lTVFGLixIzXVipSiY6VO/Jt8LfEJMTYXDaNJieQkcW6GTWlZgfgBfxXCPGV+WVz6VyE5CPjwDVcwUFBanpNwYIZry1ZEu7eVTWVzsyNyBv8EvoL/Wr3I79XfkNkuHoVAgJUx56qVS275p1n3iE8Kpwlh5bYVjiNJoeQUVbwQFT2r8aGJB8ZB65jsTa3MFJfwjQT6fJlqFzZdjIZzXch3xGXGMew+sMMk+GTTyAhQTVzt5QWZVtQ+4nazNw7k4F1Bzp1UxONxhLSVaxSykV2ksOliYqPepi4BM4fY71wQSnJZ56xbH1J05TeS5ecV7HGJcYREBRA2wptqeZbzRAZLlyAefNgwACo8NgsqbQRQvD2M28z4PcB/H32b0OTrjQaRyCjGOt3QoiaaZzLI4QYKITobRvRXIfIeNdyBe/Yod6bNrVsvdlideY464rQFVyNuMqIhiMMk2HKFPU+blzmr+1ZsydP5HmCGXt0JwKNJqMYawAwXghxTAixQgjxjRBigRDiH2AXkA/VREKTDSLjIh+xWJ3dFbxjB+TLl3HGqZnkrmBnRErJzL0zqVKkCm0rtjVEhtOnYeFCGDr03wHzmcHLw4s36r/B+lPrOXHzhPUF1GhyEBnVsR6QUr4M1Ecp2X9Qo98GSylrSylnJZ+hqskaj8VYndxi3blTuYE9LOzdnjs3FCrkvBbr7ku7Cb4SzIiGI3ATtpzkmDaTJqmfxwcfZH2P1/xfw8vdSzeM0Lg8Fv0WSykjpJRbpZRLTU0j9COpFYmKj3rEFezMMda7d+HwYcvdwGZKlnRei3XmnpkU9C5oWInN8eOweDEMGwbFi2d9n2J5ivHf2v9l4YGFumGExqUx5vFY8wiu5ArevRukzLxiLVHCOS3WC/cusOrYKob6DX3k/4A9mThRNWofPTr7e41qMor4pHhttWpcGq1YHQBXSl7asQOsP09gAAAgAElEQVTc3aFBg8xd56wW6+zA2QAMa2BMic2hQ/DLLzBiBPj6Zn+/ioUr0r16d74J+oa7MXezv6FGkwPJlGIVQuTOeJUms0TFR7lMHevff6s2hnkzObSnRAm4fh3inMg7fi/mHnND5vJS9ZcoXSALGUNWYMIENejg3Xett+foJqN5EPdAtznUuCwWKVYhRGMhxFHguOlzbSHENzaVzEWQUj7WeclZY6x376rGEK1bZ/7akiWVC/nqVevLZRRzgudwP/Y+o5tYwQebBUJCYPVqpVQLF7bevnWL16VdxXbM3DuT6Pho622s0eQQLLVYZwBtgVsAUsqDQLOMLhJCtBNCnBBChAkhxqRyfqQQ4qgQ4pAQ4m8hRJlk5/oJIU6ZXv0slDPHEZMQg0Q+GmN1Ulfwli2QlAStstA/wNwkwlncwTEJMczcO5PW5VtTt3hdQ2T48EOlUN9+2/p7v9/0fW5E3mDB/gXW31yjcXAsdgVLKS+mOJSY3nohhDuqRKc9UB3oKYSonmLZfsBfSvk0qh52munawsAEoCGqV/EEIUQhS2XNSaScxQrgJtzwcPNwOlfw5s2QJ4/lHZeS42xNIn46+BPXIq4ZZq3u2gUbNsCoUcoVbG2eLf0sjUo2YtquaU7nedFoMsJSxXpRCNEYkEKIXEKI/wHHMrimARAmpTwjpYwDlgGdky+QUm6RUkaZPu4BTHYJbYFNUsrbUso7wCagnYWy5ihSjowz4+Xu5XQW66ZNqj+wp2fmr3UmizUxKZHPd31OveL1aFmupSEyfPghFCsGw200S10Iwfjm47lw7wIL9y+0zU00GgfFUsX6GjAMKAFcBuqYPqdHCSC5lXvJdCwtBgEbMnOtEGKoECJYCBEcHh6egTiOScoh52Y83T2d6kn//Hk4dSpr8VVQU3B8fJzDYl19fDWnbp9iVJNRhjSs37IF/u//4P33lQfBVrSt0JZGJRsx5Z8pTveQqNGkh6UNIm5KKXtLKZ+QUhaTUvaRUt6ylhBCiD6AP/B5Zq6TUs6TUvpLKf19rVErYACpuYJBxVmdyRW8dq16b98+a9cLoazWnK5YpZRM3TmVCoUq0K1aNwPur6zVp56C116z7b2EEHzU4iMu3b/E/P3zbXszjcaBSLepnBDia9IZGyelfCudyy8DpZJ9Lmk6lvIerYCxQPNk7REvAy1SXLs1PVlzKmZXcEqL1cvduRTr6tVQrRpUqZL1PUqUyPmu4D/D/iToShDzOs7D3c3d7vf/6y/VUvKbb8Db2/b3a1W+FU1LN+Xjfz5mYN2BeHvY4aYajcFkZLEGAyGAN+AHnDK96gAZRcqCgEpCiHJCCE+gB6rP8EOEEHWBuUAnKeWNZKc2Am2EEIVMSUttTMecDrMr+LEYq4fzxFjv3IGtW6FLl+ztk9MtViklE7dNpEyBMvSrY/9EdynV5JoyZWDQIPvcUwjBpBaTuPLgCvNC5tnnphqNwWQ0j/UHACHE60BTKWWC6fMcVEP+9K5NEEIMRylEd2CBlDJUCDEJCJZSrkG5fvMCK0yxpgtSyk5SyttCiMko5QwwSUp5O8vfpQOTlivYmWKs69ZBYmL2FWuJEnDliirZccuBPcP+DPuTwMuBzOs472Gtsj1ZuxaCg2H+/KwlkGWV58o9R4uyLfh0x6cM9hv82EOkRuNsWPrnqRCQPCk/r+lYukgp10spK0spK0gpPzYdG29SqkgpW5nitnVMr07Jrl0gpaxoejltWmFayUvO5ApevVrF9Pz9s7dPyZIQHw85MU/NaGs1KUnFVitWhP8a0Ot/8nOTuRZxjRm79bxWjfNjqWL9DNgvhFgkhPgB2Ad8YjuxXIeHMdbUkpecwBV8756yWLt2zb6VmZPnspqt1bHPjjXEWv31V9UXeOJEy8f1WZOmpZvSuUpnpu6cyo3IGxlfoNHkYCzNCl6IatbwG7AKaGR2E2uyh9kVnGodqxNYrCtXQkyMdawkcy1rTouzGm2tJiaqnsDVq0OPHna//UOmtppKVHwUk7ZNMk4IjcYOZMaGiAWuAneAykKIDFsaajImreQlZ4mx/vijygTOrhsYcq5i/f3E74Zaq0uXwrFj8NFHarKQUVQpWoWh9YYyN2QuJ2+dNE4QjcbGWNqEfzCwHZWI9JHpfaLtxHIdIuMj8fbwfqz0whlcwWfPwvbtylq1Rh+EYsWUYshJruCEpATe//t9qhSpwoC6A+x+//h45f6tXVu5441mQvMJeHt4M2bzY63DNRqnwVKLdQRQHzgvpXwOqAvoYYtWIOXIODPO4ApesEAp1N69rbOfu7tKgspJFuuC/Qs4fvM4n7X6DA83+wc3f/gBTp+GSZMcI5P6ibxPMKrxKH47/hv/nE+3sECjybFY+qsWI6WMARBCeEkpjwPZKPXXmEk55NxMTrdYY2Jg7lx48UVVN2ktclKTiMi4SCZsnUDjUo3pXKVzxhdYmZgY5f5t2FD9HByFkY1GUip/Kd7c8CYJSQlGi6PRWB1LFeslIURBYDWwSQjxO3DedmK5DpFxkY+V2gB4uuXsGOsvv6iymLfS682VBXJSk4gZe2ZwLeIan7f+3JCewHPmqH+rTz6xjiveWuTxzMOXbb/k4PWDzAmeY7Q4Go3VsTQr+D9SyrtSyonAh8B8Ukyq0WSNqPiotC3WHOoKlhJmzYIaNaCllYe3lCihlIVMs9GmYxAeGc60ndPoUrULjUs1tvv9HzyAjz9Ws2+t/TOwBt2qdaN1+daM+79xXI+4brQ4Go1VsTR56WEDNCnlNlODh8k2k8qFiIyPTDvGmkNdwRs3wv79MGKE9S2lkiUhMlLVxzoy4/5vHFHxUXz6/KeG3H/mTLh5UylXR0QIwdftvyYqPooxf+tEJo1zYakruJsQ4mEKihAiAMiZ42QcjLRcwTnVYpVS1UyWKQP9bFCyWco01uHixfTXGUnwlWC+2/cdbzV8i6pFq9r9/rduwRdfqBaSDRrY/fYWU6VoFUY2GsmiA4vYdXGX0eJoNFbDYsUK9BdC9DR1XkqQUtqpjbdzk1bykrmOVTq6zzMFGzZAYCCMHWubfrRly6r3c+esv7c1SJJJvLnhTYrlKcaE5hMMkWHqVOUKnjLFkNtninHNxlEqfymGrh2aYz00Gk1K0lWsQojCQojCgA8wGBgFPAA+Mh3XZJOo+KjULVZ3LwDik+LtLVKWSUiAMWOU8rOFtQqOr1h/OvgTey7tYWqrqRTwLmD3+1++DF9/DX36qBi3o5PXMy9zOs4hNDyUT/7RXVI1zkFGhXUhqHmsItl7B9NLAuVtKp0LEBkXSW6PVGKsHkqxxibEGtKtJyvMnQuHD6s2hraanlKsmJoj6oiK9V7MPUZtHsUzJZ+hb+2+hsgwZYpqYfjRR4bcPku8UOkF+jzdh092fELXal2p/WRto0XSaLJFuharlLKclLJ8infzSytVKxAZn0aM1WSx5pQ4a3i4mp7y/PO27fAjhLJaHVGxjt8ynvDIcGa3n42bsH83htOn4fvvYcgQKFfO7rfPFjPbzqSwT2EGrhmoa1s1OR5Ls4KHmepYzZ8LCSHesJ1YroGUMs1yG7OVmhNqWaWEN96AiAhVZmPrmsmyZeG8g1VR7764m68Dv+aN+m9Q76l6hsgwYQLkyqWGmec0iuQuQsALAey7uo/Pd35utDgaTbaw9LF6iJTyYQtDKeUdYIhtRHIdohOigcdnscKjrmBHZ9ky5f6dNMk+cT1Hs1hjE2IZtGYQJfOXNKy85vBh+Pln1ZCjeHFDRMg2L1V/iZeqv8SErRPYf3W/0eJoNFnGUsXqLpK1jhFCuAM5I/DnwKQ12QZyjis4LAxefx2eeQb+9z/73LNsWVVS8uCBfe6XEZ/88wnHbh5jbse55PPKZ4gM48ZB/vwwapQht7caczrMwTePL71W9Xr4+6HR5DQsVax/Ar8IIZ4XQjwPLDUd02QD8yzWtDovgWNbrJGRKp7q7q5Gk9lrgLY5M9gR3MGHrh/ikx2f0OfpPrSv1N4QGf75B9asUUq1cA7P1S+Suwg/dvmREzdPMHLjSKPF0WiyhKWKdTSwBXjd9PobVXqTLkKIdkKIE0KIMCHEY+1VhBDNhBD7hBAJQoiXUpxLFEIcML3WWChnjiIqPgpI3RXs6DHWxERV0hEaqlzBZmVnDxyl5CY+MZ6Bvw+kkHchZrSdYYgMSUnw7ruq1ePbbxsigtV5vvzzvNf4Pebtm8dvx34zWhyNJtNYZGNIKZOEEPOBHagymxNSysT0rjG5iwOA1sAlIEgIsUZKeTTZsgtAfyA1J2K0lLKOJfLlVMyurlQtVgd2BZuTlVavhq++gtat7Xt/R1GsE7dOJORqCCu7r6Ro7qKGyLB8OQQFwcKFkPvxiEKOZXLLyfx99m8Grx2MX3E/yhS04ogkjcbGWJoV3AI4BcwGvgFOCiGaZXBZAyBMSnlGShkHLCNF434p5Tkp5SEgKbOCOwNmV3CqMVYHdQUnJcGbb8K8efD+++pre+MItazbz2/n0x2fMrDOQLpV72aIDLGx6mfw9NPQ15iyWZvh6e7J0m5LSUhKoNvybsQkxBgtkkZjMZa6gqcDbaSUzaWUzYC2QEa+rxJA8o6ul0zHLMVbCBEshNgjhOiS2gIhxFDTmuDw8PBMbO0YPLRYc0gda3w8DB4MAQEqUcmoBu/mWtazZ425/92Yu/T9rS/lC5VnVvtZxgiB+jmcO6f6Aru7GyaGzahUpBI//ecnQq6G8OZ6A57gNJosYqlizSWlPGH+IKU8CeSyjUgPKSOl9Ad6ATOFEBVSLpBSzpNS+ksp/X19c95MgIcx1hxQx3rnDrzwgnI5TpgA06YZO+OzQgXVEMHeSCkZtn4Yl+9fZknXJeT1zGt/IVA/jylToG1b+7vi7UmnKp0Y++xYvt//Pd/v+95ocTQai7BUsQYLIb4XQrQwvb4DgjO45jJQKtnnkqZjFiGlvGx6PwNsBepaem1OIae4gnfuhDp1YOtWpVgnTjR+cHblynDqlHJN25Pv933Pz4d/Znzz8TQs2dC+N0/Gxx/D3bvqAcfZ+ajFR7Sp0IZh64ex59Ieo8XRaDLEUsX6OnAUeMv0Omo6lh5BQCUhRDkhhCfQA7Aou9fU2cnL9HVRoInpnk6Fo7uCExPhk0+geXNVSrNzJ/Tvb5g4j1C5MkRFwZUr9rtn8JVghm8YTpsKbRj77Fj73TgFZ8+qRvv9+6v4qrPj7ubOz11/pmT+knRe1plzd88ZLZJGky4WKVYpZayU8kspZVfTa4aUMt2/+FLKBGA4sBE4BiyXUoYKISYJIToBCCHqCyEuAd2BuUKIUNPl1VBW8kFUmc9nKbKJnYL0XMFGW6x79kCTJmr8W/fusG+fY832rFxZvZ84kf46a3E7+jYvLX+JJ/I8wZKuS3B3My6o+cEHKqY6ebJhItidIrmLsK7XOuIS4+jwcwfuxtzN+CKNxiDSLbcRQhxGldekipQy3edlKeV6YH2KY+OTfR2EchGnvG4XUCu9vZ0BsyvYJ5fPY+eMirFeuqRGvy1Zolrj/fQT9O5tvOs3JVWqqPeTJ1Xjf1uSJJPos6oPVx5cYcfAHYaV1oBqBrFsGYwfr2pXXYmqRauy6uVVtFnchu4rurO+13pyuds61UOjyTwZ1bF2tIsULkpkXCQ+Hj6pTkKxtys4Kkpll06dqlzAY8cqBZvXmNycDHnqKVW3efKk7e817v/GsSFsAwEvBNCghHFme2Ki6gVcqhSMHm2YGIbyXLnn+O7F7xjw+wCGrB3Cgs4LDJkkpNGkR7qKVUr5WNM4U8zzlpQyTUtWkzq3om7xv03/45Uar9C2Qts0R8aB/VzBUioLaPRouHgRXn5ZJcSUcfB6fCGUO9jWruBFBxbx6Y5PGeo3lNf9M0orsC3z58OBA+rn5UzNIDJL/zr9uXDvAhO2TqCgd0FmtJ2BcDSXisalycgV/AzwGXAbmAz8BBQF3IQQ/5VS6n7BmWDNiTUsOrCIRQcWMcRvCHGJcanGV+FfV7AtLdagIBgxAnbvBj8/5f599lmb3c7qVK4MISG223/buW0MXTuU58s9z+wXZhv6x/vOHeVFaNZMPfy4Oh82+5Db0beZtXcWhbwLMaHFBKNF0mgekpEPZTbwCarp/v8Bg6WUTwLNAGPmYzk4STLt+o99V/eR1zMvL1V/iSWHlxARF5Gmxeom3PBw87BJjPX6dejXTyUjnT0LCxYoJZuTlCqoOOvZsxBngzD0qVun6Lq8K+ULlWdF9xWGx/I++ghu37bPvNucgBCCL9t+Sf86/Zm4bSIz98w0WiSN5iEZKVYPKeVfUsoVwDUp5R4AKeVx24uW89h8ZjP5P83PrahbqZ7fd20fdZ+sS5vybYiKjyI0PDTVGlYz3h7eVm3lJqVyJ1atqtyJ77+vYpQDBoBbDgxTVa6s6lit3SjiyoMrtFncBoFgXa91FPIpZN0bZJKjR2H2bBgyRNUTaxRuwo3vXvyObtW68c7Gd5ix25hBCBpNSjL6c5rc/IpOcU7HWFOw5ewWIuMjuR55/bFziUmJHLh2AL/iftQsVhOA4zePp+kKBlWGY62ZlKdPw3PPqZaETz8Nhw6pGtV8xowPtQq2KLm5E32HtovbEh4ZzobeG6hQ+LGGX3ZFSuWuz5dPdVrSPIqHmwdLuy3lpeovMfKvkUzb6QIdMzQOT0ZZwbWFEPcBAfiYvsb02dumkuVAjoQfAVJPODp56yRR8VH4FfejRrEaD4+n5QoGyOeVj4j4iGzLtXq1cv0KoZrnDxqUMy3UlFSrpt6PHIEuqXaTzhyRcZF0XNqRk7dOsq7XOuqXqJ/9TbPJ77/D5s3KBVzUuCofhyaXey6WdluKh5sHozePJi4xjrHPjtUJTRrDyCgr2Albe9uOw9cPA6knHO27ug8Av+J+5PfKT6n8pbh4/2K6Fmtez7w8iH2QLZk+/1wNwK5fH1ascPxs38yQL5/qGXzwYPb3io6Ppuvyruy5tIflLy2nVflW2d80m0REqPKaGjXgdWMTkh0eDzcPfvrPT3i4efDhlg+5FXWL6W2n61IcjSFYNI9VkzERcRGcvavGraSWcLTv6j68PbypWrQqADWL1eTi/YvpxljzeuYlIi5rFquUqkPPZ5+pLNIffwQvryxt5dDUrp19xRodH03nZZ3ZfGYz8zvNN2wMXEomTFAlUDt2QC7dByFDPNw8+KHLDxT2LszMvTO5HnmdRV0WPcyw12jshX6csxJHw//tuJiaK3jftX3UfqI2Hm7qWcYcZ83IYs2qYp02TSnV116DpUudU6mCUqxhYRCZxVB0VHwUnZZ1YvOZzSzovIABdQdYV8Assn+/cv8OHapaS2osw024MbPdTD57/jOWHllKh587cD/2fsYXajRWRCtWK2F2A8PjruAkmcS+q/vwK+738FgNXxVnTTfG6pkvS4p1+XLVNalnTzWz0xniqWnx9NPKOj98OOO1KYmIi6DT0k78feZvFnZeSP86/a0uX1ZITIRXX4UiRdTDkSZzCCEY3XQ0izovYuu5rTSa34gzd84YLZbGhXDiP7n25ciNIw+/Tmmxnrh5gvux92lY4t8xY5ZarA/iMhdjDQtTmb+NG8OiRc6tVEFZrKCynDPDzaibtPyhJVvPbeWHLj/Qr04/6wuXRebMUXXFM2ZAIWMrfXI0/er0468+f3H1wVUafNeAbee2GS2SxkVw8j+79uNI+BEKeBUAHrdY917eC/BIn9nqvtXxze2bbjlHZl3BCQnKSvXwUO5fTxcILZUtC/nzZy7OeuHeBZouaMrhG4dZ9coq+tbuazP5MsuVK6q+uHVr9bPUZI/nyj1H4JBAfPP40uqnVnwT9A26G6vG1mjFagWklBy6fuihqzelxRp4OZD8XvmpUrTKw2M+uXy48u4Vetfqnea+ZsVq6R+Cb7+F4GBl8ZQunYVvJAcihHIHW6pYD10/RJMFTbgWcY2/+vxFpyqdbCtgJhkxAuLj1c9SV4tYh4qFK7Jn0B7aVmjLsPXD6PNbnyznLmg0lqAVqxXYfWk3NyJv0Lp8ayB1i7X+U/UfS/33cPNIt9Yur2deEpISLGpreOMGfPihsnS6d8/CN5GDMWcGJyamv279qfU0WdCEJJnEtv7beLaMY/Vw/OMPWLkSxo1TZUQa61HAuwBreq7h45Yfs+zIMhp81+CRhEONxppoxWoFFuxfQJ5ceej9tLI+kyvC6PhoDl0/lKVxY/k8VVskS56uJ05UmbFffeV6lk7Dhqrm82g6fye/3vs1Ly59kUqFKxE4OJDaT9a2n4AWcOeOygCuWRPee89oaZwTN+HGB89+wKa+m7gVfYt68+oxO3B2uv29NZqsoBVrNomIi+CX0F94pcYrFPEpAjzqCj5w7QAJSQlZUqx5PdUw1IwSmK5cUT2ABw1SfYBdjcaN1fuuXY+fi0mI4dW1r/LWn2/RsXJHtg/YTon8jjch/O23lddh0SLXiI0bSctyLTn42kFalmvJmxvepP2S9ly+f9losTROhFas2WR56HIi4iIY5Dco1VFvgZcDAR7JCLYUs2LNyGL98kvlBh01KtO3cArKl4cnnoCdOx89fu7uOZouaMq8ffMY02QMq15e9fDf1JFYu1Y18PjgA6hXz2hpXIMn8z7JHz3/YE6HOey4sINa39Zieehyo8XSOAk2VaxCiHZCiBNCiDAhxJhUzjcTQuwTQiQIIV5Kca6fEOKU6eU4tRDJiIqPYvL2ydQqVotGJRupmCniEYv18I3D+Ob2pXi+4pne3xLFevu2Slbq2VMpGFdECGW1JrdYN5zagN9cP8Juh/F7j9/5tNWnuLs5XofO27eVC/jpp1VsVWM/hBC86v8q+1/dT6UilXhl5Su8vOJlrjy4YrRomhyOzRSrEMIdCADaA9WBnkKI6imWXQD6Az+nuLYwMAFoCDQAJgghHK6ib8r2KZy7e+7hEGwhBF4eXo9YrCdunXjYxjCz5PPKOMb6448qturqcbnGjdUEnytXE5mwZQIdfu5A6QKlCRka4nCZv8l58024eVO7gI2kcpHK7By4kynPTWHNiTVUC6hGQGAAiUkZZMNpNGlgS4u1ARAmpTwjpYwDlgGdky+QUp6TUh7i0fF0AG2BTVLK21LKO8AmoJ0NZc00QZeD+GLXF/Sr3Y9mZZo9PO7l7vWIxXri5gmqFKmS2hYZkpHFap6v2qCBsnhcmcaNgdw36fDzC0zaPol+dfqxe9Buw8e+pcfixfDzzzB+PNSta7Q0ro2Hmwdjm43l8OuHaVCiAcM3DKfxgsYcuHbAaNE0ORBbKtYSwMVkny+ZjlntWiHEUCFEsBAiODw8PMuCZpZzd8/x4tIXKZG/BF+0+eKRc8kt1tvRtwmPCs+yxfoweSmNCTdBQWpk2qBBWdreqUh8MhBe9ePwg2189+J3LOi0AJ9cPkaLlSZnz8Ibb0DTpiq2qnEMKhWpxF99/mJJ1yWcu3sO/3n+DFs3jPBI+/190eR8cnTykpRynpTSX0rp7+vra5d7ht0Oo81PbYhJiGF9r/UUzf3okEwv938V64mbagJ38sYQmSEji3X+fPDxgR49srS9UyCl5Nugb3l+SVO8Pd0p+/dOBvsNduhZnAkJ0Lu3aje5eDG4O17o16URQtCrVi+ODzvOq/VeZW7IXCp+XZFpO6cRkxBjtHiaHIAtFetloFSyzyVNx2x9rc3YfGYzz3z/DLejb7O+93qq+VZ7bI2Xh9fDOtYTt0yK1Qau4Ph41Wy/WzfV0s8ViYqPot/qfryx/g1aV2jNmMIhnN5Rj4sXM77WSKZMgd27VdKZM83HdTYK+RQioEMAh18/TLMyzRi9eTRVZ1dl2ZFluvZVky62VKxBQCUhRDkhhCfQA1hj4bUbgTZCiEKmpKU2pmNZIj4xngX7FxAVH5Wl6+9E3+HN9W/S+qfW+ObxZc/gPTQu1TjVtZ7ung9jrMdvHieXWy7KFSqXpfvmzpUbgUhVsW7dCnfvul6XJTPn756n6YKmLD60mI9afMTanmvp9kJhAP7802Dh0mHzZpg0Cfr2dW1PQ06imm811vZcy+a+mynkU4iev/ak7ty6/HbsN61gNaliM8UqpUwAhqMU4jFguZQyVAgxSQjRCUAIUV8IcQnoDswVQoSarr0NTEYp5yBgkulYlth9aTeD1gyi3+p+mWrAfSPyBlO2T6HCVxUICArgzQZvEjI0hIqFK6Z5zSOu4FsnqFi44sMZrJnFTbiRxzNPqg0iVq2CPHlUC0NXY+u5rfh/58/pO6dZ23Mt45uPx024UaMGlCwJGzYYLWHqXLqkyqKqV1e9gDU5i+fLP0/wkGAW/2cxMQkxdF3elXrz6rH6+Grd2F/zCDaNsUop10spK0spK0gpPzYdGy+lXGP6OkhKWVJKmUdKWURKWSPZtQuklBVNr4XZkSMyTk3BXnl0Je9sfCfNRIQkmcTxm8eZFzKPDj93oOSXJflwy4c0LtWY/a/u56v2X5E7V+507+Xl4fWIxZrVxCUzqU24SUqC33+H9u1VjNWVWLh/Ia1+bEXR3EUJHBxIh8odHp4TQv2bbN4McRm3V7YrcXHw8ssQEwO//qoeijQ5D3c3d3o/3ZvQN0L5scuPRMRF8J9f/oPfPD+WHVlGQlKC0SJqHICsmVI5jOiEaACeLf0ss/bOIiAogCpFqlCqQCk83T2Jjo/meuR1wm6HPXQXly1YlhENRzDYb3Cmko/MFmtCUgKnb5+mS5Uu2ZI9NcW6dy9cvQr/+U+2ts5xTN81nf9t+h+ty7dmRfcVFPAu8NiaF16A776Dbdscy5p/7z0VV12+HKpkLeSucSA83DzoW7svPWv1ZMmhJXyy4xN6/tqTMZvH8M4z7zDIb5BDdvnS2AeXUKzmTL7vO31PbEIsy44sIzQ8lMsPLhOfGI+3hzelC5SmZdmW1HqiFk1LN6VS4bcHVt8AAB4VSURBVEpZyiz18vDiXsw9Lt2/RHxSfLpuY0vI55nvMcW6Zo2audqhQxoXORlSSsb93zg+2fEJ3at356f//ISXh1eqa9u1U8lcS5Y4jmJdtEgNRxgxwnVj4s6Kh5sH/er0o2/tvqw9sZbPd33O2xvfZuK2ibxW7zVe83+NMgV1hpqr4RKKNTpeWazeHt5ULlKZWk/Ustm9zBarufY0NasqM6RmsW7aBI0aQYHsbZ0jkFLyzsZ3mLV3FkP8hvBth2/TbU3o7Q0vvQQrVqg4ptGu8u3bVcvCVq3g88+NlUVjO9yEG52rdqZz1c7svribL3Z/wbRd05i2axodK3fkDX+VuZ5ydKTGOXGJn7LZYvX28Lb5vczlNpHxKq6bXXdQXs+8jyQv3boF+/Y5jjVmayZvn8ysvbMY0XAEczvOtajfb+/e8OCBam5vJGFhyl1fvrxS9LlyGSuPxj40KtWIX1/+lTNvnWFMkzHsubSHdkvaUfnrykzbOY2rD64aLaLGxriEYjXHWH08bG++mMttzFZmnlzZy1JJabH+/bdqZegKivWHAz8wYesE+tXux5dtv7TYNd+8OZQoAT/9ZGMB0+HOHXjxRfX1H39AwYLGyaIxhjIFy/Dx8x9z8Z2LLO22lKfyPcXozaMpOaMk7Ze0Z9mRZQ+9aRrnwiUUq10tVpMr2JyJnMcze4o1ZYx10yblAvb3z9a2Ds/eS3sZ+sdQWpZryXcvfpcpF5q7O/TvD+vWqcb89iYiQsW/T59WZVEVsxdm1+RwPN096VGzB9sHbOfE8BO83/R9Qm+E0vPXnhSfXpyha4ey88JOXbLjRLiEYo2Oj8ZduJPL3fa+OHMTfrMytIYr2LyXlEqxtmypkpeclVtRt+i2vBsl8pVg+UvLs/RzGzZMuV5nzLCBgOkQE6Pcv3v3wtKlynrWaMxULlKZKS2ncO7tc/z937/pXLUzSw4voenCppSZWYaRG0ey6+Iu3Xgih+PEf54VX38N/3crBje8GT1a9WlNSFBJLfnyqQxSX1/VWq50aSheXPVwzSrmJvzmGKu1XMFSSs6eFZw/79wDzaWUvPrHq9yIvMGewXsokrtIlvYpXlzFWhcuVJ2OChe2sqCpkJCgGkBs3qwygbt1s/09NTkTN+FGy3ItaVmuJQEvBPDbsd9YcXQFAUEBzNgzgxL5StCtWje61+hO41KNddJTDsPpFeusWXC6SgzU8OGrr5Wl5+6uLIvY2MfX+/hA7drg5wf16kGbNqqbj6WYLVZruYLzeuYlISmB2MRYtm9XruxmzTK4KAfz06Gf+PXYr0xtNRW/4n7Z2mvkSKVYv/xS9ee1JXFxqk3h6tWqtKZfP9veT+M85PXMS9/afelbuy/3Y++z9sRaVh5bydyQuXwV+BXF8hTjhUov0LFSR1pXaE1+LxdtDp6DcHrFeuIEDP4jms1nvLmYIk8gNhbu34cbN+DCBTh/Xq3fv19NHfnmG7Xu6adVzKxPH9WOLj3MFqs1k5dANeLfscObQoUyliGncjPqJiM3jqRxqca82+jdbO9Xsyb06gXTp8PgwVC2bPZlTI3ISGWdbtwIX3yhhpdrNFkhv1d+ej/dm95P9+ZB7APWnVrHmhNrWH18NYsOLCKXWy6alWlGx8od6VCpAxULV3ToSU6uitMrVnd3iE2ISTUj2MtLuYF9faFGjUfPJSXBsWOwfr1Kgpk2DT79FFq0UHM0u3RJvXzCy101LrgbcxdPd89sx3XzeeUDlGL955+iNG2aPVe1IzNq0yjuxd6zuKzGEj77DH77TbnPly+3ypaPcPs2dOyoYqrz58PAgda/h8Y1yeeVjx41e9CjZg8SkhLYfXE3f5z8gz9O/cE7G9/hnY3vUKZAGZ4v9zzPl3+eluVa8mTeJ40WW4OrJC8lRGc6I9jNTSnb995Tk2SuXlV/pM+eVT1fq1aFn39WCjg55o5At2NuZ9taBZUVDHDm8j1OnlSDsZ2RwMuBLDywkJHPjKRmsZpW27dUKXj/fVVHunKl1bYF1INXkyYQEqL21kpVYys83Dx4tsyzTG09ldA3Qjnz1hlmt5+NX3E/Vh1fRe9VvSk+vTg1v6nJiA0jWHNiDfdi7hkttsviEoo1JiEm26U2vr4werQqofj9d5X41Lu3yvo8evTfdZ7ungDcjr5tlV6h5kHq24JvAvDss9ne0uGQUjJq0yiK5SnGuGbjrL7/mDHQsCEMGmS98ptly6B+fWWx/vWX6/Vt1hhLuULlGNZgGKteWcXN924SNCSIz57/jKfyPcV3+76j87LOFJ5WmHrz6vHOn++w6tiqNIePaKyPSyjW6PhofHJZpzmEuzt06qS6H82fD6GhKtHp229VOYzZFXw7+na2E5cAiuUpBkDgkRt4e6uEKmdj3al1bDu/jQnNJzx0fVuTXLmUInRzU27b69ezvldUFAwfrrJ/69RR/w90SY3GSNzd3PF/yp/RTUfzV9+/uDP6Dlv6bWHss2PJ55mPOSFz6La8G+O3jDdaVJfB6WOsoCzW7PbsTYmbm3L9dfj/9u48TqrqSuD473RDN9gtNEuzyaqAC8ERBAQGWUJUBKQRNwJxQRMHIzHqoKODAeTDaMBIIlvQREQdVBLjQgSDCwaJ2i2yCMgOsoqyyDJAA930mT/uK7poqnqtelVWn+/n0x+qXt16dXi8rsO9775z+8GwYe6666efQvf7XGLdf2x/RIaCM9MyAVizbS+dOkFKSoV3GVcKtIBRC0fRqnYrftH+F1H7nObNC5fa69XLLYbetGnp36/q3jNiBGzZAg88ABMmWJlCE39Sq6TSs3lPejbvCcCJ/BMs3b2UmqmVoLh4nKgcPdb83KiVM6xf35WsGz/ezSSe9FRhjzUSQ8G1q9cmSZLYeWAvnTtXeHdxZ+76uaz8biWje4yOegGP7t3dZLQdO1xvc84clzCLowoLF7oSkn37utu1/vlPdwuPJVXzQ5BaJZWuTbrSpl6bkhubiKgUiTUS11iLk5QEo0a5pco2rInsUHCSJFGzal0Kqu2hU6fy7+foyaN8d6QCY6BRoKqMWzSOlrVbMvhHg335zB493O1UF1wAgwe7Yfxp09xtVidOuMlo33/vkueYMW4CW+/esGqVuz911Sob+jXGFK/SDAVH6hprcYYMgVUnUvntdsgryKN6csUTK0Bqfiak7eWKK8r+3nX71jH8neF8suMT8gvy6dioI491f4wBFw6ISGwV8e6md1n+7XJeyHqBKkn+nYotW7pFx195xc30HjGi8DWRwl5sUpKbLDZypPu3rRb9UtPGmAQQ1W8zEekDPAMkA39W1d8WeT0VeAm4HNgP3KKqW0WkObAWWO81zVbV4eWNIzcvl2rJ/nwr9u6Rym+9VVVWLU1Hb3Zf1hVRcCSTlIy9ZaoABfDhlg+54S83kJKcwsguI6mRWoOXV75M1mtZjOkxhjE9xsT05vJncp6h0bmNGNp2qO+fXaUK3Hab+9m4ET75BHbudBWUMjJcEY6OHaFO+SoqGmMqsaglVhFJBqYBVwE7gSUiMldVg25O4S7ggKq2FJHBwATgFu+1zap6WSRi8avHCoW32wBsWJ3Gn/7kFrquiKN76pHSYEWZ3rPlwBaun3M9TWs2Zd6QeTTLaAbAA10e4J559/D4oscRhDE9x1QsuHJau3ct721+j/G9xvuyOEJxWrVyP8YYEwnRvMbaCdikqltU9STwGpBVpE0W8KL3+HWgt0ShC1WeAhHlFbjdBqD5eWncfz+sW1f+/e3bB0f3ZHKqWunvQcsvyGfoG0NJkqQzkiq4pfNmDpjJHZfdwdhFY5m9cnb5g6uAqZ9PJTU5lbsvr+D/OowxJs5EM7GeB+wIer7T2xayjarmA4eAwOBbCxFZLiKLRKTcZRHyC/LJL8j3ZZFzKKy8BHBTVjrnnONqDOfnl29/S5YARzPJ5QB5p/JK9Z7JOZPJ3pnNjP4zzkiqASLCs/2fpUezHtz9zt1s3L+xfMGV08HjB3nxyxf5adufnr6dyBhjEkW8zgreDTRV1XbAg8ArInLWkg4icreIfCEiX+zdG7pH5+ci53Bmj7VhnTSmT3cl7wIF/csqJwfkmCsSse/YvhLbHz5xmCcWP8HVF1xd7EzblOQUZg+aTWpyKre+eSv5BeXM/OXwwvIXOJp3lF91smr1xpjEE83EugtoEvS8sbctZBsRqQLUBPar6glV3Q+gqkuBzUDroh+gqs+pagdV7ZCZGbrn43tiDeqxpqWkcdNNbum53/ymfBV/cnKgcW33d9t7rOTh4Kc/fZr9uft54sdPlNj2vBrnMb3fdHJ25TDhXxPKHlw5nCo4xdQlU+nWtFuFl4Uzxph4FM3EugRoJSItRCQFGAzMLdJmLhBYufJGYKGqqohkepOfEJHzgVbAlvIEkZvn1orza/JScI81PSUdEbfYem6uu9e1LFTh88+h7QVeYi2h1ufB4weZlD2JGy+5kcsbla72YWD1jLGLxrJs97KyBVgO8zfOZ8uBLdzX6b6of5YxxsRC1BKrd810BLAAd+vMX1T1KxEZJyKBmyifB+qIyCbckO8j3vbuwEoRWYGb1DRcVb8vTxwx7bF6JQ1bt4Z77oFZs2DTptLva/NmV6zgijZuKHjP0T3Ftn9h+QscOXmER7s9WqaYp/WdRr20etz65q2nj1e0TP58Mo1rNGbgRQOj+jnGGBMrUb3GqqrzVbW1ql6gqv/jbRutqnO9x8dV9SZVbamqnVR1i7f9b6raRlUvU9X2qvr38saQm+/1WP2avJR85lBwwKOPujq/48aVfl85Oe7PHh1LHgou0AKmLZlG1yZdyzzEWrt6bZ4f8Dxr9q5h3KIyBFhGa/au4YMtH/DLDr+M+S02xhgTLfE6eSli/O6xBt/HGlwruEEDuPder+zhhtLtKycH0tKgaztXL7i4oeB3N77L5gObyz3E2qdlH4ZdNoyJn0xk6TdLy7WPkkzJmUJqciq/uDx6xfaNMSbWKk1ijUWBiKKr24wc6Sr+TJ5cun3l5ECHDlC1ShJ1z6lb7FDwc8ueo0F6AwZdPKhccQNMumYS9dLqMeztYZw8dbLc+wnlQO4BXlr5EkPbDj29xqwxxiSihE+sgclLfvVYReR0ci1ahL9+fbeO56xZcPBg8fs5cQJWrOB04f3MczLDDgXvP7af+RvnM7Tt0AoNsWZUy2BG/xms2rOKJxc/We79hDJz+UyO5R3jV1fYLTbGmMSW8In1dI/Vp2usUHidNdSycffdB0ePwsyZxe9j+XJXt7ZLF/c8My0zbI/19TWvu2pLEai5O+DCAQxpO4Txi8ez8ruVFd4fFN5i071Zdy5rEJEqlcYYE7cSPrEGJi/51WOFwpnBoRY6b98eunVzBSOKWws0O9v9GVjRpnGNxuw4vCNk29mrZnNR3YsilrQm95lM7eq1Gfb2sIgUjnhnwztsPbjVbrExxlQKCZ9Y/Z68BIU91nDXdX/+c3crzSefhN9HdjY0bQqNGrnnLTJasOPQjrOufW4/tJ3F2xcztO3QiK1UU+ecOkzrO41lu5fxu09/V+H9Tf58Mk1qNCHroqKloo0xJvEkfGL1u0AEuB5rWtU0kiT04b3hBjfb96WXwu8jOxs6dy583iKjBYqy/dD2M9rNXe9qbtzc5uYKxx3sxktu5IaLb2DMP8ewdu/acu9nxbcrWPj1QkZ0GuHrmqvGGBMrCZ9YY9FjTUlOOWviUrD0dJdc58xxFZmK2r0btm0rklhrtQDg6wNfn9F23sZ5tKrditZ1zqr4WGHT+k4jPSWdO+feyamCU+Xax6TPJpGekm6r2BhjKo2ET6x+F4gANxQc6vpqsNtug8OHYW7RIo8UFoYo2mMF+PpgYWI9evIoH339Ef1a9atwzKHUT6/P5D5upZxncp4p8/t3Hd7Fq6tf5a52d5FRLSMKERpjTPxJ+MR6PP84gpxxf2m0pVZJDTkjOFivXtCwIbz++tmvZWdD1arQrl3htsY1GlMlqcoZPdYPv/6QE6dO0K91dBIrwJC2Q7iu9XWMWjiK9fvWl+m9Uz+fSoEW8Osrfh2l6IwxJv4kfGLNzXOLnEdh/fSwUpNTix0KBkhKgoEDYf78s4eDs7NdUq0WNHqdnJRM05pNz+ixztswj/SUdLo36x7J8M8QWLu1epXqDHt7WKmHhI+cPMKMpTMYdPGg08PYxhhTGSR8Yj2ef9zX66sAD3Z5kIe6PlRiu+uvh2PH4P33C7fl57vFzYOHgQNaZLQ4I7H+Y/M/uOr8q6LeG294bkOmXDuFz3Z+xu+zf1+q98xaMYuDxw/yYOcHoxqbMcbEm0qRWP2cEQww8KKBpSot2LMnZGTAG28Ublu92iXbsInVGwrednAb2w9tp1fzXhGKunhD2g4h68IsHlv4GOv2rSu2bX5BPn/I/gNdGnehS5MuvsRnjDHxIuETa25+ru891tKqWhWuuw7+/nfIy3PbAoUhQibWWi3Ye2wvR04e4V/b/wVAt6bdfIlVRJjRfwZpKWnc8dYdxQ4Jv7rqVTYf2FyqXrsxxiSahE+ssRgKLousLLfmamAmcHY21KsHzZuf3TYwM3jrwa0s3r6YGqk1uLT+pb7F2iC9AVOvnUrOrhye/uzpkG3yC/J5fNHjtGvQztZcNcZUSgmfWHPzc3291aasevd2E5kWLHDPA4UhQs21Cr6XdfH2xXRt0pXkpGQfo4XBPxrMoIsHMfqj0SELR7z85ctsPrCZx3s+7uuEMWOMiRcJn1jjvceakeHqAS9YAPv3w/r1oYeBAS6scyEpySlM/HQia/auoVsTf4aBg4kI0/tOJy0ljeHzhqNBBY/zTuUx7uNxdGjUgf6t+/semzHGxIOET6y5ebm+T14qq2uugS++gFGjCp+HUqt6Lab3nX76+uqVza70KcIz1U+vz8SfTOTjbR/z4pcvnt4+a8Usth7cyrie46y3aoyptBI+scZ7jxVcIlWFZ5+Fn/3MrYATzl3t72Jkl5FknpNJx0Yd/QuyiGHthvHvTf6dke+NZP+x/ZzIP8H4xePp3LgzfVr2iVlcxhgTa1FNrCLSR0TWi8gmEXkkxOupIjLHez1HRJoHvfaot329iITpw5Us3q+xAnTsCLVquRrCEyaU3P6pq5/im//8JqY98SRJ4o/9/sihE4d4+P2HeWDBA2w/tN16q8aYSi9qy42ISDIwDbgK2AksEZG5qromqNldwAFVbSkig4EJwC0icgkwGGgDNAI+EJHWqlrmSvA/hB5rcjJMmeJWvAksE1eSeFgppm39tjzY+UEmfjoRgIe6PsRVF1wV46iMMSa2ovnt3AnYpKpbAETkNSALCE6sWcBY7/HrwFRx3Z0s4DVVPQF8LSKbvP19VtYglt297AfRgxo6NNYRlM/oHqOZu2Eubeu15cneT8Y6HGOMibloJtbzgB1Bz3cCV4Rro6r5InIIqONtzy7y3vOKfoCI3A0E1iM7IiLhqsTXBfaV9S/gox98fOtYx1/5q0/hhPSDP4YxZvFVXLzHGC6+Zn4HkuhiP55YAar6HPBcSe1E5AtV7eBDSOVi8VVcvMdo8VVMvMcH8R9jvMeXSKI5eWkX0CToeWNvW8g2IlIFqAnsL+V7jTHGmLgTzcS6BGglIi1EJAU3Ganost5zgdu9xzcCC9VVHJgLDPZmDbcAWgGfRzFWY4wxJiKiNhTsXTMdASwAkoGZqvqViIwDvlDVucDzwMve5KTvcckXr91fcBOd8oF7yzMjOEiJw8UxZvFVXLzHaPFVTLzHB/EfY7zHlzAkuCSdMcYYYyom4SsvGWOMMX6yxGqMMcZEUEIl1oqUUPQhtiYi8pGIrBGRr0Tk1yHa9BSRQyKywvsZ7Vd83udvFZFV3md/EeJ1EZHJ3vFbKSLFVDWOeGwXBh2XFSJyWETuL9LG9+MnIjNFZI+IrA7aVltE3heRjd6ftcK893avzUYRuT1UmyjF95SIrPP+Dd8UkYww7y32fIhifGNFZFfQv2PfMO8t9vc9yjHOCYpvq4isCPPeqB7DcN8r8XQOVkqqmhA/uAlSm4HzgRTgS+CSIm1+CczwHg8G5vgYX0Ogvff4XGBDiPh6Au/E8BhuBeoW83pf4F1AgM5ATgz/rb8FmsX6+AHdgfbA6qBtE4FHvMePABNCvK82sMX7s5b3uJZP8V0NVPEeTwgVX2nOhyjGNxYYWYpzoNjf92jGWOT1p4HRsTiG4b5X4ukcrIw/idRjPV1CUVVPAoESisGygMA6Z68DvcWneoequltVl3mP/w9YS4hqUnEuC3hJnWwgQ0QaxiCO3sBmVd0Wg88+g6p+jJvRHiz4PHsRGBjirdcA76vq96p6AHgfiPiyQKHiU9X3VDXfe5qNu088JsIcv9Ioze97RBQXo/f9cTPwajQ+uyTFfK/EzTlYGSVSYg1VQrFo4jqjhCIQKKHoK28Iuh2QE+LlLiLypYi8KyJtfA0MFHhPRJaKKxdZVGmOsR8GE/6LLJbHL6C+qu72Hn8L1A/RJl6O5Z24UYhQSjofommEN1Q9M8wwZrwcvyuB71R1Y5jXfTuGRb5XfkjnYMJJpMT6gyAi6cDfgPtV9XCRl5fhhjf/DZgCvOVzeN1UtT1wLXCviHT3+fNLJK7YyAAIWZg41sfvLOrG3OLynjYRGYW7T3x2mCaxOh/+CFwAXAbsxg21xqufUnxv1ZdjWNz3Sjyfg4kqkRJrRUoo+kJEquJO/tmq+kbR11X1sKoe8R7PB6qKSF2/4lPVXd6fe4A3ccNtweKh1OS1wDJV/a7oC7E+fkG+CwyRe3/uCdEmpsdSRO4A+gNDvS/es5TifIgKVf1OVU+pagHwpzCfG/Nz0fsOGQTMCdfGj2MY5nsl7s/BRJZIibUiJRSjzrsW8zywVlUnhWnTIHDNV0Q64f59fEn8IpImIucGHuMmuKwu0mwucJs4nYFDQcNNfgnbQ4jl8Ssi+Dy7HXg7RJsFwNUiUssb6rza2xZ1ItIHeBgYoKrHwrQpzfkQrfiCr9tfH+ZzS/P7Hm0/Adap6s5QL/pxDIv5XonrczDhxXr2VCR/cLNWN+BmC47yto3DfYEAVMMNIW7C1R4+38fYuuGGY1YCK7yfvsBwYLjXZgTwFW6GYzbQ1cf4zvc+90svhsDxC45PcIvXbwZWAR18/vdNwyXKmkHbYnr8cEl+N5CHu0Z1F+66/YfARuADoLbXtgPw56D33umdi5uAYT7Gtwl3bS1wHgZmyjcC5hd3PvgU38ve+bUSlyAaFo3Pe37W77tfMXrbZwXOvaC2vh7DYr5X4uYcrIw/VtLQGGOMiaBEGgo2xhhjYs4SqzHGGBNBlliNMcaYCLLEaowxxkSQJVZjjDEmgiyxGmOMMRFkidUYY4yJIEusxnhE5D9EZLe3duYmEXnLq+pTtF11EVkkIskiUl9EXhGRLV6h9c9E5HofYk0RkY+9snrGmDhiidWYQm2B/1bVy4DWwI+AS0O0uxN4AyjAFfr/WFXPV9XLcaX1or4Mm7ql0j4Ebon2ZxljysYSqzGFLgWWe49b4ko4bgjRbiiu9uqPgZOqOiPwgqpuU9Upgeder3epiHwVWDZMRJqLyOqgNiNFZKz3OE1E5nlL360WkVtCbfPe+pYXizEmjtgwkjGF2gAveauFnAf01yJLcHlDw+er6lYRGYBbqq44d6rq9yJSHVgiIn8roX0f4BtV7ed9Xs0w28AVdO9Yhr+fMcYH1mM1BhCRJsAeVb1UVS/GFfT/TYimdYGDYfYxzetVLgnafJ+IBBYFaAK0KiGUVcBVIjJBRK5U1UNhtqGqp4CTgRVUjDHxwRKrMU5bYE3Q8y+BeiHa5eJWSQK3Ykn7wAuqei/QG8gEEJGeuKXFuqhbfH259958zvzdqxa0jw3ePlcB40VkdKhtQe9NBY6X8e9qjIkiS6zGOJcCa+H0Gpe345bbOoOqHgCSRaQasBCoJiL3BDU5J+hxTeCAqh4TkYuAzt7274B6IlJHRFJxC47jfXYj4Jiq/i/wFNA+1DavbR1gn6rmVfyvb4yJFLvGaozTFughIv1ws31zgJFh2r4HdFPVD0RkIPB7EXkY2AscBf7La/cPYLiIrAXW44aDUdU8ERmHWxN4F7CuSBxPiUgBbv3Pe8JsA+gFzKvw39wYE1G2HqsxZSQi7YEHVPXWGMfxBvCIN1RsjIkTNhRsTBmp6jLgIxFJjlUM3uzktyypGhN/rMdqjDHGRJD1WI0xxpgIssRqjDHGRJAlVmOMMSaCLLEaY4wxEWSJ1RhjjIkgS6zGGGNMBP0/qtFluyEsKhMAAAAASUVORK5CYII=",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1126bcb70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xlabel(r\"$B$ (Gauss)\")\n",
    "plt.ylabel(r\"Blockade (MHz)\")\n",
    "plt.xlim(-0.4, 20.4)\n",
    "plt.ylim(0, 0.4)\n",
    "\n",
    "if sys.platform != \"win32\":\n",
    "    with Pool() as pool:\n",
    "        energies1 = pool.map(partial(getEnergies, angle=0), bfields)\n",
    "        energies2 = pool.map(partial(getEnergies, angle=np.pi / 2), bfields)\n",
    "else:\n",
    "    energies1 = list(map(partial(getEnergies, angle=0), bfields))\n",
    "    energies2 = list(map(partial(getEnergies, angle=np.pi / 2), bfields))\n",
    "\n",
    "plt.plot(bfields, energies1, \"b-\", label=r\"$\\theta = 0$\")\n",
    "plt.plot(bfields, energies2, \"g-\", label=r\"$\\theta = \\pi/2$\")\n",
    "plt.legend(loc=2, bbox_to_anchor=(1.02, 1), borderaxespad=0);"
   ]
  }
 ],
 "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.10.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}