Tutorial for Mag forward problem:

TOODs:

- Change code that we can implement SimPEG's solver class (Now use jacobi-bicgstab)
- Put linear forward problem

simpegPF/notebooks/tutorials/Tutorial_1_Mag forward modeling.ipynb

@@ -0,0 +1,443 @@
+{
+ "metadata": {
+  "name": "",
+  "signature": "sha256:550ca3df316c026a67f55f9de3baa38df81944304bfe8103f9a2d9e8a43c20c7"
+ },
+ "nbformat": 3,
+ "nbformat_minor": 0,
+ "worksheets": [
+  {
+   "cells": [
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "from SimPEG import *\n",
+      "from simpegPF.MagAnalytics import spheremodel, MagSphereAnalFun, CongruousMagBC, MagSphereAnalFunA\n",
+      "from simpegPF.Magnetics import MagneticsDiffSecondary, MagneticsDiffSecondaryInv, BaseMag\n",
+      "%pylab inline"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "output_type": "stream",
+       "stream": "stdout",
+       "text": [
+        "Populating the interactive namespace from numpy and matplotlib\n"
+       ]
+      }
+     ],
+     "prompt_number": 15
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "import matplotlib\n",
+      "matplotlib.rcParams.update({'font.size': 16, 'text.usetex': True})"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 16
+    },
+    {
+     "cell_type": "heading",
+     "level": 1,
+     "metadata": {},
+     "source": [
+      "Forward problem: Magnetics"
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "This is a tutorial for Mag forward problem using simpegPF package. We first start with analytic solution for susceptible sphere in a whole space. Then we solve steady-state Maxwell's equatoins for Mag problem (<a href=\"http://simpegpf.readthedocs.org/en/latest/api_PF.html\">See Doc</a>) using simpegPF package. "
+     ]
+    },
+    {
+     "cell_type": "heading",
+     "level": 2,
+     "metadata": {},
+     "source": [
+      "Step1: Discretize the earth"
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "We use TensorMesh class in SimPEG to discretize the 3D earth (<a href=\"http://docs.simpeg.xyz/en/latest/api_MeshCode.html?highlight=tensormesh#module-SimPEG.Mesh.TensorMesh\">See Doc</a>). Let's visualize discretized mesh on section views:"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "cs = 12.5\n",
+      "ncx, ncy, ncz, npad = 41, 41, 40, 5\n",
+      "hx = [(cs,npad,-1.4), (cs,ncx), (cs,npad,1.4)]\n",
+      "hy = [(cs,npad,-1.4), (cs,ncy), (cs,npad,1.4)]\n",
+      "hz = [(cs,npad,-1.4), (cs,ncz), (cs,npad,1.4)]\n",
+      "mesh = Mesh.TensorMesh([hx, hy, hz], 'CCC')\n",
+      "fig, ax = plt.subplots(1,2, figsize=(12, 5))\n",
+      "dat0 = mesh.plotSlice(np.zeros(mesh.nC), grid=True, ax=ax[0]); ax[0].set_title('XY plane')\n",
+      "dat1 = mesh.plotSlice(np.zeros(mesh.nC), grid=True, normal='X', ax=ax[1]); ax[1].set_title('YZ plane')"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "metadata": {},
+       "output_type": "pyout",
+       "prompt_number": 17,
+       "text": [
+        "<matplotlib.text.Text at 0x3a99ad0>"
+       ]
+      },
+      {
+       "metadata": {},
+       "output_type": "display_data",
+       "png": "iVBORw0KGgoAAAANSUhEUgAAAucAAAFeCAYAAAAmKrHoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3c+SHMd57uGXOlwDo9ENDEBfgAmLXp0Im8LA3ntI61yA\nBpT3hAjtanUEmNqbILR0hI8Ag3sSFGluxX/aGwTnBjQYYC/iLLLKXZ1dWV9mVWVXVvfviUAAld+X\n7yQBIidRXdMjAQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAJjJgaRnkn6Q9Hqg57iuf1lf\n/6q+/sTI/rru+9n4ZUqSHks6nygLAJZq6B78Q+SPSxOtkz0bAAa6LrchPwnUn0n6i6Sj1thX9ZyT\nwJybdf3jaZYoyW30f54wDwCWasgefEfSbwI/HtT9f5lwjezZADDCB3Ib8x1v/F49/gtv/PV6PHRX\npDnQT3UHRmKjB4DG1Hvwd3X/P41f2v9gzwaAkZ7IbfZX6uvmcZbQ3e876j7Q363H3514fWz0ALAy\n1R7c3IT5t+mWJok9GwBGa+7EfFVfN3dS+u68nGv9QH+1vv7vyI/ZfGI51ur5xCdyd/J9XRt9M695\nbr6Ze9nreyj333W5/nWz7k/U/az9gdwnrO9aubci/5sAYFvG7sHN4y/s2QBQqOZOTLPBWS9xNs+r\nN1+Y9FjuQP/XkR+v+cfAD5J+L/fs4yfqfgbe3+jf0uqTym+8uf4XSj2o857UH+cXWv23+i8LN18k\n+xe5O0nvtnK/EgCUY8we3Bzk/a8p6sOeDQAzaA7msV/M2Xwh0ZCXRpt3B/CfaW824dPWmL/RP1T3\nnf2H9dyuMX9tzce57vX+WZufrG51rAkA5jZ0Dx7ynDl7NgBsWXMnpe/dW3yXtXrJ8c9K+yLQZk6o\n1r7rEfv84tfafMeBZqP319Y8V99+u8e+T24/yH0iBIBSDNmDmz3xN4kfiz0bxfnR3AsAMnso6aXc\nFxRdVfdzhL7nWm3wv5H0IvFjhl52/Far5yhDrso9M3lX7hNB3/u1v4xY29X653fU/f6/L7X5bCQA\nzCl1D74p9xaMX0v69YCPx56Norw69wKAjO7KbZJ35Tbsv5HbRB9K+oMx93n980W21W36ldxLnM/k\n7oz8Xu7ly9ty6x7jXv3D94q2+98IADFi9+Brcjddnmn90ZBtYM9GFhzOsauuyT2f951Wd1Leqa8f\nyt0Ned49dbSrPeOhR2sO5Db5x5L+0au9MmItT1sZf+qo39TqO6UCwJIcyN1oeSnpbaW/ytlgz0ZR\neKwFu6p5nOXt1tj3kt6T21TvZ/zYr2nzC3buyr0U+fvAnMP652+98QO5ZxJfjljPp3Ibuv9S6125\nO04/HpENAHN5KLev/qukz0bksGejKNw5xy66K3dn/K427zy8L+nncm+BdSLpUYaPfyH3cuQNuWcZ\nb8i93PqdpN96vc0dlqf1j1/V11/J3bW5LfeFUa/IvRLwoVZ3/GPvzrwj9yzm15L+01vTQ437pAYA\nc7ip1WMsr8jt9yH/V/2vlLJnA0BG12R/A4orst8F4FTuq+39t9eyNO+V+zOt3jv3v9X9lfefaP0r\n/6/IPbd4Xv/4WO69fa/IvbTafu/eB9p8NwDJ3bH5i9a/8l9yd4Ca99lt1jT1dzwFgKlYe3DztoJ/\nUfcXTv7Qqh/1fBz2bADYcc1GDwAoH3s2isMz5wAAAEAhOJwDAAAAheBwDgAAABRizHtx7qIxb30E\nAHP6QtLfz72ILWPPBrBUwT2bw/m6l1I1Yvrnkt4M1CqNy06Z/7ncn3lffyjPH7eu+3qa348xGayD\ndQz9uH3j7frfKfz3NjVvzPy+/SM2e+/2dPbsSf/+hP5OpmSwDtYxdB1946k9Q3pTM/Lt2TzWAgAA\nABSCwzkAAABQCA7nkzqaewG1o7kXUDuaewG1o7kXUDuaewG1o7kXUDuaewG1o7kXUDuaewF76Gju\nBdSO5l5A7WjuBdSO5l5A7WjuBdSO5l5A7WjuBRTmKFsyh/NJXZl7ATXWsY51rGMd61jH/irl95x1\nrGMd61hHmfL9fryaLRkAAOcDSb/0xm5JeirpsL6+n1gHgJ3EnXMAQE53Jf20Y+xrSY/kDt2vSTpJ\nqAPAzlri4fyDjrFbchv3af0jtQ4AmN5Vdb8X+amkz1rXjyW9k1AHgJ21tMM5d2AAYDmuyx2s2651\n9D2TdBxZB4CdtqTDOXdgAGA5rkt6oM1vsnEo6dwbu6h/vhRRB4CdtqTDOXdgAGA5DiQ9D4wfemPN\nYfwwog4AO20ph3PuwADAcpzIPUrY5aJjrDl0n0fUAWCnLeWtFHPdgXkx1QIBAJLcm/92HbAb53J7\nc1tz/SKi7vm89esj8V7MAMr0vaSzqM4lHM65AwMAy3FN7muEmscK35A7XL8rt5d/o829+VCrxxat\nuufNsesFgC24ovWbB18EO0s/nG/5DozEXRgAyxB/F2bL/JspN+UO679tjX2o9Rsvx5LuJdQBYGeV\nfjjf8h0YibswAJYh/i7MjE4lvSW30Hfl3s72uaTbWn3/iauSnkj6qDXPqgPAzir9cM4dGABYrvv1\njy7vG3OtOgDspKW8W4u0eQfmcj1+W+7AfiJ3p6XrDkxfHQAAAChC6XfO27gDAwAAgJ22pDvnAAAA\nwE5b0p3zLakKzk6db/WH6v64dT1kDhlk5M6wxmPrY/unng8A2GX+d9zcdy/zfeKsNC47db7VH6r7\n49b1kDlkkJE7wxqPrY/tn3q+lb13e/rLuRcAACN07tmvbnsVAABMp8qYOyY7db7VH6r74119Vg8Z\nZJSQ0Tee2jOkN2dGKLcbz5wDAAAAheBwDgAAABSCwzkAAABQCA7nAAAAQCE4nAMAAACF4HAOAAAA\nFILDOQAAAFAIDucAAABAITicAwAAAIXgcA4AAAAUgsM5AAAAUAgO5wAAAEAhXp17AeWpCs5OnW/1\nh+r+uHU9ZA4ZZOTOsMZj62P7p54PANhlr8y9gMK8zPeJs9K47NT5Vn+o7o9b10PmkEFG7gxrPLY+\ntn/q+Vb23u3pL+deAACM0Llnv7rtVQAAMJ0qY+6Y7NT5Vn+o7o939Vk9ZJBRQkbfeGrPkN6cGaHc\nbhzOAQBTO5B0KulC0mv12G2v55akp5IO6+v7iXUA2ElLOJyzyQPAsvxa0nut66/k9vFm770r6WNJ\nn9XXdySdSHoUWQeAnbWEd2v5taT35Tb125KO5Tb5xl1JX8tt2vflDvAnCXUAwLROJP2idf1U0o3W\n9alWB29JeizpnYQ6AOysJRzO2eQBYFmOJf2udf2apD/Wv77W0f+snhNTB4CdtoTHWo4lnbWuX5P0\nH/Wv2eQBoDxnrV9fk/SDpN/W14eSzr3+i/rnSxH1F5OtEgAKtITD+Vnr12zyALAMlyX9s6S3Jd1s\njR9o9fU/jWafPoyos28D2GlLOJxLbPIAsDTP5b7O577c1/18UP/6oqO32afPI+qez1u/PpJ0Zcha\nASCz77V+vzlsKYfzLW3yEhs9gGWI3+hncKD1/fcDSffk9u3zuu73S+6GiVX3vDlupQCwFVe0fqb8\nIti5hMP5Fjd5af03K/wbN0y15flWf6juj1vXQ+aQQUbuDGs8tj62f+r5xTuW9IncXtvss813wbsk\n6Rtt3jg5lPtifUXUAWCnlf6tnrs2+ZtyB/Rm7Fzrj64cy72v+T/W11a97WW+T5yVxmWnzrf6Q3V/\n3LoeMocMMnJnWOOx9bH9U8+3sovY0y/LvS/5v7TGHsp9vdDP6+s7kr7U6n3L78i9m8tHkfXGyykX\nDgBb1rlnv7rtVST6Uu4uefsu9w25jb4Z+1Dr35ziuJ6jyDoAYDrP5fbdW/X1TyQ9kfueFY3bdf1E\n0tW6/lFCvaWaat0duWOyU+db/aG6P97VZ/WQQUYJGX3jqT1DenNmhHK7lX443/ImDwCYwLf1jz7v\nj6wDwE4q/XAusckDAABgTyzhO4QCAAAAe4HDOQAAAFAIDucAAABAIUp4262S8LZcAJZs3/Z09mwA\nS7bIt1KcQZUxd0x26nyrP1T3x63rIXPIICN3hjUeWx/bP/V8K3sfVRlzx2Snzrf6Q3V/vKvP6iGD\njBIy+sZTe4b05swI5XbjsRYAAACgEBzOAQAAgEJwOAcAAAAK8ercCyhPVXB26nyrP1T3x63rIXPI\nICN3hjUeWx/bP/V8AMAu27ev7Le8zPeJs9K47NT5Vn+o7o9b10PmkEFG7gxrPLY+tn/q+Vb23u3p\nvFsLgCXj3VoAALumypg7Jjt1vtUfqvvjXX1WDxlklJDRN57aM6Q3Z0YotxvPnAMAAACF4M75hqrg\n7NT5Vn+o7o9b10PmkEFG7gxrPLY+tn/q+QCAXbZvzydaeOZ80S+JkUGGLzQeWx/bP/V8K3vv9nSe\nOQewZDxzDgDYNVXG3DHZqfOt/lDdH+/qs3rIIKOEjL7x1J4hvTkzQrndeOYcAAAAKAR3zjdUBWen\nzrf6Q3V/3LoeMocMMnJnWOOx9bH9U88HAOyyfXs+0cIz54t+SYwMMnyh8dj62P6p51vZe7en88w5\ngCXjmXMAwNbcqn9+Q9KXkt7vqD+VdFhf30+s16pxqwyqRmanzrf6Q3V/vKvP6iGDjBIy+sZTe4b0\n5swI5XZbyuF8S5s8AGACdyTdbl1/Vf/c7N13JX0s6bNW/4mkR5F1ANhZSzicb3mTr8avOGhsdup8\nqz9U98et6yFzyCAjd4Y1Hlsf2z/1/OJdlvRnb+ye3F7c7Nunkt5r1R/X148i6wCws0p/PvGypJta\nv1N+KrfJN3fBz1u/lqTrcpv4P0TW23jmfNEviZFBhi80Hlsf2z/1fCu7iD39qqQn9c9n9dLine truncated
+       "text": [
+        "<matplotlib.figure.Figure at 0x3b6b310>"
+       ]
+      }
+     ],
+     "prompt_number": 17
+    },
+    {
+     "cell_type": "heading",
+     "level": 2,
+     "metadata": {},
+     "source": [
+      "Step2: Compose suceptibility model: susceptible sphere in whole space"
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "- $\\mu = \\mu_0(1+\\chi)$\n",
+      "- $\\mu$: magnetic permeability\n",
+      "- $\\mu_0$: magnetic permeability of vacuum space\n",
+      "- $\\chi$: magnetic susceptibility"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "from scipy.constants import mu_0\n",
+      "mu0 = 4*np.pi*1e-7\n",
+      "chibkg = 0.     # Background susceptibility\n",
+      "chiblk = 0.01   # Susceptibility for a sphere\n",
+      "chi = np.ones(mesh.nC)*chibkg"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 19
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "sph_ind = spheremodel(mesh, 0, 0, -100, 80) # A sphere is located at (0, 0, 0) and radius of the sphere is 100 m\n",
+      "chi[sph_ind] = chiblk # Assign susceptibility value for the sphere\n",
+      "mu = (1.+chi)*mu0"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 20
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "fig, ax  =  plt.subplots(1,2, figsize=(12, 7))\n",
+      "indz = int(np.argmin(abs(mesh.vectorCCz-(-100.)))); indx = int(np.argmin(abs(mesh.vectorCCx-0.)))\n",
+      "dat0 = mesh.plotSlice(chi, grid=True, ind=indz, ax=ax[0]); ax[0].set_title(('XY plane at z=%5.2f')%(mesh.vectorCCz[indz]))\n",
+      "dat1 = mesh.plotSlice(chi, grid=True, normal='X', ind=indx, ax=ax[1]); ax[1].set_title(('YZ plane at x=%5.2f')%(mesh.vectorCCx[indx]))\n",
+      "cb0 = plt.colorbar(dat0[0], orientation='horizontal',  ax=ax[0], ticks = linspace(0, 0.01, 5)); \n",
+      "cb1 = plt.colorbar(dat1[0], orientation='horizontal',  ax=ax[1], ticks = linspace(0, 0.01, 5)); \n",
+      "cb0.set_label(\"Suceptibility\")\n",
+      "cb1.set_label(\"Suceptibility\")"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "metadata": {},
+       "output_type": "display_data",
+       "png": "iVBORw0KGgoAAAANSUhEUgAAAv0AAAGjCAYAAAC/qqoOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3bGPHEfa5/mfJL7ArcPut8fas5qt+QNGnJnXuQXeldh8\nxzlrmhqts1hHTc1rHHCGKFFeGgeouRr/JcUx97ArcinjgF1AIkcaYYEzRpQ0zmGBG4rqde7OUbPJ\nA27v8IrkGZGpyorKyIjIzMiMrPp+gEZ3xvPEU9HV1ZHRmVnZEgAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAABstG1JjyQ9k/SKI2e/jH9Zbr9Tbn/qqf1Vmfda/2FKku5KOhmoFsLclnTLk3NN\n0rcyP5tP5X4dSdLlMveZpAeSrkSMZa8cT/V6bev/rOXjbMRjAnPXdb5u+x1K8fvE/D6OmPm6b/++\njwUggQtaLKKaPJL0VNJure1+2efA0edyGf9kmCFKMjuF7wesN7XqOfrZ1ANx2JMZ30ctOV/Vct6X\neQ09k3lN2W5o8cfj+1q8ho4Cx1It9j+S9LbMTqRpMbNde5x/aPgANk2X+fpI5ve06eNWmf90wDEy\nv6cXM1/37d/3sQAkdF3NC7Bqofam1f5K2e46MlP9oTDkUVV2CuO4JHOEpr7IblKN/+1a25bMa8J+\nXVR/QNh/BFYL9y3PmKrX4a+t9iOtLmbOO3KBTTX0fP1tmT/k7xjze1ox83Xf/n0fC8AIqr/Ez5Xb\n1WU9rqP11YLL/kPhmlZ/4YewrjuF3E55Vkfwnqp90f+Vmo/0Va+L+hGdatG+a+VeUNji4ZGaf/Zb\nZd3rtbZLymtnC+RgqPm6+l0e+qwZ83taMfN13/59HwvACKqjQffL7epoTtvRnxMt/6FQHdH9S+Bj\nVjuhfS2u6Xyg5UVcpWmnUPWrX+d9XatHjm/LfF9b5dfVuF3XGW7L7Ny6XH8eMqa7Wr4utu35arue\nNuaa+FjVgtq16H+k5j8Iqz8W37RyQ18TTR6oeZFRXcpTH0e1iNmq5QDoP19XC1nm92Hm9+r5t983\nVZ2tHPIPq5j5um//vo8FYCTVX+LVZOg7Alu9H6C6rvquzB8KoUdZqz8y6tf+VZd82O8xsHcK1RHd\nv2hxranrOu9bZb0H5eO8qcX3ap9urN7c/FRm0q1fP35f7ULHdK72+G+q/fn6tfVxUH4fTzXcm6Sb\nVAtq16LfFat2ZO+35A61EK/epFg/Snm7bKs+Vz/jpsUCsEn6zNfV77X9/q42zO/++f2KVi9R/Fb9\nDpI0iZmv+/bv+1gARlQt+EPfhFu9qavLad9qUWb/5V9NmIe1NnuncFvNZyKqxV5Tmz22ptONt8vH\n2bVyq8n5UG4xY+p6zWfVL/XE2bbor2JtR9//wdq+LrNIr46QVTvKcw01QlTPw/dafr6rN5B9KbND\nf1OL16jrzerApug6X3e5jp/5PWx+vy/zx8mWzJnKmANnIWLm6779+z4WgBFVf4nHLJCqN+g0LcB8\nqj6uWP3IS+g1n03XE1aTsj226nRj/Yh526TUdCo2RNOYuiz6q9O+fe+KtG19uHJ8R2zaJvaPrNzq\nZ/22zMKhfiQu5gj8nhanz7/X6vN3Rc3XJzctNIBN02W+rubP2AMNzO9h8/u5Mrc6YDH0++Fi5uu+\n/fs+FoARVZPX+1ocnQ1RHSWJnazaFrBfaXkn0LRT2JOZXK9p+TrKpp1C0xuL7J1CfYHa9PG0Zbyx\nY4pd9FenpV2nfdvGXX38TM3fY9M1lkNd3lPVabos4KCMhdy2U1pczvO0fOyYPzB971EANkXMfF3N\nU1/6Ehswv4fP74dlft/53SV0vh6if9/HwkjOTD0ATOqazJuerkl6T9LPZSau25L+4On7uPx8mmx0\nq96RWSw+kjk685HMxHlVZtx93Cg/bC+o/XtMOabqKNBFR/x8QI0/l5/3G2p30XSWoGqrduDV83VP\n0rGVe6f8HHKHi9syfyR8K+l1Lb6XUNVrlDf2YtOFztfnZQ78PNL4d1zZtPn9p+Xnn8gcoHhsxWPm\nd5eQ+Xqo/n0fCyNg0b+5zssc/flWZsEvSW+V27dlTj/ak9BQ9lraXZcYbctMvncl/cqKvdBjLA9r\nNZom0MtyH/FKNSbJ/Ax2ZRb8x44cV3uTz/oNR5JZxDf9AfJGLV75WmZn5uJbfLwjs+C/Xavf5Jyk\nfy3p32nxB0Wl2uE8FACfbZmDPc9l/sh+0rEO87vfvsz+94bMfve2pL+zco57PkbMfN23f9/HApBY\n9SatpuujQ65zrE5lxt6KqzotaV9n3XTv6Prp3+o0oX1ZyLYW30td6OlfaXE3Bvvo87WG3LrYMYXe\nx7npDjVN+p7+tfku76kuzan/7Krv1T5FXZ26to8WVs9pyO3iYm4T+ECr7xOo3rjI/fux6ULm6+rS\nlT6XYjC/++f36rLN6o+Nprv5SP3n95j5um//vo8FIKFqsnNN7r5/4S71W/RXbyq7JbPArXY29uRQ\n3ee5Uv0zsSOZW6m9o8V9oJ/JTJ7Vwq/p7gpS807hXMuYfNeDx4ypmhhvyf3cVuN7UOZcsj7qi+if\nBXzECHnTVfXaOJL5Xqv3hDTtOKvnoLqLT/Wc2tfQXqrlSYud7QMtTsvbH/WdS/UHxkk5rmtaXBrF\nXSMA/3xdxavf7WstH21vwmd+b5/fJenfa/X9TvW7+VSGmN9D52t7Do7tH5sLYCTV3WDa/vqu7izQ\ndqeHQ5lf6C6L/o9kJoJqkviLmhdnn2r5WsBzMhNq9a+9P5GZ+M5pcR/73TL3ltxHgpomoi0t7v1c\njSnkTW8xY9rS4qiT65RytYCt/juu/dH3Lj5tQhb9WzI7hmon+qXck3qV+6CW2/ScVjvL6nGrHbfr\nOXiq1TFe0PLP4UvF3WoQWGe++bo62uz6nav/7u22PA7ze/v8Xi2u7ceu9rlDz++h87U9B8f2j80F\nsCFCjq4AAOaH+R3I2ItTDwAAAABAWiz6AQAAgDXHoh8AAAAAAAAAgDnr+w8m1s3zqQcAAB19Iemf\nTz2IkTFnA5ir0edsFv3LnktFj+6fS3rVESvUr3ZM/89lXktt+a56drtvuy2nej761GAcjKPr47a1\n1+N/K/fvbWy9Pv3b5o/Q2hs3pzNnD/r74/qdjKnBOBhH13G0tcfmdMmNrTG/OZtr+gEAAIA1x6If\nAAAAWHMs+ge1O/UASrtTD6C0O/UASrtTD6C0O/UASrtTD6C0O/UASrtTD6C0O/UANtDu1AMo7U49\ngNLu1AMo7U49gNLu1AMo7U49gNLu1APIzO7UA4h2ZuoBrJdzUw+gxDiWMY5ljGMZ4xjBdUm/tdqu\nSHooaafcvhkZH0AuzznjWMY4ljGOPM3v+Zjjoj/TnQcAoME1Sb9oaPtE0mfl9pGkA0l3AuMAgEhz\nu7zHtfP4SmZncFPSyzI7h9A4ACCNPTXfVvNQiwW9JN2V9FZEHAAQaU6LfnYeADAvF2Tm3LrzDXmP\nJO0HxgEAHcxp0c/OAwDm44KkW1q9D/WOpBOr7bT8fDYgDgDoYC6LfnYeADAv25IeO9p3rLZqnt4J\niAMAOpjLop+dBwDMR9ubbk8b2qr5+CQgDgDoYA6LfnYeADAf59Q891ZOZA7I1FXbTwLiAIAOcr9l\nJzsPAJiX8zI3XqjeU/VLmXn3bZkDOF9rdV7f0eI9W7645fPa17ua472zAWyC7yQdTzqC3Bf9I+88\nJHYgAOZh+h2Ig31m9rLMPP67WtuHWj6Luy/pRkS85tVegwWAcZzT8pryi9FHkPuif+Sdh8QOBMA8\nTL8DCXAo6ZLMQN+W+V8pjyVdlfmniQcyc/oDSR/X+vniAIBIuS/669h5AMC83JT7P6B/4OnriwMA\nIsxp0c/OAwAAAOhgDnfvAQAAANADi34AAABgzc3p8p6RFBnXju3vy3fF7Xbfdpc+1KBG6hq+9tB4\n3/yh+wMAEO+FqQeQmefpdsiF+tWO7e/Ld8Xtdt92lz7UoEbqGr720Hjf/KH7+2pv3Jz+fOoBAEAP\no87ZZ8Z8MAAAhlUkrNundmx/X74rbrc35flyqEGNHGq0tcfmdMlNWcNVd1xc0w8AAACsORb9AAAA\nwJpj0Q8AAACsORb9AAAAwJpj0Q8AAACsORb9AAAAwJpj0Q8AAACsORb9AAAAwJpj0Q8AAACsORb9\nAAAAwJpj0Q8AAACsORb9AAAAwJp7YeoBZOb51AMAgB42bU5nzgYwZ6PO2WfGfLB5KBLW7VM7tr8v\n3xW3233bXfpQgxqpa/jaQ+N984fu76u9iYqEdfvUju3vy3fF7famPF8ONaiRQ4229ticLrkpa7jq\njuvM6I8IAFh325IOJZ1Kerlsu2rlXJH0UNJOuX0zMg4AiDCHRT87DwCYl/ckvVvbvi8zj1dz7zVJ\nn0j6rNw+knQg6U5gHAAQaQ5v5H1P0gcyO4urkvZldh6Va5K+ktkZ3JT5w+AgIg4AGNaBpDdr2w8l\nXaxtH2qxoJeku5LeiogDACLNYdHPzgMA5mVf0u9r2y9L+lP59fmG/Edln5A4AKCDOVzesy/puLb9\nsqR/W37NzgMA8nNc+/q8pGeSfldu70g6sfJPy89nA+JPBhslAGyQOSz6j2tfs/MAgHnYkvQbSa9L\nulxr39bi/VWVap7eCYgzbwNAB3O4vEcyO49DmTdzDbnzAACk8VjmfVR/J3OpT/VerNOG3Go+PgmI\nAwA6mMORfmmx87gp86bc6+XX7DwAID/bWp5/r0u6ITNvn5RxO18yR/F9ccvnta93JZ3rMFwASO07\nLV+8Mr45LPpH3HlIaf9ZQt/asf19+a643e7b7tKHGtRIXcPXHhrvmz90/+ztS/pUZq6t5tLine truncated
+       "text": [
+        "<matplotlib.figure.Figure at 0x42450d0>"
+       ]
+      }
+     ],
+     "prompt_number": 21
+    },
+    {
+     "cell_type": "heading",
+     "level": 2,
+     "metadata": {},
+     "source": [
+      "Step3: Set up an airborne MAG survey"
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "We have discretized 3D earth and generated suceptibility model, which means that we have discretized earth and physical property distribution. We can compute magnetic fields everywhere in our domain by solving parial differental equation (PDE), but our measurements are confined to finite locations. Therefore, we need to project computed fields $\\mathbf{u}$, which is defined everywhere in our domain to certain locations where we have receiving points. For instance in airborne mag survey these are the points where a plane or helicopter measure earth magnetic fields. This projection can be expressed as:\n",
+      "\n",
+      "$$ \\mathbf{d} = P(\\mathbf{u})$$\n",
+      "\n",
+      "where $P(\\cdot)$ is a projection from computed field to the measured data, and $d$ is the measure data. "
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "Let's assume that we have a survey area: 400 m $\\times$ 400 m. We have 21 lines of airborne MAG survey, and we measure magnetic fields for every 20 m on each line. A pilot for this helicopter is really talented so that the flight height is constant for 30 m above the surface.  "
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "xr = np.linspace(-200, 200, 21)\n",
+      "yr = np.linspace(-200, 200, 21)\n",
+      "X, Y = np.meshgrid(xr, yr)\n",
+      "Z = np.ones((size(xr), size(yr)))*30."
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 22
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "fig, ax  =  plt.subplots(1,1, figsize=(5, 5))\n",
+      "indz = int(np.argmin(abs(mesh.vectorCCz-(0.))));\n",
+      "dat0 = mesh.plotSlice(chi, grid=True, ind=indz, ax=ax); ax.set_title(('XY plane at z=%5.2f')%(mesh.vectorCCz[indz]))\n",
+      "ax.plot(X.flatten(), Y.flatten(), 'w.', ms=5)"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "metadata": {},
+       "output_type": "pyout",
+       "prompt_number": 24,
+       "text": [
+        "[<matplotlib.lines.Line2D at 0x56e54d0>]"
+       ]
+      },
+      {
+       "metadata": {},
+       "output_type": "display_data",
+       "png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFeCAYAAAC7EcWRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnc9vHEeW57+c7ls3RDZ9XWBIav+AFscewAPM9tiivHs2\n5WnMbYE1ZR/3QNpqYLFTAwNjcWWdRz/6Pmtp5bskj73GAD3AyLL73rIsLLBHU5QwDRtYW7mHyHBl\nRkZlRP6qyCp+PkCCFfFevHiVrHoZ9TIyQgIAAAAAAAAAAAAAAAAAAAAAAAAACLIm6Ymk55LOzNDZ\nyeX38/I7efluwPaDXO/V7m5Kku5JOurJFtRzQdJXMv+/h5IOGrTdknRL089VXfvnNcepNo4DLBpn\nNf2i+Hgi6QdJG4W6z/M2uzPaXMjld/pxUZIJwN/0aC819hz9MrUjDtc0veC+r+n/+lJE2y1NA++H\nkvZlLtS+C/ZaoZ9/8BwAJ4ar8n/J7JfxTaf+TF4/a0Rqg3afoxgC8PBsyX/htEF0NdDefl5ed+ov\nqXrB3p6hC3AieSjzhdjMyzb1MGsUa79UbtA+zOv3e/ZvWQPwrNRPCmwA3XDqz8pcUEPB8on8/6PV\n3O7VQt15je8CBJAMO6r9PC9/pfAo9kjloG1HUH+I7NMG8B1Nc7wPVf6iWnwB2LYr5huvqjpSuyXz\nvlbz19bvu/IHwDWZYNQmDxrj0z2V851156suT9okNxvDk4AvIR7Knz6w6YbixdxeqFcLOgAnGjuq\ntYEnNOKx+WOb37snE7RjRzU24Nuc4fua/tx1c9JuALYjqD/k7Ypt3Xzjzdzew7yfNzV9r24axd6Y\n/EEmmBTzmJ+rnlifNgv9v6n68/W6c+zm7+MH9XeD02L/D5a+gqK9cVv8VXQrr7N/7f/CdwEFODHY\n4Bt7A+1mrm9/vja5gWK/eG6O2QanvUKdG4BvyT9Ct19oX53rm+3nrKP7jao/ww88Prk08altDti2\ne79huxB2lHpVJmDaEby9eGzOblqL9fcblc+LnSVzXyYwv6npZ2nWDWGApcamEJp8CVY1/UnvfslC\n2DazZMURZ2wO+IFMECxiA6Drm81zF0eSdReR5zJBoik+n9oEYHvjquvskjXnkMr/+29kguLrKv9S\naDIy3dI01fKNqu/zQP77BL6LL8CJwAaK91W9aVKHHR02vfFWF0weqBxwfQF4SyaQHaqcV/UFYLdO\nqgbgYhDyHT/U+NvUp6YB2KZGZuVo6/y2xy/lf49vanqjzJ1yKJm0R+xUNGmacvhBJqXR5KJs/fgw\npLis/DS1A5CEQ5kbUoeSfiPpz2SCxC1J/xRo+zT/ezyYd1XekQkIT2RGpR/KjJouyvjdhWv54bKi\n+vc4pE/2J/u5GfLtCBu/z//uOPWfS3qWv/5Y0mNHfjv/GzNj45ZMwP5K0huFPmOxnyVuysGJwf60\nLY6uNhX/09OO5txcboi6GQBPNH36TiqPgH131S02F10kdgRsfZo18r+g2UGoqU9NRsDW/7qbbm5a\nwXeE+FyzbzTGjErfidTb1DRQu9jzyMMYcGKwU858ebqYvGeXAOzL9/nmEhcDsP0Z7f4kXtP0vRRp\nEoDtrAU30B56dIs09Sl2HrBvBoGP2BREHXuq3pSUpu899P9tMo3N3mdwL+72YsX8YDgR2C/XrLvq\noceOpW4B2N7AuykTbGze1P0iu2tB2AdHLslM/3pH03nEz2UuHvbL7ZuFIPkD8GaNT6GRXROfbF71\npmafW+vfw1znvHMUA+UvI44YrK92NoR97+7I/nxBT5pegB5qmsJxj+KF1gb7I5nzdahpmoXRL5wI\nfKkHF5uKqJvhsCczwmsTgD+UCYA20P9B/i/gXZVvwm3KBK+j/LgjE2Q2NZ0nu5Hr3tTsEbDvp/2q\npnOHrU8xNxib+LSq6Wj7vmsoxwapH+Qf0fa51oZlVSao2vdup4m52AuIvSjZi8UsX39Q9QJ2VuXz\ndV88mgwwN0703W4AH3+S2gEAgJMKARgAIBEEYACARKykdmBkZKkdAICl5DNJf+VWEoDLZNKkZdNP\nJb3iqZ+ovc3YtiE9n9yta1u277uLvZR9t/VtHn3P8iVF36H+u9oO1VvZr+T/njWx07btrO94rM1q\nvCUFAQCQCAIwAEAiCMC9sZHagURs0PeJ6jt1/8vVNwG4N9quYb3opHzf9H3y+l+uvgnAAACJIAAD\nACRiERdkvyrpbafuQNIjSet5+UZDOQDA3Fm0EfChpBc9dQ9kVvK/Iem0ykv+heQAAElYpAC8Jf+T\nanuSPimU70l6q4EcACAJixSAz8oEzyK+vbGeaLoPVkgOAJCMRQnAdjFn91G+dZV3TpCmGymeipAD\nACRjUQLwmqY7qLr1606dDbjrEXIAgGQsQgDe1XSrbBfftuE2sB5FyAEAkjH2aWib8gdRy5GqW3Db\n8rMIuYdPC683lP6pIwBYPL6W9DioNfYAvC0z+8HeTHtJJoDuy4yKv1A1QK9rerMuJPfQdrk5AADL\npsqDt8+8WmMPwG7q4YJMQP6gUHdd5TTFjsy22LFyAIAkLEIO2LIn6bzMZWVfZkttSbooE5R3ZZ54\neyjpo0K7kBwAIAljHwEXuaHZjxBfDrQNyQEA5s4ijYABAJYK9oQrw6acADAUlXi7SCmIOTEZwF5b\nm7FtQ3o+uVuXsjwmX8bs25h86du3UH1I1kW3z7Z1NquQggAASAQBGAAgEQRgAIBEEIABABJBAAYA\nSAQBGAAgEQRgAIBEEIABABJBAAYASAQBGAAgEQRgAIBEEIABABJBAAYASATLUZZhOUoAGAqWowwz\nGcBeW5uxbUN6Prlbl7I8Jl/G7NuYfOnbt1B9SNZFt8+2dTarkIIAAEgEARgAIBGLkIJYk9kR+VjS\n6bzuoqNzIOmRpPW87G7eGZIDAMydRRgB/0ZmV+MbMoF3RyYgWw4lPZB0O9c5LbMFfawcACAJixCA\ndyW9WSg/knSuUN6T9EmhfE/SWw3kAABJWIQUxI6kx4XyaUn/mL/e9ug/ydvEyAEAkrEII+DHhdfb\nkp5L+iAvr0s6cvSP87+nIuQAAMlYhBGwJK1K+mtJb0i6UKhf0/TGmsUG3PUI+bN+3QQAiGdRAvBT\nmRtoN2RuqF3NXx97dG3APYqQe/i08HpD0mZTXwHgxPO1yj/eF5c1p7wnk4aQpimJItsN5C4ZBwcH\nx0BHhbGPgHck3ZUJwjZdYJ+nPiXpC1VHuesyMx0UIfcwaetrjb22NmPbhvR8crcuZXlMvozZtzH5\n0rdvofqQrItun23rbFYZ+024+5KuqZyrPSfpVqHuusrzenfyNoqUAwAkYewj4KcyAfQgL78g6aHM\nwxmWi7l8V9JWLv+ogRwAIAljD8CS9GV+1HG5oxwAYO6MPQUBALC0EIABABJBAAYASARbEpXxztUD\nAOgBtiQKMxnAXlubsW1Dej65W5eyPCZfxuzbmHzp27dQfUjWRbfPtnU2q5CCAABIBAEYACARBGAA\ngERwE64MN+EAYCi4CRdmMoC9tjZj24b0fHK3LmV5TL6M2bcx+dK3b6H6kKyLbp9t62xWIQUBAJAI\nAjAAQCLIAZchBwwAQ0EOOMxkAHttbca2Den55G5dyvKYfBmzb2PypW/fQvUhWRfdPtvW2axCCgIA\nIBEEYACARJADLkMOGACGghxwmMkA9trajG0b0vPJ3bqU5TH5MmbfxuRL376F6kOyLrp9tq2zWYUU\nBABAIhZlBGw35XxJZqdkd4+3A0mPZLacl6QbDeUAAODhklP+XNOALEmHkl519HcbyItkHBwcHAMd\nFcZ+E25V0gWVR7x7MkHVjmaPCq8l6aykdyW9FikvkpEDTlEeky9j9m1MvvTtW6g+JOui22fbOpvV\neDv2HPALMsF2o1D3RNJa/nrb0+aJpJ1IOQBAMsaeA34kE0QfF+rOSbqXv16XGeEWOc7/noqQP+vL\n0ZPGysqK9vf/QlmW6cqVf1GW2Xrl9fu6ckWFer9+n7Zm6Zdlma5cWXHaHCjLXh5h365+3TlsaquZ\nfqh/OBmsyQTUjbx8XtUAuybpea4TkrukzhEtzHFwcJB999132bfffpsdHBy0ru/TVso+lv39hWQc\nUUeFsY+AXW7K3FB7nJePPTrF3HBI7uFXhdcbkjabeVhhomXMAWfZvynLfpK/3pFJ00+UZS/n9d/n\n9T+r1e/T1ix9I9sptMkKfb+cl39aaFPUT9m3q193DpvaaqYf6p8csMvXKv9w/8yrtUgB+FJ+/L5Q\nd6RpPthiy88i5B5e6eLjieHKlStaWXntx5+j0/p/0crKirLsrq5c+XlQv09bs/TLskxXrlyR9LeF\n+o+VZTsj7NvVrzuHTW010w/1Dy6bKg/e/AF4UdhVeSrZmcJrdyS7I+lOA3mR1D9RODg4lveosAgj\n4B2ZtMHHMqPXdUm/lvRlLr8uE6BvF/SvFdqH5A6TXpwu22trM7ZtSM8nd+tSlsfky5h9G5MvffsW\nqg/Juuj22bbOZpWxB+A1SXfz18Wgeavw+qLMgxm7krYkPZT0UQM5AEASxh6AjxU3V9l9NLmpHBrC\nNDSmoTENDfomdY5oYY5lmUK1DH2n7oMj+qgw9hFwAiYD2GtrM7ZtSM8nd+ualZmGNo++XX2moZED\nBhDT0JiGxjS0PiAAQyuyLNPly7/z1Cuv/0Dl0ZNfv09bs/TLMqn4a9DUX9Z0xDqmvl39unPY1FYz\n/VD/0I6xr4Y2b7x5GgCAHmBLojCTAey1tRnbNqTnk7t1zcorK3+n/f07hTviRj69i25+ik/vrvv1\n+7Q1S9/IJoU7+P9JWfa3uS1pf/+PP/4UN22K+in7dvXrzmFTW830Q/2TA46xWYUADK3Y39Line truncated
+       "text": [
+        "<matplotlib.figure.Figure at 0x482ae10>"
+       ]
+      }
+     ],
+     "prompt_number": 24
+    },
+    {
+     "cell_type": "heading",
+     "level": 2,
+     "metadata": {},
+     "source": [
+      "Step4: Analytic solution"
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "We have an analytic solution when we have a sphere in a whole-space. simpegPF provides this function so that you can compute magnetic field on your receiving locations. Another input you need to put is direction of the earth magnetic fieds, and the strength, you can easily get this information from <a href=\"http://www.ngdc.noaa.gov/geomag-web/\">NOAA's website</a>. We assume that we have veritcal earth fields and the strength is 1. "
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "Bxra, Byra, Bzra = MagSphereAnalFunA(X, Y, Z, 80., 0., 0., -100, chiblk, np.array([0., 0., 1.]), flag)\n",
+      "Bxra = np.reshape(Bxra, (size(xr), size(yr)), order='F')\n",
+      "Byra = np.reshape(Byra, (size(xr), size(yr)), order='F')\n",
+      "Bzra = np.reshape(Bzra, (size(xr), size(yr)), order='F')"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 25
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "fig, ax = plt.subplots(1,3, figsize=(18, 4))\n",
+      "dat0=ax[0].contourf(X, Y, Bxra, 30); ax[0].set_title('Bx'); cb0 = plt.colorbar(dat0, ax=ax[0])\n",
+      "dat1=ax[1].contourf(X, Y, Byra, 30); ax[1].set_title('By'); cb1 = plt.colorbar(dat0, ax=ax[1])\n",
+      "dat2=ax[2].contourf(X, Y, Bzra, 30); ax[2].set_title('Bz'); cb2 = plt.colorbar(dat0, ax=ax[2])\n",
+      "for i in range(3):\n",
+      "    ax[i].plot(X.flatten(), Y.flatten(), 'k.', ms=3)"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "metadata": {},
+       "output_type": "display_data",
+       "png": "iVBORw0KGgoAAAANSUhEUgAABEUAAAESCAYAAAAbqGtQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvX+wZVV17/s9pDndwRdpW0W90Nz+AQasWHUBqcRVVtGV\nRq3wiI0POoHyn+Q9GtuAeIXLjyS3uloqiQ1GuGITmh+5puo+H15bVCweRhFfd5W1MUHw5WlFlB/d\n2hJFQ9ONkXQ3be/3x1rznLnnnj/G/LXWXGuPT9WuffaeY8019zp7jzXXd40xJsAwDMMwDMMwDMMw\nDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMw\nDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMNI3AzgmOHxLQAXdzc0hmGYmYF9McMwTLewH2aYwljS9QCY\nmWMngGel18sB/GHz/g0APtbFoBiGYWYM9sUMwzDdwn6YYRhmxhCq+O9q2k4EsL9pZxiGYfLBvphh\nGKZb2A8zTGEc1/UAGAbAQQCfbf7+T10OhGEYZoZhX8wwDNMt7IcZpgNYFGFKYU55vQa1Sv5Z5f2z\nm/fvbGNQDMMwM4bqi69H7XPP0tg+jvqOJsMwDJMO1Q8D5hokxwBc197QGIZhmBhEqOB6TdtyAC8C\neEp5/7pmG7ng1DMaO4ZhGIaGry8+sbHfodgK4ZoFaoZhGD9C5sT/m/K4GMDTAH4FfRoOwzAecKFV\npm3eD+Bd0uvlAP4A9clho2L7MdQFp+4B8DUAfwZgFYBzso+SYRhm2FB98UEATzRtm6X3L2me78o4\nRoZhmCHjMyf+vPL6CtTi9M0Avp5rgAzDMExabMuP/Qp1CLZuCbLVjc3jzfN/aWOwDMMwAyXEF4uo\nPTmFhqP2GIZhwgidEwtEKvlX8g6TYRiGSY2t0vZqAF9t2jdp2jc1bTwBZxiGiSPEF6spNGJCziI1\nwzCMPzFzYlN6DcMwDNMDbCcAYHHS/S3LtvsbO4ZhGCaMUF/8LQAvNH/fhfpu5qtzDJBhGGbgxMyJ\nn0Hti1dlGRnDzCi8+gxTCgebZ1X0OB916PZdqNXxnW0OimEYZsYw+eK7ALwGdWHAP0Bd5+mlFsfF\nMAwzK5j88E7UYshGAHtbHA/DDB4WRZhSEEX7npDeEyLI4wA+AOAG1CKJLc+SYRiGCUfni4HF5dHv\nQj1R5wKrDMMwedD54etRz39vABdWZZjk5Fx9RqyZfS6Ax1CvJKK2PwtgRfP6Hs92pp9sBvBu5b01\nqB39GLWzF9yLOjxbVOBWV6M5CIZhbLAfZkz4+GKg9rdfQy1Mv4jplRAYhjHDvpjRQfXD5wPYhvo7\nsAeLoongRQCP5BsmwzChbFNefwuLJwSgzqWT8+i2YfLuv6ud6R/bsFhVW620/QKA/4nJ/MhLoC/k\nJ1aj4YrbDGOH/TCjw9cXy4ii13dmHyXDDAf2xYyKrx/eZLHnOTHDFMqJmHT2QP1j3i+93q+0r0dd\naZnazjAMw5hhP8zk4HrUE/D/1PVAGKYnsC9mGIaZUdagnjStkt4Td/2Beik/1cGf7dHOMAzD2GE/\nzOTgGfAykAzjA/tihmGYHpCj0OqzqB32Xum9dwJ4uPl7BaYd/IHm+dWEdoZhGMYO+2EmJTej/u6s\nbv5mGIYG+2KGYRgGQL2CyH4squSXYNrBL8eiku5qZxiGYfxgP8zE8DTqPPePdj0Qhuk57IsZhmEK\nJOfqM4LPoi4Qtbd5fUBjI6pp7ye0L3Dc0l8bHzv8qwRDZBimBQ4B+PWQDU8Axi/TzV/Eos9garL5\nYQD4taXHjX91mKO5Z4AbmwfTb9gXd0dGX/xr47oOJ8MwPYD9cGHkFkW2NY//V3pvP2qVW0a8fonQ\nvsCxw7/CB8e3WAewDys9hkvje1s/hzO3qqth5SXVPldin5f9P2z9Kn5767ui9+vD97Z+Du/a+tut\n7hMAvrr1H5z79T1+Lj639Xu4ZOuZSfvser+m39z1c59cFtrnywBuJ9peDbwmdD8DJasfBoBfHT6G\nW8YfBJD+N2IixXc4ZKx/u/Vf8H9s/Q9R+5VZNRFVr+fjW1/GtVtPSLZPKl3sl7LPvRlukIf8X2Pn\nF235f3Wc7Is7I7Mv/hXo/52UPATgghnYZ1f75c86zP1ezX64MHLUFBFcjLo69teb12c1z09gWvle\ngcX8Slc7iX1YmUUQ6Tt8XPxZuXDU2rnY6zt8rIqiUz9cMl1/R1dhL0kQYaYp5dh1/R2i0pdxDhz2\nxQzDMAWTSxQ5H7XTfhy1or0GwB9K7Xdjco318wHc5dFuhC/6afBxosGTyXBYTOqc1v1wX/7XXY6z\nlAv6IVDCsezLd57plM7mxAzDMAyNHOkzy7G4frrstHdKf9+Iet32i1GfHJ4G8HmP9gnavrh/3bq3\ntLq/nPsUx840sTt53dos+7Wxdt3Jre9Tt982JrtvWfe67PsoYb984dA6rfvhtgn9Dsd+F89a9xtB\n28VcvL993fHB28bQxX5D9imObUxaTej/Fai/UyHzkK78P9Mqg/fFwOkzss+u9sufdbj7ZUpirusB\nRDK+aHxf12MYDKVcuHY9jq73P1Qum/siEO5zxh75kzH7YcIY3ze+qOsxOOnit911JMOskqPuiIvS\noy/F+K6f+yTAvniIjLupKcIwjD9XA+yHi6KN1WeYnrAPK2deEJj1z88wQ6Xt33apYsjJP38heZ/P\nvf61yfuMZRX2ti6MhEaMMAzDMAzTLSyKMBPMsjAyq5+bYYbOrAoiOQQQn/10LZakSKvxpWRhpOSx\nMQzDMEyXsCjCTDGLwsisfd4+0v7CxcwQaPO33bUY0pYIQkUdT1ciSRdRI0OGfTHDxPLWjH1/J2Pf\nTCmwH04PiyKMFlcB1iExC5+RYWaRIQsipQkgFLoUSdoURkqOyODzHcO0SU7xw2efLJQwjIuZEUV2\nV1sAAOeNbkpun7Pvru11USM7q+0AgI2jq5x9+9gK+3kcwVWjjST77VVdwD3U3jVB3FLtBgDcNDrP\n2beP7azZC9uOuQ7As6iXRgSAeyLtY9uXo15V4LHG5lsAvh0x3uJp8ztJufjbXD0JANgxOoPUv87e\nJIZsqA4CAB4YnUjqm2ovhIXqgvr16CFS91ntw/uuP4tLHEl1LE3pNCm+BypCGCnJD8v2HcO+uHNu\nbZ6vyWCfs++S7T9FsL2yeb6D2Hcqe5NQwt+DbuxvdZvkh/2wwnGpO2SGR6l3vGJYiX18x2y2uBnA\n4wDuR+1I16Je3jDUPrZ9OYCvoT4B3N+8/tOI8TISbf2224oOOfnnLyw8hkjbn6+t/xufY7SwL2YG\nwFulxwnNo2/IY5c/DzMDsB/W0PclenhJ3hZpa4KXez88Ue2G2CV5HyYavrN+UvezH4vqMgCsB3AD\ngHcZunHZx7bfhVoNv1eyORHAwcDxdk0xS/IOSRAZqghCJXd6TVvpNKXdWGBfPFhfzEvyZmHWhQJO\nvclD3JK87IfT++GZSZ9h4hlCAda+j58J4mzNey8COD/QPrYdADYB+KhiI5y/73iZlskpiMy6ECIj\njkUucaSt1WlKrjHSMuyLmR4w6yKIino8WCTpOeyHDbAownjRZ2Gkr+NmolmBWmWWOdA8vxrAS572\nse2va/5eC+Ccxn45gI8FjpdpyP0bZzGkG9oQR3h1mlZgX8wUCgshdORjxQJJD2E/bIBrijDe8B0v\npmcsx2TYHbDoYNX3Kfax7Wuav8dYzI8EgG2B42XQX0FkyLVCUpPzWOVOh2JRHgD7YqY4uI5GHFyL\npIewHzbAkSJMEH2LGOnTWBk9pjXZRwAetW96QPOecKSq+kyxj20X+/yW1P5I8/rGgPEymckpiLTK\ns5n7X+M2SUGuyJHcESNDSaNhX2wcL9ML+AI+D+K4cvRIG7AfNo43GBZFmMHDgsiwqZqH4LZpk/2o\nlWYZ8VoXdueyj20/IP0tkEMBfcc78+T8jecQRFoTQ3KLILb9tSCQnPzzF3onjAwZ9sVM2bAY0g6c\nXtMl7IfD4fQZJpgh3PFiZoInMK00rwBgKt7tso9tf7ZpXy21yw7ed7wzDQsiEs8qjy5paSw5Umpy\nptLMuEjPvpjpAE7v6A4+9gXCftgAiyJMFKULIzM+AWUWuRuTa5qfj3oJMMEapd1lH9v+UUxWzv5D\nANd7bM9kJvWFcZZ6GCWJIC4yj5WFkd7AvphpCb4gLwf+PxQG+2ENoesjl8L4ovF9XY9h5kk9wUvV\nX2kTz1TjKV2IMnHZ3BeBiDXZf0w0PKV+0u3nOtSXY2tQL+clr4e+CcAlAN5NtE/VLhgD+GvP7Uti\nfN/4otZ3mus3nkMQSUrpAogPGdJsUqbU5Eql6dKPsy8erC8eA7d3PYYC4AvwsuGUmpqrAfbDarug\nEz+cUxS5BMDbUBdJUd9fA2An6g+1CcDnAOyRbMQHF4VU7oEesiiyu9oCADhvdFNy+5C+j2Apzhjt\nINk/WW0GgCl70wVBzs9qsjeNZWe1HQCwcXQVqe+d1XbM4wiuGm0k2W+vdgKA1l43pi3VbgDATaPz\nSP372G+pdmMpjmDH6AxS35urJwEg2t40uc75WX3tt1S78dSjLwLdngBmkTb8MEAURVJ+x3S/7xS/\nKZsgsqE6CAB4YHSis+8N1UHMv3IUo4dIQ0F1Qf1stZfEkOryxp44Lchpn6xvg0BCOjaK/ZHjl5D+\nT4D7/6oKIz7fM5utzne34bfZF3dCK3Niuihya/N8TQb7nH3b7E1iyJXN8x3E/nPalzSW3PYuW1Uc\n6ep704X9rUA912A/XBA5Cq2uB3A2gHcCeEbTvgL1MjvbUOcIXY5J538zgK8A+HrzehvqkJLine truncated
+       "text": [
+        "<matplotlib.figure.Figure at 0x56f2150>"
+       ]
+      }
+     ],
+     "prompt_number": 26
+    },
+    {
+     "cell_type": "heading",
+     "level": 2,
+     "metadata": {},
+     "source": [
+      "Step5: Solve PDE"
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "Note that data for this case magnetic fields projected to earth field direction, which is typical data type for airborne mag survey. "
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "First, set survey class"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "survey = BaseMag.BaseMagSurvey() # survey class for mag problem\n",
+      "Inc = 90.\n",
+      "Dec = 0.\n",
+      "Btot = 1\n",
+      "survey.setBackgroundField(Inc, Dec, Btot) # set inclination, declination, and strength of magnetic field\n",
+      "rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)]\n",
+      "survey.rxLoc = rxLoc # set receiver locations"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 27
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "Second, set problem class then pair with survey"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "prob = MagneticsDiffSecondary(mesh)\n",
+      "prob.pair(survey)"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 28
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "Third, run forward modeling"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "data = survey.dpred(mu)"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 29
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "Now we viualize computed solution and compare with analytic solutions! Note that you may always need to make sure that your numerical solution is reasonable enough. "
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "fig, ax = plt.subplots(1,3, figsize=(18, 4))\n",
+      "vmin = Bzra.min()\n",
+      "vmax = Bzra.max()\n",
+      "residual = data.reshape((xr.size, yr.size), order='F')-Bzra\n",
+      "dat0=ax[0].contourf(X, Y, Bzra, 30, vmin=vmin, vmax=vmax)\n",
+      "dat1=ax[1].contourf(X, Y, data.reshape((xr.size, yr.size), order='F'), 30, vmin=vmin, vmax=vmax)\n",
+      "dat2=ax[2].contourf(X, Y, residual, 30)\n",
+      "cb0 = plt.colorbar(dat0, ax=ax[0])\n",
+      "cb1 = plt.colorbar(dat0, ax=ax[1])\n",
+      "cb2 = plt.colorbar(dat2, ax=ax[2])\n",
+      "ax[0].set_title('Bz (analytic)')\n",
+      "ax[1].set_title('Bz (simpegPF)')\n",
+      "ax[2].set_title('Residual')"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "metadata": {},
+       "output_type": "pyout",
+       "prompt_number": 30,
+       "text": [
+        "<matplotlib.text.Text at 0x4858f10>"
+       ]
+      },
+      {
+       "metadata": {},
+       "output_type": "display_data",
+       "png": "iVBORw0KGgoAAAANSUhEUgAABFQAAAETCAYAAAAYmFIjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvVuMZdd55/cjk4GTzIQsVoMDCE0m1dUkIgUeZbpFIdG8\nhHY3JSQBAiRsSsnLPKlI+SEJIrtFKg9B6cXqFmULmQQwb28BohGp1kOABBC7KdMvkWPePJ7A0oB9\nqQHVEGClb0Q8jjMwKw9rrzprr7Pue+3r+X7AQZ2z1zr77HOq6r+/89/f9y0QBEEQBEEQBEEQBEEQ\nBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQ\nBEEQBEEQBEEQBEEQBEEQBEEQBEEQhCMuAp94bu8CT/f42u8Bv9nj/nO5DNyuuK+rmc/Zal7/RKVj\nEARhPOamrSWaNXXO4v8dXANeAh5MnG/eQPRaEARBEISO/KtjH4BQlTeA68bjLeArzfbngRcrv943\nUIHpTyrvtyuHBc95FhWYnwb+1NhP7r7uAt8GXga+WHAcgiBMj7loa4lmzYX3gCvG4y3gKZR2nwUe\ni8x3IXotCILJaZRZ7uM6Si9qa76Lyyiz19Y2m12Ukd7Hucg+ns8B2z2+hiAIwmjoq6iuq5kPoq7A\nfeIY64rvNcfkMnCr4HnPot7P369wDFvNvs5U2JcgCOMh2jo+OuPk257x15vxPWv+7yTuX/RaEATN\naZQevIPSHPP2B6w0//UBjuVN4MOEebvkaV4ppfG1ICye+8c+AKF37rES/hpmgeZZ4A7Ty07pyn0V\n9nEX+CHqaoEgCMtEtHUaaKPlbOHzRa8FQbC5AnzTuv0WKjvjOnAOONXzMXwReLzn1xAEoQJiqGwG\nLpMgVFt+PmGfz+N26M+iXOw7zb6usl7jDipV/t1m+xusXP83cZ+kUvdrotPmXft7j1WvlcvNvvT2\nD43tdj+C3eZ47zTP9x3vD5pjfiBwfIIgzJshtfUbqJ4hn6C053XWe3/YmuXS2avAhWb8LErztKba\nx/dJM1frr37+S7jZQqXD6+N07VNz2trnBVZXWv9zz3Nc3DNeuxTRa0EQUvlh81N6LwmCAIihsjRc\nwf0W8GVUgPunxvZz1u0ZlOt+iAqwQ2yhTiSXre3nUAbDDirgvtjs81lUMG9y2OxHB/PfAL6DCmrf\n6rBfk5ebn89Z23dRJsgPmsfPNq+t7z9jHafmNCrw/83mOF5p9vUe680p9Xt4KnB8giDMg7G19Q2U\n4XAVpZVXmn279ndo3dca9X83z73b/HwPZcp8H2XibKO01S59OYvS39soTdXaa5vNW8AN4KvN/G80\ncy+y3pPgbLPtdLPPHxLXcx+nm5/Xg7PCiF4LgpCKL5M5x1AuMchhZW6bFxZdvOd4LrjL20suWAqC\nICwOXef/enNf315GCeQt4inpWmR9deom55q5O9b2N4C/Yf0q3xus9xnQ2/7A2n6B9Xr2nP3aNZ7v\nsr7qj85cMT8T10nmMu361WvNvs3j0H0U3mGd26y/P0EQ5sNUtPUTVgawZs/x+rZmaY38qrHtwWbb\n31jPPcV6Hb7OrDGfDyud3jO2vdEcj33s5x1zQ1pqZ6jonigXWOcs6vdgvhc9/13avzPz5vqiIHot\nCILuoRLS62usx6RbrLToD1A6+iYrLTLRuvzjZp7uA2XHqrae63PJreY1LjTP0caMqd3v4e6/Yse6\n+pzzIateMfq433Qcj/RQEQRhsYSW9vwblOCGlvfUJ5AfJ76eNiRSea85DhN9QrFNEh0MpzRkdO3X\nFnwdzJtlOddYP9HEDJVQ068zrH/hAHUStU9IgiDMh6loqw62YzX7LkPF1khQGmibwLo5q/lFQgfv\nLuwvCi6D3BzThpD+TFxaqo0al6ES+j18NXP+juO1Ra8FQdD65DJkTSPdjlFzDOVSg1y/tm1EX6Xc\nUOlywVIQhAZZNnlZnGW9keEJ1EngDVTpy6vW+BYq3fka8KXE1zkZGNttjuMk6sSkM01cS3keAh8n\nvmbOfk1eQZ0InwO+1jz3BOqLSw6htHK7RElzB3XcgiDMm7G19Znmdd5DadAVVHB7KXG/LuyroT58\nS4h+wCpFXevcc6yXWMKqxNOc+37Ga4F7GeRrzbYDx/xvAN8N7M9G9FoQBM1pVnGfzRXWyy2fRp0P\nDqztL6Ji0KdonyOeQhnkHzSPX2X9HGJyFmWePEs7br7Hyuwp4RnP9l3i8bUgCA3SQ2X53GAlmK5A\nV/cwqVE7/g2UU/5tlNv9A5Qov0K31XO67PceKnA3P4PD5rkl3M2cLyckQVgmQ2rrJeCh5vWuoIJr\n3WR2So0RX2b1RcS8PcHqM9ou3Ldr1Y3XcJsppYheC4IAyqS433F7jpX+akxD2ZUVd0i7zPAZVj0E\ndb+SUKYjrMwdl+kc680VYxdl1Fxk1Z+l7xWMBGFRSIbKZqBXQXCttLODCvgPMvZ3zbFtC5WufZn1\nq7FdzJQa+325uZ1BNZG8QnpmjEZnpnyO9SvVz6JOsF+2tuvl9QRBWCZDaOuDqMyVayhjRWelnEHp\n4kXWtacmvqyNXVZND7XO3Ue7Qa/mWVSJ0QGrUiOXln6++Ci7I3otCEKMV1F661qmXceaNvfRvhin\nDfKzqHPEWZRG3kXp4g3HPnIv5qXyDVSMfQfVy+UHqPKjF5pjEgQhAclQ2QzONT/NFOtvoBzx51kP\namPooHPH2KavOn7QnsoW6mRReuWvxn71EqQvo76chFIjfSbN+6j3/U3aX562UF9oXFeJd5EAXRCW\nzBDaepJVPb+JvirZd1bFSdr1/7Bq7Gr2AbiCCsDtK5sXUVdgH2oef4D6cuDS0jEDeNFrQRBS0Drx\ngPVYG8r27fOsSh4fZJVtcglViv4YyljR8WToNV2m8xMZx/454755wfIY8Fussv66XAgVhI2jzwwV\nvVTY51FXpF50jF9n9YXZrh2MjQvrfI31LI5dVHB/iArwQRkRF1Cf7w1WXwo0d/D3BYFVHftpVldf\nrzc33Zvk3ea1X0Clpd+H+p2+wuqqbopg19jvPVZp8neAHzleRzfa+ibqS4K+Cmzu6znUiedG83q3\nm20PsPpsNVuoE6e9/KkgDInocB3G1FZt5j6L+j28izI5vty8tm0Q2/rXNTC+27zGU81rP4XKjrlG\nu0fJcyiT5z3UMsjm3Ddom0t7rHrCXDK2vYP6DPu6GutD9FroG9Hi5bHNKttZG8ov074AeBH1u9MZ\nLdogfwV1XtHEDPIrKF28iLpIqOPdLdbjT3A3Md9idd7Qxw/1L4QKglAJe3nDd2mvxX6RdofsC7Tr\nB2PjQhu9MsLfsF67eQtlEOwY8/cC81NXpNB1nyYnUEJ/u7n9GNVJ/EQz31xd4XXcq0+cbbabv/+c\n/b6Juwu5fs++VSgeZLVUnE5Jf5P1Lumnmu36WN7BvSKRXorO7pwuCEMhOtydqWjrg802vazwLVY6\naGJrlk9nrzqOxbfKzw9QfwfvslpeM6Sjr7NadeJD3Kv5wEpLteZ+FfX3Za+2plftSVl22pzve10X\notdCn4gWz4eUZZP1amzm7+AEK31+vZmj+5HYK/pctebp1YPs2Pcyq7JKWOnjbdTfwEWUsa212dQ8\nvbrQuyij5xvN87Tmao3Vx3IBpYN6nt5+nlUm4WXSm5kLgtCRB2mfKEAFmeY/of0PeYb2coWxcWF8\n7N/p1NEnQPsLSF+8QfpSqYJQG9Hh+TIlbXV9GajBadb7zsBKp3d6eM0QotdCX4gWz4sUQ+UMK8PY\nJNVQLjXI9Wvri3qmue0ykS8Yr/EhyrQ+hTJu9GvVuGApCEIP7LIeEOmrP6DEyj45nM4YF6bDJ6yW\nL54611g/MfXJJ7gzVwRhCESH581UtLUvQ+UO7i+E7zGsTmtEr4W+EC0WBEFYOH00pb1Ou/4bVB21\nrk3eZv3koOulH0gYF6bD8/gbaE0FvQzcCYY7Vp3qmduQUhBqITo8b+agrV34XVRpzuuoL5fPosyU\nv4+7H0CfiF4LfSJaLAiCIHRmC3Uy2Gken2P95KBrt3cSxoVp8S7TvrJ3FZWimFp/3xX7710QpoDo\n8PyYgrb2laECquxB1/7rVPOh36/otTA0osWCIAhCNm+y3mDOPjnolMgHEsYFQRCEPESHBUEQxke0\nWBAEYWH0uWwyqIZIF1DrsGtus1qPXaMff5wwbvDQoSrFFgRhBlwDHit54r8Oh3+VPv0Oq+UAhd51\nGESLBWFWiBaPQ89a/G8fwj+vcZyCIPRPsQ7/G3D4L9Kn+3Q4dyn2rku711ga/hzwBPCCtX0LlfV6\nF7U0Odaccygj+g3U57EH/BC44XiNIvo0VJ5GOfG6LvkUaq3z91nVf2q2WdWTxsYN7gD7NY41kz8E\nfmMDXnOs15X3uszX3T8Zn+Pmr4DfS5z72/BQ6esskAF0GMbR4k3635H3uszXHeu9ihaPwABa/M+B\nj6ocbB6/D3x9A15zrNeV97rM1320WIf/BfCPEuf+124dvogqs9V6pJdiv+TZTWx+3+NnUL2onkIZ\nUTbfpN1/7V2UaaJNmW1WhvZd1IpX1cwU6KcpLagUxW1Uk7ktlCv0FWP8Fdrrt59FrcOeOi4IgiCE\nER0WBEEYH9FiQRCmxB7tRuyXgec6zO97/C3gRZTBfJ/j+J5GmSSa6yjzRXPISnu3gR859tGJPjJU\ntlgth2gK/hvG/RdQqT1Po97cVdpvLjYuCIIg+BEdFgRBGB/RYkEQpsRpx7Y7KKO2ZH7f4ymcpb2S\n2kng+9acj3GWrNehD0PlLmmZLy92HB+RnQ15zbFed4zXHOt1x3jNMV9XGAjR4UW97hivOdLine truncated
+       "text": [
+        "<matplotlib.figure.Figure at 0x6054210>"
+       ]
+      }
+     ],
+     "prompt_number": 30
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\" /></a><br />This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">Creative Commons Attribution 4.0 International License</a>."
+     ]
+    }
+   ],
+   "metadata": {}
+  }
+ ]
+}