formelsammlung/scripts/crystal_lattices-Copy1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "eaed683c-c6f1-45e4-aaee-ae4e57209f5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy as scp\n",
    "from scipy.spatial import Voronoi, voronoi_plot_2d\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "efd84ecd-9fb3-4d2f-9f7a-47cd1ff09eea",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Lattice:\n",
    "    def __init__(self, *vecs):\n",
    "        # if the vecs were put in an iterable\n",
    "        if len(vecs) == 1:\n",
    "            vecs = vecs[0]\n",
    "        if len(vecs)  == 3:\n",
    "            pass\n",
    "        elif len(vecs) == 2:\n",
    "            pass\n",
    "        else: raise ValueError(\"Vecs must contain either 2 or 3 vectors\")\n",
    "        self.dim = len(vecs)\n",
    "        self.vecs = list(vecs)\n",
    "        for i, v in enumerate(self.vecs):\n",
    "            if type(v) != np.ndarray:\n",
    "                self.vecs[i] = np.array(v)\n",
    "            if self.vecs[i].shape != (self.dim,):\n",
    "                raise ValueError(f\"Got {self.dim} vectors, therefore all vectors must be {self.dim} dimensional but vector {i+1} has shape {self.vecs[i].shape}\")\n",
    "        self.vecs = np.array(self.vecs)\n",
    "        self.vec_lengths = np.array([np.linalg.norm(v) for v in self.vecs])\n",
    "        self.center = np.zeros(self.dim)\n",
    "\n",
    "    def get_point(self, *ns):\n",
    "        if len(ns) != len(self.vecs): raise ValueError(f\"Requires one index for each lattice vector {len(self.vecs)}, but got only {ns}\")\n",
    "        point = self.center.copy()\n",
    "        for i, n in enumerate(ns):\n",
    "            point += n * self.vecs[i]\n",
    "        return point\n",
    "\n",
    "    \n",
    "    def get_points_around_center(self, n):\n",
    "        import itertools\n",
    "        ns = [i for i in range(-i, i+1)]\n",
    "        for n in itertools.product(*[])\n",
    "        for i in range(self.dim):\n",
    "            \n",
    "            \n",
    "        for i in range(-n, n+1):\n",
    "            for j "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2855d08e-70d8-4ef7-ba19-2c06316caf15",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<scipy.spatial._qhull.Voronoi object at 0x7e3efffb4a70>\n",
      "1.7320508075688774 [0.70710678 1.58113883 0.70710678 0.70710678 1.58113883 0.70710678]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABHkUlEQVR4nO3deVwV9f7H8dc57KLgDogoSkQapLhiZWYmqaWmlpaooFi3rO61bovaYt5+ZlZ2rSzLsnLLpbTSbLMuaV2tLFfUDBeEkCXUgBAVOd/fHyQ3xNwS5nB4Px+PecDMd86cDwzfc97MfGeOzRhjEBEREbGA3eoCREREpOZSEBERERHLKIiIiIiIZRRERERExDIKIiIiImIZBRERERGxjIKIiIiIWEZBRERERCzjbnUBp+NwONi/fz916tTBZrNZXY6IiIicBWMMBQUFNGnSBLv99Mc8nDqI7N+/n5CQEKvLEBERkfOQnp5O06ZNT7uOUweROnXqAKU/iJ+fn8XViJS3a9cuFi5cyMKFC8nIyCAsLIxXX32Vjh07Wl2aiIil8vPzCQkJKXsfPx2bM3/WTH5+Pv7+/uTl5SmIiNNyOBwkJSUxZ84cpk6dSlBQEK+99hre3t4MHDgQX19fq0sUEalS5/L+rcGqIn+R3W6nR48ezJ07l6CgIAA+//xzRowYQWBgIImJiaxZswYnzvwiIpZREBGpBIsXL2bPnj3cf//9JCUl0a1bNzZv3gxAUVGRxdWJiDgPnZoRqWQOh4P169fTqVMnbDYb0dHR+Pn5kZCQwE033XRW51BFRKoTnZoRcSJ2u53OnTtjs9lwOBzcf//9eHp6kpiYSGBgICNGjCA3N9fqMkVELKEgIlKF7HY7cXFxrFq1in379vHwww+zZ88e/P39AZgzZw4//fSTxVWKiFQdnZoRcRKHDx+mefPm5Obm0qVLF+Lj4xkyZAh169a1ujQRkXOiUzMi1VCtWrVIT09n0aJF1K1blzFjxhAaGlo2uNWJ/2cQETlvCiIiTsTb25shQ4bw0UcfkZ6ezltvvYWPjw9HjhwhKiqKBx98kG3btlldpojIBaMgIuKkmjRpwo033ghAYWEh3bt3Z/bs2URGRtKxY0deeuklHSURkWpPQUSkGmjQoAEvvvgi+/fvZ+nSpQQFBbF48WJsNhvGGD799FOKi4utLlNE5JxpsKpINXX8+HHc3d354Ycf6NChA40bNyYuLo74+HjatGljdXkiUoNpsKpIDeDuXvqZle3bt2fTpk3ExcUxf/582rZtS79+/SyuTkTk7CiIiLiANm3a8Nxzz5GRkcHy5csZOnQoAPv27aN///4sW7aMY8eOWVyliEhF7lYXICIXjoeHB3379i2bz8nJITMzk0GDBlG/fn2GDh3KyJEjadeunYVVioj8j46IiLiwjh078t1337Ft2zZGjx7N0qVLeeGFF4DSG6hlZmZaXKGI1HQKIiI1QOvWrZk6dSppaWn8+9//Bko/Ibhp06b06dOHxYsXc+TIEYurFJGaSEFEpAZxd3enXr16AAwYMICZM2eSl5fHLbfcQmBgYNnREhGRqqIgIlJD1a1bl9tvv53//ve/7Ny5k7vvvpvmzZsDsHbtWp588knS09MtrlJEXJ3uIyIiFbzyyivcd999HDlyhB49ehAfH8+AAQPw9fW1ujQRqQZ0HxER+UvuuOMOsrOzmT17NseOHWP48OG8+OKLABQVFenW8iJyweiIiIic0Z49e/Dz86Nhw4aMHz+eJUuWMGLECEaMGEGLFi2sLk9EnIyOiIjIBdWyZUsaNmwIQP/+/enWrRvPPvssLVu25Oqrr2bt2rUWVygi1ZWCiIick5iYGN544w2ysrKYN28eHh4eZbeb/+KLL/jPf/6Dw+GwuEoRqS4qNYhMmTKFjh07UqdOHRo3bsyNN97Izp07K/MpRaSK+Pr6MmzYMFatWkWnTp0AePXVV+nRowctWrTgkUceISUlxeIqRcTZVWoQWb16NXfddRfffPMNq1at4vjx48TGxlJYWFiZTysiFlm8eDH//e9/6dWrFzNmzODiiy/miy++ANAAV5ELJDOviLW7c8nMK7K6lAuiSger/vLLLzRu3JjVq1dz1VVXnXF9DVYVqb6KiopYsWIF/fv3x8vLi7i4OBwOBwkJCVx77bW4ublZXaJItbN4fRrjl23FYcBugykDoxjSsZnVZVXgtINV8/LyAKhfv/4p248ePUp+fn65SUSqJx8fHwYPHoyXlxdQ+rk3W7ZsoVevXjRr1oyHHnqInJwci6sUqT4y84rKQgiAw8CEZcnV/shIlQURYwz33XcfV155JZGRkadcZ8qUKfj7+5dNISEhVVWeiFSysWPHkpyczHfffceAAQOYO3cuNpsNgK+++ooDBw5YXKGIc9ubW1gWQk4oMYbU3MPWFHSBVNmpmbvuuouVK1fy9ddf07Rp01Ouc/ToUY4ePVo2n5+fT0hIiE7NiLigkpIS3NzcKCkpITg4mIMHD9KvXz8SEhK47rrr8PDwsLpEEaeSmVfEFU/9p1wYcbPZ+Hpcd4L8fawr7BSc7tTMPffcw/Lly0lKSvrTEALg5eWFn59fuUlEXNOJMSJubm5s3ryZqVOnsmvXLvr27UvTpk3Jzc21uEIR5xLk78OUgVG4/X4k0c1m48mBkU4XQs5VpR4RMcZwzz338N577/Hll18SHh5+To/XYFWRmmfTpk2sWrWKBx54AGMMQ4YM4fLLL2fo0KE0btzY6vJELJeZV0Rq7mFCG9Zy2hByLu/flRpExowZw9tvv80HH3xARERE2XJ/f398fM78y1MQEanZCgoKSEhIYMWKFRhj6NOnDwkJCfTt27fsJmoi4nycJoicGIh2sjfffJOEhIQzPl5BREQAcnNzWbRoEXPmzGH//v2kpaXh5ubGTz/9RHh4+J++1oiINZwmiPxVCiIicrIDBw7QoEEDfv75Z5o1a0br1q1JSEggLi6OoKAgq8sTEZxwsKqIyIXSoEEDAAIDA/noo4+IjIzkkUceoWnTptx00026g6tINaMgIiLVkru7O7169WLRokVkZWXx8ssv06ZNG2w2G4WFhdx999188803CiYiTk6nZkTE5WzevJm+ffuSnp7OxRdfTEJCAsOHDz/t7QNE5MLRqRkRqdHatGlDamoqn3/+OZ06deKJJ54gLi4OAIfDweHD1ftOlCKuREFERFyS3W6nR48ezJs3j6ysLF577TUA1qxZQ2BgIImJiaxZs0anbkQspiAiIi7Pz8+Piy++GIDw8HD++c9/kpSURLdu3QgLC2PmzJkWVyhScymIiEiNEhwczMSJE9m1axerV6/m6quvprCwEIDU1FTefPNNCgoKLK5SpObQYFURkd+98cYbjB49Gh8fHwYNGkR8fDzdu3fHbtf/bCLnQoNVRUTOw6hRo0hNTWXChAl88803XHvttTz44IMAGksiUkl0RERE5BSMMaxbt44GDRoQERHBrFmzmDNnDgkJCQwePBh/f3+rSxRxWjoiIiLyF9lsNi6//PKyD+wMCwvDz8+PO+64g8DAQG699VbWr19vcZUi1Z+CiIjIWejRowcff/wx6enpTJo0iS1btrBnzx4Atm/fzo4dOyyuUKR60qkZEZHzYIzBGIPdbichIYE5c+bQsWNHEhISuOWWW6hfv77VJYpYRqdmREQqmc1mK7ua5tVXX+Xdd98lMDCQv//97wQFBfHBBx9YXKFI9eBudQEiItWdl5cXgwYNYtCgQWRnZ/P222/TuXNnACZPnszBgweJj4/nsssus7hSEeejIyIiIhdQQEAA9957L4GBgUDpZ9vMmzePNm3aEB0dzfTp0zl48KDFVYo4DwUREZFK9Oijj5KRkcEHH3xAy5YtefDBB0lLSwNgz549HDt2zOIKRaylwaoiIlXo0KFD1KtXDyj9lOCMjAyGDh1KfHw87dq1w2azWVyhyF+nwaoiIk7qRAgBePvttxk1ahTvvPMOHTp0ICoqir1791pYnUjVUxAREbHIpZdeytNPP016ejofffQRMTExNG3aFIB//etfLFmyhCNHjlhcpUjl0qkZEREnc/z4ca655hq++uor6tatyy233EJ8fDydO3fWqRupFnRqRkSkGnN3d2fNmjX8+OOPjBkzhg8//JBu3brx66+/AlBQUGBtgSIXkIKIiIiTioiIYPLkyaSmprJ+/Xrq1avHkSNHCA0NJTY2lgULFnD48GGryxT5SxREREScnJubW7mboT3zzDMcOXKEYcOGERgYyOjRoykuLrawQpHzpyAiIlKNeHt7M2rUKNasWcPu3bu57777+O233/Dw8MAYw/Tp00lNTbW6TJGzpsGqIiIuYvfu3bRp04bCwkKuvvpq4uPjuemmm6hdu7bVpUkNo8GqIiI1UFhYGNnZ2cydOxe73c6oUaPo3LkzJ/7fdOL/O6UGUxAREXEhvr6+DB8+nC+++ILU1FRefvllbDYbqamptGjRgkcffZRdu3adfiOOEtj7FWx9t/Sro6RqipcaSUFERMRFNWvWjG7dupXNx8bG8sILLxAeHs6VV17J/PnzKz5o+3KYHglzboCliaVfp0eWLhepBAoiIiI1QGhoKLNmzSIrK4uFCxdSu3ZtvvzyS6D0viSffvopJVvfgyUjIH9/+QfnZ5YuVxiRSqDBqiIiNZTD4cBut7N48WJuueUWmvi5MzzKTnwbD1o1cjtpbRv4NYGxW8F+cptIeefy/u1eRTWJiIiTsdtLD4oPHjyYFl6/MmfyPcz64RhT/3uMezp5Mv5KTzJ/+8P/qplp1Pv6XVpcNYQjR46wffv2Ctts164dADt37qSwsLBcW2hoKPXr1+eXX34hPT29XFudOnUIDw+npKSEzZs3V9huVFQUHh4e7N69m7y8vHJtwcHBBAQEcOjQoQofGujj40OrVq0A2LhxY4UBu61atcLHx4d9+/Zx4MCBcm0BAQEEBwdTUFBASkpKuTYPDw+ioqIA2Lp1a4X7uISHh1OnTh0yMjLIzs4u19agQQOaN29OUVERO3bsKNdms9mIjo4GYMeOHRQVFZVrb9GiBfXq1SM7O5uMjIxybf7+/oSFhVFcXMzWrVs5WZs2bXBzcyMlJaXC3XlDQkJo1KgRBw8erHD5t6+vLxEREQBs2LChwnZbt26Nt7c3e/fu5dChQwD89ttvFdb7U8aJ5eXlGcDk5eVZXYqIiGvb8o4xE/3MkYfrmHdu9jHv3ORtujW3G6DcFHd9V2OMMSkpKRXa/viWEhMTU6Ft3rx5xhhjZsyYUaEtNjbWGPO/1/2Tp5ycHGOMMX379q3QNm3aNGOMMUuWLKnQFh0dXVaTp6dnhfbk5GRjjDGJiYkV2saNG2eMMSYpKalCW3BwcNl2g4ODK7QnJSUZY4wZN25chbbExERjjDHJyckV2jw9Pcu2Gx0dXaF9yZIlxhhjpk2bVqGtb9++xhhjcnJyTvk7PPFeGhsbW6FtxowZxhhj5s2bV6EtJiamrKZTbTclJcUYY0xcXNyfPufp6NSMiIiUXh0z54ay2Q2ZJbSfVcj8gd60avi/UzH14mbriAg6IvJHf3ZEpFu3bmf1/q0gIiIipZfoTo8sHZiKKQsiP9zuS7sgNzRGRM6FbmgmIiLnxu4Gvab+PmM7qfH3+V5PKYTIBacgIiIipVr3g8FzwS+o/HK/JqXLW/ezpi5xaQoiIiLyP637wdhkuGF66fwN00tPxyiESCVREBERkfLsbhBcOuiU4HY6HSOVSkFERERELKMgIiIiIpZREBERERHLKIiIiIiIZRRERERExDIKIiIiImIZBRERERGxjIKIiIiIWEZBRERERCyjICIiIiKWURARERERyyiIiIiIiGUURERERMQyCiIiIiJiGQURERERsYyCiIiIiFim2gSRhIQEnnvuObKzs60uRURERC6QahFEjh49SmFhIePHjyc4OJi+ffvy7rvvUlxcbHVpIiIi8hdUiyDi5eXFO++8Q2ZmJi+88AI5OTncfvvtOBwOAHbv3o0xxuIqRURE5FxViyByQv369RkzZgzffvstO3fuxMvLiwMHDtCqVSsuvfRSpk6dSkZGhtVlioiIyFmqVkHkjxo1agRA3bp1+fDDD2nbti2PP/44zZo1o2/fvmVHS0RERMR5VdsgcoKbmxuxsbG8/fbbZGVl8eqrrxIVFYXdbufYsWP84x//YO3atTp1IyIi4oTcrS7gQvL392f06NFl83v37uX999/nhRdeIDw8nBEjRjBixAiaNWtmYZUiIiJyQrU/InI6ERER7N27ly+++IIuXbowZcoUBgwYAIAxhsOHD1tcoYiISM3m0kEEwG63c8011zBnzhyys7OZP38+ABs3biQgIICRI0eyevVqjSkRERGxgMsHkT+qXbs2rVq1AqBJkyY88MADrFmzhquvvpqwsDCef/55iysUERGpWWpUEPmjwMBAHnvsMXbt2sWaNWvo0aMHv/76KwBZWVm88cYbFBQUWFukiIiIi6vUILJmzRr69u1LkyZNsNlsvP/++5X5dOfFZrPRtWtXXn/9dSZOnAjAl19+yejRowkICGD48OF8/vnnlJSUWFypiIiI66nUIFJYWEibNm2YMWNGZT7NBXfLLbewb98+Hn30UdavX0/Pnj256667AHQZsMgFkplXxNrduWTmFVldiki14mp9p1Iv3+3duze9e/euzKeoNCEhIYwfP55x48bx7bff4uvrC8DixYt54YUXiI+PZ8iQIdStW9faQkWqocXr0xi/bCsOA3YbTBkYxZCOuqxe5Excse841RiRo0ePkp+fX26yms1mIyYmhqioKKB0kGu9evUYM2YMgYGBDBkyhHXr1llcpUj1kZlXVPZCCuAwMGFZssv8dydSWVy17zhVEJkyZQr+/v5lU0hIiNUlVXDVVVexcuVKfv75Z/7v//6P7du3k5ycDMCuXbvKvheRU9ubW1j2QnpCiTGk5uq+PiKn46p9x6mCyPjx48nLyyub0tPTrS7pTwUFBXH//fezZcsWRo0aBcD06dOJioqiQ4cOvPjii+Tm5lpcpYjzadHQF7ut/DI3m43QhrWsKUikmnDVvuNUQcTLyws/P79yk7Oz2Wy4ubkB8Nxzz7Fs2TKaNm3KfffdR5MmTVi4cKHFFYo4lyB/H6YMjMLNVvqK6maz8eTASIL8fSyuTMS5uWrfcanPmrGap6cnAwYMYMCAAeTk5LBw4UJiYmIAeP7559m7dy8JCQm0bdvW2kJFLDakYzOuurgRqbmHCW1Yq9q/kIpUFVfsO5UaRH777Td27dpVNr937142bdpE/fr1Xf6D5xo3bsw//vGPsvljx46xcOFCnn/+eS677DISEhIYPnw4DRs2tLBKEesE+fu4xIuoSFVztb5Tqadmvv/+e6Kjo4mOjgbgvvvuIzo6mscee6wyn9YpPfDAA/z888+sWLGC8PBwxo0bx44dOwDYt28fR48etbhCERGRqlepR0Suvvpq3QDsDzw8PLjhhhu44YYbOHjwIPXq1QNg2LBhbNu2jVtvvZWEhAQ6dOiAzWY7w9ZERESqP6carFqT1K9fvyxszJo1i9tvv53333+fTp06ERkZWXa0RERExJUpiDiBVq1a8dRTT5GWlsYnn3xCx44dad68OQBPP/00ixYtoqioet+wRkRE5FQURJyIm5sb1113HW+99Ra1atXCGMNnn33GrbfeSlBQEH/7299Yt26dTneJiIjLUBBxYjabjc8//5yffvqJe+65h48//pgrr7ySrKwsoPSqJBERkepMQaQaCA8P54knniA1NZXvv/+eoKAgjh8/TkREBD169GDu3LkUFhZaXaaIiMg5UxCpRux2e9ml0A6Hg8mTJ+NwOIiPjycwMJCRI0dy+HD1/swBERGpWRREqilPT08SEhJISkpiz549PPDAA+Tm5uLjU3qTm5deeondu3dbXKWIiMjpKYi4gBYtWvDYY4+xYsUKbDYbmZmZjB8/nosuuoirrrqKN954g/z8fKvLFBERqUBBxAUFBQWRlZXF/Pnz8fLyYvTo0Vx22WVlV9voqhsREXEWCiIuqlatWsTFxbFq1Sr27dvHa6+9hs1mIycnh7CwMB5++GF++uknq8sUEZEaTkGkBggJCaFnz54AFBcXc9111/Hyyy8TERFBly5deOONNyyuUEREaioFkRomODiYmTNnkpmZyeLFi6lfvz6ffvopAEePHuXjjz/m+PHjFlcpIiI1RaV+6J04L29vbwYPHszgwYNxOBwAfPHFF1x//fUEBgYyfPhw4uPjufTSSy2uVEREXJmOiAh2e+mfQe/evfn++++5+eabeeONN4iMjGTUqFEWVyciIq5MQUTK2Gw22rdvzwsvvMD+/ftZtmwZgwcPBuC7775j0KBBLF++nOLiYosrFRERV6FTM3JKnp6eDBgwoGz+t99+Y+/evfTv359GjRoRFxdHYmIikZGRFlYpIiLVnY6IyFm55ppr2LBhA5s3b2b48OG8/fbbLFy4EICDBw+Sk5NjcYUiIlIdKYjIObnsssuYNm0aP//8M+PGjQNg5syZNGnShH79+rF06VKOHj1qcZUiIlJdKIjIefHw8KBOnToA3HHHHTz//PNkZWVx00030aRJE+bNm2dxhSIiUh0oiMhf1qBBA+666y6+++47tm3bxm233UbLli0BWLFiBU8//TT79++3uEoREXFGCiJyQbVu3ZqnnnqKK664AoCtW7cyceJEQkJC6N27N4sWLaKoqMjiKkVExFkoiEilmjBhApmZmcycOZP8/HxuvfVW3nnnHaD0Shx9AJ+ISM2mICKVrm7dutx+++3897//5aeffmLQoEEAjBkzhksuuYTJkyeTlpZmcZUiImIFBRGpUuHh4fj6+gKQmJhI586defLJJwkNDeXaa69l69atFlcoIiJVSUFELNOtWzfmzp1LVlYWs2fPpqSkBH9/fwCWLVvGmjVryj4HR0REXJOCiFiuTp06jBw5kqSkJJo1awbA888/T7du3bjooot4/PHH2bNnj8VViohIZVAQEaeUlJTE6tWr6d69O8899xxhYWEkJycDaICriIgLURARp2S327nqqquYPXs2WVlZvPvuu1x66aUA9OjRg+HDh/P555+f/tSNowT2fgVb3y396iipoupFqjlHCWRsKP0+Y4P6jlQqBRFxerVq1WLQoEHYbDYcDgfXXnst3333HT179iQ0NJSHH36YvLy88g/avhymR8KcG2BpYunX6ZGly0Xkz53oOx+OLZ3/cKz6jlQqBRGpVux2OxMmTODHH39k3bp19OnTh0WLFuHt7Q3Ap59+yq/fvA1LRkD+SXdzzc8sXa4XVJFT275cfUeqnLvVBYicD5vNRkxMDDExMTgcDux2O7/99hsDBgzAUXyEGyPcSGjrSc+WbrjZbb8/ygA2+GQcXHI92N2s/BFEnIujBD55iNJ+cjL1Hak8CiJS7dntpQf2ateuza7/LGDBo7fy1uZiei84TDN/G6vja3HwyB8ekJkGK9+E4HbWFCzijDI2wM70stkduSXlvgKQmUa9r9+lxVVDOHLkCNu3b6+wmXbtSvvVzp07KSwsLNcWGhpK/fr1+eWXX0hPTy/XVqdOHcLDwykpKWHz5s0VthsVFYWHhwe7d++ucCo2ODiYgIAADh06xN69e8u1+fj40KpVKwA2btxYYbB7q1at8PHxYd++fRw4cKBcW0BAAMHBwRQUFJCSklKuzcPDg6ioKKD0oyyKi4vLtYeHh1OnTh0yMjLIzs4u19agQQOaN29OUVERO3bsKNdms9mIjo4GYMeOHRU+EqNFixbUq1eP7OxsMjIyyrX5+/sTFhZGcXHxKe/J1KZNG9zc3EhJSaGgoKBcW0hICI0aNeLgwYOkpqaWa/P19SUiIgKADRs2VNhu69at8fb2Zu/evRw6dAgovXP2WTNOLC8vzwAmLy/P6lKkutjyjjET/YzjsTpm/W2+Zuq1XqZbc7uh9F86TZo0ncPk71VxWdz1XY0xxqSkpJzyMSfExMRUaJs3b54xxpgZM2ZUaIuNjTXG/O91/+QpJyfHGGNM3759K7RNmzbNGGPMkiVLKrRFR0eX1eTp6VmhPTk52RhjTGJiYoW2cePGGWOMSUpKqtAWHBxctt3g4OAK7UlJScYYY8aNG1ehLTEx0RhjTHJycoU2T0/Psu1GR0dXaF+yZIkxxphp06ZVaOvbt68xxpicnJxT/g5PvJfGxsZWaJsxY4Yxxph58+ZVaIuJiSmr6VTbTUlJMcYYExcX96fPeTq23zfslPLz8/H39ycvLw8/Pz+ry5HqYO9XpQNTf7chs4T2swqZP9CbVg3/cDj5huk6IiLyRxkb/jdAldIjIcOWHanQd+rFzdYREXRE5I/+7IhIt27dzur9W0FEXIujpHSEf34mYMqCyA+3+9IuyA2wgV8TGLtV57lF/kh9Ry6gc3n/1lUz4lrsbtBr6u8ztpMaf5/v9ZReSEVOpr4jFlEQEdfTuh8Mngt+QeWX+zUpXd66nzV1iTg79R2xgIKIuKbW/WBsculYECj9OnarXkhFzkR9R6qYgoi4Lrvb/wakBrfTIWWRs6W+I1VIQUREREQsoyAiIiIillEQEREREcsoiIiIiIhlFERERETEMgoiIiIiYhkFEREREbGMgoiIiIhYRkFERERELKMgIiIiIpZREBERERHLKIiIiIiIZRRERERExDIKIiIiImIZBRERERGxjIKIiIiIWEZBRERERCyjICIiIiKWURARERERyyiIiIiIiGUURERERMQyCiIiIiJiGQURERERsYyCiIiIiFhGQUREREQsoyAiIiIillEQEREREctUSRB5+eWXadGiBd7e3rRv356vvvqqKp5WREREnFylB5HFixczduxYHn74YTZu3EjXrl3p3bs3aWlplf3UVSYzr4i1u3PJzCuyuhSRakP9RuT8uFrfca/sJ3juuedITExk9OjRAEyfPp1PP/2UmTNnMmXKlMp++kq3eH0a45dtxWHAboMpA6MY0rGZ1WWJODX1G5Hz44p9p1KPiBw7dowffviB2NjYcstjY2NZu3ZthfWPHj1Kfn5+ucmZZeYVlf1BADgMTFiW7DIpVaQyqN+InB9X7TuVGkRyc3MpKSkhICCg3PKAgACysrIqrD9lyhT8/f3LppCQkMos7y/bm1tY9gdxQokxpOYetqYgkWpA/Ubk/Lhq36mSwao2m63cvDGmwjKA8ePHk5eXVzalp6dXRXnnrUVDX+wn/RhuNhuhDWtZU5BINaB+I3J+XLXvVGoQadiwIW5ubhWOfuTk5FQ4SgLg5eWFn59fucmZBfn7MGVgFG6/hyo3m40nB0YS5O9jcWUizkv9RuT8uGrfqdTBqp6enrRv355Vq1YxYMCAsuWrVq2if//+lfnUVWZIx2ZcdXEjUnMPE9qwVrX/gxCpCuo3IufHFftOpV81c9999zF8+HA6dOhAly5dmDVrFmlpadxxxx2V/dRVJsjfxyX+GESqkvqNyPlxtb5T6UFkyJAhHDhwgH/9619kZmYSGRnJRx99RPPmzSv7qUVERMTJVXoQARgzZgxjxoypiqcSERGRakSfNSMiIiKWURARERERyyiIiIiIiGUURERERMQyCiIiIiJiGQURERERsYyCiIiIiFhGQUREREQsoyAiIiIillEQEREREcsoiIiIiIhlFERERETEMgoiIiIiYhkFEREREbGMgoiIiIhYRkFERERELKMgIiIiIpZREBERERHLKIiIiIiIZRRERERExDIKIiIiImIZBRERERGxjIKIiIiIWEZBRERERCyjICIiIiKWURARERERyyiIiIiIiGUURMR1OUogY0Pp9xkbSudF5MzUd6QKKYiIa9q+HKZHwodjS+c/HFs6v325lVWJOD/1HaliCiLierYvhyUjIH9/+eX5maXL9YIqcmrqO2IBBRFxLY4S+OQhwJyi8fdln4zToWaRk6nviEXcrS7gbGzatInatWuXzderV48WLVpw5MgRtm/fXmH9du3aAbBz504KCwvLtYWGhlK/fn1++eUX0tPTy7XVqVOH8PBwSkpK2Lx5c4XtRkVF4eHhwe7du8nLyyvXFhwcTEBAAIcOHWLv3r3l2nx8fGjVqhUAGzduxJjyHb1Vq1b4+Piwb98+Dhw4UK4tICCA4OBgCgoKSElJKdfm4eFBVFQUAFu3bqW4uLhce3h4OHXq1CEjI4Ps7OxybQ0aNKB58+YUFRWxY8eOcm02m43o6GgAduzYQVFRUbn2Fi1aUK9ePbKzs8nIyCjX5u/vT1hYGMXFxWzdupWTtWnTBjc3N1JSUigoKCjXFhISQqNGjTh48CCpqanl2nx9fYmIiABgw4YNFbbbunVrvL292fv1uxza+b/9uiO3pNxXADLTYOWbENyuwnZEaqyMDfCHvrPip9LXk7wjJ16vDORnwL610KKrBQWKyzJOLC8vz1AaxctNcXFxxhhjUlJSTtl+QkxMTIW2efPmGWOMmTFjRoW22NjY0z5vTk6OMcaYvn37VmibNm2aMcaYJUuWVGiLjo4uq8nT07NCe3JysjHGmMTExApt48aNM8YYk5SUVKEtODi4bLvBwcEV2pOSkowxxowbN65CW2JiojHGmOTk5Aptnp6eZduNjo6u0L5kyRJjjDHTpk2r0Na3b19jjDE5OTmn/B3m5eUZY4yJjY2t0DZjxgxjjDHz5s2r0BYTE1NW06m2m5KSYowxJu76rhXa/L0qrq9Jk6azmzzdMDe1djcrbvUxxY/WMWbLO+fzci41zIn30ROv+adjM+akf8+dSH5+Pv7+/qxevVpHRHRE5OyOiKxZzKEFiWXLd+SWMGzZEeYP9KZVQ7f/PeCG6ToiIvJHGRv+N0D1d+42wxd7S3hrczFbsh0829OLf766iuKmMXh4eFhTp1QLJ96/8/Ly8PPzO+261SKInM0PIgKUnr+eHlk6uA7DhswS2s8q5IfbfWkX5AbYwK8JjN0KdrczbU2k5jip75xsU5aDJkFBNH50B+MmPMxnn31GfHw8Q4cOpVGjRlVfrzi1c3n/1mBVcS12N+g19fcZ20mNv8/3ekohRORkZ+g7bQPdaHzzs2B3o0ePHoSGhvLAAw/QpEkTbrzxRjZu3FjVFYuLUBAR19O6HwyeC35B5Zf7NSld3rqfNXWJOLuz7Ds9e/Zk2bJl7N+/n3//+9/8/PPPZaeG16xZw4YNGyqcghb5Mzo1I67LUcKGlW/Svt9t/LD8NdpdP1JHQkTOhqOk9OqY37KhdgA0v/ys+06fPn34+OOPiYyMJD4+nri4OIKCgs78QHEpOjUjAqUvnCcGpAa3UwgROVt2t9JLdKNuKv16Dn1n+fLlrFy5ktatW/PII4/QtGlTvv7660osVqq7anEfERERqR7c3d3p06cPffr04dChQ7z77rt07NgRgFGjRuHl5UV8fDydO3fGZjt5LIrURDoiIiIilaJevXrcdttteHl5AdCsWTM+/PBDunTpQqtWrZgyZQqHDh2yuEqxmoKIiIhUiccff5zU1FQ+++wzOnTowNSpUzl+/DgA33//PYcPH7a4QrGCgoiIiFQZNzc3evbsyfz588nOzqZRo0YcP36c3r17ExgYyOjRo/n666911U0NoiAiIiKWOHHKxt3dnW+++YZ7772Xzz//nK5duxIeHq7TNjWEgoiIiFguLCyMSZMmsWfPHpKSkoiLi6NevXoYYxg5ciRvvfUWv/32m9VlSiVQEBEREadht9u5+uqrmTRpEgB5eXns27ePkSNHEhgYSHx8PElJSTgcDosrlQtFQURERJxW3bp1+c9//kNqairjxo3jv//9L8OGDSsbQ7J//36LK5S/SkFEREScXvPmzXnkkUdISUnhm2++wc3NjfT0dEJCQrjyyit5/fXXK3wqulQPCiIiIlJt2Gw2QkJCAGjUqBELFiygdu3a/O1vfyMwMJDbbrvN4grlXCmIiIhIteTt7c0tt9zCJ598QlpaGo8//jiNGzcGoLCwkIcffpgdO3ZYXKWciYKIiIhUe8HBwTz00ENMnjwZgG3btjFz5kxat25N586dmTlzpi4HdlIKIiIi4nI6depEZmYm77zzDo0bN+aee+7h5ptvBsDhcJTd0VWspyAiIiIuycvLi5tuuokVK1bw888/M336dADWrFlD06ZN+ec//8mWLVusLVIURERExPUFBgYSGRkJQEhICLfccgtz586lTZs2tGvXjgULFlhcYc2lICIiIjVKWFgY06dPJyMjg/fff5/Q0FDS09MB2LdvH++99x7Hjh2zuMqaw93qAkRERKzg6elJ//796d+/f9myFStWcM8999CgQQOGDh1KfHw87dq1w2azWVipa9MRERERkd/dfffdbN26lZEjR/LOO+/QoUMHJk6caHVZLk1BRERE5A8iIyN55plnSE9PZ+XKlQwZMgSA2bNnc/311/POO+9w5MgRi6t0HQoiIiIip+Du7k6fPn249NJLgdI7uR48eJDBgwcTFBTEnXfeybZt2yyusvpTEBERETkL/fr1Y926dfz444/ceeedfPjhhyQnJwOQkpLCzz//bHGF1ZOCiIiIyDmIiIjgySefJDU1lUGDBgEwceJEmjVrRmxsLAsWLODw4cMWV1l9KIiIiIicBzc3N9zdSy8+feWVV3jttdc4cuQIw4YNIzAwkE8//dTiCqsHXb4rIiLyF/n5+ZGYmEhiYiK7d+8uu1kawLPPPktRUREjRoygefPmFlfqfHRERERE5AIKCwtj0qRJBAYGApCdnc3UqVMJDQ2le/fuzJkzh8LCQourdB4KIiIiIpXomWeeISsrizlz5mC320lISGDXrl0A7N+/H4fDYXGF1lIQERERqWS1a9dmxIgRfPHFF/z888+0adMGYwy9evWiZcuWPPbYY2XhpKap1CAyefJkLr/8cmrVqkXdunUr86lERESqheDg4LLvZ86cSWxsLM8//zzh4eFceeWVNe4y4EoNIseOHePmm2/mzjvvrMynsVxmXhFrd+eSmVdkdSki1Yb6jdR0NpuNK664glmzZpGVlcXbb79NcHAwAQEBQOkg188++4ySkhKLK61clXrVzKRJkwB46623KvNpLLV4fRrjl23FYcBugykDoxjSsZnVZYk4NfUbkfJ8fHy49dZbufXWW4HSf+TnzJlDcnIywcHBDB8+nPj4eC655BKLK73wnGqMyNGjR8nPzy83ObPMvKKyF1MAh4EJy5L1H57IaajfiJyZp6cnW7Zs4dtvv6V///68+uqrtG3btux98dixYxZXeOE4VRCZMmUK/v7+ZVNISIjVJZ3W3tzCshfTE0qMITVXd9QT+TPqNyJnx2az0alTJ1566SUyMzNJSkrCz8+PoqIimjVrxuDBg1m5ciXHjx+3utS/5JyDyOOPP47NZjvt9P33359XMePHjycvL69sSk9PP6/tVJUWDX2x28ovc7PZCG1Yy5qCRKoB9RuRc+fl5UWXLl0AcDgcPPjgg/z444/ccMMNhISE8MADD1TbsSTnHETuvvtuduzYcdopMjLyvIrx8vLCz8+v3OTMgvx9mDIwCjdb6auqm83GkwMjCfL3sbgyEeelfiPy1/j6+nLfffexefNmNmzYwJAhQ/jpp59wc3PDGMPs2bPJzc21usyzds6DVRs2bEjDhg0ro5ZqaUjHZlx1cSNScw8T2rCWXkxFzoL6jchfZ7PZiI6OJjo6umzZzp07ueOOO7jzzju5/vrriY+Pp0+fPnh6elpY6elV6lUzaWlpHDx4kLS0NEpKSti0aRMAF110EbVr167Mp65SQf4+eiEVOUfqNyIX3iWXXML+/ftZuHAhc+bMYcCAAbRv3/68h0xUhUoNIo899hhz5swpmz+R2pKSkrj66qsr86lFRERqpEaNGvH3v/+dv//97yQnJ5ORkQFAamoq/fv3Z8SIEcTFxZV9Fo7VKvWqmbfeegtjTIVJIURERKTyRUZGct111wGlt8iIiIhgwoQJNG3alBtuuIHly5dbXKGTXb4rIiIilSMiIoIlS5aQlZXFjBkzyM3N5b333gOgoKCAb7/9FmPMGbZy4SmIiIiI1CD16tXjjjvu4JtvvmHWrFkAvP/++8TExNCqVSueeuqpstM5VUFBREREpIby8PAAYOjQoXz22We0b9+eSZMmERISwoQJE6qkBgURERGRGs7NzY2ePXuyYMECsrKymDVrFt27dwfgyy+/5LbbbuPrr7+ulFM3CiIiIiJSxt/fn9GjR9OzZ08ADh48yOeff07Xrl25+OKLeeKJJ0hLS7tgz6cgIiIiIn9q4MCB7N69m6SkJK644gqmTp3K22+/DZSGlN9+++0vbV9BRERERE7Lbrdz9dVX89Zbb5GVlcWYMWMAmDp1KoGBgSQkJJCUlITD4Tj3bV/oYkVERMR11a5du+yz4MaMGcO4ceP4+uuvueaaa2jZsiUrV648p+0piIiIiMh5ad68OY888ggpKSl8/fXX9OzZk6CgoHPaRqXe4l1ERERcn81m44orruCKK64AID8//6wfqyMiIiIiYhkFEREREbGMgoiIiIhYRkFERERELKMgIiIiIpZREBERERHLKIiIiIiIZRRERERExDIKIiIiImIZBRERERGxjIKIiIiIWEZBRERERCyjICIiIiKWURARERERyyiIiIiIiGUURERERMQyCiIiIiJiGQURcV2OEsjYUPp9xobSeRERcSoKIuKati+H6ZHw4djS+Q/Hls5vX25lVSIichIFEXE925fDkhGQv7/88vzM0uUKIyIiTkNBRFyLowQ+eQgwp2j8fdkn43SaRkTESbhbXcDZ2LRpE7Vr1y6br1evHi1atODIkSNs3769wvrt2rUDYOfOnRQWFpZrCw0NpX79+vzyyy+kp6eXa6tTpw7h4eGUlJSwefPmCtuNiorCw8OD3bt3k5eXV64tODiYgIAADh06xN69e8u1+fj40KpVKwA2btyIMeXfJFu1aoWPjw/79u3jwIED5doCAgIIDg6moKCAlJSUcm0eHh5ERUUBsHXrVoqLi8u1h4eHU6dOHTIyMsjOzi7X1qBBA5o3b05RURE7duwo12az2YiOjgZgx44dFBUVlWtv0aIF9erVIzs7m4yMjHJt/v7+hIWFUVxczNatWzlZmzZtcHNzIyUlhYKCgnJtISEhNGrUiIMHD5KamlquzdfXl4iICAA2bNhQYbutW7fG29ubvV+/y6Gd/9uvW3OOA7D7oIN2QW6AgfwM2LcWWnStsB0REalixonl5eUZSv+NLTfFxcUZY4xJSUk5ZfsJMTExFdrmzZtnjDFmxowZFdpiY2NP+7w5OTnGGGP69u1boW3atGnGGGOWLFlSoS06OrqsJk9PzwrtycnJxhhjEhMTK7SNGzfOGGNMUlJShbbg4OCy7QYHB1doT0pKMsYYM27cuAptiYmJxhhjkpOTK7R5enqWbTc6OrpC+5IlS4wxxkybNq1CW9++fY0xxuTk5Jzyd5iXl2eMMSY2NrZC24wZM4wxxsybN69CW0xMTFlNp9puSkqKMcaYuOu7Vmhzt5d+7dDEbmb09ja5D9Q2Zss75/MnKSIiZ+HE++iJ1/zTsRlz0r/nTiQ/Px9/f39Wr16tIyI6InJ2R0TWLObQgsRybfW9YUOWgzmbi/ko5Tjtg9z45qskTOiVlJSU4O5eLQ4MiohUGyfev/Py8vDz8zvtutUiiJzNDyIClI79mB5ZOjD1FONEcgoNGSUNiX42hR82bqJPnz7ExcURHx9PmzZtqr5eEREXdC7v3xqsKq7F7ga9pv4+Yzup0UZjXzvRo6eD3Y0GDRoQFxfH/Pnzadu2LW3btuXNN9+s4oJFRGo2BRFxPa37weC54BdUfrlfk9LlrfsBpafpnnvuOTIyMli+fDkXXXQRu3fvBiArK4tly5Zx7Nixqq5eRKRG0akZcV2OktKrY37LhtoB0Pzy0iMmZ+HNN99k1KhR1K9fn6FDhxIfH0/79u2x2U4+yiIiIifTqRkRKA0dLbpC1E2lX88yhACMHDmSbdu2MXr0aJYuXUrHjh25//77K7FYEZGaSUFE5E+0bt2aqVOnkpaWxscff8ywYcMAWLx4Mb1792bRokUVrigSEZFzoyAicgbu7u706tWr7JJmPz8/8vLyuPXWWwkKCuJvf/sbW7ZssbhKEZHqSUFE5Bz17t2btWvXsnPnTu6++24+/vhjvvvuOwBSU1NJS0uzuEIRkepDg1VF/iKHw0FJSQkeHh7cfvvtvP7661xzzTUkJCQwYMAAfH19rS5RRKRKabCqSBWy2+14eHgAMG3aNGbPnk1xcTHDhw8nMDCQ999/39oCRUScmI6IiFSSPXv2MG/ePIYPH07Lli15+eWXycnJYcSIEbRs2dLq8kREKo2OiIg4gZYtWzJx4sSy0JGZmcm0adMICwujW7duvPnmmxU+b0dEpKZREBGpIk888QRZWVnMmzcPDw8PEhMT2bhxI1B6J1eHw2FxhSIiVU+nZkQskp6eTnBwMHa7nSuvvJK0tDRGjBhBfHw84eHhVpcnInLedGpGpBoICQnBbi/tgs8++yx9+vRhxowZXHzxxVx++eVln3sjIuLKFEREnEBMTAyvvPIKmZmZLFq0iMDAQIKCSj+078UXX+STTz6hpKTE4ipFRC48nZoRcWIOh4NOnTrxww8/0KRJE4YNG0Z8fDytW7e2ujQRkT+lUzMiLsJut7N+/XrWr1/PwIEDef3114mKiuKXX34BoLi42OIKRUT+GgURESdns9no0KEDL774Ivv37+fLL7+kUaNGHD9+nPDwcAYNGsSKFSsUSkSkWlIQEalGvLy86Nq1KwDHjx9n7Nix7N69m379+tG0aVPuvfdejh49anGVIiJnT0FEpJry9vZm7NixbNq0iY0bNzJ06FC2bNmCp6cnAHPmzCEnJ8fiKkVETk+DVUVcUFpaGhdddBHGGHr37k1CQgLXX389Xl5eVpcmIjWABquK1HDNmjUjMzOT6dOnk5mZyaBBg2jbti1O/H+HiNRQCiIiLqpBgwbcddddrF+/nuTkZKZOnYrNZiM7O5u2bdvyzDPPkJmZaXWZIlLDKYiI1ACXXnop/fr1A+Dw4cO0atWKRx99lKZNm9K7d2+WLl1qcYUiUlMpiIjUMC1atGDhwoVkZWUxc+ZM8vPzWbRoEQBHjhxh3bp1OoUjIlVGg1VFhGPHjuHp6cl7773HwIEDufjii4mPj2f48OGEhIRYXZ6IVDMarCoi5+TEJb/9+/fn888/p3PnzkyePJnmzZszduxYa4sTEZemICIiZex2Oz169GDu3LlkZWUxe/ZsrrnmGgDWr19PYmIia9as0akbEblgFERE5JTq1KnDyJEjywa55uTkkJSURLdu3QgLC2PSpEns3bvX4ipFpLpTEBGRs3L99deza9cuVq9eTffu3Xn22Wd5/fXXAcjLy6OgoMDiCkWkOlIQEZGzZrfbueqqq5g9ezZZWVncf//9AMyYMYPAwEBGjBjBF198gcPhsLhSEakuKi2IpKamkpiYSIsWLfDx8SEsLIyJEydy7NixynpKEalCvr6+1KtXD4ARI0bw8MMP880333DttdcSGhqqe5OIyFmptCDy448/4nA4ePXVV9m2bRv//ve/eeWVV5gwYUJlPaVlMvOKWLs7l8y8IqtLEbFESEgIEyZMYOfOnaxdu5bevXvTuHFjAD755BNmzZrFr7/+am2RIuKUqvQ+Is888wwzZ85kz549Z7V+dbiPyOL1aYxfthWHAbsNpgyMYkjHZlaXJeI0JkyYwNSpU/H09OTGG28kPj6enj174ubmZnVpIlJJnPY+Inl5edSvX/9P248ePUp+fn65yZll5hWVhRAAh4EJy5J1ZETkD5588knS09OZNGkSW7ZsoXfv3ixevBhAp2pFpOqCyO7du3nxxRe54447/nSdKVOm4O/vXzY5+x0d9+YWloWQE0qMITX3sDUFiTipJk2a8OCDD5KcnMz69evp378/AKNHj6ZTp0689NJLHDx40OIqRcQK5xxEHn/8cWw222mn77//vtxj9u/fT69evbj55psZPXr0n257/Pjx5OXllU3p6enn/hNVoRYNfbHbyi9zs9kIbVjLmoJEnJzNZqNDhw74+voCcPPNNxMYGMg//vEPgoKCuOmmm/jpp58srlJEqtI5jxHJzc0lNzf3tOuEhobi7e0NlIaQ7t2707lzZ9566y3s9rPPPtVljMiEZcmUGIObzcaTAyM1RkTkHGVnZ/P2228zb948PvjgA0JCQvjoo49o2rQpl112mdXlicg5Opf370odrJqRkUH37t1p37498+fPP+fBadUhiEDpWJHU3MOENqxFkL+P1eWIuIT27duzYcMG2rZtS0JCAkOHDqVRo0ZWlyUiZ8Epgsj+/fvp1q0bzZo1Y+7cueVCSGBg4Flto7oEERG58IqLi/n444+ZM2cOK1aswBjDtm3buPjii60uTUTO4Fzev90rq4jPPvuMXbt2sWvXLpo2bVquTR+YJSJn4uHhQb9+/ejXrx8HDhxgxYoVhIeHA9CnTx/CwsJISEigXbt22Gy2M2xNRJxVpV01k5CQgDHmlJOIyLlo0KABCQkJ2Gw2jh8/TmRkJEuXLqVDhw5ERUXxzDPPUFhYaHWZInIe9FkzIlKtuLu78/TTT5OWlsZHH31EZGQk06dPx9299ADv119/zZEjRyyuUkTOVpXeWfVcaYyIiJyN4uJiPDw8yMvLo3HjxtSqVYshQ4aQkJBA586ddepGpIo57Z1VRUQqg4eHBwD+/v5s2bKFMWPGsHLlSrp06cJll11GcXGxxRWKyJ9REBERlxIREcHkyZNJTU1l1apVjBo1Cg8PD44dO8agQYNYsGABhw/r7scizkJBRERckpubG9deey333nsvUHrTtNzcXIYNG0ZgYCCJiYl89dVXGkAvYjEFERGpEUJCQli9ejW7d+/mvvvuIykpiX/84x/YbDaMMaSlpVldokiNpMGqIlIjORwOcnJyCAwM5Pvvv6djx45069aNhIQEbrrpJmrXrm11iSLVlgarioicgd1uL7vLc6tWrZg3bx7u7u6MGjWKgIAA7r//fosrFKkZFEREpMbz9fVl2LBhfP7556SmpjJhwgSCg4OB0s/MeuSRR0hJSbG4ShHXpFMzIiKn8emnnzJkyBDy8vK4/PLLiY+PZ/DgwdStW9fq0kSclk7NiIhcINdddx2ZmZksXLgQPz8/7rzzTu655x6g9EZqJSUlFlcoUr0piIiInIGPjw+33HILH3/8Menp6Tz++OMAvPPOOzRr1oyHHnqI7du3W1ukSDWlICIicg6aNGlCWFgYAG3btmXAgAG89tprXHrppXTq1IkPPvjA4gpFqhcFERGR89S6dWtmzJhBZmYm7777LgEBAfz6668AbN26lQ8//FC3lxc5AwUREZG/yMvLi0GDBrFixQri4+MBWLJkCX379qVp06bcd999bNmyxeIqRZyTgoiISCV44okn2LhxI0OHDmX+/Pm0adOGl156yeqyRJyOLt8VEalkxcXFfPzxx0RHRxMSEsL//d//sX79ehISErj++uvx9PS0ukSRC0qX74qIOBEPDw/69etHSEgIAC1btmT//v0MHDiQJk2acM8997Br1y6LqxSxhoKIiEgVGzp0KOvXryc5OZlRo0axdOlSsrKyANiyZQuZmZkWVyhSdXRqRkTEYsePH8fNzQ2bzUb37t1Zs2YN1113HQkJCfTr1w9vb2+rSxQ5Jzo1IyJSjbi7u2Oz2QBYtmwZL7/8Mr/++itDhgwhKCiI9evXW1yhSOXRERERESe1c+dOFixYwPjx4/Hx8eGhhx6ibt26DB8+nKZNm1pdnsif0hEREREXEBERwb/+9S98fHwwxnDw4EGeeOIJmjVrRs+ePVmwYAFFRUVWlynylzj1EZG8vDzq1q1Lenq6joiIiFD6n+YHH3zAggUL+Pbbb9mxYweBgYGkpaUREhJSdopHxEr5+fmEhITw66+/4u/vf9p13auopvNSUFAAUHbJm4iIlBcREWF1CSJ/qqCg4IxBxKmPiDgcDvbv30+dOnWcOuWfSH46cuN8tG+cl/aNc9P+cV7VYd8YYygoKKBJkybY7acfBeLUR0Tsdnu1GpDl5+fntH8UNZ32jfPSvnFu2j/Oy9n3zZmOhJygwaoiIiJiGQURERERsYyCyAXg5eXFxIkT8fLysroUOYn2jfPSvnFu2j/Oy9X2jVMPVhURERHXpiMiIiIiYhkFEREREbGMgoiIiIhYRkFERERELKMgch4mT57M5ZdfTq1atahbt+5ZPcYYw+OPP06TJk3w8fHh6quvZtu2bZVbaA116NAhhg8fjr+/P/7+/gwfPpxff/31tI9JSEjAZrOVm2JiYqqmYBf28ssv06JFC7y9vWnfvj1fffXVaddfvXo17du3x9vbm5YtW/LKK69UUaU1z7nsmy+//LJC/7DZbPz4449VWHHNsGbNGvr27UuTJk2w2Wy8//77Z3xMde83CiLn4dixY9x8883ceeedZ/2Yp59+mueee44ZM2awfv16AgMD6dmzZ9nn6ciFM3ToUDZt2sQnn3zCJ598wqZNmxg+fPgZH9erVy8yMzPLpo8++qgKqnVdixcvZuzYsTz88MNs3LiRrl270rt3b9LS0k65/t69e+nTpw9du3Zl48aNTJgwgb///e8sXbq0iit3fee6b07YuXNnuT4SHh5eRRXXHIWFhbRp04YZM2ac1fou0W+MnLc333zT+Pv7n3E9h8NhAgMDzVNPPVW27MiRI8bf39+88sorlVhhzbN9+3YDmG+++aZs2bp16wxgfvzxxz99XHx8vOnfv38VVFhzdOrUydxxxx3lll1yySVm3Lhxp1z/wQcfNJdcckm5ZX/7299MTExMpdVYU53rvklKSjKAOXToUBVUJycA5r333jvtOq7Qb3REpArs3buXrKwsYmNjy5Z5eXnRrVs31q5da2FlrmfdunX4+/vTuXPnsmUxMTH4+/uf8Xf95Zdf0rhxYy6++GJuu+02cnJyKrtcl3Xs2DF++OGHcn/zALGxsX+6H9atW1dh/euuu47vv/+e4uLiSqu1pjmffXNCdHQ0QUFB9OjRg6SkpMosU86SK/QbBZEqkJWVBUBAQEC55QEBAWVtcmFkZWXRuHHjCssbN2582t917969WbBgAf/5z3+YNm0a69ev55prruHo0aOVWa7Lys3NpaSk5Jz+5rOysk65/vHjx8nNza20Wmua89k3QUFBzJo1i6VLl7Js2TIiIiLo0aMHa9asqYqS5TRcod849afvVqXHH3+cSZMmnXad9evX06FDh/N+DpvNVm7eGFNhmZza2e4fqPh7hjP/rocMGVL2fWRkJB06dKB58+asXLmSgQMHnmfVcq5/86da/1TL5a87l30TERFBRERE2XyXLl1IT0/n2Wef5aqrrqrUOuXMqnu/URD53d13380tt9xy2nVCQ0PPa9uBgYFAaXINCgoqW56Tk1Mhycqpne3+2bJlC9nZ2RXafvnll3P6XQcFBdG8eXNSUlLOuVaBhg0b4ubmVuE/7NP9zQcGBp5yfXd3dxo0aFBptdY057NvTiUmJob58+df6PLkHLlCv1EQ+V3Dhg1p2LBhpWy7RYsWBAYGsmrVKqKjo4HS87SrV69m6tSplfKcruZs90+XLl3Iy8vju+++o1OnTgB8++235OXlcfnll5/18x04cID09PRywVHOnqenJ+3bt2fVqlUMGDCgbPmqVavo37//KR/TpUsXVqxYUW7ZZ599RocOHfDw8KjUemuS89k3p7Jx40b1DyfgEv3GypGy1dW+ffvMxo0bzaRJk0zt2rXNxo0bzcaNG01BQUHZOhEREWbZsmVl80899ZTx9/c3y5YtM1u3bjW33nqrCQoKMvn5+Vb8CC6tV69e5rLLLjPr1q0z69atM1FRUeaGG24ot84f909BQYH55z//adauXWv27t1rkpKSTJcuXUxwcLD2z1+waNEi4+HhYWbPnm22b99uxo4da3x9fU1qaqoxxphx48aZ4cOHl62/Z88eU6tWLXPvvfea7du3m9mzZxsPDw/z7rvvWvUjuKxz3Tf//ve/zXvvvWd++uknk5ycbMaNG2cAs3TpUqt+BJdVUFBQ9p4CmOeee85s3LjR7Nu3zxjjmv1GQeQ8xMfHG6DClJSUVLYOYN58882yeYfDYSZOnGgCAwONl5eXueqqq8zWrVurvvga4MCBAyYuLs7UqVPH1KlTx8TFxVW47PCP++fw4cMmNjbWNGrUyHh4eJhmzZqZ+Ph4k5aWVvXFu5iXXnrJNG/e3Hh6epp27dqZ1atXl7XFx8ebbt26lVv/yy+/NNHR0cbT09OEhoaamTNnVnHFNce57JupU6easLAw4+3tberVq2euvPJKs3LlSguqdn0nLpU+eYqPjzfGuGa/sRnz+6gWERERkSqmy3dFRETEMgoiIiIiYhkFEREREbGMgoiIiIhYRkFERERELKMgIiIiIpZREBERERHLKIiIiIiIZRRERERExDIKIiIiImIZBRERERGxjIKIiIiIWOb/AVqYDa5btuMtAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/fUlEQVR4nO3deVyU9d7/8fcAgqgwpsgmCNqvo6kthuWWuYZ7iyfTW1OsTo+oPGnmnXq6K7W7bDl17LSZ57awbPHk0klckk6ilrRoam6ZFSkhhCug6YDw/f0xMacRRFBmhgtez8djHjLf+V4zn68X18yb6/pe19iMMUYAAAAW4efrAgAAAKqD8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACwlwNcF1LTS0lIdOHBAISEhstlsvi4HAABUgTFGhYWFio6Olp9f5ftW6lx4OXDggGJjY31dBgAAOA9ZWVmKiYmptE+dCy8hISGSnIMPDQ31cTUAAKAqCgoKFBsb6/ocr0ydCy9lh4pCQ0MJLwAAWExVpnwwYRcAAFgK4QUAAFgK4QUAAFgK4QUAAFgK4QUAAFgK4QUAAFgK4QUAAFgK4QUAAFhKnbtIHQAA8IySEmnDBiknR4qKknr2lPz9vV+HR/e8zJ49W1dffbVCQkIUHh6um266SXv27DnncuvWrVNCQoIaNmyoNm3aaO7cuZ4sEwAAnMPSpVJ8vNSnjzR6tPPf+Hhnu7d5NLysW7dO9913nz7//HOlpaXp9OnTSkxM1IkTJ866TGZmpgYPHqyePXtqy5Yt+stf/qL7779fS5Ys8WSpAADgLJYulW65Rfr5Z/f27Gxnu7cDjM0YY7z1YgcPHlR4eLjWrVun6667rsI+U6dO1Ycffqjdu3e72pKTk7Vt2zZlZGSc8zUKCgpkt9uVn5/PdxsBAHCBSkqce1jODC5lbDYpJkbKzLywQ0jV+fz26oTd/Px8SVKzZs3O2icjI0OJiYlubQMGDNCmTZtUXFxcrr/D4VBBQYHbDQAA1IwNG84eXCTJGCkry9nPW7wWXowxmjx5sq699lp17NjxrP1yc3MVERHh1hYREaHTp0/r0KFD5frPnj1bdrvddYuNja3x2gEAqK9ycmq2X03wWniZMGGCvvnmG7377rvn7Hvm12GXHdmq6Guyp0+frvz8fNctKyurZgoGAACKiqrZfjXBK6dK//nPf9aHH36o9evXKyYmptK+kZGRys3NdWvLy8tTQECAmjdvXq5/UFCQgoKCarReAADg1LOnc05LdrbzENGZyua89OzpvZo8uufFGKMJEyZo6dKl+uSTT9S6detzLtOtWzelpaW5ta1Zs0adO3dWgwYNPFUqAACogL+/9MILzp/PPABSdn/OHO9e78Wj4eW+++7TwoUL9c477ygkJES5ubnKzc3VyZMnXX2mT5+ucePGue4nJydr3759mjx5snbv3q3XX39d8+fP15QpUzxZKgAAOIvhw6XFi6WWLd3bY2Kc7cOHe7cej54qXdEcFUl64403NH78eEnS+PHj9dNPPyk9Pd31+Lp16/TAAw9o586dio6O1tSpU5WcnFyl1+RUaQAAPMOTV9itzue3V6/z4g2EFwAArKfWXucFAADgQhFeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApRBeAACApXg0vKxfv17Dhg1TdHS0bDabPvjgg0r7p6eny2azlbt9++23niwTAABYSIAnn/zEiRO64oordPvtt+uPf/xjlZfbs2ePQkNDXfdbtGjhifIAAIAFeTS8DBo0SIMGDar2cuHh4WratGnNFwQAACyvVs556dSpk6KiotSvXz+tXbu20r4Oh0MFBQVuNwAAUHfVqvASFRWlefPmacmSJVq6dKnatm2rfv36af369WddZvbs2bLb7a5bbGysFysGAADeZjPGGK+8kM2mZcuW6aabbqrWcsOGDZPNZtOHH35Y4eMOh0MOh8N1v6CgQLGxscrPz3ebNwMAAGqvgoIC2e32Kn1+16o9LxXp2rWr9u7de9bHg4KCFBoa6nYDAAB1V60PL1u2bFFUVJSvywAAALWER882On78uL7//nvX/czMTG3dulXNmjVTq1atNH36dGVnZ+vNN9+UJM2ZM0fx8fHq0KGDioqKtHDhQi1ZskRLlizxZJkAAMBCPBpeNm3apD59+rjuT548WZKUlJSklJQU5eTkaP/+/a7Hi4qKNGXKFGVnZys4OFgdOnTQihUrNHjwYE+WCQAALMRrE3a9pToTfgAAQO1QpybsAgAA/B7hBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWEqArwsAgKooKZE2bJBycqSoKKlnT8nf39dVAfVLbdkOPbrnZf369Ro2bJiio6Nls9n0wQcfnHOZdevWKSEhQQ0bNlSbNm00d+5cT5YIwAKWLpXi46U+faTRo53/xsc72wF4R23aDj0aXk6cOKErrrhCL730UpX6Z2ZmavDgwerZs6e2bNmiv/zlL7r//vu1ZMkST5YJoBZbulS65Rbp55/d27Ozne0EGMDzatt2aDPGGK+8kM2mZcuW6aabbjprn6lTp+rDDz/U7t27XW3Jycnatm2bMjIyqvQ6BQUFstvtys/PV2ho6IWWDcCHSkqcf9md+YZZxmaTYmKkzEwOIQGe4q3tsDqf37VqzktGRoYSExPd2gYMGKD58+eruLhYDRo0KLeMw+GQw+Fw3S8oKPB4nQC8Y8OGs79hSpIxUlaWs1/v3l4rC6gT5s2bp4kTJ6pFixZq1aqV4uLi1LJlS8XExLj9u2dPpH7++exxwRfbYa0KL7m5uYqIiHBri4iI0OnTp3Xo0CFFRUWVW2b27NmaOXOmt0oE4EU5OTXbD8B/NG/eXKdOnVJWVpaysrKUkZEhf39/lZSUqLS01NXPZvOT1ExS499uX0hqUu75vLkd1rpTpW02m9v9sqNaZ7aXmT59uvLz8123rKwsj9cIwDvK/73ys6Tuko6cox+Acxk6dKjb4ZnS0lIVFxe7BRdJMqZU0iFJ+yT9IKni2Sbe3A5rVXiJjIxUbm6uW1teXp4CAgLUvHnzCpcJCgpSaGio2w1A3dCzp/NY+n/+dpkqKUPSDEnO9thYZz8A1RMUFKTbbrtNAQFVOQhjk3SRpJ2SQtwf8cF2WKvCS7du3ZSWlubWtmbNGnXu3LnC+S4A6jZ/f+mFF8ru5Ula/NvPr0tyzm+bM4fJusD5Gj9+vE6fPl1pHz8/PzVs2ERSumy2i90eK/vDwtvboUfDy/Hjx7V161Zt3bpVkvNU6K1bt2r//v2SnId8xo0b5+qfnJysffv2afLkydq9e7def/11zZ8/X1OmTPFkmQBqseHDpcWLpSZNnpVU9ib7q+z2F7V4sfNxANVXWlqq/Px8hYSEnLWPn5+fAgMD9cknH2nJksvVsqX74zEx8sl26NFTpdPT09WnT59y7UlJSUpJSdH48eP1008/KT093fXYunXr9MADD2jnzp2Kjo7W1KlTlZycXOXX5FRpoO7Jy8tTq1ZxcjhOudpCQkL1889ZbOdANe3du1cLFizQm2++qaysLIWFhenw4cM6Mw7YbDb5+/tr5cqVuv766yV59gq71fn89tp1XryF8ALUPf/93/+tv/3tbyopKXG1+fn5adasWXr44Yd9WBlgDceOHdM///lPpaSkKCMjQ3a7XaNGjVJSUpJatWql2NjYCsPL4sWLNdxLu1UIL4QXoM7Iy8tTXFycTp06Ve6x0NBQZWWx9wWoSElJidLS0pSSkqIPPvhAxcXFGjBggJKSknTDDTcoODjY1TcxMVGffPKJ2x8I8+fP1x133OG1eqvz+V2rJuwCwJmeffZZFRcXV/jY8ePH9eKLL3q5IqB227lzpx566CHFxsZq0KBB2rFjhx5//HFlZWVp5cqVGjlypFtwkaQ77rjDLbg899xzXg0u1cWeFwC1VmV7Xcqw9wWQDh8+rPfee08pKSnatGmTmjVrptGjRyspKUkJCQlnvVZamZMnTyo8PFzHjx/Xww8/rP/93//1UuX/YdmvBwCA36tsr0uZsr0vzH1BfVNcXKzVq1crJSVFy5cvlzFGgwcP1pIlSzRkyBAFBQVV+bmCg4P15JNP6ujRo3rkkUc8WHXNYM8LgFqpKntdyrD3BfXJtm3btGDBAr399tvKy8vTlVdeqaSkJI0ePVrh4eG+Lu+8MecFgOU999xzVQoukvNN76WXXjqv18nNzdX9999f7pLoQG2Sl5enOXPm6Morr9SVV16phQsXasyYMdq6dau2bNmiSZMmWTq4VBeHjQDUSqGhobr4YvereZ46dUrZ2dmKjY1VYGCg22NNmpT/oriKzJs3T++8845b2549e9S2bVvdd999F1Y0UIMcDodWrFihlJQUrVq1SjabTTfccIMef/xxDRw4sF5feZ7DRgAsIyMjQ927d9eOHTvUoUOHGnnOI0eOaNy4cVq4cKGaNm1aI88JnC9jjDZv3qyUlBS9++67OnLkiDp37qzx48dr1KhRZ/2ev7qACbsAUEVNmjTRkiVLqjW5EahpOTk5WrhwoVJSUrRr1y5FRUXpT3/6k5KSktS+fXtfl1frEF4A1GtnHn4CvOXUqVP617/+pZSUFK1Zs0YNGjTQzTffrOeee079+/ev4rc9109M2AVQ78yaNUuXXXaZGjdurIiICN1zzz3nPCUbqAnGGGVkZGjcuHGuS/QXFBTo1VdfVW5urt59910NHDiQ4HIO/O8AqFeMMSopKdFrr72mli1bateuXRo3bpwuv/xy3XPPPb4uD3XUU089pddee02//PKLTp48qaCgIDVo0EDTpk3TzJkzfV2e5RBeANQrNpvN7cMiLi5O119/vb799ltlZWVp7NixysvLU0BAgB555BGNGDHCh9XCyk6cOKFly5YpJSVFn3zyiYKDg3XLLbcoKSlJV1xxhcaPH68HHnjA12VaEuEFQL2yb98+Pfvss0pPT1d2draKi4t16tQpzZ49WwEBAa5raeTl5emqq67S4MGD1bhxY1+XDYsoLS3Vp59+qpSUFL3//vs6fvy4evXqpfnz5+uWW25RSEiIJKmoqIiJ4heA8AKg3jh06JCuueYa9enTR88//7xatmyp0tJSde7cWVdeeaWioqIUFRUlSQoPD1ezZs105MgRwgvOKTMzU2+++aYWLFigzMxMtW7dWlOmTNHYsWPVpk2bcv2ZKH5hCC8A6o2VK1fq9OnTevfdd11fVPfyyy+rqKhIV155pVvfTZs2qbS0VLGxsT6oFFZQWFioxYsXKyUlRevXr1eTJk106623KiUlRddee638/NzPiZk1a5bef/99/fjjj2rSpImGDx+uv//97/X6YnPni/ACoN5o1qyZCgoK9OGHH6p9+/Zavny5Zs+erZYtW6pFixaufocPH9a4ceP0f//3fz6sFrVRaWmp1q5dq5SUFC1dulQnT55U37599dZbb+nmm28+6146JorXLMILgHpjyJAhuvPOOzV27FgFBwfrtttu06233qp9+/a5+jgcDt18882aPn26unfv7sNqUZt89913WrBggd566y1lZWXpkksu0cMPP6zbbrtNrVq1OufylU0Ul6Sbb75Z6enp6tevnxYvXuyxcdQVhBcA9YbNZtPcuXM1d+7cCh83xmj8+PHq27evxo4d6+XqUNscO3ZMixYt0oIFC5SRkSG73a6RI0dq/Pjx6tq1q+vQY1VUNlFcku6//37dcccdWrBggaeGU6dwkToA+M1nn32mRYsW6YMPPnB9e+/27dt9XRa8qKSkRKtXr9aoUaMUGRmpe++9V3a7Xe+9955ycnL02muvqVu3btUKLmUTxQ8dOqTnn39en376qTIyMuTv7++aa9WnTx/XmUg4N/a8AMBvrr32WpWWlvq6DPjAzp07tWDBAi1cuFA5OTlq3769Hn/8cY0ZM0bR0dEX9NzVmSiOqiG8AADqpcOHD+u9995TSkqKNm3apGbNmmn06NFKSkpSQkJCtfauVKaqE8VRdYQXAEC9UVxcrNWrVyslJUXLly9XaWmpBg8erCVLlmjIkCEeuWhcVSaKo3oILwCAOm/btm1KSUnRO++8o7y8PF1xxRV65plnNHr0aIWHh3v0tc81URzVR3gBANRJeXl5evvtt7VgwQJt27ZNLVq00JgxY5SUlFTr5poMGDBAX3/9tU6cOKGYmBgtW7ZMV199ta/LqrUILwCAOsPhcGjFihVKSUnRqlWrZLPZNGzYMD3++OMaOHBgrb2a7UcffeTrEiyF8AIAsDRjjDZv3qyUlBS9++67OnLkiDp37qw5c+Zo1KhRat68ua9LRA0jvAAALCknJ0cLFy5USkqKdu3apaioKN15551KSkpShw4dfF0ePIjwAgCwjFOnTulf//qXUlJStGbNGjVo0EA33XSTnnvuOfXv318BAXys1QesZQBArWaM0eeff66UlBQtWrRI+fn56tatm1555RWNHDlSTZs29XWJ8DLCC+qP0hLp4AbpZI4UHCW16Cn5+fu6KqD+qOY2uH//fr311ltasGCB9u7dq5iYGN13330aN26c2rZt68XCUdt45buNXnnlFbVu3VoNGzZUQkKCNmzYcNa+6enpstls5W5l37wJnJespdKH8dK/+0gbRzv//TDe2Q7A86q4DZ44cUILFy5U//79FR8fryeeeEJdunRRWlqafvrpJz3xxBMEF3h+z8uiRYs0adIkvfLKK+rRo4dee+01DRo0SLt27ar0a8T37Nmj0NBQ130uoYzzlrVU2nCLJOPe/mu2s73nYil2uE9KA+qFc2yDpT3+qU/3hSslJUXvv/++jh8/ruuuu07z58/XLbfcwhcWohyPh5fnn39ed955p/70pz9JkubMmaOPPvpIr776quurwCsSHh7OcUxcuNISafNElXvTlH5rs0mbJ0ktb+QQEuAJlWyDmXlGb26QFkwepcxfStS6dWs9+OCDGjdunNq0aeP9WmEZHg0vRUVF2rx5s6ZNm+bWnpiYqI0bN1a6bKdOnXTq1Cm1b99e//M//6M+ffpU2M/hcMjhcLjuFxQUXHjhqDsObpB+/dmtqei0n06b333hWv4Baf/HUnhPLxeH6jp16pQk6eTJk/r11199XA2qJG+DlJ8jyfnHQeFJow83GS38zGj9t1KThtKILiVKmTtH197wZ/n5eWU2AyzOo+Hl0KFDKikpUUREhFt7RESEcnNzK1wmKipK8+bNU0JCghwOh9566y3169dP6enpuu6668r1nz17tmbOnOmR+lEHnMxxu1t02k+bDjTTiaIzfvUzvpDCeNOs7Xbt2iVJ2rJli44ePerjalAlh75QaWaYtv1UpDVbTurT3Q4VnTbqfan0ZrI0/GqpcUNJV4ZLBBdUkVfONjrza8WNMWf9qvG2bdu6Tcbq1q2bsrKy9Ne//rXC8DJ9+nRNnjzZdb+goECxsbE1VDksLzjK7e5pY9OJogAF+pcq0L/0Pw/YI6UmTbxcHKorODhYktSoUSM1YX3VellZWVr5r8+1evVB/ZJfqtjmfrqtVyP1uayRbr7iiBo1KPlP5zO2VaAyHg0vYWFh8vf3L7eXJS8vr9zemMp07dpVCxcurPCxoKAgj3yFOeqIFj2lRjHOiYG/O+Ye6F+qoIDfwktQCym8E3NeLCAwMND1L9t97VRYWKi0tDQtX75c33zzjZo0aaLEy0M0tFOpLo8LUFGJv44X/X5bszm30RYctkXVeTS8BAYGKiEhQWlpabr55ptd7Wlpabrxxhur/DxbtmxRVBSpHOfBz19KeOG3Mx0q3tunS+4luAAXoKSkRF988YWWL1+u9PR0nT59Wl27dtWTTz6pXr16qWHBJmlHRYf3f9smE+awDaJaPH7YaPLkyRo7dqw6d+6sbt26ad68edq/f7+Sk5MlOQ/7ZGdn680335TkPBspPj5eHTp0UFFRkRYuXKglS5ZoyZIlni4VdVXscOfp0Jsn/jZx8DdBLZzBJfxa39UGWNgPP/yg1NRUrVy5UocOHVKbNm2UnJyswYMHu1/eouG1UsfHpL0vS6eP/Ke9UYwzuHCpAlSTx8PLyJEjdfjwYc2aNUs5OTnq2LGjVq5cqbi4OEnOL9bav3+/q39RUZGmTJmi7OxsBQcHq0OHDlqxYoUGDx7s6VJRl8UOd54Ovf9j5+RceySHioDzcOzYMa1Zs0bLly/Xrl27ZLfbNWDAAA0dOlTt27c/63xGhV8rhXWT8rZI+blSty5Sq/5sgzgvNmNMRRfAsKyCggLZ7Xbl5+e7XeQOkKRff/1Vn332mZo0acKcCQv65ptvdPvtt+uf//ynLr74Yl+XU28UFxcrIyNDy5cv1/r162WMUY8ePTR06FD17NnTNRepKhwOh44fP64ePXqoUaNGHqwaVlOdz2++2wgAUKHvvvtOqampWrVqlY4cOaI//OEPuv/++zVw4EA1b97c1+WhHiO8AABcjhw5otWrV2v58uX67rvvdNFFF2nQoEEaOnQo3ymEWoPwAgD1XHFxsTZs2KDU1FR9+umnstlsuu6665ScnKzu3burQYMGvi4RcEN4AYB6yBij3bt3KzU1VatXr1Z+fr7at2+vBx98UAMGDOC75VCrEV4AoB45ePCgVq1apdTUVP3www8KCwvTjTfeqKFDhzIJGpZBeAGAOs7hcGjdunVKTU1VRkaGAgIC1Lt3b02cOFFdunRRQAAfBbAWfmMBoA4yxmj79u1KTU3VmjVrVFhYqMsvv1zTpk3T9ddfz6UkYGmEFwCoQ3Jzc7Vy5UqlpqZq3759ioiI0IgRIzRkyBDFx8f7ujygRhBeAMDiTp48qbVr1yo1NVVffvmlAgMD1a9fP02dOlWdO3eWvz9XsUXdQngBAAsyxmjLli1KTU3Vxx9/rBMnTuiqq67SI488on79+qlJkya+LhHwGMILAFhIdna2VqxYodTUVGVnZ6tly5YaM2aMhgwZopiYGF+XB3gF4QUAarkTJ07o3//+t1JTU7V582Y1atRI/fv312OPPaZOnTrJz8/P1yUCXkV4AYBaqLS0VJs2bVJqaqr+/e9/y+Fw6Oqrr9asWbPUt29fBQcH+7pEwGcILwBQi+zfv1+pqalasWKFcnNzFRsbqzvuuEODBw9WVFSUr8sDagXCCwD4WGFhodLS0pSamqpt27apcePGGjBggIYOHarLL79cNpvN1yUCtQrhBQB8oKSkRF988YVSU1OVnp6u4uJide3aVU8++aR69eqlhg0b+rpEoNYivACAF/34449KTU3VypUrdfDgQbVp00Z33323Bg8erBYtWvi6PMASCC8A4GH5+fn66KOPlJqaqp07dyo0NFQDBw7U0KFD1b59ew4LAdVEeAEADyguLlZGRoZSU1O1fv16lZaWqkePHnrmmWfUs2dPBQYG+rpEwLIILwBQg/bu3avly5dr1apVOnLkiC655BL9+c9/1sCBA9W8eXNflwfUCYQXALhAR48e1apVq5Samqo9e/booosu0qBBgzR06FC1bdvW1+UBdQ7hBQDOQ3FxsTZs2KDU1FR9+umnstlsuu6663T33Xere/fuatCgga9LBOoswgsAVJExRt9++62WL1+u1atXKz8/X5deeqkefPBBJSYm6qKLLvJ1iUC9QHgBgHM4ePCg67DQDz/8oObNm+vGG2/U0KFDdfHFF/u6PKDeIbwAQAUcDofWrVun1NRUZWRkKCAgQL1799bEiRPVpUsXBQTw9gn4ClsfAPzGGKMdO3Zo+fLlWrNmjQoLC3X55Zdr2rRpuv766xUaGurrEgGI8AIA+uWXX7RixQqlpqZq3759ioiI0IgRIzRkyBDFx8f7ujwAZyC8AKiXTp48qbVr1yo1NVVffvmlAgMD1a9fP02dOlWdO3eWv7+/r0sEcBaEFwD1hjFGW7du1fLly/Xxxx/rxIkT6tSpkx555BH169dPTZo08XWJAKqA8AKgzjtw4IBSU1O1YsUK/fzzz4qOjtaYMWM0ZMgQxcTE+Lo8ANVEeKmikhJpwwYpJ0eKipJ69pTYqwx4T0mJtHev8+c9e6T4+Mq3wV9//VUff/yxUlNTtXnzZgUHB6t///569NFH1alTJ/n5+XmlbqAuqS2fhV7Zel955RW1bt1aDRs2VEJCgjZs2FBp/3Xr1ikhIUENGzZUmzZtNHfuXG+UeVZLlzrfKPv0kUaPdv4bH+9sB+B5n34q3Xab9OKLzvvPPuu8/+mn7v1KS0v15Zdf6tFHH1ViYqJmzZolPz8/zZo1S2lpaZoxY4YSEhIILsB5qE2fhR7f87Jo0SJNmjRJr7zyinr06KHXXntNgwYN0q5du9SqVaty/TMzMzV48GDdddddWrhwoT777DPde++9atGihf74xz96utxyli6VbrlFMsa9PTvb2b54sTR8uNfLAuqNTz+VZs4s337okLP9scekuLgsLV++XCtWrFBubq5iY2N1++23a/DgwYqKivJ+0UAdU9s+C23GnFlKzerSpYuuuuoqvfrqq662Sy+9VDfddJNmz55drv/UqVP14Ycfavfu3a625ORkbdu2TRkZGed8vYKCAtntduXn51/wNRlKSpyp8uefK37cZpNiYqTMTA4hWcWvv/6qzz77TE2aNFFQUJCvy8E5lJQ497AcOuS8X1z8jQoLb1do6D/l5xeuoqI0GZOqkye3qXHjxhowYICGDh2qyy+/XDabzbfFo0IOh0PHjx9Xjx491KhRI1+Xgyrw1mdhdT6/PbrnpaioSJs3b9a0adPc2hMTE7Vx48YKl8nIyFBiYqJb24ABAzR//nwVFxeX+7Izh8Mhh8Phul9QUFBD1TuP651tZUnOBJqV5ezXu3eNvSyA3+zY8Z/g8nu//vqCTp/eJKlYDRp00d13P6lx43qpYcOGXq8RqOtq42ehRw/8Hjp0SCUlJYqIiHBrj4iIUG5uboXL5ObmVtj/9OnTOlTBu9js2bNlt9tdt9jY2BqrPyenZvsBqJ7Dh89sKZUklZTsV3Dw3WradKVCQl7S//t/AwgugIfUxs9Cr8xaO3P3rTGm0l26FfWvqF2Spk+frvz8fNctKyurBip2quqhcg6pA57RvLn7/QYNrpTd/qGaNl2m4OAk+fm1qLAfgJpTGz8LPXrYKCwsTP7+/uX2suTl5ZXbu1ImMjKywv4BAQFqXsE7VFBQkMfmLvTs6TyOl51dfpKS9J/jfD17euTlgXqvY0cpLMz90JG/f0u3Pi1aOPsB8Iza+Fno0T0vgYGBSkhIUFpamlt7WlqaunfvXuEy3bp1K9d/zZo16ty5c7n5Lp7m7y+98ILz5zN3+pTdnzOHybqAp/j7S/fdV3mfe+9lGwQ8qTZ+Fnr8sNHkyZP1f//3f3r99de1e/duPfDAA9q/f7+Sk5MlOQ/7jBs3ztU/OTlZ+/bt0+TJk7V79269/vrrmj9/vqZMmeLpUis0fLjzFLCW7n/sKSaG06QBb7j2Wufp0GFh7u0tWjjbr73WN3UB9Ult+yz0+HVeRo4cqcOHD2vWrFnKyclRx44dtXLlSsXFxUmScnJytH//flf/1q1ba+XKlXrggQf08ssvKzo6Wn//+999co2XMsOHSzfeWDuuKgjUR9deK3Xr5jz76PBh5xyXjh3ZBgFvqk2fhR6/zou31eR1XlD3cJ0XwLe4zgvOpjqf31wjGwAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBQAAWArhBfVHaYmUt0E6lCEd3e68D8B7Skuc296hDOe2yDaI8+TR8HL06FGNHTtWdrtddrtdY8eO1bFjxypdZvz48bLZbG63rl27erJM1AdZS6UP46V1Q6S9r0rfPCxl3CblferryoD6Ie9T5zb3zcPObXDdEOc2mbXU15XBgjwaXkaPHq2tW7dq9erVWr16tbZu3aqxY8eec7mBAwcqJyfHdVu5cqUny0Rdl7VU2nCL9OvP7u2OQ9KOmQQYwNPyPnVua45D7u2/Zju3TQIMqinAU0+8e/durV69Wp9//rm6dOkiSfrHP/6hbt26ac+ePWrbtu1Zlw0KClJkZKSnSkN9UloibZ4oyZy9z95XpLBukp+/18oC6o3SEmnvy2d50EiySZsnSS1vZBtElXksvGRkZMhut7uCiyR17dpVdrtdGzdurDS8pKenKzw8XE2bNlWvXr30xBNPKDw8vMK+DodDDofDdb+goKDmBgHrO7ih/B4XSUUlv9vpePqwlLdFuugyLxYG1BNHt0snjqhsR7/btidJMtKvWc5tNaK3t6uDRXksvOTm5lYYOMLDw5Wbm3vW5QYNGqQRI0YoLi5OmZmZeuSRR9S3b19t3rxZQUFB5frPnj1bM2fOrNHaUYeczHG7G2Azahx4WieKAtzfRPNzpQatvVwcUA/k50pF7h81jQNPK8B2xt7QM7ZVoDLVDi8zZsw4Z1j46quvJEk2m63cY8aYCtvLjBw50vVzx44d1blzZ8XFxWnFihUaPnx4uf7Tp0/X5MmTXfcLCgoUGxt7znGgngiOcrsbGFCqztFHdNqc8TvYrYsU3sOLhQH1RF6pVDTLrSnAZhQYUOre74xtFahMtcPLhAkTNGrUqEr7xMfH65tvvtEvv/xS7rGDBw8qIiKiyq8XFRWluLg47d27t8LHg4KCKtwjA0iSWvSUGsU4Jwb+Nu8lMKBUga4ONufjrfpzvB3whFb9JXuU2zbo7rdtsEVPb1cGC6t2eAkLC1NYWNg5+3Xr1k35+fn68ssvdc0110iSvvjiC+Xn56t79+5Vfr3Dhw8rKytLUVGkcpwHP38p4QXnGQ2yyf3N87e9LwlzCC6Ap7ANwgM8dqr0pZdeqoEDB+quu+7S559/rs8//1x33XWXhg4d6jZZt127dlq2bJkk6fjx45oyZYoyMjL0008/KT09XcOGDVNYWJhuvvlmT5WKui52uNRzsdSopXt7oxhne2z5w5EAahDbIGqYxybsStLbb7+t+++/X4mJiZKkG264QS+99JJbnz179ig/P1+S5O/vr+3bt+vNN9/UsWPHFBUVpT59+mjRokUKCQnxZKmo62KHO0/FPLjBOTEwOMq5m5q/9gDvYBtEDbIZYyq5AIb1FBQUyG63Kz8/X6Ghob4uBwAAVEF1Pr/5biMAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGAphBcAAGApAb4uwCpKSqQNG6ScHCkqSurZU/L393VVQP3BNgigjEf3vDzxxBPq3r27GjVqpKZNm1ZpGWOMZsyYoejoaAUHB6t3797auXOnJ8s8p6VLpfh4qU8fafRo57/x8c52AJ7HNgjg9zwaXoqKijRixAjdc889VV7mmWee0fPPP6+XXnpJX331lSIjI3X99dersLDQg5We3dKl0i23SD//7N6ene1s580T8Cy2QQBnshljjKdfJCUlRZMmTdKxY8cq7WeMUXR0tCZNmqSpU6dKkhwOhyIiIvT000/r7rvvPudrFRQUyG63Kz8/X6GhoRdUd0mJ86+7M980y9hsUkyMlJnJ7mvAE9gGgfqjOp/ftWrCbmZmpnJzc5WYmOhqCwoKUq9evbRx48YKl3E4HCooKHC71ZQNG87+pilJxkhZWc5+AGoe2yCAitSq8JKbmytJioiIcGuPiIhwPXam2bNny263u26xsbE1Vk9OTs32A1A9bIMAKlLt8DJjxgzZbLZKb5s2bbqgomw2m9t9Y0y5tjLTp09Xfn6+65aVlXVBr/17UVE12w9A9bANAqhItU+VnjBhgkaNGlVpn/j4+PMqJjIyUpJzD0zU796N8vLyyu2NKRMUFKSgoKDzer1z6dnTeTw9O9u5e/pMZcfbe/b0yMsD9R7bIICKVDu8hIWFKSwszBO1qHXr1oqMjFRaWpo6deokyXnG0rp16/T000975DUr4+8vvfCC84wGm839zbNsR9CcOUwUBDyFbRBARTw652X//v3aunWr9u/fr5KSEm3dulVbt27V8ePHXX3atWunZcuWSXIeLpo0aZKefPJJLVu2TDt27ND48ePVqFEjjR492pOlntXw4dLixVLLlu7tMTHO9uHDfVIWUG+wDQI4k0evsPvoo49qwYIFrvtle1PWrl2r3r17S5L27Nmj/Px8V5+HHnpIJ0+e1L333qujR4+qS5cuWrNmjUJCQjxZaqWGD5duvJGrewK+wjYI4Pe8cp0Xb6rJ67wAAADvsOx1XgAAAM6F8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACyF8AIAACzFo+HliSeeUPfu3dWoUSM1bdq0SsuMHz9eNpvN7da1a1dPlgkAACzEo+GlqKhII0aM0D333FOt5QYOHKicnBzXbeXKlR6qEAAAWE2AJ5985syZkqSUlJRqLRcUFKTIyEgPVAQAAKyuVs55SU9PV3h4uP7whz/orrvuUl5e3ln7OhwOFRQUuN0AAEDdVevCy6BBg/T222/rk08+0XPPPaevvvpKffv2lcPhqLD/7NmzZbfbXbfY2FgvVwwAALyp2uFlxowZ5SbUnnnbtGnTeRc0cuRIDRkyRB07dtSwYcO0atUqfffdd1qxYkWF/adPn678/HzXLSsr67xfGwAA1H7VnvMyYcIEjRo1qtI+8fHx51tPOVFRUYqLi9PevXsrfDwoKEhBQUE19noAAKB2q3Z4CQsLU1hYmCdqqdDhw4eVlZWlqKgor70mAACovTw652X//v3aunWr9u/fr5KSEm3dulVbt27V8ePHXX3atWunZcuWSZKOHz+uKVOmKCMjQz/99JPS09M1bNgwhYWF6eabb/ZkqQAAwCI8eqr0o48+qgULFrjud+rUSZK0du1a9e7dW5K0Z88e5efnS5L8/f21fft2vfnmmzp27JiioqLUp08fLVq0SCEhIZ4sFQAAWITNGGN8XURNKigokN1uV35+vkJDQ31dDgAAqILqfH7XulOlAQAAKkN4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlhLg6wKsoqRE2rBBysmRoqKknj0lf39fVwUAQP3jsT0vP/30k+688061bt1awcHBuvjii/XYY4+pqKio0uWMMZoxY4aio6MVHBys3r17a+fOnZ4qs0qWLpXi46U+faTRo53/xsc72wEAgHd5LLx8++23Ki0t1WuvvaadO3fqb3/7m+bOnau//OUvlS73zDPP6Pnnn9dLL72kr776SpGRkbr++utVWFjoqVIrtXSpdMst0s8/u7dnZzvbCTAAAHiXzRhjvPVizz77rF599VX9+OOPFT5ujFF0dLQmTZqkqVOnSpIcDociIiL09NNP6+677z7naxQUFMhutys/P1+hoaEXVG9JiXMPy5nBpYzNJsXESJmZHEICAOBCVOfz26sTdvPz89WsWbOzPp6Zmanc3FwlJia62oKCgtSrVy9t3LixwmUcDocKCgrcbjVlw4azBxdJMkbKynL2AwAA3uG18PLDDz/oxRdfVHJy8ln75ObmSpIiIiLc2iMiIlyPnWn27Nmy2+2uW2xsbI3VnJNTs/0AAMCFq3Z4mTFjhmw2W6W3TZs2uS1z4MABDRw4UCNGjNCf/vSnc76GzWZzu2+MKddWZvr06crPz3fdsrKyqjuks4qKqtl+AADgwlX7VOkJEyZo1KhRlfaJj493/XzgwAH16dNH3bp107x58ypdLjIyUpJzD0zU7xJBXl5eub0xZYKCghQUFFTF6qunZ0/nnJbsbOchojOVzXnp2dMjLw8AACpQ7fASFhamsLCwKvXNzs5Wnz59lJCQoDfeeEN+fpXv6GndurUiIyOVlpamTp06SZKKioq0bt06Pf3009Ut9YL5+0svvOA8q8hmcw8wZTuC5sxhsi4AAN7ksTkvBw4cUO/evRUbG6u//vWvOnjwoHJzc8vNXWnXrp2WLVsmyXm4aNKkSXryySe1bNky7dixQ+PHj1ejRo00evRoT5VaqeHDpcWLpZYt3dtjYpztw4f7pCwAAOotj11hd82aNfr+++/1/fffKyYmxu2x35+dvWfPHuXn57vuP/TQQzp58qTuvfdeHT16VF26dNGaNWsUEhLiqVLPafhw6cYbucIuAAC1gVev8+INNXmdFwAA4B219jovAAAAF4rwAgAALIXwAgAALIXwAgAALIXwAgAALIXwAgAALIXwAgAALIXwAgAALIXwAgAALMVjXw/gK2UXDC4oKPBxJQAAoKrKPrercuH/OhdeCgsLJUmxsbE+rgQAAFRXYWGh7HZ7pX3q3HcblZaW6sCBAwoJCZHNZqvR5y4oKFBsbKyysrLq5Pcm1fXxSXV/jIzP+ur6GOv6+KS6P0ZPjc8Yo8LCQkVHR8vPr/JZLXVuz4ufn1+5b7GuaaGhoXXyF7JMXR+fVPfHyPisr66Psa6PT6r7Y/TE+M61x6UME3YBAIClEF4AAIClEF6qISgoSI899piCgoJ8XYpH1PXxSXV/jIzP+ur6GOv6+KS6P8baML46N2EXAADUbex5AQAAlkJ4AQAAlkJ4AQAAlkJ4AQAAlkJ4+Z0nnnhC3bt3V6NGjdS0adMqLWOM0YwZMxQdHa3g4GD17t1bO3fudOvjcDj05z//WWFhYWrcuLFuuOEG/fzzzx4YwbkdPXpUY8eOld1ul91u19ixY3Xs2LFKl7HZbBXenn32WVef3r17l3t81KhRHh5NeeczvvHjx5ervWvXrm59ass6rO74iouLNXXqVF122WVq3LixoqOjNW7cOB04cMCtny/X3yuvvKLWrVurYcOGSkhI0IYNGyrtv27dOiUkJKhhw4Zq06aN5s6dW67PkiVL1L59ewUFBal9+/ZatmyZp8o/p+qMb+nSpbr++uvVokULhYaGqlu3bvroo4/c+qSkpFS4PZ46dcrTQzmr6owxPT29wvq//fZbt35WXYcVvZ/YbDZ16NDB1ac2rcP169dr2LBhio6Ols1m0wcffHDOZWrFNmjg8uijj5rnn3/eTJ482djt9iot89RTT5mQkBCzZMkSs337djNy5EgTFRVlCgoKXH2Sk5NNy5YtTVpamvn6669Nnz59zBVXXGFOnz7toZGc3cCBA03Hjh3Nxo0bzcaNG03Hjh3N0KFDK10mJyfH7fb6668bm81mfvjhB1efXr16mbvuusut37Fjxzw9nHLOZ3xJSUlm4MCBbrUfPnzYrU9tWYfVHd+xY8dM//79zaJFi8y3335rMjIyTJcuXUxCQoJbP1+tv/fee880aNDA/OMf/zC7du0yEydONI0bNzb79u2rsP+PP/5oGjVqZCZOnGh27dpl/vGPf5gGDRqYxYsXu/ps3LjR+Pv7myeffNLs3r3bPPnkkyYgIMB8/vnnHh/Pmao7vokTJ5qnn37afPnll+a7774z06dPNw0aNDBff/21q88bb7xhQkNDy22XvlLdMa5du9ZIMnv27HGr//fbkpXX4bFjx9zGlZWVZZo1a2Yee+wxV5/atA5XrlxpHn74YbNkyRIjySxbtqzS/rVlGyS8VOCNN96oUngpLS01kZGR5qmnnnK1nTp1ytjtdjN37lxjjPMXuUGDBua9995z9cnOzjZ+fn5m9erVNV57ZXbt2mUkuf0CZWRkGEnm22+/rfLz3HjjjaZv375ubb169TITJ06sqVLPy/mOLykpydx4441nfby2rMOaWn9ffvmlkeT25uur9XfNNdeY5ORkt7Z27dqZadOmVdj/oYceMu3atXNru/vuu03Xrl1d92+99VYzcOBAtz4DBgwwo0aNqqGqq66646tI+/btzcyZM133q/r+5C3VHWNZeDl69OhZn7MurcNly5YZm81mfvrpJ1dbbVuHZaoSXmrLNshhowuQmZmp3NxcJSYmutqCgoLUq1cvbdy4UZK0efNmFRcXu/WJjo5Wx44dXX28JSMjQ3a7XV26dHG1de3aVXa7vcq1/PLLL1qxYoXuvPPOco+9/fbbCgsLU4cOHTRlyhTXN3x7y4WMLz09XeHh4frDH/6gu+66S3l5ea7Hass6rIn1J0n5+fmy2WzlDo16e/0VFRVp8+bNbv+vkpSYmHjW8WRkZJTrP2DAAG3atEnFxcWV9vH29nY+4ztTaWmpCgsL1axZM7f248ePKy4uTjExMRo6dKi2bNlSY3VXx4WMsVOnToqKilK/fv20du1at8fq0jqcP3+++vfvr7i4OLf22rIOq6u2bIN17osZvSk3N1eSFBER4dYeERGhffv2ufoEBgbqoosuKtenbHlvyc3NVXh4eLn28PDwKteyYMEChYSEaPjw4W7tY8aMUevWrRUZGakdO3Zo+vTp2rZtm9LS0mqk9qo43/ENGjRII0aMUFxcnDIzM/XII4+ob9++2rx5s4KCgmrNOqyJ9Xfq1ClNmzZNo0ePdvtCNV+sv0OHDqmkpKTC7eds48nNza2w/+nTp3Xo0CFFRUWdtY+3t7fzGd+ZnnvuOZ04cUK33nqrq61du3ZKSUnRZZddpoKCAr3wwgvq0aOHtm3bpksuuaRGx3Au5zPGqKgozZs3TwkJCXI4HHrrrbfUr18/paen67rrrpN09vVstXWYk5OjVatW6Z133nFrr03rsLpqyzZY58PLjBkzNHPmzEr7fPXVV+rcufN5v4bNZnO7b4wp13amqvSpqqqOUSpfa3Vref311zVmzBg1bNjQrf2uu+5y/dyxY0ddcskl6ty5s77++mtdddVVVXrus/H0+EaOHOn6uWPHjurcubPi4uK0YsWKciGtOs9bVd5af8XFxRo1apRKS0v1yiuvuD3myfV3LtXdfirqf2b7+WyTnnK+tbz77ruaMWOG/vWvf7mF1q5du7pNKO/Ro4euuuoqvfjii/r73/9ec4VXQ3XG2LZtW7Vt29Z1v1u3bsrKytJf//pXV3ip7nN62vnWkpKSoqZNm+qmm25ya6+N67A6asM2WOfDy4QJE8551kR8fPx5PXdkZKQkZxKNiopytefl5blSZ2RkpIqKinT06FG3v9zz8vLUvXv383rdM1V1jN98841++eWXco8dPHiwXEquyIYNG7Rnzx4tWrTonH2vuuoqNWjQQHv37r3gDz9vja9MVFSU4uLitHfvXkmeX4feGF9xcbFuvfVWZWZm6pNPPjnn19jX5Po7m7CwMPn7+5f7a+z328+ZIiMjK+wfEBCg5s2bV9qnOr8DNeF8xldm0aJFuvPOO/X++++rf//+lfb18/PT1Vdf7fp99aYLGePvde3aVQsXLnTdrwvr0Bij119/XWPHjlVgYGClfX25Dqur1myDNTZ7pg6p7oTdp59+2tXmcDgqnLC7aNEiV58DBw74dMLuF1984Wr7/PPPqzzhMykpqdxZKmezfft2I8msW7fuvOutrgsdX5lDhw6ZoKAgs2DBAmNM7VmH5zu+oqIic9NNN5kOHTqYvLy8Kr2Wt9bfNddcY+655x63tksvvbTSCbuXXnqpW1tycnK5yYKDBg1y6zNw4ECfTfaszviMMeadd94xDRs2POfEyTKlpaWmc+fO5vbbb7+QUs/b+YzxTH/84x9Nnz59XPetvg6N+c/E5O3bt5/zNXy9DsuoihN2a8M2SHj5nX379pktW7aYmTNnmiZNmpgtW7aYLVu2mMLCQleftm3bmqVLl7ruP/XUU8Zut5ulS5ea7du3m//6r/+q8FTpmJgY8/HHH5uvv/7a9O3b16enSl9++eUmIyPDZGRkmMsuu6zcqbZnjtEYY/Lz802jRo3Mq6++Wu45v//+ezNz5kzz1VdfmczMTLNixQrTrl0706lTJ5+cSlyd8RUWFpoHH3zQbNy40WRmZpq1a9eabt26mZYtW9bKdVjd8RUXF5sbbrjBxMTEmK1bt7qdlulwOIwxvl1/Zaehzp8/3+zatctMmjTJNG7c2HVmxrRp08zYsWNd/ctO03zggQfMrl27zPz588udpvnZZ58Zf39/89RTT5ndu3ebp556yuen2VZ1fO+8844JCAgwL7/88llPW58xY4ZZvXq1+eGHH8yWLVvM7bffbgICAtxCrTdVd4x/+9vfzLJly8x3331nduzYYaZNm2YkmSVLlrj6WHkdlrnttttMly5dKnzO2rQOCwsLXZ91kszzzz9vtmzZ4jobsbZug4SX30lKSjKSyt3Wrl3r6iPJvPHGG677paWl5rHHHjORkZEmKCjIXHfddeWS9smTJ82ECRNMs2bNTHBwsBk6dKjZv3+/l0bl7vDhw2bMmDEmJCTEhISEmDFjxpQ7ZfHMMRpjzGuvvWaCg4MrvPbH/v37zXXXXWeaNWtmAgMDzcUXX2zuv//+ctdK8Ybqju/XX381iYmJpkWLFqZBgwamVatWJikpqdz6qS3rsLrjy8zMrPB3+ve/175efy+//LKJi4szgYGB5qqrrnLb25OUlGR69erl1j89Pd106tTJBAYGmvj4+AoD9fvvv2/atm1rGjRoYNq1a+f2weht1Rlfr169KlxXSUlJrj6TJk0yrVq1MoGBgaZFixYmMTHRbNy40YsjKq86Y3z66afNxRdfbBo2bGguuugic+2115oVK1aUe06rrkNjnHtrg4ODzbx58yp8vtq0Dsv2EJ3td662boM2Y36baQMAAGABXOcFAABYCuEFAABYCuEFAABYCuEFAABYCuEFAABYCuEFAABYCuEFAABYCuEFAABYCuEFAABYCuEFAABYCuEFAABYCuEFAABYyv8HCZEXIrMjqX0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sphere_point = lambda rad: np.array([np.cos(rad), np.sin(rad)])\n",
    "\n",
    "square_lattice = Lattice([1, 0], [0, 1])\n",
    "tilted_lattice = Lattice([1, 1], [0, 1])\n",
    "honeycomb_lattice = Lattice(sphere_point(0), sphere_point(np.pi * 2/3))\n",
    "\n",
    "def get_elementary_cell(lattice: Lattice):\n",
    "    points = []\n",
    "    for n1 in range(-1, 2):\n",
    "        for n2 in range(-1, 2):\n",
    "            p = lattice.get_point(n1, n2)\n",
    "            p_len = np.linalg.norm(p)\n",
    "            # if (p_len <= lattice.vec_lengths).any():\n",
    "            points.append(p)\n",
    "    return points\n",
    "\n",
    "def get_n_points_around_center(n)\n",
    "\n",
    "points = get_elementary_cell(tilted_lattice)\n",
    "vor = Voronoi(points)\n",
    "\n",
    "print(vor)\n",
    "fig = voronoi_plot_2d(vor)\n",
    "\n",
    "def get_orthogonal_2D(vec):\n",
    "    return np.array((vec[1], -vec[0]))\n",
    "def plot_voronoi(lattice, v: Voronoi, subplot_kw={}):\n",
    "    fig, ax = plt.subplots(**subplot_kw) \n",
    "    lattice_points = v.points\n",
    "    lattice_vec_norm = np.sqrt(np.sum(lattice.vec_lengths**2))\n",
    "    print(lattice_vec_norm, np.linalg.norm(v.vertices, axis=1))\n",
    "    cell_vertices = v.vertices[np.linalg.norm(v.vertices, axis=1) <= 0.5 * lattice_vec_norm]\n",
    "    \n",
    "    x, y = zip(*lattice_points)\n",
    "    ax.scatter(x, y, color=\"blue\")\n",
    "    x, y = zip(*cell_vertices)\n",
    "    arrowprops = dict(arrowstyle=\"-|>,head_width=0.4,head_length=0.8\", color=\"black\", shrinkA=0,shrinkB=0)\n",
    "    for i, vec in enumerate(lattice.vecs):\n",
    "        ax.annotate(f\"\", xy=lattice.vecs[i], xytext=lattice.center, arrowprops=arrowprops)\n",
    "        # add name of vector at a perpendicular offset starting at half length\n",
    "        ax.annotate(r\"$\\vec{a}_\"+f\"{i+1}$\", xy=0.7*lattice.vecs[i], xytext=0.7*lattice.vecs[i] + 0.06*get_orthogonal_2D(lattice.vecs[i]))\n",
    "        \n",
    "    ax.scatter(x, y, color=\"orange\")\n",
    "    ax.fill(x, y, color=\"#4444\")\n",
    "    return fig\n",
    "    \n",
    "\n",
    "plot_voronoi(tilted_lattice, vor)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2faef81-5f2a-4950-b986-e891ddfa7da8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6b6ccf4-34c7-419d-8bcc-f3ce2870ba25",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "conda",
   "language": "python",
   "name": "conda"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}