{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1586 particles were reflected out of 1000000 (0.159%).\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAF9CAYAAAAZcuPcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd5gUxdPHv317gQvAcXfknEVJIkGUoCCIBMmIKBhQUFExoGLkpwgGVFAUFBUFEQQEJCkKIkmRKBIkipJzhiNcqPePct8NM7NxNtfnefq53Z6e7t7bmdqe6gqKiCAIgiBELnGhnoAgCILgHyLIBUEQIhwR5IIgCBGOCHJBEIQIRwS5IAhChCOCXBAEIcIJmSBXSo1XSh1VSm02qb8FSqnTSql5TvVfK6W2K6U2/zdmghnjCYIghAuhXJF/CaCNif2NANBbp/5rAFcBqAUgGcADJo4pCIIQckImyIloGYCT9nVKqcr/razXKaWWK6Wu8qK/nwGc06n/nv4DwGoAZfyduyAIQjgRbjrycQAeI6LrAAwCMMasjv9TqfQGsMCsPgVBEMKB+FBPwIpSKg3ADQCmK6Ws1Un/HesC4DWd0w4Q0a0eDjEGwDIiWu7vXAVBEMKJsBHk4KeD00RU1/kAEc0EMNPXjpVSQwAUBdDf9+kJgiCEJ2GjWiGiswD+UUp1BwDF1PG3X6XUAwBuBXAnEeX7258gCEK4oUIV/VApNQXATQCyABwBMATAYgBjAZQEkADgGyLSU6no9bccbJ2SBuAEgL5E9KNSKhfAHtg2Qmd62qcgCEIkEDJBLgiCIJhD2KhWBEEQBN8QQS4IghDhhMRqJSsriypUqBCKoQVBECKWdevWHSeios71IRHkFSpUwNq1a0MxtCAIQsSilNqjVy+qFUEQhAjHb0GulCqglFqtlPpTKbVFKfWqGRMTBEEQPMMM1cplAC2I6Px/8UxWKKV+IKLfTehbEARBcIPfgvy/qILn/3ub8F8R43RBEIQgYYqOXCllUUptAHAUwEIiWqXTpp9Saq1Sau2xY8fMGFYQBEGASYKciPL+C3ZVBkBDpVRNnTbjiKg+EdUvWlRjPSMIgiD4iKlWK0R0GsASmJv5RxAEQXCBGVYrRZVS6f+9TgZwC4Bt/vYrCIIgeIYZVislAUxQSlnAPwzTiGiem3MEQRAEkzDDamUjgGtNmIsgeM+ePcCUKQAR0LMnULFiqGckCEEnnDIECYJ3/PIL0LYtcOkSv3/1VWDuXKBVq9DOSxCCjLjoC5HLs8/ahDgAXL4MPPNM6OYjCCFCBLkQuWza5FmdIEQ5IsiFyKVBA21d/frBn4cghBgR5ELk8s47QKFCtvcFCwLvvRe6+QhCiJDNTiFyadQI2L0bmDGDrVa6dgWyskI9K0EIOiLIhcgmMxPo1y/UsxCEkCKqFUEQhAhHBLkgCEKEI4JcEAQhwhFBLgiCEOGIIBcEQYhwRJALgiBEOCLIBUEQIhwR5IIgCBGOCHJBEIQIRwS5IAhChCOCXBAEIcIRQS4IghDhiCAXBEGIcESQC4IgRDgSxlYIDw4dAsaNA/bvB1q3Brp1A5QK/jxycoApU4DffgOuugq4/37H5BWCEIaIIBdCz969nLbt6FF+/9lnQP/+wMcfA/n5QFwQHxw7dwbmz7e9HzcOWLWKsw8JQpgiqhUh9IwaZRPiVsaNA5o3B+LjgVKlgNGjAz+P5csdhTgAbN0KTJgQ+LEFwQ9EkAuhZ+dObR0RsGwZ/z10CHj8cU7pFkh27NCv3749sOMKgp+IIBdCT/PmnrWbODGw82jSRL++adPAjisIfiKCXAg9jzwCNGtme2+kE7dYAjuP6tWB//3PcZO1e3dO6iwIYYxsdgqhJyUFSEiwvc/PZ2FK5NiuVKnAz2XIEKBXL2DlSqBGDd6EFYQwRwS5EBz27eOVdunS2mOrVgE//+xY5yzEAWD9+sDMzZmqVbkYcfQo8PLLwC+/ABUrAi+9ZI765coVYNIkNn2sUQN44AGgcGHj9jt2sKmkxQLcdRfPRYhJRJALgeXIEaBnT2DJEn7fpg0weTJQpIitzcGDnvUVTDNEI4iAW28FNmzg9zt38mdbuxaoVcu/fjt0AH76yVb32WfA6tX6po9z57LKJyeH3w8bBvzwA3DTTb7PQYhYwuDOEKKahx+2CXEAWLAAeOopxzbNmwPJye77uv9+U6fmkn/+AX7/3SYoraxYYRPiVq5cYXNJf1iyxFGIA8C2bcYbvM8+6zi3S5eA55/3bw5CxCKCXAgc+fnAnDna+u++c3yfkQF8+aXNg9LIo7NSJeOxdu4EPvmEfyjy832aLgDg8mWgRw8eq3FjoFw5ti+3cvas/nnnzvk+JqBvgmlUn5fHQt6ZzZv9m4MQsYggF9xz7hwwdCi7zj/6KK9WPSEuDkhP19ZnZGjrunQB1q0DFi9myxE9FizQrx85ki1OHnoIuO02NiO8cMGzOTrzwQfA9Om294cPA3feCeTm8vsWLYDMTO153bv7Np4VIx27vTWPFYsFqFdPW9+okX9zECIXIgp6ue6660iIEPLziRo3JmItLpeiRYkOHPDs/KFDHc8FiD74wLHNV18RlSjBx8qVIxo8WHsOQPThh9r+Dx8mSkjQtn3nHd8+70036Y/9++9ES5YQrVtH9OuvRNWrc32hQkRvvOHbWM68+ipRXJxtzLvuIsrL02+7YgVRwYK2thkZRH/84dhm506i3r2JrrmG/+7aZc48hZABYC3pyFRFetYBAaZ+/fq0du3aoI8r+MDcucDtt2vrX3uNLTfcQcT644kTeSX5wANAnz686v3wQ+DYMVYT2F+HCQlAtWrAli22uvLlgT//1FpxfP890K6ddtxq1dgj8/x54IsvWK9dtSo/JeTlcVAuPcuUevWAP/5wrFMKKF6cV+cAcOONrDK6eJFX5wUKuP8/GHHhAjBvHq/4O3QAjh9n3fzVVwN167o+9+RJYNYs/r927uz4vzl9moN+HTliqytRgv/XrixhhLBGKbWOiOprDuhJ90AXWZFHACdOEN1+O5FS+ivU/v1973vKFP0+7ctrr3Fp3ZromWeIDh3S7+vvv437GDCAKD5e/1h8PNGMGdr+MjK0bfX+B48+6vvnt7J5M1Hx4rY+09N55X/gANHZs/71PXas/uceN87/eQshAwYrchHkgj533OFa0M6e7XvfDRu6F+QAUVIS/y1Zkujzz437q1DBs/6cS3o60Q8/sPqIiCgnx/Nzq1YlWruWqE0bVgu1b8//k759iZo3J3rlFaIzZ1z/H269Vdtvaqrtsw8YQJSb69v/+I039Of99tu+9SeEBSLIBebgQaLz5123yckxXslaLERPPunb2FeuEA0caLzKd1dWrtTv97XXtG29GaNBA6LbbuMnkKJFPTvnxhtZP+5qzBtu0M511Squj4931IcblZEjfftfb92q7T8ujmj7dt/684YDB4jWr+frSDAVEeSxzl9/EV13HX/lyclETz9tvJGWl+e4kWYtJUt6vsmpx8svGwusggX5R8KVUCtfnujyZW2/p08T1avnKFCzsnz7sfCkxMcT9ejhWdtly2zzPHGCqHBh78a68Ubf/98TJtj+D0WLEk2c6HtfnpCbS3T//bYfkNKliZYuDeyYMYYI8lgmP99mZWFfxo41PufZZ7XtR4zQ7/u331hguVMDVKyoL6xq1PB8Bf3pp/p95+YSzZlD9NFHbK2ht0r3tRQowNY0NWoQ3X030SefGD+xOJfp021z/Owz78du08b99+uKS5fYWkXvB9Bs9PTyxYsHZ+wYwUiQix15LLBpk35M7W+/NT5n+HDgjTfYeqJ2bbavHjTIsc3eveyWfsMNbO9cvbpxTG/AMTCWlcREPo/Is89i75wDsNXH9u3sBNShA9u6p6QAL7zAJTOTvUZr1WL3dW8Cb1nt3S9d4s+6dSt7U86aZbMrd0VyMtudu6NZM+C55/Tt0/Vs7r0hKQmoXJn/z4Fm7lxt3ZEjwJo1gR871tGT7oEusiIPMrt366/2OnXyr9+uXbV9tmxJNHky0bXXEpUqRfTAA2wXXqECUZEi2vb33EOUkuL5CnXYMNv4I0faVEBZWUSVK/NrpfizWfcCRo/2bSVu7c++xMURVark/lyLhWjSJMf/18mTWtVKUhKvmE+c0O/nqqv8+46Cyb336n+GHTtCPbOoAaJaiXHatnW8uZQi+ukn//rU06N7WlJSiJ54gq1R9I5nZBD17OlYV7Ei0bFjPPaKFe7HeOopVv3Ym/h5WpRy1LvbF3cWPdYyejTRn38SLVxIlJ3N8161ivXe8fFEdesS/fgj1589q+/YVKEC0bZtvn0/hw7x/7BQIVatffGFf9+3O9avt1kambVYEBwQQR7rnDvHG5zVqxM1bcr6ZH+55hqt4ClQwHNh+f77XPSOPfQQ672nTiVq0sS2aq9Th1f4aWnu+y9RgnXEvvzQtGrF3qHO9VdfTXTkCP9110fZsrbXGRks0F1x333GfTVrxk8jW7d6/v3omXma8b27Ys0a/vFo1oxNIC9eDOx4MYYIcsF8vvlGu0l51VWeC8t69Vjt42ytopTN3fzTT30TxNbSogWvgL05p3RpflrJySHq18+2Uq5dm2jLFp7Xu+96P5fMTEdBvGkT0QsvsDXP9u0s9AYPZpWUUR8WC6uu3LFpk/75HTuafx0IQSNgghxAWQC/ANgKYAuAge7OEUEeYVy6xHFOunUjeu45ov37bcfmzWNnmB492LnmoYc8F2ylSnEfkyaxkAPYSWfMGFv/pUv7J8gBfnLQU1sYCVvr69q1WcCePEn0zz+O/5O33/Z9PjffzKtvezvvpCT2ujx4kOitt9z/39zZaG/YoH9u+/amXBJCaAikIC8JoN5/rwsC2AHgalfniCCPMNq1cxQGqakstN97z+aJmJ7OOtj//c9zgVa0qG2MS5d4tXrhguPYzjrXYJeUFFbjFClC1LkzP4XccgsLU18dm1wVpfQ3hZ3LwYPuv7fatbXnTZtm6qUhBJegqVYAzAbQylUbEeQRxJo1nguhuDii55/3TnBNmeJ6/GLFQivIw7GULeuZ6/6//7LHqlK8XzBqlDnXhBAyjAS5qXbkSqkKAK4FsErnWD+l1Fql1Npjx46ZOawQSPbs8bxtfj7wzjve9e+cWWfrVs5ZOWgQR+s7elT/vIwMID4GMhVWqACkpdneJyQAo0ZxxEMjduzgGOq33AKkpnLUyEOHgIEDAz5dIUToSXdfCoA0AOsAdHHXVlbkIWT1arbG+O47z1Z1R44QJSYGbnVZrBjRkCEcQuDmm0O/2g23Urw4r6jj4zlGy44d7Lnavz//v4YMcQzOdfq01twyM5MDfHXowKqiGjVYRSREHAhkPHKlVAKAeQB+JKL33LWXeORB5PRpXvX+9RcnOV640HbshhuARYvc58scPx7o14/jeAeKpCROsya45rHHOHn1iRO2uqwsFtnFinHcdPscqfZtjh+3vVcKWLrUODNRsDl7Fhg9Gli5kr2Jn3jCOy/cGCFg8cgBKAATAYzy9BxZkQeJM2f0Y6zYF6uFyPnzvJozYs0aR4sOX0qpUu7nI8V1SU83r6/77w/OdeiO3Fyi+vUd51a2LHu7Cg4ggDryGwH0BtBCKbXhv9LWhH4Ff/nqK/0YK/asWsVZezIygCJFgLZtbZlw7KlfH9i/n9t6ir0eNy6O9bvu5hPJuNJb+0Kczu1plJjaF6xPWCtXApMmcTwZPS5fBrKzzRvXmQULAOcn9H37OKuU4BF+C3IiWkFEiohqE1Hd/8r3ZkxO8JO//3bf5sAB4PPPgStXeC30ww9AyZIs0P/917FtgQJAy5aej28viPLzvds4jUSKFze3v/x8bd1115nX/513Au3bs4qtd2+gUiVOZG3lyhXgkUc4gXbBgpxOzl6lYxb79+vX79tn/lhRikQ/jGaaN3d9vG5dYNcu/WM//MDRBMlpD8U+B6Q7cnI8bxsN2OugzSYxkfXjBw7431dKCucx3bsXmD/fVp+XBzz7rG2MV18Fxo7l6I/5+cB333n3ROYprVrpP33ceqv5Y0UrevqWQBfRkQeJ/HzHiHRxcfz+5Zc5ycDFi0S1arnWo65da+vv7rtDryOO5WKWA5I1WFqfPvrHrTHU9VLoxcW5zzDljl27iJYvd4zD8tFHnPAE4DAEzzzj3xhRCgx05DFgiBvDKMUrr6efZquVhg3ZLtmeRx4BHn7YfV+LFrEeVQgdRL6dl5zMq/DSpYHnn+cVMABUrarf/s8/eTWsZ82UkOD7XkBODtCnD/DNN/w+M5MtcFq35uuwZ09gwwaOT1+6tG9jxCp60j3QRVbkAWLtWqKvv9bGBXHHmDGOkfqcV4AWi7nWElKCW55+Wv97P3rUOHF1pUpEr7+ure/fX78vT/Jz6kW6zMri8AyCR0CiHwaY/HyOfxGKtFa5uY45JOPiiIYP9+zcuXM5jRnAcVOKFQusA5CU4JYePdiBaMYMW+RGe44fJ3rxRX21jX1y6YIFOWCa8/U9YwZR1arcpkED4wTZRBwuQG+Oy5cbn3PpEqeQ692bf1iOH/fsuo5SRJAHksWLbRdzZiZHCjSLBQs4kUH37kTz59vq16zh+NaXLnG8EuebQynWRbri8GFt/HClOMlxqAWQFP9L4cKcZ9U+TPDtt7OO3P7aWLTIs/7ef9/x+tm4URuCOD3d0dPUnr599ft1dZ22auXYtlIljkYZo4ggDxSnT+tnyrHPnu4rX36p7XfUKMf42kWLclQ+vRukenXOwGOEUTLgjIzQCyEp/pcSJYw3SJXiNHuXLxMNHepZf02a8NOfdVVuFCDtq6/0r7dNm7Rp/Xr2NL4+ly7V7/+993y+pSIdI0Eum53+smABcO6ctn7aNP/dn4cN09a98gq7M1s5dkybkNjK9u1A374cMMli4eS4xYsDTz0FNGli7HJ/4YJ/8xbCAz3HLitEwIQJ7Kbv7C9gxMqVQKFC7CDUuTNQrpx+O6NgZjVrsgPa+++zieNttwEPPWQ8npEfhJHJbAwjgtxf7CPTeVLvDXqOEvZC3Mrx42yNYnRDDh3qGMdkzhy+Cf/5hy1biBzbS8yT2MFTIQ7wD7/Vw/Pbb4Gbb2b79itXbG2KFWP/AyNq1gQ+/dSz8Zo1078+b77Z8znHCOIQ5C+tWwNVqjjWJScD99/v+rxLl/gmMloV//STvmdf0aLaupQUYNkyYMQINg9zxlkw5+WxEAe0N4kgeMqSJcDUqRy+ITGRQ+ZmZABff21O/5UrA2+/7Wju2Ls30KWLOf1HE3r6lkCXqNKRExHt3cv6xkqViNq0IfrtN9ft33/fZs5XtqzjJiYRZ8nRyxJTsiQnI3a2KnnhBdu5vXuHXjcrJTZKfDzR2bNEjz6qPTZ2rHn317597KS0aZN5fUYoCGQYW2+J6TC2K1ZodefJyewunZXF7xctsjlt2NO5MzBzJjtNfPwxcOoU0LUrP35+9hmQmwt06sSqk0WL+JwKFThmRSBD0AqxyV138XVXpAg/YdpzzTXA5s2hmVcUYxTGVlQrwWbmTG3dxYsc28SKUfClkiX5b926LMinTuWN1h49WBWzeDHw+OOsatm9G3jtNdfqG0FISfE+omJSEtC/P8e5/+UXrRAHWHW3bBmvzyOFvDxeHEUgIsiDTeHC7utr1QLatHE8npICDBigPW/ECG3dnDmsR/zf/3yephAjpKV5L2wvX+YQs598Agwfrt8mO5uDttWv7xgxccECoEUL4KqrgCef1N+8DwXjxwNly7KOv2ZNY0uwcEVP3xLoEnU6ck/Zt4+obVutPrFKFaIrVxzbZmcTDRtG1Lw56703bNDvU2y+pQSixMXxnkxqqut2nngBP/44X6tLl3K/9sdatgzoLecRK1dq7e0LFTJ2bAohCEbyZcEFREC7dsD3TqHa27bl3X89axPrZZWfr2/BAgD16pk+VUHAww8DO3eyL0OlSsbt7E0PjfjlF/47dqz2Ov75Z+D663kl3KoVsG6d73P2lalT+T6z5+xZ7b0axogdeSBZvpzVGzt2cB7CjRu1bYoUcYz0tm0bO2pMnWozEQRYt756NT/22WOkqhEET6lalfdS7OPHf/QRMGsW53l1h8Xieh/Gap5rpEZZtYr/LlrEiUt27bJt/AcDo5y1KSnBm4OfyIo8UGzbxiuMxYvZseenn/Tb5ebaXv/8M1CnDvDmm45CHOAN0dGjtecbrdQFwVN27tRPAuKJEFcKmDHDWPAmJ3PoXADo3t19f2fOsFd0MLn3Xs5+ZU+FCtp9qjBGBHmg+OILzzwkrVl8vv+eLxxXj6pHj9peX7nC1iqzZvk3T0Hwh1dfBTp25DAQy5Zx7s25c1k4PvkksH490KABt73nHmDwYNtK12qF5Ywn6hozqVaNF1otWgClSvF99fPP7OQUIYhqJVB4mqw2M5MFdOfOjqtzPdq3t70eMwaYPt33+QmCL7RqxSZ6JUsCQ4bYcojGx9v8I2bO5Fgqu3fz/k/16rxyVwp44w3gpZd45X3lCh+zF9yJiUC3bkH/WGjalIV3hCKCPFDccQfw4Yeu28TFsXnWvHmuVyEWCwe/uu8+W92CBebMUxC8YeFCFt4TJvD+DsAbhYsWAb//zosS++t+2TIW6kOH2upSU7nk57O6xd6lPykJOHkSKFMmOJ8nShBBHiiaNGE725df5ou7Rg1ebfz5p63NCy+wRcD69cb9jBjBHnTOj6GlSgVm3oLgjnXrgH79eKP93DmOsrhsmXH7Dz7gqJ32lllXrnD0w8WLHdueO8f3xbx5gZl7lCIu+oFi6VJg0yb2wqxRg1Uo+fnAjz/yrvxNN7HjD8CecSVK8OOmPVdfDWzZwq/HjeMfhsuXWbC3acMroxB8f4LgNUuXcjRDKxMnss5cjxIlWOcuaDBy0dcYlgejRLVDUF4eUdeujs4FDzzg/rzt24kqVrSdU60a12VnEz38sNbJ4p57Qu80IkWKN6V2bb4/iIgee8y4XWYmJ5wwIzlLlAEJmhUk5szhXXxnfv0VuOEG9+cfPcqPl5Urc4CtTp0cXZyt6MVpFoRwp1s3vr7/+MPRCkuPuDhg/vyIMgMMNBI0K1hYnRs8rbfn6FG2Hy9ZkoX0fffpC3FAhLgQvqSns6emHt9+y+pFd0IcYFVk377Ab7+ZO78oRAS52Th7XjrXz5nDbvk33cQuy9nZrCMfPJh36q+/nj09P/xQUloJkcnp02x54ikWi3EExoMH2XBg/Hhz5haliGrFbK5cYTfjFStsdW3b8i78tGlAz57ac+LjtTbkCQlcPLVHF4RIZeZMVp80buxo1WVP6dIcsz8utteeoloJJHl5bCdbrRpQuzYL7okTgWeeYeE9ezavOIYM0T9fzxEoJwe49dbAzlsQgkFcHIfL1RPCGRl8nScnc1Jwo/gmBw7IosYFsiI3g+ee49yC9gwbxvawVvLy2NnBmyQPixZxoKEZM4Dz5/kHQRCiiZkz2at5wwagUSNjx7h69UITGTHMkBV5oCDibD3OjB3r+H7nTu+EeJkyHPuhc2dg0iRO6SYI0YbVg/O114yFeEYGh6QQDBFB7i9E+qmu7B8D8/L4Ik1K8rzfhg0dN4BatGBduiBEC8WLs+67USPj4G9DhrBuvFGj4M4twhBB7i9xcRxXxZk77+S/c+ZwSMw6dbzbqGnY0PF96dKc6i3GN3uECKJyZeNjKSnAp59yCIvVq/XbZGRwCNzU1MDML4qQJZ4ZvPceB+b/9Vc2perZk3Xmhw9zSExrONuLF/lvmzaug15lZAAPPcSvL1/mH4PlyznYv8QfF8KduDgOL7Fvn3Gb7GwOdetqYXLttd49xcYwIsj9ZccOoHVrYM8efl+8OG9+pqQAkyfrxySvUwe4cME4wWv79uxANHEixymP0MzeQoySnw9s3uy+3cmT2oQO9hit1AUN8pzuL088YRPiAAf7efRRfm2Uhq1wYX6sNHr0TE1lk6yvvxYhLkQ3ly4Zr8olwqfHiCD3lyVLtHVLl/ImqJ53W2oq0KcPB9TfsYNX3XXq8LHSpXl3fsqUgE5ZEMKGwoVtMc6dGTyY/27cyKkTBUNEkPuL3qramnXcPpi+lcxMW2zxuDigd2+2ob10iXWKnTqxi7MgxAJnzwLHjnGKuG+/BW6+GbjqKg5xW7UqL3Lq1OFQ0M2aAcePh3rGYYkIcn959VXHR0Ol2Cb20iX2RnNm717g7ru19UlJrEaZPNk47oQgRBtEvLH/9dd8zyxbxqvvCRM4n+3Gjba2y5ezt7SgQTw7zWD1ar7w8vJYbWINV9uokfGGzapVNhPDEyc4XOegQbw6EYRYJDHRfeLlzEy+xyZM4AVU37785GuffSiKMfLsFKsVM2jYUGv3DbDHZ8OG+rFUNm/mYx9/zBumetYtghBLuBPiVkaOtL1+6y1bUucYRlQr9hBxJu0RI8zJqH3qlL4QB1iIb9oEPPywCHFB8JQLF7R1EuJWVuT/DxFn9J4xw1bXrRtHL3Slsz51itucPw906QJUrGg7NmeO8Xn5+TYzRUEQmIQEjvwJsOrE6gBXrBg7xOnl+RQnOVmR/z8//ugoxAHeRV+40PicrVvZjPChh1i/Xa2aYx9ZWcbnzp1rHHtZEGKVhASgVStbsnIrR49yInI9QW6UxDmGEEFuZc0a/XpXm7Ivv+y4OZmby/ruvDzg77/5QjRyMS5WjL1ABUGwkZ3Niye9FIc//AC8+y7wyCPsj1GwIDBwIDB8ePDnGWaYolpRSo0H0B7AUSIyyHUW5lidcpypXdv4nPXrtXX793Oatqeesq0onBMllynDgbby820xVQRBcE2pUpyA4qOPuAj/j1kr8i8BRHaq63bttNm6b7uNs/0Yce212rpSpTj0pv1jIRE7NzRuzJubv/4KFCoE9O/PduPly5vzGQQhWrDGKbdisbCd+Qcf6G94mkF+vi2wXaRBRKYUABUAbPak7XXXXUdhSW4u0cyZRC+/TDRrFr/Pzyd65x2i6tWJKlUiGjKE6MoVbr95M1FmJhGLaqL4eKK337a9dy7330904YLjmHv2GLeXIiVWy4wZROPGEd16K1FGhuOxBuXAYAcAACAASURBVA1s96BZjBpFVLQo99+kCdHWreb2bxIA1hLpyF+9Sl+KO0EOoB+AtQDWlitXLigf2hTeeEN7kT35pO348eNEL75I9MwzRDt3Ep09S5SW5voirViRaNo0Pn/48NDfNFKkhFNJTCSqWZOoTBlePOm1mTjRvHt89mxt/5Uq8UIuzAi5ILcvYbsi16NCBe2XnJZGlJdH9M8/RHXr2upbtiQ6eZJ/3d1drBYL0caNRB9+GPobR4qUSCutWpl3j3fvrj/G8uXmjWESRoJcrFbcoaczu3yZ9WkPPMABr6z8/DPHIh84EHj/fdf95uXZ2riKySwIgpYzZ8zrKzHRu/owRAS5O6wp2+zp1o3difW8P+fO5b8dO7rPsfn55+wUpJfzUxAEY8zM4fngg1qnv7p19cNuhCmmCHKl1BQAKwFUV0rtV0r1NaPfsGD4cOD++/nX2WJh788xY/h9erq2fbFi/Pezz4zd8wVB8J2CBTmK6DPPAP/8439/zZuz89+117IT3113cWauCEKiH3rKxYusOUtJsdW99hqbGtrz1VdsanjzzZFryiQI4Ujr1uwotG6drS49nSOMVq3q2DY7m23OoywktFH0Q1GteEpysqMQB4BXXuGAPS1bsr357Nn8a96tmwhxQTCT1q35SfiPPxzrT592jIa4ejVnG0pNBcqV44VVDCBBs/zlvvu4WPn3X/buFATBHDIyOJzFhg36AbKs6pXz59mJz5picf9+jsNSrZq5OvUwRFbkZqOnNxcEwTVGCZgBFsxffQU8+yyvtJ2xWn398IM2Ty4Re09HOSLIzSY9nVcQgiB4jiehaHfv1nfPX7CAdedG5oJGgeuiCBHkgcCahkoQhMBz6RLrztu0Yb24PYmJwL33Bnb8/fu1TwJBRqRNIGjfnnNyJieHeiaCEP3ExXFegKQk9u1o3JjjmqelAb17ay1azGL3buD664GyZVmHf999nqerM5noEuTnz7O3ZUoKqziefTa4ttwXLgAvvMDOBC+9xJ6b7dvzBSUIQmB4/HHWoVeowE48K1dylqHz59npbsCAwIzbqxcv2ACWM19+Cbz5ZmDGcoee336gS8BirdxzjzZewssvB2YsPdq1cxzbYiFasYLoiy9CH5tCipRoLLVrE735pus2CQkcA8lMDhzQH6tWLXPHcQJRH2vlyhVgyhRt/YQJwRl/1y5g/nzHurw84NVXecUgCILvGDn2NGnCydJdkZMDnDtn7nxSUvRDcITIai16BLlS+l+2L55d69bxxkmJEqwa2bTJ/TmnTunXL1xo/kUkCNGM3j3bqBGrLe0pVowXSUb3npWrr+aELoULAw0asJWLv6SnA336aOtDtWjTW6YHugRMtdK3r/ZR59VXvevj0CGiQoUc+8jMdP9olptLVKxY6B81pUiJxpKQwGGfO3bkJC+9ehEdPszF3bkFC2r72rTJdzlj5fJloqFDierUIWrWjOjbb/3v0w0IdDxyb0rABHl2NtGjjxKlpxMVL876cVfB4c+fJ1q5ki8GKyNH6l8M48bp93HyJGcKyskhats29Be8FCnRWFJSiJKTHeuee47vwRo1vO9v0CDz5E4QMRLk0aNaAdjcb/RoftQ6fJiDWlks+m2nTAFKl2ZTpbJlOY44YBwjRa9+8GCgZEmgZk3eMffEqUEQBO8pVEh7D44axbFWxo/nqIXeEGWRSaNLkHvK4cPsJGANTp+TA7z9Nm9Wduum3cRISuJkEhkZrL9LTeWEyW+9xfUAcOAAsHhxUD+GIMQMJ05o6y5fZjvu7Gw29/WUuDgObhdFxKYgX7RI33B//nx2Hpg2DahYkeuqVgVefplt0q2bKtnZHA/ZmRA5AwhC1JOTo1+/fTsbJIwda3zuk0/agmZVrAhMmgTU10SCjWhiM/phiRL69SVL8t/OnYFOnYCzZ9l76557gjc3QRD0qVCBo4s6c/EicOyY/jnLlgFNm/LrnBz2+IxCYnNF3qIFxyy2JykJuPVW23ulgOnTuZ0n5oeCIASWJ5+0PSk7o2fi266dTYgDxkL8xAl+yo5gIkeQX7oEfPMNP0Lt2+dfX3FxQOXKjnWXL7NzwcKFtvfPPcd73IIghJ7p0/V15YCj6kUpTvQyfrzr/v79F2jWjDdKMzOBJ55gJ75IRM+UJdDFa/PDffuIKla0mQ4lJBBNneqL9Q6zfz9RXJy+WVLt2txm9+7Qm1xJkSLF+6Jnz334MJsb29Oggfbcd9/1Xa4EAUS0+eFrrzkmWc3JAR57zPfNxRMnjE0FN25ki5T69bWp3QDgzjuBNWs4OE9qqn6ge0EQQof1qRrgAFpZWbwvVrAgq1rmz+dsQ2vWaM+dOjV48zSRyBDkq1dr644e1d/48ISaNY11bQBbpJw8yXoz+7jiZcsCw4dzdMPPPmN3X71A94IghI6CBW2vW7a0qWOIgBUr2MqlcWP9UAARGqk0MgR5zZrausKFgTJlfOsvLo5NDEuVcqw3issyYgSbLG3bBmzdyiv22rWBJUt8G18QhMDRtSv/nTPH2MHv0iX9+z1QIW8DTGSYH770ki2dk5Vhw/RVH55Svz5vmi5bBvz2G/+Kf/ghsGOHYzsioG9foEgR9iLr3l1W4YIQapTie9OeIkWAoUPZSQhwb2qYnw907Ahs2cLOfk88AXTpEpj5BpjIEORXXQX89RcHjz95ku28zTDoj4sDbrqJC8C/3lZXfSvt2/MFArAjkQhxQQg9RCwD4uKA1q2Bnj3ZEs2aiBngCKYFC7qOPjpkCHDttYGfb4CJDEEOcMjKp58O7BiDBvGqe8wYFthdu/Lro0d51R7ivHyCINhRtiwwc6bx8exs9gMxUoHefntUCHEgkgR5MIiL483MYcP4sctiAT79FHj0UXG/F4Rwo3Bh18dff11fiDdowCv38+dZFfPkkxG7yWlFBLkeSrEQP3JEhLgghCuHD7s+Pm+efv369TbHn8WLgR9/ZGuWCCYyrFZCxYoVIsQFIVxZsACYO5dfr1zJyc7tI5BmZuqf5+y9+euvwPLlgZljkJAVuR7Z2cDEiRKWVhDCnY4d2VZ80SLHuhkzOF/A0qWe9XPwYGDmFyREkDtz6RJ7f61fH+qZCILgDiJHIQ4As2fzSr1TJ+CDD4CnnnKdSCIpiX8MIpjYU63k5wNffsn24I89xvGM7fnmG2MhftNNvHkydixfPM6BtwRBCA/ef5+tUvbvd602KVyY5YG3GYbCDb0ALIEuAcvZ6Qn9+jkGyUlLI1q6lGjHDj4+eLBxMJ6kJKLjx7ndsWOhDw4kRYoU96VsWaJ69bT13bppA2mFOYjooFlmceAAx0ix5/x5XmlXqwbUqsW2qUZcvsyen6dPs6NB8eIBna4gCCawbx9wyy1A9eq2ultuYVkQJUHvYkuQ79+vH/WQiP9u3gx8/LFxPr/kZODBB9nTs0IF9igTBCE8iHMhzq5c4ThJ69ezOnXhQvd26BFEbAnyunXd68I2beJNkkcecTRfUorLli38/vBhDhlQp07g5isIgmeUKWMcmtp6HGBPzmrVgjOnIKLIuhoNIvXr16e1a9cGfVwAwA8/cEzxM2f0jzsH42nUiNtbLLw5KghC+FGyJHDokOs2NWrw/V++fHDmFACUUuuISBNoKrZW5ABw222sK1+wAPjoI20oS+cftlWr+IsvVCh4cxQEwTucQ1LrsXUrMHBg4OcSAmLTjjw11ZZouVIlNlU6eRIoWpSzhzjTpw9HUIuLc/34JghCaOjdG1i3zn27n34K/FxCQOytyJ1p04Yft1atAh5/XL+NNQymCHFBCD8GDWK1iScUKsTxyuvXB0aN0j6BRyixuSI3olUrzr79/fehnokgCJ5w3XVA//7A//7HiSRycly3P3KEC8Ar+JMnOSdwhBN7m53uyM0FvvuOE7NOmqQfgyElhU2X3G2uCIIQePSyBdmTkMA25EePcrEnPZ2FuVGaxzBDNjs9JT4e6NYNeOstdhDSIzubd8kFQQg97hajOTnsI3LsmPbYuXPaaIgRiAhyV7hKxLp5c/DmIQiC/+gJ/I4defEW4YggNyI/HyhdmjOI6D12SZxyQYg8ypa1eYC2asWe3J4weTI7FcXFAVWqGCetCBGm/BQppdoAeB+ABcBnRPSmGf2GjI0bOZv233/z+5o1OfmzWK0IQmTz2mucUD0vz/NYSVOmOIbt+PtvXslv3Ahcc01g5uklfq/IlVIWAB8BuA3A1QDuVEpd7W+/IaVXL5sQB1iN8tBDwPjxHKtcEITQkpDg/TnduwN3381hOrwJePfMM9q6/Hw2hggTzFCtNASwi4h2E9EVAN8A6GhCv6Hh339t8VTsWb4cuO8+4MUXgz4lQRCcKFSIE0J4wu238+p52jTf9OFGuUHDyNLFDEFeGsA+u/f7/6uLTIoUARITtfUlSvDfG2+M6FgNghAVnDjBYaVLeyBqbrjB2ALNE5o00dYpxav7MMEMQa73s6TZHlZK9VNKrVVKrT2mZwZkJkeOAG+/DTz7LPDbb96dW7gwq1HsiYuzPV7ddBOwZ48p0xRCzzFkai9WIXI4cMB1pq5ixYB77/VvjI8+coyaGh8PjB4NXB1GGmS9bBPeFACNAfxo9/55AM+7OiegGYK2bSPKzHTMBDJihHd95OURjR1LdMstnEVkyRKuz80lUir0GU+kmFLOI4UG4/WQz0OKn6VUKf16pYhmzjRHrly8SPTdd0SzZvHrEAGDDEF+e3YqpeIB7ADQEsABAGsA9CIiHUUzE1DPznvvBSZMcKxLS2MPzYIF/es7J0df7SJELLmwIB6R7xASsyQmujYF7tiRPbWjhIB5dhJRLoBHAfwIYCuAaa6EeMDR26g8fx7Yu9f/vhMSODmFEDWIEI9w3Plz2FufRTGmOAQR0fdEVI2IKhPRMDP69Jnrr9fWZWWxEb8ZLFrkm+mTIAjBp0WLUM8gKESfZ+eLLzqmckpI4I0JT02V3JGZCTRoYE5fgiAEjquuAl55JdSzCAqRH2TAmRIlOO/mvHlsotS2rWcmSt4weDDbpgqCEL5MneqYdzeKib4VOcAbIF26cMZ7s4U4AHToEBUxjAUhqrFYQj0DDgXw1lsc5qNuXTZlDADRKcgDARHwySesc+vQgQPaDxlifLFkZISV55cgxBzLlgFz5njvS2ImL77IT/BbtgB//gk8+ijw3numDyOJJfTYtYtNGC9eBHr25LRQgwfzL6s91aoBp09rg9ULghAcrr8e+P13/WOpqcCFC7Z2P/zAiSSCRX4+j2dNFWmlYkVg926fujQyP9QYlgejBNQhyF+WLydKTnZ0KvjsM6L4+NA7PkiRIsVWChUi2rKFqGdP7TH7e9haBg0KrizJydGXG0WL+twlDByCYle1QqRf/8orvBK3b/fcc5wCThCE8KBUKWDfPnaTnzIFWLmS98QefJD10Pb3sJV33gF69ODUbsEgPh7o3Flb36OH+WPpSfdAl5CtyPPyiF58kahIEaLERKK77yY6dcqxTcmSoV9pSJEixXX58UeikSOJBg8mWrHC8R4+eZIoKcn43E6dgidzjh8nat+en+wtFqI77yQ6f97n7mCwItdUBKOETJC/9Zb2S+3SxbFN167aNpUqhf7ClSJFCpeKFbULruHDbffw8uVE9esbnx8XR3TmTHBlz4kTRKdP+92NkSCPLdWKcwwWgOMwnDlje//mm5wOykrBgsAXXwC9ezuel5ERmDkKguCaKlWAQ4cc64YOBU6dYrVK06aAK2MKiyX4pokZGRxZNUBEn0OQK+J0freUcjQTrFIF2LEDmD+f9Wzt2/POc9OmnO7pt9+AGjWA778HvvoqeHMXBIHRS3x+8SKwbZtnnpy9erFFSxQRW4K8b1/gyScd67p142wj9hQoAHTt6linFHDrrVxWrBAhLgihwnk1DvCTc1aW8UZm0aK8Cu/ZExg+PLDzCwGxJcgHDuRf7o8+4oiId9wBvPuu9/2sWWP+3ARB8J0ePfhpunRpTjbhTIMG/JQdpcSWjlwp4Pnngf372ZHnk084Vrm3+JM2ShAE87n2Wr6/X39d/3gwHYFCQGwJcrO46ioJZSsI4cTevUB2Nq+8b7rJ8Vh8PLvGRzGxpVoxi08/5WxBgiCEB8uXAyVLAmfPcgL1du3YYah0ac6327hxqGcYUESQ+8LBg6GegSAI9qxcaXt96hTHVdm+3byEMmGOqFZ8oU2bUM9AEASA1SZ6q+38fGDu3ODPJ0SIIPeFLl2AAQNCPQtBiF0KFeINzKws4MgR/TYxklQCiFVBfukSZw/yFaWADz8EypQxb06CIHjO2bNseXb4sH5I2HLltL4gUUzsCfKXXmLngKwsjlG8fbt35+/aZbtwatQwf36CIPhHqVLstBdl3puuiC1B/sUXwLBh7AwEAKtWAZ06cSgddxw+DNx4I1C1KlC5MtC8OfD444GdryAIWtxl3ipb1jFeUgwQW4J86lRt3bZtnILJHQMGOKaMWraM3f1vvNG8+QmC4JqUFPeJz++4IzhzCSNiS5CnpHhXb4/eDviuXcC//wLJyX5NSxAED8nO5tV2u3bGbdat4yim/uyDRRixJcgfflj7WNayJefedEeRIvr1Bw7oZwERBCEwLFkCzJplHF7j6685FEeDBjEjzGNLkLdqBcyYwXanFSqw2+6MGa7POX8e6NjRdYLlEiVMnaYgCC7YvBmoXp29Nl3xzz/AuHHBmVOIiT3Pzs6dvVtBv/IKMGeO8fH4eBb0H30EXL7s//wEQXDPP/94Fu9o587AzyUMiK0VuS+4EuJZWZx1aOBAEeKCEGw8iXfUvHng5xEGRLcgP30aeOopoGZN4LbbWLfmLVlZ2rqUFOCPPzgc7o4dwIYNfk9VEAQfcE65aL8H1rEjZwMyYu9e4Jdf2Lko0tFL5BnoErTky02bOiZdTUggWrPGuz6mTdMmb336aT527JjrbN1SpEgJXGnalOjQIaLXXyd6+GGiOXOI9uwh+vprotWrtfdyfj7/XbiQqGhRWz/JyUQTJ/ona4IEDJIvayqCUYIiyP/4Q//Lv/9+7/uaM4eoTRuiJk2IRo8mysvj+l9+8e0CTE0N/U0gRUqkl3vvJcrNdX//HjxI1KkTkcVClJlJFBen7SsxkejIEe9lQ5AxEuTRu9l55ox39a7o0IGLM9dcAyQmAleueNffhQvez0EQBEe+/BKoX98WwO7SJb4fnZOsd+tmc+YzMke8coXVLBHqTBS9OvLGjTnQvDNmBtLJyuKg9YIghIY33wS2bgWaNWPHvFKlgNGjbcd37XL0yHaFO3PGMCZ6BXliIvDdd2xvCvAG5QsvAHfe6X/fR48CbdvyL/+wYVxXvHjE/poLQsSyfz9w882cIQjgkLaPPw7Mm8fvnVfnRtxwA9CkSWDmGASiV5ADQMOGHEvl33/5C7YKXX/p3JkzkNhz5AhbxcRHr7ZKEMISvXjkvXtzTJZDh7Q5PAEOQa0UxzR/5BFg4ULvxty+nUMB5Of7NGXT0VOcB7oEzWolEOzd63oDplGj0G8CSZESDaVBA//7SEggWrSI6O67iQoWJKpUiWjMGN/v/9OniVq2tPVfpQrRli3myRc3IOY2OwOFu0e1rVv5EW3FiuDMRxCijYQEdvZZs4afcMuWBU6e5LrsbG37zEzjTcycHGDyZOCrr8yZ28svAz//bHu/axdwzz081xAS3aqVQFCggOvjZ8+KEBcEf7D32MzNBfbtAwoX1hfiN94ILF3Ke1+FCun3d/KkeXP7/ntt3dq1rmMxBQER5N7ii/miIAi+k5vLXph6fPklmwFPngwcO6YfwK5LF/PmomcJl5YGFCxo3hg+IILcWypVYpd/b3noIeMVgyAI3lOkiGPe3MREYOZMzuAFAElJwKBBwN13mzfms89qQ2EPHBjynASiI/eFb74BundnfbhS/Ct98KBx+1tu4ce/aIjpIAjhgFLA8OGs6iQCfvwRWL0aqFWL78t//+XcvOnp5o7boQPryMeM4RDXd9wB3HuvuWP4gOKN0OBSv359Wrt2bdDHNZ2tW1l/t2MH8MYbWmH+2GNAnz5s6yrJJwTBHAoWZNPC48eBcuU4VePq1bbjLVuyebAnYW6DBRHw00+8KVqnDvuhWCxed6OUWkdE9Z3rZUXuD+PHA++8Y3y8cWN2IY6GHy1BCCbFigF9+wIffKANaUHEWYCM+Pln4PPPgQcf9ElYBoQ77gCmT7e9b9OGnZZMmp/oyH3lr79cC3EAGDmS/7Zp47mHmSAIbAVy//1sNmivf65cmVUa7nj4YaBiRWD27MDN0VMWL3YU4gCwYIGpcxPp4iv2j3JGWH9tK1TglFNJSQGdkiBEFStWcB6Bgwc5R+fvv+sHrzNi3z6gRw/OqxtK1q/Xr//jD9OG8EuQK6W6K6W2KKXylVIavU1U44nlyuDBttd9+wKdOrluHx8fXno9QQgl993HK/BPP+U8nUT8dOsNV66EflVer55+/bXXmjaEvyvyzQC6AFhmwlwig3XrOBhXxYra7CPWFXhyMvDqq0DTpo4p4Fq0cN13bq65cxWEcObOO7WmfM4cPMgmfy+/zHtOy5dzxFFvdMshtvFGixZs5WZPmzacwcgs9Pz2vS0AlgCo72n7iIy1cvEi0W232WIsFChANGkSJ5145hmizz8nys4mOn+eYy/UqsXtChYkGjKE+7h8WZu1SIqUWC1PPMGZfbw5x2Ih2reP6MABoltucd++VCm+J0NNfj7RggWczWjOHFtyGi9BIDMExYQgf+897UWSlkZ05oxju7179dO/TZjAxz/+OPQ3kBQp4VIyM70/5/vv+V4aN07/eMuWHMyqVy+inTuDKycCjJEgd6taUUotUkpt1ilePRcopfoppdYqpdYeO3bM2weH0PPLL9q68+e1m569ejmqU6xMncoZTJ5/PjDzE4RIxCjYlRHx8UDt2vy6d29tDPHHHwcWLQJ27mQTxSpVzJlnmOPWjpyIbjFjICIaB2AcwA5BZvQZVCpV0q/PyuK/RLy5aRQwKyWF9X2nTgVmfoIQC7z0ki2TT4ECnANg3jyOD960KevRYxBxCPKUgQOBDz8E8vIc6ydNAurWBSZOBN5+2/j8hx7icJypqZKzUxC8pVgxYO5cThZjj8Vi7qZhhOKv+WFnpdR+AI0BzFdK/WjOtMKQ+HitEAdsmYJmzDA+t39/dhv+4w8R4oLgC0eP6mcCEgD4uSInolkAZpk0l/CmSBE2K7x40bG+VCn+W7iw/nmjRvFqHnAMSC8Igj5KsarSGaNQtoJ4dnpMWhoHwbLHYgGee45fP/KINl9ny5Y2IQ6w7bkgCK7RE+JKAa1aBX8uEYIIcm94803giy+A9u2Bu+4Cli3jELUAb7IsWAC0bs2B7gcN4tjI9nTurB/4XhBiFaWAdu0c44o7k5QEvPceUK1a8OYVYUgY22Bz5gybKC5bxiv4evV45z1csnELQrCxWNgIQC9e/3ffATfcwLHFBcMwtrIiDzaFCwPz5wPnzrEp4s8/c5jbsmVDPTNBCA15efpCvHVroGpVjvXfrx/bh3vCuHFsinjTTeYlXQ539LyEAl0i0rPTW9avJxo1imj+fM/ccSdMCL2XnRQpgSzXX09UtKhnbZs0IZoxQ+slPWqU6/vo9de1fX34oTn3dBiAQLroe1uiXpA//7zjhdSsGcdq0ePUKaIxYzgGi1Khv9mkSAlUiYvzrJ3FwvdGq1baY4UKEe3aZXzv6bn8V65s/j0eIowEuahWzGbnTt4UtWfZMmDCBG3b33/nWOWPPMJR3YiCMkVBCAme7gNdfTX/3bVLe+zsWVa39OypDYWRlwecPq09x9swABGICHKzWbtWXyDrJaJ4/HHe/BQEwcadd/Lf5s31jxNx7CLnDF0Wi37iiRjIlyuC3GyMEk7UquX4Pi+PE7EKguDIxx/z6n3YMKBGDeN28+Zp68aOBZo1s72/9Vbg3XfNn2OYIYLcbGrV4mxA9tSubavbvBnYto1XD9WrB39+ghBOFCigrdu7l9UqpUoBmzYB06bp57wtVkxbV6IEsHQp8O+/nOptwQL2yo5yJGhWIPj0U6BbN76gqlblR8Vjx3ilsGEDt7nxRuDFFznBrGQGEmIBiwVITwfKl+cUbo0a2aIX2pOcbHOcs1g4u84PP7Aznn1fTz5pPFb58qZPP5wRQR4IlOJUTvb5Bfv3twlxAPj1Vw6Nu3EjP0p+8IFnfVss+sG7BCHcSU8Hjh93rGvWjH0p7GMYPfUUUKiQY7tx4/hpd/ZsIDOTQ1/Yq1BiHPHsDAYrV7J3mjNFigAnT/IGaYMGnvVlFFBIEMKdihWB3bv5+v3gA+Dzz1kXftttnHTl1CnemOzaNdQzDVuMPDtlRR5o/vnHONhP8eL8t04dfpQ8fNh9fyLEhUjln39Y5z13LqsVrWzZwsmVR48O3dwiHNnsDDRffmkcg9waOTEhAZgyRRs9URDCGV+u1x07gI8+0tbr1QkeI4I80Jw/r1//xBPAvffa3t90E/Dss8GYkSCYQ24um9vWrMnB33r1YkuslBT99vHxHCVU756QhCt+IYI80HTvrq1LTeVHSWdefJFDelrRM80ShHBi82Yu69cDkycDf/8N9OkDjBnjGAguLg4YMYJNCnv00Pajd584c+oU8NZbbOn1+efAlSvmfY5IR89vP9Al6mOtOPPJJ0RZWba4Dz/95Lr99u1EK1cS7d9PVLNm6GNkSJHibenXj6/lX38lmjiRaM8e2/V99ixRjx4cUyUujqhzZ4455IpTp4iqVHEco21b/+7LCAQGsVbEaiVY5OSw6VWJEmx54glEbNkibvxCJJKezkkhTp3ipBCvv+6YKPnMGb7G09Pd9/Xee8DTT2vrV6xgn4wYQeKRh5qEBKBkSc+FOAD89ZcIcSFyOX2aEyZfucLql65dbb4U2dnsP9GrF8cc+ucfw2l4oQAAGThJREFU133pBdACOEidIII8rHHlWtyvH5CVFby5CLHH00/ru8b7Sl4eMHEiv27XDhg8mD02R48GGjYEDhwwPlfP+UcpcQr6DxHk4UypUhz0R49x47RecoJgJr/8wiEmzISIvZqXLHGsP36cr2kjunfn0LVW4uI4c1ClSubOL0IRw+Vw5957gR9/DPUshFhk/Xpz+7NYgN69jdUkrlbkFgv7Wjz3HAeda9w45uKpuEIEebhz1VWhnoEQK3gT/iEzkzcx3SWLKFiQbcSrVeOwtPXqsVliYqLWfNDo6dOeunW5CA6IaiXcmTtXPD6FwFO6NGer8pSsLOD994FvvwWuvVZ7PCGB46qcPcu68a1bgS5d+FjRosD48UBaGr+PiwMefpgjhgq+oWeTGOgSc3bkRuTlEU2eTNSnD9ELLxDt3et4fMoUrX1uairR0KFEKSmhtxWWIqVMGaKZM4mSkx3rn3vO/fV/9izRL7842pgLLoHk7Awhc+dyYKwCBYBbbmFTLAB46CE2v5o4ERg+nFc2u3fbzps2TdvXhQvc365dvOoRhFCyfz+bFdqHoQV41e2OggU5NEW5cgGZWiwhgjzQbNzIj5QbN3Ky2J9/5miI27cDn33m2PbECWDUKNv71FT9Plev5kdab13409KAu+/mWM+irhHMgkhb98or4gMRRESQB5qvvtJmADp8mJPH6t0Af/9te92/v7ED0ZIl7lPF2Qt6i4UjzH31Fd9g9mFEBcEd3tqTZ2c7Pl0KAUUEeaDRE9YAB9m3bvbYc/PNttdNmgCffKJ/fvnyHMfZFZcvs5XA229z+NA+fWxzev9993MXhORkVpO4s05xpnBhtlQRgoII8kBz1128GranaFHWK44b57hqvuUWYMAAx7YPPgjcc49jXYkSnAfUWS/pDBHbAg8dyjr1+fM5sP/Eiew+LQiuuOYajhF07Jh35ykFjBxprBoUTEeCZgWDWbOAF17gVXGTJrwattrCnjjBSZrLlGE3ZT3y81kVs3Ahr8T792dhXqcO6949wd5GuHx5YM8e/z+XELlUrcrOZvPns6D++2/vV91GTJ7MCw3BdIyCZokgj2Q2buQch1ZdZFyceTejEN2kpAC//w7ccQfbeJtFo0aco9ab4HCCx0j0w2ikdm1g+XKgcmV+740QdxWQS4h+srP5ic5bIe7K2um22/ipUYR40BFBHukMGeJo6QKw+7M7+vQBZs4E7ruPw4lK3IrYw9un8SpVeBU/ciSHZE5LA2rUAJ55Bli1Cvj+e7YNF4KOCPJIxzmKHMAxLGbPZpfnG27glbs9GRnAwIG8gsrL43jQojMX3DF+PHDddZxv9uBB1oVXr87XmztTWCGgiFdIpFO5sjaaXHo6ByC6/XZ+n5sLTJgA/PQTe9ENGMBxNe64Q997VBCc6dIFaNrU9r5TJ14sWBk7li2krrkm+HMTZLMz4lmxgs0WL1+21Y0YAQwa5P7chASts5IQnRQuzOaqniYsbteO439v2QJcfz3QoYPNKeivv/QFdsOGrGIJFCdP8j5QDCdUMdrslBV5pNOkCbBuHbv7nzvHGcpbt3Zsk53NMTEqVnSMz5KXF9y5CqHDmi/T07gm5cpxOAc9Fi3Sr9++3be5uePcOd7LmTWL9fq33ca+EJmZgRkvAhEdeTRwzTW8AfXZZ1ohPno0ZxqqXp3jQH/7re1YlSravlJSWIcuRA9K8T7J1197fo6rTUuj9GoVK3o3L08ZNAiYMYNX40S8qersOBfjiCCPZlav5o1Ma/CiI0c42qI1E8t33zmGCYiPB6ZPB/bt0xfyQmRiDclw6pRn7ePjeQVsRN26jqEkAPZedpWqzR+mT9fWWQW7AEAEeXQza5a2LicHmDePX199NQv3yZN5NX/wINC2La/K//qLvVGdwwvo4Ym5oxBaduzgDUojlOLv8ZprgAULHDNTnTnDK+By5YD69YFvvgEWLwY+/ZRX5716sVNagwaBmXuhQtq6tDRzE0NHOKIjj2aMnH7sVScpKfru1AkJtuwu7vB0A03wjsKFPQsFW7w4Z6P/809+Ajt3TtsmNZXzXI4aBTz7rPY769+fLU/0uOMOW97Yffv4eklLAx54gEugeewx7eb9o48GftxIQi/bRKCLZAgKEocOEWVkOGZuqVSJ6NIlz84vVcq7bDGJiaHPWBNrJTWVaO5c23dWp45+u7fesrU5dYqzUqWlcaaphx4iys7WvwZ279bvr107369LXxgzhui664jq1iV6913OrhWDwCBDkJgfRjtbt7LFwpYtHAfj5Zc5QJdeu48+4ljpeXmcqXzbNu/Gslg4suPhw/w+I4Mff48f9/9zCI5UrcpZpVq14pW7lZ9/ZnNUexISODCWfTt7Ll5ky6cyZbR5O41MDW++mdUrQlAxMj90u3p2VQCMALANwEYAswCke3KerMjDjPXrzc0B2ro10W+/2VZNFSuGfuUa6UUpoqef5pXwa68RnT5t/H1OnEhUvDiRxUJUvTrR6tXGbefNsz21KcUr9dxcxza1a2vnM26c/9ed4DUwWJH7JMD//2SgNYD4/16/BeAtT84TQR5GLFzIN703QqVAAX4cr12bKC5Oezwujmj/ftsY770XekF4zTWhn4M/5fbb+X+5YQNRly5ENWoQ3X8/0b59vn/3Fy4QFS6sHctZSP/9N1GzZnysUCGil14iys/3fVzBZwIiyB06AjoD+NqTtiLIw4TffuNVmzshkpREVLUqUePGREWLcl3BgkRDhxI1bap/zp9/2sbJziZq2TL0wtDMp45AF/sf16ZNeb9j3z4WpPbtKlUiunzZt+9/4UL9sTt31m9/5ozvYwmmYCTIzbTfuR/AD0YHlVL9lFJrlVJrj3mbcUQIDGPHemaVMmQIxz7ftcuWLebcOda3V6qkf85PP7FYAFhHm53t3dwsFq2+1l88mUPhwhyjxmraVqMGm9gF28Ty2DFg82bg33+BZcs4kciXX7IlkT27d7M/gC+UKuVdfaFCYmoaprgV5EqpRUqpzTqlo12bFwHkAjB0HSOicURUn4jqFy1a1JzZC/7hLBSsWIVYYiKbpQ0axFmM9H6A9+7V7+OZZzgDDQCMGcPJBrwhL4+FWLA5cwaYM8fmbLJ3LwtSZ3M9iwVIStKeb5ZtMxHHLZk0iWPOW+emx9tv214vWMCbnTVrAoMHAxcuGI9x9dXstm9PWpqY9kUiest0bwqAewCsBJDi6TmiWgkTJk3SPlYXLcqP0AcP8l8rv/2m/xjetq1rFUGFClp1gJ7Ko1w5oqys0Ks0PC1Vqnje1mIhKl+eNxM9PadIEcf3AwcafwcA0ebNRIsXa/cs2rd3fQ1cusTmfLfeSvTgg0RbtgT0khP8AwHa7GwD4C8ARb05TwR5GPHii2yLDBDVrEm0bp1+u/x8tuO1FxKJiUTTp+tvmHlTChQguniR6IUXQi+gPS3O9vlGpU8f1m3fcovrdtddR1SyJNt2X3+9fputW9mOWu/Y998Tde2qf+zvv4N7TQkBw0iQ+/sc+CGAggAWKqU2KKU+9rM/Idi8/jq76e/bB2zaBNSrp99OKX5sf+QRoFo1jnfeqxdHW/TE+9AVly4Bb73F/bvDKI1YcrJx0Kb0dN/nZjQHTwNELV3K6hGjiIEAe2ZOm8YhEs6dc3SPt2fDBuDhh7X1qamcQMToezBSoQnRg550D3SRFXmIuXKFV9IjRrANuS9s3Rr4Va+elUmbNq7P0VNfrFxJ1KOHOXNKSCD67juiOXP876tcOaJRo4iOH7f9X+fN43q99tu2EeXkEN1xh60uNZVo8mQ+d8wY7TlVqoipYBSBQJsfelNEkIeQs2eJ6tVzvNmHDPG+n3HjAivEr7uOdfR9+7L5Y4ECRP36sSnjV1/p268DbI5nVRUlJRENG8bzzc0lGjuWqFEj/+Z15522/8H06exwk5DAx4zmBBDFx2vrPv/c8X+6bp1+O4B15PZs3cr/m6wsHrddO6I9e4ieeII/N8B2/hs3ev/dCmGLCHKBefttrZCwWIj27vWunyVLPBN8Zcu6FnBGpWpV21g5OVpvw+Rk/fOuvZa9Hn/7zXGla2X5ct+FeGIi0V9/cT+//84/GlYh3asXj2dkV5+QQHTzzfy6cGGiV17Rzu3RR/XPffNNbduJE7XtGjTgY2fOeP99ChGBkSCXOJCxxpo12rq8PNa/ekPz5kCbNo519rHNAbYfnz7dOAqjK9LTbenr4uO14XR79tQ/r1cvtgVv3Fg/g8wNN2hjh8TFAV98weO44u672a48N5cTW+/ezfW5uRwKeNo04MYb9c8tV45jk5w7x7FnXn1V2yYnR//cqVO1dV99pa1bs4bj4xQqxElEhJhBBHmsUaeOti4uDqhVy/u+5sxhAfjgg8B77wE7d7LzUIMGbJ+8eDEH6tLboDTatLSyZg2HZDVi5EiOnW7FYuGN2CefdN1vXBxvqvbsyYK+Xj1OUtCtm2vnqMREYOBAfr1uHafOc2b2bJ5DgQLaY/v2AUeP8o+d0Q9Gr1769X/8wYmN7TGKE+9J/PgDBzhImhA96C3TA11EtRJCTp7kOB32j+SDBgV2zAEDtGoAZz29XklKYrNEV5w6xeZ9nobmdeaTT9ikr1QpfXVNejoHAVu6lNufPk3UoYP+fHv14jY33qh//Ndf3c/HSIc/dy5v2g4YQPT44/oqsmbNXPd98SJR9+62DeEaNcRuPMKA6MiF/+fCBaIvvuBNzmXLAj/e+fMs5OLjWV/esaNr5xZrUYro3LnAzWvkSNfjV6qkHd/eYsRZB/7779zmpZe0x9PSHB2sjNDTfaemEn36qaNFjsXCOvVq1fjHpndvomPHXPf9yivavuvU8e1/Z8Tmzfxd16tH9NhjREeOmNt/jCOCXAgcv/5KdN99RD17Es2ebdzu7FnH8KtvvGFspQGwpUpOjnnz3LXLtllJZNusdDX+hQu29tnZ+kHGkpN5E9XKyZOOoV/j4lgQe0JeHlvnWDeIixQhmjVLPxRw3brefX6jpBP//ONdP0bs3at1Drv6anO/wxhHBLkQGObN01qlvPuu5+fv28d22XoxrwFzzOdOnCBq0cLWZ716RP/+61lIgBMnbP1cvGgz7bMv9hY2VnJy2Nb8k094LG/Zt49/ILOzWbjrzS011bs+9SJQJiTwD48Z/O9/+vOcP9+c/gVDQS6bnYJ/vP66Npv58OFsyeEJZcrwxujVV2uPJSYaR+Lzhmefdcxms3490K8f0LWr6/OuvdYxv2mBAkDv3tp29nkrf/+dc1p27cqelv36AeXLez/nMmXYwiY5mTdo9axhmjf3rs+nn9ZuMvft65tVkR4nT+rXnzhhTv+CMXrSPdBFVuRRhFFeT1cZbPT44w+bI4+1PPOMOXMsVkw7P6WIjh41jk8CsK7bmYsXiZ56iqhECQ6ENXy4zXNy1ChtH88/b85n2LCBx7T2W64c0fbt3vfz448cSKtJE0744Wyf7w+LF2s/f4EC/H8WTAGiWhECwj33aG/ehg1962vbNvZMvOsuomnTzJujXnagIkVsQuzzz/UF+cKFno9x5oy+vj8pybwN24sX2at1wAAWwgcOmNOvmbzxhu0HuXhx1u8LpiGCXAgMhw45RuQrW9YxO1A48OWXWgE7dKjteE6Oow4dYMsab2KUfPON8cp+zx5zPsf69Y7hbQsUIPrpJ3P6NpOzZ3lT+cqVUM8k6jAS5G5c2QTBDSVKsM551SrOxt60qXsPyWBzzz3s7TluHHuL3n03cN99tuPx8cCPP3KmnT//BOrXBzp0cO+0ZE9ysn59ZiZ7dZrBc88Bp07Z3l+6BDzxBLBlizn9m0XBguwBKwSNMLvjhIhEKeD660M9C9d06sTFmV27gBde4AxGNWqw63zjxt7336YNC2znjEmffOLbfO05dIg9Z5cs0R776y/+ATX6IRFiAhHkQuxy8SJw8802d/v9+4Fff+VcmZ7GG7eSmMiCdvBg4Jdf2OrknXeAFi38m+OZM/zDsmeP/vFq1USICyLIhSCTn89BoBYt4hVs//6sngkFc+dqY6ZkZwMTJ3LMGG+pWFE/wJU/TJpkLMQTEoB33zV3PCEiETtyIbg88AAHhxo/Hvjf/9hW+8AB88fJzWWh7IqLF/Xrp0xxf26wOHhQv75DB4502L59cOcjhCUiyIXgsXMnR0u05/BhYPRo88YgAl55hTcZ09KA1q21emsr7dtzmjRntm8HXnvNvDn5g3OoYID3JN55B6hUKfjzEcISEeRC8Ni1S79+507zxvj0U2DoUM5TSQQsXMghavXIzNSP6w0As2aZNyd/aNqUf5gSEvh9cjLw/vusGxeE/xBBLgSPhg31Y3V762ruiq+/1tatWWP8Y9G6tf6c7F3zQ82rr3I886VLWaf/2GOhnpEQZoggF4JHZibw4Ye21SXAgrRfP/PGSEryrj41lRNjOPPEE+bNyQyKFweaNQuvHxghbBCrFSG49O0LtGvHq8vy5c23P+/fn9Up9lhtvI0YORKoUIEtTtLSgAEDgC5dzJ2XIAQQxV6fwaV+/fq0du3aoI8rxAgTJrBZ3vHjHFnxzTfZs1MQIhyl1Doiqu9cLytyIfq45x4ughAjiCAXYoedO9kJKCODLVnS0kI9I0EwBRHkQmwwYQJw//22JBivvAKsWGFeQCtBCCFitSJEP5cuAU895ZjJaN8+zm4UCThnYBIEJ0SQC9HP7t36acjWrQv+XLzh9GkOuZuczOnYXngByMsL9ayEMEQEuRD9VKigb7VSt27Qp+IVffuyg9OVKyzU33gDGDEi1LMSwhAR5EL0k5ICvPWWY6KIEiWAF18M3ZzccfYsJ7pwZuLE4M9FCHtks1OIDfr356z0s2ezh+mddwLp6aGelTFKAXFxWv24xRKa+QhhjQhyIXaoVYtLJFCwINC9O4fUtadv39DMRwhrRJALQrgybhw/NUydyuqhRx4BBg4M9ayEMERc9AVBECIEIxd92ewUBEGIcESQC4IgRDgiyAVBECIcEeSCIAgRjghyQRCECEcEuSAIQoQjglwQBCHCEUEuCIIQ4YggFwRBiHBEkAuCIEQ4IsgFQRAinJDEWlFKHQOwJ+gDB5YsAMdDPYkAEc2fDYjuzxfNnw2I7s+n99nKE1FR54YhEeTRiFJqrV4wm2ggmj8bEN2fL5o/GxDdn8+bzyaqFUEQhAhHBLkgCEKEI4LcPMaFegIBJJo/GxDdny+aPxsQ3Z/P488mOnJBEIQIR1bkgiAIEY4IchNRSo1QSm1TSm1USs1SSoVxmnbvUEp1V0ptUUrlK6WiwkpAKdVGKbVdKbVLKTU41PMxE6XUeKXUUaXU5lDPxWyUUmWVUr8opbb+d01GVSJTpVQBpdRqpdSf/32+V92dI4LcXBYCqElEtQHsAPB8iOdjJpsBdAGwLNQTMQOllOX/2rmb0DrKMIrj/6PGUqg7C0orKCLWojQiSKsbqYVGEUVBsIsi2IILFYQKLroodCt0VasIQjfBolhRqjbGL5qFXyhBlCiIIkZFQQlWBEv0uJgRS7i5d5JM+jLX84NA5uu9zxBymMz75gGeBG4DNgO7JG0uW1WrjgJjpYtYJfPAPtvXAFuBh4bsZ/cnsN32FmAUGJO0td8FCfIW2X7D9ny9+T6wsWQ9bbI9Y/vL0nW06EbgK9tf2z4DHAPuKlxTa2yfAn4tXcdqsP2j7U/q708DM8CGslW1x5Xf682R+qvvZGaCfPU8ALxeuohY1Abgu7O2ZxmiMPi/kHQ5cD3wQdlK2iXpfEnTwM/ApO2+93fBuSlreEh6E7ikx6H9tl+uz9lP9eff+LmsbaWa3NsQUY99WcLVIZLWAS8Cj9r+rXQ9bbL9FzBaz7O9JOla24vOdyTIl8j2jn7HJd0P3AHc6o6t7Rx0b0NmFrjsrO2NwA+FaoklkjRCFeLjto+Xrme12J6T9C7VfMeiQZ5XKy2SNAY8Dtxp+4/S9URfHwFXSbpC0oXAfcArhWuKBiQJeBaYsX2odD1tk7T+3xVvktYCO4Av+l2TIG/XYeAiYFLStKSnSxfUFkl3S5oFtgGvSpooXdNK1JPSDwMTVJNlz9v+vGxV7ZH0HPAecLWkWUl7StfUopuB3cD2+vdsWtLtpYtq0aXAO5I+pXrgmLR9ot8F+c/OiIiOyxN5RETHJcgjIjouQR4R0XEJ8oiIjkuQR0Q00HYjMkknJc1JOrFg/3jdzO2z+jNHBo2VII+IaOYo7TYie4JqGeVC48Am4DpgLbB30EAJ8oiIBno1IpN0Zf1k/bGkKUmbljDeW8DpHvtfqxtnGfiQBs33EuQREcv3DPCI7RuAx4AjbQ1cv1LZDZwcdG56rURELEPdtOsm4IWqawAAa+pj9wAHe1z2ve2dDT/iCHDK9tSgExPkERHLcx4wZ3t04YG6kdeym3lJOgCsBx5sWkhERCxR3Tr3G0n3QtXMS9KWlY4raS+wE9hl++9G16TXSkTEYHUjsluAi4GfgAPA28BTVI2uRoBjtnu9Uuk13hTV6pR1wC/AHtsTkuaBb/lvIvT4oDET5BERHZdXKxERHZcgj4jouAR5RETHJcgjIjouQR4R0XEJ8oiIjkuQR0R0XII8IqLj/gHrMhCFqNFOMgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from numpy import random,sqrt,pi,count_nonzero\n", "\n", "# Constants\n", "Z = 79 # atomic charge of nucleus\n", "e = 1.602e-19 # unit of elementary charge\n", "E = 7.7e6*e # kinetic energy\n", "epsilon0 = 8.854e-12 # permittivity of free space\n", "a0 = 5.292e-11 # Bohr radius\n", "cutoff = Z*e*e/(2*pi*epsilon0*E) # scatter if b less than this\n", "sigma = 0.01*a0 # standard deviation\n", "N = 1000000 # number of incident particles\n", "\n", "rng = random.default_rng() \n", "\n", "x = rng.normal(0, sigma, N)\n", "y = rng.normal(0, sigma, N)\n", "r = sqrt(x*x+y*y) # calculate r from x and y\n", "\n", "particle_scatter = count_nonzero(r < cutoff)\n", "scatter_percent = 100*particle_scatter/N\n", "\n", "print(\"{} particles were reflected out of {} ({:.3f}%).\".format(particle_scatter, N, scatter_percent))\n", "\n", "from numpy import where\n", "import matplotlib.pyplot as plt\n", "\n", "plt.figure(figsize=(6,6))\n", "col = where(r