mirror of
https://github.com/wassname/simpeg.git
synced 2026-10-01 12:31:30 +08:00
Change functions in the header and start a development branch.
This commit is contained in:
D Fournier committed
1 parent
80c839390c
commit
7d516d2879
1 file changed
+462
-433
@@ -1,443 +1,472 @@
|
||||
{
|
||||
"metadata": {
|
||||
"name": "",
|
||||
"signature": "sha256:550ca3df316c026a67f55f9de3baa38df81944304bfe8103f9a2d9e8a43c20c7"
|
||||
},
|
||||
"nbformat": 3,
|
||||
"nbformat_minor": 0,
|
||||
"worksheets": [
|
||||
"cells": [
|
||||
{
|
||||
"cells": [
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"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"
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Populating the interactive namespace from numpy and matplotlib\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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>."
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"WARNING: pylab import has clobbered these variables: ['linalg']\n",
|
||||
"`%matplotlib` prevents importing * from pylab and numpy\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {}
|
||||
"source": [
|
||||
"from SimPEG import *\n",
|
||||
"from simpegPF.MagAnalytics import spheremodel, MagSphereAnaFun, CongruousMagBC, MagSphereAnaFunA\n",
|
||||
"from simpegPF.Magnetics import MagneticsDiffSecondary, MagneticsDiffSecondaryInv, BaseMag\n",
|
||||
"%pylab inline"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import matplotlib\n",
|
||||
"matplotlib.rcParams.update({'font.size': 16, 'text.usetex': True})"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"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": "markdown",
|
||||
"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",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<matplotlib.text.Text at 0x15eb8588>"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAucAAAFeCAYAAAAmKrHoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3c9yHMeV7/Hf0WhNwvALEJDvfkRLs7oRE7JAe39BjeYB\nDMreGyK1q9U1aXk//OPlRMyYNLWXKEuhrSWKfgCKwgNcgyD3o3MXXU0VClWVlVWV3VmN7yeCIVae\nkz8kQSqZrC40zN0FAAAAYP1eW/cCAAAAACxwOAcAAAAyweEcAAAAyASHcwAAACATHM4BAACATHA4\nBwAAADLB4Rwbxcy2zOy5mf1gZm+29OyV9a/L6w/L688C2Y/Lvl9MtNZHZnY8RRYAzNXQPbj8eZ8f\nFyZaJ3s2VsJ4n3NsGjN7V9IjSc/c/WcN9eeSLkh6w92PyrFvJF2W9J67P2yYc03SbUmP3P1XE63z\nkaTL7v7TKfIAYK6G7MFmdlNS2yHmDUlXJbm7/9NEa2TPxkpwOMdGMrPbkq5J+oO736iM35F0IOma\nu/+pMv6mpMeSTtx9uyFveaD/ibu/nGiNbPQAoOn3YDP7TtIlLQ77n0y0RvZsrASPtWAjuftvJD2T\n9KGZ7UiLx1m0OJg/qh7My/4nkv4gaau8G/OKmd2SdFHS9akO5gCAH025B5c3YXYk3Z3qYA6sEodz\nbLL3yv8+KP97R4uXQN9rai7vsJ/o9IF+V9KhpO/c/Y+hD1g+33izfK79kZkdm9nT8k5+UGXe8rn5\np2Z228wu1voemNk3Znax/Pnx8pnNpmfty2fx75jZd5Xcwz5rAoBVmGgPvqbFTZjv3P23PfrZs5Ed\nDufYWJU7MZfLlzh3tHiJs+vOy/Lgfqfy39YDfYs9SZ9JOi7nP5N0zcyedk0ys6vlvEtaPFt5azlX\nP/4DY8klbWnxMvAPkj7U4te6J+mvtdwtSd9L+nWZ/2GZe6t8zhMAcjF4Dy4P8rfL/isRH5M9G3lx\nd37wY6N/SPpOi83w057998v+O+V//yPiY/1Q/vh1bfxmOX5QGXsk6R+V6weS/kfShdrcB5J+aBqr\nr63ycd6t9f5D0qVa72F9TfzgBz/4se4fQ/fgcq//H0n/J+JjsWfzI7sf3DnHRivvpOyUl2/0nHag\nxUurB5KeS7oe+WGf+9ln2pdflPpB2yR3f8/d/8nP3tnfVfs7EtTX9vkyrjK2L+m+l+9MU/l4H5c/\njbnDBACpRe/BZvZAi73+Dx7/nDl7NrLC4Ryb7oEWm94tSbt9niN09xeSfl9e/r5h4w1pe9nxiX78\nh0IjM9s1s2tmdqt8jvEHSY3v175Yavfayn+cSNIHTe//q8Xn5mJHBACsVOweXD5nvi/psbt/NOBD\nsmcjK6+vewFAKuVX+L8p6Za7f2RmP9fiOcIH7v7XwPQX5X9Pki6ywsw+1OIlzudavKz7Zy3uHN3Q\n4hnGMe7ox2c4T31YrfDXCAA99dqDzeyyFs97P5f0bupF1T42ezaS4HCOjVRu2Muv8F/eSflAi2cS\nH5jZTnl3JoXdjvHGLzAqvwDophq+yZGZ2dCFuPuzcrq5+98bPu41SV8PzQeAdSn3zb+q/ILRAa9y\nLrFnIys81oJNtXyc5dVX+Lv791o877cl6V7Cj/2GmR1UByrv0/vnljnLb7rxpDZvS4uv5h/z3cI+\n1+IVg1MvtZZrui3pJyOyAWBdHmixr/7B3b8YkcOejaxw5xwbp9zAdrR4nOXUnQd3/9jM3pd01cz2\nveHbRE/gRNIdM7uixbOMV7R4ubXpfXqtXNczM1t+0ySV83a1eHn0ePHLskMtvqnGi+rcHj7Q4u27\nHpvZX2prejDyLzUAWLnyDvK7P17arY72/xt4pZQ9G1nhcI6N0vI4S917WjzectfMHrW8FOoafufj\nkRbPCv5B0tXyY93xs98Qo/4xrmjxhavXyh9fS/qFFs9ePtLiJdQH5XXX+k6Nu/v35Tf0uCfpcmVN\nHzb8xQMAOQjtwRcrfR8Gcv5DPz7D3oQ9G1kx9zGvvACoKr+a/oG7v7/utQAAurFnI0c8cw4AAABk\ngsM5AAAAkAkO5wAAAEAmeOa8wsz4ZACYLXcf/P7Kc8SeDWDO2vZs3q3ljGLE3C8lvdOROyY7Zv6X\nkr4K9Lfl1cdD1109y8/HmAzWwTqGftyu8Wr9X9X+/21s3pj5XftH3+zzqBgx9zzs2X16ltdt/0/G\nZLAO1jF0HV3jsT1DemMz0u3ZPNYCAAAAZILDOQAAAJAJDueTurTuBZQurXsBpUvrXkDp0roXULq0\n7gWULq17AaVL615A6dK6F1C6tO4FnEOX1r2A0qV1L6B0ad0LKF1a9wJKl9a9gNKldS+gdGndC8jM\npWTJHM4ntbPuBZRYx2ms4zTWcRrrOL9y+ZyzjtNYx2msI0/pPh+vJ0sGAECSmd1299/Uxg4lPZO0\nLUnufi+mDgCbijvnAIBkzOyWpLcaxh67+8Py0P2Gme33rQPAJpvd4dzMbjeMHZrZvpkdmNlBbB0A\nMD0z25XU9F7kB+7+ReX6kaQPIuoAsLFmdTjnDgwAzMq7WhysXzGzyw19zyXt9akDwKabzeGcOzAA\nMB9m9q6k+5Lq3wFvW9JxbeyknHOhRx0ANtpsDufiDgwAzMmWu79oGlf5RZ4Vy8P4do86AGy0WRzO\nuQMDAPNhZvvu/rClfNIwtjx0H/eoA8BGm8tbKW65+wuz+tl89B2Yl1MuEgDOOzPbUfMBe+lYi725\nakuS3P2lmXXWz8Z9Wfn5JfFezADy9L2ko16d2R/OuQMDALNyWdJu5bHCtyVtmdnvJD1092/NrL43\nb6t8bDFUP+udqdYNAAnt6PTNg69aO7M+nK/+DozEXRgA89D/Lswq1W+mmNk1Sbvu/sfK8N3ajZc9\nSXci6gCwsbI+nGvld2Ak7sIAmIf+d2HWpfy+Elcl7ZT79j13f+HuN5bff0LSrqSn7v7Jcl6oDgCb\nLOvDOXdgAGC+yu8tca+l9nFgbmcdADbVLN6tRTp7B8bMLkqLOyxa3F3fN7NDNdyB6aoDAAAAucj6\nznkVd2AAAACw6WZz5xwAAADYdLO5c746RcbZsfND/W31+njoesgcMshInREa71sf2z/1fADAJjN3\nX/casmFmnu4vzkLjsmPnh/rb6vXx0PWQOWSQkTojNN63PrZ/6vnd2e5+5ju1bbLFng0A89S2Z7++\n6oUAADCdImHumOzY+aH+tnp9vKkv1EMGGTlkdI3H9gzpTZnRltuMZ84BAACATHA4BwAAADLB4RwA\nAADIBIdzAAAAIBMczgEAAIBMcDgHAAAAMsHhHAAAAMgEh3MAAAAgExzOAQAAgExwOAcAAAAyweEc\nAAAAyASHcwAAACATr697AfkpMs6OnR/qb6vXx0PXQ+aQQUbqjNB43/rY/qnnAwA2mbn7uteQDTPz\ndH9xFhqXHTs/1N9Wr4+HrofMIYOM1Bmh8b71sf1Tz+/OdndLFJ6lxZ4NAPPUtme/vuqFAAAwnSJh\n7pjs2Pmh/rZ6fbypL9RDBhk5ZHSNx/YM6U2Z0ZbbjMM5AGBSZrYl6UDSiaQ3JMndb9R6DiU9k7Rd\n1u/F1AFgU2V/OGeTB4DZ+cjdry8vzOwbMztY7r1mdkvSp+7+RXl908z23f1hnzoAbLI5vFvLR+7+\nsbvfKw/le2Z2sCyWm/hjd39YbvxvmNl+3zoAYHL7ZvbryvUzSVcq1wfLg3fpkaQPIuoAsLHmcDhn\nkweAedlz9z9Vrt+Q9DdJMrPLDf3PJe31qQPApsv+sRYtNvmjyvUbkv5LYpMHgBxV9+xyH/7B3f9Y\nDm1LOq5NOSl7L4Tq7v4yxZoBIBfZH87Z5AFgfszsoqR/k/SepGuV0pbKr/+pWO7T2z3q7NsANlr2\nh3OJTR4A5sbdX0i6J+memT02s9vl1/2cNLQv9+njHvWaLys/vyRpZ+CKASCl7yUd9eqcxeF8dZu8\nxEYPYB76b/SrZmZb7l7df29LuqPFPn6sxY2Tqi1JcveXZtZZP/vR3plm0QCQ1I5Onym/au3M/nC+\n2k1eOv3Jav/EDVOseH6ov61eHw9dD5lDBhmpM0Ljfetj+6eenzcz25P0Wbl3L/dZK2sX3P1bM6vf\nONnW4ov1FaoDwKYz93y/+/Fyk5f0apM3s2taHNC3lgdwd9+uzTl091+V15312sfzdH9xFhqXHTs/\n1N9Wr4+HrofMIYOM1Bmh8b71sf1Tz+/ObvtW0KtUPoZ4091/Wxl7oMXXC71fXt+U9HXlfc1vSvqb\nu3/Sp17JzfcvMAAIaNuzX1/1QiJ9LelO7S73FUkPKmN3a9+cYk+LO+vqWQcATMTdX5jZ3fKbv0nS\nTyU9dfePKj03zOyw/J4Tu2X9k77104pEv5JiZHbs/FB/W70+3tQX6iGDjBwyusZje4b0psxoy22W\n9eF89Zs8AGAsd38i6Umg5+MxdQDYVFkfziU2eQAAAJwfc/gOoQAAAMC5wOEcAAAAyASHcwAAACAT\nWb+V4qrxtlwA5iyHt1JcJfZsAHM217dSXIMiYe6Y7Nj5of62en08dD1kDhlkpM4Ijfetj+2fen4o\n+zwqEuaOyY6dH+pvq9fHm/pCPWSQkUNG13hsz5DelBltuc14rAUAAADIBIdzAAAAIBMczgEAAIBM\nvL7uBeSnyDg7dn6ov61eHw9dD5lDBhmpM0Ljfetj+6eeDwDYZLxbS8XiK/+LROmFxmXHzg/1t9Xr\n46HrIXPIICN1Rmi8b31s/9Tzu7N5txYAmA/erQUAsIGKhLljsmPnh/rb6vXxpr5QDxlk5JDRNR7b\nM6Q3ZUZbbjOeOQcAAAAywZ3zM4qMs2Pnh/rb6vXx0PWQOWSQkTojNN63PrZ/6vkAgLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x15cca3c8>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"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')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"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",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"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"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"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"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAv0AAAGjCAYAAAC/qqoOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3U2THMeV7vnnkOgtUV36AkBRsxchsVdj1iay0No3QPF+\nAAKU9gJB7mIxMwJE7QUQXF6z2wIa3JOgqNZWJEitRyBYdz8qvGybxJlFeAJRURHh4fEeWf+fWRoQ\nfo6f9MrK8vSK9IwydxcAAACA7fXS3AMAAAAAMC4W/QAAAMCWY9EPAAAAbDkW/QAAAMCWY9EPAAAA\nbDkW/QAAAMCWY9F/wpjZjpk9MrNnZvZaTc5+iH8Zjt8Lx59Fat8PeW8MNNZ7ZnY4RC20Y2Z3zOx2\nJOe6mX1rZodm9lnd8yjkXg65z8zsgZldSRjLXhjPo1j/EK+7vdL2PoG16zpfR36GBv95Yn6fRsp8\n3bd/3/vC+Izr9J88ZvampHuSHrr7jyvijyS9IulVdz8IbV9JOifpLXe/W9HnsqQbku65+y8GGuc9\nSefc/UdD1Jtb4TE65+5/m3s8ZWa2J+mBpDvu/nZNzn1Jr0m6I+mhpLck7Uk67+5/KuXelHRJ0n1J\nn0s6r/w59Dt3f7/FWO5LOh3u60tJ/yZpX9Ln7v5vhdwdSYch/6tyLXf/dexrB7ZJl/nazK5JqlsQ\nvCrpoiR395cHGiPz+8hS5uu+/fveFybi7txO4E355PRM0rVS+83Q/k6p/bXQflhT75GkHyS9MuAY\n70n6x9yP1YBfz+XwGP5k7rGUxnVR0vXwPXwm6Y+R8f+m0HZa+YL7sJS7F3I/LbV/FtpPR8a0eR7+\ne6n9Wmi/UGg7V5XLjdtJvQ09X0v6NuQP9jPG/D7ZeKLzdd/+fe+L23Q3tvecUO7+K+W/jb9nZmel\nfFuP8jOz99z941L+N5J+J2knnBF6zsyuK/8Bv+ruT6cY/8rZ3AMouS3pN8rf3WnyrvIzfb/fNLj7\nE0kfKX9evFnIvVroU3Rd+dnEN9Xsl5Ieufsnpfbfhn/PF9r2wr8PIzWBE2HI+Tq8Y3dW0kcVP484\nbinze8p83bd/3/vCRFj0n2xvhX/vhH9vKl+QvVWV7PmWjMc6+ovCnqQrkr4t/sDXCXtCr4XPDdwL\ne/8emNmNNgMu9Cvu875hZqdLeXfM7CszOx3+f7jZ51q1z9Dyzzrc7LL/vM2YwlvZm6/xvpn9PfIY\n1d1a74lvy91f8vwt+91I6p7ybTplm7azhbZfKn9OHJTu60/u/nKLxcM/lP8yUrZ5QS3e1+vh3++k\n59t9gBNtoPn6svITQd96i21yzO+SGuZ3yz+n9MxKn5sys3Oh/Q9txtRSynzdt3/f+8JEWPSfYIWz\nQefM7FvlP5hvRc7+bH4huFn4t/YXhRr7yrd5HIb+DyVdNrMHTZ3M7GLod0b5BHt901cvfnHZcEk7\nyvd5P5P0nvKvdV9SeS/ijvIF4zuh/nuh7vWwN3aIMV0O97/5f9PjdbF0eyvU9PD1jCV2huq08kVE\n2eYM+6ul3K+fF05ciLv7j2sWGZfDv/cKbZsz/R+b2TNJh2ERcGyxAJwwnefr8AvCjZB/PpJexPxe\nf+LsofJ3QS+a2YVC6I5a/mKVIGW+7tu/731hKnPvL+I2/035fs1j+68b8m+H/M2+6z8k3NczVX9m\nYLNX+1Kh7cieT+UT47F9qKH9WVVbeWyF+3mzlPsPSWdKuVfKY6r4elLG1GnPZ6Hfb0d+HuyoZk9/\nIXbse12OFY5vKH+B3XxW4JnyF9CzHce3eRz+UXy89eKF/0vl25TeKTxHH4z5mHHjtvRb1/laHfbx\nM7+3m9+VX3DgUPli+Xq4j8E+C5AyX/ft3/e+uE1740z/CRfO5mzeemv72/gl5b/VX1K+oLvanH7M\nIz/+mYHN1VzKe8CLOW95vjWk/E7EnuqvOlEe2+btxmL+BUm3/fhWlA/Df2vPcnUcU2tmdk4vrrLx\nQY86O8VbhxKxrT/FnM2/l5U//v+38ncsNmfi7qecgQ9viW/ePn+k/AW9+Hj/h6T33P11d/+9u3/s\n7r8M97dnZpfa3hewhZLnazO7o/x14Xeevo+f+T3uLeUL4i+U//Jx1Ye94k/KfN23f9/7woRY9OOO\n8snruvIFUnTvpecf0Nl8oPK3FRNiTN1bqt8osvcvLAAvW3494HthO0fdtYA9NrbwS48kvVu1h175\nY9O4QE0cU2thcf4n5W/7HrsMamTv/+b2k/A1HhZvZvZOylg8f1taqp68N22bnM21t13ST8NC/JPw\nwr95sWv1C4yZvaf8MqJvKH+uni2/OLr7h169P3nzHN1vc1/ANkqdr8M+/guS7nc80cD8HuHu3yn/\nBeg11Xy+ou38XlM/Zb7u1b/vfWFap+YeAOZj+VUcXpN03d0/MLOfKt97ecfj19V9Ev6t2sc3irAA\nvKb8bNVtSX9Ufvbqfb3Y693VTb3Y93rkbtXwNY48ps22lbozUediBTYLZMuvzFSu3UXVuwSbtn+E\n+3xsZlJ+Pf2D0njuhlj0RTOcbbygfJvBW6lnwtz9SbgvPtiLk67VfF14Z/GR4lfYGtQJnN83fyPn\nR2Z2OvxyVtR6fm8Qna8H7N/3vjABFv0nVJjcN1dx2JzNeVf5AuuOmZ2tmISGstfQXvlhr3DW+5oq\n/viXhZVdF+7+MHS3qgk0nPX6csoxhf53lH947Hx54VxQ136Mu3/RZzzB5g9slb1diG98Lanpj+7E\nFh/vKV/w1/6hsJB3Vvk2nv/w0h8hKmxj4iwTEFF4Z9EVv6BDE+b3iHAS5oryX0TeVf4u5r+V0g76\n3IfS5uu+/fveFybC9p6Ta7Ot5/lVBsJbjleV/3Z+a8T7frW8z9peXDv6jzV9Nm8TflPqt6N8+0af\n/ZWfK3+H48jZ5zCmG5L+eeAxNb5gFBa8VyOL9cPYre7t345uhvE9/96Fr/Wi8l8e/1bKPWel6zOH\nx1Q6evWdKh+EmrULfun5c/aC8itxlN+m39xX1Rk+AEfdUT4H/67nSQLm9wah/x3l26d+rfw1d790\nNR+p//yeMl/37d/3vjARzvSfQGGyO6t8W8+x/dFm9rbCJcXKZ08H8ljSTTM7r3z/53nlbyVX7W20\nMK6HZrb5Y2IK/faUv816mH9ZdkX5H5B5UuzbwrvKt7vcN7P/LI3pTt0LYIcxbd7i/MDM/lj12IYz\nQNeUn53+zvJLxhU9Kmy9GuLt39bC1pyvlX/vXlX+Nb6t8I5EKfeWmV2VdM/MPlL+9Wwe0yN//C18\njbeVP06/CvtwT0v6h+V/GKjKV+6++cX0XeUvOt+F+zLlL8qvSbrJCw7QLJzxfvPF4fNfzqv8P5F3\ngZnfa+b34GPlfwjxrXA/m9fcW2b2eeHr6zW/p8zX5Tk4tX9KLmbW9jI/3LbjpnwieSbp7w05Z1Vx\nacRSziXllxl7J/H+nyk/2/OG8kn0maS/q/pyX5/p6CXdziqfmDZnOj6V9JPQ/iCM50zIvS3ph4qa\n+yHvjVL76dDnQWFMv2nx9aSM6XT4mp5J+rLhcX0W+j2ruLW6rGrH50btJTtLj9ON8LVuLpP5RiT3\nQSH32GOq/Ez98/sN36Omx+CH8hiVv4AXvw9fKuFSg9y4bfMtNl/rxeUr637mij97Zxruh/m9eX6/\nGOK/qbifwef3tvN1eQ5O7Z+ay22+m4VvFjCJcOWDxn3aAID1YX4Hlo09/QAAAMCWY9EPAAAAbDkW\n/QAAAMCWY08/AAAAsOW4ZGeBmfEbEIDVcvdefzRobZizAazZ1HM2i/5jsh59/yzp5w11+9RO6f9n\nSX+J5NfVK7fHjptyNo9HnxqMg3F0vd+m9mL8X1X/c5tar0//pvmjbe2TKOvR9yTM2W1yNsd1P5Mp\nNRgH4+g6jqb21Jwuuak11jdns6cfAAAA2HIs+gEAAIAtx6J/UGfmHkBwZu4BBGfmHkBwZu4BBGfm\nHkBwZu4BBGfmHkBwZu4BBGfmHsAJdGbuAQRn5h5AcGbuAQRn5h5AcGbuAQRn5h5AcGbuASzMmbkH\nkOzU3APYLmfnHkDAOI5iHEcxjqMYx9jM7Ia7/6rUdkXSQ0m7kuTut1Liw1jKY844jmIcRzGOZVrf\n47G6Rf9yXzwAAGVmdl3SzyraPnX3L8LxNTO74O5328QBAOlWtb2n4cXjvrvfDYv5V83sQts4AGAc\nZrYnqeqympc2C/rgnqR3E+IAgESrWfTz4gEAq/Om8jn3OTM7V5H3SNJ+mzgAoJvVLPrFiwcArIaZ\nvSnptqTyH5/ZlXRYansc+rzSIg4A6GAVi35ePABgdXbc/UlVu8Lnqwo28/RuizgAoINVLPrFiwcA\nrEbkQ7ePK9o28/FhizgAoIPFL/p58QCA9TCzs6qeezcOlZ+QKdqRJHd/2iIOAOhg0Zfs7PviYWa8\neADAtM5J2it8pup1STtm9htJd939azMrz+u7Cp/ZisWP+3Ph/2e0xmtnAzgJvpN0MOsIFr3o1+Qv\nHhIvIADWYf4XkCrld2bN7LKkPXf/faH5o9K7uPuSbibEC34+yLgBYFxndXRN+ZfJR7DoRf/0Lx4S\nLyAA1mH+F5AYM7sk6aKks+FkzS13f+Lu75vZlfA3U/YkPXD3Tzb9YnEAQLpFL/qLePEAgHUJfxCx\n8i+gu/uHkb6NcQBAmtUs+nnxAAAAALpZ/NV7AAAAAPTDoh8AAADYcqvZ3jOdbMG1U/vH8uvi5fbY\ncZc+1KDG2DVi7W3jffOH7g8AQDpz97nHsBhm5uO9IGfqVzu1fyy/Ll5ujx136UMNaoxdI9beNt43\nf+j+zbXd3UYqvkj5nA0A6zT1nH1qyjsDAGBY2Yh1+9RO7R/Lr4uX26vyYjnUoMYSaLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x161cc710>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"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\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"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",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"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."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[<matplotlib.lines.Line2D at 0x1794ee10>]"
|
||||
]
|
||||
},
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFeCAYAAAC7EcWRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnU1sXEeS5//R3bduiNX0dYERqd67xbYH6AFmB7ZK6j03\n5TbmtsCakue2B9KSgMFODYydJlvWeSTR9xlLLd8lqm0YA/QAra++j76wwN7WJCXMwAbWVuzh5ROT\nWVkv831VZhX/P+CBlRGREVGPZFRWvnz5RFVBCCFk+vwgdQKEEHJUYQEmhJBEsAATQkgiWIAJISQR\nLMCEEJIIFmBCCEkEC/ARQUQGIrInIq9E5OQEm6HR3zPtj0z7TsD3A2P3bke57ojIbhe+SDUick5E\nnpjf32MR2ajRd1lEblp/VxP7G/2k41h372i2EK4DPjqIyCkAOwCequrPPPo9AMcAnFDV50Z2H8AK\ngPdU9ZanzzkAVwHsqOovO8pzB8CKqr7Rhb/UWOdoRVX/lDqfEhG5BmANwAMAdwGcRvG7/q2qXgz0\nXTb9FgDcBHAPwBkAQwB3VfWMZTsAsGvs77u+VPVvung/M4mq8jhCB4pC8ArApiO/ZuQfOPKTRr47\nwd8egO8BHOswxx0AX6c+Vx2+n3PmHL6ZOhcrp2WT021HfsfIFwL9y7+XXznyTSNftWQrPlseyimI\no4aqfgjgKYCPRGQJKKYeUIyEdlT1U8f+EYDfAhiIyKatE5EtFCOgC6r6chr5zziSOgGLC+bneUe+\nBUABnAr0/zWAPVX93JH/xvw8bcmWzc+ndZOce1J/AvCY/oGDUe19036CwCgWxVfIVwCWTLscQf1b\nZMxXKEZHQxQj3F0AjwFc9diOjYCtfnvG12MUo/kFx+4miq+55VfjMu87AE56Yg1QjOaeWH43It9T\nMCejf2UdE8+XY+ceUTnV+BvYi/3dTej/GMA/Tjifh0bWKIr661E1gEHq/4FcDo6AjyB6MKpdEZEn\nAJZQzPFWjWLfMz+vWT/VkscwRFEId03/pwDOicjjqk4ictb0O46iwG2VfVEUWRtFUQQeoPin/wjF\nex0C+L3jdwDgGYAPjP+PjN8tM/fdRU7nTPzyddX5Ousc7xmfat5PlywAeFg2zLmIRlV/pv6523Pm\n544lK0fAn4rIKwC7IrIrIldFZKFO3Lkj9ScAj3QHDkZ9tyPtbxj7cv5vbARU0bccyblzzOWc4Zol\nOzQCRlHQxkboRv7KJ3Nzs+Kccmy/BnDcsd1wc/K8nzo5NZoDtvr9puPfezlKvYriQ6ccwZffFJYa\n+i3z/do+Lzj4MLwHYB3FB175t/Q45f9A6iN5AjwS/eIPphCi/wlQjJp2ff9kEX1fYcKFNVjTIaYd\ndRHO/GN/78jKAuwWxqGRv+vE9X6IGN2NBufVl1PtAoyDC1dRH44Vfgb24fndf22K4q+sD6ldBC7C\nef6Wdix/bzr6DQDrnn5jH75H7UieAI9Ev3hTKFBcNHkFz1zshH7l6HDsHyrQb2IxMbnYI17fHPCy\nKWRbODyv6ivA33tiHCrAThHyHd+Hil+NnGoVYFMsJ87RBvIujzcnvMcPUHyQlu/xuON7FZ5VMhW5\nfmT5+gz1PpTLPD5L/f+Q6vgRyJHDrF44CWBLVS+JyM9RzMXeVNXfB7q/MD/3e03SQkQ+QjFa2kPx\n1fUzFKs2LuJgzrEp13Awr30oLCreY885lV/ZT0/Qr4QcqFlvbFa42NxX1ZciAhTrdZ87/W4Znfdm\nHRsRuYmiYD9BcQ2h1hpnVX1hYtWaf54nWICPGCKygmIU+0RVLxnxeRT/RDdFZElVX0x00I7lCrn3\nQpy5OLQJz40eYv57m6CqT0138RUOc/PEvWnmZPrfRHFh77RbHC0mycdQ1S8mqB4CqLrRpfID1nwA\nrQK4qarvV9gtobgI+c/q3MhjXfg7ssvTuAri6HETzuoFVX2GYl3oAMB2j7FPiMiaLbDWEn82oc+i\n+fnI6TdAMa3Q5lbOuyhG/odGeyanqwB+2nFOlcXZKmoXKgonUMzRVh4i8mZVLBSj/hVzd6Sdw5Z5\nuTPe5RCXUHyITyy+wOu/rVUUK0vcFQ9lLN83kCMBR8BHCPPPtYRi6uHQqE9VL4vI+wDOisiqO1rp\niH0A10TkNIq1uqdRLPh/oqqfuOmavJ6KSHnjCEy/ZRRf9XeLtyUbAK5bI/fYUeh5FF/3H4jI75yc\nbk4qgg1y+tp0vSQin/nOrZkq2EQxGnxmlrnZ7FnTQ9FTEBX6bRG5AGBHRK6buOV7P3RDjsnlhnk/\nH5rbkBcAfG1uZ/ZxX1XLD/PzKIrsMxNLUHxQnQRwre7UxVyRehKax3QOHFxVr7oRYAmBFQ4o5jm/\nh7OcLCL+KxSj3HdRFKxXAP4N/sX8d3D4otwSigJQjvBuo7jItIRi6uL1xSRjN+ki3PewVkEY+YLp\n89jKKXiBsWZOCzi4xfdexXktL2b5Lqq1Wg0xIeYCipF++d7v+d47Di7MfWady6pcv4dzYQ1FYbfP\n1z3w1mRuxkOmg1mAXzlfSMhRg3PAhBCSCBZgQghJBAswIYQkgnPAFiLCk0EI6QVVHVudw2VoY4wa\n9vsSwDsT/DX1Gds3ZOfTu7Km7fJ9t/GXMnbT3KYRe1IuKWKH4rf1HZKXur+C//+sjp+mfSf9j8f6\nHIdTEIQQkggWYEIISQQLcGccT51AIo4z9pGKnTr+fMVmAe6MpdQJJCLl+2bsoxd/vmKzABNCSCJY\ngAkhJBEztwxNRK5q8Wh1W7aBYjenRaDY6amOnhBCUjBTI2CzneJbHtkDVb1lCusJEVmN1RNCSCpm\npgCbPUh9d6qt6eF9W3dQ7D8aqyeEkCTMTAGG2SjaFpjH67jsodivNKgnhJCUzEQBNo9NuYHxJx0s\notjc2Wbf9DkWoSeEkGTMRAEGMFD/gyIHOHg+V0lZcBcj9IQQkozsC3Dg+WS+J7eWhXU3Qk8IIcnI\nehmaeaR11eOxd1GMcm0GAKCqL0WkUu93+aX1+jjS33VECJk9ngF4HrTKugCjeJDksnUx7W0AAxFZ\nB3BLVR+KiFugF2Eu1oX0fppuN0cIISVLODx4+8prlXUBdqceROQcgGU9/Ajz6840xRDFI7Bj9YQQ\nkoTs54BLRGQNwFkASyKyLiILAKCqF1GMklfNHW+PVfXzsl9ITwghqch6BGxj7mLz3kKsqpcDfSv1\nhBCSgpkZARNCyLzBh3Ja8KGchJC+4EM5oxj14K+pz9i+ITuf3pWlbOeUS8655ZRL17mF5CFdG9su\n+1b5HIdTEIQQkggWYEIISQQLMCGEJIIFmBBCEsECTAghiWABJoSQRLAAE0JIIliACSEkESzAhBCS\nCBZgQghJBAswIYQkggWYEEISwQJMCCGJ4HaUFtyOkhDSF9yOMopRD/6a+oztG7Lz6V1ZynZOueSc\nW065dJ1bSB7StbHtsm+Vz3E4BUEIIYlgASaEkERkPwUhIgMAawD2AZwAXj/p2LbZAPAUwKLRb9fR\nE0JICmZhBHxJVS+r6rYpvEPziHoAgIhsAXigqrdMYT0hIquxekIIScUsFOBVEfnAaj8FcNpqr6nq\nF1Z7B8D5GnpCCElC9lMQAIaq+txqnwDwTwAgIise+z0Awxg9IYSkJPsRsF18TUF9paqfGNEigF2n\ny76xPRahJ4SQZMzCCBgisgDg1wDeA3DOUg1gLqxZlAV3MUL/sttMCSEknpkowKr6AsA2gG0ReSAi\nV80FtX2PeVlwdyP0Hr60Xh8HsNQgY0LI0eYZgOdBq+xvRRaRgaruW+01ANdU9QdmSuK+qv7A0r+W\nhfSeWHmfDELIzDJztyKLyBDAHVOEy+kCMbpjqvpQRNxR7iKKlQ4I6f2MOsjc9dfUZ2zfkJ1P78pS\ntnPKJefccsql69xC8pCujW2Xfat8jpP7Rbh7KEa79lztaQA3Ldl1Z13vEMA1qx3SE0JIErIeAavq\nCxG5bu5kA4A3ADxW1UuWzUUR2TBFdtnoP4/VE0JIKrIuwACgqo8APArYXG6jJ4SQFOQ+BUEIIXML\nCzAhhCSCBZgQQhKR/TrgacJ1wISQvpi5dcBpGPXgr6nP2L4hO5/elaVs55RLzrnllEvXuYXkIV0b\n2y77Vvkch1MQhBCSCBZgQghJBAswIYQkghfhLHgRjhDSF7wIF8WoB39Nfcb2Ddn59K4sZTunXHLO\nLadcus4tJA/p2th22bfK5zicgiCEkESwABNCSCI4B2zBOWBCSF9wDjiKUQ/+mvqM7Ruy8+ldWcp2\nTrnknFtOuXSdW0ge0rWx7bJvlc9xOAVBCCGJYAEmhJBEcA7YgnPAhJC+4BxwFKMe/DX1Gds3ZOfT\nu7KU7ZxyyTm3nHLpOreQPKRrY9tl3yqf43AKghBCEjETI2DroZxvA7jnPuPN6J+ieOQ8VHW7jp4Q\nQpKgqlkfADad9n0AG1Z7C8C7tj2A1Vi941t58ODBo4/DV3OyvggnIgsAztkjXhFZA7ClqoumvVu+\nNu1TAC6o6pkYvRNPOQecop1TLjnnllMuXecWkod0bWy77DvZp+8iXO5zwG8A2BKR45ZsD8AAAERk\nxdNnD8AwRk8IISnJeg5YVZ+KyIqqPrfEpwHsmNeLAHadbvsAICLHQnpVfdl50kcEEcH6+l9AVXHl\nyr+i/CIlAiNfx5UrsOR++y59TbI/rFNcuSJOnw2o/iLD2K591Tms66uefSg+aUjqOd6a88EDFAX1\nuGmfBbDrsXkF4HhIzzng5sfGxoZ+++23+s033+jGxkZjeZe+UsaY9/cX0vEIH76alvUI2MMNFBfU\nnpv2vsemnO/djdB7+Cvr9XEASzVTdBlhHueAVf8dqj80r4cALqOY5/qFkX9n5D+utLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x16250b00>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"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)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"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",
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"Bxra, Byra, Bzra = MagSphereAnaFunA(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')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABEUAAAESCAYAAAAbqGtQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvX+wHUd17/ttl33kmBcsDmDIteWrHzaxqVD1bHAlTKXK\nqkhAxeWLzJOU4OKfm/csW8TGXMzzj9zcUglX7kUyQX4xcpB/5JKqeykThAFTLhMw5klV1DaJsbkp\nqGDwDwmEA4ZYP0ziyLKe+/0xM+f0md0/Vvd0z/TMXp+qU/uc6TU9vefsvabnO2utFlJKMAzDMAzD\nMAzDMAzDzBqn9D0AhmEYhmEYhmEYhmGYPmBRhGEYhmEYhmEYhmGYmYRFEYZhGIZhGIZhGIZhZhIW\nRRiGYRiGYRiGYRiGmUlYFGEYhmEYhmEYhmEYZiZhUYRhGIZhGIZhGIZhmJmERRGmE4QQO4UQrxp+\nviOE2Nj3GBmGYcYO+2KGYZh+YT/MMPlxat8DYGaOvQCeVf5eDuAPAewVQtwspfxEP8NiGIaZKdgX\nMwzD9Av7YYbJBCGl7HsMzAwghNgJ4EYA66WU32y0nQngAIDlUkqOXmIYhkkE+2KGYZh+YT/MMPnB\nXzamd6SUxwB8HgCEEP97z8NhGIaZSdgXMwzD9Av7YYbpBxZFmFwQS/4QYnWVW/n5xvaLq+2f7nZ4\nDMMwM0HTF99U+dyLpgyFeFwIcbi7oTEMw8wEYmqDuQbJq0KIG/sYJMOMCRZFmK7ROfrlAP4AwDNS\nyv8FAFLKZwHcDGBTo+DU3srug10MlmEYZqSQfDGAu6rXaxq2qwFcBOBvUg6SYRhmxFD9MABsavxs\nRlmPRAJ4PP1QGWbccKFVpmuuEUK8W/m7dv6vonTwC0gpPyGE+EMA9wghvgHgPwNYCeDtHY2VYRhm\nrJB8sZTymBDiiaptq2K/qXq9CwzDMEwIPnPiL6p/CyGuBrAawM5mXRKGYfxhUYTpmk2abRLAMQBr\nAPyvRttmAM8A+CbKp5I3NZRzhmEYxh8fX/w3AHYKIS6SUn632nYNpp9kMgzDMHR858QAylRyAHsA\nPCyl/JN0w2OY2YHTZ5iuWS+lPEX9AXAegO+gXIJsi2ospTyAcvJ9EcoJ+J93P2SGYZjR4eOL765e\nrwEWJuSrwFEiDMMwbfCaEwML6TWPoJwTv6fj8TLMaGFRhOmdSviowwSv0ZicV72+vlqqjGEYhomM\nyRdXqyE80WiTWBRLGIZhmAgQ5sSPo0yveVdng2KYGYBFESYLqkk3ACwRPYQQ61Gu5X4XylzLvR0P\njWEYZmYw+WKUPvh1Qoh1KHPevyGlfLHTwTEMw8wAljnxXpS19TZLKQ92PCyGGTUsijBZIISo8yqf\nULbVIsjj1WozNwNY31iNhmEYhomEzhdX1Muj34Vyos6pMwzDMAkwzIlvArARwM1cWJVh4pOs0Kqy\nZvYlAB6TUn5C0/4sgHkAkFLe49PODJatQohmDuRqlI5eohQ+au4F8FpUYYTN1WgUJZ1hGA3shxkL\nPr64XoXmGwDWAzjSXAmBYRgz7IsZAyQ/XEVN70D5GTigiCY1R6SUj6QeLMOMmSSiiBBih5TyFuXv\n7wghUF8EhBA7AXytVjqFEDuEEBullPdT2plBIqvXjZhel/0IyoiQm+twwMrh/x8oV5s5qNjWq9F8\nHgAXmGIYA+yHGQNevrjBXpSiyOc1bQzDaGBfzGjw9cOrlFddGvnDKIuvMgwTiJBSuq18OiwLYV6t\nquBV9eSdUsr56u/D9e/V3+tQfvnfTWlnGIZhzLAfZlJQhW/vAHAxL8XLMG7YFzMMwwyDFDVFXg9g\npxBipbLtCMoimfVSfk2OoHz65GxnGIZhnLAfZlJwDcplIFkQYRga7IsZhmEGQHRRREr5LMqnSAeV\nze9CGdoFlPmQhxu7HQUAIcRrCe0MwzCMBfbDTEyEEDuFEA+jDN3e2fd4GGYosC9mGIYZBklqiqhP\nkaoVRDYDqNXu5agKRSnUDn+e0M5LADIMwzhgP8xEZCOA16EM+b+378EwzJBgX8wwDJM/yVafUfg8\ngN9TVPKjGpva4R8mtC8ghIhbEIVhmKRIKZsFxUj4ftdDjzNikvlhgH3xDHGLEOIWtxmTO+yLe4Pn\nxAzDAGA/nBtJRREhxA4AOxr5x4dR5VIqLAcAKeWLQghre/MYH5K3WcdwCCs8R+3mB9u/gAu3N1fD\nSkusY67AIS/7v9v+dfz29m5ref1g+xfw7u2/3ekxAeDr2//OeVzf8+fiC9t/gE3bL4zaZ9/HNX3n\nbhKfatXvHUS761sdZXx04YcB4Db5IQDxvyMmYnyGQ8b6V9v/Cf/X9n/X6rgqK3HQafPJ7S/ho9vP\niHZMKn0cl3LMg1gZ/bgh/9e284uu/H9znOyL+6EbX0z978TkIQCXzcAx+zouv9dxHredh2Q/HJ8U\nhVYBAEKIjQC+riwhdhEASCmfwLTyPY8qv9LVTuUQViQRRIYOnxd/ViyctW5u9oYOn6t86NsP50zf\nn9GVOEgSRJhpcjl3fX+GqAxlnGOGfTHDMEzeJBFFhBDrUTrtx4UQy4UQqwH8oWJyd3WBqFkP4C6P\ndiN800+DzxMNnkyGw2JSv/Thh4fyv+5znLnc0I+BHM7lUD7zTH/0OSdmGIZhaERPn6mKSH29+lN1\n2nvrX6SUtwghbqyc/GoAT0spv0htb9L1zf0b1r610+OlPGZ97kwTu7PXrklyXBtr1p7d+TF1x+1i\nsvvWtW9Ifowcjss3Dt3Shx/umtDPcNvP4kVrfz1ovzY37+9ce1rwvm3o47ghx6zPbZu0mtD/K1B+\npkLmIX35f6Y7ZsEXA+fPyDH7Oi6/1/Eel8kJIeVw6zIJIeQV8r6+hzEacrlx7XscfR9/rFwpvtyq\nqJRP/iQXleoWIYS8T17R9zCc9PHd7juSYVZJUXfERe7Rl/X4bhKfYl88Qsrii33UFGEYxp/r2Q9n\nRherzzAD4RBWzLwgMOvvn2HGStff7VzFkLN/+UL0Pp974+uj99mWlTjYuTASGjHCMAzDMEy/sCjC\nLGGWhZFZfd8MM3ZmVRBJIYD4HKdvsSRGWo0vOQsjOY+NYRiGYfqERRFmilkURmbt/Q6R7hcuZsZA\nl9/tvsWQrkQQKs3x9CWS9BE1MmbYFzNMW96WsO/vJeybyQX2w/FhUYTR4irAOiZm4T0yzCwyZkEk\nNwGEQp8iSZfCSM4RGXy9Y5guSSl++ByThRKGcTEzosj+YhsA4NLJrdHtU/bdt70uamRvsRsAsHly\nnbNvH9vafg4ncN1kM8l+d1EWcA+1d00QtxX7AQC3Ti519u1jO2v2tW2fCCFuBPAsyqURIaW8p419\nhPblAG4B8Fhl8x0p5XdDxzsEuvxMUm7+thZPAgD2TC4g9a+zN4khG4pjAIAHJmeS+qba18JCcVn5\n9+QhUvdJ7cP7Lt+LSxyJdS5N6TQxPgdNamEkJz+s2vcJ++Ic2FW93pDAPmXfOdt/hmB7bfV6J7Hv\nWPYmoYQ/B/3Y73KbJIb98DSnxO6QGR+5PvFqwwoc4idmM4QQYieAx6WU91eOdE21vGGQfYT25QC+\nIaW8RUp5P4DlAP4kdLzMUrr6bncVHXL2L19Y+BkjXb+/rv5vfI2Zhn0xMw7epvycUf0MDXXs6vth\nxg77YcP75CV5GSpdTfBSH4cnqv3Qdkneh4m278L08mNCiMNSynnl73UAbpZSvttwPKt9hPa7ADwm\npbxXsTlTSnksZLx9k9OSvGMSRMYqglBJnV7TVTpNbg8W2BeP0xfzkrypmHWhgFNv0tBuSV72w/H9\n8MykzzDtGUMB1qGPn/FHCHGxZvMRAOtD7Nu2V2wB8HHVQHH+XuNluielIDLrQohKfS5SiSNdrU6T\nc42RLmFfzAyDWRdBmjTPB4skQ4b9sBkWRRgvhiyMDHXcTGvmARxubDsKAEKI10opX/Sxj9D+hmrb\nGiHE2yv75VLKTwSOl6lI/R1nMaQfuhBHeHWaTmBfzGQKCyF01HPFAskAYT9sgGuKMN7wEy9mYCxH\nVZhJoXawze0U+7btq6vfpZIfCSHEjsDxMhiuIDLmWiGxSXmuUqdDsSgPgH0xkx1cR6MdXItkgLAf\nNsCRIkwQQ4sYGdJYGT2mNdknAB6173pUs612pE31mWLftr0+5neU9keqv28JGC+TmJSCSKc8m7j/\n1W6TGKSKHEkdMTKWNBr2xcbxMoOAb+DTUJ9Xjh7pAvbDxvEGw6IIM3pYEBk3RfVTc/u0yWGUSrPK\ncgAwhN1Z7YUQbduPao6thhL6jnfmSfkdTyGIdCaGpBZBbMfrQCA5+5cvDE4YGTPsi5m8YTGkGzi9\npk/YD4fD6TNMMGN44sWMHynlE5hWmucBaIt3u+wjtD8L4KgQYpXSvuDgfcc767AgovBs46dPOhpL\nipSalKk0syzSsy9m+oHTO/qDz31usB82w6II04rchZFZnoAyS7i7sab5egB31X8IIVY32q32Edo/\njqWVs/8QwE0e+zOJiX1jnKQeRk4iiIvEY2VhZDCwL2Y6gm/I84H/D5nBfliDkFLG7rMzhBDyCnlf\n38OYeWJP8GL1l9vEM9Z4cheiTFwpvtxqTfafEm3PwfSa7FUfN6K8HVsN4EhjPfQtADZJKd9DsY/Y\nXiOllH/us39OCCHkffKKzo+b6jueQhCJSu4CiA8J0mxiptSkSqXp04+zLx6nLxZCSOCOvoeRAXwD\nnjecUlNyPfvh6faaXvxwMlFECLEJwDuklLdotq8GsBflOsNbAHxBSnlAsanf+DwA1JVoNccgiyL7\ni20AgEsnt0a3D+n7BJbhgskekv2TxVYAmLI33RCkfK8me9NY9ha7AQCbJ9eR+t5b7Line truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x17ede9b0>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"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)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"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",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"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"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Second, set problem class then pair with survey"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"prob = MagneticsDiffSecondary(mesh)\n",
|
||||
"prob.pair(survey)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Third, run forward modeling"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"data = survey.dpred(mu)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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",
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"ename": "NameError",
|
||||
"evalue": "name 'Bzra' is not defined",
|
||||
"output_type": "error",
|
||||
"traceback": [
|
||||
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
|
||||
"\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)",
|
||||
"\u001b[1;32m<ipython-input-11-0e249253afbe>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[0mfig\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0max\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfigsize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m18\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m4\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mvmin\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mBzra\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 3\u001b[0m \u001b[0mvmax\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mBzra\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 4\u001b[0m \u001b[0mresidual\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mxr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msize\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0myr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msize\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0morder\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'F'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mBzra\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[0mdat0\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0max\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcontourf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mY\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mBzra\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m30\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvmin\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mvmin\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvmax\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mvmax\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;31mNameError\u001b[0m: name 'Bzra' is not defined"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABBwAAAEHCAYAAAANslX6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFm1JREFUeJzt3UGSG9eZLeDzOzSWytUbaNErEJ81r2jS9gIotzfQtPzm\nZjd7ZM6aDHkBluj565CDntuyFTV+ouVegEx5AyZpLeD9b1AADUGoQgF1oUQR3xeBiELiZuEmK3mK\ncZh5Ud0dAAAAgJG+M/UEAAAAgNePwgEAAAAYTuEAAAAADKdwAAAAAIZTOAAAAADDKRwAAACA4d64\nzKCqei/J97v7/iXG3kvyLMlxknT34yvNEIAkshhganIYYDMXFg5VdSvJzSQ/SPKXdd+sqh4l+V13\nfzp7/rCq7nT3kxGTBThEshhgWnIYYDvV3esHVT1MctTdP1sz7nl3Hy88v5XkP7r7h1eeKcCBk8UA\n05LDAJsZtoZDVd1csflFktuj3gOAi8ligGnJYYB/GLlo5HGS50vbXiZJVb058H0AOJ8sBpiWHAaY\nGVk4HGW2KM6CedgubwdgN2QxwLTkMMDMyMLh5Ypt81BdbnkB2A1ZDDAtOQwwc6mPxbyk5zlrdBcd\nJUl3f7Vqh6pav2IlwES6u6aewxZkMfDakMMA07tKFg+7wqG7P883G93jJJ+s2e8gHr/4xS8mn4Nj\ndayO9fKP66plsfPYsb7Wj0M61uuq5bBz2LG+1o9DOtbuq2fxZQuHlY1GVd2oqjsLmz5aen47yYfb\nTg6Ar5HFANOSwwAbuPCWiqp6J2cBeSfJd6vqL0n+0N1/ng25leS9JE+SpLvvV9W9WcDeSPJFd/92\nZ7MHOACyGGBachhgOxcWDrMQ/XOSD855/XGSx0vbVo49dCcnJ1NP4VvjWF9Ph3Ss+0YWj3NI57Fj\nfT0d0rHuEzk8ziGdw4719XRIxzpCjbgvY+s3r+op3x/gPFWVvp6LlW1MFgP7SA4DTO+qWTzyYzEB\nAAAAkigcAAAAgB1QOAAAAADDKRwAAACA4RQOAAAAwHAKBwAAAGA4hQMAAAAwnMIBAAAAGE7hAAAA\nAAyncAAAAACGUzgAAAAAwykcAAAAgOEUDgAAAMBwCgcAAABgOIUDAAAAMJzCAQAAABhO4QAAAAAM\np3AAAAAAhlM4AAAAAMMpHAAAAIDhFA4AAADAcAoHAAAAYDiFAwAAADCcwgEAAAAYTuEAAAAADKdw\nAAAAAIZTOAAAAADDKRwAAACA4RQOAAAAwHAKBwAAAGA4hQMAAAAwnMIBAAAAGO6NywyqqntJniU5\nTpLufnyJ8S9nT4+6+4OrTBIAWQwwNTkMsJnq7osHVD1K8rvu/nT2/GGSz7r7yTnj7y2GaVW9k+T2\nqoCtql73/gBTqKp0d009jzlZDBwaOQwwvatm8WVuqbg7D9aZT5K8f8H4nyw+6e4/J3l3i7kB8A+y\nGGBachhgQxcWDlV1c8XmF0luX7Db86r6uKremn2PO0n+e/spAhw2WQwwLTkMsJ11VzgcJ3m+tO1l\nklTVm+fs836Sm0m+nN23lu7+7VUmCXDgZDHAtOQwwBbWFQ5HmS2Ks2AetsvbkyTd/WWSD5M8TfIo\nLh0DuCpZDDAtOQywhXWfUvFyxbZ5qC63vEleLajzf7r7g6q6leQ3VXWju/911fgHDx68+vrk5CQn\nJyfr5gww3OnpaU5PT6eexnlkMfDak8MPXn0th4GpjM7iCz+lYna/2tPu/s5F25Ze+9fuvr+w7a0k\nX3b3N9pfK/IC+2qfVkeXxcAhksMA09vpp1R09+f5ZqN7nLNVeVf5bpK/LX2Pvyf5w7YTBDh0shhg\nWnIYYDuX+VjMj2ar6s7dztn9aEmSqroxf727/5jkB4s7V9VRkmcD5gpwyGQxwLTkMMCGLryl4tWg\ns5V1nyW5keRFd/964bW7Sd7r7h/Nnr+ds1V5/5ZZE9zdj8/5vi4fA/bSPl3KOyeLgUMihwGmd9Us\nvlThsCvCFdhX+/gP3V2RxcA+ksMA09vpGg4AAAAA21A4AAAAAMMpHAAAAIDhFA4AAADAcAoHAAAA\nYDiFAwAAADCcwgEAAAAYTuEAAAAADKdwAAAAAIZTOAAAAADDKRwAAACA4RQOAAAAwHAKBwAAAGA4\nhQMAAAAwnMIBAAAAGE7hAAAAAAyncAAAAACGUzgAAAAAwykcAAAAgOEUDgAAAMBwCgcAAABgOIUD\nAAAAMJzCAQAAABhO4QAAAAAMp3AAAAAAhlM4AAAAAMMpHAAAAIDhFA4AAADAcAoHAAAAYDiFAwAA\nADCcwgEAAAAYTuEAAAAADKdwAAAAAIZ74zKDqupekmdJjpOkux+vGX+U5H6Sz2b7PO3uP19tqgCH\nTRYDTEsOA2xm7RUOVfUoyZ+6+8ksVL9XVXcuGH+U5A/dfb+7nyQ5SvKfw2YMcIBkMcC05DDA5qq7\nLx5Q9by7jxee30ryH939w3PGf5jks+7+9cK2t7r77yvG9rr3B5hCVaW7a+p5zMli4NDIYYDpXTWL\nL7zCoapurtj8IsntC3a7m+QPixtWBSsAlyOLAaYlhwG2s24Nh+Mkz5e2vUySqnqzu79afKGqbsy+\n/F5V/a/Z/kfd/cGIyQIcKFkMMC05DLCFdWs4HGW2KM6Cedgub0+Sebj2wv1tqaqH208R4ODJYoBp\nyWGALay7wuHlim3zUF1ueRe3PV3Y9sfZ8/ur3uDBgwevvj45OcnJycmaKQGMd3p6mtPT06mncR5Z\nDLz25PCDV1/LYWAqo7P4wkUjZ/erPe3u71y0beG1G0m+WBp/I8kXObuMbPlyMwvkAHtpnxYrk8XA\nIZLDANPb6aKR3f15vtnoHif55Jzxz5K8rKq3FzYfzV77atU+AFxMFgNMSw4DbGfdGg5J8tHSZwzf\nTvLh/ElV3Vh6/b/y9RV7f5Lk3680SwBkMcC05DDAhi68peLVoKp7SZ7lbAGcF0ufJ3w3yXvd/aOl\n8XPd3b885/u6fAzYS/t0Ke+cLAYOiRwGmN5Vs/hShcOuCFdgX+3jP3R3RRYD+0gOA0xvp2s4AAAA\nAGxD4QAAAAAMp3AAAAAAhlM4AAAAAMMpHAAAAIDhFA4AAADAcAoHAAAAYDiFAwAAADCcwgEAAAAY\nTuEAAAAADKdwAAAAAIZTOAAAAADDKRwAAACA4RQOAAAAwHAKBwAAAGA4hQMAAAAwnMIBAAAAGE7h\nAAAAAAyncAAAAACGUzgAAAAAwykcAAAAgOEUDgAAAMBwCgcAAABgOIUDAAAAMJzCAQAAABhO4QAA\nAAAMp3AAAAAAhlM4AAAAAMMpHAAAAIDhFA4AAADAcAoHAAAAYLg3LjOoqu4leZbkOEm6+/Fl36Cq\nftXdP9tuegDMyWKAaclhgM2svcKhqh4l+VN3P5mF6veq6s5lvvls3+9fcY4AB08WA0xLDgNs7jK3\nVNzt7k8Xnn+S5P11O1XVjSS97cQA+BpZDDAtOQywoQsLh6q6uWLziyS3L/G9b+UsiAG4AlkMMC05\nDLCddVc4HCd5vrTtZZJU1Zvn7VRVt5J8nKSuNDsAElkMMDU5DLCFdYXDUWaL4iyYh+3y9q/t191/\n33pWACySxQDTksMAW1hXOLxcsW0eqsstb5Kkqu5095MrzQqARbIYYFpyGGAL6wqH5zlrdBcdJUl3\nf7U8uKrezupABmB7shhgWnIYYAtvXPRid39eVctheZzzF765meTGwsI67yY5qqqfJ3nS3V8u7/Dg\nwYNXX5+cnOTk5ORyMwcY6PT0NKenp1NPYyVZDBwCOfzg1ddyGJjK6Cyu7os/paeqHib5bH5J2Oz5\n/+3u386e30jyzqpLxqrqp0l+2t0rP3e4qnrd+wNMoarS3XuzyJcsBg6NHAaY3lWzeN0tFenu+zlr\naO9U1b0kX8yDdeZWkp+umNjdJO8lebuqfl5Vb207SYBDJ4sBpiWHATa39gqHnb65NhfYU/v2P2u7\nJIuBfSSHAaa38yscAAAAADalcAAAAACGUzgAAAAAwykcAAAAgOEUDgAAAMBwCgcAAABgOIUDAAAA\nMJzCAQAAABhO4QAAAAAMp3AAAAAAhlM4AAAAAMMpHAAAAIDhFA4AAADAcAoHAAAAYDiFAwAAADCc\nwgEAAAAYTuEAAAAADKdwAAAAAIZTOAAAAADDKRwAAACA4RQOAAAAwHAKBwAAAGA4hQMAAAAwnMIB\nAAAAGE7hAAAAAAyncAAAAACGUzgAAAAAwykcAAAAgOEUDgAAAMBwCgcAAABgOIUDAAAAMJzCAQAA\nABhO4QAAAAAMp3AAAAAAhnvjMoOq6l6SZ0mOk6S7H19ifJK8m+Sz7v7gKpMEQBYDTE0OA2xmbeFQ\nVY+S/K67P509f1hVd7r7yTnjH3b3/YXnT6sqAhZge7IYYFpyGGBzl7ml4u48WGc+SfL+qoFV9VaS\nvy1t/jDJf243PQBmZDHAtOQwwIYuLByq6uaKzS+S3D5nl39K8qiq/nlp/NE2kwNAFgNMTQ4DbGfd\nFQ7HSZ4vbXuZJFX15vLg7n6W5GZ3/3Vh8w9y1gADsB1ZDDAtOQywhXVrOBxltijOgnnYHif5anmH\n7v6f+ddVdZTkx0lWtcIAXI4sBpiWHAbYwrrC4eWKbfOwXW55V/k4yb8stbtf8+DBg1dfn5yc5OTk\n5BLfFmCs09PTnJ6eTj2N88hi4LUnhx+8+loOA1MZncXV3ee/eHa/2tPu/s5F287Z92GS3y8trrM8\npi96f4CpzFYSr6nnkchi4DDJYYDpXTWLLwzI7v4832x0j7Pm/rOqupOFYK2qd7adIMChk8UA05LD\nANu5zMdifjQLy7nbOftYnyRJVd1YfL2qbucsgP9UVUdVdSPJT0ZNGOBAyWKAaclhgLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x17684eb8>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"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')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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": {
|
||||
"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.10"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 0
|
||||
}
|
||||
Reference in new issue
Block a user