MT Script - 3D.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:5ff50687aa6b5eff8d9f077fee29dfb14e6eaa12563cbe8511b1dddc2ba76e23"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from SimPEG import *\n",
      "from scipy.constants import mu_0\n",
      "\n",
      "\n",
      "def omega(freq):\n",
      "    \"\"\"Change frequency to angular frequency, omega\"\"\"\n",
      "    return 2.*np.pi*freq"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M = Mesh.TensorMesh([[(100.,32)],[(100.,34)],[(100.,18)]], x0='CCC')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print M"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "  ---- 3-D TensorMesh ----  \n",
        "   x0: -1600.00\n",
        "   y0: -1700.00\n",
        "   z0: -900.00\n",
        "  nCx: 32\n",
        "  nCy: 34\n",
        "  nCz: 18\n",
        "   hx: 32*100.00\n",
        "   hy: 34*100.00\n",
        "   hz: 18*100.00\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Setup the model\n",
      "conds = [1e-2,1]\n",
      "sig = Utils.ModelBuilder.defineBlock(M.gridCC,[0,0,-300],[500,500,-100],conds)\n",
      "sig[M.gridCC[:,2]>0] = 1e-8\n",
      "sigBG = np.zeros(M.nC) + conds[0]\n",
      "sigBG[M.gridCC[:,2]>0] = 1e-8\n",
      "colorbar(M.plotImage(log10(sig)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "<matplotlib.colorbar.Colorbar instance at 0x7f55faf27290>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEKCAYAAAAPVd6lAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FeX1+PHPySrIJqhhF43Ysi8iqHxb0IKNWFG+btgW\nqeCKijuCqAS/v1rF1h0XqlixIKJVBBEVF6AuLCIg+1IIJCFBIYSlLIFwfn/MJNyEuTc3yV3jeb9e\n9+XcZ56Ze06CJ3OfeWZGVBVjjDHxKyHaARhjjKkeK+TGGBPnrJAbY0ycs0JujDFxzgq5McbEOSvk\nxhgT56yQG2NMnLNCbmKeiPxeRL4Tkb0isk1EPhKRnu66s0TkHRH5SUQKRWS5iNwtIvZv2/xs2D92\nE9NE5B7gaeD/AacCLYDxQH8RSQcWAluA9qraALgKOBuoG52IjYk8sSs7TawSkfpADvAnVf2Xx/p/\nAvVV9dKIB2dMDLEjchPLzgNOAN73s/43wLuRC8eY2GSF3MSyRsAOVT0aYH1eBOMxJiZZITexbCdw\ncoATlzuBphGMx5iYZIXcxLJvgUPAAD/rPwOuiFw4xsQmK+QmZqnqbuARYLyIXCYitUUkWUQuFpEn\ngDHA+SIyTkTSAETkTBF50z1RaszPghVyE9NU9SngHuAh4EdgKzAMeF9VN+GcEG0FrBKRQpyTn4uB\nvVEJ2JgosOmHxhgTYiKSATwDJAKvquoTHn2eAy4G9uNMsV1a1c+zI3JjjAkhEUkEXgAygLbAtSLS\nplyffsCZqtoauAl4qTqfaYXcGGNCqzuwUVWzVPUwMBW4rFyf/sAbAKq6EGhQcp6nKqyQG2NMaDUD\nsn3e57htFfVpXtUPtEJujDGhFeyJR6nidsdJquqGkSQidkbWGBM0VS1fJCulsjWn3Ofl4tzcrUQL\nnCNuAvRp7rZVSVwUckdmtAMAvgQuCMmeVMcAIDI2ajFUz/FxVD2n0MUQDiV5+ZOZmUlmZmbY44j1\nGGIlDpFq1fBS/y/Ifg8d3/Qd0FpEWgHbgGuAa8v1mQHcDkwVkXOBQlXdXsVQ46mQG2NM5CRXcTtV\nPSIitwOf4Ew/fE1V14jIze76V1T1IxHpJyIbgf8C11cnVivkxhjjoTrFUVVnA7PLtb1S7v3t1fiI\nMqyQV0qraAdAbMQAsRFHq2gHAEDv3r2jHUJMxACxE0co1Ip2AJVghbxSTo92AMRGDBAbccRCDLFR\nvGIhBoidOEKhqkMr0WCF3BhjPMRTcYynWI0xJmLsiNwYY+JcPBXHeIrVGGMixo7IjTEmzlkhN8aY\nOGfTD40xJs7FU3GMp1iNMSZibGjFGGPiXDwVx3iK1RhjIsaOyI0xJs7FU3GMp1iNMSZi4umI3B71\nZowxHmoF+aoMEWkoInNEZL2IfCoiDfz0GyUiq0RkhYhMEZHUQPu1Qm6MMR6Sg3xV0khgjqqeBXzu\nvi/DfbLQjUBXVe2A83CKgYF2aoXcGGM8JAX5qqT+wBvu8hvA5R599gCHgdoikgTUpoLneQYs5CIy\nUUS2i8gKn7YnRWSNiCwXkfdEpL7PulEiskFE1orIRT7tZ7tfETaIyLM+7aki8rbbvkBETgsUjzHG\nREpyUnCvSkrzeTbndiCtfAdVLQD+BmzFeeZnoap+FminFR2Rvw5klGv7FGinqp2A9cAoABFpi/OQ\n0bbuNi/KsaegvgQMVdXWOA8lLdnnUGCn2/408EQF8RhjTEQkJXm/vhV4Uo+9ynPHwFd4vPr79lNV\nBY7bg4ikA3fhPAKrKVBHRP4QMNZAK1X13+54jW/bHJ+3C4Er3OXLgLdU9TCQ5T5UtIeIbAHqquoi\nt98knK8TH+N8zSh5RPm/gBcCxWOMMZGSnOjdfmEiXOjz/vH/ll2vqn397dMd4Wisqvki0gT40aNb\nN+AbVd3pbvMecD4w2d9+qztGPgT4yF1uCuT4rMsBmnm057rtuP/NBufJ08BuEWlYzZiMMaba/B2R\nl39V0gxgsLs8GJju0WctcK6I1HJHNfoAqwPttMqFXERGA0WqOqWq+zDGmFiVnBrcq5IeB/qKyHqc\nA/vHAUSkqYjMAlDV5TgjF98BP7jbTQi00ypdECQifwL6Ab/xac4FWvi8b45zJJ7rLpdvL9mmJbDN\nPTtb3x3o9/Clz3IrYuXBu8aY6Jo7dy5z584N/Y7DcLmkW9/6eLRvAy7xeT8OGBfsfisdqnui8n6g\nl6oe9Fk1A5giIk/hDJm0BhapqorIHhHpASwCBgHP+WwzGFgAXIkzr9KPCyobqjHmZ6B379707t27\n9P3YsWNDs+M4uu49YKgi8hbQCzhZRLJxTkyOAlKAOe6klG9VdZiqrhaRaThjOUeAYe5ZWYBhwD9w\nLoT6SFU/dttfA94UkQ3ATiqY9G6MMRFTUwq5ql7r0TwxQP/HgMc82pcAHTzaDwFXVxymMcZEmJ9Z\nK7Eojv7mGGNMBMVRdYyjUI0xJoIqPyMlaqyQG2OMlziqjnEUqjHGRFAcVcc4CtUYYyLITnYaY0yc\ni6PqGEehGmNMBMVRdYyjUI0xJoLiqDrGUajGGBNBNv3QGGPiXBxVR3tmpzHGeEkM8lUJInKViKwS\nkWIR6RqgXwMRedd9rOZqETk30H6tkHvIz7+Xrl2bADBv3p8YOLB96bq2bU9h2rQrWbfudo4ceZgJ\nEy49bvtevU6juPiR417XX985YjmUV92cABIShAce6Mnatbdx4MBo8vPvZfz4fhGJ35/q5vX665d5\n/q6OHHmYRo1qRSwPE4PC8/TlFcAAYH4F/Z7FucFgG6AjsKaiUI2P9PSTqF07maVL80hOTqBbt6Z8\n9dXW0vW1aiWRlbWbDz5Yxz33nMexGzwer0uXV8jL21v6fs+eQ2GN3Z9Q5fSPf1xGjx7NGTFiDsuW\n5VO3biqtWjWIVBrHCUVew4fPZsSIY08vFBGmT7+GffuK2LnzQETyMDEqPPcjXwvOvzN/3Afa/0pV\nB7vbHAF2B9qvFfJyevZsycKFuajCOec0Y+fO/eTk7Cldv2RJHkuW5AEwdGiXgPvasWM/P/20P6zx\nBiMUOfXu3YqBA9vTsePLrF27o7R95UqvRw5GRijy2ru3iL17i0rft27dkB49mnPVVe+EN3gT+6J3\nQdDpwE8i8jrQCVgC3KmqfouJFXLXrl0PoKqkpiaRkCAUFIwgOTmR1NRECgpGoAqNGgX9wA4Avvrq\nemrXTmbjxgJeeWUJb775Q8UbhVAoc7riijZs2rSLiy5KZ+bMa0lJSeTbb7O57745ZYpnJITjd1Xi\n5pu7kZ+/j+nT14Y4ahN3/FTHudtgbp7/zURkDtDYY9WDqjozyE/uCtyuqotF5BlgJPBIJUP9+enY\n8SVEhAULhnLLLbNYtiyfqVOvYMqUlXzwQeX+p962bS/Dhs3iu++2cfSo0q9fayZMuJQzz2zImDFz\nw5OAh1DmlJ5+Ei1b1uePf+zA0KEzKCoq5s9/vpAvvriO9u1foqioOExZHC+UeflKSUlk8OBOvPLK\nEo4e9T9kZn4mTvBu7n2G8yoxdmnZ9arat5qfnAPkqOpi9/27OIXcLyvkruzsPXTocCrJyYnMnLmO\nOnVS6Ny5Mf37T2XHjsoNj2zYUMCGDccePbp0aT6JiQncd995jB07L2JFIpQ5JSQIqalJXHfd9NKh\nlWuueZe8vHvp1691RI9gQ5mXryuvbMtJJ53AhAlLQhitiVvhH1rxHChX1XwRyRaRs1R1Pc4zPlcF\n2pEVcmDlyltp2bI+SUkJJCcnsnv3yNLCtWnTcADatBlPbu7eCvbk38KFOZx4YgqnnFKb7dv/G6rQ\n/Qp1Tnl5+1DVMuPjO3bsZ8eO/bRsWT8sOXgJ5+/qllvO5pNP/sPWrQHPK5mfizBURxEZgPPM4pOB\nWSKyVFUvFpGmwN9VteQBzHcAk0UkBfgPcH2EQ40/GRmTSUlJZOLE/syevZFp01YxZkwvDh0q5vHH\nvwKcQlYdXbs2Yf/+w9U6YqyMUOc0f/4WrruuE2ed1Yj163cC0LBhLU4+uTZZWYVhycFLuH5Xbdqc\nTM+eLRkw4O1Qh2ziVXhmrbwPvO/Rvg24xOf9cuCcYPdrhRzIydlDQoLQsWMaN930IZs3F9KhQxqZ\nmXPZvLlskUpKSqBdu1MAqFs3lUaNatGpUxpFRcWsWeMcrd5117ls2VLI6tU/oQq//W06o0f/ihde\nWExxcWSGVUKd01tvrWT06F8xcWJ/hg//mMOHi3niiT5s2LCT2bM3RCSncORV4uabu7Ft215mzlwX\nsVxMjLPb2MafLl0ac+hQMevX76RevVTatTuF+fO3HNevWbO6fP/9zQCoKl27NmHAgDZkZRWSnv4c\nAImJwmOP/YYWLepx+PBRNmzYyfDhHzNx4tLj9hcvOR08eIQ+fd7kmWd+y9y5g9m//zBz52bRp8+b\nHD58NG7zAjjhhCQGDerI888vIsBlAebnJo6qowS6oCVWiIhCZrTDCCnVMQCIjI1yJKFTE3OCY3mZ\n+CAiqKr/K26C24fqo0H2fYRqf151xdHfHGOMiSC7+6ExxsS5OKqOAYdWRGQizpnUH1W1g9vWEHgb\nOA3IAq5W1UJ33ShgCFAMDFfVT932s4F/4Eyx/0hV73TbU4FJOFcx7QSuUdXjBjtFRONhCMgYE30h\nG1p5Nsi+d0Z/aKWiux++DmSUaxsJzFHVs4DP3feISFvgGqCtu82LcuzOMC8BQ1W1NdBaREr2ORTY\n6bY/DTxRzXyMMSY0wnAb23AJ+OVBVf8tIq3KNfcHernLbwBzcYr5ZcBbqnoYyBKRjUAPEdkC1FXV\nRe42k4DLgY/dfZWcSfoX8IK/WGrqCbSalFdNzAlqZl41MScI8YnpOBpaqUqoaaq63V3eDqS5y02B\nBT79coBmwGF3uUSu247732xwbtUoIrtFpKGqFmCMMdFUwwt5KVVVZ2qgMcbUMDEybBKMqhTy7SLS\n2L2xSxOg5IbUuUALn37NcY7Ec93l8u0l27QEtolIElDf/9H4lz7LrXBu2WuMMZtx5l1AZmYIjyv9\n3P0wFlXlUW8zgMHu8mBguk/7QBFJEZHTgdbAIlXNB/aISA/35Ocg4AOPfV2Jc/LUjwt8XlbEjTEl\nTqekNmRmZoZut+F51FtYBAxDRN7CObF5sohk49zY/HFgmogMxZ1+CKCqq0VkGrAaOAIM85kzOAxn\n+mEtnOmHH7vtrwFvisgGnOmHA0OXmjHGVENNGVpR1Wv9rOrjp/9jwGMe7UuADh7th3D/EBhjTEwJ\nz21snwR+BxTh3p5WVT3vmywiicB3OA+Z8H4iuqsqQyvGGFPzhWdo5VOgnap2AtYDowL0vRNnhKPC\ngX8r5MYY4yUMFwSp6hxVLbld6ELKTgQpJSLNgX7Aq/h5kpCvGBmqN8aYGBP+WStDgLf8rHsauB+o\nF8yOrJAbY4yXKp7sFJE5QGOPVQ+q6ky3z2igSFWneGz/O5z7Wy0Vkd7BfKYVcmOM8eKnOs5dAnO/\n97+ZqvYNtFsR+RPOsMlv/HQ5H+gvIv1wvhfUE5FJqnpdJUM1xpifOT/VsXcP51Vi7GvB79K9YeD9\nQC9VPejVR1UfBB50+/cC7gtUxMFOdhpjjLfwzFp5HqgDzBGRpSLyIoCINBWRWX62qXDWih2RG2OM\nlzBcEOTesturfRvOsx/Kt88D5lW0XyvkxhjjJY6qYxyFaowxEWTP7DTGmDgXR9UxjkI1xpgIiqPq\nGEehGmNMBMVRdYyjUI0xJnK0ptzG1hhjfq6K46g6xlGoxhgTOVbIjTEmzh1KTQmyZ1FY4wiGFXJj\njPFQnBg/g+RWyI0xxkNxHD200wq5McZ4OGKF3Bhj4ltxHJVHu42tMcZ4KCYxqFdliMj/ichyEVkm\nIp+LSAuPPi1E5EsRWSUiK0VkeEX7tUJujDEewlHIgXGq2klVOwPTgTEefQ4Dd6tqO+Bc4DYRaRNo\np/Hz3cEYYyLoEMFOPwyequ71eVsH2OHRJx/Id5f3icgaoCmwxt9+q3xELiKj3EP/FSIyRURSRaSh\niMwRkfUi8qmINCjXf4OIrBWRi3zaz3b3sUFEnq1qPMYYE0rFJAX1qiwR+bOIbAUGA49X0LcV0AVY\nGKhflQq5u/Mbga6q2gHnWRoDgZHAHFU9C/jcfY+ItAWuAdoCGcCLIiLu7l4ChrpPzmjtPtPOGGOi\nqqpDK+7B7AqP16UAqjpaVVsC/wCe9vf5IlIHeBe4U1X3BYq1qkMre3DGcWqLSDFQG9gGjAJ6uX3e\nAObiFPPLgLdU9TCQJSIbgR4isgWoq6qL3G0mAZcDH1cxLmOMCQl/49/fzf0v383d73c7Ve0b5EdM\nAT7yWiEiycC/gH+q6vSKdlSlQq6qBSLyN2ArcAD4RFXniEiaqm53u20H0tzlpsACn13kAM1w/hjk\n+LTnuu3GGBNV/uaRd+5dj86965W+nzD2uGFuv0SktapucN9eBiz16CPAa8BqVX0mmP1WdWglHbgL\naIVTpOuIyB99+6iqEsTTn40xJhaFaYz8L+4wyzKgN3AvgIg0FZFZbp+ewB+BC0RkqfsKOORc1aGV\nbsA3qrrTDeI94DwgX0Qaq2q+iDQBfnT75wK+8yWb4xyJ57rLvu253h/5pc9yK+D0KoZujKlZNgNZ\nAGRmhu7YMRyX6KvqlX7atwGXuMtfUcmD7KrOWlkLnCsitdyvAX2A1cBMnDOxuP8tGduZAQwUkRQR\nOR1oDSxyp9nsEZEe7n4G+WxTzgU+LyvixpgSp1NSGzIzM0O21yJSgnrFgqqOkS8XkUnAd8BR4Htg\nAlAXmCYiQ3H+RF7t9l8tItNwiv0RYJg79AIwDOfsbS3gI1W1E53GmKj7WdxrRVXHAePKNRfgHJ17\n9X8MeMyjfQnQoapxGGNMOMTTvVbiJ1JjjImgeLqNrd1rpZz8/Hvp2rUJAPPm/YmBA9uXrmvb9hSm\nTbuSdetu58iRh5kw4VLPfbRu3ZCPP/4D+/aN4scf7+PFFy+hVq3o/s2sbl5paSfyz38OYMWKWykq\neohPP/3jcX2iobp59e//C2bN+j3btt3Dvn2jWLHiVu64o3vE4vdS3Zw6d27Ml18OJi/vXg4cGE1W\n1p0899zF1KuXGrEcygvF/1cl0tJOJC/vXoqLH6FJkzphizlM91oJCyvkPtLTT6J27WSWLs0jOTmB\nbt2a8tVXW0vX16qVRFbWbh59dB7Ll2/n2DD/MSeemMznn19HUVEx5533Gldf/S4ZGem89lr/SKZS\nRijySk1NYufOA/ztb9/y2Web8OgScaHIq1ev0/j662wuv/xt2rV7kSef/Ia//OU33H//+ZFMpVQo\ncjp48AgTJy6lb983OfPM5xg6dAYXXXQGb7xxeSRTKRWKnEqIwOTJ/8vChTl++4TKERKDesUCG1rx\n0bNnSxYuzEUVzjmnGTt37icnZ0/p+iVL8liyJA+AoUO7eO7j97/vQKNGtfn9799j3z7nWX633fYR\nH374e0aN+pwtW3aHP5FyQpHX1q27ufNO5zx0r16n0axZ3fAHXoFQ5HXvvZ+WeT9p0nK6dm3C1Ve3\n48knvwlf8H6EIqe1a3ewdu2xi1Ryc/fy4ovfMWZML8/+4RaKnEo8/HAvDh48wtNPL+DSS38R1riL\niN43mMqyQg7s2vUAqkpqahIJCUJBwQiSkxNJTU2koGAEqtCoUfnzut569mzBN99klxZxgDlzNnH0\nqHL++S0iWshDmVcsCXdeJ510QpnfXySEM6fmzetx5ZVtmD17Q8WdQyjUOfXu3YobbuhCly6v0L79\nqWGM3BErwybBsEIOdOz4EiLCggVDueWWWSxbls/UqVcwZcpKPvhgbaX21aRJXfLzy97f5siRoxQU\nHKBJk8gexYYyr1gSzrx69TqNgQPbM2DA2yGKNjjhyOnrr4fQuXNjTjghiU8+2cjQoTNCHHVgoczp\n1FNP5M03B3Ddde+zc+eBMEVcVqwMmwTDxsiB7Ow91K+fSnJyIjNnrmPXrgN07tyYqVNXkp29h+zs\nPRXvxBVofC/SQplXLAlXXj16NOP9969hzJi5fPRRZI9ew5HT1Ve/Q5cur3DFFdM47bQGvP2250WF\nYRPKnCZP/l8mTVrOl19mlWk/dhPV0AvXbWzDITaiiKKVK2+lZcv6JCUlkJycyO7dI0lIEFJTk9i0\nyXnCUps248nN3VvBnhx5efto0aJembakpAQaNqxFXl5w+wiFUOcVK8KVV69epzFjxrU89ti/eeKJ\nr8MRul/hyqmk//r1O8nL28s33wzll788ucz4ebiEOqcLLzydXr1OKz0JXVLAs7Lu5NVXlzJs2KxA\nm1eJDa3EkYyMyaSkJDJxYn9mz97ItGmrGDOmF4cOFfP4418BTnEO1tdfZ/PssxnUqZNSOs7at+8Z\nJCQIX3+dHZYcvIQ6r/Ki9cUjHHn169eaadOu5KGHvuSZZxZUvEGIhft3BZCY6Hz5TkqKzJfwUOfU\nvv2LZd53796MiRMv46KL/smaNT+FNPYSVsjjSE7OHhIShI4d07jppg/ZvLmQDh3SyMycy+bNhWX6\nJiUl0K7dKQDUrZtKo0a16NQpjaKiYtascY5ypkxZwcMP/5opU/6X0aO/oFGj2owf34+pU1eydWvk\nTnSGOi+ATp2cuxI3bFiLunVT6NgxDRFYvnw7kRLqvK68si2TJ/8vjz32b6ZMWUFa2okAFBcrO3b4\nv+d0LOc0dGgXdu06yOrVP3Hw4BHatz+VJ57ow5Il21i58sfjPj8ecvL9dwjOmDnAunU72L79v2HJ\nwQp5nOnSpTGHDhWzfv1O6tVLpV27U5g/f8tx/Zo1q8v3398MOGPhXbs2YcCANmRlFZKe/hwA+/cf\npk+fSTz//MV8++1QDhw4wjvvrOaeez6JaE4Q2ryA0j4l/ZYuvRlVJSnp/8KfjI9Q5jVsWDcSE4VH\nHunFI48cm55XPvdwC2VOR44cZfToX5GefhJJSQlkZ+/hvffWRHw6Zaj//ZUX7vNRh+Jo+qHE0sk5\nf0REITPaYYSUqvPwbJGxUY4kdGpiTlAz86qJOYGTl4igqtU6CyoiOk7vCKrvCHm+0p8nIvcCTwIn\nq2qBx/oM4Bmcx2i+qqpPBNqfHZEbY4yHcA2tiEgLoC9w/NcTZ30i8ALODQhzgcUiMkNV1/jbp00/\nNMYYD2G8RP8pYESA9d2Bjaqa5T7neCrOY+H8ipsj8pKvgjVNTcyrJuYENTOvmphTqIRjjriIXAbk\nqOoPAebANwN8p7jlAD0C7TduCrkxxkSSv6GVrLlb2DLXc1QEABGZAzT2WDUaGAVc5Nvdo1+lT1zG\nTSEfG8YruKJhjHuSuSblVRNzgpqZV0lONfFkZ6j4K+Qtep9Bi95nlL6fP/arcjFoX6/tRKQ9znPp\nlrtH482BJSLSXVV954WWf8ZxC5yjcr/ippAbY0wkHQrx8zhVdSWQVvJeRDYDZ3vMWvkOaC0irYBt\nwDXAtYH2bSc7jTHGQwTutVI6hCIiTUVkFoCqHgFuBz7Bec7x24FmrIAdkRtjjKdwX9mpqmf4LG8D\nLvF5PxuYHey+rJAbY4wHu0TfGGPiXDzdj9wKuTHGeIiVe40HI34iNcaYCIqnoZUqz1oRkQYi8q6I\nrBGR1SLSQ0QaisgcEVkvIp+KSAOf/qNEZIOIrBWRi3zazxaRFe66Z6ubkDHGhEIRKUG9YkF1ph8+\nC3ykqm2AjsBaYCQwR1XPAj533yMibXHmQrYFMoAX5dj1qS8BQ1W1Nc7cyYxqxGSMMSERxnuthFyV\nCrmI1Ad+paoTwZn3qKq7gf7AG263N4DL3eXLgLdU9bCqZgEbgR4i0gSoq6qL3H6TfLYxxpio+Tk8\ns/N04CcReR3oBCwB7gLSVLXkcTHbOXYVU1PA9xlaOTg3hjlM2UtPc912Y4yJqngaI69qIU8CugK3\nq+piEXkGdxilhKqq80CI0PjSZ7kVzl8SY4yBzUAWAJmZoXtQTjwV8qqOkefg3Ipxsfv+XZzCni8i\njQHcYZOSG8GUvwlMc3cfue6yb3uu1wde4POyIm6MOeZ0SqpDZmZmyPZa48fIVTUfyBaRs9ymPsAq\nYCYw2G0bDEx3l2cAA0UkRUROB1oDi9z97HFnvAgwyGcbY4yJmp/DGDnAHcBkEUkB/gNcj/N8uWki\nMhTnu87VAKq6WkSm4dwA5ggwTI89LHQY8A+gFs4smI+rEZMxxoRErEwtDEaVC7mqLgfO8VjVx0//\nx4DHPNqXAB2qGocxxoRDrAybBMNuY2uMMR7CObQiIveKyFERaehn/SgRWeVeLDlFRFID7c8KuTHG\neCgmMahXZYlIC6Av4Pm8OPeBEjcCXVW1A86Q9cBA+4yNkXoTVmNeCLx+7O2RiSPUampeJjaEcfrh\nU8AI4AM/6/fgXGNTW0SKgdr4mc1Xwgq5McZ4CEchF5HLcKZu/yB+ngGrqgUi8jdgK3AA+ERVPwu0\nXyvkxhjj4RABh6X9EpE5QGOPVaOBUcBFvt09tk/HuVK+FbAbeEdE/qCqk/19phVyY4zx4O+IfP/c\nxeyf+53f7VS1r1e7iLTHuXppuXs03hxYIiLdVfVHn67dgG9Udae73XvA+YAVcmOMqQx/hTy197mk\n9j639H3B2JeD2p+qruTY/acQkc3A2apaUK7rWuBhEakFHMSZ0r2IAGzWijHGeIjAJfqlN4YRkaYi\nMgtKr9GZBHwH/OB2mRBoR3ZEbowxHsJ9+b2qnuGzvA24xOf9OGBcsPuyQm6MMR7i6e6HVsiNMcaD\nFXITU2rqhTE1NS8TGw4V/QxummWMMTVZ8ZH4KY/xE6kxxkRQ8REbWjHGmLhmhdwYY+LckcNWyI0x\nJq4dLY6f8hg/kRpjTCTZ0IoxxsS5g/FTHuMnUmOMiaQj0Q4geFbIjTHGixVyY4yJc3FUyO02tsYY\n4+VwkK9KEJFMEckRkaXuK8NPvwYi8q6IrBGR1SJyrle/ElbIq+je/HyadO0KwJ/mzaP9wGMPue40\neDCPFBcf92p1wQXRCjcogXICSKpVi9/85S8M37SJ0QcPcnd2Nr9+6KFohFopgfIa/OWXnr+rUXv3\nRivcoFQxLH+hAAAbJElEQVT0u+p+++0MW7WKUfv2cU9uLpe9/jq1TzklGqFWSn7+vXTt2gSAefP+\nxMCB7UvXJSYK999/PmvW3Mb+/Q+ybt3t3Hprt/AFUxzkq3IUeEpVu7ivj/30exb4SFXbAB2BNYF2\nWq2hFRFJxLn5eY6qXioiDYG3gdOALOBqVS10+44ChuCkPlxVP3Xbzwb+AZzgBn5ndWKKhJPS00mu\nXZu8pUtJSE6mabdubP3qqzJ9jhYX81TTpuDzgNWDu3ZFOtSgVZSTJCTw+1mzSKlThw9vuokd69ZR\nu1Ejap98chSjrlhFeb09YAAJycml7yUhgRsXL+Y/H/v7/yv6Ksqp/cCBXPS3v/HhLbew6bPPqN+i\nBZe8/DIDJk1i8sUXRzHywNLTT6J27WSWLs0jOTmBbt2a8tVXW0vXjx17ATfe2JUbb5zJ8uX5nH9+\nCyZMuJSiomJee21p6AMK39CK91OXS1aK1Ad+paqDAVT1CM6zO/2q7hj5ncBqoK77fiQwR1XHicgD\n7vuRItIWuAZoCzQDPhOR1qqqwEvAUFVdJCIfiUhGgL9SMaFlz57kLlwIqjQ75xz279zJnpyc4/rt\n37EjCtFVTUU5dbruOpp07cpz6ekc2LkTgD3Z2dEKN2gV5XWwsLBM/zP69KFes2Z893Jwj++Khopy\natajB9t/+IFlr78OOL+n7ydMoPfYsdEKOSg9e7Zk4cJcVOGcc5qxc+d+cnL2lK4fPLgTf/3rN8yY\nsQ6ALVt20717M0aP/lV4CvnB0O/SdYeIXIdzEHxvycGuj9OBn0TkdaATsAS4U1X3+9thlQu5iDQH\n+gF/Bu5xm/sDvdzlN4C5OMX8MuAtVT0MZInIRqCHiGwB6qpqyfPoJgGXAzFZyB/YtQtVJSk1FUlI\nYERBAYnJySSmpjKioABUGdeoEQAJiYncsXEjybVqsWPdOr7961/Z8NFHUc7geMHm1OaKK8hdtIjz\n7r6bjoMGUXz4MJs//5zPRo6MyW8alfld+Tr7llvI+/578r7/PgpRBxZsThtnz6bL0KGc9utfs2X+\nfE5MS6PtVVex/sMPo52Cp127HkBVSU1NIiFBKCgYQXJyIqmpiRQUjEAVGjUaR2pqIocOlR3LOHjw\nCKed1oDmzeuVKfohUcUjchGZAzT2WDUa58D1Uff9/wF/A4aW65cEdAVuV9XFIvIMTh19xN9nVueI\n/GngfqCeT1uaqm53l7dz7EGjTYEFPv1ycI7MD7vLJXLd9pj0UseOiAhDFyxg1i23kL9sGVdMncrK\nKVNY+8EHpf12rF3LB9dfT/7y5SSlptLu6qu5duZMZtxwQ+lRUqwINqeT0tNp0KoVWlzMtCuvJKVO\nHX779NMMnD6df/TqFeAToiPYvHzVadyYX1x6KR/ddluEow1OsDn959NP+eSuu/jjJ58gCQkkJCWx\n/sMPmXHDDVGM3r+OHV9CRFiwYCi33DKLZcvymTr1CqZMWckHH6wt7Td79kaGD+/O559vYtWqn+je\nvRlDhnRBVWnatG7kCvmKubByrt/NVLVvMLsXkVeBmR6rcnCGqxe779/FKeR+VamQi8jvgB9VdamI\n9Pbqo6oqIuq1riq+9FluhfPdI9L2ZGdzaocOJCYns27mTFLq1KFx585M7d+/zDBK7sKFzlffkveL\nFnFCw4b0fOCBmCvkweYkCc558XcHDuTQbme4bsaQIdy4eDFpnTqxffnyqMTvT7B5+eoyZAiHDxxg\nxZQpEY42OMHmdNall/Lbp5/mk7vvZsu//0295s3p++STXDZxIu8PGhTFDLxlZ++hQ4dTSU5OZObM\nddSpk0Lnzo3p338qO3YcG024886PefnlS1i27BZUldzcvbz66veMHNmTo0cXA+vIzAxZyfFfyNv0\ndl4lpgY/ZCUiTVQ1z307AFhRvo+q5otItoicparrgT7AqkD7reoR+flAfxHph3OSsp6IvAlsF5HG\nbiBNgB/d/rlAC5/tm+P81cl1l33bc70+MNrzPW5duZL6LVuSkJREYnIyI3fvRhISSEpNZfimTQCM\nb9OGvbme4ZO7cCEdfv/7SIZcocrktC8vj8Tk5NIiDvDT6tUANDjttJgq5FX6XYnQ9cYbWTF5Mof3\n+x2KjJrK5PSrBx/kh3/+s3Sc/6dVqyjat4/r58/ny0ceoXDz5mimUsbKlbfSsmV9kpISSE5OZPfu\nkSQkCKmpSWzaNByANm3Gk5u7l8LCgwwc+C8SE9/j1FNPJC9vnztrRdi0KR1oSmbmGMaG6lxAJacW\nBukJEemMM3tlM3AzgIg0Bf6uqiUPYL4DmCwiKcB/gOsD7bRKhVxVHwQedAPoBdynqoNEZBwwGHjC\n/e90d5MZwBQReQpn6KQ1sMg9at8jIj2ARcAg4LmqxBRukzMySExJof/EiWycPZtV06bRa8wYig8d\n4qvHHwdgX16e3+2bdO3K7q1b/a6PhsrktGX+fHqOGEFK3boUuVPzGv3iFwAUZmVFJX5/qvK7OjMj\ng/otW7LklVeiEXKFKpWTCFpcdixZjx51VwWcMBFxGRmTSUlJZOLE/syevZFp01YxZkwvDh0q5vHH\nnZk4eXn7ymxTXKylbdde255587IoKDgQ+uAqP7WwQqp6nZ/2bcAlPu+XA+cEu99QzSMv+T7zONBX\nRNYDF7rvUdXVwDScGS6zgWHujBWAYcCrwAZgY6zOWNmTk0NhVhZpHTuy9v33Kdy8mbQOHVj/4YcU\nbt5M4ebNpf+z9BozhjMzMjgpPZ1T2ral1yOP0GXIEBY89VSUsyirMjktfvFFDu/fz4BJkzilbVua\nnnMOl/7972TNncv2H36IciZlVSavEmfffDO5ixbFXC4lKpPT2vfeo/OQIXQcNIgGrVrR8n/+h4uf\nf5785cvZ5R69x4qcnD1kZRXSsWMa77+/ls2bC+nQIY0PP1zP5s2FbN5cyNGjTqk4++wmXHllW844\n4yTOPbc577xzFR07pjF8eJhKxpEgXzGg2pfoq+o8YJ67XIAznuPV7zHgMY/2JUCH6sYRCY27dKH4\n0CF2rl9Par16nNKuHVvmzz+uX2rduvQbP546jRtz+MABdqxZwztXXcXa6dM99hpdweb03+3beePC\nC/ntU09x4+LFHCgoYMOsWcx54IEoRF2xYPMCqNu0Ka379ePDm26KcJSVE2xOX48bh6ryP6NGUf+l\nlzhYWMjmL77g81GjohB1xbp0acyhQ8WsX7+TevVSadfuFObP33Jcv9TUJB555NekpzekqKiYefOy\nOP/8iaxe/VN4Agvf9MOQk2MHxrFLRDQz2kGE2Bj35z42xr7qVkdNzAlqZl4lOYnE9tzyylIdg4ig\nqtX6ZYmIMj7I2nhb9T+vuuymWcYY4yVGhk2CYYXcGGO8WCE3xpg4F57ph2ERN2Pk8RCnMSb6QjZG\n/ucga85oGyM3xpjYFEezVuKmkNekGQNQs2dC1KScoGbmVRNzgmN5hYSNkRtjTJyLozFyK+TGGOMl\nDJfoh4sVcmOM8WJDK8YYE+eskBtjTJyLozHyUN390BhjapZDQb4qSUTuEJE1IrJSRJ4I0C9RRJaK\niNdThMqwI3JjjPEShqEVEbkA59nGHVX1sIicEqB7+Yfb+2VH5MYY4+VwkK/KuRX4i/sgelTV8x68\nPg+3fxWocLK/FXJjjPFSHOSrcloDvxaRBSIyV0S6+elX8nD7o37Wl2FDK8YY48Xf0MqOubBzrt/N\nRGQO0Nhj1WicmnuSqp4rIufgPDntjHLbV/hw+/KskBtjjBd/hbxBb+dVYn3Zh3Ooal9/uxSRW4H3\n3H6LReSoiDRS1Z0+3bwebj/J3/M+wYZWjDHGW3jGyKfjPM8YETkLSClXxFHVB1W1haqeDgwEvghU\nxMEKuTHGeAvP9MOJwBkisgJ4C7gOQESaisgsP9tUeCcwG1oxxhgvYZh+6M5WGeTRvg24xKO99OH2\ngVghN8YYL3F0ZacVcmOM8WJ3PzTGmDgXRzfNqtLJThFpISJfisgq934Bw932hiIyR0TWi8inItLA\nZ5tRIrJBRNaKyEU+7WeLyAp33bPVT8kYY0LgSJCvGFDVWSuHgbtVtR1wLnCbiLQBRgJzVPUs4HP3\nPSLSFrgGaAtkAC+KlD5j6iVgqKq2BlqLSEaVszHGmFAJz/TDsKhSIVfVfFVd5i7vA9YAzXBuBvOG\n2+0N4HJ3+TLgLVU9rKpZwEagh4g0Aeqq6iK33ySfbYwxJnrCdPfDcKj2GLmItAK6AAuBNFXd7q7a\nDqS5y02BBT6b5eAU/sPucolct90YY6IrRoZNglGtQi4idYB/AXeq6l7xeSK3qqqIhOyR1l/6LLcC\nTg/Vjo0xcW0zkOUua2Zm6HYcI8MmwajylZ0ikoxTxN9U1elu83YRaeyubwL86LbnAi18Nm+OcySe\n6y77tud6fd4FPi8r4saYEqdzrDZkhrKQh+fuh2FR1VkrArwGrFbVZ3xWzQAGu8uDce4rUNI+UERS\nROR0nFs5LlLVfGCPiPRw9znIZxtjjImeOJq1UtWhlZ7AH4EfRGSp2zYKeByYJiJDcb7tXA2gqqtF\nZBrO0y6OAMNUtWTYZRjwD6AW8JGqflzFmIwxJnRipEgHo0qFXFW/wv/RfB8/2zwGPObRvgToUJU4\njDEmbOJojNyu7DTGGC81/YjcGGNM5YnIVOAX7tsGQKGqdinXpwXONTWn4tzCdoKqPhdov1bIjTEm\nQlR1YMmyiPwVKPToVnLl/DJ3ivcSEZmjqmv87dcKuTHGRJg7S+9qnFmTZbiz+fLd5X0isgbnokor\n5MYYUzlhPdv5K2C7qv4nUKdyV877ZYXcGGM8+TvbOd99eROROUBjj1UPqupMd/laYEqgT3eHVd7F\nuXJ+X6C+VsiNMcaTvyPy89xXibKzqlW1b6C9ikgSMADoGqBPyZXz//S5ct4vK+TGGOPpQLh23AdY\n4z6n8zgBrpz3q8r3WjHGmJotbDckvwZ4y7dBRJqKyCz3bcmV8xeIyFL3FfA5DXZEbowxnsJzRZCq\nXu/Rtg24xF0OdOW8JyvkxhjjKX6u0bdCbowxnuLnGn0r5MYY4yl+jsjtZGcV3JufT5OuzsyhP82b\nR/uBA8usb9a9O0O+/poH9+/nntxcLvzzn8Hn6UmxKlBep7Rty5XTpnH7unU8fOQIl06YEK0wKy1Q\nXp2vv57rvviC+378kZG7d3Pj4sW0v/baaIUatEA5pV90EUO++Yb7fvyRB/fv544NG7jg0UdJSIr9\n47aK/t8qcXKbNozat4+HiorCGM2BIF/RF/u/2RhzUno6ybVrk7d0KQnJyTTt1o2tX31Vur5e8+YM\nmjOH1e+8w4yhQ2l01ln0nzgREeHzBx+MYuSBVZRXUq1a7M7KYt0HH3DePfdw7Hbysa2ivFpdcAFr\n33+fOffdx4GCAn45YAADJk3i6JEjrH7nnShG7l9FOR3cvZsFTz/NjytXUrR3L026duV3EyaQUrcu\nn9x9dxQjD6yivEok1arFVdOmsfnzzzkzI+BkjmqyoZUaq2XPnuQuXAiqNDvnHPbv3MmenGPPj+52\n660cLCxkxg03ALBj7Vq+fPhh+o4bx7xHH+XIwYPRCj2givLKW7KEvCVLAOgydGi0wqy0ivKaft11\nZfovePppTuvVi3ZXXx2zhbyinHIXLnTWu/bk5NCqd29O69UrGuEGraK8SvQbP54t8+eTu3AhZ158\ncRgjip+hFSvkQXpg1y5UlaTUVCQhgREFBSQmJ5OYmsqIggJQZVyjRrTo2ZP/fPppmW3/88kn9Hvh\nBRp36ULOt99GKQNvweYVb6qTV62TTmLXpk0RjrhiVc2p0S9+QXpGBmvfey8KUVesMnl1HDSIpmef\nzd/POScCQ2B2RF7jvNSxIyLC0AULmHXLLeQvW8YVU6eycsoU1n7wQWm/Oo0bs/Xf/y6z7b78fADq\nNmkS0ZiDEWxe8aaqeXX4wx9o1qMHs4cPj2C0walsTndnZ1P75JNJTElh2euv88VDD0Uh6ooFm9fJ\nv/wlF/31r/yjd2+Kwzo2XiJ+jsjtZGeQ9mRnk1q/PonJyaybOZMDu3bRuHNnVk6dyp7sbPZkZ0c7\nxCqxvI75Rf/+XDphAjOGDGH78uVRiDqwyuY0sWdPXunShfcHDSL9t78l49lnoxR5YMHklZiSwlXv\nvMMXDz3EjjV+7+YaYvHz9GU7Ig/CrStXUr9lSxKSkkhMTmbk7t1IQgJJqakMd7+Cj2/Thr25uezL\nyzvuyPvEtDQA9ublRTz2QCqTVzypSl7trrmGy15/nZk33MCKKQFvShcVVclp99atgHOe5mhxMf87\neTKfjxrF4f37o5KDl2DzSkhK4pS2bek3fjz9xo8HQESQhAQeKiriy4cf5usnnghxdPFzRG6FPAiT\nMzJITEmh/8SJbJw9m1XTptFrzBiKDx3iq8cfB2CfW6Szv/6ajoMGldn+zIwMiv77X/KXLo147IFU\nJq94Utm8ut5wAxnPPcf0665j9bvvRivsgKr7u0pITARA3P/GiqDzEuHF9u3LbPvLyy+n99ixvNyp\nE//98ccwRBcbUwuDYYU8CHtycpCEBNI6duTDm26icPNm0jp0YG5mJoWbN5fpu/illzjn9tu59O9/\nZ8HTT3NSejoXPPooi55/PuZmrFQmr4SkJE5p1w6A1Lp1qdWoEWmdOlFcVBTBr7rBqUxe5951F33G\njeOj225jy7//XfrtqbioiIO7dkUjfE+Vyem8e+7hpzVrKNiwAVWlabdu9HniCdZOn07R3r1RysBb\nZfIq/+9sb/funu2hY0fkNU7jLl0oPnSInevXk1qvHqe0a8eW+cffXH5vbi7/vOgiLnrqKW787jsO\nFhay5JVXYvZEU7B51W3WjJu//x4AVaVJ1660GTCAwqwsnktPj3TYFQo2r+7DhyMJCfzu5Zf53csv\nl7ZnzZ3LpN/8JpIhVyjYnBKSkug7bhwNWrVCjx6lMCuLRS+8wIJngrojasQFm5ensF7PEBvj38GQ\nWLiww71F4zNAIvCqqj5Rbr1mRiOwMBrj/tzHxsEVn8GqiTlBzcyrJuYETl4igqpWKzERUXgxyN7D\ngv48EekOvAAk4/ylGKaqiz36BayJ5UV91oqIJOIklgG0Ba4VkTbRjcrb5oq7hF0sxACxEUcsxACx\nEUcsxACxE0dohGXWyjjgYVXtAjzivi+jKjUx6oUc6A5sVNUsVT0MTAUui3JMnrKiHQCxEQPERhxZ\n0Q7AlRXtAIiNGCB24giNsDxYIg+o7y43ALymhFW6JsbCGHkzwHcCbA7QI0qxGGOMKyxj5COBr0Tk\nrzgH0ud59Kl0TYyFQh7UIP2YGBjL18xMxmRmhnSflc0rHDFURaA4IvW7ivTPwl9esfA7qWoMof5d\nxcLPInSqNv1QROYAjT1WjQaGA8NV9X0RuQqYCJR/WHOlfylRP9kpIucCmaqa4b4fBRz1Hdx3TjwY\nY0xwQnOyM/SfJyJ7VLWeuyxAoarWL9enwppYXiwckX8HtBaRVsA2nAeTlrkbTnV/KcYYUxlhrDkb\nRaSXqs4DLgTWe/SpsCaWF/VCrqpHROR24BOcqTavqWpsXWFijDGhcRMwXkRSccZubgIQkabA31X1\nkqrUxKgPrRhjjKmeWJh+GJCIZIjIWhHZICIPhHjfLUTkSxFZJSIrRWS4295QROaIyHoR+VREGvhs\nM8qNZa2IXOTTfraIrHDXVfo2cyKSKCJLRWRmFGNoICLvisgaEVktIj0iHYe7z1Xu9lNEJDUSMYjI\nRBHZLiIrfNpC9rluHm+77QtE5LQgY3jS/X0sF5H3RKS+z7qQx+AvDp9194rIURFpGOmfhdt+h/vz\nWCkivufRwvKziBuqGrMvnK8VG4FWOFdCLQPahHD/jYHO7nIdYB3QBmeS/gi3/QHgcXe5rRtDshvT\nRo59q1kEdHeXPwIyKhnLPcBkYIb7PhoxvAEMcZeTcOa7RiwOdz+bgFT3/dvA4EjEAPwK6AKs8GkL\n2ecCw4AX3eVrgKlBxtAXSHCXHw93DP7icNtbAB/jXPfTMAo/iwuAOUCy+/6UcP8s4uUV9QACBufM\nsfzY5/1IYGQYP2860AdYC6S5bY2Bte7yKOABn/4fA+cCTYA1Pu0DgZcr8bnNgc/cf6gz3bZIx1Af\n2OTRHrE4gIY4f0xPwvlDMhOnkEUkBrcIrAhH7m6fHu5yEvBTMDGUWzcA+Ge4Y/AXB/AO0JGyhTxi\nPwtgGnChR7+w/izi4RXrQyteE+ObheODxDlD3AVYiPM/73Z31XYgzV1u6sZQPp7y7bmVjPNp4H7g\nqE9bpGM4HfhJRF4Xke9F5O8icmIk41DVAuBvwFacs/WFqjonkjGUE8rPLf23rKpHgN2+wxNBGoJz\nVBnxGETkMiBHVX8otyqScbQGfu0OhcwVkW5RiCEmxXohj8iZWBGpA/wLuFNVy9znU50/2WGLQ0R+\nB/yoqksBzylP4Y7BlQR0xfm62RX4L843oIjFISLpwF04R2JNgToi8sdIxuBPtD63hIiMBopUNeJP\nvRCR2sCDwBjf5kjHgfNv9CRVPRfnwGdaFGKISbFeyHNxxuVKtKDsX9hqE5FknCL+pqpOd5u3i0hj\nd30ToOSu9eXjae7Gk+su+7YH+1id84H+IrIZeAu4UETejHAMuPvI0WN3YnsXp7DnRzCObsA3qrrT\nPUp6D2d4LZIx+ArF7yDHZ5uW7r6SgPruN5AKicifgH7AH3yaIxlDOs4f1+Xuv9PmwBIRSYtwHDk4\n/yZw/50eFZGTIxxDTIr1Ql46MV5EUnBOSswI1c5FRIDXgNWq6nuz5hk4J9lw/zvdp32giKSIyOk4\nX/UWqWo+sEecWR4CDPLZJiBVfVBVW6jq6ThjeF+o6qBIxuDGkQ9ki8hZblMfYBXOOHWk4lgLnCsi\ntdxt+wCrIxyDr1D8Dj7w2NeVwOfBBCDO7UzvBy5TVd8nk0QsBlVdoappqnq6++80B+jqDjtFLA6c\nn/+F7s/lLCBFVXdEOIbYFO1B+opewMU4J8A2AqNCvO//wRmXXgYsdV8ZOCfdPsO56upToIHPNg+6\nsawFfuvTfjawwl33XBXj6cWxWSsRjwHoBCwGluMc+dSPdBzACJw/ICtwZtEkRyIGnG9D24AinLHT\n60P5uUAqzlDABmAB0CqIGIa4/bf4/Pt8MZwxlIvjUMnPotz6TbgnOyPwsyiNwf238Ka7zyVA73D/\nLOLlZRcEGWNMnIv1oRVjjDEVsEJujDFxzgq5McbEOSvkxhgT56yQG2NMnLNCbowxcc4KuTHGxDkr\n5MYYE+f+P/pjeYweur3BAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f55de17de90>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get the mass matrix \n",
      "# The model\n",
      "Msig = M.getEdgeInnerProduct(sig)\n",
      "MsigBG = M.getEdgeInnerProduct(sigBG)\n",
      "\n",
      "Mmu = M.getFaceInnerProduct(mu_0, invProp=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "freq = 1e1\n",
      "C = M.edgeCurl\n",
      "A = C.T*Mmu*C - 1j*omega(freq)*Msig\n",
      "ABG = C.T*Mmu*C - 1j*omega(freq)*MsigBG"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Need to solve x and y polarizations of the source.\n",
      "from simpegMT.Utils import get1DEfields\n",
      "# Get a 1d solution for a halfspace background\n",
      "mesh1d = simpeg.Mesh.TensorMesh([M.hz],np.array([M.x0[2]]))\n",
      "e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq)\n",
      "# Setup x (east) polarization (_x)\n",
      "ex_x = np.zeros(M.vnEx,dtype=complex)\n",
      "ey_x = np.zeros((M.nEy,1),dtype=complex)\n",
      "ez_x = np.zeros((M.nEz,1),dtype=complex)\n",
      "# Assign the source to ex_x\n",
      "for i in arange(M.vnEx[0]):\n",
      "    for j in arange(M.vnEx[2]):\n",
      "        ex_x[i,j,:] = e0_1d\n",
      "eBG_x = np.vstack((simpeg.Utils.mkvc(M.r(ex_x,'Ex','Ex','V'),2),ey_x,ez_x))\n",
      "rhs_x = ABG.dot(eBG_x)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 50
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Setup y (north) polarization (_y)\n",
      "ex_y = np.zeros((M.vnEx,1))\n",
      "ey_y = np.zeros(M.vnEy)\n",
      "ez_y = np.zeros((M.vnEz,1))\n",
      "# Assign the source to ex_x\n",
      "for i in arange(M.vnEy[0]):\n",
      "    for j in arange(M.vnEy[2]):\n",
      "        ey_y[i,j,:] = e0_1d\n",
      "eBG_y = np.vstack((ex_y,simpeg.Utils.mkvc(M.r(ey_y,'Ey','Ey','V'),2),ez_y))\n",
      "rhs_y = ABG.dot(eBG_y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Solve the systems for each polarization\n",
      "Ainv = simpeg.SolverCG(A)\n",
      "e_x = Ainv*rhs_x\n",
      "\n",
      "e_y = Ainv*rhs_x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 28
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "NameError",
       "evalue": "name 'e0_1d' is not defined",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-36-4c5017fec386>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0marange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mM\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvnEx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      6\u001b[0m     \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0marange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mM\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvnEx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 7\u001b[1;33m         \u001b[0mex_x\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0me0_1d\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[1;31mNameError\u001b[0m: name 'e0_1d' is not defined"
       ]
      }
     ],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get a 1d solution for a halfspace background\n",
      "mesh1d = simpeg.Mesh.TensorMesh(M.hz,M.x0[2])\n",
      "e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq)\n",
      "# Setup x (east) polarization (_x)\n",
      "ex_x = np.zeros(M.vnEx)\n",
      "ey_x = np.zeros(M.vnEy)\n",
      "ez_x = np.zeros(M.vnEz)\n",
      "# Assign the source to ex_x\n",
      "for i in arange(M.vnEx[0]):\n",
      "    for j in arange(M.vnEx[2]):\n",
      "        ex_x[i,j,:] = e0_1d\n",
      "eBG_x = np.vstack(M.r(ex_x,'Ex','Ex','V'),M.r(ex_y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "SyntaxError",
       "evalue": "invalid syntax (<ipython-input-34-f14ef29654fd>, line 12)",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-34-f14ef29654fd>\"\u001b[1;36m, line \u001b[1;32m12\u001b[0m\n\u001b[1;33m    eBG_x = np.vstack(M.r(ex_x,'Ex','Ex','V'),M.r(ex_y)\u001b[0m\n\u001b[1;37m                                                       ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
       ]
      }
     ],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}