Skip to content

Instantly share code, notes, and snippets.

@alexalemi
Created August 5, 2014 23:44
Show Gist options
  • Save alexalemi/b030fb320079d45bad78 to your computer and use it in GitHub Desktop.
Save alexalemi/b030fb320079d45bad78 to your computer and use it in GitHub Desktop.
illustration of the modified method of particular solutions
Display the source blob
Display the rendered blob
Raw
{
"metadata": {
"name": "",
"signature": "sha256:fa3430cc9cf68c4a2083d30e9135f07addd24446eb5600004d33ce7b062a9fc6"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Demonstration of the Modified Method of Particular Solutions\n",
"\n",
"This notebook was prepared as part of an answer to this [physics stackexchange question](http://physics.stackexchange.com/questions/128382/boundary-element-method-or-boundary-integral-method-computational-aspects).\n",
"\n",
"It it demonstrating an example application of the method described in the paper:\n",
" *Reviving the Method of Particular Solutions* by Timo Betcke and Lloyd N. Trefethen from 2005. [[doi]](http://dx.doi.org/10.1137/S0036144503437336) [[pdf]](http://eprints.ma.man.ac.uk/589/01/covered/MIMS_ep2006_365.pdf)\n",
" \n",
"The idea behind the method is to choose some parametric basis of solutions to some part of the broader Helmholtz equation, and then to find the eigenvalues for which there exist a linear combination of this basis that satisfies the boundary conditions.\n",
"\n",
"Schematically, if we have some parametric basis: $ \\phi_k(x;\\lambda)$ we see a solution of the form $ \\sum_k c_k \\phi_k(x;\\lambda) $\n",
"and we do this by choose a set of points on the boundary, for which we will ensure that our solution vanishes, and a set of points on the interior for which the solution doesn't vanish.\n",
"\n",
"The details are worked out in the paper, but this means you choose a set of points $x_i$ on both the boundary and interior, then we define $\\sigma(\\lambda)$ in the following way: \n",
"\n",
" 1. Form the matrix\n",
"$$ A_{ik}(\\lambda) = \\phi_k(x_i; \\lambda) $$\n",
" 2. takes its QR decomposition\n",
"$$ Q_{il}R_{lk} = A_{ik}(\\lambda) $$\n",
" 3. then take the smallest singular value of $Q$ for those points on the boundary."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# just some standard imports\n",
"import numpy as np\n",
"import scipy as sp\n",
"import matplotlib.pyplot as plt\n",
"from scipy.special import jv as jn\n",
"import matplotlib as mpl\n",
"mpl.cm.register_cmap(name='cubehelix2', data=mpl._cm.cubehelix(gamma=1.0, s=0.3, r=-0.5, h=1.5))\n",
"%matplotlib inline"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In particular, in this example we will be interested in finding the eigenvalues and eigenfunctions for an L shaped domain, let's choose our points, both on the interior and boundary randomly."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"N = 36\n",
"NP = 2*N\n",
"k = (np.arange(N)+1.0)\n",
"alpha = 2/3.0\n",
"t1 = 1.5 * np.pi * (np.arange(NP)+0.5)/NP\n",
"r1 = 1./np.maximum( np.abs(np.sin(t1)), np.abs(np.cos(t1)))\n",
"t2 = 1.5 * np.pi*sp.rand(NP)\n",
"r2 = sp.rand(NP)/np.maximum( np.abs(np.sin(t2)), np.abs(np.cos(t2)) )\n",
"t = np.r_[t1,t2]\n",
"r = np.r_[r1,r2]"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Just so we can see what's going on, these are the points we've sampled:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"xs = r*np.cos(t)\n",
"ys = r*np.sin(t)\n",
"plt.scatter(xs,ys);\n",
"plt.axis('equal');\n",
"plt.axis('off');"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XFX5x/HPTDJJZrkzSdrSfWeHFrAsBcsii8giiOxQ\ndhUBEQEF2UFFUX/KIpvKoojADxHZf7KUIgjIJsgOLZSl2FK6JJkl+5zfH88NcxPS0iXJTZrv+/XK\nC+beO/c8mUnP3Dn3nOcBERERERERERERERERERERERERERERERERERERERERERERERERERERERER\nERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERER\nERERERERERERERERERERERERERERERERERERERERERERERHp/6ZA5mVI1EPmn8DYzx5ScQp4H0Nq\nMSR+BkQDO+OQvgmSyyD9IfD14BPB+y0kl0L6IyibGdi3LWTegkQdZB4Ehvnbt4DMq348/wBG2eby\nY8FbAMklkLoCKAfGQ+Yp/9iXgI2BJGRutfOm3wf2BirBu86PYz5EDgYSkL7ZP+4DYF+gFjL3+9vm\nADsA47q0sQmwKWT+42970o7hC5B5zd82GxgJbASZF/1tTwMTgBrI3BdoY8cuv99SSP0GiFnc6ev9\n13Y+RA7yX6NyOya51J5Tfmzgdd0KMq/7bc4Chnf/Xlc/AYzp/D5HDrJ2ksvAu97a72RF77WI9CM1\nkFgK1xVhvoMLWiE1FygrHVJ2KIzJw4sO3nQwNQfxH5b2p2+CPQvwnoPHHKTzwDa2z/sN7JSHdx08\n6aA2D+wMjIF4Fv7i4EMHJ7VA+llgqHV6f/DjOasVvDeAvWBYHp5xMNfB9Dwkfwmp9+AnbXbsNUVI\nLIbU3+BrjfC+g4cdpAoQvwN2KcA8B084qM5D4iHYu9HiftSBV4Dki/CNZovpbw7iOUh8CD/y27i2\nCIklEF8Kv/VjvKgNUu9DvA5ucrbtjFbw3oT4J3C1f9xP/eO8J0pt3OUgkQOO6fz7bZuH1C/A+z3s\n6sf9uB83M2zftnk79hlnz2UPYATEG+AWP47TWsB7GYh89r0+v+t7PcPO/7iz9nYtQPp3nf9c0n8q\nvdezu7zXItKf7ApfqAPn7KfY0QFPLB1Scxfc4ErHPOKg9qXS/uQy+8fesf+sdoheaPvS/4VXAvt+\n5qDyMuBQ2KuhtL3dQXkr8HXYrks8XiMkboHLAud51kH1B7BOtrTNOZhSB7FGWBjYdkorVGTtg6lj\n20VFqPA72I5tP2iDSBFaAtu+nodUU5c2cjAp2znGoY0wraHztmQzrFff+bkjsxBt79zGATmI/QMu\nDWx7zkH1u+AtgrcD2y8oQuwS2/dcYPulDlK/B/aDL9V3jiPeDAy193pal9e2Jo99m8DOe0GxtP9t\nZ+0HJeo6v9c/DLzX0p9EP/8QWcvVw4JyaPEfLgPy5UBD6ZDWJfBOsfR4HuCWlR6X5WxbhzktUKyz\n/4/Wd943txXaltr534tAx2k/AnDAEphfBq3+9sVAcxm0LoJ32jrHQB1kY+A3RROwsBzKGj8bT7Tx\ns3FEmj97XKQdPvQfF7FvHM1ReByYCzQCH0ehrgya/ePqgFwZLIqW4v4EaInC0nJ7DkA9UB+DaGup\nDee3UVwKc7v8fpF6iOa6ef2WQaShy/Y2aFlir+sHUeg41UKgNQLkbV/wvV4KFALvdVudnT8YQzRL\nJyt6r0UkbAlgc2y8NQLp+2BaDs4rwqQcpC7tcvxEiC+Db7bA91ptqIKtAvv3s6/nZ7XDAY2Q/ACo\n9vftDqm8XRUf2gSJBdjYejmkn4Yd8nBuEUbmIH4mEIX0w7B1zraPz0HqZ8AoSHwCRzbD6W2QyAM7\nQuoqWNePffMcpO8EDrHhhbPbYd+CP/RwgA3PnNkOBzdB8iOIHA3pgsW9XyOk5kHFmTDcj2mXPCRf\ngWQWxjsY4mBkK6TuAu9u2MJvd3LO7gGkH4HpftzjcpD4KaTvgM3849bLgXcNVJ4C6wTa8F4AxkFi\nkf1+32+DZB7YCdgbvLxdIR/kx02t7Uvm7dgjm+21YSRQBunHYDv//GNykLjQfy8ikL5/Be91rZ3/\noCZrz8tb+518vfRe798Iyff997oG2AIYsop/i9JLImEHIH1uKiQehWExWFwBXA75c4CZEJkM7kXg\nLuySMmgMcDhEysHdAbzVZf82EN3Dv4r7A6XLaYAtILoPFHPAH7HLcYAK4BiIjAL3NPB3f3s5cCRE\nxoN7AbjH3z7c314F7i7gFexveH+ITAU3B/gzdsn9RYh8GdxS4AYgC0yD6N5QzPoxLgW2hejuftw3\nYpfXu0FkBriFkDkMTt8WziuDAvDFZnjpRP/3OBwi64L7D3AnNnbdNe4ocBhE1gf3CnCH/9ruCpHt\nrQ1uxL52DAeOgEgc3N3Ay/7vvSVE9+oSt/9eRvYF1wj8CfjY3x4DjoLIOHDPAvcF3osyv41Jy3mv\na4GjIepB8X7geT6ry3sd/RJU3Awj22BBDFqPhfbbunmeiPQe7127WekcfOLsipkvhR1V/+Utgjld\n7xl0/WYzmNVAVR6e91+f/ziIF+g8Q0dCoDH3wSUKuQlwmP+NbSiwZxSb2ifdKnsLbm+3/28E/pKH\n5ldDDal/mWBX7NP8h1OB8S3ApPBCEhmUUh/Abf5V1jJnY7LsFnZU/dgESM6HifU2syT9VzpNEx30\nhkBVAV72/6Ze77hyHxF2YCKDzTRILIMN6myKYerysAMaAKqwS9MN0H2qbsQOsQ59ozr7b8VRYUck\nMlilgenA5LADkbXGcGA7dMUuIiLSezR2KLJyxgC7AutgK5C6ThUVEQndPhC/CsrOATJhBzMA7GgL\nt3ath/FZWwikmWYBZcC3IHENcCK2TkFE+lbV6TA6D790cEgTpN4BUiEFMxL4LnAqMD6kGFZCaj7c\n788GaXawaRY4MOyo+okIeHfAVnn4lbNkZun70I1nkT4VgYqCZRLsWJTzpSxwRAixTLQMhTMb4dgm\ny2TYX+fbl7dCNvCandQMnB52VP3EZMgUoOC/Nk0OhuSx1Msi0kciUNalozq8ABzf96F4N8F57aU4\nfl2E6vs+/3lhyDxvaX2LzjIiDssD24cdVT8xBUZn7bXpeC8n1QNbhh2YyCCTvgv2bYRXHfzZ+UnA\nQlhNWPt3yzne0SE84GDIs8BwK6Ax5FlI/gqI931snzEWUm9DvAliLVB1ai+3VwOpa2DIY5Z87DMF\nM/qTCki9C+e0wmsOLmqF5IfY2gAR6UMJSN9oVZGqX8Tmu4eg4luWJXGOfzW8WR6qzrMsg6e2wP85\n2LsA6QfDie8zIli+ht7uaCsh9SYc1wx3Oti9AN6D9O8x7NFQ/aDl7q+ehVWlEpFBKgLxC2ysvSrr\nX6XvDpsHCk00dxSaGPZ5JxugIlB2JNQ+aGUB2QjYETZsKA1zNDtINmE3n0VEBqTdYEp9qWMrOKhq\nZq3NEV55quWrv8XBz9r9m8qHwkaBzr3FQaqRT+vIioh0rwrS10D1e1DzDKV0fv1BHFJz4JvNVlt1\n57zdI1hbeQvg34H7Die3QeRHVlzkxGa4z8E+jeA9SvjDMhG/uMkfIXY+nafPjoDqe6H6fah+ABgd\nTogig1r6dtijAC84q4saz/JpDc1+oRZSV8OQWTZsQ0XYAfUeb6HlP+/o3E9r8+uRDgXvRhjyr/5z\nUzl5BWycgyudVdtKvYrdf4jZB/LpLVZA/exWK1quG6oifcmfClkf6FAOzRPKVEiBqnNh/Rzc6+A3\nRX/m0gZhR9WNuM31XxIorD21AdgT2BRGBaZCFh1MaEBTIUOnZcKDTlmbFW1O+48/cZQqOK+q9bEs\ngK8BS3oiusGl6WJ4bxkcPdNK1jWezWfLF/YHMcu20PE3E8G/DVIFNEHBLxZehRUIz0dZ/b8pEVk9\nVWfDuDxc4eDoZr/Asbfq5/H+x4pLT6nzh3Z26OlIpT9JPwaHNdlw3mXttrqYoVjR7Xvgi3kbstkp\nD+mHCP8egcigdCB4N0Dsx1hB5FU1A4bnSl/T/+4gvgT9g16beZD+E1TPg+p/ABsG9pVD9LuQ/gNE\nT8UKdIvIAHQMHJgtjdsXHZS1A4mwAxMRGaxikPw51L4K1Q+zesm6toLaPMz3O/fbHCT/29OByoBR\nC5mbofY1SN/G2rvoTKQ/8/4IO+bhMQeXFyFez2rNS06cAVVNMKoBEkuAL/R0pDIglIP3KnyrGR53\ncHIzpN5irZ7CKtL/RGxK29LAVMgD88A3V/N8Q7HUrv1gHraEZBMY0WUq5PgG9GEfOlWTGXQi7ZAP\nPM45oG01T7YYeJ21e9pbBTb3XEMN3WuFlmjpT6gdaIqy+n9TIrJ6Ej+1hTM3Oji5BRILgJqwo+qn\nNoDkAhietRw3yYvDDqgfikB6Fny5AH9wsFcBvH+iC0eRPheB8mOh5k5IXAkMDzug/iv9uq0cdQ4+\ndjAyhxXJls4qofI8qLkbqi5CqQdEQjMM+DKwediB9G9lbZAP3J84oRk4Leyo+qmJwO7A5LADERms\nvgiJBtiqDobmbeGJFh91z3sP/tfv2LMOJueAfcKOqv+p/CYkC7BNHaQKUHVS2BGJDELJ/8I9foeV\nczAxC+wRdlQhSWLz/Je3Sncrmyq6eZ0VfdYHYTeGQVUjvOX/Tb3rINGI0v6K9KkoRIpWAKJjqOG4\nAvCdsAMLwQ7WcY9usM6p8lvLOa4W2AmY0nehDShbwMRA9SznYMM6YNuwAxMZZLy34Cr/JuFHDobl\ngBlhR9XHYhCvs5w4zlkdVy+PZbmUVZOBqhw84b+W/3IQz2NrIESkD21oQzNDc1DZDImzwg4oBGMg\nU+h8tblTHbAvtuLyWqgsWG3XqnPRUMzn2d06+GFZ+2/0q2EHJDJY1QD7AZuFHUhIKqEyB0+70jeY\nTB7YxOayb5uHBf4V/cQ8lB0RdsADwDjgAGzWjIiE4AsQXwbr14NXAO+ysAMKyV42fDDFn+GRPAuI\nQM2Hlnen44r+Bgc1fw072P6t7GCIF2CDOvtvTB+GIn0v9T7c6ndcSx2MGcwLc4ZjRUb8q83k5bBO\nG1wd6NxPbYXkVSs4x15QfglwMoNz8U4tVBVKtWBfc9bBMyLswEQGk4E8WyYKiZ9DZd7GwxO/prTE\nfSpU3w1DHofYt1m9MfKEJVV7zMEwB99wsJ+DRB0wqvunxM+CUTm4qAi7FcB7gXCzIUaAvYATgW36\nqE3NlhHpH7x5cJM/W2axs86JncOO6vNVngKb5OADBx86mJqHqu8D61mZv0uLcKe/0Kjqh6vRQA1U\ntkCbg3kOLncwuRE4cjnHl0F5i8XjHLQ72KyB8BY5RSB9K0zOwtEFy7df9d0+aNe/cn/Bfx1e7rhy\nV1oLkT62GSQW29VWshFSP1+Nc+wD8aug7FyguqcD7N6Q2XBH4OrwHgdDnoToBfC9ttL2lxx4H69e\nG+l/wpFNdo4rgnVCuxO39AStgZi+2gAcvpq/4JqabrlvCn4s8xzEmumT6ljl+9v9i0n11tHHDu39\nNkWkO3tA9PfAuazyP/7K02BUHn7prGBy6l0g1QsxdpG5FS4MdOI/a4fMXyF6PpzSWtr+ogNv4eo2\nYu1kPoDqfwIbddk/1uqIDpkNVafbh8FxzdaR3ub8QuFj1+z3XG37wI5dhkfSBZY7pNSjIsBxwO+B\n41FGSJEwxI63nDLnFuErBauis9I3AiM25v12oAPZJQcc1YsBd5hks3wOLsDhjbYIifWAydap/rII\ntzuYmIPKH/Rgu5XASGAYJBbBWW32DWKLHCSug8w9kFwCmdfpu3Hu7oyGeA4ednZP5dJ2u3neFx2t\ndyNslIMLHEzJgXcbWhsg0tcqczajoaNqznZZ4LCVfHIEylqhPtC5H1EATujFgINGAif5P8Er0k0g\ncwcMmQXlx9FjHUvscKhosimjFfWwZ2Dh0yJnN2BX2HmOgfjldrXfJ2PxO0Nykd00T78BrNsHbU4A\nrxEa/Ncl76AmjxU4EZE+EoVolzS2R65i55y+E/ZphFcd3OLsanGtTPO6gc1/f7VjVlER9iyWXrcl\nwc69ElL/A7UvQfW9wCRgBMSXwGmtcJWDEXkoX17+mp7Wl8MiU2FsQ+fhoPXqga37MAYRgfSDMLPJ\nZnnc5yCxqjlVqiF1M6Q/guqXWHunvB0Me9R3vlKPO7ioDe51sE0eUr+3Q9N/gV0LMMvBxW12w5oL\n4ejm0vOfW5N7Af1ZJSQ/gl+0w3wHl7dD8mP65EauiARlIP03y4jovcPKT4OMQOJiu1qtaIH0Myw/\nVe7aYDoMz0Gd3zm/4CDW6A//PAXxc4EyoOKzRT2+nAVuh9MCN4DnOEguXcUYBsq49STIPGN1AjIv\noARsIqGYAN4cqwkaa4bKlR2S2d9uVn7sbC74N5shc1dvBho+72rLmrlzPSTzEP16NweV232IpYHO\nffsG4DTLa36zgycdbOkgNR+IfX67ZYfah2+0HTKP0/8zLM6A+FKIN0FVPfClsAMSGYTSL9vQQdHB\n3I6bXysxPlp5Gfws0IG97cBb1Ovhhm8rLMnaChJipa6GzXPwJwcnNkPyQ8Cz4uPTHGzlbCbJ1jng\n6M9pb0ubwviCg0YHJ7VA5rEe+l16g2cfRA/4fxePdEwJVdF1kT4UhUg7NAc66WMbsdknn+cUq3Df\n7j/vhiJUv9jbAQ8QUag4GWrvheSVWI1aLE3CJ4HX+vttwNmfc67TrFZrx3ManK2E7be6ST+wQR0w\nPezARAaZ5Cd24885aHKwQRbLY/55qsB7HjZqgF3r7WpNBbZXLPMPOK3FPhDfczAsz+cPWcyEbXKl\nD9F/OnvP+q0REG8spWH4r7OVz4wLOzCRwWY3myGzez2Mz0L6LlZ+6lwM+AqWt1u5Qz7fcEj/225C\nl7dA5fdW4jkV4D0DW2bhmAKk8qzch2+IEqdb8ZOvNNgwX/KcsCMSGYyS4N0B3jLIvMPaO5WxP0lh\nM2tWVgVwCDZcNgUrqnIm8G3A6/Ho1twYSP0LqpZB6nlgfNgBiQxC6bvha422SvXWjkVIk8KOqh+a\nhA2hjAwxhghwMaRa4eQ22Cvv5/JJhxhTVxWQmgfntsEbDn7cBsn5QDzswEQGk6jNyc4Gbn4dXsCu\nCAeK4dhYfy9ewSbOtNWpW9RZtsPo13qvrRVJ/gqqi/Bg4P3arxE4NZx4ujUFxjTY7KuOGCfVY7OM\nRKSP+Im/5gT+Ie6UBQZIWbT49+zm3fiOG7ozeqGRDW0q4keBlaWVefr+SjRji8VGOcs62fF+nVuE\n6E/6OJYVmWzj7R2phpscDMkDG4cdmMgg05Gy9+cODmqC1FwgGXZUK2FTK2L9nt+J/J+zLJE9nkdl\nL/hiXeepfbU5YEIPt/N5RtuHzLEODnSw0MEzDqqbgO37OJYViVj6hS39NNDT85C+h4GzulZkrVEB\nibugehmk32Xg3FA90Gb4BDvdRBOfzinvMZMgWbDxY+fgIeevuqzs4XY+TxS8N6yG6wEOPAfJIkSP\n7+M4VsY6UPkcpJZB5b8J9z6FyGDl3Qw7F2xJ/FVFiDcQXoGJVbGZTbP7r9/pzu7odFdlFspKqjga\nKhthVIM//LNjz7exUkZCZhYkl0HmJWzmTH9TDt7rcFIzPO3g1Bb/22CYtWRFBp2IzbleFrj6PSgP\nrGoq2gpC+dqdPMcWyGxU58/y2TWwcwuo/ivUPgiRA3ugsWqsEpOyG67YJjAiW7qhWnQwoQGYFnZg\nIoNJBCoCqwmdg6+sTL6TDuMh/SpEi1DVENIskvHAdsCQwLZNrLP/ddESdY3IQ9mxK3Gu9bH55DPo\n/GE1BSpmQ+0SqHm7mw+LCJaPZ1eUQ2U9+0bVkdKixcE6OfrntwyRtVniRzA5B79zcEKL5eJe2SLX\n3puWdKzdv7mXLLD89K614P0ehjwNqUvp1Zu2VZfBeYFCGrMdVM8JHBCBqrMhscyGoZK/gshBFv8e\nDVZYOn29HccUS+1b4yDtIOkg0wZ0fJCVQfoBGJWDaXV+Ee3B3JFFwHsQds7DdQ52L0D6MVRHVaTP\nxSH5mM08ySwEvrySz0vakE5wPvN+WeDIbo6thNQblhb4Pn9utvcEvTaUE78CLgx07o87qJ5b2l92\ntKUrfsOfVrh5HiparZi2czbvf2QemAGp62Csg038fDAfO5jhwOs43/Ewrbl0pfr7ImT+0zu/14Ax\nDhJvQnUB4nPQojiRMKT/CnsX4CUHN3WkZ11BOttPRW1I5xW/U2t2MCkL7NbNsdvB5MDCllZnc6F7\nbTrh5pYv57cO/uZgbB5igUyXtffa79rR+f/VQVWx88yb3R1wLWRuhjEO/hjYN8tBbT1QBRWL4PzA\nvg8dJOp76fcaCGJ2A/XMVnjZwfmtfmHulS26LiI9wC9w3RDonA4rACs5va5sJqTzcGge1staKoNu\nr8a3hUmBzr2ltzt3gOlQ8yAMeRLKj+kcV+p6ODtQFemKIiRb4Vr/8X8c1DqoLELlQ5Bog9MCr9Fl\nzioNcSBMboRNndVQLTo4x/kFNQarKTA62/kb3cR6YMuwAxMZTCJQUYB3Av8Qd87R/dDK8kwFvgHs\nyfLHVSvAew2OarIr6a8WID2b8Ba2TLBKQUc02VBRPAtlD0DCH1dPO7jNwT4ODnaQaoJEO+zr4NCi\nPzNnE+BY2D8HZ/rzzkf7Y/J9vsCpP/FXqDYGvtEN7Xi9RKTvVJ0BY/JwqYOZzZCcR+/kaam2CkVD\n/gGJXxB+IqmRWE6W72PDUF+DcmdDCa3ObhJPd3C3/5N5DftGcxKldQATrKP/i7NKSfs1g/doKL9N\n/xGxtNHT8/YNZ0Ye0n9HK1RF+lwteC9CRRukskBPzAkfoKoegJEOLvTH3Lf3rzyfdFAzZzlP+iJk\n3rXVol47xArLqa26ssqAw4BzsFz53dkTamfbD3utQVu9ZWOo+siGs6oWYCmKRaRvZWbb0MRiB486\nK/w8qJM87QnRW6xg+N8c/MvBxjmbOtmtMkgsglv8YYgXOl7DCavRdgTS98JmOfhBG4zO2VRVMsDX\nsemXX4fqPPzZ2U913mLuN+KQ/BiuLFqR8N8VIbEYy2EvIn0kCtH20vioc3BUATgx7MD6ga9BzWtQ\n3WDTJMuaoPxuyLxqhcAzt2N51P2EXsGZNjvWsXrVkraFMVnLpOicJQeLtUBiIcxogG0boLqx88yd\nmxzUPtyzv/oa2QzGNXR+PdarA7YJO7DBTgsNBpciVOThbf+hA95sB5aGGFNfGAaZxyDWDKmFdH/l\nexfUPwhtHpxSDjMrIbYP/HQTeGEY7LkPZO4AlkBTBF7zn7YMeDUGfLAacdXA2PZSTrJ1gHgUTh4C\nT3jwlAejKqEYeEo72BvXXyyDJTGo8x9mgUUVrP1/UyL9TdmRtoDp+FaYkQPvBdb6JE+ZJ+HkFqhz\ntno1mQc2/OxxNQ1wu3+TdU8H4xxs6+ybTpOzQifEIHa41TbdrR6G5sG7dDUDG2aJyW5zNkx2URvU\nNMH9gavg8x2k223153XOpqJ2u7YgRKmrYFIOTm+F9XLg3RB2RCKD0XhIfQCVbRBrhdh3wg6ol5Xb\nUFRLcG7/cpKlVWfhXgfDHVztj79/xcFMZytbY02UZoGsB+zPms/n3hLSc6CyYHPp41faUv68s5Wz\n2+Wh/AaouQ9q7gd2WcP2ekME2Ac4C9gPzZQRCUP6ZfhJmy06meMs6RNbhx1VL/KrT3WsrG13sHkW\n65i7SFwDExwcEvggyDqIOUtPUPWDPoi3wgqYl7dCeRt4fwbK+6BdERnAohBpL+VFcQ6ObcTmcq+K\ncmyRyroMiKu08qNtlsnJzbBNDryngVh3B0JsNuwYeH3mOygvAnv3bcwkCH9tgIgMHMlPbApkR73L\nDbOs2kyPdSw75IisJYry7qf7jrK/mQ6cjtWLXVG8HiTfh2NaLD3B5ALE+1PNUhGRbn0Z4nnYLgtj\ncn5+mFWYNZW5E05psWGdJgfb56H89F6LNhxDoPISyPwZIjMZEN9ORGSwmwKJJTCmAPEWSF68ak+v\nfgeeCwxbXOsgc0vvhNrjyrF0tEPDDkREpId578ANfrrbRQ5G5ICdlnPwV8G7ESp+wadFjzP3wBl+\nXvcWB7vmofyMvol9jYyD1LtQm7PVqKlL0RW5iKwlohApdp4WeFwB6GY6ZOx4GJ6Hyx2c1GpL7hkG\njLROcmK9lVPzHmVAzJNPvGZpfaMOvuAseRphlAkUEekNqQ8tq6FztqhnbA6rxjQE6+SvBA6H1CeW\nN+XTueFNWFZFsCWVW2Ll5QbCKufxlt73UT8D5C+dpevlwrADExHpKVtBVR1MyEKqEbwrgfFQuRRS\nDvZ2sFnR8pnPC3Tup7Rii1QGogNgl/bS71J01tlzQtiBiYj0lGmQqIMNcpBuAu8KSP8FJju4OdD5\nrVuEbVvgWWcZEBN5Bm4Bhh1sGKYjYdpcB7F2VApORNYeqQ8sl4lzsMzZdMiaN2yY4u3AlfqPHSSe\ns9kxNc8D24cd+RqIQPoOy3tyZJNlWiz/dthBiYj0lOXcUK14AMa2wbH+mPR8B2NaWL00tv1VBFtl\negKwVcixiIj0NO9d+IM/FfITByNzwK6QvM2qC5X5y+2rzgs7UhERWXlTIb4YRmYh3gTJSwL7YkAS\nK/0mIgOYss0NPm2WPMxFwUWAxsC+Vv9HREQGFm+O1bl0DhY4GJ4Ddgg7KhERWX2rsEK116wDHAjs\nRam+nIiIrJnUfLjT79gbHIzPAbsv5+BpUDMLhrwI8bNY89Wom0BiGezSAFMbwHsZG+MXEZE1NN3q\ndm5WZ1WYvN/RfQKt9SGeg2uL8LCDqXlIXNLNcaug+im4ulhaKLVvI5T9cM3OKSIiALWQmQWVTeB9\nDOzR/WGRs+C7raXhm7ccJJf4O8cB29m5VkXmQys+3XHOSx2kfreK8UeADYBpaIWpiEiHzKPwjWZL\n9zvLQTIPbNzNgT+A45pLHfFLDlKLIHEmJBphozq7smfXVWj7VpjZZAulPnGwfg44fBWCL7P6ojV5\ny0qZnA8Hws0uAAALj0lEQVRMWIXni4islaIQbSvlWHEOjioAJ3Zz7GiIL4Nz2uFGB2PzEPsFZAq2\ngtU5eMxBVQMrP6U2A+nHIdYKZa2Q/BWrllP9OPhCDgp++xe3QfUTq/B8EZG1VVUW/hNIELZNFjh0\nOQdPhNR1UHM3lB0G7A+71Zc+GJyzzJKMWMUgPLqfKfMVqH4L0gvA++1nj6m8FH4WaHuus9TEIiKD\nXuwoG9Y4rRV2zIH3Iis/JXFjSBfgPb9zfdBZ+uAeWdG6BaTycJ+D1x3sUoD09V2OOQa2yEG+I7lZ\nG2Qe74G2RUTWCidBZBZwM1C9ak+NnwLxRphQD/EGemwBVORc+H5b6ar8PWepiTuJgne7DQ2Na7AM\nl4zvmfZFRAa0im/BkDycU4TdC+C9wqrPOhkObIENr/SUU+HgplLn/pQDb0E3x0WAdYHN0SIoEZEO\nlTl4LTDmvl0WOCzsqIBaSHwERzbDj4s2dBTpD3GJiPR7/myZQuCm5JEF+k+5uSEQORsq/gfYKexg\nREQGkMxDcEQTfOjg/o7yeRuEHZWIiKyZDGTuspuV6XeBXcIOSER6nvK5Dz5NUHwPKheBWwAsDjsg\nERFZY+lbbA75Uw6uKfrTGceGHZWIiKy+CJS3Ql3ghuoheeBbYQcmIj1rTfNzy4ATbYP6wOOlDmgJ\nKxoREekR8QtgUh6udXB8i59ZMRN2VCIismYiEDkcMrdA1S+BYWEHJCI9T7NlBp8yiK8L5etDJAVN\nNYAyK4qIDGyp38H0PPzdwSXtlrOdkWFHJSIiqy8C5S2wODBbZv888M1ebncolj1yci+3IyI+zZYZ\ndCKu8+SYZge092KDu0B8Hky5xzJQpn7Ui22JiAxWyV/BJnm41cGZrZBYBAzppcaiUFUPs/1vCR87\nSzfMtF5qT0Rk0IpAxXeg9u/g3QCMCewbBuwP7AFU9EBbtVDV3Lks314NLL+sn4iIrIEaYHtgw8C2\nje3m6pcaYEoDeP8BkmvYTgTiS+DuQHWlTB6YuobnFRGRLraGeD1sWmfl6ryrgQhUPwVXFUtFPPZt\nhLIf9kB70yGxDMY0QFUTxL/XA+dcHRGoOgOq34WaNyFyUEhxiIj0htR8uMMf/57tYHQO+DJkPoSX\nA8MnlzqbNtkjEsCm2KyZjmLae0H6DxC/lM5DQ72k8nRYPwdP+4W9a/PA7r3frohI74tCpAjXO6hx\nsLmDlAN+C5lbrYhHq7OpkhvkgJk92PYukFhs7ccXwLBG+I2DU1tt6Ka359rXvgaPBj68rnSQubl3\n2xQR6TPJeVDt4A2/k3vWQVUjMAbST0BFi2WOTP4aK0bdE0ZDIgePOGhzMNrBvwId7TEtwJk91NZy\n1D4L/xto8/x2SF7bu22KiPSdo2HTYucZLKMbgE38/Wmgqofb3Ad2qCu1N87BW4H2z2iHsgsCx38V\nah+A6r8B2/RQDLtBqgA/dXBmu5/Hfv0eOreISOjGQKJQ6lyfd1CVxzr13rI1jMhB3m/zBAfTnBUM\nubWjjqs/gyZyoM2F/4OzYZtEHtiyh+KYDvErofJXwLo9dE4Rkf4i9m2IN8F6DRDPQ/RrvdxgBNI3\nw8QszMxDTR4Sf4eaOVDzHLBj6dDaF+GewFX9Lx14f+zl+ETWOsoKOfhsBrGfQG0LfFgB0d9B8a5e\nbtNBwxHQsDvMGw+8CDwLhe6OjXT+sywHIkqTISKyYt48+KM/5v6Jg5E5YOewoyopOwKG52265o0O\nUnlgux46eTmkroDkEvAWQPmxPXReEZFQ+VMhWwLDHscVgO+EHVhnkUNgyBNQ8wiwU8+dN/lzS3c8\n18EzDoblgT177vwiIqFJfQC3+R37UgdjcsCuYUfVN2resamfHR9slzlIXRd2VCK9QWOZg05uPzim\nDjaoh/GNUHc98EjIQcUprVztTQ0wL/DwnTZoXdIH7YqI9Ik0Nn887OIZtZD+J5S1WRGR+Fm93N5O\nNrXy9DY4shkSn6AqVCIiPS1zL3yj2VIefOBgVA7Yq5cbnQqR84DvA8N7uS0RkT61JXAaljumJ/K2\nr6bkEpgXGAO/qAjll4QXj4jIgFV2mOVUP6EZtslB+hkgFk4smTdKN3fbHezWD2fuiIgMCPF6+Heg\nQ52WBQ4OKZjpEM/CVxtgswbwXqDn89qIiKz1ohBth8bAUMjRBeDEEGMaiw0P7UOoQ0QiIgNa5gk4\nqQXqHTzRsQJ007CjEhGRNTMUMrMh1gzJRcC+YQckIiIiIiIiMjj1xZJvkd70JeyG7HrAa0B7uOGI\niISjCrxrofo9qHmGnqtyFIKq71pmxzPbYUYO0k+jGgUiMjil/wJ7FuAFP196PAtMCDuq1RCF8mZ4\nJzBnf9Ms8NWwAxMR6WsRKGu1aZAd89wPywPHhx3YaqiyhGNtgd/la1ngiLADE+kPlPJ30Clrg6WB\nx4sc0BhWNGugCZIvwPda4WPgXuChCPBEyHGJiISh6mwYl4fLHRzVDMn3AQ+oBMaw4jwztcAwINLN\nviiWPtfrZl8EGAFkutkXB0bz2Zv7VX48XcfQq7FsjhGLJfMwVOXA+4DORUeiwCgg1U0c6W6OC8Zd\n5seU6CZeEZF+6wDwrofyHwE1EN0PKguQKUB8KbBtl+PLwftfqGyBRBOkZwPJwP6R4L0OqUaItUDy\np4F9Q8F7EZJNUNECqd/w6YdDxbFQ0QTpAiQW8OlK2bKDS/EkFgNbAVHwrrNzJJr8hGfVwAyIL7Nj\nK/MQ3RsYC6k54PnxJM4H1gHvZYuxogWSvwZGg/dWKe7Ej4ANIfmRxVTRDJUn9OxLLyLSN8ZYAYvn\n/XHre5x1llSWDqk8A7bLQ85Z/dUDGsG7prQ/MxvOaIWig48djMthuWKAzD1wUrPd8FzqYMMccDiw\nKaTz8Jbf7g1FSM4HJkIyDy/52+9wEF8CkW/B1DzUORtnP6YJ0rdbIrT/84/9l4N4HtL/hgvbLJ7/\n+kXA08/BKS22bbGDdXOQfgPO8eNe6JccTC2Ea/wC4nMdVOeBL/TVmyEi0lN2h63rSjclnYNhOTpV\naaq9F24K7J/toPbl0v5EPcwP7D+3CJGLbJ+3EN4I7PuFg/hvgJmwb0Npe9FBrBU4GHboEk+mAN5f\n4MrAthccePNhVEPnYzevsyvuxYFtp7daDp25gW0XO6hqgQWBbT9oh2jRYunYdnAOOLbv3g6RnqEb\nqvI+vFkBi/yHbwANZZQ2AE1vwcPN4PzHs9qgfW5pf2w+zPJ3tgKPFMC9b4/L3ivtawceaoTGudbu\ns1Fo8M/xLBBttgBeiUFHadNXgMYIFF6BB5ug6G9/tAiRd2FJDOb42xYA71RA5UJ41N/WAsxqguii\nUhxtwEMFiC0tHdcKPNpow0Qd92TzwNPOYhURGXBSF0FNHnaqsyvc2JFdDsjYmPrGDbBlvY1JMyaw\nfwtI1MGMOhifhfQjlG6EbgyJJbBtHUxuAO8pbMgnYouphuVgxzobGmJvP55LoDZv25N5KDsESFi+\n9/UbYJt6SCwC1oWKb0CqYLFnCpA8B9gG4g2wfR2MzUL6PmCqDTdtVwcTs+A9juWTr7e4x2Uh/SCw\nByT8mIbnwLuR7m8gi/Rr+qOVDlOASdgS/rnd7K8Etsc67SeBbJf9w7Ci2/X+/mJgXy12kzYP/BO7\ndO6wJTZb5SXgg8D2zYDxwKvAu/62mB9DFfAUUOdvXw/YGHjHPx5sRs3WwDL/2CIwBJgO5Pw42oF1\n/OOCcY8FtgAWAs9R+soiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIi\nIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIi\nIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIi\nIiIiIjIQ/D9rSKytSyA4SwAAAABJRU5ErkJggg==\n",
"text": [
"<matplotlib.figure.Figure at 0x1095a23d0>"
]
}
],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we will define $\\sigma(\\lambda)$ as we described above.\n",
"\n",
"For our basis of representations, we will use the functions\n",
"\n",
"$$ \\phi_k(x; \\lambda) = j_{\\alpha k} ( \\sqrt{\\lambda} r ) \\sin(\\alpha \\theta k ) $$\n",
"\n",
"since these are the solution to the helmholtz equation on a semiinfinite wedge. If we use these functions at the interior angle of our $L$ shape, this will guarentee we satisfy the boundary conditions along the interior walls of the $L$, which will help us out quite a bit with the accuracy. Here, since our interior angle is $3\\pi/2$, we take $\\alpha = 2/3$"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def mat(lam):\n",
" return np.sin(alpha*t[:,None]*k[None,:])*jn(alpha*k[None,:], np.sqrt(lam)*r[:,None])\n",
" \n",
"def sigma(lam):\n",
" A = mat(lam)\n",
" Q,R = sp.linalg.qr(A, mode='economic')\n",
" s = sp.linalg.svdvals(Q[:NP]).min()\n",
" return s"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 4
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To see how this function behaves, let's look at it as a function of $\\lambda$."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"lamvec = np.r_[0.2:40:0.2]\n",
"plt.plot(lamvec, map(sigma, lamvec));\n",
"plt.xlabel(r\"$\\lambda$\")\n",
"plt.ylabel(r\"$\\sigma(\\lambda)$\");"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAERCAYAAABl3+CQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYVOXZx/HvFpCOEgxNVAQ0AokUYZHmihpBREBRgmiM\nhEiMmFhiiHkTszGWWGJMosGuiRrRJIpgAVFZQUAFBKRLEUNRVLCAiLAw7x/3ThjG2d0pZ85zzszv\nc117uTtz5sztYXfu87T7AREREREREREREREREREREREREREREcnQAGAlsBoYn+D5psBUYBGwFPiB\nb5GJiIhzRcAa4EigFpYMjo07pgy4qfL7psBWoNif8EREJFahg/fsgSWK9cAeYCIwJO6Y94FGld83\nwhJFhU/xiYhIDBd36a2ADTE/bwRK4o65D3gF2Aw0BM71JzQREYnnokURSeKYX2FdUi2BzsBdWMIQ\nERGfuWhRbAJax/zcGmtVxOoF3FD5/VrgXeAYYH7sQW3bto2sXbs2S2GKiOSstUC7ZA920aKYD7TH\nBrNrAyOAyXHHrAROqfy+GZYk1sWfaO3atUQikcB//fa3v3Ueg+JUjIpTcUa/gLapfGi7aFFUAOOA\nadgMqAeAFcDYyufvAW4EHgIWY8nsF8A23yMVERFnU05fqPyKdU/M9x8Dg/0LR0REquKi6ynvlJaW\nug4hKYrTO2GIERSn18ISZ6oKXAeQoUhlf5uIiCSpoKAAUvj8V4tCRESqpUQhIiLVUqIQEZFqhT5R\nDB8O110Hq1a5jkREJDeFPlFccQV8/DGUlkLnznD77bBNKy5ERDyTM7Oe9u6FWbPggQdgyhQYNgzG\njYNu3RxHKCISMHk766moyFoVjzwCq1fDscfC0KFwxhmwYIHr6EREwitnWhSJ7NoF998PN90E3bvD\nrbdC+/Y+RiciEkB526JIpE4d635aswZ694YTToDrr4fdu11HJiISHjndooj33ntw6aWwbh089BCU\nxG+XJCKSB9SiqMYRR9hAd1kZnHkm3HUXqAKIiEj18qpFEWvtWjj7bOjYEe69F+rX9zgyEZGAUosi\nSW3bwty5ULs29OwJGzbU/BoRkXyUt4kCoG5dePBBuPBCG+xetsx1RCIiweNq46LAKCiAn/8cmjeH\nU06BF1+Eb3/bdVQiIsGR94ki6vzzoVYtOPVUeOkl6NTJdUQiIsGgRBFjxAioqICBA+G112yWlIhI\nvnM1RjEAWAmsBsYneP7nwMLKryVABXCwH4GNGgVXXw2nnQaffOLHO4qIBJuL6bFFwCrgFGATMA8Y\nCayo4vgzgMsrj4+Xta1Qr7oKFi2CqVOtS0pEJFeEYXpsD2ANsB7YA0wEhlRz/HnA49kP60C33AL1\n6lkZcxGRfOYiUbQCYlctbKx8LJF6wGnAf7IdVLyiInj0UZg2Df75T7/fXUQkOFwkilT6igYDrwGf\nZimWajVuDP/+N/zsZ7Ciqo4xEZEc52LW0yagdczPrbFWRSLfo4Zup7Kysv99X1paSmlpaWbRxTnu\nOKs4O2oUvP66reQWEQmT8vJyysvL0369i8HsYmww+2RgM/AmiQezGwPrgMOAL6s4V9YGsw98E9sE\n6VvfgptvzvrbiYhkVRgGsyuAccA0YDnwBJYkxlZ+RQ2tPKaqJOGbggLbAOkf/7BWhYhIPsnb6rHp\nePJJK1G+cCEcdJBvbysi4qkwtChC65xzrPvphhtcRyIi4h+1KFK0aZMNcL/5Jhx1lK9vLSLiCbUo\nsqxVK1uE9/Ofu45ERMQfalGkYdcu6NDBdsY7JVFhERGRAFOLwgd16sDtt9tCvD17XEcjIpJdShRp\nGjIEWraECRNcRyIikl3qesrA8uVw4omwejUc7EsRdBGRzKnryUcdOsCgQfDnP7uOREQke9SiyNDa\ntVBSYq2KQw5xGoqISFLUovBZ27Zw5pnwpz+5jkREJDvUovDAunXQowe88w40aeI6GhGR6qlF4cBR\nR8GwYfDHP7qORETEe2pReGT9eujWDVatgqZNXUcjIlI1tSgcOfJIa1X87W+uIxER8ZZaFB5atsxK\neqxfrzLkIhJcalE41LEjfOc78Hi1m7eKiISLEoXHrrjCpsoGqKEjIpIRJQqPnXYaVFTAjBmuIxER\n8YYShccKCuDyy626rIhILnCVKAYAK4HVwPgqjikFFgJLgXJfovLI+efbDnirVrmOREQkcy5mPRUB\nq4BTgE3APGAksCLmmIOB2cBpwEagKfBxgnMFatZTrN/8Bj75BO6803UkIiIHCsOspx7AGmA9sAeY\nCAyJO+Y84D9YkoDESSLQLr7YZj/t3Ok6EhGRzLhIFK2ADTE/b6x8LFZ7oAkwA5gPXOBPaN5p3Rp6\n9oT//Md1JCIimSl28J7J9BXVAroCJwP1gLnA69iYxgHKysr+931paSmlpaVexOiJMWPgjjvggtCl\nORHJJeXl5ZSXl6f9ehdjFD2BMmxAG+AaYB9wc8wx44G6lccB3A9MBf4dd67AjlGA7afdujXMnAlH\nH+06GhERE4YxivlY19KRQG1gBDA57phngD7YwHc9oARY7l+I3qhVCy68EB54wHUkIiLpc5EoKoBx\nwDTsw/8JbMbT2MovsKmzU4G3gTeA+whhogAYPRr+/ndrXYiIhJGKAvqgXz8r7TFsmOtIRETC0fWU\nd8aMgfvvdx2FiEh61KLwwc6d0KoVrFgBzZu7jkZE8p1aFAFUrx4MHgxPPuk6EhGR1ClR+GTkSO1T\nISLhpK4nn+zZAy1bWrHANm1cRyMi+UxdTwFVqxYMHw4TJ7qOREQkNUoUPlL3k4iEkRKFj/r0sdLj\ny5a5jkREJHlKFD4qLIQRI9SqEJFwUaLw2ciRNk4RkjF4ERElCr917Woti/nzXUciIpIcJQqfFRTA\nOefAU0+5jkREJDlKFA4MHQqTJrmOQkQkOUoUDnTrBp9/DqtWuY5ERKRmShQOFBbCkCHwzDOuIxER\nqZkShSPqfhKRsFCtJ0d274ZmzVR6XET8p1pPIVG7NgwYAFOmuI5ERKR6ShQOqftJRMLAVaIYAKwE\nVgPjEzxfCnwGLKz8+rVvkflo4ECYNQu2b3cdiYhI1VwkiiLgTixZdABGAscmOO5VoEvl1/W+Reej\nRo2gVy+YNs11JCIiVXORKHoAa4D1wB5gIjAkwXFhH2hPypAh6n4SkWBzkShaARtift5Y+VisCNAL\nWAw8j7U8ctLpp8OLL8K+fa4jERFJrNjBeyYzn/UtoDWwExgITAKOTnRgWVnZ/74vLS2ltLQ04wD9\ndMQR0LQpvPUWHH+862hEJBeVl5dTXl6e9utddO/0BMqwMQqAa4B9wM3VvOZdoBuwLe7x0K6jiHXl\nldCkCfw6J4fsRSRowrCOYj7QHjgSqA2MACbHHdOM/f8TPSq/j08SOWPAAJg61XUUIiKJueh6qgDG\nAdOwGVAPACuAsZXP3wMMBy6pPHYn8D3/w/RPv37w9tu2Teohh7iORkTkQGGfWZQTXU9gg9qjR8Pw\n4a4jEZFcF4auJ0lA3U8iElRqUQTEO+9A//6wYYPtgiciki1qUYRU+/ZWKHDZMteRiIgcSIkiIAoK\n1P0kIsGkRBEgShQiEkRh7w3PmTEKsCqyLVvCBx9A/fquoxGRXKUxihBr2BA6d4Y5c1xHIiKynxJF\nwJx0ErzyiusoRET2U6IImP79lShEJFg0RhEwu3bBoYfCxo3QuLHraEQkF2mMIuTq1IEePWyLVBGR\nIFCiCCB1P4lIkChRBFD//jBjhusoRESMEkUAHX88rF0LW7e6jkTC7Isv4NFHYexY+30SSZcSRQDV\nqgV9+kAGOxdKntu9G/r2hUcesT1OevaEZ55xHZWElYuNiyQJ0e6ns892HYmE0XXXQatWMHmy1RE7\n6ywYNAjatYOOHV1HJ2GTyvTYusBI4NtYgqmH7XW9HXgD+Fflz37KuemxUQsWwAUXwPLlriORsFmw\nwJLCokXQvPn+xx96CG67DebNg3r13MUn7qU6PTbZA08BOgDPAfG9nQXAd4CTgZeBxcm+uQdyNlHs\n3WvrKZYtgxYtXEcjYXLGGbZj4k9+cuDjkQiMGmUtjVtvdRObBEM2EkUd4DBgTRLHdgT83FEhZxMF\nwNChMGIEjBzpOhIJi0WLrDWxdq2tyYm3ZQt06gSvvgodOvgfnwRDNhbc7aLqJFEbGBrzc7JJYgCw\nElgNjK/muO5ABXBWkufNKf36aeGdpObGG+GqqxInCYBmzeA3v4HLLrMWhkgy0pn11Agbq3gCeBVr\nRaSiCLgTSxYdKs91bBXH3QxMJfylRtLSt68ShSRv3TqbAHHxxdUf95OfWCn7F17wJy4Jv2RnPTUB\nhgNDgPrADuA3wMI03rMH1kJZX/nzxMrzrog77jLg31irIi917gzr18Mnn9gUR5Hq3H03XHghNGhQ\n/XHFxdbyuOYa2yyrUJPkpQbJ/orcic1oGg2UYkljBKm3JgBaARtift5Y+Vj8MUOACZU/52UjuVYt\nKCmB2bNdRyJB9+WXNqvpkkuSO/7MM23m08SJ2Y1LckOyLYpLgU9ift4F/B/wR2AB8EgK75nMh/4d\nwC8rjy2gmq6nsrKy/31fWlpKaWlpCqEEX58+8NprNpNFpCr/+pet6G/bNrnjCwrghhsssYwYAUVF\n2Y1P3CovL6c8gxW8XvT9D8e6iJLVEyjDxigArsFaKzfHHLMuJramwE7gR8DkuHPl9KwngJdfhmuv\nVatCqterF4wfD0OGJP+aSAROOAGuvloLO/NNNqbHHgQ0BD5O4tjDgf/WcEwxsApbd7EZeBMb0I4f\no4h6CJgCPJXguZxPFDt22EyVjz+GunVdRyNBtGKFreTfsMHGH1IxaZK1LN5801oZkh+yMT32K6wV\ncB62OjuRQ4CLgSOSOF8FMA6YBizHZk+tAMZWfkmMBg2s5MK8ea4jkaB6+GFbxZ9qkgAbq9ixQ2Xt\npXqp3EO0AC4CvoktwqsF7MW6hTYC9wGfeR1gDXK+RQFw5ZXQtCn86leuI5GgqaiA1q3tg/7YRJPM\nk/Dww/DYYzB9uqehSYBlq4RHrEvYPxvJtbxIFE8/Dffeq3nv8nXPPQfXXw9z56Z/jt27bRD86adt\nQFxynx9boZ4LDANOwloVkmW9e9sHwd69riORoHnsMet2ykTt2raa++abaz5W8lM6ieJYbNrqUVjL\n4tueRiRf881vWhXQJUtcRyJB8sUX8PzzcM45mZ9rzBjb/2TduszPJbknnUQxD5gEPACMwQa6Jcui\n6ylEoiZPtg2JDj0083M1aGCruicEpVNZAiWdRHFP5VcHoAF5umrab6r7JPEef9zbysKXXGID219+\n6d05JTekO3P6cOD72LTY+6l6DUS25cVgNsDq1XDyyfDfmlapSF7Ytg3atLG1E40aeXfeQYNg+HC4\n6CLvzinB48dgNtiiuuuBq3CXJPJKu3awcyds2uQ6EgmCp56CU0/1NkkAXHop3HWXSpDLgVQ3MiQK\nCqw/+vXXXUciQeB1t1PUgAFWrfjNN70/t4SXEkWIKFEIwPvvw1tv2XanXisstLGKu+7y/twSXkoU\nIXLCCUoUAk8+aaU3slX7a/RomDIFPvooO+eX8FGiCJHu3WHhQltJK/nriSfge9/L3vmbNIFhw2x/\nCxFQogiVRo1spsvbb7uORFzZvBlWrrQZcNk0Zgw8+KAGtcUoUYSMup/y26RJNjZRu3Z23+eEE+y/\nmdSQktyhRBEyPXvqjzefPf00nHVW9t+noMDGKh58MPvvJcGnRBEyYWhRbN0KY8fCsmWuI8kt27bZ\ntNXTTvPn/b7/ffjPf2y/CslvShQhc8wx9kH84YeuI0ls8WLo3BlWrbLFW+rj9s6zz9pOdvXr+/N+\nzZtDv362H7fkNyWKkCkshJKS4LYqyspso6WXX7aFW/9OZTd1qdZTT/nT7RRL3U8C6dd6Coq8qfUU\nq6zMpsjeeKPrSA60eTN06gTvvQcNG1rZ6tGjYc0aS3CSvi++gBYt7Noecoh/77tnj+2gN3MmHH20\nf+8r2eVXrSdxKKgD2g88AOeea0kCoLTUFoUFtfUTJlOn2r+7n0kCoFYt2xhJrYr85ipRDABWAquB\n8QmeHwIsBhYCC4D+/oUWfCUlMH++7ZccFHv3wv33w8UXH/j4Oeeo+8kLLrqdoi66CB59VDss5jMX\niaIIuBNLFh2AkdiuebFeAo4DugA/AO71Mb7AO+QQaNUKli93Hcl+s2fDwQdD164HPh5NFPv2uYkr\nF+zebfulDxni5v07dLBdFl991c37i3suEkUPYA2wHtgDTMRaELG+iPm+AfCxL5GFSPfuMG+e6yj2\ne+YZK/sQr2NH2z1N1UjTN2uWjQ+0aOEuhgsugEcecff+4paLRNEK2BDz88bKx+INxfa6eAH4qQ9x\nhUqQEkUkYiuGhw5N/Ly6nzLz7LMweLDbGL73Pfs33rnTbRziRrGD90x2mtKkyq++wCPAMYkOKisr\n+9/3paWllJaWZhZdSHTvDn//u+sozLJl1n993HGJnx882BZv3Xabv3HlgkjEKrm6TrQtWkCPHhbL\niBFuY5HUlZeXU15envbrXUyP7QmUYWMUANcA+4Cbq3nNWqzLamvc43k5PRZsX+OmTW3xXZ06bmO5\n4QbYsgX+8pfEz+/bZ4u35s2DI47wN7awW7XKCgBu2GBlNVx65BErcT5lits4JHNhmB47H2gPHAnU\nBkYAk+OOacv+/4no8Gh8kshrdevaKu1Fi1xHApMnVz/QWlhoZSemTvUvplzx7LNwxhnukwTYGNSs\nWdqnIh+5SBQVwDhgGrAceAIbixhb+QVwNrAEmx77ZyCL1ffDKwjjFB9/DCtWQN++1R83YIDN3JHU\nRBNFEDRoAIMGWatC8ksA7lMykrddT2DrFmbOhH/8w10MEyfCP/9prYrqfPQRtGtn/812iexc8emn\ncPjh8MEHUK+e62jMCy/AddcFc8GnJC8MXU/ikR493Lcopk1LrprpoYdaV9mcOdmPKVdMm2ZF+YKS\nJABOPRXefRfWrnUdifhJiSLEOnSAjRvhs8/cvH8kAi++mHzZ61NOsWKBkpznngtOt1NUcTGcfbYq\nynplzx5rpf3mN/DLX8LDD8P27a6j+jolihArLrYpqQsWuHn/pUttxlW7dskd378/vPJKdmPKFXv3\n2gfI6ae7juTrzj3X9u2WzEyfbn+/v/ud3XQ1bGgLV9u1g7vuClaJfhfrKMRD0QHt/g6qYaXSmgDo\n1cv2q9i+fX/hQEnsjTegZUsbowiaPn1sOvQ776iibDoiEbjlFrjzTvjb374+q23pUlt3tHAhTJhg\nhRldU4si5FyOU8yYkVqCqlfPEttrr2UvplREIrYY8J57glfw7vnng9maACgqguHDNfspXb/4hU0A\nmTvXfv/ipz536mSTVDZvhvPPD0adNCWKkHM1Rbaiwj7wTzwxtdf17x+ccYoFC2DJEnjsseCtNp46\nFQYOdB1F1c49V4kiHbffbjcB5eVw2GFVH9eggVUM3rjRxi9cU6IIubZtrStnyxZ/33fhQtvQ5tBD\nU3tdkMYpJk2yGkYvvGAzjIKyN/SHH9pmTyec4DqSqvXqZXt4r1jhOpLwmDwZ/vQnuwlIZl+ROnXs\nd/Txxy1puKREEXIFBW5aFTNmwEknpf667t3tQ3BrANbZT5pkK8rr17e4glJGe/p02/QpCH3TVSks\ntGKPalUk5913YcwYmy3WunXyrzv0UOumuuQS64pyRYkiBxx/vG1k5KcZM+zDLFW1a0Pv3u4/lFev\ntlXlJSX282mnWasiCKZNs5XsQafup+RUVFjL9Ze/tF0KU9WzJ/zkJ7aBlKuZUEoUOaBbN3jrLf/e\nb88e26go1fGJqCB0Pz3zjLUmont5f/e7NovLtX37Up9N5kpJiXV7Ll3qOpJgu/lmaNwYrrgi/XP8\n3/9Z9/Ljj3sXVyqUKHJA167+rqVYsADatIFvfCO91598svsB7eefP3Ax23HHwSefwHvvuYsJbPpw\n48Z2fYOusFCtiposWQJ33GHldjIp7FhcbLPzrrrKfk/9pkSRA444AnbtsppAfigvT298Iuq44+zu\nyFWf644dNqYT+/9QWGjlKVy3KpItiRIU0UQRpMVhQbF3r41L3HCDN+thSkpsc7CYLXh8o0SRAwoK\nrFXhV/dTuuMTUUVF9voZM7yKKDXl5TZ43aDBgY8HoUts6tRwJYru3eGrr+Dtt11HEjz33WctgTFj\nvDvnddfZdG6/a20pUeSIbt386X7as8cWCvXrl9l5XH4oV3XXftJJlrxc3R1v327/hmHapLGgQN1P\niXz4IVx7ra2sLvTwU/bQQ+Hyy23Mwk9KFDnCrxbFvHm2dqNJk8zOc+KJtvrUharu2tu0sbnrK1f6\nHxNYkiopsem6YRKt/aTup/1++1s47zz4zne8P/cVV9gGUn5OiVeiyBF+tSgyHZ+I6tjR1lL4PU6x\nbp3duVf1BxxtVbgwdWo4psXG69rVkkQQdlsMguXLbY/za6/Nzvnr17dENH68f8lZiSJHHHWUfQBm\ne5vKTMcnogoLbVe8WbMyP1cqXnrJBq2r6g5wmSimT7dpumFTUABnnQVPP+06kmD4xS/gmmsyb3VX\nZ/RoeP99/7YXVqLIEQUF0KVLdrufKiqsqmmfPt6cr18//7ufXn7ZpudWpbTUFgP63Y2yfj18/rkV\nhAujs85yX2YiCF5+2cqaXHppdt+nuBhuuskW8flRNFCJIodku/vp7bet/IBXd0r9+vnboti3zwbQ\nq0sUhx9uJdCXLfMvLtifwLwc+PRTSYnVfnrnHdeRuLN3r61zuPlmOOig7L/fkCE2g3DKlOy/l8tf\nywHASmA1MD7B86OAxcDbwGwgC8NCuSXbA9qzZ1sxOK906WJ30tu2eXfO6ixZYsXYaqq146L76aWX\nbAfAsCostDn++dz99OijNn5w9tn+vF9BAfz613D99dlvAbtKFEXAnViy6ACMBI6NO2Yd0A9LEL8H\n7vUzwDDKdoti9myr0+SV4mKrkOrX/hSvvJLc/hl+J4pkWjphkM/dT7t32wDzLbdktgI7VUOHwpdf\nZr9OmatE0QNYA6wH9gATgSFxx8wFortBvwFUU71dwLZQ3LYte3foXicK8HecoqbxiaiTTrJxCr82\njFm6FBo1shX2YXbiiVYZeONG15H478EH4dhjvf/7qElhoa2p+P3vs9uqcJUoWgEbYn7eWPlYVX4I\nPJ/ViHJAYSF07pyd7qf//tdW4Ca7P3ay/EoUe/bYeEgyU3tbtoSmTf1bbRz2bqeoWrWsftakSa4j\n8deuXVam43e/c/P+555rsx2zWZHZ1Z7ZqeS+k4DRQMJcXRZT+KS0tJTSMC1rzYJo95PXHzzR1oTX\nzeru3W3eebb30Z43z6YQN22a3PHR7qfOnbMXU9RLL8EPf5j99/HDsGHw17/CuHGuI/HPfffZ70mP\nHm7ev6jIpuNef33VU9fLy8spLy/3MyxP9ARiZwBfQ+IB7e9gXVRV3cdG5ECPPBKJnHuu9+cdNy4S\nufVW788biUQi/fpFItOmZefcUdddF4lcdVXyxz/+eCQyeHD24on66qtIpGHDSGTr1uy/lx+++CIS\nadQoEvnoI9eR+GPnzkikRYtIZMECt3Hs3h2JHHFEJDJ3bnLHk9rNurOup/lAe+BIoDYwApgcd8zh\nwFPA+ViykCRka0A7G+MTUX50PyU7PhFVWmpdVXv3Zi0kwNalHH10dhdn+alePVvQ6MeUzSC4+26b\nGty1q9s4atWCK6+EP/4xO+d3lSgqgHHANGA58ASwAhhb+QVwLXAIMAFYCLzpf5jhc/TRVsL700+9\nO+f27TY/Plt/DH37ZjdR7NxpOwD27Zv8a5o3h29+M/ub8uTK+ESsYcPyY5rszp02y8nV2ES80aOt\nuzQblWVdrqN4ATgG61a6qfKxeyq/AMYA3wC6VH456gEMl6Ii2+9h4ULvzvnGG9YHm61FRL16WSto\n167snH/2bIs/vqx4Tfxo6eRiohg0yGqCbd/uOpLsevBB26Y0G4X/0tGgAfzoR7ZRktdCug5UquN1\n91M2u53AfsG/9a3srQFJtdspKtuJYvt229HO7ymV2Xbwwfb/5FcdIhd274Zbb7VB5CC57DLbr8Lr\nKfJKFDnI6xXa2U4UYOefPTs753711fQq3kYTRbbmp8+aZbO+6tbNzvldGjYstxffPfaYdfO6mulU\nlZYtrbTH3Xd7e14lihzkZYti717revKydEci2UoUX3xhpTtKSlJ/7RFH2P4U2apf5FXJ9iAaMgRe\neMHW3uSavXvhD3+AX/3KdSSJXXWVTVH28torUeSgb33LVsd+/nnm51q61AZ2k11/kK7evWHOHO/v\n3ufOtZpS6d61n3hi9hYyeVWyPYiaNbNKuK5KtmfTU0/ZLLWg/tt16mRjco895t05lShyUHGxDbB5\nsZHMG2/YgF22HXaYFVTz+u595szMtm3N1jjFZ5/ZTnrptHTCYsgQeOYZ11F4KxKBG2+01oSfNZ1S\ndcUV8Oc/e3fjpUSRo7p29ab76Y03/Pswy0b306uvZp4osrE/xaxZ1r/tRzlqV848EyZPzq0tUqdN\ns66nQYNcR1K9U0+1AXevyvgrUeSorl29mSIb5kSxa5cly0zGV9q3tzpR773nXVxg4xNB7brwyjHH\n2Iw2P/Zy98ttt8HVVwd/35CCAiuj8pe/eHO+gP/vSrq6dMk8UXz+Obz7rn/zxL1OFPPmQYcOmdWQ\nKijITvdTLg9kx4q2KnLBokXWXThihOtIkvP979sY0YYNNR9bEyWKHNWxo63Q/PLL9M8xf74NitWq\n5V1c1enUCT74wLt9vzPtdoryekD7009h1SqbGpvrcmmc4vbbbZ1C7dquI0lOw4Zw/vkwYULm51Ki\nyFEHHWTzvDMpQeFntxPYqvKSEpv95IWZM+1DPlNetyhmzbIJArk8PhF1wgmwaZP3XXd+27QJnn0W\nLr7YdSSpGTcO7r8/sxtGUKLIaV26ZNY/7HeiAO+6n/bsgddfhz59Mj9Xx46wdSts3pz5uSC3p8XG\nKyqygd/X5ynVAAAQKUlEQVSwFwn861/hggtsK90wad8ejj8eJk7M7DxKFDksk3GKSMRNoujTx5tE\n8dZb0KaNN3/YhYVWUNCrGST5Mj4RFfZxiu3b7a788stdR5Keyy6zRJfJ7DMlihyWycynDRtsK1C/\nt+csKbH6R5kWCMx0/UQ8r7qftm2D1avtLi9ffPe71rr77LOajw2iBx+0xN6mjetI0nPaabBjR2Zd\nukoUOey442yMoqIi9ddGWxN+LyqqX9+bAoFejU9EeTWgPWuW9duHZUDUCw0aWIssjEUCKyqsGutV\nV7mOJH2FhTZW8de/ZnAO78KRoGnYEFq1shk2qXLR7RSV6TjF3r3w2mup7T9Rk86dbd/wjz/O7Dwz\nZuRXt1NUWLufnn7aCu35UZ0gm37wA3jxRRuUT4cSRY5Ld5wizIliyRKrNdSsmXcxFRfbwr3XXsvs\nPPmw0C6RM86wIoF79riOJHmRiO0Y9/Ofu44kc40awXnnpV9VVokix6WTKPbssde4muefaYFAr7ud\nojIdp9i2Ddaty6/xiahWraBtW+8mBPhhzhxrQZ55putIvHHNNfDDH6b3WiWKHJfOFNmlS+Hww6Fx\n4+zEVJNWrTIrEOjVQrt40bpP6Xr1VWuV+LWAMWjC1v30l7/AT39qU3xzQatWcOSR6b3WZaIYAKwE\nVgPjEzz/LWAusAsI8VCSW126WOmBVO7OXXY7RaXb/RSJeD/jKap7dxvvSXeLz3ybFhtvyJDwFAnc\nuBGmT7e+fXGXKIqAO7Fk0QEYCRwbd8xW4DLgNn9Dyy3f/Kbdna9fn/xr5s1zX16iT5/0xgNWrrRB\n/NatvY/poINsyvHcuem9Pp8W2iXy7W/blOtly1xHUrO774ZRo6xvX9wlih7AGmA9sAeYCAyJO+Yj\nYH7l85KBVMcpFixw34+ebosiW91OUekmsI8/tmTdtavnIYVGQYF1PwW99tOuXXDffTalVIyrRNEK\niK1puLHyMcmCVBLFl1/a2IBfFWOr0rEjbNmSeoHAmTO9nRYbL90V2q+9lt/jE1FhGKd44gn7mznm\nGNeRBIerRBGCXsrckUqiWLzYFrzVqZPdmGoSLRD4+uupvW7OHG/qO1WlVy/rmtu9O7XXZWvcJGz6\n9bMbkfffdx1JYpGILUy77DLXkQRLsaP33QTE9iK3xloVKSsrK/vf96WlpZTmcydwFVKZ+bRgAXTr\nlt14ktWrl33wDx6c3PGbNsEXX1jV3Gxp3BjatbPrmcoirFmzrEx1vqtd20pKPPss/OhH6Z9n+nRr\nbZ53nnexgY0/ffopDBzo7XldKy8vp7y8PO3Xu9r1tRhYBZwMbAbexAa0VyQ4tgzYDvwxwXORSBim\nUDgWicA3vgErVtS8CO2ii+wDcOxYf2KrzvTpcP31yU9JffJJ21A+233gl11m04evvjq547dvhxYt\nbJzCdUstCB5/3P6dnn02vddv324bUlVU2LiPl+XaR460lmxYCwAmq8Bq8yT9+e+q66kCGAdMA5YD\nT2BJYmzlF0BzbBzjCuDXwH+BBr5HmgMKCpLvfgrCQHZUSYnFk+xq3tmzM9v2NFl9+6Y2oD13rg1i\nK0mYgQOtK27nzvRe//vfQ//+NovqiSe8i2vzZtsT+6KLvDtnrnC5juIF4BigHXBT5WP3VH4BfIB1\nSTUGDgEOB3b4HGPOSCZR7NwJa9bYTnNB0KiRreZdtCi54+fMsdlS2Rad+bRvX3LHa3ziQAcfbN2b\nL7+c+ms3bYIHHoBbbrG7/jvu8G5dxt13W4vC1ULTINPK7DyRTKJYtMia9EHaeS06TlGTL76A5cv9\naQ21bGkfditXJnf8rFnZnYkVRmeckd5mRlOmWIukWTMYMMD+3b3Yv+Srr+DeezUltipKFHkimUQR\npG6nqF69klvg9uabNqXXr+6dPn2Smyb71Vd2Xf3oEguTwYNtjCLV1sDkyftrLxUWwo9/nH6hu1j/\n+pd1ZR0bv+xXACWKvHHMMdYHW93mMfPnB2fGU1SyLQq/up2ikl1PMW+eTTdu2DD7MYXJ0UfbNUml\nDtmOHXbNTztt/2MXXmgJZ+vWzOKZMAEuvTSzc+QyJYo8UVRkd9yLF1d9TBBbFEcdZStlN2yo/rjZ\ns/1PFMkMaGt8ompnnJHazKcXX7RNn2LHEJo0sdbJP/6Rfhxvvw3vvWfxSGJKFHmkuu6nHTusBHbH\njv7GVJOCgpq7n/bts+f97N45+mgb/K8pgSlRVG3w4NTGKSZPTrymZuxYuOee9Ae1J0ywNR3FrlaV\nhYASRR6pLlEsWmSznYK4RWdN3U/Ll9s6ES83KqpJQUHNdZ8qKiyBZXOleJj17m03J5s313zs3r3w\n3HOJE0Xv3jZekc5eIdu3w8SJMGZM6q/NJ0oUeaS6RBHEbqeomhKF3+MTUTUNaC9eDIcdBk2b+hdT\nmNSqZeMNzz1X87Gvv26zzRLtp1BQsL9VkapHH7U1Ga1Uaa5aShR5pFMnWL3a+vzjBXEgO6pbNytN\nXdUCLb/HJ6JqGqfQtNiaJTtOETvbKZHvfx+efz61Pc0jEet2uuSS5F+Tr5Qo8kidOlanaOnSrz8X\n5BZF3bo2dXH+/MTP+7UiO16XLlZC4pNPEj+v8YmaDRxo+3R8+WX1x9WUKA45xDZG+vvfk3/vOXPs\npql//+Rfk6+UKPJM165f737avt1mfXTo4CamZFTV/bRli02NdBF7cTH06JF4wVckohZFMpo0sYQ7\nY0bVx7zzjk3rrqnFO3asLZpLdlB7wgRbh1GoT8Ea6RLlmUTjFAsX2h17kPdKqCpRzJljUyZd/bFX\n1f20ciU0aJCdnfZyTU2rtKdMsUHsmv6NTzjBJmMkUyT1o4+sy0tbnSZHiSLPJEoUQe52ioomivi7\nRVfjE1FVLbxTayJ5Na3SrmpabLyCArj44uQGtR96CIYOtRaN1EyJIs907gxLlth0w6j584OfKFq2\ntDv01asPfHzOHLflMUpKbHZTfB+7xieSd8wxVl8s0WLQrVvtxubkk5M71wUXwNSp8OGHVR+zb58l\nEw1iJ0+JIs80agTNm1u/b1SQZzzFiu9+2rXLVtX26OEupvr1bZHivHkHPj5rlhJFsgoKql5898IL\nNthct25y5zr4YBg2DB5+uOpjpk+3vwOXvzdho0SRh2J3vPv8cyvdHIZiaPGJYv58i7t+fXcxwdfX\nU7z3niWx9u3dxRQ20e6neDXNdkpk7Fi4776qy8BHp8QWuNq2LYSUKPJQ7Mynt96C444LR/mC+ETh\nalpsvPgB7WhrQh9EyevTx1q5H3yw/7GvvrL6ToMGpXaukhKoVy/xTKoNG6xb0OstVHOdEkUeih3Q\nDku3E1hRw/fesz2Nwd2K7Hi9e1upjui4z8yZGshOVe3a8N3v2qK5qFdftWnPqZZmqW6l9n33WZJo\noL0yU6JEkYeiiSISsT7gk05yHVFyiouhe3cr5xCJuB/Ijjr0UNsTe8kS+1kD2emJnyabTrdT1KhR\nNhaxZcv+x/bsgfvv1yB2OpQo8lCzZrZK+623rEURW98/6KLdT++8Y90Lhx3mOiLTty+88orNtvng\nA1uXIqk5/XS7hrt22Y1AdP1EOho3hrPOsmmwUc88Y5UJglYhOQxcJYoBwEpgNTC+imP+Uvn8YqCL\nT3HljS5doKzMmvv16rmOJnnRRBGUbqeoMWPgtttsQLZ3b9v/Q1LzjW9Y92J5uc1mKy7ObMV9/KC2\n6jqlz0WiKALuxJJFB2AkED/n5nSgHdAeuBiY4GeAXitPZqmoz7p0sQ+14cP3PxbEOOP17Alz5pQz\nc2awEkWPHtZ18rOfWesiDNcSghfnqFF2DW+7zbqdohMC0omze3ebBvvyy7BqldU4O+ssb+ONF7Tr\n6RUXiaIHsAZYD+wBJgJD4o45E4iW93oDOBjwcbcBbwXxl6drV1vkdPrp+x8LYpzxmjSB+vXLefLJ\nYCUKgBtvtPn+/fuH41pC8OL88Y/h2mttfCL2Qz2dOAsK4MorYcQIO9fo0fY7n01Bu55ecZEoWgGx\n+4JtrHyspmMC0hudG/r3h7/9LZx7ObdubXV/OnVyHcmBmja1WVlayJWZUaOsFpMXM8cuuMAmGVx+\nOVxxRebny1cuZs8nu2Fh/Cz0NDc6lEQOPtjusMLo8MNtplEQ134ku4JYquflToutWtlWp5I+F0uC\negJl2BgFwDXAPuDmmGPuBsqxbimwge8TgZjJboB1YbXNUpwiIrlqLTYOHFjFWJBHArWBRSQezI4u\nvekJvO5XcCIiEgwDgVVYi+CaysfGVn5F3Vn5/GKgq6/RiYiIiIhIbktmwV4QrAfeBhYCb7oN5QAP\nYuM9S2IeawJMB94BXsSmJLuWKM4ybBbcwsqvAV9/me9aAzOAZcBS4KeVjwftmlYVZxnBuaZ1sCnx\ni4DlwE2VjwftWlYVZxnBuZaxirB4okVSgnY9PVeEdUkdCdQi8RhHULyL/YMETV9stXvsB/AtwC8q\nvx8P/MHvoBJIFOdvgSvdhFOl5kDnyu8bYN2qxxK8a1pVnEG7ptFaAcXY+GQfgnctIXGcQbuWUVcC\njwGTK39O6XqGsdZTMgv2giSIxaZnAZ/EPRa7yPHvwFBfI0osUZwQvGv6AXbDArADWIGtBQraNa0q\nTgjWNd1Z+d/a2I3hJwTvWkLiOCFY1xJsDdrpwP3sjy2l6xnGRJHMgr2giAAvAfOBoM/kbsb+6cdb\nCPZK+MuwSQ4PELwm85FYK+gNgn1Nj8TijM4oDNI1LcQS2hb2d5UF8VomihOCdS0B/gRcjS1DiErp\neoYxUYRp4V1v7I9xIHAp1pUSBhGCe50nAG2wLpT3gT+6DecADYD/AD8Dtsc9F6Rr2gD4NxbnDoJ3\nTfdVxnIY0A+IL4QflGsZH2cpwbuWZwAfYuMTVbV0aryeYUwUm7BBuajWWKsiiN6v/O9HwNNYt1lQ\nbcH6sAFaYL9cQfQh+3+x7yc417QWliQeASZVPhbEaxqN81H2xxnUa/oZ8BzQjWBey6honMcTvGvZ\nC+tmehd4HOiP/Y6mdD3DmCjmY1Vlj8T6Bkewf4AmSOoB0UpK9YHvcuCgbNBMBi6s/P5C9n+IBE2L\nmO+HEYxrWoB1MywH7oh5PGjXtKo4g3RNm7K/u6YucCp2Nxy0a1lVnM1jjnF9LQF+hd1MtwG+B7wC\nXEDwrmdWJFqwFzRtsP7LRdhUxCDF+TiwGdiNjfdchM3OeolgTZeLj3M08A9syvFi7Jc7CH3VfbBu\niEUcOC0yaNc0UZwDCdY1/TbwFhbj21jfOgTvWlYVZ5CuZbwT2X9THbTrKSIiIiIiIiIiIiIiIiIi\nIiIiIiIiIiIiIiIiIiL56BDgnwRzLxKRlBW5DkAkB+3CkkVHYIHjWEREJKCaYRVFRUIvjNVjRcJg\nC1Y1uJHrQEQypUQhkh11sE2BBrkOREREgqcI26y+C/CE41hERCSA7sAGsgHmYhtsiYSWZj2JeGs4\ntvXtS5U/N8c2hVnpLCIRERERERERERERERERERERERERERERERERERERydz/A5g4znqxZqTpAAAA\nAElFTkSuQmCC\n",
"text": [
"<matplotlib.figure.Figure at 0x10c31ea10>"
]
}
],
"prompt_number": 8
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can clearly see those points where it is vanishing, those correspond to our eigenvalues. Let's use a scalar minimizer to find the first 6 values."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from scipy.optimize import minimize_scalar\n",
"lams = [ minimize_scalar(sigma,(z-0.1,z+0.1), tol=1e-16).x for z in [9.6, 15.2, 19.7, 29.5, 31.9, 41.5] ]\n",
"print lams"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[9.639723844634462, 15.197251925760987, 19.739208802179888, 29.521481117131952, 31.912635948136323, 41.474509890573152]\n"
]
}
],
"prompt_number": 10
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Great, these agree very well with the results in the paper. Now that we know what the eigenvalues are to good numerical precision, let's determine the eigenfunctions. Following the method, these are determined by\n",
"\n",
"$$ R(\\lambda) c = \\hat v $$\n",
"\n",
"where $R$ is the R from the QR decomposition of $A(\\lambda)$ as above, $c$ is the vector of coefficients defining our solution, and $\\hat v$ is the right singular vector associated with the smallest eigenvalue of the boundary part of $Q$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can in fact plot the lowest 6 eigenfunctions and give their associated eigenvalues."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"x, y = np.ogrid[-1:1:200j, -1:1:200j]\n",
"rr = np.sqrt( x**2 + y**2 )\n",
"tt = (sp.arctan2(y,x)+2*np.pi)%(2*np.pi)\n",
"\n",
"fig,axs = plt.subplots(2,3, figsize=(10,5))\n",
"\n",
"for pk,ax in enumerate(axs.flat):\n",
" Q,R = sp.linalg.qr(mat(lams[pk]), mode='economic')\n",
" u,s,vt = sp.linalg.svd( Q[:NP,:], full_matrices=False )\n",
" c = sp.linalg.solve(R, vt[-1])\n",
" ans = (np.sin(alpha*tt[:,:,None]*k[None,None,:])*jn(alpha*k[None,None,:], np.sqrt(lams[pk])*rr[:,:,None])).dot(\n",
" c)\n",
" ans[100:,:100] = np.nan\n",
" #ax.imshow( ans[:,::-1],cmap='cubehelix2', extent=[-1,1,-1,1]);\n",
" ax.contourf( ans[:,::-1], 10, cmap='cubehelix2')\n",
" ax.set_title(r\"$\\lambda_{} = {}$\".format(pk, lams[pk]))\n",
" ax.axis('off');\n",
" ax.axis('equal');\n",
"plt.tight_layout()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAtEAAAFiCAYAAAAwf4RBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuwJndd5/HPGWKSCbmREEISw5wJR0miJiYSFhHxLqRY\nvLCgK8ISXV23vG65XvCGaHlB19Wq9bK7lqWs1rpW6eKtEEu01oXaKAYzBmUCciAzxCQOmDtmkpDM\n2T/66Zw+fX73/v26f/30+1U1lZnn9G0m/Xv683yfb/9aAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAwBQ2pj4ARPl6SZdK+oSkD0j6PcMyZ0v6Hkl3SjpX0s9K2pH05aufPVfSP0n6Jcc2\nD0i6X9Kpznb/RNJXS3qNpEskvUDS70r6rdXPbdt/haRPlnSmpOOS3to73hdI+iJJP+l5PfbvBdTg\nMyW9VtJ3dV77kJox8YCk75b06xHr5hqvrmMI3a8Stu97PwCmFHPut1zjL+a6K9nH01dIunq1j7sk\n/cbqddt4sr1u245tv7bXgSK+XdJPFNr2Z0h6V+fP71AzQPp+VdKh1e/ft/r9+ZIeXS2/Iene1eu2\nbR5WM3g2V8t9h6SrJG1J+rbVss9U88ax6dj+5dr7ZvQrat5QWgck/bGkN/b+DqbXY/5eQIiS41WS\nvlPNxevXeq9/o6TnSDotct1c49V1DDH7jd2+7/0AcJlivIZcd23jL/a6axtP50n6687yf7H6uW08\nmV5/umM7pv0edryOjgNTH8Ca+XlJXyXp4gLbfpmkOzp//qikz+ktc4WaT7fHV3/+0tXvH5D0WWoG\n9I6aC9uGY5uPqflkfEzSQ2o+Ld8u6dPUVIOl5lP1tqTnO7Z/kaQvlnT6ap1/lvR4Z3+vlvSn2v+N\nSP/12L8XEKLkeJWab0t+3/D645I+IumJyHVzjVfXMcTsN3b7z5T7/QBwmWK8hlx3beMv9rprGk83\nSHqJpKOd5W+T9AWyjyfTdfcTju1cbdjv8y2vf5awh6sSgng7kn5T0usk/UzgOleoqdrY/KWagf2w\npE/qvH6mmk+7f9Z57QvVDNzXqfkU/LCkt6x+9r7Vf18s6c/VDPiQbX6TpJ9b/f6PJN24+v2Gmq95\nth3bP6bmg9otkn5ZzVdc7UXzIklPSvqYmk/Jcrwe+/cCQpQcry3Th7obJJ2hpi3p7yX9gWVb/XVz\njlfXMYTu979Hbv+I7O8HgM8U4zVkzN3d+X13/Elx113TePqgpM9Wc/1rPSDpUyT9tszj6VbL659s\n2c5/suz3fYbXtwUUdoWkvyuw3WdJ+is1J/M5q328obfMD/T2/S41g6T1Skn/S83XNCHbvED2N6t/\nqf29Yf3tS00Afoekk2reYFrftNrv6yX9sOf12L8XEKrUeG29XvvbOb6y8/u/UfPBMGTdnOPVdQyx\n+43Zvu39AAgx9ngNOfdbtvEXe92V9o6n71cTdFs/qt22Ftt4Mr3+fY7tmPYb8vriUYnO7yJJZ6lp\nxP+rjNv9qKSvU/Op+h5Jf7t6reuh1eutj6hpffjg6s9vVfOp9IikL1Hzqdi1za9W87VU3/mSblJz\nA0ZXf/unS/r81e+/WM2b09+qqTS/W01lofvJ/4WW11P+XkCIUuO1ZapEdytf96sZI6YLVH9d33tA\nzHh1HUPsfkO3f1Tm94ObDccMmIw9XkOuuy3b+Iu97vbH00NqAnrroKQTagpJn6/94+mfLK8/LOlC\nw3Zs+/W9DhGic3uZmhP7x9QMkg9K+ndqBsh7tbepvxXzddNR7X499EZJP9Rb9n2SPrfz51NqvtZ5\nuZpPs58j6eOr43mVmk/Nrm1+gfbPHLCh5lPzN6y2dUjSp1u2v6PmKyep6XF+vZqvtU6qeSN86Wqd\ng2ruYn6O5fWUvxfg0x+vH1AzI8zztH+2mFbs18M7vZ+/VtKXqentlJqWJVtvdH9dKc94/VzPMcTs\nN3T7T6qZMcD0fkCIRoj+eL1PzU1610j6QzVtDH1Dx6vkv+62+uMv5bprGk8f0u69BlIThI+oGWOm\n8fSk5fW/7W3nmdr9NzPt97jjdSC710j66dXvz5H0D2oG0AvUfFj5nwO3v6nmRgCp6Z/6nc7Pnqvm\nZD9TzZtC6+bVz14m6cdXr21ot5Lr2qbUDNQv6r327WpuLni2mr/b5zm2/0pJX9NZ90Y1Nzd0vUl7\n2zlMr8f+vQAf03g9U81FwnQ+prpJe78efrGar1ql5gPjHav/Srvj2LbupvKMV9cxxO43Zvsh7weA\niev6eo6aXukcblL4ud8fr/3xl3LdNY2ns7T3m9jb1LSEmMbT51lef4ljO7b9ul7HCjMZ5PFCNRXn\nr++89gtqBuRVauY2frt2m/RTfJKkH9Tu1zg/quZrUqn5NPlv1Qzil0l6kZpK7e3aDe/fLOlpakLC\nB9XcxODaptTcPPEtkt6/+vOLJf1f7Z43O2qqx3dZti810/08Xc0dwg9I+h+d7X+Vmk+5O5LerN1P\nz93Xf1LNm0zM3wtwsY3X90j6P2rG7Y9k2M+3qjmXL1dzI+zPqflq9mvVfC19SM28q+9eLd8dx6Z1\nTyrfeLUdQ8x+U7bvej8ATGzj9RY158/Vas63Hxi4n9gx1x2v0v7xJ8Vdd13j6XWrbRxQU5lur3+2\n8WR73bQd234PO44HGMUvqpmaTWrulAdQt9yVaABl/YD2fpMCjIZ5osv6gJo5Lc9UU4ECUDe+nQPm\n48sk/RdJl019IFimp019AGtuW9K/knSlmnlS/2HawwHgcLaarzs/W+678AFM7yvVtEV8uZrWhXe5\nFwcAAAAAAAAAAJibRff/Pe/ZrzDNCblIH/jHP1z0uYD6Xf/mt+2cvcVtHK13vupGxiyq9ZLfeTvX\n1w7G63pa7BWJAL0X/x6oHQEaAFATrkoAAABAJEI0AAAAEIkQDQAAAEQ6beoDAAAAGOrSQ08EL3v3\nceIPhuMsAgAA1YsJyUO3RchGCM4SAIuW88JswsUYSFN6bMbsm3EME84KAIsy9oWZizEQZ8rwbNMe\nE+MXXZwNIzp42bVRy5+867ZCRwIsTy0XZi7GgFktY9SF8YsuzoICYsNyzHYI1kCcWi/MXIyBXbWO\nU5tLDz3B2AUhOpdcwTlmPwRqwG4uF2UuxliyIeP06ovzPUn76In4h/YydsH//YHGCs++fROogV1z\nCdAtqtJAmJzB2bbdmEBNkF42/s8nmjI8mxy87FqCNKA8ATr1Qp1SzerigowliRmrpcKza1+h45lx\nu1z8X09QW4ButcdFmAbi5LpA97czNFQD6yo0QI8Znm37ZhzDhsd+R6o1QHfN4RiBEmKr0FdfvFH0\nIp2y/bm1ogClTBmgu0KOg3G7TFSiI8wpnNLegaWp9avh7v74ehgIG6tDx+j1556358+3PvTgoO1d\nffEGFWnsw7t0oDkFaAB2U389zIUYcIsdo/3AHLNMTLj2jV8+/C4P/7cD5ArQG5uHo5bfOXbHoP1R\njQb2quHrYYI0YBczRkPCc+g2QsM04xddhOjCYoOzbd3UQE2QxhKM8fVwTiEXYqpagFmO8Gzb5tC2\nDywL79AeKVXoIcE5ZJtDK9TA0qQE6NgLdezFl4oWlsj1gTdknJYI0P3t+8aya+zy4XdZmJ0jsxIB\neug+6OfGksV+Pdz+ipWyXk3VcaB2pQP02PvB/BGiHWLC58bmYY0RoKfaH1CroZWtVq4LZ2oIN2Ha\nLCyFb6yOHWx9++MDMCRCdBZThtnQfVONBsxyht6U7XIxBtymqgxTkYYPIdoiNHTWUA2u4RiA2tTQ\nXznWPoA5c41Vxg9qRvf7ADWF143Nw9xwCEQY8+Lsu1mJmwyxBGO0J91w1pZ3mVse2Q7enmvsMm5B\nJTpRtgC9dWHzawS0dGDdpF6Up6hupe6TvmjA74aztoICdOyygAuV6DG5wrLpZ9v3Rm2eajTQmOPX\nw1S1gL1Cx2pqIL7hrK2gqnTItHddTHO3HPxfNvBVbKOq0EOqzN11AwM1QRooo3uhjvk6uCv2YgzA\nbWhFOTRIAyaE6EijBWjTtiIr0yY8wRBLFluFtl2g+69zEQbSDJmdhpYMTI2e6BJK9TkHbremGx6B\nseWYMi62ZzJm+VrbSYCajDlOCONIRYiOEBROx7hJcKQbEYF1Urq/cui6EnNGA7WyvX8wZpeNEN0z\naAaLMcOtZ19Uo4F4OSpSIdugGg0MQ/UYNSBE5zJFdZiKNJBNzosyF3hgXhizSEGIDuSs7E4ZZh37\nth0z80VjaXyVXy6gAIBYzM4xgdO3znD+/PHtx+I2mGnmDmBOTA8hqak/0Td1FtPdAcC8UYkOkKsK\nffrWGd4A3V0uZFkfeqMBt9qq0KYPAjy1EADqQ4gewZBAHLwu/dFAdWoL6ADMmOsdKQjRHdG9wgHB\nNUc1Oed2gKVx9UMTcoF5IvSiBoToVJ4Anasdo79NJ6rRAABkZ7t/4eiJnZGPBDUhRHuk9BSXrBqn\nbJu+aGA6VLsBO1sIDbnpNlc1mqo2UhGiZ8gZpAOr0Uxzh6WrIdzy0BVgWgRoDEGITuEIqmP1LtMj\nDQBYuiEhmACNoQjRGaUE2/MPnXrqV7b9GUI+LR0AgDkJnUf9lke2owNxzPKx87nffZxHcCwF/6cd\njMHTUoWOCdC2wNx//YHj/s84p2+dEf9wFgAARnb38dOKznneBmNXqxbVZ+REiB5RbLW5XT4kTO/D\nUwwBADNw9MSO9Wmjtz70YPS9AzmDsqsKzcwcoJ0jA18VOrVdo7v+kP0DAIA4sW0cWB5C9ErQbBUJ\n8zAPCc8ltgMAwJwQZlErQvRAripw7uDrqmgbj6MX+rm5EABQI19rxNhB2rc/WjkgEaKthgZOKscA\nAOQzVpAesh9m5lgWQvQMRVWjHXjgCjAtvqbG0rhCZkh199aHHiw2bkK3TRUaLUJ0KEM/tC20jlGF\nDt5HQh83AAA1yxmmY7ZFgEYX3ztkFhpuL7/Evdyd9/g/35x/6NS+6e+YNxoAMFeu6e5MuuE3Zio8\nvgVCDoToRKnTyvnCc3+5kDAdY2PzsHaO3ZF1mwDsYuesNVW66LPEOvE9dCU2SLdKB2NfFZpxujy0\nc2Tkq0KHBuj+Oq71TPvcF/Bp6QAAzEhNbRNHT+xUdTyoByF6JCkBOnR9ZgIB4vH4X2A6IVXbGsJr\n6P6pQi8TITpEQCXXFWSHBuiU7fAUQ2B6rqBOTyYQZoowXUOAR/0I0Qa+OaJjAmquAO3bHtVoAMCc\nxFZvxwi2KfugCr1c/J+focsvOWW84dA0W8dTti6Utu+VxM2FwBhoFwH8fDcZmvRDbspNiKbtpCBA\nLxv/9zOwVYFzV6H72/bN3MF0d0DTNmGb+uqWR7Z1w1lb2ffpC9C2Vg6+PsYSpQTprqnGDQEatHMU\nEhOgr7zgaXt+DdlHbFsHTy0EAEzt7uOnzSaUzulYURYhOlLOG/ZsoTkmUPvC+p7jZao7YJ/cbRep\nVWgbLtZYkpoDas3HhmkQon08wdNU+fUF25iKc2x12nZMfb6bJ4ElyRWkh2yHVg5gVy2BtT2OGo4F\n9SFEjyw2EIesF1WNBtaE6aJmC6Ih1d+hQTpkfaa1A+J0Q2zpMDvmvrAeOENGlBqg++u//74n9/2s\nf6OhdaaOziwdAPZKudEwRxWbKjQQjnCLWlCJzqzkjBytlDBONRoIExqKb3lkOypAp1ShCQsAUC/e\noZU+Q0XMTBhDq9Cm7fUr0rHVaOaLxpK4prrry32zoStAU4UGgHmiEh0htZqbO0C7ttuvhHeDvu34\nmeYOKCe1D5oqNADUjXfpjMZo5egLqUgbGarRz3v2KyiJrXzgH/8w7RFYmNzREzvWJ5jFVKNz8AXo\nIVXo69/8NsZrx61veDljFtV6ye+8nfHa8c5X3bgW45UQXVhMFfq6c56x589HHr4/yzF02zpMTzHc\n2Dysg6vfn7zrtiz7nLvnPfsVOwRpDDEkQFOFxtx9fHv8olLNzt7ii/91xDt1zxTzJ/fDs+l1V6CO\n7Y9+Kkh3ZupoK9Ku1g4CNtbBGNXo0gGagAIA0yNET8gWnl3L2sL00LYOafcDhO1mw9DeacI2alcy\nSJeeC5oADQB1IEQn6s/MYeqHLnFD4XXnPCM5SFvbOtqnMvbCdCt2Bo+QsE3QRkmuvuiSQgL0kCo0\nARoA6kGInkhMFdq0bkyQ7uoHaUnWMN0aGqpNaBtBDncfP02XHnoiad028OaoSIdWnwnQAEpIfR+0\n4b6MMPwrubShMrMhAbq/DVOY7gdp3/zR+242tITplqlvPOd806aATbBGrNBq9JDWjpjWjakC9Knb\nx3/U+IGrxpsBBcuV89ye0zmbOzCH7oNgvR//IjNnq0qnBGlJ5jAteR8VPmawJlAjt5iqdErP81gB\neorAbGI6jjmFFNSn9Lld+zk7RnCOOQYCdYN/hUJs/dA5qtCmbaYGaUn+MC1FBepWqWBNoEYpJW4K\nHCNA1xKeXfrHWFNAQZ2mPq+7+5/ifK0hONsQqBtMXLgmbOG8H+ZNN0CaHl9++tYZ9ic0bl2491eE\njc3DT/3K4eBl1/LERVhN+Ujtoyd2igfoU7c/OHnQSDXnY0dZNZ4bYx9TzQG6b07HmttyPz6soZiK\ntCRvVVpyVKa7EqrUUt4bFtsgTWV6eXw3F04xU4cvvA8N0CUu5hvbx4yv72xtZt9XV/t3oTKN2oKz\nyRjna45QGvueN7TgcOmhJxZZkV7e3ziRtSoboUQrh2kfIUFaMs8j3e+VbgWFaSk5UEv+eapDHLzs\nWoI0JlVrgLaF5Bzr5QjahOllyxGgY87xoefsqdsfLHKupgToHEUC0zZig/USg/Sy/rYzc83BKyRJ\n7z354aj1cgRpaX9VWtr7YaJUoB4apqlKo2+ManTIBWfMAJ0amlP09zUkoJQKJ6jX2B8MTeumnLO5\nz9XYAF36Pa3dfkyYXlqQXs7fdEba8Gz7c0iojg3SkqLCtBRRnZa80+aZ5AjTBGm0Sgbp0gG61vAc\ncgw1hBPUKzZAlzq/U8/ZXOdqTIAeuz0tNkwvKUhzY2EGppv1UvUDc+oyUvjNhi3b3+P8Q6eMNx+2\nnDch9g24GTEFNx0uQ+gbdu4bDX03D7bGCNAb28eqCNB9tR4X5mXM86jWc/bqizcmeRJrLfuvESG6\nIqHhuF02ZPlcQVqqI0wDQ4UG31zbKB2ga73g98Ue5xxuMsMwtZ/fofsdeq6GVKFrCq8hx7KUGTuW\nUW9fY9ccvMLb3hHT2iHZ2ztatpsPWyXbPFJaPGjrgEkbgkMvTinBe4wAneLRY0eS1muduXld8rob\n28eCvy6nrWPZavhw2B5D6VlqXFIDdOiTWFPmx7/64o1JpxCtBSE6gzvvOZC1pSNWiSAtmW86bPn6\npaXyYZogjS7fVHc2pS4EJQN0bLgYGpp924sN1TFBGusp9zneCjnXUz4E+s7ZUh/4YgN0aHC2rRMT\nqAnStHOsjaGtHa72Dl+Lh090m0cg2jtQo7uPn1ZNgH702JHsATrXfsb6qhzzk/IhMeYcjF2+pFxt\nD9efe15SgB66HVfIX0JLByE6UFA1dWJDgrRk75OWhvVKtwjSKG3qO8JD9j9GgJ4qIMTut4av6zFf\nOc7zmG2Mfb6GVKFzheextrtuCNHYIzVIS+FV6SCRNx2GYKYOlFRDgK6lulZzMMH0XOf6VOd56DZd\nx5fzW5PQAF1ayD5quulxbIToNTO0Gi3NL0hTjUbX2NXokPYNaZwAXZtcQZqWDrRKn+c1jiOTMavE\nVKTtCNFrqIYg7QvTBGmsg9DAXnOAPvGxm52/hppLKMH0fOf6WOeSbz+5vj1J7RmeItSm7nPd+6IJ\n0S4RT9ebI1+QdskxG0lwkM6Ilo5lKF2NDq0+S/5HebuEhIqUYBETknOE6qFfkwNL+jDmao+Ysirs\n2vdSWzqY4m5ERx6+f1BwHZtr+jvJPQWe5J9PWmqCtPemza0L1/4DDfJLnfIuZLuhfAF6SG9obKjI\nUVXubufii14Utd6jx44Mml8aGNM6na83nLVlfP2WR7ZHPpL1QyW6EFf4tPHN9Rxj6KPBW662DmnE\nHumAtg5aOtCXsyIdU32W6gnQudoyxtou1pvtvHed80O/bUk9V137neKbk9gq9A1nbVkDdMjPcxzD\nuiNE94Q+wMNXYcWukCANlBIbfk3rxq4/JEDnNEbIjdnHWL2mWC5XYJ7DB78cbRGx4TglTJsssaWD\nJJiJq62hy/TUwK6c1ehchlajQ0zRH41lCQ3EqcG5NTRA56pCjxkWxtgXM3QsU8z5Hnoe5vzg1zf1\neTokDOcI0ktDiK5QjiAds40x+rSzVKMztXRwcyG6Qbn/a4ghNxFK8wzQsftc0g1iGEfK+V57Rbov\npI0iRwgO2QYtHbsI0ZJO3nVbke2m9EW3aqtIj1GNBuYsJEAPqVLVHKBr2Dfma6o2Hs5XDEWIjpDr\n0d++lo5WapCuLYCHoqUDc5UjQOcIEoQCrJOQD45Dz/l1GTM5WzFo6whHiB6gf3OhqS/aVo2OCdIx\noXiuARqYqzECNC0QwHRs448bYcE80TNRMhyHBvqhQuaN9mLOaFRkaA90TrVU1E587GbvHNKuOXg3\nto9pZ2uzwJFh3dRyzmO5qET7RAa2EtVoAPUJDdBUoYG6EcaRihAdqd8XPbSyOnWQDt3/kJskgXVT\nUwVaIgQAwBQI0QXEVKOl6YP0WLI8oIZWDkwsJkCPcTMhAGAa9ESP6P33PWmdKq4N0mPM2dzdXwiq\n0EAjZ4AOQSsHUJ6vhx+woRJtEPro75apwmp7gqEvkI5RlS6xj9AnNgJzVVsLR61CAontpkJJ3FSI\nYIRfTI3kE6LXQjB0vuiQIN3+yillm2NWoXPNww3kFhugp370L4B8bB/6avrAd8sj21MfwiIRojOJ\nqUZL4eF0aKAesn7oMYZUobP0QwMTKFWBph8aaNjCqOsbi9bQajTV7P0I5OHoiU70+PZjQU/Yu/Oe\nA9ZHYrt6pE3GugExpvqcK0AHVaEDbioMacUp9Zh3rJ+5tHBcfNGLqpihg0CCKaSe/+t2vt7yyPbg\npw0SoONQHgwVEOBsYdFXka7pxr3cARqYq9QAnbOVI6QS15o6EITuP+bvBISKPf/ndr7e+lA9LWI1\nHcvUSEEWIRVNU/U0JUhL04fp2P2HBuhsVWhgRHOpQNdijAB/4Krziu8D9SnxQbLk+Tr0PD16YmfQ\n+kMqyVSh49HOsXLyrtt08LJrs2zrgeMHdP6h/RdhV2tHqw2yMW0eqVJCe0z1OWsfNPNDYyRzDtBT\ntHXEBBJfIKrpRi2UceCq84reeNs9H7tjYepvasbUhuHQ1o5c4XnoB4A5IkTH2L5X2rpwz0u23mhX\nkJYUHKZbOUL10Ep3iQCdqxdaip+aEOirMUCfuXld1HzRbVgYI0wvKZigrJ2tTevNtrFjoDX0/Kxt\nKsZbH3pQ158bXunuhuN+oE4NzrRy7EWIziA2SEvhYbo1ZatHqeozbRyoydgB2hUacigZplPCSS29\npZin1CA9ZH81ig3SrRzVZgL0foRoh51jd2hj8/DeFw3VaBdXkJbiw/RYUm4azFp9bmWuQjMzB0xq\nrEB3DQkQucL0kKpeSCChlQO1mDpAHz2xo6sv3pj0GGItsZVDIkSniWjrkPxBWqojTKfOtlGs+kwb\nB0YwZYCOqUYPrcSZQrAtWOds05g6kGA+fONhjGr0mB/47j5+mi499ET0eqnV6CFSq9B3H1/vmLne\nf7sMjNVoC1+QlhQcplslQ/XQKepibxwsEaCBIWqvQPe1F/hcQaJ0T3NogPaFEmbmQCv3GDBtO1XO\n89RXjR4zSNPGYUeI7oiaocPS1tEGxaFhuhUSdPtBu+T8zSkzbkT3PkcE6JgqNK0c6KolQKf0Ro/d\nHxorJozQxrE8rhk6QsdD7jGQ6wPfmMYI0iEBeqmtHBIhOoi1Gu3oj/Y90bAbRkMDtU3ph56kTlVX\nMjxLtHEgXS0BupUapKUyFblUsZW8mgIJ6hETpKX0MVD7+RrSG92G3NxhOrT6vOQALRGih/MEacle\nlW7lDNQ5DJnfOXnGjcIBmio0WqUDdOo8uKmzddQQplO+Bg8NJLRywKd//rnGQmrLRqkAndoX3Zcr\nTOds3Vj3fmiJEL2PraXD2RvtmbEjNExL5gBbKljnfBjKWOFZogKNdLVVoPuGTHsXEyRyGNI/SgUa\nvg+b7Tky5INlLiHna6kPe7EzdfRDsC9UDwnNS69CS4ToKN4gLWUL011Zn/yXyaA5nhNvGkwNz1Sh\nIY0boIc8lW1IeOgyBYmUYD1FIGlRhUbp+dRD9j+1IVPelbopkADdIERH8s7WETCPdD+AxobqKWR5\nMMqAGTcI0Bhiigr00Mcb5wrTXVNONxcbRgjQ6y90jEwVpMdqOQpp6WhD69TzR4eG5yW0ckiEaCPf\nLB1BQVoKfihLTaE6+1MEB05VN6R1gwANadoWjqFBWioTpsdUQyUP8zfWOKj9w96UD2Kh+rwfITpR\n0PzR3QAZ8ZTDmCDrC9yjP1o7w/zOOXqeCdCQ6uiBzhGkpb0X99oD9dDgTBV6OWLHR6lxMPWHvZgb\nDMeuSseG56VUoSVCtFXInNFt2At6GEtioPYZPSR3ZX4gSq4bBgnQkOoI0K02FOYI05L5gr8ufaME\n6OUZMqONiWss1Hyuxs7U0Q23uQN1atV5SQFaIkQ7hT58JSpMS+bwmTFYZ1fo6YG5Z9kgPKNVU4Du\nyh2mu1zhIEfALl2pIzwvW65vbKR5n6upU97ZQq8vXOds0VhagJYI0V4xTzGMDtNdIUE1Z9Ce4LHa\npaamIzyjq9YA3dW9CJcI1H1Tf1XtQnhGq+SHzFzGOF9zzR0tjdfHvMQALRGig0Q9Dlz7w2JSqDaZ\nIPimGmMuZ8Iz+uYQoPvGDtS1IDzDpsYwPfb52obSXGG6lKWG59ay//YR2sAWE6ZbpkCZLVhPZIoH\nnhCa4TLHAN1nulDXFCSGIDQjVv+cGXMs1HK+1hqmlx6eW/wrROoGuZRA3YoJoaUDd61PACQ0I9Q6\nBGgb18X6D9e/AAAgAElEQVS8xoBdS/jA+rGdW0PGwVzO135oHTtUE5rN+FcZIFeg9qk15OZCWMYQ\n6xygfeYSAICSljgObKF2aLgmLMfhXysTWxAsGa7ngICMkpYcoAGgjxA8Lv61C0sJkVMHb4Iv5oAA\nDQCYEiG6QoRYwI0ADQCY2oGpDwAAYhCg0Xf2FpcyAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAYPE2pj4ARHmNpEskvUDS70r6rdXrXy/pUkmfkPQBSb9nWPdD\nkj5Z0gOSvlvSr3u2KUmfKem1kr7LsL0XSPoiST+5+vMrVts/U9JxSW8N2E7o6wck3S/pVGeZd0j6\nqoDjBKby5ZLOlvRcSf8k6Zc6P/Ods64xbVrXNo7PlvQ9ku6UdK6kn5W04zg22/uE7z0mx/sBUIP+\nuSyFn7PddX3XLdO+bOPvKyRdvdrWXZJ+Y/W6bVzaxp9tedv2be8fQBHfLuknCm17S9K3rX7/TDWD\n87Ckz5D0rs5y71AzcPq+UdJzJJ3m2ebm6s/fqWbg/ZphWwck/bGkN67+fLn2vrn8iqSne7YT8/ph\nNSFhU9IhSd8h6aqA4wRcSo7X8yU9qmYsbki6V825K/nPWdeYNq3rGse/2tnv+1a/dx2b6X3C9x6T\n4/0A8Ck5Xlv9c1kKP2f767quW7Z9mcbfeZL+uvPnv1Azzm3j0jT+zrYsf4Zl+xeufm96/0DHgakP\nYM38vJpPmRcX2PanqflEKDWVo21Jz5f0Mkl3dJb7qKTPMaz/uKSPSHoiYJtS84nz9y3H8mpJf6rd\nbzIukvTFkk5f/fmf1XzSdW0n5vXH1HxiPibpodW2bw84TsCl5Hh9QNJnqQmrO2ouiu148Z2zrjFt\nWtc2jq9QU3U6vvrZl65+7zo20/uE7z0mx/sB4FNyvLb657IUfs7213Vdt2z7Mo2/l0g62vnzbZK+\nQPZx+UztH3+PW5Z/sWX7Xyj7+wc6TvMvggg7kn5T0usk/UzgOleo+fRp85dqBvAfSbpx9dqGmq9u\nP6jmE+MndZY/U82n3T/rbecGNZ86z5X095L+wLLN7c46pnafiyQ9Kelj2q0u3armA9ktkn5Z0p+o\nGbSu7cS8fnfn998k6ecCtwO4lByvUlO5kZoL1Z+ruZi2XOfsw3KP6f66tnH8hWoC8+vUVJ8flvQW\nz7GZ3idcx5Pz/QBwKT1eTedyy3fOmtZ1Xbds+zKNv7a9o/WApE+RdJ/2jsuDkq6U9Isyjz/bOH7C\nsv1zZX//wAohOr+3qDnxQwf5hyV9X8Byn5D0d6vfv1zSeyT9jZqB+vVqBvnZkp6nZvD0/ZmaXkmt\n1nunmgFi2mbL1Pv0SjUD89/0Xn/z6u/xM5L+Q+9nth6q2NcvUPMp+7HA5QGft6jMeG29Uk3F6T/2\nXneds2+Ve0z317W9N7xc0qdL+tern71L0v9T8+Hbdmym9wnX8eR8PwB83qJy49V2Lkv+c9a1rum6\nZVveNP7a9qvW42qC969o77j8VEl/tVrGNP5M4/g9q5/1t3/26veu9w+Ido4SLpJ0lpobBko4X9JN\nam5ykJqvZL5Ozaftz5f0t6vX+rpfR92/Wta2zVb/0/cLJb1bzRtK92efutrel6i5oeEHJb3IsZ3U\n179ae78O8y0P+JQer29VMzbfrt0+Zcl9zvrGtG3d/jh+aLVu6yNqvpJ1HZvpfcJ2PP9Ced8PAJ9S\n49V2bWu5zlnfuv3rlmt50/h7uLfcQTVVaNu4tI0/2/K27T8s9/sHRCU6t5ep+Rrkx9ScrCfUDPYt\nNV+p/LVhnZivmzYkvUHSN0j6uJom/+Nq+pnar2ffKOmHett4raQv0+5dwU/Xbs+VbZvS/k/fN6h5\nA3upmt6rg2ru8t+S9NurZf5U0uvVfE18s2U7Snz9C7R7t3LI8oBLf7y2VZz/rP2V41boeH25pO9X\nM04+ruZi9SrtVtB856xrTJvWNY3j90n63M4yp9QUTmzHdo+a8Wx6nzAdz43K+34AuPTH6z9K+t9q\n2ofepOb87Qsdr6Zr25epqXpL7nPWt27/umW7jp4j83X6Q9q9V0lqqtq3rn7fH5dvXG3DNv5M4/hT\netu/cLX9EzK/fwBFvEbST69+f46kf1i99oVqLlBfk2Ef367mhqBnqwnnn6emgnTb6udXSfqdzvLP\nVXNxffHqOKRm8N6x+q9tm62bZL8j+U2Sfnj1+1dq79/vRjU3K/i2E/v6ETXTAYUuD9iYxuuZasZM\n/36CFC+T9OOr329ofxXnJu0/Z9vxuin7mLataxrHZ6gJCa2bV/uwHZvtfcJ3PFKe9wPAxjRer1Rz\nPuf2Ju2ey62bZB+vvnVt163+8rbx93TtrQjfJulZso9L2/izLW/b/pkyv38A2b1QzVQwXb+g5hPg\nYTX9Saf3V4r0YjU3Ipxa/XpS0mVqbhT4EUnfrObGhWd01rlV0nWr33+tmt6on1PzNaxrm5L0rWr6\nse5QM8jP7Wz3q1bb/ms1HxCkZvqe71/99/WdZW3biX1dasLNlb1/F9fygIltvN6k5iKWK+B9s5qp\n535GzY1FLds5245X15g2resaxy+T9KNqqndfG3BspvcJ1/FIed4PABvXeP33q/9enWlfpnPZN15N\n676687rpumXbl2n8Sc3NfT+opnrcjmPXuDSNP9fypu1L9vcPYFQvVPn5LQEM80I1N9RQJQXqd0BN\nJXhDzBqBidDfUtZPqfmE/JiaO2EB1GtTzdeuz1HTxwigXt8i6Ro1LUxkGUziaVMfwJr7ZzVPD3qR\npP+mZk5IAHX6OzXTxb1azVO7tt2LA5jQP6np0X2xpP8q6cFpDwcAAAAAAAAAACC3xU56f/2b38Zc\npT23vuHliz0fUL+X/snbdh44Tutj693feCPjFdW67nt/j2tsx5Gf+grG6xriigQAAABEIkQDmAWq\n0ACAmnBVAgAAACIRogEAAIBIp019AOvs1O35p608cNV52bcJAMBibN8bvuzWheWOA7NHiM6gRFiO\n2RfBGsjj8e3HktY7feuMzEcCIIuYwBy6PsEaK4ToRGMGZ5/usRCoAb/UsJyyPQI2MKKhoTlmH4Tp\nxSNER6opPJu0x0eYBvbKHZxT9kugBgoZIzzb9kmYXixCdKAc4Xlj+1j0Ojtbm0n7IkwDjanCs0l7\nLIRpIKMpAnR//wTpRSJEB0gJ0CmBOWQ7saH61O0PEqSxSDWF5z7CNJBJYoDeOXZH8LIbm4fDjoMg\nvTiEaI+YAJ0rOIfuI7VKDay7mgN01+PbjxGkgVSRATomOJvWCwrTWBRCdAZjhGfXfn1hmmo0lmRw\ngB55+iuCNFBWani2bccapqlGLw4PW3EIqUJPFaBjj6H2GyKBSW3fu/srdb0BfZlzqZwD1Qgcb7kC\ndOltYp6oRA+QGqAfPXbEu8yZm9dFHwvtHVi66DCa+4Yk7tYHJNn7/cf8wFgy7O4cu0PGijTV6EUh\nRFv4KrcxATokNLvWCQ3UBGkgUOm7+RMupLR1YAm65zjfwGDuCNEJQgN0Snh2bSckTLuCNL3RWGfB\nF+SxpsOiKg04tYG6RJiOqUKfvOs24+sHL7vWuw9uNlw2eqILyRWg+9sssV1gMaaeT9aDyhyW6PSt\nM+K/hcnw4fTkXbdZA3TIzwFCdKSQKnTpoOvbfg03OwLVmSpAVx7cgVqM2c4UE45dyxor3oz5xSBE\nGwyZyWKsSjEVaWBGuKgCQZKq0iOgIg0TQnRGNQVbqtFYEm8bBCEWmJWgIJ3Y0kEgRi6E6Ai1BdOU\n0M580QCAORgSpEvc8GcL38wbvVyEaAAAUKUSrR2+WTeAUIRoAABQLW+f9IjTSNIKgi5CdCZT9UPX\n1IcNAEApsVVp5nBGaYToTGIf0w1gQXjgCpCFNUhH9kbT0oEcCNEzR3gHKhcRoGuc2guoTWyQBkoh\nREewPU4bAACMJyZIU41GKadNfQDr5MzN6+hRBrCLKjTgdP6hU86fP3DcXus7fesM/xzxQEFUog0O\nXHVe8rpjtlfQygFUjK+WASdfgG6XcS1n/PBJNRojIURHCmnpqCHc2o5zyAcEYLbGDLRbF0bvjyo0\nliYkQPeXt61DkMZUCNGFlA7SNQR1oBYlHxEcJWEfBGgsTWyA7q9rWn+sIE3gRhchOkHoDYalgi4B\nGqhMYvWZAI2lGRKgfdsJDdI2BGTEIkRb+NoeYoJ0rtAbui1mEQEsclaj2+BM9RmYRHCQ7nE9hMUV\npAnZ6GN2jgF2tja1sX0saNl++I2ZxSMmhBOgAY829G7fm75uIsIzlsxXhb78EvvP77zHXPNrt9md\nxWPfrB1bF+4b7xubh7Vz7A7jNg9edi2P90YQQrTDgavO06nbH3Qu04bW0DDdmqIlg5sKsc6ip7vq\nBmJboM5UuSY8A3au8GxaxhSozz90KnuQBnxo5/AIDZ41VIBrOAZglrqtGQPaNProewbcQgJ06Dr9\nSve+sVfo5mJXewjWGyE6QEyQbn+NaYp9AjWqIbC2wbmGYwFqYWrlSAnQ3XVN68cGaQIwhiBEB4pt\nhRgrUIdun1YOoByCMxBnSID2bSd2BhCCNFLREx2hDaK+Pum+ftCN7Z/2bc+HAI0lGetRwKUCc64p\nwIC5uvKCpxlff/99Txpfv/ySU9YbD6Xh/dGADSE6QWqYbtF6AZRVKkjnCs4EZSyZrQptC8/9n5vC\ndD9I577REDAhRA/QrfCmBurSqEJjqXIF6RzBmdCMJQs5/30B2rRsP0xHB+kMaAVZNkJ0Jv2wOnWo\nJjwDuwE45sJJaAbGFROg++vFBuk9qEZjIEJ0IbYQWzpcE56B/ca44Y/gDPj1WzlSA3R3fV+Q7hql\nraPQVHqoDyF6ZCEhNyZoE5qBaRCagbKuO+cZe/585OH7jcuZgnSXsxptQUUaIRYbos/eagbUx7fr\nuxASjIF6EZ6BPGxV6H54Nr3eD9T9ID20Gh2CfmgsNkS32jCdS42hHMAwuYKz79HFwNLZArRpuZgg\nTTUaJSw+ROcWEsoJ2kD9hgTnkAdJmJYhWGPJQgN0d3lfkLYJrUYTpOFCiJ6AL2gTsoHppITn3E9f\nI0xjHeUaJ12mIN3f55BqNOBCiK6QK2QTsIEyYsNziUDQ3TZBGuus3w8dW4Xur9sN0qHV6FCmajT9\n0JAI0bNDwAbyqik89/dDkAbKyXWDIZaLEJ3o0kNPDN7G3cfz/vOH3iRJ2MbSTdmyEbtPgjTWnasK\nfc3BK/b8+b0nP2zdhq2tI0dLB73RMCFEB8oRmkO3mTtc9+WekQSYi7HCc+gDJHJ+5Qysm36Abl+z\nBemu3C0dgAkh2qNEeI7dZ+lQDay71Jk2YgJ0ypPX2nVcF3uq0Ziz1LFnCtDdn4UEaaA00pnDFAHa\npHscBGogXOlp6oY+sri/HSpngDtAd5fpB2nfTB1AbiQyi9gAffXFG9H7OHpiJ3odAjXgVzo8S/kC\ndH+bpiBNNRrrpMTYSbXv5sJAO8fuEDN0gBRmEBqgU4Kzb/2YYN0eJ2Ea2DWX6rNr+1SksRRDprYb\nU9RNhdv3NjN9YO2RvhIMDc+h2w4N1JceeoIgjcXL9Whul5oqaACAaZG8enxV6JIB2ravkDBNVRpL\nlSs8+6rQMQE6pLrm6t2kGg2MI6WVo0VLB2iyizBmgO7vN3TftdwMCYxhrAAd6rpznhH89bRv2X5o\nn2KeagCRrRxYFMqWgUJC7PXnnhe93VsfejDqGEKr0lSkse7GaN9o+arQOR9ZDGA+qEYvG5XoDlsV\n1xegrz/3vKQA3V03dP2pquFATdYlQPu2QQ82MCLDI7+pQsOFEO3hCq1DwvOQ7YUEado6gDC1tEnM\nZZYCIKd1+BaGoL1chOhEOcNzyrapSGPJHjhex1tX7uDr214tgR8YwnfTbMjTCE3L5A7kg8KxoaqN\n9VPHlWhmSgboMfcBzFktQTq3fpCmpQPIp/u+MWRmjj6q0cu0nlehTKau9vqC9NTHB0xtyiBN+wXg\nlzpGXdXokEp1t9od9LTPXuWYUIwQTOGwEtpDHFIhvuGsLe8ytzyyHby/mBk8gKV54PiBUW80THHN\nwSv2veYKAv0ZO5g3Gkv03pMf3jN2XGOmht5qZupYHirRmYUE6Ha50GUBuD1w/EByxSuoSpXomoNX\nGAO072eSvdJNXzTWjSsAv/fkh5/6FSrkA6erlWNIFXrPuvRFrz0q0RamVglXFTo1EN9w1pa3Kj1W\nNfr6N78t7DnjBmdvrd/nsXe+6kb6ZWZmDlXpvmsOXhEUELrV6MsvOaWX/kn6eE1Vcx/6u7+R8bpU\nrhDe/ZBsPX8Dw+7Ju27b8+eDl10btJ4kXfe9vzf6eK3ZkZ/6irUYr/W+Iy5IDRXpj2/PK3gANilV\n6ZLV6BC2irTrJkMq0sB4M3L0A3T7mul1qtHLQYjOoIYQPAQBGusoV9XU9tVw7ot3apDu/gJq1v2w\n2h9XucaT7YZC66wcASHXFJR9PydILwMhOlDpKeemCuIEaKyzmKr01NVoKS1IdxGqMWexQbq//NCb\nb01VaF+Adi23L0h3f2Et0BM9UK1V6LuP+//XEqCxFEN7pd9/35Ojzdds65E2zdjRHptNP0jX8EEB\ncDny8P3e6SNDwnaOKnRogO4u3++Tts7YYdr31oVR+4syJLiXPK6ZI0QvVM0Beugjy0M+QGB52gup\nK0zfec+BqApuyAU/RWiQlsLCdKv9uxGmUQPbh9P2HO+PLVd4LlGFThEVpPtqrVDbjotwTYieA9vM\nHEdPmG/29YXIGgP00ODs2hahGl2+qrQtSI9ZjZbigrS0t8XDFygI05iD0PaO/vkeVIX2iK1C99c1\nBWlJWqt5pLvheqGBmnfQSoQ+fGWomgL0pYeeeOrXGPsBWr4+6Zhw6ZvjdghXj7SrAn7lBU8LCvz0\nTWNs/bE1tILsCtBOnQBY4umEthC+c+yO9Xwa4kJ7vQnRM5VSha4tQC9hn6iX76ZD08V4rJk6ulIe\nyNIiSGMOUoO0b70pqtCh25kyTLf79v1KsrAgzffclYt5yAoBOnzftHig5WrvMLV2xLZ19B9dnML1\nQJY2SNuCfOgNiLR2oJSQG3tjx5XpfE6pQsf66H3v2fPnZ13wfOfyptaOrpJtHkNDumn94N7uhbR3\n8K5ZgdhWDlsV2maMAB26jyEB+uqLN/b9AnLIUZEu2dYhuSvSUp6qtEvqvNuPbz/m/YVlsYXdkIr0\n++97MihAh1ah+0HRVT3uB2jbaylyVKazVJID9+G1kIo05biJuQJ0jir0XCvQoQG5v1zoB4xLDz1B\nNRp7uGbvCK1Il5qto+V7RLjtpsNW99Hhfbmr0THhuF329K0zsu0f85TS2uEK0PskhrshYdlXje6K\nqUxP2Vu9ljdKJqASPdBYNwS2TCExR4A+dfuDe35NZWiFOWZ9eqTXX0q103YBHlqRzlGNloZXpEsb\nUl2mMr0cuT6w+baTqwrtEhKwY7ddrGc5M+dxLKAaTYieUGwVOneAdoXm3IE6JLDmbM+g3WPZTGEs\npoUgJkjHqCFIl5ymL1cAJkgvw9DxZFq/RBXax9cXPQcn77rtqV+xagn0UyBEB3K1VsRWo295ZDtb\nG4eJL0DHhuMxKtOlAi9BelliKpm+QG2bvSNkiq7S/dFSfRXp3MGXIL1eUnvqTe6850BQgI6pQpcW\n2s4xlm5o7gdn18+wFyE6k9Ag7Vsu5sEqpiq0K0APqSyXDNKlg65r+7R0rI8hocsXpvtqCdKpbNXo\n2qa7I0ivP1sgdi1v4g3QVKElKSkYh6yz1Go0ITojV0D2VZ+l8gF6qNRtuILqWJViKtIIMbcgXUs1\nmrCLoVxhuv1ZaIBO4QuJz7rg+U8F5vb3oQE6tgq9sXnY+SvV0Kpy0vpr3hfN9AQRbn3oQV1/7nnO\nZVJvNCwVoHNXkE/d/qAOXOX+NyjB9e8+tP0F85czxNlmijDNdduftSN2xo4cc0iPwRVSSgfox7cf\nY9aONREyX3Rsn7Tp3PRVoYdUTWMrz6EBOiYct8tOUf11zTSyc+wOLW22DirRFrFzMQ8xlwCdW0h1\n+Ppzz/N+cAlZxrU/WjoWon0sbeDjaU3h0NQnXUNFGliioAA9oZAAPaS6PFVgpU96FyF6YlMH6I3t\nY0/9CmXbdu45qUOCcX/52HUwf96Lpi00B4RqW6/0ugbp2CrgWIGlpmCEYXLdYBi8nQnaCQ5edq03\nQA9tzUjZTm03N64DQvRK6IM3crYOlAjQITcP2oJzSqBO5atCDwnDrnXpjV6QwIpzyPJjBelawnQR\n/Q8ta94rCbshQdo2a44U9mGrdAvEWOHZtF2MjxCdYGiQvvWhB4sFaJvYgOxbrmSrCNVkhLJeNIcE\nNEvAM1WlcwdpqZ6qdGtIYJHkDswE6cWKDdKu8CxZzseI82toldZXfS4Vnvv78MlVjba1dCxtlg5C\ntIOrLzo1SLvWKxmgU4xRke7LFaAJ4guWK5g5wnRXqSAdEqZrC9x7hFabB/SoY958wTh0manOjTY4\nZwvPWxfafwUaM0iD2TkGCZmto13Op0SAzhGCN7aPaWdrc/B2pnD9uecZ/+2vvnhj1BtHURdTpcR5\n4dm+d99FrD9jRH/WgaGzdrT6Ifmag1dUEZyD+tCBQEPaO3J+G3Xwsmu9N81lnW0jIhzvW9bx9wuZ\nvSPk7wo/QrTH0RM7zj7aHD3StQbo7rZMQTr3dHdUjzGY48LiuqB0f2a8+LXb7VzIxgrSXTUEaK/U\nAG34sAK4lKhAD63SesNzrnO83Y4nTJcK0lSzG7RzTOjoiZ3sAbrUjYE5t8nNfRhbTJ/ezrE77Mv3\nLli5Wjt87R2pbNs1HYdJzkc1B6GCjUDOAD3BeeRs20hoywjm2a4v1BOGh6ES3XH38dOMcwb7qtEp\nbO0EQwO0z6PHjlh/dubmdd71fT6+fUpnb9Xz2czW0oH5C32sb+qNLu16+y5CvYppjoq0FF+Vzi10\nervawguWp5b++BxV55gHCTn/3o7KtK+9ow3SoVVpXw/4ktSTdhbCVn2WygXoR48deeqXi+/nU9xo\nOBYeuLKectwpbtxGgYq05L/hMEap6nZxBHE41NCP771Z0FEdPn3rjD2/YgSt69g3Ven8qEQHylGN\ndt3MVjJAx2iXD61KT/UYcMDHF6BNVRfX42yl3kWoYEVaUrGqdGgrRzTCLwqauvo89EbBEo+ub7dp\n/LfZutBalfb1Sbf675GE7P0I0RHaEBwbpn0zQZQI0LHheWrcVIixuL6y7P7MdMHYOXaHxgjS0rD2\njtgqtKmVo8ZHKvf/fbEMUXOSe/hCZH9Zr4HBufv+4OK6P8EaphODdIvQ7EeI7rH1RXeFhOmQKdRs\nT0msIUA/euyIsRo95ynvsMZ6FwrbBSLmTvSTd902WpCWlKUqPVkLB1VoFBD1oS3iHHT1CEfN6Wzh\nCs+hodm3nilUn751hjlIS/v+fUKmwYu1tH5oiRA9yJC5hmsO0N1t5bjZcGqmmwuZK3p5UqZyatfp\nh+ncQVryV6VbpkAdGp6LtXLkxnR3izbWNx5JoS8hPKcGZ5d2m/0wXaoqDTNCtEFINXro9k1qCtAx\nxu6LvuGsLUnSLY9sj7ZPzNvQhwqYqtJjB+lWarXZFqAHt3JQhUYGg4LzGOdggfDcH/8uttlzXGE6\nJkhLw6rSzun91hizc4xsygB94mM37/vlY9qub5YO09+ly1QBZho65FKqqmIK4vv2NXDWDmlG1WJg\ngMe3H9vzq1oBM230nX/olDVAX37Jqad+xfCtZ9qf8fgGzN6B/ahEW+SuRtvCs1Q+QLvC8omP3ayL\nL3qRc/2hcv9bUoFGDFcV+qP3vcf4+rMueL51W0Mr0n2pFekYQ6vQQKzRgnGpKnTGynNsYPZpt9cf\nv1NVpZccvgnRDrnCX44AbTMkQHeXKR2kgdrYAnT3Z6YwnRKku/ptHZI9SEvmGw5j5Khs08qBxfC0\nH9gqzya+8Bw6tm1j2BWmhwRpKSxMLzk8twjRHkOCtCs8S/62hy5TFdoVoEPCc395W5Ce+w2GplYR\nbipcNleA7i8XGqRdfP3RkjlIS+lh2heeQ59QCCxCQniWzAHaFZ5TPhR31zGNa1OYzhGkpYFBec37\noSV6ooP4wrBp+dQAHdrGkTNAA7CzBe5+m8jQ/mjJHWzff9+TwVXl1OozrRyYlVzfhCRWn2MC9JUX\nPC1Li5ZrO/19m47R+GHA0ScNNyrRgbqhuF+Zjg3ZQwN0KbR1YJYcF1JTP3RoFbq/jq1PuiumrUOK\nq0i3crRnxFShq77pCxhihOpzSHC2zQPvmomn3W7//eDyS0552zuMFWnJW5XGfpQdErSV5pCKc19M\ngLapoQrdD/cxx2/DDB2omSl8x06dFxpIS7ZaZNs2F1vMVUDldWiAdlWMrzvnGXt+2YQsY9qH6XiC\nKtJSvor0QirbhOgRxfRAS+NWocdWW09y7IchoFWirUMqE6Rd26SVA2svsG0hR4A28QViF9e6RYL0\nQkLwULxrjsQVoOdWhQ4R+4GhFKrbyCmlFSTEmEEaWJyIUFgyQOdgC9OhQbrPNf1mcpBeUADnHXoE\nKQF63arQoZVeQi+iVfiGHVuNlsoH6ZQq9GT90BX+P0WFfOdJZEW1RIAeUn12SQ3StoeyWMVWpRc2\ndvkOu7BaKrJLRCDHwcuuHfzI7xCxU95J/oewdLUBOOWhDVSzsdb6N8MlhrhSAdrnmoNXWH/23pMf\ndq573TnP2Hfz4ZUXPM17s2Hw9Hdd7b+r7V6IhYXnFu+uBfkCdGwVOuTR3kOMPTOHrS+6dPitrR8b\n4wqZZaOIgBvxfH3Jd95zICoUhyw7p17o0A8dWJi2WjqjAH3NwSucATp0mdCKdF90RbrV/bce+O++\nDvY/WZMAAAXKSURBVObz7jkz616B3tnaLLr9oUGaKjSmFvLY3Ji2ji5XmG5/VrwCXerCueALMsYX\n+wjvvpS5n33BeOjyJqHfYvFBNQ4huoCQAJ1jSri+IZVk17pTPa0wNQgToOEzWTU6UGh1uBuYRwnO\nwBqJDYyhQdRVhU4NxK71clajJYJ0DN5xM5tjBXqsNg7TzYW+1orYQOxbnlaO9dd/TK2tVzklSA8O\n35XNrVxdKwdVaIzEFRRDq9BzEzLtXYsgHaayd9B5Cw3QqVXokIrwxRe9KDgUhyw7VRW669aHHvSG\n45BlgL6YUOxaNvamwi7bzTzVBVwTQi8WImY6u3VBkPZjdo5MaqtAt+HYNHd0bY/2PnpiR1dfvOFd\nbmhItlWhedDK+nPN0tGGY9sc0FO2fpgeC55z21VJfIIcEGuJVWiUQXrIYMwAfebmdVGzdAwJzLYq\ntOmmwgNXnZe8Hyk8SANGvWmuNjYP77uxzzfdXWpYtlWh+20lqUoG6Sz6U4ylbgMYAR/G9jNNedfy\nTn23cJWVIuZnigr0GC0WpfYxVdWXXmhIw9ouSm5v7ItU9ir0kGmuCNBYEN/cz7nX60uZa54PHnaE\n6AFKBeiQ6eNKBmnXtmOmtov99ykVdAnQy2SrBOcKvq7tWKvQiYGxutYLm9gnmw18FDMQo5bzKDYQ\n5wrQQ9Tyb1cb2jkS1dAD3YbdXA9hqeEmwtxtHb4ATT/0fO37mtHQVmBq65B2A3DK0wxzV7ND5Wrr\niAnkSV/lUlnGTMWOr/ff96Tx5sIjD9/vfdhKNxibpq+LCc79pxa6MA1mXiSIBGME6J2tTeuTC/uG\nhunQ8JzrASt3Hz9Nlx56wvrzNvgODdNUoCHZg7S0PxDbQnVocHb2QRMuk1ABwxyFBOnWkEqzLUD3\nH/2dA/3R+xGiI+UI0AeuOi9omruYIC2Zw3A/WKdWm10BeuhNhTYpYTomOFOFXkOWm9xcQbprSJV5\nyI2EoUFxaDU6pS2ECyfWQcgYSxlftmq0FBekU8RUoHPh/WAvUkSEKVo42vAaE6a7crRolHjEt68a\n3UVFGTbGN3RHkJbCHscdIyg8Z6xCTzFbx1QXTqrQqMWd9xyw3pTnC9KS+ymGsXzh2VaFppUjP0J0\noNwBOrQa3RoaplOEhOchVeiYIJ0bVeg155h2rRt6UwN1VNW5kvmPh96cOHaQJkCjNqlBWtoffGNC\ndUzFuUQbRx/V6F0kiQClKtCxQVraG2xLBOqYqnOONo4pgjQBer1Y39AD5i/ONZezUcEe6NhqdK7Z\nPdpgW/ICSnhGbmOdU22ADXmSYe5WDF94pgpdBmnCo3QLRxtEUx4Fbgq8scE6tVUjZx/0mEGaAL0w\nbZAd+jCQlH0GGHJxn/IhLCXCNOEZNfCNqzaMuuZbjgnTQ4RWnX0Bmnsl0pEoHMbsgR4SprtK9C/3\nlbiRsA23pcI04Xm9ed/Qu8E2d6BOqDiPGRhLzjHd/3vEXFQJzahVyAfUmDDdNSRYx7ZqhFSfZzMH\nfaVIFhZTzQPdDahDA3VupWbh6MoZpgnOyxJcHfWF3n7IztiWkTs41vZIcIIx1kXo2AoJ012le5Zj\n2jbmdp9EjUgZBjU8SEXaH1rHDtVjhGYbVwDuB2zCMroGtxoU6GUuGS5dF3uqTEC6dvzEhOlWyuO1\nU6X0O/PekAfpo6eWAG3iC7UpIXvKoJyK0IwQ3eC67rNKmII0F0kgz3iMCdMtX7ANDdk5bwjkPSE/\n0khHzQE6xBwDMTCGJbQZdIM0F0sgv/64GtJKNdZsGaXfC5be0kGIXpl7gAYAwjMwHtd4G+teBcb8\ntAjRIkADAIB8CLfLsPj/ywRoYB5qmoUC9Tj/0CnODQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAcvn/3LI3Fjx7Ly4AAAAASUVORK5CYII=\n",
"text": [
"<matplotlib.figure.Figure at 0x10c13c510>"
]
}
],
"prompt_number": 11
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 63
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 63
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 63
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment