qpms/lepaper/figs/ewald_branchcuts.ipynb

152 lines
43 KiB
Plaintext

{
"cells": [
{
"cell_type": "code",
"execution_count": 123,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt, gridspec\n",
"from qpms import EpsMuGenerator, lorentz_drude, unitcell, eV, hbar, BaseSpec, CTMatrix, Particle, pi\n",
"from matplotlib.colors import BoundaryNorm, Normalize\n",
"eh = eV/hbar\n",
"ccycle=plt.rcParams['axes.prop_cycle'].by_key()['color']"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"# I only need to build a valid unitcell to use its empty_lattice_modes method\n",
"px = 580e-9\n",
"py = px\n",
"\n",
"a1 = np.array([px, 0])\n",
"a2 = np.array([0, py])\n",
"\n",
"bspec = BaseSpec(lMax=3)\n",
"tm = CTMatrix(bspec, None)\n",
"emg = EpsMuGenerator(lorentz_drude['Au'])\n",
"pp = Particle((0,0), tm, r=50e-9, material=emg)\n",
"u = unitcell(a1, a2, [pp])\n"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [],
"source": [
"fn = 10\n",
"kxs = np.linspace(0, pi/px)\n",
"emptyfreqs = np.empty(kxs.shape+(fn,))\n",
"for i, kx in enumerate(kxs):\n",
" emptyfreqs[i,:] = u.empty_lattice_modes((kx,0),3*eh)[:fn]\n",
"emptykappas = u.omega2k(emptyfreqs)\n",
"mykx = 0.2*pi/px"
]
},
{
"cell_type": "code",
"execution_count": 132,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEOCAYAAACetPCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd8VFX6h58zM+m9FxISSgiEntAJSrEgirJrpQkqAq7guurqT11ABHVV1oqgFEWKBSsqIiiK1AQIECAJkAAJCaSQRnqZzPn9cZMQMJA2k0m5z+dzDJl777nvxHPmO+e+531fIaVERUVFRUWlvmjMbYCKioqKSutCFQ4VFRUVlQahCoeKioqKSoNQhUNFRUVFpUGowqGioqKi0iBU4VBRUVFRaRDNJhxCCGshxH4hRLQQIkYIsbCWc6yEEF8KIRKEEJFCiMDmsk9FRUWdpyr1ozlXHKXAaCllX6AfMFYIMeSqcx4BcqSUXYG3gdeb0T4VFRV1nqrUg2YTDqlQUPmrRWW7OvrwLuDTyn9/DYwRQohmMlFFpd2jzlOV+qBrzpsJIbRAFNAV+EBKGXnVKR2AZAAppV4IcQlwAzKv6mcmMBPAzs4urHv37qY2vRmQUF4C5cWgLwZ9GVSUgUEPsgKkAa4b5S9AVDUNCC1oqpquslmAtvJn5b9LSw3kZBfi4maPlbVFk96BPrcUQ2H5X17X2Fmgc7ZqUt8qTScqKipTSulR13ktdZ5Wja9zGADoiAYpDSRfSqOoooSAgACcnJyadI+2jpSSc+fOkZmZCVCv8VAbzSocUsoKoJ8Qwhn4TgjRS0p5vBH9rABWAAwYMEAePHjQyJa2UAwGyDsPF09AVgLkJEHeBSjMgKJsKL0EZUWgL6kUHANgAMoq21+Rlf+RZKHRVgpMlbhoLTBoLCjDghIsKZJWlEod5cKCcrSUY4FeailHg15q8c96GB12f72JAL/XRpju76JSL4QQSfU5r6XO05Tnd4GEORQCsLRyrJVWlNHlzTFN6rs9kJiYyJQpU8jMzOS+++5j48aN9RoPtdGswlGFlDJXCPEHMBaoOSDPA/5AihBCBzgBWWYwsWWi0YCzv9KCbq77fIMBCtIhNxnyUqAgAwozKCvIJis7h9z8fIqLi8nPL0JIPZaATl+ODj06KtChx4IKrEQJluixpRxHKrCoPK4VV66AUni8djvUdGitkpY0Tw0GA9Igqe2JmJXW0pS3bhN89tlnPPbYYwBs2LCBSZMmsXHjxkb312zCIYTwAMorB6MNcDN/dar9AEwD9gH3AL9LNQtjgzEYDJxIz2dPfCZHknM5l11KZoETl4ptKSn3x1DbX1RKtAaBpU6LlYUGW0sd9lZaHK0tcLa1wMXOEg97KxytLbCx1GJtocFaK7DVGbDV6LHSgsfa82hr61t9+t1qaGnz1GAwsGTJEl5++WV+m/oJfk7eTOOqx57q+LomeXl5PP7446xfv55hw4axYcMGAgMDm9xvc644fIBPK5+faoCNUsqfhBAvAwellD8Aq4F1QogEIBt4oBnta3UYDAaiUy6xJyGT6ORcTl8sJD2/hKLSir98ybfQCmwtdXg6WOHpaI2/iw1dPO3p7u1IiJc9s8cuYcjI7vzfovsabU/2YANFEWl/ed12sHej+1RpdlrMPF2xYgXPPfccubm5aLVaosvO0AEvBl71saWOr9rZu3cvkydPJjk5mYULF/LCCy+g0xnnI7/ZhENKeRToX8vr82v8uwS4t7lsak0YDAYOJ19iy/FUDpzN5kxmIfkl+ivOEYC9tY6unvZ09bSnn78zw7u6E+LjgEZz/Q10ocO6cmD3KQwGQ53nXgvXCUEAFEWmKY+nhDKpq15Xafm0hHm6ceNG5s6dS0ZGBkII7r77bj7++GMcHR3J/j6e6IjzAAQJrTq+akGv1/Pqq6/y8ssv4+/vz65duxg6dKhR72EWH4dK3STnFPFj9AX2JmRyIq2A7MLSKx4x2VpqCfZyoLuPA/38nRnR1Z2uXg6Nvt/A8G7s3R5L7JFz9AoNbHQ/rhOC1Ims0ii2bNnCrFmzSE5ORgjBLbfcwqeffoq39+UVheuEID5MVzZvfTnLuB+GbYEqB/iePXuYMmUKH3zwAY6Ojka/jyocLYT0vBI2RCSxPS6DhIsFlOoN1ccstAI/Fxv6+DkzKtiDW3v6YG9t3P91Q0d256M3LPn56wNNEg4VlYayY8cOZsyYwenTpwEYPnw469ato1OnTma2rHVR0wG+fv16Jk+ebLJ7tX7haKW+84ISPV9FJbPlWCoxqXkUllZUH3O2sWBgoAvhXT0Y38+HDs62JrfH2c2esGFdidhxAilr372iomJMIiMjmT59OidOnAAgLCyMTz/9lJ49e5rZstbF1Q7w9evXm1x0W79wpB6BNzqDTz/ocSf0vges7M1tVa3sP5vFql1nOZCYTU7R5UA5W0stYQHOjO3lzX1h/jjZmmd74cAR3dizPZaYw0nqqkPFZBw5coQHH3yQY8eOAdCrVy8++eQTBgwYYGbLWh/79u1j8uTJJCUl8dJLL/Hiiy8azQF+PVq/cOisoTQfTm9X2k//BGsn8AyBrjdD3wfAqYNZTDMYDGyNTWfdviQOncuhpFx5/GShFXT3dmB0d08mDu6Iv4vpVxT1YfCNSmTv/l2nVOFQMToxMTFMmzaNqKgoAIKDg1m9ejXDhw83s2Wtj9oc4MOGDWu2+4vWHiZRHZGafgKi18OZHZCZoKTtqEJnDa6doOMw6Pk3CBiuBNOZAL3ewJdRyXx5IJnYC3noKz3adlZaBgW68kh4J8KDGhXl3yw8ft8HVFQY+PCbueY2RcXICCGipJRm+Vpvb28vCwuViO8uXbqwYsUKRo8e3ai+opKyAQgLcDWafa2Jmg7wyZMn88EHHzQq1UpTxkPrX3FU4dUdbll8+fe8VIj+HBJ+hYy4y+3gakCAvQd49oQuo6Dn35Vo7EZiMBj4ITqVj3ae5mRafvXuJxdbC27o5sGsGzoT4ts6cugMuiGYz1fsICsjDzdP4+/GUGmfFBYW4uPjw+rVq7ntttua1Fd7FQyAzz//nNmzZyOlNLkD/Hq0nRVHXejL4MRPcGIznI+CSylgqJGQT2sJDr7g3UcRk5A7wc79ul3Gpl7if1tPsTshs3oXlJeDFbf09Gb2yM7N4tQ2NnHRyfxr6kfMnXcnt987yNzmqBgRc644hFDy0/Tv35+PP/6Yfv36Nbqv9rjiyMvLY86cOaxbt46hQ4eyYcOGJjvAmzIe2o9w1EZOIhz7Bs7+qaxGijIrEwNWorUCRx9lZdJpBASPI9/Gl3e3x/P9ofNkFiqJA20ttYzp4ckztwQT4FZLkr9WhJSS+254lV6hgSx41zzfZlRMgzmFIzg4WBoMBhISEgAIDQ3lk08+oU+fPg3u6/6P9gHtJ44jIiKCSZMmkZSUxLx58/jPf/5jFAe4+qiqsbgEwg1PK62K84fhxI+QtFfxleSeg5xEtsee541NcEr6IxFoMNDfsYC5g10ZfeMo0LWNRGtCCLr19CMxPk3dlqtiNBwcHDh48CDbt29n9uzZHDp0iL59+9KvXz9Wr15NaGiouU1scZjbAX492rdw1EaH/kpDcXS//0c8n+45Q265AZD4aXKYqN3Oo2ITlmUG2IXSdNZg7w3uXcG3PwSOgI5DW6Wg9B4QSNTeeFKTs/Ht6GZuc1TaEGPGjCE+Pr5aQI4cOUJYWBi9e/dm9erVDBw40NwmtgiM5QA3Fapw1EJ6Xgkv/RDDr7Hp6A0SjYARQe4suCOkMq3HVKX+Rfw2SNwNqUfhUjJcOge5iZDwG+x8U+lMawW2buDcETx7QIcw5bGXS6AZ3+H1GXxDMGve+5WIP0/w96nqVkkV41MlIL///juzZ8/m2LFjDBo0iJ49e7Jy5Uqj51ZqTbQUB/j1UIWjBnsTMlm8OZbY1HxA8V1MHuDPc2ODsbW86k9l66rEiPS9KjFoThKc/gNS9kNGrOKEL8yA/AuQHAFRn1SeKMDCRhEVJz9w7QLevcA3DHz6gIW16d/wNejUzRt3L0cO7olXhUPFpIwePZpTp06xe/duHn30UWJiYhg2bBjBwcF8+OGHjBw50twmNhumcICbClU4gG0xaSz4IYbUSyUA+LnY8MSYIO4b0Igtui4BMGC60mpSWqCsTlL2Q9pxyE1SiizlpyqrlXP7rjxfowULO7BxBQdvZbuwaxfwCAav3uDWxWSxKABBIR04eTylSdlyVVTqS3h4OHFxcURGRvLoo49y7NgxRo0aRefOnXn//fcZN25c9bnzx4eY0VLTEBERweTJk0lMTGzWCPDG0nItawa2Hk/jpR8vC0b/js68MqGXaWIurOwheKzSrqasSNkinHJAKQubk6gISnGOIiq5icpq5WqEVvGtWNmDtTPYuYG9Fzh2UFYxzh3BpQu4BjZ4BdN3UGf2/RFHYnw6nYN96n1d9vfxalp1lUYzePBgjh49ypEjR3jkkUc4dOgQt99+O35+fixZsoSbrUJxqhxfKW1gfOn1el577TUWLlyIv78/O3fubBWR9O1SOLYcT2XhDzGk5ZUCMCDAhTfv6UMnDzPluLK0Vfwena5Rl1tfBmnHlEdfmSchOwnyz0NhFpTkQnGuUhY2s46t1RqtUktcZ60IiYWdcm9Le7ByqGyOYOXAcEtI6RxN0d5PoaKnIk5W9mDhAJY2lX1U/tToQKNRRKNmISdJ9e+teXKrND/9+vUjKiqKuLg4Hn74YSIiIvh84SoGjPs/juqUFfBAqWvV46umA3zSpEksW7asRTnAr0e7Eo6rBWNgoAtL7u3b8mMvdJbgF6a061FeApmnIPuM8igs7wLkp0FRlrJ6KS2AskKoKIGiHDBcBFlRa1cewJx+QNJOqEdJe4MUlJSuBjz/cqwoMq1VTmwV89OjRw/27dvH2bNnKXo3DiudJZ+ipC6pqgTYGsdXa3CAX492IRxHU3KZ89khzmUr+atajWA0FAtrxbHu89egqvyScvadziI2NY8zFwtJySniYn4pucXlFJXpsTIUY0cptqIEO0qwpgxrUYY1ZdhRgq0oxUGUYKMpx06UYSEqsKIcHQYs0KPVSMLwqL38c+uOMVVpAXTq1IkU65Raj0mDRK/Xt2ifQBV5eXnMnTuXtWvXtngH+PVo+X/pJpBdUMaczw+x93QWAP38nXn3gX5tTzBqkFlQwq5TWRxIzCIuLZ+UnGJyCsuqky3WRCsEtpZafJxs8HBwxs/FBj9nW7ydrPFxsubAz0f44/N9/Lh7Pja2VnXeO+X5XbWLhBpDqGIMBLWOr/N56XS1s2Pq1Km888472Nu3zLIKNVOgL1iwwGgR4OagdVpdBwaDgQU/xLIhMgmDBG9HK965vz9DurStYLary8vmFJZRcVUKGUutBnd7Szq62dHd24EePo709HWku7cjlrrr75aySPRhl76CUzHn6Tuwc5322A72vtLHUeN1FZWmUtv4kkhiSUaj0bB69WrWrFnD+PHjWb58+RUlZ81JRUUFr776arUDvCVFgDeWNiccn+8/x6KfYikqq8DaQsNzY7vz0PDWtxS8moISPb8cT+WPkxkcS7lEal4J5RWXRcJCK/BysqaTmy29OjgyuLM7Qzq7/jX+pAF06aHspjp7Mq1ewlH1nFndVaViCqrHUYRSYhYBdoN9eHjCU0z/4EkWLVrEW2+9xffff8+mTZsYMWIEK1asIDg42Gw2JyUlMWXKFHbv3t0iI8AbS5tJcnj6YgEPfbyfcznFaATcE+bHqxN6o6vjW3VLpUxv4Icj5/nh6AWOJl8it/hyJl8hwM3Oiu7e9oR39eD2vj4mKQYlpWTKzW8Q2NWbVz6cZvT+VZoXcyY5bFIy0qs4fbEAgC7X2AW5YsUKFixYQFqasjrp3bs37733XrMHE9Z0gC9btowpU6Y06/3rot0nOVzwQwxr9yYigbAAZ5ZPDsPT0XyR143BYDCw41Qm30QlcyAxh4z80upjVjoNPX0cGdrVjdt7+9DXz6lZgvKEEHTv7c+JY8lqwkOVFsO1BKOKmTNnMnPmTH7++WeefPLJ6mBCPz8/Fi9ezLRppv0SdHUE+Pr16+ncue4Ve2ui1QtHXGoemXsTsbXU8vZ9/bi1V8t4rlkfCkr0fLovke8Pn+f0xYLqAlBajaCTux03dPNg8uCOdPNyMJuN3Xp1YM/2WLIzC3DzMJ8dKipV/BabDsBNIV7XPW/cuHGMGzeOQ4cO8dhjj7F//36mT5/OP//5T5544gnmz59vdOd0zQjw+fPnM2/evFbrAL8erf4d6Q2Sm3p4smxyWJ3O3pZARl4JK3ad4ZfjaaTkXC5v6+VoxdDO7jww0L9FOfG7VEaNnzqewtBRPcxsjYoKrNx1BqhbOKoIDQ0lMjKSlJQUHnvsMbZs2cKiRYt4/fXXuf/++3nvvfdwdnZukk01HeB+fn6tJgK8sbR64ejkbseqaS07FXNyThHvbY/njxMZZBYoxZ8EEOhmyx19fJkxohPOti0z/XqnbsoKLiH2giocKq0aPz8/fvzxR4qKinjmmWdYs2YN69atY/369dx4440sXbqUnj17NrjfxMREpk6d2uYc4Nej1QuHvVXLfAtlegMrdp7ms/3JXMhVVhZaIQjxceCeMH8mDfLHugk7npoLN09HHJ1tyc7MN7cpKipGwdbWlmXLlrF06VLeeecdXn/9dXbs2EGvXr3o1q0bb7zxBnfddVe9+vriiy+YPXs2BoOhVUaAN5Zme7YjhPAXQvwhhIgVQsQIIf5ZyzkjhRCXhBBHKtv85rLPWGyPS2fCB3voPm8LS7adIjW3mGAvB964pw/xr4zl53/ewMPhnVqFaFTh4+fKqZjz5jZDxcS0lzlahUaj4amnniI9PZ3NmzcTEhLCqVOnmDBhAm5ubrz88svo9fpar83Ly2PatGlMnDiRkJAQoqOj241oQPOuOPTA01LKQ0IIByBKCPGrlDL2qvN2SSnvaEa7mkx6XgmLN8eyPS6DojIl95O7vSV/69+BJ8YE4WBtYWYLm0bnYG+2/3SE0pJyrFr5e1G5Lm12jtZFlSM9Pj6eOXPm8Ntvv7FgwQIWL17MhAkTeO+996oDCms6wFt7BHhjabZ3K6VMBVIr/50vhIgDOgBXD8pWQ8TpLBZtjiXmQh6gbJu9uYcnT98STHcfRzNbZzy6hviy5ZuDXDiXVe3zUGl7tJY5+vb9/UzWd1BQEFu3bqWgoIDnn3+eTz75hK+++oqvvvqKrl27Ymdnx/Hjx1tVCnRTYBaZFEIEAv2ByFoODxVCRAMXgGeklDHNaFqdGAwGVu8+y0c7z5JZoMRadHC2Yc7orkwc1NHM1pmG4F5+ACSdzlCFo53Qkueor7ONye9hb2/P+++/z7vvvsuqVatYuHAhCQkJ1cemTZvWrsvbNrtwCCHsgW+AJ6WUeVcdPgQESCkLhBDjgO+Bv+SrEELMBGYCdOzYPB/WuUVlvPxjLJuPpVKqNyBQ6njMHx9CH7+mbeVr6Ti7KgFXKYmZZrZEpTkwxhyt7Mck8/TH6AsAjO/ra7Q+r4VGo8HR0ZHCQiWVu62tLYWFhSxcuJDXXnuNO++8k7fffhs/Pz+T29KSaNbAByGEBcqA3CCl/Pbq41LKPCllQeW/fwYshBDutZy3Qko5QEo5wMPDw6Q25xaVMWvtQUIX/cq3hxUH8b1hfhyefzNfPzaszYsGgLuXI86udlxMu2RuU1RMjLHmaOVxk8zT9RFJrI+oR5GYJnK1AzwhIYHCwkKys7OZNWsWFhYWfP311/j7+9OrVy82bdpkcptaCs25q0oAq4E4KeVb1zjHu/I8hBCDKu3Lai4ba5JbVMbMSsHYGpuOvZWOF8d1J+7lW3nz3r4tNu7CVHj4OJMQd8HcZqiYkNY2R01JREQE/fv3Z/369SxYsICdO3fSpUsXAJydnfnwww8pKChgzZo1dO7cmZiYGCZMmICzszP/+te/KCgoMPM7MC3NueIYDkwFRtfYyjdOCDFbCDG78px7gOOVz0/fAx6QzZyFsaZgbKsUjPl3hHD0pVt59IYuzZIjqiXSoaMbudkFtPakmCrXpVXMUVNSUVHBokWLCA8Pp6Kigp07d/LSSy9dc9fUtGnTOH36NCdOnGDs2LEUFhbyzjvv4OTkRHh4OHv27Gnmd9A8NOeuqt3UUdJHSrkUWNqgjg21lz5tKEVlep7eGM3WmDQMEhytdTx5UzceDm/9KdmNQWBXT3ZsOUpOVgGu7mrOqraIyeZoK6FmCvSG1gAPDg5my5Yt6PV6XnnlFZYvX86ePXsIDw/Hw8ODWbNm8eKLL2Jt3bqSr16L1v/1Oe0YfHwb7H4b0o5DA7/8GAwG3vntFH0XbmPL8TTsrXQsGK+sMFTRuIxrZYLDrIyrfaUqKq2fL774gr59+xIdHc26devYsGFDo9KG6HS66pTue/fuZcSIEWRlZbF48WLs7e0ZOXIk+/btM8E7aF5av3A4eEFZPvz2Enw4HN4KgU1zIHYTFOde99JtMWmELf6Nd36LB+CJ0V05Mv/mNlH4ydh066nsGjl35qKZLVFp7yyfEsbyKWFG6SsvL48HH3yw2gF+5MgRo9XNGDp0KDt37iQ/P59nn30WZ2dn/vzzT4YNG4anpyfPPfccRUVFRrlXc9MGhMMHZu+Gp+LgzvfBL0wRjY0PwhudYfUtsON1SN4PFUr6gNMXC7jl7T+ZuS6KnKJybgnx4vC8W3jqluB268OoCydXpU57siocKmbG1c4SV7umb06pcoBv2LCh2gFuiroZtra2vP7662RmZrJr1y7Cw8PJysrijTfewMHBgUGDBrFlyxaj39eUtJ04eUdfCH1QaRXlkHIAErbD6e2w4zXY8Sp6Cxf+qXmWzZc6AoJgLweWTQmtszCMCri42ePm6UCm+qhKxcx8dTAZgHsH+Dfq+oqKCl577TVeeumlZk+BHh4ezq5duygpKeHVV19l5cqVHDhwgHHjxuHg4MA999zDq6++2mLqpV+Ltvn1WmsBAcNgzDyYuQOePcOOoWvoV/Qumy8F4EwBqyzeZGvFDLrsfBKiPoWs0w32j7Q3vHxdiFeTHaqYma+jUvg6KqVR1yYlJTFy5EjmzZvH/fffT3R0tFnShlhbW/Pyyy+TmprK4cOHue222ygtLeWTTz7Bx8eHoKAg3n333WsmWTQ3bVM4alBUpmfq5/FM/8OSwgod9w/059AzA7lpwnQIHAFn/oQfn4D3QxX/yNePwMGPIeOEKiRX4evvSnFxmbnNUFFpFMZygBubfv368fPPP1NcXMyHH35IcHAwCQkJPPnkk9jY2DBy5Eh27NhhbjOvoO08qqqFTUfO89zXRynRG/B2tOLjhwYS4lM5UNy7QNh0RRwy4yFxJyTugcRdcPxr5RwbV+g4FDoOBv8h4NsPdFZmez/mxsHZlsz0PPLzinFwNH2+IBUVY5Cfn8+cOXNYu3YtQ4YMYcOGDS2yBrhGo2HWrFnMmjWLzMxM5s2bx8aNG/nzzz8ZNWoUDg4O3HnnnSxatIhOncy7gadNCkdBiZ5pn0QSlZSLRsCsGzrx/LiQ2k8WAjy6KW3gDEVIss9A0h44FwFJe+HkZuVcrSX49AP/QdAhDPwGgpOf0kc7oFOQF4YKA3k5hapwqLQKIiMjmTRpUqurAe7u7s7y5ctZvnw5kZGRzJ8/nx07drBhwwY2bNiAr68vU6ZM4cUXX8TRsfkzcbf8v2ADiTidxUNrDlBcXkGgmy3rHh6Mv5tt/TsQAty6KC30QeW1ggxIjlR2ZiXvh/0roaIyBsrOUxGRDqHgG6qsSuxqTd3T6vH1V2qhn4o5T4eAtvkeVdoGVzvA//zzT8LDw81tVqMYPHgwW7duxWAw8Nlnn/G///2P6Oho3njjDd588026devGzJkzmTNnDpaWzZMKSbT2bAEDBgyQBw8eBODlH2P4eE8iUMcqo6noyyD9OKQchAuH4HwUZJ66fNzRTxEQ7z7g0we8e4Njh1a/Mkm/kMO0sf/jyZcmMPbvA8xtjkoDEUJESSnN8j+u5jxtKsWVxdJsLLW1Hq8ZAT5x4kSWL1/eInwZxqSkpIQlS5bw8ccfc/bsWUB51NW3b1/mzp3LtGnT6gwtaMp4aBPC8dvOvdy9fB+nLxZga6llzUMDGdTJrXkNKcmD1CNw4YjyM/UoZCUAlX9fGxfw6gVePcEzpLJ1B6vWk76jsKCE+294jXH3DuQfz7epAnDtgrYiHNfjyy+/ZNasWRgMBpYtW2a0YL6WTGZmJq+88gpffvklqampgBLBHhYWxty5c5k4cWKtItKuhaNbr76SCf+lTG+gn78znz06GNuWUs+7tADSYyDtqJIaJf04ZMRBeY1oUSd/8AgG92DFz+IWBO5BYOfRIlco08YuoVdoIP9+9R5zm6LSQNqKcKzblwjA1KGB1a/l5+czd+5cPv300xbtADc1SUlJLFq0iE2bNpGZqdTP0el0DBgwgDlz5lwhIk0ZDy3kE7bxJGUV4auv4JlRAcy5tbe5zbkSK3tlR1bHwZdfMxggN1HZ7psRCxdPwsU4ZUeXvrjGtU6Kn8W1c2XrBC6BSrP3BjNFuAshKC4qNcu9VVQAfjqqfKuuEo7IyEgmT57M2bNnmTdvHvPnz28VDnBTEBAQwKpVq1i1ahXx8fEsWrSIzZs3ExERQUREBNOnT6dv37488sgjTbpPq//rajCw2fJ5Qvadg+M+0G0sDHsC3Frotw2N5rIYdB93+XWDAS4lQ1a8sj0467Ty75T9EPMtSMPlc7WWykrF2V/56eQPTh2U6HkHX3D0AStHk6xYOgd7X7cuR+HhDPK2JlKRW4rW2QrHWwOx6+9pdDtU2ieFhzMoO5ePrDBw4b/7+aP4MNNeebzVO8BNQVBQEGvXrgUgPj6ehQsX8t133xEVFUVUVFST+m71wtHd25GQMZPg2FfKB23UJ0qzclSixwfOgKCbzW1m3Wg04BKgtK43XXlMXwa555SVSvZZRWBykpSfp7ZCYcZf+7OwBXtPZXVi76E8+rLzAFt3sHUFWzfF72LjAjbOYOlQr1WMvZMNBkPtjzcLD2eQ+208slwRuYo02EUrAAAgAElEQVTcUnK/VRJIquKh0lSqx1eFMr4MuaUMLO/Mq9NfYOaSp3B2bvvVOBvKuXPn+N///semTZtISrpcNdHW1rZJCRZbvXBotVq48d9Kq9BD9BdwZD2cPwSnflGaRgeePaDn3xUhsW7+fc9NQmcJ7l2VVhv6Usi7APmpl3/mp0FBuvIzM155FFacQ7Wz/i8IRWytHJRmaVfZ7MHCBqmzpkJrQ7j2PF4dckn8Kg+ps6RC6NBLDXo0OBzuj7b8ynoDstxA3tZEVThUmkze1sTqLyVV2FpYM7XrHapo1OD3339n6dKl7Nixg5ycHEB5xBwUFMS9997Lv/71L9zd3RFNeCLR6oXjCrQ6CJ2iNFAC+CI+VKLC044pbftCJfai0wgIe0j52drRWSk+ENc6okkr9Ip4FGdDYSYU5yCLcyjKy6LwUhbF+TmUF+ZSUZKHobAQTXkGuookLA0lWMpSrCijF3rCAsrQxuz9S/cpZT/Ufttc1Sei0nSuNY7a+/jKy8vjo48+4ssvv+TYsWOUlSlpgSwtLRk0aBDTp0/nkUceMWqMR9sSjqvpOERpoHxQ7l8JcZuUb+DHv1GaRqfsYup2Gwx4WPEbtEHK9AaSsopJyNATn6EhPsOW+PQKErO0lJS7AkHV51rqNPg6WePlbo2rnSVONhY42VrgbGNJakIav22M5JkXx9GxgxOWGgMWQo+lRmC5LgVR8NeKjFrn9pumRcU4REREoCvMxNvOnaXYXXGsPY6vHTt28OGHH/Lnn3+SlpZW/bqrqysjR47k8ccfZ/To0Sa7f9sWjprYucOo55VmMMDpP+DQp0pKkYw4pe1+S3k0490but8B/acoz/9bGaX6CmIv5BGdnEt0yiWOnb9EYmYh+hq+CT8XG4I87Rne1Z0Ozjb4Otvg62yNr7MNbnaW11zGxjjp2P/eFhzzihlwY/crjhXebn2FjwNAWGhwvDXQJO9Tpe1TUVHBq6++ysKFC5k+7F7mj/gHmhrfTdrL+EpKSmLZsmVs3ryZkydPVmfN1el09OrViwkTJjB37lw8PZvnkXD7EY6aaDQQNEZpoMRbHF6nFIBKPQrn9ilt24tg7QQ+faH7eOj7QIv0j1zILSbiTBZHknOJTs4lNjWP8gpFJDwdrOjj58ytPb3o6mlPkKcDnT3sGh3r4uXrAkB56V/TPVf5MdRdVSrGICkpialTp7Jr1y4eeOABlixfjsXZMpZvisVQomeqs2ObHV+ZmZmsWrWKTZs2cezYMQoLC6uPeXl5ER4ezqxZsxgzZoxZis+1T+G4Git7GPKY0gByk5XVyKltkHkSzu5U2pZ/Kw5kr57Q9WboN0nZ+trM5BaVse90FrsTMtl7OouzmcqgsrPU0tvPiUfCO9PP34l+/i54O1nX0VvD0GiVlcilnMJaj9v192yTE1mlefniiy+YPXs2BoOBtWvXMmXKFGUV3B8O7D8NWPHsrEHmNtNoZGRksGbNGn788UeOHj1KXt7lgmkODg7ccMMN3H///Tz44IPY25u/8JwqHLXh7A+j/6M0gKwzyookYTtknbq8Ivn9ZdBZKzEZAcOh5wToOMzowXkGg+RISi6/xqazK/4iMRfykFIRisGd3Zg8uCPDurgT7O2AVmPaaHMXN3vcPR1JTKhlC7CKShO5OgJ8/fr1dOnSxdxmGZ24uDjWrVvHr7/+Slxc3BUrCjs7OwYNGsT48eOZMWNGi6wGqApHfXDrDDctUBpAXipEfw4Jv0J6rBIBnhELB1YCQomX8AqBTjdCr78r0d4NpLzCQMSZLLbGpPFrbDrpeaXoNILQABeeHNON8CA3+vg5Y6Ft3mWqEAIrG4tmvadK+6BmCvR58+Yxb948LCxa/1grKyvjhx9+4PvvvycyMpKkpCTKy8urjzs4ODBkyBBuv/12pk+fjp+fnxmtrR+qcDQGRx8Y8ZTSAMpL4MRmpW7H+UNwKQXO7FDa9oXKzi17b0VMAoZDj/FKOpGr0FcY2Bl/kZ+iU/ktLp28Ej02Flpu7ObBrb28GB3shZOt+SeSEIK8S40PHlJRqUlNB3hrjwA3GAzs3buXb7/9lt27dxMfH09ubm71cSEEnp6e9OvXjzvvvJNJkya1yhgUVTiMgYU19L5baVXkJEHMd4pvJCOmMjgvBeK3wW8LQGiVCG7XzpxwHM43JQP4PsmSiwVlONlYcHOIN7f29GJEkMc100ebi6AQX47sP2NuM1TaADUd4BMnTmTZsmV1fpBaW7SM+VBWVsb27dvZtm0b+/btIyEhgezsbGomjrWzs6N3796Eh4dz9913M2rUKLM4s41N6xcOfQsN/nEJgPAnlVbFxXhlVZK0h+y0c2y61IVvsodxXHZCh57RmgPcbbWPUfapWBZ1hOTeIIZCpxta1LZga9vmKRaj0rapmQL9Cgd4HXz6cPM7xWNiYti6dSt79uwhJiaGlJSUK/wSADY2NnTr1o3Q0FDGjh3LhAkTzFKdrzlo/cKREQsrRkKvu6Hn35RSri0VjyBiyx/mk9RRbMq5QJneQE9vOxZ0yOVO7V7csg5CzjkozIb8JKX+ecQy5VqNVtnR5eANLp2UFCo+/ZUAR3uPZn8r8hr5qlRU6qKl1gAvKSkhIiKCnTt3cvjwYU6ePMmFCxfIy8u7YhUhhMDJyYm+ffsSGhrKzTffzPjx41vEbqfmovULh2MH5ee2/yjNb5CyuynkrhYjIhUGya+x6Xyy5yyRZ7OxsdByb5gfU4YE0MOn6hvJhCsvKsyEM38opWozYpVHX0VZl4MVT/5c42QBFjZKzImdp5Ip1yVQqe3hVVk0yojxJ65uDuRmF5KVkYeb55X9Zn8fT1FkmpISS4DtYG9cJwTV3pFKu6OpNcCzv49nWYSSrG+6sGrw+EpJSSEiIoLDhw8TFxfHmTNnSE1NJTc3tzpVRxVCCOzt7enWrRs9evRg0KBB3HbbbfTp06dNPG5qCs0mHEIIf2At4IXysbJCSvnuVecI4F1gHFAETJdSHrpux/aeMHOHkoY85juI+R62vqA0v4GKgPQY36idTU2lqEzPZ5HnWLM3kZScYjo42/D8bd15YGDHup3cdu7Q+16l1cRggJyzSsR72lElfUreeUVoirIUX0padC0dCtBaKAJj6aCIjK2Lki3X3lMRHAdvsPeqTM/uDTautW4t9uuk1BsvKiylZp3F7O/jKYq4nP4ASfXvqni0fEw2RytZvHgxL730Eh06dGiUA7xqfB1ECT6dLq2qx5ftWH9OnDhBXFwcp06d4uzZs5w7d460tDQyMzPJz8+npKSk1n5tbGzw8PDAz8+PHj16MHDgQEaNGkVwcHC7F4hr0WwVAIUQPoCPlPKQEMIBiAImSClja5wzDpiLMigHA+9KKQfX2mEltVYWy0yA2O8g9gflwxWU+t89xiupRDx7mLS6XlGZnvURSXz05xmyCssYFOjKQ8MDuTnEC52pt88aDJCTqJSvzYyH7DOQfwEKMqAkF0oLQV8CBj3XzpR7GQkY0KIXOsqFBWXSghJpQZFeS4XOCj069GjRo8Gt5D/oqKW2swC/19pAMslWTl0V30w1RwEcHBxkQUEBDzzwAMuXL692gBsMBgoKCsjLyyM3N5fc3Fyys7PJyckhKyuLzMxMsrOzycrK4qUOD+Nkbc8cFN9CVc6qC3npDF5+b633FUJgbW2Nk5MTHh4e+Pv7ExwcTO/evRk2bBhBQUHtVhxaZelYIcQmYKmU8tcar30E7JBSfl75+0lgpJQy9Vr91FmSMvssxP2otJT9ymsunaD77RA8DvwHK1l1jUBRmZ51+5JYsVMRjBFB7jx5UxBhAa5G6d/olJeQk3KSoydPcuHCeQpysykvzsNQXoSoKEMny7GgAltKsBFl2FCKFeVYiXIsKccCPToqsKACLRXoylYhqF2Q/f6rCoe5aegHhbHmaOV5Tf6gSXp2Bxqh+YtwSCmZGjkfHx8f/P396dq1KyEhIfTp06fNOqeNQasrHSuECAT6A5FXHeoAJNf4PaXytSsGpRBiJjAToGPHjte/mWsnGP6E0vLTFN/AiZ9h/wrYt1QpZBR0C3S7FbqMadTupVJ9BWv3JrH8z9Nkt1DByC8pZ+vxNCITszmZls/5nGJyi8upMEjADuhWfa6lVoO9tQ4XW0u8Ha1wd7DCzc4STwdrvJ2sKcm4xIp5X/PsC7czdkJY9XUpz++qfRHT8kqnq9RBU+doZR/V8xSUQDd/f38sLS3R6XTodDosLS2xsrLCxsYGKysr7O3tcXJywsnJCTc3N1xcXPDx8cHb2xuxPrvW8SU0gh07djTxHas0hGYXDiGEPfAN8KSUMq+u82tDSrkCWAHKiqPeFzp4K6nTBzwMJXlw+nc4uUWJrTj6pRJb0XEoBN2kiIlnyHUfaUkp2RqTzqs/x3Euu6hSMLoRFuDSmLdlNPR6A3tOZ7ItNp2opBzOZRdRVHZlunMrnQYvR2sC3WwJ8XUk1N+Fvh2d6OBsW2f/OVlOrCkpQ39VokPbwd5X+jhqvK7SejDGHIUr56m/v7/MyMggJyeHtWvXctNNN9Vx9V/JHqz4OJyu+iaijq/mp1mFQwhhgTIgN0gpv63llPNAzYIYfpWvGR9rR2X3Vc8JYKiAlANKGdb4bfDbS0pz8IWuo5WVSOeRSsBeJcfPX2LRT7FEns2mm5c9ax8exA3dmn9bLCjPiXecyuSbqGQOJOZwMb/0ii9mDtY6+vo5MSDQlZHBHgwMcMG6kdlxr0eVA1zdVdV6MdUc9fLy4scff2TixIncfPPNPPPMM7zyyisNKi5UNY5eUceX2WnOXVUCWA3ESSnfusZpPwBzhBBfoDjeLtX17NQoaLSXiz7dtEApv5rwm9LifoTD6wEBvv3J9LuJN9MHsPGUHhdbSxZP6MUDA/1N7/S+ioT0fFbvPsue05mk5BRTFVahFYJAd1tCO7pwU4gXY7p7YalrPttcJwSpE7mVYuo52q9fP6Kionj66adZsmQJ27dv57PPPqN79+51X1yJOr5aBs254hgOTAWOCSGOVL72AtARQEr5IfAzym6NBJStfg81o32XcfSF0AeVVqGHC4eRCdv5PjqVhTu7UkgJM3S/MafDBZxKh8D5EeAbqtQGNyEHE7NZsfMM+85kkV+iPCYSgLeTNYM7uXHfAD+GdHY1+S4Ry8rVStr5HJPeR6XZMfkctbW1Zfny5dx222088sgjhIaG8vbbbzNz5sx618B+/ZcTADw3tv6Co2Jcmk04pJS7qcNNKpUtXo83j0X1RKvjgkMvXkyU/JF2kf5+DrwxsIigHAc4mw5/LIY/AJ0N+A9Ukhh2HKrEkFjW7S+oi70JmXzwRwJR53IoqaysZ6nTMCjQhYmDAhjfxwddM64oAOwcrOkQ4K4KRxujOefonXfeydGjR5k2bRqzZ89my5YtrFq1Cnd39zqvPZSkjjtz0/ojx02IwSDZsP8c//05DoOE+XeEMG1YYGXNi7HKSUXZkLQHEvdA0m7Y8V9AKhlxffqC/xBFUPwGKRHd9SA1t5gl206yNSadgkoHtJ2lllHBHjwS3onwIPP4UmqibWaxUml7+Pj48Msvv/DOO+/wf//3f/Tp06fRjnOV5kUVjmuQnlfCv748wt7TWYR3dee1v/fG37WWFYStqxJY2GO88ntxrpIm5Nw+OBcBB1dDxAfKMQdf8AuDDmHKoy3ffkoEN8pOqFW7z7I+MomUnGIALLSCG7q58/TN3ejrb96dWioqpkCj0fDUU08xevToKxznixcvxsrKytzmqVwDVThq4Y+TGTy9MZrisgr++/fe3D/Qv97PX7Fxhm63KA1AXwbpxyDloCIo56MUhzsAggynPrxcNomtuX6US+UeQZ72zBjRiXvD/NptVKtK+8IYjnOV5kMVjhqU6Q0s2XaSFTvP0N3bgaWT+tPV06FpneoslRVGhzAYPEt5rTCL/YejWLQzj2Pp1oDAhhLu1e7mOe0XOElnONkTcnsq9c09uoN7N6XuRwtBIwT5uWoxJxXjUeU4Hzt27HUd5z5OLWcetFdU4agkObuIOZ8fJjo5lylDOvKf20NMUjBmy/FUXv4xltRL5YANPk7WzBndlYkhtmgyXCGtE6Qfh/QYpYKgobLEpNCAcwB4BIN7kJL51q2rUknQ3sukubdqo2uIL/t3nmzWe6q0D+666y4GDhxY7Tj/5ZdfWLlyZbXj/J0H+pvZQhVVOIDIM1nMXh+F3iBZPjmU23r7GP0eW46nsvCHGNLylMJT/Ts689L4kCt9Fw6joMuoy79XlENWgpJG/eIJuHgSMk/B6T+gokYBKws7JfuvS6BSQMo5AJz9lbTyTv5KWhUjC4vNNYo5FR7OIG9rIhW5pWidrXC8NRC7/p5GvbdK28fX15etW7fy9ttv8/zzz1c7zoe69VHHVwug3QvHxoPJvPjdMTq62rJ62kAC3e2M2v/PR1NZ+FMM6ZWCMSjQhTfv7UuAWz3uo7VQMvl69rjydUMFXEpWUslnnVbSrGefUdqZP6D8qkdIOmtw8LmcLt3e63IqdTt3JbW6rYuSRt3KsdZU6vWh8HAGud/GIyu3DVfklpL7bTyAOrlVGoxGo+Hpp59m9OjRTJo0iaVPvE6X8S+wTCg7Df+Zizq+zES7FY4Kg+T1X06wYucZRgS5s3RSKE42ddTIaACxFy7xjw2HSMxSPsQbJBh1odFeXmF0HXPlMSmVuhy5SXDpvCIweReUVpAOFw4rKdbLCmrvW2gU8bB2BCsnsHIAK3uwtENa2FKusaJEWtKr/AIO3QtI+GZhdVp1h+iBaMuvfP4syw3kbU1UJ7ZKo+nfvz9RUVGceWknFkJHPJdX2+r4Mg/tUjgKS/X884sj/BaXztQhASwYH2K0lCG5RWXM/fwwu+IzAejj58T7E/sbRzDqgxDKKsLOXXHIX4vSAijKrC4AJQszyc+9SF72RUrzsygvysVQlIfILUBbnoGFoRgrWYI1pVhTznDKudHXAMd+r+4ypeyHWm9VkdtC68KrtBpsbW1x1tQ+h9Tx1fy0O+HIKynnoU8OcPhcDgvv7Mm0YYFG6ddgMLB4cxyf7k2iQko8Hax4676+LSJY72ou5pdy7Hwh8ekVxGdYEJ9hz+kMKCh1Bi7nAXK3t8TXWXHgezhY4WxjiZONBU62FuzdfIRT++N5/YNJWAoDFhqJ5dokyK/4y/20zup+fJWmo3W2ql0kHNvdx5jZaVd/8ezCMh78OJKTafl8MMl4TvCopBxmfHqAnKJyrHUaXri1O4+M6GyUvptKcVkFR1NyiU7JJTr5EkeSczmfW1x93MPBiiBPe+4O7UBXLwc6u9vRwdkGbyfr6+4qKzsQz4nkbLytbfDwVoIYC8dZXOHjABAWGhxvDTTZ+1NpPzjeGqj4NMovv1ZUXsJrm1cyeVChGnHejLQb4cjIL2HKqkiSsopYMXUAo7o3/ZmoXm/gn18eZvMxpQbFbb28efeB/s2ajfZqyisMRCfnsichiz0JmRxOzqG8Qkmd6+diQ/+Ozjw0PJA+fs4EeznUXfv8GgR0Uf5+BXnF1cJR9ZxZ3fWiYgqqxlHAd8cwlFWgdbZC392KiJ+P8ektt1RHnDckVbtK42gXwnEht5jJqyJJzyvhk4cGMqxL3YnU6mLHyQzmfHaYglI9LrYWrHxwAAMCzVPxLyOvhG2x6fx+IoPIM1kUllUgBPTydeLh8E4MCnSlr78z7vbGe2R0rUh6u/6eqlComAy7/p682//yhhAfIOqWKJ555hnefPNNtm/fzoYNG9SIcxPT5oUjs6CUyasiySwoZd0jg5pcztVgMPDUxqN8f+Q8Arh/oD+v/a1Xs6cGScwsZGtMGltj0jicnIuUEOBmy99COzC8iztDu7jhfI1YCxWVtoStrS3Lli3j1ltvrY44f+edd3j00UfrnypIpUHUWziEEPZSymvs4WyZFJTqeeiTA6ReKmbDjMFNFo3zuUXcs3wfqZdKcLa14LNHBxPi42Qka+smq6CUTUcu8HVUCrGpSkXPXh0ceeqmbtzS05tuXvbqRFFp8zz/7VEAXvt7nyterxlxPmvWrOpU7W5ubuYws03TkBXHYSHEXcApKaUeQAjRTUp5yjSmNY1SfQWz10URm5rHygfDmiwaGw8m8/y3x6gwSG7s5sHqBwc0Sx2MMr2B309k8M2hFP44kYHeIOnj58S8O0K4tacXfi5Nr/mhotKaOHOx8JrHrhVxPmbMmGteo9JwGiIcrsAHQA8hxEUgDggFuprCsKZgMEie3hjN7oRM/ndvX0Z392pCXwbmfH6Yn4+loRWCV//Wm0mDOxrR2trJyCthfUQSn+0/R2ZBGR4OVjwS3om7w/zo5tXExIsqKm2YqojzMWPGMHHiRG666aZG1ThXuTYNEY5zUspRAEIIP6AHUHz9S8zD67+c4Kejqbwwrjt3h/k1up+iMj13Ld1DfEYB7vaWfPfYcPzdTPsNPzo5l0/2nGXzsVT0BsmY7p5MHhzAiCD3Zq9rrqLSmlFTtZuOhgiHkxBiGBAjpUwBUkxkU5P46egFPtp5hqlDAph5Q5dG95OcVcQdS3dzqbicvv5OfDNrmMkeTUkp2RWfyfu/x3MgMQd7Kx1ThgQwbWig0XNnqai0J2pL1a46zptOQ4TDHngG6CmEsEJ5VHVcSvlvk1jWCE6m5fPs10cJC3Bh3h0hje5nV/xFHl5zgPIKyX0D/Hjjnr5GtPIyUkp2xmfy7m+nOHQuF18na+bfEcK9A/xwsDZe3iwVlbZEiK9jg6+52nFelapddZw3joYIxxAp5RkAIYQ1EAL0NIlVjSCvpJzZ66Ows9KxbHJoo4PwNh05z5NfHAEwakqSq9mTkMmSbSc5XCkYr/ytF/eE+WGlM34NEBWVtsSC8Y372FEd58ajzk9XIcRQIYSoEg0AKWWJlPKQlHKdac2rH1JKnvoymuTsIpZNDsXLsXEVwr48cI5/fnEEjYC1jwwyiWicuVjAjE8PMHlVJBl5pbz6t97s+PcoJg8OUEVDRcXEVDnO9+/fj6OjIzfffDPPPvssZWVl5jatVVGfFceDwAdCiFPAL8AvUso005rVMD7bf47f4tKZf0cIAxsZvf3x7rO8/FMsOo3gi5lDjB4FfqmonHe3x7N2XyLWFlr+77buPDQ8UBULFZUG8uQXh4GmVQKscpw/9dRTasR5I6hTOKSUjwEIIboDtwFrhBBOwB8oQrJHSvnXlKjNRJnewCub4wjv6s5DwwMb1ccHfyTw5taTWGgF3/1jOL06GC+oT0rJN4fO88rmWHKLy3lgYEeeurkbHg5qxlgVlcaQeqnEKP3Y2try4YcfMnbsWGbMmKE6zhtAvR0BUsoTUsq3pZRjgVuB3cC9QKSpjKsPKTnFaITg9Xv6NOp/9obIJN7cehIrnYYtT4wwqmgkZxfx4Mf7eearaLp42LN57ghe+3tvVTRUVFoQEyZM4OjRowwfPpxZs2Zx9913k5WVZW6zWjSN3V/6FrBDSjkXeNqI9jSYwjI98+7oQQdnmwZfuy0mjf98dxydRvDj3HC6GimwrsIgWbPnLLe+s5NDSTksuqsnG2cNbdRuEBUVFdNT5ThfsmQJP/30E3369GH79u3mNqvF0ljhWACsFkKsAwYY0Z4G42Cl474B/g2+7mBiNrPXRyEEbJgx2GjR2GmXSpi0MoKXfoxlYKAr2566kalDA9Fo1KWvikpLpspxHhkZqTrO66CxwrEIOAlIYGN9LhBCfCyEyBBCHL/G8ZFCiEtCiCOVbX59+u3gYtPgR1RnLxYwcWUEBglLJ4YyuLNx9nL/cSKDce/t4tj5S7x5Tx/WPDSwUSshFRVzYap5akxCA1wIDXAxWf9VNc5nzpzJm2++ydChQzl58qTJ7tcaaWxa9WellJlCCDvgXWBGPa5ZAywF1l7nnF1SyjsaYohFA9NwlOkN/G35XsorJAvv7Mm4Pk2vAlimN/Dm1hOs3HWWHj6OLJ3Uny4e9k3uV0XFDKzBBPPUmDw31vQ7n67lOJ8xY4bqOKcBKw4hxG1CiEghxElgmRBiiJSyEJhVn+ullDuB7EbaaTQeWLGP3KJy/t6/g1HiNC7mlzJxZQQrd53lwaEBfPePYapoqLRaWso8bSlUOc6HDRvGzJkzVcd5JQ35ur4MeAoYAqwAlgghJhp5K+5QIUS0EGKLEOKa4aFCiJlCiINCiIMXL16sd+ev/RzHoXO5dPGw4637+zXZ2BNpeUz4YA8xFy6xdFJ/Xr6r13XrdKuotBFMOk/rYva6KGavizJaf3VRm+P8999/b7b7t0QaIhwZUso9UsocKeVvKFtyXzSiLYeAACllX+B94PtrnSilXCGlHCClHODh4VGvzneczOCjnWewtdTy/ePDm2zs9rh07l62F73BwFezhnFHH98m96mi0gow6TytDzlFZeQUNa/D+mrH+U033dSuHecNEY6zQojFQoiqhPblgN5Yhkgp86oqDEopfwYshBBNLw6Okh79sfWHECg7qJqaQHBdRBIz1h6kk4cdmx4Pp7df81UBVFExJ6acp60B1XGu0BDhMAB/A5KFELuBBGCHECLIGIYIIbxFpddJCDGo0jajPEx89NODFJdX8MiITvTv2LTdGB/9eZp53x9nTHdPNs4airdT4/Jiqai0Rkw5T1sLVY7z7777jqSkJEJDQ1m5ciVSSnOb1mzUe1eVlHISQGVK9V5A38q2UgjRWUp53bJ4QojPgZGAuxAiBSUWxKKy7w+Be4DHhBB6lAJRD0gj/J/4KfoCe05n4etsw39ub3yqdSkl7/wWz7vb47mjjw9v39+vwTu6VFRaOuaap62RCRMmMGjQIKZNm8bMmTPZsmVLu0nVXm/hEELcJaXcJKUsBaKAKCGERkppqM/1UsqJdRxfirIN0GgUlOh5+qtoJcjvkUGN7kdKyX+3nOCjnWe4J8yP1+/ug5gKXLEAACAASURBVFYN6FNpg5hjnjaU4V1bzpOxKsf5W2+9xQsvvNBuUrU35CvzTCHEYAAhhFYI8TBwwjRmGYc5nx+iVG/gsRu70KkJW2SX/p5QXVXwDVU0ACgqKjW3CSrtlCfGBPHEGKM8ITcKGo2GZ555hoiICBwcHNpFxHlDhGMS8LoQ4t9APHAXSsr1Fsmp9Hx2nLyIu50lzzYhYOizyHP879dT/L1/Bxbe2VNNHVJJQtwFLK10dAhsOd/+VFTMSWhoKFFRUTz66KNt3nHeEOEIR9l++zjwsJTyLillhGnMajr/2HAIgCX3Nb7s6y/HU/nP98cYFezB6/f0UUWjBtIgsbaxxNKysckHVFQax7SP9zPt4/3mNqNW7Ozs+Oijj9q847whwnE38A7gDnwkhFgqhJhtGrOaxpbjqSRkFNDDx4GRwZ6N6uPwuRye+OIIff2d+WByqOoIV2kL9Da3AcagpLyCknKzlQCqF1UR50OHDm2TEecNqcfxsJRyIOAAjEcp5ORlKsOawovfHUcAH04Ja9T1F/NLeWz9Ibwcrfh42kBs1W/VKm2DY+Y2oD3h6+vLtm3bePPNN9tcqvb61BzXCCE6CyG0AFLhFJCIsruqRbHpyHmyC8sID3InwM2uwdfrKwzM+ewQ/9/efYdHVaUPHP+eSe8EkpCQUJPQIZSANJWmIIpIVxGVbtfd36LIrrq6qwK6ViwggiKLiwJSLPSidAJICRASIEAglTTSy5zfH3cSEkjPZEpyPs+ThzBz5865w7m8c+97zntSsvL48rGeeLrYV/6iBiglKQMbW3UVpigVqa+J86qc+cuAY2jDb0cKIZYLIS6grfw3sU5bVwPzN2nJqHfG1OyqfN5vZzl4MZl5Y7vQqZmaEV6eyxcTadvR39zNUBSr0KNHD44ePVpvZpxXJXDkAMFSym6AB9rkoKlAcynl5DpsW7UdOH+da6nZhAR40NzTudqv3xWRwJI9F3mib0tGdw+ogxbWL3YO6haeYnpDOvgwpEPNcpfmVDTjfN26dVafOK9K4FgIPAggpVwBLAH2SClj67JhNfH6Bm3tmbdHV/9qIy07nzlrThLs48qrIzoYu2n1il6vJz0lExs1YEAxg5l3BTLzrkBzN6PGRo0aZfWJ86qc+XnAW0KIG0KIPOB1IEsIcU0IYTFXHEkZOZyLz6BlE2c6+1f/FtNbG0+TmJHLfyaEqNLolUhLySI1OZNgdatKUWqkrMS5NZVqr0rgeAf4P8BTSmkPuBr+vAttTodF+HBrJABT+7eu9mu3n4lnzdEYnhkYSNeARsZuWr3l6FS7KsOKUhMTF+1n4qL95m5Grd2aOB86dCivvPKKVSTOqxI4TgPXpJQFQoi70WaQg5b7sJgJgL+cjMVGJ5jUu8Jai7fJyS/knxvDadfUjecHW04ZA0t2Iy0LAKEmRCpKrZVMnC9YsIB+/fpZfOK8KoHjdcBXCPE84CKlXGp4vCvwcp21rBpOxKSSmpVPz5ae2FZziOg3+6K5kpzNaw90xF4NL62SCxFxAAR1UItXKYoxlCzVfvHiRXr06MGSJUssNnFe6f+Uhnkbq6SUnxoWbil6/FcppUVcU328TbtN9fzgoGq9LvFGLgt3RDGkvQ8DglXNpepydnEwdxMUpV556KGHOHnyJH379mXGjBmMGzfOIhPn9eIrdtilFBxsddwZXL3lKT/ado6c/ELm3q9GUSmKYhlKJs43btxISEiIxSXOrT5wFOgladn5BPtUr2x6bFo2P4Rd4eHezQmsRcl1RVFM64GufjzQ1c/czahTJRPnrq6uFpc4t/rAkZypfZDDO/tW63VL91xEL2GWFY8HV5SGaHLfVkzu28rczTCJolLtJRPn586dM3ezrD9wpGfnA/Bo75ZVfk1aVj4rD17m/i5+NG9c/Rnm1iwhNpUP3ljLB2+sJS0l09zNUZRqy84rJDvPsqvjGpOLi0upxHn37t3Nnji3+sCRW6DH3dGWxq5VL0a44uAlMvMKmXV3mzpsmWV688X/suWno2z56SgvPPIFhQU1PwELC6u0arCiGNWTyw7x5DLLXI+jLpWVOE9OTjZLW6w+cOilpEU1rhqklKw6fIW+bZo0yCKGLdrcHECQmpxZo//8vX21zy3qzDWjtUtRlMoVJc4XLFjAxo0bzTbj3OoDB0CnapQYOXIphcvJWYzt2TCLGE6cfjcTpt3F2CcG8NbCydg7VH/2d7MWTQBIva5udSmKqel0OmbPnl0qcT5nzhyTJs7rRXnTfoFNqrztmqNXcbKzqXYyvb5oFdSUqS/eW6t9eHg649nElcjTV43UKkVRqqsocf7Xv/6V+fPns23bNlauXEnbtm3r/L3rxRXHwLZVm7+Rk1/IzyeuMbyzL66qJHiN6XQ6GjVxJSk+3dxNUZQGrWiN87Vr15o0cV4vAoeHc9US44ejk7mRU1Dvx4CbQttO/ly5mGjuZigN0LieAYxroLeayzN69GhOnDhBnz59TJI4t/rAoRNVL7S3JzIJOxtB32rc2lLK5tnElfTULFKuZ5i7KUoDMz60OeNDm5u7GRbH39+frVu3lkqc79y5s07eqx4Ejqpv+3tkEj1beuJsr25T1VZwJ20tjktR8WZuidLQJGfmFU/8VUq7NXE+ZMiQOkmcW33gsKli5Ei8kcuZ2PRq17NSytbWEDiioxLM3BKloXl6xRGeXnHE3M2waEWJ8+nTpzN//nyjr3FussAhhFgqhEgQQpwq53khhPhECBElhDghhOhRlf1WNXAcu5wCQJ826jaVMRQt4pSTlWvmlijGVFfnqWJ6Li4uLF68mLVr1xIdHW3UUu2mvOL4BhhewfP3AcGGn5nAF1XZqcjPhsTKI+nZuBsIAe193aqyW6USru5OeDV1J/K0mgRYz3xDHZynivnUReLcZIFDSvk7UFFrRwHLDet/HAAaCSEqH/4k9bDoLji4GCqIpGfj0mnZ2BkXNQzXKIQQuHk4W+xCM0rN1Nl5qpjVrYnzdu3a1Wp/lpTj8AeulPh7jOGx2wghZgohwoQQYfnCHlrdCb/NhhVjIT22zJ2fib1Be19347e6AXNwtCM2JsXczVBMq0bnaWKiGrptbjqdjhdffJFHHnmEpKSk2u3LSG0yKSnlYillqJQy1M7OHib9CCPeh0v74Iu+cHp9qe0LCvVEX88kuKlad8OY2nbyJyE21dzNUCxUyfPU29t4g1Ie69OSx/pUvRq2oomIiKBfv34sX76cKVOm1GpflhQ4rgIlB2cHGB6rkESCENB7Bjz1B3i2gh8eh5+ehhxtZnNyZh5Sgo+bWurUmHTVGQut1Bc1Ok+NaWRIM0aGqPXuq0pKyZIlS+jRowcXL15k7dq1LF26tFb7tKTAsQF43DBqow+QJqUs+75TCQX6EvfYvYJh2la4azac+B982R8u7SMxQxv5460Ch9EV5Bei16vy6g1Ijc5TY7qWms211GxTvqXVSk5OZty4ccyYMYO+ffty4sQJRo8eXev9mixTLIT4HhgIeAkhYoA3ADsAKeWXwK/ACCAKyAKqdC2l19+SnLWxg8H/gOB7Ye0MWDaCpA6vA+3wclWBw5iat/YmNyefy+cTaRXc1NzNUYygrs5TY/rLqj8BWDWrr6nf2qrs3LmTyZMnk5CQwHvvvcdf//pXdDrjXCuYLHBIKR+p5HkJPFvd/d4aN4o17w1P7YHNc0k9vBVoh2fOFaBxdd9CKUdjb21oc0EtFoNSLEtdnaeK6eTl5fHaa6/x3nvv0bZtWzZs2ECPHsadbmNJt6pqRF/RcFAHN3jwU/R9ngPAZtUk2P85qFsriqLUQ0UJ8AULFjBjxgyOHDli9KAB9SBwACRnVFKHpVmI9meLPrD5VfjuIUhTa0koilI/lJUAX7RoES4uLnXyfvUicOw4W8VCew98BCM/hpgwbdjuydV12zBFUZQ6lpyczPjx442eAK9IvQgc+y9cr/B5gTZsVA/Q80l4eg94tYM102D1VMgyz4LviqJU34w72zDjzjbmboZF2LFjB127dmXDhg0sWLCALVu24O9f5nxMo6oXgeN0bMUr0XkYCvKlZudrDzRuA1N+00ZfnV4PX/SD86Zf8F1RlOob2rEpQzs27FF8eXl5vPLKKwwdOhQXFxcOHDjA7NmzjTZqqjJWHzh0QnAlOavCbYqG4SbdKFHJ1cZWm+8xfZuWRP9uNPz6MuRVvC9FUczrfGIG5xMb7gJi586dK06AT58+naNHj9ZJArwiVh84HO10ZOQWEp+eU+42Xm7a0rJJZSXRm3WHWb/DHU/BoUWw+G64erSumqsoSi3NXXuSuWtPmrsZJleUAO/evXtxAnzx4sV1lgCviNUHjqLbUN/tv1TuNk1ctCuOhBvlBBc7J7hvPkxeB3mZ8PU9sHsBFBYYvb2KoijVZY4EeEWsPnA0dtGuJrafKX9klb2tDv9GTlxIzKx4Z4GD4Om90Gk07Hwblt4LSVHGbK6iKEq17Ny50ywJ8IpYfeDQCUFjF3uiKrnn2cHPjTOVJNEBcPKEsUtg3DJIvgBfDoBDX1W41oeiKIqx5eXlMWfOHIYMGWKWBHhFzN8CI+jTpjH5hZLNp+LK3aa9rzsXkjLJya9ieYzOY+Dp/dCqP/z6Ny15riYNllJQ1c9SUZRqKUqAz58/32wJ8IrUi8Dx0tC2AHy+u/zbSh383CnUS6ISqjEaw90PJq2G+/8DVw7enDSorj4AOBd+FZ1O0NTf09xNURqQ5wcH8/zgYHM3o05IKfn666+LE+Br1qwxWwK8IvUicLRt6oaXqz0nY9LIKyi7DlVIcw8ADl6s5mQ/IaDXdK1goldbw6TBKWrSIJCfV4CTiwNu7k7mborSgAwI9mJAsJe5m2F0RQnw6dOn06dPH06cOMGYMWPM3awy1YvAAfBQd3/0EpbsuVDm8wGezrTxcmFPZA2XsGwSCFM2weDX4MzP8HkfOLelFi1WFKUmwq+lEX4tzdzNMKpbE+Bbt241ewK8IvUmcDw/KAioeFjugGAvDlxIJremZcBtbOGuv8GMHeDcBFaOh40vQm7DnIyUl6eGKyum99bG07y18bS5m2EUtybA9+/fbzEJ8IpYduuqwcPZnq4BHsSm5XDgfNm1q+4M9iY7v5Aj0Sm1ezO/rjBzF/R/EY58q5UsubSvdvu0QpHhV2ne2nhrSStKQ1IyAT5t2jSOHj1Kz549zd2sKqk3gQPg3w91BuCNDeFlPt83sAkOtjp+PWWElS5tHeCet2DKr1oeZNkI2Px3yC9/Bnt9U1Cgp7GXq7mboShWpawE+FdffWVxCfCK1KvA0TWgEa29nImIv0FU/I3bnnd1sGVYJ182Ho+t+e2qW7XsB0/thdApsH+hVrLk2jHj7NuCFRbqSU5Mx9bOxtxNURSrYU0J8IpYfeAovGXt2Nce6ATAqz+VXctmTA9/0rLz2Xk2wXiNcHCFBz6Ex9ZATjosGQo734XCfOO9h4XJysghLSWLtp0CzN0URbEKu3btKk6Az58/3+IT4BWx+sARm1b61tDg9j74ujtwODqF02WMvBgQ5IW3mwOrj9TBZL6gofDMPug8FnbPgyVDIL5+JPFudSM9GwB7B5MtW68oALw8vB0vD29n7mZUWV5eHq+++iqDBw8uToC//PLLFp8Ar4j1ttwgJSvvthUA3x+vLRX77Mrbq9za2ugY2yOAHWfjuXS9ktpVNeHkCWMWw4TvtJnmi++GPR+Cvn7Nsr4Qoc3SD+rQzMwtURqani0b07NlY3M3o0oiIyPp378/8+bNs7oEeEWsPnA42towZ81JUrNulkwfEOxNF393LiZlse7Y7VcWU/q3wlan46s/yp7zYRQdH4RnDkDbYbDtn7B0WL0qmJibo92Gc3K2N3NLlIbmyKVkjlyy7Am4RQnwbt26ceHCBatMgFfE6gNHQGMnkjPzePOWcd2fP9YTAby2/hR6fenZ5E3dHRnd3Z8fw2JIysilzrh6a1ceY5ZAUiR82R/2fw76sme3W5PI8KvY2dng19w6vvkp9ceCTREs2BRh7maUKzk5mXHjxll9ArwiVh84nOxseHZQED8du8qvJ28Os23u6cz9Xf24kVPAP9aduu11M+9uQ16hnmV7L9ZtA4WAruO1q4/Wd8PmV+HbByC5jt+3jhUW6nFyccDRSV1xKEqRkjPArT0BXhGrDxwAzw0OIqR5I15efaJUEcMPJnTD1cGGlYeucCImtdRrAr1dGdHZj6V7oitcPdBo3P3g0VUw6nOIOwVf9NfKtVvQ1UdBfiHffbadL+b9wtKPNnM98fYhzUXs7G3JzMghs7zFsRSlASlKgJecAW7tCfCK1IujsrPR8cWkHjjY6pj1XRgZuVopDHtbHV88piWinlx2+LZbVq8Mb0+hXvLeZhNd9goB3SfBM/uhZV+tXPvyByEl2jTvX4nI09f476KdrF+5nx+W/sEvPxwsd9u2nfwpLNBz6bwRhzUrihUqmgFeMgEeGhpq7mbVqXoROACaNXLi00e7E309i9k/HkcaSp/fGezN8E5NSc7MY/bqE6Ve06KJM1P6t2LN0RhOXTVh0TQPf61c+4OfQuxx+LyfRVx95GTfHGAghKBraOtyt23bSbv8Pheu1ihRGqaSM8DrYwK8IiYNHEKI4UKICCFElBBiThnPPymESBRC/Gn4mV6d/fcL9OLV+9rz26k4vtx9c8TUwkd64OZoy5qjV9kVUfob8rODg2jsbM9bG0+j15twnQ0hoMfj2tVHizu0q49vR2qrDppJ9z6BfL3xJeYvmcridS/Q7Y7Acrf19vPA2dWBM8cvm7CFSl2r63PUGF4f2ZHXR3Y09duWUnIG+B133FEvE+AVMVngEELYAJ8B9wEdgUeEEGX966+SUnYz/Cyp7vtMG9CaB7r6sWDzWX4+cQ0AW1sd30zpjQBmLj9SKqfh7mjHK8Pbcyg6mRUHy6+sW2c8AuCxtfDgQog7oeU+DnxhtqsP/5ZehPRuU2nxQjs7W7rfEUjY3kgKjVW+RTErU52jtdWpmQedmnmY+m2LFSXA169fz7x589i6dSsBAQ2rgoIprzh6A1FSygtSyjzgf8AoY7+JEIL3x4cQ2tKTv6z6k71RSQD0bOnJ3BHtySvUM2rhnlL5jvGhAdzd1pt3fz1LdFIdTAqsvNHQY7I28qrVANg0B5bdpw3htWC972pH5o0cTv+prjrqCZOco7W1JzKJPZFJJn/fWxPgBw4c4JVXXsHGpuHVazNl4PAHrpT4e4zhsVuNFUKcEEKsFkI0r8kbOdrZsOTxXrTxcmXm8jBOxmj5ixl3BTKonTdx6blM+zaseHshBPPGdsHWRjB79XHT3rIqycMfHv0BHvoSEs9qVx97PoRCy1z3okdfbQ2UsydjzNwSxUhMdo7Wxqc7Ivl0h2m/VJVMgBetAV4fZoDXlKUlxzcCraSUXYGtwLdlbSSEmCmECBNChCUmlr2in4ezHcun9aaRsz1PLjvERcOVxNdPhOLr7sDOiEQ+2HJzNJWfhxNvjOzE4egUvvz9vLGPq+qEgG6PwLOHIPgebdb5ksEQV3bRRnNq7O2Gh6ezSpA3LFU6R6Fq56mlu7UE+tq1ay1yDXBTM2XguAqU/HYSYHismJTyupSyaCr3EqDMkC6lXCylDJVShnp7l38vvqm7I99N640EHltykOikTHQ6HRueG4CTnQ2f7IjivyXyGmN7+PNAVz/e3xxhlkvhUtyawsQVMP5bSL8GiwfC9rcsar0PGxsdrYKaEnXmmrmbohiH0c5Rw7ZVOk8tVXJyMhMmTCg1A3z06NHmbpZFMGXgOAwECyFaCyHsgYeBDSU3EEL4lfjrg8CZynZaWcmQNt6uLJ/am6y8AiYs2k9k/A183B3Z+Fx/7GwEf//pFJtPxRW9P/PHdiXIx5Xnvz9KTEpWNQ/RyISATg9pVx9dJsAf/9HKlkTvNW+7SmjT3o+UpIa5dG49VCfnqDXatWsXISEhrFu3rjgBXh9ngNeUyQKHlLIAeA7YjNbZfpBShgsh3hJCPGjY7AUhRLgQ4jjwAvBkZfuNTcth1vKw2yb3ldTZ34NVs/oigYmLDxB+LY2gpm78d/od6AQ8/d8jhEVrRdNcHGxZNDmUgkLJUyuOkJNvASOGnBvD6C+00VeFefDNCNjwAmSnVv7aOubgYEdOdh7pqWYOskqt1dU5ak2K1gAfPHgwzs7OxQnw+joDvMaklFb949wsWLZ85WfZ819b5YWEG7IiFxIzZN93tskub2ySRy8lSyml/PXkNdnylZ9l4Ku/yMMXrxdvuzU8TrZ85Wc5c/lhWVCor3C/JpWbIeWmuVL+s5GU7wVLeXKNlHrztW/v9nA5rMvf5dkTV8zWBqVqgDBppvO0Z8+eRjuOqIQbMqqSc70mIiIiZGhoqATkjBkzZEZGhtHfw5LUpj9YfRjt4OfOsI5NScrIZcgHu1m+P7rcbVt7ufDDU33xdLHnsSUH2XE2nvs6+zF/bBcK9JKJi/azJ1JL4g3t2JQ3RnZkc3g8/1h3sngmutnZu8Cwt2HGTnDzhdVTYOUESDHDHBRQS8cqJhfo7Uqgt/HWupdSsnTpUrp378758+dZs2aNSoBXwuoDR26+nkWPh/LhhG7ohOD19eE8tuQAeQVl37oK8HTmx1l9aeXlwvRvw1jyxwUmhDbn44e7oZfw+NJDbAnXch5T+rfm2UGBfH/oCv/Zcs6Uh1W5Zt1g+g4Y9q6W8/jsDvjjAyjIq/y1daBATQJUTGTb6Xi2nY6vfMMqKJoBPm3atAY5A7ymrD5wnEu4wevrTzGqmx9/vDyIAE8n9kRdp9fb28qtP+Xj7siPT/Xl3o6+/PuXM8z96SQjuvjx1ePaAJFZ3x0pXgDqb/e24+FezVm4M4rF5hymWxYbW+j7DDx7EIKGwPY3YdGdEL3HZE1o1rwJoK3PoSim8NUfF4yyCFtRAnz9+vXFJdAb2gzwmrL6wKETguX7L9H77e0k3Mjh99kDGd8zgLTsfB74dA9/XfVnmYlzZ3tbPp/Uo/iK4vGvDxHaqjHLp/ZGJwQvrfqT9zZHIITg3w915v4ufrzz61k+2R5pObetijRqDg//Fx5ZBflZ8M39sGYG3Iir87f29tVKP+TlWeYkRUW5Vck1wIsS4C+//HKDnAFeU1YfODr6uXNPBx+SMvMY9dk+nl15jDdHdWLFtN64Odqy9thVuv9rK/uibp+XodMJZg9rzwcTQjhyKYWRC/fg5mjH+uf642Rnw2c7o3hy2SF0Aj5+uBtjevjzwdZzzN8UYXnBA6DdcHjmINw1G06vg09DYe8nZrt9pSiWpr6uAW5qVh84hICvnujFyul30MjJjt9OxRHy5haOXk7l2D/u4aFu/qRlF/DokoNM//ZwmbmPMT0C+H5mH/R6GPvFPvafv84fLw/Cz8ORXRGJDP7PbnLy9bw/LoRJd7Tgy93neWNDOIXmKk1SEXtnGPwPre5Vy36w9TX4oh9EbjV3yxTFbG5NgK9evbrBlECvC1YfOIr0C/Li6GtDeWZgIFLCB1vP0fvd7YwM8WP9s/3wcrFn25kEur21pXjCX0k9W3ryywsDGNzeh7d/PcPs1cdZ/0x/erduTPT1LPq8u51T19L490OdmXlXG5bvv1Rq0SiL0yQQJv2g1b5Cwn/HwYqxkFAv52spSrlSUlKYMGEC06ZNo3fv3pw4cYKxY8eau1lWrd4EDgCdTsfLw9tz7PV7GNLeh+TMPKZ9G8acNSdZNasPT/RrRXZeIbNWHOG+j37nYmLpGc+NnO1ZNLknbz7Yib1R1xn52R5eGhrME/1akZFbwKiFe/lo2znmjujAW6M6sTMikXFf7DP/DPOKtB0GT++He9+GK4e1q4+NL5ok/6EodeHDid34cGK3Km27e/duunbtWmoGuEqA1169ChxF3Bzt+PrJXmz7y120berKmbgbDP3gd+JSs1n3bP/ixwb/ZzezloeRVSKxK4TgiX6t+OnZfrjY2/LoVwcpKNTzxaQe2Nvq+Hh7FA988gdjugfwzZReXE3N5qHP9nLkUooZj7gStvbQ7zl44Rj0mgHHVsAn3WH7vyxi9rmiVEezRk40a+RU4Tb5+fnMnTuXQYMGlZoBrhLgxmH1geNaajZRCWXXSgpq6saWv9zNosd64uFsx+bT8Yz+fC/tfN34cEJI8WMhb27hw60RpUZfdWrmwc8vDGD6gNZ8f+gyb/18mvfHdaVdUzdOXUun19vbAPjpmf64ONgycdF+vvr9gvlKsleFSxMYsQCeOwxth8Mf78PHXeH39yEn3dytU5Qq2Xj8GhuPl19YsygB/u677zJt2jSOHDmiEuBGZvWBIzkzj6Ef7ObxpYfYGZFQ5n/cwzr7cvQfQ5kzvD3O9jZsPB7L//14nP6BTXjq7jYAfLw9ih7/3sZvp2KLX+dsb8s/HujImqf74eZoy/P/+5MOfm483rcF2fmFTP76EG//cpofZvVhSActNzLt28Ncr6Twotk1bgPjl8GsP6BFX9jxL/ioC+yaD1nJ5m6dolRoxYFLrDhwe6UEKSXLli2je/fuREVFFSfAXV2NN8tc0Vh94Gjv585f72nLmdh0piw7zNAPtbIjmbckrXU6HU8NDOTEG/fy6n3tcba35ZeTcSz+/QKD2/kwqK03qVn5PL3iKH3f3V4qgd69hSc/P38nLw4J5peTsfx8Io5nBgbi7WrPzohEBr63mwe6+vGvUZ3Ye/46Iz75g33nzVyWvSr8usKjq7TyJS36wq534MPOsOlVSL5o7tYpSpWlpKQwceJEpk6dqhLgJmD1gcNWJ3hhSDB7XxnMRxO74eZgy+vrw+nz7nb+uSGc8GulZ4/rdDpm3R3IiTfu4e8jtACy+XQ8uyMTubutFyEBHlrF3RVH6PPOzSsQe1sdf7mnLRufH0CQjyuf7zqPp4s9wzv7klNQyPPf/8mG49f4dmpoeKKXqAAAEr5JREFUcW7ktXWnLHfUVUn+PeDR/8HT+6DDA3BosZYDWfkwnNsMelVORLFcu3fvJiQkhJ9++kklwE1EWOREtmoIDQ2VYWFhpR47ejmFZXuj2XwqjrxCPe193RjXM4BR3fzxdnMota1er+frPRf5ZHsUNwz/ybf3dQMpORuv5U6aujvwxgOdGNFVW4pASsnm8Dje/vUMV5Kz6dumMZeuZ3EtLQcbIXi0TwtsgG8PXMLP3ZG3x3RhUDufuv8wjCX9Ghz+Go5+C5mJ4OYHncdC5zHg1x1KlJjOzclnVO83mfrSvUyYeheZxxJI3xxNYWouNo0ccB/WCpfuVnTs9ZgQ4oiUMtQc713WeVoTmccSeGz1n8hCPZ97eLIp7RDTF7xIUFAQK1euJDTULIdnlWrTH+pl4CiSmpXHxuPXWH30KsevpGKjEwxs6824ngEM7uCDg+3NERZ6vZ7vD19h4Y4oYtO0VfaauNjj4mDD5eRsAHzcHHh6YCBP9G2JTqcjt6CQZXujWbgjipz8QkKae3AyJo28Qomrgw1T+rdm06k4IhMyGN3dn7kjOtwWuCxaQR6c+w2O/0+bQKjPB1dfbUnbVndCiz7kOvox6o63mPrSvdzfvT2payOR+TcHGQg7HY3GBKvgYQGsPXBkHksgdW0kz+bfAGAhLmTl5/Bb9iGmLXhR5TKqSQWOKnTIqIQbrD5ylZ+OxRCfnksjZzseDGnG2B4BdA3wQAhRvO2xyyn86+fTHLucigTsbXS4O9mSlKGV7nC01fFgt2b8fUQHPJztSbyRy4fbzvHD4SsIAc09nblgWOO8uacTfQKbsO7YVRxsbXhmUCBT+7fG0drKkWenQMQmOLcJzu+EXO0WoHRw5/RVFxwCgnHIfA6R73zbS20aOeA3p7epW6zcwtoDR+y8QxSm5pKK9sWkkeFOu+pfNaMCRzU6ZKFesicqidVHYtgSHkdugZ5gH1ce6u7PsE5NCfJxK942NSuPtzae5ueTseQV6BGAh5MdN3LyKZQggNBWnrwxshOd/T24kpzF57ui+DEsBiG0bYuCTcvGzni5OXDkUgoBnk7MHdGB+zr7lgpYVkNfCAmnuRG5l/MnD5J39RS+Nino8r6mvLRZwLw7TdtG5TbWHjhi5vxR7nOqf1WfChw17JBp2fn8ejKW1UdiiifwtfF2YVgnX4Z18iXEcCWi1+v5eu9Flvx+kfgb2lBbW51AJyCvUPv8AjydmNq/NY/3aUlseg6f7zrP6iNXkBKc7W1Iz9HyJ808HLHRCa6kZNOzpSd/GdqW/kFNrCqAFBTq+WLXeT7bFUVuvh6npDQGd/Dhb9mu2Gbk37a9+kZoGaw5cOzatQuX1cn4uXrzK9qXsRHYA6p/1ZQKHEZIusWmZbP1dDybw+M4cCGZQr3E192Rezs1ZVgnX3q3boydjY7ziRm8vzmCnREJ5Bju5dvqBAWG+SM2OkGoISAENHZi0e4LrD4SQ3Z+IU52OrINr2niYk+hXk9qdgE9W3ry0tBgBgR5WXwAycgtYObyMPadv86ILr68cHcgL414X+U4rIA1Bo78/Hxef/115s+fz8yBk5jbdyYvFGq3gRfiovpXLdSmP9gauzHWys/Dicf7tuLxvq1Izcpj+5kENofH8UPYFZbvv4SHkx1DOvgwrJMvH0zohpO9Db+diuXznVGcunZz1rWUkoMXk3n4qwO4OdoyvLMvG57tz46IBJbvv8TV1GzsdILrmYZ8iZ2Os7HpTP76ED1aNOLZQUEMaueDTmd5AaRQL3l6xREOXkzm/fEhjOsZQG7OzSuMopNXjapSjCEyMpJHH32UsLAwpk6dyvsff4yIzEIYRlWp/mU+KnCUoZGzPWN7BjC2ZwDZeYXsPpfIlvA4tp2OZ+3RqzjY6ujVqjH9gprw9uguBHq5sPiPC/wQFlM8IgvgRk4BP4bF8GNYDM09nRjZ1Y9Ab1d+PBLDoehkhNCWvpVo+ZLwa+lM+zaMVk2ceaJfK8b1DMDN0c5sn8OtVh68xB+RScwb04VxPcseJ+/S3UedyEqtFM0Af+GFF7C3t2f16tU3J/N1d8X+kJaH9Julbk+ZiwoclXCyt2F4Z1+Gd/Ylv1DPwQvJbD8bz76o6yzYFMECInB3tKVvYBOeHhhIB183NofHsSk8npiU7OL9XEnJ5svfteUum7o5cF9nX5ztbdh+NoHULO1be65hrZCrKdm8ufE07206y8TeLZh0R0uCfMw71FBKyaLfL9CrlScTezU3a1uU+islJYVZs2bx448/MnDgQL777js1mc8CqcBRDXY2OgYEezEg2AuAxBu57DufxN6oJPZGXWdzeDwAvu6O9AtqwjMDG3MlJYvtZ+KJSsigqIxW/I1cfjOUNHF3sCG0pSc6HYRFp6CXkG/YMCtfz7K90SzbG01Xf3fG92rBg12b4eFs+quQS9eziEnJ5umBgRafh1Gs065du5g8eTJxcXG8++67zJ49W1WztVAqcNSCt5sDo7r5M6qbP1JKLidnsScqiX1R19l5NoG1R68C0LKJM/d38UMPRMSlc/l6VvForPTcQsIMI7psBXi52iF0OuLTSxdKPHE1nRNXT/HP9ae4t2NTxvRszp3BXiabD5JoKNwY4Hn7PA1FqY2SCfCgoCD2799f4Qzwb6aoW1TmpgKHkQghaNnEhZZNXJh0R0v0esmZuHT2RiVx9FIqYZdSivMftjpBkI8rNkKQnJlLUkYeEiiQkFBiOKutDuxtBFn5N0e+FUr4LTye38LjsdUJ7m7rzciQZgxq74OHU91diXg6a0MfE9JzKtlSUaouMjKSSZMmcfjwYaZPn86HH35Y6QxwJ3t1FWJuKnDUEZ1O0KmZB52aeRQ/Fp+ew59XUjl+JZXjMamcuJJWXB/Lyc4GVwdbsvIKyMzTigoW6Cke5gtaAr3k4OkCvWT72QS2n00AIMjbhTHd/bmrnQ8d/dyNOjKrjZcL3m4O/HIylvGhKseh1I6Ukm+++Ybnn3/+9gR4Jb7bHw3A5L6t6qx9SsVU4DChpu6OxZMLAfR6yYWkDP68ksafV1I4fiWNM7E3h/YKwMZGUFgokZQOGmWJSsxkwZZzLNhyDhsBQT6uPNitGfd29CXQ27VWgUSnEzzZrxXvbY5gc3hc8TEoSnWVTIAPGjSI5cuXVysB/vMJrWK1ChzmowKHGel0giAfN4J83IqHt+bkFxJ9PZPI+AwiEzKISrhBVEIGFxIzS119VKZQQkR8Bu9tPsd7m88B4O5oS+dmbozo2oyhHXxp6u5QrUT39Dtbszk8jpf+9ydfTu7J3W29q3fASoO3e/duHnvsMeLi4pg3bx5/+9vfVALcCpk0cAghhgMfAzbAEinlvFuedwCWAz2B68BEKWW0Kdtobo52NrT3dae9r3upx/ML9Vy6nkWUIZhEJmRw+lo6F5OqHlDScwrYdyGFfRdS+Me68OLHGznbEuzlSp9AL4Z18qGdnwd2NrfXnHKwteHrJ3rx+NJDPLnsEE/2a8XMfq1qdbyKZamrc1RKydy5c5k3bx5BQUHs27ePXr16Gf8AFJMwWeAQQtgAnwH3ADHAYSHEBinl6RKbTQNSpJRBQoiHgfnARFO10ZLZ2egI8nE1zOe4eZtIr5ckZeZyLTWH2NRsrqZmcyExk/BraURfz+JGTj6VxZXUrAIOX07l8OVUPt0ZVeY2NgLsbASujrbYA+62FA8VtusVxJaEG3SJuY6/pwvXYlNpH6RuZVmbujxHz549y9GjR5k2bRofffSRKoFu5Ux5xdEbiJJSXgAQQvwPGAWU7JSjgH8afl8NLBRCCGntBbXqkE4n8HFzxMfNkW7NG5W5TW5BIXFpOVxLzeFychbhV9MIj00jJiWbtOx8cgv0VPYJF0ooLJDklFHEMN/bm22psG3hgVKPR8+7v8bHpZhFnZ2jubm51UqAK5bNlIHDH7hS4u8xwB3lbSOlLBBCpAFNgFILeAshZgIzDX/NFUKcqpMWm54XtxyrFfMS8+vFsdSnf5N2lTxvtHMUbj9Px40bZ9Tz9IenjLm3KmtI/aFcVpkcl1IuBhYDCCHCzFXx09jUsVie+nIcoB2LKd+vPp6n9eU4oHb9oexVd+rGVaDkBIAAw2NlbiOEsAU80BJwiqLUPXWOKlViysBxGAgWQrQWQtgDDwMbbtlmA/CE4fdxwA6V31AUk1HnqFIlJrtVZbgf+hywGW2o31IpZbgQ4i0gTEq5Afga+E4IEQUko3Xcyiyus0abnjoWy1NfjgMqOZY6PEcrfW8rUl+OA2pxLFa/AqCiKIpiWqa8VaUoiqLUAypwKIqiKNVidYFDCDFeCBEuhNALIcodFieEGC6EiBBCRAkh5piyjVUlhGgshNgqhIg0/OlZznaFQog/DT+3JivNprLPWAjhIIRYZXj+oBCilelbWTVVOJYnhRCJJf4dppujnZURQiwVQiSUN7dJaD4xHOcJIUQPI7636g8Wps76g5TSqn6ADmgTV3YBoeVsYwOcB9oA9sBxoKO5215GOxcAcwy/zwHml7NdhrnbWpPPGHgG+NLw+8PAKnO3uxbH8iSw0NxtrcKx3AX0AE6V8/wI4De04st9gIOqP6j+UN3+YHVXHFLKM1LKiEo2Ky6dIKXMA4pKJ1iaUcC3ht+/BR4yY1uqqyqfccnjWw0MEZa57qy19JdKSSl/RxvtVJ5RwHKpOQA0EkL4GeGtVX+wQHXVH6wucFRRWaUT/M3Uloo0lVLGGn6PA5qWs52jECJMCHFACGEpwaUqn3Gp8hRAUXkKS1PV/jLWcDm/WghhratZ1dW5ofqDdapRf7DIkiNCiG2ULAF709+llOtN3Z7aqOhYSv5FSimFEOWNjW4ppbwqhGgD7BBCnJRSnjd2W5UKbQS+l1LmCiFmoX1zHmzmNinm06D7g0UGDinl0FruoiqlE0yiomMRQsQLIfyklLGGy8OEcvZx1fDnBSHELqA72j1Yc6pOeYoYCy9PUemxSClLtnsJWn7KGtXVuaH6g3WqUX+or7eqqlI6wRKULN/wBHDb1ZQQwlNoi+cghPAC+lO6zLW51KfyFJUeyy33fR8Ezpiwfca0AXjcMJqmD5BW4nZpbaj+YJ1q1h/MnfWvwSiB0Wj34XKBeGCz4fFmwK+3jBY4h/bN/O/mbnc5x9IE2A5EAtuAxobHQ9FWXwPoB5xEG9lxEphm7nZX9BkDbwEPGn53BH4EooBDQBtzt7kWx/IuEG74d9gJtDd3m8s5ju+BWCDfcJ5MA54CnjI8L9AWazpv6E9ljkxU/UH1h4p+VMkRRVEUpVrq660qRVEUpY6owKEoiqJUiwociqIoSrWowKEoiqJUiwociqIoSrWowKEoiqJUiwociqIoSrWowKGUSQjRRAgxWQgxSgjRwjCz1MWI+/9SCNHfWPtTlOoQQgwRQnxXB/ttEP1aBQ7lNoY6QkMBP2A6sAmttLSbEd+mD3DAiPtTlOoIAY7VwX4bRL9WgUO5jZSyQEq5Skq5QEo5EugKtJVSxgkhlgkhRgohGgkhfhVCjAYQQnxvWN3tkBDikhDi/vL2L4ToAJyTUhYKIQINK6lFG1ZSSxZCnBdCuJvocJWGKQQ4ZliV8BshxDtlrQ1Sk34NtKrvfVoFDqUqAoDOQggfoAuQglaQ8d9Syp8M24QAF6SUvYFJwBsV7O8+tKsYpFYefg8wWUrZDTgBPCSlTK+TI1EUTVe0atSbgW1Syrmy7PpL1e7XDaFPq8ChVEpKGQ1EAOPRlu5diVZQch+AEMIR8AbeNLzkNFDm+ukGwzAEDoNOQNGayB0M76UodUIIYYe2LOz3wKtSyhXlbFebfl2v+7QKHEpVrQNcgWto6y0/ZTgBAToDkVLKHMPfe6BVDUUI4W24vRUghFgqhPAAGkkprxmedwIcpZQphlXUkqS2XKei1JUOaKXTC4BCKLOf2lFOvy5rWyGEM4Z+3RD6tAocym2EEDohxLNCiDeEEH6Gk6IZWnnmrVLKHWjfph43vCQEaCGEcDSMvHoT+BBASpkIXAb+A7wADEArQ12kIzfXMuiA9a5roFiPEGAf2joby4QQTW/tp1LKfMrp1+VsO4ib/bre92kVOJSyNEWrze+JduWwCHgNCOTm5fc7wKuGEVghwFrgINo3uS+klHsBhBCuaLcFCqSUGZTIbxiUvKTPBnoIIdrX3aEpCiHAKSnlOeAV4AchhBul+2nRdrf16zL6NJTu1/W+T6v1OJRaE0LsBmZKKSNuedwWWIz2TW0C2sn3AXCH4VuaophdWf1USrmrrH5dwbZHaUD9WgUOpdaEEDFACyml3txtURRjUf26fCpwKIqiKNWichyKoihKtajAoSiKolSLChyKoihKtajAoSiKolSLChyKoihKtajAoSiKolSLChyKoihKtfw/QR2C41IYuUAAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"f = plt.figure()\n",
"gspec = f.add_gridspec(ncols=2, nrows=1, width_ratios=[2,1],\n",
" # height_ratios=heights\n",
" )\n",
"ax = f.add_subplot(gspec[0,1])\n",
"for i in range(fn):\n",
" ax.plot(kxs * px / pi, emptykappas[:,i] * px / pi, c='k')\n",
" ax.plot(mykx * px / pi, \n",
" u.omega2k(u.empty_lattice_modes((mykx,0),3*eh)[i]) * px / pi,\n",
" marker='o', c=ccycle[6])\n",
"ax.axvline(mykx * px / pi, ls='--')\n",
"ax.set_xlim([0,1])\n",
"ax.set_ylim([0,3])\n",
"ax.set_ylabel('$\\\\kappa p_x / \\\\pi$')\n",
"ax.set_xlabel('$k_x p_x / \\\\pi$')\n",
"\n",
"\n",
"# branch cut image\n",
"mybetas = u.omega2k(u.empty_lattice_modes((mykx, 0), 3*eh))[:fn]\n",
"angles = np.linspace(pi,2*pi)\n",
"# branch cuts due to the 'lilgamma' function\n",
"gammacuts[:,:] = mybetas[:,None]*0.5*(1 + np.exp(1j*angles[None,:]))\n",
"# branch cuts due to incomplete Γ\n",
"gci = np.linspace(0, -1) * pi / px\n",
"Gammacuts = 1j*gci[None,:] + (np.sqrt(mybetas[:, None]**2+gci[None,:]**2))\n",
"\n",
"gr = np.linspace(0, 3, num=500) * pi / px\n",
"gi = np.linspace(-1.0, 0,num=400) * pi / px\n",
"gr, gi = np.meshgrid(gr, gi)\n",
"g = gr + 1j*gi\n",
"eps = np.empty(g.shape, dtype=complex)\n",
"for i in np.ndindex(g.shape):\n",
" eps[i] = emg(u.k2omega(g[i])).n#em['eps']\n",
"\n",
"ax = f.add_subplot(gspec[0,0])\n",
"for i in range(fn):\n",
" ax.plot(gammacuts[i].imag * px / pi, gammacuts[i].real * px / pi, c=ccycle[0])\n",
" ax.plot(Gammacuts[i].imag * px / pi, Gammacuts[i].real * px / pi, c=ccycle[1])\n",
" ax.plot(0, \n",
" u.omega2k(u.empty_lattice_modes((mykx,0),3*eh)[i]) * px / pi,\n",
" marker='o', c=ccycle[6])\n",
"ax.contour(g.imag * px / pi, g.real * px / pi, eps.imag)\n",
"#cax = ax.pcolormesh(g.imag * px / pi, g.real * px / pi, eps.imag, cmap = plt.get_cmap('plasma'), norm = Normalize(vmin=-5, vmax=5, clip=True))\n",
"ax.set_xlim([-1,1])\n",
"ax.set_ylim([0,3])\n",
"ax.set_ylabel('$\\\\Re\\\\kappa p_x / \\\\pi$')\n",
"ax.set_xlabel('$\\\\Im\\\\kappa p_x / \\\\pi$')\n",
"f.savefig(\"ewald_branchcuts.pdf\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [default]",
"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.5.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}