Fix missing import in finiterectlat-scatter.py

[qpms.git] / examples / size_vs_lMax / AuSphere / _mode_r_dep.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import re\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.constants import hbar, e as eV, c\n",
    "eh = eV/hbar\n",
    "import glob\n",
    "def ri(z): return (z.real, z.imag)\n",
    "#m = re.compile(r\"([^_]+)_r([0-9.]+)nm_\")\n",
    "#removek = re.compile(r\"(k\\([^)]+\\)um-1_)\")\n",
    "remover = re.compile(r\"r[0-9.]+nm_\")\n",
    "\n",
    "\n",
    "markerdict = {\n",
    "    4: \"3\",\n",
    "    -4: \"4\",\n",
    "    3: \"^\",\n",
    "    -3: \"v\",\n",
    "    -2: 'x',\n",
    "    2: '+',\n",
    "    1: 's',\n",
    "    -1: 'd',\n",
    "}\n",
    "\n",
    "\n",
    "D4h_irlabels = {\"B2''\":\"$B_2''$\",\n",
    "                \"B2'\":\"$B_2'$\",\n",
    "                \"A1''\":\"$A_1''$\",\n",
    "                \"A1'\":\"$A_1'$\",\n",
    "                \"A2''\":\"$A_2''$\",\n",
    "                \"B1''\":\"$B_1''$\",\n",
    "                \"A2'\":\"$A_2'$\", \n",
    "                \"B1'\":\"$B_1'$\",\n",
    "                \"E'\":\"$E'$\",\n",
    "                \"E''\":\"$E''$\",}\n",
    "\n",
    "\n",
    "prop_cycle = plt.rcParams['axes.prop_cycle']\n",
    "colors = prop_cycle.by_key()['color']\n",
    "colordict = {i: colors[(i+1)] for i in range(-4,8)}\n",
    "\n",
    "def markerfun(b):\n",
    "    if b in markerdict.keys():\n",
    "        return markerdict[b]\n",
    "    else: return 'X'\n",
    "\n",
    "def colorfun(b):\n",
    "    if (b+1) in colordict.keys():\n",
    "        return colordict[b+1]\n",
    "    else: return colordict[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sph_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L1_cn410.npz 203\n",
      "sph_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L2_cn410.npz 251\n",
      "sph_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L3_cn410.npz 251\n",
      "sph_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L4_cn410.npz 251\n",
      "sph_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L5_cn410.npz 251\n"
     ]
    }
   ],
   "source": [
    "allfiles=glob.glob('*sph*k(0_0)*.npz')\n",
    "allgraphs=dict()\n",
    "for f in allfiles:\n",
    "    base = remover.sub('', f)\n",
    "    if base in allgraphs.keys():\n",
    "        allgraphs[base] += 1\n",
    "    else:\n",
    "        allgraphs[base] = 1\n",
    "for k in sorted(allgraphs.keys()):\n",
    "    print(k, allgraphs[k])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['sph_r281nm_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L4_cn410.npz', 'sph_r121nm_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L3_cn410.npz', 'sph_r67nm_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L3_cn410.npz', 'sph_r208nm_p580nmx580nm_mAu_n1.52_b-2_k(0_0)um-1_L3_cn410.npz'] collected_XXXsph_rXXX_p580nmx580nm_mAu_n1.52_bX2_k(0_0)um-1_LX_cnXXX.pdf\n"
     ]
    }
   ],
   "source": [
    "projectors = dict()\n",
    "projectors_list = dict()\n",
    "lMaxes = [lMax for lMax in range(1,6)]\n",
    "for lMax in lMaxes:\n",
    "    proj = np.load('projectors_D4h_lMax%d.npz' % lMax)\n",
    "    irlabels = sorted(proj.keys())\n",
    "    proj =  {f: proj[f] for f in irlabels}\n",
    "    proj_list = [proj[irlabels[i]] for i in range(len(proj))]\n",
    "    projectors[lMax] = proj\n",
    "    projectors_list[lMax] = proj_list\n",
    "globpattern = '*sph_r*_p580nmx580nm_mAu_n1.52_b?2_k(0_0)um-1_L?_cn???.npz'\n",
    "filenames=glob.glob(globpattern)\n",
    "plotfilename = 'collected_' + globpattern.replace('*', 'XXX').replace('?', 'X').replace('npz','pdf')\n",
    "print(filenames[:4], plotfilename)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#projectors\n",
    "#glob.glob('cyl_r100nm*L3*3100.npz')\n",
    "#glob.glob('sph_r100*m5*.npz')\n",
    "#dat['meta'][()],list(dat.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inpure result detected [1.         1.         0.92521572 1.         0.99999627 0.99990293\n",
      " 0.99946049] [1.59712906e-05 3.60193407e-05 2.48341492e-01 1.21848930e-03\n",
      " 3.81805601e-03 2.42649228e-02 2.99534246e-02]\n",
      "inpure result detected [1.         1.         0.99999998 0.99999999 0.99999996 0.96608887\n",
      " 0.99999852 0.99999397 0.99998951 0.99999912 0.99982435] [2.66223026e-04 2.12357147e-05 3.54211968e-05 1.06651057e-04\n",
      " 2.79595790e-04 2.41939163e-01 2.17645058e-03 3.41541473e-03\n",
      " 1.14507609e-02 1.49639498e-02 2.33483138e-02]\n",
      "inpure result detected [1.         0.99999999 1.         0.97991334 0.99999996 0.9999989\n",
      " 0.99999983 0.99999966 0.99999322 0.99999721 0.99999653] [3.28735741e-04 2.66532534e-05 2.47011478e-05 1.45012420e-01\n",
      " 2.44785416e-04 7.05405359e-04 1.60203586e-03 1.71245137e-03\n",
      " 1.03244480e-02 9.18732728e-03 1.18651583e-02]\n",
      "inpure result detected [1.         1.         0.99999998 0.99999961 0.93267685 0.99999964\n",
      " 0.99999822 0.99921774 0.99995547 0.99997301] [5.22490396e-04 3.01556792e-05 4.88795563e-05 6.29703960e-04\n",
      " 2.34414238e-01 3.72766210e-03 4.72444059e-03 7.62106094e-02\n",
      " 6.32796684e-02 5.63231562e-02]\n"
     ]
    }
   ],
   "source": [
    "plotdata = {}\n",
    "for file in filenames:\n",
    "    dat = np.load(file, allow_pickle=True)\n",
    "    kx = dat['meta'][()]['k'][0]\n",
    "    radius = dat['meta'][()]['radius']\n",
    "    b = dat['meta'][()]['band_index']\n",
    "    eigvals = dat['eigval']\n",
    "    lMax = dat['meta'][()]['lMax']\n",
    "    residuals = dat['residuals']\n",
    "    ef =dat['empty_freqs']\n",
    "    eigvecs = dat['eigvec']\n",
    "    irweights = []\n",
    "    #for proj in projectors_list[lMax]:\n",
    "    #    try:\n",
    "    #        irweights.append(np.linalg.norm(np.tensordot(proj, eigvecs, axes=(-1, -1)), axis=0,ord=2) if len(proj) != 0 else np.zeros((len(eigvecs),)))\n",
    "    #    except ValueError as err:\n",
    "    #        print(proj, len(proj))\n",
    "    #        raise err\n",
    "    irweights = np.array(irweights)\n",
    "    #print(irweights)\n",
    "    irweights = np.array([np.linalg.norm(np.tensordot(proj, eigvecs, axes=(-1, -1)), axis=0,ord=2) if len(proj) != 0 else np.zeros((len(eigvecs),)) for proj in projectors_list[lMax]]).T\n",
    "    irclass = np.argmax(irweights, axis=-1)\n",
    "    purities = np.amax(irweights, axis=-1)\n",
    "    if (np.any(purities < 0.98)):\n",
    "        print(\"inpure result detected\", purities, residuals)\n",
    "    #print(purities)\n",
    "    \n",
    "    #for i in range(len(residuals)): \n",
    "    #    if residuals[i] < 0.01:\n",
    "    #        vec = eigvecs[i]\n",
    "    #        for irlabel, proj in projectors.items():\n",
    "    #            print(irlabel, np.linalg.norm(np.dot(proj, vec))) #maybe some conj() here?\n",
    "    #        print('--->', irlabels[irclass[i]])\n",
    "\n",
    "    \n",
    "    plotdata[(lMax,radius)] = (eigvals, residuals, b, ef, irclass,)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "collected_XXXsph_rXXX_p580nmx580nm_mAu_n1.52_bX2_k(0_0)um-1_LX_cnXXX.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbsAAADOCAYAAACn8azYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd8VFXagJ/3ZgIkkAChgyAgoYckgAio6H4UsWBBIURRcXctWNcuyKorCioWFATdFaU3UVewo66iAiItgVBCaEoHgSSQhDCZ8/1x7qSQ6ZlJArkPP34zSebee86de8rbRSmFhYWFhYXFuYxR0Q2wsLCwsLAINdZiZ2FhYWFxzmMtdhYWFhYW5zzWYmdhYWFhcc5jq+gGWFhYWFhUDtasWdPQZrO9B3Tm7BSGHMBGu93+927duh0q/gdrsbOwsLCwAMBms73XuHHjDg0aNDhmGMZZ56rvcDjk8OHDHQ8cOPAecG3xv52NK7eFhYWFRWjo3KBBg6yzcaEDMAxDNWjQIBMtmZb8WwW0x8LCwsKicmKcrQudE7P9pdY2a7GzsLCwsCgTkZGRiRXdBm9Yi52FhYWFRdCx2+0lfnY4HBQUFFRQa6rYYicif4hIQkW3oyKp6vfA6n/V7j9Y9yCUfPbZZ1EXXXRR20GDBrVq165dp61bt1Zr3bp1p+HDh7fo1KlTx+3bt1f7+OOPoxMSEtp37Nixw5VXXtk6MzPTAGjWrFncyJEjm8XFxXWIi4vrsHHjxurBbFuVWexEpC7QGNhc0W0BEJH7RWS1iJwSkenldM1Kcw9EpLqITBOR3SKSLSLrROTKEF+z0vQfQERmi8h+EckSkXQR+XuIr1ep+u9ERGJFJE9EZpfDtSrVPRCRH8y+nzD/b63oNpWV1NTUmhMmTNi7ffv2NIBdu3bVuOOOO/7cvHnzpqioKMe4ceOaLFu2LH3Tpk2bu3btmjN27NhGzmOjo6MLNmzYsPnuu+8+9MADDzQPZruqUuhBHJCulDpV0Q0x2Qe8AFwBRJTTNSvTPbABfwCXAb8DVwELRSROKbUrRNesTP0HGA/8TSl1SkTaAz+IyDql1JoQXa+y9d/J28Bv5XStyngP7ldKvVfRjQgWXbp0Odm+fft8589NmjTJ79u370mAH374oeb27dtr9OjRoz3A6dOnpVu3biecn7399tuPAtx5551Hx4wZYy12AdIF2FDRjXCilPoYQES6A+eV02UrzT1QSp0Eniv2q89EZCfQDdgVostWmv4DKKXSiv9o/r8ACNViV6n6DyAiw4DjwHKgTTlcstLdg3ONyMhIh7uflVJccsklWUuWLNnp6ljDKFI2ikhQvUKrjBoT/ZCnhurkIvKZiBx38/+zUF3XTyrtPRCRRkBbIM3T58pIpeu/iEwRkRxgC7Af+CJU7aOS9V9EooHngUdD1SYXVKp7YDJeRI6IyC8icnmo2lYZuPzyy0+uXr26ltMel52dbaSmphba5mbOnBkDMG3atLqJiYkng3ntqiTZxQGLQ3VypdQ1oTp3EKmU90BEwoE5wAyl1JbgtqoEla7/Sql7ReQBoBdwORBK9Vpl6/9YYJpS6g8RCUWTXFHZ7sGTwCYgHxgGLBGRBKXU9qA3rhLQtGlT+7vvvrtr2LBhrfPz8wXg2Wef3dulS5dTAKdOnZIuXbq0dzgcMn/+/B3BvLZUheKtokdSFtBZKbW7ottTHBF5AThPKTUixNeplPdARAxgLhANXKeUOh2i61TK/hdHRN4BNiml3grBuStV/01vyDlAolIqX0SeA9oopYaH8JqV6h64QkS+Aj5XSk2qiOunpKTsio+PP1IR127WrFnc6tWrNzdp0sTu/dOeSUlJqR8fH9+y+O+qihqzFeBQSu0WkctFZLGIfCwi60XkehH5WkRWiUg9EWkkIv8TkWUislBEwkSkv4hMAxCRD0Sk/5kXEJEvi3lUnfn/y3LvcWkq3T0wJ59pQCPgxlAtdJW1/y6woW12oaCy9f9yoCXwu4gcAB4DbhSRtSHqf2W8B65QQLmJuVWJqrLYnWmUjlBKDQamALcrpa4A5gPXAceA/kqpPmiPyb8opZYCuSLyNpBr/lwCpdSVSqlabv6XcqkXEZuI1ADCgDARqSEioVQrV7p7AEwFOgCDlFK5Qe1taSpV/0WkoYgME5Fa5kR6BZAMfB+KzlPJ+g/8G72wJ5j/3wE+R3snh4pKdQ9EpI6IXOEc+yJyC9AH+DoUna/s7N27d0MwpDp3VBWbXRwljdLO9/vOeN8MiAGmikgM0BRYZ/79LbQTQesgtWkM8Gyxn4cD/6Kkh2IwqVT3QETOB+5G26gOSJHN5m6l1Jyynt8Flar/6B38SPQkbwC7gX8opT4NwrldUan6r5TKAXKcP4vICSBPKXW4rOf2QKW6B0A4OvyoPVBgnvd6pdRZH2tXGakSi51SauyZv3LzXoBbgG+UUlNFZBJFKoVXgXuAl4GkILTpOUK3sLm6XqW6B6bNpNzUNZWw/4fRMYblQmXrv4v2PRfM87m5RqW6B+YzcGFZzmHhO1VFjekP3wH3isinQAMAERkJ/KKU+jewW0SSK7KB5UBVvwdW/6t2/8G6B+ccVcIb08LCwsLCOxXpjRlMqrI3poWFhYVFFcZa7CwsLCwszhpmzpxZJ5DjKmyxE5H3ReSQiGz08rkLRaRARG4qr7ZZWFhYWFQcEyZMqN+gQYMu7du379iiRYvO/fr1uyAvL0/S0tKqL126NDqQc1akZDcdGOjpAyIShvZ6qpJxJxYWFhZVkdTU1MgxY8bs3bJly6adO3duTE9Pj1i1alXEkiVLoh955JFDgZyzwkIPlFLLRKSll489AHyEH+659evXVy1bejvt2ceaNWuOKKUa+PLZc/EeVPX+g+/3oLL0PwcdHxwZpGnmbHwGjnGKugSnBqk//T/b2bx5c8TIkSMPA6SlpVVXStG5c+e8VatWRSYmJuYFcs5KG2cnIs2AG4D/w4/FrmXLlqxevTpk7aooRMTnXH7n4j2o6v0H3+9BZej/ev7kGnSCkVn05Ut+J4U/aUQkE+mJLYCp52x7Bhawg5v5ge6cTydqs4GjADxBPN1piPgZZupP/4PBX++n+cbNRHr7nFKQk4sRGYHDWz7vzh3IeX8yf3g7Z0ZGRsSIESNa2e12OXjwYLVFixZti4mJcTz22GMBe4pW2sUOmAg8qZQqEC93UETuAu4CaNGiRTk0rfJRKe+BwwEfHoYhDcAIrca8Uva/HKls/Y8nhiX0A+Br9vBPUgr/9hV7iacujahBQ6qTRiYNiOBiGpFMG4wArSuV7R4MoSWf0op5lCzd9l/2cBkNqEM16lOdI+RzlDz60oynSCCMsApqcWDk5GJk7CCiTWtya0bi8H6EZzIyMsJjYmJOp6enbwKYPHlyvbFjxza58sort5XlvJV5sesOzDcXuvrAVSJiV0r998wPmkGe/wbo3r17lQwcDNk9cDjgw8Uw5Fr/F6wFB2H4ZnB0gOQmQWuSK8r1GVAKUlIgPh7KrzSNR8q1/w4HfPghDBni9pkQhETqA9CFutgpIIU/WcsxtnOS7ZQuVTaVbYzkFxoRTjYF9KYho+jqsxRU2eYBA4NZ9KEzdbArO9kpx/giPpNNks13lM6K9j+O8AIpRGNQh+rUIozBXFBhC6AvEhjox2HlaiJ6dic3GHva1atXR8bGxhaqKrt165YzadKkRmU9b6Vd7JRSrZzvRWQ68Jmrhc4ixHy4GIaNhPFvwaqvIDwc1q/Xf0tIcD/Z2+0wfX9w2qAUNSEiOCcL7PqkpEF8J93fdetg4ED46ivo2tXH409Q0LkWr+4SHmsFYf7OXZXpHixYAMOH61ku2XsSkTDCGEM3lILV2ad5JeoHHKIKJbufOcz3HEIB2SiyyQfgE/bzXz6nLgatiIKa1Suu/wFu+sIIYzSJkJKNunEDI6ZfxJQuezkQ/QciFEp2P7CfPykgD8jDwSF0XvS1pDCHDG4plyLugWEY0LsHQUvknpqaGtGuXbtcAIfDwbRp0+r16dMna9asWXU+//zz2ocPH7bdd999hwcPHpzlz3krbLETkXnoMh/1RWQPOilyOIBS6p2KapfFGQy5Fl58E1I2Q/crYPprcM01kJcH337rfrJ/eBt8cxwG1IGkMm7K1q4lVlcxrxhS0mDw36BnV5g5CbZsgSNH9Ksvi13KCRicyszR7XmqSQwOhFH+zF0FBfDII7QOXfkf76xNhf5JsHSB/tnhgOXLISnJ58k/JRuS1oRzd/P+PN666DCFYhX7uZ0fycJOTYQd5ONAJ6w8ioOjZEK7xh1C0jdfmPcJ3PoA2Avglhv9Pz6+Fukz4uifU4tTq2qz9MKOdC0WLXaa09zMD/zAHxgYnAaOmRrBzZxkDCnQtE6ZpZuzgY0bN0b8+uuvUUuXLq0jIqpr164np06duicqKspx6623Hj98+HDYfffdd95Zs9gppXzOK6dCXNjUwgOGARdfCBu26P+btsMDD8CYMfD115CYWFq6czjAoXTq3Nsals1eV1AATz9NWEVqIeI7Qbc4mPdfqBMNl5oS7YoVMGyY9/7F14IXWpF0XxrjxiWyrlEUBQU+SncOh5ai5s8nCzKD0p9A+PpHOJYJ94+GHz+B+++HKVOgd2+fpDuA+Ci4szmMSocWEZDcTP9eEC6iKVsoOs9pTnMLP7CNY2znBNnOD1YUK9do6XbUOBgyCKpV8+94EdpeHMX96TBmO7y/GxKiix6dcML5kKLyeAUU8CJr+YJ9NKUa3WjImH3HDwaxR5WWxYsX7/T099GjRzd58MEH/a6OYWVQqQgcDm3PcpTZlls+15k0Djq10+9fngJ9+0Lt2vDKK1qldyYLDsKUA3pbbpTB1uCc6L/+mmw/J/oU1qMIktlGBDKz9fslS+HGG+Hee2HyZK3S8+X4XXlEHLfzl03H+fAQvOSrqX3+fP2/Rg2OgP+eaAUF8PIu/VoWWpor08q18PAz0KuX/tmP3LoicH51/Vj88qfnxzKccBbSn3UM5TgjmEsf2P3njsA7UEbeeB6aN4U/9kH3gQGNKRFoVUv3/+29MHeP+8+GEcYzXMhKruNjruRpugXe9nMEh8PByJEjm1199dWZl1xySY73I0pyTi52CqUnO0dB+Swq/uBwaBXfLWnaU5Fi7Q3W5Iw52X94CG5O09cryz2w2WDUg/r9xi0w82N47DE4flyr8s7Eea17G5VNhblggf4fEcFe2OvrYTnkMJhrSGG9Xwu+UrB+g5v5+7NZ0KwR7NkPVw+Hnj3175cv9+3etqwBAqdjIwHF914me0B/YPFi/f7vf+ckvttFCp+pV3bBUzvglTJ6rQ+7Afr30e9/+hXatIF69WD3br+erfZREBUGU/bAgn2+HWNgkEwbOHLimD9NLhxTPj4DHr9/mw1efFK/37AZ5nzsT1MKL5B0MJuetfQFZu5VHFfBHffnMuPGjWu4bNmy6IULF9Z95ZVX/I43PCcXuxTWM5hr2L1glZ7sH0qvHAuewwEPbIG39sLdjbVLPkXtTWF9UC5TONkP2Qv3NNbXm+fdWcTjYE++Hvpeqt9v2QHnn6/fr1hR+t46f+4Z7ZMK0+Vi73DAL7/oxowZ49dE7ySdrah5+2BYmk/9X5cKAwbr11JUqwbXX6Xff/sTUL1IleeLdDesMTKyMbHLjoGC7zNhnoedPaAlugULYOhQmDjR+zWKUfhMxW/Xo/x/mdppyAMOByz42M1QMQy4fYh+v2ELTP8EHn5Yq7N96b9JYm14ynQ925njYrIvKICXJ5VZEi2x4Zl3AJLT9KsHUjbC4FvhlTfd3INbboLOpun4lcn+zynrszGuTeWOY1pLsNu+nl8cg8gsPu6V0vbR+Z9UjjmrEjFmzJhDaWlpm+fOnfv7E088UT5qTBGZLCK9Azm2PHDgIJccjtbJBgcweZ/e2VU08/Zr9R7oycNcCBSK0+QHdYeXTz5bjXSU08qxwrstd/0GuCZJv5bCMOCOofr9dz+Bo7q23U2dCgsXFn3O4YDZprZtVWn3cpfXZR2DGMB6iqlEFyzQ546K0p6PfhBJJOOYwKM8wOGVpuZr5iGvk8fWbfDnUf3qkonPQydzsvt5lX+qPMMAEZ786HcGp2sBZfpeD00qLtUlJPjtvukcA46BdaFfbVh6TEv4HvjwvzD8bpjwlrsNz2Doe4l+n54BLVoUOar4ODGLQKtIrcqbf3gdy9UVZBb/3l+dAk+Ngz7XeV2cvaFQLOVrHMuP6gsu9ywYxneGu26H0WNh4ScuPmAY8MQD+v3GrdppxV/sp/nb3N30jLSTnh3P13uuJZstRWM/JU07AiXfC/MCkB4t3BKoZLcNeE1EdonIyyKSEMxGlZUMtpFFFhlXHNGqNKj4XZLDATPNBTcuAt7Q7ngKxTbSsREeVPu7IDzGQ+y+6E/9i4tq+dTEnDwPtyp5MPQzJ7uZpvt5ZKRWaTmZtx++OQb9axf20RvpbOVP/iSdrfoXSsHOnbohTz2lnWAC4ChH+d8be7VH6DfHvG542sVC/Xr61SU2G1x+sX7/y2/+q/IuqkUYMLiJAwG+PQbz3ClnFyzQsWxJSfDEE97PfQbOMfCd8R3qtkbaWahnlMdjhlwPz4+CVye7kW4NA+5I1ivWdz+73/B4IakpjIuF8DA7p1QmWWpr0R8fuxcuSoTla+Chf/p8zjOJJJKRPMC/GMPXV5uZVMTw+D2JQL/LoWYktGnt5kO3FFvwf1nlX6MSouDyeoQtPMxDC3YBioMnm5DmeKhIuovvBP1MDcqMDyt+3jqHCGixU0q9qZTqBVwGHAU+EJHNIvKMiFSci7jJTQzlfWYTa7RH9a6jB/qKrIp9cObt1674/aJg7YV64kSrm57iMV7iNeIJzp4hkkgmMBFQKFuY/pZnH/a6U962HTIzYen/3OzsDQPamgvYtz/D9Plw8iTMnFl0b50SZNvIwj56I5a21KY2sc7ogrVr4aWXoE4dLdUFELjdlrbUpg5tbW21R6gADs+qscQu8NUifTm3wtrE56FTrKnK+9g/VZ55P5J/OUK/2vp+fXrQxbUKCuDTT/UfrrsugKA8PQae5QWmMomUOptMySbT4xgwDGjZAo4e8yDdJl8P/Uzb3axFrjc8LlAKjqzT/0VgYEPYnRXPpzvvZ+cJiiSbsDDt+QpQxmQcfelPNNE0GtgJ7msCU/Z7VWVm7ICcXPj2Bw9joL25G/pplV/Sp0MJ22Ma4wDaXFyTcAx++/NqcgvCwdl/Ebj2Cv1+6TJYYIUWB4sy2eyUUruVUi8rpRKBm9G5LDcHpWVlwMCgLe24iWtJiUkz3Z/2V5wqMz8fRpnetO2jSiwCCoWd08TSNqiSXVvaIRgcT6oG/erqhfbhDI/HDL0BXhwDk//jRpUJ8ObYIs9MqmmvxLff1pO9w1E0mfogSYLufwbbiCBSp4lyOODVVyE7G269NWCpTjCoQQTb2IZanlk02Xs6xrz9g4Z56H9x6S5jh1blgZZEvagz1dBG/JnUEnl7H/es1M/i4kOw9vgZH3zllSJbXVKSx3O6w8CgHwOwcxplmGEgUw76JN3WizHb626yv92strX0R5i+oPSGxwVHU+DLq+Czy+DwWu12/38x1ZiR/gqfHP2R46qY3ap3D/2avrNMqkxBtMbECNOdUcDyM292SYbeAC88DW/928MzMPF5iOugnbX+8YzP7dk+H356J4Y/qUW3lZnc0ECxMzueaZuexl58H5Z8A9w3Qj+QlmQXNMq02IlIuIgMEpE5wJdAOhBAxGXwKbSDXRGjd3WCjv2qCAalwR/50LyaG9Ve8NulB7qNbUYGang9/cseNT0eYxjQ/y96Z+t2jDk9M0Xg4h46zgq03WbeAZh6UN9rH6U6Ldk+ynhe1ZLt/PlaJTZ0KLzxRsDpuOJJ4CVe5SkeZXfvTN0mL2os0HNifj5sSfewdvUy3cCX/gS/H4LoaPdhGMU4kmLw9ZKm5Cob8TXtVANOAV8VL1hSUKAz1IjAtQGkaCuFwBX1YGRj/aOX/id2gTfHw4NPuVFlglZnD7hMvzeql9zwuCEmHloPAXs2rH5aNyO+lgDCodwGbMkudrOTb9Dn//YnHeYQIM7FbhvpqJ7m5suLKtc5Bk6d8jIGLjUXZOX7YpS9G+zKYHuzZjje3se1u44DBXy976+8vaNYhLlhmIuzguXnXkLziiJQB5X+IvI+sAedePUL4AKlVFJlSeklCOFU07u6XtHmelIBi53dDsrctr14vs+LQFmJJ4GXeZ1neZrU+lv1N73qpNfJLj0DMrP0q1uSr4d5U/RrUlLRZPeL6VBxbxOfQw60ZGvXkq1DaacMhyMgp4ziCGJKt8Lx5Op6w+ODdC+i/z80ysPOPnkwzJ0Cc9+GJx7UNrXsbB1k70G6y0yHuJwdRGCnVWo21zfUv//xaLGv5ZVXihb7YcP873jxviCEOyUbpzrfi/bAubfwqMo0DPh8Njz4Ny3l9OpVZGd1038R6PkGNBkA+76GFf+AljWhfvWdxNX/gZ0nzzj/8MH6fY/AVft6DLzGs4xmt22vz2Ng23Y47kmdD/DmC7r/b77ge3uegKZXCCf3hmEAw3JzGNxkOzVtR9l1ejMFxTfjTul2264yO+pYaALdNo4GVgAdlFKDlFJzlFK+ud6VE12I5yVeowvxpucjIc+875KHM2BpFgyoC7c0K/EnBw6+5RvCsAU9OYQg3MgQnmcccVf8Be5t6pNXartYiKmr37sd6IYBSdcXeZQ6J7u1uYCCXrV9vtdONa5ClXS1D8Ap40wKnwEjEXrW9km6T4iDieOLNtYuMQwtfSQP1gvyFVd4DrI3aT0UoobqHfyBdWFcU09PukuP6SBzZTelOsMIilRXYgygfN7wtW2ju9PWkxnOZtMTvc0G7dvrA1591WP/w8Kg7e36ffo06HsA+tU7Tv2IPXy0P6LkGrTKVGvOdhcL4Z2iMTCeFkkX6k2YD2PAqc6f5EmdX7z/PhIWBt1fhD1GXfYTRfbik9wQ047W0eu4+PwHWHi4mCo3+Xot3S79sUzSrUURgTqo/EUp9R/gmIgMF5FnAESkhYj0CGoLA2QDqTzDKD7iQ1RSQ5htptUrTx14fj78z/SGdJE2axEL+RdjuJcHg+acUpxUUniKR0k1NkDPWj47aUx6GR5/xsNAP5P27aH6JfDr80A6GL4t3ArFdywlk0wyCrbAjBn6D2WU6pwU9p8UCif75ce9euS1i/VTe5qYqIPsMzN1cmg3q6RhwPlzG7O/Z3MardxL39cO0t+0j72wU7HugaeDJtWBcwyM1mPA+ezNPORVUjAMiKiuJRyfEqQkJurNSWamV+n2gmHQKgkcuZD2Gvwlpi7/WvUZiw+2Zt6+Yse98TwMuFyrMhd86kMjXFNyDET5NAYMAwb8n6mUCbIyqEFXaN83l6+JZ2dKTdrXVGw5NoDxaxfy4+EGRdczDBhu2kZ7+JB/tQoxc+bMOt4/VZqyijpTgF5QmNQuG3i7jOcMCl2I51+M4ykeJcVI1Q/PbVsKs5aEHIcDuq+BtFM61MBFiZs2xBJNNH3pH5K0fyXi95ZnmZO953g7Eb2jzz/tx0BPTITe3YELocNWGNrQp8NSWM9UJvGceoGbHl4B33yjpaQgSHVOFGZwuXMBnnIg+I5KItCqlX596SWP0o2EGTS+N5JT2Mg74GDsBTr7eev09XR55xWIiNALZxC0EHoMvMiTPEJK0j6tXfjmmFdHpYQ4eO0FeOyfPm54RPT3FhXlVbo1DIh71EGdaifZuVAR9UlL7AW1aF5rA58cOFS0D7HZ4LYbTaenwAPMAxkDhceGwOohAq1ejKYVh8jNUNR/5k9ubXoQewHM/KMx64r7UBlnvFYhJkyYUL9BgwZd2rdv39H5f9WqVRFpaWnVly5dGh3IOct6Gy9SSt0H5AEopY4BfmZIDQ2CEEvbogf9xnpalXdjvfJpwLz9sMFM+vF4c5eTl9OJJLT5bc1z965t/qx8km79kmxEoP71wJOwJbeoBJAXnPa6fptbYHwwQzt6vPhiUKQ6KOmkkpK03+eYSxGoFq7tlj5PeElJMHKktt15kO7AlJyws+MjWPcs1AhT7GrWhnVtuuF4eoxvlRR8wDkGtEemwHBdX86bo1LAG54nnoCsLNcp5IrRYPthBttXcWnf/Zx6Fh6K/JNRiUNZceIY80vEHRrmgxik2d45BirYy7FBd4PG/QwuYhv7lykurX+Ifs1mkKvC+LJ4KIoRph+WlWurnFdmampq5OjRo/dt2bJlk/N/jx49cpcsWRL9yCOPHPJ+htKU9Sk6LSJhmIYAEWlAWYNjQsWiIzr76qKAq7r7jsOh1UVg2uqalvpIqILJi1PokUk6Krkx3N8E3vHufi4CtjA/Jnu7HY5GAvtBzdC2G18Hp8NBs6cmwYkTcNttQZvoofhkb9eTfe+6PjlpJMTBq2Ph0TF+qHINAy42QxK8SHdOG3KBhBM21cErCz5l0C8f81O7i/i6a2BxhT5hs/nspAEBbHhamXnAFi/2mO5LhjbAeLElp1bmYM+E839qxjtpk9mX045Pi6ewNNAP4MrVZbLbFXpkJjeG+5vCVN8clcJtfm54fG2TQOytig00Iy2tAfUebsmfue0Jk9O8uEMVSXdJ18G9t8PkD6pcvN3mzZsjunfvXirZs81mU4mJiXmujvFGWV0D3wI+ARqKyIvATcCYMp4zyJgjtm0E1A2DnTl64ITSWWXeATOLSB34PM7ltVJYz6ggB5OfiZZsXmMUj9HWaEdCz4Yweb/XicM52T/0lN7hd433cqGHM+C7bOjXAWKu03anQYPg5pvdHqJQbFNbiVueQ50vf4UaNWDEiBBN9Kro1QcnDRGIvQByPbmfuyIpSYdgTJkCH3ygbY+unrOkRvDLcWq+nUtd1nLhB//m5ho/cdnLP3JT60QGqpK3wZk7lJr+F291eiVvI52EpHhkeaZ20uhd22v1eKW0R2ZiFx+/lqQkHQy/YIGu4j5qlOvPGQYMiKHTa6mo25pA53082nkko377hIWHErl+n5B8HtoJavkamDIdR+9uUJ+6/vbf9Rgg+GPAT8LChM7sZYutOXun7++KAAAgAElEQVQ/rMuIa+PoObAv76a9xebsRF3rzjCgZ3d4e3rFiBB/3dycjScivX5OKchxGEQaDq8PSudaObzfwWsF9IyMjIi//e1vLQ1z/Pz9738/9Nhjjx157LHHApZWAg096CUiopSaAzwBjAf2A9crpT4MtDHBptD1GoHEaHjsfHhmd2jtdg4HzDSTDreLcOut5bQlBDuYvDhO9/viv/HH/fyIpzyRxeluzsHD6+sJ3uGASZM8OkKkqHU8deIenn0yn3A78M9/BlWqc1IYgoIULTwrs31yP8/K8uJ+fiaGAXfcAbVq6fI/L73k+mDDQJSiF1tpHvYTvzCFL2Omsr5VAmO3S0m7DbCetQzgcqQN7hKZucU52T/LaNNJI9onr1QxHxWPIRgu+sW11+qD16/3nMxZQViBIubIMX4bFE/GiMUcy2+OADtPmrfNMKBnIg6Hg3+0m4i0oJWPLSnqR6kx4ONx/o4BfzGEMKBDe50Y/MTahjSO3IkIzNhT/NY53YIrcWWEHIdBRm4EOY6gSBAZGRnhMTExp9PT0wtVmGVZ5JwEKtndDrwtIunAV8BHSinPeXgqgHgS+IglACgBefQ82HcqtHa7efvhm0zoFw1v+D03BZ2S7uf7ixLiJjf2KN365H4OetGYe1S/X50LE5+AJUt0NYQ+feCnn1za4NSmTcT9nkX8apBhw3QOzBBIdSX6nwT8chze3qcnfRfqZSdDb4Bdv8Nb78IVfbV04xOJiVodO3ky/OtfMGAAdO9e+nMXRhLGa3Qv+Ib8mjGk/HkbzXcKf7TWQeaJtU3TEg7e4FWyyUIdx68SN1Dkfg8QRxdw7NPPgA+SzcTxWrLxS403bJj+/ufP18/XnDmun7NtuXDMTrOF6bS5qDaORV24ttqfzLkZXtwOAxthVvIW1sfnM6Pdz5Dpf//hzDAkU30585CWsD2EDvg8BgIhqRHy8zHipmwlp0cN0ia243CPeeyMimcn8MoOGBVLkWfzyjUhaIQXfJDAAL1xWpkVQc/oXF89sT2xevXqyAsuuCAgVaUnAg09uEcp1RV4DqgLTBeRFSIyTkT6mHa8CscpMd3IIK0GCrXdzm6H8WbdMA+5IZ32ukKJI4SUcL83DJ/TRvnsfj7vgA4U619HZ4cJC4O33tIlcVas0NLa6dNFn3c4YM4cWg99iodegsxLO8Ps2SFTK5d0vxed8cKHDPiF7ucFfk72IrocT69eOvTkjjv0hO9cXJSC1avhhRHA/YRJN5o/1oGjMUL9fWDY4eWdFEp3C5nPIhZwI0NhN75NPmdQ4hlYZRahne25CkRhCIa/X4thaG/S6tWLCs+6YmhDuK8JhsCl92dTr9cJWqyqxwWbIFfB10f0rSrAAQbYMFCZeM715YYS/U9q5LNXqt8hGP5gVsEwgNYcQlpvpWXWHBrY9Pyxx7S2kHQ93H8HTJlBfcRvNW65YAj0rh2UhQ4gNTU1om3btpVjsXOilNqilHpDKTUQ+D/gZ2AI8GswGhcMSrget42AeuFeIobLwD+2QZr5HfV2HwpSKkVWCCnR/6RGOrAWfNrZ++R+vlyrYYitXrS4d+umbVY2G6SmaseFe+6Bm27Si8Dw4URv2kudk8Ler6YFzfvSFSXc71mvbVVCMe9UzwQkbIaFwbJl2m61caNOltyzp74H/frpFGvbfwSeBDWeFm3Oo3cfGPkyPDQGWi6Drw6CXRUwhxkoFF3K8JzoZ8AM3H8jVk/23x73zVEpECeNrl3hzjv1+xkzXKszDaMws1GYKOL/s43z4r/g0dFw07vw0lZYm6mYHDOD1x87wT2fd4AjgUl2oO9BYQiKj16pfodg+Esv7UFfv5uiYbdwTt47iYe/PAwopuw1i9s67XahmK8qKRs3boxYtGhRPWfIQYcOHTpmZmaWeTdcJgcVEXkT6ICe7lKAuUqpB8raqGBSwmaTGA1vtYVRGdChli65ESwcDkg3k8j0j9ZqQjeUSJEVYsmulM2qV7ROm+XtOF/cz+122HbKXDzqljw4ORliY/UC9/vv8O67hX9WQMpNF3DT3Cw+toWXqX/eKJ42DNDu3ApYeQKSvTsqBTzH2Gxagrv0Uli5En77Tf930rkzXPIsvBuGIfCXOZCWC02Wwa1TYUF1OPjIRL6t9Q39GMCjPMGTjA6oKSVs1zabTnDwzTHfNjxj9WTfLtYPVa5Tuk1P17GTt9wCc+d6vNf1Owrnv3cPm/cvpduK9hxuCu9fsIe5A75HTuazrMd25kBzP7pdiDN13jOMph3tiTfMItdevvuAQjD8IbkJrMwibMp+LvlXexan7sa4bBbV8rtyupqN5cd1OSTDYQeVQwyqnOKmKpbFixfvDMV5y+qNuRn4DB0X2xGYLSLvKKUml7llQaKEvl6ksDo4XTzv6vxm3gH41gxWvb2JD2q58tmplbBbohBn6rSV2T5N9h4lm4fSYelxXbvuzMVdBC68EDIytGRTty4cOQIipIyIZ/BVU3lJ3gq5ZAsu7HbLM2GKd4/E4u7nPnskFsdmg59/hvHjISVF1747fBiaNNGLwYKD8M5WQBEWBvcugrtHwsrzwXZKse07A66FW+R2wghc+i1ptwV/HJUCnuzDwuD22/Vi50wQPWdOSSm+2LMYnRxHl+hX+X1AW2b0gMs/Uyyvt5+8fxXwwf216LTVoD55vmUrOLMfZ9otjYPBGwPeUEqXrFq6FB5/vHT/e9aGt/dTpxW07tua04/fxxP2MH54Fd7e7eCqXz/jytdeADKJAt/UERYuKdNip5R6p9iPX4jIZOA3tHNvpcCZNgx0jS9JPQlPbYNdufB4i+DYipwemAoY4N2luzwpbrf8iCUkJMVrJ43J3p00nJO928HunAHb1HB/H8PDS2XDV6zFzqRykWyhyGYTS1sSjEStvpqM17RZCXHw+osw+nlo31b/7DdhYbrenQ+Eh8PlT8Pu5x103p1C3pp+JB96kRv/OowyrHWFdkswx4AzBGNFpldHJWeAfUAT/rBh2m47ebJ+Bnbv1g5LTnV3UiOdvm3yPnJ7HyM9eRTXjIGVUwaT1XkFvRZ15G9f38LgNbOASI7r2pkBUeIZCOYY8EZKCvTvD8eOwbx58N572mGp8IT6uxAUHSemsOfJH2n6RiztP/6ca49/TYvU73Hs0UH6x2pX0hjms4SgeAWIyD0i8irwJuBbLp5ywpk27Eke1Tab+FpwV1MYvSN4aaPm7NUemP9XCz6P9zh5lEcwuatrFtrtDKOoLIkXJ42EOFhi+heU2tnb7ZCRq4WDiwOxm5efDaJE/0EHVYMuaOvFSWPI9TDuGejSyYcL5eRA4++h83KYsdt7kJ6LUIjOdWHl/QuY/vFAZs8biO2jgXyaaHDgNw/n8UIpu6XzulO8B1fHd9Z2q/jOPlwoPx/6r9b34M71uk9vvgn33Wf2c6W25w0dqp2WDEMvNkDk6abEP9qQ/T/cwePfXciTvwyhXscHCF/7DN/zEcvr/4ftELB6q4TdMlhjwBU5OdDiJ7grRY+R+Hgt1XXpou3XvXvreMTiTksABQ5qvD6N8/s/Tu24Z3joy6nE/bGFqFMnWXjzRdz5fhi7LqiKicOCR7Bu3hdoleZ56Ji7SkORzcb5C4EW1XWQZjBS8Njt8LTpgWkL95oFPZUUnuVpXub1clHhwRl2O/DZScNjMVNnNYf+dT3aJysDpfr/RhsYUEcXtPUy2ads1JlUUjZ6+JDdDnenQs2VcBDtpDRiOzT/CfI8OJUlNdIZPaYUZeKPr+2gTZMN5BS0of5n75D4TSLHN8BnvaEW9QPyxitlt0xqpFOn+RCC4FP/HQ69uNf9WavyDwLvHYXoZTDjd62yfcGsELBhA3z4IbRurR12Xn8a+AN5dBj1Xv+V7gNz6fhpCrPu/pOpt/Vizv25LG95DZuP3ERz2vqy5XB7D8KLbzB71i756u44T2OgOHl50PFnqLUS/jgN//kTYpbB7L06HGX1am3Dttt1/4cP1wvftLeArfCPCUQ+8TYXDbTTN+Ublp9cyGu93uXOp+9k+MQ85o4ooM5hsWr9lIGyFm9dKCIdlFK/K6WmAYOAF4PTtOBRymYRzJI/D6XrwqzNbLDE+1jsQjzz+chUJ5WPZFe6/6aTxsyDXlV5SsFpV9nfu5uJFW6u5/d9LLX4hJgz7ZbYbDDcNP/YPScZVsoHJ52uq+DfxcJZnD43+wqg3nIt8bjCMHQm/mIbrw9lPpttq7H/tpQmLa5i2sPC6TjADnU573yfOuyCUiWvnA5FK7K8F7RFB1a7vQez/tCL+5nJnXKBv+7Qmo/Ro7VK8957oVkz2LNHOy2tnQcsgyMryevelIyna/Duu00ZP1qxu15PvhjYgf2DjtOW/2AQWSMo/QedD89AS/c+jAGPz0B+PtRfAZvzSyossoHb0uG+zVqdvWmT7n8PszDMr7/C/14BFsGxa1Bh/8fev4cR+dc6hCce4FTrJaTHZ2CcbsioLy6l5R8hdFuuApR1tp8NLBCRDSIyHVhIJcyNWbLUC9pJ5d6mZkb1MjTXboefzNCf6+rr2LJKSIlYM2cIwoA6WvXqJdYI3KgwX9un368ulb7O87nM2eAjlpSrZAvF4i0Bfs0u+erpeHdrssMBV6yDDab0FgGc7AlZvaGJeVAO0GSF+wWvmLNIATrUwGiwlIEdF7MpJ4w1fWHjLEgcD4fYttV7b11Tagw4Y6K8FLQt7pHpUrLJy4O7tuv31YC7G0Bub7irfpGd8a7tuvR39+66yO/OnTBkCNx9N3RLBjpCx/sJ/2Ub/3kukUfu2sfg6kOZdUF37m8B0wfWxXggi7qkB6z7djkG+vkWbwcengG7Hdr8CifNpjULg8zu2qThnF3fOag3BDab7v/y5TB2bLF7EAFshGceJezfs6kx7Qv+M+cJamaHc/cdz9Nh1uv8ZcdDhU+KRWCUNc5usVKqC3AL8D3aM/OKYDQspHz0p1YdDd9aNrvdwxlFE52PcVsprGcw1xRNuuWAS5uNU7LpUcvr8aUG+j+2wYYc6FxDqwT9IIX13Mggfd5yHLslbDZQGOMEeLXbuY01m7MXvjcXy/PCIasPREbqPJ97LoN7TfXu0QJd7snVdQzRPgorM3nc8TDf8g0D5ApuMYZwwqE9yN7ZC+m3wilO5pat/8XslkmNdPV2L6nDnB6ZLqV7ux2aroR89EKX2RveidP9f7cLnOgN1dF/b1BswQ8P1/lT33kHVr0P918Fm27g47SxtGcFbQnnYeMxbmlm8GZHmJ0g9L6oDWHkB/zAuB4Dpme2D2PApVTncMAVa7TaEuCehvD7pbp6x3c94GTvohow9+0o6r/Tacl5Dx5+F0ggu1t1NslT3ER31rX+L7v7fMW+ltk0W9+JkxsK3YosAiQoNjulVKpSaqZSarpSKqDyC6EkngQ+5rMiSWJIA7jHnIgCLXmfnw//M53D7m3sswdmqUm3HChls4FiaYg8q7FKVUCw2+EH06jfp45flZqhYvoPLmw2yU309/b2fp3izQ1uJRu7HUaZ/hJ1Ddh+Ucl7YRgwqT30NSfSDbkwb1/pC5h2OzV1H1kLt1GTWjzHi3SsFUZtU9Om0PnLy0Ip1bFhwEW1TLud93pxpSZ7hwOuTIFj5rPzXhu9yBWnRg043hsigRMKuruoXmC63ztEseXEQV4FkriTRLoW/jmpKRjJ17Gbgh3+9tuJyzFgM0XP2Z7V+W4rIMzZC9+bzk5xEfB2x5Iq/Ro19AYgEi35tVnp5jpaZos63o7vaMdqUmkgNsb0u5lrG3XihvnCng+uZy/he10cXOVwFm/1t4hroImgDRFpfWZaMBHpJiLX+HiO90XkkIi4NH2LyHUikioi60VktYhcEkhbQT/o8SSUHOjOh3L2kcBUmYPStCNCXKSe1HywWxV5Yoa6hl1pXNotQRcz9TbZF88i8eh23e/4SHizbYCtKf8Nqsv+O0XWFe4diN1KNv/YBnvNievNNq5V2IYBX3XVqi2AJ7eXnuwMA9WzFl/3/ZU5Q77iVkaQSFcSa8PjrbRQVB2YsAOoUdPvqgdOStktoZhXagBjYN5++NbMada/DtxynuvP1agBt5t1BDfkuV7wUfx4+Vpeung6R6jBQO4oPT4MgyOogDOogItnwEd1vstMKg4HTNij38eEwepurueAGjXg3+Y4+cOun5szMaX7wws28a79R2oA9/EsyUYyf5kjNOwF2A0MWgZUoftspHjx1hYtWnTu16/fBXl5eeIs3hpIEddAJbsPgHXAGhEZJCIzRWQHsAodtusL04GBHv7+HRCvlEoA/gq8F2BbgaIyKYUD/Y020C9K6+w9TPYuyc+HveZW+/FmPjto6DRhoS3r445SNpukRjoYHHya7AsN9K9dAA82g9XdA5Lqyjvswokz3rLQZgNasgHdMW9OGsUXuvx8WGxWzujvOU4Lmw3Gt9bv9xbAletLXWsXu+n7bTf6rEzkVkYgiC4PV1MP0GsbQmYBULN2QBWawY3d8o02eqFaekwnRXB3rJxRzNbhgOnmotWpOnzRxfMYeKsdNDP/7mLBV4bBpcsSeGjCjfzNUZ3WZQkq9EBpu6UBN5upw7q7r2TjcsMzZ29RceaJrT3b65ObQl/zq/vhaOkNT1IjHAPq0HBJGM8/8lfiieAurkIQwsLg6mXQoBdUI9K7vvUcITU1NXLMmDF7t2zZsmnnzp0b09PTI1atWhXhLN4aSBHXQBe7PCDWXIhqA5ejF6TzlFK3+nICpdQyPASJKqVOKFX4aNWkjOJAKikkc2PRg26zQaz5gP/iR35Zh0OrY9JOaZvVLc38aodAaXVKOVDKZmMYcLupyt16wqs6t9BuZ7PBm+38XuigfGr4ucMZb/mUM94Sivow9YDvk32BQ9vf/rBDcxt84aZmXXFuaVa0sfg2q8TmSqH4vMW3/HXWSwwquK5QfQeQ1ATmJsDjraGODYhu0CiQvhe/VolnwGaDC6rpkfWz+3jtwmK2Tslmzl741rRVXlbX+7Ngs8F407a7twAeKulnsy5pD8M/HsuNC/7CbQv+U9jWYFOq/1DkYPXaXp+8MgG92fmHqVGNi/Q+BxgG3GGaOdJOlZLulCHMuv975t68lOfHzuBuxlC72Piw2aDXJHBQUGVCDzZv3hxx4YUX5gCkpaVVV0rRuXPnPGfx1kCKuAaaQWUycC3wnlJqtoi0Bn5WSgX1yxCRG9Bxew2Bqz187i7gLoAWLVq4/EwX4pnHR8VSJlE0SWXk+V7QdcHBIqeUS2v75XZfOm1T8PB2D1y6+yc3gRkHYGmmDqF4u6Pb8wcjN2Aoa/j50v9Sdc2SGsH7+7Q6bvp+t9lEilcuv2jbQVo5d/RX+zDRgz7nF/HQ9TctDcw4oO+9YbCedTzf41VOd8vn5v8kIZeoQnuq016lFHx3EXTdn+G2spovY8DlM+D0Zv85W0/2LvpTKNmcNhf7CabpqHk131XZtzSFWYd0erklR+GNfKhWDQcOJhqvs/iqn3hg/GDaREWzlpvoxkclJnxfCGgMvNEGfjymv5d/bIPJHbxfaNBG7XRU14DVXX2bA5KbwAf74bss2HKyxHyznnXMrTuX5dNTmTDrH9hGvEe2cVWJ/jfoCofYlu79QsHjr6k035iN1+KtunYrRqSBw2vt1ihy3u/ivXpHRkZGxIgRI1rZ7XY5ePBgtUWLFm2LiYlxOOvaBVLfLlDJLh94XkSyRSQfeAbIEZF9IuKTZOcLSqlPlFLtgeuBsR4+92+lVHelVPcGDRr4foGJsaYaJ9Pjzr4Qu11PigZa1/+Wf0UhN5DKs4xmA6l+HecL3u6BS5uNYehUXwBbT7pV5ZUpZVQ54csz4NJuN8KUbr/NdKvOdk729gJTveD8/if5MDE6sdngUlOVtVRfS6FYytecqHaSe34YzsBHYl16B4vo+nbkuffGDKj/oMdAp2qwMbeUxOWqHSw9Bptztc02o4fvEr5haHXneeHae/Fqbfz6kAUsYgGTJjxN75VxRH/VnHaMIwpfs04XEdAYsNngYlM7uMV9Ud8STipOT2Z3tlpXFJfuvs8qnG8Uim/4imW91zF21kP8/W9X0OnrZ0v1X6Rs3rihJMeBkXGSiBxHcBweixdv3bFjR9rLL7/8+9ixY8ucgzFQyW4c8CjwoVLKLiI1lFJ5ItIGHXs3q6wNK45SapmIXCAi9ZVSARWjc7q8f8QSEkjUv7TZILYGLEXn6PNkewEt/XxjJj7+0gf11Rm4lC7LiVI5Mp334OK6Og7ou2w92btQyZyZMqkyL3qeKJEf0dn/4tLtiky3KinnZJc5sBG0MrRHr79JCd5spz14N5+Gn4+x/paDvM7L1JE6DFVDMRwSnKw+biidI9OMq7i8HqTth5+yPEp34TaQAXXh/I6B9d9mg2vraaeo/adQ9tPssmmP1poXN0EQsppuZasaRS1p57dk5w23Y+CSevDOYfj+hF6EXMwDJfKkzmhKwvzwoqTyvpLcRNs6v82G6XshuTHrjfW8xksYhkHbxp04Gb+TTZe8SATxQe+/v/gigYFZu/U4ET3rkBuMknarV6+OjI2NLVRRduvWLWfSpEllUuFD4JLdJmCfudBdDiSbv88DVpa1UQAi0kZET6si0hUdsfJnWc5ZWNOquM7eGR/3U6aHwN/CRunXtpF+D3SFIpUUuhBf7va64m0oZbNIbqyrqoNWs7iYbH1OmeTl2uVVsNZTG0r13zD0hgfgx+Nu7TaFOSK7GFr9GUj2HZsN/s90iPgli632NDLJ5CEeI+F4JzNkOHSeqqVizZxMjNWu805VngsKNzyGgRoaYP9Bqz07VYe0PHaN/oE3mEA0tWl3yUVwf1Oi/hlNx69KSzbBwv0YiNK3fvpet2OgME9qXIDPgGHACHMz9a3eXG5lC9lkcycjGXhsAFGpF9Dxp3+FrP+hwBDoXTc4Cx3o4q3t2rXLBXA4HEybNq1enz59sjZt2lRt6NCh5w8cOLB1QO0MsD3PAI1F5AEgUin1gfn7OOAJX04gIvOAFUA7EdkjIn8zE0rfY37kRmCjiKwH3gaSijms+I2zptWzPF3kpAJ6t+Uc6G1+9WyknhirPREnxvp9/YoIJj8TlzYLw4C2Zrmj77JgjivXcB9SJnmhPAvWusNtmrI320GnGrAxT0vvLvApR6QvTIyFuBqoDXk4JvyOgUFLWiErs4tSuIVIunMZawZ6Eb7E6S14zOUYCMaGp/Bal+lUZcf2HyRLZfIwj5NgdIOeUYhDOHL0C7JDoOoHD2NgRDO92fg2WzvguCAoz0ByY206ARwrMvmVFRgY9OJiDARRwpFjX4as/2cDGzdujFiwYEH9Dh06dOzcuXOHvLw84/XXX9/XsWPH/IULF+4O9LwBLXZKswA4DNwnIgtE5J9o6csn90SlVLJSqolSKlwpdZ5SappS6h1n2SCl1MtKqU5KqQSlVC+l1M+BtNWJINzEUOa7clJ53IwR+uO0jiNzRxk8ESsqmLo4bh1kJsZCczOh45R9IalUqQvWni63sj6ucGmzAf19Nqmu3/943O1i4zVHpC/YbNCnLqA4sf0Qd9lHMoSkogri32QGrxqHC9w+AxebIVxpeW6lu7JueArpXRsHDg4d+J3qOeF6sUcwa93QlKSQSTZu+5/cGDqaEv5M1x7tQem/YcCt2ua3LjKFdxyTuYt79TNgSoptVt6HchRU6FxRkSxevHjnwYMHUzdv3rxp06ZNm2fPnv17VFRUmXeAgQaVOz0jXwLeAN4CDqFThc0sa6PKnVuawez2cH8THUcWZCoyvqw4peKMnNhskHER9I6C/fmQcqLUsW6zSJxFuIw1c/J5nHa6SMuDDw+XOtZrjkh/mBjLwa45dF/VgeFPX4KBUVRBHEJqt3P7DCQ3gaFmIeym7p0ugmKvTW5Cym0H6Pu/7sx48nmGOHRhVUzptmDmbrIdodGAuIy3BL3QjGqhpbsRrs1DIjrTWZnHgC0MBwpjaSa1TkbSk976GUhqhNzflOpT7OxZ8ApZZ35HFmUiUDXmBPN1o1LqW6XUL0qpd5VS9yqlLg1W44KNW1WiYWij9KQOAUltvlz3KR7lpQpU4RWnlN0StFfZz93hv/G65t8ZuMwi4TcV79ni0mYDuv9re8D8Ti4dD8pUtfsMHDaDmaP/R0rCNmJ21Cy9uPlQiSBQ3PbfMGBuZ3ipNTzZ0v3xwQhBMYR9Q+28e/enXDS3NcY8U5J9Ixb6RVP/m3hs83JCItm4jLd0ktwU5nXSry4I2oYnqREr/7mPuNTW/POnh7RUB4VVMJRDUT/6qrPKblceHDhwIOzmm29ukZaWFjlq1Ci/64oFutg1E5FRwGYRGSciXmMxqjJahWevUBUeaDXeS7zGs4wuvbMHPaMnRLncvpdlsq/INGln4rG8kOHd8SAYks2HLOBfN7xJ5MkatFxUsyjkwY+iqoHisf9hYXqhc1NJJliSzXrW8fcrH6bVH01peCxGe8EC2GxIO50/ImfGb2Q51gV+ETe4jLd04uX7D9aGRxnC4RpHeHjSJDqvb62lumJtyE7czqY+z1Rpu50rGjduXDB37tzf//jjj43jx4/3IVasJIEuduuBE0A0cBmwS0TWisgHIvJwgOcMOaUSQpcrFa/7E4QbGcLzjCcugF1joPF2FZkm7UzKEthf6H5fhgVPodjNLpQoEkjQS45zsk9qBCPNoqpe6uwFilu7pQ8kxMHrL8DTY8vmpJHOVrKMLPZdmw9IyRRavaJRQNRvzYma56drv4+U9Rkoa8xpCuv5z6UL+PzqlbRKbVjSISipEVFPDqJdzZcsyS7IBOqgcplSapJS6m6l1MVKqYbAdeh6dsHXAwYJQehCPKmklHvVgYp0uS+OW5uNDxS633f277jK4JzipCz9PzPeMBDWs47JTOQ5eZHYxmZQ+tacoowaYeaQXOW9zl4geLRbeju2uDHhWQ0AACAASURBVPt9gDXDlfkvmtp0G3E1ElcD0vLhYdMpJrkJ+YMKqH68LvkrdgZ2ES9UxBgojgMHay7ZyqJn3qDlgpolk1AbBtlJB9lqjLEkuyATlIh3AKXUH0qpL5VSE7x/uuIolSOzHKhYibIkbm02PlA21+uKl2yduIy39IFguN9vZQtHOEwLzsd4s73OxPJ9VpHa8g0zvOUN/8NbfKXingG92D/CA9iwYdjC4XEzrVePKP1qGFRL6sDxhG1U69UqsIt4oSL7r1B8yzdkSRYHhpzWm48eNUt8JoouAWeRsXBP0Ba7s4WKyGJSqsRQBVJWKVMp/93vK5Nk6zbe0kfK6n4uxf5hs8GIM7IglSG8xZ82BPp9KOWmkKuPbCOdLDJ5kEf05s9mFu2zFdkJs5MPs+6nl8hOLu0VGwwqYgw4SSWFd5hMbWrTIbut2feS33U2qWxltCXZBZlKq3IMFc6Fp6rilDIDWexLpExqq38O9TWDjTPesi3tKsRmM4QkDAxuxHS3dzpE+Jt6qgyU5fsoq93SmabsRoboxcZF/6ONRLrV+i/RIXpeKmIMOOlCPB/zGQAtk7qAcaTUdx9NPN34KGT9r6pUucWuqlOWxV4Eht4A7WL9s1lUtg1GWdqTEAefLQjcZmNgFLmaQ5EHYDlS2fsvSEjzQlbEGCh+7cKcnAYuv/tQ97+qYi12Fn4h4v9u9lzC6n/V7j9Y9+BspcrZ7CwsLCwszl5mzpxZp/irr1iLnYWFhYVFpWLChAn1GzRo0KV9+/YdW7Ro0blfv34X5OXlSVpaWvWlS5dGO1/9Oae12FlYWFhYVCpSU1Mjx4wZs3fLli2bdu7cuTE9PT1i1apVEUuWLIl+5JFHDjlf/TnnOWezW7NmzRER8bUMRH0goGKwFYDPZdH9uAdVvf9wDt6Dqt5/sMbA2c7mzZsjRo4ceRggLS2tulKKzp07561atSoyMTEx77vvvquVmJiY5+08xZEylIg76xGR1Uqp7hXdDl8IRVurev9Ded5QYD0DVv9D3daUlJRd8fHxRwCW/ZXmRzfiNe+xUlCQgxEWicNbSEpMZ3L6vO+9AnqdOnUSGjZsmG+32+XgwYPVFi1atO3KK68sXY7FfT/qx8fHtyz+O0uNaWFhYWERMAU5GFkZRBTkBGc9ycjICI+JiTmdnp6+aceOHWkvv/zy72PHjm3i/UjPnHNqTAsLCwuLsuOLBAagHHBoJRENe5IrQVjuVq9eHRkbG1uoouzWrVvOpEmTyhyMWtUlu39XdAP8IBRtrer9D+V5Q4H1DJwd5wwVlbKtYkCj3sFZ6ABSU1Mj2rVrlwvgcDiYNm1avT59+mTNmjWrzrBhw87v27fvBR9//LFfnphQxW12FhYWFhZFFLfZVRTXXnttq19//TUqJibGLiKqa9euJ6dOnbonKirKAXD48OGw++6777yFCxe6dUByZbOz1JgWFhYWFpWGxYsXe6ztNHr06CYPPvig31nCq7oa08LCwsLiLMDhcDBy5MhmV199deYll1yS4+/xlmRnYWFhYVHpGTduXMNly5ZFZ2ZmhqWnp9d44okn/JLuLJudhYWFhQVQOWx2wcCKs7OwsLCwqJKcc2rM+vXrq5YtW1Z0M4LOmjVrjiilfKrweS7eg6ref/D9HlTG/h8lj50UJcCIIpxYol1WCz9KHvvI5RQFAEQTThuiWbtm7Vn7DCgUm8kkFzsAkdhoRA1iqFHqszmcZg85ZHMagGiqEUu0X2PAojTn3GLXsmVLVq9eXdHNCDp+5Do8J+9BVe8/+H4PKlv/13KYi/kccBBHHd7nYrrR0OVC58SOnatZyjfsJwuwUwfkxrPyGSiggEv5glwOE47wHAk8RTyGB8WaAwdz2cYsdrCZY7zHFSRKfZ/7b1Gac26xs7CwqFykk8UpHFxBEz6jPzYfph0bNr7kCsaznudYTyrH4YKGrcqhuUGlgAKuYSkrOEw48AEXczOxHhd60BXdh9OOW2hLCkeJJ6Z8GnwOYy12FhYWQUehWM//t3fm4VUU6cL/vZ1DSMISloRVkGFfDdHoIO4C4zaKzL2IUVFH/fx0xkEFnAH1+nmv4z6O4jp37lURlc1BruB2B0VFlIhhCSSsYdUAWSAkQNaTru+P6kMOIcnpk+ScbPV7nvN0nz7d1VV1uqvqfeut9z0MwGT6IQiT6VerNFMVC4uHSeQKevMCW1gYE9kuVPkNBTY2U1nF5xzgCnrxJGdzNvEBOzp/BGE0XUOYy9aD6ewMBkODk8YRruULAD5mAlPoX6d0BCGJ7synGwt35+5qyDyGmg/YyyL2kMwveJeLiSCisbPUqjGdncFgaHAS6MJyxp/cry+CwInS4nonFEYm0w+4NGiJ1hAaAv4DIvKqiIwNR2YMBkPLQBASiSORuKDUdi0JC4sp9DcdXQiYN29ep2CvcfMv7AReEJG9IvKsiIwOPmsGg8FgMATHrbfe2rdXr16j/I9lZGS0XbFiRdBRDwJ2dkqpOUqp84FLgCPA2yKyVUQeE5HBwd7QYDAYDIZAbNu2LXLNmjUdysvLJT8//2RftXz58o7Tp0/PCTY91/K1UmqfUupZpVQicBMwCdga7A0NBoPBYAjE7Nmzez300EMHBg4cWLJ+/fpo33GPx6MSExNLaru2OlwbqIhIG+BK4EZgHPAN8O/B3tBgMBgMzYA77uhDenpMwPOUgqIii5gYGwkwPztyZBFvvRUwAnpqamrU9u3bo5csWbJ39erVHTZt2hQ9bty4EwAzZ86sk+/OgJ2diEwAkoFrgLXAQuBupdSJutzQYDAYDC2IoiKLzMxoBg4spl07uyGSnDVrVu/HHnvsgGVZDBs2rDg9PT068FW140ayexiYD8xUSh2p7w0NBoPB0AxwIYEBYNuQkhLNmDHFWPW3PF25cmW71atXx27dujVmxowZlJaWWkOGDAk6fl1V3BioXKaU+i8gX0RuEZHHAESkr4icV98MGAwGg6EZY1kwdmyDdHQAjzzySO9FixbtzMrK2pyVlbV5w4YNGTt37owG2LJlS+QNN9xw5pVXXhm0l4Jgcvc6cD5apQlwDHgt2BsaDAaDwVAdS5cu7VhWViYTJ0485jvWp08fb3FxsZWdnR0xfPjwssWLF9fJIXYwnd0vlVK/B0oAlFL5QGRdbhosInKliGwXkUwRmRWOexoMBoMhvEyaNKlw3bp126seP3bs2Mbu3btX1CftYDq7chGJABSAiMQDDTIZWRvOPV8DrgKGA8kiMjzU9zUYDAZDyyGYzu5lYCnQTUSeBFYDT4UkV6dyHpCplNqtlCpDW4NODMN9DQaDwdCEOHToUMRNN93UNyMjI2b27Nk9grnW9To7pdT7IrIOvcZOgOuVUuFYVN4b8LcK+hn4ZRjuazAYDIYmRI8ePSrmz5+/vy7XBhX1QCm1DdhWlxvVg+pWKapTThC5G7gboG/fvnjxMpP7qcDmDM5gJrNafHiNqnXwHE/xE/s5Qj5zeY82tAmYhtcL98+CnFzo2hki2sCcp8HT2LExlKId1LrOpmr55/Men7CMrsQxlgu4geSADnm9XnhgNnSLg/StEB8Hv70ZzhlNwLWybqiogE3PQ69xcHQ7HFgNpXkQHQ8nciDCgsguUJIHlkBkVyjOBU8ERHZVdCK6s9vye/EynWkc4TAjGMUfmR3wHbBtWLAEvk+Bg9nQxgPXXgU3TabehnZKQc462P42nMiGCIEev/ay33qemMgzkW+mUHooAonQdVB6WL/kUXG6jhBoT1yN5a+uDhbwPsv5CFC8w/yA74BSkLoBnnsZunSCI/nQLR7mPNMw70B5OXx1C7TpDCW5oKSCigvfpbztfjxfzCa2WxuKsnVdxyXC6FkQ0bKbrbAiSqnAZwEikgQ8ApyJ7iQFUEqps0KXPRCR84HHlVJXON9no2/8dHXnJyUlqUtSL+B1Xq48xnkkcjY96Y2FxQSuIJGzG90bu69BOkweIzmr1gZJRNYppZLcpNsn6QyVl5p18ntvzuBKriaPPAToTg/+wpzTIkbfPwte/s9T0xo5DAYPALF0A1BQqJg5bgtJa95ADmRBZCRcey3cdFPQLaKNzSLms4z/IQKLt73v0GbaDMjNha5dIS9Pt0BA4ocfejcoFbjHBgYkDVAHUnefcuxsziORRHLJwYOH6TzE2SSd8gxUV35L4NILdR3k5EH3btA7vpxZG24homNHOHoU4uPh5ZcDtogbn4XUWeBpB95qXTIooJwupFPAEGwiiOQwFXQkmj08w7lsUKUBH9rq3oFhjGAQQ+hGPHnkkUAiD1UZBC76EJLvOlnlJzk3EfqeoTv/vMMQ38XLy0XT8BzJhS5dID8fHnoIkpJqHBUc3gjLLwbvscpjEeeuofPiGymc9hIVu/tTsXMwlEUBFXQnhXIi6cw+dvEbIII/k8RelerqpR2QNEAdTN2DcsbF/enPWSQiQGe6UsBR3q4yCNy4GS68Ck5U+W9GDIVhg/VgoFs85B2BP06DpMTgBkErb4LdCyq/exI20OHlaajj7Tj+xL+hCjtQkdUDCroykMV0Yi1d2+bwU8LdFPS+iKlL++fkqb3d3d8xeNLS0vYmJCTUyUNJUyItLS0uISGhn/+xYDq77cBDwGb8DFOUUnUyA3WLiHiAHWj1aRbwI3CTUiqjuvOTkpJUSmoKM7mf7Wzna7487ZzOdOYzVpJAeAI4aEnzAbrTA4XNIQ6SSy472MYWKovxBM8wkz9Vm0YwnV1SUpK6IfU3/GTv54ujn7G38/7T5OPfMY0XmHNqPv0kO9uGjz7T0sipKDrb+aw8cTmjK9IqD597LvTrB5066cbPtiEuDo4cwWuX837SVkpLCrBVBW2jYslpX8TWXvksv6IEbxuwFfz7c+2Y/nj1jnmSgFSlXDUtSUlJanrqA3zCMvaxj3WsPe2cdrTjS1af8gxUley2bIeMavQYvyldzAfFU06VE0eOhKFDnRaxG2Q7Q/TOnXWHCFR06srOb2NZ9etPsEvbUrKjB56KXrSJiaHs+AnalZcw9KdUzjq4ly1ndqfCEs7avo+NgxNI2JHGiKjI4/uKSzu4KX9KagrTmcauip2sqvgKb5uK056Bqs9bVcnuwCH4Yd3p6ffwZrHr+ABiKK08aFkwaZJu/ZXSneCRI/p7586ow0c4ejiWvTnns7HDFk7EHsfTfwdthxTT4bdriYiIwLt5GKU/dKXzl8cY98leLAQPbdjaoy0fXzqQ5z77oeDw0WOuQrskJSWpGakPspyP2EAqe9lz2jk3kMw7zD/5vapkt2sPfLVa10tV2sXA6s9h9KjTfzsFpWD9evjsM8rXp7Hx+0s40quEnZ22gUdoe9YB2gw4RKeb01FSgbe8gt6/+iVj1u0jm+F0Zxu7uJgVfW7i2dxp7CveE9IRuunsABFZrZS6MBQZc3Hvq4GXgAjgLaXUkzWdm5SUpFJTUwGwvWVs+LelzHviG5RHNZpkN4P7K0fZCiLLPJS18YKAiMVlXM6FXNJgkt3JOliUTfmNaTy08L/x/ktn8jy1S3b+KAXrNsLc93XDV6Nkd+gQpKTUmp8SD3w5AUZlwHcXwNjv4d6/w8g0mLQEIsuh41EYuBskIgIuuwwGDKiU7OLj6f2f//lzllJ9gio/UOEt5YPXX+Sj36XS1RNXq2RXFV/nn50D8V1rkOz27oUvvzxdHKqBJb+Bc1PhsyuhrA2cuR/iDkPvLPi5Nzz2JBS1g23DwRsBnQ5DcUeILIJDV1Jhb1ABFWr+5VdP72L3377n/Ue/IePOfLpZNUt2VbFtmP8P+OhjLXPWKNl98QXs3l1jOv4oYMcgOByvy5zVG3YNh96FIN3gnE9g03mQ2xN2doPb/9qW9OHCgb4l3JUGJZvcD3h8deD1lrDohb/w8Yz1iMeqUbI7La/OO/DWe5CbV0fJbuNGuPxyPQB0yr9tCOwcrP93Jfpg11g40EHXzZAdMOLrtnw0pZR/ebcdn0woZshui7vSvWwqdVf+umI6O0BExqEXlH8JlUM6pdSHDZjHeuP/kHP/dng5C/pGws7ztMqtEfCX7OJ/jOKiqfFsHbiH1c/8TPLIezibcwJ2vHXq7GwbJvwIK0/APfHwRqBhaB2wbViwAJYu1aP7ICS7Yinh4sLRXLo+FsvyaHXgHXfAOeec1orUqfwAf9gKrx6EUVGw/ryGn4BUCtatg7ffhpycWiU7unalPCeLd0ZvwPK0pehELm08UUTFxFF0IgdPRCT7Lx/IzgviEMuiK13JJRcPEXQhjlfb/22LOq5GBFX++QfhZseObN5gmHpGw5YfnFHB/VoNcPhwtZIdR7SnwYr4rnwW8zVZ7Y7SNqYTJUX5iEBkVCwlRUeI6XaC8igoE0VeaXv67o3grGXH8bap4NdZZQU5+cq1ZHeyDu7bCq8dhAmd4PPR9Z+ADAY/yY60NJRdwe5+ii9ivkMsITKqI2Ulx6hQ5cS060ZRyS7adPMibdpRti+CNlYMl75dwC9+asfA2FzvnqPuVPl1pSV3dsG8+b8FhgJtqFRjKqBJdXan8MIAWJoD+8sgaT1sPC+8D7qDBw8v8SqUlcF1KXDIy5Auw7h+eBLBrf4IEsvSViYAn+TrRqmhG3vLgptv1p8a8AC3Nexd3TMmVnd2m0vggZ3w6rCGTV9Ez1UlueqHaQPcVcdbvXrib8VBX3Rjd3hmjy7/c/vh5l4N/w54PPCaO2dKEcCvg0nbtiEtGVasIoac0sAXVINvQL/iKCw4pOsgXIjowds55+ivwADn45q82fD6XNoVeMpCkMNWQzBPfYJSKkkpdZtS6rfO546Q5awh8HjgyX56f3MRLDjYeHmxbThvPRzyQg8PfJMYno7345HQ2wM/eeH+HaG/X1MjuQeMd6a5th6vfgKmJWNZ8FBfvZ9e0rjvQF2wLBj0CwDaIm3rlMacwTDcufS7/IbKWfiY8wSMGkIUBA63Y6iRYFrblGbpueTm3jAhVu+vKWycPNi2VqelFcGISNg3BtqEVBtRSWQkXB+v9789qqW71oRlwe299f7KY43T2CsFGwr1x+W0QYNycy+Y0FHvf18Q/vvXlzl/hl9dTHukY52u93jgUmfVwreFjfMOKAUbj9Xt//d44E/3NXyemjHz5s3r5L91QzCd3YXARsdH5SYR2Swim4LNZNixLLjNWWivVOOM7N8/AK9n6/3L4sI/d/jSID1n5VPltTaSe+j5GoB3DoX/GUg7DldvhLHr4Mld4b+/ZcEgRyj4tqD5DXg8Hrj9RqT6NbfumDMYRkZBenHjvANpx+FfNuttXbAsVJX1xa2BW2+9tW+vXr1OMTbIyMhou2LFio6+rdu0AnZ2InK+iAg6Svkg4FfAtWjV+7XBZb2RSO4J9/WEvx2CRdnhvbdtw/OOA5g+beDFgeG9P+jG4kJHut12onWq8m5zlietKNDzNuEkoT0Mj4ISBf+2H7p8DQNXwmVr9DxuOJgzGEZEwebixlFnl5XB5T/A8G8hbiV0Xwnj1rov/5SJ7KPCnclndXg8cFEjvgMJ7eG9QdpgLNYp/zv73edjyvXsp+L09RMtmG3btkWuWbOmQ3l5ueTn55/sq5YvX95x+vTpOb6t2/TcSHa3AeuAZ4HLgFKl1D7fJ+gSNAaWBec7I/vvjob3QV9wUDcwI9pC5i8bzx3JBU75vyxsfvM2DUFyTxjnzN2Fe96muBi+8RvRFwC7gK+LoX9KeCStU1R5YVZnl5VB9+/gqxOwtRwOAznAyuPQ02WHb1nkoer3x/negZWF4R/wFBfDBZshDyhEl//2TLhvi7v2qCHK38yYPXt2r4ceeujAwIEDS9avX3/Sg5LH41GJiYklvq3b9AK2vEqpewBEZCg68sBcEYkFvgI+B75TStUr9EJYGBoDMaJNkC/opBu/UOP1wrOOG7fLujba0gdAl3fuAfjimJ63ubm3u+tsGz5YBpOvaxRL1gbDsmBIO/jyWKUqLxwDj5ISiE2Bmt6QLMdw6LUwTIe/NAhW5YfOMrU6vF7doR91NHCd0CaZJcAJ4EgFJK2Djb8M/fPl/w68c0Crt8PxTPuegeqUkG/k6IF4KJaE1JP/yx19MkgPaBSjUBRTZEUTYwdaQjWCkUX/SeAI6KmpqVHbt2+PXrJkyd7Vq1d32LRpU/S4ceNOAMycOTPPf+sW1/+0UmqbUupFpdSVwOXoqAeTgR+CuWGjkdgRZvf1OTkL/f1sG65O01LdqKjGUV/6Y1kwuJ3e31HsXrr9YBnc9Ht48LHmr/6cMxhGRWurxAczQ38/rxd6/wBedAN/7DwY50wxDIuEns7r91Fe+KQ7nzo7HP+lbcNVabpDBz1nlnsx5F0OhZfC5Y6kvbk4PNoGf2OlLwrDM6VRUgKdvq98Bv67nzZSg8oZyPdca+KaJMUUWbvIjC6mqMFGDrNmzer92GOPHbAsi2HDhhWnp6fX6hvXDXUa2iqlioFPnU/zQASujIM5BypdGjWEd9+aWJSt54cALurUBLwpoxv7ncXOvNVBd9Ld5Otg9Y/w6tsw5hxInhT6fIYKjwdm9ILbd0FSvd+d2vENdo44It2bA6B9e/htL61K3loGd3WBfx7X60Bn7II5Q0KbJ9BajTeyIbMk9NLt+1nwhfMO+C/qt21YnANP94c7tkFGaXDahvqQ3ANWH4G/5QQ3wFAK0jIgYYT7dsM32ClFh7kuGAtRUbC+BDIOaUlvZDQsH1mHgoQeNxIYaB+3a0mJPo8xxYGcrbth5cqV7VavXh27devWmBkzZlBaWmoNGTKkqL7pus6ZiMwRkX+KyP+KyHMiEh7Hkg1JYkeYMwhm7tJmwKHCtmG1o14fH6s7maaAxwODnEZ+7kHXcwWMSdTn2k1fWx2QVGdd9gsH3Td2SsHG9OA0AgsOVg52JnSCqY6nM3/L0B+KYHsSTOutHSCEg+SeehnCioLgrBKDrQPbhud/1vt92pzqveaDXEjeApekwZRuWsIZ69qCvH5YVqXqcl6Oewl3Yzr8+la9dcsDOysHO38foDs6cIyFnHV/F8c27vRGA2BhMYaxDdLRATzyyCO9Fy1atDMrK2tzVlbW5g0bNmTs3LkzGuDdd9/tdOONN545bty4AR9++GFQS1GCyd1W4HlgDnp69T0RaV6LP0RgSAwoG/73cOhUOYuy4fVD+iW+o1fTkOp8jO2k8xWUGsdyRrPNeM7Ox4sDtSpzc5F7VWZaBvzmTnjuVXfPjG3rJQ6gJZpPz6psYC1Lf/ctBZnuSHThekYsCwZGV+bTLRvT4ddT3Tf2Cw5o9SRoxw7+5ZscD1d01vN2H+TBu0P0ICBcnHwHCty/A7YNRUGo/71e+Oqw3r+8XeVgB3RdXOIYC32UFz6L3GbA0qVLO5aVlcnEiRNPSiN9+vTxFhcXW9nZ2RFTp049unDhwn0LFy7cu3DhwlpDPlUlmDm7vymlViilPlVK/QXtiP7/BnOzJsHoDvCHM+CRPaHR2du2VpMA/K4nTAlpRI7gSe4B9zh5cq3GsZ0RfTOfswPd0DzkNDwVFe4ar4QR8H9ugtlPw6L/CXz+gkN+KuzY0zsyjwd6OCP7j4+Ef92bzyoxs9j9vW0bikrcd/ZPO4ZZI6JOV09aFnw8yhl0FGtnD+E0fkruAfc6navb8m/PhKMFeuuG+3fAFqcTu6PP6eWbMxj6ROr5zKT1zX8+vIGYNGlS4bp167ZXPX7s2LGN3bt3P6laevjhh3tOmzYtN5i0g37CROQeEfkLWsJrJJck9UAE+rbV7XYoHrAFBysXkF/QqelZMPqrcd7NdanKjNDXpGxoGS+lx/H0/0a2O8MIEejrdJBuiu9b2jC8bc3zcB+Pcty4lYd/kXNyT+1VaEWh+3tv3w0FhXobiAUHIcOxCP9TNQ09OMYyjhbq6yodfkkJjE7RW/Sc0AIyIa59UCP5GvHPj1tVplh6wLfGRcfk9cIqx/n3hE7VS60eDzx5pt5vbFeGzQjbtrn33nt7X3PNNQUXXnhhUPN4dWmJP0WrNM8Aqg2g2mz4vqBhG2/bhrmO+mp8bNOT6nyc7zQyK466k26nTITf3QavvuVOsmnqTOmu/x8Iwn2WS+nW64XVzhjw0s41qycjI+G6rno/3IucLaty7nZ7kcvG3pmrW/Nj7efbNrzjPFMTOtVudOKTMDNKtTs924Z39kHs99q13piN2NhMYw03sQrO7No/cEZdEqwqc8pEuO+38PrcwO/A/Tu0xe+ItqeqsKtyc+9Kv621PIcVVPAsaTX+3pp46qmnuq1atarj4sWLOz/33HPxwVwbjIHKYhEZppTar5R6E+09pca4ck0ay9IP+uuHGnZEteCQfnkEuL1n05PqfCT3hN/10Pl0Y3RiWdoSE3SU1eaOZcHtzkBEcNnRuJy3fGCns9wkOrBhkv9C//edyPJVpJqQMWew7vC/KHD3Dky53hnwvF17Y7/goB5EgfZaU9s7kNxT1xPAvFy4Y522lC0D2kJFyiiS+YrXcLRaR4uPuCqbG5J76HcAgngHkhyXg7WcZ9va4hn0vFxtc7H+y4GqmbtTKFLJ4SI+YRbroFenJjp6Dh+PPvpoTkZGxtb58+fv/+Mf/xgyNeZ7wCLHJ+ZcYDHNdRKnTiN7F/jUV+M7hnfCPVh8DZAC5ma3DNVksFiOKvM1lwMeC93QpaTWXF+2DdsdTykXdQxsdOLf2D+wG44fh/g1J6WakOLxwCDHOvC7o4HPdzvg8b1P4zsGdtxgWZB6DrQXKALecWwSIqH86LmcHfW/LEY7afoVPWFXTsO5y7KsyiUEz+53OXenHOm+lvL7LzkaGxs4yTmDtbVqlhd+vfnk4QoquI/vGMsnrCGPscTDgaNh9nXYsgjGQGWZUuos4GZgJfAxcEWoMhZSLAtuc17EnSUNYyBg23rtEminu01VqvPhexHdqnF85UlpIZPpU7oHFw1jyvVw3x3w5bFb/gAAGMVJREFU+juw6KPqz1lwUHvnmBDrbs2cr7HvEqFN1DusheNKB3JJCcPKHt9/mlnq8j8Vp4OoYZ2ZbVf+dLtL7ySRkZB7vl6HBngj4b5C+EXUMjZxlDZY/J4hfMqvXOQvSF4c5BjJuHQyYDmFS1lXc315HSnx3u7uvDR5PPDrLnp/y3Eqyop5kvWMZTmvs4NyFGfRiW+4KnBahlpx4wg6RkS6OftxwAVopz8LlVLNd+l/cg8t3a04Wn/HuLYN07brjmNCp6azrq42glVlTpnoNPZza27sq+L16mjxTdHLvmXB1G563000DDfrDX1SzaBo90sJIiPhhq6nHvvb4Mo1WaHkpUH6ef3iqEtfkQEkmwUHtaQMlZKzm1Sj2vJj4UBen1HI0GPZvNZ2L1mU0odovuMqXmEsEbhPzzUeD8xw5hSTXISKm3K9M29Xw4DHtmGeoyUQXA941cuDKTxTsT6uiM8mreBRtZG1aI3tBHqyjuvw1M3/h8GPgP+GUqoImCMi64EDwL1APPCqiITBwWSIsKxKNY7bSfqaWJStfW4qYEgQDV1jYll6pKpwp8YJdnG5bcO1m+HlLO0dpCni+5/eOOSysa9l3s62tRu2uvDKULirM/TxwFv9wxdJ2+OBAZH6GVjtYjoskGTjk5AndArKOCuNI1zW9mt+/5cj7IosZiDtmUxfdvGvnEt3AvlbrBepjkHfC1ku34Gza34HFmVryR5grHvD0TRPIRP+kc0lqw7xH7OOElkC4+nOe1zI51xhOroGImAtikgfYBSQAJzhi3QgIjHAG+ioCM2TCzrD37KdSACH6t7I+F6Syzto1Uhz4cVB8E1B5QLrV4YGuCCIxeULDsLn+VqlFy7vIMEypTu8dUBL5HMPBnYMbKHLX90pCw7qdOriMcfjgf9KDO6ahiLCkZhWHwvsPmzK9fB9qpZsxiZB8m8qf7NtvW4RYGp8UGr8BLrwFVfwHJuwEN7jEtoQpuDGLw6Eb/K1UZEr59i1vAO+duCebkHN2SfQhdfOvowPX0/hvI+gs7Tloo8vx4oMg3TfinDzRCYCu5VSFVVC+nSkuc7Z+fB33ZRSR0MVrxeedyzphndoHlKdD48HZjjCuStfkS4m6OFUDyKDY5punViWnlsCd3OXU66H917V+1UlG58Kc3Azkex9+AL7phfDgwHW3FkWjD3X+VJF2lpwUPubhKDLLwjn0p0PmMAixoevowOd1wva6/0KF+r2mgY8/u2AxxNUZy8ISVZPnrr7Gq7f0IVLVkRhTdzq+nqDO9z8Ix3R8eyqcjewoGGzE2YsSzth/VUsjOlQN1XmtO1aMhrRtvEjG9SFYHxFWlKp/qwNfw8i5wflvi78JPfUxgRQaVxQE74G7JY/nDpnY9uwo95+ahsHjwce6qv3K+zA78Dka/W81eQqcZtPWmE24fWlNXGhM2e6qzzwOzD5Oqf81516/MFM3Q6Miql7OxAZCQfOhys7w0cj6paGoUbcdHa+ReRVyQOWNGx2GoGl+fDPApi6I/iAjl4vLHf83/WKal4jeh8vDtSe1zcXBTbUmXI9vPcKILU3iqv9fAKGI25gfTjFo4wLbxq2XfnxseCg9jU6vkPzME6qiidCC2pvZAd+Bz5YrtfafbC88pjXq9eWNfX1pTURjLHako91+Wc8fuoz4NOMzOhZv3YgMhI+S2z2zqFDybx58+rkNdyNgcoRqu/U3kDHW64XIvK4iGSJyEbnc7Xfb7NFJFNEtotIaFSmk+P1miAFPL0nOMvBBzPhZ69eJ/PxqJBkL+R4PNrzOugGq7bG3rIAgVvuq9ki0+uF75y1ZsM6No+Gzyd9urJKVKdu/T2GDGnfPAc8U7pXRnH3DVRqwq443UDj/h1akh8f27TXl9aEZWkH8QCrCmpvAyZfd/rieq9Xa0agUlNiqBfPP/98XHx8/FlDhw4d3rdv35Hjx48fUFJSIhkZGW1XrFhRJ3WRq5ZIKVUhIlNE5GMRWSQi/wZcBzTUDOqLSqnRzudTABEZDtwIjACuBF4XkYa3P7YsPRoF7bbIreVgWRl85ViwPd2/eY/EfGvuVriYt6qusfPn/u16sn9klJ4Pag4k94TL2ztWiQEae5+fUJ9pvb/HkAAqW6VgY2F4YgcHhWXpjhr0QKW2xr5q+W0bdpzQ+4Oim8fgpjpeGqQ1HOnFtUt31c1bPrBTa0ZGRjfPqYwmyKZNm2IeffTRrG3btm3Zs2dP+o4dO6LXrl0bvXz58o7Tp0+v05I3N+vsrnF2nwVeAl5Gh/i5AphXl5u6ZCJ6LV+pUmoPkAmcF5I7JfeE94bCfT3dWQ7atvZUnlGidfRNXVUXiOSecI/jZi7gvF2Vxs4f24btzsi2Om//TRXLgmGOZLOrvHbptuqclW+uakJgjyFpx+D6tfBARqXhYpPhpUF6gOKzSqyJKRMdIx1nbaJvIT248xjSVPF4tNcbCDwa8X8GvF742vGcdHHgZ14pWH8UFma1DN8MoWLr1q3R5557bhFARkZGW6UUI0eOLPF4PCoxMbFOvvTctEbPA58A6UqpL5xj39XlZrVwn4jcCqQCM5RS+UBvIMXvnJ+dYw2PZellB26XHiw4qEdyAA/1br6jWR+WBR7HAu6To3DLGTVHY558rV5nVdVAAbQK8MtCPeC9oGEc1IeNlwbBrhK9XOKD3JqNLJZ8rE3vAV54vNIwxYXXnIQOcG03eHk/ZJfA/LOb0KPj8UByd3hkH/SqRUthWYDSRjoA3zvLVcbHNv9Bn28pUqBn94NlWo055mz4fpgz6I1yNV+bdgyu/BHyyrXXtZvOaKC8h4A7+LZPOvkBV9srFEVUWDFE2IHWRI6kc9FbXBQwAnpmZmb07bff/guv1yvZ2dmR//jHP3Z26dLFnjlzZl4QRTgFN69abxGZDWwVkaec9XVBISJfiEh6NZ+J6Lm/AcBo4CDwgu+yapKqdsglIneLSKqIpObmBuUbNHj8zerHd6jdq3sYqXcdvDgQkuP1Eoy04zWfV52BAugR7tN79X4jzN3Uu/weD3ySAAtH6HncmvCfs/nDe9owZVx7Vw2dCJzfRT/Yi3Lg/vSGG903yDvwp37wTH+9rQ0bLaJ4K2BnkWOYEsDpsx+hUufWuw6Se8CCEYGfXduRar12pQXvQ31daTISOsCUHrohe2Vv03QuFCxFVFiZFEYXUdEgQ7fMzMw2Xbp0Kd+xY8eW3bt3Zzz77LP7n3jiiXqPpNxIdhuB40A/YCRwl4j8DKQBm5RSLwZKQCk13k1mROS/0D43QUtyfuF9OQPtwaW69P8O/B0gKSkpdDMitq3npHxm9YPbNZmheb3rwOOB90fqji6hfW13qrJ1uH97ZQyz21z6RWxAGuQZsKzAZvOWpSU6W8FXg4BSrdJ1qbK9sResOQqv7odXf4YcL8xPqFzbXVcapPwREYE7OqhcazbXCysdrylBuAdLOwb/sg4+SNR95eQGelzqXQdu/n/feZYFKYWVawvXHoOpgS8VgZeGw7pCWFMAF/8A346p//8fCtxIYAA2ihRyosfQrdhqAG83qampMYMGDTqpqjznnHOKXnnllXqvZwn4hiqlLql6zPGqMtL51AsR6amU8rmdnwSkO/vLgPki8legFzAIWFvf+9WL9w/Aqwf1SLa5+MAMBhEdyb02plzvO1l3/palh6ff+IJVuvB239zxqbEYApwBU92/h5YFc4ZDbhksOgSLD8G+IvjDmZDcXDTik6+D1evhm75AadDuwRI6wJJzIKMAbt0M3/aBl0c0k7KDnrcEeLMYUNqZdBCekyIi4OVhcEGK7vAuSoFVv2w+U9xVsRDG0r3BzFA3bdoUPWTIkGLQwVrffPPNrhdffHHhli1bIh9//PGehYWFEZ9//rmLKMJV81kHlFI/KaU+U0o9X5frq/CcEzZoE3AZ8KBzjwx0GKEtwOfA75VSjTut/54zihsZDZ+Nbr5PZ304ubD6vkrT62nbtSUr6GgSzabVqiNegDHAfidAaXBu5iwL3k+AGxw/1D8Uwi2b4YLv4M/bmqDxSlWWfAyvl0LGHhjqqT1AaTWIwOiOsL9Ea0Rf+wmmNaBKN+RYFpRVwJdLgQq4qFPQbcE5neHtUdAG3eGd813LUGk2BOnp6dGLFi2KGzZs2PCRI0cOKykpsf76178eGD58eNnixYv3BU6hehq9tVZK1Sj8K6WepCkFiF0+EiZmaO8GLb1BrxWp3JaVwSJnfmRUVMuX6gBWDwLOBs4LHKC0BiIiYMHZcN3P8Mo++OEYpDifhQchuRfMGtg01VtMvBIiS6EsDtocrvOg74/94esj8M/D8PrP2mjj/dFNtMz+eL0wPR/4HjqfA3MuDzoJEUg+Q08G/HYzbDoBZ38Hb43SHWFN9mGtgWXLljVc3EI/WnOLHTzGu4HGZ3o98Qro8T0csaGzBalJLX8Q4PU6xjlvQO8N9ercLQtu7gvfjYUn+sOEznr0mVEMj+6Cq39sglKebcN5G6CsC0SmQUrdfT1ERMCnSVoDqtBLPJtkmavywE440h94FF68sc6dvYi2xnx7lP7fN5+A81Ng/s9NcC1mC0BUC6tVEckF3Iq6cWi3Z82BIUqpABNqmiDqoE7lj0M69yXiFzbdKyJo5wHYR/buPAryg00rCEJRfgiyDs6ke584OnaDfPZxeHceqmHLHNUumvi+fWkbo62Eio8VkLXDF1nUVR2EsvxxxHY+k2794QSQwz4qGqYOup3Zh45xWrF74ugRDu6qbnTfJN4B/QzEdlOcUPs5uCcPu/7lj2oXTbd+/YiMikEpRdb2rZScqDoP5rr8dSUtLW1vQkJCk2wTDx06FDF9+vTe3377bcdbbrkl7+mnn67R3VFaWlpcQkJCP/9jja7GbGiUUrXYjZ+KiKQqpZJCmZ+GQkRS3Z7rtg5ae/l96ba0Omjt5QfzDrREevToUTF//vz9db2+heucDAaDwWAwnZ3BYDAYWgGtvbP7e2NnIAhCkdfWXv5QphsKzDPQPNIMFc0pr02OFmegYjAYDIa60ZQNVIKhOgOV1i7ZGQwGg6EV0KI7OxF5S0RyRCTd71gXEVkhIjudbWfnuIjIy06w2E0icnYY8xklImtFJE1EMkTk353jvxCRH5y8LhKRSOd4W+d7pvN7vxrSbRbld+7f4HXQ2svvnNcs6qC1l9+5f0jqwKBp0Z0dMBcd+NWfWcCXSqlBwJfOd4Cr0P43BwF3o6MxhItS4HKlVAI6+sOVIjIGHUPwRSev+cCdzvl3AvlKqYHAi8551TGX5lF+CE0dzKV1lx+aTx209vJD6OrAQAvv7JRSq4AjVQ5PBJyAZLwDXO93fJ7SpACdRCQsvq+ce/ri6rRxPgq4HPhHDXn1leEfwDiR0x0MNZfyQ2jqoLWX30m3WdRBay8/hK4OWhrz5s3r5L91S4vu7Gqguy/KgrN13PHSG/APaRG6YLHVICIRIrIRHQV+BbALOKqU8rmH9c/Pybw6vxcAXV3eqkmWH8JWB629/NBE66C1lx/CWgdNmueffz4uPj7+rKFDhw73fdauXRudkZHRdsWKFR1922DSbHEeVOqB62CxocCJ6DBaRDoBS4FhteQnFHlt1PJDo9dBay9/qNJ0TWsvPzSJOmgSbNq0Kebhhx8+8OCDD55iGfrMM8/ET58+PWf58uUdp0+fnhNMmq1Rssv2qSacra/CXAeLDSVKqaPA1+gYMp1ExDcg8c/Pybw6v8dyuqqmJpp0+SHkddDayw9NvA5ae/khLHXQpNm6dWt0UlJSUdXjHo9HJSYmlvi2waTZGiW7ZcBtwDPO9iO/4/eJyELgl0CBqgwqG1JEJB4oV0odFZFoYDx6svkr4F+BhdXk9TZgjfP7SuV+wWSTKz+EtQ5ae/n9r20yddDayw9hr4PA3PFgH9K3xQQ8TykoKraIibYDxiYaObSIt14MGAE9MzMz+s477+xnOVFU7rrrrpyZM2fmzZw5Mw/Atw2GFt3ZicgC4FIgTkR+Bv4f+gFfLCJ3AvuByc7pnwJXA5lAEfDbMGa1J/COiESgpe3FSqmPRWQLsFBE/gxsAN50zn8TeFdEMtEjuRurS7QZlR9CUAetvfzQrOqgtZcfQlQHIaeo2CJzbzQD+xXTLqbeIXgzMzPbdOnSpXzHjh1bGiJ7PowHFYPBYDAAdfSgYtuQsj6aMWcXN0Q8y4ULF8a++eabcStWrNhV1zRaRYgfg8FgMIQRy4KxSVVj79WZTZs2RQ8ePDio+Tg3mM7OYDAYDE2G9PT06B9++KHDl19+GQsgIqSkpGyLjY2tl4rUdHYGg8FgaDIsW7asuij19aY1Lj0wGAwGQyvDdHYGg8FgaPGYzs7QqIhIhYhsFJF0EflARAKv6zn1+gf8rxGRTx3vEzWd/7iIzKxPnl3ma66I/Kuz/98iMryB0v1cRMLqwspgaAmYzs7Q2BQrpUYrpUYCZcA9bi901iM9AJzs7JRSVzveJxocPy8WQaGUukspVe81Q85C4y5Kqaz6pmUwtDZMZ2doSnwLDAQQkf8RkXWi43rd7TtBRI6LyH+IyA/AI0Av4CsR+cr5fa+IxDn7t4qOS5YmIu9WvZmIDHAkpXUi8q2IDK3mnMdF5O8i8k9gnoj0c85d73zGOueJiLwqIltE5BMqnQsjIl+LSJIv/37H/1VE5jr7kx3pNk1EVtVQP5eiXUhVzePXIvKs6FhoO0TkIuf47U49LheRPSJyn4hMF5ENIpIiIl1q/CcMhhaGscY0NAkcqekq4HPn0B1KqSOONPOjiCxRSh0G2gHpSqnHnOvuAC5TSuVVSW8EujO8QCmVV0PD/nfgHqXUThH5JfA6OpxKVc4BLlRKFTsq0wlKqRIRGQQsAJKAScAQYBTQHdgCvBVEFTwGXKGUyqpFDXsV8D81/OZRSp0nIlejvYSMd46PBBKBKLRnkD8ppRJF5EXgVuClIPJoMDRbTGdnaGyiRYc0AS3Z+VwhTRORSc5+H3RAzcNABbDERbqXA//wdYJKqVMc5IpIe2As8IFU+vNrW0Nay5RSvkWzbYBXRWS0k5fBzvGLgQWO1/oDIrLSRR79+Q6YKyKLgQ9rOOcCoKb5Rt8164B+fse/UkodA46JSAGw3Dm+GTgryDwaDM0W09kZGptipdRo/wMicilaMjlfKVUkIl+jJROAEqdDCYRQe7gTCx0nbHQt5/g44bf/IJANJDhp+Ht6cON7z/+cqJMHlbrHkS6vATaKyGhHkgVARPoDPymlympIt9TZVnDqe13qt2/7fbcx77+hmTJv3rxOt956a1Bz82bOztAUiQXynY5uKDrMSU0cAzpUc/xL4AYR6QpQVY2plCoE9ojIZOd3EZEEl3k7qJSygalAhHN8FXCj6OCbPYHLarg+W0SGiYiFVn3i3H+AUuoHRz2bx6mhZuBUFa/B0Cq49dZb+/bq1WuU/7G6BG4F09kZmiafAx4R2QQ8AaTUcu7fgc98Bio+lFIZwJPANyKSBvy1mmtvBu50fs8AJrrI2+vAbSKSglZh+qS+pcBOtHrwDeCbGq6fBXwMrAT8w8c8LyKbRSQd3XGmVbnuSkxnZ2hFbNu2LXLNmjUdysvLJT8//2RfVZfArWCiHhgMTR4RaQt8p5RKauy8GFo2dYp6ECImTZrU76qrrip4++234//85z9njRs37gTAX/7yl7hA8exM1AODoRmilCpFW3waDGEjjTv6FJIe0MmDQmFTZFnE2ELtwVs7MrIogbcCBm9NTU2N2r59e/SSJUv2rl69usOmTZuifZ1dXQK3glFjGgwGg6Ee2BRZJ8iMtilqsP5k1qxZvR977LEDlmUxbNiw4vT09Oj6pmkkO4PBYDCchhsJDEBhk09KdGfGFEsDyE8rV65st3r16titW7fGzJgxg9LSUmvIkCFF9U3XSHYGg8FgqDOCRRfGNkhHB/DII4/0XrRo0c6srKzNWVlZmzds2JCxc+fOaIB3332304033njmuHHjBnz44YdBWWSazs5gMBgMTYKlS5d2LCsrk4kTJx7zHevTp4+3uLjYys7Ojpg6derRhQsX7lu4cOHehQsXdg4mbaPGNBgMBkOTYNKkSYWTJk0qrHr82LFjG/2/P/zwwz2nTZuWG0zaRrIzGAwGQ7PAtm3uvffe3tdcc03BhRdeGNQ8npHsDAaDwdAseOqpp7qtWrWqY0FBQcSOHTui/vjHP7qW7syicoPBYDAATWtReX2oblG5UWMaDAaDocVjOjuDwWAwtHhMZ2cwGAyGFo/p7AwGg8Hgw7Ztu3YHl00cJ/921eOmszMYDAaDj/Tc3NzY5trh2bYtubm5sUB61d/M0gODwWAwAOD1eu86dOjQfx86dGgkzVMYsoF0r9d7V9UfzNIDg8FgMLR4mmPPbTAYDAZDUJjOzmAwGAwtHtPZGQwGg6HFYzo7g8FgMLR4TGdnMBgMhhaP6ewMBoPB0OIxnZ3BYDAYWjz/Hx6x3+1+Kh/3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x216 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(6,3))\n",
    "axesR = {}\n",
    "axesI = {}\n",
    "for i, lMax in enumerate(lMaxes):\n",
    "    axesR[lMax] = fig.add_subplot(2,len(lMaxes),i+1)\n",
    "    axesR[lMax].set_xlim([50,300])\n",
    "    axesR[lMax].set_ylim([1.25, ef[1]/eh])\n",
    "    axesI[lMax] = fig.add_subplot(2,len(lMaxes),len(lMaxes)+i+1)\n",
    "    axesI[lMax].set_xlim([50,300])\n",
    "    axesI[lMax].set_ylim([-60, 30])\n",
    "    axesR[lMax].set_title('$l_\\max = %d $' % lMax)    \n",
    "    axesR[lMax].tick_params(labelbottom=False) \n",
    "    if i == len(lMaxes)//2:\n",
    "        axesI[lMax].set_xlabel(\"Particle radius / nm\")\n",
    "    if i == 0:\n",
    "        axesR[lMax].set_ylabel('$\\hbar \\Re \\omega / \\mathrm{eV}$')\n",
    "        axesI[lMax].set_ylabel('$\\hbar \\Im \\omega / \\mathrm{meV}$')\n",
    "    else:\n",
    "        axesR[lMax].tick_params(labelleft=False) \n",
    "        axesI[lMax].tick_params(labelleft=False) \n",
    "    axesR[lMax].set_xticks([100,300])\n",
    "    axesI[lMax].set_xticks([100,300])\n",
    "\n",
    "res_thr = 0.00015\n",
    "\n",
    "ir_labeled=set()\n",
    "if True:\n",
    "  for (lMax, radius), (eigvals, residuals, b, ef, irclass) in plotdata.items():\n",
    "    for i, (e, res, iri) in enumerate(zip(eigvals, residuals, irclass)):\n",
    "        #if i == 0:\n",
    "        if res < res_thr:# and e.real < 2.14e15:\n",
    "            if iri in ir_labeled: \n",
    "                label=None\n",
    "            else:\n",
    "                ir_labeled.add(iri)\n",
    "                label=D4h_irlabels[irlabels[iri]]\n",
    "            axesR[lMax].plot(radius*1e9, e.real/eh, \n",
    "                             marker='.',\n",
    "                             #marker=markerfun(b),\n",
    "                    ms=1, #c=colorfun(b)\n",
    "                    c=matplotlib.cm.hsv(iri/9),\n",
    "                    #c = colorfun(iri),\n",
    "                    label=label,\n",
    "                   )\n",
    "            axesI[lMax].plot(radius*1e9, e.imag/eh*1000, \n",
    "                    #marker='x', \n",
    "                    #c=colorfun(b), \n",
    "                    c=matplotlib.cm.hsv(iri/9),#colorfun(iri),\n",
    "                    marker='.', #markerfun(b),\n",
    "                             ms=1,\n",
    "                    #label=label\n",
    "                            )\n",
    "fig.legend(title=\"Irrep\", loc=\"upper right\")\n",
    "plt.subplots_adjust(bottom=0.3, right=0.82,left=0.10)\n",
    "#fig.suptitle('$l_\\mathrm{max}=%d$, residual threshold = %g' % (lMax, res_thr) )\n",
    "fig.savefig(plotfilename)\n",
    "fig.savefig(plotfilename.replace('pdf', 'png'))\n",
    "print(plotfilename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.        , 1.40635433, 1.98888536, 2.81270865, 3.14470387,\n",
       "       3.97777072, 4.21906298, 4.44728287])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ef / eh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1207"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(plotdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5, 3, 5, 0, 7, 8, 8, 9, 9])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "irclass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}