You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
218 lines
38 KiB
218 lines
38 KiB
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 67,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import scipy.constants as sc\n",
|
|
"import matplotlib.pyplot as plt"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 97,
|
|
"metadata": {
|
|
"scrolled": true
|
|
},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"pole positions\n",
|
|
"[[[ 0.01414214]\n",
|
|
" [ 0.01414214]\n",
|
|
" [ 0. ]]\n",
|
|
"\n",
|
|
" [[-0.01414214]\n",
|
|
" [ 0.01414214]\n",
|
|
" [ 0. ]]\n",
|
|
"\n",
|
|
" [[-0.01414214]\n",
|
|
" [-0.01414214]\n",
|
|
" [ 0. ]]\n",
|
|
"\n",
|
|
" [[ 0.01414214]\n",
|
|
" [-0.01414214]\n",
|
|
" [ 0. ]]]\n",
|
|
"pseudo charges\n",
|
|
"[0. 0. 0. 0.]\n",
|
|
"simulation time = 9.999999999999934e-07\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAADuCAYAAAAUXsqNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABeLUlEQVR4nO29eXxb1bU9vq4ky/Mkj4mdxHPiOLETJyYQKC9Aw0wSvgUaOpHX0tBCC7QMpS1Q6ASlpS2PvL7XiaHlB7SE0NAkDRBoHrRAQibPTjzPtmYPmofz+8M5hytZw73SvY5stD6ffMC2dHUl3XXPPnuvvTZHCEEcccSxMKE41ycQRxxxyIc4weOIYwEjTvA44ljAiBM8jjgWMOIEjyOOBYw4weOIYwFDFebv8RpaHHHID06uA8dX8DjiWMCIEzyOOBYw4gSPI44FjDjB44hjASNO8DjiWMCIEzyOOBYw4gSPI44FjDjB44hjASNO8DjiWMCIEzyOOBYw4gSPI44FjDjB44hjASNO8DjiWMCIEzyOOBYw4gSPI44FjHD94HHIAK/XC5vNBqVSCaVSCZVKBY6TrSU4jk8wuDC+6HHDBwlBCIHb7Ybb7YbT6QT/s1cqlUhISIBKpYJSqYwT/pMF2b7sOMHnCIQQOJ1OeL1ecBwHl8vl8zf6jyJO+E8U4gSfz7Db7XC5XD6hOJ/g/ogT/hOHOMHnI2hI3t/fD4/Hg5KSEvZ7p9MpmKSU7HT1B+KEX2CQ7cuLJ9lkgtfrhcvlgtfrhUKhgMfjifhYHMeB4zgoFDNFD0r2trY25ObmIiMjAyqViv2LEz4OijjBJQYhBB6Ph4XgCoUCHMf5hNsjIyPo6upCYmIisrOzkZ2djYyMDEbgcKCEpzcPegNxuVyM2CqVCgkJCSxTHyf8JxNxgksIQghcLhc8Hg8jITBDchqut7e3w+v1Yv369fB4PJiYmMDY2BjOnDkDtVrNCJ+enh6WlPTGEWiF93g8cLvd7LGU8CqVit104lj4iBNcIni9Xlb64pObwm634+jRo1i6dCmKiopY0q2goAAFBQXsMSaTCUNDQ5ienkZSUhKys7ORlZWFtLQ0waT0f31/wnMc5xPSxwm/cBEneJTg17b5qyj/7waDAVqtFg0NDUhLSwt6rKSkJCxatAiLFi0CIQQ2mw1msxkDAwOYnp5GSkoKW+FTUlIEn2MgwrvdbraNiBN+4SKeRY8C/rVtf1K4XC60tLTA4/EgPT0dy5cvZ3/jr/ZCX8tqtcJkMsFkMsFqtcLr9SI3NxdFRUVITk6OmJQ0aedwODA6OorS0tI44ecW8Sx6rIEmtYKF5GazGa2trSgvL4dSqYTJZIrq9TiOQ2pqKlJTU1FcXAxCCNra2uD1etHV1QW73Y60tDS2wiclJYk6tlKpBADYbDYmxOGv8PySXJzw8wdxgosEIQQWi4VlrwOF5L29vdDpdFi7di1SUlKg1+sRJlISDY7joFarkZubi+zsbHi9XkxPT8NkMqGjowNOpxMZGRmM8Gq1WtSxKeHpe3I6nXA4HOxmRglPxTtxwscm4gQXAVrbbmxsRE1NDZKTk33+7nA40NzcjPT0dDQ0NDDy+5fJ5IBCoUBGRgYyMjKwbNkyeL1eTE5OwmQyYWRkBG63G5mZmcjKykJ2djYSEhJmHSPYliEU4elrJyQksLJcnPCxgzjBBSBQbdsfBoMBHR0dqKqqQl5ens/f5oLg/lAoFMjKykJWVhaAmS0FJfzg4CAIIcjMzGRZepVK+KXAJzx9X06nE06nk702vw4vtL4fh/SIEzwMAtW2FQoFvF4vgJlVvbu7G2azGevWrQu495WT4EKPq1QqWbgOAG63GxMTEzCZTOjr6wMApKamwul0wuPx+KzYocCXzvLPhxLearXC4/EgJycnTvhzgDjBQyBYbZsS1mazobm5GRqNBuvXrw8alspF8GjCYJVKhZycHOTk5ACYIfzo6CgmJiZw4sQJKJVKFs5nZGRETHiLxQK73Y709PT4Cn8OECd4APBD8kC1bYVCAb1ej+HhYVRXV0Oj0YQ83rkI0cVCpVIhOzsbU1NTWLlyJZxOJ8xmM7RaLbq6utjfqcpODCkVCoXPCk/38HzC+2fp45AGcYL7IVxt2+v1YmJiAlarFQ0NDYKy04EIHouk5yfZ1Go18vPzkZ+fD2AmgUgTdlNTU4Jltf6Ju0CiG0IIHA4HS9rRTrm42030iBOch3ByU4vFgubmZqhUKlRVVQkuPcUimcUiMTERhYWFKCwsBOArq52amkJycjIjfGpqKvvswol5AhHe6/XCbrez38VbYyNHnOAILzcFZjrA+vr6UFNTg5GREVGEjYUkm9BjCSVPIFktTdhZLBakpqYiKysLLpcrYEkuGIQQnn5HycnJccKHwSee4Py+7UCrNr8D7LzzzoNKpcLo6GhMEDxWLmyO45CSkoKUlBQUFRX5yGp1Oh1cLhcmJibYCu+vHwh3bH/CT09Po7+/H9XV1QDiK3wofGIJThNpAwMD4DgOixYtmnVhTE1Nobm5mXWA8ds/aZlMCOZLiC5mBfeHx+PBqZ42nBnuQVJCIs6vrkdxcTG8Xi9UKhXS09NhMplw5swZOBwOpKensyy9WFktANbnTld4KrGlf4sTfgafSILzQ3Kv1+tjhUT/Pjg4iOHhYdTW1s7qABNL2PlC8GhwqqcNxzubUJCdB6fLiQMfvYP/d+FVIIRAqVQiPT0d6enpWLp0KZPVGo3GiGS11OgCCO52Eyf8DD5xBPcPyZVKpY+dksvlQmtrKxISEnDeeecFrP9SAweh+CQQvGukFwVZuVCrEqBWJWDCOgWd2QA1md2YwpfVAvCR1Q4PD8Pj8fio7Pz38HyC+0MI4T9J9lafGIIHq21T6yPAtwOMZosDgf8cIfAnOCEEIyMjUCgUohtB/HGukmz+SE5Mhs1hhzph5r14CUGCKgHE5Q57TL6strS0lDndmEwmDAwMzJLV+kdcofBJd7v5RBA8VG2b+pn19vZifHycdYCFQjQhutPpRHNzM5KSkqBUKjE8PAyv18v2o1lZWaJVY7GADcvX4MBH/8SUbRperxeLcwtQnFuIwYFB0eepVCqh0WiYgMhfVut2u5GQkACj0YjMzEzBnxcQ2u1mdHQU2dnZSEtLWzCEX/AED1fb9ng8GBkZQWFhIc477zxBKqpIk2w0QqisrGQtnhzHwe12w2w2w2g0oqenh+nGNRqNaNVYNIhmBS/IzsP1F14J/YQRKqUKxbmFUClVUR2Twl9WOzo6CqPRCL1ej+7ubh9ZbWZmpqjPi39NTExMICsra0HZWy1YggupbRsMBnR2ds5yWwmHSPbUNpsNHR0dLELgDz5QqVTIzc1Fbm4ugJlV3mg0MtUY9WbTaDRISUmJ2QssKzUDWakZPr+TguD+4DiOJewA+MhqOzs7kZCQEJGs1uv1zpLKznd7qwVJcCFyU9oBtmLFChgMBlHHF7OCu91uZtvU0NDAwslQF4VarWaqMb6IpKenB1arFenp6cjOzobb7ZZ8Dy415CC4/x48nKyWb08dyrwyUPIuUEjv73bj3zgTS4RfcAT3v9v6f9h2ux1NTU2sA2xyclJUuE2PK+Q5tI6+bNky2O12n72imCSRv4iElpgMBgP0ej2MRiM0Gk3AjLNYSH1xykFwQkjIVdlfVktvkIODgyFltULaZAOZXwSyt+rp6UFpaSlSU1OjfbtRYcEQnBCCiYkJcByHxMTEgBcVDeH4HWCRhNtCnjM0NISBgQHU1tYiNTUV/f39s843EtDwND09HV6vFykpKVCr1TAajRgYGAAAZGVlQaPRiGrzlAtyreBiDCqSk5ORnJyMxYsXB5XVZmdn+wyOEIpAhHc6nXj44Yfx2GOPidr6yYEFQXBa2x4cHERmZibzGef//fTp07DZbLM6wMQmzMI9x+PxoK2tDYQQJm2VE7TURo0cXC4XzGYzdDoda/PUaDSChinItdrKHaKLQaCIyGKxwGQyweFw4NixYz7mlWJktfT4SqWS3TjONeY1wf2tlKh0kQ/aAVZYWIgVK1YEFF1EEqIHWoEtFguamppQXFyM4uLic7IXS0hIQF5eHrONcjgcMBqNrOsrJSWFET4aq2WhkIvgUlUWOI5DWloa0tLSMDY2hvXr1zPzSr6slhI+MTFR0HEtFgvS09MlOcdoMG8JHshKSalU+pCV3wGWmZkZ8DiRruD+BB8dHUVvby9WrVrFFFpC3kO0F3+4UD8xMdGn68tqtcJoNDKrZXrxajSaebOCh9uDR3Nc/haIymqnpqZgMpnQ1tYGl8vlI7oJJlKyWq2ihlPIhXlJ8GC1bSpaCdQBFgyRruB8T7aOjg44HA40NDREneQSex5iH0+91ZcsWeJz8ba0tMDhcCAhIQEZGRmijRiDIZpwOtQx50oboFAokJmZiczMTJSUlDDDD7PZjKGhoaCyWkKIWCPLJQD+BKAAMwNHfkcIecrvMRyApwBcDcAKYAch5ESo484rgoerbSsUClgsFp8ZYEJkkpGG6DabDY2NjSgsLER1dXVMlUeEwP/i1Wq10Ol0MJvN6OvrYxJSmrCLhFRyrLZy3DQAYTdMfs4jmKx2z549UCgUmJqaEhzNAXADuIcQcoLjuHQAxzmOe4sQ0sZ7zFUAKs/+2wDgf87+NyjmDcHD1bapUkyv12P9+vUhZ4DxIbZxhD7HYrHgxIkTqKmpYdbE8x00AVVaWgrgYwEJnX5K68kajcanvBQK8yVEj7Sq4S+rdblc6Ovrw/79+3HFFVcgJycH+/btE/L6owBGz/7/FMdx7QCKAPAJvhXAn8jMyX7IcVwWx3GLzj43IOYFwcONCaIdYC6XC4sXLxZMbkB844jX62X11I0bN0bVKBLr8BeQ+JeXhGSbYz3JJvUxExIS8JnPfAb/9V//hQ8++AA2m030MTiOKwGwFsARvz8VARjk/Tx09nfzk+BC5KbRzgATc/FRkUxaWhpycnJigtxz2U3mX0+2WCwwGo0s25yRkcEy9Py96HwguMfjkeyYTqeTvf8IymxpAF4FcDchZDLac4lZgoezUiKEoK+vD1qtlum7DQaDT2+3lKCTS1asWAG1Wo3e3t6ojifFRX8u9/z88hLNNk9OTsJoNLLJKVlZWWxbJSXkStxJJQqanp6OqAbOcVwCZsj9/xFC9gR4yDCAJbyfi8/+LihijuC0/DU6Oor8/PyA5KYtl2lpaT4zwPzLZFKdT09PDwwGA5tcMj09veANHMTCf1QS7ZAbHx9HU1PTLMFNNKulHHtwKVfw6elpUdtEgGXI/wignRDyyyAPex3ANziOexkzybWJUPtvIMYIztf19vb2BjRdCDUDLJKMeCg4nU40NTUhPT0d69evZxeA1K8TDWLF8MEftEOuv78fa9asgdvt9mkASUpKYoQX2yEn1x5cqhXcarWKJjiACwF8EUAzx3Gnzv7uewCWAgAh5H8BHMBMiawLM2Wy/wx30JghOL+2HUiRJmQGmJTEo8KGyspKlmSiiES/LgfkyE5LDXrT4DeA0BIj7X/nd8hpNJqwarFYTrIBMyo2sSIXQsi/AIT8Qs9mz+8Qc9xzTvBgVkp8+HeABbuwpQjRCSHo7+/H2NgY6uvrAyZJIimt+WNqagoWiwUajeacN4TwMRdKNr4evLi4GIQQH7UYHXVMM/T+ghE59uBiBi6GQyQhulw4pwQPV9sGAneABQNVskUKl8uFlpYWJCYmhnR3EVta44Pv2Jqeno6+vj62P9VoNCH7lecrwr0fjuN8ZpvzxSP9/f3gOI4JbjIzM2UTz5zLPbhcOGcED2elRAhBe3t7wA6wYIgmRJ+YmEBraytKS0uxaNGisK8TyQrudrvR1tYGjuOwbt06dlHRhpCBgQF2cVDCC21ukAJylLQiQSDxCN+xxWazYWBgQNIbopRJNqvVGhOdZMA5ILiQ2rbFYoHVakVRUVHADrBg8LdAFno+VChTV1cn6IuJZAUnhOCjjz5inWbUjx2Y3RBC68s0XKWrVyBDxljIBcgN/w65I0eOQK1Wsxsi7eeOpkNOyiQbFQHFAuaU4OFq28DHHWDJyclYunSpqC9LbPKLrqgejwcbNmwQPGFD7OvodDpYLBY0NDSwvu1Qx+bXl2m4ShNSdIwvNSCUEvPlZqFQKMJ2yNEMvVAxktQhergocK4wJwT379sORG7/DrDjx4+L/tDF3Aymp6fR1NSEpUuXwuFwiH4dIWSgNXSj0Yj09HQxjQcM/uEqDecHBwdhNBqRlJQEQgg0Go2oEUDBEAshuhgE65AzGo1siAJ//x6sw8vj8UimTIywTCYLZCe4f0ge6AIKNAOMhttyOKLQKGH16tVIT0/H+Pi4aBvkcHC5XGhubkZqairWrVuHjz76yOemECmR+OH8wMAAa4/t6Ohgvcp09Yql7Pxcgd8hR7u9zGYzTCYTent7fbrB+B1yUpfJPhF7cI/HA5PJxDqPAiXSgs0Ai2Q/HQ5erxft7e1wuVw+feJSC1foDausrIyJdaQorfmD4zg2xtc/nO/t7WXhPPVXD3dTiZUkm5RQKpU+nupOpxMmk8mnQ06j0cBms0m26i74PTgNyamkdOPGjbMeE24GmNQEt1qtaGpqYmTgX8hSEpw6uwQaWiiXkIQiUDhvMpmYXVNqair7uxTh/HyEWq1GQUEB8+2jghu6yo+Pj7ObYqSf0YImOL+2HSzkMZvNaGtr81nh/CGlrnx8fBxdXV1Be7elIDg1drTb7QFdZORawUPBXz1Gs/P+4Tx1b5kPK7jUn2FycjKKioowPT2NwsJCKBQKmEwmNvWUL7gR6tZDlXmxAEkJLqS2TTvA1qxZE1LOF41ohb6+1+vFmTNnWAY7WBIlWoI7HA40NjYiNzc3ZFnvXGapQ2XnqXuLSqVCampqTBNdLrsmqmRLS0vz8WOjghvaIUfJHmomWqTdZHJAMoLTlRtAwC8gWAdYMEQaolOyulwuNDY2Ii8vD8uXLw95wUazulJ55YoVK0KWruRYwaOBfzjvdDrR1dUFk8mEo0ePxmw4L5ddU6Abh78lNe2QozPRguU4xDqqchz3DIBrAWgJIasC/H0TgL0AaI/yHkLID4UcWzKCU9FKoItYr9fj9OnTATvAgiEagut0OnR3dwuSt9Jzj0S40t/fj9HR0aCa9UDPkRpSHVOtVrOGj8LCwrDh/LmCXI6qQrTo/jPk/HMcKSkpOHXqFLxer9ib4nMAdmHGdDEY3iOEXCvmoIDEIbp/fZiG7D09PUE7wIIhUlWa3W5Hf38/1q9fL1jmKTZE93g8sNvtmJqa8pk3FgpydKDJFUYLCefp6h4qOy/HDU2uED2S4wbqkHvnnXcwPDyMdevWoaGhAb///e/DHpcQ8u5ZmybJIdutmHaAKRQK1NbWig7zxBKc9m5zHIdVq1aJ0nCLIbjVakVjYyOUSiVWrZoVTQVFrLSYhkKw8wsUzvOHKQQL5+eLXRMQvRaddsjdeeed2L17N44dO4aOjg4pz/UCjuMaAYwAuJcQ0irkSbIQnN8BNjAwEFECS6FQ+IzYDQW6D66qqsLIyEhELqlCzlGn0+HMmTNYtWoVWltbRV3A84HggLCowH/6KT+cdzqdPlbLchB8rvbg0UClUolaAMLgBIBlhJBpjuOuBvA3zFgnhz8Pqc4ACNwBNjw8zIapi4FSqYTdbg/7ejQrT/fB4+PjokP7cAQnhKC7uxsmk4m9L0rYaAgeq5lqMQgVzvf29sJms6Gvr0+w2CYc5NqDB1NZioVMNtGTvP8/wHHcbziOyyWE6MM9V1KCd3d3Izk52adUFGmyLNzzqBQ0OTnZJysv9TBBl8vFnFTXrVvHXkfsiizXCh5rlk38cJ5+dklJSZKJbeZyqkkkcDgckrf4chxXCGCcEEI4jjsPgAKAoKH2khK8oqJiFlHkIPjExARaWlpQXl4+SygjJcEDSU79nyNU7y1Xki2Ww35ydnwPP5ynnV+nT5+Gw+HwaYUVkp2PdYJHUgPnOO4lAJsA5HIcNwTgBwASAObFdgOAr3Mc5wZgA7CdCPziJSV4IKKoVKqIQ3T/Y/G162vWrAn4QUaigAtUJqMNKYEkp/Q50azgBoMB4+PjMVF6opA6vPQ/XqDOLyoRpdl52gobLJyXaw8uFSJxcyGE3Bzm77swU0YTDdmvqmjq2fznud1utLa2QqlUBtSuB3ue0NfiDxOkq0uowYViIwVKcH7eoLi4mM0Bo6FtTk6O4LFAsY5wZOSX24CPG0GGh4dZXZkKSajOQI49uJRRUCy5uQAxTHD+82jv9rJly1BUVBTyedGE6LS0l5eXF9ZJJpIV3OPxoLm5GUqlEuvWrYPH45k1x7uvr49pmenFP1cTS+VewcOB3wjCD+fp5JSsrCwoFIqYHmYYS62igAxCl1kvoFKFzYYHAiX48PAw+vv7We92OERKcKvViuPHj4eVnFKIJbjH40FHRwdKSkqwZMkS1nFH4W/bRKeEDA0NAQBb3f1D11jfg0dKnEDh/MTEBLNpMpvNbHWPthwn9VSTWOkkA+ZoBY9kDw4Ak5OTrJ1U6B5VbMRACIFOp4PBYMD5558veJaUmBsJbUMsLS3FkiVLwj6e4zgf0wKXy8UcSiYnJ5GamoqcnBxmohGrkDIioPtzm80GjUaDgoICNkjh9OnTSE5OZhGP2HlgUps9fOIILjZEp2oxqoITc5GIEch4PB60tLTA4/EgPz9f1IUhdAUfGBjAyMgIFi1aFHHolpCQ4BO6WiwWGAwGjI6OwuPxwOFwICcnB5mZmVGPBDqXIbrQYyoUirDhPN/ZJtziIKUn+icyRBdD8LGxMXR3dzO1mNgLROjKSm8iS5YsQXp6OguFpXod6h7j8XjQ0NCA3t5eScJpvrBErVbD6XQiNTWVqQeTkpKQk5MT0UomNeSSqvoTNlg4bzQa0d/fz1b/YOH8QvVEB2JoBafZa5vNhvPOOy/ixJKQMhklw6pVq5CZmYnJycmISmvBCEv7w/Pz87Fs2TKmkpJjv6xQKHw6nAIlpnJycgJaLvtjPqzgQsjo3+ZJs/PBwnmpV3Ah27C5QkzswW02G5qampCfny/KBz0QQpXJCCHo6urCxMSEjwFEJIm5YISlIhz/ZN1cCV34I4FonZlaLickJLBkndiBf5FALoKLPWa4cJ4q6txud9R6BJvNtnBD9EAIt4LTBo6VK1eG9QwXgmBkpbLJ9PR0rFu3LmpPtkDPGRkZQX9/P5tXzse5UJ3515ntdjsju81mQ0ZGBvs7tWySEnLuwSNFoHB+cHAQOp0Op06dAsdxPq2wYl9rQYfoYvbgdDU1m82ierfDIVCIPjk5yaSt1GyPj2hXcGoNRZtsAq0CsSArTUpKwuLFi7F48WLmH24wGDAwMMDMOqjPeqw2XkgtVVUoFEhOTkZubi5KSkpYxYKOORabnV/wWXT/CzmQDNThcKCpqQnZ2dkhp4XS54r5Qv1D9HCSU/qcSFdw2oeenZ0d0hoqFgjOB98/vKysDE6nE+3t7dDpdBgeHkZaWhpL1kU6EGA+ENz/mP4VC7HZ+VgaegDMQYju/wUbjUa0t7dj+fLlLDEUDDS8F0twOveL9ieHq6NH4pfGcRysViu6urpQUVExa4Z4oMfHcjeZWq1GSkoKcnNzkZWVhenpaRgMBrS0tMDr9TKNOH9YgJBzkyNEl/qYwZJsobLzAwMDADArnBe7gn/5y1/Gs88+q0VwPzYOwFMArgZgBbCDEHJC6PHnrMOBEILe3l7odDrB9k2U4GIy6rQO/tFHHyE/Px/V1dVhL4hIPNksFgvMZjPWrVsn6Audb5ZN6enpSE9PR0lJCdxut8+wABq25uTkhPwe59MKLiS55p+dd7lcLDs/NTWF1157DUajEQaDARUVFYJee8eOHXj22WevRHA/tqswY+5QCWADgP85+19BkD1EB2a+6BMnTiA1NVWQoypFJKHz5OQkJiYmUF9fL3hAn5jXobmDyclJlJWVCb5bR3ITmWsEuwGpVCo23ZMftlJDxmDTT+cLwT0eT0Q5oISEBOTn5yM/P5/lL/71r3/h0UcfxdjYGP7973+H3bdffPHFAGAM8ZCtAP50tj30Q47jsjiOW0QIGRVyjrKv4BMTE7BYLKioqAiY4AoFMSo46nI6NjbGpJxCIfQidLvdzPyhuLhYloGF5xpCoh1+2Epnf9HsvFqtZmGrlDO3KeTeg0cKjuOY6vLvf/8764WXAEUABnk/D5393bklOCEEAwMDGB0dZQkbsRBKcNpKqlKp0NDQgCNHjkRyyiFhsVjQ2NiI0tJSNvgvknZRqXGubxr+s7/oKKCenh5MTk4iKSkJycnJgiSjQhDrUQEt48VKj4AsBHe73WhpaUFCQgIaGhpw6tSpiCaFCiG4xWJhY4DDtZJGClqrX716NRsBHAuWTbGYwKKjgIqKijA8PAyLxYLJyUkmGaV797S0tIheS64QXQolm0w322EAfGlc8dnfCYLkBJ+amkJjYyNKSkqwePFiAB+r2cTuc8LJTqnklE88KUETgwaDYdboI9rfLRT+BKdNIwvF3CEYUlNT2Y3X6XSyuvv09LRPz7vQUlyshuh8SPx9vg7gGxzHvYyZ5NqE0P03IBPBpRoFHEx2GkxyKiWoOUNiYqKP2SL/3CJdwWkXm91uh9PpZFuYnJycOTN3CAQ5tOj8z02tVvv0vFOhTXNzMwghPj3vwQgX6yu42M/v5ptvBoAPENyP7QBmSmRdmCmT/aeY40tO8KKiolmklNJ4kQpLMjMzZ0lOpQK/06y4uDjgYyIN0e12OxobG7Fo0SKWdKQ156amJgBgZI80jI0VhLrgOY5DRkYGMjIyWM87LTnRnndKeH7kF8t7cJvNJrqD76WXXsJLL720KNjfz2bP74j0nGQpk816kSiMF/kEn5ycRHNzMyorK8MKSyKFwWBAR0dH0FHDFJF4stntduYao9Fo4HQ6Z9WcXS6XTxibkZHBFGWB2iSltk2WEmLI6F9yosMU2tra4Ha7WbunXEIXKQgeS1NFKeZE6BLNCu5wOAAAw8PDGBgYCOqm6g+xFwI1QxwfHxckxBFLLrPZDJ1Ohw0bNrARvYGQkJDgYzM8NTUFvV7P9OJ0dZfrQoqFdlF+zzsdpmAymaDX61l0RW96oUZQC4VUlk2xpkMHYpzgdA/e1tYGp9MZtJHDH5R8Qi8ur9cLh8OByclJwUIcocIVcnYqil6vR2FhoQ8xw50jP4ylenGDwcCMGdVqNYuOYsF22R9SrbZKpZL1vE9MTKCqqgoGgwFdXV2w2+3IzMxkPe+RfA5ShegWi0WSG46UmLMQPRKCezweDA0NoaSkRJDklEKMhp3uiVUqFaqrq0Wp7MKt4PxEXVVVFbRaLfsb1cvTz4WuIKFen5+k8nq9GBgYgMFgwKlTp3zq0ZH2es8HwwdgphRXXFzMet4nJibYjS9S+2kpzvMTvYIL9UmjMBqN6O7uRmZmJkpKSkQ9V+j+2Gw2o7W1FdXV1eju7pa0rm2323Hq1CkUFRVhyZIlMJlMzBedkjshIWEW0emNieO4kGRXKBRITU2F1+tFWVkZHA4HDAYD6/Wmq1p2drZkbiViIRfB+fDXh/Ptpy0Wi0/Pu9wViljrBQfmkOBCk2z8wQDV1dUYHx+P6PXCEXxoaAiDg4NsaGFvb6+opFmomwh1damurmZmC/SGQMnsT2Kv1+vzd2CG7PQxgcjOJ09iYqJPrzdd1Xp7e5GQkOCzugfDfFnBQyFS+2kpEGutosAchehiJKctLS1Qq9VoaGiA1WqNePRwsNejbaQul8tnQkqkk0r8MTo6it7e3oCuLjabDU6nEwkJCbM+J0pgpVLJVnaPx8P80+n7Ebq681c1u93us2cV49MWDc4FwfngBNpPS9UE9InNogvZgweSnEYyhog+L9CX5nQ60djYiNzc3Fl7erEE9388TaZNTEzM6j8nhCA5ORlpaWlsv0yTRsH2y/xVm67uHo+HvSYlvZBtRVJSEpOPUp82g8GA7u5uHxdWqSFHRBANgtlP2+12HDt2jJXiIrWf/kTvwUOF6OPj4+jq6polOY0m+x7Itqm5uRlVVVVsXFC454SCvzKtubkZSUlJqK+v97moKSmVSiUqK2dmttvtduj1enR2drIVNS8vD9nZ2QEvLP7qDoA5yYyPjyMnJwculwuEECiVSkGrO9+nje9YMjExAYVCgcLCwqg91gHp54hJqWKjpbiUlBTodDqsWbMGZrM5Kvtpq9Ua8No6lzinITohBJ2dnZiamgooOZVK4jo6Ooq+vr6QNfRIpac0mUazuhT8ZBq1TaZISkryyQLTGi+9sOjqHqwWT2ejL168mO01+au72+1mEUA4QvBdWJubm5GRkQGdTofOzk4kJyezvXsk/dJSr+By6tBVKhX73AkhrCtOjP30J6JMFgjhJKf+q16o5wl9PRrWnjlzBhaLJWwNXawhg0KhgMPhwPHjx32SaUBocgc6Dr/d0mKxQK/Xo7W1FW63Gzk5OcjNzUVmZiY4jsP09DRaWlpQWVnp04LLX90p2cXu3YGZJFRRUREzdjAYDExNRhNUQm2bpB71O1d2TRzH+dz4PB4PS1qGsp+mDTRicPDgQVx11VWnASgB/IEQ8rjfuewA8HN83EG2ixDyB6HHl4Xg/gkof6kqzTKHk5xG+mUqFAo4nU6cOHECGRkZWLt2bdhjiQ3RtVotpqamsHHjRp+7Nl1N6cUo9j1QM4Vly5bB7XbDYDBgeHgY7e3tUKvVsFqtWLVqVVCLaf+9O/8fILwMxzd2WLp06SzbppSUFHZjCtbsM59W8FCgtfVQ9tNTU1OYmpoStQf3eDy44447gBlbpiEAH3Ec9zohpM3voX8hhHxDzPuimJMVnB/+0vKUUMlpJHC73ejq6sKKFSsEu8gIJTi/ky0jI0NScvtDpVKxpBCV6hYUFKCrqwsAWEgZrCnFn+zAxzkBf5FNuC2Kv20TTVBRU0b+6k7PZaEQ3B+B7KffeustHD58GC0tLbj++uvxjW98A5mZmSGPc/ToUVRUVKC7u7sHAM62hG4F4E/wiDGnpos07BQqOY0EWq0Wo6OjWLJkiSiLKCEEd7vdaG5uRkpKCmpra3Hq1Cn2N0ocIWGwGBBC0NPTg+npaZ+yntPphF6vR29vLywWC7KyspCbmwuNRhNwj0jPyX91p2E8n/jhLnq+VpxGGtRLvKOjg7W/0jq+VIjFVlFqP/3Nb34TH3zwAX7wgx+gs7NTkKhmeHjYf8zREAIbKn6G47iLAZwB8C1CyGCAxwTEnITodrsdVqsVRUVFbFaX1KBlKrPZjJKSEtFfWjiC02TakiVLUFRUBLfbPUuZJsWqzYfX60VraysSExNnTVlVq9U+q4jZbIZer0d3dzfUajXy8vKQm5sbNANMV3eVSsVkr9RHjBI+lMiGD5VK5dMJRttfJyYm0NraitzcXEnEJbHcKgrM5E+WLl2K9evXS3K8s/g7gJcIIQ6O424D8DyAS4U+WfYVnLZfJiYmipacUoT7YunKmpycjPr6eoyNjcHpdIp6jVAEp5JW/nglvv+6HOR2uVxobGxEQUFB2GF2/qUvm80GvV6P9vZ2OJ1OaDQa5nfufzHTLYfT6UR9fT17X/4iG6FlOH776+TkJEpLS2G1WjE0NISpqSmkp6cjNzcX2dnZoqWjsbiC8yE2yVZUVITBQZ/FeJYdEyHEwPvxDwCeEHNOspouUsnpunXrcOLEiYjuwOE6w2j74LJly5hFVCQCmWB7UDpvjEpa+XA6nZiYmJBc9mi1WtHU1ITy8vKI6qrJyclYsmQJcz01Go0YHx/H6dOnkZqayvbuKpUKra2tSEpKQk1NDXsPQkQ2QlZ3QgjUajXS09OZuIS2vw4ODoLjOHbzEdIYEit78GBwuVyi3IUaGhrQ2dkJjuNKMUPs7QA+x38M52uRvAVAu5hzkoXgHo8HjY2NSExMZO2XtOQVqfFioC9Br9fj9OnTbAwwRSR+6v5lMlqjn56enpUzoGF5ZWUlM2bIzMxEXl5e0D2wUJjNZrS3t2PVqlWiSy6BoFQqZyXH6KA9i8WCzMxMLFu2LOjzA4lshEpo/W/M/PZXAD7tr7QxhK7uga6TWF/BxQp7VCoVdu3ahWuuueYNzJTJniGEtHIc90MAxwghrwO4k+O4LQDcmPFP3yHmnGQheHd3N/Lz89mKCiDivuVA001odKDT6QIOLoyE4AqFgnW88ZNp/iU2fvaZZrhpc4dOp0N3dzcSExPZHljIBBeKsbExDAwMYO3ataKeJxQ0OaZSqaDT6VBVVQWlUomBgQFMTU0xguXk5AT9nsKt7m63m4Xy4SI2//ZXurr39/cHbH+N5T14pDLaq6++GoSQKr9jPcz7/+8C+G6k5yULwZcvXy6bL5vH40FrayuUSiXWr18f8MuJ5LXoTcFms+HUqVM+IT8QWrzi39xhtVqh0+nQ2toKj8fDxCr8EhIf9IZlNptRX18vq3kDFcosX76cnS91j5mYmIBer0dfX5+PsiuUXh4ILrJxOp2i2l9pYwiAgO2vtJwnJaSWv8aah96clcmiVaUBM8mjxsZG1mMdDJGu4BaLBSdOnJjlxyY2U56SkoJly5Zh2bJlrINpcHCQrZJ5eXnIyclh7629vR0KhQJ1dXWSX8B8mEwmtqXxF2RwHIesrCxkZWWhoqJCtF4e+Hh193q9aG5uRn5+PuuMo6uvUEVdoPZXuh2anJwU1P4qBGJn3wWD1Ko9qSBbmWzWC0Xo6kITZiaTCW1tbT6Z7FDPEUtwqgc///zzfZJp0ZbB/DuYaCjf29sLlUoFh8OB/Px8lJeXy3qBaLVapscXEv5Hqpen453y8vLYTZjf607/idHL0wjJZrOxaEiq9lepVnCr1RpzOnRgjlfwSJ1VR0dHMTExIXoqqRBQvfrExAQKCgpmkVtKZRp/laRbgczMTExOTuLo0aOzdOdSYXBwEFqtFmvXro1otRKql09OTkZTUxOKi4uxaNEin+cD0evl6fcQqv01MTGRnauQLjCpkmyx2AsOxHiI7vV6YTQa2cwxoV+E0BWcP0ywsrLSxz2GEMJuSFKHzRMTEywaoXtOqgajuvP09HQWyke6J6c1bpvNhrVr10r2PgLp5QcGBqDT6dj7cblcQW8mgRJ1QpxsvN7ZY34DaQAMBgPOnDkDp9Pps7oHev9SruCx1gsOzGGILpbgDoeDldoWL14s6i4rhOD+yTSz2TwrlJQjaaLVatHb24s1a9b4rDD+arDJyUnodDqW8KJZeaFhoNfrRVtbG9RqNVavXi1b+K9SqZCWlobe3l7U19dDqVRCp9Ph5MmT4DhOkF4eCO1kQ0U2Ho8nbJ2Zb8hIJ5/q9Xp0dXWxHm/+XPP4Ci7VC4nYg9NusxUrVmBqaiqihFmo59D9PD+ZJrcyjZCZaat6vR719fUhQ2W+1RBNeOl0Opw+fRoOh8MnlA+0+tDIJCcnJ2SNWwpMTk6yz5LW7TMyMlBeXi5aLw+ELsNZrVYkJSX57N9DgV9qo+2vRqMR7e3tbJiC3W6X5HP4xBNcqVQKko9S5Rj1NLNaraJD+1DRwvDwMAYHB2ft5zmOg8Ph8KnjSgWv14szZ87A4/FEFConJSXNUqaNjo6yxg4ayickJLDIZ+nSpSgsLJTsPQQCNUSoq6sLuN+VQi9P/9vT0wNCCHJzcwFA1N4dmD3XnLa/arVatLa2sgYZjUYTkblFPEQPE6JTEthsNh/lGH+6STSvT5NpVqsV69evn6VMS0xMRHJyMj766COkpqayiy/aEgoVzVDjv2hvHP7KtKmpKeh0OtYs4nA4UFFRITu5aVZ+7dq1gggRqV6eKgpdLpfPViOYhFZoGY5ue4aHh1FdXc1GRrW2tgZtfw0Fi8UiifJQasREko26u2RnZ2P58uU+H2ik9XM+3G43GhsbkZGRgTVr1vgcnybTOI7DihUrWDcU3UdSg8S8vDzRZRC73Y6mpiYsWbLEJ6ssFfjSz9zcXLS0tKCoqAharRYDAwPIzs5GXl5e0ARTpBgZGcHIyEjEWXlAmF4+JycHvb29AICVK1fOEhcBgTPz/r3u/Mf7g/rlJSYmhm1/DTXmOBbtmoA53oMHKpNNTU2hubkZFRUVAd1doiU4bUYpKSmZRTJ+Dze9ePjdUGVlZUzwwd//5uXlhS1lTU1NsTxCuLp9tNDpdOjp6fFpiKHaAT5paCgfzbjlvr4+mEwmrF27VjINt39UMj09Db1ejw8//BAcx6GoqIh1ogX7zIM52YQbKBGozyFY+2tzczOAwN7qkQw9OHjwIO666y6cOXOmC4HtmhIB/AnAOgAGAJ8lhPSJeY1zGqKPjY2hp6dn1jxxPiIRrVDQhIp/M4qYTDlf8OHxeHwslPxVaRQ0a1tbWyt74mVoaAhjY2OzEnd8a2Z+VNLY2OiT3RY63oeW3BwOh6yKO7pX7u3tZSs8Hb4YjV7en/R8ogttf6XTX+kgBXrTodl6MVEatWt66623UF5evhKB7Zq+AsBECKngOG47gJ8B+KzgF8E5CtHpxUKH/YUK8yJdwZ1OJ86cOTMrmRZNGUypVPrc2akqraenhzWYuN1ulimPZqUMB77TS7jV1D8qoTrv7u5uWK1WFsoHk6ESQtDe3g6lUunTVioHqAV1dnY2qwDwJ5VIqZenSVWXyyV47+6vTJyamsLevXvxl7/8BWq1GgMDA/jKV74Sdrw1tWsqKysDIcQZxK5pK4BHzv7/bgC7OI7jiIjOFtkI7u/qQstkfHFJMDdVPsQSnBCC06dPM2so/oUvpTKNr0qrrKyExWJBe3s7LBYLkpOTMTw8jLy8PFED8ISC6teVSuUspxch8Nd5m0wm6HQ6ZqZIE4xqtZrpytPT0yVJEoYCbTPOz8/3saCmkFIv73K50NbWhoqKCp8oUYyTDc2BfPGLX0RbWxs2bdoEu90uSLEp0K6pCMAgABBC3BzHTQDIAaAP+wJnMacruMPhwNGjR1FaWio4nBEyZ4yCJtMyMzNn3dHlVKZ5PB50dXUhKysL69atY6t4T08PLBYLNBqNZMkuj8fDEpJS2F/xZaj8fvHGxkYQQuB0OlFQUCA7ud1uNxvWKPTaCKSXpzeq5OTkoHp5l8uFU6dOoaSkxGeljcbJZnp6GsXFxTj//PMj+wBkwpwR3GAwwGq14vzzz/eZXhIOQldwmkwrLS1FYWEhDAYDkyEGSqZJBYfDgaamJixevJiNXEpISPDpc6b11tOnTyMtLY1deGIz0HT0UlFRkU8rq1Sg/eJpaWkoKirCyZMnkZOTA5vNhg8//NBnhZRyphkl3NKlS0UZZfLhf6OyWq0B9fIpKSkByU2PEa7XPViDjNg6uBC7prM/LwEwxHGcCkAmZpJtgiF7iE4IQW9vLwwGA1JSUkSRGxBmvxQomUafR51a5CB3sCEE/ufPv/D8S3A0exyuMYLaOIV6LalAZ6aXl5czYQkVquh0Oib7pKF8JMIQCofDgVOnTqGsrEyysT98UQtfL0+bbjIzM+H1esPq5QHhDTJi55JRu6be3l6UlZWpEcCuCcDrAG4B8AGAGwC8I2b/Dci8gtOESWJiItatW4cPP/xQ9DHCreCDg4MYHh6elUyj3Wv8HmQpQVVcgXqrgyFYCa6jowNOp5OV4PzFFbQ5paamRvQNUiwsFguam5uxYsUKn554f6EKDeWbm5vh9XrZuYvxp6NOtVVVVbIMP6RQqVTQaDTo7+/H6tWrkZSUxMRBQvTyQPjVnbb/ijmnXbt24YorrgBmfNYC2TX9EcCfOY7rwoxd03ax750Lc0OIeJzjxMQETp48iSVLlrCEyfvvv4+NGzeKPlag53m9Xpw+fRpOpxOrVq2alUxraWmB2+1GYWHhrDJWtBgZGcHw8DBqa2ujWr34oOIKnU6HyclJ5vFGs+XB5KBSYnJyEq2trVi9erWo1YiqwHQ6nWB/Omre4X8jkQNOpxOnTp1CaWnprCiB6uX1er1gvTwfXq8Xu3btwuuvv45Dhw5FKnaRLbkhG8Hb29uh0Wh8vrz3338fF1xwgejV1J/gLpcLTU1NyMrKQllZ2axkGr2zTk5OQq/Xw2AwICkpCfn5+Sw7HAmo97rFYpl1U5EStBzU09MDs9mMrKwsFBQURB0Oh0I4XblQ8P3pjEZjQH86Oip6LiISSu6ysjK23Qh17lQvbzQaw+rlCSH43e9+h3feeQevvvpqNGXR+Udwt9s9K7Q+evRoRJ5jfIJbLBY0NjairKxslt46VDLNYrFAq9VCr9eD4zi29xXTfkmHEFRWVsqaUaYebRMTE1i9ejXrJtPr9azhQsoSHNWV19XVSX4DockunU4Hj8eDtLQ0GI1G1NXVya7ddjqdOHnypE8uQQxsNhv73P318hzH4dlnn8X+/fvx2muvRWuSuTAIfvz4cdTU1Ij+MCjB6RAF/zniYsUrDocDOp0OOp2O7X3z8/OD7h+pVl7IEIJoQQhBR0cHCCFYsWLFrGwttRrW6XSSlOCGh4cxOjqKuro6SbzJQsFoNKK1tRXp6elsaF8gJaAUoOSuqKiQJClJ9fJ6vR7PPfccjh07hsnJSfzjH//A0qVLoz38/CM4FbXwcerUKVRWVoqWb77//vsoLi7G6Ogo1qxZ47PKRGvQQDOsWq0W09PTs1Rd0Q4hEAOalMzIyBBUd+bXfk0m06zW0XCguvLa2lrZthsU1PCRbgH4SkCj0YiEhISw7aNCITW5/fHyyy/j+eefx2WXXYa3334bv/vd77B8+fJoDrkwCN7S0oKlS5eK2nd5vV4cPnwYOTk5AZNpUnqm+RNGrVbDYrFg9erVsjeM0Chh0aJFrJ4uBrQEp9VqYTAYQpbg+LrylStXyurkCoAZJNbV1QWN3mj7qE6ng8vlitifjpbd5CL3nj178Lvf/Q779++XcouxMAje3t6OwsJCwWSh87moQIafxJCa3P4YHR1FX18fNBoNzGYzEhISkJ+fj7y8PMn3qTabjUUJkewVA4Hu2/mEoWWsjo4OKJVKVFVVyW71Sx1k16xZIzgJ5V9REOpPR8ldWVkpS9lt3759eOqpp7B//36pM/8Lg+BnzpxhIXA40GRaeXk5hoaGfPbucirT+EMIVq9ezS4om80GrVYLnU4naaKLtpXyDRilBt2G6HQ6aLVapKSkoKysTJa9Lx9jY2NsFnyk+3u+P53RaGSRib8/ncPhwMmTJ2Wrqb/xxhv42c9+hv3798sRGcw/glOlEB/d3d1ITU0N6zbin0yjd+WUlBRZDRH5QwiWL18eNHSltVOdTgebzQaNRoP8/HzR4aTBYEBnZydqa2tlNwugOn26ilPCSKVI88fIyAhL3kk5qYWKg3Q6HevPz8zMRHd3N5YvXy4Lud9++2386Ec/wv79++XKwywMgvf19SEhISHkHnNgYGBWMq25uRlLly5FWlqabOR2uVxobm5GTk4Oli5dKvj4NLuq1Wp9BCrhhBKjo6MYGhpCXV2drG2lwMe14EA+bVSRJmUJbnBwEDqdDnV1dbJGCB6PB2NjY+js7IRKpWKfvdAkoxC8++67ePDBB7Fv3z45bbAWBsGHhobg8XgCOn16vV50dHTA7XZj1apVPqtna2srCgoK2AopNbnpHri0tDRsH28oEEKYXttoNCI5OZkluugFRwhBf38/TCaTzxZALlA5aEVFRdj9PT8ysVqtEZXg+PV7uZN39L0tX74cWVlZzJ/OYDBAoVCwyCRS041///vf+M53voN9+/bJ0tzDw/wjOG015GN0dBQ2mw1lZWU+v6fdRDk5ObPKQ4QQDA4OYmBgALm5ucjPzxdshCcEgYYQSAF+66VOp2MOK9PT0wCA6upq2QkQTFcuBHTohE6ng9lsDluCo5Jaq9WKmpqaOSN3sPfmcDjYzcput4v2pzt69Cjuvvtu/P3vf5dd+4CFQnCtVouJiQlUVlay301PT6OpqSmgJxu/h5sQwurVU1NTyM7ORn5+flQ91nQIQW1trew6b1pPd7vdSEhIYDerUA0O0YDeuMTqygOB794aqARHnU/dbjeqq6tlz8yHI7c/qD8dvVmF86c7ceIE7rjjDuzduxclJSXSv4HZWBgEp86Z1dXVAMDMDGtra31qiuHEK/wea7PZjIyMDOTn5wtuEOAPIaitrZVdwUXLfYWFhSguLmbNGVqtFhaLRZKbFR9S6cqDwb8ERwhBcnLynITl0Tap8E0d9foZYxRK9tTUVLS0tGDnzp149dVXfRYimTH/CA5glp/5xMQEBgcHUVNTg4GBAYyNjUWtTKOKKK1Wy/a9tKkkEHH5QwjmIkymvdXB9vf+N6toZ5LJqSv3h9frRUtLC5sZNjU1JTjJGAnk6ECjeYempibce++9cLlceOSRR/CFL3xB9s+Ph4VB8OnpaXR1dbEZVP57tWhlp3TfOz4+DoPBwOxvqThF6iEE4UANIYRekPyar8FgYN1MQsU1c6krp+ROS0tjORX/TjJagsvLy4u6UkDJXV1dLYte4PTp09ixYwfuuOMOtLa2YsmSJbj33nslf50gmJ8EdzqdPsaLU1NTOHr0KMrKylBSUhKwzVNKZRpfnOLxeOB0OrF06VLZ53UBH4fJq1evjjiLa7VaWQccIYSRJdDx+OIcuXXlgZxPA4GfZAQQcQlObnJ3dXXhC1/4Av785z+jrq5O8uMLwPwn+PT0NDPyu+iii3xfREZDRGDmxkKH0k9PT8Plcsma5BobG8PAwICkYbLT6WRksdvtPh1w3d3dc6YrD+d8Ggz+4iChWW25yd3f34/t27fjmWeewbp166I+3pe//GXs27cP+fn5aGlpmfV3QgjuuusuHDhwACkpKXjuuedQX18/vwmu0+nQ2dmJVatWobW1FRdccMHMwWUe1Qt8PISAv5K6XK5Z9d5IlGiBwE/eyVXjpgMYaHSSlJTEmivkJHgkzqeB4J/VpiW43Nxcn8+MVh5WrlwpizHE0NAQbrrpJvz2t7/Fhg3+jsWR4d1330VaWhq+9KUvBST4gQMH8PTTT+PAgQM4cuQI7rrrLhw5ckQ2gsvuqtrX1wetVov169dDrVazFX0uyM2f+sHfA/JdT6kSjU4rycrKQn5+fkhv7UCgpSKn04k1a9bISjQ6Fnd0dBTLli1DVlYWM0OUcnAiH1I4n1L4T13hD1CkJbi0tDScPn1aNteX0dFRbN++HU8//bRk5AaAiy++GH19fUH/vnfvXnzpS18Cx3E4//zzYTabwXHcIkLIqGQnwYOsBO/o6IDL5cL69eslTaaFAyWbw+EIO/WDX9Ollj1arRZnzpxBWloay8iHOgZ1e0lKSpJ98gfwsa68oKCAhckajSagaytNMkbjOCKH8ykFf4BieXk5bDYbRkZGcOrUKSQlJUGr1YIQIqm4aXx8HDfddBOefPJJfOpTn5LkmELhP/CguLgYZ86cKQIw/wi+bNkyJCYmzvpi5Gzz9Hg8aGlpQWpqKlatWiXq+HznUH5Gu7e3l3m68WWngG8ThwTOHmERSlfu79pKLYdaW1vh8XhYkktM3mGunE8pvF4vtFotGhoakJyczOyOpSrB6XQ63HjjjXjsscdwySWXSHz2sQdZCZ6amupj20SnRLS1taGgoAAajUZSggcaQhApOI5DZmYmMjMzUVFRwTzd+CtjRkYGTp8+jWXLlkUdtgoBTTgJ9UZPTk7G0qVLsXTpUpZ36O3tZfPIqLgm2HdA98Bz4XwKfGzGuGrVKiZ8onPA+CW47u5upvMXY6JpNBpx44034tFHH8Xll18u51sJCv+BB0NDQ8DsgQeSYc4mm9Ae7vr6ekxMTGB8fBydnZ1IT09Hfn5+1L3JQoYQRIPU1FSUlpaitLQUdrsdg4OD6OzsRFJSEqxWKywWi6yTROnFX11dHRHZ/KetGI1GjI6OoqOjI6A32lw6n/Jfj09uPhQKBbKzs5Gdnc0ml9ARS0D4EpzZbMaNN96I733ve7jmmmtkfz/BsGXLFuzatQvbt2/HkSNHkJmZCbn234DMWXS32w232x10v03DYCpMSU1NZXteMRnoSIYQRAOz2YyOjg6sWrUKarUaer0eWq0Wdrudld/EDAAIByl15f7ge6NRe+n09HSMj4+HHOssJSi5I31/4Upwk5OTuOGGG3DXXXfhxhtvlOEdfIybb74Zhw8fhl6vR0FBAR599FHWVfm1r30NhBB84xvfwMGDB5GSkoJnn30W69evn59lMrfbDafTKSiZRhNE4+Pj0Ov1Qfe8/pBjCEEo0AaVQP5i/gaOtN0yOzs7YrJTU4i5GHwAzNTwT58+jaSkJJ8EpFyGFNPT02hubpbs5sUvwfX29uIPf/gDDAYDbrvtNtx6660SnLEsmJ8E/9rXvoaCggJs27ZNtP8XNRDU6/VMcpqfn8/2W3M1hIAPOttKSIMKDYNpB53YhhhgJtvb398/J7pyYLbzaSB76UCjlSKF1OQOdPwvfOELSElJwdjYGK644go8+uijkr+OBJifBNfpdNi7dy/27NkDrVaLK6+8Etu2bRPdUkglmzqdDgqFArm5uTCZTEhJSZF9CAHw8c3EarXOMqMQ+nx+Q0xKSkrYrcjQ0BALk+XWlQPhnU/5vm5TU1MR6wUoKLlra2tlyV3YbDbcfPPNuOGGG7Bz504AM0nYOWwgEQOfC5jjuH8COEQI+UnUB5aT4HyYTCa8/vrr2LNnDwYGBrB582Zcf/31olsMp6am0NjYCI7joFar2couV/jq9XrR1taGhIQESVxI+fbGer2eNZTwo5O51JUDH3egCXU+5esF+H7sQnMn1GhSLnI7HA58/vOfxzXXXIPbb79dkgXg4MGDuOuuu+DxeHDrrbfigQce8Pn7wMAAbrnlFpjNZng8Hjz++OO4+uqrhR6enSDHcQ8D2AZADeAFQsjj0Zz3nBGcj8nJSezbtw979uxBZ2cnLrvsMmzbtg319fUhye4/hMDpdEKr1UKr1cLtdjOiSHXRuN1uNDU1IScnR7YGFX50Qi9ElUqF2tpa2XXlQPTOp/43LDrAIJi4Rm5yO51O3HLLLdi0aRPuvvtuScjt8XhQVVWFt956C8XFxWhoaMBLL72ElStXssfs3LkTa9euxde//nW0tbXh6quvDqlo8wMHABzHKQBcA+AwAA2AAwD+Sgh5lOM4TuzoYGAOy2R8ZGRk4HOf+xw+97nPwWq14sCBA/if//kftLS0YNOmTdi6dSs2bNjgs3qZzWY2A5yWUdRqNYqLi5mJgk6nw5kzZ+BwOBjZI20mcTgcaGxsDCgokRIpKSkoKSnB0qVL0dLSwpKSx44dYxl5qWaQ+WN4eBhjY2NYu3ZtxLp5vriGKtH44hp+BxwtZcpFbpfLhVtvvRUbN26UjNzAjH1TRUUFa4vdvn079u7d60NwjuMwOTkJYKbqEYmHGyHEy3HcPgCJhJB+juOuAfB3juPUhJDvcxyXCKCMENIu9JjnZAUPBrvdjjfffBOvvPIKTp48iYsuugjbtm1DZ2cncnJycOWVVwqSXLrdbuj1eoyPj8Nms7HOK6HJIeplNlfqLaq+y8jIYG20VJii1WrZe8jLy5OkIQaYSRjSphi5tgH0pkvHCrtcLqxYsQIFBQWS37Dcbjduu+02VFdX46GHHpL0+Lt378bBgwfxhz/8AQDw5z//GUeOHMGuXbvYY0ZHR3H55ZfDZDLBYrHg0KFDgrvTOI77fwD2EUJcHMcpzhJdSQjxcBxXAuAvAP4NYD1mwvbfCT33c7KCB0NSUhK2bNmCLVu2wOl04tChQ3jwwQeh0+lwySWXID09HRdffHHYUFKlUqGwsBCFhYWs82pgYICVrkIpuGjNOZjgQmoE0pUDsxtiDAZD1A0xFNT5tK6uTtZtQEJCAhYvXoy0tDS0tLSgoqICRqMRvb29EVUVgsHj8eCb3/wmysvLJSe3ULz00kvYsWMH7rnnHnzwwQf44he/iJaWlrCf7y233AIAuwBcwXHctwghtrMk93AcpyKE9HEcdwOAfgBPiSE3EGME50OtVsPhcOC8887Dk08+iX//+9/YvXs3HnjgAaxbtw7btm3DJZdcEjYrSmWl+fn58Hq9MBgMGBkZQUdHxyyi6HQ69PT0YM2aNdEn7VQApwBAAOIK/BCqKw8ndfV/D/yGGDFKQL7z6Vz4pwEz+Za2tjasWbOG1dL5VYVIZacUXq8X3/rWt5Cfn48f/vCHspA7kLzUXwr9xz/+EQcPHgQAXHDBBWxAQygb7ra2NnrcLwC4GsB/nSX59NkV3M1xnArAHQB+Twj5lthzj6kQfdaLnz03/pfm8Xjw3nvv4dVXX8U///lPrF69Gtu2bcOnP/1pUaSkRBkfH4fZbIZSqYTH40F9fX3UpRQuAQAHcBxAP17i6z8pWlceCFQJSAcOBvJh5z92Lp1PgY/JHWpyC99eWq/XMz/zQEMT/eH1enH//fdDpVLh17/+tWw3LLfbjaqqKrz99tsoKipCQ0MDXnzxRdTU1LDHXHXVVfjsZz+LHTt2oL29HZdddhmGh4fDfs4jIyMoKipKBlAC4GsAMgDcRwgx8ML1lYSQtkjOPaYJHg4ejwcffvghXn31VRw6dAhVVVW4/vrrsXnzZsHCCVrjNplMSE9PZ2WfgoKCiPXxnHqG3B+/BkA8AM723dAasNRe7LRVlPqw8/3oOjo6oFAo5mTgIDCz1WlvbxetwKMrn1arZc47dNwS/7y9Xi8efPBB2O12/OY3v5E9Gjlw4ADuvvtueDwefPnLX8b3v/99PPzww1i/fj22bNmCtrY2fPWrX8X09DQ4jsMTTzwRsqGFdlOeBcfN/FAB4HYAqQDuA/CfAF4ihIxHet7zmuB8eL1enDhxAq+88goOHjyI0tJSbNmyBVdffXXQZgk6TYXOIuM4LuCqWFBQIEofH5DgXgBueXXlfNBsNpXNpqWlobq6WtaGGIpIye0PmiylSbrs7Gz2fTz++OPQ6/X4/e9/PydaATnAIzktk6kA5AC4G8C3AfydEHJDNK+xYAjOh9frRXNzM1555RX84x//QGFhIbZs2YJrr72WjS72eDxoampCVlbWLANIikCilIKCgrD6+GAh+lzryqnzaUpKCpKTk6HVapnkVOqGGAqpyO0Pai/9yiuv4Ne//jVUKhUTk8xFQ4zU8Hg8UCqVsNvtSE5OTiCEsFG8HMe9B2CAEPL5aF9nQRKcD0II2tvbsXv3buzbtw9ZWVm47LLL8MYbb2DXrl2iJlfQnnCdTjfLknkW+Ek2NzA+NqMrFzMnOxrQG5hGo/ER6QRqiAnXFy4UtMtOrhsYIQRPPvkkWltb8e1vfxv79u3DZZddhk2bNkn+WnKCkttms+HKK6/Eu+++ezMh5OWzYfpyAN8ihNwmxWsteILzQQjBO++8gy996UsoLS2FWq1mZTmxtVlqyazVasFxHMtyB6rTz7Wu3OPx4NSpU7NKb/7wb4jJzMxkpSuxe1pK7jVr1kRlDxUMhBA8/fTT+Oijj/Dyyy9L8jmGk58CwF//+lc88sgj4DgOdXV1ePHFF6N6TUpup9OJ6667Dg0NDfjxj388P5tNYhFPPPEEPv3pT2Pt2rXo6+vDq6++itdeew1KpRLXXXcdtm3bhsWLF4siOx3lo9Vq4fV6mYouOTnZZ9rmXOwV6Zgksc6ndDIqbYgR05s/F+T+7W9/i8OHD2P37t2SREBC5KednZ246aab8M477yA7OxtarTaq6bNerxcKhQJerxfXXXcdampq8MQTTwAzSTYFIcQb9RvzwyeO4IFACMHw8DAjO727bt26FcuWLRNFdupfPj4+jqmpKajV6jkTzQitq4cDvzefTlih2xF/ctEWUznJ/cwzz+DAgQN47bXXJHuNDz74AI888gjeeOMNAMBjjz0GAPjud7/LHnP//fejqqoq4j5ySmj/323evBm1tbX41a9+BQA4W/OWnNxADAtd5hIcx6G4uBh33XUX7rzzToyPj2PPnj248847MTk5iWuuuQbbtm1DRUVFWLKr1WosWrQIZrMZycnJyMzMRHd3N3N7KSgokGXYAnU+LS8vDzsHPBz4+vKKigrWEEO7+CjZ7Xa7rOQGZmShr7/+Ol5//XVJXyOQu+mRI0d8HnPmzBkAwIUXXgiPx4NHHnkEV155paDj88k9PDwMQgiKi4vR2NiIK6+8Evfddx97nFzkBuIEnwWO41BYWIjbb78dt99+O3Q6Hf72t7/hgQcegE6nw1VXXYWtW7cGFYvQsT6ZmZksO7948WKW3KKmhxqNBgUFBZKYJ8jtfEobYkpKSth2pKmpCdPT01iyZAmbSiM1Xn75ZfzlL3/Bvn375qTq4A+3243Ozk4cPnwYQ0NDuPjii9Hc3CzIE4+S+7777sPg4CCsVivWrFmDH/7wh1i7di2Aj/fjciJO8DDIy8vDV7/6VXz1q1+FyWTC3r178cMf/hCDg4O4/PLLcf311zMTCJfLhaampoDJLZVKxRxCqbac2gFTh9NIrJ3m2vk0KSkJqampIIRgw4YNmJycRHd3d0RNPaGwZ88ePP/889i3b58stXsh8tPi4mJs2LABCQkJKC0tRVVVFTo7O9HQ0BDy2JS4Tz/9NPr7+/HCCy9g+/bt0Gq1Po+bi5xMfA8eISYmJlhPe1dXFy688EK8//77ePbZZ7F8+XLBxwmUyS4oKBDUSDLXzqfAjMFlZ2fnrLHP/HFK/JtWJDPP//73v+Ppp59mZU05IER+evDgQbz00kt4/vnnodfrsXbtWpw6dSqotLinp4e1lAIzEUh1dTVefPFFDA4O4sUXX4TVasXJkydx4YUX8p8az6LHMtra2nDttddi+fLlGBoawiWXXIKtW7fivPPOE3WXJoSwWeFUOkv94/2PQ40T5FbE8UFtnfzJ7Y9AM8+FNsQcPHgQP//5z7F//37ZW3XDyU8JIbjnnntw8OBBKJVKfP/738f27dsDHuvVV1/FLbfcgttuuw1PPvkkAODpp5/Gd7/7XWzbtg0vvPACAOAHP/gBpqam8Itf/IJ/44sTPJbx17/+FSUlJTjvvPNgs9nw5ptvYvfu3Th58iQ+9alPYevWrdi4caMoUwV+xxW1lKaS2enpabS3t4ds4pAalNxr164VVaYKJP2l5Tf/Wvbbb7+NH/3oRzhw4EDUicK5hNFoxLZt27B9+3b84x//QFlZGZ566ikAwF133YV3330XTz31FPbv34/33nsP+/fvZ4rKs4gTfD7C4XDg7bffxu7du3HkyBFccMEF2LZtGz71qU+JEmrQAX1arRbj4+NwOp0oKyvD4sWL59SQUSy5/UE7x/huuZTIfX19ePDBB7F///45mRIjNU6fPo2SkhKMjo7itttuQ2lpKf73f/8XAPCTn/wEk5OTsFgseOSRRwLdvBYGwY1GIz772c+ir68PJSUl+Otf/+p/JwMAPP/88/jxj38MAHjwwQdpUzzDli1b0NPTE3A8a6zC5XLh8OHDePXVV/Hee+9h/fr12LZtGzZt2iS4PZUSrbKyEmazGTqdLmSNWgro9Xp0d3dHTe5AsNls6OjowNe+9jWMjo7i9ttvx44dO3z2sfMRQ0ND+OpXv4ri4mL8/ve/R2NjI7KyskL5+i0Mgt9///3QaDR44IEH8Pjjj8NkMuFnP/uZz2OMRiPWr1+PY8eOgeM4rFu3DsePH2c3gj179mD37t1oamqaVwTnw+1241//+hd2796Nw4cPo7a2Ftu2bcNll10WtBwUzPmU1qi1Wq2PMYQU9sB6vZ4ZYMilnz9y5Ai+/e1v45lnnsHx48dhNBpx//33y/Jac4mRkRHcc889OHbsGBISEvDee++F6vuXr3+XEBLqn6SoqqoiIyMjhBBCRkZGSFVV1azHvPjii2Tnzp3s5507d5IXX3yREELI1NQUufDCC0lrayupqamR+vTOCdxuN3nvvffI3XffTVatWkVuuOEG8sILLxCtVkssFguxWCykq6uLHD58mJjNZva7QP8MBgNpb28n//d//0cOHz5M2traiF6vD/mcYP/6+vrIP//5T2IymSJ6vpB/7733HqmtrSW9vb2SfZ7/+Mc/SFVVFSkvLyePPfZY0Mft3r2bACAfffSRZK/tj8cee4wUFRWRrq6ucA8Nx8OI/81pHXx8fJzpowsLCzE+PruPPZDCaHh4ZvjiQw89hHvuuWfOEktzAaVSiYsuuggXXXQRvF4vjh8/jldeeQVPPPEEm8dNznZRhUvSJSUlsWmiDocDWq0WbW1tzN1UqKU0Hfsj58rd2NiI22+/HXv27BHV0RcKHo8Hd9xxh4++fMuWLT76cmCmAvHUU09hw4YNkrxuIBgMBhw6dAiHDh1CeXm5bK8TDpLbYHz605/GqlWrZv3bu3evz+PEzgY/deoUuru7cf3110t9yjEDhUKBhoYGPPHEEzh58iRWrFiBw4cP4+TJk/jc5z6HF154ASaTSdCxEhMTsWTJEqxbt44R9cyZMzhy5Ai6u7sxPT3NLLH4oOSWY89N0draittuuw2vvPIKKioqJDsu395YrVYze2N/PPTQQ/jOd74jm7wWAHJycrBv3z6sWLFCttcQAslX8EOHDgX9W0FBAUZHR7Fo0SKMjo4G7MwpKirC4cOH2c9DQ0PYtGkTPvjgAxw7dgwlJSVwu93QarXYtGmTz2MXElwuF9xuN5qbm6FWq9HW1obdu3dj27Zt0Gg02Lp1K6699lpB5SS1Wo2ioiIUFRXB7XazGdtUfVZQUID09HTodDr09fVh7dq1smXnOzo68JWvfAUvvfSSKEGQEAjRl584cQKDg4O45ppr8POf/1zS1/eHnDcQoZjTEH3Lli14/vnn8cADD+D555/H1q1bZz3miiuuwPe+9z22Ur355pt47LHHoNFo8PWvfx3ATEnl2muvXbDkBmZWYP4FWFNTg5qaGjz88MPo7OzE7t27cdNNNyElJQVbt27FddddJ6inXaVS+dgx6/V69Pf3w2w2gxCClStXRjwEIRw6OzuxY8cO/PnPf/ZRjM0VvF4vvv3tb+O5556b89c+V5DfN5eHBx54AG+99RYqKytx6NAh1mB/7Ngx1pKn0Wjw0EMPoaGhAQ0NDXj44YeDKpqMRiM2b96MyspKbN68OWj4+vzzz6OyshKVlZV4/vnnAcxkn6+55hqsWLECNTU1AZv9YxEcx6Gqqgrf+9738MEHH+APf/gDnE4nvvjFL+Kqq67Cb37zG4yMjAQMv/2hVCqZPj4xMRFVVVUYHx/Hhx9+iI6ODhiNRni90jQ69fX14Utf+hKeffZZ1NXVSXJMf4TTl1P136ZNm1BSUoIPP/wQW7ZswbFjx2Q5n5hAmCxcTOO+++5jmdLHHnuM3H///bMeYzAYSGlpKTEYDMRoNJLS0lJiNBqJxWIh77zzDiGEEIfDQS666CJy4MCBOT1/KeH1esnAwAD51a9+RS6++GKyceNG8thjj5HW1lYyPT0dNJPd09MzK0M/NTVFBgYGyLFjx8ihQ4fI0aNHSX9/P5mamoooW97R0UHq6urIhx9+KOtn4HK5SGlpKenp6SEOh4PU1taSlpaWoI//j//4D1mz6CIgWxZ9TldwqbF3714mgrnlllvwt7/9bdZj3njjDWzevBkajQbZ2dnYvHkzDh48iJSUFFxyySUAZvao9fX1GBoamsvTlxQcx2HJkiW4++67mfNJeno6vvnNb+LSSy/FL37xC3R2dvqs7OPj4xgYGJi151YoFMjJyUF1dTXOP/98LF68GAaDAUeOHEFLSwu0Wi08Ho+g8xoZGcH27duxa9cuWbPWwMz2Y9euXbjiiitQXV2Nm266iW1rXn/9dVlfO1Yxr6WqWVlZMJvNAGYikezsbPYzxS9+8QvY7XY8+OCDAIAf/ehHSE5Oxr333sseYzabUV9fj0OHDs17FVUg6HQ6vPbaa9izZw8MBgOuuuoqqFQquN1u3HvvvYL33ISnK9fr9WFtncbGxnDjjTfiySefnHfGiHMM2YQuMd8P/ulPfxpjY2Ozfv+Tn/jORhdbdqNwu924+eabceeddy5IcgMzPe07d+7Ezp07YTQa8f3vfx979+5FUVERnE4nrr/+etTU1IRt6+Q4DpmZmcjMzERFRQWzderr60NSUhKTzCYkJECn0+Gmm27C448/Hif3OUTME1yushvFzp07UVlZibvvvlvCs45dKJVKjI2N4fTp0/B6vdi3bx9+9rOfobu7G5s3b8a2bduwZs0aQWTn2zrRJpL9+/fj17/+NaxWK77//e9j8+bNc/TO4giEeR2i33fffcjJyWHadqPRSF0qGYxGI9atW4cTJ04AAOrr63H8+HFoNBo8+OCDaG9vxyuvvDIng/hiGdPT0zhw4AB2796N9vZ2XHrppdi6dSsaGhpE9bSbzWZ89rOfRUVFBbq6utDQ0IBf/vKXMp75gsDC0KJLDb1eTy699FJSUVFBLrvsMmIwGAghhHz00UfkK1/5CnvcH//4R1JeXk7Ky8vJM888QwghZHBwkAAgK1asIHV1daS0tJQUFBQE1TDb7XZy0003kfLycnLeeef56Kd/+tOfkvLyclJVVUUOHjwo75ueA1itVvLaa6+Rz3/+86SmpoZ8/etfJwcPHiQTExMhs+Wjo6Nk48aN5JVXXmHHcjgcUZ9POH35k08+Saqrq8nq1avJpZdeSvr6+qJ+zTmGbFn0eU1wqeB2u0lZWRnp7u5m5ZXW1lafx/z3f/83ue222wghhLz00kvkpptuIoQQ0traSmpra4ndbic9PT2krKyMuN3uOX8PcsFut5N9+/aRHTt2kJqaGnLrrbeSffv2zWp8GR8fJxdffDFrDJIKQr6bd955h1gsFkIIIb/5zW/YdzOPECe4nHj//ffJ5Zdfzn7+6U9/Sn7605/6PObyyy8n77//PiFkpt6ak5NDvF7vrMfyH7fQ4HQ6yRtvvEF27txJVq5cSW655Rby2muvkZGREXLJJZeQ5557TvLXFPLd8HHixAmyceNGyc9DZsTr4HIiVAdboMeoVCpkZmbCYDAIeu5CQUJCAi6//HL89re/RWNjI3bs2IG33noLq1evxuWXXz7LmEMKiP18//jHP+Kqq66S/DzmK2I+ix5HbEKlUmHTpk3YtGkTfvnLX87JQMVweOGFF3Ds2DH83//937k+lZhBfAWHMI9s/mPcbjcmJiaQk5Mj6LkLHYmJiZJPaqEQ+vkeOnQIP/nJT/D6669L4mazYBAmhv9EQIiGedeuXT5JthtvvJEQQkhLS4tPkq20tHRBJdnONYR8NydOnCBlZWXkzJkz5+gso0Y8ySY39u/fTyorK0lZWRn58Y9/TAgh5KGHHiJ79+4lhBBis9nIDTfcQMrLy0lDQwPp7u5mz/3xj39MysrKSFVV1ayGlXAlnmDltzfffJPU19eTVatWkfr6evL222/L9M5jH+G+m8suu4zk5+eTuro6UldXR6677rpzebqRIE7w+Yhoym8nTpwgw8PDhBBCmpubyeLFi+f25OOYS8Sz6PMRQiyE+B1xN9xwA95++20QQrB27VosXrwYwIzZg81mg8PhmPP3EMf8RpzgMiKa8hsfr776Kurr6+PJozhEI14mi3G0trbiO9/5Dt58881zfSpxzEPEV3AZEU35jT7++uuvx5/+9Kdzar0bx/xFnOAyoqGhAZ2dnejt7YXT6cTLL7+MLVu2+DyGGlECwO7du3HppZeC4ziYzWZcc801ePzxx/1Hzc57HDx4EMuXL0dFRQUef/zxWX93OBysI23Dhg3o6+ub+5NcKAiThYsjSkRafvvRj35EUlJSWOmnrq6OjI+Ps+NGWn6j6O/vJ6mpqeTnP/+5TO88MKKpLCxgxMtkcXwMKUjymc98htxwww1zTvBoGnsWMOJlsjg+RjTlNwD429/+htLS0nPiTS5VZSEOYQjn6BJHDILjuBsAXEkIufXsz18EsIEQ8g3eY1rOPmbo7M/dADYAsAN4C8BmAPcCmCaE/GI+nDshRD9X57lQEF/BP3l4BMCvCCHT5+j1hwEs4f1cfPZ3AR/DcZwKQCaA+BIeAeJ18PkJMSQZ8iPJBgA3cBz3BIAsAF6O4+yEkF2yn/UMPgJQyXFc6dlz3A7gc36PeR3ALQA+AHADgHdIPNSMCHGCz09EQ5JP0QdwHPcIZkL0uSI3CCFujuO+AeANAEoAzxBCWjmO+yGAY4SQ1wH8EcCfOY7rAmDEzPuLIwLE9+DzFBzHXQ3g1/iYJD/hk4TjuCQAfwawFmdJQgjp8TvGI5jjPXgcc4s4weOIYwEjnmSLI44FjDjB44hjASNO8DjiWMCIEzyOOBYw4gSPI44FjDjB44hjASNO8DjiWMCIEzyOOBYw/n/8h7FtidqN0AAAAABJRU5ErkJggg==\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"done!\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"class Ion:\n",
|
|
" def __init__(self,position, charge, mass, v_z):\n",
|
|
" self.position = np.array(position)\n",
|
|
" self.charge = charge\n",
|
|
" self.m = mass\n",
|
|
" self.v_z = v_z\n",
|
|
" self.F = np.array([[0],[0],[0]])\n",
|
|
"\n",
|
|
"class infiniteLengthQuadrupole:\n",
|
|
" def __init__(self,R,r_0):\n",
|
|
" # rod radius\n",
|
|
" self.R = R\n",
|
|
" # rod placement radius\n",
|
|
" self.r_0 = r_0\n",
|
|
" # rod positions\n",
|
|
" self.pPole = np.array([[[ (R + r_0)/np.sqrt(2)],\n",
|
|
" [(R + r_0)/np.sqrt(2)],\n",
|
|
" [0]],\n",
|
|
" [[-(R + r_0)/np.sqrt(2)],\n",
|
|
" [(R + r_0)/np.sqrt(2)],\n",
|
|
" [0]],\n",
|
|
" [[-(R + r_0)/np.sqrt(2)],\n",
|
|
" [-(R + r_0)/np.sqrt(2)],\n",
|
|
" [0]],\n",
|
|
" [[(R + r_0)/np.sqrt(2)],\n",
|
|
" [-(R + r_0)/np.sqrt(2)],\n",
|
|
" [0]]])\n",
|
|
" # rod pseudo charges\n",
|
|
" self.rodsPseudoQ = np.zeros(4)\n",
|
|
" # K factor for simplification ofcalculating the equations\n",
|
|
" self.K = 1 / (4 * np.pi * sc.epsilon_0)\n",
|
|
" \n",
|
|
" def check(self):\n",
|
|
" print('pole positions')\n",
|
|
" print(self.pPole)\n",
|
|
" print('pseudo charges')\n",
|
|
" print(self.rodsPseudoQ)\n",
|
|
" \n",
|
|
" def phi_0(self,U, V, f, t):\n",
|
|
" return U + V * np.sin(2 * np.pi * f * t)\n",
|
|
" \n",
|
|
" def coulombForceOnQ1(self,Q1,Q2,r1,r2):\n",
|
|
" a = self.K * Q1 * Q2\n",
|
|
" r12 = r1 - r2\n",
|
|
" mag_r12 = np.linalg.norm(r12)\n",
|
|
" # calculate force\n",
|
|
" return a * r12 / (mag_r12**3)\n",
|
|
" \n",
|
|
" \n",
|
|
" def calcNewIonPos(self,U,V,f,ion,t_s,t):\n",
|
|
" \n",
|
|
" # check if ion has left the r_0 boundary\n",
|
|
" if np.linalg.norm(ion.position) >= np.linalg.norm(self.r_0):\n",
|
|
" return None\n",
|
|
" \n",
|
|
" signLUT = [1,-1,1,-1]\n",
|
|
" index = 0;\n",
|
|
" \n",
|
|
" # reset force of ion\n",
|
|
" ion.F = np.array([[0],[0],[0]])\n",
|
|
" \n",
|
|
" for Q in self.rodsPseudoQ:\n",
|
|
"\n",
|
|
" \n",
|
|
" # calculate the pseudo charge points of the quadrupole\n",
|
|
" Q = signLUT[index] * self.phi_0(U,V,f,t)/2 * self.K * self.R\n",
|
|
" \n",
|
|
" # match pseudo charge with ion in z direction\n",
|
|
" self.pPole[index][2] = ion.position[2]\n",
|
|
" \n",
|
|
" # calculate force on ion\n",
|
|
" ion.F = ion.F + self.coulombForceOnQ1(ion.charge,\n",
|
|
" Q,\n",
|
|
" self.pPole[index],\n",
|
|
" ion.position)\n",
|
|
" \n",
|
|
" # update index\n",
|
|
" index = index + 1\n",
|
|
"\n",
|
|
" newPosition = t_s**2 * 1 / ion.m * ion.F + ion.position\n",
|
|
" # velocity component in z\n",
|
|
" newPosition[2] = newPosition[2] + ion.v_z * t_s\n",
|
|
"\n",
|
|
" return newPosition\n",
|
|
" \n",
|
|
"t_step = 1e-9\n",
|
|
"t_sim = 0\n",
|
|
"N = 1000\n",
|
|
"\n",
|
|
"# coffein atom 524.50002164 dalton (8.709527e-25 kg)\n",
|
|
"ion = Ion(np.array([[0],[0],[0]]), 1000 * sc.e, 8.709527e-25,0.01)\n",
|
|
"\n",
|
|
"# quadrupole\n",
|
|
"quad = infiniteLengthQuadrupole(5e-3,15e-3)\n",
|
|
"\n",
|
|
"quad.check()\n",
|
|
"\n",
|
|
"positions = []\n",
|
|
"\n",
|
|
"positions.append(ion.position)\n",
|
|
"\n",
|
|
"for n in range(0,N,1):\n",
|
|
" \n",
|
|
" newPos = quad.calcNewIonPos(100,500000,1.1e6,ion,t_step,t_sim)\n",
|
|
" \n",
|
|
" if(newPos is None):\n",
|
|
" print('ion has left the boundary')\n",
|
|
" exit\n",
|
|
" \n",
|
|
" ion.position = newPos\n",
|
|
" positions.append(newPos)\n",
|
|
" t_sim += t_step\n",
|
|
"\n",
|
|
"\n",
|
|
"print('simulation time = ', t_sim)\n",
|
|
"ax = plt.axes(projection='3d')\n",
|
|
"ax.scatter3D(positions[:][0], positions[:][1], positions[:][2], c=positions[:][2], cmap='Greens')\n",
|
|
"plt.show()\n",
|
|
"\n",
|
|
"print('done!')\n",
|
|
"# dummy check\n",
|
|
"if 0:\n",
|
|
" ion_1 = Ion(np.array([[0.0000005],[0],[0]]), 1000 * sc.e, 1,0)\n",
|
|
" ion_2 = Ion(np.array([[-0.0000005],[0],[0]]), 1000 * sc.e, 1,0)\n",
|
|
" print(quad.coulombForceOnQ1(ion_1.charge,ion_2.charge,ion_1.position,ion_2.position))\n",
|
|
" "
|
|
]
|
|
}
|
|
],
|
|
"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": 2
|
|
}
|