PS4_part2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fast PDE/IE course, Skoltech, Spring 2015"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2 (80 pts). High-Frequency Magnetic Resonance Imaging\n",
    "\n",
    "\n",
    "Magnetic resonance imaging (MRI) is a non-invasive medical diagnostic technique, in which the atomic nuclei in the body are excited by a radio-frequency magnetic field.\n",
    "\n",
    "In search of higher imaging resolutions, many research institutions are now developing MRI with high magnetic fields and resonant frequencies. For example, modern MRI systems built with $7$ Tesla superconducting magnets operate at $298.3$ MHz.\n",
    "\n",
    "Understanding the effects that high-frequency MRI have on the human body begins with a model of the underlying physics. To demonstrate these principles, several simplifying assumptions and approximations are made that reduce the governing equations to the Helmholtz equation in non-conductive media:\n",
    "$$\n",
    "\\Delta u({\\bf r}) + k^2({\\bf r}) = v({\\bf r}) \n",
    "$$\n",
    "where $u({\\bf r})$ is the electric ﬁeld strength in volts per meter, $v({\\bf r})$ is the electric excitation, $k({\\bf r})$ is the wave\n",
    "number, and ${\\bf r} = (x,y,z)$ is the vector deﬁning the location of the observation point. \n",
    "Note that $k({\\bf r})$ has the\n",
    "following relationship with the other constants of electromagnetics:\n",
    "$$\n",
    "k^2({\\bf r}) = \\omega^2 \\mu \\epsilon({\\bf r})\n",
    "$$\n",
    "where $\\omega = 2\\pi f$ is the angular frequency in radians per second, $\\mu$ is the magnetic permeability in henries per\n",
    "meter, and $\\epsilon({\\bf r})$ is the electric permittivity in farads per meter, as a function of the location vector ${\\bf r}$.\n",
    "\n",
    "Let us simulate the scattering of EM waves by the presence of a human head. Human tissue is slightly capacitive, but non-magnetic. We can model it as a material with magnetic permeability in the freespace $\\mu = \\mu_0 = 4\\pi \\cdot 10^{-7}$ H/m, and a location-dependent electric permittivity $\\epsilon({\\bf r}) = \\epsilon_r \\epsilon_0$. $u({\\bf r})$ can be approximately chosen as zero on boundaries. \n",
    "\n",
    "Here is the $\\epsilon_r$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x11edd4fd0>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAXcAAADvCAYAAAAJgpRkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4FFXbwOHfbEs2vZECCYQWIISu0osgCKKAiijKq58d\n",
       "O9YX2yt2VLCiYBdFEaUoCohI771GCKGFBBIC6T1b5vvjbCTGBEJIgiTPfV2b3cxO250zzz5z5swZ\n",
       "EEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIeogtwqGu1cw3AS0qqF1EaLOMVzoFfgX\n",
       "GwE4qVpAyT3L+77AfVWY78PAn8A3Z5lfJLC7CvOvLVcD3hW8Fw5cUc7wfqjtUR0cwHZgFzAX8DqP\n",
       "eZ1tW1fFBODxCt5bWwPLuxDLqE5V3Z/qNAnuFRsN/Op6Plf6Wd73B+6vwnzvQwW+/1TT/C6EMMAH\n",
       "OFXB+weAaMBaZngrIL6a1iEf6AS0B7KBe89jXmfb1tU9z541sLyaWIbmetQ0DQjk4in/tUaCe/m8\n",
       "gK7Ag8CNFYwTCewDZqCy6R/5Z0ACeAyVRe8GHnENmwg0R2WPb1RymmlAM+A3YFyZ8cvOTweMwCfA\n",
       "HmAxp6s7xgAbXeNOo/bLwO3AvLOMs4B//qhWV9Ze1nrUd1eiou9nHrAF9X3eXYn5ljd+JLCX8rfL\n",
       "s0AcsJozHy2WHClEosrfl67pvgUGobLu/cClZ1kXgOdd81gNzOT00ULpZVS0vhV9vjhgOqrsfs7p\n",
       "8gvwKuros7SSz1HRflSZ5XzG38u/B6oM7XC9PwohXG5B7dgAq4DO5YwTiQo43V3/f87pnSPH9dwF\n",
       "dehvBTxRBbQj0ISKq03Km6aD673DQEA505SdXyRgQ2WmALNcn6kNMB8V+AE+4p9HAQDXAs8AV1Ww\n",
       "juei7LzeLfVeL9QP0zTge6B3qffeKfX6MuCSaly/ku1jBOZwOus70/fj73q2or7rkv9L5lVWeeNH\n",
       "Uv52Kdnm7qjqqnjUD/yZ1r1kXm1R2esWVBkEGMbff0DLrksAKvhvByyoZGZ/qWWWXUbZ9T3T53Og\n",
       "theocrnV9dqAOiorma5EJBXvR+eynNLl/3rUD1IJH+oh04VegX+p0ZwOLj+6/t9WzniJqMwPVObx\n",
       "MDC51Pu9UHW6Ba7/56IC2PwzLLu8afoAO88wTXmHv4dRAQPUDhYJ+KECyRbXcCuQUma6KOAhoP8Z\n",
       "llciGhhYwXvTgeBy5lX6hOlJVCBZBqwEikq9V7psdgGmVmH9KmJFBbZGwBFO/5APoOLv5xHUeRiA\n",
       "CKAlsOkMyyhv/FTK3y5BqO1c6HrMp3JVGoeBWNfrWOAP1+s9rvmWty7hrnXpDvwEFLsev5xhGWXX\n",
       "90yfL4HT30sCkIZKaEJR+1BGOcs4035UmeWU/a52AZNQicOvwJoKPludJsH9nwKAy4EYTldv6MCT\n",
       "5Yxbum5U459VBzp/L3gaZ6+jrco05SkdKB2oQKWhgu4zZ5huMJAM3AasQO1EFfnT9ajIreXMy1zq\n",
       "/ThURv4mKkMszaPU69JVR2XXLwSVmTUHPj7DupRWgKpzt6KqGoZzOtMt7/vphwr83VDBdzkVt+o5\n",
       "2/jlbRf45zavjNLzcqKCdMnrkn27onUpr5ydbRkl61vRPAHyykz/GaoqLgT4ooJllN2PSv4/l+WU\n",
       "Fo/avkOBV4ClwMtnGL9Okjr3fxoJfI3KUJoCjVHZS+9yxm2MKngAN/PPDGENKusoqWIZgarfzKXi\n",
       "1iKrK5jmTHLOML/SlqI+XwPX/wGuz1BaJurwezpnDuygMvdHKnj4VzAvR6npNVSTyLKBHU7/ULZG\n",
       "/QhUtH6jUBmrWzmf5WwKUFniq651WUb5348PKuMsdK1Pt3/M6e/OdfxVqO1cUi1zNdV3ora8ddFR\n",
       "dfPXoL43L1QgPJ95VmQe6gf5EtQPaXnK7kcl5b2yyylb/sNc03yLyuDLq1at8yRz/6ebUIdzpc1x\n",
       "DS8bZOOAB1AZSSynqw5KdsxtwFecPnz8lNPVK2tR9YQLgf+Wmuf2M0xT0Q6fVmZ+H5Uzro46OfYc\n",
       "8Dvqh92Gqm8+Wmq8WajgrKF+1PxQmW0Sf68Hh7Nn7mXntQvVUqVEY07XyZamcbrety+n65LLm2dJ\n",
       "hp+Lyg5LPssC4E7+We0Ef/9udqDqgke55l3e9/MbMNb1WeM4XYVQdl4lKhpfL2d8HbXNZ6G2cypn\n",
       "ru7RK3hd0XsVrcsWVPXPLuAEquxkVXIZ5/L5bKgfzYxy3itR0X5U2eWULv+LUD/2b3H6aEaaSYpz\n",
       "Esm/uy15dWkM/F81zu8J/nlSrawOnG7h8NBZxn3d9fwY556513eermcPYDOqbry6GVA/Xs0reD+S\n",
       "+rEf1Tqpljk/NdHG+d+mJ6dP1FWHT4EbzjLOANTRUkPg2FnGXYg6R+Lk70cg4uw+QQXercBs1FFM\n",
       "dYpG1X//ARw8w3j1YT8Sol7oTcVZdltON727kdPZpRBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQ\n",
       "QgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEII\n",
       "IYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGE\n",
       "EEIIIYQQQlTZYGAfEA/89wKvixDVRcq1qNeMwAEgEjADO4A2F3KFhKgGUq7FRcVQA/O8DLUTHAFs\n",
       "wPfA8BpYjhC1Scq1uKjURHBvBCSW+j/JNUyIi5mUa3FRqYngrtfAPIW40KRci4uKqQbmeQyIKPV/\n",
       "BCrLKT1I/3sSJER1ag3s06p5pmct101AT6jmhQpRWhMgASpVtqt7BwD1gxEHDACOA5uA0cDeUuPo\n",
       "sL8GFn027wMP14Nl1rflll1mFFR/2a5Uub4BuLGaF3w2sy7AMuvbcv8tn3WkeqpU2a6JzN0OPAgs\n",
       "RrUw+Jy/7wBCXIykXIuLSk0Ed4BFrocQdYmUa3HRqIkTqv9iXevJMuvbci/UZ/2ntvVkmfVtuRfj\n",
       "Z62JOvfKuEB17qJ+qJE698rQZ1+AhYr641zq3OtZ5i6EEPWDBHchhKiDJLgLIUQdJMFdCCHqIAnu\n",
       "QghRB0lwF0KIOkiCuxBC1EES3IUQog6S4C6EEHWQBHchhKiDJLgLIUQdJMFdCCHqIAnuQghRB0lw\n",
       "F0KIOkiCuxBC1EES3IUQog6S4C6EEHWQBHchhKiDJLgLIUQdZDrP6Y8A2YADsAGXAQHALKCJ6/1R\n",
       "QOZ5LkeI2nQEKdfiIne+mbsO9AM6oXYAgPHAEtRdipe6/hfiYiLlWlz0qqNapuyduIcB012vpwMj\n",
       "qmEZQtQ2KdfiolYdmfsfwBbgbtewEOCE6/UJ1/9CXEykXIuL3vnWufcEkoEGqEPWfWXe110PIS4m\n",
       "Uq7FRe98g3uy6/kkMA9VP3kCCAVSgDAgtfxJ3y/1uqvrIURVbHQ9qk2Vy/WsUq/bAjHVuVai3tkD\n",
       "xFZx2rL1iufCAzACOYAn8DvwInAFkAa8gTrp5Mc/Tz7psP88Fi3EmURB1cv2eZXr2VVcqBCVMVI9\n",
       "Vapsn0/mHoLKakrm8y1qR9gC/ADcyekmY0JcLKRcizrhfDL38yGZu6hB55W5nw/J3EWNOpfMXa5Q\n",
       "FUKIOkiCuxBC1EES3IUQog6S4C6EEHWQBHchhKiDJLgLIUQdJMFdCCHqIAnuQghRB0lwF0KIOkiC\n",
       "uxBC1EES3IUQog6S4C6EEHWQBHchhKiDJLgLIUQdJMFdCCHqIAnuQghRB0lwF0KIOkiCuxBC1EES\n",
       "3IUQog6S4C6EEHVQZYL7F8AJYHepYQHAEtRdrn8H/Eq99zQQD+wDBlXPagpR7aRcizqtMsH9S2Bw\n",
       "mWHjUTtBFLDU9T9ANHCj63kw8FEllyFEbZNyLeq0yhTQ1UBGmWHDgOmu19OBEa7Xw4GZgA04AhwA\n",
       "LjvvtRSi+km5FnVaVbOPENQhLa7nENfrhkBSqfGSgEZVXIYQtU3KtagzquPQUnc9zvS+EBcbKdfi\n",
       "omaq4nQngFAgBQgDUl3DjwERpcYLdw0rx/ulXnd1PeobG2BAC8pBdysE3Yy7ORMT+YBWZkwrRTY/\n",
       "QMfNnEVRXiCc8gPsrnGNtb72/x4bXY/zdt7lelap122BmOpYq4uMDjhwlUj30//YKog2ZjvgBLsJ\n",
       "NB3MhaDr4ERObOwBYqs4rXb2UQCIBH4B2rn+fxNIA95AnXTycz1HA9+h6iMbAX8ALfhnlqOrBgn1\n",
       "WQaMaQXhBjpdsZc7en/CHOcNdHLfSgNOko/H38ZOIpy1Bb3xMuTyecI9LGgwmOdveANsuZCfC1uO\n",
       "oTaDJJTqfGilynYk1VyuZ5/fil/0HECgEdo3BXMS+N0F5AG5qLZI6fwzD0kHioB+QEMo/AE2L4Um\n",
       "brCuSGWglQ1Udd1I9VSpr6MyI80E+gJBqMzmf8DPwA9AY9QJplFApmv8Z4A7UCnlI8DicuZZz4N7\n",
       "Lpf1O4nP3CI+t9/LtQfnE5yewrJVgygutDC242SsCwv/GtsJpA8K5psf7sYQ4uTmm77G4ihiZ6cO\n",
       "3Bv2ISv29ue78UNg4XHA+4J9qn+PSgX3GinX9Tm4l/zSDRgO/p0gfx84t0CGFWJ3Q8B1LVmy4ypw\n",
       "eJ2eSNNpF/kFXitSaBwO5pZgzIfwcCAY0g/DgiXg7pAsHqo/uNeEehjcNcABxiLcmoWjz3PnxugZ\n",
       "zJl3C4VH3Ikq2sdv6d3xzsohex7sP6WmKNlhWoWBz39ALwaDBziL4YPx43ll+0uMajaTsU0/4P7k\n",
       "z9jTyIrkOZXO3KtbvQ3uDsBigkHR4OEGTk/ITob1xW7o7/VjaNwMWGaERQcAT1TKYlcT39Ec2tpV\n",
       "Rq+baKXtZdOKAThW5eL3OGip8Od82H246vXIdYUE938lBwQEYWkJj699n58PXs/h1S0osHvyWten\n",
       "6DNqCo2OFeDjB6ZCyE4H3QRB/oARkk5AkQ4NvCApV224wxEG3H/uz02FP3BNt7kM039hdNcZsCUB\n",
       "sFzgz3shSXCvTTrqhETXoeCIg3w3OJ4AIQ9CF++tHPrKHfoFg9kIx91hkQFGaxDkgFgjrCyCdnYo\n",
       "AHbZoZGN8B65fDdmNL1e3IZ2OeRvgyXL//o5qLckuP8rmfAd50vU8DhCI1I4PjucvgF/8GjCm2TN\n",
       "zSR0tBFHHx9s7c1wEhrOP4l9qY6xNdj2glsGYAX6QPIUsHiDf2s4uhnyWmrMvuZZPjA9wuhrZjDl\n",
       "7dHwxSnqb4CX4F5bnKjA3r4TeJgh1xvmO9qwwXssH2/qzdgbVhLdZR/3HvoYMiCvpwfDiufy1ICJ\n",
       "rP7xCo6MCWVHVhemBt7LdPOtfN7sAZ6Pe464T6P54ZPL6DIsjTs9XmHEwl8JaAQ/LgJrPY7w5xLc\n",
       "pRqrVuhgCiUr05+ofntYnduDFqH7efK1Z2iwMpO2XcBxtw/7+zaj0ZupNFqbinaljrkzGMyw/w9U\n",
       "0ws3IBvC+kNgA3AkQhNPaJSp0/b1V3jsmons3tmeKe8/i6mxE7XrCVGz/IGsAtCDISkOvur3Ct/4\n",
       "/Icxh9cz4pE53DvrY9aO7cLtb0zFIyGfGZffwpW2FWQP0UjMCOWNxk/wpueTjHN7lweOT+aNLc8x\n",
       "r+sIOv2Swdbf2nD/nlkU92+AWwSM6aPamImzk+BeGyxOGGLlypsXMH/ujWTuDmHKnvtJOwFuPYEB\n",
       "EPxaBt3XbFdtMd6BhBahqo3GPohph2qgFwX8Cfiggv1xcHiCbxPomg+PT5zMyuX92evZmgaeaRfs\n",
       "44r6wQl0RB0fhjWF/Rvggz/eZsl7V9P+vS18vuEBCgPdmTDzaS6ZvZv9GW14+rYXCfs5jThnc+bY\n",
       "b+V/zV7j8rRVPOt4lR9NI3nR/gKjmnyH7WMPRrT4Ee70hm1WZh69g43fA50v7Ge+mEhwrw2+OQyY\n",
       "9BsHT7TE8aWRuUeHkjsjm7ZjgTWoFtNeYJzrhGQgCposToED4AgA+2hUU7F0oBuQA3iCuQeYIkDr\n",
       "CBEvQu4SGNP1K1bNH8DASfEq7ReihjiAlm3AtyfE/wE5r3RgZp+bee/7R3nU8jZHAppw5aE/eLjw\n",
       "Qyx9i5nk9xhhRanQHgIL0/nPwek0jk1mltsouu3awXUZ8/nZMoy3Cp7hypd/ISZrL472BibtfIhf\n",
       "c4fScTDYUiCmpdodxJlJnXtt8M+i3/J0IlKSyA/15v9GjOLqywAfKJgP1s9QDeuCUGen3IFiVFq0\n",
       "H+jlGnYQ1YehEXVmSQdWANlAEcQnwNuP3860X77g9TWP8nSn/8LBdKC+BXmpc69pJZfvXuIBTe+B\n",
       "44n+NBqSRHTkLjZd2pfwUwms+mMQkd3iSWnZgJZ7EqEQVa4zUUefp1CJigPogorYrvK/tG1P7jF+\n",
       "TG66D2lrwnDcF0v6+x0xPwn5UbBiZf28bO9c6tzre8uiWqDTar+ZlbMGYDykk5vtjnEc7LmnOTHz\n",
       "DmINAbahLpexAV6o41wjkA/Z/+eJVc/DvBwoAv3/NDY1bE9o4kmafHocevDXpXwt0yCGNGgNoaQw\n",
       "btpM3h3YBXWxpRDVxwm0bg5hRtj3IXzd6W6atU9i28y+mA/aOdQmGu+OORiKHET+mayCtgGVpASg\n",
       "eudJhdjbWrDE5woGFCyj3af7IQHsYw30tG9i094+ONON4AlPjnsdazAYvcBrLPQ5BCsS61/aci6k\n",
       "WqZGGYB8EhNaEHhDCuN9JrLrewemZtBsS6KqP/RBFfySNKQDZFzlhX4VkAU+Y/MwXwV7R8Ou22BD\n",
       "S51LZu9iX4sWcBUqG9oE9lWQ9h00b5AGI2D6/juY02YYPNcZlRoJUX0cgDMN3CKgVVd43XoXLVvs\n",
       "xe32YgxeTvxTMjA57RjQMet2dXQZB+wAbJDYP5Rru/5Aj8e38Kj1Sdo3WEhIeDxFMSZM05y4x9r4\n",
       "6pKbeL7oZQaeWMRWxyWk5YViNoJ9MoQ2BndJTc9IgnuNysF3SASOwwZM/sWQC9EtQW8C7rE22Mnp\n",
       "LeAGhUPM7BrckswgPxxPo6phLoW8z6y06gseFujQEIxZOl2Pb1bXUPoC4eqqvsAhkBgVBZOcHNrv\n",
       "S8SGJDqOWId0SSCqk47qLjPaB7R80AzQ8AoPXk19HraiGgWUboV7ir+6/9ka1ZFNfdsRcTCFeQmj\n",
       "mP3ICMIedsetBYzs8wuGJs6/OuS59/AXHG0cSlC3NL6//iaC30wh8TAcOQ7EQLFb7X7ui40E9xqV\n",
       "yXUj51J0xErK8khoDZ5PgbYDCqItsBl1fNsB0CH+skiSDWH4f5RE4q9ABCQ9FsSmjh0wvAgtJoDH\n",
       "TUAu+H1doOrjHUBz0K4DYqClfT+0g6M3N6JDzk66sxZpEimqk46K347ekL8bnOnwxrjH6bJzt+pq\n",
       "zRPV+P0g6txQGpAE2y5tR/frlxF3pJX6EVgNA39bweEhzUj/uj1TUh7H7OX8qwmOV1ohCwtGsCTv\n",
       "StqyD1MTOGJzBfcOYKuvl3FUkgT3GqVjxEHHoZtxt+QxcvknOG0GnM0MFPhZyRlvUVfl/ake0QPj\n",
       "6fPWBtwSodF4yL3Ogm43cPn2DeQ1dVf18kZU3Xw+qt27ETgBhYvhg+vvY5rlPkgywF2NyHDzJ5lg\n",
       "JHMXNcEaDkdN8MIrL7Nkai91sj8amIdqAdYSDrpFMrLxLG4ZNYNxXd6j/7QZ9Dm+XpVfK+AOboeK\n",
       "+D2gP9px/XSfG64OxpJCGjEl+h6Wxvej2A0u8YJW7up9u0NK9plIrVUNMzicvBDxP6aEPETRkkQK\n",
       "bnLHYTcSkJuJ4WtUZnMA9COQlQBHXy6m4RdB7BzZCjN2NHQsnsWEvJ2JswsY1gOBQBtUf6CtYOaj\n",
       "o9C3m3jxpZdIIwhaAk7IaeJD9wY7+dX7Kuw5RdTP9gWiJtgBosEvxMSct6J4+PXfVdVLNqrlS2NI\n",
       "8/djgv+zHJ7SlF2zm9Lkjp0cmPkwNENV27jDylE9cITrfJA1lvnBQxkf/Ba+vlm843iY27d8hWNm\n",
       "Lp7dk+CGFaTqoGmoqJULg1rAmm3Sk1JFJLjXAm9TDkuThnDZO+BsVQRWMMxBNQPLAXs8mK6EgEII\n",
       "cELCSANWCrBjotmpBEJ6Z7IvFVr3BcPLDm7vNpUIcyK7He2Y99lN6B8bVRupFqhD4gzAG17v8QQ6\n",
       "Gu//9wlSnzMirYNFtXKHUF87sS/chJaoq4jvAYSDI9VA30WruX38VO58+lOeevotXkt5Gda7xvGG\n",
       "9x6+l3Ft34HxZjJv9MLn1WK0aB2K4fU9Ezi6TMehw+2+C9nvVK0os3RoPhBoAWEtoO0DsKdI0pby\n",
       "SLVMDXNz2HAaDTBJ4+2McRje1DF8rqt6yHxgKZi6otIPdzWN1z51gbUDEwEHMtmdCumftuXxb19H\n",
       "vx++COzNi52fZe5Tt6CvMNK7xzKm9bwNd+98OAnkQushu/k/21esowdDn50Lff053XutEOfH6Prj\n",
       "cIBjps644W+o1l+NAB0e7voWMfdspSM7MKDzFk9iCir6q+cvPRTe2Pg0tDvFjLZjMH9ZxLfzdX6b\n",
       "CHvehT2/66TaodmDoHVSaUkeqqsDazqwHLBBsy6q9bBUz/yTZO41ypelqy/h7ts+wpJWyLZWnVT6\n",
       "kYr6WdVQPYaXXJCkA7mgu1ouulHI+m5dCNTWcV/4N+y6Iwb2JADusD0bYrzRbnCw3t6L1ff3V1m7\n",
       "ATDC3MY3EpFzhE0BHdm0oQ+eEyDvmdawPhk5kBXnQ0N1gB8WCx7eYJsD73b6r7rZRhGkBftync8c\n",
       "PMmj0JWx7CcKEk30c2xQ89gMS/YPZcmnvenVbRYLs1UVfC6Q6oS2BjB7QcpUOFasAngBEBYGlsdQ\n",
       "R73x4O4F/ULgt5PgcEq2Wpp8FzXKmz+/82cEcxk270e+nXErjf0TyE7zUlm7kdNN0B3AXMj5HsjS\n",
       "AI1cvOn35Tri0UjeFwaWHKAALCdBy4Y40NcYsZ+0cOOHX+HbOhOGO9B8nWQGWkkMCKap4zBjO79H\n",
       "3iZvgqYV4BbRGDTpeklUnYY6X3rkEzB3B80b8p8F5kFRthtrG16GN9kU4YYbRfyacw33dPiGWWtv\n",
       "UtdlAPmDzMQsXc+VSYsoLFLtuWyaeuQASTqkZYOXL3S5CkL8IdjoqnPPRv0KtALdA45lwHU3gVWT\n",
       "KzpKk+Be40wkvdGCLQVdIB0S5zTm2hd+YOU0SJzK6SzeCAwCr1D+arkYxnHSl6tzp8UHLBAaADHN\n",
       "CYlvDCPbQIxO6+d3wh74YcJ/yIr1g9lGmj0XS5AllUKs9NmwgK4d3oU3k0l7uwE37poOLSRzF+fH\n",
       "BBSdALbB0Ux1z9NNn0G3aS8zvNN07lgwA4OrIA/2XkT4pEScDV3lLh3mXXYdHIwH3DlaBO2B4U2h\n",
       "V6i6QW3bcPCKgFNpsH85pGVBlhncbkId/WbAwnHw03zIs8Gfv8CAziq7F4oE91oQHHGIvKP+dLhy\n",
       "Oxhg3Yre6C06ciQDsmeCfTsqo7GAdimkNVP3Qm06NZm4b2DDI3dhjcnn/nveR/vYyInrg+CQE9Zo\n",
       "JG5oAvmgJxhV87EM4BcL+2xt+M5xM8fWt8fSowsB3XQCM09xYHNryJTNLs6PBtiC4EQyhHcAuw38\n",
       "8mHEwRyINxB7bwfe/PQZlmVfgSd5GP3sXOs7DzwhuWcwj33wFP5eweTPz6bjQGh9OxxJggPJKik/\n",
       "eQLIBk8D+NjBUwOnA+y7gTTQV0GkB7SJBM0CDdtD2gnpjqA02ctrnEbSE+ng5aDXPUvp8fky7J+7\n",
       "cf2Tv9H+D1/25sOhn1F9bWhAdzjYMBIdDf3zXLK+GcI9Bz5h7MAPGZ47H/0hYMtBtD+zYATk/e6n\n",
       "Wh+UnJDtDUNvnEeiOZz79KlMHPcYj3z+LpE/HuFUxyDWDWurTroKcR4MwNEUCEwBqwN874FGBrhl\n",
       "+8t8sfYxTKGF/HrPUD5650FCnSloZhi77WOwwhvXPk5aTgcO3xNBxGfJNPgvFByDA8XqYlbNCLZi\n",
       "SM+CZDtk2SDDAdHNwNQelQSNhugroPUQaN8etJ2wMUld+iEU6RWyNnjaYFw0fbst5TfDVezzb02X\n",
       "xC2EFKay6mA0wauzydgGTa4GWsHS53sQQRKBXkeZ8ct9PKq/w8mUMH7U+3PfmP9CkxZc8/Fyfim6\n",
       "Ds1RjG90xl8dMxUWelK4zBPjCTur7+zJpT9tI6b5buJ/jcL57W4ocqfuNxyTXiFrgx0YYIWsImj+\n",
       "PqQO8sX3tiwSk2DPrKFcu30uvGSj9buHOekMZbNHZ5r+kUiMbQ9aw3ymvXoZ3f8HjmsN/BzjxAq0\n",
       "8Ydm90H8GxAaAZhU0dYASxbsyYGmvcD3Qw3tGOSP0jmQA/uK6kfWLndi+rfJM8Kneaxb14tJlkdp\n",
       "80ccDoOFvc2j+az5M2wf1Zc8T9j6LaR8BY2Lkmj2+XGynHB5zlp2BXckcHcGTqcBrLHMWDiWqScf\n",
       "ZEzRVxw41YaMmaFkfB9KxnehFPzuzdS+d/Ng1Ee8cGoCjlaQf8AT5xdHociTuh/YRW0xAcsKQAuC\n",
       "g9MgID0b868aLebDVRsXcHipG5dPWsy+cY1IuyWNFzwncDgmkkPrI4lod5Smd4BzFDjMGoVA93YQ\n",
       "uR70SGjxKHivVA+fFerZbQGEXgfFSaD9BJoOu9Jhfz0J7OeqMsH9C1TLp92lhk1AVSRsdz2GlHrv\n",
       "aSAe2AcMqpa1rAsKndh1C69OeoEDlzeDQ+C7JpuJcePpm76S6Dsg0gucuWBKgTV3diFotjsP9p9E\n",
       "w6gEvnh9NG7FNtodi+Z6fmbRqP58c/B2mh0/pK72K3kUw9iZn/Nu/DhCg1O4+Ypv+M+jn4GvtAQu\n",
       "Q8p1NQkohrS9YHhNRyvQ0RuDdjsEfWtiQf+buGbTH2gBJo4kRhKZf4SWE+N599DjhFjBNBncdzho\n",
       "4AMBb4FhFWjhwJOo/NSBamDgAMIgZAqEvgOGz3S0PJ3Idn81wBFlVCa4fwkMLjNMB94GOrkei1zD\n",
       "o4EbXc+DgY8quYw6zgDZcXR220ThzVYGJv+hriY1AE2h8C4TdIfA6eB/pwV9i0Y8LXFb6SBobwY9\n",
       "F2yhz6+p06i6AAAee0lEQVTrOWEJ4UvPO7i/5WSu/XSh6mOmvJTFBwiFBQnDucy+iev0uapqSJQm\n",
       "5boaaMC2TOjcD7YuAlaCZgdTMXgts2PdYuM7+/9x//qf2V7YGS0LvDunc6JJA45/BM4QcISA0Q/V\n",
       "nYYOxIBWXsTWwZSP6jbYB7CAW6g0f6xIZQroalQbjLLKq/cZDsxEhZ0jwAHgsqquXN3iy9YX/fDL\n",
       "PUnygYa8FPG0GpwH5mQHuc098Denc/fTH5F8fQAmHLitstHUcJC4n9vgmZyPP2lc6jORmM1xBKZm\n",
       "nLGGJbe9O1e1+4nHDkzh2p9+geNyqqkMKdfVwIA6P795O3QIgozHVT/vGFB9HBWD5/R8Pnjicb65\n",
       "9WboAJkvhLDTowspxVBcBI5QA6lHIXUq0M4144pqlW3gOAncDOTDifXq4ifxT+eTfTyE6pH8c1R/\n",
       "hQANUYe1JZJQFyQLAAxkPnQIVuq8uH4Cp270hiIwZuh4eOUzccCj9GcZPb/YjsVWDP+BApsVssFm\n",
       "NmPAicl0lGsWLjhjidbt0L/pEjxWF7E7IIZj7zZB8ptKk3J9jkxAQjpsTYbDKWD/AIrcNfRAwAFa\n",
       "C3WhU0GKJz2iVhD95jY6s41CwFms5uEETP2guPUZQpIBOAQJ0wF/0CfDyWw5hKpIVb+XqUBTVM/L\n",
       "ycDkM4wrlb2lOX0gvoBGbsfZoXf8q9uBIjcj986Zzh0zvoOJqKYIRdC9YBN6kMY7bR/Gmm0HjDis\n",
       "xjN20W43GsmYFswdzb8g9HAyPtas2vlsFz8p11VkBhKA1j3gxH6wW1yHlUVAc+B/sJK+ZD3ZgBuY\n",
       "zZON3iIPiP8CLNudaKgLpw2Os9x7IA28OwG9gW6yEc6kqn3LpJZ6/Rnwi+v1MVQ3/SXCXcPK8X6p\n",
       "111dj3rCqVHoYUFzRWjdAe7bHYTGJHLik3Ae/vJt3Mx5DLEspWXqIa5/5Hsu9d6Cw2qh+bOBNCg6\n",
       "pQ5bS65sLc0OeR28yJ3jxTcJY/ivNgmHyUjd3g1ct/k5f+ddrmeVet2Wv24qVC9ogNHgqlHRQXMC\n",
       "EaDHw7Er4KbBL9Oxbywxu/ayw9mFSy4BswWwQiN3sHZQLWAwo65CLVs1YwcSIHsVBLluFF/XO8fa\n",
       "A8RWcdqqfjdhqMwG4FpOtziYD3yHOinVCFXrtqn8WTxcxUVf7HQI1Yk278PpisxaPBT7mbDZrFAM\n",
       "77/9GKNmTmfj3R3JSQ/AKz8Lzd9GblAgc+54kgRnKIXLLYTZTmDYzN8DvAX8tmYx5sUvyU70JaV5\n",
       "IPl7LdTt4F42OZhS1Rmdd7m+sapLrgN0wNgWzAdcA5yoSw5iVRcB/JRKzB/TaDMAvpl4E16eUOQB\n",
       "FEB4ILj7oNoixVF+eyQN6A/+7UDLAraofmjqshj+niD8eA7TVia4zwT6AkFAIvACqv+3jqjteRi4\n",
       "1zXun8APrmc7cD91O6pUgRNaOOlp2EARbiow50NOcy/s4SYO3deMIdrPZGp+/G4ehHdALtM/uJfj\n",
       "A8KZ5rifS+/dysJeQ3C0zeJ420DC96b9vUMNDXCDt358jowoPwoauOOlJZJFMX+/sWW9J+W6uhnA\n",
       "1AVM8acH6U7QDkHUh9AiAQ6/D/jBda8sID8M1n8PV9wKEaPBsREcgw1YVjrVKevW/P1UkbtqAx8w\n",
       "E1XdcxxylsmVGxWpTHAfXc6wL84w/muuhyiXCVYY+bDvvQxgoQoRVjCZHHRybiLfzY3cI77sbtKe\n",
       "3kGr+Tr2djgCS+YO5dFXJpNxVyBPPPUa18QupnnnvXiHr8J3Vy5FvSxsNXfA0RSesUxkjHMGqZYG\n",
       "vDjxaRyJp5DA/g9SrquZ5oSir6HA7XQ3AFox6G5AA2CLOtlf0jA9bjZ4A1lPgMdtkL4CAjOcOJ9X\n",
       "7d2LO2k4dSPmbXZ0C5jmg7YfGAH2ZbBklnRefSZ1vcrq38kORQYLGq57RmbB7u6t+G3rMH4NHUxG\n",
       "bAD5273I/Z8XMS12E/BiNhmFATSyJKJH2VkzriWfjXuQtHh/to2MwdJT59HUSWxe24sH2kym76E5\n",
       "PP766+QlesOuYiTJFLXBAOglwdzFnqZhctMhAvL2QXYhFKeBpSF06oKKQIVgC4fZB+COXuDmOv9v\n",
       "+UEn/7id3A/A2gZyOkCBB6Tfqdqw5iEtZc5EgvuF4AC7zYSbs0jF3UL4zP1O2u59mm/9xlB8yA10\n",
       "KMDKQx7vcWnjzaQZAwjanMfUBg/DL2mM6zmZxY6r+E/QZxyLa8zml3vBw/Bxq1ux5yQBWa6H5Dai\n",
       "FjkAE9gMFsCO8yRwCNgNBa7qGqcNCECddnMCNnC2VndZ2vA19PUB+wwwXQ1HvlY36sjcB9vWqdoY\n",
       "zfWQ6pgzkx++WmeEkynYJ9t4PmICdrMR/OCrmfexv2Mz8ppZaHHLn2ADT/J47tRr+EwqoO2H8Vzb\n",
       "8AfMvQvAPYfjR0OhZTGLTg7lld7j4S64vO9v2HPSUQe7JbuAELVDA1ZsBWcs+C3OR9fA2AV4FPgW\n",
       "Qgao7nx1gCzIvh92PwAZE8E4GRqZ4aQTSIfsPMAMrW6BsADIjVCnlsyojFQC+9lJcL8gCiHHSczW\n",
       "fThtqs26MwPe9hrHzabvyDoVAG4qc78j5DN6Xr2Ozm22c+h4FOnXhnDjc8c4viGA60LnkX44hA5F\n",
       "e7hn0BTWTLoEVfSlGkZcGHlArifqIjsTYNFV+6OBUORUwVkD8AafO6FpEzAawZQGXaa607QhZMRB\n",
       "QBTYs8HYHUxNIDteOgc7V1Itc0E4CB97nE/33cf4xpM51DGSnrM3sL9JFHMfH439kAnCwIyNDXRD\n",
       "b1/Mz9k34HmkAPMGG9eM/oniZAszD93KzAHzKNYMvJLwIj8uDieDthf6w4l6zAQ0MYA2HdgItjtN\n",
       "GIvsZE2CvWnqRKsOqn7lEHiNBKJB7wpu1iJCZoD/A8ApMJ0CPEHrD/u3q+oZUXmSuV8QhXQ1b4LO\n",
       "OkTBtat/Yv/NTRlhm8fbrz8AVwKNwUo+G3f2YO39A4meHcvOfdEYTDqeDfMJeuQEZl87t6b+yFLj\n",
       "FVzvN5OstKYX+oOJek4DvNoCNkj7DvTn7BACvpOha081jskE5MHO32HXG5A3E7YEt8e4Sycs0jWj\n",
       "xqgrCn4F20Jp61UVEtwviFDmzLsZ00sOjEcdeGfkscFyGZet28UnO++HFdDkkXisFBIdGEvfQcsw\n",
       "F9iZcvQxnLqBLN2bQRGLyfM3kxNoISE3ktXX9sJ5RA7ExIXjABpYYPls0L0gYADkxkFhqBF2QdxG\n",
       "aOAO5lDACmEx0Kgl5K0DXdcgGYzRqPS+GMiHwkyYHSfBvSokGtQ6G1zakD6vrmHjM92xbzQx/erR\n",
       "tNL2kd/PHffYAnCD5Z592Mwl+IVn4heewaX6FgxLDNiOm3l66mRSNjdk0ncPcLX+K5NTnoKEkv4I\n",
       "hLgwNODSMbDhC3AWaBS9Bw2+0lX/Gt5gN4FlaDB7pnsT89xBgq9DZejFkE0a7ILj70LDjqj+3HeB\n",
       "fhiMZ+luRpRPgnut0yDQTHJaKBO+fY6Quak03JZMToInt5yYzpb9XSEc1hf0Itk7mGmx93NoThSj\n",
       "h3/DjI/upOAtKzlveqF/f5wVM/sSpiVz7PogOLEPdTNVIWqfDngCJEKLTmCbYMCa78B5LdjnwuEp\n",
       "UFAAfgez8JuaBQdg9yLIN0PL773x0nMhC+KK4fAB6JkGadth2SFJWapKvrdaowEFGDvA9fN+JH5j\n",
       "NFd5LSC3jQc5rTx4b+hYlq4fCDaN6LE7cTcVEkoKN7j/yGWNNrLNvRN0hCRzGLmrdxLQzs7d+mcc\n",
       "cLYg+KYsVFfjQtQ+V8emDO0N6dvB2ATwAy0JDOsg833IPQWRvmC43IkhTYfm4BEJvhr4J+WQ4fTH\n",
       "sV01rGnZDiiCwBan+1wW506Ce63Jwv+5UCbumMLi767hgR7v0nR5Eq2C9uPTOJvnX52E092I96g0\n",
       "Rjb/ngBrKp8676Z1WBxrhvXjl/BrIAWiTh2EOAsj+i9m8NEVJDkjSH3HD1feJEStKsnYb3gAHP6Q\n",
       "kg7BT4D75w6cN0DOvXAgTlWhN+gKliY2LMXF6Hug+QvQ+gE4ebk/bTYeZvtudYOl4McADzi1RzWt\n",
       "FFUj1TK1QgcicG9m55uJd5Pb1JNWl+6j//fLSJ0ZApqmrsqIhlu7f4UP2ayiH6s/H8iWg714q/94\n",
       "di/rSOrdQXCzHfwa4v7YOvLsnnx/3fVw/CiyKcWFYAfaWCB7GxyMhfbTIPM1SFgHBZkqe/QAml0O\n",
       "dEf1X7sG4vaBeRdETgXNx0nRBvizAPo3ATKBLbBht/pRkAuWqkYiQi3QvE4S9Wk2AbZC1q/tBbvh\n",
       "4XXTVCEOBMzgFl3ImIc+4+WfXmRBx6EEtviT2wd/ys55neiwdg/Lxvdg6tx74MCftHy/mGG+v+A/\n",
       "NRMWHEQ2o7gQioGmkRDZDXZ8r3oSiLtHBXQ76jomownCImH5e1exLrsHr094DsNQCHSDbB3ytxtp\n",
       "kJ7F4mehoRnCbwPnHNi5UF2RKoG96iQq1DgdzwHtaXrtH/w28FrVwWwxkI369osBd+j22HJesz3D\n",
       "Bzfcywu/v0nH0I1cFbKALwe/jblNPPu79GP2nVdjHW/l9ftvYPK2J+AP5GJUcUE4gEYG6NYT1i04\n",
       "3QtkAafvAGkHGvmC2/seDNizkmveWEh8PET5Q4NBYH7Chw2B7Wl59Rp0Ewy4S3URfOpXONsNmcTZ\n",
       "SZ17jfPgqjlz+e2R4ep2Y2GoS+38XA8zNL8jjgm8yOrsvuTgA+kaO17sxmtvvMzuuLa0GgWv5r8C\n",
       "JPLAy2/TPPUIS94bqu5UoEtuI2qfEwjrAIWbwJqp6srdgWCD6ibAG2hsBd9xBvZ3i8R03IHBCHn5\n",
       "sG82ZK+FXUGt6PvrBlJXg5c/aHfBsS2Q75S69uogwb1GZdBv6j5+eO82yDUw5r9fqHrHzhA+8ghk\n",
       "QtioRO7tPIVdtGO6+61cyWK2Xt0Oi18xzRrE09ztKB+Mf4iNj1vxviKS35ZdzYmQIDgO7ChCpf5C\n",
       "1B4H6u5AjRvCjgMQ3hiaDYMiA4T3AbMRjGYIvAqMwU46P/4n7gmFcBfE3K7eM3SAxkeOc+oxO8dy\n",
       "oEczcC4DmxNOou53KN3enR8J7jXGCWFtWOXoD9uAm2Hl2sthO5AHSZ9E4jssnYcGv00EiehoLHr5\n",
       "WsYde5djniHc9dwHzGp/I8eyG/Jwx/8DrREeL+WxvcWlPLHvA2gJPq8Wgr/9An9OUd84gWZBkPo7\n",
       "BJrArwlsm69uLPvnCnUvjuZ9wRAIxX1MFG2BnLlg3wCWGIi6DUyvmbDedIxFB6BnMDhagcECPhvB\n",
       "xwTdrKpaR1Sd1LnXFKuTzvP2YDtpYfeiLvAHJB5pqiokj0Hjqw8xavi35OGJHROX2rdg9iom9plL\n",
       "GNN+NhGWRI5f2phUgwG0YppuzuDEWxE4xhpJ+zIAgGce/x/P/TQC+9rwC/tZRb1hB64MAY8BsOc7\n",
       "CDRA7Gp14vO4q548aggQCujgvshOcRTEzwHTd+C/Emw6NPOxM2sjDL4UTMVwdC0094fEHPAPgohL\n",
       "wL5QeoI8H5K51xSLg3jf1mSkB/LJr2Nw9ytQe4ARiICmQw7gRyYvLHyD2Tk3EmeK4qXnniL0luNk\n",
       "B/sRe7AdPyWMZN0zfSE0nFkpN/HEPRNxX5NHcnw4WKEYC5pnwIX+pKIe0QF/f0hfD81ugYiuapjm\n",
       "eq9VGJgtsOOTVhwMbIq+DiyXQOf+0MgDwtwhzAZpX6j28X4D4NCSjjQxwIGP1InUYHfAoPIgOa9a\n",
       "dRLca4QTzIV0b7SGGS1uYV1KbwpPWKEp6qxTc4gJ3MWAzOWsG9yJn5++gSfnf8BA/XfG9Z6E1s6O\n",
       "33UZeJ/KgridND2SzSXP7qTF6gRYalQ/EPEwYfFEtKdaIVenitriBAiEgBvA/zrI2q5OphpRl2tY\n",
       "2kFmH18O25viNj+Z2HnADuBOCGwHpqZgbQe7M6GBN2y+vzsfNbiL/JOQ7SrG8cmgp0JUoDQGOx9n\n",
       "C+4RwHLUpQd7UDfGAnWTrCXAfuB3/n6V8NNAPLAPGFSdK3vxsMB3bfj9+WtYYB3KV7/fqyokkwGH\n",
       "zsibvsWPTJIdIWyiK7hD2s8NGJ70C29lPgWrjXh3OkX+eyfxvzOKgql+2F7VmND5edXixhNoAM5P\n",
       "jBTr7qh2ldJq5hxJ2T5HxUBMMOTEQsFPYNjhuicqqp69TXMwNgbHSSOdf47FklxIjgNO/ApMATIg\n",
       "cyck7wGDDp26Q0GhkUe//wBjA9Wc0gw4bFB4QJoKnK+zBXcb6iZZbYFuwANAG2A8ageIApa6/geI\n",
       "Bm50PQ8GPqrEMuogA1f1XQxBOgYPp6qozATs4HN7Om39dtONDfweOIjmhgN8/9YwYm7bQcJHrTBs\n",
       "c9LjoaUc6+CDI7eYlk8lkLKjIW7r7DTLieelF5/Es0euiuVm6NppFR0mm1BtDMQ5kLJ9jnTAIwQM\n",
       "noA/OE+p6phCIMgKbm0BD7BfA02mJxL8FFw2UrV9Tz8IjsaqimbVcYhsBG4twX/IGkI+jMNzBLTt\n",
       "r/IWM2BzQJdrIBipmqmqsxXOFNRBFajbH+5FdaE/DJjuGj4dGOF6PRyYidpxjgAHgMuqb3UvFjZ+\n",
       "2zGceU8MIXhnBiQBDnDrUsgtvaZjxMEGujHty0e5/tlFJKZH8nSXV3APKGCvKRrb3FSchzOwvuqD\n",
       "M9YACfBih+f4dvRQRn77Cy36xkE+4AsbJ/XB3a3wQn/gi5GU7SooTgBLa7C+C2nTVD/rFiC0Nar7\n",
       "3gIIfCid1b9C7mowRkFYD/B6ARx3G5j2C1gNEN4BVn0InsXw5cJHeeHV8Tj7a7j7QbAZdmdB/j5p\n",
       "Dnk+ziXziAQ6ARuBEOCEa/gJ1/8ADVGhrEQSaoepRwpgYBPMfnb8DJnkRLhDKyAFml5ygEBO4Uku\n",
       "PmSpZDsZnnzpfW55fTaRt8QR69uGTe8NAZwUmCKId7Zhz+cteST3A2y3hNLJtoP4E60x3VAARdDo\n",
       "3sNsCusBoW1RLZBFFUQiZfuMdNQX0CIYtDAgFaxWNdwMeEQDBWAfC4WuVHv/Qtj/BuSshZzBnti/\n",
       "chII9AkHbStcOsGCT2wDEr0b8sjGj0l+UyckBFJtEBMI7lnQoaHUu1dVZYO7FzAHeAR1XWRpJT1+\n",
       "VqSebZuTDHvyFzy/yeWyTTtZk9kL1gAFoKGjAcGk0okdTL/vFjxHqMv7DCl2lm8ayLepY2BVHHRu\n",
       "hKFLEf8LfZ62DxzgycGv0mvYGooWuJM/1RP7SiugMyBsMROuGw/XGJATq1UiZbsSnECbIMjrAKZh\n",
       "QCwcy4csXNm1E2gHpq3gdQN0vwLCfOCkAwL6GPA0FrF3AwRpYPUEbQQ89cJbGApMTHn5SVKHZ3A8\n",
       "Gw4eUrdXtfiAYSD4XAv2elX5VX0q087djCr83wA/uYadQLVkTUFdUJ/qGn4MdaKqRLhrWDneL/W6\n",
       "q+tRFzjw9cxCM+sYTjoIiEiHJoAXREQexUoB8UQxY/UdPNZ7IsuG9WbOJTfx9odPcMI9lMXrB4Mp\n",
       "kf4/rqM4x43HZn8I2fBh2CN8/+AWbrV9DW10OKWBXWNN6uX8MPlWQo5mkOZlw57rfqG/gAtgo+tx\n",
       "zqq9bM8q9bot6krOukKzgMEKuhsQoBFk1ElzQGgoqhWYP2S9BJ7DwTQYwq4Cx7ugO5y4fe3kxEkI\n",
       "94GT6XD4wW6kEMp9we8T+nwKd/zvIcLMoHlAbhbk50HxYth7EJoHw+GU+tnmfQ/qjH9VnC24a8Dn\n",
       "wJ/Au6WGzwduA95wPf9Uavh3wNuoQ9aWwKbyZ/1w+YMver7M/q0nHsYCtK3QNWAbPxy9Ba6CyyOX\n",
       "4E02K+nL4reHsix+IJaIYoJ7HcXha6bnkvVMn3Aj0+8bzYLfb8ZktrF8UE8uX78WcxC0bLwPcsHn\n",
       "3jQKi60Uf+jJoY9acdfE98ly9+OnrSPgklPUv+qZssnBlMpMVCNl+8ZzW/GLhgE4mgItV4BjtAG9\n",
       "s0a+00GhEbwHon4G98KBFLC+A018wCMcgoywezm0HwQDbgc9Fgp2QvdFW7DFmPEgHzsmOr0Klilg\n",
       "7AVhe2DnXmjlB13Gg3mDavO+L6X+tQmL4e8Jwo/nMO3ZDnh6AmOAy1EXzm9HtRSYCAxENRfr7/of\n",
       "1I7yg+t5EXA/9ejQVfGmYKaJW5//BFNXB1e0Wggtwdo8gyTCsWDjSn6j4Z1Hsf1qIe9zLw4Pj0ZP\n",
       "NJC3x4s1+X1J+qYV+gqN4m/dWB7UBzqDoye0Me2DYLCdcMekO9Spv0TQfrawKr0v/4t6hXrWgON8\n",
       "SNk+Bxpw2AkFQWD00zEanYQOVcGbIsAXaAtd+qj+Yf7MhM17YN9hCDWA4zAU7IfDCbC/CE4ts9OI\n",
       "YzTiGIGkcWJYmGpVVgBukdDWD5KPgNZZg4PgHiPdEZyrC3UyWlf7Tl2loXIMJ6ev3SuJA0bXaycq\n",
       "EJeNDyXX+5We7kybqfR4Bupf1l6eKLgwZVuffQEWWpvKNkss+ZL1Mv+fTUlpLT2f8po8li790nIG\n",
       "RqqnSn0V0rdMjdD5+8nN0tvCVmpYRYFYr+B12XmVZT/L+0Kcn4qOC8+11JWkJGebh5TmqpPgXmMq\n",
       "KpZV2Q1qYlwhLiwprTVLKmiFEKIOkuD+/+3dT2gcZRjH8W/R9qARpIdWrYX0oGAV1Es9VNGTNhf/\n",
       "nKLgRUEExZZ6sNqLOYkU/IMevFSxKlZEsdiDoIKKCkYUY21rtSkNWK2NB6FKDyrGwzvLTtbZ7G4y\n",
       "7zvbN98PLDvZdPLMm/nl7c7O7LOSlCEnd0nKkJO7JGXIyV2SMuTkLkkZcnKXpAw5uUtShpzcJSlD\n",
       "Tu6SlCEnd0nKkJO7JGXIyV2SMuTkLkkZcnKXpAw5uUtShpzcJSlDTu6SlCEnd0nKUK/JfT3wEXAI\n",
       "OAhsLR6fAE4A3xS3sdI6jwFHgSPAzTVuq1Qns62s9fqA7L+B7cAUMAJ8DXxA+ODyp4tb2UZgvLhf\n",
       "B3wIXA78W98mS7Uw28par2fuvxLCD/An8D0h2FD94eW3AXsJfzgzwDSwaclbKdXPbCtrg7zmPgpc\n",
       "C3xRfP0Q8C3wInBh8dglhEPalhO0/2CkYTWK2VZm+p3cR4C3gG2EZzkvABuAa4CTwFMLrDu3lA2U\n",
       "IjPbylKv19wBVgJvA68B+4rHZkvf3w3sL5Z/Jpyoarm0eKzCc6Xl64pbbJOJ6jRdc7nVfRk4vZgV\n",
       "a8/288CaYvlK4KrFbNWADiaqs5zrNjXW/cCZRa7b65n7CsKh6WHg2dLjF5eW7wC+K5bfBe4EVhGe\n",
       "/VwGfFn9o7eWbqkmg8lEdZquudzqnmZ+nvoSJdtrCGddx0k3GRxKVGc5121qrGdo52l8wHV7PXPf\n",
       "DNwNHCBcFgawE7iLcNg6BxwH7i++dxh4s7j/B3gAD101nMy2stZrcv+M6mf37y2wzhPFTRpmZltZ\n",
       "q7rkK4Up4OqGait/nwA3NVD3Y+DGBupq+Wgq25IkSZIkSRpeWwiNl44COyLWmaF9JUTrkrXVhP4h\n",
       "PwLv03734VK8BJyifclcrzp1NJ+qqjlB/IZX3ZptxRzv2dLgK1WuIU22m8h1t7oTxN3XTeR6oboT\n",
       "DFe2+3IOoSfHKOENJFPAFZFqHSfsnLJdwCPF8g7gyRrq3EB463o5jN3qbCSMeSXhdzDN4touV9V8\n",
       "HHi44t/WVRPgIsJlghDe2fkDYf/FHG+3minG26+UuYY02W4i193qxt7XTeR6obq1jDd1P/dNhA2a\n",
       "ITRgeoPQkCmWzquBbgX2FMt7gNtrqPEp8HufdepqPlVVE+I3vOrWbCvmeM+GBl+pcw3xs91ErrvV\n",
       "hbj7uolcL1QXahhv6sl9HfBT6euYzZfmCG1ZvwLuKx5bSzjko7hfG6l2tzqxm0+lbHg1SniGNUm6\n",
       "8bZqDluDr5S5huay3VSuId2+HiV9rst1a8t26sk95Tv6NhN+WWPAg4TDvc5tSbE9verUtQ0pG16N\n",
       "EHqybAP+qPjZMcY7zA2+Ur9TdRiynSrXkG5fN5HrVt3as516cu9svrSe+f8T1elkcf8b8A7h8OUU\n",
       "4XUuCD1EZivWq0O3OgM0VhvYLO0A7qZ9uFZ3zVazrVdpN9uKPd5uDb5SjLcfKXMNzWW7iVxDmn3d\n",
       "RK7LdYc12307FzhGOARZRbwTT+cBFxTL5wOfE84s76J9JcOj1HNCFcJ4Ok88VdVpnRBpNZ86xuLf\n",
       "JdxZs9zwajvweoSaK4BXgGc6Ho853m41U4y3X6lyDWmzPUr6XFfVjb2vm8j1QnWHKdsDGSOcFZ4m\n",
       "XNYTwwbCL2GKcIlRq85qwmuVdV4KuRf4BfiL8LrrPT3q7CSM/QhwS0017yWE5ADhdbp9zH/NtY6a\n",
       "ANcTPlZuivZlWluIO96qmmOkGe8gUuQa0mW7iVxX1U2R7SZy3a3uMGZbkiRJkiRJkiRJkiRJkiRJ\n",
       "kiRJkv7vP0kQ9iYYaCk/AAAAAElFTkSuQmCC\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11d43dc50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.io\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "mat = scipy.io.loadmat('MRI_DATA.mat')\n",
    "e_r = mat['e_r']\n",
    "plt.subplot(121)\n",
    "plt.imshow(np.real(e_r))\n",
    "plt.subplot(122)\n",
    "plt.imshow(np.imag(e_r))\n",
    "plt.suptitle('A plot of the $\\epsilon_r = \\epsilon(r)/\\epsilon_0$. Real and imaginary parts')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An Integral Equation formulation\n",
    "* (10 pts) Show that the Helmholtz equation can be represented in the integral equation formulation\n",
    "$$\n",
    "u({\\bf r}) + \\int \\chi ({\\bf r'}) u({\\bf r'}) G({\\bf r}-{\\bf r'}) d{\\bf r'} = \n",
    "\\int v({\\bf r})G({\\bf r}-{\\bf r'}) d{\\bf r'},\n",
    "$$\n",
    "where \n",
    "$$\n",
    "G({\\bf r}-{\\bf r'}) = \\frac{e^{ik_0\\|{\\bf r}-{\\bf r'}\\|}}{4\\pi\\|{\\bf r}-{\\bf r'}\\|}\n",
    "$$\n",
    "is the Green’s function to the 3D Helmholtz operator and $\\chi({\\bf r}) = k^2({\\bf r}) - k_0^2$, where $k^2_0 = \\omega^2 \\mu_0 \\epsilon_0$. \n",
    "What kind of integral equation is this?\n",
    "**Note:** if this task seems to be too challaging for you, skip it until the whole problem set will be solved\n",
    "\n",
    "\n",
    "* (5 pts) Explain why the itegral formulation is more preferrable here as compared to the differential formulation (think about the condition number after discretization and about the influence of the domain size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discretization\n",
    "Let us consider system of linear equations obtained from a Galerkin discretization scheme of the integral formulation:\n",
    "$$\n",
    "A_h u_h = f_h\n",
    "$$\n",
    "where $u_h\\in \\mathbb{R}^{N\\times 1}$ is the vector of unknowns, $f_h\\in \\mathbb{R}^{N\\times 1}$ is the vector of excitations and $A_h\\in \\mathbb{R}^{N\\times N}$ is the dense coupling coeﬃcient matrix.\n",
    "To model the head, use the MRI data for $\\epsilon_r ({\\bf r}) = \\epsilon ({\\bf r})/\\epsilon_0$\n",
    "```MRI_DATA.mat```.\n",
    "\n",
    "Consider the grid created by importing ``` MRI_DATA.mat```, which contains $N = 256^2$ ﬂat, square and\n",
    "constant-value \"pixels\", with each side having a length $h = 1/256$ meters. Each pixel can be used as the\n",
    "basis and testing functions in a Galerkin discretization.\n",
    "\n",
    "* (5 pts) Write down an integral that can be used to evaluate $A_{ij}$, the element in the $i$-th row and $j$-th column of the coeﬃcient matrix $A$. Specify $A_{ij}$ where $\\chi = 0$ (out of head regions). \n",
    "\n",
    "Let $f({\\bf r})$ denote the excitation on the right-hand side:\n",
    "$$\n",
    "f({\\bf r}) = \\int v({\\bf r})G({\\bf r}-{\\bf r'}) d{\\bf r'}.\n",
    "$$\n",
    "Let us set $v({\\bf r})$ to an impulse (delta-function) located at $x = 0.6$ m and $y = 0.7$ m.\n",
    "\n",
    "* (5 pts) What is the analytical form of $f({\\bf r})$? Write down an integral that can be used to evaluate $f_i$, the element in the i-th row of the excitation matrix. \n",
    "\n",
    "It is worth appreciating the ease with which integral equations can handle impulses and singularities.\n",
    "Notice that the location where $f({\\bf r})$ becomes singular is outside of our region of interest (i.e. the head)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You are given two functions:\n",
    "\n",
    "1. ```DEMCEM_ST``` (Athanasios G. Polimeridis, Direct Evaluation Method in Computational ElectroMagnetics http://web.mit.edu/thanos_p/www/DEMCEM.html). It calculates\n",
    "$$\n",
    "\\text{DEMCEM_ST} = \\int_0^h \\int_0^h \\int_0^h \\int_0^h \n",
    "G\\left(k, \\sqrt{(x-x')^2 + (y - y')^2} \\ \\right) dx' dy' dx \\, dy.\n",
    "$$\n",
    "```\n",
    "k0 = 2*pi*f/299792458;\n",
    "st = DEMCEM_ST(k0,h);\n",
    "```\n",
    "\n",
    "2. ```nwspgr```. It implements sparse grid quadrature, and was authored by Florian Heiss, Viktor Winschel at http://sparse-grids.de/. \n",
    "It generates multidimensional quadrature points via the sparse grid method. As the grid is really sparse, the calculations are done very fast. \n",
    "\n",
    "The ```nwspgr``` is a MATLAB function. Fortunately, you need fixed weights and quadrature points for all integrals. Demo code how to generate weights and nodes is shown below:\n",
    "\n",
    "```\n",
    "% Quadrature for some 4-dimensional function \"func\"\n",
    "[nodes, weights] = nwspgr('KPU', 4, 3); \n",
    "% 4-dimensional quadrature nodes and weights of order 3\n",
    "```\n",
    "Given nodes and weights you can simply write in Python\n",
    "````\n",
    "integral = weights.dot(func(nodes[:, 1], nodes[:, 2],\n",
    "                            nodes[:, 3], nodes[:, 4]))\n",
    "```\n",
    "to calculate a 4-dimensional integral of a function ```func```.\n",
    "\n",
    "* (3 pt) Describe which matrix elements of $A$ can be computed via ```DEMCEM_ST```\n",
    "\n",
    "\n",
    "* (12 pts) Using the coeﬃcients generated by ```DEMCEM_ST``` and ```nwspgr```, produce code that will compute any\n",
    "element $A_{ij}$ in the coeﬃcient matrix $A$ for given $i$ and $j$. Measure the time it takes to ﬁll a single\n",
    "column of the coeﬃcient matrix $A$. Notice that $A$ has a translation-invariant structure; estimate\n",
    "the time it would take the ﬁll the entire matrix $A$ in case this structure goes unnoticed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Fast Integral Equation Solver\n",
    "Without leveraging the translation-invariant structure, forming the coeﬃcient matrix along requires the\n",
    "evaluation of $N^2$ 4-dimensional integrals using quadrature. Then, to invert the matrix directly with Gaussian elimination is of arithmetic complexity $\\mathcal{O}(N^3)$.\n",
    "Fortunately, the coupling matrix $A$ has a block-Toeplitz structure that can be exploited to\n",
    "solve the linear system with $\\mathcal{O}(N \\log N)$ complexity by embedding $A$ into a block-circulant matrix with circulant blocks.\n",
    "\n",
    "\n",
    "* (10 pts) Check that $A= h^2  I + T  diag(\\chi)$, where $T$  has block-Toeplitz with Toeplitz blocks structure. Describe the steps required to implement the matvec by $A$ using the FFT. Include the details of each operation in your description.\n",
    "\n",
    "\n",
    "* (15 pts) Implement a fast matrix-vector (matvec) product function that will quickly evaluate the product $Ax$ for\n",
    "a given $x$. \n",
    "** Note:** to verify that your function is indeed a convolution  you can apply it to a few columns of\n",
    "the identity matrix, and check that the result is translation-invariant within the machine\n",
    "precision.\n",
    "\n",
    "\n",
    "* (15 pts) Embed the forward matrix-vector product function from above into ```GMRES``` to \n",
    "solve the integral equation without explicitly forming the dense coupling matrix.  Solve the scattering problem using ```GMRES``` for $f = 21.3$ MHz and $f = 298.3$ MHz, corresponding\n",
    "to the older $0.5$ T permanent magnet MRI and the modern $7$ T superconductor magnet MRI\n",
    "respectively. Produce iteration-residual plots at both frequencies, and the scattered image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 160%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}