{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rydberg Pair Potentials Near Surfaces"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial is based on results that were published by J. Block and S. Scheel, \"van der Waals interaction potential between Rydberg atoms near surfaces\" [Phys. Rev. A 96, 062509 (2017)](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.96.062509). We will reproduce the pair potentials shown in Figure 4. The final result is that for states around the $|70p_{3/2};70p_{3/2}\\rangle$-asymptote of Rubidium the strength of the pair interaction is reduced when placing the atoms in front of a perfect mirror (perfectly conducting plate) compared to the vacuum interaction."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described in the [introduction](introduction.ipynb), we start our code with some preparations and load the necessary modules."
   ]
  },
  {
   "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",
    "from itertools import product\n",
    "\n",
    "# Operating system interfaces\n",
    "import os\n",
    "import sys\n",
    "\n",
    "# Parallel computing\n",
    "from multiprocessing import Pool\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": [
    "The plate lies in the $xy$-plane with the surface at $z = 0$. The atoms lie in the $xz$-plane with $z>0$.\n",
    "\n",
    "We can set the angle between the interatomic axis and the z-axis `theta` and the center of mass distance from the surface `distance_surface`. `distance_atom` defines the interatomic distances for which the pair potential is plotted. The units of the respective quantities are given as comments.\n",
    "\n",
    "Be careful: `theta = np.pi/2` corresponds to horizontal alignment of the two atoms with respect to the surface. For different angles, large interatomic distances `distance_atom` might lead to one of the atoms being placed inside the plate. Make sure that `distance_surface` is larger than `distance_atom*np.cos(theta)/2`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.pi / 2  # rad\n",
    "distance_atom = np.linspace(6, 1.5, 50)  # µm\n",
    "distance_surface = 1  # µm"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define the state that we are interested in using pairinteraction's `StateOne` class . As shown in Figures 4 and 5 of [Phys. Rev. A 96, 062509 (2017)](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.96.062509) we expect large changes for the $C_6$ coefficient of the $|69s_{1/2},m_j=1/2;72s_{1/2},m_j=1/2\\rangle$ pair state, so this provides a good example. \n",
    "\n",
    "We set up the one-atom system using restrictions of energy, main quantum number n and angular momentum l. This is done by means of the `restrict...` functions in `SystemOne`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "state_one1 = pi.StateOne(\"Rb\", 69, 0, 0.5, 0.5)\n",
    "state_one2 = pi.StateOne(\"Rb\", 72, 0, 0.5, 0.5)\n",
    "\n",
    "# Set up one-atom system\n",
    "system_one = pi.SystemOne(state_one1.getSpecies(), cache)\n",
    "system_one.restrictEnergy(\n",
    "    min(state_one1.getEnergy(), state_one2.getEnergy()) - 50,\n",
    "    max(state_one1.getEnergy(), state_one2.getEnergy()) + 50,\n",
    ")\n",
    "system_one.restrictN(\n",
    "    min(state_one1.getN(), state_one2.getN()) - 2,\n",
    "    max(state_one1.getN(), state_one2.getN()) + 2,\n",
    ")\n",
    "system_one.restrictL(\n",
    "    min(state_one1.getL(), state_one2.getL()) - 2,\n",
    "    max(state_one1.getL(), state_one2.getL()) + 2,\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pair state `state_two` is created from the one atom states `state_one1` and `state_one2` using the `StateTwo` class.\n",
    "\n",
    "From the previously set up system_one we define system_two using `SystemTwo` class. This class also contains methods `set..` to set angle, distance, surface distance and to `enableGreenTensor` in order implement a surface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Set up pair state\n",
    "state_two = pi.StateTwo(state_one1, state_one2)\n",
    "\n",
    "\n",
    "# Set up two-atom system\n",
    "system_two = pi.SystemTwo(system_one, system_one, cache)\n",
    "system_two.restrictEnergy(state_two.getEnergy() - 5, state_two.getEnergy() + 5)\n",
    "\n",
    "system_two.setAngle(theta)\n",
    "system_two.enableGreenTensor(True)\n",
    "system_two.setDistance(distance_atom[0])\n",
    "system_two.setSurfaceDistance(distance_surface)\n",
    "system_two.buildInteraction()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we diagonalize the system for the given interatomic distances in distance_atom and compare the free space system to a system at `distance_surface` away from the perfect mirror. The energy is calculated with respect to a Rubidium $|70p_{3/2},m_j=3/2;70p_{3/2},m_j=3/2\\rangle$ two atom state, defined in `energyzero`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Diagonalize the two-atom system for different surface and interatomic distances\n",
    "\n",
    "\n",
    "def getDiagonalizedSystems(distances):\n",
    "    system_two.setSurfaceDistance(distances[0])\n",
    "    system_two.setDistance(distances[1])\n",
    "    system_two.diagonalize(1e-3)\n",
    "    return system_two.getHamiltonian().diagonal()\n",
    "\n",
    "\n",
    "if sys.platform != \"win32\":\n",
    "    with Pool() as pool:\n",
    "        energies = pool.map(\n",
    "            getDiagonalizedSystems, product([1e12, distance_surface], distance_atom)\n",
    "        )\n",
    "else:\n",
    "    energies = list(\n",
    "        map(getDiagonalizedSystems, product([1e12, distance_surface], distance_atom))\n",
    "    )\n",
    "\n",
    "energyzero = pi.StateTwo(\n",
    "    [\"Rb\", \"Rb\"], [70, 70], [1, 1], [1.5, 1.5], [1.5, 1.5]\n",
    ").getEnergy()\n",
    "y = np.array(energies).reshape(2, -1) - energyzero\n",
    "x = np.repeat(distance_atom, system_two.getNumBasisvectors())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEOCAYAAACjJpHCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzsvXl4XMWV9/+pbqlt2UgIGU8MGCRAHscGg4wdQWeR5RjkEEgwOMubyCHJhMjOm5cXZomM4ZfJMAwyS+bHJIQk0hMCYTCQCQaSycDYwRNBkhY2O2aNjRcwSzA2xmy2Fp/3j3u7feuqb3Xr9qKWVJ/nuc/tc+su1S11fbvqnDqlRASLxWKxWMISGe4KWCwWi2VkY4XEYrFYLDlhhcRisVgsOWGFxGKxWCw5YYXEYrFYLDlhhcRisVgsOTEsQqKU+rxS6hml1AGl1FzDedVKqTuVUs8rpZ5TSsWLWU+LxWKxZGa4eiRPA+cBD2Y47wfAf4vIh4GTgecKXTGLxWKxDI2y4XioiDwHoJQKPEcpdSjQBHzNvaYX6C1C9SwWi8UyBIZFSLLkWGAncJNS6mTgUeAiEXkv3clKqTagDWDixIlzPvzhD2f9oJ07d7Jnzx6qq6uZPHlyzhXfunUru3fvTtk1NTUce+yxOd/XYrFYCsWjjz76poiEagBVoVKkKKXuB6akKbpMRH7tntMN/IOIPJLm+rnAQ8DHRGS9UuoHwF4R+W6mZ8+dO1ceeWTQLdPS1dXF0qVLU3ZnZydtbW1ZXRvEIYccwnvvHdS7iRMn8u677+Z0T4vFYikkSqlHRSTQZ22iYD4SETldRE5Ms/06y1vsAHaIyHrXvhM4Jd/1XL16tdEOQyQSMdoWi8UymijZFk5EXgdeVkpNdw8tAJ7N93MaGhqMdhimTp1qtFmyBCZNcvYWi8Uywhmu8N9zlVI7gDjwX0qpNe7xI5VS93pOvRBYpZR6CmgAOvJdl+rq6pTTXylFdXV1zve8+OKLg+0lS2DVKti929lbMbFYLCOcgvlIhpOh+Eh6enpYsGABvb29xGIx1q1bRzye+3SVrq4uVq9ezeLFi3Wfy7hx0OsJPovFYP/+nJ9nsVgsuZCLj2TMCwk4YtLd3U1zc3NeRMRILIb09aEAAVR5uS4sFovFMgzkIiSlHP5bNOLxeOEFxOUPRx7Jx7dvRzz2J4ryZIvFYikMJetsz4VNmzbR1dU13NVIy+mvvca/A28C/+7aFovFMpIZlT2SvXv3puaG5DonJN+ICF/12OWjcGjRYrGMLUZljyRJPuaE5JsvfOELRttisVhGGqOyR5Jk8eLFw12FQdx6660A3HfffZx55pkp22KxWEYqozJq69BDD5Vrr7225Ia1LBaLpVQpyRQpw8m0adOGJiJdXbBwobO3WCwWy5AY1UNbWdHVBcmkjWvXOnvbk7FYLJasGZU9kiHhd8iXoIPeYrFYShkrJH6HfDEc9HV1EIk4e4vFYsknpqH6oLanro45MCfsI62QtLVBezvU1zt7/7BWvv0ndXXI9u2ICLJ9uxUTi2W00NMDK1c6ey8zZ0I06uz9GH5U7pk0iQNKsWfSpEFlr0+cyIBSvD5xol7Q1YUsXYqsXYssXaq3W0Ftj3s8J0Rk1G1z5syRrEkkRCoqRKJRZ59IHCzr7BSBg1tnZ/b3DWAA5IB7vwMgA5DzPS0WyxDp7BRpaUn/na6tFVHK2WdblkjIAfc7fQAOtiMzZujHZ8zQ7qWVee75Vk2NVvZWTU2q7LUJE7Sy1yZMSJW9HYlo7cvbkUiqLKjtSR6f49ih2lzbI+nu5sC+fTAwwIH9+6G7+2BZAfwnW9y9+GyLxRIC09o+Qb/2kwE2a9c6+2x+tWcoe6+lBQDlsweee047nrQBDri9AOWzASrdpbqVzwaY/P77WlnSBqg4cEB7q147qO3JRxs05oXknj172CdCH7DvwAHu2bPnYKF/kas8LHp1SmUlLwIHgBdd22KxMPSxfdPaPiZB+M53Uo2puHYSU+NuKhvnW0o7ab/hexte2/SjcpuvbFuWZY+5ayuJz4bgtuejkyfzIjkStitTyttQhrZaWlrkNJBLQE4DaWlpOVi4aJHWFZRFi7K+bxCdnZ3i/p0FkM48DJdZLCOGxkaRsjJn78U0jBw0BBSL6d/PWCx1iWkIuc9X1ucp26KU9qwtSqXKNnmOHwDZ5Lluo69so1u2YMIE6XOf3weywDMMVVlZKZtA+t17VVZWpsqmTJmilU2ZMiVV1traqpW1trZ6PsZO+RPIPpA/+dqXoLYnkUikjknINnfYG/1CbFVVVVk30KaG/fmyMu0f7vmysqzumc0zW1paBtVxayQi/SBbPeOaFktJEsbH0Nioi4VXTKqqdFGoqkoVBYnCPt/xfVkKQsLX6Cc818ViMa2RjnnEKVPD/6QrFk96Gv7Ozk7th2o2DXumMhFHTGpqajQR8V6brn0xlbliskOskHje1BB/7Qd9uKt9/6irPf9w+War6yRLblZMLMNOkFhk6D1oZV4x8TmCxfM/3u/7rvV7vmvbolHtu7EtGhURkVt8gnBLloIwefJk7Vf75MmTU2VTpkzRGnBvTyBswx+mYc9UVgiAR8QKyWAhmTlzZk4fbFN5udYtbSovz+l+JkxfJIuloKQTDJNYVFXpZZ7eg7g9gZRYeHoC+8aN03sQ48alyjb4xGKDKxYiwaJQW1srvwB5A+QXILUe0TIJwrJly7SyZcuWed62uSdQSg1/vrFCEiAkM7yhdgaC/gG6W1u1f+7uNN3IfGF7JJaCMtTehUksDD0LU9jq31dWamV/7xkamjBhgtZLmODxJcyYMUNr3JPf65aWFu24179pEoREIiGxWEyUUhKLxSThDfmXkS8IYbFCEiAk2fwjmP7hnpwyRfvCPOn5VVMIrI/EkhOJhEhHhz4XSiRc78IgFq/6nN+v+noCQU7iCRMmyAUg94Fc4BOLxsZG7XvY6PGfBAlGLr2HRCIhHR0dg0RkLGOFJEBI0jmi/LS0tGj/3N5fNRf5Jv5c5PnHt1hKikTCaeyTjb63gczQu9DKXMHIJBbe3oNfLLyNe65iIRLe/2AZGiNOSIDPA8/ghDTPNZz3t+55TwO3A+OzvL8AUuPpVgdhGr5qbW3VRCYbYbJYCkpQ+Gx9vS4I9fUHywLEInk/LbTWvW8piUWy3ApGYRmJQjIDmA50BwkJcBSwFahw7f8Avpbl/bPukSS/SOL7IiVpaWmRiooKfX6JxVJIgvwZpvDZaFQv8zirg8RCxPmx5BWM5HfGisXYY8QJSerhmYXkZaAGZ92U3wItWd43ax/JloYGTUi2NDSkyhKJhIwbN06UUjJu3Dg7nmopPCZ/hqlnUVOjl3l646ZJahMnTtQa94kTJ4qIWRCsWIxORqWQuOUXAe8CO4FVGe7VBjzibgLZhf9e1Ngo+3BCfPeBXOT5dWUKE7RYciKo1zFlii4I3gCPykq9zBP1JJ2deq/Dc19TKGxQWaFCZC2lSy5CUrAVEpVS9wNT0hRdJiK/zuL6w4BzgGOBPcCvlFJLROTWdOeLSBfQ5V4r7rGM9Zz5jW/QvGEDzTiq9vVvfCPjNRZLTphW5fQk5xtkf//7B69L2slbAg8Di4HVwEdwflkB7N27V7ul116wYAGrVq3S7Ewkl7FevXo1ixcvHrSsdVtb29CWuraMfMIqUD42zENbnwdu9NjnAz/O8r4CSHt7e1ZKbEobUF5eLoCUl5cXfGirtrZWlFLaxCrLCKe93XF8e/8XTb0OU5k4wSEP19QMmtNk6nW0trZqZV7fYdAwVaa5FpbRB6NxaAs4FSdiawJOws1fABdmed/BQ1FBMfYGivllqq2t1eLvrZiMAtrbdVFIikksph/3pO8w+UhMQ0om53jY1B52rsXYYsQJCXAusAPYD/wFWOMePxK413Pe5cDzOOG//w6My/L+AkhTU5PzCSUSIuPGOekaxo3LWkw6OjokGo0KINFoVDo6OrK6LgymzKIitrdS0gT5O6qrdVGornaOF6DXEdY57lTf+jQsI1BICr0lvzB1dXXOJ7Rsmf7FzdJpnkgkpKKiQqLRqFRUVBT0l5kp11Ztba3WEFgxKSFMUVYTJuhlyV5CAXoduTjHLRYRKySBQpL6xbZokf7F9a0rYkrJ3N7eLvX19Vn7W0JjWHZTKaUNeylPIjxLkQgTZWXIhLu5vV1erK+Xzb7/q7C9jkwRhrbXYcmEFZIAIUklbTT0SEyOyKL/kgtYx8EmdBxmTL2OEP6ORCIhZWVlAkhZWZnW0w3b67DOcUuuWCEJEJJUw59IiJSXO2+3vFzzkdTU1GhfTm9alUxjy0XDkJrbUgRyiLJK15NpamrS/q9Svjwx/7DJ1OuwznFLLuQiJKNyzfZoNEpra6sey67Uwc1DfX09pwGXAKe5dpLFixdr5/rtonHMMal1opVrWwpA0Jrhb70VbF9+uV7ms3tmzWJlczM9s2aljm3cuFE7x2s3NTVpZV77/PPPJxaLoZQiFotx/vnna+fG43FWrFhBPB5P+/YsloIRVoFKecPfI+noOJiLKBp1bJe/mTFDPnD9Dx+A/I1vDZOSGVsOWr7Ukh9Mw1fTpull06YNvjaN/yRoCMsUPJGpF2x7HZZCgR3aSi8kqS9hIiFSUeGISEWFNrR1zxFHaMNG9xxxRJi/gWWkYxqiChn1FzSE1d7erh33BnLYCCvLcJGLkIzKoa0kqaGoeBzWrYMrrnD2nq5/dXW1do3ftowyli+HadOcvRdTapLzz4dYzBkWjcUc20NPTw8rV66kp6dHOx40hGVKWdLW1kZnZyctLS10dnbaVCOWkUFYBSrlraqqKutfck91dmpJG5+yvwBHL0EzzUUyO80DMiOYIrCChrBsMlBLKYLtkehMmzYt619yL774otG2jCL8TnSvncFpTjwOK1ZovVmASy65hP7+fgD6+/u55JJLUmVf/OIXtXOTdianucUy0hiVQvLEE0+wZMmSrM7d19VFDOeDiLl2KdLV1cXChQvpKtH6lRRB0Ve9vcF2Wxt0dkJLi7PP8oeIKQIraAgrHo/T3d3NlVdeSXd3t42ysox8wnZlSnkjTQx+ELdXV2vO9tuT+ZBKiEwOWJuHy4Mp+sow0zwTQRkOTBFYdgjLMpLARm2lF5KYd5ZxAHe3t0uv6yPpBbnbnwolRNbgfFNVVaU1SFVVVakym4fLh8nXYRIZA6YoK1OZnW1uGUnkIiQFW9iqFHA+G4e6ujpeeukljjnmGLZt25Y6vuj440meVebaKXp6YP58ZwgkFoPf/37QGHkxePfddwPtl156SSvz26Oanh7o7obm5oN/F9PkweRw1erVsHhx1sNXt9122yD76quvBswRWMkhrO7ubpqbm+0QlmX0ElaBSnnDN7Rl+tX+5JQp2tDWk775A1pqkmEammhsbNSSNjZ6lgMes+uYJBIiZWVOz6Ks7GCPMdPkQeMt00/2mzZtmvb/M81zTzt8ZRktYKO2dJIpUm691VmV1/Sr/efuL0jx2QBbHnpIu85vF4v1f/kLx+MEBBzv2km2uceSZduKX73h4ZJLwI2Wor/fsQH8S8VmsXQsOHNB5s2bx6WXXsq8efO0OSGHHXaYdq7XthFYFssojdpqaGhIiQjAMb7cVF77zXPPpQ1Yg7PG9Zvnnpsqe/yVV7Tr/HbReOklFKQ2vMJoKhvN+KKlUnaGyYNBXHPNNfT19QHQ19fHNddckyprbm7WzvXaNgLLYmF0+0iSbNu2LdBHcuuttzLx7rv52fvvM2HCBN7zCNDDM2Zw1s6dlAN9rj0saRuPOQa2b9ftbMpGA+n8IDAo+WbKjsed89NdY+CFF14ItE1+EOeRcSsgljHNmBASQBMPL3V1dbz//vsAvP/++9TV1aXOPeeqq/jkRz/KPOAB4F+vuqoodR3Etm1QV+f0No45xrGzKRvp9PRAU5MzdFVWBg8+eFAYPvQhPY3Jhz508HU8PuSgiHHjxhlti8USzKgc2hoKJv/JDTfckPKdiGsPG9u2wYED6YXCVDaSCfKDAFx8sX6u3w4gKC9WLBYLtK0fxGIxMyp7JJs2baKrqyurNCnHHHMM2z1DQ17/yc7f/IZ1ODPee4Fzf/Ob/FfWEkyQHwRChfL29PTQ1NREf38/ZWVlPPjgg6khqebmZjZs2JA6N50fxIbxWizpGZU9kr1797J06dKs0olceumlgfY5hx5KDEdty117pLFkyRImTZqUdcqYYaGnB1audPZegvwgSdraYM2arOeDmPJiZeMHsYtGWSzpGZVCkuQHP/jBkM/x2nft3k0vjqO9z7VLkaA8XEuWLGHVqlXs3r2bVatWlaaYJP0gl17q7L1i4vV7pLOHiCkvlsViyYGwE1By2YBrgeeBp4C7geqA8z4FvABsBi7J9v4ngHSAzPCtdpiOGTNmaBPKvNdEIhE5DeQSkNNAIpFI5lk9RcaUhysWi2ll2aSMKTpNTfoEQs/65WFTmgRNLDRNTLXpTCxjHUZari2gBShzX18NXJ3mnCjwInAcjpviSWBmNvef485Ef8y3TGk6TA2xSWRKBVMervLycq2svLx8GGsawGGH6WJx2GF6ecAytkEkEonU+y4vL9cEwZQXK3mtXcbWMlbJRUiGZWhLRNaKiBuOw0PA1DSnNQKbRWSLiPQCdwDnDOU5s594IuM5phXp5s2bp53rt+vq6ohEItTV1Q2lWnnFlIfrC1/4ApuAfmCTaw8bQX4Qf5it3x6iH8Q0sTDTapjWD2KxhKMUorb+BvhlmuNHAS977B3AqUE3UUq14UxOZw7uLO9IdjrZ1tY25CVN6+rqUtFe27dv1+afFJO5c+dq0UZz585Nvb71j39MhS8f79rDQk8PzJsHfX1QXg4PPHBwnsdf/zW8/vrBc//6r3N61KuvvhpoT5o0SSvz2xaLJRwF65Eope5XSj2dZjvHc85lOD+YV+X6PBHpEpG5InKwJZ0yJad7nn/++USjUcDJ3+WdP1AqWXfXr1/Py+Xl9AMvl5ezfv16b6VKI33KNdc4IgLO3tNLYOZM/Vy/PURMvY5du3YRcX9cRCIRdu3aldOzLBaLQ8GEREROF5ET02y/BlBKfQ04G2h1x+f8vAIc7bGnusey59lnw1XeZePGjQwMDAAwMDCgRfmY8ncVlbo6pvb1EQWm9vU5s9wPVko/d7jq6E926bVD5sYKYsuWLYF2c3Mz48aNIxqNMm7cuEE5tCwWSziGxUeilPoU0A58VkTeDzjtYWCaUupYpVQM+F/A0GYE+pdWHSKrV68OtLdt20ZtbS1KKWpra4dlWAsY3Mvw2tu2QW2t00jX1g7fzPd9+4LtZG6sK6909kPwT6SbpX7qqfrop9eOx+OsW7eOK664gnXr1llfiMWSL8J66XPZcMJ5XwaecLefusePBO71nPdp4M840VuXZXv/OemWUq2tFVFqSMurZlri1kRjY6OUlZVpa4cUhByWj21tbZWampqsliTOiqDVJE0hvqEflT46q6OjQ/ubdXR05Pwsi2UswEgL/y30NgfktQkTDn5COTS2nZ2d0tLSMmQR8TZmRRGTAJEMqn9ra6tWx5zFJJEQKS93Pt/ycl1MQs4HMbFo0SKt/osWLXIfFV78LZaxjBUS/5vyN45K6Q2ZUjl83JmJRCJaYzZcExmLOllx0SL9M3YbdhFxeimRiHM8EnHsHAkS646OjtTnH4lEbI/EYsmSXIRkVKdIue+++5wXhXI6d3XBwoXO3sPEiRONdrEw+Xic/xsC7SFjcqg3NzvzQ6JRZ58HJ3dQdJZ1qFssxWdUC8mZZ57pvMjkdK6rc+acDGViYVcXLF0Ka9c6e4+YfP/739dO9dvFwrRYU+jJikETCzM51NetgyuucPZ5cHI/64vIS9rWoW6xDANhuzKlvOEOd2SV6iKs/6SxUb/O5wcJ41vJN2VlZdrwT1lZ2cHC2lo54KaSOZDt+zb5QQrgUBdx0prU19cPSmfS1NSkvbemPD3PYhmrYH0k6YUkq8YlrP/E5BPIQG1trSiltKSBhcDo9A/zvk3vuQAOdVNurEQioflCbH4siyU3chGSUTm0NQU4jSzThIf1n5x5Jt5WjuQwWgaSqVVEJJVapVCsX7+eK2tr+W+luLK2Vp/1nuF9p13HxJd+RLN37TqYkiYScewcueuuu4x2WVkZSinKykoh04/FMnYZlUJyFM4a6x9JLtNqIqT/ZPvjjzOAk3pkwLWzoaipVbq6uHT7dhaKcOn27XpQgOF9B65jMm2afn+vXQCH+nnnnRdod3d3MzAwgIgwMDBAd3d3zs+zWCwhCduVKeUtmUb+19FoTl09k//kF8uWyXsgvSDvgfxi2bIsb1mrDdcUdHirqkqvvyfFvEiwHycWi2nrsKRCg2fO1O83c6b+vKAJiTkwY8YMiUQig1L4JxIJqaiokGg0KhUVFXZoy2LJEayPJL2QPAOpDymU89vgR0gkEvLxaFQuAfl4NDqkhqxYPpLU3I3k5pnPYppj8vFoVPaDDIDsd9+fiIjMmKHfr8Drs2SaNGnXD7FY8ocVkgAh2VJdLSI5zHY29EgKMYM675FehsiylpYWrf4tnkXANkyd6kRyuZ/jhqlTnYL2dv1+vkiqsAQJQk1NjVbHmpqavDzPYrEMpqBCAswF/hZnedx/Br4AHBb2gcXYkkKS/MVsajQzEpB+pKWlRRv+GdI901Cw1B6NjSJlZYPCk41DbI2NmpCkru3oONhLUyovM9RNKxrmPY2LxWIJpCBCAnwdeAxYDVwKXAD8H+CHwKPAL4Bjwj64kFsqaaMb/ptLIx3US7i7vV0G3IZ2AOTuHH+dm5bMLQTGOSaVlXrPo7LSOZ5IiFRUiESjzj4PQ0pBObOStLS0SEVFRc5CbbFYzBRKSL4NVBjKG4AFYR9cyC0lJB4HeJhhI6MATZ6s/2qfPFm/eIiO52Ln52psbJQOkBdAOvDNMfGKSHJz+V5Li/xzRYV8L08Nu2mui3WoWyzFo9BDW0enOTYl7AOLsc0BkVgs51/MxiEx04Q+0wzwAIqeMbi9XZ/Z7u1RxWL6e3Ojtgox1GQS646ODolGowJINBq1CRgtlgKSi5BkM49kq1LqdqXUBM+xe7O4bvg46qihLZIUkD9q8uTJwfaHP6zfw2ublpYNYP369TQ2NlJWVkZjY6M+ebAQ3HWXvgyvd7JfQ0NqrXdxbYBf/epX2i38dhhmzZqlLWc8a9asVFlzczOxWIxoNEosFrMJGC2WEiUbIdkI/AH4o1LqePeYKlyV8sCUKUMTkaYmuPRSZ+8Rk507d2qnavazz8KMGc5kxRkz9GV9TTPADaxfv56+vr60ItLV1cXChQvp8mUaDs1xxwXaj7uZdMVvi2iX+O0w3HLLLdpyxrfcckuqzCZgtFhGCJm6LMBj7v5jwLPAZ5LHSnWLRqPZD7sYkg2GnsdQU6PfM8ew1ZwiutrbRerrB4fqTpyo13HixFTRjfX10uce7wO5sb5eRHIb2rLJFy2W0oYC+0ge97w+Aqd38n7YBxZjG1JDd9hheoN62GGpIpOPxOgINjirwzBlyhStHlOmTMnuQtO8D0Mdz50yRZu1f67nea9NmCD9+FagzFiN4OSLRfcNWSyWtOQiJNkMbX3a03t5DZgPfGqIHZ9hIbWwlYlx4wLtBtc3kM7u7u6mt7eXgYEBent79VxPvkWXBtlDZPfu3UY7EF+SQ82eMUMv89j/+eabLAD+EVjg2gDU1THl/feJAlPef39Q/rG0iR4xJ1/8xje+oZX5bYvFUvoEColS6u+UUn8HfCn52rX/L84kxZLnTG9G3iVLYNIkZ++lsjLQ/vOf/6wVeW2jI/heXyyC3w5YWTGImpoaox14P1/SQ83+zGf0Mo99yimn8BBwFfCQawPw0ku6E96TcDIw0SPm5IuzZs1KZe8tKyvTnO0Wi2WEENRVAb7n2V712d8L2wUqxjbIR9Laqg/heMsmT9bLPPNBMg27GHM9Bc0jCbFuh9FH0tmph/F6yxIJZ/IgOHtvXerr9Xq4fhDvey8rK9Pe81s1Ndqz3vL4fjKtAd/a2io1NTWDhhttiK/FUhpQ6FxbePwkI2GbM2eO/gmZnN8NDXpZQ0OqqCBpSwLmaLgPFGlpSSsu3a2t8nBNjXT7GuK9M2dqEyP3ejPydnToQuJtpE3iGlCXlpYW2QTSD7LJ5zNKpjlJbuXl5akykz/JTjq0WEqDXIQk2/VIco/z9KCUulYp9bxS6iml1N1KqUFOBKXU0Uqp3yulnlVKPaOUuij0A/2LTnnt007Tyzx2W1sbP2pt5dqaGn7U2kpbW1voKqTo7U1vG9aAp6uLeatWMXf3buatWqWVbX7vPe12mt3cDLGYs0ZILKavEXLCCXo9vHZAXfbs2cM0oAyY5tpJTGvAm/xJNsTXYhkFZKM25DncF2gBytzXVwNXpznnCOAU93Ul8GdgZjb3n3PUUYOHlFpbnZ6I/5d3IuH0CpQaNBv+qc5O+cD9Bf4ByFP56JEccojeEzjkEOe4aQ14Q9mlhx+uDTddevjhg99fuiE2U96sgHVMJk6cqPU6JnpChk1rwCcSCe062+uwWEoPcuiRBK5RqpTayMGeSL1S6qlkkaM/clIO4rXWYz4EfC7NOa8Br7mv31FKPYez+OGz/nMH8corMG8ePPDAwYmJt96a/tx4HC6+2IloOu88bSJj3403Mg7nDUdcm1x7JWvXwkc/qtsA48fr53ltQ1lkzx4GcHoJ/a6tEY+nn5wZj8O6dU4GgOZm/Zx339XPde3Kykre8/R4Kj2BCQe2b0/NkhfXTnZ30znbX3vttcF1slgsIxLT0NbZOJMPPwPM8LxOHs8XfwMY43SVUnXAbCD7vCFZpiahq8s5b/NmZ+8ZNjriyCO1U/12KOJxSCSgo8PZZzOUM3NmoP14VRW9QB/Q69oay5c7S+IuX56+LitWDK7D3Llp7csvv1w77LVfcSOvxGcDvP7669p1fttisYxwwnZlMm3A/cDTabZzPOdcBtwNKMN9DsFJW39ehue1AY8Aj6Sy/9bVZe7PmZaPTSSkv7xcBkD6s0y+GJoJE/R6eCf8GYbfOjvT7qFZAAAgAElEQVQ7tXVRtIAAU2LGTASsYxLk9K+srNQc8ZXJ1PNidsRbLJbSgAKlkf8G8B2P/QqwF3gHWBb2gZ77fQ3oASYYzikH1gB/N5R7p4TEM07/fGOj7Ckrk+f9M6czLR9bgHXI02JYFjdTPYJS5L/1V3+lRXS99Vd/lVsdDaHLpoWyOjs7NZHx1zNdqLHFYikuhRKSh4FJHvtxdz8eeCDsA917fArH1zHZcI4CbgH+baj3TwmJu0jS8+6Kf8lNE5MCLR87ZEzO9pDc4XOA35HrGvFTpuh19KROGT9+vCYk48ePT5W9P2WKVo/3PdfZFCkWS2mQi5CYfCRKRHZ57F+5Q2H7gArDddnwI5xIrN8ppZ5QSv0UQCl1pFIqOQ38Y8BXgE+65zyhlPp0wP0GE4lAezsAUx57DDiYsjhpA7B3r36d3y4W69dDYyOUlTl7fwbgIc6GB1g7YQJ9OC10n2vnhD81i8dOpoJPZ49zfSLKZwM88sgj2nV+22KxlD4mIdHmdohIB4BSKgIcnstDRaReRI4WkQZ3W+Yef1VEPu2+/qOIKBE5yXNe9uugHDgAGzcC8Lqb4iPpCE7ajuFz/A6nI3j9eidIIJ2IBM0xMTB182YiOH9k5do54U/N4rG//e1va0Ve+9XycuDg55+0ASZOnKhd57eD8ndZLJYSIqirAvwY+Jc0x/8F+GnYLlAxttTQlmfm9YO1tfKmUvKgf3hn0SJ9uMa3ZnhJYBhSMnFWTY2WxfesHNPZG9OxiMivWlrkgYoK+ZVvGd7a2lrNR5Kt/6QQKzJaLJb0UKChre8AxyulNiulVrvbZqAe+IdCilveWLwYcBaF2rp9OwdE2Lp9u7441JQp+jV+2xQ+a6Cnp4eVK1fS41t1MRSGISUTn125Usvi+9mVK7XyoS6WtfzFF2nDiX5oc23Pzfjc2rU0ffABn1u7Vus1vfLKK9qM+FdeeSVV9oUVKzgep9d0vGsnKcSKjBaLpQBkUhrgOA7OITk+rGIVc5tTVaX9Wr6nslL7JX2PJzTVFFob1hGf9/xRIXskIsERXWHyiNXX12vX1HsTPVZVaRFiydnwImaH+kDyfHc/AKmyTGHDQYkgLRbL0KFAUVt1xgvdYfewDy7k5k/auDsS0Rqr3dmG1k6dqjfgU6dm+FM4dHR0iFJKAIlEIrlntM2UMTgo/YsB06JdQZgWqBpQShcEpbRrr6ytlf9WSq70DS2+XF6uifzLHrFobW3Vhr28gmGHvSyW/FIoIfkVsBo4HzgB+CvgGOCTwBVAAjgj7IMLufmFxBj+a8KwDK8JU4MbmqDMwJmy+AYQNuw2qBfwXHW19hk/V12t1z1ACCdOnKiJhTd/Vy5p6y0Wy9AoiJA492UmcCXQDbwAPA7cBiwBxod9aKG3QWnkxTAh0USItUNEwv3aD01ZmV7HsrKsLqupqdHqWJOFI940ZNfY2Ch/AtkH8ie/MBmG5kzLCIcd9qqtrRWllObUt1gsZnIREmMaeRF5VkQuE5FmEZkuIrNF5Msicqs480lGDNPXr+fQvj6m+0NrTTz+uNkOYLHr5A+y80p/v9kOoL6+3minw5QO/hvf+AYfw5mt+jF8S+YaggUuv/xyLsBJtnYBev6ulyLOv6f4bAhOW19XV8f27dsREbZv306dbzlgi8WSf7Jdj8QyBNra2mhvb6e+vp729vb8rGMShClrsIFq3zryfjtdRJdpeWHjkrmG+SdtQBew0N17P6lpkQgvAgeAF107yY/uu0+L9vrRfU7ez5c8y/+ms+28FIulAITtypTylm5oK1TOLFNEl/GyIq76lymyLMC3YoraMpUFLS9sXDLXNEQYsPaJiNmP0+8b9up3h71Mc1asg95iCYZCL7U70rZBQpJIHPQllJUNXUyGKEBFX4e8vd1Zcz2diBh8PEGhwXV1dVqDW5dFFuWM4hkULJAhWeWTEybIPpAnvdmQRWSzxwl/AGSzKyS5OOitb8UylimokAB3AWcBkbAPKfY2SEhCRl+FpWTWIQ85/yRsRFRQb8VIyJUhp0yZovU8kk76oJ6KiOOg917jd9B737MVE8tYo9BCcjqwCmeY+ipgetiHFWvzT0iUo47SG6Sjjsrl886KUI1qvonF9PedpSCUlZVpjWpZlpFgoQlY+8TUWwlKrbLJ11PZ5BES05wVpZR2P+WbB2N7K5bRTlGGtoBDgWXAyzhzSL4OlId9cCG3VK6tpJiUSqr4YhOyR2IKDQ4aDisIlZV6/b0ZCQLWiG9sbJS/4IQK/yVLv4qIyBafAG3xlGXqrViRsYwGCi4kwCTgIpwVCH8DfBG4HugO++BCbumSNgb6EUYzmebBBPgtFi1apDWci9xElmHSqhSs/krpZckehCGxpKm3YhIZ3HOTvRWGIDJ20S7LSKHQQ1t34yxCtQI4wlcW+sGF3Ab1SMYyQU5uQyOdSCRSw1tlZWWp4bmiTrTMVP/aWr3+yQY8QxTYfre3st/XWzGJjKnMNCRmijqzvRhLqVFoIZkf9ubDtQ3ykYwG8t2jclctTP4C9w97pfPxDIuQmKitdXoi3sbYFAUWMBwmYhYZU28lKHrMqUpEE5mIWxc7VGYpRQotJOel2RYAfxX2oYXe0s4jKQBF8xcUwMfTF41qjWNfNJrxmkMOOURrAA855JCc65F3TFFgQcNhIkaR2eImpEz5TzzXmUQmqCeTaagsaB5MsjxIZGw2ZEsuFFpI/gvYjZPAcTWwC1iLk5niK2EfXMitGEKSi79gyBFdIbMQm/iLT0j+4heSNENKmaK5iuqINxEUBRY0HCZiFBnT3JQ3fGLxhkcUgkQm7DBaLpMtrQBZMlFoIVkLfMhjfwhnbaMa4OmwDy7kVgwhCTvMk0gkUgkHy8vLsxOTAsyDuaGhQWuwbmhoOFgY4D/JFM1VVEd8WNINhyWPG0RGGwb0iMwH48drZR+MH58q2+8Thf2uKJh6MWF6OCIi0WhUE5mo54dBJgEKStWf/HunS+hph99GH4UWkmd9tkoeAx4P++BCbqXcIwmKiDKSSIhEo86fKxod2sz8ADo7O+UCkPtALvDXv65Ob1Tdme2mupec/yQMYUTGt9aN5pMJEKBenyD0egQhSHxEzNmQw/Zygspqamo0gfGKSabhN2+ZH1NZGOGyvan8UWgh+THwW+Cr7vYb99hE4PdhH1zIrZR9JGHXAQmVK8xAR0eHRCIRgTSLb1VU6A1nRYVbhfTRXCIjqEcSliCRMflkggTI0MMxiUWfr6wvS5EJU5ZvYcpUVlNTowU7+IXLW5YUk9bWVu24X0y8ZX72uWX7hih2YcpM10yYMEEAmeBLAZT8TML0CMOKa6GFRAGLgevc7XOACvvAYmzFEpIwFGTRqxAY07gceqjeAB56qHZdOv/OmE6IGOSTEUkvQKYejsHpb+r9mMQizFBavoUpU5mpJxZUFuYaEUdEvGVeMcm3SJqumTBhgiZ2XjEJI6wi4cV1P8gc51j+hQSIFqLXAVwLPA88hTNPpTpDHR4Hfpvt/UtZSLzL8CqlCp7Qsb29Xerr69MKVqDT3xQlFjCvY+LEiZqQeFc6dC4rEUd8qRDUwzGVZej9BAqQoSwokCDfwpSpLIwAFVvs8l3HsEKY77Lk8YIJiTgN+Trg0LAPCLhnC1Dmvr4auNpw7t/hrMqYtZAcddRRRclxFSafVjETOubU+0k3b8UwidG00uGoH/YqJkPt/WRR9lZNjQygR6OFFaawZWH8RmEb1OSzJM2z8i2S+RbPQpQljxdaSH4NvATcCPwwuYV9YJr7nwusCiib6grZJ4ciJEOKiApJqOgrz7XFSOhYX1+vNeD19fW53dCQu8skFqPCET8WCSlMYct63cbN27AnMQ3JBPlBAssSCelzy/pA9zvmWyRDimfYgIxS7pF8Nd0W9oFp7v+fwJKAsjuBOUDzUITkBJAOsoyIykSAkztU9FWRybvfIkM24aDhq9ABBhZLoTAFr+RbJEOKZxhhDVuWq4/EWRvVgIj8QilVARwjIi9kOj+JUup+YEqaostE5NfuOZfhLLu9Ks31ZwNviMijSqnmLJ7XhrtS6xzgEuBnDz6YbXXT09MD8+dDby/EYvD730M8DsCrr76qneq3w7B8+XLuuusuzjvvPK6++uqsr+vq6mL16tUsXrxYW9b3hBNO0M7z20OmpgZef123PaT+AD42b95stIPqb7EUjHg89V0exLZtwdeFKTNcU+78YE675rmpLOaWxfJUFhPhUaUeDaxoJjIpDfAZ4AVgq2s3AL8Jq1ye+34N6AEmBJSvBHYA24DXgfeBW7O59xy3q7bXs95EKJYt03+BL1uWKsr3uH9Yf0amZXHD+mPSOulN2XgNZaYeifWfWCylAQUO/30UZy2Sxz3HcprRDnwKJ6Pw5CzPb2YoUVvJccVcc0EtWqQ3jr7hq3ymCJ86darWoE71pUEJir7KNGwUxh9jFLWgbLwG/4nJR2L9JxZLaVBoIXnI3XuF5KmwD3Sv34yzQNYT7vZT9/iRwL1pzg8lJDknNzT0SPLtf2hqatLu1+RJg2Jq2E3XhSWUk97gPzH1OmyPxGIpDQotJDcCX8aZ8zENZ0Grn4Z9YDG2OePGDUlE7m5vlxvr6+Vu/zWJhNMgKuXsPb/qTXmnwpBIJOQqkBdArgKtB2Fq2JctW6aVLfOIXVhCiWSG1Ri7W1vl4Zoa6fbdK9Oz7PwTi6U4FFpIJgBXAg/jrJB4JTA+7AOLsaWbkBg0NHR3e7t8gBPf/QEMFpOAoZy8RyK1t+thgp56mBrbRCIhsVhMlFISi8XyElLc0dGhPS+rSZMh/SemiYy2t2KxFI+CCslI3PxCYhoaur26Wpukc3t19cELEwknz1Q06uwLuciTIVV8poY93/NSQjvpQ/hPTBMZrf/EYikeuQhJusgyDaXUXyulupRSa5VS/5PcMl03rGzaBF1dKfO2227Tir32/v37tTLN7u52Qn8HBpx9d3eqqKGhQbvObw+Z444LtJubm6moqCAajVJRUUFzc7N2ajweZ8WKFcSDwhmHSDweZ926dVxxxRWsW7cu+/u2tcGaNc7ey+7dgfbll1+uFXntPXv2aGV+u6uri4ULF9Ll+VtbLJbik3EeCfAr4KfAz4CBwlYnT+zdC0uXOq/b2jjuuOPYsWNHqvg4TyP91oEDgPOT12sD0NzszB9JziPxNODV1dUopRARlFJUV1fnVuerroJPfMIRrWjUsV2SDXt3dzfNzc15EwwT8Xg8f88xzD9pa2tj+oMPMvG++3jvzDOZ5xGhZ555RruN1+7q6mKp+zdeu3Zt6l4Wi2UYyNRlAR4N290Zrm1Ocghl5kwRcYZqotGogLPgj3eo5hlfSu9nvMuvOhennQFbkJxZeU4VXzKE9J/YYS+LpXhQYGf7PwH/GzgCZ1XEGqAm7AOLsaWEZMaM1IcU5Ed4JuncTgoJZP3BFytn1qggyH9SVaULSVWV55JwYcM20stiGTqFFpKtabYtYR9YjG1Oul++6bLZish148Zp0VLXjRs39L+AJTyRiC4k3hUGJThsOChqzkZ6WSzhKKiQjMRtTlXVYBHxNlbZRm1ZCo9pjQ3DsFdypcbkVlZWJiJ2yMtiCUsuQhIYtaWUave8/ryvrGNonpgiM22aHjl01116ucc+6aSTtCK/bSkw69dDYyOUlTn79esPln3nO/q5HvuUU07RipK2jfSyWIaBIIUBHkv3Op1datugCYmtrfovW+8wSSIh+5SSfpB9So0+R/dIJsOw13/V1MgbIP/lySpgJzhaLOGgQPNIVMDrdHZp40+f7rGX3HADzSL8f0CzCEtuuKG4dbMEM3dusL1kCZ/evZvJwKd374YlSwCorKzULvHaq1ev1sr8tu2tWCwhCVIYRlOPxDBDPdNa40aCIpFGCKb13EuGoKVlA5JEdnZ2ygUg94Fc4Ot12HT2FkswFMLZjjP5cC/wDs7iU3s9dl/YBxZjS5drK2iORmVlpdaAVFZWZvepm+ZGjAByWs+9FCgv1z//5NoznZ16zrIs83plctLbkGLLaKcgQjKSt7RCEsCMGTO0BmSGZ+6JkZYWvSEbYdFBeV/PvdgE+b0M81KmTJmi9Va8Exxtb8Uy1slFSDLm2hrtXHzxxUY7kMWLzXaJc+qppxrtkufWW6G11Um30trq2ADvvquf57HvWLCALmAh0OXaSUzpWKxvxWIxk02urVFNMj/TkNcMb2uDF190QonPO29wosISJ+/ruQ8HSfHwMncubNig2y7z/vM/EZxIEXHtJJWVlXzpvfdYDKwGfutx0ptCim3OL4sFO7QVGoMDPyMBs+yLSUFyhZUKQQ56Qzhxd2ur5lvxzqS3vhXLWADrIxkGIenocFZOTDZI2Sz+JGKcZV9sxlyuMNMs+qoqLcOB37fyJ5B9IH/y+VZqa2s1IamtrU2V2XxglpGEFZLhEJKwgmBYwMpSBEL0Vl6trdV6K696xKKsrEx+AfIGyC84mKpFJDh7cSbnvRUZy3BghWQ4hMT069ZEU5N+XVNTYetpyQ7T3zMS0XsrHpH5r5oaTWS8s+xjsZgWJRZz57qY0uNbkbEMF1ZIchQS05czcPinoUFveBoasntYIuH4VcDZj5VhpZFAUG+lslL/W3vnGsViusi4YiEi8veVlZrI/L17XZDAiFiRsQwfVkhyEJLOzk45DeQSkNN8X85EIiHl5eUCSHl5uS4mYXskzo1H5wJWoxXT5NOgiZEi8sH48ZrIfDB+vIgEC4yIIzJef4wVGUuxsEKSg5B8qa5OPgDpB/kA5Et1damyRYsWaV/ORYsWHbywhJzmliIQlA7HlBA0YEgsSGBERDZEo5rIbIhGU2WxWEw2uf+rm3wiU1VVpf2vVmW5QJiIMxmzrKxMm4SZvM6Kz9hhxAkJcC3wPPAUcDdQHXBeNXCne+5zQDyb+w9FSEzrkZhmO8uyZXoDsmyZdt8xFxE1lmltFamp0UVEJLjXavC59PuWfu73LP28xS1Lbls8ZZFIRBOZiOeeVVVVWplXZMIuEBYkPslrrQCNPEaikLQAZe7rq4GrA877BXCB+zoWJDj+bShC8mJDg/bFfdHj6zDmozIIyaieo2EZGun8LqZhUXfYKyUynmGvAfRloQcgVfZyebkmMi97htg2eY4fcMUkSZAAZRIfr49nKOlkamtrRSmlhUl7r7UCNHyMOCHRKgDnAqvSHD8UZ1lfNdR7DsnZvmyZ/sX1CEJHR4copQQQpZR0eOeKJBIHx8fLyzV/h/e6SCSiX2exiAQ79k3+GF8YsngbY19PRjy9FZMAvegTmRfdMpP4fNNX9k1PmUmAamtrtTL/nBtvmV9MltfUyH+DLPdExSXvGSRMYUVrrAraSBeS/wSWpDneAGwAbgYeB34GTDTcpw14BHjkmGOOyf7TSyScSBulnL1HEIw9C4OQtLe3aw78EZdZ1zK8mJYnqK11/lf9jWNtrS5A3nKDAAWJjEl83vKVveUpMwmQqSxI0EQcEfGWJcXEJExhRSuToG2NRKQfZKtvkTVTmemaIIFMvocwIhlGQN2e5NtSakIC3A88nWY7x3POZa6PZFCvA5iLk77+VNf+AXBFNs8e8jwSwxc30NdhGNq6qLFRc+BfNJSIrjwzVn9djUmCRMZUFiQyBvHp9wlJv6fRNwmQ6TpT2W5f2e4sek1hRctUttX1bSU3rzAElZmuCRJI5+MPJ5JhBDR5fI7zXktLSDI+GL4G9AATAsqnANs89ieA/8rm3umEJLBBDZszyyAkv6uv1/7xfzdMKdpt+nNLVph6OemONzbqIuP9oWQafjOUmXw8fT4h6XMb97DCVMyyMAIpEl4kwwho8viIExLgU8CzwOQM5/0BmO6+/ifg2mzu7xcSY4Pa0aFPEMzWn2EYEvvD5MnaP8gfJk/O7p55xhh1ZrHkQpCPRyRcz0gcMen3iUjyWWmFqwCiZSrLd48kSCBFhkfsRqKQbAZeBp5wt5+6x48E7vWc1+D6PZ4C7gEOy+b+fiGZOXOm1qDOnDnzYGEuWXwDJhY+09Sk/cGeGaY0KE1NTdr7brLpWCwjlSDhyrdoZSjLq4+kAD27MAKaPD7ihKTQm19IMq6CmO+Z5omE9EWjMgDSN4xpUJYtW6a972W+uS4Wi2WYKUDPLoyAvlxeLqfkICRjYoXEjKsgxuOwYoWzzxNl0SgRpSiLRvN2z6Fy/vnnE4vFUEoRi8U4//zzh60uFoslDevXQ1+fs/ezbRscOODs81Q2tbeXqAhTe3sHHX8MHh36G3AYEyskhl4FMQM9PT10d3fT3NxM3CtC3d0wMOAMbg0MOHYeRSpb4vE43d3d6etosVgseUKJ44sYVcydO1ceeeSRrM8PFIQM1yxYsIDe3l5isRjr1q07eG1PD71NTUT6+zlQVkbswQeHRUgsFoslW5RSj4rI3MxnDmZMDG2Bs7b2woUL6erq0o739PQwf/58LrvsMubPn09PT09W9+vu7qa3t5eBgQF6e3vp7u5OlS2/5x7m9ffzXWBefz/L77knj+/EYrFYSosxISRdXV0sXbqUtWvXsnTpUk1MbrnlFvbv34+IsH//fm655Zas7tnc3EwsFiMajRKLxWhubk6V3XzzzTwEXAU85NoWi8UyWhkTQvKDH/zAaIchHo9z4YUXcuyxx3LhhRdqQ2IVFRXauX7bYrFYRhNjQkj8fiCvXVVVpZX57SC6urq45ppr2Lx5M9dcc43Wy7n00ku1c/22xWKxjCbGhJBcfPHFXADcB1yAHv7729/+VjvXb9PTAytXOnsPq1evDrTb2tpobW2lpqaG1tbWvEWJWSwWSykyNsJ/gW+6rxcCylNm6q3Q0wMLFkBvL8RisG5dKvpq8eLFrF27NnXq4sWLU6+7urpYtWoVAKtWraKpqcmKicViGbWMiR7JU9/7HnBQQJI2ZJis2N3tiMjAgLP3RGa1tbXR3t5OfX097e3tmlCYeisWi8Uy2hgTQvLzvXsBJ0+I1waYNWsWUXf2eTQaZdasWQcvbG6mPxKhH+iPRMATmdXT08P111/P1q1buf7667WwYW/vJJ1tsVgso4kxISRvnnsubcAanGGuN889N1V2yy23MDAwAMDAwIAW/tu1cSPf6uvjfuBbfX10bdyYKjPNI2lra6Ozs5OWlhY6OzvtsJbFYhnVjAkhufXWW/mgtZXWmho+aG3l1ltvzeq6Z2+8keuBM4DrXTuJaR4JOGKyZs0aKyIWi2XUMyaEBBwx2bVr1yARMSU2/Ob48YwDosA4104Sj8dZt24dV1xxhZ4eZRhYvnw506ZNY/ny5cNWB4vFMnYZE1FbAA8sWcLE++7jvTPPZJ5HTEyJDU+YORN58EHN9hKPx4c9EeLy5cu55pprAFL7q6++ejirZLFYxhhjImnjA0uW0OSG4wI82NqqiUkgPT2Og72vD8rL85bFN0ySyCCOPvpoduzYkbKnTp3Kyy+/nGsVLRbLGCOXpI1jokcSdZMmKpzIrWi2SRTjce65+GJ233UXNeedx6I8iUhg1uAQHHfccZqQHHfccTnX0WKxWIbCmPCRrK2sBA6G/ybtTHR1dXHuNdfwjc2bOdeXBiUTPT09rFy5clA2YVO0VxiuuuoqLXz5qquuyul+FovFMlTGRI9k6uWX07Z0KYuB1cBHLr88q+tWr17NaUAz0O3a2URhmXodzc3NRCIRBgYGiEQig6K9hko8HucPf/iDXbzKUjL09fWxY8cO9u3bN9xVsaRh/PjxTJ06lfLy8rzdc0wISbLxv26IKyR+q6GBlrVriQG9wNqGhqyuS9frSDbwGzdupK+vD3C+cBs3bsy58S8Fp7/FkmTHjh1UVlZSV1eHUirzBZaiISLs2rWLHTt2cOyxx+btvmNiaAvCzetYVF3NeKUoA8ZHIiyqrs7quubmZtqU4r+BNqW0XodNn2IZ7ezbt49JkyZZESlBlFJMmjQp773FMSMkoWhuJjJ+PESjRMaN01KkmIhv3MgN/f20ADf09xP3zIi36VMsYwErIqVLIf42wzK0pZS6FvgMzojRi8DXRWRPmvP+FifzuwAb3fOKN/AajzsZf7u7HRHJdvho9WotwzCrV4PbE0r2iFYPcZjNYrFYSpXh6pH8DjhRRE4C/gys8J+glDoK+L/AXBE5EWeC+f8qai3BEY8VK4Y2f8Tfy/DZNn2KxVJYfvjDHzJjxgxaW1uHuypjgmHpkYjIWo/5EPC5gFPLgAqlVB8wAXi10HXLC7NmQVkZ9Pc7e29GYYvFMoh8TtIF+PGPf8z999/P1KlTteP9/f2UlY2JGKOiUgo+kr/BWbxQQ0ReAb4PvAS8BrztEyANpVSbUuoRpdQjO3fuLFhls6K7G+nvB3D2Oc4VsVhGM8lw+e9+97ssWLBg0NyrobJs2TK2bNnCmWeeyXXXXcc//dM/8ZWvfIWPfexjfOUrX2FgYIDvfOc7fOQjH+Gkk06is7Mzde21116bOv49z7pFSQYGBvja177GiSeeyKxZs7juuusAJ8DmoosuoqGhgRNPPJENGzYAsGHDBuLxOLNnz+ajH/0oL7zwQuo+//AP/8CJJ57ISSedxPXXXw/Ao48+yrx585gzZw4LFy7ktddey+mzKBoiUpANuB94Os12juecy4C7cVO1+K4/DPgfYDJQDtwDLMnm2XPmzJHh5PnGRjkAqe35xsZhrY/FUkyeffbZIZ3f0dEh0WhUAIlGo9LR0ZFzHWpra2Xnzp0iIvK9731PTjnlFHn//fdFRKSzs1OuuOIKERHZt2+fzJkzR7Zs2SJr1qyRb37zm3LgwAEZGBiQs846Sx544AHtvo888oicfvrpKfutt94SEZF58+bJBRdcICIiDzzwgJxwwgkiIvL2229LX1+fiIj87ne/k/POO09ERH784x/L4sWLU2W7du2S3t5eicfj8sYbb4iIyB133AjeEfoAABUvSURBVCFf//rXc/4s0pHubwQ8IiHb+4L18UTkdFO5UuprwNnAAvdN+Dkd2CoiO93z7wI+CmSXA95PV5fj9F68OOX4zpmenrSO+MOefBI4mJIlaR+8LL/deItlJJNckiE5gTfXSbrp+OxnP0tFRQUAa9eu5amnnuLOO+8E4O2332bTpk2sXbuWtWvXMnv2bADeffddNm3aRFNTU+o+xx13HFu2bOHCCy/krLPOoqWlJVX2pS99CYCmpib27t3Lnj17eOedd/jqV7/Kpk2bUEql5pDdf//9LFu2LDXMVlNTw9NPP83TTz/NGWecATi9liOOOCLvn0UhGK6orU8B7cA8EXk/4LSXgNOUUhOAD4AFwCMB55rp6oKlS53XyXXWcxUTw3ruavJk8OS/UpMney4Ll2vLio9ltJJckqGQ/98TJ05MvRYRrr/+ehYuXKids2bNGlasWMHSZFuRhsMOO4wnn3ySNWvW8NOf/pT/+I//4Oc//zkwOKxWKcV3v/td5s+fz9133822bduMIikinHDCCTkP7Q0Hw+Uj+RFQCfxOKfWEUuqnAEqpI5VS9wKIyHrgTuAxnNDfCJB9sisvP/iB2Q6DYT33yWefrZ3qtbu7u9m3bx8DAwPs378/q1xbPT09zJs3j0svvZR58+aNyH80i8VEPB5nxYoVRfmRtHDhQn7yk5+kegd//vOfee+991i4cCE///nPeffddwF45ZVXeOONN7Rr33zzTQ4cOMDixYv5l3/5Fx577LFU2S9/+UsA/vjHP3LooYdy6KGH8vbbb3PUUUcBcPPNN6fOPeOMM+js7KTf9aXu3r2b6dOns3PnztT3u6+vj2eeeaYwH0KeGRYhEZF6ETlaRBrcbZl7/FUR+bTnvO+JyIdF5EQR+YqI7A/5QLNtYMmSJUyaNIklS5boBc3NTk8kGnX23l8a55/PgfJyDgAHysvBs1jWnj17kj4gDhw4wJ49B6fPBCV6vOaaa7S0Ksl1RywWy9C54IILmDlzJqeccgonnngiS5cupb+/n5aWFr785S8Tj8eZNWsWn/vc53jnnXe0a1955RWam5tpaGhgyZIlrFy5MlU2fvx4Zs+ezbJly7jRXU21vb2dFStWMHv27JRoJOtwzDHHcNJJJ3HyySdz2223EYvFuPPOO1m+fDknn3wyDQ0NJBKJ4nwouRLWuVLK2yBne2eniCMfztbZOcjRlI7W1lbBcXMIIK2trfoJiYRIR4ez1w4nZF4sJpcqJfNiMUl4yltaWrR7trS0pK6pqKiQaDQqFRUV2jUzZszQrpkxY0ZW9bdYhoOhOttHA/PmzZOHH354uKuRNSPG2V5SJP0hQ3S233fffUabeDztRMXu7m7+ODDAAyJEBwa0pI2TPf4Sr50c8hKR1JBX8prp06fz3HPPpa6ZPn16VvW3WCyWYlAK80iKQ1sbrFkzJCf7mWeeabTp6YGVK529h2QUSjQaHRSF4p/jkrRNQ17t7e2p6I6ysjLa29uzfg8Wi6XwdHd3M3duqMUFRwVjo0cSklvd5Xjvu+8+zjzzzJQNOOIxf/7BqK3f/z7VOzFFoSxevJi1a9dqNsATTzyhPdtrx+NxHnzwQRu1ZbFYShIrJBm4NWht91tugf2u73//fsf2NPBDXSMkSGDC3s9isViKhRWSIpOM5vDabW1tNiuwxWIZsVghCcv558PPfw59feAL8TVx5JFHBtpeQbFYLJaRwthxtuebeNyZhHjllc4+y2En6zi3WEqfnTt3cuqppzJ79mz+8Ic/DOnaJ554gnvvvTen5//jP/4j999/f073KCa2R5ILAeG/5kuCHec2DYplzBKQt2446O/vZ926dcyaNYuf/exnQ77+iSee4JFHHuHTn/505pMD+Od//ue0xwcGBohGo1pds0mLX/D0+WEnoJTyNtzZfzPR2dkpLS0t0umZGGmakGixjCSGPCExkRCpqBCJRp19jv/7W7dulenTp8uXv/xl+fCHPyyLFy+W9957T0Sc7L1NTU1yyimnSEtLi7z66qsi4kwovOiii2TOnDny/e9/X44++mg5/PDD5eSTT5b3339f1qxZI6eddprMnj1bPve5z8k777wjIiIbNmyQeDwuJ510knzkIx+RPXv2aNfecccdWt1uuukmOeecc+T000+X2tpauf766+Vf//VfpaGhQU499VTZtWuXiIh89atflV/96lci4mQybm9vl9mzZ8vtt98+qK5bt26V+fPny6xZs+STn/ykbN++PXWPpUuXSmNjo/zt3/6tVo98T0gc9ka/EFvRhCRgZruJzs5ObZZ6UkwKkUrbYhkOhiwkHR2OiICzz/F/f+vWrQLIH//4RxER+frXvy7XXnutMU37vHnz5Fvf+lbqHjfddJN8+9vfFhGRnTt3yic+8Ql59913RUTkqquukssvv1z2798vxx57rGzYsEFEDqaM917r56abbpLjjz9e9u7dK2+88YZUVVXJT37yExERufjii+W6664TkcFCcvXVV6fu4a/r2WefLTfffLOIiNx4441yzjnnpO5x1llnSX9//6B62JntpYIh+6+J1atXD7Lb2tqKkkrbYilJknnrkt+lPPzvH3300XzsYx8DnHx5P/zhD/nUpz5lTNP+xS9+Me29HnroIZ599tnU/Xp7e4nH47zwwgscccQRfOQjHwGgqqoqq7rNnz+fyspKKisrOfTQQ/nMZz4DwKxZs3jqqafSXuOvm9fu6enhrrvuAuArX/mK5nf9/Oc/rw2FFQorJGFJl/03CyEJmi9SjFTaFktJEo87P8Ty6CNJl9JdxJym3Ztq3ouIcMYZZ3D77bdrxzdu3BiqbuPGjUu9jkQiKTsSiWiJHU11C6prpusKhY3aCktzM/2RCP1AfySS9a+otrY2Ojs7aWlpobOzUwv3LWYqbYulpIjHYcWKvDnaX3rppZRg3HbbbXz84x8Pnab9tNNO409/+hObN28G4L333uPPf/4z06dP57XXXuPhhx8G4J133qG/v5/KyspBWYMLyUc/+lHuuOMOAFatWsUnPvGJoj07yagUktdff73ga3Z0bdzIJ/r6+C7wib4+uobw66StrY01a9bYOSMWS4GYPn06N9xwAzNmzOCtt97iW9/6Vug07ZMnT+bmm2/mS1/6EieddBLxeJznn3+eWCzGL3/5Sy688EJOPvlkzjjjDPbt28f8+fN59tlnaWhoSK1RUkiuv/56brrpJk466ST+/d//nR/kY72lIaIcH8voQiklFRUVWa8+GIaFCxdqQ1QtLS2sWbPm4AklFM5osRST5557jhkzZgzb87dt28bZZ5/N008/PWx1KHXS/Y2UUo+KSKjMk6OyRwKOQyyb1QfD4s+FpdlJR/x3v+vs7YqGFotlFDNqhWRYI58My/BaLJbCUldXZ3sjRWZUCklVVRX/9m//VlCndbow3hSmZXiBrq4uFi5cSFdXuCXoLZZSZzQOmY8WCvG3GZXhv3v37uXiiy9m1qxZBRMTY9p3QzhjV1cXS5cuBUhdb53ultHE+PHj2bVrF5MmTRoUhmsZXkSEXbt2MX78+Lzed1QKCRz0kRRKSDKmfQ/IwxU0IdFiGS1MnTqVHTt2DFoN1FIajB8/nqlTp+b1nqNWSCKRSPY+kq6uIa/nDuHSvjc0NGg9mYaGhiFdb7GUOuXl5Rx77LHDXQ1LERk2H4lS6gql1FNKqSeUUmuVUkcGnPdVpdQmd/tqtvfv6+vLbuZpVxcsXQpr1zr7fPktAtZzr66uZiXwArDStQ9WxfpOLBbLCCRskq5cN6DK8/r/Aj9Nc04NsMXdH+a+PiyLewsgLS0tBzOSBSVYbGlxksUlN+81YTFkM325tVUOQGp7ubVVRIKTOVosFksxIIekjcPWIxGRvR5zIk4D6mch8DsR2S0ibwG/Az6V7TNSDvCeHmhqgksvdfbeXoJvPsggOwzd3c467gMDzt4T/jt1/XoAlM9OtwSvxWKxjASGdWa7UupK4HzgbWC+iOz0lf8DMF5E/sW1vwt8ICLfT3OvNiDpsJgNvAy8CfBhmD4RDkme+x68+7wzusSH4PCpUJss2wHb/+JeF5aAewK8eTzUVcOkZNke2PUibAOOB6o9t9kDvJhLPbLgcHJ8rwXA1il7SrFetk7ZUYp1mi4ilWEuLKizXSl1PzAlTdFlIvJrEbkMuEwptQL4P8D3wj5LRLqALve5j0jIqf6FwtYpO2ydsqcU62XrlB2lWqew1xZUSETk9CxPXQXcy2AheQVo9thTge6cK2axWCyWvDGcUVvTPOY5wPNpTlsDtCilDlNKHQa0uMcsFovFUiIM5zySq5RS04EDwHZgGYBSai6wTEQuEJHdSqkrgIfda/5ZRHZnce9SjJ+1dcoOW6fsKcV62Tplx6iq06hMI2+xWCyW4jEqkzZaLBaLpXhYIbFYLBZLToxYIVFK/Vwp9YZSKu3CA0qpZqXU224KlieUUv9YhDodrZT6vVLqWaXUM0qpi9Kco5RSP1RKbXZTxJxSAnUq6mellBqvlNqglHrSrdPlac4Zp5T6pfs5rVdK1ZVAnb6mlNrp+ZwuKGSdPM+NKqUeV0r9Nk1ZUT+nLOs0XJ/TNqXURveZg0JZi/3dy7JOw9FOVSul7lRKPa+Uek4pFfeVD/1zCjslfrg3oAk4BXg6oLwZ+G2R63QEcIr7uhL4MzDTd86ngftwJrefBqwvgToV9bNy3/sh7utyYD1wmu+c/42bNgf4X8AvS6BOXwN+VMz/Kfe5fwfclu5vVOzPKcs6DdfntA043FBe1O9elnUajnbqF8AF7usYUJ3r5zRieyQi8iCQTQRX0RCR10TkMff1O8BzwFG+084BbhGHh4BqpdQRw1ynouK+93dds9zd/FEf5+D8wwPcCSxQqnCLW2RZp6KjlJoKnAX8LOCUon5OWdapVCnqd68UUUodivMj/EYAEekVkT2+04b8OY1YIcmSuDtUcZ9S6oRiPtgdYpiN88vWy1E46VuS7KBIDbuhTlDkz8odGnkCeAMnn1rg5yQi/ThpdCZRQLKoE8Bit7t/p1Lq6ELWx+XfgHacMPl0FP1zyqJOUPzPCRzhX6uUelQ5KZP8DMd3L1OdoLjfvWOBncBN7tDkz5RSE33nDPlzGs1C8hhQKyInA9cD9xTrwUqpQ4DVwMWiJ6ccNjLUqeiflYgMiEgDTraCRqXUiYV+Zh7q9J9AnYichJNA9Bf+e+QTpdTZwBsi8mghnzMUsqxTUT8nDx8XkVOAM4FvK6WaivRcE5nqVOzvXhmOS+AnIjIbeA+4JNebjlohEZG9yaEKEbkXKFdKHV7o5yqlynEa7FUicleaU14BvL/QprrHhq1Ow/VZuc/bA/yewVmdU5+TUqoMOBTYNZx1EpFdIrLfNX8GzClwVT4GfFYptQ24A/ikUupW3znF/pwy1mkYPqfkc19x928AdwONvlOK/t3LVKdh+O7tAHZ4ett34giLlyF/TqNWSJRSU5JjxUqpRpz3WtCGyH3ejcBzIvL/B5z2G+B8NzLiNOBtEXltOOtU7M9KKTVZKVXtvq4AzmBwipzfAMmFzD4H/I+4nsDhqpNvnPizOP6mgiEiK0RkqojU4TjS/0dElvhOK+rnlE2div05uc+cqJSqTL7GSafkj+gs9ncvY52K/d0TkdeBl5WTVQRgAfCs77Qhf04jdqldpdTtOBEPhyulduAkfCwHEPl/7d1biFV1FMfx78/L6DwpakURJYZgiCn60sUguzz00kuFQmraQ2gXoXop8akLFPPSg5QhPth0IQWzyArJ0qiQQktNKhNCqJCyLBzUdKbVw3+daXfScWTPOOL5fWAz++z/vvw5w9lr/v99Zq1YRflQLZXUDRwD5g3mByzdACwA9uRcO8By4IpKv96lfCtiP3AUWHwe9Olcv1eXAmslDad8cNZFxDuSnqQU13mbEvw6Je2nfKli3iD2p799WibpDqA7+7RokPt0SkP8PvWnT0PxPl0CvJn35BHAaxHxvqQlMGSfvf70aSjuUw8Dr0pqoxQLXFz3fXKKFDMzq+WCndoyM7Nzw4HEzMxqcSAxM7NaHEjMzKwWBxIzM6vFgcTMzGpxIDEzs1ocSKwlSOpRqfewNxPkPSZpWLZ91sdxYyU9cO56+r/rt0valv8oWfdcbZI+zlQqZgPGgcRaxbGImBERUynpT26nZEMgIq7v47ixlJofQ+U+YENE9NQ9UUScALYAc2v3yqzCgcRaTibQux94KPMJdUFvbqRNOWL5WtJc4FngqhzNdOR+GzMt+F5lanBJE1Wqza3O7ZszZxfZvlAlrfouSZ25bb5KVcavJL10mlHHPcBblfNslTQl18erUiFU0npJKyV9IumApNmSOiXtk7Qmd9uY5zQbOM2Vrrx4uRAXoOsU2/6g5EPqytd3Aqsr7WOAiTRV4QTG5c92ShK+8blfNzAj29YB83N9KqUy5YTG8cDVlHTrI3PbC8DCpuu0AQebtv0IDMv1OcDrlbZvgUdzfTnwHSWH2AjgIDAKGA78OtS/Dy8X1uIRidm/9gC3SXpO0o0R8edp9lsmaRewnZJue3Ju/yEiGokxd1CCC8DNwPqIOAQQEb9Tsq7OAr7IZJq3AJOarjOBEuwAkHQl8FNENApKXQPszrbRlGm457MtgDVRKmR2Az3AiShTZCcaWWnNBoIfullLkjSJcnP9pbEtIvZJmknJfPq0pC3Ay03H3QTcClwXEUclbQVGZ/NflV17KCOW03YBWBsRT/Sxz7HKuQGmk4EjzQLeyPWpwM5KkJkOvJh9vhz4OSIaGVpHAcf7uK7ZWfGIxFqOpIuAVcDKys0VSZcBRyPiFaCDUvDnCFD9630McDiDyBTg2n5c8kPgbknj8zrjKA+975J0cWNbjjh6RcRhYHiONgBmkIFF0mRKbe092TYN2FU5vHe0QiUAZR8ORcTJfvTbrF88IrFW0Z5TSCMpzzI6geZCX9OADkl/AyeBpRHxm6RP86H2e8AKYImkbyjPILaf6cIRsVfSM8A2ST3AlxGxSNIKSj3vYXm9B4EDTYdvBmYDH1ACwvGcVttNKUh0L/BU9v1z6J3mas9ABP8NKnOATWfqs9nZcD0Ss/NYTrU9EhELJH0PzIyIIzXOtwF4PCL2DVgnreV5asvsPBYRO4GPJI0pL2sFkTZgo4OIDTSPSMzMrBaPSMzMrBYHEjMzq8WBxMzManEgMTOzWhxIzMysFgcSMzOrxYHEzMxq+QdHQh6b0iiiSgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10591fa90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot pair potentials\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.set_xlabel(r\"Distance ($\\mu m$)\")\n",
    "ax.set_ylabel(r\"Energy (GHz)\")\n",
    "ax.set_xlim(np.min(distance_atom), np.max(distance_atom))\n",
    "ax.set_ylim(-3, -1.6)\n",
    "ax.plot(x, y[0], \"ko\", ms=3, label=\"free space\")\n",
    "ax.plot(x, y[1], \"ro\", ms=3, label=\"perfect mirror\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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 (main, Mar 10 2023, 10:55:28) [GCC 11.3.0]"
  },
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}