radial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from icecream import ic\n",
    "from scipy.sparse import csr_matrix\n",
    "from scipy.sparse.linalg import LinearOperator\n",
    "from scipy.linalg import eigh\n",
    "\n",
    "%matplotlib widget\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def radial_fdm_laplacian(r_max, n, left_bc = 'neumann_simple'):\n",
    "    \"\"\"Generate laplacian in 2d polar coordinates, radial part.\n",
    "    \n",
    "    Discretize [0, r_max] with n+2 points, and apply\n",
    "    neumann conditions on r=0 and Dirichlet on r=r_max, \n",
    "    resulting in an nxn matrix for the inner grid points.\n",
    "    \n",
    "    Args:\n",
    "        r_max (float): right endpoint\n",
    "        n (int): number of inner grid points\n",
    "        left_bc (str, optional): 'neumann' for fancy neumann BC, 'neumann_simple' for simple Neumann, else dirichlet.\n",
    "    \n",
    "    Returns:\n",
    "        ndarray: Laplacian as CRS matrix\n",
    "        r: radial grid with endpoints\n",
    "        G: transition matrix from inner to global grid\n",
    "    \n",
    "    \"\"\"\n",
    "    r = np.linspace(0, r_max, n+2)\n",
    "    L = np.zeros((n+2,n+2))\n",
    "    h = r[1]-r[0]\n",
    "    \n",
    "    for i in range(n+2):\n",
    "        if i < n+1:\n",
    "            L[i,i+1] = (r[i]+h/2)\n",
    "        if i > 0:\n",
    "            L[i,i-1] = (r[i]-h/2)\n",
    "        if i > 0 and i < n+1:\n",
    "            L[i,i] = -2*r[i]\n",
    "        elif i == 0:\n",
    "            L[i,i] =  -(r[i]+h/2)\n",
    "        else:\n",
    "            L[i,i] =  -(r[i]+h/2)\n",
    "#        L[i,i] = -2*r[i] #-(r[i] + h/2)/r[i] - (r[i] - h/2)/r[i]\n",
    "    # for i in range(n+2):\n",
    "    #     if i < n+1:\n",
    "    #         L[i,i+1] = (r[i]+h/2)\n",
    "    #     if i > 0:\n",
    "    #         L[i,i-1] = (r[i]-h/2)\n",
    "    #     L[i,i] = -2 #-(r[i] + h/2)/r[i] - (r[i] - h/2)/r[i]\n",
    "    \n",
    "    G = np.zeros((n+2,n))\n",
    "    G[1:-1,:] = np.eye(n)\n",
    "    if left_bc == 'neumann_simple':\n",
    "        G[0,0] = 1.0\n",
    "\n",
    "    elif left_bc == 'neumann':\n",
    "        G[0,0] = 1.0\n",
    "        G[0,0] = 4/3\n",
    "        G[0,1] = -1/3\n",
    "    else:\n",
    "        pass\n",
    "    \n",
    "    L1 =  np.diag(r[1:-1]**(-.5)) @ (G.T @ (L @ G) / h**2) @ np.diag(r[1:-1]**(-.5))\n",
    "    return L1, r, G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ic| E[:]: array([  1.00120895,   3.00107011,   4.99589775,   6.98568789,\n",
      "                   8.97043104,  10.95011503,  12.924726  ,  14.89424888,\n",
      "                  16.85866766,  18.81796555,  20.7721251 ,  22.72112827,\n",
      "                  24.66495647,  26.60359063,  28.53701119,  30.46519815,\n",
      "                  32.38813114,  34.30579014,  36.21816192,  38.12528535,\n",
      "                  40.02750257,  41.92649879,  43.82825563,  45.74772789,\n",
      "                  47.70981883,  49.74078322,  51.85719885,  54.06375164,\n",
      "                  56.35788363,  58.73433294,  61.18739255,  63.71166991,\n",
      "                  66.30222998,  68.95454703,  71.66441814,  74.42788387,\n",
      "                  77.24116522,  80.10061606,  83.00268781,  85.94390331,\n",
      "                  88.92083753,  91.93010329,  94.96834077,  98.03220978,\n",
      "                 101.11838423, 104.22354809, 107.34439272, 110.47761507,\n",
      "                 113.61991663, 116.76800298, 119.91858373, 123.06837278,\n",
      "                 126.21408885, 129.35245611, 132.4802049 , 135.5940725 ,\n",
      "                 138.69080376, 141.76715172, 144.81987797, 147.84575278,\n",
      "                 150.84155489, 153.80407083, 156.73009369, 159.61642119,\n",
      "                 162.45985278, 165.25718572, 168.00520955, 170.70069886,\n",
      "                 173.34040354, 175.92103595, 178.43925386, 180.89163781,\n",
      "                 183.27466062, 185.58464603, 187.81771149, 189.96968724,\n",
      "                 192.03599882, 194.01148972, 195.89014081, 197.66459593,\n",
      "                 199.32528064, 200.8585028 , 202.24146609, 203.43681304,\n",
      "                 204.48865106, 205.64357306, 207.00651728, 208.55193706,\n",
      "                 210.26404835, 212.13959348, 214.18229076, 216.4013307 ,\n",
      "                 218.81150347, 221.43443692, 224.30120643, 227.45732028,\n",
      "                 230.97269189, 234.96406026, 239.65662729, 245.63276015])\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3b858c237ad24a3181acb33246173b2c",
       "version_major": 2,
       "version_minor": 0
      },
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABV9UlEQVR4nO3dd3xUVf7/8dedSSWkkIQ0SCB0pAgEpShixe5aVsGCumtdZS24rn0tv3VZ67q7LvbyXSvuWnEVxUbvTaT3hBRCAmRCSJ25vz+GCUYChGT6fT8fj3nMZObeez7JsjvvPeeecwzTNE1ERERExDJsgS5ARERERPxLAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYiICXUAoc7lcFBUVER8fj2EYgS5HREREWsA0TSorK8nKysJms2ZfmAJgGxQVFZGdnR3oMkRERKQVCgoK6Ny5c6DLCAgFwDaIj48H3P+AEhISAlyNiIiItITD4SA7O7vxe9yKFADbwDPsm5CQoAAoIiISYqx8+5Y1B75FRERELEwBUERERMRiFABFRERELEYBUERERMRiFABFRERELEYBUERERMRiFABFRERELEYBUERERMRiFABFRERELEYBUERERMRiFABFRERELEYBUERERMRiFACD0OoiBxdNnsO3a3ZgmmagyxEREZEwowAYhF6csYll+Xu47v8Wc84/ZvPFymJcLgVBERER8Q4FwCD00HnHcNPobsRF2VlT7OCWd5Yy5rmZTF1RFOjSREREJAwYpsYYW83hcJCYmEhFRQUJCQlev/7uqjremLuVN+ZsobKmAYBrR3blofOOwW4zvN6eiIiIFfj6+zsUqAcwSG2p2IJh38fEM3ox595TmXBKDwDenLuVG/69mL21DQGuUEREREKVAmAQenLRk1zwyQX8Z/1/AEiIieQPZ/Zm8pVDiI6w8d3aUi59cR7FFdUBrlRERERCkQJgEOqb3BeAKeum0OA60NN3zoBMptw0gtT20awpdvCr5+ewqqgiUGWKiIhIiFIADEJndj2T5JhkduzbwfcF3zf5bFB2Ep/cOpJe6e0prazlmtcXUbhHPYEiIiLScgqAQSjKHsUlPS8B4N017x70eecO7fjv70bSJyOesr21XPfmIt0TKCIiIi2mABikLut9GXbDzuIdi1m/e/1BnyfERPLatceR2j6atSWV3PH+MpxaK1BERERaQAEwSGXEZXBqzqkAvLf2vWaP6ZQUyytX5xEVYeObNaX89cs1/ixRREREQpQCYBC7os8VAPxv8/+oqG1+ssfgnA48c+mxALwyawvvLcz3W30iIiISmhQAg1heeh49O/SkuqGaTzZ+csjjzj82iztP7wXAQ5/8xPKCPf4pUEREREJSWAXAyZMnk5ubS0xMDHl5ecyaNatF582ZM4eIiAgGDRrk2wKPkmEYjb2A7699H5fpOuSxt53Wg3MHZNLgMrnrg+XU1Dv9VaaIiIiEmLAJgFOmTOGOO+7ggQceYNmyZYwaNYqzzz6b/PzDD4lWVFRw9dVXc9ppp/mp0qNzbrdziY+KZ/ve7cwunH3I4wzD4PGL+tMxPppNO6t45ut1fqxSREREQknYBMBnn32W6667juuvv56+ffvy3HPPkZ2dzQsvvHDY82666SauuOIKRowY4adKj05sRCwX97gYgHfXHrwkzM8ltYvirxcPAODV2VtYtHWXz+sTERGR0BMWAbCuro4lS5YwZsyYJu+PGTOGuXPnHvK8N954g02bNvHwww+3qJ3a2locDkeThz+M7TMWA4M5hXMo3lt82GNP65vOr/M6Y5rwh/+sYF+d1gcUERGRpsIiAJaVleF0OklPT2/yfnp6OiUlJc2es2HDBu69917eeecdIiIiWtTOpEmTSExMbHxkZ2e3ufaWyI7PZkCqu2dvYcnCIx7/p/OPITMxhm3l+3jiy7W+Lk9ERERCTFgEQA/DMJr8bJrmQe8BOJ1OrrjiCh599FF69erV4uvfd999VFRUND4KCgraXHNL5WXkAbB4x+IjHpsQE8kTlwwE4P/mbWPuxjKf1iYiIiKhJSwCYGpqKna7/aDevtLS0oN6BQEqKytZvHgxEyZMICIigoiICB577DFWrFhBREQE3333XbPtREdHk5CQ0OThL8elHwfA4pIjB0CAk3p15IphOQDc9/FKahs0K1hERETcwiIARkVFkZeXx/Tp05u8P336dEaOHHnQ8QkJCaxcuZLly5c3Pm6++WZ69+7N8uXLGTZsmL9Kb7HBaYOxGTa2791OSVXzw9q/dP85fUmLj2Zb+T7emLPVtwWKiIhIyAiLAAgwceJEXn31VV5//XXWrFnDnXfeSX5+PjfffDPgHr69+uqrAbDZbPTv37/JIy0tjZiYGPr3709cXFwgf5VmtY9qT9/kvkDLhoEB2kdHcM9ZfQD457cbKK2s8Vl9IiIiEjrCJgCOHTuW5557jscee4xBgwYxc+ZMvvjiC7p06QJAcXHxEdcEDHZD04cCLR8GBrhocCeOzU6iqs7JU9O0NqCIiIiAYZqmGegiQpXD4SAxMZGKigq/3A/4Q8EP/P6739M1oStTL5ra4vOW5u/m4snu5XA+vfUEjs1O8k2BIiIiIcDf39/BKGx6AK1gSPoQDAy2Orayc9/Olp+X04GLB3cC4NGpq1DmFxERsTYFwBCSEJVAn2T3PX1Ldiw5qnP/eFYf2kXZWZq/h89WFPmiPBEREQkRCoAhJi/dvR7gopJFR3VeRmIMt57SA4BJX6zVDiEiIiIWpgAYYoZm7J8I0sKZwD933Ym5dO4QS4mjhtdnb/F2aSIiIhIiFABDTF6auwdwc8VmyqvLj+rcmEg7d5/ZG4CXZm6mYl+91+sTERGR4KcAGGKSYpLo2aEncPT3AQKcPzCL3unxVNY08Mqszd4uT0REREKAAmAIalwPsBXDwDabwZ1nuPc/fn3OFsr31nq1NhEREQl+CoAh6LiM/fsCtyIAApzZL50BnRLZV+fkxRmbvFmaiIiIhAAFwBDkmQm8YfcG9tTsOerzDcPgrjHuXsB/z9vGDoe2iBMREbESBcAQlByTTPfE7gAsKT36+wABRvfqyHFdO1Db4OL57zZ6szwREREJcgqAIapxOZij2Bf459y9gO4Zwe8vyqdg1z6v1SYiIiLBTQEwRA1JGwLAyrKVrb7G8G4pjOqZSr3T5O/fbvBWaSIiIhLkFABDlGdLuPW71+MyXa2+zsT9M4I/XlZIfrl6AUVERKxAATBEdUnoQow9huqGagoqC1p9ncE5HRjVMxWny+SlmZoRLCIiYgUKgCHKbrM3Lgi9dtfaNl1rwv49gv+zeLtmBIuIiFiAAmAI69XBPXy7bte6Nl1nWLcUjuvagTqni1dmancQERGRcKcAGMI89wGu2922AAhw6/5ewHcW5LOrqq7N1xMREZHgpQAYwjwBsK1DwOBeF7B/pwSq6528OWdLm68nIiIiwUsBMIR57gEs3VfK7prdbbqWYRjcerK7F/DNuVuprKlvc30iIiISnBQAQ1hcZBw58TmAd4aBz+yXQfeOcThqGnhr/rY2X09ERESCkwJgiOud7N7No60TQQBsNoNb9vcCvjZrC9V1zjZfU0RERIKPAmCI8+Z9gAAXDMqic4dYyqvqmLIo3yvXFBERkeCiABjienfY3wPohSFggEi7jZtO6gbAq7O30OBs/S4jIiIiEpwUAEOcZwh4y54t1DprvXLNX+dlkxwXxfbd1Xz5U4lXrikiIiLBQwEwxKW3SycpOokGs4FNe7yzlVtslJ2rR3QB4OWZmzFN0yvXFRERkeCgABjiDMM4MAzshYkgHleP6EpMpI2VhRXM21zuteuKiIhI4CkAhoHGmcBeug8QIDkuikvzsgF3L6CIiIiEDwXAMODtmcAe14/KxWbAD+t2sq6k0qvXFhERkcBRAAwDvTr0AmD9rvVevV+vS0ocZ/XPANQLKCIiEk4UAMNAt8RuRNoiqayvpKiqyKvXvvGk7gB8tqKQkooar15bREREAkMBMAxE2iPpkeTewcPbw8CDspM4PjeZeqfJG3O2ePXaIiIiEhgKgGHi58PA3uZZGPrdBflU1tR7/foiIiLiXwqAYcJXE0EATumdRo+09lTWNvDB4u1ev76IiIj4lwJgmPDFUjAeNpvBdSfmAvC6tocTEREJeQqAYcIzBFy4txBHncPr179ocCeS46Io3FPNV6t2eP36IiIi4j8KgGEiMTqRrLgswDf3AcZE2rlquHt7uFdna0kYERGRUKYAGEa6J7mXbNlc4ZuANn54F6LsNpbl72HJtt0+aUNERER8TwEwjOQmuu/T21Lhm+VaOsZHc+Fgdy/jq7PUCygiIhKqFADDSLdE93ItvuoBBLh+lLuNr1aVkF++z2ftiIiIiO8oAIaRbkm+D4C90uM5qVdHXCa8MVcLQ4uIiIQiBcAw4ukBLKkqYV+973rnrt+/JMwHiwqoqNbC0CIiIqFGATCMJEYnkhyTDMAWh+9650b1TKVXenuq6py8vzDfZ+2IiIiIbygAhpnG+wD3+G4Y2DAMrj/R3c7/zd1KvRaGFhERCSkKgGHGEwB9NRPY44JBWaS2j6KoooYvfyrxaVsiIiLiXQqAYcazFIwvJ4KAe2Ho8cO7Au4lYUzT9Gl7IiIi4j0KgGHGH0vBeFw1PIeoCBs/bq9gsRaGFhERCRkKgGHGsxRMgaOAepdvZ+imtI/mkiGdAC0MLSIiEkoUAMNMert02kW0o8FsoMBR4PP2fnuCe8j569U72FZe5fP2REREpO0UAMOMYRg+3xLu53qmx3Ny746YJrwxZ6vP2xMREZG2UwAMQ/68DxBoXBLmg8UFVOzTwtAiIiLBTgEwDPljS7ifO6FHCn0y4tlX5+S9RVoYWkREJNgpAIah3AT/LAXjYRgG1+3fHu7NOVupa9DC0CIiIsFMATAM5SYduAfQZfonjLkXho6mxFHD/1YW+aVNERERaR0FwDCUHZ9NhBFBdUM1O6p2+KXN6Ag7147sAsArM7doYWgREZEgpgAYhiJtkeQk5AD+GwYGuHJYF2Ij7awudjB3U7nf2hUREZGjowAYpvw9ExigQ1wUlw3tDMDLM7UwtIiISLBSAAxT/lwL8Od+e2IuNgNmrN/JupJKv7YtIiIiLaMAGKY8AdCfPYAAXVLiOKt/BgCvaHs4ERGRoKQAGKY8awH6uwcQ4IZR7rY/XV7IDkeN39sXERGRw1MADFOetQB31exiT80ev7Y9OKcDx3XtQL3T5M25W/3atoiIiByZAmCYahfZjsy4TMD/w8AA1+/vBXxn/jaqahv83r6IiIgcmgJgGAvETGCP0/umk5sah6OmgSmLCvzevoiIiByaAmAYC9REEAC77cD2cK/N3kK9U9vDiYiIBAsFwDDmmQgSiAAI8Ou8zqS2j6JwTzVTV2h7OBERkWChABjGPBNBtlZsDUj7MZF2fnOCu4YXZ2zC5dL2cCIiIsEgrALg5MmTyc3NJSYmhry8PGbNmnXIYz/66CPOOOMMOnbsSEJCAiNGjOCrr77yY7W+1zWxKwBFe4uoc9YFpIarhnehfXQE63fs5ft1pQGpQURERJoKmwA4ZcoU7rjjDh544AGWLVvGqFGjOPvss8nPz2/2+JkzZ3LGGWfwxRdfsGTJEk455RTOP/98li1b5ufKfSclJoX2ke0xMSmoDMxEjMTYSK4c7t6X+IUfNgWkBhEREWnKME0zLMblhg0bxpAhQ3jhhRca3+vbty8XXnghkyZNatE1+vXrx9ixY/nTn/7UouMdDgeJiYlUVFSQkJDQqrp9bdzn41hVvornTnmO03JOC0gNpY4aTnzie+qcLv5z8wiO65ockDpEREQgNL6/fS0segDr6upYsmQJY8aMafL+mDFjmDt3bouu4XK5qKysJDk5vMJJl4QuQODuAwRIS4jhkrxOgHoBRUREgkFEoAvwhrKyMpxOJ+np6U3eT09Pp6SkpEXXeOaZZ6iqquKyyy475DG1tbXU1tY2/uxwOFpXsB91TegKwDbHtoDWceNJ3Xl/UQHfrS1lbYmDPhlt/H9cNQ7YOhucdcDPOrFTekJG/7ZdW0REJMyFRQD0MAyjyc+maR70XnPee+89HnnkET799FPS0tIOedykSZN49NFH21ynP3kmggQ6AOamxnFO/0z+t7KYl2Zs5m9jB7XuQg21sOg1mPU07Ctv/piuo2DErdDzTLCFRSe3iIiIV4XFt2Nqaip2u/2g3r7S0tKDegV/acqUKVx33XV88MEHnH766Yc99r777qOioqLxUVAQ/DtcNA4BO7YGthDg5tHdAfhsRREFu/Yd3ckuJyx/D/45FL66zx3+knIgZ6T70eUEyB4GtgjYOgveGwf/Os4dFhtqj3x9ERERCwmLHsCoqCjy8vKYPn06F110UeP706dP51e/+tUhz3vvvff47W9/y3vvvce55557xHaio6OJjo72Ss3+4gmAu2p24ahzkBAVuJtdB3ROZFTPVGZtKGPyD5uYdPGAlp24dye8fTGU/Oj+OT4TTr4XBl0F9l/8E64ohIUvweI3oXwj/G8iLH8Xxr4NCZle/X1ERERCVVj0AAJMnDiRV199lddff501a9Zw5513kp+fz8033wy4e++uvvrqxuPfe+89rr76ap555hmGDx9OSUkJJSUlVFRUBOpX8Im4yDg6xnYEYFtFYIeBAW47rScA/11SQOGe6iOfUFcF717mDn8xiXD6I/D7pZB37cHhDyCxE5zxGExcDWc9ATFJULgYXh4NBQu9+auIiIiErLAJgGPHjuW5557jscceY9CgQcycOZMvvviCLl3cPWDFxcVN1gR86aWXaGho4NZbbyUzM7PxcfvttwfqV/AZz32AwTAMfFzXZEZ0S6HeafLikWYEu5zw4fVQtBRiO8D138KJd0JUuyM3FN0eht8MN34PacfA3h3wxjmw5P+884uIiIiEsLBZBzAQQmUdoUfnPcp/1/+XmwbexITBEwJdDvM2lXP5K/OJstuY+cdTyEiMOfgg04Qv/gCLXgV7NFzzGeQMb12DtXvhk5thzVT3z8fdAGc/qQkiIiIWFSrf376kb0ALCJalYDyGd0vm+K7J1DldvDjjEL2Ac//pDn8YcPHLrQ9/4O4NvPTfcMqD7p8XvQJf3u0OmSIiIhakAGgBnokgwRIADcNovBfwvYX5lDpqmh6w6hOY/pD79ZmPQ78L296ozQaj74aLXwEMd7j8/vG2X1dERCQEKQBawM+XggmWEf8TeqQwJCeJ2gYXL8/cfOCDmgr3zF2AYTfD8Fu82/DAy+Dcp92vZz4Fc5/37vVFRERCgAKgBXSO74zdsFPdUE3pvtJAlwM07QV8e8E2yvbuX6tv9t/ca/yl9IQxf4YWLOR91I67Hk7d38P49QOw9C3vtyEiIhLEFAAtINIWSef4zkDwDAMDjO7VkWM7J1JT7+KVmZthTwHMm+z+8IzHwB7pu8ZH3QUj9k+ImXobrPncd22JiIgEGQVAiwimHUE8ft4L+O9526j56mFw1kKXE6H32b5u3N3DOHg8mC746EbYsdq3bYqIiAQJBUCLCLaJIB6n9kljcE4SPRvWE7PmQ/ebZ/po6PeXDAPOew5yT4L6Knj/Cqje7ft2RUREAkwB0CI8S8EEUw8guHsB7x7Tiwci3wGgqvclkDXYfwXYI+DXb0JiDuzeAh/e4F6AWkREJIwpAFpEsK0F+HMj6xcwzLaWGjOSZ1xj/V9AXAqMfQsiYmDjdPhhkv9rEBER8aNmNlOVcOQZAt5euZ16Vz2RNh9OsDgaznqY/icAXnOezRs/NTBuRyW90uP9W0fWIDj/H/Dxje7lYTKPhb7nt/h00zSZVzSPLY4tOOocVNZV4qh14DSdjOo0ilNyTiE2ItZ39YuIiBwFBUCLSGuXRmxELNUN1RRWFjbuDxxwP30EuzZBu1Q2pN+AuaaSZ75ex0vjh/q/lmPHQvFymD8ZPr4ZOvaB1J5HPG1V+SqeXPgkS0uXNvv555s/Jy4yjjFdxnB+9/PJS8/DZqjzXUREAkcB0CIMw6BLQhfW7lrLVsfW4AmAi193Pw+7iQl9B/PZ2pl8tWoHywv2MCg7yf/1nPEYFP8I22bDf38D130Dkc3sVQyUVZfxz2X/5OMNH2NiEmOPYVTnUSRGJ5IQlUB8VDx76/Yybes0CvcW8vHGj/l448f0SOrBpFGT6JPcx8+/nIiIiJsCoIV0TejK2l1rg+c+wB2roWA+GHYYPJ4eCfFcPKQz/12ynae/Wsfb1w/zf032SLjkVXjxBChZ6d6S7pynDjrs882f8+f5f6aqvgqAc3LP4c68O8mIyzjo2NuG3May0mVM3TSVr7Z+xcY9G7n8f5dz2+DbuKbfNeoNFBERv9M3j4UE3VqAS95wP/c+GxIyAbj9tJ5E2g1mbyxj9oaywNSVkAkXveR+vfBlWDO1ycezts/iwdkPUlVfRb+Ufrx19ls8cdITzYY/AJthIy89j0dGPsL/Lv4fp2afSoOrgWeXPMv1X19PSVWJr38jERGRJhQALSSo1gKs2wcrprhfD/1N49vZye24ari7zj//bzVOV4D2Lu55Boz8vfv1p7fCnnwA1pSv4a4Zd+E0nVzQ/QLePfddBqUNavFlk2OSee6U53h05KPERsSyqGQRF392MbO2z/LBLyEiItI8BUALaVwKpiIIAuCqj6C2ApK6QLdTm3x0+2k9SYyNZG1JJVMWFQSoQODUP0GnPKipgA+vp9iRz63f3kp1QzXDMofxyIhHWjV8axgGF/e8mP+c/x8GpA6gsq6S276/ja+3fu2DX0JERORgCoAW0iXR3bNWWl3aeO9awHgmf+RdC7am/wyT2kVxx+nu2bfPfL0OR029n4vbLyIKfv06RCdQuX0ht0y9nJ3VO+mR1INnT36WyDbuVdwloQv/d/b/cXbXs2lwNXD3zLv5dOOnXipeRETk0BQALSQhKoHkmGQgwMPAxT9C4RKwRcDgq5o95KrhXejWMY7yqjr+9d1GPxf4Mx26Un/e35iYnsrGBgcdoxKZfNpkEqISvHL5SFskk0ZN4uKeF+MyXTw450HeX/u+V64tIiJyKAqAFtO4JVzF1sAV4Zn80fd8aJ/W7CGRdhsPnXsMAG/M2cq28sD1WL5rq2J+bCyxLhf/2lFOpi3aq9e32+w8POJhrux7JQCPL3ic13963attiIiI/JwCoMUEfCJIbSX8+IH7dd5vDnvoyb07MqpnKnVOF5O+WOuH4g62q2YXL61wzwi+pzaSvnuK4LPfg+ndySk2w8Y9x93DDQNuAOBvS/6mnkAREfEZBUCL8SwAHbClYFb+F+r2QkoPyD3psIcahsFD5x2DzYBpq0qYt6ncT0UeMHn5ZCrrK+mT3IcLz38dbJGw9vMDvZheZBgGtw25jd8d+zsA/rLgL5oYIiIiPqEAaDEB7wFc+m/3c961YBhHPLxXejxXDnPX/OjUVTQ4XT4srqkNuzfwn/X/AeCPx/0Re6chcPoj7g+n3Q+lvumV/N2xv+OyXpdhYnLvrHtZVLLIJ+2IiIh1KQBaTOM9gI6tmF4exjyiPflQtBQwYODYFp925xm9GpeFeWPOVp+V93OmafLkoidxmS7O6HIGx2Uc5/5g+C3Q/VRoqIYPr4P6Gq+3bRgG9w+7n9NzTqfeVc9t393Gul3rvN6OiIhYlwKgxWTHZ2MzbFTVV1Fe4+ch1bX/cz/njDjk5I/mJMdFcf857n1zn52+nu279/miuiZmbJ/B/OL5RNoiuTPvzgMf2Gxw4YvQLhV2/ATfPOyT9u02O3896a/kpeext34vN39zM9srt/ukLRERsR4FQIuJskeRFZcFwJaKLf5t3LOlWt/zjvrUy4Zmc3xuMtX1Tv706Sqf9l7WO+t5evHTAIw/ZjzZ8dlND4hPhwtfcL9e8CKs+9IndUTbo/nHqf+gZ4eelFWXceu3t+Koc/ikLRERsRYFQAvyTATx632AVWWQP8/9us/RB0DDMPjLRf2JtBt8t7aUL3/y3f657619j22ObaTEpDTOyj1IrzEw/Fb3609uAUeRT2pJiErgxdNfJK1dGpsrNvOHH/5Ag6vBJ22JiIh1KABaUEDWAlz3BZguyBgIHbq06hI90uL53ejuADzy2Sqf7BBS66zl1ZWvAvD7wb+nfVT7Qx98+sOQeSxU74IPbwCX0+v1AKS1S+P5U58nNiKWecXz+OvCv/r//k0REQkrCoAW1LgnsD97ANd87n7ue36bLnPLKT3ITY2jtLKWp7/y/sSIL7d8ye7a3WTEZfCrHr86/MER0fDrNyCqPWybDbOe8Xo9Hn1T+vLXUX/FwGDKuim8u/Zdn7UlIiLhTwHQgjx7AvttLcAaB2z+3v26jQEwJtLO4xf2B+Ct+dtYsm13W6trZJom765xB6uxvccSYYs48kkp3eHc/cHvh0mwbZ7X6vmlU3NOZWLeRACeXPQkM7fP9FlbIiIS3hQALcjTA7i9cjv1Lu8Pox5k43Rw1kFyd+jYp82XG9kjlUuGdMY04Q//WcG+Ou/cE7di5wrW7FpDlC2KS3pe0vITjx0HA8e5h7g/vA6qfDe7+pp+1zTuG3z3jLtZv3u9z9oSEZHwpQBoQWnt0oiNiKXBbKBor28mLzTROPx7XosWf26JP513DJmJMWwpq+Lx/63xyjXfWfMOAOd2O5cOMR2O7uRzn3bvbuIohI+u99n9gIZh8OCwBzku4zj2Nezjtu9uo7za/zukiIhIaFMAtCCbYSMnPgfww0SQ+hrYsH87sz5tG/79ucR2kTx96bEAvLMgn+/W7mjT9Ur3lfLNtm8AuKLvFUd/geh4uOzfEBELm76DmU+1qZ7DibRH8uzoZ8mOz6ZwbyF3/nAndc46n7UnIiLhRwHQovy2J/CWGe69f+MzoVOeVy99Qo9UfntCLgB//O9KyvfWtvpaH6z7gAazgSFpQ+iT3Mph6vR+cP5z7tc//BU2ftvqeo4kKSaJ5097nvjIeJaVLuOxeY9pZrCIiLSYAqBFefYE9nkA9Cz+3Odc9y4aXvbHs3rTM609ZXtrue+jla0KQXXOusY9f1vV+/dzx45z73OMCR9eDxW+272jW2I3nh79NHbDzqebPuXNVW/6rC0REQkvCoAW5ZelYFxO9/p/0KrFn1siJtLOc+MGEWk3+Hr1Dv6z+OgD11dbv2JXzS7S2qVxas6pbS/qrCcOrA/4n99Ag++GZ0d2Gsk9x98DwN+W/I3v8r/zWVsiIhI+FAAtqjEAVvgwAObPg33lEJMEXU/0WTP9shKZeEZvAB6ZuooNOyqP6nzP0i/jeo8j0hbZ9oIiY9z3A8YkwvaFMO2etl/zMC7vczlje4/FxOTeWfeyqnyVT9sTEZHQpwBoUZ61AEurS6mqr/JNIxvdkyroOQbsXghWh3HjSd0Y2T2FfXVObnprSYt3CVm5cyU/lf/kXvql11Es/XIkHbrCxa8ABix+HRa+4r1rN+Oe4+/hhKwTqG6oZsK3E/wzu1tEREKWAqBFJUQlkByTDPhwGHjzD+7n7l4YVj0Cu83gH5cPJjMxhs1lVfzhgxW4XEe+H/DTTZ8CcGbXMxv/Hl7T60w4/RH36y/vOfD38IFIWyRPj36anh16UlZdxq3f3kpl3dH1hIqIiHUoAFqYT/cErt4NRcvdr7uN9v71m5HaPpoXrsojym7j69U7eGHGpsMeX++q5+ut7iVqzuvmm3sUOeH2/YtEO+GDa6D88DW1Rfuo9kw+bTJpsWls3LORO3+4k3qnHxb6FhGRkKMAaGGemcA+6QHcMgswIbU3JGR5//qHMCg7iUd/1Q+Ap79ex8z1Ow957ILiBeyu3U1yTDLHZx7vm4IMA87/O3QaCjV74L1xUFPhm7aAjLgMnj/teWIjYllQvIDH5mt5GBEROZgCoIV51gLc4tji/Yt7hju7nez9ax/B5cfnMHZoNqYJt72/jIJd+5o97sstXwIwpsuYlu3721qRMTDuXUjoBGXr3TODfdgz1zelL0+PfhqbYeOTjZ/w/PLnfdaWiIiEJgVAC/NpD2BjAPTP8O8vPfqrfgzsnMieffVc+8ZCdlc1XYqlpqGGb/PdCzWf0+0c3xcUn+4OgRGxsOlb+HQCuFw+a+6kzifx4PAHAXj5x5f596p/+6wtEREJPQqAFpab4N5FY5tjm3eHCfcUwK5NYNh8uvzL4cRE2nnxqjwyEmLYtLOK6/5vEdV1B/bnnVU4i6r6KjLjMjm247H+KSprEFz6Jhh2+PF9mP4Q+HB49tJel3Lb4NsAeGrxU3y68VOftSUiIqFFAdDCOsd3xmbYqKqvoqy6zHsX3jLD/dwpz70WXoBkJcXy7+uOJyEmgqX5e/j9e8tocLp73TzDv2flnoXN8ON/DXqfBb/6l/v1vOdhzt992tz1A67n6mOuBuDhuQ9roWgREQEUAC0tyh5Fp/adAC9vCRfA+/9+qVd6PK9ecxxRETa+WbODhz79icraSmYUuEPqObl+GP79pUGXw5g/u19/8zAsfctnTRmGwR+G/oELe1yI03Ry94y7WVi80GftiYhIaFAAtDiv7wlsmrB5fw9gEARAgONzk/nn5YOxGfDewgLumfYeda46chNz6d2hd2CKGvl79xIxAFNvg9W+G541DIOHRzzMqdmnUueqY8J3E1hUsshn7YmISPBTALQ4r28JV7oGqkrdkx06H+eda3rBmf0y+H8X9gfg++1fAXB217MxDCNwRZ3+KAy6CkyXe2bwyv/6rKkIWwRPjn6ycbeQW765hblFc33WnoiIBDcFQItrXAzaWz2AnuHfLiMhIto71/SSK4d14fdnZGCP2whAQUGfwK6RZxhwwT/g2MvdC0V/dAMsf9dnzUXbo/n7qX9ndOfR1Dhr+P23v2fm9pk+a09ERIKXAqDFefYE9noADJLh31/q1GkDhuHCWd2Jd+fs49Gpq1u0ZZzP2Ozwq8kw5Gp3T+Ant8DiN3zWXLQ9mr+d/DdOyzmNOlcdt39/uyaGiIhYkAKgxXl6ALdXbm/7tmHOetg2x/06SAOgZ/bv6TlnAvDm3K3c99FKnAENgTY47+9w/I2ACZ/fAQte8llzkfZInhr9FGd2PZMGVwN3/XAXUzdN9Vl7IiISfBQALS69XTpxkXE4TWfbF4QuXAJ1e6FdCqT3906BXlRWXcay0mUA3D/6cp659FhsBkxZXMCt7yxlX11D4Iqz2eDsJ2HEBPfPX/4RvnoAXM7Dn9dKkbZI/jrqr5zX7TwazAbun30/k5dP1rZxIiIWoQBocYZh0C2xGwCbKza37WKe4d/ck9yBJsjM3D4TE5N+Kf3IiMvgkrzO/PPyIUTaDaatKuGSF+axfXfz28b5hWG4l4c55QH3z/Oeh/evhNpKnzQXYYvg8RMf57r+1wHwwooXuH/2/dQ5645wpoiIhLrg+5YWv/MEwE0Vm9p2oSBb/uWXfij4AYCTs09ufO/cgZm8d8NwUttHsabYwa+en8PCLbsCUh/gDoGj/wiXvAYRMbD+S3j9LPfuKj5gM2zckXcHD494GLth5/PNn3Pj9BupqK3wSXsiIhIcFACFbknuALhlz5bWX6S+GrbvX1suNzD7/x5OTUMN84rmAXBK9ilNPhvaNZlPJ5xIv6wEyqvquPLV+by7ID+ww6EDfg3X/g/i0mDHT/DKqZA/32fN/brXr5l8+mTaR7ZnyY4lXP6/y1ldvtpn7YmISGApAIp3hoALl4KrHuIzoUNX7xTmRQuKF1DjrCEjLoNeHXod9HmnpFj+e/NIzh2YSb3T5P6PVzLhvWXs2RfA4dDOQ+GG7yB9gHttxTfOhh/+Ck7f3Ks4Mmsk/z7732TFZVFQWcBVX1zFO2ve0X2BIiJhSAFQ6J7YHYAtFVtwtnbSQcH+3qnsYe5hzCDzfcH3AJzc+eRDLv4cG2Xn+csH88ezemO3Gfzvx2LG/G0mM9bv9GepTSVlw2+nwcBx7mVifpgEb54Lu720cPcv9OzQkw/O/4BTsk+h3lXPXxf+lTu+v0NDwiIiYUYBUMhqn0WULYo6Vx1Fe4tadxHP8GTOcO8V5iUu09W44PEvh39/yTAMbjm5Bx/fMpLuHeMorazlmtcX8tAnPwVulnB0e7j4Jbj4FYiKd4ftF0+EH//j3nrPyxKjE/n7KX/n3uPvJdIWyXcF33Hp1Eu1h7CISBhRABTsNju5iblAK4eBXS4oWOB+nT3Mi5V5x+ry1eys3klcZBxDM4a26JyBnZP4322juHZkVwDemr+NM56dyec/FgVuSHTgZfC72e4t9mod8NH18M6vobyNk3eaYRgGV/a9krfPeZuc+ByKq4q57uvreGD2A+yu2e319kRExL8UAAVo40zgsnVQUwGR7SBjgJcrazvP8O/IrJFE2aNafF5MpJ1HLujH29cNo1NSLIV7qpnw7jLGvjSfnwoDNCTaoSv8ZhqcfB/Yo2DjNzB5OHz3Z6jz/hI2x6Qcw5TzpjC291gMDD7b9BkXfHIBn2z8RPcGioiEMAVAASA3aX8P4J5W9AB6hn87DwV7pBer8g7P8i9HGv49lBN7pvLNxNHceXovYiJtLNy6i/Ofn809//0xMOsG2iPg5Hvhd/Og+6ngrIOZT8G/hsFPH3p98ej2Ue15cPiDvH3O2/Tq0Is9tXt4aM5DXDvtWpaXLvdqWyIi4h8KgAI0nQhy1BqHf4Pv/r+ivUWs370em2FjVKdRrb5ObJSd20/vyXd3ncwFx2Zhmu4dRE5+6gfu+mAFG0t9s1jzYaX2gKs+gsvegsRsqMiH//4WJo+Alf/1ehAc2HEg75/3PhPzJhJjj2Fp6VLGfzme333zO1aVr/JqWyIi4lsKgAI0HQI+6qG9fPf6euQE3/1/nt6/wWmDSYpJavP1spJi+cflg/nwdyM4sUcqDS6TD5du54y/zeTmt5awZNsu/w6NGgYccwHcugBOvh9iEt1D8h9e5+4RXDEFGry3lE2kLZLf9P8NUy+ayiU9L8Fu2JldOJtxn4/jju/v4Keyn7zWloiI+I5h6kaeVnM4HCQmJlJRUUFCQkKgy2mTemc9x71zHE7TyTe//ob0uPSWnVi5A57pBRhw7zZ3AAkiN359I/OK53FX3l1c2/9ar19/RcEeJv+wka9W7Wh8r2dae8Yel81FgzuR0j7a620eVk0FLHjZvY1czR73e3EdYdCVkHcNJHfzanP5jnxeXPEin2/+HBP3/5T0T+nPuD7jOCv3LKLtfv79RURaIJy+v1tLAbANwu0f0Pkfn89Wx1ZeOuMlRmaNbNlJqz+FD66G9P7wuzm+LfAo7a3by6gpo2hwNTD1wql0Tezqs7Y27KjkpZmb+fzHImrqXQBE2g1O75vOOQMyObl3R+Jj/Hh/ZI0DFr4Mi16FyuID73c72R0Ge53p1bC+ac8mXlv5GtO2TqPeVQ9AUnQSF/a4kLNzz6Zvct9Drr8oIuJv4fb93RoKgG0Qbv+Abv/udr4r+I57j7+XK/te2bKTpt0P8/8FQ6+D8571bYFH6autX/GHGX+ga0JXpl401S9tOmrqmbqiiA8WFbBi+4GZwpF2gxHdUxlzTDqn9kkjKynWL/XgbID102DJG7DxW9jfS4c9yh0G+54Pvc+BuFSvNFdeXc7HGz/mg3UfUFx1IHhmx2czpssYxnQdozAoIgEXbt/fraEA2Abh9g/oH0v/wSsrX+GyXpfx0IiHWnbSy6dA0VK4+FUYeKlvCzxKD8x+gM82fcY1x1zDH477g9/bX1vi4JNlRUxfXcKmnVVNPstJbsfwbskM75bCsG4pZCXG+D4U7d4Gy95299qWrfvZBwZkDoSuo9yPLiPa3DvodDmZsX0Gn2/+nFnbZ1HjrGn8LC02jWGZwxofGXEZbWpLRORohdv3d2uEVQCcPHkyTz31FMXFxfTr14/nnnuOUaMOPfNzxowZTJw4kVWrVpGVlcUf//hHbr755ha3F27/gKZumsr9s+8nLz2PN89688gn1O2Dv2aDqwHuWAlJOT6vsaVcpotTPziV8ppyXh3zKsMyAztBZWPpXqav3sHXq0tYUbAH1y/+W5faPpr+nRIY0CmRflmJHJOZQKcOsdhtPgqFO9fBms9gzVQoXtH0M8MG6f0g81jIHOR+ZPSHyNb1Wu6r38fMwpl8vfXrg8IgQE58Dv1T+3NMyjH0S+lH35S+xEXGte73EhFpgXD7/m6NsAmAU6ZMYfz48UyePJkTTjiBl156iVdffZXVq1eTk3NwMNmyZQv9+/fnhhtu4KabbmLOnDnccsstvPfee1xyySUtajPc/gGtKl/FuM/HkRyTzIyxM458wpZZ8H/nQXwWTFwdVHsAr921lkunXkpsRCyzx80+qgWgfa2ypp7FW3czf0s58zfv4qfCCpy/TIRAVISNrint6Jbanm4d4+jcoR2dOsTSKSmGzMRY4qIjvFOQoxi2zoats9yPXc2sBWnY3AE/pSek9oSU7u4JJQmdIbETRLUssNU6a1leupwFxQtYULyAn8p/wmW6mjaFQaf2neia2JWuCV3JTcyla0JXMttnkt4uPaj+sxSR0BRu39+tETYBcNiwYQwZMoQXXnih8b2+ffty4YUXMmnSpIOOv+eee/jss89Ys2ZN43s333wzK1asYN68eS1qM9z+Ae2r38ewd909ZTPHzqRDTIfDnzDzKfcOFP0ugkvf9H2BR+HVla/y96V/5+TOJ/PP0/4Z6HIOq7rOyZoSB6sKK/ip0MFPRRVsKN1LXYPrsOclxESQGh9NavtoOraPJrV9FEntokiMjSSpXSSJse5HXHQE7aMjaBdlJy46gugI2+GHmysK3cP6RcvdvYPFy6Fq5+F/idgO7jDYPs096zgu1f1olwIxSe4hZc8jOt4dGCNiqKzfy4qdK1hdvprV5atZVb6KkqqSwzaVGptKRrsMOrbrSHJMcuOjQ0wHEqISiI+Kb3y0j2xPTEQMNkMrXonIAeH2/d0aXupCCKy6ujqWLFnCvffe2+T9MWPGMHfu3GbPmTdvHmPGjGny3plnnslrr71GfX09kZHBt6OFr7WLbEdWXBZFVUVsrthMXkze4U/ID94FoOcUumckn9DphABXcmSxUXaG5HRgSM6BwO10mRTtqWbTzr1s3lnFlrIqCvdUU7SnmsI91VTWNODY/9j8i/sLj8RmQGykndgoO9ERdmIibURF2ImKsBFttxEVYSPSnk6E/Wwi7ecQ0cmgg7mbzPoC0uq3k1aXT1ptAUl1xSTW7SDatQ+qd7sfO47cvocLG9H2dgy2xzDQFo3TFoXTHkOZPZGtETbyI0wK7C6225wU2eopM+qpM0zKqssoqy6D8pa3FYWdaOxEE0EUdiKxE2HYiMROJDbs2IjAjt0wsGPHhoEdm/vZcD8bGNj2PzyvwcAG+98xwGD/T573PA689rxqGsGbvmfgnd704OmTF2m9gV1O5pwTrw50GWEnLAJgWVkZTqeT9PSma9elp6dTUtJ8b0JJSUmzxzc0NFBWVkZmZuZB59TW1lJbW9v4s8Ph8EL1waVbUrcDATD9MAHQ5YKChe7XQbYA9N66vY1blIVCAGyO3WaQndyO7OR2nNz74M8dNfXsqKhh595ayvbWUb63lrK9tezZV8+e6noc1fVU7H9U1Tqpqm2gut69M4jLhKo6J1V1R7tTSMf9j8FN3o1nH5lGOVlGOSk4SDEqSDEcpBoOOlBJvLGPBPaRaFSRwD7aGe7/DtlwEeXcS5Rzb5PrpQJ9mmndBHbbbJRE2CmOiKDcbmeX3cYum/t5t91Opc1Gpc3Y/2zDtb+nsw4ndTippK7pBZtrRESCimPDbgVAHwiLAOjxy2Et0zQPO9TV3PHNve8xadIkHn300TZWGdy6JXZjduHsI+8JvHMN1FZAZBykD/BPcS20oGQBDWYDXRK6kB2fHehyfCIhJpKEmEh6pse3+Byny6SqroGaOifV9e5HTb2L6jondU4XdQ37H04n9Q0m9S4XDU6TeqeLBpeJ85cP08TlMnGZJk4XuEwT0zSpNiHfNNkGuFwmpgkm7meXCYbpJNJVQ5SzmmizmkhXNRGuWiJddUSYdUS4aokw67GbDdjMBiLMBuxmHXbTiQ0nNtOJ3WwgASeJpkl304mtwYVR78K2/35CAxcGLurNBuoMF3WG0/3ASb3hogEXDYaLelw0GCZOoMFw4cR0PwwTFyYuwImJy3D/Di7cGdFluJe9NvE8/+y1ceBnj1/mSrPZ12azxx7qvYOou0/CVHbH5v4vobRVWATA1NRU7Hb7Qb19paWlB/XyeWRkZDR7fEREBCkpKc2ec9999zFx4sTGnx0OB9nZ4RUwPFvCba44QgDMn+9+7jwU7MH1z6hx+DcrNHv/fMVuMxqDo4iIWFtY3BkdFRVFXl4e06dPb/L+9OnTGTmy+R0tRowYcdDxX3/9NUOHDj3k/X/R0dEkJCQ0eYSb7kndgRYEwMIl7ufs431c0dExTbMxAJ7Y6cQAVyMiIhKcwiIAAkycOJFXX32V119/nTVr1nDnnXeSn5/fuK7ffffdx9VXH7iH4Oabb2bbtm1MnDiRNWvW8Prrr/Paa6/xhz/4f8HgYJKbmAtASVUJVfWHmVzgCYCdhvqhqpbb4thCUVURUbYohmYEV20iIiLBIrjG7tpg7NixlJeX89hjj1FcXEz//v354osv6NKlCwDFxcXk5+c3Hp+bm8sXX3zBnXfeyb/+9S+ysrL4xz/+0eI1AMNVYnQiKTEplNeUs6ViC/1T+x98UI3DvZAwQKch/i3wCDy9f0MzhhIb4aft1kREREJM2ARAgFtuuYVbbrml2c/efPPNg94bPXo0S5cu9XFVoad7UnfKS8rZXLG5+QBYvBwwITHHve5bEJldOBvQ/X8iIiKHEzZDwOI9nokg63etb/6AxuHfwc1/HiDVDdUsLlkM6P4/ERGRw1EAlIP0SXZPuV+3e13zBzQGwCMsFO1ni0sWU+eqIzMus/FeRhERETmYAqAcpDEA7lpHszsFFi5zPwdZAJxTdGD3j8NudSYiImJxCoBykO5J3bEbdnbX7mZn9S/2gK0sAcd2MGyQOSgg9R1K4/IvWRr+FRERORwFQDlITERM4xDq2l1rm35YuH/STMc+EN3ez5Ud2vbK7Wx1bCXCiOD4zOBam1BERCTYKABKs3onuzegXbfrF/cBNt7/F1zLv8wtmgvAwI4DiY9q+fZoIiIiVqQAKM3q08F9H+DBPYD7A2BWcAVAz/Ivmv0rIiJyZAqA0qzGHsCfzwR2uaBo/xBwEE0AqXfVs7BkIQAjOzW/9Z+IiIgcoAAozfIEwHxH/oEt4XZthpoKsEdDer8AVtfUitIVVNVXkRyTTN/kvoEuR0REJOgpAEqzkmOSSWuXhonJht0b3G96ev8yjwV7ZOCK+wXP8i8jskZgM/RPWkRE5Ej0bSmH5FkPsPE+wCBdANqz/Iu2fxMREWkZBUA5pN4d3MPAwRwAy6vLWbNrDeDuARQREZEjUwCUQ/r5jiA01EHxj+4PgmgJGM/yL32S+5AamxrgakREREKDAqAckicAbtizgYaSH8FZCzGJkNwtwJUd4AmAGv4VERFpOQVAOaTO8Z1pF9GOWmct27Z+636zUx4EyT67LtN1IAB2UgAUERFpKQVAOSSbYWtcDmZt0SL3m0F0/9/aXWvZVbOLdhHtGNRxUKDLERERCRkKgHJYnokg6yo2ud8IogDo6f07PvN4IoNoWRoREZFgpwAoh9W4FIxzr/uNINoCzrP9m+7/ExEROToKgHJYjTOBoyIx47MgPj3AFbntrdvLitIVgAKgiIjI0VIAlMPqntQdOwa77HZ2ZgbP9m8LSxbSYDaQHZ9NdkJ2oMsREREJKQqAclgxETHk2mIAWNshK8DVHKDdP0RERFpPAVCOqHdtHQDromMCXImbaZrMLJwJwKjOowJcjYiISOhRAJTDq91LH0c5AGtdVQEuxm397vWUVJUQY4/h+IzjA12OiIhIyFEAlMMrWUnvuloA1lVuC3AxbjO3u3v/hmUOIyYiOHolRUREQokCoBxe8XL61NUDsM2xDUedI8AFwYztMwA4qfNJAa5EREQkNCkAyuEVLaeDy0WXiHgAlpcuD2g5u2p28ePOHwEFQBERkdZSAJTDK14OwJAO7vUAl+5YGsBi3Is/m5j0Se5DRlxGQGsREREJVQqAcmh1VVC2HoAhOaMBWLJjSSArarz/T71/IiIiracAKIdWshJMF7TPIC/nVAB+Kv+JmoaagJRT76pvXP9PAVBERKT1FADl0IqWu5+zBtE5vjMdYzvS4GpgZdnKgJSzbMcy9tbvJTkmmf4p/QNSg4iISDhQAJRD23//H5mDMAyDIelDgMDdB+iZ/XtipxOx2+wBqUFERCQcKADKof2sBxBgSNr+AFgamADouf9vdOfRAWlfREQkXCgASvPq9kHZOvfrzEEA5KXnAe6lYBpcDX4tZ5tjG1sdW4kwIhiZNdKvbYuIiIQbBUBp3o6f9k8ASYeETAB6JPUgPjKefQ37WLd7nV/L8fT+5WXk0T6qvV/bFhERCTcKgNI8z/Dv/t4/ALvNzqA098/+vg/Qc/+fhn9FRETaTgFQmueZALL//j+PQEwE2Vu3lyUl7vUHFQBFRETaTgFQmtdMDyAcuA9waelSTNP0SynfF3xPg9lAbmIuOQk5fmlTREQknCkAysHqq2HnWvfrzGObfNQvpR9Rtih21exiq2OrX8r5csuXAJzV9Sy/tCciIhLuFADlYCU/gemEuI6QkNXkoyh7FAM6DgD8Mwy8p2YP84rmAXB27tk+b09ERMQKFADlYD9bABrDOOhjf64H+PW2r2kwG+ib3JfcxFyftyciImIFCoBysENMAPHw3Ae4ZMcSn5cybes0AM7K1fCviIiItygAysGKV7iffzEBxOPYjsdiM2wU7i2kpKrEZ2XsqNrB4pLFAJzdVcO/IiIi3qIAKE011ELpGvfrX0wA8Wgf1Z7eHXoDvr0P8KutX2FiMjhtMJntM33WjoiIiNUoAEpTpavB1QCxyZDY+ZCHeYaB5xfP91kpjcO/mv0rIiLiVQqA0lTj8O+xzU4A8Tgl+xQAvs3/lnpnvdfLKHAUsLJsJTbDxpiuY7x+fREREStTAJSmfh4ADyMvPY+OsR1x1DmYUzTH62V8udW99t+wjGGkxqZ6/foiIiJWpgAoTbUwANptds7seiYAX2z5wutleBZ/1tp/IiIi3qcAKAc4692LQMMRAyDAObnnAPBDwQ/sq9/ntTI27N7Axj0bibRFclqX07x2XREREXFTAJQDytaDsxaiE6DDkRdd7p/an+z4bKobqvmh4AevleHp/Tux04kkRCV47boiIiLipgAoB3iGfzMGgu3I/zQMw2icoesJbW1V01DDRxs+Ag70MIqIiIh3KQDKAS28/+/nPCFtdtFsKmor2lzCpxs/pbymnIy4DA3/ioiI+IgCoBzQigDYo0MPenXoRYOrgW+2fdOm5htcDbyx6g0Aru13LZG2yDZdT0RERJqnAChuLhcU/+h+fRQBEA7M1G3rMPBXW7+icG8hHaI7cHHPi9t0LRERETk0BUBxK98I9VUQEQupPY/qVE8AXFiykJ37draqedM0ee2n1wC4su+VxEbEtuo6IiIicmQKgOLWOAFkANjsR3Vqp/adOLbjsZiYfLX1q1Y1P3P7TDbs3kC7iHaM6zOuVdcQERGRllEAFLfi5e7noxz+9fD0ArZ2UWhP799lvS8jMTqxVdcQERGRllEAFLdWTAD5uTO7nonNsLGybCVLdiw5qnOX7FjCstJlRNoiGX/M+Fa1LyIiIi2nAChgmq2eAOKRGpvKRT0uAuChOQ8d1c4gr658FYBf9fgVae3SWtW+iIiItJwCoMDurVBbAfYo6Nin1Ze5a+hdZMRlUFBZwN+W/K1F58wtmsvswtnYDBu/6febVrctIiIiLacAKAeGf9P7QURUqy8THxXPYyMfA+D9de8zv3j+YY9fXb6aO7+/E4ALe1xITkJOq9sWERGRllMAlDbf//dzI7JGMLb3WAD+NOdP7K3b2+xx+Y58fvfN79jXsI9hGcN4YNgDbW5bREREWkYBULwaAAEm5k2kc/vOFFcV8/Tipw/6vKy6jBun38iuml30Te7Lc6c8R5S99T2PIiIicnQiAl2ABJhpej0Atotsx/874f/x269+y4cbPiTaHk33pO5kxmXSsV1HHpz9IIV7C8mOz2by6ZNpH9XeK+2KiIhIyygAWp2jCPaVgWGHtH5eu+zQjKFcdcxVvLX6Ld5d++5BnyfHJPPS6S+RGpvqtTZFRESkZRQArc6zAHRaX4iM8eqlJ+ZNpEdSD9buWkvx3mIKqwop2ltEtD2aF09/keyEbK+2JyIiIi2jAGh1Rcvdz5mDvH7pCFsEF/e8uMl7pmkCYBiG19sTERGRlgmLSSC7d+9m/PjxJCYmkpiYyPjx49mzZ88hj6+vr+eee+5hwIABxMXFkZWVxdVXX01RUZH/ig4Wnh7ArEF+ac4wDIU/ERGRAAuLAHjFFVewfPlypk2bxrRp01i+fDnjxx96S7F9+/axdOlSHnroIZYuXcpHH33E+vXrueCCC/xYdRAwTSha5n7tgx5AERERCU6G6RmTC1Fr1qzhmGOOYf78+QwbNgyA+fPnM2LECNauXUvv3r1bdJ1FixZx/PHHs23bNnJyWrYgscPhIDExkYqKChISElr9OwRMRSH87Rj3BJD7CyEyNtAViYiI+FzIf397Qcj3AM6bN4/ExMTG8AcwfPhwEhMTmTt3bouvU1FRgWEYJCUlHfKY2tpaHA5Hk0dI8wz/duyj8CciImIhIR8AS0pKSEtLO+j9tLQ0SkpKWnSNmpoa7r33Xq644orD/j+BSZMmNd5nmJiYSHZ2iM9i9UwA8dP9fyIiIhIcgjYAPvLII40TBg71WLx4MdD8jFLTNFs02aC+vp5x48bhcrmYPHnyYY+97777qKioaHwUFBS07pcLFo0TQAYHtAwRERHxr6BdBmbChAmMGzfusMd07dqVH3/8kR07dhz02c6dO0lPTz/s+fX19Vx22WVs2bKF77777oj3AURHRxMdHX3k4kOBafp0CRgREREJXkEbAFNTU0lNPfIuESNGjKCiooKFCxdy/PHHA7BgwQIqKioYOXLkIc/zhL8NGzbw/fffk5KS4rXaQ4KjCKpK3RNAMvoHuhoRERHxo6AdAm6pvn37ctZZZ3HDDTcwf/585s+fzw033MB5553XZAZwnz59+PjjjwFoaGjg17/+NYsXL+add97B6XRSUlJCSUkJdXV1gfpV/EsTQERERCwr5AMgwDvvvMOAAQMYM2YMY8aMYeDAgbz11ltNjlm3bh0VFRUAbN++nc8++4zt27czaNAgMjMzGx9HM3M4pGkCiIiIiGUF7RDw0UhOTubtt98+7DE/X+6wa9euhPjyh23n6QHU/X8iIiKWExY9gHKUfj4BRD2AIiIilqMAaEWNE0BskK4JICIiIlajAGhFjRNA+kJUu4CWIiIiIv6nAGhFGv4VERGxNAVAK9IEEBEREUtTALQaTQARERGxPAVAq6ks1gQQERERi1MAtJqiZe7njn00AURERMSiFACtxjP8q/v/RERELEsB0Go8E0B0/5+IiIhlKQBaiWlC4RL36055ga1FREREAkYB0Ep2b4V95WCL1AQQERERC1MAtBJP71/GAIiMCWwtIiIiEjAKgFbiCYCdhwa2DhEREQkoBUAr2b7Y/dxJAVBERMTKFACtoqEOile4X2sCiIiIiKUpAFpF6Spw1kJMEqR0D3Q1IiIiEkAKgFbROPybB4YR2FpEREQkoBQArULr/4mIiMh+CoBWoRnAIiIisp8CoBVU74Gy9e7X6gEUERGxPAVAKyha6n5O6gJxqYGtRURERAJOAdAKNPwrIiIiP6MAaAXbPRNAFABFREREATD8mSYU/mwJGBEREbE8BcBwV1EAVTvBFgGZAwNdjYiIiAQBBcBw51kAOr0/RMYGthYREREJCgqA4U4LQIuIiMgvKACGO80AFhERkV9QAAxnznooWu5+rRnAIiIisp8CYDgrXQ0N1RCdCCk9Al2NiIiIBAkFwHCWv8D93DkPbPqPWkRERNyUCsLZtjnu5y4jA1uHiIiIBBUFwHBlmrBtrvt1lxMCW4uIiIgEFQXAcLVrM1SVgj0KsoYEuhoREREJIgqA4coz/NtpKETGBLYWERERCSoKgOGqcfhX9/+JiIhIUwqA4UoTQEREROQQFADDUcV22JMPhh2yjw90NSIiIhJkFADD0bZ57ufMgRAdH9haREREJOgoAIajxuFfLf8iIiIiB1MADEeaACIiIiKHoQAYbqrKoGyd+3XOiMDWIiIiIkFJATDc5O+//y/tGGiXHNhaREREJCgpAIYbz/Cvev9ERETkEBQAw43W/xMREZEjUAAMJzUOKFnpfq0AKCIiIoegABhOChaC6YIOuZCQFehqREREJEgpAIYTrf8nIiIiLaAAGE4aA6AmgIiIiMihKQCGi+rdsH2x+3XuSYGtRURERIKaAmC42PgtmE7o2BeScgJdjYiIiAQxBcBwseFr93OvMYGtQ0RERIKeAmA4cDlhw3T3655nBrYWERERCXoKgOFg+2Ko3gUxiZA9LNDViIiISJBTAAwHG75yP3c/DewRga1FREREgp4CYDhY77n/T8O/IiIicmQKgKGuohB2rAQM6HF6oKsRERGREKAAGOo8s387D4W41MDWIiIiIiFBATDUeQKgZv+KiIhICykAhrL6Gtj8g/u11v8TERGRFlIADGXbZkP9PmifARkDA12NiIiIhAgFwFDmmf3b8wwwjMDWIiIiIiFDATBUmeaB9f+0/IuIiIgcBQXAUFW2AXZvBVskdDs50NWIiIhICAmLALh7927Gjx9PYmIiiYmJjB8/nj179rT4/JtuugnDMHjuued8VqPXrf7E/dz1RIiOD2gpIiIiElrCIgBeccUVLF++nGnTpjFt2jSWL1/O+PHjW3TuJ598woIFC8jKyvJxlV7kcsLSf7tfH3t5YGsRERGRkBPyG8euWbOGadOmMX/+fIYNGwbAK6+8wogRI1i3bh29e/c+5LmFhYVMmDCBr776inPPPddfJbfdpu+gogBikuCYXwW6GhEREQkxId8DOG/ePBITExvDH8Dw4cNJTExk7ty5hzzP5XIxfvx47r77bvr16+ePUr1nyZvu50FXQGRMQEsRERGR0BPyPYAlJSWkpaUd9H5aWholJSWHPO+JJ54gIiKC2267rcVt1dbWUltb2/izw+E4umK9wVEM6750vx5yjf/bFxERkZAXtD2AjzzyCIZhHPaxePFiAIxm1sAzTbPZ9wGWLFnC3//+d958881DHtOcSZMmNU40SUxMJDs7u3W/XFssextMJ+SMgLQ+/m9fREREQl7Q9gBOmDCBcePGHfaYrl278uOPP7Jjx46DPtu5cyfp6enNnjdr1ixKS0vJyclpfM/pdHLXXXfx3HPPsXXr1mbPu++++5g4cWLjzw6Hw78h8OeTP/Ku9V+7IiIiElaCNgCmpqaSmpp6xONGjBhBRUUFCxcu5PjjjwdgwYIFVFRUMHLkyGbPGT9+PKeffnqT984880zGjx/Pb37zm0O2FR0dTXR09FH8Fl626XuoyNfkDxEREWmToA2ALdW3b1/OOussbrjhBl566SUAbrzxRs4777wmM4D79OnDpEmTuOiii0hJSSElJaXJdSIjI8nIyDjsrOGAW/KG+/nYyyEyNrC1iIiISMgK2nsAj8Y777zDgAEDGDNmDGPGjGHgwIG89dZbTY5Zt24dFRUVAarQCypLDkz+yNPkDxEREWm9kO8BBEhOTubtt98+7DGmaR7280Pd9xc0lr3lnvyRPRzS+ga6GhEREQlhYdEDGPZqK2Hx/uFfTf4QERGRNlIADAXT7gVHISTmQL8LA12NiIiIhDgFwGC35nP32n8YcNGLmvwhIiIibaYAGMwqd8DU/TuVnHAbdD0hsPWIiIhIWFAADFamCZ9NgH3lkD4ATnkg0BWJiIhImFAADFZL3oANX4M9Gi5+GSICuAC1iIiIhBUFwGBUthG+2t/jd/rDkH5MYOsRERGRsKIAGIy+fQTq90HuSTDsd4GuRkRERMJMWCwEHXYueB7apcBJd4NNGV1ERES8SwEwGMUmwfl/D3QVIiIiEqbUvSQiIiJiMQqAIiIiIhajACgiIiJiMQqAIiIiIhajACgiIiJiMQqAIiIiIhajACgiIiJiMQqAIiIiIhajACgiIiJiMQqAIiIiIhajACgiIiJiMQqAIiIiIhajACgiIiJiMRGBLiCUmaYJgMPhCHAlIiIi0lKe723P97gVKQC2QWVlJQDZ2dkBrkRERESOVmVlJYmJiYEuIyAM08rxt41cLhdFRUXEx8djGEarr+NwOMjOzqagoICEhAQvVijN0d/bv/T39i/9vf1Lf2//8tbf2zRNKisrycrKwmaz5t1w6gFsA5vNRufOnb12vYSEBP0PiB/p7+1f+nv7l/7e/qW/t3954+9t1Z4/D2vGXhERERELUwAUERERsRgFwCAQHR3Nww8/THR0dKBLsQT9vf1Lf2//0t/bv/T39i/9vb1Hk0BERERELEY9gCIiIiIWowAoIiIiYjEKgCIiIiIWowAoIiIiYjEKgAE2efJkcnNziYmJIS8vj1mzZgW6pLA1adIkjjvuOOLj40lLS+PCCy9k3bp1gS7LEiZNmoRhGNxxxx2BLiWsFRYWctVVV5GSkkK7du0YNGgQS5YsCXRZYamhoYEHH3yQ3NxcYmNj6datG4899hgulyvQpYWFmTNncv7555OVlYVhGHzyySdNPjdNk0ceeYSsrCxiY2M5+eSTWbVqVWCKDVEKgAE0ZcoU7rjjDh544AGWLVvGqFGjOPvss8nPzw90aWFpxowZ3HrrrcyfP5/p06fT0NDAmDFjqKqqCnRpYW3RokW8/PLLDBw4MNClhLXdu3dzwgknEBkZyZdffsnq1at55plnSEpKCnRpYemJJ57gxRdf5Pnnn2fNmjU8+eSTPPXUU/zzn/8MdGlhoaqqimOPPZbnn3++2c+ffPJJnn32WZ5//nkWLVpERkYGZ5xxBpWVlX6uNHRpGZgAGjZsGEOGDOGFF15ofK9v375ceOGFTJo0KYCVWcPOnTtJS0tjxowZnHTSSYEuJyzt3buXIUOGMHnyZP785z8zaNAgnnvuuUCXFZbuvfde5syZo1EEPznvvPNIT0/ntddea3zvkksuoV27drz11lsBrCz8GIbBxx9/zIUXXgi4e/+ysrK44447uOeeewCora0lPT2dJ554gptuuimA1YYO9QAGSF1dHUuWLGHMmDFN3h8zZgxz584NUFXWUlFRAUBycnKAKwlft956K+eeey6nn356oEsJe5999hlDhw7l0ksvJS0tjcGDB/PKK68EuqywdeKJJ/Ltt9+yfv16AFasWMHs2bM555xzAlxZ+NuyZQslJSVNvj+jo6MZPXq0vj+PQkSgC7CqsrIynE4n6enpTd5PT0+npKQkQFVZh2maTJw4kRNPPJH+/fsHupyw9P7777N06VIWLVoU6FIsYfPmzbzwwgtMnDiR+++/n4ULF3LbbbcRHR3N1VdfHejyws4999xDRUUFffr0wW6343Q6efzxx7n88ssDXVrY83xHNvf9uW3btkCUFJIUAAPMMIwmP5umedB74n0TJkzgxx9/ZPbs2YEuJSwVFBRw++238/XXXxMTExPocizB5XIxdOhQ/vKXvwAwePBgVq1axQsvvKAA6ANTpkzh7bff5t1336Vfv34sX76cO+64g6ysLK655ppAl2cJ+v5sGwXAAElNTcVutx/U21daWnrQ/6sR7/r973/PZ599xsyZM+ncuXOgywlLS5YsobS0lLy8vMb3nE4nM2fO5Pnnn6e2tha73R7ACsNPZmYmxxxzTJP3+vbty4cffhigisLb3Xffzb333su4ceMAGDBgANu2bWPSpEkKgD6WkZEBuHsCMzMzG9/X9+fR0T2AARIVFUVeXh7Tp09v8v706dMZOXJkgKoKb6ZpMmHCBD766CO+++47cnNzA11S2DrttNNYuXIly5cvb3wMHTqUK6+8kuXLlyv8+cAJJ5xw0LJG69evp0uXLgGqKLzt27cPm63pV6jdbtcyMH6Qm5tLRkZGk+/Puro6ZsyYoe/Po6AewACaOHEi48ePZ+jQoYwYMYKXX36Z/Px8br755kCXFpZuvfVW3n33XT799FPi4+Mbe18TExOJjY0NcHXhJT4+/qB7K+Pi4khJSdE9lz5y5513MnLkSP7yl79w2WWXsXDhQl5++WVefvnlQJcWls4//3wef/xxcnJy6NevH8uWLePZZ5/lt7/9baBLCwt79+5l48aNjT9v2bKF5cuXk5ycTE5ODnfccQd/+ctf6NmzJz179uQvf/kL7dq144orrghg1SHGlID617/+ZXbp0sWMiooyhwwZYs6YMSPQJYUtoNnHG2+8EejSLGH06NHm7bffHugywtrUqVPN/v37m9HR0WafPn3Ml19+OdAlhS2Hw2HefvvtZk5OjhkTE2N269bNfOCBB8za2tpAlxYWvv/++2b/9/qaa64xTdM0XS6X+fDDD5sZGRlmdHS0edJJJ5krV64MbNEhRusAioiIiFiM7gEUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsZj/D+KX+P2Tg+sHAAAAAElFTkSuQmCC",
      "text/html": [
       "\n",
       "            <div style=\"display: inline-block;\">\n",
       "                <div class=\"jupyter-widgets widget-label\" style=\"text-align: center;\">\n",
       "                    Figure\n",
       "                </div>\n",
       "                <img src='data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABV9UlEQVR4nO3dd3xUVf7/8dedSSWkkIQ0SCB0pAgEpShixe5aVsGCumtdZS24rn0tv3VZ67q7LvbyXSvuWnEVxUbvTaT3hBRCAmRCSJ25vz+GCUYChGT6fT8fj3nMZObeez7JsjvvPeeecwzTNE1ERERExDJsgS5ARERERPxLAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYhQARURERCxGAVBERETEYiICXUAoc7lcFBUVER8fj2EYgS5HREREWsA0TSorK8nKysJms2ZfmAJgGxQVFZGdnR3oMkRERKQVCgoK6Ny5c6DLCAgFwDaIj48H3P+AEhISAlyNiIiItITD4SA7O7vxe9yKFADbwDPsm5CQoAAoIiISYqx8+5Y1B75FRERELEwBUERERMRiFABFRERELEYBUERERMRiFABFRERELEYBUERERMRiFABFRERELEYBUERERMRiFABFRERELEYBUERERMRiFABFRERELEYBUERERMRiFACD0OoiBxdNnsO3a3ZgmmagyxEREZEwowAYhF6csYll+Xu47v8Wc84/ZvPFymJcLgVBERER8Q4FwCD00HnHcNPobsRF2VlT7OCWd5Yy5rmZTF1RFOjSREREJAwYpsYYW83hcJCYmEhFRQUJCQlev/7uqjremLuVN+ZsobKmAYBrR3blofOOwW4zvN6eiIiIFfj6+zsUqAcwSG2p2IJh38fEM3ox595TmXBKDwDenLuVG/69mL21DQGuUEREREKVAmAQenLRk1zwyQX8Z/1/AEiIieQPZ/Zm8pVDiI6w8d3aUi59cR7FFdUBrlRERERCkQJgEOqb3BeAKeum0OA60NN3zoBMptw0gtT20awpdvCr5+ewqqgiUGWKiIhIiFIADEJndj2T5JhkduzbwfcF3zf5bFB2Ep/cOpJe6e0prazlmtcXUbhHPYEiIiLScgqAQSjKHsUlPS8B4N017x70eecO7fjv70bSJyOesr21XPfmIt0TKCIiIi2mABikLut9GXbDzuIdi1m/e/1BnyfERPLatceR2j6atSWV3PH+MpxaK1BERERaQAEwSGXEZXBqzqkAvLf2vWaP6ZQUyytX5xEVYeObNaX89cs1/ixRREREQpQCYBC7os8VAPxv8/+oqG1+ssfgnA48c+mxALwyawvvLcz3W30iIiISmhQAg1heeh49O/SkuqGaTzZ+csjjzj82iztP7wXAQ5/8xPKCPf4pUEREREJSWAXAyZMnk5ubS0xMDHl5ecyaNatF582ZM4eIiAgGDRrk2wKPkmEYjb2A7699H5fpOuSxt53Wg3MHZNLgMrnrg+XU1Dv9VaaIiIiEmLAJgFOmTOGOO+7ggQceYNmyZYwaNYqzzz6b/PzDD4lWVFRw9dVXc9ppp/mp0qNzbrdziY+KZ/ve7cwunH3I4wzD4PGL+tMxPppNO6t45ut1fqxSREREQknYBMBnn32W6667juuvv56+ffvy3HPPkZ2dzQsvvHDY82666SauuOIKRowY4adKj05sRCwX97gYgHfXHrwkzM8ltYvirxcPAODV2VtYtHWXz+sTERGR0BMWAbCuro4lS5YwZsyYJu+PGTOGuXPnHvK8N954g02bNvHwww+3qJ3a2locDkeThz+M7TMWA4M5hXMo3lt82GNP65vOr/M6Y5rwh/+sYF+d1gcUERGRpsIiAJaVleF0OklPT2/yfnp6OiUlJc2es2HDBu69917eeecdIiIiWtTOpEmTSExMbHxkZ2e3ufaWyI7PZkCqu2dvYcnCIx7/p/OPITMxhm3l+3jiy7W+Lk9ERERCTFgEQA/DMJr8bJrmQe8BOJ1OrrjiCh599FF69erV4uvfd999VFRUND4KCgraXHNL5WXkAbB4x+IjHpsQE8kTlwwE4P/mbWPuxjKf1iYiIiKhJSwCYGpqKna7/aDevtLS0oN6BQEqKytZvHgxEyZMICIigoiICB577DFWrFhBREQE3333XbPtREdHk5CQ0OThL8elHwfA4pIjB0CAk3p15IphOQDc9/FKahs0K1hERETcwiIARkVFkZeXx/Tp05u8P336dEaOHHnQ8QkJCaxcuZLly5c3Pm6++WZ69+7N8uXLGTZsmL9Kb7HBaYOxGTa2791OSVXzw9q/dP85fUmLj2Zb+T7emLPVtwWKiIhIyAiLAAgwceJEXn31VV5//XXWrFnDnXfeSX5+PjfffDPgHr69+uqrAbDZbPTv37/JIy0tjZiYGPr3709cXFwgf5VmtY9qT9/kvkDLhoEB2kdHcM9ZfQD457cbKK2s8Vl9IiIiEjrCJgCOHTuW5557jscee4xBgwYxc+ZMvvjiC7p06QJAcXHxEdcEDHZD04cCLR8GBrhocCeOzU6iqs7JU9O0NqCIiIiAYZqmGegiQpXD4SAxMZGKigq/3A/4Q8EP/P6739M1oStTL5ra4vOW5u/m4snu5XA+vfUEjs1O8k2BIiIiIcDf39/BKGx6AK1gSPoQDAy2Orayc9/Olp+X04GLB3cC4NGpq1DmFxERsTYFwBCSEJVAn2T3PX1Ldiw5qnP/eFYf2kXZWZq/h89WFPmiPBEREQkRCoAhJi/dvR7gopJFR3VeRmIMt57SA4BJX6zVDiEiIiIWpgAYYoZm7J8I0sKZwD933Ym5dO4QS4mjhtdnb/F2aSIiIhIiFABDTF6auwdwc8VmyqvLj+rcmEg7d5/ZG4CXZm6mYl+91+sTERGR4KcAGGKSYpLo2aEncPT3AQKcPzCL3unxVNY08Mqszd4uT0REREKAAmAIalwPsBXDwDabwZ1nuPc/fn3OFsr31nq1NhEREQl+CoAh6LiM/fsCtyIAApzZL50BnRLZV+fkxRmbvFmaiIiIhAAFwBDkmQm8YfcG9tTsOerzDcPgrjHuXsB/z9vGDoe2iBMREbESBcAQlByTTPfE7gAsKT36+wABRvfqyHFdO1Db4OL57zZ6szwREREJcgqAIapxOZij2Bf459y9gO4Zwe8vyqdg1z6v1SYiIiLBTQEwRA1JGwLAyrKVrb7G8G4pjOqZSr3T5O/fbvBWaSIiIhLkFABDlGdLuPW71+MyXa2+zsT9M4I/XlZIfrl6AUVERKxAATBEdUnoQow9huqGagoqC1p9ncE5HRjVMxWny+SlmZoRLCIiYgUKgCHKbrM3Lgi9dtfaNl1rwv49gv+zeLtmBIuIiFiAAmAI69XBPXy7bte6Nl1nWLcUjuvagTqni1dmancQERGRcKcAGMI89wGu2922AAhw6/5ewHcW5LOrqq7N1xMREZHgpQAYwjwBsK1DwOBeF7B/pwSq6528OWdLm68nIiIiwUsBMIR57gEs3VfK7prdbbqWYRjcerK7F/DNuVuprKlvc30iIiISnBQAQ1hcZBw58TmAd4aBz+yXQfeOcThqGnhr/rY2X09ERESCkwJgiOud7N7No60TQQBsNoNb9vcCvjZrC9V1zjZfU0RERIKPAmCI8+Z9gAAXDMqic4dYyqvqmLIo3yvXFBERkeCiABjienfY3wPohSFggEi7jZtO6gbAq7O30OBs/S4jIiIiEpwUAEOcZwh4y54t1DprvXLNX+dlkxwXxfbd1Xz5U4lXrikiIiLBQwEwxKW3SycpOokGs4FNe7yzlVtslJ2rR3QB4OWZmzFN0yvXFRERkeCgABjiDMM4MAzshYkgHleP6EpMpI2VhRXM21zuteuKiIhI4CkAhoHGmcBeug8QIDkuikvzsgF3L6CIiIiEDwXAMODtmcAe14/KxWbAD+t2sq6k0qvXFhERkcBRAAwDvTr0AmD9rvVevV+vS0ocZ/XPANQLKCIiEk4UAMNAt8RuRNoiqayvpKiqyKvXvvGk7gB8tqKQkooar15bREREAkMBMAxE2iPpkeTewcPbw8CDspM4PjeZeqfJG3O2ePXaIiIiEhgKgGHi58PA3uZZGPrdBflU1tR7/foiIiLiXwqAYcJXE0EATumdRo+09lTWNvDB4u1ev76IiIj4lwJgmPDFUjAeNpvBdSfmAvC6tocTEREJeQqAYcIzBFy4txBHncPr179ocCeS46Io3FPNV6t2eP36IiIi4j8KgGEiMTqRrLgswDf3AcZE2rlquHt7uFdna0kYERGRUKYAGEa6J7mXbNlc4ZuANn54F6LsNpbl72HJtt0+aUNERER8TwEwjOQmuu/T21Lhm+VaOsZHc+Fgdy/jq7PUCygiIhKqFADDSLdE93ItvuoBBLh+lLuNr1aVkF++z2ftiIiIiO8oAIaRbkm+D4C90uM5qVdHXCa8MVcLQ4uIiIQiBcAw4ukBLKkqYV+973rnrt+/JMwHiwqoqNbC0CIiIqFGATCMJEYnkhyTDMAWh+9650b1TKVXenuq6py8vzDfZ+2IiIiIbygAhpnG+wD3+G4Y2DAMrj/R3c7/zd1KvRaGFhERCSkKgGHGEwB9NRPY44JBWaS2j6KoooYvfyrxaVsiIiLiXQqAYcazFIwvJ4KAe2Ho8cO7Au4lYUzT9Gl7IiIi4j0KgGHGH0vBeFw1PIeoCBs/bq9gsRaGFhERCRkKgGHGsxRMgaOAepdvZ+imtI/mkiGdAC0MLSIiEkoUAMNMert02kW0o8FsoMBR4PP2fnuCe8j569U72FZe5fP2REREpO0UAMOMYRg+3xLu53qmx3Ny746YJrwxZ6vP2xMREZG2UwAMQ/68DxBoXBLmg8UFVOzTwtAiIiLBTgEwDPljS7ifO6FHCn0y4tlX5+S9RVoYWkREJNgpAIah3AT/LAXjYRgG1+3fHu7NOVupa9DC0CIiIsFMATAM5SYduAfQZfonjLkXho6mxFHD/1YW+aVNERERaR0FwDCUHZ9NhBFBdUM1O6p2+KXN6Ag7147sAsArM7doYWgREZEgpgAYhiJtkeQk5AD+GwYGuHJYF2Ij7awudjB3U7nf2hUREZGjowAYpvw9ExigQ1wUlw3tDMDLM7UwtIiISLBSAAxT/lwL8Od+e2IuNgNmrN/JupJKv7YtIiIiLaMAGKY8AdCfPYAAXVLiOKt/BgCvaHs4ERGRoKQAGKY8awH6uwcQ4IZR7rY/XV7IDkeN39sXERGRw1MADFOetQB31exiT80ev7Y9OKcDx3XtQL3T5M25W/3atoiIiByZAmCYahfZjsy4TMD/w8AA1+/vBXxn/jaqahv83r6IiIgcmgJgGAvETGCP0/umk5sah6OmgSmLCvzevoiIiByaAmAYC9REEAC77cD2cK/N3kK9U9vDiYiIBAsFwDDmmQgSiAAI8Ou8zqS2j6JwTzVTV2h7OBERkWChABjGPBNBtlZsDUj7MZF2fnOCu4YXZ2zC5dL2cCIiIsEgrALg5MmTyc3NJSYmhry8PGbNmnXIYz/66CPOOOMMOnbsSEJCAiNGjOCrr77yY7W+1zWxKwBFe4uoc9YFpIarhnehfXQE63fs5ft1pQGpQURERJoKmwA4ZcoU7rjjDh544AGWLVvGqFGjOPvss8nPz2/2+JkzZ3LGGWfwxRdfsGTJEk455RTOP/98li1b5ufKfSclJoX2ke0xMSmoDMxEjMTYSK4c7t6X+IUfNgWkBhEREWnKME0zLMblhg0bxpAhQ3jhhRca3+vbty8XXnghkyZNatE1+vXrx9ixY/nTn/7UouMdDgeJiYlUVFSQkJDQqrp9bdzn41hVvornTnmO03JOC0gNpY4aTnzie+qcLv5z8wiO65ockDpEREQgNL6/fS0segDr6upYsmQJY8aMafL+mDFjmDt3bouu4XK5qKysJDk5vMJJl4QuQODuAwRIS4jhkrxOgHoBRUREgkFEoAvwhrKyMpxOJ+np6U3eT09Pp6SkpEXXeOaZZ6iqquKyyy475DG1tbXU1tY2/uxwOFpXsB91TegKwDbHtoDWceNJ3Xl/UQHfrS1lbYmDPhlt/H9cNQ7YOhucdcDPOrFTekJG/7ZdW0REJMyFRQD0MAyjyc+maR70XnPee+89HnnkET799FPS0tIOedykSZN49NFH21ynP3kmggQ6AOamxnFO/0z+t7KYl2Zs5m9jB7XuQg21sOg1mPU07Ctv/piuo2DErdDzTLCFRSe3iIiIV4XFt2Nqaip2u/2g3r7S0tKDegV/acqUKVx33XV88MEHnH766Yc99r777qOioqLxUVAQ/DtcNA4BO7YGthDg5tHdAfhsRREFu/Yd3ckuJyx/D/45FL66zx3+knIgZ6T70eUEyB4GtgjYOgveGwf/Os4dFhtqj3x9ERERCwmLHsCoqCjy8vKYPn06F110UeP706dP51e/+tUhz3vvvff47W9/y3vvvce55557xHaio6OJjo72Ss3+4gmAu2p24ahzkBAVuJtdB3ROZFTPVGZtKGPyD5uYdPGAlp24dye8fTGU/Oj+OT4TTr4XBl0F9l/8E64ohIUvweI3oXwj/G8iLH8Xxr4NCZle/X1ERERCVVj0AAJMnDiRV199lddff501a9Zw5513kp+fz8033wy4e++uvvrqxuPfe+89rr76ap555hmGDx9OSUkJJSUlVFRUBOpX8Im4yDg6xnYEYFtFYIeBAW47rScA/11SQOGe6iOfUFcF717mDn8xiXD6I/D7pZB37cHhDyCxE5zxGExcDWc9ATFJULgYXh4NBQu9+auIiIiErLAJgGPHjuW5557jscceY9CgQcycOZMvvviCLl3cPWDFxcVN1gR86aWXaGho4NZbbyUzM7PxcfvttwfqV/AZz32AwTAMfFzXZEZ0S6HeafLikWYEu5zw4fVQtBRiO8D138KJd0JUuyM3FN0eht8MN34PacfA3h3wxjmw5P+884uIiIiEsLBZBzAQQmUdoUfnPcp/1/+XmwbexITBEwJdDvM2lXP5K/OJstuY+cdTyEiMOfgg04Qv/gCLXgV7NFzzGeQMb12DtXvhk5thzVT3z8fdAGc/qQkiIiIWFSrf376kb0ALCJalYDyGd0vm+K7J1DldvDjjEL2Ac//pDn8YcPHLrQ9/4O4NvPTfcMqD7p8XvQJf3u0OmSIiIhakAGgBnokgwRIADcNovBfwvYX5lDpqmh6w6hOY/pD79ZmPQ78L296ozQaj74aLXwEMd7j8/vG2X1dERCQEKQBawM+XggmWEf8TeqQwJCeJ2gYXL8/cfOCDmgr3zF2AYTfD8Fu82/DAy+Dcp92vZz4Fc5/37vVFRERCgAKgBXSO74zdsFPdUE3pvtJAlwM07QV8e8E2yvbuX6tv9t/ca/yl9IQxf4YWLOR91I67Hk7d38P49QOw9C3vtyEiIhLEFAAtINIWSef4zkDwDAMDjO7VkWM7J1JT7+KVmZthTwHMm+z+8IzHwB7pu8ZH3QUj9k+ImXobrPncd22JiIgEGQVAiwimHUE8ft4L+O9526j56mFw1kKXE6H32b5u3N3DOHg8mC746EbYsdq3bYqIiAQJBUCLCLaJIB6n9kljcE4SPRvWE7PmQ/ebZ/po6PeXDAPOew5yT4L6Knj/Cqje7ft2RUREAkwB0CI8S8EEUw8guHsB7x7Tiwci3wGgqvclkDXYfwXYI+DXb0JiDuzeAh/e4F6AWkREJIwpAFpEsK0F+HMj6xcwzLaWGjOSZ1xj/V9AXAqMfQsiYmDjdPhhkv9rEBER8aNmNlOVcOQZAt5euZ16Vz2RNh9OsDgaznqY/icAXnOezRs/NTBuRyW90uP9W0fWIDj/H/Dxje7lYTKPhb7nt/h00zSZVzSPLY4tOOocVNZV4qh14DSdjOo0ilNyTiE2ItZ39YuIiBwFBUCLSGuXRmxELNUN1RRWFjbuDxxwP30EuzZBu1Q2pN+AuaaSZ75ex0vjh/q/lmPHQvFymD8ZPr4ZOvaB1J5HPG1V+SqeXPgkS0uXNvv555s/Jy4yjjFdxnB+9/PJS8/DZqjzXUREAkcB0CIMw6BLQhfW7lrLVsfW4AmAi193Pw+7iQl9B/PZ2pl8tWoHywv2MCg7yf/1nPEYFP8I22bDf38D130Dkc3sVQyUVZfxz2X/5OMNH2NiEmOPYVTnUSRGJ5IQlUB8VDx76/Yybes0CvcW8vHGj/l448f0SOrBpFGT6JPcx8+/nIiIiJsCoIV0TejK2l1rg+c+wB2roWA+GHYYPJ4eCfFcPKQz/12ynae/Wsfb1w/zf032SLjkVXjxBChZ6d6S7pynDjrs882f8+f5f6aqvgqAc3LP4c68O8mIyzjo2NuG3May0mVM3TSVr7Z+xcY9G7n8f5dz2+DbuKbfNeoNFBERv9M3j4UE3VqAS95wP/c+GxIyAbj9tJ5E2g1mbyxj9oaywNSVkAkXveR+vfBlWDO1ycezts/iwdkPUlVfRb+Ufrx19ls8cdITzYY/AJthIy89j0dGPsL/Lv4fp2afSoOrgWeXPMv1X19PSVWJr38jERGRJhQALSSo1gKs2wcrprhfD/1N49vZye24ari7zj//bzVOV4D2Lu55Boz8vfv1p7fCnnwA1pSv4a4Zd+E0nVzQ/QLePfddBqUNavFlk2OSee6U53h05KPERsSyqGQRF392MbO2z/LBLyEiItI8BUALaVwKpiIIAuCqj6C2ApK6QLdTm3x0+2k9SYyNZG1JJVMWFQSoQODUP0GnPKipgA+vp9iRz63f3kp1QzXDMofxyIhHWjV8axgGF/e8mP+c/x8GpA6gsq6S276/ja+3fu2DX0JERORgCoAW0iXR3bNWWl3aeO9awHgmf+RdC7am/wyT2kVxx+nu2bfPfL0OR029n4vbLyIKfv06RCdQuX0ht0y9nJ3VO+mR1INnT36WyDbuVdwloQv/d/b/cXbXs2lwNXD3zLv5dOOnXipeRETk0BQALSQhKoHkmGQgwMPAxT9C4RKwRcDgq5o95KrhXejWMY7yqjr+9d1GPxf4Mx26Un/e35iYnsrGBgcdoxKZfNpkEqISvHL5SFskk0ZN4uKeF+MyXTw450HeX/u+V64tIiJyKAqAFtO4JVzF1sAV4Zn80fd8aJ/W7CGRdhsPnXsMAG/M2cq28sD1WL5rq2J+bCyxLhf/2lFOpi3aq9e32+w8POJhrux7JQCPL3ic13963attiIiI/JwCoMUEfCJIbSX8+IH7dd5vDnvoyb07MqpnKnVOF5O+WOuH4g62q2YXL61wzwi+pzaSvnuK4LPfg+ndySk2w8Y9x93DDQNuAOBvS/6mnkAREfEZBUCL8SwAHbClYFb+F+r2QkoPyD3psIcahsFD5x2DzYBpq0qYt6ncT0UeMHn5ZCrrK+mT3IcLz38dbJGw9vMDvZheZBgGtw25jd8d+zsA/rLgL5oYIiIiPqEAaDEB7wFc+m/3c961YBhHPLxXejxXDnPX/OjUVTQ4XT4srqkNuzfwn/X/AeCPx/0Re6chcPoj7g+n3Q+lvumV/N2xv+OyXpdhYnLvrHtZVLLIJ+2IiIh1KQBaTOM9gI6tmF4exjyiPflQtBQwYODYFp925xm9GpeFeWPOVp+V93OmafLkoidxmS7O6HIGx2Uc5/5g+C3Q/VRoqIYPr4P6Gq+3bRgG9w+7n9NzTqfeVc9t393Gul3rvN6OiIhYlwKgxWTHZ2MzbFTVV1Fe4+ch1bX/cz/njDjk5I/mJMdFcf857n1zn52+nu279/miuiZmbJ/B/OL5RNoiuTPvzgMf2Gxw4YvQLhV2/ATfPOyT9u02O3896a/kpeext34vN39zM9srt/ukLRERsR4FQIuJskeRFZcFwJaKLf5t3LOlWt/zjvrUy4Zmc3xuMtX1Tv706Sqf9l7WO+t5evHTAIw/ZjzZ8dlND4hPhwtfcL9e8CKs+9IndUTbo/nHqf+gZ4eelFWXceu3t+Koc/ikLRERsRYFQAvyTATx632AVWWQP8/9us/RB0DDMPjLRf2JtBt8t7aUL3/y3f657619j22ObaTEpDTOyj1IrzEw/Fb3609uAUeRT2pJiErgxdNfJK1dGpsrNvOHH/5Ag6vBJ22JiIh1KABaUEDWAlz3BZguyBgIHbq06hI90uL53ejuADzy2Sqf7BBS66zl1ZWvAvD7wb+nfVT7Qx98+sOQeSxU74IPbwCX0+v1AKS1S+P5U58nNiKWecXz+OvCv/r//k0REQkrCoAW1LgnsD97ANd87n7ue36bLnPLKT3ITY2jtLKWp7/y/sSIL7d8ye7a3WTEZfCrHr86/MER0fDrNyCqPWybDbOe8Xo9Hn1T+vLXUX/FwGDKuim8u/Zdn7UlIiLhTwHQgjx7AvttLcAaB2z+3v26jQEwJtLO4xf2B+Ct+dtYsm13W6trZJom765xB6uxvccSYYs48kkp3eHc/cHvh0mwbZ7X6vmlU3NOZWLeRACeXPQkM7fP9FlbIiIS3hQALcjTA7i9cjv1Lu8Pox5k43Rw1kFyd+jYp82XG9kjlUuGdMY04Q//WcG+Ou/cE7di5wrW7FpDlC2KS3pe0vITjx0HA8e5h7g/vA6qfDe7+pp+1zTuG3z3jLtZv3u9z9oSEZHwpQBoQWnt0oiNiKXBbKBor28mLzTROPx7XosWf26JP513DJmJMWwpq+Lx/63xyjXfWfMOAOd2O5cOMR2O7uRzn3bvbuIohI+u99n9gIZh8OCwBzku4zj2Nezjtu9uo7za/zukiIhIaFMAtCCbYSMnPgfww0SQ+hrYsH87sz5tG/79ucR2kTx96bEAvLMgn+/W7mjT9Ur3lfLNtm8AuKLvFUd/geh4uOzfEBELm76DmU+1qZ7DibRH8uzoZ8mOz6ZwbyF3/nAndc46n7UnIiLhRwHQovy2J/CWGe69f+MzoVOeVy99Qo9UfntCLgB//O9KyvfWtvpaH6z7gAazgSFpQ+iT3Mph6vR+cP5z7tc//BU2ftvqeo4kKSaJ5097nvjIeJaVLuOxeY9pZrCIiLSYAqBFefYE9nkA9Cz+3Odc9y4aXvbHs3rTM609ZXtrue+jla0KQXXOusY9f1vV+/dzx45z73OMCR9eDxW+272jW2I3nh79NHbDzqebPuXNVW/6rC0REQkvCoAW5ZelYFxO9/p/0KrFn1siJtLOc+MGEWk3+Hr1Dv6z+OgD11dbv2JXzS7S2qVxas6pbS/qrCcOrA/4n99Ag++GZ0d2Gsk9x98DwN+W/I3v8r/zWVsiIhI+FAAtqjEAVvgwAObPg33lEJMEXU/0WTP9shKZeEZvAB6ZuooNOyqP6nzP0i/jeo8j0hbZ9oIiY9z3A8YkwvaFMO2etl/zMC7vczlje4/FxOTeWfeyqnyVT9sTEZHQpwBoUZ61AEurS6mqr/JNIxvdkyroOQbsXghWh3HjSd0Y2T2FfXVObnprSYt3CVm5cyU/lf/kXvql11Es/XIkHbrCxa8ABix+HRa+4r1rN+Oe4+/hhKwTqG6oZsK3E/wzu1tEREKWAqBFJUQlkByTDPhwGHjzD+7n7l4YVj0Cu83gH5cPJjMxhs1lVfzhgxW4XEe+H/DTTZ8CcGbXMxv/Hl7T60w4/RH36y/vOfD38IFIWyRPj36anh16UlZdxq3f3kpl3dH1hIqIiHUoAFqYT/cErt4NRcvdr7uN9v71m5HaPpoXrsojym7j69U7eGHGpsMeX++q5+ut7iVqzuvmm3sUOeH2/YtEO+GDa6D88DW1Rfuo9kw+bTJpsWls3LORO3+4k3qnHxb6FhGRkKMAaGGemcA+6QHcMgswIbU3JGR5//qHMCg7iUd/1Q+Ap79ex8z1Ow957ILiBeyu3U1yTDLHZx7vm4IMA87/O3QaCjV74L1xUFPhm7aAjLgMnj/teWIjYllQvIDH5mt5GBEROZgCoIV51gLc4tji/Yt7hju7nez9ax/B5cfnMHZoNqYJt72/jIJd+5o97sstXwIwpsuYlu3721qRMTDuXUjoBGXr3TODfdgz1zelL0+PfhqbYeOTjZ/w/PLnfdaWiIiEJgVAC/NpD2BjAPTP8O8vPfqrfgzsnMieffVc+8ZCdlc1XYqlpqGGb/PdCzWf0+0c3xcUn+4OgRGxsOlb+HQCuFw+a+6kzifx4PAHAXj5x5f596p/+6wtEREJPQqAFpab4N5FY5tjm3eHCfcUwK5NYNh8uvzL4cRE2nnxqjwyEmLYtLOK6/5vEdV1B/bnnVU4i6r6KjLjMjm247H+KSprEFz6Jhh2+PF9mP4Q+HB49tJel3Lb4NsAeGrxU3y68VOftSUiIqFFAdDCOsd3xmbYqKqvoqy6zHsX3jLD/dwpz70WXoBkJcXy7+uOJyEmgqX5e/j9e8tocLp73TzDv2flnoXN8ON/DXqfBb/6l/v1vOdhzt992tz1A67n6mOuBuDhuQ9roWgREQEUAC0tyh5Fp/adAC9vCRfA+/9+qVd6PK9ecxxRETa+WbODhz79icraSmYUuEPqObl+GP79pUGXw5g/u19/8zAsfctnTRmGwR+G/oELe1yI03Ry94y7WVi80GftiYhIaFAAtDiv7wlsmrB5fw9gEARAgONzk/nn5YOxGfDewgLumfYeda46chNz6d2hd2CKGvl79xIxAFNvg9W+G541DIOHRzzMqdmnUueqY8J3E1hUsshn7YmISPBTALQ4r28JV7oGqkrdkx06H+eda3rBmf0y+H8X9gfg++1fAXB217MxDCNwRZ3+KAy6CkyXe2bwyv/6rKkIWwRPjn6ycbeQW765hblFc33WnoiIBDcFQItrXAzaWz2AnuHfLiMhIto71/SSK4d14fdnZGCP2whAQUGfwK6RZxhwwT/g2MvdC0V/dAMsf9dnzUXbo/n7qX9ndOfR1Dhr+P23v2fm9pk+a09ERIKXAqDFefYE9noADJLh31/q1GkDhuHCWd2Jd+fs49Gpq1u0ZZzP2Ozwq8kw5Gp3T+Ant8DiN3zWXLQ9mr+d/DdOyzmNOlcdt39/uyaGiIhYkAKgxXl6ALdXbm/7tmHOetg2x/06SAOgZ/bv6TlnAvDm3K3c99FKnAENgTY47+9w/I2ACZ/fAQte8llzkfZInhr9FGd2PZMGVwN3/XAXUzdN9Vl7IiISfBQALS69XTpxkXE4TWfbF4QuXAJ1e6FdCqT3906BXlRWXcay0mUA3D/6cp659FhsBkxZXMCt7yxlX11D4Iqz2eDsJ2HEBPfPX/4RvnoAXM7Dn9dKkbZI/jrqr5zX7TwazAbun30/k5dP1rZxIiIWoQBocYZh0C2xGwCbKza37WKe4d/ck9yBJsjM3D4TE5N+Kf3IiMvgkrzO/PPyIUTaDaatKuGSF+axfXfz28b5hWG4l4c55QH3z/Oeh/evhNpKnzQXYYvg8RMf57r+1wHwwooXuH/2/dQ5645wpoiIhLrg+5YWv/MEwE0Vm9p2oSBb/uWXfij4AYCTs09ufO/cgZm8d8NwUttHsabYwa+en8PCLbsCUh/gDoGj/wiXvAYRMbD+S3j9LPfuKj5gM2zckXcHD494GLth5/PNn3Pj9BupqK3wSXsiIhIcFACFbknuALhlz5bWX6S+GrbvX1suNzD7/x5OTUMN84rmAXBK9ilNPhvaNZlPJ5xIv6wEyqvquPLV+by7ID+ww6EDfg3X/g/i0mDHT/DKqZA/32fN/brXr5l8+mTaR7ZnyY4lXP6/y1ldvtpn7YmISGApAIp3hoALl4KrHuIzoUNX7xTmRQuKF1DjrCEjLoNeHXod9HmnpFj+e/NIzh2YSb3T5P6PVzLhvWXs2RfA4dDOQ+GG7yB9gHttxTfOhh/+Ck7f3Ks4Mmsk/z7732TFZVFQWcBVX1zFO2ve0X2BIiJhSAFQ6J7YHYAtFVtwtnbSQcH+3qnsYe5hzCDzfcH3AJzc+eRDLv4cG2Xn+csH88ezemO3Gfzvx2LG/G0mM9bv9GepTSVlw2+nwcBx7mVifpgEb54Lu720cPcv9OzQkw/O/4BTsk+h3lXPXxf+lTu+v0NDwiIiYUYBUMhqn0WULYo6Vx1Fe4tadxHP8GTOcO8V5iUu09W44PEvh39/yTAMbjm5Bx/fMpLuHeMorazlmtcX8tAnPwVulnB0e7j4Jbj4FYiKd4ftF0+EH//j3nrPyxKjE/n7KX/n3uPvJdIWyXcF33Hp1Eu1h7CISBhRABTsNju5iblAK4eBXS4oWOB+nT3Mi5V5x+ry1eys3klcZBxDM4a26JyBnZP4322juHZkVwDemr+NM56dyec/FgVuSHTgZfC72e4t9mod8NH18M6vobyNk3eaYRgGV/a9krfPeZuc+ByKq4q57uvreGD2A+yu2e319kRExL8UAAVo40zgsnVQUwGR7SBjgJcrazvP8O/IrJFE2aNafF5MpJ1HLujH29cNo1NSLIV7qpnw7jLGvjSfnwoDNCTaoSv8ZhqcfB/Yo2DjNzB5OHz3Z6jz/hI2x6Qcw5TzpjC291gMDD7b9BkXfHIBn2z8RPcGioiEMAVAASA3aX8P4J5W9AB6hn87DwV7pBer8g7P8i9HGv49lBN7pvLNxNHceXovYiJtLNy6i/Ofn809//0xMOsG2iPg5Hvhd/Og+6ngrIOZT8G/hsFPH3p98ej2Ue15cPiDvH3O2/Tq0Is9tXt4aM5DXDvtWpaXLvdqWyIi4h8KgAI0nQhy1BqHf4Pv/r+ivUWs370em2FjVKdRrb5ObJSd20/vyXd3ncwFx2Zhmu4dRE5+6gfu+mAFG0t9s1jzYaX2gKs+gsvegsRsqMiH//4WJo+Alf/1ehAc2HEg75/3PhPzJhJjj2Fp6VLGfzme333zO1aVr/JqWyIi4lsKgAI0HQI+6qG9fPf6euQE3/1/nt6/wWmDSYpJavP1spJi+cflg/nwdyM4sUcqDS6TD5du54y/zeTmt5awZNsu/w6NGgYccwHcugBOvh9iEt1D8h9e5+4RXDEFGry3lE2kLZLf9P8NUy+ayiU9L8Fu2JldOJtxn4/jju/v4Keyn7zWloiI+I5h6kaeVnM4HCQmJlJRUUFCQkKgy2mTemc9x71zHE7TyTe//ob0uPSWnVi5A57pBRhw7zZ3AAkiN359I/OK53FX3l1c2/9ar19/RcEeJv+wka9W7Wh8r2dae8Yel81FgzuR0j7a620eVk0FLHjZvY1czR73e3EdYdCVkHcNJHfzanP5jnxeXPEin2/+HBP3/5T0T+nPuD7jOCv3LKLtfv79RURaIJy+v1tLAbANwu0f0Pkfn89Wx1ZeOuMlRmaNbNlJqz+FD66G9P7wuzm+LfAo7a3by6gpo2hwNTD1wql0Tezqs7Y27KjkpZmb+fzHImrqXQBE2g1O75vOOQMyObl3R+Jj/Hh/ZI0DFr4Mi16FyuID73c72R0Ge53p1bC+ac8mXlv5GtO2TqPeVQ9AUnQSF/a4kLNzz6Zvct9Drr8oIuJv4fb93RoKgG0Qbv+Abv/udr4r+I57j7+XK/te2bKTpt0P8/8FQ6+D8571bYFH6autX/GHGX+ga0JXpl401S9tOmrqmbqiiA8WFbBi+4GZwpF2gxHdUxlzTDqn9kkjKynWL/XgbID102DJG7DxW9jfS4c9yh0G+54Pvc+BuFSvNFdeXc7HGz/mg3UfUFx1IHhmx2czpssYxnQdozAoIgEXbt/fraEA2Abh9g/oH0v/wSsrX+GyXpfx0IiHWnbSy6dA0VK4+FUYeKlvCzxKD8x+gM82fcY1x1zDH477g9/bX1vi4JNlRUxfXcKmnVVNPstJbsfwbskM75bCsG4pZCXG+D4U7d4Gy95299qWrfvZBwZkDoSuo9yPLiPa3DvodDmZsX0Gn2/+nFnbZ1HjrGn8LC02jWGZwxofGXEZbWpLRORohdv3d2uEVQCcPHkyTz31FMXFxfTr14/nnnuOUaMOPfNzxowZTJw4kVWrVpGVlcUf//hHbr755ha3F27/gKZumsr9s+8nLz2PN89688gn1O2Dv2aDqwHuWAlJOT6vsaVcpotTPziV8ppyXh3zKsMyAztBZWPpXqav3sHXq0tYUbAH1y/+W5faPpr+nRIY0CmRflmJHJOZQKcOsdhtPgqFO9fBms9gzVQoXtH0M8MG6f0g81jIHOR+ZPSHyNb1Wu6r38fMwpl8vfXrg8IgQE58Dv1T+3NMyjH0S+lH35S+xEXGte73EhFpgXD7/m6NsAmAU6ZMYfz48UyePJkTTjiBl156iVdffZXVq1eTk3NwMNmyZQv9+/fnhhtu4KabbmLOnDnccsstvPfee1xyySUtajPc/gGtKl/FuM/HkRyTzIyxM458wpZZ8H/nQXwWTFwdVHsAr921lkunXkpsRCyzx80+qgWgfa2ypp7FW3czf0s58zfv4qfCCpy/TIRAVISNrint6Jbanm4d4+jcoR2dOsTSKSmGzMRY4qIjvFOQoxi2zoats9yPXc2sBWnY3AE/pSek9oSU7u4JJQmdIbETRLUssNU6a1leupwFxQtYULyAn8p/wmW6mjaFQaf2neia2JWuCV3JTcyla0JXMttnkt4uPaj+sxSR0BRu39+tETYBcNiwYQwZMoQXXnih8b2+ffty4YUXMmnSpIOOv+eee/jss89Ys2ZN43s333wzK1asYN68eS1qM9z+Ae2r38ewd909ZTPHzqRDTIfDnzDzKfcOFP0ugkvf9H2BR+HVla/y96V/5+TOJ/PP0/4Z6HIOq7rOyZoSB6sKK/ip0MFPRRVsKN1LXYPrsOclxESQGh9NavtoOraPJrV9FEntokiMjSSpXSSJse5HXHQE7aMjaBdlJy46gugI2+GHmysK3cP6RcvdvYPFy6Fq5+F/idgO7jDYPs096zgu1f1olwIxSe4hZc8jOt4dGCNiqKzfy4qdK1hdvprV5atZVb6KkqqSwzaVGptKRrsMOrbrSHJMcuOjQ0wHEqISiI+Kb3y0j2xPTEQMNkMrXonIAeH2/d0aXupCCKy6ujqWLFnCvffe2+T9MWPGMHfu3GbPmTdvHmPGjGny3plnnslrr71GfX09kZHBt6OFr7WLbEdWXBZFVUVsrthMXkze4U/ID94FoOcUumckn9DphABXcmSxUXaG5HRgSM6BwO10mRTtqWbTzr1s3lnFlrIqCvdUU7SnmsI91VTWNODY/9j8i/sLj8RmQGykndgoO9ERdmIibURF2ImKsBFttxEVYSPSnk6E/Wwi7ecQ0cmgg7mbzPoC0uq3k1aXT1ptAUl1xSTW7SDatQ+qd7sfO47cvocLG9H2dgy2xzDQFo3TFoXTHkOZPZGtETbyI0wK7C6225wU2eopM+qpM0zKqssoqy6D8pa3FYWdaOxEE0EUdiKxE2HYiMROJDbs2IjAjt0wsGPHhoEdm/vZcD8bGNj2PzyvwcAG+98xwGD/T573PA689rxqGsGbvmfgnd704OmTF2m9gV1O5pwTrw50GWEnLAJgWVkZTqeT9PSma9elp6dTUtJ8b0JJSUmzxzc0NFBWVkZmZuZB59TW1lJbW9v4s8Ph8EL1waVbUrcDATD9MAHQ5YKChe7XQbYA9N66vY1blIVCAGyO3WaQndyO7OR2nNz74M8dNfXsqKhh595ayvbWUb63lrK9tezZV8+e6noc1fVU7H9U1Tqpqm2gut69M4jLhKo6J1V1R7tTSMf9j8FN3o1nH5lGOVlGOSk4SDEqSDEcpBoOOlBJvLGPBPaRaFSRwD7aGe7/DtlwEeXcS5Rzb5PrpQJ9mmndBHbbbJRE2CmOiKDcbmeX3cYum/t5t91Opc1Gpc3Y/2zDtb+nsw4ndTippK7pBZtrRESCimPDbgVAHwiLAOjxy2Et0zQPO9TV3PHNve8xadIkHn300TZWGdy6JXZjduHsI+8JvHMN1FZAZBykD/BPcS20oGQBDWYDXRK6kB2fHehyfCIhJpKEmEh6pse3+Byny6SqroGaOifV9e5HTb2L6jondU4XdQ37H04n9Q0m9S4XDU6TeqeLBpeJ85cP08TlMnGZJk4XuEwT0zSpNiHfNNkGuFwmpgkm7meXCYbpJNJVQ5SzmmizmkhXNRGuWiJddUSYdUS4aokw67GbDdjMBiLMBuxmHXbTiQ0nNtOJ3WwgASeJpkl304mtwYVR78K2/35CAxcGLurNBuoMF3WG0/3ASb3hogEXDYaLelw0GCZOoMFw4cR0PwwTFyYuwImJy3D/Di7cGdFluJe9NvE8/+y1ceBnj1/mSrPZ12azxx7qvYOou0/CVHbH5v4vobRVWATA1NRU7Hb7Qb19paWlB/XyeWRkZDR7fEREBCkpKc2ec9999zFx4sTGnx0OB9nZ4RUwPFvCba44QgDMn+9+7jwU7MH1z6hx+DcrNHv/fMVuMxqDo4iIWFtY3BkdFRVFXl4e06dPb/L+9OnTGTmy+R0tRowYcdDxX3/9NUOHDj3k/X/R0dEkJCQ0eYSb7kndgRYEwMIl7ufs431c0dExTbMxAJ7Y6cQAVyMiIhKcwiIAAkycOJFXX32V119/nTVr1nDnnXeSn5/fuK7ffffdx9VXH7iH4Oabb2bbtm1MnDiRNWvW8Prrr/Paa6/xhz/4f8HgYJKbmAtASVUJVfWHmVzgCYCdhvqhqpbb4thCUVURUbYohmYEV20iIiLBIrjG7tpg7NixlJeX89hjj1FcXEz//v354osv6NKlCwDFxcXk5+c3Hp+bm8sXX3zBnXfeyb/+9S+ysrL4xz/+0eI1AMNVYnQiKTEplNeUs6ViC/1T+x98UI3DvZAwQKch/i3wCDy9f0MzhhIb4aft1kREREJM2ARAgFtuuYVbbrml2c/efPPNg94bPXo0S5cu9XFVoad7UnfKS8rZXLG5+QBYvBwwITHHve5bEJldOBvQ/X8iIiKHEzZDwOI9nokg63etb/6AxuHfwc1/HiDVDdUsLlkM6P4/ERGRw1EAlIP0SXZPuV+3e13zBzQGwCMsFO1ni0sWU+eqIzMus/FeRhERETmYAqAcpDEA7lpHszsFFi5zPwdZAJxTdGD3j8NudSYiImJxCoBykO5J3bEbdnbX7mZn9S/2gK0sAcd2MGyQOSgg9R1K4/IvWRr+FRERORwFQDlITERM4xDq2l1rm35YuH/STMc+EN3ez5Ud2vbK7Wx1bCXCiOD4zOBam1BERCTYKABKs3onuzegXbfrF/cBNt7/F1zLv8wtmgvAwI4DiY9q+fZoIiIiVqQAKM3q08F9H+DBPYD7A2BWcAVAz/Ivmv0rIiJyZAqA0qzGHsCfzwR2uaBo/xBwEE0AqXfVs7BkIQAjOzW/9Z+IiIgcoAAozfIEwHxH/oEt4XZthpoKsEdDer8AVtfUitIVVNVXkRyTTN/kvoEuR0REJOgpAEqzkmOSSWuXhonJht0b3G96ev8yjwV7ZOCK+wXP8i8jskZgM/RPWkRE5Ej0bSmH5FkPsPE+wCBdANqz/Iu2fxMREWkZBUA5pN4d3MPAwRwAy6vLWbNrDeDuARQREZEjUwCUQ/r5jiA01EHxj+4PgmgJGM/yL32S+5AamxrgakREREKDAqAckicAbtizgYaSH8FZCzGJkNwtwJUd4AmAGv4VERFpOQVAOaTO8Z1pF9GOWmct27Z+636zUx4EyT67LtN1IAB2UgAUERFpKQVAOSSbYWtcDmZt0SL3m0F0/9/aXWvZVbOLdhHtGNRxUKDLERERCRkKgHJYnokg6yo2ud8IogDo6f07PvN4IoNoWRoREZFgpwAoh9W4FIxzr/uNINoCzrP9m+7/ExEROToKgHJYjTOBoyIx47MgPj3AFbntrdvLitIVgAKgiIjI0VIAlMPqntQdOwa77HZ2ZgbP9m8LSxbSYDaQHZ9NdkJ2oMsREREJKQqAclgxETHk2mIAWNshK8DVHKDdP0RERFpPAVCOqHdtHQDromMCXImbaZrMLJwJwKjOowJcjYiISOhRAJTDq91LH0c5AGtdVQEuxm397vWUVJUQY4/h+IzjA12OiIhIyFEAlMMrWUnvuloA1lVuC3AxbjO3u3v/hmUOIyYiOHolRUREQokCoBxe8XL61NUDsM2xDUedI8AFwYztMwA4qfNJAa5EREQkNCkAyuEVLaeDy0WXiHgAlpcuD2g5u2p28ePOHwEFQBERkdZSAJTDK14OwJAO7vUAl+5YGsBi3Is/m5j0Se5DRlxGQGsREREJVQqAcmh1VVC2HoAhOaMBWLJjSSArarz/T71/IiIiracAKIdWshJMF7TPIC/nVAB+Kv+JmoaagJRT76pvXP9PAVBERKT1FADl0IqWu5+zBtE5vjMdYzvS4GpgZdnKgJSzbMcy9tbvJTkmmf4p/QNSg4iISDhQAJRD23//H5mDMAyDIelDgMDdB+iZ/XtipxOx2+wBqUFERCQcKADKof2sBxBgSNr+AFgamADouf9vdOfRAWlfREQkXCgASvPq9kHZOvfrzEEA5KXnAe6lYBpcDX4tZ5tjG1sdW4kwIhiZNdKvbYuIiIQbBUBp3o6f9k8ASYeETAB6JPUgPjKefQ37WLd7nV/L8fT+5WXk0T6qvV/bFhERCTcKgNI8z/Dv/t4/ALvNzqA098/+vg/Qc/+fhn9FRETaTgFQmueZALL//j+PQEwE2Vu3lyUl7vUHFQBFRETaTgFQmtdMDyAcuA9waelSTNP0SynfF3xPg9lAbmIuOQk5fmlTREQknCkAysHqq2HnWvfrzGObfNQvpR9Rtih21exiq2OrX8r5csuXAJzV9Sy/tCciIhLuFADlYCU/gemEuI6QkNXkoyh7FAM6DgD8Mwy8p2YP84rmAXB27tk+b09ERMQKFADlYD9bABrDOOhjf64H+PW2r2kwG+ib3JfcxFyftyciImIFCoBysENMAPHw3Ae4ZMcSn5cybes0AM7K1fCviIiItygAysGKV7iffzEBxOPYjsdiM2wU7i2kpKrEZ2XsqNrB4pLFAJzdVcO/IiIi3qIAKE011ELpGvfrX0wA8Wgf1Z7eHXoDvr0P8KutX2FiMjhtMJntM33WjoiIiNUoAEpTpavB1QCxyZDY+ZCHeYaB5xfP91kpjcO/mv0rIiLiVQqA0lTj8O+xzU4A8Tgl+xQAvs3/lnpnvdfLKHAUsLJsJTbDxpiuY7x+fREREStTAJSmfh4ADyMvPY+OsR1x1DmYUzTH62V8udW99t+wjGGkxqZ6/foiIiJWpgAoTbUwANptds7seiYAX2z5wutleBZ/1tp/IiIi3qcAKAc4692LQMMRAyDAObnnAPBDwQ/sq9/ntTI27N7Axj0bibRFclqX07x2XREREXFTAJQDytaDsxaiE6DDkRdd7p/an+z4bKobqvmh4AevleHp/Tux04kkRCV47boiIiLipgAoB3iGfzMGgu3I/zQMw2icoesJbW1V01DDRxs+Ag70MIqIiIh3KQDKAS28/+/nPCFtdtFsKmor2lzCpxs/pbymnIy4DA3/ioiI+IgCoBzQigDYo0MPenXoRYOrgW+2fdOm5htcDbyx6g0Aru13LZG2yDZdT0RERJqnAChuLhcU/+h+fRQBEA7M1G3rMPBXW7+icG8hHaI7cHHPi9t0LRERETk0BUBxK98I9VUQEQupPY/qVE8AXFiykJ37draqedM0ee2n1wC4su+VxEbEtuo6IiIicmQKgOLWOAFkANjsR3Vqp/adOLbjsZiYfLX1q1Y1P3P7TDbs3kC7iHaM6zOuVdcQERGRllEAFLfi5e7noxz+9fD0ArZ2UWhP799lvS8jMTqxVdcQERGRllEAFLdWTAD5uTO7nonNsLGybCVLdiw5qnOX7FjCstJlRNoiGX/M+Fa1LyIiIi2nAChgmq2eAOKRGpvKRT0uAuChOQ8d1c4gr658FYBf9fgVae3SWtW+iIiItJwCoMDurVBbAfYo6Nin1Ze5a+hdZMRlUFBZwN+W/K1F58wtmsvswtnYDBu/6febVrctIiIiLacAKAeGf9P7QURUqy8THxXPYyMfA+D9de8zv3j+YY9fXb6aO7+/E4ALe1xITkJOq9sWERGRllMAlDbf//dzI7JGMLb3WAD+NOdP7K3b2+xx+Y58fvfN79jXsI9hGcN4YNgDbW5bREREWkYBULwaAAEm5k2kc/vOFFcV8/Tipw/6vKy6jBun38iuml30Te7Lc6c8R5S99T2PIiIicnQiAl2ABJhpej0Atotsx/874f/x269+y4cbPiTaHk33pO5kxmXSsV1HHpz9IIV7C8mOz2by6ZNpH9XeK+2KiIhIyygAWp2jCPaVgWGHtH5eu+zQjKFcdcxVvLX6Ld5d++5BnyfHJPPS6S+RGpvqtTZFRESkZRQArc6zAHRaX4iM8eqlJ+ZNpEdSD9buWkvx3mIKqwop2ltEtD2aF09/keyEbK+2JyIiIi2jAGh1Rcvdz5mDvH7pCFsEF/e8uMl7pmkCYBiG19sTERGRlgmLSSC7d+9m/PjxJCYmkpiYyPjx49mzZ88hj6+vr+eee+5hwIABxMXFkZWVxdVXX01RUZH/ig4Wnh7ArEF+ac4wDIU/ERGRAAuLAHjFFVewfPlypk2bxrRp01i+fDnjxx96S7F9+/axdOlSHnroIZYuXcpHH33E+vXrueCCC/xYdRAwTSha5n7tgx5AERERCU6G6RmTC1Fr1qzhmGOOYf78+QwbNgyA+fPnM2LECNauXUvv3r1bdJ1FixZx/PHHs23bNnJyWrYgscPhIDExkYqKChISElr9OwRMRSH87Rj3BJD7CyEyNtAViYiI+FzIf397Qcj3AM6bN4/ExMTG8AcwfPhwEhMTmTt3bouvU1FRgWEYJCUlHfKY2tpaHA5Hk0dI8wz/duyj8CciImIhIR8AS0pKSEtLO+j9tLQ0SkpKWnSNmpoa7r33Xq644orD/j+BSZMmNd5nmJiYSHZ2iM9i9UwA8dP9fyIiIhIcgjYAPvLII40TBg71WLx4MdD8jFLTNFs02aC+vp5x48bhcrmYPHnyYY+97777qKioaHwUFBS07pcLFo0TQAYHtAwRERHxr6BdBmbChAmMGzfusMd07dqVH3/8kR07dhz02c6dO0lPTz/s+fX19Vx22WVs2bKF77777oj3AURHRxMdHX3k4kOBafp0CRgREREJXkEbAFNTU0lNPfIuESNGjKCiooKFCxdy/PHHA7BgwQIqKioYOXLkIc/zhL8NGzbw/fffk5KS4rXaQ4KjCKpK3RNAMvoHuhoRERHxo6AdAm6pvn37ctZZZ3HDDTcwf/585s+fzw033MB5553XZAZwnz59+PjjjwFoaGjg17/+NYsXL+add97B6XRSUlJCSUkJdXV1gfpV/EsTQERERCwr5AMgwDvvvMOAAQMYM2YMY8aMYeDAgbz11ltNjlm3bh0VFRUAbN++nc8++4zt27czaNAgMjMzGx9HM3M4pGkCiIiIiGUF7RDw0UhOTubtt98+7DE/X+6wa9euhPjyh23n6QHU/X8iIiKWExY9gHKUfj4BRD2AIiIilqMAaEWNE0BskK4JICIiIlajAGhFjRNA+kJUu4CWIiIiIv6nAGhFGv4VERGxNAVAK9IEEBEREUtTALQaTQARERGxPAVAq6ks1gQQERERi1MAtJqiZe7njn00AURERMSiFACtxjP8q/v/RERELEsB0Go8E0B0/5+IiIhlKQBaiWlC4RL36055ga1FREREAkYB0Ep2b4V95WCL1AQQERERC1MAtBJP71/GAIiMCWwtIiIiEjAKgFbiCYCdhwa2DhEREQkoBUAr2b7Y/dxJAVBERMTKFACtoqEOile4X2sCiIiIiKUpAFpF6Spw1kJMEqR0D3Q1IiIiEkAKgFbROPybB4YR2FpEREQkoBQArULr/4mIiMh+CoBWoRnAIiIisp8CoBVU74Gy9e7X6gEUERGxPAVAKyha6n5O6gJxqYGtRURERAJOAdAKNPwrIiIiP6MAaAXbPRNAFABFREREATD8mSYU/mwJGBEREbE8BcBwV1EAVTvBFgGZAwNdjYiIiAQBBcBw51kAOr0/RMYGthYREREJCgqA4U4LQIuIiMgvKACGO80AFhERkV9QAAxnznooWu5+rRnAIiIisp8CYDgrXQ0N1RCdCCk9Al2NiIiIBAkFwHCWv8D93DkPbPqPWkRERNyUCsLZtjnu5y4jA1uHiIiIBBUFwHBlmrBtrvt1lxMCW4uIiIgEFQXAcLVrM1SVgj0KsoYEuhoREREJIgqA4coz/NtpKETGBLYWERERCSoKgOGqcfhX9/+JiIhIUwqA4UoTQEREROQQFADDUcV22JMPhh2yjw90NSIiIhJkFADD0bZ57ufMgRAdH9haREREJOgoAIajxuFfLf8iIiIiB1MADEeaACIiIiKHoQAYbqrKoGyd+3XOiMDWIiIiIkFJATDc5O+//y/tGGiXHNhaREREJCgpAIYbz/Cvev9ERETkEBQAw43W/xMREZEjUAAMJzUOKFnpfq0AKCIiIoegABhOChaC6YIOuZCQFehqREREJEgpAIYTrf8nIiIiLaAAGE4aA6AmgIiIiMihKQCGi+rdsH2x+3XuSYGtRURERIKaAmC42PgtmE7o2BeScgJdjYiIiAQxBcBwseFr93OvMYGtQ0RERIKeAmA4cDlhw3T3655nBrYWERERCXoKgOFg+2Ko3gUxiZA9LNDViIiISJBTAAwHG75yP3c/DewRga1FREREgp4CYDhY77n/T8O/IiIicmQKgKGuohB2rAQM6HF6oKsRERGREKAAGOo8s387D4W41MDWIiIiIiFBATDUeQKgZv+KiIhICykAhrL6Gtj8g/u11v8TERGRFlIADGXbZkP9PmifARkDA12NiIiIhAgFwFDmmf3b8wwwjMDWIiIiIiFDATBUmeaB9f+0/IuIiIgcBQXAUFW2AXZvBVskdDs50NWIiIhICAmLALh7927Gjx9PYmIiiYmJjB8/nj179rT4/JtuugnDMHjuued8VqPXrf7E/dz1RIiOD2gpIiIiElrCIgBeccUVLF++nGnTpjFt2jSWL1/O+PHjW3TuJ598woIFC8jKyvJxlV7kcsLSf7tfH3t5YGsRERGRkBPyG8euWbOGadOmMX/+fIYNGwbAK6+8wogRI1i3bh29e/c+5LmFhYVMmDCBr776inPPPddfJbfdpu+gogBikuCYXwW6GhEREQkxId8DOG/ePBITExvDH8Dw4cNJTExk7ty5hzzP5XIxfvx47r77bvr16+ePUr1nyZvu50FXQGRMQEsRERGR0BPyPYAlJSWkpaUd9H5aWholJSWHPO+JJ54gIiKC2267rcVt1dbWUltb2/izw+E4umK9wVEM6750vx5yjf/bFxERkZAXtD2AjzzyCIZhHPaxePFiAIxm1sAzTbPZ9wGWLFnC3//+d958881DHtOcSZMmNU40SUxMJDs7u3W/XFssextMJ+SMgLQ+/m9fREREQl7Q9gBOmDCBcePGHfaYrl278uOPP7Jjx46DPtu5cyfp6enNnjdr1ixKS0vJyclpfM/pdHLXXXfx3HPPsXXr1mbPu++++5g4cWLjzw6Hw78h8OeTP/Ku9V+7IiIiElaCNgCmpqaSmpp6xONGjBhBRUUFCxcu5PjjjwdgwYIFVFRUMHLkyGbPGT9+PKeffnqT984880zGjx/Pb37zm0O2FR0dTXR09FH8Fl626XuoyNfkDxEREWmToA2ALdW3b1/OOussbrjhBl566SUAbrzxRs4777wmM4D79OnDpEmTuOiii0hJSSElJaXJdSIjI8nIyDjsrOGAW/KG+/nYyyEyNrC1iIiISMgK2nsAj8Y777zDgAEDGDNmDGPGjGHgwIG89dZbTY5Zt24dFRUVAarQCypLDkz+yNPkDxEREWm9kO8BBEhOTubtt98+7DGmaR7280Pd9xc0lr3lnvyRPRzS+ga6GhEREQlhYdEDGPZqK2Hx/uFfTf4QERGRNlIADAXT7gVHISTmQL8LA12NiIiIhDgFwGC35nP32n8YcNGLmvwhIiIibaYAGMwqd8DU/TuVnHAbdD0hsPWIiIhIWFAADFamCZ9NgH3lkD4ATnkg0BWJiIhImFAADFZL3oANX4M9Gi5+GSICuAC1iIiIhBUFwGBUthG+2t/jd/rDkH5MYOsRERGRsKIAGIy+fQTq90HuSTDsd4GuRkRERMJMWCwEHXYueB7apcBJd4NNGV1ERES8SwEwGMUmwfl/D3QVIiIiEqbUvSQiIiJiMQqAIiIiIhajACgiIiJiMQqAIiIiIhajACgiIiJiMQqAIiIiIhajACgiIiJiMQqAIiIiIhajACgiIiJiMQqAIiIiIhajACgiIiJiMQqAIiIiIhajACgiIiJiMRGBLiCUmaYJgMPhCHAlIiIi0lKe723P97gVKQC2QWVlJQDZ2dkBrkRERESOVmVlJYmJiYEuIyAM08rxt41cLhdFRUXEx8djGEarr+NwOMjOzqagoICEhAQvVijN0d/bv/T39i/9vf1Lf2//8tbf2zRNKisrycrKwmaz5t1w6gFsA5vNRufOnb12vYSEBP0PiB/p7+1f+nv7l/7e/qW/t3954+9t1Z4/D2vGXhERERELUwAUERERsRgFwCAQHR3Nww8/THR0dKBLsQT9vf1Lf2//0t/bv/T39i/9vb1Hk0BERERELEY9gCIiIiIWowAoIiIiYjEKgCIiIiIWowAoIiIiYjEKgAE2efJkcnNziYmJIS8vj1mzZgW6pLA1adIkjjvuOOLj40lLS+PCCy9k3bp1gS7LEiZNmoRhGNxxxx2BLiWsFRYWctVVV5GSkkK7du0YNGgQS5YsCXRZYamhoYEHH3yQ3NxcYmNj6datG4899hgulyvQpYWFmTNncv7555OVlYVhGHzyySdNPjdNk0ceeYSsrCxiY2M5+eSTWbVqVWCKDVEKgAE0ZcoU7rjjDh544AGWLVvGqFGjOPvss8nPzw90aWFpxowZ3HrrrcyfP5/p06fT0NDAmDFjqKqqCnRpYW3RokW8/PLLDBw4MNClhLXdu3dzwgknEBkZyZdffsnq1at55plnSEpKCnRpYemJJ57gxRdf5Pnnn2fNmjU8+eSTPPXUU/zzn/8MdGlhoaqqimOPPZbnn3++2c+ffPJJnn32WZ5//nkWLVpERkYGZ5xxBpWVlX6uNHRpGZgAGjZsGEOGDOGFF15ofK9v375ceOGFTJo0KYCVWcPOnTtJS0tjxowZnHTSSYEuJyzt3buXIUOGMHnyZP785z8zaNAgnnvuuUCXFZbuvfde5syZo1EEPznvvPNIT0/ntddea3zvkksuoV27drz11lsBrCz8GIbBxx9/zIUXXgi4e/+ysrK44447uOeeewCora0lPT2dJ554gptuuimA1YYO9QAGSF1dHUuWLGHMmDFN3h8zZgxz584NUFXWUlFRAUBycnKAKwlft956K+eeey6nn356oEsJe5999hlDhw7l0ksvJS0tjcGDB/PKK68EuqywdeKJJ/Ltt9+yfv16AFasWMHs2bM555xzAlxZ+NuyZQslJSVNvj+jo6MZPXq0vj+PQkSgC7CqsrIynE4n6enpTd5PT0+npKQkQFVZh2maTJw4kRNPPJH+/fsHupyw9P7777N06VIWLVoU6FIsYfPmzbzwwgtMnDiR+++/n4ULF3LbbbcRHR3N1VdfHejyws4999xDRUUFffr0wW6343Q6efzxx7n88ssDXVrY83xHNvf9uW3btkCUFJIUAAPMMIwmP5umedB74n0TJkzgxx9/ZPbs2YEuJSwVFBRw++238/XXXxMTExPocizB5XIxdOhQ/vKXvwAwePBgVq1axQsvvKAA6ANTpkzh7bff5t1336Vfv34sX76cO+64g6ysLK655ppAl2cJ+v5sGwXAAElNTcVutx/U21daWnrQ/6sR7/r973/PZ599xsyZM+ncuXOgywlLS5YsobS0lLy8vMb3nE4nM2fO5Pnnn6e2tha73R7ACsNPZmYmxxxzTJP3+vbty4cffhigisLb3Xffzb333su4ceMAGDBgANu2bWPSpEkKgD6WkZEBuHsCMzMzG9/X9+fR0T2AARIVFUVeXh7Tp09v8v706dMZOXJkgKoKb6ZpMmHCBD766CO+++47cnNzA11S2DrttNNYuXIly5cvb3wMHTqUK6+8kuXLlyv8+cAJJ5xw0LJG69evp0uXLgGqKLzt27cPm63pV6jdbtcyMH6Qm5tLRkZGk+/Puro6ZsyYoe/Po6AewACaOHEi48ePZ+jQoYwYMYKXX36Z/Px8br755kCXFpZuvfVW3n33XT799FPi4+Mbe18TExOJjY0NcHXhJT4+/qB7K+Pi4khJSdE9lz5y5513MnLkSP7yl79w2WWXsXDhQl5++WVefvnlQJcWls4//3wef/xxcnJy6NevH8uWLePZZ5/lt7/9baBLCwt79+5l48aNjT9v2bKF5cuXk5ycTE5ODnfccQd/+ctf6NmzJz179uQvf/kL7dq144orrghg1SHGlID617/+ZXbp0sWMiooyhwwZYs6YMSPQJYUtoNnHG2+8EejSLGH06NHm7bffHugywtrUqVPN/v37m9HR0WafPn3Ml19+OdAlhS2Hw2HefvvtZk5OjhkTE2N269bNfOCBB8za2tpAlxYWvv/++2b/9/qaa64xTdM0XS6X+fDDD5sZGRlmdHS0edJJJ5krV64MbNEhRusAioiIiFiM7gEUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsRgFQBERERGLUQAUERERsZj/D+KX+P2Tg+sHAAAAAElFTkSuQmCC' width=640.0/>\n",
       "            </div>\n",
       "        "
      ],
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = 0\n",
    "L, r, G = radial_fdm_laplacian(10, 100, left_bc='neumann_simple' if m == 0 else 'dirichlet')\n",
    "R = np.diag(r[1:-1])\n",
    "Rm2 = np.diag(1/r[1:-1]**2)\n",
    "Rm5 = np.diag(1/r[1:-1]**.5)\n",
    "\n",
    "E, U = eigh(-.5*L  + 0.5*m**2*Rm2 + .5*R**2)\n",
    "ic(E[:])\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(r[1:-1], Rm5 @ U[:,0])\n",
    "plt.plot(r[1:-1], Rm5 @ U[:,1])\n",
    "plt.plot(r[1:-1], Rm5 @ U[:,2])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "class polar_fdm_2d:\n",
    "    \"\"\"Class for 2d polar FDM SE solver.\n",
    "    \n",
    "    Attributes:\n",
    "        r_max (float): right endpoint\n",
    "        n_r (int): number of inner grid points\n",
    "        m_max (int): maximum angular momentum\n",
    "        L_neumann (ndarray): radial laplacian with Neumann BC\n",
    "        L_0 (ndarray): radial laplacian with Dirichlet BC\n",
    "        r0 (ndarray): radial grid with endpoints\n",
    "        r_inner (ndarray): inner radial grid points\n",
    "        G_neumann (ndarray): transition matrix from inner to global grid with Neumann BC\n",
    "        G_0 (ndarray): transition matrix from inner to global grid with Dirichlet BC\n",
    "        T_m (list of ndarray): Kinetic energy for each |m|\n",
    "        V_m (list of ndarray): list of potential matrices V(r) for each m.\n",
    "        V_m_max (int): number of potential matrices\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, r_max, n_r, m_max):\n",
    "        \"\"\"Initialize the polar FDM grid.\n",
    "        \n",
    "        Args:\n",
    "            r_max (float): right endpoint\n",
    "            n_r (int): number of inner grid points\n",
    "            m_max (int): maximum angular momentum\n",
    "        \"\"\"\n",
    "        self.r_max = r_max\n",
    "        self.n_r = n_r\n",
    "        self.m_max = m_max\n",
    "\n",
    "        # Set up radial laplacian for Neumann and Dirichlet BC\n",
    "        self.L_neumann, self.r0, self.G_neumann = radial_fdm_laplacian(r_max, n_r, left_bc = 'neumann')\n",
    "        self.L_neumann = csr_matrix(self.L_neumann)\n",
    "        self.L_0, _, self.G_0 = radial_fdm_laplacian(r_max, n_r, left_bc = 'dirichlet')\n",
    "        self.L_0 = csr_matrix(self.L_0)\n",
    "        self.r_inner = self.r0[1:-1]        \n",
    "        \n",
    "        # Set up Kinetic energy for each |m|\n",
    "        self.T_m = []\n",
    "        for m in range(m_max+1):\n",
    "            if m == 0:\n",
    "                self.T_m.append( -0.5*self.L_neumann  )\n",
    "            else:\n",
    "                # THE POWER IS PROBABLY WRONG HERE\n",
    "                self.T_m.append( -0.5*self.L_0 + np.diag(.5*m**2/self.r_inner**2) )\n",
    "                \n",
    "                \n",
    "        \n",
    "        # Set wavefunction shape\n",
    "        self.shape = (n_r, 2*m_max + 1)\n",
    "        # Set number of degrees of freedom\n",
    "        self.n_dof = np.prod(self.shape)\n",
    "\n",
    "            \n",
    "\n",
    "    def set_potential(self, V_m):\n",
    "        \"\"\" Set the scalar potential.\n",
    "        \n",
    "        \n",
    "        Args:\n",
    "            V_m (list of ndarray): list of potential matrices V(r) for each m.\n",
    "        \"\"\"\n",
    "        self.V_m = V_m\n",
    "        self.V_m_max = (len(V_m) - 1)//2\n",
    "        ic(self.V_m_max)\n",
    "\n",
    "\n",
    "    def apply_weight(self, psi):\n",
    "        \"\"\"RHS of the SE with weight function, generalized EVP\"\"\"\n",
    "        \n",
    "        result = np.zeros(self.shape, dtype=np.complex128)\n",
    "        for m in range(-self.m_max, self.m_max+1):\n",
    "            result[:,m+self.m_max] = self.r_inner * psi[:,m+self.m_max] \n",
    "            \n",
    "        return result\n",
    "        \n",
    "    def apply_hamiltonian(self, psi, gen_evp=False):\n",
    "        \"\"\"Apply the Hamiltonian to a wavefunction.\n",
    "        \"\"\"\n",
    "        \n",
    "        \n",
    "        assert(psi.shape == self.shape)\n",
    "        result = np.zeros(self.shape, dtype=np.complex128)\n",
    "        \n",
    "        for m in range(-self.m_max, self.m_max+1):\n",
    "            m_index = m + self.m_max\n",
    "            #ic(m, m_index)\n",
    "            result[:,m_index] = self.T_m[abs(m)] @ psi[:,m_index] \n",
    "            for m2 in range(-self.V_m_max, self.V_m_max+1):\n",
    "                m2_index = self.V_m_max + m2\n",
    "                m3 = m - m2\n",
    "                if m3 >= -self.m_max and m3 <= self.m_max:\n",
    "                    m3_index = m3 + self.m_max\n",
    "                    #ic(m_index, m2_index, m3_index)\n",
    "                    result[:,m_index] += self.V_m[m2_index] * psi[:,m3_index]\n",
    "            if gen_evp:\n",
    "                result[:,m_index] *= self.r_inner\n",
    "        \n",
    "        return result\n",
    "    \n",
    "    def get_dense_matrix(self):\n",
    "        \"\"\"Compute dense matrix representation of H.\"\"\"\n",
    "        def apply_hamiltonian(psi_vec):\n",
    "            return self.apply_hamiltonian(psi_vec.reshape(self.shape)).flatten()\n",
    "\n",
    "        H = LinearOperator((self.n_dof, self.n_dof), matvec=apply_hamiltonian, dtype=np.complex128)\n",
    "        \n",
    "        H_mat = np.zeros((self.n_dof, self.n_dof), dtype=np.complex128)\n",
    "        R_mat = np.zeros((self.n_dof, self.n_dof), dtype=np.complex128)\n",
    "        \n",
    "        ic(self.n_dof)\n",
    "        for i in range(self.n_dof):\n",
    "            e = np.zeros(self.n_dof)\n",
    "            e[i] = 1\n",
    "            H_mat[:,i] = self.apply_hamiltonian(e.reshape(self.shape), gen_evp=False).flatten()\n",
    "            R_mat[:,i] = self.apply_weight(e.reshape(self.shape)).flatten()\n",
    "            \n",
    "        return H_mat, R_mat\n",
    "\n",
    "        \n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ic| self.V_m_max:"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 1\n",
      "ic| self.n_dof: 2100\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7b7289234f8449478f45fc412530373f",
       "version_major": 2,
       "version_minor": 0
      },
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkxklEQVR4nO3de4yV5Z3A8d/IZYpkOHHEmcOsIyENurawJItdGLZb7yMkI7VuVlqaCSasl1UxEzB23WYjbhpoTRb7B6tru6Zeahf/EbdJDStGxRJALSup4iU0UoUtI0qGM2DZGaXv/uFytofhMjOc+/v5JCflvO9zju/75HT4znPec2hIkiQJAABS46xKHwAAAOUlAAEAUkYAAgCkjAAEAEgZAQgAkDICEAAgZQQgAEDKCEAAgJQRgAAAKSMAAQBSRgACAKSMAAQASBkBCACQMgIQACBlBCAAQMoIQACAlBGAAAApIwABAFJGAAIApIwABABIGQEIAJAyAhAAIGUEIABAyghAAICUEYAAACkjAAEAUkYAAgCkjAAEAEgZAQgAkDICEAAgZQQgAEDKCEAAgJQRgAAAKSMAAQBSRgAyKg8++GBMmzYtvvCFL8Ts2bPjl7/8ZaUPqeqtXLkyGhoaCm7ZbDa/P0mSWLlyZbS1tcWECRPisssui507dxY8x8DAQCxbtiwmT54cEydOjIULF8bevXvLfSpV4eWXX45rr7022traoqGhIZ555pmC/cWaz76+vuju7o5MJhOZTCa6u7vj4MGDJT67yjvd/N54441DXs9z584tGGN+T2716tXxla98JZqamqKlpSWuu+66ePfddwvGeA1TSgKQEXvqqaeip6cnvvvd78brr78ef/VXfxULFiyIDz74oNKHVvW+/OUvx759+/K3N954I7/v/vvvjzVr1sTatWvjtddei2w2G1dffXUcOnQoP6anpyfWr18f69ati82bN8fhw4ejq6srjh49WonTqahPPvkkZs2aFWvXrj3h/mLN5+LFi2PHjh2xYcOG2LBhQ+zYsSO6u7tLfn6Vdrr5jYiYP39+wev52WefLdhvfk9u06ZNcfvtt8e2bdti48aN8dlnn0VnZ2d88skn+TFew5RUAiP0F3/xF8mtt95asO1P//RPk7//+7+v0BHVhnvvvTeZNWvWCff94Q9/SLLZbPL9738/v+1//ud/kkwmk/zrv/5rkiRJcvDgwWTcuHHJunXr8mP++7//OznrrLOSDRs2lPTYq11EJOvXr8/fL9Z8vvXWW0lEJNu2bcuP2bp1axIRyTvvvFPis6oex89vkiTJkiVLkq9//esnfYz5HZn9+/cnEZFs2rQpSRKvYUrPCiAjMjg4GNu3b4/Ozs6C7Z2dnbFly5YKHVXt2LVrV7S1tcW0adPim9/8Zrz33nsREbF79+7o7e0tmNfGxsa49NJL8/O6ffv2+PTTTwvGtLW1xYwZM8z9cYo1n1u3bo1MJhNz5szJj5k7d25kMhlzHhEvvfRStLS0xIUXXhg33XRT7N+/P7/P/I5MLpeLiIjm5uaI8Bqm9AQgI/Lxxx/H0aNHo7W1tWB7a2tr9Pb2VuioasOcOXPi8ccfj//8z/+MH//4x9Hb2xvz5s2LAwcO5OfuVPPa29sb48ePj3POOeekY/hcseazt7c3Wlpahjx/S0tL6ud8wYIF8eSTT8YLL7wQ//zP/xyvvfZaXHHFFTEwMBAR5nckkiSJ5cuXx1e/+tWYMWNGRHgNU3pjK30A1KaGhoaC+0mSDNlGoQULFuT/PHPmzOjo6IgvfvGL8dhjj+Uvnh/NvJr7kyvGfJ5ovDmPWLRoUf7PM2bMiEsuuSSmTp0av/jFL+L6668/6ePM71B33HFH/PrXv47NmzcP2ec1TKlYAWREJk+eHGPGjBnym+P+/fuH/KbKqU2cODFmzpwZu3btyn8a+FTzms1mY3BwMPr6+k46hs8Vaz6z2Wx8+OGHQ57/o48+MufHmTJlSkydOjV27doVEeZ3uJYtWxY///nP48UXX4zzzz8/v91rmFITgIzI+PHjY/bs2bFx48aC7Rs3box58+ZV6Khq08DAQLz99tsxZcqUmDZtWmSz2YJ5HRwcjE2bNuXndfbs2TFu3LiCMfv27Ys333zT3B+nWPPZ0dERuVwuXn311fyYV155JXK5nDk/zoEDB2LPnj0xZcqUiDC/p5MkSdxxxx3x9NNPxwsvvBDTpk0r2O81TMlV5KMn1LR169Yl48aNSx555JHkrbfeSnp6epKJEycmv/3tbyt9aFVtxYoVyUsvvZS89957ybZt25Kurq6kqakpP2/f//73k0wmkzz99NPJG2+8kXzrW99KpkyZkvT39+ef49Zbb03OP//85Pnnn0/+67/+K7niiiuSWbNmJZ999lmlTqtiDh06lLz++uvJ66+/nkREsmbNmuT1119P3n///SRJijef8+fPT/7sz/4s2bp1a7J169Zk5syZSVdXV9nPt9xONb+HDh1KVqxYkWzZsiXZvXt38uKLLyYdHR3Jn/zJn5jfYfq7v/u7JJPJJC+99FKyb9++/O33v/99fozXMKUkABmVf/mXf0mmTp2ajB8/PvnzP//z/FcXcHKLFi1KpkyZkowbNy5pa2tLrr/++mTnzp35/X/4wx+Se++9N8lms0ljY2Pyta99LXnjjTcKnuPIkSPJHXfckTQ3NycTJkxIurq6kg8++KDcp1IVXnzxxSQihtyWLFmSJEnx5vPAgQPJt7/97aSpqSlpampKvv3tbyd9fX1lOsvKOdX8/v73v086OzuT8847Lxk3blxywQUXJEuWLBkyd+b35E40txGR/OQnP8mP8RqmlBqSJEnKveoIAEDluAYQACBlBCAAQMoIQACAlBGAAAApIwABAFJGAAIApIwABABIGQGYYg8++GBMmzYtvvCFL8Ts2bPjl7/85bAfOzAwECtXroyBgYESHmF6md/SMr+lZX5Ly/xSDL4IOqWeeuqp6O7ujgcffDD+8i//Mh5++OH4t3/7t3jrrbfiggsuOO3j+/v7I5PJRC6Xi0mTJpXhiNPF/JaW+S0t81ta5pdisAKYUmvWrImlS5fG3/7t38bFF18cP/zhD6O9vT0eeuihSh8aAFBiAjCFBgcHY/v27dHZ2VmwvbOzM7Zs2VKhowIAymVspQ+A8vv444/j6NGj0draWrC9tbU1ent7T/iYgYGBgutNDh48GBERuVyuZMeZZv39/QX/S3GZ39Iyv6VVrPlNkiQOHToUbW1tcdZZ1oPSRgCmWENDQ8H9JEmGbDtm9erVcd999w3ZPpzrBRm99vb2Sh9CXTO/pWV+S6tY87tnz544//zzi/Jc1A4BmEKTJ0+OMWPGDFnt279//5BVwWPuueeeWL58ef5+LpeLCy64IPbs2RPt7e1WAgFqTH9/f7S3t0dTU1OlD4UKEIApNH78+Jg9e3Zs3LgxvvGNb+S3b9y4Mb7+9a+f8DGNjY3R2Ng4ZPukSZPyK4c+UA5Qe072zg/1TQCm1PLly6O7uzsuueSS6OjoiB/96EfxwQcfxK233jqq5xOBAFA7BGBKLVq0KA4cOBD/9E//FPv27YsZM2bEs88+G1OnTh31c4pAAKgNvgiaUTnVF5GKQIDq5wul083nvim6U32aGACoPAFISYhAAKheApCSEYEAUJ0EICUlAgGg+ghASk4EAkB1EYCUhQgEgOohACkbEQgA1UEAUlYiEAAqTwBSdiIQACpLAFIRIhAAKkcAUjEiEAAqQwBSUSIQAMpPAFJxIhAAyksAUhVEIACUjwCkaohAACgPAUhVEYEAUHoCkKojAgGgtAQgVUkEAkDpCECqlggEgNIQgFQ1EQgAxScAqXoiEACKSwBSE0QgABSPAKRmiEAAKA4BSE0RgQBw5gQgNUcEAsCZEYDUJBEIAKMnAKlZIhAARkcAUtNEIACMnACk5olAABgZAUhdEIEAMHwCkLohAgFgeAQgdUUEAsDpCUDqjggEgFMTgNQlEQgAJycAqVsiEABOTABS10QgAAwlAKl7IhAACglAUkEEAsD/E4CkhggEgM8JQFJFBAKAACSFRCAAaScASSURCECaCUBSSwQCkFYCkFQTgQCkkQAk9UQgAGkjACFEIADpIgDh/4hAANJCAMIfEYEApIEAhOOIQADqnQCEExCBANQzAQgnIQIBqFcCEE5BBAJQjwQgnIYIBKDeCEAYBhEIQD0RgDBMIhCAeiEAYQREIAD1QADCCIlAAGqdAIRREIEA1DIBCKMkAgGoVQIQzoAIBKAWCUA4QyIQgFojAKEIRCAAtUQAQpGIQABqhQCEIhKBANQCAQhFJgIBqHYCEEpABAJQzQQglIgIBKBaCUAoIREIQDUSgFBiIhCAaiMAoQxEIADVRABCmYhAAKqFAIQyEoEAVAMBCGUmAgGoNAEIFSACAagkAQgVIgIBqBQBCBUkAgGoBAEIFSYCASg3AVhjVq5cGQ0NDQW3bDab358kSaxcuTLa2tpiwoQJcdlll8XOnTsLnmNgYCCWLVsWkydPjokTJ8bChQtj79695T4V/ogIBKCcBGAN+vKXvxz79u3L39544438vvvvvz/WrFkTa9eujddeey2y2WxcffXVcejQofyYnp6eWL9+faxbty42b94chw8fjq6urjh69GglTof/IwIBKJexlT4ARm7s2LEFq37HJEkSP/zhD+O73/1uXH/99RER8dhjj0Vra2v87Gc/i1tuuSVyuVw88sgj8cQTT8RVV10VERE//elPo729PZ5//vm45pprynouFDoWgUmSVPpQAKhjVgBr0K5du6KtrS2mTZsW3/zmN+O9996LiIjdu3dHb29vdHZ25sc2NjbGpZdeGlu2bImIiO3bt8enn35aMKatrS1mzJiRH0NlWQkEoNSsANaYOXPmxOOPPx4XXnhhfPjhh/G9730v5s2bFzt37oze3t6IiGhtbS14TGtra7z//vsREdHb2xvjx4+Pc845Z8iYY48/kYGBgRgYGMjf7+/vL9YpcQJWAgEoJQFYYxYsWJD/88yZM6OjoyO++MUvxmOPPRZz586NiBiyejScFaXTjVm9enXcd999Z3DkjJQIBKBUvAVc4yZOnBgzZ86MXbt25a8LPH4lb//+/flVwWw2G4ODg9HX13fSMSdyzz33RC6Xy9/27NlT5DPhRLwdDEApCMAaNzAwEG+//XZMmTIlpk2bFtlsNjZu3JjfPzg4GJs2bYp58+ZFRMTs2bNj3LhxBWP27dsXb775Zn7MiTQ2NsakSZMKbpSHCASg2LwFXGPuuuuuuPbaa+OCCy6I/fv3x/e+973o7++PJUuWRENDQ/T09MSqVati+vTpMX369Fi1alWcffbZsXjx4oiIyGQysXTp0lixYkWce+650dzcHHfddVfMnDkz/6lgqo+3gwEoJgFYY/bu3Rvf+ta34uOPP47zzjsv5s6dG9u2bYupU6dGRMTdd98dR44cidtuuy36+vpizpw58dxzz0VTU1P+OR544IEYO3Zs3HDDDXHkyJG48sor49FHH40xY8ZU6rQYBhEIQLE0JP42YRT6+/sjk8lELpfzdnCZiUCgGPwcTzfXAEKNcU0gAGdKAEINEoEAnAkBCDVKBAIwWgIQapgIBGA0BCDUOBEIwEgJQKgDIhCAkRCAUCdEIADDJQChjohAAIZDAEKdEYEAnI4AhDokAgE4FQEIdUoEAnAyAhDqmAgE4EQEINQ5EQjA8QQgpIAIBOCPCUBICREIwDECEFJEBAIQIQAhdUQgAAIQUkgEAqSbAISUEoEA6SUAIcVEIEA6CUBIOREIkD4CEBCBACkjAIGIEIEAaSIAgTwRCJAOAhAoIAIB6p8ABIYQgQD1TQACJyQCAeqXAAROSgQC1CcBCJySCASoPwIQOC0RCFBfBCAwLCIQoH4IQGDYRCBAfRCAwIiIQIDaJwCBEROBALVNAAKjIgIBapcABEZNBALUJgEInBERCFB7BCBwxkQgQG0RgEBRiECA2iEAgaIRgQC1QQACRSUCAaqfAASKTgQCVDcBCJSECASoXgIQKBkRCFCdBCBQUiIQoPoIQKDkRCBAdRGAQFmIQIDqIQCBshGBANVBAAJlJQIBKk8AAmUnAgEqSwACFSECASpHAAIVIwIBKkMAAhUlAgHKTwACFScCAcpLAAJVQQQClI8ABKqGCAQoDwEIVBURCFB6AhCoOiIQoLQEIFCVRCBA6QhAoGqJQIDSEIBAVROBAMUnAIGqJwIBiksAAjVBBAIUjwAEaoYIBCgOAQjUFBEIcOYEIFBzRCDAmRGAQE0SgQCjJwCBmiUCAUZHAAI1TQQCjJwABGqeCAQYGQEI1AURCDB8AhCoGyIQYHgEIFBXRCDA6QlAoO6IQIBTE4BAXRKBACcnAIG6JQIBTkwAAnVNBAIMJQCrzMsvvxzXXntttLW1RUNDQzzzzDMF+5MkiZUrV0ZbW1tMmDAhLrvssti5c2fBmIGBgVi2bFlMnjw5Jk6cGAsXLoy9e/cWjOnr64vu7u7IZDKRyWSiu7s7Dh48WOKzg8oQgQCFBGCV+eSTT2LWrFmxdu3aE+6///77Y82aNbF27dp47bXXIpvNxtVXXx2HDh3Kj+np6Yn169fHunXrYvPmzXH48OHo6uqKo0eP5scsXrw4duzYERs2bIgNGzbEjh07oru7u+TnB5UiAgH+X0OSJEmlD4ITa2hoiPXr18d1110XEZ//BdbW1hY9PT3xne98JyI+X+1rbW2NH/zgB3HLLbdELpeL8847L5544olYtGhRRET87ne/i/b29nj22Wfjmmuuibfffju+9KUvxbZt22LOnDkREbFt27bo6OiId955Jy666KLTHlt/f39kMpnI5XIxadKk0kwAlEBDQ0P4sQd+jqedFcAasnv37ujt7Y3Ozs78tsbGxrj00ktjy5YtERGxffv2+PTTTwvGtLW1xYwZM/Jjtm7dGplMJh9/ERFz586NTCaTHwP1ykogQMTYSh8Aw9fb2xsREa2trQXbW1tb4/3338+PGT9+fJxzzjlDxhx7fG9vb7S0tAx5/paWlvyY4w0MDMTAwED+fn9//+hPBCrsWARaCQTSygpgDTp+9WI4KxrHjznR+FM9z+rVq/MfGMlkMtHe3j6KI4fqYSUQSDMBWEOy2WxExJBVuv379+dXBbPZbAwODkZfX98px3z44YdDnv+jjz4asrp4zD333BO5XC5/27NnzxmfD1SaCATSSgDWkGnTpkU2m42NGzfmtw0ODsamTZti3rx5ERExe/bsGDduXMGYffv2xZtvvpkf09HREblcLl599dX8mFdeeSVyuVx+zPEaGxtj0qRJBTeoByIQSCPXAFaZw4cPx29+85v8/d27d8eOHTuiubk5Lrjggujp6YlVq1bF9OnTY/r06bFq1ao4++yzY/HixRERkclkYunSpbFixYo499xzo7m5Oe66666YOXNmXHXVVRERcfHFF8f8+fPjpptuiocffjgiIm6++ebo6uoa1ieAod64JhBIGwFYZX71q1/F5Zdfnr+/fPnyiIhYsmRJPProo3H33XfHkSNH4rbbbou+vr6YM2dOPPfcc9HU1JR/zAMPPBBjx46NG264IY4cORJXXnllPProozFmzJj8mCeffDLuvPPO/KeFFy5ceNLvHoQ0EIFAmvgeQEbF90dRr0QgaeHneLq5BhDgj7gmEEgDAQhwHBEI1DsBCHACIhCoZwIQ4CREIFCvBCDAKYhAoB4JQIDTEIFAvRGAAMMgAoF6IgABhkkEAvVCAAKMgAgE6oEABBghEQjUOgEIMAoiEKhlAhBglEQgUKsEIMAZEIFALRKAAGdIBAK1RgACFIEIBGqJAAQoEhEI1AoBCFBEIhCoBQIQoMhEIFDtBCBACYhAoJoJQIASEYFAtRKAACUkAoFqJAABSkwEAtVGAAKUgQgEqokABCgTEQhUCwEIUEYiEKgGAhCgzEQgUGkCEKACRCBQSQIQoEJEIFApAhCggkQgUAkCEKDCRCBQbgIQoAqIQKCcBCBAlRCBQLkIQIAqIgKBchCAAFVGBAKlJgABqpAIBEpJAAJUKREIlIoABKhiIhAoBQEIUOVEIFBsAhCgBohAoJgEIECNEIFAsQhAgBoiAoFiEIAANUYEAmdKAALUIBEInAkBCFCjRCAwWgIQoIaJQGA0BCBAjROBwEgJQIA6IAKBkRCAAHVCBALDJQAB6ogIBIZDAALUGREInI4ABKhDIhA4FQEIUKdEIHAyAhCgjolA4EQEIECdE4HA8QQgQAqIQOCPCUCAlBCBwDECECBFRCAQIQABUkcEAgIQIIVEIKSbAARIKREI6SUAAVJMBEI6CUCAlBOBkD4CEAARCCkjAAGICBEIaSIAAcgTgZAOAhCAAiIQ6p8ABGAIEQj1TQACcEIiEOqXAATgpEQg1CcBCMApiUCoPwIQgNMSgVBfBCAAwyICoX4IQACGTQRCfRCAAIyICITaJwABGDERCLVNAAIwKiIQapcABGDURCDUJgFYZV5++eW49tpro62tLRoaGuKZZ54p2H/jjTdGQ0NDwW3u3LkFYwYGBmLZsmUxefLkmDhxYixcuDD27t1bMKavry+6u7sjk8lEJpOJ7u7uOHjwYInPDqhHIhBqjwCsMp988knMmjUr1q5de9Ix8+fPj3379uVvzz77bMH+np6eWL9+faxbty42b94chw8fjq6urjh69Gh+zOLFi2PHjh2xYcOG2LBhQ+zYsSO6u7tLdl5AfROBUFvGVvoAKLRgwYJYsGDBKcc0NjZGNps94b5cLhePPPJIPPHEE3HVVVdFRMRPf/rTaG9vj+effz6uueaaePvtt2PDhg2xbdu2mDNnTkRE/PjHP46Ojo54991346KLLiruSQGpcCwCkySp9KEAp2EFsAa99NJL0dLSEhdeeGHcdNNNsX///vy+7du3x6effhqdnZ35bW1tbTFjxozYsmVLRERs3bo1MplMPv4iIubOnRuZTCY/BmA0rARCbbACWGMWLFgQf/M3fxNTp06N3bt3xz/+4z/GFVdcEdu3b4/Gxsbo7e2N8ePHxznnnFPwuNbW1ujt7Y2IiN7e3mhpaRny3C0tLfkxxxsYGIiBgYH8/f7+/iKeFVBPrARC9ROANWbRokX5P8+YMSMuueSSmDp1avziF7+I66+//qSPO/638hP9hn6q39xXr14d99133xkcOZAmIhCqm7eAa9yUKVNi6tSpsWvXroiIyGazMTg4GH19fQXj9u/fH62trfkxH3744ZDn+uijj/JjjnfPPfdELpfL3/bs2VPkMwHqjbeDoXoJwBp34MCB2LNnT0yZMiUiImbPnh3jxo2LjRs35sfs27cv3nzzzZg3b15ERHR0dEQul4tXX301P+aVV16JXC6XH3O8xsbGmDRpUsEN4HREIFQnbwFXmcOHD8dvfvOb/P3du3fHjh07orm5OZqbm2PlypXx13/91zFlypT47W9/G//wD/8QkydPjm984xsREZHJZGLp0qWxYsWKOPfcc6O5uTnuuuuumDlzZv5TwRdffHHMnz8/brrppnj44YcjIuLmm2+Orq4unwAGis7bwVB9BGCV+dWvfhWXX355/v7y5csjImLJkiXx0EMPxRtvvBGPP/54HDx4MKZMmRKXX355PPXUU9HU1JR/zAMPPBBjx46NG264IY4cORJXXnllPProozFmzJj8mCeffDLuvPPO/KeFFy5ceMrvHgQ4EyIQqktD4v+NjEJ/f39kMpnI5XLeDgaGTQRWDz/H0801gACUjWsCoToIQADKSgRC5QlAAMpOBEJlCUAAKkIEQuUIQAAqRgRCZQhAACpKBEL5CUAAKk4EQnkJQACqggiE8hGAAFQNEQjlIQABqCoiEEpPAAJQdUQglJYABKAqiUAoHQEIQNUSgVAaAhCAqiYCofgEIABVTwRCcQlAAGqCCITiEYAA1AwRCMUhAAGoKSIQzpwABKDmiEA4MwIQgJokAmH0BCAANUsEwugIQABqmgiEkROAANQ8EQgjIwABqAsiEIZPAAJQN0QgDI8ABKCuiEA4PQEIQN0RgXBqAhCAuiQC4eQEIAB1SwTCiQlAAOqaCIShBCAAdU8EQiEBCEAqiED4fwIQgNQQgfA5AQhAqohAEIAApJAIJO0EIACpJAJJMwEIQGqJQNJKAAKQaiKQNBKAAKSeCCRtBCAAhAgkXQQgAPwfEUhaCEAA+CMikDQQgABwHBFIvROAAHACIpB6JgAB4CREIPVKAALAKYhA6pEABIDTEIHUGwEIAMMgAqknAhAAhkkEUi8EIACMgAikHghAABghEUitE4AAMAoikFomAAFglEQgtUoAAsAZEIHUIgEIAGdIBFJrBCAAFIEIpJYIQAAoEhFIrRCAAFBEIpBaIAABoMhEINVOAAJACYhAqpkABIASEYFUKwEIACUkAqlGAhAASkwEUm0EIACUgQikmghAACgTEUi1EIAAUEYikGogAAGgzEQglSYAAaACRCCVJAABoEJEIJUiAAGggkQglSAAAaDCRCDlJgABoAqIQMpJAAJAlRCBlIsABIAqIgIpBwEIAFVGBFJqAhAAqpAIpJQEIABUKRFIqQjAKrJ69er4yle+Ek1NTdHS0hLXXXddvPvuuwVjkiSJlStXRltbW0yYMCEuu+yy2LlzZ8GYgYGBWLZsWUyePDkmTpwYCxcujL179xaM6evri+7u7shkMpHJZKK7uzsOHjxY6lMEYIREIKUgAKvIpk2b4vbbb49t27bFxo0b47PPPovOzs745JNP8mPuv//+WLNmTaxduzZee+21yGazcfXVV8ehQ4fyY3p6emL9+vWxbt262Lx5cxw+fDi6urri6NGj+TGLFy+OHTt2xIYNG2LDhg2xY8eO6O7uLuv5AjA8IpBia0iSJKn0QXBiH330UbS0tMSmTZvia1/7WiRJEm1tbdHT0xPf+c53IuLz1b7W1tb4wQ9+ELfcckvkcrk477zz4oknnohFixZFRMTvfve7aG9vj2effTauueaaePvtt+NLX/pSbNu2LebMmRMREdu2bYuOjo5455134qKLLjrtsfX390cmk4lcLheTJk0q3SQAkNfQ0BDF+mvbz/F0swJYxXK5XERENDc3R0TE7t27o7e3Nzo7O/NjGhsb49JLL40tW7ZERMT27dvj008/LRjT1tYWM2bMyI/ZunVrZDKZfPxFRMydOzcymUx+DADVx0ogxTK20gfAiSVJEsuXL4+vfvWrMWPGjIiI6O3tjYiI1tbWgrGtra3x/vvv58eMHz8+zjnnnCFjjj2+t7c3Wlpahvw3W1pa8mOONzAwEAMDA/n7/f39ozwzAM7EsQj0Bh5nwgpglbrjjjvi17/+dfz7v//7kH3H//Y3nN8Ijx9zovGnep7Vq1fnPzCSyWSivb19OKcBQAlYCeRMCcAqtGzZsvj5z38eL774Ypx//vn57dlsNiJiyCrd/v3786uC2Ww2BgcHo6+v75RjPvzwwyH/3Y8++mjI6uIx99xzT+Ryufxtz549oz9BAM6YCORMCMAqkiRJ3HHHHfH000/HCy+8ENOmTSvYP23atMhms7Fx48b8tsHBwdi0aVPMmzcvIiJmz54d48aNKxizb9++ePPNN/NjOjo6IpfLxauvvpof88orr0Qul8uPOV5jY2NMmjSp4AZAZYlARss1gFXk9ttvj5/97GfxH//xH9HU1JRf6ctkMjFhwoRoaGiInp6eWLVqVUyfPj2mT58eq1atirPPPjsWL16cH7t06dJYsWJFnHvuudHc3Bx33XVXzJw5M6666qqIiLj44otj/vz5cdNNN8XDDz8cERE333xzdHV1DesTwABUD9cEMhoCsIo89NBDERFx2WWXFWz/yU9+EjfeeGNERNx9991x5MiRuO2226Kvry/mzJkTzz33XDQ1NeXHP/DAAzF27Ni44YYb4siRI3HllVfGo48+GmPGjMmPefLJJ+POO+/Mf1p44cKFsXbt2tKeIAAlIQIZKd8DyKj4/iiA6jOSCPRzPN1cAwgAdcI1gQyXAASAOiICGQ4BCAB1RgRyOgIQAOqQCORUBCAA1CkRyMkIQACoYyKQExGAAFDnRCDHE4AAkAIikD8mAAEgJUQgxwhAAEgREUiEfwuYUTr2Tw319/dX+EgAGKlcLheZTCYiwr8fnFICkFE5cOBARES0t7dX+EgAOBOHDh3KxyDpIQAZlebm5oiI+OCDD/zgKIH+/v5ob2+PPXv2+EfaS8D8lpb5La1izW+SJHHo0KFoa2sr4tFRKwQgo3LWWZ9fPprJZPyAL6FJkyaZ3xIyv6VlfkurGPPrF/j08iEQAICUEYAAACkjABmVxsbGuPfee6OxsbHSh1KXzG9pmd/SMr+lZX4phobE578BAFLFCiAAQMoIQACAlBGAAAApIwABAFJGAAIApIwABABIGQEIAJAyAhAAIGUEIABAyghAAICUEYAAACkjAAEAUkYAAgCkjAAEAEgZAQgAkDICEAAgZQQgAEDKCEAAgJQRgAAAKSMAAQBSRgACAKSMAAQASBkBCACQMgIQACBlBCAAQMoIQACAlBGAAAApIwABAFJGAAIApIwABABIGQEIAJAyAhAAIGUEIABAyghAAICUEYAAACkjAAEAUkYAAgCkjAAEAEgZAQgAkDICEAAgZQQgAEDKCEAAgJQRgAAAKSMAAQBSRgACAKSMAAQASJn/BXJL0WQ6w21rAAAAAElFTkSuQmCC",
      "text/html": [
       "\n",
       "            <div style=\"display: inline-block;\">\n",
       "                <div class=\"jupyter-widgets widget-label\" style=\"text-align: center;\">\n",
       "                    Figure\n",
       "                </div>\n",
       "                <img src='data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkxklEQVR4nO3de4yV5Z3A8d/IZYpkOHHEmcOsIyENurawJItdGLZb7yMkI7VuVlqaCSasl1UxEzB23WYjbhpoTRb7B6tru6Zeahf/EbdJDStGxRJALSup4iU0UoUtI0qGM2DZGaXv/uFytofhMjOc+/v5JCflvO9zju/75HT4znPec2hIkiQJAABS46xKHwAAAOUlAAEAUkYAAgCkjAAEAEgZAQgAkDICEAAgZQQgAEDKCEAAgJQRgAAAKSMAAQBSRgACAKSMAAQASBkBCACQMgIQACBlBCAAQMoIQACAlBGAAAApIwABAFJGAAIApIwABABIGQEIAJAyAhAAIGUEIABAyghAAICUEYAAACkjAAEAUkYAAgCkjAAEAEgZAQgAkDICEAAgZQQgAEDKCEAAgJQRgAAAKSMAAQBSRgAyKg8++GBMmzYtvvCFL8Ts2bPjl7/8ZaUPqeqtXLkyGhoaCm7ZbDa/P0mSWLlyZbS1tcWECRPisssui507dxY8x8DAQCxbtiwmT54cEydOjIULF8bevXvLfSpV4eWXX45rr7022traoqGhIZ555pmC/cWaz76+vuju7o5MJhOZTCa6u7vj4MGDJT67yjvd/N54441DXs9z584tGGN+T2716tXxla98JZqamqKlpSWuu+66ePfddwvGeA1TSgKQEXvqqaeip6cnvvvd78brr78ef/VXfxULFiyIDz74oNKHVvW+/OUvx759+/K3N954I7/v/vvvjzVr1sTatWvjtddei2w2G1dffXUcOnQoP6anpyfWr18f69ati82bN8fhw4ejq6srjh49WonTqahPPvkkZs2aFWvXrj3h/mLN5+LFi2PHjh2xYcOG2LBhQ+zYsSO6u7tLfn6Vdrr5jYiYP39+wev52WefLdhvfk9u06ZNcfvtt8e2bdti48aN8dlnn0VnZ2d88skn+TFew5RUAiP0F3/xF8mtt95asO1P//RPk7//+7+v0BHVhnvvvTeZNWvWCff94Q9/SLLZbPL9738/v+1//ud/kkwmk/zrv/5rkiRJcvDgwWTcuHHJunXr8mP++7//OznrrLOSDRs2lPTYq11EJOvXr8/fL9Z8vvXWW0lEJNu2bcuP2bp1axIRyTvvvFPis6oex89vkiTJkiVLkq9//esnfYz5HZn9+/cnEZFs2rQpSRKvYUrPCiAjMjg4GNu3b4/Ozs6C7Z2dnbFly5YKHVXt2LVrV7S1tcW0adPim9/8Zrz33nsREbF79+7o7e0tmNfGxsa49NJL8/O6ffv2+PTTTwvGtLW1xYwZM8z9cYo1n1u3bo1MJhNz5szJj5k7d25kMhlzHhEvvfRStLS0xIUXXhg33XRT7N+/P7/P/I5MLpeLiIjm5uaI8Bqm9AQgI/Lxxx/H0aNHo7W1tWB7a2tr9Pb2VuioasOcOXPi8ccfj//8z/+MH//4x9Hb2xvz5s2LAwcO5OfuVPPa29sb48ePj3POOeekY/hcseazt7c3Wlpahjx/S0tL6ud8wYIF8eSTT8YLL7wQ//zP/xyvvfZaXHHFFTEwMBAR5nckkiSJ5cuXx1e/+tWYMWNGRHgNU3pjK30A1KaGhoaC+0mSDNlGoQULFuT/PHPmzOjo6IgvfvGL8dhjj+Uvnh/NvJr7kyvGfJ5ovDmPWLRoUf7PM2bMiEsuuSSmTp0av/jFL+L6668/6ePM71B33HFH/PrXv47NmzcP2ec1TKlYAWREJk+eHGPGjBnym+P+/fuH/KbKqU2cODFmzpwZu3btyn8a+FTzms1mY3BwMPr6+k46hs8Vaz6z2Wx8+OGHQ57/o48+MufHmTJlSkydOjV27doVEeZ3uJYtWxY///nP48UXX4zzzz8/v91rmFITgIzI+PHjY/bs2bFx48aC7Rs3box58+ZV6Khq08DAQLz99tsxZcqUmDZtWmSz2YJ5HRwcjE2bNuXndfbs2TFu3LiCMfv27Ys333zT3B+nWPPZ0dERuVwuXn311fyYV155JXK5nDk/zoEDB2LPnj0xZcqUiDC/p5MkSdxxxx3x9NNPxwsvvBDTpk0r2O81TMlV5KMn1LR169Yl48aNSx555JHkrbfeSnp6epKJEycmv/3tbyt9aFVtxYoVyUsvvZS89957ybZt25Kurq6kqakpP2/f//73k0wmkzz99NPJG2+8kXzrW99KpkyZkvT39+ef49Zbb03OP//85Pnnn0/+67/+K7niiiuSWbNmJZ999lmlTqtiDh06lLz++uvJ66+/nkREsmbNmuT1119P3n///SRJijef8+fPT/7sz/4s2bp1a7J169Zk5syZSVdXV9nPt9xONb+HDh1KVqxYkWzZsiXZvXt38uKLLyYdHR3Jn/zJn5jfYfq7v/u7JJPJJC+99FKyb9++/O33v/99fozXMKUkABmVf/mXf0mmTp2ajB8/PvnzP//z/FcXcHKLFi1KpkyZkowbNy5pa2tLrr/++mTnzp35/X/4wx+Se++9N8lms0ljY2Pyta99LXnjjTcKnuPIkSPJHXfckTQ3NycTJkxIurq6kg8++KDcp1IVXnzxxSQihtyWLFmSJEnx5vPAgQPJt7/97aSpqSlpampKvv3tbyd9fX1lOsvKOdX8/v73v086OzuT8847Lxk3blxywQUXJEuWLBkyd+b35E40txGR/OQnP8mP8RqmlBqSJEnKveoIAEDluAYQACBlBCAAQMoIQACAlBGAAAApIwABAFJGAAIApIwABABIGQGYYg8++GBMmzYtvvCFL8Ts2bPjl7/85bAfOzAwECtXroyBgYESHmF6md/SMr+lZX5Ly/xSDL4IOqWeeuqp6O7ujgcffDD+8i//Mh5++OH4t3/7t3jrrbfiggsuOO3j+/v7I5PJRC6Xi0mTJpXhiNPF/JaW+S0t81ta5pdisAKYUmvWrImlS5fG3/7t38bFF18cP/zhD6O9vT0eeuihSh8aAFBiAjCFBgcHY/v27dHZ2VmwvbOzM7Zs2VKhowIAymVspQ+A8vv444/j6NGj0draWrC9tbU1ent7T/iYgYGBgutNDh48GBERuVyuZMeZZv39/QX/S3GZ39Iyv6VVrPlNkiQOHToUbW1tcdZZ1oPSRgCmWENDQ8H9JEmGbDtm9erVcd999w3ZPpzrBRm99vb2Sh9CXTO/pWV+S6tY87tnz544//zzi/Jc1A4BmEKTJ0+OMWPGDFnt279//5BVwWPuueeeWL58ef5+LpeLCy64IPbs2RPt7e1WAgFqTH9/f7S3t0dTU1OlD4UKEIApNH78+Jg9e3Zs3LgxvvGNb+S3b9y4Mb7+9a+f8DGNjY3R2Ng4ZPukSZPyK4c+UA5Qe072zg/1TQCm1PLly6O7uzsuueSS6OjoiB/96EfxwQcfxK233jqq5xOBAFA7BGBKLVq0KA4cOBD/9E//FPv27YsZM2bEs88+G1OnTh31c4pAAKgNvgiaUTnVF5GKQIDq5wul083nvim6U32aGACoPAFISYhAAKheApCSEYEAUJ0EICUlAgGg+ghASk4EAkB1EYCUhQgEgOohACkbEQgA1UEAUlYiEAAqTwBSdiIQACpLAFIRIhAAKkcAUjEiEAAqQwBSUSIQAMpPAFJxIhAAyksAUhVEIACUjwCkaohAACgPAUhVEYEAUHoCkKojAgGgtAQgVUkEAkDpCECqlggEgNIQgFQ1EQgAxScAqXoiEACKSwBSE0QgABSPAKRmiEAAKA4BSE0RgQBw5gQgNUcEAsCZEYDUJBEIAKMnAKlZIhAARkcAUtNEIACMnACk5olAABgZAUhdEIEAMHwCkLohAgFgeAQgdUUEAsDpCUDqjggEgFMTgNQlEQgAJycAqVsiEABOTABS10QgAAwlAKl7IhAACglAUkEEAsD/E4CkhggEgM8JQFJFBAKAACSFRCAAaScASSURCECaCUBSSwQCkFYCkFQTgQCkkQAk9UQgAGkjACFEIADpIgDh/4hAANJCAMIfEYEApIEAhOOIQADqnQCEExCBANQzAQgnIQIBqFcCEE5BBAJQjwQgnIYIBKDeCEAYBhEIQD0RgDBMIhCAeiEAYQREIAD1QADCCIlAAGqdAIRREIEA1DIBCKMkAgGoVQIQzoAIBKAWCUA4QyIQgFojAKEIRCAAtUQAQpGIQABqhQCEIhKBANQCAQhFJgIBqHYCEEpABAJQzQQglIgIBKBaCUAoIREIQDUSgFBiIhCAaiMAoQxEIADVRABCmYhAAKqFAIQyEoEAVAMBCGUmAgGoNAEIFSACAagkAQgVIgIBqBQBCBUkAgGoBAEIFSYCASg3AVhjVq5cGQ0NDQW3bDab358kSaxcuTLa2tpiwoQJcdlll8XOnTsLnmNgYCCWLVsWkydPjokTJ8bChQtj79695T4V/ogIBKCcBGAN+vKXvxz79u3L39544438vvvvvz/WrFkTa9eujddeey2y2WxcffXVcejQofyYnp6eWL9+faxbty42b94chw8fjq6urjh69GglTof/IwIBKJexlT4ARm7s2LEFq37HJEkSP/zhD+O73/1uXH/99RER8dhjj0Vra2v87Gc/i1tuuSVyuVw88sgj8cQTT8RVV10VERE//elPo729PZ5//vm45pprynouFDoWgUmSVPpQAKhjVgBr0K5du6KtrS2mTZsW3/zmN+O9996LiIjdu3dHb29vdHZ25sc2NjbGpZdeGlu2bImIiO3bt8enn35aMKatrS1mzJiRH0NlWQkEoNSsANaYOXPmxOOPPx4XXnhhfPjhh/G9730v5s2bFzt37oze3t6IiGhtbS14TGtra7z//vsREdHb2xvjx4+Pc845Z8iYY48/kYGBgRgYGMjf7+/vL9YpcQJWAgEoJQFYYxYsWJD/88yZM6OjoyO++MUvxmOPPRZz586NiBiyejScFaXTjVm9enXcd999Z3DkjJQIBKBUvAVc4yZOnBgzZ86MXbt25a8LPH4lb//+/flVwWw2G4ODg9HX13fSMSdyzz33RC6Xy9/27NlT5DPhRLwdDEApCMAaNzAwEG+//XZMmTIlpk2bFtlsNjZu3JjfPzg4GJs2bYp58+ZFRMTs2bNj3LhxBWP27dsXb775Zn7MiTQ2NsakSZMKbpSHCASg2LwFXGPuuuuuuPbaa+OCCy6I/fv3x/e+973o7++PJUuWRENDQ/T09MSqVati+vTpMX369Fi1alWcffbZsXjx4oiIyGQysXTp0lixYkWce+650dzcHHfddVfMnDkz/6lgqo+3gwEoJgFYY/bu3Rvf+ta34uOPP47zzjsv5s6dG9u2bYupU6dGRMTdd98dR44cidtuuy36+vpizpw58dxzz0VTU1P+OR544IEYO3Zs3HDDDXHkyJG48sor49FHH40xY8ZU6rQYBhEIQLE0JP42YRT6+/sjk8lELpfzdnCZiUCgGPwcTzfXAEKNcU0gAGdKAEINEoEAnAkBCDVKBAIwWgIQapgIBGA0BCDUOBEIwEgJQKgDIhCAkRCAUCdEIADDJQChjohAAIZDAEKdEYEAnI4AhDokAgE4FQEIdUoEAnAyAhDqmAgE4EQEINQ5EQjA8QQgpIAIBOCPCUBICREIwDECEFJEBAIQIQAhdUQgAAIQUkgEAqSbAISUEoEA6SUAIcVEIEA6CUBIOREIkD4CEBCBACkjAIGIEIEAaSIAgTwRCJAOAhAoIAIB6p8ABIYQgQD1TQACJyQCAeqXAAROSgQC1CcBCJySCASoPwIQOC0RCFBfBCAwLCIQoH4IQGDYRCBAfRCAwIiIQIDaJwCBEROBALVNAAKjIgIBapcABEZNBALUJgEInBERCFB7BCBwxkQgQG0RgEBRiECA2iEAgaIRgQC1QQACRSUCAaqfAASKTgQCVDcBCJSECASoXgIQKBkRCFCdBCBQUiIQoPoIQKDkRCBAdRGAQFmIQIDqIQCBshGBANVBAAJlJQIBKk8AAmUnAgEqSwACFSECASpHAAIVIwIBKkMAAhUlAgHKTwACFScCAcpLAAJVQQQClI8ABKqGCAQoDwEIVBURCFB6AhCoOiIQoLQEIFCVRCBA6QhAoGqJQIDSEIBAVROBAMUnAIGqJwIBiksAAjVBBAIUjwAEaoYIBCgOAQjUFBEIcOYEIFBzRCDAmRGAQE0SgQCjJwCBmiUCAUZHAAI1TQQCjJwABGqeCAQYGQEI1AURCDB8AhCoGyIQYHgEIFBXRCDA6QlAoO6IQIBTE4BAXRKBACcnAIG6JQIBTkwAAnVNBAIMJQCrzMsvvxzXXntttLW1RUNDQzzzzDMF+5MkiZUrV0ZbW1tMmDAhLrvssti5c2fBmIGBgVi2bFlMnjw5Jk6cGAsXLoy9e/cWjOnr64vu7u7IZDKRyWSiu7s7Dh48WOKzg8oQgQCFBGCV+eSTT2LWrFmxdu3aE+6///77Y82aNbF27dp47bXXIpvNxtVXXx2HDh3Kj+np6Yn169fHunXrYvPmzXH48OHo6uqKo0eP5scsXrw4duzYERs2bIgNGzbEjh07oru7u+TnB5UiAgH+X0OSJEmlD4ITa2hoiPXr18d1110XEZ//BdbW1hY9PT3xne98JyI+X+1rbW2NH/zgB3HLLbdELpeL8847L5544olYtGhRRET87ne/i/b29nj22Wfjmmuuibfffju+9KUvxbZt22LOnDkREbFt27bo6OiId955Jy666KLTHlt/f39kMpnI5XIxadKk0kwAlEBDQ0P4sQd+jqedFcAasnv37ujt7Y3Ozs78tsbGxrj00ktjy5YtERGxffv2+PTTTwvGtLW1xYwZM/Jjtm7dGplMJh9/ERFz586NTCaTHwP1ykogQMTYSh8Aw9fb2xsREa2trQXbW1tb4/3338+PGT9+fJxzzjlDxhx7fG9vb7S0tAx5/paWlvyY4w0MDMTAwED+fn9//+hPBCrsWARaCQTSygpgDTp+9WI4KxrHjznR+FM9z+rVq/MfGMlkMtHe3j6KI4fqYSUQSDMBWEOy2WxExJBVuv379+dXBbPZbAwODkZfX98px3z44YdDnv+jjz4asrp4zD333BO5XC5/27NnzxmfD1SaCATSSgDWkGnTpkU2m42NGzfmtw0ODsamTZti3rx5ERExe/bsGDduXMGYffv2xZtvvpkf09HREblcLl599dX8mFdeeSVyuVx+zPEaGxtj0qRJBTeoByIQSCPXAFaZw4cPx29+85v8/d27d8eOHTuiubk5Lrjggujp6YlVq1bF9OnTY/r06bFq1ao4++yzY/HixRERkclkYunSpbFixYo499xzo7m5Oe66666YOXNmXHXVVRERcfHFF8f8+fPjpptuiocffjgiIm6++ebo6uoa1ieAod64JhBIGwFYZX71q1/F5Zdfnr+/fPnyiIhYsmRJPProo3H33XfHkSNH4rbbbou+vr6YM2dOPPfcc9HU1JR/zAMPPBBjx46NG264IY4cORJXXnllPProozFmzJj8mCeffDLuvPPO/KeFFy5ceNLvHoQ0EIFAmvgeQEbF90dRr0QgaeHneLq5BhDgj7gmEEgDAQhwHBEI1DsBCHACIhCoZwIQ4CREIFCvBCDAKYhAoB4JQIDTEIFAvRGAAMMgAoF6IgABhkkEAvVCAAKMgAgE6oEABBghEQjUOgEIMAoiEKhlAhBglEQgUKsEIMAZEIFALRKAAGdIBAK1RgACFIEIBGqJAAQoEhEI1AoBCFBEIhCoBQIQoMhEIFDtBCBACYhAoJoJQIASEYFAtRKAACUkAoFqJAABSkwEAtVGAAKUgQgEqokABCgTEQhUCwEIUEYiEKgGAhCgzEQgUGkCEKACRCBQSQIQoEJEIFApAhCggkQgUAkCEKDCRCBQbgIQoAqIQKCcBCBAlRCBQLkIQIAqIgKBchCAAFVGBAKlJgABqpAIBEpJAAJUKREIlIoABKhiIhAoBQEIUOVEIFBsAhCgBohAoJgEIECNEIFAsQhAgBoiAoFiEIAANUYEAmdKAALUIBEInAkBCFCjRCAwWgIQoIaJQGA0BCBAjROBwEgJQIA6IAKBkRCAAHVCBALDJQAB6ogIBIZDAALUGREInI4ABKhDIhA4FQEIUKdEIHAyAhCgjolA4EQEIECdE4HA8QQgQAqIQOCPCUCAlBCBwDECECBFRCAQIQABUkcEAgIQIIVEIKSbAARIKREI6SUAAVJMBEI6CUCAlBOBkD4CEAARCCkjAAGICBEIaSIAAcgTgZAOAhCAAiIQ6p8ABGAIEQj1TQACcEIiEOqXAATgpEQg1CcBCMApiUCoPwIQgNMSgVBfBCAAwyICoX4IQACGTQRCfRCAAIyICITaJwABGDERCLVNAAIwKiIQapcABGDURCDUJgFYZV5++eW49tpro62tLRoaGuKZZ54p2H/jjTdGQ0NDwW3u3LkFYwYGBmLZsmUxefLkmDhxYixcuDD27t1bMKavry+6u7sjk8lEJpOJ7u7uOHjwYInPDqhHIhBqjwCsMp988knMmjUr1q5de9Ix8+fPj3379uVvzz77bMH+np6eWL9+faxbty42b94chw8fjq6urjh69Gh+zOLFi2PHjh2xYcOG2LBhQ+zYsSO6u7tLdl5AfROBUFvGVvoAKLRgwYJYsGDBKcc0NjZGNps94b5cLhePPPJIPPHEE3HVVVdFRMRPf/rTaG9vj+effz6uueaaePvtt2PDhg2xbdu2mDNnTkRE/PjHP46Ojo54991346KLLiruSQGpcCwCkySp9KEAp2EFsAa99NJL0dLSEhdeeGHcdNNNsX///vy+7du3x6effhqdnZ35bW1tbTFjxozYsmVLRERs3bo1MplMPv4iIubOnRuZTCY/BmA0rARCbbACWGMWLFgQf/M3fxNTp06N3bt3xz/+4z/GFVdcEdu3b4/Gxsbo7e2N8ePHxznnnFPwuNbW1ujt7Y2IiN7e3mhpaRny3C0tLfkxxxsYGIiBgYH8/f7+/iKeFVBPrARC9ROANWbRokX5P8+YMSMuueSSmDp1avziF7+I66+//qSPO/638hP9hn6q39xXr14d99133xkcOZAmIhCqm7eAa9yUKVNi6tSpsWvXroiIyGazMTg4GH19fQXj9u/fH62trfkxH3744ZDn+uijj/JjjnfPPfdELpfL3/bs2VPkMwHqjbeDoXoJwBp34MCB2LNnT0yZMiUiImbPnh3jxo2LjRs35sfs27cv3nzzzZg3b15ERHR0dEQul4tXX301P+aVV16JXC6XH3O8xsbGmDRpUsEN4HREIFQnbwFXmcOHD8dvfvOb/P3du3fHjh07orm5OZqbm2PlypXx13/91zFlypT47W9/G//wD/8QkydPjm984xsREZHJZGLp0qWxYsWKOPfcc6O5uTnuuuuumDlzZv5TwRdffHHMnz8/brrppnj44YcjIuLmm2+Orq4unwAGis7bwVB9BGCV+dWvfhWXX355/v7y5csjImLJkiXx0EMPxRtvvBGPP/54HDx4MKZMmRKXX355PPXUU9HU1JR/zAMPPBBjx46NG264IY4cORJXXnllPProozFmzJj8mCeffDLuvPPO/KeFFy5ceMrvHgQ4EyIQqktD4v+NjEJ/f39kMpnI5XLeDgaGTQRWDz/H0801gACUjWsCoToIQADKSgRC5QlAAMpOBEJlCUAAKkIEQuUIQAAqRgRCZQhAACpKBEL5CUAAKk4EQnkJQACqggiE8hGAAFQNEQjlIQABqCoiEEpPAAJQdUQglJYABKAqiUAoHQEIQNUSgVAaAhCAqiYCofgEIABVTwRCcQlAAGqCCITiEYAA1AwRCMUhAAGoKSIQzpwABKDmiEA4MwIQgJokAmH0BCAANUsEwugIQABqmgiEkROAANQ8EQgjIwABqAsiEIZPAAJQN0QgDI8ABKCuiEA4PQEIQN0RgXBqAhCAuiQC4eQEIAB1SwTCiQlAAOqaCIShBCAAdU8EQiEBCEAqiED4fwIQgNQQgfA5AQhAqohAEIAApJAIJO0EIACpJAJJMwEIQGqJQNJKAAKQaiKQNBKAAKSeCCRtBCAAhAgkXQQgAPwfEUhaCEAA+CMikDQQgABwHBFIvROAAHACIpB6JgAB4CREIPVKAALAKYhA6pEABIDTEIHUGwEIAMMgAqknAhAAhkkEUi8EIACMgAikHghAABghEUitE4AAMAoikFomAAFglEQgtUoAAsAZEIHUIgEIAGdIBFJrBCAAFIEIpJYIQAAoEhFIrRCAAFBEIpBaIAABoMhEINVOAAJACYhAqpkABIASEYFUKwEIACUkAqlGAhAASkwEUm0EIACUgQikmghAACgTEUi1EIAAUEYikGogAAGgzEQglSYAAaACRCCVJAABoEJEIJUiAAGggkQglSAAAaDCRCDlJgABoAqIQMpJAAJAlRCBlIsABIAqIgIpBwEIAFVGBFJqAhAAqpAIpJQEIABUKRFIqQjAKrJ69er4yle+Ek1NTdHS0hLXXXddvPvuuwVjkiSJlStXRltbW0yYMCEuu+yy2LlzZ8GYgYGBWLZsWUyePDkmTpwYCxcujL179xaM6evri+7u7shkMpHJZKK7uzsOHjxY6lMEYIREIKUgAKvIpk2b4vbbb49t27bFxo0b47PPPovOzs745JNP8mPuv//+WLNmTaxduzZee+21yGazcfXVV8ehQ4fyY3p6emL9+vWxbt262Lx5cxw+fDi6urri6NGj+TGLFy+OHTt2xIYNG2LDhg2xY8eO6O7uLuv5AjA8IpBia0iSJKn0QXBiH330UbS0tMSmTZvia1/7WiRJEm1tbdHT0xPf+c53IuLz1b7W1tb4wQ9+ELfcckvkcrk477zz4oknnohFixZFRMTvfve7aG9vj2effTauueaaePvtt+NLX/pSbNu2LebMmRMREdu2bYuOjo5455134qKLLjrtsfX390cmk4lcLheTJk0q3SQAkNfQ0BDF+mvbz/F0swJYxXK5XERENDc3R0TE7t27o7e3Nzo7O/NjGhsb49JLL40tW7ZERMT27dvj008/LRjT1tYWM2bMyI/ZunVrZDKZfPxFRMydOzcymUx+DADVx0ogxTK20gfAiSVJEsuXL4+vfvWrMWPGjIiI6O3tjYiI1tbWgrGtra3x/vvv58eMHz8+zjnnnCFjjj2+t7c3Wlpahvw3W1pa8mOONzAwEAMDA/n7/f39ozwzAM7EsQj0Bh5nwgpglbrjjjvi17/+dfz7v//7kH3H//Y3nN8Ijx9zovGnep7Vq1fnPzCSyWSivb19OKcBQAlYCeRMCcAqtGzZsvj5z38eL774Ypx//vn57dlsNiJiyCrd/v3786uC2Ww2BgcHo6+v75RjPvzwwyH/3Y8++mjI6uIx99xzT+Ryufxtz549oz9BAM6YCORMCMAqkiRJ3HHHHfH000/HCy+8ENOmTSvYP23atMhms7Fx48b8tsHBwdi0aVPMmzcvIiJmz54d48aNKxizb9++ePPNN/NjOjo6IpfLxauvvpof88orr0Qul8uPOV5jY2NMmjSp4AZAZYlARss1gFXk9ttvj5/97GfxH//xH9HU1JRf6ctkMjFhwoRoaGiInp6eWLVqVUyfPj2mT58eq1atirPPPjsWL16cH7t06dJYsWJFnHvuudHc3Bx33XVXzJw5M6666qqIiLj44otj/vz5cdNNN8XDDz8cERE333xzdHV1DesTwABUD9cEMhoCsIo89NBDERFx2WWXFWz/yU9+EjfeeGNERNx9991x5MiRuO2226Kvry/mzJkTzz33XDQ1NeXHP/DAAzF27Ni44YYb4siRI3HllVfGo48+GmPGjMmPefLJJ+POO+/Mf1p44cKFsXbt2tKeIAAlIQIZKd8DyKj4/iiA6jOSCPRzPN1cAwgAdcI1gQyXAASAOiICGQ4BCAB1RgRyOgIQAOqQCORUBCAA1CkRyMkIQACoYyKQExGAAFDnRCDHE4AAkAIikD8mAAEgJUQgxwhAAEgREUiEfwuYUTr2Tw319/dX+EgAGKlcLheZTCYiwr8fnFICkFE5cOBARES0t7dX+EgAOBOHDh3KxyDpIQAZlebm5oiI+OCDD/zgKIH+/v5ob2+PPXv2+EfaS8D8lpb5La1izW+SJHHo0KFoa2sr4tFRKwQgo3LWWZ9fPprJZPyAL6FJkyaZ3xIyv6VlfkurGPPrF/j08iEQAICUEYAAACkjABmVxsbGuPfee6OxsbHSh1KXzG9pmd/SMr+lZX4phobE578BAFLFCiAAQMoIQACAlBGAAAApIwABAFJGAAIApIwABABIGQEIAJAyAhAAIGUEIABAyghAAICUEYAAACkjAAEAUkYAAgCkjAAEAEgZAQgAkDICEAAgZQQgAEDKCEAAgJQRgAAAKSMAAQBSRgACAKSMAAQASBkBCACQMgIQACBlBCAAQMoIQACAlBGAAAApIwABAFJGAAIApIwABABIGQEIAJAyAhAAIGUEIABAyghAAICUEYAAACkjAAEAUkYAAgCkjAAEAEgZAQgAkDICEAAgZQQgAEDKCEAAgJQRgAAAKSMAAQBSRgACAKSMAAQASJn/BXJL0WQ6w21rAAAAAElFTkSuQmCC' width=640.0/>\n",
       "            </div>\n",
       "        "
      ],
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "solver = polar_fdm_2d(r_max = 10, n_r = 100, m_max = 10)\n",
    "\n",
    "# Set up potential V(x, y) = (x-x0)**2/2 + (y-y0)**2/2.\n",
    "# This gives an harmonic oscillator with spectrum 1, 2, 2, 3, 3, 3, 4, 4, 4, 4, ...\n",
    "#\n",
    "# V(x, y) = (x-x0)**2/2 + (y-y0)**2/2 = r**2/2 - x*x0 - y*y0 + x0**2/2 + y0**2/2\n",
    "# = r**2/2 - r*cos(theta)*x0 - r*sin(theta)*y0 + x0**2/2 + y0**2/2\n",
    "# = r**2/2 - x0*r*(exp(i*theta) + exp(-i*theta))/2 - x0*r*(exp(i*tjheta)-exp(-i*theta))/2i + x0**2/2 + y0**2/2 \n",
    "# = r**2/2 - r*(x0 - i*y0)/2 * exp(i*theta) - r*(x0 + i*y0)/2 * exp(-i*theta) + x0**2/2 + y0**2/2\n",
    "\n",
    "# For all alpha, this is just a shifted harmonic oscillator\n",
    "# and should produce evenly spaced energy levels with correct degeneracies 1, 2, 3. ...\n",
    "x0 = 0.2\n",
    "y0 = 0.1\n",
    "alpha = x0 - 1j*y0\n",
    "V_m = []\n",
    "# m = -1\n",
    "V_m.append(-alpha*0.5*solver.r_inner)\n",
    "V_m.append(0.5*solver.r_inner**2 + 0.5*np.abs(alpha)**2)\n",
    "V_m.append(-alpha.conjugate()*0.5*solver.r_inner)\n",
    "solver.set_potential(V_m)\n",
    "\n",
    "# Compute dense matrix version of Hamiltonian.\n",
    "H_mat, R_mat = solver.get_dense_matrix()\n",
    "\n",
    "# Plot sparsity pattern. Zoom in for details.\n",
    "plt.figure()\n",
    "plt.spy(H_mat)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ic| E[:11]: array([1.00109977, 1.9988004 , 1.99908408, 2.99868241, 2.99871346,\n",
      "                   3.00051106, 3.994058  , 3.9950714 , 3.99871128, 3.99871137,\n",
      "                   4.99324552])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([1.00109977, 1.9988004 , 1.99908408, 2.99868241, 2.99871346,\n",
       "       3.00051106, 3.994058  , 3.9950714 , 3.99871128, 3.99871137,\n",
       "       4.99324552])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compute the first eigenvalues and the differences between them.\n",
    "# For a HO this should give approximately, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1.\n",
    "E, U = eigh(H_mat)\n",
    "i = np.argsort(E)\n",
    "E = E[i]\n",
    "U = U[:,i]\n",
    "ic(E[:11])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[48], line 16\u001b[0m\n\u001b[1;32m     14\u001b[0m start \u001b[38;5;241m=\u001b[39m time()\n\u001b[1;32m     15\u001b[0m \u001b[38;5;66;03m#result = eigs(H, 10, R, which='LM', sigma=sigma, tol=1e-6, maxiter=1000)\u001b[39;00m\n\u001b[0;32m---> 16\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[43meigs\u001b[49m\u001b[43m(\u001b[49m\u001b[43mH\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m10\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mwhich\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mLM\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msigma\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msigma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtol\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1e-6\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmaxiter\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1000\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m     17\u001b[0m ic(time() \u001b[38;5;241m-\u001b[39m start)\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/scipy/sparse/linalg/_eigen/arpack/arpack.py:1352\u001b[0m, in \u001b[0;36meigs\u001b[0;34m(A, k, M, sigma, which, v0, ncv, maxiter, tol, return_eigenvectors, Minv, OPinv, OPpart)\u001b[0m\n\u001b[1;32m   1350\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m _ARPACK_LOCK:\n\u001b[1;32m   1351\u001b[0m     \u001b[38;5;28;01mwhile\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m params\u001b[38;5;241m.\u001b[39mconverged:\n\u001b[0;32m-> 1352\u001b[0m         \u001b[43mparams\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43miterate\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1354\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m params\u001b[38;5;241m.\u001b[39mextract(return_eigenvectors)\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/scipy/sparse/linalg/_eigen/arpack/arpack.py:746\u001b[0m, in \u001b[0;36m_UnsymmetricArpackParams.iterate\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    744\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m    745\u001b[0m         Bxslice \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mslice\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mipntr[\u001b[38;5;241m2\u001b[39m] \u001b[38;5;241m-\u001b[39m \u001b[38;5;241m1\u001b[39m, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mipntr[\u001b[38;5;241m2\u001b[39m] \u001b[38;5;241m-\u001b[39m \u001b[38;5;241m1\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn)\n\u001b[0;32m--> 746\u001b[0m         \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mworkd[yslice] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mOPa\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mworkd\u001b[49m\u001b[43m[\u001b[49m\u001b[43mBxslice\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    747\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mido \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m2\u001b[39m:\n\u001b[1;32m    748\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mworkd[yslice] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mB(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mworkd[xslice])\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py:234\u001b[0m, in \u001b[0;36mLinearOperator.matvec\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m    231\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m x\u001b[38;5;241m.\u001b[39mshape \u001b[38;5;241m!=\u001b[39m (N,) \u001b[38;5;129;01mand\u001b[39;00m x\u001b[38;5;241m.\u001b[39mshape \u001b[38;5;241m!=\u001b[39m (N,\u001b[38;5;241m1\u001b[39m):\n\u001b[1;32m    232\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mdimension mismatch\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m--> 234\u001b[0m y \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_matvec\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    236\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(x, np\u001b[38;5;241m.\u001b[39mmatrix):\n\u001b[1;32m    237\u001b[0m     y \u001b[38;5;241m=\u001b[39m asmatrix(y)\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/scipy/sparse/linalg/_eigen/arpack/arpack.py:1027\u001b[0m, in \u001b[0;36mIterOpInv._matvec\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m   1026\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_matvec\u001b[39m(\u001b[38;5;28mself\u001b[39m, x):\n\u001b[0;32m-> 1027\u001b[0m     b, info \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mifunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mOP\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtol\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtol\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1028\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m info \u001b[38;5;241m!=\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[1;32m   1029\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mError in inverting [A-sigma*M]: function \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m   1030\u001b[0m                          \u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m did not converge (info = \u001b[39m\u001b[38;5;132;01m%i\u001b[39;00m\u001b[38;5;124m).\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m   1031\u001b[0m                          \u001b[38;5;241m%\u001b[39m (\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mifunc\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m, info))\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/scipy/sparse/linalg/_eigen/arpack/arpack.py:953\u001b[0m, in \u001b[0;36mgmres_loose\u001b[0;34m(A, b, tol)\u001b[0m\n\u001b[1;32m    951\u001b[0m b \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39masarray(b)\n\u001b[1;32m    952\u001b[0m min_tol \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1000\u001b[39m \u001b[38;5;241m*\u001b[39m np\u001b[38;5;241m.\u001b[39msqrt(b\u001b[38;5;241m.\u001b[39msize) \u001b[38;5;241m*\u001b[39m np\u001b[38;5;241m.\u001b[39mfinfo(b\u001b[38;5;241m.\u001b[39mdtype)\u001b[38;5;241m.\u001b[39meps\n\u001b[0;32m--> 953\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mgmres\u001b[49m\u001b[43m(\u001b[49m\u001b[43mA\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mb\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtol\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mmax\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mtol\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmin_tol\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43matol\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m<decorator-gen-5>:2\u001b[0m, in \u001b[0;36mgmres\u001b[0;34m(A, b, x0, tol, restart, maxiter, M, callback, restrt, atol, callback_type)\u001b[0m\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/scipy/_lib/_threadsafety.py:44\u001b[0m, in \u001b[0;36mReentrancyLock.decorate.<locals>.caller\u001b[0;34m(func, *a, **kw)\u001b[0m\n\u001b[1;32m     42\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mcaller\u001b[39m(func, \u001b[38;5;241m*\u001b[39ma, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkw):\n\u001b[1;32m     43\u001b[0m     \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28mself\u001b[39m:\n\u001b[0;32m---> 44\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkw\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/scipy/sparse/linalg/_isolve/iterative.py:694\u001b[0m, in \u001b[0;36mgmres\u001b[0;34m(A, b, x0, tol, restart, maxiter, M, callback, restrt, atol, callback_type)\u001b[0m\n\u001b[1;32m    692\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m (ijob \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m3\u001b[39m):\n\u001b[1;32m    693\u001b[0m     work[slice2] \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m=\u001b[39m sclr2\n\u001b[0;32m--> 694\u001b[0m     work[slice2] \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m sclr1\u001b[38;5;241m*\u001b[39m\u001b[43mmatvec\u001b[49m\u001b[43m(\u001b[49m\u001b[43mwork\u001b[49m\u001b[43m[\u001b[49m\u001b[43mslice1\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    695\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m resid_ready:\n\u001b[1;32m    696\u001b[0m         \u001b[38;5;28;01mif\u001b[39;00m callback_type \u001b[38;5;129;01min\u001b[39;00m (\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpr_norm\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlegacy\u001b[39m\u001b[38;5;124m'\u001b[39m):\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py:234\u001b[0m, in \u001b[0;36mLinearOperator.matvec\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m    231\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m x\u001b[38;5;241m.\u001b[39mshape \u001b[38;5;241m!=\u001b[39m (N,) \u001b[38;5;129;01mand\u001b[39;00m x\u001b[38;5;241m.\u001b[39mshape \u001b[38;5;241m!=\u001b[39m (N,\u001b[38;5;241m1\u001b[39m):\n\u001b[1;32m    232\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mdimension mismatch\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m--> 234\u001b[0m y \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_matvec\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    236\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(x, np\u001b[38;5;241m.\u001b[39mmatrix):\n\u001b[1;32m    237\u001b[0m     y \u001b[38;5;241m=\u001b[39m asmatrix(y)\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py:591\u001b[0m, in \u001b[0;36m_CustomLinearOperator._matvec\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m    590\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_matvec\u001b[39m(\u001b[38;5;28mself\u001b[39m, x):\n\u001b[0;32m--> 591\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m__matvec_impl\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/scipy/sparse/linalg/_eigen/arpack/arpack.py:1004\u001b[0m, in \u001b[0;36mIterOpInv.__init__.<locals>.mult_func_M_None\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m   1003\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mmult_func_M_None\u001b[39m(x):\n\u001b[0;32m-> 1004\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mA\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmatvec\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;241m-\u001b[39m sigma \u001b[38;5;241m*\u001b[39m x\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py:234\u001b[0m, in \u001b[0;36mLinearOperator.matvec\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m    231\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m x\u001b[38;5;241m.\u001b[39mshape \u001b[38;5;241m!=\u001b[39m (N,) \u001b[38;5;129;01mand\u001b[39;00m x\u001b[38;5;241m.\u001b[39mshape \u001b[38;5;241m!=\u001b[39m (N,\u001b[38;5;241m1\u001b[39m):\n\u001b[1;32m    232\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mdimension mismatch\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m--> 234\u001b[0m y \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_matvec\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    236\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(x, np\u001b[38;5;241m.\u001b[39mmatrix):\n\u001b[1;32m    237\u001b[0m     y \u001b[38;5;241m=\u001b[39m asmatrix(y)\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py:591\u001b[0m, in \u001b[0;36m_CustomLinearOperator._matvec\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m    590\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_matvec\u001b[39m(\u001b[38;5;28mself\u001b[39m, x):\n\u001b[0;32m--> 591\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m__matvec_impl\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\n",
      "Cell \u001b[0;32mIn[48], line 5\u001b[0m, in \u001b[0;36mapply_hamiltonian\u001b[0;34m(psi_vec, sigma)\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mapply_hamiltonian\u001b[39m(psi_vec, sigma \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0\u001b[39m):\n\u001b[0;32m----> 5\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43msolver\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mapply_hamiltonian\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpsi_vec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mreshape\u001b[49m\u001b[43m(\u001b[49m\u001b[43msolver\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mshape\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mgen_evp\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m\u001b[38;5;241m.\u001b[39mflatten() \u001b[38;5;241m-\u001b[39m sigma\u001b[38;5;241m*\u001b[39mpsi_vec\n",
      "Cell \u001b[0;32mIn[45], line 88\u001b[0m, in \u001b[0;36mpolar_fdm_2d.apply_hamiltonian\u001b[0;34m(self, psi, gen_evp)\u001b[0m\n\u001b[1;32m     86\u001b[0m m_index \u001b[38;5;241m=\u001b[39m m \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mm_max\n\u001b[1;32m     87\u001b[0m \u001b[38;5;66;03m#ic(m, m_index)\u001b[39;00m\n\u001b[0;32m---> 88\u001b[0m result[:,m_index] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mT_m\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;28;43mabs\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mm\u001b[49m\u001b[43m)\u001b[49m\u001b[43m]\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m@\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mpsi\u001b[49m\u001b[43m[\u001b[49m\u001b[43m:\u001b[49m\u001b[43m,\u001b[49m\u001b[43mm_index\u001b[49m\u001b[43m]\u001b[49m \n\u001b[1;32m     89\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m m2 \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(\u001b[38;5;241m-\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mV_m_max, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mV_m_max\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m):\n\u001b[1;32m     90\u001b[0m     m2_index \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mV_m_max \u001b[38;5;241m+\u001b[39m m2\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.10/site-packages/numpy/matrixlib/defmatrix.py:168\u001b[0m, in \u001b[0;36mmatrix.__array_finalize__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    163\u001b[0m     ret \u001b[38;5;241m=\u001b[39m N\u001b[38;5;241m.\u001b[39mndarray\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__new__\u001b[39m(subtype, shape, arr\u001b[38;5;241m.\u001b[39mdtype,\n\u001b[1;32m    164\u001b[0m                             buffer\u001b[38;5;241m=\u001b[39marr,\n\u001b[1;32m    165\u001b[0m                             order\u001b[38;5;241m=\u001b[39morder)\n\u001b[1;32m    166\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m ret\n\u001b[0;32m--> 168\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__array_finalize__\u001b[39m(\u001b[38;5;28mself\u001b[39m, obj):\n\u001b[1;32m    169\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_getitem \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mFalse\u001b[39;00m\n\u001b[1;32m    170\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m (\u001b[38;5;28misinstance\u001b[39m(obj, matrix) \u001b[38;5;129;01mand\u001b[39;00m obj\u001b[38;5;241m.\u001b[39m_getitem): \u001b[38;5;28;01mreturn\u001b[39;00m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "from scipy.sparse.linalg import eigs, cg\n",
    "\n",
    "sigma = 1.0\n",
    "def apply_hamiltonian(psi_vec, sigma = 0):\n",
    "    return solver.apply_hamiltonian(psi_vec.reshape(solver.shape), gen_evp=False).flatten() - sigma*psi_vec\n",
    "#def apply_weight(psi_vec):\n",
    "#    return solver.apply_weight(psi_vec.reshape(solver.shape)).flatten()\n",
    "\n",
    "H = LinearOperator((solver.n_dof, solver.n_dof), matvec=apply_hamiltonian, dtype=np.complex128)\n",
    "#R = LinearOperator((solver.n_dof, solver.n_dof), matvec=apply_weight, dtype=np.complex128)\n",
    "\n",
    "\n",
    "from time import time\n",
    "start = time()\n",
    "#result = eigs(H, 10, R, which='LM', sigma=sigma, tol=1e-6, maxiter=1000)\n",
    "result = eigs(H, 10, which='LM', sigma=sigma, tol=1e-6, maxiter=1000)\n",
    "ic(time() - start)\n",
    "\n",
    "#psi_vec = np.random.rand(solver.n_dof)\n",
    "\n",
    "# for k in range(10):\n",
    "#     psi_vec = cg(H, psi_vec, tol=1e-6, maxiter=1000)[0]\n",
    "#     psi_vec /= np.linalg.norm(psi_vec)\n",
    "#     E = psi_vec.conjugate().T @ (H @ psi_vec)\n",
    "#     ic(E)\n",
    "    \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'result' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[8], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m np\u001b[38;5;241m.\u001b[39mdiff(\u001b[43mresult\u001b[49m[\u001b[38;5;241m0\u001b[39m])\n",
      "\u001b[0;31mNameError\u001b[0m: name 'result' is not defined"
     ]
    }
   ],
   "source": [
    "np.diff(result[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'result' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[9], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m psi \u001b[38;5;241m=\u001b[39m \u001b[43mresult\u001b[49m[\u001b[38;5;241m1\u001b[39m][:,\u001b[38;5;241m0\u001b[39m]\u001b[38;5;241m.\u001b[39mreshape(solver\u001b[38;5;241m.\u001b[39mshape)\n\u001b[1;32m      2\u001b[0m plt\u001b[38;5;241m.\u001b[39mfigure()\n\u001b[1;32m      3\u001b[0m plt\u001b[38;5;241m.\u001b[39mplot(solver\u001b[38;5;241m.\u001b[39mr0, solver\u001b[38;5;241m.\u001b[39mG_neumann \u001b[38;5;241m@\u001b[39m psi[:,\u001b[38;5;241m0\u001b[39m]\u001b[38;5;241m.\u001b[39mreal)\n",
      "\u001b[0;31mNameError\u001b[0m: name 'result' is not defined"
     ]
    }
   ],
   "source": [
    "psi = result[1][:,0].reshape(solver.shape)\n",
    "plt.figure()\n",
    "plt.plot(solver.r0, solver.G_neumann @ psi[:,0].real)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ic| H_solver.shape: (100, 100)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(100, 100)"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_solver = solver.T_m[0] + solver.V_m[0]\n",
    "ic(H_solver.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}