{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise 1: Uniqueness" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The problem\n", "\n", "In [the first tutorial](../tutorials/1_sphere_scatterer_null_field.ipynb), we looked at two formulations for a scattering problem with a Neumann boundary condition.\n", "\n", "In this exercise, we will investigate the uniqueness so solutions to the boundary integral formulations for this problem." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The formulation\n", "\n", "In this exercise we will use the null field approach (as in [the first tutorial](../tutorials/1_sphere_scatterer_null_field.ipynb). This uses the following representation formula and boundary integral equation.\n", "\n", "### Representation formula\n", "\n", "$$\n", "p_\\text{total} = \\mathcal{D}p_\\text{total} + p_\\text{inc}\n", "$$\n", "\n", "where $\\mathcal{D}$ is the double layer potential operator.\n", "\n", "### Boundary integral equation\n", "\n", "$$\n", "(\\mathsf{D}-\\tfrac{1}{2}\\mathsf{I})p_\\text{total} = -p_\\text{inc},\n", "$$\n", "\n", "where $\\mathsf{D}$ is the double layer boundary operator, and $\\mathsf{I}$ is the identity operator." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Finding a resonance\n", "\n", "The code below plots the condition number of $\\mathsf{D}-\\tfrac{1}{2}\\mathsf{I}$ for 30 values of $k$ between 2.5 and 3.5. There is a sharp increase in the condition number near 3.2.\n", "\n", "Adjust the limits use in `np.linspace` to approximate the value of $k$ for this spike to 4 or 5 decimal places. (For example, you might start be reducing the search to between 3.0 and 3.3, so `np.linspace(3.0, 3.3, 30)`.)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEJCAYAAAB7UTvrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZxcdZnv8c9TvSdd3dm6OyEhZCEJe0AiIouDLDOjwxV0kAHFQWRkrngFR0dRr3f06swdnTvjOPNyUBFx4jJwGXYcQTESkEUwkLAGCCSQlaSz0Et673ruH+dUpdLp7lS665zqqvq+X69+VZ2l6jwnSz/nt5u7IyIiApAodAAiIjJxKCmIiEiGkoKIiGQoKYiISIaSgoiIZFQWOoDxmDFjhs+bN6/QYYiIFJWnnnpqp7s3DXesqJPCvHnzWLVqVaHDEBEpKmb2xkjHVH0kIiIZSgoiIpKhpCAiIhlKCiIikqGkICIiGUoKIiKSEVlSMLObzGyHmT2ftW+amT1gZuvC16lZx75oZq+a2ctm9kdRxSUiIiOLsqTw78AfD9n3BWCFuy8CVoTbmNkxwCXAseFnrjezighjE5ERPPXGbl7Y2lboMKRAIksK7v4wsHvI7guA5eH75cCFWftvcfded98AvAqcElVsIjKyv7n7Bb5x30uFDkMKJO42hRZ33wYQvjaH+2cDm7LO2xzuO4CZXWVmq8xsVWtra6TBipSjPXv72N7eU+gwpEAmSkOzDbNv2CXh3P0Gd1/m7suamoadukNExqGtu58dHb2FDkMKJO6ksN3MZgGErzvC/ZuBw7POmwNsjTk2kbLXP5hib98gb3X10zswWOhwpADiTgr3AJeH7y8H7s7af4mZ1ZjZfGAR8GTMsYmUvfbu/sz7He0qLZSjKLuk3gw8Diwxs81mdiXwDeA8M1sHnBdu4+4vALcCLwL3A590dz2miMSsLTspqAqpLEU2dba7XzrCoXNGOP/vgL+LKh4RObjspNDaocbmcjRRGppFZAJQSUGUFEQko01tCmVPSUFEMtINzdUVCXao+qgsKSmISEa6pDB/xmRVH5UpJQURyWjr7qeuqoI5U+tUfVSmlBREJKOtu5/GuiqaG2pUUihTSgoikpFOCk3JWnbt7WVgMFXokCRmSgoiktHW3U9DXSXNyRrcYdfevkKHJDFTUhCRjLbugaD6KFkDqFtqOVJSEJGM9u5+GuqqaG6oBVC31DKkpCAiGZmG5rCksF0lhbKjpCAiAAwMpujsDaqPZtSH1UcqKZQdJQURAaC9ZwCAxroqqisTTJtcrW6pZUhJQUSAfaOZG+uqAGhO1qihuQwpKYgIcGBSaErWaPrsMqSkICLAcCWFWlUflSElBREBDkwKLQ01tHb0kkp5IcOSmCkpiAgwfJvCQMrZ06VRzeVESUFEgH1rKTSkk0JmAJuqkMqJkoKIAEFJoaYyQW1VBcC+qS6UFMqKkoKIANDW1Z+pOoKgoRlgR7t6IJUTJQURAfZNcZHW3KCSQjlSUhAR4MCkUFtVQbK2klYlhbKipCAiwIFJAYJ2he2qPiorSgoiAoyUFDSArdwoKYgIsG8thWzBWs0qKZQTJQURYTDldITTZmdLT4rnrlHN5UJJQUQyA9eGqz7qHUhlptWW0qekICIHTHGRlu6WqtlSy4eSgoiMnBQyA9jU2FwulBREZF9SmDR8SUE9kMqHkoKIjFJS0FrN5UZJQURo7xk+KdTXVFJXVaHqozJSkKRgZn9lZi+Y2fNmdrOZ1ZrZNDN7wMzWha9TCxGbSDkaqaRgZuFYBSWFchF7UjCz2cA1wDJ3Pw6oAC4BvgCscPdFwIpwW0Ri0NbdT3XWtNnZmpMawFZOClV9VAnUmVklMAnYClwALA+PLwcuLFBsImWnfZgpLtI01UV5iT0puPsW4B+BjcA2oM3dfwW0uPu28JxtQPNwnzezq8xslZmtam1tjStskZI23LxHaU3hqGYpD4WoPppKUCqYDxwGTDazy3L9vLvf4O7L3H1ZU1NTVGGKlJXRkkJzQw2dvQN09WlUczkoRPXRucAGd291937gDuA0YLuZzQIIX3cUIDaRsjRqUtAAtrJSiKSwETjVzCaZmQHnAGuBe4DLw3MuB+4uQGwiZWm0pNCiAWxlpTLuC7r7E2Z2G/A0MACsBm4A6oFbzexKgsTxwbhjEylXQ9dnzpYpKagHUlmIPSkAuPtXgK8M2d1LUGoQkRilwmmzh66lkJYZ1azqo7KgEc0iZa6jZwD3AweupU2ZVEV1RULVR2VCSUGkzI00mjnNzIJuqao+KgtKCiJl7mBJAYKxCq0qKZQFJQWRMpdLUmjWALayoaQgUuZySgoNqj4qF0oKImUut5JCLXu6+ukbSMUVlhSIkoJImcu1+gigtVNVSKVu1KRgZhVm9n/jCkZE4tfW3U91RYLaqpF/HaSX5dzeriqkUjdqUnD3QeDkcDoKESlBbd39NNRVMdp/c81/VD5yGdG8GrjbzP4T2Jve6e53RBaViMQmWEth9F8F6ZJCqxqbS14uSWEasAs4O2ufE8xuKiJFbrTJ8NKmT64hYZoUrxwcNCm4+xVxBCIihdHW3c/0+upRz6lIGDPqNVahHBy095GZLTazFWb2fLh9gpl9OfrQRCQOuZQUQGMVykUuXVJ/AHwR6Adw92eBS6IMSkTik3NS0FrNZSGXpDDJ3Z8csk/r8omUgFTKae/JNSnUKCmUgVySwk4zW0jQuIyZXQRsizQqEYlFR+/o02Zna07WsKuzl8GUxxCZFEouvY8+SbAy2lFmtgXYAHw40qhEJBbt4WjmkRbYydbUUEvKYVdnL80NtVGHJgWSS++j9cC5ZjYZSLh7R/RhiUgccpniIi2zAluHkkIpy6X30XQz+1fgt8BKM/sXM5sefWgiErWxJQX1QCplubQp3AK0An8KXBS+/39RBiUi8TikpNCgqS7KQU4jmt3961nbf2tmF0YVkIjE51CSQlN9elI8JYVSlktJ4UEzu8TMEuHPxcB/RR2YiETvUJJCdWWCaZOrVX1U4kYsKZhZB0E3VAM+A/w0PJQAOoGvRB6diESqrbufyoQxqboip/M1VqH0jZgU3D0ZZyAiEr/0aOZcZ8dvUlIoebm0KWBmJwDzss/X1NkixS/XKS7SmpO1vLZjZ4QRSaEdNCmY2U3ACcALQHqBVk2dLVIC2sMFdnLV3FBDa2cv7p5z6UKKSy4lhVPd/ZjIIxGR2LV19zN10ujTZmdrTtbQP+js6epn2uTcPyfFI5feR4+bmZKCSAkaS/URaABbKculpLCcIDG8CfQS9EZydz8h0shEJHKHnBTCZTl3tPdy1MyoopJCyiUp3AR8BHiOfW0KIlLkUikP12c+lJLCvvmPpDTlkhQ2uvs9kUciIrHq7BsgleO02WmqPip9uSSFl8zsP4B7CaqPAHVJFSl2bV25j2ZOq6uuIFlTqfmPSlguSaGOIBn8YdY+dUkVKXJth7CWQrYmrdVc0nJZT+GKfF/UzKYANwLHESSYjwEvE8y+Og94HbjY3ffk+9oiEmg/hHmPsrUka1VSKGG5DF77EeFSnNnc/WPjuO6/APe7+0VmVg1MAr4ErHD3b5jZF4AvANeN4xoiMopDmQwvW3NDDas3vhVFSDIB5FJ99POs97XA+4GtY72gmTUA7wI+CuDufUCfmV0AnBWethxYiZKCSGQySWHSISaFZFB9pFHNpSmX6qPbs7fN7Gbg1+O45gKChXp+ZGZLgaeAa4EWd98WXnObmTUP92Ezuwq4CmDu3LnjCEOkvI25pJCspac/RUfvAA21h/ZZmfhyGdE81CJgPL+NK4G3Ad9195OAvQRVRTlx9xvcfZm7L2tqahpHGCLlra27n4qEMTnHabPTsgewSenJZY3mDjNrT78SdE0dT7XOZmCzuz8Rbt9GkCS2m9ms8JqzgB3juIaIHMShTpud1qS1mktaLtVHeV1Xwd3fNLNNZrbE3V8GzgFeDH8uB74Rvt6dz+uKyP4OdYqLtPQAtlaNai5Jua6nMBs4gv3XU3h4HNf9FPCzsOfReuAKglLLrWZ2JbAR+OA4vl9EDqLtEKfNTlP1UWnLpUvqN4E/I3iSHwx3OzDmpODua4Blwxw6Z6zfKSKHpr27n8ZDmDY7LVlTSW1VQtVHJSqXksKFwBJ312OBSAlp6+5n7vTJh/w5M6M5WatJ8UpULr2P1gPqdyZSYoI2hZxqkA/QnKxR9VGJyuVfRBewxsxWsP+EeNdEFpWIRMrdae8ZGFNDMwTtCi+/2ZHnqGQiyCUp3BP+iEiJ6OwdYDDlY08KyVp++8rOPEclE0EuXVKXxxGIiMRnrKOZ05obaujoHaC7b5C6Qxz8JhPbWEY0i0iRG3dS0GI7JUtJQaQMtXcPAIe+lkKaluUsXUoKImUoH9VHoAFspSiXwWuLgc9x4IjmsyOMS0QiNNYFdtJUfVS6cul99J/A94AfsG9Es4gUsfGWFKZOqqKqwlR9VIJySQoD7v7dyCMRkdikp82urxnb4DUzo6leA9hKUS5tCvea2dVmNsvMpqV/Io9MRCLT1t1PQ23luFZOa2qoVfVRCcrlMeHy8PVzWfucYAU1ESlCY502O1tzsoZNu7vyFJFMFLkMXpsfRyAiEp98JYWn3tiTp4hkosil91EV8AngXeGulcD33b0/wrhEJEJjXUshW3Oylt17++gbSFFdqd7tpSKXv8nvAicD14c/J4f7RKRItechKbSEYxV2dqqxuZTk0qbwdndfmrX9GzN7JqqARCR6+ag+amkMxipseaubw6bU5SMsmQByKSkMmtnC9IaZLUDjFUSKlrvnJSksaq4HYN32znyEJRNELiWFzwEPmtl6wAhGNl8RaVQiEpmuvkEGxjFtdtrsKXVMrq7gle1aV6GU5NL7aIWZLQKWECSFl7Q0p0jxGu9o5jQzY/HMpBbbKTEjJgUzO9vdf2NmHxhyaKGZ4e53RBybiEQgX0kBYElLkgde3D7u75GJY7SSwh8AvwH+2zDHHFBSEClC+UwKi1uS3PL7Tezs7GVGfc24v08Kb8Sk4O5fCd9+zd03ZB8zMw1oEylSeS0pzEwC8MqbHcw4UkmhFOTS++j2Yfbdlu9ARCQe+S4pALysxuaSMVqbwlHAsUDjkHaFBqA26sBEJBrptRTGO3gNYEZ9NdMmV6sHUgkZrU1hCXA+MIX92xU6gI9HGZSIRKetux8zSI5x2uxsZsai5npe0ViFkjFam8LdwN1m9k53fzzGmEQkQsG02VUkEmOfNjvbkplJ7nx6C+4+rqm4ZWIYrfro8+7+D8CHzOzSocfd/ZpIIxORSORjNHO2xS1JOnoH2NbWo+kuSsBo5ce14euqOAIRkXjkOymkeyC9vL1DSaEEjFZ9dG/4ujy+cEQkankvKTTv65b67iXNefteKYzRqo/uJRikNix3f18kEYlIpNq6+zmsMX9P9I2TqpjZUKtuqSVitOqjfwxfPwDMBH4abl8KvB5hTCISoXyspTDU4plJdUstEaNVHz0EYGZfd/d3ZR2618wejjwyEcm7fE2bPdSSlnp+/PguBlNORZ56NUlh5DKiuSlcQwHITHHRNN4Lm1mFma02s5+H29PM7AEzWxe+Th3vNURkf939g/QPjn/a7KEWtyTpHUixcXdXXr9X4pdLUvgrYKWZrTSzlcCDwKfzcO1r2dfDCeALwAp3XwSsCLdFJI/yOcVFtswcSKpCKnoHTQrufj+wiOCX+LXAEnf/5XguamZzgD8BbszafQGQ7um0HLhwPNcQkQNFlRSODFdhe0VrKxS9XMe5nwzMC89fGq6n8ONxXPfbwOeBZNa+FnffBuDu28xs2L5tZnYVcBXA3LlzxxGCSPlp64omKUyqrmTutEnqgVQCDlpSMLOfEPREOgN4e/izbKwXNLPzgR3u/tRYPu/uN7j7Mndf1tQ07qYNkbISVUkBgnYFVR8Vv1xKCsuAY9x9xDELh+h04H1m9l6C2VYbzOynwHYzmxWWEmYBO/J0PREJRZkUlsysZ+XLO+gbSFFdmUtzpUxEufzNPU8wTiEv3P2L7j7H3ecBlwC/cffLgHuAy8PTLgfuztc1RSQQdUlhIOVs2Lk3798t8cmlpDADeNHMngR60zsjGNH8DeBWM7sS2Ah8MM/fL1L22tPTZteOf9rsobLnQEq/l+KTy7+Mr0Z1cXdfCawM3+8CzonqWiISlBSSNZV5mzY724IZ9VQmLOiBtDTvXy8xOWhScPeHzKyFoIEZ4El3V32/SBFq6+6ncVL+q44AqisTzJ8xWT2QilwuvY8uBp4kqM65GHjCzC6KOjARyb8oprjIpjmQil8u1Uf/E3h7unRgZk3Ar4HbogxMRPIv8qTQnOQXz22ju2+QuuqKyK4j0cml91FiSHXRrhw/JyITTNRJYcnMetzh1R1as7lY5VJSuN/MfgncHG7/GXBfdCGJSFTaugeiLSm07OuBdPycxsiuI9HJpaH5c2b2AYIRzQbc4O53Rh6ZiOSVu0eylkK2I6ZPproyoXaFIjbaymtHEsxH9Ki73wHcEe5/l5ktdPfX4gpSRMavpz9F32Aq0pJCRcJY1FzPy5oYr2iN1jbwbWC4v9mu8JiIFJEoRzNnW6I5kIraaElhnrs/O3Snu68imDFVRIpIXElh8cwk29p6MteT4jJaUqgd5Vj+Vv0WkVjEWVIAWKfSQlEaLSn83sw+PnRnODfRmKa9FpHCibOkAGhkc5EarffRp4E7zezD7EsCy4Bq4P1RByYi+RVXUjissZb6mkqtwlakRkwK7r4dOM3M3g0cF+7+L3f/TSyRiUhexZUUzIxFLfW8sl0D2IpRLuMUHgQejCEWEYlQOikka6NNChC0Kzzw4vbIryP5p+kqRMpEe3c/ydpKKiKYNnuoxS1Jdu3tY2dn78FPlglFSUGkTEQ971G29CI7alcoPkoKImUizqSQPQeSFBclBZEyEWdSmFFfzbTJ1RrZXISUFETKRJxJwcxY3KI5kIqRkoJImWjr7qchhp5HacEcSJ24e2zXlPFTUhApE1GuzzycRS1JOnsH2NrWE9s1ZfyUFETKQE//IH0D0U6bPZR6IBUnJQWRMtAeDlyLcoGdoRY3h0lBjc1FRUlBpAzENcVFtsZJVcxsqFW31CKjpCBSBgqRFCCYMVUlheKipCBSBtKT082or471ukta6lm3vZPBlHogFQslBZESN5hybvzteo49rIFjZjXEeu3FLUl6B1Js3N0V63Vl7JQURErcL194k/U79/KJsxZiFv1keNnSPZA0iK14KCmIlDB35/qVrzJ/xmTec9ys2K9/ZHM9ZuqBVEyUFERK2G/X7eT5Le385bsWxDJl9lCTqis5fOok9UAqIkoKIiXs+pWv0tJQw/vfNrtgMSxuSWoAWxFRUhApUU9v3MPv1u/m42cuoKayomBxLJlZz4ade+kbSBUsBsmdkoJIibr+wdeYMqmKS0+ZW9A4FrckGUg5G3buLWgckpvYk4KZHW5mD5rZWjN7wcyuDfdPM7MHzGxd+Do17thESsXLb3bw67Xbufyd85hcc9Cl2COV6YGkdoWiUIiSwgDwWXc/GjgV+KSZHQN8AVjh7ouAFeG2iIzB9x56jUnVFXz0tHmFDoUFM+qpTJjaFYpE7EnB3be5+9Ph+w5gLTAbuABYHp62HLgw7thESsGm3V3c88xWLj1lLlMnxzuCeTjVlQkWNtXz4Ms7GBhUu8JEV9A2BTObB5wEPAG0uPs2CBIH0DzCZ64ys1Vmtqq1tTWuUEWKxg0Prydh8Bdnzi90KBlXv3shL2xt5/sPry90KHIQBUsKZlYP3A582t3bc/2cu9/g7svcfVlTU1N0AYoUodaOXm5dtYkPnDSHWY11hQ4n431LD+NPjp/Ft3/9Ci9uzfm/uxRAQZKCmVURJISfufsd4e7tZjYrPD4L2FGI2ESK2U2PbqBvMMVf/sGCQoeyHzPj6xcex5RJ1Xzm1jX0DgwWOiQZQSF6HxnwQ2Ctu38r69A9wOXh+8uBu+OOTaSYtff089PH3+C9x81iQVN9ocM5wLTJ1XzzT4/npTc7+OcH1hU6HBlBIUoKpwMfAc42szXhz3uBbwDnmdk64LxwW0Ry9JPH36Cjd4BPnLWw0KGM6OyjWrjk7Ydzw8Ovser13YUOR4YRewdmd38EGGkSlnPijEWkVPT0D/KjRzfwrsVNHDe7sdDhjOrL5x/DI6/u5LP/+Qy/uObMgo+jkP1pRLNICbh11SZ2dvZx9QQuJaTV11Tyjx9cysbdXfz9fWsLHY4MoaQgUuT6B1N8/6H1vG3uFN4xf1qhw8nJqQum8xdnzOenv9vIQ6+oa/lEoqQgUuTufWYrW97q5uqzjox9EZ3x+OwfLmFRcz2fv+0Z2rr6Cx2OhJQURIpYKuV8d+VrLGlJcvZRw473nLBqqyr41sUnsquzj7+55/lChyMhJQWRIvbrtdtZt6OTT5y1kEQBFtEZr+PnNPKpsxdx95qt/Nez2wodjqCkIFK0Xt+5l7/7xVoOn1bH+SfEv9Rmvlz97oUsndPIl+96jh0dPYUOp+wpKYgUocdf28WF1z9Ke3c//3zxiVRWFO9/5aqKBP908Yl09Q3yxdufw90LHVJZK95/SSJl6pYnN/KRHz7BjPoa7vrk6SybVxw9jkZzZHM91/3xUax4aQe3rtpU6HDKmkaNiBSJwZTz979Yy42PBIPUvvOhk2iorSp0WHnz0dPm8cCL2/navS8ybXIN5x3TUuiQypJKCiJFoKOnn4//eBU3PrKBj542j5suX1ZSCQEgkTD+6eKlHD5tEh//8SquuXk1u/f2FTqssqOkIDLBbdrdxUXffZyHXmnlby88jq++79iibkMYzWFT6rjnf5zBX527mPue38Z533qIe5/ZqnaGGJXmvyyRErHq9d1c+G+Psq2tm+VXnMJlpx5R6JAiV12Z4NpzF/HzT53JnKl1fOrm1Vz1k6fY3q6eSXFQUhCZoO54ejMf+sETNNRVcecnT+eMRTMKHVKslsxMcvsnTuNL7z2Kh19p5dxvPcStqzap1BAxJQWRCaZ/MMU/3P8Sn7n1GU4+Yip3Xn0aCyfg+ghxqKxIcNW7FnLftWdy9MwGPn/bs/z5TU+yeU9XoUMrWVbMWXfZsmW+atWqQochkhdvdfXxH09u5MePvcGb7T1cespcvnbBsVSVaPvBoUqlnJ898QZ/f99LGHDde47isnccUZQjuQvNzJ5y92XDHlNSECms9a2d/OjR17ntqc109w9y+pHT+YszFnDWkqaimuAuLpt2d/GlO5/jt+t2cmRzPR9422wuOHE2s6dMnDWpJzolBZEJxt15/LVd/PCRDax4aQfVFQkuOPEwPnbGfI6e1VDo8CY8d+fuNVv5ye/e4Kk39gDwjvnTeP9Js3nP8bNorCut7rr5pqQgMkH0Dgxyz5qt/PCRDbz0ZgfTJ1dz2alHcNmpR9CUrCl0eEVp464u7lqzhbtWb2H9zr1UVyY49+hmLjxxNmctaaa6UtVvQykpiBRQV98Aq17fw6Ov7uT2p7ews7OXJS1JrjxjPu878TBqqyoKHWJJcHee3dzGnau3cO8zW9m1t48pk6o4/4RZvG/pbJYe3khNpf6sQUlBJFZ9AynWbHqLx17byWOv7mL1pj30DzqVCePMRTO48owFnH7kdLUXRKh/MMUj63Zy5+ot/OrFN+npT1FdkeDoWUlOmDOFE+Y0cuLhU1jQVE9FGTZUKymIRGgw5Ty/pY3HXtvFY6/t5Pev76anP4UZHD+7kXcunM5pC2ew7IipWqS+ADp7B/jtK62s2fwWz2x6i+e3tNPZOwDA5OoKjpsdJIh0spgzta7kE7aSwhDuzmAquG8H3MFx0n8UB2yHn0mfm9453Dnp99nfk70z+1iwve/P34eck441CkP/0dt+x4aci2X22ZAPDD2W/l4Lv8f2nRhu77t+9jmZ78jaTn9fwg6MNy6DKWdXZy9vtvewvT143dHew5ttPWzv6GVHew+b93RnfsksaUmGSWA675g/ncZJavCcaFIpZ/3OTtZsauPZzW/xzOY21m5tp28wBUBNZYLmhhqak7U0J2uCn4ba/V+TNUydVE0iYaRSTsqdlBO+Br9fUr7vd03mffq8VPC7YzA8Pzhv3+c9/K7sz6Y8+PeYfj+9vnrMnRKUFIZYs+ktLvy3RyOISKJmBolhEkr6fTqBpJNKImH7JRfCcxK2/7mJYV47ewZo7ezNPECkJQyakjXMbKiluaGWwxprWTZvGqcumK7G4iLVN5DipTfbeWZzG5t2d7G9vYcd7b3s6OhhR0cvHT0DhQ7xAOefMIvvfOhtY/rsaEmhLMuysxpr+cx5i7OeWrOecA94Ug3P2e9pdujT8L79Bz5R237XSX/X/tv7DD3ngBPywYduHlhayT51aKlmvxJNVuloX0nLh3wuqxSV9V1BSWq4ktqBx1LhRnpfKuv9vuuln7Cyvj98qkpvB7/fnVTmqSz72L7XlDuTqiuZ2VBLS2MtLckaZjbW0tJQy4z6mrKshy5l1ZWJsPpoyrDHu/sGMwkinSx27+3LPHBUJMKHjYSRMKMi/F1QEW4nEvseRiqyHjwSifQDimUeVNLfkXl4yfqsGVSE+6ZNro7kz6Isk0JLQy3XnLOo0GGISJGoq67giOmTOWL65EKHEjl14BURkQwlBRERyVBSEBGRDCUFERHJUFIQEZEMJQUREclQUhARkQwlBRERySjqaS7MrBV4o9BxjMEMYGehg4iZ7rk8lNs9F+v9HuHuTcMdKOqkUKzMbNVI846UKt1zeSi3ey7F+1X1kYiIZCgpiIhIhpJCYdxQ6AAKQPdcHsrtnkvuftWmICIiGSopiIhIhpKCiIhkKClExMwON7MHzWytmb1gZteOcN5ZZrYmPOehuOPMp1zu2cwazexeM3smPOeKQsSaD2ZWa2ZPZt3L/x7mHDOzfzWzV83sWTMb2/qJE0SO9/zh8F6fNbPHzGxpIWLNl1zuOevct5vZoJldFGeMeeXhotH6ye8PMAt4W/g+CbwCHDPknCnAi8DccLu50HHHcM9fAr4Zvm8CdgPVhY59jPdrQH34vhb7URsAAAaASURBVAp4Ajh1yDnvBe4Lzz0VeKLQccdwz6cBU8P37ymHew6PVQC/AX4BXFTouMf6o5JCRNx9m7s/Hb7vANYCs4ec9iHgDnffGJ63I94o8yvHe3YgacFC1/UESWHirYqeAw90hptV4c/QnhsXAD8Oz/0dMMXMZsUZZz7lcs/u/pi77wk3fwfMiTHEvMvx7xngU8DtQFH/P1ZSiIGZzQNOInjCyLYYmGpmK83sKTP787hji8oo9/wd4GhgK/AccK27p2INLo/MrMLM1hD8InjA3Yfe72xgU9b2Zg5MlEUlh3vOdiVBSamoHeyezWw28H7ge4WIL5+UFCJmZvUETw+fdvf2IYcrgZOBPwH+CPhfZrY45hDz7iD3/EfAGuAw4ETgO2bWEHOIeePug+5+IsHT8ClmdtyQU2y4j0UfWXRyuGcAzOzdBEnhujjji0IO9/xt4Dp3H4w/uvxSUoiQmVUR/HL8mbvfMcwpm4H73X2vu+8EHgaKvVHuYPd8BUGVmbv7q8AG4Kg4Y4yCu78FrAT+eMihzcDhWdtzCEpJRW+Ue8bMTgBuBC5w910xhxaZUe55GXCLmb0OXARcb2YXxhtdfigpRCSsM/8hsNbdvzXCaXcDZ5pZpZlNAt5BUA9flHK8543AOeH5LcASYH08EeaXmTWZ2ZTwfR1wLvDSkNPuAf487IV0KtDm7ttiDjVvcrlnM5sL3AF8xN1fiT/K/Mrlnt19vrvPc/d5wG3A1e5+V+zB5kFloQMoYacDHwGeC+siIeh5MxfA3b/n7mvN7H7gWSAF3Ojuzxck2vw46D0DXwf+3cyeI6hauS4sJRWjWcByM6sgeMC61d1/bmb/HTL3+wuCHkivAl0EJaVilss9/w0wneBpGWDAi3sm0VzuuWRomgsREclQ9ZGIiGQoKYiISIaSgoiIZCgpiIhIhpKCiIhkKCmIiEiGkoKIiGQoKUjRMbN/NrNPZ23/0sxuzNr+JzP7TGGiy52ZzTOzvA5WNLM6M3soHGiFmZ1rZj/JOl5tZg+bmQauyrCUFKQYPUYwZz9mlgBmAMdmHT8NeLQAccUqnDpj6P/hjxHMLZWemG0psDp90N37gBXAn8UTpRQbJQUpRo8SJgWCZPA80GFmU82shmBq7tVmdlc4JfkLZnZV+sNm9k0zuzpr+6tm9lkzuyxcYWuNmX0/62l7ngWryf0g/K5fhU/k+z3pm9lfh981z8xeMrMbzex5M/tZ+MT+qJmtM7NTsu6l0syWW7BK2W3hHFjkEMv1wNPsP9kewIcJ5tRKWxr+WdSY2b+b2f8B7grPEzmAkoIUHXffCgyEE6+dBjxOsG7DOwlmq3w2fCL+mLufHO67xsymh19xC/s/KV8MrAr3nR5OkTzI/r84FwH/5u7HAm8Bf3qQMI8E/gU4gWAW2A8BZwB/TTAfVNoS4AZ3PwFoB642s6MPEssSgoV7TnL3N9I7zawaWODur2edu5RgDYBfAr929y8RJNG3HyR+KVOqV5RilS4tnAZ8i2DhmtOANoLqJQgSwfvD94cT/GLf5e6rzazZzA4jWBJ0D3A8wdoWvw8ncatj/xW0Nrh7epK/p4B5wCOjxLfB3Z8DMLMXgBXu7uFEgPOyztvk7umqrp8C1wA9B4nljXAVt6FmECQswutWhde6GfhLd38cgrUBzKzPzJLhCnkiGUoKUqzS7QrHEzz5bgI+S/C0fZOZnUUwxfE73b3LzFYCtVmfv41g3vuZBCUHA5a7+xdHuF5v1vtBgl/UA+xf2q4d4fxU1naK/f/fDZ2R0nOIZe8I+7uHxHAM8HtgWhhzthqC5COyH1UfSbF6FDgf2B2uirUbmEJQhfQ40AjsCRPCUcCpQz5/C3AJQWK4jaDx9SIzawYws2lmdsRBYtgONJvZ9LAt4/wx3MdcM3tn+P5SgtLHWGIhXBe5wszSiWEpQfK8BPiRBetXEFajtbp7/xjilRKnpCDF6jmC6pLfDdnXFq7PcD9BI+6zBGs47Ffd4u4vAElgi7tvc/cXgS8Dvwo/8wDBPPojCn+pfo2gPePnHLjATi7WApeH15wGfHcssWT5FUHbBQRJ4flwoZvrgFvDKqV3E6zzIHIAracgUkLM7CTgM+7+kVHOuQP4oru/HF9kUixUUhApIe6+Gngw3YV1qLCH0l1KCDISlRRERCRDJQUREclQUhARkQwlBRERyVBSEBGRDCUFERHJUFIQEZGM/w+TUUOAFsBd/AAAAABJRU5ErkJggg==\n", "text/plain": [ "<Figure size 432x288 with 1 Axes>" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "import bempp.api\n", "from bempp.api.operators.boundary import helmholtz, sparse\n", "from bempp.api.operators.potential import helmholtz as helmholtz_potential\n", "from bempp.api.linalg import gmres\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "grid = bempp.api.shapes.regular_sphere(3)\n", "\n", "space = bempp.api.function_space(grid, \"DP\", 0)\n", "\n", "identity = sparse.identity(space, space, space)\n", "\n", "x_data = []\n", "y_data = []\n", "for k in np.linspace(2.5, 3.5, 30):\n", " double_layer = helmholtz.double_layer(space, space, space, k)\n", " x_data.append(k)\n", " y_data.append(np.linalg.cond((double_layer - 0.5 * identity).weak_form().to_dense()))\n", " \n", "plt.plot(x_data, y_data)\n", "plt.xlabel(\"Wavenumber ($k$)\")\n", "plt.ylabel(\"Condition number\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The effect on the solution\n", "\n", "The code below has been copied from [the first tutorial](../tutorials/1_sphere_scatterer_null_field.ipynb) and the wavenumber has been changed to 3. The solution plot looks like a reasonable soluition.\n", "\n", "Change the value of the wavenumber to the resonance that you found above. Does the solution still look reasonable?" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjYAAAHwCAYAAAC17yUBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOy9e9Rmx1Xe+dT3dbdaUkuRdbFkWb4AsgHHEwzRsrkYohBgbOOJgwkZmwyJSUAhsQcyuTCekLEXCSvLk1mTCzFgFPAiziQGVrDBIc4YGCYxhhB8wSayFRNHAUuWLNGWZamtS3frq/njnHrf+vZbe9euc857f35r9Trfe651rl311LN3hRgjCCGEEEJ2gYN1F4AQQgghZCpYsSGEEELIzsCKDSGEEEJ2BlZsCCGEELIzsGJDCCGEkJ2BFRtCCCGE7Ays2BDiJIRwOoQQQwg3rbssNUIIl4QQzoUQbpxof98TQviVKfblONanQwgvXsWxCCG7Bys2ZKvp//NO/45CCI9lv/9sZduXhBA+saqyrpIY4xMxxjMxxntbtw0hfEkI4eIyyrXLhBCuCyH86xDC50MI/y2E8KfXXSZC9pET6y4AIWOIMZ5Jf4cQfg/Ad8UYV6IsECK4HcBnATwVwIsA/HwI4SMxxv+y3mIRsl9QsSE7TQjh0hDCj4QQ7gsh3BNC+D9DCCdDCNcAeCeAL8wUnmtCCF8TQviPIYTPhRDuDSH8wxBCtQEQQnhpCOH92e/3hRDem/3+QAjhJf3fb+hb9I+EEO4IIXxzP//yvhw3Z9s9vVehntL//pYQwu+EEB4KIfxaCOF5SnmOdZuFEH46hPCPQgjv6Y/76yGEZymn814Ah9l1+fL5bsMP98f+ryGEb8iOd3UI4W19N9LdIYQ3hhCK35cQwptCCG8PIfxcX5b3hxD+sLKuej+yc/zuvjyfDSH8Q7H9XwohfDyE8GAI4d+EEJ6uHOcnhPr3ZAjh9cr1KW3/FAD/A4D/Pcb4+RjjrwJ4DwBTNSSETA8rNmTX+UEAfwTAfwfgjwK4FcD3xxg/A+BbANzVd9mc6eddAPA6AFcD+Fp0/1l9l+M47wPwR0IIV4YQTgP4IgDP7f/zvQLA8wH8er/uxwF8NYA/BOD/APDTIYRrY4yfB/AuAK/O9vsqAO+JMX42hPCVAH4UwHcCuAbAP0enCniV128H8L/153Zff21KfB2AJ7Pr8tvZ/A/0x34zgJ/ItvkXAD4H4AsBvBDAnwLwHUZZvhXAP+vL8gsA3hFCOCys57kfLwXw5QC+AsB3hhBuBYAQwqsA/NV+m+sB/DaA/7tUmBjjd6XzBfD1AB4E8Iv9fn65r8yV/v2rfhdfAuCRGOPvZ7v9CIBihY0QsjxYsSG7zp8F8MYY49kY4/0AfgjGf7gxxt+KMb4/xvhkjPG/ovvP+4/VDhJjfATA7wB4MYCvBPB+AL/V//1iAL/Tr4MY48/EGO+LMR7FGP85gE+hq3QBwL/E8YrNt/fzAOAvAXhzjPGDffluB3BJtm2Nn40xfijGeKHf5wuc2yU+HmN8W4zxSXSVkmeFEK7qlZ+vA/DXYoyPxhjvA/DD6CplGr8RY3xXX5Y3AbgWXcXkGM778fdijA/HGP8bOrUpnddfAvBDMcbf7Y/zgwBeHEK4XitUCOFpAH4OwHfHGO/oy/CNMcarlH/JR3MGXcUu53MArjCuASFkCdBjQ3aWEEIAcAOAvBX9+wCK3RH9Ns8D8H+h+0/2UnTvyK9r6wv+PTpF6Fz/d0T3n/Cl/e90jL8I4PsAPLOfdQbdf+xA133xz0IIXwbgMQDPAfCv+2XPAvBnQgh/MzvmKet8BJ/O/n60P24Lcnv0+3gWgNMA/qC75AC6RpNlzL47/RFjvBhCuBfAQgSX835o5/UsAG8JIfxItvwigJsA3F841iUA3gHgn8YYf94oe4lzAK4U864E8EjjfgghI6FiQ3aW2A1d/2l0/8ElnolOIQG6iofknwL4EIAvijFeCeDvAAiF9Uqkis3X9X//e3QVmz/W/40QwnMB/BMAtwG4OsZ4FboKQOjLfAHAv0Kn2vxZAO+MMT7W7/9uAG8QisFlMcZ3OMvnpXRdLO5G9x/7U7JyXRljXFBgMp6R/ui7oG4EUIrgGnM/7gbwGnG9Lo0xflBZ/y19Gf5uPjOE8KvCf5P/e2e/2n8GcGUI4ZnZpl8G4KPOshJCJoIVG7LrvB3AG3tj8FMB/ADmPov7ATw1hJArF1cA+FyM8VxvaP3uhmP9Grr/zJ6Pzs/x2wC+FJ3/4339OmcAHAH4AwAHIYTvAXCz2M+/RNeN82rMu6GALurmfw4h3BI6zoQQ/mQI4bKGMnp4AJ15+JnVNQH0XUC/CeDvhxCuCCEchBCeE+xcNF8dQnh5COEkgO8H8Bl0FRjJmPvxFgB/O4TwxUBn8A0hfGtpxRDC96Hr0vtzfYU4P7+vz/xG8t+39Ot8Fp0n5++EEC7rfT4vQec9IoSsEFZsyK7zBgAfQ9dy/jC6boy/3y/7CDqz7u/3RtCrAfwvAL4rhHAOwI8A+BnvgWKMD/XH+u3eE3IE4IMA7uyXIcb4IXT/4X4AnYH3C/q/c94L4BCduXgWuh5j/HUA3wvgxwE8BOB30XlwWhWW2nl8Ft01+mB/XTxenFcDuAqdcvEguuumelnQ+Vj+Arrw6G8F8K29d0cy5n68HZ3J+R0hhIfR3f9vNMr/xQDuz9SYv+Y9Vs93o+tSPAvgpwD8RYZ6E7J6gmicEELIUgkhvAnAtTFGT7QZIYQ0QcWGEEIIITvD2io2fX6P3wohfCSE8NEQgpZTgxBCCCHExdq6ovpQ3Mt7U+BJdObK74sx/uZaCkQIIYSQrWdteWz6yINz/c+T/T8afgghhBAymLV6bEIIhyGED6MLL/3lGON/XGd5CCGEELLdrDXzcB/e+YIQwlUA3hlCeH5KY54IIdyGLpkZLrnk8j96ww1fgpTcNGRpug4Ojs9Lv7Vp/vfhoVh2dFSe5n8/+aS9bt7Fpy2TU/l36XeJEMq/5TQ/+SEXSrtg/e8hl6m0jbxM2mXT5pV+53gvV/639/mylsnLd2w77wUrzZMXqvQMWs+cF+3C5X9rF6j0DNYuVOGCyVPVLldp3dpzlf895jIlrOdKuywtl0n9buV/D3mupniepvhuAYsXJk3TyVsXTF4g4wVsuUzeT3/pWUzTe+/94NkY43VYES8JIZ6deJ8f7Mawe8nEu52EjRhSIcb4UAjh36FLaHWHWHY7usRkePazb4k/8AMfwKlT3bI0BYDTp8vTyy4rT/O/z5w5/vvwfJ/o9dE+a/y5c/ON5Lz0W04fe2y+zRNPHJ934UI3PX/++BQALl7spvJNkx+K/OVPL+qJE8d/nzzZTdOFShcl/7t2wYDFC5R+i+lj5+cfjHR50vTznz8+TZcp/Qbml+fxx8tT63Jpl62E/N6ly5am+XMlL52cXnppN7388vk28hLKy5eve+by/r5qz5X8nf+tXSh5wfK/tQuWyH/n/1kA8wuU5qfnK/+7dsHyZ1B7ObULCOCJi93Nq72GwPwyyFcyvX6lyyQvT1q39J+VRP4fnC7JkOfKukza6xgez7458sI80o/s8Jj4tpWeq7SO9t0C5hctzUsvXLpw1guY0L5bl1wyX0d+7OV3Sv4G5i/lFf0wXdoFy/4+9/nupqXvkfx+5ZdJe+a01xBY/Ha94Q0hH+Zl6ZzFYrKssYT5MDAbx9oqNiGE6wBc6Cs1lwL4BnQjHW838j8DD6UWivyCytZNsdmvlMUqk9YcLJWpMs2/ZbL4nkZgTdUpVVq0ZVaDUX5ztdaapyyl//Cs1p4sW+xHB5iPsCTuldUi1Sjdd01JkYUu7VvOs9bxls2zTmFd7ZqW0P5vtbZt2f/OY1V6E7UaX+lFlN+ytG36ppVeplrZrHVG3EzPoz7kk78Wpi5o5bqGEJ4B4G3oxus7AnB7jPEfi3UCgH8M4GXoxnl7TZ/EdBTrVGyehm6wv0N0Xp+fjTH+4hrLQwghhJBpuAjgr8cYPxRCuAJdJvNfjjF+LFvnpegG+n0OgBcB+LF+Oop1RkX9DroxdJpoabgulVorNldU5LzUlPQoKUPKpPUtWxfP6ofW1umnSXUoNc60RpPVmPKqJPk8qdR4ugzkPkoNxlp5rT742rbHjq3dj4b7YD4z2rNnKTW145XKVtumtE5Jlapt2zOmUe7pLdGeJ6vHbqOoKR6li+Dxy2j79UimaZl8Fj0ya4u5SpbV8BUd9s+y9sgNERpLn1spVq2UFSs2Mcb70A0bgxjjIyGEOwE8Hd2wM4lXAHhbHyX9myGEq0IIT+u3Hcwmv5KEEEII2XJCCM9GJ2TIyOenA7g7+31PP28UG2EeJoQQQsiSCGEZ0uK1IYTck3x7H+wjDh3OoBv09q/GGB+Wiwv7HZ3PbisrNkuXflt0yBaNUiOX9KX2rZmJLZ2zZibO/67FMZfmOc7ZY/zN53vW9dDSZeBRwqfwLsrjmdt67p13XWub2vPVivcZaQnLLVC7hla4t6dbqYbVnZiweui0bVdmVvY4p4fsr+VF0fqCrRD02m/PjSnN7/8+OCh3RZW6jryf241j+oKejTHeYq3QjyrwcwD+RYzxHYVV7gHwjOz3TQDuHVuwbbklhBBCCNkS+oinnwRwZ4zxHyirvQvAnwsdXwngc2P9NcCWKjYeWiqnk7aWWtQdz7qaudPjovY0LTwt64q505PEKmGZe2vhuVaDTvMLWvPkJbb2P4V3sUTVPFy6P7IZ6VHShrghtf1bIehDylY5biyo1Zq6N/Zd9iiLNaZQhjYSyyGf8Lj3a4b1Uj6HlgSAteiDAYw1Dyda1LxJWU5XVI2vAfAdAP5TP8IAAPwtAM8EgBjjWwC8G12o9yfQhXt/5xQH3tmKDSGEEELWQ4zxfSh7aPJ1IoDXTn3snanYrKwy6vVAeEKrS781dcQTlutJIa4dZyE/e2G/sryGYqMl5LMaUTVPRKlVriknpQalDOde9jPTolY1eVOseYCt0I2JNx2i8sjfLYkfeyw1zMKbNqDFduKxjNT2YR1zo1SdKQrjuWDS7Fb6XtXyLVjZQS15deaxQXFqZVuwPpnaNmtlIwqxGnamYkMIIYQQBVZs9pOZ3yHN8HgWWlrAngik1AxIY65oZbFUHo9fo6LCuKKiDMVG0hIhtKxIkXRMqdzIBqPVxa/9bgnkKC6Tj1pLR70nPGMKj408vhXhJMsmt/Uczyhji/pSa+yX5g1RhhKrVgKXfgDPA9wSwlgLIytl4KztY+g6FazPoUQb0LS0P7JcWLEhhBBCdpn1mIfXxs5WbFYW6aSt66niW9Ra4aXRva0yaGXyRLzU9lvAEz2kbaPNL7XKNR+ONfaeZ0SLKRjU2Gzx2GjI4RNKy6Qi6EGLeCrhMSJ4wkhgq2EW3uvveSat51ezw6XLnwavthSiITmb5D7WFXBzjGXJrV7DVOnYDnn1QAxc73nUWwL+9qhOsRHsbMWGEEIIIT17VLvam4pNqbKutpI8LdEx1nhrMEwtYiC1OqyOW7ltSzZhbZ/Wuslj0zf+Sy1rLadLSX2ptY7HNgLlJdSio4YM5ulpcTcpNUOep5YP1xD1xVrX23xtOa/CujUlsHSNNT+X9U1oSYFS89hYZdL2NVSl2jgs1TDhCU+r4Qk9a4iKsl4p7bMqt8lPc1UqscqedUXtz5kSQgghZOfZG8WGEEII2Vv2SLHZuoqNdm+mVjHNg1nrlsy3lSEJivNmjkDFCGpljmrpotAS9VnlFMfxpLMf062khYrnyzTTcL5t7XZOlUjNKm/LfgDYz0gL8tmQ2vjU/R0eY3BtHdHVWcLqKtJOyboftV6RFvOwp0uqhUHv0LL+M2uJh5fIG2BdKHmDPdEI2roDGGMQ3qN6xMaxdRUbQgghhDSyRzWtna/YDFIGtGRp+d9aVb6keKR5Fy6Ut8l/y+aeVG6GlKkUw1hRYaYK9661MkuNsykGNLQacut+v5ce7t3SzPTsJ6EN2QH4kgTWljvKXUugNza5ona8qc3DWpkSlqrkCnqY4jlaNTI+voWWhH2Fj05AdwNC/xxbl8l6DbRtp3oliY+dr9gQQgghe82eRUVtXcVmJaGPQ6rXQ6rtpVau95jWNlrrOVd9aj4ca7DCntgP3GqFbidkK7PlPnqGXxjCskMwNzpMVybq8/iwWvAogDWVJ3lsjPtu5WnT1rWordvisRniifGcq7Z8sgR9q/4PUFOlgcUPh0fasnJKOLEe39rndmPVmY0r0PLYnzMlhBBCyM6zdYrNUEqtHndfe6m6Xot08qgwqYWSt0Kkt0b7XSunVhbvNo51hrRAPRFDtYZWqVU+RY4vi41UX1rCMLRl6TkeEuUyJjrLWubw2HioKRye5IpjPDZy/7l1pJYA0LodTc+6Fu04tvU+RMXzkp+YHLG2tE7rfo2oq0NlkFbPgJYtwacrZ8+6ovbnTAkhhBCy8+ycYjNFRM0kHhtrfy1JETSlxoqK0mz7loqk7UM7Zkap5dvSAtW2WRebUo6lIJ+rIduWnoNazqRSk7eyTfJwWTaKlrQmWqO/NH5ii6qjHS+dlqWweIQIt+9niFev5LuTvze5pe8x+DU860OCQ+WnsqTOWOPSroxNvo8Ts3MVG0IIIYQIWLHZDjwVcSsaxx0tUWoJaa3NligQ6dfRtvOWqbYPS+XxHiv7PTTKw7vtFOTHGXJpW7fdKdJzqfgOXNtO4MMpKSoehjx7tW9CSXmUYtQQX44nv1PL/hfYBvVlajwXbmFe+Vn3vAK1/xrI6tjqig0hhBBCKtA8TAghhBCynWylYjMkUVtLt1UyLYaWcGkpvXu6l0ruM81l6wn3lpp4qSy1bT2uOaVoHoaYiseEcq+zkeLpEVxgir6WTcbqDhUXxtOVY6VskMZizWhsHUcbXDUnLdMSPlo9IB7zsyyLZLbNiYau5pbvyCYxpB9uyebhlt78tYome6TYbGXFhhBCCCFO9qwramcrNkPSp7tylGlVb8tdJpUZK+ZPay15wr1rrrV8fou502lOtq7fFOLC1O/lGJNfSy48z6MxCZug4Ggn1GJAHpCgz3PqLdn4h0QMy0j6ISHjpXQVSzUPt4SIr8NV700caUl11rryYorT8Yx6oyk4pW02WQzbJXa2YkMIIYSQHio224us4HtaaZJUeT+0qt61FlDeQr1wwd4234d3ZMaS+qKt42kta0n9DFrCWccwxCI01fG9t6ElHLRIrVluNeE9DGrmV5j6Q9ng4fIoK5o3yxNarXlrSmWRxfa8wppHyDqOVpbZN60lnYPnpbJkBk19nvoFrNFyHMfH3yNsJTRLo7UNWQ07V7EhhBBCiGCPalpbV7EpdaXK+bXt5LxRjdmWDtbatsBiS6ilBdQSCTHEYLLlL0atn9taPiSqYVBU1LbiPdmR0XZDUtJr77cn2m5Isj2p3rUM3VCylIz6TnkV5vzvIQkZ0zYXL7Zvu2xG3GBrZJlajlZrm5WzZ+bh/TlTQgghhOw8W6fYAJkHplD7rbVuWgZjdHls5PwScpsh1ni5/3wfWjNgI5InTE8pRZDW1W/huR01i4JHoGuJoGoKgVl2gpYxjDYd+dAugfy7tI21rJZOyqIWHWXtz8rHo62TyprybwFZDi6tcGPHCGhRhxND5LZlhBE5buKQR9NKZbYR7Nj/Axb7c6aEEEII2Xm2UrEZwpAG77FMnhKtWm61ZGQLqBQ+oYVYpGaZR14YYu4YugxtrVjJELWkhLyk1mCFNYGrtHzImI6ayFbcxxC1ZdnUUurmLDHayqN0eBiSD8aTi0bz1KR1TxS+si1BcNo2VtbiBbW5lnzFWteDR6ZYm8Gkx6FOej6dtUtrbb820YQeG0IIIYSQ7WRvFBtCCCFkb9kjxWarKjYxlhNwAX6ZeVRIt/y7tE6pv8E7TEI+TxZU7iOXdb0OV08ZBjz8Q7pnrFPXTrXUDeC9XFZZWkbD0H6XbrdcZoWO4uKIvgktC6WVOl7bV7FwCqtKwpah9SZYodst3wCte8fq/dOWaSbi0nGs/VsDfWrz1a4oK6uc9uBa3VYt3zZZ8LX3z2QoD8WQ7vDS++7Nu7pUNuE6r4j9OVNCCCGE7DxbpdgkWpJlWfO1tOazdT0Gu5ZEd55Q8XRwzUQ8ZP/WukOXj9z9GPKWUMoH5rlc8lbVGqgW1jZSmbFudzXcu8SYWGS5j2WpLkMMzg7TsKTlWzAkkl5i+U+1VnmL6FZatyWZX5Uh3yvtd2nZmER9Y0alHYpyo2vDJXiXbQw0DxNCCCGEbCdbqdgkrGRc2sj1kzdQZS241OpI89JgmAkrd38tT/uWUGsklKxCmn/F8rNY+5X7b7EkecpQO45lcwhQmuqezHO1pnxpmxYlRWbCbImP1n5PRMtu5Sl7/CxyXWswTK9Co/kDS+tYw8W48jEOCfceM0aAvKjWt20MG6A6aNZFjyWJHpvVsNUVG0IIIYRU2LOuqK2r2BwdDfPYtDR8Z620PkV58FS9l+2x0SIKvMduxZLDlMNakUjWEAQaLftv2Z9M1GfZpbzqixUw4omOUpv3nnAfj0oyJDtdKvCyhl1YAi1KjucSePavvbKJUqK+WtCbNRimax8nKg9uwnrYW8YIkNuWDEfaN8yT1M+jOE3ATEHtv/1T2BXJ6tm6ig0hhBBCGtmjGthWV2xKeSESsqFrNTprFoVDTxNeW26tayVokWjem6kZIIdZXfA1hohhJVtAS56IIR4erRFrqTxadNSxMtb8MlaYjJVTv/QbqJs7LPVQSg/WCzghQ6Ojhth9PHYl+Vs+c/JZLL26nttcK5Pk2C2oPdxD1JKSIU574UrfttrLuRFGlDKeb1zLd4QsF15mQgghZNc5OJj2n4MQwltDCA+EEO5Qlt8aQvhcCOHD/b83THGqW63YEEIIIaTC+szDPwXgzQDeZqzzazHGl0950K2q2MTYKZul0XITNbl5aeq59dC09K3UupxaHs5lhYjPuijqh9O6bCzVWVsn7SOX4rVuAA81k2+p3J7z8PYCAPB3QVnxvwmra8oVI4xyH4sW9m3hGeZBsuxu1p4xXmhPl9eYXmOr51Fbp+k4Q148z37kS1QqVEtuBnmc2j5WyJZm3lgLMcb3hhCeverjblXFhhBCCCED2FyDz1eFED4C4F4AfyPG+NGxO9zKio3HE+lp1WhpzOV0tHl4CHL7Ic2/UpypZEQruaVBJ+d78ndpSkq+7ZDLnC5HTYUpzfN4MGuN2WNlvlh5+CyJQHPIl7atPewWqcBWljpt3gQKzdhXqSUZXkLOk0XNl9da8C3+79Lxal7xkql4IVWF9iJ6DMGeF0PKqaWbJi+q5q73fBy0l2uD8BiOt5xrQwgfyH7fHmO8vXEfHwLwrBjjuRDCywD8PIDnjC3YVlZsCCGEEOJkOR6bszHGW8bsIMb4cPb3u0MIPxpCuDbGeHbMfre6YuNpMLYuH4yn73dIlrplUWtZe2JTe1q65BOlMGktlFqGz+YNvSF+iZqSYqkvcp2k/pw8qW+TljV5bGTzPD/RltBwuU3td+kmymUtYeQeVuSt8VC7TB51R85Pz0hpPSmktQz3YKk+C2pzTY3J/x4yOuwYBaXlw+FBk9CW5M+Ruy0dpiWzxz4RQrgBwP0xxhhCeCE65+Znxu53qys2hBBCCHGwhoZ0COHtAG5F1211D4A3AjgJADHGtwD40wD+cgjhIoDHALwqxvGpzreuYhNjW4O0JfGVGlRyotCqaTGR1B6o0rrrasU2HNfTJV4TqawIJE9UUUsDUZZJO17JdlBLzGeVSe5jnrYddcOEZdBoiaSS23p/l5ZZEU7aC2cNDaGwK34EGVSWz5NY45fWbFFNio01Foh8uEvr1l680nestm7pBdTKpi239rtkNv55XUMBY4yvrix/M7pw8EnZ9FtBCCGEEOJm6xSbnI3y2HiSG6R1LlxYUiEmxNFktAaMqykbLRFInjwzniEVpFKTvA8ySirfR1omt7HOo6oE5df24sXj87RpWg9YVEw0qbHFl2MhW9+lZ90jI9SWyTJtULNLG7IFmD9X8hmUgro1sOUQ4ax0m6toL6a1bul3zQCXnlePZythqUjyd03B8bJEFcP6tq2cPRvde3/OlBBCCCE7z9YpNkdH5fQsWounxQ6g2ho8soJcbrHq6rvHaOTJDKtgDQTZor5o21iqjLyEVtRBTT2yGoyaUtNiVZg9r5YaJhWc0n2p5a+xwmRqPpySEUS2sIcYQVqUGwfLfnU0FaQ0P81Lz4IniKxF4NJeTSvR9MK3y/Odqr2AudJck2BLL4am1LR8Mz2/tXN2qFQp/4+HrRNAtq7Aw9m6ig0hhBBCGmBX1GoIITwjhPD/hRDuDCF8NITwfesqCyGEEEJ2g3UqNhcB/PUY44dCCFcA+GAI4ZdjjB/z7sCSbTVWliCpFDO8SbRcqIrx9ODgsJ/ON9G6nLRI0tK8FsOxh9r+UiK9vItTlkUm5rO64dSQ8IuOLhytm2mqbdL91QarzE+gxXhcG4sg/+18GUvPlbXOlCw764Kn584zkoWcP/v7ROWh94RwWy9e7UXP77F2k7T8C9Zxavtcxjrbzj6cY8/azjTGeF+M8UP9348AuBPA09dVHkIIIYRsPxvhsemHNf9yAP+xsOw2ALcBwFOe8kwcHdkRhPK3lRfMm98sN5QtNBjXHsdXQDPpWRdhxIUqKTbeMO+SsKWZhUuNtqScyOhS7RLk23gMzfI8NENwqUxq4zgP3Zam4dpIh9YyzzZSqZH3Oz+RmvpiLWsxD8t10nO1EV+n4ViXYIro+Kachx4T8ZDcDNr+Sy+eptB5RqEdUjbPEApL+G5bWT2WNKpDnT3z2Kz90xFCOAPg5wD81XxArEQ/WujtAPCMZ9wyOtUyIYQQsnewYrMaQggn0VVq/kWM8R2t25daKrVkVaNHodhkhSbhUWq0jvtS1q+aYlNQLWrRn1ZDq+aBKZFUEkupqe23VCapvmhqTO7L0ZbNEhpa6ouWc6DkefIm6rP2L8nnpxOQ0plVptqojlbY+gA8lotVtZK1b0spPUULXjGsaF9qUVZqL601uqM8+SnMb9Z+rJs6RE0awD/aUQUAACAASURBVLL9V2Q4a6vYhBACgJ8EcGeM8R+sqxyEEELIzrNJjfEls07F5msAfAeA/xRC+HA/72/FGN9tbXR05Et45QngqDV0B6k7y354SieiKTSeTv6Wzn7VY9P9LEU4aV3kVgBEatmmPmrLl9NyjzQlRWuoWuUdkmhwlF+mpKDJB9caONNSc2rItPXWC1IzgnhMJOK5smwgFuv6jnuG90h4hC1tnnUp02VP/sAwJJpIe2nzebWopZICqMmqltxW8whZapJ1rmK/nldpCHtUp9gI1laxiTG+DwUvLiGEEEImhObh7SFv5UiPhcaQmne+jTrk2jY0IUud/C0RNco2yTtyeDivp9b8MqXGlGycJe+LbD1Z21jU1JZ0eXIvj1dxKpVpIWJKDpdQmufxy3gHzrT8LC0vgmxptzwjFlM1hzG9j0Y+I6vKf+W5FJq6Y26rKR6WzGopKp51Sr8963oMci1qkmf/ZKfY6ooNIYQQQhzsUYVuqyo2MR732JTG65O/PR6bUfkgtIdlqoeoxZugzZct7vzvWmRNqQyq12Z+Q2rqSCmaKCkltfw1Jc9F2kYGaXi63i2bQK3cHh/IXLFp8NhYik0t140VtaSZOizkNkOUP+u4yrpJCQyhrgTK5SVqdpBNpcUGJ5fNvpVpgWVYa7mo2n4seVXbhzXfm7+m9NJqZTVecM/ndmXZ66dkz7qi9udMCSGEELLzbJViQwghhJAB7JFis9UVG0+0aety7zqDmGLHJR1UDlaoPcAes2fpojq7PkpdUUMyoWsqdymMVnYJeS6x1o1khaBrU2uYh4VllunW270E6EnwWkLEW+LktS5Nq3tMbmuVaQJkRDpQ7m2pbbsp3QxWj13LNjNa+uw82TS1dax+vlqYt9xH6/49+2tkQn87WSFbXbEhhBBCiAMqNptLbh620DyxrcdaOUOaCNI5O8Y53RIqLBUbY1gBzXRbapxpjcFSyzsfT7JGTREqNUxrA2bK4RhK+10YSqGk2NSGUrAMxzUFx1pXkp9ILa7Yeq48L5zTcHx4uKgEDoncHfNdl6phSfBKz0ZNIRqLNhpKCdU8XPqtLWtRSbQhFizGGJqtB8ET7i3Nw4WMDBqayL0pqt8xaB4mhBBCCNlOtk6xybEasS3b17bN50/aGLN8LXIdT1NRG6yw1IrWMhpaqkLVY7NYpJpHpdQ408K+S43BIS12qR5Zio1X3bG2MR+02rUtKSxDttGakdYLk05AhnmXFCH5XHqUQFmGEYYG133o8WTa16KXt4WF8raEcGuSUykMW76kLSqP3K+Vb0Ee19pGzrNe8Aql10b7XG/8M0LFhhBCCCFk+9hKxcZTM9ZaWi3+nKUxdWiWzEbniXiqtZIHKDYzLwnmwyvUFJuWCCTLl+Oh1qgcUqaEFcjRpNgMyZfvkRprHh4LT9Rd7RmxfD8aMyVQ99jIgBePDcSa732ePIE2NQFkLPLyTT5orzyBNCptvkzbXio61rrSnOR5AeVxSmqSp6zSY+P4HA5h7WrOnnlstrJiQwghhJAGWLHZXI6OFi0kgJ5GIzGoNVPC+3AMraLL7aS3psUTIcvqUQpKLXmnYpPvP7WypZ+lRR2Rg0nOIjyyBtkUwRc1FcBax2qVz85NDqVg+Uw8yop2zzzb1JSa0vzSCyd/15q6VtKpynm02DRyNJuHx2pRs6RY6aQ2idntSOVvuZiel9Ua70SiPXMek5vXtFdbpqwbEYpFtB7bEZYwsmS2rmJDCCGEkEb2SLHZnzMlhBBCyM6z14rN2hMpjdEy821kV4H8XTqOlP9lF4iVjn9AV5Scyq4pQE9+l4pSUpbHRJXK+aVupVr0qqWEz5jAmF28H57EfJIxQypYI45r23jnO7HuVT7f2taaL58NzSRsdTt5jlNbd+ndWq4Hd8T+Sve5FkZulWmKsVkaznlMd9Nktocp2TPz8P6cKSGEEEJ2nq1SbGLsGo2liqeWZtwTtrekxuUiQ9LMa2bS0jYJzUTsCf+1QpIbtjmomIVLjSct55flT7QGyJRYUaXa/rUGo+V3VFWk0r3Twrst9UV7FjyqWy0ff15oWTYZ9j3E0Fw6duW5Kg3V4fGsttxnibzPLcMYJCw1Uc7zJA2s7WMyhrwYlqwqkcqy9VHwDrZpjc3iOA/tUZRFzv+u/d+ycabiPVJstqpiQwghhJBG9qwraqsrNlOrL0urYU9Rlfd09nq9NqV5HsWmKdy7m46J0tROz+oib/FNtbTcvfNzFga/TLQ8uJ7lnod+SOI/eZKeoT9qIwJaPq8J4mit50lr5FsKYGKKRIBWGgEL7TyacD2wDSdZW6c036vUlG6IZq6yjHGe8U9E0YZ8qjdOmSHbXbEhhBBCiAMqNptLjBtaQx7SyrTSzNciXIYc50R2u7XjlI6b5l28eHyZ0cJOaoUcWsHyRmjraN34+d/p0FqLO1+mNeA8ngiP/WAh+MN6NjTzhpWQsXb9Ldmq9pzm82teLY+C4z2uQUuQjKVmeMZ21J4Ny64mt615uUrzar89DFJyhh5AFrAlPFHbh+fFGyL9OtY9qnzaLDuZnO+xk5HlsnUVG0IIIYQ0QsVm+9g4Fcca795T2FoEUgltWamZqSkDpTJqvgmHL0fms5Gt5lIj0Io40tCUm9I62v7MYRGc811M7bHxbKM1FbXIp3x7eWFKz7Y3UstapkzzwVWB47KEdU9rKltpbMfZURTV0BN9px2/ZANpES1qxzFpeX5a5CNtXUvS0uZbspvXtGctMy6uLK41tILXh7NR/yftmXl4f86UEEIIITvPVis2HvFi8lrzlLb5Umu2Vv0f0vKaOirKoyYJxUYWyWrFauuWuvG1RsiY+97S+jd9DUOelRa0/Xu8L55BMNPJyvst89l4ylab76Rm5bBUkRZ7Rk2hsWxMHi+PtsyjQLXkvBnVSJcPd35imtKXfHhjvDbWTZQnX3ogGkIx5eCXHovbbNuGxN9rz3IPULEhhBBCCNlGtlqxIYQQQkiFPfPYbFXFJoV6l6KXEy0ezJWH3rVkgWoJy60dT8as5vO04w3ptip1RTmHVqgtA8qh3JpXcYi5cyqa9jeka0bq2mN07vQSlPrU5PNjUXu2h4StF7Y9PLQHV23xuZbm13yzljldbjPEC1s6rjeS2urWXRlWkj/tJfX0NWvT0o1oGDDT27tufQ41hsQILJU9qtjsz5kSQgghZGWEEN4aQngghHCHsjyEEH44hPCJEMLvhBC+YorjbpVis1HUWqZWfKA135tm3spAp21TqrFrzrdSRqohio3Sah0SQVr6rZ3aOltGkxx7Cvd76R5qeJQbue4Q8/CSXJQeVUTOLy2viQhjwrwtk3KidB5ynRaVaunIQliZDLUX2cpo6A3zblF5CopN7bdHwG7J5rEW1vOg/BSANwN4m7L8pQCe0/97EYAf66ejoGJDCCGEkMmJMb4XwIPGKq8A8LbY8ZsArgohPG3scbdasfFESXuXl1haBddqJmjGH6mklNZLy2rKTWmelqivdGxtH4X9a0MrJCx/w5Ax+SyWcT+bfFpWDK82f5P6xS1zW2LIi9agBMqEj57nqRYhXIpilstaopg9ak/NMmI9KkP2j4sTKIBD8Mircl3PTfTIVh7FxjmUQn7ZvEMprNXPKdlc8/DTAdyd/b6nn3ffmJ1udcWGEEIIIWvh2hDCB7Lft8cYb2/cRykT2OhqICs2Do5VdGWDZ0gUiPbb2o+l1CTSMk+Twuvlyf+eoBPZihSpbbNOrMujsgkFl9Seo3x+y/gBmnem5X0YgRWllPCoInLdWvRSCc0yYkU6DbGMjKJ07acwicjQxfx50iTZIaOFWnLuAI+NpdCUfmvzrPm1ZUtn+u/R2RjjLSP3cQ+AZ2S/bwJw78h90mNDCCGE7DSpK2rKf9PwLgB/ro+O+koAn4sxjuqGArZQsTk68l3TpdWMh0R71JoFHst97fg5ngiFlsQNE5yHHFpBFhVoaxUPQXbxe1pr3ufI6oM3jSDWPA0pT1ijOS4Dz0XR3geP2cD02JSncnnLui3eF09QV80O4lmnVCZthAOPIjTqe1JizEtaM895PgqacSqfp0xj1vtRG0qh9LjWhlKwAkv3jRDC2wHciq7b6h4AbwRwEgBijG8B8G4ALwPwCQCPAvjOKY67dRUbQgghhDSyhq7xGOOrK8sjgNdOfdydq9isbLCxKcOwPFV72XRoUW48ylCpLNqxtOZNwzl7PBGbRIuwlUgtw2A1x7XWqyeBypAPlQzzGcNUzVAt6q4Bj1/G87vmfWkpywmRddsa21GuWxIgWlSetG6KShyk1LSE8pRUYcA209US9FjLNEnNse4Q4boloFTuW9sPWR47V7EhhBBCiGATgxmWBCs2hBBCyC6zuXlslgIrNhmD7rtH4q051Kz91sKyLaSJOE+s1hLfqK1joZiHl20MHiL5aqbi/O/aZbLk7ZlIXjp5ywGqoa2jGSyBxW6e1EWw9sxh47B8qLWhFKwunHRZxnRFyW6mht6SY+vKbiqt26pkOK5+nzzd0yVq61hua60LqvS79n60OMB7PLEanq6o2j5KsCtqNbBiQwghhOw6VGy2g7G1X83XNsjYapl7a8Zfy2WmudtaBsGU80t41KQhDtoNpBb27WmdtUTsH3piebUme4uhMmGNCLll90rDE0Evl0llw1I6aibi0mXU1i0dR6o4Uo3JxdXa/qyxI2dMIWlOjecm1k7eeqeUfRxlvnntsrSMCSwZkzaCTMNWV2wIIYQQUoEem81Gq/nWas9DKDaIay2flg5Wq7mvmTqsE9IGwbTMI56Qba8PY0lNFKtotRZXCU3AslpnWgvO1aKTpghLsdFapKX4X8mqPlxDjCcTH9qzvKa8lmwgUgVJj366hdZj7RXfrGUldUdTaDyP1ag45pZvj0bLWBelZ13bn0d2c4R7a9+LljGBayJ7zlqVmz2q2OzPmRJCCCFk59k6xaaVlSWC80QbyN8tVXrLlyPRBsFs6Si2WnJDoqQmxDqsZ6zQ2niNluXJoxCl3HezcmphMqVlmrfGE1KTKEkRsmWblsnoqLy1POblWdGLt67BLz37kH6Z/LZL1SVNT55c3H+apyk0aX5R1DsvH0Yc/z1V3v+WbVsUv1o0lCWHiWlKmFl6/4eMMNOSn3QjPDZUbAghhBBCto+tVGw2ovar4RkM0zNfG2nNkwNH1sxldJRVhqmjJmY+n/bdSdXFijZoubRStJDr5qMNpFZwahVLNSaNO3nqlL6f2dAKssmd/z1mVMSaP6eEVHUsX4MHr0JjRWxNTM2K5LE8yctkFTWd2hSDYFqinlRoSqLewlAK8qHUXjLADvlLTPkRHqJOWhJaxVszRJH1WJFaIiVXzp6Zh/fnTAkhhBCy82ylYjMGT6oE2QCetX4AfyesJ9LJI0G0WO3lyclm5thR2Zbc3NB8MR5BKJ2qx/6jXa5S8FgqU2roSuVG/gbmKo5sJJ8sNcdrCU2sqA+5TfqdCpAjw3yGeCJaDC3yt/XiyTIaaMX2HEYrWsmbom3jibbTPDalbMLaNC+TXJb2L/04+WM1SorQzCOlF3TVEoTH7FRRbDyf2yGfaG15zsoGaS6xR4rN3lVsCCGEkL2CXVGEEEIIIdvJ3ig2K6usepLZeYx72m+5D2teqW9FW9fqw2kJNZcoF95zOE9vXE09L+235VmQ3RVJ7pem4bwrKp1bmjfrijpVMA9rfRGyv+GJJxYLVZtahmCtj8XT3aTFQFv7T5T239J9NQGWV1X23rZEJCe0sGzLE6vdds+6xTB2aRKWv1v6ZeTy0rIpsYzycn5DV6fnlLVEfTla0IH1DdqIgBcqNoQQQggh28feKDaTMcY83OKOrTnSPC40Ley7Ns+7/wF4TL0tfkctxLJFufEID2mqmYZzxUbOm4V9n+6a1qEUy1tTbvJtagMBWlJEy8ivQ6StMebhWnx2gZZckS3qS81gbm1bu6WlefKWWcn8pHqYlh8eFIIDNKWmZaRGTwoLTQG2FEB5MT0pB6SC44nZN8zD8jLUEvWV1klYz6AncejS2SPFhhUbQgghZJfZM/PwXldsvIPkjQ6JrqkvVnNAS9SXU0qHX9ombwZ6OpMnRLPpWA1GrXGZqyMyHLu0TkKesifFftpPWid5a9KlTNaX/NKmbc6f76bSh3MyXzktTCvXlBtgbsBIhdES9U3FmHEFPEkDvUnYMCxLgSYMlA7jsQ/VjqNF6k81pIIa5p1flJqnpvQCtozq6FV38t/aB9dzsVuG6BD7GzOUQsLjy2lRi8ly2euKDSGEELIXULHZLzwG+wVq1fNS86ClU7ZGSWFJ87T87/nxUgtoyjLlKBfR47Gp9XtbA06m36U+bembSOt6godk69iToE96bJK6c/LybPyFmlJTMmhoXhorTEYuk89PS8iY5cdp8dp4vUL5Nkto8ea71y6PZ3spUmmDV+bH0W6ZZ0iFBT/O+ewhlC+ENi2pxZ4IKnmcKT15pedW/tbGpCjMaxGePF69mqBVuiRr9dbsIazYEEIIIbsMPTarI4TwVgAvB/BAjPH56yyLSUtCAqtV09KBqw2/YKXE1xSadT7QqdUkPC9WI1BrAc2ii4zgD62xWTq2KKJVfFVckF6bfJ5UapKd5smjuUR06FVqSiE1QyKREi3KnNbcLI1FIOUvq2UtzSIOxUZrLcvlHiy7j3a/rVdKU2o8Qypo25TmyW1mQ76UQvO8hpDSNp5RHSXyhkzl96p9w4yb2eKXqX1+rW3k8rEj2EzOHlVs1n2mPwXgJWsuAyGEEEJ2hLUqNjHG94YQnr3OMhBCCCE7DbuiNpcYO0nwxIhSD8k5VsRKdAXYbrGW4QtazL21rqeSUXAKLAOfUoQWiVfrksr/1kLCrUvcZBbvSd0JMqS7FO4tQ8HTumkKAJem/qkh5mHZX5JOLO3z0Ufn28isbp//fDe95ppu+sADWOCGG7rpgw/2hb30+PL8GZr1s1XGJCg9I1rXWj9N4br57sd43tM9LHns037kOnLb0rxaT1o+TEJtUHdPMr/ZOqWHXTMLW7kTvK7Y/O9a5rl8vhbU0IIn0mOh/7CfOE5DO3VraAVPN9NGdEXtERtfsQkh3AbgNgC48spnrrk0hBBCyBZCxWZziDHeDuB2AHja025xywxTNAqKTGEelliZo7TfHvPwFC2ksfTHruX+AvTGpRWhKltc1rq1BH1p3ZIimNZJaotmIs73q5mIT5+erzsbZiEpHtq0NCpiOqGrruqm99zTTV/84m56883zba688ti6T6I77uFRF5v+xFG2/55LDrplTx50yw7RH++hh7rpww/PV77rrm76G7/RTW+6qZs+9lg3lYpRfm7axSyYP2sJ1azX0+OxHpOYT1NdLBVGMw17zMOzIRTOi4ce8MuenpepdHE9wQzAeAO794YYN9FSiRPa59VSeWph3lRn1sfGV2wIIYQQMhIqNqshhPB2ALcCuDaEcA+AN8YYf3IZx2oJ7VzAqnoPGRHN6sCtNS88HbkJTxryIZ2/tcRq2d/JH+FJXa6FWHo8NnKakuSVLq08nmxptzQkZb6wfD9J3Xn88fL8vLwnNQPFZZcdXxGY+WKe+MIvBQBc8vAfAADu/7pvAwBc/9TF8N94opMGUmhwulOxV2MuwaIYGnHy+Lq9ynN0VXf8w6QUAcCznw0A+IM/8icAANehK9MTV17X7f+uO7v1zp2bb5Okq4rHxhOeW6Jl3cQYgWCIYjNmXTXXQWleTQ7N/27JODdElqh5AT14DHLi+fGI6R4lUDtlTcmplWEl0Dy8OmKMr17n8QkhhBCyW+x8V9QkXpu8ll2renukiCFockNOy0BxK6LW/+zxT8ghCtLvfJ5cVrrkNRuAjIixsNQ+LdKlFBU1U2xOKaEvV1/dTW+8cb5RUnHSfq69tpukFmOvkoWs2R+EIiN/l9DWOZyd6/yk0zFTcePBdcfKGL+kU5fC44/Nd5TUm+TVURSbi4XrVQv2aaEUFdWSD04WWxv80krQN0ixkd4aj2JjGUFq3pqxGeeWqdA0GKWs0/B8ZrVtPPkM1+q72SPFZn/OlBBCCCE7z84oNiWvw9h9zfBUs1sGgWsxmki7fstxxmxbotZackQmSFWmFMghFRo5Nl8pkkoqNaVGbE2xkSnx5bFytEEMgcWWdQoMSuueysbAnB3noI+O6n0r504+BQBw5vJCunyU9zM7D4caMzUz787B8RT6qYyzMuWF7pWmR0536s4Vp/ubKN6LixfnFzepXZ5IlCGPe23ohNIrkO5zy5AK8rlpUWzCkfKQt5jPSttokoMlfw6JhpK0SKUNx0kqYsvIEJ48SbX8NRsZDUWPDSGEEEJ2ClZsdheHid5HS3V9TASBtV+JbKIOMY94kHJFCSUiwWoEakEZssvfyjxsNWJrLfh0WlYkVcIRELYwTXlskkUmL1+ym5w+3bUyL7mkL0vBL5MY4ptZNtUyZBcqnVu6Hk/2UViP9dcpKR+5JylPpgzMBQhP5MtCWWvRkAVK97um1JRS+IzJY6N6ayzFRnvxrNTfMrRwrGlE5l8ak0Je23f2tyZ6l9AU5dKpa15AK9XYRnhs9oi9q9gQQgghe8ceKTb7c6aEEEII2Xn2WrGpGY6LFVyvljg0c1jSM2vxqyV9NYTFed7yDemj0xL1ZX9rp2qp2zWfo6W4p24LK0GfdknTJS+p8xpWFGra9syZbnpdHwH9xVfeN1/pTDfQ5IWLobq/XUX2mKZuuNRTcc2pR2brXjxzBYB5150WoQyMG+e15pPPX7XUnaQlbZTdTfmylnDvqmm49GJojnzZh1taJl37lutWuq2tPsF08mP6ZQa8IC2ZOGph4K3blHr+VgrNw4QQQgjZKVix2T1axAxXuLdHgtDWHYIcdG7ZlGQsTU6wFBsl15cV7l0LDbcUGxn2nSfz026RbIVbRkFJqYWdeOCBbvqCF3TTNC4kLnuqup9dJxmGS6hqVZK8AFwlxs1M437mA4tKpsghJ78fpTFJ5bOgGYOtZZZiUzUNWy+GlCY8g2B6JIiWcO+pgxjy4zSoxZaq54kFqYWIl7ZZm1Kzp+zJ55QQQgjZU9gVtRsMuYeyVbZ0j01L3n8PLWqOJzV5DdnyKmxbG8gyP13NBmDlFJPLpFKTt5Ry9QaYX65039Pp5NvUrE5p2zyE+4bONoM/+Se76fW4v/vjTKfUpMEkiY9c5Un5/a473fluvuIrOs9NUm7S6AyAP3GaR82VGQ7yR11TWzQlJ/877TedV2ld1VuTDGVWNkrt5SlJplqWS49pxOObSRekJpNYtEjvDdQ8NZZYpX12rcwf+0II4SUA/jGAQwA/EWN8k1h+K4BfAPDf+lnviDH+nbHH3dmKDSGEEEJ6VqzYhBAOAfwIgG8EcA+A94cQ3hVj/JhY9ddijC+f8thbVbEJob171pMSfdT91jLPldYZMhyC1popzR9yIp5trGx02vz+bxnkZTUCNU+NNUyCXJYanbKhqh0zJ+2r1BqXZZPr5ts897nddKbinL66fEDSzCwBYO+7OdWLFv3oDEXFRrvfJYVOLkvI+5wrKtKjpQ2lMDRB34K3RlNqhmSuzLfRMs+VIqjkNtq3rTTCqPZ7SXKGRyCv+WZaBsMsbTNFjtbRrL4r6oUAPhFjvAsAQgg/DeAVAGTFZnL2p9ONEEIIIVNxbQjhA9m/28TypwO4O/t9Tz9P8lUhhI+EEP5tCOEPT1GwrVJspsZtM/FY41uW1zpnSwyp8i8rT4Sm4BTS5ctT1YIz8r9r09wrk/5OwxVYfpySr6d0Wnk3vmx9peMkXvSibvrKVy5uk3jyoGuWsxUxHhlVle7Z1VcfnwLAnXd207vvPr6uZ6QRTZSU/pl8e5nPRlNl8m00r83MVwP4lRrLfKaFFlp5bOQ61lgj8oKV5E8rx82UiIjMFjRRyTPMSkvOm5WzHPPw2RjjLdZRC/Pkf3ofAvCsGOO5EMLLAPw8gOeMLRi/tYQQQgiZmnsAPCP7fROAe/MVYowPxxjP9X+/G8DJEMK1Yw+81YrNsrsMiwP6tXhevLSkum3ZX8sFskYE1JYZaXK1fufU8pGDFwJ6Q9FqZMpswVrgSH4sLSWQjI7K100t6SylCgDgf/ymz/Z/zRfEUycXzo0sB/kI5u/s87+0uwH3339ckrlotORrYmQpW7k3J02u2Mjsymkfhwd9+R/PHlyvUmOZzzRzmzUIpnxRWqI3S+qMtv8hL0rN91fAOg1ZFEtM9xa/tI8989i8H8BzQghfAOBTAF4F4NvzFUIINwC4P8YYQwgvRCe2fGbsgbe6YkMIIYSQzSPGeDGE8DoA70EX7v3WGONHQwjf0y9/C4A/DeAvhxAuAngMwKtiHO8iZ8WGEEII2WXWlKCv7156t5j3luzvNwN489TH3bmKjXbvPGHigzJ9e0IWh5iFh4SIT4kVF18L+0Y9FFJ2SZXW0ZL7lczD0ispB8PMj2l1ReTHzY/56KPd9Iu/uJt+8zf3K6S+qcKYCHuU6LOKNZTCGOQ1Pnacg+6Ffv7zu5+f/GQ3vb/PlyiT7uV/17qg8m+FNAvXBri01nUNj9AypEJtJFlrG+uF8XZF5S+T94WYuL9myO5qCfsAPTp+6ryrk7FHH6T9OVNCCCGE7Dw7p9gMQVZkF7J2Lyvce0w1Pt922TXx2tAJBQVHiy6VHsZcxPI2Lj3mYWtIBdnSsu5/2v5jfUqpf/SPuun11/f7ROf6XPaLlCsRRVM7WSDd36uu6qZpoMyPf7ybXnNNN/UoNlKpKakvaV1pJvYk6JuZhmUyvvzvIYqN9mJYY5pIpable+WJpdcU7InjpD1xHrIonsPV9usxHq8FKjaEEEIIIdvHVio2LcMiWEMqDKrALiM7U17F1zw1U1T1rUxkJcOBRBsltKTYKI09TckB6hndrQR9tcEw83mafUkmWgOAZz6zHekXNgAAIABJREFUm77mNd30i859pPvjhpu7Mp+4HGQzSY/jJQfdQ3DJxW68ha/6qk6que++brkcHDXfVlNqSuqLNi0l6FvYj0d9qU1L6osmd5YyZNaUmjFh2db2U8sYyjdsSG7T0ql7L8tGeWw4ujchhBBCdgpWbLafScaDHNPhao2aNoaWfbSoL9a2NU9NQbGR3pdalBRQH0JhWR4bGW2VD6T41V/dTV/wgn6dm78MAHDyoNtx8rssK+qHTEAvi8SrO6Xm2v6r99hj3fRTn5qvmnw4mlJTCgj0JuYrRkVp3pqxHhtNmbEkU22dqbyA2sdgybTYHafIs2rtY11BrvvGzlZsCCGEEAJ2RZECpSr4mKq3p3kwJPfNFFgGptR8lUMslAbBjMenNeWmNE9mdLesBC0eG5mWQx73ec+bL/vb390lPblw9fXHtokHQ5IeTUNShxgd5UOqaVde2U2f/6xHAABPPHHFbFlScTSlRkY6Wcs8ik1ThFPNW5M/7No6UtIsvYBymym+dYCeKGxK2cSB53SsnDRe4am0zao/5/sKKzaEEELIrkPFhhBCCCE7Ays2u4d1T+WySYZW2CQ8J2/Fxdf208/PJX9tCAXNEJyvq+UWK3UltZiHa1GsDz7YTf/CX8jOrc/udvJEdwJPHtEkvO3MuvD64TCe/ez5sk98optqZuFSSoBaeHcpQV84UkbfHmMeLo0Foq1bGi5B64KysteVhk7I5+dorv0pKBxPFsnT8z+GliR/ZLnsTcWGEEII2UtoHt4fJknQV/udz/O4zDbRXaZJWmK+J4mVVHA8ioqmypTWtdQdbf+XXdZNv/7ru+nXvuCR+UYHXfxvUqM26dtAE/EwpJn46qsX/3788ePbJL+8J9lebYgFAMDFiqJiKTbSvVqSPzXTsHwxSi+tfOHk8hJSuZFDK4xlScEUsnjy0pZUGK/wlF++Tfys7zJ7XbEhhBBC9oJNapUtmb2r2Fh56ZoS9NWq4OusontMQto6MpQ7p5KYr0WxKTXAaj6ZUvRpS8Z42Uh99NFumhLy3XZbN/3M+Xn471W9mrM/n4T94/z5+d9pwMxPf7qbakpNrr5ooeFy/sxXA7QNj6BJmNpDb+3PilXWsmrKclhIpcbzUbAy202g+MhdLOvT3JLFgyyXvavYEEIIIXsFPTbbzVLu3dTVbE/u7RbSfmon3xIa5lk2gWJT6rvWgjGs4A+tMWuVKRX/7Nlu+h3X/tvuj8tunZfp4FJsOvTajOOSU/Prdt3pcwCAT6NT7TSlJhc8tXXS/Nnrc97wzXiS7UkvjUyyZ4UY1hL1AbqaMyRBX9oml7a8+1vDuAOahcdSYWqnk+9rI5SaParY7M+ZEkIIIWTn2TnFRrKySuqyVJ1l7T/hjHg6Ns+R82aIYqNFIli2ANnA9ew/8dSndtPXvrY/7kteCgAIF+etWA5yubsUFa4+t81NN3U/8wFRATuPjRxaYcFb4/HAlGRJLaJJ5qLxhADK+SVZQctjU6KWv8ZjOJnY9KK9qy0WniE+GY8gtVblhooNIYQQQsj2sbOKjRX9JPEk3Z2hJT7QlmvzpmTZNXHtAjny2MjfHkVFqjCaggMsNl6t/cuG7cc+1k1f/epumnKXnD6dNce3CHptfFjXJ13Dp1zVrXP2bPf70t5qZeWxkVFQhwf9cc47FBVPSu6a+mLlpPFEImlKTen7ld59aYyT3wZLkpVshBHFT0tx135qNA8TQgghZKfYo4rN/pwpIYQQQnaerVZsll4BnVo/HLK/KcrgGdBS/i715VXCvp+8sLhoiHlYKtbSO2n5HVv2/1f+Sje95lQ/hMLpzjhKo/D+Muum6h+6G27o+p5SN6U0CAOLhuLZMmuYhJqZ1+pD1ULDS9to+Q9Kx9G6oDzfoHQRWrbxMMTd2zNFj5cnCGEretb2rCtqf86UEEIIITvPVis2OUNy042qwNaq49viLPM4p2vm4cI2Mju75rH2JPWzvIwtik0qZmq8vv71acnOvAYAaCLWaLoevexyuhsDdSGpYy5oSqVmdhyPk72WSM/ajxb2DfhNw/lxWpQaeUE02aKUTVP+9kggS/5GTiE0eaLX15B7cM4eKTa79UUnhBBCyCKs2Owu1viQ6n1fpaIy5bGkT2bsg11Rd/KWSvpbtl6scfZqCa5KrSqtAWopQs99bje9445u+oIXdLG8e/cy7AlDlKukeqVH/NSp48vzR3826OWBUE5aPDbyQS4NqdAiT2rqi5U7QVN3WpBh36UP7gjfzBDWZW0k64PfckIIIWSX2TPzMCs2pEyLtOV4YWoKDqA3Kq1GpsRqxJ7rxjfEL/1SN33PD3+8L//N3XFgnPMWkkd30W/TRrpehwfzOUDZVjaLgpoiBLBk1ND8MZ6xRmrrWsexzCJJDZYKjWTjRoI8jleUKq1X8+WUlGWyGlixIYQQQnYdKjakmZbWyMSDvk2CEeE0a51VXgxPIEeidAlkq6bkw5HbWlnfJVdf3U1f85p+RjLbbGBLkoxnSpXK8tgsDHJZ88SU1vEMj9AiZXq9NaVnX76cpXVkeJiMkvK8U1PnvHEebuw2/FxsPqzYEEIIIbsMPTbEZMr8Ndve8TrgRbG6+OU6Ek8/t5YqAwAefLCbvvzl3fTCxa61eeLEbnlr9plleYlSxFMxR5BXqfGEAHr8MrWH3lqmKThWWaZmRyUPSwRLrFWs36OKzf6cKSGEEEJ2Hio2hBBCyC7Drigyu/8exXQTjcAeSrnhB5Lk+RbzcHE/sfzbUs81Suuk9Pg33NBNT6JLghZxsr7DLYfDLEzDglE4/3tIl1GLqdfbvVSi5QXU9p//x1iaZ23bgrXNjnZjkWlhxYYQQgjZdajYkEnZdpOwZMQL4mnQjWmUyW3zoqaQ3aTYxBO7r9SQiRlj7i1tI39beRBqL4/HnKwl+Sth7V9+A2SivvT7RPZfzBaoLVtQxOHsUcVmf86UEEIIITsPFRtSZsm1+xafTMsQCpJ83ME0pMIlJ9IOu3PMhx4g28nKvEMe9aU2X9uPd9shYdnbIEWsKBue5zC7JrLvm3l4f86UEEIIITsPFRvsVUV2HDKCagsuXF7EZz6zm1446hLynTxghBCZkFrU0tBoHy3aaWo1oyWEUa5jDZq7o7SoOhsRPLsF3+upWOuZhhBeEkL4eAjhEyGE16+zLIQQQshOkrqipvznOqz9f3zo+OF++e+EEL5iitNdm2ITQjgE8CMAvhHAPQDeH0J4V4zxY9Z2e1Tp3BxW1BrzNkCHNlTvuqubnnz8EQBAPHPFsB2R/aXFYzNm/63LvOu2JINKJGli6u/AzhlZjrPjp1fF+X/8SwE8p//3IgA/1k9HUa0mhBBeF0J4ytgDFXghgE/EGO+KMZ4H8NMAXrGE4xBCCCH7zeoVG8//8a8A8LbY8ZsArgohPG30qTrWuQFdTetne1lpqhCSpwO4O/t9Tz/vGCGE20IIHwghfODRR/9gokMTQgghZATXpv+b+3+3ieWe/+Nd9YBWqhWbGOPfRicT/SSA1wD4LyGEvxdC+KKRxy5VkBYsVjHG22OMt8QYb7nssutGHpIQQgjZPyLCpP8AnE3/N/f/bheH9Pwf76oHtOLSk2KMEcCn+38XATwFwL8KIfz9Ece+B8Azst83Abh3xP4IIYQQshl4/o9fSj2gah4OIXwvgD8P4CyAnwDwN2OMF0IIBwD+C4DvH3js9wN4TgjhCwB8CsCrAHx7baNtyDO1c6zIBec1hg81kH/hF3bTC6c70/DJvmHABH3EjfXwTfFxsvbf8uBr6yYngRz6wGJZwQM7HiKen16eKHRdrOH/Ts//8e8C8LoQwk+jMw1/LsZ439gDe6KirgXwyhjj7+czY4xHIYSXDz1wjPFiCOF1AN4D4BDAW2OMHx26P0IIIYQsEuPqKzba//EhhO/pl78FwLsBvAzAJwA8CuA7pzh2tWITY3yDsezOMQePMb4b3YmtFapATmSWqS24cHkRP/nJbnry4PiQCoRMQlJBpFqizS+tU0KL15g694Xcn/V+e85px0mKjEf0SrdwIxL1rZDS//F9hSb9HQG8durjMvMwIYQQsuNsQTt0MlixIWWW/BZ4GnppHasrvtZwPJE94Vdd1U2fuNjt8NSpehnIdpB8UksfDLPFY+NRaLR1rG3TC3HhQn1dz/42hVIZl1Buz2F23P6z87BiQwghhOww6/DYrBNWbFbBrlX/R7whmv2gNm/o/vOiPv54N/30p7vps27sWrzxxMnhByT7RYsKo01L+0m/S36a9P2ovTyl/UsjSPqdzB5W6E7aR3qJPCqS9ru2/YawBUUczD5VbHb4NhJCCCFk36BiU6CpZjvZCBMrJp3kBDb95Gs4OJhfC6uxqu4nlH/LxuxQS4FUbG68sVNq9uElWLr3ZE+IB50SEfIHLCkdLYpNQj7Umjpj7cdSe+Q6FrIsllKjHdMjydaYKpcPOQYVG0IIIYSQLWQfGquEEELI3kLzMLGpSaEtUqkM29w2BrwpUk0vqefaJbT8iDWVHgCuvrqb/uIvdtMXvbDvnunPI3UzkO0lHx5jyu63J4+OP6h5t+usWyrlFpBdOJZ5WAvdLq3b0sXl7bYqbSu71qb+H3FHu5OkL7t02dbpXNinis1uPmGEEEII2Uuo2ExFSytkEw3HqTpfqtZbTZAMTyMzYUW1yt/Wti2G4gcf7KY/9EPd9O/+T7/b/XHzzfWNydYxZdK+8+e7ael5PnGiNxTPZ3RTTbnJ10nqiPxdeplqbvr8BUrh3bV18+CBVF6Z/98yPdfKaLHiYRmm8jFvo+C0b11RW3iLCCGEEELKULEhZayR3WTVf/Zb96ikhpxs6OWNzCENU4mlFJ05002/6Zu66fsf/mIAwAv64p/YsWo+Q7yHk9Qe9VHPSCLLyRPKw+fxzWgPfT5PSpgy+Z61f6nQWMdJeEZ3rCUftfa/AXhzp+brWXaonNJ8zyVdFlRsCCGEEEK2kL1TbFqEiBmrbGlMeSyZfG9slb3SfC15XzSFpuSfqSkypRaS3I8nYOR3e2vN85/fTU9efKz748RpAMcja8j2M8Rrk9Y96qOhpMcmf65P9iNypMipQ+mxKfnXat6ak9kwH+mj5U0EmBewpgzlL9uQZJ2aMjQkyd6SvrNT5wjcRvbNY7N3FRtCCCFk32DFZgup3bTS8lE3esp8NutsHnhMBNo6Cy3ReetPKjNShWlJ15G2LQVctKT4kI3kN72pm/7g3xCD/2059NaUaVJu+ofl8cc75UQqNieML+fBqf441sOoRUXJ+aXtZSGS6SMvVNpe+m80rw2gRzSVvgmeMMf8eKVtWkIml/yNnCJAq6TmSXZtPORNZWcqNoQQQggpQ8VmS1j6jZq6lTCms3fMybaoMIlSf7tWhn5+S2bglkCR1BCVys3Y/f/oj3bT7/3eKwAAl/WDZJ4+vXgeZD+YqTr9Q5cGTL300m6aW180ZjYZzWsD6B6bkmoi1Ry5zczk8+TiNmmeLEvpOFaUlaTmk5nq2znCj6MpJ0PE9NJ3pOU4+1Sp2AS2umJDCCGEEBuahwkhhBCyU7BiswO0RC6qvTIlyVLqkJortjZvSmTq9mXtX/4W8y25dkhXUVLYH3/8+Pz8kieFXYtutRT95z2vm7797d30ta/t93XxwsJG2xACTtOwD8tEnOZ99qFunVOnuvktidXmUdlK+Degm4dL3UpamLdnGAate6k0pIKWoM/6rqQXT4Z5W07/KfqINoAhXVr7VLlYJztbsSGEEEIIu6J2jpXdzKlbG6kFJFtjU5+Q3F9qpVmGY03BGRCGbTXotJxfpW1kQ9daNxX3gQe66Rvf2E1f90X/tvvj1ltn28QTDrco2UqKys25cwCAe+7pDOVJsbHCuxPp+VoQWE71g2Ra8qFUbkqyZJqnhXnnhUwPuSyMZiLO15FYclUtIZ/HgDzxoMDz+3l8vy1i+hjftPx0e/dLpmPnKzaEEELIvkPFZotZys1blhqz7JBIiXVxhiwTCs5B9jQNUWxkd71UbkqNWS2ZX8nWoOUdu/rqbvrPz74UAPCyR+frXHWqvO0mQW/NOJ44P2/ZP/z4FceWXbgg19bRwoBnz21JUamFcANzxUTKkvJ3KUGfFnLukaASllqsvdBW7odalro1ZLHTxCMrel2zMSVKYtg6Kxf7VLHZ4M81IYQQQkgbO6fY1LCipRZqtJ5ogKHLl0mpBaetI0kXyOOxkYrNAI9Nfpmk2pIapFLgKtkPZIPXugSp2GfOdNPLLuum731vN/2OP/VItnI/MOYBvTa7SvLRAMBDDx1fll4Hj3Ij7SwLY1eemj+4QT64LUMqyBelpABLpSbtL61biqQcopRoL7Qn7FEz0ZX2PwJvQr2pjiPvv1WWVbFv5mEqNoQQQgjZGfZOsckZVIPVErRov/N5Hgv+OpUeDWfkVItiI/00wLyRqTXoZEO1tK7MJG+dTpqmPDm/+qvd9Nc+PPdZfO0Ln+jK2/tYnjwKx463TuitGYa8bp95cP7OPfhgN9WeRQvtWU9qT76Pk/KB1aaleVIaKuW+kUqNVHBKyqx2kpbMrb3QLRfOw5IiqGTxPOl4LFGqtC9gfmta8qtNDRUbQgghhJAtZG8Um5agn5YsozM2oQmv4Tn52qCY1n76+XlLOGVdlS0f6Z8pqS9a2o6SYqOpOdYpa6mBrr22m771rfN1v/a5neniwtXXH9vHJt9uYjPLX3Ou81L93u/NFbqkrsjgoRbFRmvR57aZE6f73DZDFBs5TYUthQtqUVGWCU0mYEnrWh4cT6SnV+oYQuGF9ySF9yzzIk9vkxSSffPY7E3FhhBCCNlX9qliw3YnIYQQQnYGKjYeSjrlmCRSHil2SWa5KlYXlHTAaeHfAELopfaKp9AyHMtEfaVuAembTEWzuhOlVCxP+WMfm//9Q/+064L61m/tft98c1+mo+NjN6xykEyahttI1yvdo4cf7ubffXfXBXX27Hzd06fL+yhFYSe00U+0sSvzZSen6Iqykvp5TcR5Aa2T1fBk4Gzddkl4Pt1WJLo2rqi2j5xBNoeJoGJDCCGEELKF7KxiM6R22pSgr7aOlVJ8DC3ONGudWtxhvq1mLDYUm8P+/GtKTSmnmJymhmgyduaXdoh5WApQch9XZFn1f//3u+mHP9xNv/T8R7o/eukmXna5fiCyGfTKQ+ilmrMPXQNgnowvT9BXe61Kz5cUNjzP+lzN6VSkwzGKTalQ8sWT61gvSm04hhKek15WKHgFjzA+RZHS6Vk5TtcwWgQAmocJIYQQsmOwYrPhDIhELv6eJEGfnD8kLXhpXAHZGTtFDGGpg1fKFy0qj6HY1BpwpYapFg2qKTh5sYd4bDyttE9+spt+//d30xf+uy8DAFzfH+cS0bglm8PMXnLUPTCPn+iUmv/wH7r513Q/iykHxlAbZQBY9N8cnppAscmlI2k+097z/DgyBHwKj01pJEht26mYfYeOyyMth2nxArac1j5VLtYJP8eEEELIDsOuqD1E3vAF+8kQFaalY3cIq8wQJ+WPmucGwEGlcVlq0Mmkfdp0qMemlCitRMkW8IIXdNO3v72bfvM3d9Mvf34aHXH+Ki0jQoqRUO2k+5i8NEl9u+qq4+tN9cHXhlAoPXfp73nEn+K1yf/2ToHFpH2eIRWkmU0bwkFuJ5eVfpeWyYjPId9ZA4+SYqlqQ/drBbIysedqYMWGEEII2XGo2Gwx2s3z5A+YZCiFUnV9SE4amUihpb97Ciwjk2MYhprHJl2K3C+j5bGR8z0emxKpeOnYct3SwJzpmCly5lOf6qa3395Nf+zvnev+OHNmvtGJk8eOJ/e1j8hcMlORrvHsucqVrX7hHXcc91qk58cT7ONBUwI9eWzktgtem/zvIYpNOllH7qkFpHLjiYqSlF6mGhO/KEN25/HRaJ5AGR2Vs87KxT5VbPb4U0sIIYSQXWPnFBtCCCGEzNlE83AI4WoAPwPg2QB+D8CfiTF+trDe7wF4BMCTAC7GGG+p7XurKzbLvlGzUYDzmcsIWczl2ilybg/JAmVlldKWGXK2twvKStAnfY8t5uESJQ9kqcwlQ3Pa/xNPHN/mZ37pKQCAV75yPu/o/PF1GAq+PBYfvfnNu/PO7qad7++H5xXVHnGZesCD1Z2RniNpOJ6ZiEtZA2XiPKsrSoZuy33IfpN8HbmsNnZAaVnJjVvra26hpUvNUSRZFI8BWMvMIY9X2pYAAF4P4P+NMb4phPD6/vf/qqz7x2OMZ5VlC/CTSwghhOw4m6bYAHgFgFv7v/8ZgH8HvWLTxF5XbLxJ/ZYW7t1iJpbNDc9TOqZ54BmTwBoEszdxppTx2mCYQzySJfNwbYSI/JjaqZVa1rUypAEzU0gxANx6aze98cZueslBiv/tNl7lgJm7hgx7PzrqrmW6/nfdNV/26KPdtDYyQItIKkWTEtIQbCkFcn9JXTp9eq6khLQwqTg15QZYTMjnGSah9hLlF0o6ZOW3zFJjli1b9Od0cNCuXGtR6y15BtM0v1wtn+0t4toQwgey37fHGG9v2P76GON9ABBjvC+E8FRlvQjgl0IIEcCPe46x1xUbQgghZNdZksfmbM3vEkL4FQA3FBb9QMNxvibGeG9f8fnlEMJ/jjG+19pgqyo2MXa14JZ+bmtIBWueG63fuLROi/dFNvO0IRbG0jKYZi3MuzBeRQjHB8OUCfry05CnXPPalJZZIdZyEE3Z2Cw1MtO8pNCkRvPp08fn58f93d/tpikR3BWPPtj98VStQUK8zHxv5x4BAJxHN2Lp2ULvu3zGWpSamoqTh3bLZ1obFDOfv+it0dc9KR/8mnID1BP0ld5Z6aWRL2IpQZ8cSVaeoDXPkypjAixLjxbOLdfNT09LwGj9V+CxKS2bdahFMcZv0JaFEO4PITytV2ueBuABZR/39tMHQgjvBPBCAGbFhlYmQgghhKyadwH48/3ffx7AL8gVQgiXhxCuSH8D+CYAd9R2vFWKTQtDaqfaGI/H2iJTJJnyNCGG0NLicRuMDGpDLUBPsldqBGpDKmiZ3vNlLYNfaqpOKUGfLEtSalKj+dJLF4/z6U9307e9rZu+/OXXAwBuuqz7fcVli53v9N3o5L6aJ8531+nhxzulJnmc0n3J760mEFSHUCkgE0HmPi8td6ZUX0qvufTWlBSbE73fJniVm3yeNhhmaRvt3Z9iMMwcz1AKGp6bNYBasS2VR/vs5tuky78kUcrFBvp73gTgZ0MIfxHAJwF8GwCEEG4E8BMxxpcBuB7AO0N34U4A+Jcxxv+ntuOdrdgQQgghZDOJMX4GwJ8ozL8XwMv6v+8C8GWt+96bik1LRb8pKkpW2z3rDkEmTVhSy2VGSfoY4LFJkQlaiyhvMNYuqTU+oFa0ki2gFhWVKzbSW6NN8zKleTfd1E2TgpMu6fOvybqSb+h8dVb6/V0iV1+kSiXvz+yanjs3W+ehR684Niv5mJLiYQ1+2mIn05DeGGD+jEhhwxIm5oNfHl+W9lEaMHPBa2NFRckXQwtHLElc0oTmUZ9lDhxr3Smxvj2Kn690H7QieqxCWrqfvEi1nDfLZhMT9C2TvanYEEIIIfsKKzZ7gjT2a8uX5rGxmhCy1SSfyrGdtbJ8nrw1Em1gvexv7VSt3B41b43VMJX7ylu+WgBHopS9WNoZZFSUnJ8vu+yy42W7//50Xk+brXtjNnZmvu2CarHDyHuXsvKm+efOXzFblok3AOyouDGviFe5AebPp4yQswIZNf+NpdioXhvrxah5bfJCyodNZi8ueW00pab0bfPIIV4G/C9dKpKWbsdjFRqyDVkNe/DZJIQQQvaXfeuKYp2SEEIIITvD3ik2Q3pcinjiAhOWsbj1eJb7rKarToWUs0tUDHyloibFWzMNlwzH6W+PKU+7DLJ7LA/l1bqitCkw705K00su6aaXX764rgwjT+bac5/vCnPyRN/dV8gINx+kNRZ/rwMtbH1WpnzYjf7CP/pot80Vp7s+nTMnj3dvnD81vyGpey+ZhWUPi6fLYFbWQi+vd0zFfHm6h7mhuHR8K0GfpyvKbSIuzUsbW5nn5EWUfcEt3xGrP2bIIL01jG5wT5ek1kVe+nTLedpxNq0rap8Um72r2BBCCCH7Bis2W0jNCDxkXzM8Ve+WVojlMpNxgdrgcp4T1bYdSi2pX6np21ML4c7nyYaiDI3NW8aaeVju07OOJ9y79jv/WybxS2pDUmdywlF/Px9+GABwJrlknyhs1J90Ui0u6Y8z93GuXrlJSo30kqYynjrVlynNAGZO4Cv6c55doDNnjv0uhdJL0601HEbLN6E2BEvJK1tKwJdTatFrY1R6FJuqibg0T74oJdVVztPCv6dmSSbi9PwfHoZjh7H8zB7RWwvvluvuU0Vi09iZig0hhBBCFtk38/DOV2zG3MzitrVYvpaYwhZkM2AZ/dRLoOat8fRdSx9F7oHRouHl8T3rlIZ5kMpM8suk+el3mgKLHhvLjzNrZMtmeSrsg/0AmknVAObKxhd+aTftR4A8e3QdAOD6py76cuKJ7gSkHydRUnc0v8w8DcJcCkl+mQcf7La5Dn/QLbiyK1P4z3d2v/N47Suv7KaVARpPnJg/65oQUVNNLKyBcYd8P+SrWiqblpDRo9hUvTaleTWvTT5PvlQlyWlI2osp4u89Q8FUblrpNLTPerokuUqsbdOS15Asl7Vc+hDCt4UQPhpCOAohmMOeE0IIIWQcR0fT/ttk1qXY3AHglQB+fFUHtFpltVa/WfXWIp5atrGaEKkZIJt4Vu74qaIXNLTEfIVWU1IC+kHMXPaimrpTGrfP05CTXgvP/jWlRk5zxUb6bpJyU/LjzI71qFBq0v1+9NFu+thj843uvbc75m/9Vve7V3Cuv6Mf9PbFL+6mN9882yQkdaQfg+CoTzt5eNQ1RZ84ymSwnksOumVPHnTLkkJmfvkLAAAgAElEQVRz+NBD3Qq5inTXXQCA637jN7rfz3tet49U7jT2wR/6Q/Nt5AVSnqeDgg/L89gOSZhmPdKl9Vr2X/J2jFJsTosH1qPYSOWm9DJJ84gngqrlf7oppAyPpDaLjjo+rEupKN6ptp98/hRBsGQca6nYxBjvBOb/2RFCCCFkOdBjsyNM7q1JtDTL5Dxt27wlJEOBEkOSJEzh7RlLfzEPD8uDYeanJRuVMl1+KbN7qTFZOPyx7aV1QCo1pagobz4bYK7eyG3T/LyRPIuGStFC2jTv5Jc+q6SgJFUmKTfve998m7R9r7Ydfu5z3e9rrunK9kA2MGeiH6DzMPl8rrji+HHzk06KTCqDjHSKhXw885Cp4+clfEYeH5ZHNdEinKw8NkOiLEsjD8iyaa37FsXmyaNu48OSYlO7QJ6XSfPa5FgJhGq0SGgt6ykRmda3Qvu8WoqN9iyWgsg2oVKxCWVYFUur2IQQfgXADYVFPxBj/IWG/dwG4DYAuPLKZ05UOkIIIYTsIkur2MQYv2Gi/dwO4HYAuPHGW+KyUrGY65YaFlYCDcCOALCUFM2r4+nTLskg2r6n7AZ05LGRRSi1nrz93JbHJrWSavltSmVKZcmjrmRUlGwcp+X5eUiFRnprcqFjplpozfJSE14zh6UTSfvM77dUTFKhkrKSCpVvk5Qgmd45nWz+DM2S6lxyvExWYhhFoZHTPGJL5iRJtAgGUoUpPSNpHZm6xcpjIzNXy9OxImvkpbAUG7nu4alCOJ/mrZHfnvxlkj4+TZIooV2w/BmZIpLTY4JaWFZWi/O/5WWxbJCJWj6bnE3w2+yTYrMBl5sQQgghZBrWFe79LSGEewB8FYB/E0J4zzrKQQghhOw6yTzMcO8lEmN8J4B3ruPYg/DE+lnLaw5HSxtN0m6L7in3tc7oMxFymbDUbU0BLw2pUFLUgfJl0jyQsguhFDWrdSfJbqfSftKydAsPDzK53tsFVeqbqHX3WF+fljBd7fkpuSNlV4RM058fL91IrXuhsI0M3R0z8GBLCgirR00+p9pgmDmy10cbDDP/W07TNrPhKkoPrtb3W3oBtSR+1jdHXphlfWtqz6lxMw/Et8ETdm+NFer9nA/JZ7hMNr0yMiUbcLkJIYQQQqZhZ8O9W/Am5TpGrdpeMsqNcTxK8uZAal6meZpJr8WsPJaKedhaVjPy5acnW9ZSucmpJejz5DmTIejyN6AbjWeqTj4Q5BDFRotJ1n5byzxShKQUA117iUrza4kei/uafiiRfPfa5XHkgTOVxXx56Tjyt3W7ZW69hSEWSoWpJeoDFkeb1fIilJbJExuDFYTgeY7FPO07Iv/Of3vCvb1h33lx18W+5bGhYkMIIYSQnWGvFZtaLXq2vEXFaFFFSs0BLc7QahnVwrtLCo6ct2QfjuZr8bSeZFd/3jCV82Rr1hIgNA9PKXRbem2spHsyFHxBCXo0ayV7w71L3pS0TCof0ucyFrn/dIE8MkbNrFJaZqzriTzWdp+wiuQZLaSG9gyWcuFpXpv8eZLz5t6a4/uYDbGQr6x9a0ovoOapsV5WeWHkPjwh3vK5KtEieYgypXQBKVXAsWWVhI8JS7HRLltJoVun12afFJu9rtgQQgghu86+dUWxYtOKpxNW/q5FJ03tfdHW9RwnMUUyLWcR5DLvFFgMHrMum2YZkC2vPEGf11vjGX9wNnxCyUBRU2rybbzelLyVW1NzSi1szUdRYkiEluXv0eYrVraWIras41FupHhRsq9o20i18MnCIyL3J5WatE0aYgEADuWOa2E/8u8cKUmUTiThUV/kumkfHqXRowQqN600GGbNz+cR1RPWM2hdQjI9rNgQQgghOw4Vmx1kZTfV41WxFBytj1o2/zxShNa0KM1rGeahxY8zK1v74Tyn4zlVua78LZUaK4+N/C0bxvn+FtSckl+mptTIXC/537Wp1QIuqTryt9bqlq1la/8SK6JK/l7SS2sJW9o61rraqyo9NiUBTYt0yiOqKmOEFoPgDqVSU/Pc1JZp6y7jHpWedW1siwZF0HPKHuthWkfmHvKMQEHFZjXsTcWGEEII2Veo2OwwLUEZqZEQMa+2hzEeGy2dpafVJJsDnnwjVvNjSHjJki39LRYkuUyeqtWglJfUSsYqU3toEVSlPDYLas7jBb+MNErIZnlJHdGa6lZzsNSs9yL3byULGuOxcTRntWehJehKzi8lUNaUGksEk8W3PDYyGkpGSZXGptTEPOm5AbIIqRbDWk16KL2gy1RuSnh8WQ0em9plaflEa8vzoizJulhl38zDaww+I4QQQgiZlr1TbAghhJB9Y58Um62s2KwzyVEVS2tsCcPWxhHwZHoaMjJgTV/1bNu6rLKJPI1SyKQsruwOsC6t1gXl6VZKv+Xy0ropQVjRPKx1K1nmSC3k1ZMMLyEvlByWoxWvK7LUl7Nkk7A23xMpLLugPEVN60rzcL5t7XZ7hlSQpuL8FqTu84Wuc/kAewIYLIf/FPfOk6agZhb2hHun99HokW/psdO6rUqpJzb6/6wdZCsrNoQQQgjxsW8em52v2KwsvM5jHpa/PaHVcn46IatlrUkRuZpUU2Y8MZFraoZYh/UkwtIuT6kR623JlYZhmJWzpTmuNeU9rdhEySisSQ5SQcl/j8kDv6IXz3MYKyxaW7cmglmm9IQnQd+QIRVqCk6+bCHsWxag9E0YQouCI9exjltTaKwBX8V0zNAKHlHdI1zvU+Vinex8xYYQQgjZd/apUrV1FZsWmwowLryu1Fg4rFXlW/wmVgeu7LCV8YKlpqoW89yS2c5qmtRYUqeyR+AacjiP+qKJVNrvYpmkt8ZqXWoSQSkmWbLqkNs1fClrhyz5WDSsxr4Wfe8tR+k4pfBfeZySsCbnyWEXrMdK/V61mEcS+UPuVeQ83yl5AtbL7LlgFV9OCvvOD6VlxrBU3JYxQz2ntmz2qWJDSxMhhBBCdoatU2xyxtZ+a+MBNmFFItWilCxTh2esAMkYNckzDIOnA3qDqQlnHnuUldRvUbFpCMOpeW1ytGUrjDxaF548bXJZLYjMWuYRqWSrXNpb8m2lt0YqQ5YNS1NwSucxa7qOeWeX9X7LC+Z5xj03Rns4DMWmZiNssSdKsT1fd12v4b6Zh7frfyRCCCGEEIOtVmymZlCN1qOKyOp/y4iNpep/a5lKRqNahFP+e0gUlGjWHBVSuEzJGEHLY0WqXaaSQLdQptLJe5QZiddjY0mPLc/TBqONjWgtsxr7LemEakgVpiTM1gS7fHu5P4/yVP0+mVKjQW0dz/Ae6aWylJva+9GighaK7rUgeVTcluiodbBPig0rNoQQQsgOw64oQgghhJAtZecUm5WNnlrTFFs0R08YZcswCZ5Qbk9Ip3asllBxhZUlTpyIIb7p2VAKljSu9Zd4hkcY0gSb8sJPpatrMbcNeEy9nt81b2rL5ZMjdefbyl4YuW4pur+lm2w2FIQcWiHhuXctw2toN6DU5Zn2K+OkSyciL5TnxtTCvbP/9byJ+axPp9b1XOp6XCdUbAghhBBClkQI4dtCCB8NIRyFEG4x1ntJCOHjIYRPhBBe79n31ik2Bwe+hsTSasi1GD9rmyGxwmnakt6+Nspjvk6LFDHiPDx+15aBBodQM/2VDMLe56jktZ7hiU1uOWkpG6xa/vJcFO19aHl5C89Vzczb4CEdFH3v8VzL2zHWpCxFvNo+8nXTvIVEfYmhH8oxL6kccNUyD2sfBS1RXz5Pmc6UVMyHV9DUl9Lj6h2fOH8Flv1t87CBis0dAF4J4Me1FUIIhwB+BMA3ArgHwPtDCO+KMX7M2vHWVWwIIYQQ4mcTzcMxxjsBINiNnRcC+ESM8a5+3Z8G8AoArNiM5ViXr1zoif+trVtCkxOspmNNqbFCuD2dykOSBQq01qdnm3XS4q2ZsQkFl9Seo/wj0+J50dZteR9GsOrBMC00kbXFy9OiPA2idO0neL/N75NUZtI68kLlx9dkNivOv3ZRs3VTsj5vWofaPGt+bRkp8nQAd2e/7wHwotpGrNgQQgghO84S2lnXhhA+kP2+PcZ4e75CCOFXANxQ2PYHYoy/4DhGSc6pNou3umLT4oUYUlNeWoPbag5ospzcpmUQzCHWfmvUt4b9p6gMrZVcTCrWYzXKhiAbjFPQEjjS5K2ZtHk+ER4FZ8jFbQhJkYkePc9TbdxEK6BG88t4sIJ9al4e61HxBM7JearHZtlYJyLLYp3IkMyJAxQbj+8uzbtwobzNfN/zv5u+E9vD2RijavoFgBjjN4w8xj0AnpH9vgnAvbWNtrpiQwghhBCbTfTYOHk/gOeEEL4AwKcAvArAt9c2YsVmKDVvjSeJQWm+ZrEfUjarTHLdllHfGnw5ni5xSS2IyNM4WyeTNIqnaGG3hGVYTUpZhpYhQawyTYhl6agNrWAF1HjSCdUovcqaIlQ6j1qU1VqfeY90Zs0D7MjP2o1okcNKio34H9Bjmax9Xumj8RFC+BYA/wTAdQD+TQjhwzHG/z6EcCOAn4gxvizGeDGE8DoA70FncX1rjPGjtX2zYkMIIYTsOJvQ6MuJMb4TwDsL8+8F8LLs97sBvLtl31tVsQmhqw17UsZ4lq+837Mll4RcVz6VLRFViZIKo23jaaI0KDZTRH+UbEU1i0qOp0t/Cpr2N6R5J18A2enfQotSY1F7tj3PuOO5GjNIpedZ0YQHz7OoUTr1WrSVFUGlUfIKrZwx0mzpW1e70aUb0fCQaGpLy+dQY8mBgM1sWsVmmVA0I4QQQsjOsFWKDSGEEELa2GLz8CC2umIzpDdm8oN6DqCtU3KbDcn6VFt3TPdSaV5Nt83+TuG5CUtJrnUDlCT5ZRgpW8zJZqLBIc9KC9r+W7qXpJRf6ufVXJFDzMMjr0GtW6ZhLERXdLE237rvmu/f0yvjiV72JBhcOOaQy66N4VA6uGcciRrSRJzP006+9EA0GI3T8Aq1oRVK1IZWyEnLLl5cXEamZ6srNoQQQgipQ8VmC9m4ELuWlm+J2rqWO9YTh6gtK60rFQDPkA2KedgKb63l4rIYYu7UxIsT2VvRYk52M8RVOOQB9zjlhzjxS9toEqn1HDsdm7GYfLTDo3BoSqD1jGhJ/CzFqGZOt9QkOb82r2V5sXAtO2wxAI8xC6dpi9t6SLh3Yb8pUV9CPuJW9o4liZOTsm9dURt06QkhhBBCxrF1ik0K+d44hph6Sp2yssNWayK2hM9aLWwtMV++/zQvSRkOX442lIInDbzWWra2SVit8JoNQIus146tzV8og/VsSCXF08lfu/5WZ3/N+FE6nqbQtXh5RpjeLA+M5sMqoT0bLZ4tjyAhrSKl17C23yGt65aBZQdRklflMo+SIpEvXv6MjDFKNazrHVqhNK82tEI+b51DK1CxIYQQQgjZQrZOsclpsSqMCV4azRTRMUOiojQPQ2leS1SUYxtvl7jV0JJ4olisVqtsLXk8PGM8Nkm1Cp77P8Zj43noZViGZxt5wTwKjebHssrkeQadWM+TfBZLYqim6gwZkFWeujXsg0WLKlXdicWQk2x5Wb37Ld2QWnRUycBUS9SXoT1ya/1/Y2L2SbHZ6ooNIYQQQmz2zTy8VRWbELpWUItVYUxjefKauCf1vSeiCTj+lNaiV0oXY0gO8YZthnSN13wT1jaetPPS+yDxKEIe348sy6zhbvlX5Dolr5NcpxaBVArlSPNSIa2HXCuDx3jgUV+cz5XHh1W671q+lxalbkiEXsJ6JjVvWMv+l/YfleeCeVWRsXKoFjmlyW/avMp5pMEwPSmbvP+3bIuSs4tsVcWGEEIIIe1QsdkTrBQeK2GEl2BMCxjAYsva09pv8ESkjMNaA+6iWA4sNsKSsGU1Aj0td4kmflnd9bXBCs3gj5bmnza11Bd5f6zQMM86ktr9t7bxzndi3at8vrWtNV828jU/i5VdthTcox1fK1PL7RmE68EduT+J96TGKkMtMrFS3jF2r3VGPpGOva7YEEIIIbvOvnls2AtICCGEkJ1h6xSbkooP2F5LYEJ50FvtHSq5y+1k14G133SSWh+b1QWiJeorrTvAPCy7njwKcuoGkNt6Bjr0oEWSWmG5WtdUSWWfmYenMNIOMerm23j7VErHmcJRWXoBnecxNOq4pfuwtj/rPifSsrV3cWcsPGotF9PzstbyOeQXzPsMlo6jmYhbuqKM/cvBMBPWYzvGSbAO9kmx2bqKDSGEEELaYMVmw/HUkP//9s42VLusrOP/63kex5FxQmMs38YSHMQQM7CJXsjM0R7ENCvBCrT8MPhBsg/BOA1kaVIyEEEGOaBQMKaCDUoz6YxYmMTkjKKmzijTUPg40TCFpUjoca4+nHuf9llnXde61tp73/vl/v/gcM5+W3vd6157n7X+63+tFelkWkze46qZ9SlSG0ux7p4BNRK7GFRq+osVRr1+LZ6/XGRni2JT6gRG813KEy4EQulzKot1jfVBvP3RBVgtSbS/XdONtVSZirx0RnTArwv9/Tki5uFo/cmJD9YEkLnFVccgLb4mVTpSYJ4cWXI/e8swdHQfJKdKlx68XJx8ae0U51i6tEKaxVbWoupshVU2bAghhBAS49DMw6tu2ESG+CPXl649tX/MyuH1gLM3RyyUN73WUgNy6UdUhYAXIrqUQs4vY4V55zpnQxSbzruTihle+pFZ2s/s8ypaNMx76DXpkgppWrlty7OVU1+s5Rdq/EUDTAuh72FHZDK8lgnzlkTRW+OpF5HVQksPQs7gVCpM7wG07pu7T80qugVyr0zrdU1VZjmsumFDCCGEkDJrbaS3sLqGzblzdT4ZT6yI3GvvpL3WSG0sKTUtUTiePyOg2FgduTTCyZsJ3VJLjhzPRcQykHod0miWfhrdvbprrEkEL7vs7H1ObAbpYpg59aUUDeUpNtZEff3C9dQci5ZFKmseuKASOFYU3JAXe2SR1e5YJIBxCJalLsfJse6zByKEzmxH1iuJqCPRL6D/hafvwfQhjfhn0m3nRWV5bHJYgntk5Zx9c2hDURTPCCGEELIZVqfY9PGCMyxqWuKjM0bCXk94qoiXYM89F70yxkzoXiewU1RSW0DXyfTsS50K0+Lh8TqzXZ7OLIbpKTY16oil0FhKTu7cliUVLFUpl4f0Wi9PI+AF7pQ+ak4RmpuIHS5yzQkRj01HySDnnRNZ46TbZ6nSXlRUjXqUk4UrGWD7WhxLqdv7YANfFyGEEELIMatWbAghhBBS5pAUm1U1bERs83BJAa/xz7pE3Ytj1aI0U7l0h4QMWy7rgeHeJWU6NRMDdpi3p26npk6v2LvslpZSaDFBRyJh3aGokmnYG/YppdH/cDWaupde6dxIngrXdqbrmgkZIyMUkbqyJKyvLGIaPnOONc6bO5buj6RTc62V6cjDFBnTtvLqmYcvnM5KjiUtmRGF5mFCCCGEkJWyKsUmpcVg12ICq1Jyhp4zBUN62unfmXO6nnUkLNfr4Fmh1J46EvEqpqQdw1SF6ftbO4Nxml/P0Gwt3nmSWL88030Rxca6psZwHKFkGo7UkUj6IzgzB/hD3fT2/chGiqJl1Q3zAWmVw0rn1IRhe/ut8G7rgazJ24FxSB+big0hhBBCNsMsio2I3Azg5wF8G8C/APgNVf165NopPDaRaNkwUzeLh/SEPW9H5D4Fb43X+SspHv2/rXDpSOespvi7+1heGy+/nlJk5qnGY1MTul0Tjl0jbZTqRmS+BS9P1rVJvfIsFx5LFEhTLOtRzTsuR1fsAkN6ihROjSpSE+adbqcPYP+DldLNPbRW/nOkHpvAo9TC3GoJPTb74S4Az1XV5wH4CoAbZ8oHIYQQQjbELIqNqt7Z27wbwC+3pJNrRZdUliYVps8Sx2lb1BdrXnZPISgoNp5qEZnzq+TL6Tw3uY5jurZjrvPX7UuXVOjSzQVlpEsqWH4fb5mHEwUqXVqhf9NSV7HfpbcUm3Tyvdx9vBU/02vSRTAjebImD/TSH9AdjggRY/tvLKx3y9AompJS4wl0QyKEzjy0NaqIJ7eV8CKcOrwvtcX308AI1rC9sqR/WVOzBPPw6wG83zooItcDuB4AnvjEZ+wrT4QQQshmYMNmBETkYwCenDl0k6p+aHfOTQCOANxqpaOqtwC4BQCe8YwX6LlzsaF9qwOZO8f63fV6TsapcyxZwYmsljdCQT2aRv8gPv1ETn2xrvWWL4gEcnTZTlUYK0oql5d0+0zkk5On7vf5vmTU3bzbV5qjxjsWucZTdYDTdcaLzEqJygoehhK4ViIWt9I13rFB82+l+4G29U+s9HMPuKWypEpO7j4teYtIdRNUMk+Y2pd6eOhM1rBR1eu84yLyOgAvB/BiVX7dhBBCyBQcmnl4rqioiwBuAPBCVf1Wazq5nkppPHtvs0YucWW9PjUF1eCxKUU0eRFIJT+OZwuoWdvR8vD061W6UGY6R40X3VWcgbh/s1I01NjXdKQFVjEzcOgB9GSL4MOYq1feOWPSZbemXrWk7xVTRHhN95/8XXqIvIe2I6KKWA9i7gG3Mp5TET251ktzinPIapjLY/NOAI8FcJccP7V3q+obZsoLIYQQsmkOqe02V1TUs+a4LyGEEHJocChq4Zw7F1tKwYqIzZ1TnJQpEpOcHvfYt+HYGzqIFFQBbzimZACOXOMtX+Ap31Y+rfvlIqAtk7CnkJeMxo+5zBkiSh3Nue8lrahWiHW/MCIG4BTrmprhsfS4d58Kpn50rNDtbn+/aEtDRd5jF9lvPZqet9scigLy27lzvfHd0tCWtyaF9cB5Q5Ol/Hsvko7AhH3/HyhSnhfkkBoKa2N1DRtCCCGE1HFIDbFVN2xqekLR481E3IVTORCnINAb16RXM1aUZo3Kkyo13sNbMoLmFJtOZek6k6VQ7n76Zsc3504uGXT7IeLWxHyWBJk7FmEM1aUlRHyyh7SdNEo+p+ikgkN6jjdNRUtR577mItaD6Z2b2y5WciON3LHug3jyZ7rdEi3g5WECIiISmYZVN2wIIYQQ4kOPzcIRiSk1kcjU9Jg5Zt3isYmEOXbUnDs1Fb3l3Ee1jkWG+qPqjhehGukgduJHGl2a6wRan8Pz/Vh56vb3lS6xKqq3Cl9J3UkrtPdhre1IupGHqqNh+YS5H4WxqFGWW+xLuTRP/j6q8MBYlTt3bkki9bw1JfoPYDrfgqXQRDxDE7P0+rq0/InIqwH8HoDnALhWVe81zvtXAN8A8F0AR6r6glLaq2vYEEIIIWT1fAHALwJ4V+DcF6nqI9GEV92wWZTHJjJ5co28MDWlQf2KLmONWtKRG+KPqiI5RaWmSC2xwhPbrHM6z003uV/u3O5YLo/no6pIxaKk7oMRUWqsa63tNH+1LMhTY33UNHCnXwcja4TmtoHymqG5c9O85QQ0MyoqEk0UWbE2pcYQ19Fl3KuDLe9IS9WZaIL7NFlvrdA5WZpio6r3AYAMXpn6LMt5oxBCCCGEnEYB3Ckin94til1klYpNpMMY6cRaNoYz1xxlejWlHtDQ5vEY6USmrB/QW4504CIWJCu9klcld6ymZ5R6bXJ5jvp9vDx5nVhTsYmoL6X5bPrXNnhcTJ+MdTy3ryWSKmHoo2R1CD0rUhoFZyk3ETyblGVXipzrRUedzMdSehA9+dPbbz14NQ9kWtg1yo3n/1kIQ2xGYzORefgqEen7Ym7ZLVp9QmQx7AA/qaoPicj34Xi1gvtV9RPeBats2BBCCCFkVh4pGXlLi2FHUNWHdr8fFpHbAFwLYDsNG5FyT6nGYzPqgpg1BpPIOHcprUhvt+sRXRj5azbu7U1QWop4yl1T6hTmjtUMo5e8Nrl7Rj5HVOUBAFxoiHCyvvvUAJKr7KVVJIeqMB2pTBKZ2XpPXpsx7ED977AU0NaSftYvY5xTdZ+WBy+STkQ5seperr6m9ymlsUeW4JdpYYGiVhERuQLAOVX9xu7vlwJ4a+k6emwIIYSQDdMNRY35MxQReZWIXALw4wBuF5GP7vY/VUTu2J32/QA+KSKfA/ApALer6kdKaa9KsSGEEELI+lHV2wDcltn/EICX7f5+EMAP16a96oaNN6xk+Sq9dEwpOTLOYB33zq3RNPelIzYMFQwxyUWKydseY0kF7ytN0y2FpOeOWfc5lakaE3HJTZqe16fkpPX2pQ/ZZOO6fjai57SMoJXeBZGRYCtSPzKDgvc5Su+yU19Bi+s9xRqHzaVjXetNPpp+2MgQ2ExE3nE175E5WGCxTgaHogghhBCyGVan2Jw71+Y/9CaxsnyaZ0Im+3+XZIUcNdeWei+5/WOaLz33opEFr/PXMjdhTfo16abKTbpdYwT2Zqa3PnNWqLPkBc+JWqPUWB86wgQTaE1FSzR7yzU583BISTHuXYqsz53jphFyrsOv7JZcmaPlvRiRv6zFOidSd7yFfS3WooSsJZ9jsLqGDSGEEELicBHMFeGJCpGJr0bB8s3kjCApkVnqrDQm9DSMSelhyg3bRxfDBMoRpB5pEXrqS6nj26fG1tD1EM3FMGukxoipI515zqMkEdRcMzI1yVpLG/SXweiwotRTy0ikCLyiqJk0sPRuq1JsrONAWR3J7RvybmupGwv472z5ZbwiX7B9aJOsumFDCCGEkDKH1KhaZcOmxWPj2Q7MmegjTW+rl5OjxWPTEjnV/fZUnVJ+Bz4FUz5EkYkAI9elfokW/4/X8U3PzW2ffFVppjxqIptKaeRMI2PQoigmeWiNhrKO1QhQVrF4diPrI0fu462GURLmqmh5X1nbuWNDZq/zJpacCqMCpYJT5KOT5bDKhg0hhBBCYtBjs2BE7HHqaG+sacg/18W2zskNwJamsS/t66fRfehczyj9cJHJE1qiuwpJRM71PnpaXF4ns1RcfdIebtrJzK08YQlz3sKZNdFQJ4pNVchLck7JGJJLx5KvapjIP+NhqRSR4DFru0+68kR6bvRfdbUAABEpSURBVO47tDwvngqTXuMJE1465v6jpPJZD1ek4nrKdendVvPumdM3GIz89Bgy7c8+OKSGDeexIYQQQshmWJViMwSvJ2dO4BrpbVg9oIhvxusJWas6pspN9N7e/iHHEOu41ygrwdueIrIIZnfMmsrFU1ZaZg5Nz3XTiCg0Hfvq2XZ1LHK/MdWbroACYlUNYwR51dzHK7YhAl2HO/F0S1TUEFmhpOAA80d2BkxJkVfnGJ68OaBiQwghhBCyQg5GsSGEEEIOEZqHV0DN5FiRqcrN9FMDXv/vllDIiFxbSuNMJlEeW9lYjfbW1fP2W8Xj1afSEJqn6JfSOHWsy0NLTHJoxrYEa+a5sYiE/Y5QLyPzGFrXeMdaIupLpuEaY7MX8WyZlE+WgAHKD0bN9xNZa2TIfTzGrpdAaHyxpWrWrEBBpmWVDRtCCCGExNlY/9ZldQ0bqydTaoR7vaVIB9ekJjSydC1wtgc9Viy1dZ8hbtiVUVJovJ5Wi/jV4uFeLdEPWxGzn/ueWvymJfW2xjzsKYElhSg3PUVk0r1B76moibj/d4vKMlSZmZIaiT8hMhloTRWfq5gObSiK5mFCCCGEbIbVKTYl0h6PNXeZx5lw75pYv8hsbJF4wRqJoCRBRLoJXgh65W1Lx2qJjFW3CFwRSum1WBay1MT/DolB7n6PUVB7KuyI2po711JkvO3Us9Myj6G1AG8unWge++daea3ywEQeKk/CtCr8vmWBmmch8PKvsVCWll1YmkKytPxMCRUbQgghhGyGzSk2HZGGfGR8+wxWk97rspcGZfuUpij3etxWdFSaRv8+NbNLncnL6d1e+Y2h4Izd4ygVV2teWqrGKMywxMEZrA/UsphrhpZIp44hK0/UpG+JbRG1OBLdNYrXJt2fO2alUdo3JtGVPnPmy0hlKRRm7jVpFWHEZjlnpNQhKTabbdgQQggh5PDMw6ts2LRMsV4zNn4yH0TEb5I2xXNN85q5HqIentwHsgZ9W0wqA6NXLFrmFYpEr1j0s7xvQaMpKmpIYa6BCr/akDljAHuhybGn+ymtPOHZpNxlEYy8WOk3zTOTY8kTsLTIVRXnttgfa0TvQ2pczMkqGzaEEEIIiXNIjaoVdvkIIYQQQvKsWrGpGV7KXRNWNb1JrEpOsj7WcE9u2KpFuyyF8tZI1dZ1ve1zgdozxOw5BmMZgw+pt3NCaTn0yLUtBZdck5vYLsKQJRWsR6lmSZaa6RAiIeg16Z/hEMdCIgV35gvPJzV0BYq5oceGEEIIIZuCDZsVE5kcq0hL7yZiMB5BJXFVGavL4MUfWnn0JhhMyIWztpiEl+KFXUo+JmHI282rv6VuqjfJm3FNZ+I/f/5s5akx35Ym84ssdRDBusYzv0cm5ivlxTUPl8gFLqTbS/6PmCvcSKEatJiH01dlTqxcooqzZTbXsCGEEELIaQ6pUXUwDZtcr8cLFT2FN9NSpJdTmrQsoqSk27nMdse6bmxNuLqXV8tjUzEJorWdC2EthXfn1J/0mrEf4kWqOJHvrnRsSGhvP83o6pQ1ymaGMT02OaWxtIRCLmslRShXNCVlyFM4q6Y/sGSFoQ/IlCHhrdJWNF3HY2N5abyFca3t3PWH1LiYk4Np2BBCCCGHCM3DC2cvPecxPDaRwVlvvu4hUVFW9yA3nm4N/nqRWjs6D4TsupDeMLcXnVbCWpCwn07NjP1pulPVqUWqPB2RHnxan2qI1ONSN/ZEGTwreViKR07pqPHLlKKiWtJouTaS3plrxvrHte//gDWzEkZWNbZmb6ygxp5oBbseUkNiaayuYUMIIYSQOg6pobX5hs0oM2/XKCo56aDUO474WmomUrDS95Z5aAkHcIj2WiPTzLeQ63F7vfp94n6ugnoRYqwuY+rnSuez6acfNRGM5LGxvsNcZF7acQ976zLHajw2Vhpenjo8y5J5LOLrW/J/t6hPK0erlLb7W3H8BUSqsRV02jHw1TkJhzYUtWSxnBBCCCGkitUpNlarc4zeffhm3rm5Jn8pgsrbl66+13Xt+gpO2oMu+Wf6pHmLzF+T3Of8rqc1ZKZVD8+v02XfWqxwqGJX0wu3jnkd0WL+InUkQk33cgxaFEAz+s4uwMj8NSnd9/Gd75w9VoquqxEGUhVmqOeq6Rma6vsdInuWJFknaqlq4qIRpN8apWVfj1YrS8vPlFCxIYQQQshmWJ1iQwghhJA6DkmxOZiGTU65DJnwrGM1wz3ptZ77LD0nDbVNh6S8e1txiNa+XF69c5OhgpyCbMnx3kzoQ8JlPUp5sc6L5Cnie6wyDQ+pTxHSehW5Np34MXf9EPOwlVbAYO6NTFgjEt47wRrKHMs8bNUXbw65VRJ5MAaEZZ9QM6lfpnBrXuvWazW9Jje7Bs3D+2HNjwwhhBBCyCk2q9iM2ssZGmrb0lQuxR32P6AVhuupAJZx07tvQ7h3xOuXXmPtz3XsrM6eFx0/JKq0hhbj8Shhul53s2Umw44alSfS9Y0YjNGuXkTLv6VORu4bmZSyRvEr3W8RRCSuIelav3MrmXZ48f2GYhOp6kOMxXOwhDzsiyU9EoQQQgghg9isYtPCmanJIz3SmunmLVWkfx+rR51Olta/JqrUePfxes/WsROPzels5LBCoL1e7Ng9UasHneZl6CRvoygCXncw/e4is4qV6kSkO9dlNlUI+9dbeesYMtFkJiuemmedm+7P7auJZp6qvpbut7cbRCpw99AcHcXTt9L11JdInmrOKVAjwKePYWThzH1xaB4bNmwIIYSQjcOGzQrZ25cW9UB4TXxv21Nz+tu53ojltcn1kq2J+SIKQarYXDibJUvR8KJYSqpIbjHMNIolwlw9azf4oyZKzdsH+ApOy4KW6f0ihhCr3nqLqxplkKsjEaK2jIgg4T120XyU9rWmPxljZCaSRmQivZqwxJI8HPDYeK/OlmDEQ2pULIElPUaEEEIImYBHHx33ZygicrOI3C8inxeR20TkCcZ5F0XkyyLygIi8OZL2ZhSblJqC31sEVUm5aT23I+1J10Q6efsLS0J484HUqBbRuT1yfhZrjpqcFankk8mlb/lwPOUpTSPHybGjgC+qwwrZmNpj05Er3FI4ScubcHeNoJ+mnMqC9X20YllGch/Zwqrji1JjWsg94KmkUVqTon+sZSXT9D5ePr0Hu5IaVabGHnfg3AXgRlU9EpF3ALgRwA39E0TkPIA/A/ASAJcA3CMiH1bVL3kJr/1RI4QQQohDZx5ekmKjqneqauc2vxvA0zOnXQvgAVV9UFW/DeB9AF5ZSnuVis3k45VDlJWWc3JTVJau9XpCkWvSbkXNLMWBz2x1lry5Paxzc4sVlu4b8USkefFmHI7u984JqQpjeGy876lFCYwQrSO5euXNxJ1Qs35iybOVu6ZERAmsYd8RVW0heoH00oi5mjQjkmlpu9XI1Hlsdv9erddin4gQTqp5PYD3Z/Y/DcBXe9uXAPxYKbFVNmwIIYQQEmeCBthVInJvb/sWVb2lf4KIfAzAkzPX3qSqH9qdcxOAIwC3Zs7LjTkWIyDYsCGEEEI2zETz2Dyiqi/w76vXecdF5HUAXg7gxarZkM1LAK7ubT8dwEOljM3SsBGRt+F4nOxRAA8D+HVVLWYWyCvZs1CS9HMusQrJ3Y03BPJO25K07xVeTZxj8rszd0pmLKcUpRmZUM0y/fb3pV5GTwEvmXsjefL2W8NtrlpufR8V34Nbn6yQfWvbI5e3Cxfy50SGxQrm9NP7TxfukKEcy//ax/LCrsYIXBqOyb1HalaFTUkLNfe/qksvMhacpjvEve8MdZVezUPNw97Q1qEiIhdxbBZ+oap+yzjtHgDXiMgzAXwNwGsA/Gop7bkez5tV9Xmq+nwAfwPgd2fKByGEELJ5lmYeBvBOAFcCuEtEPisifw4AIvJUEbkDAHbm4jcC+CiA+wB8QFW/WEp4FsVGVf+nt3kFAmNmq6HlG89NWuadA+S7mWkvqkYZ8lxzJYVg9/t8rwdWCvP21JfS735Hr8tmer/IIpiRkOGoehS5pqP/mSWt+ul3lVMzSnUs971b5nCva2rJFN12/4uIGstbur49ahSa0vQB3jHLC3tQRKTMUnh3ZEZF7z5juPdHDvNuiKUgGVT1Wcb+hwC8rLd9B4A7atKezWMjIm8H8FoA/w3gRXPlgxBCCNk6h9QAk7xfZ4SEA27o3Xk3ArhcVd9ipHM9gOt3m88F8IWx87ohrgLwyNyZWDgsIx+WTxmWkQ/Lp8yzVfXKfd1MRD6C4+9lTB5R1YsjpzkKkzVswhkQ+QEAt6vqcwPn3ltyYR8yLJ8yLCMflk8ZlpEPy6cMy2haZjEPi8g1vc1XALh/jnwQQgghZFvM5bH5IxF5No7Dvf8NwBtmygchhBBCNsRcUVG/1HjpLeVTDhqWTxmWkQ/LpwzLyIflU4ZlNCGze2wIIYQQQsZiLfNnEkIIIYQUWV3DRkTeJiKf381UeKeIPHXuPC0JEblZRO7fldFtIvKEufO0NETk1SLyRRF5VEQYmbBDRC6KyJdF5AERefPc+VkaIvIeEXlYRDjlRAYRuVpE/k5E7ts9X2+aO09LQkQuF5FPicjnduXz+3PnaausbihKRL6nm7lYRH4TwA+pKs3HO0TkpQA+rqpHIvIOAFDVG2bO1qIQkefg2Lj+LgC/rar3Fi7ZPCJyHsBXALwExwvP3QPgV1T1S7NmbEGIyE8D+CaAv4xMT3FoiMhTADxFVT8jIlcC+DSAX2AdOkaOF9O7QlW/KSKPAfBJAG9S1btnztrmWJ1is+nlGEZAVe/cra8BAHfjeDVU0kNV71PVL8+dj4VxLYAHVPVBVf02gPfheKFaskNVPwHgv+bOx1JR1X9X1c/s/v4Gjtf2edq8uVoOesw3d5uP2f3w/9cErK5hAxwvxyAiXwXwa+ACmh6vB/C3c2eCrIKnAfhqb/sS+E+JNCIiPwjgRwD807w5WRYicl5EPgvgYQB3qSrLZwIW2bARkY+JyBcyP68EAFW9SVWvBnArjlf+PChK5bM75yYARzguo4MjUkbkFLnVCtmbJNWIyOMBfBDAbyUK+8Gjqt9V1efjWEm/VkQ4pDkBsy2C6aGq1wVPfS+A2wFk15naKqXyEZHXAXg5gBfr2kxUI1FRh8gxlwBc3dt+OoCHZsoLWSk778gHAdyqqn89d36Wiqp+XUT+HsBFcP3D0VmkYuPB5Rh8ROQigBsAvEJVvzV3fshquAfANSLyTBG5DMBrAHx45jyRFbEzx74bwH2q+sdz52dpiMiTuihVEXkcgOvA/1+TsMaoqA8COLUcg6p+bd5cLQcReQDAYwH8527X3YwaO42IvArAnwJ4EoCvA/isqv7cvLmaHxF5GYA/AXAewHtU9e0zZ2lRiMhfAfgZHK+S/B8A3qKq7541UwtCRH4KwD8A+Gccv58B4HdU9Y75crUcROR5AP4Cx8/XOQAfUNW3zpurbbK6hg0hhBBCiMXqhqIIIYQQQizYsCGEEELIZmDDhhBCCCGbgQ0bQgghhGwGNmwIIYQQshnYsCGEEELIZmDDhhBCCCGbgQ0bQghE5EdF5PMicrmIXCEiX+Q6NoSQNcIJ+gghAAAR+QMAlwN4HIBLqvqHM2eJEEKqYcOGEAIA2K0RdQ+A/wXwE6r63ZmzRAgh1XAoihDS8b0AHg/gShwrN4QQsjqo2BBCAAAi8mEA7wPwTABPUdU3zpwlQgip5sLcGSCEzI+IvBbAkaq+V0TOA/hHEflZVf343HkjhJAaqNgQQgghZDPQY0MIIYSQzcCGDSGEEEI2Axs2hBBCCNkMbNgQQgghZDOwYUMIIYSQzcCGDSGEEEI2Axs2hBBCCNkMbNgQQgghZDP8H/mE+fRsq6wTAAAAAElFTkSuQmCC\n", "text/plain": [ "<Figure size 720x576 with 2 Axes>" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "import bempp.api\n", "from bempp.api.operators.boundary import helmholtz, sparse\n", "from bempp.api.operators.potential import helmholtz as helmholtz_potential\n", "from bempp.api.linalg import gmres\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "k = 3.\n", "\n", "grid = bempp.api.shapes.regular_sphere(3)\n", "\n", "space = bempp.api.function_space(grid, \"DP\", 0)\n", "\n", "identity = sparse.identity(space, space, space)\n", "double_layer = helmholtz.double_layer(space, space, space, k)\n", "\n", "@bempp.api.complex_callable\n", "def p_inc_callable(x, n, domain_index, result):\n", " result[0] = np.exp(1j * k * x[0])\n", "\n", "p_inc = bempp.api.GridFunction(space, fun=p_inc_callable)\n", "\n", "p_total, info = gmres(double_layer - 0.5 * identity, -p_inc, tol=1E-5)\n", "\n", "Nx = 200\n", "Ny = 200\n", "xmin, xmax, ymin, ymax = [-3, 3, -3, 3]\n", "plot_grid = np.mgrid[xmin:xmax:Nx * 1j, ymin:ymax:Ny * 1j]\n", "points = np.vstack((plot_grid[0].ravel(),\n", " plot_grid[1].ravel(),\n", " np.zeros(plot_grid[0].size)))\n", "\n", "p_inc_evaluated = np.real(np.exp(1j * k * points[0, :]))\n", "p_inc_evaluated = p_inc_evaluated.reshape((Nx, Ny))\n", "double_pot = helmholtz_potential.double_layer(space, points, k)\n", "p_s = np.real(double_pot.evaluate(p_total))\n", "p_s = p_s.reshape((Nx, Ny))\n", "\n", "vmax = max(np.abs((p_inc_evaluated + p_s).flat))\n", "\n", "fig = plt.figure(figsize=(10, 8))\n", "plt.imshow(np.real((p_inc_evaluated + p_s).T), extent=[-3, 3, -3, 3],\n", " cmap=plt.get_cmap(\"bwr\"), vmin=-vmax, vmax=vmax)\n", "plt.xlabel('x')\n", "plt.ylabel('y')\n", "plt.colorbar()\n", "plt.title(\"Total wave in the plane z=0\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Obtaining the solution for this wavenumber: the Burton–Miller formulation\n", "\n", "The Burton–Miller formulation can be used to obtain solutions to acoustic problems while avoiding spurious resonances.\n", "\n", "### Representation formula\n", "\n", "$$\n", "p_\\text{s} = \\mathcal{D}p_\\text{total},\n", "$$\n", "\n", "where $\\mathcal{D}$ is the double layer potential operator.\n", "\n", "### Boundary integral equation\n", "\n", "$$\n", "\\left(\\mathsf{D}-\\tfrac{1}{2}\\mathsf{I}+\\frac{1}{\\mathsf{i}k}\\mathsf{H}\\right)p_\\text{total}\n", "=\n", "-p_\\text{inc} + \\frac{1}{\\mathsf{i}k}\\frac{\\partial p_\\text{inc}}{\\partial \\mathbf{n}},\n", "$$\n", "where $\\mathsf{D}$ is the double layer boundary operator; $\\mathsf{H}$ is the hypersingular boundary operator; $\\mathsf{D}$ is the double layer boundary operator; and $\\mathsf{I}$ is the identity operator." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solving with Bempp\n", "\n", "Your task is to adapt and combine the example code in [the first tutorial](../tutorials/1_sphere_scatterer_null_field.ipynb) to solve the problem at the wavenumber you found above using the Burton–Miller formulation.\n", "\n", "We can create the hypersingular operator in Bempp by calling `helmholtz.hypersingular`. Complex number can be used in Python by writing (for example) `2 + 1j`, `3j`, or `1j * 3`. In order for the hypersingular operator to be defined, we must use a P1 space. The code needed to create the relevant operators is given below. Your task is to use these to implement the Burton–Miller formulation.\n", "\n", "Does the solution you obtain here look more reasonable that the solution above? You might like to adapt the previous example to use a P1 space to be sure that the resonances are still a problem with this alternative space.\n", "\n", "Hint: the normal derivative ($\\frac{\\partial p_\\text{inc}}{\\partial\\mathbf{n}}$) in this case is $\\mathrm{i}kn_0\\mathrm{e}^{\\mathrm{i}kx_0}$, where $\\mathbf{n}=(n_0,n_1,n_2)$. If you're not sure how to implement this, have a look at [tutorial 2](../tutorials/2_sphere_scatterer_direct.ipynb)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import bempp.api\n", "from bempp.api.operators.boundary import helmholtz, sparse\n", "from bempp.api.operators.potential import helmholtz as helmholtz_potential\n", "from bempp.api.linalg import gmres\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "k = 1 # Enter your value here\n", "\n", "grid = bempp.api.shapes.regular_sphere(3)\n", "\n", "space = bempp.api.function_space(grid, \"P\", 1)\n", "\n", "identity = sparse.identity(space, space, space)\n", "double_layer = helmholtz.double_layer(space, space, space, k)\n", "hypersingular = helmholtz.hypersingular(space, space, space, k)\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What next?\n", "\n", "After attempting this exercises, you should read [tutorial 2](../tutorials/2_sphere_scatterer_direct.ipynb)." ] } ], "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.8.10" } }, "nbformat": 4, "nbformat_minor": 4 }